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