*DECK DDEBDF SUBROUTINE DDEBDF (DF, NEQ, T, Y, TOUT, INFO, RTOL, ATOL, IDID, + RWORK, LRW, IWORK, LIW, RPAR, IPAR, DJAC) C***BEGIN PROLOGUE DDEBDF C***PURPOSE Solve an initial value problem in ordinary differential C equations using backward differentiation formulas. It is C intended primarily for stiff problems. C***LIBRARY SLATEC (DEPAC) C***CATEGORY I1A2 C***TYPE DOUBLE PRECISION (DEBDF-S, DDEBDF-D) C***KEYWORDS BACKWARD DIFFERENTIATION FORMULAS, DEPAC, C INITIAL VALUE PROBLEMS, ODE, C ORDINARY DIFFERENTIAL EQUATIONS, STIFF C***AUTHOR Shampine, L. F., (SNLA) C Watts, H. A., (SNLA) C***DESCRIPTION C C This is the backward differentiation code in the package of C differential equation solvers DEPAC, consisting of the codes C DDERKF, DDEABM, and DDEBDF. Design of the package was by C L. F. Shampine and H. A. Watts. It is documented in C SAND-79-2374 , DEPAC - Design of a User Oriented Package of ODE C Solvers. C DDEBDF is a driver for a modification of the code LSODE written by C A. C. Hindmarsh C Lawrence Livermore Laboratory C Livermore, California 94550 C C ********************************************************************** C ** DEPAC PACKAGE OVERVIEW ** C ********************************************************************** C C You have a choice of three differential equation solvers from C DEPAC. The following brief descriptions are meant to aid you C in choosing the most appropriate code for your problem. C C DDERKF is a fifth order Runge-Kutta code. It is the simplest of C the three choices, both algorithmically and in the use of the C code. DDERKF is primarily designed to solve non-stiff and mild- C ly stiff differential equations when derivative evaluations are C not expensive. It should generally not be used to get high C accuracy results nor answers at a great many specific points. C Because DDERKF has very low overhead costs, it will usually C result in the least expensive integration when solving C problems requiring a modest amount of accuracy and having C equations that are not costly to evaluate. DDERKF attempts to C discover when it is not suitable for the task posed. C C DDEABM is a variable order (one through twelve) Adams code. Its C complexity lies somewhere between that of DDERKF and DDEBDF. C DDEABM is primarily designed to solve non-stiff and mildly C stiff differential equations when derivative evaluations are C expensive, high accuracy results are needed or answers at C many specific points are required. DDEABM attempts to discover C when it is not suitable for the task posed. C C DDEBDF is a variable order (one through five) backward C differentiation formula code. It is the most complicated of C the three choices. DDEBDF is primarily designed to solve stiff C differential equations at crude to moderate tolerances. C If the problem is very stiff at all, DDERKF and DDEABM will be C quite inefficient compared to DDEBDF. However, DDEBDF will be C inefficient compared to DDERKF and DDEABM on non-stiff problems C because it uses much more storage, has a much larger overhead, C and the low order formulas will not give high accuracies C efficiently. C C The concept of stiffness cannot be described in a few words. C If you do not know the problem to be stiff, try either DDERKF C or DDEABM. Both of these codes will inform you of stiffness C when the cost of solving such problems becomes important. C C ********************************************************************** C ** ABSTRACT ** C ********************************************************************** C C Subroutine DDEBDF uses the backward differentiation formulas of C orders one through five to integrate a system of NEQ first order C ordinary differential equations of the form C DU/DX = DF(X,U) C when the vector Y(*) of initial values for U(*) at X=T is given. C The subroutine integrates from T to TOUT. It is easy to continue the C integration to get results at additional TOUT. This is the interval C mode of operation. It is also easy for the routine to return with C the solution at each intermediate step on the way to TOUT. This is C the intermediate-output mode of operation. C C ********************************************************************** C * Description of The Arguments To DDEBDF (An Overview) * C ********************************************************************** C C The Parameters are: C C DF -- This is the name of a subroutine which you provide to C define the differential equations. C C NEQ -- This is the number of (first order) differential C equations to be integrated. C C T -- This is a DOUBLE PRECISION value of the independent C variable. C C Y(*) -- This DOUBLE PRECISION array contains the solution C components at T. C C TOUT -- This is a DOUBLE PRECISION point at which a solution is C desired. C C INFO(*) -- The basic task of the code is to integrate the C differential equations from T to TOUT and return an C answer at TOUT. INFO(*) is an INTEGER array which is used C to communicate exactly how you want this task to be C carried out. C C RTOL, ATOL -- These DOUBLE PRECISION quantities C represent relative and absolute error tolerances which you C provide to indicate how accurately you wish the solution C to be computed. You may choose them to be both scalars C or else both vectors. C C IDID -- This scalar quantity is an indicator reporting what C the code did. You must monitor this INTEGER variable to C decide what action to take next. C C RWORK(*), LRW -- RWORK(*) is a DOUBLE PRECISION work array of C length LRW which provides the code with needed storage C space. C C IWORK(*), LIW -- IWORK(*) is an INTEGER work array of length LIW C which provides the code with needed storage space and an C across call flag. C C RPAR, IPAR -- These are DOUBLE PRECISION and INTEGER parameter C arrays which you can use for communication between your C calling program and the DF subroutine (and the DJAC C subroutine). C C DJAC -- This is the name of a subroutine which you may choose to C provide for defining the Jacobian matrix of partial C derivatives DF/DU. C C Quantities which are used as input items are C NEQ, T, Y(*), TOUT, INFO(*), C RTOL, ATOL, RWORK(1), LRW, C IWORK(1), IWORK(2), and LIW. C C Quantities which may be altered by the code are C T, Y(*), INFO(1), RTOL, ATOL, C IDID, RWORK(*) and IWORK(*). C C ********************************************************************** C * INPUT -- What To Do On The First Call To DDEBDF * C ********************************************************************** 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 DF -- Provide a subroutine of the form C DF(X,U,UPRIME,RPAR,IPAR) C to define the system of first order differential equations C which is to be solved. For the given values of X and the C vector U(*)=(U(1),U(2),...,U(NEQ)) , the subroutine must C evaluate the NEQ components of the system of differential C equations DU/DX=DF(X,U) and store the derivatives in the C array UPRIME(*), that is, UPRIME(I) = * DU(I)/DX * for C equations I=1,...,NEQ. C C Subroutine DF must not alter X or U(*). You must declare C the name DF in an external statement in your program that C calls DDEBDF. You must dimension U and UPRIME in DF. C C RPAR and IPAR are DOUBLE PRECISION and INTEGER parameter C arrays which you can use for communication between your C calling program and subroutine DF. They are not used or C altered by DDEBDF. 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 them in C your calling program and in DF as arrays of appropriate C 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 You must use a program variable for T because the code C changes its value. C C Y(*) -- Set this vector to the initial values of the NEQ solution C components at the initial point. You must dimension Y at C least NEQ in your calling program. C C TOUT -- Set it to the first point at which a solution is desired. C You can take TOUT = T, in which case the code C will evaluate the derivative of the solution at T and C return. Integration either forward in T (TOUT .GT. T) C or 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 following C each intermediate step (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 C step 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 to accommodate other members of C DEPAC or possible future extensions, though DDEBDF uses C only the first six entries. You must respond to all of C the following items which are arranged as questions. The C simplest use of the code corresponds to answering all C questions as YES ,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 stiff problems it is necessary to use the C Jacobian matrix of partial derivatives of the system C of differential equations. If you do not provide a C subroutine to evaluate it analytically (see the C description of the item DJAC 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 DJAC. Sometimes numerical differencing C is cheaper than evaluating derivatives in DJAC and C sometimes it is not - this depends on your problem. C C If your problem is linear, i.e. has the form C DU/DX = DF(X,U) = J(X)*U + G(X) for some matrix J(X) C and vector G(X), the Jacobian matrix DF/DU = J(X). C Since you must provide a subroutine to evaluate DF(X,U) C analytically, it is little extra trouble to provide C subroutine DJAC for evaluating J(X) analytically. C Furthermore, in such cases, numerical differencing is C much more expensive than analytic evaluation. 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 DJAC for evaluating the C Jacobian matrix **** C C INFO(6) -- DDEBDF will perform much better if the Jacobian C matrix 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 more cheaply, 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 Jacobian, C the code works with a full matrix of NEQ**2 elements C (stored in the conventional way). Computations with C 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 Jacobian matrix has a banded structure and C you want to provide subroutine DJAC to compute the C partial derivatives, then you must be careful to store C the elements of the Jacobian matrix in the special form C indicated in the description of DJAC. C C **** Do you want to solve the problem using a full C (dense) Jacobian 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 RTOL, ATOL -- You must assign relative (RTOL) and absolute (ATOL) C error tolerances to tell the code how accurately you want C the solution to be computed. They must be defined as C program 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 test C 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 true C solution of the initial value problem and the computed C approximation. Practically all present day codes, C including this one, control the local error at each step C and do not even attempt to control the global error C directly. Roughly speaking, they produce a solution Y(T) C which satisfies the differential equations with a C residual R(T), DY(T)/DT = DF(T,Y(T)) + R(T) , C and, almost always, R(T) is bounded by the error C tolerances. Usually, but not always, the true accuracy of C the computed Y is comparable to the error tolerances. This C code will usually, but not always, deliver a more accurate C solution if you reduce the tolerances and integrate again. C By comparing two such solutions you can get a fairly C reliable idea of the true error in the solution at the C 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 abso- C lute error test on that component. A mixed test with non- C zero RTOL and ATOL corresponds roughly to a relative error C test when the solution component is much bigger than ATOL C and to an absolute error test when the solution component C is smaller than the threshold ATOL. C C Proper selection of the absolute error control parameters C ATOL requires you to have some idea of the scale of the C solution components. To acquire this information may mean C that you will have to solve the problem more than once. In C the absence of scale information, you should ask for some C relative accuracy in all the components (by setting RTOL C values non-zero) and perhaps impose extremely small C absolute error tolerances to protect against the danger of C a solution component becoming zero. 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 DOUBLE PRECISION work array of length C LRW in your calling program. C C RWORK(1) -- If you have set INFO(4)=0, you can ignore this C optional input parameter. Otherwise you must define a C stopping point TSTOP by setting RWORK(1) = TSTOP. 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.) C C LRW -- Set it to the declared length of the RWORK array. C You must have C LRW .GE. 250+10*NEQ+NEQ**2 C for the full (dense) Jacobian case (when INFO(6)=0), or C LRW .GE. 250+10*NEQ+(2*ML+MU+1)*NEQ C for the banded Jacobian case (when INFO(6)=1). C C IWORK(*) -- Dimension this INTEGER work array of length LIW in C your calling program. C C IWORK(1), IWORK(2) -- If you have set INFO(6)=0, you can ignore C these optional input parameters. Otherwise you must define C the half-bandwidths ML (lower) and MU (upper) of the C Jacobian matrix by setting IWORK(1) = ML and C IWORK(2) = MU. (The code will work with a full matrix C of NEQ**2 elements unless it is told that the problem has C a banded Jacobian, in which case the code will work with C a matrix containing at most (2*ML+MU+1)*NEQ elements.) C C LIW -- Set it to the declared length of the IWORK array. C You must have LIW .GE. 56+NEQ. C C RPAR, IPAR -- These are parameter arrays, of DOUBLE PRECISION and C INTEGER type, respectively. You can use them for C communication between your program that calls DDEBDF and C the DF subroutine (and the DJAC subroutine). They are not C used or altered by DDEBDF. If you do not need RPAR or C IPAR, ignore these parameters by treating them as dummy C arguments. If you do choose to use them, dimension them in C your calling program and in DF (and in DJAC) as arrays of C appropriate length. C C DJAC -- If you have set INFO(5)=0, you can ignore this parameter C by treating it as a dummy argument. (For some compilers C you may have to write a dummy subroutine named DJAC in C order to avoid problems associated with missing external C routine names.) Otherwise, you must provide a subroutine C of the form C DJAC(X,U,PD,NROWPD,RPAR,IPAR) C to define the Jacobian matrix of partial derivatives DF/DU C of the system of differential equations DU/DX = DF(X,U). C For the given values of X and the vector C U(*)=(U(1),U(2),...,U(NEQ)), the subroutine must evaluate C the non-zero partial derivatives DF(I)/DU(J) for each C differential equation I=1,...,NEQ and each solution C component J=1,...,NEQ , and store these values in the C matrix PD. The elements of PD are set to zero before each C call to DJAC so only non-zero elements need to be defined. C C Subroutine DJAC must not alter X, U(*), or NROWPD. You C must declare the name DJAC in an external statement in C your program that calls DDEBDF. NROWPD is the row C dimension of the PD matrix and is assigned by the code. C Therefore you must dimension PD in DJAC according to C DIMENSION PD(NROWPD,1) C You must also dimension U in DJAC. C C The way you must store the elements into the PD matrix C depends on the structure of the Jacobian which you C indicated by INFO(6). C *** INFO(6)=0 -- Full (Dense) Jacobian *** 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) = * DF(I)/DU(J) * C *** INFO(6)=1 -- Banded Jacobian with ML Lower and MU C Upper Diagonal Bands (refer to INFO(6) description of C ML and MU) *** 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) = * DF(I)/DU(J) * C C RPAR and IPAR are DOUBLE PRECISION and INTEGER parameter C arrays which you can use for communication between your C calling program and your Jacobian subroutine DJAC. They C are not altered by DDEBDF. If you do not need RPAR or C IPAR, ignore these parameters by treating them as dummy C arguments. If you do choose to use them, dimension them C in your calling program and in DJAC as arrays of C appropriate length. C C ********************************************************************** C * OUTPUT -- After any return from DDEBDF * C ********************************************************************** 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 You may also be interested in the approximate derivative C of the solution at T. It is contained in C RWORK(21),...,RWORK(20+NEQ). 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 TOUT was successfully C completed (T=TOUT) by stepping exactly to TOUT. 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 C *** Task Interrupted *** C Reported by negative values of IDID C C IDID = -1 -- A large amount of work has been expended. C (500 steps attempted) 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 = -4,-5 -- Not applicable for this code but used C by other members of DEPAC. C C IDID = -6 -- DDEBDF had repeated convergence test failures C on the last attempted step. C C IDID = -7 -- DDEBDF had repeated error test failures on C the last attempted step. C C IDID = -8,..,-32 -- Not applicable for this code but C used by other members of DEPAC or possible C future extensions. 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 C occurs 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 be C appropriate for continuing the integration. However, the C reported solution at T was obtained using the input values C 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(11)--which contains the step size H to be C attempted on the next step. C C RWORK(12)--If the tolerances have been increased by the C code (IDID = -2) , they were multiplied by the C value in RWORK(12). C C RWORK(13)--which contains the current value of the C independent variable, i.e. the farthest point C integration has reached. This will be C different from T only when interpolation has C been performed (IDID=3). C C RWORK(20+I)--which contains the approximate derivative C of the solution component Y(I). In DDEBDF, it C is never obtained by calling subroutine DF to C evaluate the differential equation using T and C Y(*), except at the initial point of C integration. C C ********************************************************************** C ** INPUT -- What To Do To Continue The Integration ** C ** (calls after the first) ** C ********************************************************************** 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(*), RWORK(*), IWORK(*) or C the differential equation in subroutine DF. Any such alteration C constitutes a new problem and must be treated as such, i.e. C 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 The parameter INFO(1) is used by the code to indicate the C beginning of a new problem and to indicate whether integration C is to be continued. You must input the value INFO(1) = 0 C when starting a new problem. You must input the value C INFO(1) = 1 if you wish to continue after an interrupted task. C Do not set INFO(1) = 0 on a continuation call unless you C want the code to restart at the current T. 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 reset INFO(1) = 1 C If C IDID = -1, the code has attempted 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 but used C by other members of DEPAC. C C IDID = -6, repeated convergence test failures occurred C on the last attempted step in DDEBDF. An inaccu- C rate Jacobian may be the problem. If you are C absolutely certain you want to continue, restart C the integration at the current T by setting C INFO(1)=0 and call the code again. C C IDID = -7, repeated error test failures occurred on the C last attempted step in DDEBDF. A singularity in C the solution may be present. You should re- C examine the problem being solved. If you are C absolutely certain you want to continue, restart C the integration at the current T by setting C INFO(1)=0 and call the code again. C C IDID = -8,..,-32 --- cannot occur with this code but C used by other members of DDEPAC or possible future C extensions. C C *** Following a Terminated Task *** C If C IDID = -33, you cannot continue the solution of this C problem. An attempt to do so will result in your C run being terminated. C C ********************************************************************** C C ***** Warning ***** C C If DDEBDF is to be used in an overlay situation, you must save and C restore certain items used internally by DDEBDF (values in the C common block DDEBD1). This can be accomplished as follows. C C To save the necessary values upon return from DDEBDF, simply call C DSVCO(RWORK(22+NEQ),IWORK(21+NEQ)). C C To restore the necessary values before the next call to DDEBDF, C simply call DRSCO(RWORK(22+NEQ),IWORK(21+NEQ)). C C***REFERENCES L. F. Shampine and H. A. Watts, DEPAC - design of a user C oriented package of ODE solvers, Report SAND79-2374, C Sandia Laboratories, 1979. C***ROUTINES CALLED DLSOD, XERMSG C***COMMON BLOCKS DDEBD1 C***REVISION HISTORY (YYMMDD) C 820301 DATE WRITTEN C 890831 Modified array declarations. (WRB) C 891024 Changed references from DVNORM to DHVNRM. (WRB) C 891024 REVISION DATE from Version 3.2 C 891214 Prologue converted to Version 4.0 format. (BAB) C 900326 Removed duplicate information from DESCRIPTION section. C (WRB) C 900510 Convert XERRWV calls to XERMSG calls, make Prologue comments C consistent with DEBDF. (RWC) C 920501 Reformatted the REFERENCES section. (WRB) C***END PROLOGUE DDEBDF INTEGER IACOR, IBAND, IBEGIN, ICOMI, ICOMR, IDELSN, IDID, IER, 1 IEWT, IINOUT, IINTEG, IJAC, ILRW, INFO, INIT, 2 IOWNS, IPAR, IQUIT, ISAVF, ITOL, ITSTAR, ITSTOP, IWM, 3 IWORK, IYH, IYPOUT, JSTART, KFLAG, KSTEPS, L, LIW, LRW, 4 MAXORD, METH, MITER, ML, MU, N, NEQ, NFE, NJE, NQ, NQU, 5 NST DOUBLE PRECISION ATOL, EL0, H, HMIN, HMXI, HU, ROWNS, RPAR, 1 RTOL, RWORK, T, TN, TOLD, TOUT, UROUND, Y LOGICAL INTOUT CHARACTER*8 XERN1, XERN2 CHARACTER*16 XERN3 C DIMENSION Y(*),INFO(15),RTOL(*),ATOL(*),RWORK(*),IWORK(*), 1 RPAR(*),IPAR(*) C COMMON /DDEBD1/ TOLD,ROWNS(210),EL0,H,HMIN,HMXI,HU,TN,UROUND, 1 IQUIT,INIT,IYH,IEWT,IACOR,ISAVF,IWM,KSTEPS,IBEGIN, 2 ITOL,IINTEG,ITSTOP,IJAC,IBAND,IOWNS(6),IER,JSTART, 3 KFLAG,L,METH,MITER,MAXORD,N,NQ,NST,NFE,NJE,NQU C EXTERNAL DF, DJAC C C CHECK FOR AN APPARENT INFINITE LOOP C C***FIRST EXECUTABLE STATEMENT DDEBDF IF (INFO(1) .EQ. 0) IWORK(LIW) = 0 C IF (IWORK(LIW).GE. 5) THEN IF (T .EQ. RWORK(21+NEQ)) THEN WRITE (XERN3, '(1PE15.6)') T CALL XERMSG ('SLATEC', 'DDEBDF', * 'AN APPARENT INFINITE LOOP HAS BEEN DETECTED.$$' // * 'YOU HAVE MADE REPEATED CALLS AT T = ' // XERN3 // * ' AND THE INTEGRATION HAS NOT ADVANCED. CHECK THE ' // * 'WAY YOU HAVE SET PARAMETERS FOR THE CALL TO THE ' // * 'CODE, PARTICULARLY INFO(1).', 13, 2) RETURN ENDIF ENDIF C IDID = 0 C C CHECK VALIDITY OF INFO PARAMETERS C IF (INFO(1) .NE. 0 .AND. INFO(1) .NE. 1) THEN WRITE (XERN1, '(I8)') INFO(1) CALL XERMSG ('SLATEC', 'DDEBDF', 'INFO(1) MUST BE SET TO 0 ' // * 'FOR THE START OF A NEW PROBLEM, AND MUST BE SET TO 1 ' // * 'FOLLOWING AN INTERRUPTED TASK. YOU ARE ATTEMPTING TO ' // * 'CONTINUE THE INTEGRATION ILLEGALLY BY CALLING THE ' // * 'CODE WITH INFO(1) = ' // XERN1, 3, 1) IDID = -33 ENDIF C IF (INFO(2) .NE. 0 .AND. INFO(2) .NE. 1) THEN WRITE (XERN1, '(I8)') INFO(2) CALL XERMSG ('SLATEC', 'DDEBDF', 'INFO(2) MUST BE 0 OR 1 ' // * 'INDICATING SCALAR AND VECTOR ERROR TOLERANCES, ' // * 'RESPECTIVELY. YOU HAVE CALLED THE CODE WITH INFO(2) = ' // * XERN1, 4, 1) IDID = -33 ENDIF C IF (INFO(3) .NE. 0 .AND. INFO(3) .NE. 1) THEN WRITE (XERN1, '(I8)') INFO(3) CALL XERMSG ('SLATEC', 'DDEBDF', 'INFO(3) MUST BE 0 OR 1 ' // * 'INDICATING THE INTERVAL OR INTERMEDIATE-OUTPUT MODE OF ' // * 'INTEGRATION, RESPECTIVELY. YOU HAVE CALLED THE CODE ' // * 'WITH INFO(3) = ' // XERN1, 5, 1) IDID = -33 ENDIF C IF (INFO(4) .NE. 0 .AND. INFO(4) .NE. 1) THEN WRITE (XERN1, '(I8)') INFO(4) CALL XERMSG ('SLATEC', 'DDEBDF', 'INFO(4) MUST BE 0 OR 1 ' // * 'INDICATING WHETHER OR NOT THE INTEGRATION INTERVAL IS ' // * 'TO BE RESTRICTED BY A POINT TSTOP. YOU HAVE CALLED ' // * 'THE CODE WITH INFO(4) = ' // XERN1, 14, 1) IDID = -33 ENDIF C IF (INFO(5) .NE. 0 .AND. INFO(5) .NE. 1) THEN WRITE (XERN1, '(I8)') INFO(5) CALL XERMSG ('SLATEC', 'DDEBDF', 'INFO(5) MUST BE 0 OR 1 ' // * 'INDICATING WHETHER THE CODE IS TOLD TO FORM THE ' // * 'JACOBIAN MATRIX BY NUMERICAL DIFFERENCING OR YOU ' // * 'PROVIDE A SUBROUTINE TO EVALUATE IT ANALYTICALLY. ' // * 'YOU HAVE CALLED THE CODE WITH INFO(5) = ' // XERN1, 15, 1) IDID = -33 ENDIF C IF (INFO(6) .NE. 0 .AND. INFO(6) .NE. 1) THEN WRITE (XERN1, '(I8)') INFO(6) CALL XERMSG ('SLATEC', 'DDEBDF', 'INFO(6) MUST BE 0 OR 1 ' // * 'INDICATING WHETHER THE CODE IS TOLD TO TREAT THE ' // * 'JACOBIAN AS A FULL (DENSE) MATRIX OR AS HAVING A ' // * 'SPECIAL BANDED STRUCTURE. YOU HAVE CALLED THE CODE ' // * 'WITH INFO(6) = ' // XERN1, 16, 1) IDID = -33 ENDIF C ILRW = NEQ IF (INFO(6) .NE. 0) THEN C C CHECK BANDWIDTH PARAMETERS C ML = IWORK(1) MU = IWORK(2) ILRW = 2*ML + MU + 1 C IF (ML.LT.0 .OR. ML.GE.NEQ .OR. MU.LT.0 .OR. MU.GE.NEQ) THEN WRITE (XERN1, '(I8)') ML WRITE (XERN2, '(I8)') MU CALL XERMSG ('SLATEC', 'DDEBDF', 'YOU HAVE SET INFO(6) ' // * '= 1, TELLING THE CODE THAT THE JACOBIAN MATRIX HAS ' // * 'A SPECIAL BANDED STRUCTURE. HOWEVER, THE LOWER ' // * '(UPPER) BANDWIDTHS ML (MU) VIOLATE THE CONSTRAINTS ' // * 'ML,MU .GE. 0 AND ML,MU .LT. NEQ. YOU HAVE CALLED ' // * 'THE CODE WITH ML = ' // XERN1 // ' AND MU = ' // XERN2, * 17, 1) IDID = -33 ENDIF ENDIF C C CHECK LRW AND LIW FOR SUFFICIENT STORAGE ALLOCATION C IF (LRW .LT. 250 + (10 + ILRW)*NEQ) THEN WRITE (XERN1, '(I8)') LRW IF (INFO(6) .EQ. 0) THEN CALL XERMSG ('SLATEC', 'DDEBDF', 'LENGTH OF ARRAY RWORK ' // * 'MUST BE AT LEAST 250 + 10*NEQ + NEQ*NEQ.$$' // * 'YOU HAVE CALLED THE CODE WITH LRW = ' // XERN1, 1, 1) ELSE CALL XERMSG ('SLATEC', 'DDEBDF', 'LENGTH OF ARRAY RWORK ' // * 'MUST BE AT LEAST 250 + 10*NEQ + (2*ML+MU+1)*NEQ.$$' // * 'YOU HAVE CALLED THE CODE WITH LRW = ' // XERN1, 18, 1) ENDIF IDID = -33 ENDIF C IF (LIW .LT. 56 + NEQ) THEN WRITE (XERN1, '(I8)') LIW CALL XERMSG ('SLATEC', 'DDEBDF', 'LENGTH OF ARRAY IWORK ' // * 'BE AT LEAST 56 + NEQ. YOU HAVE CALLED THE CODE WITH ' // * 'LIW = ' // XERN1, 2, 1) IDID = -33 ENDIF C C COMPUTE THE INDICES FOR THE ARRAYS TO BE STORED IN THE WORK C ARRAY AND RESTORE COMMON BLOCK DATA C ICOMI = 21 + NEQ IINOUT = ICOMI + 33 C IYPOUT = 21 ITSTAR = 21 + NEQ ICOMR = 22 + NEQ C IF (INFO(1) .NE. 0) INTOUT = IWORK(IINOUT) .NE. (-1) C CALL DRSCO(RWORK(ICOMR),IWORK(ICOMI)) C IYH = ICOMR + 218 IEWT = IYH + 6*NEQ ISAVF = IEWT + NEQ IACOR = ISAVF + NEQ IWM = IACOR + NEQ IDELSN = IWM + 2 + ILRW*NEQ C IBEGIN = INFO(1) ITOL = INFO(2) IINTEG = INFO(3) ITSTOP = INFO(4) IJAC = INFO(5) IBAND = INFO(6) RWORK(ITSTAR) = T C CALL DLSOD(DF,NEQ,T,Y,TOUT,RTOL,ATOL,IDID,RWORK(IYPOUT), 1 RWORK(IYH),RWORK(IYH),RWORK(IEWT),RWORK(ISAVF), 2 RWORK(IACOR),RWORK(IWM),IWORK(1),DJAC,INTOUT, 3 RWORK(1),RWORK(12),RWORK(IDELSN),RPAR,IPAR) C IWORK(IINOUT) = -1 IF (INTOUT) IWORK(IINOUT) = 1 C IF (IDID .NE. (-2)) IWORK(LIW) = IWORK(LIW) + 1 IF (T .NE. RWORK(ITSTAR)) IWORK(LIW) = 0 C CALL DSVCO(RWORK(ICOMR),IWORK(ICOMI)) RWORK(11) = H RWORK(13) = TN INFO(1) = IBEGIN C RETURN END