/[ascend]/trunk/lsod/lsode.f
ViewVC logotype

Contents of /trunk/lsod/lsode.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (show annotations) (download)
Fri Oct 29 20:54:12 2004 UTC (15 years, 7 months ago) by aw0a
File size: 118049 byte(s)
Setting up web subdirectory in repository
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,nq)
1946 nqnyh = nq*nyh