dlsodkr Subroutine

public subroutine dlsodkr(f, Neq, Y, T, Tout, Itol, Rtol, Atol, Itask, Istate, Iopt, Rwork, Lrw, Iwork, Liw, jac, psol, Mf, g, Ng, Jroot)

Synopsis

DLSODKR: Livermore Solver for Ordinary Differential equations, with preconditioned Krylov iteration methods for the Newton correction linear systems, and with Rootfinding.

DLSODKR solves the initial value problem for stiff or nonstiff systems of first order ODEs,

     dy/dt = f(t,y),  or, in component form,
     dy(i)/dt = f(i) = f(i,t,y(1),y(2),...,y(NEQ)) (i = 1,...,NEQ).

At the same time, it locates the roots of any of a set of functions

     g(i) = g(i,t,y(1),...,y(NEQ))  (i = 1,...,ng).

This version is in double precision.

Introduction.

This is a modification of the DLSODE package, and differs from it in five ways: 1. It uses various preconditioned Krylov subspace iteration methods for the linear algebraic systems that arise in the case of stiff systems. See the introductory notes below. 2. It does automatic switching between functional (fixpoint) iteration and Newton iteration in the corrector iteration. 3. It finds the root of at least one of a set of constraint functions g(i) of the independent and dependent variables. It finds only those roots for which some g(i), as a function of t, changes sign in the interval of integration. It then returns the solution at the root, if that occurs sooner than the specified stop condition, and otherwise returns the solution according the specified stop condition. 4. It supplies to JAC an input flag, JOK, which indicates whether JAC may (optionally) bypass the evaluation of Jacobian matrix data and instead process saved data (with the current value of scalar hl0). 5. It contains a new subroutine that calculates the initial step size to be attempted.

Introduction to the Krylov methods in DLSODKR:

The linear systems that must be solved have the form

   A * x  = b,  where  A = identity - hl0 * (df/dy) .

Here hl0 is a scalar, and df/dy is the Jacobian matrix of partial derivatives of f (NEQ by NEQ).

The particular Krylov method is chosen by setting the second digit, MITER, in the method flag MF. Currently, the values of MITER have the following meanings:

The user can apply preconditioning to the linear system A*x = b, by means of arbitrary matrices (the preconditioners).

In the case of SPIOM and SPIGMR, one can apply left and right preconditioners P1 and P2, and the basic iterative method is then applied to the matrix (P1-inverse)*A*(P2-inverse) instead of to the matrix A. The product P1*P2 should be an approximation to matrix A such that linear systems with P1 or P2 are easier to solve than with A. Preconditioning from the left only or right only means using P2 = identity or P1 = identity, respectively.

In the case of the PCG and PCGS methods, there is only one preconditioner matrix P (but it can be the product of more than one). It should approximate the matrix A but allow for relatively easy solution of linear systems with coefficient matrix P. For PCG, P should be positive definite symmetric, or nearly so, and for PCGS, the scaled preconditioner D-inverse * P * D should be symmetric or nearly so.

If the Jacobian J = df/dy splits in a natural way into a sum J = J1 + J2, then one possible choice of preconditioners is P1 = identity - hl0 * J1 and P2 = identity - hl0 * J2 provided each of these is easy to solve (or approximately solve).

Summary of Usage.

Communication between the user and the DLSODKR package, for normal situations, is summarized here. This summary describes only a subset of the full set of options available. See the full description for details, including optional communication, nonstandard options, and instructions for special situations. See also the demonstration program distributed with this solver.

A. First provide a subroutine of the form:

               SUBROUTINE F (NEQ, T, Y, YDOT)
               DOUBLE PRECISION T, Y(*), YDOT(*)

which supplies the vector function f by loading YDOT(i) with f(i).

B. Provide a subroutine of the form:

               SUBROUTINE G (NEQ, T, Y, NG, GOUT)
               DOUBLE PRECISION T, Y(*), GOUT(NG)

which supplies the vector function g by loading GOUT(i) with g(i), the i-th constraint function whose root is sought.

C. Next determine (or guess) whether or not the problem is stiff. Stiffness occurs when the Jacobian matrix df/dy has an eigenvalue whose real part is negative and large in magnitude, compared to the reciprocal of the t span of interest. If the problem is nonstiff, use a method flag MF = 10. If it is stiff, MF should be between 21 and 24, or possibly 29. MF = 22 is generally the best choice. Use 23 or 24 only if symmetry is present. Use MF = 29 if the complete linear system solution is to be provided by the user. The following four parameters must also be set.

  IWORK(1) = LWP  = length of real array WP for preconditioning.
  IWORK(2) = LIWP = length of integer array IWP for preconditioning.
  IWORK(3) = JPRE = preconditioner type flag:
                  = 0 for no preconditioning (P1 = P2 = P = identity)
                  = 1 for left-only preconditioning (P2 = identity)
                  = 2 for right-only preconditioning (P1 = identity)
                  = 3 for two-sided preconditioning (and PCG or PCGS)
  IWORK(4) = JACFLG = flag for whether JAC is called.
                    = 0 if JAC is not to be called,
                    = 1 if JAC is to be called.

