# sdriv3.f

          SUBROUTINE SDRIV3 (N, T, Y, F, NSTATE, TOUT, NTASK, NROOT, EPS,
         8   EWT, IERROR, MINT, MITER, IMPL, ML, MU, MXORD, HMAX, WORK,
         8   LENW, IWORK, LENIW, JACOBN, FA, NDE, MXSTEP, G, USERS, IERFLG)
    C***BEGIN PROLOGUE  SDRIV3
    C***PURPOSE  The function of SDRIV3 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.  Other important options are available.  SDRIV3
    C            uses single precision arithmetic.
    C***LIBRARY   SLATEC (SDRIVE)
    C***CATEGORY  I1A2, I1A1B
    C***TYPE      SINGLE PRECISION (SDRIV3-S, DDRIV3-D, CDRIV3-C)
    C***KEYWORDS  GEAR'S METHOD, INITIAL VALUE PROBLEMS, ODE,
    C             ORDINARY DIFFERENTIAL EQUATIONS, SDRIVE, SINGLE PRECISION,
    C             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.  ABSTRACT  .......................................................
    C
    C    The primary function of SDRIV3 is to solve N ordinary differential
    C    equations of the form dY(I)/dT = F(Y(I),T), given the initial
    C    conditions Y(I) = YI.  The program has options to allow the
    C    solution of both stiff and non-stiff differential equations.  In
    C    addition, SDRIV3 may be used to solve:
    C      1. The initial value problem, A*dY(I)/dT = F(Y(I),T), where A is
    C         a non-singular matrix depending on Y and T.
    C      2. The hybrid differential/algebraic initial value problem,
    C         A*dY(I)/dT = F(Y(I),T), where A is a vector (whose values may
    C         depend upon Y and T) some of whose components will be zero
    C         corresponding to those equations which are algebraic rather
    C         than differential.
    C    SDRIV3 is to be called once for each output point of T.
    C
    C  II.  PARAMETERS  ....................................................
    C
    C    The user should use parameter names in the call sequence of SDRIV3
    C    for those quantities whose value may be altered by SDRIV3.  The
    C    parameters in the call sequence are:
    C
    C    N      = (Input) The number of dependent functions whose solution
    C             is desired.  N must not be altered during a problem.
    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,
    C             JACOBN, FA, USERS, and G.  Thus parameters required by
    C             those routines can be stored in this array in components
    C             N+1 and above.  (Note: Changes by the user to the first
    C             N components of this array will take effect only after a
    C             restart, i.e., after setting NSTATE to 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                   REAL 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 SDRIV3.
    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  SDRIV3.  However, if the user would like
    C             to abort the calculation, i.e., return control to the
    C             program which calls SDRIV3, he should set N to zero.
    C             SDRIV3 will signal this by returning a value of NSTATE
    C             equal to 6 .  Altering the value of N in F has no effect
    C             on the value of N in the call sequence of SDRIV3.
    C
    C    NSTATE = An integer describing the status of integration.  The
    C             meaning of NSTATE is as follows:
    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 NSTATE should be tested by the
    C                  user, but must not be altered.  (As a convenience to
    C                  the user who may wish to put out the initial
    C                  conditions, SDRIV3 can be called with NSTATE=1, and
    C                  TOUT=T.  In this case the program will return with
    C                  NSTATE unchanged, i.e., NSTATE=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                  MXSTEP steps without reaching TOUT.  The user can
    C                  continue the integration by simply calling SDRIV3
    C                  again.
    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 SDRIV3
    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 SDRIV3 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)(Unsuccessful) N has been set to zero in
    C                  SUBROUTINE JACOBN.  See description of JACOBN below.
    C               9  (Output)(Unsuccessful) N has been set to zero in
    C                  SUBROUTINE FA.  See description of FA below.
    C              10  (Output)(Unsuccessful) N has been set to zero in
    C                  SUBROUTINE USERS.  See description of USERS below.
    C              11  (Output)(Successful) For NTASK = 2 or 3, 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 SDRIV3 again.
    C              12  (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 NSTATE to 1 before
    C                  calling SDRIV3 again.  Otherwise the program will
    C                  terminate the run.
    C
    C    TOUT   = (Input) The point at which the solution is desired.  The
    C             position of TOUT relative to T on the first call
    C             determines the direction of integration.
    C
    C    NTASK  = (Input) An index specifying the manner of returning the
    C             solution, according to the following:
    C               NTASK = 1  Means SDRIV3 will integrate past TOUT and
    C                          interpolate the solution.  This is the most
    C                          efficient mode.
    C               NTASK = 2  Means SDRIV3 will return the solution after
    C                          each internal integration step, or at TOUT,
    C                          whichever comes first.  In the latter case,
    C                          the program integrates exactly to TOUT.
    C               NTASK = 3  Means SDRIV3 will adjust its internal step to
    C                          reach TOUT exactly (useful if a singularity
    C                          exists beyond TOUT.)
    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.)  SDRIV3 attempts
    C             to find the value of T at which one of the equations
    C             changes sign.  SDRIV3 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 SDRIV3 increases as EPS
    C             decreases.
    C
    C    EWT    = (Input) Problem zero, i.e., the smallest, nonzero,
    C             physically meaningful value for the solution.  (Array,
    C             possibly of length one.  See following description of
    C             IERROR.)  Setting EWT smaller than necessary can adversely
    C             affect the running time.
    C
    C    IERROR = (Input) Error control indicator.  A value of 3 is
    C             suggested for most problems.  Other choices and detailed
    C             explanations of EWT and IERROR are given below for those
    C             who may need extra flexibility.
    C
    C             These last three input quantities EPS, EWT and IERROR
    C             control the accuracy of the computed solution.  EWT and
    C             IERROR are used internally to compute an array YWT.  One
    C             step error estimates divided by YWT(I) are kept less than
    C             EPS in root mean square norm.
    C                 IERROR (Set by the user) =
    C                 1  Means YWT(I) = 1. (Absolute error control)
    C                                   EWT is ignored.
    C                 2  Means YWT(I) = ABS(Y(I)),  (Relative error control)
    C                                   EWT is ignored.
    C                 3  Means YWT(I) = MAX(ABS(Y(I)), EWT(1)).
    C                 4  Means YWT(I) = MAX(ABS(Y(I)), EWT(I)).
    C                    This choice is useful when the solution components
    C                    have differing scales.
    C                 5  Means YWT(I) = EWT(I).
    C             If IERROR is 3, EWT need only be dimensioned one.
    C             If IERROR is 4 or 5, the user must dimension EWT at least
    C             N, and set its values.
    C
    C    MINT   = (Input) The integration method indicator.
    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.  When using the Adams methods, the
    C                         program uses a value of MITER=0; when using
    C                         the Gear methods, the program uses the value
    C                         of MITER provided by the user.  Only a value
    C                         of IMPL = 0 and a value of MITER = 1, 2, 4, or
    C                         5 is allowed for this option.  The user may
    C                         not alter the value of MINT or MITER without
    C                         restarting, i.e., setting NSTATE to 1.
    C
    C    MITER  = (Input) The iteration method indicator.
    C               MITER = 0  Means functional iteration.  This value is
    C                          suggested for non-stiff problems.
    C               MITER = 1  Means chord method with analytic Jacobian.
    C                          In this case, the user supplies subroutine
    C                          JACOBN (see description below).
    C               MITER = 2  Means chord method with Jacobian calculated
    C                          internally by finite differences.
    C               MITER = 3  Means chord method with corrections computed
    C                          by the user-written routine USERS (see
    C                          description of USERS below.)  This option
    C                          allows all matrix algebra and storage
    C                          decisions to be made by the user.  When using
    C                          a value of MITER = 3, the subroutine FA is
    C                          not required, even if IMPL is not 0.  For
    C                          further information on using this option, see
    C                          Section IV-E below.
    C               MITER = 4  Means the same as MITER = 1 but the A and
    C                          Jacobian matrices are assumed to be banded.
    C               MITER = 5  Means the same as MITER = 2 but the A and
    C                          Jacobian matrices are assumed to be banded.
    C
    C    IMPL   = (Input) The implicit method indicator.
    C               IMPL = 0    Means solving dY(I)/dT = F(Y(I),T).
    C               IMPL = 1    Means solving A*dY(I)/dT = F(Y(I),T), non-
    C                           singular A (see description of FA below.)
    C                           Only MINT = 1 or 2, and MITER = 1, 2, 3, 4,
    C                           or 5 are allowed for this option.
    C               IMPL = 2,3  Means solving certain systems of hybrid
    C                           differential/algebraic equations (see
    C                           description of FA below.)  Only MINT = 2 and
    C                           MITER = 1, 2, 3, 4, or 5, are allowed for
    C                           this option.
    C               The value of IMPL must not be changed during a problem.
    C
    C    ML     = (Input) The lower half-bandwidth in the case of a banded
    C             A or Jacobian matrix.  (I.e., maximum(R-C) for nonzero
    C             A(R,C).)
    C
    C    MU     = (Input) The upper half-bandwidth in the case of a banded
    C             A or Jacobian matrix.  (I.e., maximum(C-R).)
    C
    C    MXORD  = (Input) The maximum order desired. This is .LE. 12 for
    C             the Adams methods and .LE. 5 for the Gear methods.  Normal
    C             value is 12 and 5, respectively.  If MINT is 3, the
    C             maximum order used will be MIN(MXORD, 12) when using the
    C             Adams methods, and MIN(MXORD, 5) when using the Gear
    C             methods.  MXORD must not be altered during a problem.
    C
    C    HMAX   = (Input) The maximum magnitude of the step size that will
    C             be used for the problem.  This is useful for ensuring that
    C             important details are not missed.  If this is not the
    C             case, a large value, such as the interval length, is
    C             suggested.
    C
    C    WORK
    C    LENW   = (Input)
    C             WORK is an array of LENW real 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                       REAL WORK(...)
    C             The following table gives the required minimum value for
    C             the length of WORK, depending on the value of IMPL and
    C             MITER.  LENW should be set to the value used.  The
    C             contents of WORK should not be disturbed between calls to
    C             SDRIV3.
    C
    C      IMPL =   0            1               2             3
    C              ---------------------------------------------------------
    C MITER =  0   (MXORD+4)*N   Not allowed     Not allowed   Not allowed
    C              + 2*NROOT
    C              + 250
    C
    C         1,2  N*N +         2*N*N +         N*N +         N*(N + NDE)
    C              (MXORD+5)*N   (MXORD+5)*N     (MXORD+6)*N   + (MXORD+5)*N
    C              + 2*NROOT     + 2*NROOT       + 2*NROOT     + 2*NROOT
    C              + 250         + 250           + 250         + 250
    C
    C          3   (MXORD+4)*N   (MXORD+4)*N     (MXORD+4)*N   (MXORD+4)*N
    C              + 2*NROOT     + 2*NROOT       + 2*NROOT     + 2*NROOT
    C              + 250         + 250           + 250         + 250
    C
    C         4,5  (2*ML+MU+1)   2*(2*ML+MU+1)   (2*ML+MU+1)   (2*ML+MU+1)*
    C              *N +          *N +            *N +          (N+NDE) +
    C              (MXORD+5)*N   (MXORD+5)*N     (MXORD+6)*N   + (MXORD+5)*N
    C              + 2*NROOT     + 2*NROOT       + 2*NROOT     + 2*NROOT
    C              + 250         + 250           + 250         + 250
    C              ---------------------------------------------------------
    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 MITER is 0 or 3, or
    C               N+50    if MITER is 1, 2, 4, or 5, or MINT is 3,
    C             and LENIW should be set to the value used.  The contents
    C             of IWORK should not be disturbed between calls to SDRIV3.
    C
    C    JACOBN = A subroutine supplied by the user, if MITER is 1 or 4.
    C             If this is the case, the name must be declared EXTERNAL in
    C             the user's calling program.  Given a system of N
    C             differential equations, it is meaningful to speak about
    C             the partial derivative of the I-th right hand side with
    C             respect to the J-th dependent variable.  In general there
    C             are N*N such quantities.  Often however the equations can
    C             be ordered so that the I-th differential equation only
    C             involves dependent variables with index near I, e.g., I+1,
    C             I-2.  Such a system is called banded.  If, for all I, the
    C             I-th equation depends on at most the variables
    C               Y(I-ML), Y(I-ML+1), ... , Y(I), Y(I+1), ... , Y(I+MU)
    C             then we call ML+MU+1 the bandwidth of the system.  In a
    C             banded system many of the partial derivatives above are
    C             automatically zero.  For the cases MITER = 1, 2, 4, and 5,
    C             some of these partials are needed.  For the cases
    C             MITER = 2 and 5 the necessary derivatives are
    C             approximated numerically by SDRIV3, and we only ask the
    C             user to tell SDRIV3 the value of ML and MU if the system
    C             is banded.  For the cases MITER = 1 and 4 the user must
    C             derive these partials algebraically and encode them in
    C             subroutine JACOBN.  By computing these derivatives the
    C             user can often save 20-30 per cent of the computing time.
    C             Usually, however, the accuracy is not much affected and
    C             most users will probably forego this option.  The optional
    C             user-written subroutine JACOBN has the form:
    C                   SUBROUTINE JACOBN (N, T, Y, DFDY, MATDIM, ML, MU)
    C                   REAL Y(*), DFDY(MATDIM,*)
    C                     .
    C                     .
    C                     Calculate values of DFDY
    C                     .
    C                     .
    C                   END (Sample)
    C             Here Y is a vector of length at least N.  The actual
    C             length of Y is determined by the user's declaration in the
    C             program which calls SDRIV3.  Thus the dimensioning of Y in
    C             JACOBN, while required by FORTRAN convention, does not
    C             actually allocate any storage.  When this subroutine is
    C             called, the first N components of Y are intermediate
    C             approximations to the solution components.  The user
    C             should not alter these values.  If the system is not
    C             banded (MITER=1), the partials of the I-th equation with
    C             respect to the J-th dependent function are to be stored in
    C             DFDY(I,J).  Thus partials of the I-th equation are stored
    C             in the I-th row of DFDY.  If the system is banded
    C             (MITER=4), then the partials of the I-th equation with
    C             respect to Y(J) are to be stored in DFDY(K,J), where
    C             K=I-J+MU+1 .  Normally a return from JACOBN passes control
    C             back to SDRIV3.  However, if the user would like to abort
    C             the calculation, i.e., return control to the program which
    C             calls SDRIV3, he should set N to zero.  SDRIV3 will signal
    C             this by returning a value of NSTATE equal to +8(-8).
    C             Altering the value of N in JACOBN has no effect on the
    C             value of N in the call sequence of SDRIV3.
    C
    C    FA     = A subroutine supplied by the user if IMPL is not zero, and
    C             MITER is not 3.  If so, the name must be declared EXTERNAL
    C             in the user's calling program.  This subroutine computes
    C             the array A, where A*dY(I)/dT = F(Y(I),T).
    C             There are three cases:
    C
    C               IMPL=1.
    C               Subroutine FA is of the form:
    C                   SUBROUTINE FA (N, T, Y, A, MATDIM, ML, MU, NDE)
    C                   REAL Y(*), A(MATDIM,*)
    C                     .
    C                     .
    C                     Calculate ALL values of A
    C                     .
    C                     .
    C                   END (Sample)
    C               In this case A is assumed to be a nonsingular matrix,
    C               with the same structure as DFDY (see JACOBN description
    C               above).  Programming considerations prevent complete
    C               generality.  If MITER is 1 or 2, A is assumed to be full
    C               and the user must compute and store all values of
    C               A(I,J), I,J=1, ... ,N.  If MITER is 4 or 5, A is assumed
    C               to be banded with lower and upper half bandwidth ML and
    C               MU.  The left hand side of the I-th equation is a linear
    C               combination of dY(I-ML)/dT, dY(I-ML+1)/dT, ... ,
    C               dY(I)/dT, ... , dY(I+MU-1)/dT, dY(I+MU)/dT.  Thus in the
    C               I-th equation, the coefficient of dY(J)/dT is to be
    C               stored in A(K,J), where K=I-J+MU+1.
    C               NOTE: The array A will be altered between calls to FA.
    C
    C               IMPL=2.
    C               Subroutine FA is of the form:
    C                   SUBROUTINE FA (N, T, Y, A, MATDIM, ML, MU, NDE)
    C                   REAL Y(*), A(*)
    C                     .
    C                     .
    C                     Calculate non-zero values of A(1),...,A(NDE)
    C                     .
    C                     .
    C                   END (Sample)
    C               In this case it is assumed that the system is ordered by
    C               the user so that the differential equations appear
    C               first, and the algebraic equations appear last.  The
    C               algebraic equations must be written in the form:
    C               0 = F(Y(I),T).  When using this option it is up to the
    C               user to provide initial values for the Y(I) that satisfy
    C               the algebraic equations as well as possible.  It is
    C               further assumed that A is a vector of length NDE.  All
    C               of the components of A, which may depend on T, Y(I),
    C               etc., must be set by the user to non-zero values.
    C
    C               IMPL=3.
    C               Subroutine FA is of the form:
    C                   SUBROUTINE FA (N, T, Y, A, MATDIM, ML, MU, NDE)
    C                   REAL Y(*), A(MATDIM,*)
    C                     .
    C                     .
    C                     Calculate ALL values of A
    C                     .
    C                     .
    C                   END (Sample)
    C               In this case A is assumed to be a nonsingular NDE by NDE
    C               matrix with the same structure as DFDY (see JACOBN
    C               description above).  Programming considerations prevent
    C               complete generality.  If MITER is 1 or 2, A is assumed
    C               to be full and the user must compute and store all
    C               values of A(I,J), I,J=1, ... ,NDE.  If MITER is 4 or 5,
    C               A is assumed to be banded with lower and upper half
    C               bandwidths ML and MU.  The left hand side of the I-th
    C               equation is a linear combination of dY(I-ML)/dT,
    C               dY(I-ML+1)/dT, ... , dY(I)/dT, ... , dY(I+MU-1)/dT,
    C               dY(I+MU)/dT.  Thus in the I-th equation, the coefficient
    C               of dY(J)/dT is to be stored in A(K,J), where K=I-J+MU+1.
    C               It is assumed that the system is ordered by the user so
    C               that the differential equations appear first, and the
    C               algebraic equations appear last.  The algebraic
    C               equations must be written in the form 0 = F(Y(I),T).
    C               When using this option it is up to the user to provide
    C               initial values for the Y(I) that satisfy the algebraic
    C               equations as well as possible.
    C               NOTE: For IMPL = 3, the array A will be altered between
    C               calls to FA.
    C             Here Y is a vector of length at least N.  The actual
    C             length of Y is determined by the user's declaration in the
    C             program which calls SDRIV3.  Thus the dimensioning of Y in
    C             FA, while required by FORTRAN convention, does not
    C             actually allocate any storage.  When this subroutine is
    C             called, the first N components of Y are intermediate
    C             approximations to the solution components.  The user
    C             should not alter these values.  FA is always called
    C             immediately after calling F, with the same values of T
    C             and Y.  Normally a return from FA passes control back to
    C             SDRIV3.  However, if the user would like to abort the
    C             calculation, i.e., return control to the program which
    C             calls SDRIV3, he should set N to zero.  SDRIV3 will signal
    C             this by returning a value of NSTATE equal to +9(-9).
    C             Altering the value of N in FA has no effect on the value
    C             of N in the call sequence of SDRIV3.
    C
    C    NDE    = (Input) The number of differential equations.  This is
    C             required only for IMPL = 2 or 3, with NDE .LT. N.
    C
    C    MXSTEP = (Input) The maximum number of internal steps allowed on
    C             one call to SDRIV3.
    C
    C    G      = A real 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 obtain
    C             the value of each of the NROOT equations for which a root
    C             is desired.  G is of the form:
    C                   REAL FUNCTION G (N, T, Y, IROOT)
    C                   REAL 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 SDRIV3.  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  SDRIV3.  However, if the user would like
    C             to abort the calculation, i.e., return control to the
    C             program which calls SDRIV3, he should set N to zero.
    C             SDRIV3 will signal this by returning a value of NSTATE
    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 SDRIV3.
    C
    C    USERS  = A subroutine supplied by the user, if MITER is 3.
    C             If this is the case, the name must be declared EXTERNAL in
    C             the user's calling program.  The routine USERS is called
    C             by SDRIV3 when certain linear systems must be solved.  The
    C             user may choose any method to form, store and solve these
    C             systems in order to obtain the solution result that is
    C             returned to SDRIV3.  In particular, this allows sparse
    C             matrix methods to be used.  The call sequence for this
    C             routine is:
    C
    C                SUBROUTINE USERS (Y, YH, YWT, SAVE1, SAVE2, T, H, EL,
    C               8                  IMPL, N, NDE, IFLAG)
    C                REAL Y(*), YH(*), YWT(*), SAVE1(*),
    C               8     SAVE2(*), T, H, EL
    C
    C             The input variable IFLAG indicates what action is to be
    C             taken.  Subroutine USERS should perform the following
    C             operations, depending on the value of IFLAG and IMPL.
    C
    C               IFLAG = 0
    C                 IMPL = 0.  USERS is not called.
    C                 IMPL = 1, 2 or 3.  Solve the system A*X = SAVE2,
    C                   returning the result in SAVE2.  The array SAVE1 can
    C                   be used as a work array.  For IMPL = 1, there are N
    C                   components to the system, and for IMPL = 2 or 3,
    C                   there are NDE components to the system.
    C
    C               IFLAG = 1
    C                 IMPL = 0.  Compute, decompose and store the matrix
    C                   (I - H*EL*J), where I is the identity matrix and J
    C                   is the Jacobian matrix of the right hand side.  The
    C                   array SAVE1 can be used as a work array.
    C                 IMPL = 1, 2 or 3. Compute, decompose and store the
    C                   matrix (A - H*EL*J).  The array SAVE1 can be used as
    C                   a work array.
    C
    C               IFLAG = 2
    C                 IMPL = 0.   Solve the system
    C                     (I - H*EL*J)*X = H*SAVE2 - YH - SAVE1,
    C                   returning the result in SAVE2.
    C                 IMPL = 1, 2 or 3.  Solve the system
    C                   (A - H*EL*J)*X = H*SAVE2 - A*(YH + SAVE1)
    C                   returning the result in SAVE2.
    C                 The array SAVE1 should not be altered.
    C             If IFLAG is 0 and IMPL is 1 or 2 and the matrix A is
    C             singular, or if IFLAG is 1 and one of the matrices
    C             (I - H*EL*J), (A - H*EL*J) is singular, the INTEGER
    C             variable IFLAG is to be set to -1 before RETURNing.
    C             Normally a return from USERS passes control back to
    C             SDRIV3.  However, if the user would like to abort the
    C             calculation, i.e., return control to the program which
    C             calls SDRIV3, he should set N to zero.  SDRIV3 will signal
    C             this by returning a value of NSTATE equal to +10(-10).
    C             Altering the value of N in USERS has no effect on the
    C             value of N in the call sequence of SDRIV3.
    C
    C    IERFLG = An error flag.  The error number associated with a
    C             diagnostic message (see Section III-A below) is the same
    C             as the corresponding value of IERFLG.  The meaning of
    C             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 NTASK = 2 or 3, 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              24  (Recoverable) MITER is less than 0 or greater than
    C                  5 .
    C              25  (Recoverable) IMPL is less than 0 or greater than 3 .
    C              26  (Recoverable) The value of NSTATE is less than 1 or
    C                  greater than 12 .
    C              27  (Recoverable) EPS is less than zero.
    C              28  (Recoverable) MXORD is not positive.
    C              29  (Recoverable) For MINT = 3, either MITER = 0 or 3, or
    C                  IMPL = 0 .
    C              30  (Recoverable) For MITER = 0, IMPL is not 0 .
    C              31  (Recoverable) For MINT = 1, IMPL is 2 or 3 .
    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              43  (Recoverable)  For IMPL greater than 0, the matrix A
    C                  is singular.
    C             999  (Fatal) The value of NSTATE is 12 .
    C
    C  III.  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  IV.  REMARKS  .......................................................
    C
    C    A. Other routines used:
    C         SDNTP, SDZRO, SDSTP, SDNTL, SDPST, SDCOR, SDCST,
    C         SDPSC, and SDSCL;
    C         SGEFA, SGESL, SGBFA, SGBSL, and SNRM2 (from LINPACK)
    C         R1MACH (from the Bell Laboratories Machine Constants Package)
    C         XERMSG (from the SLATEC Common Math Library)
    C       The last seven routines above, not having been written by the
    C       present authors, are not explicitly part of this package.
    C
    C    B. On any return from SDRIV3 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    C. 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 SDRIV3.
    C
    C    D. Changing parameters during an integration.
    C       The value of NROOT, EPS, EWT, IERROR, MINT, MITER, or HMAX may
    C       be altered by the user between calls to SDRIV3.  For example, if
    C       too much accuracy has been requested (the program returns with
    C       NSTATE = 4 and an increased value of EPS) the user may wish to
    C       increase EPS further.  In general, prudence is necessary when
    C       making changes in parameters since such changes are not
    C       implemented until the next integration step, which is not
    C       necessarily the next call to SDRIV3.  This can happen if the
    C       program has already integrated to a point which is beyond the
    C       new point TOUT.
    C
    C    E. As the price for complete control of matrix algebra, the SDRIV3
    C       USERS option puts all responsibility for Jacobian matrix
    C       evaluation on the user.  It is often useful to approximate
    C       numerically all or part of the Jacobian matrix.  However this
    C       must be done carefully.  The FORTRAN sequence below illustrates
    C       the method we recommend.  It can be inserted directly into
    C       subroutine USERS to approximate Jacobian elements in rows I1
    C       to I2 and columns J1 to J2.
    C              REAL DFDY(N,N), EPSJ, H, R, R1MACH,
    C             8     SAVE1(N), SAVE2(N), T, UROUND, Y(N), YJ, YWT(N)
    C              UROUND = R1MACH(4)
    C              EPSJ = SQRT(UROUND)
    C              DO 30 J = J1,J2
    C                R = EPSJ*MAX(ABS(YWT(J)), ABS(Y(J)))
    C                IF (R .EQ. 0.E0) R = YWT(J)
    C                YJ = Y(J)
    C                Y(J) = Y(J) + R
    C                CALL F (N, T, Y, SAVE1)
    C                IF (N .EQ. 0) RETURN
    C                Y(J) = YJ
    C                DO 20 I = I1,I2
    C         20       DFDY(I,J) = (SAVE1(I) - SAVE2(I))/R
    C         30     CONTINUE
    C       Many problems give rise to structured sparse Jacobians, e.g.,
    C       block banded.  It is possible to approximate them with fewer
    C       function evaluations than the above procedure uses; see Curtis,
    C       Powell and Reid, J. Inst. Maths Applics, (1974), Vol. 13,
    C       pp. 117-119.
    C
    C    F. When any of the routines JACOBN, FA, G, or USERS, is not
    C       required, difficulties associated with unsatisfied externals can
    C       be avoided by using the name of the routine which calculates the
    C       right hand side of the differential equations in place of the
    C       corresponding name in the call sequence of SDRIV3.
    C
    C***REFERENCES  C. W. Gear, Numerical Initial Value Problems in
    C                 Ordinary Differential Equations, Prentice-Hall, 1971.
    C***ROUTINES CALLED  R1MACH, SDNTP, SDSTP, SDZRO, SGBFA, SGBSL, SGEFA,
    C                    SGESL, SNRM2, XERMSG
    C***REVISION HISTORY  (YYMMDD)
    C   790601  DATE WRITTEN
    C   900329  Initial submission to SLATEC.
    C***END PROLOGUE  SDRIV3