88 #include "rheolef/solver_option.h"
93 template <
class Matrix,
class Vector,
class Vector2,
class Preconditioner>
94 int cg (
const Matrix &A,
Vector &x,
const Vector2 &Mb,
const Preconditioner &
M,
100 typedef typename Vector::float_type Real;
101 std::string label = (sopt.label !=
"" ? sopt.label :
"cg");
103 Real norm2_b =
dot(Mb,
b);
104 if (sopt.absolute_stopping || norm2_b == Real(0)) norm2_b = 1;
107 if (sopt.p_err) (*sopt.p_err) <<
"[" << label <<
"] #iteration residue" << std::endl;
109 for (sopt.n_iter = 0; sopt.n_iter <= sopt.max_iter; sopt.n_iter++) {
111 Real prev_norm2_r = norm2_r;
112 norm2_r =
dot(Mr, r);
113 sopt.residue = sqrt(norm2_r/norm2_b);
114 if (sopt.p_err) (*sopt.p_err) <<
"[" << label <<
"] " << sopt.n_iter <<
" " << sopt.residue << std::endl;
115 if (sopt.residue <= sopt.tol)
return 0;
116 if (sopt.n_iter == 0) {
119 Real
beta = norm2_r/prev_norm2_r;
see the solver_option page for the full documentation
Float alpha[pmax+1][pmax+1]
This file is part of Rheolef.
rheolef::std enable_if ::type dot const Expr1 expr1, const Expr2 expr2 dot(const Expr1 &expr1, const Expr2 &expr2)
dot(x,y): see the expression page for the full documentation
int cg(const Matrix &A, Vector &x, const Vector2 &Mb, const Preconditioner &M, const solver_option &sopt=solver_option())
DATA::size_type size_type