*DECK DPCGS
SUBROUTINE DPCGS (NEQ, TN, Y, SAVF, R, WGHT, N, MAXL, DELTA, HL0,
1 JPRE, MNEWT, F, PSOL, NPSL, X, P, W, Z, LPCG, WP, IWP, WK, IFLAG)
EXTERNAL F, PSOL
INTEGER NEQ, N, MAXL, JPRE, MNEWT, NPSL, LPCG, IWP, IFLAG
DOUBLE PRECISION TN,Y,SAVF,R,WGHT,DELTA,HL0,X,P,W,Z,WP,WK
DIMENSION NEQ(*), Y(*), SAVF(*), R(*), WGHT(*), X(*), P(*), W(*),
1 Z(*), WP(*), IWP(*), WK(*)
C-----------------------------------------------------------------------
C This routine computes the solution to the system A*x = b using a
C scaled preconditioned version of the Conjugate Gradient algorithm.
C It is assumed here that the scaled matrix D**-1 * A * D and the
C scaled preconditioner D**-1 * M * D are close to being
C symmetric positive definite.
C-----------------------------------------------------------------------
C
C On entry
C
C NEQ = problem size, passed to F and PSOL (NEQ(1) = N).
C
C TN = current value of t.
C
C Y = array containing current dependent variable vector.
C
C SAVF = array containing current value of f(t,y).
C
C R = the right hand side of the system A*x = b.
C
C WGHT = array of length N containing scale factors.
C 1/WGHT(i) are the diagonal elements of the diagonal
C scaling matrix D.
C
C N = the order of the matrix A, and the lengths
C of the vectors Y, SAVF, R, WGHT, P, W, Z, WK, and X.
C
C MAXL = the maximum allowable number of iterates.
C
C DELTA = tolerance on residuals b - A*x in weighted RMS-norm.
C
C HL0 = current value of (step size h) * (coefficient l0).
C
C JPRE = preconditioner type flag.
C
C MNEWT = Newton iteration counter (.ge. 0).
C
C WK = real work array used by routine DATP.
C
C WP = real work array used by preconditioner PSOL.
C
C IWP = integer work array used by preconditioner PSOL.
C
C On return
C
C X = the final computed approximation to the solution
C of the system A*x = b.
C
C LPCG = the number of iterations performed, and current
C order of the upper Hessenberg matrix HES.
C
C NPSL = the number of calls to PSOL.
C
C IFLAG = integer error flag:
C 0 means convergence in LPCG iterations, LPCG .le. MAXL.
C 1 means the convergence test did not pass in MAXL
C iterations, but the residual norm is .lt. 1,
C or .lt. norm(b) if MNEWT = 0, and so X is computed.
C 2 means the convergence test did not pass in MAXL
C iterations, residual .gt. 1, and X is undefined.
C 3 means there was a recoverable error in PSOL
C caused by the preconditioner being out of date.
C 4 means there was a zero denominator in the algorithm.
C the scaled matrix or scaled preconditioner is not
C sufficiently close to being symmetric pos. definite.
C -1 means there was a nonrecoverable error in PSOL.
C
C-----------------------------------------------------------------------
INTEGER I, IER
DOUBLE PRECISION ALPHA, BETA, BNRM, PTW, RNRM, DVNORM, ZTR, ZTR0
C
IFLAG = 0
NPSL = 0
LPCG = 0
DO 10 I = 1,N
10 X(I) = 0.0D0
BNRM = DVNORM (N, R, WGHT)
C Test for immediate return with X = 0 or X = b. -----------------------
IF (BNRM .GT. DELTA) GO TO 20
IF (MNEWT .GT. 0) RETURN
CALL DCOPY (N, R, 1, X, 1)
RETURN
C
20 ZTR = 0.0D0
C Loop point for PCG iterations. ---------------------------------------
30 CONTINUE
LPCG = LPCG + 1
CALL DCOPY (N, R, 1, Z, 1)
IER = 0
IF (JPRE .EQ. 0) GO TO 40
CALL PSOL (NEQ, TN, Y, SAVF, WK, HL0, WP, IWP, Z, 3, IER)
NPSL = NPSL + 1
IF (IER .NE. 0) GO TO 100
40 CONTINUE
ZTR0 = ZTR
ZTR = 0.0D0
DO 45 I = 1,N
45 ZTR = ZTR + Z(I)*R(I)*WGHT(I)**2
IF (LPCG .NE. 1) GO TO 50
CALL DCOPY (N, Z, 1, P, 1)
GO TO 70
50 CONTINUE
IF (ZTR0 .EQ. 0.0D0) GO TO 200
BETA = ZTR/ZTR0
DO 60 I = 1,N
60 P(I) = Z(I) + BETA*P(I)
70 CONTINUE
C-----------------------------------------------------------------------
C Call DATP to compute A*p and return the answer in W.
C-----------------------------------------------------------------------
CALL DATP (NEQ, Y, SAVF, P, WGHT, HL0, WK, F, W)
C
PTW = 0.0D0
DO 80 I = 1,N
80 PTW = PTW + P(I)*W(I)*WGHT(I)**2
IF (PTW .EQ. 0.0D0) GO TO 200
ALPHA = ZTR/PTW
CALL DAXPY (N, ALPHA, P, 1, X, 1)
ALPHA = -ALPHA
CALL DAXPY (N, ALPHA, W, 1, R, 1)
RNRM = DVNORM (N, R, WGHT)
IF (RNRM .LE. DELTA) RETURN
IF (LPCG .LT. MAXL) GO TO 30
IFLAG = 2
IF (RNRM .LE. 1.0D0) IFLAG = 1
IF (RNRM .LE. BNRM .AND. MNEWT .EQ. 0) IFLAG = 1
RETURN
C-----------------------------------------------------------------------
C This block handles error returns from PSOL.
C-----------------------------------------------------------------------
100 CONTINUE
IF (IER .LT. 0) IFLAG = -1
IF (IER .GT. 0) IFLAG = 3
RETURN
C-----------------------------------------------------------------------
C This block handles division by zero errors.
C-----------------------------------------------------------------------
200 CONTINUE
IFLAG = 4
RETURN
C----------------------- End of Subroutine DPCGS -----------------------
END