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