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;
68 template<
typename _MatrixType,
typename _OrderingType>
69 class SparseQR :
public SparseSolverBase<SparseQR<_MatrixType,_OrderingType> >
72 typedef SparseSolverBase<SparseQR<_MatrixType,_OrderingType> > Base;
73 using Base::m_isInitialized;
75 using Base::_solve_impl;
76 typedef _MatrixType MatrixType;
77 typedef _OrderingType OrderingType;
78 typedef typename MatrixType::Scalar Scalar;
79 typedef typename MatrixType::RealScalar RealScalar;
80 typedef typename MatrixType::StorageIndex StorageIndex;
81 typedef SparseMatrix<Scalar,ColMajor,StorageIndex> QRMatrixType;
82 typedef Matrix<StorageIndex, Dynamic, 1> IndexVector;
83 typedef Matrix<Scalar, Dynamic, 1> ScalarVector;
84 typedef PermutationMatrix<Dynamic, Dynamic, StorageIndex> PermutationType;
86 SparseQR () : m_analysisIsok(false), m_lastError(
""), m_useDefaultThreshold(true),m_isQSorted(false),m_isEtreeOk(false)
95 explicit SparseQR(
const MatrixType& mat) : m_analysisIsok(false), m_lastError(
""), m_useDefaultThreshold(true),m_isQSorted(false),m_isEtreeOk(false)
116 inline Index
rows()
const {
return m_pmat.
rows(); }
124 const QRMatrixType&
matrixR()
const {
return m_R; }
132 eigen_assert(m_isInitialized &&
"The factorization should be called first, use compute()");
133 return m_nonzeropivots;
154 SparseQRMatrixQReturnType<SparseQR>
matrixQ()
const
155 {
return SparseQRMatrixQReturnType<SparseQR>(*this); }
162 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
163 return m_outputPerm_c;
172 template<
typename Rhs,
typename Dest>
175 eigen_assert(m_isInitialized &&
"The factorization should be called first, use compute()");
176 eigen_assert(this->
rows() == B.rows() &&
"SparseQR::solve() : invalid number of rows in the right hand side matrix");
181 typename Dest::PlainObject y, b;
182 y = this->
matrixQ().transpose() * B;
186 y.
resize((std::max<Index>)(
cols(),y.rows()),y.cols());
188 y.bottomRows(y.rows()-
rank).setZero();
205 m_useDefaultThreshold =
false;
206 m_threshold = threshold;
213 template<
typename Rhs>
216 eigen_assert(m_isInitialized &&
"The factorization should be called first, use compute()");
217 eigen_assert(this->
rows() == B.rows() &&
"SparseQR::solve() : invalid number of rows in the right hand side matrix");
220 template<
typename Rhs>
223 eigen_assert(m_isInitialized &&
"The factorization should be called first, use compute()");
224 eigen_assert(this->
rows() == B.
rows() &&
"SparseQR::solve() : invalid number of rows in the right hand side matrix");
238 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
244 inline void _sort_matrix_Q()
246 if(this->m_isQSorted)
return;
250 this->m_isQSorted =
true;
256 bool m_factorizationIsok;
258 std::string m_lastError;
262 ScalarVector m_hcoeffs;
263 PermutationType m_perm_c;
264 PermutationType m_pivotperm;
265 PermutationType m_outputPerm_c;
266 RealScalar m_threshold;
267 bool m_useDefaultThreshold;
268 Index m_nonzeropivots;
270 IndexVector m_firstRowElt;
274 template <
typename,
typename >
friend struct SparseQR_QProduct;
287 template <
typename MatrixType,
typename OrderingType>
290 eigen_assert(mat.isCompressed() &&
"SparseQR requires a sparse matrix in compressed mode. Call .makeCompressed() before passing it to SparseQR");
292 typename internal::conditional<MatrixType::IsRowMajor,QRMatrixType,const MatrixType&>::type matCpy(mat);
295 ord(matCpy, m_perm_c);
296 Index n = mat.cols();
297 Index m = mat.rows();
298 Index diagSize = (std::min)(m,n);
300 if (!m_perm_c.size())
303 m_perm_c.indices().setLinSpaced(n, 0,StorageIndex(n-1));
307 m_outputPerm_c = m_perm_c.inverse();
308 internal::coletree(matCpy, m_etree, m_firstRowElt, m_outputPerm_c.indices().data());
312 m_Q.resize(m, diagSize);
315 m_R.reserve(2*mat.nonZeros());
316 m_Q.reserve(2*mat.nonZeros());
317 m_hcoeffs.resize(diagSize);
318 m_analysisIsok =
true;
328 template <
typename MatrixType,
typename OrderingType>
333 eigen_assert(m_analysisIsok &&
"analyzePattern() should be called before this step");
334 StorageIndex m = StorageIndex(mat.rows());
335 StorageIndex n = StorageIndex(mat.cols());
336 StorageIndex diagSize = (std::min)(m,n);
337 IndexVector mark((std::max)(m,n)); mark.
setConstant(-1);
338 IndexVector Ridx(n), Qidx(m);
339 Index nzcolR, nzcolQ;
340 ScalarVector tval(m);
341 RealScalar pivotThreshold = m_threshold;
348 m_outputPerm_c = m_perm_c.inverse();
349 internal::coletree(m_pmat, m_etree, m_firstRowElt, m_outputPerm_c.indices().data());
360 IndexVector originalOuterIndicesCpy;
361 const StorageIndex *originalOuterIndices = mat.outerIndexPtr();
362 if(MatrixType::IsRowMajor)
364 originalOuterIndicesCpy = IndexVector::Map(m_pmat.outerIndexPtr(),n+1);
365 originalOuterIndices = originalOuterIndicesCpy.
data();
368 for (
int i = 0; i < n; i++)
370 Index p = m_perm_c.size() ? m_perm_c.indices()(i) : i;
371 m_pmat.outerIndexPtr()[p] = originalOuterIndices[i];
372 m_pmat.innerNonZeroPtr()[p] = originalOuterIndices[i+1] - originalOuterIndices[i];
380 if(m_useDefaultThreshold)
382 RealScalar max2Norm = 0.0;
383 for (
int j = 0; j < n; j++) max2Norm = numext::maxi(max2Norm, m_pmat.col(j).norm());
384 if(max2Norm==RealScalar(0))
385 max2Norm = RealScalar(1);
390 m_pivotperm.setIdentity(n);
392 StorageIndex nonzeroCol = 0;
396 for (StorageIndex col = 0; col < n; ++col)
400 mark(nonzeroCol) = col;
401 Qidx(0) = nonzeroCol;
402 nzcolR = 0; nzcolQ = 1;
403 bool found_diag = nonzeroCol>=m;
410 for (
typename QRMatrixType::InnerIterator itp(m_pmat, col); itp || !found_diag; ++itp)
412 StorageIndex curIdx = nonzeroCol;
413 if(itp) curIdx = StorageIndex(itp.row());
414 if(curIdx == nonzeroCol) found_diag =
true;
417 StorageIndex st = m_firstRowElt(curIdx);
420 m_lastError =
"Empty row found during numerical factorization";
427 for (; mark(st) != col; st = m_etree(st))
435 Index nt = nzcolR-bi;
436 for(Index i = 0; i < nt/2; i++) std::swap(Ridx(bi+i), Ridx(nzcolR-i-1));
439 if(itp) tval(curIdx) = itp.value();
440 else tval(curIdx) = Scalar(0);
443 if(curIdx > nonzeroCol && mark(curIdx) != col )
445 Qidx(nzcolQ) = curIdx;
452 for (Index i = nzcolR-1; i >= 0; i--)
454 Index curIdx = Ridx(i);
460 tdot = m_Q.col(curIdx).dot(tval);
462 tdot *= m_hcoeffs(curIdx);
466 for (
typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq)
467 tval(itq.row()) -= itq.value() * tdot;
470 if(m_etree(Ridx(i)) == nonzeroCol)
472 for (
typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq)
474 StorageIndex iQ = StorageIndex(itq.row());
484 Scalar tau = RealScalar(0);
487 if(nonzeroCol < diagSize)
491 Scalar c0 = nzcolQ ? tval(Qidx(0)) : Scalar(0);
494 RealScalar sqrNorm = 0.;
495 for (Index itq = 1; itq < nzcolQ; ++itq) sqrNorm += numext::abs2(tval(Qidx(itq)));
496 if(sqrNorm == RealScalar(0) && numext::imag(c0) == RealScalar(0))
498 beta = numext::real(c0);
504 beta = sqrt(numext::abs2(c0) + sqrNorm);
505 if(numext::real(c0) >= RealScalar(0))
508 for (Index itq = 1; itq < nzcolQ; ++itq)
509 tval(Qidx(itq)) /= (c0 - beta);
510 tau = numext::conj((beta-c0) / beta);
516 for (Index i = nzcolR-1; i >= 0; i--)
518 Index curIdx = Ridx(i);
519 if(curIdx < nonzeroCol)
521 m_R.insertBackByOuterInnerUnordered(col, curIdx) = tval(curIdx);
522 tval(curIdx) = Scalar(0.);
526 if(nonzeroCol < diagSize && abs(beta) >= pivotThreshold)
528 m_R.insertBackByOuterInner(col, nonzeroCol) = beta;
530 m_hcoeffs(nonzeroCol) = tau;
532 for (Index itq = 0; itq < nzcolQ; ++itq)
534 Index iQ = Qidx(itq);
535 m_Q.insertBackByOuterInnerUnordered(nonzeroCol,iQ) = tval(iQ);
536 tval(iQ) = Scalar(0.);
539 if(nonzeroCol<diagSize)
540 m_Q.startVec(nonzeroCol);
545 for (Index j = nonzeroCol; j < n-1; j++)
546 std::swap(m_pivotperm.indices()(j), m_pivotperm.indices()[j+1]);
549 internal::coletree(m_pmat, m_etree, m_firstRowElt, m_pivotperm.indices().data());
554 m_hcoeffs.tail(diagSize-nonzeroCol).setZero();
558 m_Q.makeCompressed();
560 m_R.makeCompressed();
563 m_nonzeropivots = nonzeroCol;
568 QRMatrixType tempR(m_R);
569 m_R = tempR * m_pivotperm;
572 m_outputPerm_c = m_outputPerm_c * m_pivotperm;
575 m_isInitialized =
true;
576 m_factorizationIsok =
true;
580 template <
typename SparseQRType,
typename Derived>
581 struct SparseQR_QProduct : ReturnByValue<SparseQR_QProduct<SparseQRType, Derived> >
583 typedef typename SparseQRType::QRMatrixType MatrixType;
584 typedef typename SparseQRType::Scalar Scalar;
586 SparseQR_QProduct(
const SparseQRType& qr,
const Derived& other,
bool transpose) :
587 m_qr(qr),m_other(other),m_transpose(transpose) {}
588 inline Index rows()
const {
return m_transpose ? m_qr.rows() : m_qr.cols(); }
589 inline Index cols()
const {
return m_other.cols(); }
592 template<
typename DesType>
593 void evalTo(DesType& res)
const
595 Index m = m_qr.rows();
596 Index n = m_qr.cols();
597 Index diagSize = (std::min)(m,n);
601 eigen_assert(m_qr.m_Q.rows() == m_other.rows() &&
"Non conforming object sizes");
603 for(Index j = 0; j < res.cols(); j++){
604 for (Index k = 0; k < diagSize; k++)
606 Scalar tau = Scalar(0);
607 tau = m_qr.m_Q.col(k).dot(res.col(j));
608 if(tau==Scalar(0))
continue;
609 tau = tau * m_qr.m_hcoeffs(k);
610 res.col(j) -= tau * m_qr.m_Q.col(k);
616 eigen_assert(m_qr.m_Q.rows() == m_other.rows() &&
"Non conforming object sizes");
618 for(Index j = 0; j < res.cols(); j++)
620 for (Index k = diagSize-1; k >=0; k--)
622 Scalar tau = Scalar(0);
623 tau = m_qr.m_Q.col(k).dot(res.col(j));
624 if(tau==Scalar(0))
continue;
625 tau = tau * m_qr.m_hcoeffs(k);
626 res.col(j) -= tau * m_qr.m_Q.col(k);
632 const SparseQRType& m_qr;
633 const Derived& m_other;
637 template<
typename SparseQRType>
638 struct SparseQRMatrixQReturnType :
public EigenBase<SparseQRMatrixQReturnType<SparseQRType> >
640 typedef typename SparseQRType::Scalar Scalar;
641 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
643 RowsAtCompileTime = Dynamic,
644 ColsAtCompileTime = Dynamic
646 explicit SparseQRMatrixQReturnType(
const SparseQRType& qr) : m_qr(qr) {}
647 template<
typename Derived>
648 SparseQR_QProduct<SparseQRType, Derived> operator*(
const MatrixBase<Derived>& other)
650 return SparseQR_QProduct<SparseQRType,Derived>(m_qr,other.derived(),
false);
652 SparseQRMatrixQTransposeReturnType<SparseQRType> adjoint()
const
654 return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr);
656 inline Index rows()
const {
return m_qr.rows(); }
657 inline Index cols()
const {
return (std::min)(m_qr.rows(),m_qr.cols()); }
659 SparseQRMatrixQTransposeReturnType<SparseQRType> transpose()
const
661 return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr);
663 const SparseQRType& m_qr;
666 template<
typename SparseQRType>
667 struct SparseQRMatrixQTransposeReturnType
669 explicit SparseQRMatrixQTransposeReturnType(
const SparseQRType& qr) : m_qr(qr) {}
670 template<
typename Derived>
671 SparseQR_QProduct<SparseQRType,Derived> operator*(
const MatrixBase<Derived>& other)
673 return SparseQR_QProduct<SparseQRType,Derived>(m_qr,other.derived(),
true);
675 const SparseQRType& m_qr;
680 template<
typename SparseQRType>
681 struct evaluator_traits<SparseQRMatrixQReturnType<SparseQRType> >
683 typedef typename SparseQRType::MatrixType MatrixType;
684 typedef typename storage_kind_to_evaluator_kind<typename MatrixType::StorageKind>::Kind Kind;
685 typedef SparseShape Shape;
686 static const int AssumeAliasing = 0;
689 template<
typename DstXprType,
typename SparseQRType>
690 struct Assignment<DstXprType, SparseQRMatrixQReturnType<SparseQRType>,
internal::assign_op<typename DstXprType::Scalar>, Sparse2Sparse>
692 typedef SparseQRMatrixQReturnType<SparseQRType> SrcXprType;
693 typedef typename DstXprType::Scalar Scalar;
694 typedef typename DstXprType::StorageIndex StorageIndex;
695 static void run(DstXprType &dst,
const SrcXprType &src,
const internal::assign_op<Scalar> &)
697 typename DstXprType::PlainObject idMat(src.m_qr.rows(), src.m_qr.rows());
700 const_cast<SparseQRType *
>(&src.m_qr)->_sort_matrix_Q();
701 dst = SparseQR_QProduct<SparseQRType, DstXprType>(src.m_qr, idMat,
false);
705 template<
typename DstXprType,
typename SparseQRType>
706 struct Assignment<DstXprType, SparseQRMatrixQReturnType<SparseQRType>,
internal::assign_op<typename DstXprType::Scalar>, Sparse2Dense>
708 typedef SparseQRMatrixQReturnType<SparseQRType> SrcXprType;
709 typedef typename DstXprType::Scalar Scalar;
710 typedef typename DstXprType::StorageIndex StorageIndex;
711 static void run(DstXprType &dst,
const SrcXprType &src,
const internal::assign_op<Scalar> &)
713 dst = src.m_qr.matrixQ() * DstXprType::Identity(src.m_qr.rows(), src.m_qr.rows());
Index cols() const
Definition: SparseMatrix.h:132
A versatible sparse matrix representation.
Definition: SparseMatrix.h:92
RowsBlockXpr topRows(Index n)
Definition: DenseBase.h:400
Index rank() const
Definition: SparseQR.h:130
Index size() const
Definition: PermutationMatrix.h:110
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:236
const QRMatrixType & matrixR() const
Definition: SparseQR.h:124
void resize(Index newSize)
Definition: DenseBase.h:252
void factorize(const MatrixType &mat)
Performs the numerical QR factorization of the input matrix.
Definition: SparseQR.h:329
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: SparseMatrixBase.h:26
Index cols() const
Definition: SparseQR.h:120
Definition: Constants.h:431
Definition: Constants.h:424
void analyzePattern(const MatrixType &mat)
Preprocessing step of a QR factorization.
Definition: SparseQR.h:288
Block< Derived > topLeftCorner(Index cRows, Index cCols)
Definition: SparseMatrixBase.h:164
Definition: Eigen_Colamd.h:54
Index rows() const
Definition: SparseQR.h:116
const PermutationType & colsPermutation() const
Definition: SparseQR.h:160
const Solve< SparseQR, Rhs > solve(const MatrixBase< Rhs > &B) const
Definition: SparseQR.h:214
SparseQRMatrixQReturnType< SparseQR > matrixQ() const
Definition: SparseQR.h:154
Index rows() const
Definition: SparseMatrixBase.h:152
Pseudo expression representing a solving operation.
Definition: Solve.h:63
const Scalar * data() const
Definition: PlainObjectBase.h:228
ComputationInfo
Definition: Constants.h:422
SparseQR(const MatrixType &mat)
Definition: SparseQR.h:95
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
Index rows() const
Definition: SparseMatrix.h:130
void setPivotThreshold(const RealScalar &threshold)
Definition: SparseQR.h:203
std::string lastErrorMessage() const
Definition: SparseQR.h:169
Derived & setConstant(Index size, const Scalar &value)
Definition: CwiseNullaryOp.h:352
void compute(const MatrixType &mat)
Definition: SparseQR.h:106