177 typedef typename Vector::float_type Real;
178 std::string label = (sopt.label !=
"" ? sopt.label :
"gmres");
179 Size m = sopt.krylov_dimension;
181 SmallVector s(m+1), cs(m+1), sn(m+1);
183 Real norm_b =
norm(
M.solve(b));
186 if (sopt.p_err) (*sopt.p_err) <<
"[" << label <<
"] # norm_b=" << norm_b << std::endl
187 <<
"[" << label <<
"] #iteration residue" << std::endl;
188 if (sopt.absolute_stopping || norm_b == Real(0)) norm_b = 1;
190 sopt.residue =
norm(r)/norm_b;
191 if (sopt.residue <= sopt.tol)
return 0;
192 std::vector<Vector> v (m+1);
193 for (sopt.n_iter = 1; sopt.n_iter <= sopt.max_iter; ) {
195 for (Size i = 0; i < m+1; i++) s(i) = 0;
197 for (Size i = 0; i < m && sopt.n_iter <= sopt.max_iter; i++, sopt.n_iter++) {
198 w =
M.solve(
A * v[i]);
199 for (Size k = 0; k <= i; k++) {
200 H(k, i) =
dot(w, v[k]);
205 for (Size k = 0; k < i; k++) {
211 sopt.residue = abs(s(i+1))/norm_b;
212 if (sopt.p_err) (*sopt.p_err) <<
"[" << label <<
"] " << sopt.n_iter <<
" " << sopt.residue << std::endl;
213 if (sopt.residue <= sopt.tol) {
219 r =
M.solve(b -
A * x);
221 sopt.residue = beta/norm_b;
222 if (sopt.p_err) (*sopt.p_err) <<
"[" << label <<
"] " << sopt.n_iter <<
" " << sopt.residue << std::endl;
223 if (sopt.residue < sopt.tol)
return 0;
int gmres(const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M, SmallMatrix &H, const SmallVector &V, const solver_option &sopt=solver_option())
field residue(Float p, const field &uh)