Solves a real symmetric definite linear system using the Jacobi-preconditioned conjugate gradient method with reverse communication.
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.
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.
Generic: CALL JCGRC (IDO, DIAGNL, X, P, R, Z [,…])
Specific: The specific interface names are S_JCGRC and D_JPCGRC.
Single: CALL JCGRC (IDO, N, DIAGNL, X, P, R, Z, RELERR, ITMAX)
Double: The double precision name is DJCGRC.
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.
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
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.
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
|
Visual Numerics, Inc. PHONE: 713.784.3131 FAX:713.781.9260 |