ViennaCL - The Vienna Computing Library  1.7.0
Free open-source GPU-accelerated linear algebra and solver library.
qr.hpp
Go to the documentation of this file.
1 #ifndef VIENNACL_LINALG_QR_HPP
2 #define VIENNACL_LINALG_QR_HPP
3 
4 /* =========================================================================
5  Copyright (c) 2010-2015, Institute for Microelectronics,
6  Institute for Analysis and Scientific Computing,
7  TU Wien.
8  Portions of this software are copyright by UChicago Argonne, LLC.
9 
10  -----------------
11  ViennaCL - The Vienna Computing Library
12  -----------------
13 
14  Project Head: Karl Rupp rupp@iue.tuwien.ac.at
15 
16  (A list of authors and contributors can be found in the manual)
17 
18  License: MIT (X11), see file LICENSE in the base directory
19 ============================================================================= */
20 
25 #include <utility>
26 #include <iostream>
27 #include <fstream>
28 #include <string>
29 #include <algorithm>
30 #include <vector>
31 #include <math.h>
32 #include <cmath>
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"
39 
40 #include "viennacl/matrix.hpp"
41 #include "viennacl/matrix_proxy.hpp"
42 #include "viennacl/linalg/prod.hpp"
43 #include "viennacl/range.hpp"
44 
45 namespace viennacl
46 {
47  namespace linalg
48  {
49  namespace detail
50  {
51 
52  template<typename MatrixType, typename VectorType>
53  typename MatrixType::value_type setup_householder_vector_ublas(MatrixType const & A, VectorType & v, MatrixType & matrix_1x1, vcl_size_t j)
54  {
55  using boost::numeric::ublas::range;
56  using boost::numeric::ublas::project;
57 
58  typedef typename MatrixType::value_type ScalarType;
59 
60  //compute norm of column below diagonal:
61  matrix_1x1 = boost::numeric::ublas::prod( trans(project(A, range(j+1, A.size1()), range(j, j+1))),
62  project(A, range(j+1, A.size1()), range(j, j+1))
63  );
64  ScalarType sigma = matrix_1x1(0,0);
65  ScalarType beta = 0;
66  ScalarType A_jj = A(j,j);
67 
68  assert( sigma >= 0.0 && bool("sigma must be non-negative!"));
69 
70  //get v from A:
71  v(j,0) = 1.0;
72  project(v, range(j+1, A.size1()), range(0,1)) = project(A, range(j+1, A.size1()), range(j,j+1));
73 
74  if (sigma <= 0)
75  return 0;
76  else
77  {
78  ScalarType mu = std::sqrt(sigma + A_jj*A_jj);
79 
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);
82 
83  //divide v by its diagonal element v[j]
84  project(v, range(j+1, A.size1()), range(0,1)) /= v1;
85  }
86 
87  return beta;
88  }
89 
90 
91  template<typename MatrixType, typename VectorType>
92  typename viennacl::result_of::cpu_value_type< typename MatrixType::value_type >::type
93  setup_householder_vector_viennacl(MatrixType const & A, VectorType & v, MatrixType & matrix_1x1, vcl_size_t j)
94  {
95  using viennacl::range;
96  using viennacl::project;
97 
98  typedef typename viennacl::result_of::cpu_value_type< typename MatrixType::value_type >::type ScalarType;
99 
100  //compute norm of column below diagonal:
101  matrix_1x1 = viennacl::linalg::prod( trans(project(A, range(j+1, A.size1()), range(j, j+1))),
102  project(A, range(j+1, A.size1()), range(j, j+1))
103  );
104  ScalarType sigma = matrix_1x1(0,0);
105  ScalarType beta = 0;
106  ScalarType A_jj = A(j,j);
107 
108  assert( sigma >= 0.0 && bool("sigma must be non-negative!"));
109 
110  //get v from A:
111  v(j,0) = 1.0;
112  project(v, range(j+1, A.size1()), range(0,1)) = project(A, range(j+1, A.size1()), range(j,j+1));
113 
114  if (sigma == 0)
115  return 0;
116  else
117  {
118  ScalarType mu = std::sqrt(sigma + A_jj*A_jj);
119 
120  ScalarType v1 = (A_jj <= 0) ? (A_jj - mu) : (-sigma / (A_jj + mu));
121 
122  beta = 2.0 * v1 * v1 / (sigma + v1 * v1);
123 
124  //divide v by its diagonal element v[j]
125  project(v, range(j+1, A.size1()), range(0,1)) /= v1;
126  }
127 
128  return beta;
129  }
130 
131 
132  // Apply (I - beta v v^T) to the k-th column of A, where v is the reflector starting at j-th row/column
133  template<typename MatrixType, typename VectorType, typename ScalarType>
134  void householder_reflect(MatrixType & A, VectorType & v, ScalarType beta, vcl_size_t j, vcl_size_t k)
135  {
136  ScalarType v_in_col = A(j,k);
137  for (vcl_size_t i=j+1; i<A.size1(); ++i)
138  v_in_col += v[i] * A(i,k);
139 
140  //assert(v[j] == 1.0);
141 
142  for (vcl_size_t i=j; i<A.size1(); ++i)
143  A(i,k) -= beta * v_in_col * v[i];
144  }
145 
146  template<typename MatrixType, typename VectorType, typename ScalarType>
147  void householder_reflect_ublas(MatrixType & A, VectorType & v, MatrixType & matrix_1x1, ScalarType beta, vcl_size_t j, vcl_size_t k)
148  {
149  using boost::numeric::ublas::range;
150  using boost::numeric::ublas::project;
151 
152  ScalarType v_in_col = A(j,k);
153  matrix_1x1 = boost::numeric::ublas::prod(trans(project(v, range(j+1, A.size1()), range(0, 1))),
154  project(A, range(j+1, A.size1()), range(k,k+1)));
155  v_in_col += matrix_1x1(0,0);
156 
157  project(A, range(j, A.size1()), range(k, k+1)) -= (beta * v_in_col) * project(v, range(j, A.size1()), range(0, 1));
158  }
159 
160  template<typename MatrixType, typename VectorType, typename ScalarType>
161  void householder_reflect_viennacl(MatrixType & A, VectorType & v, MatrixType & matrix_1x1, ScalarType beta, vcl_size_t j, vcl_size_t k)
162  {
163  using viennacl::range;
164  using viennacl::project;
165 
166  ScalarType v_in_col = A(j,k);
167 
168  matrix_1x1 = viennacl::linalg::prod(trans(project(v, range(j+1, A.size1()), range(0, 1))),
169  project(A, range(j+1, A.size1()), range(k,k+1)));
170  v_in_col += matrix_1x1(0,0);
171 
172  if ( beta * v_in_col != 0.0)
173  {
174  VectorType temp = project(v, range(j, A.size1()), range(0, 1));
175  project(v, range(j, A.size1()), range(0, 1)) *= (beta * v_in_col);
176  project(A, range(j, A.size1()), range(k, k+1)) -= project(v, range(j, A.size1()), range(0, 1));
177  project(v, range(j, A.size1()), range(0, 1)) = temp;
178  }
179  }
180 
181 
182  // Apply (I - beta v v^T) to A, where v is the reflector starting at j-th row/column
183  template<typename MatrixType, typename VectorType, typename ScalarType>
184  void householder_reflect(MatrixType & A, VectorType & v, ScalarType beta, vcl_size_t j)
185  {
186  vcl_size_t column_end = A.size2();
187 
188  for (vcl_size_t k=j; k<column_end; ++k) //over columns
189  householder_reflect(A, v, beta, j, k);
190  }
191 
192 
193  template<typename MatrixType, typename VectorType>
194  void write_householder_to_A(MatrixType & A, VectorType const & v, vcl_size_t j)
195  {
196  for (vcl_size_t i=j+1; i<A.size1(); ++i)
197  A(i,j) = v[i];
198  }
199 
200  template<typename MatrixType, typename VectorType>
201  void write_householder_to_A_ublas(MatrixType & A, VectorType const & v, vcl_size_t j)
202  {
203  using boost::numeric::ublas::range;
204  using boost::numeric::ublas::project;
205 
206  //VectorType temp = project(v, range(j+1, A.size1()));
207  project( A, range(j+1, A.size1()), range(j, j+1) ) = project(v, range(j+1, A.size1()), range(0, 1) );;
208  }
209 
210  template<typename MatrixType, typename VectorType>
211  void write_householder_to_A_viennacl(MatrixType & A, VectorType const & v, vcl_size_t j)
212  {
213  using viennacl::range;
214  using viennacl::project;
215 
216  //VectorType temp = project(v, range(j+1, A.size1()));
217  project( A, range(j+1, A.size1()), range(j, j+1) ) = project(v, range(j+1, A.size1()), range(0, 1) );;
218  }
219 
220 
221 
227  template<typename MatrixType>
228  std::vector<typename MatrixType::value_type> inplace_qr_ublas(MatrixType & A, vcl_size_t block_size = 32)
229  {
230  typedef typename MatrixType::value_type ScalarType;
231  typedef boost::numeric::ublas::matrix_range<MatrixType> MatrixRange;
232 
233  using boost::numeric::ublas::range;
234  using boost::numeric::ublas::project;
235 
236  std::vector<ScalarType> betas(A.size2());
237  MatrixType v(A.size1(), 1);
238  MatrixType matrix_1x1(1,1);
239 
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);
242 
243  //run over A in a block-wise manner:
244  for (vcl_size_t j = 0; j < std::min(A.size1(), A.size2()); j += block_size)
245  {
246  vcl_size_t effective_block_size = std::min(std::min(A.size1(), A.size2()), j+block_size) - j;
247 
248  //determine Householder vectors:
249  for (vcl_size_t k = 0; k < effective_block_size; ++k)
250  {
251  betas[j+k] = detail::setup_householder_vector_ublas(A, v, matrix_1x1, j+k);
252 
253  for (vcl_size_t l = k; l < effective_block_size; ++l)
254  detail::householder_reflect_ublas(A, v, matrix_1x1, betas[j+k], j+k, j+l);
255 
256  detail::write_householder_to_A_ublas(A, v, j+k);
257  }
258 
259  //
260  // Setup Y:
261  //
262  Y.clear(); Y.resize(A.size1(), block_size);
263  for (vcl_size_t k = 0; k < effective_block_size; ++k)
264  {
265  //write Householder to Y:
266  Y(j+k,k) = 1.0;
267  project(Y, range(j+k+1, A.size1()), range(k, k+1)) = project(A, range(j+k+1, A.size1()), range(j+k, j+k+1));
268  }
269 
270  //
271  // Setup W:
272  //
273 
274  //first vector:
275  W.clear(); W.resize(A.size1(), block_size);
276  W(j, 0) = -betas[j];
277  project(W, range(j+1, A.size1()), range(0, 1)) = -betas[j] * project(A, range(j+1, A.size1()), range(j, j+1));
278 
279 
280  //k-th column of W is given by -beta * (Id + W*Y^T) v_k, where W and Y have k-1 columns
281  for (vcl_size_t k = 1; k < effective_block_size; ++k)
282  {
283  MatrixRange Y_old = project(Y, range(j, A.size1()), range(0, k));
284  MatrixRange v_k = project(Y, range(j, A.size1()), range(k, k+1));
285  MatrixRange W_old = project(W, range(j, A.size1()), range(0, k));
286  MatrixRange z = project(W, range(j, A.size1()), range(k, k+1));
287 
288  MatrixType YT_prod_v = boost::numeric::ublas::prod(boost::numeric::ublas::trans(Y_old), v_k);
289  z = - betas[j+k] * (v_k + prod(W_old, YT_prod_v));
290  }
291 
292  //
293  //apply (I+WY^T)^T = I + Y W^T to the remaining columns of A:
294  //
295 
296  if (A.size2() - j - effective_block_size > 0)
297  {
298 
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));
301  MatrixType temp = boost::numeric::ublas::prod(trans(W_part), A_part);
302 
303  A_part += prod(project(Y, range(j, A.size1()), range(0, effective_block_size)),
304  temp);
305  }
306  }
307 
308  return betas;
309  }
310 
311 
320  template<typename MatrixType>
321  std::vector< typename viennacl::result_of::cpu_value_type< typename MatrixType::value_type >::type >
322  inplace_qr_viennacl(MatrixType & A, vcl_size_t block_size = 16)
323  {
324  typedef typename viennacl::result_of::cpu_value_type< typename MatrixType::value_type >::type ScalarType;
325  typedef viennacl::matrix_range<MatrixType> MatrixRange;
326 
327  using viennacl::range;
328  using viennacl::project;
329 
330  std::vector<ScalarType> betas(A.size2());
331  MatrixType v(A.size1(), 1);
332  MatrixType matrix_1x1(1,1);
333 
334  MatrixType Y(A.size1(), block_size); Y.clear();
335  MatrixType W(A.size1(), block_size); W.clear();
336 
337  MatrixType YT_prod_v(block_size, 1);
338  MatrixType z(A.size1(), 1);
339 
340  //run over A in a block-wise manner:
341  for (vcl_size_t j = 0; j < std::min(A.size1(), A.size2()); j += block_size)
342  {
343  vcl_size_t effective_block_size = std::min(std::min(A.size1(), A.size2()), j+block_size) - j;
344 
345  //determine Householder vectors:
346  for (vcl_size_t k = 0; k < effective_block_size; ++k)
347  {
348  betas[j+k] = detail::setup_householder_vector_viennacl(A, v, matrix_1x1, j+k);
349  for (vcl_size_t l = k; l < effective_block_size; ++l)
350  detail::householder_reflect_viennacl(A, v, matrix_1x1, betas[j+k], j+k, j+l);
351 
352  detail::write_householder_to_A_viennacl(A, v, j+k);
353  }
354 
355  //
356  // Setup Y:
357  //
358  Y.clear();
359  for (vcl_size_t k = 0; k < effective_block_size; ++k)
360  {
361  //write Householder to Y:
362  Y(j+k,k) = 1.0;
363  project(Y, range(j+k+1, A.size1()), range(k, k+1)) = project(A, range(j+k+1, A.size1()), range(j+k, j+k+1));
364  }
365 
366  //
367  // Setup W:
368  //
369 
370  //first vector:
371  W.clear();
372  W(j, 0) = -betas[j];
373  //project(W, range(j+1, A.size1()), range(0, 1)) = -betas[j] * project(A, range(j+1, A.size1()), range(j, j+1));
374  project(W, range(j+1, A.size1()), range(0, 1)) = project(A, range(j+1, A.size1()), range(j, j+1));
375  project(W, range(j+1, A.size1()), range(0, 1)) *= -betas[j];
376 
377 
378  //k-th column of W is given by -beta * (Id + W*Y^T) v_k, where W and Y have k-1 columns
379  for (vcl_size_t k = 1; k < effective_block_size; ++k)
380  {
381  MatrixRange Y_old = project(Y, range(j, A.size1()), range(0, k));
382  MatrixRange v_k = project(Y, range(j, A.size1()), range(k, k+1));
383  MatrixRange W_old = project(W, range(j, A.size1()), range(0, k));
384 
385  project(YT_prod_v, range(0, k), range(0,1)) = prod(trans(Y_old), v_k);
386  project(z, range(j, A.size1()), range(0,1)) = prod(W_old, project(YT_prod_v, range(0, k), range(0,1)));
387  project(W, range(j, A.size1()), range(k, k+1)) = project(z, range(j, A.size1()), range(0,1));
388  project(W, range(j, A.size1()), range(k, k+1)) += v_k;
389  project(W, range(j, A.size1()), range(k, k+1)) *= - betas[j+k];
390  }
391 
392  //
393  //apply (I+WY^T)^T = I + Y W^T to the remaining columns of A:
394  //
395 
396  if (A.size2() > j + effective_block_size)
397  {
398 
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);
402 
403  A_part += prod(project(Y, range(j, A.size1()), range(0, effective_block_size)),
404  temp);
405  }
406  }
407 
408  return betas;
409  }
410 
411 
412 
413 
414 
415 
416  //MatrixType is ViennaCL-matrix
424  template<typename MatrixType>
425  std::vector< typename viennacl::result_of::cpu_value_type< typename MatrixType::value_type >::type >
426  inplace_qr_hybrid(MatrixType & A, vcl_size_t block_size = 16)
427  {
428  typedef typename viennacl::result_of::cpu_value_type< typename MatrixType::value_type >::type ScalarType;
429 
430  typedef viennacl::matrix_range<MatrixType> VCLMatrixRange;
431  typedef boost::numeric::ublas::matrix<ScalarType> UblasMatrixType;
432  typedef boost::numeric::ublas::matrix_range<UblasMatrixType> UblasMatrixRange;
433 
434  std::vector<ScalarType> betas(A.size2());
435  UblasMatrixType v(A.size1(), 1);
436  UblasMatrixType matrix_1x1(1,1);
437 
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);
440 
441  UblasMatrixType ublasA(A.size1(), A.size1());
442 
443  MatrixType vclW(ublasW.size1(), ublasW.size2());
444  MatrixType vclY(ublasY.size1(), ublasY.size2());
445 
446 
447  //run over A in a block-wise manner:
448  for (vcl_size_t j = 0; j < std::min(A.size1(), A.size2()); j += block_size)
449  {
450  vcl_size_t effective_block_size = std::min(std::min(A.size1(), A.size2()), j+block_size) - j;
451  UblasMatrixRange ublasA_part = boost::numeric::ublas::project(ublasA,
452  boost::numeric::ublas::range(0, A.size1()),
453  boost::numeric::ublas::range(j, j + effective_block_size));
454  viennacl::copy(viennacl::project(A,
455  viennacl::range(0, A.size1()),
456  viennacl::range(j, j+effective_block_size)),
457  ublasA_part
458  );
459 
460  //determine Householder vectors:
461  for (vcl_size_t k = 0; k < effective_block_size; ++k)
462  {
463  betas[j+k] = detail::setup_householder_vector_ublas(ublasA, v, matrix_1x1, j+k);
464 
465  for (vcl_size_t l = k; l < effective_block_size; ++l)
466  detail::householder_reflect_ublas(ublasA, v, matrix_1x1, betas[j+k], j+k, j+l);
467 
468  detail::write_householder_to_A_ublas(ublasA, v, j+k);
469  }
470 
471  //
472  // Setup Y:
473  //
474  ublasY.clear(); ublasY.resize(A.size1(), block_size);
475  for (vcl_size_t k = 0; k < effective_block_size; ++k)
476  {
477  //write Householder to Y:
478  ublasY(j+k,k) = 1.0;
479  boost::numeric::ublas::project(ublasY,
480  boost::numeric::ublas::range(j+k+1, A.size1()),
481  boost::numeric::ublas::range(k, k+1))
482  = boost::numeric::ublas::project(ublasA,
483  boost::numeric::ublas::range(j+k+1, A.size1()),
484  boost::numeric::ublas::range(j+k, j+k+1));
485  }
486 
487  //
488  // Setup W:
489  //
490 
491  //first vector:
492  ublasW.clear(); ublasW.resize(A.size1(), block_size);
493  ublasW(j, 0) = -betas[j];
494  boost::numeric::ublas::project(ublasW,
495  boost::numeric::ublas::range(j+1, A.size1()),
496  boost::numeric::ublas::range(0, 1))
497  = -betas[j] * boost::numeric::ublas::project(ublasA,
498  boost::numeric::ublas::range(j+1, A.size1()),
499  boost::numeric::ublas::range(j, j+1));
500 
501 
502  //k-th column of W is given by -beta * (Id + W*Y^T) v_k, where W and Y have k-1 columns
503  for (vcl_size_t k = 1; k < effective_block_size; ++k)
504  {
505  UblasMatrixRange Y_old = boost::numeric::ublas::project(ublasY,
506  boost::numeric::ublas::range(j, A.size1()),
507  boost::numeric::ublas::range(0, k));
508  UblasMatrixRange v_k = boost::numeric::ublas::project(ublasY,
509  boost::numeric::ublas::range(j, A.size1()),
510  boost::numeric::ublas::range(k, k+1));
511  UblasMatrixRange W_old = boost::numeric::ublas::project(ublasW,
512  boost::numeric::ublas::range(j, A.size1()),
513  boost::numeric::ublas::range(0, k));
514  UblasMatrixRange z = boost::numeric::ublas::project(ublasW,
515  boost::numeric::ublas::range(j, A.size1()),
516  boost::numeric::ublas::range(k, k+1));
517 
518  UblasMatrixType YT_prod_v = boost::numeric::ublas::prod(boost::numeric::ublas::trans(Y_old), v_k);
519  z = - betas[j+k] * (v_k + prod(W_old, YT_prod_v));
520  }
521 
522 
523 
524  //
525  //apply (I+WY^T)^T = I + Y W^T to the remaining columns of A:
526  //
527 
528  VCLMatrixRange A_part = viennacl::project(A,
529  viennacl::range(0, A.size1()),
530  viennacl::range(j, j+effective_block_size));
531 
532  viennacl::copy(boost::numeric::ublas::project(ublasA,
533  boost::numeric::ublas::range(0, A.size1()),
534  boost::numeric::ublas::range(j, j+effective_block_size)),
535  A_part);
536 
537  viennacl::copy(ublasW, vclW);
538  viennacl::copy(ublasY, vclY);
539 
540  if (A.size2() > j + effective_block_size)
541  {
542 
543  VCLMatrixRange A_part2(A, viennacl::range(j, A.size1()), viennacl::range(j+effective_block_size, A.size2()));
544  VCLMatrixRange W_part(vclW, viennacl::range(j, A.size1()), viennacl::range(0, effective_block_size));
545  MatrixType temp = viennacl::linalg::prod(trans(W_part), A_part2);
546 
547  A_part2 += viennacl::linalg::prod(viennacl::project(vclY, viennacl::range(j, A.size1()), viennacl::range(0, effective_block_size)),
548  temp);
549  }
550  }
551 
552  return betas;
553  }
554 
555 
556 
557  } //namespace detail
558 
559 
560 
561 
562  //takes an inplace QR matrix A and generates Q and R explicitly
563  template<typename MatrixType, typename VectorType>
564  void recoverQ(MatrixType const & A, VectorType const & betas, MatrixType & Q, MatrixType & R)
565  {
566  typedef typename MatrixType::value_type ScalarType;
567 
568  std::vector<ScalarType> v(A.size1());
569 
570  Q.clear();
571  R.clear();
572 
573  //
574  // Recover R from upper-triangular part of A:
575  //
576  vcl_size_t i_max = std::min(R.size1(), R.size2());
577  for (vcl_size_t i=0; i<i_max; ++i)
578  for (vcl_size_t j=i; j<R.size2(); ++j)
579  R(i,j) = A(i,j);
580 
581  //
582  // Recover Q by applying all the Householder reflectors to the identity matrix:
583  //
584  for (vcl_size_t i=0; i<Q.size1(); ++i)
585  Q(i,i) = 1.0;
586 
587  vcl_size_t j_max = std::min(A.size1(), A.size2());
588  for (vcl_size_t j=0; j<j_max; ++j)
589  {
590  vcl_size_t col_index = j_max - j - 1;
591  v[col_index] = 1.0;
592  for (vcl_size_t i=col_index+1; i<A.size1(); ++i)
593  v[i] = A(i, col_index);
594 
595  if (betas[col_index] > 0 || betas[col_index] < 0)
596  detail::householder_reflect(Q, v, betas[col_index], col_index);
597  }
598  }
599 
600 
607  template<typename MatrixType, typename VectorType1, typename VectorType2>
608  void inplace_qr_apply_trans_Q(MatrixType const & A, VectorType1 const & betas, VectorType2 & b)
609  {
610  typedef typename viennacl::result_of::cpu_value_type<typename MatrixType::value_type>::type ScalarType;
611 
612  //
613  // Apply Q^T = (I - beta_m v_m v_m^T) \times ... \times (I - beta_0 v_0 v_0^T) by applying all the Householder reflectors to b:
614  //
615  for (vcl_size_t col_index=0; col_index<std::min(A.size1(), A.size2()); ++col_index)
616  {
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];
620 
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;
624  }
625  }
626 
627  template<typename T, typename F, unsigned int ALIGNMENT, typename VectorType1, unsigned int A2>
628  void inplace_qr_apply_trans_Q(viennacl::matrix<T, F, ALIGNMENT> const & A, VectorType1 const & betas, viennacl::vector<T, A2> & b)
629  {
630  boost::numeric::ublas::matrix<T> ublas_A(A.size1(), A.size2());
631  viennacl::copy(A, ublas_A);
632 
633  std::vector<T> stl_b(b.size());
634  viennacl::copy(b, stl_b);
635 
636  inplace_qr_apply_trans_Q(ublas_A, betas, stl_b);
637 
638  viennacl::copy(stl_b, b);
639  }
640 
646  template<typename T, typename F, unsigned int ALIGNMENT>
647  std::vector<T> inplace_qr(viennacl::matrix<T, F, ALIGNMENT> & A, vcl_size_t block_size = 16)
648  {
649  return detail::inplace_qr_hybrid(A, block_size);
650  }
651 
657  template<typename MatrixType>
658  std::vector<typename MatrixType::value_type> inplace_qr(MatrixType & A, vcl_size_t block_size = 16)
659  {
660  return detail::inplace_qr_ublas(A, block_size);
661  }
662 
663 
664 
665  } //linalg
666 } //viennacl
667 
668 
669 #endif
viennacl::linalg::inplace_qr
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
viennacl::linalg::detail::inplace_qr_hybrid
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
viennacl::linalg::cuda::trans
void trans(matrix_expression< const matrix_base< NumericT, SizeT, DistanceT >, const matrix_base< NumericT, SizeT, DistanceT >, op_trans > const &proxy, matrix_base< NumericT > &temp_trans)
Definition: matrix_operations.hpp:94
prod.hpp
Generic interface for matrix-vector and matrix-matrix products. See viennacl/linalg/vector_operations...
viennacl::linalg::detail::inplace_qr_ublas
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
matrix.hpp
Implementation of the dense matrix class.
viennacl::matrix
A dense matrix class.
Definition: forwards.h:375
viennacl::linalg::detail::write_householder_to_A_viennacl
void write_householder_to_A_viennacl(MatrixType &A, VectorType const &v, vcl_size_t j)
Definition: qr.hpp:211
viennacl::linalg::detail::setup_householder_vector_ublas
MatrixType::value_type setup_householder_vector_ublas(MatrixType const &A, VectorType &v, MatrixType &matrix_1x1, vcl_size_t j)
Definition: qr.hpp:53
viennacl::linalg::detail::write_householder_to_A
void write_householder_to_A(MatrixType &A, VectorType const &v, vcl_size_t j)
Definition: qr.hpp:194
viennacl::range
basic_range range
Definition: forwards.h:424
viennacl
Main namespace in ViennaCL. Holds all the basic types such as vector, matrix, etc. and defines operations upon them.
Definition: cpu_ram.hpp:34
v1
viennacl::vector< float > v1
Definition: global_variables.cpp:60
viennacl::linalg::prod
VectorT prod(std::vector< std::vector< T, A1 >, A2 > const &matrix, VectorT const &vector)
Definition: prod.hpp:102
viennacl::linalg::detail::householder_reflect_ublas
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::linalg::detail::setup_householder_vector_viennacl
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
viennacl::linalg::detail::householder_reflect
void householder_reflect(MatrixType &A, VectorType &v, ScalarType beta, vcl_size_t j, vcl_size_t k)
Definition: qr.hpp:134
detail
Definition: blas3.hpp:36
viennacl::linalg::detail::inplace_qr_viennacl
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
viennacl::linalg::recoverQ
void recoverQ(MatrixType const &A, VectorType const &betas, MatrixType &Q, MatrixType &R)
Definition: qr.hpp:564
viennacl::project
matrix_range< MatrixType > project(MatrixType const &A, viennacl::range const &r1, viennacl::range const &r2)
Definition: matrix_proxy.hpp:293
viennacl::vcl_size_t
std::size_t vcl_size_t
Definition: forwards.h:75
viennacl::matrix_base< NumericT >::size2
size_type size2() const
Returns the number of columns.
Definition: matrix_def.hpp:226
viennacl::linalg::trans
void trans(const matrix_expression< const matrix_base< NumericT, SizeT, DistanceT >, const matrix_base< NumericT, SizeT, DistanceT >, op_trans > &proxy, matrix_base< NumericT > &temp_trans)
Definition: matrix_operations.hpp:83
viennacl::result_of::cpu_value_type::type
T::ERROR_CANNOT_DEDUCE_CPU_SCALAR_TYPE_FOR_T type
Definition: result_of.hpp:271
viennacl::matrix_base< NumericT >::size1
size_type size1() const
Returns the number of rows.
Definition: matrix_def.hpp:224
viennacl::linalg::detail::write_householder_to_A_ublas
void write_householder_to_A_ublas(MatrixType &A, VectorType const &v, vcl_size_t j)
Definition: qr.hpp:201
matrix_proxy.hpp
Proxy classes for matrices.
prod
void prod(std::vector< std::map< IndexT, NumericT > > const &stl_A, std::vector< std::map< IndexT, NumericT > > const &stl_B, std::vector< std::map< IndexT, NumericT > > &stl_C)
Definition: sparse_prod.cpp:114
viennacl::copy
void copy(std::vector< NumericT > &cpu_vec, circulant_matrix< NumericT, AlignmentV > &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
viennacl::linalg::detail::householder_reflect_viennacl
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
viennacl::vector_base< NumericT >::size
size_type size() const
Returns the length of the vector (cf. std::vector)
Definition: vector_def.hpp:118
viennacl::basic_range
A range class that refers to an interval [start, stop), where 'start' is included, and 'stop' is excluded.
Definition: forwards.h:424
ScalarType
float ScalarType
Definition: fft_1d.cpp:42
range.hpp
Implementation of a range object for use with proxy objects.
viennacl::matrix_range
Class for representing non-strided submatrices of a bigger matrix A.
Definition: forwards.h:440
viennacl::linalg::detail::min
T min(const T &lhs, const T &rhs)
Minimum.
Definition: util.hpp:45
viennacl::linalg::inplace_qr_apply_trans_Q
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