solver_gmres_cg.cc

Go to the documentation of this file.

// solver wrapper: direct or iterative
//
#include "solver_gmres_cg.h"
 
#include "rheolef/cg.h"
#include "rheolef/gmres.h"
#include "rheolef/vec_expr_v2.h"
#include "rheolef/eye.h"
 
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Weffc++"
#include <Eigen/Dense>
#pragma GCC diagnostic pop
 
namespace rheolef {
using namespace std;
 
template<class T, class M>
vec<T,M>
solver_gmres_cg_rep<T,M>::solve (const vec<T,M>& b) const
{
  if (!_precond.initialized()) {
    _precond = eye_basic<T,M>();
  }
  vec<T,M> x (b.ownership(), 0);
  if (_a.is_symmetric() && _a.is_definite_positive()) {
    int status = cg (_a, x, b, _precond, base::option());
  } else {
    using namespace Eigen;
    size_type m = base::option().krylov_dimension;
    Matrix<T,Dynamic,Dynamic> h(m+1,m+1);
    Matrix<T,Dynamic,1>   dummy(m);
    int status = gmres (_a, x, b, _precond, h, dummy, base::option());
  }
  check_macro (base::option().residue <= sqrt(base::option().tol),
        "solver: precision "<<base::option().tol<<" not reached: get "<<base::option().residue
        << " after " << base::option().max_iter << " iterations");
  if (base::option().residue > base::option().tol) warning_macro (
        "solver: precision "<<base::option().tol<<" not reached: get "<<base::option().residue
        << " after " << base::option().max_iter << " iterations");
  return x;
}
template<class T, class M>
vec<T,M>
solver_gmres_cg_rep<T,M>::trans_solve (const vec<T,M>& b) const
{
  if (_a.is_symmetric()) return solve(b);
  // TODO: wrap a as a*x return a.trans_mult(x) and the same with _precond.trans_solve(b)
  error_macro ("iterative trans_solve: not yet");
  return b;
}
template <class T, class M>
typename solver_gmres_cg_rep<T,M>::determinant_type
solver_gmres_cg_rep<T,M>::det() const
{
  error_macro ("undefined determinant computation for iterative solver (HINT: use a direct method)");
  return determinant_type();
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
 
template class solver_gmres_cg_rep<Float,sequential>;
#ifdef _RHEOLEF_HAVE_MPI
template class solver_gmres_cg_rep<Float,distributed>;
#endif // _RHEOLEF_HAVE_MPI
 
} // namespace rheolef