Use JACFLG = 1 if JAC computes any nonconstant data for use in preconditioning, such as Jacobian elements. The arrays WP and IWP are work arrays under the user’s control, for use in the routines that perform preconditioning operations.

D. If the problem is stiff, you must supply two routines that deal with the preconditioning of the linear systems to be solved.

These are as follows:

     SUBROUTINE JAC (F, NEQ, T, Y, YSV, REWT, FTY,V,HL0,JOK,WP,IWP,IER)
     DOUBLE PRECISION T, Y(*), YSV(*), REWT(*), FTY(*), V(*), HL0,WP(*)
     INTEGER IWP(*)

This routine must evaluate and preprocess any parts of the Jacobian matrix df/dy involved in the preconditioners P1, P2, P. The Y and FTY arrays contain the current values of y and f(t,y), respectively, and YSV also contains the current value of y. The array V is work space of length NEQ.

JAC must multiply all computed Jacobian elements by the scalar -HL0, add the identity matrix, and do any factorization operations called for, in preparation for solving linear systems with a coefficient matrix of P1, P2, or P. The matrix P1*P2 or P should be an approximation to identity - hl0 * (df/dy). JAC should return IER = 0 if successful, and IER .ne. 0 if not. (If IER .ne. 0, a smaller time step will be tried.)

JAC may alter Y and V, but not YSV, REWT, FTY, or HL0.

The JOK argument can be ignored (or see full description below).

     SUBROUTINE PSOL (NEQ, T, Y, FTY, WK, HL0, WP, IWP, B, LR, IER)
     DOUBLE PRECISION T, Y(*), FTY(*), WK(*), HL0, WP(*), B(*)
     INTEGER IWP(*)

This routine must solve a linear system with B as right-hand side and one of the preconditioning matrices, P1, P2, or P, as coefficient matrix, and return the solution vector in B.

LR is a flag concerning left vs right preconditioning, input to PSOL. PSOL is to use P1 if LR = 1 and P2 if LR = 2.

In the case of the PCG or PCGS method, LR will be 3, and PSOL should solve the system P*x = B with the preconditioner matrix P.

In the case MF = 29 (no Krylov iteration), LR will be 0, and PSOL is to return in B the desired approximate solution to A * x = B, where A = identity - hl0 * (df/dy).

PSOL can use data generated in the JAC routine and stored in WP and IWP. WK is a work array of length NEQ.

The argument HL0 is the current value of the scalar appearing in the linear system. If the old value, at the time of the last JAC call, is needed, it must have been saved by JAC in WP.

on return, PSOL should set the error flag IER as follows:

       IER = 0 if PSOL was successful,
       IER .gt. 0 if a recoverable error occurred, meaning that the
              time step will be retried,
       IER .lt. 0 if an unrecoverable error occurred, meaning that the
              solver is to stop immediately.

E. Write a main program which calls Subroutine DLSODKR once for each point at which answers are desired. This should also provide for possible use of logical unit 6 for output of error messages by DLSODKR. On the first call to DLSODKR, supply arguments as follows:

F

name of subroutine for right-hand side vector f. This name must be declared External in calling program.

NEQ

number of first order ODEs.

Y

array of initial values, of length NEQ.

T

the initial value of the independent variable.

TOUT

first point where output is desired (.ne. T).

ITOL

1 or 2 according as ATOL (below) is a scalar or array.

RTOL

relative tolerance parameter (scalar).

ATOL

absolute tolerance parameter (scalar or array). The estimated local error in y(i) will be controlled so as to be roughly less (in magnitude) than EWT(i) = RTOL*ABS(Y(i)) + ATOL if ITOL = 1, or EWT(i) = RTOL*ABS(Y(i)) + ATOL(i) if ITOL = 2. Thus the local error test passes if, in each component, either the absolute error is less than ATOL (or ATOL(i)), or the relative error is less than RTOL. Use RTOL = 0.0 for pure absolute error control, and use ATOL = 0.0 (or ATOL(i) = 0.0) for pure relative error control. Caution: Actual (global) errors may exceed these local tolerances, so choose them conservatively.

ITASK

