The IMSL_SP_GMRES function solves a linear system Ax = b using the restarted generalized minimum residual (GMRES) method.

This routine requires an IDL Advanced Math and Stats license. For more information, contact your sales or technical support representative.

The IMSL_SP_GMRES, function based on the FORTRAN subroutine GMRESD 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 x0 and an initial residual r0 = b – Ax0. At iteration m, a correction zm is determined in the Krylov subspace:

κm(v) = span(v, Av, ..., Am–1v)

v = r0 which solves the least squares problem:

Then at iteration m, xm = x0 + zm.

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.

Example

This example finds the solution to a linear system.

The coefficient matrix is stored in coordinate format. Consider the following matrix: Let xT = (1, 2, 3, 4, 5, 6) so that Ax = (10, 7, 45, 33, –34, 31)T. The number of nonzeros in A is 15.

FUNCTION Amultp, p
; This function uses IMSL_SP_MVMUL to multiply a sparse
; matrix stored in coordinate storage mode and a vector.
; The common block holds the sparse matrix.
COMMON Gmres_ex1, nrows, ncols, a
RETURN, IMSL_SP_MVMUL(nrows, ncols, a, p)
END
 
PRO Gmres1
; This procedure defines the sparse matrix A stored in coordinate
; storage mode, and then calls IMSL_SP_GMRES to compute the
; solution to Ax = b.
COMMON Gmres_ex1, nrows, ncols, a
; Initialize sparse matrix structure variables
@imsl_init
 
A = REPLICATE(imsl_f_sp_elem, 15)
a(*).row = [0, 1, 1, 1, 2, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5]
a(*).col = [0, 1, 2, 3, 2, 0, 3, 4, 0, 3, 4, 5, 0, 1, 5]
a(*).val = [10, 10, -3, -1, 15, -2, 10, -1, -1, -5, $
1, -3, -1, -2, 6]
nrows = 6
ncols = 6
b = [10, 7, 45, 33, -34, 31]
itmax = 10
; Itmax is input/output.
x = IMSL_SP_GMRES('amultp', b, Itmax = itmax)
pm, x, title = 'Solution to Ax = b'
pm, itmax, title = 'Number of iterations'
END
; Output of this procedure:
Solution to Ax = b
  1.0000000
  2.0000000
  3.0000000
  4.0000000
  5.0000000
  6.0000000
Number of iterations
  6

Syntax

Result = IMSL_SP_GMRES(Amultp, B [, /DOUBLE] [, HH_REORTH=value] [, ITMAX=value] [, MAX_KRYLOV=value] [, PRECOND=value] [, TOLERANCE=value])

Return Value

A one-dimensional array containing the solution of the linear system Ax = b.

Arguments

Amultp

Scalar string specifying a user supplied function that computes z = Ap. The function accepts the argument p, and returns the vector Ap.

B

One-dimensional matrix containing the right-hand side.

Keywords

DOUBLE (optional)

If present and nonzero, double precision is used.

HH_REORTH (optional)

If present and nonzero, perform orthogonalization by Householder transformations, replacing the Gram-Schmidt process.

ITMAX (optional)

Initially set to the maximum number of GMRES iterations allowed. On exit, the number of iterations used is returned. Default: 1000

MAX_KRYLOV (optional)

The maximum Krylov subspace dimension, that is, the maximum allowable number of GMRES iterations allowed before restarting. Default: Minimum(N_ELEMENTS(b), 20).

PRECOND (optional)

Scalar sting specifying a user supplied function which sets z = M–1r, where M is the preconditioning matrix.

TOLERANCE (optional)

The algorithm attempts to generate x such that: where t = TOLERANCE. Default: SQRT(machine precision).

Version History

6.4

Introduced

See Also

IMSL_SP_BDPDFAC, IMSL_SP_BDPDSOL, IMSL_SP_BDSOL, IMSL_SP_CG, IMSL_SP_LUFAC, IMSL_SP_LUSOL, IMSL_SP_MVMUL, IMSL_SP_PDFAC, IMSL_SP_PDSOL