File:GMRESm.f90

From openpipeflow.org
Revision as of 01:30, 16 December 2016 by Apwillis (talk | contribs)
Jump to navigation Jump to search

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

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 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 added 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 Lapack package is also required.

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 action of multiplication by $A$,
  • an external subroutine that replaces a vector $\vec{x}$ with the solution of $M\vec{x}'=\vec{x}$ for $\vec{x}'$, where $M$ is a preconditioner matrix. This may simply be an empty subroutine if no preconditioner is required, i.e. $M=I$.

Parallel use

No special adaptations are required for parallel (MPI) use -- let each thread pass its subsection for the vector $\vec{x}$, and let the dot product function allreduce the result of the dot product.

File history

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

Date/TimeDimensionsUserComment
current03: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.