2329
|
1 SUBROUTINE DDASSL (RES, NEQ, T, Y, YPRIME, TOUT, INFO, RTOL, ATOL, |
|
2 + IDID, RWORK, LRW, IWORK, LIW, RPAR, IPAR, JAC) |
|
3 C***BEGIN PROLOGUE DDASSL |
|
4 C***PURPOSE This code solves a system of differential/algebraic |
|
5 C equations of the form G(T,Y,YPRIME) = 0. |
|
6 C***LIBRARY SLATEC (DASSL) |
|
7 C***CATEGORY I1A2 |
|
8 C***TYPE DOUBLE PRECISION (SDASSL-S, DDASSL-D) |
|
9 C***KEYWORDS DIFFERENTIAL/ALGEBRAIC, BACKWARD DIFFERENTIATION FORMULAS, |
|
10 C IMPLICIT DIFFERENTIAL SYSTEMS |
|
11 C***AUTHOR PETZOLD, LINDA R., (LLNL) |
|
12 C COMPUTING AND MATHEMATICS RESEARCH DIVISION |
|
13 C LAWRENCE LIVERMORE NATIONAL LABORATORY |
|
14 C L - 316, P.O. BOX 808, |
|
15 C LIVERMORE, CA. 94550 |
|
16 C***DESCRIPTION |
|
17 C |
|
18 C *Usage: |
|
19 C |
|
20 C EXTERNAL RES, JAC |
|
21 C INTEGER NEQ, INFO(N), IDID, LRW, LIW, IWORK(LIW), IPAR |
|
22 C DOUBLE PRECISION T, Y(NEQ), YPRIME(NEQ), TOUT, RTOL, ATOL, |
|
23 C * RWORK(LRW), RPAR |
|
24 C |
|
25 C CALL DDASSL (RES, NEQ, T, Y, YPRIME, TOUT, INFO, RTOL, ATOL, |
|
26 C * IDID, RWORK, LRW, IWORK, LIW, RPAR, IPAR, JAC) |
|
27 C |
|
28 C |
|
29 C *Arguments: |
|
30 C (In the following, all real arrays should be type DOUBLE PRECISION.) |
|
31 C |
|
32 C RES:EXT This is a subroutine which you provide to define the |
|
33 C differential/algebraic system. |
|
34 C |
|
35 C NEQ:IN This is the number of equations to be solved. |
|
36 C |
|
37 C T:INOUT This is the current value of the independent variable. |
|
38 C |
|
39 C Y(*):INOUT This array contains the solution components at T. |
|
40 C |
|
41 C YPRIME(*):INOUT This array contains the derivatives of the solution |
|
42 C components at T. |
|
43 C |
|
44 C TOUT:IN This is a point at which a solution is desired. |
|
45 C |
|
46 C INFO(N):IN The basic task of the code is to solve the system from T |
|
47 C to TOUT and return an answer at TOUT. INFO is an integer |
|
48 C array which is used to communicate exactly how you want |
|
49 C this task to be carried out. (See below for details.) |
|
50 C N must be greater than or equal to 15. |
|
51 C |
|
52 C RTOL,ATOL:INOUT These quantities represent relative and absolute |
|
53 C error tolerances which you provide to indicate how |
|
54 C accurately you wish the solution to be computed. You |
|
55 C may choose them to be both scalars or else both vectors. |
|
56 C Caution: In Fortran 77, a scalar is not the same as an |
|
57 C array of length 1. Some compilers may object |
|
58 C to using scalars for RTOL,ATOL. |
|
59 C |
|
60 C IDID:OUT This scalar quantity is an indicator reporting what the |
|
61 C code did. You must monitor this integer variable to |
|
62 C decide what action to take next. |
|
63 C |
|
64 C RWORK:WORK A real work array of length LRW which provides the |
|
65 C code with needed storage space. |
|
66 C |
|
67 C LRW:IN The length of RWORK. (See below for required length.) |
|
68 C |
|
69 C IWORK:WORK An integer work array of length LIW which probides the |
|
70 C code with needed storage space. |
|
71 C |
|
72 C LIW:IN The length of IWORK. (See below for required length.) |
|
73 C |
|
74 C RPAR,IPAR:IN These are real and integer parameter arrays which |
|
75 C you can use for communication between your calling |
|
76 C program and the RES subroutine (and the JAC subroutine) |
|
77 C |
|
78 C JAC:EXT This is the name of a subroutine which you may choose |
|
79 C to provide for defining a matrix of partial derivatives |
|
80 C described below. |
|
81 C |
|
82 C Quantities which may be altered by DDASSL are: |
|
83 C T, Y(*), YPRIME(*), INFO(1), RTOL, ATOL, |
|
84 C IDID, RWORK(*) AND IWORK(*) |
|
85 C |
|
86 C *Description |
|
87 C |
|
88 C Subroutine DDASSL uses the backward differentiation formulas of |
|
89 C orders one through five to solve a system of the above form for Y and |
|
90 C YPRIME. Values for Y and YPRIME at the initial time must be given as |
|
91 C input. These values must be consistent, (that is, if T,Y,YPRIME are |
|
92 C the given initial values, they must satisfy G(T,Y,YPRIME) = 0.). The |
|
93 C subroutine solves the system from T to TOUT. It is easy to continue |
|
94 C the solution to get results at additional TOUT. This is the interval |
|
95 C mode of operation. Intermediate results can also be obtained easily |
|
96 C by using the intermediate-output capability. |
|
97 C |
|
98 C The following detailed description is divided into subsections: |
|
99 C 1. Input required for the first call to DDASSL. |
|
100 C 2. Output after any return from DDASSL. |
|
101 C 3. What to do to continue the integration. |
|
102 C 4. Error messages. |
|
103 C |
|
104 C |
|
105 C -------- INPUT -- WHAT TO DO ON THE FIRST CALL TO DDASSL ------------ |
|
106 C |
|
107 C The first call of the code is defined to be the start of each new |
|
108 C problem. Read through the descriptions of all the following items, |
|
109 C provide sufficient storage space for designated arrays, set |
|
110 C appropriate variables for the initialization of the problem, and |
|
111 C give information about how you want the problem to be solved. |
|
112 C |
|
113 C |
|
114 C RES -- Provide a subroutine of the form |
|
115 C SUBROUTINE RES(T,Y,YPRIME,DELTA,IRES,RPAR,IPAR) |
|
116 C to define the system of differential/algebraic |
|
117 C equations which is to be solved. For the given values |
|
118 C of T,Y and YPRIME, the subroutine should |
|
119 C return the residual of the defferential/algebraic |
|
120 C system |
|
121 C DELTA = G(T,Y,YPRIME) |
|
122 C (DELTA(*) is a vector of length NEQ which is |
|
123 C output for RES.) |
|
124 C |
|
125 C Subroutine RES must not alter T,Y or YPRIME. |
|
126 C You must declare the name RES in an external |
|
127 C statement in your program that calls DDASSL. |
|
128 C You must dimension Y,YPRIME and DELTA in RES. |
|
129 C |
|
130 C IRES is an integer flag which is always equal to |
|
131 C zero on input. Subroutine RES should alter IRES |
|
132 C only if it encounters an illegal value of Y or |
|
133 C a stop condition. Set IRES = -1 if an input value |
|
134 C is illegal, and DDASSL will try to solve the problem |
|
135 C without getting IRES = -1. If IRES = -2, DDASSL |
|
136 C will return control to the calling program |
|
137 C with IDID = -11. |
|
138 C |
|
139 C RPAR and IPAR are real and integer parameter arrays which |
|
140 C you can use for communication between your calling program |
|
141 C and subroutine RES. They are not altered by DDASSL. If you |
|
142 C do not need RPAR or IPAR, ignore these parameters by treat- |
|
143 C ing them as dummy arguments. If you do choose to use them, |
|
144 C dimension them in your calling program and in RES as arrays |
|
145 C of appropriate length. |
|
146 C |
|
147 C NEQ -- Set it to the number of differential equations. |
|
148 C (NEQ .GE. 1) |
|
149 C |
|
150 C T -- Set it to the initial point of the integration. |
|
151 C T must be defined as a variable. |
|
152 C |
|
153 C Y(*) -- Set this vector to the initial values of the NEQ solution |
|
154 C components at the initial point. You must dimension Y of |
|
155 C length at least NEQ in your calling program. |
|
156 C |
|
157 C YPRIME(*) -- Set this vector to the initial values of the NEQ |
|
158 C first derivatives of the solution components at the initial |
|
159 C point. You must dimension YPRIME at least NEQ in your |
|
160 C calling program. If you do not know initial values of some |
|
161 C of the solution components, see the explanation of INFO(11). |
|
162 C |
|
163 C TOUT -- Set it to the first point at which a solution |
|
164 C is desired. You can not take TOUT = T. |
|
165 C integration either forward in T (TOUT .GT. T) or |
|
166 C backward in T (TOUT .LT. T) is permitted. |
|
167 C |
|
168 C The code advances the solution from T to TOUT using |
|
169 C step sizes which are automatically selected so as to |
|
170 C achieve the desired accuracy. If you wish, the code will |
|
171 C return with the solution and its derivative at |
|
172 C intermediate steps (intermediate-output mode) so that |
|
173 C you can monitor them, but you still must provide TOUT in |
|
174 C accord with the basic aim of the code. |
|
175 C |
|
176 C The first step taken by the code is a critical one |
|
177 C because it must reflect how fast the solution changes near |
|
178 C the initial point. The code automatically selects an |
|
179 C initial step size which is practically always suitable for |
|
180 C the problem. By using the fact that the code will not step |
|
181 C past TOUT in the first step, you could, if necessary, |
|
182 C restrict the length of the initial step size. |
|
183 C |
|
184 C For some problems it may not be permissible to integrate |
|
185 C past a point TSTOP because a discontinuity occurs there |
|
186 C or the solution or its derivative is not defined beyond |
|
187 C TSTOP. When you have declared a TSTOP point (SEE INFO(4) |
|
188 C and RWORK(1)), you have told the code not to integrate |
|
189 C past TSTOP. In this case any TOUT beyond TSTOP is invalid |
|
190 C input. |
|
191 C |
|
192 C INFO(*) -- Use the INFO array to give the code more details about |
|
193 C how you want your problem solved. This array should be |
|
194 C dimensioned of length 15, though DDASSL uses only the first |
|
195 C eleven entries. You must respond to all of the following |
|
196 C items, which are arranged as questions. The simplest use |
|
197 C of the code corresponds to answering all questions as yes, |
|
198 C i.e. setting all entries of INFO to 0. |
|
199 C |
|
200 C INFO(1) - This parameter enables the code to initialize |
|
201 C itself. You must set it to indicate the start of every |
|
202 C new problem. |
|
203 C |
|
204 C **** Is this the first call for this problem ... |
|
205 C Yes - Set INFO(1) = 0 |
|
206 C No - Not applicable here. |
|
207 C See below for continuation calls. **** |
|
208 C |
|
209 C INFO(2) - How much accuracy you want of your solution |
|
210 C is specified by the error tolerances RTOL and ATOL. |
|
211 C The simplest use is to take them both to be scalars. |
|
212 C To obtain more flexibility, they can both be vectors. |
|
213 C The code must be told your choice. |
|
214 C |
|
215 C **** Are both error tolerances RTOL, ATOL scalars ... |
|
216 C Yes - Set INFO(2) = 0 |
|
217 C and input scalars for both RTOL and ATOL |
|
218 C No - Set INFO(2) = 1 |
|
219 C and input arrays for both RTOL and ATOL **** |
|
220 C |
|
221 C INFO(3) - The code integrates from T in the direction |
|
222 C of TOUT by steps. If you wish, it will return the |
|
223 C computed solution and derivative at the next |
|
224 C intermediate step (the intermediate-output mode) or |
|
225 C TOUT, whichever comes first. This is a good way to |
|
226 C proceed if you want to see the behavior of the solution. |
|
227 C If you must have solutions at a great many specific |
|
228 C TOUT points, this code will compute them efficiently. |
|
229 C |
|
230 C **** Do you want the solution only at |
|
231 C TOUT (and not at the next intermediate step) ... |
|
232 C Yes - Set INFO(3) = 0 |
|
233 C No - Set INFO(3) = 1 **** |
|
234 C |
|
235 C INFO(4) - To handle solutions at a great many specific |
|
236 C values TOUT efficiently, this code may integrate past |
|
237 C TOUT and interpolate to obtain the result at TOUT. |
|
238 C Sometimes it is not possible to integrate beyond some |
|
239 C point TSTOP because the equation changes there or it is |
|
240 C not defined past TSTOP. Then you must tell the code |
|
241 C not to go past. |
|
242 C |
|
243 C **** Can the integration be carried out without any |
|
244 C restrictions on the independent variable T ... |
|
245 C Yes - Set INFO(4)=0 |
|
246 C No - Set INFO(4)=1 |
|
247 C and define the stopping point TSTOP by |
|
248 C setting RWORK(1)=TSTOP **** |
|
249 C |
|
250 C INFO(5) - To solve differential/algebraic problems it is |
|
251 C necessary to use a matrix of partial derivatives of the |
|
252 C system of differential equations. If you do not |
|
253 C provide a subroutine to evaluate it analytically (see |
|
254 C description of the item JAC in the call list), it will |
|
255 C be approximated by numerical differencing in this code. |
|
256 C although it is less trouble for you to have the code |
|
257 C compute partial derivatives by numerical differencing, |
|
258 C the solution will be more reliable if you provide the |
|
259 C derivatives via JAC. Sometimes numerical differencing |
|
260 C is cheaper than evaluating derivatives in JAC and |
|
261 C sometimes it is not - this depends on your problem. |
|
262 C |
|
263 C **** Do you want the code to evaluate the partial |
|
264 C derivatives automatically by numerical differences ... |
|
265 C Yes - Set INFO(5)=0 |
|
266 C No - Set INFO(5)=1 |
|
267 C and provide subroutine JAC for evaluating the |
|
268 C matrix of partial derivatives **** |
|
269 C |
|
270 C INFO(6) - DDASSL will perform much better if the matrix of |
|
271 C partial derivatives, DG/DY + CJ*DG/DYPRIME, |
|
272 C (here CJ is a scalar determined by DDASSL) |
|
273 C is banded and the code is told this. In this |
|
274 C case, the storage needed will be greatly reduced, |
|
275 C numerical differencing will be performed much cheaper, |
|
276 C and a number of important algorithms will execute much |
|
277 C faster. The differential equation is said to have |
|
278 C half-bandwidths ML (lower) and MU (upper) if equation i |
|
279 C involves only unknowns Y(J) with |
|
280 C I-ML .LE. J .LE. I+MU |
|
281 C for all I=1,2,...,NEQ. Thus, ML and MU are the widths |
|
282 C of the lower and upper parts of the band, respectively, |
|
283 C with the main diagonal being excluded. If you do not |
|
284 C indicate that the equation has a banded matrix of partial |
|
285 C derivatives, the code works with a full matrix of NEQ**2 |
|
286 C elements (stored in the conventional way). Computations |
|
287 C with banded matrices cost less time and storage than with |
|
288 C full matrices if 2*ML+MU .LT. NEQ. If you tell the |
|
289 C code that the matrix of partial derivatives has a banded |
|
290 C structure and you want to provide subroutine JAC to |
|
291 C compute the partial derivatives, then you must be careful |
|
292 C to store the elements of the matrix in the special form |
|
293 C indicated in the description of JAC. |
|
294 C |
|
295 C **** Do you want to solve the problem using a full |
|
296 C (dense) matrix (and not a special banded |
|
297 C structure) ... |
|
298 C Yes - Set INFO(6)=0 |
|
299 C No - Set INFO(6)=1 |
|
300 C and provide the lower (ML) and upper (MU) |
|
301 C bandwidths by setting |
|
302 C IWORK(1)=ML |
|
303 C IWORK(2)=MU **** |
|
304 C |
|
305 C |
|
306 C INFO(7) -- You can specify a maximum (absolute value of) |
|
307 C stepsize, so that the code |
|
308 C will avoid passing over very |
|
309 C large regions. |
|
310 C |
|
311 C **** Do you want the code to decide |
|
312 C on its own maximum stepsize? |
|
313 C Yes - Set INFO(7)=0 |
|
314 C No - Set INFO(7)=1 |
|
315 C and define HMAX by setting |
|
316 C RWORK(2)=HMAX **** |
|
317 C |
|
318 C INFO(8) -- Differential/algebraic problems |
|
319 C may occaisionally suffer from |
|
320 C severe scaling difficulties on the |
|
321 C first step. If you know a great deal |
|
322 C about the scaling of your problem, you can |
|
323 C help to alleviate this problem by |
|
324 C specifying an initial stepsize HO. |
|
325 C |
|
326 C **** Do you want the code to define |
|
327 C its own initial stepsize? |
|
328 C Yes - Set INFO(8)=0 |
|
329 C No - Set INFO(8)=1 |
|
330 C and define HO by setting |
|
331 C RWORK(3)=HO **** |
|
332 C |
|
333 C INFO(9) -- If storage is a severe problem, |
|
334 C you can save some locations by |
|
335 C restricting the maximum order MAXORD. |
|
336 C the default value is 5. for each |
|
337 C order decrease below 5, the code |
|
338 C requires NEQ fewer locations, however |
|
339 C it is likely to be slower. In any |
|
340 C case, you must have 1 .LE. MAXORD .LE. 5 |
|
341 C **** Do you want the maximum order to |
|
342 C default to 5? |
|
343 C Yes - Set INFO(9)=0 |
|
344 C No - Set INFO(9)=1 |
|
345 C and define MAXORD by setting |
|
346 C IWORK(3)=MAXORD **** |
|
347 C |
|
348 C INFO(10) --If you know that the solutions to your equations |
|
349 C will always be nonnegative, it may help to set this |
|
350 C parameter. However, it is probably best to |
|
351 C try the code without using this option first, |
|
352 C and only to use this option if that doesn't |
|
353 C work very well. |
|
354 C **** Do you want the code to solve the problem without |
|
355 C invoking any special nonnegativity constraints? |
|
356 C Yes - Set INFO(10)=0 |
|
357 C No - Set INFO(10)=1 |
|
358 C |
|
359 C INFO(11) --DDASSL normally requires the initial T, |
|
360 C Y, and YPRIME to be consistent. That is, |
|
361 C you must have G(T,Y,YPRIME) = 0 at the initial |
|
362 C time. If you do not know the initial |
|
363 C derivative precisely, you can let DDASSL try |
|
364 C to compute it. |
|
365 C **** Are the initialHE INITIAL T, Y, YPRIME consistent? |
|
366 C Yes - Set INFO(11) = 0 |
|
367 C No - Set INFO(11) = 1, |
|
368 C and set YPRIME to an initial approximation |
|
369 C to YPRIME. (If you have no idea what |
|
370 C YPRIME should be, set it to zero. Note |
|
371 C that the initial Y should be such |
|
372 C that there must exist a YPRIME so that |
|
373 C G(T,Y,YPRIME) = 0.) |
|
374 C |
|
375 C RTOL, ATOL -- You must assign relative (RTOL) and absolute (ATOL |
|
376 C error tolerances to tell the code how accurately you |
|
377 C want the solution to be computed. They must be defined |
|
378 C as variables because the code may change them. You |
|
379 C have two choices -- |
|
380 C Both RTOL and ATOL are scalars. (INFO(2)=0) |
|
381 C Both RTOL and ATOL are vectors. (INFO(2)=1) |
|
382 C in either case all components must be non-negative. |
|
383 C |
|
384 C The tolerances are used by the code in a local error |
|
385 C test at each step which requires roughly that |
|
386 C ABS(LOCAL ERROR) .LE. RTOL*ABS(Y)+ATOL |
|
387 C for each vector component. |
|
388 C (More specifically, a root-mean-square norm is used to |
|
389 C measure the size of vectors, and the error test uses the |
|
390 C magnitude of the solution at the beginning of the step.) |
|
391 C |
|
392 C The true (global) error is the difference between the |
|
393 C true solution of the initial value problem and the |
|
394 C computed approximation. Practically all present day |
|
395 C codes, including this one, control the local error at |
|
396 C each step and do not even attempt to control the global |
|
397 C error directly. |
|
398 C Usually, but not always, the true accuracy of the |
|
399 C computed Y is comparable to the error tolerances. This |
|
400 C code will usually, but not always, deliver a more |
|
401 C accurate solution if you reduce the tolerances and |
|
402 C integrate again. By comparing two such solutions you |
|
403 C can get a fairly reliable idea of the true error in the |
|
404 C solution at the bigger tolerances. |
|
405 C |
|
406 C Setting ATOL=0. results in a pure relative error test on |
|
407 C that component. Setting RTOL=0. results in a pure |
|
408 C absolute error test on that component. A mixed test |
|
409 C with non-zero RTOL and ATOL corresponds roughly to a |
|
410 C relative error test when the solution component is much |
|
411 C bigger than ATOL and to an absolute error test when the |
|
412 C solution component is smaller than the threshhold ATOL. |
|
413 C |
|
414 C The code will not attempt to compute a solution at an |
|
415 C accuracy unreasonable for the machine being used. It will |
|
416 C advise you if you ask for too much accuracy and inform |
|
417 C you as to the maximum accuracy it believes possible. |
|
418 C |
|
419 C RWORK(*) -- Dimension this real work array of length LRW in your |
|
420 C calling program. |
|
421 C |
|
422 C LRW -- Set it to the declared length of the RWORK array. |
|
423 C You must have |
|
424 C LRW .GE. 40+(MAXORD+4)*NEQ+NEQ**2 |
|
425 C for the full (dense) JACOBIAN case (when INFO(6)=0), or |
|
426 C LRW .GE. 40+(MAXORD+4)*NEQ+(2*ML+MU+1)*NEQ |
|
427 C for the banded user-defined JACOBIAN case |
|
428 C (when INFO(5)=1 and INFO(6)=1), or |
|
429 C LRW .GE. 40+(MAXORD+4)*NEQ+(2*ML+MU+1)*NEQ |
|
430 C +2*(NEQ/(ML+MU+1)+1) |
|
431 C for the banded finite-difference-generated JACOBIAN case |
|
432 C (when INFO(5)=0 and INFO(6)=1) |
|
433 C |
|
434 C IWORK(*) -- Dimension this integer work array of length LIW in |
|
435 C your calling program. |
|
436 C |
|
437 C LIW -- Set it to the declared length of the IWORK array. |
|
438 C You must have LIW .GE. 20+NEQ |
|
439 C |
|
440 C RPAR, IPAR -- These are parameter arrays, of real and integer |
|
441 C type, respectively. You can use them for communication |
|
442 C between your program that calls DDASSL and the |
|
443 C RES subroutine (and the JAC subroutine). They are not |
|
444 C altered by DDASSL. If you do not need RPAR or IPAR, |
|
445 C ignore these parameters by treating them as dummy |
|
446 C arguments. If you do choose to use them, dimension |
|
447 C them in your calling program and in RES (and in JAC) |
|
448 C as arrays of appropriate length. |
|
449 C |
|
450 C JAC -- If you have set INFO(5)=0, you can ignore this parameter |
|
451 C by treating it as a dummy argument. Otherwise, you must |
|
452 C provide a subroutine of the form |
|
453 C SUBROUTINE JAC(T,Y,YPRIME,PD,CJ,RPAR,IPAR) |
|
454 C to define the matrix of partial derivatives |
|
455 C PD=DG/DY+CJ*DG/DYPRIME |
|
456 C CJ is a scalar which is input to JAC. |
|
457 C For the given values of T,Y,YPRIME, the |
|
458 C subroutine must evaluate the non-zero partial |
|
459 C derivatives for each equation and each solution |
|
460 C component, and store these values in the |
|
461 C matrix PD. The elements of PD are set to zero |
|
462 C before each call to JAC so only non-zero elements |
|
463 C need to be defined. |
|
464 C |
|
465 C Subroutine JAC must not alter T,Y,(*),YPRIME(*), or CJ. |
|
466 C You must declare the name JAC in an EXTERNAL statement in |
|
467 C your program that calls DDASSL. You must dimension Y, |
|
468 C YPRIME and PD in JAC. |
|
469 C |
|
470 C The way you must store the elements into the PD matrix |
|
471 C depends on the structure of the matrix which you |
|
472 C indicated by INFO(6). |
|
473 C *** INFO(6)=0 -- Full (dense) matrix *** |
|
474 C Give PD a first dimension of NEQ. |
|
475 C When you evaluate the (non-zero) partial derivative |
|
476 C of equation I with respect to variable J, you must |
|
477 C store it in PD according to |
|
478 C PD(I,J) = "DG(I)/DY(J)+CJ*DG(I)/DYPRIME(J)" |
|
479 C *** INFO(6)=1 -- Banded JACOBIAN with ML lower and MU |
|
480 C upper diagonal bands (refer to INFO(6) description |
|
481 C of ML and MU) *** |
|
482 C Give PD a first dimension of 2*ML+MU+1. |
|
483 C when you evaluate the (non-zero) partial derivative |
|
484 C of equation I with respect to variable J, you must |
|
485 C store it in PD according to |
|
486 C IROW = I - J + ML + MU + 1 |
|
487 C PD(IROW,J) = "DG(I)/DY(J)+CJ*DG(I)/DYPRIME(J)" |
|
488 C |
|
489 C RPAR and IPAR are real and integer parameter arrays |
|
490 C which you can use for communication between your calling |
|
491 C program and your JACOBIAN subroutine JAC. They are not |
|
492 C altered by DDASSL. If you do not need RPAR or IPAR, |
|
493 C ignore these parameters by treating them as dummy |
|
494 C arguments. If you do choose to use them, dimension |
|
495 C them in your calling program and in JAC as arrays of |
|
496 C appropriate length. |
|
497 C |
|
498 C |
|
499 C OPTIONALLY REPLACEABLE NORM ROUTINE: |
|
500 C |
|
501 C DDASSL uses a weighted norm DDANRM to measure the size |
|
502 C of vectors such as the estimated error in each step. |
|
503 C A FUNCTION subprogram |
|
504 C DOUBLE PRECISION FUNCTION DDANRM(NEQ,V,WT,RPAR,IPAR) |
|
505 C DIMENSION V(NEQ),WT(NEQ) |
|
506 C is used to define this norm. Here, V is the vector |
|
507 C whose norm is to be computed, and WT is a vector of |
|
508 C weights. A DDANRM routine has been included with DDASSL |
|
509 C which computes the weighted root-mean-square norm |
|
510 C given by |
|
511 C DDANRM=SQRT((1/NEQ)*SUM(V(I)/WT(I))**2) |
|
512 C this norm is suitable for most problems. In some |
|
513 C special cases, it may be more convenient and/or |
|
514 C efficient to define your own norm by writing a function |
|
515 C subprogram to be called instead of DDANRM. This should, |
|
516 C however, be attempted only after careful thought and |
|
517 C consideration. |
|
518 C |
|
519 C |
|
520 C -------- OUTPUT -- AFTER ANY RETURN FROM DDASSL --------------------- |
|
521 C |
|
522 C The principal aim of the code is to return a computed solution at |
|
523 C TOUT, although it is also possible to obtain intermediate results |
|
524 C along the way. To find out whether the code achieved its goal |
|
525 C or if the integration process was interrupted before the task was |
|
526 C completed, you must check the IDID parameter. |
|
527 C |
|
528 C |
|
529 C T -- The solution was successfully advanced to the |
|
530 C output value of T. |
|
531 C |
|
532 C Y(*) -- Contains the computed solution approximation at T. |
|
533 C |
|
534 C YPRIME(*) -- Contains the computed derivative |
|
535 C approximation at T. |
|
536 C |
|
537 C IDID -- Reports what the code did. |
|
538 C |
|
539 C *** Task completed *** |
|
540 C Reported by positive values of IDID |
|
541 C |
|
542 C IDID = 1 -- A step was successfully taken in the |
|
543 C intermediate-output mode. The code has not |
|
544 C yet reached TOUT. |
|
545 C |
|
546 C IDID = 2 -- The integration to TSTOP was successfully |
|
547 C completed (T=TSTOP) by stepping exactly to TSTOP. |
|
548 C |
|
549 C IDID = 3 -- The integration to TOUT was successfully |
|
550 C completed (T=TOUT) by stepping past TOUT. |
|
551 C Y(*) is obtained by interpolation. |
|
552 C YPRIME(*) is obtained by interpolation. |
|
553 C |
|
554 C *** Task interrupted *** |
|
555 C Reported by negative values of IDID |
|
556 C |
|
557 C IDID = -1 -- A large amount of work has been expended. |
|
558 C (About 500 steps) |
|
559 C |
|
560 C IDID = -2 -- The error tolerances are too stringent. |
|
561 C |
|
562 C IDID = -3 -- The local error test cannot be satisfied |
|
563 C because you specified a zero component in ATOL |
|
564 C and the corresponding computed solution |
|
565 C component is zero. Thus, a pure relative error |
|
566 C test is impossible for this component. |
|
567 C |
|
568 C IDID = -6 -- DDASSL had repeated error test |
|
569 C failures on the last attempted step. |
|
570 C |
|
571 C IDID = -7 -- The corrector could not converge. |
|
572 C |
|
573 C IDID = -8 -- The matrix of partial derivatives |
|
574 C is singular. |
|
575 C |
|
576 C IDID = -9 -- The corrector could not converge. |
|
577 C there were repeated error test failures |
|
578 C in this step. |
|
579 C |
|
580 C IDID =-10 -- The corrector could not converge |
|
581 C because IRES was equal to minus one. |
|
582 C |
|
583 C IDID =-11 -- IRES equal to -2 was encountered |
|
584 C and control is being returned to the |
|
585 C calling program. |
|
586 C |
|
587 C IDID =-12 -- DDASSL failed to compute the initial |
|
588 C YPRIME. |
|
589 C |
|
590 C |
|
591 C |
|
592 C IDID = -13,..,-32 -- Not applicable for this code |
|
593 C |
|
594 C *** Task terminated *** |
|
595 C Reported by the value of IDID=-33 |
|
596 C |
|
597 C IDID = -33 -- The code has encountered trouble from which |
|
598 C it cannot recover. A message is printed |
|
599 C explaining the trouble and control is returned |
|
600 C to the calling program. For example, this occurs |
|
601 C when invalid input is detected. |
|
602 C |
|
603 C RTOL, ATOL -- These quantities remain unchanged except when |
|
604 C IDID = -2. In this case, the error tolerances have been |
|
605 C increased by the code to values which are estimated to |
|
606 C be appropriate for continuing the integration. However, |
|
607 C the reported solution at T was obtained using the input |
|
608 C values of RTOL and ATOL. |
|
609 C |
|
610 C RWORK, IWORK -- Contain information which is usually of no |
|
611 C interest to the user but necessary for subsequent calls. |
|
612 C However, you may find use for |
|
613 C |
|
614 C RWORK(3)--Which contains the step size H to be |
|
615 C attempted on the next step. |
|
616 C |
|
617 C RWORK(4)--Which contains the current value of the |
|
618 C independent variable, i.e., the farthest point |
|
619 C integration has reached. This will be different |
|
620 C from T only when interpolation has been |
|
621 C performed (IDID=3). |
|
622 C |
|
623 C RWORK(7)--Which contains the stepsize used |
|
624 C on the last successful step. |
|
625 C |
|
626 C IWORK(7)--Which contains the order of the method to |
|
627 C be attempted on the next step. |
|
628 C |
|
629 C IWORK(8)--Which contains the order of the method used |
|
630 C on the last step. |
|
631 C |
|
632 C IWORK(11)--Which contains the number of steps taken so |
|
633 C far. |
|
634 C |
|
635 C IWORK(12)--Which contains the number of calls to RES |
|
636 C so far. |
|
637 C |
|
638 C IWORK(13)--Which contains the number of evaluations of |
|
639 C the matrix of partial derivatives needed so |
|
640 C far. |
|
641 C |
|
642 C IWORK(14)--Which contains the total number |
|
643 C of error test failures so far. |
|
644 C |
|
645 C IWORK(15)--Which contains the total number |
|
646 C of convergence test failures so far. |
|
647 C (includes singular iteration matrix |
|
648 C failures.) |
|
649 C |
|
650 C |
|
651 C -------- INPUT -- WHAT TO DO TO CONTINUE THE INTEGRATION ------------ |
|
652 C (CALLS AFTER THE FIRST) |
|
653 C |
|
654 C This code is organized so that subsequent calls to continue the |
|
655 C integration involve little (if any) additional effort on your |
|
656 C part. You must monitor the IDID parameter in order to determine |
|
657 C what to do next. |
|
658 C |
|
659 C Recalling that the principal task of the code is to integrate |
|
660 C from T to TOUT (the interval mode), usually all you will need |
|
661 C to do is specify a new TOUT upon reaching the current TOUT. |
|
662 C |
|
663 C Do not alter any quantity not specifically permitted below, |
|
664 C in particular do not alter NEQ,T,Y(*),YPRIME(*),RWORK(*),IWORK(*) |
|
665 C or the differential equation in subroutine RES. Any such |
|
666 C alteration constitutes a new problem and must be treated as such, |
|
667 C i.e., you must start afresh. |
|
668 C |
|
669 C You cannot change from vector to scalar error control or vice |
|
670 C versa (INFO(2)), but you can change the size of the entries of |
|
671 C RTOL, ATOL. Increasing a tolerance makes the equation easier |
|
672 C to integrate. Decreasing a tolerance will make the equation |
|
673 C harder to integrate and should generally be avoided. |
|
674 C |
|
675 C You can switch from the intermediate-output mode to the |
|
676 C interval mode (INFO(3)) or vice versa at any time. |
|
677 C |
|
678 C If it has been necessary to prevent the integration from going |
|
679 C past a point TSTOP (INFO(4), RWORK(1)), keep in mind that the |
|
680 C code will not integrate to any TOUT beyond the currently |
|
681 C specified TSTOP. Once TSTOP has been reached you must change |
|
682 C the value of TSTOP or set INFO(4)=0. You may change INFO(4) |
|
683 C or TSTOP at any time but you must supply the value of TSTOP in |
|
684 C RWORK(1) whenever you set INFO(4)=1. |
|
685 C |
|
686 C Do not change INFO(5), INFO(6), IWORK(1), or IWORK(2) |
|
687 C unless you are going to restart the code. |
|
688 C |
|
689 C *** Following a completed task *** |
|
690 C If |
|
691 C IDID = 1, call the code again to continue the integration |
|
692 C another step in the direction of TOUT. |
|
693 C |
|
694 C IDID = 2 or 3, define a new TOUT and call the code again. |
|
695 C TOUT must be different from T. You cannot change |
|
696 C the direction of integration without restarting. |
|
697 C |
|
698 C *** Following an interrupted task *** |
|
699 C To show the code that you realize the task was |
|
700 C interrupted and that you want to continue, you |
|
701 C must take appropriate action and set INFO(1) = 1 |
|
702 C If |
|
703 C IDID = -1, The code has taken about 500 steps. |
|
704 C If you want to continue, set INFO(1) = 1 and |
|
705 C call the code again. An additional 500 steps |
|
706 C will be allowed. |
|
707 C |
|
708 C IDID = -2, The error tolerances RTOL, ATOL have been |
|
709 C increased to values the code estimates appropriate |
|
710 C for continuing. You may want to change them |
|
711 C yourself. If you are sure you want to continue |
|
712 C with relaxed error tolerances, set INFO(1)=1 and |
|
713 C call the code again. |
|
714 C |
|
715 C IDID = -3, A solution component is zero and you set the |
|
716 C corresponding component of ATOL to zero. If you |
|
717 C are sure you want to continue, you must first |
|
718 C alter the error criterion to use positive values |
|
719 C for those components of ATOL corresponding to zero |
|
720 C solution components, then set INFO(1)=1 and call |
|
721 C the code again. |
|
722 C |
|
723 C IDID = -4,-5 --- Cannot occur with this code. |
|
724 C |
|
725 C IDID = -6, Repeated error test failures occurred on the |
|
726 C last attempted step in DDASSL. A singularity in the |
|
727 C solution may be present. If you are absolutely |
|
728 C certain you want to continue, you should restart |
|
729 C the integration. (Provide initial values of Y and |
|
730 C YPRIME which are consistent) |
|
731 C |
|
732 C IDID = -7, Repeated convergence test failures occurred |
|
733 C on the last attempted step in DDASSL. An inaccurate |
|
734 C or ill-conditioned JACOBIAN may be the problem. If |
|
735 C you are absolutely certain you want to continue, you |
|
736 C should restart the integration. |
|
737 C |
|
738 C IDID = -8, The matrix of partial derivatives is singular. |
|
739 C Some of your equations may be redundant. |
|
740 C DDASSL cannot solve the problem as stated. |
|
741 C It is possible that the redundant equations |
|
742 C could be removed, and then DDASSL could |
|
743 C solve the problem. It is also possible |
|
744 C that a solution to your problem either |
|
745 C does not exist or is not unique. |
|
746 C |
|
747 C IDID = -9, DDASSL had multiple convergence test |
|
748 C failures, preceeded by multiple error |
|
749 C test failures, on the last attempted step. |
|
750 C It is possible that your problem |
|
751 C is ill-posed, and cannot be solved |
|
752 C using this code. Or, there may be a |
|
753 C discontinuity or a singularity in the |
|
754 C solution. If you are absolutely certain |
|
755 C you want to continue, you should restart |
|
756 C the integration. |
|
757 C |
|
758 C IDID =-10, DDASSL had multiple convergence test failures |
|
759 C because IRES was equal to minus one. |
|
760 C If you are absolutely certain you want |
|
761 C to continue, you should restart the |
|
762 C integration. |
|
763 C |
|
764 C IDID =-11, IRES=-2 was encountered, and control is being |
|
765 C returned to the calling program. |
|
766 C |
|
767 C IDID =-12, DDASSL failed to compute the initial YPRIME. |
|
768 C This could happen because the initial |
|
769 C approximation to YPRIME was not very good, or |
|
770 C if a YPRIME consistent with the initial Y |
|
771 C does not exist. The problem could also be caused |
|
772 C by an inaccurate or singular iteration matrix. |
|
773 C |
|
774 C IDID = -13,..,-32 --- Cannot occur with this code. |
|
775 C |
|
776 C |
|
777 C *** Following a terminated task *** |
|
778 C |
|
779 C If IDID= -33, you cannot continue the solution of this problem. |
|
780 C An attempt to do so will result in your |
|
781 C run being terminated. |
|
782 C |
|
783 C |
|
784 C -------- ERROR MESSAGES --------------------------------------------- |
|
785 C |
|
786 C The SLATEC error print routine XERMSG is called in the event of |
|
787 C unsuccessful completion of a task. Most of these are treated as |
|
788 C "recoverable errors", which means that (unless the user has directed |
|
789 C otherwise) control will be returned to the calling program for |
|
790 C possible action after the message has been printed. |
|
791 C |
|
792 C In the event of a negative value of IDID other than -33, an appro- |
|
793 C priate message is printed and the "error number" printed by XERMSG |
|
794 C is the value of IDID. There are quite a number of illegal input |
|
795 C errors that can lead to a returned value IDID=-33. The conditions |
|
796 C and their printed "error numbers" are as follows: |
|
797 C |
|
798 C Error number Condition |
|
799 C |
|
800 C 1 Some element of INFO vector is not zero or one. |
|
801 C 2 NEQ .le. 0 |
|
802 C 3 MAXORD not in range. |
|
803 C 4 LRW is less than the required length for RWORK. |
|
804 C 5 LIW is less than the required length for IWORK. |
|
805 C 6 Some element of RTOL is .lt. 0 |
|
806 C 7 Some element of ATOL is .lt. 0 |
|
807 C 8 All elements of RTOL and ATOL are zero. |
|
808 C 9 INFO(4)=1 and TSTOP is behind TOUT. |
|
809 C 10 HMAX .lt. 0.0 |
|
810 C 11 TOUT is behind T. |
|
811 C 12 INFO(8)=1 and H0=0.0 |
|
812 C 13 Some element of WT is .le. 0.0 |
|
813 C 14 TOUT is too close to T to start integration. |
|
814 C 15 INFO(4)=1 and TSTOP is behind T. |
|
815 C 16 --( Not used in this version )-- |
|
816 C 17 ML illegal. Either .lt. 0 or .gt. NEQ |
|
817 C 18 MU illegal. Either .lt. 0 or .gt. NEQ |
|
818 C 19 TOUT = T. |
|
819 C |
|
820 C If DDASSL is called again without any action taken to remove the |
|
821 C cause of an unsuccessful return, XERMSG will be called with a fatal |
|
822 C error flag, which will cause unconditional termination of the |
|
823 C program. There are two such fatal errors: |
|
824 C |
|
825 C Error number -998: The last step was terminated with a negative |
|
826 C value of IDID other than -33, and no appropriate action was |
|
827 C taken. |
|
828 C |
|
829 C Error number -999: The previous call was terminated because of |
|
830 C illegal input (IDID=-33) and there is illegal input in the |
|
831 C present call, as well. (Suspect infinite loop.) |
|
832 C |
|
833 C --------------------------------------------------------------------- |
|
834 C |
|
835 C***REFERENCES A DESCRIPTION OF DASSL: A DIFFERENTIAL/ALGEBRAIC |
|
836 C SYSTEM SOLVER, L. R. PETZOLD, SAND82-8637, |
|
837 C SANDIA NATIONAL LABORATORIES, SEPTEMBER 1982. |
|
838 C***ROUTINES CALLED D1MACH, DDAINI, DDANRM, DDASTP, DDATRP, DDAWTS, |
|
839 C XERMSG |
|
840 C***REVISION HISTORY (YYMMDD) |
|
841 C 830315 DATE WRITTEN |
|
842 C 880387 Code changes made. All common statements have been |
|
843 C replaced by a DATA statement, which defines pointers into |
|
844 C RWORK, and PARAMETER statements which define pointers |
|
845 C into IWORK. As well the documentation has gone through |
|
846 C grammatical changes. |
|
847 C 881005 The prologue has been changed to mixed case. |
|
848 C The subordinate routines had revision dates changed to |
|
849 C this date, although the documentation for these routines |
|
850 C is all upper case. No code changes. |
|
851 C 890511 Code changes made. The DATA statement in the declaration |
|
852 C section of DDASSL was replaced with a PARAMETER |
|
853 C statement. Also the statement S = 100.D0 was removed |
|
854 C from the top of the Newton iteration in DDASTP. |
|
855 C The subordinate routines had revision dates changed to |
|
856 C this date. |
|
857 C 890517 The revision date syntax was replaced with the revision |
|
858 C history syntax. Also the "DECK" comment was added to |
|
859 C the top of all subroutines. These changes are consistent |
|
860 C with new SLATEC guidelines. |
|
861 C The subordinate routines had revision dates changed to |
|
862 C this date. No code changes. |
|
863 C 891013 Code changes made. |
|
864 C Removed all occurrances of FLOAT or DBLE. All operations |
|
865 C are now performed with "mixed-mode" arithmetic. |
|
866 C Also, specific function names were replaced with generic |
|
867 C function names to be consistent with new SLATEC guidelines. |
|
868 C In particular: |
|
869 C Replaced DSQRT with SQRT everywhere. |
|
870 C Replaced DABS with ABS everywhere. |
|
871 C Replaced DMIN1 with MIN everywhere. |
|
872 C Replaced MIN0 with MIN everywhere. |
|
873 C Replaced DMAX1 with MAX everywhere. |
|
874 C Replaced MAX0 with MAX everywhere. |
|
875 C Replaced DSIGN with SIGN everywhere. |
|
876 C Also replaced REVISION DATE with REVISION HISTORY in all |
|
877 C subordinate routines. |
|
878 C 901004 Miscellaneous changes to prologue to complete conversion |
|
879 C to SLATEC 4.0 format. No code changes. (F.N.Fritsch) |
|
880 C 901009 Corrected GAMS classification code and converted subsidiary |
|
881 C routines to 4.0 format. No code changes. (F.N.Fritsch) |
|
882 C 901010 Converted XERRWV calls to XERMSG calls. (R.Clemens,AFWL) |
|
883 C 901019 Code changes made. |
|
884 C Merged SLATEC 4.0 changes with previous changes made |
|
885 C by C. Ulrich. Below is a history of the changes made by |
|
886 C C. Ulrich. (Changes in subsidiary routines are implied |
|
887 C by this history) |
|
888 C 891228 Bug was found and repaired inside the DDASSL |
|
889 C and DDAINI routines. DDAINI was incorrectly |
|
890 C returning the initial T with Y and YPRIME |
|
891 C computed at T+H. The routine now returns T+H |
|
892 C rather than the initial T. |
|
893 C Cosmetic changes made to DDASTP. |
|
894 C 900904 Three modifications were made to fix a bug (inside |
|
895 C DDASSL) re interpolation for continuation calls and |
|
896 C cases where TN is very close to TSTOP: |
|
897 C |
|
898 C 1) In testing for whether H is too large, just |
|
899 C compare H to (TSTOP - TN), rather than |
|
900 C (TSTOP - TN) * (1-4*UROUND), and set H to |
|
901 C TSTOP - TN. This will force DDASTP to step |
|
902 C exactly to TSTOP under certain situations |
|
903 C (i.e. when H returned from DDASTP would otherwise |
|
904 C take TN beyond TSTOP). |
|
905 C |
|
906 C 2) Inside the DDASTP loop, interpolate exactly to |
|
907 C TSTOP if TN is very close to TSTOP (rather than |
|
908 C interpolating to within roundoff of TSTOP). |
|
909 C |
|
910 C 3) Modified IDID description for IDID = 2 to say that |
|
911 C the solution is returned by stepping exactly to |
|
912 C TSTOP, rather than TOUT. (In some cases the |
|
913 C solution is actually obtained by extrapolating |
|
914 C over a distance near unit roundoff to TSTOP, |
|
915 C but this small distance is deemed acceptable in |
|
916 C these circumstances.) |
|
917 C 901026 Added explicit declarations for all variables and minor |
|
918 C cosmetic changes to prologue, removed unreferenced labels, |
|
919 C and improved XERMSG calls. (FNF) |
|
920 C 901030 Added ERROR MESSAGES section and reworked other sections to |
|
921 C be of more uniform format. (FNF) |
|
922 C 910624 Fixed minor bug related to HMAX (five lines ending in |
|
923 C statement 526 in DDASSL). (LRP) |
|
924 C |
|
925 C***END PROLOGUE DDASSL |
|
926 C |
|
927 C**End |
|
928 C |
|
929 C Declare arguments. |
|
930 C |
|
931 INTEGER NEQ, INFO(15), IDID, LRW, IWORK(*), LIW, IPAR(*) |
|
932 DOUBLE PRECISION |
|
933 * T, Y(*), YPRIME(*), TOUT, RTOL(*), ATOL(*), RWORK(*), |
|
934 * RPAR(*) |
|
935 EXTERNAL RES, JAC |
|
936 C |
|
937 C Declare externals. |
|
938 C |
|
939 EXTERNAL D1MACH, DDAINI, DDANRM, DDASTP, DDATRP, DDAWTS, XERMSG |
|
940 DOUBLE PRECISION D1MACH, DDANRM |
|
941 C |
|
942 C Declare local variables. |
|
943 C |
|
944 INTEGER I, ITEMP, LALPHA, LBETA, LCJ, LCJOLD, LCTF, LDELTA, |
|
945 * LENIW, LENPD, LENRW, LE, LETF, LGAMMA, LH, LHMAX, LHOLD, LIPVT, |
|
946 * LJCALC, LK, LKOLD, LIWM, LML, LMTYPE, LMU, LMXORD, LNJE, LNPD, |
|
947 * LNRE, LNS, LNST, LNSTL, LPD, LPHASE, LPHI, LPSI, LROUND, LS, |
|
948 * LSIGMA, LTN, LTSTOP, LWM, LWT, MBAND, MSAVE, MXORD, NPD, NTEMP, |
|
949 * NZFLG |
|
950 DOUBLE PRECISION |
|
951 * ATOLI, H, HMAX, HMIN, HO, R, RH, RTOLI, TDIST, TN, TNEXT, |
|
952 * TSTOP, UROUND, YPNORM |
|
953 LOGICAL DONE |
|
954 C Auxiliary variables for conversion of values to be included in |
|
955 C error messages. |
|
956 CHARACTER*8 XERN1, XERN2 |
|
957 CHARACTER*16 XERN3, XERN4 |
|
958 C |
|
959 C SET POINTERS INTO IWORK |
|
960 PARAMETER (LML=1, LMU=2, LMXORD=3, LMTYPE=4, LNST=11, |
|
961 * LNRE=12, LNJE=13, LETF=14, LCTF=15, LNPD=16, |
|
962 * LIPVT=21, LJCALC=5, LPHASE=6, LK=7, LKOLD=8, |
|
963 * LNS=9, LNSTL=10, LIWM=1) |
|
964 C |
|
965 C SET RELATIVE OFFSET INTO RWORK |
|
966 PARAMETER (NPD=1) |
|
967 C |
|
968 C SET POINTERS INTO RWORK |
|
969 PARAMETER (LTSTOP=1, LHMAX=2, LH=3, LTN=4, |
|
970 * LCJ=5, LCJOLD=6, LHOLD=7, LS=8, LROUND=9, |
|
971 * LALPHA=11, LBETA=17, LGAMMA=23, |
|
972 * LPSI=29, LSIGMA=35, LDELTA=41) |
|
973 C |
|
974 C***FIRST EXECUTABLE STATEMENT DDASSL |
|
975 IF(INFO(1).NE.0)GO TO 100 |
|
976 C |
|
977 C----------------------------------------------------------------------- |
|
978 C THIS BLOCK IS EXECUTED FOR THE INITIAL CALL ONLY. |
|
979 C IT CONTAINS CHECKING OF INPUTS AND INITIALIZATIONS. |
|
980 C----------------------------------------------------------------------- |
|
981 C |
|
982 C FIRST CHECK INFO ARRAY TO MAKE SURE ALL ELEMENTS OF INFO |
|
983 C ARE EITHER ZERO OR ONE. |
|
984 DO 10 I=2,11 |
|
985 IF(INFO(I).NE.0.AND.INFO(I).NE.1)GO TO 701 |
|
986 10 CONTINUE |
|
987 C |
|
988 IF(NEQ.LE.0)GO TO 702 |
|
989 C |
|
990 C CHECK AND COMPUTE MAXIMUM ORDER |
|
991 MXORD=5 |
|
992 IF(INFO(9).EQ.0)GO TO 20 |
|
993 MXORD=IWORK(LMXORD) |
|
994 IF(MXORD.LT.1.OR.MXORD.GT.5)GO TO 703 |
|
995 20 IWORK(LMXORD)=MXORD |
|
996 C |
|
997 C COMPUTE MTYPE,LENPD,LENRW.CHECK ML AND MU. |
|
998 IF(INFO(6).NE.0)GO TO 40 |
|
999 LENPD=NEQ**2 |
|
1000 LENRW=40+(IWORK(LMXORD)+4)*NEQ+LENPD |
|
1001 IF(INFO(5).NE.0)GO TO 30 |
|
1002 IWORK(LMTYPE)=2 |
|
1003 GO TO 60 |
|
1004 30 IWORK(LMTYPE)=1 |
|
1005 GO TO 60 |
|
1006 40 IF(IWORK(LML).LT.0.OR.IWORK(LML).GE.NEQ)GO TO 717 |
|
1007 IF(IWORK(LMU).LT.0.OR.IWORK(LMU).GE.NEQ)GO TO 718 |
|
1008 LENPD=(2*IWORK(LML)+IWORK(LMU)+1)*NEQ |
|
1009 IF(INFO(5).NE.0)GO TO 50 |
|
1010 IWORK(LMTYPE)=5 |
|
1011 MBAND=IWORK(LML)+IWORK(LMU)+1 |
|
1012 MSAVE=(NEQ/MBAND)+1 |
|
1013 LENRW=40+(IWORK(LMXORD)+4)*NEQ+LENPD+2*MSAVE |
|
1014 GO TO 60 |
|
1015 50 IWORK(LMTYPE)=4 |
|
1016 LENRW=40+(IWORK(LMXORD)+4)*NEQ+LENPD |
|
1017 C |
|
1018 C CHECK LENGTHS OF RWORK AND IWORK |
|
1019 60 LENIW=20+NEQ |
|
1020 IWORK(LNPD)=LENPD |
|
1021 IF(LRW.LT.LENRW)GO TO 704 |
|
1022 IF(LIW.LT.LENIW)GO TO 705 |
|
1023 C |
|
1024 C CHECK TO SEE THAT TOUT IS DIFFERENT FROM T |
|
1025 IF(TOUT .EQ. T)GO TO 719 |
|
1026 C |
|
1027 C CHECK HMAX |
|
1028 IF(INFO(7).EQ.0)GO TO 70 |
|
1029 HMAX=RWORK(LHMAX) |
|
1030 IF(HMAX.LE.0.0D0)GO TO 710 |
|
1031 70 CONTINUE |
|
1032 C |
|
1033 C INITIALIZE COUNTERS |
|
1034 IWORK(LNST)=0 |
|
1035 IWORK(LNRE)=0 |
|
1036 IWORK(LNJE)=0 |
|
1037 C |
|
1038 IWORK(LNSTL)=0 |
|
1039 IDID=1 |
|
1040 GO TO 200 |
|
1041 C |
|
1042 C----------------------------------------------------------------------- |
|
1043 C THIS BLOCK IS FOR CONTINUATION CALLS |
|
1044 C ONLY. HERE WE CHECK INFO(1),AND IF THE |
|
1045 C LAST STEP WAS INTERRUPTED WE CHECK WHETHER |
|
1046 C APPROPRIATE ACTION WAS TAKEN. |
|
1047 C----------------------------------------------------------------------- |
|
1048 C |
|
1049 100 CONTINUE |
|
1050 IF(INFO(1).EQ.1)GO TO 110 |
|
1051 IF(INFO(1).NE.-1)GO TO 701 |
|
1052 C |
|
1053 C IF WE ARE HERE, THE LAST STEP WAS INTERRUPTED |
|
1054 C BY AN ERROR CONDITION FROM DDASTP,AND |
|
1055 C APPROPRIATE ACTION WAS NOT TAKEN. THIS |
|
1056 C IS A FATAL ERROR. |
|
1057 WRITE (XERN1, '(I8)') IDID |
|
1058 CALL XERMSG ('SLATEC', 'DDASSL', |
|
1059 * 'THE LAST STEP TERMINATED WITH A NEGATIVE VALUE OF IDID = ' // |
|
1060 * XERN1 // ' AND NO APPROPRIATE ACTION WAS TAKEN. ' // |
|
1061 * 'RUN TERMINATED', -998, 2) |
|
1062 RETURN |
|
1063 110 CONTINUE |
|
1064 IWORK(LNSTL)=IWORK(LNST) |
|
1065 C |
|
1066 C----------------------------------------------------------------------- |
|
1067 C THIS BLOCK IS EXECUTED ON ALL CALLS. |
|
1068 C THE ERROR TOLERANCE PARAMETERS ARE |
|
1069 C CHECKED, AND THE WORK ARRAY POINTERS |
|
1070 C ARE SET. |
|
1071 C----------------------------------------------------------------------- |
|
1072 C |
|
1073 200 CONTINUE |
|
1074 C CHECK RTOL,ATOL |
|
1075 NZFLG=0 |
|
1076 RTOLI=RTOL(1) |
|
1077 ATOLI=ATOL(1) |
|
1078 DO 210 I=1,NEQ |
|
1079 IF(INFO(2).EQ.1)RTOLI=RTOL(I) |
|
1080 IF(INFO(2).EQ.1)ATOLI=ATOL(I) |
|
1081 IF(RTOLI.GT.0.0D0.OR.ATOLI.GT.0.0D0)NZFLG=1 |
|
1082 IF(RTOLI.LT.0.0D0)GO TO 706 |
|
1083 IF(ATOLI.LT.0.0D0)GO TO 707 |
|
1084 210 CONTINUE |
|
1085 IF(NZFLG.EQ.0)GO TO 708 |
|
1086 C |
|
1087 C SET UP RWORK STORAGE.IWORK STORAGE IS FIXED |
|
1088 C IN DATA STATEMENT. |
|
1089 LE=LDELTA+NEQ |
|
1090 LWT=LE+NEQ |
|
1091 LPHI=LWT+NEQ |
|
1092 LPD=LPHI+(IWORK(LMXORD)+1)*NEQ |
|
1093 LWM=LPD |
|
1094 NTEMP=NPD+IWORK(LNPD) |
|
1095 IF(INFO(1).EQ.1)GO TO 400 |
|
1096 C |
|
1097 C----------------------------------------------------------------------- |
|
1098 C THIS BLOCK IS EXECUTED ON THE INITIAL CALL |
|
1099 C ONLY. SET THE INITIAL STEP SIZE, AND |
|
1100 C THE ERROR WEIGHT VECTOR, AND PHI. |
|
1101 C COMPUTE INITIAL YPRIME, IF NECESSARY. |
|
1102 C----------------------------------------------------------------------- |
|
1103 C |
|
1104 TN=T |
|
1105 IDID=1 |
|
1106 C |
|
1107 C SET ERROR WEIGHT VECTOR WT |
|
1108 CALL DDAWTS(NEQ,INFO(2),RTOL,ATOL,Y,RWORK(LWT),RPAR,IPAR) |
|
1109 DO 305 I = 1,NEQ |
|
1110 IF(RWORK(LWT+I-1).LE.0.0D0) GO TO 713 |
|
1111 305 CONTINUE |
|
1112 C |
|
1113 C COMPUTE UNIT ROUNDOFF AND HMIN |
|
1114 UROUND = D1MACH(4) |
|
1115 RWORK(LROUND) = UROUND |
|
1116 HMIN = 4.0D0*UROUND*MAX(ABS(T),ABS(TOUT)) |
|
1117 C |
|
1118 C CHECK INITIAL INTERVAL TO SEE THAT IT IS LONG ENOUGH |
|
1119 TDIST = ABS(TOUT - T) |
|
1120 IF(TDIST .LT. HMIN) GO TO 714 |
|
1121 C |
|
1122 C CHECK HO, IF THIS WAS INPUT |
|
1123 IF (INFO(8) .EQ. 0) GO TO 310 |
|
1124 HO = RWORK(LH) |
|
1125 IF ((TOUT - T)*HO .LT. 0.0D0) GO TO 711 |
|
1126 IF (HO .EQ. 0.0D0) GO TO 712 |
|
1127 GO TO 320 |
|
1128 310 CONTINUE |
|
1129 C |
|
1130 C COMPUTE INITIAL STEPSIZE, TO BE USED BY EITHER |
|
1131 C DDASTP OR DDAINI, DEPENDING ON INFO(11) |
|
1132 HO = 0.001D0*TDIST |
|
1133 YPNORM = DDANRM(NEQ,YPRIME,RWORK(LWT),RPAR,IPAR) |
|
1134 IF (YPNORM .GT. 0.5D0/HO) HO = 0.5D0/YPNORM |
|
1135 HO = SIGN(HO,TOUT-T) |
|
1136 C ADJUST HO IF NECESSARY TO MEET HMAX BOUND |
|
1137 320 IF (INFO(7) .EQ. 0) GO TO 330 |
|
1138 RH = ABS(HO)/RWORK(LHMAX) |
|
1139 IF (RH .GT. 1.0D0) HO = HO/RH |
|
1140 C COMPUTE TSTOP, IF APPLICABLE |
|
1141 330 IF (INFO(4) .EQ. 0) GO TO 340 |
|
1142 TSTOP = RWORK(LTSTOP) |
|
1143 IF ((TSTOP - T)*HO .LT. 0.0D0) GO TO 715 |
|
1144 IF ((T + HO - TSTOP)*HO .GT. 0.0D0) HO = TSTOP - T |
|
1145 IF ((TSTOP - TOUT)*HO .LT. 0.0D0) GO TO 709 |
|
1146 C |
|
1147 C COMPUTE INITIAL DERIVATIVE, UPDATING TN AND Y, IF APPLICABLE |
|
1148 340 IF (INFO(11) .EQ. 0) GO TO 350 |
|
1149 CALL DDAINI(TN,Y,YPRIME,NEQ, |
|
1150 * RES,JAC,HO,RWORK(LWT),IDID,RPAR,IPAR, |
|
1151 * RWORK(LPHI),RWORK(LDELTA),RWORK(LE), |
|
1152 * RWORK(LWM),IWORK(LIWM),HMIN,RWORK(LROUND), |
|
1153 * INFO(10),NTEMP) |
|
1154 IF (IDID .LT. 0) GO TO 390 |
|
1155 C |
|
1156 C LOAD H WITH HO. STORE H IN RWORK(LH) |
|
1157 350 H = HO |
|
1158 RWORK(LH) = H |
|
1159 C |
|
1160 C LOAD Y AND H*YPRIME INTO PHI(*,1) AND PHI(*,2) |
|
1161 ITEMP = LPHI + NEQ |
|
1162 DO 370 I = 1,NEQ |
|
1163 RWORK(LPHI + I - 1) = Y(I) |
|
1164 370 RWORK(ITEMP + I - 1) = H*YPRIME(I) |
|
1165 C |
|
1166 390 GO TO 500 |
|
1167 C |
|
1168 C------------------------------------------------------- |
|
1169 C THIS BLOCK IS FOR CONTINUATION CALLS ONLY. ITS |
|
1170 C PURPOSE IS TO CHECK STOP CONDITIONS BEFORE |
|
1171 C TAKING A STEP. |
|
1172 C ADJUST H IF NECESSARY TO MEET HMAX BOUND |
|
1173 C------------------------------------------------------- |
|
1174 C |
|
1175 400 CONTINUE |
|
1176 UROUND=RWORK(LROUND) |
|
1177 DONE = .FALSE. |
|
1178 TN=RWORK(LTN) |
|
1179 H=RWORK(LH) |
|
1180 IF(INFO(7) .EQ. 0) GO TO 410 |
|
1181 RH = ABS(H)/RWORK(LHMAX) |
|
1182 IF(RH .GT. 1.0D0) H = H/RH |
|
1183 410 CONTINUE |
|
1184 IF(T .EQ. TOUT) GO TO 719 |
|
1185 IF((T - TOUT)*H .GT. 0.0D0) GO TO 711 |
|
1186 IF(INFO(4) .EQ. 1) GO TO 430 |
|
1187 IF(INFO(3) .EQ. 1) GO TO 420 |
|
1188 IF((TN-TOUT)*H.LT.0.0D0)GO TO 490 |
|
1189 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD), |
|
1190 * RWORK(LPHI),RWORK(LPSI)) |
|
1191 T=TOUT |
|
1192 IDID = 3 |
|
1193 DONE = .TRUE. |
|
1194 GO TO 490 |
|
1195 420 IF((TN-T)*H .LE. 0.0D0) GO TO 490 |
|
1196 IF((TN - TOUT)*H .GT. 0.0D0) GO TO 425 |
|
1197 CALL DDATRP(TN,TN,Y,YPRIME,NEQ,IWORK(LKOLD), |
|
1198 * RWORK(LPHI),RWORK(LPSI)) |
|
1199 T = TN |
|
1200 IDID = 1 |
|
1201 DONE = .TRUE. |
|
1202 GO TO 490 |
|
1203 425 CONTINUE |
|
1204 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD), |
|
1205 * RWORK(LPHI),RWORK(LPSI)) |
|
1206 T = TOUT |
|
1207 IDID = 3 |
|
1208 DONE = .TRUE. |
|
1209 GO TO 490 |
|
1210 430 IF(INFO(3) .EQ. 1) GO TO 440 |
|
1211 TSTOP=RWORK(LTSTOP) |
|
1212 IF((TN-TSTOP)*H.GT.0.0D0) GO TO 715 |
|
1213 IF((TSTOP-TOUT)*H.LT.0.0D0)GO TO 709 |
|
1214 IF((TN-TOUT)*H.LT.0.0D0)GO TO 450 |
|
1215 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD), |
|
1216 * RWORK(LPHI),RWORK(LPSI)) |
|
1217 T=TOUT |
|
1218 IDID = 3 |
|
1219 DONE = .TRUE. |
|
1220 GO TO 490 |
|
1221 440 TSTOP = RWORK(LTSTOP) |
|
1222 IF((TN-TSTOP)*H .GT. 0.0D0) GO TO 715 |
|
1223 IF((TSTOP-TOUT)*H .LT. 0.0D0) GO TO 709 |
|
1224 IF((TN-T)*H .LE. 0.0D0) GO TO 450 |
|
1225 IF((TN - TOUT)*H .GT. 0.0D0) GO TO 445 |
|
1226 CALL DDATRP(TN,TN,Y,YPRIME,NEQ,IWORK(LKOLD), |
|
1227 * RWORK(LPHI),RWORK(LPSI)) |
|
1228 T = TN |
|
1229 IDID = 1 |
|
1230 DONE = .TRUE. |
|
1231 GO TO 490 |
|
1232 445 CONTINUE |
|
1233 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD), |
|
1234 * RWORK(LPHI),RWORK(LPSI)) |
|
1235 T = TOUT |
|
1236 IDID = 3 |
|
1237 DONE = .TRUE. |
|
1238 GO TO 490 |
|
1239 450 CONTINUE |
|
1240 C CHECK WHETHER WE ARE WITHIN ROUNDOFF OF TSTOP |
|
1241 IF(ABS(TN-TSTOP).GT.100.0D0*UROUND* |
|
1242 * (ABS(TN)+ABS(H)))GO TO 460 |
|
1243 CALL DDATRP(TN,TSTOP,Y,YPRIME,NEQ,IWORK(LKOLD), |
|
1244 * RWORK(LPHI),RWORK(LPSI)) |
|
1245 IDID=2 |
|
1246 T=TSTOP |
|
1247 DONE = .TRUE. |
|
1248 GO TO 490 |
|
1249 460 TNEXT=TN+H |
|
1250 IF((TNEXT-TSTOP)*H.LE.0.0D0)GO TO 490 |
|
1251 H=TSTOP-TN |
|
1252 RWORK(LH)=H |
|
1253 C |
|
1254 490 IF (DONE) GO TO 580 |
|
1255 C |
|
1256 C------------------------------------------------------- |
|
1257 C THE NEXT BLOCK CONTAINS THE CALL TO THE |
|
1258 C ONE-STEP INTEGRATOR DDASTP. |
|
1259 C THIS IS A LOOPING POINT FOR THE INTEGRATION STEPS. |
|
1260 C CHECK FOR TOO MANY STEPS. |
|
1261 C UPDATE WT. |
|
1262 C CHECK FOR TOO MUCH ACCURACY REQUESTED. |
|
1263 C COMPUTE MINIMUM STEPSIZE. |
|
1264 C------------------------------------------------------- |
|
1265 C |
|
1266 500 CONTINUE |
|
1267 C CHECK FOR FAILURE TO COMPUTE INITIAL YPRIME |
|
1268 IF (IDID .EQ. -12) GO TO 527 |
|
1269 C |
|
1270 C CHECK FOR TOO MANY STEPS |
|
1271 IF((IWORK(LNST)-IWORK(LNSTL)).LT.500) |
|
1272 * GO TO 510 |
|
1273 IDID=-1 |
|
1274 GO TO 527 |
|
1275 C |
|
1276 C UPDATE WT |
|
1277 510 CALL DDAWTS(NEQ,INFO(2),RTOL,ATOL,RWORK(LPHI), |
|
1278 * RWORK(LWT),RPAR,IPAR) |
|
1279 DO 520 I=1,NEQ |
|
1280 IF(RWORK(I+LWT-1).GT.0.0D0)GO TO 520 |
|
1281 IDID=-3 |
|
1282 GO TO 527 |
|
1283 520 CONTINUE |
|
1284 C |
|
1285 C TEST FOR TOO MUCH ACCURACY REQUESTED. |
|
1286 R=DDANRM(NEQ,RWORK(LPHI),RWORK(LWT),RPAR,IPAR)* |
|
1287 * 100.0D0*UROUND |
|
1288 IF(R.LE.1.0D0)GO TO 525 |
|
1289 C MULTIPLY RTOL AND ATOL BY R AND RETURN |
|
1290 IF(INFO(2).EQ.1)GO TO 523 |
|
1291 RTOL(1)=R*RTOL(1) |
|
1292 ATOL(1)=R*ATOL(1) |
|
1293 IDID=-2 |
|
1294 GO TO 527 |
|
1295 523 DO 524 I=1,NEQ |
|
1296 RTOL(I)=R*RTOL(I) |
|
1297 524 ATOL(I)=R*ATOL(I) |
|
1298 IDID=-2 |
|
1299 GO TO 527 |
|
1300 525 CONTINUE |
|
1301 C |
|
1302 C COMPUTE MINIMUM STEPSIZE |
|
1303 HMIN=4.0D0*UROUND*MAX(ABS(TN),ABS(TOUT)) |
|
1304 C |
|
1305 C TEST H VS. HMAX |
|
1306 IF (INFO(7) .EQ. 0) GO TO 526 |
|
1307 RH = ABS(H)/RWORK(LHMAX) |
|
1308 IF (RH .GT. 1.0D0) H = H/RH |
|
1309 526 CONTINUE |
|
1310 C |
|
1311 CALL DDASTP(TN,Y,YPRIME,NEQ, |
|
1312 * RES,JAC,H,RWORK(LWT),INFO(1),IDID,RPAR,IPAR, |
|
1313 * RWORK(LPHI),RWORK(LDELTA),RWORK(LE), |
|
1314 * RWORK(LWM),IWORK(LIWM), |
|
1315 * RWORK(LALPHA),RWORK(LBETA),RWORK(LGAMMA), |
|
1316 * RWORK(LPSI),RWORK(LSIGMA), |
|
1317 * RWORK(LCJ),RWORK(LCJOLD),RWORK(LHOLD), |
|
1318 * RWORK(LS),HMIN,RWORK(LROUND), |
|
1319 * IWORK(LPHASE),IWORK(LJCALC),IWORK(LK), |
|
1320 * IWORK(LKOLD),IWORK(LNS),INFO(10),NTEMP) |
|
1321 527 IF(IDID.LT.0)GO TO 600 |
|
1322 C |
|
1323 C-------------------------------------------------------- |
|
1324 C THIS BLOCK HANDLES THE CASE OF A SUCCESSFUL RETURN |
|
1325 C FROM DDASTP (IDID=1). TEST FOR STOP CONDITIONS. |
|
1326 C-------------------------------------------------------- |
|
1327 C |
|
1328 IF(INFO(4).NE.0)GO TO 540 |
|
1329 IF(INFO(3).NE.0)GO TO 530 |
|
1330 IF((TN-TOUT)*H.LT.0.0D0)GO TO 500 |
|
1331 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ, |
|
1332 * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI)) |
|
1333 IDID=3 |
|
1334 T=TOUT |
|
1335 GO TO 580 |
|
1336 530 IF((TN-TOUT)*H.GE.0.0D0)GO TO 535 |
|
1337 T=TN |
|
1338 IDID=1 |
|
1339 GO TO 580 |
|
1340 535 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ, |
|
1341 * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI)) |
|
1342 IDID=3 |
|
1343 T=TOUT |
|
1344 GO TO 580 |
|
1345 540 IF(INFO(3).NE.0)GO TO 550 |
|
1346 IF((TN-TOUT)*H.LT.0.0D0)GO TO 542 |
|
1347 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ, |
|
1348 * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI)) |
|
1349 T=TOUT |
|
1350 IDID=3 |
|
1351 GO TO 580 |
|
1352 542 IF(ABS(TN-TSTOP).LE.100.0D0*UROUND* |
|
1353 * (ABS(TN)+ABS(H)))GO TO 545 |
|
1354 TNEXT=TN+H |
|
1355 IF((TNEXT-TSTOP)*H.LE.0.0D0)GO TO 500 |
|
1356 H=TSTOP-TN |
|
1357 GO TO 500 |
|
1358 545 CALL DDATRP(TN,TSTOP,Y,YPRIME,NEQ, |
|
1359 * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI)) |
|
1360 IDID=2 |
|
1361 T=TSTOP |
|
1362 GO TO 580 |
|
1363 550 IF((TN-TOUT)*H.GE.0.0D0)GO TO 555 |
|
1364 IF(ABS(TN-TSTOP).LE.100.0D0*UROUND*(ABS(TN)+ABS(H)))GO TO 552 |
|
1365 T=TN |
|
1366 IDID=1 |
|
1367 GO TO 580 |
|
1368 552 CALL DDATRP(TN,TSTOP,Y,YPRIME,NEQ, |
|
1369 * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI)) |
|
1370 IDID=2 |
|
1371 T=TSTOP |
|
1372 GO TO 580 |
|
1373 555 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ, |
|
1374 * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI)) |
|
1375 T=TOUT |
|
1376 IDID=3 |
|
1377 GO TO 580 |
|
1378 C |
|
1379 C-------------------------------------------------------- |
|
1380 C ALL SUCCESSFUL RETURNS FROM DDASSL ARE MADE FROM |
|
1381 C THIS BLOCK. |
|
1382 C-------------------------------------------------------- |
|
1383 C |
|
1384 580 CONTINUE |
|
1385 RWORK(LTN)=TN |
|
1386 RWORK(LH)=H |
|
1387 RETURN |
|
1388 C |
|
1389 C----------------------------------------------------------------------- |
|
1390 C THIS BLOCK HANDLES ALL UNSUCCESSFUL |
|
1391 C RETURNS OTHER THAN FOR ILLEGAL INPUT. |
|
1392 C----------------------------------------------------------------------- |
|
1393 C |
|
1394 600 CONTINUE |
|
1395 ITEMP=-IDID |
|
1396 GO TO (610,620,630,690,690,640,650,660,670,675, |
|
1397 * 680,685), ITEMP |
|
1398 C |
|
1399 C THE MAXIMUM NUMBER OF STEPS WAS TAKEN BEFORE |
|
1400 C REACHING TOUT |
|
1401 610 WRITE (XERN3, '(1P,D15.6)') TN |
|
1402 CALL XERMSG ('SLATEC', 'DDASSL', |
|
1403 * 'AT CURRENT T = ' // XERN3 // ' 500 STEPS TAKEN ON THIS ' // |
|
1404 * 'CALL BEFORE REACHING TOUT', IDID, 1) |
|
1405 GO TO 690 |
|
1406 C |
|
1407 C TOO MUCH ACCURACY FOR MACHINE PRECISION |
|
1408 620 WRITE (XERN3, '(1P,D15.6)') TN |
|
1409 CALL XERMSG ('SLATEC', 'DDASSL', |
|
1410 * 'AT T = ' // XERN3 // ' TOO MUCH ACCURACY REQUESTED FOR ' // |
|
1411 * 'PRECISION OF MACHINE. RTOL AND ATOL WERE INCREASED TO ' // |
|
1412 * 'APPROPRIATE VALUES', IDID, 1) |
|
1413 GO TO 690 |
|
1414 C |
|
1415 C WT(I) .LE. 0.0 FOR SOME I (NOT AT START OF PROBLEM) |
|
1416 630 WRITE (XERN3, '(1P,D15.6)') TN |
|
1417 CALL XERMSG ('SLATEC', 'DDASSL', |
|
1418 * 'AT T = ' // XERN3 // ' SOME ELEMENT OF WT HAS BECOME .LE. ' // |
|
1419 * '0.0', IDID, 1) |
|
1420 GO TO 690 |
|
1421 C |
|
1422 C ERROR TEST FAILED REPEATEDLY OR WITH H=HMIN |
|
1423 640 WRITE (XERN3, '(1P,D15.6)') TN |
|
1424 WRITE (XERN4, '(1P,D15.6)') H |
|
1425 CALL XERMSG ('SLATEC', 'DDASSL', |
|
1426 * 'AT T = ' // XERN3 // ' AND STEPSIZE H = ' // XERN4 // |
|
1427 * ' THE ERROR TEST FAILED REPEATEDLY OR WITH ABS(H)=HMIN', |
|
1428 * IDID, 1) |
|
1429 GO TO 690 |
|
1430 C |
|
1431 C CORRECTOR CONVERGENCE FAILED REPEATEDLY OR WITH H=HMIN |
|
1432 650 WRITE (XERN3, '(1P,D15.6)') TN |
|
1433 WRITE (XERN4, '(1P,D15.6)') H |
|
1434 CALL XERMSG ('SLATEC', 'DDASSL', |
|
1435 * 'AT T = ' // XERN3 // ' AND STEPSIZE H = ' // XERN4 // |
|
1436 * ' THE CORRECTOR FAILED TO CONVERGE REPEATEDLY OR WITH ' // |
|
1437 * 'ABS(H)=HMIN', IDID, 1) |
|
1438 GO TO 690 |
|
1439 C |
|
1440 C THE ITERATION MATRIX IS SINGULAR |
|
1441 660 WRITE (XERN3, '(1P,D15.6)') TN |
|
1442 WRITE (XERN4, '(1P,D15.6)') H |
|
1443 CALL XERMSG ('SLATEC', 'DDASSL', |
|
1444 * 'AT T = ' // XERN3 // ' AND STEPSIZE H = ' // XERN4 // |
|
1445 * ' THE ITERATION MATRIX IS SINGULAR', IDID, 1) |
|
1446 GO TO 690 |
|
1447 C |
|
1448 C CORRECTOR FAILURE PRECEEDED BY ERROR TEST FAILURES. |
|
1449 670 WRITE (XERN3, '(1P,D15.6)') TN |
|
1450 WRITE (XERN4, '(1P,D15.6)') H |
|
1451 CALL XERMSG ('SLATEC', 'DDASSL', |
|
1452 * 'AT T = ' // XERN3 // ' AND STEPSIZE H = ' // XERN4 // |
|
1453 * ' THE CORRECTOR COULD NOT CONVERGE. ALSO, THE ERROR TEST ' // |
|
1454 * 'FAILED REPEATEDLY.', IDID, 1) |
|
1455 GO TO 690 |
|
1456 C |
|
1457 C CORRECTOR FAILURE BECAUSE IRES = -1 |
|
1458 675 WRITE (XERN3, '(1P,D15.6)') TN |
|
1459 WRITE (XERN4, '(1P,D15.6)') H |
|
1460 CALL XERMSG ('SLATEC', 'DDASSL', |
|
1461 * 'AT T = ' // XERN3 // ' AND STEPSIZE H = ' // XERN4 // |
|
1462 * ' THE CORRECTOR COULD NOT CONVERGE BECAUSE IRES WAS EQUAL ' // |
|
1463 * 'TO MINUS ONE', IDID, 1) |
|
1464 GO TO 690 |
|
1465 C |
|
1466 C FAILURE BECAUSE IRES = -2 |
|
1467 680 WRITE (XERN3, '(1P,D15.6)') TN |
|
1468 WRITE (XERN4, '(1P,D15.6)') H |
|
1469 CALL XERMSG ('SLATEC', 'DDASSL', |
|
1470 * 'AT T = ' // XERN3 // ' AND STEPSIZE H = ' // XERN4 // |
|
1471 * ' IRES WAS EQUAL TO MINUS TWO', IDID, 1) |
|
1472 GO TO 690 |
|
1473 C |
|
1474 C FAILED TO COMPUTE INITIAL YPRIME |
|
1475 685 WRITE (XERN3, '(1P,D15.6)') TN |
|
1476 WRITE (XERN4, '(1P,D15.6)') HO |
|
1477 CALL XERMSG ('SLATEC', 'DDASSL', |
|
1478 * 'AT T = ' // XERN3 // ' AND STEPSIZE H = ' // XERN4 // |
|
1479 * ' THE INITIAL YPRIME COULD NOT BE COMPUTED', IDID, 1) |
|
1480 GO TO 690 |
|
1481 C |
|
1482 690 CONTINUE |
|
1483 INFO(1)=-1 |
|
1484 T=TN |
|
1485 RWORK(LTN)=TN |
|
1486 RWORK(LH)=H |
|
1487 RETURN |
|
1488 C |
|
1489 C----------------------------------------------------------------------- |
|
1490 C THIS BLOCK HANDLES ALL ERROR RETURNS DUE |
|
1491 C TO ILLEGAL INPUT, AS DETECTED BEFORE CALLING |
|
1492 C DDASTP. FIRST THE ERROR MESSAGE ROUTINE IS |
|
1493 C CALLED. IF THIS HAPPENS TWICE IN |
|
1494 C SUCCESSION, EXECUTION IS TERMINATED |
|
1495 C |
|
1496 C----------------------------------------------------------------------- |
|
1497 701 CALL XERMSG ('SLATEC', 'DDASSL', |
|
1498 * 'SOME ELEMENT OF INFO VECTOR IS NOT ZERO OR ONE', 1, 1) |
|
1499 GO TO 750 |
|
1500 C |
|
1501 702 WRITE (XERN1, '(I8)') NEQ |
|
1502 CALL XERMSG ('SLATEC', 'DDASSL', |
|
1503 * 'NEQ = ' // XERN1 // ' .LE. 0', 2, 1) |
|
1504 GO TO 750 |
|
1505 C |
|
1506 703 WRITE (XERN1, '(I8)') MXORD |
|
1507 CALL XERMSG ('SLATEC', 'DDASSL', |
|
1508 * 'MAXORD = ' // XERN1 // ' NOT IN RANGE', 3, 1) |
|
1509 GO TO 750 |
|
1510 C |
|
1511 704 WRITE (XERN1, '(I8)') LENRW |
|
1512 WRITE (XERN2, '(I8)') LRW |
|
1513 CALL XERMSG ('SLATEC', 'DDASSL', |
|
1514 * 'RWORK LENGTH NEEDED, LENRW = ' // XERN1 // |
|
1515 * ', EXCEEDS LRW = ' // XERN2, 4, 1) |
|
1516 GO TO 750 |
|
1517 C |
|
1518 705 WRITE (XERN1, '(I8)') LENIW |
|
1519 WRITE (XERN2, '(I8)') LIW |
|
1520 CALL XERMSG ('SLATEC', 'DDASSL', |
|
1521 * 'IWORK LENGTH NEEDED, LENIW = ' // XERN1 // |
|
1522 * ', EXCEEDS LIW = ' // XERN2, 5, 1) |
|
1523 GO TO 750 |
|
1524 C |
|
1525 706 CALL XERMSG ('SLATEC', 'DDASSL', |
|
1526 * 'SOME ELEMENT OF RTOL IS .LT. 0', 6, 1) |
|
1527 GO TO 750 |
|
1528 C |
|
1529 707 CALL XERMSG ('SLATEC', 'DDASSL', |
|
1530 * 'SOME ELEMENT OF ATOL IS .LT. 0', 7, 1) |
|
1531 GO TO 750 |
|
1532 C |
|
1533 708 CALL XERMSG ('SLATEC', 'DDASSL', |
|
1534 * 'ALL ELEMENTS OF RTOL AND ATOL ARE ZERO', 8, 1) |
|
1535 GO TO 750 |
|
1536 C |
|
1537 709 WRITE (XERN3, '(1P,D15.6)') TSTOP |
|
1538 WRITE (XERN4, '(1P,D15.6)') TOUT |
|
1539 CALL XERMSG ('SLATEC', 'DDASSL', |
|
1540 * 'INFO(4) = 1 AND TSTOP = ' // XERN3 // ' BEHIND TOUT = ' // |
|
1541 * XERN4, 9, 1) |
|
1542 GO TO 750 |
|
1543 C |
|
1544 710 WRITE (XERN3, '(1P,D15.6)') HMAX |
|
1545 CALL XERMSG ('SLATEC', 'DDASSL', |
|
1546 * 'HMAX = ' // XERN3 // ' .LT. 0.0', 10, 1) |
|
1547 GO TO 750 |
|
1548 C |
|
1549 711 WRITE (XERN3, '(1P,D15.6)') TOUT |
|
1550 WRITE (XERN4, '(1P,D15.6)') T |
|
1551 CALL XERMSG ('SLATEC', 'DDASSL', |
|
1552 * 'TOUT = ' // XERN3 // ' BEHIND T = ' // XERN4, 11, 1) |
|
1553 GO TO 750 |
|
1554 C |
|
1555 712 CALL XERMSG ('SLATEC', 'DDASSL', |
|
1556 * 'INFO(8)=1 AND H0=0.0', 12, 1) |
|
1557 GO TO 750 |
|
1558 C |
|
1559 713 CALL XERMSG ('SLATEC', 'DDASSL', |
|
1560 * 'SOME ELEMENT OF WT IS .LE. 0.0', 13, 1) |
|
1561 GO TO 750 |
|
1562 C |
|
1563 714 WRITE (XERN3, '(1P,D15.6)') TOUT |
|
1564 WRITE (XERN4, '(1P,D15.6)') T |
|
1565 CALL XERMSG ('SLATEC', 'DDASSL', |
|
1566 * 'TOUT = ' // XERN3 // ' TOO CLOSE TO T = ' // XERN4 // |
|
1567 * ' TO START INTEGRATION', 14, 1) |
|
1568 GO TO 750 |
|
1569 C |
|
1570 715 WRITE (XERN3, '(1P,D15.6)') TSTOP |
|
1571 WRITE (XERN4, '(1P,D15.6)') T |
|
1572 CALL XERMSG ('SLATEC', 'DDASSL', |
|
1573 * 'INFO(4)=1 AND TSTOP = ' // XERN3 // ' BEHIND T = ' // XERN4, |
|
1574 * 15, 1) |
|
1575 GO TO 750 |
|
1576 C |
|
1577 717 WRITE (XERN1, '(I8)') IWORK(LML) |
|
1578 CALL XERMSG ('SLATEC', 'DDASSL', |
|
1579 * 'ML = ' // XERN1 // ' ILLEGAL. EITHER .LT. 0 OR .GT. NEQ', |
|
1580 * 17, 1) |
|
1581 GO TO 750 |
|
1582 C |
|
1583 718 WRITE (XERN1, '(I8)') IWORK(LMU) |
|
1584 CALL XERMSG ('SLATEC', 'DDASSL', |
|
1585 * 'MU = ' // XERN1 // ' ILLEGAL. EITHER .LT. 0 OR .GT. NEQ', |
|
1586 * 18, 1) |
|
1587 GO TO 750 |
|
1588 C |
|
1589 719 WRITE (XERN3, '(1P,D15.6)') TOUT |
|
1590 CALL XERMSG ('SLATEC', 'DDASSL', |
|
1591 * 'TOUT = T = ' // XERN3, 19, 1) |
|
1592 GO TO 750 |
|
1593 C |
|
1594 750 IDID=-33 |
|
1595 IF(INFO(1).EQ.-1) THEN |
|
1596 CALL XERMSG ('SLATEC', 'DDASSL', |
|
1597 * 'REPEATED OCCURRENCES OF ILLEGAL INPUT$$' // |
|
1598 * 'RUN TERMINATED. APPARENT INFINITE LOOP', -999, 2) |
|
1599 ENDIF |
|
1600 C |
|
1601 INFO(1)=-1 |
|
1602 RETURN |
|
1603 C-----------END OF SUBROUTINE DDASSL------------------------------------ |
|
1604 END |