11 #ifndef EIGEN_SPARSE_QR_H
12 #define EIGEN_SPARSE_QR_H
14 #include "./InternalHeaderCheck.h"
18 template<
typename MatrixType,
typename OrderingType>
class SparseQR;
19 template<
typename SparseQRType>
struct SparseQRMatrixQReturnType;
20 template<
typename SparseQRType>
struct SparseQRMatrixQTransposeReturnType;
21 template<
typename SparseQRType,
typename Derived>
struct SparseQR_QProduct;
23 template <
typename SparseQRType>
struct traits<SparseQRMatrixQReturnType<SparseQRType> >
25 typedef typename SparseQRType::MatrixType ReturnType;
26 typedef typename ReturnType::StorageIndex StorageIndex;
27 typedef typename ReturnType::StorageKind StorageKind;
33 template <
typename SparseQRType>
struct traits<SparseQRMatrixQTransposeReturnType<SparseQRType> >
35 typedef typename SparseQRType::MatrixType ReturnType;
37 template <
typename SparseQRType,
typename Derived>
struct traits<SparseQR_QProduct<SparseQRType, Derived> >
39 typedef typename Derived::PlainObject ReturnType;
85 template<
typename MatrixType_,
typename OrderingType_>
90 using Base::m_isInitialized;
92 using Base::_solve_impl;
93 typedef MatrixType_ MatrixType;
94 typedef OrderingType_ OrderingType;
95 typedef typename MatrixType::Scalar Scalar;
96 typedef typename MatrixType::RealScalar RealScalar;
97 typedef typename MatrixType::StorageIndex StorageIndex;
104 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
105 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
109 SparseQR () : m_analysisIsok(
false), m_lastError(
""), m_useDefaultThreshold(
true),m_isQSorted(
false),m_isEtreeOk(
false)
118 explicit SparseQR(
const MatrixType& mat) : m_analysisIsok(false), m_lastError(
""), m_useDefaultThreshold(true),m_isQSorted(false),m_isEtreeOk(false)
166 eigen_assert(m_isInitialized &&
"The factorization should be called first, use compute()");
167 return m_nonzeropivots;
188 SparseQRMatrixQReturnType<SparseQR>
matrixQ()
const
189 {
return SparseQRMatrixQReturnType<SparseQR>(*
this); }
196 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
197 return m_outputPerm_c;
206 template<
typename Rhs,
typename Dest>
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");
215 typename Dest::PlainObject y, b;
216 y = this->
matrixQ().adjoint() * B;
220 y.resize((std::max<Index>)(
cols(),y.rows()),y.cols());
221 y.topRows(
rank) = this->
matrixR().topLeftCorner(rank,
rank).template triangularView<Upper>().solve(b.topRows(
rank));
222 y.bottomRows(y.rows()-
rank).setZero();
226 else dest = y.topRows(
cols());
239 m_useDefaultThreshold =
false;
240 m_threshold = threshold;
247 template<
typename Rhs>
250 eigen_assert(m_isInitialized &&
"The factorization should be called first, use compute()");
251 eigen_assert(this->
rows() == B.
rows() &&
"SparseQR::solve() : invalid number of rows in the right hand side matrix");
254 template<
typename Rhs>
257 eigen_assert(m_isInitialized &&
"The factorization should be called first, use compute()");
258 eigen_assert(this->
rows() == B.
rows() &&
"SparseQR::solve() : invalid number of rows in the right hand side matrix");
272 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
278 inline void _sort_matrix_Q()
280 if(this->m_isQSorted)
return;
284 this->m_isQSorted =
true;
290 bool m_factorizationIsok;
292 std::string m_lastError;
296 ScalarVector m_hcoeffs;
297 PermutationType m_perm_c;
298 PermutationType m_pivotperm;
299 PermutationType m_outputPerm_c;
300 RealScalar m_threshold;
301 bool m_useDefaultThreshold;
302 Index m_nonzeropivots;
304 IndexVector m_firstRowElt;
308 template <
typename,
typename >
friend struct SparseQR_QProduct;
321 template <
typename MatrixType,
typename OrderingType>
324 eigen_assert(mat.isCompressed() &&
"SparseQR requires a sparse matrix in compressed mode. Call .makeCompressed() before passing it to SparseQR");
326 std::conditional_t<MatrixType::IsRowMajor,QRMatrixType,const MatrixType&> matCpy(mat);
329 ord(matCpy, m_perm_c);
330 Index n = mat.cols();
331 Index m = mat.rows();
332 Index diagSize = (std::min)(m,n);
334 if (!m_perm_c.size())
337 m_perm_c.indices().setLinSpaced(n, 0,StorageIndex(n-1));
341 m_outputPerm_c = m_perm_c.inverse();
342 internal::coletree(matCpy, m_etree, m_firstRowElt, m_outputPerm_c.indices().data());
346 m_Q.resize(m, diagSize);
349 m_R.reserve(2*mat.nonZeros());
350 m_Q.reserve(2*mat.nonZeros());
351 m_hcoeffs.resize(diagSize);
352 m_analysisIsok =
true;
362 template <
typename MatrixType,
typename OrderingType>
367 eigen_assert(m_analysisIsok &&
"analyzePattern() should be called before this step");
368 StorageIndex m = StorageIndex(mat.rows());
369 StorageIndex n = StorageIndex(mat.cols());
370 StorageIndex diagSize = (std::min)(m,n);
373 Index nzcolR, nzcolQ;
375 RealScalar pivotThreshold = m_threshold;
382 m_outputPerm_c = m_perm_c.inverse();
383 internal::coletree(m_pmat, m_etree, m_firstRowElt, m_outputPerm_c.indices().data());
395 const StorageIndex *originalOuterIndices = mat.outerIndexPtr();
396 if(MatrixType::IsRowMajor)
398 originalOuterIndicesCpy = IndexVector::Map(m_pmat.outerIndexPtr(),n+1);
399 originalOuterIndices = originalOuterIndicesCpy.
data();
402 for (
int i = 0; i < n; i++)
404 Index p = m_perm_c.size() ? m_perm_c.indices()(i) : i;
405 m_pmat.outerIndexPtr()[p] = originalOuterIndices[i];
406 m_pmat.innerNonZeroPtr()[p] = originalOuterIndices[i+1] - originalOuterIndices[i];
414 if(m_useDefaultThreshold)
416 RealScalar max2Norm = 0.0;
417 for (
int j = 0; j < n; j++) max2Norm = numext::maxi(max2Norm, m_pmat.col(j).norm());
418 if(max2Norm==RealScalar(0))
419 max2Norm = RealScalar(1);
424 m_pivotperm.setIdentity(n);
426 StorageIndex nonzeroCol = 0;
430 for (StorageIndex col = 0; col < n; ++col)
434 mark(nonzeroCol) = col;
435 Qidx(0) = nonzeroCol;
436 nzcolR = 0; nzcolQ = 1;
437 bool found_diag = nonzeroCol>=m;
444 for (
typename QRMatrixType::InnerIterator itp(m_pmat, col); itp || !found_diag; ++itp)
446 StorageIndex curIdx = nonzeroCol;
447 if(itp) curIdx = StorageIndex(itp.row());
448 if(curIdx == nonzeroCol) found_diag =
true;
451 StorageIndex st = m_firstRowElt(curIdx);
454 m_lastError =
"Empty row found during numerical factorization";
461 for (; mark(st) != col; st = m_etree(st))
469 Index nt = nzcolR-bi;
470 for(
Index i = 0; i < nt/2; i++) std::swap(Ridx(bi+i), Ridx(nzcolR-i-1));
473 if(itp) tval(curIdx) = itp.value();
474 else tval(curIdx) = Scalar(0);
477 if(curIdx > nonzeroCol && mark(curIdx) != col )
479 Qidx(nzcolQ) = curIdx;
486 for (
Index i = nzcolR-1; i >= 0; i--)
488 Index curIdx = Ridx(i);
494 tdot = m_Q.col(curIdx).dot(tval);
496 tdot *= m_hcoeffs(curIdx);
500 for (
typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq)
501 tval(itq.row()) -= itq.value() * tdot;
504 if(m_etree(Ridx(i)) == nonzeroCol)
506 for (
typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq)
508 StorageIndex iQ = StorageIndex(itq.row());
518 Scalar tau = RealScalar(0);
521 if(nonzeroCol < diagSize)
525 Scalar c0 = nzcolQ ? tval(Qidx(0)) : Scalar(0);
528 RealScalar sqrNorm = 0.;
529 for (
Index itq = 1; itq < nzcolQ; ++itq) sqrNorm += numext::abs2(tval(Qidx(itq)));
530 if(sqrNorm == RealScalar(0) && numext::imag(c0) == RealScalar(0))
532 beta = numext::real(c0);
538 beta =
sqrt(numext::abs2(c0) + sqrNorm);
539 if(numext::real(c0) >= RealScalar(0))
542 for (
Index itq = 1; itq < nzcolQ; ++itq)
543 tval(Qidx(itq)) /= (c0 - beta);
544 tau = numext::conj((beta-c0) / beta);
550 for (
Index i = nzcolR-1; i >= 0; i--)
552 Index curIdx = Ridx(i);
553 if(curIdx < nonzeroCol)
555 m_R.insertBackByOuterInnerUnordered(col, curIdx) = tval(curIdx);
556 tval(curIdx) = Scalar(0.);
560 if(nonzeroCol < diagSize &&
abs(beta) >= pivotThreshold)
562 m_R.insertBackByOuterInner(col, nonzeroCol) = beta;
564 m_hcoeffs(nonzeroCol) = tau;
566 for (
Index itq = 0; itq < nzcolQ; ++itq)
568 Index iQ = Qidx(itq);
569 m_Q.insertBackByOuterInnerUnordered(nonzeroCol,iQ) = tval(iQ);
570 tval(iQ) = Scalar(0.);
573 if(nonzeroCol<diagSize)
574 m_Q.startVec(nonzeroCol);
579 for (
Index j = nonzeroCol; j < n-1; j++)
580 std::swap(m_pivotperm.indices()(j), m_pivotperm.indices()[j+1]);
583 internal::coletree(m_pmat, m_etree, m_firstRowElt, m_pivotperm.indices().data());
588 m_hcoeffs.tail(diagSize-nonzeroCol).setZero();
592 m_Q.makeCompressed();
594 m_R.makeCompressed();
597 m_nonzeropivots = nonzeroCol;
603 m_R = tempR * m_pivotperm;
606 m_outputPerm_c = m_outputPerm_c * m_pivotperm;
609 m_isInitialized =
true;
610 m_factorizationIsok =
true;
614 template <
typename SparseQRType,
typename Derived>
615 struct SparseQR_QProduct : ReturnByValue<SparseQR_QProduct<SparseQRType, Derived> >
617 typedef typename SparseQRType::QRMatrixType MatrixType;
618 typedef typename SparseQRType::Scalar Scalar;
620 SparseQR_QProduct(
const SparseQRType& qr,
const Derived& other,
bool transpose) :
621 m_qr(qr),m_other(other),m_transpose(transpose) {}
622 inline Index rows()
const {
return m_qr.matrixQ().rows(); }
623 inline Index cols()
const {
return m_other.cols(); }
626 template<
typename DesType>
627 void evalTo(DesType& res)
const
629 Index m = m_qr.rows();
630 Index n = m_qr.cols();
631 Index diagSize = (std::min)(m,n);
635 eigen_assert(m_qr.m_Q.rows() == m_other.rows() &&
"Non conforming object sizes");
637 for(
Index j = 0; j < res.cols(); j++){
638 for (
Index k = 0; k < diagSize; k++)
640 Scalar tau = Scalar(0);
641 tau = m_qr.m_Q.col(k).dot(res.col(j));
642 if(tau==Scalar(0))
continue;
643 tau = tau * m_qr.m_hcoeffs(k);
644 res.col(j) -= tau * m_qr.m_Q.col(k);
650 eigen_assert(m_qr.matrixQ().cols() == m_other.rows() &&
"Non conforming object sizes");
652 res.conservativeResize(rows(), cols());
655 for(
Index j = 0; j < res.cols(); j++)
657 Index start_k = internal::is_identity<Derived>::value ? numext::mini(j,diagSize-1) : diagSize-1;
658 for (
Index k = start_k; k >=0; k--)
660 Scalar tau = Scalar(0);
661 tau = m_qr.m_Q.col(k).dot(res.col(j));
662 if(tau==Scalar(0))
continue;
663 tau = tau * numext::conj(m_qr.m_hcoeffs(k));
664 res.col(j) -= tau * m_qr.m_Q.col(k);
670 const SparseQRType& m_qr;
671 const Derived& m_other;
675 template<
typename SparseQRType>
676 struct SparseQRMatrixQReturnType :
public EigenBase<SparseQRMatrixQReturnType<SparseQRType> >
678 typedef typename SparseQRType::Scalar Scalar;
679 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
684 explicit SparseQRMatrixQReturnType(
const SparseQRType& qr) : m_qr(qr) {}
685 template<
typename Derived>
686 SparseQR_QProduct<SparseQRType, Derived>
operator*(
const MatrixBase<Derived>& other)
688 return SparseQR_QProduct<SparseQRType,Derived>(m_qr,other.derived(),
false);
691 SparseQRMatrixQTransposeReturnType<SparseQRType> adjoint()
const
693 return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr);
695 inline Index rows()
const {
return m_qr.rows(); }
696 inline Index cols()
const {
return m_qr.rows(); }
698 SparseQRMatrixQTransposeReturnType<SparseQRType> transpose()
const
700 return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr);
702 const SparseQRType& m_qr;
706 template<
typename SparseQRType>
707 struct SparseQRMatrixQTransposeReturnType
709 explicit SparseQRMatrixQTransposeReturnType(
const SparseQRType& qr) : m_qr(qr) {}
710 template<
typename Derived>
711 SparseQR_QProduct<SparseQRType,Derived>
operator*(
const MatrixBase<Derived>& other)
713 return SparseQR_QProduct<SparseQRType,Derived>(m_qr,other.derived(),
true);
715 const SparseQRType& m_qr;
720 template<
typename SparseQRType>
721 struct evaluator_traits<SparseQRMatrixQReturnType<SparseQRType> >
723 typedef typename SparseQRType::MatrixType MatrixType;
724 typedef typename storage_kind_to_evaluator_kind<typename MatrixType::StorageKind>::Kind Kind;
725 typedef SparseShape Shape;
728 template<
typename DstXprType,
typename SparseQRType>
729 struct Assignment<DstXprType, SparseQRMatrixQReturnType<SparseQRType>, internal::assign_op<typename DstXprType::Scalar,typename DstXprType::Scalar>, Sparse2Sparse>
731 typedef SparseQRMatrixQReturnType<SparseQRType> SrcXprType;
732 typedef typename DstXprType::Scalar Scalar;
733 typedef typename DstXprType::StorageIndex StorageIndex;
734 static void run(DstXprType &dst,
const SrcXprType &src,
const internal::assign_op<Scalar,Scalar> &)
736 typename DstXprType::PlainObject idMat(src.rows(), src.cols());
739 const_cast<SparseQRType *
>(&src.m_qr)->_sort_matrix_Q();
740 dst = SparseQR_QProduct<SparseQRType, DstXprType>(src.m_qr, idMat,
false);
744 template<
typename DstXprType,
typename SparseQRType>
745 struct Assignment<DstXprType, SparseQRMatrixQReturnType<SparseQRType>, internal::assign_op<typename DstXprType::Scalar,typename DstXprType::Scalar>, Sparse2Dense>
747 typedef SparseQRMatrixQReturnType<SparseQRType> SrcXprType;
748 typedef typename DstXprType::Scalar Scalar;
749 typedef typename DstXprType::StorageIndex StorageIndex;
750 static void run(DstXprType &dst,
const SrcXprType &src,
const internal::assign_op<Scalar,Scalar> &)
752 dst = src.m_qr.matrixQ() * DstXprType::Identity(src.m_qr.rows(), src.m_qr.rows());
Derived & derived()
Definition: EigenBase.h:48
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: EigenBase.h:62
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:52
Index size() const
Definition: PermutationMatrix.h:99
Derived & setConstant(Index size, const Scalar &val)
Definition: CwiseNullaryOp.h:363
const Scalar * data() const
Definition: PlainObjectBase.h:259
Derived & setZero(Index size)
Definition: CwiseNullaryOp.h:564
Pseudo expression representing a solving operation.
Definition: Solve.h:65
Base class of any sparse matrices or sparse expressions.
Definition: SparseMatrixBase.h:30
Index rows() const
Definition: SparseMatrixBase.h:176
Index cols() const
Definition: SparseMatrix.h:142
Index rows() const
Definition: SparseMatrix.h:140
Sparse left-looking QR factorization with numerical column pivoting.
Definition: SparseQR.h:87
void setPivotThreshold(const RealScalar &threshold)
Definition: SparseQR.h:237
void analyzePattern(const MatrixType &mat)
Preprocessing step of a QR factorization.
Definition: SparseQR.h:322
void factorize(const MatrixType &mat)
Performs the numerical QR factorization of the input matrix.
Definition: SparseQR.h:363
SparseQR(const MatrixType &mat)
Definition: SparseQR.h:118
Index cols() const
Definition: SparseQR.h:143
const QRMatrixType & matrixR() const
Definition: SparseQR.h:158
Index rows() const
Definition: SparseQR.h:139
const PermutationType & colsPermutation() const
Definition: SparseQR.h:194
SparseQRMatrixQReturnType< SparseQR > matrixQ() const
Definition: SparseQR.h:188
const Solve< SparseQR, Rhs > solve(const MatrixBase< Rhs > &B) const
Definition: SparseQR.h:248
std::string lastErrorMessage() const
Definition: SparseQR.h:203
void compute(const MatrixType &mat)
Definition: SparseQR.h:129
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SparseQR.h:270
Index rank() const
Definition: SparseQR.h:164
A base class for sparse solvers.
Definition: SparseSolverBase.h:70
ComputationInfo
Definition: Constants.h:442
@ InvalidInput
Definition: Constants.h:451
@ Success
Definition: Constants.h:444
Namespace containing all symbols from the Eigen library.
Definition: Core:139
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
const Product< MatrixDerived, PermutationDerived, AliasFreeProduct > operator*(const MatrixBase< MatrixDerived > &matrix, const PermutationBase< PermutationDerived > &permutation)
Definition: PermutationMatrix.h:517
const int Dynamic
Definition: Constants.h:24
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_sqrt_op< typename Derived::Scalar >, const Derived > sqrt(const Eigen::ArrayBase< Derived > &x)
Derived & derived()
Definition: EigenBase.h:48
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:231