1 for normal computation of output values of y at t = TOUT.

ISTATE

integer flag (input and output). Set ISTATE = 1.

IOPT

0 to indicate no optional inputs used.

RWORK

real work array of length at least:

        20 + 16*NEQ + 3*NG           for MF = 10,
        45 + 17*NEQ + 3*NG + LWP     for MF = 21,
        61 + 17*NEQ + 3*NG + LWP     for MF = 22,
        20 + 15*NEQ + 3*NG + LWP     for MF = 23 or 24,
        20 + 12*NEQ + 3*NG + LWP     for MF = 29.

LRW

declared length of RWORK (in user’s dimension).

IWORK

integer work array of length at least:

        30            for MF = 10,
        35 + LIWP     for MF = 21,
        30 + LIWP     for MF = 22, 23, 24, or 29.

LIW

declared length of IWORK (in user’s dimension).

JAC,PSOL

names of subroutines for preconditioning. These names must be declared External in the calling program.

MF

method flag. Standard values are:

G

name of subroutine for constraint functions, whose roots are desired during the integration. This name must be declared External in calling program.

NG

number of constraint functions g(i). If there are none, set NG = 0, and pass a dummy name for G.

JROOT

integer array of length NG for output of root information. See next paragraph.

Note that the main program must declare arrays Y, RWORK, IWORK, JROOT, and possibly ATOL.

F. The output from the first call (or any call) is:

Y

array of computed values of y(t) vector.

T

corresponding value of independent variable (normally TOUT).

ISTATE

values and meanings:

JROOT

array showing roots found if ISTATE = 3 on return. JROOT(i) = 1 if g(i) has a root at T, or 0 otherwise.

G. To continue the integration after a successful return, proceed as follows:

1. If ISTATE = 2 on return, reset TOUT and call DLSODKR again.
2. If ISTATE = 3 on return, reset ISTATE to 2 and call DLSODKR again.

In either case, no other parameters need be reset.

Full Description of User Interface to DLSODKR.

The user interface to DLSODKR consists of the following parts.

1. The call sequence to Subroutine DLSODKR, which is a driver routine for the solver. This includes descriptions of both the call sequence arguments and of user-supplied routines. Following these descriptions is a description of optional inputs available through the call sequence, and then a description of optional outputs (in the work arrays).

2. Descriptions of other routines in the DLSODKR package that may be (optionally) called by the user. These provide the ability to alter error message handling, save and restore the internal Common, and obtain specified derivatives of the solution y(t).

3. Descriptions of Common blocks to be declared in overlay or similar environments, or to be saved when doing an interrupt of the problem and continued solution later.

4. Description of two routines in the DLSODKR package, either of which the user may replace with his/her own version, if desired. These relate to the measurement of errors.

Part 1. Call Sequence.

The call sequence parameters used for input only are F, NEQ, TOUT, ITOL, RTOL, ATOL, ITASK, IOPT, LRW, LIW, JAC, PSOL, MF, G, and NG,

that used only for output is JROOT, and those used for both input and output are Y, T, ISTATE.

The work arrays RWORK and IWORK are also used for conditional and optional inputs and optional outputs. (The term output here refers to the return from Subroutine DLSODKR to the user’s calling program.)

The legality of input parameters will be thoroughly checked on the initial call for the problem, but not checked thereafter unless a change in input parameters is flagged by ISTATE = 3 on input.

The descriptions of the call arguments are as follows.

F

the name of the user-supplied subroutine defining the ODE system. The system must be put in the first-order form dy/dt = f(t,y), where f is a vector-valued function of the scalar t and the vector y. Subroutine F is to compute the function f. It is to have the form

           SUBROUTINE F (NEQ, T, Y, YDOT)
           DOUBLE PRECISION T, Y(*), YDOT(*)

where NEQ, T, and Y are input, and the array YDOT = f(t,y) is output. Y and YDOT are arrays of length NEQ. Subroutine F should not alter Y(1),…,Y(NEQ). F must be declared External in the calling program.

Subroutine F may access user-defined quantities in NEQ(2),… and/or in Y(NEQ(1)+1),… if NEQ is an array (dimensioned in F) and/or Y has length exceeding NEQ(1). See the descriptions of NEQ and Y below.

If quantities computed in the F routine are needed externally to DLSODKR, an extra call to F should be made for this purpose, for consistent and accurate results. If only the derivative dy/dt is needed, use DINTDY instead.

NEQ

the size of the ODE system (number of first order ordinary differential equations). Used only for input. NEQ may be decreased, but not increased, during the problem. If NEQ is decreased (with ISTATE = 3 on input), the remaining components of Y should be left undisturbed, if these are to be accessed in the user-supplied routines.

