1 #ifndef VIENNACL_LINALG_QR_HPP
2 #define VIENNACL_LINALG_QR_HPP
33 #include "boost/numeric/ublas/vector.hpp"
34 #include "boost/numeric/ublas/matrix.hpp"
35 #include "boost/numeric/ublas/matrix_proxy.hpp"
36 #include "boost/numeric/ublas/vector_proxy.hpp"
37 #include "boost/numeric/ublas/io.hpp"
38 #include "boost/numeric/ublas/matrix_expression.hpp"
52 template <
typename MatrixType,
typename VectorType>
58 typedef typename MatrixType::value_type ScalarType;
64 ScalarType sigma = matrix_1x1(0,0);
66 ScalarType A_jj = A(j,j);
68 assert( sigma >= 0.0 &&
bool(
"sigma must be non-negative!"));
78 ScalarType mu = std::sqrt(sigma + A_jj*A_jj);
80 ScalarType v1 = (A_jj <= 0) ? (A_jj - mu) : (-sigma / (A_jj + mu));
81 beta =
static_cast<ScalarType
>(2.0) * v1 * v1 / (sigma + v1 * v1);
91 template <
typename MatrixType,
typename VectorType>
104 ScalarType sigma = matrix_1x1(0,0);
106 ScalarType A_jj = A(j,j);
108 assert( sigma >= 0.0 &&
bool(
"sigma must be non-negative!"));
118 ScalarType mu = std::sqrt(sigma + A_jj*A_jj);
120 ScalarType v1 = (A_jj <= 0) ? (A_jj - mu) : (-sigma / (A_jj + mu));
122 beta = 2.0 * v1 * v1 / (sigma + v1 * v1);
133 template <
typename MatrixType,
typename VectorType,
typename ScalarType>
136 ScalarType v_in_col = A(j,k);
138 v_in_col += v[i] * A(i,k);
143 A(i,k) -= beta * v_in_col * v[i];
146 template <
typename MatrixType,
typename VectorType,
typename ScalarType>
152 ScalarType v_in_col = A(j,k);
155 v_in_col += matrix_1x1(0,0);
160 template <
typename MatrixType,
typename VectorType,
typename ScalarType>
166 ScalarType v_in_col = A(j,k);
170 v_in_col += matrix_1x1(0,0);
172 if ( beta * v_in_col != 0.0)
183 template <
typename MatrixType,
typename VectorType,
typename ScalarType>
193 template <
typename MatrixType,
typename VectorType>
200 template <
typename MatrixType,
typename VectorType>
210 template <
typename MatrixType,
typename VectorType>
227 template<
typename MatrixType>
230 typedef typename MatrixType::value_type ScalarType;
231 typedef boost::numeric::ublas::matrix_range<MatrixType> MatrixRange;
236 std::vector<ScalarType> betas(A.size2());
237 MatrixType v(A.size1(), 1);
238 MatrixType matrix_1x1(1,1);
240 MatrixType Y(A.size1(), block_size); Y.clear(); Y.resize(A.size1(), block_size);
241 MatrixType W(A.size1(), block_size); W.clear(); W.resize(A.size1(), block_size);
244 for (
vcl_size_t j = 0; j < std::min(A.size1(), A.size2()); j += block_size)
246 vcl_size_t effective_block_size = std::min(std::min(A.size1(), A.size2()), j+block_size) - j;
249 for (
vcl_size_t k = 0; k < effective_block_size; ++k)
253 for (
vcl_size_t l = k; l < effective_block_size; ++l)
262 Y.clear(); Y.resize(A.size1(), block_size);
263 for (
vcl_size_t k = 0; k < effective_block_size; ++k)
275 W.clear(); W.resize(A.size1(), block_size);
281 for (
vcl_size_t k = 1; k < effective_block_size; ++k)
289 z = - betas[j+k] * (v_k +
prod(W_old, YT_prod_v));
296 if (A.size2() - j - effective_block_size > 0)
299 MatrixRange A_part(A,
range(j, A.size1()),
range(j+effective_block_size, A.size2()));
300 MatrixRange W_part(W,
range(j, A.size1()),
range(0, effective_block_size));
320 template<
typename MatrixType>
321 std::vector< typename viennacl::result_of::cpu_value_type< typename MatrixType::value_type >::type >
330 std::vector<ScalarType> betas(A.size2());
331 MatrixType v(A.size1(), 1);
332 MatrixType matrix_1x1(1,1);
334 MatrixType Y(A.size1(), block_size); Y.clear();
335 MatrixType W(A.size1(), block_size); W.clear();
337 MatrixType YT_prod_v(block_size, 1);
338 MatrixType z(A.size1(), 1);
341 for (
vcl_size_t j = 0; j < std::min(A.size1(), A.size2()); j += block_size)
343 vcl_size_t effective_block_size = std::min(std::min(A.size1(), A.size2()), j+block_size) - j;
346 for (
vcl_size_t k = 0; k < effective_block_size; ++k)
349 for (
vcl_size_t l = k; l < effective_block_size; ++l)
359 for (
vcl_size_t k = 0; k < effective_block_size; ++k)
379 for (
vcl_size_t k = 1; k < effective_block_size; ++k)
396 if (A.size2() > j + effective_block_size)
399 MatrixRange A_part(A,
range(j, A.size1()),
range(j+effective_block_size, A.size2()));
400 MatrixRange W_part(W,
range(j, A.size1()),
range(0, effective_block_size));
401 MatrixType temp =
prod(
trans(W_part), A_part);
424 template<
typename MatrixType>
425 std::vector< typename viennacl::result_of::cpu_value_type< typename MatrixType::value_type >::type >
431 typedef boost::numeric::ublas::matrix<ScalarType> UblasMatrixType;
432 typedef boost::numeric::ublas::matrix_range<UblasMatrixType> UblasMatrixRange;
434 std::vector<ScalarType> betas(A.size2());
435 UblasMatrixType v(A.size1(), 1);
436 UblasMatrixType matrix_1x1(1,1);
438 UblasMatrixType ublasW(A.size1(), block_size); ublasW.clear(); ublasW.resize(A.size1(), block_size);
439 UblasMatrixType ublasY(A.size1(), block_size); ublasY.clear(); ublasY.resize(A.size1(), block_size);
441 UblasMatrixType ublasA(A.size1(), A.size1());
443 MatrixType vclW(ublasW.size1(), ublasW.size2());
444 MatrixType vclY(ublasY.size1(), ublasY.size2());
448 for (
vcl_size_t j = 0; j < std::min(A.size1(), A.size2()); j += block_size)
450 vcl_size_t effective_block_size = std::min(std::min(A.size1(), A.size2()), j+block_size) - j;
461 for (
vcl_size_t k = 0; k < effective_block_size; ++k)
465 for (
vcl_size_t l = k; l < effective_block_size; ++l)
474 ublasY.clear(); ublasY.resize(A.size1(), block_size);
475 for (
vcl_size_t k = 0; k < effective_block_size; ++k)
492 ublasW.clear(); ublasW.resize(A.size1(), block_size);
493 ublasW(j, 0) = -betas[j];
503 for (
vcl_size_t k = 1; k < effective_block_size; ++k)
519 z = - betas[j+k] * (v_k +
prod(W_old, YT_prod_v));
540 if (A.size2() > j + effective_block_size)
563 template <
typename MatrixType,
typename VectorType>
564 void recoverQ(MatrixType
const & A, VectorType
const & betas, MatrixType & Q, MatrixType & R)
566 typedef typename MatrixType::value_type ScalarType;
568 std::vector<ScalarType> v(A.size1());
576 vcl_size_t i_max = std::min(R.size1(), R.size2());
587 vcl_size_t j_max = std::min(A.size1(), A.size2());
592 for (
vcl_size_t i=col_index+1; i<A.size1(); ++i)
593 v[i] = A(i, col_index);
595 if (betas[col_index] != 0)
607 template <
typename MatrixType,
typename VectorType1,
typename VectorType2>
615 for (
vcl_size_t col_index=0; col_index<std::min(A.size1(), A.size2()); ++col_index)
617 ScalarType v_in_b = b[col_index];
618 for (
vcl_size_t i=col_index+1; i<A.size1(); ++i)
619 v_in_b += A(i, col_index) * b[i];
621 b[col_index] -= betas[col_index] * v_in_b;
622 for (
vcl_size_t i=col_index+1; i<A.size1(); ++i)
623 b[i] -= betas[col_index] * A(i, col_index) * v_in_b;
627 template <
typename T,
typename F,
unsigned int ALIGNMENT,
typename VectorType1,
unsigned int A2>
630 boost::numeric::ublas::matrix<T> ublas_A(A.
size1(), A.
size2());
633 std::vector<T> stl_b(b.
size());
646 template<
typename T,
typename F,
unsigned int ALIGNMENT>
657 template<
typename MatrixType>
size_type size1() const
Returns the number of rows.
Definition: matrix.hpp:625
std::size_t vcl_size_t
Definition: forwards.h:58
std::vector< T > inplace_qr(viennacl::matrix< T, F, ALIGNMENT > &A, vcl_size_t block_size=16)
Overload of inplace-QR factorization of a ViennaCL matrix A.
Definition: qr.hpp:647
size_type size() const
Returns the length of the vector (cf. std::vector)
Definition: vector.hpp:859
viennacl::enable_if< viennacl::is_any_sparse_matrix< M1 >::value, matrix_expression< const M1, const M1, op_trans > >::type trans(const M1 &mat)
Returns an expression template class representing a transposed matrix.
Definition: sparse_matrix_operations.hpp:330
std::vector< typename viennacl::result_of::cpu_value_type< typename MatrixType::value_type >::type > inplace_qr_hybrid(MatrixType &A, vcl_size_t block_size=16)
Implementation of a hybrid QR factorization using uBLAS on the CPU and ViennaCL for GPUs (or multi-co...
Definition: qr.hpp:426
Generic interface for matrix-vector and matrix-matrix products. See viennacl/linalg/vector_operations...
std::vector< typename MatrixType::value_type > inplace_qr_ublas(MatrixType &A, vcl_size_t block_size=32)
Implementation of inplace-QR factorization for a general Boost.uBLAS compatible matrix A...
Definition: qr.hpp:228
Implementation of the dense matrix class.
matrix_range< MatrixType > project(MatrixType &A, viennacl::range const &r1, viennacl::range const &r2)
Definition: matrix_proxy.hpp:251
void prod(const T1 &A, bool transposed_A, const T2 &B, bool transposed_B, T3 &C, ScalarType alpha, ScalarType beta)
Definition: matrix_operations.hpp:2305
A dense matrix class.
Definition: forwards.h:293
void write_householder_to_A_viennacl(MatrixType &A, VectorType const &v, vcl_size_t j)
Definition: qr.hpp:211
size_type size2() const
Returns the number of columns.
Definition: matrix.hpp:627
MatrixType::value_type setup_householder_vector_ublas(MatrixType const &A, VectorType &v, MatrixType &matrix_1x1, vcl_size_t j)
Definition: qr.hpp:53
void write_householder_to_A(MatrixType &A, VectorType const &v, vcl_size_t j)
Definition: qr.hpp:194
VectorT prod(std::vector< std::vector< T, A1 >, A2 > const &matrix, VectorT const &vector)
Definition: prod.hpp:91
void householder_reflect_ublas(MatrixType &A, VectorType &v, MatrixType &matrix_1x1, ScalarType beta, vcl_size_t j, vcl_size_t k)
Definition: qr.hpp:147
viennacl::result_of::cpu_value_type< typename MatrixType::value_type >::type setup_householder_vector_viennacl(MatrixType const &A, VectorType &v, MatrixType &matrix_1x1, vcl_size_t j)
Definition: qr.hpp:93
void householder_reflect(MatrixType &A, VectorType &v, ScalarType beta, vcl_size_t j, vcl_size_t k)
Definition: qr.hpp:134
std::vector< typename viennacl::result_of::cpu_value_type< typename MatrixType::value_type >::type > inplace_qr_viennacl(MatrixType &A, vcl_size_t block_size=16)
Implementation of a OpenCL-only QR factorization for GPUs (or multi-core CPU). DEPRECATED! Use only i...
Definition: qr.hpp:322
void copy(std::vector< SCALARTYPE > &cpu_vec, circulant_matrix< SCALARTYPE, ALIGNMENT > &gpu_mat)
Copies a circulant matrix from the std::vector to the OpenCL device (either GPU or multi-core CPU) ...
Definition: circulant_matrix.hpp:150
void recoverQ(MatrixType const &A, VectorType const &betas, MatrixType &Q, MatrixType &R)
Definition: qr.hpp:564
A vector class representing a linear memory sequence on the GPU. Inspired by boost::numeric::ublas::v...
Definition: forwards.h:208
void write_householder_to_A_ublas(MatrixType &A, VectorType const &v, vcl_size_t j)
Definition: qr.hpp:201
Proxy classes for matrices.
T::ERROR_CANNOT_DEDUCE_CPU_SCALAR_TYPE_FOR_T type
Definition: result_of.hpp:276
void householder_reflect_viennacl(MatrixType &A, VectorType &v, MatrixType &matrix_1x1, ScalarType beta, vcl_size_t j, vcl_size_t k)
Definition: qr.hpp:161
A range class that refers to an interval [start, stop), where 'start' is included, and 'stop' is excluded.
Definition: forwards.h:339
Implementation of a range object for use with proxy objects.
basic_range range
Definition: forwards.h:339
Class for representing non-strided submatrices of a bigger matrix A.
Definition: forwards.h:355
void inplace_qr_apply_trans_Q(MatrixType const &A, VectorType1 const &betas, VectorType2 &b)
Computes Q^T b, where Q is an implicit orthogonal matrix defined via its Householder reflectors store...
Definition: qr.hpp:608