File:GMRESm.f90: Difference between revisions

From openpipeflow.org
Jump to navigation Jump to search
mNo edit summary
Line 11: Line 11:
== The code ==
== The code ==


To download, click the link above.  In addition to scalar and array variables, the routine needs to be passed
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 function that calculates dot products,
* an external subroutine that calculates the action of multiplication by $A$,
* an external subroutine that calculates the action of multiplication by $A$,

Revision as of 01:29, 16 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. 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.