63 #ifndef EIGEN_SIMPLICIAL_CHOLESKY_H
64 #define EIGEN_SIMPLICIAL_CHOLESKY_H
85 template<
typename Derived>
89 typedef typename internal::traits<Derived>::MatrixType
MatrixType;
90 enum {
UpLo = internal::traits<Derived>::UpLo };
91 typedef typename MatrixType::Scalar
Scalar;
93 typedef typename MatrixType::Index
Index;
114 Derived&
derived() {
return *
static_cast<Derived*
>(
this); }
115 const Derived&
derived()
const {
return *
static_cast<const Derived*
>(
this); }
135 template<
typename Rhs>
136 inline const internal::solve_retval<SimplicialCholeskyBase, Rhs>
141 &&
"SimplicialCholeskyBase::solve(): invalid number of rows of the right hand side matrix b");
142 return internal::solve_retval<SimplicialCholeskyBase, Rhs>(*
this, b.derived());
149 template<
typename Rhs>
150 inline const internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs>
155 &&
"SimplicialCholesky::solve(): invalid number of rows of the right hand side matrix b");
156 return internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs>(*
this, b.
derived());
185 #ifndef EIGEN_PARSED_BY_DOXYGEN
187 template<
typename Stream>
188 void dumpMemory(Stream& s)
192 s <<
" diag: " << ((total+=
m_diag.size() *
sizeof(
Scalar)) >> 20) <<
"Mb" <<
"\n";
193 s <<
" tree: " << ((total+=
m_parent.size() *
sizeof(
int)) >> 20) <<
"Mb" <<
"\n";
194 s <<
" nonzeros: " << ((total+=
m_nonZerosPerCol.size() *
sizeof(
int)) >> 20) <<
"Mb" <<
"\n";
195 s <<
" perm: " << ((total+=
m_P.
size() *
sizeof(
int)) >> 20) <<
"Mb" <<
"\n";
196 s <<
" perm^-1: " << ((total+=
m_Pinv.
size() *
sizeof(
int)) >> 20) <<
"Mb" <<
"\n";
197 s <<
" TOTAL: " << (total>> 20) <<
"Mb" <<
"\n";
201 template<
typename Rhs,
typename Dest>
202 void _solve(
const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest)
const
204 eigen_assert(
m_factorizationIsOk &&
"The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
216 derived().matrixL().solveInPlace(dest);
219 dest =
m_diag.asDiagonal().inverse() * dest;
222 derived().matrixU().solveInPlace(dest);
229 template<
typename Rhs,
typename DestScalar,
int DestOptions,
typename DestIndex>
230 void _solve_sparse(
const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest)
const
232 eigen_assert(
m_factorizationIsOk &&
"The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
236 static const int NbColsAtOnce = 4;
237 int rhsCols = b.cols();
240 for(
int k=0; k<rhsCols; k+=NbColsAtOnce)
242 int actualCols = std::min<int>(rhsCols-k, NbColsAtOnce);
243 tmp.leftCols(actualCols) = b.middleCols(k,actualCols);
244 tmp.leftCols(actualCols) =
derived().solve(tmp.leftCols(actualCols));
245 dest.middleCols(k,actualCols) = tmp.leftCols(actualCols).sparseView();
249 #endif // EIGEN_PARSED_BY_DOXYGEN
254 template<
bool DoLDLT>
258 Index size = matrix.cols();
262 factorize_preordered<DoLDLT>(ap);
265 template<
bool DoLDLT>
271 ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().
twistedBy(
m_Pinv);
272 factorize_preordered<DoLDLT>(ap);
275 template<
bool DoLDLT>
314 template<
typename _MatrixType,
int _UpLo = Lower>
class SimplicialLLT;
315 template<
typename _MatrixType,
int _UpLo = Lower>
class SimplicialLDLT;
322 typedef _MatrixType MatrixType;
323 enum { UpLo = _UpLo };
324 typedef typename MatrixType::Scalar Scalar;
325 typedef typename MatrixType::Index Index;
326 typedef SparseMatrix<Scalar, ColMajor, Index> CholMatrixType;
327 typedef SparseTriangularView<CholMatrixType, Eigen::Lower> MatrixL;
328 typedef SparseTriangularView<typename CholMatrixType::AdjointReturnType, Eigen::Upper> MatrixU;
329 static inline MatrixL getL(
const MatrixType& m) {
return m; }
330 static inline MatrixU getU(
const MatrixType& m) {
return m.adjoint(); }
333 template<
typename _MatrixType,
int _UpLo>
struct traits<SimplicialLDLT<_MatrixType,_UpLo> >
335 typedef _MatrixType MatrixType;
336 enum { UpLo = _UpLo };
337 typedef typename MatrixType::Scalar Scalar;
338 typedef typename MatrixType::Index Index;
339 typedef SparseMatrix<Scalar, ColMajor, Index> CholMatrixType;
340 typedef SparseTriangularView<CholMatrixType, Eigen::UnitLower> MatrixL;
341 typedef SparseTriangularView<typename CholMatrixType::AdjointReturnType, Eigen::UnitUpper> MatrixU;
342 static inline MatrixL getL(
const MatrixType& m) {
return m; }
343 static inline MatrixU getU(
const MatrixType& m) {
return m.adjoint(); }
346 template<
typename _MatrixType,
int _UpLo>
struct traits<SimplicialCholesky<_MatrixType,_UpLo> >
348 typedef _MatrixType MatrixType;
349 enum { UpLo = _UpLo };
368 template<
typename _MatrixType,
int _UpLo>
375 typedef typename MatrixType::Scalar
Scalar;
377 typedef typename MatrixType::Index
Index;
380 typedef internal::traits<SimplicialLLT>
Traits;
405 Base::template compute<false>(matrix);
428 Base::template factorize<false>(a);
453 template<
typename _MatrixType,
int _UpLo>
460 typedef typename MatrixType::Scalar
Scalar;
462 typedef typename MatrixType::Index
Index;
465 typedef internal::traits<SimplicialLDLT>
Traits;
496 Base::template compute<true>(matrix);
519 Base::template factorize<true>(a);
535 template<
typename _MatrixType,
int _UpLo>
542 typedef typename MatrixType::Scalar
Scalar;
544 typedef typename MatrixType::Index
Index;
547 typedef internal::traits<SimplicialCholesky>
Traits;
548 typedef internal::traits<SimplicialLDLT<MatrixType,UpLo> >
LDLTTraits;
549 typedef internal::traits<SimplicialLLT<MatrixType,UpLo> >
LLTTraits;
589 Base::template compute<true>(matrix);
591 Base::template compute<false>(matrix);
615 Base::template factorize<true>(a);
617 Base::template factorize<false>(a);
621 template<
typename Rhs,
typename Dest>
675 template<
typename Derived>
679 const Index size = a.rows();
683 C = a.template selfadjointView<UpLo>();
691 m_Pinv = m_P.inverse();
696 ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().
twistedBy(m_Pinv);
699 template<
typename Derived>
703 m_matrix.resize(size, size);
704 m_parent.resize(size);
705 m_nonZerosPerCol.resize(size);
709 for(
Index k = 0; k < size; ++k)
714 m_nonZerosPerCol[k] = 0;
715 for(
typename CholMatrixType::InnerIterator it(ap,k); it; ++it)
717 Index i = it.index();
721 for(; tags[i] != k; i = m_parent[i])
724 if (m_parent[i] == -1)
726 m_nonZerosPerCol[i]++;
734 Index* Lp = m_matrix.outerIndexPtr();
736 for(
Index k = 0; k < size; ++k)
737 Lp[k+1] = Lp[k] + m_nonZerosPerCol[k] + (doLDLT ? 0 : 1);
739 m_matrix.resizeNonZeros(Lp[size]);
741 m_isInitialized =
true;
743 m_analysisIsOk =
true;
744 m_factorizationIsOk =
false;
748 template<
typename Derived>
749 template<
bool DoLDLT>
752 eigen_assert(m_analysisIsOk &&
"You must first call analyzePattern()");
758 const Index* Lp = m_matrix.outerIndexPtr();
759 Index* Li = m_matrix.innerIndexPtr();
760 Scalar* Lx = m_matrix.valuePtr();
767 m_diag.resize(DoLDLT ? size : 0);
769 for(
Index k = 0; k < size; ++k)
775 m_nonZerosPerCol[k] = 0;
776 for(
typename MatrixType::InnerIterator it(ap,k); it; ++it)
778 Index i = it.index();
783 for(len = 0; tags[i] != k; i = m_parent[i])
789 pattern[--top] = pattern[--len];
797 for(; top < size; ++top)
799 Index i = pattern[top];
806 l_ki = yi / m_diag[i];
808 yi = l_ki = yi / Lx[Lp[i]];
810 Index p2 = Lp[i] + m_nonZerosPerCol[i];
812 for(p = Lp[i] + (DoLDLT ? 0 : 1); p < p2; ++
p)
817 ++m_nonZerosPerCol[i];
830 Index p = Lp[k] + m_nonZerosPerCol[k]++;
841 m_factorizationIsOk =
true;
846 template<
typename Derived,
typename Rhs>
848 : solve_retval_base<SimplicialCholeskyBase<Derived>, Rhs>
853 template<typename Dest>
void evalTo(Dest& dst)
const
855 dec().derived()._solve(rhs(),dst);
859 template<
typename Derived,
typename Rhs>
860 struct sparse_solve_retval<SimplicialCholeskyBase<Derived>, Rhs>
861 : sparse_solve_retval_base<SimplicialCholeskyBase<Derived>, Rhs>
863 typedef SimplicialCholeskyBase<Derived> Dec;
866 template<typename Dest>
void evalTo(Dest& dst)
const
868 dec().derived()._solve_sparse(rhs(),dst);
876 #endif // EIGEN_SIMPLICIAL_CHOLESKY_H