Normally, NEQ is a scalar, and it is generally referred to as a scalar in this user interface description. However, NEQ may be an array, with NEQ(1) set to the system size. (The DLSODKR package accesses only NEQ(1).) In either case, this parameter is passed as the NEQ argument in all calls to the user-supplied routines. Hence, if it is an array, locations NEQ(2),… may be used to store other integer data and pass it to the user-supplied routines. Each such routine must include NEQ in a Dimension statement in that case.

Y

a real array for the vector of dependent variables, of length NEQ or more. Used for both input and output on the first call (ISTATE = 1), and only for output on other calls. On the first call, Y must contain the vector of initial values. On output, Y contains the computed solution vector, evaluated at T. If desired, the Y array may be used for other purposes between calls to the solver.

This array is passed as the Y argument in all calls to F, G, JAC, and PSOL. Hence its length may exceed NEQ, and locations Y(NEQ+1),… may be used to store other real data and pass it to the user-supplied routines. (The DLSODKR package accesses only Y(1),…,Y(NEQ).)

T

the independent variable. On input, T is used only on the first call, as the initial point of the integration. On output, after each call, T is the value at which a computed solution y is evaluated (usually the same as TOUT). If a root was found, T is the computed location of the root reached first, on output. On an error return, T is the farthest point reached.

TOUT

the next value of t at which a computed solution is desired. Used only for input.

When starting the problem (ISTATE = 1), TOUT may be equal to T for one call, then should .ne. T for the next call. For the initial T, an input value of TOUT .ne. T is used in order to determine the direction of the integration (i.e. the algebraic sign of the step sizes) and the rough scale of the problem. Integration in either direction (forward or backward in t) is permitted.

If ITASK = 2 or 5 (one-step modes), TOUT is ignored after the first call (i.e. the first call with TOUT .ne. T). Otherwise, TOUT is required on every call.

If ITASK = 1, 3, or 4, the values of TOUT need not be monotone, but a value of TOUT which backs up is limited to the current internal T interval, whose endpoints are TCUR - HU and TCUR (see optional outputs, below, for TCUR and HU).

ITOL

an indicator for the type of error control. See description below under ATOL. Used only for input.

RTOL

a relative error tolerance parameter, either a scalar or an array of length NEQ. See description below under ATOL. Input only.

ATOL

an absolute error tolerance parameter, either a scalar or an array of length NEQ. Input only.

The input parameters ITOL, RTOL, and ATOL determine the error control performed by the solver. The solver will control the vector E = (E(i)) of estimated local errors in y, according to an inequality of the form

                    RMS-norm of ( E(i)/EWT(i) )   .le.   1,
        where       EWT(i) = RTOL(i)*ABS(Y(i)) + ATOL(i),

and the RMS-norm (root-mean-square norm) here is RMS-norm(v) = SQRT(sum v(i)**2 / NEQ). Here EWT = (EWT(i)) is a vector of weights which must always be positive, and the values of RTOL and ATOL should all be non-negative. The following table gives the types (scalar/array) of RTOL and ATOL, and the corresponding form of EWT(i).

When either of these parameters is a scalar, it need not be dimensioned in the user’s calling program.

If none of the above choices (with ITOL, RTOL, and ATOL fixed throughout the problem) is suitable, more general error controls can be obtained by substituting user-supplied routines for the setting of EWT and/or for the norm calculation. See Part 4 below.

If global errors are to be estimated by making a repeated run on the same problem with smaller tolerances, then all components of RTOL and ATOL (i.e. of EWT) should be scaled down uniformly.

ITASK

an index specifying the task to be performed. Input only. ITASK has the following values and meanings.

Note: If ITASK = 4 or 5 and the solver reaches TCRIT (within roundoff), it will return T = TCRIT (exactly) to indicate this (unless ITASK = 4 and TOUT comes before TCRIT, in which case answers at T = TOUT are returned first).

ISTATE

an index used for input and output to specify the the state of the calculation.

On input, the values of ISTATE are as follows.

On output, ISTATE has the following values and meanings.

Note: Since the normal output value of ISTATE is 2, it does not need to be reset for normal continuation. Also, since a negative input value of ISTATE will be regarded as illegal, a negative output value requires the user to change it, and possibly other inputs, before calling the solver again.

IOPT

an integer flag to specify whether or not any optional inputs are being used on this call. Input only. The optional inputs are listed separately below.

  IOPT = 0 means no optional inputs are being used.
           Default values will be used in all cases.
  IOPT = 1 means one or more optional inputs are being used.

RWORK

