11#ifndef EIGEN_SPARSE_QR_H
12#define EIGEN_SPARSE_QR_H
16template<
typename MatrixType,
typename OrderingType>
class SparseQR;
17template<
typename SparseQRType>
struct SparseQRMatrixQReturnType;
18template<
typename SparseQRType>
struct SparseQRMatrixQTransposeReturnType;
19template<
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;
83template<
typename _MatrixType,
typename _OrderingType>
88 using Base::m_isInitialized;
90 using Base::_solve_impl;
91 typedef _MatrixType MatrixType;
92 typedef _OrderingType OrderingType;
93 typedef typename MatrixType::Scalar Scalar;
94 typedef typename MatrixType::RealScalar RealScalar;
95 typedef typename MatrixType::StorageIndex StorageIndex;
102 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
103 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
107 SparseQR () : m_analysisIsok(
false), m_lastError(
""), m_useDefaultThreshold(
true),m_isQSorted(
false),m_isEtreeOk(
false)
116 explicit SparseQR(
const MatrixType& mat) : m_analysisIsok(false), m_lastError(
""), m_useDefaultThreshold(true),m_isQSorted(false),m_isEtreeOk(false)
164 eigen_assert(m_isInitialized &&
"The factorization should be called first, use compute()");
165 return m_nonzeropivots;
186 SparseQRMatrixQReturnType<SparseQR>
matrixQ()
const
187 {
return SparseQRMatrixQReturnType<SparseQR>(*
this); }
194 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
195 return m_outputPerm_c;
204 template<
typename Rhs,
typename Dest>
207 eigen_assert(m_isInitialized &&
"The factorization should be called first, use compute()");
208 eigen_assert(this->
rows() == B.
rows() &&
"SparseQR::solve() : invalid number of rows in the right hand side matrix");
213 typename Dest::PlainObject y, b;
214 y = this->
matrixQ().adjoint() * B;
218 y.resize((std::max<Index>)(
cols(),y.rows()),y.cols());
219 y.topRows(
rank) = this->
matrixR().topLeftCorner(rank,
rank).template triangularView<Upper>().solve(b.topRows(
rank));
220 y.bottomRows(y.rows()-
rank).setZero();
224 else dest = y.topRows(
cols());
237 m_useDefaultThreshold =
false;
238 m_threshold = threshold;
245 template<
typename Rhs>
248 eigen_assert(m_isInitialized &&
"The factorization should be called first, use compute()");
249 eigen_assert(this->
rows() == B.
rows() &&
"SparseQR::solve() : invalid number of rows in the right hand side matrix");
252 template<
typename Rhs>
255 eigen_assert(m_isInitialized &&
"The factorization should be called first, use compute()");
256 eigen_assert(this->
rows() == B.
rows() &&
"SparseQR::solve() : invalid number of rows in the right hand side matrix");
270 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
276 inline void _sort_matrix_Q()
278 if(this->m_isQSorted)
return;
282 this->m_isQSorted =
true;
288 bool m_factorizationIsok;
290 std::string m_lastError;
294 ScalarVector m_hcoeffs;
295 PermutationType m_perm_c;
296 PermutationType m_pivotperm;
297 PermutationType m_outputPerm_c;
298 RealScalar m_threshold;
299 bool m_useDefaultThreshold;
300 Index m_nonzeropivots;
302 IndexVector m_firstRowElt;
306 template <
typename,
typename >
friend struct SparseQR_QProduct;
319template <
typename MatrixType,
typename OrderingType>
322 eigen_assert(mat.isCompressed() &&
"SparseQR requires a sparse matrix in compressed mode. Call .makeCompressed() before passing it to SparseQR");
324 typename internal::conditional<MatrixType::IsRowMajor,QRMatrixType,const MatrixType&>::type matCpy(mat);
327 ord(matCpy, m_perm_c);
328 Index n = mat.cols();
329 Index m = mat.rows();
330 Index diagSize = (std::min)(m,n);
332 if (!m_perm_c.size())
335 m_perm_c.indices().setLinSpaced(n, 0,StorageIndex(n-1));
339 m_outputPerm_c = m_perm_c.inverse();
340 internal::coletree(matCpy, m_etree, m_firstRowElt, m_outputPerm_c.indices().data());
344 m_Q.resize(m, diagSize);
347 m_R.reserve(2*mat.nonZeros());
348 m_Q.reserve(2*mat.nonZeros());
349 m_hcoeffs.resize(diagSize);
350 m_analysisIsok =
true;
360template <
typename MatrixType,
typename OrderingType>
365 eigen_assert(m_analysisIsok &&
"analyzePattern() should be called before this step");
366 StorageIndex m = StorageIndex(mat.rows());
367 StorageIndex n = StorageIndex(mat.cols());
368 StorageIndex diagSize = (std::min)(m,n);
371 Index nzcolR, nzcolQ;
373 RealScalar pivotThreshold = m_threshold;
380 m_outputPerm_c = m_perm_c.inverse();
381 internal::coletree(m_pmat, m_etree, m_firstRowElt, m_outputPerm_c.indices().data());
393 const StorageIndex *originalOuterIndices = mat.outerIndexPtr();
394 if(MatrixType::IsRowMajor)
396 originalOuterIndicesCpy = IndexVector::Map(m_pmat.outerIndexPtr(),n+1);
397 originalOuterIndices = originalOuterIndicesCpy.
data();
400 for (
int i = 0; i < n; i++)
402 Index p = m_perm_c.size() ? m_perm_c.indices()(i) : i;
403 m_pmat.outerIndexPtr()[p] = originalOuterIndices[i];
404 m_pmat.innerNonZeroPtr()[p] = originalOuterIndices[i+1] - originalOuterIndices[i];
412 if(m_useDefaultThreshold)
414 RealScalar max2Norm = 0.0;
415 for (
int j = 0; j < n; j++) max2Norm = numext::maxi(max2Norm, m_pmat.col(j).norm());
416 if(max2Norm==RealScalar(0))
417 max2Norm = RealScalar(1);
422 m_pivotperm.setIdentity(n);
424 StorageIndex nonzeroCol = 0;
428 for (StorageIndex col = 0; col < n; ++col)
432 mark(nonzeroCol) = col;
433 Qidx(0) = nonzeroCol;
434 nzcolR = 0; nzcolQ = 1;
435 bool found_diag = nonzeroCol>=m;
442 for (
typename QRMatrixType::InnerIterator itp(m_pmat, col); itp || !found_diag; ++itp)
444 StorageIndex curIdx = nonzeroCol;
445 if(itp) curIdx = StorageIndex(itp.row());
446 if(curIdx == nonzeroCol) found_diag =
true;
449 StorageIndex st = m_firstRowElt(curIdx);
452 m_lastError =
"Empty row found during numerical factorization";
459 for (; mark(st) != col; st = m_etree(st))
467 Index nt = nzcolR-bi;
468 for(
Index i = 0; i < nt/2; i++) std::swap(Ridx(bi+i), Ridx(nzcolR-i-1));
471 if(itp) tval(curIdx) = itp.value();
472 else tval(curIdx) = Scalar(0);
475 if(curIdx > nonzeroCol && mark(curIdx) != col )
477 Qidx(nzcolQ) = curIdx;
484 for (
Index i = nzcolR-1; i >= 0; i--)
486 Index curIdx = Ridx(i);
492 tdot = m_Q.col(curIdx).dot(tval);
494 tdot *= m_hcoeffs(curIdx);
498 for (
typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq)
499 tval(itq.row()) -= itq.value() * tdot;
502 if(m_etree(Ridx(i)) == nonzeroCol)
504 for (
typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq)
506 StorageIndex iQ = StorageIndex(itq.row());
516 Scalar tau = RealScalar(0);
519 if(nonzeroCol < diagSize)
523 Scalar c0 = nzcolQ ? tval(Qidx(0)) : Scalar(0);
526 RealScalar sqrNorm = 0.;
527 for (
Index itq = 1; itq < nzcolQ; ++itq) sqrNorm += numext::abs2(tval(Qidx(itq)));
528 if(sqrNorm == RealScalar(0) && numext::imag(c0) == RealScalar(0))
530 beta = numext::real(c0);
536 beta = sqrt(numext::abs2(c0) + sqrNorm);
537 if(numext::real(c0) >= RealScalar(0))
540 for (
Index itq = 1; itq < nzcolQ; ++itq)
541 tval(Qidx(itq)) /= (c0 - beta);
542 tau = numext::conj((beta-c0) / beta);
548 for (
Index i = nzcolR-1; i >= 0; i--)
550 Index curIdx = Ridx(i);
551 if(curIdx < nonzeroCol)
553 m_R.insertBackByOuterInnerUnordered(col, curIdx) = tval(curIdx);
554 tval(curIdx) = Scalar(0.);
558 if(nonzeroCol < diagSize && abs(beta) >= pivotThreshold)
560 m_R.insertBackByOuterInner(col, nonzeroCol) = beta;
562 m_hcoeffs(nonzeroCol) = tau;
564 for (
Index itq = 0; itq < nzcolQ; ++itq)
566 Index iQ = Qidx(itq);
567 m_Q.insertBackByOuterInnerUnordered(nonzeroCol,iQ) = tval(iQ);
568 tval(iQ) = Scalar(0.);
571 if(nonzeroCol<diagSize)
572 m_Q.startVec(nonzeroCol);
577 for (
Index j = nonzeroCol; j < n-1; j++)
578 std::swap(m_pivotperm.indices()(j), m_pivotperm.indices()[j+1]);
581 internal::coletree(m_pmat, m_etree, m_firstRowElt, m_pivotperm.indices().data());
586 m_hcoeffs.tail(diagSize-nonzeroCol).setZero();
590 m_Q.makeCompressed();
592 m_R.makeCompressed();
595 m_nonzeropivots = nonzeroCol;
601 m_R = tempR * m_pivotperm;
604 m_outputPerm_c = m_outputPerm_c * m_pivotperm;
607 m_isInitialized =
true;
608 m_factorizationIsok =
true;
612template <
typename SparseQRType,
typename Derived>
613struct SparseQR_QProduct : ReturnByValue<SparseQR_QProduct<SparseQRType, Derived> >
615 typedef typename SparseQRType::QRMatrixType MatrixType;
616 typedef typename SparseQRType::Scalar Scalar;
618 SparseQR_QProduct(
const SparseQRType& qr,
const Derived& other,
bool transpose) :
619 m_qr(qr),m_other(other),m_transpose(transpose) {}
620 inline Index rows()
const {
return m_qr.matrixQ().rows(); }
621 inline Index cols()
const {
return m_other.cols(); }
624 template<
typename DesType>
625 void evalTo(DesType& res)
const
627 Index m = m_qr.rows();
628 Index n = m_qr.cols();
629 Index diagSize = (std::min)(m,n);
633 eigen_assert(m_qr.m_Q.rows() == m_other.rows() &&
"Non conforming object sizes");
635 for(
Index j = 0; j < res.cols(); j++){
636 for (
Index k = 0; k < diagSize; k++)
638 Scalar tau = Scalar(0);
639 tau = m_qr.m_Q.col(k).dot(res.col(j));
640 if(tau==Scalar(0))
continue;
641 tau = tau * m_qr.m_hcoeffs(k);
642 res.col(j) -= tau * m_qr.m_Q.col(k);
648 eigen_assert(m_qr.matrixQ().cols() == m_other.rows() &&
"Non conforming object sizes");
650 res.conservativeResize(rows(), cols());
653 for(
Index j = 0; j < res.cols(); j++)
655 Index start_k = internal::is_identity<Derived>::value ? numext::mini(j,diagSize-1) : diagSize-1;
656 for (
Index k = start_k; k >=0; k--)
658 Scalar tau = Scalar(0);
659 tau = m_qr.m_Q.col(k).dot(res.col(j));
660 if(tau==Scalar(0))
continue;
661 tau = tau * numext::conj(m_qr.m_hcoeffs(k));
662 res.col(j) -= tau * m_qr.m_Q.col(k);
668 const SparseQRType& m_qr;
669 const Derived& m_other;
673template<
typename SparseQRType>
674struct SparseQRMatrixQReturnType :
public EigenBase<SparseQRMatrixQReturnType<SparseQRType> >
676 typedef typename SparseQRType::Scalar Scalar;
677 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
682 explicit SparseQRMatrixQReturnType(
const SparseQRType& qr) : m_qr(qr) {}
683 template<
typename Derived>
684 SparseQR_QProduct<SparseQRType, Derived> operator*(
const MatrixBase<Derived>& other)
686 return SparseQR_QProduct<SparseQRType,Derived>(m_qr,other.derived(),
false);
689 SparseQRMatrixQTransposeReturnType<SparseQRType> adjoint()
const
691 return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr);
693 inline Index rows()
const {
return m_qr.rows(); }
694 inline Index cols()
const {
return m_qr.rows(); }
696 SparseQRMatrixQTransposeReturnType<SparseQRType> transpose()
const
698 return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr);
700 const SparseQRType& m_qr;
704template<
typename SparseQRType>
705struct SparseQRMatrixQTransposeReturnType
707 explicit SparseQRMatrixQTransposeReturnType(
const SparseQRType& qr) : m_qr(qr) {}
708 template<
typename Derived>
709 SparseQR_QProduct<SparseQRType,Derived> operator*(
const MatrixBase<Derived>& other)
711 return SparseQR_QProduct<SparseQRType,Derived>(m_qr,other.derived(),
true);
713 const SparseQRType& m_qr;
718template<
typename SparseQRType>
719struct evaluator_traits<SparseQRMatrixQReturnType<SparseQRType> >
721 typedef typename SparseQRType::MatrixType MatrixType;
722 typedef typename storage_kind_to_evaluator_kind<typename MatrixType::StorageKind>::Kind Kind;
723 typedef SparseShape Shape;
726template<
typename DstXprType,
typename SparseQRType>
727struct Assignment<DstXprType, SparseQRMatrixQReturnType<SparseQRType>, internal::assign_op<typename DstXprType::Scalar,typename DstXprType::Scalar>, Sparse2Sparse>
729 typedef SparseQRMatrixQReturnType<SparseQRType> SrcXprType;
730 typedef typename DstXprType::Scalar Scalar;
731 typedef typename DstXprType::StorageIndex StorageIndex;
732 static void run(DstXprType &dst,
const SrcXprType &src,
const internal::assign_op<Scalar,Scalar> &)
734 typename DstXprType::PlainObject idMat(src.rows(), src.cols());
737 const_cast<SparseQRType *
>(&src.m_qr)->_sort_matrix_Q();
738 dst = SparseQR_QProduct<SparseQRType, DstXprType>(src.m_qr, idMat,
false);
742template<
typename DstXprType,
typename SparseQRType>
743struct Assignment<DstXprType, SparseQRMatrixQReturnType<SparseQRType>, internal::assign_op<typename DstXprType::Scalar,typename DstXprType::Scalar>, Sparse2Dense>
745 typedef SparseQRMatrixQReturnType<SparseQRType> SrcXprType;
746 typedef typename DstXprType::Scalar Scalar;
747 typedef typename DstXprType::StorageIndex StorageIndex;
748 static void run(DstXprType &dst,
const SrcXprType &src,
const internal::assign_op<Scalar,Scalar> &)
750 dst = src.m_qr.matrixQ() * DstXprType::Identity(src.m_qr.rows(), src.m_qr.rows());
Derived & derived()
Definition: EigenBase.h:46
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: EigenBase.h:60
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:50
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:180
Index size() const
Definition: PermutationMatrix.h:97
Permutation matrix.
Definition: PermutationMatrix.h:298
Derived & setConstant(Index size, const Scalar &val)
Definition: CwiseNullaryOp.h:361
Derived & setZero(Index size)
Definition: CwiseNullaryOp.h:562
const Scalar * data() const
Definition: PlainObjectBase.h:247
Pseudo expression representing a solving operation.
Definition: Solve.h:63
Base class of any sparse matrices or sparse expressions.
Definition: SparseMatrixBase.h:28
Index rows() const
Definition: SparseMatrixBase.h:176
A versatible sparse matrix representation.
Definition: SparseMatrix.h:98
Index rows() const
Definition: SparseMatrix.h:138
Index cols() const
Definition: SparseMatrix.h:140
Sparse left-looking QR factorization with numerical column pivoting.
Definition: SparseQR.h:85
std::string lastErrorMessage() const
Definition: SparseQR.h:201
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SparseQR.h:268
void analyzePattern(const MatrixType &mat)
Preprocessing step of a QR factorization.
Definition: SparseQR.h:320
void factorize(const MatrixType &mat)
Performs the numerical QR factorization of the input matrix.
Definition: SparseQR.h:361
Index cols() const
Definition: SparseQR.h:141
const Solve< SparseQR, Rhs > solve(const MatrixBase< Rhs > &B) const
Definition: SparseQR.h:246
const PermutationType & colsPermutation() const
Definition: SparseQR.h:192
Index rank() const
Definition: SparseQR.h:162
SparseQRMatrixQReturnType< SparseQR > matrixQ() const
Definition: SparseQR.h:186
const QRMatrixType & matrixR() const
Definition: SparseQR.h:156
Index rows() const
Definition: SparseQR.h:137
SparseQR(const MatrixType &mat)
Definition: SparseQR.h:116
void setPivotThreshold(const RealScalar &threshold)
Definition: SparseQR.h:235
void compute(const MatrixType &mat)
Definition: SparseQR.h:127
A base class for sparse solvers.
Definition: SparseSolverBase.h:68
ComputationInfo
Definition: Constants.h:440
@ InvalidInput
Definition: Constants.h:449
@ Success
Definition: Constants.h:442
Namespace containing all symbols from the Eigen library.
Definition: Core:141
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
const int Dynamic
Definition: Constants.h:22
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:39
Derived & derived()
Definition: EigenBase.h:46
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:233