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::Index Index;
25 typedef typename ReturnType::StorageKind StorageKind;
27 template <
typename SparseQRType>
struct traits<SparseQRMatrixQTransposeReturnType<SparseQRType> >
29 typedef typename SparseQRType::MatrixType ReturnType;
31 template <
typename SparseQRType,
typename Derived>
struct traits<SparseQR_QProduct<SparseQRType, Derived> >
33 typedef typename Derived::PlainObject ReturnType;
63 template<
typename _MatrixType,
typename _OrderingType>
67 typedef _MatrixType MatrixType;
68 typedef _OrderingType OrderingType;
69 typedef typename MatrixType::Scalar Scalar;
70 typedef typename MatrixType::RealScalar RealScalar;
71 typedef typename MatrixType::Index Index;
72 typedef SparseMatrix<Scalar,ColMajor,Index> QRMatrixType;
73 typedef Matrix<Index, Dynamic, 1> IndexVector;
74 typedef Matrix<Scalar, Dynamic, 1> ScalarVector;
75 typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType;
77 SparseQR () : m_isInitialized(false), m_analysisIsok(false), m_lastError(
""), m_useDefaultThreshold(true),m_isQSorted(false)
80 SparseQR(
const MatrixType& mat) : m_isInitialized(false), m_analysisIsok(false), m_lastError(
""), m_useDefaultThreshold(true),m_isQSorted(false)
84 void compute(
const MatrixType& mat)
94 inline Index
rows()
const {
return m_pmat.
rows(); }
98 inline Index
cols()
const {
return m_pmat.
cols();}
110 eigen_assert(m_isInitialized &&
"The factorization should be called first, use compute()");
111 return m_nonzeropivots;
132 SparseQRMatrixQReturnType<SparseQR>
matrixQ()
const
133 {
return SparseQRMatrixQReturnType<SparseQR>(*this); }
140 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
141 return m_outputPerm_c;
150 template<
typename Rhs,
typename Dest>
153 eigen_assert(m_isInitialized &&
"The factorization should be called first, use compute()");
154 eigen_assert(this->
rows() == B.rows() &&
"SparseQR::solve() : invalid number of rows in the right hand side matrix");
159 typename Dest::PlainObject y, b;
160 y = this->
matrixQ().transpose() * B;
164 y.
resize((std::max)(
cols(),Index(y.rows())),y.cols());
166 y.bottomRows(y.rows()-
rank).setZero();
184 m_useDefaultThreshold =
false;
185 m_threshold = threshold;
192 template<
typename Rhs>
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 return internal::solve_retval<SparseQR, Rhs>(*
this, B.derived());
199 template<
typename Rhs>
202 eigen_assert(m_isInitialized &&
"The factorization should be called first, use compute()");
203 eigen_assert(this->
rows() == B.
rows() &&
"SparseQR::solve() : invalid number of rows in the right hand side matrix");
204 return internal::sparse_solve_retval<SparseQR, Rhs>(*
this, B.
derived());
217 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
222 inline void sort_matrix_Q()
224 if(this->m_isQSorted)
return;
228 this->m_isQSorted =
true;
233 bool m_isInitialized;
235 bool m_factorizationIsok;
237 std::string m_lastError;
241 ScalarVector m_hcoeffs;
242 PermutationType m_perm_c;
243 PermutationType m_pivotperm;
244 PermutationType m_outputPerm_c;
245 RealScalar m_threshold;
246 bool m_useDefaultThreshold;
247 Index m_nonzeropivots;
249 IndexVector m_firstRowElt;
252 template <
typename,
typename >
friend struct SparseQR_QProduct;
253 template <
typename >
friend struct SparseQRMatrixQReturnType;
264 template <
typename MatrixType,
typename OrderingType>
270 Index n = mat.cols();
271 Index m = mat.rows();
273 if (!m_perm_c.size())
276 m_perm_c.indices().setLinSpaced(n, 0,n-1);
280 m_outputPerm_c = m_perm_c.inverse();
287 m_R.reserve(2*mat.nonZeros());
288 m_Q.reserve(2*mat.nonZeros());
290 m_analysisIsok =
true;
300 template <
typename MatrixType,
typename OrderingType>
306 eigen_assert(m_analysisIsok &&
"analyzePattern() should be called before this step");
307 Index m = mat.rows();
308 Index n = mat.cols();
311 Index nzcolR, nzcolQ;
318 for (
int i = 0; i < n; i++)
320 Index p = m_perm_c.size() ? m_perm_c.indices()(i) : i;
321 m_pmat.outerIndexPtr()[p] = mat.outerIndexPtr()[i];
322 m_pmat.innerNonZeroPtr()[p] = mat.outerIndexPtr()[i+1] - mat.outerIndexPtr()[i];
329 if(m_useDefaultThreshold)
331 RealScalar max2Norm = 0.0;
332 for (
int j = 0; j < n; j++) max2Norm = (max)(max2Norm, m_pmat.col(j).norm());
337 m_pivotperm.setIdentity(n);
339 Index nonzeroCol = 0;
342 for (Index col = 0; col < (std::min)(n,m); ++col)
347 mark(nonzeroCol) = col;
348 Qidx(0) = nonzeroCol;
349 nzcolR = 0; nzcolQ = 1;
357 for (
typename MatrixType::InnerIterator itp(m_pmat, col); itp || !found_diag; ++itp)
359 Index curIdx = nonzeroCol ;
360 if(itp) curIdx = itp.row();
361 if(curIdx == nonzeroCol) found_diag =
true;
364 Index st = m_firstRowElt(curIdx);
367 m_lastError =
"Empty row found during numerical factorization";
374 for (; mark(st) != col; st = m_etree(st))
382 Index nt = nzcolR-bi;
383 for(Index i = 0; i < nt/2; i++) std::swap(Ridx(bi+i), Ridx(nzcolR-i-1));
386 if(itp) tval(curIdx) = itp.value();
387 else tval(curIdx) = Scalar(0);
390 if(curIdx > nonzeroCol && mark(curIdx) != col )
392 Qidx(nzcolQ) = curIdx;
399 for (Index i = nzcolR-1; i >= 0; i--)
401 Index curIdx = m_pivotperm.indices()(Ridx(i));
407 tdot = m_Q.col(curIdx).dot(tval);
409 tdot *= m_hcoeffs(curIdx);
413 for (
typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq)
414 tval(itq.row()) -= itq.value() * tdot;
417 if(m_etree(Ridx(i)) == nonzeroCol)
419 for (
typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq)
421 Index iQ = itq.row();
435 Scalar c0 = nzcolQ ? tval(Qidx(0)) : Scalar(0);
438 RealScalar sqrNorm = 0.;
439 for (Index itq = 1; itq < nzcolQ; ++itq) sqrNorm += numext::abs2(tval(Qidx(itq)));
441 if(sqrNorm == RealScalar(0) && numext::imag(c0) == RealScalar(0))
444 beta = numext::real(c0);
449 beta = std::sqrt(numext::abs2(c0) + sqrNorm);
450 if(numext::real(c0) >= RealScalar(0))
453 for (Index itq = 1; itq < nzcolQ; ++itq)
454 tval(Qidx(itq)) /= (c0 - beta);
455 tau = numext::conj((beta-c0) / beta);
460 for (Index i = nzcolR-1; i >= 0; i--)
462 Index curIdx = Ridx(i);
463 if(curIdx < nonzeroCol)
465 m_R.insertBackByOuterInnerUnordered(col, curIdx) = tval(curIdx);
466 tval(curIdx) = Scalar(0.);
470 if(abs(beta) >= m_threshold)
472 m_R.insertBackByOuterInner(col, nonzeroCol) = beta;
475 m_hcoeffs(col) = tau;
477 for (Index itq = 0; itq < nzcolQ; ++itq)
479 Index iQ = Qidx(itq);
480 m_Q.insertBackByOuterInnerUnordered(col,iQ) = tval(iQ);
481 tval(iQ) = Scalar(0.);
487 m_hcoeffs(col) = Scalar(0);
488 for (Index j = nonzeroCol; j < n-1; j++)
489 std::swap(m_pivotperm.indices()(j), m_pivotperm.indices()[j+1]);
498 m_Q.makeCompressed();
500 m_R.makeCompressed();
503 m_nonzeropivots = nonzeroCol;
508 MatrixType tempR(m_R);
509 m_R = tempR * m_pivotperm;
512 m_outputPerm_c = m_outputPerm_c * m_pivotperm;
515 m_isInitialized =
true;
516 m_factorizationIsok =
true;
522 template<
typename _MatrixType,
typename OrderingType,
typename Rhs>
523 struct solve_retval<
SparseQR<_MatrixType,OrderingType>, Rhs>
524 : solve_retval_base<SparseQR<_MatrixType,OrderingType>, Rhs>
527 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
529 template<typename Dest>
void evalTo(Dest& dst)
const
531 dec()._solve(rhs(),dst);
534 template<
typename _MatrixType,
typename OrderingType,
typename Rhs>
535 struct sparse_solve_retval<SparseQR<_MatrixType, OrderingType>, Rhs>
536 : sparse_solve_retval_base<SparseQR<_MatrixType, OrderingType>, Rhs>
538 typedef SparseQR<_MatrixType, OrderingType> Dec;
539 EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec, Rhs)
541 template<typename Dest>
void evalTo(Dest& dst)
const
543 this->defaultEvalTo(dst);
548 template <
typename SparseQRType,
typename Derived>
549 struct SparseQR_QProduct : ReturnByValue<SparseQR_QProduct<SparseQRType, Derived> >
551 typedef typename SparseQRType::QRMatrixType MatrixType;
552 typedef typename SparseQRType::Scalar Scalar;
553 typedef typename SparseQRType::Index Index;
555 SparseQR_QProduct(
const SparseQRType& qr,
const Derived& other,
bool transpose) :
556 m_qr(qr),m_other(other),m_transpose(transpose) {}
557 inline Index rows()
const {
return m_transpose ? m_qr.rows() : m_qr.cols(); }
558 inline Index cols()
const {
return m_other.cols(); }
561 template<
typename DesType>
562 void evalTo(DesType& res)
const
564 Index n = m_qr.cols();
568 eigen_assert(m_qr.m_Q.rows() == m_other.rows() &&
"Non conforming object sizes");
570 for(Index j = 0; j < res.cols(); j++){
571 for (Index k = 0; k < n; k++)
573 Scalar tau = Scalar(0);
574 tau = m_qr.m_Q.col(k).dot(res.col(j));
575 if(tau==Scalar(0))
continue;
576 tau = tau * m_qr.m_hcoeffs(k);
577 res.col(j) -= tau * m_qr.m_Q.col(k);
583 eigen_assert(m_qr.m_Q.rows() == m_other.rows() &&
"Non conforming object sizes");
585 for(Index j = 0; j < res.cols(); j++)
587 for (Index k = n-1; k >=0; k--)
589 Scalar tau = Scalar(0);
590 tau = m_qr.m_Q.col(k).dot(res.col(j));
591 if(tau==Scalar(0))
continue;
592 tau = tau * m_qr.m_hcoeffs(k);
593 res.col(j) -= tau * m_qr.m_Q.col(k);
599 const SparseQRType& m_qr;
600 const Derived& m_other;
604 template<
typename SparseQRType>
605 struct SparseQRMatrixQReturnType :
public EigenBase<SparseQRMatrixQReturnType<SparseQRType> >
607 typedef typename SparseQRType::Index Index;
608 typedef typename SparseQRType::Scalar Scalar;
609 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
610 SparseQRMatrixQReturnType(
const SparseQRType& qr) : m_qr(qr) {}
611 template<
typename Derived>
612 SparseQR_QProduct<SparseQRType, Derived>
operator*(
const MatrixBase<Derived>& other)
614 return SparseQR_QProduct<SparseQRType,Derived>(m_qr,other.derived(),
false);
616 SparseQRMatrixQTransposeReturnType<SparseQRType> adjoint()
const
618 return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr);
620 inline Index rows()
const {
return m_qr.rows(); }
621 inline Index cols()
const {
return m_qr.cols(); }
623 SparseQRMatrixQTransposeReturnType<SparseQRType> transpose()
const
625 return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr);
627 template<
typename Dest>
void evalTo(MatrixBase<Dest>& dest)
const
629 dest.derived() = m_qr.matrixQ() * Dest::Identity(m_qr.rows(), m_qr.rows());
631 template<
typename Dest>
void evalTo(SparseMatrixBase<Dest>& dest)
const
633 Dest idMat(m_qr.rows(), m_qr.rows());
636 const_cast<SparseQRType *
>(&m_qr)->sort_matrix_Q();
637 dest.derived() = SparseQR_QProduct<SparseQRType, Dest>(m_qr, idMat,
false);
640 const SparseQRType& m_qr;
643 template<
typename SparseQRType>
644 struct SparseQRMatrixQTransposeReturnType
646 SparseQRMatrixQTransposeReturnType(
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(),
true);
652 const SparseQRType& m_qr;
Index rows() const
Definition: SparseMatrix.h:119
void resize(Index newSize)
Definition: DenseBase.h:217
std::string lastErrorMessage() const
Definition: SparseQR.h:147
Index cols() const
Definition: SparseMatrix.h:121
Index rank() const
Definition: SparseQR.h:108
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:88
const internal::permut_matrix_product_retval< PermutationDerived, Derived, OnTheRight > operator*(const MatrixBase< Derived > &matrix, const PermutationBase< PermutationDerived > &permutation)
Definition: PermutationMatrix.h:510
void factorize(const MatrixType &mat)
Performs the numerical QR factorization of the input matrix.
Definition: SparseQR.h:301
SparseQRMatrixQReturnType< SparseQR > matrixQ() const
Definition: SparseQR.h:132
int coletree(const MatrixType &mat, IndexVector &parent, IndexVector &firstRowElt, typename MatrixType::Index *perm=0)
Definition: SparseColEtree.h:61
Sparse left-looking rank-revealing QR factorization.
Definition: SparseQR.h:16
Index size() const
Definition: PermutationMatrix.h:114
Base class of any sparse matrices or sparse expressions.
Definition: SparseMatrixBase.h:26
Index rows() const
Definition: SparseQR.h:94
Derived & derived()
Definition: EigenBase.h:34
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SparseQR.h:215
Index cols() const
Definition: SparseQR.h:98
Derived & setConstant(Index size, const Scalar &value)
Definition: CwiseNullaryOp.h:348
RowsBlockXpr topRows(Index n)
Definition: DenseBase.h:381
Definition: Constants.h:383
void setPivotThreshold(const RealScalar &threshold)
Definition: SparseQR.h:182
const PermutationType & colsPermutation() const
Definition: SparseQR.h:138
const internal::solve_retval< SparseQR, Rhs > solve(const MatrixBase< Rhs > &B) const
Definition: SparseQR.h:193
const QRMatrixType & matrixR() const
Definition: SparseQR.h:102
void analyzePattern(const MatrixType &mat)
Preprocessing step of a QR factorization.
Definition: SparseQR.h:265
Definition: Constants.h:376
Block< Derived > topLeftCorner(Index cRows, Index cCols)
Definition: SparseMatrixBase.h:157
Index rows() const
Definition: SparseMatrixBase.h:150
ComputationInfo
Definition: Constants.h:374
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
Derived & setZero(Index size)
Definition: CwiseNullaryOp.h:515