a real working array (double precision).

      The length of RWORK must be at least
         20 + NYH*(MAXORD+1) + 3*NEQ + 3*NG + LENLS + LWP    where
      NYH    = the initial value of NEQ,
      MAXORD = 12 (if METH = 1) or 5 (if METH = 2) (unless a
               smaller value is given as an optional input),
      LENLS = length of work space for linear system (Krylov)
              method, excluding preconditioning:
        LENLS = 0                               if MITER = 0,
        LENLS = NEQ*(MAXL+3) + MAXL**2          if MITER = 1,
        LENLS = NEQ*(MAXL+3+MIN(1,MAXL-KMP))
             + (MAXL+3)*MAXL + 1                if MITER = 2,
        LENLS = 6*NEQ                           if MITER = 3 or 4,
        LENLS = 3*NEQ                           if MITER = 9.
      (See the MF description for METH and MITER, and the
      list of optional inputs for MAXL and KMP.)
      LWP = length of real user work space for preconditioning
      (see JAC/PSOL).
      Thus if default values are used and NEQ is constant,
      this length is:
         20 + 16*NEQ + 3*NG           for MF = 10,
         45 + 24*NEQ + 3*NG + LWP     for MF = 11,
         61 + 24*NEQ + 3*NG + LWP     for MF = 12,
         20 + 22*NEQ + 3*NG + LWP     for MF = 13 or 14,
         20 + 19*NEQ + 3*NG + LWP     for MF = 19,
         20 + 9*NEQ + 3*NG            for MF = 20,
         45 + 17*NEQ + 3*NG + LWP     for MF = 21,
         61 + 17*NEQ + 3*NG + LWP     for MF = 22,
         20 + 15*NEQ + 3*NG + LWP     for MF = 23 or 24,
         20 + 12*NEQ + 3*NG + LWP     for MF = 29.

The first 20 words of RWORK are reserved for conditional and optional inputs and optional outputs.

The following word in RWORK is a conditional input:

RWORK(1) = TCRIT = critical value of t which the solver is not to overshoot. Required if ITASK is 4 or 5, and ignored otherwise. (See ITASK.)

LRW

the length of the array RWORK, as declared by the user. (This will be checked by the solver.)

IWORK

an integer work array. The length of IWORK must be at least 30 if MITER = 0 (MF = 10 or 20), 30 + MAXL + LIWP if MITER = 1 (MF = 11, 21), 30 + LIWP if MITER = 2, 3, 4, or 9. MAXL = 5 unless a different optional input value is given. LIWP = length of integer user work space for preconditioning (see conditional input list following). The first few words of IWORK are used for conditional and optional inputs and optional outputs.

  The following 4 words in IWORK are conditional inputs,
  required if MITER .ge. 1:
  IWORK(1) = LWP  = length of real array WP for use in
             preconditioning (part of RWORK array).
  IWORK(2) = LIWP = length of integer array IWP for use in
             preconditioning (part of IWORK array).
             The arrays WP and IWP are work arrays under the
             user's control, for use in the routines that
             perform preconditioning operations (JAC and PSOL).
  IWORK(3) = JPRE = preconditioner type flag:
           = 0 for no preconditioning (P1 = P2 = P = identity)
           = 1 for left-only preconditioning (P2 = identity)
           = 2 for right-only preconditioning (P1 = identity)
           = 3 for two-sided preconditioning (and PCG or PCGS)
  IWORK(4) = JACFLG = flag for whether JAC is called.
           = 0 if JAC is not to be called,
           = 1 if JAC is to be called.
             Use JACFLG = 1 if JAC computes any nonconstant
             data needed in preconditioning operations,
             such as some of the Jacobian elements.

LIW

the length of the array IWORK, as declared by the user. (This will be checked by the solver.)

Note: The work arrays must not be altered between calls to DLSODKR for the same problem, except possibly for the conditional and optional inputs, and except for the last 3*NEQ words of RWORK. The latter space is used for internal scratch space, and so is available for use by the user outside DLSODKR between calls, if desired (but not for use by any of the user-supplied routines).

JAC

the name of the user-supplied routine to compute any Jacobian elements (or approximations) involved in the matrix preconditioning operations (MITER .ge. 1). It is to have the form

        SUBROUTINE JAC (F, NEQ, T, Y, YSV, REWT, FTY, V, &
       &                HL0, JOK, WP, IWP, IER)
        DOUBLE PRECISION T, Y(*), YSV(*), REWT(*), FTY(*), V(*), &
       &                 HL0, WP(*)
        INTEGER IWP(*)

This routine must evaluate and preprocess any parts of the Jacobian matrix df/dy used in the preconditioners P1, P2, P.

