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 &)
54 template<
typename Derived>
58 using Base::m_isInitialized;
61 typedef typename internal::traits<Derived>::MatrixType MatrixType;
62 typedef typename internal::traits<Derived>::OrderingType OrderingType;
63 enum { UpLo = internal::traits<Derived>::UpLo };
64 typedef typename MatrixType::Scalar Scalar;
65 typedef typename MatrixType::RealScalar RealScalar;
66 typedef typename MatrixType::StorageIndex StorageIndex;
73 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
74 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
83 : m_info(
Success), m_shiftOffset(0), m_shiftScale(1)
87 : m_info(
Success), m_shiftOffset(0), m_shiftScale(1)
89 derived().compute(matrix);
96 Derived& derived() {
return *static_cast<Derived*>(
this); }
97 const Derived& derived()
const {
return *static_cast<const Derived*>(
this); }
99 inline Index cols()
const {
return m_matrix.
cols(); }
100 inline Index rows()
const {
return m_matrix.
rows(); }
109 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
132 Derived&
setShift(
const RealScalar& offset,
const RealScalar& scale = 1)
134 m_shiftOffset = offset;
135 m_shiftScale = scale;
139 #ifndef EIGEN_PARSED_BY_DOXYGEN 141 template<
typename Stream>
142 void dumpMemory(Stream& s)
145 s <<
" L: " << ((total+=(m_matrix.
cols()+1) *
sizeof(
int) + m_matrix.
nonZeros()*(
sizeof(int)+
sizeof(Scalar))) >> 20) <<
"Mb" <<
"\n";
146 s <<
" diag: " << ((total+=m_diag.size() *
sizeof(Scalar)) >> 20) <<
"Mb" <<
"\n";
147 s <<
" tree: " << ((total+=m_parent.size() *
sizeof(int)) >> 20) <<
"Mb" <<
"\n";
148 s <<
" nonzeros: " << ((total+=m_nonZerosPerCol.size() *
sizeof(int)) >> 20) <<
"Mb" <<
"\n";
149 s <<
" perm: " << ((total+=m_P.
size() *
sizeof(int)) >> 20) <<
"Mb" <<
"\n";
150 s <<
" perm^-1: " << ((total+=m_Pinv.
size() *
sizeof(int)) >> 20) <<
"Mb" <<
"\n";
151 s <<
" TOTAL: " << (total>> 20) <<
"Mb" <<
"\n";
155 template<
typename Rhs,
typename Dest>
156 void _solve_impl(
const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest)
const 158 eigen_assert(m_factorizationIsOk &&
"The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
159 eigen_assert(m_matrix.
rows()==b.rows());
170 derived().matrixL().solveInPlace(dest);
173 dest = m_diag.asDiagonal().inverse() * dest;
176 derived().matrixU().solveInPlace(dest);
179 dest = m_Pinv * dest;
182 template<
typename Rhs,
typename Dest>
183 void _solve_impl(
const SparseMatrixBase<Rhs> &b, SparseMatrixBase<Dest> &dest)
const 185 internal::solve_sparse_through_dense_panels(derived(), b, dest);
188 #endif // EIGEN_PARSED_BY_DOXYGEN 193 template<
bool DoLDLT>
196 eigen_assert(matrix.rows()==matrix.cols());
197 Index size = matrix.cols();
199 ConstCholMatrixPtr pmat;
200 ordering(matrix, pmat, tmp);
201 analyzePattern_preordered(*pmat, DoLDLT);
202 factorize_preordered<DoLDLT>(*pmat);
205 template<
bool DoLDLT>
206 void factorize(
const MatrixType& a)
208 eigen_assert(a.rows()==a.cols());
209 Index size = a.cols();
210 CholMatrixType tmp(size,size);
211 ConstCholMatrixPtr pmat;
216 internal::simplicial_cholesky_grab_input<CholMatrixType,MatrixType>::run(a, pmat, tmp);
220 tmp.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
224 factorize_preordered<DoLDLT>(*pmat);
227 template<
bool DoLDLT>
228 void factorize_preordered(
const CholMatrixType& a);
230 void analyzePattern(
const MatrixType& a,
bool doLDLT)
232 eigen_assert(a.rows()==a.cols());
233 Index size = a.cols();
234 CholMatrixType tmp(size,size);
235 ConstCholMatrixPtr pmat;
236 ordering(a, pmat, tmp);
237 analyzePattern_preordered(*pmat,doLDLT);
239 void analyzePattern_preordered(
const CholMatrixType& a,
bool doLDLT);
241 void ordering(
const MatrixType& a, ConstCholMatrixPtr &pmat, CholMatrixType& ap);
245 inline bool operator() (
const Index& row,
const Index& col,
const Scalar&)
const 252 bool m_factorizationIsOk;
262 RealScalar m_shiftOffset;
263 RealScalar m_shiftScale;
266 template<
typename _MatrixType,
int _UpLo = Lower,
typename _Ordering = AMDOrdering<
typename _MatrixType::StorageIndex> >
class SimplicialLLT;
267 template<
typename _MatrixType,
int _UpLo = Lower,
typename _Ordering = AMDOrdering<
typename _MatrixType::StorageIndex> >
class SimplicialLDLT;
268 template<
typename _MatrixType,
int _UpLo = Lower,
typename _Ordering = AMDOrdering<
typename _MatrixType::StorageIndex> >
class SimplicialCholesky;
272 template<
typename _MatrixType,
int _UpLo,
typename _Ordering>
struct traits<
SimplicialLLT<_MatrixType,_UpLo,_Ordering> >
274 typedef _MatrixType MatrixType;
275 typedef _Ordering OrderingType;
276 enum { UpLo = _UpLo };
277 typedef typename MatrixType::Scalar Scalar;
278 typedef typename MatrixType::StorageIndex StorageIndex;
282 static inline MatrixL getL(
const MatrixType& m) {
return MatrixL(m); }
283 static inline MatrixU getU(
const MatrixType& m) {
return MatrixU(m.adjoint()); }
286 template<
typename _MatrixType,
int _UpLo,
typename _Ordering>
struct traits<SimplicialLDLT<_MatrixType,_UpLo,_Ordering> >
288 typedef _MatrixType MatrixType;
289 typedef _Ordering OrderingType;
290 enum { UpLo = _UpLo };
291 typedef typename MatrixType::Scalar Scalar;
292 typedef typename MatrixType::StorageIndex StorageIndex;
293 typedef SparseMatrix<Scalar, ColMajor, StorageIndex> CholMatrixType;
294 typedef TriangularView<const CholMatrixType, Eigen::UnitLower> MatrixL;
295 typedef TriangularView<const typename CholMatrixType::AdjointReturnType, Eigen::UnitUpper> MatrixU;
296 static inline MatrixL getL(
const MatrixType& m) {
return MatrixL(m); }
297 static inline MatrixU getU(
const MatrixType& m) {
return MatrixU(m.adjoint()); }
300 template<
typename _MatrixType,
int _UpLo,
typename _Ordering>
struct traits<SimplicialCholesky<_MatrixType,_UpLo,_Ordering> >
302 typedef _MatrixType MatrixType;
303 typedef _Ordering OrderingType;
304 enum { UpLo = _UpLo };
329 template<
typename _MatrixType,
int _UpLo,
typename _Ordering>
330 class SimplicialLLT :
public SimplicialCholeskyBase<SimplicialLLT<_MatrixType,_UpLo,_Ordering> >
333 typedef _MatrixType MatrixType;
334 enum { UpLo = _UpLo };
335 typedef SimplicialCholeskyBase<SimplicialLLT> Base;
336 typedef typename MatrixType::Scalar Scalar;
337 typedef typename MatrixType::RealScalar RealScalar;
338 typedef typename MatrixType::StorageIndex StorageIndex;
339 typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
340 typedef Matrix<Scalar,Dynamic,1> VectorType;
341 typedef internal::traits<SimplicialLLT> Traits;
342 typedef typename Traits::MatrixL MatrixL;
343 typedef typename Traits::MatrixU MatrixU;
353 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial LLT not factorized");
354 return Traits::getL(Base::m_matrix);
359 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial LLT not factorized");
360 return Traits::getU(Base::m_matrix);
366 Base::template compute<false>(matrix);
378 Base::analyzePattern(a,
false);
389 Base::template factorize<false>(a);
395 Scalar detL = Base::m_matrix.diagonal().prod();
396 return numext::abs2(detL);
420 template<
typename _MatrixType,
int _UpLo,
typename _Ordering>
421 class SimplicialLDLT :
public SimplicialCholeskyBase<SimplicialLDLT<_MatrixType,_UpLo,_Ordering> >
424 typedef _MatrixType MatrixType;
425 enum { UpLo = _UpLo };
426 typedef SimplicialCholeskyBase<SimplicialLDLT> Base;
427 typedef typename MatrixType::Scalar Scalar;
428 typedef typename MatrixType::RealScalar RealScalar;
429 typedef typename MatrixType::StorageIndex StorageIndex;
430 typedef SparseMatrix<Scalar,ColMajor,StorageIndex> CholMatrixType;
431 typedef Matrix<Scalar,Dynamic,1> VectorType;
432 typedef internal::traits<SimplicialLDLT> Traits;
433 typedef typename Traits::MatrixL MatrixL;
434 typedef typename Traits::MatrixU MatrixU;
445 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial LDLT not factorized");
450 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial LDLT not factorized");
451 return Traits::getL(Base::m_matrix);
456 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial LDLT not factorized");
457 return Traits::getU(Base::m_matrix);
463 Base::template compute<true>(matrix);
475 Base::analyzePattern(a,
true);
486 Base::template factorize<true>(a);
492 return Base::m_diag.prod();
502 template<
typename _MatrixType,
int _UpLo,
typename _Ordering>
503 class SimplicialCholesky :
public SimplicialCholeskyBase<SimplicialCholesky<_MatrixType,_UpLo,_Ordering> >
506 typedef _MatrixType MatrixType;
507 enum { UpLo = _UpLo };
508 typedef SimplicialCholeskyBase<SimplicialCholesky> Base;
509 typedef typename MatrixType::Scalar Scalar;
510 typedef typename MatrixType::RealScalar RealScalar;
511 typedef typename MatrixType::StorageIndex StorageIndex;
512 typedef SparseMatrix<Scalar,ColMajor,StorageIndex> CholMatrixType;
513 typedef Matrix<Scalar,Dynamic,1> VectorType;
514 typedef internal::traits<SimplicialCholesky> Traits;
515 typedef internal::traits<SimplicialLDLT<MatrixType,UpLo> > LDLTTraits;
516 typedef internal::traits<SimplicialLLT<MatrixType,UpLo> > LLTTraits;
518 SimplicialCholesky() : Base(), m_LDLT(true) {}
520 explicit SimplicialCholesky(
const MatrixType& matrix)
521 : Base(), m_LDLT(true)
526 SimplicialCholesky& setMode(SimplicialCholeskyMode mode)
530 case SimplicialCholeskyLLT:
533 case SimplicialCholeskyLDLT:
543 inline const VectorType vectorD()
const {
544 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial Cholesky not factorized");
547 inline const CholMatrixType rawMatrix()
const {
548 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial Cholesky not factorized");
549 return Base::m_matrix;
556 Base::template compute<true>(matrix);
558 Base::template compute<false>(matrix);
570 Base::analyzePattern(a, m_LDLT);
582 Base::template factorize<true>(a);
584 Base::template factorize<false>(a);
588 template<
typename Rhs,
typename Dest>
591 eigen_assert(Base::m_factorizationIsOk &&
"The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
592 eigen_assert(Base::m_matrix.rows()==b.
rows());
597 if(Base::m_P.size()>0)
598 dest = Base::m_P * b;
602 if(Base::m_matrix.nonZeros()>0)
605 LDLTTraits::getL(Base::m_matrix).solveInPlace(dest);
607 LLTTraits::getL(Base::m_matrix).solveInPlace(dest);
610 if(Base::m_diag.size()>0)
611 dest = Base::m_diag.
asDiagonal().inverse() * dest;
613 if (Base::m_matrix.nonZeros()>0)
616 LDLTTraits::getU(Base::m_matrix).solveInPlace(dest);
618 LLTTraits::getU(Base::m_matrix).solveInPlace(dest);
621 if(Base::m_P.size()>0)
622 dest = Base::m_Pinv * dest;
626 template<
typename Rhs,
typename Dest>
627 void _solve_impl(
const SparseMatrixBase<Rhs> &b, SparseMatrixBase<Dest> &dest)
const 629 internal::solve_sparse_through_dense_panels(*
this, b, dest);
632 Scalar determinant()
const 636 return Base::m_diag.prod();
640 Scalar detL = Diagonal<const CholMatrixType>(Base::m_matrix).prod();
641 return numext::abs2(detL);
649 template<
typename Derived>
650 void SimplicialCholeskyBase<Derived>::ordering(
const MatrixType& a, ConstCholMatrixPtr &pmat, CholMatrixType& ap)
652 eigen_assert(a.rows()==a.cols());
653 const Index size = a.rows();
656 if(!internal::is_same<OrderingType,NaturalOrdering<Index> >::value)
660 C = a.template selfadjointView<UpLo>();
662 OrderingType ordering;
666 if(m_Pinv.size()>0) m_P = m_Pinv.inverse();
669 ap.resize(size,size);
670 ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
676 if(
int(UpLo)==
int(
Lower) || MatrixType::IsRowMajor)
679 ap.resize(size,size);
680 ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>();
683 internal::simplicial_cholesky_grab_input<CholMatrixType,MatrixType>::run(a, pmat, ap);
689 #endif // EIGEN_SIMPLICIAL_CHOLESKY_H SimplicialLDLT()
Definition: SimplicialCholesky.h:437
Index cols() const
Definition: SparseMatrix.h:138
Scalar determinant() const
Definition: SimplicialCholesky.h:490
const MatrixU matrixU() const
Definition: SimplicialCholesky.h:358
Index nonZeros() const
Definition: SparseCompressedBase.h:56
SimplicialCholeskyBase()
Definition: SimplicialCholesky.h:82
Derived & setShift(const RealScalar &offset, const RealScalar &scale=1)
Definition: SimplicialCholesky.h:132
SimplicialLLT & compute(const MatrixType &matrix)
Definition: SimplicialCholesky.h:364
const MatrixU matrixU() const
Definition: SimplicialCholesky.h:455
Index rows() const
Definition: SparseMatrix.h:136
A base class for sparse solvers.
Definition: SparseSolverBase.h:67
Namespace containing all symbols from the Eigen library.
Definition: Core:306
A direct sparse LDLT Cholesky factorizations without square root.
Definition: SimplicialCholesky.h:267
SimplicialLLT()
Definition: SimplicialCholesky.h:346
SimplicialLDLT & compute(const MatrixType &matrix)
Definition: SimplicialCholesky.h:461
Scalar determinant() const
Definition: SimplicialCholesky.h:393
const MatrixL matrixL() const
Definition: SimplicialCholesky.h:449
void factorize(const MatrixType &a)
Definition: SimplicialCholesky.h:387
Index size() const
Definition: PermutationMatrix.h:108
void analyzePattern(const MatrixType &a)
Definition: SimplicialCholesky.h:568
void compute(const MatrixType &matrix)
Definition: SimplicialCholesky.h:194
Definition: Constants.h:204
Definition: SimplicialCholesky.h:268
A base class for direct sparse Cholesky factorizations.
Definition: SimplicialCholesky.h:55
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
SimplicialLLT(const MatrixType &matrix)
Definition: SimplicialCholesky.h:348
const DiagonalWrapper< const Derived > asDiagonal() const
Definition: DiagonalMatrix.h:277
void factorize(const MatrixType &a)
Definition: SimplicialCholesky.h:579
Definition: Constants.h:206
Definition: Constants.h:432
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SimplicialCholesky.h:107
const PermutationMatrix< Dynamic, Dynamic, StorageIndex > & permutationPinv() const
Definition: SimplicialCholesky.h:120
const PermutationMatrix< Dynamic, Dynamic, StorageIndex > & permutationP() const
Definition: SimplicialCholesky.h:115
SimplicialLDLT(const MatrixType &matrix)
Definition: SimplicialCholesky.h:440
void factorize(const MatrixType &a)
Definition: SimplicialCholesky.h:484
Definition: SimplicialCholesky.h:244
const VectorType vectorD() const
Definition: SimplicialCholesky.h:444
Expression of a triangular part in a matrix.
Definition: TriangularMatrix.h:186
Index rows() const
Definition: EigenBase.h:59
void analyzePattern(const MatrixType &a)
Definition: SimplicialCholesky.h:376
A direct sparse LLT Cholesky factorizations.
Definition: SimplicialCholesky.h:266
ComputationInfo
Definition: Constants.h:430
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
const MatrixL matrixL() const
Definition: SimplicialCholesky.h:352
SimplicialCholesky & compute(const MatrixType &matrix)
Definition: SimplicialCholesky.h:553
void analyzePattern(const MatrixType &a)
Definition: SimplicialCholesky.h:473