solver_eigen.cc

Go to the documentation of this file.

// direct solver eigen, seq implementations
//
// Note : why a dis implementation based on eigen ?
// Because when dis_ext_nnz == 0, then the matrix is block diagonal.
// in that case the eigen could be better than mumps that initialize stuff
// for the distributed case.
// Is could appends e.g. for block-diagonal mass matrix "Pkd"
// This also occurs when nproc==1.
//
#include "rheolef/config.h"
#ifdef _RHEOLEF_HAVE_EIGEN
#include "solver_eigen.h"
 
namespace rheolef {
using namespace std;
 
// =========================================================================
// the class interface
// =========================================================================
template<class T, class M>
void
solver_eigen_rep<T,M>::update_values (const csr<T,M>& a)
{
  _a = a;
  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());
  if (a.nrow() != 0) {
    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();
  if (_a.is_symmetric() && _a.is_definite_positive()) {
    _ldlt_a.compute (a_tmp);
    check_macro (_ldlt_a.info() == Success, "eigen ldlt factorization failed");
    if (base::option().compute_determinant) {
      T det_a = _ldlt_a.determinant(); // TODO: wait for eigen::ldlt to support logAbsDeterminant & signDeterminant
      _det.mantissa = (det_a >= 0) ? 1 : -1;
      _det.exponant = log(fabs(det_a)) / log(T(10));
      _det.base     = 10;
    }
  } else {
    if (a.nrow() != 0) {
      _superlu_a.analyzePattern (a_tmp); 
      _superlu_a.factorize      (a_tmp);
      check_macro (_superlu_a.info() == Success, "eigen lu factorization failed");
    }
    if (base::option().compute_determinant) {
      _det.mantissa = _superlu_a.signDeterminant();
      _det.exponant = _superlu_a.logAbsDeterminant() / log(T(10));
      _det.base     = 10;
    }
  }
}
template<class T, class M>
vec<T,M>
solver_eigen_rep<T,M>::solve (const vec<T,M>& b) const
{
  if (b.dis_size() == 0) return b; // empty matrix
  vec<T,M> x (b.ownership());
  using namespace Eigen;
  Map<Matrix<T,Dynamic,1> > b_map ((T*)(b.begin().operator->()), b.size()),
                            x_map (     x.begin().operator->(),  x.size());
  if (_a.is_symmetric() && _a.is_definite_positive()) {
    x_map =    _ldlt_a.solve (b_map);
  } else {
    x_map = _superlu_a.solve (b_map);
  }
  return x;
}
template<class T, class M>
vec<T,M>
solver_eigen_rep<T,M>::trans_solve (const vec<T,M>& b) const
{
  if (_a.is_symmetric()) return solve(b);
  if (b.dis_size() == 0) return b; // empty matrix
  vec<T,M> x(b.ownership());
  fatal_macro ("eigen superlu trans_solve: not yet implemented, sorry");
#ifdef TODO
  x_map = _superlu_a.colssPermutation()*b_map; 
  _superlu_a.matrixU().transSolveInPlace (x_map);
  _superlu_a.matrixL().transSolveInPlace (x_map); // L & U trans_solve not implemented
  x_map = _superlu_a.rowsPermutation()*x_map; 
#endif // TODO
  return x;
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
 
template class solver_eigen_rep<Float,sequential>;
 
#ifdef _RHEOLEF_HAVE_MPI
template class solver_eigen_rep<Float,distributed>;
#endif // _RHEOLEF_HAVE_MPI
 
} // namespace rheolef
#endif // HAVE_EIGEN