1 # ifndef _RHEOLEF_GMRES_H
2 # define _RHEOLEF_GMRES_H
124 #include "rheolef/solver_option.h"
131 template <
class SmallMatrix,
class SmallVector,
class Vector,
class Vector2,
class Size>
132 void update (
Vector& x, Size k,
const SmallMatrix&
h,
const SmallVector& s, Vector2& v) {
135 for (
int i = k; i >= 0; i--) {
137 for (
int j = i - 1; j >= 0; j--)
138 y(j) -=
h(j,i) * y(i);
140 for (Size j = 0; j <= k; j++) {
149 }
else if (abs(dy) > abs(dx)) {
151 sn = 1.0 / sqrt( 1.0 + temp*temp );
155 cs = 1.0 / sqrt( 1.0 + temp*temp );
161 Real temp = cs * dx + sn * dy;
162 dy = -sn * dx + cs * dy;
169 template <
class Matrix,
class Vector,
class Preconditioner,
170 class SmallMatrix,
class SmallVector>
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;
230 # endif // _RHEOLEF_GMRES_H