Prusa Slicer 2.6.0
Loading...
Searching...
No Matches
Eigen::SparseQR< _MatrixType, _OrderingType > Class Template Reference

Sparse left-looking rank-revealing QR factorization. More...

#include <src/eigen/Eigen/src/SparseQR/SparseQR.h>

+ Inheritance diagram for Eigen::SparseQR< _MatrixType, _OrderingType >:
+ Collaboration diagram for Eigen::SparseQR< _MatrixType, _OrderingType >:

Public Types

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

Public Member Functions

 SparseQR ()
 
 SparseQR (const MatrixType &mat)
 
void compute (const MatrixType &mat)
 
void analyzePattern (const MatrixType &mat)
 Preprocessing step of a QR factorization.
 
void factorize (const MatrixType &mat)
 Performs the numerical QR factorization of the input matrix.
 
Index rows () const
 
Index cols () const
 
const QRMatrixTypematrixR () const
 
Index rank () const
 
SparseQRMatrixQReturnType< SparseQRmatrixQ () const
 
const PermutationTypecolsPermutation () const
 
std::string lastErrorMessage () const
 
template<typename Rhs , typename Dest >
bool _solve_impl (const MatrixBase< Rhs > &B, MatrixBase< Dest > &dest) 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
 
ComputationInfo info () const
 Reports whether previous computation was successful.
 
void _sort_matrix_Q ()
 
SparseQR< _MatrixType, _OrderingType > & derived ()
 
const SparseQR< _MatrixType, _OrderingType > & derived () const
 
void _solve_impl (const SparseMatrixBase< Rhs > &b, SparseMatrixBase< Dest > &dest) const
 

Protected Types

typedef SparseSolverBase< SparseQR< _MatrixType, _OrderingType > > Base
 

Protected Attributes

bool m_analysisIsok
 
bool m_factorizationIsok
 
ComputationInfo m_info
 
std::string m_lastError
 
QRMatrixType m_pmat
 
QRMatrixType m_R
 
QRMatrixType m_Q
 
ScalarVector m_hcoeffs
 
PermutationType m_perm_c
 
PermutationType m_pivotperm
 
PermutationType m_outputPerm_c
 
RealScalar m_threshold
 
bool m_useDefaultThreshold
 
Index m_nonzeropivots
 
IndexVector m_etree
 
IndexVector m_firstRowElt
 
bool m_isQSorted
 
bool m_isEtreeOk
 
bool m_isInitialized
 

Friends

template<typename , typename >
struct SparseQR_QProduct
 

Detailed Description

template<typename _MatrixType, typename _OrderingType>
class Eigen::SparseQR< _MatrixType, _OrderingType >

Sparse left-looking rank-revealing QR factorization.

This class implements a left-looking rank-revealing QR decomposition of sparse matrices. 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 rank-revealing 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
_MatrixTypeThe type of the sparse matrix A, must be a column-major SparseMatrix<>
_OrderingTypeThe fill-reducing ordering method. See the OrderingMethods module for the list of built-in and external ordering methods.

\implsparsesolverconcept

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.

Member Typedef Documentation

◆ Base

template<typename _MatrixType , typename _OrderingType >
typedef SparseSolverBase<SparseQR<_MatrixType,_OrderingType> > Eigen::SparseQR< _MatrixType, _OrderingType >::Base
protected

◆ IndexVector

template<typename _MatrixType , typename _OrderingType >
typedef Matrix<StorageIndex, Dynamic, 1> Eigen::SparseQR< _MatrixType, _OrderingType >::IndexVector

◆ MatrixType

template<typename _MatrixType , typename _OrderingType >
typedef _MatrixType Eigen::SparseQR< _MatrixType, _OrderingType >::MatrixType

◆ OrderingType

template<typename _MatrixType , typename _OrderingType >
typedef _OrderingType Eigen::SparseQR< _MatrixType, _OrderingType >::OrderingType

◆ PermutationType

template<typename _MatrixType , typename _OrderingType >
typedef PermutationMatrix<Dynamic, Dynamic, StorageIndex> Eigen::SparseQR< _MatrixType, _OrderingType >::PermutationType

◆ QRMatrixType

template<typename _MatrixType , typename _OrderingType >
typedef SparseMatrix<Scalar,ColMajor,StorageIndex> Eigen::SparseQR< _MatrixType, _OrderingType >::QRMatrixType

◆ RealScalar

