Eigen  3.2.91
SparseQR.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2012-2013 Desire Nuentsa <desire.nuentsa_wakam@inria.fr>
5 // Copyright (C) 2012-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
6 //
7 // This Source Code Form is subject to the terms of the Mozilla
8 // Public License v. 2.0. If a copy of the MPL was not distributed
9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
10 
11 #ifndef EIGEN_SPARSE_QR_H
12 #define EIGEN_SPARSE_QR_H
13 
14 namespace Eigen {
15 
16 template<typename MatrixType, typename OrderingType> class SparseQR;
17 template<typename SparseQRType> struct SparseQRMatrixQReturnType;
18 template<typename SparseQRType> struct SparseQRMatrixQTransposeReturnType;
19 template<typename SparseQRType, typename Derived> struct SparseQR_QProduct;
20 namespace internal {
21  template <typename SparseQRType> struct traits<SparseQRMatrixQReturnType<SparseQRType> >
22  {
23  typedef typename SparseQRType::MatrixType ReturnType;
24  typedef typename ReturnType::StorageIndex StorageIndex;
25  typedef typename ReturnType::StorageKind StorageKind;
26  enum {
27  RowsAtCompileTime = Dynamic,
28  ColsAtCompileTime = Dynamic
29  };
30  };
31  template <typename SparseQRType> struct traits<SparseQRMatrixQTransposeReturnType<SparseQRType> >
32  {
33  typedef typename SparseQRType::MatrixType ReturnType;
34  };
35  template <typename SparseQRType, typename Derived> struct traits<SparseQR_QProduct<SparseQRType, Derived> >
36  {
37  typedef typename Derived::PlainObject ReturnType;
38  };
39 } // End namespace internal
40 
68 template<typename _MatrixType, typename _OrderingType>
69 class SparseQR : public SparseSolverBase<SparseQR<_MatrixType,_OrderingType> >
70 {
71  protected:
72  typedef SparseSolverBase<SparseQR<_MatrixType,_OrderingType> > Base;
73  using Base::m_isInitialized;
74  public:
75  using Base::_solve_impl;
76  typedef _MatrixType MatrixType;
77  typedef _OrderingType OrderingType;
78  typedef typename MatrixType::Scalar Scalar;
79  typedef typename MatrixType::RealScalar RealScalar;
80  typedef typename MatrixType::StorageIndex StorageIndex;
81  typedef SparseMatrix<Scalar,ColMajor,StorageIndex> QRMatrixType;
82  typedef Matrix<StorageIndex, Dynamic, 1> IndexVector;
83  typedef Matrix<Scalar, Dynamic, 1> ScalarVector;
84  typedef PermutationMatrix<Dynamic, Dynamic, StorageIndex> PermutationType;
85  public:
86  SparseQR () : m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false),m_isEtreeOk(false)
87  { }
88 
95  explicit SparseQR(const MatrixType& mat) : m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false),m_isEtreeOk(false)
96  {
97  compute(mat);
98  }
99 
106  void compute(const MatrixType& mat)
107  {
108  analyzePattern(mat);
109  factorize(mat);
110  }
111  void analyzePattern(const MatrixType& mat);
112  void factorize(const MatrixType& mat);
113 
116  inline Index rows() const { return m_pmat.rows(); }
117 
120  inline Index cols() const { return m_pmat.cols();}
121 
124  const QRMatrixType& matrixR() const { return m_R; }
125 
130  Index rank() const
131  {
132  eigen_assert(m_isInitialized && "The factorization should be called first, use compute()");
133  return m_nonzeropivots;
134  }
135 
154  SparseQRMatrixQReturnType<SparseQR> matrixQ() const
155  { return SparseQRMatrixQReturnType<SparseQR>(*this); }
156 
160  const PermutationType& colsPermutation() const
161  {
162  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
163  return m_outputPerm_c;
164  }
165 
169  std::string lastErrorMessage() const { return m_lastError; }
170 
172  template<typename Rhs, typename Dest>
173  bool _solve_impl(const MatrixBase<Rhs> &B, MatrixBase<Dest> &dest) const
174  {
175  eigen_assert(m_isInitialized && "The factorization should be called first, use compute()");
176  eigen_assert(this->rows() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix");
177 
178  Index rank = this->rank();
179 
180  // Compute Q^T * b;
181  typename Dest::PlainObject y, b;
182  y = this->matrixQ().transpose() * B;
183  b = y;
184 
185  // Solve with the triangular matrix R
186  y.resize((std::max<Index>)(cols(),y.rows()),y.cols());
187  y.topRows(rank) = this->matrixR().topLeftCorner(rank, rank).template triangularView<Upper>().solve(b.topRows(rank));
188  y.bottomRows(y.rows()-rank).setZero();
189 
190  // Apply the column permutation
191  if (m_perm_c.size()) dest = colsPermutation() * y.topRows(cols());
192  else dest = y.topRows(cols());
193 
194  m_info = Success;
195  return true;
196  }
197 
203  void setPivotThreshold(const RealScalar& threshold)
204  {
205  m_useDefaultThreshold = false;
206  m_threshold = threshold;
207  }
208 
213  template<typename Rhs>
214  inline const Solve<SparseQR, Rhs> solve(const MatrixBase<Rhs>& B) const
215  {
216  eigen_assert(m_isInitialized && "The factorization should be called first, use compute()");
217  eigen_assert(this->rows() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix");
218  return Solve<SparseQR, Rhs>(*this, B.derived());
219  }
220  template<typename Rhs>
221  inline const Solve<SparseQR, Rhs> solve(const SparseMatrixBase<Rhs>& B) const
222  {
223  eigen_assert(m_isInitialized && "The factorization should be called first, use compute()");
224  eigen_assert(this->rows() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix");
225  return Solve<SparseQR, Rhs>(*this, B.derived());
226  }
227 
237  {
238  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
239  return m_info;
240  }
241 
242 
244  inline void _sort_matrix_Q()
245  {
246  if(this->m_isQSorted) return;
247  // The matrix Q is sorted during the transposition
249  this->m_Q = mQrm;
250  this->m_isQSorted = true;
251  }
252 
253 
254  protected:
255  bool m_analysisIsok;
256  bool m_factorizationIsok;
257  mutable ComputationInfo m_info;
258  std::string m_lastError;
259  QRMatrixType m_pmat; // Temporary matrix
260  QRMatrixType m_R; // The triangular factor matrix
261  QRMatrixType m_Q; // The orthogonal reflectors
262  ScalarVector m_hcoeffs; // The Householder coefficients
263  PermutationType m_perm_c; // Fill-reducing Column permutation
264  PermutationType m_pivotperm; // The permutation for rank revealing
265  PermutationType m_outputPerm_c; // The final column permutation
266  RealScalar m_threshold; // Threshold to determine null Householder reflections
267  bool m_useDefaultThreshold; // Use default threshold
268  Index m_nonzeropivots; // Number of non zero pivots found
269  IndexVector m_etree; // Column elimination tree
270  IndexVector m_firstRowElt; // First element in each row
271  bool m_isQSorted; // whether Q is sorted or not
272  bool m_isEtreeOk; // whether the elimination tree match the initial input matrix
273 
274  template <typename, typename > friend struct SparseQR_QProduct;
275 
276 };
277 
287 template <typename MatrixType, typename OrderingType>
289 {
290  eigen_assert(mat.isCompressed() && "SparseQR requires a sparse matrix in compressed mode. Call .makeCompressed() before passing it to SparseQR");
291  // Copy to a column major matrix if the input is rowmajor
292  typename internal::conditional<MatrixType::IsRowMajor,QRMatrixType,const MatrixType&>::type matCpy(mat);
293  // Compute the column fill reducing ordering
294  OrderingType ord;
295  ord(matCpy, m_perm_c);
296  Index n = mat.cols();
297  Index m = mat.rows();
298  Index diagSize = (std::min)(m,n);
299 
300  if (!m_perm_c.size())
301  {
302  m_perm_c.resize(n);
303  m_perm_c.indices().setLinSpaced(n, 0,StorageIndex(n-1));
304  }
305 
306  // Compute the column elimination tree of the permuted matrix
307  m_outputPerm_c = m_perm_c.inverse();
308  internal::coletree(matCpy, m_etree, m_firstRowElt, m_outputPerm_c.indices().data());
309  m_isEtreeOk = true;
310 
311  m_R.resize(m, n);
312  m_Q.resize(m, diagSize);
313 
314  // Allocate space for nonzero elements : rough estimation
315  m_R.reserve(2*mat.nonZeros()); //FIXME Get a more accurate estimation through symbolic factorization with the etree
316  m_Q.reserve(2*mat.nonZeros());
317  m_hcoeffs.resize(diagSize);
318  m_analysisIsok = true;
319 }
320 
328 template <typename MatrixType, typename OrderingType>
330 {
331  using std::abs;
332 
333  eigen_assert(m_analysisIsok && "analyzePattern() should be called before this step");
334  StorageIndex m = StorageIndex(mat.rows());
335  StorageIndex n = StorageIndex(mat.cols());
336  StorageIndex diagSize = (std::min)(m,n);
337  IndexVector mark((std::max)(m,n)); mark.setConstant(-1); // Record the visited nodes
338  IndexVector Ridx(n), Qidx(m); // Store temporarily the row indexes for the current column of R and Q
339  Index nzcolR, nzcolQ; // Number of nonzero for the current column of R and Q
340  ScalarVector tval(m); // The dense vector used to compute the current column
341  RealScalar pivotThreshold = m_threshold;
342 
343  m_R.setZero();
344  m_Q.setZero();
345  m_pmat = mat;
346  if(!m_isEtreeOk)
347  {
348  m_outputPerm_c = m_perm_c.inverse();
349  internal::coletree(m_pmat, m_etree, m_firstRowElt, m_outputPerm_c.indices().data());
350  m_isEtreeOk = true;
351  }
352 
353  m_pmat.uncompress(); // To have the innerNonZeroPtr allocated
354 
355  // Apply the fill-in reducing permutation lazily:
356  {
357  // If the input is row major, copy the original column indices,
358  // otherwise directly use the input matrix
359  //
360  IndexVector originalOuterIndicesCpy;
361  const StorageIndex *originalOuterIndices = mat.outerIndexPtr();
362  if(MatrixType::IsRowMajor)
363  {
364  originalOuterIndicesCpy = IndexVector::Map(m_pmat.outerIndexPtr(),n+1);
365  originalOuterIndices = originalOuterIndicesCpy.data();
366  }
367 
368  for (int i = 0; i < n; i++)
369  {
370  Index p = m_perm_c.size() ? m_perm_c.indices()(i) : i;
371  m_pmat.outerIndexPtr()[p] = originalOuterIndices[i];
372  m_pmat.innerNonZeroPtr()[p] = originalOuterIndices[i+1] - originalOuterIndices[i];
373  }
374  }
375 
376  /* Compute the default threshold as in MatLab, see:
377  * Tim Davis, "Algorithm 915, SuiteSparseQR: Multifrontal Multithreaded Rank-Revealing
378  * Sparse QR Factorization, ACM Trans. on Math. Soft. 38(1), 2011, Page 8:3
379  */
380  if(m_useDefaultThreshold)
381  {
382  RealScalar max2Norm = 0.0;
383  for (int j = 0; j < n; j++) max2Norm = numext::maxi(max2Norm, m_pmat.col(j).norm());
384  if(max2Norm==RealScalar(0))
385  max2Norm = RealScalar(1);
386  pivotThreshold = 20 * (m + n) * max2Norm * NumTraits<RealScalar>::epsilon();
387  }
388 
389  // Initialize the numerical permutation
390  m_pivotperm.setIdentity(n);
391 
392  StorageIndex nonzeroCol = 0; // Record the number of valid pivots
393  m_Q.startVec(0);
394 
395  // Left looking rank-revealing QR factorization: compute a column of R and Q at a time
396  for (StorageIndex col = 0; col < n; ++col)
397  {
398  mark.setConstant(-1);
399  m_R.startVec(col);
400  mark(nonzeroCol) = col;
401  Qidx(0) = nonzeroCol;
402  nzcolR = 0; nzcolQ = 1;
403  bool found_diag = nonzeroCol>=m;
404  tval.setZero();
405 
406  // Symbolic factorization: find the nonzero locations of the column k of the factors R and Q, i.e.,
407  // all the nodes (with indexes lower than rank) reachable through the column elimination tree (etree) rooted at node k.
408  // Note: if the diagonal entry does not exist, then its contribution must be explicitly added,
409  // thus the trick with found_diag that permits to do one more iteration on the diagonal element if this one has not been found.
410  for (typename QRMatrixType::InnerIterator itp(m_pmat, col); itp || !found_diag; ++itp)
411  {
412  StorageIndex curIdx = nonzeroCol;
413  if(itp) curIdx = StorageIndex(itp.row());
414  if(curIdx == nonzeroCol) found_diag = true;
415 
416  // Get the nonzeros indexes of the current column of R
417  StorageIndex st = m_firstRowElt(curIdx); // The traversal of the etree starts here
418  if (st < 0 )
419  {
420  m_lastError = "Empty row found during numerical factorization";
421  m_info = InvalidInput;
422  return;
423  }
424 
425  // Traverse the etree
426  Index bi = nzcolR;
427  for (; mark(st) != col; st = m_etree(st))
428  {
429  Ridx(nzcolR) = st; // Add this row to the list,
430  mark(st) = col; // and mark this row as visited
431  nzcolR++;
432  }
433 
434  // Reverse the list to get the topological ordering
435  Index nt = nzcolR-bi;
436  for(Index i = 0; i < nt/2; i++) std::swap(Ridx(bi+i), Ridx(nzcolR-i-1));
437 
438  // Copy the current (curIdx,pcol) value of the input matrix
439  if(itp) tval(curIdx) = itp.value();
440  else tval(curIdx) = Scalar(0);
441 
442  // Compute the pattern of Q(:,k)
443  if(curIdx > nonzeroCol && mark(curIdx) != col )
444  {
445  Qidx(nzcolQ) = curIdx; // Add this row to the pattern of Q,
446  mark(curIdx) = col; // and mark it as visited
447  nzcolQ++;
448  }
449  }
450 
451  // Browse all the indexes of R(:,col) in reverse order
452  for (Index i = nzcolR-1; i >= 0; i--)
453  {
454  Index curIdx = Ridx(i);
455 
456  // Apply the curIdx-th householder vector to the current column (temporarily stored into tval)
457  Scalar tdot(0);
458 
459  // First compute q' * tval
460  tdot = m_Q.col(curIdx).dot(tval);
461 
462  tdot *= m_hcoeffs(curIdx);
463 
464  // Then update tval = tval - q * tau
465  // FIXME: tval -= tdot * m_Q.col(curIdx) should amount to the same (need to check/add support for efficient "dense ?= sparse")
466  for (typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq)
467  tval(itq.row()) -= itq.value() * tdot;
468 
469  // Detect fill-in for the current column of Q
470  if(m_etree(Ridx(i)) == nonzeroCol)
471  {
472  for (typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq)
473  {
474  StorageIndex iQ = StorageIndex(itq.row());
475  if (mark(iQ) != col)
476  {
477  Qidx(nzcolQ++) = iQ; // Add this row to the pattern of Q,
478  mark(iQ) = col; // and mark it as visited
479  }
480  }
481  }
482  } // End update current column
483 
484  Scalar tau = RealScalar(0);
485  RealScalar beta = 0;
486 
487  if(nonzeroCol < diagSize)
488  {
489  // Compute the Householder reflection that eliminate the current column
490  // FIXME this step should call the Householder module.
491  Scalar c0 = nzcolQ ? tval(Qidx(0)) : Scalar(0);
492 
493  // First, the squared norm of Q((col+1):m, col)
494  RealScalar sqrNorm = 0.;
495  for (Index itq = 1; itq < nzcolQ; ++itq) sqrNorm += numext::abs2(tval(Qidx(itq)));
496  if(sqrNorm == RealScalar(0) && numext::imag(c0) == RealScalar(0))
497  {
498  beta = numext::real(c0);
499  tval(Qidx(0)) = 1;
500  }
501  else
502  {
503  using std::sqrt;
504  beta = sqrt(numext::abs2(c0) + sqrNorm);
505  if(numext::real(c0) >= RealScalar(0))
506  beta = -beta;
507  tval(Qidx(0)) = 1;
508  for (Index itq = 1; itq < nzcolQ; ++itq)
509  tval(Qidx(itq)) /= (c0 - beta);
510  tau = numext::conj((beta-c0) / beta);
511 
512  }
513  }
514 
515  // Insert values in R
516  for (Index i = nzcolR-1; i >= 0; i--)
517  {
518  Index curIdx = Ridx(i);
519  if(curIdx < nonzeroCol)
520  {
521  m_R.insertBackByOuterInnerUnordered(col, curIdx) = tval(curIdx);
522  tval(curIdx) = Scalar(0.);
523  }
524  }
525 
526  if(nonzeroCol < diagSize && abs(beta) >= pivotThreshold)
527  {
528  m_R.insertBackByOuterInner(col, nonzeroCol) = beta;
529  // The householder coefficient
530  m_hcoeffs(nonzeroCol) = tau;
531  // Record the householder reflections
532  for (Index itq = 0; itq < nzcolQ; ++itq)
533  {
534  Index iQ = Qidx(itq);
535  m_Q.insertBackByOuterInnerUnordered(nonzeroCol,iQ) = tval(iQ);
536  tval(iQ) = Scalar(0.);
537  }
538  nonzeroCol++;
539  if(nonzeroCol<diagSize)
540  m_Q.startVec(nonzeroCol);
541  }
542  else
543  {
544  // Zero pivot found: move implicitly this column to the end
545  for (Index j = nonzeroCol; j < n-1; j++)
546  std::swap(m_pivotperm.indices()(j), m_pivotperm.indices()[j+1]);
547 
548  // Recompute the column elimination tree
549  internal::coletree(m_pmat, m_etree, m_firstRowElt, m_pivotperm.indices().data());
550  m_isEtreeOk = false;
551  }
552  }
553 
554  m_hcoeffs.tail(diagSize-nonzeroCol).setZero();
555 
556  // Finalize the column pointers of the sparse matrices R and Q
557  m_Q.finalize();
558  m_Q.makeCompressed();
559  m_R.finalize();
560  m_R.makeCompressed();
561  m_isQSorted = false;
562 
563  m_nonzeropivots = nonzeroCol;
564 
565  if(nonzeroCol<n)
566  {
567  // Permute the triangular factor to put the 'dead' columns to the end
568  QRMatrixType tempR(m_R);
569  m_R = tempR * m_pivotperm;
570 
571  // Update the column permutation
572  m_outputPerm_c = m_outputPerm_c * m_pivotperm;
573  }
574 
575  m_isInitialized = true;
576  m_factorizationIsok = true;
577  m_info = Success;
578 }
579 
580 template <typename SparseQRType, typename Derived>
581 struct SparseQR_QProduct : ReturnByValue<SparseQR_QProduct<SparseQRType, Derived> >
582 {
583  typedef typename SparseQRType::QRMatrixType MatrixType;
584  typedef typename SparseQRType::Scalar Scalar;
585  // Get the references
586  SparseQR_QProduct(const SparseQRType& qr, const Derived& other, bool transpose) :
587  m_qr(qr),m_other(other),m_transpose(transpose) {}
588  inline Index rows() const { return m_transpose ? m_qr.rows() : m_qr.cols(); }
589  inline Index cols() const { return m_other.cols(); }
590 
591  // Assign to a vector
592  template<typename DesType>
593  void evalTo(DesType& res) const
594  {
595  Index m = m_qr.rows();
596  Index n = m_qr.cols();
597  Index diagSize = (std::min)(m,n);
598  res = m_other;
599  if (m_transpose)
600  {
601  eigen_assert(m_qr.m_Q.rows() == m_other.rows() && "Non conforming object sizes");
602  //Compute res = Q' * other column by column
603  for(Index j = 0; j < res.cols(); j++){
604  for (Index k = 0; k < diagSize; k++)
605  {
606  Scalar tau = Scalar(0);
607  tau = m_qr.m_Q.col(k).dot(res.col(j));
608  if(tau==Scalar(0)) continue;
609  tau = tau * m_qr.m_hcoeffs(k);
610  res.col(j) -= tau * m_qr.m_Q.col(k);
611  }
612  }
613  }
614  else
615  {
616  eigen_assert(m_qr.m_Q.rows() == m_other.rows() && "Non conforming object sizes");
617  // Compute res = Q * other column by column
618  for(Index j = 0; j < res.cols(); j++)
619  {
620  for (Index k = diagSize-1; k >=0; k--)
621  {
622  Scalar tau = Scalar(0);
623  tau = m_qr.m_Q.col(k).dot(res.col(j));
624  if(tau==Scalar(0)) continue;
625  tau = tau * m_qr.m_hcoeffs(k);
626  res.col(j) -= tau * m_qr.m_Q.col(k);
627  }
628  }
629  }
630  }
631 
632  const SparseQRType& m_qr;
633  const Derived& m_other;
634  bool m_transpose;
635 };
636 
637 template<typename SparseQRType>
638 struct SparseQRMatrixQReturnType : public EigenBase<SparseQRMatrixQReturnType<SparseQRType> >
639 {
640  typedef typename SparseQRType::Scalar Scalar;
641  typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
642  enum {
643  RowsAtCompileTime = Dynamic,
644  ColsAtCompileTime = Dynamic
645  };
646  explicit SparseQRMatrixQReturnType(const SparseQRType& qr) : m_qr(qr) {}
647  template<typename Derived>
648  SparseQR_QProduct<SparseQRType, Derived> operator*(const MatrixBase<Derived>& other)
649  {
650  return SparseQR_QProduct<SparseQRType,Derived>(m_qr,other.derived(),false);
651  }
652  SparseQRMatrixQTransposeReturnType<SparseQRType> adjoint() const
653  {
654  return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr);
655  }
656  inline Index rows() const { return m_qr.rows(); }
657  inline Index cols() const { return (std::min)(m_qr.rows(),m_qr.cols()); }
658  // To use for operations with the transpose of Q
659  SparseQRMatrixQTransposeReturnType<SparseQRType> transpose() const
660  {
661  return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr);
662  }
663  const SparseQRType& m_qr;
664 };
665 
666 template<typename SparseQRType>
667 struct SparseQRMatrixQTransposeReturnType
668 {
669  explicit SparseQRMatrixQTransposeReturnType(const SparseQRType& qr) : m_qr(qr) {}
670  template<typename Derived>
671  SparseQR_QProduct<SparseQRType,Derived> operator*(const MatrixBase<Derived>& other)
672  {
673  return SparseQR_QProduct<SparseQRType,Derived>(m_qr,other.derived(), true);
674  }
675  const SparseQRType& m_qr;
676 };
677 
678 namespace internal {
679 
680 template<typename SparseQRType>
681 struct evaluator_traits<SparseQRMatrixQReturnType<SparseQRType> >
682 {
683  typedef typename SparseQRType::MatrixType MatrixType;
684  typedef typename storage_kind_to_evaluator_kind<typename MatrixType::StorageKind>::Kind Kind;
685  typedef SparseShape Shape;
686  static const int AssumeAliasing = 0;
687 };
688 
689 template< typename DstXprType, typename SparseQRType>
690 struct Assignment<DstXprType, SparseQRMatrixQReturnType<SparseQRType>, internal::assign_op<typename DstXprType::Scalar>, Sparse2Sparse>
691 {
692  typedef SparseQRMatrixQReturnType<SparseQRType> SrcXprType;
693  typedef typename DstXprType::Scalar Scalar;
694  typedef typename DstXprType::StorageIndex StorageIndex;
695  static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar> &/*func*/)
696  {
697  typename DstXprType::PlainObject idMat(src.m_qr.rows(), src.m_qr.rows());
698  idMat.setIdentity();
699  // Sort the sparse householder reflectors if needed
700  const_cast<SparseQRType *>(&src.m_qr)->_sort_matrix_Q();
701  dst = SparseQR_QProduct<SparseQRType, DstXprType>(src.m_qr, idMat, false);
702  }
703 };
704 
705 template< typename DstXprType, typename SparseQRType>
706 struct Assignment<DstXprType, SparseQRMatrixQReturnType<SparseQRType>, internal::assign_op<typename DstXprType::Scalar>, Sparse2Dense>
707 {
708  typedef SparseQRMatrixQReturnType<SparseQRType> SrcXprType;
709  typedef typename DstXprType::Scalar Scalar;
710  typedef typename DstXprType::StorageIndex StorageIndex;
711  static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar> &/*func*/)
712  {
713  dst = src.m_qr.matrixQ() * DstXprType::Identity(src.m_qr.rows(), src.m_qr.rows());
714  }
715 };
716 
717 } // end namespace internal
718 
719 } // end namespace Eigen
720 
721 #endif
Index cols() const
Definition: SparseMatrix.h:132
A versatible sparse matrix representation.
Definition: SparseMatrix.h:92
RowsBlockXpr topRows(Index n)
Definition: DenseBase.h:400
Index rank() const
Definition: SparseQR.h:130
Index size() const
Definition: PermutationMatrix.h:110
Definition: LDLT.h:16
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:107
Derived & derived()
Definition: EigenBase.h:44
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:37
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SparseQR.h:236
const QRMatrixType & matrixR() const
Definition: SparseQR.h:124
void resize(Index newSize)
Definition: DenseBase.h:252
void factorize(const MatrixType &mat)
Performs the numerical QR factorization of the input matrix.
Definition: SparseQR.h:329
Derived & setZero(Index size)
Definition: CwiseNullaryOp.h:520
Sparse left-looking rank-revealing QR factorization.
Definition: SparseQR.h:16
Base class of any sparse matrices or sparse expressions.
Definition: SparseMatrixBase.h:26
Index cols() const
Definition: SparseQR.h:120
Definition: Constants.h:431
Definition: Constants.h:424
void analyzePattern(const MatrixType &mat)
Preprocessing step of a QR factorization.
Definition: SparseQR.h:288
Block< Derived > topLeftCorner(Index cRows, Index cCols)
Definition: SparseMatrixBase.h:164
Definition: Eigen_Colamd.h:54
Index rows() const
Definition: SparseQR.h:116
const PermutationType & colsPermutation() const
Definition: SparseQR.h:160
const Solve< SparseQR, Rhs > solve(const MatrixBase< Rhs > &B) const
Definition: SparseQR.h:214
SparseQRMatrixQReturnType< SparseQR > matrixQ() const
Definition: SparseQR.h:154
Index rows() const
Definition: SparseMatrixBase.h:152
Pseudo expression representing a solving operation.
Definition: Solve.h:63
const Scalar * data() const
Definition: PlainObjectBase.h:228
ComputationInfo
Definition: Constants.h:422
SparseQR(const MatrixType &mat)
Definition: SparseQR.h:95
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
Index rows() const
Definition: SparseMatrix.h:130
void setPivotThreshold(const RealScalar &threshold)
Definition: SparseQR.h:203
std::string lastErrorMessage() const
Definition: SparseQR.h:169
Derived & setConstant(Index size, const Scalar &value)
Definition: CwiseNullaryOp.h:352
void compute(const MatrixType &mat)
Definition: SparseQR.h:106