Eigen::SparseQR< MatrixType_, OrderingType_ > Class Template Reference

Sparse left-looking QR factorization with numerical column pivoting. More...

+ Inheritance diagram for Eigen::SparseQR< MatrixType_, OrderingType_ >:

Public Types

enum  {
  ColsAtCompileTime ,
  MaxColsAtCompileTime
}
 
typedef Matrix< StorageIndex, Dynamic, 1 > IndexVector
 
typedef MatrixType_ MatrixType
 
typedef OrderingType_ OrderingType
 
typedef PermutationMatrix< Dynamic, Dynamic, StorageIndexPermutationType
 
typedef SparseMatrix< Scalar, ColMajor, StorageIndexQRMatrixType
 
typedef MatrixType::RealScalar RealScalar
 
typedef MatrixType::Scalar Scalar
 
typedef Matrix< Scalar, Dynamic, 1 > ScalarVector
 
typedef MatrixType::StorageIndex StorageIndex
 

Public Member Functions

template<typename Rhs , typename Dest >
bool _solve_impl (const MatrixBase< Rhs > &B, MatrixBase< Dest > &dest) const
 
void _sort_matrix_Q ()
 
void analyzePattern (const MatrixType &mat)
 Preprocessing step of a QR factorization. More...
 
Index cols () const
 
const PermutationTypecolsPermutation () const
 
void compute (const MatrixType &mat)
 
void factorize (const MatrixType &mat)
 Performs the numerical QR factorization of the input matrix. More...
 
ComputationInfo info () const
 Reports whether previous computation was successful. More...
 
std::string lastErrorMessage () const
 
SparseQRMatrixQReturnType< SparseQRmatrixQ () const
 
const QRMatrixTypematrixR () const
 
Index rank () const
 
Index rows () const
 
void setPivotThreshold (const RealScalar &threshold)
 
template<typename Rhs >
const Solve< SparseQR, Rhs > solve (const MatrixBase< Rhs > &B) const
 
template<typename Rhs >
const Solve< SparseQR, Rhs > solve (const SparseMatrixBase< Rhs > &B) const
 
 SparseQR ()
 
 SparseQR (const MatrixType &mat)
 
- Public Member Functions inherited from Eigen::SparseSolverBase< SparseQR< MatrixType_, OrderingType_ > >
SparseQR< MatrixType_, OrderingType_ > & derived ()
 
const SparseQR< MatrixType_, OrderingType_ > & derived () const
 
const Solve< SparseQR< MatrixType_, OrderingType_ >, Rhs > solve (const MatrixBase< Rhs > &b) const
 
const Solve< SparseQR< MatrixType_, OrderingType_ >, Rhs > solve (const SparseMatrixBase< Rhs > &b) const
 
 SparseSolverBase ()
 
 SparseSolverBase (SparseSolverBase &&other)
 
 ~SparseSolverBase ()
 

Protected Types

typedef SparseSolverBase< SparseQR< MatrixType_, OrderingType_ > > Base
 

Protected Attributes

bool m_analysisIsok
 
IndexVector m_etree
 
bool m_factorizationIsok
 
IndexVector m_firstRowElt
 
ScalarVector m_hcoeffs
 
ComputationInfo m_info
 
bool m_isEtreeOk
 
bool m_isQSorted
 
std::string m_lastError
 
Index m_nonzeropivots
 
PermutationType m_outputPerm_c
 
PermutationType m_perm_c
 
PermutationType m_pivotperm
 
QRMatrixType m_pmat
 
QRMatrixType m_Q
 
QRMatrixType m_R
 
RealScalar m_threshold
 
bool m_useDefaultThreshold
 
- Protected Attributes inherited from Eigen::SparseSolverBase< SparseQR< MatrixType_, OrderingType_ > >
bool m_isInitialized
 

Detailed Description

template<typename MatrixType_, typename OrderingType_>
class Eigen::SparseQR< MatrixType_, OrderingType_ >