template<typename _MatrixType , typename _OrderingType >
typedef MatrixType::RealScalar Eigen::SparseQR< _MatrixType, _OrderingType >::RealScalar

◆ Scalar

template<typename _MatrixType , typename _OrderingType >
typedef MatrixType::Scalar Eigen::SparseQR< _MatrixType, _OrderingType >::Scalar

◆ ScalarVector

template<typename _MatrixType , typename _OrderingType >
typedef Matrix<Scalar, Dynamic, 1> Eigen::SparseQR< _MatrixType, _OrderingType >::ScalarVector

◆ StorageIndex

template<typename _MatrixType , typename _OrderingType >
typedef MatrixType::StorageIndex Eigen::SparseQR< _MatrixType, _OrderingType >::StorageIndex

Member Enumeration Documentation

◆ anonymous enum

template<typename _MatrixType , typename _OrderingType >
anonymous enum
Enumerator
ColsAtCompileTime 
MaxColsAtCompileTime 
89 {
90 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
91 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
92 };
@ MaxColsAtCompileTime
Definition SparseQR.h:91
@ ColsAtCompileTime
Definition SparseQR.h:90

Constructor & Destructor Documentation

◆ SparseQR() [1/2]

template<typename _MatrixType , typename _OrderingType >
Eigen::SparseQR< _MatrixType, _OrderingType >::SparseQR ( )
inline
96 { }
std::string m_lastError
Definition SparseQR.h:278
bool m_isEtreeOk
Definition SparseQR.h:292
bool m_useDefaultThreshold
Definition SparseQR.h:287
bool m_analysisIsok
Definition SparseQR.h:275
bool m_isQSorted
Definition SparseQR.h:291

◆ 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()
105 {
106 compute(mat);
107 }
void compute(const MatrixType &mat)
Definition SparseQR.h:115

References Eigen::SparseQR< _MatrixType, _OrderingType >::compute().

+ Here is the call graph for this function:

Member Function Documentation

