11 #ifndef EIGEN_SPARSE_QR_H
12 #define EIGEN_SPARSE_QR_H
16 template<
typename MatrixType,
typename OrderingType>
class SparseQR;
17 template<
typename SparseQRType>
struct SparseQRMatrixQReturnType;
18 template<
typename SparseQRType>
struct SparseQRMatrixQTransposeReturnType;
19 template<
typename SparseQRType,
typename Derived>
struct SparseQR_QProduct;
21 template <
typename SparseQRType>
struct traits<SparseQRMatrixQReturnType<SparseQRType> >
23 typedef typename SparseQRType::MatrixType ReturnType;
24 typedef typename ReturnType::StorageIndex StorageIndex;
25 typedef typename ReturnType::StorageKind StorageKind;
27 RowsAtCompileTime = Dynamic,
28 ColsAtCompileTime = Dynamic
31 template <
typename SparseQRType>
struct traits<SparseQRMatrixQTransposeReturnType<SparseQRType> >
33 typedef typename SparseQRType::MatrixType ReturnType;
35 template <
typename SparseQRType,
typename Derived>
struct traits<SparseQR_QProduct<SparseQRType, Derived> >
37 typedef typename Derived::PlainObject ReturnType;
70 template<
typename _MatrixType,
typename _OrderingType>
71 class SparseQR :
public SparseSolverBase<SparseQR<_MatrixType,_OrderingType> >
74 typedef SparseSolverBase<SparseQR<_MatrixType,_OrderingType> > Base;
75 using Base::m_isInitialized;
77 using Base::_solve_impl;
78 typedef _MatrixType MatrixType;
79 typedef _OrderingType OrderingType;
80 typedef typename MatrixType::Scalar Scalar;
81 typedef typename MatrixType::RealScalar RealScalar;
82 typedef typename MatrixType::StorageIndex StorageIndex;
83 typedef SparseMatrix<Scalar,ColMajor,StorageIndex> QRMatrixType;
84 typedef Matrix<StorageIndex, Dynamic, 1> IndexVector;
85 typedef Matrix<Scalar, Dynamic, 1> ScalarVector;
86 typedef PermutationMatrix<Dynamic, Dynamic, StorageIndex> PermutationType;
89 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
90 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
94 SparseQR () : m_analysisIsok(false), m_lastError(
""), m_useDefaultThreshold(true),m_isQSorted(false),m_isEtreeOk(false)
103 explicit SparseQR(
const MatrixType& mat) : m_analysisIsok(false), m_lastError(
""), m_useDefaultThreshold(true),m_isQSorted(false),m_isEtreeOk(false)
124 inline Index
rows()
const {
return m_pmat.
rows(); }
132 const QRMatrixType&
matrixR()
const {
return m_R; }
140 eigen_assert(m_isInitialized &&
"The factorization should be called first, use compute()");
141 return m_nonzeropivots;
162 SparseQRMatrixQReturnType<SparseQR>
matrixQ()
const
163 {
return SparseQRMatrixQReturnType<SparseQR>(*this); }
170 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
171 return m_outputPerm_c;
180 template<
typename Rhs,
typename Dest>
183 eigen_assert(m_isInitialized &&
"The factorization should be called first, use compute()");
184 eigen_assert(this->
rows() == B.rows() &&
"SparseQR::solve() : invalid number of rows in the right hand side matrix");
189 typename Dest::PlainObject y, b;
190 y = this->
matrixQ().transpose() * B;
194 y.
resize((std::max<Index>)(
cols(),y.rows()),y.cols());
196 y.bottomRows(y.rows()-
rank).setZero();
213 m_useDefaultThreshold =
false;
214 m_threshold = threshold;
221 template<
typename Rhs>
224 eigen_assert(m_isInitialized &&
"The factorization should be called first, use compute()");
225 eigen_assert(this->
rows() == B.rows() &&
"SparseQR::solve() : invalid number of rows in the right hand side matrix");
228 template<
typename Rhs>
231 eigen_assert(m_isInitialized &&
"The factorization should be called first, use compute()");
232 eigen_assert(this->
rows() == B.
rows() &&
"SparseQR::solve() : invalid number of rows in the right hand side matrix");
246 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
252 inline void _sort_matrix_Q()
254 if(this->m_isQSorted)
return;
258 this->m_isQSorted =
true;
264 bool m_factorizationIsok;
266 std::string m_lastError;
270 ScalarVector m_hcoeffs;
271 PermutationType m_perm_c;
272 PermutationType m_pivotperm;
273 PermutationType m_outputPerm_c;
274 RealScalar m_threshold;
275 bool m_useDefaultThreshold;
276 Index m_nonzeropivots;
278 IndexVector m_firstRowElt;
282 template <
typename,
typename >
friend struct SparseQR_QProduct;
295 template <
typename MatrixType,
typename OrderingType>
298 eigen_assert(mat.isCompressed() &&
"SparseQR requires a sparse matrix in compressed mode. Call .makeCompressed() before passing it to SparseQR");
300 typename internal::conditional<MatrixType::IsRowMajor,QRMatrixType,const MatrixType&>::type matCpy(mat);
303 ord(matCpy, m_perm_c);
304 Index n = mat.cols();
305 Index m = mat.rows();
306 Index diagSize = (std::min)(m,n);
308 if (!m_perm_c.size())
311 m_perm_c.indices().setLinSpaced(n, 0,StorageIndex(n-1));
315 m_outputPerm_c = m_perm_c.inverse();
316 internal::coletree(matCpy, m_etree, m_firstRowElt, m_outputPerm_c.indices().data());
320 m_Q.resize(m, diagSize);
323 m_R.reserve(2*mat.nonZeros());
324 m_Q.reserve(2*mat.nonZeros());
325 m_hcoeffs.resize(diagSize);
326 m_analysisIsok =
true;
336 template <
typename MatrixType,
typename OrderingType>
341 eigen_assert(m_analysisIsok &&
"analyzePattern() should be called before this step");
342 StorageIndex m = StorageIndex(mat.rows());
343 StorageIndex n = StorageIndex(mat.cols());
344 StorageIndex diagSize = (std::min)(m,n);
345 IndexVector mark((std::max)(m,n)); mark.
setConstant(-1);
346 IndexVector Ridx(n), Qidx(m);
347 Index nzcolR, nzcolQ;
348 ScalarVector tval(m);
349 RealScalar pivotThreshold = m_threshold;
356 m_outputPerm_c = m_perm_c.inverse();
357 internal::coletree(m_pmat, m_etree, m_firstRowElt, m_outputPerm_c.indices().data());
368 IndexVector originalOuterIndicesCpy;
369 const StorageIndex *originalOuterIndices = mat.outerIndexPtr();
370 if(MatrixType::IsRowMajor)
372 originalOuterIndicesCpy = IndexVector::Map(m_pmat.outerIndexPtr(),n+1);
373 originalOuterIndices = originalOuterIndicesCpy.
data();
376 for (
int i = 0; i < n; i++)
378 Index p = m_perm_c.size() ? m_perm_c.indices()(i) : i;
379 m_pmat.outerIndexPtr()[p] = originalOuterIndices[i];
380 m_pmat.innerNonZeroPtr()[p] = originalOuterIndices[i+1] - originalOuterIndices[i];
388 if(m_useDefaultThreshold)
390 RealScalar max2Norm = 0.0;
391 for (
int j = 0; j < n; j++) max2Norm = numext::maxi(max2Norm, m_pmat.col(j).norm());
392 if(max2Norm==RealScalar(0))
393 max2Norm = RealScalar(1);
398 m_pivotperm.setIdentity(n);
400 StorageIndex nonzeroCol = 0;
404 for (StorageIndex col = 0; col < n; ++col)
408 mark(nonzeroCol) = col;
409 Qidx(0) = nonzeroCol;
410 nzcolR = 0; nzcolQ = 1;
411 bool found_diag = nonzeroCol>=m;
418 for (
typename QRMatrixType::InnerIterator itp(m_pmat, col); itp || !found_diag; ++itp)
420 StorageIndex curIdx = nonzeroCol;
421 if(itp) curIdx = StorageIndex(itp.row());
422 if(curIdx == nonzeroCol) found_diag =
true;
425 StorageIndex st = m_firstRowElt(curIdx);
428 m_lastError =
"Empty row found during numerical factorization";
435 for (; mark(st) != col; st = m_etree(st))
443 Index nt = nzcolR-bi;
444 for(Index i = 0; i < nt/2; i++) std::swap(Ridx(bi+i), Ridx(nzcolR-i-1));
447 if(itp) tval(curIdx) = itp.value();
448 else tval(curIdx) = Scalar(0);
451 if(curIdx > nonzeroCol && mark(curIdx) != col )
453 Qidx(nzcolQ) = curIdx;
460 for (Index i = nzcolR-1; i >= 0; i--)
462 Index curIdx = Ridx(i);
468 tdot = m_Q.col(curIdx).dot(tval);
470 tdot *= m_hcoeffs(curIdx);
474 for (
typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq)
475 tval(itq.row()) -= itq.value() * tdot;
478 if(m_etree(Ridx(i)) == nonzeroCol)
480 for (
typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq)
482 StorageIndex iQ = StorageIndex(itq.row());
492 Scalar tau = RealScalar(0);
495 if(nonzeroCol < diagSize)
499 Scalar c0 = nzcolQ ? tval(Qidx(0)) : Scalar(0);
502 RealScalar sqrNorm = 0.;
503 for (Index itq = 1; itq < nzcolQ; ++itq) sqrNorm += numext::abs2(tval(Qidx(itq)));
504 if(sqrNorm == RealScalar(0) && numext::imag(c0) == RealScalar(0))
506 beta = numext::real(c0);
512 beta = sqrt(numext::abs2(c0) + sqrNorm);
513 if(numext::real(c0) >= RealScalar(0))
516 for (Index itq = 1; itq < nzcolQ; ++itq)
517 tval(Qidx(itq)) /= (c0 - beta);
518 tau = numext::conj((beta-c0) / beta);
524 for (Index i = nzcolR-1; i >= 0; i--)
526 Index curIdx = Ridx(i);
527 if(curIdx < nonzeroCol)
529 m_R.insertBackByOuterInnerUnordered(col, curIdx) = tval(curIdx);
530 tval(curIdx) = Scalar(0.);
534 if(nonzeroCol < diagSize && abs(beta) >= pivotThreshold)
536 m_R.insertBackByOuterInner(col, nonzeroCol) = beta;
538 m_hcoeffs(nonzeroCol) = tau;
540 for (Index itq = 0; itq < nzcolQ; ++itq)
542 Index iQ = Qidx(itq);
543 m_Q.insertBackByOuterInnerUnordered(nonzeroCol,iQ) = tval(iQ);
544 tval(iQ) = Scalar(0.);
547 if(nonzeroCol<diagSize)
548 m_Q.startVec(nonzeroCol);
553 for (Index j = nonzeroCol; j < n-1; j++)
554 std::swap(m_pivotperm.indices()(j), m_pivotperm.indices()[j+1]);
557 internal::coletree(m_pmat, m_etree, m_firstRowElt, m_pivotperm.indices().data());
562 m_hcoeffs.tail(diagSize-nonzeroCol).setZero();
566 m_Q.makeCompressed();
568 m_R.makeCompressed();
571 m_nonzeropivots = nonzeroCol;
576 QRMatrixType tempR(m_R);
577 m_R = tempR * m_pivotperm;
580 m_outputPerm_c = m_outputPerm_c * m_pivotperm;
583 m_isInitialized =
true;
584 m_factorizationIsok =
true;
588 template <
typename SparseQRType,
typename Derived>
589 struct SparseQR_QProduct : ReturnByValue<SparseQR_QProduct<SparseQRType, Derived> >
591 typedef typename SparseQRType::QRMatrixType MatrixType;
592 typedef typename SparseQRType::Scalar Scalar;
594 SparseQR_QProduct(
const SparseQRType& qr,
const Derived& other,
bool transpose) :
595 m_qr(qr),m_other(other),m_transpose(transpose) {}
596 inline Index rows()
const {
return m_transpose ? m_qr.rows() : m_qr.cols(); }
597 inline Index cols()
const {
return m_other.cols(); }
600 template<
typename DesType>
601 void evalTo(DesType& res)
const
603 Index m = m_qr.rows();
604 Index n = m_qr.cols();
605 Index diagSize = (std::min)(m,n);
609 eigen_assert(m_qr.m_Q.rows() == m_other.rows() &&
"Non conforming object sizes");
611 for(Index j = 0; j < res.cols(); j++){
612 for (Index k = 0; k < diagSize; k++)
614 Scalar tau = Scalar(0);
615 tau = m_qr.m_Q.col(k).dot(res.col(j));
616 if(tau==Scalar(0))
continue;
617 tau = tau * m_qr.m_hcoeffs(k);
618 res.col(j) -= tau * m_qr.m_Q.col(k);
624 eigen_assert(m_qr.m_Q.rows() == m_other.rows() &&
"Non conforming object sizes");
626 for(Index j = 0; j < res.cols(); j++)
628 for (Index k = diagSize-1; k >=0; k--)
630 Scalar tau = Scalar(0);
631 tau = m_qr.m_Q.col(k).dot(res.col(j));
632 if(tau==Scalar(0))
continue;
633 tau = tau * m_qr.m_hcoeffs(k);
634 res.col(j) -= tau * m_qr.m_Q.col(k);
640 const SparseQRType& m_qr;
641 const Derived& m_other;
645 template<
typename SparseQRType>
646 struct SparseQRMatrixQReturnType :
public EigenBase<SparseQRMatrixQReturnType<SparseQRType> >
648 typedef typename SparseQRType::Scalar Scalar;
649 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
651 RowsAtCompileTime = Dynamic,
652 ColsAtCompileTime = Dynamic
654 explicit SparseQRMatrixQReturnType(
const SparseQRType& qr) : m_qr(qr) {}
655 template<
typename Derived>
656 SparseQR_QProduct<SparseQRType, Derived> operator*(
const MatrixBase<Derived>& other)
658 return SparseQR_QProduct<SparseQRType,Derived>(m_qr,other.derived(),
false);
660 SparseQRMatrixQTransposeReturnType<SparseQRType> adjoint()
const
662 return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr);
664 inline Index rows()
const {
return m_qr.rows(); }
665 inline Index cols()
const {
return (std::min)(m_qr.rows(),m_qr.cols()); }
667 SparseQRMatrixQTransposeReturnType<SparseQRType> transpose()
const
669 return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr);
671 const SparseQRType& m_qr;
674 template<
typename SparseQRType>
675 struct SparseQRMatrixQTransposeReturnType
677 explicit SparseQRMatrixQTransposeReturnType(
const SparseQRType& qr) : m_qr(qr) {}
678 template<
typename Derived>
679 SparseQR_QProduct<SparseQRType,Derived> operator*(
const MatrixBase<Derived>& other)
681 return SparseQR_QProduct<SparseQRType,Derived>(m_qr,other.derived(),
true);
683 const SparseQRType& m_qr;
688 template<
typename SparseQRType>
689 struct evaluator_traits<SparseQRMatrixQReturnType<SparseQRType> >
691 typedef typename SparseQRType::MatrixType MatrixType;
692 typedef typename storage_kind_to_evaluator_kind<typename MatrixType::StorageKind>::Kind Kind;
693 typedef SparseShape Shape;
694 static const int AssumeAliasing = 0;
697 template<
typename DstXprType,
typename SparseQRType>
698 struct Assignment<DstXprType, SparseQRMatrixQReturnType<SparseQRType>,
internal::assign_op<typename DstXprType::Scalar>, Sparse2Sparse>
700 typedef SparseQRMatrixQReturnType<SparseQRType> SrcXprType;
701 typedef typename DstXprType::Scalar Scalar;
702 typedef typename DstXprType::StorageIndex StorageIndex;
703 static void run(DstXprType &dst,
const SrcXprType &src,
const internal::assign_op<Scalar> &)
705 typename DstXprType::PlainObject idMat(src.m_qr.rows(), src.m_qr.rows());
708 const_cast<SparseQRType *
>(&src.m_qr)->_sort_matrix_Q();
709 dst = SparseQR_QProduct<SparseQRType, DstXprType>(src.m_qr, idMat,
false);
713 template<
typename DstXprType,
typename SparseQRType>
714 struct Assignment<DstXprType, SparseQRMatrixQReturnType<SparseQRType>,
internal::assign_op<typename DstXprType::Scalar>, Sparse2Dense>
716 typedef SparseQRMatrixQReturnType<SparseQRType> SrcXprType;
717 typedef typename DstXprType::Scalar Scalar;
718 typedef typename DstXprType::StorageIndex StorageIndex;
719 static void run(DstXprType &dst,
const SrcXprType &src,
const internal::assign_op<Scalar> &)
721 dst = src.m_qr.matrixQ() * DstXprType::Identity(src.m_qr.rows(), src.m_qr.rows());
Index cols() const
Definition: SparseMatrix.h:133
A versatible sparse matrix representation.
Definition: SparseMatrix.h:92
RowsBlockXpr topRows(Index n)
Definition: DenseBase.h:400
Index rank() const
Definition: SparseQR.h:138
Index size() const
Definition: PermutationMatrix.h:109
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:107
Derived & derived()
Definition: EigenBase.h:44
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:37
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SparseQR.h:244
const QRMatrixType & matrixR() const
Definition: SparseQR.h:132
void resize(Index newSize)
Definition: DenseBase.h:245
void factorize(const MatrixType &mat)
Performs the numerical QR factorization of the input matrix.
Definition: SparseQR.h:337
Derived & setZero(Index size)
Definition: CwiseNullaryOp.h:520
Sparse left-looking rank-revealing QR factorization.
Definition: SparseQR.h:16
Base class of any sparse matrices or sparse expressions.
Definition: ForwardDeclarations.h:278
Index cols() const
Definition: SparseQR.h:128
Definition: Constants.h:439
Definition: Constants.h:432
void analyzePattern(const MatrixType &mat)
Preprocessing step of a QR factorization.
Definition: SparseQR.h:296
Block< Derived > topLeftCorner(Index cRows, Index cCols)
Definition: SparseMatrixBase.h:164
Definition: Eigen_Colamd.h:54
Index rows() const
Definition: SparseQR.h:124
const PermutationType & colsPermutation() const
Definition: SparseQR.h:168
const Solve< SparseQR, Rhs > solve(const MatrixBase< Rhs > &B) const
Definition: SparseQR.h:222
SparseQRMatrixQReturnType< SparseQR > matrixQ() const
Definition: SparseQR.h:162
Index rows() const
Definition: SparseMatrixBase.h:163
Pseudo expression representing a solving operation.
Definition: Solve.h:62
const Scalar * data() const
Definition: PlainObjectBase.h:228
ComputationInfo
Definition: Constants.h:430
SparseQR(const MatrixType &mat)
Definition: SparseQR.h:103
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
Index rows() const
Definition: SparseMatrix.h:131
void setPivotThreshold(const RealScalar &threshold)
Definition: SparseQR.h:211
std::string lastErrorMessage() const
Definition: SparseQR.h:177
Derived & setConstant(Index size, const Scalar &value)
Definition: CwiseNullaryOp.h:352
void compute(const MatrixType &mat)
Definition: SparseQR.h:114