Sparse left-looking QR factorization with numerical column pivoting.

This class implements a left-looking QR decomposition of sparse matrices with numerical column pivoting. When a column has a norm less than a given tolerance it is implicitly permuted to the end. The QR factorization thus obtained is given by A*P = Q*R where R is upper triangular or trapezoidal.

P is the column permutation which is the product of the fill-reducing and the numerical permutations. Use colsPermutation() to get it.

Q is the orthogonal matrix represented as products of Householder reflectors. Use matrixQ() to get an expression and matrixQ().adjoint() to get the adjoint. You can then apply it to a vector.

R is the sparse triangular or trapezoidal matrix. The later occurs when A is rank-deficient. matrixR().topLeftCorner(rank(), rank()) always returns a triangular factor of full rank.

Template Parameters
MatrixType_The type of the sparse matrix A, must be a column-major SparseMatrix<>
OrderingType_The fill-reducing ordering method. See the OrderingMethods module for the list of built-in and external ordering methods.

This class follows the sparse solver concept .

The numerical pivoting strategy and default threshold are the same as in SuiteSparse QR, and detailed in the following paper: Tim Davis, "Algorithm 915, SuiteSparseQR: Multifrontal Multithreaded Rank-Revealing Sparse QR Factorization, ACM Trans. on Math. Soft. 38(1), 2011. Even though it is qualified as "rank-revealing", this strategy might fail for some rank deficient problems. When this class is used to solve linear or least-square problems it is thus strongly recommended to check the accuracy of the computed solution. If it failed, it usually helps to increase the threshold with setPivotThreshold.

Warning
The input sparse matrix A must be in compressed mode (see SparseMatrix::makeCompressed()).
For complex matrices matrixQ().transpose() will actually return the adjoint matrix.

Definition at line 86 of file SparseQR.h.

Member Typedef Documentation

◆ Base

template<typename MatrixType_ , typename OrderingType_ >
typedef SparseSolverBase<SparseQR<MatrixType_,OrderingType_> > Eigen::SparseQR< MatrixType_, OrderingType_ >::Base
protected

Definition at line 89 of file SparseQR.h.

◆ IndexVector

template<typename MatrixType_ , typename OrderingType_ >
typedef Matrix<StorageIndex, Dynamic, 1> Eigen::SparseQR< MatrixType_, OrderingType_ >::IndexVector

Definition at line 99 of file SparseQR.h.

◆ MatrixType

template<typename MatrixType_ , typename OrderingType_ >
typedef MatrixType_ Eigen::SparseQR< MatrixType_, OrderingType_ >::MatrixType

Definition at line 93 of file SparseQR.h.

◆ OrderingType

template<typename MatrixType_ , typename OrderingType_ >
typedef OrderingType_ Eigen::SparseQR< MatrixType_, OrderingType_ >::OrderingType

Definition at line 94 of file SparseQR.h.

◆ PermutationType

template<typename MatrixType_ , typename OrderingType_ >
typedef PermutationMatrix<Dynamic, Dynamic, StorageIndex> Eigen::SparseQR< MatrixType_, OrderingType_ >::PermutationType

Definition at line 101 of file SparseQR.h.

◆ QRMatrixType

template<typename MatrixType_ , typename OrderingType_ >
typedef SparseMatrix<Scalar,ColMajor,StorageIndex> Eigen::SparseQR< MatrixType_, OrderingType_ >::QRMatrixType

Definition at line 98 of file SparseQR.h.

◆ RealScalar

template<typename MatrixType_ , typename OrderingType_ >
typedef MatrixType::RealScalar Eigen::SparseQR< MatrixType_, OrderingType_ >::RealScalar

Definition at line 96 of file SparseQR.h.

◆ Scalar

template<typename MatrixType_ , typename OrderingType_ >
typedef MatrixType::Scalar Eigen::SparseQR< MatrixType_, OrderingType_ >::Scalar

