Loading...
Searching...
No Matches
solver_cholmod.cc
Go to the documentation of this file.
// direct solver cholmod, seq implementations
// use the eigen interface, for convenience
//
// - cholmod LLt supernodal : when matrix is symmetric definite positive
// TODO:
// - cholmod simplicial could be considered when matrix is sym. indefinite
//
#include "solver_cholmod.h"
#if defined(_RHEOLEF_HAVE_CHOLMOD)
namespace rheolef {
using namespace std;
template<class T, class M>
void
solver_cholmod_rep<T,M>::update_values (const csr<T,M>& a)
{
// copy into eigen; then eigen->cholmod without copy: use pointers to arrays (ia,ja,a)
_a = a;
if (_a.nnz() == 0) {
return; // empty matrix
}
using namespace Eigen;
Matrix<int,Dynamic,1> nnz_row (a.nrow());
typename csr<T,M>::const_iterator ia = a.begin();
for (size_type i = 0, q = 0, n = a.nrow(); i < n; ++i) {
nnz_row[i] = ia[i+1] - ia[i];
}
SparseMatrix<T> a_tmp (a.nrow(),a.ncol());
a_tmp.reserve (nnz_row);
for (size_type i = 0, n = a.nrow(); i < n; ++i) {
for (typename csr<T,M>::const_data_iterator p = ia[i]; p < ia[i+1]; ++p) {
a_tmp.insert (i, (*p).first) = (*p).second;
}
}
a_tmp.makeCompressed();
_llt_a.compute (a_tmp);
check_macro (_llt_a.info() == Success, "cholmod LLt factorization failed: non-positive definite matrix");
// TODO: note: when matrix is distributed and block diagonal, each proc returns its block determinant
// => do the product
if (base::option().compute_determinant) {
T det_a = _llt_a.determinant();
_det.mantissa = 1; // sign, always positive
_det.exponant = _llt_a.logDeterminant() / log(T(10));
_det.base = 10;
}
}
template<class T, class M>
vec<T,M>
solver_cholmod_rep<T,M>::solve (const vec<T,M>& b) const
{
vec<T,M> x(b.ownership());
if (_a.nnz() == 0) return x; // empty matrix
using namespace Eigen;
Map<Matrix<T,Dynamic,1> > b_map ((T*)(&(*b.begin())), b.size()),
x_map ( &(*x.begin()), x.size());
x_map = _llt_a.solve (b_map);
return x;
}
template<class T, class M>
vec<T,M>
solver_cholmod_rep<T,M>::trans_solve (const vec<T,M>& b) const
{
return solve(b);
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
template class solver_cholmod_rep<double,sequential>;
#ifdef _RHEOLEF_HAVE_MPI
template class solver_cholmod_rep<double,distributed>;
#endif // _RHEOLEF_HAVE_MPI
} // namespace rheolef
#endif // _RHEOLEF_HAVE_CHOLMOD
rheolef::csr
see the csr page for the full documentation
Definition csr.h:317
rheolef::solver_abstract_rep::option
const solver_option & option() const
Definition solver.h:197
rheolef::solver_cholmod_rep
Definition solver_cholmod.h:40
rheolef::solver_cholmod_rep::_a
csr< T, M > _a
Definition solver_cholmod.h:64
rheolef::solver_cholmod_rep::_llt_a
Eigen::CholmodDecomposition< Eigen::SparseMatrix< double > > _llt_a
Definition solver_cholmod.h:66
rheolef::solver_cholmod_rep::update_values
void update_values(const csr< T, M > &a)
Definition solver_cholmod.cc:35
rheolef::solver_cholmod_rep::size_type
base::size_type size_type
Definition solver_cholmod.h:45
rheolef::solver_cholmod_rep::trans_solve
vec< T, M > trans_solve(const vec< T, M > &rhs) const
Definition solver_cholmod.cc:81
rheolef::solver_cholmod_rep::solve
vec< T, M > solve(const vec< T, M > &rhs) const
Definition solver_cholmod.cc:69
rheolef::solver_cholmod_rep::_det
determinant_type _det
Definition solver_cholmod.h:67
rheolef::vec
see the vec page for the full documentation
Definition vec.h:79
T
Expr1::float_type T
Definition field_expr.h:230
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
Eigen
Definition field_eigen.h:35
rheolef
This file is part of Rheolef.
Definition compiler_eigen.h:39
rheolef::solve
void solve(tiny_matrix< T > &a, tiny_vector< size_t > &piv, const tiny_vector< T > &b, tiny_vector< T > &x)
Definition tiny_lu.h:92
solver_cholmod.h
p
Definition sphere.icc:25