94int 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;
123 Real alpha = norm2_r/
dot(Mq,
p);
int cg(const Matrix &A, Vector &x, const Vector2 &Mb, const Preconditioner &M, const solver_option &sopt=solver_option())