11 #ifndef EIGEN_SPARSE_QR_H 12 #define EIGEN_SPARSE_QR_H 16 template<
typename MatrixType,
typename OrderingType>
class SparseQR;
24 typedef typename ReturnType::Index
Index;
63 template<
typename _MatrixType,
typename _OrderingType>
69 typedef typename MatrixType::Scalar
Scalar;
71 typedef typename MatrixType::Index
Index;
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)
89 void analyzePattern(
const MatrixType& mat);
90 void factorize(
const MatrixType& mat);
94 inline Index
rows()
const {
return m_pmat.rows(); }
98 inline Index
cols()
const {
return m_pmat.cols();}
102 const QRMatrixType&
matrixR()
const {
return m_R; }
110 eigen_assert(m_isInitialized &&
"The factorization should be called first, use compute()");
111 return m_nonzeropivots;
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");
156 Index rank = this->rank();
159 typename Dest::PlainObject
y, b;
160 y = this->matrixQ().transpose() * B;
164 y.topRows(rank) = this->matrixR().topLeftCorner(rank, rank).template triangularView<Upper>().solve(b.topRows(rank));
165 y.bottomRows(y.size()-rank).setZero();
168 if (m_perm_c.size()) dest.
topRows(cols()) = colsPermutation() * y.topRows(cols());
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");
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");
216 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
223 if(this->m_isQSorted)
return;
227 this->m_isQSorted =
true;
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)
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;
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();
438 for (
Index itq = 1; itq < nzcolQ; ++itq) sqrNorm +=
numext::abs2(tval(Qidx(itq)));
452 for (
Index itq = 1; itq < nzcolQ; ++itq)
453 tval(Qidx(itq)) /= (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);
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;
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>
528 template<typename Dest>
void evalTo(Dest& dst)
const 530 dec()._solve(
rhs(),dst);
533 template<
typename _MatrixType,
typename OrderingType,
typename Rhs>
540 template<typename Dest>
void evalTo(Dest& dst)
const 542 this->defaultEvalTo(dst);
547 template <
typename SparseQRType,
typename Derived>
551 typedef typename SparseQRType::Scalar
Scalar;
552 typedef typename SparseQRType::Index
Index;
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>
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);
601 template<
typename SparseQRType>
604 typedef typename SparseQRType::Index
Index;
605 typedef typename SparseQRType::Scalar
Scalar;
608 template<
typename Derived>
617 inline Index
rows()
const {
return m_qr.rows(); }
618 inline Index
cols()
const {
return m_qr.cols(); }
626 dest.derived() = m_qr.matrixQ() * Dest::Identity(m_qr.rows(), m_qr.rows());
630 Dest idMat(m_qr.rows(), m_qr.rows());
633 const_cast<SparseQRType *
>(&m_qr)->sort_matrix_Q();
640 template<
typename SparseQRType>
644 template<
typename Derived>
_OrderingType OrderingType
Matrix< Scalar, Dynamic, 1 > ScalarVector
const AutoDiffScalar< DerType > & conj(const AutoDiffScalar< DerType > &x)
IntermediateState sqrt(const Expression &arg)
Derived::PlainObject ReturnType
const SparseQRType & m_qr
RowsBlockXpr topRows(Index n)
SparseQRType::MatrixType ReturnType
SparseQRMatrixQReturnType(const SparseQRType &qr)
const internal::sparse_solve_retval< SparseQR, Rhs > solve(const SparseMatrixBase< Rhs > &B) const
void evalTo(MatrixBase< Dest > &dest) const
const SparseQRType & m_qr
SparseQR_QProduct(const SparseQRType &qr, const Derived &other, bool transpose)
bool _solve(const MatrixBase< Rhs > &B, MatrixBase< Dest > &dest) const
iterative scaling algorithm to equilibrate rows and column norms in matrices
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
SparseMatrix< Scalar, ColMajor, Index > QRMatrixType
SparseQR< _MatrixType, OrderingType > Dec
SparseQRType::Scalar Scalar
EIGEN_STRONG_INLINE const CwiseUnaryOp< internal::scalar_abs2_op< Scalar >, const Derived > abs2() const
ComputationInfo info() const
Reports whether previous computation was successful.
const ImagReturnType imag() const
const QRMatrixType & matrixR() const
SparseQRType::QRMatrixType MatrixType
void factorize(const MatrixType &mat)
Performs the numerical QR factorization of the input matrix.
SparseQR_QProduct< SparseQRType, Derived > operator*(const MatrixBase< Derived > &other)
int coletree(const MatrixType &mat, IndexVector &parent, IndexVector &firstRowElt, typename MatrixType::Index *perm=0)
EIGEN_STRONG_INLINE const CwiseUnaryOp< internal::scalar_abs_op< Scalar >, const Derived > abs() const
SparseQRType::Index Index
RealReturnType real() const
Sparse left-looking rank-revealing QR factorization.
IndexVector m_firstRowElt
Base class of any sparse matrices or sparse expressions.
PermutationMatrix< Dynamic, Dynamic, Index > PermutationType
Matrix< Scalar, Dynamic, Dynamic > DenseMatrix
PermutationType m_pivotperm
SparseQRMatrixQTransposeReturnType< SparseQRType > adjoint() const
SparseQRType::Index Index
Derived & setZero(Index size)
#define EIGEN_MAKE_SPARSE_SOLVE_HELPERS(DecompositionType, Rhs)
MatrixType::Scalar Scalar
Derived & setConstant(Index size, const Scalar &value)
void evalTo(DesType &res) const
void rhs(const real_t *x, real_t *f)
void evalTo(SparseMatrixBase< Dest > &dest) const
void analyzePattern(const MatrixType &mat)
Preprocessing step of a QR factorization.
SparseQRType::Scalar Scalar
PermutationType m_outputPerm_c
SparseQRType::MatrixType ReturnType
const Derived & derived() const
const PermutationType & colsPermutation() const
const internal::solve_retval< SparseQR, Rhs > solve(const MatrixBase< Rhs > &B) const
SparseQRMatrixQReturnType< SparseQR > matrixQ() const
SparseQR< _MatrixType, OrderingType > Dec
const SparseQRType & m_qr
SparseQRMatrixQTransposeReturnType(const SparseQRType &qr)
#define EIGEN_MAKE_SOLVE_HELPERS(DecompositionType, Rhs)
Matrix< Index, Dynamic, 1 > IndexVector
MatrixType::RealScalar RealScalar
ReturnType::StorageKind StorageKind
SparseQR(const MatrixType &mat)
Base class for all dense matrices, vectors, and expressions.
SparseQRMatrixQTransposeReturnType< SparseQRType > transpose() const
void setPivotThreshold(const RealScalar &threshold)
std::string lastErrorMessage() const
bool m_useDefaultThreshold
SparseQR_QProduct< SparseQRType, Derived > operator*(const MatrixBase< Derived > &other)
void compute(const MatrixType &mat)