1*DECK PJAC 2 SUBROUTINE PJAC (NEQ, Y, YH, NYH, EWT, FTEM, SAVF, WM, IWM, F, 3 + JAC, RPAR, IPAR) 4C***BEGIN PROLOGUE PJAC 5C***SUBSIDIARY 6C***PURPOSE Subsidiary to DEBDF 7C***LIBRARY SLATEC 8C***TYPE SINGLE PRECISION (PJAC-S, DPJAC-D) 9C***AUTHOR Watts, H. A., (SNLA) 10C***DESCRIPTION 11C 12C PJAC sets up the iteration matrix (involving the Jacobian) for the 13C integration package DEBDF. 14C 15C***SEE ALSO DEBDF 16C***ROUTINES CALLED SGBFA, SGEFA, VNWRMS 17C***COMMON BLOCKS DEBDF1 18C***REVISION HISTORY (YYMMDD) 19C 800901 DATE WRITTEN 20C 890531 Changed all specific intrinsics to generic. (WRB) 21C 891214 Prologue converted to Version 4.0 format. (BAB) 22C 900328 Added TYPE section. (WRB) 23C 910722 Updated AUTHOR section. (ALS) 24C 920422 Changed DIMENSION statement. (WRB) 25C***END PROLOGUE PJAC 26C 27CLLL. OPTIMIZE 28 INTEGER NEQ, NYH, IWM, I, I1, I2, IER, II, IOWND, IOWNS, J, J1, 29 1 JJ, JSTART, KFLAG, L, LENP, MAXORD, MBA, MBAND, MEB1, MEBAND, 30 2 METH, MITER, ML, ML3, MU, N, NFE, NJE, NQ, NQU, NST 31 EXTERNAL F, JAC 32 REAL Y, YH, EWT, FTEM, SAVF, WM, 33 1 ROWND, ROWNS, EL0, H, HMIN, HMXI, HU, TN, UROUND, 34 2 CON, DI, FAC, HL0, R, R0, SRUR, YI, YJ, YJJ, VNWRMS 35 DIMENSION Y(*), YH(NYH,*), EWT(*), FTEM(*), SAVF(*), 36 1 WM(*), IWM(*), RPAR(*), IPAR(*) 37 COMMON /DEBDF1/ ROWND, ROWNS(210), 38 1 EL0, H, HMIN, HMXI, HU, TN, UROUND, IOWND(14), IOWNS(6), 39 2 IER, JSTART, KFLAG, L, METH, MITER, MAXORD, N, NQ, NST, NFE, 40 3 NJE, NQU 41C----------------------------------------------------------------------- 42C PJAC IS CALLED BY STOD TO COMPUTE AND PROCESS THE MATRIX 43C P = I - H*EL(1)*J , WHERE J IS AN APPROXIMATION TO THE JACOBIAN. 44C HERE J IS COMPUTED BY THE USER-SUPPLIED ROUTINE JAC IF 45C MITER = 1 OR 4, OR BY FINITE DIFFERENCING IF MITER = 2, 3, OR 5. 46C IF MITER = 3, A DIAGONAL APPROXIMATION TO J IS USED. 47C J IS STORED IN WM AND REPLACED BY P. IF MITER .NE. 3, P IS THEN 48C SUBJECTED TO LU DECOMPOSITION IN PREPARATION FOR LATER SOLUTION 49C OF LINEAR SYSTEMS WITH P AS COEFFICIENT MATRIX. THIS IS DONE 50C BY SGEFA IF MITER = 1 OR 2, AND BY SGBFA IF MITER = 4 OR 5. 51C 52C IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION 53C WITH PJAC USES THE FOLLOWING.. 54C Y = ARRAY CONTAINING PREDICTED VALUES ON ENTRY. 55C FTEM = WORK ARRAY OF LENGTH N (ACOR IN STOD ). 56C SAVF = ARRAY CONTAINING F EVALUATED AT PREDICTED Y. 57C WM = REAL WORK SPACE FOR MATRICES. ON OUTPUT IT CONTAINS THE 58C INVERSE DIAGONAL MATRIX IF MITER = 3 AND THE LU DECOMPOSITION 59C OF P IF MITER IS 1, 2 , 4, OR 5. 60C STORAGE OF MATRIX ELEMENTS STARTS AT WM(3). 61C WM ALSO CONTAINS THE FOLLOWING MATRIX-RELATED DATA.. 62C WM(1) = SQRT(UROUND), USED IN NUMERICAL JACOBIAN INCREMENTS. 63C WM(2) = H*EL0, SAVED FOR LATER USE IF MITER = 3. 64C IWM = INTEGER WORK SPACE CONTAINING PIVOT INFORMATION, STARTING AT 65C IWM(21), IF MITER IS 1, 2, 4, OR 5. IWM ALSO CONTAINS THE 66C BAND PARAMETERS ML = IWM(1) AND MU = IWM(2) IF MITER IS 4 OR 5. 67C EL0 = EL(1) (INPUT). 68C IER = OUTPUT ERROR FLAG, = 0 IF NO TROUBLE, .NE. 0 IF 69C P MATRIX FOUND TO BE SINGULAR. 70C THIS ROUTINE ALSO USES THE COMMON VARIABLES EL0, H, TN, UROUND, 71C MITER, N, NFE, AND NJE. 72C----------------------------------------------------------------------- 73C***FIRST EXECUTABLE STATEMENT PJAC 74 NJE = NJE + 1 75 HL0 = H*EL0 76 GO TO (100, 200, 300, 400, 500), MITER 77C IF MITER = 1, CALL JAC AND MULTIPLY BY SCALAR. ----------------------- 78 100 LENP = N*N 79 DO 110 I = 1,LENP 80 110 WM(I+2) = 0.0E0 81 CALL JAC (TN, Y, WM(3), N, RPAR, IPAR) 82 CON = -HL0 83 DO 120 I = 1,LENP 84 120 WM(I+2) = WM(I+2)*CON 85 GO TO 240 86C IF MITER = 2, MAKE N CALLS TO F TO APPROXIMATE J. -------------------- 87 200 FAC = VNWRMS (N, SAVF, EWT) 88 R0 = 1000.0E0*ABS(H)*UROUND*N*FAC 89 IF (R0 .EQ. 0.0E0) R0 = 1.0E0 90 SRUR = WM(1) 91 J1 = 2 92 DO 230 J = 1,N 93 YJ = Y(J) 94 R = MAX(SRUR*ABS(YJ),R0*EWT(J)) 95 Y(J) = Y(J) + R 96 FAC = -HL0/R 97 CALL F (TN, Y, FTEM, RPAR, IPAR) 98 DO 220 I = 1,N 99 220 WM(I+J1) = (FTEM(I) - SAVF(I))*FAC 100 Y(J) = YJ 101 J1 = J1 + N 102 230 CONTINUE 103 NFE = NFE + N 104C ADD IDENTITY MATRIX. ------------------------------------------------- 105 240 J = 3 106 DO 250 I = 1,N 107 WM(J) = WM(J) + 1.0E0 108 250 J = J + (N + 1) 109C DO LU DECOMPOSITION ON P. -------------------------------------------- 110 CALL SGEFA (WM(3), N, N, IWM(21), IER) 111 RETURN 112C IF MITER = 3, CONSTRUCT A DIAGONAL APPROXIMATION TO J AND P. --------- 113 300 WM(2) = HL0 114 IER = 0 115 R = EL0*0.1E0 116 DO 310 I = 1,N 117 310 Y(I) = Y(I) + R*(H*SAVF(I) - YH(I,2)) 118 CALL F (TN, Y, WM(3), RPAR, IPAR) 119 NFE = NFE + 1 120 DO 320 I = 1,N 121 R0 = H*SAVF(I) - YH(I,2) 122 DI = 0.1E0*R0 - H*(WM(I+2) - SAVF(I)) 123 WM(I+2) = 1.0E0 124 IF (ABS(R0) .LT. UROUND*EWT(I)) GO TO 320 125 IF (ABS(DI) .EQ. 0.0E0) GO TO 330 126 WM(I+2) = 0.1E0*R0/DI 127 320 CONTINUE 128 RETURN 129 330 IER = -1 130 RETURN 131C IF MITER = 4, CALL JAC AND MULTIPLY BY SCALAR. ----------------------- 132 400 ML = IWM(1) 133 MU = IWM(2) 134 ML3 = 3 135 MBAND = ML + MU + 1 136 MEBAND = MBAND + ML 137 LENP = MEBAND*N 138 DO 410 I = 1,LENP 139 410 WM(I+2) = 0.0E0 140 CALL JAC (TN, Y, WM(ML3), MEBAND, RPAR, IPAR) 141 CON = -HL0 142 DO 420 I = 1,LENP 143 420 WM(I+2) = WM(I+2)*CON 144 GO TO 570 145C IF MITER = 5, MAKE MBAND CALLS TO F TO APPROXIMATE J. ---------------- 146 500 ML = IWM(1) 147 MU = IWM(2) 148 MBAND = ML + MU + 1 149 MBA = MIN(MBAND,N) 150 MEBAND = MBAND + ML 151 MEB1 = MEBAND - 1 152 SRUR = WM(1) 153 FAC = VNWRMS (N, SAVF, EWT) 154 R0 = 1000.0E0*ABS(H)*UROUND*N*FAC 155 IF (R0 .EQ. 0.0E0) R0 = 1.0E0 156 DO 560 J = 1,MBA 157 DO 530 I = J,N,MBAND 158 YI = Y(I) 159 R = MAX(SRUR*ABS(YI),R0*EWT(I)) 160 530 Y(I) = Y(I) + R 161 CALL F (TN, Y, FTEM, RPAR, IPAR) 162 DO 550 JJ = J,N,MBAND 163 Y(JJ) = YH(JJ,1) 164 YJJ = Y(JJ) 165 R = MAX(SRUR*ABS(YJJ),R0*EWT(JJ)) 166 FAC = -HL0/R 167 I1 = MAX(JJ-MU,1) 168 I2 = MIN(JJ+ML,N) 169 II = JJ*MEB1 - ML + 2 170 DO 540 I = I1,I2 171 540 WM(II+I) = (FTEM(I) - SAVF(I))*FAC 172 550 CONTINUE 173 560 CONTINUE 174 NFE = NFE + MBA 175C ADD IDENTITY MATRIX. ------------------------------------------------- 176 570 II = MBAND + 2 177 DO 580 I = 1,N 178 WM(II) = WM(II) + 1.0E0 179 580 II = II + MEBAND 180C DO LU DECOMPOSITION OF P. -------------------------------------------- 181 CALL SGBFA (WM(3), MEBAND, N, ML, MU, IWM(21), IER) 182 RETURN 183C----------------------- END OF SUBROUTINE PJAC ----------------------- 184 END 185