VERB_code_2.3
GMRES2.h
1 /*
2  * GMRES2.h
3  *
4  * Created on: Aug 31, 2010
5  * Author: subbotin
6  */
7 
8 #ifndef GMRES2_H_
9 #define GMRES2_H_
10 
11 #include <math.h>
12 
13 
14 //*****************************************************************
15 // Iterative template routine -- GMRES
16 //
17 // GMRES solves the unsymmetric linear system Ax = b using the
18 // Generalized Minimum Residual method
19 //
20 // GMRES follows the algorithm described on p. 20 of the
21 // SIAM Templates book.
22 //
23 // The return value indicates convergence within max_iter (input)
24 // iterations (0), or no convergence within max_iter iterations (1).
25 //
26 // Upon successful return, output arguments have the following values:
27 //
28 // x -- approximate solution to Ax = b
29 // max_iter -- the number of iterations performed before the
30 // tolerance was reached
31 // tol -- the residual after the final iteration
32 //
33 //*****************************************************************
34 
35 template < class Matrix, class Vector >
36 void
37 //Update(Vector &x, int k, Matrix &h, Vector &s, Vector v[])
38 Update(Vector &x, int k, Matrix &h, Vector &s, std::vector<Vector> &v)
39 {
40  Vector y(s);
41 
42  // Backsolve:
43  for (int i = k; i >= 0; i--) {
44  y(i) /= h(i,i);
45  for (int j = i - 1; j >= 0; j--)
46  y(j) -= h(j,i) * y(i);
47  }
48 
49  for (int j = 0; j <= k; j++) {
50  x += v[j] * y(j); // Need to use matrix operations here, instead of elements multiplication
51  }
52 }
53 
54 template < class Real >
55 Real
56 abs(Real x)
57 {
58  return (x > 0 ? x : -x);
59 }
60 
61 double norm(const CPPL::dcovector& V) {
62  double res = 0;
63  int i;
64  for (i = 0; i < V.l; i++) {
65  res += (V(i) * V(i));
66  }
67  return sqrt(res);
68 }
69 
70 double dot(const CPPL::dcovector& V, const CPPL::dcovector& W) {
71  double res = 0;
72  int i;
73  for (i = 0; i < V.l; i++) {
74  res += (V(i) * W(i));
75  }
76  return res;
77 }
78 
79 
80 template<class Real>
81 void GeneratePlaneRotation(Real &dx, Real &dy, Real &cs, Real &sn)
82 {
83  if (dy == 0.0) {
84  cs = 1.0;
85  sn = 0.0;
86  } else if (abs(dy) > abs(dx)) {
87  Real temp = dx / dy;
88  sn = 1.0 / sqrt( 1.0 + temp*temp );
89  cs = temp * sn;
90  } else {
91  Real temp = dy / dx;
92  cs = 1.0 / sqrt( 1.0 + temp*temp );
93  sn = temp * cs;
94  }
95 }
96 
97 
98 template<class Real>
99 void ApplyPlaneRotation(Real &dx, Real &dy, Real &cs, Real &sn)
100 {
101  Real temp = cs * dx + sn * dy;
102  dy = -sn * dx + cs * dy;
103  dx = temp;
104 }
105 
106 
107 template < class Operator, class Vector, class Preconditioner,
108  class Matrix, class Real >
109 int
110 GMRES(const Operator &A, Vector &x, const Vector &b,
111  const Preconditioner &M, Matrix &H, int &m, int &max_iter,
112  Real &tol)
113 {
114  Real resid;
115  int i, j = 1, k;
116  Vector s(m+1), cs(m+1), sn(m+1);
117  // std::vector<Real> s(m+1, 0), cs(m+1, 0), sn(m+1, 0);
118  Vector w;
119 
120  Real normb = norm(M.solve(b));
121  //Real normb = M.solve(b).norm();
122 
123  Vector r = M.solve(b - A * x);
124  Real beta = norm(r);
125  //Real beta = r.norm();
126 
127  if (normb == 0.0)
128  normb = 1;
129 
130  // initial check, if we don't really need to do anything
131  if (false && (resid = norm(r) / normb) <= tol) {
132  //if ((resid = r.norm() / normb) <= tol) {
133  tol = resid;
134  max_iter = 0;
135  return 0;
136  }
137 
138  std::vector<Vector> v(m+1, Vector());
139  //Vector *v = new Vector[m+1];
140 
141  while (j <= max_iter) {
142  v[0] = r * (1.0 / beta); // ??? r / beta
143  //s = 0.0; ????
144  for (i = 0; i < s.l; i++) s(i) = 0;
145  s(0) = beta;
146 
147  for (i = 0; i < m && j <= max_iter; i++, j++) {
148  w = M.solve(A * v[i]);
149 
150  for (k = 0; k <= i; k++) {
151  H(k, i) = dot(w, v[k]);
152  //H(k, i) = w.dot(v[k]);
153  w -= H(k, i) * v[k]; // Need to use matrix operations here, instead of elements multiplication
154  }
155  H(i+1, i) = norm(w);
156  //H(i+1, i) = w.norm();
157  v[i+1] = w * (1.0 / H(i+1, i)); // ??? w / H(i+1, i)
158 
159  for (k = 0; k < i; k++)
160  ApplyPlaneRotation(H(k,i), H(k+1,i), cs(k), sn(k));
161 
162  GeneratePlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
163  ApplyPlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
164  ApplyPlaneRotation(s(i), s(i+1), cs(i), sn(i)); //
165 
166 
167  if ((resid = abs(s(i+1)) / normb) < tol) {
168  Update(x, i, H, s, v);
169  tol = resid;
170  max_iter = j;
172  return 0;
173  }
174  }
175  Update(x, i - 1, H, s, v);
176  r = M.solve(b - A * x);
177  beta = norm(r);
178  //beta = r.norm();
179  if ((resid = beta / normb) < tol) {
180  tol = resid;
181  max_iter = j;
183  return 0;
184  }
185  }
186 
187  tol = resid;
189  return 1;
190 }
191 
192 
193 
194 #endif /* GMRES2_H_ */
static const double m
mass of electron in grams