SLATEC Routines --- DDASSL ---
*DECK DDASSL
SUBROUTINE DDASSL (RES, NEQ, T, Y, YPRIME, TOUT, INFO, RTOL, ATOL,
* IDID, RWORK, LRW, IWORK, LIW, RPAR, IPAR, JAC)
C***BEGIN PROLOGUE DDASSL
C***PURPOSE This code solves a system of differential/algebraic
C equations of the form G(T,Y,YPRIME) = 0.
C***LIBRARY SLATEC (DASSL)
C***CATEGORY I1A2
C***TYPE DOUBLE PRECISION (SDASSL-S, DDASSL-D)
C***KEYWORDS BACKWARD DIFFERENTIATION FORMULAS, DASSL,
C DIFFERENTIAL/ALGEBRAIC, IMPLICIT DIFFERENTIAL SYSTEMS
C***AUTHOR Petzold, Linda R., (LLNL)
C Computing and Mathematics Research Division
C Lawrence Livermore National Laboratory
C L - 316, P.O. Box 808,
C Livermore, CA. 94550
C***DESCRIPTION
C
C *Usage:
C
C EXTERNAL RES, JAC
C INTEGER NEQ, INFO(N), IDID, LRW, LIW, IWORK(LIW), IPAR
C DOUBLE PRECISION T, Y(NEQ), YPRIME(NEQ), TOUT, RTOL, ATOL,
C * RWORK(LRW), RPAR
C
C CALL DDASSL (RES, NEQ, T, Y, YPRIME, TOUT, INFO, RTOL, ATOL,
C * IDID, RWORK, LRW, IWORK, LIW, RPAR, IPAR, JAC)
C
C
C *Arguments:
C (In the following, all real arrays should be type DOUBLE PRECISION.)
C
C RES:EXT This is a subroutine which you provide to define the
C differential/algebraic system.
C
C NEQ:IN This is the number of equations to be solved.
C
C T:INOUT This is the current value of the independent variable.
C
C Y(*):INOUT This array contains the solution components at T.
C
C YPRIME(*):INOUT This array contains the derivatives of the solution
C components at T.
C
C TOUT:IN This is a point at which a solution is desired.
C
C INFO(N):IN The basic task of the code is to solve the system from T
C to TOUT and return an answer at TOUT. INFO is an integer
C array which is used to communicate exactly how you want
C this task to be carried out. (See below for details.)
C N must be greater than or equal to 15.
C
C RTOL,ATOL:INOUT These quantities represent relative and absolute
C error tolerances which you provide to indicate how
C accurately you wish the solution to be computed. You
C may choose them to be both scalars or else both vectors.
C Caution: In Fortran 77, a scalar is not the same as an
C array of length 1. Some compilers may object
C to using scalars for RTOL,ATOL.
C
C IDID:OUT This scalar quantity is an indicator reporting what the
C code did. You must monitor this integer variable to
C decide what action to take next.
C
C RWORK:WORK A real work array of length LRW which provides the
C code with needed storage space.
C
C LRW:IN The length of RWORK. (See below for required length.)
C
C IWORK:WORK An integer work array of length LIW which provides the
C code with needed storage space.
C
C LIW:IN The length of IWORK. (See below for required length.)
C
C RPAR,IPAR:IN These are real and integer parameter arrays which
C you can use for communication between your calling
C program and the RES subroutine (and the JAC subroutine)
C
C JAC:EXT This is the name of a subroutine which you may choose
C to provide for defining a matrix of partial derivatives
C described below.
C
C Quantities which may be altered by DDASSL are:
C T, Y(*), YPRIME(*), INFO(1), RTOL, ATOL,
C IDID, RWORK(*) AND IWORK(*)
C
C *Description
C
C Subroutine DDASSL uses the backward differentiation formulas of
C orders one through five to solve a system of the above form for Y and
C YPRIME. Values for Y and YPRIME at the initial time must be given as
C input. These values must be consistent, (that is, if T,Y,YPRIME are
C the given initial values, they must satisfy G(T,Y,YPRIME) = 0.). The
C subroutine solves the system from T to TOUT. It is easy to continue
C the solution to get results at additional TOUT. This is the interval
C mode of operation. Intermediate results can also be obtained easily
C by using the intermediate-output capability.
C
C The following detailed description is divided into subsections:
C 1. Input required for the first call to DDASSL.
C 2. Output after any return from DDASSL.
C 3. What to do to continue the integration.
C 4. Error messages.
C
C
C -------- INPUT -- WHAT TO DO ON THE FIRST CALL TO DDASSL ------------
C
C The first call of the code is defined to be the start of each new
C problem. Read through the descriptions of all the following items,
C provide sufficient storage space for designated arrays, set
C appropriate variables for the initialization of the problem, and
C give information about how you want the problem to be solved.
C
C
C RES -- Provide a subroutine of the form
C SUBROUTINE RES(T,Y,YPRIME,DELTA,IRES,RPAR,IPAR)
C to define the system of differential/algebraic
C equations which is to be solved. For the given values
C of T,Y and YPRIME, the subroutine should
C return the residual of the differential/algebraic
C system
C DELTA = G(T,Y,YPRIME)
C (DELTA(*) is a vector of length NEQ which is
C output for RES.)
C
C Subroutine RES must not alter T,Y or YPRIME.
C You must declare the name RES in an external
C statement in your program that calls DDASSL.
C You must dimension Y,YPRIME and DELTA in RES.
C
C IRES is an integer flag which is always equal to
C zero on input. Subroutine RES should alter IRES
C only if it encounters an illegal value of Y or
C a stop condition. Set IRES = -1 if an input value
C is illegal, and DDASSL will try to solve the problem
C without getting IRES = -1. If IRES = -2, DDASSL
C will return control to the calling program
C with IDID = -11.
C
C RPAR and IPAR are real and integer parameter arrays which
C you can use for communication between your calling program
C and subroutine RES. They are not altered by DDASSL. If you
C do not need RPAR or IPAR, ignore these parameters by treat-
C ing them as dummy arguments. If you do choose to use them,
C dimension them in your calling program and in RES as arrays
C of appropriate length.
C
C NEQ -- Set it to the number of differential equations.
C (NEQ .GE. 1)
C
C T -- Set it to the initial point of the integration.
C T must be defined as a variable.
C
C Y(*) -- Set this vector to the initial values of the NEQ solution
C components at the initial point. You must dimension Y of
C length at least NEQ in your calling program.
C
C YPRIME(*) -- Set this vector to the initial values of the NEQ
C first derivatives of the solution components at the initial
C point. You must dimension YPRIME at least NEQ in your
C calling program. If you do not know initial values of some
C of the solution components, see the explanation of INFO(11).
C
C TOUT -- Set it to the first point at which a solution
C is desired. You can not take TOUT = T.
C integration either forward in T (TOUT .GT. T) or
C backward in T (TOUT .LT. T) is permitted.
C
C The code advances the solution from T to TOUT using
C step sizes which are automatically selected so as to
C achieve the desired accuracy. If you wish, the code will
C return with the solution and its derivative at
C intermediate steps (intermediate-output mode) so that
C you can monitor them, but you still must provide TOUT in
C accord with the basic aim of the code.
C
C The first step taken by the code is a critical one
C because it must reflect how fast the solution changes near
C the initial point. The code automatically selects an
C initial step size which is practically always suitable for
C the problem. By using the fact that the code will not step
C past TOUT in the first step, you could, if necessary,
C restrict the length of the initial step size.
C
C For some problems it may not be permissible to integrate
C past a point TSTOP because a discontinuity occurs there
C or the solution or its derivative is not defined beyond
C TSTOP. When you have declared a TSTOP point (SEE INFO(4)
C and RWORK(1)), you have told the code not to integrate
C past TSTOP. In this case any TOUT beyond TSTOP is invalid
C input.
C
C INFO(*) -- Use the INFO array to give the code more details about
C how you want your problem solved. This array should be
C dimensioned of length 15, though DDASSL uses only the first
C eleven entries. You must respond to all of the following
C items, which are arranged as questions. The simplest use
C of the code corresponds to answering all questions as yes,
C i.e. setting all entries of INFO to 0.
C
C INFO(1) - This parameter enables the code to initialize
C itself. You must set it to indicate the start of every
C new problem.
C
C **** Is this the first call for this problem ...
C Yes - Set INFO(1) = 0
C No - Not applicable here.
C See below for continuation calls. ****
C
C INFO(2) - How much accuracy you want of your solution
C is specified by the error tolerances RTOL and ATOL.
C The simplest use is to take them both to be scalars.
C To obtain more flexibility, they can both be vectors.
C The code must be told your choice.
C
C **** Are both error tolerances RTOL, ATOL scalars ...
C Yes - Set INFO(2) = 0
C and input scalars for both RTOL and ATOL
C No - Set INFO(2) = 1
C and input arrays for both RTOL and ATOL ****
C
C INFO(3) - The code integrates from T in the direction
C of TOUT by steps. If you wish, it will return the
C computed solution and derivative at the next
C intermediate step (the intermediate-output mode) or
C TOUT, whichever comes first. This is a good way to
C proceed if you want to see the behavior of the solution.
C If you must have solutions at a great many specific
C TOUT points, this code will compute them efficiently.
C
C **** Do you want the solution only at
C TOUT (and not at the next intermediate step) ...
C Yes - Set INFO(3) = 0
C No - Set INFO(3) = 1 ****
C
C INFO(4) - To handle solutions at a great many specific
C values TOUT efficiently, this code may integrate past
C TOUT and interpolate to obtain the result at TOUT.
C Sometimes it is not possible to integrate beyond some
C point TSTOP because the equation changes there or it is
C not defined past TSTOP. Then you must tell the code
C not to go past.
C
C **** Can the integration be carried out without any
C restrictions on the independent variable T ...
C Yes - Set INFO(4)=0
C No - Set INFO(4)=1
C and define the stopping point TSTOP by
C setting RWORK(1)=TSTOP ****
C
C INFO(5) - To solve differential/algebraic problems it is
C necessary to use a matrix of partial derivatives of the
C system of differential equations. If you do not
C provide a subroutine to evaluate it analytically (see
C description of the item JAC in the call list), it will
C be approximated by numerical differencing in this code.
C although it is less trouble for you to have the code
C compute partial derivatives by numerical differencing,
C the solution will be more reliable if you provide the
C derivatives via JAC. Sometimes numerical differencing
C is cheaper than evaluating derivatives in JAC and
C sometimes it is not - this depends on your problem.
C
C **** Do you want the code to evaluate the partial
C derivatives automatically by numerical differences ...
C Yes - Set INFO(5)=0
C No - Set INFO(5)=1
C and provide subroutine JAC for evaluating the
C matrix of partial derivatives ****
C
C INFO(6) - DDASSL will perform much better if the matrix of
C partial derivatives, DG/DY + CJ*DG/DYPRIME,
C (here CJ is a scalar determined by DDASSL)
C is banded and the code is told this. In this
C case, the storage needed will be greatly reduced,
C numerical differencing will be performed much cheaper,
C and a number of important algorithms will execute much
C faster. The differential equation is said to have
C half-bandwidths ML (lower) and MU (upper) if equation i
C involves only unknowns Y(J) with
C I-ML .LE. J .LE. I+MU
C for all I=1,2,...,NEQ. Thus, ML and MU are the widths
C of the lower and upper parts of the band, respectively,
C with the main diagonal being excluded. If you do not
C indicate that the equation has a banded matrix of partial
C derivatives, the code works with a full matrix of NEQ**2
C elements (stored in the conventional way). Computations
C with banded matrices cost less time and storage than with
C full matrices if 2*ML+MU .LT. NEQ. If you tell the
C code that the matrix of partial derivatives has a banded
C structure and you want to provide subroutine JAC to
C compute the partial derivatives, then you must be careful
C to store the elements of the matrix in the special form
C indicated in the description of JAC.
C
C **** Do you want to solve the problem using a full
C (dense) matrix (and not a special banded
C structure) ...
C Yes - Set INFO(6)=0
C No - Set INFO(6)=1
C and provide the lower (ML) and upper (MU)
C bandwidths by setting
C IWORK(1)=ML
C IWORK(2)=MU ****
C
C
C INFO(7) -- You can specify a maximum (absolute value of)
C stepsize, so that the code
C will avoid passing over very
C large regions.
C
C **** Do you want the code to decide
C on its own maximum stepsize?
C Yes - Set INFO(7)=0
C No - Set INFO(7)=1
C and define HMAX by setting
C RWORK(2)=HMAX ****
C
C INFO(8) -- Differential/algebraic problems
C may occasionally suffer from
C severe scaling difficulties on the
C first step. If you know a great deal
C about the scaling of your problem, you can
C help to alleviate this problem by
C specifying an initial stepsize HO.
C
C **** Do you want the code to define
C its own initial stepsize?
C Yes - Set INFO(8)=0
C No - Set INFO(8)=1
C and define HO by setting
C RWORK(3)=HO ****
C
C INFO(9) -- If storage is a severe problem,
C you can save some locations by
C restricting the maximum order MAXORD.
C the default value is 5. for each
C order decrease below 5, the code
C requires NEQ fewer locations, however
C it is likely to be slower. In any
C case, you must have 1 .LE. MAXORD .LE. 5
C **** Do you want the maximum order to
C default to 5?
C Yes - Set INFO(9)=0
C No - Set INFO(9)=1
C and define MAXORD by setting
C IWORK(3)=MAXORD ****
C
C INFO(10) --If you know that the solutions to your equations
C will always be nonnegative, it may help to set this
C parameter. However, it is probably best to
C try the code without using this option first,
C and only to use this option if that doesn't
C work very well.
C **** Do you want the code to solve the problem without
C invoking any special nonnegativity constraints?
C Yes - Set INFO(10)=0
C No - Set INFO(10)=1
C
C INFO(11) --DDASSL normally requires the initial T,
C Y, and YPRIME to be consistent. That is,
C you must have G(T,Y,YPRIME) = 0 at the initial
C time. If you do not know the initial
C derivative precisely, you can let DDASSL try
C to compute it.
C **** Are the initial T, Y, YPRIME consistent?
C Yes - Set INFO(11) = 0
C No - Set INFO(11) = 1,
C and set YPRIME to an initial approximation
C to YPRIME. (If you have no idea what
C YPRIME should be, set it to zero. Note
C that the initial Y should be such
C that there must exist a YPRIME so that
C G(T,Y,YPRIME) = 0.)
C
C RTOL, ATOL -- You must assign relative (RTOL) and absolute (ATOL
C error tolerances to tell the code how accurately you
C want the solution to be computed. They must be defined
C as variables because the code may change them. You
C have two choices --
C Both RTOL and ATOL are scalars. (INFO(2)=0)
C Both RTOL and ATOL are vectors. (INFO(2)=1)
C in either case all components must be non-negative.
C
C The tolerances are used by the code in a local error
C test at each step which requires roughly that
C ABS(LOCAL ERROR) .LE. RTOL*ABS(Y)+ATOL
C for each vector component.
C (More specifically, a root-mean-square norm is used to
C measure the size of vectors, and the error test uses the
C magnitude of the solution at the beginning of the step.)
C
C The true (global) error is the difference between the
C true solution of the initial value problem and the
C computed approximation. Practically all present day
C codes, including this one, control the local error at
C each step and do not even attempt to control the global
C error directly.
C Usually, but not always, the true accuracy of the
C computed Y is comparable to the error tolerances. This
C code will usually, but not always, deliver a more
C accurate solution if you reduce the tolerances and
C integrate again. By comparing two such solutions you
C can get a fairly reliable idea of the true error in the
C solution at the bigger tolerances.
C
C Setting ATOL=0. results in a pure relative error test on
C that component. Setting RTOL=0. results in a pure
C absolute error test on that component. A mixed test
C with non-zero RTOL and ATOL corresponds roughly to a
C relative error test when the solution component is much
C bigger than ATOL and to an absolute error test when the
C solution component is smaller than the threshhold ATOL.
C
C The code will not attempt to compute a solution at an
C accuracy unreasonable for the machine being used. It will
C advise you if you ask for too much accuracy and inform
C you as to the maximum accuracy it believes possible.
C
C RWORK(*) -- Dimension this real work array of length LRW in your
C calling program.
C
C LRW -- Set it to the declared length of the RWORK array.
C You must have
C LRW .GE. 40+(MAXORD+4)*NEQ+NEQ**2
C for the full (dense) JACOBIAN case (when INFO(6)=0), or
C LRW .GE. 40+(MAXORD+4)*NEQ+(2*ML+MU+1)*NEQ
C for the banded user-defined JACOBIAN case
C (when INFO(5)=1 and INFO(6)=1), or
C LRW .GE. 40+(MAXORD+4)*NEQ+(2*ML+MU+1)*NEQ
C +2*(NEQ/(ML+MU+1)+1)
C for the banded finite-difference-generated JACOBIAN case
C (when INFO(5)=0 and INFO(6)=1)
C
C IWORK(*) -- Dimension this integer work array of length LIW in
C your calling program.
C
C LIW -- Set it to the declared length of the IWORK array.
C You must have LIW .GE. 20+NEQ
C
C RPAR, IPAR -- These are parameter arrays, of real and integer
C type, respectively. You can use them for communication
C between your program that calls DDASSL and the
C RES subroutine (and the JAC subroutine). They are not
C altered by DDASSL. If you do not need RPAR or IPAR,
C ignore these parameters by treating them as dummy
C arguments. If you do choose to use them, dimension
C them in your calling program and in RES (and in JAC)
C as arrays of appropriate length.
C
C JAC -- If you have set INFO(5)=0, you can ignore this parameter
C by treating it as a dummy argument. Otherwise, you must
C provide a subroutine of the form
C SUBROUTINE JAC(T,Y,YPRIME,PD,CJ,RPAR,IPAR)
C to define the matrix of partial derivatives
C PD=DG/DY+CJ*DG/DYPRIME
C CJ is a scalar which is input to JAC.
C For the given values of T,Y,YPRIME, the
C subroutine must evaluate the non-zero partial
C derivatives for each equation and each solution
C component, and store these values in the
C matrix PD. The elements of PD are set to zero
C before each call to JAC so only non-zero elements
C need to be defined.
C
C Subroutine JAC must not alter T,Y,(*),YPRIME(*), or CJ.
C You must declare the name JAC in an EXTERNAL statement in
C your program that calls DDASSL. You must dimension Y,
C YPRIME and PD in JAC.
C
C The way you must store the elements into the PD matrix
C depends on the structure of the matrix which you
C indicated by INFO(6).
C *** INFO(6)=0 -- Full (dense) matrix ***
C Give PD a first dimension of NEQ.
C When you evaluate the (non-zero) partial derivative
C of equation I with respect to variable J, you must
C store it in PD according to
C PD(I,J) = "DG(I)/DY(J)+CJ*DG(I)/DYPRIME(J)"
C *** INFO(6)=1 -- Banded JACOBIAN with ML lower and MU
C upper diagonal bands (refer to INFO(6) description
C of ML and MU) ***
C Give PD a first dimension of 2*ML+MU+1.
C when you evaluate the (non-zero) partial derivative
C of equation I with respect to variable J, you must
C store it in PD according to
C IROW = I - J + ML + MU + 1
C PD(IROW,J) = "DG(I)/DY(J)+CJ*DG(I)/DYPRIME(J)"
C
C RPAR and IPAR are real and integer parameter arrays
C which you can use for communication between your calling
C program and your JACOBIAN subroutine JAC. They are not
C altered by DDASSL. If you do not need RPAR or IPAR,
C ignore these parameters by treating them as dummy
C arguments. If you do choose to use them, dimension
C them in your calling program and in JAC as arrays of
C appropriate length.
C
C
C OPTIONALLY REPLACEABLE NORM ROUTINE:
C
C DDASSL uses a weighted norm DDANRM to measure the size
C of vectors such as the estimated error in each step.
C A FUNCTION subprogram
C DOUBLE PRECISION FUNCTION DDANRM(NEQ,V,WT,RPAR,IPAR)
C DIMENSION V(NEQ),WT(NEQ)
C is used to define this norm. Here, V is the vector
C whose norm is to be computed, and WT is a vector of
C weights. A DDANRM routine has been included with DDASSL
C which computes the weighted root-mean-square norm
C given by
C DDANRM=SQRT((1/NEQ)*SUM(V(I)/WT(I))**2)
C this norm is suitable for most problems. In some
C special cases, it may be more convenient and/or
C efficient to define your own norm by writing a function
C subprogram to be called instead of DDANRM. This should,
C however, be attempted only after careful thought and
C consideration.
C
C
C -------- OUTPUT -- AFTER ANY RETURN FROM DDASSL ---------------------
C
C The principal aim of the code is to return a computed solution at
C TOUT, although it is also possible to obtain intermediate results
C along the way. To find out whether the code achieved its goal
C or if the integration process was interrupted before the task was
C completed, you must check the IDID parameter.
C
C
C T -- The solution was successfully advanced to the
C output value of T.
C
C Y(*) -- Contains the computed solution approximation at T.
C
C YPRIME(*) -- Contains the computed derivative
C approximation at T.
C
C IDID -- Reports what the code did.
C
C *** Task completed ***
C Reported by positive values of IDID
C
C IDID = 1 -- A step was successfully taken in the
C intermediate-output mode. The code has not
C yet reached TOUT.
C
C IDID = 2 -- The integration to TSTOP was successfully
C completed (T=TSTOP) by stepping exactly to TSTOP.
C
C IDID = 3 -- The integration to TOUT was successfully
C completed (T=TOUT) by stepping past TOUT.
C Y(*) is obtained by interpolation.
C YPRIME(*) is obtained by interpolation.
C
C *** Task interrupted ***
C Reported by negative values of IDID
C
C IDID = -1 -- A large amount of work has been expended.
C (About 500 steps)
C
C IDID = -2 -- The error tolerances are too stringent.
C
C IDID = -3 -- The local error test cannot be satisfied
C because you specified a zero component in ATOL
C and the corresponding computed solution
C component is zero. Thus, a pure relative error
C test is impossible for this component.
C
C IDID = -6 -- DDASSL had repeated error test
C failures on the last attempted step.
C
C IDID = -7 -- The corrector could not converge.
C
C IDID = -8 -- The matrix of partial derivatives
C is singular.
C
C IDID = -9 -- The corrector could not converge.
C there were repeated error test failures
C in this step.
C
C IDID =-10 -- The corrector could not converge
C because IRES was equal to minus one.
C
C IDID =-11 -- IRES equal to -2 was encountered
C and control is being returned to the
C calling program.
C
C IDID =-12 -- DDASSL failed to compute the initial
C YPRIME.
C
C
C
C IDID = -13,..,-32 -- Not applicable for this code
C
C *** Task terminated ***
C Reported by the value of IDID=-33
C
C IDID = -33 -- The code has encountered trouble from which
C it cannot recover. A message is printed
C explaining the trouble and control is returned
C to the calling program. For example, this occurs
C when invalid input is detected.
C
C RTOL, ATOL -- These quantities remain unchanged except when
C IDID = -2. In this case, the error tolerances have been
C increased by the code to values which are estimated to
C be appropriate for continuing the integration. However,
C the reported solution at T was obtained using the input
C values of RTOL and ATOL.
C
C RWORK, IWORK -- Contain information which is usually of no
C interest to the user but necessary for subsequent calls.
C However, you may find use for
C
C RWORK(3)--Which contains the step size H to be
C attempted on the next step.
C
C RWORK(4)--Which contains the current value of the
C independent variable, i.e., the farthest point
C integration has reached. This will be different
C from T only when interpolation has been
C performed (IDID=3).
C
C RWORK(7)--Which contains the stepsize used
C on the last successful step.
C
C IWORK(7)--Which contains the order of the method to
C be attempted on the next step.
C
C IWORK(8)--Which contains the order of the method used
C on the last step.
C
C IWORK(11)--Which contains the number of steps taken so
C far.
C
C IWORK(12)--Which contains the number of calls to RES
C so far.
C
C IWORK(13)--Which contains the number of evaluations of
C the matrix of partial derivatives needed so
C far.
C
C IWORK(14)--Which contains the total number
C of error test failures so far.
C
C IWORK(15)--Which contains the total number
C of convergence test failures so far.
C (includes singular iteration matrix
C failures.)
C
C
C -------- INPUT -- WHAT TO DO TO CONTINUE THE INTEGRATION ------------
C (CALLS AFTER THE FIRST)
C
C This code is organized so that subsequent calls to continue the
C integration involve little (if any) additional effort on your
C part. You must monitor the IDID parameter in order to determine
C what to do next.
C
C Recalling that the principal task of the code is to integrate
C from T to TOUT (the interval mode), usually all you will need
C to do is specify a new TOUT upon reaching the current TOUT.
C
C Do not alter any quantity not specifically permitted below,
C in particular do not alter NEQ,T,Y(*),YPRIME(*),RWORK(*),IWORK(*)
C or the differential equation in subroutine RES. Any such
C alteration constitutes a new problem and must be treated as such,
C i.e., you must start afresh.
C
C You cannot change from vector to scalar error control or vice
C versa (INFO(2)), but you can change the size of the entries of
C RTOL, ATOL. Increasing a tolerance makes the equation easier
C to integrate. Decreasing a tolerance will make the equation
C harder to integrate and should generally be avoided.
C
C You can switch from the intermediate-output mode to the
C interval mode (INFO(3)) or vice versa at any time.
C
C If it has been necessary to prevent the integration from going
C past a point TSTOP (INFO(4), RWORK(1)), keep in mind that the
C code will not integrate to any TOUT beyond the currently
C specified TSTOP. Once TSTOP has been reached you must change
C the value of TSTOP or set INFO(4)=0. You may change INFO(4)
C or TSTOP at any time but you must supply the value of TSTOP in
C RWORK(1) whenever you set INFO(4)=1.
C
C Do not change INFO(5), INFO(6), IWORK(1), or IWORK(2)
C unless you are going to restart the code.
C
C *** Following a completed task ***
C If
C IDID = 1, call the code again to continue the integration
C another step in the direction of TOUT.
C
C IDID = 2 or 3, define a new TOUT and call the code again.
C TOUT must be different from T. You cannot change
C the direction of integration without restarting.
C
C *** Following an interrupted task ***
C To show the code that you realize the task was
C interrupted and that you want to continue, you
C must take appropriate action and set INFO(1) = 1
C If
C IDID = -1, The code has taken about 500 steps.
C If you want to continue, set INFO(1) = 1 and
C call the code again. An additional 500 steps
C will be allowed.
C
C IDID = -2, The error tolerances RTOL, ATOL have been
C increased to values the code estimates appropriate
C for continuing. You may want to change them
C yourself. If you are sure you want to continue
C with relaxed error tolerances, set INFO(1)=1 and
C call the code again.
C
C IDID = -3, A solution component is zero and you set the
C corresponding component of ATOL to zero. If you
C are sure you want to continue, you must first
C alter the error criterion to use positive values
C for those components of ATOL corresponding to zero
C solution components, then set INFO(1)=1 and call
C the code again.
C
C IDID = -4,-5 --- Cannot occur with this code.
C
C IDID = -6, Repeated error test failures occurred on the
C last attempted step in DDASSL. A singularity in the
C solution may be present. If you are absolutely
C certain you want to continue, you should restart
C the integration. (Provide initial values of Y and
C YPRIME which are consistent)
C
C IDID = -7, Repeated convergence test failures occurred
C on the last attempted step in DDASSL. An inaccurate
C or ill-conditioned JACOBIAN may be the problem. If
C you are absolutely certain you want to continue, you
C should restart the integration.
C
C IDID = -8, The matrix of partial derivatives is singular.
C Some of your equations may be redundant.
C DDASSL cannot solve the problem as stated.
C It is possible that the redundant equations
C could be removed, and then DDASSL could
C solve the problem. It is also possible
C that a solution to your problem either
C does not exist or is not unique.
C
C IDID = -9, DDASSL had multiple convergence test
C failures, preceded by multiple error
C test failures, on the last attempted step.
C It is possible that your problem
C is ill-posed, and cannot be solved
C using this code. Or, there may be a
C discontinuity or a singularity in the
C solution. If you are absolutely certain
C you want to continue, you should restart
C the integration.
C
C IDID =-10, DDASSL had multiple convergence test failures
C because IRES was equal to minus one.
C If you are absolutely certain you want
C to continue, you should restart the
C integration.
C
C IDID =-11, IRES=-2 was encountered, and control is being
C returned to the calling program.
C
C IDID =-12, DDASSL failed to compute the initial YPRIME.
C This could happen because the initial
C approximation to YPRIME was not very good, or
C if a YPRIME consistent with the initial Y
C does not exist. The problem could also be caused
C by an inaccurate or singular iteration matrix.
C
C IDID = -13,..,-32 --- Cannot occur with this code.
C
C
C *** Following a terminated task ***
C
C If IDID= -33, you cannot continue the solution of this problem.
C An attempt to do so will result in your
C run being terminated.
C
C
C -------- ERROR MESSAGES ---------------------------------------------
C
C The SLATEC error print routine XERMSG is called in the event of
C unsuccessful completion of a task. Most of these are treated as
C "recoverable errors", which means that (unless the user has directed
C otherwise) control will be returned to the calling program for
C possible action after the message has been printed.
C
C In the event of a negative value of IDID other than -33, an appro-
C priate message is printed and the "error number" printed by XERMSG
C is the value of IDID. There are quite a number of illegal input
C errors that can lead to a returned value IDID=-33. The conditions
C and their printed "error numbers" are as follows:
C
C Error number Condition
C
C 1 Some element of INFO vector is not zero or one.
C 2 NEQ .le. 0
C 3 MAXORD not in range.
C 4 LRW is less than the required length for RWORK.
C 5 LIW is less than the required length for IWORK.
C 6 Some element of RTOL is .lt. 0
C 7 Some element of ATOL is .lt. 0
C 8 All elements of RTOL and ATOL are zero.
C 9 INFO(4)=1 and TSTOP is behind TOUT.
C 10 HMAX .lt. 0.0
C 11 TOUT is behind T.
C 12 INFO(8)=1 and H0=0.0
C 13 Some element of WT is .le. 0.0
C 14 TOUT is too close to T to start integration.
C 15 INFO(4)=1 and TSTOP is behind T.
C 16 --( Not used in this version )--
C 17 ML illegal. Either .lt. 0 or .gt. NEQ
C 18 MU illegal. Either .lt. 0 or .gt. NEQ
C 19 TOUT = T.
C
C If DDASSL is called again without any action taken to remove the
C cause of an unsuccessful return, XERMSG will be called with a fatal
C error flag, which will cause unconditional termination of the
C program. There are two such fatal errors:
C
C Error number -998: The last step was terminated with a negative
C value of IDID other than -33, and no appropriate action was
C taken.
C
C Error number -999: The previous call was terminated because of
C illegal input (IDID=-33) and there is illegal input in the
C present call, as well. (Suspect infinite loop.)
C
C ---------------------------------------------------------------------
C
C***REFERENCES A DESCRIPTION OF DASSL: A DIFFERENTIAL/ALGEBRAIC
C SYSTEM SOLVER, L. R. PETZOLD, SAND82-8637,
C SANDIA NATIONAL LABORATORIES, SEPTEMBER 1982.
C***ROUTINES CALLED D1MACH, DDAINI, DDANRM, DDASTP, DDATRP, DDAWTS,
C XERMSG
C***REVISION HISTORY (YYMMDD)
C 830315 DATE WRITTEN
C 880387 Code changes made. All common statements have been
C replaced by a DATA statement, which defines pointers into
C RWORK, and PARAMETER statements which define pointers
C into IWORK. As well the documentation has gone through
C grammatical changes.
C 881005 The prologue has been changed to mixed case.
C The subordinate routines had revision dates changed to
C this date, although the documentation for these routines
C is all upper case. No code changes.
C 890511 Code changes made. The DATA statement in the declaration
C section of DDASSL was replaced with a PARAMETER
C statement. Also the statement S = 100.D0 was removed
C from the top of the Newton iteration in DDASTP.
C The subordinate routines had revision dates changed to
C this date.
C 890517 The revision date syntax was replaced with the revision
C history syntax. Also the "DECK" comment was added to
C the top of all subroutines. These changes are consistent
C with new SLATEC guidelines.
C The subordinate routines had revision dates changed to
C this date. No code changes.
C 891013 Code changes made.
C Removed all occurrences of FLOAT or DBLE. All operations
C are now performed with "mixed-mode" arithmetic.
C Also, specific function names were replaced with generic
C function names to be consistent with new SLATEC guidelines.
C In particular:
C Replaced DSQRT with SQRT everywhere.
C Replaced DABS with ABS everywhere.
C Replaced DMIN1 with MIN everywhere.
C Replaced MIN0 with MIN everywhere.
C Replaced DMAX1 with MAX everywhere.
C Replaced MAX0 with MAX everywhere.
C Replaced DSIGN with SIGN everywhere.
C Also replaced REVISION DATE with REVISION HISTORY in all
C subordinate routines.
C 901004 Miscellaneous changes to prologue to complete conversion
C to SLATEC 4.0 format. No code changes. (F.N.Fritsch)
C 901009 Corrected GAMS classification code and converted subsidiary
C routines to 4.0 format. No code changes. (F.N.Fritsch)
C 901010 Converted XERRWV calls to XERMSG calls. (R.Clemens, AFWL)
C 901019 Code changes made.
C Merged SLATEC 4.0 changes with previous changes made
C by C. Ulrich. Below is a history of the changes made by
C C. Ulrich. (Changes in subsidiary routines are implied
C by this history)
C 891228 Bug was found and repaired inside the DDASSL
C and DDAINI routines. DDAINI was incorrectly
C returning the initial T with Y and YPRIME
C computed at T+H. The routine now returns T+H
C rather than the initial T.
C Cosmetic changes made to DDASTP.
C 900904 Three modifications were made to fix a bug (inside
C DDASSL) re interpolation for continuation calls and
C cases where TN is very close to TSTOP:
C
C 1) In testing for whether H is too large, just
C compare H to (TSTOP - TN), rather than
C (TSTOP - TN) * (1-4*UROUND), and set H to
C TSTOP - TN. This will force DDASTP to step
C exactly to TSTOP under certain situations
C (i.e. when H returned from DDASTP would otherwise
C take TN beyond TSTOP).
C
C 2) Inside the DDASTP loop, interpolate exactly to
C TSTOP if TN is very close to TSTOP (rather than
C interpolating to within roundoff of TSTOP).
C
C 3) Modified IDID description for IDID = 2 to say
C that the solution is returned by stepping exactly
C to TSTOP, rather than TOUT. (In some cases the
C solution is actually obtained by extrapolating
C over a distance near unit roundoff to TSTOP,
C but this small distance is deemed acceptable in
C these circumstances.)
C 901026 Added explicit declarations for all variables and minor
C cosmetic changes to prologue, removed unreferenced labels,
C and improved XERMSG calls. (FNF)
C 901030 Added ERROR MESSAGES section and reworked other sections to
C be of more uniform format. (FNF)
C 910624 Fixed minor bug related to HMAX (six lines after label
C 525). (LRP)
C***END PROLOGUE DDASSL