|
4583
|
1 C----------------------------------------------------------------------- |
|
|
2 C----------------------------------------------------------------------- |
|
|
3 C THIS IS THE SEPTEMBER 1, 1986 VERSION OF ODESSA.. |
|
|
4 C AN ORDINARY DIFFERENTIAL EQUATION SOLVER WITH EXPLICIT SIMULTANEOUS |
|
|
5 C SENSITIVITY ANALYSIS. |
|
|
6 C |
|
|
7 C THIS PACKAGE IS A MODIFICATION OF THE AUGUST 13, 1981 VERSION OF |
|
|
8 C LSODE.. LIVERMORE SOLVER FOR ORDINARY DIFFERENTIAL EQUATIONS. |
|
|
9 C THIS VERSION IS IN DOUBLE PRECISION. |
|
|
10 C |
|
|
11 C ODESSA SOLVES FOR THE FIRST-ORDER SENSITIVITY COEFFICIENTS.. |
|
|
12 C DY(I)/DP, FOR A SINGLE PARAMETER, OR, |
|
|
13 C DY(I)/DP(J), FOR MULTIPLE PARAMETERS, |
|
|
14 C ASSOCIATED WITH A SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS.. |
|
|
15 C DY/DT = F(Y,T;P). |
|
|
16 C----------------------------------------------------------------------- |
|
|
17 C REFERENCES... |
|
|
18 C |
|
|
19 C 1. JORGE R. LEIS AND MARK A. KRAMER, THE SIMULTANEOUS SOLUTION AND |
|
|
20 C EXPLICIT SENSITIVITY ANALYSIS OF SYSTEMS DESCRIBED BY ORDINARY |
|
|
21 C DIFFERENTIAL EQUATIONS. SUBMITTED TO ACM TRANS. MATH. SOFTWARE, |
|
|
22 C (1985). |
|
|
23 C |
|
|
24 C 2. JORGE R. LEIS AND MARK A. KRAMER, ODESSA - AN ORDINARY DIFFERENTIA |
|
|
25 C EQUATION SOLVER WITH EXPLICIT SIMULTANEOUS SENSITIVITY ANALYSIS. |
|
|
26 C SUBMITTED TO ACM TRANS. MATH. SOFTWARE, (1985). |
|
|
27 C |
|
|
28 C 3. ALAN C. HINDMARSH, LSODE AND LSODI, TWO NEW INITIAL VALUE |
|
|
29 C ORDINARY DIFFERENTIAL EQUATION SOLVERS, ACM-SIGNUM NEWSLETTER, |
|
|
30 C VOL. 15, NO. 4 (1980), PP. 10-11. |
|
|
31 C----------------------------------------------------------------------- |
|
|
32 C PROBLEM STATEMENT.. |
|
|
33 C |
|
|
34 C THE ODESSA MODIFICATION OF THE LSODE PACKAGE PROVIDES THE OPTION TO |
|
|
35 C CALCULATE FIRST-ORDER SENSITIVITY COEFFICIENTS FOR A SYSTEM OF STIFF |
|
|
36 C OR NON-STIFF EXPLICIT ORDINARY DIFFERENTIAL EQUATIONS OF THE GENERAL |
|
|
37 C FORM : |
|
|
38 C |
|
|
39 C DY/DT = F(Y,T;P) (1) |
|
|
40 C |
|
|
41 C WHERE Y IS AN N-DIMENSIONAL DEPENDENT VARIABLE VECTOR, T IS THE |
|
|
42 C INDEPENDENT INTEGRATION VARIABLE, AND P IS AN NPAR-DIMENSIONAL |
|
|
43 C CONSTANT VECTOR. THE GOVERNING EQUATIONS FOR THE FIRST-ORDER |
|
|
44 C SENSITIVITY COEFFICIENTS ARE GIVEN BY : |
|
|
45 C |
|
|
46 C S'(T) = J(T)*S(T) + DF/DP (2) |
|
|
47 C |
|
|
48 C WHERE |
|
|
49 C |
|
|
50 C S(T) = DY(T)/DP (= SENSITIVITY FUNCTIONS) |
|
|
51 C S'(T) = D(DY(T)/DP)/DT |
|
|
52 C J(T) = DF(Y,T;P)/DY(T) (= JACOBIAN MATRIX) |
|
|
53 C AND DF/DP = DF(Y,T;P)/DP (= INHOMOGENEITY MATRIX) |
|
|
54 C |
|
|
55 C SOLUTION OF EQUATIONS (1) AND (2) PROCEEDS SIMULTANEOUSLY VIA AN |
|
|
56 C EXTENSION OF THE LSODE PACKAGE AS DESCRIBED IN [1]. |
|
|
57 C---------------------------------------------------------------------- |
|
|
58 C ACKNOWLEDGEMENT : THE FOLLOWING ODESSA PACKAGE DOCUMENTATION IS A |
|
|
59 C MODIFICATION OF THE LSODE DOCUMENTATION WHICH |
|
|
60 C ACCOMPANIES THE LSODE PACKAGE CODE. |
|
|
61 C---------------------------------------------------------------------- |
|
|
62 C SUMMARY OF USAGE. |
|
|
63 C |
|
|
64 C COMMUNICATION BETWEEN THE USER AND THE ODESSA PACKAGE, FOR NORMAL |
|
|
65 C SITUATIONS, IS SUMMARIZED HERE. THIS SUMMARY DESCRIBES ONLY A SUBSET |
|
|
66 C OF THE FULL SET OF OPTIONS AVAILABLE. SEE THE FULL DESCRIPTION FOR |
|
|
67 C DETAILS, INCLUDING OPTIONAL COMMUNICATION, NONSTANDARD OPTIONS, |
|
|
68 C AND INSTRUCTIONS FOR SPECIAL SITUATIONS. SEE ALSO THE EXAMPLE |
|
|
69 C PROBLEM (WITH PROGRAM AND OUTPUT) FOLLOWING THIS SUMMARY. |
|
|
70 C |
|
|
71 C A. FIRST PROVIDE A SUBROUTINE OF THE FORM.. |
|
|
72 C SUBROUTINE F (N, T, Y, PAR, YDOT) |
|
|
73 C DOUBLE PRECISION T, Y, PAR, YDOT |
|
|
74 C DIMENSION Y(N), YDOT(N), PAR(NPAR) |
|
|
75 C WHICH SUPPLIES THE VECTOR FUNCTION F BY LOADING YDOT(I) WITH F(I). |
|
|
76 C N IS THE NUMBER OF FIRST-ORDER DIFFERENTIAL EQUATIONS IN THE |
|
|
77 C ABOVE MODEL. NPAR IS THE NUMBER OF MODEL PARAMETERS FOR WHICH |
|
|
78 C VECTOR SENSITIVITY FUNCTIONS ARE DESIRED. YOU ARE ALSO ENCOURAGED |
|
|
79 C TO PROVIDE A SUBROUTINE OF THE FORM.. |
|
|
80 C SUBROUTINE DF (N, T, Y, PAR, DFDP, JPAR) |
|
|
81 C DOUBLE PRECISION T, Y, PAR, DFDP |
|
|
82 C DIMENSION Y(N), PAR(NPAR), DFDP(N) |
|
|
83 C GO TO (1,...,NPAR) JPAR |
|
|
84 C 1 DFDP(1) = DF(1)/DP(1) |
|
|
85 C . |
|
|
86 C DFDP(I) = DF(I)/DP(1) |
|
|
87 C . |
|
|
88 C DFDP(N) = DF(N)/DP(1) |
|
|
89 C RETURN |
|
|
90 C 2 DFDP(1) = DF(1)/DP(2) |
|
|
91 C . |
|
|
92 C DFDP(I) = DF(I)/DP(2) |
|
|
93 C . |
|
|
94 C DFDP(N) = DF(N)/DP(2) |
|
|
95 C RETURN |
|
|
96 C . . |
|
|
97 C . . |
|
|
98 C RETURN |
|
|
99 C NPAR DFDP(1) = DF(1)/DP(NPAR) |
|
|
100 C . |
|
|
101 C DFDP(I) = DF(I)/DP(NPAR) |
|
|
102 C . |
|
|
103 C DFDP(N) = DF(N)/DP(NPAR) |
|
|
104 C RETURN |
|
|
105 C END |
|
|
106 C ONLY NONZERO ELEMENTS NEED BE LOADED. IF THIS IS NOT FEASIBLE, |
|
|
107 C ODESSA WILL GENERATE THIS MATRIX INTERNALLY BY DIFFERENCE QUOTIENTS. |
|
|
108 C |
|
|
109 C B. NEXT DETERMINE (OR GUESS) WHETHER OR NOT THE PROBLEM IS STIFF. |
|
|
110 C STIFFNESS OCCURS WHEN THE JACOBIAN MATRIX DF/DY HAS AN EIGENVALUE |
|
|
111 C WHOSE REAL PART IS NEGATIVE AND LARGE IN MAGNITUDE, COMPARED TO THE |
|
|
112 C RECIPROCAL OF THE T SPAN OF INTEREST. IF THE PROBLEM IS NONSTIFF, |
|
|
113 C USE METH = 10. IF IT IS STIFF, USE METH = 20. THE USER IS REQUIRED |
|
|
114 C TO INPUT THE METHOD FLAG MF = 10*METH + MITER. THERE ARE FOUR |
|
|
115 C STANDARD CHOICES FOR MITER WHEN A SENSITIVITY ANALYSIS IS DESIRED, |
|
|
116 C AND ODESSA REQUIRES THE JACOBIAN MATRIX IN SOME FORM. |
|
|
117 C THIS MATRIX IS REGARDED EITHER AS FULL (MITER = 1 OR 2), |
|
|
118 C OR BANDED (MITER = 4 OR 5). IN THE BANDED CASE, ODESSA REQUIRES TWO |
|
|
119 C HALF-BANDWIDTH PARAMETERS ML AND MU. THESE ARE, RESPECTIVELY, THE |
|
|
120 C WIDTHS OF THE LOWER AND UPPER PARTS OF THE BAND, EXCLUDING THE MAIN |
|
|
121 C DIAGONAL. THUS THE BAND CONSISTS OF THE LOCATIONS (I,J) WITH |
|
|
122 C I-ML .LE. J .LE. I+MU, AND THE FULL BANDWIDTH IS ML+MU+1. |
|
|
123 C |
|
|
124 C C. YOU ARE ENCOURAGED TO SUPPLY THE JACOBIAN DIRECTLY (MF = 11, 14, |
|
|
125 C 21, OR 24), BUT IF THIS IS NOT FEASIBLE, ODESSA WILL COMPUTE IT |
|
|
126 C INTERNALLY BY DIFFERENCE QUOTIENTS (MF = 12, 15, 22, OR 25). IF YOU |
|
|
127 C ARE SUPPLYING THE JACOBIAN, PROVIDE A SUBROUTINE OF THE FORM.. |
|
|
128 C SUBROUTINE JAC (NEQ, T, Y, PAR, ML, MU, PD, NROWPD) |
|
|
129 C DOUBLE PRECISION T, Y, PAR, PD |
|
|
130 C DIMENSION Y(N), PD(NROWPD,N), PAR(NPAR) |
|
|
131 C WHICH SUPPLIES DF/DY BY LOADING PD AS FOLLOWS.. |
|
|
132 C FOR A FULL JACOBIAN (MF = 11, OR 21), LOAD PD(I,J) WITH DF(I)/DY(J), |
|
|
133 C THE PARTIAL DERIVATIVE OF F(I) WITH RESPECT TO Y(J). (IGNORE THE |
|
|
134 C ML AND MU ARGUMENTS IN THIS CASE.) |
|
|
135 C FOR A BANDED JACOBIAN (MF = 14, OR 24), LOAD PD(I-J+MU+1,J) WITH |
|
|
136 C DF(I)/DY(J), I.E. LOAD THE DIAGONAL LINES OF DF/DY INTO THE ROWS OF |
|
|
137 C PD FROM THE TOP DOWN. |
|
|
138 C IN EITHER CASE, ONLY NONZERO ELEMENTS NEED BE LOADED. |
|
|
139 C |
|
|
140 C D. WRITE A MAIN PROGRAM WHICH CALLS SUBROUTINE ODESSA ONCE FOR |
|
|
141 C EACH POINT AT WHICH ANSWERS ARE DESIRED. THIS SHOULD ALSO PROVIDE |
|
|
142 C FOR POSSIBLE USE OF LOGICAL UNIT 6 FOR OUTPUT OF ERROR MESSAGES BY |
|
|
143 C ODESSA. ON THE FIRST CALL TO ODESSA, SUPPLY ARGUMENTS AS FOLLOWS.. |
|
|
144 C F = NAME OF SUBROUTINE FOR RIGHT-HAND SIDE VECTOR F (MODEL). |
|
|
145 C THIS NAME MUST BE DECLARED EXTERNAL IN CALLING PROGRAM. |
|
|
146 C DF = NAME OF SUBROUTINE FOR INHOMOGENEITY MATRIX DF/DP. |
|
|
147 C IF USED (IDF = 1), THIS NAME MUST BE DECLARED EXTERNAL IN |
|
|
148 C CALLING PROGRAM. IF NOT USED (IDF = 0), PASS A DUMMY NAME. |
|
|
149 C N = NUMBER OF FIRST ORDER ODE-S IN MODEL; LOAD INTO NEQ(1). |
|
|
150 C NPAR = NUMBER OF MODEL PARAMETERS OF INTEREST; LOAD INTO NEQ(2). |
|
|
151 C Y = AN (N) BY (NPAR+1) REAL ARRAY OF INITIAL VALUES.. |
|
|
152 C Y(I,1) , I = 1,N , CONTAIN THE STATE, OR MODEL, DEPENDENT |
|
|
153 C VARIABLES, |
|
|
154 C Y(I,J) , J = 2,NPAR , CONTAIN THE DEPENDENT SENSITIVITY |
|
|
155 C COEFFICIENTS. |
|
|
156 C PAR = A REAL ARRAY OF LENGTH NPAR CONTAINING MODEL PARAMETERS |
|
|
157 C OF INTEREST. |
|
|
158 C T = THE INITIAL VALUE OF THE INDEPENDENT VARIABLE. |
|
|
159 C TOUT = FIRST POINT WHERE OUTPUT IS DESIRED (.NE. T). |
|
|
160 C ITOL = 1, 2, 3, OR 4 ACCORDING AS RTOL, ATOL (BELOW) ARE SCALARS |
|
|
161 C OR ARRAYS. |
|
|
162 C RTOL = RELATIVE TOLERANCE PARAMETER (SCALAR OR (N) BY (NPAR+1) |
|
|
163 C ARRAY). |
|
|
164 C ATOL = ABSOLUTE TOLERANCE PARAMETER (SCALAR OR (N) BY (NPAR+1) |
|
|
165 C ARRAY). |
|
|
166 C THE ESTIMATED LOCAL ERROR IN Y(I,J) WILL BE CONTROLLED SO AS |
|
|
167 C TO BE ROUGHLY LESS (IN MAGNITUDE) THAN |
|
|
168 C EWT(I,J) = RTOL*ABS(Y(I,J)) + ATOL IF ITOL = 1, |
|
|
169 C EWT(I,J) = RTOL*ABS(Y(I,J)) + ATOL(I,J) IF ITOL = 2, |
|
|
170 C EWT(I,J) = RTOL(I,J)*ABS(Y(I,J) + ATOL IF ITOL = 3, OR |
|
|
171 C EWT(I,J) = RTOL(I,J)*ABS(Y(I,J) + ATOL(I,J) IF ITOL = 4. |
|
|
172 C THUS THE LOCAL ERROR TEST PASSES IF, IN EACH COMPONENT, |
|
|
173 C EITHER THE ABSOLUTE ERROR IS LESS THAN ATOL (OR ATOL(I,J)), |
|
|
174 C OR THE RELATIVE ERROR IS LESS THAN RTOL (OR RTOL(I,J)). |
|
|
175 C USE RTOL = 0.0 FOR PURE ABSOLUTE ERROR CONTROL, AND |
|
|
176 C USE ATOL = 0.0 FOR PURE RELATIVE ERROR CONTROL. |
|
|
177 C CAUTION.. ACTUAL (GLOBAL) ERRORS MAY EXCEED THESE LOCAL |
|
|
178 C TOLERANCES, SO CHOOSE THEM CONSERVATIVELY. |
|
|
179 C ITASK = 1 FOR NORMAL COMPUTATION OF OUTPUT VALUES OF Y AT T = TOUT. |
|
|
180 C ISTATE = INTEGER FLAG (INPUT AND OUTPUT). SET ISTATE = 1. |
|
|
181 C IOPT = 0, TO INDICATE NO OPTIONAL INPUTS FOR INTEGRATION; |
|
|
182 C LOAD INTO IOPT(1). |
|
|
183 C ISOPT = 1, TO INDICATE SENSITIVITY ANALYSIS, = 0, TO INDICATE |
|
|
184 C NO SENSITIVITY ANALYSIS; LOAD INTO IOPT(2). |
|
|
185 C IDF = 1, IF SUBROUTINE DF (ABOVE) IS SUPPLIED BY THE USER, |
|
|
186 C = 0, OTHERWISE; LOAD INTO IOPT(3). |
|
|
187 C RWORK = REAL WORK ARRAY OF LENGTH AT LEAST.. |
|
|
188 C 22 + 16*N + N**2 FOR MF = 11 OR 12, |
|
|
189 C 22 + 17*N + (2*ML + MU)*N FOR MF = 14 OR 15, |
|
|
190 C 22 + 9*N + N**2 FOR MF = 21 OR 22, |
|
|
191 C 22 + 10*N + (2*ML + MU)*N FOR MF = 24 OR 25, |
|
|
192 C IF ISOPT = 0, OR.. |
|
|
193 C 22 + 15*(NPAR+1)*N + N**2 + N FOR MF = 11 OR 12, |
|
|
194 C 24 + 15*(NPAR+1)*N + (2*ML+MU+2)*N + N FOR MF = 14 OR 15, |
|
|
195 C 22 + 8*(NPAR+1)*N + N**2 + N FOR MF = 21 OR 22, |
|
|
196 C 24 + 8*(NPAR+1)*N + (2*ML+MU+2)*N + N FOR MF = 24 OR 25, |
|
|
197 C IF ISOPT = 1. |
|
|
198 C LRW = DECLARED LENGTH OF RWORK (IN USER-S DIMENSION STATEMENT). |
|
|
199 C IWORK = INTEGER WORK ARRAY OF LENGTH AT LEAST.. |
|
|
200 C 20 + N IF ISOPT = 0, |
|
|
201 C 21 + N + NPAR IF ISOPT = 1. |
|
|
202 C IF MITER = 4 OR 5, INPUT IN IWORK(1),IWORK(2) THE LOWER |
|
|
203 C AND UPPER HALF-BANDWIDTHS ML,MU (EXCLUDING MAIN DIAGONAL). |
|
|
204 C LIW = DECLARED LENGTH OF IWORK (IN USER-S DIMENSION STATEMENT). |
|
|
205 C JAC = NAME OF SUBROUTINE FOR JACOBIAN MATRIX. |
|
|
206 C IF USED, THIS NAME MUST BE DECLARED EXTERNAL IN CALLING |
|
|
207 C PROGRAM. IF NOT USED, PASS A DUMMY NAME. |
|
|
208 C MF = METHOD FLAG. STANDARD VALUES FOR ISOPT = 0 ARE.. |
|
|
209 C 10 FOR NONSTIFF (ADAMS) METHOD, NO JACOBIAN USED. |
|
|
210 C 21 FOR STIFF (BDF) METHOD, USER-SUPPLIED FULL JACOBIAN. |
|
|
211 C 22 FOR STIFF METHOD, INTERNALLY GENERATED FULL JACOBIAN. |
|
|
212 C 24 FOR STIFF METHOD, USER-SUPPLIED BANDED JACOBIAN. |
|
|
213 C 25 FOR STIFF METHOD, INTERNALLY GENERATED BANDED JACOBIAN. |
|
|
214 C IF ISOPT = 1, MF = 10 IS ILLEGAL AND CAN BE REPLACED BY.. |
|
|
215 C 11 FOR NONSTIFF METHOD, USER-SUPPLIED FULL JACOBIAN. |
|
|
216 C 12 FOR NONSTIFF METHOD, INTERNALLY GENERATED FULL JACOBIAN. |
|
|
217 C 14 FOR NONSTIFF METHOD, USER-SUPPLIED BANDED JACOBIAN. |
|
|
218 C 15 FOR NONSTIFF METHOD, INTERNALLY GENERATED BANDED JACOBIAN. |
|
|
219 C NOTE THAT THE MAIN PROGRAM MUST DECLARE ARRAYS Y, RWORK, IWORK, AND |
|
|
220 C POSSIBLY ATOL AND RTOL, AS WELL AS NEQ, IOPT, AND PAR IF ISOPT = 1. |
|
|
221 C |
|
|
222 C E. THE OUTPUT FROM THE FIRST CALL (OR ANY CALL) IS.. |
|
|
223 C Y = ARRAY OF COMPUTED VALUES OF Y(T) VECTOR. |
|
|
224 C T = CORRESPONDING VALUE OF INDEPENDENT VARIABLE (NORMALLY TOUT). |
|
|
225 C ISTATE = 2 IF ODESSA WAS SUCCESSFUL, NEGATIVE OTHERWISE. |
|
|
226 C -1 MEANS EXCESS WORK DONE ON THIS CALL (PERHAPS WRONG MF). |
|
|
227 C -2 MEANS EXCESS ACCURACY REQUESTED (TOLERANCES TOO SMALL). |
|
|
228 C -3 MEANS ILLEGAL INPUT DETECTED (SEE PRINTED MESSAGE). |
|
|
229 C -4 MEANS REPEATED ERROR TEST FAILURES (CHECK ALL INPUTS). |
|
|
230 C -5 MEANS REPEATED CONVERGENCE FAILURES (PERHAPS BAD JACOBIAN |
|
|
231 C SUPPLIED OR WRONG CHOICE OF MF OR TOLERANCES). |
|
|
232 C -6 MEANS ERROR WEIGHT BECAME ZERO DURING PROBLEM. (SOLUTION |
|
|
233 C COMPONENT I,J VANISHED, AND ATOL OR ATOL(I,J) = 0.0) |
|
|
234 C |
|
|
235 C F. TO CONTINUE THE INTEGRATION AFTER A SUCCESSFUL RETURN, SIMPLY |
|
|
236 C RESET TOUT AND CALL ODESSA AGAIN. NO OTHER PARAMETERS NEED BE RESET. |
|
|
237 C---------------------------------------------------------------------- |
|
|
238 C EXAMPLE PROBLEM. |
|
|
239 C |
|
|
240 C THE FOLLOWING IS A SIMPLE EXAMPLE PROBLEM, WITH THE CODING |
|
|
241 C NEEDED FOR ITS SOLUTION BY ODESSA. THE PROBLEM IS FROM CHEMICAL |
|
|
242 C KINETICS, AND CONSISTS OF THE FOLLOWING THREE RATE EQUATIONS.. |
|
|
243 C DY1/DT = -PAR(1)*Y1 + PAR(2)*Y2*Y3 ; PAR(1) = .04, PAR(2) = 1.E4 |
|
|
244 C DY2/DT = PAR(1)*Y1 - PAR(2)*Y2*Y3 - PAR(3)*Y2**2 ; PAR(3) = 3.E7 |
|
|
245 C DY3/DT = PAR(3)*Y2**2 |
|
|
246 C ON THE INTERVAL FROM T = 0.0 TO T = 4.E10, WITH INITIAL CONDITIONS |
|
|
247 C Y1 = 1.0, Y2 = Y3 = 0, AND S(I,J) = 0, I = 1,3, J = 1,3. |
|
|
248 C THE PROBLEM IS STIFF. |
|
|
249 C |
|
|
250 C THE FOLLOWING CODING SOLVES THIS PROBLEM WITH ODESSA, USING |
|
|
251 C MF = 21 AND PRINTING RESULTS AT T = .4, 4., ..., 4.E10. |
|
|
252 C IT USES ITOL = 4 AND ATOL MUCH SMALLER FOR Y2 THAN Y1 OR Y3, |
|
|
253 C BECAUSE Y2 HAS MUCH SMALLER VALUES. LESS STRINGENT TOLERANCES |
|
|
254 C ARE ASSIGNED FOR THE SENSITIVITIES TO ACHIEVE GREATER EFFICIENCY. |
|
|
255 C AT THE END OF THE RUN, STATISTICAL QUANTITIES OF INTEREST ARE |
|
|
256 C PRINTED (SEE OPTIONAL OUTPUTS IN THE FULL DESCRIPTION BELOW). |
|
|
257 C |
|
|
258 C DOUBLE PRECISION ATOL, RWORK, RTOL, T, TOUT, Y, PAR |
|
|
259 C EXTERNAL FEX, JEX, DFEX |
|
|
260 C DIMENSION Y(3,4), PAR(3), ATOL(3,4), RTOL(3,4), RWORK(130), |
|
|
261 C 1 IWORK(27), NEQ(2), IOPT(3) |
|
|
262 C N = 3 |
|
|
263 C NPAR = 3 |
|
|
264 C NEQ(1) = N |
|
|
265 C NEQ(2) = NPAR |
|
|
266 C NSV = NPAR+1 |
|
|
267 C DO 10 I = 1,N |
|
|
268 C DO 10 J = 1,NSV |
|
|
269 C 10 Y(I,J) = 0.0D0 |
|
|
270 C Y(1,1) = 1.0D0 |
|
|
271 C PAR(1) = 0.04D0 |
|
|
272 C PAR(2) = 1.0D4 |
|
|
273 C PAR(3) = 3.0D7 |
|
|
274 C T = 0.D0 |
|
|
275 C TOUT = .4D0 |
|
|
276 C ITOL = 4 |
|
|
277 C ATOL(1,1) = 1.D-6 |
|
|
278 C ATOL(2,1) = 1.D-10 |
|
|
279 C ATOL(3,1) = 1.D-6 |
|
|
280 C DO 20 I = 1,N |
|
|
281 C RTOL(I,1) = 1.D-4 |
|
|
282 C DO 15 J = 2,NSV |
|
|
283 C RTOL(I,J) = 1.D-3 |
|
|
284 C 15 ATOL(I,J) = 1.D2 * ATOL(I,1) |
|
|
285 C 20 CONTINUE |
|
|
286 C ITASK = 1 |
|
|
287 C ISTATE = 1 |
|
|
288 C IOPT(1) = 0 |
|
|
289 C IOPT(2) = 1 |
|
|
290 C IOPT(3) = 1 |
|
|
291 C LRW = 130 |
|
|
292 C LIW = 27 |
|
|
293 C MF = 21 |
|
|
294 C DO 60 IOUT = 1,12 |
|
|
295 C CALL ODESSA(FEX,DFEX,NEQ,Y,PAR,T,TOUT,ITOL,RTOL,ATOL, |
|
|
296 C 1 ITASK,ISTATE, IOPT,RWORK,LRW,IWORK,LIW,JEX,MF) |
|
|
297 C WRITE(6,30)T,Y(1,1),Y(2,1),Y(3,1) |
|
|
298 C 30 FORMAT(1X,7H AT T =,E12.4,6H Y =,3E14.6) |
|
|
299 C DO 50 J = 2,NSV |
|
|
300 C JPAR = J-1 |
|
|
301 C WRITE(6,40)JPAR,Y(1,J),Y(2,J),Y(3,J) |
|
|
302 C 40 FORMAT(20X,2HS(,I1,3H) =,3E14.6) |
|
|
303 C 50 CONTINUE |
|
|
304 C IF (ISTATE .LT. 0) GO TO 80 |
|
|
305 C 60 TOUT = TOUT*10.D0 |
|
|
306 C WRITE(6,70)IWORK(11),IWORK(12),IWORK(13),IWORK(19) |
|
|
307 C 70 FORMAT(1X,/,12H NO. STEPS =,I4,11H NO. F-S =,I4,11H NO. J-S =, |
|
|
308 C 1 I4,12H NO. DF-S =,I4) |
|
|
309 C STOP |
|
|
310 C 80 WRITE(6,90)ISTATE |
|
|
311 C 90 FORMAT(///22H ERROR HALT.. ISTATE =,I3) |
|
|
312 C STOP |
|
|
313 C END |
|
|
314 C |
|
|
315 C SUBROUTINE FEX (NEQ, T, Y, PAR, YDOT) |
|
|
316 C DOUBLE PRECISION T, Y, YDOT, PAR |
|
|
317 C DIMENSION Y(3), YDOT(3), PAR(3) |
|
|
318 C YDOT(1) = -PAR(1)*Y(1) + PAR(2)*Y(2)*Y(3) |
|
|
319 C YDOT(3) = PAR(3)*Y(2)*Y(2) |
|
|
320 C YDOT(2) = -YDOT(1) - YDOT(3) |
|
|
321 C RETURN |
|
|
322 C END |
|
|
323 C |
|
|
324 C SUBROUTINE JEX (NEQ, T, Y, PAR, ML, MU, PD, NRPD) |
|
|
325 C DOUBLE PRECISION PD, T, Y, PAR |
|
|
326 C DIMENSION Y(3), PD(NRPD,3), PAR(3) |
|
|
327 C PD(1,1) = -PAR(1) |
|
|
328 C PD(1,2) = PAR(2)*Y(3) |
|
|
329 C PD(1,3) = PAR(2)*Y(2) |
|
|
330 C PD(2,1) = PAR(1) |
|
|
331 C PD(2,3) = -PD(1,3) |
|
|
332 C PD(3,2) = 2.D0*PAR(3)*Y(2) |
|
|
333 C PD(2,2) = -PD(1,2) - PD(3,2) |
|
|
334 C RETURN |
|
|
335 C END |
|
|
336 C |
|
|
337 C SUBROUTINE DFEX (NEQ, T, Y, PAR, DFDP, JPAR) |
|
|
338 C DOUBLE PRECISION T, Y, PAR, DFDP |
|
|
339 C DIMENSION Y(3), PAR(3), DFDP(3) |
|
|
340 C GO TO (1,2,3), JPAR |
|
|
341 C 1 DFDP(1) = -Y(1) |
|
|
342 C DFDP(2) = Y(1) |
|
|
343 C RETURN |
|
|
344 C 2 DFDP(1) = Y(2)*Y(3) |
|
|
345 C DFDP(2) = -Y(2)*Y(3) |
|
|
346 C RETURN |
|
|
347 C 3 DFDP(2) = -Y(2)*Y(2) |
|
|
348 C DFDP(3) = Y(2)*Y(2) |
|
|
349 C RETURN |
|
|
350 C END |
|
|
351 C |
|
|
352 C THE OUTPUT OF THIS PROGRAM (ON A DATA GENERAL MV-8000 IN |
|
|
353 C DOUBLE PRECISION IS AS FOLLOWS: |
|
|
354 C |
|
|
355 C AT T = .4000E+00 Y = .985173E+00 .338641E-04 .147930E-01 |
|
|
356 C S(1) = -.355914E+00 .390261E-03 .355524E+00 |
|
|
357 C S(2) = .955150E-07 -.213065E-09 -.953019E-07 |
|
|
358 C S(3) = -.158466E-10 -.529012E-12 .163756E-10 |
|
|
359 C AT T = .4000E+01 Y = .905516E+00 .224044E-04 .944615E-01 |
|
|
360 C S(1) = -.187621E+01 .179197E-03 .187603E+01 |
|
|
361 C S(2) = .296093E-05 -.583104E-09 -.296034E-05 |
|
|
362 C S(3) = -.493267E-09 -.276246E-12 .493544E-09 |
|
|
363 C AT T = .4000E+02 Y = .715848E+00 .918628E-05 .284143E+00 |
|
|
364 C S(1) = -.424730E+01 .459360E-04 .424726E+01 |
|
|
365 C S(2) = .137294E-04 -.235815E-09 -.137291E-04 |
|
|
366 C S(3) = -.228818E-08 -.113803E-12 .228829E-08 |
|
|
367 C AT T = .4000E+03 Y = .450526E+00 .322299E-05 .549471E+00 |
|
|
368 C S(1) = -.595837E+01 .354310E-05 .595836E+01 |
|
|
369 C S(2) = .227380E-04 -.226041E-10 -.227380E-04 |
|
|
370 C S(3) = -.378971E-08 -.499501E-13 .378976E-08 |
|
|
371 C AT T = .4000E+04 Y = .183185E+00 .894131E-06 .816814E+00 |
|
|
372 C S(1) = -.475006E+01 -.599504E-05 .475007E+01 |
|
|
373 C S(2) = .188089E-04 .231330E-10 -.188089E-04 |
|
|
374 C S(3) = -.313478E-08 -.187575E-13 .313480E-08 |
|
|
375 C AT T = .4000E+05 Y = .389733E-01 .162133E-06 .961027E+00 |
|
|
376 C S(1) = -.157477E+01 -.276199E-05 .157477E+01 |
|
|
377 C S(2) = .628668E-05 .110026E-10 -.628670E-05 |
|
|
378 C S(3) = -.104776E-08 -.453588E-14 .104776E-08 |
|
|
379 C AT T = .4000E+06 Y = .493609E-02 .198411E-07 .995064E+00 |
|
|
380 C S(1) = -.236244E+00 -.458262E-06 .236244E+00 |
|
|
381 C S(2) = .944669E-06 .183193E-11 -.944671E-06 |
|
|
382 C S(3) = -.157441E-09 -.635990E-15 .157442E-09 |
|
|
383 C AT T = .4000E+07 Y = .516087E-03 .206540E-08 .999484E+00 |
|
|
384 C S(1) = -.256277E-01 -.509808E-07 .256278E-01 |
|
|
385 C S(2) = .102506E-06 .203905E-12 -.102506E-06 |
|
|
386 C S(3) = -.170825E-10 -.684002E-16 .170826E-10 |
|
|
387 C AT T = .4000E+08 Y = .519314E-04 .207736E-09 .999948E+00 |
|
|
388 C S(1) = -.259316E-02 -.518029E-08 .259316E-02 |
|
|
389 C S(2) = .103726E-07 .207209E-13 -.103726E-07 |
|
|
390 C S(3) = -.172845E-11 -.691450E-17 .172845E-11 |
|
|
391 C AT T = .4000E+09 Y = .544710E-05 .217885E-10 .999995E+00 |
|
|
392 C S(1) = -.271637E-03 -.541849E-09 .271638E-03 |
|
|
393 C S(2) = .108655E-08 .216739E-14 -.108655E-08 |
|
|
394 C S(3) = -.180902E-12 -.723615E-18 .180902E-12 |
|
|
395 C AT T = .4000E+10 Y = .446748E-06 .178699E-11 .100000E+01 |
|
|
396 C S(1) = -.322322E-04 -.842541E-10 .322323E-04 |
|
|
397 C S(2) = .128929E-09 .337016E-15 -.128929E-09 |
|
|
398 C S(3) = -.209715E-13 -.838859E-19 .209715E-13 |
|
|
399 C AT T = .4000E+11 Y = -.363960E-07 -.145584E-12 .100000E+01 |
|
|
400 C S(1) = -.164109E-06 -.429604E-11 .164113E-06 |
|
|
401 C S(2) = .656436E-12 .171842E-16 -.656451E-12 |
|
|
402 C S(3) = -.689361E-15 -.275745E-20 .689363E-15 |
|
|
403 C |
|
|
404 C NO. STEPS = 340 NO. F-S = 412 NO. J-S = 343 NO. DF-S =1023 |
|
|
405 C---------------------------------------------------------------------- |
|
|
406 C FULL DESCRIPTION OF USER INTERFACE TO ODESSA. |
|
|
407 C |
|
|
408 C THE USER INTERFACE TO ODESSA CONSISTS OF THE FOLLOWING PARTS. |
|
|
409 C |
|
|
410 C I. THE CALL SEQUENCE TO SUBROUTINE ODESSA, WHICH IS A DRIVER |
|
|
411 C ROUTINE FOR THE SOLVER. THIS INCLUDES DESCRIPTIONS OF BOTH |
|
|
412 C THE CALL SEQUENCE ARGUMENTS AND OF USER-SUPPLIED ROUTINES. |
|
|
413 C FOLLOWING THESE DESCRIPTIONS IS A DESCRIPTION OF |
|
|
414 C OPTIONAL INPUTS AVAILABLE THROUGH THE CALL SEQUENCE, AND THEN |
|
|
415 C A DESCRIPTION OF OPTIONAL OUTPUTS (IN THE WORK ARRAYS). |
|
|
416 C |
|
|
417 C II. DESCRIPTIONS OF OTHER ROUTINES IN THE ODESSA PACKAGE THAT MAY |
|
|
418 C BE (OPTIONALLY) CALLED BY THE USER. THESE PROVIDE THE ABILITY |
|
|
419 C TO ALTER ERROR MESSAGE HANDLING, SAVE AND RESTORE THE INTERNAL |
|
|
420 C COMMON, AND OBTAIN SPECIFIED DERIVATIVES OF THE SOLUTION Y(T). |
|
|
421 C |
|
|
422 C III. DESCRIPTIONS OF COMMON BLOCKS TO BE DECLARED IN OVERLAY |
|
|
423 C OR SIMILAR ENVIRONMENTS, OR TO BE SAVED WHEN DOING AN INTERRUPT |
|
|
424 C OF THE PROBLEM AND CONTINUED SOLUTION LATER. |
|
|
425 C |
|
|
426 C IV. DESCRIPTION OF TWO SUBROUTINES IN THE ODESSA PACKAGE, EITHER OF |
|
|
427 C WHICH THE USER MAY REPLACE WITH HIS OWN VERSION, IF DESIRED. |
|
|
428 C THESE RELATE TO THE MEASUREMENT OF ERRORS. |
|
|
429 C |
|
|
430 C V. GENERAL REMARKS WHICH HIGHLIGHT DIFFERENCES BETWEEN THE LSODE |
|
|
431 C PACKAGE AND THE ODESSA PACKAGE. |
|
|
432 C---------------------------------------------------------------------- |
|
|
433 C PART I. CALL SEQUENCE. |
|
|
434 C |
|
|
435 C THE CALL SEQUENCE PARAMETERS USED FOR INPUT ONLY ARE.. |
|
|
436 C F, DF, NEQ, PAR, TOUT, ITOL, RTOL, ATOL, ITASK, IOPT, LRW, LIW, |
|
|
437 C JAC, MF, |
|
|
438 C AND THOSE USED FOR BOTH INPUT AND OUTPUT ARE |
|
|
439 C Y, T, ISTATE. |
|
|
440 C THE WORK ARRAYS RWORK AND IWORK ARE ALSO USED FOR CONDITIONAL AND |
|
|
441 C OPTIONAL INPUTS AND OPTIONAL OUTPUTS. (THE TERM OUTPUT HERE REFERS |
|
|
442 C TO THE RETURN FROM SUBROUTINE ODESSA TO THE USER-S CALLING PROGRAM.) |
|
|
443 C |
|
|
444 C THE LEGALITY OF INPUT PARAMETERS WILL BE THOROUGHLY CHECKED ON THE |
|
|
445 C INITIAL CALL FOR THE PROBLEM, BUT NOT CHECKED THEREAFTER UNLESS A |
|
|
446 C CHANGE IN INPUT PARAMETERS IS FLAGGED BY ISTATE = 3 ON INPUT. |
|
|
447 C |
|
|
448 C THE DESCRIPTIONS OF THE CALL ARGUMENTS ARE AS FOLLOWS. |
|
|
449 C |
|
|
450 C F = THE NAME OF THE USER-SUPPLIED SUBROUTINE DEFINING THE |
|
|
451 C ODE MODEL. THIS SYSTEM MUST BE PUT IN THE FIRST-ORDER |
|
|
452 C FORM DY/DT = F(Y,T;P), WHERE F IS A VECTOR-VALUED FUNCTION |
|
|
453 C OF THE SCALAR T AND VECTORS Y, AND PAR. SUBROUTINE F IS TO |
|
|
454 C COMPUTE THE FUNCTION F. IT IS TO HAVE THE FORM.. |
|
|
455 C SUBROUTINE F (NEQ, T, Y, PAR, YDOT) |
|
|
456 C DOUBLE PRECISION T, Y, PAR, YDOT |
|
|
457 C DIMENSION Y(1), PAR(1), YDOT(1) |
|
|
458 C WHERE NEQ, T, Y, AND PAR ARE INPUT, AND YDOT = F(Y,T;P) |
|
|
459 C IS OUTPUT. Y AND YDOT ARE ARRAYS OF LENGTH N (= NEQ(1)). |
|
|
460 C (IN THE DIMENSION STATEMENT ABOVE, 1 IS A DUMMY |
|
|
461 C DIMENSION.. IT CAN BE REPLACED BY ANY VALUE.) |
|
|
462 C F SHOULD NOT ALTER ARRAY Y, OR PAR(1),...,PAR(NPAR). |
|
|
463 C F MUST BE DECLARED EXTERNAL IN THE CALLING PROGRAM. |
|
|
464 C |
|
|
465 C SUBROUTINE F MAY ACCESS USER-DEFINED QUANTITIES IN |
|
|
466 C NEQ(2),... AND PAR(NPAR+1),... IF NEQ IS AN ARRAY |
|
|
467 C (DIMENSIONED IN F) AND PAR HAS LENGTH EXCEEDING NPAR. |
|
|
468 C SEE THE DESCRIPTIONS OF NEQ AND PAR BELOW. |
|
|
469 C |
|
|
470 C DF = THE NAME OF THE USER-SUPPLIED ROUTINE (IDF = 1) TO COMPUTE |
|
|
471 C THE INHOMOGENEITY MATRIX, DF/DP, AS A FUNCTION OF THE SCALAR |
|
|
472 C T, AND THE VECTORS Y, AND PAR. IT IS TO HAVE THE FORM |
|
|
473 C SUBROUTINE DF (NEQ, T, Y, PAR, DFDP, JPAR) |
|
|
474 C DOUBLE PRECISION T, Y, PAR, DFDP |
|
|
475 C DIMENSION Y(1), PAR(1), DFDP(1) |
|
|
476 C GO TO (1,2,...,NPAR) JPAR |
|
|
477 C 1 DFDP(1) = DF(1)/DP(1) |
|
|
478 C . |
|
|
479 C DFDP(I) = DF(I)/DP(1) |
|
|
480 C . |
|
|
481 C DFDP(N) = DF(N)/DP(1) |
|
|
482 C RETURN |
|
|
483 C 2 DFDP(1) = DF(1)/DP(2) |
|
|
484 C . |
|
|
485 C DFDP(I) = DF(I)/DP(2) |
|
|
486 C . |
|
|
487 C DFDP(N) = DF(N)/DP(2) |
|
|
488 C . |
|
|
489 C RETURN |
|
|
490 C . . |
|
|
491 C . . |
|
|
492 C NPAR DFDP(1) = DF(1)/DP(NPAR) |
|
|
493 C . |
|
|
494 C DFDP(I) = DF(I)/DP(NPAR) |
|
|
495 C . |
|
|
496 C DFDP(N) = DF(N)/DP(NPAR) |
|
|
497 C RETURN |
|
|
498 C END |
|
|
499 C WHERE NEQ, T, Y, PAR, AND JPAR ARE INPUT AND THE VECTOR |
|
|
500 C DFDP(*,JPAR) IS TO BE LOADED WITH THE PARTIAL DERIVATIVES |
|
|
501 C DF(Y,T;PAR)/DP(JPAR) ON OUTPUT. ONLY NONZERO ELEMENTS NEED |
|
|
502 C BE LOADED. T, Y, AND PAR HAVE THE SAME MEANING AS IN |
|
|
503 C SUBROUTINE F. (IN THE DIMENSION STATEMENT ABOVE, 1 IS A DUMMY |
|
|
504 C DIMENSION.. IT CAN BE REPLACED BY ANY VALUE). |
|
|
505 C |
|
|
506 C DFDP(*,JPAR) IS PRESET TO ZERO BY THE SOLVER, SO THAT ONLY |
|
|
507 C THE NONZERO ELEMENTS NEED BE LOADED BY DF. SUBROUTINE DF |
|
|
508 C MUST BE DECLARED EXTERNAL IN THE CALLING PROGRAM IF USED. |
|
|
509 C IF IDF = 0 (OR ISOPT = 0), A DUMMY ARGUMENT CAN BE USED. |
|
|
510 C |
|
|
511 C SUBROUTINE DF MAY ACCESS USER-DEFINED QUANTITIES IN |
|
|
512 C NEQ(2),... AND PAR(NPAR+1),... IF NEQ IS AN ARRAY |
|
|
513 C (DIMENSIONED IN DF) AND PAR HAS A LENGTH EXCEEDING NPAR. |
|
|
514 C SEE THE DESCRIPTIONS OF NEQ AND PAR (BELOW). |
|
|
515 C |
|
|
516 C NEQ = THE SIZE OF THE ODE SYSTEM (NUMBER OF FIRST ORDER ORDINARY |
|
|
517 C DIFFERENTIAL EQUATIONS (N) IN THE MODEL). USED ONLY FOR |
|
|
518 C INPUT. NEQ MAY NOT BE CHANGED DURING THE PROBLEM. |
|
|
519 C |
|
|
520 C FOR ISOPT = 0, NEQ IS NORMALLY A SCALAR. HOWEVER, NEQ MAY |
|
|
521 C BE AN ARRAY, WITH NEQ(1) SET TO THE SYSTEM SIZE (N), IN WHICH |
|
|
522 C CASE THE ODESSA PACKAGE ACCESSES ONLY NEQ(1). HOWEVER, |
|
|
523 C THIS PARAMETER IS PASSED AS THE NEQ ARGUMENT IN ALL CALLS |
|
|
524 C TO F, DF, AND JAC. HENCE, IF IT IS AN ARRAY, LOCATIONS |
|
|
525 C NEQ(2),... MAY BE USED TO STORE OTHER INTEGER DATA AND PASS |
|
|
526 C IT TO F, DF, AND/OR JAC. FOR ISOPT = 1, NPAR MUST BE LOADED |
|
|
527 C INTO NEQ(2), AND IS NOT ALLOWED TO CHANGE DURING THE PROBLEM. |
|
|
528 C IN THESE CASES, SUBROUTINES F, DF, AND/OR JAC MUST INCLUDE |
|
|
529 C NEQ IN A DIMENSION STATEMENT. |
|
|
530 C |
|
|
531 C Y = A REAL ARRAY FOR THE VECTOR OF DEPENDENT VARIABLES, OF |
|
|
532 C DIMENSION (N) BY (NPAR+1). USED FOR BOTH INPUT AND |
|
|
533 C OUTPUT ON THE FIRST CALL (ISTATE = 1), AND ONLY FOR |
|
|
534 C OUTPUT ON OTHER CALLS. ON THE FIRST CALL, Y MUST CONTAIN |
|
|
535 C THE VECTORS OF INITIAL VALUES. ON OUTPUT, Y CONTAINS THE |
|
|
536 C COMPUTED SOLUTION VECTORS, EVALUATED AT T. |
|
|
537 C |
|
|
538 C PAR = A REAL ARRAY FOR THE VECTOR OF CONSTANT MODEL PARAMETERS |
|
|
539 C OF INTEREST IN THE SENSITIVITY ANALYSIS, OF LENGTH NPAR |
|
|
540 C OR MORE. PAR IS PASSED AS AN ARGUMENT IN ALL CALLS TO F, |
|
|
541 C DF, AND JAC. HENCE LOCATIONS PAR(NPAR+1),... MAY BE USED |
|
|
542 C TO STORE OTHER REAL DATA AND PASS IT TO F, DF, AND/OR JAC. |
|
|
543 C LOCATIONS PAR(1),...,PAR(NPAR) ARE USED AS INPUT ONLY, |
|
|
544 C AND MUST NOT BE CHANGED DURING THE PROBLEM. |
|
|
545 C |
|
|
546 C T = THE INDEPENDENT VARIABLE. ON INPUT, T IS USED ONLY ON THE |
|
|
547 C FIRST CALL, AS THE INITIAL POINT OF THE INTEGRATION. |
|
|
548 C ON OUTPUT, AFTER EACH CALL, T IS THE VALUE AT WHICH A |
|
|
549 C COMPUTED SOLUTION Y IS EVALUATED (USUALLY THE SAME AS TOUT). |
|
|
550 C ON AN ERROR RETURN, T IS THE FARTHEST POINT REACHED. |
|
|
551 C |
|
|
552 C TOUT = THE NEXT VALUE OF T AT WHICH A COMPUTED SOLUTION IS DESIRED. |
|
|
553 C USED ONLY FOR INPUT. |
|
|
554 C |
|
|
555 C WHEN STARTING THE PROBLEM (ISTATE = 1), TOUT MAY BE EQUAL |
|
|
556 C TO T FOR ONE CALL, THEN SHOULD .NE. T FOR THE NEXT CALL. |
|
|
557 C FOR THE INITIAL T, AN INPUT VALUE OF TOUT .NE. T IS USED |
|
|
558 C IN ORDER TO DETERMINE THE DIRECTION OF THE INTEGRATION |
|
|
559 C (I.E. THE ALGEBRAIC SIGN OF THE STEP SIZES) AND THE ROUGH |
|
|
560 C SCALE OF THE PROBLEM. INTEGRATION IN EITHER DIRECTION |
|
|
561 C (FORWARD OR BACKWARD IN T) IS PERMITTED. |
|
|
562 C |
|
|
563 C IF ITASK = 2 OR 5 (ONE-STEP MODES), TOUT IS IGNORED AFTER |
|
|
564 C THE FIRST CALL (I.E. THE FIRST CALL WITH TOUT .NE. T). |
|
|
565 C OTHERWISE, TOUT IS REQUIRED ON EVERY CALL. |
|
|
566 C |
|
|
567 C IF ITASK = 1, 3, OR 4, THE VALUES OF TOUT NEED NOT BE |
|
|
568 C MONOTONE, BUT A VALUE OF TOUT WHICH BACKS UP IS LIMITED |
|
|
569 C TO THE CURRENT INTERNAL T INTERVAL, WHOSE ENDPOINTS ARE |
|
|
570 C TCUR - HU AND TCUR (SEE OPTIONAL OUTPUTS, BELOW, FOR |
|
|
571 C TCUR AND HU). |
|
|
572 C |
|
|
573 C ITOL = AN INDICATOR FOR THE TYPE OF ERROR CONTROL. SEE |
|
|
574 C DESCRIPTION BELOW UNDER ATOL. USED ONLY FOR INPUT. |
|
|
575 C |
|
|
576 C RTOL = A RELATIVE ERROR TOLERANCE PARAMETER, EITHER A SCALAR OR |
|
|
577 C AN ARRAY OF SPACE (N) BY (NPAR+1). SEE DESCRIPTION BELOW |
|
|
578 C UNDER ATOL. INPUT ONLY. |
|
|
579 C |
|
|
580 C ATOL = AN ABSOLUTE ERROR TOLERANCE PARAMETER, EITHER A SCALAR OR |
|
|
581 C AN ARRAY OF SPACE (N) BY (NPAR+1). INPUT ONLY. |
|
|
582 C |
|
|
583 C THE INPUT PARAMETERS ITOL, RTOL, AND ATOL DETERMINE |
|
|
584 C THE ERROR CONTROL PERFORMED BY THE SOLVER. THE SOLVER WILL |
|
|
585 C CONTROL THE VECTOR E = (E(I,J)) OF ESTIMATED LOCAL ERRORS |
|
|
586 C IN Y, ACCORDING TO AN INEQUALITY OF THE FORM |
|
|
587 C RMS-NORM OF ( E(I,J)/EWT(I,J) ) .LE. 1, |
|
|
588 C WHERE EWT(I,J) = RTOL(I,J)*ABS(Y(I,J)) + ATOL(I,J), |
|
|
589 C AND THE RMS-NORM (ROOT-MEAN-SQUARE NORM) HERE IS |
|
|
590 C RMS-NORM(V) = SQRT ( (1/N) * SUM (V(I,J)**2) ); I =1,...,N. |
|
|
591 C HERE EWT = (EWT(I,J)) IS A VECTOR OF WEIGHTS WHICH MUST |
|
|
592 C ALWAYS BE POSITIVE, AND THE VALUES OF RTOL AND ATOL SHOULD |
|
|
593 C ALL BE NON-NEGATIVE. THE FOLLOWING TABLE GIVES THE TYPES |
|
|
594 C (SCALAR/ARRAY) OF RTOL AND ATOL, AND THE CORRESPONDING FORM |
|
|
595 C OF EWT(I,J). |
|
|
596 C |
|
|
597 C ITOL RTOL ATOL EWT(I,J) |
|
|
598 C 1 SCALAR SCALAR RTOL*ABS(Y(I,J)) + ATOL |
|
|
599 C 2 SCALAR ARRAY RTOL*ABS(Y(I,J)) + ATOL(I,J) |
|
|
600 C 3 ARRAY SCALAR RTOL(I,J)*ABS(Y(I,J)) + ATOL |
|
|
601 C 4 ARRAY ARRAY RTOL(I,J)*ABS(Y(I,J)) + ATOL(I,J) |
|
|
602 C |
|
|
603 C WHEN EITHER OF THESE PARAMETERS IS A SCALAR, IT NEED NOT |
|
|
604 C BE DIMENSIONED IN THE USER-S CALLING PROGRAM. |
|
|
605 C |
|
|
606 C THE TOTAL NUMBER OF ERROR TEST FAILURES DUE TO THE SENSITIVITY |
|
|
607 C ANALYSIS, AND WHICH REQUIRE AN INTEGRATION STEP TO BE |
|
|
608 C REPEATED, ARE ACCUMULATED IN THE LAST NPAR+1 LOCATIONS OF THE |
|
|
609 C INTEGER WORK ARRAY IWORK (SEE OPTIONAL OUTPUTS BELOW). |
|
|
610 C THIS INFORMATION MAY BE OF VALUE IN DETERMINING APPROPRIATE |
|
|
611 C ERROR TOLERANCES TO BE APPLIED TO THE SENSITIVITY FUNCTIONS. |
|
|
612 C |
|
|
613 C IF NONE OF THE ABOVE CHOICES (WITH ITOL, RTOL, AND ATOL |
|
|
614 C FIXED THROUGHOUT THE PROBLEM) IS SUITABLE, MORE GENERAL |
|
|
615 C ERROR CONTROLS CAN BE OBTAINED BY SUBSTITUTING |
|
|
616 C USER-SUPPLIED ROUTINES FOR THE SETTING OF EWT AND/OR FOR |
|
|
617 C THE NORM CALCULATION. SEE PART IV BELOW. |
|
|
618 C |
|
|
619 C IF GLOBAL ERRORS ARE TO BE ESTIMATED BY MAKING A REPEATED |
|
|
620 C RUN ON THE SAME PROBLEM WITH SMALLER TOLERANCES, THEN ALL |
|
|
621 C COMPONENTS OF RTOL AND ATOL (I.E. OF EWT) SHOULD BE SCALED |
|
|
622 C DOWN UNIFORMLY. |
|
|
623 C |
|
|
624 C ITASK = AN INDEX SPECIFYING THE TASK TO BE PERFORMED. |
|
|
625 C INPUT ONLY. ITASK HAS THE FOLLOWING VALUES AND MEANINGS. |
|
|
626 C 1 MEANS NORMAL COMPUTATION OF OUTPUT VALUES OF Y(T) AT |
|
|
627 C T = TOUT (BY OVERSHOOTING AND INTERPOLATING). |
|
|
628 C 2 MEANS TAKE ONE STEP ONLY AND RETURN. |
|
|
629 C 3 MEANS STOP AT THE FIRST INTERNAL MESH POINT AT OR |
|
|
630 C BEYOND T = TOUT AND RETURN. |
|
|
631 C 4 MEANS NORMAL COMPUTATION OF OUTPUT VALUES OF Y(T) AT |
|
|
632 C T = TOUT BUT WITHOUT OVERSHOOTING T = TCRIT. |
|
|
633 C TCRIT MUST BE INPUT AS RWORK(1). TCRIT MAY BE EQUAL TO |
|
|
634 C OR BEYOND TOUT, BUT NOT BEHIND IT IN THE DIRECTION OF |
|
|
635 C INTEGRATION. THIS OPTION IS USEFUL IF THE PROBLEM |
|
|
636 C HAS A SINGULARITY AT OR BEYOND T = TCRIT. |
|
|
637 C 5 MEANS TAKE ONE STEP, WITHOUT PASSING TCRIT, AND RETURN. |
|
|
638 C TCRIT MUST BE INPUT AS RWORK(1). |
|
|
639 C |
|
|
640 C NOTE.. IF ITASK = 4 OR 5 AND THE SOLVER REACHES TCRIT |
|
|
641 C (WITHIN ROUNDOFF), IT WILL RETURN T = TCRIT (EXACTLY) TO |
|
|
642 C INDICATE THIS (UNLESS ITASK = 4 AND TOUT COMES BEFORE TCRIT, |
|
|
643 C IN WHICH CASE ANSWERS AT T = TOUT ARE RETURNED FIRST). |
|
|
644 C |
|
|
645 C ISTATE = AN INDEX USED FOR INPUT AND OUTPUT TO SPECIFY THE |
|
|
646 C THE STATE OF THE CALCULATION. |
|
|
647 C |
|
|
648 C ON INPUT, THE VALUES OF ISTATE ARE AS FOLLOWS. |
|
|
649 C 1 MEANS THIS IS THE FIRST CALL FOR THE PROBLEM |
|
|
650 C (INITIALIZATIONS WILL BE DONE). SEE NOTE BELOW. |
|
|
651 C 2 MEANS THIS IS NOT THE FIRST CALL, AND THE CALCULATION |
|
|
652 C IS TO CONTINUE NORMALLY, WITH NO CHANGE IN ANY INPUT |
|
|
653 C PARAMETERS EXCEPT POSSIBLY TOUT AND ITASK. |
|
|
654 C (IF ITOL, RTOL, AND/OR ATOL ARE CHANGED BETWEEN CALLS |
|
|
655 C WITH ISTATE = 2, THE NEW VALUES WILL BE USED BUT NOT |
|
|
656 C TESTED FOR LEGALITY.) |
|
|
657 C 3 MEANS THIS IS NOT THE FIRST CALL, AND THE |
|
|
658 C CALCULATION IS TO CONTINUE NORMALLY, BUT WITH |
|
|
659 C A CHANGE IN INPUT PARAMETERS OTHER THAN |
|
|
660 C TOUT AND ITASK. CHANGES ARE ALLOWED IN |
|
|
661 C ITOL, RTOL, ATOL, IOPT, LRW, LIW, MF, ML, MU, |
|
|
662 C AND ANY OF THE OPTIONAL INPUTS EXCEPT H0. |
|
|
663 C (SEE IWORK DESCRIPTION FOR ML AND MU.) |
|
|
664 C NOTE.. A PRELIMINARY CALL WITH TOUT = T IS NOT COUNTED |
|
|
665 C AS A FIRST CALL HERE, AS NO INITIALIZATION OR CHECKING OF |
|
|
666 C INPUT IS DONE. (SUCH A CALL IS SOMETIMES USEFUL FOR THE |
|
|
667 C PURPOSE OF OUTPUTTING THE INITIAL CONDITIONS.) |
|
|
668 C THUS THE FIRST CALL FOR WHICH TOUT .NE. T REQUIRES |
|
|
669 C ISTATE = 1 ON INPUT. |
|
|
670 C |
|
|
671 C ON OUTPUT, ISTATE HAS THE FOLLOWING VALUES AND MEANINGS. |
|
|
672 C 1 MEANS NOTHING WAS DONE, AS TOUT WAS EQUAL TO T WITH |
|
|
673 C ISTATE = 1 ON INPUT. (HOWEVER, AN INTERNAL COUNTER WAS |
|
|
674 C SET TO DETECT AND PREVENT REPEATED CALLS OF THIS TYPE.) |
|
|
675 C 2 MEANS THE INTEGRATION WAS PERFORMED SUCCESSFULLY. |
|
|
676 C -1 MEANS AN EXCESSIVE AMOUNT OF WORK (MORE THAN MXSTEP |
|
|
677 C STEPS) WAS DONE ON THIS CALL, BEFORE COMPLETING THE |
|
|
678 C REQUESTED TASK, BUT THE INTEGRATION WAS OTHERWISE |
|
|
679 C SUCCESSFUL AS FAR AS T. (MXSTEP IS AN OPTIONAL INPUT |
|
|
680 C AND IS NORMALLY 500.) TO CONTINUE, THE USER MAY |
|
|
681 C SIMPLY RESET ISTATE TO A VALUE .GT. 1 AND CALL AGAIN |
|
|
682 C (THE EXCESS WORK STEP COUNTER WILL BE RESET TO 0). |
|
|
683 C IN ADDITION, THE USER MAY INCREASE MXSTEP TO AVOID |
|
|
684 C THIS ERROR RETURN (SEE BELOW ON OPTIONAL INPUTS). |
|
|
685 C -2 MEANS TOO MUCH ACCURACY WAS REQUESTED FOR THE PRECISION |
|
|
686 C OF THE MACHINE BEING USED. THIS WAS DETECTED BEFORE |
|
|
687 C COMPLETING THE REQUESTED TASK, BUT THE INTEGRATION |
|
|
688 C WAS SUCCESSFUL AS FAR AS T. TO CONTINUE, THE TOLERANCE |
|
|
689 C PARAMETERS MUST BE RESET, AND ISTATE MUST BE SET |
|
|
690 C TO 3. THE OPTIONAL OUTPUT TOLSF MAY BE USED FOR THIS |
|
|
691 C PURPOSE. (NOTE.. IF THIS CONDITION IS DETECTED BEFORE |
|
|
692 C TAKING ANY STEPS, THEN AN ILLEGAL INPUT RETURN |
|
|
693 C (ISTATE = -3) OCCURS INSTEAD.) |
|
|
694 C -3 MEANS ILLEGAL INPUT WAS DETECTED, BEFORE TAKING ANY |
|
|
695 C INTEGRATION STEPS. SEE WRITTEN MESSAGE FOR DETAILS. |
|
|
696 C NOTE.. IF THE SOLVER DETECTS AN INFINITE LOOP OF CALLS |
|
|
697 C TO THE SOLVER WITH ILLEGAL INPUT, IT WILL CAUSE |
|
|
698 C THE RUN TO STOP. |
|
|
699 C -4 MEANS THERE WERE REPEATED ERROR TEST FAILURES ON |
|
|
700 C ONE ATTEMPTED STEP, BEFORE COMPLETING THE REQUESTED |
|
|
701 C TASK, BUT THE INTEGRATION WAS SUCCESSFUL AS FAR AS T. |
|
|
702 C THE PROBLEM MAY HAVE A SINGULARITY, OR THE INPUT |
|
|
703 C MAY BE INAPPROPRIATE. |
|
|
704 C -5 MEANS THERE WERE REPEATED CONVERGENCE TEST FAILURES ON |
|
|
705 C ONE ATTEMPTED STEP, BEFORE COMPLETING THE REQUESTED |
|
|
706 C TASK, BUT THE INTEGRATION WAS SUCCESSFUL AS FAR AS T. |
|
|
707 C THIS MAY BE CAUSED BY AN INACCURATE JACOBIAN MATRIX, |
|
|
708 C IF ONE IS BEING USED. |
|
|
709 C -6 MEANS EWT(I,J) BECAME ZERO FOR SOME I,J DURING THE |
|
|
710 C INTEGRATION. PURE RELATIVE ERROR CONTROL (ATOL(I,J)=0.0) |
|
|
711 C WAS REQUESTED ON A VARIABLE WHICH HAS NOW VANISHED. |
|
|
712 C THE INTEGRATION WAS SUCCESSFUL AS FAR AS T. |
|
|
713 C |
|
|
714 C NOTE.. SINCE THE NORMAL OUTPUT VALUE OF ISTATE IS 2, |
|
|
715 C IT DOES NOT NEED TO BE RESET FOR NORMAL CONTINUATION. |
|
|
716 C ALSO, SINCE A NEGATIVE INPUT VALUE OF ISTATE WILL BE |
|
|
717 C REGARDED AS ILLEGAL, A NEGATIVE OUTPUT VALUE REQUIRES THE |
|
|
718 C USER TO CHANGE IT, AND POSSIBLY OTHER INPUTS, BEFORE |
|
|
719 C CALLING THE SOLVER AGAIN. |
|
|
720 C |
|
|
721 C IOPT = AN INTEGER ARRAY FLAG TO SPECIFY WHETHER OR NOT ANY OPTIONAL |
|
|
722 C INPUTS ARE BEING USED ON THIS CALL. INPUT ONLY. |
|
|
723 C THE OPTIONAL INPUTS ARE LISTED SEPARATELY BELOW. |
|
|
724 C IOPT(1) = 0 MEANS NO OPTIONAL INPUTS FOR THE SOLVER WILL BE |
|
|
725 C USED. DEFAULT VALUES WILL BE USED IN ALL CASES. |
|
|
726 C = 1 MEANS ONE OR MORE OPTIONAL INPUTS FOR THE |
|
|
727 C SOLVER ARE BEING USED. |
|
|
728 C NOTE : IOPT(1) IS INDEPENDENT OF ISOPT AND IDF. |
|
|
729 C IOPT(2) = 0 MEANS NO SENSITIVITY ANALYSIS WILL BE PERFORMED. |
|
|
730 C = 1 MEANS A SENSITIVITY ANALYSIS WILL BE PERFORMED. |
|
|
731 C NOTE : IOPT(2) IS RENAMED TO ISOPT IN ODESSA. |
|
|
732 C = 0 MEANS DF/DP WILL BE CALCULATED BY FINITE |
|
|
733 C DIFFERENCE WITHIN ODESSA. |
|
|
734 C IOPT(3) = 1 MEANS DF/DP WILL BE CALCULATED BY A USER-SUPPLIED |
|
|
735 C ROUTINE. |
|
|
736 C NOTE : IOPT(3) IS RENAMED TO IDF IN ODESSA. |
|
|
737 C IF IDF = 1, THE USER MUST SUPPLY A |
|
|
738 C SUBROUTINE DF (THE NAME IS ARBITRARY) AS |
|
|
739 C DESCRIBED BELOW UNDER DF. FOR IDF = 0, |
|
|
740 C A DUMMY ARGUMENT CAN BE USED. |
|
|
741 C |
|
|
742 C RWORK = A REAL WORKING ARRAY (DOUBLE PRECISION). |
|
|
743 C FOR ISOPT = 0, THE LENGTH OF RWORK MUST BE AT LEAST.. |
|
|
744 C 20 + NYH*(MAXORD + 1) + 3*NEQ + LWM |
|
|
745 C FOR ISOPT = 1, THE LENGTH OF RWORK MUST BE AT LEAST.. |
|
|
746 C 20 + NYH*(MAXORD + 1) + 2*NYH + LWM + N |
|
|
747 C WHERE.. |
|
|
748 C NYH = THE TOTAL NUMBER OF DEPENDENT VARIABLES; |
|
|
749 C (= N IF ISOPT = 0, AND N*(NPAR+1) IF ISOPT = 1). |
|
|
750 C MAXORD = 12 (IF METH = 1) OR 5 (IF METH = 2) (UNLESS A |
|
|
751 C SMALLER VALUE IS GIVEN AS AN OPTIONAL INPUT), |
|
|
752 C LWM = 0 IF MITER = 0, |
|
|
753 C LWM = N**2 + 2 IF MITER IS 1 OR 2, |
|
|
754 C LWM = N + 2 IF MITER = 3, AND |
|
|
755 C LWM = (2*ML+MU+1)*N + 2 IF MITER IS 4 OR 5. |
|
|
756 C (SEE THE MF DESCRIPTION FOR METH AND MITER.) |
|
|
757 C |
|
|
758 C THE FIRST 20 WORDS OF RWORK ARE RESERVED FOR CONDITIONAL |
|
|
759 C AND OPTIONAL INPUTS AND OPTIONAL OUTPUTS. |
|
|
760 C |
|
|
761 C THE FOLLOWING WORD IN RWORK IS A CONDITIONAL INPUT.. |
|
|
762 C RWORK(1) = TCRIT = CRITICAL VALUE OF T WHICH THE SOLVER |
|
|
763 C IS NOT TO OVERSHOOT. REQUIRED IF ITASK IS |
|
|
764 C 4 OR 5, AND IGNORED OTHERWISE. (SEE ITASK.) |
|
|
765 C |
|
|
766 C LRW = THE LENGTH OF THE ARRAY RWORK, AS DECLARED BY THE USER. |
|
|
767 C (THIS WILL BE CHECKED BY THE SOLVER.) |
|
|
768 C |
|
|
769 C IWORK = AN INTEGER WORK ARRAY. THE LENGTH MUST BE AT LEAST.. |
|
|
770 C 20 IF MITER = 0 OR 3 (MF = 10, 13, 20, 23), OR |
|
|
771 C 20 + N OTHERWISE (MF = 11, 12, 14, 15, 21, 22, 24, 25). |
|
|
772 C FOR ISOPT = 0, OR.. |
|
|
773 C 21 + N + NPAR |
|
|
774 C FOR ISOPT = 1. |
|
|
775 C THE FIRST FEW WORDS OF IWORK ARE USED FOR CONDITIONAL AND |
|
|
776 C OPTIONAL INPUTS AND OPTIONAL OUTPUTS. |
|
|
777 C |
|
|
778 C THE FOLLOWING 2 WORDS IN IWORK ARE CONDITIONAL INPUTS.. |
|
|
779 C IWORK(1) = ML THESE ARE THE LOWER AND UPPER |
|
|
780 C IWORK(2) = MU HALF-BANDWIDTHS, RESPECTIVELY, OF THE |
|
|
781 C BANDED JACOBIAN, EXCLUDING THE MAIN DIAGONAL. |
|
|
782 C THE BAND IS DEFINED BY THE MATRIX LOCATIONS |
|
|
783 C (I,J) WITH I-ML .LE. J .LE. I+MU. ML AND MU |
|
|
784 C MUST SATISFY 0 .LE. ML,MU .LE. NEQ-1. |
|
|
785 C THESE ARE REQUIRED IF MITER IS 4 OR 5, AND |
|
|
786 C IGNORED OTHERWISE. ML AND MU MAY IN FACT BE |
|
|
787 C THE BAND PARAMETERS FOR A MATRIX TO WHICH |
|
|
788 C DF/DY IS ONLY APPROXIMATELY EQUAL. |
|
|
789 * |
|
|
790 C |
|
|
791 C LIW = THE LENGTH OF THE ARRAY IWORK, AS DECLARED BY THE USER. |
|
|
792 C (THIS WILL BE CHECKED BY THE SOLVER.) |
|
|
793 C |
|
|
794 C NOTE.. THE WORK ARRAYS MUST NOT BE ALTERED BETWEEN CALLS TO ODESSA |
|
|
795 C FOR THE SAME PROBLEM, EXCEPT POSSIBLY FOR THE CONDITIONAL AND |
|
|
796 C OPTIONAL INPUTS, AND EXCEPT FOR THE LAST 2*NYH + N WORDS OF RWORK. |
|
|
797 C THE LATTER SPACE IS USED FOR INTERNAL SCRATCH SPACE, AND SO IS |
|
|
798 C AVAILABLE FOR USE BY THE USER OUTSIDE ODESSA BETWEEN CALLS, IF |
|
|
799 C DESIRED (BUT NOT FOR USE BY F, DF, OR JAC). |
|
|
800 C |
|
|
801 C JAC = THE NAME OF THE USER-SUPPLIED ROUTINE (MITER = 1 OR 4) TO |
|
|
802 C COMPUTE THE JACOBIAN MATRIX, DF/DY, AS A FUNCTION OF THE |
|
|
803 C SCALAR T AND THE VECTORS Y, AND PAR. IT IS TO HAVE THE FORM |
|
|
804 C SUBROUTINE JAC (NEQ, T, Y, PAR, ML, MU, PD, NROWPD) |
|
|
805 C DOUBLE PRECISION T, Y, PAR, PD |
|
|
806 C DIMENSION Y(1), PAR(1), PD(NROWPD,1) |
|
|
807 C WHERE NEQ, T, Y, PAR, ML, MU, AND NROWPD ARE INPUT AND THE |
|
|
808 C ARRAY PD IS TO BE LOADED WITH PARTIAL DERIVATIVES (ELEMENTS |
|
|
809 C OF THE JACOBIAN MATRIX) ON OUTPUT. PD MUST BE GIVEN A FIRST |
|
|
810 C DIMENSION OF NROWPD. T, Y, AND PAR HAVE THE SAME MEANING AS |
|
|
811 C IN SUBROUTINE F. (IN THE DIMENSION STATEMENT ABOVE, 1 IS A |
|
|
812 C DUMMY DIMENSION.. IT CAN BE REPLACED BY ANY VALUE.) |
|
|
813 C IN THE FULL MATRIX CASE (MITER = 1), ML AND MU ARE |
|
|
814 C IGNORED, AND THE JACOBIAN IS TO BE LOADED INTO PD IN |
|
|
815 C COLUMNWISE MANNER, WITH DF(I)/DY(J) LOADED INTO PD(I,J). |
|
|
816 C IN THE BAND MATRIX CASE (MITER = 4), THE ELEMENTS |
|
|
817 C WITHIN THE BAND ARE TO BE LOADED INTO PD IN COLUMNWISE |
|
|
818 C MANNER, WITH DIAGONAL LINES OF DF/DY LOADED INTO THE ROWS |
|
|
819 C OF PD. THUS DF(I)/DY(J) IS TO BE LOADED INTO PD(I-J+MU+1,J). |
|
|
820 C ML AND MU ARE THE HALF-BANDWIDTH PARAMETERS (SEE IWORK). |
|
|
821 C THE LOCATIONS IN PD IN THE TWO TRIANGULAR AREAS WHICH |
|
|
822 C CORRESPOND TO NONEXISTENT MATRIX ELEMENTS CAN BE IGNORED |
|
|
823 C OR LOADED ARBITRARILY, AS THEY ARE OVERWRITTEN BY ODESSA. |
|
|
824 C PD IS PRESET TO ZERO BY THE SOLVER, SO THAT ONLY THE |
|
|
825 C NONZERO ELEMENTS NEED BE LOADED BY JAC. EACH CALL TO JAC IS |
|
|
826 C PRECEDED BY A CALL TO F WITH THE SAME ARGUMENTS NEQ, T, Y, |
|
|
827 C AND PAR. THUS TO GAIN SOME EFFICIENCY, INTERMEDIATE |
|
|
828 C QUANTITIES SHARED BY BOTH CALCULATIONS MAY BE SAVED IN A |
|
|
829 C USER COMMON BLOCK BY F AND NOT RECOMPUTED BY JAC, IF |
|
|
830 C DESIRED. ALSO, JAC MAY ALTER THE Y ARRAY, IF DESIRED. |
|
|
831 C JAC MUST BE DECLARED EXTERNAL IN THE CALLING PROGRAM. |
|
|
832 C SUBROUTINE JAC MAY ACCESS USER-DEFINED QUANTITIES IN |
|
|
833 C NEQ(2),... AND PAR(NPAR+1),.... SEE THE DESCRIPTIONS OF |
|
|
834 C NEQ (ABOVE) AND PAR (BELOW). |
|
|
835 C |
|
|
836 C MF = THE METHOD FLAG. USED ONLY FOR INPUT. THE LEGAL VALUES OF |
|
|
837 C MF ARE 10, 11, 12, 13, 14, 15, 20, 21, 22, 23, 24, AND 25. |
|
|
838 C MF HAS DECIMAL DIGITS METH AND MITER.. MF = 10*METH + MITER. |
|
|
839 C METH INDICATES THE BASIC LINEAR MULTISTEP METHOD.. |
|
|
840 C METH = 1 MEANS THE IMPLICIT ADAMS METHOD. |
|
|
841 * |
|
|
842 C METH = 2 MEANS THE METHOD BASED ON BACKWARD |
|
|
843 C DIFFERENTIATION FORMULAS (BDF-S). |
|
|
844 C MITER INDICATES THE CORRECTOR ITERATION METHOD.. |
|
|
845 C MITER = 0 MEANS FUNCTIONAL ITERATION (NO JACOBIAN MATRIX |
|
|
846 C IS INVOLVED). |
|
|
847 C MITER = 1 MEANS CHORD ITERATION WITH A USER-SUPPLIED |
|
|
848 C FULL (NEQ BY NEQ) JACOBIAN. |
|
|
849 C MITER = 2 MEANS CHORD ITERATION WITH AN INTERNALLY |
|
|
850 C GENERATED (DIFFERENCE QUOTIENT) FULL JACOBIAN |
|
|
851 C (USING NEQ EXTRA CALLS TO F PER DF/DY VALUE). |
|
|
852 C MITER = 3 MEANS CHORD ITERATION WITH AN INTERNALLY |
|
|
853 C GENERATED DIAGONAL JACOBIAN APPROXIMATION. |
|
|
854 C (USING 1 EXTRA CALL TO F PER DF/DY EVALUATION). |
|
|
855 C MITER = 4 MEANS CHORD ITERATION WITH A USER-SUPPLIED |
|
|
856 C BANDED JACOBIAN. |
|
|
857 C MITER = 5 MEANS CHORD ITERATION WITH AN INTERNALLY |
|
|
858 C GENERATED BANDED JACOBIAN (USING ML+MU+1 EXTRA |
|
|
859 C CALLS TO F PER DF/DY EVALUATION). |
|
|
860 C IF MITER = 1 OR 4, THE USER MUST SUPPLY A SUBROUTINE JAC |
|
|
861 C (THE NAME IS ARBITRARY) AS DESCRIBED ABOVE UNDER JAC. |
|
|
862 C FOR OTHER VALUES OF MITER, A DUMMY ARGUMENT CAN BE USED. |
|
|
863 C |
|
|
864 C IF A SENSITIVITY ANLYSIS IS DESIRED (ISOPT = 1), MITER = 0 |
|
|
865 C AND 3 ARE DISALLOWED. IN THESE CASES, THE USER IS RECOMMENDED |
|
|
866 C TO SUPPLY AN ANALYTICAL JACOBIAN (MITER = 1 OR 4) AND AN |
|
|
867 C ANALYTICAL INHOMOGENEITY MATRIX (IDF = 1). |
|
|
868 C---------------------------------------------------------------------- |
|
|
869 C OPTIONAL INPUTS. |
|
|
870 C |
|
|
871 C THE FOLLOWING IS A LIST OF THE OPTIONAL INPUTS PROVIDED FOR IN THE |
|
|
872 C CALL SEQUENCE. (SEE ALSO PART II.) FOR EACH SUCH INPUT VARIABLE, |
|
|
873 C THIS TABLE LISTS ITS NAME AS USED IN THIS DOCUMENTATION, ITS |
|
|
874 C LOCATION IN THE CALL SEQUENCE, ITS MEANING, AND THE DEFAULT VALUE. |
|
|
875 C THE USE OF ANY OF THESE INPUTS REQUIRES IOPT(1) = 1, AND IN THAT |
|
|
876 C CASE ALL OF THESE INPUTS ARE EXAMINED. A VALUE OF ZERO FOR ANY |
|
|
877 C OF THESE OPTIONAL INPUTS WILL CAUSE THE DEFAULT VALUE TO BE USED. |
|
|
878 C THUS TO USE A SUBSET OF THE OPTIONAL INPUTS, SIMPLY PRELOAD |
|
|
879 C LOCATIONS 5 TO 10 IN RWORK AND IWORK TO 0.0 AND 0 RESPECTIVELY, AND |
|
|
880 C THEN SET THOSE OF INTEREST TO NONZERO VALUES. |
|
|
881 C |
|
|
882 C NAME LOCATION MEANING AND DEFAULT VALUE |
|
|
883 C |
|
|
884 C H0 RWORK(5) THE STEP SIZE TO BE ATTEMPTED ON THE FIRST STEP. |
|
|
885 C THE DEFAULT VALUE IS DETERMINED BY THE SOLVER. |
|
|
886 C |
|
|
887 C HMAX RWORK(6) THE MAXIMUM ABSOLUTE STEP SIZE ALLOWED. |
|
|
888 C THE DEFAULT VALUE IS INFINITE. |
|
|
889 C |
|
|
890 C HMIN RWORK(7) THE MINIMUM ABSOLUTE STEP SIZE ALLOWED. |
|
|
891 C THE DEFAULT VALUE IS 0. (THIS LOWER BOUND IS NOT |
|
|
892 C ENFORCED ON THE FINAL STEP BEFORE REACHING TCRIT |
|
|
893 C WHEN ITASK = 4 OR 5.) |
|
|
894 C |
|
|
895 C MAXORD IWORK(5) THE MAXIMUM ORDER TO BE ALLOWED. THE DEFAULT |
|
|
896 C VALUE IS 12 IF METH = 1, AND 5 IF METH = 2. |
|
|
897 C IF MAXORD EXCEEDS THE DEFAULT VALUE, IT WILL |
|
|
898 C BE REDUCED TO THE DEFAULT VALUE. |
|
|
899 C IF MAXORD IS CHANGED DURING THE PROBLEM, IT MAY |
|
|
900 C CAUSE THE CURRENT ORDER TO BE REDUCED. |
|
|
901 C |
|
|
902 C MXSTEP IWORK(6) MAXIMUM NUMBER OF (INTERNALLY DEFINED) STEPS |
|
|
903 C ALLOWED DURING ONE CALL TO THE SOLVER. |
|
|
904 C THE DEFAULT VALUE IS 500. |
|
|
905 C |
|
|
906 C MXHNIL IWORK(7) MAXIMUM NUMBER OF MESSAGES PRINTED (PER PROBLEM) |
|
|
907 C WARNING THAT T + H = T ON A STEP (H = STEP SIZE). |
|
|
908 C THIS MUST BE POSITIVE TO RESULT IN A NON-DEFAULT |
|
|
909 C VALUE. THE DEFAULT VALUE IS 10. |
|
|
910 C---------------------------------------------------------------------- |
|
|
911 C OPTIONAL OUTPUTS. |
|
|
912 C |
|
|
913 C AS OPTIONAL ADDITIONAL OUTPUT FROM ODESSA, THE VARIABLES LISTED |
|
|
914 C BELOW ARE QUANTITIES RELATED TO THE PERFORMANCE OF ODESSA |
|
|
915 C WHICH ARE AVAILABLE TO THE USER. THESE ARE COMMUNICATED BY WAY OF |
|
|
916 C THE WORK ARRAYS, BUT ALSO HAVE INTERNAL MNEMONIC NAMES AS SHOWN. |
|
|
917 C EXCEPT WHERE STATED OTHERWISE, ALL OF THESE OUTPUTS ARE DEFINED |
|
|
918 C ON ANY SUCCESSFUL RETURN FROM ODESSA, AND ON ANY RETURN WITH |
|
|
919 C ISTATE = -1, -2, -4, -5, OR -6. ON AN ILLEGAL INPUT RETURN |
|
|
920 C (ISTATE = -3), THEY WILL BE UNCHANGED FROM THEIR EXISTING VALUES |
|
|
921 C (IF ANY), EXCEPT POSSIBLY FOR TOLSF, LENRW, AND LENIW. |
|
|
922 C ON ANY ERROR RETURN, OUTPUTS RELEVANT TO THE ERROR WILL BE DEFINED, |
|
|
923 C AS NOTED BELOW. |
|
|
924 C |
|
|
925 C NAME LOCATION MEANING |
|
|
926 C |
|
|
927 C HU RWORK(11) THE STEP SIZE IN T LAST USED (SUCCESSFULLY). |
|
|
928 C |
|
|
929 C HCUR RWORK(12) THE STEP SIZE TO BE ATTEMPTED ON THE NEXT STEP. |
|
|
930 C |
|
|
931 C TCUR RWORK(13) THE CURRENT VALUE OF THE INDEPENDENT VARIABLE |
|
|
932 C WHICH THE SOLVER HAS ACTUALLY REACHED, I.E. THE |
|
|
933 C CURRENT INTERNAL MESH POINT IN T. ON OUTPUT, TCUR |
|
|
934 C WILL ALWAYS BE AT LEAST AS FAR AS THE ARGUMENT |
|
|
935 C T, BUT MAY BE FARTHER (IF INTERPOLATION WAS DONE). |
|
|
936 C |
|
|
937 C TOLSF RWORK(14) A TOLERANCE SCALE FACTOR, GREATER THAN 1.0, |
|
|
938 C COMPUTED WHEN A REQUEST FOR TOO MUCH ACCURACY WAS |
|
|
939 C DETECTED (ISTATE = -3 IF DETECTED AT THE START OF |
|
|
940 C THE PROBLEM, ISTATE = -2 OTHERWISE). IF ITOL IS |
|
|
941 C LEFT UNALTERED BUT RTOL AND ATOL ARE UNIFORMLY |
|
|
942 C SCALED UP BY A FACTOR OF TOLSF FOR THE NEXT CALL, |
|
|
943 C THEN THE SOLVER IS DEEMED LIKELY TO SUCCEED. |
|
|
944 C (THE USER MAY ALSO IGNORE TOLSF AND ALTER THE |
|
|
945 C TOLERANCE PARAMETERS IN ANY OTHER WAY APPROPRIATE.) |
|
|
946 C |
|
|
947 C NST IWORK(11) THE NUMBER OF STEPS TAKEN FOR THE PROBLEM SO FAR. |
|
|
948 C |
|
|
949 C NFE IWORK(12) THE NUMBER OF F EVALUATIONS FOR THE PROBLEM SO FAR. |
|
|
950 C |
|
|
951 C NJE IWORK(13) THE NUMBER OF JACOBIAN EVALUATIONS (AND OF MATRIX |
|
|
952 C LU DECOMPOSITIONS IF ISOPT = 0) FOR THE PROBLEM SO |
|
|
953 C FAR. IF ISOPT = 1, THE NUMBER OF LU DECOMPOSITIONS |
|
|
954 C IS EQUAL TO NJE - NSPE (SEE BELOW). |
|
|
955 C |
|
|
956 C NQU IWORK(14) THE METHOD ORDER LAST USED (SUCCESSFULLY). |
|
|
957 C |
|
|
958 C NQCUR IWORK(15) THE ORDER TO BE ATTEMPTED ON THE NEXT STEP. |
|
|
959 C |
|
|
960 C IMXER IWORK(16) THE INDEX OF THE COMPONENT OF LARGEST MAGNITUDE IN |
|
|
961 C THE WEIGHTED LOCAL ERROR VECTOR (E(I,J)/EWT(I,J)), |
|
|
962 C ON AN ERROR RETURN WITH ISTATE = -4 OR -5. |
|
|
963 C |
|
|
964 C LENRW IWORK(17) THE LENGTH OF RWORK ACTUALLY REQUIRED. |
|
|
965 C THIS IS DEFINED ON NORMAL RETURNS AND ON AN ILLEGAL |
|
|
966 C INPUT RETURN FOR INSUFFICIENT STORAGE. |
|
|
967 C |
|
|
968 C LENIW IWORK(18) THE LENGTH OF IWORK ACTUALLY REQUIRED. |
|
|
969 C THIS IS DEFINED ON NORMAL RETURNS AND ON AN ILLEGAL |
|
|
970 C INPUT RETURN FOR INSUFFICIENT STORAGE. |
|
|
971 C |
|
|
972 C NDFE IWORK(19) THE NUMBER OF DF/DP (VECTOR) EVALUATIONS. |
|
|
973 C |
|
|
974 C NSPE IWORK(20) THE NUMBER OF CALLS TO SUBROUTINE ODESSA_SPRIME. EACH CALL |
|
|
975 C TO ODESSA_SPRIME REQUIRES A JACOBIAN EVALUATION, BUT NOT |
|
|
976 C AN LU DECOMPOSITION. |
|
|
977 C |
|
|
978 C THE FOLLOWING ARRAYS ARE SEGMENTS OF THE RWORK AND IWORK ARRAYS |
|
|
979 C WHICH MAY ALSO BE OF INTEREST TO THE USER AS OPTIONAL OUTPUTS. |
|
|
980 C FOR EACH ARRAY, THE TABLE BELOW GIVES ITS INTERNAL NAME, ITS BASE |
|
|
981 C ADDRESS IN RWORK OR IWORK, AND ITS DESCRIPTION. |
|
|
982 C |
|
|
983 C NAME BASE ADDRESS DESCRIPTION |
|
|
984 C |
|
|
985 C YH 21 IN RWORK THE NORDSIECK HISTORY ARRAY, OF SIZE NYH BY |
|
|
986 C (NQCUR + 1). FOR J = 0,1,...,NQCUR, COLUMN J+1 |
|
|
987 C OF YH CONTAINS HCUR**J/FACTORIAL(J) TIMES |
|
|
988 C THE J-TH DERIVATIVE OF THE INTERPOLATING |
|
|
989 C POLYNOMIAL CURRENTLY REPRESENTING THE SOLUTION, |
|
|
990 C EVALUATED AT T = TCUR. |
|
|
991 C |
|
|
992 C ACOR LENRW-NYH+1 ARRAY OF SIZE NYH USED FOR THE ACCUMULATED |
|
|
993 C IN RWORK CORRECTIONS ON EACH STEP, SCALED ON OUTPUT |
|
|
994 C TO REPRESENT THE ESTIMATED LOCAL ERROR IN Y |
|
|
995 C ON THE LAST STEP. THIS IS THE VECTOR E IN |
|
|
996 C THE DESCRIPTION OF THE ERROR CONTROL. |
|
|
997 C IT IS DEFINED ONLY ON A SUCCESSFUL RETURN |
|
|
998 C FROM ODESSA. |
|
|
999 C NRS LENIW-NPAR ARRAY OF SIZE NPAR+1, USED TO STORE THE |
|
|
1000 C IN IWORK ACCUMULATED NUMBER OF REPEATED STEPS DUE TO |
|
|
1001 C THE SENSITIVITY ANALYSIS.. |
|
|
1002 C NRS(1) = TOTAL NUMBER OF REPEATED STEPS, |
|
|
1003 C NRS(2),... = NUMBER OF REPEATED STEPS DUE TO |
|
|
1004 C MODEL PARAMETER 1,... |
|
|
1005 C |
|
|
1006 C---------------------------------------------------------------------- |
|
|
1007 C PART II. OTHER ROUTINES CALLABLE. |
|
|
1008 C |
|
|
1009 C THE FOLLOWING ARE OPTIONAL CALLS WHICH THE USER MAY MAKE TO |
|
|
1010 C GAIN ADDITIONAL CAPABILITIES IN CONJUNCTION WITH ODESSA. |
|
|
1011 C |
|
|
1012 C CALL ODESSA_SVCOM (RSAV, ISAV) STORE IN RSAV AND ISAV THE CONTENTS |
|
|
1013 C OF THE INTERNAL COMMON BLOCKS USED BY |
|
|
1014 C ODESSA (SEE PART III BELOW). |
|
|
1015 C RSAV MUST BE A REAL ARRAY OF LENGTH 222 |
|
|
1016 C OR MORE, AND ISAV MUST BE AN INTEGER |
|
|
1017 C ARRAY OF LENGTH 54 OR MORE. |
|
|
1018 C |
|
|
1019 C CALL ODESSA_RSCOM (RSAV, ISAV) RESTORE, FROM RSAV AND ISAV, THE CONTENTS |
|
|
1020 C OF THE INTERNAL COMMON BLOCKS USED BY |
|
|
1021 C ODESSA. PRESUMES A PRIOR CALL TO ODESSA_SVCOM |
|
|
1022 C WITH THE SAME ARGUMENTS. |
|
|
1023 C |
|
|
1024 C ODESSA_SVCOM AND ODESSA_RSCOM ARE USEFUL IF |
|
|
1025 C INTERRUPTING A RUN AND RESTARTING |
|
|
1026 C LATER, OR ALTERNATING BETWEEN TWO OR |
|
|
1027 C MORE PROBLEMS SOLVED WITH ODESSA. |
|
|
1028 C |
|
|
1029 C CALL ODESSA_INTDY(,,,,,) PROVIDE DERIVATIVES OF Y, OF VARIOUS |
|
|
1030 C (SEE BELOW) ORDERS, AT A SPECIFIED POINT T, IF |
|
|
1031 C DESIRED. IT MAY BE CALLED ONLY AFTER |
|
|
1032 C A SUCCESSFUL RETURN FROM ODESSA. |
|
|
1033 C |
|
|
1034 C THE DETAILED INSTRUCTIONS FOR USING ODESSA_INTDY ARE AS FOLLOWS. |
|
|
1035 C THE FORM OF THE CALL IS.. |
|
|
1036 C |
|
|
1037 C CALL ODESSA_INTDY (T, K, RWORK(21), NYH, DKY, IFLAG) |
|
|
1038 C |
|
|
1039 C THE INPUT PARAMETERS ARE.. |
|
|
1040 C |
|
|
1041 C T = VALUE OF INDEPENDENT VARIABLE WHERE ANSWERS ARE DESIRED |
|
|
1042 C (NORMALLY THE SAME AS THE T LAST RETURNED BY ODESSA). |
|
|
1043 C FOR VALID RESULTS, T MUST LIE BETWEEN TCUR - HU AND TCUR. |
|
|
1044 C (SEE OPTIONAL OUTPUTS FOR TCUR AND HU.) |
|
|
1045 C K = INTEGER ORDER OF THE DERIVATIVE DESIRED. K MUST SATISFY |
|
|
1046 C 0 .LE. K .LE. NQCUR, WHERE NQCUR IS THE CURRENT ORDER |
|
|
1047 C (SEE OPTIONAL OUTPUTS). THE CAPABILITY CORRESPONDING |
|
|
1048 C TO K = 0, I.E. COMPUTING Y(T), IS ALREADY PROVIDED |
|
|
1049 C BY ODESSA DIRECTLY. SINCE NQCUR .GE. 1, THE FIRST |
|
|
1050 C DERIVATIVE DY/DT IS ALWAYS AVAILABLE WITH ODESSA_INTDY. |
|
|
1051 C RWORK(21) = THE BASE ADDRESS OF THE HISTORY ARRAY YH. |
|
|
1052 C NYH = COLUMN LENGTH OF YH, EQUAL TO THE TOTAL NUMBER OF |
|
|
1053 C DEPENDENT VARIABLES. IF ISOPT = 0, NYH = N. IF ISOPT = 1, |
|
|
1054 C NYH = N * (NPAR + 1). |
|
|
1055 C |
|
|
1056 C THE OUTPUT PARAMETERS ARE.. |
|
|
1057 C |
|
|
1058 C DKY = A REAL ARRAY OF LENGTH NYH CONTAINING THE COMPUTED VALUE |
|
|
1059 C OF THE K-TH DERIVATIVE OF Y(T). |
|
|
1060 C IFLAG = INTEGER FLAG, RETURNED AS 0 IF K AND T WERE LEGAL, |
|
|
1061 C -1 IF K WAS ILLEGAL, AND -2 IF T WAS ILLEGAL. |
|
|
1062 C ON AN ERROR RETURN, A MESSAGE IS ALSO WRITTEN. |
|
|
1063 C---------------------------------------------------------------------- |
|
|
1064 C PART III. COMMON BLOCKS. |
|
|
1065 C |
|
|
1066 C IF ODESSA IS TO BE USED IN AN OVERLAY SITUATION, THE USER |
|
|
1067 C MUST DECLARE, IN THE PRIMARY OVERLAY, THE VARIABLES IN.. |
|
|
1068 C (1) THE CALL SEQUENCE TO ODESSA, |
|
|
1069 C (2) THE THREE INTERNAL COMMON BLOCKS |
|
|
1070 C /ODE001/ OF LENGTH 258 (219 DOUBLE PRECISION WORDS |
|
|
1071 C FOLLOWED BY 39 INTEGER WORDS), |
|
|
1072 C /ODE002/ OF LENGTH 14 (3 DOUBLE PRECISION WORDS FOLLOWED |
|
|
1073 C BY 11 INTEGER WORDS), |
|
|
1074 C |
|
|
1075 C IF ODESSA IS USED ON A SYSTEM IN WHICH THE CONTENTS OF INTERNAL |
|
|
1076 C COMMON BLOCKS ARE NOT PRESERVED BETWEEN CALLS, THE USER SHOULD |
|
|
1077 C DECLARE THE ABOVE THREE COMMON BLOCKS IN HIS MAIN PROGRAM TO INSURE |
|
|
1078 C THAT THEIR CONTENTS ARE PRESERVED. |
|
|
1079 C |
|
|
1080 C IF THE SOLUTION OF A GIVEN PROBLEM BY ODESSA IS TO BE INTERRUPTED |
|
|
1081 C AND THEN LATER CONTINUED, SUCH AS WHEN RESTARTING AN INTERRUPTED RUN |
|
|
1082 C OR ALTERNATING BETWEEN TWO OR MORE PROBLEMS, THE USER SHOULD SAVE, |
|
|
1083 C FOLLOWING THE RETURN FROM THE LAST ODESSA CALL PRIOR TO THE |
|
|
1084 C INTERRUPTION, THE CONTENTS OF THE CALL SEQUENCE VARIABLES AND THE |
|
|
1085 C INTERNAL COMMON BLOCKS, AND LATER RESTORE THESE VALUES BEFORE THE |
|
|
1086 C NEXT ODESSA CALL FOR THAT PROBLEM. TO SAVE AND RESTORE THE COMMON |
|
|
1087 C BLOCKS, USE SUBROUTINES ODESSA_SVCOM AND ODESSA_RSCOM (SEE PART II ABOVE). |
|
|
1088 C |
|
|
1089 C---------------------------------------------------------------------- |
|
|
1090 C PART IV. OPTIONALLY REPLACEABLE SOLVER ROUTINES. |
|
|
1091 C |
|
|
1092 C BELOW ARE DESCRIPTIONS OF TWO ROUTINES IN THE ODESSA PACKAGE WHICH |
|
|
1093 C RELATE TO THE MEASUREMENT OF ERRORS. EITHER ROUTINE CAN BE |
|
|
1094 C REPLACED BY A USER-SUPPLIED VERSION, IF DESIRED. HOWEVER, SINCE SUCH |
|
|
1095 C A REPLACEMENT MAY HAVE A MAJOR IMPACT ON PERFORMANCE, IT SHOULD BE |
|
|
1096 C DONE ONLY WHEN ABSOLUTELY NECESSARY, AND ONLY WITH GREAT CAUTION. |
|
|
1097 C (NOTE.. THE MEANS BY WHICH THE PACKAGE VERSION OF A ROUTINE IS |
|
|
1098 C SUPERSEDED BY THE USER-S VERSION MAY BE SYSTEM-DEPENDENT.) |
|
|
1099 C |
|
|
1100 C (A) ODESSA_EWSET. |
|
|
1101 C THE FOLLOWING SUBROUTINE IS CALLED JUST BEFORE EACH INTERNAL |
|
|
1102 C INTEGRATION STEP, AND SETS THE ARRAY OF ERROR WEIGHTS, EWT, AS |
|
|
1103 C DESCRIBED UNDER ITOL/RTOL/ATOL ABOVE.. |
|
|
1104 C SUBROUTINE ODESSA_EWSET (NYH, ITOL, RTOL, ATOL, YCUR, EWT) |
|
|
1105 C WHERE NEQ, ITOL, RTOL, AND ATOL ARE AS IN THE ODESSA CALL SEQUENCE, |
|
|
1106 C YCUR CONTAINS THE CURRENT DEPENDENT VARIABLE VECTOR, AND |
|
|
1107 C EWT IS THE ARRAY OF WEIGHTS SET BY ODESSA_EWSET. |
|
|
1108 C |
|
|
1109 C IF THE USER SUPPLIES THIS SUBROUTINE, IT MUST RETURN IN EWT(I) |
|
|
1110 C (I = 1,...,NYH) A POSITIVE QUANTITY SUITABLE FOR COMPARING ERRORS |
|
|
1111 C IN Y(I) TO. THE EWT ARRAY RETURNED BY ODESSA_EWSET IS PASSED TO THE |
|
|
1112 C ODESSA_VNORM ROUTINE (SEE BELOW), AND ALSO USED BY ODESSA IN THE COMPUTATION |
|
|
1113 C OF THE OPTIONAL OUTPUT IMXER, THE DIAGONAL JACOBIAN APPROXIMATION, |
|
|
1114 C AND THE INCREMENTS FOR DIFFERENCE QUOTIENT JACOBIANS. |
|
|
1115 C |
|
|
1116 C IN THE USER-SUPPLIED VERSION OF ODESSA_EWSET, IT MAY BE DESIRABLE TO USE |
|
|
1117 C THE CURRENT VALUES OF DERIVATIVES OF Y. DERIVATIVES UP TO ORDER NQ |
|
|
1118 C ARE AVAILABLE FROM THE HISTORY ARRAY YH, DESCRIBED ABOVE UNDER |
|
|
1119 C OPTIONAL OUTPUTS. IN ODESSA_EWSET, YH IS IDENTICAL TO THE YCUR ARRAY, |
|
|
1120 C EXTENDED TO NQ + 1 COLUMNS WITH A COLUMN LENGTH OF NYH AND SCALE |
|
|
1121 C FACTORS OF H**J/FACTORIAL(J). ON THE FIRST CALL FOR THE PROBLEM, |
|
|
1122 C GIVEN BY NST = 0, NQ IS 1 AND H IS TEMPORARILY SET TO 1.0. |
|
|
1123 C THE QUANTITIES NQ, NYH, H, AND NST CAN BE OBTAINED BY INCLUDING |
|
|
1124 C IN ODESSA_EWSET THE STATEMENTS.. |
|
|
1125 C DOUBLE PRECISION H, RLS |
|
|
1126 C COMMON /ODE001/ RLS(219),ILS(39) |
|
|
1127 C NQ = ILS(35) |
|
|
1128 C NYH = ILS(14) |
|
|
1129 C NST = ILS(36) |
|
|
1130 C H = RLS(213) |
|
|
1131 C THUS, FOR EXAMPLE, THE CURRENT VALUE OF DY/DT CAN BE OBTAINED AS |
|
|
1132 C YCUR(NYH+I)/H (I=1,...,N) (AND THE DIVISION BY H IS |
|
|
1133 C UNNECESSARY WHEN NST = 0). |
|
|
1134 C |
|
|
1135 C (B) ODESSA_VNORM. |
|
|
1136 C THE FOLLOWING IS A REAL FUNCTION ROUTINE WHICH COMPUTES THE WEIGHTED |
|
|
1137 C ROOT-MEAN-SQUARE NORM OF A VECTOR V.. |
|
|
1138 C D = ODESSA_VNORM (LV, V, W) |
|
|
1139 C WHERE.. |
|
|
1140 C LV = THE LENGTH OF THE VECTOR, |
|
|
1141 C V = REAL ARRAY OF LENGTH N CONTAINING THE VECTOR, |
|
|
1142 C W = REAL ARRAY OF LENGTH N CONTAINING WEIGHTS, |
|
|
1143 C D = SQRT( (1/N) * SUM(V(I)*W(I))**2 ). |
|
|
1144 C ODESSA_VNORM IS CALLED WITH LV = N AND WITH W(I) = 1.0/EWT(I), WHERE |
|
|
1145 C EWT IS AS SET BY SUBROUTINE ODESSA_EWSET. |
|
|
1146 C |
|
|
1147 C IF THE USER SUPPLIES THIS FUNCTION, IT SHOULD RETURN A NON-NEGATIVE |
|
|
1148 C VALUE OF ODESSA_VNORM SUITABLE FOR USE IN THE ERROR CONTROL IN ODESSA. |
|
|
1149 C NONE OF THE ARGUMENTS SHOULD BE ALTERED BY ODESSA_VNORM. |
|
|
1150 C FOR EXAMPLE, A USER-SUPPLIED ODESSA_VNORM ROUTINE MIGHT.. |
|
|
1151 C -SUBSTITUTE A MAX-NORM OF (V(I)*W(I)) FOR THE RMS-NORM, OR |
|
|
1152 C -IGNORE SOME COMPONENTS OF V IN THE NORM, WITH THE EFFECT OF |
|
|
1153 C SUPPRESSING THE ERROR CONTROL ON THOSE COMPONENTS OF Y. |
|
|
1154 C---------------------------------------------------------------------- |
|
|
1155 C OTHER ROUTINES IN THE ODESSA PACKAGE. |
|
|
1156 C |
|
|
1157 C IN ADDITION TO SUBROUTINE ODESSA, THE ODESSA PACKAGE INCLUDES THE |
|
|
1158 C FOLLOWING SUBROUTINES AND FUNCTION ROUTINES.. |
|
|
1159 C ODESSA_INTDY COMPUTES AN INTERPOLATED VALUE OF THE Y VECTOR AT T = TOUT. |
|
|
1160 C ODESSA_STODE IS THE CORE INTEGRATOR, WHICH DOES ONE STEP OF THE |
|
|
1161 C INTEGRATION AND THE ASSOCIATED ERROR CONTROL. |
|
|
1162 C ODESSA_STESA MANAGES THE SOLUTION OF THE SENSITIVITY FUNCTIONS. |
|
|
1163 C ODESSA_CFODE SETS ALL METHOD COEFFICIENTS AND TEST CONSTANTS. |
|
|
1164 C ODESSA_PREPJ COMPUTES AND PREPROCESSES THE JACOBIAN MATRIX J = DF/DY |
|
|
1165 C AND THE NEWTON ITERATION MATRIX P = I - H*L0*J. |
|
|
1166 C IT IS ALSO CALLED BY ODESSA_SPRIME (WITH JOPT = 1) TO JUST |
|
|
1167 C COMPUTE THE JACOBIAN MATRIX. |
|
|
1168 C ODESSA_PREPDF COMPUTES THE INHOMOGENEITY MATRIX DF/DP. |
|
|
1169 C ODESSA_SPRIME DEFINES THE SYSTEM OF SENSITIVITY EQUATIONS. |
|
|
1170 C ODESSA_SOLSY MANAGES SOLUTION OF LINEAR SYSTEM IN CHORD ITERATION. |
|
|
1171 C ODESSA_EWSET SETS THE ERROR WEIGHT VECTOR EWT BEFORE EACH STEP. |
|
|
1172 C ODESSA_VNORM COMPUTES THE WEIGHTED R.M.S. NORM OF A VECTOR. |
|
|
1173 C ODESSA_SVCOM AND ODESSA_RSCOM ARE USER-CALLABLE ROUTINES TO SAVE AND RESTORE, |
|
|
1174 C RESPECTIVELY, THE CONTENTS OF THE INTERNAL COMMON BLOCKS. |
|
|
1175 C DGETRF AND DGETRS ARE ROUTINES FROM LAPACK FOR SOLVING FULL |
|
|
1176 C SYSTEMS OF LINEAR ALGEBRAIC EQUATIONS. |
|
|
1177 C DGBTRF AND DGBTRS ARE ROUTINES FROM LAPACK FOR SOLVING BANDED |
|
|
1178 C LINEAR SYSTEMS. |
|
|
1179 C DAXPY, DSCAL, IDAMAX, AND DDOT ARE BASIC LINEAR ALGEBRA MODULES |
|
|
1180 C (BLAS) USED BY THE ABOVE LINPACK ROUTINES. |
|
|
1181 C D1MACH COMPUTES THE UNIT ROUNDOFF IN A MACHINE-INDEPENDENT MANNER. |
|
|
1182 C XERRWD, XSETUN, AND ODESSA_XSETF HANDLE THE PRINTING OF ALL ERROR |
|
|
1183 C MESSAGES AND WARNINGS. |
|
|
1184 C NOTE.. ODESSA_VNORM, IDAMAX, DDOT, AND D1MACH ARE FUNCTION ROUTINES. |
|
|
1185 C ALL THE OTHERS ARE SUBROUTINES. |
|
|
1186 C |
|
|
1187 C THE FORTRAN GENERIC INTRINSIC FUNCTIONS USED BY ODESSA ARE.. |
|
|
1188 C ABS, MAX, MIN, REAL, MOD, SIGN, SQRT, AND WRITE |
|
|
1189 C |
|
|
1190 C A BLOCK DATA SUBPROGRAM IS ALSO INCLUDED WITH THE PACKAGE, |
|
|
1191 C FOR LOADING SOME OF THE VARIABLES IN INTERNAL COMMON. |
|
|
1192 C |
|
|
1193 C---------------------------------------------------------------------- |
|
|
1194 C PART V. GENERAL REMARKS |
|
|
1195 C |
|
|
1196 C THIS SECTION HIGHLIGHTS THE BASIC DIFFERENCES BETWEEN THE ORIGINAL |
|
|
1197 C LSODE PACKAGE AND THE ODESSA MODIFICATION. THIS IS PROVIDED AS A |
|
|
1198 C SERVICE TO EXPERIENCED LSODE USERS TO EXPEDITE FAMILIARIZATION WITH |
|
|
1199 C ODESSA. |
|
|
1200 C |
|
|
1201 C (A). ORIGINAL SUBROUTINES AND FUNCTIONS. |
|
|
1202 C |
|
|
1203 C OF THE ORIGINAL 22 SUBROUTINES AND FUNCTIONS USED IN THE LSODE |
|
|
1204 C PACKAGE, ALL ARE USED BY ODESSA, WITH THE FOLLOWING HAVING BEEN |
|
|
1205 C MODIFIED.. |
|
|
1206 C |
|
|
1207 C LSODE THE ORIGINAL DRIVER SUBROUTINE FOR THE LSODE PACKAGE IS |
|
|
1208 C EXTENSIVELY MODIFIED AND RENAMED ODESSA, WHICH NOW |
|
|
1209 C CONTAINS A CALL TO ODESSA_SPRIME TO ESTABLISH INITIAL CONDITIONS |
|
|
1210 C FOR THE SENSITIVITY CALCULATIONS. |
|
|
1211 C |
|
|
1212 C ODESSA_STODE THE ONE STEP INTEGRATOR IS SLIGHTLY MODIFIED AND RETAINS |
|
|
1213 C ITS ORIGINAL NAME. IT NOW CONTAINS THE CALL TO ODESSA_STESA, |
|
|
1214 C AND ALSO CALLS ODESSA_SPRIME IF KFLAG .LE. -3. |
|
|
1215 C |
|
|
1216 C ODESSA_PREPJ ALSO NAMED ODESSA_PREPJ IN ODESSA IS SLIGHTLY MODIFIED TO ALLOW |
|
|
1217 C FOR THE CALCULATION OF JACOBIAN WITH NO PREPROCESSING |
|
|
1218 C (JOPT = 1). |
|
|
1219 C |
|
|
1220 C (B). NEW SUBROUTINES. |
|
|
1221 C |
|
|
1222 C IN ADDITION TO THE CHANGES NOTED ABOVE, THREE NEW SUBROUTINES |
|
|
1223 C HAVE BEEN INTRODUCED (SEE ODESSA_STESA, ODESSA_SPRIME, AND ODESSA_PREPDF AS DESCRIBED |
|
|
1224 C IN PART IV. ABOVE). |
|
|
1225 C |
|
|
1226 C (C). COMMON BLOCKS. |
|
|
1227 C |
|
|
1228 C /LS0001/ RETAINS THE SAME LENGTH AND IS RENAMED /ODE001/; |
|
|
1229 C HOWEVER THE REAL ARRAY ROWNS(209) IS SHORTENED TO A |
|
|
1230 C LENGTH OF (173) REAL WORDS, ALLOWING THE REMOVAL OF |
|
|
1231 C TESCO(3,12) WHICH IS NOW PASSED FROM ODESSA_STODE TO ODESSA_STESA. |
|
|
1232 C IN ADDITION, THE INTEGER ARRAY IOWNS(6) IS SHORTENED |
|
|
1233 C TO A LENGTH OF (4) INTEGER WORDS, ALLOWING THE REMOVAL |
|
|
1234 C OF IALTH AND LMAX WHICH ARE NOW PASSED FROM ODESSA_STODE TO |
|
|
1235 C ODESSA_STESA. |
|
|
1236 C |
|
|
1237 C /ODE002/ ADDED COMMON BLOCK FOR VARIABLES IMPORTANT TO |
|
|
1238 C SENSITIVITY ANALYSIS (SEE PART III. ABOVE). A BLOCK |
|
|
1239 C DATA PROGRAM IS NOT REQUIRED FOR THIS COMMON BLOCK. |
|
|
1240 C |
|
|
1241 C ODESSA_SVCOM,ODESSA_RSCOM THESE TWO SUBROUTINES ARE MODIFIED TO HANDLE |
|
|
1242 C COMMON BLOCK /ODE002/ AS WELL. |
|
|
1243 C |
|
|
1244 C (D). OPTIONAL INPUTS. |
|
|
1245 C |
|
|
1246 C THE FULL SET OF OPTIONAL INPUTS AVAILABLE IN LSODE IS ALSO |
|
|
1247 C AVAILABLE IN ODESSA, WITH THE EXCEPTION THAT THE NUMBER OF ODE'S |
|
|
1248 C IN THE MODEL (NEQ(1)), MAY NOT BE CHANGED DURING THE PROBLEM. |
|
|
1249 C IN ODESSA, NYH NOW REFERS TO THE TOTAL NUMBER OF FIRST-ORDER |
|
|
1250 C ODE'S (MODEL AND SENSITIVITY EQUATIONS) WHICH IS EQUAL TO |
|
|
1251 C NEQ(1) IF ISOPT = 0, OR NEQ(1)*(NEQ(2)+1) IF ISOPT = 1. |
|
|
1252 C NEQ(1), NEQ(2), AND NYH ARE NOT ALLOWED TO CHANGE DURING |
|
|
1253 C THE COURSE OF AN INTEGRATION. |
|
|
1254 C |
|
|
1255 C (E). OPTIONAL OUTPUTS. |
|
|
1256 C |
|
|
1257 C THE FULL SET OF OPTIONAL OUTPUTS AVAILABLE IN LSODE IS ALSO |
|
|
1258 C AVAILABLE IN ODESSA. IN ADDITION, IWORK(19) AND IWORK(20) ARE |
|
|
1259 C LOADED WITH NDFE AND NSPE, RESPECTIVELY, UPON OUTPUT. THE TOTAL |
|
|
1260 C NUMBER OF LU DECOMPOSITIONS OF THE PROCESSED JACOBIAN IS EQUAL |
|
|
1261 C TO NJE - NSPE. |
|
|
1262 C----------------------------------------------------------------------- |
|
|
1263 SUBROUTINE DODESSA (F, DF, NEQ, Y, PAR, T, TOUT, ITOL, RTOL, ATOL, |
|
|
1264 1 ITASK, ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JAC, MF) |
|
|
1265 IMPLICIT DOUBLE PRECISION (A-H,O-Z) |
|
|
1266 LOGICAL IHIT |
|
|
1267 EXTERNAL F, DF, JAC, ODESSA_PREPJ, ODESSA_SOLSY, ODESSA_PREPDF |
|
|
1268 DIMENSION NEQ(*), Y(*), PAR(*), RTOL(*), ATOL(*), IOPT(*), |
|
|
1269 1 RWORK(LRW), IWORK(LIW), MORD(2) |
|
|
1270 C----------------------------------------------------------------------- |
|
|
1271 C THIS IS THE SEPTEMBER 1, 1986 VERSION OF ODESSA.. |
|
|
1272 C AN ORDINARY DIFFERENTIAL EQUATION SOLVER WITH EXPLICIT SIMULTANEOUS |
|
|
1273 C SENSITIVITY ANALYSIS. |
|
|
1274 C |
|
|
1275 C THIS PACKAGE IS A MODIFICATION OF THE AUGUST 13, 1981 VERSION OF |
|
|
1276 C LSODE.. LIVERMORE SOLVER FOR ORDINARY DIFFERENTIAL EQUATIONS. |
|
|
1277 C THIS VERSION IS IN DOUBLE PRECISION. |
|
|
1278 C |
|
|
1279 C ODESSA SOLVES FOR THE FIRST-ORDER SENSITIVITY COEFFICIENTS.. |
|
|
1280 C DY(I)/DP, FOR A SINGLE PARAMETER, OR, |
|
|
1281 C DY(I)/DP(J), FOR MULTIPLE PARAMETERS, |
|
|
1282 C ASSOCIATED WITH A SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS.. |
|
|
1283 C DY(T)/DT = F(Y,T;P). |
|
|
1284 C----------------------------------------------------------------------- |
|
|
1285 C REFERENCES... |
|
|
1286 C |
|
|
1287 C 1. JORGE R. LEIS AND MARK A. KRAMER, THE SIMULTANEOUS SOLUTION AND |
|
|
1288 C EXPLICIT SENSITIVITY ANALYSIS OF SYSTEMS DESCRIBED BY ORDINARY |
|
|
1289 C DIFFERENTIAL EQUATIONS, SUBMITTED TO ACM TRANS. MATH. SOFTWARE, |
|
|
1290 C (1985). |
|
|
1291 C |
|
|
1292 C 2. JORGE R. LEIS AND MARK A. KRAMER, ODESSA - AN ORDINARY |
|
|
1293 C DIFFERENTIAL EQUATION SOLVER WITH EXPLICIT SIMULTANEOUS |
|
|
1294 C SENSITIVITY ANALYSIS, SUBMITTED TO ACM TRANS. MATH. SOFTWARE. |
|
|
1295 C (1985). |
|
|
1296 C |
|
|
1297 C 3. ALAN C. HINDMARSH, LSODE AND LSODI, TWO NEW INITIAL VALUE |
|
|
1298 C ORDINARY DIFFERENTIAL EQUATION SOLVERS, ACM-SIGNUM NEWSLETTER, |
|
|
1299 C VOL. 15, NO. 4 (1980), PP. 10-11. |
|
|
1300 C----------------------------------------------------------------------- |
|
|
1301 C THE FOLLOWING INTERNAL COMMON BLOCKS CONTAIN |
|
|
1302 C (A) VARIABLES WHICH ARE LOCAL TO ANY SUBROUTINE BUT WHOSE VALUES MUST |
|
|
1303 C BE PRESERVED BETWEEN CALLS TO THE ROUTINE (OWN VARIABLES), AND |
|
|
1304 C (B) VARIABLES WHICH ARE COMMUNICATED BETWEEN SUBROUTINES. |
|
|
1305 C THE STRUCTURE OF THE BLOCKS ARE AS FOLLOWS.. ALL REAL VARIABLES ARE |
|
|
1306 C LISTED FIRST, FOLLOWED BY ALL INTEGERS. WITHIN EACH TYPE, THE |
|
|
1307 C VARIABLES ARE GROUPED WITH THOSE LOCAL TO SUBROUTINE ODESSA FIRST, |
|
|
1308 C THEN THOSE LOCAL TO SUBROUTINE ODESSA_STODE, AND FINALLY THOSE USED |
|
|
1309 C FOR COMMUNICATION. THE BLOCKS ARE DECLARED IN SUBROUTINES ODESSA |
|
|
1310 C ODESSA_INTDY, ODESSA_STODE, ODESSA_STESA, ODESSA_PREPJ, ODESSA_PREPDF, |
|
|
1311 C AND ODESSA_SOLSY. GROUPS OF VARIABLES ARE REPLACED BY DUMMY ARRAYS IN |
|
|
1312 C THE COMMON DECLARATIONS IN ROUTINES WHERE THOSE VARIABLES ARE NOT USED. |
|
|
1313 C----------------------------------------------------------------------- |
|
|
1314 COMMON /ODE001/ TRET, ROWNS(173), |
|
|
1315 1 TESCO(3,12), CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, |
|
|
1316 2 ILLIN, INIT, LYH, LEWT, LACOR, LSAVF, LWM, LIWM, |
|
|
1317 3 MXSTEP, MXHNIL, NHNIL, NTREP, NSLAST, NYH, IOWNS(4), |
|
|
1318 4 IALTH, LMAX, ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, |
|
|
1319 5 MITER, MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU |
|
|
1320 COMMON /ODE002/ DUPS, DSMS, DDNS, |
|
|
1321 1 NPAR, LDFDP, LNRS, |
|
|
1322 2 ISOPT, NSV, NDFE, NSPE, IDF, IERSP, JOPT, KFLAGS |
|
|
1323 PARAMETER (ZERO=0.0D0,ONE=1.0D0,TWO=2.0D0,FOUR=4.0D0) |
|
|
1324 DATA MORD(1),MORD(2)/12,5/, MXSTP0/500/, MXHNL0/10/ |
|
|
1325 C----------------------------------------------------------------------- |
|
|
1326 C BLOCK A. |
|
|
1327 C THIS CODE BLOCK IS EXECUTED ON EVERY CALL. |
|
|
1328 C IT TESTS ISTATE AND ITASK FOR LEGALITY AND BRANCHES APPROPIATELY. |
|
|
1329 C IF ISTATE .GT. 1 BUT THE FLAG INIT SHOWS THAT INITIALIZATION HAS |
|
|
1330 C NOT YET BEEN DONE, AN ERROR RETURN OCCURS. |
|
|
1331 C IF ISTATE = 1 AND TOUT = T, JUMP TO BLOCK G AND RETURN IMMEDIATELY. |
|
|
1332 C----------------------------------------------------------------------- |
|
|
1333 IF (ISTATE .LT. 1 .OR. ISTATE .GT. 3) GO TO 601 |
|
|
1334 IF (ITASK .LT. 1 .OR. ITASK .GT. 5) GO TO 602 |
|
|
1335 IF (ISTATE .EQ. 1) GO TO 10 |
|
|
1336 IF (INIT .EQ. 0) GO TO 603 |
|
|
1337 IF (ISTATE .EQ. 2) GO TO 200 |
|
|
1338 GO TO 20 |
|
|
1339 10 INIT = 0 |
|
|
1340 IF (TOUT .EQ. T) GO TO 430 |
|
|
1341 20 NTREP = 0 |
|
|
1342 C----------------------------------------------------------------------- |
|
|
1343 C BLOCK B. |
|
|
1344 C THE NEXT CODE BLOCK IS EXECUTED FOR THE INITIAL CALL (ISTATE = 1), |
|
|
1345 C OR FOR A CONTINUATION CALL WITH PARAMETER CHANGES (ISTATE = 3). |
|
|
1346 C IT CONTAINS CHECKING OF ALL INPUTS AND VARIOUS INITIALIZATIONS. |
|
|
1347 C |
|
|
1348 C FIRST CHECK LEGALITY OF THE NON-OPTIONAL INPUTS NEQ, ITOL, IOPT, |
|
|
1349 C MF, ML, AND MU. |
|
|
1350 C----------------------------------------------------------------------- |
|
|
1351 IF (NEQ(1) .LE. 0) GO TO 604 |
|
|
1352 IF (ISTATE .EQ. 1) GO TO 25 |
|
|
1353 IF (NEQ(1) .NE. N) GO TO 605 |
|
|
1354 25 N = NEQ(1) |
|
|
1355 IF (ITOL .LT. 1 .OR. ITOL .GT. 4) GO TO 606 |
|
|
1356 DO 26 I = 1,3 |
|
|
1357 26 IF (IOPT(I) .LT. 0 .OR. IOPT(I) .GT. 1) GO TO 607 |
|
|
1358 ISOPT = IOPT(2) |
|
|
1359 IDF = IOPT(3) |
|
|
1360 NYH = N |
|
|
1361 NSV = 1 |
|
|
1362 METH = MF/10 |
|
|
1363 MITER = MF - 10*METH |
|
|
1364 IF (METH .LT. 1 .OR. METH .GT. 2) GO TO 608 |
|
|
1365 IF (MITER .LT. 0 .OR. MITER .GT. 5) GO TO 608 |
|
|
1366 IF (MITER .LE. 3) GO TO 30 |
|
|
1367 ML = IWORK(1) |
|
|
1368 MU = IWORK(2) |
|
|
1369 IF (ML .LT. 0 .OR. ML .GE. N) GO TO 609 |
|
|
1370 IF (MU .LT. 0 .OR. MU .GE. N) GO TO 610 |
|
|
1371 30 IF (ISOPT .EQ. 0) GO TO 32 |
|
|
1372 C CHECK LEGALITY OF THE NON-OPTIONAL INPUTS ISOPT, NPAR. |
|
|
1373 C COMPUTE NUMBER OF SOLUTION VECTORS AND TOTAL NUMBER OF EQUATIONS. |
|
|
1374 IF (NEQ(2) .LE. 0) GO TO 628 |
|
|
1375 IF (ISTATE .EQ. 1) GO TO 31 |
|
|
1376 IF (NEQ(2) .NE. NPAR) GO TO 629 |
|
|
1377 31 NPAR = NEQ(2) |
|
|
1378 NSV = NPAR + 1 |
|
|
1379 NYH = NSV * N |
|
|
1380 IF (MITER .EQ. 0 .OR. MITER .EQ. 3) GO TO 630 |
|
|
1381 C NEXT PROCESS AND CHECK THE OPTIONAL INPUTS. -------------------------- |
|
|
1382 32 IF (IOPT(1) .EQ. 1) GO TO 40 |
|
|
1383 MAXORD = MORD(METH) |
|
|
1384 MXSTEP = MXSTP0 |
|
|
1385 MXHNIL = MXHNL0 |
|
|
1386 IF (ISTATE .EQ. 1) H0 = ZERO |
|
|
1387 HMXI = ZERO |
|
|
1388 HMIN = ZERO |
|
|
1389 GO TO 60 |
|
|
1390 40 MAXORD = IWORK(5) |
|
|
1391 IF (MAXORD .LT. 0) GO TO 611 |
|
|
1392 IF (MAXORD .EQ. 0) MAXORD = 100 |
|
|
1393 MAXORD = MIN(MAXORD,MORD(METH)) |
|
|
1394 MXSTEP = IWORK(6) |
|
|
1395 IF (MXSTEP .LT. 0) GO TO 612 |
|
|
1396 IF (MXSTEP .EQ. 0) MXSTEP = MXSTP0 |
|
|
1397 MXHNIL = IWORK(7) |
|
|
1398 IF (MXHNIL .LT. 0) GO TO 613 |
|
|
1399 IF (MXHNIL .EQ. 0) MXHNIL = MXHNL0 |
|
|
1400 IF (ISTATE .NE. 1) GO TO 50 |
|
|
1401 H0 = RWORK(5) |
|
|
1402 IF ((TOUT - T)*H0 .LT. ZERO) GO TO 614 |
|
|
1403 50 HMAX = RWORK(6) |
|
|
1404 IF (HMAX .LT. ZERO) GO TO 615 |
|
|
1405 HMXI = ZERO |
|
|
1406 IF (HMAX .GT. ZERO) HMXI = ONE/HMAX |
|
|
1407 HMIN = RWORK(7) |
|
|
1408 IF (HMIN .LT. ZERO) GO TO 616 |
|
|
1409 C----------------------------------------------------------------------- |
|
|
1410 C SET WORK ARRAY POINTERS AND CHECK LENGTHS LRW AND LIW. |
|
|
1411 C POINTERS TO SEGMENTS OF RWORK AND IWORK ARE NAMED BY PREFIXING L TO |
|
|
1412 C THE NAME OF THE SEGMENT. E.G., THE SEGMENT YH STARTS AT RWORK(LYH). |
|
|
1413 C SEGMENTS OF RWORK (IN ORDER) ARE DENOTED YH, WM, EWT, SAVF, ACOR. |
|
|
1414 C WORK SPACE FOR DFDP IS CONTAINED IN ACOR. |
|
|
1415 C----------------------------------------------------------------------- |
|
|
1416 60 LYH = 21 |
|
|
1417 LWM = LYH + (MAXORD + 1)*NYH |
|
|
1418 IF (MITER .EQ. 0) LENWM = 0 |
|
|
1419 IF (MITER .EQ. 1 .OR. MITER .EQ. 2) LENWM = N*N + 2 |
|
|
1420 IF (MITER .EQ. 3) LENWM = N + 2 |
|
|
1421 IF (MITER .GE. 4) LENWM = (2*ML + MU + 1)*N + 2 |
|
|
1422 LEWT = LWM + LENWM |
|
|
1423 LSAVF = LEWT + NYH |
|
|
1424 LACOR = LSAVF + N |
|
|
1425 LDFDP = LACOR + N |
|
|
1426 LENRW = LACOR + NYH - 1 |
|
|
1427 IWORK(17) = LENRW |
|
|
1428 LIWM = 1 |
|
|
1429 LENIW = 20 + N |
|
|
1430 IF (MITER .EQ. 0 .OR. MITER .EQ. 3) LENIW = 20 |
|
|
1431 LNRS = LENIW + LIWM |
|
|
1432 IF (ISOPT .EQ. 1) LENIW = LNRS + NPAR |
|
|
1433 IWORK(18) = LENIW |
|
|
1434 IF (LENRW .GT. LRW) GO TO 617 |
|
|
1435 IF (LENIW .GT. LIW) GO TO 618 |
|
|
1436 C CHECK RTOL AND ATOL FOR LEGALITY. ------------------------------------ |
|
|
1437 RTOLI = RTOL(1) |
|
|
1438 ATOLI = ATOL(1) |
|
|
1439 DO 70 I = 1,NYH |
|
|
1440 IF (ITOL .GE. 3) RTOLI = RTOL(I) |
|
|
1441 IF (ITOL .EQ. 2 .OR. ITOL .EQ. 4) ATOLI = ATOL(I) |
|
|
1442 IF (RTOLI .LT. ZERO) GO TO 619 |
|
|
1443 IF (ATOLI .LT. ZERO) GO TO 620 |
|
|
1444 70 CONTINUE |
|
|
1445 IF (ISTATE .EQ. 1) GO TO 100 |
|
|
1446 C IF ISTATE = 3, SET FLAG TO SIGNAL PARAMETER CHANGES TO ODESSA_STODE. - |
|
|
1447 JSTART = -1 |
|
|
1448 IF (NQ .LE. MAXORD) GO TO 90 |
|
|
1449 C MAXORD WAS REDUCED BELOW NQ. COPY YH(*,MAXORD+2) INTO SAVF. --------- |
|
|
1450 DO 80 I = 1,N |
|
|
1451 80 RWORK(I+LSAVF-1) = RWORK(I+LWM-1) |
|
|
1452 C RELOAD WM(1) = RWORK(LWM), SINCE LWM MAY HAVE CHANGED. --------------- |
|
|
1453 90 IF (MITER .GT. 0) RWORK(LWM) = DSQRT(UROUND) |
|
|
1454 GO TO 200 |
|
|
1455 C----------------------------------------------------------------------- |
|
|
1456 C BLOCK C. |
|
|
1457 C THE NEXT BLOCK IS FOR THE INITIAL CALL ONLY (ISTATE = 1). |
|
|
1458 C IT CONTAINS ALL REMAINING INITIALIZATIONS, THE INITIAL CALL TO F, |
|
|
1459 C THE INITIAL CALL TO ODESSA_SPRIME IF ISOPT = 1, |
|
|
1460 C AND THE CALCULATION OF THE INITIAL STEP SIZE. |
|
|
1461 C THE ERROR WEIGHTS IN EWT ARE INVERTED AFTER BEING LOADED. |
|
|
1462 C----------------------------------------------------------------------- |
|
|
1463 100 UROUND = D1MACH(4) |
|
|
1464 TN = T |
|
|
1465 IF (ITASK .NE. 4 .AND. ITASK .NE. 5) GO TO 105 |
|
|
1466 TCRIT = RWORK(1) |
|
|
1467 IF ((TCRIT - TOUT)*(TOUT - T) .LT. ZERO) GO TO 625 |
|
|
1468 IF (H0 .NE. ZERO .AND. (T + H0 - TCRIT)*H0 .GT. ZERO) |
|
|
1469 1 H0 = TCRIT - T |
|
|
1470 105 JSTART = 0 |
|
|
1471 IF (MITER .GT. 0) RWORK(LWM) = DSQRT(UROUND) |
|
|
1472 NHNIL = 0 |
|
|
1473 NST = 0 |
|
|
1474 NJE = 0 |
|
|
1475 NSLAST = 0 |
|
|
1476 HU = ZERO |
|
|
1477 NQU = 0 |
|
|
1478 CCMAX = 0.3D0 |
|
|
1479 MAXCOR = 3 |
|
|
1480 IF (ISOPT .EQ. 1) MAXCOR = 4 |
|
|
1481 MSBP = 20 |
|
|
1482 MXNCF = 10 |
|
|
1483 C INITIAL CALL TO F. (LF0 POINTS TO YH(1,2) AND LOADS IN VALUES). |
|
|
1484 LF0 = LYH + NYH |
|
|
1485 CALL F (NEQ, T, Y, PAR, RWORK(LF0)) |
|
|
1486 NFE = 1 |
|
|
1487 DUPS = ZERO |
|
|
1488 DSMS = ZERO |
|
|
1489 DDNS = ZERO |
|
|
1490 NDFE = 0 |
|
|
1491 NSPE = 0 |
|
|
1492 IF (ISOPT .EQ. 0) GO TO 114 |
|
|
1493 C INITIALIZE COUNTS FOR REPEATED STEPS DUE TO SENSITIVITY ANALYSIS. |
|
|
1494 DO 110 J = 1,NSV |
|
|
1495 110 IWORK(J + LNRS - 1) = 0 |
|
|
1496 C LOAD THE INITIAL VALUE VECTOR IN YH. --------------------------------- |
|
|
1497 114 DO 115 I = 1,NYH |
|
|
1498 115 RWORK(I+LYH-1) = Y(I) |
|
|
1499 C LOAD AND INVERT THE EWT ARRAY. (H IS TEMPORARILY SET TO ONE.) ------- |
|
|
1500 NQ = 1 |
|
|
1501 H = ONE |
|
|
1502 CALL ODESSA_EWSET (NYH, ITOL, RTOL, ATOL, RWORK(LYH), RWORK(LEWT)) |
|
|
1503 DO 120 I = 1,NYH |
|
|
1504 IF (RWORK(I+LEWT-1) .LE. ZERO) GO TO 621 |
|
|
1505 120 RWORK(I+LEWT-1) = ONE/RWORK(I+LEWT-1) |
|
|
1506 IF (ISOPT .EQ. 0) GO TO 125 |
|
|
1507 C CALL ODESSA_SPRIME TO LOAD FIRST-ORDER SENSITIVITY DERIVATIVES INTO |
|
|
1508 C REMAINING YH(*,2) POSITIONS. |
|
|
1509 CALL ODESSA_SPRIME (NEQ, Y, RWORK(LYH), NYH, N, NSV, RWORK(LWM), |
|
|
1510 1 IWORK(LIWM), RWORK(LEWT), RWORK(LF0), RWORK(LACOR), |
|
|
1511 2 RWORK(LDFDP), PAR, F, JAC, DF, ODESSA_PREPJ, ODESSA_PREPDF) |
|
|
1512 IF (IERSP .EQ. -1) GO TO 631 |
|
|
1513 IF (IERSP .EQ. -2) GO TO 632 |
|
|
1514 C----------------------------------------------------------------------- |
|
|
1515 C THE CODING BELOW COMPUTES THE STEP SIZE, H0, TO BE ATTEMPTED ON THE |
|
|
1516 C FIRST STEP, UNLESS THE USER HAS SUPPLIED A VALUE FOR THIS. |
|
|
1517 C FIRST CHECK THAT TOUT - T DIFFERS SIGNIFICANTLY FROM ZERO. |
|
|
1518 C A SCALAR TOLERANCE QUANTITY TOL IS COMPUTED, AS MAX(RTOL(I)) |
|
|
1519 C IF THIS IS POSITIVE, OR MAX(ATOL(I)/ABS(Y(I))) OTHERWISE, ADJUSTED |
|
|
1520 C SO AS TO BE BETWEEN 100*UROUND AND 1.0E-3. ONLY THE ORIGINAL |
|
|
1521 C SOLUTION VECTOR IS CONSIDERED IN THIS CALCULATION (ISOPT = 0 OR 1). |
|
|
1522 C THEN THE COMPUTED VALUE H0 IS GIVEN BY.. |
|
|
1523 C NEQ |
|
|
1524 C H0**2 = TOL / ( W0**-2 + (1/NEQ) * SUM ( F(I)/YWT(I) )**2 ) |
|
|
1525 C 1 |
|
|
1526 C WHERE W0 = MAX ( ABS(T), ABS(TOUT) ), |
|
|
1527 C F(I) = I-TH COMPONENT OF INITIAL VALUE OF F, |
|
|
1528 C YWT(I) = EWT(I)/TOL (A WEIGHT FOR Y(I)). |
|
|
1529 C THE SIGN OF H0 IS INFERRED FROM THE INITIAL VALUES OF TOUT AND T. |
|
|
1530 C----------------------------------------------------------------------- |
|
|
1531 125 IF (H0 .NE. ZERO) GO TO 180 |
|
|
1532 TDIST = DABS(TOUT - T) |
|
|
1533 W0 = DMAX1(DABS(T),DABS(TOUT)) |
|
|
1534 IF (TDIST .LT. TWO*UROUND*W0) GO TO 622 |
|
|
1535 TOL = RTOL(1) |
|
|
1536 IF (ITOL .LE. 2) GO TO 140 |
|
|
1537 DO 130 I = 1,N |
|
|
1538 130 TOL = DMAX1(TOL,RTOL(I)) |
|
|
1539 140 IF (TOL .GT. ZERO) GO TO 160 |
|
|
1540 ATOLI = ATOL(1) |
|
|
1541 DO 150 I = 1,N |
|
|
1542 IF (ITOL .EQ. 2 .OR. ITOL .EQ. 4) ATOLI = ATOL(I) |
|
|
1543 AYI = DABS(Y(I)) |
|
|
1544 IF (AYI .NE. ZERO) TOL = DMAX1(TOL,ATOLI/AYI) |
|
|
1545 150 CONTINUE |
|
|
1546 160 TOL = DMAX1(TOL,100.0D0*UROUND) |
|
|
1547 TOL = DMIN1(TOL,0.001D0) |
|
|
1548 SUM = ODESSA_VNORM (N, RWORK(LF0), RWORK(LEWT)) |
|
|
1549 SUM = ONE/(TOL*W0*W0) + TOL*SUM**2 |
|
|
1550 H0 = ONE/DSQRT(SUM) |
|
|
1551 H0 = MIN(H0,TDIST) |
|
|
1552 H0 = DSIGN(H0,TOUT-T) |
|
|
1553 C ADJUST H0 IF NECESSARY TO MEET HMAX BOUND. --------------------------- |
|
|
1554 180 RH = DABS(H0)*HMXI |
|
|
1555 IF (RH .GT. ONE) H0 = H0/RH |
|
|
1556 C LOAD H WITH H0 AND SCALE YH(*,2) BY H0. ------------------------------ |
|
|
1557 H = H0 |
|
|
1558 DO 190 I = 1,NYH |
|
|
1559 190 RWORK(I+LF0-1) = H0*RWORK(I+LF0-1) |
|
|
1560 GO TO 270 |
|
|
1561 C----------------------------------------------------------------------- |
|
|
1562 C BLOCK D. |
|
|
1563 C THE NEXT CODE BLOCK IS FOR CONTINUATION CALLS ONLY (ISTATE = 2 OR 3) |
|
|
1564 C AND IS TO CHECK STOP CONDITIONS BEFORE TAKING A STEP. |
|
|
1565 C----------------------------------------------------------------------- |
|
|
1566 200 NSLAST = NST |
|
|
1567 GO TO (210, 250, 220, 230, 240), ITASK |
|
|
1568 210 IF ((TN - TOUT)*H .LT. ZERO) GO TO 250 |
|
|
1569 CALL ODESSA_INTDY (TOUT, 0, RWORK(LYH), NYH, Y, IFLAG) |
|
|
1570 IF (IFLAG .NE. 0) GO TO 627 |
|
|
1571 T = TOUT |
|
|
1572 GO TO 420 |
|
|
1573 220 TP = TN - HU*(ONE + 100.0D0*UROUND) |
|
|
1574 IF ((TP - TOUT)*H .GT. ZERO) GO TO 623 |
|
|
1575 IF ((TN - TOUT)*H .LT. ZERO) GO TO 250 |
|
|
1576 GO TO 400 |
|
|
1577 230 TCRIT = RWORK(1) |
|
|
1578 IF ((TN - TCRIT)*H .GT. ZERO) GO TO 624 |
|
|
1579 IF ((TCRIT - TOUT)*H .LT. ZERO) GO TO 625 |
|
|
1580 IF ((TN - TOUT)*H .LT. ZERO) GO TO 245 |
|
|
1581 CALL ODESSA_INTDY (TOUT, 0, RWORK(LYH), NYH, Y, IFLAG) |
|
|
1582 IF (IFLAG .NE. 0) GO TO 627 |
|
|
1583 T = TOUT |
|
|
1584 GO TO 420 |
|
|
1585 240 TCRIT = RWORK(1) |
|
|
1586 IF ((TN - TCRIT)*H .GT. ZERO) GO TO 624 |
|
|
1587 245 HMX = DABS(TN) + DABS(H) |
|
|
1588 IHIT = DABS(TN - TCRIT) .LE. 100.0D0*UROUND*HMX |
|
|
1589 IF (IHIT) GO TO 400 |
|
|
1590 TNEXT = TN + H*(ONE + FOUR*UROUND) |
|
|
1591 IF ((TNEXT - TCRIT)*H .LE. ZERO) GO TO 250 |
|
|
1592 H = (TCRIT - TN)*(ONE - FOUR*UROUND) |
|
|
1593 IF (ISTATE .EQ. 2) JSTART = -2 |
|
|
1594 C----------------------------------------------------------------------- |
|
|
1595 C BLOCK E. |
|
|
1596 C THE NEXT BLOCK IS NORMALLY EXECUTED FOR ALL CALLS AND CONTAINS |
|
|
1597 C THE CALL TO THE ONE-STEP CORE INTEGRATOR ODESSA_STODE. |
|
|
1598 C |
|
|
1599 C THIS IS A LOOPING POINT FOR THE INTEGRATION STEPS. |
|
|
1600 C |
|
|
1601 C FIRST CHECK FOR TOO MANY STEPS BEING TAKEN, UPDATE EWT (IF NOT AT |
|
|
1602 C START OF PROBLEM), CHECK FOR TOO MUCH ACCURACY BEING REQUESTED, AND |
|
|
1603 C CHECK FOR H BELOW THE ROUNDOFF LEVEL IN T. |
|
|
1604 C TOLSF IS CALCULATED CONSIDERING ALL SOLUTION VECTORS. |
|
|
1605 C----------------------------------------------------------------------- |
|
|
1606 250 CONTINUE |
|
|
1607 IF ((NST-NSLAST) .GE. MXSTEP) GO TO 500 |
|
|
1608 CALL ODESSA_EWSET (NYH, ITOL, RTOL, ATOL, RWORK(LYH), RWORK(LEWT)) |
|
|
1609 DO 260 I = 1,NYH |
|
|
1610 IF (RWORK(I+LEWT-1) .LE. ZERO) GO TO 510 |
|
|
1611 260 RWORK(I+LEWT-1) = ONE/RWORK(I+LEWT-1) |
|
|
1612 270 TOLSF = UROUND*ODESSA_VNORM (NYH, RWORK(LYH), RWORK(LEWT)) |
|
|
1613 IF (TOLSF .LE. ONE) GO TO 280 |
|
|
1614 TOLSF = TOLSF*2.0D0 |
|
|
1615 IF (NST .EQ. 0) GO TO 626 |
|
|
1616 GO TO 520 |
|
|
1617 280 IF (ODESSA_ADDX(TN,H) .NE. TN) GO TO 290 |
|
|
1618 NHNIL = NHNIL + 1 |
|
|
1619 IF (NHNIL .GT. MXHNIL) GO TO 290 |
|
|
1620 CALL XERRWD ('ODESSA - WARNING..INTERNAL T (=R1) AND H (=R2) ARE', |
|
|
1621 1 50, 101, 1, 0, 0, 0, 0, ZERO, ZERO) |
|
|
1622 CALL XERRWD |
|
|
1623 1 ('SUCH THAT IN THE MACHINE, T + H = T ON THE NEXT STEP', |
|
|
1624 1 52, 101, 1, 0, 0, 0, 0, ZERO, ZERO) |
|
|
1625 CALL XERRWD ('(H = STEP SIZE). SOLVER WILL CONTINUE ANYWAY', |
|
|
1626 1 44, 101, 1, 0, 0, 0, 2, TN, H) |
|
|
1627 IF (NHNIL .LT. MXHNIL) GO TO 290 |
|
|
1628 CALL XERRWD ('ODESSA - ABOVE WARNING HAS BEEN ISSUED I1 TIMES.', |
|
|
1629 1 48, 102, 1, 0, 0, 0, 0, ZERO, ZERO) |
|
|
1630 CALL XERRWD ('IT WILL NOT BE ISSUED AGAIN FOR THIS PROBLEM', |
|
|
1631 1 44, 102, 1, 1, MXHNIL, 0, 0, ZERO,ZERO) |
|
|
1632 290 CONTINUE |
|
|
1633 C----------------------------------------------------------------------- |
|
|
1634 C CALL ODESSA_STODE(NEQ,Y,YH,NYH,YH,WM,IWM,EWT,SAVF,ACOR,PAR,NRS, |
|
|
1635 C 1 F,JAC,DF,ODESSA_PREPJ,ODESSA_PREPDF,ODESSA_SOLSY) |
|
|
1636 C----------------------------------------------------------------------- |
|
|
1637 CALL ODESSA_STODE (NEQ, Y, RWORK(LYH), NYH, RWORK(LYH), |
|
|
1638 1 RWORK(LWM), IWORK(LIWM), RWORK(LEWT), RWORK(LSAVF), |
|
|
1639 2 RWORK(LACOR), PAR, IWORK(LNRS), F, JAC, DF, ODESSA_PREPJ, |
|
|
1640 3 ODESSA_PREPDF, ODESSA_SOLSY) |
|
|
1641 KGO = 1 - KFLAG |
|
|
1642 GO TO (300, 530, 540, 633), KGO |
|
|
1643 C----------------------------------------------------------------------- |
|
|
1644 C BLOCK F. |
|
|
1645 C THE FOLLOWING BLOCK HANDLES THE CASE OF A SUCCESSFUL RETURN FROM THE |
|
|
1646 C CORE INTEGRATOR (KFLAG = 0). TEST FOR STOP CONDITIONS. |
|
|
1647 C----------------------------------------------------------------------- |
|
|
1648 300 INIT = 1 |
|
|
1649 GO TO (310, 400, 330, 340, 350), ITASK |
|
|
1650 C ITASK = 1. IF TOUT HAS BEEN REACHED, INTERPOLATE. ------------------- |
|
|
1651 310 IF ((TN - TOUT)*H .LT. ZERO) GO TO 250 |
|
|
1652 CALL ODESSA_INTDY (TOUT, 0, RWORK(LYH), NYH, Y, IFLAG) |
|
|
1653 T = TOUT |
|
|
1654 GO TO 420 |
|
|
1655 C ITASK = 3. JUMP TO EXIT IF TOUT WAS REACHED. ------------------------ |
|
|
1656 330 IF ((TN - TOUT)*H .GE. ZERO) GO TO 400 |
|
|
1657 GO TO 250 |
|
|
1658 C ITASK = 4. SEE IF TOUT OR TCRIT WAS REACHED. ADJUST H IF NECESSARY. |
|
|
1659 340 IF ((TN - TOUT)*H .LT. ZERO) GO TO 345 |
|
|
1660 CALL ODESSA_INTDY (TOUT, 0, RWORK(LYH), NYH, Y, IFLAG) |
|
|
1661 T = TOUT |
|
|
1662 GO TO 420 |
|
|
1663 345 HMX = DABS(TN) + DABS(H) |
|
|
1664 IHIT = DABS(TN - TCRIT) .LE. 100.0D0*UROUND*HMX |
|
|
1665 IF (IHIT) GO TO 400 |
|
|
1666 TNEXT = TN + H*(ONE + FOUR*UROUND) |
|
|
1667 IF ((TNEXT - TCRIT)*H .LE. ZERO) GO TO 250 |
|
|
1668 H = (TCRIT - TN)*(ONE - FOUR*UROUND) |
|
|
1669 JSTART = -2 |
|
|
1670 GO TO 250 |
|
|
1671 C ITASK = 5. SEE IF TCRIT WAS REACHED AND JUMP TO EXIT. --------------- |
|
|
1672 350 HMX = DABS(TN) + DABS(H) |
|
|
1673 IHIT = DABS(TN - TCRIT) .LE. 100.0D0*UROUND*HMX |
|
|
1674 C----------------------------------------------------------------------- |
|
|
1675 C BLOCK G. |
|
|
1676 C THE FOLLOWING BLOCK HANDLES ALL SUCCESSFUL RETURNS FROM ODESSA. |
|
|
1677 C IF ITASK .NE. 1, Y IS LOADED FROM YH AND T IS SET ACCORDINGLY. |
|
|
1678 C ISTATE IS SET TO 2, THE ILLEGAL INPUT COUNTER IS ZEROED, AND THE |
|
|
1679 C OPTIONAL OUTPUTS ARE LOADED INTO THE WORK ARRAYS BEFORE RETURNING. |
|
|
1680 C IF ISTATE = 1 AND TOUT = T, THERE IS A RETURN WITH NO ACTION TAKEN, |
|
|
1681 C EXCEPT THAT IF THIS HAS HAPPENED REPEATEDLY, THE RUN IS TERMINATED. |
|
|
1682 C----------------------------------------------------------------------- |
|
|
1683 400 DO 410 I = 1,NYH |
|
|
1684 410 Y(I) = RWORK(I+LYH-1) |
|
|
1685 T = TN |
|
|
1686 IF (ITASK .NE. 4 .AND. ITASK .NE. 5) GO TO 420 |
|
|
1687 IF (IHIT) T = TCRIT |
|
|
1688 420 ISTATE = 2 |
|
|
1689 ILLIN = 0 |
|
|
1690 RWORK(11) = HU |
|
|
1691 RWORK(12) = H |
|
|
1692 RWORK(13) = TN |
|
|
1693 IWORK(11) = NST |
|
|
1694 IWORK(12) = NFE |
|
|
1695 IWORK(13) = NJE |
|
|
1696 IWORK(14) = NQU |
|
|
1697 IWORK(15) = NQ |
|
|
1698 IF (ISOPT .EQ. 0) RETURN |
|
|
1699 IWORK(19) = NDFE |
|
|
1700 IWORK(20) = NSPE |
|
|
1701 RETURN |
|
|
1702 430 NTREP = NTREP + 1 |
|
|
1703 IF (NTREP .LT. 5) RETURN |
|
|
1704 CALL XERRWD ('ODESSA -- REPEATED CALLS WITH ISTATE = 1 AND |
|
|
1705 1 TOUT = T (=R1)', 59, 301, 1, 0, 0, 0, 1, T, ZERO) |
|
|
1706 GO TO 800 |
|
|
1707 C----------------------------------------------------------------------- |
|
|
1708 C BLOCK H. |
|
|
1709 C THE FOLLOWING BLOCK HANDLES ALL UNSUCCESSFUL RETURNS OTHER THAN |
|
|
1710 C THOSE FOR ILLEGAL INPUT. FIRST THE ERROR MESSAGE ROUTINE IS CALLED. |
|
|
1711 C IF THERE WAS AN ERROR TEST OR CONVERGENCE TEST FAILURE, IMXER IS SET. |
|
|
1712 C THEN Y IS LOADED FROM YH, T IS SET TO TN, AND THE ILLEGAL INPUT |
|
|
1713 C COUNTER ILLIN IS SET TO 0. THE OPTIONAL OUTPUTS ARE LOADED INTO |
|
|
1714 C THE WORK ARRAYS BEFORE RETURNING. |
|
|
1715 C----------------------------------------------------------------------- |
|
|
1716 C THE MAXIMUM NUMBER OF STEPS WAS TAKEN BEFORE REACHING TOUT. ---------- |
|
|
1717 500 CALL XERRWD ('ODESSA - AT CURRENT T (=R1), MXSTEP (=I1) STEPS', |
|
|
1718 1 47, 201, 1, 0, 0, 0, 0, ZERO,ZERO) |
|
|
1719 CALL XERRWD ('TAKEN ON THIS CALL BEFORE REACHING TOUT', |
|
|
1720 1 39, 201, 1, 1, MXSTEP, 0, 1, TN, ZERO) |
|
|
1721 ISTATE = -1 |
|
|
1722 GO TO 580 |
|
|
1723 C EWT(I) .LE. 0.0 FOR SOME I (NOT AT START OF PROBLEM). ---------------- |
|
|
1724 510 EWTI = RWORK(LEWT+I-1) |
|
|
1725 CALL XERRWD ('ODESSA - AT T (=R1), EWT(I1) HAS BECOME R2 .LE. 0.', |
|
|
1726 1 50, 202, 1, 1, I, 0, 2, TN, EWTI) |
|
|
1727 ISTATE = -6 |
|
|
1728 GO TO 580 |
|
|
1729 C TOO MUCH ACCURACY REQUESTED FOR MACHINE PRECISION. ------------------- |
|
|
1730 520 CALL XERRWD ('ODESSA - AT T (=R1), TOO MUCH ACCURACY REQUESTED', |
|
|
1731 1 48, 203, 1, 0, 0, 0, 0, ZERO,ZERO) |
|
|
1732 CALL XERRWD ('FOR PRECISION OF MACHINE.. SEE TOLSF (=R2)', |
|
|
1733 1 43, 203, 1, 0, 0, 0, 2, TN, TOLSF) |
|
|
1734 RWORK(14) = TOLSF |
|
|
1735 ISTATE = -2 |
|
|
1736 GO TO 580 |
|
|
1737 C KFLAG = -1. ERROR TEST FAILED REPEATEDLY OR WITH ABS(H) = HMIN. ----- |
|
|
1738 530 CALL XERRWD ('ODESSA - AT T(=R1) AND STEP SIZE H(=R2), THE ERROR', |
|
|
1739 1 50, 204, 1, 0, 0, 0, 0, ZERO, ZERO) |
|
|
1740 CALL XERRWD ('TEST FAILED REPEATEDLY OR WITH ABS(H) = HMIN', |
|
|
1741 1 44, 204, 1, 0, 0, 0, 2, TN, H) |
|
|
1742 ISTATE = -4 |
|
|
1743 GO TO 560 |
|
|
1744 C KFLAG = -2. CONVERGENCE FAILED REPEATEDLY OR WITH ABS(H) = HMIN. ---- |
|
|
1745 540 CALL XERRWD ('ODESSA - AT T (=R1) AND STEP SIZE H (=R2), THE', |
|
|
1746 1 46, 205, 1, 0, 0, 0, 0, ZERO,ZERO) |
|
|
1747 CALL XERRWD ('CORRECTOR CONVERGENCE FAILED REPEATEDLY', |
|
|
1748 1 39, 205, 1, 0, 0, 0, 0, ZERO,ZERO) |
|
|
1749 CALL XERRWD ('OR WITH ABS(H) = HMIN', |
|
|
1750 1 21, 0, 1, 0, 0, 0, 2, TN, H) |
|
|
1751 ISTATE = -5 |
|
|
1752 C COMPUTE IMXER IF RELEVANT. ------------------------------------------- |
|
|
1753 560 BIG = ZERO |
|
|
1754 IMXER = 1 |
|
|
1755 DO 570 I = 1,NYH |
|
|
1756 SIZE = DABS(RWORK(I+LACOR-1)*RWORK(I+LEWT-1)) |
|
|
1757 IF (BIG .GE. SIZE) GO TO 570 |
|
|
1758 BIG = SIZE |
|
|
1759 IMXER = I |
|
|
1760 570 CONTINUE |
|
|
1761 IWORK(16) = IMXER |
|
|
1762 C SET Y VECTOR, T, ILLIN, AND OPTIONAL OUTPUTS. ------------------------ |
|
|
1763 580 DO 590 I = 1,NYH |
|
|
1764 590 Y(I) = RWORK(I+LYH-1) |
|
|
1765 T = TN |
|
|
1766 ILLIN = 0 |
|
|
1767 RWORK(11) = HU |
|
|
1768 RWORK(12) = H |
|
|
1769 RWORK(13) = TN |
|
|
1770 IWORK(11) = NST |
|
|
1771 IWORK(12) = NFE |
|
|
1772 IWORK(13) = NJE |
|
|
1773 IWORK(14) = NQU |
|
|
1774 IWORK(15) = NQ |
|
|
1775 IF (ISOPT .EQ. 0) RETURN |
|
|
1776 IWORK(19) = NDFE |
|
|
1777 IWORK(20) = NSPE |
|
|
1778 RETURN |
|
|
1779 C----------------------------------------------------------------------- |
|
|
1780 C BLOCK I. |
|
|
1781 C THE FOLLOWING BLOCK HANDLES ALL ERROR RETURNS DUE TO ILLEGAL INPUT |
|
|
1782 C (ISTATE = -3), AS DETECTED BEFORE CALLING THE CORE INTEGRATOR. |
|
|
1783 C FIRST THE ERROR MESSAGE ROUTINE IS CALLED. THEN IF THERE HAVE BEEN |
|
|
1784 C 5 CONSECUTIVE SUCH RETURNS JUST BEFORE THIS CALL TO THE SOLVER, |
|
|
1785 C THE RUN IS HALTED. |
|
|
1786 C----------------------------------------------------------------------- |
|
|
1787 601 CALL XERRWD ('ODESSA - ISTATE (=I1) ILLEGAL', |
|
|
1788 1 29, 1, 1, 1, ISTATE, 0, 0, ZERO,ZERO) |
|
|
1789 GO TO 700 |
|
|
1790 602 CALL XERRWD ('ODESSA - ITASK (=I1) ILLEGAL', |
|
|
1791 1 28, 2, 1, 1, ITASK, 0, 0, ZERO,ZERO) |
|
|
1792 GO TO 700 |
|
|
1793 603 CALL XERRWD ('ODESSA - ISTATE .GT. 1 BUT ODESSA NOT INITIALIZED', |
|
|
1794 1 49, 3, 1, 0, 0, 0, 0, ZERO,ZERO) |
|
|
1795 GO TO 700 |
|
|
1796 604 CALL XERRWD ('ODESSA - NEQ (=I1) .LT. 1', |
|
|
1797 1 25, 4, 1, 1, NEQ(1), 0, 0, ZERO,ZERO) |
|
|
1798 GO TO 700 |
|
|
1799 605 CALL XERRWD ('ODESSA - ISTATE = 3 AND NEQ CHANGED. (I1 TO I2)', |
|
|
1800 1 48, 5, 1, 2, N, NEQ(1), 0, ZERO,ZERO) |
|
|
1801 GO TO 700 |
|
|
1802 606 CALL XERRWD ('ODESSA - ITOL (=I1) ILLEGAL', |
|
|
1803 1 27, 6, 1, 1, ITOL, 0, 0, ZERO,ZERO) |
|
|
1804 GO TO 700 |
|
|
1805 607 CALL XERRWD ('ODESSA - IOPT (=I1) ILLEGAL', |
|
|
1806 1 27, 7, 1, 1, IOPT, 0, 0, ZERO,ZERO) |
|
|
1807 GO TO 700 |
|
|
1808 608 CALL XERRWD('ODESSA - MF (=I1) ILLEGAL', |
|
|
1809 1 25, 8, 1, 1, MF, 0, 0, ZERO,ZERO) |
|
|
1810 GO TO 700 |
|
|
1811 609 CALL XERRWD('ODESSA - ML (=I1) ILLEGAL.. .LT.0 OR .GE.NEQ (=I2)', |
|
|
1812 1 50, 9, 1, 2, ML, NEQ(1), 0, ZERO,ZERO) |
|
|
1813 GO TO 700 |
|
|
1814 610 CALL XERRWD('ODESSA - MU (=I1) ILLEGAL.. .LT.0 OR .GE.NEQ (=I2)', |
|
|
1815 1 50, 10, 1, 2, MU, NEQ(1), 0, ZERO,ZERO) |
|
|
1816 GO TO 700 |
|
|
1817 611 CALL XERRWD('ODESSA - MAXORD (=I1) .LT. 0', |
|
|
1818 1 28, 11, 1, 1, MAXORD, 0, 0, ZERO,ZERO) |
|
|
1819 GO TO 700 |
|
|
1820 612 CALL XERRWD('ODESSA - MXSTEP (=I1) .LT. 0', |
|
|
1821 1 28, 12, 1, 1, MXSTEP, 0, 0, ZERO,ZERO) |
|
|
1822 GO TO 700 |
|
|
1823 613 CALL XERRWD('ODESSA - MXHNIL (=I1) .LT. 0', |
|
|
1824 1 28, 13, 1, 1, MXHNIL, 0, 0, ZERO,ZERO) |
|
|
1825 GO TO 700 |
|
|
1826 614 CALL XERRWD('ODESSA - TOUT (=R1) BEHIND T (=R2)', |
|
|
1827 1 34, 14, 1, 0, 0, 0, 2, TOUT, T) |
|
|
1828 CALL XERRWD('INTEGRATION DIRECTION IS GIVEN BY H0 (=R1)', |
|
|
1829 1 42, 14, 1, 0, 0, 0, 1, H0, ZERO) |
|
|
1830 GO TO 700 |
|
|
1831 615 CALL XERRWD('ODESSA - HMAX (=R1) .LT. 0.0', |
|
|
1832 1 28, 15, 1, 0, 0, 0, 1, HMAX, ZERO) |
|
|
1833 GO TO 700 |
|
|
1834 616 CALL XERRWD('ODESSA - HMIN (=R1) .LT. 0.0', |
|
|
1835 1 28, 16, 1, 0, 0, 0, 1, HMIN, ZERO) |
|
|
1836 GO TO 700 |
|
|
1837 617 CALL XERRWD('ODESSA - RWORK LENGTH NEEDED, LENRW (=I1), EXCEEDS |
|
|
1838 1 LRW (=I2)', 60, 17, 1, 2, LENRW, LRW, 0, ZERO,ZERO) |
|
|
1839 GO TO 700 |
|
|
1840 618 CALL XERRWD('ODESSA - IWORK LENGTH NEEDED, LENIW (=I1), EXCEEDS |
|
|
1841 1 LIW (=I2)', 60, 18, 1, 2, LENIW, LIW, 0, ZERO,ZERO) |
|
|
1842 GO TO 700 |
|
|
1843 619 CALL XERRWD('ODESSA - RTOL(I1) IS R1 .LT. 0.0', |
|
|
1844 1 32, 19, 1, 1, I, 0, 1, RTOLI, ZREO) |
|
|
1845 GO TO 700 |
|
|
1846 620 CALL XERRWD('ODESSA - ATOL(I1) IS R1 .LT. 0.0', |
|
|
1847 1 32, 20, 1, 1, I, 0, 1, ATOLI, ZERO) |
|
|
1848 GO TO 700 |
|
|
1849 * |
|
|
1850 621 EWTI = RWORK(LEWT+I-1) |
|
|
1851 CALL XERRWD('ODESSA - EWT(I1) IS R1 .LE. 0.0', |
|
|
1852 1 31, 21, 1, 1, I, 0, 1, EWTI, ZERO) |
|
|
1853 GO TO 700 |
|
|
1854 622 CALL XERRWD('ODESSA - TOUT (=R1) TOO CLOSE TO T(=R2) TO START |
|
|
1855 1 INTEGRATION', 60, 22, 1, 0, 0, 0, 2, TOUT, T) |
|
|
1856 GO TO 700 |
|
|
1857 623 CALL XERRWD('ODESSA - ITASK = I1 AND TOUT (=R1) BEHIND TCUR - HU |
|
|
1858 1 (= R2)', 58, 23, 1, 1, ITASK, 0, 2, TOUT, TP) |
|
|
1859 GO TO 700 |
|
|
1860 624 CALL XERRWD('ODESSA - ITASK = 4 OR 5 AND TCRIT (=R1) BEHIND TCUR |
|
|
1861 1 (=R2)', 57, 24, 1, 0, 0, 0, 2, TCRIT, TN) |
|
|
1862 GO TO 700 |
|
|
1863 625 CALL XERRWD('ODESSA - ITASK = 4 OR 5 AND TCRIT (=R1) BEHIND TOUT |
|
|
1864 1 (=R2)', 57, 25, 1, 0, 0, 0, 2, TCRIT, TOUT) |
|
|
1865 GO TO 700 |
|
|
1866 626 CALL XERRWD('ODESSA - AT START OF PROBLEM, TOO MUCH ACCURACY', |
|
|
1867 1 47, 26, 1, 0, 0, 0, 0, ZERO,ZERO) |
|
|
1868 CALL XERRWD('REQUESTED FOR PRECISION OF MACHINE. SEE TOLSF (=R1)', |
|
|
1869 1 51, 26, 1, 0, 0, 0, 1, TOLSF, ZERO) |
|
|
1870 RWORK(14) = TOLSF |
|
|
1871 GO TO 700 |
|
|
1872 627 CALL XERRWD |
|
|
1873 1 ('ODESSA - TROUBLE FROM ODESSA_INTDY. ITASK = I1, TOUT = R1', |
|
|
1874 1 57, 27, 1, 1, ITASK, 0, 1, TOUT, ZERO) |
|
|
1875 GO TO 700 |
|
|
1876 C ERROR STATEMENTS ASSOCIATED WITH SENSITIVITY ANALYSIS. |
|
|
1877 628 CALL XERRWD('ODESSA - NPAR (=I1) .LT. 1', |
|
|
1878 1 26, 28, 1, 1, NPAR, 0, 0, ZERO,ZERO) |
|
|
1879 GO TO 700 |
|
|
1880 629 CALL XERRWD('ODESSA - ISTATE = 3 AND NPAR CHANGED (I1 TO I2)', |
|
|
1881 1 47, 29, 1, 2, NP, NPAR, 0, ZERO,ZERO) |
|
|
1882 GO TO 700 |
|
|
1883 630 CALL XERRWD('ODESSA - MITER (=I1) ILLEGAL', |
|
|
1884 1 28, 30, 1, 1, MITER, 0, 0, ZERO,ZERO) |
|
|
1885 GO TO 700 |
|
|
1886 631 CALL XERRWD('ODESSA - TROUBLE IN ODESSA_SPRIME (IERPJ)', |
|
|
1887 1 41, 31, 1, 0, 0, 0, 0, ZERO,ZERO) |
|
|
1888 GO TO 700 |
|
|
1889 632 CALL XERRWD('ODESSA - TROUBLE IN ODESSA_SPRIME (MITER)', |
|
|
1890 1 41, 32, 1, 0, 0, 0, 0, ZERO,ZERO) |
|
|
1891 GO TO 700 |
|
|
1892 633 CALL XERRWD('ODESSA - FATAL ERROR IN ODESSA_STODE (KFLAG = -3)', |
|
|
1893 1 49, 33, 2, 0, 0, 0, 0, ZERO,ZERO) |
|
|
1894 GO TO 801 |
|
|
1895 C |
|
|
1896 700 IF (ILLIN .EQ. 5) GO TO 710 |
|
|
1897 ILLIN = ILLIN + 1 |
|
|
1898 ISTATE = -3 |
|
|
1899 RETURN |
|
|
1900 710 CALL XERRWD('ODESSA - REPEATED OCCURRENCES OF ILLEGAL INPUT', |
|
|
1901 1 46, 302, 1, 0, 0, 0, 0, ZERO,ZERO) |
|
|
1902 C |
|
|
1903 800 CALL XERRWD('ODESSA - RUN ABORTED.. APPARENT INFINITE LOOP', |
|
|
1904 1 45, 303, 2, 0, 0, 0, 0, ZERO,ZERO) |
|
|
1905 RETURN |
|
|
1906 801 CALL XERRWD('ODESSA - RUN ABORTED', |
|
|
1907 1 20, 304, 2, 0, 0, 0, 0, ZERO,ZERO) |
|
|
1908 RETURN |
|
|
1909 C-------------------- END OF SUBROUTINE ODESSA ------------------------- |
|
|
1910 END |
