23 template <
class Matrix,
class Vector >
25 Update(Vector &x,
int k, Matrix &h, Vector &s, Vector v[])
30 for (
int i = k; i >= 0; i--) {
32 for (
int j = i - 1; j >= 0; j--)
33 y(j) -= h(j,i) * y(i);
36 for (
int j = 0; j <= k; j++)
41 template <
class Real >
45 return (x > 0 ? x : -x);
49 template <
class Operator,
class Vector,
class Preconditioner,
50 class Matrix,
class Real >
52 GMRES(
const Operator &A, Vector &x,
const Vector &b,
53 const Preconditioner &M, Matrix &H,
int &m,
int &max_iter,
58 Vector s(m+1), cs(m+1), sn(m+1), w;
60 Real normb =
norm(M.solve(b));
61 Vector r = M.solve(b - A * x);
67 if ((resid =
norm(r) / normb) <= tol) {
73 Vector *v =
new Vector[m+1];
75 while (j <= max_iter) {
76 v[0] = r * (1.0 / beta);
80 for (i = 0; i < m && j <= max_iter; i++, j++) {
81 w = M.solve(A * v[i]);
82 for (k = 0; k <= i; k++) {
83 H(k, i) =
dot(w, v[k]);
87 v[i+1] = w * (1.0 / H(i+1, i));
89 for (k = 0; k < i; k++)
96 if ((resid =
abs(s(i+1)) / normb) < tol) {
104 Update(x, m - 1, H, s, v);
105 r = M.solve(b - A * x);
107 if ((resid = beta / normb) < tol) {
130 }
else if (
abs(dy) >
abs(dx)) {
132 sn = 1.0 / sqrt( 1.0 + temp*temp );
136 cs = 1.0 / sqrt( 1.0 + temp*temp );
145 Real temp = cs * dx + sn * dy;
146 dy = -sn * dx + cs * dy;
void Update(Vector &x, int k, Matrix &h, Vector &s, Vector v[])
void ApplyPlaneRotation(Real &dx, Real &dy, Real &cs, Real &sn)
void GeneratePlaneRotation(Real &dx, Real &dy, Real &cs, Real &sn)
int GMRES(const Operator &A, Vector &x, const Vector &b, const Preconditioner &M, Matrix &H, int &m, int &max_iter, Real &tol)
A::value_type dot(const A &a, const B &b)
Dot product of array or expression.
A::value_type norm(const A &a)
Euklidean norm of array or expression.