Loading...
Searching...
No Matches
gmres.h
Go to the documentation of this file.
# ifndef _RHEOLEF_GMRES_H
# define _RHEOLEF_GMRES_H
//
// This file is part of Rheolef.
//
// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
//
// Rheolef is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// Rheolef is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with Rheolef; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// =========================================================================
// AUTHOR: Pierre Saramito
// DATE: 12 march 1997
namespace rheolef {
} // namespace rheolef
// ========== [ method implementation ] ====================================
#include "rheolef/solver_option.h"
#include <cmath>
namespace rheolef {
namespace details {
// [verbatim_gmres_cont]
template <class SmallMatrix, class SmallVector, class Vector, class Vector2, class Size>
void update (Vector& x, Size k, const SmallMatrix& h, const SmallVector& s, Vector2& v) {
SmallVector y = s;
// back solve:
for (int i = k; i >= 0; i--) {
y(i) /= h(i,i);
for (int j = i - 1; j >= 0; j--)
y(j) -= h(j,i) * y(i);
}
for (Size j = 0; j <= k; j++) {
x += v[j] * y(j);
}
}
template<class Real>
void generate_plane_rotation (const Real& dx, const Real& dy, Real& cs, Real& sn) {
if (dy == Real(0)) {
cs = 1.0;
sn = 0.0;
} else if (abs(dy) > abs(dx)) {
Real temp = dx / dy;
sn = 1.0 / sqrt( 1.0 + temp*temp );
cs = temp * sn;
} else {
Real temp = dy / dx;
cs = 1.0 / sqrt( 1.0 + temp*temp );
sn = temp * cs;
}
}
template<class Real>
void apply_plane_rotation (Real& dx, Real& dy, const Real& cs, const Real& sn) {
Real temp = cs * dx + sn * dy;
dy = -sn * dx + cs * dy;
dx = temp;
}
// [verbatim_gmres_cont]
} // namespace details
// [verbatim_gmres_synopsis]
template <class Matrix, class Vector, class Preconditioner,
class SmallMatrix, class SmallVector>
int gmres (const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M,
SmallMatrix &H, const SmallVector& V, const solver_option& sopt = solver_option())
// [verbatim_gmres_synopsis]
// [verbatim_gmres]
{
typedef typename Vector::size_type Size;
typedef typename Vector::float_type Real;
std::string label = (sopt.label != "" ? sopt.label : "gmres");
Size m = sopt.krylov_dimension;
Vector w;
SmallVector s(m+1), cs(m+1), sn(m+1);
Real residue;
Real norm_b = norm(M.solve(b));
Vector r = M.solve(b - A * x);
Real beta = norm(r);
if (sopt.p_err) (*sopt.p_err) << "[" << label << "] # norm_b=" << norm_b << std::endl
<< "[" << label << "] #iteration residue" << std::endl;
if (sopt.absolute_stopping || norm_b == Real(0)) norm_b = 1;
sopt.n_iter = 0;
sopt.residue = norm(r)/norm_b;
if (sopt.residue <= sopt.tol) return 0;
std::vector<Vector> v (m+1);
for (sopt.n_iter = 1; sopt.n_iter <= sopt.max_iter; ) {
v[0] = r/beta;
for (Size i = 0; i < m+1; i++) s(i) = 0; // std::numeric_limits<Float>::max();
s(0) = beta;
for (Size i = 0; i < m && sopt.n_iter <= sopt.max_iter; i++, sopt.n_iter++) {
w = M.solve(A * v[i]);
for (Size k = 0; k <= i; k++) {
H(k, i) = dot(w, v[k]);
w -= H(k, i) * v[k];
}
H(i+1, i) = norm(w);
v[i+1] = w/H(i+1,i);
for (Size k = 0; k < i; k++) {
details::apply_plane_rotation (H(k,i), H(k+1,i), cs(k), sn(k));
}
details::generate_plane_rotation (H(i,i), H(i+1,i), cs(i), sn(i));
details::apply_plane_rotation (H(i,i), H(i+1,i), cs(i), sn(i));
details::apply_plane_rotation (s(i), s(i+1), cs(i), sn(i));
sopt.residue = abs(s(i+1))/norm_b;
if (sopt.p_err) (*sopt.p_err) << "[" << label << "] " << sopt.n_iter << " " << sopt.residue << std::endl;
if (sopt.residue <= sopt.tol) {
details::update (x, i, H, s, v);
return 0;
}
}
details::update (x, m - 1, H, s, v);
r = M.solve(b - A * x);
beta = norm(r);
sopt.residue = beta/norm_b;
if (sopt.p_err) (*sopt.p_err) << "[" << label << "] " << sopt.n_iter << " " << sopt.residue << std::endl;
if (sopt.residue < sopt.tol) return 0;
}
return 1;
}
// [verbatim_gmres]
}// namespace rheolef
# endif // _RHEOLEF_GMRES_H
rheolef::solver_option
see the solver_option page for the full documentation
Definition solver_option.h:155
rheolef::details
Definition compiler_eigen.h:39
rheolef::details::update
void update(Vector &x, Size k, const SmallMatrix &h, const SmallVector &s, Vector2 &v)
Definition gmres.h:132
rheolef::details::apply_plane_rotation
void apply_plane_rotation(Real &dx, Real &dy, const Real &cs, const Real &sn)
Definition gmres.h:160
rheolef::details::generate_plane_rotation
void generate_plane_rotation(const Real &dx, const Real &dy, Real &cs, Real &sn)
Definition gmres.h:145
rheolef
This file is part of Rheolef.
Definition compiler_eigen.h:39
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::norm
T norm(const vec< T, M > &x)
norm(x): see the expression page for the full documentation
Definition vec.h:387
rheolef::dot
rheolef::std::enable_if< details::is_vec_expr_v2_arg< Expr1 >::value &&details::is_vec_expr_v2_arg< Expr2 >::value, promote< Expr1::float_type, Expr2::float_type >::type > type dot const Expr1 expr1, const Expr2 expr2 dot(const Expr1 &expr1, const Expr2 &expr2)
Definition vec_expr_v2.h:415
residue
field residue(Float p, const field &uh)
Definition p_laplacian_post.cc:35
A
Definition diffusion_isotropic.h:25
rheolef::Vector
Definition Vector.h:176
rheolef::Vector::size_type
DATA::size_type size_type
Definition Vector.h:188
M
Expr1::memory_type M
Definition vec_expr_v2.h:416