Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
solver_gmres_cg.cc
Go to the documentation of this file.
1
21// 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
35namespace rheolef {
36using namespace std;
37
38template<class T, class M>
39vec<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}
63template<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}
72template <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
determinant_type det() const
vec< T, M > trans_solve(const vec< T, M > &rhs) const
vec< T, M > solve(const vec< T, M > &rhs) const
see the vec page for the full documentation
Definition vec.h:79
#define error_macro(message)
Definition dis_macros.h:49
#define warning_macro(message)
Definition dis_macros.h:53
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
This file is part of Rheolef.
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
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
int cg(const Matrix &A, Vector &x, const Vector2 &Mb, const Preconditioner &M, const solver_option &sopt=solver_option())
Definition cg.h:94
STL namespace.
field residue(Float p, const field &uh)