Forgot to load lapack in a few examples
[maxima.git] / share / odepack / fortran / dpcgs.f
blob64e86c96b06c4f96b170d618945fd9d0efad4405
1 *DECK DPCGS
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)
4 EXTERNAL F, PSOL
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-----------------------------------------------------------------------
17 C On entry
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
31 C scaling matrix D.
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.
52 C On return
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-----------------------------------------------------------------------
77 INTEGER I, IER
78 DOUBLE PRECISION ALPHA, BETA, BNRM, PTW, RNRM, DVNORM, ZTR, ZTR0
80 IFLAG = 0
81 NPSL = 0
82 LPCG = 0
83 DO 10 I = 1,N
84 10 X(I) = 0.0D0
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)
90 RETURN
92 20 ZTR = 0.0D0
93 C Loop point for PCG iterations. ---------------------------------------
94 30 CONTINUE
95 LPCG = LPCG + 1
96 CALL DCOPY (N, R, 1, Z, 1)
97 IER = 0
98 IF (JPRE .EQ. 0) GO TO 40
99 CALL PSOL (NEQ, TN, Y, SAVF, WK, HL0, WP, IWP, Z, 3, IER)
100 NPSL = NPSL + 1
101 IF (IER .NE. 0) GO TO 100
102 40 CONTINUE
103 ZTR0 = ZTR
104 ZTR = 0.0D0
105 DO 45 I = 1,N
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)
109 GO TO 70
110 50 CONTINUE
111 IF (ZTR0 .EQ. 0.0D0) GO TO 200
112 BETA = ZTR/ZTR0
113 DO 60 I = 1,N
114 60 P(I) = Z(I) + BETA*P(I)
115 70 CONTINUE
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)
121 PTW = 0.0D0
122 DO 80 I = 1,N
123 80 PTW = PTW + P(I)*W(I)*WGHT(I)**2
124 IF (PTW .EQ. 0.0D0) GO TO 200
125 ALPHA = ZTR/PTW
126 CALL DAXPY (N, ALPHA, P, 1, X, 1)
127 ALPHA = -ALPHA
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
132 IFLAG = 2
133 IF (RNRM .LE. 1.0D0) IFLAG = 1
134 IF (RNRM .LE. BNRM .AND. MNEWT .EQ. 0) IFLAG = 1
135 RETURN
136 C-----------------------------------------------------------------------
137 C This block handles error returns from PSOL.
138 C-----------------------------------------------------------------------
139 100 CONTINUE
140 IF (IER .LT. 0) IFLAG = -1
141 IF (IER .GT. 0) IFLAG = 3
142 RETURN
143 C-----------------------------------------------------------------------
144 C This block handles division by zero errors.
145 C-----------------------------------------------------------------------
146 200 CONTINUE
147 IFLAG = 4
148 RETURN
149 C----------------------- End of Subroutine DPCGS -----------------------