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;
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;
71 template<
typename _MatrixType,
typename _OrderingType>
72 class SparseQR :
public SparseSolverBase<SparseQR<_MatrixType,_OrderingType> >
75 typedef SparseSolverBase<SparseQR<_MatrixType,_OrderingType> > Base;
76 using Base::m_isInitialized;
78 using Base::_solve_impl;
79 typedef _MatrixType MatrixType;
80 typedef _OrderingType OrderingType;
81 typedef typename MatrixType::Scalar Scalar;
82 typedef typename MatrixType::RealScalar RealScalar;
83 typedef typename MatrixType::StorageIndex StorageIndex;
84 typedef SparseMatrix<Scalar,ColMajor,StorageIndex> QRMatrixType;
85 typedef Matrix<StorageIndex, Dynamic, 1> IndexVector;
86 typedef Matrix<Scalar, Dynamic, 1> ScalarVector;
87 typedef PermutationMatrix<Dynamic, Dynamic, StorageIndex> PermutationType;
90 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
91 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
95 SparseQR () : m_analysisIsok(false), m_lastError(
""), m_useDefaultThreshold(true),m_isQSorted(false),m_isEtreeOk(false)
104 explicit SparseQR(
const MatrixType& mat) : m_analysisIsok(false), m_lastError(
""), m_useDefaultThreshold(true),m_isQSorted(false),m_isEtreeOk(false)
152 eigen_assert(m_isInitialized &&
"The factorization should be called first, use compute()");
153 return m_nonzeropivots;
174 SparseQRMatrixQReturnType<SparseQR>
matrixQ()
const 175 {
return SparseQRMatrixQReturnType<SparseQR>(*
this); }
182 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
183 return m_outputPerm_c;
192 template<
typename Rhs,
typename Dest>
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");
201 typename Dest::PlainObject y, b;
202 y = this->
matrixQ().adjoint() * B;
206 y.resize((std::max<Index>)(
cols(),y.rows()),y.cols());
208 y.bottomRows(y.rows()-
rank).setZero();
212 else dest = y.topRows(
cols());
225 m_useDefaultThreshold =
false;
226 m_threshold = threshold;
233 template<
typename Rhs>
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");
240 template<
typename Rhs>
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");
258 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
264 inline void _sort_matrix_Q()
266 if(this->m_isQSorted)
return;
270 this->m_isQSorted =
true;
276 bool m_factorizationIsok;
278 std::string m_lastError;
282 ScalarVector m_hcoeffs;
283 PermutationType m_perm_c;
284 PermutationType m_pivotperm;
285 PermutationType m_outputPerm_c;
286 RealScalar m_threshold;
287 bool m_useDefaultThreshold;
288 Index m_nonzeropivots;
290 IndexVector m_firstRowElt;
294 template <
typename,
typename >
friend struct SparseQR_QProduct;
307 template <
typename MatrixType,
typename OrderingType>
310 eigen_assert(mat.isCompressed() &&
"SparseQR requires a sparse matrix in compressed mode. Call .makeCompressed() before passing it to SparseQR");
312 typename internal::conditional<MatrixType::IsRowMajor,QRMatrixType,const MatrixType&>::type matCpy(mat);
315 ord(matCpy, m_perm_c);
316 Index n = mat.cols();
317 Index m = mat.rows();
318 Index diagSize = (std::min)(m,n);
320 if (!m_perm_c.size())
323 m_perm_c.indices().setLinSpaced(n, 0,StorageIndex(n-1));
327 m_outputPerm_c = m_perm_c.inverse();
328 internal::coletree(matCpy, m_etree, m_firstRowElt, m_outputPerm_c.indices().data());
332 m_Q.resize(m, diagSize);
335 m_R.reserve(2*mat.nonZeros());
336 m_Q.reserve(2*mat.nonZeros());
337 m_hcoeffs.resize(diagSize);
338 m_analysisIsok =
true;
348 template <
typename MatrixType,
typename OrderingType>
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);
359 Index nzcolR, nzcolQ;
361 RealScalar pivotThreshold = m_threshold;
368 m_outputPerm_c = m_perm_c.inverse();
369 internal::coletree(m_pmat, m_etree, m_firstRowElt, m_outputPerm_c.indices().data());
381 const StorageIndex *originalOuterIndices = mat.outerIndexPtr();
382 if(MatrixType::IsRowMajor)
384 originalOuterIndicesCpy = IndexVector::Map(m_pmat.outerIndexPtr(),n+1);
385 originalOuterIndices = originalOuterIndicesCpy.
data();
388 for (
int i = 0; i < n; i++)
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];
400 if(m_useDefaultThreshold)
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);
410 m_pivotperm.setIdentity(n);
412 StorageIndex nonzeroCol = 0;
416 for (StorageIndex col = 0; col < n; ++col)
420 mark(nonzeroCol) = col;
421 Qidx(0) = nonzeroCol;
422 nzcolR = 0; nzcolQ = 1;
423 bool found_diag = nonzeroCol>=m;
430 for (
typename QRMatrixType::InnerIterator itp(m_pmat, col); itp || !found_diag; ++itp)
432 StorageIndex curIdx = nonzeroCol;
433 if(itp) curIdx = StorageIndex(itp.row());
434 if(curIdx == nonzeroCol) found_diag =
true;
437 StorageIndex st = m_firstRowElt(curIdx);
440 m_lastError =
"Empty row found during numerical factorization";
447 for (; mark(st) != col; st = m_etree(st))
455 Index nt = nzcolR-bi;
456 for(
Index i = 0; i < nt/2; i++) std::swap(Ridx(bi+i), Ridx(nzcolR-i-1));
459 if(itp) tval(curIdx) = itp.value();
460 else tval(curIdx) = Scalar(0);
463 if(curIdx > nonzeroCol && mark(curIdx) != col )
465 Qidx(nzcolQ) = curIdx;
472 for (
Index i = nzcolR-1; i >= 0; i--)
474 Index curIdx = Ridx(i);
480 tdot = m_Q.col(curIdx).dot(tval);
482 tdot *= m_hcoeffs(curIdx);
486 for (
typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq)
487 tval(itq.row()) -= itq.value() * tdot;
490 if(m_etree(Ridx(i)) == nonzeroCol)
492 for (
typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq)
494 StorageIndex iQ = StorageIndex(itq.row());
504 Scalar tau = RealScalar(0);
507 if(nonzeroCol < diagSize)
511 Scalar c0 = nzcolQ ? tval(Qidx(0)) : Scalar(0);
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))
518 beta = numext::real(c0);
524 beta =
sqrt(numext::abs2(c0) + sqrNorm);
525 if(numext::real(c0) >= RealScalar(0))
528 for (
Index itq = 1; itq < nzcolQ; ++itq)
529 tval(Qidx(itq)) /= (c0 - beta);
530 tau = numext::conj((beta-c0) / beta);
536 for (
Index i = nzcolR-1; i >= 0; i--)
538 Index curIdx = Ridx(i);
539 if(curIdx < nonzeroCol)
541 m_R.insertBackByOuterInnerUnordered(col, curIdx) = tval(curIdx);
542 tval(curIdx) = Scalar(0.);
546 if(nonzeroCol < diagSize &&
abs(beta) >= pivotThreshold)
548 m_R.insertBackByOuterInner(col, nonzeroCol) = beta;
550 m_hcoeffs(nonzeroCol) = tau;
552 for (
Index itq = 0; itq < nzcolQ; ++itq)
554 Index iQ = Qidx(itq);
555 m_Q.insertBackByOuterInnerUnordered(nonzeroCol,iQ) = tval(iQ);
556 tval(iQ) = Scalar(0.);
559 if(nonzeroCol<diagSize)
560 m_Q.startVec(nonzeroCol);
565 for (
Index j = nonzeroCol; j < n-1; j++)
566 std::swap(m_pivotperm.indices()(j), m_pivotperm.indices()[j+1]);
569 internal::coletree(m_pmat, m_etree, m_firstRowElt, m_pivotperm.indices().data());
574 m_hcoeffs.tail(diagSize-nonzeroCol).setZero();
578 m_Q.makeCompressed();
580 m_R.makeCompressed();
583 m_nonzeropivots = nonzeroCol;
589 m_R = tempR * m_pivotperm;
592 m_outputPerm_c = m_outputPerm_c * m_pivotperm;
595 m_isInitialized =
true;
596 m_factorizationIsok =
true;
600 template <
typename SparseQRType,
typename Derived>
601 struct SparseQR_QProduct : ReturnByValue<SparseQR_QProduct<SparseQRType, Derived> >
603 typedef typename SparseQRType::QRMatrixType MatrixType;
604 typedef typename SparseQRType::Scalar Scalar;
606 SparseQR_QProduct(
const SparseQRType& qr,
const Derived& other,
bool transpose) :
607 m_qr(qr),m_other(other),m_transpose(transpose) {}
608 inline Index rows()
const {
return m_qr.matrixQ().rows(); }
609 inline Index cols()
const {
return m_other.cols(); }
612 template<
typename DesType>
613 void evalTo(DesType& res)
const 615 Index m = m_qr.rows();
616 Index n = m_qr.cols();
617 Index diagSize = (std::min)(m,n);
621 eigen_assert(m_qr.m_Q.rows() == m_other.rows() &&
"Non conforming object sizes");
623 for(
Index j = 0; j < res.cols(); j++){
624 for (
Index k = 0; k < diagSize; k++)
626 Scalar tau = Scalar(0);
627 tau = m_qr.m_Q.col(k).dot(res.col(j));
628 if(tau==Scalar(0))
continue;
629 tau = tau * m_qr.m_hcoeffs(k);
630 res.col(j) -= tau * m_qr.m_Q.col(k);
636 eigen_assert(m_qr.matrixQ().cols() == m_other.rows() &&
"Non conforming object sizes");
638 res.conservativeResize(rows(), cols());
641 for(
Index j = 0; j < res.cols(); j++)
643 for (
Index k = diagSize-1; k >=0; k--)
645 Scalar tau = Scalar(0);
646 tau = m_qr.m_Q.col(k).dot(res.col(j));
647 if(tau==Scalar(0))
continue;
648 tau = tau * numext::conj(m_qr.m_hcoeffs(k));
649 res.col(j) -= tau * m_qr.m_Q.col(k);
655 const SparseQRType& m_qr;
656 const Derived& m_other;
660 template<
typename SparseQRType>
661 struct SparseQRMatrixQReturnType :
public EigenBase<SparseQRMatrixQReturnType<SparseQRType> >
663 typedef typename SparseQRType::Scalar Scalar;
664 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
669 explicit SparseQRMatrixQReturnType(
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(),
false);
676 SparseQRMatrixQTransposeReturnType<SparseQRType> adjoint()
const 678 return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr);
680 inline Index rows()
const {
return m_qr.rows(); }
681 inline Index cols()
const {
return m_qr.rows(); }
683 SparseQRMatrixQTransposeReturnType<SparseQRType> transpose()
const 685 return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr);
687 const SparseQRType& m_qr;
691 template<
typename SparseQRType>
692 struct SparseQRMatrixQTransposeReturnType
694 explicit SparseQRMatrixQTransposeReturnType(
const SparseQRType& qr) : m_qr(qr) {}
695 template<
typename Derived>
696 SparseQR_QProduct<SparseQRType,Derived>
operator*(
const MatrixBase<Derived>& other)
698 return SparseQR_QProduct<SparseQRType,Derived>(m_qr,other.derived(),
true);
700 const SparseQRType& m_qr;
705 template<
typename SparseQRType>
706 struct evaluator_traits<SparseQRMatrixQReturnType<SparseQRType> >
708 typedef typename SparseQRType::MatrixType MatrixType;
709 typedef typename storage_kind_to_evaluator_kind<typename MatrixType::StorageKind>::Kind Kind;
710 typedef SparseShape Shape;
713 template<
typename DstXprType,
typename SparseQRType>
714 struct Assignment<DstXprType, SparseQRMatrixQReturnType<SparseQRType>, internal::assign_op<typename DstXprType::Scalar,typename DstXprType::Scalar>, Sparse2Sparse>
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,Scalar> &)
721 typename DstXprType::PlainObject idMat(src.rows(), src.cols());
724 const_cast<SparseQRType *>(&src.m_qr)->_sort_matrix_Q();
725 dst = SparseQR_QProduct<SparseQRType, DstXprType>(src.m_qr, idMat,
false);
729 template<
typename DstXprType,
typename SparseQRType>
730 struct Assignment<DstXprType, SparseQRMatrixQReturnType<SparseQRType>, internal::assign_op<typename DstXprType::Scalar,typename DstXprType::Scalar>, Sparse2Dense>
732 typedef SparseQRMatrixQReturnType<SparseQRType> SrcXprType;
733 typedef typename DstXprType::Scalar Scalar;
734 typedef typename DstXprType::StorageIndex StorageIndex;
735 static void run(DstXprType &dst,
const SrcXprType &src,
const internal::assign_op<Scalar,Scalar> &)
737 dst = src.m_qr.matrixQ() * DstXprType::Identity(src.m_qr.rows(), src.m_qr.rows());
Index cols() const
Definition: SparseMatrix.h:138
Index cols() const
Definition: SparseQR.h:129
SparseQRMatrixQReturnType< SparseQR > matrixQ() const
Definition: SparseQR.h:174
Index rows() const
Definition: SparseMatrix.h:136
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_sqrt_op< typename Derived::Scalar >, const Derived > sqrt(const Eigen::ArrayBase< Derived > &x)
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SparseQR.h:256
Namespace containing all symbols from the Eigen library.
Definition: Core:306
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:150
Derived & derived()
Definition: EigenBase.h:45
Index rows() const
Definition: SparseQR.h:125
const QRMatrixType & matrixR() const
Definition: SparseQR.h:144
void factorize(const MatrixType &mat)
Performs the numerical QR factorization of the input matrix.
Definition: SparseQR.h:349
Index size() const
Definition: PermutationMatrix.h:108
Derived & setZero(Index size)
Definition: CwiseNullaryOp.h:515
Sparse left-looking rank-revealing QR factorization.
Definition: SparseQR.h:16
Index rows() const
Definition: SparseMatrixBase.h:171
Base class of any sparse matrices or sparse expressions.
Definition: ForwardDeclarations.h:281
Derived & derived()
Definition: EigenBase.h:45
const Product< MatrixDerived, PermutationDerived, AliasFreeProduct > operator*(const MatrixBase< MatrixDerived > &matrix, const PermutationBase< PermutationDerived > &permutation)
Definition: PermutationMatrix.h:543
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
BlockXpr topLeftCorner(Index cRows, Index cCols)
Definition: SparseMatrixBase.h:179
Derived & setConstant(Index size, const Scalar &val)
Definition: CwiseNullaryOp.h:341
Definition: Constants.h:439
Definition: Constants.h:432
const Scalar * data() const
Definition: PlainObjectBase.h:255
void analyzePattern(const MatrixType &mat)
Preprocessing step of a QR factorization.
Definition: SparseQR.h:308
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
Index rank() const
Definition: SparseQR.h:150
Index rows() const
Definition: EigenBase.h:59
const int Dynamic
Definition: Constants.h:21
Pseudo expression representing a solving operation.
Definition: Solve.h:62
const PermutationType & colsPermutation() const
Definition: SparseQR.h:180
std::string lastErrorMessage() const
Definition: SparseQR.h:189
ComputationInfo
Definition: Constants.h:430
SparseQR(const MatrixType &mat)
Definition: SparseQR.h:104
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
void setPivotThreshold(const RealScalar &threshold)
Definition: SparseQR.h:223
const Solve< SparseQR, Rhs > solve(const MatrixBase< Rhs > &B) const
Definition: SparseQR.h:234
void compute(const MatrixType &mat)
Definition: SparseQR.h:115