00001 SUBROUTINE DDASSL (RES, NEQ, T, Y, YPRIME, TOUT, INFO, RTOL, ATOL, 00002 + IDID, RWORK, LRW, IWORK, LIW, RPAR, IPAR, JAC) 00003 C***BEGIN PROLOGUE DDASSL 00004 C***PURPOSE This code solves a system of differential/algebraic 00005 C equations of the form G(T,Y,YPRIME) = 0. 00006 C***LIBRARY SLATEC (DASSL) 00007 C***CATEGORY I1A2 00008 C***TYPE DOUBLE PRECISION (SDASSL-S, DDASSL-D) 00009 C***KEYWORDS DIFFERENTIAL/ALGEBRAIC, BACKWARD DIFFERENTIATION FORMULAS, 00010 C IMPLICIT DIFFERENTIAL SYSTEMS 00011 C***AUTHOR PETZOLD, LINDA R., (LLNL) 00012 C COMPUTING AND MATHEMATICS RESEARCH DIVISION 00013 C LAWRENCE LIVERMORE NATIONAL LABORATORY 00014 C L - 316, P.O. BOX 808, 00015 C LIVERMORE, CA. 94550 00016 C***DESCRIPTION 00017 C 00018 C *Usage: 00019 C 00020 C EXTERNAL RES, JAC 00021 C INTEGER NEQ, INFO(N), IDID, LRW, LIW, IWORK(LIW), IPAR 00022 C DOUBLE PRECISION T, Y(NEQ), YPRIME(NEQ), TOUT, RTOL, ATOL, 00023 C * RWORK(LRW), RPAR 00024 C 00025 C CALL DDASSL (RES, NEQ, T, Y, YPRIME, TOUT, INFO, RTOL, ATOL, 00026 C * IDID, RWORK, LRW, IWORK, LIW, RPAR, IPAR, JAC) 00027 C 00028 C 00029 C *Arguments: 00030 C (In the following, all real arrays should be type DOUBLE PRECISION.) 00031 C 00032 C RES:EXT This is a subroutine which you provide to define the 00033 C differential/algebraic system. 00034 C 00035 C NEQ:IN This is the number of equations to be solved. 00036 C 00037 C T:INOUT This is the current value of the independent variable. 00038 C 00039 C Y(*):INOUT This array contains the solution components at T. 00040 C 00041 C YPRIME(*):INOUT This array contains the derivatives of the solution 00042 C components at T. 00043 C 00044 C TOUT:IN This is a point at which a solution is desired. 00045 C 00046 C INFO(N):IN The basic task of the code is to solve the system from T 00047 C to TOUT and return an answer at TOUT. INFO is an integer 00048 C array which is used to communicate exactly how you want 00049 C this task to be carried out. (See below for details.) 00050 C N must be greater than or equal to 15. 00051 C 00052 C RTOL,ATOL:INOUT These quantities represent relative and absolute 00053 C error tolerances which you provide to indicate how 00054 C accurately you wish the solution to be computed. You 00055 C may choose them to be both scalars or else both vectors. 00056 C Caution: In Fortran 77, a scalar is not the same as an 00057 C array of length 1. Some compilers may object 00058 C to using scalars for RTOL,ATOL. 00059 C 00060 C IDID:OUT This scalar quantity is an indicator reporting what the 00061 C code did. You must monitor this integer variable to 00062 C decide what action to take next. 00063 C 00064 C RWORK:WORK A real work array of length LRW which provides the 00065 C code with needed storage space. 00066 C 00067 C LRW:IN The length of RWORK. (See below for required length.) 00068 C 00069 C IWORK:WORK An integer work array of length LIW which probides the 00070 C code with needed storage space. 00071 C 00072 C LIW:IN The length of IWORK. (See below for required length.) 00073 C 00074 C RPAR,IPAR:IN These are real and integer parameter arrays which 00075 C you can use for communication between your calling 00076 C program and the RES subroutine (and the JAC subroutine) 00077 C 00078 C JAC:EXT This is the name of a subroutine which you may choose 00079 C to provide for defining a matrix of partial derivatives 00080 C described below. 00081 C 00082 C Quantities which may be altered by DDASSL are: 00083 C T, Y(*), YPRIME(*), INFO(1), RTOL, ATOL, 00084 C IDID, RWORK(*) AND IWORK(*) 00085 C 00086 C *Description 00087 C 00088 C Subroutine DDASSL uses the backward differentiation formulas of 00089 C orders one through five to solve a system of the above form for Y and 00090 C YPRIME. Values for Y and YPRIME at the initial time must be given as 00091 C input. These values must be consistent, (that is, if T,Y,YPRIME are 00092 C the given initial values, they must satisfy G(T,Y,YPRIME) = 0.). The 00093 C subroutine solves the system from T to TOUT. It is easy to continue 00094 C the solution to get results at additional TOUT. This is the interval 00095 C mode of operation. Intermediate results can also be obtained easily 00096 C by using the intermediate-output capability. 00097 C 00098 C The following detailed description is divided into subsections: 00099 C 1. Input required for the first call to DDASSL. 00100 C 2. Output after any return from DDASSL. 00101 C 3. What to do to continue the integration. 00102 C 4. Error messages. 00103 C 00104 C 00105 C -------- INPUT -- WHAT TO DO ON THE FIRST CALL TO DDASSL ------------ 00106 C 00107 C The first call of the code is defined to be the start of each new 00108 C problem. Read through the descriptions of all the following items, 00109 C provide sufficient storage space for designated arrays, set 00110 C appropriate variables for the initialization of the problem, and 00111 C give information about how you want the problem to be solved. 00112 C 00113 C 00114 C RES -- Provide a subroutine of the form 00115 C SUBROUTINE RES(T,Y,YPRIME,DELTA,IRES,RPAR,IPAR) 00116 C to define the system of differential/algebraic 00117 C equations which is to be solved. For the given values 00118 C of T,Y and YPRIME, the subroutine should 00119 C return the residual of the defferential/algebraic 00120 C system 00121 C DELTA = G(T,Y,YPRIME) 00122 C (DELTA(*) is a vector of length NEQ which is 00123 C output for RES.) 00124 C 00125 C Subroutine RES must not alter T,Y or YPRIME. 00126 C You must declare the name RES in an external 00127 C statement in your program that calls DDASSL. 00128 C You must dimension Y,YPRIME and DELTA in RES. 00129 C 00130 C IRES is an integer flag which is always equal to 00131 C zero on input. Subroutine RES should alter IRES 00132 C only if it encounters an illegal value of Y or 00133 C a stop condition. Set IRES = -1 if an input value 00134 C is illegal, and DDASSL will try to solve the problem 00135 C without getting IRES = -1. If IRES = -2, DDASSL 00136 C will return control to the calling program 00137 C with IDID = -11. 00138 C 00139 C RPAR and IPAR are real and integer parameter arrays which 00140 C you can use for communication between your calling program 00141 C and subroutine RES. They are not altered by DDASSL. If you 00142 C do not need RPAR or IPAR, ignore these parameters by treat- 00143 C ing them as dummy arguments. If you do choose to use them, 00144 C dimension them in your calling program and in RES as arrays 00145 C of appropriate length. 00146 C 00147 C NEQ -- Set it to the number of differential equations. 00148 C (NEQ .GE. 1) 00149 C 00150 C T -- Set it to the initial point of the integration. 00151 C T must be defined as a variable. 00152 C 00153 C Y(*) -- Set this vector to the initial values of the NEQ solution 00154 C components at the initial point. You must dimension Y of 00155 C length at least NEQ in your calling program. 00156 C 00157 C YPRIME(*) -- Set this vector to the initial values of the NEQ 00158 C first derivatives of the solution components at the initial 00159 C point. You must dimension YPRIME at least NEQ in your 00160 C calling program. If you do not know initial values of some 00161 C of the solution components, see the explanation of INFO(11). 00162 C 00163 C TOUT -- Set it to the first point at which a solution 00164 C is desired. You can not take TOUT = T. 00165 C integration either forward in T (TOUT .GT. T) or 00166 C backward in T (TOUT .LT. T) is permitted. 00167 C 00168 C The code advances the solution from T to TOUT using 00169 C step sizes which are automatically selected so as to 00170 C achieve the desired accuracy. If you wish, the code will 00171 C return with the solution and its derivative at 00172 C intermediate steps (intermediate-output mode) so that 00173 C you can monitor them, but you still must provide TOUT in 00174 C accord with the basic aim of the code. 00175 C 00176 C The first step taken by the code is a critical one 00177 C because it must reflect how fast the solution changes near 00178 C the initial point. The code automatically selects an 00179 C initial step size which is practically always suitable for 00180 C the problem. By using the fact that the code will not step 00181 C past TOUT in the first step, you could, if necessary, 00182 C restrict the length of the initial step size. 00183 C 00184 C For some problems it may not be permissible to integrate 00185 C past a point TSTOP because a discontinuity occurs there 00186 C or the solution or its derivative is not defined beyond 00187 C TSTOP. When you have declared a TSTOP point (SEE INFO(4) 00188 C and RWORK(1)), you have told the code not to integrate 00189 C past TSTOP. In this case any TOUT beyond TSTOP is invalid 00190 C input. 00191 C 00192 C INFO(*) -- Use the INFO array to give the code more details about 00193 C how you want your problem solved. This array should be 00194 C dimensioned of length 15, though DDASSL uses only the first 00195 C eleven entries. You must respond to all of the following 00196 C items, which are arranged as questions. The simplest use 00197 C of the code corresponds to answering all questions as yes, 00198 C i.e. setting all entries of INFO to 0. 00199 C 00200 C INFO(1) - This parameter enables the code to initialize 00201 C itself. You must set it to indicate the start of every 00202 C new problem. 00203 C 00204 C **** Is this the first call for this problem ... 00205 C Yes - Set INFO(1) = 0 00206 C No - Not applicable here. 00207 C See below for continuation calls. **** 00208 C 00209 C INFO(2) - How much accuracy you want of your solution 00210 C is specified by the error tolerances RTOL and ATOL. 00211 C The simplest use is to take them both to be scalars. 00212 C To obtain more flexibility, they can both be vectors. 00213 C The code must be told your choice. 00214 C 00215 C **** Are both error tolerances RTOL, ATOL scalars ... 00216 C Yes - Set INFO(2) = 0 00217 C and input scalars for both RTOL and ATOL 00218 C No - Set INFO(2) = 1 00219 C and input arrays for both RTOL and ATOL **** 00220 C 00221 C INFO(3) - The code integrates from T in the direction 00222 C of TOUT by steps. If you wish, it will return the 00223 C computed solution and derivative at the next 00224 C intermediate step (the intermediate-output mode) or 00225 C TOUT, whichever comes first. This is a good way to 00226 C proceed if you want to see the behavior of the solution. 00227 C If you must have solutions at a great many specific 00228 C TOUT points, this code will compute them efficiently. 00229 C 00230 C **** Do you want the solution only at 00231 C TOUT (and not at the next intermediate step) ... 00232 C Yes - Set INFO(3) = 0 00233 C No - Set INFO(3) = 1 **** 00234 C 00235 C INFO(4) - To handle solutions at a great many specific 00236 C values TOUT efficiently, this code may integrate past 00237 C TOUT and interpolate to obtain the result at TOUT. 00238 C Sometimes it is not possible to integrate beyond some 00239 C point TSTOP because the equation changes there or it is 00240 C not defined past TSTOP. Then you must tell the code 00241 C not to go past. 00242 C 00243 C **** Can the integration be carried out without any 00244 C restrictions on the independent variable T ... 00245 C Yes - Set INFO(4)=0 00246 C No - Set INFO(4)=1 00247 C and define the stopping point TSTOP by 00248 C setting RWORK(1)=TSTOP **** 00249 C 00250 C INFO(5) - To solve differential/algebraic problems it is 00251 C necessary to use a matrix of partial derivatives of the 00252 C system of differential equations. If you do not 00253 C provide a subroutine to evaluate it analytically (see 00254 C description of the item JAC in the call list), it will 00255 C be approximated by numerical differencing in this code. 00256 C although it is less trouble for you to have the code 00257 C compute partial derivatives by numerical differencing, 00258 C the solution will be more reliable if you provide the 00259 C derivatives via JAC. Sometimes numerical differencing 00260 C is cheaper than evaluating derivatives in JAC and 00261 C sometimes it is not - this depends on your problem. 00262 C 00263 C **** Do you want the code to evaluate the partial 00264 C derivatives automatically by numerical differences ... 00265 C Yes - Set INFO(5)=0 00266 C No - Set INFO(5)=1 00267 C and provide subroutine JAC for evaluating the 00268 C matrix of partial derivatives **** 00269 C 00270 C INFO(6) - DDASSL will perform much better if the matrix of 00271 C partial derivatives, DG/DY + CJ*DG/DYPRIME, 00272 C (here CJ is a scalar determined by DDASSL) 00273 C is banded and the code is told this. In this 00274 C case, the storage needed will be greatly reduced, 00275 C numerical differencing will be performed much cheaper, 00276 C and a number of important algorithms will execute much 00277 C faster. The differential equation is said to have 00278 C half-bandwidths ML (lower) and MU (upper) if equation i 00279 C involves only unknowns Y(J) with 00280 C I-ML .LE. J .LE. I+MU 00281 C for all I=1,2,...,NEQ. Thus, ML and MU are the widths 00282 C of the lower and upper parts of the band, respectively, 00283 C with the main diagonal being excluded. If you do not 00284 C indicate that the equation has a banded matrix of partial 00285 C derivatives, the code works with a full matrix of NEQ**2 00286 C elements (stored in the conventional way). Computations 00287 C with banded matrices cost less time and storage than with 00288 C full matrices if 2*ML+MU .LT. NEQ. If you tell the 00289 C code that the matrix of partial derivatives has a banded 00290 C structure and you want to provide subroutine JAC to 00291 C compute the partial derivatives, then you must be careful 00292 C to store the elements of the matrix in the special form 00293 C indicated in the description of JAC. 00294 C 00295 C **** Do you want to solve the problem using a full 00296 C (dense) matrix (and not a special banded 00297 C structure) ... 00298 C Yes - Set INFO(6)=0 00299 C No - Set INFO(6)=1 00300 C and provide the lower (ML) and upper (MU) 00301 C bandwidths by setting 00302 C IWORK(1)=ML 00303 C IWORK(2)=MU **** 00304 C 00305 C 00306 C INFO(7) -- You can specify a maximum (absolute value of) 00307 C stepsize, so that the code 00308 C will avoid passing over very 00309 C large regions. 00310 C 00311 C **** Do you want the code to decide 00312 C on its own maximum stepsize? 00313 C Yes - Set INFO(7)=0 00314 C No - Set INFO(7)=1 00315 C and define HMAX by setting 00316 C RWORK(2)=HMAX **** 00317 C 00318 C INFO(8) -- Differential/algebraic problems 00319 C may occaisionally suffer from 00320 C severe scaling difficulties on the 00321 C first step. If you know a great deal 00322 C about the scaling of your problem, you can 00323 C help to alleviate this problem by 00324 C specifying an initial stepsize HO. 00325 C 00326 C **** Do you want the code to define 00327 C its own initial stepsize? 00328 C Yes - Set INFO(8)=0 00329 C No - Set INFO(8)=1 00330 C and define HO by setting 00331 C RWORK(3)=HO **** 00332 C 00333 C INFO(9) -- If storage is a severe problem, 00334 C you can save some locations by 00335 C restricting the maximum order MAXORD. 00336 C the default value is 5. for each 00337 C order decrease below 5, the code 00338 C requires NEQ fewer locations, however 00339 C it is likely to be slower. In any 00340 C case, you must have 1 .LE. MAXORD .LE. 5 00341 C **** Do you want the maximum order to 00342 C default to 5? 00343 C Yes - Set INFO(9)=0 00344 C No - Set INFO(9)=1 00345 C and define MAXORD by setting 00346 C IWORK(3)=MAXORD **** 00347 C 00348 C INFO(10) --If you know that the solutions to your equations 00349 C will always be nonnegative, it may help to set this 00350 C parameter. However, it is probably best to 00351 C try the code without using this option first, 00352 C and only to use this option if that doesn't 00353 C work very well. 00354 C **** Do you want the code to solve the problem without 00355 C invoking any special nonnegativity constraints? 00356 C Yes - Set INFO(10)=0 00357 C No - Set INFO(10)=1 00358 C 00359 C INFO(11) --DDASSL normally requires the initial T, 00360 C Y, and YPRIME to be consistent. That is, 00361 C you must have G(T,Y,YPRIME) = 0 at the initial 00362 C time. If you do not know the initial 00363 C derivative precisely, you can let DDASSL try 00364 C to compute it. 00365 C **** Are the initialHE INITIAL T, Y, YPRIME consistent? 00366 C Yes - Set INFO(11) = 0 00367 C No - Set INFO(11) = 1, 00368 C and set YPRIME to an initial approximation 00369 C to YPRIME. (If you have no idea what 00370 C YPRIME should be, set it to zero. Note 00371 C that the initial Y should be such 00372 C that there must exist a YPRIME so that 00373 C G(T,Y,YPRIME) = 0.) 00374 C 00375 C RTOL, ATOL -- You must assign relative (RTOL) and absolute (ATOL 00376 C error tolerances to tell the code how accurately you 00377 C want the solution to be computed. They must be defined 00378 C as variables because the code may change them. You 00379 C have two choices -- 00380 C Both RTOL and ATOL are scalars. (INFO(2)=0) 00381 C Both RTOL and ATOL are vectors. (INFO(2)=1) 00382 C in either case all components must be non-negative. 00383 C 00384 C The tolerances are used by the code in a local error 00385 C test at each step which requires roughly that 00386 C ABS(LOCAL ERROR) .LE. RTOL*ABS(Y)+ATOL 00387 C for each vector component. 00388 C (More specifically, a root-mean-square norm is used to 00389 C measure the size of vectors, and the error test uses the 00390 C magnitude of the solution at the beginning of the step.) 00391 C 00392 C The true (global) error is the difference between the 00393 C true solution of the initial value problem and the 00394 C computed approximation. Practically all present day 00395 C codes, including this one, control the local error at 00396 C each step and do not even attempt to control the global 00397 C error directly. 00398 C Usually, but not always, the true accuracy of the 00399 C computed Y is comparable to the error tolerances. This 00400 C code will usually, but not always, deliver a more 00401 C accurate solution if you reduce the tolerances and 00402 C integrate again. By comparing two such solutions you 00403 C can get a fairly reliable idea of the true error in the 00404 C solution at the bigger tolerances. 00405 C 00406 C Setting ATOL=0. results in a pure relative error test on 00407 C that component. Setting RTOL=0. results in a pure 00408 C absolute error test on that component. A mixed test 00409 C with non-zero RTOL and ATOL corresponds roughly to a 00410 C relative error test when the solution component is much 00411 C bigger than ATOL and to an absolute error test when the 00412 C solution component is smaller than the threshhold ATOL. 00413 C 00414 C The code will not attempt to compute a solution at an 00415 C accuracy unreasonable for the machine being used. It will 00416 C advise you if you ask for too much accuracy and inform 00417 C you as to the maximum accuracy it believes possible. 00418 C 00419 C RWORK(*) -- Dimension this real work array of length LRW in your 00420 C calling program. 00421 C 00422 C LRW -- Set it to the declared length of the RWORK array. 00423 C You must have 00424 C LRW .GE. 40+(MAXORD+4)*NEQ+NEQ**2 00425 C for the full (dense) JACOBIAN case (when INFO(6)=0), or 00426 C LRW .GE. 40+(MAXORD+4)*NEQ+(2*ML+MU+1)*NEQ 00427 C for the banded user-defined JACOBIAN case 00428 C (when INFO(5)=1 and INFO(6)=1), or 00429 C LRW .GE. 40+(MAXORD+4)*NEQ+(2*ML+MU+1)*NEQ 00430 C +2*(NEQ/(ML+MU+1)+1) 00431 C for the banded finite-difference-generated JACOBIAN case 00432 C (when INFO(5)=0 and INFO(6)=1) 00433 C 00434 C IWORK(*) -- Dimension this integer work array of length LIW in 00435 C your calling program. 00436 C 00437 C LIW -- Set it to the declared length of the IWORK array. 00438 C You must have LIW .GE. 21+NEQ 00439 C 00440 C RPAR, IPAR -- These are parameter arrays, of real and integer 00441 C type, respectively. You can use them for communication 00442 C between your program that calls DDASSL and the 00443 C RES subroutine (and the JAC subroutine). They are not 00444 C altered by DDASSL. If you do not need RPAR or IPAR, 00445 C ignore these parameters by treating them as dummy 00446 C arguments. If you do choose to use them, dimension 00447 C them in your calling program and in RES (and in JAC) 00448 C as arrays of appropriate length. 00449 C 00450 C JAC -- If you have set INFO(5)=0, you can ignore this parameter 00451 C by treating it as a dummy argument. Otherwise, you must 00452 C provide a subroutine of the form 00453 C SUBROUTINE JAC(T,Y,YPRIME,PD,CJ,RPAR,IPAR) 00454 C to define the matrix of partial derivatives 00455 C PD=DG/DY+CJ*DG/DYPRIME 00456 C CJ is a scalar which is input to JAC. 00457 C For the given values of T,Y,YPRIME, the 00458 C subroutine must evaluate the non-zero partial 00459 C derivatives for each equation and each solution 00460 C component, and store these values in the 00461 C matrix PD. The elements of PD are set to zero 00462 C before each call to JAC so only non-zero elements 00463 C need to be defined. 00464 C 00465 C Subroutine JAC must not alter T,Y,(*),YPRIME(*), or CJ. 00466 C You must declare the name JAC in an EXTERNAL statement in 00467 C your program that calls DDASSL. You must dimension Y, 00468 C YPRIME and PD in JAC. 00469 C 00470 C The way you must store the elements into the PD matrix 00471 C depends on the structure of the matrix which you 00472 C indicated by INFO(6). 00473 C *** INFO(6)=0 -- Full (dense) matrix *** 00474 C Give PD a first dimension of NEQ. 00475 C When you evaluate the (non-zero) partial derivative 00476 C of equation I with respect to variable J, you must 00477 C store it in PD according to 00478 C PD(I,J) = "DG(I)/DY(J)+CJ*DG(I)/DYPRIME(J)" 00479 C *** INFO(6)=1 -- Banded JACOBIAN with ML lower and MU 00480 C upper diagonal bands (refer to INFO(6) description 00481 C of ML and MU) *** 00482 C Give PD a first dimension of 2*ML+MU+1. 00483 C when you evaluate the (non-zero) partial derivative 00484 C of equation I with respect to variable J, you must 00485 C store it in PD according to 00486 C IROW = I - J + ML + MU + 1 00487 C PD(IROW,J) = "DG(I)/DY(J)+CJ*DG(I)/DYPRIME(J)" 00488 C 00489 C RPAR and IPAR are real and integer parameter arrays 00490 C which you can use for communication between your calling 00491 C program and your JACOBIAN subroutine JAC. They are not 00492 C altered by DDASSL. If you do not need RPAR or IPAR, 00493 C ignore these parameters by treating them as dummy 00494 C arguments. If you do choose to use them, dimension 00495 C them in your calling program and in JAC as arrays of 00496 C appropriate length. 00497 C 00498 C 00499 C OPTIONALLY REPLACEABLE NORM ROUTINE: 00500 C 00501 C DDASSL uses a weighted norm DDANRM to measure the size 00502 C of vectors such as the estimated error in each step. 00503 C A FUNCTION subprogram 00504 C DOUBLE PRECISION FUNCTION DDANRM(NEQ,V,WT,RPAR,IPAR) 00505 C DIMENSION V(NEQ),WT(NEQ) 00506 C is used to define this norm. Here, V is the vector 00507 C whose norm is to be computed, and WT is a vector of 00508 C weights. A DDANRM routine has been included with DDASSL 00509 C which computes the weighted root-mean-square norm 00510 C given by 00511 C DDANRM=SQRT((1/NEQ)*SUM(V(I)/WT(I))**2) 00512 C this norm is suitable for most problems. In some 00513 C special cases, it may be more convenient and/or 00514 C efficient to define your own norm by writing a function 00515 C subprogram to be called instead of DDANRM. This should, 00516 C however, be attempted only after careful thought and 00517 C consideration. 00518 C 00519 C 00520 C -------- OUTPUT -- AFTER ANY RETURN FROM DDASSL --------------------- 00521 C 00522 C The principal aim of the code is to return a computed solution at 00523 C TOUT, although it is also possible to obtain intermediate results 00524 C along the way. To find out whether the code achieved its goal 00525 C or if the integration process was interrupted before the task was 00526 C completed, you must check the IDID parameter. 00527 C 00528 C 00529 C T -- The solution was successfully advanced to the 00530 C output value of T. 00531 C 00532 C Y(*) -- Contains the computed solution approximation at T. 00533 C 00534 C YPRIME(*) -- Contains the computed derivative 00535 C approximation at T. 00536 C 00537 C IDID -- Reports what the code did. 00538 C 00539 C *** Task completed *** 00540 C Reported by positive values of IDID 00541 C 00542 C IDID = 1 -- A step was successfully taken in the 00543 C intermediate-output mode. The code has not 00544 C yet reached TOUT. 00545 C 00546 C IDID = 2 -- The integration to TSTOP was successfully 00547 C completed (T=TSTOP) by stepping exactly to TSTOP. 00548 C 00549 C IDID = 3 -- The integration to TOUT was successfully 00550 C completed (T=TOUT) by stepping past TOUT. 00551 C Y(*) is obtained by interpolation. 00552 C YPRIME(*) is obtained by interpolation. 00553 C 00554 C *** Task interrupted *** 00555 C Reported by negative values of IDID 00556 C 00557 C IDID = -1 -- A large amount of work has been expended. 00558 C (About 500 steps) 00559 C 00560 C IDID = -2 -- The error tolerances are too stringent. 00561 C 00562 C IDID = -3 -- The local error test cannot be satisfied 00563 C because you specified a zero component in ATOL 00564 C and the corresponding computed solution 00565 C component is zero. Thus, a pure relative error 00566 C test is impossible for this component. 00567 C 00568 C IDID = -6 -- DDASSL had repeated error test 00569 C failures on the last attempted step. 00570 C 00571 C IDID = -7 -- The corrector could not converge. 00572 C 00573 C IDID = -8 -- The matrix of partial derivatives 00574 C is singular. 00575 C 00576 C IDID = -9 -- The corrector could not converge. 00577 C there were repeated error test failures 00578 C in this step. 00579 C 00580 C IDID =-10 -- The corrector could not converge 00581 C because IRES was equal to minus one. 00582 C 00583 C IDID =-11 -- IRES equal to -2 was encountered 00584 C and control is being returned to the 00585 C calling program. 00586 C 00587 C IDID =-12 -- DDASSL failed to compute the initial 00588 C YPRIME. 00589 C 00590 C 00591 C 00592 C IDID = -13,..,-32 -- Not applicable for this code 00593 C 00594 C *** Task terminated *** 00595 C Reported by the value of IDID=-33 00596 C 00597 C IDID = -33 -- The code has encountered trouble from which 00598 C it cannot recover. A message is printed 00599 C explaining the trouble and control is returned 00600 C to the calling program. For example, this occurs 00601 C when invalid input is detected. 00602 C 00603 C RTOL, ATOL -- These quantities remain unchanged except when 00604 C IDID = -2. In this case, the error tolerances have been 00605 C increased by the code to values which are estimated to 00606 C be appropriate for continuing the integration. However, 00607 C the reported solution at T was obtained using the input 00608 C values of RTOL and ATOL. 00609 C 00610 C RWORK, IWORK -- Contain information which is usually of no 00611 C interest to the user but necessary for subsequent calls. 00612 C However, you may find use for 00613 C 00614 C RWORK(3)--Which contains the step size H to be 00615 C attempted on the next step. 00616 C 00617 C RWORK(4)--Which contains the current value of the 00618 C independent variable, i.e., the farthest point 00619 C integration has reached. This will be different 00620 C from T only when interpolation has been 00621 C performed (IDID=3). 00622 C 00623 C RWORK(7)--Which contains the stepsize used 00624 C on the last successful step. 00625 C 00626 C IWORK(7)--Which contains the order of the method to 00627 C be attempted on the next step. 00628 C 00629 C IWORK(8)--Which contains the order of the method used 00630 C on the last step. 00631 C 00632 C IWORK(11)--Which contains the number of steps taken so 00633 C far. 00634 C 00635 C IWORK(12)--Which contains the number of calls to RES 00636 C so far. 00637 C 00638 C IWORK(13)--Which contains the number of evaluations of 00639 C the matrix of partial derivatives needed so 00640 C far. 00641 C 00642 C IWORK(14)--Which contains the total number 00643 C of error test failures so far. 00644 C 00645 C IWORK(15)--Which contains the total number 00646 C of convergence test failures so far. 00647 C (includes singular iteration matrix 00648 C failures.) 00649 C 00650 C 00651 C -------- INPUT -- WHAT TO DO TO CONTINUE THE INTEGRATION ------------ 00652 C (CALLS AFTER THE FIRST) 00653 C 00654 C This code is organized so that subsequent calls to continue the 00655 C integration involve little (if any) additional effort on your 00656 C part. You must monitor the IDID parameter in order to determine 00657 C what to do next. 00658 C 00659 C Recalling that the principal task of the code is to integrate 00660 C from T to TOUT (the interval mode), usually all you will need 00661 C to do is specify a new TOUT upon reaching the current TOUT. 00662 C 00663 C Do not alter any quantity not specifically permitted below, 00664 C in particular do not alter NEQ,T,Y(*),YPRIME(*),RWORK(*),IWORK(*) 00665 C or the differential equation in subroutine RES. Any such 00666 C alteration constitutes a new problem and must be treated as such, 00667 C i.e., you must start afresh. 00668 C 00669 C You cannot change from vector to scalar error control or vice 00670 C versa (INFO(2)), but you can change the size of the entries of 00671 C RTOL, ATOL. Increasing a tolerance makes the equation easier 00672 C to integrate. Decreasing a tolerance will make the equation 00673 C harder to integrate and should generally be avoided. 00674 C 00675 C You can switch from the intermediate-output mode to the 00676 C interval mode (INFO(3)) or vice versa at any time. 00677 C 00678 C If it has been necessary to prevent the integration from going 00679 C past a point TSTOP (INFO(4), RWORK(1)), keep in mind that the 00680 C code will not integrate to any TOUT beyond the currently 00681 C specified TSTOP. Once TSTOP has been reached you must change 00682 C the value of TSTOP or set INFO(4)=0. You may change INFO(4) 00683 C or TSTOP at any time but you must supply the value of TSTOP in 00684 C RWORK(1) whenever you set INFO(4)=1. 00685 C 00686 C Do not change INFO(5), INFO(6), IWORK(1), or IWORK(2) 00687 C unless you are going to restart the code. 00688 C 00689 C *** Following a completed task *** 00690 C If 00691 C IDID = 1, call the code again to continue the integration 00692 C another step in the direction of TOUT. 00693 C 00694 C IDID = 2 or 3, define a new TOUT and call the code again. 00695 C TOUT must be different from T. You cannot change 00696 C the direction of integration without restarting. 00697 C 00698 C *** Following an interrupted task *** 00699 C To show the code that you realize the task was 00700 C interrupted and that you want to continue, you 00701 C must take appropriate action and set INFO(1) = 1 00702 C If 00703 C IDID = -1, The code has taken about 500 steps. 00704 C If you want to continue, set INFO(1) = 1 and 00705 C call the code again. An additional 500 steps 00706 C will be allowed. 00707 C 00708 C IDID = -2, The error tolerances RTOL, ATOL have been 00709 C increased to values the code estimates appropriate 00710 C for continuing. You may want to change them 00711 C yourself. If you are sure you want to continue 00712 C with relaxed error tolerances, set INFO(1)=1 and 00713 C call the code again. 00714 C 00715 C IDID = -3, A solution component is zero and you set the 00716 C corresponding component of ATOL to zero. If you 00717 C are sure you want to continue, you must first 00718 C alter the error criterion to use positive values 00719 C for those components of ATOL corresponding to zero 00720 C solution components, then set INFO(1)=1 and call 00721 C the code again. 00722 C 00723 C IDID = -4,-5 --- Cannot occur with this code. 00724 C 00725 C IDID = -6, Repeated error test failures occurred on the 00726 C last attempted step in DDASSL. A singularity in the 00727 C solution may be present. If you are absolutely 00728 C certain you want to continue, you should restart 00729 C the integration. (Provide initial values of Y and 00730 C YPRIME which are consistent) 00731 C 00732 C IDID = -7, Repeated convergence test failures occurred 00733 C on the last attempted step in DDASSL. An inaccurate 00734 C or ill-conditioned JACOBIAN may be the problem. If 00735 C you are absolutely certain you want to continue, you 00736 C should restart the integration. 00737 C 00738 C IDID = -8, The matrix of partial derivatives is singular. 00739 C Some of your equations may be redundant. 00740 C DDASSL cannot solve the problem as stated. 00741 C It is possible that the redundant equations 00742 C could be removed, and then DDASSL could 00743 C solve the problem. It is also possible 00744 C that a solution to your problem either 00745 C does not exist or is not unique. 00746 C 00747 C IDID = -9, DDASSL had multiple convergence test 00748 C failures, preceeded by multiple error 00749 C test failures, on the last attempted step. 00750 C It is possible that your problem 00751 C is ill-posed, and cannot be solved 00752 C using this code. Or, there may be a 00753 C discontinuity or a singularity in the 00754 C solution. If you are absolutely certain 00755 C you want to continue, you should restart 00756 C the integration. 00757 C 00758 C IDID =-10, DDASSL had multiple convergence test failures 00759 C because IRES was equal to minus one. 00760 C If you are absolutely certain you want 00761 C to continue, you should restart the 00762 C integration. 00763 C 00764 C IDID =-11, IRES=-2 was encountered, and control is being 00765 C returned to the calling program. 00766 C 00767 C IDID =-12, DDASSL failed to compute the initial YPRIME. 00768 C This could happen because the initial 00769 C approximation to YPRIME was not very good, or 00770 C if a YPRIME consistent with the initial Y 00771 C does not exist. The problem could also be caused 00772 C by an inaccurate or singular iteration matrix. 00773 C 00774 C IDID = -13,..,-32 --- Cannot occur with this code. 00775 C 00776 C 00777 C *** Following a terminated task *** 00778 C 00779 C If IDID= -33, you cannot continue the solution of this problem. 00780 C An attempt to do so will result in your 00781 C run being terminated. 00782 C 00783 C 00784 C -------- ERROR MESSAGES --------------------------------------------- 00785 C 00786 C The SLATEC error print routine XERMSG is called in the event of 00787 C unsuccessful completion of a task. Most of these are treated as 00788 C "recoverable errors", which means that (unless the user has directed 00789 C otherwise) control will be returned to the calling program for 00790 C possible action after the message has been printed. 00791 C 00792 C In the event of a negative value of IDID other than -33, an appro- 00793 C priate message is printed and the "error number" printed by XERMSG 00794 C is the value of IDID. There are quite a number of illegal input 00795 C errors that can lead to a returned value IDID=-33. The conditions 00796 C and their printed "error numbers" are as follows: 00797 C 00798 C Error number Condition 00799 C 00800 C 1 Some element of INFO vector is not zero or one. 00801 C 2 NEQ .le. 0 00802 C 3 MAXORD not in range. 00803 C 4 LRW is less than the required length for RWORK. 00804 C 5 LIW is less than the required length for IWORK. 00805 C 6 Some element of RTOL is .lt. 0 00806 C 7 Some element of ATOL is .lt. 0 00807 C 8 All elements of RTOL and ATOL are zero. 00808 C 9 INFO(4)=1 and TSTOP is behind TOUT. 00809 C 10 HMAX .lt. 0.0 00810 C 11 TOUT is behind T. 00811 C 12 INFO(8)=1 and H0=0.0 00812 C 13 Some element of WT is .le. 0.0 00813 C 14 TOUT is too close to T to start integration. 00814 C 15 INFO(4)=1 and TSTOP is behind T. 00815 C 16 --( Not used in this version )-- 00816 C 17 ML illegal. Either .lt. 0 or .gt. NEQ 00817 C 18 MU illegal. Either .lt. 0 or .gt. NEQ 00818 C 19 TOUT = T. 00819 C 00820 C If DDASSL is called again without any action taken to remove the 00821 C cause of an unsuccessful return, XERMSG will be called with a fatal 00822 C error flag, which will cause unconditional termination of the 00823 C program. There are two such fatal errors: 00824 C 00825 C Error number -998: The last step was terminated with a negative 00826 C value of IDID other than -33, and no appropriate action was 00827 C taken. 00828 C 00829 C Error number -999: The previous call was terminated because of 00830 C illegal input (IDID=-33) and there is illegal input in the 00831 C present call, as well. (Suspect infinite loop.) 00832 C 00833 C --------------------------------------------------------------------- 00834 C 00835 C***REFERENCES A DESCRIPTION OF DASSL: A DIFFERENTIAL/ALGEBRAIC 00836 C SYSTEM SOLVER, L. R. PETZOLD, SAND82-8637, 00837 C SANDIA NATIONAL LABORATORIES, SEPTEMBER 1982. 00838 C***ROUTINES CALLED D1MACH, DDAINI, DDANRM, DDASTP, DDATRP, DDAWTS, 00839 C XERMSG 00840 C***REVISION HISTORY (YYMMDD) 00841 C 830315 DATE WRITTEN 00842 C 880387 Code changes made. All common statements have been 00843 C replaced by a DATA statement, which defines pointers into 00844 C RWORK, and PARAMETER statements which define pointers 00845 C into IWORK. As well the documentation has gone through 00846 C grammatical changes. 00847 C 881005 The prologue has been changed to mixed case. 00848 C The subordinate routines had revision dates changed to 00849 C this date, although the documentation for these routines 00850 C is all upper case. No code changes. 00851 C 890511 Code changes made. The DATA statement in the declaration 00852 C section of DDASSL was replaced with a PARAMETER 00853 C statement. Also the statement S = 100.D0 was removed 00854 C from the top of the Newton iteration in DDASTP. 00855 C The subordinate routines had revision dates changed to 00856 C this date. 00857 C 890517 The revision date syntax was replaced with the revision 00858 C history syntax. Also the "DECK" comment was added to 00859 C the top of all subroutines. These changes are consistent 00860 C with new SLATEC guidelines. 00861 C The subordinate routines had revision dates changed to 00862 C this date. No code changes. 00863 C 891013 Code changes made. 00864 C Removed all occurrances of FLOAT or DBLE. All operations 00865 C are now performed with "mixed-mode" arithmetic. 00866 C Also, specific function names were replaced with generic 00867 C function names to be consistent with new SLATEC guidelines. 00868 C In particular: 00869 C Replaced DSQRT with SQRT everywhere. 00870 C Replaced DABS with ABS everywhere. 00871 C Replaced DMIN1 with MIN everywhere. 00872 C Replaced MIN0 with MIN everywhere. 00873 C Replaced DMAX1 with MAX everywhere. 00874 C Replaced MAX0 with MAX everywhere. 00875 C Replaced DSIGN with SIGN everywhere. 00876 C Also replaced REVISION DATE with REVISION HISTORY in all 00877 C subordinate routines. 00878 C 901004 Miscellaneous changes to prologue to complete conversion 00879 C to SLATEC 4.0 format. No code changes. (F.N.Fritsch) 00880 C 901009 Corrected GAMS classification code and converted subsidiary 00881 C routines to 4.0 format. No code changes. (F.N.Fritsch) 00882 C 901010 Converted XERRWV calls to XERMSG calls. (R.Clemens,AFWL) 00883 C 901019 Code changes made. 00884 C Merged SLATEC 4.0 changes with previous changes made 00885 C by C. Ulrich. Below is a history of the changes made by 00886 C C. Ulrich. (Changes in subsidiary routines are implied 00887 C by this history) 00888 C 891228 Bug was found and repaired inside the DDASSL 00889 C and DDAINI routines. DDAINI was incorrectly 00890 C returning the initial T with Y and YPRIME 00891 C computed at T+H. The routine now returns T+H 00892 C rather than the initial T. 00893 C Cosmetic changes made to DDASTP. 00894 C 900904 Three modifications were made to fix a bug (inside 00895 C DDASSL) re interpolation for continuation calls and 00896 C cases where TN is very close to TSTOP: 00897 C 00898 C 1) In testing for whether H is too large, just 00899 C compare H to (TSTOP - TN), rather than 00900 C (TSTOP - TN) * (1-4*UROUND), and set H to 00901 C TSTOP - TN. This will force DDASTP to step 00902 C exactly to TSTOP under certain situations 00903 C (i.e. when H returned from DDASTP would otherwise 00904 C take TN beyond TSTOP). 00905 C 00906 C 2) Inside the DDASTP loop, interpolate exactly to 00907 C TSTOP if TN is very close to TSTOP (rather than 00908 C interpolating to within roundoff of TSTOP). 00909 C 00910 C 3) Modified IDID description for IDID = 2 to say that 00911 C the solution is returned by stepping exactly to 00912 C TSTOP, rather than TOUT. (In some cases the 00913 C solution is actually obtained by extrapolating 00914 C over a distance near unit roundoff to TSTOP, 00915 C but this small distance is deemed acceptable in 00916 C these circumstances.) 00917 C 901026 Added explicit declarations for all variables and minor 00918 C cosmetic changes to prologue, removed unreferenced labels, 00919 C and improved XERMSG calls. (FNF) 00920 C 901030 Added ERROR MESSAGES section and reworked other sections to 00921 C be of more uniform format. (FNF) 00922 C 910624 Fixed minor bug related to HMAX (five lines ending in 00923 C statement 526 in DDASSL). (LRP) 00924 C 00925 C***END PROLOGUE DDASSL 00926 C 00927 C**End 00928 C 00929 C Declare arguments. 00930 C 00931 INTEGER NEQ, INFO(15), IDID, LRW, IWORK(*), LIW, IPAR(*) 00932 DOUBLE PRECISION 00933 * T, Y(*), YPRIME(*), TOUT, RTOL(*), ATOL(*), RWORK(*), 00934 * RPAR(*) 00935 EXTERNAL RES, JAC 00936 C 00937 C Declare externals. 00938 C 00939 EXTERNAL D1MACH, DDAINI, DDANRM, DDASTP, DDATRP, DDAWTS, XERMSG 00940 DOUBLE PRECISION D1MACH, DDANRM 00941 C 00942 C Declare local variables. 00943 C 00944 INTEGER I, ITEMP, LALPHA, LBETA, LCJ, LCJOLD, LCTF, LDELTA, 00945 * LENIW, LENPD, LENRW, LE, LETF, LGAMMA, LH, LHMAX, LHOLD, 00946 * LMXSTP, LIPVT, 00947 * LJCALC, LK, LKOLD, LIWM, LML, LMTYPE, LMU, LMXORD, LNJE, LNPD, 00948 * LNRE, LNS, LNST, LNSTL, LPD, LPHASE, LPHI, LPSI, LROUND, LS, 00949 * LSIGMA, LTN, LTSTOP, LWM, LWT, MBAND, MSAVE, MXORD, NPD, NTEMP, 00950 * NZFLG 00951 DOUBLE PRECISION 00952 * ATOLI, H, HMAX, HMIN, HO, R, RH, RTOLI, TDIST, TN, TNEXT, 00953 * TSTOP, UROUND, YPNORM 00954 LOGICAL DONE 00955 C Auxiliary variables for conversion of values to be included in 00956 C error messages. 00957 CHARACTER*8 XERN1, XERN2 00958 CHARACTER*16 XERN3, XERN4 00959 C 00960 C SET POINTERS INTO IWORK 00961 PARAMETER (LML=1, LMU=2, LMXORD=3, LMTYPE=4, LNST=11, 00962 * LNRE=12, LNJE=13, LETF=14, LCTF=15, LNPD=16, LMXSTP=21, 00963 * LIPVT=22, LJCALC=5, LPHASE=6, LK=7, LKOLD=8, 00964 * LNS=9, LNSTL=10, LIWM=1) 00965 C 00966 C SET RELATIVE OFFSET INTO RWORK 00967 PARAMETER (NPD=1) 00968 C 00969 C SET POINTERS INTO RWORK 00970 PARAMETER (LTSTOP=1, LHMAX=2, LH=3, LTN=4, 00971 * LCJ=5, LCJOLD=6, LHOLD=7, LS=8, LROUND=9, 00972 * LALPHA=11, LBETA=17, LGAMMA=23, 00973 * LPSI=29, LSIGMA=35, LDELTA=41) 00974 C 00975 C***FIRST EXECUTABLE STATEMENT DDASSL 00976 IF(INFO(1).NE.0)GO TO 100 00977 C 00978 C----------------------------------------------------------------------- 00979 C THIS BLOCK IS EXECUTED FOR THE INITIAL CALL ONLY. 00980 C IT CONTAINS CHECKING OF INPUTS AND INITIALIZATIONS. 00981 C----------------------------------------------------------------------- 00982 C 00983 C FIRST CHECK INFO ARRAY TO MAKE SURE ALL ELEMENTS OF INFO 00984 C ARE EITHER ZERO OR ONE. 00985 DO 10 I=2,11 00986 IF(INFO(I).NE.0.AND.INFO(I).NE.1)GO TO 701 00987 10 CONTINUE 00988 C 00989 IF(NEQ.LE.0)GO TO 702 00990 C 00991 C CHECK AND COMPUTE MAXIMUM ORDER 00992 MXORD=5 00993 IF(INFO(9).EQ.0)GO TO 20 00994 MXORD=IWORK(LMXORD) 00995 IF(MXORD.LT.1.OR.MXORD.GT.5)GO TO 703 00996 20 IWORK(LMXORD)=MXORD 00997 C 00998 C COMPUTE MTYPE,LENPD,LENRW.CHECK ML AND MU. 00999 IF(INFO(6).NE.0)GO TO 40 01000 LENPD=NEQ**2 01001 LENRW=40+(IWORK(LMXORD)+4)*NEQ+LENPD 01002 IF(INFO(5).NE.0)GO TO 30 01003 IWORK(LMTYPE)=2 01004 GO TO 60 01005 30 IWORK(LMTYPE)=1 01006 GO TO 60 01007 40 IF(IWORK(LML).LT.0.OR.IWORK(LML).GE.NEQ)GO TO 717 01008 IF(IWORK(LMU).LT.0.OR.IWORK(LMU).GE.NEQ)GO TO 718 01009 LENPD=(2*IWORK(LML)+IWORK(LMU)+1)*NEQ 01010 IF(INFO(5).NE.0)GO TO 50 01011 IWORK(LMTYPE)=5 01012 MBAND=IWORK(LML)+IWORK(LMU)+1 01013 MSAVE=(NEQ/MBAND)+1 01014 LENRW=40+(IWORK(LMXORD)+4)*NEQ+LENPD+2*MSAVE 01015 GO TO 60 01016 50 IWORK(LMTYPE)=4 01017 LENRW=40+(IWORK(LMXORD)+4)*NEQ+LENPD 01018 C 01019 C CHECK LENGTHS OF RWORK AND IWORK 01020 60 LENIW=21+NEQ 01021 IWORK(LNPD)=LENPD 01022 IF(LRW.LT.LENRW)GO TO 704 01023 IF(LIW.LT.LENIW)GO TO 705 01024 C 01025 C CHECK TO SEE THAT TOUT IS DIFFERENT FROM T 01026 IF(TOUT .EQ. T)GO TO 719 01027 C 01028 C CHECK HMAX 01029 IF(INFO(7).EQ.0)GO TO 70 01030 HMAX=RWORK(LHMAX) 01031 IF(HMAX.LE.0.0D0)GO TO 710 01032 70 CONTINUE 01033 C 01034 C CHECK AND COMPUTE MAXIMUM STEPS 01035 MXSTP=500 01036 IF(INFO(12).EQ.0)GO TO 80 01037 MXSTP=IWORK(LMXSTP) 01038 IF(MXSTP.LT.0)GO TO 716 01039 80 IWORK(LMXSTP)=MXSTP 01040 C 01041 C INITIALIZE COUNTERS 01042 IWORK(LNST)=0 01043 IWORK(LNRE)=0 01044 IWORK(LNJE)=0 01045 C 01046 IWORK(LNSTL)=0 01047 IDID=1 01048 GO TO 200 01049 C 01050 C----------------------------------------------------------------------- 01051 C THIS BLOCK IS FOR CONTINUATION CALLS 01052 C ONLY. HERE WE CHECK INFO(1),AND IF THE 01053 C LAST STEP WAS INTERRUPTED WE CHECK WHETHER 01054 C APPROPRIATE ACTION WAS TAKEN. 01055 C----------------------------------------------------------------------- 01056 C 01057 100 CONTINUE 01058 IF(INFO(1).EQ.1)GO TO 110 01059 IF(INFO(1).NE.-1)GO TO 701 01060 C 01061 C IF WE ARE HERE, THE LAST STEP WAS INTERRUPTED 01062 C BY AN ERROR CONDITION FROM DDASTP,AND 01063 C APPROPRIATE ACTION WAS NOT TAKEN. THIS 01064 C IS A FATAL ERROR. 01065 WRITE (XERN1, '(I8)') IDID 01066 CALL XERMSG ('SLATEC', 'DDASSL', 01067 * 'THE LAST STEP TERMINATED WITH A NEGATIVE VALUE OF IDID = ' // 01068 * XERN1 // ' AND NO APPROPRIATE ACTION WAS TAKEN. ' // 01069 * 'RUN TERMINATED', -998, 2) 01070 RETURN 01071 110 CONTINUE 01072 IWORK(LNSTL)=IWORK(LNST) 01073 C 01074 C----------------------------------------------------------------------- 01075 C THIS BLOCK IS EXECUTED ON ALL CALLS. 01076 C THE ERROR TOLERANCE PARAMETERS ARE 01077 C CHECKED, AND THE WORK ARRAY POINTERS 01078 C ARE SET. 01079 C----------------------------------------------------------------------- 01080 C 01081 200 CONTINUE 01082 C CHECK RTOL,ATOL 01083 NZFLG=0 01084 RTOLI=RTOL(1) 01085 ATOLI=ATOL(1) 01086 DO 210 I=1,NEQ 01087 IF(INFO(2).EQ.1)RTOLI=RTOL(I) 01088 IF(INFO(2).EQ.1)ATOLI=ATOL(I) 01089 IF(RTOLI.GT.0.0D0.OR.ATOLI.GT.0.0D0)NZFLG=1 01090 IF(RTOLI.LT.0.0D0)GO TO 706 01091 IF(ATOLI.LT.0.0D0)GO TO 707 01092 210 CONTINUE 01093 IF(NZFLG.EQ.0)GO TO 708 01094 C 01095 C SET UP RWORK STORAGE.IWORK STORAGE IS FIXED 01096 C IN DATA STATEMENT. 01097 LE=LDELTA+NEQ 01098 LWT=LE+NEQ 01099 LPHI=LWT+NEQ 01100 LPD=LPHI+(IWORK(LMXORD)+1)*NEQ 01101 LWM=LPD 01102 NTEMP=NPD+IWORK(LNPD) 01103 IF(INFO(1).EQ.1)GO TO 400 01104 C 01105 C----------------------------------------------------------------------- 01106 C THIS BLOCK IS EXECUTED ON THE INITIAL CALL 01107 C ONLY. SET THE INITIAL STEP SIZE, AND 01108 C THE ERROR WEIGHT VECTOR, AND PHI. 01109 C COMPUTE INITIAL YPRIME, IF NECESSARY. 01110 C----------------------------------------------------------------------- 01111 C 01112 TN=T 01113 IDID=1 01114 C 01115 C SET ERROR WEIGHT VECTOR WT 01116 CALL DDAWTS(NEQ,INFO(2),RTOL,ATOL,Y,RWORK(LWT),RPAR,IPAR) 01117 DO 305 I = 1,NEQ 01118 IF(RWORK(LWT+I-1).LE.0.0D0) GO TO 713 01119 305 CONTINUE 01120 C 01121 C COMPUTE UNIT ROUNDOFF AND HMIN 01122 UROUND = D1MACH(4) 01123 RWORK(LROUND) = UROUND 01124 HMIN = 4.0D0*UROUND*MAX(ABS(T),ABS(TOUT)) 01125 C 01126 C CHECK INITIAL INTERVAL TO SEE THAT IT IS LONG ENOUGH 01127 TDIST = ABS(TOUT - T) 01128 IF(TDIST .LT. HMIN) GO TO 714 01129 C 01130 C CHECK HO, IF THIS WAS INPUT 01131 IF (INFO(8) .EQ. 0) GO TO 310 01132 HO = RWORK(LH) 01133 IF ((TOUT - T)*HO .LT. 0.0D0) GO TO 711 01134 IF (HO .EQ. 0.0D0) GO TO 712 01135 GO TO 320 01136 310 CONTINUE 01137 C 01138 C COMPUTE INITIAL STEPSIZE, TO BE USED BY EITHER 01139 C DDASTP OR DDAINI, DEPENDING ON INFO(11) 01140 HO = 0.001D0*TDIST 01141 YPNORM = DDANRM(NEQ,YPRIME,RWORK(LWT),RPAR,IPAR) 01142 IF (YPNORM .GT. 0.5D0/HO) HO = 0.5D0/YPNORM 01143 HO = SIGN(HO,TOUT-T) 01144 C ADJUST HO IF NECESSARY TO MEET HMAX BOUND 01145 320 IF (INFO(7) .EQ. 0) GO TO 330 01146 RH = ABS(HO)/RWORK(LHMAX) 01147 IF (RH .GT. 1.0D0) HO = HO/RH 01148 C COMPUTE TSTOP, IF APPLICABLE 01149 330 IF (INFO(4) .EQ. 0) GO TO 340 01150 TSTOP = RWORK(LTSTOP) 01151 IF ((TSTOP - T)*HO .LT. 0.0D0) GO TO 715 01152 IF ((T + HO - TSTOP)*HO .GT. 0.0D0) HO = TSTOP - T 01153 IF ((TSTOP - TOUT)*HO .LT. 0.0D0) GO TO 709 01154 C 01155 C COMPUTE INITIAL DERIVATIVE, UPDATING TN AND Y, IF APPLICABLE 01156 340 IF (INFO(11) .EQ. 0) GO TO 350 01157 CALL DDAINI(TN,Y,YPRIME,NEQ, 01158 * RES,JAC,HO,RWORK(LWT),IDID,RPAR,IPAR, 01159 * RWORK(LPHI),RWORK(LDELTA),RWORK(LE), 01160 * RWORK(LWM),IWORK(LIWM),HMIN,RWORK(LROUND), 01161 * INFO(10),NTEMP) 01162 IF (IDID .LT. 0) GO TO 390 01163 C 01164 C LOAD H WITH HO. STORE H IN RWORK(LH) 01165 350 H = HO 01166 RWORK(LH) = H 01167 C 01168 C LOAD Y AND H*YPRIME INTO PHI(*,1) AND PHI(*,2) 01169 ITEMP = LPHI + NEQ 01170 DO 370 I = 1,NEQ 01171 RWORK(LPHI + I - 1) = Y(I) 01172 370 RWORK(ITEMP + I - 1) = H*YPRIME(I) 01173 C 01174 390 GO TO 500 01175 C 01176 C------------------------------------------------------- 01177 C THIS BLOCK IS FOR CONTINUATION CALLS ONLY. ITS 01178 C PURPOSE IS TO CHECK STOP CONDITIONS BEFORE 01179 C TAKING A STEP. 01180 C ADJUST H IF NECESSARY TO MEET HMAX BOUND 01181 C------------------------------------------------------- 01182 C 01183 400 CONTINUE 01184 UROUND=RWORK(LROUND) 01185 DONE = .FALSE. 01186 TN=RWORK(LTN) 01187 H=RWORK(LH) 01188 IF(INFO(7) .EQ. 0) GO TO 410 01189 RH = ABS(H)/RWORK(LHMAX) 01190 IF(RH .GT. 1.0D0) H = H/RH 01191 410 CONTINUE 01192 IF(T .EQ. TOUT) GO TO 719 01193 IF((T - TOUT)*H .GT. 0.0D0) GO TO 711 01194 IF(INFO(4) .EQ. 1) GO TO 430 01195 IF(INFO(3) .EQ. 1) GO TO 420 01196 IF((TN-TOUT)*H.LT.0.0D0)GO TO 490 01197 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD), 01198 * RWORK(LPHI),RWORK(LPSI)) 01199 T=TOUT 01200 IDID = 3 01201 DONE = .TRUE. 01202 GO TO 490 01203 420 IF((TN-T)*H .LE. 0.0D0) GO TO 490 01204 IF((TN - TOUT)*H .GT. 0.0D0) GO TO 425 01205 CALL DDATRP(TN,TN,Y,YPRIME,NEQ,IWORK(LKOLD), 01206 * RWORK(LPHI),RWORK(LPSI)) 01207 T = TN 01208 IDID = 1 01209 DONE = .TRUE. 01210 GO TO 490 01211 425 CONTINUE 01212 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD), 01213 * RWORK(LPHI),RWORK(LPSI)) 01214 T = TOUT 01215 IDID = 3 01216 DONE = .TRUE. 01217 GO TO 490 01218 430 IF(INFO(3) .EQ. 1) GO TO 440 01219 TSTOP=RWORK(LTSTOP) 01220 IF((TN-TSTOP)*H.GT.0.0D0) GO TO 715 01221 IF((TSTOP-TOUT)*H.LT.0.0D0)GO TO 709 01222 IF((TN-TOUT)*H.LT.0.0D0)GO TO 450 01223 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD), 01224 * RWORK(LPHI),RWORK(LPSI)) 01225 T=TOUT 01226 IDID = 3 01227 DONE = .TRUE. 01228 GO TO 490 01229 440 TSTOP = RWORK(LTSTOP) 01230 IF((TN-TSTOP)*H .GT. 0.0D0) GO TO 715 01231 IF((TSTOP-TOUT)*H .LT. 0.0D0) GO TO 709 01232 IF((TN-T)*H .LE. 0.0D0) GO TO 450 01233 IF((TN - TOUT)*H .GT. 0.0D0) GO TO 445 01234 CALL DDATRP(TN,TN,Y,YPRIME,NEQ,IWORK(LKOLD), 01235 * RWORK(LPHI),RWORK(LPSI)) 01236 T = TN 01237 IDID = 1 01238 DONE = .TRUE. 01239 GO TO 490 01240 445 CONTINUE 01241 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD), 01242 * RWORK(LPHI),RWORK(LPSI)) 01243 T = TOUT 01244 IDID = 3 01245 DONE = .TRUE. 01246 GO TO 490 01247 450 CONTINUE 01248 C CHECK WHETHER WE ARE WITHIN ROUNDOFF OF TSTOP 01249 IF(ABS(TN-TSTOP).GT.100.0D0*UROUND* 01250 * (ABS(TN)+ABS(H)))GO TO 460 01251 CALL DDATRP(TN,TSTOP,Y,YPRIME,NEQ,IWORK(LKOLD), 01252 * RWORK(LPHI),RWORK(LPSI)) 01253 IDID=2 01254 T=TSTOP 01255 DONE = .TRUE. 01256 GO TO 490 01257 460 TNEXT=TN+H 01258 IF((TNEXT-TSTOP)*H.LE.0.0D0)GO TO 490 01259 H=TSTOP-TN 01260 RWORK(LH)=H 01261 C 01262 490 IF (DONE) GO TO 580 01263 C 01264 C------------------------------------------------------- 01265 C THE NEXT BLOCK CONTAINS THE CALL TO THE 01266 C ONE-STEP INTEGRATOR DDASTP. 01267 C THIS IS A LOOPING POINT FOR THE INTEGRATION STEPS. 01268 C CHECK FOR TOO MANY STEPS. 01269 C UPDATE WT. 01270 C CHECK FOR TOO MUCH ACCURACY REQUESTED. 01271 C COMPUTE MINIMUM STEPSIZE. 01272 C------------------------------------------------------- 01273 C 01274 500 CONTINUE 01275 C CHECK FOR FAILURE TO COMPUTE INITIAL YPRIME 01276 IF (IDID .EQ. -12) GO TO 527 01277 C 01278 C CHECK FOR TOO MANY STEPS 01279 IF((IWORK(LNST)-IWORK(LNSTL)).LT.IWORK(LMXSTP)) 01280 * GO TO 510 01281 IDID=-1 01282 GO TO 527 01283 C 01284 C UPDATE WT 01285 510 CALL DDAWTS(NEQ,INFO(2),RTOL,ATOL,RWORK(LPHI), 01286 * RWORK(LWT),RPAR,IPAR) 01287 DO 520 I=1,NEQ 01288 IF(RWORK(I+LWT-1).GT.0.0D0)GO TO 520 01289 IDID=-3 01290 GO TO 527 01291 520 CONTINUE 01292 C 01293 C TEST FOR TOO MUCH ACCURACY REQUESTED. 01294 R=DDANRM(NEQ,RWORK(LPHI),RWORK(LWT),RPAR,IPAR)* 01295 * 100.0D0*UROUND 01296 IF(R.LE.1.0D0)GO TO 525 01297 C MULTIPLY RTOL AND ATOL BY R AND RETURN 01298 IF(INFO(2).EQ.1)GO TO 523 01299 RTOL(1)=R*RTOL(1) 01300 ATOL(1)=R*ATOL(1) 01301 IDID=-2 01302 GO TO 527 01303 523 DO 524 I=1,NEQ 01304 RTOL(I)=R*RTOL(I) 01305 524 ATOL(I)=R*ATOL(I) 01306 IDID=-2 01307 GO TO 527 01308 525 CONTINUE 01309 C 01310 C COMPUTE MINIMUM STEPSIZE 01311 HMIN=4.0D0*UROUND*MAX(ABS(TN),ABS(TOUT)) 01312 C 01313 C TEST H VS. HMAX 01314 IF (INFO(7) .EQ. 0) GO TO 526 01315 RH = ABS(H)/RWORK(LHMAX) 01316 IF (RH .GT. 1.0D0) H = H/RH 01317 526 CONTINUE 01318 C 01319 CALL DDASTP(TN,Y,YPRIME,NEQ, 01320 * RES,JAC,H,RWORK(LWT),INFO(1),IDID,RPAR,IPAR, 01321 * RWORK(LPHI),RWORK(LDELTA),RWORK(LE), 01322 * RWORK(LWM),IWORK(LIWM), 01323 * RWORK(LALPHA),RWORK(LBETA),RWORK(LGAMMA), 01324 * RWORK(LPSI),RWORK(LSIGMA), 01325 * RWORK(LCJ),RWORK(LCJOLD),RWORK(LHOLD), 01326 * RWORK(LS),HMIN,RWORK(LROUND), 01327 * IWORK(LPHASE),IWORK(LJCALC),IWORK(LK), 01328 * IWORK(LKOLD),IWORK(LNS),INFO(10),NTEMP) 01329 527 IF(IDID.LT.0)GO TO 600 01330 C 01331 C-------------------------------------------------------- 01332 C THIS BLOCK HANDLES THE CASE OF A SUCCESSFUL RETURN 01333 C FROM DDASTP (IDID=1). TEST FOR STOP CONDITIONS. 01334 C-------------------------------------------------------- 01335 C 01336 IF(INFO(4).NE.0)GO TO 540 01337 IF(INFO(3).NE.0)GO TO 530 01338 IF((TN-TOUT)*H.LT.0.0D0)GO TO 500 01339 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ, 01340 * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI)) 01341 IDID=3 01342 T=TOUT 01343 GO TO 580 01344 530 IF((TN-TOUT)*H.GE.0.0D0)GO TO 535 01345 T=TN 01346 IDID=1 01347 GO TO 580 01348 535 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ, 01349 * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI)) 01350 IDID=3 01351 T=TOUT 01352 GO TO 580 01353 540 IF(INFO(3).NE.0)GO TO 550 01354 IF((TN-TOUT)*H.LT.0.0D0)GO TO 542 01355 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ, 01356 * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI)) 01357 T=TOUT 01358 IDID=3 01359 GO TO 580 01360 542 IF(ABS(TN-TSTOP).LE.100.0D0*UROUND* 01361 * (ABS(TN)+ABS(H)))GO TO 545 01362 TNEXT=TN+H 01363 IF((TNEXT-TSTOP)*H.LE.0.0D0)GO TO 500 01364 H=TSTOP-TN 01365 GO TO 500 01366 545 CALL DDATRP(TN,TSTOP,Y,YPRIME,NEQ, 01367 * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI)) 01368 IDID=2 01369 T=TSTOP 01370 GO TO 580 01371 550 IF((TN-TOUT)*H.GE.0.0D0)GO TO 555 01372 IF(ABS(TN-TSTOP).LE.100.0D0*UROUND*(ABS(TN)+ABS(H)))GO TO 552 01373 T=TN 01374 IDID=1 01375 GO TO 580 01376 552 CALL DDATRP(TN,TSTOP,Y,YPRIME,NEQ, 01377 * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI)) 01378 IDID=2 01379 T=TSTOP 01380 GO TO 580 01381 555 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ, 01382 * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI)) 01383 T=TOUT 01384 IDID=3 01385 GO TO 580 01386 C 01387 C-------------------------------------------------------- 01388 C ALL SUCCESSFUL RETURNS FROM DDASSL ARE MADE FROM 01389 C THIS BLOCK. 01390 C-------------------------------------------------------- 01391 C 01392 580 CONTINUE 01393 RWORK(LTN)=TN 01394 RWORK(LH)=H 01395 RETURN 01396 C 01397 C----------------------------------------------------------------------- 01398 C THIS BLOCK HANDLES ALL UNSUCCESSFUL 01399 C RETURNS OTHER THAN FOR ILLEGAL INPUT. 01400 C----------------------------------------------------------------------- 01401 C 01402 600 CONTINUE 01403 ITEMP=-IDID 01404 GO TO (610,620,630,690,690,640,650,660,670,675, 01405 * 680,685), ITEMP 01406 C 01407 C THE MAXIMUM NUMBER OF STEPS WAS TAKEN BEFORE 01408 C REACHING TOUT 01409 610 WRITE (XERN3, '(1P,D15.6)') TN 01410 CALL XERMSG ('SLATEC', 'DDASSL', 01411 * 'AT CURRENT T = ' // XERN3 // ' 500 STEPS TAKEN ON THIS ' // 01412 * 'CALL BEFORE REACHING TOUT', IDID, 1) 01413 GO TO 690 01414 C 01415 C TOO MUCH ACCURACY FOR MACHINE PRECISION 01416 620 WRITE (XERN3, '(1P,D15.6)') TN 01417 CALL XERMSG ('SLATEC', 'DDASSL', 01418 * 'AT T = ' // XERN3 // ' TOO MUCH ACCURACY REQUESTED FOR ' // 01419 * 'PRECISION OF MACHINE. RTOL AND ATOL WERE INCREASED TO ' // 01420 * 'APPROPRIATE VALUES', IDID, 1) 01421 GO TO 690 01422 C 01423 C WT(I) .LE. 0.0 FOR SOME I (NOT AT START OF PROBLEM) 01424 630 WRITE (XERN3, '(1P,D15.6)') TN 01425 CALL XERMSG ('SLATEC', 'DDASSL', 01426 * 'AT T = ' // XERN3 // ' SOME ELEMENT OF WT HAS BECOME .LE. ' // 01427 * '0.0', IDID, 1) 01428 GO TO 690 01429 C 01430 C ERROR TEST FAILED REPEATEDLY OR WITH H=HMIN 01431 640 WRITE (XERN3, '(1P,D15.6)') TN 01432 WRITE (XERN4, '(1P,D15.6)') H 01433 CALL XERMSG ('SLATEC', 'DDASSL', 01434 * 'AT T = ' // XERN3 // ' AND STEPSIZE H = ' // XERN4 // 01435 * ' THE ERROR TEST FAILED REPEATEDLY OR WITH ABS(H)=HMIN', 01436 * IDID, 1) 01437 GO TO 690 01438 C 01439 C CORRECTOR CONVERGENCE FAILED REPEATEDLY OR WITH H=HMIN 01440 650 WRITE (XERN3, '(1P,D15.6)') TN 01441 WRITE (XERN4, '(1P,D15.6)') H 01442 CALL XERMSG ('SLATEC', 'DDASSL', 01443 * 'AT T = ' // XERN3 // ' AND STEPSIZE H = ' // XERN4 // 01444 * ' THE CORRECTOR FAILED TO CONVERGE REPEATEDLY OR WITH ' // 01445 * 'ABS(H)=HMIN', IDID, 1) 01446 GO TO 690 01447 C 01448 C THE ITERATION MATRIX IS SINGULAR 01449 660 WRITE (XERN3, '(1P,D15.6)') TN 01450 WRITE (XERN4, '(1P,D15.6)') H 01451 CALL XERMSG ('SLATEC', 'DDASSL', 01452 * 'AT T = ' // XERN3 // ' AND STEPSIZE H = ' // XERN4 // 01453 * ' THE ITERATION MATRIX IS SINGULAR', IDID, 1) 01454 GO TO 690 01455 C 01456 C CORRECTOR FAILURE PRECEEDED BY ERROR TEST FAILURES. 01457 670 WRITE (XERN3, '(1P,D15.6)') TN 01458 WRITE (XERN4, '(1P,D15.6)') H 01459 CALL XERMSG ('SLATEC', 'DDASSL', 01460 * 'AT T = ' // XERN3 // ' AND STEPSIZE H = ' // XERN4 // 01461 * ' THE CORRECTOR COULD NOT CONVERGE. ALSO, THE ERROR TEST ' // 01462 * 'FAILED REPEATEDLY.', IDID, 1) 01463 GO TO 690 01464 C 01465 C CORRECTOR FAILURE BECAUSE IRES = -1 01466 675 WRITE (XERN3, '(1P,D15.6)') TN 01467 WRITE (XERN4, '(1P,D15.6)') H 01468 CALL XERMSG ('SLATEC', 'DDASSL', 01469 * 'AT T = ' // XERN3 // ' AND STEPSIZE H = ' // XERN4 // 01470 * ' THE CORRECTOR COULD NOT CONVERGE BECAUSE IRES WAS EQUAL ' // 01471 * 'TO MINUS ONE', IDID, 1) 01472 GO TO 690 01473 C 01474 C FAILURE BECAUSE IRES = -2 01475 680 WRITE (XERN3, '(1P,D15.6)') TN 01476 WRITE (XERN4, '(1P,D15.6)') H 01477 CALL XERMSG ('SLATEC', 'DDASSL', 01478 * 'AT T = ' // XERN3 // ' AND STEPSIZE H = ' // XERN4 // 01479 * ' IRES WAS EQUAL TO MINUS TWO', IDID, 1) 01480 GO TO 690 01481 C 01482 C FAILED TO COMPUTE INITIAL YPRIME 01483 685 WRITE (XERN3, '(1P,D15.6)') TN 01484 WRITE (XERN4, '(1P,D15.6)') HO 01485 CALL XERMSG ('SLATEC', 'DDASSL', 01486 * 'AT T = ' // XERN3 // ' AND STEPSIZE H = ' // XERN4 // 01487 * ' THE INITIAL YPRIME COULD NOT BE COMPUTED', IDID, 1) 01488 GO TO 690 01489 C 01490 690 CONTINUE 01491 INFO(1)=-1 01492 T=TN 01493 RWORK(LTN)=TN 01494 RWORK(LH)=H 01495 RETURN 01496 C 01497 C----------------------------------------------------------------------- 01498 C THIS BLOCK HANDLES ALL ERROR RETURNS DUE 01499 C TO ILLEGAL INPUT, AS DETECTED BEFORE CALLING 01500 C DDASTP. FIRST THE ERROR MESSAGE ROUTINE IS 01501 C CALLED. IF THIS HAPPENS TWICE IN 01502 C SUCCESSION, EXECUTION IS TERMINATED 01503 C 01504 C----------------------------------------------------------------------- 01505 701 CALL XERMSG ('SLATEC', 'DDASSL', 01506 * 'SOME ELEMENT OF INFO VECTOR IS NOT ZERO OR ONE', 1, 1) 01507 GO TO 750 01508 C 01509 702 WRITE (XERN1, '(I8)') NEQ 01510 CALL XERMSG ('SLATEC', 'DDASSL', 01511 * 'NEQ = ' // XERN1 // ' .LE. 0', 2, 1) 01512 GO TO 750 01513 C 01514 703 WRITE (XERN1, '(I8)') MXORD 01515 CALL XERMSG ('SLATEC', 'DDASSL', 01516 * 'MAXORD = ' // XERN1 // ' NOT IN RANGE', 3, 1) 01517 GO TO 750 01518 C 01519 704 WRITE (XERN1, '(I8)') LENRW 01520 WRITE (XERN2, '(I8)') LRW 01521 CALL XERMSG ('SLATEC', 'DDASSL', 01522 * 'RWORK LENGTH NEEDED, LENRW = ' // XERN1 // 01523 * ', EXCEEDS LRW = ' // XERN2, 4, 1) 01524 GO TO 750 01525 C 01526 705 WRITE (XERN1, '(I8)') LENIW 01527 WRITE (XERN2, '(I8)') LIW 01528 CALL XERMSG ('SLATEC', 'DDASSL', 01529 * 'IWORK LENGTH NEEDED, LENIW = ' // XERN1 // 01530 * ', EXCEEDS LIW = ' // XERN2, 5, 1) 01531 GO TO 750 01532 C 01533 706 CALL XERMSG ('SLATEC', 'DDASSL', 01534 * 'SOME ELEMENT OF RTOL IS .LT. 0', 6, 1) 01535 GO TO 750 01536 C 01537 707 CALL XERMSG ('SLATEC', 'DDASSL', 01538 * 'SOME ELEMENT OF ATOL IS .LT. 0', 7, 1) 01539 GO TO 750 01540 C 01541 708 CALL XERMSG ('SLATEC', 'DDASSL', 01542 * 'ALL ELEMENTS OF RTOL AND ATOL ARE ZERO', 8, 1) 01543 GO TO 750 01544 C 01545 709 WRITE (XERN3, '(1P,D15.6)') TSTOP 01546 WRITE (XERN4, '(1P,D15.6)') TOUT 01547 CALL XERMSG ('SLATEC', 'DDASSL', 01548 * 'INFO(4) = 1 AND TSTOP = ' // XERN3 // ' BEHIND TOUT = ' // 01549 * XERN4, 9, 1) 01550 GO TO 750 01551 C 01552 710 WRITE (XERN3, '(1P,D15.6)') HMAX 01553 CALL XERMSG ('SLATEC', 'DDASSL', 01554 * 'HMAX = ' // XERN3 // ' .LT. 0.0', 10, 1) 01555 GO TO 750 01556 C 01557 711 WRITE (XERN3, '(1P,D15.6)') TOUT 01558 WRITE (XERN4, '(1P,D15.6)') T 01559 CALL XERMSG ('SLATEC', 'DDASSL', 01560 * 'TOUT = ' // XERN3 // ' BEHIND T = ' // XERN4, 11, 1) 01561 GO TO 750 01562 C 01563 712 CALL XERMSG ('SLATEC', 'DDASSL', 01564 * 'INFO(8)=1 AND H0=0.0', 12, 1) 01565 GO TO 750 01566 C 01567 713 CALL XERMSG ('SLATEC', 'DDASSL', 01568 * 'SOME ELEMENT OF WT IS .LE. 0.0', 13, 1) 01569 GO TO 750 01570 C 01571 714 WRITE (XERN3, '(1P,D15.6)') TOUT 01572 WRITE (XERN4, '(1P,D15.6)') T 01573 CALL XERMSG ('SLATEC', 'DDASSL', 01574 * 'TOUT = ' // XERN3 // ' TOO CLOSE TO T = ' // XERN4 // 01575 * ' TO START INTEGRATION', 14, 1) 01576 GO TO 750 01577 C 01578 715 WRITE (XERN3, '(1P,D15.6)') TSTOP 01579 WRITE (XERN4, '(1P,D15.6)') T 01580 CALL XERMSG ('SLATEC', 'DDASSL', 01581 * 'INFO(4)=1 AND TSTOP = ' // XERN3 // ' BEHIND T = ' // XERN4, 01582 * 15, 1) 01583 GO TO 750 01584 C 01585 716 WRITE (XERN1, '(I8)') MXSTP 01586 CALL XERMSG ('SLATEC', 'DDASSL', 01587 * 'INFO(12)=1 AND MXSTP = ' // XERN1 // ' ILLEGAL.', 3, 1) 01588 GO TO 750 01589 C 01590 717 WRITE (XERN1, '(I8)') IWORK(LML) 01591 CALL XERMSG ('SLATEC', 'DDASSL', 01592 * 'ML = ' // XERN1 // ' ILLEGAL. EITHER .LT. 0 OR .GT. NEQ', 01593 * 17, 1) 01594 GO TO 750 01595 C 01596 718 WRITE (XERN1, '(I8)') IWORK(LMU) 01597 CALL XERMSG ('SLATEC', 'DDASSL', 01598 * 'MU = ' // XERN1 // ' ILLEGAL. EITHER .LT. 0 OR .GT. NEQ', 01599 * 18, 1) 01600 GO TO 750 01601 C 01602 719 WRITE (XERN3, '(1P,D15.6)') TOUT 01603 CALL XERMSG ('SLATEC', 'DDASSL', 01604 * 'TOUT = T = ' // XERN3, 19, 1) 01605 GO TO 750 01606 C 01607 750 IDID=-33 01608 IF(INFO(1).EQ.-1) THEN 01609 CALL XERMSG ('SLATEC', 'DDASSL', 01610 * 'REPEATED OCCURRENCES OF ILLEGAL INPUT$$' // 01611 * 'RUN TERMINATED. APPARENT INFINITE LOOP', -999, 2) 01612 ENDIF 01613 C 01614 INFO(1)=-1 01615 RETURN 01616 C-----------END OF SUBROUTINE DDASSL------------------------------------ 01617 END
1.7.1
