|
3911
|
1 C Work performed under the auspices of the U.S. Department of Energy |
|
|
2 C by Lawrence Livermore National Laboratory under contract number |
|
|
3 C W-7405-Eng-48. |
|
|
4 C |
|
|
5 SUBROUTINE DDASPK (RES, NEQ, T, Y, YPRIME, TOUT, INFO, RTOL, ATOL, |
|
|
6 * IDID, RWORK, LRW, IWORK, LIW, RPAR, IPAR, JAC, PSOL) |
|
|
7 C |
|
|
8 C***BEGIN PROLOGUE DDASPK |
|
|
9 C***DATE WRITTEN 890101 (YYMMDD) |
|
|
10 C***REVISION DATE 910624 |
|
|
11 C***REVISION DATE 920929 (CJ in RES call, RES counter fix.) |
|
|
12 C***REVISION DATE 921215 (Warnings on poor iteration performance) |
|
|
13 C***REVISION DATE 921216 (NRMAX as optional input) |
|
|
14 C***REVISION DATE 930315 (Name change: DDINI to DDINIT) |
|
|
15 C***REVISION DATE 940822 (Replaced initial condition calculation) |
|
|
16 C***REVISION DATE 941101 (Added linesearch in I.C. calculations) |
|
|
17 C***REVISION DATE 941220 (Misc. corrections throughout) |
|
|
18 C***REVISION DATE 950125 (Added DINVWT routine) |
|
|
19 C***REVISION DATE 950714 (Misc. corrections throughout) |
|
|
20 C***REVISION DATE 950802 (Default NRMAX = 5, based on tests.) |
|
|
21 C***REVISION DATE 950808 (Optional error test added.) |
|
|
22 C***REVISION DATE 950814 (Added I.C. constraints and INFO(14)) |
|
|
23 C***REVISION DATE 950828 (Various minor corrections.) |
|
|
24 C***REVISION DATE 951006 (Corrected WT scaling in DFNRMK.) |
|
|
25 C***REVISION DATE 960129 (Corrected RL bug in DLINSD, DLINSK.) |
|
|
26 C***REVISION DATE 960301 (Added NONNEG to SAVE statement.) |
|
|
27 C***CATEGORY NO. I1A2 |
|
|
28 C***KEYWORDS DIFFERENTIAL/ALGEBRAIC, BACKWARD DIFFERENTIATION FORMULAS, |
|
|
29 C IMPLICIT DIFFERENTIAL SYSTEMS, KRYLOV ITERATION |
|
|
30 C***AUTHORS Linda R. Petzold, Peter N. Brown, Alan C. Hindmarsh, and |
|
|
31 C Clement W. Ulrich |
|
|
32 C Center for Computational Sciences & Engineering, L-316 |
|
|
33 C Lawrence Livermore National Laboratory |
|
|
34 C P.O. Box 808, |
|
|
35 C Livermore, CA 94551 |
|
|
36 C***PURPOSE This code solves a system of differential/algebraic |
|
|
37 C equations of the form |
|
|
38 C G(t,y,y') = 0 , |
|
|
39 C using a combination of Backward Differentiation Formula |
|
|
40 C (BDF) methods and a choice of two linear system solution |
|
|
41 C methods: direct (dense or band) or Krylov (iterative). |
|
|
42 C This version is in double precision. |
|
|
43 C----------------------------------------------------------------------- |
|
|
44 C***DESCRIPTION |
|
|
45 C |
|
|
46 C *Usage: |
|
|
47 C |
|
|
48 C IMPLICIT DOUBLE PRECISION(A-H,O-Z) |
|
|
49 C INTEGER NEQ, INFO(N), IDID, LRW, LIW, IWORK(LIW), IPAR(*) |
|
|
50 C DOUBLE PRECISION T, Y(*), YPRIME(*), TOUT, RTOL(*), ATOL(*), |
|
|
51 C RWORK(LRW), RPAR(*) |
|
|
52 C EXTERNAL RES, JAC, PSOL |
|
|
53 C |
|
|
54 C CALL DDASPK (RES, NEQ, T, Y, YPRIME, TOUT, INFO, RTOL, ATOL, |
|
|
55 C * IDID, RWORK, LRW, IWORK, LIW, RPAR, IPAR, JAC, PSOL) |
|
|
56 C |
|
|
57 C Quantities which may be altered by the code are: |
|
|
58 C T, Y(*), YPRIME(*), INFO(1), RTOL, ATOL, IDID, RWORK(*), IWORK(*) |
|
|
59 C |
|
|
60 C |
|
|
61 C *Arguments: |
|
|
62 C |
|
|
63 C RES:EXT This is the name of a subroutine which you |
|
|
64 C provide to define the residual function G(t,y,y') |
|
|
65 C of the differential/algebraic system. |
|
|
66 C |
|
|
67 C NEQ:IN This is the number of equations in the system. |
|
|
68 C |
|
|
69 C T:INOUT This is the current value of the independent |
|
|
70 C variable. |
|
|
71 C |
|
|
72 C Y(*):INOUT This array contains the solution components at T. |
|
|
73 C |
|
|
74 C YPRIME(*):INOUT This array contains the derivatives of the solution |
|
|
75 C components at T. |
|
|
76 C |
|
|
77 C TOUT:IN This is a point at which a solution is desired. |
|
|
78 C |
|
|
79 C INFO(N):IN This is an integer array used to communicate details |
|
|
80 C of how the solution is to be carried out, such as |
|
|
81 C tolerance type, matrix structure, step size and |
|
|
82 C order limits, and choice of nonlinear system method. |
|
|
83 C N must be at least 20. |
|
|
84 C |
|
|
85 C RTOL,ATOL:INOUT These quantities represent absolute and relative |
|
|
86 C error tolerances (on local error) which you provide |
|
|
87 C to indicate how accurately you wish the solution to |
|
|
88 C be computed. You may choose them to be both scalars |
|
|
89 C or else both arrays of length NEQ. |
|
|
90 C |
|
|
91 C IDID:OUT This integer scalar is an indicator reporting what |
|
|
92 C the code did. You must monitor this variable to |
|
|
93 C decide what action to take next. |
|
|
94 C |
|
|
95 C RWORK:WORK A real work array of length LRW which provides the |
|
|
96 C code with needed storage space. |
|
|
97 C |
|
|
98 C LRW:IN The length of RWORK. |
|
|
99 C |
|
|
100 C IWORK:WORK An integer work array of length LIW which provides |
|
|
101 C the code with needed storage space. |
|
|
102 C |
|
|
103 C LIW:IN The length of IWORK. |
|
|
104 C |
|
|
105 C RPAR,IPAR:IN These are real and integer parameter arrays which |
|
|
106 C you can use for communication between your calling |
|
|
107 C program and the RES, JAC, and PSOL subroutines. |
|
|
108 C |
|
|
109 C JAC:EXT This is the name of a subroutine which you may |
|
|
110 C provide (optionally) for calculating Jacobian |
|
|
111 C (partial derivative) data involved in solving linear |
|
|
112 C systems within DDASPK. |
|
|
113 C |
|
|
114 C PSOL:EXT This is the name of a subroutine which you must |
|
|
115 C provide for solving linear systems if you selected |
|
|
116 C a Krylov method. The purpose of PSOL is to solve |
|
|
117 C linear systems involving a left preconditioner P. |
|
|
118 C |
|
|
119 C *Overview |
|
|
120 C |
|
|
121 C The DDASPK solver uses the backward differentiation formulas of |
|
|
122 C orders one through five to solve a system of the form G(t,y,y') = 0 |
|
|
123 C for y = Y and y' = YPRIME. Values for Y and YPRIME at the initial |
|
|
124 C time must be given as input. These values should be consistent, |
|
|
125 C that is, if T, Y, YPRIME are the given initial values, they should |
|
|
126 C satisfy G(T,Y,YPRIME) = 0. However, if consistent values are not |
|
|
127 C known, in many cases you can have DDASPK solve for them -- see INFO(11). |
|
|
128 C (This and other options are described in more detail below.) |
|
|
129 C |
|
|
130 C Normally, DDASPK solves the system from T to TOUT. It is easy to |
|
|
131 C continue the solution to get results at additional TOUT. This is |
|
|
132 C the interval mode of operation. Intermediate results can also be |
|
|
133 C obtained easily by specifying INFO(3). |
|
|
134 C |
|
|
135 C On each step taken by DDASPK, a sequence of nonlinear algebraic |
|
|
136 C systems arises. These are solved by one of two types of |
|
|
137 C methods: |
|
|
138 C * a Newton iteration with a direct method for the linear |
|
|
139 C systems involved (INFO(12) = 0), or |
|
|
140 C * a Newton iteration with a preconditioned Krylov iterative |
|
|
141 C method for the linear systems involved (INFO(12) = 1). |
|
|
142 C |
|
|
143 C The direct method choices are dense and band matrix solvers, |
|
|
144 C with either a user-supplied or an internal difference quotient |
|
|
145 C Jacobian matrix, as specified by INFO(5) and INFO(6). |
|
|
146 C In the band case, INFO(6) = 1, you must supply half-bandwidths |
|
|
147 C in IWORK(1) and IWORK(2). |
|
|
148 C |
|
|
149 C The Krylov method is the Generalized Minimum Residual (GMRES) |
|
|
150 C method, in either complete or incomplete form, and with |
|
|
151 C scaling and preconditioning. The method is implemented |
|
|
152 C in an algorithm called SPIGMR. Certain options in the Krylov |
|
|
153 C method case are specified by INFO(13) and INFO(15). |
|
|
154 C |
|
|
155 C If the Krylov method is chosen, you may supply a pair of routines, |
|
|
156 C JAC and PSOL, to apply preconditioning to the linear system. |
|
|
157 C If the system is A*x = b, the matrix is A = dG/dY + CJ*dG/dYPRIME |
|
|
158 C (of order NEQ). This system can then be preconditioned in the form |
|
|
159 C (P-inverse)*A*x = (P-inverse)*b, with left preconditioner P. |
|
|
160 C (DDASPK does not allow right preconditioning.) |
|
|
161 C Then the Krylov method is applied to this altered, but equivalent, |
|
|
162 C linear system, hopefully with much better performance than without |
|
|
163 C preconditioning. (In addition, a diagonal scaling matrix based on |
|
|
164 C the tolerances is also introduced into the altered system.) |
|
|
165 C |
|
|
166 C The JAC routine evaluates any data needed for solving systems |
|
|
167 C with coefficient matrix P, and PSOL carries out that solution. |
|
|
168 C In any case, in order to improve convergence, you should try to |
|
|
169 C make P approximate the matrix A as much as possible, while keeping |
|
|
170 C the system P*x = b reasonably easy and inexpensive to solve for x, |
|
|
171 C given a vector b. |
|
|
172 C |
|
|
173 C |
|
|
174 C *Description |
|
|
175 C |
|
|
176 C------INPUT - WHAT TO DO ON THE FIRST CALL TO DDASPK------------------- |
|
|
177 C |
|
|
178 C |
|
|
179 C The first call of the code is defined to be the start of each new |
|
|
180 C problem. Read through the descriptions of all the following items, |
|
|
181 C provide sufficient storage space for designated arrays, set |
|
|
182 C appropriate variables for the initialization of the problem, and |
|
|
183 C give information about how you want the problem to be solved. |
|
|
184 C |
|
|
185 C |
|
|
186 C RES -- Provide a subroutine of the form |
|
|
187 C |
|
|
188 C SUBROUTINE RES (T, Y, YPRIME, CJ, DELTA, IRES, RPAR, IPAR) |
|
|
189 C |
|
|
190 C to define the system of differential/algebraic |
|
|
191 C equations which is to be solved. For the given values |
|
|
192 C of T, Y and YPRIME, the subroutine should return |
|
|
193 C the residual of the differential/algebraic system |
|
|
194 C DELTA = G(T,Y,YPRIME) |
|
|
195 C DELTA is a vector of length NEQ which is output from RES. |
|
|
196 C |
|
|
197 C Subroutine RES must not alter T, Y, YPRIME, or CJ. |
|
|
198 C You must declare the name RES in an EXTERNAL |
|
|
199 C statement in your program that calls DDASPK. |
|
|
200 C You must dimension Y, YPRIME, and DELTA in RES. |
|
|
201 C |
|
|
202 C The input argument CJ can be ignored, or used to rescale |
|
|
203 C constraint equations in the system (see Ref. 2, p. 145). |
|
|
204 C Note: In this respect, DDASPK is not downward-compatible |
|
|
205 C with DDASSL, which does not have the RES argument CJ. |
|
|
206 C |
|
|
207 C IRES is an integer flag which is always equal to zero |
|
|
208 C on input. Subroutine RES should alter IRES only if it |
|
|
209 C encounters an illegal value of Y or a stop condition. |
|
|
210 C Set IRES = -1 if an input value is illegal, and DDASPK |
|
|
211 C will try to solve the problem without getting IRES = -1. |
|
|
212 C If IRES = -2, DDASPK will return control to the calling |
|
|
213 C program with IDID = -11. |
|
|
214 C |
|
|
215 C RPAR and IPAR are real and integer parameter arrays which |
|
|
216 C you can use for communication between your calling program |
|
|
217 C and subroutine RES. They are not altered by DDASPK. If you |
|
|
218 C do not need RPAR or IPAR, ignore these parameters by treat- |
|
|
219 C ing them as dummy arguments. If you do choose to use them, |
|
|
220 C dimension them in your calling program and in RES as arrays |
|
|
221 C of appropriate length. |
|
|
222 C |
|
|
223 C NEQ -- Set it to the number of equations in the system (NEQ .GE. 1). |
|
|
224 C |
|
|
225 C T -- Set it to the initial point of the integration. (T must be |
|
|
226 C a variable.) |
|
|
227 C |
|
|
228 C Y(*) -- Set this array to the initial values of the NEQ solution |
|
|
229 C components at the initial point. You must dimension Y of |
|
|
230 C length at least NEQ in your calling program. |
|
|
231 C |
|
|
232 C YPRIME(*) -- Set this array to the initial values of the NEQ first |
|
|
233 C derivatives of the solution components at the initial |
|
|
234 C point. You must dimension YPRIME at least NEQ in your |
|
|
235 C calling program. |
|
|
236 C |
|
|
237 C TOUT - Set it to the first point at which a solution is desired. |
|
|
238 C You cannot take TOUT = T. Integration either forward in T |
|
|
239 C (TOUT .GT. T) or backward in T (TOUT .LT. T) is permitted. |
|
|
240 C |
|
|
241 C The code advances the solution from T to TOUT using step |
|
|
242 C sizes which are automatically selected so as to achieve the |
|
|
243 C desired accuracy. If you wish, the code will return with the |
|
|
244 C solution and its derivative at intermediate steps (the |
|
|
245 C intermediate-output mode) so that you can monitor them, |
|
|
246 C but you still must provide TOUT in accord with the basic |
|
|
247 C aim of the code. |
|
|
248 C |
|
|
249 C The first step taken by the code is a critical one because |
|
|
250 C it must reflect how fast the solution changes near the |
|
|
251 C initial point. The code automatically selects an initial |
|
|
252 C step size which is practically always suitable for the |
|
|
253 C problem. By using the fact that the code will not step past |
|
|
254 C TOUT in the first step, you could, if necessary, restrict the |
|
|
255 C length of the initial step. |
|
|
256 C |
|
|
257 C For some problems it may not be permissible to integrate |
|
|
258 C past a point TSTOP, because a discontinuity occurs there |
|
|
259 C or the solution or its derivative is not defined beyond |
|
|
260 C TSTOP. When you have declared a TSTOP point (see INFO(4) |
|
|
261 C and RWORK(1)), you have told the code not to integrate past |
|
|
262 C TSTOP. In this case any tout beyond TSTOP is invalid input. |
|
|
263 C |
|
|
264 C INFO(*) - Use the INFO array to give the code more details about |
|
|
265 C how you want your problem solved. This array should be |
|
|
266 C dimensioned of length 20, though DDASPK uses only the |
|
|
267 C first 15 entries. You must respond to all of the following |
|
|
268 C items, which are arranged as questions. The simplest use |
|
|
269 C of DDASPK corresponds to setting all entries of INFO to 0. |
|
|
270 C |
|
|
271 C INFO(1) - This parameter enables the code to initialize itself. |
|
|
272 C You must set it to indicate the start of every new |
|
|
273 C problem. |
|
|
274 C |
|
|
275 C **** Is this the first call for this problem ... |
|
|
276 C yes - set INFO(1) = 0 |
|
|
277 C no - not applicable here. |
|
|
278 C See below for continuation calls. **** |
|
|
279 C |
|
|
280 C INFO(2) - How much accuracy you want of your solution |
|
|
281 C is specified by the error tolerances RTOL and ATOL. |
|
|
282 C The simplest use is to take them both to be scalars. |
|
|
283 C To obtain more flexibility, they can both be arrays. |
|
|
284 C The code must be told your choice. |
|
|
285 C |
|
|
286 C **** Are both error tolerances RTOL, ATOL scalars ... |
|
|
287 C yes - set INFO(2) = 0 |
|
|
288 C and input scalars for both RTOL and ATOL |
|
|
289 C no - set INFO(2) = 1 |
|
|
290 C and input arrays for both RTOL and ATOL **** |
|
|
291 C |
|
|
292 C INFO(3) - The code integrates from T in the direction of TOUT |
|
|
293 C by steps. If you wish, it will return the computed |
|
|
294 C solution and derivative at the next intermediate step |
|
|
295 C (the intermediate-output mode) or TOUT, whichever comes |
|
|
296 C first. This is a good way to proceed if you want to |
|
|
297 C see the behavior of the solution. If you must have |
|
|
298 C solutions at a great many specific TOUT points, this |
|
|
299 C code will compute them efficiently. |
|
|
300 C |
|
|
301 C **** Do you want the solution only at |
|
|
302 C TOUT (and not at the next intermediate step) ... |
|
|
303 C yes - set INFO(3) = 0 |
|
|
304 C no - set INFO(3) = 1 **** |
|
|
305 C |
|
|
306 C INFO(4) - To handle solutions at a great many specific |
|
|
307 C values TOUT efficiently, this code may integrate past |
|
|
308 C TOUT and interpolate to obtain the result at TOUT. |
|
|
309 C Sometimes it is not possible to integrate beyond some |
|
|
310 C point TSTOP because the equation changes there or it is |
|
|
311 C not defined past TSTOP. Then you must tell the code |
|
|
312 C this stop condition. |
|
|
313 C |
|
|
314 C **** Can the integration be carried out without any |
|
|
315 C restrictions on the independent variable T ... |
|
|
316 C yes - set INFO(4) = 0 |
|
|
317 C no - set INFO(4) = 1 |
|
|
318 C and define the stopping point TSTOP by |
|
|
319 C setting RWORK(1) = TSTOP **** |
|
|
320 C |
|
|
321 C INFO(5) - used only when INFO(12) = 0 (direct methods). |
|
|
322 C To solve differential/algebraic systems you may wish |
|
|
323 C to use a matrix of partial derivatives of the |
|
|
324 C system of differential equations. If you do not |
|
|
325 C provide a subroutine to evaluate it analytically (see |
|
|
326 C description of the item JAC in the call list), it will |
|
|
327 C be approximated by numerical differencing in this code. |
|
|
328 C Although it is less trouble for you to have the code |
|
|
329 C compute partial derivatives by numerical differencing, |
|
|
330 C the solution will be more reliable if you provide the |
|
|
331 C derivatives via JAC. Usually numerical differencing is |
|
|
332 C more costly than evaluating derivatives in JAC, but |
|
|
333 C sometimes it is not - this depends on your problem. |
|
|
334 C |
|
|
335 C **** Do you want the code to evaluate the partial deriv- |
|
|
336 C atives automatically by numerical differences ... |
|
|
337 C yes - set INFO(5) = 0 |
|
|
338 C no - set INFO(5) = 1 |
|
|
339 C and provide subroutine JAC for evaluating the |
|
|
340 C matrix of partial derivatives **** |
|
|
341 C |
|
|
342 C INFO(6) - used only when INFO(12) = 0 (direct methods). |
|
|
343 C DDASPK will perform much better if the matrix of |
|
|
344 C partial derivatives, dG/dY + CJ*dG/dYPRIME (here CJ is |
|
|
345 C a scalar determined by DDASPK), is banded and the code |
|
|
346 C is told this. In this case, the storage needed will be |
|
|
347 C greatly reduced, numerical differencing will be performed |
|
|
348 C much cheaper, and a number of important algorithms will |
|
|
349 C execute much faster. The differential equation is said |
|
|
350 C to have half-bandwidths ML (lower) and MU (upper) if |
|
|
351 C equation i involves only unknowns Y(j) with |
|
|
352 C i-ML .le. j .le. i+MU . |
|
|
353 C For all i=1,2,...,NEQ. Thus, ML and MU are the widths |
|
|
354 C of the lower and upper parts of the band, respectively, |
|
|
355 C with the main diagonal being excluded. If you do not |
|
|
356 C indicate that the equation has a banded matrix of partial |
|
|
357 C derivatives the code works with a full matrix of NEQ**2 |
|
|
358 C elements (stored in the conventional way). Computations |
|
|
359 C with banded matrices cost less time and storage than with |
|
|
360 C full matrices if 2*ML+MU .lt. NEQ. If you tell the |
|
|
361 C code that the matrix of partial derivatives has a banded |
|
|
362 C structure and you want to provide subroutine JAC to |
|
|
363 C compute the partial derivatives, then you must be careful |
|
|
364 C to store the elements of the matrix in the special form |
|
|
365 C indicated in the description of JAC. |
|
|
366 C |
|
|
367 C **** Do you want to solve the problem using a full (dense) |
|
|
368 C matrix (and not a special banded structure) ... |
|
|
369 C yes - set INFO(6) = 0 |
|
|
370 C no - set INFO(6) = 1 |
|
|
371 C and provide the lower (ML) and upper (MU) |
|
|
372 C bandwidths by setting |
|
|
373 C IWORK(1)=ML |
|
|
374 C IWORK(2)=MU **** |
|
|
375 C |
|
|
376 C INFO(7) - You can specify a maximum (absolute value of) |
|
|
377 C stepsize, so that the code will avoid passing over very |
|
|
378 C large regions. |
|
|
379 C |
|
|
380 C **** Do you want the code to decide on its own the maximum |
|
|
381 C stepsize ... |
|
|
382 C yes - set INFO(7) = 0 |
|
|
383 C no - set INFO(7) = 1 |
|
|
384 C and define HMAX by setting |
|
|
385 C RWORK(2) = HMAX **** |
|
|
386 C |
|
|
387 C INFO(8) - Differential/algebraic problems may occasionally |
|
|
388 C suffer from severe scaling difficulties on the first |
|
|
389 C step. If you know a great deal about the scaling of |
|
|
390 C your problem, you can help to alleviate this problem |
|
|
391 C by specifying an initial stepsize H0. |
|
|
392 C |
|
|
393 C **** Do you want the code to define its own initial |
|
|
394 C stepsize ... |
|
|
395 C yes - set INFO(8) = 0 |
|
|
396 C no - set INFO(8) = 1 |
|
|
397 C and define H0 by setting |
|
|
398 C RWORK(3) = H0 **** |
|
|
399 C |
|
|
400 C INFO(9) - If storage is a severe problem, you can save some |
|
|
401 C storage by restricting the maximum method order MAXORD. |
|
|
402 C The default value is 5. For each order decrease below 5, |
|
|
403 C the code requires NEQ fewer locations, but it is likely |
|
|
404 C to be slower. In any case, you must have |
|
|
405 C 1 .le. MAXORD .le. 5. |
|
|
406 C **** Do you want the maximum order to default to 5 ... |
|
|
407 C yes - set INFO(9) = 0 |
|
|
408 C no - set INFO(9) = 1 |
|
|
409 C and define MAXORD by setting |
|
|
410 C IWORK(3) = MAXORD **** |
|
|
411 C |
|
|
412 C INFO(10) - If you know that certain components of the |
|
|
413 C solutions to your equations are always nonnegative |
|
|
414 C (or nonpositive), it may help to set this |
|
|
415 C parameter. There are three options that are |
|
|
416 C available: |
|
|
417 C 1. To have constraint checking only in the initial |
|
|
418 C condition calculation. |
|
|
419 C 2. To enforce nonnegativity in Y during the integration. |
|
|
420 C 3. To enforce both options 1 and 2. |
|
|
421 C |
|
|
422 C When selecting option 2 or 3, it is probably best to try the |
|
|
423 C code without using this option first, and only use |
|
|
424 C this option if that does not work very well. |
|
|
425 C |
|
|
426 C **** Do you want the code to solve the problem without |
|
|
427 C invoking any special inequality constraints ... |
|
|
428 C yes - set INFO(10) = 0 |
|
|
429 C no - set INFO(10) = 1 to have option 1 enforced |
|
|
430 C no - set INFO(10) = 2 to have option 2 enforced |
|
|
431 C no - set INFO(10) = 3 to have option 3 enforced **** |
|
|
432 C |
|
|
433 C If you have specified INFO(10) = 1 or 3, then you |
|
|
434 C will also need to identify how each component of Y |
|
|
435 C in the initial condition calculation is constrained. |
|
|
436 C You must set: |
|
|
437 C IWORK(40+I) = +1 if Y(I) must be .GE. 0, |
|
|
438 C IWORK(40+I) = +2 if Y(I) must be .GT. 0, |
|
|
439 C IWORK(40+I) = -1 if Y(I) must be .LE. 0, while |
|
|
440 C IWORK(40+I) = -2 if Y(I) must be .LT. 0, while |
|
|
441 C IWORK(40+I) = 0 if Y(I) is not constrained. |
|
|
442 C |
|
|
443 C INFO(11) - DDASPK normally requires the initial T, Y, and |
|
|
444 C YPRIME to be consistent. That is, you must have |
|
|
445 C G(T,Y,YPRIME) = 0 at the initial T. If you do not know |
|
|
446 C the initial conditions precisely, in some cases |
|
|
447 C DDASPK may be able to compute it. |
|
|
448 C |
|
|
449 C Denoting the differential variables in Y by Y_d |
|
|
450 C and the algebraic variables by Y_a, DDASPK can solve |
|
|
451 C one of two initialization problems: |
|
|
452 C 1. Given Y_d, calculate Y_a and Y'_d, or |
|
|
453 C 2. Given Y', calculate Y. |
|
|
454 C In either case, initial values for the given |
|
|
455 C components are input, and initial guesses for |
|
|
456 C the unknown components must also be provided as input. |
|
|
457 C |
|
|
458 C **** Are the initial T, Y, YPRIME consistent ... |
|
|
459 C |
|
|
460 C yes - set INFO(11) = 0 |
|
|
461 C no - set INFO(11) = 1 to calculate option 1 above, |
|
|
462 C or set INFO(11) = 2 to calculate option 2 **** |
|
|
463 C |
|
|
464 C If you have specified INFO(11) = 1, then you |
|
|
465 C will also need to identify which are the |
|
|
466 C differential and which are the algebraic |
|
|
467 C components (algebraic components are components |
|
|
468 C whose derivatives do not appear explicitly |
|
|
469 C in the function G(T,Y,YPRIME)). You must set: |
|
|
470 C IWORK(LID+I) = +1 if Y(I) is a differential variable |
|
|
471 C IWORK(LID+I) = -1 if Y(I) is an algebraic variable, |
|
|
472 C where LID = 40 if INFO(10) = 0 or 2 and LID = 40+NEQ |
|
|
473 C if INFO(10) = 1 or 3. |
|
|
474 C |
|
|
475 C INFO(12) - Except for the addition of the RES argument CJ, |
|
|
476 C DDASPK by default is downward-compatible with DDASSL, |
|
|
477 C which uses only direct (dense or band) methods to solve |
|
|
478 C the linear systems involved. You must set INFO(12) to |
|
|
479 C indicate whether you want the direct methods or the |
|
|
480 C Krylov iterative method. |
|
|
481 C **** Do you want DDASPK to use standard direct methods |
|
|
482 C (dense or band) or the Krylov (iterative) method ... |
|
|
483 C direct methods - set INFO(12) = 0. |
|
|
484 C Krylov method - set INFO(12) = 1, |
|
|
485 C and check the settings of INFO(13) and INFO(15). |
|
|
486 C |
|
|
487 C INFO(13) - used when INFO(12) = 1 (Krylov methods). |
|
|
488 C DDASPK uses scalars MAXL, KMP, NRMAX, and EPLI for the |
|
|
489 C iterative solution of linear systems. INFO(13) allows |
|
|
490 C you to override the default values of these parameters. |
|
|
491 C These parameters and their defaults are as follows: |
|
|
492 C MAXL = maximum number of iterations in the SPIGMR |
|
|
493 C algorithm (MAXL .le. NEQ). The default is |
|
|
494 C MAXL = MIN(5,NEQ). |
|
|
495 C KMP = number of vectors on which orthogonalization is |
|
|
496 C done in the SPIGMR algorithm. The default is |
|
|
497 C KMP = MAXL, which corresponds to complete GMRES |
|
|
498 C iteration, as opposed to the incomplete form. |
|
|
499 C NRMAX = maximum number of restarts of the SPIGMR |
|
|
500 C algorithm per nonlinear iteration. The default is |
|
|
501 C NRMAX = 5. |
|
|
502 C EPLI = convergence test constant in SPIGMR algorithm. |
|
|
503 C The default is EPLI = 0.05. |
|
|
504 C Note that the length of RWORK depends on both MAXL |
|
|
505 C and KMP. See the definition of LRW below. |
|
|
506 C **** Are MAXL, KMP, and EPLI to be given their |
|
|
507 C default values ... |
|
|
508 C yes - set INFO(13) = 0 |
|
|
509 C no - set INFO(13) = 1, |
|
|
510 C and set all of the following: |
|
|
511 C IWORK(24) = MAXL (1 .le. MAXL .le. NEQ) |
|
|
512 C IWORK(25) = KMP (1 .le. KMP .le. MAXL) |
|
|
513 C IWORK(26) = NRMAX (NRMAX .ge. 0) |
|
|
514 C RWORK(10) = EPLI (0 .lt. EPLI .lt. 1.0) **** |
|
|
515 C |
|
|
516 C INFO(14) - used with INFO(11) > 0 (initial condition |
|
|
517 C calculation is requested). In this case, you may |
|
|
518 C request control to be returned to the calling program |
|
|
519 C immediately after the initial condition calculation, |
|
|
520 C before proceeding to the integration of the system |
|
|
521 C (e.g. to examine the computed Y and YPRIME). |
|
|
522 C If this is done, and if the initialization succeeded |
|
|
523 C (IDID = 4), you should reset INFO(11) to 0 for the |
|
|
524 C next call, to prevent the solver from repeating the |
|
|
525 C initialization (and to avoid an infinite loop). |
|
|
526 C **** Do you want to proceed to the integration after |
|
|
527 C the initial condition calculation is done ... |
|
|
528 C yes - set INFO(14) = 0 |
|
|
529 C no - set INFO(14) = 1 **** |
|
|
530 C |
|
|
531 C INFO(15) - used when INFO(12) = 1 (Krylov methods). |
|
|
532 C When using preconditioning in the Krylov method, |
|
|
533 C you must supply a subroutine, PSOL, which solves the |
|
|
534 C associated linear systems using P. |
|
|
535 C The usage of DDASPK is simpler if PSOL can carry out |
|
|
536 C the solution without any prior calculation of data. |
|
|
537 C However, if some partial derivative data is to be |
|
|
538 C calculated in advance and used repeatedly in PSOL, |
|
|
539 C then you must supply a JAC routine to do this, |
|
|
540 C and set INFO(15) to indicate that JAC is to be called |
|
|
541 C for this purpose. For example, P might be an |
|
|
542 C approximation to a part of the matrix A which can be |
|
|
543 C calculated and LU-factored for repeated solutions of |
|
|
544 C the preconditioner system. The arrays WP and IWP |
|
|
545 C (described under JAC and PSOL) can be used to |
|
|
546 C communicate data between JAC and PSOL. |
|
|
547 C **** Does PSOL operate with no prior preparation ... |
|
|
548 C yes - set INFO(15) = 0 (no JAC routine) |
|
|
549 C no - set INFO(15) = 1 |
|
|
550 C and supply a JAC routine to evaluate and |
|
|
551 C preprocess any required Jacobian data. **** |
|
|
552 C |
|
|
553 C INFO(16) - option to exclude algebraic variables from |
|
|
554 C the error test. |
|
|
555 C **** Do you wish to control errors locally on |
|
|
556 C all the variables... |
|
|
557 C yes - set INFO(16) = 0 |
|
|
558 C no - set INFO(16) = 1 |
|
|
559 C If you have specified INFO(16) = 1, then you |
|
|
560 C will also need to identify which are the |
|
|
561 C differential and which are the algebraic |
|
|
562 C components (algebraic components are components |
|
|
563 C whose derivatives do not appear explicitly |
|
|
564 C in the function G(T,Y,YPRIME)). You must set: |
|
|
565 C IWORK(LID+I) = +1 if Y(I) is a differential |
|
|
566 C variable, and |
|
|
567 C IWORK(LID+I) = -1 if Y(I) is an algebraic |
|
|
568 C variable, |
|
|
569 C where LID = 40 if INFO(10) = 0 or 2 and |
|
|
570 C LID = 40 + NEQ if INFO(10) = 1 or 3. |
|
|
571 C |
|
|
572 C INFO(17) - used when INFO(11) > 0 (DDASPK is to do an |
|
|
573 C initial condition calculation). |
|
|
574 C DDASPK uses several heuristic control quantities in the |
|
|
575 C initial condition calculation. They have default values, |
|
|
576 C but can also be set by the user using INFO(17). |
|
|
577 C These parameters and their defaults are as follows: |
|
|
578 C MXNIT = maximum number of Newton iterations |
|
|
579 C per Jacobian or preconditioner evaluation. |
|
|
580 C The default is: |
|
|
581 C MXNIT = 5 in the direct case (INFO(12) = 0), and |
|
|
582 C MXNIT = 15 in the Krylov case (INFO(12) = 1). |
|
|
583 C MXNJ = maximum number of Jacobian or preconditioner |
|
|
584 C evaluations. The default is: |
|
|
585 C MXNJ = 6 in the direct case (INFO(12) = 0), and |
|
|
586 C MXNJ = 2 in the Krylov case (INFO(12) = 1). |
|
|
587 C MXNH = maximum number of values of the artificial |
|
|
588 C stepsize parameter H to be tried if INFO(11) = 1. |
|
|
589 C The default is MXNH = 5. |
|
|
590 C NOTE: the maximum number of Newton iterations |
|
|
591 C allowed in all is MXNIT*MXNJ*MXNH if INFO(11) = 1, |
|
|
592 C and MXNIT*MXNJ if INFO(11) = 2. |
|
|
593 C LSOFF = flag to turn off the linesearch algorithm |
|
|
594 C (LSOFF = 0 means linesearch is on, LSOFF = 1 means |
|
|
595 C it is turned off). The default is LSOFF = 0. |
|
|
596 C STPTOL = minimum scaled step in linesearch algorithm. |
|
|
597 C The default is STPTOL = (unit roundoff)**(2/3). |
|
|
598 C EPINIT = swing factor in the Newton iteration convergence |
|
|
599 C test. The test is applied to the residual vector, |
|
|
600 C premultiplied by the approximate Jacobian (in the |
|
|
601 C direct case) or the preconditioner (in the Krylov |
|
|
602 C case). For convergence, the weighted RMS norm of |
|
|
603 C this vector (scaled by the error weights) must be |
|
|
604 C less than EPINIT*EPCON, where EPCON = .33 is the |
|
|
605 C analogous test constant used in the time steps. |
|
|
606 C The default is EPINIT = .01. |
|
|
607 C **** Are the initial condition heuristic controls to be |
|
|
608 C given their default values... |
|
|
609 C yes - set INFO(17) = 0 |
|
|
610 C no - set INFO(17) = 1, |
|
|
611 C and set all of the following: |
|
|
612 C IWORK(32) = MXNIT (.GT. 0) |
|
|
613 C IWORK(33) = MXNJ (.GT. 0) |
|
|
614 C IWORK(34) = MXNH (.GT. 0) |
|
|
615 C IWORK(35) = LSOFF ( = 0 or 1) |
|
|
616 C RWORK(14) = STPTOL (.GT. 0.0) |
|
|
617 C RWORK(15) = EPINIT (.GT. 0.0) **** |
|
|
618 C |
|
|
619 C INFO(18) - option to get extra printing in initial condition |
|
|
620 C calculation. |
|
|
621 C **** Do you wish to have extra printing... |
|
|
622 C no - set INFO(18) = 0 |
|
|
623 C yes - set INFO(18) = 1 for minimal printing, or |
|
|
624 C set INFO(18) = 2 for full printing. |
|
|
625 C If you have specified INFO(18) .ge. 1, data |
|
|
626 C will be printed with the error handler routines. |
|
|
627 C To print to a non-default unit number L, include |
|
|
628 C the line CALL XSETUN(L) in your program. **** |
|
|
629 C |
|
|
630 C RTOL, ATOL -- You must assign relative (RTOL) and absolute (ATOL) |
|
|
631 C error tolerances to tell the code how accurately you |
|
|
632 C want the solution to be computed. They must be defined |
|
|
633 C as variables because the code may change them. |
|
|
634 C you have two choices -- |
|
|
635 C Both RTOL and ATOL are scalars (INFO(2) = 0), or |
|
|
636 C both RTOL and ATOL are vectors (INFO(2) = 1). |
|
|
637 C In either case all components must be non-negative. |
|
|
638 C |
|
|
639 C The tolerances are used by the code in a local error |
|
|
640 C test at each step which requires roughly that |
|
|
641 C abs(local error in Y(i)) .le. EWT(i) , |
|
|
642 C where EWT(i) = RTOL*abs(Y(i)) + ATOL is an error weight |
|
|
643 C quantity, for each vector component. |
|
|
644 C (More specifically, a root-mean-square norm is used to |
|
|
645 C measure the size of vectors, and the error test uses the |
|
|
646 C magnitude of the solution at the beginning of the step.) |
|
|
647 C |
|
|
648 C The true (global) error is the difference between the |
|
|
649 C true solution of the initial value problem and the |
|
|
650 C computed approximation. Practically all present day |
|
|
651 C codes, including this one, control the local error at |
|
|
652 C each step and do not even attempt to control the global |
|
|
653 C error directly. |
|
|
654 C |
|
|
655 C Usually, but not always, the true accuracy of |
|
|
656 C the computed Y is comparable to the error tolerances. |
|
|
657 C This code will usually, but not always, deliver a more |
|
|
658 C accurate solution if you reduce the tolerances and |
|
|
659 C integrate again. By comparing two such solutions you |
|
|
660 C can get a fairly reliable idea of the true error in the |
|
|
661 C solution at the larger tolerances. |
|
|
662 C |
|
|
663 C Setting ATOL = 0. results in a pure relative error test |
|
|
664 C on that component. Setting RTOL = 0. results in a pure |
|
|
665 C absolute error test on that component. A mixed test |
|
|
666 C with non-zero RTOL and ATOL corresponds roughly to a |
|
|
667 C relative error test when the solution component is |
|
|
668 C much bigger than ATOL and to an absolute error test |
|
|
669 C when the solution component is smaller than the |
|
|
670 C threshold ATOL. |
|
|
671 C |
|
|
672 C The code will not attempt to compute a solution at an |
|
|
673 C accuracy unreasonable for the machine being used. It |
|
|
674 C will advise you if you ask for too much accuracy and |
|
|
675 C inform you as to the maximum accuracy it believes |
|
|
676 C possible. |
|
|
677 C |
|
|
678 C RWORK(*) -- a real work array, which should be dimensioned in your |
|
|
679 C calling program with a length equal to the value of |
|
|
680 C LRW (or greater). |
|
|
681 C |
|
|
682 C LRW -- Set it to the declared length of the RWORK array. The |
|
|
683 C minimum length depends on the options you have selected, |
|
|
684 C given by a base value plus additional storage as described |
|
|
685 C below. |
|
|
686 C |
|
|
687 C If INFO(12) = 0 (standard direct method), the base value is |
|
|
688 C base = 50 + max(MAXORD+4,7)*NEQ. |
|
|
689 C The default value is MAXORD = 5 (see INFO(9)). With the |
|
|
690 C default MAXORD, base = 50 + 9*NEQ. |
|
|
691 C Additional storage must be added to the base value for |
|
|
692 C any or all of the following options: |
|
|
693 C if INFO(6) = 0 (dense matrix), add NEQ**2 |
|
|
694 C if INFO(6) = 1 (banded matrix), then |
|
|
695 C if INFO(5) = 0, add (2*ML+MU+1)*NEQ + 2*(NEQ/(ML+MU+1)+1), |
|
|
696 C if INFO(5) = 1, add (2*ML+MU+1)*NEQ, |
|
|
697 C if INFO(16) = 1, add NEQ. |
|
|
698 C |
|
|
699 C If INFO(12) = 1 (Krylov method), the base value is |
|
|
700 C base = 50 + (MAXORD+5)*NEQ + (MAXL+3+MIN0(1,MAXL-KMP))*NEQ + |
|
|
701 C + (MAXL+3)*MAXL + 1 + LENWP. |
|
|
702 C See PSOL for description of LENWP. The default values are: |
|
|
703 C MAXORD = 5 (see INFO(9)), MAXL = min(5,NEQ) and KMP = MAXL |
|
|
704 C (see INFO(13)). |
|
|
705 C With the default values for MAXORD, MAXL and KMP, |
|
|
706 C base = 91 + 18*NEQ + LENWP. |
|
|
707 C Additional storage must be added to the base value for |
|
|
708 C any or all of the following options: |
|
|
709 C if INFO(16) = 1, add NEQ. |
|
|
710 C |
|
|
711 C |
|
|
712 C IWORK(*) -- an integer work array, which should be dimensioned in |
|
|
713 C your calling program with a length equal to the value |
|
|
714 C of LIW (or greater). |
|
|
715 C |
|
|
716 C LIW -- Set it to the declared length of the IWORK array. The |
|
|
717 C minimum length depends on the options you have selected, |
|
|
718 C given by a base value plus additional storage as described |
|
|
719 C below. |
|
|
720 C |
|
|
721 C If INFO(12) = 0 (standard direct method), the base value is |
|
|
722 C base = 40 + NEQ. |
|
|
723 C IF INFO(10) = 1 or 3, add NEQ to the base value. |
|
|
724 C If INFO(11) = 1 or INFO(16) =1, add NEQ to the base value. |
|
|
725 C |
|
|
726 C If INFO(12) = 1 (Krylov method), the base value is |
|
|
727 C base = 40 + LENIWP. |
|
|
728 C See PSOL for description of LENIWP. |
|
|
729 C IF INFO(10) = 1 or 3, add NEQ to the base value. |
|
|
730 C If INFO(11) = 1 or INFO(16) = 1, add NEQ to the base value. |
|
|
731 C |
|
|
732 C |
|
|
733 C RPAR, IPAR -- These are arrays of double precision and integer type, |
|
|
734 C respectively, which are available for you to use |
|
|
735 C for communication between your program that calls |
|
|
736 C DDASPK and the RES subroutine (and the JAC and PSOL |
|
|
737 C subroutines). They are not altered by DDASPK. |
|
|
738 C If you do not need RPAR or IPAR, ignore these |
|
|
739 C parameters by treating them as dummy arguments. |
|
|
740 C If you do choose to use them, dimension them in |
|
|
741 C your calling program and in RES (and in JAC and PSOL) |
|
|
742 C as arrays of appropriate length. |
|
|
743 C |
|
|
744 C JAC -- This is the name of a routine that you may supply |
|
|
745 C (optionally) that relates to the Jacobian matrix of the |
|
|
746 C nonlinear system that the code must solve at each T step. |
|
|
747 C The role of JAC (and its call sequence) depends on whether |
|
|
748 C a direct (INFO(12) = 0) or Krylov (INFO(12) = 1) method |
|
|
749 C is selected. |
|
|
750 C |
|
|
751 C **** INFO(12) = 0 (direct methods): |
|
|
752 C If you are letting the code generate partial derivatives |
|
|
753 C numerically (INFO(5) = 0), then JAC can be absent |
|
|
754 C (or perhaps a dummy routine to satisfy the loader). |
|
|
755 C Otherwise you must supply a JAC routine to compute |
|
|
756 C the matrix A = dG/dY + CJ*dG/dYPRIME. It must have |
|
|
757 C the form |
|
|
758 C |
|
|
759 C SUBROUTINE JAC (T, Y, YPRIME, PD, CJ, RPAR, IPAR) |
|
|
760 C |
|
|
761 C The JAC routine must dimension Y, YPRIME, and PD (and RPAR |
|
|
762 C and IPAR if used). CJ is a scalar which is input to JAC. |
|
|
763 C For the given values of T, Y, and YPRIME, the JAC routine |
|
|
764 C must evaluate the nonzero elements of the matrix A, and |
|
|
765 C store these values in the array PD. The elements of PD are |
|
|
766 C set to zero before each call to JAC, so that only nonzero |
|
|
767 C elements need to be defined. |
|
|
768 C The way you store the elements into the PD array depends |
|
|
769 C on the structure of the matrix indicated by INFO(6). |
|
|
770 C *** INFO(6) = 0 (full or dense matrix) *** |
|
|
771 C Give PD a first dimension of NEQ. When you evaluate the |
|
|
772 C nonzero partial derivatives of equation i (i.e. of G(i)) |
|
|
773 C with respect to component j (of Y and YPRIME), you must |
|
|
774 C store the element in PD according to |
|
|
775 C PD(i,j) = dG(i)/dY(j) + CJ*dG(i)/dYPRIME(j). |
|
|
776 C *** INFO(6) = 1 (banded matrix with half-bandwidths ML, MU |
|
|
777 C as described under INFO(6)) *** |
|
|
778 C Give PD a first dimension of 2*ML+MU+1. When you |
|
|
779 C evaluate the nonzero partial derivatives of equation i |
|
|
780 C (i.e. of G(i)) with respect to component j (of Y and |
|
|
781 C YPRIME), you must store the element in PD according to |
|
|
782 C IROW = i - j + ML + MU + 1 |
|
|
783 C PD(IROW,j) = dG(i)/dY(j) + CJ*dG(i)/dYPRIME(j). |
|
|
784 C |
|
|
785 C **** INFO(12) = 1 (Krylov method): |
|
|
786 C If you are not calculating Jacobian data in advance for use |
|
|
787 C in PSOL (INFO(15) = 0), JAC can be absent (or perhaps a |
|
|
788 C dummy routine to satisfy the loader). Otherwise, you may |
|
|
789 C supply a JAC routine to compute and preprocess any parts of |
|
|
790 C of the Jacobian matrix A = dG/dY + CJ*dG/dYPRIME that are |
|
|
791 C involved in the preconditioner matrix P. |
|
|
792 C It is to have the form |
|
|
793 C |
|
|
794 C SUBROUTINE JAC (RES, IRES, NEQ, T, Y, YPRIME, REWT, SAVR, |
|
|
795 C WK, H, CJ, WP, IWP, IER, RPAR, IPAR) |
|
|
796 C |
|
|
797 C The JAC routine must dimension Y, YPRIME, REWT, SAVR, WK, |
|
|
798 C and (if used) WP, IWP, RPAR, and IPAR. |
|
|
799 C The Y, YPRIME, and SAVR arrays contain the current values |
|
|
800 C of Y, YPRIME, and the residual G, respectively. |
|
|
801 C The array WK is work space of length NEQ. |
|
|
802 C H is the step size. CJ is a scalar, input to JAC, that is |
|
|
803 C normally proportional to 1/H. REWT is an array of |
|
|
804 C reciprocal error weights, 1/EWT(i), where EWT(i) is |
|
|
805 C RTOL*abs(Y(i)) + ATOL (unless you supplied routine DDAWTS |
|
|
806 C instead), for use in JAC if needed. For example, if JAC |
|
|
807 C computes difference quotient approximations to partial |
|
|
808 C derivatives, the REWT array may be useful in setting the |
|
|
809 C increments used. The JAC routine should do any |
|
|
810 C factorization operations called for, in preparation for |
|
|
811 C solving linear systems in PSOL. The matrix P should |
|
|
812 C be an approximation to the Jacobian, |
|
|
813 C A = dG/dY + CJ*dG/dYPRIME. |
|
|
814 C |
|
|
815 C WP and IWP are real and integer work arrays which you may |
|
|
816 C use for communication between your JAC routine and your |
|
|
817 C PSOL routine. These may be used to store elements of the |
|
|
818 C preconditioner P, or related matrix data (such as factored |
|
|
819 C forms). They are not altered by DDASPK. |
|
|
820 C If you do not need WP or IWP, ignore these parameters by |
|
|
821 C treating them as dummy arguments. If you do use them, |
|
|
822 C dimension them appropriately in your JAC and PSOL routines. |
|
|
823 C See the PSOL description for instructions on setting |
|
|
824 C the lengths of WP and IWP. |
|
|
825 C |
|
|
826 C On return, JAC should set the error flag IER as follows.. |
|
|
827 C IER = 0 if JAC was successful, |
|
|
828 C IER .ne. 0 if JAC was unsuccessful (e.g. if Y or YPRIME |
|
|
829 C was illegal, or a singular matrix is found). |
|
|
830 C (If IER .ne. 0, a smaller stepsize will be tried.) |
|
|
831 C IER = 0 on entry to JAC, so need be reset only on a failure. |
|
|
832 C If RES is used within JAC, then a nonzero value of IRES will |
|
|
833 C override any nonzero value of IER (see the RES description). |
|
|
834 C |
|
|
835 C Regardless of the method type, subroutine JAC must not |
|
|
836 C alter T, Y(*), YPRIME(*), H, CJ, or REWT(*). |
|
|
837 C You must declare the name JAC in an EXTERNAL statement in |
|
|
838 C your program that calls DDASPK. |
|
|
839 C |
|
|
840 C PSOL -- This is the name of a routine you must supply if you have |
|
|
841 C selected a Krylov method (INFO(12) = 1) with preconditioning. |
|
|
842 C In the direct case (INFO(12) = 0), PSOL can be absent |
|
|
843 C (a dummy routine may have to be supplied to satisfy the |
|
|
844 C loader). Otherwise, you must provide a PSOL routine to |
|
|
845 C solve linear systems arising from preconditioning. |
|
|
846 C When supplied with INFO(12) = 1, the PSOL routine is to |
|
|
847 C have the form |
|
|
848 C |
|
|
849 C SUBROUTINE PSOL (NEQ, T, Y, YPRIME, SAVR, WK, CJ, WGHT, |
|
|
850 C WP, IWP, B, EPLIN, IER, RPAR, IPAR) |
|
|
851 C |
|
|
852 C The PSOL routine must solve linear systems of the form |
|
|
853 C P*x = b where P is the left preconditioner matrix. |
|
|
854 C |
|
|
855 C The right-hand side vector b is in the B array on input, and |
|
|
856 C PSOL must return the solution vector x in B. |
|
|
857 C The Y, YPRIME, and SAVR arrays contain the current values |
|
|
858 C of Y, YPRIME, and the residual G, respectively. |
|
|
859 C |
|
|
860 C Work space required by JAC and/or PSOL, and space for data to |
|
|
861 C be communicated from JAC to PSOL is made available in the form |
|
|
862 C of arrays WP and IWP, which are parts of the RWORK and IWORK |
|
|
863 C arrays, respectively. The lengths of these real and integer |
|
|
864 C work spaces WP and IWP must be supplied in LENWP and LENIWP, |
|
|
865 C respectively, as follows.. |
|
|
866 C IWORK(27) = LENWP = length of real work space WP |
|
|
867 C IWORK(28) = LENIWP = length of integer work space IWP. |
|
|
868 C |
|
|
869 C WK is a work array of length NEQ for use by PSOL. |
|
|
870 C CJ is a scalar, input to PSOL, that is normally proportional |
|
|
871 C to 1/H (H = stepsize). If the old value of CJ |
|
|
872 C (at the time of the last JAC call) is needed, it must have |
|
|
873 C been saved by JAC in WP. |
|
|
874 C |
|
|
875 C WGHT is an array of weights, to be used if PSOL uses an |
|
|
876 C iterative method and performs a convergence test. (In terms |
|
|
877 C of the argument REWT to JAC, WGHT is REWT/sqrt(NEQ).) |
|
|
878 C If PSOL uses an iterative method, it should use EPLIN |
|
|
879 C (a heuristic parameter) as the bound on the weighted norm of |
|
|
880 C the residual for the computed solution. Specifically, the |
|
|
881 C residual vector R should satisfy |
|
|
882 C SQRT (SUM ( (R(i)*WGHT(i))**2 ) ) .le. EPLIN |
|
|
883 C |
|
|
884 C PSOL must not alter NEQ, T, Y, YPRIME, SAVR, CJ, WGHT, EPLIN. |
|
|
885 C |
|
|
886 C On return, PSOL should set the error flag IER as follows.. |
|
|
887 C IER = 0 if PSOL was successful, |
|
|
888 C IER .lt. 0 if an unrecoverable error occurred, meaning |
|
|
889 C control will be passed to the calling routine, |
|
|
890 C IER .gt. 0 if a recoverable error occurred, meaning that |
|
|
891 C the step will be retried with the same step size |
|
|
892 C but with a call to JAC to update necessary data, |
|
|
893 C unless the Jacobian data is current, in which case |
|
|
894 C the step will be retried with a smaller step size. |
|
|
895 C IER = 0 on entry to PSOL so need be reset only on a failure. |
|
|
896 C |
|
|
897 C You must declare the name PSOL in an EXTERNAL statement in |
|
|
898 C your program that calls DDASPK. |
|
|
899 C |
|
|
900 C |
|
|
901 C OPTIONALLY REPLACEABLE SUBROUTINE: |
|
|
902 C |
|
|
903 C DDASPK uses a weighted root-mean-square norm to measure the |
|
|
904 C size of various error vectors. The weights used in this norm |
|
|
905 C are set in the following subroutine: |
|
|
906 C |
|
|
907 C SUBROUTINE DDAWTS (NEQ, IWT, RTOL, ATOL, Y, EWT, RPAR, IPAR) |
|
|
908 C DIMENSION RTOL(*), ATOL(*), Y(*), EWT(*), RPAR(*), IPAR(*) |
|
|
909 C |
|
|
910 C A DDAWTS routine has been included with DDASPK which sets the |
|
|
911 C weights according to |
|
|
912 C EWT(I) = RTOL*ABS(Y(I)) + ATOL |
|
|
913 C in the case of scalar tolerances (IWT = 0) or |
|
|
914 C EWT(I) = RTOL(I)*ABS(Y(I)) + ATOL(I) |
|
|
915 C in the case of array tolerances (IWT = 1). (IWT is INFO(2).) |
|
|
916 C In some special cases, it may be appropriate for you to define |
|
|
917 C your own error weights by writing a subroutine DDAWTS to be |
|
|
918 C called instead of the version supplied. However, this should |
|
|
919 C be attempted only after careful thought and consideration. |
|
|
920 C If you supply this routine, you may use the tolerances and Y |
|
|
921 C as appropriate, but do not overwrite these variables. You |
|
|
922 C may also use RPAR and IPAR to communicate data as appropriate. |
|
|
923 C ***Note: Aside from the values of the weights, the choice of |
|
|
924 C norm used in DDASPK (weighted root-mean-square) is not subject |
|
|
925 C to replacement by the user. In this respect, DDASPK is not |
|
|
926 C downward-compatible with the original DDASSL solver (in which |
|
|
927 C the norm routine was optionally user-replaceable). |
|
|
928 C |
|
|
929 C |
|
|
930 C------OUTPUT - AFTER ANY RETURN FROM DDASPK---------------------------- |
|
|
931 C |
|
|
932 C The principal aim of the code is to return a computed solution at |
|
|
933 C T = TOUT, although it is also possible to obtain intermediate |
|
|
934 C results along the way. To find out whether the code achieved its |
|
|
935 C goal or if the integration process was interrupted before the task |
|
|
936 C was completed, you must check the IDID parameter. |
|
|
937 C |
|
|
938 C |
|
|
939 C T -- The output value of T is the point to which the solution |
|
|
940 C was successfully advanced. |
|
|
941 C |
|
|
942 C Y(*) -- contains the computed solution approximation at T. |
|
|
943 C |
|
|
944 C YPRIME(*) -- contains the computed derivative approximation at T. |
|
|
945 C |
|
|
946 C IDID -- reports what the code did, described as follows: |
|
|
947 C |
|
|
948 C *** TASK COMPLETED *** |
|
|
949 C Reported by positive values of IDID |
|
|
950 C |
|
|
951 C IDID = 1 -- a step was successfully taken in the |
|
|
952 C intermediate-output mode. The code has not |
|
|
953 C yet reached TOUT. |
|
|
954 C |
|
|
955 C IDID = 2 -- the integration to TSTOP was successfully |
|
|
956 C completed (T = TSTOP) by stepping exactly to TSTOP. |
|
|
957 C |
|
|
958 C IDID = 3 -- the integration to TOUT was successfully |
|
|
959 C completed (T = TOUT) by stepping past TOUT. |
|
|
960 C Y(*) and YPRIME(*) are obtained by interpolation. |
|
|
961 C |
|
|
962 C IDID = 4 -- the initial condition calculation, with |
|
|
963 C INFO(11) > 0, was successful, and INFO(14) = 1. |
|
|
964 C No integration steps were taken, and the solution |
|
|
965 C is not considered to have been started. |
|
|
966 C |
|
|
967 C *** TASK INTERRUPTED *** |
|
|
968 C Reported by negative values of IDID |
|
|
969 C |
|
|
970 C IDID = -1 -- a large amount of work has been expended |
|
|
971 C (about 500 steps). |
|
|
972 C |
|
|
973 C IDID = -2 -- the error tolerances are too stringent. |
|
|
974 C |
|
|
975 C IDID = -3 -- the local error test cannot be satisfied |
|
|
976 C because you specified a zero component in ATOL |
|
|
977 C and the corresponding computed solution component |
|
|
978 C is zero. Thus, a pure relative error test is |
|
|
979 C impossible for this component. |
|
|
980 C |
|
|
981 C IDID = -5 -- there were repeated failures in the evaluation |
|
|
982 C or processing of the preconditioner (in JAC). |
|
|
983 C |
|
|
984 C IDID = -6 -- DDASPK had repeated error test failures on the |
|
|
985 C last attempted step. |
|
|
986 C |
|
|
987 C IDID = -7 -- the nonlinear system solver in the time integration |
|
|
988 C could not converge. |
|
|
989 C |
|
|
990 C IDID = -8 -- the matrix of partial derivatives appears |
|
|
991 C to be singular (direct method). |
|
|
992 C |
|
|
993 C IDID = -9 -- the nonlinear system solver in the time integration |
|
|
994 C failed to achieve convergence, and there were repeated |
|
|
995 C error test failures in this step. |
|
|
996 C |
|
|
997 C IDID =-10 -- the nonlinear system solver in the time integration |
|
|
998 C failed to achieve convergence because IRES was equal |
|
|
999 C to -1. |
|
|
1000 C |
|
|
1001 C IDID =-11 -- IRES = -2 was encountered and control is |
|
|
1002 C being returned to the calling program. |
|
|
1003 C |
|
|
1004 C IDID =-12 -- DDASPK failed to compute the initial Y, YPRIME. |
|
|
1005 C |
|
|
1006 C IDID =-13 -- unrecoverable error encountered inside user's |
|
|
1007 C PSOL routine, and control is being returned to |
|
|
1008 C the calling program. |
|
|
1009 C |
|
|
1010 C IDID =-14 -- the Krylov linear system solver could not |
|
|
1011 C achieve convergence. |
|
|
1012 C |
|
|
1013 C IDID =-15,..,-32 -- Not applicable for this code. |
|
|
1014 C |
|
|
1015 C *** TASK TERMINATED *** |
|
|
1016 C reported by the value of IDID=-33 |
|
|
1017 C |
|
|
1018 C IDID = -33 -- the code has encountered trouble from which |
|
|
1019 C it cannot recover. A message is printed |
|
|
1020 C explaining the trouble and control is returned |
|
|
1021 C to the calling program. For example, this occurs |
|
|
1022 C when invalid input is detected. |
|
|
1023 C |
|
|
1024 C RTOL, ATOL -- these quantities remain unchanged except when |
|
|
1025 C IDID = -2. In this case, the error tolerances have been |
|
|
1026 C increased by the code to values which are estimated to |
|
|
1027 C be appropriate for continuing the integration. However, |
|
|
1028 C the reported solution at T was obtained using the input |
|
|
1029 C values of RTOL and ATOL. |
|
|
1030 C |
|
|
1031 C RWORK, IWORK -- contain information which is usually of no interest |
|
|
1032 C to the user but necessary for subsequent calls. |
|
|
1033 C However, you may be interested in the performance data |
|
|
1034 C listed below. These quantities are accessed in RWORK |
|
|
1035 C and IWORK but have internal mnemonic names, as follows.. |
|
|
1036 C |
|
|
1037 C RWORK(3)--contains H, the step size h to be attempted |
|
|
1038 C on the next step. |
|
|
1039 C |
|
|
1040 C RWORK(4)--contains TN, the current value of the |
|
|
1041 C independent variable, i.e. the farthest point |
|
|
1042 C integration has reached. This will differ |
|
|
1043 C from T if interpolation has been performed |
|
|
1044 C (IDID = 3). |
|
|
1045 C |
|
|
1046 C RWORK(7)--contains HOLD, the stepsize used on the last |
|
|
1047 C successful step. If INFO(11) = INFO(14) = 1, |
|
|
1048 C this contains the value of H used in the |
|
|
1049 C initial condition calculation. |
|
|
1050 C |
|
|
1051 C IWORK(7)--contains K, the order of the method to be |
|
|
1052 C attempted on the next step. |
|
|
1053 C |
|
|
1054 C IWORK(8)--contains KOLD, the order of the method used |
|
|
1055 C on the last step. |
|
|
1056 C |
|
|
1057 C IWORK(11)--contains NST, the number of steps (in T) |
|
|
1058 C taken so far. |
|
|
1059 C |
|
|
1060 C IWORK(12)--contains NRE, the number of calls to RES |
|
|
1061 C so far. |
|
|
1062 C |
|
|
1063 C IWORK(13)--contains NJE, the number of calls to JAC so |
|
|
1064 C far (Jacobian or preconditioner evaluations). |
|
|
1065 C |
|
|
1066 C IWORK(14)--contains NETF, the total number of error test |
|
|
1067 C failures so far. |
|
|
1068 C |
|
|
1069 C IWORK(15)--contains NCFN, the total number of nonlinear |
|
|
1070 C convergence failures so far (includes counts |
|
|
1071 C of singular iteration matrix or singular |
|
|
1072 C preconditioners). |
|
|
1073 C |
|
|
1074 C IWORK(16)--contains NCFL, the number of convergence |
|
|
1075 C failures of the linear iteration so far. |
|
|
1076 C |
|
|
1077 C IWORK(17)--contains LENIW, the length of IWORK actually |
|
|
1078 C required. This is defined on normal returns |
|
|
1079 C and on an illegal input return for |
|
|
1080 C insufficient storage. |
|
|
1081 C |
|
|
1082 C IWORK(18)--contains LENRW, the length of RWORK actually |
|
|
1083 C required. This is defined on normal returns |
|
|
1084 C and on an illegal input return for |
|
|
1085 C insufficient storage. |
|
|
1086 C |
|
|
1087 C IWORK(19)--contains NNI, the total number of nonlinear |
|
|
1088 C iterations so far (each of which calls a |
|
|
1089 C linear solver). |
|
|
1090 C |
|
|
1091 C IWORK(20)--contains NLI, the total number of linear |
|
|
1092 C (Krylov) iterations so far. |
|
|
1093 C |
|
|
1094 C IWORK(21)--contains NPS, the number of PSOL calls so |
|
|
1095 C far, for preconditioning solve operations or |
|
|
1096 C for solutions with the user-supplied method. |
|
|
1097 C |
|
|
1098 C Note: The various counters in IWORK do not include |
|
|
1099 C counts during a call made with INFO(11) > 0 and |
|
|
1100 C INFO(14) = 1. |
|
|
1101 C |
|
|
1102 C |
|
|
1103 C------INPUT - WHAT TO DO TO CONTINUE THE INTEGRATION ----------------- |
|
|
1104 C (CALLS AFTER THE FIRST) |
|
|
1105 C |
|
|
1106 C This code is organized so that subsequent calls to continue the |
|
|
1107 C integration involve little (if any) additional effort on your |
|
|
1108 C part. You must monitor the IDID parameter in order to determine |
|
|
1109 C what to do next. |
|
|
1110 C |
|
|
1111 C Recalling that the principal task of the code is to integrate |
|
|
1112 C from T to TOUT (the interval mode), usually all you will need |
|
|
1113 C to do is specify a new TOUT upon reaching the current TOUT. |
|
|
1114 C |
|
|
1115 C Do not alter any quantity not specifically permitted below. In |
|
|
1116 C particular do not alter NEQ, T, Y(*), YPRIME(*), RWORK(*), |
|
|
1117 C IWORK(*), or the differential equation in subroutine RES. Any |
|
|
1118 C such alteration constitutes a new problem and must be treated |
|
|
1119 C as such, i.e. you must start afresh. |
|
|
1120 C |
|
|
1121 C You cannot change from array to scalar error control or vice |
|
|
1122 C versa (INFO(2)), but you can change the size of the entries of |
|
|
1123 C RTOL or ATOL. Increasing a tolerance makes the equation easier |
|
|
1124 C to integrate. Decreasing a tolerance will make the equation |
|
|
1125 C harder to integrate and should generally be avoided. |
|
|
1126 C |
|
|
1127 C You can switch from the intermediate-output mode to the |
|
|
1128 C interval mode (INFO(3)) or vice versa at any time. |
|
|
1129 C |
|
|
1130 C If it has been necessary to prevent the integration from going |
|
|
1131 C past a point TSTOP (INFO(4), RWORK(1)), keep in mind that the |
|
|
1132 C code will not integrate to any TOUT beyond the currently |
|
|
1133 C specified TSTOP. Once TSTOP has been reached, you must change |
|
|
1134 C the value of TSTOP or set INFO(4) = 0. You may change INFO(4) |
|
|
1135 C or TSTOP at any time but you must supply the value of TSTOP in |
|
|
1136 C RWORK(1) whenever you set INFO(4) = 1. |
|
|
1137 C |
|
|
1138 C Do not change INFO(5), INFO(6), INFO(12-17) or their associated |
|
|
1139 C IWORK/RWORK locations unless you are going to restart the code. |
|
|
1140 C |
|
|
1141 C *** FOLLOWING A COMPLETED TASK *** |
|
|
1142 C |
|
|
1143 C If.. |
|
|
1144 C IDID = 1, call the code again to continue the integration |
|
|
1145 C another step in the direction of TOUT. |
|
|
1146 C |
|
|
1147 C IDID = 2 or 3, define a new TOUT and call the code again. |
|
|
1148 C TOUT must be different from T. You cannot change |
|
|
1149 C the direction of integration without restarting. |
|
|
1150 C |
|
|
1151 C IDID = 4, reset INFO(11) = 0 and call the code again to begin |
|
|
1152 C the integration. (If you leave INFO(11) > 0 and |
|
|
1153 C INFO(14) = 1, you may generate an infinite loop.) |
|
|
1154 C In this situation, the next call to DASPK is |
|
|
1155 C considered to be the first call for the problem, |
|
|
1156 C in that all initializations are done. |
|
|
1157 C |
|
|
1158 C *** FOLLOWING AN INTERRUPTED TASK *** |
|
|
1159 C |
|
|
1160 C To show the code that you realize the task was interrupted and |
|
|
1161 C that you want to continue, you must take appropriate action and |
|
|
1162 C set INFO(1) = 1. |
|
|
1163 C |
|
|
1164 C If.. |
|
|
1165 C IDID = -1, the code has taken about 500 steps. If you want to |
|
|
1166 C continue, set INFO(1) = 1 and call the code again. |
|
|
1167 C An additional 500 steps will be allowed. |
|
|
1168 C |
|
|
1169 C |
|
|
1170 C IDID = -2, the error tolerances RTOL, ATOL have been increased |
|
|
1171 C to values the code estimates appropriate for |
|
|
1172 C continuing. You may want to change them yourself. |
|
|
1173 C If you are sure you want to continue with relaxed |
|
|
1174 C error tolerances, set INFO(1) = 1 and call the code |
|
|
1175 C again. |
|
|
1176 C |
|
|
1177 C IDID = -3, a solution component is zero and you set the |
|
|
1178 C corresponding component of ATOL to zero. If you |
|
|
1179 C are sure you want to continue, you must first alter |
|
|
1180 C the error criterion to use positive values of ATOL |
|
|
1181 C for those components corresponding to zero solution |
|
|
1182 C components, then set INFO(1) = 1 and call the code |
|
|
1183 C again. |
|
|
1184 C |
|
|
1185 C IDID = -4 --- cannot occur with this code. |
|
|
1186 C |
|
|
1187 C IDID = -5, your JAC routine failed with the Krylov method. Check |
|
|
1188 C for errors in JAC and restart the integration. |
|
|
1189 C |
|
|
1190 C IDID = -6, repeated error test failures occurred on the last |
|
|
1191 C attempted step in DDASPK. A singularity in the |
|
|
1192 C solution may be present. If you are absolutely |
|
|
1193 C certain you want to continue, you should restart |
|
|
1194 C the integration. (Provide initial values of Y and |
|
|
1195 C YPRIME which are consistent.) |
|
|
1196 C |
|
|
1197 C IDID = -7, repeated convergence test failures occurred on the last |
|
|
1198 C attempted step in DDASPK. An inaccurate or ill- |
|
|
1199 C conditioned Jacobian or preconditioner may be the |
|
|
1200 C problem. If you are absolutely certain you want |
|
|
1201 C to continue, you should restart the integration. |
|
|
1202 C |
|
|
1203 C |
|
|
1204 C IDID = -8, the matrix of partial derivatives is singular, with |
|
|
1205 C the use of direct methods. Some of your equations |
|
|
1206 C may be redundant. DDASPK cannot solve the problem |
|
|
1207 C as stated. It is possible that the redundant |
|
|
1208 C equations could be removed, and then DDASPK could |
|
|
1209 C solve the problem. It is also possible that a |
|
|
1210 C solution to your problem either does not exist |
|
|
1211 C or is not unique. |
|
|
1212 C |
|
|
1213 C IDID = -9, DDASPK had multiple convergence test failures, preceded |
|
|
1214 C by multiple error test failures, on the last |
|
|
1215 C attempted step. It is possible that your problem is |
|
|
1216 C ill-posed and cannot be solved using this code. Or, |
|
|
1217 C there may be a discontinuity or a singularity in the |
|
|
1218 C solution. If you are absolutely certain you want to |
|
|
1219 C continue, you should restart the integration. |
|
|
1220 C |
|
|
1221 C IDID = -10, DDASPK had multiple convergence test failures |
|
|
1222 C because IRES was equal to -1. If you are |
|
|
1223 C absolutely certain you want to continue, you |
|
|
1224 C should restart the integration. |
|
|
1225 C |
|
|
1226 C IDID = -11, there was an unrecoverable error (IRES = -2) from RES |
|
|
1227 C inside the nonlinear system solver. Determine the |
|
|
1228 C cause before trying again. |
|
|
1229 C |
|
|
1230 C IDID = -12, DDASPK failed to compute the initial Y and YPRIME |
|
|
1231 C vectors. This could happen because the initial |
|
|
1232 C approximation to Y or YPRIME was not very good, or |
|
|
1233 C because no consistent values of these vectors exist. |
|
|
1234 C The problem could also be caused by an inaccurate or |
|
|
1235 C singular iteration matrix, or a poor preconditioner. |
|
|
1236 C |
|
|
1237 C IDID = -13, there was an unrecoverable error encountered inside |
|
|
1238 C your PSOL routine. Determine the cause before |
|
|
1239 C trying again. |
|
|
1240 C |
|
|
1241 C IDID = -14, the Krylov linear system solver failed to achieve |
|
|
1242 C convergence. This may be due to ill-conditioning |
|
|
1243 C in the iteration matrix, or a singularity in the |
|
|
1244 C preconditioner (if one is being used). |
|
|
1245 C Another possibility is that there is a better |
|
|
1246 C choice of Krylov parameters (see INFO(13)). |
|
|
1247 C Possibly the failure is caused by redundant equations |
|
|
1248 C in the system, or by inconsistent equations. |
|
|
1249 C In that case, reformulate the system to make it |
|
|
1250 C consistent and non-redundant. |
|
|
1251 C |
|
|
1252 C IDID = -15,..,-32 --- Cannot occur with this code. |
|
|
1253 C |
|
|
1254 C *** FOLLOWING A TERMINATED TASK *** |
|
|
1255 C |
|
|
1256 C If IDID = -33, you cannot continue the solution of this problem. |
|
|
1257 C An attempt to do so will result in your run being |
|
|
1258 C terminated. |
|
|
1259 C |
|
|
1260 C --------------------------------------------------------------------- |
|
|
1261 C |
|
|
1262 C***REFERENCES |
|
|
1263 C 1. L. R. Petzold, A Description of DASSL: A Differential/Algebraic |
|
|
1264 C System Solver, in Scientific Computing, R. S. Stepleman et al. |
|
|
1265 C (Eds.), North-Holland, Amsterdam, 1983, pp. 65-68. |
|
|
1266 C 2. K. E. Brenan, S. L. Campbell, and L. R. Petzold, Numerical |
|
|
1267 C Solution of Initial-Value Problems in Differential-Algebraic |
|
|
1268 C Equations, Elsevier, New York, 1989. |
|
|
1269 C 3. P. N. Brown and A. C. Hindmarsh, Reduced Storage Matrix Methods |
|
|
1270 C in Stiff ODE Systems, J. Applied Mathematics and Computation, |
|
|
1271 C 31 (1989), pp. 40-91. |
|
|
1272 C 4. P. N. Brown, A. C. Hindmarsh, and L. R. Petzold, Using Krylov |
|
|
1273 C Methods in the Solution of Large-Scale Differential-Algebraic |
|
|
1274 C Systems, SIAM J. Sci. Comp., 15 (1994), pp. 1467-1488. |
|
|
1275 C 5. P. N. Brown, A. C. Hindmarsh, and L. R. Petzold, Consistent |
|
|
1276 C Initial Condition Calculation for Differential-Algebraic |
|
|
1277 C Systems, LLNL Report UCRL-JC-122175, August 1995; submitted to |
|
|
1278 C SIAM J. Sci. Comp. |
|
|
1279 C |
|
|
1280 C***ROUTINES CALLED |
|
|
1281 C |
|
|
1282 C The following are all the subordinate routines used by DDASPK. |
|
|
1283 C |
|
|
1284 C DDASIC computes consistent initial conditions. |
|
|
1285 C DYYPNW updates Y and YPRIME in linesearch for initial condition |
|
|
1286 C calculation. |
|
|
1287 C DDSTP carries out one step of the integration. |
|
|
1288 C DCNSTR/DCNST0 check the current solution for constraint violations. |
|
|
1289 C DDAWTS sets error weight quantities. |
|
|
1290 C DINVWT tests and inverts the error weights. |
|
|
1291 C DDATRP performs interpolation to get an output solution. |
|
|
1292 C DDWNRM computes the weighted root-mean-square norm of a vector. |
|
|
1293 C D1MACH provides the unit roundoff of the computer. |
|
|
1294 C XERRWD/XSETF/XSETUN/IXSAV is a package to handle error messages. |
|
|
1295 C DDASID nonlinear equation driver to initialize Y and YPRIME using |
|
|
1296 C direct linear system solver methods. Interfaces to Newton |
|
|
1297 C solver (direct case). |
|
|
1298 C DNSID solves the nonlinear system for unknown initial values by |
|
|
1299 C modified Newton iteration and direct linear system methods. |
|
|
1300 C DLINSD carries out linesearch algorithm for initial condition |
|
|
1301 C calculation (direct case). |
|
|
1302 C DFNRMD calculates weighted norm of preconditioned residual in |
|
|
1303 C initial condition calculation (direct case). |
|
|
1304 C DNEDD nonlinear equation driver for direct linear system solver |
|
|
1305 C methods. Interfaces to Newton solver (direct case). |
|
|
1306 C DMATD assembles the iteration matrix (direct case). |
|
|
1307 C DNSD solves the associated nonlinear system by modified |
|
|
1308 C Newton iteration and direct linear system methods. |
|
|
1309 C DSLVD interfaces to linear system solver (direct case). |
|
|
1310 C DDASIK nonlinear equation driver to initialize Y and YPRIME using |
|
|
1311 C Krylov iterative linear system methods. Interfaces to |
|
|
1312 C Newton solver (Krylov case). |
|
|
1313 C DNSIK solves the nonlinear system for unknown initial values by |
|
|
1314 C Newton iteration and Krylov iterative linear system methods. |
|
|
1315 C DLINSK carries out linesearch algorithm for initial condition |
|
|
1316 C calculation (Krylov case). |
|
|
1317 C DFNRMK calculates weighted norm of preconditioned residual in |
|
|
1318 C initial condition calculation (Krylov case). |
|
|
1319 C DNEDK nonlinear equation driver for iterative linear system solver |
|
|
1320 C methods. Interfaces to Newton solver (Krylov case). |
|
|
1321 C DNSK solves the associated nonlinear system by Inexact Newton |
|
|
1322 C iteration and (linear) Krylov iteration. |
|
|
1323 C DSLVK interfaces to linear system solver (Krylov case). |
|
|
1324 C DSPIGM solves a linear system by SPIGMR algorithm. |
|
|
1325 C DATV computes matrix-vector product in Krylov algorithm. |
|
|
1326 C DORTH performs orthogonalization of Krylov basis vectors. |
|
|
1327 C DHEQR performs QR factorization of Hessenberg matrix. |
|
|
1328 C DHELS finds least-squares solution of Hessenberg linear system. |
|
4329
|
1329 C DGETRF, DGETRS, DGBTRF, DGBTRS are LAPACK routines for solving |
|
3911
|
1330 C linear systems (dense or band direct methods). |
|
|
1331 C DAXPY, DCOPY, DDOT, DNRM2, DSCAL are Basic Linear Algebra (BLAS) |
|
|
1332 C routines. |
|
|
1333 C |
|
|
1334 C The routines called directly by DDASPK are: |
|
|
1335 C DCNST0, DDAWTS, DINVWT, D1MACH, DDWNRM, DDASIC, DDATRP, DDSTP, |
|
|
1336 C XERRWD |
|
|
1337 C |
|
|
1338 C***END PROLOGUE DDASPK |
|
|
1339 C |
|
|
1340 C |
|
|
1341 IMPLICIT DOUBLE PRECISION(A-H,O-Z) |
|
|
1342 LOGICAL DONE, LAVL, LCFN, LCFL, LWARN |
|
|
1343 DIMENSION Y(*),YPRIME(*) |
|
|
1344 DIMENSION INFO(20) |
|
|
1345 DIMENSION RWORK(LRW),IWORK(LIW) |
|
|
1346 DIMENSION RTOL(*),ATOL(*) |
|
|
1347 DIMENSION RPAR(*),IPAR(*) |
|
|
1348 CHARACTER MSG*80 |
|
|
1349 EXTERNAL RES, JAC, PSOL, DDASID, DDASIK, DNEDD, DNEDK |
|
|
1350 C |
|
|
1351 C Set pointers into IWORK. |
|
|
1352 C |
|
|
1353 PARAMETER (LML=1, LMU=2, LMTYPE=4, |
|
|
1354 * LIWM=1, LMXORD=3, LJCALC=5, LPHASE=6, LK=7, LKOLD=8, |
|
|
1355 * LNS=9, LNSTL=10, LNST=11, LNRE=12, LNJE=13, LETF=14, LNCFN=15, |
|
|
1356 * LNCFL=16, LNIW=17, LNRW=18, LNNI=19, LNLI=20, LNPS=21, |
|
|
1357 * LNPD=22, LMITER=23, LMAXL=24, LKMP=25, LNRMAX=26, LLNWP=27, |
|
|
1358 * LLNIWP=28, LLOCWP=29, LLCIWP=30, LKPRIN=31, |
|
|
1359 * LMXNIT=32, LMXNJ=33, LMXNH=34, LLSOFF=35, LICNS=41) |
|
|
1360 C |
|
|
1361 C Set pointers into RWORK. |
|
|
1362 C |
|
|
1363 PARAMETER (LTSTOP=1, LHMAX=2, LH=3, LTN=4, LCJ=5, LCJOLD=6, |
|
|
1364 * LHOLD=7, LS=8, LROUND=9, LEPLI=10, LSQRN=11, LRSQRN=12, |
|
|
1365 * LEPCON=13, LSTOL=14, LEPIN=15, |
|
|
1366 * LALPHA=21, LBETA=27, LGAMMA=33, LPSI=39, LSIGMA=45, LDELTA=51) |
|
|
1367 C |
|
|
1368 SAVE LID, LENID, NONNEG |
|
|
1369 C |
|
|
1370 C |
|
|
1371 C***FIRST EXECUTABLE STATEMENT DDASPK |
|
|
1372 C |
|
|
1373 C |
|
|
1374 IF(INFO(1).NE.0) GO TO 100 |
|
|
1375 C |
|
|
1376 C----------------------------------------------------------------------- |
|
|
1377 C This block is executed for the initial call only. |
|
|
1378 C It contains checking of inputs and initializations. |
|
|
1379 C----------------------------------------------------------------------- |
|
|
1380 C |
|
|
1381 C First check INFO array to make sure all elements of INFO |
|
|
1382 C Are within the proper range. (INFO(1) is checked later, because |
|
|
1383 C it must be tested on every call.) ITEMP holds the location |
|
|
1384 C within INFO which may be out of range. |
|
|
1385 C |
|
|
1386 DO 10 I=2,9 |
|
|
1387 ITEMP = I |
|
|
1388 IF (INFO(I) .NE. 0 .AND. INFO(I) .NE. 1) GO TO 701 |
|
|
1389 10 CONTINUE |
|
|
1390 ITEMP = 10 |
|
|
1391 IF(INFO(10).LT.0 .OR. INFO(10).GT.3) GO TO 701 |
|
|
1392 ITEMP = 11 |
|
|
1393 IF(INFO(11).LT.0 .OR. INFO(11).GT.2) GO TO 701 |
|
|
1394 DO 15 I=12,17 |
|
|
1395 ITEMP = I |
|
|
1396 IF (INFO(I) .NE. 0 .AND. INFO(I) .NE. 1) GO TO 701 |
|
|
1397 15 CONTINUE |
|
|
1398 ITEMP = 18 |
|
|
1399 IF(INFO(18).LT.0 .OR. INFO(18).GT.2) GO TO 701 |
|
|
1400 |
|
|
1401 C |
|
|
1402 C Check NEQ to see if it is positive. |
|
|
1403 C |
|
|
1404 IF (NEQ .LE. 0) GO TO 702 |
|
|
1405 C |
|
|
1406 C Check and compute maximum order. |
|
|
1407 C |
|
|
1408 MXORD=5 |
|
|
1409 IF (INFO(9) .NE. 0) THEN |
|
|
1410 MXORD=IWORK(LMXORD) |
|
|
1411 IF (MXORD .LT. 1 .OR. MXORD .GT. 5) GO TO 703 |
|
|
1412 ENDIF |
|
|
1413 IWORK(LMXORD)=MXORD |
|
|
1414 C |
|
|
1415 C Set and/or check inputs for constraint checking (INFO(10) .NE. 0). |
|
|
1416 C Set values for ICNFLG, NONNEG, and pointer LID. |
|
|
1417 C |
|
|
1418 ICNFLG = 0 |
|
|
1419 NONNEG = 0 |
|
|
1420 LID = LICNS |
|
|
1421 IF (INFO(10) .EQ. 0) GO TO 20 |
|
|
1422 IF (INFO(10) .EQ. 1) THEN |
|
|
1423 ICNFLG = 1 |
|
|
1424 NONNEG = 0 |
|
|
1425 LID = LICNS + NEQ |
|
|
1426 ELSEIF (INFO(10) .EQ. 2) THEN |
|
|
1427 ICNFLG = 0 |
|
|
1428 NONNEG = 1 |
|
|
1429 ELSE |
|
|
1430 ICNFLG = 1 |
|
|
1431 NONNEG = 1 |
|
|
1432 LID = LICNS + NEQ |
|
|
1433 ENDIF |
|
|
1434 C |
|
|
1435 20 CONTINUE |
|
|
1436 C |
|
|
1437 C Set and/or check inputs for Krylov solver (INFO(12) .NE. 0). |
|
|
1438 C If indicated, set default values for MAXL, KMP, NRMAX, and EPLI. |
|
|
1439 C Otherwise, verify inputs required for iterative solver. |
|
|
1440 C |
|
|
1441 IF (INFO(12) .EQ. 0) GO TO 25 |
|
|
1442 C |
|
|
1443 IWORK(LMITER) = INFO(12) |
|
|
1444 IF (INFO(13) .EQ. 0) THEN |
|
|
1445 IWORK(LMAXL) = MIN(5,NEQ) |
|
|
1446 IWORK(LKMP) = IWORK(LMAXL) |
|
|
1447 IWORK(LNRMAX) = 5 |
|
|
1448 RWORK(LEPLI) = 0.05D0 |
|
|
1449 ELSE |
|
|
1450 IF(IWORK(LMAXL) .LT. 1 .OR. IWORK(LMAXL) .GT. NEQ) GO TO 720 |
|
|
1451 IF(IWORK(LKMP) .LT. 1 .OR. IWORK(LKMP) .GT. IWORK(LMAXL)) |
|
|
1452 1 GO TO 721 |
|
|
1453 IF(IWORK(LNRMAX) .LT. 0) GO TO 722 |
|
|
1454 IF(RWORK(LEPLI).LE.0.0D0 .OR. RWORK(LEPLI).GE.1.0D0)GO TO 723 |
|
|
1455 ENDIF |
|
|
1456 C |
|
|
1457 25 CONTINUE |
|
|
1458 C |
|
|
1459 C Set and/or check controls for the initial condition calculation |
|
|
1460 C (INFO(11) .GT. 0). If indicated, set default values. |
|
|
1461 C Otherwise, verify inputs required for iterative solver. |
|
|
1462 C |
|
|
1463 IF (INFO(11) .EQ. 0) GO TO 30 |
|
|
1464 IF (INFO(17) .EQ. 0) THEN |
|
|
1465 IWORK(LMXNIT) = 5 |
|
|
1466 IF (INFO(12) .GT. 0) IWORK(LMXNIT) = 15 |
|
|
1467 IWORK(LMXNJ) = 6 |
|
|
1468 IF (INFO(12) .GT. 0) IWORK(LMXNJ) = 2 |
|
|
1469 IWORK(LMXNH) = 5 |
|
|
1470 IWORK(LLSOFF) = 0 |
|
|
1471 RWORK(LEPIN) = 0.01D0 |
|
|
1472 ELSE |
|
|
1473 IF (IWORK(LMXNIT) .LE. 0) GO TO 725 |
|
|
1474 IF (IWORK(LMXNJ) .LE. 0) GO TO 725 |
|
|
1475 IF (IWORK(LMXNH) .LE. 0) GO TO 725 |
|
|
1476 LSOFF = IWORK(LLSOFF) |
|
|
1477 IF (LSOFF .LT. 0 .OR. LSOFF .GT. 1) GO TO 725 |
|
|
1478 IF (RWORK(LEPIN) .LE. 0.0D0) GO TO 725 |
|
|
1479 ENDIF |
|
|
1480 C |
|
|
1481 30 CONTINUE |
|
|
1482 C |
|
|
1483 C Below is the computation and checking of the work array lengths |
|
|
1484 C LENIW and LENRW, using direct methods (INFO(12) = 0) or |
|
|
1485 C the Krylov methods (INFO(12) = 1). |
|
|
1486 C |
|
|
1487 LENIC = 0 |
|
|
1488 IF (INFO(10) .EQ. 1 .OR. INFO(10) .EQ. 3) LENIC = NEQ |
|
|
1489 LENID = 0 |
|
|
1490 IF (INFO(11) .EQ. 1 .OR. INFO(16) .EQ. 1) LENID = NEQ |
|
|
1491 IF (INFO(12) .EQ. 0) THEN |
|
|
1492 C |
|
|
1493 C Compute MTYPE, etc. Check ML and MU. |
|
|
1494 C |
|
|
1495 NCPHI = MAX(MXORD + 1, 4) |
|
|
1496 IF(INFO(6).EQ.0) THEN |
|
|
1497 LENPD = NEQ**2 |
|
|
1498 LENRW = 50 + (NCPHI+3)*NEQ + LENPD |
|
|
1499 IF(INFO(5).EQ.0) THEN |
|
|
1500 IWORK(LMTYPE)=2 |
|
|
1501 ELSE |
|
|
1502 IWORK(LMTYPE)=1 |
|
|
1503 ENDIF |
|
|
1504 ELSE |
|
|
1505 IF(IWORK(LML).LT.0.OR.IWORK(LML).GE.NEQ)GO TO 717 |
|
|
1506 IF(IWORK(LMU).LT.0.OR.IWORK(LMU).GE.NEQ)GO TO 718 |
|
|
1507 LENPD=(2*IWORK(LML)+IWORK(LMU)+1)*NEQ |
|
|
1508 IF(INFO(5).EQ.0) THEN |
|
|
1509 IWORK(LMTYPE)=5 |
|
|
1510 MBAND=IWORK(LML)+IWORK(LMU)+1 |
|
|
1511 MSAVE=(NEQ/MBAND)+1 |
|
|
1512 LENRW = 50 + (NCPHI+3)*NEQ + LENPD + 2*MSAVE |
|
|
1513 ELSE |
|
|
1514 IWORK(LMTYPE)=4 |
|
|
1515 LENRW = 50 + (NCPHI+3)*NEQ + LENPD |
|
|
1516 ENDIF |
|
|
1517 ENDIF |
|
|
1518 C |
|
|
1519 C Compute LENIW, LENWP, LENIWP. |
|
|
1520 C |
|
|
1521 LENIW = 40 + LENIC + LENID + NEQ |
|
|
1522 LENWP = 0 |
|
|
1523 LENIWP = 0 |
|
|
1524 C |
|
|
1525 ELSE IF (INFO(12) .EQ. 1) THEN |
|
|
1526 MAXL = IWORK(LMAXL) |
|
|
1527 LENWP = IWORK(LLNWP) |
|
|
1528 LENIWP = IWORK(LLNIWP) |
|
|
1529 LENPD = (MAXL+3+MIN0(1,MAXL-IWORK(LKMP)))*NEQ |
|
|
1530 1 + (MAXL+3)*MAXL + 1 + LENWP |
|
|
1531 LENRW = 50 + (IWORK(LMXORD)+5)*NEQ + LENPD |
|
|
1532 LENIW = 40 + LENIC + LENID + LENIWP |
|
|
1533 C |
|
|
1534 ENDIF |
|
|
1535 IF(INFO(16) .NE. 0) LENRW = LENRW + NEQ |
|
|
1536 C |
|
|
1537 C Check lengths of RWORK and IWORK. |
|
|
1538 C |
|
|
1539 IWORK(LNIW)=LENIW |
|
|
1540 IWORK(LNRW)=LENRW |
|
|
1541 IWORK(LNPD)=LENPD |
|
|
1542 IWORK(LLOCWP) = LENPD-LENWP+1 |
|
|
1543 IF(LRW.LT.LENRW)GO TO 704 |
|
|
1544 IF(LIW.LT.LENIW)GO TO 705 |
|
|
1545 C |
|
|
1546 C Check ICNSTR for legality. |
|
|
1547 C |
|
|
1548 IF (LENIC .GT. 0) THEN |
|
|
1549 DO 40 I = 1,NEQ |
|
|
1550 ICI = IWORK(LICNS-1+I) |
|
|
1551 IF (ICI .LT. -2 .OR. ICI .GT. 2) GO TO 726 |
|
|
1552 40 CONTINUE |
|
|
1553 ENDIF |
|
|
1554 C |
|
|
1555 C Check Y for consistency with constraints. |
|
|
1556 C |
|
|
1557 IF (LENIC .GT. 0) THEN |
|
|
1558 CALL DCNST0(NEQ,Y,IWORK(LICNS),IRET) |
|
|
1559 IF (IRET .NE. 0) GO TO 727 |
|
|
1560 ENDIF |
|
|
1561 C |
|
|
1562 C Check ID for legality. |
|
|
1563 C |
|
|
1564 IF (LENID .GT. 0) THEN |
|
|
1565 DO 50 I = 1,NEQ |
|
|
1566 IDI = IWORK(LID-1+I) |
|
|
1567 IF (IDI .NE. 1 .AND. IDI .NE. -1) GO TO 724 |
|
|
1568 50 CONTINUE |
|
|
1569 ENDIF |
|
|
1570 C |
|
|
1571 C Check to see that TOUT is different from T. |
|
|
1572 C |
|
|
1573 IF(TOUT .EQ. T)GO TO 719 |
|
|
1574 C |
|
|
1575 C Check HMAX. |
|
|
1576 C |
|
|
1577 IF(INFO(7) .NE. 0) THEN |
|
|
1578 HMAX = RWORK(LHMAX) |
|
|
1579 IF (HMAX .LE. 0.0D0) GO TO 710 |
|
|
1580 ENDIF |
|
|
1581 C |
|
|
1582 C Initialize counters and other flags. |
|
|
1583 C |
|
|
1584 IWORK(LNST)=0 |
|
|
1585 IWORK(LNRE)=0 |
|
|
1586 IWORK(LNJE)=0 |
|
|
1587 IWORK(LETF)=0 |
|
|
1588 IWORK(LNCFN)=0 |
|
|
1589 IWORK(LNNI)=0 |
|
|
1590 IWORK(LNLI)=0 |
|
|
1591 IWORK(LNPS)=0 |
|
|
1592 IWORK(LNCFL)=0 |
|
|
1593 IWORK(LKPRIN)=INFO(18) |
|
|
1594 IDID=1 |
|
|
1595 GO TO 200 |
|
|
1596 C |
|
|
1597 C----------------------------------------------------------------------- |
|
|
1598 C This block is for continuation calls only. |
|
|
1599 C Here we check INFO(1), and if the last step was interrupted, |
|
|
1600 C we check whether appropriate action was taken. |
|
|
1601 C----------------------------------------------------------------------- |
|
|
1602 C |
|
|
1603 100 CONTINUE |
|
|
1604 IF(INFO(1).EQ.1)GO TO 110 |
|
|
1605 ITEMP = 1 |
|
|
1606 IF(INFO(1).NE.-1)GO TO 701 |
|
|
1607 C |
|
|
1608 C If we are here, the last step was interrupted by an error |
|
|
1609 C condition from DDSTP, and appropriate action was not taken. |
|
|
1610 C This is a fatal error. |
|
|
1611 C |
|
|
1612 MSG = 'DASPK-- THE LAST STEP TERMINATED WITH A NEGATIVE' |
|
|
1613 CALL XERRWD(MSG,49,201,0,0,0,0,0,0.0D0,0.0D0) |
|
|
1614 MSG = 'DASPK-- VALUE (=I1) OF IDID AND NO APPROPRIATE' |
|
|
1615 CALL XERRWD(MSG,47,202,0,1,IDID,0,0,0.0D0,0.0D0) |
|
|
1616 MSG = 'DASPK-- ACTION WAS TAKEN. RUN TERMINATED' |
|
|
1617 CALL XERRWD(MSG,41,203,1,0,0,0,0,0.0D0,0.0D0) |
|
|
1618 RETURN |
|
|
1619 110 CONTINUE |
|
|
1620 C |
|
|
1621 C----------------------------------------------------------------------- |
|
|
1622 C This block is executed on all calls. |
|
|
1623 C |
|
|
1624 C Counters are saved for later checks of performance. |
|
|
1625 C Then the error tolerance parameters are checked, and the |
|
|
1626 C work array pointers are set. |
|
|
1627 C----------------------------------------------------------------------- |
|
|
1628 C |
|
|
1629 200 CONTINUE |
|
|
1630 C |
|
|
1631 C Save counters for use later. |
|
|
1632 C |
|
|
1633 IWORK(LNSTL)=IWORK(LNST) |
|
|
1634 NLI0 = IWORK(LNLI) |
|
|
1635 NNI0 = IWORK(LNNI) |
|
|
1636 NCFN0 = IWORK(LNCFN) |
|
|
1637 NCFL0 = IWORK(LNCFL) |
|
|
1638 NWARN = 0 |
|
|
1639 C |
|
|
1640 C Check RTOL and ATOL. |
|
|
1641 C |
|
|
1642 NZFLG = 0 |
|
|
1643 RTOLI = RTOL(1) |
|
|
1644 ATOLI = ATOL(1) |
|
|
1645 DO 210 I=1,NEQ |
|
|
1646 IF (INFO(2) .EQ. 1) RTOLI = RTOL(I) |
|
|
1647 IF (INFO(2) .EQ. 1) ATOLI = ATOL(I) |
|
|
1648 IF (RTOLI .GT. 0.0D0 .OR. ATOLI .GT. 0.0D0) NZFLG = 1 |
|
|
1649 IF (RTOLI .LT. 0.0D0) GO TO 706 |
|
|
1650 IF (ATOLI .LT. 0.0D0) GO TO 707 |
|
|
1651 210 CONTINUE |
|
|
1652 IF (NZFLG .EQ. 0) GO TO 708 |
|
|
1653 C |
|
|
1654 C Set pointers to RWORK and IWORK segments. |
|
|
1655 C For direct methods, SAVR is not used. |
|
|
1656 C |
|
|
1657 IWORK(LLCIWP) = LID + LENID |
|
|
1658 LSAVR = LDELTA |
|
|
1659 IF (INFO(12) .NE. 0) LSAVR = LDELTA + NEQ |
|
|
1660 LE = LSAVR + NEQ |
|
|
1661 LWT = LE + NEQ |
|
|
1662 LVT = LWT |
|
|
1663 IF (INFO(16) .NE. 0) LVT = LWT + NEQ |
|
|
1664 LPHI = LVT + NEQ |
|
|
1665 LWM = LPHI + (IWORK(LMXORD)+1)*NEQ |
|
|
1666 IF (INFO(1) .EQ. 1) GO TO 400 |
|
|
1667 C |
|
|
1668 C----------------------------------------------------------------------- |
|
|
1669 C This block is executed on the initial call only. |
|
|
1670 C Set the initial step size, the error weight vector, and PHI. |
|
|
1671 C Compute unknown initial components of Y and YPRIME, if requested. |
|
|
1672 C----------------------------------------------------------------------- |
|
|
1673 C |
|
|
1674 300 CONTINUE |
|
|
1675 TN=T |
|
|
1676 IDID=1 |
|
|
1677 C |
|
|
1678 C Set error weight array WT and altered weight array VT. |
|
|
1679 C |
|
|
1680 CALL DDAWTS(NEQ,INFO(2),RTOL,ATOL,Y,RWORK(LWT),RPAR,IPAR) |
|
|
1681 CALL DINVWT(NEQ,RWORK(LWT),IER) |
|
|
1682 IF (IER .NE. 0) GO TO 713 |
|
|
1683 IF (INFO(16) .NE. 0) THEN |
|
|
1684 DO 305 I = 1, NEQ |
|
|
1685 305 RWORK(LVT+I-1) = MAX(IWORK(LID+I-1),0)*RWORK(LWT+I-1) |
|
|
1686 ENDIF |
|
|
1687 C |
|
|
1688 C Compute unit roundoff and HMIN. |
|
|
1689 C |
|
|
1690 UROUND = D1MACH(4) |
|
|
1691 RWORK(LROUND) = UROUND |
|
|
1692 HMIN = 4.0D0*UROUND*MAX(ABS(T),ABS(TOUT)) |
|
|
1693 C |
|
|
1694 C Set/check STPTOL control for initial condition calculation. |
|
|
1695 C |
|
|
1696 IF (INFO(11) .NE. 0) THEN |
|
|
1697 IF( INFO(17) .EQ. 0) THEN |
|
|
1698 RWORK(LSTOL) = UROUND**.6667D0 |
|
|
1699 ELSE |
|
|
1700 IF (RWORK(LSTOL) .LE. 0.0D0) GO TO 725 |
|
|
1701 ENDIF |
|
|
1702 ENDIF |
|
|
1703 C |
|
|
1704 C Compute EPCON and square root of NEQ and its reciprocal, used |
|
|
1705 C inside iterative solver. |
|
|
1706 C |
|
|
1707 RWORK(LEPCON) = 0.33D0 |
|
|
1708 FLOATN = NEQ |
|
|
1709 RWORK(LSQRN) = SQRT(FLOATN) |
|
|
1710 RWORK(LRSQRN) = 1.D0/RWORK(LSQRN) |
|
|
1711 C |
|
|
1712 C Check initial interval to see that it is long enough. |
|
|
1713 C |
|
|
1714 TDIST = ABS(TOUT - T) |
|
|
1715 IF(TDIST .LT. HMIN) GO TO 714 |
|
|
1716 C |
|
|
1717 C Check H0, if this was input. |
|
|
1718 C |
|
|
1719 IF (INFO(8) .EQ. 0) GO TO 310 |
|
|
1720 H0 = RWORK(LH) |
|
|
1721 IF ((TOUT - T)*H0 .LT. 0.0D0) GO TO 711 |
|
|
1722 IF (H0 .EQ. 0.0D0) GO TO 712 |
|
|
1723 GO TO 320 |
|
|
1724 310 CONTINUE |
|
|
1725 C |
|
|
1726 C Compute initial stepsize, to be used by either |
|
|
1727 C DDSTP or DDASIC, depending on INFO(11). |
|
|
1728 C |
|
|
1729 H0 = 0.001D0*TDIST |
|
|
1730 YPNORM = DDWNRM(NEQ,YPRIME,RWORK(LVT),RPAR,IPAR) |
|
|
1731 IF (YPNORM .GT. 0.5D0/H0) H0 = 0.5D0/YPNORM |
|
|
1732 H0 = SIGN(H0,TOUT-T) |
|
|
1733 C |
|
|
1734 C Adjust H0 if necessary to meet HMAX bound. |
|
|
1735 C |
|
|
1736 320 IF (INFO(7) .EQ. 0) GO TO 330 |
|
|
1737 RH = ABS(H0)/RWORK(LHMAX) |
|
|
1738 IF (RH .GT. 1.0D0) H0 = H0/RH |
|
|
1739 C |
|
|
1740 C Check against TSTOP, if applicable. |
|
|
1741 C |
|
|
1742 330 IF (INFO(4) .EQ. 0) GO TO 340 |
|
|
1743 TSTOP = RWORK(LTSTOP) |
|
|
1744 IF ((TSTOP - T)*H0 .LT. 0.0D0) GO TO 715 |
|
|
1745 IF ((T + H0 - TSTOP)*H0 .GT. 0.0D0) H0 = TSTOP - T |
|
|
1746 IF ((TSTOP - TOUT)*H0 .LT. 0.0D0) GO TO 709 |
|
|
1747 C |
|
|
1748 340 IF (INFO(11) .EQ. 0) GO TO 370 |
|
|
1749 C |
|
|
1750 C Compute unknown components of initial Y and YPRIME, depending |
|
|
1751 C on INFO(11) and INFO(12). INFO(12) represents the nonlinear |
|
|
1752 C solver type (direct/Krylov). Pass the name of the specific |
|
|
1753 C nonlinear solver, depending on INFO(12). The location of the work |
|
|
1754 C arrays SAVR, YIC, YPIC, PWK also differ in the two cases. |
|
|
1755 C |
|
|
1756 NWT = 1 |
|
|
1757 EPCONI = RWORK(LEPIN)*RWORK(LEPCON) |
|
|
1758 350 IF (INFO(12) .EQ. 0) THEN |
|
|
1759 LYIC = LPHI + 2*NEQ |
|
|
1760 LYPIC = LYIC + NEQ |
|
|
1761 LPWK = LYPIC |
|
|
1762 CALL DDASIC(TN,Y,YPRIME,NEQ,INFO(11),IWORK(LID), |
|
|
1763 * RES,JAC,PSOL,H0,RWORK(LWT),NWT,IDID,RPAR,IPAR, |
|
|
1764 * RWORK(LPHI),RWORK(LSAVR),RWORK(LDELTA),RWORK(LE), |
|
|
1765 * RWORK(LYIC),RWORK(LYPIC),RWORK(LPWK),RWORK(LWM),IWORK(LIWM), |
|
|
1766 * HMIN,RWORK(LROUND),RWORK(LEPLI),RWORK(LSQRN),RWORK(LRSQRN), |
|
|
1767 * EPCONI,RWORK(LSTOL),INFO(15),ICNFLG,IWORK(LICNS),DDASID) |
|
|
1768 ELSE IF (INFO(12) .EQ. 1) THEN |
|
|
1769 LYIC = LWM |
|
|
1770 LYPIC = LYIC + NEQ |
|
|
1771 LPWK = LYPIC + NEQ |
|
|
1772 CALL DDASIC(TN,Y,YPRIME,NEQ,INFO(11),IWORK(LID), |
|
|
1773 * RES,JAC,PSOL,H0,RWORK(LWT),NWT,IDID,RPAR,IPAR, |
|
|
1774 * RWORK(LPHI),RWORK(LSAVR),RWORK(LDELTA),RWORK(LE), |
|
|
1775 * RWORK(LYIC),RWORK(LYPIC),RWORK(LPWK),RWORK(LWM),IWORK(LIWM), |
|
|
1776 * HMIN,RWORK(LROUND),RWORK(LEPLI),RWORK(LSQRN),RWORK(LRSQRN), |
|
|
1777 * EPCONI,RWORK(LSTOL),INFO(15),ICNFLG,IWORK(LICNS),DDASIK) |
|
|
1778 ENDIF |
|
|
1779 C |
|
|
1780 IF (IDID .LT. 0) GO TO 600 |
|
|
1781 C |
|
|
1782 C DDASIC was successful. If this was the first call to DDASIC, |
|
|
1783 C update the WT array (with the current Y) and call it again. |
|
|
1784 C |
|
|
1785 IF (NWT .EQ. 2) GO TO 355 |
|
|
1786 NWT = 2 |
|
|
1787 CALL DDAWTS(NEQ,INFO(2),RTOL,ATOL,Y,RWORK(LWT),RPAR,IPAR) |
|
|
1788 CALL DINVWT(NEQ,RWORK(LWT),IER) |
|
|
1789 IF (IER .NE. 0) GO TO 713 |
|
|
1790 GO TO 350 |
|
|
1791 C |
|
|
1792 C If INFO(14) = 1, return now with IDID = 4. |
|
|
1793 C |
|
|
1794 355 IF (INFO(14) .EQ. 1) THEN |
|
|
1795 IDID = 4 |
|
|
1796 H = H0 |
|
|
1797 IF (INFO(11) .EQ. 1) RWORK(LHOLD) = H0 |
|
|
1798 GO TO 590 |
|
|
1799 ENDIF |
|
|
1800 C |
|
|
1801 C Update the WT and VT arrays one more time, with the new Y. |
|
|
1802 C |
|
|
1803 CALL DDAWTS(NEQ,INFO(2),RTOL,ATOL,Y,RWORK(LWT),RPAR,IPAR) |
|
|
1804 CALL DINVWT(NEQ,RWORK(LWT),IER) |
|
|
1805 IF (IER .NE. 0) GO TO 713 |
|
|
1806 IF (INFO(16) .NE. 0) THEN |
|
|
1807 DO 357 I = 1, NEQ |
|
|
1808 357 RWORK(LVT+I-1) = MAX(IWORK(LID+I-1),0)*RWORK(LWT+I-1) |
|
|
1809 ENDIF |
|
|
1810 C |
|
|
1811 C Reset the initial stepsize to be used by DDSTP. |
|
|
1812 C Use H0, if this was input. Otherwise, recompute H0, |
|
|
1813 C and adjust it if necessary to meet HMAX bound. |
|
|
1814 C |
|
|
1815 IF (INFO(8) .NE. 0) THEN |
|
|
1816 H0 = RWORK(LH) |
|
|
1817 GO TO 360 |
|
|
1818 ENDIF |
|
|
1819 C |
|
|
1820 H0 = 0.001D0*TDIST |
|
|
1821 YPNORM = DDWNRM(NEQ,YPRIME,RWORK(LVT),RPAR,IPAR) |
|
|
1822 IF (YPNORM .GT. 0.5D0/H0) H0 = 0.5D0/YPNORM |
|
|
1823 H0 = SIGN(H0,TOUT-T) |
|
|
1824 C |
|
|
1825 360 IF (INFO(7) .NE. 0) THEN |
|
|
1826 RH = ABS(H0)/RWORK(LHMAX) |
|
|
1827 IF (RH .GT. 1.0D0) H0 = H0/RH |
|
|
1828 ENDIF |
|
|
1829 C |
|
|
1830 C Check against TSTOP, if applicable. |
|
|
1831 C |
|
|
1832 IF (INFO(4) .NE. 0) THEN |
|
|
1833 TSTOP = RWORK(LTSTOP) |
|
|
1834 IF ((T + H0 - TSTOP)*H0 .GT. 0.0D0) H0 = TSTOP - T |
|
|
1835 ENDIF |
|
|
1836 C |
|
|
1837 C Load H and RWORK(LH) with H0. |
|
|
1838 C |
|
|
1839 370 H = H0 |
|
|
1840 RWORK(LH) = H |
|
|
1841 C |
|
|
1842 C Load Y and H*YPRIME into PHI(*,1) and PHI(*,2). |
|
|
1843 C |
|
|
1844 ITEMP = LPHI + NEQ |
|
|
1845 DO 380 I = 1,NEQ |
|
|
1846 RWORK(LPHI + I - 1) = Y(I) |
|
|
1847 380 RWORK(ITEMP + I - 1) = H*YPRIME(I) |
|
|
1848 C |
|
|
1849 GO TO 500 |
|
|
1850 C |
|
|
1851 C----------------------------------------------------------------------- |
|
|
1852 C This block is for continuation calls only. |
|
|
1853 C Its purpose is to check stop conditions before taking a step. |
|
|
1854 C Adjust H if necessary to meet HMAX bound. |
|
|
1855 C----------------------------------------------------------------------- |
|
|
1856 C |
|
|
1857 400 CONTINUE |
|
|
1858 UROUND=RWORK(LROUND) |
|
|
1859 DONE = .FALSE. |
|
|
1860 TN=RWORK(LTN) |
|
|
1861 H=RWORK(LH) |
|
|
1862 IF(INFO(7) .EQ. 0) GO TO 410 |
|
|
1863 RH = ABS(H)/RWORK(LHMAX) |
|
|
1864 IF(RH .GT. 1.0D0) H = H/RH |
|
|
1865 410 CONTINUE |
|
|
1866 IF(T .EQ. TOUT) GO TO 719 |
|
|
1867 IF((T - TOUT)*H .GT. 0.0D0) GO TO 711 |
|
|
1868 IF(INFO(4) .EQ. 1) GO TO 430 |
|
|
1869 IF(INFO(3) .EQ. 1) GO TO 420 |
|
|
1870 IF((TN-TOUT)*H.LT.0.0D0)GO TO 490 |
|
|
1871 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD), |
|
|
1872 * RWORK(LPHI),RWORK(LPSI)) |
|
|
1873 T=TOUT |
|
|
1874 IDID = 3 |
|
|
1875 DONE = .TRUE. |
|
|
1876 GO TO 490 |
|
|
1877 420 IF((TN-T)*H .LE. 0.0D0) GO TO 490 |
|
|
1878 IF((TN - TOUT)*H .GT. 0.0D0) GO TO 425 |
|
|
1879 CALL DDATRP(TN,TN,Y,YPRIME,NEQ,IWORK(LKOLD), |
|
|
1880 * RWORK(LPHI),RWORK(LPSI)) |
|
|
1881 T = TN |
|
|
1882 IDID = 1 |
|
|
1883 DONE = .TRUE. |
|
|
1884 GO TO 490 |
|
|
1885 425 CONTINUE |
|
|
1886 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD), |
|
|
1887 * RWORK(LPHI),RWORK(LPSI)) |
|
|
1888 T = TOUT |
|
|
1889 IDID = 3 |
|
|
1890 DONE = .TRUE. |
|
|
1891 GO TO 490 |
|
|
1892 430 IF(INFO(3) .EQ. 1) GO TO 440 |
|
|
1893 TSTOP=RWORK(LTSTOP) |
|
|
1894 IF((TN-TSTOP)*H.GT.0.0D0) GO TO 715 |
|
|
1895 IF((TSTOP-TOUT)*H.LT.0.0D0)GO TO 709 |
|
|
1896 IF((TN-TOUT)*H.LT.0.0D0)GO TO 450 |
|
|
1897 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD), |
|
|
1898 * RWORK(LPHI),RWORK(LPSI)) |
|
|
1899 T=TOUT |
|
|
1900 IDID = 3 |
|
|
1901 DONE = .TRUE. |
|
|
1902 GO TO 490 |
|
|
1903 440 TSTOP = RWORK(LTSTOP) |
|
|
1904 IF((TN-TSTOP)*H .GT. 0.0D0) GO TO 715 |
|
|
1905 IF((TSTOP-TOUT)*H .LT. 0.0D0) GO TO 709 |
|
|
1906 IF((TN-T)*H .LE. 0.0D0) GO TO 450 |
|
|
1907 IF((TN - TOUT)*H .GT. 0.0D0) GO TO 445 |
|
|
1908 CALL DDATRP(TN,TN,Y,YPRIME,NEQ,IWORK(LKOLD), |
|
|
1909 * RWORK(LPHI),RWORK(LPSI)) |
|
|
1910 T = TN |
|
|
1911 IDID = 1 |
|
|
1912 DONE = .TRUE. |
|
|
1913 GO TO 490 |
|
|
1914 445 CONTINUE |
|
|
1915 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD), |
|
|
1916 * RWORK(LPHI),RWORK(LPSI)) |
|
|
1917 T = TOUT |
|
|
1918 IDID = 3 |
|
|
1919 DONE = .TRUE. |
|
|
1920 GO TO 490 |
|
|
1921 450 CONTINUE |
|
|
1922 C |
|
|
1923 C Check whether we are within roundoff of TSTOP. |
|
|
1924 C |
|
|
1925 IF(ABS(TN-TSTOP).GT.100.0D0*UROUND* |
|
|
1926 * (ABS(TN)+ABS(H)))GO TO 460 |
|
|
1927 CALL DDATRP(TN,TSTOP,Y,YPRIME,NEQ,IWORK(LKOLD), |
|
|
1928 * RWORK(LPHI),RWORK(LPSI)) |
|
|
1929 IDID=2 |
|
|
1930 T=TSTOP |
|
|
1931 DONE = .TRUE. |
|
|
1932 GO TO 490 |
|
|
1933 460 TNEXT=TN+H |
|
|
1934 IF((TNEXT-TSTOP)*H.LE.0.0D0)GO TO 490 |
|
|
1935 H=TSTOP-TN |
|
|
1936 RWORK(LH)=H |
|
|
1937 C |
|
|
1938 490 IF (DONE) GO TO 590 |
|
|
1939 C |
|
|
1940 C----------------------------------------------------------------------- |
|
|
1941 C The next block contains the call to the one-step integrator DDSTP. |
|
|
1942 C This is a looping point for the integration steps. |
|
|
1943 C Check for too many steps. |
|
|
1944 C Check for poor Newton/Krylov performance. |
|
|
1945 C Update WT. Check for too much accuracy requested. |
|
|
1946 C Compute minimum stepsize. |
|
|
1947 C----------------------------------------------------------------------- |
|
|
1948 C |
|
|
1949 500 CONTINUE |
|
|
1950 C |
|
|
1951 C Check for too many steps. |
|
|
1952 C |
|
|
1953 IF((IWORK(LNST)-IWORK(LNSTL)).LT.500) GO TO 505 |
|
|
1954 IDID=-1 |
|
|
1955 GO TO 527 |
|
|
1956 C |
|
|
1957 C Check for poor Newton/Krylov performance. |
|
|
1958 C |
|
|
1959 505 IF (INFO(12) .EQ. 0) GO TO 510 |
|
|
1960 NSTD = IWORK(LNST) - IWORK(LNSTL) |
|
|
1961 NNID = IWORK(LNNI) - NNI0 |
|
|
1962 IF (NSTD .LT. 10 .OR. NNID .EQ. 0) GO TO 510 |
|
|
1963 AVLIN = REAL(IWORK(LNLI) - NLI0)/REAL(NNID) |
|
|
1964 RCFN = REAL(IWORK(LNCFN) - NCFN0)/REAL(NSTD) |
|
|
1965 RCFL = REAL(IWORK(LNCFL) - NCFL0)/REAL(NNID) |
|
|
1966 FMAXL = IWORK(LMAXL) |
|
|
1967 LAVL = AVLIN .GT. FMAXL |
|
|
1968 LCFN = RCFN .GT. 0.9D0 |
|
|
1969 LCFL = RCFL .GT. 0.9D0 |
|
|
1970 LWARN = LAVL .OR. LCFN .OR. LCFL |
|
|
1971 IF (.NOT.LWARN) GO TO 510 |
|
|
1972 NWARN = NWARN + 1 |
|
|
1973 IF (NWARN .GT. 10) GO TO 510 |
|
|
1974 IF (LAVL) THEN |
|
|
1975 MSG = 'DASPK-- Warning. Poor iterative algorithm performance ' |
|
|
1976 CALL XERRWD (MSG, 56, 501, 0, 0, 0, 0, 0, 0.0D0, 0.0D0) |
|
|
1977 MSG = ' at T = R1. Average no. of linear iterations = R2 ' |
|
|
1978 CALL XERRWD (MSG, 56, 501, 0, 0, 0, 0, 2, TN, AVLIN) |
|
|
1979 ENDIF |
|
|
1980 IF (LCFN) THEN |
|
|
1981 MSG = 'DASPK-- Warning. Poor iterative algorithm performance ' |
|
|
1982 CALL XERRWD (MSG, 56, 502, 0, 0, 0, 0, 0, 0.0D0, 0.0D0) |
|
|
1983 MSG = ' at T = R1. Nonlinear convergence failure rate = R2' |
|
|
1984 CALL XERRWD (MSG, 56, 502, 0, 0, 0, 0, 2, TN, RCFN) |
|
|
1985 ENDIF |
|
|
1986 IF (LCFL) THEN |
|
|
1987 MSG = 'DASPK-- Warning. Poor iterative algorithm performance ' |
|
|
1988 CALL XERRWD (MSG, 56, 503, 0, 0, 0, 0, 0, 0.0D0, 0.0D0) |
|
|
1989 MSG = ' at T = R1. Linear convergence failure rate = R2 ' |
|
|
1990 CALL XERRWD (MSG, 56, 503, 0, 0, 0, 0, 2, TN, RCFL) |
|
|
1991 ENDIF |
|
|
1992 C |
|
|
1993 C Update WT and VT, if this is not the first call. |
|
|
1994 C |
|
|
1995 510 CALL DDAWTS(NEQ,INFO(2),RTOL,ATOL,RWORK(LPHI),RWORK(LWT), |
|
|
1996 * RPAR,IPAR) |
|
|
1997 CALL DINVWT(NEQ,RWORK(LWT),IER) |
|
|
1998 IF (IER .NE. 0) THEN |
|
|
1999 IDID = -3 |
|
|
2000 GO TO 527 |
|
|
2001 ENDIF |
|
|
2002 IF (INFO(16) .NE. 0) THEN |
|
|
2003 DO 515 I = 1, NEQ |
|
|
2004 515 RWORK(LVT+I-1) = MAX(IWORK(LID+I-1),0)*RWORK(LWT+I-1) |
|
|
2005 ENDIF |
|
|
2006 C |
|
|
2007 C Test for too much accuracy requested. |
|
|
2008 C |
|
|
2009 R = DDWNRM(NEQ,RWORK(LPHI),RWORK(LWT),RPAR,IPAR)*100.0D0*UROUND |
|
|
2010 IF (R .LE. 1.0D0) GO TO 525 |
|
|
2011 C |
|
|
2012 C Multiply RTOL and ATOL by R and return. |
|
|
2013 C |
|
|
2014 IF(INFO(2).EQ.1)GO TO 523 |
|
|
2015 RTOL(1)=R*RTOL(1) |
|
|
2016 ATOL(1)=R*ATOL(1) |
|
|
2017 IDID=-2 |
|
|
2018 GO TO 527 |
|
|
2019 523 DO 524 I=1,NEQ |
|
|
2020 RTOL(I)=R*RTOL(I) |
|
|
2021 524 ATOL(I)=R*ATOL(I) |
|
|
2022 IDID=-2 |
|
|
2023 GO TO 527 |
|
|
2024 525 CONTINUE |
|
|
2025 C |
|
|
2026 C Compute minimum stepsize. |
|
|
2027 C |
|
|
2028 HMIN=4.0D0*UROUND*MAX(ABS(TN),ABS(TOUT)) |
|
|
2029 C |
|
|
2030 C Test H vs. HMAX |
|
|
2031 IF (INFO(7) .NE. 0) THEN |
|
|
2032 RH = ABS(H)/RWORK(LHMAX) |
|
|
2033 IF (RH .GT. 1.0D0) H = H/RH |
|
|
2034 ENDIF |
|
|
2035 C |
|
|
2036 C Call the one-step integrator. |
|
|
2037 C Note that INFO(12) represents the nonlinear solver type. |
|
|
2038 C Pass the required nonlinear solver, depending upon INFO(12). |
|
|
2039 C |
|
|
2040 IF (INFO(12) .EQ. 0) THEN |
|
|
2041 CALL DDSTP(TN,Y,YPRIME,NEQ, |
|
|
2042 * RES,JAC,PSOL,H,RWORK(LWT),RWORK(LVT),INFO(1),IDID,RPAR,IPAR, |
|
|
2043 * RWORK(LPHI),RWORK(LSAVR),RWORK(LDELTA),RWORK(LE), |
|
|
2044 * RWORK(LWM),IWORK(LIWM), |
|
|
2045 * RWORK(LALPHA),RWORK(LBETA),RWORK(LGAMMA), |
|
|
2046 * RWORK(LPSI),RWORK(LSIGMA), |
|
|
2047 * RWORK(LCJ),RWORK(LCJOLD),RWORK(LHOLD),RWORK(LS),HMIN, |
|
|
2048 * RWORK(LROUND), RWORK(LEPLI),RWORK(LSQRN),RWORK(LRSQRN), |
|
|
2049 * RWORK(LEPCON), IWORK(LPHASE),IWORK(LJCALC),INFO(15), |
|
|
2050 * IWORK(LK), IWORK(LKOLD),IWORK(LNS),NONNEG,INFO(12), |
|
|
2051 * DNEDD) |
|
|
2052 ELSE IF (INFO(12) .EQ. 1) THEN |
|
|
2053 CALL DDSTP(TN,Y,YPRIME,NEQ, |
|
|
2054 * RES,JAC,PSOL,H,RWORK(LWT),RWORK(LVT),INFO(1),IDID,RPAR,IPAR, |
|
|
2055 * RWORK(LPHI),RWORK(LSAVR),RWORK(LDELTA),RWORK(LE), |
|
|
2056 * RWORK(LWM),IWORK(LIWM), |
|
|
2057 * RWORK(LALPHA),RWORK(LBETA),RWORK(LGAMMA), |
|
|
2058 * RWORK(LPSI),RWORK(LSIGMA), |
|
|
2059 * RWORK(LCJ),RWORK(LCJOLD),RWORK(LHOLD),RWORK(LS),HMIN, |
|
|
2060 * RWORK(LROUND), RWORK(LEPLI),RWORK(LSQRN),RWORK(LRSQRN), |
|
|
2061 * RWORK(LEPCON), IWORK(LPHASE),IWORK(LJCALC),INFO(15), |
|
|
2062 * IWORK(LK), IWORK(LKOLD),IWORK(LNS),NONNEG,INFO(12), |
|
|
2063 * DNEDK) |
|
|
2064 ENDIF |
|
|
2065 C |
|
|
2066 527 IF(IDID.LT.0)GO TO 600 |
|
|
2067 C |
|
|
2068 C----------------------------------------------------------------------- |
|
|
2069 C This block handles the case of a successful return from DDSTP |
|
|
2070 C (IDID=1). Test for stop conditions. |
|
|
2071 C----------------------------------------------------------------------- |
|
|
2072 C |
|
|
2073 IF(INFO(4).NE.0)GO TO 540 |
|
|
2074 IF(INFO(3).NE.0)GO TO 530 |
|
|
2075 IF((TN-TOUT)*H.LT.0.0D0)GO TO 500 |
|
|
2076 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ, |
|
|
2077 * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI)) |
|
|
2078 IDID=3 |
|
|
2079 T=TOUT |
|
|
2080 GO TO 580 |
|
|
2081 530 IF((TN-TOUT)*H.GE.0.0D0)GO TO 535 |
|
|
2082 T=TN |
|
|
2083 IDID=1 |
|
|
2084 GO TO 580 |
|
|
2085 535 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ, |
|
|
2086 * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI)) |
|
|
2087 IDID=3 |
|
|
2088 T=TOUT |
|
|
2089 GO TO 580 |
|
|
2090 540 IF(INFO(3).NE.0)GO TO 550 |
|
|
2091 IF((TN-TOUT)*H.LT.0.0D0)GO TO 542 |
|
|
2092 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ, |
|
|
2093 * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI)) |
|
|
2094 T=TOUT |
|
|
2095 IDID=3 |
|
|
2096 GO TO 580 |
|
|
2097 542 IF(ABS(TN-TSTOP).LE.100.0D0*UROUND* |
|
|
2098 * (ABS(TN)+ABS(H)))GO TO 545 |
|
|
2099 TNEXT=TN+H |
|
|
2100 IF((TNEXT-TSTOP)*H.LE.0.0D0)GO TO 500 |
|
|
2101 H=TSTOP-TN |
|
|
2102 GO TO 500 |
|
|
2103 545 CALL DDATRP(TN,TSTOP,Y,YPRIME,NEQ, |
|
|
2104 * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI)) |
|
|
2105 IDID=2 |
|
|
2106 T=TSTOP |
|
|
2107 GO TO 580 |
|
|
2108 550 IF((TN-TOUT)*H.GE.0.0D0)GO TO 555 |
|
|
2109 IF(ABS(TN-TSTOP).LE.100.0D0*UROUND*(ABS(TN)+ABS(H)))GO TO 552 |
|
|
2110 T=TN |
|
|
2111 IDID=1 |
|
|
2112 GO TO 580 |
|
|
2113 552 CALL DDATRP(TN,TSTOP,Y,YPRIME,NEQ, |
|
|
2114 * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI)) |
|
|
2115 IDID=2 |
|
|
2116 T=TSTOP |
|
|
2117 GO TO 580 |
|
|
2118 555 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ, |
|
|
2119 * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI)) |
|
|
2120 T=TOUT |
|
|
2121 IDID=3 |
|
|
2122 580 CONTINUE |
|
|
2123 C |
|
|
2124 C----------------------------------------------------------------------- |
|
|
2125 C All successful returns from DDASPK are made from this block. |
|
|
2126 C----------------------------------------------------------------------- |
|
|
2127 C |
|
|
2128 590 CONTINUE |
|
|
2129 RWORK(LTN)=TN |
|
|
2130 RWORK(LH)=H |
|
|
2131 RETURN |
|
|
2132 C |
|
|
2133 C----------------------------------------------------------------------- |
|
|
2134 C This block handles all unsuccessful returns other than for |
|
|
2135 C illegal input. |
|
|
2136 C----------------------------------------------------------------------- |
|
|
2137 C |
|
|
2138 600 CONTINUE |
|
|
2139 ITEMP = -IDID |
|
|
2140 GO TO (610,620,630,700,655,640,650,660,670,675, |
|
|
2141 * 680,685,690,695), ITEMP |
|
|
2142 C |
|
|
2143 C The maximum number of steps was taken before |
|
|
2144 C reaching tout. |
|
|
2145 C |
|
|
2146 610 MSG = 'DASPK-- AT CURRENT T (=R1) 500 STEPS' |
|
|
2147 CALL XERRWD(MSG,38,610,0,0,0,0,1,TN,0.0D0) |
|
|
2148 MSG = 'DASPK-- TAKEN ON THIS CALL BEFORE REACHING TOUT' |
|
|
2149 CALL XERRWD(MSG,48,611,0,0,0,0,0,0.0D0,0.0D0) |
|
|
2150 GO TO 700 |
|
|
2151 C |
|
|
2152 C Too much accuracy for machine precision. |
|
|
2153 C |
|
|
2154 620 MSG = 'DASPK-- AT T (=R1) TOO MUCH ACCURACY REQUESTED' |
|
|
2155 CALL XERRWD(MSG,47,620,0,0,0,0,1,TN,0.0D0) |
|
|
2156 MSG = 'DASPK-- FOR PRECISION OF MACHINE. RTOL AND ATOL' |
|
|
2157 CALL XERRWD(MSG,48,621,0,0,0,0,0,0.0D0,0.0D0) |
|
|
2158 MSG = 'DASPK-- WERE INCREASED TO APPROPRIATE VALUES' |
|
|
2159 CALL XERRWD(MSG,45,622,0,0,0,0,0,0.0D0,0.0D0) |
|
|
2160 GO TO 700 |
|
|
2161 C |
|
|
2162 C WT(I) .LE. 0.0D0 for some I (not at start of problem). |
|
|
2163 C |
|
|
2164 630 MSG = 'DASPK-- AT T (=R1) SOME ELEMENT OF WT' |
|
|
2165 CALL XERRWD(MSG,38,630,0,0,0,0,1,TN,0.0D0) |
|
|
2166 MSG = 'DASPK-- HAS BECOME .LE. 0.0' |
|
|
2167 CALL XERRWD(MSG,28,631,0,0,0,0,0,0.0D0,0.0D0) |
|
|
2168 GO TO 700 |
|
|
2169 C |
|
|
2170 C Error test failed repeatedly or with H=HMIN. |
|
|
2171 C |
|
|
2172 640 MSG = 'DASPK-- AT T (=R1) AND STEPSIZE H (=R2) THE' |
|
|
2173 CALL XERRWD(MSG,44,640,0,0,0,0,2,TN,H) |
|
|
2174 MSG='DASPK-- ERROR TEST FAILED REPEATEDLY OR WITH ABS(H)=HMIN' |
|
|
2175 CALL XERRWD(MSG,57,641,0,0,0,0,0,0.0D0,0.0D0) |
|
|
2176 GO TO 700 |
|
|
2177 C |
|
|
2178 C Nonlinear solver failed to converge repeatedly or with H=HMIN. |
|
|
2179 C |
|
|
2180 650 MSG = 'DASPK-- AT T (=R1) AND STEPSIZE H (=R2) THE' |
|
|
2181 CALL XERRWD(MSG,44,650,0,0,0,0,2,TN,H) |
|
|
2182 MSG = 'DASPK-- NONLINEAR SOLVER FAILED TO CONVERGE' |
|
|
2183 CALL XERRWD(MSG,44,651,0,0,0,0,0,0.0D0,0.0D0) |
|
|
2184 MSG = 'DASPK-- REPEATEDLY OR WITH ABS(H)=HMIN' |
|
|
2185 CALL XERRWD(MSG,40,652,0,0,0,0,0,0.0D0,0.0D0) |
|
|
2186 GO TO 700 |
|
|
2187 C |
|
|
2188 C The preconditioner had repeated failures. |
|
|
2189 C |
|
|
2190 655 MSG = 'DASPK-- AT T (=R1) AND STEPSIZE H (=R2) THE' |
|
|
2191 CALL XERRWD(MSG,44,655,0,0,0,0,2,TN,H) |
|
|
2192 MSG = 'DASPK-- PRECONDITIONER HAD REPEATED FAILURES.' |
|
|
2193 CALL XERRWD(MSG,46,656,0,0,0,0,0,0.0D0,0.0D0) |
|
|
2194 GO TO 700 |
|
|
2195 C |
|
|
2196 C The iteration matrix is singular. |
|
|
2197 C |
|
|
2198 660 MSG = 'DASPK-- AT T (=R1) AND STEPSIZE H (=R2) THE' |
|
|
2199 CALL XERRWD(MSG,44,660,0,0,0,0,2,TN,H) |
|
|
2200 MSG = 'DASPK-- ITERATION MATRIX IS SINGULAR.' |
|
|
2201 CALL XERRWD(MSG,38,661,0,0,0,0,0,0.0D0,0.0D0) |
|
|
2202 GO TO 700 |
|
|
2203 C |
|
|
2204 C Nonlinear system failure preceded by error test failures. |
|
|
2205 C |
|
|
2206 670 MSG = 'DASPK-- AT T (=R1) AND STEPSIZE H (=R2) THE' |
|
|
2207 CALL XERRWD(MSG,44,670,0,0,0,0,2,TN,H) |
|
|
2208 MSG = 'DASPK-- NONLINEAR SOLVER COULD NOT CONVERGE.' |
|
|
2209 CALL XERRWD(MSG,45,671,0,0,0,0,0,0.0D0,0.0D0) |
|
|
2210 MSG = 'DASPK-- ALSO, THE ERROR TEST FAILED REPEATEDLY.' |
|
|
2211 CALL XERRWD(MSG,49,672,0,0,0,0,0,0.0D0,0.0D0) |
|
|
2212 GO TO 700 |
|
|
2213 C |
|
|
2214 C Nonlinear system failure because IRES = -1. |
|
|
2215 C |
|
|
2216 675 MSG = 'DASPK-- AT T (=R1) AND STEPSIZE H (=R2) THE' |
|
|
2217 CALL XERRWD(MSG,44,675,0,0,0,0,2,TN,H) |
|
|
2218 MSG = 'DASPK-- NONLINEAR SYSTEM SOLVER COULD NOT CONVERGE' |
|
|
2219 CALL XERRWD(MSG,51,676,0,0,0,0,0,0.0D0,0.0D0) |
|
|
2220 MSG = 'DASPK-- BECAUSE IRES WAS EQUAL TO MINUS ONE' |
|
|
2221 CALL XERRWD(MSG,44,677,0,0,0,0,0,0.0D0,0.0D0) |
|
|
2222 GO TO 700 |
|
|
2223 C |
|
|
2224 C Failure because IRES = -2. |
|
|
2225 C |
|
|
2226 680 MSG = 'DASPK-- AT T (=R1) AND STEPSIZE H (=R2)' |
|
|
2227 CALL XERRWD(MSG,40,680,0,0,0,0,2,TN,H) |
|
|
2228 MSG = 'DASPK-- IRES WAS EQUAL TO MINUS TWO' |
|
|
2229 CALL XERRWD(MSG,36,681,0,0,0,0,0,0.0D0,0.0D0) |
|
|
2230 GO TO 700 |
|
|
2231 C |
|
|
2232 C Failed to compute initial YPRIME. |
|
|
2233 C |
|
|
2234 685 MSG = 'DASPK-- AT T (=R1) AND STEPSIZE H (=R2) THE' |
|
|
2235 CALL XERRWD(MSG,44,685,0,0,0,0,0,0.0D0,0.0D0) |
|
|
2236 MSG = 'DASPK-- INITIAL (Y,YPRIME) COULD NOT BE COMPUTED' |
|
|
2237 CALL XERRWD(MSG,49,686,0,0,0,0,2,TN,H0) |
|
|
2238 GO TO 700 |
|
|
2239 C |
|
|
2240 C Failure because IER was negative from PSOL. |
|
|
2241 C |
|
|
2242 690 MSG = 'DASPK-- AT T (=R1) AND STEPSIZE H (=R2)' |
|
|
2243 CALL XERRWD(MSG,40,690,0,0,0,0,2,TN,H) |
|
|
2244 MSG = 'DASPK-- IER WAS NEGATIVE FROM PSOL' |
|
|
2245 CALL XERRWD(MSG,35,691,0,0,0,0,0,0.0D0,0.0D0) |
|
|
2246 GO TO 700 |
|
|
2247 C |
|
|
2248 C Failure because the linear system solver could not converge. |
|
|
2249 C |
|
|
2250 695 MSG = 'DASPK-- AT T (=R1) AND STEPSIZE H (=R2) THE' |
|
|
2251 CALL XERRWD(MSG,44,695,0,0,0,0,2,TN,H) |
|
|
2252 MSG = 'DASPK-- LINEAR SYSTEM SOLVER COULD NOT CONVERGE.' |
|
|
2253 CALL XERRWD(MSG,50,696,0,0,0,0,0,0.0D0,0.0D0) |
|
|
2254 GO TO 700 |
|
|
2255 C |
|
|
2256 C |
|
|
2257 700 CONTINUE |
|
|
2258 INFO(1)=-1 |
|
|
2259 T=TN |
|
|
2260 RWORK(LTN)=TN |
|
|
2261 RWORK(LH)=H |
|
|
2262 RETURN |
|
|
2263 C |
|
|
2264 C----------------------------------------------------------------------- |
|
|
2265 C This block handles all error returns due to illegal input, |
|
|
2266 C as detected before calling DDSTP. |
|
|
2267 C First the error message routine is called. If this happens |
|
|
2268 C twice in succession, execution is terminated. |
|
|
2269 C----------------------------------------------------------------------- |
|
|
2270 C |
|
|
2271 701 MSG = 'DASPK-- ELEMENT (=I1) OF INFO VECTOR IS NOT VALID' |
|
|
2272 CALL XERRWD(MSG,50,1,0,1,ITEMP,0,0,0.0D0,0.0D0) |
|
|
2273 GO TO 750 |
|
|
2274 702 MSG = 'DASPK-- NEQ (=I1) .LE. 0' |
|
|
2275 CALL XERRWD(MSG,25,2,0,1,NEQ,0,0,0.0D0,0.0D0) |
|
|
2276 GO TO 750 |
|
|
2277 703 MSG = 'DASPK-- MAXORD (=I1) NOT IN RANGE' |
|
|
2278 CALL XERRWD(MSG,34,3,0,1,MXORD,0,0,0.0D0,0.0D0) |
|
|
2279 GO TO 750 |
|
|
2280 704 MSG='DASPK-- RWORK LENGTH NEEDED, LENRW (=I1), EXCEEDS LRW (=I2)' |
|
|
2281 CALL XERRWD(MSG,60,4,0,2,LENRW,LRW,0,0.0D0,0.0D0) |
|
|
2282 GO TO 750 |
|
|
2283 705 MSG='DASPK-- IWORK LENGTH NEEDED, LENIW (=I1), EXCEEDS LIW (=I2)' |
|
|
2284 CALL XERRWD(MSG,60,5,0,2,LENIW,LIW,0,0.0D0,0.0D0) |
|
|
2285 GO TO 750 |
|
|
2286 706 MSG = 'DASPK-- SOME ELEMENT OF RTOL IS .LT. 0' |
|
|
2287 CALL XERRWD(MSG,39,6,0,0,0,0,0,0.0D0,0.0D0) |
|
|
2288 GO TO 750 |
|
|
2289 707 MSG = 'DASPK-- SOME ELEMENT OF ATOL IS .LT. 0' |
|
|
2290 CALL XERRWD(MSG,39,7,0,0,0,0,0,0.0D0,0.0D0) |
|
|
2291 GO TO 750 |
|
|
2292 708 MSG = 'DASPK-- ALL ELEMENTS OF RTOL AND ATOL ARE ZERO' |
|
|
2293 CALL XERRWD(MSG,47,8,0,0,0,0,0,0.0D0,0.0D0) |
|
|
2294 GO TO 750 |
|
|
2295 709 MSG='DASPK-- INFO(4) = 1 AND TSTOP (=R1) BEHIND TOUT (=R2)' |
|
|
2296 CALL XERRWD(MSG,54,9,0,0,0,0,2,TSTOP,TOUT) |
|
|
2297 GO TO 750 |
|
|
2298 710 MSG = 'DASPK-- HMAX (=R1) .LT. 0.0' |
|
|
2299 CALL XERRWD(MSG,28,10,0,0,0,0,1,HMAX,0.0D0) |
|
|
2300 GO TO 750 |
|
|
2301 711 MSG = 'DASPK-- TOUT (=R1) BEHIND T (=R2)' |
|
|
2302 CALL XERRWD(MSG,34,11,0,0,0,0,2,TOUT,T) |
|
|
2303 GO TO 750 |
|
|
2304 712 MSG = 'DASPK-- INFO(8)=1 AND H0=0.0' |
|
|
2305 CALL XERRWD(MSG,29,12,0,0,0,0,0,0.0D0,0.0D0) |
|
|
2306 GO TO 750 |
|
|
2307 713 MSG = 'DASPK-- SOME ELEMENT OF WT IS .LE. 0.0' |
|
|
2308 CALL XERRWD(MSG,39,13,0,0,0,0,0,0.0D0,0.0D0) |
|
|
2309 GO TO 750 |
|
|
2310 714 MSG='DASPK-- TOUT (=R1) TOO CLOSE TO T (=R2) TO START INTEGRATION' |
|
|
2311 CALL XERRWD(MSG,60,14,0,0,0,0,2,TOUT,T) |
|
|
2312 GO TO 750 |
|
|
2313 715 MSG = 'DASPK-- INFO(4)=1 AND TSTOP (=R1) BEHIND T (=R2)' |
|
|
2314 CALL XERRWD(MSG,49,15,0,0,0,0,2,TSTOP,T) |
|
|
2315 GO TO 750 |
|
|
2316 717 MSG = 'DASPK-- ML (=I1) ILLEGAL. EITHER .LT. 0 OR .GT. NEQ' |
|
|
2317 CALL XERRWD(MSG,52,17,0,1,IWORK(LML),0,0,0.0D0,0.0D0) |
|
|
2318 GO TO 750 |
|
|
2319 718 MSG = 'DASPK-- MU (=I1) ILLEGAL. EITHER .LT. 0 OR .GT. NEQ' |
|
|
2320 CALL XERRWD(MSG,52,18,0,1,IWORK(LMU),0,0,0.0D0,0.0D0) |
|
|
2321 GO TO 750 |
|
|
2322 719 MSG = 'DASPK-- TOUT (=R1) IS EQUAL TO T (=R2)' |
|
|
2323 CALL XERRWD(MSG,39,19,0,0,0,0,2,TOUT,T) |
|
|
2324 GO TO 750 |
|
|
2325 720 MSG = 'DASPK-- MAXL (=I1) ILLEGAL. EITHER .LT. 1 OR .GT. NEQ' |
|
|
2326 CALL XERRWD(MSG,54,20,0,1,IWORK(LMAXL),0,0,0.0D0,0.0D0) |
|
|
2327 GO TO 750 |
|
|
2328 721 MSG = 'DASPK-- KMP (=I1) ILLEGAL. EITHER .LT. 1 OR .GT. MAXL' |
|
|
2329 CALL XERRWD(MSG,54,21,0,1,IWORK(LKMP),0,0,0.0D0,0.0D0) |
|
|
2330 GO TO 750 |
|
|
2331 722 MSG = 'DASPK-- NRMAX (=I1) ILLEGAL. .LT. 0' |
|
|
2332 CALL XERRWD(MSG,36,22,0,1,IWORK(LNRMAX),0,0,0.0D0,0.0D0) |
|
|
2333 GO TO 750 |
|
|
2334 723 MSG = 'DASPK-- EPLI (=R1) ILLEGAL. EITHER .LE. 0.D0 OR .GE. 1.D0' |
|
|
2335 CALL XERRWD(MSG,58,23,0,0,0,0,1,RWORK(LEPLI),0.0D0) |
|
|
2336 GO TO 750 |
|
|
2337 724 MSG = 'DASPK-- ILLEGAL IWORK VALUE FOR INFO(11) .NE. 0' |
|
|
2338 CALL XERRWD(MSG,48,24,0,0,0,0,0,0.0D0,0.0D0) |
|
|
2339 GO TO 750 |
|
|
2340 725 MSG = 'DASPK-- ONE OF THE INPUTS FOR INFO(17) = 1 IS ILLEGAL' |
|
|
2341 CALL XERRWD(MSG,54,25,0,0,0,0,0,0.0D0,0.0D0) |
|
|
2342 GO TO 750 |
|
|
2343 726 MSG = 'DASPK-- ILLEGAL IWORK VALUE FOR INFO(10) .NE. 0' |
|
|
2344 CALL XERRWD(MSG,48,26,0,0,0,0,0,0.0D0,0.0D0) |
|
|
2345 GO TO 750 |
|
|
2346 727 MSG = 'DASPK-- Y(I) AND IWORK(40+I) (I=I1) INCONSISTENT' |
|
|
2347 CALL XERRWD(MSG,49,27,0,1,IRET,0,0,0.0D0,0.0D0) |
|
|
2348 GO TO 750 |
|
|
2349 750 IF(INFO(1).EQ.-1) GO TO 760 |
|
|
2350 INFO(1)=-1 |
|
|
2351 IDID=-33 |
|
|
2352 RETURN |
|
|
2353 760 MSG = 'DASPK-- REPEATED OCCURRENCES OF ILLEGAL INPUT' |
|
|
2354 CALL XERRWD(MSG,46,701,0,0,0,0,0,0.0D0,0.0D0) |
|
|
2355 770 MSG = 'DASPK-- RUN TERMINATED. APPARENT INFINITE LOOP' |
|
|
2356 CALL XERRWD(MSG,47,702,1,0,0,0,0,0.0D0,0.0D0) |
|
|
2357 RETURN |
|
|
2358 C |
|
|
2359 C------END OF SUBROUTINE DDASPK----------------------------------------- |
|
|
2360 END |
