File:GMRESm.f90

$ \renewcommand{\vec}[1]{ {\bf #1} } \newcommand{\bnabla}{ \vec{\nabla} } \newcommand{\Rey}{Re} \def\vechat#1{ \hat{ \vec{#1} } } \def\mat#1{#1} $

The GMRES(m) Method

This implements the classic GMRES(m) method for solving the system \[A\vec{x}=\vec{b}\] for $\vec{x}$. This implementation introduces the extra constraint $|\vec{x}|<\delta$, which may be ignored by putting $\delta$ large.

The main advantage of the method is that it only requires calculations of multiplies by $A$ for a given $\vec{x}$ -- it does not need to know $A$ itself. This means that $A$ need not even be stored, and could correspond to a very complex linear 'action' on $\vec{x}$, e.g. a time integral with initial condition $\vec{x}$. The method seeks eigenvectors in $\mathrm{span}\{\vec{x},\,A\vec{x},\,A^2\vec{x},...\}$, but uses Gram-Schmidt orthogonalisation to improve the suitability of the basis. The set of orthogonalised vectors is called the Krylov-subspace, and m is the maximum number of vectors stored.

Whereas m is traditionally a small number, e.g. 3 or 4, the extra constraint renders restarts difficult. If the constraint is important, then m should be chosen sufficiently large to solve within m iterations.

The code

To download, click the link above. The routine needs to be passed an external functions that calculates dot products, and two external subroutines, one that calculates the action of multiplication by $A$, and another that replaces a vector with the action of a preconditioner matrix. The latter may simply be an empty subroutine if no preconditioner is required.