GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-06-04 18:54:00 +08:00
Code Issues Packages Projects Releases Wiki Activity

extend SimplicialCholesky for sparse rhs, and add determinant

Browse Source
This commit is contained in:
Gael Guennebaud 2011-10-11 11:31:12 +02:00
parent 5dc8458293
commit 4f237f035c
1 changed files with 57 additions and 18 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

75
unsupported/Eigen/src/SparseExtra/SimplicialCholesky.h
View File

@ -150,15 +150,15 @@ class SimplicialCholeskyBase
* *
* \sa compute() * \sa compute()
*/ */
// template<typename Rhs> template<typename Rhs>
// inline const internal::sparse_solve_retval<CholmodDecomposition, Rhs> inline const internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs>
// solve(const SparseMatrixBase<Rhs>& b) const solve(const SparseMatrixBase<Rhs>& b) const
// { {
// eigen_assert(m_isInitialized && "SimplicialCholesky is not initialized."); eigen_assert(m_isInitialized && "Simplicial LLt or LDLt is not initialized.");
// eigen_assert(rows()==b.rows() eigen_assert(rows()==b.rows()
// && "SimplicialCholesky::solve(): invalid number of rows of the right hand side matrix b"); && "SimplicialCholesky::solve(): invalid number of rows of the right hand side matrix b");
// return internal::sparse_solve_retval<SimplicialCholesky, Rhs>(*this, b.derived()); return internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs>(*this, b.derived());
// } }
/** \returns the permutation P /** \returns the permutation P
* \sa permutationPinv() */ * \sa permutationPinv() */
@ -214,13 +214,26 @@ class SimplicialCholeskyBase
} }
/** \internal */ /** \internal */
/* template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex>
template<typename RhsScalar, int RhsOptions, typename RhsIndex, typename DestScalar, int DestOptions, typename DestIndex> void _solve_sparse(const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const
void _solve(const SparseMatrix<RhsScalar,RhsOptions,RhsIndex> &b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const {
{ eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
// TODO eigen_assert(m_matrix.rows()==b.rows());
}
*/ // we process the sparse rhs per block of NbColsAtOnce columns temporarily stored into a dense matrix.
static const int NbColsAtOnce = 4;
int rhsCols = b.cols();
int size = b.rows();
Eigen::Matrix<DestScalar,Dynamic,Dynamic> tmp(size,rhsCols);
for(int k=0; k<rhsCols; k+=NbColsAtOnce)
{
int actualCols = std::min<int>(rhsCols-k, NbColsAtOnce);
tmp.leftCols(actualCols) = b.middleCols(k,actualCols);
tmp.leftCols(actualCols) = derived().solve(tmp.leftCols(actualCols));
dest.middleCols(k,actualCols) = tmp.leftCols(actualCols).sparseView();
}
}
#endif // EIGEN_PARSED_BY_DOXYGEN #endif // EIGEN_PARSED_BY_DOXYGEN
protected: protected:
@ -382,6 +395,11 @@ public:
Base::template factorize<false>(a); Base::template factorize<false>(a);
} }
Scalar determinant() const
{
Scalar detL = Diagonal<const CholMatrixType>(Base::m_matrix).prod();
return internal::abs2(detL);
}
}; };
/** \class SimplicialLDLt /** \class SimplicialLDLt
@ -453,6 +471,10 @@ public:
Base::template factorize<true>(a); Base::template factorize<true>(a);
} }
Scalar determinant() const
{
return Base::m_diag.prod();
}
}; };
/** \class SimplicialCholesky /** \class SimplicialCholesky
@ -477,7 +499,10 @@ public:
public: public:
SimplicialCholesky() : Base(), m_LDLt(true) {} SimplicialCholesky() : Base(), m_LDLt(true) {}
SimplicialCholesky(const MatrixType& matrix) SimplicialCholesky(const MatrixType& matrix)
: Base(matrix), m_LDLt(true) {} : Base(), m_LDLt(true)
{
Base::compute(matrix);
}
SimplicialCholesky& setMode(SimplicialCholeskyMode mode) SimplicialCholesky& setMode(SimplicialCholeskyMode mode)
{ {
@ -567,6 +592,20 @@ public:
if(Base::m_P.size()>0) if(Base::m_P.size()>0)
dest = Base::m_P * dest; dest = Base::m_P * dest;
} }
Scalar determinant() const
{
if(m_LDLt)
{
return Base::m_diag.prod();
}
else
{
Scalar detL = Diagonal<const CholMatrixType>(Base::m_matrix).prod();
return internal::abs2(detL);
}
}
protected: protected:
bool m_LDLt; bool m_LDLt;
}; };
@ -754,7 +793,7 @@ struct sparse_solve_retval<SimplicialCholeskyBase<Derived>, Rhs>
template<typename Dest> void evalTo(Dest& dst) const template<typename Dest> void evalTo(Dest& dst) const
{ {
dec().derived()._solve(rhs(),dst); dec().derived()._solve_sparse(rhs(),dst);
} }
}; };