The Y and FTY arrays contain the current values of y and f(t,y), respectively, and the YSV array also contains the current y vector. The array V is work space of length NEQ for use by JAC. REWT is the array of reciprocal error weights (1/EWT). JAC must multiply all computed Jacobian elements by the scalar -HL0, add the identity matrix, and do any factorization operations called for, in preparation for solving linear systems with a coefficient matrix of P1, P2, or P. The matrix P1*P2 or P should be an & approximation to identity - hl0 * (df/dy). JAC should return IER = 0 if successful, and IER .ne. 0 if not.

(If IER .ne. 0, a smaller time step will be tried.) The arrays WP (of length LWP) and IWP (of length LIWP) are for use by JAC and PSOL for work space and for storage of data needed for the solution of the preconditioner linear systems. Their lengths and contents are under the user’s control.

The argument JOK is an input flag for optional use by JAC in deciding whether to recompute Jacobian elements or use saved values. If JOK = -1, then JAC must compute any relevant Jacobian elements (or approximations) used in the preconditioners. Optionally, JAC may also save these elements for later reuse. If JOK = 1, the integrator has made a judgement (based on the convergence history and the value of HL0) that JAC need not recompute Jacobian elements, but instead use saved values, and the current value of HL0, to reconstruct the preconditioner matrices, followed by any required factorizations. This may be cost-effective if Jacobian elements are costly and storage is available.

JAC may alter Y and V, but not YSV, REWT, FTY, or HL0. JAC must be declared External in the calling program.

Subroutine JAC may access user-defined quantities in NEQ(2),… and/or in Y(NEQ(1)+1),… if NEQ is an array (dimensioned in JAC) and/or Y has length exceeding NEQ(1). See the descriptions of NEQ and Y above.

PSOL

the name of the user-supplied routine for the solution of preconditioner linear systems. It is to have the form

         SUBROUTINE PSOL (NEQ, T, Y, FTY, WK,HL0, WP,IWP, B, LR,IER)
         DOUBLE PRECISION T, Y(*), FTY(*), WK(*), HL0, WP(*), B(*)
         INTEGER IWP(*)

This routine must solve a linear system with B as right-hand side and one of the preconditioning matrices, P1, P2, or P, as coefficient matrix, and return the solution vector in B. LR is a flag concerning left vs right preconditioning, input to PSOL. PSOL is to use P1 if LR = 1 and P2 if LR = 2. In the case of the PCG or PCGS method, LR will be 3, and PSOL should solve the system P*x = B with the preconditioner P. In the case MITER = 9 (no Krylov iteration), LR will be 0, and PSOL is to return in B the desired approximate solution to A * x = B, where A = identity - hl0 * (df/dy). PSOL can use data generated in the JAC routine and stored in WP and IWP.

The Y and FTY arrays contain the current values of y and f(t,y), respectively. The array WK is work space of length NEQ for use by PSOL.

The argument HL0 is the current value of the scalar appearing in the linear system. If the old value, as of the last JAC call, is needed, it must have been saved by JAC in WP.

On return, PSOL should set the error flag IER as follows:

        IER = 0 if PSOL was successful,
        IER .gt. 0 on a recoverable error, meaning that the
               time step will be retried,
        IER .lt. 0 on an unrecoverable error, meaning that the
               solver is to stop immediately.
      PSOL may not alter Y, FTY, or HL0.
      PSOL must be declared External in the calling program.
           Subroutine PSOL may access user-defined quantities in
      NEQ(2),... and Y(NEQ(1)+1),... if NEQ is an array
      (dimensioned in PSOL) and/or Y has length exceeding NEQ(1).
      See the descriptions of NEQ and Y above.

MF

the method flag. Used only for input. The legal values of MF are 10, 11, 12, 13, 14, 19, 20, 21, 22, 23, 24, and 29. MF has decimal digits METH and MITER: MF = 10*METH + MITER. METH indicates the basic linear multistep method:

Differentiation Formulas (BDFs).

MITER indicates the corrector iteration method: MITER | description ----- | ---------------------------------------------------- 0 | means functional iteration (no linear system | is involved). 1 | means Newton iteration with Scaled Preconditioned | Incomplete Orthogonalization Method (SPIOM) | for the linear systems. 2 | means Newton iteration with Scaled Preconditioned | Incomplete Generalized Minimal Residual method | (SPIGMR) for the linear systems. 3 | means Newton iteration with Preconditioned | Conjugate Gradient method (PCG) | for the linear systems. 4 | means Newton iteration with scaled preconditioned | Conjugate Gradient method (PCGS) | for the linear systems. 9 | means Newton iteration with only the | user-supplied PSOL routine called (no Krylov | iteration) for the linear systems. | JPRE is ignored, and PSOL is called with LR = 0.