Definition at line 95 of file SparseQR.h.

◆ ScalarVector

template<typename MatrixType_ , typename OrderingType_ >
typedef Matrix<Scalar, Dynamic, 1> Eigen::SparseQR< MatrixType_, OrderingType_ >::ScalarVector

Definition at line 100 of file SparseQR.h.

◆ StorageIndex

template<typename MatrixType_ , typename OrderingType_ >
typedef MatrixType::StorageIndex Eigen::SparseQR< MatrixType_, OrderingType_ >::StorageIndex

Definition at line 97 of file SparseQR.h.

Member Enumeration Documentation

◆ anonymous enum

template<typename MatrixType_ , typename OrderingType_ >
anonymous enum
Enumerator
ColsAtCompileTime 
MaxColsAtCompileTime 

Definition at line 103 of file SparseQR.h.

103  {
104  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
105  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
106  };
@ MaxColsAtCompileTime
Definition: SparseQR.h:105

Constructor & Destructor Documentation

◆ SparseQR() [1/2]

template<typename MatrixType_ , typename OrderingType_ >
Eigen::SparseQR< MatrixType_, OrderingType_ >::SparseQR ( )
inline

Definition at line 109 of file SparseQR.h.

109  : m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false),m_isEtreeOk(false)
110  { }
bool m_isQSorted
Definition: SparseQR.h:305
bool m_isEtreeOk
Definition: SparseQR.h:306
std::string m_lastError
Definition: SparseQR.h:292
bool m_useDefaultThreshold
Definition: SparseQR.h:301
bool m_analysisIsok
Definition: SparseQR.h:289

◆ SparseQR() [2/2]

template<typename MatrixType_ , typename OrderingType_ >
Eigen::SparseQR< MatrixType_, OrderingType_ >::SparseQR ( const MatrixType mat)
inlineexplicit

Construct a QR factorization of the matrix mat.

Warning
The matrix mat must be in compressed mode (see SparseMatrix::makeCompressed()).
See also
compute()

Definition at line 118 of file SparseQR.h.

118  : m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false),m_isEtreeOk(false)
119  {
120  compute(mat);
121  }
void compute(const MatrixType &mat)
Definition: SparseQR.h:129

Member Function Documentation

◆ _solve_impl()

template<typename MatrixType_ , typename OrderingType_ >
template<typename Rhs , typename Dest >
bool Eigen::SparseQR< MatrixType_, OrderingType_ >::_solve_impl ( const MatrixBase< Rhs > &  B,
MatrixBase< Dest > &  dest 
) const
inline

Definition at line 207 of file SparseQR.h.

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  }
Array< int, 3, 1 > b
MatrixXf B
#define eigen_assert(x)
Definition: Macros.h:902
PermutationType m_perm_c
Definition: SparseQR.h:297
Index cols() const
Definition: SparseQR.h:143
const QRMatrixType & matrixR() const
Definition: SparseQR.h:158
ComputationInfo m_info
Definition: SparseQR.h:291
Index rows() const
Definition: SparseQR.h:139
const PermutationType & colsPermutation() const
Definition: SparseQR.h:194
SparseQRMatrixQReturnType< SparseQR > matrixQ() const
Definition: SparseQR.h:188
Index rank() const
Definition: SparseQR.h:164
@ Success
Definition: Constants.h:446
const Scalar & y
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:82

◆ _sort_matrix_Q()

template<typename MatrixType_ , typename OrderingType_ >
void Eigen::SparseQR< MatrixType_, OrderingType_ >::_sort_matrix_Q ( )
inline

Definition at line 278 of file SparseQR.h.

279  {
280  if(this->m_isQSorted) return;
281  // The matrix Q is sorted during the transposition
282  SparseMatrix<Scalar, RowMajor, Index> mQrm(this->m_Q);
283  this->m_Q = mQrm;
284  this->m_isQSorted = true;
285  }
QRMatrixType m_Q
Definition: SparseQR.h:295

