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.

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.

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 = 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.

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).

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 = N.

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

Example 1

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

Output

Solution

1 1.001

2 -4.000

3 7.000

Example 2

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

Output

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