.TH "gmres" 5rheolef "Sat Mar 13 2021" "Version 7.1" "rheolef" \" -*- nroff -*-
.ad l
.nh
.SH NAME
gmres \- generalized minimum residual algorithm (rheolef-7\&.1)
.PP
 
.SH "SYNOPSIS"
.PP
.PP
.nf
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())
.fi
.PP
 
.SH "EXAMPLE"
.PP
.PP
.nf
    solver_option sopt;
    sopt.tol = 1e-7;
    sopt.max_iter = 100;
    size_t m = sopt.krylov_dimension = 6;
    Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> H(m+1,m+1);
    Eigen::Matrix<T,Eigen::Dynamic,1>              V(m);
    int status = gmres (A, x, b, ilut(a), H, V, sopt);
.fi
.PP
 
.SH "DESCRIPTION"
.PP
This function solves the \fIunsymmetric\fP linear system \fCA*x=b\fP with the generalized minimum residual algorithm\&. The \fCgmres\fP function follows the algorithm described on p\&. 20 in 
.PP
.nf
    R. Barrett, M. Berry, T. F. Chan, J. Demmel, J. Donato,
    J. Dongarra, V. Eijkhout, R. Pozo, C. Romine and H. Van der Vorst,
    Templates for the solution of linear systems: building blocks for iterative methods, 
    SIAM, 1994.

.fi
.PP
 The fourth argument of \fCgmres\fP is a preconditionner: here, the \fBilut(5)\fP one is used, for simplicity\&.
.PP
Next, \fCH\fP specifies a matrix to hold the coefficients of the upper Hessenberg matrix constructed by the \fCgmres\fP iterations, \fCm\fP specifies the number of iterations for each restart\&. We have used here the \fCeigen\fP dense matrix and vector types for the \fCH\fP and \fCV\fP vectors, with sizes related to the Krylov space dimension \fCm\fP\&. Finally, the \fBsolver_option(4)\fP variable \fCsopt\fP transmits the stopping criterion with \fCsopt\&.tol\fP and \fCsopt\&.max_iter\fP\&.
.PP
On return, the \fCsopt\&.residue\fP and \fCsopt\&.n_iter\fP indicate the reached residue and the number of iterations effectively performed\&. The return \fCstatus\fP is zero when the prescribed tolerance \fCtol\fP has been obtained, and non-zero otherwise\&. Also, the \fCx\fP variable contains the approximate solution\&. See also the \fBsolver_option(4)\fP for more controls upon the stopping criterion\&.
.SH "REMARKS"
.PP
\fCgmres\fP requires two matrices as input: \fCA\fP and \fCH\fP\&. The \fCA\fP matrix, which will typically be \fIsparse\fP, corresponds to the matrix involved in the linear system A*x=b\&. Conversely, the \fCH\fP matrix, which will typically be \fIdense\fP, corresponds to the upper Hessenberg matrix that is constructed during the \fCgmres\fP iterations\&. Within \fCgmres\fP, \fCH\fP is used in a different way than \fCA\fP, so its class must supply different functionality\&. That is, \fCA\fP is only accessed though its matrix-vector and transpose-matrix-vector multiplication functions\&. On the other hand, \fCgmres\fP solves a dense upper triangular linear system of equations on \fCH\fP\&. Therefore, the class to which \fCH\fP belongs must provide \fCH(i,j)\fP accessors\&.
.PP
It is important to remember that we use the convention that indices are 0-based\&. That is \fCH(0,0)\fP is the first component of the matrix\&. 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\&.
.SH "IMPLEMENTATION"
.PP
This documentation has been generated from file linalg/lib/gmres\&.h
.PP
The present template implementation is inspired from the \fCIML++ 1\&.2\fP iterative method library, http://math.nist.gov/iml++
.PP
.PP
.nf
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())
.fi
.PP
.PP
.nf
{
  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;
}
.fi
.PP
.PP
.nf
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;
}
.fi
.PP

.SH AUTHOR
Pierre  Saramito  <Pierre.Saramito@imag.fr>
.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.