10 #ifndef EIGEN_SIMPLICIAL_CHOLESKY_H 11 #define EIGEN_SIMPLICIAL_CHOLESKY_H 15 enum SimplicialCholeskyMode {
16 SimplicialCholeskyLLT,
17 SimplicialCholeskyLDLT
21 template<
typename CholMatrixType,
typename InputMatrixType>
22 struct simplicial_cholesky_grab_input {
23 typedef CholMatrixType
const * ConstCholMatrixPtr;
24 static void run(
const InputMatrixType& input, ConstCholMatrixPtr &pmat, CholMatrixType &tmp)
31 template<
typename MatrixType>
32 struct simplicial_cholesky_grab_input<MatrixType,MatrixType> {
33 typedef MatrixType
const * ConstMatrixPtr;
34 static void run(
const MatrixType& input, ConstMatrixPtr &pmat, MatrixType &)
56 template<
typename Derived>
60 using Base::m_isInitialized;
63 typedef typename internal::traits<Derived>::MatrixType MatrixType;
64 typedef typename internal::traits<Derived>::OrderingType OrderingType;
65 enum { UpLo = internal::traits<Derived>::UpLo };
66 typedef typename MatrixType::Scalar Scalar;
67 typedef typename MatrixType::RealScalar RealScalar;
68 typedef typename MatrixType::StorageIndex StorageIndex;
70 typedef CholMatrixType
const * ConstCholMatrixPtr;
75 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
76 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
85 : m_info(
Success), m_shiftOffset(0), m_shiftScale(1)
89 : m_info(
Success), m_shiftOffset(0), m_shiftScale(1)
91 derived().compute(matrix);
98 Derived& derived() {
return *
static_cast<Derived*
>(
this); }
99 const Derived& derived()
const {
return *
static_cast<const Derived*
>(
this); }
101 inline Index cols()
const {
return m_matrix.cols(); }
102 inline Index rows()
const {
return m_matrix.rows(); }
111 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
134 Derived&
setShift(
const RealScalar& offset,
const RealScalar& scale = 1)
136 m_shiftOffset = offset;
137 m_shiftScale = scale;
141 #ifndef EIGEN_PARSED_BY_DOXYGEN 143 template<
typename Stream>
144 void dumpMemory(Stream& s)
147 s <<
" L: " << ((total+=(m_matrix.cols()+1) *
sizeof(
int) + m_matrix.nonZeros()*(
sizeof(int)+
sizeof(Scalar))) >> 20) <<
"Mb" <<
"\n";
148 s <<
" diag: " << ((total+=m_diag.size() *
sizeof(Scalar)) >> 20) <<
"Mb" <<
"\n";
149 s <<
" tree: " << ((total+=m_parent.size() *
sizeof(int)) >> 20) <<
"Mb" <<
"\n";
150 s <<
" nonzeros: " << ((total+=m_nonZerosPerCol.size() *
sizeof(int)) >> 20) <<
"Mb" <<
"\n";
151 s <<
" perm: " << ((total+=m_P.size() *
sizeof(int)) >> 20) <<
"Mb" <<
"\n";
152 s <<
" perm^-1: " << ((total+=m_Pinv.size() *
sizeof(int)) >> 20) <<
"Mb" <<
"\n";
153 s <<
" TOTAL: " << (total>> 20) <<
"Mb" <<
"\n";
157 template<
typename Rhs,
typename Dest>
160 eigen_assert(m_factorizationIsOk &&
"The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
161 eigen_assert(m_matrix.rows()==b.rows());
171 if(m_matrix.nonZeros()>0)
172 derived().matrixL().solveInPlace(dest);
177 if (m_matrix.nonZeros()>0)
178 derived().matrixU().solveInPlace(dest);
181 dest = m_Pinv * dest;
184 template<
typename Rhs,
typename Dest>
187 internal::solve_sparse_through_dense_panels(derived(), b, dest);
190 #endif // EIGEN_PARSED_BY_DOXYGEN 195 template<
bool DoLDLT>
198 eigen_assert(matrix.rows()==matrix.cols());
199 Index size = matrix.cols();
200 CholMatrixType tmp(size,size);
201 ConstCholMatrixPtr pmat;
202 ordering(matrix, pmat, tmp);
203 analyzePattern_preordered(*pmat, DoLDLT);
204 factorize_preordered<DoLDLT>(*pmat);
207 template<
bool DoLDLT>
208 void factorize(
const MatrixType& a)
210 eigen_assert(a.rows()==a.cols());
211 Index size = a.cols();
212 CholMatrixType tmp(size,size);
213 ConstCholMatrixPtr pmat;
218 internal::simplicial_cholesky_grab_input<CholMatrixType,MatrixType>::run(a, pmat, tmp);
222 tmp.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().
twistedBy(m_P);
226 factorize_preordered<DoLDLT>(*pmat);
229 template<
bool DoLDLT>
230 void factorize_preordered(
const CholMatrixType& a);
232 void analyzePattern(
const MatrixType& a,
bool doLDLT)
234 eigen_assert(a.rows()==a.cols());
235 Index size = a.cols();
236 CholMatrixType tmp(size,size);
237 ConstCholMatrixPtr pmat;
238 ordering(a, pmat, tmp);
239 analyzePattern_preordered(*pmat,doLDLT);
241 void analyzePattern_preordered(
const CholMatrixType& a,
bool doLDLT);
243 void ordering(
const MatrixType& a, ConstCholMatrixPtr &pmat, CholMatrixType& ap);
247 inline bool operator() (
const Index& row,
const Index& col,
const Scalar&)
const 254 bool m_factorizationIsOk;
257 CholMatrixType m_matrix;
260 VectorI m_nonZerosPerCol;
264 RealScalar m_shiftOffset;
265 RealScalar m_shiftScale;
268 template<
typename _MatrixType,
int _UpLo = Lower,
typename _Ordering = AMDOrdering<
typename _MatrixType::StorageIndex> >
class SimplicialLLT;
269 template<
typename _MatrixType,
int _UpLo = Lower,
typename _Ordering = AMDOrdering<
typename _MatrixType::StorageIndex> >
class SimplicialLDLT;
270 template<
typename _MatrixType,
int _UpLo = Lower,
typename _Ordering = AMDOrdering<
typename _MatrixType::StorageIndex> >
class SimplicialCholesky;
274 template<
typename _MatrixType,
int _UpLo,
typename _Ordering>
struct traits<
SimplicialLLT<_MatrixType,_UpLo,_Ordering> >
276 typedef _MatrixType MatrixType;
277 typedef _Ordering OrderingType;
278 enum { UpLo = _UpLo };
279 typedef typename MatrixType::Scalar Scalar;
280 typedef typename MatrixType::StorageIndex StorageIndex;
284 static inline MatrixL getL(
const MatrixType& m) {
return MatrixL(m); }
285 static inline MatrixU getU(
const MatrixType& m) {
return MatrixU(m.adjoint()); }
288 template<
typename _MatrixType,
int _UpLo,
typename _Ordering>
struct traits<
SimplicialLDLT<_MatrixType,_UpLo,_Ordering> >
290 typedef _MatrixType MatrixType;
291 typedef _Ordering OrderingType;
292 enum { UpLo = _UpLo };
293 typedef typename MatrixType::Scalar Scalar;
294 typedef typename MatrixType::StorageIndex StorageIndex;
298 static inline MatrixL getL(
const MatrixType& m) {
return MatrixL(m); }
299 static inline MatrixU getU(
const MatrixType& m) {
return MatrixU(m.adjoint()); }
302 template<
typename _MatrixType,
int _UpLo,
typename _Ordering>
struct traits<
SimplicialCholesky<_MatrixType,_UpLo,_Ordering> >
304 typedef _MatrixType MatrixType;
305 typedef _Ordering OrderingType;
306 enum { UpLo = _UpLo };
331 template<
typename _MatrixType,
int _UpLo,
typename _Ordering>
335 typedef _MatrixType MatrixType;
336 enum { UpLo = _UpLo };
338 typedef typename MatrixType::Scalar Scalar;
339 typedef typename MatrixType::RealScalar RealScalar;
340 typedef typename MatrixType::StorageIndex StorageIndex;
343 typedef internal::traits<SimplicialLLT> Traits;
344 typedef typename Traits::MatrixL MatrixL;
345 typedef typename Traits::MatrixU MatrixU;
355 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial LLT not factorized");
356 return Traits::getL(Base::m_matrix);
361 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial LLT not factorized");
362 return Traits::getU(Base::m_matrix);
368 Base::template compute<false>(matrix);
380 Base::analyzePattern(a,
false);
391 Base::template factorize<false>(a);
397 Scalar detL = Base::m_matrix.diagonal().prod();
398 return numext::abs2(detL);
422 template<
typename _MatrixType,
int _UpLo,
typename _Ordering>
426 typedef _MatrixType MatrixType;
427 enum { UpLo = _UpLo };
429 typedef typename MatrixType::Scalar Scalar;
430 typedef typename MatrixType::RealScalar RealScalar;
431 typedef typename MatrixType::StorageIndex StorageIndex;
434 typedef internal::traits<SimplicialLDLT> Traits;
435 typedef typename Traits::MatrixL MatrixL;
436 typedef typename Traits::MatrixU MatrixU;
447 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial LDLT not factorized");
452 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial LDLT not factorized");
453 return Traits::getL(Base::m_matrix);
458 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial LDLT not factorized");
459 return Traits::getU(Base::m_matrix);
465 Base::template compute<true>(matrix);
477 Base::analyzePattern(a,
true);
488 Base::template factorize<true>(a);
494 return Base::m_diag.prod();
504 template<
typename _MatrixType,
int _UpLo,
typename _Ordering>
508 typedef _MatrixType MatrixType;
509 enum { UpLo = _UpLo };
511 typedef typename MatrixType::Scalar Scalar;
512 typedef typename MatrixType::RealScalar RealScalar;
513 typedef typename MatrixType::StorageIndex StorageIndex;
516 typedef internal::traits<SimplicialCholesky> Traits;
517 typedef internal::traits<SimplicialLDLT<MatrixType,UpLo> > LDLTTraits;
518 typedef internal::traits<SimplicialLLT<MatrixType,UpLo> > LLTTraits;
523 : Base(), m_LDLT(
true)
532 case SimplicialCholeskyLLT:
535 case SimplicialCholeskyLDLT:
545 inline const VectorType vectorD()
const {
546 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial Cholesky not factorized");
549 inline const CholMatrixType rawMatrix()
const {
550 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial Cholesky not factorized");
551 return Base::m_matrix;
558 Base::template compute<true>(matrix);
560 Base::template compute<false>(matrix);
572 Base::analyzePattern(a, m_LDLT);
584 Base::template factorize<true>(a);
586 Base::template factorize<false>(a);
590 template<
typename Rhs,
typename Dest>
593 eigen_assert(Base::m_factorizationIsOk &&
"The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
594 eigen_assert(Base::m_matrix.rows()==b.rows());
599 if(Base::m_P.size()>0)
600 dest = Base::m_P * b;
604 if(Base::m_matrix.nonZeros()>0)
607 LDLTTraits::getL(Base::m_matrix).solveInPlace(dest);
609 LLTTraits::getL(Base::m_matrix).solveInPlace(dest);
612 if(Base::m_diag.size()>0)
613 dest = Base::m_diag.
asDiagonal().inverse() * dest;
615 if (Base::m_matrix.nonZeros()>0)
618 LDLTTraits::getU(Base::m_matrix).solveInPlace(dest);
620 LLTTraits::getU(Base::m_matrix).solveInPlace(dest);
623 if(Base::m_P.size()>0)
624 dest = Base::m_Pinv * dest;
628 template<
typename Rhs,
typename Dest>
631 internal::solve_sparse_through_dense_panels(*
this, b, dest);
634 Scalar determinant()
const 638 return Base::m_diag.prod();
643 return numext::abs2(detL);
651 template<
typename Derived>
654 eigen_assert(a.rows()==a.cols());
655 const Index size = a.rows();
662 C = a.template selfadjointView<UpLo>();
664 OrderingType ordering;
668 if(m_Pinv.size()>0) m_P = m_Pinv.inverse();
672 ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().
twistedBy(m_P);
678 if(
int(UpLo)==
int(
Lower) || MatrixType::IsRowMajor)
682 ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>();
685 internal::simplicial_cholesky_grab_input<CholMatrixType,MatrixType>::run(a, pmat, ap);
691 #endif // EIGEN_SIMPLICIAL_CHOLESKY_H SimplicialLDLT()
Definition: SimplicialCholesky.h:439
const MatrixL matrixL() const
Definition: SimplicialCholesky.h:354
SimplicialCholeskyBase()
Definition: SimplicialCholesky.h:84
Derived & setShift(const RealScalar &offset, const RealScalar &scale=1)
Definition: SimplicialCholesky.h:134
SimplicialLLT & compute(const MatrixType &matrix)
Definition: SimplicialCholesky.h:366
A base class for sparse solvers.
Definition: SparseSolverBase.h:53
Definition: Ordering.h:94
A direct sparse LDLT Cholesky factorizations without square root.
Definition: SimplicialCholesky.h:269
SimplicialLLT()
Definition: SimplicialCholesky.h:348
SimplicialLDLT & compute(const MatrixType &matrix)
Definition: SimplicialCholesky.h:463
void factorize(const MatrixType &a)
Definition: SimplicialCholesky.h:389
void analyzePattern(const MatrixType &a)
Definition: SimplicialCholesky.h:570
void compute(const MatrixType &matrix)
Definition: SimplicialCholesky.h:196
Definition: Constants.h:204
Base class of any sparse matrices or sparse expressions.
Definition: ForwardDeclarations.h:278
Scalar determinant() const
Definition: SimplicialCholesky.h:492
const PermutationMatrix< Dynamic, Dynamic, StorageIndex > & permutationP() const
Definition: SimplicialCholesky.h:117
Definition: SimplicialCholesky.h:270
A direct sparse Cholesky factorizations.
Definition: SimplicialCholesky.h:57
SimplicialLLT(const MatrixType &matrix)
Definition: SimplicialCholesky.h:350
const MatrixL matrixL() const
Definition: SimplicialCholesky.h:451
void resize(Index rows, Index cols)
Definition: SparseMatrix.h:611
const MatrixU matrixU() const
Definition: SimplicialCholesky.h:360
void factorize(const MatrixType &a)
Definition: SimplicialCholesky.h:581
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SimplicialCholesky.h:109
const MatrixU matrixU() const
Definition: SimplicialCholesky.h:457
SparseSymmetricPermutationProduct< Derived, Upper|Lower > twistedBy(const PermutationMatrix< Dynamic, Dynamic, StorageIndex > &perm) const
Definition: SparseMatrixBase.h:309
Definition: Constants.h:432
const PermutationMatrix< Dynamic, Dynamic, StorageIndex > & permutationPinv() const
Definition: SimplicialCholesky.h:122
Definition: Eigen_Colamd.h:54
SimplicialLDLT(const MatrixType &matrix)
Definition: SimplicialCholesky.h:442
Scalar determinant() const
Definition: SimplicialCholesky.h:395
void factorize(const MatrixType &a)
Definition: SimplicialCholesky.h:486
Definition: SimplicialCholesky.h:246
Expression of a triangular part in a matrix.
Definition: TriangularMatrix.h:186
void analyzePattern(const MatrixType &a)
Definition: SimplicialCholesky.h:378
Expression of a diagonal/subdiagonal/superdiagonal in a matrix.
Definition: Diagonal.h:63
A direct sparse LLT Cholesky factorizations.
Definition: SimplicialCholesky.h:268
const DiagonalWrapper< const Derived > asDiagonal() const
Definition: DiagonalMatrix.h:278
Definition: Constants.h:206
ComputationInfo
Definition: Constants.h:430
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
const VectorType vectorD() const
Definition: SimplicialCholesky.h:446
SimplicialCholesky & compute(const MatrixType &matrix)
Definition: SimplicialCholesky.h:555
void analyzePattern(const MatrixType &a)
Definition: SimplicialCholesky.h:475