subroutine lsodar (f, neq, y, t, tout, itol, rtol, atol, itask,
1 istate, iopt, rwork, lrw, iwork, liw, jac, jt,
2 g, ng, jroot)
external f, jac, g
integer neq, itol, itask, istate, iopt, lrw, iwork, liw, jt,
1 ng, jroot
double precision y, t, tout, rtol, atol, rwork
dimension neq(1), y(1), rtol(1), atol(1), rwork(lrw), iwork(liw),
1 jroot(ng)
c-----------------------------------------------------------------------
c this is the may 7, 1982 version of
c lsodar.. livermore solver for ordinary differential equations, with
c automatic method switching for stiff and nonstiff problems,
c and with root-finding.
c
c this version is in double precision.
c
c lsodar solves the initial value problem for stiff or nonstiff
c systems of first order ode-s,
c dy/dt = f(t,y) , or, in component form,
c dy(i)/dt = f(i) = f(i,t,y(1),y(2),...,y(neq)) (i = 1,...,neq).
c at the same time, it locates the roots of any of a set of functions
c g(i) = g(i,t,y(1),...,y(neq)) (i = 1,...,ng).
c
c this a variant version of the lsode package. it differs from lsode
c in two ways..
c (a) it switches automatically between stiff and nonstiff methods.
c this means that the user does not have to determine whether the
c problem is stiff or not, and the solver will automatically choose the
c appropriate method. it always starts with the nonstiff method.
c (b) it finds the root of at least one of a set of constraint
c functions g(i) of the independent and dependent variables.
c it finds only those roots for which some g(i), as a function
c of t, changes sign in the interval of integration.
c it then returns the solution at the root, if that occurs
c sooner than the specified stop condition, and otherwise returns
c the solution according the specified stop condition.
c
c authors..
c linda r. petzold
c applied mathematics division 8331
c sandia national laboratories
c livermore, ca 94550
c and
c alan c. hindmarsh,
c mathematics and statistics division, l-316
c lawrence livermore national laboratory
c livermore, ca 94550.
c
c references..
c 1. alan c. hindmarsh, lsode and lsodi, two new initial value
c ordinary differential equation solvers,
c acm-signum newsletter, vol. 15, no. 4 (1980), pp. 10-11.
c 2. linda r. petzold, automatic selection of methods for solving
c stiff and nonstiff systems of ordinary differential equations,
c siam j. sci. stat. comput. 4 (1983), pp. 136-148.
c 3. kathie l. hiebert and lawrence f. shampine, implicitly defined
c output points for solutions of ode-s, sandia report sand80-0180,
c february, 1980.
c-----------------------------------------------------------------------
c summary of usage.
c
c communication between the user and the lsodar package, for normal
c situations, is summarized here. this summary describes only a subset
c of the full set of options available. see the full description for
c details, including alternative treatment of the jacobian matrix,
c optional inputs and outputs, nonstandard options, and
c instructions for special situations. see also the example
c problem (with program and output) following this summary.
c
c a. first provide a subroutine of the form..
c subroutine f (neq, t, y, ydot)
c dimension y(neq), ydot(neq)
c which supplies the vector function f by loading ydot(i) with f(i).
c
c b. provide a subroutine of the form..
c subroutine g (neq, t, y, ng, gout)
c dimension y(neq), gout(ng)
c which supplies the vector function g by loading gout(i) with
c g(i), the i-th constraint function whose root is sought.
c
c c. write a main program which calls subroutine lsodar once for
c each point at which answers are desired. this should also provide
c for possible use of logical unit 6 for output of error messages by
c lsodar. on the first call to lsodar, supply arguments as follows..
c f = name of subroutine for right-hand side vector f.
c this name must be declared external in calling program.
c neq = number of first order ode-s.
c y = array of initial values, of length neq.
c t = the initial value of the independent variable.
c tout = first point where output is desired (.ne. t).
c itol = 1 or 2 according as atol (below) is a scalar or array.
c rtol = relative tolerance parameter (scalar).
c atol = absolute tolerance parameter (scalar or array).
c the estimated local error in y(i) will be controlled so as
c to be less than
c ewt(i) = rtol*abs(y(i)) + atol if itol = 1, or
c ewt(i) = rtol*abs(y(i)) + atol(i) if itol = 2.
c thus the local error test passes if, in each component,
c either the absolute error is less than atol (or atol(i)),
c or the relative error is less than rtol.
c use rtol = 0.0 for pure absolute error control, and
c use atol = 0.0 (or atol(i) = 0.0) for pure relative error
c control. caution.. actual (global) errors may exceed these
c local tolerances, so choose them conservatively.
c itask = 1 for normal computation of output values of y at t = tout.
c istate = integer flag (input and output). set istate = 1.
c iopt = 0 to indicate no optional inputs used.
c rwork = real work array of length at least..
c 22 + neq * max(16, neq + 9) + 3*ng.
c see also paragraph f below.
c lrw = declared length of rwork (in user-s dimension).
c iwork = integer work array of length at least 20 + neq.
c liw = declared length of iwork (in user-s dimension).
c jac = name of subroutine for jacobian matrix.
c use a dummy name. see also paragraph f below.
c jt = jacobian type indicator. set jt = 2.
c see also paragraph f below.
c g = name of subroutine for constraint functions, whose
c roots are desired during the integration.
c this name must be declared external in calling program.
c ng = number of constraint functions g(i). if there are none,
c set ng = 0, and pass a dummy name for g.
c jroot = integer array of length ng for output of root information.
c see next paragraph.
c note that the main program must declare arrays y, rwork, iwork,
c jroot, and possibly atol.
c
c d. the output from the first call (or any call) is..
c y = array of computed values of y(t) vector.
c t = corresponding value of independent variable. this is
c tout if istate = 2, or the root location if istate = 3,
c or the farthest point reached if lsodar was unsuccessful.
c istate = 2 or 3 if lsodar was successful, negative otherwise.
c 2 means no root was found, and tout was reached as desired.
c 3 means a root was found prior to reaching tout.
c -1 means excess work done on this call (perhaps wrong jt).
c -2 means excess accuracy requested (tolerances too small).
c -3 means illegal input detected (see printed message).
c -4 means repeated error test failures (check all inputs).
c -5 means repeated convergence failures (perhaps bad jacobian
c supplied or wrong choice of jt or tolerances).
c -6 means error weight became zero during problem. (solution
c component i vanished, and atol or atol(i) = 0.)
c -7 means work space insufficient to finish (see messages).
c jroot = array showing roots found if istate = 3 on return.
c jroot(i) = 1 if g(i) has a root at t, or 0 otherwise.
c
c e. to continue the integration after a successful return, proceed
c as follows..
c (a) if istate = 2 on return, reset tout and call lsodar again.
c (b) if istate = 3 on return, reset istate to 2 and call lsodar again.
c in either case, no other parameters need be reset.
c
c f. note.. if and when lsodar regards the problem as stiff, and
c switches methods accordingly, it must make use of the neq by neq
c jacobian matrix, j = df/dy. for the sake of simplicity, the
c inputs to lsodar recommended in paragraph c above cause lsodar to
c treat j as a full matrix, and to approximate it internally by
c difference quotients. alternatively, j can be treated as a band
c matrix (with great potential reduction in the size of the rwork
c array). also, in either the full or banded case, the user can supply
c j in closed form, with a routine whose name is passed as the jac
c argument. these alternatives are described in the paragraphs on
c rwork, jac, and jt in the full description of the call sequence below.
c
c-----------------------------------------------------------------------
c example problem.
c
c the following is a simple example problem, with the coding
c needed for its solution by lsodar. the problem is from chemical
c kinetics, and consists of the following three rate equations..
c dy1/dt = -.04*y1 + 1.e4*y2*y3
c dy2/dt = .04*y1 - 1.e4*y2*y3 - 3.e7*y2**2
c dy3/dt = 3.e7*y2**2
c on the interval from t = 0.0 to t = 4.e10, with initial conditions
c y1 = 1.0, y2 = y3 = 0. the problem is stiff.
c in addition, we want to find the values of t, y1, y2, and y3 at which
c (1) y1 reaches the value 1.e-4, and
c (2) y3 reaches the value 1.e-2.
c
c the following coding solves this problem with lsodar,
c printing results at t = .4, 4., ..., 4.e10, and at the computed
c roots. it uses itol = 2 and atol much smaller for y2 than y1 or y3
c because y2 has much smaller values.
c at the end of the run, statistical quantities of interest are
c printed (see optional outputs in the full description below).
c
c external fex, gex
c double precision atol, rtol, rwork, t, tout, y
c dimension y(3), atol(3), rwork(76), iwork(23), jroot(2)
c neq = 3
c y(1) = 1.0d0
c y(2) = 0.0d0
c y(3) = 0.0d0
c t = 0.0d0
c tout = 0.4d0
c itol = 2
c rtol = 1.0d-4
c atol(1) = 1.0d-6
c atol(2) = 1.0d-10
c atol(3) = 1.0d-6
c itask = 1
c istate = 1
c iopt = 0
c lrw = 76
c liw = 23
c jt = 2
c ng = 2
c do 40 iout = 1,12
c 10 call lsodar(fex,neq,y,t,tout,itol,rtol,atol,itask,istate,
c 1 iopt,rwork,lrw,iwork,liw,jdum,jt,gex,ng,jroot)
c write(6,20)t,y(1),y(2),y(3)
c 20 format(7h at t =,e12.4,6h y =,3e14.6)
c if (istate .lt. 0) go to 80
c if (istate .eq. 2) go to 40
c write(6,30)jroot(1),jroot(2)
c 30 format(5x,35h the above line is a root, jroot =,2i5)
c istate = 2
c go to 10
c 40 tout = tout*10.0d0
c write(6,60)iwork(11),iwork(12),iwork(13),iwork(10),
c 1 iwork(19),rwork(15)
c 60 format(/12h no. steps =,i4,11h no. f-s =,i4,11h no. j-s =,i4,
c 1 11h no. g-s =,i4/
c 2 19h method last used =,i2,25h last switch was at t =,e12.4)
c stop
c 80 write(6,90)istate
c 90 format(///22h error halt.. istate =,i3)
c stop
c end
c
c subroutine fex (neq, t, y, ydot)
c double precision t, y, ydot
c dimension y(3), ydot(3)
c ydot(1) = -0.04d0*y(1) + 1.0d4*y(2)*y(3)
c ydot(3) = 3.0d7*y(2)*y(2)
c ydot(2) = -ydot(1) - ydot(3)
c return
c end
c
c subroutine gex (neq, t, y, ng, gout)
c double precision t, y, gout
c dimension y(3), gout(2)
c gout(1) = y(1) - 1.0d-4
c gout(2) = y(3) - 1.0d-2
c return
c end
c
c the output of this program (on a cdc-7600 in single precision)
c is as follows..
c
c at t = 2.6400e-01 y = 9.899653e-01 3.470563e-05 1.000000e-02
c the above line is a root, jroot = 0 1
c at t = 4.0000e-01 y = 9.851712e-01 3.386380e-05 1.479493e-02
c at t = 4.0000e+00 y = 9.055333e-01 2.240655e-05 9.444430e-02
c at t = 4.0000e+01 y = 7.158403e-01 9.186334e-06 2.841505e-01
c at t = 4.0000e+02 y = 4.505250e-01 3.222964e-06 5.494717e-01
c at t = 4.0000e+03 y = 1.831975e-01 8.941774e-07 8.168016e-01
c at t = 4.0000e+04 y = 3.898730e-02 1.621940e-07 9.610125e-01
c at t = 4.0000e+05 y = 4.936363e-03 1.984221e-08 9.950636e-01
c at t = 4.0000e+06 y = 5.161831e-04 2.065786e-09 9.994838e-01
c at t = 2.0745e+07 y = 1.000000e-04 4.000395e-10 9.999000e-01
c the above line is a root, jroot = 1 0
c at t = 4.0000e+07 y = 5.179817e-05 2.072032e-10 9.999482e-01
c at t = 4.0000e+08 y = 5.283401e-06 2.113371e-11 9.999947e-01
c at t = 4.0000e+09 y = 4.659031e-07 1.863613e-12 9.999995e-01
c at t = 4.0000e+10 y = 1.404280e-08 5.617126e-14 1.000000e+00
c
c no. steps = 361 no. f-s = 693 no. j-s = 64 no. g-s = 390
c method last used = 2 last switch was at t = 6.0092e-03
c-----------------------------------------------------------------------
c full description of user interface to lsodar.
c
c the user interface to lsodar consists of the following parts.
c
c i. the call sequence to subroutine lsodar, which is a driver
c routine for the solver. this includes descriptions of both
c the call sequence arguments and of user-supplied routines.
c following these descriptions is a description of
c optional inputs available through the call sequence, and then
c a description of optional outputs (in the work arrays).
c
c ii. descriptions of other routines in the lsodar package that may be
c (optionally) called by the user. these provide the ability to
c alter error message handling, save and restore the internal
c common, and obtain specified derivatives of the solution y(t).
c
c iii. descriptions of common blocks to be declared in overlay
c or similar environments, or to be saved when doing an interrupt
c of the problem and continued solution later.
c
c iv. description of a subroutine in the lsodar package,
c which the user may replace with his own version, if desired.
c this relates to the measurement of errors.
c
c-----------------------------------------------------------------------
c part i. call sequence.
c
c the call sequence parameters used for input only are
c f, neq, tout, itol, rtol, atol, itask, iopt, lrw, liw, jac,
c jt, g, and ng,
c that used only for output is jroot,
c and those used for both input and output are
c y, t, istate.
c the work arrays rwork and iwork are also used for conditional and
c optional inputs and optional outputs. (the term output here refers
c to the return from subroutine lsodar to the user-s calling program.)
c
c the legality of input parameters will be thoroughly checked on the
c initial call for the problem, but not checked thereafter unless a
c change in input parameters is flagged by istate = 3 on input.
c
c the descriptions of the call arguments are as follows.
c
c f = the name of the user-supplied subroutine defining the
c ode system. the system must be put in the first-order
c form dy/dt = f(t,y), where f is a vector-valued function
c of the scalar t and the vector y. subroutine f is to
c compute the function f. it is to have the form
c subroutine f (neq, t, y, ydot)
c dimension y(1), ydot(1)
c where neq, t, and y are input, and the array ydot = f(t,y)
c is output. y and ydot are arrays of length neq.
c (in the dimension statement above, 1 is a dummy
c dimension.. it can be replaced by any value.)
c subroutine f should not alter y(1),...,y(neq).
c f must be declared external in the calling program.
c
c subroutine f may access user-defined quantities in
c neq(2),... and y(neq(1)+1),... if neq is an array
c (dimensioned in f) and y has length exceeding neq(1).
c see the descriptions of neq and y below.
c
c neq = the size of the ode system (number of first order
c ordinary differential equations). used only for input.
c neq may be decreased, but not increased, during the problem.
c if neq is decreased (with istate = 3 on input), the
c remaining components of y should be left undisturbed, if
c these are to be accessed in f and/or jac.
c
c normally, neq is a scalar, and it is generally referred to
c as a scalar in this user interface description. however,
c neq may be an array, with neq(1) set to the system size.
c (the lsodar package accesses only neq(1).) in either case,
c this parameter is passed as the neq argument in all calls
c to f, jac, and g. hence, if it is an array, locations
c neq(2),... may be used to store other integer data and pass
c it to f, jac, and g. each such subroutine must include
c neq in a dimension statement in that case.
c
c y = a real array for the vector of dependent variables, of
c length neq or more. used for both input and output on the
c first call (istate = 1), and only for output on other calls.
c on the first call, y must contain the vector of initial
c values. on output, y contains the computed solution vector,
c evaluated at t. if desired, the y array may be used
c for other purposes between calls to the solver.
c
c this array is passed as the y argument in all calls to f,
c jac, and g. hence its length may exceed neq, and locations
c y(neq+1),... may be used to store other real data and
c pass it to f, jac, and g. (the lsodar package accesses only
c y(1),...,y(neq).)
c
c t = the independent variable. on input, t is used only on the
c first call, as the initial point of the integration.
c on output, after each call, t is the value at which a
c computed solution y is evaluated (usually the same as tout).
c if a root was found, t is the computed location of the
c root reached first, on output.
c on an error return, t is the farthest point reached.
c
c tout = the next value of t at which a computed solution is desired.
c used only for input.
c
c when starting the problem (istate = 1), tout may be equal
c to t for one call, then should .ne. t for the next call.
c for the initial t, an input value of tout .ne. t is used
c in order to determine the direction of the integration
c (i.e. the algebraic sign of the step sizes) and the rough
c scale of the problem. integration in either direction
c (forward or backward in t) is permitted.
c
c if itask = 2 or 5 (one-step modes), tout is ignored after
c the first call (i.e. the first call with tout .ne. t).
c otherwise, tout is required on every call.
c
c if itask = 1, 3, or 4, the values of tout need not be
c monotone, but a value of tout which backs up is limited
c to the current internal t interval, whose endpoints are
c tcur - hu and tcur (see optional outputs, below, for
c tcur and hu).
c
c itol = an indicator for the type of error control. see
c description below under atol. used only for input.
c
c rtol = a relative error tolerance parameter, either a scalar or
c an array of length neq. see description below under atol.
c input only.
c
c atol = an absolute error tolerance parameter, either a scalar or
c an array of length neq. input only.
c
c the input parameters itol, rtol, and atol determine
c the error control performed by the solver. the solver will
c control the vector e = (e(i)) of estimated local errors
c in y, according to an inequality of the form
c max-norm of ( e(i)/ewt(i) ) .le. 1,
c where ewt = (ewt(i)) is a vector of positive error weights.
c the values of rtol and atol should all be non-negative.
c the following table gives the types (scalar/array) of
c rtol and atol, and the corresponding form of ewt(i).
c
c itol rtol atol ewt(i)
c 1 scalar scalar rtol*abs(y(i)) + atol
c 2 scalar array rtol*abs(y(i)) + atol(i)
c 3 array scalar rtol(i)*abs(y(i)) + atol
c 4 array array rtol(i)*abs(y(i)) + atol(i)
c
c when either of these parameters is a scalar, it need not
c be dimensioned in the user-s calling program.
c
c if none of the above choices (with itol, rtol, and atol
c fixed throughout the problem) is suitable, more general
c error controls can be obtained by substituting a
c user-supplied routine for the setting of ewt.
c see part iv below.
c
c if global errors are to be estimated by making a repeated
c run on the same problem with smaller tolerances, then all
c components of rtol and atol (i.e. of ewt) should be scaled
c down uniformly.
c
c itask = an index specifying the task to be performed.
c input only. itask has the following values and meanings.
c 1 means normal computation of output values of y(t) at
c t = tout (by overshooting and interpolating).
c 2 means take one step only and return.
c 3 means stop at the first internal mesh point at or
c beyond t = tout and return.
c 4 means normal computation of output values of y(t) at
c t = tout but without overshooting t = tcrit.
c tcrit must be input as rwork(1). tcrit may be equal to
c or beyond tout, but not behind it in the direction of
c integration. this option is useful if the problem
c has a singularity at or beyond t = tcrit.
c 5 means take one step, without passing tcrit, and return.
c tcrit must be input as rwork(1).
c
c note.. if itask = 4 or 5 and the solver reaches tcrit
c (within roundoff), it will return t = tcrit (exactly) to
c indicate this (unless itask = 4 and tout comes before tcrit,
c in which case answers at t = tout are returned first).
c
c istate = an index used for input and output to specify the
c the state of the calculation.
c
c on input, the values of istate are as follows.
c 1 means this is the first call for the problem
c (initializations will be done). see note below.
c 2 means this is not the first call, and the calculation
c is to continue normally, with no change in any input
c parameters except possibly tout and itask.
c (if itol, rtol, and/or atol are changed between calls
c with istate = 2, the new values will be used but not
c tested for legality.)
c 3 means this is not the first call, and the
c calculation is to continue normally, but with
c a change in input parameters other than
c tout and itask. changes are allowed in
c neq, itol, rtol, atol, iopt, lrw, liw, jt, ml, mu,
c and any optional inputs except h0, mxordn, and mxords.
c (see iwork description for ml and mu.)
c in addition, immediately following a return with
c istate = 3 (root found), ng and g may be changed.
c (but changing ng from 0 to .gt. 0 is not allowed.)
c note.. a preliminary call with tout = t is not counted
c as a first call here, as no initialization or checking of
c input is done. (such a call is sometimes useful for the
c purpose of outputting the initial conditions.)
c thus the first call for which tout .ne. t requires
c istate = 1 on input.
c
c on output, istate has the following values and meanings.
c 1 means nothing was done, as tout was equal to t with
c istate = 1 on input. (however, an internal counter was
c set to detect and prevent repeated calls of this type.)
c 2 means the integration was performed successfully, and
c no roots were found.
c 3 means the integration was successful, and one or more
c roots were found before satisfying the stop condition
c specified by itask. see jroot.
c -1 means an excessive amount of work (more than mxstep
c steps) was done on this call, before completing the
c requested task, but the integration was otherwise
c successful as far as t. (mxstep is an optional input
c and is normally 500.) to continue, the user may
c simply reset istate to a value .gt. 1 and call again
c (the excess work step counter will be reset to 0).
c in addition, the user may increase mxstep to avoid
c this error return (see below on optional inputs).
c -2 means too much accuracy was requested for the precision
c of the machine being used. this was detected before
c completing the requested task, but the integration
c was successful as far as t. to continue, the tolerance
c parameters must be reset, and istate must be set
c to 3. the optional output tolsf may be used for this
c purpose. (note.. if this condition is detected before
c taking any steps, then an illegal input return
c (istate = -3) occurs instead.)
c -3 means illegal input was detected, before taking any
c integration steps. see written message for details.
c note.. if the solver detects an infinite loop of calls
c to the solver with illegal input, it will cause
c the run to stop.
c -4 means there were repeated error test failures on
c one attempted step, before completing the requested
c task, but the integration was successful as far as t.
c the problem may have a singularity, or the input
c may be inappropriate.
c -5 means there were repeated convergence test failures on
c one attempted step, before completing the requested
c task, but the integration was successful as far as t.
c this may be caused by an inaccurate jacobian matrix,
c if one is being used.
c -6 means ewt(i) became zero for some i during the
c integration. pure relative error control (atol(i)=0.0)
c was requested on a variable which has now vanished.
c the integration was successful as far as t.
c -7 means the length of rwork and/or iwork was too small to
c proceed, but the integration was successful as far as t.
c this happens when lsodar chooses to switch methods
c but lrw and/or liw is too small for the new method.
c
c note.. since the normal output value of istate is 2,
c it does not need to be reset for normal continuation.
c also, since a negative input value of istate will be
c regarded as illegal, a negative output value requires the
c user to change it, and possibly other inputs, before
c calling the solver again.
c
c iopt = an integer flag to specify whether or not any optional
c inputs are being used on this call. input only.
c the optional inputs are listed separately below.
c iopt = 0 means no optional inputs are being used.
c default values will be used in all cases.
c iopt = 1 means one or more optional inputs are being used.
c
c rwork = a real array (double precision) for work space, and (in the
c first 20 words) for conditional and optional inputs and
c optional outputs.
c as lsodar switches automatically between stiff and nonstiff
c methods, the required length of rwork can change during the
c problem. thus the rwork array passed to lsodar can either
c have a static (fixed) length large enough for both methods,
c or have a dynamic (changing) length altered by the calling
c program in response to output from lsodar.
c
c --- fixed length case ---
c if the rwork length is to be fixed, it should be at least
c max (lrn, lrs),
c where lrn and lrs are the rwork lengths required when the
c current method is nonstiff or stiff, respectively.
c
c the separate rwork length requirements lrn and lrs are
c as follows..
c if neq is constant and the maximum method orders have
c their default values, then
c lrn = 20 + 16*neq + 3*ng,
c lrs = 22 + 9*neq + neq**2 + 3*ng (jt = 1 or 2),
c lrs = 22 + 10*neq + (2*ml+mu)*neq + 3*ng (jt = 4 or 5).
c under any other conditions, lrn and lrs are given by..
c lrn = 20 + nyh*(mxordn+1) + 3*neq + 3*ng,
c lrs = 20 + nyh*(mxords+1) + 3*neq + lmat + 3*ng,
c where
c nyh = the initial value of neq,
c mxordn = 12, unless a smaller value is given as an
c optional input,
c mxords = 5, unless a smaller value is given as an
c optional input,
c lmat = length of matrix work space..
c lmat = neq**2 + 2 if jt = 1 or 2,
c lmat = (2*ml + mu + 1)*neq + 2 if jt = 4 or 5.
c
c --- dynamic length case ---
c if the length of rwork is to be dynamic, then it should
c be at least lrn or lrs, as defined above, depending on the
c current method. initially, it must be at least lrn (since
c lsodar starts with the nonstiff method). on any return
c from lsodar, the optional output mcur indicates the current
c method. if mcur differs from the value it had on the
c previous return, or if there has only been one call to
c lsodar and mcur is now 2, then lsodar has switched
c methods during the last call, and the length of rwork
c should be reset (to lrn if mcur = 1, or to lrs if
c mcur = 2). (an increase in the rwork length is required
c if lsodar returned istate = -7, but not otherwise.)
c after resetting the length, call lsodar with istate = 3
c to signal that change.
c
c lrw = the length of the array rwork, as declared by the user.
c (this will be checked by the solver.)
c
c iwork = an integer array for work space.
c as lsodar switches automatically between stiff and nonstiff
c methods, the required length of iwork can change during
c problem, between
c lis = 20 + neq and lin = 20,
c respectively. thus the iwork array passed to lsodar can
c either have a fixed length of at least 20 + neq, or have a
c dynamic length of at least lin or lis, depending on the
c current method. the comments on dynamic length under
c rwork above apply here. initially, this length need
c only be at least lin = 20.
c
c the first few words of iwork are used for conditional and
c optional inputs and optional outputs.
c
c the following 2 words in iwork are conditional inputs..
c iwork(1) = ml these are the lower and upper
c iwork(2) = mu half-bandwidths, respectively, of the
c banded jacobian, excluding the main diagonal.
c the band is defined by the matrix locations
c (i,j) with i-ml .le. j .le. i+mu. ml and mu
c must satisfy 0 .le. ml,mu .le. neq-1.
c these are required if jt is 4 or 5, and
c ignored otherwise. ml and mu may in fact be
c the band parameters for a matrix to which
c df/dy is only approximately equal.
c
c liw = the length of the array iwork, as declared by the user.
c (this will be checked by the solver.)
c
c note.. the base addresses of the work arrays must not be
c altered between calls to lsodar for the same problem.
c the contents of the work arrays must not be altered
c between calls, except possibly for the conditional and
c optional inputs, and except for the last 3*neq words of rwork.
c the latter space is used for internal scratch space, and so is
c available for use by the user outside lsodar between calls, if
c desired (but not for use by f, jac, or g).
c
c jac = the name of the user-supplied routine to compute the
c jacobian matrix, df/dy, if jt = 1 or 4. the jac routine
c is optional, but if the problem is expected to be stiff much
c of the time, you are encouraged to supply jac, for the sake
c of efficiency. (alternatively, set jt = 2 or 5 to have
c lsodar compute df/dy internally by difference quotients.)
c if and when lsodar uses df/dy, if treats this neq by neq
c matrix either as full (jt = 1 or 2), or as banded (jt =
c 4 or 5) with half-bandwidths ml and mu (discussed under
c iwork above). in either case, if jt = 1 or 4, the jac
c routine must compute df/dy as a function of the scalar t
c and the vector y. it is to have the form
c subroutine jac (neq, t, y, ml, mu, pd, nrowpd)
c dimension y(1), pd(nrowpd,1)
c where neq, t, y, ml, mu, and nrowpd are input and the array
c pd is to be loaded with partial derivatives (elements of
c the jacobian matrix) on output. pd must be given a first
c dimension of nrowpd. t and y have the same meaning as in
c subroutine f. (in the dimension statement above, 1 is a
c dummy dimension.. it can be replaced by any value.)
c in the full matrix case (jt = 1), ml and mu are
c ignored, and the jacobian is to be loaded into pd in
c columnwise manner, with df(i)/dy(j) loaded into pd(i,j).
c in the band matrix case (jt = 4), the elements
c within the band are to be loaded into pd in columnwise
c manner, with diagonal lines of df/dy loaded into the rows
c of pd. thus df(i)/dy(j) is to be loaded into pd(i-j+mu+1,j).
c ml and mu are the half-bandwidth parameters (see iwork).
c the locations in pd in the two triangular areas which
c correspond to nonexistent matrix elements can be ignored
c or loaded arbitrarily, as they are overwritten by lsodar.
c jac need not provide df/dy exactly. a crude
c approximation (possibly with a smaller bandwidth) will do.
c in either case, pd is preset to zero by the solver,
c so that only the nonzero elements need be loaded by jac.
c each call to jac is preceded by a call to f with the same
c arguments neq, t, and y. thus to gain some efficiency,
c intermediate quantities shared by both calculations may be
c saved in a user common block by f and not recomputed by jac,
c if desired. also, jac may alter the y array, if desired.
c jac must be declared external in the calling program.
c subroutine jac may access user-defined quantities in
c neq(2),... and y(neq(1)+1),... if neq is an array
c (dimensioned in jac) and y has length exceeding neq(1).
c see the descriptions of neq and y above.
c
c jt = jacobian type indicator. used only for input.
c jt specifies how the jacobian matrix df/dy will be
c treated, if and when lsodar requires this matrix.
c jt has the following values and meanings..
c 1 means a user-supplied full (neq by neq) jacobian.
c 2 means an internally generated (difference quotient) full
c jacobian (using neq extra calls to f per df/dy value).
c 4 means a user-supplied banded jacobian.
c 5 means an internally generated banded jacobian (using
c ml+mu+1 extra calls to f per df/dy evaluation).
c if jt = 1 or 4, the user must supply a subroutine jac
c (the name is arbitrary) as described above under jac.
c if jt = 2 or 5, a dummy argument can be used.
c
c g = the name of subroutine for constraint functions, whose
c roots are desired during the integration. it is to have
c the form
c subroutine g (neq, t, y, ng, gout)
c dimension y(neq), gout(ng)
c where neq, t, y, and ng are input, and the array gout
c is output. neq, t, and y have the same meaning as in
c the f routine, and gout is an array of length ng.
c for i = 1,...,ng, this routine is to load into gout(i)
c the value at (t,y) of the i-th constraint function g(i).
c lsodar will find roots of the g(i) of odd multiplicity
c (i.e. sign changes) as they occur during the integration.
c g must be declared external in the calling program.
c
c caution.. because of numerical errors in the functions
c g(i) due to roundoff and integration error, lsodar may
c return false roots, or return the same root at two or more
c nearly equal values of t. if such false roots are
c suspected, the user should consider smaller error tolerances
c and/or higher precision in the evaluation of the g(i).
c
c if a root of some g(i) defines the end of the problem,
c the input to lsodar should nevertheless allow integration
c to a point slightly past that root, so that lsodar can
c locate the root by interpolation.
c
c subroutine g may access user-defined quantities in
c neq(2),... and y(neq(1)+1),... if neq is an array
c (dimensioned in g) and y has length exceeding neq(1).
c see the descriptions of neq and y above.
c
c ng = number of constraint functions g(i). if there are none,
c set ng = 0, and pass a dummy name for g.
c
c jroot = integer array of length ng. used only for output.
c on a return with istate = 3 (one or more roots found),
c jroot(i) = 1 if g(i) has a root at t, or jroot(i) = 0 if not.
c-----------------------------------------------------------------------
c optional inputs.
c
c the following is a list of the optional inputs provided for in the
c call sequence. (see also part ii.) for each such input variable,
c this table lists its name as used in this documentation, its
c location in the call sequence, its meaning, and the default value.
c the use of any of these inputs requires iopt = 1, and in that
c case all of these inputs are examined. a value of zero for any
c of these optional inputs will cause the default value to be used.
c thus to use a subset of the optional inputs, simply preload
c locations 5 to 10 in rwork and iwork to 0.0 and 0 respectively, and
c then set those of interest to nonzero values.
c
c name location meaning and default value
c
c h0 rwork(5) the step size to be attempted on the first step.
c the default value is determined by the solver.
c
c hmax rwork(6) the maximum absolute step size allowed.
c the default value is infinite.
c
c hmin rwork(7) the minimum absolute step size allowed.
c the default value is 0. (this lower bound is not
c enforced on the final step before reaching tcrit
c when itask = 4 or 5.)
c
c ixpr iwork(5) flag to generate extra printing at method switches.
c ixpr = 0 means no extra printing (the default).
c ixpr = 1 means print data on each switch.
c t, h, and nst will be printed on the same logical
c unit as used for error messages.
c
c mxstep iwork(6) maximum number of (internally defined) steps
c allowed during one call to the solver.
c the default value is 500.
c
c mxhnil iwork(7) maximum number of messages printed (per problem)
c warning that t + h = t on a step (h = step size).
c this must be positive to result in a non-default
c value. the default value is 10.
c
c mxordn iwork(8) the maximum order to be allowed for the nonstiff
c (adams) method. the default value is 12.
c if mxordn exceeds the default value, it will
c be reduced to the default value.
c mxordn is held constant during the problem.
c
c mxords iwork(9) the maximum order to be allowed for the stiff
c (bdf) method. the default value is 5.
c if mxords exceeds the default value, it will
c be reduced to the default value.
c mxords is held constant during the problem.
c-----------------------------------------------------------------------
c optional outputs.
c
c as optional additional output from lsodar, the variables listed
c below are quantities related to the performance of lsodar
c which are available to the user. these are communicated by way of
c the work arrays, but also have internal mnemonic names as shown.
c except where stated otherwise, all of these outputs are defined
c on any successful return from lsodar, and on any return with
c istate = -1, -2, -4, -5, or -6. on an illegal input return
c (istate = -3), they will be unchanged from their existing values
c (if any), except possibly for tolsf, lenrw, and leniw.
c on any error return, outputs relevant to the error will be defined,
c as noted below.
c
c name location meaning
c
c hu rwork(11) the step size in t last used (successfully).
c
c hcur rwork(12) the step size to be attempted on the next step.
c
c tcur rwork(13) the current value of the independent variable
c which the solver has actually reached, i.e. the
c current internal mesh point in t. on output, tcur
c will always be at least as far as the argument
c t, but may be farther (if interpolation was done).
c
c tolsf rwork(14) a tolerance scale factor, greater than 1.0,
c computed when a request for too much accuracy was
c detected (istate = -3 if detected at the start of
c the problem, istate = -2 otherwise). if itol is
c left unaltered but rtol and atol are uniformly
c scaled up by a factor of tolsf for the next call,
c then the solver is deemed likely to succeed.
c (the user may also ignore tolsf and alter the
c tolerance parameters in any other way appropriate.)
c
c tsw rwork(15) the value of t at the time of the last method
c switch, if any.
c
c nge iwork(10) the number of g evaluations for the problem so far.
c
c nst iwork(11) the number of steps taken for the problem so far.
c
c nfe iwork(12) the number of f evaluations for the problem so far.
c
c nje iwork(13) the number of jacobian evaluations (and of matrix
c lu decompositions) for the problem so far.
c
c nqu iwork(14) the method order last used (successfully).
c
c nqcur iwork(15) the order to be attempted on the next step.
c
c imxer iwork(16) the index of the component of largest magnitude in
c the weighted local error vector ( e(i)/ewt(i) ),
c on an error return with istate = -4 or -5.
c
c lenrw iwork(17) the length of rwork actually required, assuming
c that the length of rwork is to be fixed for the
c rest of the problem, and that switching may occur.
c this is defined on normal returns and on an illegal
c input return for insufficient storage.
c
c leniw iwork(18) the length of iwork actually required, assuming
c that the length of iwork is to be fixed for the
c rest of the problem, and that switching may occur.
c this is defined on normal returns and on an illegal
c input return for insufficient storage.
c
c mused iwork(19) the method indicator for the last successful step..
c 1 means adams (nonstiff), 2 means bdf (stiff).
c
c mcur iwork(20) the current method indicator..
c 1 means adams (nonstiff), 2 means bdf (stiff).
c this is the method to be attempted
c on the next step. thus it differs from mused
c only if a method switch has just been made.
c
c the following two arrays are segments of the rwork array which
c may also be of interest to the user as optional outputs.
c for each array, the table below gives its internal name,
c its base address in rwork, and its description.
c
c name base address description
c
c yh 21 + 3*ng the nordsieck history array, of size nyh by
c (nqcur + 1), where nyh is the initial value
c of neq. for j = 0,1,...,nqcur, column j+1
c of yh contains hcur**j/factorial(j) times
c the j-th derivative of the interpolating
c polynomial currently representing the solution,
c evaluated at t = tcur.
c
c acor lacor array of size neq used for the accumulated
c (from common corrections on each step, scaled on output
c as noted) to represent the estimated local error in y
c on the last step. this is the vector e in
c the description of the error control. it is
c defined only on a successful return from
c lsodar. the base address lacor is obtained by
c including in the user-s program the
c following 3 lines..
c double precision rls
c common /ls0001/ rls(219), ils(39)
c lacor = ils(5)
c
c-----------------------------------------------------------------------
c part ii. other routines callable.
c
c the following are optional calls which the user may make to
c gain additional capabilities in conjunction with lsodar.
c (the routines xsetun and xsetf are designed to conform to the
c slatec error handling package.)
c
c form of call function
c call xsetun(lun) set the logical unit number, lun, for
c output of messages from lsodar, if
c the default is not desired.
c the default value of lun is 6.
c
c call xsetf(mflag) set a flag to control the printing of
c messages by lsodar.
c mflag = 0 means do not print. (danger..
c this risks losing valuable information.)
c mflag = 1 means print (the default).
c
c either of the above calls may be made at
c any time and will take effect immediately.
c
c call svcar (rsav, isav) store in rsav and isav the contents
c of the internal common blocks used by
c lsodar (see part iii below).
c rsav must be a real array of length 246
c or more, and isav must be an integer
c array of length 59 or more.
c
c call rscar (rsav, isav) restore, from rsav and isav, the contents
c of the internal common blocks used by
c lsodar. presumes a prior call to svcar
c with the same arguments.
c
c svcar and rscar are useful if
c interrupting a run and restarting
c later, or alternating between two or
c more problems solved with lsodar.
c
c call intdy(,,,,,) provide derivatives of y, of various
c (see below) orders, at a specified point t, if
c desired. it may be called only after
c a successful return from lsodar.
c
c the detailed instructions for using intdy are as follows.
c the form of the call is..
c
c call intdy (t, k, rwork(lyh), nyh, dky, iflag)
c
c the input parameters are..
c
c t = value of independent variable where answers are desired
c (normally the same as the t last returned by lsodar).
c for valid results, t must lie between tcur - hu and tcur.
c (see optional outputs for tcur and hu.)
c k = integer order of the derivative desired. k must satisfy
c 0 .le. k .le. nqcur, where nqcur is the current order
c (see optional outputs). the capability corresponding
c to k = 0, i.e. computing y(t), is already provided
c by lsodar directly. since nqcur .ge. 1, the first
c derivative dy/dt is always available with intdy.
c lyh = 21 + 3*ng = base address in rwork of the history array yh.
c nyh = column length of yh, equal to the initial value of neq.
c
c the output parameters are..
c
c dky = a real array of length neq containing the computed value
c of the k-th derivative of y(t).
c iflag = integer flag, returned as 0 if k and t were legal,
c -1 if k was illegal, and -2 if t was illegal.
c on an error return, a message is also written.
c-----------------------------------------------------------------------
c part iii. common blocks.
c
c if lsodar is to be used in an overlay situation, the user
c must declare, in the primary overlay, the variables in..
c (1) the call sequence to lsodar,
c (2) the four internal common blocks
c /ls0001/ of length 258 (219 double precision words
c followed by 39 integer words),
c /lsa001/ of length 31 (22 double precision words
c followed by 9 integer words),
c /lsr001/ of length 14 (5 double precision words
c followed by 9 integer words),
c /eh0001/ of length 2 (integer words).
c
c if lsodar is used on a system in which the contents of internal
c common blocks are not preserved between calls, the user should
c declare the above common blocks in his main program to insure
c that their contents are preserved.
c
c if the solution of a given problem by lsodar is to be interrupted
c and then later continued, such as when restarting an interrupted run
c or alternating between two or more problems, the user should save,
c following the return from the last lsodar call prior to the
c interruption, the contents of the call sequence variables and the
c internal common blocks, and later restore these values before the
c next lsodar call for that problem. to save and restore the common
c blocks, use subroutines svcar and rscar (see part ii above).
c
c-----------------------------------------------------------------------
c part iv. optionally replaceable solver routines.
c
c below is a description of a routine in the lsodar package which
c relates to the measurement of errors, and can be
c replaced by a user-supplied version, if desired. however, since such
c a replacement may have a major impact on performance, it should be
c done only when absolutely necessary, and only with great caution.
c (note.. the means by which the package version of a routine is
c superseded by the user-s version may be system-dependent.)
c
c (a) ewset.
c the following subroutine is called just before each internal
c integration step, and sets the array of error weights, ewt, as
c described under itol/rtol/atol above..
c subroutine ewset (neq, itol, rtol, atol, ycur, ewt)
c where neq, itol, rtol, and atol are as in the lsodar call sequence,
c ycur contains the current dependent variable vector, and
c ewt is the array of weights set by ewset.
c
c if the user supplies this subroutine, it must return in ewt(i)
c (i = 1,...,neq) a positive quantity suitable for comparing errors
c in y(i) to. the ewt array returned by ewset is passed to the
c vmnorm routine, and also used by lsodar in the computation
c of the optional output imxer, and the increments for difference
c quotient jacobians.
c
c in the user-supplied version of ewset, it may be desirable to use
c the current values of derivatives of y. derivatives up to order nq
c are available from the history array yh, described above under
c optional outputs. in ewset, yh is identical to the ycur array,
c extended to nq + 1 columns with a column length of nyh and scale
c factors of h**j/factorial(j). on the first call for the problem,
c given by nst = 0, nq is 1 and h is temporarily set to 1.0.
c the quantities nq, nyh, h, and nst can be obtained by including
c in ewset the statements..
c double precision h, rls
c common /ls0001/ rls(219),ils(39)
c nq = ils(35)
c nyh = ils(14)
c nst = ils(36)
c h = rls(213)
c thus, for example, the current value of dy/dt can be obtained as
c ycur(nyh+i)/h (i=1,...,neq) (and the division by h is
c unnecessary when nst = 0).
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c other routines in the lsodar package.
c
c in addition to subroutine lsodar, the lsodar package includes the
c following subroutines and function routines..
c rchek does preliminary checking for roots, and serves as an
c interface between subroutine lsodar and subroutine roots.
c roots finds the leftmost root of a set of functions.
c intdy computes an interpolated value of the y vector at t = tout.
c stoda is the core integrator, which does one step of the
c integration and the associated error control.
c cfode sets all method coefficients and test constants.
c prja computes and preprocesses the jacobian matrix j = df/dy
c and the newton iteration matrix p = i - h*l0*j.
c solsy manages solution of linear system in chord iteration.
c ewset sets the error weight vector ewt before each step.
c vmnorm computes the weighted max-norm of a vector.
c fnorm computes the norm of a full matrix consistent with the
c weighted max-norm on vectors.
c bnorm computes the norm of a band matrix consistent with the
c weighted max-norm on vectors.
c svcar and rscar are user-callable routines to save and restore,
c respectively, the contents of the internal common blocks.
c dgefa and dgesl are routines from linpack for solving full
c systems of linear algebraic equations.
c dgbfa and dgbsl are routines from linpack for solving banded
c linear systems.
c daxpy, dscal, idamax, ddot, and dcopy are basic linear algebra
c modules (blas) used by the above linpack routines.
c d1mach computes the unit roundoff in a machine-independent manner.
c xerrwx, xsetun, and xsetf handle the printing of all error
c messages and warnings. xerrwx is machine-dependent.
c note.. vmnorm, fnorm, bnorm, idamax, ddot, and d1mach are function
c routines. all the others are subroutines.
c
c the intrinsic and external routines used by lsodar are..
c dabs, dmax1, dmin1, dfloat, max0, min0, mod, dsign, dsqrt, and write.
c
c a block data subprogram is also included with the package,
c for loading some of the variables in internal common.
c
c-----------------------------------------------------------------------
c the following card is for optimized compilation on lll compilers.
clll. optimize
c-----------------------------------------------------------------------
external prja, solsy
integer illin, init, lyh, lewt, lacor, lsavf, lwm, liwm,
1 mxstep, mxhnil, nhnil, ntrep, nslast, nyh, iowns
integer icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
1 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
integer insufr, insufi, ixpr, iowns2, jtyp, mused, mxordn, mxords
integer lg0, lg1, lgx, iownr3, irfnd, itaskc, ngc, nge
integer i, i1, i2, iflag, imxer, kgo, lf0,
1 leniw, lenrw, lenwm, ml, mord, mu, mxhnl0, mxstp0
integer len1, len1c, len1n, len1s, len2, leniwc,
1 lenrwc, lenrwn, lenrws
integer irfp, irt, lenyh, lyhnew
double precision tret, rowns,
1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
double precision tsw, rowns2, pdnorm
double precision rownr3, t0, tlast, toutc
double precision atoli, ayi, big, ewti, h0, hmax, hmx, rh, rtoli,
1 tcrit, tdist, tnext, tol, tolsf, tp, size, sum, w0,
2 d1mach, vmnorm
dimension mord(2)
logical ihit
c-----------------------------------------------------------------------
c the following three internal common blocks contain
c (a) variables which are local to any subroutine but whose values must
c be preserved between calls to the routine (own variables), and
c (b) variables which are communicated between subroutines.
c the structure of each block is as follows.. all real variables are
c listed first, followed by all integers. within each type, the
c variables are grouped with those local to subroutine lsodar first,
c then those local to subroutine roots or subroutine stoda
c (no other routines have own variables), and finally those used
c for communication. the block ls0001 is declared in subroutines
c lsodar, intdy, stoda, prja, and solsy. the block lsa001 is declared
c in subroutines lsodar, stoda, and prja. the block lsr001 is declared
c in subroutines lsodar, rchek, and roots. groups of variables are
c replaced by dummy arrays in the common declarations in routines
c where those variables are not used.
c-----------------------------------------------------------------------
common /ls0001/ tret, rowns(209),
1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
2 illin, init, lyh, lewt, lacor, lsavf, lwm, liwm,
3 mxstep, mxhnil, nhnil, ntrep, nslast, nyh, iowns(6),
4 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
common /lsa001/ tsw, rowns2(20), pdnorm,
1 insufr, insufi, ixpr, iowns2(2), jtyp, mused, mxordn, mxords
common /lsr001/ rownr3(2), t0, tlast, toutc,
1 lg0, lg1, lgx, iownr3(2), irfnd, itaskc, ngc, nge
c
data mord(1),mord(2)/12,5/, mxstp0/500/, mxhnl0/10/
c-----------------------------------------------------------------------
c block a.
c this code block is executed on every call.
c it tests istate and itask for legality and branches appropriately.
c if istate .gt. 1 but the flag init shows that initialization has
c not yet been done, an error return occurs.
c if istate = 1 and tout = t, jump to block g and return immediately.
c-----------------------------------------------------------------------
if (istate .lt. 1 .or. istate .gt. 3) go to 601
if (itask .lt. 1 .or. itask .gt. 5) go to 602
itaskc = itask
if (istate .eq. 1) go to 10
if (init .eq. 0) go to 603
if (istate .eq. 2) go to 200
go to 20
10 init = 0
if (tout .eq. t) go to 430
20 ntrep = 0
c-----------------------------------------------------------------------
c block b.
c the next code block is executed for the initial call (istate = 1),
c or for a continuation call with parameter changes (istate = 3).
c it contains checking of all inputs and various initializations.
c
c first check legality of the non-optional inputs neq, itol, iopt,
c jt, ml, mu, and ng.
c-----------------------------------------------------------------------
if (neq(1) .le. 0) go to 604
if (istate .eq. 1) go to 25
if (neq(1) .gt. n) go to 605
25 n = neq(1)
if (itol .lt. 1 .or. itol .gt. 4) go to 606
if (iopt .lt. 0 .or. iopt .gt. 1) go to 607
if (jt .eq. 3 .or. jt .lt. 1 .or. jt .gt. 5) go to 608
jtyp = jt
if (jt .le. 2) go to 30
ml = iwork(1)
mu = iwork(2)
if (ml .lt. 0 .or. ml .ge. n) go to 609
if (mu .lt. 0 .or. mu .ge. n) go to 610
30 continue
if (ng .lt. 0) go to 630
if (istate .eq. 1) go to 35
if (irfnd .eq. 0 .and. ng .ne. ngc) go to 631
35 ngc = ng
c next process and check the optional inputs. --------------------------
if (iopt .eq. 1) go to 40
ixpr = 0
mxstep = mxstp0
mxhnil = mxhnl0
hmxi = 0.0d0
hmin = 0.0d0
if (istate .ne. 1) go to 60
h0 = 0.0d0
mxordn = mord(1)
mxords = mord(2)
go to 60
40 ixpr = iwork(5)
if (ixpr .lt. 0 .or. ixpr .gt. 1) go to 611
mxstep = iwork(6)
if (mxstep .lt. 0) go to 612
if (mxstep .eq. 0) mxstep = mxstp0
mxhnil = iwork(7)
if (mxhnil .lt. 0) go to 613
if (mxhnil .eq. 0) mxhnil = mxhnl0
if (istate .ne. 1) go to 50
h0 = rwork(5)
mxordn = iwork(8)
if (mxordn .lt. 0) go to 628
if (mxordn .eq. 0) mxordn = 100
mxordn = min0(mxordn,mord(1))
mxords = iwork(9)
if (mxords .lt. 0) go to 629
if (mxords .eq. 0) mxords = 100
mxords = min0(mxords,mord(2))
if ((tout - t)*h0 .lt. 0.0d0) go to 614
50 hmax = rwork(6)
if (hmax .lt. 0.0d0) go to 615
hmxi = 0.0d0
if (hmax .gt. 0.0d0) hmxi = 1.0d0/hmax
hmin = rwork(7)
if (hmin .lt. 0.0d0) go to 616
c-----------------------------------------------------------------------
c set work array pointers and check lengths lrw and liw.
c if istate = 1, meth is initialized to 1 here to facilitate the
c checking of work space lengths.
c pointers to segments of rwork and iwork are named by prefixing l to
c the name of the segment. e.g., the segment yh starts at rwork(lyh).
c segments of rwork (in order) are denoted g0, g1, gx, yh, wm,
c ewt, savf, acor.
c if the lengths provided are insufficient for the current method,
c an error return occurs. this is treated as illegal input on the
c first call, but as a problem interruption with istate = -7 on a
c continuation call. if the lengths are sufficient for the current
c method but not for both methods, a warning message is sent.
c-----------------------------------------------------------------------
60 if (istate .eq. 1) meth = 1
if (istate .eq. 1) nyh = n
lg0 = 21
lg1 = lg0 + ng
lgx = lg1 + ng
lyhnew = lgx + ng
if (istate .eq. 1) lyh = lyhnew
if (lyhnew .eq. lyh) go to 62
c if istate = 3 and ng was changed, shift yh to its new location. ------
lenyh = l*nyh
if (lrw .lt. lyhnew-1+lenyh) go to 62
i1 = 1
if (lyhnew .gt. lyh) i1 = -1
call dcopy (lenyh, rwork(lyh), i1, rwork(lyhnew), i1)
lyh = lyhnew
62 continue
len1n = lyhnew - 1 + (mxordn + 1)*nyh
len1s = lyhnew - 1 + (mxords + 1)*nyh
lwm = len1s + 1
if (jt .le. 2) lenwm = n*n + 2
if (jt .ge. 4) lenwm = (2*ml + mu + 1)*n + 2
len1s = len1s + lenwm
len1c = len1n
if (meth .eq. 2) len1c = len1s
len1 = max0(len1n,len1s)
len2 = 3*n
lenrw = len1 + len2
lenrwn = len1n + len2
lenrws = len1s + len2
lenrwc = len1c + len2
iwork(17) = lenrw
liwm = 1
leniw = 20 + n
leniwc = 20
if (meth .eq. 2) leniwc = leniw
iwork(18) = leniw
if (istate .eq. 1 .and. lrw .lt. lenrwc) go to 617
if (istate .eq. 1 .and. liw .lt. leniwc) go to 618
if (istate .eq. 3 .and. lrw .lt. lenrwc) go to 550
if (istate .eq. 3 .and. liw .lt. leniwc) go to 555
lewt = len1 + 1
insufr = 0
if (lrw .ge. lenrw) go to 65
insufr = 2
lewt = len1c + 1
call xerrwx(
1 60hlsodar- warning.. rwork length is sufficient for now, but ,
1 60, 103, 1, 0, 0, 0, 0, 0.0d0, 0.0d0)
call xerrwx(
1 60h may not be later. integration will proceed anyway. ,
1 60, 103, 1, 0, 0, 0, 0, 0.0d0, 0.0d0)
call xerrwx(
1 50h length needed is lenrw = i1, while lrw = i2.,
1 50, 103, 1, 2, lenrw, lrw, 0, 0.0d0, 0.0d0)
65 lsavf = lewt + n
lacor = lsavf + n
insufi = 0
if (liw .ge. leniw) go to 70
insufi = 2
call xerrwx(
1 60hlsodar- warning.. iwork length is sufficient for now, but ,
1 60, 104, 1, 0, 0, 0, 0, 0.0d0, 0.0d0)
call xerrwx(
1 60h may not be later. integration will proceed anyway. ,
1 60, 104, 1, 0, 0, 0, 0, 0.0d0, 0.0d0)
call xerrwx(
1 50h length needed is leniw = i1, while liw = i2.,
1 50, 104, 1, 2, leniw, liw, 0, 0.0d0, 0.0d0)
70 continue
c check rtol and atol for legality. ------------------------------------
rtoli = rtol(1)
atoli = atol(1)
do 75 i = 1,n
if (itol .ge. 3) rtoli = rtol(i)
if (itol .eq. 2 .or. itol .eq. 4) atoli = atol(i)
if (rtoli .lt. 0.0d0) go to 619
if (atoli .lt. 0.0d0) go to 620
75 continue
if (istate .eq. 1) go to 100
c if istate = 3, set flag to signal parameter changes to stoda. --------
jstart = -1
if (n .eq. nyh) go to 200
c neq was reduced. zero part of yh to avoid undefined references. -----
i1 = lyh + l*nyh
i2 = lyh + (maxord + 1)*nyh - 1
if (i1 .gt. i2) go to 200
do 95 i = i1,i2
95 rwork(i) = 0.0d0
go to 200
c-----------------------------------------------------------------------
c block c.
c the next block is for the initial call only (istate = 1).
c it contains all remaining initializations, the initial call to f,
c and the calculation of the initial step size.
c the error weights in ewt are inverted after being loaded.
c-----------------------------------------------------------------------
100 uround = d1mach(4)
tn = t
tsw = t
maxord = mxordn
if (itask .ne. 4 .and. itask .ne. 5) go to 110
tcrit = rwork(1)
if ((tcrit - tout)*(tout - t) .lt. 0.0d0) go to 625
if (h0 .ne. 0.0d0 .and. (t + h0 - tcrit)*h0 .gt. 0.0d0)
1 h0 = tcrit - t
110 jstart = 0
nhnil = 0
nst = 0
nje = 0
nslast = 0
hu = 0.0d0
nqu = 0
mused = 0
miter = 0
ccmax = 0.3d0
maxcor = 3
msbp = 20
mxncf = 10
c initial call to f. (lf0 points to yh(*,2).) -------------------------
lf0 = lyh + nyh
call f (neq, t, y, rwork(lf0))
nfe = 1
c load the initial value vector in yh. ---------------------------------
do 115 i = 1,n
115 rwork(i+lyh-1) = y(i)
c load and invert the ewt array. (h is temporarily set to 1.0.) -------
nq = 1
h = 1.0d0
call ewset (n, itol, rtol, atol, rwork(lyh), rwork(lewt))
do 120 i = 1,n
if (rwork(i+lewt-1) .le. 0.0d0) go to 621
120 rwork(i+lewt-1) = 1.0d0/rwork(i+lewt-1)
c-----------------------------------------------------------------------
c the coding below computes the step size, h0, to be attempted on the
c first step, unless the user has supplied a value for this.
c first check that tout - t differs significantly from zero.
c a scalar tolerance quantity tol is computed, as max(rtol(i))
c if this is positive, or max(atol(i)/abs(y(i))) otherwise, adjusted
c so as to be between 100*uround and 1.0e-3.
c then the computed value h0 is given by..
c
c h0**(-2) = 1./(tol * w0**2) + tol * (norm(f))**2
c
c where w0 = max ( abs(t), abs(tout) ),
c f = the initial value of the vector f(t,y), and
c norm() = the weighted vector norm used throughout, given by
c the vmnorm function routine, and weighted by the
c tolerances initially loaded into the ewt array.
c the sign of h0 is inferred from the initial values of tout and t.
c abs(h0) is made .le. abs(tout-t) in any case.
c-----------------------------------------------------------------------
if (h0 .ne. 0.0d0) go to 180
tdist = dabs(tout - t)
w0 = dmax1(dabs(t),dabs(tout))
if (tdist .lt. 2.0d0*uround*w0) go to 622
tol = rtol(1)
if (itol .le. 2) go to 140
do 130 i = 1,n
130 tol = dmax1(tol,rtol(i))
140 if (tol .gt. 0.0d0) go to 160
atoli = atol(1)
do 150 i = 1,n
if (itol .eq. 2 .or. itol .eq. 4) atoli = atol(i)
ayi = dabs(y(i))
if (ayi .ne. 0.0d0) tol = dmax1(tol,atoli/ayi)
150 continue
160 tol = dmax1(tol,100.0d0*uround)
tol = dmin1(tol,0.001d0)
sum = vmnorm (n, rwork(lf0), rwork(lewt))
sum = 1.0d0/(tol*w0*w0) + tol*sum**2
h0 = 1.0d0/dsqrt(sum)
h0 = dmin1(h0,tdist)
h0 = dsign(h0,tout-t)
c adjust h0 if necessary to meet hmax bound. ---------------------------
180 rh = dabs(h0)*hmxi
if (rh .gt. 1.0d0) h0 = h0/rh
c load h with h0 and scale yh(*,2) by h0. ------------------------------
h = h0
do 190 i = 1,n
190 rwork(i+lf0-1) = h0*rwork(i+lf0-1)
c
c check for a zero of g at t. ------------------------------------------
irfnd = 0
toutc = tout
if (ngc .eq. 0) go to 270
call rchek (1, g, neq, y, rwork(lyh), nyh,
1 rwork(lg0), rwork(lg1), rwork(lgx), jroot, irt)
if (irt .eq. 0) go to 270
go to 632
c-----------------------------------------------------------------------
c block d.
c the next code block is for continuation calls only (istate = 2 or 3)
c and is to check stop conditions before taking a step.
c first, rchek is called to check for a root within the last step
c taken, other than the last root found there, if any.
c if itask = 2 or 5, and y(tn) has not yet been returned to the user
c because of an intervening root, return through block g.
c-----------------------------------------------------------------------
200 nslast = nst
c
irfp = irfnd
if (ngc .eq. 0) go to 205
if (itask .eq. 1 .or. itask .eq. 4) toutc = tout
call rchek (2, g, neq, y, rwork(lyh), nyh,
1 rwork(lg0), rwork(lg1), rwork(lgx), jroot, irt)
if (irt .ne. 1) go to 205
irfnd = 1
istate = 3
t = t0
go to 425
205 continue
irfnd = 0
if (irfp .eq. 1 .and. tlast .ne. tn .and. itask .eq. 2) go to 400
c
go to (210, 250, 220, 230, 240), itask
210 if ((tn - tout)*h .lt. 0.0d0) go to 250
call intdy (tout, 0, rwork(lyh), nyh, y, iflag)
if (iflag .ne. 0) go to 627
t = tout
go to 420
220 tp = tn - hu*(1.0d0 + 100.0d0*uround)
if ((tp - tout)*h .gt. 0.0d0) go to 623
if ((tn - tout)*h .lt. 0.0d0) go to 250
t = tn
go to 400
230 tcrit = rwork(1)
if ((tn - tcrit)*h .gt. 0.0d0) go to 624
if ((tcrit - tout)*h .lt. 0.0d0) go to 625
if ((tn - tout)*h .lt. 0.0d0) go to 245
call intdy (tout, 0, rwork(lyh), nyh, y, iflag)
if (iflag .ne. 0) go to 627
t = tout
go to 420
240 tcrit = rwork(1)
if ((tn - tcrit)*h .gt. 0.0d0) go to 624
245 hmx = dabs(tn) + dabs(h)
ihit = dabs(tn - tcrit) .le. 100.0d0*uround*hmx
if (ihit) t = tcrit
if (irfp .eq. 1 .and. tlast .ne. tn .and. itask .eq. 5) go to 400
if (ihit) go to 400
tnext = tn + h*(1.0d0 + 4.0d0*uround)
if ((tnext - tcrit)*h .le. 0.0d0) go to 250
h = (tcrit - tn)*(1.0d0 - 4.0d0*uround)
if (istate .eq. 2) jstart = -2
c-----------------------------------------------------------------------
c block e.
c the next block is normally executed for all calls and contains
c the call to the one-step core integrator stoda.
c
c this is a looping point for the integration steps.
c
c first check for too many steps being taken, update ewt (if not at
c start of problem), check for too much accuracy being requested, and
c check for h below the roundoff level in t.
c-----------------------------------------------------------------------
250 continue
if (meth .eq. mused) go to 255
if (insufr .eq. 1) go to 550
if (insufi .eq. 1) go to 555
255 if ((nst-nslast) .ge. mxstep) go to 500
call ewset (n, itol, rtol, atol, rwork(lyh), rwork(lewt))
do 260 i = 1,n
if (rwork(i+lewt-1) .le. 0.0d0) go to 510
260 rwork(i+lewt-1) = 1.0d0/rwork(i+lewt-1)
270 tolsf = uround*vmnorm (n, rwork(lyh), rwork(lewt))
if (tolsf .le. 0.01d0) go to 280
tolsf = tolsf*200.0d0
if (nst .eq. 0) go to 626
go to 520
280 if ((tn + h) .ne. tn) go to 290
nhnil = nhnil + 1
if (nhnil .gt. mxhnil) go to 290
call xerrwx(50hlsodar- warning..internal t (=r1) and h (=r2) are,
1 50, 101, 1, 0, 0, 0, 0, 0.0d0, 0.0d0)
call xerrwx(
1 60h such that in the machine, t + h = t on the next step ,
1 60, 101, 1, 0, 0, 0, 0, 0.0d0, 0.0d0)
call xerrwx(50h (h = step size). solver will continue anyway,
1 50, 101, 1, 0, 0, 0, 2, tn, h)
if (nhnil .lt. mxhnil) go to 290
call xerrwx(50hlsodar- above warning has been issued i1 times. ,
1 50, 102, 1, 0, 0, 0, 0, 0.0d0, 0.0d0)
call xerrwx(50h it will not be issued again for this problem,
1 50, 102, 1, 1, mxhnil, 0, 0, 0.0d0, 0.0d0)
290 continue
c-----------------------------------------------------------------------
c call stoda(neq,y,yh,nyh,yh,ewt,savf,acor,wm,iwm,f,jac,prja,solsy)
c-----------------------------------------------------------------------
call stoda (neq, y, rwork(lyh), nyh, rwork(lyh), rwork(lewt),
1 rwork(lsavf), rwork(lacor), rwork(lwm), iwork(liwm),
2 f, jac, prja, solsy)
kgo = 1 - kflag
go to (300, 530, 540), kgo
c-----------------------------------------------------------------------
c block f.
c the following block handles the case of a successful return from the
c core integrator (kflag = 0).
c if a method switch was just made, record tsw, reset maxord,
c set jstart to -1 to signal stoda to complete the switch,
c and do extra printing of data if ixpr = 1.
c then call rchek to check for a root within the last step.
c then, if no root was found, check for stop conditions.
c-----------------------------------------------------------------------
300 init = 1
if (meth .eq. mused) go to 310
tsw = tn
maxord = mxordn
if (meth .eq. 2) maxord = mxords
if (meth .eq. 2) rwork(lwm) = dsqrt(uround)
insufr = min0(insufr,1)
insufi = min0(insufi,1)
jstart = -1
if (ixpr .eq. 0) go to 310
if (meth .eq. 2) call xerrwx(
1 60hlsodar- a switch to the bdf (stiff) method has occurred ,
1 60, 105, 1, 0, 0, 0, 0, 0.0d0, 0.0d0)
if (meth .eq. 1) call xerrwx(
1 60hlsodar- a switch to the adams (nonstiff) method has occurred,
1 60, 106, 1, 0, 0, 0, 0, 0.0d0, 0.0d0)
call xerrwx(
1 60h at t = r1, tentative step size h = r2, step nst = i1 ,
1 60, 107, 1, 1, nst, 0, 2, tn, h)
310 continue
c
if (ngc .eq. 0) go to 315
call rchek (3, g, neq, y, rwork(lyh), nyh,
1 rwork(lg0), rwork(lg1), rwork(lgx), jroot, irt)
if (irt .ne. 1) go to 315
irfnd = 1
istate = 3
t = t0
go to 425
315 continue
c
go to (320, 400, 330, 340, 350), itask
c itask = 1. if tout has been reached, interpolate. -------------------
320 if ((tn - tout)*h .lt. 0.0d0) go to 250
call intdy (tout, 0, rwork(lyh), nyh, y, iflag)
t = tout
go to 420
c itask = 3. jump to exit if tout was reached. ------------------------
330 if ((tn - tout)*h .ge. 0.0d0) go to 400
go to 250
c itask = 4. see if tout or tcrit was reached. adjust h if necessary.
340 if ((tn - tout)*h .lt. 0.0d0) go to 345
call intdy (tout, 0, rwork(lyh), nyh, y, iflag)
t = tout
go to 420
345 hmx = dabs(tn) + dabs(h)
ihit = dabs(tn - tcrit) .le. 100.0d0*uround*hmx
if (ihit) go to 400
tnext = tn + h*(1.0d0 + 4.0d0*uround)
if ((tnext - tcrit)*h .le. 0.0d0) go to 250
h = (tcrit - tn)*(1.0d0 - 4.0d0*uround)
jstart = -2
go to 250
c itask = 5. see if tcrit was reached and jump to exit. ---------------
350 hmx = dabs(tn) + dabs(h)
ihit = dabs(tn - tcrit) .le. 100.0d0*uround*hmx
c-----------------------------------------------------------------------
c block g.
c the following block handles all successful returns from lsodar.
c if itask .ne. 1, y is loaded from yh and t is set accordingly.
c istate is set to 2, the illegal input counter is zeroed, and the
c optional outputs are loaded into the work arrays before returning.
c if istate = 1 and tout = t, there is a return with no action taken,
c except that if this has happened repeatedly, the run is terminated.
c-----------------------------------------------------------------------
400 do 410 i = 1,n
410 y(i) = rwork(i+lyh-1)
t = tn
if (itask .ne. 4 .and. itask .ne. 5) go to 420
if (ihit) t = tcrit
420 istate = 2
425 continue
illin = 0
rwork(11) = hu
rwork(12) = h
rwork(13) = tn
rwork(15) = tsw
iwork(11) = nst
iwork(12) = nfe
iwork(13) = nje
iwork(14) = nqu
iwork(15) = nq
iwork(19) = mused
iwork(20) = meth
iwork(10) = nge
tlast = t
return
c
430 ntrep = ntrep + 1
if (ntrep .lt. 5) return
call xerrwx(
1 60hlsodar- repeated calls with istate = 1 and tout = t (=r1) ,
1 60, 301, 1, 0, 0, 0, 1, t, 0.0d0)
go to 800
c-----------------------------------------------------------------------
c block h.
c the following block handles all unsuccessful returns other than
c those for illegal input. first the error message routine is called.
c if there was an error test or convergence test failure, imxer is set.
c then y is loaded from yh, t is set to tn, and the illegal input
c counter illin is set to 0. the optional outputs are loaded into
c the work arrays before returning.
c-----------------------------------------------------------------------
c the maximum number of steps was taken before reaching tout. ----------
500 call xerrwx(50hlsodar- at current t (=r1), mxstep (=i1) steps ,
1 50, 201, 1, 0, 0, 0, 0, 0.0d0, 0.0d0)
call xerrwx(50h taken on this call before reaching tout ,
1 50, 201, 1, 1, mxstep, 0, 1, tn, 0.0d0)
istate = -1
go to 580
c ewt(i) .le. 0.0 for some i (not at start of problem). ----------------
510 ewti = rwork(lewt+i-1)
call xerrwx(50hlsodar- at t (=r1), ewt(i1) has become r2 .le. 0.,
1 50, 202, 1, 1, i, 0, 2, tn, ewti)
istate = -6
go to 580
c too much accuracy requested for machine precision. -------------------
520 call xerrwx(50hlsodar- at t (=r1), too much accuracy requested ,
1 50, 203, 1, 0, 0, 0, 0, 0.0d0, 0.0d0)
call xerrwx(50h for precision of machine.. see tolsf (=r2) ,
1 50, 203, 1, 0, 0, 0, 2, tn, tolsf)
rwork(14) = tolsf
istate = -2
go to 580
c kflag = -1. error test failed repeatedly or with abs(h) = hmin. -----
530 call xerrwx(50hlsodar- at t(=r1) and step size h(=r2), the error,
1 50, 204, 1, 0, 0, 0, 0, 0.0d0, 0.0d0)
call xerrwx(50h test failed repeatedly or with abs(h) = hmin,
1 50, 204, 1, 0, 0, 0, 2, tn, h)
istate = -4
go to 560
c kflag = -2. convergence failed repeatedly or with abs(h) = hmin. ----
540 call xerrwx(50hlsodar- at t (=r1) and step size h (=r2), the ,
1 50, 205, 1, 0, 0, 0, 0, 0.0d0, 0.0d0)
call xerrwx(50h corrector convergence failed repeatedly ,
1 50, 205, 1, 0, 0, 0, 0, 0.0d0, 0.0d0)
call xerrwx(30h or with abs(h) = hmin ,
1 30, 205, 1, 0, 0, 0, 2, tn, h)
istate = -5
go to 560
c rwork length too small to proceed. -----------------------------------
550 call xerrwx(50hlsodar- at current t(=r1), rwork length too small,
1 50, 206, 1, 0, 0, 0, 0, 0.0d0, 0.0d0)
call xerrwx(
1 60h to proceed. the integration was otherwise successful.,
1 60, 206, 1, 0, 0, 0, 1, tn, 0.0d0)
istate = -7
go to 580
c iwork length too small to proceed. -----------------------------------
555 call xerrwx(50hlsodar- at current t(=r1), iwork length too small,
1 50, 207, 1, 0, 0, 0, 0, 0.0d0, 0.0d0)
call xerrwx(
1 60h to proceed. the integration was otherwise successful.,
1 60, 207, 1, 0, 0, 0, 1, tn, 0.0d0)
istate = -7
go to 580
c compute imxer if relevant. -------------------------------------------
560 big = 0.0d0
imxer = 1
do 570 i = 1,n
size = dabs(rwork(i+lacor-1)*rwork(i+lewt-1))
if (big .ge. size) go to 570
big = size
imxer = i
570 continue
iwork(16) = imxer
c set y vector, t, illin, and optional outputs. ------------------------
580 do 590 i = 1,n
590 y(i) = rwork(i+lyh-1)
t = tn
illin = 0
rwork(11) = hu
rwork(12) = h
rwork(13) = tn
rwork(15) = tsw
iwork(11) = nst
iwork(12) = nfe
iwork(13) = nje
iwork(14) = nqu
iwork(15) = nq
iwork(19) = mused
iwork(20) = meth
iwork(10) = nge
tlast = t
return
c-----------------------------------------------------------------------
c block i.
c the following block handles all error returns due to illegal input
c (istate = -3), as detected before calling the core integrator.
c first the error message routine is called. then if there have been
c 5 consecutive such returns just before this call to the solver,
c the run is halted.
c-----------------------------------------------------------------------
601 call xerrwx(30hlsodar- istate (=i1) illegal ,
1 30, 1, 1, 1, istate, 0, 0, 0.0d0, 0.0d0)
go to 700
602 call xerrwx(30hlsodar- itask (=i1) illegal ,
1 30, 2, 1, 1, itask, 0, 0, 0.0d0, 0.0d0)
go to 700
603 call xerrwx(50hlsodar- istate .gt. 1 but lsodar not initialized ,
1 50, 3, 1, 0, 0, 0, 0, 0.0d0, 0.0d0)
go to 700
604 call xerrwx(30hlsodar- neq (=i1) .lt. 1 ,
1 30, 4, 1, 1, neq(1), 0, 0, 0.0d0, 0.0d0)
go to 700
605 call xerrwx(50hlsodar- istate = 3 and neq increased (i1 to i2) ,
1 50, 5, 1, 2, n, neq(1), 0, 0.0d0, 0.0d0)
go to 700
606 call xerrwx(30hlsodar- itol (=i1) illegal ,
1 30, 6, 1, 1, itol, 0, 0, 0.0d0, 0.0d0)
go to 700
607 call xerrwx(30hlsodar- iopt (=i1) illegal ,
1 30, 7, 1, 1, iopt, 0, 0, 0.0d0, 0.0d0)
go to 700
608 call xerrwx(30hlsodar- jt (=i1) illegal ,
1 30, 8, 1, 1, jt, 0, 0, 0.0d0, 0.0d0)
go to 700
609 call xerrwx(50hlsodar- ml (=i1) illegal.. .lt.0 or .ge.neq (=i2),
1 50, 9, 1, 2, ml, neq(1), 0, 0.0d0, 0.0d0)
go to 700
610 call xerrwx(50hlsodar- mu (=i1) illegal.. .lt.0 or .ge.neq (=i2),
1 50, 10, 1, 2, mu, neq(1), 0, 0.0d0, 0.0d0)
go to 700
611 call xerrwx(30hlsodar- ixpr (=i1) illegal ,
1 30, 11, 1, 1, ixpr, 0, 0, 0.0d0, 0.0d0)
go to 700
612 call xerrwx(30hlsodar- mxstep (=i1) .lt. 0 ,
1 30, 12, 1, 1, mxstep, 0, 0, 0.0d0, 0.0d0)
go to 700
613 call xerrwx(30hlsodar- mxhnil (=i1) .lt. 0 ,
1 30, 13, 1, 1, mxhnil, 0, 0, 0.0d0, 0.0d0)
go to 700
614 call xerrwx(40hlsodar- tout (=r1) behind t (=r2) ,
1 40, 14, 1, 0, 0, 0, 2, tout, t)
call xerrwx(50h integration direction is given by h0 (=r1) ,
1 50, 14, 1, 0, 0, 0, 1, h0, 0.0d0)
go to 700
615 call xerrwx(30hlsodar- hmax (=r1) .lt. 0.0 ,
1 30, 15, 1, 0, 0, 0, 1, hmax, 0.0d0)
go to 700
616 call xerrwx(30hlsodar- hmin (=r1) .lt. 0.0 ,
1 30, 16, 1, 0, 0, 0, 1, hmin, 0.0d0)
go to 700
617 call xerrwx(
1 60hlsodar- rwork length needed, lenrw (=i1), exceeds lrw (=i2),
1 60, 17, 1, 2, lenrw, lrw, 0, 0.0d0, 0.0d0)
go to 700
618 call xerrwx(
1 60hlsodar- iwork length needed, leniw (=i1), exceeds liw (=i2),
1 60, 18, 1, 2, leniw, liw, 0, 0.0d0, 0.0d0)
go to 700
619 call xerrwx(40hlsodar- rtol(i1) is r1 .lt. 0.0 ,
1 40, 19, 1, 1, i, 0, 1, rtoli, 0.0d0)
go to 700
620 call xerrwx(40hlsodar- atol(i1) is r1 .lt. 0.0 ,
1 40, 20, 1, 1, i, 0, 1, atoli, 0.0d0)
go to 700
621 ewti = rwork(lewt+i-1)
call xerrwx(40hlsodar- ewt(i1) is r1 .le. 0.0 ,
1 40, 21, 1, 1, i, 0, 1, ewti, 0.0d0)
go to 700
622 call xerrwx(
1 60hlsodar- tout (=r1) too close to t(=r2) to start integration,
1 60, 22, 1, 0, 0, 0, 2, tout, t)
go to 700
623 call xerrwx(
1 60hlsodar- itask = i1 and tout (=r1) behind tcur - hu (= r2) ,
1 60, 23, 1, 1, itask, 0, 2, tout, tp)
go to 700
624 call xerrwx(
1 60hlsodar- itask = 4 or 5 and tcrit (=r1) behind tcur (=r2) ,
1 60, 24, 1, 0, 0, 0, 2, tcrit, tn)
go to 700
625 call xerrwx(
1 60hlsodar- itask = 4 or 5 and tcrit (=r1) behind tout (=r2) ,
1 60, 25, 1, 0, 0, 0, 2, tcrit, tout)
go to 700
626 call xerrwx(50hlsodar- at start of problem, too much accuracy ,
1 50, 26, 1, 0, 0, 0, 0, 0.0d0, 0.0d0)
call xerrwx(
1 60h requested for precision of machine.. see tolsf (=r1) ,
1 60, 26, 1, 0, 0, 0, 1, tolsf, 0.0d0)
rwork(14) = tolsf
go to 700
627 call xerrwx(50hlsodar- trouble from intdy. itask = i1, tout = r1,
1 50, 27, 1, 1, itask, 0, 1, tout, 0.0d0)
go to 700
628 call xerrwx(30hlsodar- mxordn (=i1) .lt. 0 ,
1 30, 28, 1, 1, mxordn, 0, 0, 0.0d0, 0.0d0)
go to 700
629 call xerrwx(30hlsodar- mxords (=i1) .lt. 0 ,
1 30, 29, 1, 1, mxords, 0, 0, 0.0d0, 0.0d0)
go to 700
630 call xerrwx(30hlsodar- ng (=i1) .lt. 0 ,
1 30, 30, 1, 1, ng, 0, 0, 0.0d0, 0.0d0)
go to 700
631 call xerrwx(50hlsodar- ng changed (from i1 to i2) illegally, ,
1 50, 31, 1, 0, 0, 0, 0, 0.0d0, 0.0d0)
call xerrwx(50h i.e. not immediately after a root was found ,
1 50, 31, 1, 2, ngc, ng, 0, 0.0d0, 0.0d0)
go to 700
632 call xerrwx(50hlsodar- one or more components of g has a root ,
1 50, 32, 1, 0, 0, 0, 0, 0.0d0, 0.0d0)
call xerrwx(40h too near to the initial point ,
1 40, 32, 1, 0, 0, 0, 0, 0.0d0, 0.0d0)
c
700 if (illin .eq. 5) go to 710
illin = illin + 1
tlast = t
istate = -3
return
710 call xerrwx(50hlsodar- repeated occurrences of illegal input ,
1 50, 302, 1, 0, 0, 0, 0, 0.0d0, 0.0d0)
c
800 call xerrwx(50hlsodar- run aborted.. apparent infinite loop ,
1 50, 303, 2, 0, 0, 0, 0, 0.0d0, 0.0d0)
return
c----------------------- end of subroutine lsodar ----------------------
end
block data
c-----------------------------------------------------------------------
c this data subprogram loads variables into the internal common
c blocks used by lsode and its variants. the variables are
c defined as follows..
c illin = counter for the number of consecutive times the package
c was called with illegal input. the run is stopped when
c illin reaches 5.
c ntrep = counter for the number of consecutive times the package
c was called with istate = 1 and tout = t. the run is
c stopped when ntrep reaches 5.
c mesflg = flag to control printing of error messages. 1 means print,
c 0 means no printing.
c lunit = default value of logical unit number for printing of error
c messages.
c-----------------------------------------------------------------------
integer illin, iduma, ntrep, idumb, iowns, icomm, mesflg, lunit
double precision rownd, rowns, rcomm
common /ls0001/ rownd, rowns(209), rcomm(9),
1 illin, iduma(10), ntrep, idumb(2), iowns(6), icomm(19)
common /eh0001/ mesflg, lunit
data illin/0/, ntrep/0/
data mesflg/1/, lunit/6/
c
c----------------------- end of block data -----------------------------
end