PCGRC
Solves a real symmetric definite linear system using a 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 = AP, where A is the matrix, and call PCGRC again. If the routine returns with IDO = 2, then set Z to the solution of the system MZ = R, where M is the preconditioning matrix, and call PCGRC again. If the routine returns with IDO = 3, then the iteration has converged and X contains the solution.
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 = 2, it contains the solution of MZ = R, where M is the preconditioning matrix. 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 = N.
FORTRAN 90 Interface
Generic: CALL PCGRC (IDO, X, P, R, Z [, …])
Specific: The specific interface names are S_PCGRC and D_PCGRC.
FORTRAN 77 Interface
Single: CALL PCGRC (IDO, N, X, P, R, Z, RELERR, ITMAX)
Double: The double precision name is DPCGRC.
Description
Routine PCGRC solves the symmetric definite linear system Ax = b using the preconditioned 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).
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.
The number of iterations needed depends on the matrix and the error tolerance RELERR. As a rough guide, ITMAX = N1/2 is often sufficient when N >> 1. See the references for further information.
Let M be the preconditioning matrix, let b, p, r, x and z be vectors and let τ be the desired relative error. Then the algorithm used is as follows.
λ = −1
p0 = x0
r1 = b − Ap
For k = 1, …, itmax
zk = M-1rk
If k = 1 then
βk = 1
pk = zk
Else
End if
zk = Ap
αk = zk-1Trk-1/zkTpk
xk = xk + αkpk
rk = rk − αkzk
If (∥zk∥2 ≤ τ(1 − λ)∥xk∥2) Then
Recompute λ
If (∥zk∥2 ≤ τ(1 − λ)∥xk∥2) Exit
End if
End loop
Here λ is an estimate of λmax(G), the largest eigenvalue of the iteration matrix G = I − M-1 A. The stopping criterion is based on the result (Hageman and Young, 1981, pages 148−151)
Where
It is known that
where the Tn are the symmetric, tridiagonal matrices
with
and
The largest eigenvalue of Tk is found using the routine EVASB. Usually this eigenvalue computation is needed for only a few of the iterations.
Comments
1. Workspace may be explicitly provided, if desired, by use of P2GRC/DP2GRC. The reference is:
CALL P2GRC (IDO, N, 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 PCGRC/DPCGRC returning with IDO = 3.
WK — Workspace of length 5 * ITMAX.
IWK — Workspace of length ITMAX.
2. Informational errors
|
Type |
Code |
Description |
|
4 |
1 |
The preconditioning matrix is singular. |
|
4 |
2 |
The preconditioning matrix is not definite. |
|
4 |
3 |
The linear system is not definite. |
|
4 |
4 |
The linear system is singular. |
|
4 |
5 |
No convergence after ITMAX iterations. |
Examples
In this example, the solution to a linear system is found. The coefficient matrix A is stored as a full matrix. The preconditioning matrix is the diagonal of A. This is called the Jacobi preconditioner. It is also used by the IMSL routine JCGRC.
USE PCGRC_INT
USE MURRV_INT
USE WRRRN_INT
USE SCOPY_INT
INTEGER LDA, N
PARAMETER (N=3, LDA=N)
!
INTEGER IDO, ITMAX, J
REAL A(LDA,N), B(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)
!
ITMAX = 100
IDO = 0
10 CALL PCGRC (IDO, X, P, R, Z, ITMAX=ITMAX)
IF (IDO .EQ. 1) THEN
! Set z = Ap
CALL MURRV (A, P, Z)
GO TO 10
ELSE IF (IDO .EQ. 2) THEN
! Use diagonal of A as the
! preconditioning matrix M
! and set z = inv(M)*r
DO 20 J=1, N
Z(J) = R(J)/A(J,J)
20 CONTINUE
GO TO 10
END IF
! Print the solution
CALL WRRRN (’Solution’, X)
!
END
Solution
1 1.001
2 -4.000
3 7.000
In this example, a more complicated preconditioner is used to find the solution of a linear system which occurs in a finite-difference solution of Laplace’s equation on a 4 × 4 grid. The matrix is
The preconditioning matrix M is the symmetric tridiagonal part of A,
Note that M, called PRECND in the program, is factored once.
USE IMSL_LIBRARIES
INTEGER LDA, LDPRE, N, NCODA, NCOPRE
PARAMETER (N=9, NCODA=3, NCOPRE=1, LDA=2*NCODA+1,&
LDPRE=NCOPRE+1)
!
INTEGER IDO, ITMAX
REAL A(LDA,N), P(N), PRECND(LDPRE,N), PREFAC(LDPRE,N),&
R(N), RCOND, RELERR, X(N), Z(N)
! Set A in band form
DATA A/3*0.0, 4.0, -1.0, 0.0, -1.0, 2*0.0, -1.0, 4.0, -1.0, 0.0,&
-1.0, 2*0.0, -1.0, 4.0, -1.0, 0.0, -1.0, -1.0, 0.0, -1.0,&
4.0, -1.0, 0.0, -1.0, -1.0, 0.0, -1.0, 4.0, -1.0, 0.0,&
-1.0, -1.0, 0.0, -1.0, 4.0, -1.0, 0.0, -1.0, -1.0, 0.0,&
-1.0, 4.0, -1.0, 2*0.0, -1.0, 0.0, -1.0, 4.0, -1.0, 2*0.0,&
-1.0, 0.0, -1.0, 4.0, 3*0.0/
! Set PRECND in band symmetric form
DATA PRECND/0.0, 4.0, -1.0, 4.0, -1.0, 4.0, -1.0, 4.0, -1.0, 4.0,&
-1.0, 4.0, -1.0, 4.0, -1.0, 4.0, -1.0, 4.0/
! Right side is (1, ..., 1)
R = 1.0E0
! Initial guess for X is 0
X = 0.0E0
! Factor the preconditioning matrix
CALL LFCQS (PRECND, NCOPRE, PREFAC, RCOND)
!
ITMAX = 100
RELERR = 1.0E-4
IDO = 0
10 CALL PCGRC (IDO, X, P, R, Z, RELERR=RELERR, ITMAX=ITMAX)
IF (IDO .EQ. 1) THEN
! Set z = Ap
CALL MURBV (A, NCODA, NCODA, P, Z)
GO TO 10
ELSE IF (IDO .EQ. 2) THEN
! Solve PRECND*z = r for r
CALL LSLQS (PREFAC, NCOPRE, R, Z)
GO TO 10
END IF
! Print the solution
CALL WRRRN (’Solution’, X)
!
END
Solution
1 0.955
2 1.241
3 1.349
4 1.578
5 1.660
6 1.578
7 1.349
8 1.241
9 0.955
