Rheolef  7.1
an efficient C++ finite element environment
solver_gmres_cg.cc
Go to the documentation of this file.
1 // solver wrapper: direct or iterative
22 //
23 #include "solver_gmres_cg.h"
24 
25 #include "rheolef/cg.h"
26 #include "rheolef/gmres.h"
27 #include "rheolef/vec_expr_v2.h"
28 #include "rheolef/eye.h"
29 
30 #pragma GCC diagnostic push
31 #pragma GCC diagnostic ignored "-Weffc++"
32 #include <Eigen/Dense>
33 #pragma GCC diagnostic pop
34 
35 namespace rheolef {
36 using namespace std;
37 
38 template<class T, class M>
39 vec<T,M>
41 {
42  if (!_precond.initialized()) {
43  _precond = eye_basic<T,M>();
44  }
45  vec<T,M> x (b.ownership(), 0);
46  if (_a.is_symmetric() && _a.is_definite_positive()) {
47  int status = cg (_a, x, b, _precond, base::option());
48  } else {
49  using namespace Eigen;
50  size_type m = base::option().krylov_dimension;
51  Matrix<T,Dynamic,Dynamic> h(m+1,m+1);
52  Matrix<T,Dynamic,1> dummy(m);
53  int status = gmres (_a, x, b, _precond, h, dummy, base::option());
54  }
55  check_macro (base::option().residue <= sqrt(base::option().tol),
56  "solver: precision "<<base::option().tol<<" not reached: get "<<base::option().residue
57  << " after " << base::option().max_iter << " iterations");
58  if (base::option().residue > base::option().tol) warning_macro (
59  "solver: precision "<<base::option().tol<<" not reached: get "<<base::option().residue
60  << " after " << base::option().max_iter << " iterations");
61  return x;
62 }
63 template<class T, class M>
66 {
67  if (_a.is_symmetric()) return solve(b);
68  // TODO: wrap a as a*x return a.trans_mult(x) and the same with _precond.trans_solve(b)
69  error_macro ("iterative trans_solve: not yet");
70  return b;
71 }
72 template <class T, class M>
75 {
76  error_macro ("undefined determinant computation for iterative solver (HINT: use a direct method)");
77  return determinant_type();
78 }
79 // ----------------------------------------------------------------------------
80 // instanciation in library
81 // ----------------------------------------------------------------------------
82 
84 #ifdef _RHEOLEF_HAVE_MPI
86 #endif // _RHEOLEF_HAVE_MPI
87 
88 } // namespace rheolef
Eigen
Definition: field_eigen.h:35
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
warning_macro
#define warning_macro(message)
Definition: dis_macros.h:53
mkgeo_contraction.h
h
Definition: mkgeo_contraction.sh:179
residue
field residue(Float p, const field &uh)
Definition: p_laplacian_post.cc:35
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::dummy
static iorheo::force_initialization dummy
Definition: iorheo.cc:147
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
rheolef::solver_abstract_rep::determinant_type
Definition: solver.h:215
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::solver_gmres_cg_rep::det
determinant_type det() const
Definition: solver_gmres_cg.cc:74
mkgeo_sector.m
m
Definition: mkgeo_sector.sh:118
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
error_macro
#define error_macro(message)
Definition: dis_macros.h:49
rheolef::solver_abstract_rep::size_type
csr< T, M >::size_type size_type
Definition: solver.h:193
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
mkgeo_ball.b
b
Definition: mkgeo_ball.sh:152
rheolef::solver_gmres_cg_rep
Definition: solver_gmres_cg.h:27
solver_gmres_cg.h
mkgeo_contraction.status
status
Definition: mkgeo_contraction.sh:290