Loading...
Searching...
No Matches
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
rheolef::solver_abstract_rep::option
const solver_option & option() const
Definition solver.h:197
rheolef::solver_gmres_cg_rep
Definition solver_gmres_cg.h:27
rheolef::solver_gmres_cg_rep::_precond
solver_basic< T, M > _precond
Definition solver_gmres_cg.h:52
rheolef::solver_gmres_cg_rep::_a
csr< T, M > _a
Definition solver_gmres_cg.h:51
rheolef::solver_gmres_cg_rep::det
determinant_type det() const
Definition solver_gmres_cg.cc:74
rheolef::solver_gmres_cg_rep::size_type
base::size_type size_type
Definition solver_gmres_cg.h:32
rheolef::solver_gmres_cg_rep::determinant_type
base::determinant_type determinant_type
Definition solver_gmres_cg.h:33
rheolef::solver_gmres_cg_rep::trans_solve
vec< T, M > trans_solve(const vec< T, M > &rhs) const
Definition solver_gmres_cg.cc:65
rheolef::solver_gmres_cg_rep::solve
vec< T, M > solve(const vec< T, M > &rhs) const
Definition solver_gmres_cg.cc:40
rheolef::vec
see the vec page for the full documentation
Definition vec.h:79
error_macro
#define error_macro(message)
Definition dis_macros.h:49
warning_macro
#define warning_macro(message)
Definition dis_macros.h:53
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
rheolef::gmres
int gmres(const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M, SmallMatrix &H, const SmallVector &V, const solver_option &sopt=solver_option())
Definition gmres.h:171
rheolef::eye_basic
solver_basic< T, M > eye_basic()
Definition eye.h:93
rheolef::cg
int cg(const Matrix &A, Vector &x, const Vector2 &Mb, const Preconditioner &M, const solver_option &sopt=solver_option())
Definition cg.h:94
residue
field residue(Float p, const field &uh)
Definition p_laplacian_post.cc:35
solver_gmres_cg.h