◆ analyzePattern()

template<typename MatrixType , typename OrderingType >
void Eigen::SparseQR< MatrixType, OrderingType >::analyzePattern ( const MatrixType mat)

Preprocessing step of a QR factorization.

Warning
The matrix mat must be in compressed mode (see SparseMatrix::makeCompressed()).

In this step, the fill-reducing permutation is computed and applied to the columns of A and the column elimination tree is computed as well. Only the sparsity pattern of mat is exploited.

Note
In this step it is assumed that there is no empty row in the matrix mat.

Definition at line 322 of file SparseQR.h.

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
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 }
Matrix3f m
int n
void resize(Index newSize)
InverseReturnType inverse() const
const IndicesType & indices() const
constexpr void resize(Index rows, Index cols)
void resize(Index rows, Index cols)
Definition: SparseMatrix.h:739
void reserve(Index reserveSize)
Definition: SparseMatrix.h:300
PermutationType m_outputPerm_c
Definition: SparseQR.h:299
IndexVector m_firstRowElt
Definition: SparseQR.h:304
ScalarVector m_hcoeffs
Definition: SparseQR.h:296
OrderingType_ OrderingType
Definition: SparseQR.h:94
QRMatrixType m_R
Definition: SparseQR.h:294
IndexVector m_etree
Definition: SparseQR.h:303
MatrixType::StorageIndex StorageIndex
Definition: SparseQR.h:97
bfloat16() min(const bfloat16 &a, const bfloat16 &b)
Definition: BFloat16.h:684
int coletree(const MatrixType &mat, IndexVector &parent, IndexVector &firstRowElt, typename MatrixType::StorageIndex *perm=0)

◆ cols()

template<typename MatrixType_ , typename OrderingType_ >
Index Eigen::SparseQR< MatrixType_, OrderingType_ >::cols ( void  ) const
inline
Returns
the number of columns of the represented matrix.

Definition at line 143 of file SparseQR.h.

143 { return m_pmat.cols();}
Index cols() const
Definition: SparseMatrix.h:167
QRMatrixType m_pmat
Definition: SparseQR.h:293

◆ colsPermutation()

template<typename MatrixType_ , typename OrderingType_ >
const PermutationType& Eigen::SparseQR< MatrixType_, OrderingType_ >::colsPermutation ( ) const
inline
Returns
a const reference to the column permutation P that was applied to A such that A*P = Q*R It is the combination of the fill-in reducing permutation and numerical column pivoting.

Definition at line 194 of file SparseQR.h.

195  {
196  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
197  return m_outputPerm_c;
198  }

◆ compute()

template<typename MatrixType_ , typename OrderingType_ >
void Eigen::SparseQR< MatrixType_, OrderingType_ >::compute ( const MatrixType mat)
inline

Computes the QR factorization of the sparse matrix mat.

Warning
The matrix mat must be in compressed mode (see SparseMatrix::makeCompressed()).
See also
analyzePattern(), factorize()

Definition at line 129 of file SparseQR.h.

130  {
132  factorize(mat);
133  }
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

◆ factorize()

template<typename MatrixType , typename OrderingType >
void Eigen::SparseQR< MatrixType, OrderingType >::factorize ( const MatrixType mat)

Performs the numerical QR factorization of the input matrix.

The function SparseQR::analyzePattern(const MatrixType&) must have been called beforehand with a matrix having the same sparsity pattern than mat.

Parameters
matThe sparse column-major matrix

