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) {
122 Real beta = sqrt(beta2);
130 Vector u_old (x.ownership(), 0.);
131 Vector Mv_old (x.ownership(), 0.);
132 Vector w (x.ownership(), 0.);
133 Vector w_old (x.ownership(), 0.);
134 Vector w_old2 (x.ownership(), 0.);
135 for (sopt.n_iter = 1; sopt.n_iter <= sopt.max_iter; sopt.n_iter++) {
139 Real alpha =
dot(Mr,
u);
140 Mr = Mr - alpha*Mv - beta*Mv_old;
141 z = z - alpha*
u - beta*u_old;
145 sopt.residue = norm_r/norm_b;
148 Real beta_old = beta;
155 Real rho0 = c_old*alpha - c_old2*s_old*beta_old;
156 Real rho2 = s_old*alpha + c_old2*c_old*beta_old;
157 Real rho1 = sqrt(sqr(rho0) + sqr(beta));
158 Real rho3 = s_old2 * beta_old;
165 w = (
u - rho2*w_old - rho3*w_old2)/rho1;
174 sopt.residue = norm_r/norm_b;
175 if (sopt.p_err) (*sopt.p_err) <<
"[" << label <<
"] " << sopt.n_iter <<
" " << sopt.residue << std::endl;
176 if (sopt.residue <= sopt.tol)
return 0;
int minres(const Matrix &A, Vector &x, const Vector &Mb, const Preconditioner &M, const solver_option &sopt=solver_option())