GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
ddassl.f
Go to the documentation of this file.
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. 21+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,
946  * lmxstp, lipvt,
947  * ljcalc, lk, lkold, liwm, lml, lmtype, lmu, lmxord, lnje, lnpd,
948  * lnre, lns, lnst, lnstl, lpd, lphase, lphi, lpsi, lround, ls,
949  * lsigma, ltn, ltstop, lwm, lwt, mband, msave, mxord, npd, ntemp,
950  * nzflg
951  double precision
952  * atoli, h, hmax, hmin, ho, r, rh, rtoli, tdist, tn, tnext,
953  * tstop, uround, ypnorm
954  LOGICAL DONE
955 C Auxiliary variables for conversion of values to be included in
956 C error messages.
957  CHARACTER*8 XERN1, XERN2
958  CHARACTER*16 XERN3, XERN4
959 C
960 C SET POINTERS INTO IWORK
961  parameter(lml=1, lmu=2, lmxord=3, lmtype=4, lnst=11,
962  * lnre=12, lnje=13, letf=14, lctf=15, lnpd=16, lmxstp=21,
963  * lipvt=22, ljcalc=5, lphase=6, lk=7, lkold=8,
964  * lns=9, lnstl=10, liwm=1)
965 C
966 C SET RELATIVE OFFSET INTO RWORK
967  parameter(npd=1)
968 C
969 C SET POINTERS INTO RWORK
970  parameter(ltstop=1, lhmax=2, lh=3, ltn=4,
971  * lcj=5, lcjold=6, lhold=7, ls=8, lround=9,
972  * lalpha=11, lbeta=17, lgamma=23,
973  * lpsi=29, lsigma=35, ldelta=41)
974 C
975 C***FIRST EXECUTABLE STATEMENT DDASSL
976  IF(info(1).NE.0)GO TO 100
977 C
978 C-----------------------------------------------------------------------
979 C THIS BLOCK IS EXECUTED FOR THE INITIAL CALL ONLY.
980 C IT CONTAINS CHECKING OF INPUTS AND INITIALIZATIONS.
981 C-----------------------------------------------------------------------
982 C
983 C FIRST CHECK INFO ARRAY TO MAKE SURE ALL ELEMENTS OF INFO
984 C ARE EITHER ZERO OR ONE.
985  DO 10 i=2,11
986  IF(info(i).NE.0.AND.info(i).NE.1)GO TO 701
987 10 CONTINUE
988 C
989  IF(neq.LE.0)GO TO 702
990 C
991 C CHECK AND COMPUTE MAXIMUM ORDER
992  mxord=5
993  IF(info(9).EQ.0)GO TO 20
994  mxord=iwork(lmxord)
995  IF(mxord.LT.1.OR.mxord.GT.5)GO TO 703
996 20 iwork(lmxord)=mxord
997 C
998 C COMPUTE MTYPE,LENPD,LENRW.CHECK ML AND MU.
999  IF(info(6).NE.0)GO TO 40
1000  lenpd=neq**2
1001  lenrw=40+(iwork(lmxord)+4)*neq+lenpd
1002  IF(info(5).NE.0)GO TO 30
1003  iwork(lmtype)=2
1004  GO TO 60
1005 30 iwork(lmtype)=1
1006  GO TO 60
1007 40 IF(iwork(lml).LT.0.OR.iwork(lml).GE.neq)GO TO 717
1008  IF(iwork(lmu).LT.0.OR.iwork(lmu).GE.neq)GO TO 718
1009  lenpd=(2*iwork(lml)+iwork(lmu)+1)*neq
1010  IF(info(5).NE.0)GO TO 50
1011  iwork(lmtype)=5
1012  mband=iwork(lml)+iwork(lmu)+1
1013  msave=(neq/mband)+1
1014  lenrw=40+(iwork(lmxord)+4)*neq+lenpd+2*msave
1015  GO TO 60
1016 50 iwork(lmtype)=4
1017  lenrw=40+(iwork(lmxord)+4)*neq+lenpd
1018 C
1019 C CHECK LENGTHS OF RWORK AND IWORK
1020 60 leniw=21+neq
1021  iwork(lnpd)=lenpd
1022  IF(lrw.LT.lenrw)GO TO 704
1023  IF(liw.LT.leniw)GO TO 705
1024 C
1025 C CHECK TO SEE THAT TOUT IS DIFFERENT FROM T
1026  IF(tout .EQ. t)GO TO 719
1027 C
1028 C CHECK HMAX
1029  IF(info(7).EQ.0)GO TO 70
1030  hmax=rwork(lhmax)
1031  IF(hmax.LE.0.0d0)GO TO 710
1032 70 CONTINUE
1033 C
1034 C CHECK AND COMPUTE MAXIMUM STEPS
1035  mxstp=500
1036  IF(info(12).EQ.0)GO TO 80
1037  mxstp=iwork(lmxstp)
1038  IF(mxstp.LT.0)GO TO 716
1039 80 iwork(lmxstp)=mxstp
1040 C
1041 C INITIALIZE COUNTERS
1042  iwork(lnst)=0
1043  iwork(lnre)=0
1044  iwork(lnje)=0
1045 C
1046  iwork(lnstl)=0
1047  idid=1
1048  GO TO 200
1049 C
1050 C-----------------------------------------------------------------------
1051 C THIS BLOCK IS FOR CONTINUATION CALLS
1052 C ONLY. HERE WE CHECK INFO(1),AND IF THE
1053 C LAST STEP WAS INTERRUPTED WE CHECK WHETHER
1054 C APPROPRIATE ACTION WAS TAKEN.
1055 C-----------------------------------------------------------------------
1056 C
1057 100 CONTINUE
1058  IF(info(1).EQ.1)GO TO 110
1059  IF(info(1).NE.-1)GO TO 701
1060 C
1061 C IF WE ARE HERE, THE LAST STEP WAS INTERRUPTED
1062 C BY AN ERROR CONDITION FROM DDASTP,AND
1063 C APPROPRIATE ACTION WAS NOT TAKEN. THIS
1064 C IS A FATAL ERROR.
1065  WRITE (xern1, '(I8)') idid
1066  CALL xermsg ('SLATEC', 'DDASSL',
1067  * 'THE LAST STEP TERMINATED WITH A NEGATIVE VALUE OF IDID = ' //
1068  * xern1 // ' AND NO APPROPRIATE ACTION WAS TAKEN. ' //
1069  * 'RUN TERMINATED', -998, 2)
1070  RETURN
1071 110 CONTINUE
1072  iwork(lnstl)=iwork(lnst)
1073 C
1074 C-----------------------------------------------------------------------
1075 C THIS BLOCK IS EXECUTED ON ALL CALLS.
1076 C THE ERROR TOLERANCE PARAMETERS ARE
1077 C CHECKED, AND THE WORK ARRAY POINTERS
1078 C ARE SET.
1079 C-----------------------------------------------------------------------
1080 C
1081 200 CONTINUE
1082 C CHECK RTOL,ATOL
1083  nzflg=0
1084  rtoli=rtol(1)
1085  atoli=atol(1)
1086  DO 210 i=1,neq
1087  IF(info(2).EQ.1)rtoli=rtol(i)
1088  IF(info(2).EQ.1)atoli=atol(i)
1089  IF(rtoli.GT.0.0d0.OR.atoli.GT.0.0d0)nzflg=1
1090  IF(rtoli.LT.0.0d0)GO TO 706
1091  IF(atoli.LT.0.0d0)GO TO 707
1092 210 CONTINUE
1093  IF(nzflg.EQ.0)GO TO 708
1094 C
1095 C SET UP RWORK STORAGE.IWORK STORAGE IS FIXED
1096 C IN DATA STATEMENT.
1097  le=ldelta+neq
1098  lwt=le+neq
1099  lphi=lwt+neq
1100  lpd=lphi+(iwork(lmxord)+1)*neq
1101  lwm=lpd
1102  ntemp=npd+iwork(lnpd)
1103  IF(info(1).EQ.1)GO TO 400
1104 C
1105 C-----------------------------------------------------------------------
1106 C THIS BLOCK IS EXECUTED ON THE INITIAL CALL
1107 C ONLY. SET THE INITIAL STEP SIZE, AND
1108 C THE ERROR WEIGHT VECTOR, AND PHI.
1109 C COMPUTE INITIAL YPRIME, IF NECESSARY.
1110 C-----------------------------------------------------------------------
1111 C
1112  tn=t
1113  idid=1
1114 C
1115 C SET ERROR WEIGHT VECTOR WT
1116  CALL ddawts(neq,info(2),rtol,atol,y,rwork(lwt),rpar,ipar)
1117  DO 305 i = 1,neq
1118  IF(rwork(lwt+i-1).LE.0.0d0) GO TO 713
1119 305 CONTINUE
1120 C
1121 C COMPUTE UNIT ROUNDOFF AND HMIN
1122  uround = d1mach(4)
1123  rwork(lround) = uround
1124  hmin = 4.0d0*uround*max(abs(t),abs(tout))
1125 C
1126 C CHECK INITIAL INTERVAL TO SEE THAT IT IS LONG ENOUGH
1127  tdist = abs(tout - t)
1128  IF(tdist .LT. hmin) GO TO 714
1129 C
1130 C CHECK HO, IF THIS WAS INPUT
1131  IF (info(8) .EQ. 0) GO TO 310
1132  ho = rwork(lh)
1133  IF ((tout - t)*ho .LT. 0.0d0) GO TO 711
1134  IF (ho .EQ. 0.0d0) GO TO 712
1135  GO TO 320
1136 310 CONTINUE
1137 C
1138 C COMPUTE INITIAL STEPSIZE, TO BE USED BY EITHER
1139 C DDASTP OR DDAINI, DEPENDING ON INFO(11)
1140  ho = 0.001d0*tdist
1141  ypnorm = ddanrm(neq,yprime,rwork(lwt),rpar,ipar)
1142  IF (ypnorm .GT. 0.5d0/ho) ho = 0.5d0/ypnorm
1143  ho = sign(ho,tout-t)
1144 C ADJUST HO IF NECESSARY TO MEET HMAX BOUND
1145 320 IF (info(7) .EQ. 0) GO TO 330
1146  rh = abs(ho)/rwork(lhmax)
1147  IF (rh .GT. 1.0d0) ho = ho/rh
1148 C COMPUTE TSTOP, IF APPLICABLE
1149 330 IF (info(4) .EQ. 0) GO TO 340
1150  tstop = rwork(ltstop)
1151  IF ((tstop - t)*ho .LT. 0.0d0) GO TO 715
1152  IF ((t + ho - tstop)*ho .GT. 0.0d0) ho = tstop - t
1153  IF ((tstop - tout)*ho .LT. 0.0d0) GO TO 709
1154 C
1155 C COMPUTE INITIAL DERIVATIVE, UPDATING TN AND Y, IF APPLICABLE
1156 340 IF (info(11) .EQ. 0) GO TO 350
1157  CALL ddaini(tn,y,yprime,neq,
1158  * res,jac,ho,rwork(lwt),idid,rpar,ipar,
1159  * rwork(lphi),rwork(ldelta),rwork(le),
1160  * rwork(lwm),iwork(liwm),hmin,rwork(lround),
1161  * info(10),ntemp)
1162  IF (idid .LT. 0) GO TO 390
1163 C
1164 C LOAD H WITH HO. STORE H IN RWORK(LH)
1165 350 h = ho
1166  rwork(lh) = h
1167 C
1168 C LOAD Y AND H*YPRIME INTO PHI(*,1) AND PHI(*,2)
1169  itemp = lphi + neq
1170  DO 370 i = 1,neq
1171  rwork(lphi + i - 1) = y(i)
1172 370 rwork(itemp + i - 1) = h*yprime(i)
1173 C
1174 390 GO TO 500
1175 C
1176 C-------------------------------------------------------
1177 C THIS BLOCK IS FOR CONTINUATION CALLS ONLY. ITS
1178 C PURPOSE IS TO CHECK STOP CONDITIONS BEFORE
1179 C TAKING A STEP.
1180 C ADJUST H IF NECESSARY TO MEET HMAX BOUND
1181 C-------------------------------------------------------
1182 C
1183 400 CONTINUE
1184  uround=rwork(lround)
1185  done = .false.
1186  tn=rwork(ltn)
1187  h=rwork(lh)
1188  IF(info(7) .EQ. 0) GO TO 410
1189  rh = abs(h)/rwork(lhmax)
1190  IF(rh .GT. 1.0d0) h = h/rh
1191 410 CONTINUE
1192  IF(t .EQ. tout) GO TO 719
1193  IF((t - tout)*h .GT. 0.0d0) GO TO 711
1194  IF(info(4) .EQ. 1) GO TO 430
1195  IF(info(3) .EQ. 1) GO TO 420
1196  IF((tn-tout)*h.LT.0.0d0)GO TO 490
1197  CALL ddatrp(tn,tout,y,yprime,neq,iwork(lkold),
1198  * rwork(lphi),rwork(lpsi))
1199  t=tout
1200  idid = 3
1201  done = .true.
1202  GO TO 490
1203 420 IF((tn-t)*h .LE. 0.0d0) GO TO 490
1204  IF((tn - tout)*h .GT. 0.0d0) GO TO 425
1205  CALL ddatrp(tn,tn,y,yprime,neq,iwork(lkold),
1206  * rwork(lphi),rwork(lpsi))
1207  t = tn
1208  idid = 1
1209  done = .true.
1210  GO TO 490
1211 425 CONTINUE
1212  CALL ddatrp(tn,tout,y,yprime,neq,iwork(lkold),
1213  * rwork(lphi),rwork(lpsi))
1214  t = tout
1215  idid = 3
1216  done = .true.
1217  GO TO 490
1218 430 IF(info(3) .EQ. 1) GO TO 440
1219  tstop=rwork(ltstop)
1220  IF((tn-tstop)*h.GT.0.0d0) GO TO 715
1221  IF((tstop-tout)*h.LT.0.0d0)GO TO 709
1222  IF((tn-tout)*h.LT.0.0d0)GO TO 450
1223  CALL ddatrp(tn,tout,y,yprime,neq,iwork(lkold),
1224  * rwork(lphi),rwork(lpsi))
1225  t=tout
1226  idid = 3
1227  done = .true.
1228  GO TO 490
1229 440 tstop = rwork(ltstop)
1230  IF((tn-tstop)*h .GT. 0.0d0) GO TO 715
1231  IF((tstop-tout)*h .LT. 0.0d0) GO TO 709
1232  IF((tn-t)*h .LE. 0.0d0) GO TO 450
1233  IF((tn - tout)*h .GT. 0.0d0) GO TO 445
1234  CALL ddatrp(tn,tn,y,yprime,neq,iwork(lkold),
1235  * rwork(lphi),rwork(lpsi))
1236  t = tn
1237  idid = 1
1238  done = .true.
1239  GO TO 490
1240 445 CONTINUE
1241  CALL ddatrp(tn,tout,y,yprime,neq,iwork(lkold),
1242  * rwork(lphi),rwork(lpsi))
1243  t = tout
1244  idid = 3
1245  done = .true.
1246  GO TO 490
1247 450 CONTINUE
1248 C CHECK WHETHER WE ARE WITHIN ROUNDOFF OF TSTOP
1249  IF(abs(tn-tstop).GT.100.0d0*uround*
1250  * (abs(tn)+abs(h)))GO TO 460
1251  CALL ddatrp(tn,tstop,y,yprime,neq,iwork(lkold),
1252  * rwork(lphi),rwork(lpsi))
1253  idid=2
1254  t=tstop
1255  done = .true.
1256  GO TO 490
1257 460 tnext=tn+h
1258  IF((tnext-tstop)*h.LE.0.0d0)GO TO 490
1259  h=tstop-tn
1260  rwork(lh)=h
1261 C
1262 490 IF (done) GO TO 580
1263 C
1264 C-------------------------------------------------------
1265 C THE NEXT BLOCK CONTAINS THE CALL TO THE
1266 C ONE-STEP INTEGRATOR DDASTP.
1267 C THIS IS A LOOPING POINT FOR THE INTEGRATION STEPS.
1268 C CHECK FOR TOO MANY STEPS.
1269 C UPDATE WT.
1270 C CHECK FOR TOO MUCH ACCURACY REQUESTED.
1271 C COMPUTE MINIMUM STEPSIZE.
1272 C-------------------------------------------------------
1273 C
1274 500 CONTINUE
1275 C CHECK FOR FAILURE TO COMPUTE INITIAL YPRIME
1276  IF (idid .EQ. -12) GO TO 527
1277 C
1278 C CHECK FOR TOO MANY STEPS
1279  IF((iwork(lnst)-iwork(lnstl)).LT.iwork(lmxstp))
1280  * GO TO 510
1281  idid=-1
1282  GO TO 527
1283 C
1284 C UPDATE WT
1285 510 CALL ddawts(neq,info(2),rtol,atol,rwork(lphi),
1286  * rwork(lwt),rpar,ipar)
1287  DO 520 i=1,neq
1288  IF(rwork(i+lwt-1).GT.0.0d0)GO TO 520
1289  idid=-3
1290  GO TO 527
1291 520 CONTINUE
1292 C
1293 C TEST FOR TOO MUCH ACCURACY REQUESTED.
1294  r=ddanrm(neq,rwork(lphi),rwork(lwt),rpar,ipar)*
1295  * 100.0d0*uround
1296  IF(r.LE.1.0d0)GO TO 525
1297 C MULTIPLY RTOL AND ATOL BY R AND RETURN
1298  IF(info(2).EQ.1)GO TO 523
1299  rtol(1)=r*rtol(1)
1300  atol(1)=r*atol(1)
1301  idid=-2
1302  GO TO 527
1303 523 DO 524 i=1,neq
1304  rtol(i)=r*rtol(i)
1305 524 atol(i)=r*atol(i)
1306  idid=-2
1307  GO TO 527
1308 525 CONTINUE
1309 C
1310 C COMPUTE MINIMUM STEPSIZE
1311  hmin=4.0d0*uround*max(abs(tn),abs(tout))
1312 C
1313 C TEST H VS. HMAX
1314  IF (info(7) .EQ. 0) GO TO 526
1315  rh = abs(h)/rwork(lhmax)
1316  IF (rh .GT. 1.0d0) h = h/rh
1317 526 CONTINUE
1318 C
1319  CALL ddastp(tn,y,yprime,neq,
1320  * res,jac,h,rwork(lwt),info(1),idid,rpar,ipar,
1321  * rwork(lphi),rwork(ldelta),rwork(le),
1322  * rwork(lwm),iwork(liwm),
1323  * rwork(lalpha),rwork(lbeta),rwork(lgamma),
1324  * rwork(lpsi),rwork(lsigma),
1325  * rwork(lcj),rwork(lcjold),rwork(lhold),
1326  * rwork(ls),hmin,rwork(lround),
1327  * iwork(lphase),iwork(ljcalc),iwork(lk),
1328  * iwork(lkold),iwork(lns),info(10),ntemp)
1329 527 IF(idid.LT.0)GO TO 600
1330 C
1331 C--------------------------------------------------------
1332 C THIS BLOCK HANDLES THE CASE OF A SUCCESSFUL RETURN
1333 C FROM DDASTP (IDID=1). TEST FOR STOP CONDITIONS.
1334 C--------------------------------------------------------
1335 C
1336  IF(info(4).NE.0)GO TO 540
1337  IF(info(3).NE.0)GO TO 530
1338  IF((tn-tout)*h.LT.0.0d0)GO TO 500
1339  CALL ddatrp(tn,tout,y,yprime,neq,
1340  * iwork(lkold),rwork(lphi),rwork(lpsi))
1341  idid=3
1342  t=tout
1343  GO TO 580
1344 530 IF((tn-tout)*h.GE.0.0d0)GO TO 535
1345  t=tn
1346  idid=1
1347  GO TO 580
1348 535 CALL ddatrp(tn,tout,y,yprime,neq,
1349  * iwork(lkold),rwork(lphi),rwork(lpsi))
1350  idid=3
1351  t=tout
1352  GO TO 580
1353 540 IF(info(3).NE.0)GO TO 550
1354  IF((tn-tout)*h.LT.0.0d0)GO TO 542
1355  CALL ddatrp(tn,tout,y,yprime,neq,
1356  * iwork(lkold),rwork(lphi),rwork(lpsi))
1357  t=tout
1358  idid=3
1359  GO TO 580
1360 542 IF(abs(tn-tstop).LE.100.0d0*uround*
1361  * (abs(tn)+abs(h)))GO TO 545
1362  tnext=tn+h
1363  IF((tnext-tstop)*h.LE.0.0d0)GO TO 500
1364  h=tstop-tn
1365  GO TO 500
1366 545 CALL ddatrp(tn,tstop,y,yprime,neq,
1367  * iwork(lkold),rwork(lphi),rwork(lpsi))
1368  idid=2
1369  t=tstop
1370  GO TO 580
1371 550 IF((tn-tout)*h.GE.0.0d0)GO TO 555
1372  IF(abs(tn-tstop).LE.100.0d0*uround*(abs(tn)+abs(h)))GO TO 552
1373  t=tn
1374  idid=1
1375  GO TO 580
1376 552 CALL ddatrp(tn,tstop,y,yprime,neq,
1377  * iwork(lkold),rwork(lphi),rwork(lpsi))
1378  idid=2
1379  t=tstop
1380  GO TO 580
1381 555 CALL ddatrp(tn,tout,y,yprime,neq,
1382  * iwork(lkold),rwork(lphi),rwork(lpsi))
1383  t=tout
1384  idid=3
1385  GO TO 580
1386 C
1387 C--------------------------------------------------------
1388 C ALL SUCCESSFUL RETURNS FROM DDASSL ARE MADE FROM
1389 C THIS BLOCK.
1390 C--------------------------------------------------------
1391 C
1392 580 CONTINUE
1393  rwork(ltn)=tn
1394  rwork(lh)=h
1395  RETURN
1396 C
1397 C-----------------------------------------------------------------------
1398 C THIS BLOCK HANDLES ALL UNSUCCESSFUL
1399 C RETURNS OTHER THAN FOR ILLEGAL INPUT.
1400 C-----------------------------------------------------------------------
1401 C
1402 600 CONTINUE
1403  itemp=-idid
1404  GO TO (610,620,630,690,690,640,650,660,670,675,
1405  * 680,685), itemp
1406 C
1407 C THE MAXIMUM NUMBER OF STEPS WAS TAKEN BEFORE
1408 C REACHING TOUT
1409 610 WRITE (xern3, '(1P,D15.6)') tn
1410  CALL xermsg ('SLATEC', 'DDASSL',
1411  * 'AT CURRENT T = ' // xern3 // ' 500 STEPS TAKEN ON THIS ' //
1412  * 'CALL BEFORE REACHING TOUT', idid, 1)
1413  GO TO 690
1414 C
1415 C TOO MUCH ACCURACY FOR MACHINE PRECISION
1416 620 WRITE (xern3, '(1P,D15.6)') tn
1417  CALL xermsg ('SLATEC', 'DDASSL',
1418  * 'AT T = ' // xern3 // ' TOO MUCH ACCURACY REQUESTED FOR ' //
1419  * 'PRECISION OF MACHINE. RTOL AND ATOL WERE INCREASED TO ' //
1420  * 'APPROPRIATE VALUES', idid, 1)
1421  GO TO 690
1422 C
1423 C WT(I) .LE. 0.0 FOR SOME I (NOT AT START OF PROBLEM)
1424 630 WRITE (xern3, '(1P,D15.6)') tn
1425  CALL xermsg ('SLATEC', 'DDASSL',
1426  * 'AT T = ' // xern3 // .LE.' SOME ELEMENT OF WT HAS BECOME ' //
1427  * '0.0', idid, 1)
1428  GO TO 690
1429 C
1430 C ERROR TEST FAILED REPEATEDLY OR WITH H=HMIN
1431 640 WRITE (xern3, '(1P,D15.6)') tn
1432  WRITE (xern4, '(1P,D15.6)') h
1433  CALL xermsg ('SLATEC', 'DDASSL',
1434  * 'AT T = ' // xern3 // ' AND STEPSIZE H = ' // xern4 //
1435  * ' THE ERROR TEST FAILED REPEATEDLY OR WITH ABS(H)=HMIN',
1436  * idid, 1)
1437  GO TO 690
1438 C
1439 C CORRECTOR CONVERGENCE FAILED REPEATEDLY OR WITH H=HMIN
1440 650 WRITE (xern3, '(1P,D15.6)') tn
1441  WRITE (xern4, '(1P,D15.6)') h
1442  CALL xermsg ('SLATEC', 'DDASSL',
1443  * 'AT T = ' // xern3 // ' AND STEPSIZE H = ' // xern4 //
1444  * ' THE CORRECTOR FAILED TO CONVERGE REPEATEDLY OR WITH ' //
1445  * 'ABS(H)=HMIN', idid, 1)
1446  GO TO 690
1447 C
1448 C THE ITERATION MATRIX IS SINGULAR
1449 660 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 ITERATION MATRIX IS SINGULAR', idid, 1)
1454  GO TO 690
1455 C
1456 C CORRECTOR FAILURE PRECEEDED BY ERROR TEST FAILURES.
1457 670 WRITE (xern3, '(1P,D15.6)') tn
1458  WRITE (xern4, '(1P,D15.6)') h
1459  CALL xermsg ('SLATEC', 'DDASSL',
1460  * 'AT T = ' // xern3 // ' AND STEPSIZE H = ' // xern4 //
1461  * ' THE CORRECTOR COULD NOT CONVERGE. ALSO, THE ERROR TEST ' //
1462  * 'FAILED REPEATEDLY.', idid, 1)
1463  GO TO 690
1464 C
1465 C CORRECTOR FAILURE BECAUSE IRES = -1
1466 675 WRITE (xern3, '(1P,D15.6)') tn
1467  WRITE (xern4, '(1P,D15.6)') h
1468  CALL xermsg ('SLATEC', 'DDASSL',
1469  * 'AT T = ' // xern3 // ' AND STEPSIZE H = ' // xern4 //
1470  * ' THE CORRECTOR COULD NOT CONVERGE BECAUSE IRES WAS EQUAL ' //
1471  * 'TO MINUS ONE', idid, 1)
1472  GO TO 690
1473 C
1474 C FAILURE BECAUSE IRES = -2
1475 680 WRITE (xern3, '(1P,D15.6)') tn
1476  WRITE (xern4, '(1P,D15.6)') h
1477  CALL xermsg ('SLATEC', 'DDASSL',
1478  * 'AT T = ' // xern3 // ' AND STEPSIZE H = ' // xern4 //
1479  * ' IRES WAS EQUAL TO MINUS TWO', idid, 1)
1480  GO TO 690
1481 C
1482 C FAILED TO COMPUTE INITIAL YPRIME
1483 685 WRITE (xern3, '(1P,D15.6)') tn
1484  WRITE (xern4, '(1P,D15.6)') ho
1485  CALL xermsg ('SLATEC', 'DDASSL',
1486  * 'AT T = ' // xern3 // ' AND STEPSIZE H = ' // xern4 //
1487  * ' THE INITIAL YPRIME COULD NOT BE COMPUTED', idid, 1)
1488  GO TO 690
1489 C
1490 690 CONTINUE
1491  info(1)=-1
1492  t=tn
1493  rwork(ltn)=tn
1494  rwork(lh)=h
1495  RETURN
1496 C
1497 C-----------------------------------------------------------------------
1498 C THIS BLOCK HANDLES ALL ERROR RETURNS DUE
1499 C TO ILLEGAL INPUT, AS DETECTED BEFORE CALLING
1500 C DDASTP. FIRST THE ERROR MESSAGE ROUTINE IS
1501 C CALLED. IF THIS HAPPENS TWICE IN
1502 C SUCCESSION, EXECUTION IS TERMINATED
1503 C
1504 C-----------------------------------------------------------------------
1505 701 CALL xermsg ('SLATEC', 'DDASSL',
1506  * 'SOME ELEMENT OF INFO VECTOR IS NOT ZERO OR ONE', 1, 1)
1507  GO TO 750
1508 C
1509 702 WRITE (xern1, '(I8)') neq
1510  CALL xermsg ('SLATEC', 'DDASSL',
1511  * 'NEQ = ' // xern1 // .LE.' 0', 2, 1)
1512  GO TO 750
1513 C
1514 703 WRITE (xern1, '(I8)') mxord
1515  CALL xermsg ('SLATEC', 'DDASSL',
1516  * 'MAXORD = ' // xern1 // ' NOT IN RANGE', 3, 1)
1517  GO TO 750
1518 C
1519 704 WRITE (xern1, '(I8)') lenrw
1520  WRITE (xern2, '(I8)') lrw
1521  CALL xermsg ('SLATEC', 'DDASSL',
1522  * 'RWORK LENGTH NEEDED, LENRW = ' // xern1 //
1523  * ', EXCEEDS LRW = ' // xern2, 4, 1)
1524  GO TO 750
1525 C
1526 705 WRITE (xern1, '(I8)') leniw
1527  WRITE (xern2, '(I8)') liw
1528  CALL xermsg ('SLATEC', 'DDASSL',
1529  * 'IWORK LENGTH NEEDED, LENIW = ' // xern1 //
1530  * ', EXCEEDS LIW = ' // xern2, 5, 1)
1531  GO TO 750
1532 C
1533 706 CALL xermsg ('SLATEC', 'DDASSL',
1534  * .LT.'SOME ELEMENT OF RTOL IS 0', 6, 1)
1535  GO TO 750
1536 C
1537 707 CALL xermsg ('SLATEC', 'DDASSL',
1538  * .LT.'SOME ELEMENT OF ATOL IS 0', 7, 1)
1539  GO TO 750
1540 C
1541 708 CALL xermsg ('SLATEC', 'DDASSL',
1542  * 'ALL ELEMENTS OF RTOL AND ATOL ARE ZERO', 8, 1)
1543  GO TO 750
1544 C
1545 709 WRITE (xern3, '(1P,D15.6)') tstop
1546  WRITE (xern4, '(1P,D15.6)') tout
1547  CALL xermsg ('SLATEC', 'DDASSL',
1548  * 'INFO(4) = 1 AND TSTOP = ' // xern3 // ' BEHIND TOUT = ' //
1549  * xern4, 9, 1)
1550  GO TO 750
1551 C
1552 710 WRITE (xern3, '(1P,D15.6)') hmax
1553  CALL xermsg ('SLATEC', 'DDASSL',
1554  * 'HMAX = ' // xern3 // .LT.' 0.0', 10, 1)
1555  GO TO 750
1556 C
1557 711 WRITE (xern3, '(1P,D15.6)') tout
1558  WRITE (xern4, '(1P,D15.6)') t
1559  CALL xermsg ('SLATEC', 'DDASSL',
1560  * 'TOUT = ' // xern3 // ' BEHIND T = ' // xern4, 11, 1)
1561  GO TO 750
1562 C
1563 712 CALL xermsg ('SLATEC', 'DDASSL',
1564  * 'INFO(8)=1 AND H0=0.0', 12, 1)
1565  GO TO 750
1566 C
1567 713 CALL xermsg ('SLATEC', 'DDASSL',
1568  * .LE.'SOME ELEMENT OF WT IS 0.0', 13, 1)
1569  GO TO 750
1570 C
1571 714 WRITE (xern3, '(1P,D15.6)') tout
1572  WRITE (xern4, '(1P,D15.6)') t
1573  CALL xermsg ('SLATEC', 'DDASSL',
1574  * 'TOUT = ' // xern3 // ' TOO CLOSE TO T = ' // xern4 //
1575  * ' TO START INTEGRATION', 14, 1)
1576  GO TO 750
1577 C
1578 715 WRITE (xern3, '(1P,D15.6)') tstop
1579  WRITE (xern4, '(1P,D15.6)') t
1580  CALL xermsg ('SLATEC', 'DDASSL',
1581  * 'INFO(4)=1 AND TSTOP = ' // xern3 // ' BEHIND T = ' // xern4,
1582  * 15, 1)
1583  GO TO 750
1584 C
1585 716 WRITE (xern1, '(I8)') mxstp
1586  CALL xermsg ('SLATEC', 'DDASSL',
1587  * 'INFO(12)=1 AND MXSTP = ' // xern1 // ' ILLEGAL.', 3, 1)
1588  GO TO 750
1589 C
1590 717 WRITE (xern1, '(I8)') iwork(lml)
1591  CALL xermsg ('SLATEC', 'DDASSL',
1592  * 'ML = ' // xern1 // .LT..GT.' ILLEGAL. EITHER 0 OR NEQ',
1593  * 17, 1)
1594  GO TO 750
1595 C
1596 718 WRITE (xern1, '(I8)') iwork(lmu)
1597  CALL xermsg ('SLATEC', 'DDASSL',
1598  * 'MU = ' // xern1 // .LT..GT.' ILLEGAL. EITHER 0 OR NEQ',
1599  * 18, 1)
1600  GO TO 750
1601 C
1602 719 WRITE (xern3, '(1P,D15.6)') tout
1603  CALL xermsg ('SLATEC', 'DDASSL',
1604  * 'TOUT = T = ' // xern3, 19, 1)
1605  GO TO 750
1606 C
1607 750 idid=-33
1608  IF(info(1).EQ.-1) THEN
1609  CALL xermsg ('SLATEC', 'DDASSL',
1610  * 'REPEATED OCCURRENCES OF ILLEGAL INPUT$$' //
1611  * 'RUN TERMINATED. APPARENT INFINITE LOOP', -999, 2)
1612  ENDIF
1613 C
1614  info(1)=-1
1615  RETURN
1616 C-----------END OF SUBROUTINE DDASSL------------------------------------
1617  END
charNDArray max(char d, const charNDArray &m)
Definition: chNDArray.cc:230
double precision function d1mach(i)
Definition: d1mach.f:23
subroutine ddaini(X, Y, YPRIME, NEQ, RES, JAC, H, WT, IDID, RPAR, IPAR, PHI, DELTA, E, WM, IWM, HMIN, UROUND, NONNEG, NTEMP)
Definition: ddaini.f:3
double precision function ddanrm(NEQ, V, WT, RPAR, IPAR)
Definition: ddanrm.f:2
subroutine ddassl(RES, NEQ, T, Y, YPRIME, TOUT, INFO, RTOL, ATOL, IDID, RWORK, LRW, IWORK, LIW, RPAR, IPAR, JAC)
Definition: ddassl.f:3
subroutine ddastp(X, Y, YPRIME, NEQ, RES, JAC, H, WT, JSTART, IDID, RPAR, IPAR, PHI, DELTA, E, WM, IWM, ALPHA, BETA, GAMMA, PSI, SIGMA, CJ, CJOLD, HOLD, S, HMIN, UROUND, IPHASE, JCALC, K, KOLD, NS, NONNEG, NTEMP)
Definition: ddastp.f:5
subroutine ddatrp(X, XOUT, YOUT, YPOUT, NEQ, KOLD, PHI, PSI)
Definition: ddatrp.f:2
subroutine ddawts(NEQ, IWT, RTOL, ATOL, Y, WT, RPAR, IPAR)
Definition: ddawts.f:2
static T abs(T x)
Definition: pr-output.cc:1678
subroutine xermsg(LIBRAR, SUBROU, MESSG, NERR, LEVEL)
Definition: xermsg.f:3