2 SUBROUTINE DPCGS
(NEQ
, TN
, Y
, SAVF
, R
, WGHT
, N
, MAXL
, DELTA
, HL0
,
3 1 JPRE
, MNEWT
, F
, PSOL
, NPSL
, X
, P
, W
, Z
, LPCG
, WP
, IWP
, WK
, IFLAG
)
5 INTEGER NEQ
, N
, MAXL
, JPRE
, MNEWT
, NPSL
, LPCG
, IWP
, IFLAG
6 DOUBLE PRECISION TN
,Y
,SAVF
,R
,WGHT
,DELTA
,HL0
,X
,P
,W
,Z
,WP
,WK
7 DIMENSION NEQ
(*), Y
(*), SAVF
(*), R
(*), WGHT
(*), X
(*), P
(*), W
(*),
8 1 Z
(*), WP
(*), IWP
(*), WK
(*)
9 C-----------------------------------------------------------------------
10 C This routine computes the solution to the system A*x = b using a
11 C scaled preconditioned version of the Conjugate Gradient algorithm.
12 C It is assumed here that the scaled matrix D**-1 * A * D and the
13 C scaled preconditioner D**-1 * M * D are close to being
14 C symmetric positive definite.
15 C-----------------------------------------------------------------------
19 C NEQ = problem size, passed to F and PSOL (NEQ(1) = N).
21 C TN = current value of t.
23 C Y = array containing current dependent variable vector.
25 C SAVF = array containing current value of f(t,y).
27 C R = the right hand side of the system A*x = b.
29 C WGHT = array of length N containing scale factors.
30 C 1/WGHT(i) are the diagonal elements of the diagonal
33 C N = the order of the matrix A, and the lengths
34 C of the vectors Y, SAVF, R, WGHT, P, W, Z, WK, and X.
36 C MAXL = the maximum allowable number of iterates.
38 C DELTA = tolerance on residuals b - A*x in weighted RMS-norm.
40 C HL0 = current value of (step size h) * (coefficient l0).
42 C JPRE = preconditioner type flag.
44 C MNEWT = Newton iteration counter (.ge. 0).
46 C WK = real work array used by routine DATP.
48 C WP = real work array used by preconditioner PSOL.
50 C IWP = integer work array used by preconditioner PSOL.
54 C X = the final computed approximation to the solution
55 C of the system A*x = b.
57 C LPCG = the number of iterations performed, and current
58 C order of the upper Hessenberg matrix HES.
60 C NPSL = the number of calls to PSOL.
62 C IFLAG = integer error flag:
63 C 0 means convergence in LPCG iterations, LPCG .le. MAXL.
64 C 1 means the convergence test did not pass in MAXL
65 C iterations, but the residual norm is .lt. 1,
66 C or .lt. norm(b) if MNEWT = 0, and so X is computed.
67 C 2 means the convergence test did not pass in MAXL
68 C iterations, residual .gt. 1, and X is undefined.
69 C 3 means there was a recoverable error in PSOL
70 C caused by the preconditioner being out of date.
71 C 4 means there was a zero denominator in the algorithm.
72 C the scaled matrix or scaled preconditioner is not
73 C sufficiently close to being symmetric pos. definite.
74 C -1 means there was a nonrecoverable error in PSOL.
76 C-----------------------------------------------------------------------
78 DOUBLE PRECISION ALPHA
, BETA
, BNRM
, PTW
, RNRM
, DVNORM
, ZTR
, ZTR0
85 BNRM
= DVNORM
(N
, R
, WGHT
)
86 C Test for immediate return with X = 0 or X = b. -----------------------
87 IF (BNRM
.GT
. DELTA
) GO TO 20
88 IF (MNEWT
.GT
. 0) RETURN
89 CALL DCOPY
(N
, R
, 1, X
, 1)
93 C Loop point for PCG iterations. ---------------------------------------
96 CALL DCOPY
(N
, R
, 1, Z
, 1)
98 IF (JPRE
.EQ
. 0) GO TO 40
99 CALL PSOL
(NEQ
, TN
, Y
, SAVF
, WK
, HL0
, WP
, IWP
, Z
, 3, IER
)
101 IF (IER
.NE
. 0) GO TO 100
106 45 ZTR
= ZTR
+ Z
(I
)*R
(I
)*WGHT
(I
)**2
107 IF (LPCG
.NE
. 1) GO TO 50
108 CALL DCOPY
(N
, Z
, 1, P
, 1)
111 IF (ZTR0
.EQ
. 0.0D0
) GO TO 200
114 60 P
(I
) = Z
(I
) + BETA*P
(I
)
116 C-----------------------------------------------------------------------
117 C Call DATP to compute A*p and return the answer in W.
118 C-----------------------------------------------------------------------
119 CALL DATP
(NEQ
, Y
, SAVF
, P
, WGHT
, HL0
, WK
, F
, W
)
123 80 PTW
= PTW
+ P
(I
)*W
(I
)*WGHT
(I
)**2
124 IF (PTW
.EQ
. 0.0D0
) GO TO 200
126 CALL DAXPY
(N
, ALPHA
, P
, 1, X
, 1)
128 CALL DAXPY
(N
, ALPHA
, W
, 1, R
, 1)
129 RNRM
= DVNORM
(N
, R
, WGHT
)
130 IF (RNRM
.LE
. DELTA
) RETURN
131 IF (LPCG
.LT
. MAXL
) GO TO 30
133 IF (RNRM
.LE
. 1.0D0
) IFLAG
= 1
134 IF (RNRM
.LE
. BNRM
.AND
. MNEWT
.EQ
. 0) IFLAG
= 1
136 C-----------------------------------------------------------------------
137 C This block handles error returns from PSOL.
138 C-----------------------------------------------------------------------
140 IF (IER
.LT
. 0) IFLAG
= -1
141 IF (IER
.GT
. 0) IFLAG
= 3
143 C-----------------------------------------------------------------------
144 C This block handles division by zero errors.
145 C-----------------------------------------------------------------------
149 C----------------------- End of Subroutine DPCGS -----------------------