Math-Lsoda
view release on metacpan or search on metacpan
C END
C
C The output from this program (on a Cray-1 in single precision)
C is as follows.
C
C At t = 4.0000e-01 y = 9.851726e-01 3.386406e-05 1.479357e-02
C At t = 4.0000e+00 y = 9.055142e-01 2.240418e-05 9.446344e-02
C At t = 4.0000e+01 y = 7.158050e-01 9.184616e-06 2.841858e-01
C At t = 4.0000e+02 y = 4.504846e-01 3.222434e-06 5.495122e-01
C At t = 4.0000e+03 y = 1.831701e-01 8.940379e-07 8.168290e-01
C At t = 4.0000e+04 y = 3.897016e-02 1.621193e-07 9.610297e-01
C At t = 4.0000e+05 y = 4.935213e-03 1.983756e-08 9.950648e-01
C At t = 4.0000e+06 y = 5.159269e-04 2.064759e-09 9.994841e-01
C At t = 4.0000e+07 y = 5.306413e-05 2.122677e-10 9.999469e-01
C At t = 4.0000e+08 y = 5.494530e-06 2.197825e-11 9.999945e-01
C At t = 4.0000e+09 y = 5.129458e-07 2.051784e-12 9.999995e-01
C At t = 4.0000e+10 y = -7.170603e-08 -2.868241e-13 1.000000e+00
C
C No. steps = 330, No. f-s = 405, No. J-s = 69
C
C *Accuracy:
C The accuracy of the solution depends on the choice of tolerances
C RTOL and ATOL. Actual (global) errors may exceed these local
C tolerances, so choose them conservatively.
C
C *Cautions:
C The work arrays should not be altered between calls to DLSODE for
C the same problem, except possibly for the conditional and optional
C inputs.
C
C *Portability:
C Since NEQ is dimensioned inside DLSODE, some compilers may object
C to a call to DLSODE with NEQ a scalar variable. In this event,
C use DIMENSION NEQ(1). Similar remarks apply to RTOL and ATOL.
C
C Note to Cray users:
C For maximum efficiency, use the CFT77 compiler. Appropriate
C compiler optimization directives have been inserted for CFT77.
C
C *Reference:
C Alan C. Hindmarsh, "ODEPACK, A Systematized Collection of ODE
C Solvers," in Scientific Computing, R. S. Stepleman, et al., Eds.
C (North-Holland, Amsterdam, 1983), pp. 55-64.
C
C *Long Description:
C The following complete description of the user interface to
C DLSODE consists of four parts:
C
C 1. The call sequence to subroutine DLSODE, which is a driver
C routine for the solver. This includes descriptions of both
C the call sequence arguments and user-supplied routines.
C Following these descriptions is a description of optional
C inputs available through the call sequence, and then a
C description of optional outputs in the work arrays.
C
C 2. Descriptions of other routines in the DLSODE package that may
C be (optionally) called by the user. These provide the ability
C to alter error message handling, save and restore the internal
C COMMON, and obtain specified derivatives of the solution y(t).
C
C 3. Descriptions of COMMON block to be declared in overlay or
C similar environments, or to be saved when doing an interrupt
C of the problem and continued solution later.
C
C 4. Description of two routines in the DLSODE package, either of
C which the user may replace with his own version, if desired.
C These relate to the measurement of errors.
C
C
C Part 1. Call Sequence
C ----------------------
C
C Arguments
C ---------
C The call sequence parameters used for input only are
C
C F, NEQ, TOUT, ITOL, RTOL, ATOL, ITASK, IOPT, LRW, LIW, JAC, MF,
C
C and those used for both input and output are
C
C Y, T, ISTATE.
C
C The work arrays RWORK and IWORK are also used for conditional and
C optional inputs and optional outputs. (The term output here
C refers to the return from subroutine DLSODE to the user's calling
C program.)
C
C The legality of input parameters will be thoroughly checked on the
C initial call for the problem, but not checked thereafter unless a
C change in input parameters is flagged by ISTATE = 3 on input.
C
C The descriptions of the call arguments are as follows.
C
C F The name of the user-supplied subroutine defining the ODE
C system. The system must be put in the first-order form
C dy/dt = f(t,y), where f is a vector-valued function of
C the scalar t and the vector y. Subroutine F is to compute
C the function f. It is to have the form
C
C SUBROUTINE F (NEQ, T, Y, YDOT)
C DOUBLE PRECISION T, Y(*), YDOT(*)
C
C where NEQ, T, and Y are input, and the array YDOT =
C f(T,Y) is output. Y and YDOT are arrays of length NEQ.
C Subroutine F should not alter Y(1),...,Y(NEQ). F must be
C declared EXTERNAL in the calling program.
C
C Subroutine F may access user-defined quantities in
C NEQ(2),... and/or in Y(NEQ(1)+1),..., if NEQ is an array
C (dimensioned in F) and/or Y has length exceeding NEQ(1).
C See the descriptions of NEQ and Y below.
C
C If quantities computed in the F routine are needed
C externally to DLSODE, an extra call to F should be made
C for this purpose, for consistent and accurate results.
C If only the derivative dy/dt is needed, use DINTDY
C instead.
C
C NEQ The size of the ODE system (number of first-order
C ordinary differential equations). Used only for input.
C NEQ may be decreased, but not increased, during the
C immediately.
C CALL XSETF(MFLAG) Set a flag to control the printing of
C messages by DLSODE. MFLAG = 0 means do
C not print. (Danger: this risks losing
C valuable information.) MFLAG = 1 means
C print (the default). This call may be
C made at any time and will take effect
C immediately.
C CALL DSRCOM(RSAV,ISAV,JOB) Saves and restores the contents of the
C internal COMMON blocks used by DLSODE
C (see Part 3 below). RSAV must be a
C real array of length 218 or more, and
C ISAV must be an integer array of length
C 37 or more. JOB = 1 means save COMMON
C into RSAV/ISAV. JOB = 2 means restore
C COMMON from same. DSRCOM is useful if
C one is interrupting a run and restarting
C later, or alternating between two or
C more problems solved with DLSODE.
C CALL DINTDY(,,,,,) Provide derivatives of y, of various
C (see below) orders, at a specified point t, if
C desired. It may be called only after a
C successful return from DLSODE. Detailed
C instructions follow.
C
C Detailed instructions for using DINTDY
C --------------------------------------
C The form of the CALL is:
C
C CALL DINTDY (T, K, RWORK(21), NYH, DKY, IFLAG)
C
C The input parameters are:
C
C T Value of independent variable where answers are
C desired (normally the same as the T last returned by
C DLSODE). For valid results, T must lie between
C TCUR - HU and TCUR. (See "Optional Outputs" above
C for TCUR and HU.)
C K Integer order of the derivative desired. K must
C satisfy 0 <= K <= NQCUR, where NQCUR is the current
C order (see "Optional Outputs"). The capability
C corresponding to K = 0, i.e., computing y(t), is
C already provided by DLSODE directly. Since
C NQCUR >= 1, the first derivative dy/dt is always
C available with DINTDY.
C RWORK(21) The base address of the history array YH.
C NYH Column length of YH, equal to the initial value of NEQ.
C
C The output parameters are:
C
C DKY Real array of length NEQ containing the computed value
C of the Kth derivative of y(t).
C IFLAG Integer flag, returned as 0 if K and T were legal,
C -1 if K was illegal, and -2 if T was illegal.
C On an error return, a message is also written.
C
C
C Part 3. Common Blocks
C ----------------------
C
C If DLSODE is to be used in an overlay situation, the user must
C declare, in the primary overlay, the variables in:
C (1) the call sequence to DLSODE,
C (2) the internal COMMON block /DLS001/, of length 255
C (218 double precision words followed by 37 integer words).
C
C If DLSODE is used on a system in which the contents of internal
C COMMON blocks are not preserved between calls, the user should
C declare the above COMMON block in his main program to insure that
C its contents are preserved.
C
C If the solution of a given problem by DLSODE is to be interrupted
C and then later continued, as when restarting an interrupted run or
C alternating between two or more problems, the user should save,
C following the return from the last DLSODE call prior to the
C interruption, the contents of the call sequence variables and the
C internal COMMON block, and later restore these values before the
C next DLSODE call for that problem. In addition, if XSETUN and/or
C XSETF was called for non-default handling of error messages, then
C these calls must be repeated. To save and restore the COMMON
C block, use subroutine DSRCOM (see Part 2 above).
C
C
C Part 4. Optionally Replaceable Solver Routines
C -----------------------------------------------
C
C Below are descriptions of two routines in the DLSODE package which
C relate to the measurement of errors. Either routine can be
C replaced by a user-supplied version, if desired. However, since
C such a replacement may have a major impact on performance, it
C should be done only when absolutely necessary, and only with great
C caution. (Note: The means by which the package version of a
C routine is superseded by the user's version may be system-
C dependent.)
C
C DEWSET
C ------
C The following subroutine is called just before each internal
C integration step, and sets the array of error weights, EWT, as
C described under ITOL/RTOL/ATOL above:
C
C SUBROUTINE DEWSET (NEQ, ITOL, RTOL, ATOL, YCUR, EWT)
C
C where NEQ, ITOL, RTOL, and ATOL are as in the DLSODE call
C sequence, YCUR contains the current dependent variable vector,
C and EWT is the array of weights set by DEWSET.
C
C If the user supplies this subroutine, it must return in EWT(i)
C (i = 1,...,NEQ) a positive quantity suitable for comparing errors
C in Y(i) to. The EWT array returned by DEWSET is passed to the
C DVNORM routine (see below), and also used by DLSODE in the
C computation of the optional output IMXER, the diagonal Jacobian
C approximation, and the increments for difference quotient
C Jacobians.
C
C In the user-supplied version of DEWSET, it may be desirable to use
C the current values of derivatives of y. Derivatives up to order NQ
C are available from the history array YH, described above under
C optional outputs. In DEWSET, YH is identical to the YCUR array,
C extended to NQ + 1 columns with a column length of NYH and scale
C factors of H**j/factorial(j). On the first call for the problem,
C given by NST = 0, NQ is 1 and H is temporarily set to 1.0.
C PDJ(10) = -RK17*Y(10)
C PDJ(12) = -RK15*Y(2) - RK17*Y(10)
C RETURN
C END
C
C The output of this program (on a Cray-1 in single precision)
C is as follows:
C
C
C At t = 1.000e-01 No. steps = 12 Last step = 1.515e-02
C Y array = 9.90050e-01 6.28228e-03 3.65313e-03 7.51934e-07
C 1.12167e-09 1.18458e-09 1.77291e-12 3.26476e-07
C 5.46720e-08 9.99500e-06 4.48483e-08 2.76398e-06
C
C
C At t = 1.000e+00 No. steps = 33 Last step = 7.880e-02
C Y array = 9.04837e-01 9.13105e-03 8.20622e-02 2.49177e-05
C 1.85055e-06 1.96797e-06 1.46157e-07 2.39557e-05
C 3.26306e-05 7.21621e-04 5.06433e-05 3.05010e-03
C
C
C At t = 1.000e+01 No. steps = 48 Last step = 1.239e+00
C Y array = 3.67876e-01 3.68958e-03 3.65133e-01 4.48325e-05
C 6.10798e-05 4.33148e-05 5.90211e-05 1.18449e-04
C 3.15235e-03 3.56531e-03 4.15520e-03 2.48741e-01
C
C
C At t = 1.000e+02 No. steps = 91 Last step = 3.764e+00
C Y array = 4.44981e-05 4.42666e-07 4.47273e-04 -3.53257e-11
C 2.81577e-08 -9.67741e-11 2.77615e-07 1.45322e-07
C 1.56230e-02 4.37394e-06 1.60104e-02 9.52246e-01
C
C
C At t = 1.000e+03 No. steps = 111 Last step = 4.156e+02
C Y array = -2.65492e-13 2.60539e-14 -8.59563e-12 6.29355e-14
C -1.78066e-13 5.71471e-13 -1.47561e-12 4.58078e-15
C 1.56314e-02 1.37878e-13 1.60184e-02 9.52719e-01
C
C
C Required RWORK size = 442 IWORK size = 30
C No. steps = 111 No. f-s = 142 No. J-s = 2 No. LU-s = 20
C No. of nonzeros in J = 44 No. of nonzeros in LU = 50
C
C-----------------------------------------------------------------------
C Full Description of User Interface to DLSODES.
C
C The user interface to DLSODES consists of the following parts.
C
C 1. The call sequence to Subroutine DLSODES, which is a driver
C routine for the solver. This includes descriptions of both
C the call sequence arguments and of user-supplied routines.
C Following these descriptions is a description of
C optional inputs available through the call sequence, and then
C a description of optional outputs (in the work arrays).
C
C 2. Descriptions of other routines in the DLSODES package that may be
C (optionally) called by the user. These provide the ability to
C alter error message handling, save and restore the internal
C Common, and obtain specified derivatives of the solution y(t).
C
C 3. Descriptions of Common blocks to be declared in overlay
C or similar environments, or to be saved when doing an interrupt
C of the problem and continued solution later.
C
C 4. Description of two routines in the DLSODES package, either of
C which the user may replace with his/her own version, if desired.
C These relate to the measurement of errors.
C
C-----------------------------------------------------------------------
C Part 1. Call Sequence.
C
C The call sequence parameters used for input only are
C F, NEQ, TOUT, ITOL, RTOL, ATOL, ITASK, IOPT, LRW, LIW, JAC, MF,
C and those used for both input and output are
C Y, T, ISTATE.
C The work arrays RWORK and IWORK are also used for conditional and
C optional inputs and optional outputs. (The term output here refers
C to the return from Subroutine DLSODES to the user's calling program.)
C
C The legality of input parameters will be thoroughly checked on the
C initial call for the problem, but not checked thereafter unless a
C change in input parameters is flagged by ISTATE = 3 on input.
C
C The descriptions of the call arguments are as follows.
C
C F = the name of the user-supplied subroutine defining the
C ODE system. The system must be put in the first-order
C form dy/dt = f(t,y), where f is a vector-valued function
C of the scalar t and the vector y. Subroutine F is to
C compute the function f. It is to have the form
C SUBROUTINE F (NEQ, T, Y, YDOT)
C DOUBLE PRECISION T, Y(*), YDOT(*)
C where NEQ, T, and Y are input, and the array YDOT = f(t,y)
C is output. Y and YDOT are arrays of length NEQ.
C Subroutine F should not alter y(1),...,y(NEQ).
C F must be declared External in the calling program.
C
C Subroutine F may access user-defined quantities in
C NEQ(2),... and/or in Y(NEQ(1)+1),... if NEQ is an array
C (dimensioned in F) and/or Y has length exceeding NEQ(1).
C See the descriptions of NEQ and Y below.
C
C If quantities computed in the F routine are needed
C externally to DLSODES, an extra call to F should be made
C for this purpose, for consistent and accurate results.
C If only the derivative dy/dt is needed, use DINTDY instead.
C
C NEQ = the size of the ODE system (number of first order
C ordinary differential equations). Used only for input.
C NEQ may be decreased, but not increased, during the problem.
C If NEQ is decreased (with ISTATE = 3 on input), the
C remaining components of Y should be left undisturbed, if
C these are to be accessed in F and/or JAC.
C
C Normally, NEQ is a scalar, and it is generally referred to
C as a scalar in this user interface description. However,
C NEQ may be an array, with NEQ(1) set to the system size.
C (The DLSODES package accesses only NEQ(1).) In either case,
C this parameter is passed as the NEQ argument in all calls
C to F and JAC. Hence, if it is an array, locations
C NEQ(2),... may be used to store other integer data and pass
C
C CALL XSETF(MFLAG) Set a flag to control the printing of
C messages by DLSODES.
C MFLAG = 0 means do not print. (Danger:
C This risks losing valuable information.)
C MFLAG = 1 means print (the default).
C
C Either of the above calls may be made at
C any time and will take effect immediately.
C
C CALL DSRCMS(RSAV,ISAV,JOB) saves and restores the contents of
C the internal Common blocks used by
C DLSODES (see Part 3 below).
C RSAV must be a real array of length 224
C or more, and ISAV must be an integer
C array of length 71 or more.
C JOB=1 means save Common into RSAV/ISAV.
C JOB=2 means restore Common from RSAV/ISAV.
C DSRCMS is useful if one is
C interrupting a run and restarting
C later, or alternating between two or
C more problems solved with DLSODES.
C
C CALL DINTDY(,,,,,) Provide derivatives of y, of various
C (see below) orders, at a specified point t, if
C desired. It may be called only after
C a successful return from DLSODES.
C
C The detailed instructions for using DINTDY are as follows.
C The form of the call is:
C
C LYH = IWORK(22)
C CALL DINTDY (T, K, RWORK(LYH), NYH, DKY, IFLAG)
C
C The input parameters are:
C
C T = value of independent variable where answers are desired
C (normally the same as the T last returned by DLSODES).
C For valid results, T must lie between TCUR - HU and TCUR.
C (See optional outputs for TCUR and HU.)
C K = integer order of the derivative desired. K must satisfy
C 0 .le. K .le. NQCUR, where NQCUR is the current order
C (See optional outputs). The capability corresponding
C to K = 0, i.e. computing y(T), is already provided
C by DLSODES directly. Since NQCUR .ge. 1, the first
C derivative dy/dt is always available with DINTDY.
C LYH = the base address of the history array YH, obtained
C as an optional output as shown above.
C NYH = column length of YH, equal to the initial value of NEQ.
C
C The output parameters are:
C
C DKY = a real array of length NEQ containing the computed value
C of the K-th derivative of y(t).
C IFLAG = integer flag, returned as 0 if K and T were legal,
C -1 if K was illegal, and -2 if T was illegal.
C On an error return, a message is also written.
C-----------------------------------------------------------------------
C Part 3. Common Blocks.
C
C If DLSODES is to be used in an overlay situation, the user
C must declare, in the primary overlay, the variables in:
C (1) the call sequence to DLSODES, and
C (2) the two internal Common blocks
C /DLS001/ of length 255 (218 double precision words
C followed by 37 integer words),
C /DLSS01/ of length 40 (6 double precision words
C followed by 34 integer words),
C
C If DLSODES is used on a system in which the contents of internal
C Common blocks are not preserved between calls, the user should
C declare the above Common blocks in the calling program to insure
C that their contents are preserved.
C
C If the solution of a given problem by DLSODES is to be interrupted
C and then later continued, such as when restarting an interrupted run
C or alternating between two or more problems, the user should save,
C following the return from the last DLSODES call prior to the
C interruption, the contents of the call sequence variables and the
C internal Common blocks, and later restore these values before the
C next DLSODES call for that problem. To save and restore the Common
C blocks, use Subroutine DSRCMS (see Part 2 above).
C
C-----------------------------------------------------------------------
C Part 4. Optionally Replaceable Solver Routines.
C
C Below are descriptions of two routines in the DLSODES package which
C relate to the measurement of errors. Either routine can be
C replaced by a user-supplied version, if desired. However, since such
C a replacement may have a major impact on performance, it should be
C done only when absolutely necessary, and only with great caution.
C (Note: The means by which the package version of a routine is
C superseded by the user's version may be system-dependent.)
C
C (a) DEWSET.
C The following subroutine is called just before each internal
C integration step, and sets the array of error weights, EWT, as
C described under ITOL/RTOL/ATOL above:
C Subroutine DEWSET (NEQ, ITOL, RTOL, ATOL, YCUR, EWT)
C where NEQ, ITOL, RTOL, and ATOL are as in the DLSODES call sequence,
C YCUR contains the current dependent variable vector, and
C EWT is the array of weights set by DEWSET.
C
C If the user supplies this subroutine, it must return in EWT(i)
C (i = 1,...,NEQ) a positive quantity suitable for comparing errors
C in y(i) to. The EWT array returned by DEWSET is passed to the DVNORM
C routine (see below), and also used by DLSODES in the computation
C of the optional output IMXER, the diagonal Jacobian approximation,
C and the increments for difference quotient Jacobians.
C
C In the user-supplied version of DEWSET, it may be desirable to use
C the current values of derivatives of y. Derivatives up to order NQ
C are available from the history array YH, described above under
C optional outputs. In DEWSET, YH is identical to the YCUR array,
C extended to NQ + 1 columns with a column length of NYH and scale
C factors of H**j/factorial(j). On the first call for the problem,
C given by NST = 0, NQ is 1 and H is temporarily set to 1.0.
C NYH is the initial value of NEQ. The quantities NQ, H, and NST
C can be obtained by including in DEWSET the statements:
C DOUBLE PRECISION RLS
C COMMON /DLS001/ RLS(218),ILS(37)
C NQ = ILS(33)
C DO 40 IOUT = 1,12
C CALL DLSODA(FEX,NEQ,Y,T,TOUT,ITOL,RTOL,ATOL,ITASK,ISTATE,
C 1 IOPT,RWORK,LRW,IWORK,LIW,JDUM,JT)
C WRITE(6,20)T,Y(1),Y(2),Y(3)
C 20 FORMAT(' At t =',D12.4,' Y =',3D14.6)
C IF (ISTATE .LT. 0) GO TO 80
C 40 TOUT = TOUT*10.
C WRITE(6,60)IWORK(11),IWORK(12),IWORK(13),IWORK(19),RWORK(15)
C 60 FORMAT(/' No. steps =',I4,' No. f-s =',I4,' No. J-s =',I4/
C 1 ' Method last used =',I2,' Last switch was at t =',D12.4)
C STOP
C 80 WRITE(6,90)ISTATE
C 90 FORMAT(///' Error halt.. ISTATE =',I3)
C STOP
C END
C
C SUBROUTINE FEX (NEQ, T, Y, YDOT)
C DOUBLE PRECISION T, Y, YDOT
C DIMENSION Y(3), YDOT(3)
C YDOT(1) = -.04*Y(1) + 1.D4*Y(2)*Y(3)
C YDOT(3) = 3.D7*Y(2)*Y(2)
C YDOT(2) = -YDOT(1) - YDOT(3)
C RETURN
C END
C
C The output of this program (on a CDC-7600 in single precision)
C is as follows:
C
C At t = 4.0000e-01 y = 9.851712e-01 3.386380e-05 1.479493e-02
C At t = 4.0000e+00 Y = 9.055333e-01 2.240655e-05 9.444430e-02
C At t = 4.0000e+01 Y = 7.158403e-01 9.186334e-06 2.841505e-01
C At t = 4.0000e+02 Y = 4.505250e-01 3.222964e-06 5.494717e-01
C At t = 4.0000e+03 Y = 1.831975e-01 8.941774e-07 8.168016e-01
C At t = 4.0000e+04 Y = 3.898730e-02 1.621940e-07 9.610125e-01
C At t = 4.0000e+05 Y = 4.936363e-03 1.984221e-08 9.950636e-01
C At t = 4.0000e+06 Y = 5.161831e-04 2.065786e-09 9.994838e-01
C At t = 4.0000e+07 Y = 5.179817e-05 2.072032e-10 9.999482e-01
C At t = 4.0000e+08 Y = 5.283401e-06 2.113371e-11 9.999947e-01
C At t = 4.0000e+09 Y = 4.659031e-07 1.863613e-12 9.999995e-01
C At t = 4.0000e+10 Y = 1.404280e-08 5.617126e-14 1.000000e+00
C
C No. steps = 361 No. f-s = 693 No. J-s = 64
C Method last used = 2 Last switch was at t = 6.0092e-03
C-----------------------------------------------------------------------
C Full description of user interface to DLSODA.
C
C The user interface to DLSODA consists of the following parts.
C
C 1. The call sequence to Subroutine DLSODA, which is a driver
C routine for the solver. This includes descriptions of both
C the call sequence arguments and of user-supplied routines.
C following these descriptions is a description of
C optional inputs available through the call sequence, and then
C a description of optional outputs (in the work arrays).
C
C 2. Descriptions of other routines in the DLSODA package that may be
C (optionally) called by the user. These provide the ability to
C alter error message handling, save and restore the internal
C Common, and obtain specified derivatives of the solution y(t).
C
C 3. Descriptions of Common blocks to be declared in overlay
C or similar environments, or to be saved when doing an interrupt
C of the problem and continued solution later.
C
C 4. Description of a subroutine in the DLSODA package,
C which the user may replace with his/her own version, if desired.
C this relates to the measurement of errors.
C
C-----------------------------------------------------------------------
C Part 1. Call Sequence.
C
C The call sequence parameters used for input only are
C F, NEQ, TOUT, ITOL, RTOL, ATOL, ITASK, IOPT, LRW, LIW, JAC, JT,
C and those used for both input and output are
C Y, T, ISTATE.
C The work arrays RWORK and IWORK are also used for conditional and
C optional inputs and optional outputs. (The term output here refers
C to the return from Subroutine DLSODA to the user's calling program.)
C
C The legality of input parameters will be thoroughly checked on the
C initial call for the problem, but not checked thereafter unless a
C change in input parameters is flagged by ISTATE = 3 on input.
C
C The descriptions of the call arguments are as follows.
C
C F = the name of the user-supplied subroutine defining the
C ODE system. The system must be put in the first-order
C form dy/dt = f(t,y), where f is a vector-valued function
C of the scalar t and the vector y. Subroutine F is to
C compute the function f. It is to have the form
C SUBROUTINE F (NEQ, T, Y, YDOT)
C DOUBLE PRECISION T, Y(*), YDOT(*)
C where NEQ, T, and Y are input, and the array YDOT = f(t,y)
C is output. Y and YDOT are arrays of length NEQ.
C Subroutine F should not alter Y(1),...,Y(NEQ).
C F must be declared External in the calling program.
C
C Subroutine F may access user-defined quantities in
C NEQ(2),... and/or in Y(NEQ(1)+1),... if NEQ is an array
C (dimensioned in F) and/or Y has length exceeding NEQ(1).
C See the descriptions of NEQ and Y below.
C
C If quantities computed in the F routine are needed
C externally to DLSODA, an extra call to F should be made
C for this purpose, for consistent and accurate results.
C If only the derivative dy/dt is needed, use DINTDY instead.
C
C NEQ = the size of the ODE system (number of first order
C ordinary differential equations). Used only for input.
C NEQ may be decreased, but not increased, during the problem.
C If NEQ is decreased (with ISTATE = 3 on input), the
C remaining components of Y should be left undisturbed, if
C these are to be accessed in F and/or JAC.
C
C Normally, NEQ is a scalar, and it is generally referred to
C as a scalar in this user interface description. However,
C NEQ may be an array, with NEQ(1) set to the system size.
C (The DLSODA package accesses only NEQ(1).) In either case,
C this parameter is passed as the NEQ argument in all calls
C to F and JAC. Hence, if it is an array, locations
C NEQ(2),... may be used to store other integer data and pass
C the default is not desired.
C The default value of LUN is 6.
C
C CALL XSETF(MFLAG) set a flag to control the printing of
C messages by DLSODA.
C MFLAG = 0 means do not print. (Danger:
C This risks losing valuable information.)
C MFLAG = 1 means print (the default).
C
C Either of the above calls may be made at
C any time and will take effect immediately.
C
C CALL DSRCMA(RSAV,ISAV,JOB) saves and restores the contents of
C the internal Common blocks used by
C DLSODA (see Part 3 below).
C RSAV must be a real array of length 240
C or more, and ISAV must be an integer
C array of length 46 or more.
C JOB=1 means save Common into RSAV/ISAV.
C JOB=2 means restore Common from RSAV/ISAV.
C DSRCMA is useful if one is
C interrupting a run and restarting
C later, or alternating between two or
C more problems solved with DLSODA.
C
C CALL DINTDY(,,,,,) provide derivatives of y, of various
C (see below) orders, at a specified point t, if
C desired. It may be called only after
C a successful return from DLSODA.
C
C The detailed instructions for using DINTDY are as follows.
C The form of the call is:
C
C CALL DINTDY (T, K, RWORK(21), NYH, DKY, IFLAG)
C
C The input parameters are:
C
C T = value of independent variable where answers are desired
C (normally the same as the T last returned by DLSODA).
C For valid results, T must lie between TCUR - HU and TCUR.
C (See optional outputs for TCUR and HU.)
C K = integer order of the derivative desired. K must satisfy
C 0 .le. K .le. NQCUR, where NQCUR is the current order
C (see optional outputs). The capability corresponding
C to K = 0, i.e. computing y(T), is already provided
C by DLSODA directly. Since NQCUR .ge. 1, the first
C derivative dy/dt is always available with DINTDY.
C RWORK(21) = the base address of the history array YH.
C NYH = column length of YH, equal to the initial value of NEQ.
C
C The output parameters are:
C
C DKY = a real array of length NEQ containing the computed value
C of the K-th derivative of y(t).
C IFLAG = integer flag, returned as 0 if K and T were legal,
C -1 if K was illegal, and -2 if T was illegal.
C On an error return, a message is also written.
C-----------------------------------------------------------------------
C Part 3. Common Blocks.
C
C If DLSODA is to be used in an overlay situation, the user
C must declare, in the primary overlay, the variables in:
C (1) the call sequence to DLSODA, and
C (2) the two internal Common blocks
C /DLS001/ of length 255 (218 double precision words
C followed by 37 integer words),
C /DLSA01/ of length 31 (22 double precision words
C followed by 9 integer words).
C
C If DLSODA is used on a system in which the contents of internal
C Common blocks are not preserved between calls, the user should
C declare the above Common blocks in the calling program to insure
C that their contents are preserved.
C
C If the solution of a given problem by DLSODA is to be interrupted
C and then later continued, such as when restarting an interrupted run
C or alternating between two or more problems, the user should save,
C following the return from the last DLSODA call prior to the
C interruption, the contents of the call sequence variables and the
C internal Common blocks, and later restore these values before the
C next DLSODA call for that problem. To save and restore the Common
C blocks, use Subroutine DSRCMA (see Part 2 above).
C
C-----------------------------------------------------------------------
C Part 4. Optionally Replaceable Solver Routines.
C
C Below is a description of a routine in the DLSODA package which
C relates to the measurement of errors, and can be
C replaced by a user-supplied version, if desired. However, since such
C a replacement may have a major impact on performance, it should be
C done only when absolutely necessary, and only with great caution.
C (Note: The means by which the package version of a routine is
C superseded by the user's version may be system-dependent.)
C
C (a) DEWSET.
C The following subroutine is called just before each internal
C integration step, and sets the array of error weights, EWT, as
C described under ITOL/RTOL/ATOL above:
C Subroutine DEWSET (NEQ, ITOL, RTOL, ATOL, YCUR, EWT)
C where NEQ, ITOL, RTOL, and ATOL are as in the DLSODA call sequence,
C YCUR contains the current dependent variable vector, and
C EWT is the array of weights set by DEWSET.
C
C If the user supplies this subroutine, it must return in EWT(i)
C (i = 1,...,NEQ) a positive quantity suitable for comparing errors
C in y(i) to. The EWT array returned by DEWSET is passed to the
C DMNORM routine, and also used by DLSODA in the computation
C of the optional output IMXER, and the increments for difference
C quotient Jacobians.
C
C In the user-supplied version of DEWSET, it may be desirable to use
C the current values of derivatives of y. Derivatives up to order NQ
C are available from the history array YH, described above under
C optional outputs. In DEWSET, YH is identical to the YCUR array,
C extended to NQ + 1 columns with a column length of NYH and scale
C factors of H**j/factorial(j). On the first call for the problem,
C given by NST = 0, NQ is 1 and H is temporarily set to 1.0.
C NYH is the initial value of NEQ. The quantities NQ, H, and NST
C can be obtained by including in DEWSET the statements:
C DOUBLE PRECISION RLS
C COMMON /DLS001/ RLS(218),ILS(37)
C NQ = ILS(33)
C STOP
C END
C
C SUBROUTINE FEX (NEQ, T, Y, YDOT)
C DOUBLE PRECISION T, Y, YDOT
C DIMENSION Y(3), YDOT(3)
C YDOT(1) = -.04*Y(1) + 1.D4*Y(2)*Y(3)
C YDOT(3) = 3.D7*Y(2)*Y(2)
C YDOT(2) = -YDOT(1) - YDOT(3)
C RETURN
C END
C
C SUBROUTINE GEX (NEQ, T, Y, NG, GOUT)
C DOUBLE PRECISION T, Y, GOUT
C DIMENSION Y(3), GOUT(2)
C GOUT(1) = Y(1) - 1.D-4
C GOUT(2) = Y(3) - 1.D-2
C RETURN
C END
C
C The output of this program (on a CDC-7600 in single precision)
C is as follows:
C
C At t = 2.6400e-01 y = 9.899653e-01 3.470563e-05 1.000000e-02
C The above line is a root, JROOT = 0 1
C At t = 4.0000e-01 Y = 9.851712e-01 3.386380e-05 1.479493e-02
C At t = 4.0000e+00 Y = 9.055333e-01 2.240655e-05 9.444430e-02
C At t = 4.0000e+01 Y = 7.158403e-01 9.186334e-06 2.841505e-01
C At t = 4.0000e+02 Y = 4.505250e-01 3.222964e-06 5.494717e-01
C At t = 4.0000e+03 Y = 1.831975e-01 8.941774e-07 8.168016e-01
C At t = 4.0000e+04 Y = 3.898730e-02 1.621940e-07 9.610125e-01
C At t = 4.0000e+05 Y = 4.936363e-03 1.984221e-08 9.950636e-01
C At t = 4.0000e+06 Y = 5.161831e-04 2.065786e-09 9.994838e-01
C At t = 2.0745e+07 Y = 1.000000e-04 4.000395e-10 9.999000e-01
C The above line is a root, JROOT = 1 0
C At t = 4.0000e+07 Y = 5.179817e-05 2.072032e-10 9.999482e-01
C At t = 4.0000e+08 Y = 5.283401e-06 2.113371e-11 9.999947e-01
C At t = 4.0000e+09 Y = 4.659031e-07 1.863613e-12 9.999995e-01
C At t = 4.0000e+10 Y = 1.404280e-08 5.617126e-14 1.000000e+00
C
C No. steps = 361 No. f-s = 693 No. J-s = 64 No. g-s = 390
C Method last used = 2 Last switch was at t = 6.0092e-03
C
C-----------------------------------------------------------------------
C Full Description of User Interface to DLSODAR.
C
C The user interface to DLSODAR consists of the following parts.
C
C 1. The call sequence to Subroutine DLSODAR, which is a driver
C routine for the solver. This includes descriptions of both
C the call sequence arguments and of user-supplied routines.
C Following these descriptions is a description of
C optional inputs available through the call sequence, and then
C a description of optional outputs (in the work arrays).
C
C 2. Descriptions of other routines in the DLSODAR package that may be
C (optionally) called by the user. These provide the ability to
C alter error message handling, save and restore the internal
C Common, and obtain specified derivatives of the solution y(t).
C
C 3. Descriptions of Common blocks to be declared in overlay
C or similar environments, or to be saved when doing an interrupt
C of the problem and continued solution later.
C
C 4. Description of a subroutine in the DLSODAR package,
C which the user may replace with his/her own version, if desired.
C this relates to the measurement of errors.
C
C-----------------------------------------------------------------------
C Part 1. Call Sequence.
C
C The call sequence parameters used for input only are
C F, NEQ, TOUT, ITOL, RTOL, ATOL, ITASK, IOPT, LRW, LIW, JAC,
C JT, G, and NG,
C that used only for output is JROOT,
C and those used for both input and output are
C Y, T, ISTATE.
C The work arrays RWORK and IWORK are also used for conditional and
C optional inputs and optional outputs. (The term output here refers
C to the return from Subroutine DLSODAR to the user's calling program.)
C
C The legality of input parameters will be thoroughly checked on the
C initial call for the problem, but not checked thereafter unless a
C change in input parameters is flagged by ISTATE = 3 on input.
C
C The descriptions of the call arguments are as follows.
C
C F = the name of the user-supplied subroutine defining the
C ODE system. The system must be put in the first-order
C form dy/dt = f(t,y), where f is a vector-valued function
C of the scalar t and the vector y. Subroutine F is to
C compute the function f. It is to have the form
C SUBROUTINE F (NEQ, T, Y, YDOT)
C DOUBLE PRECISION T, Y(*), YDOT(*)
C where NEQ, T, and Y are input, and the array YDOT = f(t,y)
C is output. Y and YDOT are arrays of length NEQ.
C Subroutine F should not alter Y(1),...,Y(NEQ).
C F must be declared External in the calling program.
C
C Subroutine F may access user-defined quantities in
C NEQ(2),... and/or in Y(NEQ(1)+1),... if NEQ is an array
C (dimensioned in F) and/or Y has length exceeding NEQ(1).
C See the descriptions of NEQ and Y below.
C
C If quantities computed in the F routine are needed
C externally to DLSODAR, an extra call to F should be made
C for this purpose, for consistent and accurate results.
C If only the derivative dy/dt is needed, use DINTDY instead.
C
C NEQ = the size of the ODE system (number of first order
C ordinary differential equations). Used only for input.
C NEQ may be decreased, but not increased, during the problem.
C If NEQ is decreased (with ISTATE = 3 on input), the
C remaining components of Y should be left undisturbed, if
C these are to be accessed in F and/or JAC.
C
C Normally, NEQ is a scalar, and it is generally referred to
C as a scalar in this user interface description. However,
C NEQ may be an array, with NEQ(1) set to the system size.
C (The DLSODAR package accesses only NEQ(1).) In either case,
C this parameter is passed as the NEQ argument in all calls
C The default value of LUN is 6.
C
C CALL XSETF(MFLAG) Set a flag to control the printing of
C messages by DLSODAR.
C MFLAG = 0 means do not print. (Danger:
C This risks losing valuable information.)
C MFLAG = 1 means print (the default).
C
C Either of the above calls may be made at
C any time and will take effect immediately.
C
C CALL DSRCAR(RSAV,ISAV,JOB) saves and restores the contents of
C the internal Common blocks used by
C DLSODAR (see Part 3 below).
C RSAV must be a real array of length 245
C or more, and ISAV must be an integer
C array of length 55 or more.
C JOB=1 means save Common into RSAV/ISAV.
C JOB=2 means restore Common from RSAV/ISAV.
C DSRCAR is useful if one is
C interrupting a run and restarting
C later, or alternating between two or
C more problems solved with DLSODAR.
C
C CALL DINTDY(,,,,,) Provide derivatives of y, of various
C (see below) orders, at a specified point t, if
C desired. It may be called only after
C a successful return from DLSODAR.
C
C The detailed instructions for using DINTDY are as follows.
C The form of the call is:
C
C LYH = 21 + 3*NG
C CALL DINTDY (T, K, RWORK(LYH), NYH, DKY, IFLAG)
C
C The input parameters are:
C
C T = value of independent variable where answers are desired
C (normally the same as the T last returned by DLSODAR).
C For valid results, T must lie between TCUR - HU and TCUR.
C (See optional outputs for TCUR and HU.)
C K = integer order of the derivative desired. K must satisfy
C 0 .le. K .le. NQCUR, where NQCUR is the current order
C (see optional outputs). The capability corresponding
C to K = 0, i.e. computing y(t), is already provided
C by DLSODAR directly. Since NQCUR .ge. 1, the first
C derivative dy/dt is always available with DINTDY.
C LYH = 21 + 3*NG = base address in RWORK of the history array YH.
C NYH = column length of YH, equal to the initial value of NEQ.
C
C The output parameters are:
C
C DKY = a real array of length NEQ containing the computed value
C of the K-th derivative of y(t).
C IFLAG = integer flag, returned as 0 if K and T were legal,
C -1 if K was illegal, and -2 if T was illegal.
C On an error return, a message is also written.
C-----------------------------------------------------------------------
C Part 3. Common Blocks.
C
C If DLSODAR is to be used in an overlay situation, the user
C must declare, in the primary overlay, the variables in:
C (1) the call sequence to DLSODAR, and
C (2) the three internal Common blocks
C /DLS001/ of length 255 (218 double precision words
C followed by 37 integer words),
C /DLSA01/ of length 31 (22 double precision words
C followed by 9 integer words).
C /DLSR01/ of length 7 (3 double precision words
C followed by 4 integer words).
C
C If DLSODAR is used on a system in which the contents of internal
C Common blocks are not preserved between calls, the user should
C declare the above Common blocks in the calling program to insure
C that their contents are preserved.
C
C If the solution of a given problem by DLSODAR is to be interrupted
C and then later continued, such as when restarting an interrupted run
C or alternating between two or more problems, the user should save,
C following the return from the last DLSODAR call prior to the
C interruption, the contents of the call sequence variables and the
C internal Common blocks, and later restore these values before the
C next DLSODAR call for that problem. To save and restore the Common
C blocks, use Subroutine DSRCAR (see Part 2 above).
C
C-----------------------------------------------------------------------
C Part 4. Optionally Replaceable Solver Routines.
C
C Below is a description of a routine in the DLSODAR package which
C relates to the measurement of errors, and can be
C replaced by a user-supplied version, if desired. However, since such
C a replacement may have a major impact on performance, it should be
C done only when absolutely necessary, and only with great caution.
C (Note: The means by which the package version of a routine is
C superseded by the user's version may be system-dependent.)
C
C (a) DEWSET.
C The following subroutine is called just before each internal
C integration step, and sets the array of error weights, EWT, as
C described under ITOL/RTOL/ATOL above:
C Subroutine DEWSET (NEQ, ITOL, RTOL, ATOL, YCUR, EWT)
C where NEQ, ITOL, RTOL, and ATOL are as in the DLSODAR call sequence,
C YCUR contains the current dependent variable vector, and
C EWT is the array of weights set by DEWSET.
C
C If the user supplies this subroutine, it must return in EWT(i)
C (i = 1,...,NEQ) a positive quantity suitable for comparing errors
C in y(i) to. The EWT array returned by DEWSET is passed to the
C DMNORM routine, and also used by DLSODAR in the computation
C of the optional output IMXER, and the increments for difference
C quotient Jacobians.
C
C In the user-supplied version of DEWSET, it may be desirable to use
C the current values of derivatives of y. Derivatives up to order NQ
C are available from the history array YH, described above under
C optional outputs. In DEWSET, YH is identical to the YCUR array,
C extended to NQ + 1 columns with a column length of NYH and scale
C factors of H**j/factorial(j). On the first call for the problem,
C given by NST = 0, NQ is 1 and H is temporarily set to 1.0.
C NYH is the initial value of NEQ. The quantities NQ, H, and NST
C can be obtained by including in DEWSET the statements:
C DOUBLE PRECISION RLS
C RWORK = real work array of length at least:
C 20 + 16*NEQ for MF = 10,
C 45 + 17*NEQ + LWP for MF = 21,
C 61 + 17*NEQ + LWP for MF = 22,
C 20 + 15*NEQ + LWP for MF = 23 or 24,
C 20 + 12*NEQ + LWP for MF = 29.
C LRW = declared length of RWORK (in user's dimension).
C IWORK = integer work array of length at least:
C 30 for MF = 10,
C 35 + LIWP for MF = 21,
C 30 + LIWP for MF = 22, 23, 24, or 29.
C LIW = declared length of IWORK (in user's dimension).
C JAC,PSOL = names of subroutines for preconditioning.
C These names must be declared External in the calling program.
C MF = method flag. Standard values are:
C 10 for nonstiff (Adams) method.
C 21 for stiff (BDF) method, with preconditioned SIOM.
C 22 for stiff method, with preconditioned GMRES method.
C 23 for stiff method, with preconditioned CG method.
C 24 for stiff method, with scaled preconditioned CG method.
C 29 for stiff method, with user's PSOL routine only.
C Note that the main program must declare arrays Y, RWORK, IWORK,
C and possibly ATOL.
C
C E. The output from the first call (or any call) is:
C Y = array of computed values of y(t) vector.
C T = corresponding value of independent variable (normally TOUT).
C ISTATE = 2 if DLSODPK was successful, negative otherwise.
C -1 means excess work done on this call (perhaps wrong MF).
C -2 means excess accuracy requested (tolerances too small).
C -3 means illegal input detected (see printed message).
C -4 means repeated error test failures (check all inputs).
C -5 means repeated convergence failures (perhaps bad JAC
C or PSOL routine supplied or wrong choice of MF or
C tolerances, or this solver is inappropriate).
C -6 means error weight became zero during problem. (Solution
C component i vanished, and ATOL or ATOL(i) = 0.)
C -7 means an unrecoverable error occurred in PSOL.
C
C F. To continue the integration after a successful return, simply
C reset TOUT and call DLSODPK again. No other parameters need be reset.
C
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C Full Description of User Interface to DLSODPK.
C
C The user interface to DLSODPK consists of the following parts.
C
C 1. The call sequence to Subroutine DLSODPK, which is a driver
C routine for the solver. This includes descriptions of both
C the call sequence arguments and of user-supplied routines.
C Following these descriptions is a description of
C optional inputs available through the call sequence, and then
C a description of optional outputs (in the work arrays).
C
C 2. Descriptions of other routines in the DLSODPK package that may be
C (optionally) called by the user. These provide the ability to
C alter error message handling, save and restore the internal
C Common, and obtain specified derivatives of the solution y(t).
C
C 3. Descriptions of Common blocks to be declared in overlay
C or similar environments, or to be saved when doing an interrupt
C of the problem and continued solution later.
C
C 4. Description of two routines in the DLSODPK package, either of
C which the user may replace with his/her own version, if desired.
C These relate to the measurement of errors.
C
C-----------------------------------------------------------------------
C Part 1. Call Sequence.
C
C The call sequence parameters used for input only are
C F, NEQ, TOUT, ITOL, RTOL, ATOL, ITASK, IOPT, LRW, LIW, JAC, PSOL, MF,
C and those used for both input and output are
C Y, T, ISTATE.
C The work arrays RWORK and IWORK are also used for conditional and
C optional inputs and optional outputs. (The term output here refers
C to the return from Subroutine DLSODPK to the user's calling program.)
C
C The legality of input parameters will be thoroughly checked on the
C initial call for the problem, but not checked thereafter unless a
C change in input parameters is flagged by ISTATE = 3 on input.
C
C The descriptions of the call arguments are as follows.
C
C F = the name of the user-supplied subroutine defining the
C ODE system. The system must be put in the first-order
C form dy/dt = f(t,y), where f is a vector-valued function
C of the scalar t and the vector y. Subroutine F is to
C compute the function f. It is to have the form
C SUBROUTINE F (NEQ, T, Y, YDOT)
C DOUBLE PRECISION T, Y(*), YDOT(*)
C where NEQ, T, and Y are input, and the array YDOT = f(t,y)
C is output. Y and YDOT are arrays of length NEQ.
C Subroutine F should not alter Y(1),...,Y(NEQ).
C F must be declared External in the calling program.
C
C Subroutine F may access user-defined quantities in
C NEQ(2),... and/or in Y(NEQ(1)+1),... if NEQ is an array
C (dimensioned in F) and/or Y has length exceeding NEQ(1).
C See the descriptions of NEQ and Y below.
C
C If quantities computed in the F routine are needed
C externally to DLSODPK, an extra call to F should be made
C for this purpose, for consistent and accurate results.
C If only the derivative dy/dt is needed, use DINTDY instead.
C
C NEQ = the size of the ODE system (number of first order
C ordinary differential equations). Used only for input.
C NEQ may be decreased, but not increased, during the problem.
C If NEQ is decreased (with ISTATE = 3 on input), the
C remaining components of Y should be left undisturbed, if
C these are to be accessed in the user-supplied subroutines.
C
C Normally, NEQ is a scalar, and it is generally referred to
C as a scalar in this user interface description. However,
C NEQ may be an array, with NEQ(1) set to the system size.
C (The DLSODPK package accesses only NEQ(1).) In either case,
C this parameter is passed as the NEQ argument in all calls
C to F, JAC, and PSOL. Hence, if it is an array, locations
C NEQ(2),... may be used to store other integer data and pass
C the default is not desired.
C The default value of lun is 6.
C
C CALL XSETF(MFLAG) Set a flag to control the printing of
C messages by DLSODPK.
C MFLAG = 0 means do not print. (Danger:
C This risks losing valuable information.)
C MFLAG = 1 means print (the default).
C
C Either of the above calls may be made at
C any time and will take effect immediately.
C
C CALL DSRCPK(RSAV,ISAV,JOB) saves and restores the contents of
C the internal Common blocks used by
C DLSODPK (see Part 3 below).
C RSAV must be a real array of length 222
C or more, and ISAV must be an integer
C array of length 50 or more.
C JOB=1 means save Common into RSAV/ISAV.
C JOB=2 means restore Common from RSAV/ISAV.
C DSRCPK is useful if one is
C interrupting a run and restarting
C later, or alternating between two or
C more problems solved with DLSODPK.
C
C CALL DINTDY(,,,,,) Provide derivatives of y, of various
C (See below) orders, at a specified point t, if
C desired. It may be called only after
C a successful return from DLSODPK.
C
C The detailed instructions for using DINTDY are as follows.
C The form of the call is:
C
C CALL DINTDY (T, K, RWORK(21), NYH, DKY, IFLAG)
C
C The input parameters are:
C
C T = value of independent variable where answers are desired
C (normally the same as the T last returned by DLSODPK).
C for valid results, T must lie between TCUR - HU and TCUR.
C (See optional outputs for TCUR and HU.)
C K = integer order of the derivative desired. K must satisfy
C 0 .le. K .le. NQCUR, where NQCUR is the current order
C (see optional outputs). The capability corresponding
C to K = 0, i.e. computing y(T), is already provided
C by DLSODPK directly. Since NQCUR .ge. 1, the first
C derivative dy/dt is always available with DINTDY.
C RWORK(21) = the base address of the history array YH.
C NYH = column length of YH, equal to the initial value of NEQ.
C
C The output parameters are:
C
C DKY = a real array of length NEQ containing the computed value
C of the K-th derivative of y(t).
C IFLAG = integer flag, returned as 0 if K and T were legal,
C -1 if K was illegal, and -2 if T was illegal.
C On an error return, a message is also written.
C-----------------------------------------------------------------------
C Part 3. Common Blocks.
C
C If DLSODPK is to be used in an overlay situation, the user
C must declare, in the primary overlay, the variables in:
C (1) the call sequence to DLSODPK, and
C (2) the two internal Common blocks
C /DLS001/ of length 255 (218 double precision words
C followed by 37 integer words),
C /DLPK01/ of length 17 (4 double precision words
C followed by 13 integer words).
C
C If DLSODPK is used on a system in which the contents of internal
C Common blocks are not preserved between calls, the user should
C declare the above Common blocks in the calling program to insure
C that their contents are preserved.
C
C If the solution of a given problem by DLSODPK is to be interrupted
C and then later continued, such as when restarting an interrupted run
C or alternating between two or more problems, the user should save,
C following the return from the last DLSODPK call prior to the
C interruption, the contents of the call sequence variables and the
C internal Common blocks, and later restore these values before the
C next DLSODPK call for that problem. To save and restore the Common
C blocks, use Subroutine DSRCPK (see Part 2 above).
C
C-----------------------------------------------------------------------
C Part 4. Optionally Replaceable Solver Routines.
C
C below are descriptions of two routines in the DLSODPK package which
C relate to the measurement of errors. Either routine can be
C replaced by a user-supplied version, if desired. However, since such
C a replacement may have a major impact on performance, it should be
C done only when absolutely necessary, and only with great caution.
C (Note: The means by which the package version of a routine is
C superseded by the user's version may be system-dependent.)
C
C (a) DEWSET.
C The following subroutine is called just before each internal
C integration step, and sets the array of error weights, EWT, as
C described under ITOL/RTOL/ATOL above:
C SUBROUTINE DEWSET (NEQ, ITOL, RTOL, ATOL, YCUR, EWT)
C where NEQ, ITOL, RTOL, and ATOL are as in the DLSODPK call sequence,
C YCUR contains the current dependent variable vector, and
C EWT is the array of weights set by DEWSET.
C
C If the user supplies this subroutine, it must return in EWT(i)
C (i = 1,...,NEQ) a positive quantity suitable for comparing errors
C in y(i) to. The EWT array returned by DEWSET is passed to the DVNORM
C routine (see below), and also used by DLSODPK in the computation
C of the optional output IMXER, the diagonal Jacobian approximation,
C and the increments for difference quotient Jacobians.
C
C In the user-supplied version of DEWSET, it may be desirable to use
C the current values of derivatives of y. Derivatives up to order NQ
C are available from the history array YH, described above under
C optional outputs. In DEWSET, YH is identical to the YCUR array,
C extended to NQ + 1 columns with a column length of NYH and scale
C factors of H**j/factorial(j). On the first call for the problem,
C given by NST = 0, NQ is 1 and H is temporarily set to 1.0.
C NYH is the initial value of NEQ. The quantities NQ, H, and NST
C can be obtained by including in DEWSET the statements:
C DOUBLE PRECISION RLS
C COMMON /DLS001/ RLS(218),ILS(37)
C NQ = ILS(33)
C MF = method flag. Standard values are:
C 10 for nonstiff (Adams) method.
C 21 for stiff (BDF) method, with preconditioned SIOM.
C 22 for stiff method, with preconditioned GMRES method.
C 23 for stiff method, with preconditioned CG method.
C 24 for stiff method, with scaled preconditioned CG method.
C 29 for stiff method, with user's PSOL routine only.
C G = name of subroutine for constraint functions, whose
C roots are desired during the integration.
C This name must be declared External in calling program.
C NG = number of constraint functions g(i). If there are none,
C set NG = 0, and pass a dummy name for G.
C JROOT = integer array of length NG for output of root information.
C See next paragraph.
C Note that the main program must declare arrays Y, RWORK, IWORK,
C JROOT, and possibly ATOL.
C
C F. The output from the first call (or any call) is:
C Y = array of computed values of y(t) vector.
C T = corresponding value of independent variable (normally TOUT).
C ISTATE = 2 or 3 if DLSODKR was successful, negative otherwise.
C 2 means no root was found, and TOUT was reached as desired.
C 3 means a root was found prior to reaching TOUT.
C -1 means excess work done on this call (perhaps wrong MF).
C -2 means excess accuracy requested (tolerances too small).
C -3 means illegal input detected (see printed message).
C -4 means repeated error test failures (check all inputs).
C -5 means repeated convergence failures (perhaps bad JAC
C or PSOL routine supplied or wrong choice of MF or
C tolerances, or this solver is inappropriate).
C -6 means error weight became zero during problem. (Solution
C component i vanished, and ATOL or ATOL(i) = 0.)
C -7 means an unrecoverable error occurred in PSOL.
C JROOT = array showing roots found if ISTATE = 3 on return.
C JROOT(i) = 1 if g(i) has a root at T, or 0 otherwise.
C
C G. To continue the integration after a successful return, proceed
C as follows:
C (a) If ISTATE = 2 on return, reset TOUT and call DLSODKR again.
C (b) If ISTATE = 3 on return, reset ISTATE to 2 and call DLSODKR again.
C In either case, no other parameters need be reset.
C
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C Full Description of User Interface to DLSODKR.
C
C The user interface to DLSODKR consists of the following parts.
C
C 1. The call sequence to Subroutine DLSODKR, which is a driver
C routine for the solver. This includes descriptions of both
C the call sequence arguments and of user-supplied routines.
C Following these descriptions is a description of
C optional inputs available through the call sequence, and then
C a description of optional outputs (in the work arrays).
C
C 2. Descriptions of other routines in the DLSODKR package that may be
C (optionally) called by the user. These provide the ability to
C alter error message handling, save and restore the internal
C Common, and obtain specified derivatives of the solution y(t).
C
C 3. Descriptions of Common blocks to be declared in overlay
C or similar environments, or to be saved when doing an interrupt
C of the problem and continued solution later.
C
C 4. Description of two routines in the DLSODKR package, either of
C which the user may replace with his/her own version, if desired.
C These relate to the measurement of errors.
C
C-----------------------------------------------------------------------
C Part 1. Call Sequence.
C
C The call sequence parameters used for input only are
C F, NEQ, TOUT, ITOL, RTOL, ATOL, ITASK, IOPT, LRW, LIW, JAC, PSOL, MF,
C G, and NG,
C that used only for output is JROOT,
C and those used for both input and output are
C Y, T, ISTATE.
C The work arrays RWORK and IWORK are also used for conditional and
C optional inputs and optional outputs. (The term output here refers
C to the return from Subroutine DLSODKR to the user's calling program.)
C
C The legality of input parameters will be thoroughly checked on the
C initial call for the problem, but not checked thereafter unless a
C change in input parameters is flagged by ISTATE = 3 on input.
C
C The descriptions of the call arguments are as follows.
C
C F = the name of the user-supplied subroutine defining the
C ODE system. The system must be put in the first-order
C form dy/dt = f(t,y), where f is a vector-valued function
C of the scalar t and the vector y. Subroutine F is to
C compute the function f. It is to have the form
C SUBROUTINE F (NEQ, T, Y, YDOT)
C DOUBLE PRECISION T, Y(*), YDOT(*)
C where NEQ, T, and Y are input, and the array YDOT = f(t,y)
C is output. Y and YDOT are arrays of length NEQ.
C Subroutine F should not alter Y(1),...,Y(NEQ).
C F must be declared External in the calling program.
C
C Subroutine F may access user-defined quantities in
C NEQ(2),... and/or in Y(NEQ(1)+1),... if NEQ is an array
C (dimensioned in F) and/or Y has length exceeding NEQ(1).
C See the descriptions of NEQ and Y below.
C
C If quantities computed in the F routine are needed
C externally to DLSODKR, an extra call to F should be made
C for this purpose, for consistent and accurate results.
C If only the derivative dy/dt is needed, use DINTDY instead.
C
C NEQ = the size of the ODE system (number of first order
C ordinary differential equations). Used only for input.
C NEQ may be decreased, but not increased, during the problem.
C If NEQ is decreased (with ISTATE = 3 on input), the
C remaining components of Y should be left undisturbed, if
C these are to be accessed in the user-supplied routines.
C
C Normally, NEQ is a scalar, and it is generally referred to
C as a scalar in this user interface description. However,
C NEQ may be an array, with NEQ(1) set to the system size.
C (The DLSODKR package accesses only NEQ(1).) In either case,
C this parameter is passed as the NEQ argument in all calls
C The default value of LUN is 6.
C
C CALL XSETF(MFLAG) Set a flag to control the printing of
C messages by DLSODKR.
C MFLAG = 0 means do not print. (Danger:
C This risks losing valuable information.)
C MFLAG = 1 means print (the default).
C
C Either of the above calls may be made at
C any time and will take effect immediately.
C
C CALL DSRCKR(RSAV,ISAV,JOB) saves and restores the contents of
C the internal Common blocks used by
C DLSODKR (see Part 3 below).
C RSAV must be a real array of length 228
C or more, and ISAV must be an integer
C array of length 63 or more.
C JOB=1 means save Common into RSAV/ISAV.
C JOB=2 means restore Common from RSAV/ISAV.
C DSRCKR is useful if one is
C interrupting a run and restarting
C later, or alternating between two or
C more problems solved with DLSODKR.
C
C CALL DINTDY(,,,,,) Provide derivatives of y, of various
C (see below) orders, at a specified point t, if
C desired. It may be called only after
C a successful return from DLSODKR.
C
C The detailed instructions for using DINTDY are as follows.
C The form of the call is:
C
C LYH = 21 + 3*NG
C CALL DINTDY (T, K, RWORK(LYH), NYH, DKY, IFLAG)
C
C The input parameters are:
C
C T = value of independent variable where answers are desired
C (normally the same as the T last returned by DLSODKR).
C For valid results, T must lie between TCUR - HU and TCUR.
C (See optional outputs for TCUR and HU.)
C K = integer order of the derivative desired. K must satisfy
C 0 .le. K .le. NQCUR, where NQCUR is the current order
C (see optional outputs). The capability corresponding
C to K = 0, i.e. computing y(T), is already provided
C by DLSODKR directly. Since NQCUR .ge. 1, the first
C derivative dy/dt is always available with DINTDY.
C LYH = 21 + 3*NG = base address in RWORK of the history array YH.
C NYH = column length of YH, equal to the initial value of NEQ.
C
C The output parameters are:
C
C DKY = a real array of length NEQ containing the computed value
C of the K-th derivative of y(t).
C IFLAG = integer flag, returned as 0 if K and T were legal,
C -1 if K was illegal, and -2 if T was illegal.
C On an error return, a message is also written.
C-----------------------------------------------------------------------
C Part 3. Common Blocks.
C
C If DLSODKR is to be used in an overlay situation, the user
C must declare, in the primary overlay, the variables in:
C (1) the call sequence to DLSODKR, and
C (2) the four internal Common blocks
C /DLS001/ of length 255 (218 double precision words
C followed by 37 integer words),
C /DLS002/ of length 5 (1 double precision word
C followed by 4 integer words),
C /DLPK01/ of length 17 (4 double precision words
C followed by 13 integer words),
C /DLSR01/ of length 14 (5 double precision words
C followed by 9 integer words).
C
C If DLSODKR is used on a system in which the contents of internal
C Common blocks are not preserved between calls, the user should
C declare the above Common blocks in the calling program to insure
C that their contents are preserved.
C
C If the solution of a given problem by DLSODKR is to be interrupted
C and then later continued, such as when restarting an interrupted run
C or alternating between two or more problems, the user should save,
C following the return from the last DLSODKR call prior to the
C interruption, the contents of the call sequence variables and the
C internal Common blocks, and later restore these values before the
C next DLSODKR call for that problem. To save and restore the Common
C blocks, use Subroutine DSRCKR (see Part 2 above).
C
C-----------------------------------------------------------------------
C Part 4. Optionally Replaceable Solver Routines.
C
C Below are descriptions of two routines in the DLSODKR package which
C relate to the measurement of errors. Either routine can be
C replaced by a user-supplied version, if desired. However, since such
C a replacement may have a major impact on performance, it should be
C done only when absolutely necessary, and only with great caution.
C (Note: The means by which the package version of a routine is
C superseded by the user's version may be system-dependent.)
C
C (a) DEWSET.
C The following subroutine is called just before each internal
C integration step, and sets the array of error weights, EWT, as
C described under ITOL/RTOL/ATOL above:
C SUBROUTINE DEWSET (NEQ, ITOL, RTOL, ATOL, YCUR, EWT)
C where NEQ, ITOL, RTOL, and ATOL are as in the DLSODKR call sequence,
C YCUR contains the current dependent variable vector, and
C EWT is the array of weights set by DEWSET.
C
C If the user supplies this subroutine, it must return in EWT(i)
C (i = 1,...,NEQ) a positive quantity suitable for comparing errors
C in y(i) to. The EWT array returned by DEWSET is passed to the DVNORM
C routine (see below), and also used by DLSODKR in the computation
C of the optional output IMXER, the diagonal Jacobian approximation,
C and the increments for difference quotient Jacobians.
C
C In the user-supplied version of DEWSET, it may be desirable to use
C the current values of derivatives of y. Derivatives up to order NQ
C are available from the history array YH, described above under
C optional outputs. In DEWSET, YH is identical to the YCUR array,
C extended to NQ + 1 columns with a column length of NYH and scale
C factors of H**j/factorial(j). On the first call for the problem,
C given by NST = 0, NQ is 1 and H is temporarily set to 1.0.
C NYH is the initial value of NEQ. The quantities NQ, H, and NST
C END
C
C SUBROUTINE APLUSP(NEQ, T, Y, ML, MU, P, NROWP)
C DOUBLE PRECISION T, Y, P
C DIMENSION Y(3), P(NROWP,3)
C P(1,1) = P(1,1) + 1.
C P(2,2) = P(2,2) + 1.
C RETURN
C END
C
C SUBROUTINE DGBYDY(NEQ, T, Y, S, ML, MU, P, NROWP)
C DOUBLE PRECISION T, Y, S, P
C DIMENSION Y(3), S(3), P(NROWP,3)
C P(1,1) = -.04
C P(1,2) = 1.D4*Y(3)
C P(1,3) = 1.D4*Y(2)
C P(2,1) = .04
C P(2,2) = -1.D4*Y(3) - 6.D7*Y(2)
C P(2,3) = -1.D4*Y(2)
C P(3,1) = 1.
C P(3,2) = 1.
C P(3,3) = 1.
C RETURN
C END
C
C The output of this program (on a CDC-7600 in single precision)
C is as follows:
C
C At t = 4.0000e-01 Y = 9.851726e-01 3.386406e-05 1.479357e-02
C At t = 4.0000e+00 Y = 9.055142e-01 2.240418e-05 9.446344e-02
C At t = 4.0000e+01 Y = 7.158050e-01 9.184616e-06 2.841858e-01
C At t = 4.0000e+02 Y = 4.504846e-01 3.222434e-06 5.495122e-01
C At t = 4.0000e+03 Y = 1.831701e-01 8.940379e-07 8.168290e-01
C At t = 4.0000e+04 Y = 3.897016e-02 1.621193e-07 9.610297e-01
C At t = 4.0000e+05 Y = 4.935213e-03 1.983756e-08 9.950648e-01
C At t = 4.0000e+06 Y = 5.159269e-04 2.064759e-09 9.994841e-01
C At t = 4.0000e+07 Y = 5.306413e-05 2.122677e-10 9.999469e-01
C At t = 4.0000e+08 Y = 5.494532e-06 2.197826e-11 9.999945e-01
C At t = 4.0000e+09 Y = 5.129457e-07 2.051784e-12 9.999995e-01
C At t = 4.0000e+10 Y = -7.170472e-08 -2.868188e-13 1.000000e+00
C
C No. steps = 330 No. r-s = 404 No. J-s = 69
C
C-----------------------------------------------------------------------
C Full Description of User Interface to DLSODI.
C
C The user interface to DLSODI consists of the following parts.
C
C 1. The call sequence to Subroutine DLSODI, which is a driver
C routine for the solver. This includes descriptions of both
C the call sequence arguments and of user-supplied routines.
C Following these descriptions is a description of
C optional inputs available through the call sequence, and then
C a description of optional outputs (in the work arrays).
C
C 2. Descriptions of other routines in the DLSODI package that may be
C (optionally) called by the user. These provide the ability to
C alter error message handling, save and restore the internal
C Common, and obtain specified derivatives of the solution y(t).
C
C 3. Descriptions of Common blocks to be declared in overlay
C or similar environments, or to be saved when doing an interrupt
C of the problem and continued solution later.
C
C 4. Description of two routines in the DLSODI package, either of
C which the user may replace with his/her own version, if desired.
C These relate to the measurement of errors.
C
C-----------------------------------------------------------------------
C Part 1. Call Sequence.
C
C The call sequence parameters used for input only are
C RES, ADDA, JAC, NEQ, TOUT, ITOL, RTOL, ATOL, ITASK,
C IOPT, LRW, LIW, MF,
C and those used for both input and output are
C Y, T, ISTATE, YDOTI.
C The work arrays RWORK and IWORK are also used for conditional and
C optional inputs and optional outputs. (The term output here refers
C to the return from Subroutine DLSODI to the user's calling program.)
C
C The legality of input parameters will be thoroughly checked on the
C initial call for the problem, but not checked thereafter unless a
C change in input parameters is flagged by ISTATE = 3 on input.
C
C The descriptions of the call arguments are as follows.
C
C RES = the name of the user-supplied subroutine which supplies
C the residual vector for the ODE system, defined by
C r = g(t,y) - A(t,y) * s
C as a function of the scalar t and the vectors
C s and y (s approximates dy/dt). This subroutine
C is to have the form
C SUBROUTINE RES (NEQ, T, Y, S, R, IRES)
C DOUBLE PRECISION T, Y(*), S(*), R(*)
C where NEQ, T, Y, S, and IRES are input, and R and
C IRES are output. Y, S, and R are arrays of length NEQ.
C On input, IRES indicates how DLSODI will use the
C returned array R, as follows:
C IRES = 1 means that DLSODI needs the full residual,
C r = g - A*s, exactly.
C IRES = -1 means that DLSODI is using R only to compute
C the Jacobian dr/dy by difference quotients.
C The RES routine can ignore IRES, or it can omit some terms
C if IRES = -1. If A does not depend on y, then RES can
C just return R = g when IRES = -1. If g - A*s contains other
C additive terms that are independent of y, these can also be
C dropped, if done consistently, when IRES = -1.
C The subroutine should set the flag IRES if it
C encounters a halt condition or illegal input.
C Otherwise, it should not reset IRES. On output,
C IRES = 1 or -1 represents a normal return, and
C DLSODI continues integrating the ODE. Leave IRES
C unchanged from its input value.
C IRES = 2 tells DLSODI to immediately return control
C to the calling program, with ISTATE = 3. This lets
C the calling program change parameters of the problem,
C if necessary.
C IRES = 3 represents an error condition (for example, an
C illegal value of y). DLSODI tries to integrate the system
C without getting IRES = 3 from RES. If it cannot, DLSODI
C returns with ISTATE = -7 or -1.
C the default is not desired.
C The default value of LUN is 6.
C
C CALL XSETF(MFLAG) Set a flag to control the printing of
C messages by DLSODI.
C MFLAG = 0 means do not print. (Danger:
C This risks losing valuable information.)
C MFLAG = 1 means print (the default).
C
C Either of the above calls may be made at
C any time and will take effect immediately.
C
C CALL DSRCOM(RSAV,ISAV,JOB) saves and restores the contents of
C the internal Common blocks used by
C DLSODI (see Part 3 below).
C RSAV must be a real array of length 218
C or more, and ISAV must be an integer
C array of length 37 or more.
C JOB=1 means save Common into RSAV/ISAV.
C JOB=2 means restore Common from RSAV/ISAV.
C DSRCOM is useful if one is
C interrupting a run and restarting
C later, or alternating between two or
C more problems solved with DLSODI.
C
C CALL DINTDY(,,,,,) Provide derivatives of y, of various
C (see below) orders, at a specified point t, if
C desired. It may be called only after
C a successful return from DLSODI.
C
C The detailed instructions for using DINTDY are as follows.
C The form of the call is:
C
C CALL DINTDY (T, K, RWORK(21), NYH, DKY, IFLAG)
C
C The input parameters are:
C
C T = value of independent variable where answers are desired
C (normally the same as the T last returned by DLSODI).
C For valid results, T must lie between TCUR - HU and TCUR.
C (See optional outputs for TCUR and HU.)
C K = integer order of the derivative desired. K must satisfy
C 0 .le. K .le. NQCUR, where NQCUR is the current order
C (see optional outputs). The capability corresponding
C to K = 0, i.e. computing y(T), is already provided
C by DLSODI directly. Since NQCUR .ge. 1, the first
C derivative dy/dt is always available with DINTDY.
C RWORK(21) = the base address of the history array YH.
C NYH = column length of YH, equal to the initial value of NEQ.
C
C The output parameters are:
C
C DKY = a real array of length NEQ containing the computed value
C of the K-th derivative of y(t).
C IFLAG = integer flag, returned as 0 if K and T were legal,
C -1 if K was illegal, and -2 if T was illegal.
C On an error return, a message is also written.
C-----------------------------------------------------------------------
C Part 3. Common Blocks.
C
C If DLSODI is to be used in an overlay situation, the user
C must declare, in the primary overlay, the variables in:
C (1) the call sequence to DLSODI, and
C (2) the internal Common block
C /DLS001/ of length 255 (218 double precision words
C followed by 37 integer words),
C
C If DLSODI is used on a system in which the contents of internal
C Common blocks are not preserved between calls, the user should
C declare the above Common block in the calling program to insure
C that their contents are preserved.
C
C If the solution of a given problem by DLSODI is to be interrupted
C and then later continued, such as when restarting an interrupted run
C or alternating between two or more problems, the user should save,
C following the return from the last DLSODI call prior to the
C interruption, the contents of the call sequence variables and the
C internal Common blocks, and later restore these values before the
C next DLSODI call for that problem. To save and restore the Common
C blocks, use Subroutine DSRCOM (see Part 2 above).
C
C-----------------------------------------------------------------------
C Part 4. Optionally Replaceable Solver Routines.
C
C Below are descriptions of two routines in the DLSODI package which
C relate to the measurement of errors. Either routine can be
C replaced by a user-supplied version, if desired. However, since such
C a replacement may have a major impact on performance, it should be
C done only when absolutely necessary, and only with great caution.
C (Note: The means by which the package version of a routine is
C superseded by the user's version may be system-dependent.)
C
C (a) DEWSET.
C The following subroutine is called just before each internal
C integration step, and sets the array of error weights, EWT, as
C described under ITOL/RTOL/ATOL above:
C SUBROUTINE DEWSET (NEQ, ITOL, RTOL, ATOL, YCUR, EWT)
C where NEQ, ITOL, RTOL, and ATOL are as in the DLSODI call sequence,
C YCUR contains the current dependent variable vector, and
C EWT is the array of weights set by DEWSET.
C
C If the user supplies this subroutine, it must return in EWT(i)
C (i = 1,...,NEQ) a positive quantity suitable for comparing errors
C in y(i) to. The EWT array returned by DEWSET is passed to the DVNORM
C routine (see below), and also used by DLSODI in the computation
C of the optional output IMXER, the diagonal Jacobian approximation,
C and the increments for difference quotient Jacobians.
C
C In the user-supplied version of DEWSET, it may be desirable to use
C the current values of derivatives of y. Derivatives up to order NQ
C are available from the history array YH, described above under
C optional outputs. In DEWSET, YH is identical to the YCUR array,
C extended to NQ + 1 columns with a column length of NYH and scale
C factors of H**j/factorial(j). On the first call for the problem,
C given by NST = 0, NQ is 1 and H is temporarily set to 1.0.
C NYH is the initial value of NEQ. The quantities NQ, H, and NST
C can be obtained by including in DEWSET the statements:
C DOUBLE PRECISION RLS
C COMMON /DLS001/ RLS(218),ILS(37)
C NQ = ILS(33)
C NST = ILS(34)
C H = RLS(212)
C 10 CONTINUE
C PA(1,1,N) = PA(1,1,N) + 1.0
C RETURN
C END
C
C SUBROUTINE JACBT (N, T, Y, S, MB, NB, PA, PB, PC)
C DOUBLE PRECISION T, Y, S, PA, PB, PC, ETA, DELX, EODSQ
C DIMENSION Y(N), S(N), PA(MB,MB,NB),PB(MB,MB,NB),PC(MB,MB,NB)
C DATA ETA/0.05/, DELX/0.05/
C EODSQ = ETA/DELX**2
C PA(1,1,1) = EODSQ
C PB(1,1,1) = -2.0*EODSQ
C PC(1,1,1) = EODSQ
C DO 10 K = 2,N
C PA(1,1,K) = -2.0*EODSQ
C PB(1,1,K) = -Y(K+1)*(0.5/DELX) + EODSQ
C PC(1,1,K) = Y(K-1)*(0.5/DELX) + EODSQ
C 10 CONTINUE
C PB(1,1,N) = EODSQ
C PC(1,1,N) = -2.0*EODSQ
C PA(1,1,N) = EODSQ
C RETURN
C END
C
C The output of this program (on a CDC-7600 in single precision)
C is as follows:
C
C At t = 0.10 No. steps = 35 No. r-s = 45 No. J-s = 9
C At t = 0.20 No. steps = 43 No. r-s = 54 No. J-s = 10
C At t = 0.30 No. steps = 48 No. r-s = 60 No. J-s = 11
C At t = 0.40 No. steps = 51 No. r-s = 64 No. J-s = 12
C
C Final solution values..
C 1.2747e-02 1.1997e-02 1.5560e-02 2.3767e-02 3.7224e-02
C 5.6646e-02 8.2645e-02 1.1557e-01 1.5541e-01 2.0177e-01
C 2.5397e-01 3.1104e-01 3.7189e-01 4.3530e-01 5.0000e-01
C 5.6472e-01 6.2816e-01 6.8903e-01 7.4612e-01 7.9829e-01
C 8.4460e-01 8.8438e-01 9.1727e-01 9.4330e-01 9.6281e-01
C 9.7632e-01 9.8426e-01 9.8648e-01 9.8162e-01 9.6617e-01
C 9.3374e-01 8.7535e-01 7.8236e-01 6.5321e-01 5.0003e-01
C 3.4709e-01 2.1876e-01 1.2771e-01 7.3671e-02 5.0642e-02
C 5.4496e-02
C
C-----------------------------------------------------------------------
C Full Description of User Interface to DLSOIBT.
C
C The user interface to DLSOIBT consists of the following parts.
C
C 1. The call sequence to Subroutine DLSOIBT, which is a driver
C routine for the solver. This includes descriptions of both
C the call sequence arguments and of user-supplied routines.
C Following these descriptions is a description of
C optional inputs available through the call sequence, and then
C a description of optional outputs (in the work arrays).
C
C 2. Descriptions of other routines in the DLSOIBT package that may be
C (optionally) called by the user. These provide the ability to
C alter error message handling, save and restore the internal
C Common, and obtain specified derivatives of the solution y(t).
C
C 3. Descriptions of Common blocks to be declared in overlay
C or similar environments, or to be saved when doing an interrupt
C of the problem and continued solution later.
C
C 4. Description of two routines in the DLSOIBT package, either of
C which the user may replace with his/her own version, if desired.
C These relate to the measurement of errors.
C
C-----------------------------------------------------------------------
C Part 1. Call Sequence.
C
C The call sequence parameters used for input only are
C RES, ADDA, JAC, NEQ, TOUT, ITOL, RTOL, ATOL, ITASK,
C IOPT, LRW, LIW, MF,
C and those used for both input and output are
C Y, T, ISTATE, YDOTI.
C The work arrays RWORK and IWORK are also used for additional and
C optional inputs and optional outputs. (The term output here refers
C to the return from Subroutine DLSOIBT to the user's calling program.)
C
C The legality of input parameters will be thoroughly checked on the
C initial call for the problem, but not checked thereafter unless a
C change in input parameters is flagged by ISTATE = 3 on input.
C
C The descriptions of the call arguments are as follows.
C
C RES = the name of the user-supplied subroutine which supplies
C the residual vector for the ODE system, defined by
C r = g(t,y) - A(t,y) * s
C as a function of the scalar t and the vectors
C s and y (s approximates dy/dt). This subroutine
C is to have the form
C SUBROUTINE RES (NEQ, T, Y, S, R, IRES)
C DOUBLE PRECISION T, Y(*), S(*), R(*)
C where NEQ, T, Y, S, and IRES are input, and R and
C IRES are output. Y, S, and R are arrays of length NEQ.
C On input, IRES indicates how DLSOIBT will use the
C returned array R, as follows:
C IRES = 1 means that DLSOIBT needs the full residual,
C r = g - A*s, exactly.
C IRES = -1 means that DLSOIBT is using R only to compute
C the Jacobian dr/dy by difference quotients.
C The RES routine can ignore IRES, or it can omit some terms
C if IRES = -1. If A does not depend on y, then RES can
C just return R = g when IRES = -1. If g - A*s contains other
C additive terms that are independent of y, these can also be
C dropped, if done consistently, when IRES = -1.
C The subroutine should set the flag IRES if it
C encounters a halt condition or illegal input.
C Otherwise, it should not reset IRES. On output,
C IRES = 1 or -1 represents a normal return, and
C DLSOIBT continues integrating the ODE. Leave IRES
C unchanged from its input value.
C IRES = 2 tells DLSOIBT to immediately return control
C to the calling program, with ISTATE = 3. This lets
C the calling program change parameters of the problem
C if necessary.
C IRES = 3 represents an error condition (for example, an
C illegal value of y). DLSOIBT tries to integrate the system
C without getting IRES = 3 from RES. If it cannot, DLSOIBT
C returns with ISTATE = -7 or -1.
C the default is not desired.
C The default value of LUN is 6.
C
C CALL XSETF(MFLAG) Set a flag to control the printing of
C messages by DLSOIBT.
C MFLAG = 0 means do not print. (Danger:
C This risks losing valuable information.)
C MFLAG = 1 means print (the default).
C
C Either of the above calls may be made at
C any time and will take effect immediately.
C
C CALL DSRCOM(RSAV,ISAV,JOB) saves and restores the contents of
C the internal Common blocks used by
C DLSOIBT (see Part 3 below).
C RSAV must be a real array of length 218
C or more, and ISAV must be an integer
C array of length 37 or more.
C JOB=1 means save Common into RSAV/ISAV.
C JOB=2 means restore Common from RSAV/ISAV.
C DSRCOM is useful if one is
C interrupting a run and restarting
C later, or alternating between two or
C more problems solved with DLSOIBT.
C
C CALL DINTDY(,,,,,) Provide derivatives of y, of various
C (see below) orders, at a specified point t, if
C desired. It may be called only after
C a successful return from DLSOIBT.
C
C The detailed instructions for using DINTDY are as follows.
C The form of the call is:
C
C CALL DINTDY (T, K, RWORK(21), NYH, DKY, IFLAG)
C
C The input parameters are:
C
C T = value of independent variable where answers are desired
C (normally the same as the t last returned by DLSOIBT).
C For valid results, T must lie between TCUR - HU and TCUR.
C (See optional outputs for TCUR and HU.)
C K = integer order of the derivative desired. K must satisfy
C 0 .le. K .le. NQCUR, where NQCUR is the current order
C (see optional outputs). The capability corresponding
C to K = 0, i.e. computing y(t), is already provided
C by DLSOIBT directly. Since NQCUR .ge. 1, the first
C derivative dy/dt is always available with DINTDY.
C RWORK(21) = the base address of the history array YH.
C NYH = column length of YH, equal to the initial value of NEQ.
C
C The output parameters are:
C
C DKY = a real array of length NEQ containing the computed value
C of the K-th derivative of y(t).
C IFLAG = integer flag, returned as 0 if K and T were legal,
C -1 if K was illegal, and -2 if T was illegal.
C On an error return, a message is also written.
C-----------------------------------------------------------------------
C Part 3. Common Blocks.
C
C If DLSOIBT is to be used in an overlay situation, the user
C must declare, in the primary overlay, the variables in:
C (1) the call sequence to DLSOIBT, and
C (2) the internal Common block
C /DLS001/ of length 255 (218 double precision words
C followed by 37 integer words),
C
C If DLSOIBT is used on a system in which the contents of internal
C Common blocks are not preserved between calls, the user should
C declare the above Common block in the calling program to insure
C that their contents are preserved.
C
C If the solution of a given problem by DLSOIBT is to be interrupted
C and then later continued, such as when restarting an interrupted run
C or alternating between two or more problems, the user should save,
C following the return from the last DLSOIBT call prior to the
C interruption, the contents of the call sequence variables and the
C internal Common blocks, and later restore these values before the
C next DLSOIBT call for that problem. To save and restore the Common
C blocks, use Subroutine DSRCOM (see Part 2 above).
C
C-----------------------------------------------------------------------
C Part 4. Optionally Replaceable Solver Routines.
C
C Below are descriptions of two routines in the DLSOIBT package which
C relate to the measurement of errors. Either routine can be
C replaced by a user-supplied version, if desired. However, since such
C a replacement may have a major impact on performance, it should be
C done only when absolutely necessary, and only with great caution.
C (Note: The means by which the package version of a routine is
C superseded by the user's version may be system-dependent.)
C
C (a) DEWSET.
C The following subroutine is called just before each internal
C integration step, and sets the array of error weights, EWT, as
C described under ITOL/RTOL/ATOL above:
C SUBROUTINE DEWSET (NEQ, ITOL, RTOL, ATOL, YCUR, EWT)
C where NEQ, ITOL, RTOL, and ATOL are as in the DLSOIBT call sequence,
C YCUR contains the current dependent variable vector, and
C EWT is the array of weights set by DEWSET.
C
C If the user supplies this subroutine, it must return in EWT(i)
C (i = 1,...,NEQ) a positive quantity suitable for comparing errors
C in y(i) to. The EWT array returned by DEWSET is passed to the DVNORM
C routine (see below), and also used by DLSOIBT in the computation
C of the optional output IMXER, the diagonal Jacobian approximation,
C and the increments for difference quotient Jacobians.
C
C In the user-supplied version of DEWSET, it may be desirable to use
C the current values of derivatives of y. Derivatives up to order NQ
C are available from the history array YH, described above under
C optional outputs. In DEWSET, YH is identical to the YCUR array,
C extended to NQ + 1 columns with a column length of NYH and scale
C factors of H**j/factorial(j). On the first call for the problem,
C given by NST = 0, NQ is 1 and H is temporarily set to 1.0.
C NYH is the initial value of NEQ. The quantities NQ, H, and NST
C can be obtained by including in DEWSET the statements:
C DOUBLE PRECISION RLS
C COMMON /DLS001/ RLS(218),ILS(37)
C NQ = ILS(33)
C NST = ILS(34)
C H = RLS(212)
C IF (J .EQ. 1) JM1 = N
C P(J) = P(J) + (2.0/3.0)
C P(JP1) = P(JP1) + (1.0/6.0)
C P(JM1) = P(JM1) + (1.0/6.0)
C RETURN
C END
C
C SUBROUTINE JACSP (N, T, Y, S, J, IP, JP, PDJ)
C DOUBLE PRECISION T, Y, S, PDJ, R4D, EODSQ
C DIMENSION Y(N), S(N), IP(*), JP(*), PDJ(N)
C COMMON /TEST1/ R4D, EODSQ, NM1
C JM1 = J - 1
C JP1 = J + 1
C IF (J .EQ. 1) JM1 = N
C IF (J .EQ. N) JP1 = 1
C PDJ(JM1) = -2.0*R4D*Y(J) + EODSQ
C PDJ(J) = -2.0*EODSQ
C PDJ(JP1) = 2.0*R4D*Y(J) + EODSQ
C RETURN
C END
C
C The output of this program (on a CDC-7600 in single precision)
C is as follows:
C
C At t = 0.10 No. steps = 15 Last step = 1.6863e-02
C At t = 0.20 No. steps = 19 Last step = 2.4101e-02
C At t = 0.30 No. steps = 22 Last step = 4.3143e-02
C At t = 0.40 No. steps = 24 Last step = 5.7819e-02
C
C Final solution values..
C 1.8371e-02 1.3578e-02 1.5864e-02 2.3805e-02 3.7245e-02
C 5.6630e-02 8.2538e-02 1.1538e-01 1.5522e-01 2.0172e-01
C 2.5414e-01 3.1150e-01 3.7259e-01 4.3608e-01 5.0060e-01
C 5.6482e-01 6.2751e-01 6.8758e-01 7.4415e-01 7.9646e-01
C 8.4363e-01 8.8462e-01 9.1853e-01 9.4500e-01 9.6433e-01
C 9.7730e-01 9.8464e-01 9.8645e-01 9.8138e-01 9.6584e-01
C 9.3336e-01 8.7497e-01 7.8213e-01 6.5315e-01 4.9997e-01
C 3.4672e-01 2.1758e-01 1.2461e-01 6.6208e-02 3.3784e-02
C
C Required RW size = 1409 IW size = 30
C No. steps = 24 No. r-s = 33 No. J-s = 8
C No. of nonzeros in P matrix = 120 No. of nonzeros in LU = 194
C
C-----------------------------------------------------------------------
C Full Description of User Interface to DLSODIS.
C
C The user interface to DLSODIS consists of the following parts.
C
C 1. The call sequence to Subroutine DLSODIS, which is a driver
C routine for the solver. This includes descriptions of both
C the call sequence arguments and of user-supplied routines.
C Following these descriptions is a description of
C optional inputs available through the call sequence, and then
C a description of optional outputs (in the work arrays).
C
C 2. Descriptions of other routines in the DLSODIS package that may be
C (optionally) called by the user. These provide the ability to
C alter error message handling, save and restore the internal
C Common, and obtain specified derivatives of the solution y(t).
C
C 3. Descriptions of Common blocks to be declared in overlay
C or similar environments, or to be saved when doing an interrupt
C of the problem and continued solution later.
C
C 4. Description of two routines in the DLSODIS package, either of
C which the user may replace with his/her own version, if desired.
C These relate to the measurement of errors.
C
C-----------------------------------------------------------------------
C Part 1. Call Sequence.
C
C The call sequence parameters used for input only are
C RES, ADDA, JAC, NEQ, TOUT, ITOL, RTOL, ATOL, ITASK,
C IOPT, LRW, LIW, MF,
C and those used for both input and output are
C Y, T, ISTATE, YDOTI.
C The work arrays RWORK and IWORK are also used for conditional and
C optional inputs and optional outputs. (The term output here refers
C to the return from Subroutine DLSODIS to the user's calling program.)
C
C The legality of input parameters will be thoroughly checked on the
C initial call for the problem, but not checked thereafter unless a
C change in input parameters is flagged by ISTATE = 3 on input.
C
C The descriptions of the call arguments are as follows.
C
C RES = the name of the user-supplied subroutine which supplies
C the residual vector for the ODE system, defined by
C r = g(t,y) - A(t,y) * s
C as a function of the scalar t and the vectors
C s and y (s approximates dy/dt). This subroutine
C is to have the form
C SUBROUTINE RES (NEQ, T, Y, S, R, IRES)
C DOUBLE PRECISION T, Y(*), S(*), R(*)
C where NEQ, T, Y, S, and IRES are input, and R and
C IRES are output. Y, S, and R are arrays of length NEQ.
C On input, IRES indicates how DLSODIS will use the
C returned array R, as follows:
C IRES = 1 means that DLSODIS needs the full residual,
C r = g - A*s, exactly.
C IRES = -1 means that DLSODIS is using R only to compute
C the Jacobian dr/dy by difference quotients.
C The RES routine can ignore IRES, or it can omit some terms
C if IRES = -1. If A does not depend on y, then RES can
C just return R = g when IRES = -1. If g - A*s contains other
C additive terms that are independent of y, these can also be
C dropped, if done consistently, when IRES = -1.
C The subroutine should set the flag IRES if it
C encounters a halt condition or illegal input.
C Otherwise, it should not reset IRES. On output,
C IRES = 1 or -1 represents a normal return, and
C DLSODIS continues integrating the ODE. Leave IRES
C unchanged from its input value.
C IRES = 2 tells DLSODIS to immediately return control
C to the calling program, with ISTATE = 3. This lets
C the calling program change parameters of the problem
C if necessary.
C IRES = 3 represents an error condition (for example, an
C illegal value of y). DLSODIS tries to integrate the system
C without getting IRES = 3 from RES. If it cannot, DLSODIS
C returns with ISTATE = -7 or -1.
C
C CALL XSETF(MFLAG) Set a flag to control the printing of
C messages by DLSODIS.
C MFLAG = 0 means do not print. (Danger:
C This risks losing valuable information.)
C MFLAG = 1 means print (the default).
C
C Either of the above calls may be made at
C any time and will take effect immediately.
C
C CALL DSRCMS(RSAV,ISAV,JOB) saves and restores the contents of
C the internal Common blocks used by
C DLSODIS (see Part 3 below).
C RSAV must be a real array of length 224
C or more, and ISAV must be an integer
C array of length 71 or more.
C JOB=1 means save Common into RSAV/ISAV.
C JOB=2 means restore Common from RSAV/ISAV.
C DSRCMS is useful if one is
C interrupting a run and restarting
C later, or alternating between two or
C more problems solved with DLSODIS.
C
C CALL DINTDY(,,,,,) Provide derivatives of y, of various
C (see below) orders, at a specified point t, if
C desired. It may be called only after
C a successful return from DLSODIS.
C
C The detailed instructions for using DINTDY are as follows.
C The form of the call is:
C
C LYH = IWORK(22)
C CALL DINTDY (T, K, RWORK(LYH), NYH, DKY, IFLAG)
C
C The input parameters are:
C
C T = value of independent variable where answers are desired
C (normally the same as the T last returned by DLSODIS).
C For valid results, T must lie between TCUR - HU and TCUR.
C (See optional outputs for TCUR and HU.)
C K = integer order of the derivative desired. K must satisfy
C 0 .le. K .le. NQCUR, where NQCUR is the current order
C (see optional outputs). The capability corresponding
C to K = 0, i.e. computing y(t), is already provided
C by DLSODIS directly. Since NQCUR .ge. 1, the first
C derivative dy/dt is always available with DINTDY.
C LYH = the base address of the history array YH, obtained
C as an optional output as shown above.
C NYH = column length of YH, equal to the initial value of NEQ.
C
C The output parameters are:
C
C DKY = a real array of length NEQ containing the computed value
C of the K-th derivative of y(t).
C IFLAG = integer flag, returned as 0 if K and T were legal,
C -1 if K was illegal, and -2 if T was illegal.
C On an error return, a message is also written.
C-----------------------------------------------------------------------
C Part 3. Common Blocks.
C
C If DLSODIS is to be used in an overlay situation, the user
C must declare, in the primary overlay, the variables in:
C (1) the call sequence to DLSODIS, and
C (2) the two internal Common blocks
C /DLS001/ of length 255 (218 double precision words
C followed by 37 integer words),
C /DLSS01/ of length 40 (6 double precision words
C followed by 34 integer words).
C
C If DLSODIS is used on a system in which the contents of internal
C Common blocks are not preserved between calls, the user should
C declare the above Common blocks in the calling program to insure
C that their contents are preserved.
C
C If the solution of a given problem by DLSODIS is to be interrupted
C and then later continued, such as when restarting an interrupted run
C or alternating between two or more problems, the user should save,
C following the return from the last DLSODIS call prior to the
C interruption, the contents of the call sequence variables and the
C internal Common blocks, and later restore these values before the
C next DLSODIS call for that problem. To save and restore the Common
C blocks, use Subroutines DSRCMS (see Part 2 above).
C
C-----------------------------------------------------------------------
C Part 4. Optionally Replaceable Solver Routines.
C
C Below are descriptions of two routines in the DLSODIS package which
C relate to the measurement of errors. Either routine can be
C replaced by a user-supplied version, if desired. However, since such
C a replacement may have a major impact on performance, it should be
C done only when absolutely necessary, and only with great caution.
C (Note: The means by which the package version of a routine is
C superseded by the user's version may be system-dependent.)
C
C (a) DEWSET.
C The following subroutine is called just before each internal
C integration step, and sets the array of error weights, EWT, as
C described under ITOL/RTOL/ATOL above:
C SUBROUTINE DEWSET (NEQ, ITOL, RTOL, ATOL, YCUR, EWT)
C where NEQ, ITOL, RTOL, and ATOL are as in the DLSODIS call sequence,
C YCUR contains the current dependent variable vector, and
C EWT is the array of weights set by DEWSET.
C
C If the user supplies this subroutine, it must return in EWT(i)
C (i = 1,...,NEQ) a positive quantity suitable for comparing errors
C in y(i) to. The EWT array returned by DEWSET is passed to the DVNORM
C routine (see below), and also used by DLSODIS in the computation
C of the optional output IMXER, and the increments for difference
C quotient Jacobians.
C
C In the user-supplied version of DEWSET, it may be desirable to use
C the current values of derivatives of y. Derivatives up to order NQ
C are available from the history array YH, described above under
C optional outputs. In DEWSET, YH is identical to the YCUR array,
C extended to NQ + 1 columns with a column length of NYH and scale
C factors of H**j/factorial(j). On the first call for the problem,
C given by NST = 0, NQ is 1 and H is temporarily set to 1.0.
C NYH is the initial value of NEQ. The quantities NQ, H, and NST
C can be obtained by including in DEWSET the statements:
C DOUBLE PRECISION RLS
C COMMON /DLS001/ RLS(218),ILS(37)
C NQ = ILS(33)
( run in 2.990 seconds using v1.01-cache-2.11-cpan-39bf76dae61 )