1 # ifndef _RHEOLEF_MINRES_H
2 # define _RHEOLEF_MINRES_H
94 #include "rheolef/solver_option.h"
99 template <
class Matrix,
class Vector,
class Preconditioner>
107 typedef typename Vector::float_type Real;
108 std::string label = (sopt.label !=
"" ? sopt.label :
"minres");
110 Real norm_b = sqrt(fabs(
dot(Mb,
b)));
111 if (sopt.absolute_stopping || norm_b == Real(0.)) norm_b = 1;
114 Real beta2 =
dot(Mr, z);
115 Real norm_r = sqrt(fabs(beta2));
116 sopt.residue = norm_r/norm_b;
117 if (sopt.p_err) (*sopt.p_err) <<
"[" << label <<
"] #iteration residue" << std::endl
118 <<
"[" << label <<
"] 0 " << sopt.residue << std::endl;
119 if (beta2 < 0 || sopt.residue <= sopt.tol) {
123 Real
beta = sqrt(beta2);
131 Vector u_old (x.ownership(), 0.);
132 Vector Mv_old (x.ownership(), 0.);
133 Vector w (x.ownership(), 0.);
134 Vector w_old (x.ownership(), 0.);
135 Vector w_old2 (x.ownership(), 0.);
136 for (sopt.n_iter = 1; sopt.n_iter <= sopt.max_iter; sopt.n_iter++) {
146 sopt.residue = norm_r/norm_b;
149 Real beta_old =
beta;
156 Real rho0 = c_old*
alpha - c_old2*s_old*beta_old;
157 Real rho2 = s_old*
alpha + c_old2*c_old*beta_old;
158 Real rho1 = sqrt(sqr(rho0) + sqr(
beta));
159 Real rho3 = s_old2 * beta_old;
166 w = (
u - rho2*w_old - rho3*w_old2)/rho1;
175 sopt.residue = norm_r/norm_b;
176 if (sopt.p_err) (*sopt.p_err) <<
"[" << label <<
"] " << sopt.n_iter <<
" " << sopt.residue << std::endl;
177 if (sopt.residue <= sopt.tol)
return 0;
184 # endif // _RHEOLEF_MINRES_H