JCGRC

Solves a real symmetric definite linear system using the Jacobi-preconditioned conjugate gradient method with reverse communication.

Required Arguments

IDO — Flag indicating task to be done. (Input/Output)

On the initial call IDO must be 0. If the routine returns with IDO = 1, then set

Z = A * P, where A is the matrix, and call JCGRC again. If the routine returns with IDO = 2, then the iteration has converged and X contains the solution.

On the initial call IDO must be 0. If the routine returns with IDO = 1, then set

Z = A * P, where A is the matrix, and call JCGRC again. If the routine returns with IDO = 2, then the iteration has converged and X contains the solution.

DIAGNL — Vector of length N containing the diagonal of the matrix. (Input)

Its elements must be all strictly positive or all strictly negative.

Its elements must be all strictly positive or all strictly negative.

X — Array of length N containing the solution. (Input/Output)

On input, X contains the initial guess of the solution. On output, X contains the solution to the system.

On input, X contains the initial guess of the solution. On output, X contains the solution to the system.

P — Array of length N. (Output)

Its use is described under IDO.

Its use is described under IDO.

R — Array of length N. (Input/Output)

On initial input, it contains the right-hand side of the linear system. On output, it contains the residual.

On initial input, it contains the right-hand side of the linear system. On output, it contains the residual.

Z — Array of length N. (Input)

When IDO = 1, it contains AP, where A is the linear system. When IDO = 0, it is ignored. Its use is described under IDO.

When IDO = 1, it contains AP, where A is the linear system. When IDO = 0, it is ignored. Its use is described under IDO.

Optional Arguments

N — Order of the linear system. (Input)

Default: N = size (X,1).

Default: N = size (X,1).

RELERR — Relative error desired. (Input)

Default: RELERR = 1.e-5 for single precision and 1.d-10 for double precision.

Default: RELERR = 1.e-5 for single precision and 1.d-10 for double precision.

ITMAX — Maximum number of iterations allowed. (Input)

Default: ITMAX = 100.

Default: ITMAX = 100.

FORTRAN 90 Interface

Generic: CALL JCGRC (IDO, DIAGNL, X, P, R, Z [, …])

Specific: The specific interface names are S_JCGRC and D_JPCGRC.

FORTRAN 77 Interface

Single: CALL JCGRC (IDO, N, DIAGNL, X, P, R, Z, RELERR, ITMAX)

Double: The double precision name is DJCGRC.

Description

Routine JCGRC solves the symmetric definite linear system Ax = b using the Jacobi conjugate gradient method. This method is described in detail by Golub and Van Loan (1983, Chapter 10), and in Hageman and Young (1981, Chapter 7).

This routine is a special case of the routine PCGRC, with the diagonal of the matrix A used as the preconditioning matrix. For details of the algorithm see PCGRC.

The number of iterations needed depends on the matrix and the error tolerance RELERR. As a rough guide, ITMAX = N is often sufficient when N ≫ 1. See the references for further information.

Comments

1. Workspace may be explicitly provided, if desired, by use of J2GRC/DJ2GRC. The reference is:

CALL J2GRC (IDO, N, DIAGNL, X, P, R, Z, RELERR, ITMAX, TRI, WK, IWK)

The additional arguments are as follows:

TRI — Workspace of length 2 * ITMAX containing a tridiagonal matrix (in band symmetric form) whose largest eigenvalue is approximately the same as the largest eigenvalue of the iteration matrix. The workspace arrays TRI, WK and IWK should not be changed between the initial call with IDO = 0 and JCGRC/DJCGRC returning with IDO = 2.

WK — Workspace of length 5 * ITMAX.

IWK — Workspace of length ITMAX.

2. Informational errors

Type | Code | Description |
---|---|---|

4 | 1 | The diagonal contains a zero. |

4 | 2 | The diagonal elements have different signs. |

4 | 3 | No convergence after ITMAX iterations. |

4 | 4 | The linear system is not definite. |

4 | 5 | The linear system is singular. |

Example

In this example, the solution to a linear system is found. The coefficient matrix A is stored as a full matrix.

USE IMSL_LIBRARIES

INTEGER LDA, N

PARAMETER (LDA=3, N=3)

!

INTEGER IDO, ITMAX

REAL A(LDA,N), B(N), DIAGNL(N), P(N), R(N), X(N), &

Z(N)

! ( 1, -3, 2 )

! A = ( -3, 10, -5 )

! ( 2, -5, 6 )

DATA A/1.0, -3.0, 2.0, -3.0, 10.0, -5.0, 2.0, -5.0, 6.0/

! B = ( 27.0, -78.0, 64.0 )

DATA B/27.0, -78.0, 64.0/

! Set R to right side

CALL SCOPY (N, B, 1, R, 1)

! Initial guess for X is B

CALL SCOPY (N, B, 1, X, 1)

! Copy diagonal of A to DIAGNL

CALL SCOPY (N, A(:, 1), LDA+1, DIAGNL, 1)

! Set parameters

ITMAX = 100

IDO = 0

10 CALL JCGRC (IDO, DIAGNL, X, P, R, Z, ITMAX=ITMAX)

IF (IDO .EQ. 1) THEN

! Set z = Ap

CALL MURRV (A, P, Z)

GO TO 10

END IF

! Print the solution

CALL WRRRN (’Solution’, X)

!

END

Output

Solution

1 1.001

2 -4.000

3 7.000