10#ifndef EIGEN_SIMPLICIAL_CHOLESKY_H
11#define EIGEN_SIMPLICIAL_CHOLESKY_H
13#include "./InternalHeaderCheck.h"
17enum SimplicialCholeskyMode {
18 SimplicialCholeskyLLT,
19 SimplicialCholeskyLDLT
23 template<
typename CholMatrixType,
typename InputMatrixType>
24 struct simplicial_cholesky_grab_input {
25 typedef CholMatrixType
const * ConstCholMatrixPtr;
26 static void run(
const InputMatrixType& input, ConstCholMatrixPtr &pmat, CholMatrixType &tmp)
33 template<
typename MatrixType>
34 struct simplicial_cholesky_grab_input<MatrixType,MatrixType> {
35 typedef MatrixType
const * ConstMatrixPtr;
36 static void run(
const MatrixType& input, ConstMatrixPtr &pmat, MatrixType &)
56template<
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;
75 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
76 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
86 m_factorizationIsOk(false),
87 m_analysisIsOk(false),
94 m_factorizationIsOk(false),
95 m_analysisIsOk(false),
99 derived().compute(matrix);
102 ~SimplicialCholeskyBase()
106 Derived& derived() {
return *
static_cast<Derived*
>(
this); }
107 const Derived& derived()
const {
return *
static_cast<const Derived*
>(
this); }
109 inline Index cols()
const {
return m_matrix.
cols(); }
110 inline Index rows()
const {
return m_matrix.
rows(); }
119 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
142 Derived&
setShift(
const RealScalar& offset,
const RealScalar& scale = 1)
144 m_shiftOffset = offset;
145 m_shiftScale = scale;
149#ifndef EIGEN_PARSED_BY_DOXYGEN
151 template<
typename Stream>
152 void dumpMemory(Stream& s)
155 s <<
" L: " << ((total+=(m_matrix.
cols()+1) *
sizeof(
int) + m_matrix.
nonZeros()*(
sizeof(int)+
sizeof(Scalar))) >> 20) <<
"Mb" <<
"\n";
156 s <<
" diag: " << ((total+=m_diag.size() *
sizeof(Scalar)) >> 20) <<
"Mb" <<
"\n";
157 s <<
" tree: " << ((total+=m_parent.size() *
sizeof(int)) >> 20) <<
"Mb" <<
"\n";
158 s <<
" nonzeros: " << ((total+=m_nonZerosPerCol.size() *
sizeof(int)) >> 20) <<
"Mb" <<
"\n";
159 s <<
" perm: " << ((total+=m_P.
size() *
sizeof(int)) >> 20) <<
"Mb" <<
"\n";
160 s <<
" perm^-1: " << ((total+=m_Pinv.
size() *
sizeof(int)) >> 20) <<
"Mb" <<
"\n";
161 s <<
" TOTAL: " << (total>> 20) <<
"Mb" <<
"\n";
165 template<
typename Rhs,
typename Dest>
166 void _solve_impl(
const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest)
const
168 eigen_assert(m_factorizationIsOk &&
"The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
169 eigen_assert(m_matrix.
rows()==b.rows());
180 derived().matrixL().solveInPlace(dest);
183 dest = m_diag.asDiagonal().inverse() * dest;
186 derived().matrixU().solveInPlace(dest);
189 dest = m_Pinv * dest;
192 template<
typename Rhs,
typename Dest>
193 void _solve_impl(
const SparseMatrixBase<Rhs> &b, SparseMatrixBase<Dest> &dest)
const
195 internal::solve_sparse_through_dense_panels(derived(), b, dest);
203 template<
bool DoLDLT>
206 eigen_assert(matrix.rows()==matrix.cols());
207 Index size = matrix.cols();
209 ConstCholMatrixPtr pmat;
210 ordering(matrix, pmat, tmp);
211 analyzePattern_preordered(*pmat, DoLDLT);
212 factorize_preordered<DoLDLT>(*pmat);
215 template<
bool DoLDLT>
216 void factorize(
const MatrixType& a)
218 eigen_assert(a.rows()==a.cols());
219 Index size = a.cols();
220 CholMatrixType tmp(size,size);
221 ConstCholMatrixPtr pmat;
226 internal::simplicial_cholesky_grab_input<CholMatrixType,MatrixType>::run(a, pmat, tmp);
230 tmp.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
234 factorize_preordered<DoLDLT>(*pmat);
237 template<
bool DoLDLT>
238 void factorize_preordered(
const CholMatrixType& a);
240 void analyzePattern(
const MatrixType& a,
bool doLDLT)
242 eigen_assert(a.rows()==a.cols());
243 Index size = a.cols();
244 CholMatrixType tmp(size,size);
245 ConstCholMatrixPtr pmat;
246 ordering(a, pmat, tmp);
247 analyzePattern_preordered(*pmat,doLDLT);
249 void analyzePattern_preordered(
const CholMatrixType& a,
bool doLDLT);
251 void ordering(
const MatrixType& a, ConstCholMatrixPtr &pmat, CholMatrixType& ap);
255 inline bool operator() (
const Index& row,
const Index& col,
const Scalar&)
const
262 bool m_factorizationIsOk;
272 RealScalar m_shiftOffset;
273 RealScalar m_shiftScale;
276template<
typename MatrixType_,
int UpLo_ = Lower,
typename Ordering_ = AMDOrdering<
typename MatrixType_::StorageIndex> >
class SimplicialLLT;
277template<
typename MatrixType_,
int UpLo_ = Lower,
typename Ordering_ = AMDOrdering<
typename MatrixType_::StorageIndex> >
class SimplicialLDLT;
278template<
typename MatrixType_,
int UpLo_ = Lower,
typename Ordering_ = AMDOrdering<
typename MatrixType_::StorageIndex> >
class SimplicialCholesky;
282template<
typename MatrixType_,
int UpLo_,
typename Ordering_>
struct traits<
SimplicialLLT<MatrixType_,UpLo_,Ordering_> >
284 typedef MatrixType_ MatrixType;
285 typedef Ordering_ OrderingType;
286 enum { UpLo = UpLo_ };
287 typedef typename MatrixType::Scalar Scalar;
288 typedef typename MatrixType::StorageIndex StorageIndex;
289 typedef SparseMatrix<Scalar, ColMajor, StorageIndex> CholMatrixType;
290 typedef TriangularView<const CholMatrixType, Eigen::Lower> MatrixL;
291 typedef TriangularView<const typename CholMatrixType::AdjointReturnType, Eigen::Upper> MatrixU;
292 static inline MatrixL getL(
const CholMatrixType& m) {
return MatrixL(m); }
293 static inline MatrixU getU(
const CholMatrixType& m) {
return MatrixU(m.adjoint()); }
296template<
typename MatrixType_,
int UpLo_,
typename Ordering_>
struct traits<SimplicialLDLT<MatrixType_,UpLo_,Ordering_> >
298 typedef MatrixType_ MatrixType;
299 typedef Ordering_ OrderingType;
300 enum { UpLo = UpLo_ };
301 typedef typename MatrixType::Scalar Scalar;
302 typedef typename MatrixType::StorageIndex StorageIndex;
303 typedef SparseMatrix<Scalar, ColMajor, StorageIndex> CholMatrixType;
304 typedef TriangularView<const CholMatrixType, Eigen::UnitLower> MatrixL;
305 typedef TriangularView<const typename CholMatrixType::AdjointReturnType, Eigen::UnitUpper> MatrixU;
306 static inline MatrixL getL(
const CholMatrixType& m) {
return MatrixL(m); }
307 static inline MatrixU getU(
const CholMatrixType& m) {
return MatrixU(m.adjoint()); }
310template<
typename MatrixType_,
int UpLo_,
typename Ordering_>
struct traits<SimplicialCholesky<MatrixType_,UpLo_,Ordering_> >
312 typedef MatrixType_ MatrixType;
313 typedef Ordering_ OrderingType;
314 enum { UpLo = UpLo_ };
339template<
typename MatrixType_,
int UpLo_,
typename Ordering_>
343 typedef MatrixType_ MatrixType;
344 enum { UpLo = UpLo_ };
346 typedef typename MatrixType::Scalar Scalar;
347 typedef typename MatrixType::RealScalar RealScalar;
348 typedef typename MatrixType::StorageIndex StorageIndex;
351 typedef internal::traits<SimplicialLLT> Traits;
352 typedef typename Traits::MatrixL MatrixL;
353 typedef typename Traits::MatrixU MatrixU;
363 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial LLT not factorized");
364 return Traits::getL(Base::m_matrix);
369 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial LLT not factorized");
370 return Traits::getU(Base::m_matrix);
376 Base::template compute<false>(matrix);
388 Base::analyzePattern(a,
false);
399 Base::template factorize<false>(a);
405 Scalar detL = Base::m_matrix.diagonal().prod();
406 return numext::abs2(detL);
430template<
typename MatrixType_,
int UpLo_,
typename Ordering_>
434 typedef MatrixType_ MatrixType;
435 enum { UpLo = UpLo_ };
437 typedef typename MatrixType::Scalar Scalar;
438 typedef typename MatrixType::RealScalar RealScalar;
439 typedef typename MatrixType::StorageIndex StorageIndex;
442 typedef internal::traits<SimplicialLDLT> Traits;
443 typedef typename Traits::MatrixL MatrixL;
444 typedef typename Traits::MatrixU MatrixU;
455 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial LDLT not factorized");
460 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial LDLT not factorized");
461 return Traits::getL(Base::m_matrix);
466 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial LDLT not factorized");
467 return Traits::getU(Base::m_matrix);
473 Base::template compute<true>(matrix);
485 Base::analyzePattern(a,
true);
496 Base::template factorize<true>(a);
502 return Base::m_diag.prod();
512template<
typename MatrixType_,
int UpLo_,
typename Ordering_>
516 typedef MatrixType_ MatrixType;
517 enum { UpLo = UpLo_ };
519 typedef typename MatrixType::Scalar Scalar;
520 typedef typename MatrixType::RealScalar RealScalar;
521 typedef typename MatrixType::StorageIndex StorageIndex;
524 typedef internal::traits<SimplicialCholesky> Traits;
525 typedef internal::traits<SimplicialLDLT<MatrixType,UpLo> > LDLTTraits;
526 typedef internal::traits<SimplicialLLT<MatrixType,UpLo> > LLTTraits;
531 :
Base(), m_LDLT(
true)
540 case SimplicialCholeskyLLT:
543 case SimplicialCholeskyLDLT:
554 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial Cholesky not factorized");
558 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial Cholesky not factorized");
559 return Base::m_matrix;
566 Base::template compute<true>(matrix);
568 Base::template compute<false>(matrix);
580 Base::analyzePattern(a, m_LDLT);
592 Base::template factorize<true>(a);
594 Base::template factorize<false>(a);
598 template<
typename Rhs,
typename Dest>
601 eigen_assert(Base::m_factorizationIsOk &&
"The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
602 eigen_assert(Base::m_matrix.rows()==b.
rows());
607 if(Base::m_P.size()>0)
608 dest = Base::m_P * b;
612 if(Base::m_matrix.nonZeros()>0)
615 LDLTTraits::getL(Base::m_matrix).solveInPlace(dest);
617 LLTTraits::getL(Base::m_matrix).solveInPlace(dest);
620 if(Base::m_diag.size()>0)
621 dest = Base::m_diag.real().
asDiagonal().inverse() * dest;
623 if (Base::m_matrix.nonZeros()>0)
626 LDLTTraits::getU(Base::m_matrix).solveInPlace(dest);
628 LLTTraits::getU(Base::m_matrix).solveInPlace(dest);
631 if(Base::m_P.size()>0)
632 dest = Base::m_Pinv * dest;
636 template<
typename Rhs,
typename Dest>
637 void _solve_impl(
const SparseMatrixBase<Rhs> &b, SparseMatrixBase<Dest> &dest)
const
639 internal::solve_sparse_through_dense_panels(*
this, b, dest);
642 Scalar determinant()
const
646 return Base::m_diag.prod();
650 Scalar detL = Diagonal<const CholMatrixType>(Base::m_matrix).prod();
651 return numext::abs2(detL);
659template<
typename Derived>
660void SimplicialCholeskyBase<Derived>::ordering(
const MatrixType& a, ConstCholMatrixPtr &pmat, CholMatrixType& ap)
662 eigen_assert(a.rows()==a.cols());
663 const Index size = a.rows();
666 if(!internal::is_same<OrderingType,NaturalOrdering<Index> >::value)
670 C = a.template selfadjointView<UpLo>();
672 OrderingType ordering;
676 if(m_Pinv.size()>0) m_P = m_Pinv.inverse();
679 ap.resize(size,size);
680 ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
686 if(
int(UpLo)==
int(
Lower) || MatrixType::IsRowMajor)
689 ap.resize(size,size);
690 ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>();
693 internal::simplicial_cholesky_grab_input<CholMatrixType,MatrixType>::run(a, pmat, ap);
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: EigenBase.h:62
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:52
const DiagonalWrapper< const Derived > asDiagonal() const
Definition: DiagonalMatrix.h:325
Index size() const
Definition: PermutationMatrix.h:99
A base class for direct sparse Cholesky factorizations.
Definition: SimplicialCholesky.h:58
SimplicialCholeskyBase()
Definition: SimplicialCholesky.h:84
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SimplicialCholesky.h:117
const PermutationMatrix< Dynamic, Dynamic, StorageIndex > & permutationP() const
Definition: SimplicialCholesky.h:125
Derived & setShift(const RealScalar &offset, const RealScalar &scale=1)
Definition: SimplicialCholesky.h:142
void compute(const MatrixType &matrix)
Definition: SimplicialCholesky.h:204
const PermutationMatrix< Dynamic, Dynamic, StorageIndex > & permutationPinv() const
Definition: SimplicialCholesky.h:130
Definition: SimplicialCholesky.h:514
void factorize(const MatrixType &a)
Definition: SimplicialCholesky.h:589
void analyzePattern(const MatrixType &a)
Definition: SimplicialCholesky.h:578
SimplicialCholesky & compute(const MatrixType &matrix)
Definition: SimplicialCholesky.h:563
A direct sparse LDLT Cholesky factorizations without square root.
Definition: SimplicialCholesky.h:432
Scalar determinant() const
Definition: SimplicialCholesky.h:500
const MatrixL matrixL() const
Definition: SimplicialCholesky.h:459
SimplicialLDLT(const MatrixType &matrix)
Definition: SimplicialCholesky.h:450
const MatrixU matrixU() const
Definition: SimplicialCholesky.h:465
void analyzePattern(const MatrixType &a)
Definition: SimplicialCholesky.h:483
void factorize(const MatrixType &a)
Definition: SimplicialCholesky.h:494
const VectorType vectorD() const
Definition: SimplicialCholesky.h:454
SimplicialLDLT & compute(const MatrixType &matrix)
Definition: SimplicialCholesky.h:471
SimplicialLDLT()
Definition: SimplicialCholesky.h:447
A direct sparse LLT Cholesky factorizations.
Definition: SimplicialCholesky.h:341
const MatrixL matrixL() const
Definition: SimplicialCholesky.h:362
SimplicialLLT & compute(const MatrixType &matrix)
Definition: SimplicialCholesky.h:374
Scalar determinant() const
Definition: SimplicialCholesky.h:403
const MatrixU matrixU() const
Definition: SimplicialCholesky.h:368
SimplicialLLT()
Definition: SimplicialCholesky.h:356
void analyzePattern(const MatrixType &a)
Definition: SimplicialCholesky.h:386
void factorize(const MatrixType &a)
Definition: SimplicialCholesky.h:397
SimplicialLLT(const MatrixType &matrix)
Definition: SimplicialCholesky.h:358
Index nonZeros() const
Definition: SparseCompressedBase.h:58
Index cols() const
Definition: SparseMatrix.h:142
Index rows() const
Definition: SparseMatrix.h:140
A base class for sparse solvers.
Definition: SparseSolverBase.h:70
ComputationInfo
Definition: Constants.h:442
@ Lower
Definition: Constants.h:211
@ Upper
Definition: Constants.h:213
@ Success
Definition: Constants.h:444
Namespace containing all symbols from the Eigen library.
Definition: B01_Experimental.dox:1
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
Definition: SimplicialCholesky.h:254