rheolef  6.6
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  odiststream* p_derr = 0, std::string label = "pgmres")
140 {
141  Vector w;
142  SmallVector s(m+1), cs(m+1), sn(m+1);
143  Size i;
144  Size j = 1;
145  Size k;
146  Real residue;
147  Real norm_b = norm(M.solve(b));
148  Vector r = M.solve(b - A * x);
149  Real beta = norm(r);
150  if (p_derr) (*p_derr) << "[" << label << "] # norm_b=" << norm_b << std::endl;
151  if (p_derr) (*p_derr) << "[" << label << "] #iteration residue" << std::endl;
152  if (norm_b == Real(0)) {
153  norm_b = 1;
154  }
155  residue = norm(r);
156  if (residue <= tol*max(1.,norm_b)) {
157  tol = residue/norm_b;
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  residue = abs(s(i+1));
181  if (p_derr) (*p_derr) << "[" << label << "] " << j << " " << residue/norm_b << std::endl;
182  if (residue < tol*max(1.,norm_b)) {
183  Update(x, i, H, s, v);
184  tol = residue/norm_b;
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  residue = beta;
194  if (p_derr) (*p_derr) << "[" << label << "] " << j << " " << residue/norm_b << std::endl;
195  if (residue < tol*max(1.,norm_b)) {
196  tol = residue/norm_b;
197  max_iter = j;
198  delete [] v;
199  return 0;
200  }
201  }
202  tol = residue/norm_b;
203  delete [] v;
204  return 1;
205 }
206 @endcode
207 }// namespace rheolef
208 # endif // _SKIT_GMRES_H
idiststream, odiststream - distributed interface for large data streams
Definition: diststream.h:68
int pgmres(const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M, SmallMatrix &H, const SmallVector &, const Size &m, Size &max_iter, Real &tol, odiststream *p_derr=0, std::string label="pgmres")
Definition: pgmres.h:137
void GeneratePlaneRotation(Real &dx, Real &dy, Real &cs, Real &sn)
Definition: pgmres.h:112
irheostream, orheostream - large data streams
Definition: compiler.h:7
field residue(Float p, const field &uh)
T norm(const vec< T, M > &x)
Definition: vec_expr.h:326
field_nonlinear_expr< field_nonlinear_expr_uf< details::sqrt_,field_expr_terminal_field< T, M > > > sqrt(const field_basic< T, M > &x)
void ApplyPlaneRotation(Real &dx, Real &dy, Real &cs, Real &sn)
Definition: pgmres.h:128
T dot(const vec< T, M > &x, const int &y)
boost::proto::result_of::make_expr< boost::proto::tag::function, vec_domain, const vec_detail::max_, const int &, const vec< T, M > & >::type max(const int &x, const vec< T, M > &y)
boost::proto::result_of::make_expr< boost::proto::tag::function, vec_domain, const vec_detail::abs_, const vec< T, M > & >::type abs(const vec< T, M > &x)
void Update(Vector &x, Size k, SmallMatrix &h, SmallVector &s, Vector v[])
Definition: pgmres.h:99