Provided by: librheolef-dev_7.0-3_amd64 
NAME
gmres -- generalized minimum residual method
SYNOPSIS
template <class Matrix, class Vector, class Preconditioner,
class SmallMatrix, class SmallVector, class Real, class Size>
int gmres (const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M,
SmallMatrix &H, const SmallVector& dummy,
const solver_option& sopt);
EXAMPLE
The simplest call to gmres has the folling form:
solver_option sopt;
sopt.tol = 1e-7;
sopt.max_iter = 100;
size_t m = sopt.krylov_dimension;
boost::numeric::ublas::matrix<T> H(m+1,m+1);
vec<T,sequential> dummy(m,0.0);
int status = gmres (a, x, b, ic0(a), H, dummy, sopt);
DESCRIPTION
gmres solves the unsymmetric linear system Ax = b using the generalized minimum residual
method.
The return value indicates convergence within max_iter (input) iterations (0), or no
convergence within max_iter iterations (1). Upon successful return, output arguments have
the following values:
x approximate solution to Ax = b
sopt.n_iter
the number of iterations performed before the tolerance was reached
sopt.residue
the residual after the final iteration See also the solver_option(2). In addition,
M specifies a preconditioner, H specifies a matrix to hold the coefficients of the
upper Hessenberg matrix constructed by the gmres iterations, m specifies the number
of iterations for each restart.
gmres requires two matrices as input, A and H. The matrix A, which will typically be a
sparse matrix) corresponds to the matrix in the linear system Ax=b. The matrix H, which
will be typically a dense matrix, corresponds to the upper Hessenberg matrix H that is
constructed during the gmres iterations. Within gmres, H is used in a different way than
A, so its class must supply different functionality. That is, A is only accessed though
its matrix-vector and transpose-matrix-vector multiplication functions. On the other
hand, gmres solves a dense upper triangular linear system of equations on H. Therefore,
the class to which H belongs must provide H(i,j) operator for element access.
NOTE
It is important to remember that we use the convention that indices are 0-based. That is
H(0,0) is the first component of the matrix H. Also, the type of the matrix must be
compatible with the type of single vector entry. That is, operations such as H(i,j)*x(j)
must be able to be carried out.
gmres is an iterative template routine.
gmres follows the algorithm described on p. 20 in Templates for the solution of linear
systems: building blocks for iterative methods, 2nd Edition, R. Barrett, M. Berry, T. F.
Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo, C. Romine, H. Van der
Vorst, SIAM, 1994, ftp.netlib.org/templates/templates.ps.
The present implementation is inspired from IML++ 1.2 iterative method library,
http://math.nist.gov/iml++.
IMPLEMENTATION
namespace details {
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;
}
} // namespace details
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&, const solver_option& sopt = solver_option())
{
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 * (1.0 / beta); // ??? r / beta
s = 0.0;
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 * (1.0 / H(i+1, i)); // ??? 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;
}
SEE ALSO
solver_option(2)
COPYRIGHT
Copyright (C) 2000-2018 Pierre Saramito <Pierre.Saramito@imag.fr> GPLv3+: GNU GPL version
3 or later <http://gnu.org/licenses/gpl.html>. This is free software: you are free to
change and redistribute it. There is NO WARRANTY, to the extent permitted by law.