00001 C Work performed under the auspices of the U.S. Department of Energy 00002 C by Lawrence Livermore National Laboratory under contract number 00003 C W-7405-Eng-48. 00004 C 00005 SUBROUTINE DDASPK (RES, NEQ, T, Y, YPRIME, TOUT, INFO, RTOL, ATOL, 00006 * IDID, RWORK, LRW, IWORK, LIW, RPAR, IPAR, JAC, PSOL) 00007 C 00008 C***BEGIN PROLOGUE DDASPK 00009 C***DATE WRITTEN 890101 (YYMMDD) 00010 C***REVISION DATE 910624 00011 C***REVISION DATE 920929 (CJ in RES call, RES counter fix.) 00012 C***REVISION DATE 921215 (Warnings on poor iteration performance) 00013 C***REVISION DATE 921216 (NRMAX as optional input) 00014 C***REVISION DATE 930315 (Name change: DDINI to DDINIT) 00015 C***REVISION DATE 940822 (Replaced initial condition calculation) 00016 C***REVISION DATE 941101 (Added linesearch in I.C. calculations) 00017 C***REVISION DATE 941220 (Misc. corrections throughout) 00018 C***REVISION DATE 950125 (Added DINVWT routine) 00019 C***REVISION DATE 950714 (Misc. corrections throughout) 00020 C***REVISION DATE 950802 (Default NRMAX = 5, based on tests.) 00021 C***REVISION DATE 950808 (Optional error test added.) 00022 C***REVISION DATE 950814 (Added I.C. constraints and INFO(14)) 00023 C***REVISION DATE 950828 (Various minor corrections.) 00024 C***REVISION DATE 951006 (Corrected WT scaling in DFNRMK.) 00025 C***REVISION DATE 960129 (Corrected RL bug in DLINSD, DLINSK.) 00026 C***REVISION DATE 960301 (Added NONNEG to SAVE statement.) 00027 C***CATEGORY NO. I1A2 00028 C***KEYWORDS DIFFERENTIAL/ALGEBRAIC, BACKWARD DIFFERENTIATION FORMULAS, 00029 C IMPLICIT DIFFERENTIAL SYSTEMS, KRYLOV ITERATION 00030 C***AUTHORS Linda R. Petzold, Peter N. Brown, Alan C. Hindmarsh, and 00031 C Clement W. Ulrich 00032 C Center for Computational Sciences & Engineering, L-316 00033 C Lawrence Livermore National Laboratory 00034 C P.O. Box 808, 00035 C Livermore, CA 94551 00036 C***PURPOSE This code solves a system of differential/algebraic 00037 C equations of the form 00038 C G(t,y,y') = 0 , 00039 C using a combination of Backward Differentiation Formula 00040 C (BDF) methods and a choice of two linear system solution 00041 C methods: direct (dense or band) or Krylov (iterative). 00042 C This version is in double precision. 00043 C----------------------------------------------------------------------- 00044 C***DESCRIPTION 00045 C 00046 C *Usage: 00047 C 00048 C IMPLICIT DOUBLE PRECISION(A-H,O-Z) 00049 C INTEGER NEQ, INFO(N), IDID, LRW, LIW, IWORK(LIW), IPAR(*) 00050 C DOUBLE PRECISION T, Y(*), YPRIME(*), TOUT, RTOL(*), ATOL(*), 00051 C RWORK(LRW), RPAR(*) 00052 C EXTERNAL RES, JAC, PSOL 00053 C 00054 C CALL DDASPK (RES, NEQ, T, Y, YPRIME, TOUT, INFO, RTOL, ATOL, 00055 C * IDID, RWORK, LRW, IWORK, LIW, RPAR, IPAR, JAC, PSOL) 00056 C 00057 C Quantities which may be altered by the code are: 00058 C T, Y(*), YPRIME(*), INFO(*), RTOL, ATOL, IDID, RWORK(*), IWORK(*) 00059 C 00060 C 00061 C *Arguments: 00062 C 00063 C RES:EXT This is the name of a subroutine which you 00064 C provide to define the residual function G(t,y,y') 00065 C of the differential/algebraic system. 00066 C 00067 C NEQ:IN This is the number of equations in the system. 00068 C 00069 C T:INOUT This is the current value of the independent 00070 C variable. 00071 C 00072 C Y(*):INOUT This array contains the solution components at T. 00073 C 00074 C YPRIME(*):INOUT This array contains the derivatives of the solution 00075 C components at T. 00076 C 00077 C TOUT:IN This is a point at which a solution is desired. 00078 C 00079 C INFO(N):IN This is an integer array used to communicate details 00080 C of how the solution is to be carried out, such as 00081 C tolerance type, matrix structure, step size and 00082 C order limits, and choice of nonlinear system method. 00083 C N must be at least 20. 00084 C 00085 C RTOL,ATOL:INOUT These quantities represent absolute and relative 00086 C error tolerances (on local error) which you provide 00087 C to indicate how accurately you wish the solution to 00088 C be computed. You may choose them to be both scalars 00089 C or else both arrays of length NEQ. 00090 C 00091 C IDID:OUT This integer scalar is an indicator reporting what 00092 C the code did. You must monitor this variable to 00093 C decide what action to take next. 00094 C 00095 C RWORK:WORK A real work array of length LRW which provides the 00096 C code with needed storage space. 00097 C 00098 C LRW:IN The length of RWORK. 00099 C 00100 C IWORK:WORK An integer work array of length LIW which provides 00101 C the code with needed storage space. 00102 C 00103 C LIW:IN The length of IWORK. 00104 C 00105 C RPAR,IPAR:IN These are real and integer parameter arrays which 00106 C you can use for communication between your calling 00107 C program and the RES, JAC, and PSOL subroutines. 00108 C 00109 C JAC:EXT This is the name of a subroutine which you may 00110 C provide (optionally) for calculating Jacobian 00111 C (partial derivative) data involved in solving linear 00112 C systems within DDASPK. 00113 C 00114 C PSOL:EXT This is the name of a subroutine which you must 00115 C provide for solving linear systems if you selected 00116 C a Krylov method. The purpose of PSOL is to solve 00117 C linear systems involving a left preconditioner P. 00118 C 00119 C *Overview 00120 C 00121 C The DDASPK solver uses the backward differentiation formulas of 00122 C orders one through five to solve a system of the form G(t,y,y') = 0 00123 C for y = Y and y' = YPRIME. Values for Y and YPRIME at the initial 00124 C time must be given as input. These values should be consistent, 00125 C that is, if T, Y, YPRIME are the given initial values, they should 00126 C satisfy G(T,Y,YPRIME) = 0. However, if consistent values are not 00127 C known, in many cases you can have DDASPK solve for them -- see INFO(11). 00128 C (This and other options are described in more detail below.) 00129 C 00130 C Normally, DDASPK solves the system from T to TOUT. It is easy to 00131 C continue the solution to get results at additional TOUT. This is 00132 C the interval mode of operation. Intermediate results can also be 00133 C obtained easily by specifying INFO(3). 00134 C 00135 C On each step taken by DDASPK, a sequence of nonlinear algebraic 00136 C systems arises. These are solved by one of two types of 00137 C methods: 00138 C * a Newton iteration with a direct method for the linear 00139 C systems involved (INFO(12) = 0), or 00140 C * a Newton iteration with a preconditioned Krylov iterative 00141 C method for the linear systems involved (INFO(12) = 1). 00142 C 00143 C The direct method choices are dense and band matrix solvers, 00144 C with either a user-supplied or an internal difference quotient 00145 C Jacobian matrix, as specified by INFO(5) and INFO(6). 00146 C In the band case, INFO(6) = 1, you must supply half-bandwidths 00147 C in IWORK(1) and IWORK(2). 00148 C 00149 C The Krylov method is the Generalized Minimum Residual (GMRES) 00150 C method, in either complete or incomplete form, and with 00151 C scaling and preconditioning. The method is implemented 00152 C in an algorithm called SPIGMR. Certain options in the Krylov 00153 C method case are specified by INFO(13) and INFO(15). 00154 C 00155 C If the Krylov method is chosen, you may supply a pair of routines, 00156 C JAC and PSOL, to apply preconditioning to the linear system. 00157 C If the system is A*x = b, the matrix is A = dG/dY + CJ*dG/dYPRIME 00158 C (of order NEQ). This system can then be preconditioned in the form 00159 C (P-inverse)*A*x = (P-inverse)*b, with left preconditioner P. 00160 C (DDASPK does not allow right preconditioning.) 00161 C Then the Krylov method is applied to this altered, but equivalent, 00162 C linear system, hopefully with much better performance than without 00163 C preconditioning. (In addition, a diagonal scaling matrix based on 00164 C the tolerances is also introduced into the altered system.) 00165 C 00166 C The JAC routine evaluates any data needed for solving systems 00167 C with coefficient matrix P, and PSOL carries out that solution. 00168 C In any case, in order to improve convergence, you should try to 00169 C make P approximate the matrix A as much as possible, while keeping 00170 C the system P*x = b reasonably easy and inexpensive to solve for x, 00171 C given a vector b. 00172 C 00173 C 00174 C *Description 00175 C 00176 C------INPUT - WHAT TO DO ON THE FIRST CALL TO DDASPK------------------- 00177 C 00178 C 00179 C The first call of the code is defined to be the start of each new 00180 C problem. Read through the descriptions of all the following items, 00181 C provide sufficient storage space for designated arrays, set 00182 C appropriate variables for the initialization of the problem, and 00183 C give information about how you want the problem to be solved. 00184 C 00185 C 00186 C RES -- Provide a subroutine of the form 00187 C 00188 C SUBROUTINE RES (T, Y, YPRIME, CJ, DELTA, IRES, RPAR, IPAR) 00189 C 00190 C to define the system of differential/algebraic 00191 C equations which is to be solved. For the given values 00192 C of T, Y and YPRIME, the subroutine should return 00193 C the residual of the differential/algebraic system 00194 C DELTA = G(T,Y,YPRIME) 00195 C DELTA is a vector of length NEQ which is output from RES. 00196 C 00197 C Subroutine RES must not alter T, Y, YPRIME, or CJ. 00198 C You must declare the name RES in an EXTERNAL 00199 C statement in your program that calls DDASPK. 00200 C You must dimension Y, YPRIME, and DELTA in RES. 00201 C 00202 C The input argument CJ can be ignored, or used to rescale 00203 C constraint equations in the system (see Ref. 2, p. 145). 00204 C Note: In this respect, DDASPK is not downward-compatible 00205 C with DDASSL, which does not have the RES argument CJ. 00206 C 00207 C IRES is an integer flag which is always equal to zero 00208 C on input. Subroutine RES should alter IRES only if it 00209 C encounters an illegal value of Y or a stop condition. 00210 C Set IRES = -1 if an input value is illegal, and DDASPK 00211 C will try to solve the problem without getting IRES = -1. 00212 C If IRES = -2, DDASPK will return control to the calling 00213 C program with IDID = -11. 00214 C 00215 C RPAR and IPAR are real and integer parameter arrays which 00216 C you can use for communication between your calling program 00217 C and subroutine RES. They are not altered by DDASPK. If you 00218 C do not need RPAR or IPAR, ignore these parameters by treat- 00219 C ing them as dummy arguments. If you do choose to use them, 00220 C dimension them in your calling program and in RES as arrays 00221 C of appropriate length. 00222 C 00223 C NEQ -- Set it to the number of equations in the system (NEQ .GE. 1). 00224 C 00225 C T -- Set it to the initial point of the integration. (T must be 00226 C a variable.) 00227 C 00228 C Y(*) -- Set this array to the initial values of the NEQ solution 00229 C components at the initial point. You must dimension Y of 00230 C length at least NEQ in your calling program. 00231 C 00232 C YPRIME(*) -- Set this array to the initial values of the NEQ first 00233 C derivatives of the solution components at the initial 00234 C point. You must dimension YPRIME at least NEQ in your 00235 C calling program. 00236 C 00237 C TOUT - Set it to the first point at which a solution is desired. 00238 C You cannot take TOUT = T. Integration either forward in T 00239 C (TOUT .GT. T) or backward in T (TOUT .LT. T) is permitted. 00240 C 00241 C The code advances the solution from T to TOUT using step 00242 C sizes which are automatically selected so as to achieve the 00243 C desired accuracy. If you wish, the code will return with the 00244 C solution and its derivative at intermediate steps (the 00245 C intermediate-output mode) so that you can monitor them, 00246 C but you still must provide TOUT in accord with the basic 00247 C aim of the code. 00248 C 00249 C The first step taken by the code is a critical one because 00250 C it must reflect how fast the solution changes near the 00251 C initial point. The code automatically selects an initial 00252 C step size which is practically always suitable for the 00253 C problem. By using the fact that the code will not step past 00254 C TOUT in the first step, you could, if necessary, restrict the 00255 C length of the initial step. 00256 C 00257 C For some problems it may not be permissible to integrate 00258 C past a point TSTOP, because a discontinuity occurs there 00259 C or the solution or its derivative is not defined beyond 00260 C TSTOP. When you have declared a TSTOP point (see INFO(4) 00261 C and RWORK(1)), you have told the code not to integrate past 00262 C TSTOP. In this case any tout beyond TSTOP is invalid input. 00263 C 00264 C INFO(*) - Use the INFO array to give the code more details about 00265 C how you want your problem solved. This array should be 00266 C dimensioned of length 20, though DDASPK uses only the 00267 C first 15 entries. You must respond to all of the following 00268 C items, which are arranged as questions. The simplest use 00269 C of DDASPK corresponds to setting all entries of INFO to 0. 00270 C 00271 C INFO(1) - This parameter enables the code to initialize itself. 00272 C You must set it to indicate the start of every new 00273 C problem. 00274 C 00275 C **** Is this the first call for this problem ... 00276 C yes - set INFO(1) = 0 00277 C no - not applicable here. 00278 C See below for continuation calls. **** 00279 C 00280 C INFO(2) - How much accuracy you want of your solution 00281 C is specified by the error tolerances RTOL and ATOL. 00282 C The simplest use is to take them both to be scalars. 00283 C To obtain more flexibility, they can both be arrays. 00284 C The code must be told your choice. 00285 C 00286 C **** Are both error tolerances RTOL, ATOL scalars ... 00287 C yes - set INFO(2) = 0 00288 C and input scalars for both RTOL and ATOL 00289 C no - set INFO(2) = 1 00290 C and input arrays for both RTOL and ATOL **** 00291 C 00292 C INFO(3) - The code integrates from T in the direction of TOUT 00293 C by steps. If you wish, it will return the computed 00294 C solution and derivative at the next intermediate step 00295 C (the intermediate-output mode) or TOUT, whichever comes 00296 C first. This is a good way to proceed if you want to 00297 C see the behavior of the solution. If you must have 00298 C solutions at a great many specific TOUT points, this 00299 C code will compute them efficiently. 00300 C 00301 C **** Do you want the solution only at 00302 C TOUT (and not at the next intermediate step) ... 00303 C yes - set INFO(3) = 0 00304 C no - set INFO(3) = 1 **** 00305 C 00306 C INFO(4) - To handle solutions at a great many specific 00307 C values TOUT efficiently, this code may integrate past 00308 C TOUT and interpolate to obtain the result at TOUT. 00309 C Sometimes it is not possible to integrate beyond some 00310 C point TSTOP because the equation changes there or it is 00311 C not defined past TSTOP. Then you must tell the code 00312 C this stop condition. 00313 C 00314 C **** Can the integration be carried out without any 00315 C restrictions on the independent variable T ... 00316 C yes - set INFO(4) = 0 00317 C no - set INFO(4) = 1 00318 C and define the stopping point TSTOP by 00319 C setting RWORK(1) = TSTOP **** 00320 C 00321 C INFO(5) - used only when INFO(12) = 0 (direct methods). 00322 C To solve differential/algebraic systems you may wish 00323 C to use a matrix of partial derivatives of the 00324 C system of differential equations. If you do not 00325 C provide a subroutine to evaluate it analytically (see 00326 C description of the item JAC in the call list), it will 00327 C be approximated by numerical differencing in this code. 00328 C Although it is less trouble for you to have the code 00329 C compute partial derivatives by numerical differencing, 00330 C the solution will be more reliable if you provide the 00331 C derivatives via JAC. Usually numerical differencing is 00332 C more costly than evaluating derivatives in JAC, but 00333 C sometimes it is not - this depends on your problem. 00334 C 00335 C **** Do you want the code to evaluate the partial deriv- 00336 C atives automatically by numerical differences ... 00337 C yes - set INFO(5) = 0 00338 C no - set INFO(5) = 1 00339 C and provide subroutine JAC for evaluating the 00340 C matrix of partial derivatives **** 00341 C 00342 C INFO(6) - used only when INFO(12) = 0 (direct methods). 00343 C DDASPK will perform much better if the matrix of 00344 C partial derivatives, dG/dY + CJ*dG/dYPRIME (here CJ is 00345 C a scalar determined by DDASPK), is banded and the code 00346 C is told this. In this case, the storage needed will be 00347 C greatly reduced, numerical differencing will be performed 00348 C much cheaper, and a number of important algorithms will 00349 C execute much faster. The differential equation is said 00350 C to have half-bandwidths ML (lower) and MU (upper) if 00351 C equation i involves only unknowns Y(j) with 00352 C i-ML .le. j .le. i+MU . 00353 C For all i=1,2,...,NEQ. Thus, ML and MU are the widths 00354 C of the lower and upper parts of the band, respectively, 00355 C with the main diagonal being excluded. If you do not 00356 C indicate that the equation has a banded matrix of partial 00357 C derivatives the code works with a full matrix of NEQ**2 00358 C elements (stored in the conventional way). Computations 00359 C with banded matrices cost less time and storage than with 00360 C full matrices if 2*ML+MU .lt. NEQ. If you tell the 00361 C code that the matrix of partial derivatives has a banded 00362 C structure and you want to provide subroutine JAC to 00363 C compute the partial derivatives, then you must be careful 00364 C to store the elements of the matrix in the special form 00365 C indicated in the description of JAC. 00366 C 00367 C **** Do you want to solve the problem using a full (dense) 00368 C matrix (and not a special banded structure) ... 00369 C yes - set INFO(6) = 0 00370 C no - set INFO(6) = 1 00371 C and provide the lower (ML) and upper (MU) 00372 C bandwidths by setting 00373 C IWORK(1)=ML 00374 C IWORK(2)=MU **** 00375 C 00376 C INFO(7) - You can specify a maximum (absolute value of) 00377 C stepsize, so that the code will avoid passing over very 00378 C large regions. 00379 C 00380 C **** Do you want the code to decide on its own the maximum 00381 C stepsize ... 00382 C yes - set INFO(7) = 0 00383 C no - set INFO(7) = 1 00384 C and define HMAX by setting 00385 C RWORK(2) = HMAX **** 00386 C 00387 C INFO(8) - Differential/algebraic problems may occasionally 00388 C suffer from severe scaling difficulties on the first 00389 C step. If you know a great deal about the scaling of 00390 C your problem, you can help to alleviate this problem 00391 C by specifying an initial stepsize H0. 00392 C 00393 C **** Do you want the code to define its own initial 00394 C stepsize ... 00395 C yes - set INFO(8) = 0 00396 C no - set INFO(8) = 1 00397 C and define H0 by setting 00398 C RWORK(3) = H0 **** 00399 C 00400 C INFO(9) - If storage is a severe problem, you can save some 00401 C storage by restricting the maximum method order MAXORD. 00402 C The default value is 5. For each order decrease below 5, 00403 C the code requires NEQ fewer locations, but it is likely 00404 C to be slower. In any case, you must have 00405 C 1 .le. MAXORD .le. 5. 00406 C **** Do you want the maximum order to default to 5 ... 00407 C yes - set INFO(9) = 0 00408 C no - set INFO(9) = 1 00409 C and define MAXORD by setting 00410 C IWORK(3) = MAXORD **** 00411 C 00412 C INFO(10) - If you know that certain components of the 00413 C solutions to your equations are always nonnegative 00414 C (or nonpositive), it may help to set this 00415 C parameter. There are three options that are 00416 C available: 00417 C 1. To have constraint checking only in the initial 00418 C condition calculation. 00419 C 2. To enforce nonnegativity in Y during the integration. 00420 C 3. To enforce both options 1 and 2. 00421 C 00422 C When selecting option 2 or 3, it is probably best to try the 00423 C code without using this option first, and only use 00424 C this option if that does not work very well. 00425 C 00426 C **** Do you want the code to solve the problem without 00427 C invoking any special inequality constraints ... 00428 C yes - set INFO(10) = 0 00429 C no - set INFO(10) = 1 to have option 1 enforced 00430 C no - set INFO(10) = 2 to have option 2 enforced 00431 C no - set INFO(10) = 3 to have option 3 enforced **** 00432 C 00433 C If you have specified INFO(10) = 1 or 3, then you 00434 C will also need to identify how each component of Y 00435 C in the initial condition calculation is constrained. 00436 C You must set: 00437 C IWORK(40+I) = +1 if Y(I) must be .GE. 0, 00438 C IWORK(40+I) = +2 if Y(I) must be .GT. 0, 00439 C IWORK(40+I) = -1 if Y(I) must be .LE. 0, while 00440 C IWORK(40+I) = -2 if Y(I) must be .LT. 0, while 00441 C IWORK(40+I) = 0 if Y(I) is not constrained. 00442 C 00443 C INFO(11) - DDASPK normally requires the initial T, Y, and 00444 C YPRIME to be consistent. That is, you must have 00445 C G(T,Y,YPRIME) = 0 at the initial T. If you do not know 00446 C the initial conditions precisely, in some cases 00447 C DDASPK may be able to compute it. 00448 C 00449 C Denoting the differential variables in Y by Y_d 00450 C and the algebraic variables by Y_a, DDASPK can solve 00451 C one of two initialization problems: 00452 C 1. Given Y_d, calculate Y_a and Y'_d, or 00453 C 2. Given Y', calculate Y. 00454 C In either case, initial values for the given 00455 C components are input, and initial guesses for 00456 C the unknown components must also be provided as input. 00457 C 00458 C **** Are the initial T, Y, YPRIME consistent ... 00459 C 00460 C yes - set INFO(11) = 0 00461 C no - set INFO(11) = 1 to calculate option 1 above, 00462 C or set INFO(11) = 2 to calculate option 2 **** 00463 C 00464 C If you have specified INFO(11) = 1, then you 00465 C will also need to identify which are the 00466 C differential and which are the algebraic 00467 C components (algebraic components are components 00468 C whose derivatives do not appear explicitly 00469 C in the function G(T,Y,YPRIME)). You must set: 00470 C IWORK(LID+I) = +1 if Y(I) is a differential variable 00471 C IWORK(LID+I) = -1 if Y(I) is an algebraic variable, 00472 C where LID = 40 if INFO(10) = 0 or 2 and LID = 40+NEQ 00473 C if INFO(10) = 1 or 3. 00474 C 00475 C INFO(12) - Except for the addition of the RES argument CJ, 00476 C DDASPK by default is downward-compatible with DDASSL, 00477 C which uses only direct (dense or band) methods to solve 00478 C the linear systems involved. You must set INFO(12) to 00479 C indicate whether you want the direct methods or the 00480 C Krylov iterative method. 00481 C **** Do you want DDASPK to use standard direct methods 00482 C (dense or band) or the Krylov (iterative) method ... 00483 C direct methods - set INFO(12) = 0. 00484 C Krylov method - set INFO(12) = 1, 00485 C and check the settings of INFO(13) and INFO(15). 00486 C 00487 C INFO(13) - used when INFO(12) = 1 (Krylov methods). 00488 C DDASPK uses scalars MAXL, KMP, NRMAX, and EPLI for the 00489 C iterative solution of linear systems. INFO(13) allows 00490 C you to override the default values of these parameters. 00491 C These parameters and their defaults are as follows: 00492 C MAXL = maximum number of iterations in the SPIGMR 00493 C algorithm (MAXL .le. NEQ). The default is 00494 C MAXL = MIN(5,NEQ). 00495 C KMP = number of vectors on which orthogonalization is 00496 C done in the SPIGMR algorithm. The default is 00497 C KMP = MAXL, which corresponds to complete GMRES 00498 C iteration, as opposed to the incomplete form. 00499 C NRMAX = maximum number of restarts of the SPIGMR 00500 C algorithm per nonlinear iteration. The default is 00501 C NRMAX = 5. 00502 C EPLI = convergence test constant in SPIGMR algorithm. 00503 C The default is EPLI = 0.05. 00504 C Note that the length of RWORK depends on both MAXL 00505 C and KMP. See the definition of LRW below. 00506 C **** Are MAXL, KMP, and EPLI to be given their 00507 C default values ... 00508 C yes - set INFO(13) = 0 00509 C no - set INFO(13) = 1, 00510 C and set all of the following: 00511 C IWORK(24) = MAXL (1 .le. MAXL .le. NEQ) 00512 C IWORK(25) = KMP (1 .le. KMP .le. MAXL) 00513 C IWORK(26) = NRMAX (NRMAX .ge. 0) 00514 C RWORK(10) = EPLI (0 .lt. EPLI .lt. 1.0) **** 00515 C 00516 C INFO(14) - used with INFO(11) > 0 (initial condition 00517 C calculation is requested). In this case, you may 00518 C request control to be returned to the calling program 00519 C immediately after the initial condition calculation, 00520 C before proceeding to the integration of the system 00521 C (e.g. to examine the computed Y and YPRIME). 00522 C If this is done, and if the initialization succeeded 00523 C (IDID = 4), you should reset INFO(11) to 0 for the 00524 C next call, to prevent the solver from repeating the 00525 C initialization (and to avoid an infinite loop). 00526 C **** Do you want to proceed to the integration after 00527 C the initial condition calculation is done ... 00528 C yes - set INFO(14) = 0 00529 C no - set INFO(14) = 1 **** 00530 C 00531 C INFO(15) - used when INFO(12) = 1 (Krylov methods). 00532 C When using preconditioning in the Krylov method, 00533 C you must supply a subroutine, PSOL, which solves the 00534 C associated linear systems using P. 00535 C The usage of DDASPK is simpler if PSOL can carry out 00536 C the solution without any prior calculation of data. 00537 C However, if some partial derivative data is to be 00538 C calculated in advance and used repeatedly in PSOL, 00539 C then you must supply a JAC routine to do this, 00540 C and set INFO(15) to indicate that JAC is to be called 00541 C for this purpose. For example, P might be an 00542 C approximation to a part of the matrix A which can be 00543 C calculated and LU-factored for repeated solutions of 00544 C the preconditioner system. The arrays WP and IWP 00545 C (described under JAC and PSOL) can be used to 00546 C communicate data between JAC and PSOL. 00547 C **** Does PSOL operate with no prior preparation ... 00548 C yes - set INFO(15) = 0 (no JAC routine) 00549 C no - set INFO(15) = 1 00550 C and supply a JAC routine to evaluate and 00551 C preprocess any required Jacobian data. **** 00552 C 00553 C INFO(16) - option to exclude algebraic variables from 00554 C the error test. 00555 C **** Do you wish to control errors locally on 00556 C all the variables... 00557 C yes - set INFO(16) = 0 00558 C no - set INFO(16) = 1 00559 C If you have specified INFO(16) = 1, then you 00560 C will also need to identify which are the 00561 C differential and which are the algebraic 00562 C components (algebraic components are components 00563 C whose derivatives do not appear explicitly 00564 C in the function G(T,Y,YPRIME)). You must set: 00565 C IWORK(LID+I) = +1 if Y(I) is a differential 00566 C variable, and 00567 C IWORK(LID+I) = -1 if Y(I) is an algebraic 00568 C variable, 00569 C where LID = 40 if INFO(10) = 0 or 2 and 00570 C LID = 40 + NEQ if INFO(10) = 1 or 3. 00571 C 00572 C INFO(17) - used when INFO(11) > 0 (DDASPK is to do an 00573 C initial condition calculation). 00574 C DDASPK uses several heuristic control quantities in the 00575 C initial condition calculation. They have default values, 00576 C but can also be set by the user using INFO(17). 00577 C These parameters and their defaults are as follows: 00578 C MXNIT = maximum number of Newton iterations 00579 C per Jacobian or preconditioner evaluation. 00580 C The default is: 00581 C MXNIT = 5 in the direct case (INFO(12) = 0), and 00582 C MXNIT = 15 in the Krylov case (INFO(12) = 1). 00583 C MXNJ = maximum number of Jacobian or preconditioner 00584 C evaluations. The default is: 00585 C MXNJ = 6 in the direct case (INFO(12) = 0), and 00586 C MXNJ = 2 in the Krylov case (INFO(12) = 1). 00587 C MXNH = maximum number of values of the artificial 00588 C stepsize parameter H to be tried if INFO(11) = 1. 00589 C The default is MXNH = 5. 00590 C NOTE: the maximum number of Newton iterations 00591 C allowed in all is MXNIT*MXNJ*MXNH if INFO(11) = 1, 00592 C and MXNIT*MXNJ if INFO(11) = 2. 00593 C LSOFF = flag to turn off the linesearch algorithm 00594 C (LSOFF = 0 means linesearch is on, LSOFF = 1 means 00595 C it is turned off). The default is LSOFF = 0. 00596 C STPTOL = minimum scaled step in linesearch algorithm. 00597 C The default is STPTOL = (unit roundoff)**(2/3). 00598 C EPINIT = swing factor in the Newton iteration convergence 00599 C test. The test is applied to the residual vector, 00600 C premultiplied by the approximate Jacobian (in the 00601 C direct case) or the preconditioner (in the Krylov 00602 C case). For convergence, the weighted RMS norm of 00603 C this vector (scaled by the error weights) must be 00604 C less than EPINIT*EPCON, where EPCON = .33 is the 00605 C analogous test constant used in the time steps. 00606 C The default is EPINIT = .01. 00607 C **** Are the initial condition heuristic controls to be 00608 C given their default values... 00609 C yes - set INFO(17) = 0 00610 C no - set INFO(17) = 1, 00611 C and set all of the following: 00612 C IWORK(32) = MXNIT (.GT. 0) 00613 C IWORK(33) = MXNJ (.GT. 0) 00614 C IWORK(34) = MXNH (.GT. 0) 00615 C IWORK(35) = LSOFF ( = 0 or 1) 00616 C RWORK(14) = STPTOL (.GT. 0.0) 00617 C RWORK(15) = EPINIT (.GT. 0.0) **** 00618 C 00619 C INFO(18) - option to get extra printing in initial condition 00620 C calculation. 00621 C **** Do you wish to have extra printing... 00622 C no - set INFO(18) = 0 00623 C yes - set INFO(18) = 1 for minimal printing, or 00624 C set INFO(18) = 2 for full printing. 00625 C If you have specified INFO(18) .ge. 1, data 00626 C will be printed with the error handler routines. 00627 C To print to a non-default unit number L, include 00628 C the line CALL XSETUN(L) in your program. **** 00629 C 00630 C RTOL, ATOL -- You must assign relative (RTOL) and absolute (ATOL) 00631 C error tolerances to tell the code how accurately you 00632 C want the solution to be computed. They must be defined 00633 C as variables because the code may change them. 00634 C you have two choices -- 00635 C Both RTOL and ATOL are scalars (INFO(2) = 0), or 00636 C both RTOL and ATOL are vectors (INFO(2) = 1). 00637 C In either case all components must be non-negative. 00638 C 00639 C The tolerances are used by the code in a local error 00640 C test at each step which requires roughly that 00641 C abs(local error in Y(i)) .le. EWT(i) , 00642 C where EWT(i) = RTOL*abs(Y(i)) + ATOL is an error weight 00643 C quantity, for each vector component. 00644 C (More specifically, a root-mean-square norm is used to 00645 C measure the size of vectors, and the error test uses the 00646 C magnitude of the solution at the beginning of the step.) 00647 C 00648 C The true (global) error is the difference between the 00649 C true solution of the initial value problem and the 00650 C computed approximation. Practically all present day 00651 C codes, including this one, control the local error at 00652 C each step and do not even attempt to control the global 00653 C error directly. 00654 C 00655 C Usually, but not always, the true accuracy of 00656 C the computed Y is comparable to the error tolerances. 00657 C This code will usually, but not always, deliver a more 00658 C accurate solution if you reduce the tolerances and 00659 C integrate again. By comparing two such solutions you 00660 C can get a fairly reliable idea of the true error in the 00661 C solution at the larger tolerances. 00662 C 00663 C Setting ATOL = 0. results in a pure relative error test 00664 C on that component. Setting RTOL = 0. results in a pure 00665 C absolute error test on that component. A mixed test 00666 C with non-zero RTOL and ATOL corresponds roughly to a 00667 C relative error test when the solution component is 00668 C much bigger than ATOL and to an absolute error test 00669 C when the solution component is smaller than the 00670 C threshold ATOL. 00671 C 00672 C The code will not attempt to compute a solution at an 00673 C accuracy unreasonable for the machine being used. It 00674 C will advise you if you ask for too much accuracy and 00675 C inform you as to the maximum accuracy it believes 00676 C possible. 00677 C 00678 C RWORK(*) -- a real work array, which should be dimensioned in your 00679 C calling program with a length equal to the value of 00680 C LRW (or greater). 00681 C 00682 C LRW -- Set it to the declared length of the RWORK array. The 00683 C minimum length depends on the options you have selected, 00684 C given by a base value plus additional storage as described 00685 C below. 00686 C 00687 C If INFO(12) = 0 (standard direct method), the base value is 00688 C base = 50 + max(MAXORD+4,7)*NEQ. 00689 C The default value is MAXORD = 5 (see INFO(9)). With the 00690 C default MAXORD, base = 50 + 9*NEQ. 00691 C Additional storage must be added to the base value for 00692 C any or all of the following options: 00693 C if INFO(6) = 0 (dense matrix), add NEQ**2 00694 C if INFO(6) = 1 (banded matrix), then 00695 C if INFO(5) = 0, add (2*ML+MU+1)*NEQ + 2*(NEQ/(ML+MU+1)+1), 00696 C if INFO(5) = 1, add (2*ML+MU+1)*NEQ, 00697 C if INFO(16) = 1, add NEQ. 00698 C 00699 C If INFO(12) = 1 (Krylov method), the base value is 00700 C base = 50 + (MAXORD+5)*NEQ + (MAXL+3+MIN0(1,MAXL-KMP))*NEQ + 00701 C + (MAXL+3)*MAXL + 1 + LENWP. 00702 C See PSOL for description of LENWP. The default values are: 00703 C MAXORD = 5 (see INFO(9)), MAXL = min(5,NEQ) and KMP = MAXL 00704 C (see INFO(13)). 00705 C With the default values for MAXORD, MAXL and KMP, 00706 C base = 91 + 18*NEQ + LENWP. 00707 C Additional storage must be added to the base value for 00708 C any or all of the following options: 00709 C if INFO(16) = 1, add NEQ. 00710 C 00711 C 00712 C IWORK(*) -- an integer work array, which should be dimensioned in 00713 C your calling program with a length equal to the value 00714 C of LIW (or greater). 00715 C 00716 C LIW -- Set it to the declared length of the IWORK array. The 00717 C minimum length depends on the options you have selected, 00718 C given by a base value plus additional storage as described 00719 C below. 00720 C 00721 C If INFO(12) = 0 (standard direct method), the base value is 00722 C base = 40 + NEQ. 00723 C IF INFO(10) = 1 or 3, add NEQ to the base value. 00724 C If INFO(11) = 1 or INFO(16) =1, add NEQ to the base value. 00725 C 00726 C If INFO(12) = 1 (Krylov method), the base value is 00727 C base = 40 + LENIWP. 00728 C See PSOL for description of LENIWP. 00729 C IF INFO(10) = 1 or 3, add NEQ to the base value. 00730 C If INFO(11) = 1 or INFO(16) = 1, add NEQ to the base value. 00731 C 00732 C 00733 C RPAR, IPAR -- These are arrays of double precision and integer type, 00734 C respectively, which are available for you to use 00735 C for communication between your program that calls 00736 C DDASPK and the RES subroutine (and the JAC and PSOL 00737 C subroutines). They are not altered by DDASPK. 00738 C If you do not need RPAR or IPAR, ignore these 00739 C parameters by treating them as dummy arguments. 00740 C If you do choose to use them, dimension them in 00741 C your calling program and in RES (and in JAC and PSOL) 00742 C as arrays of appropriate length. 00743 C 00744 C JAC -- This is the name of a routine that you may supply 00745 C (optionally) that relates to the Jacobian matrix of the 00746 C nonlinear system that the code must solve at each T step. 00747 C The role of JAC (and its call sequence) depends on whether 00748 C a direct (INFO(12) = 0) or Krylov (INFO(12) = 1) method 00749 C is selected. 00750 C 00751 C **** INFO(12) = 0 (direct methods): 00752 C If you are letting the code generate partial derivatives 00753 C numerically (INFO(5) = 0), then JAC can be absent 00754 C (or perhaps a dummy routine to satisfy the loader). 00755 C Otherwise you must supply a JAC routine to compute 00756 C the matrix A = dG/dY + CJ*dG/dYPRIME. It must have 00757 C the form 00758 C 00759 C SUBROUTINE JAC (T, Y, YPRIME, PD, CJ, RPAR, IPAR) 00760 C 00761 C The JAC routine must dimension Y, YPRIME, and PD (and RPAR 00762 C and IPAR if used). CJ is a scalar which is input to JAC. 00763 C For the given values of T, Y, and YPRIME, the JAC routine 00764 C must evaluate the nonzero elements of the matrix A, and 00765 C store these values in the array PD. The elements of PD are 00766 C set to zero before each call to JAC, so that only nonzero 00767 C elements need to be defined. 00768 C The way you store the elements into the PD array depends 00769 C on the structure of the matrix indicated by INFO(6). 00770 C *** INFO(6) = 0 (full or dense matrix) *** 00771 C Give PD a first dimension of NEQ. When you evaluate the 00772 C nonzero partial derivatives of equation i (i.e. of G(i)) 00773 C with respect to component j (of Y and YPRIME), you must 00774 C store the element in PD according to 00775 C PD(i,j) = dG(i)/dY(j) + CJ*dG(i)/dYPRIME(j). 00776 C *** INFO(6) = 1 (banded matrix with half-bandwidths ML, MU 00777 C as described under INFO(6)) *** 00778 C Give PD a first dimension of 2*ML+MU+1. When you 00779 C evaluate the nonzero partial derivatives of equation i 00780 C (i.e. of G(i)) with respect to component j (of Y and 00781 C YPRIME), you must store the element in PD according to 00782 C IROW = i - j + ML + MU + 1 00783 C PD(IROW,j) = dG(i)/dY(j) + CJ*dG(i)/dYPRIME(j). 00784 C 00785 C **** INFO(12) = 1 (Krylov method): 00786 C If you are not calculating Jacobian data in advance for use 00787 C in PSOL (INFO(15) = 0), JAC can be absent (or perhaps a 00788 C dummy routine to satisfy the loader). Otherwise, you may 00789 C supply a JAC routine to compute and preprocess any parts of 00790 C of the Jacobian matrix A = dG/dY + CJ*dG/dYPRIME that are 00791 C involved in the preconditioner matrix P. 00792 C It is to have the form 00793 C 00794 C SUBROUTINE JAC (RES, IRES, NEQ, T, Y, YPRIME, REWT, SAVR, 00795 C WK, H, CJ, WP, IWP, IER, RPAR, IPAR) 00796 C 00797 C The JAC routine must dimension Y, YPRIME, REWT, SAVR, WK, 00798 C and (if used) WP, IWP, RPAR, and IPAR. 00799 C The Y, YPRIME, and SAVR arrays contain the current values 00800 C of Y, YPRIME, and the residual G, respectively. 00801 C The array WK is work space of length NEQ. 00802 C H is the step size. CJ is a scalar, input to JAC, that is 00803 C normally proportional to 1/H. REWT is an array of 00804 C reciprocal error weights, 1/EWT(i), where EWT(i) is 00805 C RTOL*abs(Y(i)) + ATOL (unless you supplied routine DDAWTS 00806 C instead), for use in JAC if needed. For example, if JAC 00807 C computes difference quotient approximations to partial 00808 C derivatives, the REWT array may be useful in setting the 00809 C increments used. The JAC routine should do any 00810 C factorization operations called for, in preparation for 00811 C solving linear systems in PSOL. The matrix P should 00812 C be an approximation to the Jacobian, 00813 C A = dG/dY + CJ*dG/dYPRIME. 00814 C 00815 C WP and IWP are real and integer work arrays which you may 00816 C use for communication between your JAC routine and your 00817 C PSOL routine. These may be used to store elements of the 00818 C preconditioner P, or related matrix data (such as factored 00819 C forms). They are not altered by DDASPK. 00820 C If you do not need WP or IWP, ignore these parameters by 00821 C treating them as dummy arguments. If you do use them, 00822 C dimension them appropriately in your JAC and PSOL routines. 00823 C See the PSOL description for instructions on setting 00824 C the lengths of WP and IWP. 00825 C 00826 C On return, JAC should set the error flag IER as follows.. 00827 C IER = 0 if JAC was successful, 00828 C IER .ne. 0 if JAC was unsuccessful (e.g. if Y or YPRIME 00829 C was illegal, or a singular matrix is found). 00830 C (If IER .ne. 0, a smaller stepsize will be tried.) 00831 C IER = 0 on entry to JAC, so need be reset only on a failure. 00832 C If RES is used within JAC, then a nonzero value of IRES will 00833 C override any nonzero value of IER (see the RES description). 00834 C 00835 C Regardless of the method type, subroutine JAC must not 00836 C alter T, Y(*), YPRIME(*), H, CJ, or REWT(*). 00837 C You must declare the name JAC in an EXTERNAL statement in 00838 C your program that calls DDASPK. 00839 C 00840 C PSOL -- This is the name of a routine you must supply if you have 00841 C selected a Krylov method (INFO(12) = 1) with preconditioning. 00842 C In the direct case (INFO(12) = 0), PSOL can be absent 00843 C (a dummy routine may have to be supplied to satisfy the 00844 C loader). Otherwise, you must provide a PSOL routine to 00845 C solve linear systems arising from preconditioning. 00846 C When supplied with INFO(12) = 1, the PSOL routine is to 00847 C have the form 00848 C 00849 C SUBROUTINE PSOL (NEQ, T, Y, YPRIME, SAVR, WK, CJ, WGHT, 00850 C WP, IWP, B, EPLIN, IER, RPAR, IPAR) 00851 C 00852 C The PSOL routine must solve linear systems of the form 00853 C P*x = b where P is the left preconditioner matrix. 00854 C 00855 C The right-hand side vector b is in the B array on input, and 00856 C PSOL must return the solution vector x in B. 00857 C The Y, YPRIME, and SAVR arrays contain the current values 00858 C of Y, YPRIME, and the residual G, respectively. 00859 C 00860 C Work space required by JAC and/or PSOL, and space for data to 00861 C be communicated from JAC to PSOL is made available in the form 00862 C of arrays WP and IWP, which are parts of the RWORK and IWORK 00863 C arrays, respectively. The lengths of these real and integer 00864 C work spaces WP and IWP must be supplied in LENWP and LENIWP, 00865 C respectively, as follows.. 00866 C IWORK(27) = LENWP = length of real work space WP 00867 C IWORK(28) = LENIWP = length of integer work space IWP. 00868 C 00869 C WK is a work array of length NEQ for use by PSOL. 00870 C CJ is a scalar, input to PSOL, that is normally proportional 00871 C to 1/H (H = stepsize). If the old value of CJ 00872 C (at the time of the last JAC call) is needed, it must have 00873 C been saved by JAC in WP. 00874 C 00875 C WGHT is an array of weights, to be used if PSOL uses an 00876 C iterative method and performs a convergence test. (In terms 00877 C of the argument REWT to JAC, WGHT is REWT/sqrt(NEQ).) 00878 C If PSOL uses an iterative method, it should use EPLIN 00879 C (a heuristic parameter) as the bound on the weighted norm of 00880 C the residual for the computed solution. Specifically, the 00881 C residual vector R should satisfy 00882 C SQRT (SUM ( (R(i)*WGHT(i))**2 ) ) .le. EPLIN 00883 C 00884 C PSOL must not alter NEQ, T, Y, YPRIME, SAVR, CJ, WGHT, EPLIN. 00885 C 00886 C On return, PSOL should set the error flag IER as follows.. 00887 C IER = 0 if PSOL was successful, 00888 C IER .lt. 0 if an unrecoverable error occurred, meaning 00889 C control will be passed to the calling routine, 00890 C IER .gt. 0 if a recoverable error occurred, meaning that 00891 C the step will be retried with the same step size 00892 C but with a call to JAC to update necessary data, 00893 C unless the Jacobian data is current, in which case 00894 C the step will be retried with a smaller step size. 00895 C IER = 0 on entry to PSOL so need be reset only on a failure. 00896 C 00897 C You must declare the name PSOL in an EXTERNAL statement in 00898 C your program that calls DDASPK. 00899 C 00900 C 00901 C OPTIONALLY REPLACEABLE SUBROUTINE: 00902 C 00903 C DDASPK uses a weighted root-mean-square norm to measure the 00904 C size of various error vectors. The weights used in this norm 00905 C are set in the following subroutine: 00906 C 00907 C SUBROUTINE DDAWTS (NEQ, IWT, RTOL, ATOL, Y, EWT, RPAR, IPAR) 00908 C DIMENSION RTOL(*), ATOL(*), Y(*), EWT(*), RPAR(*), IPAR(*) 00909 C 00910 C A DDAWTS routine has been included with DDASPK which sets the 00911 C weights according to 00912 C EWT(I) = RTOL*ABS(Y(I)) + ATOL 00913 C in the case of scalar tolerances (IWT = 0) or 00914 C EWT(I) = RTOL(I)*ABS(Y(I)) + ATOL(I) 00915 C in the case of array tolerances (IWT = 1). (IWT is INFO(2).) 00916 C In some special cases, it may be appropriate for you to define 00917 C your own error weights by writing a subroutine DDAWTS to be 00918 C called instead of the version supplied. However, this should 00919 C be attempted only after careful thought and consideration. 00920 C If you supply this routine, you may use the tolerances and Y 00921 C as appropriate, but do not overwrite these variables. You 00922 C may also use RPAR and IPAR to communicate data as appropriate. 00923 C ***Note: Aside from the values of the weights, the choice of 00924 C norm used in DDASPK (weighted root-mean-square) is not subject 00925 C to replacement by the user. In this respect, DDASPK is not 00926 C downward-compatible with the original DDASSL solver (in which 00927 C the norm routine was optionally user-replaceable). 00928 C 00929 C 00930 C------OUTPUT - AFTER ANY RETURN FROM DDASPK---------------------------- 00931 C 00932 C The principal aim of the code is to return a computed solution at 00933 C T = TOUT, although it is also possible to obtain intermediate 00934 C results along the way. To find out whether the code achieved its 00935 C goal or if the integration process was interrupted before the task 00936 C was completed, you must check the IDID parameter. 00937 C 00938 C 00939 C T -- The output value of T is the point to which the solution 00940 C was successfully advanced. 00941 C 00942 C Y(*) -- contains the computed solution approximation at T. 00943 C 00944 C YPRIME(*) -- contains the computed derivative approximation at T. 00945 C 00946 C IDID -- reports what the code did, described as follows: 00947 C 00948 C *** TASK COMPLETED *** 00949 C Reported by positive values of IDID 00950 C 00951 C IDID = 1 -- a step was successfully taken in the 00952 C intermediate-output mode. The code has not 00953 C yet reached TOUT. 00954 C 00955 C IDID = 2 -- the integration to TSTOP was successfully 00956 C completed (T = TSTOP) by stepping exactly to TSTOP. 00957 C 00958 C IDID = 3 -- the integration to TOUT was successfully 00959 C completed (T = TOUT) by stepping past TOUT. 00960 C Y(*) and YPRIME(*) are obtained by interpolation. 00961 C 00962 C IDID = 4 -- the initial condition calculation, with 00963 C INFO(11) > 0, was successful, and INFO(14) = 1. 00964 C No integration steps were taken, and the solution 00965 C is not considered to have been started. 00966 C 00967 C *** TASK INTERRUPTED *** 00968 C Reported by negative values of IDID 00969 C 00970 C IDID = -1 -- a large amount of work has been expended 00971 C (about 500 steps). 00972 C 00973 C IDID = -2 -- the error tolerances are too stringent. 00974 C 00975 C IDID = -3 -- the local error test cannot be satisfied 00976 C because you specified a zero component in ATOL 00977 C and the corresponding computed solution component 00978 C is zero. Thus, a pure relative error test is 00979 C impossible for this component. 00980 C 00981 C IDID = -5 -- there were repeated failures in the evaluation 00982 C or processing of the preconditioner (in JAC). 00983 C 00984 C IDID = -6 -- DDASPK had repeated error test failures on the 00985 C last attempted step. 00986 C 00987 C IDID = -7 -- the nonlinear system solver in the time integration 00988 C could not converge. 00989 C 00990 C IDID = -8 -- the matrix of partial derivatives appears 00991 C to be singular (direct method). 00992 C 00993 C IDID = -9 -- the nonlinear system solver in the time integration 00994 C failed to achieve convergence, and there were repeated 00995 C error test failures in this step. 00996 C 00997 C IDID =-10 -- the nonlinear system solver in the time integration 00998 C failed to achieve convergence because IRES was equal 00999 C to -1. 01000 C 01001 C IDID =-11 -- IRES = -2 was encountered and control is 01002 C being returned to the calling program. 01003 C 01004 C IDID =-12 -- DDASPK failed to compute the initial Y, YPRIME. 01005 C 01006 C IDID =-13 -- unrecoverable error encountered inside user's 01007 C PSOL routine, and control is being returned to 01008 C the calling program. 01009 C 01010 C IDID =-14 -- the Krylov linear system solver could not 01011 C achieve convergence. 01012 C 01013 C IDID =-15,..,-32 -- Not applicable for this code. 01014 C 01015 C *** TASK TERMINATED *** 01016 C reported by the value of IDID=-33 01017 C 01018 C IDID = -33 -- the code has encountered trouble from which 01019 C it cannot recover. A message is printed 01020 C explaining the trouble and control is returned 01021 C to the calling program. For example, this occurs 01022 C when invalid input is detected. 01023 C 01024 C RTOL, ATOL -- these quantities remain unchanged except when 01025 C IDID = -2. In this case, the error tolerances have been 01026 C increased by the code to values which are estimated to 01027 C be appropriate for continuing the integration. However, 01028 C the reported solution at T was obtained using the input 01029 C values of RTOL and ATOL. 01030 C 01031 C RWORK, IWORK -- contain information which is usually of no interest 01032 C to the user but necessary for subsequent calls. 01033 C However, you may be interested in the performance data 01034 C listed below. These quantities are accessed in RWORK 01035 C and IWORK but have internal mnemonic names, as follows.. 01036 C 01037 C RWORK(3)--contains H, the step size h to be attempted 01038 C on the next step. 01039 C 01040 C RWORK(4)--contains TN, the current value of the 01041 C independent variable, i.e. the farthest point 01042 C integration has reached. This will differ 01043 C from T if interpolation has been performed 01044 C (IDID = 3). 01045 C 01046 C RWORK(7)--contains HOLD, the stepsize used on the last 01047 C successful step. If INFO(11) = INFO(14) = 1, 01048 C this contains the value of H used in the 01049 C initial condition calculation. 01050 C 01051 C IWORK(7)--contains K, the order of the method to be 01052 C attempted on the next step. 01053 C 01054 C IWORK(8)--contains KOLD, the order of the method used 01055 C on the last step. 01056 C 01057 C IWORK(11)--contains NST, the number of steps (in T) 01058 C taken so far. 01059 C 01060 C IWORK(12)--contains NRE, the number of calls to RES 01061 C so far. 01062 C 01063 C IWORK(13)--contains NJE, the number of calls to JAC so 01064 C far (Jacobian or preconditioner evaluations). 01065 C 01066 C IWORK(14)--contains NETF, the total number of error test 01067 C failures so far. 01068 C 01069 C IWORK(15)--contains NCFN, the total number of nonlinear 01070 C convergence failures so far (includes counts 01071 C of singular iteration matrix or singular 01072 C preconditioners). 01073 C 01074 C IWORK(16)--contains NCFL, the number of convergence 01075 C failures of the linear iteration so far. 01076 C 01077 C IWORK(17)--contains LENIW, the length of IWORK actually 01078 C required. This is defined on normal returns 01079 C and on an illegal input return for 01080 C insufficient storage. 01081 C 01082 C IWORK(18)--contains LENRW, the length of RWORK actually 01083 C required. This is defined on normal returns 01084 C and on an illegal input return for 01085 C insufficient storage. 01086 C 01087 C IWORK(19)--contains NNI, the total number of nonlinear 01088 C iterations so far (each of which calls a 01089 C linear solver). 01090 C 01091 C IWORK(20)--contains NLI, the total number of linear 01092 C (Krylov) iterations so far. 01093 C 01094 C IWORK(21)--contains NPS, the number of PSOL calls so 01095 C far, for preconditioning solve operations or 01096 C for solutions with the user-supplied method. 01097 C 01098 C Note: The various counters in IWORK do not include 01099 C counts during a call made with INFO(11) > 0 and 01100 C INFO(14) = 1. 01101 C 01102 C 01103 C------INPUT - WHAT TO DO TO CONTINUE THE INTEGRATION ----------------- 01104 C (CALLS AFTER THE FIRST) 01105 C 01106 C This code is organized so that subsequent calls to continue the 01107 C integration involve little (if any) additional effort on your 01108 C part. You must monitor the IDID parameter in order to determine 01109 C what to do next. 01110 C 01111 C Recalling that the principal task of the code is to integrate 01112 C from T to TOUT (the interval mode), usually all you will need 01113 C to do is specify a new TOUT upon reaching the current TOUT. 01114 C 01115 C Do not alter any quantity not specifically permitted below. In 01116 C particular do not alter NEQ, T, Y(*), YPRIME(*), RWORK(*), 01117 C IWORK(*), or the differential equation in subroutine RES. Any 01118 C such alteration constitutes a new problem and must be treated 01119 C as such, i.e. you must start afresh. 01120 C 01121 C You cannot change from array to scalar error control or vice 01122 C versa (INFO(2)), but you can change the size of the entries of 01123 C RTOL or ATOL. Increasing a tolerance makes the equation easier 01124 C to integrate. Decreasing a tolerance will make the equation 01125 C harder to integrate and should generally be avoided. 01126 C 01127 C You can switch from the intermediate-output mode to the 01128 C interval mode (INFO(3)) or vice versa at any time. 01129 C 01130 C If it has been necessary to prevent the integration from going 01131 C past a point TSTOP (INFO(4), RWORK(1)), keep in mind that the 01132 C code will not integrate to any TOUT beyond the currently 01133 C specified TSTOP. Once TSTOP has been reached, you must change 01134 C the value of TSTOP or set INFO(4) = 0. You may change INFO(4) 01135 C or TSTOP at any time but you must supply the value of TSTOP in 01136 C RWORK(1) whenever you set INFO(4) = 1. 01137 C 01138 C Do not change INFO(5), INFO(6), INFO(12-17) or their associated 01139 C IWORK/RWORK locations unless you are going to restart the code. 01140 C 01141 C *** FOLLOWING A COMPLETED TASK *** 01142 C 01143 C If.. 01144 C IDID = 1, call the code again to continue the integration 01145 C another step in the direction of TOUT. 01146 C 01147 C IDID = 2 or 3, define a new TOUT and call the code again. 01148 C TOUT must be different from T. You cannot change 01149 C the direction of integration without restarting. 01150 C 01151 C IDID = 4, reset INFO(11) = 0 and call the code again to begin 01152 C the integration. (If you leave INFO(11) > 0 and 01153 C INFO(14) = 1, you may generate an infinite loop.) 01154 C In this situation, the next call to DASPK is 01155 C considered to be the first call for the problem, 01156 C in that all initializations are done. 01157 C 01158 C *** FOLLOWING AN INTERRUPTED TASK *** 01159 C 01160 C To show the code that you realize the task was interrupted and 01161 C that you want to continue, you must take appropriate action and 01162 C set INFO(1) = 1. 01163 C 01164 C If.. 01165 C IDID = -1, the code has taken about 500 steps. If you want to 01166 C continue, set INFO(1) = 1 and call the code again. 01167 C An additional 500 steps will be allowed. 01168 C 01169 C 01170 C IDID = -2, the error tolerances RTOL, ATOL have been increased 01171 C to values the code estimates appropriate for 01172 C continuing. You may want to change them yourself. 01173 C If you are sure you want to continue with relaxed 01174 C error tolerances, set INFO(1) = 1 and call the code 01175 C again. 01176 C 01177 C IDID = -3, a solution component is zero and you set the 01178 C corresponding component of ATOL to zero. If you 01179 C are sure you want to continue, you must first alter 01180 C the error criterion to use positive values of ATOL 01181 C for those components corresponding to zero solution 01182 C components, then set INFO(1) = 1 and call the code 01183 C again. 01184 C 01185 C IDID = -4 --- cannot occur with this code. 01186 C 01187 C IDID = -5, your JAC routine failed with the Krylov method. Check 01188 C for errors in JAC and restart the integration. 01189 C 01190 C IDID = -6, repeated error test failures occurred on the last 01191 C attempted step in DDASPK. A singularity in the 01192 C solution may be present. If you are absolutely 01193 C certain you want to continue, you should restart 01194 C the integration. (Provide initial values of Y and 01195 C YPRIME which are consistent.) 01196 C 01197 C IDID = -7, repeated convergence test failures occurred on the last 01198 C attempted step in DDASPK. An inaccurate or ill- 01199 C conditioned Jacobian or preconditioner may be the 01200 C problem. If you are absolutely certain you want 01201 C to continue, you should restart the integration. 01202 C 01203 C 01204 C IDID = -8, the matrix of partial derivatives is singular, with 01205 C the use of direct methods. Some of your equations 01206 C may be redundant. DDASPK cannot solve the problem 01207 C as stated. It is possible that the redundant 01208 C equations could be removed, and then DDASPK could 01209 C solve the problem. It is also possible that a 01210 C solution to your problem either does not exist 01211 C or is not unique. 01212 C 01213 C IDID = -9, DDASPK had multiple convergence test failures, preceded 01214 C by multiple error test failures, on the last 01215 C attempted step. It is possible that your problem is 01216 C ill-posed and cannot be solved using this code. Or, 01217 C there may be a discontinuity or a singularity in the 01218 C solution. If you are absolutely certain you want to 01219 C continue, you should restart the integration. 01220 C 01221 C IDID = -10, DDASPK had multiple convergence test failures 01222 C because IRES was equal to -1. If you are 01223 C absolutely certain you want to continue, you 01224 C should restart the integration. 01225 C 01226 C IDID = -11, there was an unrecoverable error (IRES = -2) from RES 01227 C inside the nonlinear system solver. Determine the 01228 C cause before trying again. 01229 C 01230 C IDID = -12, DDASPK failed to compute the initial Y and YPRIME 01231 C vectors. This could happen because the initial 01232 C approximation to Y or YPRIME was not very good, or 01233 C because no consistent values of these vectors exist. 01234 C The problem could also be caused by an inaccurate or 01235 C singular iteration matrix, or a poor preconditioner. 01236 C 01237 C IDID = -13, there was an unrecoverable error encountered inside 01238 C your PSOL routine. Determine the cause before 01239 C trying again. 01240 C 01241 C IDID = -14, the Krylov linear system solver failed to achieve 01242 C convergence. This may be due to ill-conditioning 01243 C in the iteration matrix, or a singularity in the 01244 C preconditioner (if one is being used). 01245 C Another possibility is that there is a better 01246 C choice of Krylov parameters (see INFO(13)). 01247 C Possibly the failure is caused by redundant equations 01248 C in the system, or by inconsistent equations. 01249 C In that case, reformulate the system to make it 01250 C consistent and non-redundant. 01251 C 01252 C IDID = -15,..,-32 --- Cannot occur with this code. 01253 C 01254 C *** FOLLOWING A TERMINATED TASK *** 01255 C 01256 C If IDID = -33, you cannot continue the solution of this problem. 01257 C An attempt to do so will result in your run being 01258 C terminated. 01259 C 01260 C --------------------------------------------------------------------- 01261 C 01262 C***REFERENCES 01263 C 1. L. R. Petzold, A Description of DASSL: A Differential/Algebraic 01264 C System Solver, in Scientific Computing, R. S. Stepleman et al. 01265 C (Eds.), North-Holland, Amsterdam, 1983, pp. 65-68. 01266 C 2. K. E. Brenan, S. L. Campbell, and L. R. Petzold, Numerical 01267 C Solution of Initial-Value Problems in Differential-Algebraic 01268 C Equations, Elsevier, New York, 1989. 01269 C 3. P. N. Brown and A. C. Hindmarsh, Reduced Storage Matrix Methods 01270 C in Stiff ODE Systems, J. Applied Mathematics and Computation, 01271 C 31 (1989), pp. 40-91. 01272 C 4. P. N. Brown, A. C. Hindmarsh, and L. R. Petzold, Using Krylov 01273 C Methods in the Solution of Large-Scale Differential-Algebraic 01274 C Systems, SIAM J. Sci. Comp., 15 (1994), pp. 1467-1488. 01275 C 5. P. N. Brown, A. C. Hindmarsh, and L. R. Petzold, Consistent 01276 C Initial Condition Calculation for Differential-Algebraic 01277 C Systems, LLNL Report UCRL-JC-122175, August 1995; submitted to 01278 C SIAM J. Sci. Comp. 01279 C 01280 C***ROUTINES CALLED 01281 C 01282 C The following are all the subordinate routines used by DDASPK. 01283 C 01284 C DDASIC computes consistent initial conditions. 01285 C DYYPNW updates Y and YPRIME in linesearch for initial condition 01286 C calculation. 01287 C DDSTP carries out one step of the integration. 01288 C DCNSTR/DCNST0 check the current solution for constraint violations. 01289 C DDAWTS sets error weight quantities. 01290 C DINVWT tests and inverts the error weights. 01291 C DDATRP performs interpolation to get an output solution. 01292 C DDWNRM computes the weighted root-mean-square norm of a vector. 01293 C D1MACH provides the unit roundoff of the computer. 01294 C XERRWD/XSETF/XSETUN/IXSAV is a package to handle error messages. 01295 C DDASID nonlinear equation driver to initialize Y and YPRIME using 01296 C direct linear system solver methods. Interfaces to Newton 01297 C solver (direct case). 01298 C DNSID solves the nonlinear system for unknown initial values by 01299 C modified Newton iteration and direct linear system methods. 01300 C DLINSD carries out linesearch algorithm for initial condition 01301 C calculation (direct case). 01302 C DFNRMD calculates weighted norm of preconditioned residual in 01303 C initial condition calculation (direct case). 01304 C DNEDD nonlinear equation driver for direct linear system solver 01305 C methods. Interfaces to Newton solver (direct case). 01306 C DMATD assembles the iteration matrix (direct case). 01307 C DNSD solves the associated nonlinear system by modified 01308 C Newton iteration and direct linear system methods. 01309 C DSLVD interfaces to linear system solver (direct case). 01310 C DDASIK nonlinear equation driver to initialize Y and YPRIME using 01311 C Krylov iterative linear system methods. Interfaces to 01312 C Newton solver (Krylov case). 01313 C DNSIK solves the nonlinear system for unknown initial values by 01314 C Newton iteration and Krylov iterative linear system methods. 01315 C DLINSK carries out linesearch algorithm for initial condition 01316 C calculation (Krylov case). 01317 C DFNRMK calculates weighted norm of preconditioned residual in 01318 C initial condition calculation (Krylov case). 01319 C DNEDK nonlinear equation driver for iterative linear system solver 01320 C methods. Interfaces to Newton solver (Krylov case). 01321 C DNSK solves the associated nonlinear system by Inexact Newton 01322 C iteration and (linear) Krylov iteration. 01323 C DSLVK interfaces to linear system solver (Krylov case). 01324 C DSPIGM solves a linear system by SPIGMR algorithm. 01325 C DATV computes matrix-vector product in Krylov algorithm. 01326 C DORTH performs orthogonalization of Krylov basis vectors. 01327 C DHEQR performs QR factorization of Hessenberg matrix. 01328 C DHELS finds least-squares solution of Hessenberg linear system. 01329 C DGETRF, DGETRS, DGBTRF, DGBTRS are LAPACK routines for solving 01330 C linear systems (dense or band direct methods). 01331 C DAXPY, DCOPY, DDOT, DNRM2, DSCAL are Basic Linear Algebra (BLAS) 01332 C routines. 01333 C 01334 C The routines called directly by DDASPK are: 01335 C DCNST0, DDAWTS, DINVWT, D1MACH, DDWNRM, DDASIC, DDATRP, DDSTP, 01336 C XERRWD 01337 C 01338 C***END PROLOGUE DDASPK 01339 C 01340 C 01341 IMPLICIT DOUBLE PRECISION(A-H,O-Z) 01342 LOGICAL DONE, LAVL, LCFN, LCFL, LWARN 01343 DIMENSION Y(*),YPRIME(*) 01344 DIMENSION INFO(20) 01345 DIMENSION RWORK(LRW),IWORK(LIW) 01346 DIMENSION RTOL(*),ATOL(*) 01347 DIMENSION RPAR(*),IPAR(*) 01348 CHARACTER MSG*80 01349 EXTERNAL RES, JAC, PSOL, DDASID, DDASIK, DNEDD, DNEDK 01350 C 01351 C Set pointers into IWORK. 01352 C 01353 PARAMETER (LML=1, LMU=2, LMTYPE=4, 01354 * LIWM=1, LMXORD=3, LJCALC=5, LPHASE=6, LK=7, LKOLD=8, 01355 * LNS=9, LNSTL=10, LNST=11, LNRE=12, LNJE=13, LETF=14, LNCFN=15, 01356 * LNCFL=16, LNIW=17, LNRW=18, LNNI=19, LNLI=20, LNPS=21, 01357 * LNPD=22, LMITER=23, LMAXL=24, LKMP=25, LNRMAX=26, LLNWP=27, 01358 * LLNIWP=28, LLOCWP=29, LLCIWP=30, LKPRIN=31, 01359 * LMXNIT=32, LMXNJ=33, LMXNH=34, LLSOFF=35, LICNS=41) 01360 C 01361 C Set pointers into RWORK. 01362 C 01363 PARAMETER (LTSTOP=1, LHMAX=2, LH=3, LTN=4, LCJ=5, LCJOLD=6, 01364 * LHOLD=7, LS=8, LROUND=9, LEPLI=10, LSQRN=11, LRSQRN=12, 01365 * LEPCON=13, LSTOL=14, LEPIN=15, 01366 * LALPHA=21, LBETA=27, LGAMMA=33, LPSI=39, LSIGMA=45, LDELTA=51) 01367 C 01368 SAVE LID, LENID, NONNEG 01369 C 01370 C 01371 C***FIRST EXECUTABLE STATEMENT DDASPK 01372 C 01373 C 01374 IF(INFO(1).NE.0) GO TO 100 01375 C 01376 C----------------------------------------------------------------------- 01377 C This block is executed for the initial call only. 01378 C It contains checking of inputs and initializations. 01379 C----------------------------------------------------------------------- 01380 C 01381 C First check INFO array to make sure all elements of INFO 01382 C Are within the proper range. (INFO(1) is checked later, because 01383 C it must be tested on every call.) ITEMP holds the location 01384 C within INFO which may be out of range. 01385 C 01386 DO 10 I=2,9 01387 ITEMP = I 01388 IF (INFO(I) .NE. 0 .AND. INFO(I) .NE. 1) GO TO 701 01389 10 CONTINUE 01390 ITEMP = 10 01391 IF(INFO(10).LT.0 .OR. INFO(10).GT.3) GO TO 701 01392 ITEMP = 11 01393 IF(INFO(11).LT.0 .OR. INFO(11).GT.2) GO TO 701 01394 DO 15 I=12,17 01395 ITEMP = I 01396 IF (INFO(I) .NE. 0 .AND. INFO(I) .NE. 1) GO TO 701 01397 15 CONTINUE 01398 ITEMP = 18 01399 IF(INFO(18).LT.0 .OR. INFO(18).GT.2) GO TO 701 01400 01401 C 01402 C Check NEQ to see if it is positive. 01403 C 01404 IF (NEQ .LE. 0) GO TO 702 01405 C 01406 C Check and compute maximum order. 01407 C 01408 MXORD=5 01409 IF (INFO(9) .NE. 0) THEN 01410 MXORD=IWORK(LMXORD) 01411 IF (MXORD .LT. 1 .OR. MXORD .GT. 5) GO TO 703 01412 ENDIF 01413 IWORK(LMXORD)=MXORD 01414 C 01415 C Set and/or check inputs for constraint checking (INFO(10) .NE. 0). 01416 C Set values for ICNFLG, NONNEG, and pointer LID. 01417 C 01418 ICNFLG = 0 01419 NONNEG = 0 01420 LID = LICNS 01421 IF (INFO(10) .EQ. 0) GO TO 20 01422 IF (INFO(10) .EQ. 1) THEN 01423 ICNFLG = 1 01424 NONNEG = 0 01425 LID = LICNS + NEQ 01426 ELSEIF (INFO(10) .EQ. 2) THEN 01427 ICNFLG = 0 01428 NONNEG = 1 01429 ELSE 01430 ICNFLG = 1 01431 NONNEG = 1 01432 LID = LICNS + NEQ 01433 ENDIF 01434 C 01435 20 CONTINUE 01436 C 01437 C Set and/or check inputs for Krylov solver (INFO(12) .NE. 0). 01438 C If indicated, set default values for MAXL, KMP, NRMAX, and EPLI. 01439 C Otherwise, verify inputs required for iterative solver. 01440 C 01441 IF (INFO(12) .EQ. 0) GO TO 25 01442 C 01443 IWORK(LMITER) = INFO(12) 01444 IF (INFO(13) .EQ. 0) THEN 01445 IWORK(LMAXL) = MIN(5,NEQ) 01446 IWORK(LKMP) = IWORK(LMAXL) 01447 IWORK(LNRMAX) = 5 01448 RWORK(LEPLI) = 0.05D0 01449 ELSE 01450 IF(IWORK(LMAXL) .LT. 1 .OR. IWORK(LMAXL) .GT. NEQ) GO TO 720 01451 IF(IWORK(LKMP) .LT. 1 .OR. IWORK(LKMP) .GT. IWORK(LMAXL)) 01452 1 GO TO 721 01453 IF(IWORK(LNRMAX) .LT. 0) GO TO 722 01454 IF(RWORK(LEPLI).LE.0.0D0 .OR. RWORK(LEPLI).GE.1.0D0)GO TO 723 01455 ENDIF 01456 C 01457 25 CONTINUE 01458 C 01459 C Set and/or check controls for the initial condition calculation 01460 C (INFO(11) .GT. 0). If indicated, set default values. 01461 C Otherwise, verify inputs required for iterative solver. 01462 C 01463 IF (INFO(11) .EQ. 0) GO TO 30 01464 IF (INFO(17) .EQ. 0) THEN 01465 IWORK(LMXNIT) = 5 01466 IF (INFO(12) .GT. 0) IWORK(LMXNIT) = 15 01467 IWORK(LMXNJ) = 6 01468 IF (INFO(12) .GT. 0) IWORK(LMXNJ) = 2 01469 IWORK(LMXNH) = 5 01470 IWORK(LLSOFF) = 0 01471 RWORK(LEPIN) = 0.01D0 01472 ELSE 01473 IF (IWORK(LMXNIT) .LE. 0) GO TO 725 01474 IF (IWORK(LMXNJ) .LE. 0) GO TO 725 01475 IF (IWORK(LMXNH) .LE. 0) GO TO 725 01476 LSOFF = IWORK(LLSOFF) 01477 IF (LSOFF .LT. 0 .OR. LSOFF .GT. 1) GO TO 725 01478 IF (RWORK(LEPIN) .LE. 0.0D0) GO TO 725 01479 ENDIF 01480 C 01481 30 CONTINUE 01482 C 01483 C Below is the computation and checking of the work array lengths 01484 C LENIW and LENRW, using direct methods (INFO(12) = 0) or 01485 C the Krylov methods (INFO(12) = 1). 01486 C 01487 LENIC = 0 01488 IF (INFO(10) .EQ. 1 .OR. INFO(10) .EQ. 3) LENIC = NEQ 01489 LENID = 0 01490 IF (INFO(11) .EQ. 1 .OR. INFO(16) .EQ. 1) LENID = NEQ 01491 IF (INFO(12) .EQ. 0) THEN 01492 C 01493 C Compute MTYPE, etc. Check ML and MU. 01494 C 01495 NCPHI = MAX(MXORD + 1, 4) 01496 IF(INFO(6).EQ.0) THEN 01497 LENPD = NEQ**2 01498 LENRW = 50 + (NCPHI+3)*NEQ + LENPD 01499 IF(INFO(5).EQ.0) THEN 01500 IWORK(LMTYPE)=2 01501 ELSE 01502 IWORK(LMTYPE)=1 01503 ENDIF 01504 ELSE 01505 IF(IWORK(LML).LT.0.OR.IWORK(LML).GE.NEQ)GO TO 717 01506 IF(IWORK(LMU).LT.0.OR.IWORK(LMU).GE.NEQ)GO TO 718 01507 LENPD=(2*IWORK(LML)+IWORK(LMU)+1)*NEQ 01508 IF(INFO(5).EQ.0) THEN 01509 IWORK(LMTYPE)=5 01510 MBAND=IWORK(LML)+IWORK(LMU)+1 01511 MSAVE=(NEQ/MBAND)+1 01512 LENRW = 50 + (NCPHI+3)*NEQ + LENPD + 2*MSAVE 01513 ELSE 01514 IWORK(LMTYPE)=4 01515 LENRW = 50 + (NCPHI+3)*NEQ + LENPD 01516 ENDIF 01517 ENDIF 01518 C 01519 C Compute LENIW, LENWP, LENIWP. 01520 C 01521 LENIW = 40 + LENIC + LENID + NEQ 01522 LENWP = 0 01523 LENIWP = 0 01524 C 01525 ELSE IF (INFO(12) .EQ. 1) THEN 01526 MAXL = IWORK(LMAXL) 01527 LENWP = IWORK(LLNWP) 01528 LENIWP = IWORK(LLNIWP) 01529 LENPD = (MAXL+3+MIN0(1,MAXL-IWORK(LKMP)))*NEQ 01530 1 + (MAXL+3)*MAXL + 1 + LENWP 01531 LENRW = 50 + (IWORK(LMXORD)+5)*NEQ + LENPD 01532 LENIW = 40 + LENIC + LENID + LENIWP 01533 C 01534 ENDIF 01535 IF(INFO(16) .NE. 0) LENRW = LENRW + NEQ 01536 C 01537 C Check lengths of RWORK and IWORK. 01538 C 01539 IWORK(LNIW)=LENIW 01540 IWORK(LNRW)=LENRW 01541 IWORK(LNPD)=LENPD 01542 IWORK(LLOCWP) = LENPD-LENWP+1 01543 IF(LRW.LT.LENRW)GO TO 704 01544 IF(LIW.LT.LENIW)GO TO 705 01545 C 01546 C Check ICNSTR for legality. 01547 C 01548 IF (LENIC .GT. 0) THEN 01549 DO 40 I = 1,NEQ 01550 ICI = IWORK(LICNS-1+I) 01551 IF (ICI .LT. -2 .OR. ICI .GT. 2) GO TO 726 01552 40 CONTINUE 01553 ENDIF 01554 C 01555 C Check Y for consistency with constraints. 01556 C 01557 IF (LENIC .GT. 0) THEN 01558 CALL DCNST0(NEQ,Y,IWORK(LICNS),IRET) 01559 IF (IRET .NE. 0) GO TO 727 01560 ENDIF 01561 C 01562 C Check ID for legality. 01563 C 01564 IF (LENID .GT. 0) THEN 01565 DO 50 I = 1,NEQ 01566 IDI = IWORK(LID-1+I) 01567 IF (IDI .NE. 1 .AND. IDI .NE. -1) GO TO 724 01568 50 CONTINUE 01569 ENDIF 01570 C 01571 C Check to see that TOUT is different from T. 01572 C 01573 IF(TOUT .EQ. T)GO TO 719 01574 C 01575 C Check HMAX. 01576 C 01577 IF(INFO(7) .NE. 0) THEN 01578 HMAX = RWORK(LHMAX) 01579 IF (HMAX .LE. 0.0D0) GO TO 710 01580 ENDIF 01581 C 01582 C Initialize counters and other flags. 01583 C 01584 IWORK(LNST)=0 01585 IWORK(LNRE)=0 01586 IWORK(LNJE)=0 01587 IWORK(LETF)=0 01588 IWORK(LNCFN)=0 01589 IWORK(LNNI)=0 01590 IWORK(LNLI)=0 01591 IWORK(LNPS)=0 01592 IWORK(LNCFL)=0 01593 IWORK(LKPRIN)=INFO(18) 01594 IDID=1 01595 GO TO 200 01596 C 01597 C----------------------------------------------------------------------- 01598 C This block is for continuation calls only. 01599 C Here we check INFO(1), and if the last step was interrupted, 01600 C we check whether appropriate action was taken. 01601 C----------------------------------------------------------------------- 01602 C 01603 100 CONTINUE 01604 IF(INFO(1).EQ.1)GO TO 110 01605 ITEMP = 1 01606 IF(INFO(1).NE.-1)GO TO 701 01607 C 01608 C If we are here, the last step was interrupted by an error 01609 C condition from DDSTP, and appropriate action was not taken. 01610 C This is a fatal error. 01611 C 01612 MSG = 'DASPK-- THE LAST STEP TERMINATED WITH A NEGATIVE' 01613 CALL XERRWD(MSG,49,201,0,0,0,0,0,0.0D0,0.0D0) 01614 MSG = 'DASPK-- VALUE (=I1) OF IDID AND NO APPROPRIATE' 01615 CALL XERRWD(MSG,47,202,0,1,IDID,0,0,0.0D0,0.0D0) 01616 MSG = 'DASPK-- ACTION WAS TAKEN. RUN TERMINATED' 01617 CALL XERRWD(MSG,41,203,1,0,0,0,0,0.0D0,0.0D0) 01618 RETURN 01619 110 CONTINUE 01620 C 01621 C----------------------------------------------------------------------- 01622 C This block is executed on all calls. 01623 C 01624 C Counters are saved for later checks of performance. 01625 C Then the error tolerance parameters are checked, and the 01626 C work array pointers are set. 01627 C----------------------------------------------------------------------- 01628 C 01629 200 CONTINUE 01630 C 01631 C Save counters for use later. 01632 C 01633 IWORK(LNSTL)=IWORK(LNST) 01634 NLI0 = IWORK(LNLI) 01635 NNI0 = IWORK(LNNI) 01636 NCFN0 = IWORK(LNCFN) 01637 NCFL0 = IWORK(LNCFL) 01638 NWARN = 0 01639 C 01640 C Check RTOL and ATOL. 01641 C 01642 NZFLG = 0 01643 RTOLI = RTOL(1) 01644 ATOLI = ATOL(1) 01645 DO 210 I=1,NEQ 01646 IF (INFO(2) .EQ. 1) RTOLI = RTOL(I) 01647 IF (INFO(2) .EQ. 1) ATOLI = ATOL(I) 01648 IF (RTOLI .GT. 0.0D0 .OR. ATOLI .GT. 0.0D0) NZFLG = 1 01649 IF (RTOLI .LT. 0.0D0) GO TO 706 01650 IF (ATOLI .LT. 0.0D0) GO TO 707 01651 210 CONTINUE 01652 IF (NZFLG .EQ. 0) GO TO 708 01653 C 01654 C Set pointers to RWORK and IWORK segments. 01655 C For direct methods, SAVR is not used. 01656 C 01657 IWORK(LLCIWP) = LID + LENID 01658 LSAVR = LDELTA 01659 IF (INFO(12) .NE. 0) LSAVR = LDELTA + NEQ 01660 LE = LSAVR + NEQ 01661 LWT = LE + NEQ 01662 LVT = LWT 01663 IF (INFO(16) .NE. 0) LVT = LWT + NEQ 01664 LPHI = LVT + NEQ 01665 LWM = LPHI + (IWORK(LMXORD)+1)*NEQ 01666 IF (INFO(1) .EQ. 1) GO TO 400 01667 C 01668 C----------------------------------------------------------------------- 01669 C This block is executed on the initial call only. 01670 C Set the initial step size, the error weight vector, and PHI. 01671 C Compute unknown initial components of Y and YPRIME, if requested. 01672 C----------------------------------------------------------------------- 01673 C 01674 300 CONTINUE 01675 TN=T 01676 IDID=1 01677 C 01678 C Set error weight array WT and altered weight array VT. 01679 C 01680 CALL DDAWTS(NEQ,INFO(2),RTOL,ATOL,Y,RWORK(LWT),RPAR,IPAR) 01681 CALL DINVWT(NEQ,RWORK(LWT),IER) 01682 IF (IER .NE. 0) GO TO 713 01683 IF (INFO(16) .NE. 0) THEN 01684 DO 305 I = 1, NEQ 01685 305 RWORK(LVT+I-1) = MAX(IWORK(LID+I-1),0)*RWORK(LWT+I-1) 01686 ENDIF 01687 C 01688 C Compute unit roundoff and HMIN. 01689 C 01690 UROUND = D1MACH(4) 01691 RWORK(LROUND) = UROUND 01692 HMIN = 4.0D0*UROUND*MAX(ABS(T),ABS(TOUT)) 01693 C 01694 C Set/check STPTOL control for initial condition calculation. 01695 C 01696 IF (INFO(11) .NE. 0) THEN 01697 IF( INFO(17) .EQ. 0) THEN 01698 RWORK(LSTOL) = UROUND**.6667D0 01699 ELSE 01700 IF (RWORK(LSTOL) .LE. 0.0D0) GO TO 725 01701 ENDIF 01702 ENDIF 01703 C 01704 C Compute EPCON and square root of NEQ and its reciprocal, used 01705 C inside iterative solver. 01706 C 01707 RWORK(LEPCON) = 0.33D0 01708 FLOATN = NEQ 01709 RWORK(LSQRN) = SQRT(FLOATN) 01710 RWORK(LRSQRN) = 1.D0/RWORK(LSQRN) 01711 C 01712 C Check initial interval to see that it is long enough. 01713 C 01714 TDIST = ABS(TOUT - T) 01715 IF(TDIST .LT. HMIN) GO TO 714 01716 C 01717 C Check H0, if this was input. 01718 C 01719 IF (INFO(8) .EQ. 0) GO TO 310 01720 H0 = RWORK(LH) 01721 IF ((TOUT - T)*H0 .LT. 0.0D0) GO TO 711 01722 IF (H0 .EQ. 0.0D0) GO TO 712 01723 GO TO 320 01724 310 CONTINUE 01725 C 01726 C Compute initial stepsize, to be used by either 01727 C DDSTP or DDASIC, depending on INFO(11). 01728 C 01729 H0 = 0.001D0*TDIST 01730 YPNORM = DDWNRM(NEQ,YPRIME,RWORK(LVT),RPAR,IPAR) 01731 IF (YPNORM .GT. 0.5D0/H0) H0 = 0.5D0/YPNORM 01732 H0 = SIGN(H0,TOUT-T) 01733 C 01734 C Adjust H0 if necessary to meet HMAX bound. 01735 C 01736 320 IF (INFO(7) .EQ. 0) GO TO 330 01737 RH = ABS(H0)/RWORK(LHMAX) 01738 IF (RH .GT. 1.0D0) H0 = H0/RH 01739 C 01740 C Check against TSTOP, if applicable. 01741 C 01742 330 IF (INFO(4) .EQ. 0) GO TO 340 01743 TSTOP = RWORK(LTSTOP) 01744 IF ((TSTOP - T)*H0 .LT. 0.0D0) GO TO 715 01745 IF ((T + H0 - TSTOP)*H0 .GT. 0.0D0) H0 = TSTOP - T 01746 IF ((TSTOP - TOUT)*H0 .LT. 0.0D0) GO TO 709 01747 C 01748 340 IF (INFO(11) .EQ. 0) GO TO 370 01749 C 01750 C Compute unknown components of initial Y and YPRIME, depending 01751 C on INFO(11) and INFO(12). INFO(12) represents the nonlinear 01752 C solver type (direct/Krylov). Pass the name of the specific 01753 C nonlinear solver, depending on INFO(12). The location of the work 01754 C arrays SAVR, YIC, YPIC, PWK also differ in the two cases. 01755 C 01756 NWT = 1 01757 EPCONI = RWORK(LEPIN)*RWORK(LEPCON) 01758 350 IF (INFO(12) .EQ. 0) THEN 01759 LYIC = LPHI + 2*NEQ 01760 LYPIC = LYIC + NEQ 01761 LPWK = LYPIC 01762 CALL DDASIC(TN,Y,YPRIME,NEQ,INFO(11),IWORK(LID), 01763 * RES,JAC,PSOL,H0,RWORK(LWT),NWT,IDID,RPAR,IPAR, 01764 * RWORK(LPHI),RWORK(LSAVR),RWORK(LDELTA),RWORK(LE), 01765 * RWORK(LYIC),RWORK(LYPIC),RWORK(LPWK),RWORK(LWM),IWORK(LIWM), 01766 * HMIN,RWORK(LROUND),RWORK(LEPLI),RWORK(LSQRN),RWORK(LRSQRN), 01767 * EPCONI,RWORK(LSTOL),INFO(15),ICNFLG,IWORK(LICNS),DDASID) 01768 ELSE IF (INFO(12) .EQ. 1) THEN 01769 LYIC = LWM 01770 LYPIC = LYIC + NEQ 01771 LPWK = LYPIC + NEQ 01772 CALL DDASIC(TN,Y,YPRIME,NEQ,INFO(11),IWORK(LID), 01773 * RES,JAC,PSOL,H0,RWORK(LWT),NWT,IDID,RPAR,IPAR, 01774 * RWORK(LPHI),RWORK(LSAVR),RWORK(LDELTA),RWORK(LE), 01775 * RWORK(LYIC),RWORK(LYPIC),RWORK(LPWK),RWORK(LWM),IWORK(LIWM), 01776 * HMIN,RWORK(LROUND),RWORK(LEPLI),RWORK(LSQRN),RWORK(LRSQRN), 01777 * EPCONI,RWORK(LSTOL),INFO(15),ICNFLG,IWORK(LICNS),DDASIK) 01778 ENDIF 01779 C 01780 IF (IDID .LT. 0) GO TO 600 01781 C 01782 C DDASIC was successful. If this was the first call to DDASIC, 01783 C update the WT array (with the current Y) and call it again. 01784 C 01785 IF (NWT .EQ. 2) GO TO 355 01786 NWT = 2 01787 CALL DDAWTS(NEQ,INFO(2),RTOL,ATOL,Y,RWORK(LWT),RPAR,IPAR) 01788 CALL DINVWT(NEQ,RWORK(LWT),IER) 01789 IF (IER .NE. 0) GO TO 713 01790 GO TO 350 01791 C 01792 C If INFO(14) = 1, return now with IDID = 4. 01793 C 01794 355 IF (INFO(14) .EQ. 1) THEN 01795 IDID = 4 01796 H = H0 01797 IF (INFO(11) .EQ. 1) RWORK(LHOLD) = H0 01798 GO TO 590 01799 ENDIF 01800 C 01801 C Update the WT and VT arrays one more time, with the new Y. 01802 C 01803 CALL DDAWTS(NEQ,INFO(2),RTOL,ATOL,Y,RWORK(LWT),RPAR,IPAR) 01804 CALL DINVWT(NEQ,RWORK(LWT),IER) 01805 IF (IER .NE. 0) GO TO 713 01806 IF (INFO(16) .NE. 0) THEN 01807 DO 357 I = 1, NEQ 01808 357 RWORK(LVT+I-1) = MAX(IWORK(LID+I-1),0)*RWORK(LWT+I-1) 01809 ENDIF 01810 C 01811 C Reset the initial stepsize to be used by DDSTP. 01812 C Use H0, if this was input. Otherwise, recompute H0, 01813 C and adjust it if necessary to meet HMAX bound. 01814 C 01815 IF (INFO(8) .NE. 0) THEN 01816 H0 = RWORK(LH) 01817 GO TO 360 01818 ENDIF 01819 C 01820 H0 = 0.001D0*TDIST 01821 YPNORM = DDWNRM(NEQ,YPRIME,RWORK(LVT),RPAR,IPAR) 01822 IF (YPNORM .GT. 0.5D0/H0) H0 = 0.5D0/YPNORM 01823 H0 = SIGN(H0,TOUT-T) 01824 C 01825 360 IF (INFO(7) .NE. 0) THEN 01826 RH = ABS(H0)/RWORK(LHMAX) 01827 IF (RH .GT. 1.0D0) H0 = H0/RH 01828 ENDIF 01829 C 01830 C Check against TSTOP, if applicable. 01831 C 01832 IF (INFO(4) .NE. 0) THEN 01833 TSTOP = RWORK(LTSTOP) 01834 IF ((T + H0 - TSTOP)*H0 .GT. 0.0D0) H0 = TSTOP - T 01835 ENDIF 01836 C 01837 C Load H and RWORK(LH) with H0. 01838 C 01839 370 H = H0 01840 RWORK(LH) = H 01841 C 01842 C Load Y and H*YPRIME into PHI(*,1) and PHI(*,2). 01843 C 01844 ITEMP = LPHI + NEQ 01845 DO 380 I = 1,NEQ 01846 RWORK(LPHI + I - 1) = Y(I) 01847 380 RWORK(ITEMP + I - 1) = H*YPRIME(I) 01848 C 01849 GO TO 500 01850 C 01851 C----------------------------------------------------------------------- 01852 C This block is for continuation calls only. 01853 C Its purpose is to check stop conditions before taking a step. 01854 C Adjust H if necessary to meet HMAX bound. 01855 C----------------------------------------------------------------------- 01856 C 01857 400 CONTINUE 01858 UROUND=RWORK(LROUND) 01859 DONE = .FALSE. 01860 TN=RWORK(LTN) 01861 H=RWORK(LH) 01862 IF(INFO(7) .EQ. 0) GO TO 410 01863 RH = ABS(H)/RWORK(LHMAX) 01864 IF(RH .GT. 1.0D0) H = H/RH 01865 410 CONTINUE 01866 IF(T .EQ. TOUT) GO TO 719 01867 IF((T - TOUT)*H .GT. 0.0D0) GO TO 711 01868 IF(INFO(4) .EQ. 1) GO TO 430 01869 IF(INFO(3) .EQ. 1) GO TO 420 01870 IF((TN-TOUT)*H.LT.0.0D0)GO TO 490 01871 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD), 01872 * RWORK(LPHI),RWORK(LPSI)) 01873 T=TOUT 01874 IDID = 3 01875 DONE = .TRUE. 01876 GO TO 490 01877 420 IF((TN-T)*H .LE. 0.0D0) GO TO 490 01878 IF((TN - TOUT)*H .GT. 0.0D0) GO TO 425 01879 CALL DDATRP(TN,TN,Y,YPRIME,NEQ,IWORK(LKOLD), 01880 * RWORK(LPHI),RWORK(LPSI)) 01881 T = TN 01882 IDID = 1 01883 DONE = .TRUE. 01884 GO TO 490 01885 425 CONTINUE 01886 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD), 01887 * RWORK(LPHI),RWORK(LPSI)) 01888 T = TOUT 01889 IDID = 3 01890 DONE = .TRUE. 01891 GO TO 490 01892 430 IF(INFO(3) .EQ. 1) GO TO 440 01893 TSTOP=RWORK(LTSTOP) 01894 IF((TN-TSTOP)*H.GT.0.0D0) GO TO 715 01895 IF((TSTOP-TOUT)*H.LT.0.0D0)GO TO 709 01896 IF((TN-TOUT)*H.LT.0.0D0)GO TO 450 01897 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD), 01898 * RWORK(LPHI),RWORK(LPSI)) 01899 T=TOUT 01900 IDID = 3 01901 DONE = .TRUE. 01902 GO TO 490 01903 440 TSTOP = RWORK(LTSTOP) 01904 IF((TN-TSTOP)*H .GT. 0.0D0) GO TO 715 01905 IF((TSTOP-TOUT)*H .LT. 0.0D0) GO TO 709 01906 IF((TN-T)*H .LE. 0.0D0) GO TO 450 01907 IF((TN - TOUT)*H .GT. 0.0D0) GO TO 445 01908 CALL DDATRP(TN,TN,Y,YPRIME,NEQ,IWORK(LKOLD), 01909 * RWORK(LPHI),RWORK(LPSI)) 01910 T = TN 01911 IDID = 1 01912 DONE = .TRUE. 01913 GO TO 490 01914 445 CONTINUE 01915 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD), 01916 * RWORK(LPHI),RWORK(LPSI)) 01917 T = TOUT 01918 IDID = 3 01919 DONE = .TRUE. 01920 GO TO 490 01921 450 CONTINUE 01922 C 01923 C Check whether we are within roundoff of TSTOP. 01924 C 01925 IF(ABS(TN-TSTOP).GT.100.0D0*UROUND* 01926 * (ABS(TN)+ABS(H)))GO TO 460 01927 CALL DDATRP(TN,TSTOP,Y,YPRIME,NEQ,IWORK(LKOLD), 01928 * RWORK(LPHI),RWORK(LPSI)) 01929 IDID=2 01930 T=TSTOP 01931 DONE = .TRUE. 01932 GO TO 490 01933 460 TNEXT=TN+H 01934 IF((TNEXT-TSTOP)*H.LE.0.0D0)GO TO 490 01935 H=TSTOP-TN 01936 RWORK(LH)=H 01937 C 01938 490 IF (DONE) GO TO 590 01939 C 01940 C----------------------------------------------------------------------- 01941 C The next block contains the call to the one-step integrator DDSTP. 01942 C This is a looping point for the integration steps. 01943 C Check for too many steps. 01944 C Check for poor Newton/Krylov performance. 01945 C Update WT. Check for too much accuracy requested. 01946 C Compute minimum stepsize. 01947 C----------------------------------------------------------------------- 01948 C 01949 500 CONTINUE 01950 C 01951 C Check for too many steps. 01952 C 01953 IF((IWORK(LNST)-IWORK(LNSTL)).LT.500) GO TO 505 01954 IDID=-1 01955 GO TO 527 01956 C 01957 C Check for poor Newton/Krylov performance. 01958 C 01959 505 IF (INFO(12) .EQ. 0) GO TO 510 01960 NSTD = IWORK(LNST) - IWORK(LNSTL) 01961 NNID = IWORK(LNNI) - NNI0 01962 IF (NSTD .LT. 10 .OR. NNID .EQ. 0) GO TO 510 01963 AVLIN = REAL(IWORK(LNLI) - NLI0)/REAL(NNID) 01964 RCFN = REAL(IWORK(LNCFN) - NCFN0)/REAL(NSTD) 01965 RCFL = REAL(IWORK(LNCFL) - NCFL0)/REAL(NNID) 01966 FMAXL = IWORK(LMAXL) 01967 LAVL = AVLIN .GT. FMAXL 01968 LCFN = RCFN .GT. 0.9D0 01969 LCFL = RCFL .GT. 0.9D0 01970 LWARN = LAVL .OR. LCFN .OR. LCFL 01971 IF (.NOT.LWARN) GO TO 510 01972 NWARN = NWARN + 1 01973 IF (NWARN .GT. 10) GO TO 510 01974 IF (LAVL) THEN 01975 MSG = 'DASPK-- Warning. Poor iterative algorithm performance ' 01976 CALL XERRWD (MSG, 56, 501, 0, 0, 0, 0, 0, 0.0D0, 0.0D0) 01977 MSG = ' at T = R1. Average no. of linear iterations = R2 ' 01978 CALL XERRWD (MSG, 56, 501, 0, 0, 0, 0, 2, TN, AVLIN) 01979 ENDIF 01980 IF (LCFN) THEN 01981 MSG = 'DASPK-- Warning. Poor iterative algorithm performance ' 01982 CALL XERRWD (MSG, 56, 502, 0, 0, 0, 0, 0, 0.0D0, 0.0D0) 01983 MSG = ' at T = R1. Nonlinear convergence failure rate = R2' 01984 CALL XERRWD (MSG, 56, 502, 0, 0, 0, 0, 2, TN, RCFN) 01985 ENDIF 01986 IF (LCFL) THEN 01987 MSG = 'DASPK-- Warning. Poor iterative algorithm performance ' 01988 CALL XERRWD (MSG, 56, 503, 0, 0, 0, 0, 0, 0.0D0, 0.0D0) 01989 MSG = ' at T = R1. Linear convergence failure rate = R2 ' 01990 CALL XERRWD (MSG, 56, 503, 0, 0, 0, 0, 2, TN, RCFL) 01991 ENDIF 01992 C 01993 C Update WT and VT, if this is not the first call. 01994 C 01995 510 CALL DDAWTS(NEQ,INFO(2),RTOL,ATOL,RWORK(LPHI),RWORK(LWT), 01996 * RPAR,IPAR) 01997 CALL DINVWT(NEQ,RWORK(LWT),IER) 01998 IF (IER .NE. 0) THEN 01999 IDID = -3 02000 GO TO 527 02001 ENDIF 02002 IF (INFO(16) .NE. 0) THEN 02003 DO 515 I = 1, NEQ 02004 515 RWORK(LVT+I-1) = MAX(IWORK(LID+I-1),0)*RWORK(LWT+I-1) 02005 ENDIF 02006 C 02007 C Test for too much accuracy requested. 02008 C 02009 R = DDWNRM(NEQ,RWORK(LPHI),RWORK(LWT),RPAR,IPAR)*100.0D0*UROUND 02010 IF (R .LE. 1.0D0) GO TO 525 02011 C 02012 C Multiply RTOL and ATOL by R and return. 02013 C 02014 IF(INFO(2).EQ.1)GO TO 523 02015 RTOL(1)=R*RTOL(1) 02016 ATOL(1)=R*ATOL(1) 02017 IDID=-2 02018 GO TO 527 02019 523 DO 524 I=1,NEQ 02020 RTOL(I)=R*RTOL(I) 02021 524 ATOL(I)=R*ATOL(I) 02022 IDID=-2 02023 GO TO 527 02024 525 CONTINUE 02025 C 02026 C Compute minimum stepsize. 02027 C 02028 HMIN=4.0D0*UROUND*MAX(ABS(TN),ABS(TOUT)) 02029 C 02030 C Test H vs. HMAX 02031 IF (INFO(7) .NE. 0) THEN 02032 RH = ABS(H)/RWORK(LHMAX) 02033 IF (RH .GT. 1.0D0) H = H/RH 02034 ENDIF 02035 C 02036 C Call the one-step integrator. 02037 C Note that INFO(12) represents the nonlinear solver type. 02038 C Pass the required nonlinear solver, depending upon INFO(12). 02039 C 02040 IF (INFO(12) .EQ. 0) THEN 02041 CALL DDSTP(TN,Y,YPRIME,NEQ, 02042 * RES,JAC,PSOL,H,RWORK(LWT),RWORK(LVT),INFO(1),IDID,RPAR,IPAR, 02043 * RWORK(LPHI),RWORK(LSAVR),RWORK(LDELTA),RWORK(LE), 02044 * RWORK(LWM),IWORK(LIWM), 02045 * RWORK(LALPHA),RWORK(LBETA),RWORK(LGAMMA), 02046 * RWORK(LPSI),RWORK(LSIGMA), 02047 * RWORK(LCJ),RWORK(LCJOLD),RWORK(LHOLD),RWORK(LS),HMIN, 02048 * RWORK(LROUND), RWORK(LEPLI),RWORK(LSQRN),RWORK(LRSQRN), 02049 * RWORK(LEPCON), IWORK(LPHASE),IWORK(LJCALC),INFO(15), 02050 * IWORK(LK), IWORK(LKOLD),IWORK(LNS),NONNEG,INFO(12), 02051 * DNEDD) 02052 ELSE IF (INFO(12) .EQ. 1) THEN 02053 CALL DDSTP(TN,Y,YPRIME,NEQ, 02054 * RES,JAC,PSOL,H,RWORK(LWT),RWORK(LVT),INFO(1),IDID,RPAR,IPAR, 02055 * RWORK(LPHI),RWORK(LSAVR),RWORK(LDELTA),RWORK(LE), 02056 * RWORK(LWM),IWORK(LIWM), 02057 * RWORK(LALPHA),RWORK(LBETA),RWORK(LGAMMA), 02058 * RWORK(LPSI),RWORK(LSIGMA), 02059 * RWORK(LCJ),RWORK(LCJOLD),RWORK(LHOLD),RWORK(LS),HMIN, 02060 * RWORK(LROUND), RWORK(LEPLI),RWORK(LSQRN),RWORK(LRSQRN), 02061 * RWORK(LEPCON), IWORK(LPHASE),IWORK(LJCALC),INFO(15), 02062 * IWORK(LK), IWORK(LKOLD),IWORK(LNS),NONNEG,INFO(12), 02063 * DNEDK) 02064 ENDIF 02065 C 02066 527 IF(IDID.LT.0)GO TO 600 02067 C 02068 C----------------------------------------------------------------------- 02069 C This block handles the case of a successful return from DDSTP 02070 C (IDID=1). Test for stop conditions. 02071 C----------------------------------------------------------------------- 02072 C 02073 IF(INFO(4).NE.0)GO TO 540 02074 IF(INFO(3).NE.0)GO TO 530 02075 IF((TN-TOUT)*H.LT.0.0D0)GO TO 500 02076 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ, 02077 * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI)) 02078 IDID=3 02079 T=TOUT 02080 GO TO 580 02081 530 IF((TN-TOUT)*H.GE.0.0D0)GO TO 535 02082 T=TN 02083 IDID=1 02084 GO TO 580 02085 535 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ, 02086 * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI)) 02087 IDID=3 02088 T=TOUT 02089 GO TO 580 02090 540 IF(INFO(3).NE.0)GO TO 550 02091 IF((TN-TOUT)*H.LT.0.0D0)GO TO 542 02092 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ, 02093 * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI)) 02094 T=TOUT 02095 IDID=3 02096 GO TO 580 02097 542 IF(ABS(TN-TSTOP).LE.100.0D0*UROUND* 02098 * (ABS(TN)+ABS(H)))GO TO 545 02099 TNEXT=TN+H 02100 IF((TNEXT-TSTOP)*H.LE.0.0D0)GO TO 500 02101 H=TSTOP-TN 02102 GO TO 500 02103 545 CALL DDATRP(TN,TSTOP,Y,YPRIME,NEQ, 02104 * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI)) 02105 IDID=2 02106 T=TSTOP 02107 GO TO 580 02108 550 IF((TN-TOUT)*H.GE.0.0D0)GO TO 555 02109 IF(ABS(TN-TSTOP).LE.100.0D0*UROUND*(ABS(TN)+ABS(H)))GO TO 552 02110 T=TN 02111 IDID=1 02112 GO TO 580 02113 552 CALL DDATRP(TN,TSTOP,Y,YPRIME,NEQ, 02114 * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI)) 02115 IDID=2 02116 T=TSTOP 02117 GO TO 580 02118 555 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ, 02119 * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI)) 02120 T=TOUT 02121 IDID=3 02122 580 CONTINUE 02123 C 02124 C----------------------------------------------------------------------- 02125 C All successful returns from DDASPK are made from this block. 02126 C----------------------------------------------------------------------- 02127 C 02128 590 CONTINUE 02129 RWORK(LTN)=TN 02130 RWORK(LH)=H 02131 RETURN 02132 C 02133 C----------------------------------------------------------------------- 02134 C This block handles all unsuccessful returns other than for 02135 C illegal input. 02136 C----------------------------------------------------------------------- 02137 C 02138 600 CONTINUE 02139 ITEMP = -IDID 02140 GO TO (610,620,630,700,655,640,650,660,670,675, 02141 * 680,685,690,695), ITEMP 02142 C 02143 C The maximum number of steps was taken before 02144 C reaching tout. 02145 C 02146 610 MSG = 'DASPK-- AT CURRENT T (=R1) 500 STEPS' 02147 CALL XERRWD(MSG,38,610,0,0,0,0,1,TN,0.0D0) 02148 MSG = 'DASPK-- TAKEN ON THIS CALL BEFORE REACHING TOUT' 02149 CALL XERRWD(MSG,48,611,0,0,0,0,0,0.0D0,0.0D0) 02150 GO TO 700 02151 C 02152 C Too much accuracy for machine precision. 02153 C 02154 620 MSG = 'DASPK-- AT T (=R1) TOO MUCH ACCURACY REQUESTED' 02155 CALL XERRWD(MSG,47,620,0,0,0,0,1,TN,0.0D0) 02156 MSG = 'DASPK-- FOR PRECISION OF MACHINE. RTOL AND ATOL' 02157 CALL XERRWD(MSG,48,621,0,0,0,0,0,0.0D0,0.0D0) 02158 MSG = 'DASPK-- WERE INCREASED TO APPROPRIATE VALUES' 02159 CALL XERRWD(MSG,45,622,0,0,0,0,0,0.0D0,0.0D0) 02160 GO TO 700 02161 C 02162 C WT(I) .LE. 0.0D0 for some I (not at start of problem). 02163 C 02164 630 MSG = 'DASPK-- AT T (=R1) SOME ELEMENT OF WT' 02165 CALL XERRWD(MSG,38,630,0,0,0,0,1,TN,0.0D0) 02166 MSG = 'DASPK-- HAS BECOME .LE. 0.0' 02167 CALL XERRWD(MSG,28,631,0,0,0,0,0,0.0D0,0.0D0) 02168 GO TO 700 02169 C 02170 C Error test failed repeatedly or with H=HMIN. 02171 C 02172 640 MSG = 'DASPK-- AT T (=R1) AND STEPSIZE H (=R2) THE' 02173 CALL XERRWD(MSG,44,640,0,0,0,0,2,TN,H) 02174 MSG='DASPK-- ERROR TEST FAILED REPEATEDLY OR WITH ABS(H)=HMIN' 02175 CALL XERRWD(MSG,57,641,0,0,0,0,0,0.0D0,0.0D0) 02176 GO TO 700 02177 C 02178 C Nonlinear solver failed to converge repeatedly or with H=HMIN. 02179 C 02180 650 MSG = 'DASPK-- AT T (=R1) AND STEPSIZE H (=R2) THE' 02181 CALL XERRWD(MSG,44,650,0,0,0,0,2,TN,H) 02182 MSG = 'DASPK-- NONLINEAR SOLVER FAILED TO CONVERGE' 02183 CALL XERRWD(MSG,44,651,0,0,0,0,0,0.0D0,0.0D0) 02184 MSG = 'DASPK-- REPEATEDLY OR WITH ABS(H)=HMIN' 02185 CALL XERRWD(MSG,40,652,0,0,0,0,0,0.0D0,0.0D0) 02186 GO TO 700 02187 C 02188 C The preconditioner had repeated failures. 02189 C 02190 655 MSG = 'DASPK-- AT T (=R1) AND STEPSIZE H (=R2) THE' 02191 CALL XERRWD(MSG,44,655,0,0,0,0,2,TN,H) 02192 MSG = 'DASPK-- PRECONDITIONER HAD REPEATED FAILURES.' 02193 CALL XERRWD(MSG,46,656,0,0,0,0,0,0.0D0,0.0D0) 02194 GO TO 700 02195 C 02196 C The iteration matrix is singular. 02197 C 02198 660 MSG = 'DASPK-- AT T (=R1) AND STEPSIZE H (=R2) THE' 02199 CALL XERRWD(MSG,44,660,0,0,0,0,2,TN,H) 02200 MSG = 'DASPK-- ITERATION MATRIX IS SINGULAR.' 02201 CALL XERRWD(MSG,38,661,0,0,0,0,0,0.0D0,0.0D0) 02202 GO TO 700 02203 C 02204 C Nonlinear system failure preceded by error test failures. 02205 C 02206 670 MSG = 'DASPK-- AT T (=R1) AND STEPSIZE H (=R2) THE' 02207 CALL XERRWD(MSG,44,670,0,0,0,0,2,TN,H) 02208 MSG = 'DASPK-- NONLINEAR SOLVER COULD NOT CONVERGE.' 02209 CALL XERRWD(MSG,45,671,0,0,0,0,0,0.0D0,0.0D0) 02210 MSG = 'DASPK-- ALSO, THE ERROR TEST FAILED REPEATEDLY.' 02211 CALL XERRWD(MSG,49,672,0,0,0,0,0,0.0D0,0.0D0) 02212 GO TO 700 02213 C 02214 C Nonlinear system failure because IRES = -1. 02215 C 02216 675 MSG = 'DASPK-- AT T (=R1) AND STEPSIZE H (=R2) THE' 02217 CALL XERRWD(MSG,44,675,0,0,0,0,2,TN,H) 02218 MSG = 'DASPK-- NONLINEAR SYSTEM SOLVER COULD NOT CONVERGE' 02219 CALL XERRWD(MSG,51,676,0,0,0,0,0,0.0D0,0.0D0) 02220 MSG = 'DASPK-- BECAUSE IRES WAS EQUAL TO MINUS ONE' 02221 CALL XERRWD(MSG,44,677,0,0,0,0,0,0.0D0,0.0D0) 02222 GO TO 700 02223 C 02224 C Failure because IRES = -2. 02225 C 02226 680 MSG = 'DASPK-- AT T (=R1) AND STEPSIZE H (=R2)' 02227 CALL XERRWD(MSG,40,680,0,0,0,0,2,TN,H) 02228 MSG = 'DASPK-- IRES WAS EQUAL TO MINUS TWO' 02229 CALL XERRWD(MSG,36,681,0,0,0,0,0,0.0D0,0.0D0) 02230 GO TO 700 02231 C 02232 C Failed to compute initial YPRIME. 02233 C 02234 685 MSG = 'DASPK-- AT T (=R1) AND STEPSIZE H (=R2) THE' 02235 CALL XERRWD(MSG,44,685,0,0,0,0,0,0.0D0,0.0D0) 02236 MSG = 'DASPK-- INITIAL (Y,YPRIME) COULD NOT BE COMPUTED' 02237 CALL XERRWD(MSG,49,686,0,0,0,0,2,TN,H0) 02238 GO TO 700 02239 C 02240 C Failure because IER was negative from PSOL. 02241 C 02242 690 MSG = 'DASPK-- AT T (=R1) AND STEPSIZE H (=R2)' 02243 CALL XERRWD(MSG,40,690,0,0,0,0,2,TN,H) 02244 MSG = 'DASPK-- IER WAS NEGATIVE FROM PSOL' 02245 CALL XERRWD(MSG,35,691,0,0,0,0,0,0.0D0,0.0D0) 02246 GO TO 700 02247 C 02248 C Failure because the linear system solver could not converge. 02249 C 02250 695 MSG = 'DASPK-- AT T (=R1) AND STEPSIZE H (=R2) THE' 02251 CALL XERRWD(MSG,44,695,0,0,0,0,2,TN,H) 02252 MSG = 'DASPK-- LINEAR SYSTEM SOLVER COULD NOT CONVERGE.' 02253 CALL XERRWD(MSG,50,696,0,0,0,0,0,0.0D0,0.0D0) 02254 GO TO 700 02255 C 02256 C 02257 700 CONTINUE 02258 INFO(1)=-1 02259 T=TN 02260 RWORK(LTN)=TN 02261 RWORK(LH)=H 02262 RETURN 02263 C 02264 C----------------------------------------------------------------------- 02265 C This block handles all error returns due to illegal input, 02266 C as detected before calling DDSTP. 02267 C First the error message routine is called. If this happens 02268 C twice in succession, execution is terminated. 02269 C----------------------------------------------------------------------- 02270 C 02271 701 MSG = 'DASPK-- ELEMENT (=I1) OF INFO VECTOR IS NOT VALID' 02272 CALL XERRWD(MSG,50,1,0,1,ITEMP,0,0,0.0D0,0.0D0) 02273 GO TO 750 02274 702 MSG = 'DASPK-- NEQ (=I1) .LE. 0' 02275 CALL XERRWD(MSG,25,2,0,1,NEQ,0,0,0.0D0,0.0D0) 02276 GO TO 750 02277 703 MSG = 'DASPK-- MAXORD (=I1) NOT IN RANGE' 02278 CALL XERRWD(MSG,34,3,0,1,MXORD,0,0,0.0D0,0.0D0) 02279 GO TO 750 02280 704 MSG='DASPK-- RWORK LENGTH NEEDED, LENRW (=I1), EXCEEDS LRW (=I2)' 02281 CALL XERRWD(MSG,60,4,0,2,LENRW,LRW,0,0.0D0,0.0D0) 02282 GO TO 750 02283 705 MSG='DASPK-- IWORK LENGTH NEEDED, LENIW (=I1), EXCEEDS LIW (=I2)' 02284 CALL XERRWD(MSG,60,5,0,2,LENIW,LIW,0,0.0D0,0.0D0) 02285 GO TO 750 02286 706 MSG = 'DASPK-- SOME ELEMENT OF RTOL IS .LT. 0' 02287 CALL XERRWD(MSG,39,6,0,0,0,0,0,0.0D0,0.0D0) 02288 GO TO 750 02289 707 MSG = 'DASPK-- SOME ELEMENT OF ATOL IS .LT. 0' 02290 CALL XERRWD(MSG,39,7,0,0,0,0,0,0.0D0,0.0D0) 02291 GO TO 750 02292 708 MSG = 'DASPK-- ALL ELEMENTS OF RTOL AND ATOL ARE ZERO' 02293 CALL XERRWD(MSG,47,8,0,0,0,0,0,0.0D0,0.0D0) 02294 GO TO 750 02295 709 MSG='DASPK-- INFO(4) = 1 AND TSTOP (=R1) BEHIND TOUT (=R2)' 02296 CALL XERRWD(MSG,54,9,0,0,0,0,2,TSTOP,TOUT) 02297 GO TO 750 02298 710 MSG = 'DASPK-- HMAX (=R1) .LT. 0.0' 02299 CALL XERRWD(MSG,28,10,0,0,0,0,1,HMAX,0.0D0) 02300 GO TO 750 02301 711 MSG = 'DASPK-- TOUT (=R1) BEHIND T (=R2)' 02302 CALL XERRWD(MSG,34,11,0,0,0,0,2,TOUT,T) 02303 GO TO 750 02304 712 MSG = 'DASPK-- INFO(8)=1 AND H0=0.0' 02305 CALL XERRWD(MSG,29,12,0,0,0,0,0,0.0D0,0.0D0) 02306 GO TO 750 02307 713 MSG = 'DASPK-- SOME ELEMENT OF WT IS .LE. 0.0' 02308 CALL XERRWD(MSG,39,13,0,0,0,0,0,0.0D0,0.0D0) 02309 GO TO 750 02310 714 MSG='DASPK-- TOUT (=R1) TOO CLOSE TO T (=R2) TO START INTEGRATION' 02311 CALL XERRWD(MSG,60,14,0,0,0,0,2,TOUT,T) 02312 GO TO 750 02313 715 MSG = 'DASPK-- INFO(4)=1 AND TSTOP (=R1) BEHIND T (=R2)' 02314 CALL XERRWD(MSG,49,15,0,0,0,0,2,TSTOP,T) 02315 GO TO 750 02316 717 MSG = 'DASPK-- ML (=I1) ILLEGAL. EITHER .LT. 0 OR .GT. NEQ' 02317 CALL XERRWD(MSG,52,17,0,1,IWORK(LML),0,0,0.0D0,0.0D0) 02318 GO TO 750 02319 718 MSG = 'DASPK-- MU (=I1) ILLEGAL. EITHER .LT. 0 OR .GT. NEQ' 02320 CALL XERRWD(MSG,52,18,0,1,IWORK(LMU),0,0,0.0D0,0.0D0) 02321 GO TO 750 02322 719 MSG = 'DASPK-- TOUT (=R1) IS EQUAL TO T (=R2)' 02323 CALL XERRWD(MSG,39,19,0,0,0,0,2,TOUT,T) 02324 GO TO 750 02325 720 MSG = 'DASPK-- MAXL (=I1) ILLEGAL. EITHER .LT. 1 OR .GT. NEQ' 02326 CALL XERRWD(MSG,54,20,0,1,IWORK(LMAXL),0,0,0.0D0,0.0D0) 02327 GO TO 750 02328 721 MSG = 'DASPK-- KMP (=I1) ILLEGAL. EITHER .LT. 1 OR .GT. MAXL' 02329 CALL XERRWD(MSG,54,21,0,1,IWORK(LKMP),0,0,0.0D0,0.0D0) 02330 GO TO 750 02331 722 MSG = 'DASPK-- NRMAX (=I1) ILLEGAL. .LT. 0' 02332 CALL XERRWD(MSG,36,22,0,1,IWORK(LNRMAX),0,0,0.0D0,0.0D0) 02333 GO TO 750 02334 723 MSG = 'DASPK-- EPLI (=R1) ILLEGAL. EITHER .LE. 0.D0 OR .GE. 1.D0' 02335 CALL XERRWD(MSG,58,23,0,0,0,0,1,RWORK(LEPLI),0.0D0) 02336 GO TO 750 02337 724 MSG = 'DASPK-- ILLEGAL IWORK VALUE FOR INFO(11) .NE. 0' 02338 CALL XERRWD(MSG,48,24,0,0,0,0,0,0.0D0,0.0D0) 02339 GO TO 750 02340 725 MSG = 'DASPK-- ONE OF THE INPUTS FOR INFO(17) = 1 IS ILLEGAL' 02341 CALL XERRWD(MSG,54,25,0,0,0,0,0,0.0D0,0.0D0) 02342 GO TO 750 02343 726 MSG = 'DASPK-- ILLEGAL IWORK VALUE FOR INFO(10) .NE. 0' 02344 CALL XERRWD(MSG,48,26,0,0,0,0,0,0.0D0,0.0D0) 02345 GO TO 750 02346 727 MSG = 'DASPK-- Y(I) AND IWORK(40+I) (I=I1) INCONSISTENT' 02347 CALL XERRWD(MSG,49,27,0,1,IRET,0,0,0.0D0,0.0D0) 02348 GO TO 750 02349 750 IF(INFO(1).EQ.-1) GO TO 760 02350 INFO(1)=-1 02351 IDID=-33 02352 RETURN 02353 760 MSG = 'DASPK-- REPEATED OCCURRENCES OF ILLEGAL INPUT' 02354 CALL XERRWD(MSG,46,701,0,0,0,0,0,0.0D0,0.0D0) 02355 770 MSG = 'DASPK-- RUN TERMINATED. APPARENT INFINITE LOOP' 02356 CALL XERRWD(MSG,47,702,1,0,0,0,0,0.0D0,0.0D0) 02357 RETURN 02358 C 02359 C------END OF SUBROUTINE DDASPK----------------------------------------- 02360 END
1.7.1