Definition at line 363 of file SparseQR.h.

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  {
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  */
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
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";
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++)
581 
582  // Recompute the column elimination tree
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();
593  m_R.finalize();
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
607  }
608 
609  m_isInitialized = true;
610  m_factorizationIsok = true;
611  m_info = Success;
612 }
const AbsReturnType abs() const
const SqrtReturnType sqrt() const
ColXpr col(Index i)
This is the const version of col().
const ImagReturnType imag() const
RealReturnType real() const
float * p
static ConstMapType Map(const Scalar *data)
Derived & setZero(Index size)
Scalar dot(const MatrixBase< OtherDerived > &other) const
RealScalar norm() const
Definition: SparseDot.h:86
void startVec(Index outer)
Definition: SparseMatrix.h:447
Scalar & insertBackByOuterInner(Index outer, Index inner)
Definition: SparseMatrix.h:425
const StorageIndex * innerNonZeroPtr() const
Definition: SparseMatrix.h:204
const StorageIndex * outerIndexPtr() const
Definition: SparseMatrix.h:195
Scalar & insertBackByOuterInnerUnordered(Index outer, Index inner)
Definition: SparseMatrix.h:437
MatrixType::Scalar Scalar
Definition: SparseQR.h:95
RealScalar m_threshold
Definition: SparseQR.h:300
bool m_factorizationIsok
Definition: SparseQR.h:290
Matrix< StorageIndex, Dynamic, 1 > IndexVector
Definition: SparseQR.h:99
Matrix< Scalar, Dynamic, 1 > ScalarVector
Definition: SparseQR.h:100
PermutationType m_pivotperm
Definition: SparseQR.h:298
MatrixType::RealScalar RealScalar
Definition: SparseQR.h:96
SparseMatrix< Scalar, ColMajor, StorageIndex > QRMatrixType
Definition: SparseQR.h:98
Index m_nonzeropivots
Definition: SparseQR.h:302
@ InvalidInput
Definition: Constants.h:453
bfloat16() max(const bfloat16 &a, const bfloat16 &b)
Definition: BFloat16.h:690
void swap(scoped_array< T > &a, scoped_array< T > &b)
Definition: Memory.h:788
bool abs2(bool x)
EIGEN_ALWAYS_INLINE T maxi(const T &x, const T &y)
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_conjugate_op< typename Derived::Scalar >, const Derived > conj(const Eigen::ArrayBase< Derived > &x)
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)
std::ptrdiff_t j

◆ info()

template<typename MatrixType_ , typename OrderingType_ >
ComputationInfo Eigen::SparseQR< MatrixType_, OrderingType_ >::info ( ) const
inline

Reports whether previous computation was successful.

Returns
Success if computation was successful, NumericalIssue if the QR factorization reports a numerical problem InvalidInput if the input matrix is invalid
See also
iparm()

Definition at line 270 of file SparseQR.h.

271  {
272  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
273  return m_info;
274  }

◆ lastErrorMessage()

template<typename MatrixType_ , typename OrderingType_ >
std::string Eigen::SparseQR< MatrixType_, OrderingType_ >::lastErrorMessage ( ) const
inline
Returns
A string describing the type of error. This method is provided to ease debugging, not to handle errors.

Definition at line 203 of file SparseQR.h.

203 { return m_lastError; }

◆ matrixQ()

template<typename MatrixType_ , typename OrderingType_ >
SparseQRMatrixQReturnType<SparseQR> Eigen::SparseQR< MatrixType_, OrderingType_ >::matrixQ ( void  ) const
inline
Returns
an expression of the matrix Q as products of sparse Householder reflectors. The common usage of this function is to apply it to a dense matrix or vector
VectorXd B1, B2;
// Initialize B1
B2 = matrixQ() * B1;
Matrix< double, Dynamic, 1 > VectorXd
DynamicĂ—1 vector of type double.
Definition: Matrix.h:502

To get a plain SparseMatrix representation of Q:

SparseMatrix<double> Q;
Q = SparseQR<SparseMatrix<double> >(A).matrixQ();
MatrixXcf A
MatrixXf Q

Internally, this call simply performs a sparse product between the matrix Q and a sparse identity matrix. However, due to the fact that the sparse reflectors are stored unsorted, two transpositions are needed to sort them before performing the product.

