Rheolef  7.1
an efficient C++ finite element environment
gmres.h
Go to the documentation of this file.
1 # ifndef _RHEOLEF_GMRES_H
2 # define _RHEOLEF_GMRES_H
3 //
4 // This file is part of Rheolef.
5 //
6 // Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
7 //
8 // Rheolef is free software; you can redistribute it and/or modify
9 // it under the terms of the GNU General Public License as published by
10 // the Free Software Foundation; either version 2 of the License, or
11 // (at your option) any later version.
12 //
13 // Rheolef is distributed in the hope that it will be useful,
14 // but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 // GNU General Public License for more details.
17 //
18 // You should have received a copy of the GNU General Public License
19 // along with Rheolef; if not, write to the Free Software
20 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
21 //
22 // =========================================================================
23 // AUTHOR: Pierre Saramito
24 // DATE: 12 march 1997
25 
26 
27 namespace rheolef {
120 } // namespace rheolef
121 
122 // ========== [ method implementation ] ====================================
123 
124 #include "rheolef/solver_option.h"
125 #include <cmath>
126 
127 namespace rheolef {
128 
129 namespace details {
130 // [verbatim_gmres_cont]
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) {
133  SmallVector y = s;
134  // back solve:
135  for (int i = k; i >= 0; i--) {
136  y(i) /= h(i,i);
137  for (int j = i - 1; j >= 0; j--)
138  y(j) -= h(j,i) * y(i);
139  }
140  for (Size j = 0; j <= k; j++) {
141  x += v[j] * y(j);
142  }
143 }
144 template<class Real>
145 void generate_plane_rotation (const Real& dx, const Real& dy, Real& cs, Real& sn) {
146  if (dy == Real(0)) {
147  cs = 1.0;
148  sn = 0.0;
149  } else if (abs(dy) > abs(dx)) {
150  Real temp = dx / dy;
151  sn = 1.0 / sqrt( 1.0 + temp*temp );
152  cs = temp * sn;
153  } else {
154  Real temp = dy / dx;
155  cs = 1.0 / sqrt( 1.0 + temp*temp );
156  sn = temp * cs;
157  }
158 }
159 template<class Real>
160 void apply_plane_rotation (Real& dx, Real& dy, const Real& cs, const Real& sn) {
161  Real temp = cs * dx + sn * dy;
162  dy = -sn * dx + cs * dy;
163  dx = temp;
164 }
165 // [verbatim_gmres_cont]
166 } // namespace details
167 
168 // [verbatim_gmres_synopsis]
169 template <class Matrix, class Vector, class Preconditioner,
170  class SmallMatrix, class SmallVector>
171 int gmres (const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M,
172  SmallMatrix &H, const SmallVector& V, const solver_option& sopt = solver_option())
173 // [verbatim_gmres_synopsis]
174 // [verbatim_gmres]
175 {
176  typedef typename Vector::size_type Size;
177  typedef typename Vector::float_type Real;
178  std::string label = (sopt.label != "" ? sopt.label : "gmres");
179  Size m = sopt.krylov_dimension;
180  Vector w;
181  SmallVector s(m+1), cs(m+1), sn(m+1);
182  Real residue;
183  Real norm_b = norm(M.solve(b));
184  Vector r = M.solve(b - A * x);
185  Real beta = norm(r);
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;
189  sopt.n_iter = 0;
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; ) {
194  v[0] = r/beta;
195  for (Size i = 0; i < m+1; i++) s(i) = 0; // std::numeric_limits<Float>::max();
196  s(0) = beta;
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]);
201  w -= H(k, i) * v[k];
202  }
203  H(i+1, i) = norm(w);
204  v[i+1] = w/H(i+1,i);
205  for (Size k = 0; k < i; k++) {
206  details::apply_plane_rotation (H(k,i), H(k+1,i), cs(k), sn(k));
207  }
208  details::generate_plane_rotation (H(i,i), H(i+1,i), cs(i), sn(i));
209  details::apply_plane_rotation (H(i,i), H(i+1,i), cs(i), sn(i));
210  details::apply_plane_rotation (s(i), s(i+1), cs(i), sn(i));
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) {
214  details::update (x, i, H, s, v);
215  return 0;
216  }
217  }
218  details::update (x, m - 1, H, s, v);
219  r = M.solve(b - A * x);
220  beta = norm(r);
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;
224  }
225  return 1;
226 }
227 // [verbatim_gmres]
228 
229 }// namespace rheolef
230 # endif // _RHEOLEF_GMRES_H
rheolef::Vector
Definition: Vector.h:175
rheolef::dot
rheolef::std enable_if ::type dot const Expr1 expr1, const Expr2 expr2 dot(const Expr1 &expr1, const Expr2 &expr2)
dot(x,y): see the expression page for the full documentation
Definition: vec_expr_v2.h:415
mkgeo_contraction.h
h
Definition: mkgeo_contraction.sh:179
residue
field residue(Float p, const field &uh)
Definition: p_laplacian_post.cc:35
rheolef::norm
T norm(const vec< T, M > &x)
norm(x): see the expression page for the full documentation
Definition: vec.h:387
rheolef::details::generate_plane_rotation
void generate_plane_rotation(const Real &dx, const Real &dy, Real &cs, Real &sn)
Definition: gmres.h:145
rheolef::gmres
int gmres(const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M, SmallMatrix &H, const SmallVector &V, const solver_option &sopt=solver_option())
Definition: gmres.h:171
rheolef::details::update
void update(Vector &x, Size k, const SmallMatrix &h, const SmallVector &s, Vector2 &v)
Definition: gmres.h:132
rheolef::Vector::size_type
DATA::size_type size_type
Definition: Vector.h:188
mkgeo_sector.m
m
Definition: mkgeo_sector.sh:118
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
mkgeo_ball.b
b
Definition: mkgeo_ball.sh:152
rheolef::details::apply_plane_rotation
void apply_plane_rotation(Real &dx, Real &dy, const Real &cs, const Real &sn)
Definition: gmres.h:160
M
Expr1::memory_type M
Definition: vec_expr_v2.h:385
rk::beta
Float beta[][pmax+1]
Definition: runge_kutta_semiimplicit.icc:60
rheolef::solver_option
see the solver_option page for the full documentation
Definition: solver_option.h:155