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