MOAB
4.9.3pre
|
A base class for direct sparse Cholesky factorizations. More...
#include <SimplicialCholesky.h>
Classes | |
struct | keep_diag |
Public Types | |
enum | { UpLo = internal::traits<Derived>::UpLo } |
enum | { ColsAtCompileTime = MatrixType::ColsAtCompileTime, MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime } |
typedef internal::traits < Derived >::MatrixType | MatrixType |
typedef internal::traits < Derived >::OrderingType | OrderingType |
typedef MatrixType::Scalar | Scalar |
typedef MatrixType::RealScalar | RealScalar |
typedef MatrixType::StorageIndex | StorageIndex |
typedef SparseMatrix< Scalar, ColMajor, StorageIndex > | CholMatrixType |
typedef CholMatrixType const * | ConstCholMatrixPtr |
typedef Matrix< Scalar, Dynamic, 1 > | VectorType |
typedef Matrix< StorageIndex, Dynamic, 1 > | VectorI |
Public Member Functions | |
SimplicialCholeskyBase () | |
SimplicialCholeskyBase (const MatrixType &matrix) | |
~SimplicialCholeskyBase () | |
Derived & | derived () |
const Derived & | derived () const |
Index | cols () const |
Index | rows () const |
ComputationInfo | info () const |
Reports whether previous computation was successful. | |
const PermutationMatrix < Dynamic, Dynamic, StorageIndex > & | permutationP () const |
const PermutationMatrix < Dynamic, Dynamic, StorageIndex > & | permutationPinv () const |
Derived & | setShift (const RealScalar &offset, const RealScalar &scale=1) |
template<typename Stream > | |
void | dumpMemory (Stream &s) |
template<typename Rhs , typename Dest > | |
void | _solve_impl (const MatrixBase< Rhs > &b, MatrixBase< Dest > &dest) const |
template<typename Rhs , typename Dest > | |
void | _solve_impl (const SparseMatrixBase< Rhs > &b, SparseMatrixBase< Dest > &dest) const |
Protected Member Functions | |
template<bool DoLDLT> | |
void | compute (const MatrixType &matrix) |
template<bool DoLDLT> | |
void | factorize (const MatrixType &a) |
template<bool DoLDLT> | |
void | factorize_preordered (const CholMatrixType &a) |
void | analyzePattern (const MatrixType &a, bool doLDLT) |
void | analyzePattern_preordered (const CholMatrixType &a, bool doLDLT) |
void | ordering (const MatrixType &a, ConstCholMatrixPtr &pmat, CholMatrixType &ap) |
Protected Attributes | |
ComputationInfo | m_info |
bool | m_factorizationIsOk |
bool | m_analysisIsOk |
CholMatrixType | m_matrix |
VectorType | m_diag |
VectorI | m_parent |
VectorI | m_nonZerosPerCol |
PermutationMatrix< Dynamic, Dynamic, StorageIndex > | m_P |
PermutationMatrix< Dynamic, Dynamic, StorageIndex > | m_Pinv |
RealScalar | m_shiftOffset |
RealScalar | m_shiftScale |
Private Types | |
typedef SparseSolverBase< Derived > | Base |
A base class for direct sparse Cholesky factorizations.
This is a base class for LL^T and LDL^T Cholesky factorizations of sparse matrices that are selfadjoint and positive definite. These factorizations allow for solving A.X = B where X and B can be either dense or sparse.
In order to reduce the fill-in, a symmetric permutation P is applied prior to the factorization such that the factorized matrix is P A P^-1.
Derived | the type of the derived class, that is the actual factorization type. |
Definition at line 55 of file SimplicialCholesky.h.
typedef SparseSolverBase<Derived> Eigen::SimplicialCholeskyBase< Derived >::Base [private] |
Reimplemented in Eigen::SimplicialCholesky< _MatrixType, _UpLo, _Ordering >, Eigen::SimplicialLDLT< _MatrixType, _UpLo, _Ordering >, and Eigen::SimplicialLLT< _MatrixType, _UpLo, _Ordering >.
Definition at line 57 of file SimplicialCholesky.h.
typedef SparseMatrix<Scalar,ColMajor,StorageIndex> Eigen::SimplicialCholeskyBase< Derived >::CholMatrixType |
Reimplemented in Eigen::SimplicialCholesky< _MatrixType, _UpLo, _Ordering >, Eigen::SimplicialLDLT< _MatrixType, _UpLo, _Ordering >, and Eigen::SimplicialLLT< _MatrixType, _UpLo, _Ordering >.
Definition at line 67 of file SimplicialCholesky.h.
typedef CholMatrixType const* Eigen::SimplicialCholeskyBase< Derived >::ConstCholMatrixPtr |
Definition at line 68 of file SimplicialCholesky.h.
typedef internal::traits<Derived>::MatrixType Eigen::SimplicialCholeskyBase< Derived >::MatrixType |
Reimplemented in Eigen::SimplicialCholesky< _MatrixType, _UpLo, _Ordering >, Eigen::SimplicialLDLT< _MatrixType, _UpLo, _Ordering >, and Eigen::SimplicialLLT< _MatrixType, _UpLo, _Ordering >.
Definition at line 61 of file SimplicialCholesky.h.
typedef internal::traits<Derived>::OrderingType Eigen::SimplicialCholeskyBase< Derived >::OrderingType |
Definition at line 62 of file SimplicialCholesky.h.
typedef MatrixType::RealScalar Eigen::SimplicialCholeskyBase< Derived >::RealScalar |
Reimplemented in Eigen::SimplicialCholesky< _MatrixType, _UpLo, _Ordering >, Eigen::SimplicialLDLT< _MatrixType, _UpLo, _Ordering >, and Eigen::SimplicialLLT< _MatrixType, _UpLo, _Ordering >.
Definition at line 65 of file SimplicialCholesky.h.
typedef MatrixType::Scalar Eigen::SimplicialCholeskyBase< Derived >::Scalar |
Reimplemented in Eigen::SimplicialCholesky< _MatrixType, _UpLo, _Ordering >, Eigen::SimplicialLDLT< _MatrixType, _UpLo, _Ordering >, and Eigen::SimplicialLLT< _MatrixType, _UpLo, _Ordering >.
Definition at line 64 of file SimplicialCholesky.h.
typedef MatrixType::StorageIndex Eigen::SimplicialCholeskyBase< Derived >::StorageIndex |
Reimplemented in Eigen::SimplicialCholesky< _MatrixType, _UpLo, _Ordering >, Eigen::SimplicialLDLT< _MatrixType, _UpLo, _Ordering >, and Eigen::SimplicialLLT< _MatrixType, _UpLo, _Ordering >.
Definition at line 66 of file SimplicialCholesky.h.
typedef Matrix<StorageIndex,Dynamic,1> Eigen::SimplicialCholeskyBase< Derived >::VectorI |
Definition at line 70 of file SimplicialCholesky.h.
typedef Matrix<Scalar,Dynamic,1> Eigen::SimplicialCholeskyBase< Derived >::VectorType |
Reimplemented in Eigen::SimplicialCholesky< _MatrixType, _UpLo, _Ordering >, Eigen::SimplicialLDLT< _MatrixType, _UpLo, _Ordering >, and Eigen::SimplicialLLT< _MatrixType, _UpLo, _Ordering >.
Definition at line 69 of file SimplicialCholesky.h.
anonymous enum |
Definition at line 63 of file SimplicialCholesky.h.
{ UpLo = internal::traits<Derived>::UpLo };
anonymous enum |
Definition at line 72 of file SimplicialCholesky.h.
{ ColsAtCompileTime = MatrixType::ColsAtCompileTime, MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime };
Eigen::SimplicialCholeskyBase< Derived >::SimplicialCholeskyBase | ( | ) | [inline] |
Default constructor
Definition at line 82 of file SimplicialCholesky.h.
: m_info(Success), m_shiftOffset(0), m_shiftScale(1) {}
Eigen::SimplicialCholeskyBase< Derived >::SimplicialCholeskyBase | ( | const MatrixType & | matrix | ) | [inline, explicit] |
Definition at line 86 of file SimplicialCholesky.h.
: m_info(Success), m_shiftOffset(0), m_shiftScale(1) { derived().compute(matrix); }
Eigen::SimplicialCholeskyBase< Derived >::~SimplicialCholeskyBase | ( | ) | [inline] |
Definition at line 92 of file SimplicialCholesky.h.
{ }
void Eigen::SimplicialCholeskyBase< Derived >::_solve_impl | ( | const MatrixBase< Rhs > & | b, |
MatrixBase< Dest > & | dest | ||
) | const [inline] |
Reimplemented in Eigen::SimplicialCholesky< _MatrixType, _UpLo, _Ordering >.
Definition at line 156 of file SimplicialCholesky.h.
{ eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()"); eigen_assert(m_matrix.rows()==b.rows()); if(m_info!=Success) return; if(m_P.size()>0) dest = m_P * b; else dest = b; if(m_matrix.nonZeros()>0) // otherwise L==I derived().matrixL().solveInPlace(dest); if(m_diag.size()>0) dest = m_diag.asDiagonal().inverse() * dest; if (m_matrix.nonZeros()>0) // otherwise U==I derived().matrixU().solveInPlace(dest); if(m_P.size()>0) dest = m_Pinv * dest; }
void Eigen::SimplicialCholeskyBase< Derived >::_solve_impl | ( | const SparseMatrixBase< Rhs > & | b, |
SparseMatrixBase< Dest > & | dest | ||
) | const [inline] |
default implementation of solving with a sparse rhs
Reimplemented from Eigen::SparseSolverBase< Derived >.
Reimplemented in Eigen::SimplicialCholesky< _MatrixType, _UpLo, _Ordering >.
Definition at line 183 of file SimplicialCholesky.h.
{ internal::solve_sparse_through_dense_panels(derived(), b, dest); }
void Eigen::SimplicialCholeskyBase< Derived >::analyzePattern | ( | const MatrixType & | a, |
bool | doLDLT | ||
) | [inline, protected] |
Definition at line 230 of file SimplicialCholesky.h.
{ eigen_assert(a.rows()==a.cols()); Index size = a.cols(); CholMatrixType tmp(size,size); ConstCholMatrixPtr pmat; ordering(a, pmat, tmp); analyzePattern_preordered(*pmat,doLDLT); }
void Eigen::SimplicialCholeskyBase< Derived >::analyzePattern_preordered | ( | const CholMatrixType & | a, |
bool | doLDLT | ||
) | [protected] |
Definition at line 51 of file SimplicialCholesky_impl.h.
{ const StorageIndex size = StorageIndex(ap.rows()); m_matrix.resize(size, size); m_parent.resize(size); m_nonZerosPerCol.resize(size); ei_declare_aligned_stack_constructed_variable(StorageIndex, tags, size, 0); for(StorageIndex k = 0; k < size; ++k) { /* L(k,:) pattern: all nodes reachable in etree from nz in A(0:k-1,k) */ m_parent[k] = -1; /* parent of k is not yet known */ tags[k] = k; /* mark node k as visited */ m_nonZerosPerCol[k] = 0; /* count of nonzeros in column k of L */ for(typename CholMatrixType::InnerIterator it(ap,k); it; ++it) { StorageIndex i = it.index(); if(i < k) { /* follow path from i to root of etree, stop at flagged node */ for(; tags[i] != k; i = m_parent[i]) { /* find parent of i if not yet determined */ if (m_parent[i] == -1) m_parent[i] = k; m_nonZerosPerCol[i]++; /* L (k,i) is nonzero */ tags[i] = k; /* mark i as visited */ } } } } /* construct Lp index array from m_nonZerosPerCol column counts */ StorageIndex* Lp = m_matrix.outerIndexPtr(); Lp[0] = 0; for(StorageIndex k = 0; k < size; ++k) Lp[k+1] = Lp[k] + m_nonZerosPerCol[k] + (doLDLT ? 0 : 1); m_matrix.resizeNonZeros(Lp[size]); m_isInitialized = true; m_info = Success; m_analysisIsOk = true; m_factorizationIsOk = false; }
Index Eigen::SimplicialCholeskyBase< Derived >::cols | ( | void | ) | const [inline] |
Definition at line 99 of file SimplicialCholesky.h.
void Eigen::SimplicialCholeskyBase< Derived >::compute | ( | const MatrixType & | matrix | ) | [inline, protected] |
Computes the sparse Cholesky decomposition of matrix
Reimplemented in Eigen::SimplicialCholesky< _MatrixType, _UpLo, _Ordering >, Eigen::SimplicialLDLT< _MatrixType, _UpLo, _Ordering >, and Eigen::SimplicialLLT< _MatrixType, _UpLo, _Ordering >.
Definition at line 194 of file SimplicialCholesky.h.
{ eigen_assert(matrix.rows()==matrix.cols()); Index size = matrix.cols(); CholMatrixType tmp(size,size); ConstCholMatrixPtr pmat; ordering(matrix, pmat, tmp); analyzePattern_preordered(*pmat, DoLDLT); factorize_preordered<DoLDLT>(*pmat); }
Derived& Eigen::SimplicialCholeskyBase< Derived >::derived | ( | ) | [inline] |
Reimplemented from Eigen::SparseSolverBase< Derived >.
Definition at line 96 of file SimplicialCholesky.h.
{ return *static_cast<Derived*>(this); }
const Derived& Eigen::SimplicialCholeskyBase< Derived >::derived | ( | ) | const [inline] |
Reimplemented from Eigen::SparseSolverBase< Derived >.
Definition at line 97 of file SimplicialCholesky.h.
{ return *static_cast<const Derived*>(this); }
void Eigen::SimplicialCholeskyBase< Derived >::dumpMemory | ( | Stream & | s | ) | [inline] |
Definition at line 142 of file SimplicialCholesky.h.
{ int total = 0; s << " L: " << ((total+=(m_matrix.cols()+1) * sizeof(int) + m_matrix.nonZeros()*(sizeof(int)+sizeof(Scalar))) >> 20) << "Mb" << "\n"; s << " diag: " << ((total+=m_diag.size() * sizeof(Scalar)) >> 20) << "Mb" << "\n"; s << " tree: " << ((total+=m_parent.size() * sizeof(int)) >> 20) << "Mb" << "\n"; s << " nonzeros: " << ((total+=m_nonZerosPerCol.size() * sizeof(int)) >> 20) << "Mb" << "\n"; s << " perm: " << ((total+=m_P.size() * sizeof(int)) >> 20) << "Mb" << "\n"; s << " perm^-1: " << ((total+=m_Pinv.size() * sizeof(int)) >> 20) << "Mb" << "\n"; s << " TOTAL: " << (total>> 20) << "Mb" << "\n"; }
void Eigen::SimplicialCholeskyBase< Derived >::factorize | ( | const MatrixType & | a | ) | [inline, protected] |
Reimplemented in Eigen::SimplicialCholesky< _MatrixType, _UpLo, _Ordering >, Eigen::SimplicialLDLT< _MatrixType, _UpLo, _Ordering >, and Eigen::SimplicialLLT< _MatrixType, _UpLo, _Ordering >.
Definition at line 206 of file SimplicialCholesky.h.
{ eigen_assert(a.rows()==a.cols()); Index size = a.cols(); CholMatrixType tmp(size,size); ConstCholMatrixPtr pmat; if(m_P.size()==0 && (UpLo&Upper)==Upper) { // If there is no ordering, try to directly use the input matrix without any copy internal::simplicial_cholesky_grab_input<CholMatrixType,MatrixType>::run(a, pmat, tmp); } else { tmp.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P); pmat = &tmp; } factorize_preordered<DoLDLT>(*pmat); }
void Eigen::SimplicialCholeskyBase< Derived >::factorize_preordered | ( | const CholMatrixType & | a | ) | [protected] |
Definition at line 101 of file SimplicialCholesky_impl.h.
{ using std::sqrt; eigen_assert(m_analysisIsOk && "You must first call analyzePattern()"); eigen_assert(ap.rows()==ap.cols()); eigen_assert(m_parent.size()==ap.rows()); eigen_assert(m_nonZerosPerCol.size()==ap.rows()); const StorageIndex size = StorageIndex(ap.rows()); const StorageIndex* Lp = m_matrix.outerIndexPtr(); StorageIndex* Li = m_matrix.innerIndexPtr(); Scalar* Lx = m_matrix.valuePtr(); ei_declare_aligned_stack_constructed_variable(Scalar, y, size, 0); ei_declare_aligned_stack_constructed_variable(StorageIndex, pattern, size, 0); ei_declare_aligned_stack_constructed_variable(StorageIndex, tags, size, 0); bool ok = true; m_diag.resize(DoLDLT ? size : 0); for(StorageIndex k = 0; k < size; ++k) { // compute nonzero pattern of kth row of L, in topological order y[k] = 0.0; // Y(0:k) is now all zero StorageIndex top = size; // stack for pattern is empty tags[k] = k; // mark node k as visited m_nonZerosPerCol[k] = 0; // count of nonzeros in column k of L for(typename CholMatrixType::InnerIterator it(ap,k); it; ++it) { StorageIndex i = it.index(); if(i <= k) { y[i] += numext::conj(it.value()); /* scatter A(i,k) into Y (sum duplicates) */ Index len; for(len = 0; tags[i] != k; i = m_parent[i]) { pattern[len++] = i; /* L(k,i) is nonzero */ tags[i] = k; /* mark i as visited */ } while(len > 0) pattern[--top] = pattern[--len]; } } /* compute numerical values kth row of L (a sparse triangular solve) */ RealScalar d = numext::real(y[k]) * m_shiftScale + m_shiftOffset; // get D(k,k), apply the shift function, and clear Y(k) y[k] = 0.0; for(; top < size; ++top) { Index i = pattern[top]; /* pattern[top:n-1] is pattern of L(:,k) */ Scalar yi = y[i]; /* get and clear Y(i) */ y[i] = 0.0; /* the nonzero entry L(k,i) */ Scalar l_ki; if(DoLDLT) l_ki = yi / m_diag[i]; else yi = l_ki = yi / Lx[Lp[i]]; Index p2 = Lp[i] + m_nonZerosPerCol[i]; Index p; for(p = Lp[i] + (DoLDLT ? 0 : 1); p < p2; ++p) y[Li[p]] -= numext::conj(Lx[p]) * yi; d -= numext::real(l_ki * numext::conj(yi)); Li[p] = k; /* store L(k,i) in column form of L */ Lx[p] = l_ki; ++m_nonZerosPerCol[i]; /* increment count of nonzeros in col i */ } if(DoLDLT) { m_diag[k] = d; if(d == RealScalar(0)) { ok = false; /* failure, D(k,k) is zero */ break; } } else { Index p = Lp[k] + m_nonZerosPerCol[k]++; Li[p] = k ; /* store L(k,k) = sqrt (d) in column k */ if(d <= RealScalar(0)) { ok = false; /* failure, matrix is not positive definite */ break; } Lx[p] = sqrt(d) ; } } m_info = ok ? Success : NumericalIssue; m_factorizationIsOk = true; }
ComputationInfo Eigen::SimplicialCholeskyBase< Derived >::info | ( | ) | const [inline] |
Reports whether previous computation was successful.
Success
if computation was succesful, NumericalIssue
if the matrix.appears to be negative. Definition at line 107 of file SimplicialCholesky.h.
{ eigen_assert(m_isInitialized && "Decomposition is not initialized."); return m_info; }
void Eigen::SimplicialCholeskyBase< Derived >::ordering | ( | const MatrixType & | a, |
ConstCholMatrixPtr & | pmat, | ||
CholMatrixType & | ap | ||
) | [protected] |
Definition at line 650 of file SimplicialCholesky.h.
{ eigen_assert(a.rows()==a.cols()); const Index size = a.rows(); pmat = ≈ // Note that ordering methods compute the inverse permutation if(!internal::is_same<OrderingType,NaturalOrdering<Index> >::value) { { CholMatrixType C; C = a.template selfadjointView<UpLo>(); OrderingType ordering; ordering(C,m_Pinv); } if(m_Pinv.size()>0) m_P = m_Pinv.inverse(); else m_P.resize(0); ap.resize(size,size); ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P); } else { m_Pinv.resize(0); m_P.resize(0); if(int(UpLo)==int(Lower) || MatrixType::IsRowMajor) { // we have to transpose the lower part to to the upper one ap.resize(size,size); ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>(); } else internal::simplicial_cholesky_grab_input<CholMatrixType,MatrixType>::run(a, pmat, ap); } }
const PermutationMatrix<Dynamic,Dynamic,StorageIndex>& Eigen::SimplicialCholeskyBase< Derived >::permutationP | ( | ) | const [inline] |
Definition at line 115 of file SimplicialCholesky.h.
{ return m_P; }
const PermutationMatrix<Dynamic,Dynamic,StorageIndex>& Eigen::SimplicialCholeskyBase< Derived >::permutationPinv | ( | ) | const [inline] |
Definition at line 120 of file SimplicialCholesky.h.
{ return m_Pinv; }
Index Eigen::SimplicialCholeskyBase< Derived >::rows | ( | void | ) | const [inline] |
Definition at line 100 of file SimplicialCholesky.h.
Derived& Eigen::SimplicialCholeskyBase< Derived >::setShift | ( | const RealScalar & | offset, |
const RealScalar & | scale = 1 |
||
) | [inline] |
Sets the shift parameters that will be used to adjust the diagonal coefficients during the numerical factorization.
During the numerical factorization, the diagonal coefficients are transformed by the following linear model:
d_ii
= offset + scale * d_ii
The default is the identity transformation with offset=0, and scale=1.
*this
. Definition at line 132 of file SimplicialCholesky.h.
{ m_shiftOffset = offset; m_shiftScale = scale; return derived(); }
bool Eigen::SimplicialCholeskyBase< Derived >::m_analysisIsOk [protected] |
Definition at line 253 of file SimplicialCholesky.h.
VectorType Eigen::SimplicialCholeskyBase< Derived >::m_diag [protected] |
Definition at line 256 of file SimplicialCholesky.h.
bool Eigen::SimplicialCholeskyBase< Derived >::m_factorizationIsOk [protected] |
Definition at line 252 of file SimplicialCholesky.h.
ComputationInfo Eigen::SimplicialCholeskyBase< Derived >::m_info [mutable, protected] |
Definition at line 251 of file SimplicialCholesky.h.
CholMatrixType Eigen::SimplicialCholeskyBase< Derived >::m_matrix [protected] |
Definition at line 255 of file SimplicialCholesky.h.
VectorI Eigen::SimplicialCholeskyBase< Derived >::m_nonZerosPerCol [protected] |
Definition at line 258 of file SimplicialCholesky.h.
PermutationMatrix<Dynamic,Dynamic,StorageIndex> Eigen::SimplicialCholeskyBase< Derived >::m_P [protected] |
Definition at line 259 of file SimplicialCholesky.h.
VectorI Eigen::SimplicialCholeskyBase< Derived >::m_parent [protected] |
Definition at line 257 of file SimplicialCholesky.h.
PermutationMatrix<Dynamic,Dynamic,StorageIndex> Eigen::SimplicialCholeskyBase< Derived >::m_Pinv [protected] |
Definition at line 260 of file SimplicialCholesky.h.
RealScalar Eigen::SimplicialCholeskyBase< Derived >::m_shiftOffset [protected] |
Definition at line 262 of file SimplicialCholesky.h.
RealScalar Eigen::SimplicialCholeskyBase< Derived >::m_shiftScale [protected] |
Definition at line 263 of file SimplicialCholesky.h.