linSolGenMinResidual¶

Solves a linear system \(Ax = b\) using the restarted generalized minimum residual (GMRES) method.

Synopsis¶

linSolGenMinResidual (amultp, b)

Required Arguments¶

amultp (p, z): User-supplied function which computes \(z = Ap\).
float b (Input): Vector of length n containing the right-hand side.

Return Value¶

The solution x of the linear system \(Ax = b\). If no solution was computed, then None is returned.

Optional Arguments¶

maxIter (Input/Output)

An integer, initially set to the maximum number of GMRES iterations allowed. On exit, the number of iterations used is returned.

Default: maxIter = 1000

relErr (Input)

The relative error desired.

Default: relErr = sqrt(machine(4))

precond (r, z) (Input)

User supplied function which sets \(z = M^{-1}r\), where M is the preconditioning matrix.

maxKrylovSubspaceDim (Input)

The maximum Krylov subspace dimension, i.e., the maximum allowable number of GMRES iterations allowed before restarting.

Default: maxKrylovSubspaceDim = iMin(n, 20)

householderReorthog

Perform orthogonalization by Householder transformations, replacing the Gram-Schmidt process.

Description¶

The function linSolGenMinResidual, based on the FORTRAN subroutine GMRES by H.F. Walker, solves the linear system \(Ax = b\) using the GMRES method. This method is described in detail by Saad and Schultz (1986) and Walker (1988).

The GMRES method begins with an approximate solution \(x_0\) and an initial residual \(r_0 = b - Ax_0\). At iteration m, a correction \(z_m\) is determined in the Krylov subspace

\[\kappa^m (v) = \texttt{span} (v, Av, \ldots, A^{m-1}v)\]

\(v = r_0\) which solves the least-squares problem

\[\left(\mathit{z} \in \min_{\kappa_m} \left(r_0\right)\right) \phantom{.....} \| b - A \left(x_0 + z\right) \|_2\]

Then at iteration m, \(x_m = x_0 + z_m\).

Orthogonalization by Householder transformations requires less storage but more arithmetic than Gram-Schmidt. However, Walker (1988) reports numerical experiments which suggest the Householder approach is more stable, especially as the limits of residual reduction are reached.

Examples¶

Example 1¶

As an example, consider the following matrix:

\[\begin{split}A = \begin{bmatrix} 10 & 0 & 0 & 0 & 0 & 0 \\ 0 & 10 & -3 & -1 & 0 & 0 \\ 0 & 0 & 15 & 0 & 0 & 0 \\ -2 & 0 & 0 & 10 & -1 & 0 \\ -1 & 0 & 0 & -5 & 1 & -3 \\ -1 & -2 & 0 & 0 & 0 & 6 \\ \end{bmatrix}\end{split}\]

Let \(x^T = (1, 2, 3, 4, 5, 6)\) so that \(Ax = (10, 7, 45, 33, -34, 31)^T\). The function matMulRectCoordinate is used to form the product Ax.

from numpy import *
from pyimsl.math.linSolGenMinResidual import linSolGenMinResidual
from pyimsl.math.matMulRectCoordinate import matMulRectCoordinate
from pyimsl.math.writeMatrix import writeMatrix


def amultp(p, z):
    a = [[0, 0, 10.0],
         [1, 1, 10.0],
         [1, 2, -3.0],
         [1, 3, -1.0],
         [2, 2, 15.0],
         [3, 0, -2.0],
         [3, 3, 10.0],
         [3, 4, -1.0],
         [4, 0, -1.0],
         [4, 3, -5.0],
         [4, 4, 1.0],
         [4, 5, -3.0],
         [5, 0, -1.0],
         [5, 1, -2.0],
         [5, 5, 6.0]]

    n = 6
    nz = 15
    b = empty((n), dtype=double)
    for i in range(0, n):
        b[i] = p[i]
    ztmp = matMulRectCoordinate("A*x",
                                aMatrix={'nrowa': n, 'ncola': n, 'a': a},
                                xVector=b)
    for j in range(0, len(ztmp)):
        z[j] = ztmp[j]


b = array([10.0, 7.0, 45.0, 33.0, -34.0, 31.0])
x = linSolGenMinResidual(amultp, b)
writeMatrix("Solution, x, to Ax = b", x)

Output¶

 
                           Solution, x, to Ax = b
          1            2            3            4            5            6
          1            2            3            4            5            6

Example 2¶

In this example, the same system given in the first example is solved. This time a preconditioner is provided. The preconditioned matrix is chosen as the diagonal of A.

from __future__ import print_function
from numpy import *
from pyimsl.math.linSolGenMinResidual import linSolGenMinResidual
from pyimsl.math.matMulRectCoordinate import matMulRectCoordinate
from pyimsl.math.writeMatrix import writeMatrix


def amultp(p, z):
    a = [[0, 0, 10.0],
         [1, 1, 10.0],
         [1, 2, -3.0],
         [1, 3, -1.0],
         [2, 2, 15.0],
         [3, 0, -2.0],
         [3, 3, 10.0],
         [3, 4, -1.0],
         [4, 0, -1.0],
         [4, 3, -5.0],
         [4, 4, 1.0],
         [4, 5, -3.0],
         [5, 0, -1.0],
         [5, 1, -2.0],
         [5, 5, 6.0]]

    n = 6
    nz = 15
    b = empty((n), dtype=double)
    for i in range(0, n):
        b[i] = p[i]
    ztmp = matMulRectCoordinate("A*x",
                                aMatrix={'nrowa': n, 'ncola': n, 'a': a},
                                xVector=b)
    for j in range(0, len(ztmp)):
        z[j] = ztmp[j]

# Solve Mz = r


def precond(r, z):
    diagonal_inverse = array([0.1, 0.1, 1.0 / 15.0, 0.1, 1.0, 1.0 / 6.0])
    n = 6
    for i in range(0, n):
        z[i] = diagonal_inverse[i] * r[i]


b = array([10.0, 7.0, 45.0, 33.0, -34.0, 31.0])
n = 6
maxit = [1000]
x = linSolGenMinResidual(amultp, b,
                         maxIter=maxit,
                         precond=precond)

writeMatrix("Solution, x, to Ax = b", x)
print("\nNumber of iterations taken = %d\n" % maxit[0])

Output¶

Number of iterations taken = 5

 
                           Solution, x, to Ax = b
          1            2            3            4            5            6
          1            2            3            4            5            6

Fatal Errors¶

IMSL_STOP_USER_FCN

Request from user supplied function to stop algorithm.

User flag = “#”.