Parent Directory | Revision Log
Moving integrators to own directory, about to make them self-contained shared libraries.
1 | subroutine lsode (f, neq, y, t, tout, itol, rtol, atol, itask, |
2 | 1 istate, iopt, rwork, lrw, iwork, liw, jac, mf) |
3 | external f, jac, xascwv |
4 | integer neq, itol, itask, istate, iopt, lrw, iwork, liw, mf |
5 | double precision y, t, tout, rtol, atol, rwork |
6 | dimension neq(1), y(1), rtol(1), atol(1), rwork(lrw), iwork(liw) |
7 | c----------------------------------------------------------------------- |
8 | c this is the march 30, 1987 version of |
9 | c lsode.. livermore solver for ordinary differential equations. |
10 | c this version is in double precision. |
11 | c |
12 | c lsode solves the initial value problem for stiff or nonstiff |
13 | c systems of first order ode-s, |
14 | c dy/dt = f(t,y) , or, in component form, |
15 | c dy(i)/dt = f(i) = f(i,t,y(1),y(2),...,y(neq)) (i = 1,...,neq). |
16 | c lsode is a package based on the gear and gearb packages, and on the |
17 | c october 23, 1978 version of the tentative odepack user interface |
18 | c standard, with minor modifications. |
19 | c----------------------------------------------------------------------- |
20 | c reference.. |
21 | c alan c. hindmarsh, odepack, a systematized collection of ode |
22 | c solvers, in scientific computing, r. s. stepleman et al. (eds.), |
23 | c north-holland, amsterdam, 1983, pp. 55-64. |
24 | c----------------------------------------------------------------------- |
25 | c author and contact.. alan c. hindmarsh, |
26 | c computing and mathematics research div., l-316 |
27 | c lawrence livermore national laboratory |
28 | c livermore, ca 94550. |
29 | c----------------------------------------------------------------------- |
30 | c summary of usage. |
31 | c |
32 | c communication between the user and the lsode package, for normal |
33 | c situations, is summarized here. this summary describes only a subset |
34 | c of the full set of options available. see the full description for |
35 | c details, including optional communication, nonstandard options, |
36 | c and instructions for special situations. see also the example |
37 | c problem (with program and output) following this summary. |
38 | c |
39 | c a. first provide a subroutine of the form.. |
40 | c subroutine f (neq, t, y, ydot) |
41 | c dimension y(neq), ydot(neq) |
42 | c which supplies the vector function f by loading ydot(i) with f(i). |
43 | c |
44 | c b. next determine (or guess) whether or not the problem is stiff. |
45 | c stiffness occurs when the jacobian matrix df/dy has an eigenvalue |
46 | c whose real part is negative and large in magnitude, compared to the |
47 | c reciprocal of the t span of interest. if the problem is nonstiff, |
48 | c use a method flag mf = 10. if it is stiff, there are four standard |
49 | c choices for mf, and lsode requires the jacobian matrix in some form. |
50 | c this matrix is regarded either as full (mf = 21 or 22), |
51 | c or banded (mf = 24 or 25). in the banded case, lsode requires two |
52 | c half-bandwidth parameters ml and mu. these are, respectively, the |
53 | c widths of the lower and upper parts of the band, excluding the main |
54 | c diagonal. thus the band consists of the locations (i,j) with |
55 | c i-ml .le. j .le. i+mu, and the full bandwidth is ml+mu+1. |
56 | c |
57 | c c. if the problem is stiff, you are encouraged to supply the jacobian |
58 | c directly (mf = 21 or 24), but if this is not feasible, lsode will |
59 | c compute it internally by difference quotients (mf = 22 or 25). |
60 | c if you are supplying the jacobian, provide a subroutine of the form.. |
61 | c subroutine jac (neq, t, y, ml, mu, pd, nrowpd) |
62 | c dimension y(neq), pd(nrowpd,neq) |
63 | c which supplies df/dy by loading pd as follows.. |
64 | c for a full jacobian (mf = 21), load pd(i,j) with df(i)/dy(j), |
65 | c the partial derivative of f(i) with respect to y(j). (ignore the |
66 | c ml and mu arguments in this case.) |
67 | c for a banded jacobian (mf = 24), load pd(i-j+mu+1,j) with |
68 | c df(i)/dy(j), i.e. load the diagonal lines of df/dy into the rows of |
69 | c pd from the top down. |
70 | c in either case, only nonzero elements need be loaded. |
71 | c |
72 | c d. write a main program which calls subroutine lsode once for |
73 | c each point at which answers are desired. this should also provide |
74 | c for possible use of logical unit 6 for output of error messages |
75 | c by lsode. on the first call to lsode, supply arguments as follows.. |
76 | c f = name of subroutine for right-hand side vector f. |
77 | c this name must be declared external in calling program. |
78 | c neq = number of first order ode-s. |
79 | c y = array of initial values, of length neq. |
80 | c t = the initial value of the independent variable. |
81 | c tout = first point where output is desired (.ne. t). |
82 | c itol = 1 or 2 according as atol (below) is a scalar or array. |
83 | c rtol = relative tolerance parameter (scalar). |
84 | c atol = absolute tolerance parameter (scalar or array). |
85 | c the estimated local error in y(i) will be controlled so as |
86 | c to be roughly less (in magnitude) than |
87 | c ewt(i) = rtol*abs(y(i)) + atol if itol = 1, or |
88 | c ewt(i) = rtol*abs(y(i)) + atol(i) if itol = 2. |
89 | c thus the local error test passes if, in each component, |
90 | c either the absolute error is less than atol (or atol(i)), |
91 | c or the relative error is less than rtol. |
92 | c use rtol = 0.0 for pure absolute error control, and |
93 | c use atol = 0.0 (or atol(i) = 0.0) for pure relative error |
94 | c control. caution.. actual (global) errors may exceed these |
95 | c local tolerances, so choose them conservatively. |
96 | c itask = 1 for normal computation of output values of y at t = tout. |
97 | c istate = integer flag (input and output). set istate = 1. |
98 | c iopt = 0 to indicate no optional inputs used. |
99 | c rwork = real work array of length at least.. |
100 | c 20 + 16*neq for mf = 10, |
101 | c 22 + 9*neq + neq**2 for mf = 21 or 22, |
102 | c 22 + 10*neq + (2*ml + mu)*neq for mf = 24 or 25. |
103 | c lrw = declared length of rwork (in user-s dimension). |
104 | c iwork = integer work array of length at least.. |
105 | c 20 for mf = 10, |
106 | c 20 + neq for mf = 21, 22, 24, or 25. |
107 | c if mf = 24 or 25, input in iwork(1),iwork(2) the lower |
108 | c and upper half-bandwidths ml,mu. |
109 | c liw = declared length of iwork (in user-s dimension). |
110 | c jac = name of subroutine for jacobian matrix (mf = 21 or 24). |
111 | c if used, this name must be declared external in calling |
112 | c program. if not used, pass a dummy name. |
113 | c mf = method flag. standard values are.. |
114 | c 10 for nonstiff (adams) method, no jacobian used. |
115 | c 21 for stiff (bdf) method, user-supplied full jacobian. |
116 | c 22 for stiff method, internally generated full jacobian. |
117 | c 24 for stiff method, user-supplied banded jacobian. |
118 | c 25 for stiff method, internally generated banded jacobian. |
119 | c note that the main program must declare arrays y, rwork, iwork, |
120 | c and possibly atol. |
121 | c |
122 | c e. the output from the first call (or any call) is.. |
123 | c y = array of computed values of y(t) vector. |
124 | c t = corresponding value of independent variable (normally tout). |
125 | c istate = 2 if lsode was successful, negative otherwise. |
126 | c -1 means excess work done on this call (perhaps wrong mf). |
127 | c -2 means excess accuracy requested (tolerances too small). |
128 | c -3 means illegal input detected (see printed message). |
129 | c -4 means repeated error test failures (check all inputs). |
130 | c -5 means repeated convergence failures (perhaps bad jacobian |
131 | c supplied or wrong choice of mf or tolerances). |
132 | c -6 means error weight became zero during problem. (solution |
133 | c component i vanished, and atol or atol(i) = 0.) |
134 | c |
135 | c f. to continue the integration after a successful return, simply |
136 | c reset tout and call lsode again. no other parameters need be reset. |
137 | c |
138 | c----------------------------------------------------------------------- |
139 | c example problem. |
140 | c |
141 | c the following is a simple example problem, with the coding |
142 | c needed for its solution by lsode. the problem is from chemical |
143 | c kinetics, and consists of the following three rate equations.. |
144 | c dy1/dt = -.04*y1 + 1.e4*y2*y3 |
145 | c dy2/dt = .04*y1 - 1.e4*y2*y3 - 3.e7*y2**2 |
146 | c dy3/dt = 3.e7*y2**2 |
147 | c on the interval from t = 0.0 to t = 4.e10, with initial conditions |
148 | c y1 = 1.0, y2 = y3 = 0. the problem is stiff. |
149 | c |
150 | c the following coding solves this problem with lsode, using mf = 21 |
151 | c and printing results at t = .4, 4., ..., 4.e10. it uses |
152 | c itol = 2 and atol much smaller for y2 than y1 or y3 because |
153 | c y2 has much smaller values. |
154 | c at the end of the run, statistical quantities of interest are |
155 | c printed (see optional outputs in the full description below). |
156 | c |
157 | c external fex, jex |
158 | c double precision atol, rtol, rwork, t, tout, y |
159 | c dimension y(3), atol(3), rwork(58), iwork(23) |
160 | c neq = 3 |
161 | c y(1) = 1.d0 |
162 | c y(2) = 0.d0 |
163 | c y(3) = 0.d0 |
164 | c t = 0.d0 |
165 | c tout = .4d0 |
166 | c itol = 2 |
167 | c rtol = 1.d-4 |
168 | c atol(1) = 1.d-6 |
169 | c atol(2) = 1.d-10 |
170 | c atol(3) = 1.d-6 |
171 | c itask = 1 |
172 | c istate = 1 |
173 | c iopt = 0 |
174 | c lrw = 58 |
175 | c liw = 23 |
176 | c mf = 21 |
177 | c do 40 iout = 1,12 |
178 | c call lsode(fex,neq,y,t,tout,itol,rtol,atol,itask,istate, |
179 | c 1 iopt,rwork,lrw,iwork,liw,jex,mf) |
180 | c write(6,20)t,y(1),y(2),y(3) |
181 | c 20 format(7h at t =,e12.4,6h y =,3e14.6) |
182 | c if (istate .lt. 0) go to 80 |
183 | c 40 tout = tout*10.d0 |
184 | c write(6,60)iwork(11),iwork(12),iwork(13) |
185 | c 60 format(/12h no. steps =,i4,11h no. f-s =,i4,11h no. j-s =,i4) |
186 | c stop |
187 | c 80 write(6,90)istate |
188 | c 90 format(///22h error halt.. istate =,i3) |
189 | c stop |
190 | c end |
191 | c |
192 | c subroutine fex (neq, t, y, ydot) |
193 | c double precision t, y, ydot |
194 | c dimension y(3), ydot(3) |
195 | c ydot(1) = -.04d0*y(1) + 1.d4*y(2)*y(3) |
196 | c ydot(3) = 3.d7*y(2)*y(2) |
197 | c ydot(2) = -ydot(1) - ydot(3) |
198 | c return |
199 | c end |
200 | c |
201 | c subroutine jex (neq, t, y, ml, mu, pd, nrpd) |
202 | c double precision pd, t, y |
203 | c dimension y(3), pd(nrpd,3) |
204 | c pd(1,1) = -.04d0 |
205 | c pd(1,2) = 1.d4*y(3) |
206 | c pd(1,3) = 1.d4*y(2) |
207 | c pd(2,1) = .04d0 |
208 | c pd(2,3) = -pd(1,3) |
209 | c pd(3,2) = 6.d7*y(2) |
210 | c pd(2,2) = -pd(1,2) - pd(3,2) |
211 | c return |
212 | c end |
213 | c |
214 | c the output of this program (on a cdc-7600 in single precision) |
215 | c is as follows.. |
216 | c |
217 | c at t = 4.0000e-01 y = 9.851726e-01 3.386406e-05 1.479357e-02 |
218 | c at t = 4.0000e+00 y = 9.055142e-01 2.240418e-05 9.446344e-02 |
219 | c at t = 4.0000e+01 y = 7.158050e-01 9.184616e-06 2.841858e-01 |
220 | c at t = 4.0000e+02 y = 4.504846e-01 3.222434e-06 5.495122e-01 |
221 | c at t = 4.0000e+03 y = 1.831701e-01 8.940379e-07 8.168290e-01 |
222 | c at t = 4.0000e+04 y = 3.897016e-02 1.621193e-07 9.610297e-01 |
223 | c at t = 4.0000e+05 y = 4.935213e-03 1.983756e-08 9.950648e-01 |
224 | c at t = 4.0000e+06 y = 5.159269e-04 2.064759e-09 9.994841e-01 |
225 | c at t = 4.0000e+07 y = 5.306413e-05 2.122677e-10 9.999469e-01 |
226 | c at t = 4.0000e+08 y = 5.494529e-06 2.197824e-11 9.999945e-01 |
227 | c at t = 4.0000e+09 y = 5.129458e-07 2.051784e-12 9.999995e-01 |
228 | c at t = 4.0000e+10 y = -7.170586e-08 -2.868234e-13 1.000000e+00 |
229 | c |
230 | c no. steps = 330 no. f-s = 405 no. j-s = 69 |
231 | c----------------------------------------------------------------------- |
232 | c full description of user interface to lsode. |
233 | c |
234 | c the user interface to lsode consists of the following parts. |
235 | c |
236 | c i. the call sequence to subroutine lsode, which is a driver |
237 | c routine for the solver. this includes descriptions of both |
238 | c the call sequence arguments and of user-supplied routines. |
239 | c following these descriptions is a description of |
240 | c optional inputs available through the call sequence, and then |
241 | c a description of optional outputs (in the work arrays). |
242 | c |
243 | c ii. descriptions of other routines in the lsode package that may be |
244 | c (optionally) called by the user. these provide the ability to |
245 | c alter error message handling, save and restore the internal |
246 | c common, and obtain specified derivatives of the solution y(t). |
247 | c |
248 | c iii. descriptions of common blocks to be declared in overlay |
249 | c or similar environments, or to be saved when doing an interrupt |
250 | c of the problem and continued solution later. |
251 | c |
252 | c iv. description of two routines in the lsode package, either of |
253 | c which the user may replace with his own version, if desired. |
254 | c these relate to the measurement of errors. |
255 | c |
256 | c----------------------------------------------------------------------- |
257 | c part i. call sequence. |
258 | c |
259 | c the call sequence parameters used for input only are |
260 | c f, neq, tout, itol, rtol, atol, itask, iopt, lrw, liw, jac, mf, |
261 | c and those used for both input and output are |
262 | c y, t, istate. |
263 | c the work arrays rwork and iwork are also used for conditional and |
264 | c optional inputs and optional outputs. (the term output here refers |
265 | c to the return from subroutine lsode to the user-s calling program.) |
266 | c |
267 | c the legality of input parameters will be thoroughly checked on the |
268 | c initial call for the problem, but not checked thereafter unless a |
269 | c change in input parameters is flagged by istate = 3 on input. |
270 | c |
271 | c the descriptions of the call arguments are as follows. |
272 | c |
273 | c f = the name of the user-supplied subroutine defining the |
274 | c ode system. the system must be put in the first-order |
275 | c form dy/dt = f(t,y), where f is a vector-valued function |
276 | c of the scalar t and the vector y. subroutine f is to |
277 | c compute the function f. it is to have the form |
278 | c subroutine f (neq, t, y, ydot) |
279 | c dimension y(1), ydot(1) |
280 | c where neq, t, and y are input, and the array ydot = f(t,y) |
281 | c is output. y and ydot are arrays of length neq. |
282 | c (in the dimension statement above, 1 is a dummy |
283 | c dimension.. it can be replaced by any value.) |
284 | c subroutine f should not alter y(1),...,y(neq). |
285 | c f must be declared external in the calling program. |
286 | c |
287 | c subroutine f may access user-defined quantities in |
288 | c neq(2),... and/or in y(neq(1)+1),... if neq is an array |
289 | c (dimensioned in f) and/or y has length exceeding neq(1). |
290 | c see the descriptions of neq and y below. |
291 | c |
292 | c if quantities computed in the f routine are needed |
293 | c externally to lsode, an extra call to f should be made |
294 | c for this purpose, for consistent and accurate results. |
295 | c if only the derivative dy/dt is needed, use intdy instead. |
296 | c |
297 | c neq = the size of the ode system (number of first order |
298 | c ordinary differential equations). used only for input. |
299 | c neq may be decreased, but not increased, during the problem. |
300 | c if neq is decreased (with istate = 3 on input), the |
301 | c remaining components of y should be left undisturbed, if |
302 | c these are to be accessed in f and/or jac. |
303 | c |
304 | c normally, neq is a scalar, and it is generally referred to |
305 | c as a scalar in this user interface description. however, |
306 | c neq may be an array, with neq(1) set to the system size. |
307 | c (the lsode package accesses only neq(1).) in either case, |
308 | c this parameter is passed as the neq argument in all calls |
309 | c to f and jac. hence, if it is an array, locations |
310 | c neq(2),... may be used to store other integer data and pass |
311 | c it to f and/or jac. subroutines f and/or jac must include |
312 | c neq in a dimension statement in that case. |
313 | c |
314 | c y = a real array for the vector of dependent variables, of |
315 | c length neq or more. used for both input and output on the |
316 | c first call (istate = 1), and only for output on other calls. |
317 | c on the first call, y must contain the vector of initial |
318 | c values. on output, y contains the computed solution vector, |
319 | c evaluated at t. if desired, the y array may be used |
320 | c for other purposes between calls to the solver. |
321 | c |
322 | c this array is passed as the y argument in all calls to |
323 | c f and jac. hence its length may exceed neq, and locations |
324 | c y(neq+1),... may be used to store other real data and |
325 | c pass it to f and/or jac. (the lsode package accesses only |
326 | c y(1),...,y(neq).) |
327 | c |
328 | c t = the independent variable. on input, t is used only on the |
329 | c first call, as the initial point of the integration. |
330 | c on output, after each call, t is the value at which a |
331 | c computed solution y is evaluated (usually the same as tout). |
332 | c on an error return, t is the farthest point reached. |
333 | c |
334 | c tout = the next value of t at which a computed solution is desired. |
335 | c used only for input. |
336 | c |
337 | c when starting the problem (istate = 1), tout may be equal |
338 | c to t for one call, then should .ne. t for the next call. |
339 | c for the initial t, an input value of tout .ne. t is used |
340 | c in order to determine the direction of the integration |
341 | c (i.e. the algebraic sign of the step sizes) and the rough |
342 | c scale of the problem. integration in either direction |
343 | c (forward or backward in t) is permitted. |
344 | c |
345 | c if itask = 2 or 5 (one-step modes), tout is ignored after |
346 | c the first call (i.e. the first call with tout .ne. t). |
347 | c otherwise, tout is required on every call. |
348 | c |
349 | c if itask = 1, 3, or 4, the values of tout need not be |
350 | c monotone, but a value of tout which backs up is limited |
351 | c to the current internal t interval, whose endpoints are |
352 | c tcur - hu and tcur (see optional outputs, below, for |
353 | c tcur and hu). |
354 | c |
355 | c itol = an indicator for the type of error control. see |
356 | c description below under atol. used only for input. |
357 | c |
358 | c rtol = a relative error tolerance parameter, either a scalar or |
359 | c an array of length neq. see description below under atol. |
360 | c input only. |
361 | c |
362 | c atol = an absolute error tolerance parameter, either a scalar or |
363 | c an array of length neq. input only. |
364 | c |
365 | c the input parameters itol, rtol, and atol determine |
366 | c the error control performed by the solver. the solver will |
367 | c control the vector e = (e(i)) of estimated local errors |
368 | c in y, according to an inequality of the form |
369 | c rms-norm of ( e(i)/ewt(i) ) .le. 1, |
370 | c where ewt(i) = rtol(i)*abs(y(i)) + atol(i), |
371 | c and the rms-norm (root-mean-square norm) here is |
372 | c rms-norm(v) = sqrt(sum v(i)**2 / neq). here ewt = (ewt(i)) |
373 | c is a vector of weights which must always be positive, and |
374 | c the values of rtol and atol should all be non-negative. |
375 | c the following table gives the types (scalar/array) of |
376 | c rtol and atol, and the corresponding form of ewt(i). |
377 | c |
378 | c itol rtol atol ewt(i) |
379 | c 1 scalar scalar rtol*abs(y(i)) + atol |
380 | c 2 scalar array rtol*abs(y(i)) + atol(i) |
381 | c 3 array scalar rtol(i)*abs(y(i)) + atol |
382 | c 4 array array rtol(i)*abs(y(i)) + atol(i) |
383 | c |
384 | c when either of these parameters is a scalar, it need not |
385 | c be dimensioned in the user-s calling program. |
386 | c |
387 | c if none of the above choices (with itol, rtol, and atol |
388 | c fixed throughout the problem) is suitable, more general |
389 | c error controls can be obtained by substituting |
390 | c user-supplied routines for the setting of ewt and/or for |
391 | c the norm calculation. see part iv below. |
392 | c |
393 | c if global errors are to be estimated by making a repeated |
394 | c run on the same problem with smaller tolerances, then all |
395 | c components of rtol and atol (i.e. of ewt) should be scaled |
396 | c down uniformly. |
397 | c |
398 | c itask = an index specifying the task to be performed. |
399 | c input only. itask has the following values and meanings. |
400 | c 1 means normal computation of output values of y(t) at |
401 | c t = tout (by overshooting and interpolating). |
402 | c 2 means take one step only and return. |
403 | c 3 means stop at the first internal mesh point at or |
404 | c beyond t = tout and return. |
405 | c 4 means normal computation of output values of y(t) at |
406 | c t = tout but without overshooting t = tcrit. |
407 | c tcrit must be input as rwork(1). tcrit may be equal to |
408 | c or beyond tout, but not behind it in the direction of |
409 | c integration. this option is useful if the problem |
410 | c has a singularity at or beyond t = tcrit. |
411 | c 5 means take one step, without passing tcrit, and return. |
412 | c tcrit must be input as rwork(1). |
413 | c |
414 | c note.. if itask = 4 or 5 and the solver reaches tcrit |
415 | c (within roundoff), it will return t = tcrit (exactly) to |
416 | c indicate this (unless itask = 4 and tout comes before tcrit, |
417 | c in which case answers at t = tout are returned first). |
418 | c |
419 | c istate = an index used for input and output to specify the |
420 | c the state of the calculation. |
421 | c |
422 | c on input, the values of istate are as follows. |
423 | c 1 means this is the first call for the problem |
424 | c (initializations will be done). see note below. |
425 | c 2 means this is not the first call, and the calculation |
426 | c is to continue normally, with no change in any input |
427 | c parameters except possibly tout and itask. |
428 | c (if itol, rtol, and/or atol are changed between calls |
429 | c with istate = 2, the new values will be used but not |
430 | c tested for legality.) |
431 | c 3 means this is not the first call, and the |
432 | c calculation is to continue normally, but with |
433 | c a change in input parameters other than |
434 | c tout and itask. changes are allowed in |
435 | c neq, itol, rtol, atol, iopt, lrw, liw, mf, ml, mu, |
436 | c and any of the optional inputs except h0. |
437 | c (see iwork description for ml and mu.) |
438 | c note.. a preliminary call with tout = t is not counted |
439 | c as a first call here, as no initialization or checking of |
440 | c input is done. (such a call is sometimes useful for the |
441 | c purpose of outputting the initial conditions.) |
442 | c thus the first call for which tout .ne. t requires |
443 | c istate = 1 on input. |
444 | c |
445 | c on output, istate has the following values and meanings. |
446 | c 1 means nothing was done, as tout was equal to t with |
447 | c istate = 1 on input. (however, an internal counter was |
448 | c set to detect and prevent repeated calls of this type.) |
449 | c 2 means the integration was performed successfully. |
450 | c -1 means an excessive amount of work (more than mxstep |
451 | c steps) was done on this call, before completing the |
452 | c requested task, but the integration was otherwise |
453 | c successful as far as t. (mxstep is an optional input |
454 | c and is normally 500.) to continue, the user may |
455 | c simply reset istate to a value .gt. 1 and call again |
456 | c (the excess work step counter will be reset to 0). |
457 | c in addition, the user may increase mxstep to avoid |
458 | c this error return (see below on optional inputs). |
459 | c -2 means too much accuracy was requested for the precision |
460 | c of the machine being used. this was detected before |
461 | c completing the requested task, but the integration |
462 | c was successful as far as t. to continue, the tolerance |
463 | c parameters must be reset, and istate must be set |
464 | c to 3. the optional output tolsf may be used for this |
465 | c purpose. (note.. if this condition is detected before |
466 | c taking any steps, then an illegal input return |
467 | c (istate = -3) occurs instead.) |
468 | c -3 means illegal input was detected, before taking any |
469 | c integration steps. see written message for details. |
470 | c note.. if the solver detects an infinite loop of calls |
471 | c to the solver with illegal input, it will cause |
472 | c the run to stop. |
473 | c -4 means there were repeated error test failures on |
474 | c one attempted step, before completing the requested |
475 | c task, but the integration was successful as far as t. |
476 | c the problem may have a singularity, or the input |
477 | c may be inappropriate. |
478 | c -5 means there were repeated convergence test failures on |
479 | c one attempted step, before completing the requested |
480 | c task, but the integration was successful as far as t. |
481 | c this may be caused by an inaccurate jacobian matrix, |
482 | c if one is being used. |
483 | c -6 means ewt(i) became zero for some i during the |
484 | c integration. pure relative error control (atol(i)=0.0) |
485 | c was requested on a variable which has now vanished. |
486 | c the integration was successful as far as t. |
487 | c |
488 | c note.. since the normal output value of istate is 2, |
489 | c it does not need to be reset for normal continuation. |
490 | c also, since a negative input value of istate will be |
491 | c regarded as illegal, a negative output value requires the |
492 | c user to change it, and possibly other inputs, before |
493 | c calling the solver again. |
494 | c |
495 | c iopt = an integer flag to specify whether or not any optional |
496 | c inputs are being used on this call. input only. |
497 | c the optional inputs are listed separately below. |
498 | c iopt = 0 means no optional inputs are being used. |
499 | c default values will be used in all cases. |
500 | c iopt = 1 means one or more optional inputs are being used. |
501 | c |
502 | c rwork = a real working array (double precision). |
503 | c the length of rwork must be at least |
504 | c 20 + nyh*(maxord + 1) + 3*neq + lwm where |
505 | c nyh = the initial value of neq, |
506 | c maxord = 12 (if meth = 1) or 5 (if meth = 2) (unless a |
507 | c smaller value is given as an optional input), |
508 | c lwm = 0 if miter = 0, |
509 | c lwm = neq**2 + 2 if miter is 1 or 2, |
510 | c lwm = neq + 2 if miter = 3, and |
511 | c lwm = (2*ml+mu+1)*neq + 2 if miter is 4 or 5. |
512 | c (see the mf description for meth and miter.) |
513 | c thus if maxord has its default value and neq is constant, |
514 | c this length is.. |
515 | c 20 + 16*neq for mf = 10, |
516 | c 22 + 16*neq + neq**2 for mf = 11 or 12, |
517 | c 22 + 17*neq for mf = 13, |
518 | c 22 + 17*neq + (2*ml+mu)*neq for mf = 14 or 15, |
519 | c 20 + 9*neq for mf = 20, |
520 | c 22 + 9*neq + neq**2 for mf = 21 or 22, |
521 | c 22 + 10*neq for mf = 23, |
522 | c 22 + 10*neq + (2*ml+mu)*neq for mf = 24 or 25. |
523 | c the first 20 words of rwork are reserved for conditional |
524 | c and optional inputs and optional outputs. |
525 | c |
526 | c the following word in rwork is a conditional input.. |
527 | c rwork(1) = tcrit = critical value of t which the solver |
528 | c is not to overshoot. required if itask is |
529 | c 4 or 5, and ignored otherwise. (see itask.) |
530 | c |
531 | c lrw = the length of the array rwork, as declared by the user. |
532 | c (this will be checked by the solver.) |
533 | c |
534 | c iwork = an integer work array. the length of iwork must be at least |
535 | c 20 if miter = 0 or 3 (mf = 10, 13, 20, 23), or |
536 | c 20 + neq otherwise (mf = 11, 12, 14, 15, 21, 22, 24, 25). |
537 | c the first few words of iwork are used for conditional and |
538 | c optional inputs and optional outputs. |
539 | c |
540 | c the following 2 words in iwork are conditional inputs.. |
541 | c iwork(1) = ml these are the lower and upper |
542 | c iwork(2) = mu half-bandwidths, respectively, of the |
543 | c banded jacobian, excluding the main diagonal. |
544 | c the band is defined by the matrix locations |
545 | c (i,j) with i-ml .le. j .le. i+mu. ml and mu |
546 | c must satisfy 0 .le. ml,mu .le. neq-1. |
547 | c these are required if miter is 4 or 5, and |
548 | c ignored otherwise. ml and mu may in fact be |
549 | c the band parameters for a matrix to which |
550 | c df/dy is only approximately equal. |
551 | c |
552 | c liw = the length of the array iwork, as declared by the user. |
553 | c (this will be checked by the solver.) |
554 | c |
555 | c note.. the work arrays must not be altered between calls to lsode |
556 | c for the same problem, except possibly for the conditional and |
557 | c optional inputs, and except for the last 3*neq words of rwork. |
558 | c the latter space is used for internal scratch space, and so is |
559 | c available for use by the user outside lsode between calls, if |
560 | c desired (but not for use by f or jac). |
561 | c |
562 | c jac = the name of the user-supplied routine (miter = 1 or 4) to |
563 | c compute the jacobian matrix, df/dy, as a function of |
564 | c the scalar t and the vector y. it is to have the form |
565 | c subroutine jac (neq, t, y, ml, mu, pd, nrowpd) |
566 | c dimension y(1), pd(nrowpd,1) |
567 | c where neq, t, y, ml, mu, and nrowpd are input and the array |
568 | c pd is to be loaded with partial derivatives (elements of |
569 | c the jacobian matrix) on output. pd must be given a first |
570 | c dimension of nrowpd. t and y have the same meaning as in |
571 | c subroutine f. (in the dimension statement above, 1 is a |
572 | c dummy dimension.. it can be replaced by any value.) |
573 | c in the full matrix case (miter = 1), ml and mu are |
574 | c ignored, and the jacobian is to be loaded into pd in |
575 | c columnwise manner, with df(i)/dy(j) loaded into pd(i,j). |
576 | c in the band matrix case (miter = 4), the elements |
577 | c within the band are to be loaded into pd in columnwise |
578 | c manner, with diagonal lines of df/dy loaded into the rows |
579 | c of pd. thus df(i)/dy(j) is to be loaded into pd(i-j+mu+1,j). |
580 | c ml and mu are the half-bandwidth parameters (see iwork). |
581 | c the locations in pd in the two triangular areas which |
582 | c correspond to nonexistent matrix elements can be ignored |
583 | c or loaded arbitrarily, as they are overwritten by lsode. |
584 | c jac need not provide df/dy exactly. a crude |
585 | c approximation (possibly with a smaller bandwidth) will do. |
586 | c in either case, pd is preset to zero by the solver, |
587 | c so that only the nonzero elements need be loaded by jac. |
588 | c each call to jac is preceded by a call to f with the same |
589 | c arguments neq, t, and y. thus to gain some efficiency, |
590 | c intermediate quantities shared by both calculations may be |
591 | c saved in a user common block by f and not recomputed by jac, |
592 | c if desired. also, jac may alter the y array, if desired. |
593 | c jac must be declared external in the calling program. |
594 | c subroutine jac may access user-defined quantities in |
595 | c neq(2),... and/or in y(neq(1)+1),... if neq is an array |
596 | c (dimensioned in jac) and/or y has length exceeding neq(1). |
597 | c see the descriptions of neq and y above. |
598 | c |
599 | c mf = the method flag. used only for input. the legal values of |
600 | c mf are 10, 11, 12, 13, 14, 15, 20, 21, 22, 23, 24, and 25. |
601 | c mf has decimal digits meth and miter.. mf = 10*meth + miter. |
602 | c meth indicates the basic linear multistep method.. |
603 | c meth = 1 means the implicit adams method. |
604 | c meth = 2 means the method based on backward |
605 | c differentiation formulas (bdf-s). |
606 | c miter indicates the corrector iteration method.. |
607 | c miter = 0 means functional iteration (no jacobian matrix |
608 | c is involved). |
609 | c miter = 1 means chord iteration with a user-supplied |
610 | c full (neq by neq) jacobian. |
611 | c miter = 2 means chord iteration with an internally |
612 | c generated (difference quotient) full jacobian |
613 | c (using neq extra calls to f per df/dy value). |
614 | c miter = 3 means chord iteration with an internally |
615 | c generated diagonal jacobian approximation. |
616 | c (using 1 extra call to f per df/dy evaluation). |
617 | c miter = 4 means chord iteration with a user-supplied |
618 | c banded jacobian. |
619 | c miter = 5 means chord iteration with an internally |
620 | c generated banded jacobian (using ml+mu+1 extra |
621 | c calls to f per df/dy evaluation). |
622 | c if miter = 1 or 4, the user must supply a subroutine jac |
623 | c (the name is arbitrary) as described above under jac. |
624 | c for other values of miter, a dummy argument can be used. |
625 | c----------------------------------------------------------------------- |
626 | c optional inputs. |
627 | c |
628 | c the following is a list of the optional inputs provided for in the |
629 | c call sequence. (see also part ii.) for each such input variable, |
630 | c this table lists its name as used in this documentation, its |
631 | c location in the call sequence, its meaning, and the default value. |
632 | c the use of any of these inputs requires iopt = 1, and in that |
633 | c case all of these inputs are examined. a value of zero for any |
634 | c of these optional inputs will cause the default value to be used. |
635 | c thus to use a subset of the optional inputs, simply preload |
636 | c locations 5 to 10 in rwork and iwork to 0.0 and 0 respectively, and |
637 | c then set those of interest to nonzero values. |
638 | c |
639 | c name location meaning and default value |
640 | c |
641 | c h0 rwork(5) the step size to be attempted on the first step. |
642 | c the default value is determined by the solver. |
643 | c |
644 | c hmax rwork(6) the maximum absolute step size allowed. |
645 | c the default value is infinite. |
646 | c |
647 | c hmin rwork(7) the minimum absolute step size allowed. |
648 | c the default value is 0. (this lower bound is not |
649 | c enforced on the final step before reaching tcrit |
650 | c when itask = 4 or 5.) |
651 | c |
652 | c maxord iwork(5) the maximum order to be allowed. the default |
653 | c value is 12 if meth = 1, and 5 if meth = 2. |
654 | c if maxord exceeds the default value, it will |
655 | c be reduced to the default value. |
656 | c if maxord is changed during the problem, it may |
657 | c cause the current order to be reduced. |
658 | c |
659 | c mxstep iwork(6) maximum number of (internally defined) steps |
660 | c allowed during one call to the solver. |
661 | c the default value is 500. |
662 | c |
663 | c mxhnil iwork(7) maximum number of messages printed (per problem) |
664 | c warning that t + h = t on a step (h = step size). |
665 | c this must be positive to result in a non-default |
666 | c value. the default value is 10. |
667 | c----------------------------------------------------------------------- |
668 | c optional outputs. |
669 | c |
670 | c as optional additional output from lsode, the variables listed |
671 | c below are quantities related to the performance of lsode |
672 | c which are available to the user. these are communicated by way of |
673 | c the work arrays, but also have internal mnemonic names as shown. |
674 | c except where stated otherwise, all of these outputs are defined |
675 | c on any successful return from lsode, and on any return with |
676 | c istate = -1, -2, -4, -5, or -6. on an illegal input return |
677 | c (istate = -3), they will be unchanged from their existing values |
678 | c (if any), except possibly for tolsf, lenrw, and leniw. |
679 | c on any error return, outputs relevant to the error will be defined, |
680 | c as noted below. |
681 | c |
682 | c name location meaning |
683 | c |
684 | c hu rwork(11) the step size in t last used (successfully). |
685 | c |
686 | c hcur rwork(12) the step size to be attempted on the next step. |
687 | c |
688 | c tcur rwork(13) the current value of the independent variable |
689 | c which the solver has actually reached, i.e. the |
690 | c current internal mesh point in t. on output, tcur |
691 | c will always be at least as far as the argument |
692 | c t, but may be farther (if interpolation was done). |
693 | c |
694 | c tolsf rwork(14) a tolerance scale factor, greater than 1.0, |
695 | c computed when a request for too much accuracy was |
696 | c detected (istate = -3 if detected at the start of |
697 | c the problem, istate = -2 otherwise). if itol is |
698 | c left unaltered but rtol and atol are uniformly |
699 | c scaled up by a factor of tolsf for the next call, |
700 | c then the solver is deemed likely to succeed. |
701 | c (the user may also ignore tolsf and alter the |
702 | c tolerance parameters in any other way appropriate.) |
703 | c |
704 | c nst iwork(11) the number of steps taken for the problem so far. |
705 | c |
706 | c nfe iwork(12) the number of f evaluations for the problem so far. |
707 | c |
708 | c nje iwork(13) the number of jacobian evaluations (and of matrix |
709 | c lu decompositions) for the problem so far. |
710 | c |
711 | c nqu iwork(14) the method order last used (successfully). |
712 | c |
713 | c nqcur iwork(15) the order to be attempted on the next step. |
714 | c |
715 | c imxer iwork(16) the index of the component of largest magnitude in |
716 | c the weighted local error vector ( e(i)/ewt(i) ), |
717 | c on an error return with istate = -4 or -5. |
718 | c |
719 | c lenrw iwork(17) the length of rwork actually required. |
720 | c this is defined on normal returns and on an illegal |
721 | c input return for insufficient storage. |
722 | c |
723 | c leniw iwork(18) the length of iwork actually required. |
724 | c this is defined on normal returns and on an illegal |
725 | c input return for insufficient storage. |
726 | c |
727 | c the following two arrays are segments of the rwork array which |
728 | c may also be of interest to the user as optional outputs. |
729 | c for each array, the table below gives its internal name, |
730 | c its base address in rwork, and its description. |
731 | c |
732 | c name base address description |
733 | c |
734 | c yh 21 the nordsieck history array, of size nyh by |
735 | c (nqcur + 1), where nyh is the initial value |
736 | c of neq. for j = 0,1,...,nqcur, column j+1 |
737 | c of yh contains hcur**j/factorial(j) times |
738 | c the j-th derivative of the interpolating |
739 | c polynomial currently representing the solution, |
740 | c evaluated at t = tcur. |
741 | c |
742 | c acor lenrw-neq+1 array of size neq used for the accumulated |
743 | c corrections on each step, scaled on output |
744 | c to represent the estimated local error in y |
745 | c on the last step. this is the vector e in |
746 | c the description of the error control. it is |
747 | c defined only on a successful return from lsode. |
748 | c |
749 | c----------------------------------------------------------------------- |
750 | c part ii. other routines callable. |
751 | c |
752 | c the following are optional calls which the user may make to |
753 | c gain additional capabilities in conjunction with lsode. |
754 | c (the routines xsetun and xsetf are designed to conform to the |
755 | c slatec error handling package.) |
756 | c |
757 | c form of call function |
758 | c call xsetun(lun) set the logical unit number, lun, for |
759 | c output of messages from lsode, if |
760 | c the default is not desired. |
761 | c the default value of lun is 6. |
762 | c |
763 | c call xsetf(mflag) set a flag to control the printing of |
764 | c messages by lsode. |
765 | c mflag = 0 means do not print. (danger.. |
766 | c this risks losing valuable information.) |
767 | c mflag = 1 means print (the default). |
768 | c |
769 | c either of the above calls may be made at |
770 | c any time and will take effect immediately. |
771 | c |
772 | c call srcom(rsav,isav,job) saves and restores the contents of |
773 | c the internal common blocks used by |
774 | c lsode (see part iii below). |
775 | c rsav must be a real array of length 218 |
776 | c or more, and isav must be an integer |
777 | c array of length 41 or more. |
778 | c job=1 means save common into rsav/isav. |
779 | c job=2 means restore common from rsav/isav. |
780 | c srcom is useful if one is |
781 | c interrupting a run and restarting |
782 | c later, or alternating between two or |
783 | c more problems solved with lsode. |
784 | c |
785 | c call intdy(,,,,,) provide derivatives of y, of various |
786 | c (see below) orders, at a specified point t, if |
787 | c desired. it may be called only after |
788 | c a successful return from lsode. |
789 | c |
790 | c the detailed instructions for using intdy are as follows. |
791 | c the form of the call is.. |
792 | c |
793 | c call intdy (t, k, rwork(21), nyh, dky, iflag) |
794 | c |
795 | c the input parameters are.. |
796 | c |
797 | c t = value of independent variable where answers are desired |
798 | c (normally the same as the t last returned by lsode). |
799 | c for valid results, t must lie between tcur - hu and tcur. |
800 | c (see optional outputs for tcur and hu.) |
801 | c k = integer order of the derivative desired. k must satisfy |
802 | c 0 .le. k .le. nqcur, where nqcur is the current order |
803 | c (see optional outputs). the capability corresponding |
804 | c to k = 0, i.e. computing y(t), is already provided |
805 | c by lsode directly. since nqcur .ge. 1, the first |
806 | c derivative dy/dt is always available with intdy. |
807 | c rwork(21) = the base address of the history array yh. |
808 | c nyh = column length of yh, equal to the initial value of neq. |
809 | c |
810 | c the output parameters are.. |
811 | c |
812 | c dky = a real array of length neq containing the computed value |
813 | c of the k-th derivative of y(t). |
814 | c iflag = integer flag, returned as 0 if k and t were legal, |
815 | c -1 if k was illegal, and -2 if t was illegal. |
816 | c on an error return, a message is also written. |
817 | c----------------------------------------------------------------------- |
818 | c part iii. common blocks. |
819 | c |
820 | c if lsode is to be used in an overlay situation, the user |
821 | c must declare, in the primary overlay, the variables in.. |
822 | c (1) the call sequence to lsode, |
823 | c (2) the two internal common blocks |
824 | c /ls0001/ of length 257 (218 double precision words |
825 | c followed by 39 integer words), |
826 | c /eh0001/ of length 2 (integer words). |
827 | c |
828 | c if lsode is used on a system in which the contents of internal |
829 | c common blocks are not preserved between calls, the user should |
830 | c declare the above two common blocks in his main program to insure |
831 | c that their contents are preserved. |
832 | c |
833 | c if the solution of a given problem by lsode is to be interrupted |
834 | c and then later continued, such as when restarting an interrupted run |
835 | c or alternating between two or more problems, the user should save, |
836 | c following the return from the last lsode call prior to the |
837 | c interruption, the contents of the call sequence variables and the |
838 | c internal common blocks, and later restore these values before the |
839 | c next lsode call for that problem. to save and restore the common |
840 | c blocks, use subroutine srcom (see part ii above). |
841 | c |
842 | c----------------------------------------------------------------------- |
843 | c part iv. optionally replaceable solver routines. |
844 | c |
845 | c below are descriptions of two routines in the lsode package which |
846 | c relate to the measurement of errors. either routine can be |
847 | c replaced by a user-supplied version, if desired. however, since such |
848 | c a replacement may have a major impact on performance, it should be |
849 | c done only when absolutely necessary, and only with great caution. |
850 | c (note.. the means by which the package version of a routine is |
851 | c superseded by the user-s version may be system-dependent.) |
852 | c |
853 | c (a) ewset. |
854 | c the following subroutine is called just before each internal |
855 | c integration step, and sets the array of error weights, ewt, as |
856 | c described under itol/rtol/atol above.. |
857 | c subroutine ewset (neq, itol, rtol, atol, ycur, ewt) |
858 | c where neq, itol, rtol, and atol are as in the lsode call sequence, |
859 | c ycur contains the current dependent variable vector, and |
860 | c ewt is the array of weights set by ewset. |
861 | c |
862 | c if the user supplies this subroutine, it must return in ewt(i) |
863 | c (i = 1,...,neq) a positive quantity suitable for comparing errors |
864 | c in y(i) to. the ewt array returned by ewset is passed to the |
865 | c vnorm routine (see below), and also used by lsode in the computation |
866 | c of the optional output imxer, the diagonal jacobian approximation, |
867 | c and the increments for difference quotient jacobians. |
868 | c |
869 | c in the user-supplied version of ewset, it may be desirable to use |
870 | c the current values of derivatives of y. derivatives up to order nq |
871 | c are available from the history array yh, described above under |
872 | c optional outputs. in ewset, yh is identical to the ycur array, |
873 | c extended to nq + 1 columns with a column length of nyh and scale |
874 | c factors of h**j/factorial(j). on the first call for the problem, |
875 | c given by nst = 0, nq is 1 and h is temporarily set to 1.0. |
876 | c the quantities nq, nyh, h, and nst can be obtained by including |
877 | c in ewset the statements.. |
878 | c double precision h, rls |
879 | c common /ls0001/ rls(218),ils(39) |
880 | c nq = ils(35) |
881 | c nyh = ils(14) |
882 | c nst = ils(36) |
883 | c h = rls(212) |
884 | c thus, for example, the current value of dy/dt can be obtained as |
885 | c ycur(nyh+i)/h (i=1,...,neq) (and the division by h is |
886 | c unnecessary when nst = 0). |
887 | c |
888 | c (b) vnorm. |
889 | c the following is a real function routine which computes the weighted |
890 | c root-mean-square norm of a vector v.. |
891 | c d = vnorm (n, v, w) |
892 | c where.. |
893 | c n = the length of the vector, |
894 | c v = real array of length n containing the vector, |
895 | c w = real array of length n containing weights, |
896 | c d = sqrt( (1/n) * sum(v(i)*w(i))**2 ). |
897 | c vnorm is called with n = neq and with w(i) = 1.0/ewt(i), where |
898 | c ewt is as set by subroutine ewset. |
899 | c |
900 | c if the user supplies this function, it should return a non-negative |
901 | c value of vnorm suitable for use in the error control in lsode. |
902 | c none of the arguments should be altered by vnorm. |
903 | c for example, a user-supplied vnorm routine might.. |
904 | c -substitute a max-norm of (v(i)*w(i)) for the rms-norm, or |
905 | c -ignore some components of v in the norm, with the effect of |
906 | c suppressing the error control on those components of y. |
907 | c----------------------------------------------------------------------- |
908 | c----------------------------------------------------------------------- |
909 | c other routines in the lsode package. |
910 | c |
911 | c in addition to subroutine lsode, the lsode package includes the |
912 | c following subroutines and function routines.. |
913 | c intdy computes an interpolated value of the y vector at t = tout. |
914 | c stode is the core integrator, which does one step of the |
915 | c integration and the associated error control. |
916 | c cfode sets all method coefficients and test constants. |
917 | c prepj computes and preprocesses the jacobian matrix j = df/dy |
918 | c and the newton iteration matrix p = i - h*l0*j. |
919 | c solsy manages solution of linear system in chord iteration. |
920 | c ewset sets the error weight vector ewt before each step. |
921 | c vnorm computes the weighted r.m.s. norm of a vector. |
922 | c srcom is a user-callable routine to save and restore |
923 | c the contents of the internal common blocks. |
924 | c dgefa and dgesl are routines from linpack for solving full |
925 | c systems of linear algebraic equations. |
926 | c dgbfa and dgbsl are routines from linpack for solving banded |
927 | c linear systems. |
928 | c daxpy, dscal, idamax, and ddot are basic linear algebra modules |
929 | c (blas) used by the above linpack routines. |
930 | c d1mach computes the unit roundoff in a machine-independent manner. |
931 | c XERRWV, xsetun, and xsetf handle the printing of all error |
932 | c messages and warnings. XERRWV is machine-dependent. |
933 | C ASCEND C XERRWV has been replaced by XASCWV from interface/Lsode.c. |
934 | c note.. vnorm, idamax, ddot, and d1mach are function routines. |
935 | c all the others are subroutines. |
936 | c |
937 | c the intrinsic and external routines used by lsode are.. |
938 | c dabs, dmax1, dmin1, dfloat, max0, min0, mod, dsign, dsqrt, and write. |
939 | c |
940 | c a block data subprogram is also included with the package, |
941 | c for loading some of the variables in internal common. |
942 | c |
943 | c----------------------------------------------------------------------- |
944 | c the following card is for optimized compilation on llnl compilers. |
945 | clll. optimize |
946 | c----------------------------------------------------------------------- |
947 | cascend changes |
948 | Ckaa external aftime |
949 | external prepj, solsy |
950 | integer illin, init, lyh, lewt, lacor, lsavf, lwm, liwm, |
951 | 1 mxstep, mxhnil, nhnil, ntrep, nslast, nyh, iowns |
952 | integer icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, |
953 | 1 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu |
954 | integer i, i1, i2, iflag, imxer, kgo, lf0, |
955 | 1 leniw, lenrw, lenwm, ml, mord, mu, mxhnl0, mxstp0 |
956 | double precision rowns, |
957 | 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround |
958 | double precision atoli, ayi, big, ewti, h0, hmax, hmx, rh, rtoli, |
959 | 1 tcrit, tdist, tnext, tol, tolsf, tp, size, sum, w0, |
960 | 2 d1mach, vnorm |
961 | dimension mord(2) |
962 | logical ihit |
963 | c----------------------------------------------------------------------- |
964 | c the following internal common block contains |
965 | c (a) variables which are local to any subroutine but whose values must |
966 | c be preserved between calls to the routine (own variables), and |
967 | c (b) variables which are communicated between subroutines. |
968 | c the structure of the block is as follows.. all real variables are |
969 | c listed first, followed by all integers. within each type, the |
970 | c variables are grouped with those local to subroutine lsode first, |
971 | c then those local to subroutine stode, and finally those used |
972 | c for communication. the block is declared in subroutines |
973 | c lsode, intdy, stode, prepj, and solsy. groups of variables are |
974 | c replaced by dummy arrays in the common declarations in routines |
975 | c where those variables are not used. |
976 | c----------------------------------------------------------------------- |
977 | common /ls0001/ rowns(209), |
978 | 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround, |
979 | 2 illin, init, lyh, lewt, lacor, lsavf, lwm, liwm, |
980 | 3 mxstep, mxhnil, nhnil, ntrep, nslast, nyh, iowns(6), |
981 | 4 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, |
982 | 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu |
983 | c |
984 | data mord(1),mord(2)/12,5/, mxstp0/500/, mxhnl0/10/ |
985 | |
986 | C ascend debugging options |
987 | C type *, 'iopt',iopt |
988 | C type *, 'h0',rwork(5) |
989 | C type *, 'hmax',rwork(6) |
990 | C type *, 'hmin',rwork(7) |
991 | C type *, 'maxs',iwork(6) |
992 | c----------------------------------------------------------------------- |
993 | c block a. |
994 | c this code block is executed on every call. |
995 | c it tests istate and itask for legality and branches appropriately. |
996 | c if istate .gt. 1 but the flag init shows that initialization has |
997 | c not yet been done, an error return occurs. |
998 | c if istate = 1 and tout = t, jump to block g and return immediately. |
999 | c----------------------------------------------------------------------- |
1000 | if (istate .lt. 1 .or. istate .gt. 3) go to 601 |
1001 | if (itask .lt. 1 .or. itask .gt. 5) go to 602 |
1002 | if (istate .eq. 1) go to 10 |
1003 | if (init .eq. 0) go to 603 |
1004 | if (istate .eq. 2) go to 200 |
1005 | go to 20 |
1006 | 10 init = 0 |
1007 | if (tout .eq. t) go to 430 |
1008 | 20 ntrep = 0 |
1009 | c----------------------------------------------------------------------- |
1010 | c block b. |
1011 | c the next code block is executed for the initial call (istate = 1), |
1012 | c or for a continuation call with parameter changes (istate = 3). |
1013 | c it contains checking of all inputs and various initializations. |
1014 | c |
1015 | c first check legality of the non-optional inputs neq, itol, iopt, |
1016 | c mf, ml, and mu. |
1017 | c----------------------------------------------------------------------- |
1018 | if (neq(1) .le. 0) go to 604 |
1019 | if (istate .eq. 1) go to 25 |
1020 | if (neq(1) .gt. n) go to 605 |
1021 | 25 n = neq(1) |
1022 | if (itol .lt. 1 .or. itol .gt. 4) go to 606 |
1023 | if (iopt .lt. 0 .or. iopt .gt. 1) go to 607 |
1024 | meth = mf/10 |
1025 | miter = mf - 10*meth |
1026 | if (meth .lt. 1 .or. meth .gt. 2) go to 608 |
1027 | if (miter .lt. 0 .or. miter .gt. 5) go to 608 |
1028 | if (miter .le. 3) go to 30 |
1029 | ml = iwork(1) |
1030 | mu = iwork(2) |
1031 | if (ml .lt. 0 .or. ml .ge. n) go to 609 |
1032 | if (mu .lt. 0 .or. mu .ge. n) go to 610 |
1033 | 30 continue |
1034 | c next process and check the optional inputs. -------------------------- |
1035 | if (iopt .eq. 1) go to 40 |
1036 | maxord = mord(meth) |
1037 | mxstep = mxstp0 |
1038 | mxhnil = mxhnl0 |
1039 | if (istate .eq. 1) h0 = 0.0d0 |
1040 | hmxi = 0.0d0 |
1041 | hmin = 0.0d0 |
1042 | go to 60 |
1043 | 40 maxord = iwork(5) |
1044 | if (maxord .lt. 0) go to 611 |
1045 | if (maxord .eq. 0) maxord = 100 |
1046 | maxord = min0(maxord,mord(meth)) |
1047 | mxstep = iwork(6) |
1048 | if (mxstep .lt. 0) go to 612 |
1049 | if (mxstep .eq. 0) mxstep = mxstp0 |
1050 | mxhnil = iwork(7) |
1051 | if (mxhnil .lt. 0) go to 613 |
1052 | if (mxhnil .eq. 0) mxhnil = mxhnl0 |
1053 | if (istate .ne. 1) go to 50 |
1054 | h0 = rwork(5) |
1055 | if ((tout - t)*h0 .lt. 0.0d0) go to 614 |
1056 | 50 hmax = rwork(6) |
1057 | if (hmax .lt. 0.0d0) go to 615 |
1058 | hmxi = 0.0d0 |
1059 | if (hmax .gt. 0.0d0) hmxi = 1.0d0/hmax |
1060 | hmin = rwork(7) |
1061 | if (hmin .lt. 0.0d0) go to 616 |
1062 | c----------------------------------------------------------------------- |
1063 | c set work array pointers and check lengths lrw and liw. |
1064 | c pointers to segments of rwork and iwork are named by prefixing l to |
1065 | c the name of the segment. e.g., the segment yh starts at rwork(lyh). |
1066 | c segments of rwork (in order) are denoted yh, wm, ewt, savf, acor. |
1067 | c----------------------------------------------------------------------- |
1068 | 60 lyh = 21 |
1069 | if (istate .eq. 1) nyh = n |
1070 | lwm = lyh + (maxord + 1)*nyh |
1071 | if (miter .eq. 0) lenwm = 0 |
1072 | if (miter .eq. 1 .or. miter .eq. 2) lenwm = n*n + 2 |
1073 | if (miter .eq. 3) lenwm = n + 2 |
1074 | if (miter .ge. 4) lenwm = (2*ml + mu + 1)*n + 2 |
1075 | lewt = lwm + lenwm |
1076 | lsavf = lewt + n |
1077 | lacor = lsavf + n |
1078 | lenrw = lacor + n - 1 |
1079 | iwork(17) = lenrw |
1080 | liwm = 1 |
1081 | leniw = 20 + n |
1082 | if (miter .eq. 0 .or. miter .eq. 3) leniw = 20 |
1083 | iwork(18) = leniw |
1084 | if (lenrw .gt. lrw) go to 617 |
1085 | if (leniw .gt. liw) go to 618 |
1086 | c check rtol and atol for legality. ------------------------------------ |
1087 | rtoli = rtol(1) |
1088 | atoli = atol(1) |
1089 | do 70 i = 1,n |
1090 | if (itol .ge. 3) rtoli = rtol(i) |
1091 | if (itol .eq. 2 .or. itol .eq. 4) atoli = atol(i) |
1092 | if (rtoli .lt. 0.0d0) go to 619 |
1093 | if (atoli .lt. 0.0d0) go to 620 |
1094 | 70 continue |
1095 | if (istate .eq. 1) go to 100 |
1096 | c if istate = 3, set flag to signal parameter changes to stode. -------- |
1097 | jstart = -1 |
1098 | if (nq .le. maxord) go to 90 |
1099 | c maxord was reduced below nq. copy yh(*,maxord+2) into savf. --------- |
1100 | do 80 i = 1,n |
1101 | 80 rwork(i+lsavf-1) = rwork(i+lwm-1) |
1102 | c reload wm(1) = rwork(lwm), since lwm may have changed. --------------- |
1103 | 90 if (miter .gt. 0) rwork(lwm) = dsqrt(uround) |
1104 | if (n .eq. nyh) go to 200 |
1105 | c neq was reduced. zero part of yh to avoid undefined references. ----- |
1106 | i1 = lyh + l*nyh |
1107 | i2 = lyh + (maxord + 1)*nyh - 1 |
1108 | if (i1 .gt. i2) go to 200 |
1109 | do 95 i = i1,i2 |
1110 | 95 rwork(i) = 0.0d0 |
1111 | go to 200 |
1112 | c----------------------------------------------------------------------- |
1113 | c block c. |
1114 | c the next block is for the initial call only (istate = 1). |
1115 | c it contains all remaining initializations, the initial call to f, |
1116 | c and the calculation of the initial step size. |
1117 | c the error weights in ewt are inverted after being loaded. |
1118 | c----------------------------------------------------------------------- |
1119 | 100 uround = d1mach(4) |
1120 | tn = t |
1121 | if (itask .ne. 4 .and. itask .ne. 5) go to 110 |
1122 | tcrit = rwork(1) |
1123 | if ((tcrit - tout)*(tout - t) .lt. 0.0d0) go to 625 |
1124 | if (h0 .ne. 0.0d0 .and. (t + h0 - tcrit)*h0 .gt. 0.0d0) |
1125 | 1 h0 = tcrit - t |
1126 | 110 jstart = 0 |
1127 | if (miter .gt. 0) rwork(lwm) = dsqrt(uround) |
1128 | nhnil = 0 |
1129 | nst = 0 |
1130 | nje = 0 |
1131 | nslast = 0 |
1132 | hu = 0.0d0 |
1133 | nqu = 0 |
1134 | ccmax = 0.3d0 |
1135 | maxcor = 3 |
1136 | msbp = 20 |
1137 | mxncf = 10 |
1138 | c initial call to f. (lf0 points to yh(*,2).) ------------------------- |
1139 | lf0 = lyh + nyh |
1140 | call f (neq, t, y, rwork(lf0)) |
1141 | nfe = 1 |
1142 | c load the initial value vector in yh. --------------------------------- |
1143 | do 115 i = 1,n |
1144 | 115 rwork(i+lyh-1) = y(i) |
1145 | c load and invert the ewt array. (h is temporarily set to 1.0.) ------- |
1146 | nq = 1 |
1147 | h = 1.0d0 |
1148 | call ewset (n, itol, rtol, atol, rwork(lyh), rwork(lewt)) |
1149 | do 120 i = 1,n |
1150 | if (rwork(i+lewt-1) .le. 0.0d0) go to 621 |
1151 | 120 rwork(i+lewt-1) = 1.0d0/rwork(i+lewt-1) |
1152 | c----------------------------------------------------------------------- |
1153 | c the coding below computes the step size, h0, to be attempted on the |
1154 | c first step, unless the user has supplied a value for this. |
1155 | c first check that tout - t differs significantly from zero. |
1156 | c a scalar tolerance quantity tol is computed, as max(rtol(i)) |
1157 | c if this is positive, or max(atol(i)/abs(y(i))) otherwise, adjusted |
1158 | c so as to be between 100*uround and 1.0e-3. |
1159 | c then the computed value h0 is given by.. |
1160 | c neq |
1161 | c h0**2 = tol / ( w0**-2 + (1/neq) * sum ( f(i)/ywt(i) )**2 ) |
1162 | c 1 |
1163 | c where w0 = max ( abs(t), abs(tout) ), |
1164 | c f(i) = i-th component of initial value of f, |
1165 | c ywt(i) = ewt(i)/tol (a weight for y(i)). |
1166 | c the sign of h0 is inferred from the initial values of tout and t. |
1167 | c----------------------------------------------------------------------- |
1168 | if (h0 .ne. 0.0d0) go to 180 |
1169 | tdist = dabs(tout - t) |
1170 | w0 = dmax1(dabs(t),dabs(tout)) |
1171 | if (tdist .lt. 2.0d0*uround*w0) go to 622 |
1172 | tol = rtol(1) |
1173 | if (itol .le. 2) go to 140 |
1174 | do 130 i = 1,n |
1175 | 130 tol = dmax1(tol,rtol(i)) |
1176 | 140 if (tol .gt. 0.0d0) go to 160 |
1177 | atoli = atol(1) |
1178 | do 150 i = 1,n |
1179 | if (itol .eq. 2 .or. itol .eq. 4) atoli = atol(i) |
1180 | ayi = dabs(y(i)) |
1181 | if (ayi .ne. 0.0d0) tol = dmax1(tol,atoli/ayi) |
1182 | 150 continue |
1183 | 160 tol = dmax1(tol,100.0d0*uround) |
1184 | tol = dmin1(tol,0.001d0) |
1185 | sum = vnorm (n, rwork(lf0), rwork(lewt)) |
1186 | sum = 1.0d0/(tol*w0*w0) + tol*sum**2 |
1187 | h0 = 1.0d0/dsqrt(sum) |
1188 | h0 = dmin1(h0,tdist) |
1189 | h0 = dsign(h0,tout-t) |
1190 | c adjust h0 if necessary to meet hmax bound. --------------------------- |
1191 | 180 rh = dabs(h0)*hmxi |
1192 | if (rh .gt. 1.0d0) h0 = h0/rh |
1193 | c load h with h0 and scale yh(*,2) by h0. ------------------------------ |
1194 | h = h0 |
1195 | do 190 i = 1,n |
1196 | 190 rwork(i+lf0-1) = h0*rwork(i+lf0-1) |
1197 | go to 270 |
1198 | c----------------------------------------------------------------------- |
1199 | c block d. |
1200 | c the next code block is for continuation calls only (istate = 2 or 3) |
1201 | c and is to check stop conditions before taking a step. |
1202 | c----------------------------------------------------------------------- |
1203 | 200 nslast = nst |
1204 | go to (210, 250, 220, 230, 240), itask |
1205 | 210 if ((tn - tout)*h .lt. 0.0d0) go to 250 |
1206 | call intdy (tout, 0, rwork(lyh), nyh, y, iflag) |
1207 | if (iflag .ne. 0) go to 627 |
1208 | t = tout |
1209 | go to 420 |
1210 | 220 tp = tn - hu*(1.0d0 + 100.0d0*uround) |
1211 | if ((tp - tout)*h .gt. 0.0d0) go to 623 |
1212 | if ((tn - tout)*h .lt. 0.0d0) go to 250 |
1213 | go to 400 |
1214 | 230 tcrit = rwork(1) |
1215 | if ((tn - tcrit)*h .gt. 0.0d0) go to 624 |
1216 | if ((tcrit - tout)*h .lt. 0.0d0) go to 625 |
1217 | if ((tn - tout)*h .lt. 0.0d0) go to 245 |
1218 | call intdy (tout, 0, rwork(lyh), nyh, y, iflag) |
1219 | if (iflag .ne. 0) go to 627 |
1220 | t = tout |
1221 | go to 420 |
1222 | 240 tcrit = rwork(1) |
1223 | if ((tn - tcrit)*h .gt. 0.0d0) go to 624 |
1224 | 245 hmx = dabs(tn) + dabs(h) |
1225 | ihit = dabs(tn - tcrit) .le. 100.0d0*uround*hmx |
1226 | if (ihit) go to 400 |
1227 | tnext = tn + h*(1.0d0 + 4.0d0*uround) |
1228 | if ((tnext - tcrit)*h .le. 0.0d0) go to 250 |
1229 | h = (tcrit - tn)*(1.0d0 - 4.0d0*uround) |
1230 | if (istate .eq. 2) jstart = -2 |
1231 | c----------------------------------------------------------------------- |
1232 | c block e. |
1233 | c the next block is normally executed for all calls and contains |
1234 | c the call to the one-step core integrator stode. |
1235 | c |
1236 | c this is a looping point for the integration steps. |
1237 | c |
1238 | c first check for too many steps being taken, update ewt (if not at |
1239 | c start of problem), check for too much accuracy being requested, and |
1240 | c check for h below the roundoff level in t. |
1241 | c----------------------------------------------------------------------- |
1242 | 250 continue |
1243 | if ((nst-nslast) .ge. mxstep) go to 500 |
1244 | call ewset (n, itol, rtol, atol, rwork(lyh), rwork(lewt)) |
1245 | do 260 i = 1,n |
1246 | if (rwork(i+lewt-1) .le. 0.0d0) go to 510 |
1247 | 260 rwork(i+lewt-1) = 1.0d0/rwork(i+lewt-1) |
1248 | 270 tolsf = uround*vnorm (n, rwork(lyh), rwork(lewt)) |
1249 | if (tolsf .le. 1.0d0) go to 280 |
1250 | tolsf = tolsf*2.0d0 |
1251 | if (nst .eq. 0) go to 626 |
1252 | go to 520 |
1253 | 280 if ((tn + h) .ne. tn) go to 290 |
1254 | nhnil = nhnil + 1 |
1255 | if (nhnil .gt. mxhnil) go to 290 |
1256 | call xascwv(50hlsode-- warning..internal t (=r1) and h (=r2) are, |
1257 | 1 50, 101, 0, 0, 0, 0, 0, 0.0d0, 0.0d0) |
1258 | call xascwv( |
1259 | 1 60h such that in the machine, t + h = t on the next step , |
1260 | 1 60, 101, 0, 0, 0, 0, 0, 0.0d0, 0.0d0) |
1261 | call xascwv(50h (h = step size). solver will continue anyway, |
1262 | 1 50, 101, 0, 0, 0, 0, 2, tn, h) |
1263 | if (nhnil .lt. mxhnil) go to 290 |
1264 | call xascwv(50hlsode-- above warning has been issued i1 times. , |
1265 | 1 50, 102, 0, 0, 0, 0, 0, 0.0d0, 0.0d0) |
1266 | call xascwv(50h it will not be issued again for this problem, |
1267 | 1 50, 102, 0, 1, mxhnil, 0, 0, 0.0d0, 0.0d0) |
1268 | 290 continue |
1269 | c----------------------------------------------------------------------- |
1270 | c call stode(neq,y,yh,nyh,yh,ewt,savf,acor,wm,iwm,f,jac,prepj,solsy) |
1271 | c----------------------------------------------------------------------- |
1272 | c boyd |
1273 | call stode (neq, y, rwork(lyh), nyh, rwork(lyh), rwork(lewt), |
1274 | 1 rwork(lsavf), rwork(lacor), rwork(lwm), iwork(liwm), |
1275 | 2 f, jac, prepj, solsy) |
1276 | kgo = 1 - kflag |
1277 | go to (300, 530, 540), kgo |
1278 | c----------------------------------------------------------------------- |
1279 | c block f. |
1280 | c the following block handles the case of a successful return from the |
1281 | c core integrator (kflag = 0). test for stop conditions. |
1282 | c----------------------------------------------------------------------- |
1283 | 300 init = 1 |
1284 | go to (310, 400, 330, 340, 350), itask |
1285 | c itask = 1. if tout has been reached, interpolate. ------------------- |
1286 | 310 if ((tn - tout)*h .lt. 0.0d0) go to 250 |
1287 | call intdy (tout, 0, rwork(lyh), nyh, y, iflag) |
1288 | t = tout |
1289 | go to 420 |
1290 | c itask = 3. jump to exit if tout was reached. ------------------------ |
1291 | 330 if ((tn - tout)*h .ge. 0.0d0) go to 400 |
1292 | go to 250 |
1293 | c itask = 4. see if tout or tcrit was reached. adjust h if necessary. |
1294 | 340 if ((tn - tout)*h .lt. 0.0d0) go to 345 |
1295 | call intdy (tout, 0, rwork(lyh), nyh, y, iflag) |
1296 | t = tout |
1297 | go to 420 |
1298 | 345 hmx = dabs(tn) + dabs(h) |
1299 | ihit = dabs(tn - tcrit) .le. 100.0d0*uround*hmx |
1300 | if (ihit) go to 400 |
1301 | tnext = tn + h*(1.0d0 + 4.0d0*uround) |
1302 | if ((tnext - tcrit)*h .le. 0.0d0) go to 250 |
1303 | h = (tcrit - tn)*(1.0d0 - 4.0d0*uround) |
1304 | jstart = -2 |
1305 | go to 250 |
1306 | c itask = 5. see if tcrit was reached and jump to exit. --------------- |
1307 | 350 hmx = dabs(tn) + dabs(h) |
1308 | ihit = dabs(tn - tcrit) .le. 100.0d0*uround*hmx |
1309 | c----------------------------------------------------------------------- |
1310 | c block g. |
1311 | c the following block handles all successful returns from lsode. |
1312 | c if itask .ne. 1, y is loaded from yh and t is set accordingly. |
1313 | c istate is set to 2, the illegal input counter is zeroed, and the |
1314 | c optional outputs are loaded into the work arrays before returning. |
1315 | c if istate = 1 and tout = t, there is a return with no action taken, |
1316 | c except that if this has happened repeatedly, the run is terminated. |
1317 | c----------------------------------------------------------------------- |
1318 | 400 do 410 i = 1,n |
1319 | 410 y(i) = rwork(i+lyh-1) |
1320 | t = tn |
1321 | if (itask .ne. 4 .and. itask .ne. 5) go to 420 |
1322 | if (ihit) t = tcrit |
1323 | 420 istate = 2 |
1324 | illin = 0 |
1325 | rwork(11) = hu |
1326 | rwork(12) = h |
1327 | rwork(13) = tn |
1328 | iwork(11) = nst |
1329 | iwork(12) = nfe |
1330 | iwork(13) = nje |
1331 | iwork(14) = nqu |
1332 | iwork(15) = nq |
1333 | return |
1334 | c |
1335 | 430 ntrep = ntrep + 1 |
1336 | if (ntrep .lt. 5) return |
1337 | call xascwv( |
1338 | 1 60hlsode-- repeated calls with istate = 1 and tout = t (=r1) , |
1339 | 1 60, 301, 0, 0, 0, 0, 1, t, 0.0d0) |
1340 | go to 800 |
1341 | c----------------------------------------------------------------------- |
1342 | c block h. |
1343 | c the following block handles all unsuccessful returns other than |
1344 | c those for illegal input. first the error message routine is called. |
1345 | c if there was an error test or convergence test failure, imxer is set. |
1346 | c then y is loaded from yh, t is set to tn, and the illegal input |
1347 | c counter illin is set to 0. the optional outputs are loaded into |
1348 | c the work arrays before returning. |
1349 | c----------------------------------------------------------------------- |
1350 | c the maximum number of steps was taken before reaching tout. ---------- |
1351 | 500 call xascwv(50hlsode-- at current t (=r1), mxstep (=i1) steps , |
1352 | 1 50, 201, 0, 0, 0, 0, 0, 0.0d0, 0.0d0) |
1353 | call xascwv(50h taken on this call before reaching tout , |
1354 | 1 50, 201, 0, 1, mxstep, 0, 1, tn, 0.0d0) |
1355 | istate = -1 |
1356 | go to 580 |
1357 | c ewt(i) .le. 0.0 for some i (not at start of problem). ---------------- |
1358 | 510 ewti = rwork(lewt+i-1) |
1359 | call xascwv(50hlsode-- at t (=r1), ewt(i1) has become r2 .le. 0., |
1360 | 1 50, 202, 0, 1, i, 0, 2, tn, ewti) |
1361 | istate = -6 |
1362 | go to 580 |
1363 | c too much accuracy requested for machine precision. ------------------- |
1364 | 520 call xascwv(50hlsode-- at t (=r1), too much accuracy requested , |
1365 | 1 50, 203, 0, 0, 0, 0, 0, 0.0d0, 0.0d0) |
1366 | call xascwv(50h for precision of machine.. see tolsf (=r2) , |
1367 | 1 50, 203, 0, 0, 0, 0, 2, tn, tolsf) |
1368 | rwork(14) = tolsf |
1369 | istate = -2 |
1370 | go to 580 |
1371 | c kflag = -1. error test failed repeatedly or with abs(h) = hmin. ----- |
1372 | 530 call xascwv(50hlsode-- at t(=r1) and step size h(=r2), the error, |
1373 | 1 50, 204, 0, 0, 0, 0, 0, 0.0d0, 0.0d0) |
1374 | call xascwv(50h test failed repeatedly or with abs(h) = hmin, |
1375 | 1 50, 204, 0, 0, 0, 0, 2, tn, h) |
1376 | istate = -4 |
1377 | go to 560 |
1378 | c kflag = -2. convergence failed repeatedly or with abs(h) = hmin. ---- |
1379 | 540 call xascwv(50hlsode-- at t (=r1) and step size h (=r2), the , |
1380 | 1 50, 205, 0, 0, 0, 0, 0, 0.0d0, 0.0d0) |
1381 | call xascwv(50h corrector convergence failed repeatedly , |
1382 | 1 50, 205, 0, 0, 0, 0, 0, 0.0d0, 0.0d0) |
1383 | call xascwv(30h or with abs(h) = hmin , |
1384 | 1 30, 205, 0, 0, 0, 0, 2, tn, h) |
1385 | istate = -5 |
1386 | c compute imxer if relevant. ------------------------------------------- |
1387 | 560 big = 0.0d0 |
1388 | imxer = 1 |
1389 | do 570 i = 1,n |
1390 | size = dabs(rwork(i+lacor-1)*rwork(i+lewt-1)) |
1391 | if (big .ge. size) go to 570 |
1392 | big = size |
1393 | imxer = i |
1394 | 570 continue |
1395 | iwork(16) = imxer |
1396 | c set y vector, t, illin, and optional outputs. ------------------------ |
1397 | 580 do 590 i = 1,n |
1398 | 590 y(i) = rwork(i+lyh-1) |
1399 | t = tn |
1400 | illin = 0 |
1401 | rwork(11) = hu |
1402 | rwork(12) = h |
1403 | rwork(13) = tn |
1404 | iwork(11) = nst |
1405 | iwork(12) = nfe |
1406 | iwork(13) = nje |
1407 | iwork(14) = nqu |
1408 | iwork(15) = nq |
1409 | return |
1410 | c----------------------------------------------------------------------- |
1411 | c block i. |
1412 | c the following block handles all error returns due to illegal input |
1413 | c (istate = -3), as detected before calling the core integrator. |
1414 | c first the error message routine is called. then if there have been |
1415 | c 5 consecutive such returns just before this call to the solver, |
1416 | c the run is halted. |
1417 | c----------------------------------------------------------------------- |
1418 | 601 call xascwv(30hlsode-- istate (=i1) illegal , |
1419 | 1 30, 1, 0, 1, istate, 0, 0, 0.0d0, 0.0d0) |
1420 | go to 700 |
1421 | 602 call xascwv(30hlsode-- itask (=i1) illegal , |
1422 | 1 30, 2, 0, 1, itask, 0, 0, 0.0d0, 0.0d0) |
1423 | go to 700 |
1424 | 603 call xascwv(50hlsode-- istate .gt. 1 but lsode not initialized , |
1425 | 1 50, 3, 0, 0, 0, 0, 0, 0.0d0, 0.0d0) |
1426 | go to 700 |
1427 | 604 call xascwv(30hlsode-- neq (=i1) .lt. 1 , |
1428 | 1 30, 4, 0, 1, neq(1), 0, 0, 0.0d0, 0.0d0) |
1429 | go to 700 |
1430 | 605 call xascwv(50hlsode-- istate = 3 and neq increased (i1 to i2) , |
1431 | 1 50, 5, 0, 2, n, neq(1), 0, 0.0d0, 0.0d0) |
1432 | go to 700 |
1433 | 606 call xascwv(30hlsode-- itol (=i1) illegal , |
1434 | 1 30, 6, 0, 1, itol, 0, 0, 0.0d0, 0.0d0) |
1435 | go to 700 |
1436 | 607 call xascwv(30hlsode-- iopt (=i1) illegal , |
1437 | 1 30, 7, 0, 1, iopt, 0, 0, 0.0d0, 0.0d0) |
1438 | go to 700 |
1439 | 608 call xascwv(30hlsode-- mf (=i1) illegal , |
1440 | 1 30, 8, 0, 1, mf, 0, 0, 0.0d0, 0.0d0) |
1441 | go to 700 |
1442 | 609 call xascwv(50hlsode-- ml (=i1) illegal.. .lt.0 or .ge.neq (=i2), |
1443 | 1 50, 9, 0, 2, ml, neq(1), 0, 0.0d0, 0.0d0) |
1444 | go to 700 |
1445 | 610 call xascwv(50hlsode-- mu (=i1) illegal.. .lt.0 or .ge.neq (=i2), |
1446 | 1 50, 10, 0, 2, mu, neq(1), 0, 0.0d0, 0.0d0) |
1447 | go to 700 |
1448 | 611 call xascwv(30hlsode-- maxord (=i1) .lt. 0 , |
1449 | 1 30, 11, 0, 1, maxord, 0, 0, 0.0d0, 0.0d0) |
1450 | go to 700 |
1451 | 612 call xascwv(30hlsode-- mxstep (=i1) .lt. 0 , |
1452 | 1 30, 12, 0, 1, mxstep, 0, 0, 0.0d0, 0.0d0) |
1453 | go to 700 |
1454 | 613 call xascwv(30hlsode-- mxhnil (=i1) .lt. 0 , |
1455 | 1 30, 13, 0, 1, mxhnil, 0, 0, 0.0d0, 0.0d0) |
1456 | go to 700 |
1457 | 614 call xascwv(40hlsode-- tout (=r1) behind t (=r2) , |
1458 | 1 40, 14, 0, 0, 0, 0, 2, tout, t) |
1459 | call xascwv(50h integration direction is given by h0 (=r1) , |
1460 | 1 50, 14, 0, 0, 0, 0, 1, h0, 0.0d0) |
1461 | go to 700 |
1462 | 615 call xascwv(30hlsode-- hmax (=r1) .lt. 0.0 , |
1463 | 1 30, 15, 0, 0, 0, 0, 1, hmax, 0.0d0) |
1464 | go to 700 |
1465 | 616 call xascwv(30hlsode-- hmin (=r1) .lt. 0.0 , |
1466 | 1 30, 16, 0, 0, 0, 0, 1, hmin, 0.0d0) |
1467 | go to 700 |
1468 | 617 call xascwv( |
1469 | 1 60hlsode-- rwork length needed, lenrw (=i1), exceeds lrw (=i2), |
1470 | 1 60, 17, 0, 2, lenrw, lrw, 0, 0.0d0, 0.0d0) |
1471 | go to 700 |
1472 | 618 call xascwv( |
1473 | 1 60hlsode-- iwork length needed, leniw (=i1), exceeds liw (=i2), |
1474 | 1 60, 18, 0, 2, leniw, liw, 0, 0.0d0, 0.0d0) |
1475 | go to 700 |
1476 | 619 call xascwv(40hlsode-- rtol(i1) is r1 .lt. 0.0 , |
1477 | 1 40, 19, 0, 1, i, 0, 1, rtoli, 0.0d0) |
1478 | go to 700 |
1479 | 620 call xascwv(40hlsode-- atol(i1) is r1 .lt. 0.0 , |
1480 | 1 40, 20, 0, 1, i, 0, 1, atoli, 0.0d0) |
1481 | go to 700 |
1482 | 621 ewti = rwork(lewt+i-1) |
1483 | call xascwv(40hlsode-- ewt(i1) is r1 .le. 0.0 , |
1484 | 1 40, 21, 0, 1, i, 0, 1, ewti, 0.0d0) |
1485 | go to 700 |
1486 | 622 call xascwv( |
1487 | 1 60hlsode-- tout (=r1) too close to t(=r2) to start integration, |
1488 | 1 60, 22, 0, 0, 0, 0, 2, tout, t) |
1489 | go to 700 |
1490 | 623 call xascwv( |
1491 | 1 60hlsode-- itask = i1 and tout (=r1) behind tcur - hu (= r2) , |
1492 | 1 60, 23, 0, 1, itask, 0, 2, tout, tp) |
1493 | go to 700 |
1494 | 624 call xascwv( |
1495 | 1 60hlsode-- itask = 4 or 5 and tcrit (=r1) behind tcur (=r2) , |
1496 | 1 60, 24, 0, 0, 0, 0, 2, tcrit, tn) |
1497 | go to 700 |
1498 | 625 call xascwv( |
1499 | 1 60hlsode-- itask = 4 or 5 and tcrit (=r1) behind tout (=r2) , |
1500 | 1 60, 25, 0, 0, 0, 0, 2, tcrit, tout) |
1501 | go to 700 |
1502 | 626 call xascwv(50hlsode-- at start of problem, too much accuracy , |
1503 | 1 50, 26, 0, 0, 0, 0, 0, 0.0d0, 0.0d0) |
1504 | call xascwv( |
1505 | 1 60h requested for precision of machine.. see tolsf (=r1) , |
1506 | 1 60, 26, 0, 0, 0, 0, 1, tolsf, 0.0d0) |
1507 | rwork(14) = tolsf |
1508 | go to 700 |
1509 | 627 call xascwv(50hlsode-- trouble from intdy. itask = i1, tout = r1, |
1510 | 1 50, 27, 0, 1, itask, 0, 1, tout, 0.0d0) |
1511 | c |
1512 | 700 if (illin .eq. 5) go to 710 |
1513 | illin = illin + 1 |
1514 | istate = -3 |
1515 | return |
1516 | 710 call xascwv(50hlsode-- repeated occurrences of illegal input , |
1517 | 1 50, 302, 0, 0, 0, 0, 0, 0.0d0, 0.0d0) |
1518 | c |
1519 | 800 call xascwv(50hlsode-- run aborted.. apparent infinite loop , |
1520 | 1 50, 303, 2, 0, 0, 0, 0, 0.0d0, 0.0d0) |
1521 | return |
1522 | c----------------------- end of subroutine lsode ----------------------- |
1523 | end |
1524 | subroutine solsy (wm, iwm, x, tem) |
1525 | clll. optimize |
1526 | integer iwm |
1527 | integer iownd, iowns, |
1528 | 1 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, |
1529 | 2 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu |
1530 | integer i, meband, ml, mu |
1531 | double precision wm, x, tem |
1532 | double precision rowns, |
1533 | 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround |
1534 | double precision di, hl0, phl0, r |
1535 | dimension wm(1), iwm(1), x(1), tem(1) |
1536 | common /ls0001/ rowns(209), |
1537 | 2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround, |
1538 | 3 iownd(14), iowns(6), |
1539 | 4 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, |
1540 | 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu |
1541 | c----------------------------------------------------------------------- |
1542 | c this routine manages the solution of the linear system arising from |
1543 | c a chord iteration. it is called if miter .ne. 0. |
1544 | c if miter is 1 or 2, it calls dgesl to accomplish this. |
1545 | c if miter = 3 it updates the coefficient h*el0 in the diagonal |
1546 | c matrix, and then computes the solution. |
1547 | c if miter is 4 or 5, it calls dgbsl. |
1548 | c communication with solsy uses the following variables.. |
1549 | c wm = real work space containing the inverse diagonal matrix if |
1550 | c miter = 3 and the lu decomposition of the matrix otherwise. |
1551 | c storage of matrix elements starts at wm(3). |
1552 | c wm also contains the following matrix-related data.. |
1553 | c wm(1) = sqrt(uround) (not used here), |
1554 | c wm(2) = hl0, the previous value of h*el0, used if miter = 3. |
1555 | c iwm = integer work space containing pivot information, starting at |
1556 | c iwm(21), if miter is 1, 2, 4, or 5. iwm also contains band |
1557 | c parameters ml = iwm(1) and mu = iwm(2) if miter is 4 or 5. |
1558 | c x = the right-hand side vector on input, and the solution vector |
1559 | c on output, of length n. |
1560 | c tem = vector of work space of length n, not used in this version. |
1561 | c iersl = output flag (in common). iersl = 0 if no trouble occurred. |
1562 | c iersl = 1 if a singular matrix arose with miter = 3. |
1563 | c this routine also uses the common variables el0, h, miter, and n. |
1564 | c----------------------------------------------------------------------- |
1565 | iersl = 0 |
1566 | go to (100, 100, 300, 400, 400), miter |
1567 | 100 call dgesl (wm(3), n, n, iwm(21), x, 0) |
1568 | return |
1569 | c |
1570 | 300 phl0 = wm(2) |
1571 | hl0 = h*el0 |
1572 | wm(2) = hl0 |
1573 | if (hl0 .eq. phl0) go to 330 |
1574 | r = hl0/phl0 |
1575 | do 320 i = 1,n |
1576 | di = 1.0d0 - r*(1.0d0 - 1.0d0/wm(i+2)) |
1577 | if (dabs(di) .eq. 0.0d0) go to 390 |
1578 | 320 wm(i+2) = 1.0d0/di |
1579 | 330 do 340 i = 1,n |
1580 | 340 x(i) = wm(i+2)*x(i) |
1581 | return |
1582 | 390 iersl = 1 |
1583 | return |
1584 | c |
1585 | 400 ml = iwm(1) |
1586 | mu = iwm(2) |
1587 | meband = 2*ml + mu + 1 |
1588 | call dgbsl (wm(3), meband, n, ml, mu, iwm(21), x, 0) |
1589 | return |
1590 | c----------------------- end of subroutine solsy ----------------------- |
1591 | end |
1592 | subroutine prepj (neq, y, yh, nyh, ewt, ftem, savf, wm, iwm, |
1593 | 1 f, jac) |
1594 | clll. optimize |
1595 | cascend changes |
1596 | external f, jac |
1597 | integer neq, nyh, iwm |
1598 | integer iownd, iowns, |
1599 | 1 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, |
1600 | 2 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu |
1601 | integer i, i1, i2, ier, ii, j, j1, jj, lenp, |
1602 | 1 mba, mband, meb1, meband, ml, ml3, mu, np1 |
1603 | double precision y, yh, ewt, ftem, savf, wm |
1604 | double precision rowns, |
1605 | 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround |
1606 | double precision con, di, fac, hl0, r, r0, srur, yi, yj, yjj, |
1607 | 1 vnorm |
1608 | Ckaa double precision time1, time2 |
1609 | dimension neq(1), y(1), yh(nyh,1), ewt(1), ftem(1), savf(1), |
1610 | 1 wm(1), iwm(1) |
1611 | common /ls0001/ rowns(209), |
1612 | 2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround, |
1613 | 3 iownd(14), iowns(6), |
1614 | 4 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, |
1615 | 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu |
1616 | c----------------------------------------------------------------------- |
1617 | c prepj is called by stode to compute and process the matrix |
1618 | c p = i - h*el(1)*j , where j is an approximation to the jacobian. |
1619 | c here j is computed by the user-supplied routine jac if |
1620 | c miter = 1 or 4, or by finite differencing if miter = 2, 3, or 5. |
1621 | c if miter = 3, a diagonal approximation to j is used. |
1622 | c j is stored in wm and replaced by p. if miter .ne. 3, p is then |
1623 | c subjected to lu decomposition in preparation for later solution |
1624 | c of linear systems with p as coefficient matrix. this is done |
1625 | c by dgefa if miter = 1 or 2, and by dgbfa if miter = 4 or 5. |
1626 | c |
1627 | c in addition to variables described previously, communication |
1628 | c with prepj uses the following.. |
1629 | c y = array containing predicted values on entry. |
1630 | c ftem = work array of length n (acor in stode). |
1631 | c savf = array containing f evaluated at predicted y. |
1632 | c wm = real work space for matrices. on output it contains the |
1633 | c inverse diagonal matrix if miter = 3 and the lu decomposition |
1634 | c of p if miter is 1, 2 , 4, or 5. |
1635 | c storage of matrix elements starts at wm(3). |
1636 | c wm also contains the following matrix-related data.. |
1637 | c wm(1) = sqrt(uround), used in numerical jacobian increments. |
1638 | c wm(2) = h*el0, saved for later use if miter = 3. |
1639 | c iwm = integer work space containing pivot information, starting at |
1640 | c iwm(21), if miter is 1, 2, 4, or 5. iwm also contains band |
1641 | c parameters ml = iwm(1) and mu = iwm(2) if miter is 4 or 5. |
1642 | c el0 = el(1) (input). |
1643 | c ierpj = output error flag, = 0 if no trouble, .gt. 0 if |
1644 | c p matrix found to be singular. |
1645 | c jcur = output flag = 1 to indicate that the jacobian matrix |
1646 | c (or approximation) is now current. |
1647 | c this routine also uses the common variables el0, h, tn, uround, |
1648 | c miter, n, nfe, and nje. |
1649 | c----------------------------------------------------------------------- |
1650 | nje = nje + 1 |
1651 | ierpj = 0 |
1652 | jcur = 1 |
1653 | hl0 = h*el0 |
1654 | go to (100, 200, 300, 400, 500), miter |
1655 | c if miter = 1, call jac and multiply by scalar. ----------------------- |
1656 | 100 lenp = n*n |
1657 | do 110 i = 1,lenp |
1658 | 110 wm(i+2) = 0.0d0 |
1659 | call jac (neq, tn, y, 0, 0, wm(3), n) |
1660 | con = -hl0 |
1661 | do 120 i = 1,lenp |
1662 | 120 wm(i+2) = wm(i+2)*con |
1663 | go to 240 |
1664 | c if miter = 2, make n calls to f to approximate j. -------------------- |
1665 | 200 fac = vnorm (n, savf, ewt) |
1666 | r0 = 1000.0d0*dabs(h)*uround*dfloat(n)*fac |
1667 | if (r0 .eq. 0.0d0) r0 = 1.0d0 |
1668 | srur = wm(1) |
1669 | j1 = 2 |
1670 | do 230 j = 1,n |
1671 | yj = y(j) |
1672 | r = dmax1(srur*dabs(yj),r0/ewt(j)) |
1673 | y(j) = y(j) + r |
1674 | fac = -hl0/r |
1675 | call f (neq, tn, y, ftem) |
1676 | do 220 i = 1,n |
1677 | 220 wm(i+j1) = (ftem(i) - savf(i))*fac |
1678 | y(j) = yj |
1679 | j1 = j1 + n |
1680 | 230 continue |
1681 | nfe = nfe + n |
1682 | c add identity matrix. ------------------------------------------------- |
1683 | 240 j = 3 |
1684 | np1 = n + 1 |
1685 | do 250 i = 1,n |
1686 | wm(j) = wm(j) + 1.0d0 |
1687 | 250 j = j + np1 |
1688 | c do lu decomposition on p. -------------------------------------------- |
1689 | Ckaa call aftime(time1) |
1690 | call dgefa (wm(3), n, n, iwm(21), ier) |
1691 | Ckaa call aftime(time2) |
1692 | Ckaa time2 = time2 - time1 |
1693 | Ckaa write (6,*) "Time for decomposition" |
1694 | Ckaa write (6,700) time2 |
1695 | if (ier .ne. 0) ierpj = 1 |
1696 | return |
1697 | c if miter = 3, construct a diagonal approximation to j and p. --------- |
1698 | 300 wm(2) = hl0 |
1699 | r = el0*0.1d0 |
1700 | do 310 i = 1,n |
1701 | 310 y(i) = y(i) + r*(h*savf(i) - yh(i,2)) |
1702 | call f (neq, tn, y, wm(3)) |
1703 | nfe = nfe + 1 |
1704 | do 320 i = 1,n |
1705 | r0 = h*savf(i) - yh(i,2) |
1706 | di = 0.1d0*r0 - h*(wm(i+2) - savf(i)) |
1707 | wm(i+2) = 1.0d0 |
1708 | if (dabs(r0) .lt. uround/ewt(i)) go to 320 |
1709 | if (dabs(di) .eq. 0.0d0) go to 330 |
1710 | wm(i+2) = 0.1d0*r0/di |
1711 | 320 continue |
1712 | return |
1713 | 330 ierpj = 1 |
1714 | return |
1715 | c if miter = 4, call jac and multiply by scalar. ----------------------- |
1716 | 400 ml = iwm(1) |
1717 | mu = iwm(2) |
1718 | ml3 = ml + 3 |
1719 | mband = ml + mu + 1 |
1720 | meband = mband + ml |
1721 | lenp = meband*n |
1722 | do 410 i = 1,lenp |
1723 | 410 wm(i+2) = 0.0d0 |
1724 | call jac (neq, tn, y, ml, mu, wm(ml3), meband) |
1725 | con = -hl0 |
1726 | do 420 i = 1,lenp |
1727 | 420 wm(i+2) = wm(i+2)*con |
1728 | go to 570 |
1729 | c if miter = 5, make mband calls to f to approximate j. ---------------- |
1730 | 500 ml = iwm(1) |
1731 | mu = iwm(2) |
1732 | mband = ml + mu + 1 |
1733 | mba = min0(mband,n) |
1734 | meband = mband + ml |
1735 | meb1 = meband - 1 |
1736 | srur = wm(1) |
1737 | fac = vnorm (n, savf, ewt) |
1738 | r0 = 1000.0d0*dabs(h)*uround*dfloat(n)*fac |
1739 | if (r0 .eq. 0.0d0) r0 = 1.0d0 |
1740 | do 560 j = 1,mba |
1741 | do 530 i = j,n,mband |
1742 | yi = y(i) |
1743 | r = dmax1(srur*dabs(yi),r0/ewt(i)) |
1744 | 530 y(i) = y(i) + r |
1745 | call f (neq, tn, y, ftem) |
1746 | do 550 jj = j,n,mband |
1747 | y(jj) = yh(jj,1) |
1748 | yjj = y(jj) |
1749 | r = dmax1(srur*dabs(yjj),r0/ewt(jj)) |
1750 | fac = -hl0/r |
1751 | i1 = max0(jj-mu,1) |
1752 | i2 = min0(jj+ml,n) |
1753 | ii = jj*meb1 - ml + 2 |
1754 | do 540 i = i1,i2 |
1755 | 540 wm(ii+i) = (ftem(i) - savf(i))*fac |
1756 | 550 continue |
1757 | 560 continue |
1758 | nfe = nfe + mba |
1759 | c add identity matrix. ------------------------------------------------- |
1760 | 570 ii = mband + 2 |
1761 | do 580 i = 1,n |
1762 | wm(ii) = wm(ii) + 1.0d0 |
1763 | 580 ii = ii + meband |
1764 | c do lu decomposition of p. -------------------------------------------- |
1765 | call dgbfa (wm(3), meband, n, ml, mu, iwm(21), ier) |
1766 | if (ier .ne. 0) ierpj = 1 |
1767 | return |
1768 | 700 format(d21.13) |
1769 | c----------------------- end of subroutine prepj ----------------------- |
1770 | end |
1771 | subroutine stode (neq, y, yh, nyh, yh1, ewt, savf, acor, |
1772 | 1 wm, iwm, f, jac, pjac, slvs) |
1773 | clll. optimize |
1774 | external f, jac, pjac, slvs |
1775 | integer neq, nyh, iwm |
1776 | integer iownd, ialth, ipup, lmax, meo, nqnyh, nslp, |
1777 | 1 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, |
1778 | 2 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu |
1779 | integer i, i1, iredo, iret, j, jb, m, ncf, newq |
1780 | double precision y, yh, yh1, ewt, savf, acor, wm |
1781 | double precision conit, crate, el, elco, hold, rmax, tesco, |
1782 | 2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround |
1783 | double precision dcon, ddn, del, delp, dsm, dup, exdn, exsm, exup, |
1784 | 1 r, rh, rhdn, rhsm, rhup, told, vnorm |
1785 | dimension neq(1), y(1), yh(nyh,1), yh1(1), ewt(1), savf(1), |
1786 | 1 acor(1), wm(1), iwm(1) |
1787 | common /ls0001/ conit, crate, el(13), elco(13,12), |
1788 | 1 hold, rmax, tesco(3,12), |
1789 | 2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround, iownd(14), |
1790 | 3 ialth, ipup, lmax, meo, nqnyh, nslp, |
1791 | 4 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, |
1792 | 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu |
1793 | c----------------------------------------------------------------------- |
1794 | c stode performs one step of the integration of an initial value |
1795 | c problem for a system of ordinary differential equations. |
1796 | c note.. stode is independent of the value of the iteration method |
1797 | c indicator miter, when this is .ne. 0, and hence is independent |
1798 | c of the type of chord method used, or the jacobian structure. |
1799 | c communication with stode is done with the following variables.. |
1800 | c |
1801 | c neq = integer array containing problem size in neq(1), and |
1802 | c passed as the neq argument in all calls to f and jac. |
1803 | c y = an array of length .ge. n used as the y argument in |
1804 | c all calls to f and jac. |
1805 | c yh = an nyh by lmax array containing the dependent variables |
1806 | c and their approximate scaled derivatives, where |
1807 | c lmax = maxord + 1. yh(i,j+1) contains the approximate |
1808 | c j-th derivative of y(i), scaled by h**j/factorial(j) |
1809 | c (j = 0,1,...,nq). on entry for the first step, the first |
1810 | c two columns of yh must be set from the initial values. |
1811 | c nyh = a constant integer .ge. n, the first dimension of yh. |
1812 | c yh1 = a one-dimensional array occupying the same space as yh. |
1813 | c ewt = an array of length n containing multiplicative weights |
1814 | c for local error measurements. local errors in y(i) are |
1815 | c compared to 1.0/ewt(i) in various error tests. |
1816 | c savf = an array of working storage, of length n. |
1817 | c also used for input of yh(*,maxord+2) when jstart = -1 |
1818 | c and maxord .lt. the current order nq. |
1819 | c acor = a work array of length n, used for the accumulated |
1820 | c corrections. on a successful return, acor(i) contains |
1821 | c the estimated one-step local error in y(i). |
1822 | c wm,iwm = real and integer work arrays associated with matrix |
1823 | c operations in chord iteration (miter .ne. 0). |
1824 | c pjac = name of routine to evaluate and preprocess jacobian matrix |
1825 | c and p = i - h*el0*jac, if a chord method is being used. |
1826 | c slvs = name of routine to solve linear system in chord iteration. |
1827 | c ccmax = maximum relative change in h*el0 before pjac is called. |
1828 | c h = the step size to be attempted on the next step. |
1829 | c h is altered by the error control algorithm during the |
1830 | c problem. h can be either positive or negative, but its |
1831 | c sign must remain constant throughout the problem. |
1832 | c hmin = the minimum absolute value of the step size h to be used. |
1833 | c hmxi = inverse of the maximum absolute value of h to be used. |
1834 | c hmxi = 0.0 is allowed and corresponds to an infinite hmax. |
1835 | c hmin and hmxi may be changed at any time, but will not |
1836 | c take effect until the next change of h is considered. |
1837 | c tn = the independent variable. tn is updated on each step taken. |
1838 | c jstart = an integer used for input only, with the following |
1839 | c values and meanings.. |
1840 | c 0 perform the first step. |
1841 | c .gt.0 take a new step continuing from the last. |
1842 | c -1 take the next step with a new value of h, maxord, |
1843 | c n, meth, miter, and/or matrix parameters. |
1844 | c -2 take the next step with a new value of h, |
1845 | c but with other inputs unchanged. |
1846 | c on return, jstart is set to 1 to facilitate continuation. |
1847 | c kflag = a completion code with the following meanings.. |
1848 | c 0 the step was succesful. |
1849 | c -1 the requested error could not be achieved. |
1850 | c -2 corrector convergence could not be achieved. |
1851 | c -3 fatal error in pjac or slvs. |
1852 | c a return with kflag = -1 or -2 means either |
1853 | c abs(h) = hmin or 10 consecutive failures occurred. |
1854 | c on a return with kflag negative, the values of tn and |
1855 | c the yh array are as of the beginning of the last |
1856 | c step, and h is the last step size attempted. |
1857 | c maxord = the maximum order of integration method to be allowed. |
1858 | c maxcor = the maximum number of corrector iterations allowed. |
1859 | c msbp = maximum number of steps between pjac calls (miter .gt. 0). |
1860 | c mxncf = maximum number of convergence failures allowed. |
1861 | c meth/miter = the method flags. see description in driver. |
1862 | c n = the number of first-order differential equations. |
1863 | c----------------------------------------------------------------------- |
1864 | kflag = 0 |
1865 | told = tn |
1866 | ncf = 0 |
1867 | ierpj = 0 |
1868 | iersl = 0 |
1869 | jcur = 0 |
1870 | icf = 0 |
1871 | delp = 0.0d0 |
1872 | if (jstart .gt. 0) go to 200 |
1873 | if (jstart .eq. -1) go to 100 |
1874 | if (jstart .eq. -2) go to 160 |
1875 | c----------------------------------------------------------------------- |
1876 | c on the first call, the order is set to 1, and other variables are |
1877 | c initialized. rmax is the maximum ratio by which h can be increased |
1878 | c in a single step. it is initially 1.e4 to compensate for the small |
1879 | c initial h, but then is normally equal to 10. if a failure |
1880 | c occurs (in corrector convergence or error test), rmax is set at 2 |
1881 | c for the next increase. |
1882 | c----------------------------------------------------------------------- |
1883 | lmax = maxord + 1 |
1884 | nq = 1 |
1885 | l = 2 |
1886 | ialth = 2 |
1887 | rmax = 10000.0d0 |
1888 | rc = 0.0d0 |
1889 | el0 = 1.0d0 |
1890 | crate = 0.7d0 |
1891 | hold = h |
1892 | meo = meth |
1893 | nslp = 0 |
1894 | ipup = miter |
1895 | iret = 3 |
1896 | go to 140 |
1897 | c----------------------------------------------------------------------- |
1898 | c the following block handles preliminaries needed when jstart = -1. |
1899 | c ipup is set to miter to force a matrix update. |
1900 | c if an order increase is about to be considered (ialth = 1), |
1901 | c ialth is reset to 2 to postpone consideration one more step. |
1902 | c if the caller has changed meth, cfode is called to reset |
1903 | c the coefficients of the method. |
1904 | c if the caller has changed maxord to a value less than the current |
1905 | c order nq, nq is reduced to maxord, and a new h chosen accordingly. |
1906 | c if h is to be changed, yh must be rescaled. |
1907 | c if h or meth is being changed, ialth is reset to l = nq + 1 |
1908 | c to prevent further changes in h for that many steps. |
1909 | c----------------------------------------------------------------------- |
1910 | 100 ipup = miter |
1911 | lmax = maxord + 1 |
1912 | if (ialth .eq. 1) ialth = 2 |
1913 | if (meth .eq. meo) go to 110 |
1914 | call cfode (meth, elco, tesco) |
1915 | meo = meth |
1916 | if (nq .gt. maxord) go to 120 |
1917 | ialth = l |
1918 | iret = 1 |
1919 | go to 150 |
1920 | 110 if (nq .le. maxord) go to 160 |
1921 | 120 nq = maxord |
1922 | l = lmax |
1923 | do 125 i = 1,l |
1924 | 125 el(i) = elco(i,nq) |
1925 | nqnyh = nq*nyh |
1926 | rc = rc*el(1)/el0 |
1927 | el0 = el(1) |
1928 | conit = 0.5d0/dfloat(nq+2) |
1929 | ddn = vnorm (n, savf, ewt)/tesco(1,l) |
1930 | exdn = 1.0d0/dfloat(l) |
1931 | rhdn = 1.0d0/(1.3d0*ddn**exdn + 0.0000013d0) |
1932 | rh = dmin1(rhdn,1.0d0) |
1933 | iredo = 3 |
1934 | if (h .eq. hold) go to 170 |
1935 | rh = dmin1(rh,dabs(h/hold)) |
1936 | h = hold |
1937 | go to 175 |
1938 | c----------------------------------------------------------------------- |
1939 | c cfode is called to get all the integration coefficients for the |
1940 | c current meth. then the el vector and related constants are reset |
1941 | c whenever the order nq is changed, or at the start of the problem. |
1942 | c----------------------------------------------------------------------- |
1943 | 140 call cfode (meth, elco, tesco) |
1944 | 150 do 155 i = 1,l |
1945 | 155 el(i) = elco(i,n |