IMSL_SP

IMSL_SP_CG

The IMSL_SP_CG function solves a real symmetric definite linear system using a

conjugate gradient method. A preconditioner can be supplied by using keywords.

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

The IMSL_SP_CG function solves the symmetric definite linear system Ax = b using the conjugate gradient method with optional preconditioning. This method is described in detail by Golub and Van Loan (1983, chapter 10), and in Hageman and Young (1981, chapter 7).

The preconditioning matrix M, is a matrix that approximates A, and for which the linear system Mz = r is easy to solve. These two properties are in conflict; balancing them is a topic of much current research. In the default usage of IMSL_SP_CG, M = I. If the keyword JACOBI is selected, M is set to the diagonal of A.

The number of iterations needed depends on the matrix and the error tolerance. As a rough guide:

Itmax = √n

for

n >> 1

See the academic references for details.

Let M be the preconditioning matrix, let b, p, r, x, and z be vectors and let t be the desired relative error. Then the algorithm used is as follows:

Example

This example finds the solution to a linear system. The coefficient matrix is stored as a full matrix.

FUNCTION Amultp, p

; Since A is in dense form, we use the # operator to perform the

; matrix-vector product. The common block us used to hold A.

  COMMON Cg_comm1, a

  RETURN, a#p

END

Pro CG_EX1

  COMMON Cg_comm1, a

  a = TRANSPOSE([[ 1, -3, 2], [-3, 10, -5], [ 2, -5, 6]])

  b = [27, -78, 64]

  x = IMSL_SP_CG('amultp', b)

  ; Use IMSL_SP_CG to compute the solution, then output

  ; the result.

  PM, x, title = 'Solution to Ax = b'

END

; Output of this procedure:

Solution to Ax = b

  1.0000000

  -4.0000000

  7.0000000

Syntax

Result = IMSL_SP_CG(Amultp, B [, /DOUBLE] [, ITMAX=value] [, JACOBI=vector] [, PRECOND=value] [, REL_ERR=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 which 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.

ITMAX (optional)

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

JACOBI (optional)

If present, use the Jacobi preconditioner, that is, M = diag(A). The supplied vector Jacobi should be set so that JACOBI(i) = A_i,i.

PRECOND (optional)

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

REL_ERR (optional)

Initially set to relative error desired. On exit, the computed relative error is returned. Default: SQRT(machine precision)

Version History

6.4	Introduced

Module	Math&Stats

Version	9.0