GNU Octave  4.0.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
ddasrt.f
Go to the documentation of this file.
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?
390 C Yes - Set INFO(12) = 0
391 C No - Set INFO(12) = 1,
392 C and define the maximum number of steps
393 C by setting IWORK(21)=MXSTEP
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
445 C LRW .GE. 50+(MAXORD+4)*NEQ+NEQ**2+3*NG
446 C for the full (dense) JACOBIAN case (when INFO(6)=0), or
447 C LRW .GE. 50+(MAXORD+4)*NEQ+(2*ML+MU+1)*NEQ+3*NG
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
451 C +2*(NEQ/(ML+MU+1)+1)+3*NG
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.
459 C you must have LIW .GE. 21+NEQ
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,
862 C XERRWD,D1MACH
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,
880  * lirfnd=18, lmxstp=21, lipvt=22, ljcalc=5, lphase=6, lk=7,
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
919  lenrw=50+(iwork(lmxord)+4)*neq+lenpd+3*ng
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
932  lenrw=50+(iwork(lmxord)+4)*neq+lenpd+2*msave+3*ng
933  go to 60
934 50 iwork(lmtype)=4
935  lenrw=50+(iwork(lmxord)+4)*neq+lenpd+3*ng
936 C
937 C CHECK LENGTHS OF RWORK AND IWORK
938 60 leniw=21+neq
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)
958  IF(mxstp.LT.0)go to 716
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'
986  CALL xerrwd(msg,49,201,0,0,0,0,0,0.0d0,0.0d0)
987  msg = 'DASRT-- VALUE (=I1) OF IDID AND NO APPROPRIATE'
988  CALL xerrwd(msg,47,202,0,1,idid,0,0,0.0d0,0.0d0)
989  msg = 'DASRT-- ACTION WAS TAKEN. RUN TERMINATED'
990  CALL xerrwd(msg,41,203,1,0,0,0,0,0.0d0,0.0d0)
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'
1396  CALL xerrwd(msg,38,610,0,0,0,0,1,tn,0.0d0)
1397  msg = 'DASRT-- TAKEN ON THIS CALL BEFORE REACHING TOUT'
1398  CALL xerrwd(msg,48,611,0,0,0,0,0,0.0d0,0.0d0)
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'
1403  CALL xerrwd(msg,47,620,0,0,0,0,1,tn,0.0d0)
1404  msg = 'DASRT-- FOR PRECISION OF MACHINE. RTOL AND ATOL'
1405  CALL xerrwd(msg,48,621,0,0,0,0,0,0.0d0,0.0d0)
1406  msg = 'DASRT-- WERE INCREASED TO APPROPRIATE VALUES'
1407  CALL xerrwd(msg,45,622,0,0,0,0,0,0.0d0,0.0d0)
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'
1412  CALL xerrwd(msg,38,630,0,0,0,0,1,tn,0.0d0)
1413  msg = .LE.'DASRT-- HAS BECOME 0.0'
1414  CALL xerrwd(msg,28,631,0,0,0,0,0,0.0d0,0.0d0)
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'
1419  CALL xerrwd(msg,44,640,0,0,0,0,2,tn,h)
1420  msg='DASRT-- ERROR TEST FAILED REPEATEDLY OR WITH ABS(H)=HMIN'
1421  CALL xerrwd(msg,57,641,0,0,0,0,0,0.0d0,0.0d0)
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'
1426  CALL xerrwd(msg,44,650,0,0,0,0,2,tn,h)
1427  msg = 'DASRT-- CORRECTOR FAILED TO CONVERGE REPEATEDLY'
1428  CALL xerrwd(msg,48,651,0,0,0,0,0,0.0d0,0.0d0)
1429  msg = 'DASRT-- OR WITH ABS(H)=HMIN'
1430  CALL xerrwd(msg,28,652,0,0,0,0,0,0.0d0,0.0d0)
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'
1435  CALL xerrwd(msg,44,660,0,0,0,0,2,tn,h)
1436  msg = 'DASRT-- ITERATION MATRIX IS SINGULAR'
1437  CALL xerrwd(msg,37,661,0,0,0,0,0,0.0d0,0.0d0)
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'
1442  CALL xerrwd(msg,44,670,0,0,0,0,2,tn,h)
1443  msg = 'DASRT-- CORRECTOR COULD NOT CONVERGE. ALSO, THE'
1444  CALL xerrwd(msg,49,671,0,0,0,0,0,0.0d0,0.0d0)
1445  msg = 'DASRT-- ERROR TEST FAILED REPEATEDLY.'
1446  CALL xerrwd(msg,38,672,0,0,0,0,0,0.0d0,0.0d0)
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'
1451  CALL xerrwd(msg,44,675,0,0,0,0,2,tn,h)
1452  msg = 'DASRT-- CORRECTOR COULD NOT CONVERGE BECAUSE'
1453  CALL xerrwd(msg,45,676,0,0,0,0,0,0.0d0,0.0d0)
1454  msg = 'DASRT-- IRES WAS EQUAL TO MINUS ONE'
1455  CALL xerrwd(msg,36,677,0,0,0,0,0,0.0d0,0.0d0)
1456  go to 690
1457 C
1458 C FAILURE BECAUSE IRES = -2
1459 680 msg = 'DASRT-- AT T (=R1) AND STEPSIZE H (=R2)'
1460  CALL xerrwd(msg,40,680,0,0,0,0,2,tn,h)
1461  msg = 'DASRT-- IRES WAS EQUAL TO MINUS TWO'
1462  CALL xerrwd(msg,36,681,0,0,0,0,0,0.0d0,0.0d0)
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'
1467  CALL xerrwd(msg,44,685,0,0,0,0,2,tn,ho)
1468  msg = 'DASRT-- INITIAL YPRIME COULD NOT BE COMPUTED'
1469  CALL xerrwd(msg,45,686,0,0,0,0,0,0.0d0,0.0d0)
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'
1486  CALL xerrwd(msg,55,1,0,0,0,0,0,0.0d0,0.0d0)
1487  go to 750
1488 702 msg = .LE.'DASRT-- NEQ (=I1) 0'
1489  CALL xerrwd(msg,25,2,0,1,neq,0,0,0.0d0,0.0d0)
1490  go to 750
1491 703 msg = 'DASRT-- MAXORD (=I1) NOT IN RANGE'
1492  CALL xerrwd(msg,34,3,0,1,mxord,0,0,0.0d0,0.0d0)
1493  go to 750
1494 704 msg='DASRT-- RWORK LENGTH NEEDED, LENRW (=I1), EXCEEDS LRW (=I2)'
1495  CALL xerrwd(msg,60,4,0,2,lenrw,lrw,0,0.0d0,0.0d0)
1496  go to 750
1497 705 msg='DASRT-- IWORK LENGTH NEEDED, LENIW (=I1), EXCEEDS LIW (=I2)'
1498  CALL xerrwd(msg,60,5,0,2,leniw,liw,0,0.0d0,0.0d0)
1499  go to 750
1500 706 msg = .LT.'DASRT-- SOME ELEMENT OF RTOL IS 0'
1501  CALL xerrwd(msg,39,6,0,0,0,0,0,0.0d0,0.0d0)
1502  go to 750
1503 707 msg = .LT.'DASRT-- SOME ELEMENT OF ATOL IS 0'
1504  CALL xerrwd(msg,39,7,0,0,0,0,0,0.0d0,0.0d0)
1505  go to 750
1506 708 msg = 'DASRT-- ALL ELEMENTS OF RTOL AND ATOL ARE ZERO'
1507  CALL xerrwd(msg,47,8,0,0,0,0,0,0.0d0,0.0d0)
1508  go to 750
1509 709 msg='DASRT-- INFO(4) = 1 AND TSTOP (=R1) BEHIND TOUT (=R2)'
1510  CALL xerrwd(msg,54,9,0,0,0,0,2,tstop,tout)
1511  go to 750
1512 710 msg = .LT.'DASRT-- HMAX (=R1) 0.0'
1513  CALL xerrwd(msg,28,10,0,0,0,0,1,hmax,0.0d0)
1514  go to 750
1515 711 msg = 'DASRT-- TOUT (=R1) BEHIND T (=R2)'
1516  CALL xerrwd(msg,34,11,0,0,0,0,2,tout,t)
1517  go to 750
1518 712 msg = 'DASRT-- INFO(8)=1 AND H0=0.0'
1519  CALL xerrwd(msg,29,12,0,0,0,0,0,0.0d0,0.0d0)
1520  go to 750
1521 713 msg = .LE.'DASRT-- SOME ELEMENT OF WT IS 0.0'
1522  CALL xerrwd(msg,39,13,0,0,0,0,0,0.0d0,0.0d0)
1523  go to 750
1524 714 msg='DASRT-- TOUT (=R1) TOO CLOSE TO T (=R2) TO START INTEGRATION'
1525  CALL xerrwd(msg,60,14,0,0,0,0,2,tout,t)
1526  go to 750
1527 715 msg = 'DASRT-- INFO(4)=1 AND TSTOP (=R1) BEHIND T (=R2)'
1528  CALL xerrwd(msg,49,15,0,0,0,0,2,tstop,t)
1529  go to 750
1530 716 msg = .LT.'DASRT-- INFO(12)=1 AND MXSTP (=I1) 0'
1531  CALL xerrwd(msg,42,16,0,1,iwork(lmxstp),0,0,0.0d0,0.0d0)
1532  go to 750
1533 717 msg = .LT..GT.'DASRT-- ML (=I1) ILLEGAL. EITHER 0 OR NEQ'
1534  CALL xerrwd(msg,52,17,0,1,iwork(lml),0,0,0.0d0,0.0d0)
1535  go to 750
1536 718 msg = .LT..GT.'DASRT-- MU (=I1) ILLEGAL. EITHER 0 OR NEQ'
1537  CALL xerrwd(msg,52,18,0,1,iwork(lmu),0,0,0.0d0,0.0d0)
1538  go to 750
1539 719 msg = 'DASRT-- TOUT (=R1) IS EQUAL TO T (=R2)'
1540  CALL xerrwd(msg,39,19,0,0,0,0,2,tout,t)
1541  go to 750
1542 730 msg = .LT.'DASRT-- NG (=I1) 0'
1543  CALL xerrwd(msg,24,30,1,1,ng,0,0,0.0d0,0.0d0)
1544  go to 750
1545 732 msg = 'DASRT-- ONE OR MORE COMPONENTS OF G HAS A ROOT'
1546  CALL xerrwd(msg,47,32,1,0,0,0,0,0.0d0,0.0d0)
1547  msg = ' TOO NEAR TO THE INITIAL POINT'
1548  CALL xerrwd(msg,38,32,1,0,0,0,0,0.0d0,0.0d0)
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'
1554  CALL xerrwd(msg,46,801,0,0,0,0,0,0.0d0,0.0d0)
1555 770 msg = 'DASRT-- RUN TERMINATED. APPARENT INFINITE LOOP'
1556  CALL xerrwd(msg,47,802,1,0,0,0,0,0.0d0,0.0d0)
1557  RETURN
1558 C-----------END OF SUBROUTINE DDASRT------------------------------------
1559  END
subroutine ddasrt(RES, NEQ, T, Y, YPRIME, TOUT, INFO, RTOL, ATOL, IDID, RWORK, LRW, IWORK, LIW, RPAR, IPAR, JAC, G, NG, JROOT)
Definition: ddasrt.f:1
std::string dimension(void) const
subroutine xerrwd(MSG, NMES, NERR, LEVEL, NI, I1, I2, NR, R1, R2)
Definition: xerrwd.f:3
subroutine ddatrp(X, XOUT, YOUT, YPOUT, NEQ, KOLD, PHI, PSI)
Definition: ddatrp.f:1
double precision function d1mach(i)
Definition: d1mach.f:1
subroutine ddawts(NEQ, IWT, RTOL, ATOL, Y, WT, RPAR, IPAR)
Definition: ddawts.f:1
double precision function ddanrm(NEQ, V, WT, RPAR, IPAR)
Definition: ddanrm.f:1
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:1
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:1
T abs(T x)
Definition: pr-output.cc:3062
octave_value lgamma(void) const
Definition: ov.h:1189
subroutine drchek(JOB, G, NG, NEQ, TN, TOUT, Y, YP, PHI, PSI, KOLD, G0, G1, GX, JROOT, IRT, UROUND, INFO3, RWORK, IWORK, RPAR, IPAR)
Definition: drchek.f:1