# IMSP_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.

Note: This routine requires an IDL Analyst license. For more information, contact your Exelis VIS 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

If present and nonzero, double precision is used.

### ITMAX

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

### JACOBI

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

### PRECOND

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

### REL_ERR

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

## Version History

6.4 | Introduced

##