10#ifndef EIGEN_SIMPLICIAL_CHOLESKY_H
11#define EIGEN_SIMPLICIAL_CHOLESKY_H
15enum 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 &)
54template<
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
84 m_factorizationIsOk(false),
85 m_analysisIsOk(false),
92 m_factorizationIsOk(false),
93 m_analysisIsOk(false),
97 derived().compute(matrix);
100 ~SimplicialCholeskyBase()
104 Derived& derived() {
return *
static_cast<Derived*
>(
this); }
105 const Derived& derived()
const {
return *
static_cast<const Derived*
>(
this); }
107 inline Index cols()
const {
return m_matrix.
cols(); }
108 inline Index rows()
const {
return m_matrix.
rows(); }
117 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
140 Derived&
setShift(
const RealScalar& offset,
const RealScalar& scale = 1)
142 m_shiftOffset = offset;
143 m_shiftScale = scale;
147#ifndef EIGEN_PARSED_BY_DOXYGEN
149 template<
typename Stream>
150 void dumpMemory(Stream& s)
153 s <<
" L: " << ((total+=(m_matrix.
cols()+1) *
sizeof(
int) + m_matrix.
nonZeros()*(
sizeof(int)+
sizeof(Scalar))) >> 20) <<
"Mb" <<
"\n";
154 s <<
" diag: " << ((total+=m_diag.size() *
sizeof(Scalar)) >> 20) <<
"Mb" <<
"\n";
155 s <<
" tree: " << ((total+=m_parent.size() *
sizeof(int)) >> 20) <<
"Mb" <<
"\n";
156 s <<
" nonzeros: " << ((total+=m_nonZerosPerCol.size() *
sizeof(int)) >> 20) <<
"Mb" <<
"\n";
157 s <<
" perm: " << ((total+=m_P.
size() *
sizeof(int)) >> 20) <<
"Mb" <<
"\n";
158 s <<
" perm^-1: " << ((total+=m_Pinv.
size() *
sizeof(int)) >> 20) <<
"Mb" <<
"\n";
159 s <<
" TOTAL: " << (total>> 20) <<
"Mb" <<
"\n";
163 template<
typename Rhs,
typename Dest>
164 void _solve_impl(
const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest)
const
166 eigen_assert(m_factorizationIsOk &&
"The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
167 eigen_assert(m_matrix.
rows()==b.rows());
178 derived().matrixL().solveInPlace(dest);
181 dest = m_diag.asDiagonal().inverse() * dest;
184 derived().matrixU().solveInPlace(dest);
187 dest = m_Pinv * dest;
190 template<
typename Rhs,
typename Dest>
191 void _solve_impl(
const SparseMatrixBase<Rhs> &b, SparseMatrixBase<Dest> &dest)
const
193 internal::solve_sparse_through_dense_panels(derived(), b, dest);
201 template<
bool DoLDLT>
204 eigen_assert(matrix.rows()==matrix.cols());
205 Index size = matrix.cols();
207 ConstCholMatrixPtr pmat;
208 ordering(matrix, pmat, tmp);
209 analyzePattern_preordered(*pmat, DoLDLT);
210 factorize_preordered<DoLDLT>(*pmat);
213 template<
bool DoLDLT>
214 void factorize(
const MatrixType& a)
216 eigen_assert(a.rows()==a.cols());
217 Index size = a.cols();
218 CholMatrixType tmp(size,size);
219 ConstCholMatrixPtr pmat;
224 internal::simplicial_cholesky_grab_input<CholMatrixType,MatrixType>::run(a, pmat, tmp);
228 tmp.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
232 factorize_preordered<DoLDLT>(*pmat);
235 template<
bool DoLDLT>
236 void factorize_preordered(
const CholMatrixType& a);
238 void analyzePattern(
const MatrixType& a,
bool doLDLT)
240 eigen_assert(a.rows()==a.cols());
241 Index size = a.cols();
242 CholMatrixType tmp(size,size);
243 ConstCholMatrixPtr pmat;
244 ordering(a, pmat, tmp);
245 analyzePattern_preordered(*pmat,doLDLT);
247 void analyzePattern_preordered(
const CholMatrixType& a,
bool doLDLT);
249 void ordering(
const MatrixType& a, ConstCholMatrixPtr &pmat, CholMatrixType& ap);
253 inline bool operator() (
const Index& row,
const Index& col,
const Scalar&)
const
260 bool m_factorizationIsOk;
270 RealScalar m_shiftOffset;
271 RealScalar m_shiftScale;
274template<
typename _MatrixType,
int _UpLo = Lower,
typename _Ordering = AMDOrdering<
typename _MatrixType::StorageIndex> >
class SimplicialLLT;
275template<
typename _MatrixType,
int _UpLo = Lower,
typename _Ordering = AMDOrdering<
typename _MatrixType::StorageIndex> >
class SimplicialLDLT;
276template<
typename _MatrixType,
int _UpLo = Lower,
typename _Ordering = AMDOrdering<
typename _MatrixType::StorageIndex> >
class SimplicialCholesky;
280template<
typename _MatrixType,
int _UpLo,
typename _Ordering>
struct traits<
SimplicialLLT<_MatrixType,_UpLo,_Ordering> >
282 typedef _MatrixType MatrixType;
283 typedef _Ordering OrderingType;
284 enum { UpLo = _UpLo };
285 typedef typename MatrixType::Scalar Scalar;
286 typedef typename MatrixType::StorageIndex StorageIndex;
287 typedef SparseMatrix<Scalar, ColMajor, StorageIndex> CholMatrixType;
288 typedef TriangularView<const CholMatrixType, Eigen::Lower> MatrixL;
289 typedef TriangularView<const typename CholMatrixType::AdjointReturnType, Eigen::Upper> MatrixU;
290 static inline MatrixL getL(
const CholMatrixType& m) {
return MatrixL(m); }
291 static inline MatrixU getU(
const CholMatrixType& m) {
return MatrixU(m.adjoint()); }
294template<
typename _MatrixType,
int _UpLo,
typename _Ordering>
struct traits<SimplicialLDLT<_MatrixType,_UpLo,_Ordering> >
296 typedef _MatrixType MatrixType;
297 typedef _Ordering OrderingType;
298 enum { UpLo = _UpLo };
299 typedef typename MatrixType::Scalar Scalar;
300 typedef typename MatrixType::StorageIndex StorageIndex;
301 typedef SparseMatrix<Scalar, ColMajor, StorageIndex> CholMatrixType;
302 typedef TriangularView<const CholMatrixType, Eigen::UnitLower> MatrixL;
303 typedef TriangularView<const typename CholMatrixType::AdjointReturnType, Eigen::UnitUpper> MatrixU;
304 static inline MatrixL getL(
const CholMatrixType& m) {
return MatrixL(m); }
305 static inline MatrixU getU(
const CholMatrixType& m) {
return MatrixU(m.adjoint()); }
308template<
typename _MatrixType,
int _UpLo,
typename _Ordering>
struct traits<SimplicialCholesky<_MatrixType,_UpLo,_Ordering> >
310 typedef _MatrixType MatrixType;
311 typedef _Ordering OrderingType;
312 enum { UpLo = _UpLo };
337template<
typename _MatrixType,
int _UpLo,
typename _Ordering>
341 typedef _MatrixType MatrixType;
342 enum { UpLo = _UpLo };
344 typedef typename MatrixType::Scalar Scalar;
345 typedef typename MatrixType::RealScalar RealScalar;
346 typedef typename MatrixType::StorageIndex StorageIndex;
349 typedef internal::traits<SimplicialLLT> Traits;
350 typedef typename Traits::MatrixL MatrixL;
351 typedef typename Traits::MatrixU MatrixU;
361 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial LLT not factorized");
362 return Traits::getL(Base::m_matrix);
367 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial LLT not factorized");
368 return Traits::getU(Base::m_matrix);
374 Base::template compute<false>(matrix);
386 Base::analyzePattern(a,
false);
397 Base::template factorize<false>(a);
403 Scalar detL = Base::m_matrix.diagonal().prod();
404 return numext::abs2(detL);
428template<
typename _MatrixType,
int _UpLo,
typename _Ordering>
432 typedef _MatrixType MatrixType;
433 enum { UpLo = _UpLo };
435 typedef typename MatrixType::Scalar Scalar;
436 typedef typename MatrixType::RealScalar RealScalar;
437 typedef typename MatrixType::StorageIndex StorageIndex;
440 typedef internal::traits<SimplicialLDLT> Traits;
441 typedef typename Traits::MatrixL MatrixL;
442 typedef typename Traits::MatrixU MatrixU;
453 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial LDLT not factorized");
458 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial LDLT not factorized");
459 return Traits::getL(Base::m_matrix);
464 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial LDLT not factorized");
465 return Traits::getU(Base::m_matrix);
471 Base::template compute<true>(matrix);
483 Base::analyzePattern(a,
true);
494 Base::template factorize<true>(a);
500 return Base::m_diag.prod();
510template<
typename _MatrixType,
int _UpLo,
typename _Ordering>
514 typedef _MatrixType MatrixType;
515 enum { UpLo = _UpLo };
517 typedef typename MatrixType::Scalar Scalar;
518 typedef typename MatrixType::RealScalar RealScalar;
519 typedef typename MatrixType::StorageIndex StorageIndex;
522 typedef internal::traits<SimplicialCholesky> Traits;
523 typedef internal::traits<SimplicialLDLT<MatrixType,UpLo> > LDLTTraits;
524 typedef internal::traits<SimplicialLLT<MatrixType,UpLo> > LLTTraits;
529 :
Base(), m_LDLT(
true)
538 case SimplicialCholeskyLLT:
541 case SimplicialCholeskyLDLT:
552 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial Cholesky not factorized");
556 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial Cholesky not factorized");
557 return Base::m_matrix;
564 Base::template compute<true>(matrix);
566 Base::template compute<false>(matrix);
578 Base::analyzePattern(a, m_LDLT);
590 Base::template factorize<true>(a);
592 Base::template factorize<false>(a);
596 template<
typename Rhs,
typename Dest>
599 eigen_assert(Base::m_factorizationIsOk &&
"The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
600 eigen_assert(Base::m_matrix.rows()==b.
rows());
605 if(Base::m_P.size()>0)
606 dest = Base::m_P * b;
610 if(Base::m_matrix.nonZeros()>0)
613 LDLTTraits::getL(Base::m_matrix).solveInPlace(dest);
615 LLTTraits::getL(Base::m_matrix).solveInPlace(dest);
618 if(Base::m_diag.size()>0)
619 dest = Base::m_diag.real().
asDiagonal().inverse() * dest;
621 if (Base::m_matrix.nonZeros()>0)
624 LDLTTraits::getU(Base::m_matrix).solveInPlace(dest);
626 LLTTraits::getU(Base::m_matrix).solveInPlace(dest);
629 if(Base::m_P.size()>0)
630 dest = Base::m_Pinv * dest;
634 template<
typename Rhs,
typename Dest>
635 void _solve_impl(
const SparseMatrixBase<Rhs> &b, SparseMatrixBase<Dest> &dest)
const
637 internal::solve_sparse_through_dense_panels(*
this, b, dest);
640 Scalar determinant()
const
644 return Base::m_diag.prod();
648 Scalar detL = Diagonal<const CholMatrixType>(Base::m_matrix).prod();
649 return numext::abs2(detL);
657template<
typename Derived>
658void SimplicialCholeskyBase<Derived>::ordering(
const MatrixType& a, ConstCholMatrixPtr &pmat, CholMatrixType& ap)
660 eigen_assert(a.rows()==a.cols());
661 const Index size = a.rows();
664 if(!internal::is_same<OrderingType,NaturalOrdering<Index> >::value)
668 C = a.template selfadjointView<UpLo>();
670 OrderingType ordering;
674 if(m_Pinv.size()>0) m_P = m_Pinv.inverse();
677 ap.resize(size,size);
678 ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
684 if(
int(UpLo)==
int(Lower) || MatrixType::IsRowMajor)
687 ap.resize(size,size);
688 ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>();
691 internal::simplicial_cholesky_grab_input<CholMatrixType,MatrixType>::run(a, pmat, ap);
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: EigenBase.h:60
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:50
const DiagonalWrapper< const Derived > asDiagonal() const
Definition: DiagonalMatrix.h:325
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:180
Index size() const
Definition: PermutationMatrix.h:97
Permutation matrix.
Definition: PermutationMatrix.h:298
A base class for direct sparse Cholesky factorizations.
Definition: SimplicialCholesky.h:56
SimplicialCholeskyBase()
Definition: SimplicialCholesky.h:82
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SimplicialCholesky.h:115
const PermutationMatrix< Dynamic, Dynamic, StorageIndex > & permutationP() const
Definition: SimplicialCholesky.h:123
Derived & setShift(const RealScalar &offset, const RealScalar &scale=1)
Definition: SimplicialCholesky.h:140
void compute(const MatrixType &matrix)
Definition: SimplicialCholesky.h:202
const PermutationMatrix< Dynamic, Dynamic, StorageIndex > & permutationPinv() const
Definition: SimplicialCholesky.h:128
Definition: SimplicialCholesky.h:512
SimplicialCholesky & compute(const MatrixType &matrix)
Definition: SimplicialCholesky.h:561
void analyzePattern(const MatrixType &a)
Definition: SimplicialCholesky.h:576
void factorize(const MatrixType &a)
Definition: SimplicialCholesky.h:587
A direct sparse LDLT Cholesky factorizations without square root.
Definition: SimplicialCholesky.h:430
SimplicialLDLT(const MatrixType &matrix)
Definition: SimplicialCholesky.h:448
SimplicialLDLT()
Definition: SimplicialCholesky.h:445
SimplicialLDLT & compute(const MatrixType &matrix)
Definition: SimplicialCholesky.h:469
void factorize(const MatrixType &a)
Definition: SimplicialCholesky.h:492
Scalar determinant() const
Definition: SimplicialCholesky.h:498
void analyzePattern(const MatrixType &a)
Definition: SimplicialCholesky.h:481
const VectorType vectorD() const
Definition: SimplicialCholesky.h:452
const MatrixL matrixL() const
Definition: SimplicialCholesky.h:457
const MatrixU matrixU() const
Definition: SimplicialCholesky.h:463
A direct sparse LLT Cholesky factorizations.
Definition: SimplicialCholesky.h:339
const MatrixU matrixU() const
Definition: SimplicialCholesky.h:366
SimplicialLLT(const MatrixType &matrix)
Definition: SimplicialCholesky.h:356
SimplicialLLT & compute(const MatrixType &matrix)
Definition: SimplicialCholesky.h:372
void factorize(const MatrixType &a)
Definition: SimplicialCholesky.h:395
Scalar determinant() const
Definition: SimplicialCholesky.h:401
SimplicialLLT()
Definition: SimplicialCholesky.h:354
void analyzePattern(const MatrixType &a)
Definition: SimplicialCholesky.h:384
const MatrixL matrixL() const
Definition: SimplicialCholesky.h:360
A versatible sparse matrix representation.
Definition: SparseMatrix.h:98
Index nonZeros() const
Definition: SparseCompressedBase.h:56
Index rows() const
Definition: SparseMatrix.h:138
Index cols() const
Definition: SparseMatrix.h:140
A base class for sparse solvers.
Definition: SparseSolverBase.h:68
ComputationInfo
Definition: Constants.h:440
@ Upper
Definition: Constants.h:211
@ Success
Definition: Constants.h:442
Namespace containing all symbols from the Eigen library.
Definition: Core:141
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
Definition: SimplicialCholesky.h:252