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;
80 : m_info(
Success), m_shiftOffset(0), m_shiftScale(1)
84 : m_info(
Success), m_shiftOffset(0), m_shiftScale(1)
86 derived().compute(matrix);
93 Derived& derived() {
return *
static_cast<Derived*
>(
this); }
94 const Derived& derived()
const {
return *
static_cast<const Derived*
>(
this); }
96 inline Index cols()
const {
return m_matrix.
cols(); }
97 inline Index rows()
const {
return m_matrix.
rows(); }
106 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
129 Derived&
setShift(
const RealScalar& offset,
const RealScalar& scale = 1)
131 m_shiftOffset = offset;
132 m_shiftScale = scale;
136 #ifndef EIGEN_PARSED_BY_DOXYGEN
138 template<
typename Stream>
139 void dumpMemory(Stream& s)
142 s <<
" L: " << ((total+=(m_matrix.
cols()+1) *
sizeof(
int) + m_matrix.nonZeros()*(
sizeof(int)+
sizeof(Scalar))) >> 20) <<
"Mb" <<
"\n";
143 s <<
" diag: " << ((total+=m_diag.size() *
sizeof(Scalar)) >> 20) <<
"Mb" <<
"\n";
144 s <<
" tree: " << ((total+=m_parent.size() *
sizeof(int)) >> 20) <<
"Mb" <<
"\n";
145 s <<
" nonzeros: " << ((total+=m_nonZerosPerCol.size() *
sizeof(int)) >> 20) <<
"Mb" <<
"\n";
146 s <<
" perm: " << ((total+=m_P.
size() *
sizeof(int)) >> 20) <<
"Mb" <<
"\n";
147 s <<
" perm^-1: " << ((total+=m_Pinv.
size() *
sizeof(int)) >> 20) <<
"Mb" <<
"\n";
148 s <<
" TOTAL: " << (total>> 20) <<
"Mb" <<
"\n";
152 template<
typename Rhs,
typename Dest>
153 void _solve_impl(
const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest)
const
155 eigen_assert(m_factorizationIsOk &&
"The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
156 eigen_assert(m_matrix.
rows()==b.rows());
166 if(m_matrix.nonZeros()>0)
167 derived().matrixL().solveInPlace(dest);
170 dest = m_diag.asDiagonal().inverse() * dest;
172 if (m_matrix.nonZeros()>0)
173 derived().matrixU().solveInPlace(dest);
176 dest = m_Pinv * dest;
179 template<
typename Rhs,
typename Dest>
180 void _solve_impl(
const SparseMatrixBase<Rhs> &b, SparseMatrixBase<Dest> &dest)
const
182 internal::solve_sparse_through_dense_panels(derived(), b, dest);
185 #endif // EIGEN_PARSED_BY_DOXYGEN
190 template<
bool DoLDLT>
193 eigen_assert(matrix.rows()==matrix.cols());
194 Index size = matrix.cols();
195 CholMatrixType tmp(size,size);
196 ConstCholMatrixPtr pmat;
197 ordering(matrix, pmat, tmp);
198 analyzePattern_preordered(*pmat, DoLDLT);
199 factorize_preordered<DoLDLT>(*pmat);
202 template<
bool DoLDLT>
203 void factorize(
const MatrixType& a)
205 eigen_assert(a.rows()==a.cols());
206 Index size = a.cols();
207 CholMatrixType tmp(size,size);
208 ConstCholMatrixPtr pmat;
213 internal::simplicial_cholesky_grab_input<CholMatrixType,MatrixType>::run(a, pmat, tmp);
217 tmp.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
221 factorize_preordered<DoLDLT>(*pmat);
224 template<
bool DoLDLT>
225 void factorize_preordered(
const CholMatrixType& a);
227 void analyzePattern(
const MatrixType& a,
bool doLDLT)
229 eigen_assert(a.rows()==a.cols());
230 Index size = a.cols();
231 CholMatrixType tmp(size,size);
232 ConstCholMatrixPtr pmat;
233 ordering(a, pmat, tmp);
234 analyzePattern_preordered(*pmat,doLDLT);
236 void analyzePattern_preordered(
const CholMatrixType& a,
bool doLDLT);
238 void ordering(
const MatrixType& a, ConstCholMatrixPtr &pmat, CholMatrixType& ap);
242 inline bool operator() (
const Index& row,
const Index& col,
const Scalar&)
const
249 bool m_factorizationIsOk;
252 CholMatrixType m_matrix;
255 VectorI m_nonZerosPerCol;
259 RealScalar m_shiftOffset;
260 RealScalar m_shiftScale;
263 template<
typename _MatrixType,
int _UpLo = Lower,
typename _Ordering = AMDOrdering<
typename _MatrixType::StorageIndex> >
class SimplicialLLT;
264 template<
typename _MatrixType,
int _UpLo = Lower,
typename _Ordering = AMDOrdering<
typename _MatrixType::StorageIndex> >
class SimplicialLDLT;
265 template<
typename _MatrixType,
int _UpLo = Lower,
typename _Ordering = AMDOrdering<
typename _MatrixType::StorageIndex> >
class SimplicialCholesky;
269 template<
typename _MatrixType,
int _UpLo,
typename _Ordering>
struct traits<
SimplicialLLT<_MatrixType,_UpLo,_Ordering> >
271 typedef _MatrixType MatrixType;
272 typedef _Ordering OrderingType;
273 enum { UpLo = _UpLo };
274 typedef typename MatrixType::Scalar Scalar;
275 typedef typename MatrixType::StorageIndex StorageIndex;
279 static inline MatrixL getL(
const MatrixType& m) {
return MatrixL(m); }
280 static inline MatrixU getU(
const MatrixType& m) {
return MatrixU(m.adjoint()); }
283 template<
typename _MatrixType,
int _UpLo,
typename _Ordering>
struct traits<SimplicialLDLT<_MatrixType,_UpLo,_Ordering> >
285 typedef _MatrixType MatrixType;
286 typedef _Ordering OrderingType;
287 enum { UpLo = _UpLo };
288 typedef typename MatrixType::Scalar Scalar;
289 typedef typename MatrixType::StorageIndex StorageIndex;
290 typedef SparseMatrix<Scalar, ColMajor, StorageIndex> CholMatrixType;
291 typedef TriangularView<const CholMatrixType, Eigen::UnitLower> MatrixL;
292 typedef TriangularView<const typename CholMatrixType::AdjointReturnType, Eigen::UnitUpper> MatrixU;
293 static inline MatrixL getL(
const MatrixType& m) {
return MatrixL(m); }
294 static inline MatrixU getU(
const MatrixType& m) {
return MatrixU(m.adjoint()); }
297 template<
typename _MatrixType,
int _UpLo,
typename _Ordering>
struct traits<SimplicialCholesky<_MatrixType,_UpLo,_Ordering> >
299 typedef _MatrixType MatrixType;
300 typedef _Ordering OrderingType;
301 enum { UpLo = _UpLo };
324 template<
typename _MatrixType,
int _UpLo,
typename _Ordering>
325 class SimplicialLLT :
public SimplicialCholeskyBase<SimplicialLLT<_MatrixType,_UpLo,_Ordering> >
328 typedef _MatrixType MatrixType;
329 enum { UpLo = _UpLo };
330 typedef SimplicialCholeskyBase<SimplicialLLT> Base;
331 typedef typename MatrixType::Scalar Scalar;
332 typedef typename MatrixType::RealScalar RealScalar;
333 typedef typename MatrixType::StorageIndex StorageIndex;
334 typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
335 typedef Matrix<Scalar,Dynamic,1> VectorType;
336 typedef internal::traits<SimplicialLLT> Traits;
337 typedef typename Traits::MatrixL MatrixL;
338 typedef typename Traits::MatrixU MatrixU;
348 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial LLT not factorized");
349 return Traits::getL(Base::m_matrix);
354 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial LLT not factorized");
355 return Traits::getU(Base::m_matrix);
361 Base::template compute<false>(matrix);
373 Base::analyzePattern(a,
false);
384 Base::template factorize<false>(a);
390 Scalar detL = Base::m_matrix.diagonal().prod();
391 return numext::abs2(detL);
413 template<
typename _MatrixType,
int _UpLo,
typename _Ordering>
414 class SimplicialLDLT :
public SimplicialCholeskyBase<SimplicialLDLT<_MatrixType,_UpLo,_Ordering> >
417 typedef _MatrixType MatrixType;
418 enum { UpLo = _UpLo };
419 typedef SimplicialCholeskyBase<SimplicialLDLT> Base;
420 typedef typename MatrixType::Scalar Scalar;
421 typedef typename MatrixType::RealScalar RealScalar;
422 typedef typename MatrixType::StorageIndex StorageIndex;
423 typedef SparseMatrix<Scalar,ColMajor,StorageIndex> CholMatrixType;
424 typedef Matrix<Scalar,Dynamic,1> VectorType;
425 typedef internal::traits<SimplicialLDLT> Traits;
426 typedef typename Traits::MatrixL MatrixL;
427 typedef typename Traits::MatrixU MatrixU;
438 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial LDLT not factorized");
443 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial LDLT not factorized");
444 return Traits::getL(Base::m_matrix);
449 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial LDLT not factorized");
450 return Traits::getU(Base::m_matrix);
456 Base::template compute<true>(matrix);
468 Base::analyzePattern(a,
true);
479 Base::template factorize<true>(a);
485 return Base::m_diag.prod();
495 template<
typename _MatrixType,
int _UpLo,
typename _Ordering>
496 class SimplicialCholesky :
public SimplicialCholeskyBase<SimplicialCholesky<_MatrixType,_UpLo,_Ordering> >
499 typedef _MatrixType MatrixType;
500 enum { UpLo = _UpLo };
501 typedef SimplicialCholeskyBase<SimplicialCholesky> Base;
502 typedef typename MatrixType::Scalar Scalar;
503 typedef typename MatrixType::RealScalar RealScalar;
504 typedef typename MatrixType::StorageIndex StorageIndex;
505 typedef SparseMatrix<Scalar,ColMajor,StorageIndex> CholMatrixType;
506 typedef Matrix<Scalar,Dynamic,1> VectorType;
507 typedef internal::traits<SimplicialCholesky> Traits;
508 typedef internal::traits<SimplicialLDLT<MatrixType,UpLo> > LDLTTraits;
509 typedef internal::traits<SimplicialLLT<MatrixType,UpLo> > LLTTraits;
511 SimplicialCholesky() : Base(), m_LDLT(true) {}
513 explicit SimplicialCholesky(
const MatrixType& matrix)
514 : Base(), m_LDLT(true)
519 SimplicialCholesky& setMode(SimplicialCholeskyMode mode)
523 case SimplicialCholeskyLLT:
526 case SimplicialCholeskyLDLT:
536 inline const VectorType vectorD()
const {
537 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial Cholesky not factorized");
540 inline const CholMatrixType rawMatrix()
const {
541 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial Cholesky not factorized");
542 return Base::m_matrix;
549 Base::template compute<true>(matrix);
551 Base::template compute<false>(matrix);
563 Base::analyzePattern(a, m_LDLT);
575 Base::template factorize<true>(a);
577 Base::template factorize<false>(a);
581 template<
typename Rhs,
typename Dest>
584 eigen_assert(Base::m_factorizationIsOk &&
"The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
585 eigen_assert(Base::m_matrix.rows()==b.rows());
590 if(Base::m_P.size()>0)
591 dest = Base::m_P * b;
595 if(Base::m_matrix.nonZeros()>0)
598 LDLTTraits::getL(Base::m_matrix).solveInPlace(dest);
600 LLTTraits::getL(Base::m_matrix).solveInPlace(dest);
603 if(Base::m_diag.size()>0)
604 dest = Base::m_diag.
asDiagonal().inverse() * dest;
606 if (Base::m_matrix.nonZeros()>0)
609 LDLTTraits::getU(Base::m_matrix).solveInPlace(dest);
611 LLTTraits::getU(Base::m_matrix).solveInPlace(dest);
614 if(Base::m_P.size()>0)
615 dest = Base::m_Pinv * dest;
619 template<
typename Rhs,
typename Dest>
620 void _solve_impl(
const SparseMatrixBase<Rhs> &b, SparseMatrixBase<Dest> &dest)
const
622 internal::solve_sparse_through_dense_panels(*
this, b, dest);
625 Scalar determinant()
const
629 return Base::m_diag.prod();
633 Scalar detL = Diagonal<const CholMatrixType>(Base::m_matrix).prod();
634 return numext::abs2(detL);
642 template<
typename Derived>
643 void SimplicialCholeskyBase<Derived>::ordering(
const MatrixType& a, ConstCholMatrixPtr &pmat, CholMatrixType& ap)
645 eigen_assert(a.rows()==a.cols());
646 const Index size = a.rows();
649 if(!internal::is_same<OrderingType,NaturalOrdering<Index> >::value)
653 C = a.template selfadjointView<UpLo>();
655 OrderingType ordering;
659 if(m_Pinv.size()>0) m_P = m_Pinv.inverse();
662 ap.resize(size,size);
663 ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
669 if(
int(UpLo)==
int(
Lower) || MatrixType::IsRowMajor)
672 ap.resize(size,size);
673 ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>();
676 internal::simplicial_cholesky_grab_input<CholMatrixType,MatrixType>::run(a, pmat, ap);
682 #endif // EIGEN_SIMPLICIAL_CHOLESKY_H
SimplicialLDLT()
Definition: SimplicialCholesky.h:430
const MatrixL matrixL() const
Definition: SimplicialCholesky.h:347
Index cols() const
Definition: SparseMatrix.h:132
SimplicialCholeskyBase()
Definition: SimplicialCholesky.h:79
Derived & setShift(const RealScalar &offset, const RealScalar &scale=1)
Definition: SimplicialCholesky.h:129
SimplicialLLT & compute(const MatrixType &matrix)
Definition: SimplicialCholesky.h:359
Definition: Constants.h:196
A base class for sparse solvers.
Definition: SparseSolverBase.h:53
Index size() const
Definition: PermutationMatrix.h:110
A direct sparse LDLT Cholesky factorizations without square root.
Definition: SimplicialCholesky.h:264
SimplicialLLT()
Definition: SimplicialCholesky.h:341
SimplicialLDLT & compute(const MatrixType &matrix)
Definition: SimplicialCholesky.h:454
void factorize(const MatrixType &a)
Definition: SimplicialCholesky.h:382
void analyzePattern(const MatrixType &a)
Definition: SimplicialCholesky.h:561
Definition: Constants.h:198
void compute(const MatrixType &matrix)
Definition: SimplicialCholesky.h:191
Scalar determinant() const
Definition: SimplicialCholesky.h:483
const PermutationMatrix< Dynamic, Dynamic, StorageIndex > & permutationP() const
Definition: SimplicialCholesky.h:112
Definition: SimplicialCholesky.h:265
A direct sparse Cholesky factorizations.
Definition: SimplicialCholesky.h:57
SimplicialLLT(const MatrixType &matrix)
Definition: SimplicialCholesky.h:343
const MatrixL matrixL() const
Definition: SimplicialCholesky.h:442
const MatrixU matrixU() const
Definition: SimplicialCholesky.h:353
void factorize(const MatrixType &a)
Definition: SimplicialCholesky.h:572
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SimplicialCholesky.h:104
const MatrixU matrixU() const
Definition: SimplicialCholesky.h:448
Definition: Constants.h:424
const PermutationMatrix< Dynamic, Dynamic, StorageIndex > & permutationPinv() const
Definition: SimplicialCholesky.h:117
Definition: Eigen_Colamd.h:54
SimplicialLDLT(const MatrixType &matrix)
Definition: SimplicialCholesky.h:433
Scalar determinant() const
Definition: SimplicialCholesky.h:388
void factorize(const MatrixType &a)
Definition: SimplicialCholesky.h:477
Definition: SimplicialCholesky.h:241
Expression of a triangular part in a matrix.
Definition: TriangularMatrix.h:186
void analyzePattern(const MatrixType &a)
Definition: SimplicialCholesky.h:371
A direct sparse LLT Cholesky factorizations.
Definition: SimplicialCholesky.h:263
const DiagonalWrapper< const Derived > asDiagonal() const
Definition: DiagonalMatrix.h:278
ComputationInfo
Definition: Constants.h:422
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
const VectorType vectorD() const
Definition: SimplicialCholesky.h:437
Index rows() const
Definition: SparseMatrix.h:130
SimplicialCholesky & compute(const MatrixType &matrix)
Definition: SimplicialCholesky.h:546
void analyzePattern(const MatrixType &a)
Definition: SimplicialCholesky.h:466