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