|
3983
|
1 SUBROUTINE ODESSA_PREPJ (NEQ, Y, YH, NYH, WM, IWM, EWT, SAVF, |
|
|
2 1 FTEM, PAR, F, JAC, JOPT) |
|
|
3 IMPLICIT DOUBLE PRECISION (A-H,O-Z) |
|
|
4 DIMENSION NEQ(*), Y(*), YH(NYH,*), WM(*), IWM(*), EWT(*), |
|
|
5 1 SAVF(*), FTEM(*), PAR(*) |
|
|
6 EXTERNAL F, JAC |
|
|
7 PARAMETER (ZERO=0.0D0,ONE=1.0D0) |
|
|
8 COMMON /ODE001/ ROWND, ROWNS(173), |
|
|
9 2 RDUM1(37), EL0, H, RDUM2(4), TN, UROUND, |
|
|
10 3 IOWND(14), IOWNS(4), |
|
|
11 4 IDUM1(3), IERPJ, IDUM2, JCUR, IDUM3(4), |
|
|
12 5 MITER, IDUM4(4), N, IDUM5(2), NFE, NJE, IDUM6 |
|
|
13 C----------------------------------------------------------------------- |
|
|
14 C ODESSA_PREPJ IS CALLED BY ODESSA_STODE TO COMPUTE AND PROCESS THE MATRIX |
|
|
15 C P = I - H*EL(1)*J , WHERE J IS AN APPROXIMATION TO THE JACOBIAN. |
|
|
16 C IF ISOPT = 1, ODESSA_PREPJ IS ALSO CALLED BY ODESSA_SPRIME WITH JOPT = 1. |
|
|
17 C HERE J IS COMPUTED BY THE USER-SUPPLIED ROUTINE JAC IF |
|
|
18 C MITER = 1 OR 4, OR BY FINITE DIFFERENCING IF MITER = 2, 3, OR 5. |
|
|
19 C IF MITER = 3, A DIAGONAL APPROXIMATION TO J IS USED. |
|
|
20 C J IS STORED IN WM AND REPLACED BY P. IF MITER .NE. 3, P IS THEN |
|
|
21 C SUBJECTED TO LU DECOMPOSITION (JOPT = 0) IN PREPARATION FOR LATER |
|
|
22 C SOLUTION OF LINEAR SYSTEMS WITH P AS COEFFICIENT MATRIX. THIS IS |
|
4329
|
23 C DONE BY DGETRF IF MITER = 1 OR 2, AND BY DGBTRF IF MITER = 4 OR 5. |
|
3983
|
24 C |
|
|
25 C IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION |
|
|
26 C WITH ODESSA_PREPJ USES THE FOLLOWING.. |
|
|
27 C Y = ARRAY CONTAINING PREDICTED VALUES ON ENTRY. |
|
|
28 C FTEM = WORK ARRAY OF LENGTH N (ACOR IN ODESSA_STODE). |
|
|
29 C SAVF = ARRAY CONTAINING F EVALUATED AT PREDICTED Y. |
|
|
30 C WM = REAL WORK SPACE FOR MATRICES. ON OUTPUT IT CONTAINS THE |
|
|
31 C INVERSE DIAGONAL MATRIX IF MITER = 3 AND THE LU DECOMPOSITION |
|
|
32 C OF P IF MITER IS 1, 2 , 4, OR 5. |
|
|
33 C STORAGE OF MATRIX ELEMENTS STARTS AT WM(3). |
|
|
34 C WM ALSO CONTAINS THE FOLLOWING MATRIX-RELATED DATA.. |
|
|
35 C WM(1) = SQRT(UROUND), USED IN NUMERICAL JACOBIAN INCREMENTS. |
|
|
36 C WM(2) = H*EL0, SAVED FOR LATER USE IF MITER = 3. |
|
|
37 C IWM = INTEGER WORK SPACE CONTAINING PIVOT INFORMATION, STARTING AT |
|
|
38 C IWM(21), IF MITER IS 1, 2, 4, OR 5. IWM ALSO CONTAINS BAND |
|
|
39 C PARAMETERS ML = IWM(1) AND MU = IWM(2) IF MITER IS 4 OR 5. |
|
|
40 C EL0 = EL(1) (INPUT). |
|
|
41 C IERPJ = OUTPUT ERROR FLAG, = 0 IF NO TROUBLE, .GT. 0 IF |
|
|
42 C P MATRIX FOUND TO BE SINGULAR. |
|
|
43 C JCUR = OUTPUT FLAG = 1 TO INDICATE THAT THE JACOBIAN MATRIX |
|
|
44 C (OR APPROXIMATION) IS NOW CURRENT. |
|
|
45 C JOPT = INPUT JACOBIAN OPTION, = 1 IF JAC IS DESIRED ONLY. |
|
|
46 C THIS ROUTINE ALSO USES THE COMMON VARIABLES EL0, H, TN, UROUND, |
|
|
47 C IERPJ, JCUR, MITER, N, NFE, AND NJE. |
|
|
48 C----------------------------------------------------------------------- |
|
|
49 NJE = NJE + 1 |
|
|
50 IERPJ = 0 |
|
|
51 JCUR = 1 |
|
|
52 HL0 = H*EL0 |
|
|
53 GO TO (100, 200, 300, 400, 500), MITER |
|
|
54 C IF MITER = 1, CALL JAC AND MULTIPLY BY SCALAR. ----------------------- |
|
|
55 100 LENP = N*N |
|
|
56 DO 110 I = 1,LENP |
|
|
57 110 WM(I+2) = ZERO |
|
|
58 CALL JAC (NEQ, TN, Y, PAR, 0, 0, WM(3), N) |
|
|
59 IF (JOPT .EQ. 1) RETURN |
|
|
60 CON = -HL0 |
|
|
61 DO 120 I = 1,LENP |
|
|
62 120 WM(I+2) = WM(I+2)*CON |
|
|
63 GO TO 240 |
|
|
64 C IF MITER = 2, MAKE N CALLS TO F TO APPROXIMATE J. -------------------- |
|
|
65 200 FAC = ODESSA_VNORM (N, SAVF, EWT) |
|
3987
|
66 R0 = 1000.0D0*DABS(H)*UROUND*DBLE(N)*FAC |
|
3983
|
67 IF (R0 .EQ. ZERO) R0 = ONE |
|
|
68 SRUR = WM(1) |
|
|
69 J1 = 2 |
|
|
70 DO 230 J = 1,N |
|
|
71 YJ = Y(J) |
|
3987
|
72 R = DMAX1(SRUR*DABS(YJ),R0/EWT(J)) |
|
3983
|
73 Y(J) = Y(J) + R |
|
|
74 FAC = -HL0/R |
|
|
75 CALL F (NEQ, TN, Y, PAR, FTEM) |
|
|
76 DO 220 I = 1,N |
|
|
77 220 WM(I+J1) = (FTEM(I) - SAVF(I))*FAC |
|
|
78 Y(J) = YJ |
|
|
79 J1 = J1 + N |
|
|
80 230 CONTINUE |
|
|
81 NFE = NFE + N |
|
|
82 IF (JOPT .EQ. 1) RETURN |
|
|
83 C ADD IDENTITY MATRIX. ------------------------------------------------- |
|
|
84 240 J = 3 |
|
|
85 DO 250 I = 1,N |
|
|
86 WM(J) = WM(J) + ONE |
|
|
87 250 J = J + (N + 1) |
|
|
88 C DO LU DECOMPOSITION ON P. -------------------------------------------- |
|
4329
|
89 CALL DGETRF ( N, N, WM(3), N, IWM(21), IER) |
|
3983
|
90 IF (IER .NE. 0) IERPJ = 1 |
|
|
91 RETURN |
|
|
92 C IF MITER = 3, CONSTRUCT A DIAGONAL APPROXIMATION TO J AND P. --------- |
|
|
93 300 WM(2) = HL0 |
|
|
94 R = EL0*0.1D0 |
|
|
95 DO 310 I = 1,N |
|
|
96 310 Y(I) = Y(I) + R*(H*SAVF(I) - YH(I,2)) |
|
|
97 CALL F (NEQ, TN, Y, PAR, WM(3)) |
|
|
98 NFE = NFE + 1 |
|
|
99 DO 320 I = 1,N |
|
|
100 R0 = H*SAVF(I) - YH(I,2) |
|
|
101 DI = 0.1D0*R0 - H*(WM(I+2) - SAVF(I)) |
|
|
102 WM(I+2) = 1.0D0 |
|
3987
|
103 IF (DABS(R0) .LT. UROUND/EWT(I)) GO TO 320 |
|
|
104 IF (DABS(DI) .EQ. ZERO) GO TO 330 |
|
3983
|
105 WM(I+2) = 0.1D0*R0/DI |
|
|
106 320 CONTINUE |
|
|
107 RETURN |
|
|
108 330 IERPJ = 1 |
|
|
109 RETURN |
|
|
110 C IF MITER = 4, CALL JAC AND MULTIPLY BY SCALAR. ----------------------- |
|
|
111 400 ML = IWM(1) |
|
|
112 MU = IWM(2) |
|
|
113 ML3 = ML + 3 |
|
|
114 MBAND = ML + MU + 1 |
|
|
115 MEBAND = MBAND + ML |
|
|
116 LENP = MEBAND*N |
|
|
117 DO 410 I = 1,LENP |
|
|
118 410 WM(I+2) = ZERO |
|
|
119 CALL JAC (NEQ, TN, Y, PAR, ML, MU, WM(ML3), MEBAND) |
|
|
120 IF (JOPT .EQ. 1) RETURN |
|
|
121 CON = -HL0 |
|
|
122 DO 420 I = 1,LENP |
|
|
123 420 WM(I+2) = WM(I+2)*CON |
|
|
124 GO TO 570 |
|
|
125 C IF MITER = 5, MAKE MBAND CALLS TO F TO APPROXIMATE J. ---------------- |
|
|
126 500 ML = IWM(1) |
|
|
127 MU = IWM(2) |
|
|
128 MBAND = ML + MU + 1 |
|
3987
|
129 MBA = MIN0(MBAND,N) |
|
3983
|
130 MEBAND = MBAND + ML |
|
|
131 MEB1 = MEBAND - 1 |
|
|
132 SRUR = WM(1) |
|
|
133 FAC = ODESSA_VNORM (N, SAVF, EWT) |
|
3987
|
134 R0 = 1000.0D0*DABS(H)*UROUND*DBLE(N)*FAC |
|
3983
|
135 IF (R0 .EQ. ZERO) R0 = ONE |
|
|
136 DO 560 J = 1,MBA |
|
|
137 DO 530 I = J,N,MBAND |
|
|
138 YI = Y(I) |
|
3987
|
139 R = DMAX1(SRUR*DABS(YI),R0/EWT(I)) |
|
3983
|
140 530 Y(I) = Y(I) + R |
|
|
141 CALL F (NEQ, TN, Y, PAR, FTEM) |
|
|
142 DO 550 JJ = J,N,MBAND |
|
|
143 Y(JJ) = YH(JJ,1) |
|
|
144 YJJ = Y(JJ) |
|
3987
|
145 R = DMAX1(SRUR*DABS(YJJ),R0/EWT(JJ)) |
|
3983
|
146 FAC = -HL0/R |
|
3987
|
147 I1 = MAX0(JJ-MU,1) |
|
|
148 I2 = MIN0(JJ+ML,N) |
|
3983
|
149 II = JJ*MEB1 - ML + 2 |
|
|
150 DO 540 I = I1,I2 |
|
|
151 540 WM(II+I) = (FTEM(I) - SAVF(I))*FAC |
|
|
152 550 CONTINUE |
|
|
153 560 CONTINUE |
|
|
154 NFE = NFE + MBA |
|
|
155 IF (JOPT .EQ. 1) RETURN |
|
|
156 C ADD IDENTITY MATRIX. ------------------------------------------------- |
|
|
157 570 II = MBAND + 2 |
|
|
158 DO 580 I = 1,N |
|
|
159 WM(II) = WM(II) + ONE |
|
|
160 580 II = II + MEBAND |
|
|
161 C DO LU DECOMPOSITION OF P. -------------------------------------------- |
|
4329
|
162 CALL DGBTRF ( N, N, ML, MU, WM(3), MEBAND, IWM(21), IER) |
|
3983
|
163 IF (IER .NE. 0) IERPJ = 1 |
|
|
164 RETURN |
|
|
165 C---------------- END OF SUBROUTINE ODESSA_PREPJ ----------------------- |
|
|
166 END |
