gmres(5rheolef)

📦
version 7.2 · 4 versions
, in debian trixie, ubuntu noble, ubuntu questing

🌸 NAME 🌸

okay so bestie, you just opened gmres(5rheolef) — that's the rheolef manual, section 5, and this entire page is about gmres 💅

gmres — generalized minimum residual algorithm (rheolef-7.2) ✨

🦋 SYNOPSIS 🦋

here's the signature, verbatim, no edits, no cap:

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())

🎀 EXAMPLE 🎀

this is how it looks in the wild, obsessed:

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);

💖 DESCRIPTION 💖

so this function solves the unsymmetric linear system A*x=b with the generalized minimum residual algorithm 🌈 and it follows the algorithm described on p. 20 in this absolute classic:

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.

the fourth argument of gmres is a preconditioner, and here the ilut(5) one is used, for simplicity 💅

next up, H specifies a matrix to hold the coefficients of the upper Hessenberg matrix that gets constructed by the gmres iterations, and m specifies the number of iterations for each restart ✨ we've used the eigen dense matrix and vector types here for the H and V vectors, with sizes related to the Krylov space dimension m 🦋 and finally, the solver_option(4) variable sopt transmits the stopping criterion via sopt.tol and sopt.max_iter.

on return, sopt.residue and sopt.n_iter tell you the reached residue and the number of iterations effectively performed 👑 the return status is zero when the prescribed tolerance tol has been obtained, and non-zero otherwise. also, the x variable contains the approximate solution 🫶 see also the solver_option(4) for more controls upon the stopping criterion.

🌟 REMARKS 🌟

gmres requires two matrices as input: A and H 😤

  • the A matrix, which will typically be sparse, corresponds to the matrix involved in the linear system A*x=b.
  • conversely, the H matrix, which will typically be dense, corresponds to the upper Hessenberg matrix 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 through 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) accessors 👑

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. 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 ✨

🧚 IMPLEMENTATION 🧚

this documentation has been generated from file linalg/lib/gmres.h 📄

the present template implementation is inspired from the IML++ 1.2 iterative method library, http://math.nist.gov/iml++ 🌈

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())

{
  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;
}

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;
}

👑 AUTHOR 👑

Pierre Saramito <Pierre.Saramito@imag.fr>

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.

okay that's the whole page bestie, you ate 💅 this was gmres from the rheolef manual, Version 7.2 🎀 go forth and solve those unsymmetric systems, mwah 💖✨🦋

📦 where this page comes from