◆ _solve_impl() [1/2]

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
194 {
195 eigen_assert(m_isInitialized && "The factorization should be called first, use compute()");
196 eigen_assert(this->rows() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix");
197
198 Index rank = this->rank();
199
200 // Compute Q^* * b;
201 typename Dest::PlainObject y, b;
202 y = this->matrixQ().adjoint() * B;
203 b = y;
204
205 // Solve with the triangular matrix R
206 y.resize((std::max<Index>)(cols(),y.rows()),y.cols());
207 y.topRows(rank) = this->matrixR().topLeftCorner(rank, rank).template triangularView<Upper>().solve(b.topRows(rank));
208 y.bottomRows(y.rows()-rank).setZero();
209
210 // Apply the column permutation
211 if (m_perm_c.size()) dest = colsPermutation() * y.topRows(cols());
212 else dest = y.topRows(cols());
213
214 m_info = Success;
215 return true;
216 }
#define eigen_assert(x)
Definition Macros.h:579
Index size() const
Definition PermutationMatrix.h:108
ComputationInfo m_info
Definition SparseQR.h:277
PermutationType m_perm_c
Definition SparseQR.h:283
Index cols() const
Definition SparseQR.h:129
const PermutationType & colsPermutation() const
Definition SparseQR.h:180
Index rank() const
Definition SparseQR.h:150
SparseQRMatrixQReturnType< SparseQR > matrixQ() const
Definition SparseQR.h:174
const QRMatrixType & matrixR() const
Definition SparseQR.h:144
Index rows() const
Definition SparseQR.h:125
@ Success
Definition Constants.h:432
const Scalar & y
Definition MathFunctions.h:552
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition Meta.h:33

References Eigen::SparseQR< _MatrixType, _OrderingType >::cols(), Eigen::SparseQR< _MatrixType, _OrderingType >::colsPermutation(), eigen_assert, Eigen::SparseQR< _MatrixType, _OrderingType >::m_info, Eigen::SparseSolverBase< SparseQR< _MatrixType, _OrderingType > >::m_isInitialized, Eigen::SparseQR< _MatrixType, _OrderingType >::m_perm_c, Eigen::SparseQR< _MatrixType, _OrderingType >::matrixQ(), Eigen::SparseQR< _MatrixType, _OrderingType >::matrixR(), Eigen::SparseQR< _MatrixType, _OrderingType >::rank(), Eigen::SparseQR< _MatrixType, _OrderingType >::rows(), Eigen::PermutationBase< Derived >::size(), and Eigen::Success.

+ Here is the call graph for this function:

◆ _solve_impl() [2/2]

void Eigen::SparseSolverBase< SparseQR< _MatrixType, _OrderingType > >::_solve_impl ( const SparseMatrixBase< Rhs > &  b,
SparseMatrixBase< Dest > &  dest 
) const
inlineinherited
112 {
113 internal::solve_sparse_through_dense_panels(derived(), b.derived(), dest.derived());
114 }
SparseQR< _MatrixType, _OrderingType > & derived()
Definition SparseSolverBase.h:79

◆ _sort_matrix_Q()

template<typename _MatrixType , typename _OrderingType >
void Eigen::SparseQR< _MatrixType, _OrderingType >::_sort_matrix_Q ( )
inline
265 {
266 if(this->m_isQSorted) return;
267 // The matrix Q is sorted during the transposition
268 SparseMatrix<Scalar, RowMajor, Index> mQrm(this->m_Q);
269 this->m_Q = mQrm;
270 this->m_isQSorted = true;
271 }
QRMatrixType m_Q
Definition SparseQR.h:281

References Eigen::SparseQR< _MatrixType, _OrderingType >::m_isQSorted, and Eigen::SparseQR< _MatrixType, _OrderingType >::m_Q.

◆ 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.
309{
310 eigen_assert(mat.isCompressed() && "SparseQR requires a sparse matrix in compressed mode. Call .makeCompressed() before passing it to SparseQR");
311 // Copy to a column major matrix if the input is rowmajor
313 // Compute the column fill reducing ordering
314 OrderingType ord;
315 ord(matCpy, m_perm_c);
316 Index n = mat.cols();
317 Index m = mat.rows();
318 Index diagSize = (std::min)(m,n);
319
320 if (!m_perm_c.size())
321 {
322 m_perm_c.resize(n);
323 m_perm_c.indices().setLinSpaced(n, 0,StorageIndex(n-1));
324 }
325
326 // Compute the column elimination tree of the permuted matrix
329 m_isEtreeOk = true;
330
331 m_R.resize(m, n);
332 m_Q.resize(m, diagSize);
333
334 // Allocate space for nonzero elements : rough estimation
335 m_R.reserve(2*mat.nonZeros()); //FIXME Get a more accurate estimation through symbolic factorization with the etree
336 m_Q.reserve(2*mat.nonZeros());
337 m_hcoeffs.resize(diagSize);
338 m_analysisIsok = true;
339}
void resize(Index newSize)
Definition PermutationMatrix.h:136
InverseReturnType inverse() const
Definition PermutationMatrix.h:196
const IndicesType & indices() const
Definition PermutationMatrix.h:388
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void resize(Index rows, Index cols)
Definition PlainObjectBase.h:279
void reserve(Index reserveSize)
Definition SparseMatrix.h:262
void resize(Index rows, Index cols)
Definition SparseMatrix.h:621
IndexVector m_firstRowElt
Definition SparseQR.h:290
MatrixType::StorageIndex StorageIndex
Definition SparseQR.h:83
IndexVector m_etree
Definition SparseQR.h:289
PermutationType m_outputPerm_c
Definition SparseQR.h:285
_OrderingType OrderingType
Definition SparseQR.h:80
QRMatrixType m_R
Definition SparseQR.h:280
ScalarVector m_hcoeffs
Definition SparseQR.h:282
Then type
Definition Meta.h:58
int coletree(const MatrixType &mat, IndexVector &parent, IndexVector &firstRowElt, typename MatrixType::StorageIndex *perm=0)
Definition SparseColEtree.h:61

References Eigen::internal::coletree(), and eigen_assert.

Referenced by Eigen::SparseQR< _MatrixType, _OrderingType >::compute().

+ Here is the call graph for this function:
+ Here is the caller graph for this function:

◆ cols()

template<typename _MatrixType , typename _OrderingType >
Index Eigen::SparseQR< _MatrixType, _OrderingType >::cols ( ) const
inline
Returns
the number of columns of the represented matrix.
129{ return m_pmat.cols();}
Index cols() const
Definition SparseMatrix.h:138
QRMatrixType m_pmat
Definition SparseQR.h:279

References Eigen::SparseMatrix< _Scalar, _Options, _StorageIndex >::cols(), and Eigen::SparseQR< _MatrixType, _OrderingType >::m_pmat.

Referenced by Eigen::SparseQR< _MatrixType, _OrderingType >::_solve_impl().

+ Here is the call graph for this function:
+ Here is the caller graph for this function:

◆ 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.
181 {
182 eigen_assert(m_isInitialized && "Decomposition is not initialized.");
183 return m_outputPerm_c;
184 }

References eigen_assert, Eigen::SparseSolverBase< SparseQR< _MatrixType, _OrderingType > >::m_isInitialized, and Eigen::SparseQR< _MatrixType, _OrderingType >::m_outputPerm_c.

Referenced by Eigen::SparseQR< _MatrixType, _OrderingType >::_solve_impl().

+ Here is the caller graph for this function:

◆ 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()
116 {
117 analyzePattern(mat);
118 factorize(mat);
119 }
void analyzePattern(const MatrixType &mat)
Preprocessing step of a QR factorization.
Definition SparseQR.h:308
void factorize(const MatrixType &mat)
Performs the numerical QR factorization of the input matrix.
Definition SparseQR.h:349

References Eigen::SparseQR< _MatrixType, _OrderingType >::analyzePattern(), and Eigen::SparseQR< _MatrixType, _OrderingType >::factorize().

Referenced by Eigen::SparseQR< _MatrixType, _OrderingType >::SparseQR().

+ Here is the call graph for this function:
+ Here is the caller graph for this function:

◆ derived() [1/2]

SparseQR< _MatrixType, _OrderingType > & Eigen::SparseSolverBase< SparseQR< _MatrixType, _OrderingType > >::derived ( )
inlineinherited
79{ return *static_cast<Derived*>(this); }

◆ derived() [2/2]

const SparseQR< _MatrixType, _OrderingType > & Eigen::SparseSolverBase< SparseQR< _MatrixType, _OrderingType > >::derived ( ) const
inlineinherited
80{ return *static_cast<const Derived*>(this); }

◆ 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
350{
351 using std::abs;
352
353 eigen_assert(m_analysisIsok && "analyzePattern() should be called before this step");
354 StorageIndex m = StorageIndex(mat.rows());
355 StorageIndex n = StorageIndex(mat.cols());
356 StorageIndex diagSize = (std::min)(m,n);
357 IndexVector mark((std::max)(m,n)); mark.setConstant(-1); // Record the visited nodes
358 IndexVector Ridx(n), Qidx(m); // Store temporarily the row indexes for the current column of R and Q
359 Index nzcolR, nzcolQ; // Number of nonzero for the current column of R and Q
360 ScalarVector tval(m); // The dense vector used to compute the current column
361 RealScalar pivotThreshold = m_threshold;
362
363 m_R.setZero();
364 m_Q.setZero();
365 m_pmat = mat;
366 if(!m_isEtreeOk)
367 {
370 m_isEtreeOk = true;
371 }
372
373 m_pmat.uncompress(); // To have the innerNonZeroPtr allocated
374
375 // Apply the fill-in reducing permutation lazily:
376 {
377 // If the input is row major, copy the original column indices,
378 // otherwise directly use the input matrix
379 //
380 IndexVector originalOuterIndicesCpy;
381 const StorageIndex *originalOuterIndices = mat.outerIndexPtr();
382 if(MatrixType::IsRowMajor)
383 {
384 originalOuterIndicesCpy = IndexVector::Map(m_pmat.outerIndexPtr(),n+1);
385 originalOuterIndices = originalOuterIndicesCpy.data();
386 }
387
388 for (int i = 0; i < n; i++)
389 {
390 Index p = m_perm_c.size() ? m_perm_c.indices()(i) : i;
391 m_pmat.outerIndexPtr()[p] = originalOuterIndices[i];
392 m_pmat.innerNonZeroPtr()[p] = originalOuterIndices[i+1] - originalOuterIndices[i];
393 }
394 }
395
396 /* Compute the default threshold as in MatLab, see:
397 * Tim Davis, "Algorithm 915, SuiteSparseQR: Multifrontal Multithreaded Rank-Revealing
398 * Sparse QR Factorization, ACM Trans. on Math. Soft. 38(1), 2011, Page 8:3
399 */
401 {
402 RealScalar max2Norm = 0.0;
403 for (int j = 0; j < n; j++) max2Norm = numext::maxi(max2Norm, m_pmat.col(j).norm());
404 if(max2Norm==RealScalar(0))
405 max2Norm = RealScalar(1);
406 pivotThreshold = 20 * (m + n) * max2Norm * NumTraits<RealScalar>::epsilon();
407 }
408
409 // Initialize the numerical permutation
411
412 StorageIndex nonzeroCol = 0; // Record the number of valid pivots
413 m_Q.startVec(0);
414
415 // Left looking rank-revealing QR factorization: compute a column of R and Q at a time
416 for (StorageIndex col = 0; col < n; ++col)
417 {
418 mark.setConstant(-1);
420 mark(nonzeroCol) = col;
421 Qidx(0) = nonzeroCol;
422 nzcolR = 0; nzcolQ = 1;
423 bool found_diag = nonzeroCol>=m;
424 tval.setZero();
425
426 // Symbolic factorization: find the nonzero locations of the column k of the factors R and Q, i.e.,
427 // all the nodes (with indexes lower than rank) reachable through the column elimination tree (etree) rooted at node k.
428 // Note: if the diagonal entry does not exist, then its contribution must be explicitly added,
429 // thus the trick with found_diag that permits to do one more iteration on the diagonal element if this one has not been found.
430 for (typename QRMatrixType::InnerIterator itp(m_pmat, col); itp || !found_diag; ++itp)
431 {
432 StorageIndex curIdx = nonzeroCol;
433 if(itp) curIdx = StorageIndex(itp.row());
434 if(curIdx == nonzeroCol) found_diag = true;
435
436 // Get the nonzeros indexes of the current column of R
437 StorageIndex st = m_firstRowElt(curIdx); // The traversal of the etree starts here
438 if (st < 0 )
439 {
440 m_lastError = "Empty row found during numerical factorization";
442 return;
443 }
444
445 // Traverse the etree
446 Index bi = nzcolR;
447 for (; mark(st) != col; st = m_etree(st))
448 {
449 Ridx(nzcolR) = st; // Add this row to the list,
450 mark(st) = col; // and mark this row as visited
451 nzcolR++;
452 }
453
454 // Reverse the list to get the topological ordering
455 Index nt = nzcolR-bi;
456 for(Index i = 0; i < nt/2; i++) std::swap(Ridx(bi+i), Ridx(nzcolR-i-1));
457
458 // Copy the current (curIdx,pcol) value of the input matrix
459 if(itp) tval(curIdx) = itp.value();
460 else tval(curIdx) = Scalar(0);
461
462 // Compute the pattern of Q(:,k)
463 if(curIdx > nonzeroCol && mark(curIdx) != col )
464 {
465 Qidx(nzcolQ) = curIdx; // Add this row to the pattern of Q,
466 mark(curIdx) = col; // and mark it as visited
467 nzcolQ++;
468 }
469 }
470
471 // Browse all the indexes of R(:,col) in reverse order
472 for (Index i = nzcolR-1; i >= 0; i--)
473 {
474 Index curIdx = Ridx(i);
475
476 // Apply the curIdx-th householder vector to the current column (temporarily stored into tval)
477 Scalar tdot(0);
478
479 // First compute q' * tval
480 tdot = m_Q.col(curIdx).dot(tval);
481
482 tdot *= m_hcoeffs(curIdx);
483
484 // Then update tval = tval - q * tau
485 // FIXME: tval -= tdot * m_Q.col(curIdx) should amount to the same (need to check/add support for efficient "dense ?= sparse")
486 for (typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq)
487 tval(itq.row()) -= itq.value() * tdot;
488
489 // Detect fill-in for the current column of Q
490 if(m_etree(Ridx(i)) == nonzeroCol)
491 {
492 for (typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq)
493 {
494 StorageIndex iQ = StorageIndex(itq.row());
495 if (mark(iQ) != col)
496 {
497 Qidx(nzcolQ++) = iQ; // Add this row to the pattern of Q,
498 mark(iQ) = col; // and mark it as visited
499 }
500 }
501 }
502 } // End update current column
503
504 Scalar tau = RealScalar(0);
505 RealScalar beta = 0;
506
507 if(nonzeroCol < diagSize)
508 {
509 // Compute the Householder reflection that eliminate the current column
510 // FIXME this step should call the Householder module.
511 Scalar c0 = nzcolQ ? tval(Qidx(0)) : Scalar(0);
512
513 // First, the squared norm of Q((col+1):m, col)
514 RealScalar sqrNorm = 0.;
515 for (Index itq = 1; itq < nzcolQ; ++itq) sqrNorm += numext::abs2(tval(Qidx(itq)));
516 if(sqrNorm == RealScalar(0) && numext::imag(c0) == RealScalar(0))
517 {
518 beta = numext::real(c0);
519 tval(Qidx(0)) = 1;
520 }
521 else
522 {
523 using std::sqrt;
524 beta = sqrt(numext::abs2(c0) + sqrNorm);
525 if(numext::real(c0) >= RealScalar(0))
526 beta = -beta;
527 tval(Qidx(0)) = 1;
528 for (Index itq = 1; itq < nzcolQ; ++itq)
529 tval(Qidx(itq)) /= (c0 - beta);
530 tau = numext::conj((beta-c0) / beta);
531
532 }
533 }
534
535 // Insert values in R
536 for (Index i = nzcolR-1; i >= 0; i--)
537 {
538 Index curIdx = Ridx(i);
539 if(curIdx < nonzeroCol)
540 {
541 m_R.insertBackByOuterInnerUnordered(col, curIdx) = tval(curIdx);
542 tval(curIdx) = Scalar(0.);
543 }
544 }
545
546 if(nonzeroCol < diagSize && abs(beta) >= pivotThreshold)
547 {
548 m_R.insertBackByOuterInner(col, nonzeroCol) = beta;
549 // The householder coefficient
550 m_hcoeffs(nonzeroCol) = tau;
551 // Record the householder reflections
552 for (Index itq = 0; itq < nzcolQ; ++itq)
553 {
554 Index iQ = Qidx(itq);
555 m_Q.insertBackByOuterInnerUnordered(nonzeroCol,iQ) = tval(iQ);
556 tval(iQ) = Scalar(0.);
557 }
558 nonzeroCol++;
559 if(nonzeroCol<diagSize)
560 m_Q.startVec(nonzeroCol);
561 }
562 else
563 {
564 // Zero pivot found: move implicitly this column to the end
565 for (Index j = nonzeroCol; j < n-1; j++)
566 std::swap(m_pivotperm.indices()(j), m_pivotperm.indices()[j+1]);
567
568 // Recompute the column elimination tree
570 m_isEtreeOk = false;
571 }
572 }
573
574 m_hcoeffs.tail(diagSize-nonzeroCol).setZero();
575
576 // Finalize the column pointers of the sparse matrices R and Q
577 m_Q.finalize();
579 m_R.finalize();
581 m_isQSorted = false;
582
583 m_nonzeropivots = nonzeroCol;
584
585 if(nonzeroCol<n)
586 {
587 // Permute the triangular factor to put the 'dead' columns to the end
588 QRMatrixType tempR(m_R);
589 m_R = tempR * m_pivotperm;
590
591 // Update the column permutation
593 }
594
595 m_isInitialized = true;
596 m_factorizationIsok = true;
597 m_info = Success;
598}
EIGEN_DEVICE_FUNC const SqrtReturnType sqrt() const
Definition ArrayCwiseUnaryOps.h:152
EIGEN_DEVICE_FUNC ColXpr col(Index i)
This is the const version of col().
Definition BlockMethods.h:838
void setIdentity()
Definition PermutationMatrix.h:142
friend class Eigen::Map
Definition PlainObjectBase.h:121
EIGEN_DEVICE_FUNC Derived & setZero(Index size)
Definition CwiseNullaryOp.h:515
Scalar dot(const MatrixBase< OtherDerived > &other) const
RealScalar norm() const
Definition SparseDot.h:84
Base::InnerIterator InnerIterator
Definition SparseMatrix.h:112
void startVec(Index outer)
Definition SparseMatrix.h:412
const StorageIndex * innerNonZeroPtr() const
Definition SparseMatrix.h:175
Scalar & insertBackByOuterInnerUnordered(Index outer, Index inner)
Definition SparseMatrix.h:402
void finalize()
Definition SparseMatrix.h:422
void makeCompressed()
Definition SparseMatrix.h:464
void uncompress()
Definition SparseMatrix.h:495
const StorageIndex * outerIndexPtr() const
Definition SparseMatrix.h:166
Scalar & insertBackByOuterInner(Index outer, Index inner)
Definition SparseMatrix.h:390
void setZero()
Definition SparseMatrix.h:251
bool m_factorizationIsok
Definition SparseQR.h:276
MatrixType::Scalar Scalar
Definition SparseQR.h:81
MatrixType::RealScalar RealScalar
Definition SparseQR.h:82
PermutationType m_pivotperm
Definition SparseQR.h:284
Index m_nonzeropivots
Definition SparseQR.h:288
RealScalar m_threshold
Definition SparseQR.h:286
Matrix< Scalar, Dynamic, 1 > ScalarVector
Definition SparseQR.h:86
SparseMatrix< Scalar, ColMajor, StorageIndex > QRMatrixType
Definition SparseQR.h:84
Matrix< StorageIndex, Dynamic, 1 > IndexVector
Definition SparseQR.h:85
@ InvalidInput
Definition Constants.h:439
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC half abs(const half &a)
Definition Half.h:445
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T maxi(const T &x, const T &y)
Definition MathFunctions.h:823

References col(), Eigen::internal::coletree(), Eigen::PlainObjectBase< Derived >::data(), eigen_assert, Eigen::InvalidInput, Eigen::numext::maxi(), Eigen::PlainObjectBase< Derived >::setConstant(), Eigen::PlainObjectBase< Derived >::setZero(), sqrt(), and Eigen::Success.

Referenced by Eigen::SparseQR< _MatrixType, _OrderingType >::compute().

+ Here is the call graph for this function:
+ Here is the caller graph for this function:

◆ 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()
257 {
258 eigen_assert(m_isInitialized && "Decomposition is not initialized.");
259 return m_info;
260 }

References eigen_assert, Eigen::SparseQR< _MatrixType, _OrderingType >::m_info, and Eigen::SparseSolverBase< SparseQR< _MatrixType, _OrderingType > >::m_isInitialized.

◆ 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.
189{ return m_lastError; }

References Eigen::SparseQR< _MatrixType, _OrderingType >::m_lastError.

◆ matrixQ()

template<typename _MatrixType , typename _OrderingType >
SparseQRMatrixQReturnType< SparseQR > Eigen::SparseQR< _MatrixType, _OrderingType >::matrixQ ( ) 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;

To get a plain SparseMatrix representation of Q:

A versatible sparse matrix representation.
Definition SparseMatrix.h:98
Sparse left-looking rank-revealing QR factorization.
Definition SparseQR.h:73

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.

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

Referenced by Eigen::SparseQR< _MatrixType, _OrderingType >::_solve_impl().

+ Here is the caller graph for this function:

◆ 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
144{ return m_R; }

References Eigen::SparseQR< _MatrixType, _OrderingType >::m_R.

Referenced by Eigen::SparseQR< _MatrixType, _OrderingType >::_solve_impl().

+ Here is the caller graph for this function:

◆ 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()
151 {
152 eigen_assert(m_isInitialized && "The factorization should be called first, use compute()");
153 return m_nonzeropivots;
154 }

References eigen_assert, Eigen::SparseSolverBase< SparseQR< _MatrixType, _OrderingType > >::m_isInitialized, and Eigen::SparseQR< _MatrixType, _OrderingType >::m_nonzeropivots.

Referenced by Eigen::SparseQR< _MatrixType, _OrderingType >::_solve_impl().

+ Here is the caller graph for this function:

◆ rows()

template<typename _MatrixType , typename _OrderingType >
Index Eigen::SparseQR< _MatrixType, _OrderingType >::rows ( ) const
inline
Returns
the number of rows of the represented matrix.
125{ return m_pmat.rows(); }
Index rows() const
Definition SparseMatrix.h:136

References Eigen::SparseQR< _MatrixType, _OrderingType >::m_pmat, and Eigen::SparseMatrix< _Scalar, _Options, _StorageIndex >::rows().

Referenced by Eigen::SparseQR< _MatrixType, _OrderingType >::_solve_impl(), Eigen::SparseQR< _MatrixType, _OrderingType >::solve(), and Eigen::SparseQR< _MatrixType, _OrderingType >::solve().

+ Here is the call graph for this function:
+ Here is the caller graph for this function:

◆ 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.

224 {
225 m_useDefaultThreshold = false;
226 m_threshold = threshold;
227 }

References Eigen::SparseQR< _MatrixType, _OrderingType >::m_threshold, and Eigen::SparseQR< _MatrixType, _OrderingType >::m_useDefaultThreshold.

◆ 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()
235 {
236 eigen_assert(m_isInitialized && "The factorization should be called first, use compute()");
237 eigen_assert(this->rows() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix");
238 return Solve<SparseQR, Rhs>(*this, B.derived());
239 }

References eigen_assert, Eigen::SparseSolverBase< SparseQR< _MatrixType, _OrderingType > >::m_isInitialized, and Eigen::SparseQR< _MatrixType, _OrderingType >::rows().

+ Here is the call graph for this function:

◆ 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
242 {
243 eigen_assert(m_isInitialized && "The factorization should be called first, use compute()");
244 eigen_assert(this->rows() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix");
245 return Solve<SparseQR, Rhs>(*this, B.derived());
246 }

References Eigen::SparseMatrixBase< Derived >::derived(), eigen_assert, Eigen::SparseSolverBase< SparseQR< _MatrixType, _OrderingType > >::m_isInitialized, Eigen::SparseMatrixBase< Derived >::rows(), and Eigen::SparseQR< _MatrixType, _OrderingType >::rows().

+ Here is the call graph for this function:

Friends And Related Symbol Documentation

◆ SparseQR_QProduct

template<typename _MatrixType , typename _OrderingType >
template<typename , typename >
friend struct SparseQR_QProduct
friend

Member Data Documentation

◆ m_analysisIsok

template<typename _MatrixType , typename _OrderingType >
bool Eigen::SparseQR< _MatrixType, _OrderingType >::m_analysisIsok
protected

◆ m_etree

template<typename _MatrixType , typename _OrderingType >
IndexVector Eigen::SparseQR< _MatrixType, _OrderingType >::m_etree
protected

◆ m_factorizationIsok

template<typename _MatrixType , typename _OrderingType >
bool Eigen::SparseQR< _MatrixType, _OrderingType >::m_factorizationIsok
protected

◆ m_firstRowElt

template<typename _MatrixType , typename _OrderingType >
IndexVector Eigen::SparseQR< _MatrixType, _OrderingType >::m_firstRowElt
protected

◆ m_hcoeffs

template<typename _MatrixType , typename _OrderingType >
ScalarVector Eigen::SparseQR< _MatrixType, _OrderingType >::m_hcoeffs
protected

◆ m_info

template<typename _MatrixType , typename _OrderingType >
ComputationInfo Eigen::SparseQR< _MatrixType, _OrderingType >::m_info
mutableprotected

◆ m_isEtreeOk

template<typename _MatrixType , typename _OrderingType >
bool Eigen::SparseQR< _MatrixType, _OrderingType >::m_isEtreeOk
protected

◆ m_isInitialized

bool Eigen::SparseSolverBase< SparseQR< _MatrixType, _OrderingType > >::m_isInitialized
mutableprotectedinherited

◆ m_isQSorted

template<typename _MatrixType , typename _OrderingType >
bool Eigen::SparseQR< _MatrixType, _OrderingType >::m_isQSorted
protected

◆ m_lastError

template<typename _MatrixType , typename _OrderingType >
std::string Eigen::SparseQR< _MatrixType, _OrderingType >::m_lastError
protected

◆ m_nonzeropivots

template<typename _MatrixType , typename _OrderingType >
Index Eigen::SparseQR< _MatrixType, _OrderingType >::m_nonzeropivots
protected

◆ m_outputPerm_c

template<typename _MatrixType , typename _OrderingType >
PermutationType Eigen::SparseQR< _MatrixType, _OrderingType >::m_outputPerm_c
protected

◆ m_perm_c

template<typename _MatrixType , typename _OrderingType >
PermutationType Eigen::SparseQR< _MatrixType, _OrderingType >::m_perm_c
protected

◆ m_pivotperm

template<typename _MatrixType , typename _OrderingType >
PermutationType Eigen::SparseQR< _MatrixType, _OrderingType >::m_pivotperm
protected

◆ m_pmat

template<typename _MatrixType , typename _OrderingType >
QRMatrixType Eigen::SparseQR< _MatrixType, _OrderingType >::m_pmat
protected

◆ m_Q

template<typename _MatrixType , typename _OrderingType >
QRMatrixType Eigen::SparseQR< _MatrixType, _OrderingType >::m_Q
protected

◆ m_R

template<typename _MatrixType , typename _OrderingType >
QRMatrixType Eigen::SparseQR< _MatrixType, _OrderingType >::m_R
protected

◆ m_threshold

template<typename _MatrixType , typename _OrderingType >
RealScalar Eigen::SparseQR< _MatrixType, _OrderingType >::m_threshold
protected

◆ m_useDefaultThreshold

template<typename _MatrixType , typename _OrderingType >
bool Eigen::SparseQR< _MatrixType, _OrderingType >::m_useDefaultThreshold
protected

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