File:GMRESm.f90: Difference between revisions

Revision as of 06:50, 13 November 2018

$ \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 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 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 added constraint renders restarts difficult. If the constraint is important, then m should be chosen sufficiently large to solve within m iterations.

GMRES is closely related to the Arnoldi method. See the remarks at File:Arnoldi.f on improved suitability via timestepping or exponentiation of the matrix.

The code

To download, click the link above. The Lapack package is also required.

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

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	04:48, 26 June 2019	(5 KB)	Apwillis (talk \| contribs)	Minor edits in comments only.
	04: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 5 pages use this file:

@@ Line 3: / Line 3: @@
 == 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.
+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 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.
+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 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 extra constraint renders restarts difficult.  If the constraint is important, then m should be chosen sufficiently large to solve within m iterations.
+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.
+GMRES is closely related to the Arnoldi method.  See the remarks at [[File:Arnoldi.f]] on improved suitability via timestepping or exponentiation of the matrix.
 == 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.
+To download, click the link above.  The Lapack package is also required.
+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 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$.
+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.