48 #include "../Core/util/NonMPL2.h"
50 #ifndef EIGEN_SIMPLICIAL_CHOLESKY_H
51 #define EIGEN_SIMPLICIAL_CHOLESKY_H
55 enum SimplicialCholeskyMode {
56 SimplicialCholeskyLLT,
57 SimplicialCholeskyLDLT
75 template<
typename Derived>
79 typedef typename internal::traits<Derived>::MatrixType MatrixType;
80 enum { UpLo = internal::traits<Derived>::UpLo };
81 typedef typename MatrixType::Scalar Scalar;
82 typedef typename MatrixType::RealScalar RealScalar;
83 typedef typename MatrixType::Index Index;
91 : m_info(
Success), m_isInitialized(false), m_shiftOffset(0), m_shiftScale(1)
95 : m_info(
Success), m_isInitialized(false), m_shiftOffset(0), m_shiftScale(1)
97 derived().compute(matrix);
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.");
125 template<
typename Rhs>
126 inline const internal::solve_retval<SimplicialCholeskyBase, Rhs>
129 eigen_assert(m_isInitialized &&
"Simplicial LLT or LDLT is not initialized.");
130 eigen_assert(rows()==b.rows()
131 &&
"SimplicialCholeskyBase::solve(): invalid number of rows of the right hand side matrix b");
132 return internal::solve_retval<SimplicialCholeskyBase, Rhs>(*
this, b.derived());
139 template<
typename Rhs>
140 inline const internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs>
143 eigen_assert(m_isInitialized &&
"Simplicial LLT or LDLT is not initialized.");
144 eigen_assert(rows()==b.
rows()
145 &&
"SimplicialCholesky::solve(): invalid number of rows of the right hand side matrix b");
146 return internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs>(*
this, b.
derived());
168 Derived&
setShift(
const RealScalar& offset,
const RealScalar& scale = 1)
170 m_shiftOffset = offset;
171 m_shiftScale = scale;
175 #ifndef EIGEN_PARSED_BY_DOXYGEN
177 template<
typename Stream>
178 void dumpMemory(Stream& s)
181 s <<
" L: " << ((total+=(m_matrix.
cols()+1) *
sizeof(
int) + m_matrix.
nonZeros()*(
sizeof(int)+
sizeof(Scalar))) >> 20) <<
"Mb" <<
"\n";
182 s <<
" diag: " << ((total+=m_diag.size() *
sizeof(Scalar)) >> 20) <<
"Mb" <<
"\n";
183 s <<
" tree: " << ((total+=m_parent.size() *
sizeof(int)) >> 20) <<
"Mb" <<
"\n";
184 s <<
" nonzeros: " << ((total+=m_nonZerosPerCol.size() *
sizeof(int)) >> 20) <<
"Mb" <<
"\n";
185 s <<
" perm: " << ((total+=m_P.
size() *
sizeof(int)) >> 20) <<
"Mb" <<
"\n";
186 s <<
" perm^-1: " << ((total+=m_Pinv.
size() *
sizeof(int)) >> 20) <<
"Mb" <<
"\n";
187 s <<
" TOTAL: " << (total>> 20) <<
"Mb" <<
"\n";
191 template<
typename Rhs,
typename Dest>
192 void _solve(
const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest)
const
194 eigen_assert(m_factorizationIsOk &&
"The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
195 eigen_assert(m_matrix.
rows()==b.rows());
206 derived().matrixL().solveInPlace(dest);
209 dest = m_diag.asDiagonal().inverse() * dest;
212 derived().matrixU().solveInPlace(dest);
215 dest = m_Pinv * dest;
219 template<
typename Rhs,
typename DestScalar,
int DestOptions,
typename DestIndex>
220 void _solve_sparse(
const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest)
const
222 eigen_assert(m_factorizationIsOk &&
"The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
223 eigen_assert(m_matrix.
rows()==b.rows());
226 static const int NbColsAtOnce = 4;
227 int rhsCols = b.cols();
230 for(
int k=0; k<rhsCols; k+=NbColsAtOnce)
232 int actualCols = std::min<int>(rhsCols-k, NbColsAtOnce);
233 tmp.leftCols(actualCols) = b.middleCols(k,actualCols);
234 tmp.leftCols(actualCols) = derived().solve(tmp.leftCols(actualCols));
235 dest.middleCols(k,actualCols) = tmp.leftCols(actualCols).sparseView();
239 #endif // EIGEN_PARSED_BY_DOXYGEN
244 template<
bool DoLDLT>
247 eigen_assert(matrix.rows()==matrix.cols());
248 Index size = matrix.cols();
250 ordering(matrix, ap);
251 analyzePattern_preordered(ap, DoLDLT);
252 factorize_preordered<DoLDLT>(ap);
255 template<
bool DoLDLT>
256 void factorize(
const MatrixType& a)
258 eigen_assert(a.rows()==a.cols());
260 CholMatrixType ap(size,size);
261 ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
262 factorize_preordered<DoLDLT>(ap);
265 template<
bool DoLDLT>
266 void factorize_preordered(
const CholMatrixType& a);
268 void analyzePattern(
const MatrixType& a,
bool doLDLT)
270 eigen_assert(a.rows()==a.cols());
272 CholMatrixType ap(size,size);
274 analyzePattern_preordered(ap,doLDLT);
276 void analyzePattern_preordered(
const CholMatrixType& a,
bool doLDLT);
278 void ordering(
const MatrixType& a, CholMatrixType& ap);
282 inline bool operator() (
const Index& row,
const Index& col,
const Scalar&)
const
289 bool m_isInitialized;
290 bool m_factorizationIsOk;
300 RealScalar m_shiftOffset;
301 RealScalar m_shiftScale;
304 template<
typename _MatrixType,
int _UpLo = Lower>
class SimplicialLLT;
305 template<
typename _MatrixType,
int _UpLo = Lower>
class SimplicialLDLT;
310 template<
typename _MatrixType,
int _UpLo>
struct traits<
SimplicialLLT<_MatrixType,_UpLo> >
312 typedef _MatrixType MatrixType;
313 enum { UpLo = _UpLo };
314 typedef typename MatrixType::Scalar Scalar;
315 typedef typename MatrixType::Index Index;
316 typedef SparseMatrix<Scalar, ColMajor, Index> CholMatrixType;
317 typedef SparseTriangularView<CholMatrixType, Eigen::Lower> MatrixL;
318 typedef SparseTriangularView<typename CholMatrixType::AdjointReturnType, Eigen::Upper> MatrixU;
319 static inline MatrixL getL(
const MatrixType& m) {
return m; }
320 static inline MatrixU getU(
const MatrixType& m) {
return m.adjoint(); }
323 template<
typename _MatrixType,
int _UpLo>
struct traits<SimplicialLDLT<_MatrixType,_UpLo> >
325 typedef _MatrixType MatrixType;
326 enum { UpLo = _UpLo };
327 typedef typename MatrixType::Scalar Scalar;
328 typedef typename MatrixType::Index Index;
329 typedef SparseMatrix<Scalar, ColMajor, Index> CholMatrixType;
330 typedef SparseTriangularView<CholMatrixType, Eigen::UnitLower> MatrixL;
331 typedef SparseTriangularView<typename CholMatrixType::AdjointReturnType, Eigen::UnitUpper> MatrixU;
332 static inline MatrixL getL(
const MatrixType& m) {
return m; }
333 static inline MatrixU getU(
const MatrixType& m) {
return m.adjoint(); }
336 template<
typename _MatrixType,
int _UpLo>
struct traits<SimplicialCholesky<_MatrixType,_UpLo> >
338 typedef _MatrixType MatrixType;
339 enum { UpLo = _UpLo };
361 template<
typename _MatrixType,
int _UpLo>
365 typedef _MatrixType MatrixType;
366 enum { UpLo = _UpLo };
368 typedef typename MatrixType::Scalar Scalar;
369 typedef typename MatrixType::RealScalar RealScalar;
370 typedef typename MatrixType::Index Index;
373 typedef internal::traits<SimplicialLLT> Traits;
374 typedef typename Traits::MatrixL MatrixL;
375 typedef typename Traits::MatrixU MatrixU;
385 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial LLT not factorized");
386 return Traits::getL(Base::m_matrix);
391 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial LLT not factorized");
392 return Traits::getU(Base::m_matrix);
398 Base::template compute<false>(matrix);
410 Base::analyzePattern(a,
false);
421 Base::template factorize<false>(a);
427 Scalar detL = Base::m_matrix.diagonal().prod();
428 return internal::abs2(detL);
449 template<
typename _MatrixType,
int _UpLo>
453 typedef _MatrixType MatrixType;
454 enum { UpLo = _UpLo };
456 typedef typename MatrixType::Scalar Scalar;
457 typedef typename MatrixType::RealScalar RealScalar;
458 typedef typename MatrixType::Index Index;
461 typedef internal::traits<SimplicialLDLT> Traits;
462 typedef typename Traits::MatrixL MatrixL;
463 typedef typename Traits::MatrixU MatrixU;
474 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial LDLT not factorized");
479 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial LDLT not factorized");
480 return Traits::getL(Base::m_matrix);
485 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial LDLT not factorized");
486 return Traits::getU(Base::m_matrix);
492 Base::template compute<true>(matrix);
504 Base::analyzePattern(a,
true);
515 Base::template factorize<true>(a);
521 return Base::m_diag.prod();
531 template<
typename _MatrixType,
int _UpLo>
535 typedef _MatrixType MatrixType;
536 enum { UpLo = _UpLo };
538 typedef typename MatrixType::Scalar Scalar;
539 typedef typename MatrixType::RealScalar RealScalar;
540 typedef typename MatrixType::Index Index;
543 typedef internal::traits<SimplicialCholesky> Traits;
544 typedef internal::traits<SimplicialLDLT<MatrixType,UpLo> > LDLTTraits;
545 typedef internal::traits<SimplicialLLT<MatrixType,UpLo> > LLTTraits;
550 :
Base(), m_LDLT(
true)
559 case SimplicialCholeskyLLT:
562 case SimplicialCholeskyLDLT:
573 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial Cholesky not factorized");
577 eigen_assert(Base::m_factorizationIsOk &&
"Simplicial Cholesky not factorized");
578 return Base::m_matrix;
585 Base::template compute<true>(matrix);
587 Base::template compute<false>(matrix);
599 Base::analyzePattern(a, m_LDLT);
611 Base::template factorize<true>(a);
613 Base::template factorize<false>(a);
617 template<
typename Rhs,
typename Dest>
620 eigen_assert(Base::m_factorizationIsOk &&
"The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
621 eigen_assert(Base::m_matrix.rows()==b.rows());
626 if(Base::m_P.size()>0)
627 dest = Base::m_P * b;
631 if(Base::m_matrix.nonZeros()>0)
634 LDLTTraits::getL(Base::m_matrix).solveInPlace(dest);
636 LLTTraits::getL(Base::m_matrix).solveInPlace(dest);
639 if(Base::m_diag.size()>0)
640 dest = Base::m_diag.
asDiagonal().inverse() * dest;
642 if (Base::m_matrix.nonZeros()>0)
645 LDLTTraits::getU(Base::m_matrix).solveInPlace(dest);
647 LLTTraits::getU(Base::m_matrix).solveInPlace(dest);
650 if(Base::m_P.size()>0)
651 dest = Base::m_Pinv * dest;
654 Scalar determinant()
const
658 return Base::m_diag.
prod();
662 Scalar detL = Diagonal<const CholMatrixType>(Base::m_matrix).prod();
663 return internal::abs2(detL);
671 template<
typename Derived>
672 void SimplicialCholeskyBase<Derived>::ordering(
const MatrixType& a, CholMatrixType& ap)
674 eigen_assert(a.rows()==a.cols());
675 const Index size = a.rows();
680 C = a.template selfadjointView<UpLo>();
684 internal::minimum_degree_ordering(C, m_Pinv);
688 m_P = m_Pinv.inverse();
692 ap.resize(size,size);
693 ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
696 template<
typename Derived>
697 void SimplicialCholeskyBase<Derived>::analyzePattern_preordered(
const CholMatrixType& ap,
bool doLDLT)
699 const Index size = ap.rows();
700 m_matrix.resize(size, size);
701 m_parent.resize(size);
702 m_nonZerosPerCol.resize(size);
704 ei_declare_aligned_stack_constructed_variable(Index, tags, size, 0);
706 for(Index k = 0; k < size; ++k)
711 m_nonZerosPerCol[k] = 0;
712 for(
typename CholMatrixType::InnerIterator it(ap,k); it; ++it)
714 Index i = it.index();
718 for(; tags[i] != k; i = m_parent[i])
721 if (m_parent[i] == -1)
723 m_nonZerosPerCol[i]++;
731 Index* Lp = m_matrix.outerIndexPtr();
733 for(Index k = 0; k < size; ++k)
734 Lp[k+1] = Lp[k] + m_nonZerosPerCol[k] + (doLDLT ? 0 : 1);
736 m_matrix.resizeNonZeros(Lp[size]);
738 m_isInitialized =
true;
740 m_analysisIsOk =
true;
741 m_factorizationIsOk =
false;
745 template<
typename Derived>
746 template<
bool DoLDLT>
747 void SimplicialCholeskyBase<Derived>::factorize_preordered(
const CholMatrixType& ap)
749 eigen_assert(m_analysisIsOk &&
"You must first call analyzePattern()");
750 eigen_assert(ap.rows()==ap.cols());
751 const Index size = ap.rows();
752 eigen_assert(m_parent.size()==size);
753 eigen_assert(m_nonZerosPerCol.size()==size);
755 const Index* Lp = m_matrix.outerIndexPtr();
756 Index* Li = m_matrix.innerIndexPtr();
757 Scalar* Lx = m_matrix.valuePtr();
759 ei_declare_aligned_stack_constructed_variable(Scalar, y, size, 0);
760 ei_declare_aligned_stack_constructed_variable(Index, pattern, size, 0);
761 ei_declare_aligned_stack_constructed_variable(Index, tags, size, 0);
764 m_diag.resize(DoLDLT ? size : 0);
766 for(Index k = 0; k < size; ++k)
772 m_nonZerosPerCol[k] = 0;
773 for(
typename MatrixType::InnerIterator it(ap,k); it; ++it)
775 Index i = it.index();
778 y[i] += internal::conj(it.value());
780 for(len = 0; tags[i] != k; i = m_parent[i])
786 pattern[--top] = pattern[--len];
792 RealScalar d = internal::real(y[k]) * m_shiftScale + m_shiftOffset;
794 for(; top < size; ++top)
796 Index i = pattern[top];
803 l_ki = yi / m_diag[i];
805 yi = l_ki = yi / Lx[Lp[i]];
807 Index p2 = Lp[i] + m_nonZerosPerCol[i];
809 for(p = Lp[i] + (DoLDLT ? 0 : 1); p < p2; ++p)
810 y[Li[p]] -= internal::conj(Lx[p]) * yi;
811 d -= internal::real(l_ki * internal::conj(yi));
814 ++m_nonZerosPerCol[i];
819 if(d == RealScalar(0))
827 Index p = Lp[k] + m_nonZerosPerCol[k]++;
829 if(d <= RealScalar(0)) {
833 Lx[p] = internal::sqrt(d) ;
838 m_factorizationIsOk =
true;
843 template<
typename Derived,
typename Rhs>
844 struct solve_retval<SimplicialCholeskyBase<Derived>, Rhs>
845 : solve_retval_base<SimplicialCholeskyBase<Derived>, Rhs>
847 typedef SimplicialCholeskyBase<Derived> Dec;
848 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
850 template<typename Dest>
void evalTo(Dest& dst)
const
852 dec().derived()._solve(rhs(),dst);
856 template<
typename Derived,
typename Rhs>
857 struct sparse_solve_retval<SimplicialCholeskyBase<Derived>, Rhs>
858 : sparse_solve_retval_base<SimplicialCholeskyBase<Derived>, Rhs>
860 typedef SimplicialCholeskyBase<Derived> Dec;
861 EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
863 template<typename Dest>
void evalTo(Dest& dst)
const
865 dec().derived()._solve_sparse(rhs(),dst);
873 #endif // EIGEN_SIMPLICIAL_CHOLESKY_H