See comments in the introduction about the choice of MITER. If MITER .ge. 1, the user must supply routines JAC and PSOL (the names are arbitrary) as described above. For MITER = 0, a dummy argument can be used.

G

the name of subroutine for constraint functions, whose roots are desired during the integration. It is to have the form

          SUBROUTINE G (NEQ, T, Y, NG, GOUT)
          DOUBLE PRECISION T, Y(*), GOUT(NG)

where NEQ, T, Y, and NG are input, and the array GOUT is output. NEQ, T, and Y have the same meaning as in the F routine, and GOUT is an array of length NG. For i = 1,…,NG, this routine is to load into GOUT(i) the value at (t,y) of the i-th constraint function g(i). DLSODKR will find roots of the g(i) of odd multiplicity (i.e. sign changes) as they occur during the integration. G must be declared External in the calling program.

Caution: Because of numerical errors in the functions g(i) due to roundoff and integration error, DLSODKR may return false roots, or return the same root at two or more nearly equal values of t. If such false roots are suspected, the user should consider smaller error tolerances and/or higher precision in the evaluation of the g(i).

If a root of some g(i) defines the end of the problem, the input to DLSODKR should nevertheless allow integration to a point slightly past that root, so that DLSODKR can locate the root by interpolation.

Subroutine G may access user-defined quantities in NEQ(2),… and Y(NEQ(1)+1),… if NEQ is an array (dimensioned in G) and/or Y has length exceeding NEQ(1). See the descriptions of NEQ and Y above.

NG

number of constraint functions g(i). If there are none, set NG = 0, and pass a dummy name for G.

JROOT

integer array of length NG. Used only for output. On a return with ISTATE = 3 (one or more roots found), JROOT(i) = 1 if g(i) has a root at t, or JROOT(i) = 0 if not.

Optional Inputs.

The following is a list of the optional inputs provided for in the call sequence. (See also Part 2.) For each such input variable, this table lists its name as used in this documentation, its location in the call sequence, its meaning, and the default value. The use of any of these inputs requires IOPT = 1, and in that case all of these inputs are examined. A value of zero for any of these optional inputs will cause the default value to be used. Thus to use a subset of the optional inputs, simply preload locations 5 to 10 in RWORK and IWORK to 0.0 and 0 respectively, and then set those of interest to nonzero values.

Optional Outputs.

As optional additional output from DLSODKR, the variables listed below are quantities related to the performance of DLSODKR which are available to the user. These are communicated by way of the work arrays, but also have internal mnemonic names as shown.

Except where stated otherwise, all of these outputs are defined on any successful return from DLSODKR, and on any return with ISTATE = -1, -2, -4, -5, -6, or -7. On an illegal input return (ISTATE = -3), they will be unchanged from their existing values (if any), except possibly for TOLSF, LENRW, and LENIW.

On any error return, outputs relevant to the error will be defined, as noted below.

The following two arrays are segments of the RWORK array which may also be of interest to the user as optional outputs. For each array, the table below gives its internal name, its base address in RWORK, and its description.

Part 2. Other Routines Callable.

The following are optional calls which the user may make to gain additional capabilities in conjunction with DLSODKR. (The routines XSETUN and XSETF are designed to conform to the SLATEC error handling package.)

The detailed instructions for using DINTDY are as follows. The form of the call is:

   LYH = 21 + 3*NG
   CALL DINTDY (T, K, RWORK(LYH), NYH, DKY, IFLAG)

The input parameters are:

T

value of independent variable where answers are desired (normally the same as the T last returned by DLSODKR). For valid results, T must lie between TCUR - HU and TCUR. (See optional outputs for TCUR and HU.) K

integer order of the derivative desired. K must satisfy 0 .le. K .le. NQCUR, where NQCUR is the current order (see optional outputs). The capability corresponding to K = 0, i.e. computing y(T), is already provided by DLSODKR directly. Since NQCUR .ge. 1, the first derivative dy/dt is always available with DINTDY. LYH

21 + 3*NG = base address in RWORK of the history array YH. NYH

column length of YH, equal to the initial value of NEQ.

The output parameters are:

DKY

a real array of length NEQ containing the computed value of the K-th derivative of y(t).

IFLAG

integer flag, returned as 0 if K and T were legal, -1 if K was illegal, and -2 if T was illegal. On an error return, a message is also written.

Part 3. Common Blocks.

