Original WRF subgrid support version from John Michalakes without fire
[wrffire.git] / wrfv2_fire / chem / KPP / kpp / kpp-2.1 / int / atm_odessa_ddm.f
blob55bfd2e731bebdf55bf716c98792e071994f2e87
1 SUBROUTINE INTEGRATE( NSENSIT, Y, TIN, TOUT )
3 INCLUDE 'KPP_ROOT_Parameters.h'
4 INCLUDE 'KPP_ROOT_Global.h'
6 C TIN - Start Time
7 REAL*8 TIN
8 C TOUT - End Time
9 REAL*8 TOUT
10 C Concentrations and Sensitivities
11 REAL*8 Y(NVAR,NSENSIT+1), PARAMS(NSENSIT)
12 C --- Note: Y contains: (1:NVAR) concentrations, followed by
13 C --- (1:NVAR) sensitivities w.r.t. first parameter, followed by
14 C --- etc., followed by
15 C --- (1:NVAR) sensitivities w.r.t. NSENSIT's parameter
16 INTEGER i
18 INTEGER LIW, LRW
19 C PARAMETER (LRW = 22 + 8*(NSENSIT+1)*NVAR + NVAR**2 + NVAR)
20 C PARAMETER (LIW = 21 + NVAR + NSENSIT)
21 C REAL*8 RWORK(LRW)
22 C INTEGER IWORK(LIW)
23 C Note: the following dynamic allocation is not standard F77 and may not work on
24 C some systems. Declare LRW, LIW parameters as above with some upper bound
25 C used for NSENSIT
26 REAL*8 RWORK(22 + 8*(NSENSIT+1)*NVAR + NVAR**2 + NVAR)
27 INTEGER IWORK(21 + NVAR + NSENSIT)
29 INTEGER IOPT(3), NEQ(2)
31 EXTERNAL FUNC_CHEM,JAC,DFUNC_CHEMDPAR
33 MF = 21 ! --- BDF plus user-supplied Jacobian
35 LRW = 22 + 8*(NSENSIT+1)*NVAR + NVAR**2 + NVAR
36 LIW = 21 + NVAR + NSENSIT
38 NEQ(1) = NVAR ! --- No. of Variables
39 NEQ(2) = NSENSIT ! --- No of parameters
41 ITOL=1 ! --- 1=Scalar Tolerances; 4 = VECTOR TOLERANCES
42 ITASK=1 ! --- Normal Output
43 ISTATE=1
44 IOPT(1)=1 ! --- 0= No optional parameters, 1=Optional parameters
45 IOPT(2)=1 ! --- 1=Perform sensitivity analysis; 0 if not
46 IOPT(3)=0 ! --- DFUNC_CHEMDPAR supplied by the user;
47 ! --- 0 if finite differences are to be used
48 C --- Set optional parameters
49 DO 10 i=1,LRW
50 RWORK(i) = 0.0D0
51 10 CONTINUE
52 DO 20 i=1,LIW
53 IWORK(i) = 0
54 20 CONTINUE
56 RWORK(5) = STEPMIN ! THE STEP SIZE TO BE ATTEMPTED ON THE FIRST STEP.
57 RWORK(6) = STEPMAX ! THE MAXIMUM ABSOLUTE STEP SIZE ALLOWED.
58 RWORK(7) = 0.0D0 ! THE MINIMUM ABSOLUTE STEP SIZE ALLOWED.
59 IWORK(6) = 5000 ! MAXIMUM NUMBER OF (INTERNALLY DEFINED) STEPS
61 CALL ODESSA( FUNC_CHEM,DFUNC_CHEMDPAR,NEQ,Y,PARAMS,TIN,TOUT,
62 & ITOL,RTOL,ATOL,
63 1 ITASK,ISTATE,IOPT,RWORK,LRW,IWORK,LIW,
64 & JAC,MF)
66 IF (ISTATE.LT.0) THEN
67 print *,'ODESSA: Unsucessfull exit at T=',
68 & TIN,' (ISTATE=',ISTATE,')'
69 ENDIF
71 RETURN
72 END
76 SUBROUTINE FUNC_CHEM (N, T, V, PARAMS, FCT)
77 INCLUDE 'KPP_ROOT_Parameters.h'
78 INCLUDE 'KPP_ROOT_Global.h'
80 DIMENSION V(NVAR), PARAMS(*), FCT(NVAR)
81 TOLD = TIME
82 TIME = T
83 CALL Update_SUN()
84 CALL Update_RCONST()
85 TIME = TOLD
86 CALL Fun(V, FIX, RCONST, FCT)
87 RETURN
88 END
90 SUBROUTINE DFUNC_CHEMDPAR (N, T, V, PARAMS, DFCT, JPAR)
91 INCLUDE 'KPP_ROOT_Parameters.h'
92 INCLUDE 'KPP_ROOT_Global.h'
94 DIMENSION V(NVAR), PARAMS(*), DFCT(NVAR)
95 INTEGER JPAR, i
96 C This setting is required for sensitivities w.r.t. initial conditions
97 DO i=1,NVAR
98 DFCT(i) = 0.d0
99 END DO
100 RETURN
103 SUBROUTINE JAC (N, T, V, PARAMS, ML, MU, JV, NROWPD)
104 INCLUDE 'KPP_ROOT_Parameters.h'
105 INCLUDE 'KPP_ROOT_Global.h'
107 REAL*8 V(NVAR), PARAMS(*), JV(NVAR,NVAR)
108 INTEGER ML, MU, NROWPD
109 TOLD = TIME
110 TIME = T
111 CALL Update_SUN()
112 CALL Update_RCONST()
113 TIME = TOLD
114 DO i=1,NVAR
115 DO j=1,NVAR
116 JV(i,j) = 0.D0
117 END DO
118 END DO
119 CALL Jac(V, FIX, RCONST, JV)
120 RETURN
124 C ALGORITHM 658, COLLECTED ALGORITHMS FROM ACM.
125 C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
126 C VOL. 14, NO. 1, P.61.
127 C-----------------------------------------------------------------------
128 C THIS IS THE SEPTEMBER 1, 1986 VERSION OF ODESSA..
129 C AN ORDINARY DIFFERENTIAL EQUATION KppSolveR WITH EXPLICIT SIMULTANEOUS
130 C SENSITIVITY ANALYSIS.
132 C THIS PACKAGE IS A MODIFICATION OF THE AUGUST 13, 1981 VERSION OF
133 C LSODE.. LIVERMORE KppSolveR FOR ORDINARY DIFFERENTIAL EQUATIONS.
134 C THIS VERSION IS IN DOUBLE PRECISION.
136 C ODESSA KppSolveS FOR THE FIRST-ORDER SENSITIVITY COEFFICIENTS..
137 C DY(I)/DP, FOR A SINGLE PARAMETER, OR,
138 C DY(I)/DP(J), FOR MULTIPLE PARAMETERS,
139 C ASSOCIATED WITH A SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS..
140 C DY/DT = F(Y,T;P).
141 C-----------------------------------------------------------------------
142 C REFERENCES...
144 C 1. JORGE R. LEIS AND MARK A. KRAMER, THE SIMULTANEOUS SOLUTION AND
145 C EXPLICIT SENSITIVITY ANALYSIS OF SYSTEMS DESCRIBED BY ORDINARY
146 C DIFFERENTIAL EQUATIONS. SUBMITTED TO ACM TRANS. MATH. SOFTWARE,
147 C (1985).
149 C 2. JORGE R. LEIS AND MARK A. KRAMER, ODESSA - AN ORDINARY DIFFERENTIA
150 C EQUATION KppSolveR WITH EXPLICIT SIMULTANEOUS SENSITIVITY ANALYSIS.
151 C SUBMITTED TO ACM TRANS. MATH. SOFTWARE, (1985).
153 C 3. ALAN C. HINDMARSH, LSODE AND LSODI, TWO NEW INITIAL VALUE
154 C ORDINARY DIFFERENTIAL EQUATION KppSolveRS, ACM-SIGNUM NEWSLETTER,
155 C VOL. 15, NO. 4 (1980), PP. 10-11.
156 C-----------------------------------------------------------------------
157 C PROBLEM STATEMENT..
159 C THE ODESSA MODIFICATION OF THE LSODE PACKAGE PROVIDES THE OPTION TO
160 C CALCULATE FIRST-ORDER SENSITIVITY COEFFICIENTS FOR A SYSTEM OF STIFF
161 C OR NON-STIFF EXPLICIT ORDINARY DIFFERENTIAL EQUATIONS OF THE GENERAL
162 C FORM :
164 C DY/DT = F(Y,T;P) (1)
166 C WHERE Y IS AN N-DIMENSIONAL DEPENDENT VARIABLE VECTOR, T IS THE
167 C INDEPENDENT INTEGRATION VARIABLE, AND P IS AN NPAR-DIMENSIONAL
168 C CONSTANT VECTOR. THE GOVERNING EQUATIONS FOR THE FIRST-ORDER
169 C SENSITIVITY COEFFICIENTS ARE GIVEN BY :
171 C S'(T) = J(T)*S(T) + DF/DP (2)
173 C WHERE
175 C S(T) = DY(T)/DP (= SENSITIVITY FUNCTIONS)
176 C S'(T) = D(DY(T)/DP)/DT
177 C J(T) = DF(Y,T;P)/DY(T) (= JACOBIAN MATRIX)
178 C AND DF/DP = DF(Y,T;P)/DP (= INHOMOGENEITY MATRIX)
180 C SOLUTION OF EQUATIONS (1) AND (2) PROCEEDS SIMULTANEOUSLY VIA AN
181 C EXTENSION OF THE LSODE PACKAGE AS DESCRIBED IN [1].
182 C----------------------------------------------------------------------
183 C ACKNOWLEDGEMENT : THE FOLLOWING ODESSA PACKAGE DOCUMENTATION IS A
184 C MODIFICATION OF THE LSODE DOCUMENTATION WHICH
185 C ACCOMPANIES THE LSODE PACKAGE CODE.
186 C----------------------------------------------------------------------
187 C SUMMARY OF USAGE.
189 C COMMUNICATION BETWEEN THE USER AND THE ODESSA PACKAGE, FOR NORMAL
190 C SITUATIONS, IS SUMMARIZED HERE. THIS SUMMARY DESCRIBES ONLY A SUBSET
191 C OF THE FULL SET OF OPTIONS AVAILABLE. SEE THE FULL DESCRIPTION FOR
192 C DETAILS, INCLUDING OPTIONAL COMMUNICATION, NONSTANDARD OPTIONS,
193 C AND INSTRUCTIONS FOR SPECIAL SITUATIONS. SEE ALSO THE EXAMPLE
194 C PROBLEM (WITH PROGRAM AND OUTPUT) FOLLOWING THIS SUMMARY.
196 C A. FIRST PROVIDE A SUBROUTINE OF THE FORM..
197 C SUBROUTINE F (N, T, Y, PAR, YDOT)
198 C DOUBLE PRECISION T, Y, PAR, YDOT
199 C DIMENSION Y(N), YDOT(N), PAR(NPAR)
200 C WHICH SUPPLIES THE VECTOR FUNCTION F BY LOADING YDOT(I) WITH F(I).
201 C N IS THE NUMBER OF FIRST-ORDER DIFFERENTIAL EQUATIONS IN THE
202 C ABOVE MODEL. NPAR IS THE NUMBER OF MODEL PARAMETERS FOR WHICH
203 C VECTOR SENSITIVITY FUNCTIONS ARE DESIRED. YOU ARE ALSO ENCOURAGED
204 C TO PROVIDE A SUBROUTINE OF THE FORM..
205 C SUBROUTINE DF (N, T, Y, PAR, DFDP, JPAR)
206 C DOUBLE PRECISION T, Y, PAR, DFDP
207 C DIMENSION Y(N), PAR(NPAR), DFDP(N)
208 C GO TO (1,...,NPAR) JPAR
209 C 1 DFDP(1) = DF(1)/DP(1)
211 C DFDP(I) = DF(I)/DP(1)
213 C DFDP(N) = DF(N)/DP(1)
214 C RETURN
215 C 2 DFDP(1) = DF(1)/DP(2)
217 C DFDP(I) = DF(I)/DP(2)
219 C DFDP(N) = DF(N)/DP(2)
220 C RETURN
221 C . .
222 C . .
223 C RETURN
224 C NPAR DFDP(1) = DF(1)/DP(NPAR)
226 C DFDP(I) = DF(I)/DP(NPAR)
228 C DFDP(N) = DF(N)/DP(NPAR)
229 C RETURN
230 C END
231 C ONLY NONZERO ELEMENTS NEED BE LOADED. IF THIS IS NOT FEASIBLE,
232 C ODESSA WILL GENERATE THIS MATRIX INTERNALLY BY DIFFERENCE QUOTIENTS.
234 C B. NEXT DETERMINE (OR GUESS) WHETHER OR NOT THE PROBLEM IS STIFF.
235 C STIFFNESS OCCURS WHEN THE JACOBIAN MATRIX DF/DY HAS AN EIGENVALUE
236 C WHOSE REAL PART IS NEGATIVE AND LARGE IN MAGNITUDE, COMPARED TO THE
237 C RECIPROCAL OF THE T SPAN OF INTEREST. IF THE PROBLEM IS NONSTIFF,
238 C USE METH = 10. IF IT IS STIFF, USE METH = 20. THE USER IS REQUIRED
239 C TO INPUT THE METHOD FLAG MF = 10*METH + MITER. THERE ARE FOUR
240 C STANDARD CHOICES FOR MITER WHEN A SENSITIVITY ANALYSIS IS DESIRED,
241 C AND ODESSA REQUIRES THE JACOBIAN MATRIX IN SOME FORM.
242 C THIS MATRIX IS REGARDED EITHER AS FULL (MITER = 1 OR 2),
243 C OR BANDED (MITER = 4 OR 5). IN THE BANDED CASE, ODESSA REQUIRES TWO
244 C HALF-BANDWIDTH PARAMETERS ML AND MU. THESE ARE, RESPECTIVELY, THE
245 C WIDTHS OF THE LOWER AND UPPER PARTS OF THE BAND, EXCLUDING THE MAIN
246 C DIAGONAL. THUS THE BAND CONSISTS OF THE LOCATIONS (I,J) WITH
247 C I-ML .LE. J .LE. I+MU, AND THE FULL BANDWIDTH IS ML+MU+1.
249 C C. YOU ARE ENCOURAGED TO SUPPLY THE JACOBIAN DIRECTLY (MF = 11, 14,
250 C 21, OR 24), BUT IF THIS IS NOT FEASIBLE, ODESSA WILL COMPUTE IT
251 C INTERNALLY BY DIFFERENCE QUOTIENTS (MF = 12, 15, 22, OR 25). IF YOU
252 C ARE SUPPLYING THE JACOBIAN, PROVIDE A SUBROUTINE OF THE FORM..
253 C SUBROUTINE JAC (NEQ, T, Y, PAR, ML, MU, PD, NROWPD)
254 C DOUBLE PRECISION T, Y, PAR, PD
255 C DIMENSION Y(N), PD(NROWPD,N), PAR(NPAR)
256 C WHICH SUPPLIES DF/DY BY LOADING PD AS FOLLOWS..
257 C FOR A FULL JACOBIAN (MF = 11, OR 21), LOAD PD(I,J) WITH DF(I)/DY(J),
258 C THE PARTIAL DERIVATIVE OF F(I) WITH RESPECT TO Y(J). (IGNORE THE
259 C ML AND MU ARGUMENTS IN THIS CASE.)
260 C FOR A BANDED JACOBIAN (MF = 14, OR 24), LOAD PD(I-J+MU+1,J) WITH
261 C DF(I)/DY(J), I.E. LOAD THE DIAGONAL LINES OF DF/DY INTO THE ROWS OF
262 C PD FROM THE TOP DOWN.
263 C IN EITHER CASE, ONLY NONZERO ELEMENTS NEED BE LOADED.
265 C D. WRITE A MAIN PROGRAM WHICH CALLS SUBROUTINE ODESSA ONCE FOR
266 C EACH POINT AT WHICH ANSWERS ARE DESIRED. THIS SHOULD ALSO PROVIDE
267 C FOR POSSIBLE USE OF LOGICAL UNIT 6 FOR OUTPUT OF ERROR MESSAGES BY
268 C ODESSA. ON THE FIRST CALL TO ODESSA, SUPPLY ARGUMENTS AS FOLLOWS..
269 C F = NAME OF SUBROUTINE FOR RIGHT-HAND SIDE VECTOR F (MODEL).
270 C THIS NAME MUST BE DECLARED EXTERNAL IN CALLING PROGRAM.
271 C DF = NAME OF SUBROUTINE FOR INHOMOGENEITY MATRIX DF/DP.
272 C IF USED (IDF = 1), THIS NAME MUST BE DECLARED EXTERNAL IN
273 C CALLING PROGRAM. IF NOT USED (IDF = 0), PASS A DUMMY NAME.
274 C N = NUMBER OF FIRST ORDER ODE-S IN MODEL; LOAD INTO NEQ(1).
275 C NPAR = NUMBER OF MODEL PARAMETERS OF INTEREST; LOAD INTO NEQ(2).
276 C Y = AN (N) BY (NPAR+1) REAL ARRAY OF INITIAL VALUES..
277 C Y(I,1) , I = 1,N , CONTAIN THE STATE, OR MODEL, DEPENDENT
278 C VARIABLES,
279 C Y(I,J) , J = 2,NPAR , CONTAIN THE DEPENDENT SENSITIVITY
280 C COEFFICIENTS.
281 C PAR = A REAL ARRAY OF LENGTH NPAR CONTAINING MODEL PARAMETERS
282 C OF INTEREST.
283 C T = THE INITIAL VALUE OF THE INDEPENDENT VARIABLE.
284 C TOUT = FIRST POINT WHERE OUTPUT IS DESIRED (.NE. T).
285 C ITOL = 1, 2, 3, OR 4 ACCORDING AS RTOL, ATOL (BELOW) ARE SCALARS
286 C OR ARRAYS.
287 C RTOL = RELATIVE TOLERANCE PARAMETER (SCALAR OR (N) BY (NPAR+1)
288 C ARRAY).
289 C ATOL = ABSOLUTE TOLERANCE PARAMETER (SCALAR OR (N) BY (NPAR+1)
290 C ARRAY).
291 C THE ESTIMATED LOCAL ERROR IN Y(I,J) WILL BE CONTROLLED SO AS
292 C TO BE ROUGHLY LESS (IN MAGNITUDE) THAN
293 C EWT(I,J) = RTOL*ABS(Y(I,J)) + ATOL IF ITOL = 1,
294 C EWT(I,J) = RTOL*ABS(Y(I,J)) + ATOL(I,J) IF ITOL = 2,
295 C EWT(I,J) = RTOL(I,J)*ABS(Y(I,J) + ATOL IF ITOL = 3, OR
296 C EWT(I,J) = RTOL(I,J)*ABS(Y(I,J) + ATOL(I,J) IF ITOL = 4.
297 C THUS THE LOCAL ERROR TEST PASSES IF, IN EACH COMPONENT,
298 C EITHER THE ABSOLUTE ERROR IS LESS THAN ATOL (OR ATOL(I,J)),
299 C OR THE RELATIVE ERROR IS LESS THAN RTOL (OR RTOL(I,J)).
300 C USE RTOL = 0.0 FOR PURE ABSOLUTE ERROR CONTROL, AND
301 C USE ATOL = 0.0 FOR PURE RELATIVE ERROR CONTROL.
302 C CAUTION.. ACTUAL (GLOBAL) ERRORS MAY EXCEED THESE LOCAL
303 C TOLERANCES, SO CHOOSE THEM CONSERVATIVELY.
304 C ITASK = 1 FOR NORMAL COMPUTATION OF OUTPUT VALUES OF Y AT T = TOUT.
305 C ISTATE = INTEGER FLAG (INPUT AND OUTPUT). SET ISTATE = 1.
306 C IOPT = 0, TO INDICATE NO OPTIONAL INPUTS FOR INTEGRATION;
307 C LOAD INTO IOPT(1).
308 C ISOPT = 1, TO INDICATE SENSITIVITY ANALYSIS, = 0, TO INDICATE
309 C NO SENSITIVITY ANALYSIS; LOAD INTO IOPT(2).
310 C IDF = 1, IF SUBROUTINE DF (ABOVE) IS SUPPLIED BY THE USER,
311 C = 0, OTHERWISE; LOAD INTO IOPT(3).
312 C RWORK = REAL WORK ARRAY OF LENGTH AT LEAST..
313 C 22 + 16*N + N**2 FOR MF = 11 OR 12,
314 C 22 + 17*N + (2*ML + MU)*N FOR MF = 14 OR 15,
315 C 22 + 9*N + N**2 FOR MF = 21 OR 22,
316 C 22 + 10*N + (2*ML + MU)*N FOR MF = 24 OR 25,
317 C IF ISOPT = 0, OR..
318 C 22 + 15*(NPAR+1)*N + N**2 + N FOR MF = 11 OR 12,
319 C 24 + 15*(NPAR+1)*N + (2*ML+MU+2)*N + N FOR MF = 14 OR 15,
320 C 22 + 8*(NPAR+1)*N + N**2 + N FOR MF = 21 OR 22,
321 C 24 + 8*(NPAR+1)*N + (2*ML+MU+2)*N + N FOR MF = 24 OR 25,
322 C IF ISOPT = 1.
323 C LRW = DECLARED LENGTH OF RWORK (IN USER-S DIMENSION STATEMENT).
324 C IWORK = INTEGER WORK ARRAY OF LENGTH AT LEAST..
325 C 20 + N IF ISOPT = 0,
326 C 21 + N + NPAR IF ISOPT = 1.
327 C IF MITER = 4 OR 5, INPUT IN IWORK(1),IWORK(2) THE LOWER
328 C AND UPPER HALF-BANDWIDTHS ML,MU (EXCLUDING MAIN DIAGONAL).
329 C LIW = DECLARED LENGTH OF IWORK (IN USER-S DIMENSION STATEMENT).
330 C JAC = NAME OF SUBROUTINE FOR JACOBIAN MATRIX.
331 C IF USED, THIS NAME MUST BE DECLARED EXTERNAL IN CALLING
332 C PROGRAM. IF NOT USED, PASS A DUMMY NAME.
333 C MF = METHOD FLAG. STANDARD VALUES FOR ISOPT = 0 ARE..
334 C 10 FOR NONSTIFF (ADAMS) METHOD, NO JACOBIAN USED.
335 C 21 FOR STIFF (BDF) METHOD, USER-SUPPLIED FULL JACOBIAN.
336 C 22 FOR STIFF METHOD, INTERNALLY GENERATED FULL JACOBIAN.
337 C 24 FOR STIFF METHOD, USER-SUPPLIED BANDED JACOBIAN.
338 C 25 FOR STIFF METHOD, INTERNALLY GENERATED BANDED JACOBIAN.
339 C IF ISOPT = 1, MF = 10 IS ILLEGAL AND CAN BE REPLACED BY..
340 C 11 FOR NONSTIFF METHOD, USER-SUPPLIED FULL JACOBIAN.
341 C 12 FOR NONSTIFF METHOD, INTERNALLY GENERATED FULL JACOBIAN.
342 C 14 FOR NONSTIFF METHOD, USER-SUPPLIED BANDED JACOBIAN.
343 C 15 FOR NONSTIFF METHOD, INTERNALLY GENERATED BANDED JACOBIAN.
344 C NOTE THAT THE MAIN PROGRAM MUST DECLARE ARRAYS Y, RWORK, IWORK, AND
345 C POSSIBLY ATOL AND RTOL, AS WELL AS NEQ, IOPT, AND PAR IF ISOPT = 1.
347 C E. THE OUTPUT FROM THE FIRST CALL (OR ANY CALL) IS..
348 C Y = ARRAY OF COMPUTED VALUES OF Y(T) VECTOR.
349 C T = CORRESPONDING VALUE OF INDEPENDENT VARIABLE (NORMALLY TOUT).
350 C ISTATE = 2 IF ODESSA WAS SUCCESSFUL, NEGATIVE OTHERWISE.
351 C -1 MEANS EXCESS WORK DONE ON THIS CALL (PERHAPS WRONG MF).
352 C -2 MEANS EXCESS ACCURACY REQUESTED (TOLERANCES TOO SMALL).
353 C -3 MEANS ILLEGAL INPUT DETECTED (SEE PRINTED MESSAGE).
354 C -4 MEANS REPEATED ERROR TEST FAILURES (CHECK ALL INPUTS).
355 C -5 MEANS REPEATED CONVERGENCE FAILURES (PERHAPS BAD JACOBIAN
356 C SUPPLIED OR WRONG CHOICE OF MF OR TOLERANCES).
357 C -6 MEANS ERROR WEIGHT BECAME ZERO DURING PROBLEM. (SOLUTION
358 C COMPONENT I,J VANISHED, AND ATOL OR ATOL(I,J) = 0.0)
360 C F. TO CONTINUE THE INTEGRATION AFTER A SUCCESSFUL RETURN, SIMPLY
361 C RESET TOUT AND CALL ODESSA AGAIN. NO OTHER PARAMETERS NEED BE RESET.
362 C----------------------------------------------------------------------
363 C EXAMPLE PROBLEM.
365 C THE FOLLOWING IS A SIMPLE EXAMPLE PROBLEM, WITH THE CODING
366 C NEEDED FOR ITS SOLUTION BY ODESSA. THE PROBLEM IS FROM CHEMICAL
367 C KINETICS, AND CONSISTS OF THE FOLLOWING THREE RATE EQUATIONS..
368 C DY1/DT = -PAR(1)*Y1 + PAR(2)*Y2*Y3 ; PAR(1) = .04, PAR(2) = 1.E4
369 C DY2/DT = PAR(1)*Y1 - PAR(2)*Y2*Y3 - PAR(3)*Y2**2 ; PAR(3) = 3.E7
370 C DY3/DT = PAR(3)*Y2**2
371 C ON THE INTERVAL FROM T = 0.0 TO T = 4.E10, WITH INITIAL CONDITIONS
372 C Y1 = 1.0, Y2 = Y3 = 0, AND S(I,J) = 0, I = 1,3, J = 1,3.
373 C THE PROBLEM IS STIFF.
375 C THE FOLLOWING CODING KppSolveS THIS PROBLEM WITH ODESSA, USING
376 C MF = 21 AND PRINTING RESULTS AT T = .4, 4., ..., 4.E10.
377 C IT USES ITOL = 4 AND ATOL MUCH SMALLER FOR Y2 THAN Y1 OR Y3,
378 C BECAUSE Y2 HAS MUCH SMALLER VALUES. LESS STRINGENT TOLERANCES
379 C ARE ASSIGNED FOR THE SENSITIVITIES TO ACHIEVE GREATER EFFICIENCY.
380 C AT THE END OF THE RUN, STATISTICAL QUANTITIES OF INTEREST ARE
381 C PRINTED (SEE OPTIONAL OUTPUTS IN THE FULL DESCRIPTION BELOW).
383 C DOUBLE PRECISION ATOL, RWORK, RTOL, T, TOUT, Y, PAR
384 C EXTERNAL FEX, JEX, DFEX
385 C DIMENSION Y(3,4), PAR(3), ATOL(3,4), RTOL(3,4), RWORK(130),
386 C 1 IWORK(27), NEQ(2), IOPT(3)
387 C N = 3
388 C NPAR = 3
389 C NEQ(1) = N
390 C NEQ(2) = NPAR
391 C NSV = NPAR+1
392 C DO 10 I = 1,N
393 C DO 10 J = 1,NSV
394 C 10 Y(I,J) = 0.0D0
395 C Y(1,1) = 1.0D0
396 C PAR(1) = 0.04D0
397 C PAR(2) = 1.0D4
398 C PAR(3) = 3.0D7
399 C T = 0.D0
400 C TOUT = .4D0
401 C ITOL = 4
402 C ATOL(1,1) = 1.D-6
403 C ATOL(2,1) = 1.D-10
404 C ATOL(3,1) = 1.D-6
405 C DO 20 I = 1,N
406 C RTOL(I,1) = 1.D-4
407 C DO 15 J = 2,NSV
408 C RTOL(I,J) = 1.D-3
409 C 15 ATOL(I,J) = 1.D2 * ATOL(I,1)
410 C 20 CONTINUE
411 C ITASK = 1
412 C ISTATE = 1
413 C IOPT(1) = 0
414 C IOPT(2) = 1
415 C IOPT(3) = 1
416 C LRW = 130
417 C LIW = 27
418 C MF = 21
419 C DO 60 IOUT = 1,12
420 C CALL ODESSA(FEX,DFEX,NEQ,Y,PAR,T,TOUT,ITOL,RTOL,ATOL,
421 C 1 ITASK,ISTATE, IOPT,RWORK,LRW,IWORK,LIW,JEX,MF)
422 C WRITE(6,30)T,Y(1,1),Y(2,1),Y(3,1)
423 C 30 FORMAT(1X,7H AT T =,E12.4,6H Y =,3E14.6)
424 C DO 50 J = 2,NSV
425 C JPAR = J-1
426 C WRITE(6,40)JPAR,Y(1,J),Y(2,J),Y(3,J)
427 C 40 FORMAT(20X,2HS(,I1,3H) =,3E14.6)
428 C 50 CONTINUE
429 C IF (ISTATE .LT. 0) GO TO 80
430 C 60 TOUT = TOUT*10.D0
431 C WRITE(6,70)IWORK(11),IWORK(12),IWORK(13),IWORK(19)
432 C 70 FORMAT(1X,/,12H NO. STEPS =,I4,11H NO. F-S =,I4,11H NO. J-S =,
433 C 1 I4,12H NO. DF-S =,I4)
434 C STOP
435 C 80 WRITE(6,90)ISTATE
436 C 90 FORMAT(///22H ERROR HALT.. ISTATE =,I3)
437 C STOP
438 C END
440 C SUBROUTINE FEX (NEQ, T, Y, PAR, YDOT)
441 C DOUBLE PRECISION T, Y, YDOT, PAR
442 C DIMENSION Y(3), YDOT(3), PAR(3)
443 C YDOT(1) = -PAR(1)*Y(1) + PAR(2)*Y(2)*Y(3)
444 C YDOT(3) = PAR(3)*Y(2)*Y(2)
445 C YDOT(2) = -YDOT(1) - YDOT(3)
446 C RETURN
447 C END
449 C SUBROUTINE JEX (NEQ, T, Y, PAR, ML, MU, PD, NRPD)
450 C DOUBLE PRECISION PD, T, Y, PAR
451 C DIMENSION Y(3), PD(NRPD,3), PAR(3)
452 C PD(1,1) = -PAR(1)
453 C PD(1,2) = PAR(2)*Y(3)
454 C PD(1,3) = PAR(2)*Y(2)
455 C PD(2,1) = PAR(1)
456 C PD(2,3) = -PD(1,3)
457 C PD(3,2) = 2.D0*PAR(3)*Y(2)
458 C PD(2,2) = -PD(1,2) - PD(3,2)
459 C RETURN
460 C END
462 C SUBROUTINE DFEX (NEQ, T, Y, PAR, DFDP, JPAR)
463 C DOUBLE PRECISION T, Y, PAR, DFDP
464 C DIMENSION Y(3), PAR(3), DFDP(3)
465 C GO TO (1,2,3), JPAR
466 C 1 DFDP(1) = -Y(1)
467 C DFDP(2) = Y(1)
468 C RETURN
469 C 2 DFDP(1) = Y(2)*Y(3)
470 C DFDP(2) = -Y(2)*Y(3)
471 C RETURN
472 C 3 DFDP(2) = -Y(2)*Y(2)
473 C DFDP(3) = Y(2)*Y(2)
474 C RETURN
475 C END
477 C THE OUTPUT OF THIS PROGRAM (ON A DATA GENERAL MV-8000 IN
478 C DOUBLE PRECISION IS AS FOLLOWS:
480 C AT T = .4000E+00 Y = .985173E+00 .338641E-04 .147930E-01
481 C S(1) = -.355914E+00 .390261E-03 .355524E+00
482 C S(2) = .955150E-07 -.213065E-09 -.953019E-07
483 C S(3) = -.158466E-10 -.529012E-12 .163756E-10
484 C AT T = .4000E+01 Y = .905516E+00 .224044E-04 .944615E-01
485 C S(1) = -.187621E+01 .179197E-03 .187603E+01
486 C S(2) = .296093E-05 -.583104E-09 -.296034E-05
487 C S(3) = -.493267E-09 -.276246E-12 .493544E-09
488 C AT T = .4000E+02 Y = .715848E+00 .918628E-05 .284143E+00
489 C S(1) = -.424730E+01 .459360E-04 .424726E+01
490 C S(2) = .137294E-04 -.235815E-09 -.137291E-04
491 C S(3) = -.228818E-08 -.113803E-12 .228829E-08
492 C AT T = .4000E+03 Y = .450526E+00 .322299E-05 .549471E+00
493 C S(1) = -.595837E+01 .354310E-05 .595836E+01
494 C S(2) = .227380E-04 -.226041E-10 -.227380E-04
495 C S(3) = -.378971E-08 -.499501E-13 .378976E-08
496 C AT T = .4000E+04 Y = .183185E+00 .894131E-06 .816814E+00
497 C S(1) = -.475006E+01 -.599504E-05 .475007E+01
498 C S(2) = .188089E-04 .231330E-10 -.188089E-04
499 C S(3) = -.313478E-08 -.187575E-13 .313480E-08
500 C AT T = .4000E+05 Y = .389733E-01 .162133E-06 .961027E+00
501 C S(1) = -.157477E+01 -.276199E-05 .157477E+01
502 C S(2) = .628668E-05 .110026E-10 -.628670E-05
503 C S(3) = -.104776E-08 -.453588E-14 .104776E-08
504 C AT T = .4000E+06 Y = .493609E-02 .198411E-07 .995064E+00
505 C S(1) = -.236244E+00 -.458262E-06 .236244E+00
506 C S(2) = .944669E-06 .183193E-11 -.944671E-06
507 C S(3) = -.157441E-09 -.635990E-15 .157442E-09
508 C AT T = .4000E+07 Y = .516087E-03 .206540E-08 .999484E+00
509 C S(1) = -.256277E-01 -.509808E-07 .256278E-01
510 C S(2) = .102506E-06 .203905E-12 -.102506E-06
511 C S(3) = -.170825E-10 -.684002E-16 .170826E-10
512 C AT T = .4000E+08 Y = .519314E-04 .207736E-09 .999948E+00
513 C S(1) = -.259316E-02 -.518029E-08 .259316E-02
514 C S(2) = .103726E-07 .207209E-13 -.103726E-07
515 C S(3) = -.172845E-11 -.691450E-17 .172845E-11
516 C AT T = .4000E+09 Y = .544710E-05 .217885E-10 .999995E+00
517 C S(1) = -.271637E-03 -.541849E-09 .271638E-03
518 C S(2) = .108655E-08 .216739E-14 -.108655E-08
519 C S(3) = -.180902E-12 -.723615E-18 .180902E-12
520 C AT T = .4000E+10 Y = .446748E-06 .178699E-11 .100000E+01
521 C S(1) = -.322322E-04 -.842541E-10 .322323E-04
522 C S(2) = .128929E-09 .337016E-15 -.128929E-09
523 C S(3) = -.209715E-13 -.838859E-19 .209715E-13
524 C AT T = .4000E+11 Y = -.363960E-07 -.145584E-12 .100000E+01
525 C S(1) = -.164109E-06 -.429604E-11 .164113E-06
526 C S(2) = .656436E-12 .171842E-16 -.656451E-12
527 C S(3) = -.689361E-15 -.275745E-20 .689363E-15
529 C NO. STEPS = 340 NO. F-S = 412 NO. J-S = 343 NO. DF-S =1023
530 C----------------------------------------------------------------------
531 C FULL DESCRIPTION OF USER INTERFACE TO ODESSA.
533 C THE USER INTERFACE TO ODESSA CONSISTS OF THE FOLLOWING PARTS.
535 C I. THE CALL SEQUENCE TO SUBROUTINE ODESSA, WHICH IS A DRIVER
536 C ROUTINE FOR THE KppSolveR. THIS INCLUDES DESCRIPTIONS OF BOTH
537 C THE CALL SEQUENCE ARGUMENTS AND OF USER-SUPPLIED ROUTINES.
538 C FOLLOWING THESE DESCRIPTIONS IS A DESCRIPTION OF
539 C OPTIONAL INPUTS AVAILABLE THROUGH THE CALL SEQUENCE, AND THEN
540 C A DESCRIPTION OF OPTIONAL OUTPUTS (IN THE WORK ARRAYS).
542 C II. DESCRIPTIONS OF OTHER ROUTINES IN THE ODESSA PACKAGE THAT MAY
543 C BE (OPTIONALLY) CALLED BY THE USER. THESE PROVIDE THE ABILITY
544 C TO ALTER ERROR MESSAGE HANDLING, SAVE AND RESTORE THE INTERNAL
545 C COMMON, AND OBTAIN SPECIFIED DERIVATIVES OF THE SOLUTION Y(T).
547 C III. DESCRIPTIONS OF COMMON BLOCKS TO BE DECLARED IN OVERLAY
548 C OR SIMILAR ENVIRONMENTS, OR TO BE SAVED WHEN DOING AN INTERRUPT
549 C OF THE PROBLEM AND CONTINUED SOLUTION LATER.
551 C IV. DESCRIPTION OF TWO SUBROUTINES IN THE ODESSA PACKAGE, EITHER OF
552 C WHICH THE USER MAY REPLACE WITH HIS OWN VERSION, IF DESIRED.
553 C THESE RELATE TO THE MEASUREMENT OF ERRORS.
555 C V. GENERAL REMARKS WHICH HIGHLIGHT DIFFERENCES BETWEEN THE LSODE
556 C PACKAGE AND THE ODESSA PACKAGE.
557 C----------------------------------------------------------------------
558 C PART I. CALL SEQUENCE.
560 C THE CALL SEQUENCE PARAMETERS USED FOR INPUT ONLY ARE..
561 C F, DF, NEQ, PAR, TOUT, ITOL, RTOL, ATOL, ITASK, IOPT, LRW, LIW,
562 C JAC, MF,
563 C AND THOSE USED FOR BOTH INPUT AND OUTPUT ARE
564 C Y, T, ISTATE.
565 C THE WORK ARRAYS RWORK AND IWORK ARE ALSO USED FOR CONDITIONAL AND
566 C OPTIONAL INPUTS AND OPTIONAL OUTPUTS. (THE TERM OUTPUT HERE REFERS
567 C TO THE RETURN FROM SUBROUTINE ODESSA TO THE USER-S CALLING PROGRAM.)
569 C THE LEGALITY OF INPUT PARAMETERS WILL BE THOROUGHLY CHECKED ON THE
570 C INITIAL CALL FOR THE PROBLEM, BUT NOT CHECKED THEREAFTER UNLESS A
571 C CHANGE IN INPUT PARAMETERS IS FLAGGED BY ISTATE = 3 ON INPUT.
573 C THE DESCRIPTIONS OF THE CALL ARGUMENTS ARE AS FOLLOWS.
575 C F = THE NAME OF THE USER-SUPPLIED SUBROUTINE DEFINING THE
576 C ODE MODEL. THIS SYSTEM MUST BE PUT IN THE FIRST-ORDER
577 C FORM DY/DT = F(Y,T;P), WHERE F IS A VECTOR-VALUED FUNCTION
578 C OF THE SCALAR T AND VECTORS Y, AND PAR. SUBROUTINE F IS TO
579 C COMPUTE THE FUNCTION F. IT IS TO HAVE THE FORM..
580 C SUBROUTINE F (NEQ, T, Y, PAR, YDOT)
581 C DOUBLE PRECISION T, Y, PAR, YDOT
582 C DIMENSION Y(1), PAR(1), YDOT(1)
583 C WHERE NEQ, T, Y, AND PAR ARE INPUT, AND YDOT = F(Y,T;P)
584 C IS OUTPUT. Y AND YDOT ARE ARRAYS OF LENGTH N (= NEQ(1)).
585 C (IN THE DIMENSION STATEMENT ABOVE, 1 IS A DUMMY
586 C DIMENSION.. IT CAN BE REPLACED BY ANY VALUE.)
587 C F SHOULD NOT ALTER ARRAY Y, OR PAR(1),...,PAR(NPAR).
588 C F MUST BE DECLARED EXTERNAL IN THE CALLING PROGRAM.
590 C SUBROUTINE F MAY ACCESS USER-DEFINED QUANTITIES IN
591 C NEQ(2),... AND PAR(NPAR+1),... IF NEQ IS AN ARRAY
592 C (DIMENSIONED IN F) AND PAR HAS LENGTH EXCEEDING NPAR.
593 C SEE THE DESCRIPTIONS OF NEQ AND PAR BELOW.
595 C DF = THE NAME OF THE USER-SUPPLIED ROUTINE (IDF = 1) TO COMPUTE
596 C THE INHOMOGENEITY MATRIX, DF/DP, AS A FUNCTION OF THE SCALAR
597 C T, AND THE VECTORS Y, AND PAR. IT IS TO HAVE THE FORM
598 C SUBROUTINE DF (NEQ, T, Y, PAR, DFDP, JPAR)
599 C DOUBLE PRECISION T, Y, PAR, DFDP
600 C DIMENSION Y(1), PAR(1), DFDP(1)
601 C GO TO (1,2,...,NPAR) JPAR
602 C 1 DFDP(1) = DF(1)/DP(1)
604 C DFDP(I) = DF(I)/DP(1)
606 C DFDP(N) = DF(N)/DP(1)
607 C RETURN
608 C 2 DFDP(1) = DF(1)/DP(2)
610 C DFDP(I) = DF(I)/DP(2)
612 C DFDP(N) = DF(N)/DP(2)
614 C RETURN
615 C . .
616 C . .
617 C NPAR DFDP(1) = DF(1)/DP(NPAR)
619 C DFDP(I) = DF(I)/DP(NPAR)
621 C DFDP(N) = DF(N)/DP(NPAR)
622 C RETURN
623 C END
624 C WHERE NEQ, T, Y, PAR, AND JPAR ARE INPUT AND THE VECTOR
625 C DFDP(*,JPAR) IS TO BE LOADED WITH THE PARTIAL DERIVATIVES
626 C DF(Y,T;PAR)/DP(JPAR) ON OUTPUT. ONLY NONZERO ELEMENTS NEED
627 C BE LOADED. T, Y, AND PAR HAVE THE SAME MEANING AS IN
628 C SUBROUTINE F. (IN THE DIMENSION STATEMENT ABOVE, 1 IS A DUMMY
629 C DIMENSION.. IT CAN BE REPLACED BY ANY VALUE).
631 C DFDP(*,JPAR) IS PRESET TO ZERO BY THE KppSolveR, SO THAT ONLY
632 C THE NONZERO ELEMENTS NEED BE LOADED BY DF. SUBROUTINE DF
633 C MUST BE DECLARED EXTERNAL IN THE CALLING PROGRAM IF USED.
634 C IF IDF = 0 (OR ISOPT = 0), A DUMMY ARGUMENT CAN BE USED.
636 C SUBROUTINE DF MAY ACCESS USER-DEFINED QUANTITIES IN
637 C NEQ(2),... AND PAR(NPAR+1),... IF NEQ IS AN ARRAY
638 C (DIMENSIONED IN DF) AND PAR HAS A LENGTH EXCEEDING NPAR.
639 C SEE THE DESCRIPTIONS OF NEQ AND PAR (BELOW).
641 C NEQ = THE SIZE OF THE ODE SYSTEM (NUMBER OF FIRST ORDER ORDINARY
642 C DIFFERENTIAL EQUATIONS (N) IN THE MODEL). USED ONLY FOR
643 C INPUT. NEQ MAY NOT BE CHANGED DURING THE PROBLEM.
645 C FOR ISOPT = 0, NEQ IS NORMALLY A SCALAR. HOWEVER, NEQ MAY
646 C BE AN ARRAY, WITH NEQ(1) SET TO THE SYSTEM SIZE (N), IN WHICH
647 C CASE THE ODESSA PACKAGE ACCESSES ONLY NEQ(1). HOWEVER,
648 C THIS PARAMETER IS PASSED AS THE NEQ ARGUMENT IN ALL CALLS
649 C TO F, DF, AND JAC. HENCE, IF IT IS AN ARRAY, LOCATIONS
650 C NEQ(2),... MAY BE USED TO STORE OTHER INTEGER DATA AND PASS
651 C IT TO F, DF, AND/OR JAC. FOR ISOPT = 1, NPAR MUST BE LOADED
652 C INTO NEQ(2), AND IS NOT ALLOWED TO CHANGE DURING THE PROBLEM.
653 C IN THESE CASES, SUBROUTINES F, DF, AND/OR JAC MUST INCLUDE
654 C NEQ IN A DIMENSION STATEMENT.
656 C Y = A REAL ARRAY FOR THE VECTOR OF DEPENDENT VARIABLES, OF
657 C DIMENSION (N) BY (NPAR+1). USED FOR BOTH INPUT AND
658 C OUTPUT ON THE FIRST CALL (ISTATE = 1), AND ONLY FOR
659 C OUTPUT ON OTHER CALLS. ON THE FIRST CALL, Y MUST CONTAIN
660 C THE VECTORS OF INITIAL VALUES. ON OUTPUT, Y CONTAINS THE
661 C COMPUTED SOLUTION VECTORS, EVALUATED AT T.
663 C PAR = A REAL ARRAY FOR THE VECTOR OF CONSTANT MODEL PARAMETERS
664 C OF INTEREST IN THE SENSITIVITY ANALYSIS, OF LENGTH NPAR
665 C OR MORE. PAR IS PASSED AS AN ARGUMENT IN ALL CALLS TO F,
666 C DF, AND JAC. HENCE LOCATIONS PAR(NPAR+1),... MAY BE USED
667 C TO STORE OTHER REAL DATA AND PASS IT TO F, DF, AND/OR JAC.
668 C LOCATIONS PAR(1),...,PAR(NPAR) ARE USED AS INPUT ONLY,
669 C AND MUST NOT BE CHANGED DURING THE PROBLEM.
671 C T = THE INDEPENDENT VARIABLE. ON INPUT, T IS USED ONLY ON THE
672 C FIRST CALL, AS THE INITIAL POINT OF THE INTEGRATION.
673 C ON OUTPUT, AFTER EACH CALL, T IS THE VALUE AT WHICH A
674 C COMPUTED SOLUTION Y IS EVALUATED (USUALLY THE SAME AS TOUT).
675 C ON AN ERROR RETURN, T IS THE FARTHEST POINT REACHED.
677 C TOUT = THE NEXT VALUE OF T AT WHICH A COMPUTED SOLUTION IS DESIRED.
678 C USED ONLY FOR INPUT.
680 C WHEN STARTING THE PROBLEM (ISTATE = 1), TOUT MAY BE EQUAL
681 C TO T FOR ONE CALL, THEN SHOULD .NE. T FOR THE NEXT CALL.
682 C FOR THE INITIAL T, AN INPUT VALUE OF TOUT .NE. T IS USED
683 C IN ORDER TO DETERMINE THE DIRECTION OF THE INTEGRATION
684 C (I.E. THE ALGEBRAIC SIGN OF THE STEP SIZES) AND THE ROUGH
685 C SCALE OF THE PROBLEM. INTEGRATION IN EITHER DIRECTION
686 C (FORWARD OR BACKWARD IN T) IS PERMITTED.
688 C IF ITASK = 2 OR 5 (ONE-STEP MODES), TOUT IS IGNORED AFTER
689 C THE FIRST CALL (I.E. THE FIRST CALL WITH TOUT .NE. T).
690 C OTHERWISE, TOUT IS REQUIRED ON EVERY CALL.
692 C IF ITASK = 1, 3, OR 4, THE VALUES OF TOUT NEED NOT BE
693 C MONOTONE, BUT A VALUE OF TOUT WHICH BACKS UP IS LIMITED
694 C TO THE CURRENT INTERNAL T INTERVAL, WHOSE ENDPOINTS ARE
695 C TCUR - HU AND TCUR (SEE OPTIONAL OUTPUTS, BELOW, FOR
696 C TCUR AND HU).
698 C ITOL = AN INDICATOR FOR THE TYPE OF ERROR CONTROL. SEE
699 C DESCRIPTION BELOW UNDER ATOL. USED ONLY FOR INPUT.
701 C RTOL = A RELATIVE ERROR TOLERANCE PARAMETER, EITHER A SCALAR OR
702 C AN ARRAY OF SPACE (N) BY (NPAR+1). SEE DESCRIPTION BELOW
703 C UNDER ATOL. INPUT ONLY.
705 C ATOL = AN ABSOLUTE ERROR TOLERANCE PARAMETER, EITHER A SCALAR OR
706 C AN ARRAY OF SPACE (N) BY (NPAR+1). INPUT ONLY.
708 C THE INPUT PARAMETERS ITOL, RTOL, AND ATOL DETERMINE
709 C THE ERROR CONTROL PERFORMED BY THE KppSolveR. THE KppSolveR WILL
710 C CONTROL THE VECTOR E = (E(I,J)) OF ESTIMATED LOCAL ERRORS
711 C IN Y, ACCORDING TO AN INEQUALITY OF THE FORM
712 C RMS-NORM OF ( E(I,J)/EWT(I,J) ) .LE. 1,
713 C WHERE EWT(I,J) = RTOL(I,J)*ABS(Y(I,J)) + ATOL(I,J),
714 C AND THE RMS-NORM (ROOT-MEAN-SQUARE NORM) HERE IS
715 C RMS-NORM(V) = SQRT ( (1/N) * SUM (V(I,J)**2) ); I =1,...,N.
716 C HERE EWT = (EWT(I,J)) IS A VECTOR OF WEIGHTS WHICH MUST
717 C ALWAYS BE POSITIVE, AND THE VALUES OF RTOL AND ATOL SHOULD
718 C ALL BE NON-NEGATIVE. THE FOLLOWING TABLE GIVES THE TYPES
719 C (SCALAR/ARRAY) OF RTOL AND ATOL, AND THE CORRESPONDING FORM
720 C OF EWT(I,J).
722 C ITOL RTOL ATOL EWT(I,J)
723 C 1 SCALAR SCALAR RTOL*ABS(Y(I,J)) + ATOL
724 C 2 SCALAR ARRAY RTOL*ABS(Y(I,J)) + ATOL(I,J)
725 C 3 ARRAY SCALAR RTOL(I,J)*ABS(Y(I,J)) + ATOL
726 C 4 ARRAY ARRAY RTOL(I,J)*ABS(Y(I,J)) + ATOL(I,J)
728 C WHEN EITHER OF THESE PARAMETERS IS A SCALAR, IT NEED NOT
729 C BE DIMENSIONED IN THE USER-S CALLING PROGRAM.
731 C THE TOTAL NUMBER OF ERROR TEST FAILURES DUE TO THE SENSITIVITY
732 C ANALYSIS, AND WHICH REQUIRE AN INTEGRATION STEP TO BE
733 C REPEATED, ARE ACCUMULATED IN THE LAST NPAR+1 LOCATIONS OF THE
734 C INTEGER WORK ARRAY IWORK (SEE OPTIONAL OUTPUTS BELOW).
735 C THIS INFORMATION MAY BE OF VALUE IN DETERMINING APPROPRIATE
736 C ERROR TOLERANCES TO BE APPLIED TO THE SENSITIVITY FUNCTIONS.
738 C IF NONE OF THE ABOVE CHOICES (WITH ITOL, RTOL, AND ATOL
739 C FIXED THROUGHOUT THE PROBLEM) IS SUITABLE, MORE GENERAL
740 C ERROR CONTROLS CAN BE OBTAINED BY SUBSTITUTING
741 C USER-SUPPLIED ROUTINES FOR THE SETTING OF EWT AND/OR FOR
742 C THE NORM CALCULATION. SEE PART IV BELOW.
744 C IF GLOBAL ERRORS ARE TO BE ESTIMATED BY MAKING A REPEATED
745 C RUN ON THE SAME PROBLEM WITH SMALLER TOLERANCES, THEN ALL
746 C COMPONENTS OF RTOL AND ATOL (I.E. OF EWT) SHOULD BE SCALED
747 C DOWN UNIFORMLY.
749 C ITASK = AN INDEX SPECIFYING THE TASK TO BE PERFORMED.
750 C INPUT ONLY. ITASK HAS THE FOLLOWING VALUES AND MEANINGS.
751 C 1 MEANS NORMAL COMPUTATION OF OUTPUT VALUES OF Y(T) AT
752 C T = TOUT (BY OVERSHOOTING AND INTERPOLATING).
753 C 2 MEANS TAKE ONE STEP ONLY AND RETURN.
754 C 3 MEANS STOP AT THE FIRST INTERNAL MESH POINT AT OR
755 C BEYOND T = TOUT AND RETURN.
756 C 4 MEANS NORMAL COMPUTATION OF OUTPUT VALUES OF Y(T) AT
757 C T = TOUT BUT WITHOUT OVERSHOOTING T = TCRIT.
758 C TCRIT MUST BE INPUT AS RWORK(1). TCRIT MAY BE EQUAL TO
759 C OR BEYOND TOUT, BUT NOT BEHIND IT IN THE DIRECTION OF
760 C INTEGRATION. THIS OPTION IS USEFUL IF THE PROBLEM
761 C HAS A SINGULARITY AT OR BEYOND T = TCRIT.
762 C 5 MEANS TAKE ONE STEP, WITHOUT PASSING TCRIT, AND RETURN.
763 C TCRIT MUST BE INPUT AS RWORK(1).
765 C NOTE.. IF ITASK = 4 OR 5 AND THE KppSolveR REACHES TCRIT
766 C (WITHIN ROUNDOFF), IT WILL RETURN T = TCRIT (EXACTLY) TO
767 C INDICATE THIS (UNLESS ITASK = 4 AND TOUT COMES BEFORE TCRIT,
768 C IN WHICH CASE ANSWERS AT T = TOUT ARE RETURNED FIRST).
770 C ISTATE = AN INDEX USED FOR INPUT AND OUTPUT TO SPECIFY THE
771 C THE STATE OF THE CALCULATION.
773 C ON INPUT, THE VALUES OF ISTATE ARE AS FOLLOWS.
774 C 1 MEANS THIS IS THE FIRST CALL FOR THE PROBLEM
775 C (INITIALIZATIONS WILL BE DONE). SEE NOTE BELOW.
776 C 2 MEANS THIS IS NOT THE FIRST CALL, AND THE CALCULATION
777 C IS TO CONTINUE NORMALLY, WITH NO CHANGE IN ANY INPUT
778 C PARAMETERS EXCEPT POSSIBLY TOUT AND ITASK.
779 C (IF ITOL, RTOL, AND/OR ATOL ARE CHANGED BETWEEN CALLS
780 C WITH ISTATE = 2, THE NEW VALUES WILL BE USED BUT NOT
781 C TESTED FOR LEGALITY.)
782 C 3 MEANS THIS IS NOT THE FIRST CALL, AND THE
783 C CALCULATION IS TO CONTINUE NORMALLY, BUT WITH
784 C A CHANGE IN INPUT PARAMETERS OTHER THAN
785 C TOUT AND ITASK. CHANGES ARE ALLOWED IN
786 C ITOL, RTOL, ATOL, IOPT, LRW, LIW, MF, ML, MU,
787 C AND ANY OF THE OPTIONAL INPUTS EXCEPT H0.
788 C (SEE IWORK DESCRIPTION FOR ML AND MU.)
789 C NOTE.. A PRELIMINARY CALL WITH TOUT = T IS NOT COUNTED
790 C AS A FIRST CALL HERE, AS NO INITIALIZATION OR CHECKING OF
791 C INPUT IS DONE. (SUCH A CALL IS SOMETIMES USEFUL FOR THE
792 C PURPOSE OF OUTPUTTING THE INITIAL CONDITIONS.)
793 C THUS THE FIRST CALL FOR WHICH TOUT .NE. T REQUIRES
794 C ISTATE = 1 ON INPUT.
796 C ON OUTPUT, ISTATE HAS THE FOLLOWING VALUES AND MEANINGS.
797 C 1 MEANS NOTHING WAS DONE, AS TOUT WAS EQUAL TO T WITH
798 C ISTATE = 1 ON INPUT. (HOWEVER, AN INTERNAL COUNTER WAS
799 C SET TO DETECT AND PREVENT REPEATED CALLS OF THIS TYPE.)
800 C 2 MEANS THE INTEGRATION WAS PERFORMED SUCCESSFULLY.
801 C -1 MEANS AN EXCESSIVE AMOUNT OF WORK (MORE THAN MXSTEP
802 C STEPS) WAS DONE ON THIS CALL, BEFORE COMPLETING THE
803 C REQUESTED TASK, BUT THE INTEGRATION WAS OTHERWISE
804 C SUCCESSFUL AS FAR AS T. (MXSTEP IS AN OPTIONAL INPUT
805 C AND IS NORMALLY 500.) TO CONTINUE, THE USER MAY
806 C SIMPLY RESET ISTATE TO A VALUE .GT. 1 AND CALL AGAIN
807 C (THE EXCESS WORK STEP COUNTER WILL BE RESET TO 0).
808 C IN ADDITION, THE USER MAY INCREASE MXSTEP TO AVOID
809 C THIS ERROR RETURN (SEE BELOW ON OPTIONAL INPUTS).
810 C -2 MEANS TOO MUCH ACCURACY WAS REQUESTED FOR THE PRECISION
811 C OF THE MACHINE BEING USED. THIS WAS DETECTED BEFORE
812 C COMPLETING THE REQUESTED TASK, BUT THE INTEGRATION
813 C WAS SUCCESSFUL AS FAR AS T. TO CONTINUE, THE TOLERANCE
814 C PARAMETERS MUST BE RESET, AND ISTATE MUST BE SET
815 C TO 3. THE OPTIONAL OUTPUT TOLSF MAY BE USED FOR THIS
816 C PURPOSE. (NOTE.. IF THIS CONDITION IS DETECTED BEFORE
817 C TAKING ANY STEPS, THEN AN ILLEGAL INPUT RETURN
818 C (ISTATE = -3) OCCURS INSTEAD.)
819 C -3 MEANS ILLEGAL INPUT WAS DETECTED, BEFORE TAKING ANY
820 C INTEGRATION STEPS. SEE WRITTEN MESSAGE FOR DETAILS.
821 C NOTE.. IF THE KppSolveR DETECTS AN INFINITE LOOP OF CALLS
822 C TO THE KppSolveR WITH ILLEGAL INPUT, IT WILL CAUSE
823 C THE RUN TO STOP.
824 C -4 MEANS THERE WERE REPEATED ERROR TEST FAILURES ON
825 C ONE ATTEMPTED STEP, BEFORE COMPLETING THE REQUESTED
826 C TASK, BUT THE INTEGRATION WAS SUCCESSFUL AS FAR AS T.
827 C THE PROBLEM MAY HAVE A SINGULARITY, OR THE INPUT
828 C MAY BE INAPPROPRIATE.
829 C -5 MEANS THERE WERE REPEATED CONVERGENCE TEST FAILURES ON
830 C ONE ATTEMPTED STEP, BEFORE COMPLETING THE REQUESTED
831 C TASK, BUT THE INTEGRATION WAS SUCCESSFUL AS FAR AS T.
832 C THIS MAY BE CAUSED BY AN INACCURATE JACOBIAN MATRIX,
833 C IF ONE IS BEING USED.
834 C -6 MEANS EWT(I,J) BECAME ZERO FOR SOME I,J DURING THE
835 C INTEGRATION. PURE RELATIVE ERROR CONTROL (ATOL(I,J)=0.0)
836 C WAS REQUESTED ON A VARIABLE WHICH HAS NOW VANISHED.
837 C THE INTEGRATION WAS SUCCESSFUL AS FAR AS T.
839 C NOTE.. SINCE THE NORMAL OUTPUT VALUE OF ISTATE IS 2,
840 C IT DOES NOT NEED TO BE RESET FOR NORMAL CONTINUATION.
841 C ALSO, SINCE A NEGATIVE INPUT VALUE OF ISTATE WILL BE
842 C REGARDED AS ILLEGAL, A NEGATIVE OUTPUT VALUE REQUIRES THE
843 C USER TO CHANGE IT, AND POSSIBLY OTHER INPUTS, BEFORE
844 C CALLING THE KppSolveR AGAIN.
846 C IOPT = AN INTEGER ARRAY FLAG TO SPECIFY WHETHER OR NOT ANY OPTIONAL
847 C INPUTS ARE BEING USED ON THIS CALL. INPUT ONLY.
848 C THE OPTIONAL INPUTS ARE LISTED SEPARATELY BELOW.
849 C IOPT(1) = 0 MEANS NO OPTIONAL INPUTS FOR THE KppSolveR WILL BE
850 C USED. DEFAULT VALUES WILL BE USED IN ALL CASES.
851 C = 1 MEANS ONE OR MORE OPTIONAL INPUTS FOR THE
852 C KppSolveR ARE BEING USED.
853 C NOTE : IOPT(1) IS INDEPENDENT OF ISOPT AND IDF.
854 C IOPT(2) = 0 MEANS NO SENSITIVITY ANALYSIS WILL BE PERFORMED.
855 C = 1 MEANS A SENSITIVITY ANALYSIS WILL BE PERFORMED.
856 C NOTE : IOPT(2) IS RENAMED TO ISOPT IN ODESSA.
857 C = 0 MEANS DF/DP WILL BE CALCULATED BY FINITE
858 C DIFFERENCE WITHIN ODESSA.
859 C IOPT(3) = 1 MEANS DF/DP WILL BE CALCULATED BY A USER-SUPPLIED
860 C ROUTINE.
861 C NOTE : IOPT(3) IS RENAMED TO IDF IN ODESSA.
862 C IF IDF = 1, THE USER MUST SUPPLY A
863 C SUBROUTINE DF (THE NAME IS ARBITRARY) AS
864 C DESCRIBED BELOW UNDER DF. FOR IDF = 0,
865 C A DUMMY ARGUMENT CAN BE USED.
867 C RWORK = A REAL WORKING ARRAY (DOUBLE PRECISION).
868 C FOR ISOPT = 0, THE LENGTH OF RWORK MUST BE AT LEAST..
869 C 20 + NYH*(MAXORD + 1) + 3*NEQ + LWM
870 C FOR ISOPT = 1, THE LENGTH OF RWORK MUST BE AT LEAST..
871 C 20 + NYH*(MAXORD + 1) + 2*NYH + LWM + N
872 C WHERE..
873 C NYH = THE TOTAL NUMBER OF DEPENDENT VARIABLES;
874 C (= N IF ISOPT = 0, AND N*(NPAR+1) IF ISOPT = 1).
875 C MAXORD = 12 (IF METH = 1) OR 5 (IF METH = 2) (UNLESS A
876 C SMALLER VALUE IS GIVEN AS AN OPTIONAL INPUT),
877 C LWM = 0 IF MITER = 0,
878 C LWM = N**2 + 2 IF MITER IS 1 OR 2,
879 C LWM = N + 2 IF MITER = 3, AND
880 C LWM = (2*ML+MU+1)*N + 2 IF MITER IS 4 OR 5.
881 C (SEE THE MF DESCRIPTION FOR METH AND MITER.)
883 C THE FIRST 20 WORDS OF RWORK ARE RESERVED FOR CONDITIONAL
884 C AND OPTIONAL INPUTS AND OPTIONAL OUTPUTS.
886 C THE FOLLOWING WORD IN RWORK IS A CONDITIONAL INPUT..
887 C RWORK(1) = TCRIT = CRITICAL VALUE OF T WHICH THE KppSolveR
888 C IS NOT TO OVERSHOOT. REQUIRED IF ITASK IS
889 C 4 OR 5, AND IGNORED OTHERWISE. (SEE ITASK.)
891 C LRW = THE LENGTH OF THE ARRAY RWORK, AS DECLARED BY THE USER.
892 C (THIS WILL BE CHECKED BY THE KppSolveR.)
894 C IWORK = AN INTEGER WORK ARRAY. THE LENGTH MUST BE AT LEAST..
895 C 20 IF MITER = 0 OR 3 (MF = 10, 13, 20, 23), OR
896 C 20 + N OTHERWISE (MF = 11, 12, 14, 15, 21, 22, 24, 25).
897 C FOR ISOPT = 0, OR..
898 C 21 + N + NPAR
899 C FOR ISOPT = 1.
900 C THE FIRST FEW WORDS OF IWORK ARE USED FOR CONDITIONAL AND
901 C OPTIONAL INPUTS AND OPTIONAL OUTPUTS.
903 C THE FOLLOWING 2 WORDS IN IWORK ARE CONDITIONAL INPUTS..
904 C IWORK(1) = ML THESE ARE THE LOWER AND UPPER
905 C IWORK(2) = MU HALF-BANDWIDTHS, RESPECTIVELY, OF THE
906 C BANDED JACOBIAN, EXCLUDING THE MAIN DIAGONAL.
907 C THE BAND IS DEFINED BY THE MATRIX LOCATIONS
908 C (I,J) WITH I-ML .LE. J .LE. I+MU. ML AND MU
909 C MUST SATISFY 0 .LE. ML,MU .LE. NEQ-1.
910 C THESE ARE REQUIRED IF MITER IS 4 OR 5, AND
911 C IGNORED OTHERWISE. ML AND MU MAY IN FACT BE
912 C THE BAND PARAMETERS FOR A MATRIX TO WHICH
913 C DF/DY IS ONLY APPROXIMATELY EQUAL.
916 C LIW = THE LENGTH OF THE ARRAY IWORK, AS DECLARED BY THE USER.
917 C (THIS WILL BE CHECKED BY THE KppSolveR.)
919 C NOTE.. THE WORK ARRAYS MUST NOT BE ALTERED BETWEEN CALLS TO ODESSA
920 C FOR THE SAME PROBLEM, EXCEPT POSSIBLY FOR THE CONDITIONAL AND
921 C OPTIONAL INPUTS, AND EXCEPT FOR THE LAST 2*NYH + N WORDS OF RWORK.
922 C THE LATTER SPACE IS USED FOR INTERNAL SCRATCH SPACE, AND SO IS
923 C AVAILABLE FOR USE BY THE USER OUTSIDE ODESSA BETWEEN CALLS, IF
924 C DESIRED (BUT NOT FOR USE BY F, DF, OR JAC).
926 C JAC = THE NAME OF THE USER-SUPPLIED ROUTINE (MITER = 1 OR 4) TO
927 C COMPUTE THE JACOBIAN MATRIX, DF/DY, AS A FUNCTION OF THE
928 C SCALAR T AND THE VECTORS Y, AND PAR. IT IS TO HAVE THE FORM
929 C SUBROUTINE JAC (NEQ, T, Y, PAR, ML, MU, PD, NROWPD)
930 C DOUBLE PRECISION T, Y, PAR, PD
931 C DIMENSION Y(1), PAR(1), PD(NROWPD,1)
932 C WHERE NEQ, T, Y, PAR, ML, MU, AND NROWPD ARE INPUT AND THE
933 C ARRAY PD IS TO BE LOADED WITH PARTIAL DERIVATIVES (ELEMENTS
934 C OF THE JACOBIAN MATRIX) ON OUTPUT. PD MUST BE GIVEN A FIRST
935 C DIMENSION OF NROWPD. T, Y, AND PAR HAVE THE SAME MEANING AS
936 C IN SUBROUTINE F. (IN THE DIMENSION STATEMENT ABOVE, 1 IS A
937 C DUMMY DIMENSION.. IT CAN BE REPLACED BY ANY VALUE.)
938 C IN THE FULL MATRIX CASE (MITER = 1), ML AND MU ARE
939 C IGNORED, AND THE JACOBIAN IS TO BE LOADED INTO PD IN
940 C COLUMNWISE MANNER, WITH DF(I)/DY(J) LOADED INTO PD(I,J).
941 C IN THE BAND MATRIX CASE (MITER = 4), THE ELEMENTS
942 C WITHIN THE BAND ARE TO BE LOADED INTO PD IN COLUMNWISE
943 C MANNER, WITH DIAGONAL LINES OF DF/DY LOADED INTO THE ROWS
944 C OF PD. THUS DF(I)/DY(J) IS TO BE LOADED INTO PD(I-J+MU+1,J).
945 C ML AND MU ARE THE HALF-BANDWIDTH PARAMETERS (SEE IWORK).
946 C THE LOCATIONS IN PD IN THE TWO TRIANGULAR AREAS WHICH
947 C CORRESPOND TO NONEXISTENT MATRIX ELEMENTS CAN BE IGNORED
948 C OR LOADED ARBITRARILY, AS THEY ARE OVERWRITTEN BY ODESSA.
949 C PD IS PRESET TO ZERO BY THE KppSolveR, SO THAT ONLY THE
950 C NONZERO ELEMENTS NEED BE LOADED BY JAC. EACH CALL TO JAC IS
951 C PRECEDED BY A CALL TO F WITH THE SAME ARGUMENTS NEQ, T, Y,
952 C AND PAR. THUS TO GAIN SOME EFFICIENCY, INTERMEDIATE
953 C QUANTITIES SHARED BY BOTH CALCULATIONS MAY BE SAVED IN A
954 C USER COMMON BLOCK BY F AND NOT RECOMPUTED BY JAC, IF
955 C DESIRED. ALSO, JAC MAY ALTER THE Y ARRAY, IF DESIRED.
956 C JAC MUST BE DECLARED EXTERNAL IN THE CALLING PROGRAM.
957 C SUBROUTINE JAC MAY ACCESS USER-DEFINED QUANTITIES IN
958 C NEQ(2),... AND PAR(NPAR+1),.... SEE THE DESCRIPTIONS OF
959 C NEQ (ABOVE) AND PAR (BELOW).
961 C MF = THE METHOD FLAG. USED ONLY FOR INPUT. THE LEGAL VALUES OF
962 C MF ARE 10, 11, 12, 13, 14, 15, 20, 21, 22, 23, 24, AND 25.
963 C MF HAS DECIMAL DIGITS METH AND MITER.. MF = 10*METH + MITER.
964 C METH INDICATES THE BASIC LINEAR MULTISTEP METHOD..
965 C METH = 1 MEANS THE IMPLICIT ADAMS METHOD.
967 C METH = 2 MEANS THE METHOD BASED ON BACKWARD
968 C DIFFERENTIATION FORMULAS (BDF-S).
969 C MITER INDICATES THE CORRECTOR ITERATION METHOD..
970 C MITER = 0 MEANS FUNCTIONAL ITERATION (NO JACOBIAN MATRIX
971 C IS INVOLVED).
972 C MITER = 1 MEANS CHORD ITERATION WITH A USER-SUPPLIED
973 C FULL (NEQ BY NEQ) JACOBIAN.
974 C MITER = 2 MEANS CHORD ITERATION WITH AN INTERNALLY
975 C GENERATED (DIFFERENCE QUOTIENT) FULL JACOBIAN
976 C (USING NEQ EXTRA CALLS TO F PER DF/DY VALUE).
977 C MITER = 3 MEANS CHORD ITERATION WITH AN INTERNALLY
978 C GENERATED DIAGONAL JACOBIAN APPROXIMATION.
979 C (USING 1 EXTRA CALL TO F PER DF/DY EVALUATION).
980 C MITER = 4 MEANS CHORD ITERATION WITH A USER-SUPPLIED
981 C BANDED JACOBIAN.
982 C MITER = 5 MEANS CHORD ITERATION WITH AN INTERNALLY
983 C GENERATED BANDED JACOBIAN (USING ML+MU+1 EXTRA
984 C CALLS TO F PER DF/DY EVALUATION).
985 C IF MITER = 1 OR 4, THE USER MUST SUPPLY A SUBROUTINE JAC
986 C (THE NAME IS ARBITRARY) AS DESCRIBED ABOVE UNDER JAC.
987 C FOR OTHER VALUES OF MITER, A DUMMY ARGUMENT CAN BE USED.
989 C IF A SENSITIVITY ANLYSIS IS DESIRED (ISOPT = 1), MITER = 0
990 C AND 3 ARE DISALLOWED. IN THESE CASES, THE USER IS RECOMMENDED
991 C TO SUPPLY AN ANALYTICAL JACOBIAN (MITER = 1 OR 4) AND AN
992 C ANALYTICAL INHOMOGENEITY MATRIX (IDF = 1).
993 C----------------------------------------------------------------------
994 C OPTIONAL INPUTS.
996 C THE FOLLOWING IS A LIST OF THE OPTIONAL INPUTS PROVIDED FOR IN THE
997 C CALL SEQUENCE. (SEE ALSO PART II.) FOR EACH SUCH INPUT VARIABLE,
998 C THIS TABLE LISTS ITS NAME AS USED IN THIS DOCUMENTATION, ITS
999 C LOCATION IN THE CALL SEQUENCE, ITS MEANING, AND THE DEFAULT VALUE.
1000 C THE USE OF ANY OF THESE INPUTS REQUIRES IOPT(1) = 1, AND IN THAT
1001 C CASE ALL OF THESE INPUTS ARE EXAMINED. A VALUE OF ZERO FOR ANY
1002 C OF THESE OPTIONAL INPUTS WILL CAUSE THE DEFAULT VALUE TO BE USED.
1003 C THUS TO USE A SUBSET OF THE OPTIONAL INPUTS, SIMPLY PRELOAD
1004 C LOCATIONS 5 TO 10 IN RWORK AND IWORK TO 0.0 AND 0 RESPECTIVELY, AND
1005 C THEN SET THOSE OF INTEREST TO NONZERO VALUES.
1007 C NAME LOCATION MEANING AND DEFAULT VALUE
1009 C H0 RWORK(5) THE STEP SIZE TO BE ATTEMPTED ON THE FIRST STEP.
1010 C THE DEFAULT VALUE IS DETERMINED BY THE KppSolveR.
1012 C HMAX RWORK(6) THE MAXIMUM ABSOLUTE STEP SIZE ALLOWED.
1013 C THE DEFAULT VALUE IS INFINITE.
1015 C HMIN RWORK(7) THE MINIMUM ABSOLUTE STEP SIZE ALLOWED.
1016 C THE DEFAULT VALUE IS 0. (THIS LOWER BOUND IS NOT
1017 C ENFORCED ON THE FINAL STEP BEFORE REACHING TCRIT
1018 C WHEN ITASK = 4 OR 5.)
1020 C MAXORD IWORK(5) THE MAXIMUM ORDER TO BE ALLOWED. THE DEFAULT
1021 C VALUE IS 12 IF METH = 1, AND 5 IF METH = 2.
1022 C IF MAXORD EXCEEDS THE DEFAULT VALUE, IT WILL
1023 C BE REDUCED TO THE DEFAULT VALUE.
1024 C IF MAXORD IS CHANGED DURING THE PROBLEM, IT MAY
1025 C CAUSE THE CURRENT ORDER TO BE REDUCED.
1027 C MXSTEP IWORK(6) MAXIMUM NUMBER OF (INTERNALLY DEFINED) STEPS
1028 C ALLOWED DURING ONE CALL TO THE KppSolveR.
1029 C THE DEFAULT VALUE IS 500.
1031 C MXHNIL IWORK(7) MAXIMUM NUMBER OF MESSAGES PRINTED (PER PROBLEM)
1032 C WARNING THAT T + H = T ON A STEP (H = STEP SIZE).
1033 C THIS MUST BE POSITIVE TO RESULT IN A NON-DEFAULT
1034 C VALUE. THE DEFAULT VALUE IS 10.
1035 C----------------------------------------------------------------------
1036 C OPTIONAL OUTPUTS.
1038 C AS OPTIONAL ADDITIONAL OUTPUT FROM ODESSA, THE VARIABLES LISTED
1039 C BELOW ARE QUANTITIES RELATED TO THE PERFORMANCE OF ODESSA
1040 C WHICH ARE AVAILABLE TO THE USER. THESE ARE COMMUNICATED BY WAY OF
1041 C THE WORK ARRAYS, BUT ALSO HAVE INTERNAL MNEMONIC NAMES AS SHOWN.
1042 C EXCEPT WHERE STATED OTHERWISE, ALL OF THESE OUTPUTS ARE DEFINED
1043 C ON ANY SUCCESSFUL RETURN FROM ODESSA, AND ON ANY RETURN WITH
1044 C ISTATE = -1, -2, -4, -5, OR -6. ON AN ILLEGAL INPUT RETURN
1045 C (ISTATE = -3), THEY WILL BE UNCHANGED FROM THEIR EXISTING VALUES
1046 C (IF ANY), EXCEPT POSSIBLY FOR TOLSF, LENRW, AND LENIW.
1047 C ON ANY ERROR RETURN, OUTPUTS RELEVANT TO THE ERROR WILL BE DEFINED,
1048 C AS NOTED BELOW.
1050 C NAME LOCATION MEANING
1052 C HU RWORK(11) THE STEP SIZE IN T LAST USED (SUCCESSFULLY).
1054 C HCUR RWORK(12) THE STEP SIZE TO BE ATTEMPTED ON THE NEXT STEP.
1056 C TCUR RWORK(13) THE CURRENT VALUE OF THE INDEPENDENT VARIABLE
1057 C WHICH THE KppSolveR HAS ACTUALLY REACHED, I.E. THE
1058 C CURRENT INTERNAL MESH POINT IN T. ON OUTPUT, TCUR
1059 C WILL ALWAYS BE AT LEAST AS FAR AS THE ARGUMENT
1060 C T, BUT MAY BE FARTHER (IF INTERPOLATION WAS DONE).
1062 C TOLSF RWORK(14) A TOLERANCE SCALE FACTOR, GREATER THAN 1.0,
1063 C COMPUTED WHEN A REQUEST FOR TOO MUCH ACCURACY WAS
1064 C DETECTED (ISTATE = -3 IF DETECTED AT THE START OF
1065 C THE PROBLEM, ISTATE = -2 OTHERWISE). IF ITOL IS
1066 C LEFT UNALTERED BUT RTOL AND ATOL ARE UNIFORMLY
1067 C SCALED UP BY A FACTOR OF TOLSF FOR THE NEXT CALL,
1068 C THEN THE KppSolveR IS DEEMED LIKELY TO SUCCEED.
1069 C (THE USER MAY ALSO IGNORE TOLSF AND ALTER THE
1070 C TOLERANCE PARAMETERS IN ANY OTHER WAY APPROPRIATE.)
1072 C NST IWORK(11) THE NUMBER OF STEPS TAKEN FOR THE PROBLEM SO FAR.
1074 C NFE IWORK(12) THE NUMBER OF F EVALUATIONS FOR THE PROBLEM SO FAR.
1076 C NJE IWORK(13) THE NUMBER OF JACOBIAN EVALUATIONS (AND OF MATRIX
1077 C LU DECOMPOSITIONS IF ISOPT = 0) FOR THE PROBLEM SO
1078 C FAR. IF ISOPT = 1, THE NUMBER OF LU DECOMPOSITIONS
1079 C IS EQUAL TO NJE - NSPE (SEE BELOW).
1081 C NQU IWORK(14) THE METHOD ORDER LAST USED (SUCCESSFULLY).
1083 C NQCUR IWORK(15) THE ORDER TO BE ATTEMPTED ON THE NEXT STEP.
1085 C IMXER IWORK(16) THE INDEX OF THE COMPONENT OF LARGEST MAGNITUDE IN
1086 C THE WEIGHTED LOCAL ERROR VECTOR (E(I,J)/EWT(I,J)),
1087 C ON AN ERROR RETURN WITH ISTATE = -4 OR -5.
1089 C LENRW IWORK(17) THE LENGTH OF RWORK ACTUALLY REQUIRED.
1090 C THIS IS DEFINED ON NORMAL RETURNS AND ON AN ILLEGAL
1091 C INPUT RETURN FOR INSUFFICIENT STORAGE.
1093 C LENIW IWORK(18) THE LENGTH OF IWORK ACTUALLY REQUIRED.
1094 C THIS IS DEFINED ON NORMAL RETURNS AND ON AN ILLEGAL
1095 C INPUT RETURN FOR INSUFFICIENT STORAGE.
1097 C NDFE IWORK(19) THE NUMBER OF DF/DP (VECTOR) EVALUATIONS.
1099 C NSPE IWORK(20) THE NUMBER OF CALLS TO SUBROUTINE SPRIME. EACH CALL
1100 C TO SPRIME REQUIRES A JACOBIAN EVALUATION, BUT NOT
1101 C AN LU DECOMPOSITION.
1103 C THE FOLLOWING ARRAYS ARE SEGMENTS OF THE RWORK AND IWORK ARRAYS
1104 C WHICH MAY ALSO BE OF INTEREST TO THE USER AS OPTIONAL OUTPUTS.
1105 C FOR EACH ARRAY, THE TABLE BELOW GIVES ITS INTERNAL NAME, ITS BASE
1106 C ADDRESS IN RWORK OR IWORK, AND ITS DESCRIPTION.
1108 C NAME BASE ADDRESS DESCRIPTION
1110 C YH 21 IN RWORK THE NORDSIECK HISTORY ARRAY, OF SIZE NYH BY
1111 C (NQCUR + 1). FOR J = 0,1,...,NQCUR, COLUMN J+1
1112 C OF YH CONTAINS HCUR**J/FACTORIAL(J) TIMES
1113 C THE J-TH DERIVATIVE OF THE INTERPOLATING
1114 C POLYNOMIAL CURRENTLY REPRESENTING THE SOLUTION,
1115 C EVALUATED AT T = TCUR.
1117 C ACOR LENRW-NYH+1 ARRAY OF SIZE NYH USED FOR THE ACCUMULATED
1118 C IN RWORK CORRECTIONS ON EACH STEP, SCALED ON OUTPUT
1119 C TO REPRESENT THE ESTIMATED LOCAL ERROR IN Y
1120 C ON THE LAST STEP. THIS IS THE VECTOR E IN
1121 C THE DESCRIPTION OF THE ERROR CONTROL.
1122 C IT IS DEFINED ONLY ON A SUCCESSFUL RETURN
1123 C FROM ODESSA.
1124 C NRS LENIW-NPAR ARRAY OF SIZE NPAR+1, USED TO STORE THE
1125 C IN IWORK ACCUMULATED NUMBER OF REPEATED STEPS DUE TO
1126 C THE SENSITIVITY ANALYSIS..
1127 C NRS(1) = TOTAL NUMBER OF REPEATED STEPS,
1128 C NRS(2),... = NUMBER OF REPEATED STEPS DUE TO
1129 C MODEL PARAMETER 1,...
1131 C----------------------------------------------------------------------
1132 C PART II. OTHER ROUTINES CALLABLE.
1134 C THE FOLLOWING ARE OPTIONAL CALLS WHICH THE USER MAY MAKE TO
1135 C GAIN ADDITIONAL CAPABILITIES IN CONJUNCTION WITH ODESSA.
1136 C (THE ROUTINES XSETUN AND XSETF ARE DESIGNED TO CONFORM TO THE
1137 C SLATEC ERROR HANDLING PACKAGE.)
1139 C FORM OF CALL FUNCTION
1140 C CALL XSETUN(LUN) SET THE LOGICAL UNIT NUMBER, LUN, FOR
1141 C OUTPUT OF MESSAGES FROM ODESSA, IF
1142 C THE DEFAULT IS NOT DESIRED.
1143 C THE DEFAULT VALUE OF LUN IS 6.
1145 C CALL XSETF(MFLAG) SET A FLAG TO CONTROL THE PRINTING OF
1146 C MESSAGES BY ODESSA..
1147 C MFLAG = 0 MEANS DO NOT PRINT. (DANGER..
1148 C THIS RISKS LOSING VALUABLE INFORMATION.)
1149 C MFLAG = 1 MEANS PRINT (THE DEFAULT).
1151 C EITHER OF THE ABOVE CALLS MAY BE MADE AT
1152 C ANY TIME AND WILL TAKE EFFECT IMMEDIATELY.
1154 C CALL SVCOM (RSAV, ISAV) STORE IN RSAV AND ISAV THE CONTENTS
1155 C OF THE INTERNAL COMMON BLOCKS USED BY
1156 C ODESSA (SEE PART III BELOW).
1157 C RSAV MUST BE A REAL ARRAY OF LENGTH 222
1158 C OR MORE, AND ISAV MUST BE AN INTEGER
1159 C ARRAY OF LENGTH 54 OR MORE.
1161 C CALL RSCOM (RSAV, ISAV) RESTORE, FROM RSAV AND ISAV, THE CONTENTS
1162 C OF THE INTERNAL COMMON BLOCKS USED BY
1163 C ODESSA. PRESUMES A PRIOR CALL TO SVCOM
1164 C WITH THE SAME ARGUMENTS.
1166 C SVCOM AND RSCOM ARE USEFUL IF
1167 C INTERRUPTING A RUN AND RESTARTING
1168 C LATER, OR ALTERNATING BETWEEN TWO OR
1169 C MORE PROBLEMS KppSolveD WITH ODESSA.
1171 C CALL INTDY(,,,,,) PROVIDE DERIVATIVES OF Y, OF VARIOUS
1172 C (SEE BELOW) ORDERS, AT A SPECIFIED POINT T, IF
1173 C DESIRED. IT MAY BE CALLED ONLY AFTER
1174 C A SUCCESSFUL RETURN FROM ODESSA.
1176 C THE DETAILED INSTRUCTIONS FOR USING INTDY ARE AS FOLLOWS.
1177 C THE FORM OF THE CALL IS..
1179 C CALL INTDY (T, K, RWORK(21), NYH, DKY, IFLAG)
1181 C THE INPUT PARAMETERS ARE..
1183 C T = VALUE OF INDEPENDENT VARIABLE WHERE ANSWERS ARE DESIRED
1184 C (NORMALLY THE SAME AS THE T LAST RETURNED BY ODESSA).
1185 C FOR VALID RESULTS, T MUST LIE BETWEEN TCUR - HU AND TCUR.
1186 C (SEE OPTIONAL OUTPUTS FOR TCUR AND HU.)
1187 C K = INTEGER ORDER OF THE DERIVATIVE DESIRED. K MUST SATISFY
1188 C 0 .LE. K .LE. NQCUR, WHERE NQCUR IS THE CURRENT ORDER
1189 C (SEE OPTIONAL OUTPUTS). THE CAPABILITY CORRESPONDING
1190 C TO K = 0, I.E. COMPUTING Y(T), IS ALREADY PROVIDED
1191 C BY ODESSA DIRECTLY. SINCE NQCUR .GE. 1, THE FIRST
1192 C DERIVATIVE DY/DT IS ALWAYS AVAILABLE WITH INTDY.
1193 C RWORK(21) = THE BASE ADDRESS OF THE HISTORY ARRAY YH.
1194 C NYH = COLUMN LENGTH OF YH, EQUAL TO THE TOTAL NUMBER OF
1195 C DEPENDENT VARIABLES. IF ISOPT = 0, NYH = N. IF ISOPT = 1,
1196 C NYH = N * (NPAR + 1).
1198 C THE OUTPUT PARAMETERS ARE..
1200 C DKY = A REAL ARRAY OF LENGTH NYH CONTAINING THE COMPUTED VALUE
1201 C OF THE K-TH DERIVATIVE OF Y(T).
1202 C IFLAG = INTEGER FLAG, RETURNED AS 0 IF K AND T WERE LEGAL,
1203 C -1 IF K WAS ILLEGAL, AND -2 IF T WAS ILLEGAL.
1204 C ON AN ERROR RETURN, A MESSAGE IS ALSO WRITTEN.
1205 C----------------------------------------------------------------------
1206 C PART III. COMMON BLOCKS.
1208 C IF ODESSA IS TO BE USED IN AN OVERLAY SITUATION, THE USER
1209 C MUST DECLARE, IN THE PRIMARY OVERLAY, THE VARIABLES IN..
1210 C (1) THE CALL SEQUENCE TO ODESSA,
1211 C (2) THE THREE INTERNAL COMMON BLOCKS
1212 C /ODE001/ OF LENGTH 258 (219 DOUBLE PRECISION WORDS
1213 C FOLLOWED BY 39 INTEGER WORDS),
1214 C /ODE002/ OF LENGTH 14 (3 DOUBLE PRECISION WORDS FOLLOWED
1215 C BY 11 INTEGER WORDS),
1216 C /EH0001/ OF LENGTH 2 (INTEGER WORDS).
1218 C IF ODESSA IS USED ON A SYSTEM IN WHICH THE CONTENTS OF INTERNAL
1219 C COMMON BLOCKS ARE NOT PRESERVED BETWEEN CALLS, THE USER SHOULD
1220 C DECLARE THE ABOVE THREE COMMON BLOCKS IN HIS MAIN PROGRAM TO INSURE
1221 C THAT THEIR CONTENTS ARE PRESERVED.
1223 C IF THE SOLUTION OF A GIVEN PROBLEM BY ODESSA IS TO BE INTERRUPTED
1224 C AND THEN LATER CONTINUED, SUCH AS WHEN RESTARTING AN INTERRUPTED RUN
1225 C OR ALTERNATING BETWEEN TWO OR MORE PROBLEMS, THE USER SHOULD SAVE,
1226 C FOLLOWING THE RETURN FROM THE LAST ODESSA CALL PRIOR TO THE
1227 C INTERRUPTION, THE CONTENTS OF THE CALL SEQUENCE VARIABLES AND THE
1228 C INTERNAL COMMON BLOCKS, AND LATER RESTORE THESE VALUES BEFORE THE
1229 C NEXT ODESSA CALL FOR THAT PROBLEM. TO SAVE AND RESTORE THE COMMON
1230 C BLOCKS, USE SUBROUTINES SVCOM AND RSCOM (SEE PART II ABOVE).
1232 C----------------------------------------------------------------------
1233 C PART IV. OPTIONALLY REPLACEABLE KppSolveR ROUTINES.
1235 C BELOW ARE DESCRIPTIONS OF TWO ROUTINES IN THE ODESSA PACKAGE WHICH
1236 C RELATE TO THE MEASUREMENT OF ERRORS. EITHER ROUTINE CAN BE
1237 C REPLACED BY A USER-SUPPLIED VERSION, IF DESIRED. HOWEVER, SINCE SUCH
1238 C A REPLACEMENT MAY HAVE A MAJOR IMPACT ON PERFORMANCE, IT SHOULD BE
1239 C DONE ONLY WHEN ABSOLUTELY NECESSARY, AND ONLY WITH GREAT CAUTION.
1240 C (NOTE.. THE MEANS BY WHICH THE PACKAGE VERSION OF A ROUTINE IS
1241 C SUPERSEDED BY THE USER-S VERSION MAY BE SYSTEM-DEPENDENT.)
1243 C (A) EWSET.
1244 C THE FOLLOWING SUBROUTINE IS CALLED JUST BEFORE EACH INTERNAL
1245 C INTEGRATION STEP, AND SETS THE ARRAY OF ERROR WEIGHTS, EWT, AS
1246 C DESCRIBED UNDER ITOL/RTOL/ATOL ABOVE..
1247 C SUBROUTINE EWSET (NYH, ITOL, RTOL, ATOL, YCUR, EWT)
1248 C WHERE NEQ, ITOL, RTOL, AND ATOL ARE AS IN THE ODESSA CALL SEQUENCE,
1249 C YCUR CONTAINS THE CURRENT DEPENDENT VARIABLE VECTOR, AND
1250 C EWT IS THE ARRAY OF WEIGHTS SET BY EWSET.
1252 C IF THE USER SUPPLIES THIS SUBROUTINE, IT MUST RETURN IN EWT(I)
1253 C (I = 1,...,NYH) A POSITIVE QUANTITY SUITABLE FOR COMPARING ERRORS
1254 C IN Y(I) TO. THE EWT ARRAY RETURNED BY EWSET IS PASSED TO THE
1255 C VNORM ROUTINE (SEE BELOW), AND ALSO USED BY ODESSA IN THE COMPUTATION
1256 C OF THE OPTIONAL OUTPUT IMXER, THE DIAGONAL JACOBIAN APPROXIMATION,
1257 C AND THE INCREMENTS FOR DIFFERENCE QUOTIENT JACOBIANS.
1259 C IN THE USER-SUPPLIED VERSION OF EWSET, IT MAY BE DESIRABLE TO USE
1260 C THE CURRENT VALUES OF DERIVATIVES OF Y. DERIVATIVES UP TO ORDER NQ
1261 C ARE AVAILABLE FROM THE HISTORY ARRAY YH, DESCRIBED ABOVE UNDER
1262 C OPTIONAL OUTPUTS. IN EWSET, YH IS IDENTICAL TO THE YCUR ARRAY,
1263 C EXTENDED TO NQ + 1 COLUMNS WITH A COLUMN LENGTH OF NYH AND SCALE
1264 C FACTORS OF H**J/FACTORIAL(J). ON THE FIRST CALL FOR THE PROBLEM,
1265 C GIVEN BY NST = 0, NQ IS 1 AND H IS TEMPORARILY SET TO 1.0.
1266 C THE QUANTITIES NQ, NYH, H, AND NST CAN BE OBTAINED BY INCLUDING
1267 C IN EWSET THE STATEMENTS..
1268 C DOUBLE PRECISION H, RLS
1269 C COMMON /ODE001/ RLS(219),ILS(39)
1270 C NQ = ILS(35)
1271 C NYH = ILS(14)
1272 C NST = ILS(36)
1273 C H = RLS(213)
1274 C THUS, FOR EXAMPLE, THE CURRENT VALUE OF DY/DT CAN BE OBTAINED AS
1275 C YCUR(NYH+I)/H (I=1,...,N) (AND THE DIVISION BY H IS
1276 C UNNECESSARY WHEN NST = 0).
1278 C (B) VNORM.
1279 C THE FOLLOWING IS A REAL FUNCTION ROUTINE WHICH COMPUTES THE WEIGHTED
1280 C ROOT-MEAN-SQUARE NORM OF A VECTOR V..
1281 C D = VNORM (LV, V, W)
1282 C WHERE..
1283 C LV = THE LENGTH OF THE VECTOR,
1284 C V = REAL ARRAY OF LENGTH N CONTAINING THE VECTOR,
1285 C W = REAL ARRAY OF LENGTH N CONTAINING WEIGHTS,
1286 C D = SQRT( (1/N) * SUM(V(I)*W(I))**2 ).
1287 C VNORM IS CALLED WITH LV = N AND WITH W(I) = 1.0/EWT(I), WHERE
1288 C EWT IS AS SET BY SUBROUTINE EWSET.
1290 C IF THE USER SUPPLIES THIS FUNCTION, IT SHOULD RETURN A NON-NEGATIVE
1291 C VALUE OF VNORM SUITABLE FOR USE IN THE ERROR CONTROL IN ODESSA.
1292 C NONE OF THE ARGUMENTS SHOULD BE ALTERED BY VNORM.
1293 C FOR EXAMPLE, A USER-SUPPLIED VNORM ROUTINE MIGHT..
1294 C -SUBSTITUTE A MAX-NORM OF (V(I)*W(I)) FOR THE RMS-NORM, OR
1295 C -IGNORE SOME COMPONENTS OF V IN THE NORM, WITH THE EFFECT OF
1296 C SUPPRESSING THE ERROR CONTROL ON THOSE COMPONENTS OF Y.
1297 C----------------------------------------------------------------------
1298 C OTHER ROUTINES IN THE ODESSA PACKAGE.
1300 C IN ADDITION TO SUBROUTINE ODESSA, THE ODESSA PACKAGE INCLUDES THE
1301 C FOLLOWING SUBROUTINES AND FUNCTION ROUTINES..
1302 C INTDY COMPUTES AN INTERPOLATED VALUE OF THE Y VECTOR AT T = TOUT.
1303 C STODE IS THE CORE INTEGRATOR, WHICH DOES ONE STEP OF THE
1304 C INTEGRATION AND THE ASSOCIATED ERROR CONTROL.
1305 C STESA MANAGES THE SOLUTION OF THE SENSITIVITY FUNCTIONS.
1306 C CFODE SETS ALL METHOD COEFFICIENTS AND TEST CONSTANTS.
1307 C PREPJ COMPUTES AND PREPROCESSES THE JACOBIAN MATRIX J = DF/DY
1308 C AND THE NEWTON ITERATION MATRIX P = I - H*L0*J.
1309 C IT IS ALSO CALLED BY SPRIME (WITH JOPT = 1) TO JUST
1310 C COMPUTE THE JACOBIAN MATRIX.
1311 C PREPDF COMPUTES THE INHOMOGENEITY MATRIX DF/DP.
1312 C SPRIME DEFINES THE SYSTEM OF SENSITIVITY EQUATIONS.
1313 C SOLSY MANAGES SOLUTION OF LINEAR SYSTEM IN CHORD ITERATION.
1314 C EWSET SETS THE ERROR WEIGHT VECTOR EWT BEFORE EACH STEP.
1315 C VNORM COMPUTES THE WEIGHTED R.M.S. NORM OF A VECTOR.
1316 C SVCOM AND RSCOM ARE USER-CALLABLE ROUTINES TO SAVE AND RESTORE,
1317 C RESPECTIVELY, THE CONTENTS OF THE INTERNAL COMMON BLOCKS.
1318 C DGEFA AND DGESL ARE ROUTINES FROM LINPACK FOR SOLVING FULL
1319 C SYSTEMS OF LINEAR ALGEBRAIC EQUATIONS.
1320 C DGBFA AND DGBSL ARE ROUTINES FROM LINPACK FOR SOLVING BANDED
1321 C LINEAR SYSTEMS.
1322 C DAXPY, DSCAL, IDAMAX, AND DDOT ARE BASIC LINEAR ALGEBRA MODULES
1323 C (BLAS) USED BY THE ABOVE LINPACK ROUTINES.
1324 C D1MACH COMPUTES THE UNIT ROUNDOFF IN A MACHINE-INDEPENDENT MANNER.
1325 C XERR, XSETUN, AND XSETF HANDLE THE PRINTING OF ALL ERROR
1326 C MESSAGES AND WARNINGS.
1327 C NOTE.. VNORM, IDAMAX, DDOT, AND D1MACH ARE FUNCTION ROUTINES.
1328 C ALL THE OTHERS ARE SUBROUTINES.
1330 C THE FORTRAN GENERIC INTRINSIC FUNCTIONS USED BY ODESSA ARE..
1331 C ABS, MAX, MIN, REAL, MOD, SIGN, SQRT, AND WRITE
1333 C A BLOCK DATA SUBPROGRAM IS ALSO INCLUDED WITH THE PACKAGE,
1334 C FOR LOADING SOME OF THE VARIABLES IN INTERNAL COMMON.
1336 C----------------------------------------------------------------------
1337 C PART V. GENERAL REMARKS
1339 C THIS SECTION HIGHLIGHTS THE BASIC DIFFERENCES BETWEEN THE ORIGINAL
1340 C LSODE PACKAGE AND THE ODESSA MODIFICATION. THIS IS PROVIDED AS A
1341 C SERVICE TO EXPERIENCED LSODE USERS TO EXPEDITE FAMILIARIZATION WITH
1342 C ODESSA.
1344 C (A). ORIGINAL SUBROUTINES AND FUNCTIONS.
1346 C OF THE ORIGINAL 22 SUBROUTINES AND FUNCTIONS USED IN THE LSODE
1347 C PACKAGE, ALL ARE USED BY ODESSA, WITH THE FOLLOWING HAVING BEEN
1348 C MODIFIED..
1350 C LSODE THE ORIGINAL DRIVER SUBROUTINE FOR THE LSODE PACKAGE IS
1351 C EXTENSIVELY MODIFIED AND RENAMED ODESSA, WHICH NOW
1352 C CONTAINS A CALL TO SPRIME TO ESTABLISH INITIAL CONDITIONS
1353 C FOR THE SENSITIVITY CALCULATIONS.
1355 C STODE THE ONE STEP INTEGRATOR IS SLIGHTLY MODIFIED AND RETAINS
1356 C ITS ORIGINAL NAME. IT NOW CONTAINS THE CALL TO STESA,
1357 C AND ALSO CALLS SPRIME IF KFLAG .LE. -3.
1359 C PREPJ ALSO NAMED PREPJ IN ODESSA IS SLIGHTLY MODIFIED TO ALLOW
1360 C FOR THE CALCULATION OF JACOBIAN WITH NO PREPROCESSING
1361 C (JOPT = 1).
1363 C (B). NEW SUBROUTINES.
1365 C IN ADDITION TO THE CHANGES NOTED ABOVE, THREE NEW SUBROUTINES
1366 C HAVE BEEN INTRODUCED (SEE STESA, SPRIME, AND PREPDF AS DESCRIBED
1367 C IN PART IV. ABOVE).
1369 C (C). COMMON BLOCKS.
1371 C /LS0001/ RETAINS THE SAME LENGTH AND IS RENAMED /ODE001/;
1372 C HOWEVER THE REAL ARRAY ROWNS(209) IS SHORTENED TO A
1373 C LENGTH OF (173) REAL WORDS, ALLOWING THE REMOVAL OF
1374 C TESCO(3,12) WHICH IS NOW PASSED FROM STODE TO STESA.
1375 C IN ADDITION, THE INTEGER ARRAY IOWNS(6) IS SHORTENED
1376 C TO A LENGTH OF (4) INTEGER WORDS, ALLOWING THE REMOVAL
1377 C OF IALTH AND LMAX WHICH ARE NOW PASSED FROM STODE TO
1378 C STESA.
1380 C /ODE002/ ADDED COMMON BLOCK FOR VARIABLES IMPORTANT TO
1381 C SENSITIVITY ANALYSIS (SEE PART III. ABOVE). A BLOCK
1382 C DATA PROGRAM IS NOT REQUIRED FOR THIS COMMON BLOCK.
1384 C SVCOM,RSCOM THESE TWO SUBROUTINES ARE MODIFIED TO HANDLE
1385 C COMMON BLOCK /ODE002/ AS WELL.
1387 C (D). OPTIONAL INPUTS.
1389 C THE FULL SET OF OPTIONAL INPUTS AVAILABLE IN LSODE IS ALSO
1390 C AVAILABLE IN ODESSA, WITH THE EXCEPTION THAT THE NUMBER OF ODE'S
1391 C IN THE MODEL (NEQ(1)), MAY NOT BE CHANGED DURING THE PROBLEM.
1392 C IN ODESSA, NYH NOW REFERS TO THE TOTAL NUMBER OF FIRST-ORDER
1393 C ODE'S (MODEL AND SENSITIVITY EQUATIONS) WHICH IS EQUAL TO
1394 C NEQ(1) IF ISOPT = 0, OR NEQ(1)*(NEQ(2)+1) IF ISOPT = 1.
1395 C NEQ(1), NEQ(2), AND NYH ARE NOT ALLOWED TO CHANGE DURING
1396 C THE COURSE OF AN INTEGRATION.
1398 C (E). OPTIONAL OUTPUTS.
1400 C THE FULL SET OF OPTIONAL OUTPUTS AVAILABLE IN LSODE IS ALSO
1401 C AVAILABLE IN ODESSA. IN ADDITION, IWORK(19) AND IWORK(20) ARE
1402 C LOADED WITH NDFE AND NSPE, RESPECTIVELY, UPON OUTPUT. THE TOTAL
1403 C NUMBER OF LU DECOMPOSITIONS OF THE PROCESSED JACOBIAN IS EQUAL
1404 C TO NJE - NSPE.
1405 C-----------------------------------------------------------------------
1406 SUBROUTINE ODESSA (F, DF, NEQ, Y, PAR, T, TOUT, ITOL, RTOL, ATOL,
1407 1 ITASK, ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JAC, MF)
1408 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1409 LOGICAL IHIT
1410 EXTERNAL F, DF, JAC, PREPJ, SOLSY, PREPDF
1411 DIMENSION NEQ(*), Y(*), PAR(*), RTOL(*), ATOL(*), IOPT(*),
1412 1 RWORK(LRW), IWORK(LIW), MORD(2)
1413 C-----------------------------------------------------------------------
1414 C THIS IS THE SEPTEMBER 1, 1986 VERSION OF ODESSA..
1415 C AN ORDINARY DIFFERENTIAL EQUATION KppSolveR WITH EXPLICIT SIMULTANEOUS
1416 C SENSITIVITY ANALYSIS.
1418 C THIS PACKAGE IS A MODIFICATION OF THE AUGUST 13, 1981 VERSION OF
1419 C LSODE.. LIVERMORE KppSolveR FOR ORDINARY DIFFERENTIAL EQUATIONS.
1420 C THIS VERSION IS IN DOUBLE PRECISION.
1422 C ODESSA KppSolveS FOR THE FIRST-ORDER SENSITIVITY COEFFICIENTS..
1423 C DY(I)/DP, FOR A SINGLE PARAMETER, OR,
1424 C DY(I)/DP(J), FOR MULTIPLE PARAMETERS,
1425 C ASSOCIATED WITH A SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS..
1426 C DY(T)/DT = F(Y,T;P).
1427 C-----------------------------------------------------------------------
1428 C REFERENCES...
1430 C 1. JORGE R. LEIS AND MARK A. KRAMER, THE SIMULTANEOUS SOLUTION AND
1431 C EXPLICIT SENSITIVITY ANALYSIS OF SYSTEMS DESCRIBED BY ORDINARY
1432 C DIFFERENTIAL EQUATIONS, SUBMITTED TO ACM TRANS. MATH. SOFTWARE,
1433 C (1985).
1435 C 2. JORGE R. LEIS AND MARK A. KRAMER, ODESSA - AN ORDINARY
1436 C DIFFERENTIAL EQUATION KppSolveR WITH EXPLICIT SIMULTANEOUS
1437 C SENSITIVITY ANALYSIS, SUBMITTED TO ACM TRANS. MATH. SOFTWARE.
1438 C (1985).
1440 C 3. ALAN C. HINDMARSH, LSODE AND LSODI, TWO NEW INITIAL VALUE
1441 C ORDINARY DIFFERENTIAL EQUATION KppSolveRS, ACM-SIGNUM NEWSLETTER,
1442 C VOL. 15, NO. 4 (1980), PP. 10-11.
1443 C-----------------------------------------------------------------------
1444 C THE FOLLOWING INTERNAL COMMON BLOCKS CONTAIN
1445 C (A) VARIABLES WHICH ARE LOCAL TO ANY SUBROUTINE BUT WHOSE VALUES MUST
1446 C BE PRESERVED BETWEEN CALLS TO THE ROUTINE (OWN VARIABLES), AND
1447 C (B) VARIABLES WHICH ARE COMMUNICATED BETWEEN SUBROUTINES.
1448 C THE STRUCTURE OF THE BLOCKS ARE AS FOLLOWS.. ALL REAL VARIABLES ARE
1449 C LISTED FIRST, FOLLOWED BY ALL INTEGERS. WITHIN EACH TYPE, THE
1450 C VARIABLES ARE GROUPED WITH THOSE LOCAL TO SUBROUTINE ODESSA FIRST,
1451 C THEN THOSE LOCAL TO SUBROUTINE STODE, AND FINALLY THOSE USED
1452 C FOR COMMUNICATION. THE BLOCKS ARE DECLARED IN SUBROUTINES ODESSA
1453 C INTDY, STODE, STESA, PREPJ, PREPDF, AND SOLSY. GROUPS OF VARIABLES
1454 C ARE REPLACED BY DUMMY ARRAYS IN THE COMMON DECLARATIONS IN ROUTINES
1455 C WHERE THOSE VARIABLES ARE NOT USED.
1456 C-----------------------------------------------------------------------
1457 COMMON /ODE001/ TRET, ROWNS(173),
1458 1 TESCO(3,12), CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND,
1459 2 ILLIN, INIT, LYH, LEWT, LACOR, LSAVF, LWM, LIWM,
1460 3 MXSTEP, MXHNIL, NHNIL, NTREP, NSLAST, NYH, IOWNS(4),
1461 4 IALTH, LMAX, ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH,
1462 5 MITER, MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
1463 COMMON /ODE002/ DUPS, DSMS, DDNS,
1464 1 NPAR, LDFDP, LNRS,
1465 2 ISOPT, NSV, NDFE, NSPE, IDF, IERSP, JOPT, KFLAGS
1466 PARAMETER (ZERO=0.0D0,ONE=1.0D0,TWO=2.0D0,FOUR=4.0D0)
1467 DATA MORD(1),MORD(2)/12,5/, MXSTP0/500/, MXHNL0/10/
1468 C-----------------------------------------------------------------------
1469 C BLOCK A.
1470 C THIS CODE BLOCK IS EXECUTED ON EVERY CALL.
1471 C IT TESTS ISTATE AND ITASK FOR LEGALITY AND BRANCHES APPROPIATELY.
1472 C IF ISTATE .GT. 1 BUT THE FLAG INIT SHOWS THAT INITIALIZATION HAS
1473 C NOT YET BEEN DONE, AN ERROR RETURN OCCURS.
1474 C IF ISTATE = 1 AND TOUT = T, JUMP TO BLOCK G AND RETURN IMMEDIATELY.
1475 C-----------------------------------------------------------------------
1476 IF (ISTATE .LT. 1 .OR. ISTATE .GT. 3) GO TO 601
1477 IF (ITASK .LT. 1 .OR. ITASK .GT. 5) GO TO 602
1478 IF (ISTATE .EQ. 1) GO TO 10
1479 IF (INIT .EQ. 0) GO TO 603
1480 IF (ISTATE .EQ. 2) GO TO 200
1481 GO TO 20
1482 10 INIT = 0
1483 IF (TOUT .EQ. T) GO TO 430
1484 20 NTREP = 0
1485 C-----------------------------------------------------------------------
1486 C BLOCK B.
1487 C THE NEXT CODE BLOCK IS EXECUTED FOR THE INITIAL CALL (ISTATE = 1),
1488 C OR FOR A CONTINUATION CALL WITH PARAMETER CHANGES (ISTATE = 3).
1489 C IT CONTAINS CHECKING OF ALL INPUTS AND VARIOUS INITIALIZATIONS.
1491 C FIRST CHECK LEGALITY OF THE NON-OPTIONAL INPUTS NEQ, ITOL, IOPT,
1492 C MF, ML, AND MU.
1493 C-----------------------------------------------------------------------
1494 IF (NEQ(1) .LE. 0) GO TO 604
1495 IF (ISTATE .EQ. 1) GO TO 25
1496 IF (NEQ(1) .NE. N) GO TO 605
1497 25 N = NEQ(1)
1498 IF (ITOL .LT. 1 .OR. ITOL .GT. 4) GO TO 606
1499 DO 26 I = 1,3
1500 26 IF (IOPT(I) .LT. 0 .OR. IOPT(I) .GT. 1) GO TO 607
1501 ISOPT = IOPT(2)
1502 IDF = IOPT(3)
1503 NYH = N
1504 NSV = 1
1505 METH = MF/10
1506 MITER = MF - 10*METH
1507 IF (METH .LT. 1 .OR. METH .GT. 2) GO TO 608
1508 IF (MITER .LT. 0 .OR. MITER .GT. 5) GO TO 608
1509 IF (MITER .LE. 3) GO TO 30
1510 ML = IWORK(1)
1511 MU = IWORK(2)
1512 IF (ML .LT. 0 .OR. ML .GE. N) GO TO 609
1513 IF (MU .LT. 0 .OR. MU .GE. N) GO TO 610
1514 30 IF (ISOPT .EQ. 0) GO TO 32
1515 C CHECK LEGALITY OF THE NON-OPTIONAL INPUTS ISOPT, NPAR.
1516 C COMPUTE NUMBER OF SOLUTION VECTORS AND TOTAL NUMBER OF EQUATIONS.
1517 IF (NEQ(2) .LE. 0) GO TO 628
1518 IF (ISTATE .EQ. 1) GO TO 31
1519 IF (NEQ(2) .NE. NPAR) GO TO 629
1520 31 NPAR = NEQ(2)
1521 NSV = NPAR + 1
1522 NYH = NSV * N
1523 IF (MITER .EQ. 0 .OR. MITER .EQ. 3) GO TO 630
1524 C NEXT PROCESS AND CHECK THE OPTIONAL INPUTS. --------------------------
1525 32 IF (IOPT(1) .EQ. 1) GO TO 40
1526 MAXORD = MORD(METH)
1527 MXSTEP = MXSTP0
1528 MXHNIL = MXHNL0
1529 IF (ISTATE .EQ. 1) H0 = ZERO
1530 HMXI = ZERO
1531 HMIN = ZERO
1532 GO TO 60
1533 40 MAXORD = IWORK(5)
1534 IF (MAXORD .LT. 0) GO TO 611
1535 IF (MAXORD .EQ. 0) MAXORD = 100
1536 MAXORD = MIN(MAXORD,MORD(METH))
1537 MXSTEP = IWORK(6)
1538 IF (MXSTEP .LT. 0) GO TO 612
1539 IF (MXSTEP .EQ. 0) MXSTEP = MXSTP0
1540 MXHNIL = IWORK(7)
1541 IF (MXHNIL .LT. 0) GO TO 613
1542 IF (MXHNIL .EQ. 0) MXHNIL = MXHNL0
1543 IF (ISTATE .NE. 1) GO TO 50
1544 H0 = RWORK(5)
1545 IF ((TOUT - T)*H0 .LT. ZERO) GO TO 614
1546 50 HMAX = RWORK(6)
1547 IF (HMAX .LT. ZERO) GO TO 615
1548 HMXI = ZERO
1549 IF (HMAX .GT. ZERO) HMXI = ONE/HMAX
1550 HMIN = RWORK(7)
1551 IF (HMIN .LT. ZERO) GO TO 616
1552 C-----------------------------------------------------------------------
1553 C SET WORK ARRAY POINTERS AND CHECK LENGTHS LRW AND LIW.
1554 C POINTERS TO SEGMENTS OF RWORK AND IWORK ARE NAMED BY PREFIXING L TO
1555 C THE NAME OF THE SEGMENT. E.G., THE SEGMENT YH STARTS AT RWORK(LYH).
1556 C SEGMENTS OF RWORK (IN ORDER) ARE DENOTED YH, WM, EWT, SAVF, ACOR.
1557 C WORK SPACE FOR DFDP IS CONTAINED IN ACOR.
1558 C-----------------------------------------------------------------------
1559 60 LYH = 21
1560 LWM = LYH + (MAXORD + 1)*NYH
1561 IF (MITER .EQ. 0) LENWM = 0
1562 IF (MITER .EQ. 1 .OR. MITER .EQ. 2) LENWM = N*N + 2
1563 IF (MITER .EQ. 3) LENWM = N + 2
1564 IF (MITER .GE. 4) LENWM = (2*ML + MU + 1)*N + 2
1565 LEWT = LWM + LENWM
1566 LSAVF = LEWT + NYH
1567 LACOR = LSAVF + N
1568 LDFDP = LACOR + N
1569 LENRW = LACOR + NYH - 1
1570 IWORK(17) = LENRW
1571 LIWM = 1
1572 LENIW = 20 + N
1573 IF (MITER .EQ. 0 .OR. MITER .EQ. 3) LENIW = 20
1574 LNRS = LENIW + LIWM
1575 IF (ISOPT .EQ. 1) LENIW = LNRS + NPAR
1576 IWORK(18) = LENIW
1577 IF (LENRW .GT. LRW) GO TO 617
1578 IF (LENIW .GT. LIW) GO TO 618
1579 C CHECK RTOL AND ATOL FOR LEGALITY. ------------------------------------
1580 RTOLI = RTOL(1)
1581 ATOLI = ATOL(1)
1582 DO 70 I = 1,NYH
1583 IF (ITOL .GE. 3) RTOLI = RTOL(I)
1584 IF (ITOL .EQ. 2 .OR. ITOL .EQ. 4) ATOLI = ATOL(I)
1585 IF (RTOLI .LT. ZERO) GO TO 619
1586 IF (ATOLI .LT. ZERO) GO TO 620
1587 70 CONTINUE
1588 IF (ISTATE .EQ. 1) GO TO 100
1589 C IF ISTATE = 3, SET FLAG TO SIGNAL PARAMETER CHANGES TO STODE. --------
1590 JSTART = -1
1591 IF (NQ .LE. MAXORD) GO TO 90
1592 C MAXORD WAS REDUCED BELOW NQ. COPY YH(*,MAXORD+2) INTO SAVF. ---------
1593 DO 80 I = 1,N
1594 80 RWORK(I+LSAVF-1) = RWORK(I+LWM-1)
1595 C RELOAD WM(1) = RWORK(LWM), SINCE LWM MAY HAVE CHANGED. ---------------
1596 90 IF (MITER .GT. 0) RWORK(LWM) = SQRT(UROUND)
1597 GO TO 200
1598 C-----------------------------------------------------------------------
1599 C BLOCK C.
1600 C THE NEXT BLOCK IS FOR THE INITIAL CALL ONLY (ISTATE = 1).
1601 C IT CONTAINS ALL REMAINING INITIALIZATIONS, THE INITIAL CALL TO F,
1602 C THE INITIAL CALL TO SPRIME IF ISOPT = 1,
1603 C AND THE CALCULATION OF THE INITIAL STEP SIZE.
1604 C THE ERROR WEIGHTS IN EWT ARE INVERTED AFTER BEING LOADED.
1605 C-----------------------------------------------------------------------
1606 100 UROUND = D1MACH(4)
1607 TN = T
1608 IF (ITASK .NE. 4 .AND. ITASK .NE. 5) GO TO 105
1609 TCRIT = RWORK(1)
1610 IF ((TCRIT - TOUT)*(TOUT - T) .LT. ZERO) GO TO 625
1611 IF (H0 .NE. ZERO .AND. (T + H0 - TCRIT)*H0 .GT. ZERO)
1612 1 H0 = TCRIT - T
1613 105 JSTART = 0
1614 IF (MITER .GT. 0) RWORK(LWM) = SQRT(UROUND)
1615 NHNIL = 0
1616 NST = 0
1617 NJE = 0
1618 NSLAST = 0
1619 HU = ZERO
1620 NQU = 0
1621 CCMAX = 0.3D0
1622 MAXCOR = 3
1623 IF (ISOPT .EQ. 1) MAXCOR = 4
1624 MSBP = 20
1625 MXNCF = 10
1626 C INITIAL CALL TO F. (LF0 POINTS TO YH(1,2) AND LOADS IN VALUES).
1627 LF0 = LYH + NYH
1628 CALL F (NEQ, T, Y, PAR, RWORK(LF0))
1629 NFE = 1
1630 DUPS = ZERO
1631 DSMS = ZERO
1632 DDNS = ZERO
1633 NDFE = 0
1634 NSPE = 0
1635 IF (ISOPT .EQ. 0) GO TO 114
1636 C INITIALIZE COUNTS FOR REPEATED STEPS DUE TO SENSITIVITY ANALYSIS.
1637 DO 110 J = 1,NSV
1638 110 IWORK(J + LNRS - 1) = 0
1639 C LOAD THE INITIAL VALUE VECTOR IN YH. ---------------------------------
1640 114 DO 115 I = 1,NYH
1641 115 RWORK(I+LYH-1) = Y(I)
1642 C LOAD AND INVERT THE EWT ARRAY. (H IS TEMPORARILY SET TO ONE.) -------
1643 NQ = 1
1644 H = ONE
1645 CALL EWSET (NYH, ITOL, RTOL, ATOL, RWORK(LYH), RWORK(LEWT))
1646 DO 120 I = 1,NYH
1647 IF (RWORK(I+LEWT-1) .LE. ZERO) GO TO 621
1648 120 RWORK(I+LEWT-1) = ONE/RWORK(I+LEWT-1)
1649 IF (ISOPT .EQ. 0) GO TO 125
1650 C CALL SPRIME TO LOAD FIRST-ORDER SENSITIVITY DERIVATIVES INTO
1651 C REMAINING YH(*,2) POSITIONS.
1652 CALL SPRIME (NEQ, Y, RWORK(LYH), NYH, N, NSV, RWORK(LWM),
1653 1 IWORK(LIWM), RWORK(LEWT), RWORK(LF0), RWORK(LACOR),
1654 2 RWORK(LDFDP), PAR, F, JAC, DF, PREPJ, PREPDF)
1655 IF (IERSP .EQ. -1) GO TO 631
1656 IF (IERSP .EQ. -2) GO TO 632
1657 C-----------------------------------------------------------------------
1658 C THE CODING BELOW COMPUTES THE STEP SIZE, H0, TO BE ATTEMPTED ON THE
1659 C FIRST STEP, UNLESS THE USER HAS SUPPLIED A VALUE FOR THIS.
1660 C FIRST CHECK THAT TOUT - T DIFFERS SIGNIFICANTLY FROM ZERO.
1661 C A SCALAR TOLERANCE QUANTITY TOL IS COMPUTED, AS MAX(RTOL(I))
1662 C IF THIS IS POSITIVE, OR MAX(ATOL(I)/ABS(Y(I))) OTHERWISE, ADJUSTED
1663 C SO AS TO BE BETWEEN 100*UROUND AND 1.0E-3. ONLY THE ORIGINAL
1664 C SOLUTION VECTOR IS CONSIDERED IN THIS CALCULATION (ISOPT = 0 OR 1).
1665 C THEN THE COMPUTED VALUE H0 IS GIVEN BY..
1666 C NEQ
1667 C H0**2 = TOL / ( W0**-2 + (1/NEQ) * SUM ( F(I)/YWT(I) )**2 )
1669 C WHERE W0 = MAX ( ABS(T), ABS(TOUT) ),
1670 C F(I) = I-TH COMPONENT OF INITIAL VALUE OF F,
1671 C YWT(I) = EWT(I)/TOL (A WEIGHT FOR Y(I)).
1672 C THE SIGN OF H0 IS INFERRED FROM THE INITIAL VALUES OF TOUT AND T.
1673 C-----------------------------------------------------------------------
1674 125 IF (H0 .NE. ZERO) GO TO 180
1675 TDIST = ABS(TOUT - T)
1676 W0 = MAX(ABS(T),ABS(TOUT))
1677 IF (TDIST .LT. TWO*UROUND*W0) GO TO 622
1678 TOL = RTOL(1)
1679 IF (ITOL .LE. 2) GO TO 140
1680 DO 130 I = 1,N
1681 130 TOL = MAX(TOL,RTOL(I))
1682 140 IF (TOL .GT. ZERO) GO TO 160
1683 ATOLI = ATOL(1)
1684 DO 150 I = 1,N
1685 IF (ITOL .EQ. 2 .OR. ITOL .EQ. 4) ATOLI = ATOL(I)
1686 AYI = ABS(Y(I))
1687 IF (AYI .NE. ZERO) TOL = MAX(TOL,ATOLI/AYI)
1688 150 CONTINUE
1689 160 TOL = MAX(TOL,100.0D0*UROUND)
1690 TOL = MIN(TOL,0.001D0)
1691 SUM = VNORM (N, RWORK(LF0), RWORK(LEWT))
1692 SUM = ONE/(TOL*W0*W0) + TOL*SUM**2
1693 H0 = ONE/SQRT(SUM)
1694 H0 = MIN(H0,TDIST)
1695 H0 = SIGN(H0,TOUT-T)
1696 C ADJUST H0 IF NECESSARY TO MEET HMAX BOUND. ---------------------------
1697 180 RH = ABS(H0)*HMXI
1698 IF (RH .GT. ONE) H0 = H0/RH
1699 C LOAD H WITH H0 AND SCALE YH(*,2) BY H0. ------------------------------
1700 H = H0
1701 DO 190 I = 1,NYH
1702 190 RWORK(I+LF0-1) = H0*RWORK(I+LF0-1)
1703 GO TO 270
1704 C-----------------------------------------------------------------------
1705 C BLOCK D.
1706 C THE NEXT CODE BLOCK IS FOR CONTINUATION CALLS ONLY (ISTATE = 2 OR 3)
1707 C AND IS TO CHECK STOP CONDITIONS BEFORE TAKING A STEP.
1708 C-----------------------------------------------------------------------
1709 200 NSLAST = NST
1710 GO TO (210, 250, 220, 230, 240), ITASK
1711 210 IF ((TN - TOUT)*H .LT. ZERO) GO TO 250
1712 CALL INTDY (TOUT, 0, RWORK(LYH), NYH, Y, IFLAG)
1713 IF (IFLAG .NE. 0) GO TO 627
1714 T = TOUT
1715 GO TO 420
1716 220 TP = TN - HU*(ONE + 100.0D0*UROUND)
1717 IF ((TP - TOUT)*H .GT. ZERO) GO TO 623
1718 IF ((TN - TOUT)*H .LT. ZERO) GO TO 250
1719 GO TO 400
1720 230 TCRIT = RWORK(1)
1721 IF ((TN - TCRIT)*H .GT. ZERO) GO TO 624
1722 IF ((TCRIT - TOUT)*H .LT. ZERO) GO TO 625
1723 IF ((TN - TOUT)*H .LT. ZERO) GO TO 245
1724 CALL INTDY (TOUT, 0, RWORK(LYH), NYH, Y, IFLAG)
1725 IF (IFLAG .NE. 0) GO TO 627
1726 T = TOUT
1727 GO TO 420
1728 240 TCRIT = RWORK(1)
1729 IF ((TN - TCRIT)*H .GT. ZERO) GO TO 624
1730 245 HMX = ABS(TN) + ABS(H)
1731 IHIT = ABS(TN - TCRIT) .LE. 100.0D0*UROUND*HMX
1732 IF (IHIT) GO TO 400
1733 TNEXT = TN + H*(ONE + FOUR*UROUND)
1734 IF ((TNEXT - TCRIT)*H .LE. ZERO) GO TO 250
1735 H = (TCRIT - TN)*(ONE - FOUR*UROUND)
1736 IF (ISTATE .EQ. 2) JSTART = -2
1737 C-----------------------------------------------------------------------
1738 C BLOCK E.
1739 C THE NEXT BLOCK IS NORMALLY EXECUTED FOR ALL CALLS AND CONTAINS
1740 C THE CALL TO THE ONE-STEP CORE INTEGRATOR STODE.
1742 C THIS IS A LOOPING POINT FOR THE INTEGRATION STEPS.
1744 C FIRST CHECK FOR TOO MANY STEPS BEING TAKEN, UPDATE EWT (IF NOT AT
1745 C START OF PROBLEM), CHECK FOR TOO MUCH ACCURACY BEING REQUESTED, AND
1746 C CHECK FOR H BELOW THE ROUNDOFF LEVEL IN T.
1747 C TOLSF IS CALCULATED CONSIDERING ALL SOLUTION VECTORS.
1748 C-----------------------------------------------------------------------
1749 250 CONTINUE
1750 IF ((NST-NSLAST) .GE. MXSTEP) GO TO 500
1751 CALL EWSET (NYH, ITOL, RTOL, ATOL, RWORK(LYH), RWORK(LEWT))
1752 DO 260 I = 1,NYH
1753 IF (RWORK(I+LEWT-1) .LE. ZERO) GO TO 510
1754 260 RWORK(I+LEWT-1) = ONE/RWORK(I+LEWT-1)
1755 270 TOLSF = UROUND*VNORM (NYH, RWORK(LYH), RWORK(LEWT))
1756 IF (TOLSF .LE. ONE) GO TO 280
1757 TOLSF = TOLSF*2.0D0
1758 IF (NST .EQ. 0) GO TO 626
1759 GO TO 520
1760 280 IF (ADDX(TN,H) .NE. TN) GO TO 290
1761 NHNIL = NHNIL + 1
1762 IF (NHNIL .GT. MXHNIL) GO TO 290
1763 CALL XERR ('ODESSA - WARNING..INTERNAL T (=R1) AND H (=R2) ARE',
1764 1 101, 1, 0, 0, 0, 0, ZERO, ZERO)
1765 CALL XERR ('SUCH THAT IN THE MACHINE, T + H = T ON THE NEXT STEP',
1766 1 101, 1, 0, 0, 0, 0, ZERO, ZERO)
1767 CALL XERR ('(H = STEP SIZE). KppSolveR WILL CONTINUE ANYWAY',
1768 1 101, 1, 0, 0, 0, 2, TN, H)
1769 IF (NHNIL .LT. MXHNIL) GO TO 290
1770 CALL XERR ('ODESSA - ABOVE WARNING HAS BEEN ISSUED I1 TIMES.',
1771 1 102, 1, 0, 0, 0, 0, ZERO, ZERO)
1772 CALL XERR ('IT WILL NOT BE ISSUED AGAIN FOR THIS PROBLEM',
1773 1 102, 1, 1, MXHNIL, 0, 0, ZERO,ZERO)
1774 290 CONTINUE
1775 C-----------------------------------------------------------------------
1776 C CALL STODE(NEQ,Y,YH,NYH,YH,WM,IWM,EWT,SAVF,ACOR,PAR,NRS,
1777 C 1 F,JAC,DF,PREPJ,PREPDF,SOLSY)
1778 C-----------------------------------------------------------------------
1779 CALL STODE (NEQ, Y, RWORK(LYH), NYH, RWORK(LYH), RWORK(LWM),
1780 1 IWORK(LIWM), RWORK(LEWT), RWORK(LSAVF), RWORK(LACOR),
1781 2 PAR, IWORK(LNRS), F, JAC, DF, PREPJ, PREPDF, SOLSY)
1782 KGO = 1 - KFLAG
1783 GO TO (300, 530, 540, 633), KGO
1784 C-----------------------------------------------------------------------
1785 C BLOCK F.
1786 C THE FOLLOWING BLOCK HANDLES THE CASE OF A SUCCESSFUL RETURN FROM THE
1787 C CORE INTEGRATOR (KFLAG = 0). TEST FOR STOP CONDITIONS.
1788 C-----------------------------------------------------------------------
1789 300 INIT = 1
1790 GO TO (310, 400, 330, 340, 350), ITASK
1791 C ITASK = 1. IF TOUT HAS BEEN REACHED, INTERPOLATE. -------------------
1792 310 IF ((TN - TOUT)*H .LT. ZERO) GO TO 250
1793 CALL INTDY (TOUT, 0, RWORK(LYH), NYH, Y, IFLAG)
1794 T = TOUT
1795 GO TO 420
1796 C ITASK = 3. JUMP TO EXIT IF TOUT WAS REACHED. ------------------------
1797 330 IF ((TN - TOUT)*H .GE. ZERO) GO TO 400
1798 GO TO 250
1799 C ITASK = 4. SEE IF TOUT OR TCRIT WAS REACHED. ADJUST H IF NECESSARY.
1800 340 IF ((TN - TOUT)*H .LT. ZERO) GO TO 345
1801 CALL INTDY (TOUT, 0, RWORK(LYH), NYH, Y, IFLAG)
1802 T = TOUT
1803 GO TO 420
1804 345 HMX = ABS(TN) + ABS(H)
1805 IHIT = ABS(TN - TCRIT) .LE. 100.0D0*UROUND*HMX
1806 IF (IHIT) GO TO 400
1807 TNEXT = TN + H*(ONE + FOUR*UROUND)
1808 IF ((TNEXT - TCRIT)*H .LE. ZERO) GO TO 250
1809 H = (TCRIT - TN)*(ONE - FOUR*UROUND)
1810 JSTART = -2
1811 GO TO 250
1812 C ITASK = 5. SEE IF TCRIT WAS REACHED AND JUMP TO EXIT. ---------------
1813 350 HMX = ABS(TN) + ABS(H)
1814 IHIT = ABS(TN - TCRIT) .LE. 100.0D0*UROUND*HMX
1815 C-----------------------------------------------------------------------
1816 C BLOCK G.
1817 C THE FOLLOWING BLOCK HANDLES ALL SUCCESSFUL RETURNS FROM ODESSA.
1818 C IF ITASK .NE. 1, Y IS LOADED FROM YH AND T IS SET ACCORDINGLY.
1819 C ISTATE IS SET TO 2, THE ILLEGAL INPUT COUNTER IS ZEROED, AND THE
1820 C OPTIONAL OUTPUTS ARE LOADED INTO THE WORK ARRAYS BEFORE RETURNING.
1821 C IF ISTATE = 1 AND TOUT = T, THERE IS A RETURN WITH NO ACTION TAKEN,
1822 C EXCEPT THAT IF THIS HAS HAPPENED REPEATEDLY, THE RUN IS TERMINATED.
1823 C-----------------------------------------------------------------------
1824 400 DO 410 I = 1,NYH
1825 410 Y(I) = RWORK(I+LYH-1)
1826 T = TN
1827 IF (ITASK .NE. 4 .AND. ITASK .NE. 5) GO TO 420
1828 IF (IHIT) T = TCRIT
1829 420 ISTATE = 2
1830 ILLIN = 0
1831 RWORK(11) = HU
1832 RWORK(12) = H
1833 RWORK(13) = TN
1834 IWORK(11) = NST
1835 IWORK(12) = NFE
1836 IWORK(13) = NJE
1837 IWORK(14) = NQU
1838 IWORK(15) = NQ
1839 IF (ISOPT .EQ. 0) RETURN
1840 IWORK(19) = NDFE
1841 IWORK(20) = NSPE
1842 RETURN
1843 430 NTREP = NTREP + 1
1844 IF (NTREP .LT. 5) RETURN
1845 CALL XERR ('ODESSA -- REPEATED CALLS WITH ISTATE = 1 AND
1846 1TOUT = T (=R1)', 301, 1, 0, 0, 0, 1, T, ZERO)
1847 GO TO 800
1848 C-----------------------------------------------------------------------
1849 C BLOCK H.
1850 C THE FOLLOWING BLOCK HANDLES ALL UNSUCCESSFUL RETURNS OTHER THAN
1851 C THOSE FOR ILLEGAL INPUT. FIRST THE ERROR MESSAGE ROUTINE IS CALLED.
1852 C IF THERE WAS AN ERROR TEST OR CONVERGENCE TEST FAILURE, IMXER IS SET.
1853 C THEN Y IS LOADED FROM YH, T IS SET TO TN, AND THE ILLEGAL INPUT
1854 C COUNTER ILLIN IS SET TO 0. THE OPTIONAL OUTPUTS ARE LOADED INTO
1855 C THE WORK ARRAYS BEFORE RETURNING.
1856 C-----------------------------------------------------------------------
1857 C THE MAXIMUM NUMBER OF STEPS WAS TAKEN BEFORE REACHING TOUT. ----------
1858 500 CALL XERR ('ODESSA - AT CURRENT T (=R1), MXSTEP (=I1) STEPS',
1859 1 201, 1, 0, 0, 0, 0, ZERO,ZERO)
1860 CALL XERR ('TAKEN ON THIS CALL BEFORE REACHING TOUT',
1861 1 201, 1, 1, MXSTEP, 0, 1, TN, ZERO)
1862 ISTATE = -1
1863 GO TO 580
1864 C EWT(I) .LE. 0.0 FOR SOME I (NOT AT START OF PROBLEM). ----------------
1865 510 EWTI = RWORK(LEWT+I-1)
1866 CALL XERR ('ODESSA - AT T (=R1), EWT(I1) HAS BECOME R2 .LE. 0.',
1867 1 202, 1, 1, I, 0, 2, TN, EWTI)
1868 ISTATE = -6
1869 GO TO 580
1870 C TOO MUCH ACCURACY REQUESTED FOR MACHINE PRECISION. -------------------
1871 520 CALL XERR ('ODESSA - AT T (=R1), TOO MUCH ACCURACY REQUESTED',
1872 1 203, 1, 0, 0, 0, 0, ZERO,ZERO)
1873 CALL XERR ('FOR PRECISION OF MACHINE.. SEE TOLSF (=R2)',
1874 1 203, 1, 0, 0, 0, 2, TN, TOLSF)
1875 RWORK(14) = TOLSF
1876 ISTATE = -2
1877 GO TO 580
1878 C KFLAG = -1. ERROR TEST FAILED REPEATEDLY OR WITH ABS(H) = HMIN. -----
1879 530 CALL XERR ('ODESSA - AT T(=R1) AND STEP SIZE H(=R2), THE ERROR',
1880 1 204, 1, 0, 0, 0, 0, ZERO, ZERO)
1881 CALL XERR ('TEST FAILED REPEATEDLY OR WITH ABS(H) = HMIN',
1882 1 204, 1, 0, 0, 0, 2, TN, H)
1883 ISTATE = -4
1884 GO TO 560
1885 C KFLAG = -2. CONVERGENCE FAILED REPEATEDLY OR WITH ABS(H) = HMIN. ----
1886 540 CALL XERR ('ODESSA - AT T (=R1) AND STEP SIZE H (=R2), THE',
1887 1 205, 1, 0, 0, 0, 0, ZERO,ZERO)
1888 CALL XERR ('CORRECTOR CONVERGENCE FAILED REPEATEDLY',
1889 1 205, 1, 0, 0, 0, 0, ZERO,ZERO)
1890 CALL XERR ('OR WITH ABS(H) = HMIN',
1891 1 205, 1, 0, 0, 0, 2, TN, H)
1892 ISTATE = -5
1893 C COMPUTE IMXER IF RELEVANT. -------------------------------------------
1894 560 BIG = ZERO
1895 IMXER = 1
1896 DO 570 I = 1,NYH
1897 SIZE = ABS(RWORK(I+LACOR-1)*RWORK(I+LEWT-1))
1898 IF (BIG .GE. SIZE) GO TO 570
1899 BIG = SIZE
1900 IMXER = I
1901 570 CONTINUE
1902 IWORK(16) = IMXER
1903 C SET Y VECTOR, T, ILLIN, AND OPTIONAL OUTPUTS. ------------------------
1904 580 DO 590 I = 1,NYH
1905 590 Y(I) = RWORK(I+LYH-1)
1906 T = TN
1907 ILLIN = 0
1908 RWORK(11) = HU
1909 RWORK(12) = H
1910 RWORK(13) = TN
1911 IWORK(11) = NST
1912 IWORK(12) = NFE
1913 IWORK(13) = NJE
1914 IWORK(14) = NQU
1915 IWORK(15) = NQ
1916 IF (ISOPT .EQ. 0) RETURN
1917 IWORK(19) = NDFE
1918 IWORK(20) = NSPE
1919 RETURN
1920 C-----------------------------------------------------------------------
1921 C BLOCK I.
1922 C THE FOLLOWING BLOCK HANDLES ALL ERROR RETURNS DUE TO ILLEGAL INPUT
1923 C (ISTATE = -3), AS DETECTED BEFORE CALLING THE CORE INTEGRATOR.
1924 C FIRST THE ERROR MESSAGE ROUTINE IS CALLED. THEN IF THERE HAVE BEEN
1925 C 5 CONSECUTIVE SUCH RETURNS JUST BEFORE THIS CALL TO THE KppSolveR,
1926 C THE RUN IS HALTED.
1927 C-----------------------------------------------------------------------
1928 601 CALL XERR ('ODESSA - ISTATE (=I1) ILLEGAL',
1929 1 1, 1, 1, ISTATE, 0, 0, ZERO,ZERO)
1930 GO TO 700
1931 602 CALL XERR ('ODESSA - ITASK (=I1) ILLEGAL',
1932 1 2, 1, 1, ITASK, 0, 0, ZERO,ZERO)
1933 GO TO 700
1934 603 CALL XERR ('ODESSA - ISTATE .GT. 1 BUT ODESSA NOT INITIALIZED',
1935 1 3, 1, 0, 0, 0, 0, ZERO,ZERO)
1936 GO TO 700
1937 604 CALL XERR ('ODESSA - NEQ (=I1) .LT. 1',
1938 1 4, 1, 1, NEQ(1), 0, 0, ZERO,ZERO)
1939 GO TO 700
1940 605 CALL XERR ('ODESSA - ISTATE = 3 AND NEQ CHANGED. (I1 TO I2)',
1941 1 5, 1, 2, N, NEQ(1), 0, ZERO,ZERO)
1942 GO TO 700
1943 606 CALL XERR ('ODESSA - ITOL (=I1) ILLEGAL',
1944 1 6, 1, 1, ITOL, 0, 0, ZERO,ZERO)
1945 GO TO 700
1946 607 CALL XERR ('ODESSA - IOPT (=I1) ILLEGAL',
1947 1 7, 1, 1, IOPT, 0, 0, ZERO,ZERO)
1948 GO TO 700
1949 608 CALL XERR('ODESSA - MF (=I1) ILLEGAL',
1950 1 8, 1, 1, MF, 0, 0, ZERO,ZERO)
1951 GO TO 700
1952 609 CALL XERR('ODESSA - ML (=I1) ILLEGAL.. .LT.0 OR .GE.NEQ (=I2)',
1953 1 9, 1, 2, ML, NEQ(1), 0, ZERO,ZERO)
1954 GO TO 700
1955 610 CALL XERR('ODESSA - MU (=I1) ILLEGAL.. .LT.0 OR .GE.NEQ (=I2)',
1956 1 10, 1, 2, MU, NEQ(1), 0, ZERO,ZERO)
1957 GO TO 700
1958 611 CALL XERR('ODESSA - MAXORD (=I1) .LT. 0',
1959 1 11, 1, 1, MAXORD, 0, 0, ZERO,ZERO)
1960 GO TO 700
1961 612 CALL XERR('ODESSA - MXSTEP (=I1) .LT. 0',
1962 1 12, 1, 1, MXSTEP, 0, 0, ZERO,ZERO)
1963 GO TO 700
1964 613 CALL XERR('ODESSA - MXHNIL (=I1) .LT. 0',
1965 1 13, 1, 1, MXHNIL, 0, 0, ZERO,ZERO)
1966 GO TO 700
1967 614 CALL XERR('ODESSA - TOUT (=R1) BEHIND T (=R2)',
1968 1 14, 1, 0, 0, 0, 2, TOUT, T)
1969 CALL XERR('INTEGRATION DIRECTION IS GIVEN BY H0 (=R1)',
1970 1 14, 1, 0, 0, 0, 1, H0, ZERO)
1971 GO TO 700
1972 615 CALL XERR('ODESSA - HMAX (=R1) .LT. 0.0',
1973 1 15, 1, 0, 0, 0, 1, HMAX, ZERO)
1974 GO TO 700
1975 616 CALL XERR('ODESSA - HMIN (=R1) .LT. 0.0',
1976 1 16, 1, 0, 0, 0, 1, HMIN, ZERO)
1977 GO TO 700
1978 617 CALL XERR('ODESSA - RWORK LENGTH NEEDED, LENRW (=I1), EXCEEDS
1979 1 LRW (=I2)', 17, 1, 2, LENRW, LRW, 0, ZERO,ZERO)
1980 GO TO 700
1981 618 CALL XERR('ODESSA - IWORK LENGTH NEEDED, LENIW (=I1), EXCEEDS
1982 1 LIW (=I2)', 18, 1, 2, LENIW, LIW, 0, ZERO,ZERO)
1983 GO TO 700
1984 619 CALL XERR('ODESSA - RTOL(I1) IS R1 .LT. 0.0',
1985 1 19, 1, 1, I, 0, 1, RTOLI, ZREO)
1986 GO TO 700
1987 620 CALL XERR('ODESSA - ATOL(I1) IS R1 .LT. 0.0',
1988 1 20, 1, 1, I, 0, 1, ATOLI, ZERO)
1989 GO TO 700
1991 621 EWTI = RWORK(LEWT+I-1)
1992 CALL XERR('ODESSA - EWT(I1) IS R1 .LE. 0.0',
1993 1 21, 1, 1, I, 0, 1, EWTI, ZERO)
1994 GO TO 700
1995 622 CALL XERR('ODESSA - TOUT (=R1) TOO CLOSE TO T(=R2) TO START
1996 1 INTEGRATION', 22, 1, 0, 0, 0, 2, TOUT, T)
1997 GO TO 700
1998 623 CALL XERR('ODESSA - ITASK = I1 AND TOUT (=R1) BEHIND TCUR - HU
1999 1 (= R2)', 23, 1, 1, ITASK, 0, 2, TOUT, TP)
2000 GO TO 700
2001 624 CALL XERR('ODESSA - ITASK = 4 OR 5 AND TCRIT (=R1) BEHIND TCUR
2002 1 (=R2)', 24, 1, 0, 0, 0, 2, TCRIT, TN)
2003 GO TO 700
2004 625 CALL XERR('ODESSA - ITASK = 4 OR 5 AND TCRIT (=R1) BEHIND TOUT
2005 1 (=R2)', 25, 1, 0, 0, 0, 2, TCRIT, TOUT)
2006 GO TO 700
2007 626 CALL XERR('ODESSA - AT START OF PROBLEM, TOO MUCH ACCURACY',
2008 1 26, 1, 0, 0, 0, 0, ZERO,ZERO)
2009 CALL XERR('REQUESTED FOR PRECISION OF MACHINE. SEE TOLSF (=R1)',
2010 1 26, 1, 0, 0, 0, 1, TOLSF, ZERO)
2011 RWORK(14) = TOLSF
2012 GO TO 700
2013 627 CALL XERR('ODESSA - TROUBLE FROM INTDY. ITASK = I1, TOUT = R1',
2014 1 27, 1, 1, ITASK, 0, 1, TOUT, ZERO)
2015 GO TO 700
2016 C ERROR STATEMENTS ASSOCIATED WITH SENSITIVITY ANALYSIS.
2017 628 CALL XERR('ODESSA - NPAR (=I1) .LT. 1',
2018 1 28, 1, 1, NPAR, 0, 0, ZERO,ZERO)
2019 GO TO 700
2020 629 CALL XERR('ODESSA - ISTATE = 3 AND NPAR CHANGED (I1 TO I2)',
2021 1 29, 1, 2, NP, NPAR, 0, ZERO,ZERO)
2022 GO TO 700
2023 630 CALL XERR('ODESSA - MITER (=I1) ILLEGAL',
2024 1 30, 1, 1, MITER, 0, 0, ZERO,ZERO)
2025 GO TO 700
2026 631 CALL XERR('ODESSA - TROUBLE IN SPRIME (IERPJ)',
2027 1 31, 1, 0, 0, 0, 0, ZERO,ZERO)
2028 GO TO 700
2029 632 CALL XERR('ODESSA - TROUBLE IN SPRIME (MITER)',
2030 1 32, 1, 0, 0, 0, 0, ZERO,ZERO)
2031 GO TO 700
2032 633 CALL XERR('ODESSA - FATAL ERROR IN STODE (KFLAG = -3)',
2033 1 33, 2, 0, 0, 0, 0, ZERO,ZERO)
2034 GO TO 801
2036 700 IF (ILLIN .EQ. 5) GO TO 710
2037 ILLIN = ILLIN + 1
2038 ISTATE = -3
2039 RETURN
2040 710 CALL XERR('ODESSA - REPEATED OCCURRENCES OF ILLEGAL INPUT',
2041 1 302, 1, 0, 0, 0, 0, ZERO,ZERO)
2043 800 CALL XERR('ODESSA - RUN ABORTED.. APPARENT INFINITE LOOP',
2044 1 303, 2, 0, 0, 0, 0, ZERO,ZERO)
2045 RETURN
2046 801 CALL XERR('ODESSA - RUN ABORTED',
2047 1 304, 2, 0, 0, 0, 0, ZERO,ZERO)
2048 RETURN
2049 C-------------------- END OF SUBROUTINE ODESSA -------------------------
2051 DOUBLE PRECISION FUNCTION ADDX(A,B)
2052 DOUBLE PRECISION A,B
2054 C THIS FUNCTION IS NECESSARY TO FORCE OPTIMIZING COMPILERS TO
2055 C EXECUTE AND STORE A SUM, FOR SUCCESSFUL EXECUTION OF THE
2056 C TEST A + B = B.
2058 ADDX = A + B
2059 RETURN
2060 C-------------------- END OF FUNCTION SUM ------------------------------
2062 SUBROUTINE SPRIME (NEQ, Y, YH, NYH, NROW, NCOL, WM, IWM,
2063 1 EWT, SAVF, FTEM, DFDP, PAR, F, JAC, DF, PJAC, PDF)
2064 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2065 DIMENSION NEQ(*), Y(*), YH(NROW,NCOL,*), WM(*), IWM(*),
2066 1 EWT(*), SAVF(*), FTEM(*), DFDP(NROW,*), PAR(*)
2067 EXTERNAL F, JAC, DF, PJAC, PDF
2068 PARAMETER (ONE=1.0D0,ZERO=0.0D0)
2069 COMMON /ODE001/ ROWND, ROWNS(173),
2070 1 RDUM1(37),EL0, H, RDUM2(6),
2071 2 IOWND1(14), IOWNS(4),
2072 3 IDUM1(3), IERPJ, IDUM2(6),
2073 4 MITER, IDUM3(4), N, IDUM4(5)
2074 COMMON /ODE002/ RDUM3(3),
2075 1 IOWND2(3), IDUM5, NSV, IDUM6, NSPE, IDUM7, IERSP, JOPT, IDUM8
2076 C-----------------------------------------------------------------------
2077 C SPRIME IS CALLED BY ODESSA TO INITIALIZE THE YH ARRAY. IT IS ALSO
2078 C CALLED BY STODE TO REEVALUATE FIRST ORDER DERIVATIVES WHEN KFLAG
2079 C .LE. -3. SPRIME COMPUTES THE FIRST DERIVATIVES OF THE SENSITIVITY
2080 C COEFFICIENTS WITH RESPECT TO THE INDEPENDENT VARIABLE T...
2082 C SPRIME = D(DY/DP)/DT = JAC*DY/DP + DF/DP
2083 C WHERE JAC = JACOBIAN MATRIX
2084 C DY/DP = SENSITIVITY MATRIX
2085 C DF/DP = INHOMOGENEITY MATRIX
2086 C THIS ROUTINE USES THE COMMON VARIABLES EL0, H, IERPJ, MITER, N,
2087 C NSV, NSPE, IERSP, JOPT
2088 C-----------------------------------------------------------------------
2089 C CALL PREPJ WITH JOPT = 1.
2090 C IF MITER = 2 OR 5, EL0 IS TEMPORARILY SET TO -1.0 AND H IS
2091 C TEMPORARILY SET TO 1.0D0.
2092 C-----------------------------------------------------------------------
2093 NSPE = NSPE + 1
2094 JOPT = 1
2095 IF (MITER .EQ. 1 .OR. MITER .EQ. 4) GO TO 10
2096 HTEMP = H
2097 ETEMP = EL0
2098 H = ONE
2099 EL0 = -ONE
2100 10 CALL PJAC (NEQ, Y, YH, NYH, WM, IWM, EWT, SAVF, FTEM,
2101 1 PAR, F, JAC, JOPT)
2102 IF (IERPJ .NE. 0) GO TO 300
2103 JOPT = 0
2104 IF (MITER .EQ. 1 .OR. MITER .EQ. 4) GO TO 20
2105 H = HTEMP
2106 EL0 = ETEMP
2107 C-----------------------------------------------------------------------
2108 C CALL PREPDF AND LOAD DFDP(*,JPAR).
2109 C-----------------------------------------------------------------------
2110 20 DO 30 J = 2,NSV
2111 JPAR = J - 1
2112 CALL PDF (NEQ, Y, WM, SAVF, FTEM, DFDP(1,JPAR), PAR,
2113 1 F, DF, JPAR)
2114 30 CONTINUE
2115 C-----------------------------------------------------------------------
2116 C COMPUTE JAC*DY/DP AND STORE RESULTS IN YH(*,*,2).
2117 C-----------------------------------------------------------------------
2118 GO TO (40,40,310,100,100) MITER
2119 C THE JACOBIAN IS FULL.------------------------------------------------
2120 C FOR EACH ROW OF THE JACOBIAN..
2121 40 DO 70 IROW = 1,N
2122 C AND EACH COLUMN OF THE SENSITIVITY MATRIX..
2123 DO 60 J = 2,NSV
2124 SUM = ZERO
2125 C TAKE THE VECTOR DOT PRODUCT..
2126 DO 50 I = 1,N
2127 IPD = IROW + N*(I-1) + 2
2128 SUM = SUM + WM(IPD)*YH(I,J,1)
2129 50 CONTINUE
2130 YH(IROW,J,2) = SUM
2131 60 CONTINUE
2132 70 CONTINUE
2133 GO TO 200
2134 C THE JACOBIAN IS BANDED.-----------------------------------------------
2135 100 ML = IWM(1)
2136 MU = IWM(2)
2137 ICOUNT = 1
2138 MBAND = ML + MU + 1
2139 MEBAND = MBAND + ML
2140 NMU = N - MU
2141 ML1 = ML + 1
2142 C FOR EACH ROW OF THE JACOBIAN..
2143 DO 160 IROW = 1,N
2144 IF (IROW .GT. ML1) GO TO 110
2145 IPD = MBAND + IROW + 1
2146 IYH = 1
2147 LBAND = MU + IROW
2148 GO TO 120
2149 110 ICOUNT = ICOUNT + 1
2150 IPD = ICOUNT*MEBAND + 2
2151 IYH = IYH + 1
2152 LBAND = LBAND - 1
2153 IF (IROW .LE. NMU) LBAND = MBAND
2154 C AND EACH COLUMN OF THE SENSITIVITY MATRIX..
2155 120 DO 150 J = 2,NSV
2156 SUM = ZERO
2157 I1 = IPD
2158 I2 = IYH
2159 C TAKE THE VECTOR DOT PRODUCT.
2160 DO 140 I = 1,LBAND
2161 SUM = SUM + WM(I1)*YH(I2,J,1)
2162 I1 = I1 + MEBAND - 1
2163 I2 = I2 + 1
2164 140 CONTINUE
2165 YH(IROW,J,2) = SUM
2166 150 CONTINUE
2167 160 CONTINUE
2168 C-----------------------------------------------------------------------
2169 C ADD THE INHOMOGENEITY TERM, I.E., ADD DFDP(*,JPAR) TO YH(*,JPAR+1,2).
2170 C-----------------------------------------------------------------------
2171 200 DO 220 J = 2,NSV
2172 JPAR = J - 1
2173 DO 210 I = 1,N
2174 YH(I,J,2) = YH(I,J,2) + DFDP(I,JPAR)
2175 210 CONTINUE
2176 220 CONTINUE
2177 RETURN
2178 C-----------------------------------------------------------------------
2179 C ERROR RETURNS.
2180 C-----------------------------------------------------------------------
2181 300 IERSP = -1
2182 RETURN
2183 310 IERSP = -2
2184 RETURN
2185 C------------------------END OF SUBROUTINE SPRIME-----------------------
2187 SUBROUTINE PREPJ (NEQ, Y, YH, NYH, WM, IWM, EWT, SAVF, FTEM,
2188 1 PAR, F, JAC, JOPT)
2189 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2190 DIMENSION NEQ(*), Y(*), YH(NYH,*), WM(*), IWM(*), EWT(*),
2191 1 SAVF(*), FTEM(*), PAR(*)
2192 EXTERNAL F, JAC
2193 PARAMETER (ZERO=0.0D0,ONE=1.0D0)
2194 COMMON /ODE001/ ROWND, ROWNS(173),
2195 2 RDUM1(37), EL0, H, RDUM2(4), TN, UROUND,
2196 3 IOWND(14), IOWNS(4),
2197 4 IDUM1(3), IERPJ, IDUM2, JCUR, IDUM3(4),
2198 5 MITER, IDUM4(4), N, IDUM5(2), NFE, NJE, IDUM6
2199 C-----------------------------------------------------------------------
2200 C PREPJ IS CALLED BY STODE TO COMPUTE AND PROCESS THE MATRIX
2201 C P = I - H*EL(1)*J , WHERE J IS AN APPROXIMATION TO THE JACOBIAN.
2202 C IF ISOPT = 1, PREPJ IS ALSO CALLED BY SPRIME WITH JOPT = 1.
2203 C HERE J IS COMPUTED BY THE USER-SUPPLIED ROUTINE JAC IF
2204 C MITER = 1 OR 4, OR BY FINITE DIFFERENCING IF MITER = 2, 3, OR 5.
2205 C IF MITER = 3, A DIAGONAL APPROXIMATION TO J IS USED.
2206 C J IS STORED IN WM AND REPLACED BY P. IF MITER .NE. 3, P IS THEN
2207 C SUBJECTED TO LU DECOMPOSITION (JOPT = 0) IN PREPARATION FOR LATER
2208 C SOLUTION OF LINEAR SYSTEMS WITH P AS COEFFICIENT MATRIX. THIS IS
2209 C DONE BY DGEFA IF MITER = 1 OR 2, AND BY DGBFA IF MITER = 4 OR 5.
2211 C IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION
2212 C WITH PREPJ USES THE FOLLOWING..
2213 C Y = ARRAY CONTAINING PREDICTED VALUES ON ENTRY.
2214 C FTEM = WORK ARRAY OF LENGTH N (ACOR IN STODE).
2215 C SAVF = ARRAY CONTAINING F EVALUATED AT PREDICTED Y.
2216 C WM = REAL WORK SPACE FOR MATRICES. ON OUTPUT IT CONTAINS THE
2217 C INVERSE DIAGONAL MATRIX IF MITER = 3 AND THE LU DECOMPOSITION
2218 C OF P IF MITER IS 1, 2 , 4, OR 5.
2219 C STORAGE OF MATRIX ELEMENTS STARTS AT WM(3).
2220 C WM ALSO CONTAINS THE FOLLOWING MATRIX-RELATED DATA..
2221 C WM(1) = SQRT(UROUND), USED IN NUMERICAL JACOBIAN INCREMENTS.
2222 C WM(2) = H*EL0, SAVED FOR LATER USE IF MITER = 3.
2223 C IWM = INTEGER WORK SPACE CONTAINING PIVOT INFORMATION, STARTING AT
2224 C IWM(21), IF MITER IS 1, 2, 4, OR 5. IWM ALSO CONTAINS BAND
2225 C PARAMETERS ML = IWM(1) AND MU = IWM(2) IF MITER IS 4 OR 5.
2226 C EL0 = EL(1) (INPUT).
2227 C IERPJ = OUTPUT ERROR FLAG, = 0 IF NO TROUBLE, .GT. 0 IF
2228 C P MATRIX FOUND TO BE SINGULAR.
2229 C JCUR = OUTPUT FLAG = 1 TO INDICATE THAT THE JACOBIAN MATRIX
2230 C (OR APPROXIMATION) IS NOW CURRENT.
2231 C JOPT = INPUT JACOBIAN OPTION, = 1 IF JAC IS DESIRED ONLY.
2232 C THIS ROUTINE ALSO USES THE COMMON VARIABLES EL0, H, TN, UROUND,
2233 C IERPJ, JCUR, MITER, N, NFE, AND NJE.
2234 C-----------------------------------------------------------------------
2235 NJE = NJE + 1
2236 IERPJ = 0
2237 JCUR = 1
2238 HL0 = H*EL0
2239 GO TO (100, 200, 300, 400, 500), MITER
2240 C IF MITER = 1, CALL JAC AND MULTIPLY BY SCALAR. -----------------------
2241 100 LENP = N*N
2242 DO 110 I = 1,LENP
2243 110 WM(I+2) = ZERO
2244 CALL JAC (NEQ, TN, Y, PAR, 0, 0, WM(3), N)
2245 IF (JOPT .EQ. 1) RETURN
2246 CON = -HL0
2247 DO 120 I = 1,LENP
2248 120 WM(I+2) = WM(I+2)*CON
2249 GO TO 240
2250 C IF MITER = 2, MAKE N CALLS TO F TO APPROXIMATE J. --------------------
2251 200 FAC = VNORM (N, SAVF, EWT)
2252 R0 = 1000.0D0*ABS(H)*UROUND*REAL(N)*FAC
2253 IF (R0 .EQ. ZERO) R0 = ONE
2254 SRUR = WM(1)
2255 J1 = 2
2256 DO 230 J = 1,N
2257 YJ = Y(J)
2258 R = MAX(SRUR*ABS(YJ),R0/EWT(J))
2259 Y(J) = Y(J) + R
2260 FAC = -HL0/R
2261 CALL F (NEQ, TN, Y, PAR, FTEM)
2262 DO 220 I = 1,N
2263 220 WM(I+J1) = (FTEM(I) - SAVF(I))*FAC
2264 Y(J) = YJ
2265 J1 = J1 + N
2266 230 CONTINUE
2267 NFE = NFE + N
2268 IF (JOPT .EQ. 1) RETURN
2269 C ADD IDENTITY MATRIX. -------------------------------------------------
2270 240 J = 3
2271 DO 250 I = 1,N
2272 WM(J) = WM(J) + ONE
2273 250 J = J + (N + 1)
2274 C DO LU DECOMPOSITION ON P. --------------------------------------------
2275 CALL DGEFA (WM(3), N, N, IWM(21), IER)
2276 IF (IER .NE. 0) IERPJ = 1
2277 RETURN
2278 C IF MITER = 3, CONSTRUCT A DIAGONAL APPROXIMATION TO J AND P. ---------
2279 300 WM(2) = HL0
2280 R = EL0*0.1D0
2281 DO 310 I = 1,N
2282 310 Y(I) = Y(I) + R*(H*SAVF(I) - YH(I,2))
2283 CALL F (NEQ, TN, Y, PAR, WM(3))
2284 NFE = NFE + 1
2285 DO 320 I = 1,N
2286 R0 = H*SAVF(I) - YH(I,2)
2287 DI = 0.1D0*R0 - H*(WM(I+2) - SAVF(I))
2288 WM(I+2) = 1.0D0
2289 IF (ABS(R0) .LT. UROUND/EWT(I)) GO TO 320
2290 IF (ABS(DI) .EQ. ZERO) GO TO 330
2291 WM(I+2) = 0.1D0*R0/DI
2292 320 CONTINUE
2293 RETURN
2294 330 IERPJ = 1
2295 RETURN
2296 C IF MITER = 4, CALL JAC AND MULTIPLY BY SCALAR. -----------------------
2297 400 ML = IWM(1)
2298 MU = IWM(2)
2299 ML3 = ML + 3
2300 MBAND = ML + MU + 1
2301 MEBAND = MBAND + ML
2302 LENP = MEBAND*N
2303 DO 410 I = 1,LENP
2304 410 WM(I+2) = ZERO
2305 CALL JAC (NEQ, TN, Y, PAR, ML, MU, WM(ML3), MEBAND)
2306 IF (JOPT .EQ. 1) RETURN
2307 CON = -HL0
2308 DO 420 I = 1,LENP
2309 420 WM(I+2) = WM(I+2)*CON
2310 GO TO 570
2311 C IF MITER = 5, MAKE MBAND CALLS TO F TO APPROXIMATE J. ----------------
2312 500 ML = IWM(1)
2313 MU = IWM(2)
2314 MBAND = ML + MU + 1
2315 MBA = MIN(MBAND,N)
2316 MEBAND = MBAND + ML
2317 MEB1 = MEBAND - 1
2318 SRUR = WM(1)
2319 FAC = VNORM (N, SAVF, EWT)
2320 R0 = 1000.0D0*ABS(H)*UROUND*REAL(N)*FAC
2321 IF (R0 .EQ. ZERO) R0 = ONE
2322 DO 560 J = 1,MBA
2323 DO 530 I = J,N,MBAND
2324 YI = Y(I)
2325 R = MAX(SRUR*ABS(YI),R0/EWT(I))
2326 530 Y(I) = Y(I) + R
2327 CALL F (NEQ, TN, Y, PAR, FTEM)
2328 DO 550 JJ = J,N,MBAND
2329 Y(JJ) = YH(JJ,1)
2330 YJJ = Y(JJ)
2331 R = MAX(SRUR*ABS(YJJ),R0/EWT(JJ))
2332 FAC = -HL0/R
2333 I1 = MAX(JJ-MU,1)
2334 I2 = MIN(JJ+ML,N)
2335 II = JJ*MEB1 - ML + 2
2336 DO 540 I = I1,I2
2337 540 WM(II+I) = (FTEM(I) - SAVF(I))*FAC
2338 550 CONTINUE
2339 560 CONTINUE
2340 NFE = NFE + MBA
2341 IF (JOPT .EQ. 1) RETURN
2342 C ADD IDENTITY MATRIX. -------------------------------------------------
2343 570 II = MBAND + 2
2344 DO 580 I = 1,N
2345 WM(II) = WM(II) + ONE
2346 580 II = II + MEBAND
2347 C DO LU DECOMPOSITION OF P. --------------------------------------------
2348 CALL DGBFA (WM(3), MEBAND, N, ML, MU, IWM(21), IER)
2349 IF (IER .NE. 0) IERPJ = 1
2350 RETURN
2351 C----------------------- END OF SUBROUTINE PREPJ -----------------------
2353 SUBROUTINE PREPDF (NEQ, Y, SRUR, SAVF, FTEM, DFDP, PAR,
2354 1 F, DF, JPAR)
2355 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2356 EXTERNAL F, DF
2357 DIMENSION NEQ(*), Y(*), SAVF(*), FTEM(*), DFDP(*), PAR(*)
2358 COMMON /ODE001/ ROWND, ROWNS(173),
2359 1 RDUM1(43), TN, RDUM2,
2360 2 IOWND1(14), IOWNS(4),
2361 3 IDUM1(10), MITER, IDUM2(4), N, IDUM3(2), NFE, IDUM4(2)
2362 COMMON /ODE002/ RDUM3(3),
2363 1 IOWND2(3), IDUM5(2), NDFE, IDUM6, IDF, IDUM7(3)
2364 C-----------------------------------------------------------------------
2365 C PREPDF IS CALLED BY SPRIME AND STESA TO COMPUTE THE INHOMOGENEITY
2366 C VECTORS DF(I)/DP(JPAR). HERE DF/DP IS COMPUTED BY THE USER-SUPPLIED
2367 C ROUTINE DF IF IDF = 1, OR BY FINITE DIFFERENCING IF IDF = 0.
2369 C IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION WITH
2370 C PREPDF USES THE FOLLOWING..
2371 C Y = REAL ARRAY OF LENGTH NYH CONTAINING DEPENDENT VARIABLES.
2372 C PREPDF USES ONLY THE FIRST N ENTRIES OF Y(*).
2373 C SRUR = SQRT(UROUND) (= WM(1)).
2374 C SAVF = REAL ARRAY OF LENGTH N CONTAINING DERIVATIVES DY/DT.
2375 C FTEM = REAL ARRAY OF LENGTH N USED TO TEMPORARILY STORE DY/DT FOR
2376 C NUMERICAL DIFFERENTIATION.
2377 C DFDP = REAL ARRAY OF LENGTH N USED TO STORE DF(I)/DP(JPAR), I = 1,N.
2378 C PAR = REAL ARRAY OF LENGTH NPAR CONTAINING EQUATION PARAMETERS
2379 C OF INTEREST.
2380 C JPAR = INPUT PARAMETER, 2 .LE. JPAR .LE. NSV, DESIGNATING THE
2381 C APPROPRIATE SOLUTION VECTOR CORRESPONDING TO PAR(JPAR).
2382 C THIS ROUTINE ALSO USES THE COMMON VARIABLES TN, MITER, N, NFE, NDFE,
2383 C AND IDF.
2384 C-----------------------------------------------------------------------
2385 NDFE = NDFE + 1
2386 IDF1 = IDF + 1
2387 GO TO (100, 200), IDF1
2388 C IDF = 0, CALL F TO APPROXIMATE DFDP. ---------------------------------
2389 100 RPAR = PAR(JPAR)
2390 R = MAX(SRUR*ABS(RPAR),SRUR)
2391 PAR(JPAR) = RPAR + R
2392 FAC = 1.0D0/R
2393 CALL F (NEQ, TN, Y, PAR, FTEM)
2394 DO 110 I = 1,N
2395 110 DFDP(I) = (FTEM(I) - SAVF(I))*FAC
2396 PAR(JPAR) = RPAR
2397 NFE = NFE + 1
2398 RETURN
2399 C IDF = 1, CALL USER SUPPLIED DF. --------------------------------------
2400 200 DO 210 I = 1,N
2401 210 DFDP(I) = 0.0D0
2402 CALL DF (NEQ, TN, Y, PAR, DFDP, JPAR)
2403 RETURN
2404 C -------------------- END OF SUBROUTINE PREPDF ------------------------
2406 SUBROUTINE INTDY (T, K, YH, NYH, DKY, IFLAG)
2407 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2408 DIMENSION YH(NYH,1), DKY(1)
2409 COMMON /ODE001/ ROWND, ROWNS(173),
2410 2 RDUM1(38),H, RDUM2(2), HU, RDUM3, TN, UROUND,
2411 3 IOWND(14), IOWNS(4),
2412 4 IDUM1(8), L, IDUM2,
2413 5 IDUM3(5), N, NQ, IDUM4(4)
2414 C-----------------------------------------------------------------------
2415 C INTDY COMPUTES INTERPOLATED VALUES OF THE K-TH DERIVATIVE OF THE
2416 C DEPENDENT VARIABLE VECTOR Y, AND STORES IT IN DKY. THIS ROUTINE
2417 C IS CALLED WITHIN THE PACKAGE WITH K = 0 AND T = TOUT, BUT MAY
2418 C ALSO BE CALLED BY THE USER FOR ANY K UP TO THE CURRENT ORDER.
2419 C (SEE DETAILED INSTRUCTIONS IN THE USAGE DOCUMENTATION.)
2420 C-----------------------------------------------------------------------
2421 C THE COMPUTED VALUES IN DKY ARE GOTTEN BY INTERPOLATION USING THE
2422 C NORDSIECK HISTORY ARRAY YH. THIS ARRAY CORRESPONDS UNIQUELY TO A
2423 C VECTOR-VALUED POLYNOMIAL OF DEGREE NQCUR OR LESS, AND DKY IS SET
2424 C TO THE K-TH DERIVATIVE OF THIS POLYNOMIAL AT T.
2425 C THE FORMULA FOR DKY IS..
2427 C DKY(I) = SUM C(J,K) * (T - TN)**(J-K) * H**(-J) * YH(I,J+1)
2428 C J=K
2429 C WHERE C(J,K) = J*(J-1)*...*(J-K+1), Q = NQCUR, TN = TCUR, H = HCUR.
2430 C THE QUANTITIES NQ = NQCUR, L = NQ+1, N = NEQ, TN, AND H ARE
2431 C COMMUNICATED BY COMMON. THE ABOVE SUM IS DONE IN REVERSE ORDER.
2432 C IFLAG IS RETURNED NEGATIVE IF EITHER K OR T IS OUT OF BOUNDS.
2433 C-----------------------------------------------------------------------
2434 IFLAG = 0
2435 IF (K .LT. 0 .OR. K .GT. NQ) GO TO 80
2436 TP = TN - HU*(1.0D0 + 100.0D0*UROUND)
2437 IF ((T-TP)*(T-TN) .GT. 0.0D0) GO TO 90
2439 S = (T - TN)/H
2440 IC = 1
2441 IF (K .EQ. 0) GO TO 15
2442 JJ1 = L - K
2443 DO 10 JJ = JJ1,NQ
2444 10 IC = IC*JJ
2445 15 C = REAL(IC)
2446 DO 20 I = 1,NYH
2447 20 DKY(I) = C*YH(I,L)
2448 IF (K .EQ. NQ) GO TO 55
2449 JB2 = NQ - K
2450 DO 50 JB = 1,JB2
2451 J = NQ - JB
2452 JP1 = J + 1
2453 IC = 1
2454 IF (K .EQ. 0) GO TO 35
2455 JJ1 = JP1 - K
2456 DO 30 JJ = JJ1,J
2457 30 IC = IC*JJ
2458 35 C = REAL(IC)
2459 DO 40 I = 1,NYH
2460 40 DKY(I) = C*YH(I,JP1) + S*DKY(I)
2461 50 CONTINUE
2462 IF (K .EQ. 0) RETURN
2463 55 R = H**(-K)
2464 DO 60 I = 1,NYH
2465 60 DKY(I) = R*DKY(I)
2466 RETURN
2468 80 CALL XERR('INTDY-- K (=I1) ILLEGAL',
2469 1 51, 1, 1, K, 0, 0, ZERO,ZERO)
2470 IFLAG = -1
2471 RETURN
2472 90 CALL XERR ('INTDY-- T (=R1) ILLEGAL',
2473 1 52, 1, 0, 0, 0, 1, T, ZERO)
2474 CALL XERR('T NOT IN INTERVAL TCUR - HU (= R1) TO TCUR (=R2)',
2475 1 52, 1, 0, 0, 0, 2, TP, TN)
2476 IFLAG = -2
2477 RETURN
2478 C----------------------- END OF SUBROUTINE INTDY -----------------------
2480 SUBROUTINE STESA (NEQ, Y, NROW, NCOL, YH, WM, IWM, EWT, SAVF,
2481 1 ACOR, PAR, NRS, F, JAC, DF, PJAC, PDF, KppSolve)
2482 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2483 EXTERNAL F, JAC, DF, PJAC, PDF, KppSolve
2484 DIMENSION NEQ(*), Y(NROW,*), YH(NROW,NCOL,*), WM(*), IWM(*),
2485 1 EWT(NROW,*), SAVF(*), ACOR(NROW,*), PAR(*), NRS(*)
2486 PARAMETER (ONE=1.0D0,ZERO=0.0D0)
2487 COMMON /ODE001/ ROWND, ROWNS(173),
2488 1 TESCO(3,12), RDUM1, EL0, H, RDUM2(4), TN, RDUM3,
2489 2 IOWND1(14), IOWNS(4),
2490 3 IALTH, LMAX, IDUM1, IERPJ, IERSL, JCUR, IDUM2, KFLAG, L, IDUM3,
2491 4 MITER, IDUM4(4), N, NQ, IDUM5, NFE, IDUM6(2)
2492 COMMON /ODE002/ DUPS, DSMS, DDNS,
2493 1 IOWND2(3), IDUM7, NSV, IDUM8(2), IDF, IDUM9, JOPT, KFLAGS
2494 C-----------------------------------------------------------------------
2495 C STESA IS CALLED BY STODE TO PERFORM AN EXPLICIT CALCULATION FOR THE
2496 C FIRST-ORDER SENSITIVITY COEFFICIENTS DY(I)/DP(J), I = 1,N; J = 1,NPAR.
2498 C IN ADDITION TO THE VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION
2499 C WITH STESA USES THE FOLLOWING..
2500 C Y = AN NROW (=N) BY NCOL (=NSV) REAL ARRAY CONTAINING THE
2501 C CORRECTED DEPENDENT VARIABLES ON OUTPUT..
2502 C Y(I,1) , I = 1,N = STATE VARIABLES (INPUT);
2503 C Y(I,J) , I = 1,N , J = 2,NSV ,
2504 C = SENSITIVITY COEFFICIENTS, DY(I)/DP(J).
2505 C YH = AN N BY NSV BY LMAX REAL ARRAY CONTAINING THE PREDICTED
2506 C DEPENDENT VARIABLES AND THEIR APPROXIMATE SCALED DERIVATIVES.
2507 C SAVF = A REAL ARRAY OF LENGTH N USED TO STORE FIRST DERIVATIVES
2508 C OF DEPENDENT VARIABLES IF MITER = 2 OR 5.
2509 C PAR = A REAL ARRAY OF LENGTH NPAR CONTAINING THE EQUATION
2510 C PARAMETERS OF INTEREST.
2511 C NRS = AN INTEGER ARRAY OF LENGTH NPAR + 1 CONTAINING THE NUMBER
2512 C OF REPEATED STEPS (KFLAGS .LT. 0) DUE TO THE SENSITIVITY
2513 C CALCULATIONS..
2514 C NRS(1) = TOTAL NUMBER OF REPEATED STEPS
2515 C NRS(I) , I = 2,NPAR = NUMBER OF REPEATED STEPS DUE
2516 C TO PARAMETER I.
2517 C NSV = NUMBER OF SOLUTION VECTORS = NPAR + 1.
2518 C KFLAGS = LOCAL ERROR TEST FLAG, = 0 IF TEST PASSES, .LT. 0 IF TEST
2519 C FAILS, AND STEP NEEDS TO BE REPEATED. ERROR TEST IS APPLIED
2520 C TO EACH SOLUTION VECTOR INDEPENDENTLY.
2521 C DUPS, DSMS, DDNS = REAL SCALARS USED FOR COMPUTING RHUP, RHSM, RHDN,
2522 C ON RETURN TO STODE (IALTH .EQ. 1).
2523 C THIS ROUTINE ALSO USES THE COMMON VARIABLES EL0, H, TN, IALTH, LMAX,
2524 C IERPJ, IERSL, JCUR, KFLAG, L, MITER, N, NQ, NFE, AND JOPT.
2525 C-----------------------------------------------------------------------
2526 DUPS = ZERO
2527 DSMS = ZERO
2528 DDNS = ZERO
2529 HL0 = H*EL0
2530 EL0I = ONE/EL0
2531 TI2 = ONE/TESCO(2,NQ)
2532 TI3 = ONE/TESCO(3,NQ)
2533 C IF MITER = 2 OR 5 (OR IDF = 0), SUPPLY DERIVATIVES AT CORRECTED
2534 C Y(*,1) VALUES FOR NUMERICAL DIFFERENTIATION IN PJAC AND/OR PDF.
2535 IF (MITER .EQ. 2 .OR. MITER .EQ. 5 .OR. IDF .EQ. 0) GO TO 10
2536 GO TO 15
2537 10 CALL F (NEQ, TN, Y, PAR, SAVF)
2538 NFE = NFE + 1
2539 C IF JCUR = 0, UPDATE THE JACOBIAN MATRIX.
2540 C IF MITER = 5, LOAD CORRECTED Y(*,1) VALUES INTO Y(*,2).
2541 15 IF (JCUR .EQ. 1) GO TO 30
2542 IF (MITER .NE. 5) GO TO 25
2543 DO 20 I = 1,N
2544 20 Y(I,2) = Y(I,1)
2545 25 CALL PJAC (NEQ, Y, Y(1,2), N, WM, IWM, EWT, SAVF, ACOR(1,2),
2546 1 PAR, F, JAC, JOPT)
2547 IF (IERPJ .NE. 0) RETURN
2548 C-----------------------------------------------------------------------
2549 C THIS IS A LOOPING POINT FOR THE SENSITIVITY CALCULATIONS.
2550 C-----------------------------------------------------------------------
2551 C FOR EACH PARAMETER PAR(*), A SENSITIVITY SOLUTION VECTOR IS COMPUTED
2552 C USING THE SAME STEP SIZE (H) AND ORDER (NQ) AS IN STODE.
2553 C A LOCAL ERROR TEST IS APPLIED INDEPENDENTLY TO EACH SOLUTION VECTOR.
2554 C-----------------------------------------------------------------------
2555 30 DO 100 J = 2,NSV
2556 JPAR = J - 1
2557 C EVALUATE INHOMOGENEITY TERM, TEMPORARILY LOAD INTO Y(*,JPAR+1). ------
2558 CALL PDF(NEQ, Y, WM, SAVF, ACOR(1,J), Y(1,J), PAR,
2559 1 F, DF, JPAR)
2560 C-----------------------------------------------------------------------
2561 C LOAD RHS OF SENSITIVITY SOLUTION (CORRECTOR) EQUATION..
2563 C RHS = DY/DP - EL(1)*H*D(DY/DP)/DT + EL(1)*H*DF/DP
2565 C-----------------------------------------------------------------------
2566 DO 40 I = 1,N
2567 40 Y(I,J) = YH(I,J,1) - EL0*YH(I,J,2) + HL0*Y(I,J)
2568 C-----------------------------------------------------------------------
2569 C KppSolve CORRECTOR EQUATION: THE SOLUTIONS ARE LOCATED IN Y(*,JPAR+1).
2570 C THE EXPLICIT FORMULA IS..
2572 C (I - EL(1)*H*JAC) * DY/DP(CORRECTED) = RHS
2574 C-----------------------------------------------------------------------
2575 CALL KppSolve (WM, IWM, Y(1,J), DUM)
2576 IF (IERSL .NE. 0) RETURN
2577 C ESTIMATE LOCAL TRUNCATION ERROR. -------------------------------------
2578 DO 50 I = 1,N
2579 50 ACOR(I,J) = (Y(I,J) - YH(I,J,1))*EL0I
2580 ERR = VNORM(N, ACOR(1,J), EWT(1,J))*TI2
2581 IF (ERR .GT. ONE) GO TO 200
2582 C-----------------------------------------------------------------------
2583 C LOCAL ERROR TEST PASSED. SET KFLAGS TO 0 TO INDICATE THIS.
2584 C IF IALTH = 1, COMPUTE DSMS, DDNS, AND DUPS (IF L .LT. LMAX).
2585 C-----------------------------------------------------------------------
2586 KFLAGS = 0
2587 IF (IALTH .GT. 1) GO TO 100
2588 IF (L .EQ. LMAX) GO TO 70
2589 DO 60 I= 1,N
2590 60 Y(I,J) = ACOR(I,J) - YH(I,J,LMAX)
2591 DUPS = MAX(DUPS,VNORM(N,Y(1,J),EWT(1,J))*TI3)
2592 70 DSMS = MAX(DSMS,ERR)
2593 100 CONTINUE
2594 RETURN
2595 C-----------------------------------------------------------------------
2596 C THIS SECTION IS REACHED IF THE ERROR TOLERANCE FOR SENSITIVITY
2597 C SOLUTION VECTOR JPAR HAS BEEN VIOLATED. KFLAGS IS MADE NEGATIVE TO
2598 C INDICATE THIS. IF KFLAGS = -1, SET KFLAG EQUAL TO ZERO SO THAT KFLAG
2599 C IS SET TO -1 ON RETURN TO STODE BEFORE REPEATING THE STEP.
2600 C INCREMENT NRS(1) (= TOTAL NUMBER OF REPEATED STEPS DUE TO ALL
2601 C SENSITIVITY SOLUTION VECTORS) BY ONE.
2602 C INCREMENT NRS(JPAR+1) (= TOTAL NUMBER OF REPEATED STEPS DUE TO
2603 C SOLUTION VECTOR JPAR+1) BY ONE.
2604 C LOAD DSMS FOR RH CALCULATION IN STODE.
2605 C-----------------------------------------------------------------------
2606 200 KFLAGS = KFLAGS - 1
2607 IF (KFLAGS .EQ. -1) KFLAG = 0
2608 NRS(1) = NRS(1) + 1
2609 NRS(J) = NRS(J) + 1
2610 DSMS = ERR
2611 RETURN
2612 C------------------------ END OF SUBROUTINE STESA ----------------------
2614 SUBROUTINE STODE (NEQ, Y, YH, NYH, YH1, WM, IWM, EWT, SAVF, ACOR,
2615 1 PAR, NRS, F, JAC, DF, PJAC, PDF, SLVS)
2616 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2617 EXTERNAL F, JAC, DF, PJAC, PDF, SLVS
2618 DIMENSION NEQ(*), Y(*), YH(NYH,*), YH1(*), WM(*), IWM(*), EWT(*),
2619 1 SAVF(*), ACOR(*), PAR(*), NRS(*)
2620 PARAMETER (ONE=1.0D0,ZERO=0.0D0)
2621 COMMON /ODE001/ ROWND,
2622 1 CONIT, CRATE, EL(13), ELCO(13,12), HOLD, RMAX,
2623 2 TESCO(3,12), CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND,
2624 3 IOWND1(14), IPUP, MEO, NQNYH, NSLP,
2625 4 IALTH, LMAX, ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH,
2626 5 MITER, MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
2627 COMMON /ODE002/ DUPS, DSMS, DDNS,
2628 1 IOWND2(3), ISOPT, NSV, NDFE, NSPE, IDF, IERSP, JOPT, KFLAGS
2629 C-----------------------------------------------------------------------
2630 C STODE PERFORMS ONE STEP OF THE INTEGRATION OF AN INITIAL VALUE
2631 C PROBLEM FOR A SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS.
2632 C NOTE.. STODE IS INDEPENDENT OF THE VALUE OF THE ITERATION METHOD
2633 C INDICATOR MITER, WHEN THIS IS .NE. 0, AND HENCE IS INDEPENDENT
2634 C OF THE TYPE OF CHORD METHOD USED, OR THE JACOBIAN STRUCTURE.
2635 C FOR ISOPT = 1, STODE CALLS STESA FOR SENSITIVITY CALCULATIONS.
2636 C VARIABLES USED FOR COMMUNICATION WITH STESA ARE DESCRIBED IN STESA.
2637 C COMMUNICATION WITH STODE IS DONE WITH THE FOLLOWING VARIABLES..
2639 C NEQ = INTEGER ARRAY CONTAINING PROBLEM SIZE IN NEQ(1), AND
2640 C NUMBER OF PARAMETERS TO BE CONSIDERED IN THE SENSITIVITY
2641 C ANALYSIS NEQ(2) (FOR ISOPT = 1), AND PASSED AS THE
2642 C NEQ ARGUMENT IN ALL CALLS TO F, JAC, AND DF.
2643 C Y = AN ARRAY OF LENGTH .GE. N USED AS THE Y ARGUMENT IN
2644 C ALL CALLS TO F, JAC, AND DF.
2645 C YH = AN NYH BY LMAX ARRAY CONTAINING THE DEPENDENT VARIABLES
2646 C AND THEIR APPROXIMATE SCALED DERIVATIVES, WHERE
2647 C LMAX = MAXORD + 1. YH(I,J+1) CONTAINS THE APPROXIMATE
2648 C J-TH DERIVATIVE OF Y(I), SCALED BY H**J/FACTORIAL(J)
2649 C (J = 0,1,...,NQ). ON ENTRY FOR THE FIRST STEP, THE FIRST
2650 C TWO COLUMNS OF YH MUST BE SET FROM THE INITIAL VALUES.
2651 C NYH = A CONSTANT INTEGER .GE. N, THE FIRST DIMENSION OF YH.
2652 C THE TOTAL NUMBER OF FIRST-ORDER DIFFERENTIAL EQUATIONS..
2653 C NYH = N, ISOPT = 0,
2654 C NYH = N * (NPAR + 1), ISOPT = 1
2655 C YH1 = A ONE-DIMENSIONAL ARRAY OCCUPYING THE SAME SPACE AS YH.
2656 C EWT = AN ARRAY OF LENGTH NYH CONTAINING MULTIPLICATIVE WEIGHTS
2657 C FOR LOCAL ERROR MEASUREMENTS. LOCAL ERRORS IN Y(I) ARE
2658 C COMPARED TO 1.0/EWT(I) IN VARIOUS ERROR TESTS.
2659 C SAVF = AN ARRAY OF WORKING STORAGE, OF LENGTH N.
2660 C ALSO USED FOR INPUT OF YH(*,MAXORD+2) WHEN JSTART = -1
2661 C AND MAXORD .LT. THE CURRENT ORDER NQ.
2662 C ACOR = A WORK ARRAY OF LENGTH NYH, USED FOR THE ACCUMULATED
2663 C CORRECTIONS. ON A SUCCESSFUL RETURN, ACOR(I) CONTAINS
2664 C THE ESTIMATED ONE-STEP LOCAL ERROR IN Y(I).
2665 C WM,IWM = REAL AND INTEGER WORK ARRAYS ASSOCIATED WITH MATRIX
2666 C OPERATIONS IN CHORD ITERATION (MITER .NE. 0).
2667 C PJAC = NAME OF ROUTINE TO EVALUATE AND PREPROCESS JACOBIAN MATRIX
2668 C AND P = I - H*EL0*JAC, IF A CHORD METHOD IS BEING USED.
2669 C IF ISOPT = 1, PJAC CAN BE CALLED TO CALCULATE JAC BY
2670 C SETTING JOPT = 1.
2671 C SLVS = NAME OF ROUTINE TO KppSolve LINEAR SYSTEM IN CHORD ITERATION.
2672 C CCMAX = MAXIMUM RELATIVE CHANGE IN H*EL0 BEFORE PJAC IS CALLED.
2673 C H = THE STEP SIZE TO BE ATTEMPTED ON THE NEXT STEP.
2674 C H IS ALTERED BY THE ERROR CONTROL ALGORITHM DURING THE
2675 C PROBLEM. H CAN BE EITHER POSITIVE OR NEGATIVE, BUT ITS
2676 C SIGN MUST REMAIN CONSTANT THROUGHOUT THE PROBLEM.
2677 C HMIN = THE MINIMUM ABSOLUTE VALUE OF THE STEP SIZE H TO BE USED.
2678 C HMXI = INVERSE OF THE MAXIMUM ABSOLUTE VALUE OF H TO BE USED.
2679 C HMXI = 0.0 IS ALLOWED AND CORRESPONDS TO AN INFINITE HMAX.
2680 C HMIN AND HMXI MAY BE CHANGED AT ANY TIME, BUT WILL NOT
2681 C TAKE EFFECT UNTIL THE NEXT CHANGE OF H IS CONSIDERED.
2682 C TN = THE INDEPENDENT VARIABLE. TN IS UPDATED ON EACH STEP TAKEN.
2683 C JSTART = AN INTEGER USED FOR INPUT ONLY, WITH THE FOLLOWING
2684 C VALUES AND MEANINGS..
2685 C 0 PERFORM THE FIRST STEP.
2686 C .GT.0 TAKE A NEW STEP CONTINUING FROM THE LAST.
2687 C -1 TAKE THE NEXT STEP WITH A NEW VALUE OF H, MAXORD,
2688 C N, METH, OR MITER.
2689 C -2 TAKE THE NEXT STEP WITH A NEW VALUE OF H,
2690 C BUT WITH OTHER INPUTS UNCHANGED.
2691 C ON RETURN, JSTART IS SET TO 1 TO FACILITATE CONTINUATION.
2692 C KFLAG = A COMPLETION CODE WITH THE FOLLOWING MEANINGS..
2693 C 0 THE STEP WAS SUCCESFUL.
2694 C -1 THE REQUESTED ERROR COULD NOT BE ACHIEVED.
2695 C -2 CORRECTOR CONVERGENCE COULD NOT BE ACHIEVED.
2696 C -3 FATAL ERROR IN PJAC, OR SLVS, (OR STESA).
2697 C A RETURN WITH KFLAG = -1 OR -2 MEANS EITHER
2698 C ABS(H) = HMIN OR 10 CONSECUTIVE FAILURES OCCURRED.
2699 C ON A RETURN WITH KFLAG NEGATIVE, THE VALUES OF TN AND
2700 C THE YH ARRAY ARE AS OF THE BEGINNING OF THE LAST
2701 C STEP, AND H IS THE LAST STEP SIZE ATTEMPTED.
2702 C MAXORD = THE MAXIMUM ORDER OF INTEGRATION METHOD TO BE ALLOWED.
2703 C MAXCOR = THE MAXIMUM NUMBER OF CORRECTOR ITERATIONS ALLOWED.
2704 C (= 3, IF ISOPT = 0)
2705 C (= 4, IF ISOPT = 1)
2706 C MSBP = MAXIMUM NUMBER OF STEPS BETWEEN PJAC CALLS (MITER .GT. 0).
2707 C IF ISOPT = 1, PJAC IS CALLED AT LEAST ONCE EVERY STEP.
2708 C MXNCF = MAXIMUM NUMBER OF CONVERGENCE FAILURES ALLOWED.
2709 C METH/MITER = THE METHOD FLAGS. SEE DESCRIPTION IN DRIVER.
2710 C N = THE NUMBER OF FIRST-ORDER MODEL DIFFERENTIAL EQUATIONS.
2711 C-----------------------------------------------------------------------
2712 KFLAG = 0
2713 KFLAGS = 0
2714 TOLD = TN
2715 NCF = 0
2716 IERPJ = 0
2717 IERSL = 0
2718 JCUR = 0
2719 ICF = 0
2720 IF (JSTART .GT. 0) GO TO 200
2721 IF (JSTART .EQ. -1) GO TO 100
2722 IF (JSTART .EQ. -2) GO TO 160
2723 C-----------------------------------------------------------------------
2724 C ON THE FIRST CALL, THE ORDER IS SET TO 1, AND OTHER VARIABLES ARE
2725 C INITIALIZED. RMAX IS THE MAXIMUM RATIO BY WHICH H CAN BE INCREASED
2726 C IN A SINGLE STEP. IT IS INITIALLY 1.E4 TO COMPENSATE FOR THE SMALL
2727 C INITIAL H, BUT THEN IS NORMALLY EQUAL TO 10. IF A FAILURE
2728 C OCCURS (IN CORRECTOR CONVERGENCE OR ERROR TEST), RMAX IS SET AT 2
2729 C FOR THE NEXT INCREASE.
2730 C THESE COMPUTATIONS CONSIDER ONLY THE ORIGINAL SOLUTION VECTOR.
2731 C THE SENSITIVITY SOLUTION VECTORS ARE CONSIDERED IN STESA (ISOPT = 1).
2732 C-----------------------------------------------------------------------
2733 LMAX = MAXORD + 1
2734 NQ = 1
2735 L = 2
2736 IALTH = 2
2737 RMAX = 10000.0D0
2738 RC = ZERO
2739 EL0 = ONE
2740 CRATE = 0.7D0
2741 DELP = ZERO
2742 HOLD = H
2743 MEO = METH
2744 NSLP = 0
2745 IPUP = MITER
2746 IRET = 3
2747 GO TO 140
2748 C-----------------------------------------------------------------------
2749 C THE FOLLOWING BLOCK HANDLES PRELIMINARIES NEEDED WHEN JSTART = -1.
2750 C IPUP IS SET TO MITER TO FORCE A MATRIX UPDATE.
2751 C IF AN ORDER INCREASE IS ABOUT TO BE CONSIDERED (IALTH = 1),
2752 C IALTH IS RESET TO 2 TO POSTPONE CONSIDERATION ONE MORE STEP.
2753 C IF THE CALLER HAS CHANGED METH, CFODE IS CALLED TO RESET
2754 C THE COEFFICIENTS OF THE METHOD.
2755 C IF THE CALLER HAS CHANGED MAXORD TO A VALUE LESS THAN THE CURRENT
2756 C ORDER NQ, NQ IS REDUCED TO MAXORD, AND A NEW H CHOSEN ACCORDINGLY.
2757 C IF H IS TO BE CHANGED, YH MUST BE RESCALED.
2758 C IF H OR METH IS BEING CHANGED, IALTH IS RESET TO L = NQ + 1
2759 C TO PREVENT FURTHER CHANGES IN H FOR THAT MANY STEPS.
2760 C-----------------------------------------------------------------------
2761 100 IPUP = MITER
2762 LMAX = MAXORD + 1
2763 IF (IALTH .EQ. 1) IALTH = 2
2764 IF (METH .EQ. MEO) GO TO 110
2765 CALL CFODE (METH, ELCO, TESCO)
2766 MEO = METH
2767 IF (NQ .GT. MAXORD) GO TO 120
2768 IALTH = L
2769 IRET = 1
2770 GO TO 150
2771 110 IF (NQ .LE. MAXORD) GO TO 160
2772 120 NQ = MAXORD
2773 L = LMAX
2774 DO 125 I = 1,L
2775 125 EL(I) = ELCO(I,NQ)
2776 NQNYH = NQ*NYH
2777 RC = RC*EL(1)/EL0
2778 EL0 = EL(1)
2779 CONIT = 0.5D0/REAL(NQ+2)
2780 DDN = VNORM (N, SAVF, EWT)/TESCO(1,L)
2781 EXDN = ONE/REAL(L)
2782 RHDN = ONE/(1.3D0*DDN**EXDN + 0.0000013D0)
2783 RH = MIN(RHDN,ONE)
2784 IREDO = 3
2785 IF (H .EQ. HOLD) GO TO 170
2786 RH = MIN(RH,ABS(H/HOLD))
2787 H = HOLD
2788 GO TO 175
2789 C-----------------------------------------------------------------------
2790 C CFODE IS CALLED TO GET ALL THE INTEGRATION COEFFICIENTS FOR THE
2791 C CURRENT METH. THEN THE EL VECTOR AND RELATED CONSTANTS ARE RESET
2792 C WHENEVER THE ORDER NQ IS CHANGED, OR AT THE START OF THE PROBLEM.
2793 C-----------------------------------------------------------------------
2794 140 CALL CFODE (METH, ELCO, TESCO)
2795 150 DO 155 I = 1,L
2796 155 EL(I) = ELCO(I,NQ)
2797 NQNYH = NQ*NYH
2798 RC = RC*EL(1)/EL0
2799 EL0 = EL(1)
2800 CONIT = 0.5D0/REAL(NQ+2)
2801 GO TO (160, 170, 200), IRET
2802 C-----------------------------------------------------------------------
2803 C IF H IS BEING CHANGED, THE H RATIO RH IS CHECKED AGAINST
2804 C RMAX, HMIN, AND HMXI, AND THE YH ARRAY RESCALED. IALTH IS SET TO
2805 C L = NQ + 1 TO PREVENT A CHANGE OF H FOR THAT MANY STEPS, UNLESS
2806 C FORCED BY A CONVERGENCE OR ERROR TEST FAILURE.
2807 C-----------------------------------------------------------------------
2808 160 IF (H .EQ. HOLD) GO TO 200
2809 RH = H/HOLD
2810 H = HOLD
2811 IREDO = 3
2812 GO TO 175
2813 170 RH = MAX(RH,HMIN/ABS(H))
2814 175 RH = MIN(RH,RMAX)
2815 RH = RH/MAX(ONE,ABS(H)*HMXI*RH)
2816 R = ONE
2817 DO 180 J = 2,L
2818 R = R*RH
2819 DO 180 I = 1,NYH
2820 180 YH(I,J) = YH(I,J)*R
2821 H = H*RH
2822 RC = RC*RH
2823 IALTH = L
2824 IF (IREDO .EQ. 0) GO TO 690
2825 C-----------------------------------------------------------------------
2826 C THIS SECTION COMPUTES THE PREDICTED VALUES BY EFFECTIVELY
2827 C MULTIPLYING THE YH ARRAY BY THE PASCAL TRIANGLE MATRIX.
2828 C RC IS THE RATIO OF NEW TO OLD VALUES OF THE COEFFICIENT H*EL(1).
2829 C WHEN RC DIFFERS FROM 1 BY MORE THAN CCMAX, IPUP IS SET TO MITER
2830 C TO FORCE PJAC TO BE CALLED, IF A JACOBIAN IS INVOLVED.
2831 C IN ANY CASE, PJAC IS CALLED AT LEAST EVERY MSBP STEPS FOR ISOPT = 0,
2832 C AND AT LEAST ONCE EVERY STEP FOR ISOPT = 1.
2833 C-----------------------------------------------------------------------
2834 200 IF (ABS(RC-ONE) .GT. CCMAX) IPUP = MITER
2835 IF (NST .GE. NSLP+MSBP) IPUP = MITER
2836 TN = TN + H
2837 I1 = NQNYH + 1
2838 DO 215 JB = 1,NQ
2839 I1 = I1 - NYH
2840 DO 210 I = I1,NQNYH
2841 210 YH1(I) = YH1(I) + YH1(I+NYH)
2842 215 CONTINUE
2843 C-----------------------------------------------------------------------
2844 C UP TO MAXCOR CORRECTOR ITERATIONS ARE TAKEN. (= 3, FOR ISOPT = 0;
2845 C = 4, FOR ISOPT = 1). A CONVERGENCE TEST IS MADE ON THE R.M.S. NORM
2846 C OF EACH CORRECTION, WEIGHTED BY THE ERROR WEIGHT VECTOR EWT. THE SUM
2847 C OF THE CORRECTIONS IS ACCUMULATED IN THE VECTOR ACOR(I), I = 1,N.
2848 C (ACOR(I), I = N+1,NYH IS LOADED IN SUBROUTINE STESA (ISOPT = 1).)
2849 C THE YH ARRAY IS NOT ALTERED IN THE CORRECTOR LOOP.
2850 C-----------------------------------------------------------------------
2851 220 M = 0
2852 DO 230 I = 1,N
2853 230 Y(I) = YH(I,1)
2854 CALL F (NEQ, TN, Y, PAR, SAVF)
2855 NFE = NFE + 1
2856 IF (IPUP .LE. 0) GO TO 250
2857 C-----------------------------------------------------------------------
2858 C IF INDICATED, THE MATRIX P = I - H*EL(1)*J IS REEVALUATED AND
2859 C PREPROCESSED BEFORE STARTING THE CORRECTOR ITERATION. IPUP IS SET
2860 C TO 0 AS AN INDICATOR THAT THIS HAS BEEN DONE.
2861 C-----------------------------------------------------------------------
2862 IPUP = 0
2863 RC = ONE
2864 NSLP = NST
2865 CRATE = 0.7D0
2866 CALL PJAC (NEQ, Y, YH, NYH, WM, IWM, EWT, SAVF, ACOR, PAR,
2867 1 F, JAC, JOPT)
2868 IF (IERPJ .NE. 0) GO TO 430
2869 250 DO 260 I = 1,N
2870 260 ACOR(I) = ZERO
2871 270 IF (MITER .NE. 0) GO TO 350
2872 C-----------------------------------------------------------------------
2873 C IN THE CASE OF FUNCTIONAL ITERATION, UPDATE Y DIRECTLY FROM
2874 C THE RESULT OF THE LAST FUNCTION EVALUATION.
2875 C (IF ISOPT = 1, FUNCTIONAL ITERATION IS NOT ALLOWED.)
2876 C-----------------------------------------------------------------------
2877 DO 290 I = 1,N
2878 SAVF(I) = H*SAVF(I) - YH(I,2)
2879 290 Y(I) = SAVF(I) - ACOR(I)
2880 DEL = VNORM (N, Y, EWT)
2881 DO 300 I = 1,N
2882 Y(I) = YH(I,1) + EL(1)*SAVF(I)
2883 300 ACOR(I) = SAVF(I)
2884 GO TO 400
2885 C-----------------------------------------------------------------------
2886 C IN THE CASE OF THE CHORD METHOD, COMPUTE THE CORRECTOR ERROR,
2887 C AND KppSolve THE LINEAR SYSTEM WITH THAT AS RIGHT-HAND SIDE AND
2888 C P AS COEFFICIENT MATRIX.
2889 C-----------------------------------------------------------------------
2890 350 DO 360 I = 1,N
2891 360 Y(I) = H*SAVF(I) - (YH(I,2) + ACOR(I))
2892 CALL SLVS (WM, IWM, Y, SAVF)
2893 IF (IERSL .LT. 0) GO TO 430
2894 IF (IERSL .GT. 0) GO TO 410
2895 DEL = VNORM (N, Y, EWT)
2896 DO 380 I = 1,N
2897 ACOR(I) = ACOR(I) + Y(I)
2898 380 Y(I) = YH(I,1) + EL(1)*ACOR(I)
2899 C-----------------------------------------------------------------------
2900 C TEST FOR CONVERGENCE. IF M.GT.0, AN ESTIMATE OF THE CONVERGENCE
2901 C RATE CONSTANT IS STORED IN CRATE, AND THIS IS USED IN THE TEST.
2902 C-----------------------------------------------------------------------
2903 400 IF (M .NE. 0) CRATE = MAX(0.2D0*CRATE,DEL/DELP)
2904 DCON = DEL*MIN(ONE,1.5D0*CRATE)/(TESCO(2,NQ)*CONIT)
2905 IF (DCON .LE. ONE) GO TO 450
2906 M = M + 1
2907 IF (M .EQ. MAXCOR) GO TO 410
2908 IF (M .GE. 2 .AND. DEL .GT. 2.0D0*DELP) GO TO 410
2909 DELP = DEL
2910 CALL F (NEQ, TN, Y, PAR, SAVF)
2911 NFE = NFE + 1
2912 GO TO 270
2913 C-----------------------------------------------------------------------
2914 C THE CORRECTOR ITERATION FAILED TO CONVERGE IN MAXCOR TRIES.
2915 C IF MITER .NE. 0 AND THE JACOBIAN IS OUT OF DATE, PJAC IS CALLED FOR
2916 C THE NEXT TRY. OTHERWISE THE YH ARRAY IS RETRACTED TO ITS VALUES
2917 C BEFORE PREDICTION, AND H IS REDUCED, IF POSSIBLE. IF H CANNOT BE
2918 C REDUCED OR MXNCF FAILURES HAVE OCCURRED, EXIT WITH KFLAG = -2.
2919 C-----------------------------------------------------------------------
2920 410 IF (MITER .EQ. 0 .OR. JCUR .EQ. 1) GO TO 430
2921 ICF = 1
2922 IPUP = MITER
2923 GO TO 220
2924 430 ICF = 2
2925 NCF = NCF + 1
2926 RMAX = 2.0D0
2927 TN = TOLD
2928 I1 = NQNYH + 1
2929 DO 445 JB = 1,NQ
2930 I1 = I1 - NYH
2931 DO 440 I = I1,NQNYH
2932 440 YH1(I) = YH1(I) - YH1(I+NYH)
2933 445 CONTINUE
2934 IF (IERPJ .LT. 0 .OR. IERSL .LT. 0) GO TO 680
2935 IF (ABS(H) .LE. HMIN*1.00001D0) GO TO 670
2936 IF (NCF .EQ. MXNCF) GO TO 670
2937 RH = 0.25D0
2938 IPUP = MITER
2939 IREDO = 1
2940 GO TO 170
2941 C-----------------------------------------------------------------------
2942 C THE CORRECTOR HAS CONVERGED.
2943 C THE LOCAL ERROR TEST IS MADE AND CONTROL PASSES TO STATEMENT 500
2944 C IF IT FAILS. OTHERWISE, STESA IS CALLED (ISOPT = 1) TO PERFORM
2945 C SENSITIVITY CALCULATIONS AT CURRENT STEP SIZE AND ORDER.
2946 C-----------------------------------------------------------------------
2947 450 CONTINUE
2948 IF (M .EQ. 0) DSM = DEL/TESCO(2,NQ)
2949 IF (M .GT. 0) DSM = VNORM (N, ACOR, EWT)/TESCO(2,NQ)
2950 IF (DSM .GT. ONE) GO TO 500
2952 IF (ISOPT .EQ. 0) GO TO 460
2953 C-----------------------------------------------------------------------
2954 C CALL STESA TO PERFORM EXPLICIT SENSITIVITY ANALYSIS.
2955 C IF THE LOCAL ERROR TEST FAILS (WITHIN STESA) FOR ANY SOLUTION VECTOR,
2956 C KFLAGS IS SET .LT. 0 AND CONTROL PASSES TO STATEMENT 500 UPON RETURN.
2957 C IN EITHER CASE, JCUR IS SET TO ZERO TO SIGNAL THAT THE JACOBIAN MAY
2958 C NEED UPDATING LATER.
2959 C-----------------------------------------------------------------------
2960 CALL STESA (NEQ, Y, N, NSV, YH, WM, IWM, EWT, SAVF, ACOR,
2961 1 PAR, NRS, F, JAC, DF, PJAC, PDF, SLVS)
2962 IF (IERPJ .NE. 0 .OR. IERSL .NE. 0) GO TO 680
2963 IF (KFLAGS .LT. 0) GO TO 500
2964 C-----------------------------------------------------------------------
2965 C AFTER A SUCCESSFUL STEP, UPDATE THE YH ARRAY.
2966 C CONSIDER CHANGING H IF IALTH = 1. OTHERWISE DECREASE IALTH BY 1.
2967 C IF IALTH IS THEN 1 AND NQ .LT. MAXORD, THEN ACOR IS SAVED FOR
2968 C USE IN A POSSIBLE ORDER INCREASE ON THE NEXT STEP.
2969 C IF A CHANGE IN H IS CONSIDERED, AN INCREASE OR DECREASE IN ORDER
2970 C BY ONE IS CONSIDERED ALSO. A CHANGE IN H IS MADE ONLY IF IT IS BY A
2971 C FACTOR OF AT LEAST 1.1. IF NOT, IALTH IS SET TO 3 TO PREVENT
2972 C TESTING FOR THAT MANY STEPS.
2973 C-----------------------------------------------------------------------
2974 460 JCUR = 0
2975 KFLAG = 0
2976 IREDO = 0
2977 NST = NST + 1
2978 HU = H
2979 NQU = NQ
2980 DO 470 J = 1,L
2981 DO 470 I = 1,NYH
2982 470 YH(I,J) = YH(I,J) + EL(J)*ACOR(I)
2983 IALTH = IALTH - 1
2984 IF (IALTH .EQ. 0) GO TO 520
2985 IF (IALTH .GT. 1) GO TO 700
2986 IF (L .EQ. LMAX) GO TO 700
2987 DO 490 I = 1,NYH
2988 490 YH(I,LMAX) = ACOR(I)
2989 GO TO 700
2990 C-----------------------------------------------------------------------
2991 C THE ERROR TEST FAILED IN EITHER STODE OR STESA.
2992 C KFLAG KEEPS TRACK OF MULTIPLE FAILURES.
2993 C RESTORE TN AND THE YH ARRAY TO THEIR PREVIOUS VALUES, AND PREPARE
2994 C TO TRY THE STEP AGAIN. COMPUTE THE OPTIMUM STEP SIZE FOR THIS OR
2995 C ONE LOWER ORDER. AFTER 2 OR MORE FAILURES, H IS FORCED TO DECREASE
2996 C BY A FACTOR OF 0.2 OR LESS.
2997 C-----------------------------------------------------------------------
2998 500 KFLAG = KFLAG - 1
2999 JCUR = 0
3000 TN = TOLD
3001 I1 = NQNYH + 1
3002 DO 515 JB = 1,NQ
3003 I1 = I1 - NYH
3004 DO 510 I = I1,NQNYH
3005 510 YH1(I) = YH1(I) - YH1(I+NYH)
3006 515 CONTINUE
3007 RMAX = 2.0D0
3008 IF (ABS(H) .LE. HMIN*1.00001D0) GO TO 660
3009 IF (KFLAG .LE. -3) GO TO 640
3010 IREDO = 2
3011 RHUP = ZERO
3012 GO TO 540
3013 C-----------------------------------------------------------------------
3015 C REGARDLESS OF THE SUCCESS OR FAILURE OF THE STEP, FACTORS
3016 C RHDN, RHSM, AND RHUP ARE COMPUTED, BY WHICH H COULD BE MULTIPLIED
3017 C AT ORDER NQ - 1, ORDER NQ, OR ORDER NQ + 1, RESPECTIVELY.
3018 C IN THE CASE OF FAILURE, RHUP = 0.0 TO AVOID AN ORDER INCREASE.
3019 C THE LARGEST OF THESE IS DETERMINED AND THE NEW ORDER CHOSEN
3020 C ACCORDINGLY. IF THE ORDER IS TO BE INCREASED, WE COMPUTE ONE
3021 C ADDITIONAL SCALED DERIVATIVE.
3022 C FOR ISOPT = 1, DUPS AND DSMS ARE LOADED WITH THE LARGEST RMS-NORMS
3023 C OBTAINED BY CONSIDERING SEPARATELY THE SENSITIVITY SOLUTION VECTORS.
3024 C-----------------------------------------------------------------------
3025 520 RHUP = ZERO
3026 IF (L .EQ. LMAX) GO TO 540
3027 DO 530 I = 1,N
3028 530 SAVF(I) = ACOR(I) - YH(I,LMAX)
3029 DUP = VNORM (N, SAVF, EWT)/TESCO(3,NQ)
3030 DUP = MAX(DUP,DUPS)
3031 EXUP = ONE/REAL(L+1)
3032 RHUP = ONE/(1.4D0*DUP**EXUP + 0.0000014D0)
3033 540 EXSM = ONE/REAL(L)
3034 DSM = MAX(DSM,DSMS)
3035 RHSM = ONE/(1.2D0*DSM**EXSM + 0.0000012D0)
3036 RHDN = ZERO
3037 IF (NQ .EQ. 1) GO TO 560
3038 JPOINT = 1
3039 DO 550 J = 1,NSV
3040 DDN = VNORM (N, YH(JPOINT,L), EWT(JPOINT))/TESCO(1,NQ)
3041 DDNS = MAX(DDNS,DDN)
3042 JPOINT = JPOINT + N
3043 550 CONTINUE
3044 DDN = DDNS
3045 DDNS = ZERO
3046 EXDN = ONE/REAL(NQ)
3047 RHDN = ONE/(1.3D0*DDN**EXDN + 0.0000013D0)
3048 560 IF (RHSM .GE. RHUP) GO TO 570
3049 IF (RHUP .GT. RHDN) GO TO 590
3050 GO TO 580
3051 570 IF (RHSM .LT. RHDN) GO TO 580
3052 NEWQ = NQ
3053 RH = RHSM
3054 GO TO 620
3055 580 NEWQ = NQ - 1
3056 RH = RHDN
3057 IF (KFLAG .LT. 0 .AND. RH .GT. ONE) RH = ONE
3058 GO TO 620
3059 590 NEWQ = L
3060 RH = RHUP
3061 IF (RH .LT. 1.1D0) GO TO 610
3062 R = EL(L)/REAL(L)
3063 DO 600 I = 1,NYH
3064 600 YH(I,NEWQ+1) = ACOR(I)*R
3065 GO TO 630
3066 610 IALTH = 3
3067 GO TO 700
3068 620 IF ((KFLAG .EQ. 0) .AND. (RH .LT. 1.1D0)) GO TO 610
3069 IF (KFLAG .LE. -2) RH = MIN(RH,0.2D0)
3070 C-----------------------------------------------------------------------
3071 C IF THERE IS A CHANGE OF ORDER, RESET NQ, L, AND THE COEFFICIENTS.
3072 C IN ANY CASE H IS RESET ACCORDING TO RH AND THE YH ARRAY IS RESCALED.
3073 C THEN EXIT FROM 690 IF THE STEP WAS OK, OR REDO THE STEP OTHERWISE.
3074 C-----------------------------------------------------------------------
3075 IF (NEWQ .EQ. NQ) GO TO 170
3076 630 NQ = NEWQ
3077 L = NQ + 1
3078 IRET = 2
3079 GO TO 150
3080 C-----------------------------------------------------------------------
3081 C CONTROL REACHES THIS SECTION IF 3 OR MORE FAILURES HAVE OCCURED.
3082 C IF 10 FAILURES HAVE OCCURRED, EXIT WITH KFLAG = -1.
3083 C IT IS ASSUMED THAT THE DERIVATIVES THAT HAVE ACCUMULATED IN THE
3084 C YH ARRAY HAVE ERRORS OF THE WRONG ORDER. HENCE THE FIRST
3085 C DERIVATIVE IS RECOMPUTED, AND THE ORDER IS SET TO 1. THEN
3086 C H IS REDUCED BY A FACTOR OF 10, AND THE STEP IS RETRIED,
3087 C UNTIL IT SUCCEEDS OR H REACHES HMIN.
3088 C-----------------------------------------------------------------------
3089 640 IF (KFLAG .EQ. -10) GO TO 660
3090 RH = 0.1D0
3091 RH = MAX(HMIN/ABS(H),RH)
3092 H = H*RH
3093 DO 645 I = 1,NYH
3094 645 Y(I) = YH(I,1)
3095 CALL F (NEQ, TN, Y, PAR, SAVF)
3096 NFE = NFE + 1
3097 IF (ISOPT .EQ. 0) GO TO 649
3098 CALL SPRIME (NEQ, Y, YH, NYH, N, NSV, WM, IWM, EWT, SAVF, ACOR,
3099 1 ACOR(N+1), PAR, F, JAC, DF, PJAC, PDF)
3100 IF (IERSP .LT. 0) GO TO 680
3101 DO 646 I = N+1,NYH
3102 646 YH(I,2) = H*YH(I,2)
3103 649 DO 650 I = 1,N
3104 650 YH(I,2) = H*SAVF(I)
3105 IPUP = MITER
3106 IALTH = 5
3107 IF (NQ .EQ. 1) GO TO 200
3108 NQ = 1
3109 L = 2
3110 IRET = 3
3111 GO TO 150
3112 C-----------------------------------------------------------------------
3113 C ALL RETURNS ARE MADE THROUGH THIS SECTION. H IS SAVED IN HOLD
3114 C TO ALLOW THE CALLER TO CHANGE H ON THE NEXT STEP.
3115 C-----------------------------------------------------------------------
3116 660 KFLAG = -1
3117 GO TO 720
3118 670 KFLAG = -2
3119 GO TO 720
3120 680 KFLAG = -3
3121 GO TO 720
3122 690 RMAX = 10.0D0
3123 700 R = ONE/TESCO(2,NQU)
3124 DO 710 I = 1,NYH
3125 710 ACOR(I) = ACOR(I)*R
3126 720 HOLD = H
3127 JSTART = 1
3128 RETURN
3129 C----------------------- END OF SUBROUTINE STODE -----------------------
3131 SUBROUTINE CFODE (METH, ELCO, TESCO)
3132 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3133 DIMENSION ELCO(13,12), TESCO(3,12)
3134 C-----------------------------------------------------------------------
3135 C CFODE IS CALLED BY THE INTEGRATOR ROUTINE TO SET COEFFICIENTS
3136 C NEEDED THERE. THE COEFFICIENTS FOR THE CURRENT METHOD, AS
3137 C GIVEN BY THE VALUE OF METH, ARE SET FOR ALL ORDERS AND SAVED.
3138 C THE MAXIMUM ORDER ASSUMED HERE IS 12 IF METH = 1 AND 5 IF METH = 2.
3139 C (A SMALLER VALUE OF THE MAXIMUM ORDER IS ALSO ALLOWED.)
3140 C CFODE IS CALLED ONCE AT THE BEGINNING OF THE PROBLEM,
3141 C AND IS NOT CALLED AGAIN UNLESS AND UNTIL METH IS CHANGED.
3143 C THE ELCO ARRAY CONTAINS THE BASIC METHOD COEFFICIENTS.
3144 C THE COEFFICIENTS EL(I), 1 .LE. I .LE. NQ+1, FOR THE METHOD OF
3145 C ORDER NQ ARE STORED IN ELCO(I,NQ). THEY ARE GIVEN BY A GENETRATING
3146 C POLYNOMIAL, I.E.,
3147 C L(X) = EL(1) + EL(2)*X + ... + EL(NQ+1)*X**NQ.
3148 C FOR THE IMPLICIT ADAMS METHODS, L(X) IS GIVEN BY
3149 C DL/DX = (X+1)*(X+2)*...*(X+NQ-1)/FACTORIAL(NQ-1), L(-1) = 0.
3150 C FOR THE BDF METHODS, L(X) IS GIVEN BY
3151 C L(X) = (X+1)*(X+2)* ... *(X+NQ)/K,
3152 C WHERE K = FACTORIAL(NQ)*(1 + 1/2 + ... + 1/NQ).
3154 C THE TESCO ARRAY CONTAINS TEST CONSTANTS USED FOR THE
3155 C LOCAL ERROR TEST AND THE SELECTION OF STEP SIZE AND/OR ORDER.
3156 C AT ORDER NQ, TESCO(K,NQ) IS USED FOR THE SELECTION OF STEP
3157 C SIZE AT ORDER NQ - 1 IF K = 1, AT ORDER NQ IF K = 2, AND AT ORDER
3158 C NQ + 1 IF K = 3.
3159 C-----------------------------------------------------------------------
3160 DIMENSION PC(12)
3161 PARAMETER (ONE=1.0D0,ZERO=0.0D0)
3163 GO TO (100, 200), METH
3165 100 ELCO(1,1) = ONE
3166 ELCO(2,1) = ONE
3167 TESCO(1,1) = ZERO
3168 TESCO(2,1) = 2.0D0
3169 TESCO(1,2) = ONE
3170 TESCO(3,12) = ZERO
3171 PC(1) = ONE
3172 RQFAC = ONE
3173 DO 140 NQ = 2,12
3174 C-----------------------------------------------------------------------
3175 C THE PC ARRAY WILL CONTAIN THE COEFFICIENTS OF THE POLYNOMIAL
3176 C P(X) = (X+1)*(X+2)*...*(X+NQ-1).
3177 C INITIALLY, P(X) = 1.
3178 C-----------------------------------------------------------------------
3179 RQ1FAC = RQFAC
3180 RQFAC = RQFAC/REAL(NQ)
3181 NQM1 = NQ - 1
3182 FNQM1 = REAL(NQM1)
3183 NQP1 = NQ + 1
3184 C FORM COEFFICIENTS OF P(X)*(X+NQ-1). ----------------------------------
3185 PC(NQ) = ZERO
3186 DO 110 IB = 1,NQM1
3187 I = NQP1 - IB
3188 110 PC(I) = PC(I-1) + FNQM1*PC(I)
3189 PC(1) = FNQM1*PC(1)
3190 C COMPUTE INTEGRAL, -1 TO 0, OF P(X) AND X*P(X). -----------------------
3191 PINT = PC(1)
3192 XPIN = PC(1)/2.0D0
3193 TSIGN = ONE
3194 DO 120 I = 2,NQ
3195 TSIGN = -TSIGN
3196 PINT = PINT + TSIGN*PC(I)/REAL(I)
3197 120 XPIN = XPIN + TSIGN*PC(I)/REAL(I+1)
3198 C STORE COEFFICIENTS IN ELCO AND TESCO. --------------------------------
3199 ELCO(1,NQ) = PINT*RQ1FAC
3200 ELCO(2,NQ) = ONE
3201 DO 130 I = 2,NQ
3202 130 ELCO(I+1,NQ) = RQ1FAC*PC(I)/REAL(I)
3203 AGAMQ = RQFAC*XPIN
3204 RAGQ = ONE/AGAMQ
3205 TESCO(2,NQ) = RAGQ
3206 IF (NQ .LT. 12) TESCO(1,NQP1) = RAGQ*RQFAC/REAL(NQP1)
3207 TESCO(3,NQM1) = RAGQ
3208 140 CONTINUE
3209 RETURN
3211 200 PC(1) = ONE
3212 RQ1FAC = ONE
3213 DO 230 NQ = 1,5
3214 C-----------------------------------------------------------------------
3215 C THE PC ARRAY WILL CONTAIN THE COEFFICIENTS OF THE POLYNOMIAL
3216 C P(X) = (X+1)*(X+2)*...*(X+NQ).
3217 C INITIALLY, P(X) = 1.
3218 C-----------------------------------------------------------------------
3219 FNQ = REAL(NQ)
3220 NQP1 = NQ + 1
3221 C FORM COEFFICIENTS OF P(X)*(X+NQ). ------------------------------------
3222 PC(NQP1) = ZERO
3223 DO 210 IB = 1,NQ
3224 I = NQ + 2 - IB
3225 210 PC(I) = PC(I-1) + FNQ*PC(I)
3226 PC(1) = FNQ*PC(1)
3227 C STORE COEFFICIENTS IN ELCO AND TESCO. --------------------------------
3228 DO 220 I = 1,NQP1
3229 220 ELCO(I,NQ) = PC(I)/PC(2)
3230 ELCO(2,NQ) = ONE
3231 TESCO(1,NQ) = RQ1FAC
3232 TESCO(2,NQ) = REAL(NQP1)/ELCO(1,NQ)
3233 TESCO(3,NQ) = REAL(NQ+2)/ELCO(1,NQ)
3234 RQ1FAC = RQ1FAC/FNQ
3235 230 CONTINUE
3236 RETURN
3237 C----------------------- END OF SUBROUTINE CFODE -----------------------
3239 SUBROUTINE SOLSY (WM, IWM, X, TEM)
3240 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3241 DIMENSION WM(*), IWM(*), X(*), TEM(*)
3242 PARAMETER (ZERO=0.0D0,ONE=1.0D0)
3243 COMMON /ODE001/ ROWND, ROWNS(173),
3244 2 RDUM1(37), EL0, H, RDUM2(6),
3245 3 IOWND(14), IOWNS(4),
3246 4 IDUM1(4), IERSL, IDUM2(5),
3247 5 MITER, IDUM3(4), N, IDUM4(5)
3248 C-----------------------------------------------------------------------
3249 C THIS ROUTINE MANAGES THE SOLUTION OF THE LINEAR SYSTEM ARISING FROM
3250 C A CHORD ITERATION. IT IS CALLED IF MITER .NE. 0.
3251 C IF MITER IS 1 OR 2, IT CALLS DGESL TO ACCOMPLISH THIS.
3252 C IF MITER = 3 IT UPDATES THE COEFFICIENT H*EL0 IN THE DIAGONAL
3253 C MATRIX, AND THEN COMPUTES THE SOLUTION.
3254 C IF MITER IS 4 OR 5, IT CALLS DGBSL.
3255 C COMMUNICATION WITH SOLSY USES THE FOLLOWING VARIABLES..
3256 C WM = REAL WORK SPACE CONTAINING THE INVERSE DIAGONAL MATRIX IF
3257 C MITER = 3 AND THE LU DECOMPOSITION OF THE MATRIX OTHERWISE.
3258 C STORAGE OF MATRIX ELEMENTS STARTS AT WM(3).
3259 C WM ALSO CONTAINS THE FOLLOWING MATRIX-RELATED DATA..
3260 C WM(1) = SQRT(UROUND) (NOT USED HERE),
3261 C WM(2) = HL0, THE PREVIOUS VALUE OF H*EL0, USED IF MITER = 3.
3262 C IWM = INTEGER WORK SPACE CONTAINING PIVOT INFORMATION, STARTING AT
3263 C IWM(21), IF MITER IS 1, 2, 4, OR 5. IWM ALSO CONTAINS BAND
3264 C PARAMETERS ML = IWM(1) AND MU = IWM(2) IF MITER IS 4 OR 5.
3265 C X = THE RIGHT-HAND SIDE VECTOR ON INPUT, AND THE SOLUTION VECTOR
3266 C ON OUTPUT, OF LENGTH N.
3267 C TEM = VECTOR OF WORK SPACE OF LENGTH N, NOT USED IN THIS VERSION.
3268 C IERSL = OUTPUT FLAG (IN COMMON). IERSL = 0 IF NO TROUBLE OCCURRED.
3269 C IERSL = 1 IF A SINGULAR MATRIX AROSE WITH MITER = 3.
3270 C THIS ROUTINE ALSO USES THE COMMON VARIABLES EL0, H, MITER, AND N.
3271 C-----------------------------------------------------------------------
3272 IERSL = 0
3273 GO TO (100, 100, 300, 400, 400), MITER
3274 100 CALL DGESL (WM(3), N, N, IWM(21), X, 0)
3275 RETURN
3277 300 PHL0 = WM(2)
3278 HL0 = H*EL0
3279 WM(2) = HL0
3280 IF (HL0 .EQ. PHL0) GO TO 330
3281 R = HL0/PHL0
3282 DO 320 I = 1,N
3283 DI = ONE - R*(ONE - ONE/WM(I+2))
3284 IF (ABS(DI) .EQ. ZERO) GO TO 390
3285 320 WM(I+2) = ONE/DI
3286 330 DO 340 I = 1,N
3287 340 X(I) = WM(I+2)*X(I)
3288 RETURN
3289 390 IERSL = 1
3290 RETURN
3292 400 ML = IWM(1)
3293 MU = IWM(2)
3294 MEBAND = 2*ML + MU + 1
3295 CALL DGBSL (WM(3), MEBAND, N, ML, MU, IWM(21), X, 0)
3296 RETURN
3297 C----------------------- END OF SUBROUTINE SOLSY -----------------------
3299 SUBROUTINE EWSET (N, ITOL, RTOL, ATOL, YCUR, EWT)
3300 C-----------------------------------------------------------------------
3301 C THIS SUBROUTINE SETS THE ERROR WEIGHT VECTOR EWT ACCORDING TO
3302 C EWT(I) = RTOL(I)*ABS(YCUR(I)) + ATOL(I), I = 1,...,N,
3303 C WITH THE SUBSCRIPT ON RTOL AND/OR ATOL POSSIBLY REPLACED BY 1 ABOVE,
3304 C DEPENDING ON THE VALUE OF ITOL.
3305 C-----------------------------------------------------------------------
3306 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3307 DIMENSION RTOL(*), ATOL(*), YCUR(N), EWT(N)
3308 RTOLI = RTOL(1)
3309 ATOLI = ATOL(1)
3310 DO 10 I = 1,N
3311 IF (ITOL .GE. 3) RTOLI = RTOL(I)
3312 IF (ITOL .EQ. 2 .OR. ITOL .EQ. 4) ATOLI = ATOL(I)
3313 EWT(I) = RTOLI*ABS(YCUR(I)) + ATOLI
3314 10 CONTINUE
3315 RETURN
3316 C----------------------- END OF SUBROUTINE EWSET -----------------------
3318 DOUBLE PRECISION FUNCTION VNORM (N, V, W)
3319 C-----------------------------------------------------------------------
3320 C THIS FUNCTION ROUTINE COMPUTES THE WEIGHTED ROOT-MEAN-SQUARE NORM
3321 C OF THE VECTOR OF LENGTH N CONTAINED IN THE ARRAY V, WITH WEIGHTS
3322 C CONTAINED IN THE ARRAY W OF LENGTH N..
3323 C VNORM = SQRT( (1/N) * SUM( V(I)*W(I) )**2 )
3324 C PROTECTION FOR UNDERFLOW/OVERFLOW IS ACCOMPLISHED USING TWO
3325 C CONSTANTS WHICH ARE HOPEFULLY APPLICABLE FOR ALL MACHINES.
3326 C THESE ARE:
3327 C CUTLO = maximum of SQRT(U/EPS) over all known machines
3328 C CUTHI = minimum of SQRT(Z) over all known machines
3329 C WHERE
3330 C EPS = smallest number s.t. EPS + 1 .GT. 1
3331 C U = smallest positive number (underflow limit)
3332 C Z = largest number (overflow limit)
3334 C DETAILS OF THE ALGORITHM AND OF VALUES OF CUTLO AND CUTHI ARE
3335 C FOUND IN THE BLAS ROUTINE SNRM2 (SEE ALSO ALGORITHM 539, TRANS.
3336 C MATH. SOFTWARE, VOL. 5 NO. 3, 1979, 308-323.
3337 C FOR SINGLE PRECISION, THE FOLLOWING VALUES SHOULD BE UNIVERSAL:
3338 C DATA CUTLO,CUTHI /4.441E-16,1.304E19/
3339 C FOR DOUBLE PRECISION, USE
3340 C DATA CUTLO,CUTHI /8.232D-11,1.304D19/
3342 C-----------------------------------------------------------------------
3343 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3344 INTEGER NEXT,I,J,N
3345 DIMENSION V(N),W(N)
3346 DATA CUTLO,CUTHI /8.232D-11,1.304D19/
3347 DATA ZERO,ONE/0.0D0,1.0D0/
3348 C BLAS ALGORITHM
3349 NEXT = 1
3350 SUM = ZERO
3351 I = 1
3352 20 SX = V(I)*W(I)
3353 GO TO (30,40,70,80),NEXT
3354 30 IF (ABS(SX).GT.CUTLO) GO TO 110
3355 NEXT = 2
3356 XMAX = ZERO
3357 40 IF (SX.EQ.ZERO) GO TO 130
3358 IF (ABS(SX).GT.CUTLO) GO TO 110
3359 NEXT = 3
3360 GO TO 60
3361 50 I=J
3362 NEXT = 4
3363 SUM = (SUM/SX)/SX
3364 60 XMAX = ABS(SX)
3365 GO TO 90
3366 70 IF(ABS(SX).GT.CUTLO) GO TO 100
3367 80 IF(ABS(SX).LE.XMAX) GO TO 90
3368 SUM = ONE + SUM * (XMAX/SX)**2
3369 XMAX = ABS(SX)
3370 GO TO 130
3371 90 SUM = SUM + (SX/XMAX)**2
3372 GO TO 130
3373 100 SUM = (SUM*XMAX)*XMAX
3374 110 HITEST = CUTHI/REAL(N)
3375 DO 120 J = I,N
3376 SX = V(J)*W(J)
3377 IF(ABS(SX).GE.HITEST) GO TO 50
3378 SUM = SUM + SX**2
3379 120 CONTINUE
3380 VNORM = SQRT(SUM)
3381 GO TO 140
3382 130 CONTINUE
3383 I = I + 1
3384 IF (I.LE.N) GO TO 20
3385 VNORM = XMAX * SQRT(SUM)
3386 140 CONTINUE
3387 RETURN
3388 C----------------------- END OF FUNCTION VNORM -------------------------
3390 SUBROUTINE SVCOM (RSAV, ISAV)
3391 C-----------------------------------------------------------------------
3392 C THIS ROUTINE STORES IN RSAV AND ISAV THE CONTENTS OF COMMON BLOCKS
3393 C ODE001, ODE002 AND EH0001, WHICH ARE USED INTERNALLY IN THE ODESSA
3394 C PACKAGE.
3395 C RSAV = REAL ARRAY OF LENGTH 222 OR MORE.
3396 C ISAV = INTEGER ARRAY OF LENGTH 52 OR MORE.
3397 C-----------------------------------------------------------------------
3398 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3399 DIMENSION RSAV(*), ISAV(*)
3400 COMMON /ODE001/ RODE1(219), IODE1(39)
3401 COMMON /ODE002/ RODE2(3), IODE2(11)
3402 COMMON /EH0001/ IEH(2)
3403 DATA LRODE1/219/, LIODE1/39/, LRODE2/3/, LIODE2/11/
3405 DO 10 I = 1,LRODE1
3406 10 RSAV(I) = RODE1(I)
3407 DO 20 I = 1,LRODE2
3408 J = LRODE1 + I
3409 20 RSAV(J) = RODE2(I)
3410 DO 30 I = 1,LIODE1
3411 30 ISAV(I) = IODE1(I)
3412 DO 40 I = 1,LIODE2
3413 J = LIODE1 + I
3414 40 ISAV(J) = IODE2(I)
3415 ISAV(LIODE1+LIODE2+1) = IEH(1)
3416 ISAV(LIODE1+LIODE2+2) = IEH(2)
3417 RETURN
3418 C----------------------- END OF SUBROUTINE SVCOM -----------------------
3420 SUBROUTINE RSCOM (RSAV, ISAV)
3421 C-----------------------------------------------------------------------
3422 C THIS ROUTINE RESTORES FROM RSAV AND ISAV THE CONTENTS OF COMMON BLOCKS
3423 C ODE001, ODE002 AND EH0001, WHICH ARE USED INTERNALLY IN THE ODESSSA
3424 C PACKAGE. THIS PRESUMES THAT RSAV AND ISAV WERE LOADED BY MEANS
3425 C OF SUBROUTINE SVCOM OR THE EQUIVALENT.
3426 C-----------------------------------------------------------------------
3427 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3428 DIMENSION RSAV(*), ISAV(*)
3429 COMMON /ODE001/ RODE1(219), IODE1(39)
3430 COMMON /ODE002/ RODE2(3), IODE2(11)
3431 COMMON /EH0001/ IEH(2)
3432 DATA LRODE1/219/, LIODE1/39/, LRODE2/3/, LIODE2/11/
3434 DO 10 I = 1,LRODE1
3435 10 RODE1(I) = RSAV(I)
3436 DO 20 I = 1,LRODE2
3437 J = LRODE1 + I
3438 20 RODE2(I) = RSAV(J)
3439 DO 30 I = 1,LIODE1
3440 30 IODE1(I) = ISAV(I)
3441 DO 40 I = 1,LODE2
3442 J = LIODE1 + I
3443 40 IODE2(I) = ISAV(J)
3444 IEH(1) = ISAV(LIODE1+LIODE2+1)
3445 IEH(2) = ISAV(LIODE1+LIODE2+2)
3446 RETURN
3447 C----------------------- END OF SUBROUTINE RSCOM -----------------------
3449 SUBROUTINE DGEFA(A,LDA,N,IPVT,INFO)
3450 INTEGER LDA,N,IPVT(*),INFO
3451 DOUBLE PRECISION A(LDA,*)
3453 C DGEFA FACTORS A DOUBLE PRECISION MATRIX BY GAUSSIAN ELIMINATION.
3455 C DGEFA IS USUALLY CALLED BY DGECO, BUT IT CAN BE CALLED
3456 C DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED.
3457 C (TIME FOR DGECO) = (1 + 9/N)*(TIME FOR DGEFA) .
3459 C ON ENTRY
3461 C A DOUBLE PRECISION(LDA, N)
3462 C THE MATRIX TO BE FACTORED.
3464 C LDA INTEGER
3465 C THE LEADING DIMENSION OF THE ARRAY A .
3467 C N INTEGER
3468 C THE ORDER OF THE MATRIX A .
3470 C ON RETURN
3472 C A AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS
3473 C WHICH WERE USED TO OBTAIN IT.
3474 C THE FACTORIZATION CAN BE WRITTEN A = L*U WHERE
3475 C L IS A PRODUCT OF PERMUTATION AND UNIT LOWER
3476 C TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR.
3478 C IPVT INTEGER(N)
3479 C AN INTEGER VECTOR OF PIVOT INDICES.
3481 C INFO INTEGER
3482 C = 0 NORMAL VALUE.
3483 C = K IF U(K,K) .EQ. 0.0 . THIS IS NOT AN ERROR
3484 C CONDITION FOR THIS SUBROUTINE, BUT IT DOES
3485 C INDICATE THAT DGESL OR DGEDI WILL DIVIDE BY ZERO
3486 C IF CALLED. USE RCOND IN DGECO FOR A RELIABLE
3487 C INDICATION OF SINGULARITY.
3489 C LINPACK. THIS VERSION DATED 08/14/78 .
3490 C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
3492 C SUBROUTINES AND FUNCTIONS
3494 C BLAS DAXPY,DSCAL,IDAMAX
3496 C INTERNAL VARIABLES
3498 DOUBLE PRECISION T
3499 INTEGER IDAMAX,J,K,KP1,L,NM1
3502 C GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING
3504 INFO = 0
3505 NM1 = N - 1
3506 IF (NM1 .LT. 1) GO TO 70
3507 DO 60 K = 1, NM1
3508 KP1 = K + 1
3510 C FIND L = PIVOT INDEX
3512 L = IDAMAX(N-K+1,A(K,K),1) + K - 1
3513 IPVT(K) = L
3515 C ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED
3517 IF (A(L,K) .EQ. 0.0D0) GO TO 40
3519 C INTERCHANGE IF NECESSARY
3521 IF (L .EQ. K) GO TO 10
3522 T = A(L,K)
3523 A(L,K) = A(K,K)
3524 A(K,K) = T
3525 10 CONTINUE
3527 C COMPUTE MULTIPLIERS
3529 T = -1.0D0/A(K,K)
3530 CALL DSCAL(N-K,T,A(K+1,K),1)
3532 C ROW ELIMINATION WITH COLUMN INDEXING
3534 DO 30 J = KP1, N
3535 T = A(L,J)
3536 IF (L .EQ. K) GO TO 20
3537 A(L,J) = A(K,J)
3538 A(K,J) = T
3539 20 CONTINUE
3540 CALL DAXPY(N-K,T,A(K+1,K),1,A(K+1,J),1)
3541 30 CONTINUE
3542 GO TO 50
3543 40 CONTINUE
3544 INFO = K
3545 50 CONTINUE
3546 60 CONTINUE
3547 70 CONTINUE
3548 IPVT(N) = N
3549 IF (A(N,N) .EQ. 0.0D0) INFO = N
3550 RETURN
3552 SUBROUTINE DGESL(A,LDA,N,IPVT,B,JOB)
3553 INTEGER LDA,N,IPVT(*),JOB
3554 DOUBLE PRECISION A(LDA,*),B(*)
3556 C DGESL KppSolveS THE DOUBLE PRECISION SYSTEM
3557 C A * X = B OR TRANS(A) * X = B
3558 C USING THE FACTORS COMPUTED BY DGECO OR DGEFA.
3560 C ON ENTRY
3562 C A DOUBLE PRECISION(LDA, N)
3563 C THE OUTPUT FROM DGECO OR DGEFA.
3565 C LDA INTEGER
3566 C THE LEADING DIMENSION OF THE ARRAY A .
3568 C N INTEGER
3569 C THE ORDER OF THE MATRIX A .
3571 C IPVT INTEGER(N)
3572 C THE PIVOT VECTOR FROM DGECO OR DGEFA.
3574 C B DOUBLE PRECISION(N)
3575 C THE RIGHT HAND SIDE VECTOR.
3577 C JOB INTEGER
3578 C = 0 TO KppSolve A*X = B ,
3579 C = NONZERO TO KppSolve TRANS(A)*X = B WHERE
3580 C TRANS(A) IS THE TRANSPOSE.
3582 C ON RETURN
3584 C B THE SOLUTION VECTOR X .
3586 C ERROR CONDITION
3588 C A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS A
3589 C ZERO ON THE DIAGONAL. TECHNICALLY THIS INDICATES SINGULARITY
3590 C BUT IT IS OFTEN CAUSED BY IMPROPER ARGUMENTS OR IMPROPER
3591 C SETTING OF LDA . IT WILL NOT OCCUR IF THE SUBROUTINES ARE
3592 C CALLED CORRECTLY AND IF DGECO HAS SET RCOND .GT. 0.0
3593 C OR DGEFA HAS SET INFO .EQ. 0 .
3595 C TO COMPUTE INVERSE(A) * C WHERE C IS A MATRIX
3596 C WITH P COLUMNS
3597 C CALL DGECO(A,LDA,N,IPVT,RCOND,Z)
3598 C IF (RCOND IS TOO SMALL) GO TO ...
3599 C DO 10 J = 1, P
3600 C CALL DGESL(A,LDA,N,IPVT,C(1,J),0)
3601 C 10 CONTINUE
3603 C LINPACK. THIS VERSION DATED 08/14/78 .
3604 C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
3606 C SUBROUTINES AND FUNCTIONS
3608 C BLAS DAXPY,DDOT
3610 C INTERNAL VARIABLES
3612 DOUBLE PRECISION DDOT,T
3613 INTEGER K,KB,L,NM1
3615 NM1 = N - 1
3616 IF (JOB .NE. 0) GO TO 50
3618 C JOB = 0 , KppSolve A * X = B
3619 C FIRST KppSolve L*Y = B
3621 IF (NM1 .LT. 1) GO TO 30
3622 DO 20 K = 1, NM1
3623 L = IPVT(K)
3624 T = B(L)
3625 IF (L .EQ. K) GO TO 10
3626 B(L) = B(K)
3627 B(K) = T
3628 10 CONTINUE
3629 CALL DAXPY(N-K,T,A(K+1,K),1,B(K+1),1)
3630 20 CONTINUE
3631 30 CONTINUE
3633 C NOW KppSolve U*X = Y
3635 DO 40 KB = 1, N
3636 K = N + 1 - KB
3637 B(K) = B(K)/A(K,K)
3638 T = -B(K)
3639 CALL DAXPY(K-1,T,A(1,K),1,B(1),1)
3640 40 CONTINUE
3641 GO TO 100
3642 50 CONTINUE
3644 C JOB = NONZERO, KppSolve TRANS(A) * X = B
3645 C FIRST KppSolve TRANS(U)*Y = B
3647 DO 60 K = 1, N
3648 T = DDOT(K-1,A(1,K),1,B(1),1)
3649 B(K) = (B(K) - T)/A(K,K)
3650 60 CONTINUE
3652 C NOW KppSolve TRANS(L)*X = Y
3654 IF (NM1 .LT. 1) GO TO 90
3655 DO 80 KB = 1, NM1
3656 K = N - KB
3657 B(K) = B(K) + DDOT(N-K,A(K+1,K),1,B(K+1),1)
3658 L = IPVT(K)
3659 IF (L .EQ. K) GO TO 70
3660 T = B(L)
3661 B(L) = B(K)
3662 B(K) = T
3663 70 CONTINUE
3664 80 CONTINUE
3665 90 CONTINUE
3666 100 CONTINUE
3667 RETURN
3669 SUBROUTINE DGBFA(ABD,LDA,N,ML,MU,IPVT,INFO)
3670 INTEGER LDA,N,ML,MU,IPVT(*),INFO
3671 DOUBLE PRECISION ABD(LDA,*)
3673 C DGBFA FACTORS A DOUBLE PRECISION BAND MATRIX BY ELIMINATION.
3675 C DGBFA IS USUALLY CALLED BY DGBCO, BUT IT CAN BE CALLED
3676 C DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED.
3678 C ON ENTRY
3680 C ABD DOUBLE PRECISION(LDA, N)
3681 C CONTAINS THE MATRIX IN BAND STORAGE. THE COLUMNS
3682 C OF THE MATRIX ARE STORED IN THE COLUMNS OF ABD AND
3683 C THE DIAGONALS OF THE MATRIX ARE STORED IN ROWS
3684 C ML+1 THROUGH 2*ML+MU+1 OF ABD .
3685 C SEE THE COMMENTS BELOW FOR DETAILS.
3687 C LDA INTEGER
3688 C THE LEADING DIMENSION OF THE ARRAY ABD .
3689 C LDA MUST BE .GE. 2*ML + MU + 1 .
3691 C N INTEGER
3692 C THE ORDER OF THE ORIGINAL MATRIX.
3694 C ML INTEGER
3695 C NUMBER OF DIAGONALS BELOW THE MAIN DIAGONAL.
3696 C 0 .LE. ML .LT. N .
3698 C MU INTEGER
3699 C NUMBER OF DIAGONALS ABOVE THE MAIN DIAGONAL.
3700 C 0 .LE. MU .LT. N .
3701 C MORE EFFICIENT IF ML .LE. MU .
3702 C ON RETURN
3704 C ABD AN UPPER TRIANGULAR MATRIX IN BAND STORAGE AND
3705 C THE MULTIPLIERS WHICH WERE USED TO OBTAIN IT.
3706 C THE FACTORIZATION CAN BE WRITTEN A = L*U WHERE
3707 C L IS A PRODUCT OF PERMUTATION AND UNIT LOWER
3708 C TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR.
3710 C IPVT INTEGER(N)
3711 C AN INTEGER VECTOR OF PIVOT INDICES.
3713 C INFO INTEGER
3714 C = 0 NORMAL VALUE.
3715 C = K IF U(K,K) .EQ. 0.0 . THIS IS NOT AN ERROR
3716 C CONDITION FOR THIS SUBROUTINE, BUT IT DOES
3717 C INDICATE THAT DGBSL WILL DIVIDE BY ZERO IF
3718 C CALLED. USE RCOND IN DGBCO FOR A RELIABLE
3719 C INDICATION OF SINGULARITY.
3721 C BAND STORAGE
3723 C IF A IS A BAND MATRIX, THE FOLLOWING PROGRAM SEGMENT
3724 C WILL SET UP THE INPUT.
3726 C ML = (BAND WIDTH BELOW THE DIAGONAL)
3727 C MU = (BAND WIDTH ABOVE THE DIAGONAL)
3728 C M = ML + MU + 1
3729 C DO 20 J = 1, N
3730 C I1 = MAX0(1, J-MU)
3731 C I2 = MIN0(N, J+ML)
3732 C DO 10 I = I1, I2
3733 C K = I - J + M
3734 C ABD(K,J) = A(I,J)
3735 C 10 CONTINUE
3736 C 20 CONTINUE
3738 C THIS USES ROWS ML+1 THROUGH 2*ML+MU+1 OF ABD .
3739 C IN ADDITION, THE FIRST ML ROWS IN ABD ARE USED FOR
3740 C ELEMENTS GENERATED DURING THE TRIANGULARIZATION.
3741 C THE TOTAL NUMBER OF ROWS NEEDED IN ABD IS 2*ML+MU+1 .
3742 C THE ML+MU BY ML+MU UPPER LEFT TRIANGLE AND THE
3743 C ML BY ML LOWER RIGHT TRIANGLE ARE NOT REFERENCED.
3745 C LINPACK. THIS VERSION DATED 08/14/78 .
3746 C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
3748 C SUBROUTINES AND FUNCTIONS
3750 C BLAS DAXPY,DSCAL,IDAMAX
3751 C FORTRAN MAX0,MIN0
3753 C INTERNAL VARIABLES
3755 DOUBLE PRECISION T
3756 INTEGER I,IDAMAX,I0,J,JU,JZ,J0,J1,K,KP1,L,LM,M,MM,NM1
3759 M = ML + MU + 1
3760 INFO = 0
3762 C ZERO INITIAL FILL-IN COLUMNS
3764 J0 = MU + 2
3765 J1 = MIN0(N,M) - 1
3766 IF (J1 .LT. J0) GO TO 30
3767 DO 20 JZ = J0, J1
3768 I0 = M + 1 - JZ
3769 DO 10 I = I0, ML
3770 ABD(I,JZ) = 0.0D0
3771 10 CONTINUE
3772 20 CONTINUE
3773 30 CONTINUE
3774 JZ = J1
3775 JU = 0
3777 C GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING
3779 NM1 = N - 1
3780 IF (NM1 .LT. 1) GO TO 130
3781 DO 120 K = 1, NM1
3782 KP1 = K + 1
3784 C ZERO NEXT FILL-IN COLUMN
3786 JZ = JZ + 1
3787 IF (JZ .GT. N) GO TO 50
3788 IF (ML .LT. 1) GO TO 50
3789 DO 40 I = 1, ML
3790 ABD(I,JZ) = 0.0D0
3791 40 CONTINUE
3792 50 CONTINUE
3794 C FIND L = PIVOT INDEX
3796 LM = MIN0(ML,N-K)
3797 L = IDAMAX(LM+1,ABD(M,K),1) + M - 1
3798 IPVT(K) = L + K - M
3800 C ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED
3802 IF (ABD(L,K) .EQ. 0.0D0) GO TO 100
3804 C INTERCHANGE IF NECESSARY
3806 IF (L .EQ. M) GO TO 60
3807 T = ABD(L,K)
3808 ABD(L,K) = ABD(M,K)
3809 ABD(M,K) = T
3810 60 CONTINUE
3812 C COMPUTE MULTIPLIERS
3814 T = -1.0D0/ABD(M,K)
3815 CALL DSCAL(LM,T,ABD(M+1,K),1)
3817 C ROW ELIMINATION WITH COLUMN INDEXING
3819 JU = MIN0(MAX0(JU,MU+IPVT(K)),N)
3820 MM = M
3821 IF (JU .LT. KP1) GO TO 90
3822 DO 80 J = KP1, JU
3823 L = L - 1
3824 MM = MM - 1
3825 T = ABD(L,J)
3826 IF (L .EQ. MM) GO TO 70
3827 ABD(L,J) = ABD(MM,J)
3828 ABD(MM,J) = T
3829 70 CONTINUE
3830 CALL DAXPY(LM,T,ABD(M+1,K),1,ABD(MM+1,J),1)
3831 80 CONTINUE
3832 90 CONTINUE
3833 GO TO 110
3834 100 CONTINUE
3835 INFO = K
3836 110 CONTINUE
3837 120 CONTINUE
3838 130 CONTINUE
3839 IPVT(N) = N
3840 IF (ABD(M,N) .EQ. 0.0D0) INFO = N
3841 RETURN
3843 SUBROUTINE DGBSL(ABD,LDA,N,ML,MU,IPVT,B,JOB)
3844 INTEGER LDA,N,ML,MU,IPVT(*),JOB
3845 DOUBLE PRECISION ABD(LDA,*),B(*)
3847 C DGBSL KppSolveS THE DOUBLE PRECISION BAND SYSTEM
3848 C A * X = B OR TRANS(A) * X = B
3849 C USING THE FACTORS COMPUTED BY DGBCO OR DGBFA.
3851 C ON ENTRY
3853 C ABD DOUBLE PRECISION(LDA, N)
3854 C THE OUTPUT FROM DGBCO OR DGBFA.
3856 C LDA INTEGER
3857 C THE LEADING DIMENSION OF THE ARRAY ABD .
3859 C N INTEGER
3860 C THE ORDER OF THE ORIGINAL MATRIX.
3862 C ML INTEGER
3863 C NUMBER OF DIAGONALS BELOW THE MAIN DIAGONAL.
3865 C MU INTEGER
3866 C NUMBER OF DIAGONALS ABOVE THE MAIN DIAGONAL.
3868 C IPVT INTEGER(N)
3869 C THE PIVOT VECTOR FROM DGBCO OR DGBFA.
3871 C B DOUBLE PRECISION(N)
3872 C THE RIGHT HAND SIDE VECTOR.
3874 C JOB INTEGER
3875 C = 0 TO KppSolve A*X = B ,
3876 C = NONZERO TO KppSolve TRANS(A)*X = B , WHERE
3877 C TRANS(A) IS THE TRANSPOSE.
3879 C ON RETURN
3881 C B THE SOLUTION VECTOR X .
3883 C ERROR CONDITION
3885 C A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS A
3886 C ZERO ON THE DIAGONAL. TECHNICALLY THIS INDICATES SINGULARITY
3887 C BUT IT IS OFTEN CAUSED BY IMPROPER ARGUMENTS OR IMPROPER
3888 C SETTING OF LDA . IT WILL NOT OCCUR IF THE SUBROUTINES ARE
3889 C CALLED CORRECTLY AND IF DGBCO HAS SET RCOND .GT. 0.0
3890 C OR DGBFA HAS SET INFO .EQ. 0 .
3892 C TO COMPUTE INVERSE(A) * C WHERE C IS A MATRIX
3893 C WITH P COLUMNS
3894 C CALL DGBCO(ABD,LDA,N,ML,MU,IPVT,RCOND,Z)
3895 C IF (RCOND IS TOO SMALL) GO TO ...
3896 C DO 10 J = 1, P
3897 C CALL DGBSL(ABD,LDA,N,ML,MU,IPVT,C(1,J),0)
3898 C 10 CONTINUE
3900 C LINPACK. THIS VERSION DATED 08/14/78 .
3901 C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
3903 C SUBROUTINES AND FUNCTIONS
3905 C BLAS DAXPY,DDOT
3906 C FORTRAN MIN0
3908 C INTERNAL VARIABLES
3910 DOUBLE PRECISION DDOT,T
3911 INTEGER K,KB,L,LA,LB,LM,M,NM1
3913 M = MU + ML + 1
3914 NM1 = N - 1
3915 IF (JOB .NE. 0) GO TO 50
3917 C JOB = 0 , KppSolve A * X = B
3918 C FIRST KppSolve L*Y = B
3920 IF (ML .EQ. 0) GO TO 30
3921 IF (NM1 .LT. 1) GO TO 30
3922 DO 20 K = 1, NM1
3923 LM = MIN0(ML,N-K)
3924 L = IPVT(K)
3925 T = B(L)
3926 IF (L .EQ. K) GO TO 10
3927 B(L) = B(K)
3928 B(K) = T
3929 10 CONTINUE
3930 CALL DAXPY(LM,T,ABD(M+1,K),1,B(K+1),1)
3931 20 CONTINUE
3932 30 CONTINUE
3934 C NOW KppSolve U*X = Y
3936 DO 40 KB = 1, N
3937 K = N + 1 - KB
3938 B(K) = B(K)/ABD(M,K)
3939 LM = MIN0(K,M) - 1
3940 LA = M - LM
3941 LB = K - LM
3942 T = -B(K)
3943 CALL DAXPY(LM,T,ABD(LA,K),1,B(LB),1)
3944 40 CONTINUE
3945 GO TO 100
3946 50 CONTINUE
3948 C JOB = NONZERO, KppSolve TRANS(A) * X = B
3949 C FIRST KppSolve TRANS(U)*Y = B
3951 DO 60 K = 1, N
3952 LM = MIN0(K,M) - 1
3953 LA = M - LM
3954 LB = K - LM
3955 T = DDOT(LM,ABD(LA,K),1,B(LB),1)
3956 B(K) = (B(K) - T)/ABD(M,K)
3957 60 CONTINUE
3959 C NOW KppSolve TRANS(L)*X = Y
3961 IF (ML .EQ. 0) GO TO 90
3962 IF (NM1 .LT. 1) GO TO 90
3963 DO 80 KB = 1, NM1
3964 K = N - KB
3965 LM = MIN0(ML,N-K)
3966 B(K) = B(K) + DDOT(LM,ABD(M+1,K),1,B(K+1),1)
3967 L = IPVT(K)
3968 IF (L .EQ. K) GO TO 70
3969 T = B(L)
3970 B(L) = B(K)
3971 B(K) = T
3972 70 CONTINUE
3973 80 CONTINUE
3974 90 CONTINUE
3975 100 CONTINUE
3976 RETURN
3978 SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY)
3980 C CONSTANT TIMES A VECTOR PLUS A VECTOR.
3981 C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
3982 C JACK DONGARRA, LINPACK, 3/11/78.
3984 DOUBLE PRECISION DX(*),DY(*),DA
3985 INTEGER I,INCX,INCY,IX,IY,M,MP1,N
3987 IF(N.LE.0)RETURN
3988 IF (DA .EQ. 0.0D0) RETURN
3989 IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
3991 C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
3992 C NOT EQUAL TO 1
3994 IX = 1
3995 IY = 1
3996 IF(INCX.LT.0)IX = (-N+1)*INCX + 1
3997 IF(INCY.LT.0)IY = (-N+1)*INCY + 1
3998 DO 10 I = 1,N
3999 DY(IY) = DY(IY) + DA*DX(IX)
4000 IX = IX + INCX
4001 IY = IY + INCY
4002 10 CONTINUE
4003 RETURN
4005 C CODE FOR BOTH INCREMENTS EQUAL TO 1
4008 C CLEAN-UP LOOP
4010 20 M = MOD(N,4)
4011 IF( M .EQ. 0 ) GO TO 40
4012 DO 30 I = 1,M
4013 DY(I) = DY(I) + DA*DX(I)
4014 30 CONTINUE
4015 IF( N .LT. 4 ) RETURN
4016 40 MP1 = M + 1
4017 DO 50 I = MP1,N,4
4018 DY(I) = DY(I) + DA*DX(I)
4019 DY(I + 1) = DY(I + 1) + DA*DX(I + 1)
4020 DY(I + 2) = DY(I + 2) + DA*DX(I + 2)
4021 DY(I + 3) = DY(I + 3) + DA*DX(I + 3)
4022 50 CONTINUE
4023 RETURN
4025 SUBROUTINE DSCAL(N,DA,DX,INCX)
4027 C SCALES A VECTOR BY A CONSTANT.
4028 C USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE.
4029 C JACK DONGARRA, LINPACK, 3/11/78.
4031 DOUBLE PRECISION DA,DX(*)
4032 INTEGER I,INCX,M,MP1,N,NINCX
4034 IF(N.LE.0)RETURN
4035 IF(INCX.EQ.1)GO TO 20
4037 C CODE FOR INCREMENT NOT EQUAL TO 1
4040 NINCX = N*INCX
4041 DO 10 I = 1,NINCX,INCX
4042 DX(I) = DA*DX(I)
4043 10 CONTINUE
4044 RETURN
4046 C CODE FOR INCREMENT EQUAL TO 1
4049 C CLEAN-UP LOOP
4051 20 M = MOD(N,5)
4052 IF( M .EQ. 0 ) GO TO 40
4053 DO 30 I = 1,M
4054 DX(I) = DA*DX(I)
4055 30 CONTINUE
4056 IF( N .LT. 5 ) RETURN
4057 40 MP1 = M + 1
4058 DO 50 I = MP1,N,5
4059 DX(I) = DA*DX(I)
4060 DX(I + 1) = DA*DX(I + 1)
4061 DX(I + 2) = DA*DX(I + 2)
4062 DX(I + 3) = DA*DX(I + 3)
4063 DX(I + 4) = DA*DX(I + 4)
4064 50 CONTINUE
4065 RETURN
4067 DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY)
4069 C FORMS THE DOT PRODUCT OF TWO VECTORS.
4070 C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
4071 C JACK DONGARRA, LINPACK, 3/11/78.
4073 DOUBLE PRECISION DX(*),DY(*),DTEMP
4074 INTEGER I,INCX,INCY,IX,IY,M,MP1,N
4076 DDOT = 0.0D0
4077 DTEMP = 0.0D0
4078 IF(N.LE.0)RETURN
4079 IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
4081 C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
4082 C NOT EQUAL TO 1
4084 IX = 1
4085 IY = 1
4086 IF(INCX.LT.0)IX = (-N+1)*INCX + 1
4087 IF(INCY.LT.0)IY = (-N+1)*INCY + 1
4088 DO 10 I = 1,N
4089 DTEMP = DTEMP + DX(IX)*DY(IY)
4090 IX = IX + INCX
4091 IY = IY + INCY
4092 10 CONTINUE
4093 DDOT = DTEMP
4094 RETURN
4096 C CODE FOR BOTH INCREMENTS EQUAL TO 1
4099 C CLEAN-UP LOOP
4101 20 M = MOD(N,5)
4102 IF( M .EQ. 0 ) GO TO 40
4103 DO 30 I = 1,M
4104 DTEMP = DTEMP + DX(I)*DY(I)
4105 30 CONTINUE
4106 IF( N .LT. 5 ) GO TO 60
4107 40 MP1 = M + 1
4108 DO 50 I = MP1,N,5
4109 DTEMP = DTEMP + DX(I)*DY(I) + DX(I + 1)*DY(I + 1) +
4110 * DX(I + 2)*DY(I + 2) + DX(I + 3)*DY(I + 3) + DX(I + 4)*DY(I + 4)
4111 50 CONTINUE
4112 60 DDOT = DTEMP
4113 RETURN
4115 INTEGER FUNCTION IDAMAX(N,DX,INCX)
4117 C FINDS THE INDEX OF ELEMENT HAVING MAX. ABSOLUTE VALUE.
4118 C JACK DONGARRA, LINPACK, 3/11/78.
4120 DOUBLE PRECISION DX(*),DMAX
4121 INTEGER I,INCX,IX,N
4123 IDAMAX = 0
4124 IF( N .LT. 1 ) RETURN
4125 IDAMAX = 1
4126 IF(N.EQ.1)RETURN
4127 IF(INCX.EQ.1)GO TO 20
4129 C CODE FOR INCREMENT NOT EQUAL TO 1
4131 IX = 1
4132 DMAX = DABS(DX(1))
4133 IX = IX + INCX
4134 DO 10 I = 2,N
4135 IF(DABS(DX(IX)).LE.DMAX) GO TO 5
4136 IDAMAX = I
4137 DMAX = DABS(DX(IX))
4138 5 IX = IX + INCX
4139 10 CONTINUE
4140 RETURN
4142 C CODE FOR INCREMENT EQUAL TO 1
4144 20 DMAX = DABS(DX(1))
4145 DO 30 I = 2,N
4146 IF(DABS(DX(I)).LE.DMAX) GO TO 30
4147 IDAMAX = I
4148 DMAX = DABS(DX(I))
4149 30 CONTINUE
4150 RETURN
4152 DOUBLE PRECISION FUNCTION D1MACH (IDUM)
4153 INTEGER IDUM
4154 C-----------------------------------------------------------------------
4155 C THIS ROUTINE COMPUTES THE UNIT ROUNDOFF OF THE MACHINE IN DOUBLE
4156 C PRECISION. THIS IS DEFINED AS THE SMALLEST POSITIVE MACHINE NUMBER
4157 C U SUCH THAT 1.0D0 + U .NE. 1.0D0 (IN DOUBLE PRECISION).
4158 C-----------------------------------------------------------------------
4159 DOUBLE PRECISION U, COMP
4160 U = 1.0D0
4161 10 U = U*0.5D0
4162 COMP = 1.0D0 + U
4163 IF (COMP .NE. 1.0D0) GO TO 10
4164 D1MACH = U*2.0D0
4165 RETURN
4166 C----------------------- END OF FUNCTION D1MACH ------------------------
4168 SUBROUTINE XERR (MSG, NERR, IERT, NI, I1, I2, NR, R1, R2)
4169 INTEGER NERR, IERT, NI, I1, I2, NR,
4170 1 LUN, LUNIT, MESFLG
4171 DOUBLE PRECISION R1, R2
4172 CHARACTER*(*) MSG
4173 C-------------------------------------------------------------------
4175 C ALL ARGUMENTS ARE INPUT ARGUMENTS.
4177 C MSG = THE MESSAGE (CHARACTER VARIABLE)
4178 C NERR = THE ERROR NUMBER (NOT USED).
4179 C IERT = THE ERROR TYPE..
4180 C 1 MEANS RECOVERABLE (CONTROL RETURNS TO CALLER).
4181 C 2 MEANS FATAL (RUN IS ABORTED--SEE NOTE BELOW).
4182 C NI = NUMBER OF INTEGERS (0, 1, OR 2) TO BE PRINTED WITH MESSAGE.
4183 C I1,I2 = INTEGERS TO BE PRINTED, DEPENDING ON NI.
4184 C NR = NUMBER OF REALS (0, 1, OR 2) TO BE PRINTED WITH MESSAGE.
4185 C R1,R2 = REALS TO BE PRINTED, DEPENDING ON NR.
4187 C NOTES:
4188 C 1. THE DIMENSION OF MSG IS ASSUMED TO BE AT MOST 60.
4189 C (MULTI-LINE MESSAGES ARE GENERATED BY REPEATED CALLS.)
4190 C 2. IF IERT = 2, CONTROL PASSES TO THE STATEMENT STOP
4191 C TO ABORT THE RUN. THIS STATEMENT MAY BE MACHINE-DEPENDENT.
4192 C 3. R1 AND R2 ARE ASSUMED TO BE IN DOUBLE PRECISION AND ARE PRINTED
4193 C IN D21.13 FORMAT.
4194 C 4. THE COMMON BLOCK /EH0001/ BELOW IS DATA-LOADED (A MACHINE-
4195 C DEPENDENT FEATURE) WITH DEFAULT VALUES.
4196 C THIS BLOCK IS NEEDED FOR PROPER RETENTION OF PARAMETERS USED BY
4197 C THIS ROUTINE WHICH THE USER CAN RESET BY CALLING XSETF OR XSETUN.
4198 C THE VARIABLES IN THIS BLOCK ARE AS FOLLOWS..
4199 C MESFLG = PRINT CONTROL FLAG..
4200 C 1 MEANS PRINT ALL MESSAGES (THE DEFAULT).
4201 C 0 MEANS NO PRINTING.
4202 C LUNIT = LOGICAL UNIT NUMBER FOR MESSAGES.
4203 C THE DEFAULT IS 6 (MACHINE-DEPENDENT).
4204 C 5. TO CHANGE THE DEFAULT OUTPUT UNIT, CHANGE THE DATA STATEMENT
4205 C IN THE BLOCK DATA SUBPROGRAM BELOW.
4207 C FOR A DIFFERENT RUN-ABORT COMMAND, CHANGE THE STATEMENT FOLLOWING
4208 C STATEMENT 100 AT THE END.
4209 C-----------------------------------------------------------------------
4210 COMMON /EH0001/ MESFLG, LUNIT
4211 IF (MESFLG .EQ. 0) GO TO 100
4212 C GET LOGICAL UNIT NUMBER. ---------------------------------------------
4213 LUN = LUNIT
4214 C WRITE THE MESSAGE. ---------------------------------------------------
4215 WRITE (LUN, 10) MSG
4216 10 FORMAT(1X,A)
4217 C-----------------------------------------------------------------------
4218 IF (NI .EQ. 1) WRITE (LUN, 20) I1
4219 20 FORMAT(6X,'IN ABOVE MESSAGE, I1 = ',I10)
4220 IF (NI .EQ. 2) WRITE (LUN, 30) I1,I2
4221 30 FORMAT(6X,'IN ABOVE MESSAGE, I1 = ',I10,3X,'I2 = ',I10)
4222 IF (NR .EQ. 1) WRITE (LUN, 40) R1
4223 40 FORMAT(6X,'IN ABOVE MESSAGE, R1 = ',D21.13)
4224 IF (NR .EQ. 2) WRITE (LUN, 50) R1,R2
4225 50 FORMAT(6X,'IN ABOVE, R1 = ',D21.13,3X,'R2 = ',D21.13)
4226 C ABORT THE RUN IF IERT = 2. -------------------------------------------
4227 100 IF (IERT .NE. 2) RETURN
4228 STOP
4229 C----------------------- END OF SUBROUTINE XERR ----------------------
4231 SUBROUTINE XSETF (MFLAG)
4233 C THIS ROUTINE RESETS THE PRINT CONTROL FLAG MFLAG.
4235 INTEGER MFLAG, MESFLG, LUNIT
4236 COMMON /EH0001/ MESFLG, LUNIT
4238 IF (MFLAG .EQ. 0 .OR. MFLAG .EQ. 1) MESFLG = MFLAG
4239 RETURN
4240 C----------------------- END OF SUBROUTINE XSETF -----------------------
4242 SUBROUTINE XSETUN (LUN)
4244 C THIS ROUTINE RESETS THE LOGICAL UNIT NUMBER FOR MESSAGES.
4246 INTEGER LUN, MESFLG, LUNIT
4247 COMMON /EH0001/ MESFLG, LUNIT
4249 IF (LUN .GT. 0) LUNIT = LUN
4250 RETURN
4251 C----------------------- END OF SUBROUTINE XSETUN ----------------------
4253 BLOCK DATA
4254 C-----------------------------------------------------------------------
4255 C THIS DATA SUBPROGRAM LOADS VARIABLES INTO THE INTERNAL COMMON
4256 C BLOCKS USED BY ODESSA AND ITS VARIANTS. THE VARIABLES ARE
4257 C DEFINED AS FOLLOWS..
4258 C ILLIN = COUNTER FOR THE NUMBER OF CONSECUTIVE TIMES THE PACKAGE
4259 C WAS CALLED WITH ILLEGAL INPUT. THE RUN IS STOPPED WHEN
4260 C ILLIN REACHES 5.
4261 C NTREP = COUNTER FOR THE NUMBER OF CONSECUTIVE TIMES THE PACKAGE
4262 C WAS CALLED WITH ISTATE = 1 AND TOUT = T. THE RUN IS
4263 C STOPPED WHEN NTREP REACHES 5.
4264 C MESFLG = FLAG TO CONTROL PRINTING OF ERROR MESSAGES. 1 MEANS PRINT,
4265 C 0 MEANS NO PRINTING.
4266 C LUNIT = DEFAULT VALUE OF LOGICAL UNIT NUMBER FOR PRINTING OF ERROR
4267 C MESSAGES.
4268 C-----------------------------------------------------------------------
4269 INTEGER ILLIN, IDUMA, NTREP, IDUMB, IOWNS, ICOMM, MESFLG, LUNIT
4270 DOUBLE PRECISION ROWND, ROWNS, RCOMM
4271 COMMON /ODE001/ ROWND, ROWNS(173), RCOMM(45),
4272 1 ILLIN, IDUMA(10), NTREP, IDUMB(2), IOWNS(4), ICOMM(21)
4273 COMMON /EH0001/ MESFLG, LUNIT
4274 DATA ILLIN/0/, NTREP/0/
4275 DATA MESFLG/1/, LUNIT/6/
4277 C------------------------ END OF BLOCK DATA ----------------------------
4279 C-----------------------------------------------------------------------
4280 C INSTRUCTIONS FOR INSTALLING THE ODESSA PACKAGE. (see @ below.)
4282 C ODESSA is an enhanced version of the widely disseminated ODE solver
4283 C LSODE, and as such retains the same properties regarding portability.
4284 C The notes below, adapted from the installation instructions for LSODE,
4285 C are intended to facilitate the installation of the ODESSA package in
4286 C the user's software library.
4288 C 1. Both a single and a double precision version of ODESSA are
4289 C provided in this release. It is expected that most users will
4290 C utilize the double precision version, except in the case of
4291 C extended word-length computers. Most routines used by ODESSA
4292 C are named the same regardless of whether they are single or
4293 C double precision. The exceptions are the LINPAK and BLAS
4294 C routines that follow the LINPAK/BLAS naming conventions, i.e.
4295 C D--- for a double precision routine, and S--- for a single
4296 C precision routine. Thus, care should be taken if both single
4297 C and double precision versions are stored in the same library.
4299 C 2. Several routines in ODESSA have the same names as the LSODE
4300 C routines from which they were derived, although they contain
4301 C different code. These are: INTDY, STODE, PREPJ, SVCOM, and
4302 C RSCOM. If ODESSA is added to a subroutine library of which
4303 C LSODE is already a member, these routine names must be changed
4304 C in one of the two programs. Also see the note regarding BLOCK
4305 C DATA subroutines below.
4307 C 3. In many cases, ODESSA uses unaltered LSODE routines and
4308 C common library routines that may already reside on your system.
4309 C The installation of ODESSA should be done so that identical routines
4310 C are shared rather than kept as duplicate copies.
4311 C a. Normally, the user calls only subroutine ODESSA, but for optional
4312 C capabilities the user may also CALL XSETUN, XSETF, SVCOM, RSCOM,
4313 C or INTDY, as described in Part II of the Full Description in the
4314 C User Documentation (ODESSA.DOC, see below). Except for INTDY,
4315 C none of these are called from within the package.
4316 C b. Two routines, EWSET and VNORM, are optionally replaceable by the
4317 C user if the package version is unsuitable. Hence, the install-
4318 C ation of the package should be done so that the user's version
4319 C for either routine overrides the package version.
4320 C c. The function routine D1MACH is provided to compute the unit
4321 C roundoff of the machine and precision in use, in a manner com-
4322 C patible with machine parameter routines developed at Bell Lab-
4323 C oratories. If such a routine has already been installed on
4324 C your system, the version supplied here may be discarded.
4325 C d. Linear algebraic systems are solved with routines from the
4326 C LINPACK collection, in conjunction with routines from the Basic
4327 C Linear Algebra module collection (BLAS). In double precision,
4328 C the names are DGEFA, DGESL, DGBFA, and DGBSL (from LINPACK), and
4329 C DAXPY, DSCAL, IDAMAX, and DDOT (from BLAS). If these routines
4330 C have already been installed on your system, copies supplied with
4331 C ODESSA may be discarded. The single precision versions of these
4332 C routines are used in the single precision version.
4334 C 4. There are four integer variables, in the two labeled COMMON
4335 C blocks /ODE001/ and /EH0001/, which need to be loaded with DATA
4336 C statements. They can vary during execution, and are in common to
4337 C assure their retention between calls. This is legal in ANSI Fortran
4338 C only if done in a BLOCK DATA subprogram, and this package has a
4339 C BLOCK DATA for this purpose. However, BLOCK DATA subprograms can be
4340 C difficult to install in libraries, and many compilers allow such DATA
4341 C statements in subroutines. If your system allows this, the location
4342 C of the DATA statements are just after the initial type and common
4343 C declarations in subroutines ODESSA and XERR. In ODESSA, ILLIN and
4344 C NTREP are DATA-loaded as 0. In XERR, MESFLG is loaded as 1 and
4345 C LUNIT is loaded as the appropriate default logical unit number.
4347 C 5. The ODESSA package contains subscript expressions which may not
4348 C be accepted by some compilers. Subscripts of the form I + J, I - J,
4349 C etc., occur in various routines. If any of these forms are
4350 C unacceptable to your compiler, an extra line of code setting the
4351 C subscript to a dummy integer value should be added for each subscipt.
4353 C 6. User documentation is provided in a two-level structure
4354 C to accommmodate both the casual and serious user. The novice or
4355 C casual user should need to read only the Summary of Usage and the
4356 C Example Problem located at the beginning of the documentation. More
4357 C experienced users, requiring the full set of available options,
4358 C should read the Full Description which follows the Example Problem.
4360 C 7. The user documentation may need corrections in the following ways:
4361 C a. If subroutine names have been changed to avoid conflicts between
4362 C the LSODE and ODESSA packages, the corresponding name changes
4363 C should be made in the documentation.
4364 C b. In the Summary of Usage, and in the description of XSETUN under
4365 C Part II of the Full Description, the default logical unit number
4366 C should be corrected if it is not 6.
4367 C c. In the Summary of Usage, users should be instructed to execute
4368 C CALL XSETF(1) before the first CALL to ODESSA, if this is neces-
4369 C sary for proper error message handling. (see note 2(e) above.)
4370 C d. In the description of the subroutines DF and JAC in the Summary
4371 C of Usage and in Part I of the Full Description, it is stated
4372 C that dummy names may be passed if these two routines are not user
4373 C supplied. Your system may require the user to supply a dummy
4374 C subroutine instead.
4375 C e. The ODESSA package treats the arguments NEQ, RTOL, and ATOL as
4376 C arrays (possibly of length 1), while the usage documentation
4377 C states that these arguments may be either arrays or scalars.
4378 C If your system does not allow such a mismatch, then the
4379 C documentation should be changed to reflect this.
4380 C 8. A demonstration program is provided with the package for
4381 C verification.
4384 C Jorge R. Leis and Mark A. Kramer
4385 C Department of Chemical Engineering
4386 C Massachusetts Institute of Technology
4387 C Cambridge, Massachusetts 02139
4388 C U.S.A.
4390 C Current address of J.R. Leis (Jan. 1988):
4392 C Shell Development Company
4393 C Westhollow Research Center
4394 C Houston, TX
4396 C @ Adapted from 'Instructions for Installing LSODE', written by
4397 C Alan C. Hindmarsh, Mathematics & Statistics Division, L-316,
4398 C Lawrence Livermore National Laboratory, Livermore, CA. 94550