File:GMRESm.f90: Difference between revisions

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

The 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 extra constraint means that restarts are not straight forward. m should be chosen sufficiently large to solve within the specified residual error, $|Ax-b|/|b|$<res.

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.