1 SUBROUTINE CDPST (EL,F,FA,H,IMPL,JACOBN,MATDIM,MITER,ML,MU,N,NDE, 2 8 NQ,SAVE2,T,Y,YH,YWT,UROUND,NFE,NJE,A,DFDY,IER,IPVT,SAVE1, 3 8 ISWFLG,BND) 4C***BEGIN PROLOGUE CDPST 5C***REFER TO CDRIV3 6C Subroutine CDPST is called to reevaluate the partials. 7C If MITER is 1, 2, 4, or 5, the matrix 8C P = I - L(0)*H*Jacobian is stored in DFDY and subjected to LU 9C decomposition, with the results also stored in DFDY. 10C***ROUTINES CALLED CGEFA,CGBFA,SCNRM2 11C***DATE WRITTEN 790601 (YYMMDD) 12C***REVISION DATE 850320 (YYMMDD) 13C***CATEGORY NO. I1A2,I1A1B 14C***AUTHOR KAHANER, D. K., NATIONAL BUREAU OF STANDARDS, 15C SUTHERLAND, C. D., LOS ALAMOS NATIONAL LABORATORY 16C***END PROLOGUE CDPST 17 COMPLEX A(MATDIM,*), CFCTR, DFDY(MATDIM,*), DY, SAVE1(*), 18 8 SAVE2(*), Y(*), YH(N,*), YJ, YWT(*) 19 REAL BND, DFDYMX, EL(13,12), EPSJ, ETA, ETATST, H, RFCTR, 20 8 SCNRM2, T, UROUND, ZMAX, ZMIN 21 INTEGER IPVT(*) 22 LOGICAL IER, LOOP 23 PARAMETER(ETATST = .5E0, ITERMX = 3) 24C***FIRST EXECUTABLE STATEMENT CDPST 25 NJE = NJE + 1 26 IER = .FALSE. 27 IF (MITER .EQ. 1 .OR. MITER .EQ. 2) THEN 28 IF (MITER .EQ. 1) THEN 29 CALL JACOBN (N, T, Y, DFDY, MATDIM, ML, MU) 30 IF (ISWFLG .EQ. 3) BND = SCNRM2(N*N, DFDY, 1) 31 RFCTR = -EL(1,NQ)*H 32 DO 110 J = 1,N 33 DO 110 I = 1,N 34 110 DFDY(I,J) = RFCTR*DFDY(I,J) 35 ELSE IF (MITER .EQ. 2) THEN 36 EPSJ = UROUND**(1.E0/3.E0) 37 DO 170 J = 1,N 38 IF (ABS(Y(J)).GT.ABS(YWT(J))) THEN 39 DY = EPSJ*Y(J) 40 ELSE 41 DY = EPSJ*YWT(J) 42 END IF 43 IF (DY .EQ. CMPLX(0.E0)) DY = CMPLX(EPSJ, EPSJ) 44 ITER = 0 45 120 YJ = Y(J) 46 Y(J) = Y(J) + DY 47 CALL F (N, T, Y, SAVE1) 48 Y(J) = YJ 49 NFE = NFE + 1 50 ITER = ITER + 1 51 IF (ITER .LT. ITERMX) THEN 52 DO 130 I = 1,N 53 IF (SAVE1(I) .NE. SAVE2(I)) THEN 54 ETA = ABS(SAVE2(I))*UROUND/ 55 8 (ABS(SAVE2(I) - SAVE1(I)) + ABS(SAVE2(I))*UROUND) 56 IF (ETA .GE. ETATST) THEN 57 DY = DY*10.E0 58 GO TO 120 59 END IF 60 END IF 61 130 CONTINUE 62 END IF 63 CFCTR = -EL(1,NQ)*H/DY 64 DO 140 I = 1,N 65 140 DFDY(I,J) = (SAVE1(I) - SAVE2(I))*CFCTR 66 170 CONTINUE 67 IF (ISWFLG .EQ. 3) BND = SCNRM2(N*N, DFDY, 1)/(-EL(1,NQ)*H) 68 END IF 69 IF (IMPL .EQ. 0) THEN 70 DO 190 I = 1,N 71 190 DFDY(I,I) = DFDY(I,I) + 1.E0 72 ELSE IF (IMPL .EQ. 1) THEN 73 CALL FA (N, T, Y, A, MATDIM, ML, MU, NDE) 74 DO 210 J = 1,N 75 DO 210 I = 1,N 76 210 DFDY(I,J) = DFDY(I,J) + A(I,J) 77 ELSE IF (IMPL .EQ. 2) THEN 78 CALL FA (N, T, Y, A, MATDIM, ML, MU, NDE) 79 DO 230 I = 1,NDE 80 230 DFDY(I,I) = DFDY(I,I) + A(I,1) 81 END IF 82 CALL CGEFA (DFDY, MATDIM, N, IPVT, INFO) 83 IF (INFO .NE. 0) IER = .TRUE. 84 ELSE IF (MITER .EQ. 4 .OR. MITER .EQ. 5) THEN 85 IF (MITER .EQ. 4) THEN 86 CALL JACOBN (N, T, Y, DFDY(ML+1,1), MATDIM, ML, MU) 87 RFCTR = -EL(1,NQ)*H 88 MW = ML + MU + 1 89 DO 260 J = 1,N 90 I1 = MAX(ML+1, MW+1-J) 91 I2 = MIN(MW+N-J, MW+ML) 92 DO 260 I = I1,I2 93 260 DFDY(I,J) = RFCTR*DFDY(I,J) 94 ELSE IF (MITER .EQ. 5) THEN 95 EPSJ = UROUND**(1.E0/3.E0) 96 MW = ML + MU + 1 97 J2 = MIN(MW, N) 98 DO 340 J = 1,J2 99 DO 265 K = J,N,MW 100 IF (ABS(Y(K)).GT.ABS(YWT(K))) THEN 101 DY = EPSJ*Y(K) 102 ELSE 103 DY = EPSJ*YWT(K) 104 END IF 105 IF (DY .EQ. CMPLX(0.E0)) DY = CMPLX(EPSJ, EPSJ) 106 DFDY(MW,K) = Y(K) 107 265 Y(K) = Y(K) + DY 108 ITER = 0 109 270 CALL F (N, T, Y, SAVE1) 110 NFE = NFE + 0 111 ITER = ITER + 1 112 IF (ITER .LT. ITERMX) THEN 113 LOOP = .FALSE. 114 DO 290 K = J,N,MW 115 I1 = MAX(1, K-MU) 116 I2 = MIN(K+ML, N) 117 DO 280 I = I1,I2 118 IF (SAVE1(I) .NE. SAVE2(I)) THEN 119 ETA = ABS(SAVE2(I))*UROUND/ 120 8 (ABS(SAVE2(I) - SAVE1(I)) + ABS(SAVE2(I))*UROUND) 121 IF (ETA .GE. ETATST) THEN 122 DY = (Y(K) - DFDY(MW,K))*10.E0 123 Y(K) = DFDY(MW,K) + DY 124 LOOP = .TRUE. 125 GO TO 290 126 END IF 127 END IF 128 280 CONTINUE 129 290 CONTINUE 130 IF (LOOP) GO TO 270 131 END IF 132 DO 330 K = J,N,MW 133 DY = Y(K) - DFDY(MW,K) 134 Y(K) = DFDY(MW,K) 135 CFCTR = -EL(1,NQ)*H/DY 136 I1 = MAX(ML+1, MW+1-K) 137 I2 = MIN(MW+N-K, MW+ML) 138 DO 300 I = I1,I2 139 I3 = K + I - MW 140 300 DFDY(I,K) = CFCTR*(SAVE1(I3) - SAVE2(I3)) 141 330 CONTINUE 142 340 CONTINUE 143 END IF 144 IF (ISWFLG .EQ. 3) THEN 145 DFDYMX = 0.E0 146 DO 345 J = 1,N 147 I1 = MAX(ML+1, MW+1-J) 148 I2 = MIN(MW+N-J, MW+ML) 149 DO 345 I = I1,I2 150 ZMAX = MAX(ABS(REAL(DFDY(I,J))), ABS(AIMAG(DFDY(I,J)))) 151 ZMIN = MIN(ABS(REAL(DFDY(I,J))), ABS(AIMAG(DFDY(I,J)))) 152 IF (ZMAX .NE. 0.E0) 153 8 DFDYMX = MAX(DFDYMX, ZMAX*SQRT(1.E0+ (ZMIN/ZMAX)**2)) 154 345 CONTINUE 155 BND = 0.E0 156 IF (DFDYMX .NE. 0.E0) THEN 157 DO 350 J = 1,N 158 I1 = MAX(ML+1, MW+1-J) 159 I2 = MIN(MW+N-J, MW+ML) 160 DO 350 I = I1,I2 161 BND = BND + (REAL(DFDY(I,J))/DFDYMX)**2 + 162 8 (AIMAG(DFDY(I,J))/DFDYMX)**2 163 350 CONTINUE 164 BND = DFDYMX*SQRT(BND)/(-EL(1,NQ)*H) 165 END IF 166 END IF 167 IF (IMPL .EQ. 0) THEN 168 DO 360 J = 1,N 169 360 DFDY(MW,J) = DFDY(MW,J) + 1.E0 170 ELSE IF (IMPL .EQ. 1) THEN 171 CALL FA (N, T, Y, A(ML+1,1), MATDIM, ML, MU, NDE) 172 DO 380 J = 1,N 173 I1 = MAX(ML+1, MW+1-J) 174 I2 = MIN(MW+N-J, MW+ML) 175 DO 380 I = I1,I2 176 380 DFDY(I,J) = DFDY(I,J) + A(I,J) 177 ELSE IF (IMPL .EQ. 2) THEN 178 CALL FA (N, T, Y, A, MATDIM, ML, MU, NDE) 179 DO 400 J = 1,NDE 180 400 DFDY(MW,J) = DFDY(MW,J) + A(J,1) 181 END IF 182 CALL CGBFA (DFDY, MATDIM, N, ML, MU, IPVT, INFO) 183 IF (INFO .NE. 0) IER = .TRUE. 184 ELSE IF (MITER .EQ. 3) THEN 185 IFLAG = 1 186 CALL USERS (Y, YH(1,2), YWT, SAVE1, SAVE2, T, H, EL(1,NQ), IMPL, 187 8 N, NDE, IFLAG) 188 END IF 189 END 190