If DLSODKR is to be used in an overlay situation, the user must declare, in the primary overlay, the variables in:

   (1) the call sequence to DLSODKR, and
   (2) the four internal Common blocks
         /DLS001/  of length  255  (218 double precision words
                      followed by 37 integer words),
         /DLS002/  of length   5  (1 double precision word
                      followed by  4 integer words),
         /DLPK01/  of length  17  (4 double precision words
                      followed by 13 integer words),
         /DLSR01/  of length  14     (5 double precision words
                      followed by  9 integer words).

If DLSODKR is used on a system in which the contents of internal Common blocks are not preserved between calls, the user should declare the above Common blocks in the calling program to insure that their contents are preserved.

If the solution of a given problem by DLSODKR is to be interrupted and then later continued, such as when restarting an interrupted run or alternating between two or more problems, the user should save, following the return from the last DLSODKR call prior to the interruption, the contents of the call sequence variables and the internal Common blocks, and later restore these values before the next DLSODKR call for that problem. To save and restore the Common blocks, use Subroutine DSRCKR (see Part 2 above).

Part 4. Optionally Replaceable Solver Routines.

Below are descriptions of two routines in the DLSODKR package which relate to the measurement of errors. Either routine can be replaced by a user-supplied version, if desired. However, since such a replacement may have a major impact on performance, it should be done only when absolutely necessary, and only with great caution. (Note: The means by which the package version of a routine is superseded by the user’s version may be system-dependent.)

(a) DEWSET()

The following subroutine is called just before each internal integration step, and sets the array of error weights, EWT, as described under ITOL/RTOL/ATOL above:

     SUBROUTINE DEWSET (NEQ, ITOL, RTOL, ATOL, YCUR, EWT)

where NEQ, ITOL, RTOL, and ATOL are as in the DLSODKR call sequence, YCUR contains the current dependent variable vector, and EWT is the array of weights set by DEWSET.

If the user supplies this subroutine, it must return in EWT(i) (i = 1,…,NEQ) a positive quantity suitable for comparing errors in y(i) to. The EWT array returned by DEWSET is passed to the DVNORM routine (see below), and also used by DLSODKR in the computation of the optional output IMXER, the diagonal Jacobian approximation, and the increments for difference quotient Jacobians.

In the user-supplied version of DEWSET, it may be desirable to use the current values of derivatives of y. Derivatives up to order NQ are available from the history array YH, described above under optional outputs. In DEWSET, YH is identical to the YCUR array, extended to NQ + 1 columns with a column length of NYH and scale factors of H**j/factorial(j). On the first call for the problem, given by NST = 0, NQ is 1 and H is temporarily set to 1.0. NYH is the initial value of NEQ. The quantities NQ, H, and NST can be obtained by including in DEWSET the statements:

     DOUBLE PRECISION RLS
     COMMON /DLS001/ RLS(218),ILS(37)
     NQ = ILS(33)
     NST = ILS(34)
     H = RLS(212)

Thus, for example, the current value of dy/dt can be obtained as YCUR(NYH+i)/H (i=1,…,NEQ) (and the division by H is unnecessary when NST = 0).

(b) DVNORM() The following is a real function routine which computes the weighted root-mean-square norm of a vector v:

     D = DVNORM (N, V, W)
 where:
   N = the length of the vector,
   V = real array of length N containing the vector,
   W = real array of length N containing weights,
   D = SQRT( (1/N) * sum(V(i)*W(i))**2 ).

DVNORM is called with N = NEQ and with W(i) = 1.0/EWT(i), where EWT is as set by Subroutine DEWSET.

If the user supplies this function, it should return a non-negative value of DVNORM suitable for use in the error control in DLSODKR. None of the arguments should be altered by DVNORM. For example, a user-supplied DVNORM routine might:

• substitute a max-norm of (V(i)*W(i)) for the RMS-norm, or
• ignore some components of V in the norm, with the effect of suppressing the error control on those components of y.

Pedigree:

This is the 18 November 2003 version of DLSODKR is derived from the Livermore Solver for Ordinary Differential Equations package ODEPACK,

References: 1. Peter N. Brown and Alan C. Hindmarsh, Reduced Storage Matrix Methods in Stiff ODE Systems, J. Appl. Math. & Comp., 31 (1989), pp. 40-91; also L.L.N.L. Report UCRL-95088, Rev. 1, June 1987. 2. Alan C. Hindmarsh, ODEPACK, A Systematized Collection of ODE Solvers, in Scientific Computing, R. S. Stepleman et al. (Eds.), North-Holland, Amsterdam, 1983, pp. 55-64.

Authors: Alan C. Hindmarsh and Peter N. Brown Center for Applied Scientific Computing, L-561 Lawrence Livermore National Laboratory Livermore, CA 94551

Arguments