10 #ifndef EIGEN_SIMPLICIAL_CHOLESKY_H
11 #define EIGEN_SIMPLICIAL_CHOLESKY_H
15 enum SimplicialCholeskyMode {
16 SimplicialCholeskyLLT,
17 SimplicialCholeskyLDLT
35 template<
typename Derived>
39 typedef typename internal::traits<Derived>::MatrixType MatrixType;
40 enum { UpLo = internal::traits<Derived>::UpLo };
41 typedef typename MatrixType::Scalar Scalar;
42 typedef typename MatrixType::RealScalar RealScalar;
43 typedef typename MatrixType::Index Index;
51 : m_info(
Success), m_isInitialized(false), m_shiftOffset(0), m_shiftScale(1)
55 : m_info(
Success), m_isInitialized(false), m_shiftOffset(0), m_shiftScale(1)
57 derived().compute(matrix);
64 Derived& derived() {
return *
static_cast<Derived*
>(
this); }
65 const Derived& derived()
const {
return *
static_cast<const Derived*
>(
this); }
67 inline Index cols()
const {
return m_matrix.
cols(); }
68 inline Index rows()
const {
return m_matrix.
rows(); }
77 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
85 template<
typename Rhs>
86 inline const internal::solve_retval<SimplicialCholeskyBase, Rhs>
89 eigen_assert(m_isInitialized &&
"Simplicial LLT or LDLT is not initialized.");
90 eigen_assert(rows()==b.rows()
91 &&
"SimplicialCholeskyBase::solve(): invalid number of rows of the right hand side matrix b");
92 return internal::solve_retval<SimplicialCholeskyBase, Rhs>(*
this, b.derived());
99 template<
typename Rhs>
100 inline const internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs>
103 eigen_assert(m_isInitialized &&
"Simplicial LLT or LDLT is not initialized.");
104 eigen_assert(rows()==b.
rows()
105 &&
"SimplicialCholesky::solve(): invalid number of rows of the right hand side matrix b");
106 return internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs>(*
this, b.
derived());
128 Derived&
setShift(
const RealScalar& offset,
const RealScalar& scale = 1)
130 m_shiftOffset = offset;
131 m_shiftScale = scale;
135 #ifndef EIGEN_PARSED_BY_DOXYGEN
137 template<
typename Stream>
138 void dumpMemory(Stream& s)
141 s <<
" L: " << ((total+=(m_matrix.
cols()+1) *
sizeof(
int) + m_matrix.
nonZeros()*(
sizeof(int)+
sizeof(Scalar))) >> 20) <<
"Mb" <<
"\n";
142 s <<
" diag: " << ((total+=m_diag.size() *
sizeof(Scalar)) >> 20) <<
"Mb" <<
"\n";
143 s <<
" tree: " << ((total+=m_parent.size() *
sizeof(int)) >> 20) <<
"Mb" <<
"\n";
144 s <<
" nonzeros: " << ((total+=m_nonZerosPerCol.size() *
sizeof(int)) >> 20) <<
"Mb" <<
"\n";
145 s <<
" perm: " << ((total+=m_P.
size() *
sizeof(int)) >> 20) <<
"Mb" <<
"\n";
146 s <<
" perm^-1: " << ((total+=m_Pinv.
size() *
sizeof(int)) >> 20) <<
"Mb" <<
"\n";
147 s <<
" TOTAL: " << (total>> 20) <<
"Mb" <<
"\n";
151 template<
typename Rhs,
typename Dest>
152 void _solve(
const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest)
const
154 eigen_assert(m_factorizationIsOk &&
"The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
155 eigen_assert(m_matrix.
rows()==b.rows());
166 derived().matrixL().solveInPlace(dest);
169 dest = m_diag.asDiagonal().inverse() * dest;
172 derived().matrixU().solveInPlace(dest);
175 dest = m_Pinv * dest;
178 #endif // EIGEN_PARSED_BY_DOXYGEN
183 template<
bool DoLDLT>
186 eigen_assert(matrix.rows()==matrix.cols());
187 Index size = matrix.cols();
189 ordering(matrix, ap);
190 analyzePattern_preordered(ap, DoLDLT);
191 factorize_preordered<DoLDLT>(ap);
194 template<
bool DoLDLT>
195 void factorize(
const MatrixType& a)
197 eigen_assert(a.rows()==a.cols());
199 CholMatrixType ap(size,size);
200 ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
201 factorize_preordered<DoLDLT>(ap);
204 template<
bool DoLDLT>
205 void factorize_preordered(
const CholMatrixType& a);
207 void analyzePattern(
const MatrixType& a,
bool doLDLT)
209 eigen_assert(a.rows()==a.cols());
211 CholMatrixType ap(size,size);
213 analyzePattern_preordered(ap,doLDLT);
215 void analyzePattern_preordered(
const CholMatrixType& a,
bool doLDLT);
217 void ordering(
const MatrixType& a, CholMatrixType& ap);
221 inline bool operator() (
const Index& row,
const Index& col,
const Scalar&)
const
228 bool m_isInitialized;
229 bool m_factorizationIsOk;
239 RealScalar m_shiftOffset;
240 RealScalar m_shiftScale;
249 template<
typename _MatrixType,
int _UpLo>
struct traits<
SimplicialLLT<_MatrixType,_UpLo> >
251 typedef _MatrixType MatrixType;
252 enum { UpLo = _UpLo };
253 typedef typename MatrixType::Scalar Scalar;
254 typedef typename MatrixType::Index Index;
256 typedef SparseTriangularView<CholMatrixType, Eigen::Lower> MatrixL;
257 typedef SparseTriangularView<typename CholMatrixType::AdjointReturnType, Eigen::Upper> MatrixU;
258 static inline MatrixL getL(
const MatrixType& m) {
return m; }
259 static inline MatrixU getU(
const MatrixType& m) {
return m.adjoint(); }
262 template<
typename _MatrixType,
int _UpLo>
struct traits<SimplicialLDLT<_MatrixType,_UpLo> >
264 typedef _MatrixType MatrixType;
265 enum { UpLo = _UpLo };
266 typedef typename MatrixType::Scalar Scalar;
267 typedef typename MatrixType::Index Index;
268 typedef SparseMatrix<Scalar, ColMajor, Index> CholMatrixType;
269 typedef SparseTriangularView<CholMatrixType, Eigen::UnitLower> MatrixL;
270 typedef SparseTriangularView<typename CholMatrixType::AdjointReturnType, Eigen::UnitUpper> MatrixU;
271 static inline MatrixL getL(
const MatrixType& m) {
return m; }
272 static inline MatrixU getU(
const MatrixType& m) {
return m.adjoint(); }
275 template<
typename _MatrixType,
int _UpLo>
struct traits<SimplicialCholesky<_MatrixType,_UpLo> >
277 typedef _MatrixType MatrixType;
278 enum { UpLo = _UpLo };
300 template<
typename _MatrixType,
int _UpLo>
301 class SimplicialLLT :
public SimplicialCholeskyBase<SimplicialLLT<_MatrixType,_UpLo> >
304 typedef _MatrixType MatrixType;
305 enum { UpLo = _UpLo };
306 typedef SimplicialCholeskyBase<SimplicialLLT> Base;
307 typedef typename MatrixType::Scalar Scalar;
308 typedef typename MatrixType::RealScalar RealScalar;
309 typedef typename MatrixType::Index Index;
310 typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
311 typedef Matrix<Scalar,Dynamic,1> VectorType;
312 typedef internal::traits<SimplicialLLT> Traits;
313 typedef typename Traits::MatrixL MatrixL;
314 typedef typename Traits::MatrixU MatrixU;
324 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial LLT not factorized");
325 return Traits::getL(Base::m_matrix);
330 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial LLT not factorized");
331 return Traits::getU(Base::m_matrix);
337 Base::template compute<false>(matrix);
349 Base::analyzePattern(a,
false);
360 Base::template factorize<false>(a);
366 Scalar detL = Base::m_matrix.diagonal().prod();
367 return numext::abs2(detL);
388 template<
typename _MatrixType,
int _UpLo>
389 class SimplicialLDLT :
public SimplicialCholeskyBase<SimplicialLDLT<_MatrixType,_UpLo> >
392 typedef _MatrixType MatrixType;
393 enum { UpLo = _UpLo };
394 typedef SimplicialCholeskyBase<SimplicialLDLT> Base;
395 typedef typename MatrixType::Scalar Scalar;
396 typedef typename MatrixType::RealScalar RealScalar;
397 typedef typename MatrixType::Index Index;
398 typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
399 typedef Matrix<Scalar,Dynamic,1> VectorType;
400 typedef internal::traits<SimplicialLDLT> Traits;
401 typedef typename Traits::MatrixL MatrixL;
402 typedef typename Traits::MatrixU MatrixU;
413 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial LDLT not factorized");
418 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial LDLT not factorized");
419 return Traits::getL(Base::m_matrix);
424 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial LDLT not factorized");
425 return Traits::getU(Base::m_matrix);
431 Base::template compute<true>(matrix);
443 Base::analyzePattern(a,
true);
454 Base::template factorize<true>(a);
460 return Base::m_diag.prod();
470 template<
typename _MatrixType,
int _UpLo>
471 class SimplicialCholesky :
public SimplicialCholeskyBase<SimplicialCholesky<_MatrixType,_UpLo> >
474 typedef _MatrixType MatrixType;
475 enum { UpLo = _UpLo };
476 typedef SimplicialCholeskyBase<SimplicialCholesky> Base;
477 typedef typename MatrixType::Scalar Scalar;
478 typedef typename MatrixType::RealScalar RealScalar;
479 typedef typename MatrixType::Index Index;
480 typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
481 typedef Matrix<Scalar,Dynamic,1> VectorType;
482 typedef internal::traits<SimplicialCholesky> Traits;
483 typedef internal::traits<SimplicialLDLT<MatrixType,UpLo> > LDLTTraits;
484 typedef internal::traits<SimplicialLLT<MatrixType,UpLo> > LLTTraits;
486 SimplicialCholesky() : Base(), m_LDLT(true) {}
488 SimplicialCholesky(
const MatrixType& matrix)
489 : Base(), m_LDLT(true)
494 SimplicialCholesky& setMode(SimplicialCholeskyMode mode)
498 case SimplicialCholeskyLLT:
501 case SimplicialCholeskyLDLT:
511 inline const VectorType vectorD()
const {
512 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial Cholesky not factorized");
515 inline const CholMatrixType rawMatrix()
const {
516 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial Cholesky not factorized");
517 return Base::m_matrix;
524 Base::template compute<true>(matrix);
526 Base::template compute<false>(matrix);
538 Base::analyzePattern(a, m_LDLT);
550 Base::template factorize<true>(a);
552 Base::template factorize<false>(a);
556 template<
typename Rhs,
typename Dest>
559 eigen_assert(Base::m_factorizationIsOk &&
"The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
560 eigen_assert(Base::m_matrix.rows()==b.rows());
565 if(Base::m_P.size()>0)
566 dest = Base::m_P * b;
570 if(Base::m_matrix.nonZeros()>0)
573 LDLTTraits::getL(Base::m_matrix).solveInPlace(dest);
575 LLTTraits::getL(Base::m_matrix).solveInPlace(dest);
578 if(Base::m_diag.size()>0)
579 dest = Base::m_diag.
asDiagonal().inverse() * dest;
581 if (Base::m_matrix.nonZeros()>0)
584 LDLTTraits::getU(Base::m_matrix).solveInPlace(dest);
586 LLTTraits::getU(Base::m_matrix).solveInPlace(dest);
589 if(Base::m_P.size()>0)
590 dest = Base::m_Pinv * dest;
593 Scalar determinant()
const
597 return Base::m_diag.
prod();
601 Scalar detL = Diagonal<const CholMatrixType>(Base::m_matrix).prod();
602 return numext::abs2(detL);
610 template<
typename Derived>
611 void SimplicialCholeskyBase<Derived>::ordering(
const MatrixType& a, CholMatrixType& ap)
613 eigen_assert(a.rows()==a.cols());
614 const Index size = a.rows();
619 C = a.template selfadjointView<UpLo>();
623 internal::minimum_degree_ordering(C, m_Pinv);
627 m_P = m_Pinv.inverse();
631 ap.resize(size,size);
632 ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
637 template<
typename Derived,
typename Rhs>
638 struct solve_retval<SimplicialCholeskyBase<Derived>, Rhs>
639 : solve_retval_base<SimplicialCholeskyBase<Derived>, Rhs>
641 typedef SimplicialCholeskyBase<Derived> Dec;
642 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
644 template<typename Dest>
void evalTo(Dest& dst)
const
646 dec().derived()._solve(rhs(),dst);
650 template<
typename Derived,
typename Rhs>
651 struct sparse_solve_retval<SimplicialCholeskyBase<Derived>, Rhs>
652 : sparse_solve_retval_base<SimplicialCholeskyBase<Derived>, Rhs>
654 typedef SimplicialCholeskyBase<Derived> Dec;
655 EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
657 template<typename Dest>
void evalTo(Dest& dst)
const
659 this->defaultEvalTo(dst);
667 #endif // EIGEN_SIMPLICIAL_CHOLESKY_H