rheolef  6.5
pgmres.h
Go to the documentation of this file.
1 # ifndef _SKIT_GMRES_H
2 # define _SKIT_GMRES_H
3 
4 #ifdef HAVE_CONFIG_H
5 #include "rheolef/compiler.h"
6 #endif
7 #include <cmath>
8 
95 namespace rheolef {
96 
97 template <class SmallMatrix, class Vector, class SmallVector, class Size>
98 void
99 Update(Vector &x, Size k, SmallMatrix &h, SmallVector &s, Vector v[])
100 {
101  SmallVector y = s;
102  for (int i = k; i >= 0; i--) {
103  y(i) /= h(i,i);
104  for (int j = i - 1; j >= 0; j--)
105  y(j) -= h(j,i) * y(i);
106  }
107  for (Size j = 0; j <= k; j++) {
108  x += v[j] * y(j);
109  }
110 }
111 template<class Real>
112 void GeneratePlaneRotation(Real &dx, Real &dy, Real &cs, Real &sn)
113 {
114  if (dy == Real(0)) {
115  cs = 1.0;
116  sn = 0.0;
117  } else if (abs(dy) > abs(dx)) {
118  Real temp = dx / dy;
119  sn = 1.0 / ::sqrt( 1.0 + temp*temp );
120  cs = temp * sn;
121  } else {
122  Real temp = dy / dx;
123  cs = 1.0 / ::sqrt( 1.0 + temp*temp );
124  sn = temp * cs;
125  }
126 }
127 template<class Real>
128 void ApplyPlaneRotation(Real &dx, Real &dy, Real &cs, Real &sn)
129 {
130  Real temp = cs * dx + sn * dy;
131  dy = -sn * dx + cs * dy;
132  dx = temp;
133 }
134 template <class Matrix, class Vector, class Preconditioner,
135  class SmallMatrix, class SmallVector, class Real, class Size>
136 int
137 pgmres (const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M,
138  SmallMatrix &H, const SmallVector&, const Size &m, Size &max_iter, Real &tol)
139 {
140  odiststream* p_derr = &derr;
141  std::string label = "pgmres";
142  if (p_derr) (*p_derr) << "[" << label << "] #iteration residue" << std::endl;
143  Vector w;
144  SmallVector s(m+1), cs(m+1), sn(m+1);
145  Size i;
146  Size j = 1;
147  Size k;
148  Real resid;
149  Real normb = norm(M.solve(b));
150  Vector r = M.solve(b - A * x);
151  Real beta = norm(r);
152  if (normb == Real(0)) {
153  normb = 1;
154  }
155  resid = norm(r) / normb;
156  if (resid <= tol) {
157  tol = resid;
158  max_iter = 0;
159  return 0;
160  }
161  Vector *v = new Vector[m+1];
162  while (j <= max_iter) {
163  v[0] = r * (1.0 / beta); // ??? r / beta
164  s = 0.0;
165  s(0) = beta;
166  for (i = 0; i < m && j <= max_iter; i++, j++) {
167  w = M.solve(A * v[i]);
168  for (k = 0; k <= i; k++) {
169  H(k, i) = dot(w, v[k]);
170  w -= H(k, i) * v[k];
171  }
172  H(i+1, i) = norm(w);
173  v[i+1] = w * (1.0 / H(i+1, i)); // ??? w / H(i+1, i)
174  for (k = 0; k < i; k++) {
175  ApplyPlaneRotation(H(k,i), H(k+1,i), cs(k), sn(k));
176  }
177  GeneratePlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
178  ApplyPlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
179  ApplyPlaneRotation(s(i), s(i+1), cs(i), sn(i));
180  resid = abs(s(i+1)) / normb;
181  if (p_derr) (*p_derr) << "[" << label << "] " << j << " " << resid << std::endl;
182  if (resid < tol) {
183  Update(x, i, H, s, v);
184  tol = resid;
185  max_iter = j;
186  delete [] v;
187  return 0;
188  }
189  }
190  Update(x, m - 1, H, s, v);
191  r = M.solve(b - A * x);
192  beta = norm(r);
193  resid = beta / normb;
194  if (p_derr) (*p_derr) << "[" << label << "] " << j << " " << resid << std::endl;
195  if (resid < tol) {
196  tol = resid;
197  max_iter = j;
198  delete [] v;
199  return 0;
200  }
201  }
202  tol = resid;
203  delete [] v;
204  return 1;
205 }
206 @endcode
207 }// namespace rheolef
208 # endif // _SKIT_GMRES_H
209