Definition at line 188 of file SparseQR.h.

189  { return SparseQRMatrixQReturnType<SparseQR>(*this); }

◆ matrixR()

template<typename MatrixType_ , typename OrderingType_ >
const QRMatrixType& Eigen::SparseQR< MatrixType_, OrderingType_ >::matrixR ( ) const
inline
Returns
a const reference to the sparse upper triangular matrix R of the QR factorization.
Warning
The entries of the returned matrix are not sorted. This means that using it in algorithms expecting sorted entries will fail. This include random coefficient accesses (SpaseMatrix::coeff()), and coefficient-wise operations. Matrix products and triangular solves are fine though.

To sort the entries, you can assign it to a row-major matrix, and if a column-major matrix is required, you can copy it again:

SparseMatrix<double> R = qr.matrixR(); // column-major, not sorted!
SparseMatrix<double,RowMajor> Rr = qr.matrixR(); // row-major, sorted
SparseMatrix<double> Rc = Rr; // column-major, sorted
HouseholderQR< MatrixXf > qr(A)

Definition at line 158 of file SparseQR.h.

158 { return m_R; }

◆ rank()

template<typename MatrixType_ , typename OrderingType_ >
Index Eigen::SparseQR< MatrixType_, OrderingType_ >::rank ( ) const
inline
Returns
the number of non linearly dependent columns as determined by the pivoting threshold.
See also
setPivotThreshold()

Definition at line 164 of file SparseQR.h.

165  {
166  eigen_assert(m_isInitialized && "The factorization should be called first, use compute()");
167  return m_nonzeropivots;
168  }

◆ rows()

template<typename MatrixType_ , typename OrderingType_ >
Index Eigen::SparseQR< MatrixType_, OrderingType_ >::rows ( void  ) const
inline
Returns
the number of rows of the represented matrix.

Definition at line 139 of file SparseQR.h.

139 { return m_pmat.rows(); }
Index rows() const
Definition: SparseMatrix.h:165

◆ setPivotThreshold()

template<typename MatrixType_ , typename OrderingType_ >
void Eigen::SparseQR< MatrixType_, OrderingType_ >::setPivotThreshold ( const RealScalar threshold)
inline

Sets the threshold that is used to determine linearly dependent columns during the factorization.

In practice, if during the factorization the norm of the column that has to be eliminated is below this threshold, then the entire column is treated as zero, and it is moved at the end.

Definition at line 237 of file SparseQR.h.

238  {
239  m_useDefaultThreshold = false;
240  m_threshold = threshold;
241  }

◆ solve() [1/2]

template<typename MatrixType_ , typename OrderingType_ >
template<typename Rhs >
const Solve<SparseQR, Rhs> Eigen::SparseQR< MatrixType_, OrderingType_ >::solve ( const MatrixBase< Rhs > &  B) const
inline
Returns
the solution X of \( A X = B \) using the current decomposition of A.
See also
compute()

Definition at line 248 of file SparseQR.h.

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  }

◆ solve() [2/2]

template<typename MatrixType_ , typename OrderingType_ >
template<typename Rhs >
const Solve<SparseQR, Rhs> Eigen::SparseQR< MatrixType_, OrderingType_ >::solve ( const SparseMatrixBase< Rhs > &  B) const
inline

Definition at line 255 of file SparseQR.h.

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  }

Member Data Documentation

◆ m_analysisIsok

template<typename MatrixType_ , typename OrderingType_ >
bool Eigen::SparseQR< MatrixType_, OrderingType_ >::m_analysisIsok
protected

Definition at line 289 of file SparseQR.h.

◆ m_etree

template<typename MatrixType_ , typename OrderingType_ >
IndexVector Eigen::SparseQR< MatrixType_, OrderingType_ >::m_etree
protected

Definition at line 303 of file SparseQR.h.

