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;
64 template<
typename _MatrixType,
typename _OrderingType>
68 typedef _MatrixType MatrixType;
69 typedef _OrderingType OrderingType;
70 typedef typename MatrixType::Scalar Scalar;
71 typedef typename MatrixType::RealScalar RealScalar;
72 typedef typename MatrixType::Index Index;
78 SparseQR () : m_isInitialized(false), m_analysisIsok(false), m_lastError(
""), m_useDefaultThreshold(true),m_isQSorted(false)
87 SparseQR(
const MatrixType& mat) : m_isInitialized(false), m_analysisIsok(false), m_lastError(
""), m_useDefaultThreshold(true),m_isQSorted(false)
103 void analyzePattern(
const MatrixType& mat);
104 void factorize(
const MatrixType& mat);
108 inline Index
rows()
const {
return m_pmat.rows(); }
112 inline Index
cols()
const {
return m_pmat.cols();}
116 const QRMatrixType&
matrixR()
const {
return m_R; }
124 eigen_assert(m_isInitialized &&
"The factorization should be called first, use compute()");
125 return m_nonzeropivots;
146 SparseQRMatrixQReturnType<SparseQR>
matrixQ()
const 147 {
return SparseQRMatrixQReturnType<SparseQR>(*this); }
154 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
155 return m_outputPerm_c;
164 template<
typename Rhs,
typename Dest>
167 eigen_assert(m_isInitialized &&
"The factorization should be called first, use compute()");
168 eigen_assert(this->rows() == B.rows() &&
"SparseQR::solve() : invalid number of rows in the right hand side matrix");
170 Index rank = this->rank();
173 typename Dest::PlainObject y, b;
174 y = this->matrixQ().transpose() * B;
178 y.
resize((std::max)(cols(),Index(y.rows())),y.cols());
179 y.topRows(rank) = this->matrixR().topLeftCorner(rank, rank).template triangularView<Upper>().solve(b.topRows(rank));
180 y.bottomRows(y.rows()-rank).setZero();
183 if (m_perm_c.size()) dest = colsPermutation() * y.
topRows(cols());
198 m_useDefaultThreshold =
false;
199 m_threshold = threshold;
206 template<
typename Rhs>
209 eigen_assert(m_isInitialized &&
"The factorization should be called first, use compute()");
210 eigen_assert(this->rows() == B.rows() &&
"SparseQR::solve() : invalid number of rows in the right hand side matrix");
211 return internal::solve_retval<SparseQR, Rhs>(*
this, B.derived());
213 template<
typename Rhs>
214 inline const internal::sparse_solve_retval<SparseQR, Rhs> solve(
const SparseMatrixBase<Rhs>& B)
const 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");
218 return internal::sparse_solve_retval<SparseQR, Rhs>(*
this, B.
derived());
231 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
236 inline void sort_matrix_Q()
238 if(this->m_isQSorted)
return;
242 this->m_isQSorted =
true;
247 bool m_isInitialized;
249 bool m_factorizationIsok;
251 std::string m_lastError;
255 ScalarVector m_hcoeffs;
256 PermutationType m_perm_c;
257 PermutationType m_pivotperm;
258 PermutationType m_outputPerm_c;
259 RealScalar m_threshold;
260 bool m_useDefaultThreshold;
261 Index m_nonzeropivots;
263 IndexVector m_firstRowElt;
266 template <
typename,
typename >
friend struct SparseQR_QProduct;
267 template <
typename >
friend struct SparseQRMatrixQReturnType;
280 template <
typename MatrixType,
typename OrderingType>
283 eigen_assert(mat.isCompressed() &&
"SparseQR requires a sparse matrix in compressed mode. Call .makeCompressed() before passing it to SparseQR");
287 Index n = mat.cols();
288 Index m = mat.rows();
289 Index diagSize = (std::min)(m,n);
291 if (!m_perm_c.size())
294 m_perm_c.indices().setLinSpaced(n, 0,n-1);
298 m_outputPerm_c = m_perm_c.inverse();
302 m_Q.resize(m, diagSize);
305 m_R.reserve(2*mat.nonZeros());
306 m_Q.reserve(2*mat.nonZeros());
307 m_hcoeffs.resize(diagSize);
308 m_analysisIsok =
true;
318 template <
typename MatrixType,
typename OrderingType>
324 eigen_assert(m_analysisIsok &&
"analyzePattern() should be called before this step");
325 Index m = mat.rows();
326 Index n = mat.cols();
327 Index diagSize = (std::min)(m,n);
330 Index nzcolR, nzcolQ;
332 RealScalar pivotThreshold = m_threshold;
337 for (
int i = 0; i < n; i++)
339 Index p = m_perm_c.size() ? m_perm_c.indices()(i) : i;
340 m_pmat.outerIndexPtr()[p] = mat.outerIndexPtr()[i];
341 m_pmat.innerNonZeroPtr()[p] = mat.outerIndexPtr()[i+1] - mat.outerIndexPtr()[i];
348 if(m_useDefaultThreshold)
350 RealScalar max2Norm = 0.0;
351 for (
int j = 0; j < n; j++) max2Norm = (max)(max2Norm, m_pmat.col(j).norm());
356 m_pivotperm.setIdentity(n);
358 Index nonzeroCol = 0;
362 for (Index col = 0; col < n; ++col)
366 mark(nonzeroCol) = col;
367 Qidx(0) = nonzeroCol;
368 nzcolR = 0; nzcolQ = 1;
369 bool found_diag = nonzeroCol>=m;
376 for (
typename MatrixType::InnerIterator itp(m_pmat, col); itp || !found_diag; ++itp)
378 Index curIdx = nonzeroCol;
379 if(itp) curIdx = itp.row();
380 if(curIdx == nonzeroCol) found_diag =
true;
383 Index st = m_firstRowElt(curIdx);
386 m_lastError =
"Empty row found during numerical factorization";
393 for (; mark(st) != col; st = m_etree(st))
401 Index nt = nzcolR-bi;
402 for(Index i = 0; i < nt/2; i++) std::swap(Ridx(bi+i), Ridx(nzcolR-i-1));
405 if(itp) tval(curIdx) = itp.value();
406 else tval(curIdx) = Scalar(0);
409 if(curIdx > nonzeroCol && mark(curIdx) != col )
411 Qidx(nzcolQ) = curIdx;
418 for (Index i = nzcolR-1; i >= 0; i--)
420 Index curIdx = Ridx(i);
426 tdot = m_Q.col(curIdx).dot(tval);
428 tdot *= m_hcoeffs(curIdx);
432 for (
typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq)
433 tval(itq.row()) -= itq.value() * tdot;
436 if(m_etree(Ridx(i)) == nonzeroCol)
438 for (
typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq)
440 Index iQ = itq.row();
453 if(nonzeroCol < diagSize)
457 Scalar c0 = nzcolQ ? tval(Qidx(0)) : Scalar(0);
460 RealScalar sqrNorm = 0.;
461 for (Index itq = 1; itq < nzcolQ; ++itq) sqrNorm += numext::abs2(tval(Qidx(itq)));
462 if(sqrNorm == RealScalar(0) && numext::imag(c0) == RealScalar(0))
465 beta = numext::real(c0);
471 beta = sqrt(numext::abs2(c0) + sqrNorm);
472 if(numext::real(c0) >= RealScalar(0))
475 for (Index itq = 1; itq < nzcolQ; ++itq)
476 tval(Qidx(itq)) /= (c0 - beta);
477 tau = numext::conj((beta-c0) / beta);
483 for (Index i = nzcolR-1; i >= 0; i--)
485 Index curIdx = Ridx(i);
486 if(curIdx < nonzeroCol)
488 m_R.insertBackByOuterInnerUnordered(col, curIdx) = tval(curIdx);
489 tval(curIdx) = Scalar(0.);
493 if(nonzeroCol < diagSize && abs(beta) >= pivotThreshold)
495 m_R.insertBackByOuterInner(col, nonzeroCol) = beta;
497 m_hcoeffs(nonzeroCol) = tau;
499 for (Index itq = 0; itq < nzcolQ; ++itq)
501 Index iQ = Qidx(itq);
502 m_Q.insertBackByOuterInnerUnordered(nonzeroCol,iQ) = tval(iQ);
503 tval(iQ) = Scalar(0.);
506 if(nonzeroCol<diagSize)
507 m_Q.startVec(nonzeroCol);
512 for (Index j = nonzeroCol; j < n-1; j++)
513 std::swap(m_pivotperm.indices()(j), m_pivotperm.indices()[j+1]);
520 m_hcoeffs.tail(diagSize-nonzeroCol).setZero();
524 m_Q.makeCompressed();
526 m_R.makeCompressed();
529 m_nonzeropivots = nonzeroCol;
534 MatrixType tempR(m_R);
535 m_R = tempR * m_pivotperm;
538 m_outputPerm_c = m_outputPerm_c * m_pivotperm;
541 m_isInitialized =
true;
542 m_factorizationIsok =
true;
548 template<
typename _MatrixType,
typename OrderingType,
typename Rhs>
549 struct solve_retval<SparseQR<_MatrixType,OrderingType>, Rhs>
550 : solve_retval_base<SparseQR<_MatrixType,OrderingType>, Rhs>
553 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
555 template<
typename Dest>
void evalTo(Dest& dst)
const 557 dec()._solve(rhs(),dst);
560 template<
typename _MatrixType,
typename OrderingType,
typename Rhs>
561 struct sparse_solve_retval<SparseQR<_MatrixType, OrderingType>, Rhs>
562 : sparse_solve_retval_base<SparseQR<_MatrixType, OrderingType>, Rhs>
565 EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec, Rhs)
567 template<
typename Dest>
void evalTo(Dest& dst)
const 569 this->defaultEvalTo(dst);
574 template <
typename SparseQRType,
typename Derived>
575 struct SparseQR_QProduct : ReturnByValue<SparseQR_QProduct<SparseQRType, Derived> >
577 typedef typename SparseQRType::QRMatrixType MatrixType;
578 typedef typename SparseQRType::Scalar Scalar;
579 typedef typename SparseQRType::Index Index;
581 SparseQR_QProduct(
const SparseQRType& qr,
const Derived& other,
bool transpose) :
582 m_qr(qr),m_other(other),m_transpose(transpose) {}
583 inline Index rows()
const {
return m_transpose ? m_qr.rows() : m_qr.cols(); }
584 inline Index cols()
const {
return m_other.cols(); }
587 template<
typename DesType>
588 void evalTo(DesType& res)
const 590 Index m = m_qr.rows();
591 Index n = m_qr.cols();
592 Index diagSize = (std::min)(m,n);
596 eigen_assert(m_qr.m_Q.rows() == m_other.rows() &&
"Non conforming object sizes");
598 for(Index j = 0; j < res.cols(); j++){
599 for (Index k = 0; k < diagSize; k++)
601 Scalar tau = Scalar(0);
602 tau = m_qr.m_Q.col(k).dot(res.col(j));
603 if(tau==Scalar(0))
continue;
604 tau = tau * m_qr.m_hcoeffs(k);
605 res.col(j) -= tau * m_qr.m_Q.col(k);
611 eigen_assert(m_qr.m_Q.rows() == m_other.rows() &&
"Non conforming object sizes");
613 for(Index j = 0; j < res.cols(); j++)
615 for (Index k = diagSize-1; k >=0; k--)
617 Scalar tau = Scalar(0);
618 tau = m_qr.m_Q.col(k).dot(res.col(j));
619 if(tau==Scalar(0))
continue;
620 tau = tau * m_qr.m_hcoeffs(k);
621 res.col(j) -= tau * m_qr.m_Q.col(k);
627 const SparseQRType& m_qr;
628 const Derived& m_other;
632 template<
typename SparseQRType>
633 struct SparseQRMatrixQReturnType :
public EigenBase<SparseQRMatrixQReturnType<SparseQRType> >
635 typedef typename SparseQRType::Index Index;
636 typedef typename SparseQRType::Scalar Scalar;
638 SparseQRMatrixQReturnType(
const SparseQRType& qr) : m_qr(qr) {}
639 template<
typename Derived>
642 return SparseQR_QProduct<SparseQRType,Derived>(m_qr,other.derived(),
false);
644 SparseQRMatrixQTransposeReturnType<SparseQRType> adjoint()
const 646 return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr);
648 inline Index rows()
const {
return m_qr.rows(); }
649 inline Index cols()
const {
return (std::min)(m_qr.rows(),m_qr.cols()); }
651 SparseQRMatrixQTransposeReturnType<SparseQRType> transpose()
const 653 return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr);
657 dest.derived() = m_qr.matrixQ() * Dest::Identity(m_qr.rows(), m_qr.rows());
661 Dest idMat(m_qr.rows(), m_qr.rows());
664 const_cast<SparseQRType *
>(&m_qr)->sort_matrix_Q();
665 dest.
derived() = SparseQR_QProduct<SparseQRType, Dest>(m_qr, idMat,
false);
668 const SparseQRType& m_qr;
671 template<
typename SparseQRType>
672 struct SparseQRMatrixQTransposeReturnType
674 SparseQRMatrixQTransposeReturnType(
const SparseQRType& qr) : m_qr(qr) {}
675 template<
typename Derived>
678 return SparseQR_QProduct<SparseQRType,Derived>(m_qr,other.derived(),
true);
680 const SparseQRType& m_qr;
void resize(Index newSize)
Definition: DenseBase.h:217
Index rank() const
Definition: SparseQR.h:122
const QRMatrixType & matrixR() const
Definition: SparseQR.h:116
void compute(const MatrixType &mat)
Definition: SparseQR.h:98
Index cols() const
Definition: SparseQR.h:112
const PermutationType & colsPermutation() const
Definition: SparseQR.h:152
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:88
void factorize(const MatrixType &mat)
Performs the numerical QR factorization of the input matrix.
Definition: SparseQR.h:319
Index rows() const
Definition: SparseMatrixBase.h:150
Definition: EigenBase.h:26
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
Base class of any sparse matrices or sparse expressions.
Definition: SparseMatrixBase.h:26
SparseQRMatrixQReturnType< SparseQR > matrixQ() const
Definition: SparseQR.h:146
Derived & derived()
Definition: EigenBase.h:34
SparseQR(const MatrixType &mat)
Definition: SparseQR.h:87
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:196
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SparseQR.h:229
Definition: Eigen_Colamd.h:54
Derived & setIdentity()
Definition: CwiseNullaryOp.h:772
void analyzePattern(const MatrixType &mat)
Preprocessing step of a QR factorization.
Definition: SparseQR.h:281
Definition: Constants.h:376
std::string lastErrorMessage() const
Definition: SparseQR.h:161
ComputationInfo
Definition: Constants.h:374
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
Index rows() const
Definition: SparseQR.h:108
const internal::solve_retval< SparseQR, Rhs > solve(const MatrixBase< Rhs > &B) const
Definition: SparseQR.h:207
Derived & setZero(Index size)
Definition: CwiseNullaryOp.h:515