Go to the documentation of this file.
11 #ifndef EIGEN_SPARSE_QR_H
12 #define EIGEN_SPARSE_QR_H
16 template<
typename MatrixType,
typename OrderingType>
class SparseQR;
83 template<
typename _MatrixType,
typename _OrderingType>
204 template<
typename Rhs,
typename Dest>
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;
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());
245 template<
typename Rhs>
249 eigen_assert(this->
rows() == B.rows() &&
"SparseQR::solve() : invalid number of rows in the right hand side matrix");
252 template<
typename Rhs>
256 eigen_assert(this->
rows() == B.
rows() &&
"SparseQR::solve() : invalid number of rows in the right hand side matrix");
319 template <
typename MatrixType,
typename OrderingType>
322 eigen_assert(
mat.isCompressed() &&
"SparseQR requires a sparse matrix in compressed mode. Call .makeCompressed() before passing it to SparseQR");
327 ord(matCpy, m_perm_c);
332 if (!m_perm_c.size())
339 m_outputPerm_c = m_perm_c.inverse();
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;
360 template <
typename MatrixType,
typename OrderingType>
365 eigen_assert(m_analysisIsok &&
"analyzePattern() should be called before this step");
371 Index nzcolR, nzcolQ;
380 m_outputPerm_c = m_perm_c.inverse();
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)
415 for (
int j = 0;
j <
n;
j++) max2Norm =
numext::maxi(max2Norm, m_pmat.col(
j).norm());
422 m_pivotperm.setIdentity(
n);
432 mark(nonzeroCol) =
col;
433 Qidx(0) = nonzeroCol;
434 nzcolR = 0; nzcolQ = 1;
435 bool found_diag = nonzeroCol>=
m;
446 if(curIdx == nonzeroCol) found_diag =
true;
452 m_lastError =
"Empty row found during numerical factorization";
459 for (; mark(st) !=
col; st = m_etree(st))
467 Index nt = nzcolR-bi;
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--)
492 tdot = m_Q.col(curIdx).dot(tval);
494 tdot *= m_hcoeffs(curIdx);
499 tval(itq.row()) -= itq.value() * tdot;
502 if(m_etree(Ridx(
i)) == nonzeroCol)
519 if(nonzeroCol < diagSize)
527 for (
Index itq = 1; itq < nzcolQ; ++itq) sqrNorm +=
numext::abs2(tval(Qidx(itq)));
540 for (
Index itq = 1; itq < nzcolQ; ++itq)
541 tval(Qidx(itq)) /= (c0 -
beta);
548 for (
Index i = nzcolR-1;
i >= 0;
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);
571 if(nonzeroCol<diagSize)
572 m_Q.startVec(nonzeroCol);
578 std::swap(m_pivotperm.indices()(
j), m_pivotperm.indices()[
j+1]);
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;
612 template <
typename SparseQRType,
typename Derived>
624 template<
typename DesType>
636 for (
Index k = 0; k < diagSize; k++)
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);
656 for (
Index k = start_k; k >=0; k--)
659 tau =
m_qr.m_Q.col(k).dot(
res.col(
j));
660 if(tau==
Scalar(0))
continue;
662 res.col(
j) -= tau *
m_qr.m_Q.col(k);
673 template<
typename SparseQRType>
683 template<
typename Derived>
704 template<
typename SparseQRType>
708 template<
typename Derived>
718 template<
typename SparseQRType>
726 template<
typename DstXprType,
typename SparseQRType>
734 typename DstXprType::PlainObject idMat(src.
rows(), src.
cols());
737 const_cast<SparseQRType *
>(&src.
m_qr)->_sort_matrix_Q();
742 template<
typename DstXprType,
typename SparseQRType>
750 dst = src.
m_qr.matrixQ() * DstXprType::Identity(src.
m_qr.rows(), src.
m_qr.rows());
void setPivotThreshold(const RealScalar &threshold)
SparseQRType::MatrixType MatrixType
Matrix< Scalar, Dynamic, 1 > ScalarVector
Namespace containing all symbols from the Eigen library.
SparseQRMatrixQTransposeReturnType< SparseQRType > adjoint() const
PermutationMatrix< Dynamic, Dynamic, StorageIndex > PermutationType
MatrixType::StorageIndex StorageIndex
bool _solve_impl(const MatrixBase< Rhs > &B, MatrixBase< Dest > &dest) const
void _solve_impl(const SparseMatrixBase< Rhs > &b, SparseMatrixBase< Dest > &dest) const
void analyzePattern(const MatrixType &mat)
Preprocessing step of a QR factorization.
SparseQRMatrixQReturnType< SparseQR > matrixQ() const
SparseQRType::MatrixType ReturnType
SparseQRMatrixQReturnType< SparseQRType > SrcXprType
Derived::PlainObject ReturnType
SparseQR_QProduct< SparseQRType, Derived > operator*(const MatrixBase< Derived > &other)
SparseQRType::MatrixType ReturnType
const SparseQRType & m_qr
SparseQR(const MatrixType &mat)
ReturnType::StorageKind StorageKind
const PermutationType & colsPermutation() const
Matrix< StorageIndex, Dynamic, 1 > IndexVector
SparseQR_QProduct(const SparseQRType &qr, const Derived &other, bool transpose)
void evalTo(DesType &res) const
std::string lastErrorMessage() const
cout<< "Here is the matrix m:"<< endl<< m<< endl;Matrix< ptrdiff_t, 3, 1 > res
storage_kind_to_evaluator_kind< typename MatrixType::StorageKind >::Kind Kind
double beta(double a, double b)
void factorize(const MatrixType &mat)
Performs the numerical QR factorization of the input matrix.
PermutationType m_outputPerm_c
EIGEN_DEVICE_FUNC Derived & setConstant(Index size, const Scalar &val)
SparseQRMatrixQTransposeReturnType< SparseQRType > transpose() const
static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op< Scalar, Scalar > &)
SparseQRType::QRMatrixType MatrixType
SparseMatrix< Scalar, ColMajor, StorageIndex > QRMatrixType
SparseQRType::Scalar Scalar
IndexVector m_firstRowElt
const Solve< SparseQR, Rhs > solve(const MatrixBase< Rhs > &B) const
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T mini(const T &x, const T &y)
DstXprType::StorageIndex StorageIndex
void swap(GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &a, GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &b)
const SparseQRType & m_qr
const EIGEN_DEVICE_FUNC ImagReturnType imag() const
AnnoyingScalar conj(const AnnoyingScalar &x)
_OrderingType OrderingType
const Solve< SparseQR, Rhs > solve(const SparseMatrixBase< Rhs > &B) const
PermutationType m_pivotperm
const SparseQRType & m_qr
int coletree(const MatrixType &mat, IndexVector &parent, IndexVector &firstRowElt, typename MatrixType::StorageIndex *perm=0)
static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op< Scalar, Scalar > &)
NumTraits< Scalar >::Real RealScalar
SparseQRMatrixQReturnType< SparseQRType > SrcXprType
SparseQRType::Scalar Scalar
EIGEN_DEVICE_FUNC Index size() const
SparseQR_QProduct< SparseQRType, Derived > operator*(const MatrixBase< Derived > &other)
A base class for sparse solvers.
EIGEN_DEVICE_FUNC bool abs2(bool x)
Pseudo expression representing a solving operation.
HouseholderQR< MatrixXf > qr(A)
Matrix< Scalar, Dynamic, Dynamic > DenseMatrix
Base::InnerIterator InnerIterator
MatrixType::RealScalar RealScalar
DstXprType::Scalar Scalar
ReturnType::StorageIndex StorageIndex
Base class of any sparse matrices or sparse expressions.
EIGEN_DEVICE_FUNC const EIGEN_STRONG_INLINE Scalar * data() const
SparseSolverBase< SparseQR< _MatrixType, _OrderingType > > Base
const QRMatrixType & matrixR() const
ComputationInfo info() const
Reports whether previous computation was successful.
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T maxi(const T &x, const T &y)
DstXprType::StorageIndex StorageIndex
EIGEN_DEVICE_FUNC Derived & setZero(Index size)
Base class for all dense matrices, vectors, and expressions.
void compute(const MatrixType &mat)
SparseQRMatrixQReturnType(const SparseQRType &qr)
bool m_useDefaultThreshold
DstXprType::Scalar Scalar
Sparse left-looking QR factorization with numerical column pivoting.
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Jet< T, N > sqrt(const Jet< T, N > &f)
MatrixType::Scalar Scalar
SparseQRMatrixQTransposeReturnType(const SparseQRType &qr)
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
gtsam
Author(s):
autogenerated on Fri Nov 1 2024 03:36:03