# File:GMRESm.f90

GMRESm.f90 (file size: 5 KB, MIME type: text/plain)

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

FORTRAN: link above.

MATLAB / GNU Octave: File:GMRESm.m

## 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 minimises the error $\|A\vec{x}-\vec{b}\|$ subject to the additional constraint $\|\vec{x}\|<\delta$. (This constraint may be ignored by supplying $\delta<0$ in the implementation.)

The main advantage of the GMRES 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}$. For a given starting vector $\vec{x}_0$, the method seeks solutions for $\vec{x}$ in $\mathrm{span}\{\vec{x}_0,\,A\vec{x}_0,\,A^2\vec{x}_0,...\}$, but uses Gram-Schmidt orthogonalisation to improve the numerical suitability of this basis set. 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 additional constraint renders restarts difficult. If the constraint is important, then m must be chosen sufficiently large to solve to the desired accuracy within m iterations.

**Preconditioning**.

The implementations above can be supplied a preconditioner routine. (This can be avoided if combined with timestepping; see the remarks at File:Arnoldi.f).

GMRES is likely to find it easier to solve $M^{-1}A\,x=M^{-1}b$ than the original system, if $M^{-1}$ is an approximate inverse for $A$. For example, if $A$ is dominated by its diagonal elements, we might take $M$ to be the banded matrix consisting of the diagonal and the first sub- and super-diagonals of $A$. Each GMRES iteration applied to the modified system now requires a muliplication by $A$ then by $M^{-1}$. This is fine, as it is quick and easy to solve $Mx'=x$ for $x'$ for a banded matrix $M$. Like $A$, we don't need to know the matrix $M^{-1}$ itself, only the result of multiplication by these matrices.

## The code

To download, click the link above. The code uses the LAPACK package.

This constraint $\|\vec{x}\|<\delta$ may be ignored by supplying negative '`del`'.

In addition to scalar and array variables, the routine needs to be passed

- an external function that calculates dot products,
- an external subroutine that calculates the result of multiplication by $A$,
- an external subroutine that replaces a given vector $\vec{x}$ with the solution $\vec{x}'$ of the system $M\vec{x}'=\vec{x}$. This may simply be an empty subroutine if no preconditioner is required, i.e. $M=I$.

The functions above may require auxiliary data in addition to $\vec{x}$ or $\vec{\delta x}$. Place this data in a module and access via '`use mymodule`' in the function/subroutine.

## Parallel use

It is NOT necessary to edit this code for parallel (MPI) use:

- let each thread pass its subsection for the vector $\vec{x}$,
- make the dot product function
`mpi_allreduce`the result of the dot product. - to avoid multiple outputs to the terminal, set
`info=1`on rank 0 and`info=0`for the other ranks.

## File history

Click on a date/time to view the file as it appeared at that time.

Date/Time | Dimensions | User | Comment | |
---|---|---|---|---|

current | 03:48, 26 June 2019 | (5 KB) | Apwillis (talk | contribs) | Minor edits in comments only. |

03:30, 13 December 2016 | (5 KB) | Apwillis (talk | contribs) | Solve Ax=b for x, subject to constraint |x|<delta. |

You cannot overwrite this file.

## File usage

The following 6 pages link to this file: