Mercurial > octave
comparison libcruft/odepack/prepj.f @ 1644:395bb6d6c096
[project @ 1995-12-14 08:32:12 by jwe]
Initial revision
| author | jwe |
|---|---|
| date | Thu, 14 Dec 1995 08:32:21 +0000 |
| parents | |
| children | 44ed237bdc1e |
comparison
equal
deleted
inserted
replaced
| 1643:5e108d51e370 | 1644:395bb6d6c096 |
|---|---|
| 1 SUBROUTINE PREPJ (NEQ, Y, YH, NYH, EWT, FTEM, SAVF, WM, IWM, | |
| 2 1 F, JAC, IERR) | |
| 3 CLLL. OPTIMIZE | |
| 4 EXTERNAL F, JAC | |
| 5 INTEGER NEQ, NYH, IWM | |
| 6 INTEGER IOWND, IOWNS, | |
| 7 1 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, MITER, | |
| 8 2 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU | |
| 9 INTEGER I, I1, I2, IER, II, J, J1, JJ, LENP, | |
| 10 1 MBA, MBAND, MEB1, MEBAND, ML, ML3, MU, NP1 | |
| 11 DOUBLE PRECISION Y, YH, EWT, FTEM, SAVF, WM | |
| 12 DOUBLE PRECISION ROWNS, | |
| 13 1 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND | |
| 14 DOUBLE PRECISION CON, DI, FAC, HL0, R, R0, SRUR, YI, YJ, YJJ, | |
| 15 1 VNORM | |
| 16 DIMENSION NEQ(1), Y(1), YH(NYH,1), EWT(1), FTEM(1), SAVF(1), | |
| 17 1 WM(1), IWM(1) | |
| 18 COMMON /LS0001/ ROWNS(209), | |
| 19 2 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, | |
| 20 3 IOWND(14), IOWNS(6), | |
| 21 4 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, MITER, | |
| 22 5 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU | |
| 23 C----------------------------------------------------------------------- | |
| 24 C PREPJ IS CALLED BY STODE TO COMPUTE AND PROCESS THE MATRIX | |
| 25 C P = I - H*EL(1)*J , WHERE J IS AN APPROXIMATION TO THE JACOBIAN. | |
| 26 C HERE J IS COMPUTED BY THE USER-SUPPLIED ROUTINE JAC IF | |
| 27 C MITER = 1 OR 4, OR BY FINITE DIFFERENCING IF MITER = 2, 3, OR 5. | |
| 28 C IF MITER = 3, A DIAGONAL APPROXIMATION TO J IS USED. | |
| 29 C J IS STORED IN WM AND REPLACED BY P. IF MITER .NE. 3, P IS THEN | |
| 30 C SUBJECTED TO LU DECOMPOSITION IN PREPARATION FOR LATER SOLUTION | |
| 31 C OF LINEAR SYSTEMS WITH P AS COEFFICIENT MATRIX. THIS IS DONE | |
| 32 C BY DGEFA IF MITER = 1 OR 2, AND BY DGBFA IF MITER = 4 OR 5. | |
| 33 C | |
| 34 C IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION | |
| 35 C WITH PREPJ USES THE FOLLOWING.. | |
| 36 C Y = ARRAY CONTAINING PREDICTED VALUES ON ENTRY. | |
| 37 C FTEM = WORK ARRAY OF LENGTH N (ACOR IN STODE). | |
| 38 C SAVF = ARRAY CONTAINING F EVALUATED AT PREDICTED Y. | |
| 39 C WM = REAL WORK SPACE FOR MATRICES. ON OUTPUT IT CONTAINS THE | |
| 40 C INVERSE DIAGONAL MATRIX IF MITER = 3 AND THE LU DECOMPOSITION | |
| 41 C OF P IF MITER IS 1, 2 , 4, OR 5. | |
| 42 C STORAGE OF MATRIX ELEMENTS STARTS AT WM(3). | |
| 43 C WM ALSO CONTAINS THE FOLLOWING MATRIX-RELATED DATA.. | |
| 44 C WM(1) = SQRT(UROUND), USED IN NUMERICAL JACOBIAN INCREMENTS. | |
| 45 C WM(2) = H*EL0, SAVED FOR LATER USE IF MITER = 3. | |
| 46 C IWM = INTEGER WORK SPACE CONTAINING PIVOT INFORMATION, STARTING AT | |
| 47 C IWM(21), IF MITER IS 1, 2, 4, OR 5. IWM ALSO CONTAINS BAND | |
| 48 C PARAMETERS ML = IWM(1) AND MU = IWM(2) IF MITER IS 4 OR 5. | |
| 49 C EL0 = EL(1) (INPUT). | |
| 50 C IERPJ = OUTPUT ERROR FLAG, = 0 IF NO TROUBLE, .GT. 0 IF | |
| 51 C P MATRIX FOUND TO BE SINGULAR. | |
| 52 C JCUR = OUTPUT FLAG = 1 TO INDICATE THAT THE JACOBIAN MATRIX | |
| 53 C (OR APPROXIMATION) IS NOW CURRENT. | |
| 54 C THIS ROUTINE ALSO USES THE COMMON VARIABLES EL0, H, TN, UROUND, | |
| 55 C MITER, N, NFE, AND NJE. | |
| 56 C----------------------------------------------------------------------- | |
| 57 NJE = NJE + 1 | |
| 58 IERPJ = 0 | |
| 59 JCUR = 1 | |
| 60 HL0 = H*EL0 | |
| 61 GO TO (100, 200, 300, 400, 500), MITER | |
| 62 C IF MITER = 1, CALL JAC AND MULTIPLY BY SCALAR. ----------------------- | |
| 63 100 LENP = N*N | |
| 64 DO 110 I = 1,LENP | |
| 65 110 WM(I+2) = 0.0D0 | |
| 66 CALL JAC (NEQ, TN, Y, 0, 0, WM(3), N) | |
| 67 CON = -HL0 | |
| 68 DO 120 I = 1,LENP | |
| 69 120 WM(I+2) = WM(I+2)*CON | |
| 70 GO TO 240 | |
| 71 C IF MITER = 2, MAKE N CALLS TO F TO APPROXIMATE J. -------------------- | |
| 72 200 FAC = VNORM (N, SAVF, EWT) | |
| 73 R0 = 1000.0D0*DABS(H)*UROUND*DBLE(N)*FAC | |
| 74 IF (R0 .EQ. 0.0D0) R0 = 1.0D0 | |
| 75 SRUR = WM(1) | |
| 76 J1 = 2 | |
| 77 DO 230 J = 1,N | |
| 78 YJ = Y(J) | |
| 79 R = DMAX1(SRUR*DABS(YJ),R0/EWT(J)) | |
| 80 Y(J) = Y(J) + R | |
| 81 FAC = -HL0/R | |
| 82 IERR = 0 | |
| 83 CALL F (NEQ, TN, Y, FTEM, IERR) | |
| 84 IF (IERR .LT. 0) RETURN | |
| 85 DO 220 I = 1,N | |
| 86 220 WM(I+J1) = (FTEM(I) - SAVF(I))*FAC | |
| 87 Y(J) = YJ | |
| 88 J1 = J1 + N | |
| 89 230 CONTINUE | |
| 90 NFE = NFE + N | |
| 91 C ADD IDENTITY MATRIX. ------------------------------------------------- | |
| 92 240 J = 3 | |
| 93 NP1 = N + 1 | |
| 94 DO 250 I = 1,N | |
| 95 WM(J) = WM(J) + 1.0D0 | |
| 96 250 J = J + NP1 | |
| 97 C DO LU DECOMPOSITION ON P. -------------------------------------------- | |
| 98 CALL DGEFA (WM(3), N, N, IWM(21), IER) | |
| 99 IF (IER .NE. 0) IERPJ = 1 | |
| 100 RETURN | |
| 101 C IF MITER = 3, CONSTRUCT A DIAGONAL APPROXIMATION TO J AND P. --------- | |
| 102 300 WM(2) = HL0 | |
| 103 R = EL0*0.1D0 | |
| 104 DO 310 I = 1,N | |
| 105 310 Y(I) = Y(I) + R*(H*SAVF(I) - YH(I,2)) | |
| 106 IERR = 0 | |
| 107 CALL F (NEQ, TN, Y, WM(3), IERR) | |
| 108 IF (IERR .LT. 0) RETURN | |
| 109 NFE = NFE + 1 | |
| 110 DO 320 I = 1,N | |
| 111 R0 = H*SAVF(I) - YH(I,2) | |
| 112 DI = 0.1D0*R0 - H*(WM(I+2) - SAVF(I)) | |
| 113 WM(I+2) = 1.0D0 | |
| 114 IF (DABS(R0) .LT. UROUND/EWT(I)) GO TO 320 | |
| 115 IF (DABS(DI) .EQ. 0.0D0) GO TO 330 | |
| 116 WM(I+2) = 0.1D0*R0/DI | |
| 117 320 CONTINUE | |
| 118 RETURN | |
| 119 330 IERPJ = 1 | |
| 120 RETURN | |
| 121 C IF MITER = 4, CALL JAC AND MULTIPLY BY SCALAR. ----------------------- | |
| 122 400 ML = IWM(1) | |
| 123 MU = IWM(2) | |
| 124 ML3 = ML + 3 | |
| 125 MBAND = ML + MU + 1 | |
| 126 MEBAND = MBAND + ML | |
| 127 LENP = MEBAND*N | |
| 128 DO 410 I = 1,LENP | |
| 129 410 WM(I+2) = 0.0D0 | |
| 130 CALL JAC (NEQ, TN, Y, ML, MU, WM(ML3), MEBAND) | |
| 131 CON = -HL0 | |
| 132 DO 420 I = 1,LENP | |
| 133 420 WM(I+2) = WM(I+2)*CON | |
| 134 GO TO 570 | |
| 135 C IF MITER = 5, MAKE MBAND CALLS TO F TO APPROXIMATE J. ---------------- | |
| 136 500 ML = IWM(1) | |
| 137 MU = IWM(2) | |
| 138 MBAND = ML + MU + 1 | |
| 139 MBA = MIN0(MBAND,N) | |
| 140 MEBAND = MBAND + ML | |
| 141 MEB1 = MEBAND - 1 | |
| 142 SRUR = WM(1) | |
| 143 FAC = VNORM (N, SAVF, EWT) | |
| 144 R0 = 1000.0D0*DABS(H)*UROUND*DBLE(N)*FAC | |
| 145 IF (R0 .EQ. 0.0D0) R0 = 1.0D0 | |
| 146 DO 560 J = 1,MBA | |
| 147 DO 530 I = J,N,MBAND | |
| 148 YI = Y(I) | |
| 149 R = DMAX1(SRUR*DABS(YI),R0/EWT(I)) | |
| 150 530 Y(I) = Y(I) + R | |
| 151 IERR = 0 | |
| 152 CALL F (NEQ, TN, Y, FTEM, IERR) | |
| 153 IF (IERR .LT. 0) RETURN | |
| 154 DO 550 JJ = J,N,MBAND | |
| 155 Y(JJ) = YH(JJ,1) | |
| 156 YJJ = Y(JJ) | |
| 157 R = DMAX1(SRUR*DABS(YJJ),R0/EWT(JJ)) | |
| 158 FAC = -HL0/R | |
| 159 I1 = MAX0(JJ-MU,1) | |
| 160 I2 = MIN0(JJ+ML,N) | |
| 161 II = JJ*MEB1 - ML + 2 | |
| 162 DO 540 I = I1,I2 | |
| 163 540 WM(II+I) = (FTEM(I) - SAVF(I))*FAC | |
| 164 550 CONTINUE | |
| 165 560 CONTINUE | |
| 166 NFE = NFE + MBA | |
| 167 C ADD IDENTITY MATRIX. ------------------------------------------------- | |
| 168 570 II = MBAND + 2 | |
| 169 DO 580 I = 1,N | |
| 170 WM(II) = WM(II) + 1.0D0 | |
| 171 580 II = II + MEBAND | |
| 172 C DO LU DECOMPOSITION OF P. -------------------------------------------- | |
| 173 CALL DGBFA (WM(3), MEBAND, N, ML, MU, IWM(21), IER) | |
| 174 IF (IER .NE. 0) IERPJ = 1 | |
| 175 RETURN | |
| 176 C----------------------- END OF SUBROUTINE PREPJ ----------------------- | |
| 177 END |