All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
direct_solve.hpp
Go to the documentation of this file.
#ifndef VIENNACL_LINALG_HOST_BASED_DIRECT_SOLVE_HPP
#define VIENNACL_LINALG_HOST_BASED_DIRECT_SOLVE_HPP
/* =========================================================================
Copyright (c) 2010-2016, Institute for Microelectronics,
Institute for Analysis and Scientific Computing,
TU Wien.
Portions of this software are copyright by UChicago Argonne, LLC.
-----------------
ViennaCL - The Vienna Computing Library
-----------------
Project Head: Karl Rupp rupp@iue.tuwien.ac.at
(A list of authors and contributors can be found in the manual)
License: MIT (X11), see file LICENSE in the base directory
============================================================================= */
#include "viennacl/vector.hpp"
#include "viennacl/matrix.hpp"
#include "viennacl/linalg/host_based/common.hpp"
namespace viennacl
{
namespace linalg
{
namespace host_based
{
namespace detail
{
//
// Upper solve:
//
template<typename MatrixT1, typename MatrixT2>
void upper_inplace_solve_matrix(MatrixT1 & A, MatrixT2 & B, vcl_size_t A_size, vcl_size_t B_size, bool unit_diagonal)
{
typedef typename MatrixT2::value_type value_type;
for (vcl_size_t i = 0; i < A_size; ++i)
{
vcl_size_t current_row = A_size - i - 1;
for (vcl_size_t j = current_row + 1; j < A_size; ++j)
{
value_type A_element = A(current_row, j);
for (vcl_size_t k=0; k < B_size; ++k)
B(current_row, k) -= A_element * B(j, k);
}
if (!unit_diagonal)
{
value_type A_diag = A(current_row, current_row);
for (vcl_size_t k=0; k < B_size; ++k)
B(current_row, k) /= A_diag;
}
}
}
template<typename MatrixT1, typename MatrixT2>
void inplace_solve_matrix(MatrixT1 & A, MatrixT2 & B, vcl_size_t A_size, vcl_size_t B_size, viennacl::linalg::unit_upper_tag)
{
upper_inplace_solve_matrix(A, B, A_size, B_size, true);
}
template<typename MatrixT1, typename MatrixT2>
void inplace_solve_matrix(MatrixT1 & A, MatrixT2 & B, vcl_size_t A_size, vcl_size_t B_size, viennacl::linalg::upper_tag)
{
upper_inplace_solve_matrix(A, B, A_size, B_size, false);
}
//
// Lower solve:
//
template<typename MatrixT1, typename MatrixT2>
void lower_inplace_solve_matrix(MatrixT1 & A, MatrixT2 & B, vcl_size_t A_size, vcl_size_t B_size, bool unit_diagonal)
{
typedef typename MatrixT2::value_type value_type;
for (vcl_size_t i = 0; i < A_size; ++i)
{
for (vcl_size_t j = 0; j < i; ++j)
{
value_type A_element = A(i, j);
for (vcl_size_t k=0; k < B_size; ++k)
B(i, k) -= A_element * B(j, k);
}
if (!unit_diagonal)
{
value_type A_diag = A(i, i);
for (vcl_size_t k=0; k < B_size; ++k)
B(i, k) /= A_diag;
}
}
}
template<typename MatrixT1, typename MatrixT2>
void inplace_solve_matrix(MatrixT1 & A, MatrixT2 & B, vcl_size_t A_size, vcl_size_t B_size, viennacl::linalg::unit_lower_tag)
{
lower_inplace_solve_matrix(A, B, A_size, B_size, true);
}
template<typename MatrixT1, typename MatrixT2>
void inplace_solve_matrix(MatrixT1 & A, MatrixT2 & B, vcl_size_t A_size, vcl_size_t B_size, viennacl::linalg::lower_tag)
{
lower_inplace_solve_matrix(A, B, A_size, B_size, false);
}
}
//
// Note: By convention, all size checks are performed in the calling frontend. No need to double-check here.
//
template<typename NumericT, typename SolverTagT>
void inplace_solve(matrix_base<NumericT> const & A,
matrix_base<NumericT> & B,
SolverTagT)
{
typedef NumericT value_type;
value_type const * data_A = detail::extract_raw_pointer<value_type>(A);
value_type * data_B = detail::extract_raw_pointer<value_type>(B);
vcl_size_t A_start1 = viennacl::traits::start1(A);
vcl_size_t A_start2 = viennacl::traits::start2(A);
vcl_size_t A_inc1 = viennacl::traits::stride1(A);
vcl_size_t A_inc2 = viennacl::traits::stride2(A);
//vcl_size_t A_size1 = viennacl::traits::size1(A);
vcl_size_t A_size2 = viennacl::traits::size2(A);
vcl_size_t A_internal_size1 = viennacl::traits::internal_size1(A);
vcl_size_t A_internal_size2 = viennacl::traits::internal_size2(A);
vcl_size_t B_start1 = viennacl::traits::start1(B);
vcl_size_t B_start2 = viennacl::traits::start2(B);
vcl_size_t B_inc1 = viennacl::traits::stride1(B);
vcl_size_t B_inc2 = viennacl::traits::stride2(B);
//vcl_size_t B_size1 = viennacl::traits::size1(B);
vcl_size_t B_size2 = viennacl::traits::size2(B);
vcl_size_t B_internal_size1 = viennacl::traits::internal_size1(B);
vcl_size_t B_internal_size2 = viennacl::traits::internal_size2(B);
if (A.row_major() && B.row_major())
{
detail::matrix_array_wrapper<value_type const, row_major, false> wrapper_A(data_A, A_start1, A_start2, A_inc1, A_inc2, A_internal_size1, A_internal_size2);
detail::matrix_array_wrapper<value_type, row_major, false> wrapper_B(data_B, B_start1, B_start2, B_inc1, B_inc2, B_internal_size1, B_internal_size2);
detail::inplace_solve_matrix(wrapper_A, wrapper_B, A_size2, B_size2, SolverTagT());
}
else if (A.row_major() && !B.row_major())
{
detail::matrix_array_wrapper<value_type const, row_major, false> wrapper_A(data_A, A_start1, A_start2, A_inc1, A_inc2, A_internal_size1, A_internal_size2);
detail::matrix_array_wrapper<value_type, column_major, false> wrapper_B(data_B, B_start1, B_start2, B_inc1, B_inc2, B_internal_size1, B_internal_size2);
detail::inplace_solve_matrix(wrapper_A, wrapper_B, A_size2, B_size2, SolverTagT());
}
else if (!A.row_major() && B.row_major())
{
detail::matrix_array_wrapper<value_type const, column_major, false> wrapper_A(data_A, A_start1, A_start2, A_inc1, A_inc2, A_internal_size1, A_internal_size2);
detail::matrix_array_wrapper<value_type, row_major, false> wrapper_B(data_B, B_start1, B_start2, B_inc1, B_inc2, B_internal_size1, B_internal_size2);
detail::inplace_solve_matrix(wrapper_A, wrapper_B, A_size2, B_size2, SolverTagT());
}
else
{
detail::matrix_array_wrapper<value_type const, column_major, false> wrapper_A(data_A, A_start1, A_start2, A_inc1, A_inc2, A_internal_size1, A_internal_size2);
detail::matrix_array_wrapper<value_type, column_major, false> wrapper_B(data_B, B_start1, B_start2, B_inc1, B_inc2, B_internal_size1, B_internal_size2);
detail::inplace_solve_matrix(wrapper_A, wrapper_B, A_size2, B_size2, SolverTagT());
}
}
//
// Solve on vector
//
namespace detail
{
//
// Upper solve:
//
template<typename MatrixT, typename VectorT>
void upper_inplace_solve_vector(MatrixT & A, VectorT & b, vcl_size_t A_size, bool unit_diagonal)
{
typedef typename VectorT::value_type value_type;
for (vcl_size_t i = 0; i < A_size; ++i)
{
vcl_size_t current_row = A_size - i - 1;
for (vcl_size_t j = current_row + 1; j < A_size; ++j)
{
value_type A_element = A(current_row, j);
b(current_row) -= A_element * b(j);
}
if (!unit_diagonal)
b(current_row) /= A(current_row, current_row);
}
}
template<typename MatrixT, typename VectorT>
void inplace_solve_vector(MatrixT & A, VectorT & b, vcl_size_t A_size, viennacl::linalg::unit_upper_tag)
{
upper_inplace_solve_vector(A, b, A_size, true);
}
template<typename MatrixT, typename VectorT>
void inplace_solve_vector(MatrixT & A, VectorT & b, vcl_size_t A_size, viennacl::linalg::upper_tag)
{
upper_inplace_solve_vector(A, b, A_size, false);
}
//
// Lower solve:
//
template<typename MatrixT, typename VectorT>
void lower_inplace_solve_vector(MatrixT & A, VectorT & b, vcl_size_t A_size, bool unit_diagonal)
{
typedef typename VectorT::value_type value_type;
for (vcl_size_t i = 0; i < A_size; ++i)
{
for (vcl_size_t j = 0; j < i; ++j)
{
value_type A_element = A(i, j);
b(i) -= A_element * b(j);
}
if (!unit_diagonal)
b(i) /= A(i, i);
}
}
template<typename MatrixT, typename VectorT>
void inplace_solve_vector(MatrixT & A, VectorT & b, vcl_size_t A_size, viennacl::linalg::unit_lower_tag)
{
lower_inplace_solve_vector(A, b, A_size, true);
}
template<typename MatrixT, typename VectorT>
void inplace_solve_vector(MatrixT & A, VectorT & b, vcl_size_t A_size, viennacl::linalg::lower_tag)
{
lower_inplace_solve_vector(A, b, A_size, false);
}
}
template<typename NumericT, typename SolverTagT>
void inplace_solve(matrix_base<NumericT> const & mat,
vector_base<NumericT> & vec,
SolverTagT)
{
typedef NumericT value_type;
value_type const * data_A = detail::extract_raw_pointer<value_type>(mat);
value_type * data_v = detail::extract_raw_pointer<value_type>(vec);
vcl_size_t A_start1 = viennacl::traits::start1(mat);
vcl_size_t A_start2 = viennacl::traits::start2(mat);
vcl_size_t A_inc1 = viennacl::traits::stride1(mat);
vcl_size_t A_inc2 = viennacl::traits::stride2(mat);
vcl_size_t A_size2 = viennacl::traits::size2(mat);
vcl_size_t A_internal_size1 = viennacl::traits::internal_size1(mat);
vcl_size_t A_internal_size2 = viennacl::traits::internal_size2(mat);
vcl_size_t start1 = viennacl::traits::start(vec);
vcl_size_t inc1 = viennacl::traits::stride(vec);
if (mat.row_major())
{
detail::matrix_array_wrapper<value_type const, row_major, false> wrapper_A(data_A, A_start1, A_start2, A_inc1, A_inc2, A_internal_size1, A_internal_size2);
detail::vector_array_wrapper<value_type> wrapper_v(data_v, start1, inc1);
detail::inplace_solve_vector(wrapper_A, wrapper_v, A_size2, SolverTagT());
}
else
{
detail::matrix_array_wrapper<value_type const, column_major, false> wrapper_A(data_A, A_start1, A_start2, A_inc1, A_inc2, A_internal_size1, A_internal_size2);
detail::vector_array_wrapper<value_type> wrapper_v(data_v, start1, inc1);
detail::inplace_solve_vector(wrapper_A, wrapper_v, A_size2, SolverTagT());
}
}
} // namespace host_based
} // namespace linalg
} // namespace viennacl
#endif
viennacl::linalg::host_based::detail::vector_array_wrapper
Helper class for accessing a strided subvector of a larger vector.
Definition: common.hpp:50
viennacl::traits::stride1
result_of::size_type< matrix_base< NumericT > >::type stride1(matrix_base< NumericT > const &s)
Definition: stride.hpp:55
matrix.hpp
Implementation of the dense matrix class.
viennacl::traits::internal_size1
vcl_size_t internal_size1(matrix_base< NumericT > const &mat)
Helper routine for obtaining the internal number of entries per row of a ViennaCL matrix...
Definition: size.hpp:386
viennacl::linalg::host_based::detail::lower_inplace_solve_vector
void lower_inplace_solve_vector(MatrixT &A, VectorT &b, vcl_size_t A_size, bool unit_diagonal)
Definition: direct_solve.hpp:234
viennacl::linalg::lower_tag
A tag class representing a lower triangular matrix.
Definition: forwards.h:849
viennacl::linalg::host_based::inplace_solve
void inplace_solve(matrix_base< NumericT > const &A, matrix_base< NumericT > &B, SolverTagT)
Direct inplace solver for triangular systems with multiple right hand sides, i.e. A \ B (MATLAB notat...
Definition: direct_solve.hpp:130
viennacl::traits::internal_size2
vcl_size_t internal_size2(matrix_base< NumericT > const &mat)
Helper routine for obtaining the internal number of entries per column of a ViennaCL matrix...
Definition: size.hpp:394
viennacl::matrix_base< NumericT >
viennacl::traits::stride
result_of::size_type< viennacl::vector_base< T > >::type stride(viennacl::vector_base< T > const &s)
Definition: stride.hpp:45
viennacl::traits::start1
result_of::size_type< T >::type start1(T const &obj)
Definition: start.hpp:65
viennacl::linalg::host_based::detail::upper_inplace_solve_vector
void upper_inplace_solve_vector(MatrixT &A, VectorT &b, vcl_size_t A_size, bool unit_diagonal)
Definition: direct_solve.hpp:199
viennacl::traits::size2
result_of::size_type< MatrixType >::type size2(MatrixType const &mat)
Generic routine for obtaining the number of columns of a matrix (ViennaCL, uBLAS, etc...
Definition: size.hpp:201
NumericT
float NumericT
Definition: bisect.cpp:40
viennacl::linalg::host_based::detail::lower_inplace_solve_matrix
void lower_inplace_solve_matrix(MatrixT1 &A, MatrixT2 &B, vcl_size_t A_size, vcl_size_t B_size, bool unit_diagonal)
Definition: direct_solve.hpp:83
viennacl::traits::start2
result_of::size_type< T >::type start2(T const &obj)
Definition: start.hpp:84
viennacl::linalg::host_based::detail::matrix_array_wrapper
Helper array for accessing a strided submatrix embedded in a larger matrix.
Definition: common.hpp:73
viennacl::linalg::upper_tag
A tag class representing an upper triangular matrix.
Definition: forwards.h:854
viennacl::traits::start
result_of::size_type< T >::type start(T const &obj)
Definition: start.hpp:44
viennacl::vector_base< NumericT >
viennacl::vcl_size_t
std::size_t vcl_size_t
Definition: forwards.h:75
common.hpp
Common routines for single-threaded or OpenMP-enabled execution on CPU.
viennacl::traits::stride2
result_of::size_type< matrix_base< NumericT > >::type stride2(matrix_base< NumericT > const &s)
Definition: stride.hpp:65
vector.hpp
The vector type with operator-overloads and proxy classes is defined here. Linear algebra operations ...
viennacl::matrix_base::row_major
bool row_major() const
Definition: matrix_def.hpp:248
viennacl::linalg::host_based::detail::upper_inplace_solve_matrix
void upper_inplace_solve_matrix(MatrixT1 &A, MatrixT2 &B, vcl_size_t A_size, vcl_size_t B_size, bool unit_diagonal)
Definition: direct_solve.hpp:43
viennacl::linalg::host_based::detail::inplace_solve_vector
void inplace_solve_vector(MatrixT &A, VectorT &b, vcl_size_t A_size, viennacl::linalg::unit_upper_tag)
Definition: direct_solve.hpp:219
viennacl::linalg::unit_lower_tag
A tag class representing a lower triangular matrix with unit diagonal.
Definition: forwards.h:859
viennacl::linalg::host_based::detail::inplace_solve_matrix
void inplace_solve_matrix(MatrixT1 &A, MatrixT2 &B, vcl_size_t A_size, vcl_size_t B_size, viennacl::linalg::unit_upper_tag)
Definition: direct_solve.hpp:68
viennacl::linalg::unit_upper_tag
A tag class representing an upper triangular matrix with unit diagonal.
Definition: forwards.h:864