# ddriv2.f

          SUBROUTINE DDRIV2 (N, T, Y, F, TOUT, MSTATE, NROOT, EPS, EWT,
         8   MINT, WORK, LENW, IWORK, LENIW, G, IERFLG)
    C***BEGIN PROLOGUE  DDRIV2
    C***PURPOSE  The function of DDRIV2 is to solve N ordinary differential
    C            equations of the form dY(I)/dT = F(Y(I),T), given the
    C            initial conditions Y(I) = YI.  The program has options to
    C            allow the solution of both stiff and non-stiff differential
    C            equations.  DDRIV2 uses double precision arithmetic.
    C***LIBRARY   SLATEC (SDRIVE)
    C***CATEGORY  I1A2, I1A1B
    C***TYPE      DOUBLE PRECISION (SDRIV2-S, DDRIV2-D, CDRIV2-C)
    C***KEYWORDS  DOUBLE PRECISION, GEAR'S METHOD, INITIAL VALUE PROBLEMS,
    C             ODE, ORDINARY DIFFERENTIAL EQUATIONS, SDRIVE, STIFF
    C***AUTHOR  Kahaner, D. K., (NIST)
    C             National Institute of Standards and Technology
    C             Gaithersburg, MD  20899
    C           Sutherland, C. D., (LANL)
    C             Mail Stop D466
    C             Los Alamos National Laboratory
    C             Los Alamos, NM  87545
    C***DESCRIPTION
    C
    C  I.  PARAMETERS  .....................................................
    C
    C       (REMEMBER--To run DDRIV2 correctly in double precision, ALL
    C       non-integer arguments in the call sequence, including
    C       arrays, MUST be declared double precision.)
    C
    C    The user should use parameter names in the call sequence of DDRIV2
    C    for those quantities whose value may be altered by DDRIV2.  The
    C    parameters in the call sequence are:
    C
    C    N      = (Input) The number of differential equations.
    C
    C    T      = The independent variable.  On input for the first call, T
    C             is the initial point.  On output, T is the point at which
    C             the solution is given.
    C
    C    Y      = The vector of dependent variables.  Y is used as input on
    C             the first call, to set the initial values.  On output, Y
    C             is the computed solution vector.  This array Y is passed
    C             in the call sequence of the user-provided routines F and
    C             G.  Thus parameters required by F and G can be stored in
    C             this array in components N+1 and above.  (Note: Changes
    C             by the user to the first N components of this array will
    C             take effect only after a restart, i.e., after setting
    C             MSTATE to +1(-1).)
    C
    C    F      = A subroutine supplied by the user.  The name must be
    C             declared EXTERNAL in the user's calling program.  This
    C             subroutine is of the form:
    C                   SUBROUTINE F (N, T, Y, YDOT)
    C                   DOUBLE PRECISION Y(*), YDOT(*)
    C                     .
    C                     .
    C                   YDOT(1) = ...
    C                     .
    C                     .
    C                   YDOT(N) = ...
    C                   END (Sample)
    C             This computes YDOT = F(Y,T), the right hand side of the
    C             differential equations.  Here Y is a vector of length at
    C             least N.  The actual length of Y is determined by the
    C             user's declaration in the program which calls DDRIV2.
    C             Thus the dimensioning of Y in F, while required by FORTRAN
    C             convention, does not actually allocate any storage.  When
    C             this subroutine is called, the first N components of Y are
    C             intermediate approximations to the solution components.
    C             The user should not alter these values.  Here YDOT is a
    C             vector of length N.  The user should only compute YDOT(I)
    C             for I from 1 to N.  Normally a return from F passes
    C             control back to  DDRIV2.  However, if the user would like
    C             to abort the calculation, i.e., return control to the
    C             program which calls DDRIV2, he should set N to zero.
    C             DDRIV2 will signal this by returning a value of MSTATE
    C             equal to +6(-6).  Altering the value of N in F has no
    C             effect on the value of N in the call sequence of DDRIV2.
    C
    C    TOUT   = (Input) The point at which the solution is desired.
    C
    C    MSTATE = An integer describing the status of integration.  The user
    C             must initialize MSTATE to +1 or -1.  If MSTATE is
    C             positive, the routine will integrate past TOUT and
    C             interpolate the solution.  This is the most efficient
    C             mode.  If MSTATE is negative, the routine will adjust its
    C             internal step to reach TOUT exactly (useful if a
    C             singularity exists beyond TOUT.)  The meaning of the
    C             magnitude of MSTATE:
    C               1  (Input) Means the first call to the routine.  This
    C                  value must be set by the user.  On all subsequent
    C                  calls the value of MSTATE should be tested by the
    C                  user.  Unless DDRIV2 is to be reinitialized, only the
    C                  sign of MSTATE may be changed by the user.  (As a
    C                  convenience to the user who may wish to put out the
    C                  initial conditions, DDRIV2 can be called with
    C                  MSTATE=+1(-1), and TOUT=T.  In this case the program
    C                  will return with MSTATE unchanged, i.e.,
    C                  MSTATE=+1(-1).)
    C               2  (Output) Means a successful integration.  If a normal
    C                  continuation is desired (i.e., a further integration
    C                  in the same direction), simply advance TOUT and call
    C                  again.  All other parameters are automatically set.
    C               3  (Output)(Unsuccessful) Means the integrator has taken
    C                  1000 steps without reaching TOUT.  The user can
    C                  continue the integration by simply calling DDRIV2
    C                  again.  Other than an error in problem setup, the
    C                  most likely cause for this condition is trying to
    C                  integrate a stiff set of equations with the non-stiff
    C                  integrator option. (See description of MINT below.)
    C               4  (Output)(Unsuccessful) Means too much accuracy has
    C                  been requested.  EPS has been increased to a value
    C                  the program estimates is appropriate.  The user can
    C                  continue the integration by simply calling DDRIV2
    C                  again.
    C               5  (Output) A root was found at a point less than TOUT.
    C                  The user can continue the integration toward TOUT by
    C                  simply calling DDRIV2 again.
    C               6  (Output)(Unsuccessful) N has been set to zero in
    C                  SUBROUTINE F.
    C               7  (Output)(Unsuccessful) N has been set to zero in
    C                  FUNCTION G.  See description of G below.
    C               8  (Output)(Successful) For MSTATE negative, T is beyond
    C                  TOUT.  The solution was obtained by interpolation.
    C                  The user can continue the integration by simply
    C                  advancing TOUT and calling DDRIV2 again.
    C               9  (Output)(Unsuccessful) The solution could not be
    C                  obtained.  The value of IERFLG (see description
    C                  below) for a "Recoverable" situation indicates the
    C                  type of difficulty encountered: either an illegal
    C                  value for a parameter or an inability to continue the
    C                  solution.  For this condition the user should take
    C                  corrective action and reset MSTATE to +1(-1) before
    C                  calling DDRIV2 again.  Otherwise the program will
    C                  terminate the run.
    C
    C    NROOT  = (Input) The number of equations whose roots are desired.
    C             If NROOT is zero, the root search is not active.  This
    C             option is useful for obtaining output at points which are
    C             not known in advance, but depend upon the solution, e.g.,
    C             when some solution component takes on a specified value.
    C             The root search is carried out using the user-written
    C             function G (see description of G below.)  DDRIV2 attempts
    C             to find the value of T at which one of the equations
    C             changes sign.  DDRIV2 can find at most one root per
    C             equation per internal integration step, and will then
    C             return the solution either at TOUT or at a root, whichever
    C             occurs first in the direction of integration.  The initial
    C             point is never reported as a root.  The index of the
    C             equation whose root is being reported is stored in the
    C             sixth element of IWORK.
    C             NOTE: NROOT is never altered by this program.
    C
    C    EPS    = On input, the requested relative accuracy in all solution
    C             components.  EPS = 0 is allowed.  On output, the adjusted
    C             relative accuracy if the input value was too small.  The
    C             value of EPS should be set as large as is reasonable,
    C             because the amount of work done by DDRIV2 increases as
    C             EPS decreases.
    C
    C    EWT    = (Input) Problem zero, i.e., the smallest physically
    C             meaningful value for the solution.  This is used inter-
    C             nally to compute an array YWT(I) = MAX(ABS(Y(I)), EWT).
    C             One step error estimates divided by YWT(I) are kept less
    C             than EPS.  Setting EWT to zero provides pure relative
    C             error control.  However, setting EWT smaller than
    C             necessary can adversely affect the running time.
    C
    C    MINT   = (Input) The integration method flag.
    C               MINT = 1  Means the Adams methods, and is used for
    C                         non-stiff problems.
    C               MINT = 2  Means the stiff methods of Gear (i.e., the
    C                         backward differentiation formulas), and is
    C                         used for stiff problems.
    C               MINT = 3  Means the program dynamically selects the
    C                         Adams methods when the problem is non-stiff
    C                         and the Gear methods when the problem is
    C                         stiff.
    C             MINT may not be changed without restarting, i.e., setting
    C             the magnitude of MSTATE to 1.
    C
    C    WORK
    C    LENW   = (Input)
    C             WORK is an array of LENW double precision words used
    C             internally for temporary storage.  The user must allocate
    C             space for this array in the calling program by a statement
    C             such as
    C                       DOUBLE PRECISION WORK(...)
    C             The length of WORK should be at least
    C               16*N + 2*NROOT + 250         if MINT is 1, or
    C               N*N + 10*N + 2*NROOT + 250   if MINT is 2, or
    C               N*N + 17*N + 2*NROOT + 250   if MINT is 3,
    C             and LENW should be set to the value used.  The contents of
    C             WORK should not be disturbed between calls to DDRIV2.
    C
    C    IWORK
    C    LENIW  = (Input)
    C             IWORK is an integer array of length LENIW used internally
    C             for temporary storage.  The user must allocate space for
    C             this array in the calling program by a statement such as
    C                       INTEGER IWORK(...)
    C             The length of IWORK should be at least
    C               50      if MINT is 1, or
    C               N+50    if MINT is 2 or 3,
    C             and LENIW should be set to the value used.  The contents
    C             of IWORK should not be disturbed between calls to DDRIV2.
    C
    C    G      = A double precision FORTRAN function supplied by the user
    C             if NROOT is not 0.  In this case, the name must be
    C             declared EXTERNAL in the user's calling program.  G is
    C             repeatedly called with different values of IROOT to
    C             obtain the value of each of the NROOT equations for which
    C             a root is desired.  G is of the form:
    C                   DOUBLE PRECISION FUNCTION G (N, T, Y, IROOT)
    C                   DOUBLE PRECISION Y(*)
    C                   GO TO (10, ...), IROOT
    C              10   G = ...
    C                     .
    C                     .
    C                   END (Sample)
    C             Here, Y is a vector of length at least N, whose first N
    C             components are the solution components at the point T.
    C             The user should not alter these values.  The actual length
    C             of Y is determined by the user's declaration in the
    C             program which calls DDRIV2.  Thus the dimensioning of Y in
    C             G, while required by FORTRAN convention, does not actually
    C             allocate any storage.  Normally a return from G passes
    C             control back to  DDRIV2.  However, if the user would like
    C             to abort the calculation, i.e., return control to the
    C             program which calls DDRIV2, he should set N to zero.
    C             DDRIV2 will signal this by returning a value of MSTATE
    C             equal to +7(-7).  In this case, the index of the equation
    C             being evaluated is stored in the sixth element of IWORK.
    C             Altering the value of N in G has no effect on the value of
    C             N in the call sequence of DDRIV2.
    C
    C    IERFLG = An error flag.  The error number associated with a
    C             diagnostic message (see Section II-A below) is the same as
    C             the corresponding value of IERFLG.  The meaning of IERFLG:
    C               0  The routine completed successfully. (No message is
    C                  issued.)
    C               3  (Warning) The number of steps required to reach TOUT
    C                  exceeds MXSTEP.
    C               4  (Warning) The value of EPS is too small.
    C              11  (Warning) For MSTATE negative, T is beyond TOUT.
    C                  The solution was obtained by interpolation.
    C              15  (Warning) The integration step size is below the
    C                  roundoff level of T.  (The program issues this
    C                  message as a warning but does not return control to
    C                  the user.)
    C              22  (Recoverable) N is not positive.
    C              23  (Recoverable) MINT is less than 1 or greater than 3 .
    C              26  (Recoverable) The magnitude of MSTATE is either 0 or
    C                  greater than 9 .
    C              27  (Recoverable) EPS is less than zero.
    C              32  (Recoverable) Insufficient storage has been allocated
    C                  for the WORK array.
    C              33  (Recoverable) Insufficient storage has been allocated
    C                  for the IWORK array.
    C              41  (Recoverable) The integration step size has gone
    C                  to zero.
    C              42  (Recoverable) The integration step size has been
    C                  reduced about 50 times without advancing the
    C                  solution.  The problem setup may not be correct.
    C             999  (Fatal) The magnitude of MSTATE is 9 .
    C
    C  II.  OTHER COMMUNICATION TO THE USER  ...............................
    C
    C    A. The solver communicates to the user through the parameters
    C       above.  In addition it writes diagnostic messages through the
    C       standard error handling program XERMSG.  A complete description
    C       of XERMSG is given in "Guide to the SLATEC Common Mathematical
    C       Library" by Kirby W. Fong et al..  At installations which do not
    C       have this error handling package the short but serviceable
    C       routine, XERMSG, available with this package, can be used.  That
    C       program uses the file named OUTPUT to transmit messages.
    C
    C    B. The first three elements of WORK and the first five elements of
    C       IWORK will contain the following statistical data:
    C         AVGH     The average step size used.
    C         HUSED    The step size last used (successfully).
    C         AVGORD   The average order used.
    C         IMXERR   The index of the element of the solution vector that
    C                  contributed most to the last error test.
    C         NQUSED   The order last used (successfully).
    C         NSTEP    The number of steps taken since last initialization.
    C         NFE      The number of evaluations of the right hand side.
    C         NJE      The number of evaluations of the Jacobian matrix.
    C
    C  III.  REMARKS  ......................................................
    C
    C    A. On any return from DDRIV2 all information necessary to continue
    C       the calculation is contained in the call sequence parameters,
    C       including the work arrays.  Thus it is possible to suspend one
    C       problem, integrate another, and then return to the first.
    C
    C    B. If this package is to be used in an overlay situation, the user
    C       must declare in the primary overlay the variables in the call
    C       sequence to DDRIV2.
    C
    C    C. When the routine G is not required, difficulties associated with
    C       an unsatisfied external can be avoided by using the name of the
    C       routine which calculates the right hand side of the differential
    C       equations in place of G in the call sequence of DDRIV2.
    C
    C  IV.  USAGE  .........................................................
    C
    C               PROGRAM SAMPLE
    C               EXTERNAL F
    C               PARAMETER(MINT = 1, NROOT = 0, N = ...,
    C              8          LENW = 16*N + 2*NROOT + 250, LENIW = 50)
    C         C                                 N is the number of equations
    C               DOUBLE PRECISION EPS, EWT, T, TOUT, WORK(LENW), Y(N)
    C               INTEGER IWORK(LENIW)
    C               OPEN(FILE='TAPE6', UNIT=6, STATUS='NEW')
    C         C                                                Initial point
    C               T = 0.
    C         C                                       Set initial conditions
    C               DO 10 I = 1,N
    C          10     Y(I) = ...
    C               TOUT = T
    C               EWT = ...
    C               MSTATE = 1
    C               EPS = ...
    C          20   CALL DDRIV2 (N, T, Y, F, TOUT, MSTATE, NROOT, EPS, EWT,
    C              8             MINT, WORK, LENW, IWORK, LENIW, F, IERFLG)
    C         C                                 Next to last argument is not
    C         C                                    F if rootfinding is used.
    C               IF (MSTATE .GT. 2) STOP
    C               WRITE(6, 100) TOUT, (Y(I), I=1,N)
    C               TOUT = TOUT + 1.
    C               IF (TOUT .LE. 10.) GO TO 20
    C          100  FORMAT(...)
    C               END (Sample)
    C
    C***REFERENCES  C. W. Gear, Numerical Initial Value Problems in
    C                 Ordinary Differential Equations, Prentice-Hall, 1971.
    C***ROUTINES CALLED  DDRIV3, XERMSG
    C***REVISION HISTORY  (YYMMDD)
    C   790601  DATE WRITTEN
    C   900329  Initial submission to SLATEC.
    C***END PROLOGUE  DDRIV2