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.
DIAGNL — Vector of length N containing the diagonal of the matrix. (Input)
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.
P — Array of length N. (Output)
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.
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.
Optional Arguments
N — Order of the linear system. (Input)
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.
ITMAX — Maximum number of iterations allowed. (Input)
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
Solution
1 1.001
2 -4.000
3 7.000