◆ m_factorizationIsok

template<typename MatrixType_ , typename OrderingType_ >
bool Eigen::SparseQR< MatrixType_, OrderingType_ >::m_factorizationIsok
protected

Definition at line 290 of file SparseQR.h.

◆ m_firstRowElt

template<typename MatrixType_ , typename OrderingType_ >
IndexVector Eigen::SparseQR< MatrixType_, OrderingType_ >::m_firstRowElt
protected

Definition at line 304 of file SparseQR.h.

◆ m_hcoeffs

template<typename MatrixType_ , typename OrderingType_ >
ScalarVector Eigen::SparseQR< MatrixType_, OrderingType_ >::m_hcoeffs
protected

Definition at line 296 of file SparseQR.h.

◆ m_info

template<typename MatrixType_ , typename OrderingType_ >
ComputationInfo Eigen::SparseQR< MatrixType_, OrderingType_ >::m_info
mutableprotected

Definition at line 291 of file SparseQR.h.

◆ m_isEtreeOk

template<typename MatrixType_ , typename OrderingType_ >
bool Eigen::SparseQR< MatrixType_, OrderingType_ >::m_isEtreeOk
protected

Definition at line 306 of file SparseQR.h.

◆ m_isQSorted

template<typename MatrixType_ , typename OrderingType_ >
bool Eigen::SparseQR< MatrixType_, OrderingType_ >::m_isQSorted
protected

Definition at line 305 of file SparseQR.h.

◆ m_lastError

template<typename MatrixType_ , typename OrderingType_ >
std::string Eigen::SparseQR< MatrixType_, OrderingType_ >::m_lastError
protected

Definition at line 292 of file SparseQR.h.

◆ m_nonzeropivots

template<typename MatrixType_ , typename OrderingType_ >
Index Eigen::SparseQR< MatrixType_, OrderingType_ >::m_nonzeropivots
protected

Definition at line 302 of file SparseQR.h.

◆ m_outputPerm_c

template<typename MatrixType_ , typename OrderingType_ >
PermutationType Eigen::SparseQR< MatrixType_, OrderingType_ >::m_outputPerm_c
protected

Definition at line 299 of file SparseQR.h.

◆ m_perm_c

template<typename MatrixType_ , typename OrderingType_ >
PermutationType Eigen::SparseQR< MatrixType_, OrderingType_ >::m_perm_c
protected

Definition at line 297 of file SparseQR.h.

◆ m_pivotperm

template<typename MatrixType_ , typename OrderingType_ >
PermutationType Eigen::SparseQR< MatrixType_, OrderingType_ >::m_pivotperm
protected

Definition at line 298 of file SparseQR.h.

◆ m_pmat

template<typename MatrixType_ , typename OrderingType_ >
QRMatrixType Eigen::SparseQR< MatrixType_, OrderingType_ >::m_pmat
protected

Definition at line 293 of file SparseQR.h.

◆ m_Q

template<typename MatrixType_ , typename OrderingType_ >
QRMatrixType Eigen::SparseQR< MatrixType_, OrderingType_ >::m_Q
protected

Definition at line 295 of file SparseQR.h.

◆ m_R

template<typename MatrixType_ , typename OrderingType_ >
QRMatrixType Eigen::SparseQR< MatrixType_, OrderingType_ >::m_R
protected

Definition at line 294 of file SparseQR.h.

◆ m_threshold

template<typename MatrixType_ , typename OrderingType_ >
RealScalar Eigen::SparseQR< MatrixType_, OrderingType_ >::m_threshold
protected

Definition at line 300 of file SparseQR.h.

◆ m_useDefaultThreshold

template<typename MatrixType_ , typename OrderingType_ >
bool Eigen::SparseQR< MatrixType_, OrderingType_ >::m_useDefaultThreshold
protected

Definition at line 301 of file SparseQR.h.


The documentation for this class was generated from the following file: