File:GMRESm.f90: Difference between revisions

From openpipeflow.org
Jump to navigation Jump to search
(Solve Ax=b for x, subject to constraint |x|<delta.)
 
mNo edit summary
 
(29 intermediate revisions by the same user not shown)
Line 1: Line 1:
Solve Ax=b for x, subject to constraint |x|<delta.
{{latexPreamble}}
 
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 '<tt>del</tt>'.
 
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 '<tt>use mymodule</tt>' 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 <tt>mpi_allreduce</tt> the result of the dot product.
* to avoid multiple outputs to the terminal, set <tt>info=1</tt> on rank 0 and <tt>info=0</tt> for the other ranks.

Latest revision as of 06:28, 31 October 2019

$ \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/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.