.\"
.de Id
..
.de Sp
.if n .sp
.if t .sp 0.4
..
.TH gmres 4rheolef "rheolef-7.0" "rheolef-7.0" "rheolef-7.0"
.\" label: /*Class:gmres
.SH NAME
\fBgmres\fP -- generalized minimum residual method
.\" skip: @findex qmr
.\" skip: @cindex generalized minimum residual method
.\" skip: @cindex iterative solver
.SH SYNOPSIS

.\" begin_example
.Sp
.nf
  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);
.Sp
.fi
.\" end_example
.SH EXAMPLE

The simplest call to \fBgmres\fP has the folling form:
.\" begin_example
.Sp
.nf
  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);
.Sp
.fi
.\" end_example
.SH DESCRIPTION

\fBgmres\fP solves the unsymmetric linear system Ax = b
using the generalized minimum residual method.
.PP
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:
.\" begin table
.\" start item
.TP
.B x
approximate solution to Ax = b

.\" start item
.TP
.B sopt.n_iter
the number of iterations performed before the tolerance was reached

.\" start item
.TP
.B sopt.residue
the residual after the final iteration
.\" end table
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 \fBgmres\fP iterations, \fBm\fP specifies the number of iterations
for each restart.
.PP
\fBgmres\fP 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
\fBgmres\fP iterations. Within \fBgmres\fP, 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, \fBgmres\fP 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.
.PP
.SH 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.
.PP
\fBgmres\fP is an iterative template routine.
.PP
\fBgmres\fP follows the algorithm described on p. 20 in
\fITemplates for the solution of linear systems: building blocks for iterative methods\fP,
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,
\fBftp.netlib.org/templates/templates.ps\fP.
.PP
The present implementation is inspired from \fBIML++ 1.2\fP iterative method library,
\fBhttp://math.nist.gov/iml++\fP.
.PP
.\" skip start:AUTHOR:



.\" skip start:DATE:



.\" skip start:METHODS:
.\" END
.SH IMPLEMENTATION
.\" begin_example
.Sp
.nf
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;
}
.Sp
.fi
.\" end_example

.\" LENGTH = 1
.SH SEE ALSO
solver_option(2)
.SH 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.