63 #ifndef EIGEN_SIMPLICIAL_CHOLESKY_H
64 #define EIGEN_SIMPLICIAL_CHOLESKY_H
88 template<
typename Derived>
92 typedef typename internal::traits<Derived>::MatrixType
MatrixType;
93 enum {
UpLo = internal::traits<Derived>::UpLo };
94 typedef typename MatrixType::Scalar
Scalar;
96 typedef typename MatrixType::Index
Index;
117 Derived&
derived() {
return *
static_cast<Derived*
>(
this); }
118 const Derived&
derived()
const {
return *
static_cast<const Derived*
>(
this); }
138 template<
typename Rhs>
139 inline const internal::solve_retval<SimplicialCholeskyBase, Rhs>
144 &&
"SimplicialCholeskyBase::solve(): invalid number of rows of the right hand side matrix b");
145 return internal::solve_retval<SimplicialCholeskyBase, Rhs>(*
this, b.derived());
152 template<
typename Rhs>
153 inline const internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs>
158 &&
"SimplicialCholesky::solve(): invalid number of rows of the right hand side matrix b");
159 return internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs>(*
this, b.
derived());
188 #ifndef EIGEN_PARSED_BY_DOXYGEN
190 template<
typename Stream>
191 void dumpMemory(Stream& s)
195 s <<
" diag: " << ((total+=
m_diag.size() *
sizeof(
Scalar)) >> 20) <<
"Mb" <<
"\n";
196 s <<
" tree: " << ((total+=
m_parent.size() *
sizeof(
int)) >> 20) <<
"Mb" <<
"\n";
197 s <<
" nonzeros: " << ((total+=
m_nonZerosPerCol.size() *
sizeof(
int)) >> 20) <<
"Mb" <<
"\n";
198 s <<
" perm: " << ((total+=
m_P.
size() *
sizeof(
int)) >> 20) <<
"Mb" <<
"\n";
199 s <<
" perm^-1: " << ((total+=
m_Pinv.
size() *
sizeof(
int)) >> 20) <<
"Mb" <<
"\n";
200 s <<
" TOTAL: " << (total>> 20) <<
"Mb" <<
"\n";
204 template<
typename Rhs,
typename Dest>
205 void _solve(
const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest)
const
207 eigen_assert(
m_factorizationIsOk &&
"The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
219 derived().matrixL().solveInPlace(dest);
222 dest =
m_diag.asDiagonal().inverse() * dest;
225 derived().matrixU().solveInPlace(dest);
232 template<
typename Rhs,
typename DestScalar,
int DestOptions,
typename DestIndex>
233 void _solve_sparse(
const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest)
const
235 eigen_assert(
m_factorizationIsOk &&
"The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
239 static const int NbColsAtOnce = 4;
240 int rhsCols = b.cols();
243 for(
int k=0; k<rhsCols; k+=NbColsAtOnce)
245 int actualCols = std::min<int>(rhsCols-k, NbColsAtOnce);
246 tmp.leftCols(actualCols) = b.middleCols(k,actualCols);
247 tmp.leftCols(actualCols) =
derived().solve(tmp.leftCols(actualCols));
248 dest.middleCols(k,actualCols) = tmp.leftCols(actualCols).sparseView();
252 #endif // EIGEN_PARSED_BY_DOXYGEN
257 template<
bool DoLDLT>
261 Index size = matrix.cols();
265 factorize_preordered<DoLDLT>(ap);
268 template<
bool DoLDLT>
274 ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().
twistedBy(
m_P);
275 factorize_preordered<DoLDLT>(ap);
278 template<
bool DoLDLT>
317 template<
typename _MatrixType,
int _UpLo = Lower>
class SimplicialLLT;
318 template<
typename _MatrixType,
int _UpLo = Lower>
class SimplicialLDLT;
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::Lower> MatrixL;
331 typedef SparseTriangularView<typename CholMatrixType::AdjointReturnType, Eigen::Upper> 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<SimplicialLDLT<_MatrixType,_UpLo> >
338 typedef _MatrixType MatrixType;
339 enum { UpLo = _UpLo };
340 typedef typename MatrixType::Scalar Scalar;
341 typedef typename MatrixType::Index Index;
342 typedef SparseMatrix<Scalar, ColMajor, Index> CholMatrixType;
343 typedef SparseTriangularView<CholMatrixType, Eigen::UnitLower> MatrixL;
344 typedef SparseTriangularView<typename CholMatrixType::AdjointReturnType, Eigen::UnitUpper> MatrixU;
345 static inline MatrixL getL(
const MatrixType& m) {
return m; }
346 static inline MatrixU getU(
const MatrixType& m) {
return m.adjoint(); }
349 template<
typename _MatrixType,
int _UpLo>
struct traits<SimplicialCholesky<_MatrixType,_UpLo> >
351 typedef _MatrixType MatrixType;
352 enum { UpLo = _UpLo };
374 template<
typename _MatrixType,
int _UpLo>
381 typedef typename MatrixType::Scalar
Scalar;
383 typedef typename MatrixType::Index
Index;
386 typedef internal::traits<SimplicialLLT>
Traits;
411 Base::template compute<false>(matrix);
434 Base::template factorize<false>(a);
462 template<
typename _MatrixType,
int _UpLo>
469 typedef typename MatrixType::Scalar
Scalar;
471 typedef typename MatrixType::Index
Index;
474 typedef internal::traits<SimplicialLDLT>
Traits;
505 Base::template compute<true>(matrix);
528 Base::template factorize<true>(a);
544 template<
typename _MatrixType,
int _UpLo>
551 typedef typename MatrixType::Scalar
Scalar;
553 typedef typename MatrixType::Index
Index;
556 typedef internal::traits<SimplicialCholesky>
Traits;
557 typedef internal::traits<SimplicialLDLT<MatrixType,UpLo> >
LDLTTraits;
558 typedef internal::traits<SimplicialLLT<MatrixType,UpLo> >
LLTTraits;
598 Base::template compute<true>(matrix);
600 Base::template compute<false>(matrix);
624 Base::template factorize<true>(a);
626 Base::template factorize<false>(a);
630 template<
typename Rhs,
typename Dest>
684 template<
typename Derived>
688 const Index size = a.rows();
693 C = a.template selfadjointView<UpLo>();
701 m_P = m_Pinv.inverse();
706 ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().
twistedBy(m_P);
709 template<
typename Derived>
713 m_matrix.resize(size, size);
714 m_parent.resize(size);
715 m_nonZerosPerCol.resize(size);
719 for(
Index k = 0; k < size; ++k)
724 m_nonZerosPerCol[k] = 0;
725 for(
typename CholMatrixType::InnerIterator it(ap,k); it; ++it)
727 Index i = it.index();
731 for(; tags[i] != k; i = m_parent[i])
734 if (m_parent[i] == -1)
736 m_nonZerosPerCol[i]++;
744 Index* Lp = m_matrix.outerIndexPtr();
746 for(
Index k = 0; k < size; ++k)
747 Lp[k+1] = Lp[k] + m_nonZerosPerCol[k] + (doLDLT ? 0 : 1);
749 m_matrix.resizeNonZeros(Lp[size]);
751 m_isInitialized =
true;
753 m_analysisIsOk =
true;
754 m_factorizationIsOk =
false;
758 template<
typename Derived>
759 template<
bool DoLDLT>
762 eigen_assert(m_analysisIsOk &&
"You must first call analyzePattern()");
768 const Index* Lp = m_matrix.outerIndexPtr();
769 Index* Li = m_matrix.innerIndexPtr();
770 Scalar* Lx = m_matrix.valuePtr();
777 m_diag.resize(DoLDLT ? size : 0);
779 for(
Index k = 0; k < size; ++k)
785 m_nonZerosPerCol[k] = 0;
786 for(
typename MatrixType::InnerIterator it(ap,k); it; ++it)
788 Index i = it.index();
793 for(len = 0; tags[i] != k; i = m_parent[i])
799 pattern[--top] = pattern[--len];
807 for(; top < size; ++top)
809 Index i = pattern[top];
816 l_ki = yi / m_diag[i];
818 yi = l_ki = yi / Lx[Lp[i]];
820 Index p2 = Lp[i] + m_nonZerosPerCol[i];
822 for(p = Lp[i] + (DoLDLT ? 0 : 1); p < p2; ++
p)
827 ++m_nonZerosPerCol[i];
840 Index p = Lp[k] + m_nonZerosPerCol[k]++;
851 m_factorizationIsOk =
true;
856 template<
typename Derived,
typename Rhs>
858 : solve_retval_base<SimplicialCholeskyBase<Derived>, Rhs>
863 template<typename Dest>
void evalTo(Dest& dst)
const
865 dec().derived()._solve(rhs(),dst);
869 template<
typename Derived,
typename Rhs>
870 struct sparse_solve_retval<SimplicialCholeskyBase<Derived>, Rhs>
871 : sparse_solve_retval_base<SimplicialCholeskyBase<Derived>, Rhs>
873 typedef SimplicialCholeskyBase<Derived> Dec;
876 template<typename Dest>
void evalTo(Dest& dst)
const
878 dec().derived()._solve_sparse(rhs(),dst);
886 #endif // EIGEN_SIMPLICIAL_CHOLESKY_H