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