Eigen  3.4.90 (git rev 67eeba6e720c5745abc77ae6c92ce0a44aa7b7ae)
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 #include "./InternalHeaderCheck.h"
15 
16 namespace Eigen {
17 
18 template<typename MatrixType, typename OrderingType> class SparseQR;
19 template<typename SparseQRType> struct SparseQRMatrixQReturnType;
20 template<typename SparseQRType> struct SparseQRMatrixQTransposeReturnType;
21 template<typename SparseQRType, typename Derived> struct SparseQR_QProduct;
22 namespace internal {
23  template <typename SparseQRType> struct traits<SparseQRMatrixQReturnType<SparseQRType> >
24  {
25  typedef typename SparseQRType::MatrixType ReturnType;
26  typedef typename ReturnType::StorageIndex StorageIndex;
27  typedef typename ReturnType::StorageKind StorageKind;
28  enum {
29  RowsAtCompileTime = Dynamic,
30  ColsAtCompileTime = Dynamic
31  };
32  };
33  template <typename SparseQRType> struct traits<SparseQRMatrixQTransposeReturnType<SparseQRType> >
34  {
35  typedef typename SparseQRType::MatrixType ReturnType;
36  };
37  template <typename SparseQRType, typename Derived> struct traits<SparseQR_QProduct<SparseQRType, Derived> >
38  {
39  typedef typename Derived::PlainObject ReturnType;
40  };
41 } // End namespace internal
42 
85 template<typename MatrixType_, typename OrderingType_>
86 class SparseQR : public SparseSolverBase<SparseQR<MatrixType_,OrderingType_> >
87 {
88  protected:
90  using Base::m_isInitialized;
91  public:
92  using Base::_solve_impl;
93  typedef MatrixType_ MatrixType;
94  typedef OrderingType_ OrderingType;
95  typedef typename MatrixType::Scalar Scalar;
96  typedef typename MatrixType::RealScalar RealScalar;
97  typedef typename MatrixType::StorageIndex StorageIndex;
102 
103  enum {
104  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
105  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
106  };
107 
108  public:
109  SparseQR () : m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false),m_isEtreeOk(false)
110  { }
111 
118  explicit SparseQR(const MatrixType& mat) : m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false),m_isEtreeOk(false)
119  {
120  compute(mat);
121  }
122 
129  void compute(const MatrixType& mat)
130  {
131  analyzePattern(mat);
132  factorize(mat);
133  }
134  void analyzePattern(const MatrixType& mat);
135  void factorize(const MatrixType& mat);
136 
139  inline Index rows() const { return m_pmat.rows(); }
140 
143  inline Index cols() const { return m_pmat.cols();}
144 
158  const QRMatrixType& matrixR() const { return m_R; }
159 
164  Index rank() const
165  {
166  eigen_assert(m_isInitialized && "The factorization should be called first, use compute()");
167  return m_nonzeropivots;
168  }
169 
188  SparseQRMatrixQReturnType<SparseQR> matrixQ() const
189  { return SparseQRMatrixQReturnType<SparseQR>(*this); }
190 
195  {
196  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
197  return m_outputPerm_c;
198  }
199 
203  std::string lastErrorMessage() const { return m_lastError; }
204 
206  template<typename Rhs, typename Dest>
207  bool _solve_impl(const MatrixBase<Rhs> &B, MatrixBase<Dest> &dest) const
208  {
209  eigen_assert(m_isInitialized && "The factorization should be called first, use compute()");
210  eigen_assert(this->rows() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix");
211 
212  Index rank = this->rank();
213 
214  // Compute Q^* * b;
215  typename Dest::PlainObject y, b;
216  y = this->matrixQ().adjoint() * B;
217  b = y;
218 
219  // Solve with the triangular matrix R
220  y.resize((std::max<Index>)(cols(),y.rows()),y.cols());
221  y.topRows(rank) = this->matrixR().topLeftCorner(rank, rank).template triangularView<Upper>().solve(b.topRows(rank));
222  y.bottomRows(y.rows()-rank).setZero();
223 
224  // Apply the column permutation
225  if (m_perm_c.size()) dest = colsPermutation() * y.topRows(cols());
226  else dest = y.topRows(cols());
227 
228  m_info = Success;
229  return true;
230  }
231 
237  void setPivotThreshold(const RealScalar& threshold)
238  {
239  m_useDefaultThreshold = false;
240  m_threshold = threshold;
241  }
242 
247  template<typename Rhs>
248  inline const Solve<SparseQR, Rhs> solve(const MatrixBase<Rhs>& B) const
249  {
250  eigen_assert(m_isInitialized && "The factorization should be called first, use compute()");
251  eigen_assert(this->rows() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix");
252  return Solve<SparseQR, Rhs>(*this, B.derived());
253  }
254  template<typename Rhs>
255  inline const Solve<SparseQR, Rhs> solve(const SparseMatrixBase<Rhs>& B) const
256  {
257  eigen_assert(m_isInitialized && "The factorization should be called first, use compute()");
258  eigen_assert(this->rows() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix");
259  return Solve<SparseQR, Rhs>(*this, B.derived());
260  }
261 
271  {
272  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
273  return m_info;
274  }
275 
276 
278  inline void _sort_matrix_Q()
279  {
280  if(this->m_isQSorted) return;
281  // The matrix Q is sorted during the transposition
283  this->m_Q = mQrm;
284  this->m_isQSorted = true;
285  }
286 
287 
288  protected:
289  bool m_analysisIsok;
290  bool m_factorizationIsok;
291  mutable ComputationInfo m_info;
292  std::string m_lastError;
293  QRMatrixType m_pmat; // Temporary matrix
294  QRMatrixType m_R; // The triangular factor matrix
295  QRMatrixType m_Q; // The orthogonal reflectors
296  ScalarVector m_hcoeffs; // The Householder coefficients
297  PermutationType m_perm_c; // Fill-reducing Column permutation
298  PermutationType m_pivotperm; // The permutation for rank revealing
299  PermutationType m_outputPerm_c; // The final column permutation
300  RealScalar m_threshold; // Threshold to determine null Householder reflections
301  bool m_useDefaultThreshold; // Use default threshold
302  Index m_nonzeropivots; // Number of non zero pivots found
303  IndexVector m_etree; // Column elimination tree
304  IndexVector m_firstRowElt; // First element in each row
305  bool m_isQSorted; // whether Q is sorted or not
306  bool m_isEtreeOk; // whether the elimination tree match the initial input matrix
307 
308  template <typename, typename > friend struct SparseQR_QProduct;
309 
310 };
311 
321 template <typename MatrixType, typename OrderingType>
323 {
324  eigen_assert(mat.isCompressed() && "SparseQR requires a sparse matrix in compressed mode. Call .makeCompressed() before passing it to SparseQR");
325  // Copy to a column major matrix if the input is rowmajor
326  std::conditional_t<MatrixType::IsRowMajor,QRMatrixType,const MatrixType&> matCpy(mat);
327  // Compute the column fill reducing ordering
328  OrderingType ord;
329  ord(matCpy, m_perm_c);
330  Index n = mat.cols();
331  Index m = mat.rows();
332  Index diagSize = (std::min)(m,n);
333 
334  if (!m_perm_c.size())
335  {
336  m_perm_c.resize(n);
337  m_perm_c.indices().setLinSpaced(n, 0,StorageIndex(n-1));
338  }
339 
340  // Compute the column elimination tree of the permuted matrix
341  m_outputPerm_c = m_perm_c.inverse();
342  internal::coletree(matCpy, m_etree, m_firstRowElt, m_outputPerm_c.indices().data());
343  m_isEtreeOk = true;
344 
345  m_R.resize(m, n);
346  m_Q.resize(m, diagSize);
347 
348  // Allocate space for nonzero elements: rough estimation
349  m_R.reserve(2*mat.nonZeros()); //FIXME Get a more accurate estimation through symbolic factorization with the etree
350  m_Q.reserve(2*mat.nonZeros());
351  m_hcoeffs.resize(diagSize);
352  m_analysisIsok = true;
353 }
354 
362 template <typename MatrixType, typename OrderingType>
364 {
365  using std::abs;
366 
367  eigen_assert(m_analysisIsok && "analyzePattern() should be called before this step");
368  StorageIndex m = StorageIndex(mat.rows());
369  StorageIndex n = StorageIndex(mat.cols());
370  StorageIndex diagSize = (std::min)(m,n);
371  IndexVector mark((std::max)(m,n)); mark.setConstant(-1); // Record the visited nodes
372  IndexVector Ridx(n), Qidx(m); // Store temporarily the row indexes for the current column of R and Q
373  Index nzcolR, nzcolQ; // Number of nonzero for the current column of R and Q
374  ScalarVector tval(m); // The dense vector used to compute the current column
375  RealScalar pivotThreshold = m_threshold;
376 
377  m_R.setZero();
378  m_Q.setZero();
379  m_pmat = mat;
380  if(!m_isEtreeOk)
381  {
382  m_outputPerm_c = m_perm_c.inverse();
383  internal::coletree(m_pmat, m_etree, m_firstRowElt, m_outputPerm_c.indices().data());
384  m_isEtreeOk = true;
385  }
386 
387  m_pmat.uncompress(); // To have the innerNonZeroPtr allocated
388 
389  // Apply the fill-in reducing permutation lazily:
390  {
391  // If the input is row major, copy the original column indices,
392  // otherwise directly use the input matrix
393  //
394  IndexVector originalOuterIndicesCpy;
395  const StorageIndex *originalOuterIndices = mat.outerIndexPtr();
396  if(MatrixType::IsRowMajor)
397  {
398  originalOuterIndicesCpy = IndexVector::Map(m_pmat.outerIndexPtr(),n+1);
399  originalOuterIndices = originalOuterIndicesCpy.data();
400  }
401 
402  for (int i = 0; i < n; i++)
403  {
404  Index p = m_perm_c.size() ? m_perm_c.indices()(i) : i;
405  m_pmat.outerIndexPtr()[p] = originalOuterIndices[i];
406  m_pmat.innerNonZeroPtr()[p] = originalOuterIndices[i+1] - originalOuterIndices[i];
407  }
408  }
409 
410  /* Compute the default threshold as in MatLab, see:
411  * Tim Davis, "Algorithm 915, SuiteSparseQR: Multifrontal Multithreaded Rank-Revealing
412  * Sparse QR Factorization, ACM Trans. on Math. Soft. 38(1), 2011, Page 8:3
413  */
414  if(m_useDefaultThreshold)
415  {
416  RealScalar max2Norm = 0.0;
417  for (int j = 0; j < n; j++) max2Norm = numext::maxi(max2Norm, m_pmat.col(j).norm());
418  if(max2Norm==RealScalar(0))
419  max2Norm = RealScalar(1);
420  pivotThreshold = 20 * (m + n) * max2Norm * NumTraits<RealScalar>::epsilon();
421  }
422 
423  // Initialize the numerical permutation
424  m_pivotperm.setIdentity(n);
425 
426  StorageIndex nonzeroCol = 0; // Record the number of valid pivots
427  m_Q.startVec(0);
428 
429  // Left looking rank-revealing QR factorization: compute a column of R and Q at a time
430  for (StorageIndex col = 0; col < n; ++col)
431  {
432  mark.setConstant(-1);
433  m_R.startVec(col);
434  mark(nonzeroCol) = col;
435  Qidx(0) = nonzeroCol;
436  nzcolR = 0; nzcolQ = 1;
437  bool found_diag = nonzeroCol>=m;
438  tval.setZero();
439 
440  // Symbolic factorization: find the nonzero locations of the column k of the factors R and Q, i.e.,
441  // all the nodes (with indexes lower than rank) reachable through the column elimination tree (etree) rooted at node k.
442  // Note: if the diagonal entry does not exist, then its contribution must be explicitly added,
443  // thus the trick with found_diag that permits to do one more iteration on the diagonal element if this one has not been found.
444  for (typename QRMatrixType::InnerIterator itp(m_pmat, col); itp || !found_diag; ++itp)
445  {
446  StorageIndex curIdx = nonzeroCol;
447  if(itp) curIdx = StorageIndex(itp.row());
448  if(curIdx == nonzeroCol) found_diag = true;
449 
450  // Get the nonzeros indexes of the current column of R
451  StorageIndex st = m_firstRowElt(curIdx); // The traversal of the etree starts here
452  if (st < 0 )
453  {
454  m_lastError = "Empty row found during numerical factorization";
455  m_info = InvalidInput;
456  return;
457  }
458 
459  // Traverse the etree
460  Index bi = nzcolR;
461  for (; mark(st) != col; st = m_etree(st))
462  {
463  Ridx(nzcolR) = st; // Add this row to the list,
464  mark(st) = col; // and mark this row as visited
465  nzcolR++;
466  }
467 
468  // Reverse the list to get the topological ordering
469  Index nt = nzcolR-bi;
470  for(Index i = 0; i < nt/2; i++) std::swap(Ridx(bi+i), Ridx(nzcolR-i-1));
471 
472  // Copy the current (curIdx,pcol) value of the input matrix
473  if(itp) tval(curIdx) = itp.value();
474  else tval(curIdx) = Scalar(0);
475 
476  // Compute the pattern of Q(:,k)
477  if(curIdx > nonzeroCol && mark(curIdx) != col )
478  {
479  Qidx(nzcolQ) = curIdx; // Add this row to the pattern of Q,
480  mark(curIdx) = col; // and mark it as visited
481  nzcolQ++;
482  }
483  }
484 
485  // Browse all the indexes of R(:,col) in reverse order
486  for (Index i = nzcolR-1; i >= 0; i--)
487  {
488  Index curIdx = Ridx(i);
489 
490  // Apply the curIdx-th householder vector to the current column (temporarily stored into tval)
491  Scalar tdot(0);
492 
493  // First compute q' * tval
494  tdot = m_Q.col(curIdx).dot(tval);
495 
496  tdot *= m_hcoeffs(curIdx);
497 
498  // Then update tval = tval - q * tau
499  // FIXME: tval -= tdot * m_Q.col(curIdx) should amount to the same (need to check/add support for efficient "dense ?= sparse")
500  for (typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq)
501  tval(itq.row()) -= itq.value() * tdot;
502 
503  // Detect fill-in for the current column of Q
504  if(m_etree(Ridx(i)) == nonzeroCol)
505  {
506  for (typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq)
507  {
508  StorageIndex iQ = StorageIndex(itq.row());
509  if (mark(iQ) != col)
510  {
511  Qidx(nzcolQ++) = iQ; // Add this row to the pattern of Q,
512  mark(iQ) = col; // and mark it as visited
513  }
514  }
515  }
516  } // End update current column
517 
518  Scalar tau = RealScalar(0);
519  RealScalar beta = 0;
520 
521  if(nonzeroCol < diagSize)
522  {
523  // Compute the Householder reflection that eliminate the current column
524  // FIXME this step should call the Householder module.
525  Scalar c0 = nzcolQ ? tval(Qidx(0)) : Scalar(0);
526 
527  // First, the squared norm of Q((col+1):m, col)
528  RealScalar sqrNorm = 0.;
529  for (Index itq = 1; itq < nzcolQ; ++itq) sqrNorm += numext::abs2(tval(Qidx(itq)));
530  if(sqrNorm == RealScalar(0) && numext::imag(c0) == RealScalar(0))
531  {
532  beta = numext::real(c0);
533  tval(Qidx(0)) = 1;
534  }
535  else
536  {
537  using std::sqrt;
538  beta = sqrt(numext::abs2(c0) + sqrNorm);
539  if(numext::real(c0) >= RealScalar(0))
540  beta = -beta;
541  tval(Qidx(0)) = 1;
542  for (Index itq = 1; itq < nzcolQ; ++itq)
543  tval(Qidx(itq)) /= (c0 - beta);
544  tau = numext::conj((beta-c0) / beta);
545 
546  }
547  }
548 
549  // Insert values in R
550  for (Index i = nzcolR-1; i >= 0; i--)
551  {
552  Index curIdx = Ridx(i);
553  if(curIdx < nonzeroCol)
554  {
555  m_R.insertBackByOuterInnerUnordered(col, curIdx) = tval(curIdx);
556  tval(curIdx) = Scalar(0.);
557  }
558  }
559 
560  if(nonzeroCol < diagSize && abs(beta) >= pivotThreshold)
561  {
562  m_R.insertBackByOuterInner(col, nonzeroCol) = beta;
563  // The householder coefficient
564  m_hcoeffs(nonzeroCol) = tau;
565  // Record the householder reflections
566  for (Index itq = 0; itq < nzcolQ; ++itq)
567  {
568  Index iQ = Qidx(itq);
569  m_Q.insertBackByOuterInnerUnordered(nonzeroCol,iQ) = tval(iQ);
570  tval(iQ) = Scalar(0.);
571  }
572  nonzeroCol++;
573  if(nonzeroCol<diagSize)
574  m_Q.startVec(nonzeroCol);
575  }
576  else
577  {
578  // Zero pivot found: move implicitly this column to the end
579  for (Index j = nonzeroCol; j < n-1; j++)
580  std::swap(m_pivotperm.indices()(j), m_pivotperm.indices()[j+1]);
581 
582  // Recompute the column elimination tree
583  internal::coletree(m_pmat, m_etree, m_firstRowElt, m_pivotperm.indices().data());
584  m_isEtreeOk = false;
585  }
586  }
587 
588  m_hcoeffs.tail(diagSize-nonzeroCol).setZero();
589 
590  // Finalize the column pointers of the sparse matrices R and Q
591  m_Q.finalize();
592  m_Q.makeCompressed();
593  m_R.finalize();
594  m_R.makeCompressed();
595  m_isQSorted = false;
596 
597  m_nonzeropivots = nonzeroCol;
598 
599  if(nonzeroCol<n)
600  {
601  // Permute the triangular factor to put the 'dead' columns to the end
602  QRMatrixType tempR(m_R);
603  m_R = tempR * m_pivotperm;
604 
605  // Update the column permutation
606  m_outputPerm_c = m_outputPerm_c * m_pivotperm;
607  }
608 
609  m_isInitialized = true;
610  m_factorizationIsok = true;
611  m_info = Success;
612 }
613 
614 template <typename SparseQRType, typename Derived>
615 struct SparseQR_QProduct : ReturnByValue<SparseQR_QProduct<SparseQRType, Derived> >
616 {
617  typedef typename SparseQRType::QRMatrixType MatrixType;
618  typedef typename SparseQRType::Scalar Scalar;
619  // Get the references
620  SparseQR_QProduct(const SparseQRType& qr, const Derived& other, bool transpose) :
621  m_qr(qr),m_other(other),m_transpose(transpose) {}
622  inline Index rows() const { return m_qr.matrixQ().rows(); }
623  inline Index cols() const { return m_other.cols(); }
624 
625  // Assign to a vector
626  template<typename DesType>
627  void evalTo(DesType& res) const
628  {
629  Index m = m_qr.rows();
630  Index n = m_qr.cols();
631  Index diagSize = (std::min)(m,n);
632  res = m_other;
633  if (m_transpose)
634  {
635  eigen_assert(m_qr.m_Q.rows() == m_other.rows() && "Non conforming object sizes");
636  //Compute res = Q' * other column by column
637  for(Index j = 0; j < res.cols(); j++){
638  for (Index k = 0; k < diagSize; k++)
639  {
640  Scalar tau = Scalar(0);
641  tau = m_qr.m_Q.col(k).dot(res.col(j));
642  if(tau==Scalar(0)) continue;
643  tau = tau * m_qr.m_hcoeffs(k);
644  res.col(j) -= tau * m_qr.m_Q.col(k);
645  }
646  }
647  }
648  else
649  {
650  eigen_assert(m_qr.matrixQ().cols() == m_other.rows() && "Non conforming object sizes");
651 
652  res.conservativeResize(rows(), cols());
653 
654  // Compute res = Q * other column by column
655  for(Index j = 0; j < res.cols(); j++)
656  {
657  Index start_k = internal::is_identity<Derived>::value ? numext::mini(j,diagSize-1) : diagSize-1;
658  for (Index k = start_k; k >=0; k--)
659  {
660  Scalar tau = Scalar(0);
661  tau = m_qr.m_Q.col(k).dot(res.col(j));
662  if(tau==Scalar(0)) continue;
663  tau = tau * numext::conj(m_qr.m_hcoeffs(k));
664  res.col(j) -= tau * m_qr.m_Q.col(k);
665  }
666  }
667  }
668  }
669 
670  const SparseQRType& m_qr;
671  const Derived& m_other;
672  bool m_transpose; // TODO this actually means adjoint
673 };
674 
675 template<typename SparseQRType>
676 struct SparseQRMatrixQReturnType : public EigenBase<SparseQRMatrixQReturnType<SparseQRType> >
677 {
678  typedef typename SparseQRType::Scalar Scalar;
679  typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
680  enum {
681  RowsAtCompileTime = Dynamic,
682  ColsAtCompileTime = Dynamic
683  };
684  explicit SparseQRMatrixQReturnType(const SparseQRType& qr) : m_qr(qr) {}
685  template<typename Derived>
686  SparseQR_QProduct<SparseQRType, Derived> operator*(const MatrixBase<Derived>& other)
687  {
688  return SparseQR_QProduct<SparseQRType,Derived>(m_qr,other.derived(),false);
689  }
690  // To use for operations with the adjoint of Q
691  SparseQRMatrixQTransposeReturnType<SparseQRType> adjoint() const
692  {
693  return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr);
694  }
695  inline Index rows() const { return m_qr.rows(); }
696  inline Index cols() const { return m_qr.rows(); }
697  // To use for operations with the transpose of Q FIXME this is the same as adjoint at the moment
698  SparseQRMatrixQTransposeReturnType<SparseQRType> transpose() const
699  {
700  return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr);
701  }
702  const SparseQRType& m_qr;
703 };
704 
705 // TODO this actually represents the adjoint of Q
706 template<typename SparseQRType>
707 struct SparseQRMatrixQTransposeReturnType
708 {
709  explicit SparseQRMatrixQTransposeReturnType(const SparseQRType& qr) : m_qr(qr) {}
710  template<typename Derived>
711  SparseQR_QProduct<SparseQRType,Derived> operator*(const MatrixBase<Derived>& other)
712  {
713  return SparseQR_QProduct<SparseQRType,Derived>(m_qr,other.derived(), true);
714  }
715  const SparseQRType& m_qr;
716 };
717 
718 namespace internal {
719 
720 template<typename SparseQRType>
721 struct evaluator_traits<SparseQRMatrixQReturnType<SparseQRType> >
722 {
723  typedef typename SparseQRType::MatrixType MatrixType;
724  typedef typename storage_kind_to_evaluator_kind<typename MatrixType::StorageKind>::Kind Kind;
725  typedef SparseShape Shape;
726 };
727 
728 template< typename DstXprType, typename SparseQRType>
729 struct Assignment<DstXprType, SparseQRMatrixQReturnType<SparseQRType>, internal::assign_op<typename DstXprType::Scalar,typename DstXprType::Scalar>, Sparse2Sparse>
730 {
731  typedef SparseQRMatrixQReturnType<SparseQRType> SrcXprType;
732  typedef typename DstXprType::Scalar Scalar;
733  typedef typename DstXprType::StorageIndex StorageIndex;
734  static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar,Scalar> &/*func*/)
735  {
736  typename DstXprType::PlainObject idMat(src.rows(), src.cols());
737  idMat.setIdentity();
738  // Sort the sparse householder reflectors if needed
739  const_cast<SparseQRType *>(&src.m_qr)->_sort_matrix_Q();
740  dst = SparseQR_QProduct<SparseQRType, DstXprType>(src.m_qr, idMat, false);
741  }
742 };
743 
744 template< typename DstXprType, typename SparseQRType>
745 struct Assignment<DstXprType, SparseQRMatrixQReturnType<SparseQRType>, internal::assign_op<typename DstXprType::Scalar,typename DstXprType::Scalar>, Sparse2Dense>
746 {
747  typedef SparseQRMatrixQReturnType<SparseQRType> SrcXprType;
748  typedef typename DstXprType::Scalar Scalar;
749  typedef typename DstXprType::StorageIndex StorageIndex;
750  static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar,Scalar> &/*func*/)
751  {
752  dst = src.m_qr.matrixQ() * DstXprType::Identity(src.m_qr.rows(), src.m_qr.rows());
753  }
754 };
755 
756 } // end namespace internal
757 
758 } // end namespace Eigen
759 
760 #endif
Derived & derived()
Definition: EigenBase.h:48
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: EigenBase.h:62
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:52
Index size() const
Definition: PermutationMatrix.h:99
Derived & setConstant(Index size, const Scalar &val)
Definition: CwiseNullaryOp.h:363
const Scalar * data() const
Definition: PlainObjectBase.h:259
Derived & setZero(Index size)
Definition: CwiseNullaryOp.h:564
Pseudo expression representing a solving operation.
Definition: Solve.h:65
Base class of any sparse matrices or sparse expressions.
Definition: SparseMatrixBase.h:30
Index rows() const
Definition: SparseMatrixBase.h:176
Index cols() const
Definition: SparseMatrix.h:142
Index rows() const
Definition: SparseMatrix.h:140
Sparse left-looking QR factorization with numerical column pivoting.
Definition: SparseQR.h:87
void setPivotThreshold(const RealScalar &threshold)
Definition: SparseQR.h:237
void analyzePattern(const MatrixType &mat)
Preprocessing step of a QR factorization.
Definition: SparseQR.h:322
void factorize(const MatrixType &mat)
Performs the numerical QR factorization of the input matrix.
Definition: SparseQR.h:363
SparseQR(const MatrixType &mat)
Definition: SparseQR.h:118
Index cols() const
Definition: SparseQR.h:143
const QRMatrixType & matrixR() const
Definition: SparseQR.h:158
Index rows() const
Definition: SparseQR.h:139
const PermutationType & colsPermutation() const
Definition: SparseQR.h:194
SparseQRMatrixQReturnType< SparseQR > matrixQ() const
Definition: SparseQR.h:188
const Solve< SparseQR, Rhs > solve(const MatrixBase< Rhs > &B) const
Definition: SparseQR.h:248
std::string lastErrorMessage() const
Definition: SparseQR.h:203
void compute(const MatrixType &mat)
Definition: SparseQR.h:129
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SparseQR.h:270
Index rank() const
Definition: SparseQR.h:164
A base class for sparse solvers.
Definition: SparseSolverBase.h:70
ComputationInfo
Definition: Constants.h:442
@ InvalidInput
Definition: Constants.h:451
@ Success
Definition: Constants.h:444
Namespace containing all symbols from the Eigen library.
Definition: Core:139
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
const Product< MatrixDerived, PermutationDerived, AliasFreeProduct > operator*(const MatrixBase< MatrixDerived > &matrix, const PermutationBase< PermutationDerived > &permutation)
Definition: PermutationMatrix.h:517
const int Dynamic
Definition: Constants.h:24
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_sqrt_op< typename Derived::Scalar >, const Derived > sqrt(const Eigen::ArrayBase< Derived > &x)
Derived & derived()
Definition: EigenBase.h:48
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:231