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;
165 y.bottomRows(y.size()-
rank).setZero();
183 m_useDefaultThreshold =
false;
184 m_threshold = threshold;
191 template<
typename Rhs>
194 eigen_assert(m_isInitialized &&
"The factorization should be called first, use compute()");
195 eigen_assert(this->
rows() == B.rows() &&
"SparseQR::solve() : invalid number of rows in the right hand side matrix");
196 return internal::solve_retval<SparseQR, Rhs>(*
this, B.derived());
198 template<
typename Rhs>
201 eigen_assert(m_isInitialized &&
"The factorization should be called first, use compute()");
202 eigen_assert(this->
rows() == B.
rows() &&
"SparseQR::solve() : invalid number of rows in the right hand side matrix");
203 return internal::sparse_solve_retval<SparseQR, Rhs>(*
this, B.
derived());
216 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
221 inline void sort_matrix_Q()
223 if(this->m_isQSorted)
return;
227 this->m_isQSorted =
true;
232 bool m_isInitialized;
234 bool m_factorizationIsok;
236 std::string m_lastError;
240 ScalarVector m_hcoeffs;
241 PermutationType m_perm_c;
242 PermutationType m_pivotperm;
243 PermutationType m_outputPerm_c;
244 RealScalar m_threshold;
245 bool m_useDefaultThreshold;
246 Index m_nonzeropivots;
248 IndexVector m_firstRowElt;
251 template <
typename,
typename >
friend struct SparseQR_QProduct;
252 template <
typename >
friend struct SparseQRMatrixQReturnType;
263 template <
typename MatrixType,
typename OrderingType>
269 Index n = mat.cols();
270 Index m = mat.rows();
272 if (!m_perm_c.size())
275 m_perm_c.indices().setLinSpaced(n, 0,n-1);
279 m_outputPerm_c = m_perm_c.inverse();
286 m_R.reserve(2*mat.nonZeros());
287 m_Q.reserve(2*mat.nonZeros());
289 m_analysisIsok =
true;
299 template <
typename MatrixType,
typename OrderingType>
305 eigen_assert(m_analysisIsok &&
"analyzePattern() should be called before this step");
306 Index m = mat.rows();
307 Index n = mat.cols();
310 Index nzcolR, nzcolQ;
317 for (
int i = 0; i < n; i++)
319 Index p = m_perm_c.size() ? m_perm_c.indices()(i) : i;
320 m_pmat.outerIndexPtr()[p] = mat.outerIndexPtr()[i];
321 m_pmat.innerNonZeroPtr()[p] = mat.outerIndexPtr()[i+1] - mat.outerIndexPtr()[i];
328 if(m_useDefaultThreshold)
330 RealScalar max2Norm = 0.0;
331 for (
int j = 0; j < n; j++) max2Norm = (max)(max2Norm, m_pmat.col(j).norm());
336 m_pivotperm.setIdentity(n);
338 Index nonzeroCol = 0;
341 for (Index col = 0; col < n; ++col)
346 mark(nonzeroCol) = col;
347 Qidx(0) = nonzeroCol;
348 nzcolR = 0; nzcolQ = 1;
356 for (
typename MatrixType::InnerIterator itp(m_pmat, col); itp || !found_diag; ++itp)
358 Index curIdx = nonzeroCol ;
359 if(itp) curIdx = itp.row();
360 if(curIdx == nonzeroCol) found_diag =
true;
363 Index st = m_firstRowElt(curIdx);
366 m_lastError =
"Empty row found during numerical factorization";
373 for (; mark(st) != col; st = m_etree(st))
381 Index nt = nzcolR-bi;
382 for(Index i = 0; i < nt/2; i++) std::swap(Ridx(bi+i), Ridx(nzcolR-i-1));
385 if(itp) tval(curIdx) = itp.value();
386 else tval(curIdx) = Scalar(0);
389 if(curIdx > nonzeroCol && mark(curIdx) != col )
391 Qidx(nzcolQ) = curIdx;
398 for (Index i = nzcolR-1; i >= 0; i--)
400 Index curIdx = m_pivotperm.indices()(Ridx(i));
406 tdot = m_Q.col(curIdx).dot(tval);
408 tdot *= m_hcoeffs(curIdx);
412 for (
typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq)
413 tval(itq.row()) -= itq.value() * tdot;
416 if(m_etree(Ridx(i)) == nonzeroCol)
418 for (
typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq)
420 Index iQ = itq.row();
434 Scalar c0 = nzcolQ ? tval(Qidx(0)) : Scalar(0);
437 RealScalar sqrNorm = 0.;
438 for (Index itq = 1; itq < nzcolQ; ++itq) sqrNorm += numext::abs2(tval(Qidx(itq)));
440 if(sqrNorm == RealScalar(0) && numext::imag(c0) == RealScalar(0))
443 beta = numext::real(c0);
448 beta = std::sqrt(numext::abs2(c0) + sqrNorm);
449 if(numext::real(c0) >= RealScalar(0))
452 for (Index itq = 1; itq < nzcolQ; ++itq)
453 tval(Qidx(itq)) /= (c0 - beta);
454 tau = numext::conj((beta-c0) / beta);
459 for (Index i = nzcolR-1; i >= 0; i--)
461 Index curIdx = Ridx(i);
462 if(curIdx < nonzeroCol)
464 m_R.insertBackByOuterInnerUnordered(col, curIdx) = tval(curIdx);
465 tval(curIdx) = Scalar(0.);
469 if(abs(beta) >= m_threshold)
471 m_R.insertBackByOuterInner(col, nonzeroCol) = beta;
474 m_hcoeffs(col) = tau;
476 for (Index itq = 0; itq < nzcolQ; ++itq)
478 Index iQ = Qidx(itq);
479 m_Q.insertBackByOuterInnerUnordered(col,iQ) = tval(iQ);
480 tval(iQ) = Scalar(0.);
486 m_hcoeffs(col) = Scalar(0);
487 for (Index j = nonzeroCol; j < n-1; j++)
488 std::swap(m_pivotperm.indices()(j), m_pivotperm.indices()[j+1]);
497 m_Q.makeCompressed();
499 m_R.makeCompressed();
502 m_nonzeropivots = nonzeroCol;
507 MatrixType tempR(m_R);
508 m_R = tempR * m_pivotperm;
511 m_outputPerm_c = m_outputPerm_c * m_pivotperm;
514 m_isInitialized =
true;
515 m_factorizationIsok =
true;
521 template<
typename _MatrixType,
typename OrderingType,
typename Rhs>
522 struct solve_retval<
SparseQR<_MatrixType,OrderingType>, Rhs>
523 : solve_retval_base<SparseQR<_MatrixType,OrderingType>, Rhs>
526 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
528 template<typename Dest>
void evalTo(Dest& dst)
const
530 dec()._solve(rhs(),dst);
533 template<
typename _MatrixType,
typename OrderingType,
typename Rhs>
534 struct sparse_solve_retval<SparseQR<_MatrixType, OrderingType>, Rhs>
535 : sparse_solve_retval_base<SparseQR<_MatrixType, OrderingType>, Rhs>
537 typedef SparseQR<_MatrixType, OrderingType> Dec;
538 EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec, Rhs)
540 template<typename Dest>
void evalTo(Dest& dst)
const
542 this->defaultEvalTo(dst);
547 template <
typename SparseQRType,
typename Derived>
548 struct SparseQR_QProduct : ReturnByValue<SparseQR_QProduct<SparseQRType, Derived> >
550 typedef typename SparseQRType::QRMatrixType MatrixType;
551 typedef typename SparseQRType::Scalar Scalar;
552 typedef typename SparseQRType::Index Index;
554 SparseQR_QProduct(
const SparseQRType& qr,
const Derived& other,
bool transpose) :
555 m_qr(qr),m_other(other),m_transpose(transpose) {}
556 inline Index rows()
const {
return m_transpose ? m_qr.rows() : m_qr.cols(); }
557 inline Index cols()
const {
return m_other.cols(); }
560 template<
typename DesType>
561 void evalTo(DesType& res)
const
563 Index n = m_qr.cols();
567 eigen_assert(m_qr.m_Q.rows() == m_other.rows() &&
"Non conforming object sizes");
569 for(Index j = 0; j < res.cols(); j++){
570 for (Index k = 0; k < n; k++)
572 Scalar tau = Scalar(0);
573 tau = m_qr.m_Q.col(k).dot(res.col(j));
574 tau = tau * m_qr.m_hcoeffs(k);
575 res.col(j) -= tau * m_qr.m_Q.col(k);
581 eigen_assert(m_qr.m_Q.cols() == m_other.rows() &&
"Non conforming object sizes");
583 for(Index j = 0; j < res.cols(); j++)
585 for (Index k = n-1; k >=0; k--)
587 Scalar tau = Scalar(0);
588 tau = m_qr.m_Q.col(k).dot(res.col(j));
589 tau = tau * m_qr.m_hcoeffs(k);
590 res.col(j) -= tau * m_qr.m_Q.col(k);
596 const SparseQRType& m_qr;
597 const Derived& m_other;
601 template<
typename SparseQRType>
602 struct SparseQRMatrixQReturnType :
public EigenBase<SparseQRMatrixQReturnType<SparseQRType> >
604 typedef typename SparseQRType::Index Index;
605 typedef typename SparseQRType::Scalar Scalar;
606 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
607 SparseQRMatrixQReturnType(
const SparseQRType& qr) : m_qr(qr) {}
608 template<
typename Derived>
609 SparseQR_QProduct<SparseQRType, Derived>
operator*(
const MatrixBase<Derived>& other)
611 return SparseQR_QProduct<SparseQRType,Derived>(m_qr,other.derived(),
false);
613 SparseQRMatrixQTransposeReturnType<SparseQRType> adjoint()
const
615 return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr);
617 inline Index rows()
const {
return m_qr.rows(); }
618 inline Index cols()
const {
return m_qr.cols(); }
620 SparseQRMatrixQTransposeReturnType<SparseQRType> transpose()
const
622 return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr);
624 template<
typename Dest>
void evalTo(MatrixBase<Dest>& dest)
const
626 dest.derived() = m_qr.matrixQ() * Dest::Identity(m_qr.rows(), m_qr.rows());
628 template<
typename Dest>
void evalTo(SparseMatrixBase<Dest>& dest)
const
630 Dest idMat(m_qr.rows(), m_qr.rows());
633 const_cast<SparseQRType *
>(&m_qr)->sort_matrix_Q();
634 dest.derived() = SparseQR_QProduct<SparseQRType, Dest>(m_qr, idMat,
false);
637 const SparseQRType& m_qr;
640 template<
typename SparseQRType>
641 struct SparseQRMatrixQTransposeReturnType
643 SparseQRMatrixQTransposeReturnType(
const SparseQRType& qr) : m_qr(qr) {}
644 template<
typename Derived>
645 SparseQR_QProduct<SparseQRType,Derived>
operator*(
const MatrixBase<Derived>& other)
647 return SparseQR_QProduct<SparseQRType,Derived>(m_qr,other.derived(),
true);
649 const SparseQRType& m_qr;
Index rows() const
Definition: SparseMatrix.h:119
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:300
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:214
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:181
const PermutationType & colsPermutation() const
Definition: SparseQR.h:138
const internal::solve_retval< SparseQR, Rhs > solve(const MatrixBase< Rhs > &B) const
Definition: SparseQR.h:192
const QRMatrixType & matrixR() const
Definition: SparseQR.h:102
void analyzePattern(const MatrixType &mat)
Preprocessing step of a QR factorization.
Definition: SparseQR.h:264
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