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