1 #ifndef VIENNACL_LINALG_QR_HPP
2 #define VIENNACL_LINALG_QR_HPP
32 #include "boost/numeric/ublas/vector.hpp"
33 #include "boost/numeric/ublas/matrix.hpp"
34 #include "boost/numeric/ublas/matrix_proxy.hpp"
35 #include "boost/numeric/ublas/io.hpp"
36 #include "boost/numeric/ublas/matrix_expression.hpp"
47 template <
typename MatrixType,
typename VectorType>
50 typedef typename MatrixType::value_type ScalarType;
55 for (
size_t k = j+1; k<A.size1(); ++k)
56 sigma += A(k, j) * A(k, j);
59 for (
size_t k = j+1; k<A.size1(); ++k)
66 ScalarType mu = sqrt(sigma + A(j,j)*A(j,j));
74 v1 = -sigma / (A(j,j) + mu);
76 beta = 2.0 * v1 * v1 / (sigma + v1 * v1);
80 for (
size_t k = j+1; k<A.size1(); ++k)
88 template <
typename MatrixType,
typename VectorType,
typename ScalarType>
91 ScalarType v_in_col = A(j,k);
92 for (
size_t i=j+1; i<A.size1(); ++i)
93 v_in_col += v[i] * A(i,k);
95 for (
size_t i=j; i<A.size1(); ++i)
96 A(i,k) -= beta * v_in_col * v[i];
100 template <
typename MatrixType,
typename VectorType,
typename ScalarType>
103 size_t column_end = A.size2();
105 for (
size_t k=j; k<column_end; ++k)
110 template <
typename MatrixType,
typename VectorType>
113 for (
size_t i=j+1; i<A.size1(); ++i)
119 template <
typename MatrixType,
typename VectorType>
120 void recoverQ(MatrixType
const & A, VectorType
const & betas, MatrixType & Q, MatrixType & R)
122 typedef typename MatrixType::value_type ScalarType;
124 std::vector<ScalarType> v(A.size1());
132 size_t i_max = std::min(R.size1(), R.size2());
133 for (
size_t i=0; i<i_max; ++i)
134 for (
size_t j=i; j<R.size2(); ++j)
140 for (
size_t i=0; i<Q.size1(); ++i)
143 size_t j_max = std::min(A.size1(), A.size2());
144 for (
size_t j=0; j<j_max; ++j)
146 size_t col_index = j_max - j - 1;
148 for (
size_t i=col_index+1; i<A.size1(); ++i)
149 v[i] = A(i, col_index);
157 if (betas[col_index] != 0)
186 range(
size_t start,
size_t end) : start_(start), end_(end) {}
188 size_t lower()
const {
return start_; }
189 size_t upper()
const {
return end_; }
196 template <
typename MatrixType>
204 range col_range) : mat_(mat), row_range_(row_range), col_range_(col_range) {}
208 assert(row <
size1());
209 assert(col <
size2());
210 return mat_(row + row_range_.
lower(), col + col_range_.
lower());
224 template <
typename MatrixTypeA,
typename MatrixTypeB,
typename MatrixTypeC>
225 void prod_AA(MatrixTypeA
const & A, MatrixTypeB
const & B, MatrixTypeC & C)
227 assert(C.size1() == A.size1());
228 assert(A.size2() == B.size1());
229 assert(B.size2() == C.size2());
231 typedef typename MatrixTypeC::value_type ScalarType;
233 for (
size_t i=0; i<C.size1(); ++i)
235 for (
size_t j=0; j<C.size2(); ++j)
238 for (
size_t k=0; k<A.size2(); ++k)
239 val += A(i, k) * B(k, j);
246 template <
typename MatrixTypeA,
typename MatrixTypeB,
typename MatrixTypeC>
247 void prod_TA(MatrixTypeA
const & A, MatrixTypeB
const & B, MatrixTypeC & C)
249 assert(C.size1() == A.size2());
250 assert(A.size1() == B.size1());
251 assert(B.size2() == C.size2());
253 typedef typename MatrixTypeC::value_type ScalarType;
255 for (
size_t i=0; i<C.size1(); ++i)
257 for (
size_t j=0; j<C.size2(); ++j)
260 for (
size_t k=0; k<A.size1(); ++k)
261 val += A(k, i) * B(k, j);
269 template<
typename MatrixType>
270 std::vector<typename MatrixType::value_type>
inplace_qr(MatrixType & A, std::size_t block_size = 32)
272 typedef typename MatrixType::value_type ScalarType;
274 if ( A.size2() % block_size != 0 )
275 std::cout <<
"ViennaCL: Warning in inplace_qr(): Matrix columns are not divisible by block_size!" << std::endl;
277 std::vector<ScalarType> betas(A.size2());
278 std::vector<ScalarType> v(A.size1());
281 MatrixType Y(A.size1(), block_size); Y.clear();
282 MatrixType W(A.size1(), block_size); W.clear();
285 for (
size_t j = 0; j < std::min(A.size1(), A.size2()); j += block_size)
288 for (
size_t k = 0; k < block_size; ++k)
291 for (
size_t l = k; l < block_size; ++l)
300 for (
size_t k = 0; k < block_size; ++k)
304 for (
size_t l=k+1; l<A.size1(); ++l)
314 for (
size_t l=j+1; l<A.size1(); ++l)
315 W(l,0) = -betas[j] * A(l, j);
318 for (
size_t k = 1; k < block_size; ++k)
321 std::vector<ScalarType> temp(k);
322 for (
size_t l=0; l<k; ++l)
323 for (
size_t n=j; n<A.size1(); ++n)
324 temp[l] += Y(n, l) * Y(n, k);
327 for (
size_t n=0; n<A.size1(); ++n)
330 for (
size_t l=0; l<k; ++l)
331 val += temp[l] * W(n, l);
332 W(n, k) = -1.0 * betas[j+k] * (Y(n, k) + val);
340 if (A.size2() - j - block_size > 0)
344 MatrixType temp(block_size, A.size2() - j - block_size);
348 boost::numeric::ublas::matrix_range<MatrixType> A_part(A, A_rows, A_cols);
355 boost::numeric::ublas::matrix_range<MatrixType> W_part(W, A_rows, W_cols);
366 boost::numeric::ublas::matrix_range<MatrixType> Y_part(Y, A_rows, Y_cols);
372 gpu_A_part +=
prod(gpu_Y, gpu_temp);
374 MatrixType A_part_back(A_part.size1(), A_part.size2());
377 A_part = A_part_back;
386 template<
typename MatrixType>
389 typedef typename MatrixType::value_type ScalarType;
391 std::vector<ScalarType> betas(A.size2());
392 std::vector<ScalarType> v(A.size1());
394 size_t block_size = 3;
395 MatrixType Y(A.size1(), block_size); Y.clear();
396 MatrixType W(A.size1(), block_size); W.clear();
399 for (
size_t j = 0; j < std::min(A.size1(), A.size2()); j += block_size)
402 for (
size_t k = 0; k < block_size; ++k)
405 for (
size_t l = k; l < block_size; ++l)
414 for (
size_t k = 0; k < block_size; ++k)
418 for (
size_t l=k+1; l<A.size1(); ++l)
428 for (
size_t l=j+1; l<A.size1(); ++l)
429 W(l,0) = -betas[j] * A(l, j);
432 for (
size_t k = 1; k < block_size; ++k)
435 std::vector<ScalarType> temp(k);
436 for (
size_t l=0; l<k; ++l)
437 for (
size_t n=j; n<A.size1(); ++n)
438 temp[l] += Y(n, l) * Y(n, k);
441 for (
size_t n=0; n<A.size1(); ++n)
444 for (
size_t l=0; l<k; ++l)
445 val += temp[l] * W(n, l);
446 W(n, k) = -1.0 * betas[j+k] * (Y(n, k) + val);
454 if (A.size2() - j - block_size > 0)
457 MatrixType temp(block_size, A.size2() - j - block_size);
461 boost::numeric::ublas::matrix_range<MatrixType> A_part(A, A_rows, A_cols);
464 boost::numeric::ublas::matrix_range<MatrixType> W_part(W, A_rows, W_cols);
471 boost::numeric::ublas::matrix_range<MatrixType> Y_part(Y, A_rows, Y_cols);
473 A_part +=
prod(Y_part, temp);
481 template<
typename MatrixType>
484 typedef typename MatrixType::value_type ScalarType;
486 std::vector<ScalarType> betas(A.size2());
487 std::vector<ScalarType> v(A.size1());
489 size_t block_size = 5;
490 MatrixType Y(A.size1(), block_size); Y.clear();
491 MatrixType W(A.size1(), block_size); W.clear();
494 for (
size_t j = 0; j < std::min(A.size1(), A.size2()); j += block_size)
497 for (
size_t k = 0; k < block_size; ++k)
500 for (
size_t l = k; l < block_size; ++l)
509 for (
size_t k = 0; k < block_size; ++k)
513 for (
size_t l=k+1; l<A.size1(); ++l)
523 for (
size_t l=j+1; l<A.size1(); ++l)
524 W(l,0) = -betas[j] * A(l, j);
527 for (
size_t k = 1; k < block_size; ++k)
530 std::vector<ScalarType> temp(k);
531 for (
size_t l=0; l<k; ++l)
532 for (
size_t n=j; n<A.size1(); ++n)
533 temp[l] += Y(n, l) * Y(n, k);
536 for (
size_t n=0; n<A.size1(); ++n)
539 for (
size_t l=0; l<k; ++l)
540 val += temp[l] * W(n, l);
541 W(n, k) = -1.0 * betas[j+k] * (Y(n, k) + val);
549 if (A.size2() - j - block_size > 0)
552 MatrixType temp(block_size, A.size2() - j - block_size);
553 ScalarType entry = 0;
554 for (
size_t l = 0; l < temp.size2(); ++l)
556 for (
size_t k = 0; k < temp.size1(); ++k)
559 for (
size_t n = j; n < A.size1(); ++n)
560 entry += W(n, k) * A(n, j + block_size + l);
566 for (
size_t l = j+block_size; l < A.size2(); ++l)
568 for (
size_t k = j; k<A.size1(); ++k)
571 for (
size_t n=0; n<block_size; ++n)
572 val += Y(k, n) * temp(n, l-j-block_size);