Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
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
27namespace rheolef {
120} // namespace rheolef
121
122// ========== [ method implementation ] ====================================
123
124#include "rheolef/solver_option.h"
125#include <cmath>
126
127namespace rheolef {
128
129namespace details {
130// [verbatim_gmres_cont]
131template <class SmallMatrix, class SmallVector, class Vector, class Vector2, class Size>
132void 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}
144template<class Real>
145void 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}
159template<class Real>
160void 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]
169template <class Matrix, class Vector, class Preconditioner,
170 class SmallMatrix, class SmallVector>
171int 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
see the solver_option page for the full documentation
void update(Vector &x, Size k, const SmallMatrix &h, const SmallVector &s, Vector2 &v)
Definition gmres.h:132
void apply_plane_rotation(Real &dx, Real &dy, const Real &cs, const Real &sn)
Definition gmres.h:160
void generate_plane_rotation(const Real &dx, const Real &dy, Real &cs, Real &sn)
Definition gmres.h:145
This file is part of Rheolef.
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::std enable_if ::type dot const Expr1 expr1, const Expr2 expr2 dot(const Expr1 &expr1, const Expr2 &expr2)
T norm(const vec< T, M > &x)
norm(x): see the expression page for the full documentation
Definition vec.h:387
field residue(Float p, const field &uh)
DATA::size_type size_type
Definition Vector.h:188
Expr1::memory_type M