File:GMRESm.f90: Difference between revisions

Revision as of 05:05, 13 December 2016

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

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 minimises the error $|A\vec{x}-\vec{b}|$ subject to the additional constraint $|\vec{x}|<\delta$.  This 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 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 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.
 == 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.  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$.