1
SUBROUTINE DDASSL (RES, NEQ, T, Y, YPRIME, TOUT, INFO, RTOL, ATOL,
2
+ IDID, RWORK, LRW, IWORK, LIW, RPAR, IPAR, JAC)
3
C***BEGIN PROLOGUE DDASSL
4
C***PURPOSE This code solves a system of differential/algebraic
5
C equations of the form G(T,Y,YPRIME) = 0.
6
C***LIBRARY SLATEC (DASSL)
8
C***TYPE DOUBLE PRECISION (SDASSL-S, DDASSL-D)
9
C***KEYWORDS DIFFERENTIAL/ALGEBRAIC, BACKWARD DIFFERENTIATION FORMULAS,
10
C IMPLICIT DIFFERENTIAL SYSTEMS
11
C***AUTHOR PETZOLD, LINDA R., (LLNL)
12
C COMPUTING AND MATHEMATICS RESEARCH DIVISION
13
C LAWRENCE LIVERMORE NATIONAL LABORATORY
14
C L - 316, P.O. BOX 808,
15
C LIVERMORE, CA. 94550
21
C INTEGER NEQ, INFO(N), IDID, LRW, LIW, IWORK(LIW), IPAR
22
C DOUBLE PRECISION T, Y(NEQ), YPRIME(NEQ), TOUT, RTOL, ATOL,
25
C CALL DDASSL (RES, NEQ, T, Y, YPRIME, TOUT, INFO, RTOL, ATOL,
26
C * IDID, RWORK, LRW, IWORK, LIW, RPAR, IPAR, JAC)
30
C (In the following, all real arrays should be type DOUBLE PRECISION.)
32
C RES:EXT This is a subroutine which you provide to define the
33
C differential/algebraic system.
35
C NEQ:IN This is the number of equations to be solved.
37
C T:INOUT This is the current value of the independent variable.
39
C Y(*):INOUT This array contains the solution components at T.
41
C YPRIME(*):INOUT This array contains the derivatives of the solution
44
C TOUT:IN This is a point at which a solution is desired.
46
C INFO(N):IN The basic task of the code is to solve the system from T
47
C to TOUT and return an answer at TOUT. INFO is an integer
48
C array which is used to communicate exactly how you want
49
C this task to be carried out. (See below for details.)
50
C N must be greater than or equal to 15.
52
C RTOL,ATOL:INOUT These quantities represent relative and absolute
53
C error tolerances which you provide to indicate how
54
C accurately you wish the solution to be computed. You
55
C may choose them to be both scalars or else both vectors.
56
C Caution: In Fortran 77, a scalar is not the same as an
57
C array of length 1. Some compilers may object
58
C to using scalars for RTOL,ATOL.
60
C IDID:OUT This scalar quantity is an indicator reporting what the
61
C code did. You must monitor this integer variable to
62
C decide what action to take next.
64
C RWORK:WORK A real work array of length LRW which provides the
65
C code with needed storage space.
67
C LRW:IN The length of RWORK. (See below for required length.)
69
C IWORK:WORK An integer work array of length LIW which probides the
70
C code with needed storage space.
72
C LIW:IN The length of IWORK. (See below for required length.)
74
C RPAR,IPAR:IN These are real and integer parameter arrays which
75
C you can use for communication between your calling
76
C program and the RES subroutine (and the JAC subroutine)
78
C JAC:EXT This is the name of a subroutine which you may choose
79
C to provide for defining a matrix of partial derivatives
82
C Quantities which may be altered by DDASSL are:
83
C T, Y(*), YPRIME(*), INFO(1), RTOL, ATOL,
84
C IDID, RWORK(*) AND IWORK(*)
88
C Subroutine DDASSL uses the backward differentiation formulas of
89
C orders one through five to solve a system of the above form for Y and
90
C YPRIME. Values for Y and YPRIME at the initial time must be given as
91
C input. These values must be consistent, (that is, if T,Y,YPRIME are
92
C the given initial values, they must satisfy G(T,Y,YPRIME) = 0.). The
93
C subroutine solves the system from T to TOUT. It is easy to continue
94
C the solution to get results at additional TOUT. This is the interval
95
C mode of operation. Intermediate results can also be obtained easily
96
C by using the intermediate-output capability.
98
C The following detailed description is divided into subsections:
99
C 1. Input required for the first call to DDASSL.
100
C 2. Output after any return from DDASSL.
101
C 3. What to do to continue the integration.
105
C -------- INPUT -- WHAT TO DO ON THE FIRST CALL TO DDASSL ------------
107
C The first call of the code is defined to be the start of each new
108
C problem. Read through the descriptions of all the following items,
109
C provide sufficient storage space for designated arrays, set
110
C appropriate variables for the initialization of the problem, and
111
C give information about how you want the problem to be solved.
114
C RES -- Provide a subroutine of the form
115
C SUBROUTINE RES(T,Y,YPRIME,DELTA,IRES,RPAR,IPAR)
116
C to define the system of differential/algebraic
117
C equations which is to be solved. For the given values
118
C of T,Y and YPRIME, the subroutine should
119
C return the residual of the defferential/algebraic
121
C DELTA = G(T,Y,YPRIME)
122
C (DELTA(*) is a vector of length NEQ which is
125
C Subroutine RES must not alter T,Y or YPRIME.
126
C You must declare the name RES in an external
127
C statement in your program that calls DDASSL.
128
C You must dimension Y,YPRIME and DELTA in RES.
130
C IRES is an integer flag which is always equal to
131
C zero on input. Subroutine RES should alter IRES
132
C only if it encounters an illegal value of Y or
133
C a stop condition. Set IRES = -1 if an input value
134
C is illegal, and DDASSL will try to solve the problem
135
C without getting IRES = -1. If IRES = -2, DDASSL
136
C will return control to the calling program
139
C RPAR and IPAR are real and integer parameter arrays which
140
C you can use for communication between your calling program
141
C and subroutine RES. They are not altered by DDASSL. If you
142
C do not need RPAR or IPAR, ignore these parameters by treat-
143
C ing them as dummy arguments. If you do choose to use them,
144
C dimension them in your calling program and in RES as arrays
145
C of appropriate length.
147
C NEQ -- Set it to the number of differential equations.
150
C T -- Set it to the initial point of the integration.
151
C T must be defined as a variable.
153
C Y(*) -- Set this vector to the initial values of the NEQ solution
154
C components at the initial point. You must dimension Y of
155
C length at least NEQ in your calling program.
157
C YPRIME(*) -- Set this vector to the initial values of the NEQ
158
C first derivatives of the solution components at the initial
159
C point. You must dimension YPRIME at least NEQ in your
160
C calling program. If you do not know initial values of some
161
C of the solution components, see the explanation of INFO(11).
163
C TOUT -- Set it to the first point at which a solution
164
C is desired. You can not take TOUT = T.
165
C integration either forward in T (TOUT .GT. T) or
166
C backward in T (TOUT .LT. T) is permitted.
168
C The code advances the solution from T to TOUT using
169
C step sizes which are automatically selected so as to
170
C achieve the desired accuracy. If you wish, the code will
171
C return with the solution and its derivative at
172
C intermediate steps (intermediate-output mode) so that
173
C you can monitor them, but you still must provide TOUT in
174
C accord with the basic aim of the code.
176
C The first step taken by the code is a critical one
177
C because it must reflect how fast the solution changes near
178
C the initial point. The code automatically selects an
179
C initial step size which is practically always suitable for
180
C the problem. By using the fact that the code will not step
181
C past TOUT in the first step, you could, if necessary,
182
C restrict the length of the initial step size.
184
C For some problems it may not be permissible to integrate
185
C past a point TSTOP because a discontinuity occurs there
186
C or the solution or its derivative is not defined beyond
187
C TSTOP. When you have declared a TSTOP point (SEE INFO(4)
188
C and RWORK(1)), you have told the code not to integrate
189
C past TSTOP. In this case any TOUT beyond TSTOP is invalid
192
C INFO(*) -- Use the INFO array to give the code more details about
193
C how you want your problem solved. This array should be
194
C dimensioned of length 15, though DDASSL uses only the first
195
C eleven entries. You must respond to all of the following
196
C items, which are arranged as questions. The simplest use
197
C of the code corresponds to answering all questions as yes,
198
C i.e. setting all entries of INFO to 0.
200
C INFO(1) - This parameter enables the code to initialize
201
C itself. You must set it to indicate the start of every
204
C **** Is this the first call for this problem ...
205
C Yes - Set INFO(1) = 0
206
C No - Not applicable here.
207
C See below for continuation calls. ****
209
C INFO(2) - How much accuracy you want of your solution
210
C is specified by the error tolerances RTOL and ATOL.
211
C The simplest use is to take them both to be scalars.
212
C To obtain more flexibility, they can both be vectors.
213
C The code must be told your choice.
215
C **** Are both error tolerances RTOL, ATOL scalars ...
216
C Yes - Set INFO(2) = 0
217
C and input scalars for both RTOL and ATOL
218
C No - Set INFO(2) = 1
219
C and input arrays for both RTOL and ATOL ****
221
C INFO(3) - The code integrates from T in the direction
222
C of TOUT by steps. If you wish, it will return the
223
C computed solution and derivative at the next
224
C intermediate step (the intermediate-output mode) or
225
C TOUT, whichever comes first. This is a good way to
226
C proceed if you want to see the behavior of the solution.
227
C If you must have solutions at a great many specific
228
C TOUT points, this code will compute them efficiently.
230
C **** Do you want the solution only at
231
C TOUT (and not at the next intermediate step) ...
232
C Yes - Set INFO(3) = 0
233
C No - Set INFO(3) = 1 ****
235
C INFO(4) - To handle solutions at a great many specific
236
C values TOUT efficiently, this code may integrate past
237
C TOUT and interpolate to obtain the result at TOUT.
238
C Sometimes it is not possible to integrate beyond some
239
C point TSTOP because the equation changes there or it is
240
C not defined past TSTOP. Then you must tell the code
243
C **** Can the integration be carried out without any
244
C restrictions on the independent variable T ...
245
C Yes - Set INFO(4)=0
247
C and define the stopping point TSTOP by
248
C setting RWORK(1)=TSTOP ****
250
C INFO(5) - To solve differential/algebraic problems it is
251
C necessary to use a matrix of partial derivatives of the
252
C system of differential equations. If you do not
253
C provide a subroutine to evaluate it analytically (see
254
C description of the item JAC in the call list), it will
255
C be approximated by numerical differencing in this code.
256
C although it is less trouble for you to have the code
257
C compute partial derivatives by numerical differencing,
258
C the solution will be more reliable if you provide the
259
C derivatives via JAC. Sometimes numerical differencing
260
C is cheaper than evaluating derivatives in JAC and
261
C sometimes it is not - this depends on your problem.
263
C **** Do you want the code to evaluate the partial
264
C derivatives automatically by numerical differences ...
265
C Yes - Set INFO(5)=0
267
C and provide subroutine JAC for evaluating the
268
C matrix of partial derivatives ****
270
C INFO(6) - DDASSL will perform much better if the matrix of
271
C partial derivatives, DG/DY + CJ*DG/DYPRIME,
272
C (here CJ is a scalar determined by DDASSL)
273
C is banded and the code is told this. In this
274
C case, the storage needed will be greatly reduced,
275
C numerical differencing will be performed much cheaper,
276
C and a number of important algorithms will execute much
277
C faster. The differential equation is said to have
278
C half-bandwidths ML (lower) and MU (upper) if equation i
279
C involves only unknowns Y(J) with
280
C I-ML .LE. J .LE. I+MU
281
C for all I=1,2,...,NEQ. Thus, ML and MU are the widths
282
C of the lower and upper parts of the band, respectively,
283
C with the main diagonal being excluded. If you do not
284
C indicate that the equation has a banded matrix of partial
285
C derivatives, the code works with a full matrix of NEQ**2
286
C elements (stored in the conventional way). Computations
287
C with banded matrices cost less time and storage than with
288
C full matrices if 2*ML+MU .LT. NEQ. If you tell the
289
C code that the matrix of partial derivatives has a banded
290
C structure and you want to provide subroutine JAC to
291
C compute the partial derivatives, then you must be careful
292
C to store the elements of the matrix in the special form
293
C indicated in the description of JAC.
295
C **** Do you want to solve the problem using a full
296
C (dense) matrix (and not a special banded
298
C Yes - Set INFO(6)=0
300
C and provide the lower (ML) and upper (MU)
301
C bandwidths by setting
306
C INFO(7) -- You can specify a maximum (absolute value of)
307
C stepsize, so that the code
308
C will avoid passing over very
311
C **** Do you want the code to decide
312
C on its own maximum stepsize?
313
C Yes - Set INFO(7)=0
315
C and define HMAX by setting
318
C INFO(8) -- Differential/algebraic problems
319
C may occaisionally suffer from
320
C severe scaling difficulties on the
321
C first step. If you know a great deal
322
C about the scaling of your problem, you can
323
C help to alleviate this problem by
324
C specifying an initial stepsize HO.
326
C **** Do you want the code to define
327
C its own initial stepsize?
328
C Yes - Set INFO(8)=0
330
C and define HO by setting
333
C INFO(9) -- If storage is a severe problem,
334
C you can save some locations by
335
C restricting the maximum order MAXORD.
336
C the default value is 5. for each
337
C order decrease below 5, the code
338
C requires NEQ fewer locations, however
339
C it is likely to be slower. In any
340
C case, you must have 1 .LE. MAXORD .LE. 5
341
C **** Do you want the maximum order to
343
C Yes - Set INFO(9)=0
345
C and define MAXORD by setting
346
C IWORK(3)=MAXORD ****
348
C INFO(10) --If you know that the solutions to your equations
349
C will always be nonnegative, it may help to set this
350
C parameter. However, it is probably best to
351
C try the code without using this option first,
352
C and only to use this option if that doesn't
354
C **** Do you want the code to solve the problem without
355
C invoking any special nonnegativity constraints?
356
C Yes - Set INFO(10)=0
357
C No - Set INFO(10)=1
359
C INFO(11) --DDASSL normally requires the initial T,
360
C Y, and YPRIME to be consistent. That is,
361
C you must have G(T,Y,YPRIME) = 0 at the initial
362
C time. If you do not know the initial
363
C derivative precisely, you can let DDASSL try
365
C **** Are the initialHE INITIAL T, Y, YPRIME consistent?
366
C Yes - Set INFO(11) = 0
367
C No - Set INFO(11) = 1,
368
C and set YPRIME to an initial approximation
369
C to YPRIME. (If you have no idea what
370
C YPRIME should be, set it to zero. Note
371
C that the initial Y should be such
372
C that there must exist a YPRIME so that
373
C G(T,Y,YPRIME) = 0.)
375
C RTOL, ATOL -- You must assign relative (RTOL) and absolute (ATOL
376
C error tolerances to tell the code how accurately you
377
C want the solution to be computed. They must be defined
378
C as variables because the code may change them. You
379
C have two choices --
380
C Both RTOL and ATOL are scalars. (INFO(2)=0)
381
C Both RTOL and ATOL are vectors. (INFO(2)=1)
382
C in either case all components must be non-negative.
384
C The tolerances are used by the code in a local error
385
C test at each step which requires roughly that
386
C ABS(LOCAL ERROR) .LE. RTOL*ABS(Y)+ATOL
387
C for each vector component.
388
C (More specifically, a root-mean-square norm is used to
389
C measure the size of vectors, and the error test uses the
390
C magnitude of the solution at the beginning of the step.)
392
C The true (global) error is the difference between the
393
C true solution of the initial value problem and the
394
C computed approximation. Practically all present day
395
C codes, including this one, control the local error at
396
C each step and do not even attempt to control the global
398
C Usually, but not always, the true accuracy of the
399
C computed Y is comparable to the error tolerances. This
400
C code will usually, but not always, deliver a more
401
C accurate solution if you reduce the tolerances and
402
C integrate again. By comparing two such solutions you
403
C can get a fairly reliable idea of the true error in the
404
C solution at the bigger tolerances.
406
C Setting ATOL=0. results in a pure relative error test on
407
C that component. Setting RTOL=0. results in a pure
408
C absolute error test on that component. A mixed test
409
C with non-zero RTOL and ATOL corresponds roughly to a
410
C relative error test when the solution component is much
411
C bigger than ATOL and to an absolute error test when the
412
C solution component is smaller than the threshhold ATOL.
414
C The code will not attempt to compute a solution at an
415
C accuracy unreasonable for the machine being used. It will
416
C advise you if you ask for too much accuracy and inform
417
C you as to the maximum accuracy it believes possible.
419
C RWORK(*) -- Dimension this real work array of length LRW in your
422
C LRW -- Set it to the declared length of the RWORK array.
424
C LRW .GE. 40+(MAXORD+4)*NEQ+NEQ**2
425
C for the full (dense) JACOBIAN case (when INFO(6)=0), or
426
C LRW .GE. 40+(MAXORD+4)*NEQ+(2*ML+MU+1)*NEQ
427
C for the banded user-defined JACOBIAN case
428
C (when INFO(5)=1 and INFO(6)=1), or
429
C LRW .GE. 40+(MAXORD+4)*NEQ+(2*ML+MU+1)*NEQ
430
C +2*(NEQ/(ML+MU+1)+1)
431
C for the banded finite-difference-generated JACOBIAN case
432
C (when INFO(5)=0 and INFO(6)=1)
434
C IWORK(*) -- Dimension this integer work array of length LIW in
435
C your calling program.
437
C LIW -- Set it to the declared length of the IWORK array.
438
C You must have LIW .GE. 20+NEQ
440
C RPAR, IPAR -- These are parameter arrays, of real and integer
441
C type, respectively. You can use them for communication
442
C between your program that calls DDASSL and the
443
C RES subroutine (and the JAC subroutine). They are not
444
C altered by DDASSL. If you do not need RPAR or IPAR,
445
C ignore these parameters by treating them as dummy
446
C arguments. If you do choose to use them, dimension
447
C them in your calling program and in RES (and in JAC)
448
C as arrays of appropriate length.
450
C JAC -- If you have set INFO(5)=0, you can ignore this parameter
451
C by treating it as a dummy argument. Otherwise, you must
452
C provide a subroutine of the form
453
C SUBROUTINE JAC(T,Y,YPRIME,PD,CJ,RPAR,IPAR)
454
C to define the matrix of partial derivatives
455
C PD=DG/DY+CJ*DG/DYPRIME
456
C CJ is a scalar which is input to JAC.
457
C For the given values of T,Y,YPRIME, the
458
C subroutine must evaluate the non-zero partial
459
C derivatives for each equation and each solution
460
C component, and store these values in the
461
C matrix PD. The elements of PD are set to zero
462
C before each call to JAC so only non-zero elements
463
C need to be defined.
465
C Subroutine JAC must not alter T,Y,(*),YPRIME(*), or CJ.
466
C You must declare the name JAC in an EXTERNAL statement in
467
C your program that calls DDASSL. You must dimension Y,
468
C YPRIME and PD in JAC.
470
C The way you must store the elements into the PD matrix
471
C depends on the structure of the matrix which you
472
C indicated by INFO(6).
473
C *** INFO(6)=0 -- Full (dense) matrix ***
474
C Give PD a first dimension of NEQ.
475
C When you evaluate the (non-zero) partial derivative
476
C of equation I with respect to variable J, you must
477
C store it in PD according to
478
C PD(I,J) = "DG(I)/DY(J)+CJ*DG(I)/DYPRIME(J)"
479
C *** INFO(6)=1 -- Banded JACOBIAN with ML lower and MU
480
C upper diagonal bands (refer to INFO(6) description
482
C Give PD a first dimension of 2*ML+MU+1.
483
C when you evaluate the (non-zero) partial derivative
484
C of equation I with respect to variable J, you must
485
C store it in PD according to
486
C IROW = I - J + ML + MU + 1
487
C PD(IROW,J) = "DG(I)/DY(J)+CJ*DG(I)/DYPRIME(J)"
489
C RPAR and IPAR are real and integer parameter arrays
490
C which you can use for communication between your calling
491
C program and your JACOBIAN subroutine JAC. They are not
492
C altered by DDASSL. If you do not need RPAR or IPAR,
493
C ignore these parameters by treating them as dummy
494
C arguments. If you do choose to use them, dimension
495
C them in your calling program and in JAC as arrays of
496
C appropriate length.
499
C OPTIONALLY REPLACEABLE NORM ROUTINE:
501
C DDASSL uses a weighted norm DDANRM to measure the size
502
C of vectors such as the estimated error in each step.
503
C A FUNCTION subprogram
504
C DOUBLE PRECISION FUNCTION DDANRM(NEQ,V,WT,RPAR,IPAR)
505
C DIMENSION V(NEQ),WT(NEQ)
506
C is used to define this norm. Here, V is the vector
507
C whose norm is to be computed, and WT is a vector of
508
C weights. A DDANRM routine has been included with DDASSL
509
C which computes the weighted root-mean-square norm
511
C DDANRM=SQRT((1/NEQ)*SUM(V(I)/WT(I))**2)
512
C this norm is suitable for most problems. In some
513
C special cases, it may be more convenient and/or
514
C efficient to define your own norm by writing a function
515
C subprogram to be called instead of DDANRM. This should,
516
C however, be attempted only after careful thought and
520
C -------- OUTPUT -- AFTER ANY RETURN FROM DDASSL ---------------------
522
C The principal aim of the code is to return a computed solution at
523
C TOUT, although it is also possible to obtain intermediate results
524
C along the way. To find out whether the code achieved its goal
525
C or if the integration process was interrupted before the task was
526
C completed, you must check the IDID parameter.
529
C T -- The solution was successfully advanced to the
532
C Y(*) -- Contains the computed solution approximation at T.
534
C YPRIME(*) -- Contains the computed derivative
535
C approximation at T.
537
C IDID -- Reports what the code did.
539
C *** Task completed ***
540
C Reported by positive values of IDID
542
C IDID = 1 -- A step was successfully taken in the
543
C intermediate-output mode. The code has not
546
C IDID = 2 -- The integration to TSTOP was successfully
547
C completed (T=TSTOP) by stepping exactly to TSTOP.
549
C IDID = 3 -- The integration to TOUT was successfully
550
C completed (T=TOUT) by stepping past TOUT.
551
C Y(*) is obtained by interpolation.
552
C YPRIME(*) is obtained by interpolation.
554
C *** Task interrupted ***
555
C Reported by negative values of IDID
557
C IDID = -1 -- A large amount of work has been expended.
560
C IDID = -2 -- The error tolerances are too stringent.
562
C IDID = -3 -- The local error test cannot be satisfied
563
C because you specified a zero component in ATOL
564
C and the corresponding computed solution
565
C component is zero. Thus, a pure relative error
566
C test is impossible for this component.
568
C IDID = -6 -- DDASSL had repeated error test
569
C failures on the last attempted step.
571
C IDID = -7 -- The corrector could not converge.
573
C IDID = -8 -- The matrix of partial derivatives
576
C IDID = -9 -- The corrector could not converge.
577
C there were repeated error test failures
580
C IDID =-10 -- The corrector could not converge
581
C because IRES was equal to minus one.
583
C IDID =-11 -- IRES equal to -2 was encountered
584
C and control is being returned to the
587
C IDID =-12 -- DDASSL failed to compute the initial
592
C IDID = -13,..,-32 -- Not applicable for this code
594
C *** Task terminated ***
595
C Reported by the value of IDID=-33
597
C IDID = -33 -- The code has encountered trouble from which
598
C it cannot recover. A message is printed
599
C explaining the trouble and control is returned
600
C to the calling program. For example, this occurs
601
C when invalid input is detected.
603
C RTOL, ATOL -- These quantities remain unchanged except when
604
C IDID = -2. In this case, the error tolerances have been
605
C increased by the code to values which are estimated to
606
C be appropriate for continuing the integration. However,
607
C the reported solution at T was obtained using the input
608
C values of RTOL and ATOL.
610
C RWORK, IWORK -- Contain information which is usually of no
611
C interest to the user but necessary for subsequent calls.
612
C However, you may find use for
614
C RWORK(3)--Which contains the step size H to be
615
C attempted on the next step.
617
C RWORK(4)--Which contains the current value of the
618
C independent variable, i.e., the farthest point
619
C integration has reached. This will be different
620
C from T only when interpolation has been
621
C performed (IDID=3).
623
C RWORK(7)--Which contains the stepsize used
624
C on the last successful step.
626
C IWORK(7)--Which contains the order of the method to
627
C be attempted on the next step.
629
C IWORK(8)--Which contains the order of the method used
632
C IWORK(11)--Which contains the number of steps taken so
635
C IWORK(12)--Which contains the number of calls to RES
638
C IWORK(13)--Which contains the number of evaluations of
639
C the matrix of partial derivatives needed so
642
C IWORK(14)--Which contains the total number
643
C of error test failures so far.
645
C IWORK(15)--Which contains the total number
646
C of convergence test failures so far.
647
C (includes singular iteration matrix
651
C -------- INPUT -- WHAT TO DO TO CONTINUE THE INTEGRATION ------------
652
C (CALLS AFTER THE FIRST)
654
C This code is organized so that subsequent calls to continue the
655
C integration involve little (if any) additional effort on your
656
C part. You must monitor the IDID parameter in order to determine
659
C Recalling that the principal task of the code is to integrate
660
C from T to TOUT (the interval mode), usually all you will need
661
C to do is specify a new TOUT upon reaching the current TOUT.
663
C Do not alter any quantity not specifically permitted below,
664
C in particular do not alter NEQ,T,Y(*),YPRIME(*),RWORK(*),IWORK(*)
665
C or the differential equation in subroutine RES. Any such
666
C alteration constitutes a new problem and must be treated as such,
667
C i.e., you must start afresh.
669
C You cannot change from vector to scalar error control or vice
670
C versa (INFO(2)), but you can change the size of the entries of
671
C RTOL, ATOL. Increasing a tolerance makes the equation easier
672
C to integrate. Decreasing a tolerance will make the equation
673
C harder to integrate and should generally be avoided.
675
C You can switch from the intermediate-output mode to the
676
C interval mode (INFO(3)) or vice versa at any time.
678
C If it has been necessary to prevent the integration from going
679
C past a point TSTOP (INFO(4), RWORK(1)), keep in mind that the
680
C code will not integrate to any TOUT beyond the currently
681
C specified TSTOP. Once TSTOP has been reached you must change
682
C the value of TSTOP or set INFO(4)=0. You may change INFO(4)
683
C or TSTOP at any time but you must supply the value of TSTOP in
684
C RWORK(1) whenever you set INFO(4)=1.
686
C Do not change INFO(5), INFO(6), IWORK(1), or IWORK(2)
687
C unless you are going to restart the code.
689
C *** Following a completed task ***
691
C IDID = 1, call the code again to continue the integration
692
C another step in the direction of TOUT.
694
C IDID = 2 or 3, define a new TOUT and call the code again.
695
C TOUT must be different from T. You cannot change
696
C the direction of integration without restarting.
698
C *** Following an interrupted task ***
699
C To show the code that you realize the task was
700
C interrupted and that you want to continue, you
701
C must take appropriate action and set INFO(1) = 1
703
C IDID = -1, The code has taken about 500 steps.
704
C If you want to continue, set INFO(1) = 1 and
705
C call the code again. An additional 500 steps
708
C IDID = -2, The error tolerances RTOL, ATOL have been
709
C increased to values the code estimates appropriate
710
C for continuing. You may want to change them
711
C yourself. If you are sure you want to continue
712
C with relaxed error tolerances, set INFO(1)=1 and
713
C call the code again.
715
C IDID = -3, A solution component is zero and you set the
716
C corresponding component of ATOL to zero. If you
717
C are sure you want to continue, you must first
718
C alter the error criterion to use positive values
719
C for those components of ATOL corresponding to zero
720
C solution components, then set INFO(1)=1 and call
723
C IDID = -4,-5 --- Cannot occur with this code.
725
C IDID = -6, Repeated error test failures occurred on the
726
C last attempted step in DDASSL. A singularity in the
727
C solution may be present. If you are absolutely
728
C certain you want to continue, you should restart
729
C the integration. (Provide initial values of Y and
730
C YPRIME which are consistent)
732
C IDID = -7, Repeated convergence test failures occurred
733
C on the last attempted step in DDASSL. An inaccurate
734
C or ill-conditioned JACOBIAN may be the problem. If
735
C you are absolutely certain you want to continue, you
736
C should restart the integration.
738
C IDID = -8, The matrix of partial derivatives is singular.
739
C Some of your equations may be redundant.
740
C DDASSL cannot solve the problem as stated.
741
C It is possible that the redundant equations
742
C could be removed, and then DDASSL could
743
C solve the problem. It is also possible
744
C that a solution to your problem either
745
C does not exist or is not unique.
747
C IDID = -9, DDASSL had multiple convergence test
748
C failures, preceeded by multiple error
749
C test failures, on the last attempted step.
750
C It is possible that your problem
751
C is ill-posed, and cannot be solved
752
C using this code. Or, there may be a
753
C discontinuity or a singularity in the
754
C solution. If you are absolutely certain
755
C you want to continue, you should restart
758
C IDID =-10, DDASSL had multiple convergence test failures
759
C because IRES was equal to minus one.
760
C If you are absolutely certain you want
761
C to continue, you should restart the
764
C IDID =-11, IRES=-2 was encountered, and control is being
765
C returned to the calling program.
767
C IDID =-12, DDASSL failed to compute the initial YPRIME.
768
C This could happen because the initial
769
C approximation to YPRIME was not very good, or
770
C if a YPRIME consistent with the initial Y
771
C does not exist. The problem could also be caused
772
C by an inaccurate or singular iteration matrix.
774
C IDID = -13,..,-32 --- Cannot occur with this code.
777
C *** Following a terminated task ***
779
C If IDID= -33, you cannot continue the solution of this problem.
780
C An attempt to do so will result in your
781
C run being terminated.
784
C -------- ERROR MESSAGES ---------------------------------------------
786
C The SLATEC error print routine XERMSG is called in the event of
787
C unsuccessful completion of a task. Most of these are treated as
788
C "recoverable errors", which means that (unless the user has directed
789
C otherwise) control will be returned to the calling program for
790
C possible action after the message has been printed.
792
C In the event of a negative value of IDID other than -33, an appro-
793
C priate message is printed and the "error number" printed by XERMSG
794
C is the value of IDID. There are quite a number of illegal input
795
C errors that can lead to a returned value IDID=-33. The conditions
796
C and their printed "error numbers" are as follows:
798
C Error number Condition
800
C 1 Some element of INFO vector is not zero or one.
802
C 3 MAXORD not in range.
803
C 4 LRW is less than the required length for RWORK.
804
C 5 LIW is less than the required length for IWORK.
805
C 6 Some element of RTOL is .lt. 0
806
C 7 Some element of ATOL is .lt. 0
807
C 8 All elements of RTOL and ATOL are zero.
808
C 9 INFO(4)=1 and TSTOP is behind TOUT.
810
C 11 TOUT is behind T.
811
C 12 INFO(8)=1 and H0=0.0
812
C 13 Some element of WT is .le. 0.0
813
C 14 TOUT is too close to T to start integration.
814
C 15 INFO(4)=1 and TSTOP is behind T.
815
C 16 --( Not used in this version )--
816
C 17 ML illegal. Either .lt. 0 or .gt. NEQ
817
C 18 MU illegal. Either .lt. 0 or .gt. NEQ
820
C If DDASSL is called again without any action taken to remove the
821
C cause of an unsuccessful return, XERMSG will be called with a fatal
822
C error flag, which will cause unconditional termination of the
823
C program. There are two such fatal errors:
825
C Error number -998: The last step was terminated with a negative
826
C value of IDID other than -33, and no appropriate action was
829
C Error number -999: The previous call was terminated because of
830
C illegal input (IDID=-33) and there is illegal input in the
831
C present call, as well. (Suspect infinite loop.)
833
C ---------------------------------------------------------------------
835
C***REFERENCES A DESCRIPTION OF DASSL: A DIFFERENTIAL/ALGEBRAIC
836
C SYSTEM SOLVER, L. R. PETZOLD, SAND82-8637,
837
C SANDIA NATIONAL LABORATORIES, SEPTEMBER 1982.
838
C***ROUTINES CALLED D1MACH, DDAINI, DDANRM, DDASTP, DDATRP, DDAWTS,
840
C***REVISION HISTORY (YYMMDD)
841
C 830315 DATE WRITTEN
842
C 880387 Code changes made. All common statements have been
843
C replaced by a DATA statement, which defines pointers into
844
C RWORK, and PARAMETER statements which define pointers
845
C into IWORK. As well the documentation has gone through
846
C grammatical changes.
847
C 881005 The prologue has been changed to mixed case.
848
C The subordinate routines had revision dates changed to
849
C this date, although the documentation for these routines
850
C is all upper case. No code changes.
851
C 890511 Code changes made. The DATA statement in the declaration
852
C section of DDASSL was replaced with a PARAMETER
853
C statement. Also the statement S = 100.D0 was removed
854
C from the top of the Newton iteration in DDASTP.
855
C The subordinate routines had revision dates changed to
857
C 890517 The revision date syntax was replaced with the revision
858
C history syntax. Also the "DECK" comment was added to
859
C the top of all subroutines. These changes are consistent
860
C with new SLATEC guidelines.
861
C The subordinate routines had revision dates changed to
862
C this date. No code changes.
863
C 891013 Code changes made.
864
C Removed all occurrances of FLOAT or DBLE. All operations
865
C are now performed with "mixed-mode" arithmetic.
866
C Also, specific function names were replaced with generic
867
C function names to be consistent with new SLATEC guidelines.
869
C Replaced DSQRT with SQRT everywhere.
870
C Replaced DABS with ABS everywhere.
871
C Replaced DMIN1 with MIN everywhere.
872
C Replaced MIN0 with MIN everywhere.
873
C Replaced DMAX1 with MAX everywhere.
874
C Replaced MAX0 with MAX everywhere.
875
C Replaced DSIGN with SIGN everywhere.
876
C Also replaced REVISION DATE with REVISION HISTORY in all
877
C subordinate routines.
878
C 901004 Miscellaneous changes to prologue to complete conversion
879
C to SLATEC 4.0 format. No code changes. (F.N.Fritsch)
880
C 901009 Corrected GAMS classification code and converted subsidiary
881
C routines to 4.0 format. No code changes. (F.N.Fritsch)
882
C 901010 Converted XERRWV calls to XERMSG calls. (R.Clemens,AFWL)
883
C 901019 Code changes made.
884
C Merged SLATEC 4.0 changes with previous changes made
885
C by C. Ulrich. Below is a history of the changes made by
886
C C. Ulrich. (Changes in subsidiary routines are implied
888
C 891228 Bug was found and repaired inside the DDASSL
889
C and DDAINI routines. DDAINI was incorrectly
890
C returning the initial T with Y and YPRIME
891
C computed at T+H. The routine now returns T+H
892
C rather than the initial T.
893
C Cosmetic changes made to DDASTP.
894
C 900904 Three modifications were made to fix a bug (inside
895
C DDASSL) re interpolation for continuation calls and
896
C cases where TN is very close to TSTOP:
898
C 1) In testing for whether H is too large, just
899
C compare H to (TSTOP - TN), rather than
900
C (TSTOP - TN) * (1-4*UROUND), and set H to
901
C TSTOP - TN. This will force DDASTP to step
902
C exactly to TSTOP under certain situations
903
C (i.e. when H returned from DDASTP would otherwise
904
C take TN beyond TSTOP).
906
C 2) Inside the DDASTP loop, interpolate exactly to
907
C TSTOP if TN is very close to TSTOP (rather than
908
C interpolating to within roundoff of TSTOP).
910
C 3) Modified IDID description for IDID = 2 to say that
911
C the solution is returned by stepping exactly to
912
C TSTOP, rather than TOUT. (In some cases the
913
C solution is actually obtained by extrapolating
914
C over a distance near unit roundoff to TSTOP,
915
C but this small distance is deemed acceptable in
916
C these circumstances.)
917
C 901026 Added explicit declarations for all variables and minor
918
C cosmetic changes to prologue, removed unreferenced labels,
919
C and improved XERMSG calls. (FNF)
920
C 901030 Added ERROR MESSAGES section and reworked other sections to
921
C be of more uniform format. (FNF)
922
C 910624 Fixed minor bug related to HMAX (five lines ending in
923
C statement 526 in DDASSL). (LRP)
925
C***END PROLOGUE DDASSL
931
INTEGER NEQ, INFO(15), IDID, LRW, IWORK(*), LIW, IPAR(*)
933
* T, Y(*), YPRIME(*), TOUT, RTOL(*), ATOL(*), RWORK(*),
939
EXTERNAL D1MACH, DDAINI, DDANRM, DDASTP, DDATRP, DDAWTS, XERMSG
940
DOUBLE PRECISION D1MACH, DDANRM
942
C Declare local variables.
944
INTEGER I, ITEMP, LALPHA, LBETA, LCJ, LCJOLD, LCTF, LDELTA,
945
* LENIW, LENPD, LENRW, LE, LETF, LGAMMA, LH, LHMAX, LHOLD, LIPVT,
946
* LJCALC, LK, LKOLD, LIWM, LML, LMTYPE, LMU, LMXORD, LNJE, LNPD,
947
* LNRE, LNS, LNST, LNSTL, LPD, LPHASE, LPHI, LPSI, LROUND, LS,
948
* LSIGMA, LTN, LTSTOP, LWM, LWT, MBAND, MSAVE, MXORD, NPD, NTEMP,
951
* ATOLI, H, HMAX, HMIN, HO, R, RH, RTOLI, TDIST, TN, TNEXT,
952
* TSTOP, UROUND, YPNORM
954
C Auxiliary variables for conversion of values to be included in
956
CHARACTER*8 XERN1, XERN2
957
CHARACTER*16 XERN3, XERN4
959
C SET POINTERS INTO IWORK
960
PARAMETER (LML=1, LMU=2, LMXORD=3, LMTYPE=4, LNST=11,
961
* LNRE=12, LNJE=13, LETF=14, LCTF=15, LNPD=16,
962
* LIPVT=21, LJCALC=5, LPHASE=6, LK=7, LKOLD=8,
963
* LNS=9, LNSTL=10, LIWM=1)
965
C SET RELATIVE OFFSET INTO RWORK
968
C SET POINTERS INTO RWORK
969
PARAMETER (LTSTOP=1, LHMAX=2, LH=3, LTN=4,
970
* LCJ=5, LCJOLD=6, LHOLD=7, LS=8, LROUND=9,
971
* LALPHA=11, LBETA=17, LGAMMA=23,
972
* LPSI=29, LSIGMA=35, LDELTA=41)
974
C***FIRST EXECUTABLE STATEMENT DDASSL
975
IF(INFO(1).NE.0)GO TO 100
977
C-----------------------------------------------------------------------
978
C THIS BLOCK IS EXECUTED FOR THE INITIAL CALL ONLY.
979
C IT CONTAINS CHECKING OF INPUTS AND INITIALIZATIONS.
980
C-----------------------------------------------------------------------
982
C FIRST CHECK INFO ARRAY TO MAKE SURE ALL ELEMENTS OF INFO
983
C ARE EITHER ZERO OR ONE.
985
IF(INFO(I).NE.0.AND.INFO(I).NE.1)GO TO 701
988
IF(NEQ.LE.0)GO TO 702
990
C CHECK AND COMPUTE MAXIMUM ORDER
992
IF(INFO(9).EQ.0)GO TO 20
994
IF(MXORD.LT.1.OR.MXORD.GT.5)GO TO 703
995
20 IWORK(LMXORD)=MXORD
997
C COMPUTE MTYPE,LENPD,LENRW.CHECK ML AND MU.
998
IF(INFO(6).NE.0)GO TO 40
1000
LENRW=40+(IWORK(LMXORD)+4)*NEQ+LENPD
1001
IF(INFO(5).NE.0)GO TO 30
1006
40 IF(IWORK(LML).LT.0.OR.IWORK(LML).GE.NEQ)GO TO 717
1007
IF(IWORK(LMU).LT.0.OR.IWORK(LMU).GE.NEQ)GO TO 718
1008
LENPD=(2*IWORK(LML)+IWORK(LMU)+1)*NEQ
1009
IF(INFO(5).NE.0)GO TO 50
1011
MBAND=IWORK(LML)+IWORK(LMU)+1
1013
LENRW=40+(IWORK(LMXORD)+4)*NEQ+LENPD+2*MSAVE
1016
LENRW=40+(IWORK(LMXORD)+4)*NEQ+LENPD
1018
C CHECK LENGTHS OF RWORK AND IWORK
1021
IF(LRW.LT.LENRW)GO TO 704
1022
IF(LIW.LT.LENIW)GO TO 705
1024
C CHECK TO SEE THAT TOUT IS DIFFERENT FROM T
1025
IF(TOUT .EQ. T)GO TO 719
1028
IF(INFO(7).EQ.0)GO TO 70
1030
IF(HMAX.LE.0.0D0)GO TO 710
1033
C INITIALIZE COUNTERS
1042
C-----------------------------------------------------------------------
1043
C THIS BLOCK IS FOR CONTINUATION CALLS
1044
C ONLY. HERE WE CHECK INFO(1),AND IF THE
1045
C LAST STEP WAS INTERRUPTED WE CHECK WHETHER
1046
C APPROPRIATE ACTION WAS TAKEN.
1047
C-----------------------------------------------------------------------
1050
IF(INFO(1).EQ.1)GO TO 110
1051
IF(INFO(1).NE.-1)GO TO 701
1053
C IF WE ARE HERE, THE LAST STEP WAS INTERRUPTED
1054
C BY AN ERROR CONDITION FROM DDASTP,AND
1055
C APPROPRIATE ACTION WAS NOT TAKEN. THIS
1057
WRITE (XERN1, '(I8)') IDID
1058
CALL XERMSG ('SLATEC', 'DDASSL',
1059
* 'THE LAST STEP TERMINATED WITH A NEGATIVE VALUE OF IDID = ' //
1060
* XERN1 // ' AND NO APPROPRIATE ACTION WAS TAKEN. ' //
1061
* 'RUN TERMINATED', -998, 2)
1064
IWORK(LNSTL)=IWORK(LNST)
1066
C-----------------------------------------------------------------------
1067
C THIS BLOCK IS EXECUTED ON ALL CALLS.
1068
C THE ERROR TOLERANCE PARAMETERS ARE
1069
C CHECKED, AND THE WORK ARRAY POINTERS
1071
C-----------------------------------------------------------------------
1079
IF(INFO(2).EQ.1)RTOLI=RTOL(I)
1080
IF(INFO(2).EQ.1)ATOLI=ATOL(I)
1081
IF(RTOLI.GT.0.0D0.OR.ATOLI.GT.0.0D0)NZFLG=1
1082
IF(RTOLI.LT.0.0D0)GO TO 706
1083
IF(ATOLI.LT.0.0D0)GO TO 707
1085
IF(NZFLG.EQ.0)GO TO 708
1087
C SET UP RWORK STORAGE.IWORK STORAGE IS FIXED
1088
C IN DATA STATEMENT.
1092
LPD=LPHI+(IWORK(LMXORD)+1)*NEQ
1094
NTEMP=NPD+IWORK(LNPD)
1095
IF(INFO(1).EQ.1)GO TO 400
1097
C-----------------------------------------------------------------------
1098
C THIS BLOCK IS EXECUTED ON THE INITIAL CALL
1099
C ONLY. SET THE INITIAL STEP SIZE, AND
1100
C THE ERROR WEIGHT VECTOR, AND PHI.
1101
C COMPUTE INITIAL YPRIME, IF NECESSARY.
1102
C-----------------------------------------------------------------------
1107
C SET ERROR WEIGHT VECTOR WT
1108
CALL DDAWTS(NEQ,INFO(2),RTOL,ATOL,Y,RWORK(LWT),RPAR,IPAR)
1110
IF(RWORK(LWT+I-1).LE.0.0D0) GO TO 713
1113
C COMPUTE UNIT ROUNDOFF AND HMIN
1115
RWORK(LROUND) = UROUND
1116
HMIN = 4.0D0*UROUND*MAX(ABS(T),ABS(TOUT))
1118
C CHECK INITIAL INTERVAL TO SEE THAT IT IS LONG ENOUGH
1119
TDIST = ABS(TOUT - T)
1120
IF(TDIST .LT. HMIN) GO TO 714
1122
C CHECK HO, IF THIS WAS INPUT
1123
IF (INFO(8) .EQ. 0) GO TO 310
1125
IF ((TOUT - T)*HO .LT. 0.0D0) GO TO 711
1126
IF (HO .EQ. 0.0D0) GO TO 712
1130
C COMPUTE INITIAL STEPSIZE, TO BE USED BY EITHER
1131
C DDASTP OR DDAINI, DEPENDING ON INFO(11)
1133
YPNORM = DDANRM(NEQ,YPRIME,RWORK(LWT),RPAR,IPAR)
1134
IF (YPNORM .GT. 0.5D0/HO) HO = 0.5D0/YPNORM
1135
HO = SIGN(HO,TOUT-T)
1136
C ADJUST HO IF NECESSARY TO MEET HMAX BOUND
1137
320 IF (INFO(7) .EQ. 0) GO TO 330
1138
RH = ABS(HO)/RWORK(LHMAX)
1139
IF (RH .GT. 1.0D0) HO = HO/RH
1140
C COMPUTE TSTOP, IF APPLICABLE
1141
330 IF (INFO(4) .EQ. 0) GO TO 340
1142
TSTOP = RWORK(LTSTOP)
1143
IF ((TSTOP - T)*HO .LT. 0.0D0) GO TO 715
1144
IF ((T + HO - TSTOP)*HO .GT. 0.0D0) HO = TSTOP - T
1145
IF ((TSTOP - TOUT)*HO .LT. 0.0D0) GO TO 709
1147
C COMPUTE INITIAL DERIVATIVE, UPDATING TN AND Y, IF APPLICABLE
1148
340 IF (INFO(11) .EQ. 0) GO TO 350
1149
CALL DDAINI(TN,Y,YPRIME,NEQ,
1150
* RES,JAC,HO,RWORK(LWT),IDID,RPAR,IPAR,
1151
* RWORK(LPHI),RWORK(LDELTA),RWORK(LE),
1152
* RWORK(LWM),IWORK(LIWM),HMIN,RWORK(LROUND),
1154
IF (IDID .LT. 0) GO TO 390
1156
C LOAD H WITH HO. STORE H IN RWORK(LH)
1160
C LOAD Y AND H*YPRIME INTO PHI(*,1) AND PHI(*,2)
1163
RWORK(LPHI + I - 1) = Y(I)
1164
370 RWORK(ITEMP + I - 1) = H*YPRIME(I)
1168
C-------------------------------------------------------
1169
C THIS BLOCK IS FOR CONTINUATION CALLS ONLY. ITS
1170
C PURPOSE IS TO CHECK STOP CONDITIONS BEFORE
1172
C ADJUST H IF NECESSARY TO MEET HMAX BOUND
1173
C-------------------------------------------------------
1176
UROUND=RWORK(LROUND)
1180
IF(INFO(7) .EQ. 0) GO TO 410
1181
RH = ABS(H)/RWORK(LHMAX)
1182
IF(RH .GT. 1.0D0) H = H/RH
1184
IF(T .EQ. TOUT) GO TO 719
1185
IF((T - TOUT)*H .GT. 0.0D0) GO TO 711
1186
IF(INFO(4) .EQ. 1) GO TO 430
1187
IF(INFO(3) .EQ. 1) GO TO 420
1188
IF((TN-TOUT)*H.LT.0.0D0)GO TO 490
1189
CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD),
1190
* RWORK(LPHI),RWORK(LPSI))
1195
420 IF((TN-T)*H .LE. 0.0D0) GO TO 490
1196
IF((TN - TOUT)*H .GT. 0.0D0) GO TO 425
1197
CALL DDATRP(TN,TN,Y,YPRIME,NEQ,IWORK(LKOLD),
1198
* RWORK(LPHI),RWORK(LPSI))
1204
CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD),
1205
* RWORK(LPHI),RWORK(LPSI))
1210
430 IF(INFO(3) .EQ. 1) GO TO 440
1212
IF((TN-TSTOP)*H.GT.0.0D0) GO TO 715
1213
IF((TSTOP-TOUT)*H.LT.0.0D0)GO TO 709
1214
IF((TN-TOUT)*H.LT.0.0D0)GO TO 450
1215
CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD),
1216
* RWORK(LPHI),RWORK(LPSI))
1221
440 TSTOP = RWORK(LTSTOP)
1222
IF((TN-TSTOP)*H .GT. 0.0D0) GO TO 715
1223
IF((TSTOP-TOUT)*H .LT. 0.0D0) GO TO 709
1224
IF((TN-T)*H .LE. 0.0D0) GO TO 450
1225
IF((TN - TOUT)*H .GT. 0.0D0) GO TO 445
1226
CALL DDATRP(TN,TN,Y,YPRIME,NEQ,IWORK(LKOLD),
1227
* RWORK(LPHI),RWORK(LPSI))
1233
CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD),
1234
* RWORK(LPHI),RWORK(LPSI))
1240
C CHECK WHETHER WE ARE WITHIN ROUNDOFF OF TSTOP
1241
IF(ABS(TN-TSTOP).GT.100.0D0*UROUND*
1242
* (ABS(TN)+ABS(H)))GO TO 460
1243
CALL DDATRP(TN,TSTOP,Y,YPRIME,NEQ,IWORK(LKOLD),
1244
* RWORK(LPHI),RWORK(LPSI))
1250
IF((TNEXT-TSTOP)*H.LE.0.0D0)GO TO 490
1254
490 IF (DONE) GO TO 580
1256
C-------------------------------------------------------
1257
C THE NEXT BLOCK CONTAINS THE CALL TO THE
1258
C ONE-STEP INTEGRATOR DDASTP.
1259
C THIS IS A LOOPING POINT FOR THE INTEGRATION STEPS.
1260
C CHECK FOR TOO MANY STEPS.
1262
C CHECK FOR TOO MUCH ACCURACY REQUESTED.
1263
C COMPUTE MINIMUM STEPSIZE.
1264
C-------------------------------------------------------
1267
C CHECK FOR FAILURE TO COMPUTE INITIAL YPRIME
1268
IF (IDID .EQ. -12) GO TO 527
1270
C CHECK FOR TOO MANY STEPS
1271
IF((IWORK(LNST)-IWORK(LNSTL)).LT.500)
1277
510 CALL DDAWTS(NEQ,INFO(2),RTOL,ATOL,RWORK(LPHI),
1278
* RWORK(LWT),RPAR,IPAR)
1280
IF(RWORK(I+LWT-1).GT.0.0D0)GO TO 520
1285
C TEST FOR TOO MUCH ACCURACY REQUESTED.
1286
R=DDANRM(NEQ,RWORK(LPHI),RWORK(LWT),RPAR,IPAR)*
1288
IF(R.LE.1.0D0)GO TO 525
1289
C MULTIPLY RTOL AND ATOL BY R AND RETURN
1290
IF(INFO(2).EQ.1)GO TO 523
1297
524 ATOL(I)=R*ATOL(I)
1302
C COMPUTE MINIMUM STEPSIZE
1303
HMIN=4.0D0*UROUND*MAX(ABS(TN),ABS(TOUT))
1306
IF (INFO(7) .EQ. 0) GO TO 526
1307
RH = ABS(H)/RWORK(LHMAX)
1308
IF (RH .GT. 1.0D0) H = H/RH
1311
CALL DDASTP(TN,Y,YPRIME,NEQ,
1312
* RES,JAC,H,RWORK(LWT),INFO(1),IDID,RPAR,IPAR,
1313
* RWORK(LPHI),RWORK(LDELTA),RWORK(LE),
1314
* RWORK(LWM),IWORK(LIWM),
1315
* RWORK(LALPHA),RWORK(LBETA),RWORK(LGAMMA),
1316
* RWORK(LPSI),RWORK(LSIGMA),
1317
* RWORK(LCJ),RWORK(LCJOLD),RWORK(LHOLD),
1318
* RWORK(LS),HMIN,RWORK(LROUND),
1319
* IWORK(LPHASE),IWORK(LJCALC),IWORK(LK),
1320
* IWORK(LKOLD),IWORK(LNS),INFO(10),NTEMP)
1321
527 IF(IDID.LT.0)GO TO 600
1323
C--------------------------------------------------------
1324
C THIS BLOCK HANDLES THE CASE OF A SUCCESSFUL RETURN
1325
C FROM DDASTP (IDID=1). TEST FOR STOP CONDITIONS.
1326
C--------------------------------------------------------
1328
IF(INFO(4).NE.0)GO TO 540
1329
IF(INFO(3).NE.0)GO TO 530
1330
IF((TN-TOUT)*H.LT.0.0D0)GO TO 500
1331
CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,
1332
* IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
1336
530 IF((TN-TOUT)*H.GE.0.0D0)GO TO 535
1340
535 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,
1341
* IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
1345
540 IF(INFO(3).NE.0)GO TO 550
1346
IF((TN-TOUT)*H.LT.0.0D0)GO TO 542
1347
CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,
1348
* IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
1352
542 IF(ABS(TN-TSTOP).LE.100.0D0*UROUND*
1353
* (ABS(TN)+ABS(H)))GO TO 545
1355
IF((TNEXT-TSTOP)*H.LE.0.0D0)GO TO 500
1358
545 CALL DDATRP(TN,TSTOP,Y,YPRIME,NEQ,
1359
* IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
1363
550 IF((TN-TOUT)*H.GE.0.0D0)GO TO 555
1364
IF(ABS(TN-TSTOP).LE.100.0D0*UROUND*(ABS(TN)+ABS(H)))GO TO 552
1368
552 CALL DDATRP(TN,TSTOP,Y,YPRIME,NEQ,
1369
* IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
1373
555 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,
1374
* IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
1379
C--------------------------------------------------------
1380
C ALL SUCCESSFUL RETURNS FROM DDASSL ARE MADE FROM
1382
C--------------------------------------------------------
1389
C-----------------------------------------------------------------------
1390
C THIS BLOCK HANDLES ALL UNSUCCESSFUL
1391
C RETURNS OTHER THAN FOR ILLEGAL INPUT.
1392
C-----------------------------------------------------------------------
1396
GO TO (610,620,630,690,690,640,650,660,670,675,
1399
C THE MAXIMUM NUMBER OF STEPS WAS TAKEN BEFORE
1401
610 WRITE (XERN3, '(1P,D15.6)') TN
1402
CALL XERMSG ('SLATEC', 'DDASSL',
1403
* 'AT CURRENT T = ' // XERN3 // ' 500 STEPS TAKEN ON THIS ' //
1404
* 'CALL BEFORE REACHING TOUT', IDID, 1)
1407
C TOO MUCH ACCURACY FOR MACHINE PRECISION
1408
620 WRITE (XERN3, '(1P,D15.6)') TN
1409
CALL XERMSG ('SLATEC', 'DDASSL',
1410
* 'AT T = ' // XERN3 // ' TOO MUCH ACCURACY REQUESTED FOR ' //
1411
* 'PRECISION OF MACHINE. RTOL AND ATOL WERE INCREASED TO ' //
1412
* 'APPROPRIATE VALUES', IDID, 1)
1415
C WT(I) .LE. 0.0 FOR SOME I (NOT AT START OF PROBLEM)
1416
630 WRITE (XERN3, '(1P,D15.6)') TN
1417
CALL XERMSG ('SLATEC', 'DDASSL',
1418
* 'AT T = ' // XERN3 // ' SOME ELEMENT OF WT HAS BECOME .LE. ' //
1422
C ERROR TEST FAILED REPEATEDLY OR WITH H=HMIN
1423
640 WRITE (XERN3, '(1P,D15.6)') TN
1424
WRITE (XERN4, '(1P,D15.6)') H
1425
CALL XERMSG ('SLATEC', 'DDASSL',
1426
* 'AT T = ' // XERN3 // ' AND STEPSIZE H = ' // XERN4 //
1427
* ' THE ERROR TEST FAILED REPEATEDLY OR WITH ABS(H)=HMIN',
1431
C CORRECTOR CONVERGENCE FAILED REPEATEDLY OR WITH H=HMIN
1432
650 WRITE (XERN3, '(1P,D15.6)') TN
1433
WRITE (XERN4, '(1P,D15.6)') H
1434
CALL XERMSG ('SLATEC', 'DDASSL',
1435
* 'AT T = ' // XERN3 // ' AND STEPSIZE H = ' // XERN4 //
1436
* ' THE CORRECTOR FAILED TO CONVERGE REPEATEDLY OR WITH ' //
1437
* 'ABS(H)=HMIN', IDID, 1)
1440
C THE ITERATION MATRIX IS SINGULAR
1441
660 WRITE (XERN3, '(1P,D15.6)') TN
1442
WRITE (XERN4, '(1P,D15.6)') H
1443
CALL XERMSG ('SLATEC', 'DDASSL',
1444
* 'AT T = ' // XERN3 // ' AND STEPSIZE H = ' // XERN4 //
1445
* ' THE ITERATION MATRIX IS SINGULAR', IDID, 1)
1448
C CORRECTOR FAILURE PRECEEDED BY ERROR TEST FAILURES.
1449
670 WRITE (XERN3, '(1P,D15.6)') TN
1450
WRITE (XERN4, '(1P,D15.6)') H
1451
CALL XERMSG ('SLATEC', 'DDASSL',
1452
* 'AT T = ' // XERN3 // ' AND STEPSIZE H = ' // XERN4 //
1453
* ' THE CORRECTOR COULD NOT CONVERGE. ALSO, THE ERROR TEST ' //
1454
* 'FAILED REPEATEDLY.', IDID, 1)
1457
C CORRECTOR FAILURE BECAUSE IRES = -1
1458
675 WRITE (XERN3, '(1P,D15.6)') TN
1459
WRITE (XERN4, '(1P,D15.6)') H
1460
CALL XERMSG ('SLATEC', 'DDASSL',
1461
* 'AT T = ' // XERN3 // ' AND STEPSIZE H = ' // XERN4 //
1462
* ' THE CORRECTOR COULD NOT CONVERGE BECAUSE IRES WAS EQUAL ' //
1463
* 'TO MINUS ONE', IDID, 1)
1466
C FAILURE BECAUSE IRES = -2
1467
680 WRITE (XERN3, '(1P,D15.6)') TN
1468
WRITE (XERN4, '(1P,D15.6)') H
1469
CALL XERMSG ('SLATEC', 'DDASSL',
1470
* 'AT T = ' // XERN3 // ' AND STEPSIZE H = ' // XERN4 //
1471
* ' IRES WAS EQUAL TO MINUS TWO', IDID, 1)
1474
C FAILED TO COMPUTE INITIAL YPRIME
1475
685 WRITE (XERN3, '(1P,D15.6)') TN
1476
WRITE (XERN4, '(1P,D15.6)') HO
1477
CALL XERMSG ('SLATEC', 'DDASSL',
1478
* 'AT T = ' // XERN3 // ' AND STEPSIZE H = ' // XERN4 //
1479
* ' THE INITIAL YPRIME COULD NOT BE COMPUTED', IDID, 1)
1489
C-----------------------------------------------------------------------
1490
C THIS BLOCK HANDLES ALL ERROR RETURNS DUE
1491
C TO ILLEGAL INPUT, AS DETECTED BEFORE CALLING
1492
C DDASTP. FIRST THE ERROR MESSAGE ROUTINE IS
1493
C CALLED. IF THIS HAPPENS TWICE IN
1494
C SUCCESSION, EXECUTION IS TERMINATED
1496
C-----------------------------------------------------------------------
1497
701 CALL XERMSG ('SLATEC', 'DDASSL',
1498
* 'SOME ELEMENT OF INFO VECTOR IS NOT ZERO OR ONE', 1, 1)
1501
702 WRITE (XERN1, '(I8)') NEQ
1502
CALL XERMSG ('SLATEC', 'DDASSL',
1503
* 'NEQ = ' // XERN1 // ' .LE. 0', 2, 1)
1506
703 WRITE (XERN1, '(I8)') MXORD
1507
CALL XERMSG ('SLATEC', 'DDASSL',
1508
* 'MAXORD = ' // XERN1 // ' NOT IN RANGE', 3, 1)
1511
704 WRITE (XERN1, '(I8)') LENRW
1512
WRITE (XERN2, '(I8)') LRW
1513
CALL XERMSG ('SLATEC', 'DDASSL',
1514
* 'RWORK LENGTH NEEDED, LENRW = ' // XERN1 //
1515
* ', EXCEEDS LRW = ' // XERN2, 4, 1)
1518
705 WRITE (XERN1, '(I8)') LENIW
1519
WRITE (XERN2, '(I8)') LIW
1520
CALL XERMSG ('SLATEC', 'DDASSL',
1521
* 'IWORK LENGTH NEEDED, LENIW = ' // XERN1 //
1522
* ', EXCEEDS LIW = ' // XERN2, 5, 1)
1525
706 CALL XERMSG ('SLATEC', 'DDASSL',
1526
* 'SOME ELEMENT OF RTOL IS .LT. 0', 6, 1)
1529
707 CALL XERMSG ('SLATEC', 'DDASSL',
1530
* 'SOME ELEMENT OF ATOL IS .LT. 0', 7, 1)
1533
708 CALL XERMSG ('SLATEC', 'DDASSL',
1534
* 'ALL ELEMENTS OF RTOL AND ATOL ARE ZERO', 8, 1)
1537
709 WRITE (XERN3, '(1P,D15.6)') TSTOP
1538
WRITE (XERN4, '(1P,D15.6)') TOUT
1539
CALL XERMSG ('SLATEC', 'DDASSL',
1540
* 'INFO(4) = 1 AND TSTOP = ' // XERN3 // ' BEHIND TOUT = ' //
1544
710 WRITE (XERN3, '(1P,D15.6)') HMAX
1545
CALL XERMSG ('SLATEC', 'DDASSL',
1546
* 'HMAX = ' // XERN3 // ' .LT. 0.0', 10, 1)
1549
711 WRITE (XERN3, '(1P,D15.6)') TOUT
1550
WRITE (XERN4, '(1P,D15.6)') T
1551
CALL XERMSG ('SLATEC', 'DDASSL',
1552
* 'TOUT = ' // XERN3 // ' BEHIND T = ' // XERN4, 11, 1)
1555
712 CALL XERMSG ('SLATEC', 'DDASSL',
1556
* 'INFO(8)=1 AND H0=0.0', 12, 1)
1559
713 CALL XERMSG ('SLATEC', 'DDASSL',
1560
* 'SOME ELEMENT OF WT IS .LE. 0.0', 13, 1)
1563
714 WRITE (XERN3, '(1P,D15.6)') TOUT
1564
WRITE (XERN4, '(1P,D15.6)') T
1565
CALL XERMSG ('SLATEC', 'DDASSL',
1566
* 'TOUT = ' // XERN3 // ' TOO CLOSE TO T = ' // XERN4 //
1567
* ' TO START INTEGRATION', 14, 1)
1570
715 WRITE (XERN3, '(1P,D15.6)') TSTOP
1571
WRITE (XERN4, '(1P,D15.6)') T
1572
CALL XERMSG ('SLATEC', 'DDASSL',
1573
* 'INFO(4)=1 AND TSTOP = ' // XERN3 // ' BEHIND T = ' // XERN4,
1577
717 WRITE (XERN1, '(I8)') IWORK(LML)
1578
CALL XERMSG ('SLATEC', 'DDASSL',
1579
* 'ML = ' // XERN1 // ' ILLEGAL. EITHER .LT. 0 OR .GT. NEQ',
1583
718 WRITE (XERN1, '(I8)') IWORK(LMU)
1584
CALL XERMSG ('SLATEC', 'DDASSL',
1585
* 'MU = ' // XERN1 // ' ILLEGAL. EITHER .LT. 0 OR .GT. NEQ',
1589
719 WRITE (XERN3, '(1P,D15.6)') TOUT
1590
CALL XERMSG ('SLATEC', 'DDASSL',
1591
* 'TOUT = T = ' // XERN3, 19, 1)
1595
IF(INFO(1).EQ.-1) THEN
1596
CALL XERMSG ('SLATEC', 'DDASSL',
1597
* 'REPEATED OCCURRENCES OF ILLEGAL INPUT$$' //
1598
* 'RUN TERMINATED. APPARENT INFINITE LOOP', -999, 2)
1603
C-----------END OF SUBROUTINE DDASSL------------------------------------
1605
SUBROUTINE DDAWTS (NEQ, IWT, RTOL, ATOL, Y, WT, RPAR, IPAR)
1606
C***BEGIN PROLOGUE DDAWTS
1608
C***PURPOSE Set error weight vector for DDASSL.
1609
C***LIBRARY SLATEC (DASSL)
1610
C***TYPE DOUBLE PRECISION (SDAWTS-S, DDAWTS-D)
1611
C***AUTHOR PETZOLD, LINDA R., (LLNL)
1613
C-----------------------------------------------------------------------
1614
C THIS SUBROUTINE SETS THE ERROR WEIGHT VECTOR
1615
C WT ACCORDING TO WT(I)=RTOL(I)*ABS(Y(I))+ATOL(I),
1617
C RTOL AND ATOL ARE SCALARS IF IWT = 0,
1618
C AND VECTORS IF IWT = 1.
1619
C-----------------------------------------------------------------------
1620
C***ROUTINES CALLED (NONE)
1621
C***REVISION HISTORY (YYMMDD)
1622
C 830315 DATE WRITTEN
1623
C 901009 Finished conversion to SLATEC 4.0 format (F.N.Fritsch)
1624
C 901019 Merged changes made by C. Ulrich with SLATEC 4.0 format.
1625
C 901026 Added explicit declarations for all variables and minor
1626
C cosmetic changes to prologue. (FNF)
1627
C***END PROLOGUE DDAWTS
1629
INTEGER NEQ, IWT, IPAR(*)
1630
DOUBLE PRECISION RTOL(*), ATOL(*), Y(*), WT(*), RPAR(*)
1633
DOUBLE PRECISION ATOLI, RTOLI
1635
C***FIRST EXECUTABLE STATEMENT DDAWTS
1639
IF (IWT .EQ.0) GO TO 10
1642
10 WT(I)=RTOLI*ABS(Y(I))+ATOLI
1645
C-----------END OF SUBROUTINE DDAWTS------------------------------------
1647
DOUBLE PRECISION FUNCTION DDANRM (NEQ, V, WT, RPAR, IPAR)
1648
C***BEGIN PROLOGUE DDANRM
1650
C***PURPOSE Compute vector norm for DDASSL.
1651
C***LIBRARY SLATEC (DASSL)
1652
C***TYPE DOUBLE PRECISION (SDANRM-S, DDANRM-D)
1653
C***AUTHOR PETZOLD, LINDA R., (LLNL)
1655
C-----------------------------------------------------------------------
1656
C THIS FUNCTION ROUTINE COMPUTES THE WEIGHTED
1657
C ROOT-MEAN-SQUARE NORM OF THE VECTOR OF LENGTH
1658
C NEQ CONTAINED IN THE ARRAY V,WITH WEIGHTS
1659
C CONTAINED IN THE ARRAY WT OF LENGTH NEQ.
1660
C DDANRM=SQRT((1/NEQ)*SUM(V(I)/WT(I))**2)
1661
C-----------------------------------------------------------------------
1662
C***ROUTINES CALLED (NONE)
1663
C***REVISION HISTORY (YYMMDD)
1664
C 830315 DATE WRITTEN
1665
C 901009 Finished conversion to SLATEC 4.0 format (F.N.Fritsch)
1666
C 901019 Merged changes made by C. Ulrich with SLATEC 4.0 format.
1667
C 901026 Added explicit declarations for all variables and minor
1668
C cosmetic changes to prologue. (FNF)
1669
C***END PROLOGUE DDANRM
1671
INTEGER NEQ, IPAR(*)
1672
DOUBLE PRECISION V(NEQ), WT(NEQ), RPAR(*)
1675
DOUBLE PRECISION SUM, VMAX
1677
C***FIRST EXECUTABLE STATEMENT DDANRM
1681
IF(ABS(V(I)/WT(I)) .GT. VMAX) VMAX = ABS(V(I)/WT(I))
1683
IF(VMAX .LE. 0.0D0) GO TO 30
1686
20 SUM = SUM + ((V(I)/WT(I))/VMAX)**2
1687
DDANRM = VMAX*SQRT(SUM/NEQ)
1690
C------END OF FUNCTION DDANRM------
1692
SUBROUTINE DDAINI (X, Y, YPRIME, NEQ, RES, JAC, H, WT, IDID, RPAR,
1693
+ IPAR, PHI, DELTA, E, WM, IWM, HMIN, UROUND, NONNEG, NTEMP)
1694
C***BEGIN PROLOGUE DDAINI
1696
C***PURPOSE Initialization routine for DDASSL.
1697
C***LIBRARY SLATEC (DASSL)
1698
C***TYPE DOUBLE PRECISION (SDAINI-S, DDAINI-D)
1699
C***AUTHOR PETZOLD, LINDA R., (LLNL)
1701
C-----------------------------------------------------------------
1702
C DDAINI TAKES ONE STEP OF SIZE H OR SMALLER
1703
C WITH THE BACKWARD EULER METHOD, TO
1704
C FIND YPRIME. X AND Y ARE UPDATED TO BE CONSISTENT WITH THE
1705
C NEW STEP. A MODIFIED DAMPED NEWTON ITERATION IS USED TO
1706
C SOLVE THE CORRECTOR ITERATION.
1708
C THE INITIAL GUESS FOR YPRIME IS USED IN THE
1709
C PREDICTION, AND IN FORMING THE ITERATION
1710
C MATRIX, BUT IS NOT INVOLVED IN THE
1711
C ERROR TEST. THIS MAY HAVE TROUBLE
1712
C CONVERGING IF THE INITIAL GUESS IS NO
1713
C GOOD, OR IF G(X,Y,YPRIME) DEPENDS
1714
C NONLINEARLY ON YPRIME.
1716
C THE PARAMETERS REPRESENT:
1717
C X -- INDEPENDENT VARIABLE
1718
C Y -- SOLUTION VECTOR AT X
1719
C YPRIME -- DERIVATIVE OF SOLUTION VECTOR
1720
C NEQ -- NUMBER OF EQUATIONS
1721
C H -- STEPSIZE. IMDER MAY USE A STEPSIZE
1723
C WT -- VECTOR OF WEIGHTS FOR ERROR
1725
C IDID -- COMPLETION CODE WITH THE FOLLOWING MEANINGS
1726
C IDID= 1 -- YPRIME WAS FOUND SUCCESSFULLY
1727
C IDID=-12 -- DDAINI FAILED TO FIND YPRIME
1728
C RPAR,IPAR -- REAL AND INTEGER PARAMETER ARRAYS
1729
C THAT ARE NOT ALTERED BY DDAINI
1730
C PHI -- WORK SPACE FOR DDAINI
1731
C DELTA,E -- WORK SPACE FOR DDAINI
1732
C WM,IWM -- REAL AND INTEGER ARRAYS STORING
1733
C MATRIX INFORMATION
1735
C-----------------------------------------------------------------
1736
C***ROUTINES CALLED DDAJAC, DDANRM, DDASLV
1737
C***REVISION HISTORY (YYMMDD)
1738
C 830315 DATE WRITTEN
1739
C 901009 Finished conversion to SLATEC 4.0 format (F.N.Fritsch)
1740
C 901019 Merged changes made by C. Ulrich with SLATEC 4.0 format.
1741
C 901026 Added explicit declarations for all variables and minor
1742
C cosmetic changes to prologue. (FNF)
1743
C 901030 Minor corrections to declarations. (FNF)
1744
C***END PROLOGUE DDAINI
1746
INTEGER NEQ, IDID, IPAR(*), IWM(*), NONNEG, NTEMP
1748
* X, Y(*), YPRIME(*), H, WT(*), RPAR(*), PHI(NEQ,*), DELTA(*),
1749
* E(*), WM(*), HMIN, UROUND
1752
EXTERNAL DDAJAC, DDANRM, DDASLV
1753
DOUBLE PRECISION DDANRM
1755
INTEGER I, IER, IRES, JCALC, LNJE, LNRE, M, MAXIT, MJAC, NCF,
1758
* CJ, DAMP, DELNRM, ERR, OLDNRM, R, RATE, S, XOLD, YNORM
1764
DATA MAXIT/10/,MJAC/5/
1768
C---------------------------------------------------
1771
C---------------------------------------------------
1773
C***FIRST EXECUTABLE STATEMENT DDAINI
1779
YNORM=DDANRM(NEQ,Y,WT,RPAR,IPAR)
1781
C SAVE Y AND YPRIME IN PHI
1784
100 PHI(I,2)=YPRIME(I)
1787
C----------------------------------------------------
1789
C DO ONE BACKWARD EULER STEP.
1790
C----------------------------------------------------
1792
C SET UP FOR START OF CORRECTOR ITERATION
1796
C PREDICT SOLUTION AND DERIVATIVE
1798
250 Y(I)=Y(I)+H*YPRIME(I)
1806
300 IWM(LNRE)=IWM(LNRE)+1
1809
CALL RES(X,Y,YPRIME,DELTA,IRES,RPAR,IPAR)
1810
IF (IRES.LT.0) GO TO 430
1813
C EVALUATE THE ITERATION MATRIX
1814
IF (JCALC.NE.-1) GO TO 310
1815
IWM(LNJE)=IWM(LNJE)+1
1817
CALL DDAJAC(NEQ,X,Y,YPRIME,DELTA,CJ,H,
1818
* IER,WT,E,WM,IWM,RES,IRES,
1819
* UROUND,JAC,RPAR,IPAR,NTEMP)
1822
IF (IRES.LT.0) GO TO 430
1823
IF (IER.NE.0) GO TO 430
1828
C MULTIPLY RESIDUAL BY DAMPING FACTOR
1831
320 DELTA(I)=DELTA(I)*DAMP
1833
C COMPUTE A NEW ITERATE (BACK SUBSTITUTION)
1834
C STORE THE CORRECTION IN DELTA
1836
CALL DDASLV(NEQ,DELTA,WM,IWM)
1838
C UPDATE Y AND YPRIME
1841
330 YPRIME(I)=YPRIME(I)-CJ*DELTA(I)
1843
C TEST FOR CONVERGENCE OF THE ITERATION.
1845
DELNRM=DDANRM(NEQ,DELTA,WT,RPAR,IPAR)
1846
IF (DELNRM.LE.100.D0*UROUND*YNORM)
1849
IF (M.GT.0) GO TO 340
1853
340 RATE=(DELNRM/OLDNRM)**(1.0D0/M)
1854
IF (RATE.GT.0.90D0) GO TO 430
1857
350 IF (S*DELNRM .LE. 0.33D0) GO TO 400
1860
C THE CORRECTOR HAS NOT YET CONVERGED. UPDATE
1861
C M AND AND TEST WHETHER THE MAXIMUM
1862
C NUMBER OF ITERATIONS HAVE BEEN TRIED.
1863
C EVERY MJAC ITERATIONS, GET A NEW
1867
IF (M.GE.MAXIT) GO TO 430
1869
IF ((M/MJAC)*MJAC.EQ.M) JCALC=-1
1873
C THE ITERATION HAS CONVERGED.
1874
C CHECK NONNEGATIVITY CONSTRAINTS
1875
400 IF (NONNEG.EQ.0) GO TO 450
1877
410 DELTA(I)=MIN(Y(I),0.0D0)
1879
DELNRM=DDANRM(NEQ,DELTA,WT,RPAR,IPAR)
1880
IF (DELNRM.GT.0.33D0) GO TO 430
1884
420 YPRIME(I)=YPRIME(I)-CJ*DELTA(I)
1888
C EXITS FROM CORRECTOR LOOP.
1890
450 IF (.NOT.CONVGD) GO TO 600
1894
C-----------------------------------------------------
1896
C THE CORRECTOR ITERATION CONVERGED.
1898
C-----------------------------------------------------
1901
510 E(I)=Y(I)-PHI(I,1)
1902
ERR=DDANRM(NEQ,E,WT,RPAR,IPAR)
1904
IF (ERR.LE.1.0D0) RETURN
1908
C--------------------------------------------------------
1910
C THE BACKWARD EULER STEP FAILED. RESTORE X, Y
1911
C AND YPRIME TO THEIR ORIGINAL VALUES.
1912
C REDUCE STEPSIZE AND TRY AGAIN, IF
1914
C---------------------------------------------------------
1920
610 YPRIME(I)=PHI(I,2)
1922
IF (CONVGD) GO TO 640
1923
IF (IER.EQ.0) GO TO 620
1926
IF (NSF.LT.3.AND.ABS(H).GE.HMIN) GO TO 690
1929
620 IF (IRES.GT.-2) GO TO 630
1934
IF (NCF.LT.10.AND.ABS(H).GE.HMIN) GO TO 690
1939
R=0.90D0/(2.0D0*ERR+0.0001D0)
1940
R=MAX(0.1D0,MIN(0.5D0,R))
1942
IF (ABS(H).GE.HMIN.AND.NEF.LT.10) GO TO 690
1947
C-------------END OF SUBROUTINE DDAINI----------------------
1949
SUBROUTINE DDATRP (X, XOUT, YOUT, YPOUT, NEQ, KOLD, PHI, PSI)
1950
C***BEGIN PROLOGUE DDATRP
1952
C***PURPOSE Interpolation routine for DDASSL.
1953
C***LIBRARY SLATEC (DASSL)
1954
C***TYPE DOUBLE PRECISION (SDATRP-S, DDATRP-D)
1955
C***AUTHOR PETZOLD, LINDA R., (LLNL)
1957
C-----------------------------------------------------------------------
1958
C THE METHODS IN SUBROUTINE DDASTP USE POLYNOMIALS
1959
C TO APPROXIMATE THE SOLUTION. DDATRP APPROXIMATES THE
1960
C SOLUTION AND ITS DERIVATIVE AT TIME XOUT BY EVALUATING
1961
C ONE OF THESE POLYNOMIALS,AND ITS DERIVATIVE,THERE.
1962
C INFORMATION DEFINING THIS POLYNOMIAL IS PASSED FROM
1963
C DDASTP, SO DDATRP CANNOT BE USED ALONE.
1965
C THE PARAMETERS ARE:
1966
C X THE CURRENT TIME IN THE INTEGRATION.
1967
C XOUT THE TIME AT WHICH THE SOLUTION IS DESIRED
1968
C YOUT THE INTERPOLATED APPROXIMATION TO Y AT XOUT
1970
C YPOUT THE INTERPOLATED APPROXIMATION TO YPRIME AT XOUT
1972
C NEQ NUMBER OF EQUATIONS
1973
C KOLD ORDER USED ON LAST SUCCESSFUL STEP
1974
C PHI ARRAY OF SCALED DIVIDED DIFFERENCES OF Y
1975
C PSI ARRAY OF PAST STEPSIZE HISTORY
1976
C-----------------------------------------------------------------------
1977
C***ROUTINES CALLED (NONE)
1978
C***REVISION HISTORY (YYMMDD)
1979
C 830315 DATE WRITTEN
1980
C 901009 Finished conversion to SLATEC 4.0 format (F.N.Fritsch)
1981
C 901019 Merged changes made by C. Ulrich with SLATEC 4.0 format.
1982
C 901026 Added explicit declarations for all variables and minor
1983
C cosmetic changes to prologue. (FNF)
1984
C***END PROLOGUE DDATRP
1987
DOUBLE PRECISION X, XOUT, YOUT(*), YPOUT(*), PHI(NEQ,*), PSI(*)
1989
INTEGER I, J, KOLDP1
1990
DOUBLE PRECISION C, D, GAMMA, TEMP1
1992
C***FIRST EXECUTABLE STATEMENT DDATRP
2002
D=D*GAMMA+C/PSI(J-1)
2004
GAMMA=(TEMP1+PSI(J-1))/PSI(J)
2006
YOUT(I)=YOUT(I)+C*PHI(I,J)
2007
20 YPOUT(I)=YPOUT(I)+D*PHI(I,J)
2011
C------END OF SUBROUTINE DDATRP------
2013
SUBROUTINE DDASTP (X, Y, YPRIME, NEQ, RES, JAC, H, WT, JSTART,
2014
+ IDID, RPAR, IPAR, PHI, DELTA, E, WM, IWM, ALPHA, BETA, GAMMA,
2015
+ PSI, SIGMA, CJ, CJOLD, HOLD, S, HMIN, UROUND, IPHASE, JCALC,
2016
+ K, KOLD, NS, NONNEG, NTEMP)
2017
C***BEGIN PROLOGUE DDASTP
2019
C***PURPOSE Perform one step of the DDASSL integration.
2020
C***LIBRARY SLATEC (DASSL)
2021
C***TYPE DOUBLE PRECISION (SDASTP-S, DDASTP-D)
2022
C***AUTHOR PETZOLD, LINDA R., (LLNL)
2024
C-----------------------------------------------------------------------
2025
C DDASTP SOLVES A SYSTEM OF DIFFERENTIAL/
2026
C ALGEBRAIC EQUATIONS OF THE FORM
2027
C G(X,Y,YPRIME) = 0, FOR ONE STEP (NORMALLY
2030
C THE METHODS USED ARE MODIFIED DIVIDED
2031
C DIFFERENCE,FIXED LEADING COEFFICIENT
2032
C FORMS OF BACKWARD DIFFERENTIATION
2033
C FORMULAS. THE CODE ADJUSTS THE STEPSIZE
2034
C AND ORDER TO CONTROL THE LOCAL ERROR PER
2038
C THE PARAMETERS REPRESENT
2039
C X -- INDEPENDENT VARIABLE
2040
C Y -- SOLUTION VECTOR AT X
2041
C YPRIME -- DERIVATIVE OF SOLUTION VECTOR
2042
C AFTER SUCCESSFUL STEP
2043
C NEQ -- NUMBER OF EQUATIONS TO BE INTEGRATED
2044
C RES -- EXTERNAL USER-SUPPLIED SUBROUTINE
2045
C TO EVALUATE THE RESIDUAL. THE CALL IS
2046
C CALL RES(X,Y,YPRIME,DELTA,IRES,RPAR,IPAR)
2047
C X,Y,YPRIME ARE INPUT. DELTA IS OUTPUT.
2048
C ON INPUT, IRES=0. RES SHOULD ALTER IRES ONLY
2049
C IF IT ENCOUNTERS AN ILLEGAL VALUE OF Y OR A
2050
C STOP CONDITION. SET IRES=-1 IF AN INPUT VALUE
2051
C OF Y IS ILLEGAL, AND DDASTP WILL TRY TO SOLVE
2052
C THE PROBLEM WITHOUT GETTING IRES = -1. IF
2053
C IRES=-2, DDASTP RETURNS CONTROL TO THE CALLING
2054
C PROGRAM WITH IDID = -11.
2055
C JAC -- EXTERNAL USER-SUPPLIED ROUTINE TO EVALUATE
2056
C THE ITERATION MATRIX (THIS IS OPTIONAL)
2057
C THE CALL IS OF THE FORM
2058
C CALL JAC(X,Y,YPRIME,PD,CJ,RPAR,IPAR)
2059
C PD IS THE MATRIX OF PARTIAL DERIVATIVES,
2060
C PD=DG/DY+CJ*DG/DYPRIME
2061
C H -- APPROPRIATE STEP SIZE FOR NEXT STEP.
2062
C NORMALLY DETERMINED BY THE CODE
2063
C WT -- VECTOR OF WEIGHTS FOR ERROR CRITERION.
2064
C JSTART -- INTEGER VARIABLE SET 0 FOR
2065
C FIRST STEP, 1 OTHERWISE.
2066
C IDID -- COMPLETION CODE WITH THE FOLLOWING MEANINGS:
2067
C IDID= 1 -- THE STEP WAS COMPLETED SUCCESSFULLY
2068
C IDID=-6 -- THE ERROR TEST FAILED REPEATEDLY
2069
C IDID=-7 -- THE CORRECTOR COULD NOT CONVERGE
2070
C IDID=-8 -- THE ITERATION MATRIX IS SINGULAR
2071
C IDID=-9 -- THE CORRECTOR COULD NOT CONVERGE.
2072
C THERE WERE REPEATED ERROR TEST
2073
C FAILURES ON THIS STEP.
2074
C IDID=-10-- THE CORRECTOR COULD NOT CONVERGE
2075
C BECAUSE IRES WAS EQUAL TO MINUS ONE
2076
C IDID=-11-- IRES EQUAL TO -2 WAS ENCOUNTERED,
2077
C AND CONTROL IS BEING RETURNED TO
2078
C THE CALLING PROGRAM
2079
C RPAR,IPAR -- REAL AND INTEGER PARAMETER ARRAYS THAT
2080
C ARE USED FOR COMMUNICATION BETWEEN THE
2081
C CALLING PROGRAM AND EXTERNAL USER ROUTINES
2082
C THEY ARE NOT ALTERED BY DDASTP
2083
C PHI -- ARRAY OF DIVIDED DIFFERENCES USED BY
2084
C DDASTP. THE LENGTH IS NEQ*(K+1),WHERE
2085
C K IS THE MAXIMUM ORDER
2086
C DELTA,E -- WORK VECTORS FOR DDASTP OF LENGTH NEQ
2087
C WM,IWM -- REAL AND INTEGER ARRAYS STORING
2088
C MATRIX INFORMATION SUCH AS THE MATRIX
2089
C OF PARTIAL DERIVATIVES,PERMUTATION
2090
C VECTOR,AND VARIOUS OTHER INFORMATION.
2092
C THE OTHER PARAMETERS ARE INFORMATION
2093
C WHICH IS NEEDED INTERNALLY BY DDASTP TO
2094
C CONTINUE FROM STEP TO STEP.
2096
C-----------------------------------------------------------------------
2097
C***ROUTINES CALLED DDAJAC, DDANRM, DDASLV, DDATRP
2098
C***REVISION HISTORY (YYMMDD)
2099
C 830315 DATE WRITTEN
2100
C 901009 Finished conversion to SLATEC 4.0 format (F.N.Fritsch)
2101
C 901019 Merged changes made by C. Ulrich with SLATEC 4.0 format.
2102
C 901026 Added explicit declarations for all variables and minor
2103
C cosmetic changes to prologue. (FNF)
2104
C***END PROLOGUE DDASTP
2106
INTEGER NEQ, JSTART, IDID, IPAR(*), IWM(*), IPHASE, JCALC, K,
2107
* KOLD, NS, NONNEG, NTEMP
2109
* X, Y(*), YPRIME(*), H, WT(*), RPAR(*), PHI(NEQ,*), DELTA(*),
2110
* E(*), WM(*), ALPHA(*), BETA(*), GAMMA(*), PSI(*), SIGMA(*), CJ,
2111
* CJOLD, HOLD, S, HMIN, UROUND
2114
EXTERNAL DDAJAC, DDANRM, DDASLV, DDATRP
2115
DOUBLE PRECISION DDANRM
2117
INTEGER I, IER, IRES, J, J1, KDIFF, KM1, KNEW, KP1, KP2, LCTF,
2118
* LETF, LMXORD, LNJE, LNRE, LNST, M, MAXIT, NCF, NEF, NSF, NSP1
2120
* ALPHA0, ALPHAS, CJLAST, CK, DELNRM, ENORM, ERK, ERKM1,
2121
* ERKM2, ERKP1, ERR, EST, HNEW, OLDNRM, PNORM, R, RATE, TEMP1,
2122
* TEMP2, TERK, TERKM1, TERKM2, TERKP1, XOLD, XRATE
2125
PARAMETER (LMXORD=3)
2139
C-----------------------------------------------------------------------
2141
C INITIALIZE. ON THE FIRST CALL,SET
2142
C THE ORDER TO 1 AND INITIALIZE
2144
C-----------------------------------------------------------------------
2146
C INITIALIZATIONS FOR ALL CALLS
2147
C***FIRST EXECUTABLE STATEMENT DDASTP
2153
IF(JSTART .NE. 0) GO TO 120
2155
C IF THIS IS THE FIRST STEP,PERFORM
2156
C OTHER INITIALIZATIONS
2177
C-----------------------------------------------------------------------
2179
C COMPUTE COEFFICIENTS OF FORMULAS FOR
2181
C-----------------------------------------------------------------------
2187
IF(H.NE.HOLD.OR.K .NE. KOLD) NS = 0
2190
IF(KP1 .LT. NS)GO TO 230
2200
BETA(I)=BETA(I-1)*PSI(I-1)/TEMP2
2203
SIGMA(I)=(I-1)*SIGMA(I-1)*ALPHA(I)
2204
GAMMA(I)=GAMMA(I-1)+ALPHA(I-1)/H
2209
C COMPUTE ALPHAS, ALPHA0
2213
ALPHAS = ALPHAS - 1.0D0/I
2214
ALPHA0 = ALPHA0 - ALPHA(I)
2217
C COMPUTE LEADING COEFFICIENT CJ
2221
C COMPUTE VARIABLE STEPSIZE ERROR COEFFICIENT CK
2222
CK = ABS(ALPHA(KP1) + ALPHAS - ALPHA0)
2223
CK = MAX(CK,ALPHA(KP1))
2225
C DECIDE WHETHER NEW JACOBIAN IS NEEDED
2226
TEMP1 = (1.0D0 - XRATE)/(1.0D0 + XRATE)
2228
IF (CJ/CJOLD .LT. TEMP1 .OR. CJ/CJOLD .GT. TEMP2) JCALC = -1
2229
IF (CJ .NE. CJLAST) S = 100.D0
2231
C CHANGE PHI TO PHI STAR
2232
IF(KP1 .LT. NSP1) GO TO 280
2235
260 PHI(I,J)=BETA(J)*PHI(I,J)
2246
C-----------------------------------------------------------------------
2248
C PREDICT THE SOLUTION AND DERIVATIVE,
2249
C AND SOLVE THE CORRECTOR EQUATION
2250
C-----------------------------------------------------------------------
2252
C FIRST,PREDICT THE SOLUTION AND DERIVATIVE
2260
320 YPRIME(I)=YPRIME(I)+GAMMA(J)*PHI(I,J)
2262
PNORM = DDANRM (NEQ,Y,WT,RPAR,IPAR)
2266
C SOLVE THE CORRECTOR EQUATION USING A
2267
C MODIFIED NEWTON SCHEME.
2270
IWM(LNRE)=IWM(LNRE)+1
2272
CALL RES(X,Y,YPRIME,DELTA,IRES,RPAR,IPAR)
2273
IF (IRES .LT. 0) GO TO 380
2276
C IF INDICATED,REEVALUATE THE
2277
C ITERATION MATRIX PD = DG/DY + CJ*DG/DYPRIME
2278
C (WHERE G(X,Y,YPRIME)=0). SET
2279
C JCALC TO 0 AS AN INDICATOR THAT
2280
C THIS HAS BEEN DONE.
2281
IF(JCALC .NE. -1)GO TO 340
2282
IWM(LNJE)=IWM(LNJE)+1
2284
CALL DDAJAC(NEQ,X,Y,YPRIME,DELTA,CJ,H,
2285
* IER,WT,E,WM,IWM,RES,IRES,UROUND,JAC,RPAR,
2289
IF (IRES .LT. 0) GO TO 380
2290
IF(IER .NE. 0)GO TO 380
2294
C INITIALIZE THE ERROR ACCUMULATION VECTOR E.
2303
C MULTIPLY RESIDUAL BY TEMP1 TO ACCELERATE CONVERGENCE
2304
TEMP1 = 2.0D0/(1.0D0 + CJ/CJOLD)
2306
355 DELTA(I) = DELTA(I) * TEMP1
2308
C COMPUTE A NEW ITERATE (BACK-SUBSTITUTION).
2309
C STORE THE CORRECTION IN DELTA.
2310
CALL DDASLV(NEQ,DELTA,WM,IWM)
2312
C UPDATE Y,E,AND YPRIME
2316
360 YPRIME(I)=YPRIME(I)-CJ*DELTA(I)
2318
C TEST FOR CONVERGENCE OF THE ITERATION
2319
DELNRM=DDANRM(NEQ,DELTA,WT,RPAR,IPAR)
2320
IF (DELNRM .LE. 100.D0*UROUND*PNORM) GO TO 375
2321
IF (M .GT. 0) GO TO 365
2324
365 RATE = (DELNRM/OLDNRM)**(1.0D0/M)
2325
IF (RATE .GT. 0.90D0) GO TO 370
2326
S = RATE/(1.0D0 - RATE)
2327
367 IF (S*DELNRM .LE. 0.33D0) GO TO 375
2329
C THE CORRECTOR HAS NOT YET CONVERGED.
2330
C UPDATE M AND TEST WHETHER THE
2331
C MAXIMUM NUMBER OF ITERATIONS HAVE
2334
IF(M.GE.MAXIT)GO TO 370
2336
C EVALUATE THE RESIDUAL
2337
C AND GO BACK TO DO ANOTHER ITERATION
2338
IWM(LNRE)=IWM(LNRE)+1
2340
CALL RES(X,Y,YPRIME,DELTA,IRES,
2342
IF (IRES .LT. 0) GO TO 380
2346
C THE CORRECTOR FAILED TO CONVERGE IN MAXIT
2347
C ITERATIONS. IF THE ITERATION MATRIX
2348
C IS NOT CURRENT,RE-DO THE STEP WITH
2349
C A NEW ITERATION MATRIX.
2351
IF(JCALC.EQ.0)GO TO 380
2356
C THE ITERATION HAS CONVERGED. IF NONNEGATIVITY OF SOLUTION IS
2357
C REQUIRED, SET THE SOLUTION NONNEGATIVE, IF THE PERTURBATION
2358
C TO DO IT IS SMALL ENOUGH. IF THE CHANGE IS TOO LARGE, THEN
2359
C CONSIDER THE CORRECTOR ITERATION TO HAVE FAILED.
2360
375 IF(NONNEG .EQ. 0) GO TO 390
2362
377 DELTA(I) = MIN(Y(I),0.0D0)
2363
DELNRM = DDANRM(NEQ,DELTA,WT,RPAR,IPAR)
2364
IF(DELNRM .GT. 0.33D0) GO TO 380
2366
378 E(I) = E(I) - DELTA(I)
2370
C EXITS FROM BLOCK 3
2371
C NO CONVERGENCE WITH CURRENT ITERATION
2372
C MATRIX,OR SINGULAR ITERATION MATRIX
2375
IF(.NOT.CONVGD)GO TO 600
2381
C-----------------------------------------------------------------------
2383
C ESTIMATE THE ERRORS AT ORDERS K,K-1,K-2
2384
C AS IF CONSTANT STEPSIZE WAS USED. ESTIMATE
2385
C THE LOCAL ERROR AT ORDER K AND TEST
2386
C WHETHER THE CURRENT STEP IS SUCCESSFUL.
2387
C-----------------------------------------------------------------------
2389
C ESTIMATE ERRORS AT ORDERS K,K-1,K-2
2390
ENORM = DDANRM(NEQ,E,WT,RPAR,IPAR)
2391
ERK = SIGMA(K+1)*ENORM
2395
IF(K .EQ. 1)GO TO 430
2397
405 DELTA(I) = PHI(I,KP1) + E(I)
2398
ERKM1=SIGMA(K)*DDANRM(NEQ,DELTA,WT,RPAR,IPAR)
2400
IF(K .GT. 2)GO TO 410
2401
IF(TERKM1 .LE. 0.5D0*TERK)GO TO 420
2405
415 DELTA(I) = PHI(I,K) + DELTA(I)
2406
ERKM2=SIGMA(K-1)*DDANRM(NEQ,DELTA,WT,RPAR,IPAR)
2407
TERKM2 = (K-1)*ERKM2
2408
IF(MAX(TERKM1,TERKM2).GT.TERK)GO TO 430
2415
C CALCULATE THE LOCAL ERROR FOR THE CURRENT STEP
2416
C TO SEE IF THE STEP WAS SUCCESSFUL
2419
IF(ERR .GT. 1.0D0)GO TO 600
2425
C-----------------------------------------------------------------------
2427
C THE STEP IS SUCCESSFUL. DETERMINE
2428
C THE BEST ORDER AND STEPSIZE FOR
2429
C THE NEXT STEP. UPDATE THE DIFFERENCES
2430
C FOR THE NEXT STEP.
2431
C-----------------------------------------------------------------------
2433
IWM(LNST)=IWM(LNST)+1
2439
C ESTIMATE THE ERROR AT ORDER K+1 UNLESS:
2440
C ALREADY DECIDED TO LOWER ORDER, OR
2441
C ALREADY USING MAXIMUM ORDER, OR
2442
C STEPSIZE NOT CONSTANT, OR
2443
C ORDER RAISED IN PREVIOUS STEP
2444
IF(KNEW.EQ.KM1.OR.K.EQ.IWM(LMXORD))IPHASE=1
2445
IF(IPHASE .EQ. 0)GO TO 545
2446
IF(KNEW.EQ.KM1)GO TO 540
2447
IF(K.EQ.IWM(LMXORD)) GO TO 550
2448
IF(KP1.GE.NS.OR.KDIFF.EQ.1)GO TO 550
2450
510 DELTA(I)=E(I)-PHI(I,KP2)
2451
ERKP1 = (1.0D0/(K+2))*DDANRM(NEQ,DELTA,WT,RPAR,IPAR)
2452
TERKP1 = (K+2)*ERKP1
2454
IF(TERKP1.GE.0.5D0*TERK)GO TO 550
2456
520 IF(TERKM1.LE.MIN(TERK,TERKP1))GO TO 540
2457
IF(TERKP1.GE.TERK.OR.K.EQ.IWM(LMXORD))GO TO 550
2469
C IF IPHASE = 0, INCREASE ORDER BY ONE AND MULTIPLY STEPSIZE BY
2477
C DETERMINE THE APPROPRIATE STEPSIZE FOR
2481
R=(2.0D0*EST+0.0001D0)**(-1.0D0/TEMP2)
2482
IF(R .LT. 2.0D0) GO TO 555
2485
555 IF(R .GT. 1.0D0) GO TO 560
2486
R = MAX(0.5D0,MIN(0.9D0,R))
2491
C UPDATE DIFFERENCES FOR NEXT STEP
2493
IF(KOLD.EQ.IWM(LMXORD))GO TO 585
2498
590 PHI(I,KP1)=PHI(I,KP1)+E(I)
2502
595 PHI(I,J)=PHI(I,J)+PHI(I,J+1)
2509
C-----------------------------------------------------------------------
2511
C THE STEP IS UNSUCCESSFUL. RESTORE X,PSI,PHI
2512
C DETERMINE APPROPRIATE STEPSIZE FOR
2513
C CONTINUING THE INTEGRATION, OR EXIT WITH
2514
C AN ERROR FLAG IF THERE HAVE BEEN MANY
2516
C-----------------------------------------------------------------------
2521
IF(KP1.LT.NSP1)GO TO 630
2525
610 PHI(I,J)=TEMP1*PHI(I,J)
2529
640 PSI(I-1)=PSI(I)-H
2532
C TEST WHETHER FAILURE IS DUE TO CORRECTOR ITERATION
2535
IWM(LCTF)=IWM(LCTF)+1
2538
C THE NEWTON ITERATION FAILED TO CONVERGE WITH
2539
C A CURRENT ITERATION MATRIX. DETERMINE THE CAUSE
2540
C OF THE FAILURE AND TAKE APPROPRIATE ACTION.
2541
IF(IER.EQ.0)GO TO 650
2543
C THE ITERATION MATRIX IS SINGULAR. REDUCE
2544
C THE STEPSIZE BY A FACTOR OF 4. IF
2545
C THIS HAPPENS THREE TIMES IN A ROW ON
2546
C THE SAME STEP, RETURN WITH AN ERROR FLAG
2550
IF (NSF .LT. 3 .AND. ABS(H) .GE. HMIN) GO TO 690
2555
C THE NEWTON ITERATION FAILED TO CONVERGE FOR A REASON
2556
C OTHER THAN A SINGULAR ITERATION MATRIX. IF IRES = -2, THEN
2557
C RETURN. OTHERWISE, REDUCE THE STEPSIZE AND TRY AGAIN, UNLESS
2558
C TOO MANY FAILURES HAVE OCCURED.
2560
IF (IRES .GT. -2) GO TO 655
2566
IF (NCF .LT. 10 .AND. ABS(H) .GE. HMIN) GO TO 690
2568
IF (IRES .LT. 0) IDID = -10
2569
IF (NEF .GE. 3) IDID = -9
2573
C THE NEWTON SCHEME CONVERGED,AND THE CAUSE
2574
C OF THE FAILURE WAS THE ERROR ESTIMATE
2575
C EXCEEDING THE TOLERANCE.
2577
IWM(LETF)=IWM(LETF)+1
2578
IF (NEF .GT. 1) GO TO 665
2580
C ON FIRST ERROR TEST FAILURE, KEEP CURRENT ORDER OR LOWER
2581
C ORDER BY ONE. COMPUTE NEW STEPSIZE BASED ON DIFFERENCES
2585
R = 0.90D0*(2.0D0*EST+0.0001D0)**(-1.0D0/TEMP2)
2586
R = MAX(0.25D0,MIN(0.9D0,R))
2588
IF (ABS(H) .GE. HMIN) GO TO 690
2592
C ON SECOND ERROR TEST FAILURE, USE THE CURRENT ORDER OR
2593
C DECREASE ORDER BY ONE. REDUCE THE STEPSIZE BY A FACTOR OF
2595
665 IF (NEF .GT. 2) GO TO 670
2598
IF (ABS(H) .GE. HMIN) GO TO 690
2602
C ON THIRD AND SUBSEQUENT ERROR TEST FAILURES, SET THE ORDER TO
2603
C ONE AND REDUCE THE STEPSIZE BY A FACTOR OF FOUR.
2606
IF (ABS(H) .GE. HMIN) GO TO 690
2613
C FOR ALL CRASHES, RESTORE Y TO ITS LAST VALUE,
2614
C INTERPOLATE TO FIND YPRIME AT LAST X, AND RETURN
2616
CALL DDATRP(X,X,Y,YPRIME,NEQ,K,PHI,PSI)
2620
C GO BACK AND TRY THIS STEP AGAIN
2623
C------END OF SUBROUTINE DDASTP------
2625
SUBROUTINE DDAJAC (NEQ, X, Y, YPRIME, DELTA, CJ, H,
2626
+ IER, WT, E, WM, IWM, RES, IRES, UROUND, JAC, RPAR,
2628
C***BEGIN PROLOGUE DDAJAC
2630
C***PURPOSE Compute the iteration matrix for DDASSL and form the
2632
C***LIBRARY SLATEC (DASSL)
2633
C***TYPE DOUBLE PRECISION (SDAJAC-S, DDAJAC-D)
2634
C***AUTHOR PETZOLD, LINDA R., (LLNL)
2636
C-----------------------------------------------------------------------
2637
C THIS ROUTINE COMPUTES THE ITERATION MATRIX
2638
C PD=DG/DY+CJ*DG/DYPRIME (WHERE G(X,Y,YPRIME)=0).
2639
C HERE PD IS COMPUTED BY THE USER-SUPPLIED
2640
C ROUTINE JAC IF IWM(MTYPE) IS 1 OR 4, AND
2641
C IT IS COMPUTED BY NUMERICAL FINITE DIFFERENCING
2642
C IF IWM(MTYPE)IS 2 OR 5
2643
C THE PARAMETERS HAVE THE FOLLOWING MEANINGS.
2644
C Y = ARRAY CONTAINING PREDICTED VALUES
2645
C YPRIME = ARRAY CONTAINING PREDICTED DERIVATIVES
2646
C DELTA = RESIDUAL EVALUATED AT (X,Y,YPRIME)
2647
C (USED ONLY IF IWM(MTYPE)=2 OR 5)
2648
C CJ = SCALAR PARAMETER DEFINING ITERATION MATRIX
2649
C H = CURRENT STEPSIZE IN INTEGRATION
2650
C IER = VARIABLE WHICH IS .NE. 0
2651
C IF ITERATION MATRIX IS SINGULAR,
2653
C WT = VECTOR OF WEIGHTS FOR COMPUTING NORMS
2654
C E = WORK SPACE (TEMPORARY) OF LENGTH NEQ
2655
C WM = REAL WORK SPACE FOR MATRICES. ON
2656
C OUTPUT IT CONTAINS THE LU DECOMPOSITION
2657
C OF THE ITERATION MATRIX.
2658
C IWM = INTEGER WORK SPACE CONTAINING
2659
C MATRIX INFORMATION
2660
C RES = NAME OF THE EXTERNAL USER-SUPPLIED ROUTINE
2661
C TO EVALUATE THE RESIDUAL FUNCTION G(X,Y,YPRIME)
2662
C IRES = FLAG WHICH IS EQUAL TO ZERO IF NO ILLEGAL VALUES
2663
C IN RES, AND LESS THAN ZERO OTHERWISE. (IF IRES
2664
C IS LESS THAN ZERO, THE MATRIX WAS NOT COMPLETED)
2665
C IN THIS CASE (IF IRES .LT. 0), THEN IER = 0.
2666
C UROUND = THE UNIT ROUNDOFF ERROR OF THE MACHINE BEING USED.
2667
C JAC = NAME OF THE EXTERNAL USER-SUPPLIED ROUTINE
2668
C TO EVALUATE THE ITERATION MATRIX (THIS ROUTINE
2669
C IS ONLY USED IF IWM(MTYPE) IS 1 OR 4)
2670
C-----------------------------------------------------------------------
2671
C***ROUTINES CALLED DGBFA, DGEFA
2672
C***REVISION HISTORY (YYMMDD)
2673
C 830315 DATE WRITTEN
2674
C 901009 Finished conversion to SLATEC 4.0 format (F.N.Fritsch)
2675
C 901010 Modified three MAX calls to be all on one line. (FNF)
2676
C 901019 Merged changes made by C. Ulrich with SLATEC 4.0 format.
2677
C 901026 Added explicit declarations for all variables and minor
2678
C cosmetic changes to prologue. (FNF)
2679
C 901101 Corrected PURPOSE. (FNF)
2680
C***END PROLOGUE DDAJAC
2682
INTEGER NEQ, IER, IWM(*), IRES, IPAR(*), NTEMP
2684
* X, Y(*), YPRIME(*), DELTA(*), CJ, H, WT(*), E(*), WM(*),
2688
EXTERNAL DGBFA, DGEFA
2690
INTEGER I, I1, I2, II, IPSAVE, ISAVE, J, K, L, LENPD, LIPVT,
2691
* LML, LMTYPE, LMU, MBA, MBAND, MEB1, MEBAND, MSAVE, MTYPE, N,
2693
DOUBLE PRECISION DEL, DELINV, SQUR, YPSAVE, YSAVE
2698
PARAMETER (LMTYPE=4)
2699
PARAMETER (LIPVT=21)
2701
C***FIRST EXECUTABLE STATEMENT DDAJAC
2705
GO TO (100,200,300,400,500),MTYPE
2708
C DENSE USER-SUPPLIED MATRIX
2711
110 WM(NPDM1+I)=0.0D0
2712
CALL JAC(X,Y,YPRIME,WM(NPD),CJ,RPAR,IPAR)
2716
C DENSE FINITE-DIFFERENCE-GENERATED MATRIX
2721
DEL=SQUR*MAX(ABS(Y(I)),ABS(H*YPRIME(I)),ABS(WT(I)))
2722
DEL=SIGN(DEL,H*YPRIME(I))
2727
YPRIME(I)=YPRIME(I)+CJ*DEL
2728
CALL RES(X,Y,YPRIME,E,IRES,RPAR,IPAR)
2729
IF (IRES .LT. 0) RETURN
2732
220 WM(NROW+L)=(E(L)-DELTA(L))*DELINV
2739
C DO DENSE-MATRIX LU DECOMPOSITION ON PD
2740
230 CALL DGEFA(WM(NPD),NEQ,NEQ,IWM(LIPVT),IER)
2744
C DUMMY SECTION FOR IWM(MTYPE)=3
2748
C BANDED USER-SUPPLIED MATRIX
2749
400 LENPD=(2*IWM(LML)+IWM(LMU)+1)*NEQ
2751
410 WM(NPDM1+I)=0.0D0
2752
CALL JAC(X,Y,YPRIME,WM(NPD),CJ,RPAR,IPAR)
2753
MEBAND=2*IWM(LML)+IWM(LMU)+1
2757
C BANDED FINITE-DIFFERENCE-GENERATED MATRIX
2758
500 MBAND=IWM(LML)+IWM(LMU)+1
2760
MEBAND=MBAND+IWM(LML)
2768
DO 510 N=J,NEQ,MBAND
2771
WM(IPSAVE+K)=YPRIME(N)
2772
DEL=SQUR*MAX(ABS(Y(N)),ABS(H*YPRIME(N)),ABS(WT(N)))
2773
DEL=SIGN(DEL,H*YPRIME(N))
2776
510 YPRIME(N)=YPRIME(N)+CJ*DEL
2777
CALL RES(X,Y,YPRIME,E,IRES,RPAR,IPAR)
2778
IF (IRES .LT. 0) RETURN
2779
DO 530 N=J,NEQ,MBAND
2782
YPRIME(N)=WM(IPSAVE+K)
2783
DEL=SQUR*MAX(ABS(Y(N)),ABS(H*YPRIME(N)),ABS(WT(N)))
2784
DEL=SIGN(DEL,H*YPRIME(N))
2787
I1=MAX(1,(N-IWM(LMU)))
2788
I2=MIN(NEQ,(N+IWM(LML)))
2789
II=N*MEB1-IWM(LML)+NPDM1
2791
520 WM(II+I)=(E(I)-DELTA(I))*DELINV
2796
C DO LU DECOMPOSITION OF BANDED PD
2797
550 CALL DGBFA(WM(NPD),MEBAND,NEQ,
2798
* IWM(LML),IWM(LMU),IWM(LIPVT),IER)
2800
C------END OF SUBROUTINE DDAJAC------
2802
SUBROUTINE DDASLV (NEQ, DELTA, WM, IWM)
2803
C***BEGIN PROLOGUE DDASLV
2805
C***PURPOSE Linear system solver for DDASSL.
2806
C***LIBRARY SLATEC (DASSL)
2807
C***TYPE DOUBLE PRECISION (SDASLV-S, DDASLV-D)
2808
C***AUTHOR PETZOLD, LINDA R., (LLNL)
2810
C-----------------------------------------------------------------------
2811
C THIS ROUTINE MANAGES THE SOLUTION OF THE LINEAR
2812
C SYSTEM ARISING IN THE NEWTON ITERATION.
2813
C MATRICES AND REAL TEMPORARY STORAGE AND
2814
C REAL INFORMATION ARE STORED IN THE ARRAY WM.
2815
C INTEGER MATRIX INFORMATION IS STORED IN
2817
C FOR A DENSE MATRIX, THE LINPACK ROUTINE
2819
C FOR A BANDED MATRIX,THE LINPACK ROUTINE
2821
C-----------------------------------------------------------------------
2822
C***ROUTINES CALLED DGBSL, DGESL
2823
C***REVISION HISTORY (YYMMDD)
2824
C 830315 DATE WRITTEN
2825
C 901009 Finished conversion to SLATEC 4.0 format (F.N.Fritsch)
2826
C 901019 Merged changes made by C. Ulrich with SLATEC 4.0 format.
2827
C 901026 Added explicit declarations for all variables and minor
2828
C cosmetic changes to prologue. (FNF)
2829
C***END PROLOGUE DDASLV
2832
DOUBLE PRECISION DELTA(*), WM(*)
2834
EXTERNAL DGBSL, DGESL
2836
INTEGER LIPVT, LML, LMU, LMTYPE, MEBAND, MTYPE, NPD
2840
PARAMETER (LMTYPE=4)
2841
PARAMETER (LIPVT=21)
2843
C***FIRST EXECUTABLE STATEMENT DDASLV
2845
GO TO(100,100,300,400,400),MTYPE
2848
100 CALL DGESL(WM(NPD),NEQ,NEQ,IWM(LIPVT),DELTA,0)
2851
C DUMMY SECTION FOR MTYPE=3
2856
400 MEBAND=2*IWM(LML)+IWM(LMU)+1
2857
CALL DGBSL(WM(NPD),MEBAND,NEQ,IWM(LML),
2858
* IWM(LMU),IWM(LIPVT),DELTA,0)
2860
C------END OF SUBROUTINE DDASLV------
2863
SUBROUTINE XERMSG (LIBRAR, SUBROU, MESSG, NERR, LEVEL)
2864
C***BEGIN PROLOGUE XERMSG
2865
C***PURPOSE Processes error messages for SLATEC and other libraries
2869
C***KEYWORDS ERROR MESSAGE, XERROR
2870
C***AUTHOR FONG, KIRBY, (NMFECC AT LLNL)
2872
C FRITSCH, F. N., (LLNL)
2875
C XERMSG processes a diagnostic message in a manner determined by the
2876
C value of LEVEL and the current value of the library error control
2877
C flag, KONTRL. See subroutine XSETF for details.
2878
C (XSETF is inoperable in this version.).
2880
C LIBRAR A character constant (or character variable) with the name
2881
C of the library. This will be 'SLATEC' for the SLATEC
2882
C Common Math Library. The error handling package is
2883
C general enough to be used by many libraries
2884
C simultaneously, so it is desirable for the routine that
2885
C detects and reports an error to identify the library name
2886
C as well as the routine name.
2888
C SUBROU A character constant (or character variable) with the name
2889
C of the routine that detected the error. Usually it is the
2890
C name of the routine that is calling XERMSG. There are
2891
C some instances where a user callable library routine calls
2892
C lower level subsidiary routines where the error is
2893
C detected. In such cases it may be more informative to
2894
C supply the name of the routine the user called rather than
2895
C the name of the subsidiary routine that detected the
2898
C MESSG A character constant (or character variable) with the text
2899
C of the error or warning message. In the example below,
2900
C the message is a character constant that contains a
2903
C CALL XERMSG ('SLATEC', 'MMPY',
2904
C *'THE ORDER OF THE MATRIX EXCEEDS THE ROW DIMENSION',
2907
C It is possible (and is sometimes desirable) to generate a
2908
C specific message--e.g., one that contains actual numeric
2909
C values. Specific numeric values can be converted into
2910
C character strings using formatted WRITE statements into
2911
C character variables. This is called standard Fortran
2912
C internal file I/O and is exemplified in the first three
2913
C lines of the following example. You can also catenate
2914
C substrings of characters to construct the error message.
2915
C Here is an example showing the use of both writing to
2916
C an internal file and catenating character strings.
2918
C CHARACTER*5 CHARN, CHARL
2919
C WRITE (CHARN,10) N
2920
C WRITE (CHARL,10) LDA
2922
C CALL XERMSG ('SLATEC', 'MMPY', 'THE ORDER'//CHARN//
2923
C * ' OF THE MATRIX EXCEEDS ITS ROW DIMENSION OF'//
2926
C There are two subtleties worth mentioning. One is that
2927
C the // for character catenation is used to construct the
2928
C error message so that no single character constant is
2929
C continued to the next line. This avoids confusion as to
2930
C whether there are trailing blanks at the end of the line.
2931
C The second is that by catenating the parts of the message
2932
C as an actual argument rather than encoding the entire
2933
C message into one large character variable, we avoid
2934
C having to know how long the message will be in order to
2935
C declare an adequate length for that large character
2936
C variable. XERMSG calls XERPRN to print the message using
2937
C multiple lines if necessary. If the message is very long,
2938
C XERPRN will break it into pieces of 72 characters (as
2939
C requested by XERMSG) for printing on multiple lines.
2940
C Also, XERMSG asks XERPRN to prefix each line with ' * '
2941
C so that the total line length could be 76 characters.
2942
C Note also that XERPRN scans the error message backwards
2943
C to ignore trailing blanks. Another feature is that
2944
C the substring '$$' is treated as a new line sentinel
2945
C by XERPRN. If you want to construct a multiline
2946
C message without having to count out multiples of 72
2947
C characters, just use '$$' as a separator. '$$'
2948
C obviously must occur within 72 characters of the
2949
C start of each line to have its intended effect since
2950
C XERPRN is asked to wrap around at 72 characters in
2951
C addition to looking for '$$'.
2953
C NERR An integer value that is chosen by the library routine's
2954
C author. It must be in the range -9999999 to 99999999 (8
2955
C printable digits). Each distinct error should have its
2956
C own error number. These error numbers should be described
2957
C in the machine readable documentation for the routine.
2958
C The error numbers need be unique only within each routine,
2959
C so it is reasonable for each routine to start enumerating
2960
C errors from 1 and proceeding to the next integer.
2962
C LEVEL An integer value in the range 0 to 2 that indicates the
2963
C level (severity) of the error. Their meanings are
2965
C -1 A warning message. This is used if it is not clear
2966
C that there really is an error, but the user's attention
2967
C may be needed. An attempt is made to only print this
2970
C 0 A warning message. This is used if it is not clear
2971
C that there really is an error, but the user's attention
2974
C 1 A recoverable error. This is used even if the error is
2975
C so serious that the routine cannot return any useful
2976
C answer. If the user has told the error package to
2977
C return after recoverable errors, then XERMSG will
2978
C return to the Library routine which can then return to
2979
C the user's routine. The user may also permit the error
2980
C package to terminate the program upon encountering a
2981
C recoverable error.
2983
C 2 A fatal error. XERMSG will not return to its caller
2984
C after it receives a fatal error. This level should
2985
C hardly ever be used; it is much better to allow the
2986
C user a chance to recover. An example of one of the few
2987
C cases in which it is permissible to declare a level 2
2988
C error is a reverse communication Library routine that
2989
C is likely to be called repeatedly until it integrates
2990
C across some interval. If there is a serious error in
2991
C the input such that another step cannot be taken and
2992
C the Library routine is called again without the input
2993
C error having been corrected by the caller, the Library
2994
C routine will probably be called forever with improper
2995
C input. In this case, it is reasonable to declare the
2996
C error to be fatal.
2998
C Each of the arguments to XERMSG is input; none will be modified by
2999
C XERMSG. A routine may make multiple calls to XERMSG with warning
3000
C level messages; however, after a call to XERMSG with a recoverable
3001
C error, the routine should return to the user.
3003
C***REFERENCES JONES, RONDALL E. AND KAHANER, DAVID K., "XERROR, THE
3004
C SLATEC ERROR-HANDLING PACKAGE", SOFTWARE - PRACTICE
3005
C AND EXPERIENCE, VOLUME 13, NO. 3, PP. 251-257,
3007
C***ROUTINES CALLED XERHLT, XERPRN
3008
C***REVISION HISTORY (YYMMDD)
3009
C 880101 DATE WRITTEN
3010
C 880621 REVISED AS DIRECTED AT SLATEC CML MEETING OF FEBRUARY 1988.
3011
C THERE ARE TWO BASIC CHANGES.
3012
C 1. A NEW ROUTINE, XERPRN, IS USED INSTEAD OF XERPRT TO
3013
C PRINT MESSAGES. THIS ROUTINE WILL BREAK LONG MESSAGES
3014
C INTO PIECES FOR PRINTING ON MULTIPLE LINES. '$$' IS
3015
C ACCEPTED AS A NEW LINE SENTINEL. A PREFIX CAN BE
3016
C ADDED TO EACH LINE TO BE PRINTED. XERMSG USES EITHER
3017
C ' ***' OR ' * ' AND LONG MESSAGES ARE BROKEN EVERY
3018
C 72 CHARACTERS (AT MOST) SO THAT THE MAXIMUM LINE
3019
C LENGTH OUTPUT CAN NOW BE AS GREAT AS 76.
3020
C 2. THE TEXT OF ALL MESSAGES IS NOW IN UPPER CASE SINCE THE
3021
C FORTRAN STANDARD DOCUMENT DOES NOT ADMIT THE EXISTENCE
3023
C 880708 REVISED AFTER THE SLATEC CML MEETING OF JUNE 29 AND 30.
3024
C THE PRINCIPAL CHANGES ARE
3025
C 1. CLARIFY COMMENTS IN THE PROLOGUES
3026
C 2. RENAME XRPRNT TO XERPRN
3027
C 3. REWORK HANDLING OF '$$' IN XERPRN TO HANDLE BLANK LINES
3028
C SIMILAR TO THE WAY FORMAT STATEMENTS HANDLE THE /
3029
C CHARACTER FOR NEW RECORDS.
3030
C 890706 REVISED WITH THE HELP OF FRED FRITSCH AND REG CLEMENS TO
3031
C CLEAN UP THE CODING.
3032
C 890721 REVISED TO USE NEW FEATURE IN XERPRN TO COUNT CHARACTERS IN
3034
C 891013 REVISED TO CORRECT COMMENTS.
3035
C 891214 Prologue converted to Version 4.0 format. (WRB)
3036
C 900510 Changed test on NERR to be -9999999 < NERR < 99999999, but
3037
C NERR .ne. 0, and on LEVEL to be -2 < LEVEL < 3. Added
3038
C LEVEL=-1 logic, changed calls to XERSAV to XERSVE, and
3039
C XERCTL to XERCNT. (RWC)
3040
C 901011 Removed error saving features to produce a simplified
3041
C version for distribution with DASSL and other LLNL codes.
3043
C***END PROLOGUE XERMSG
3044
CHARACTER*(*) LIBRAR, SUBROU, MESSG
3046
C***FIRST EXECUTABLE STATEMENT XERMSG
3048
C WE PRINT A FATAL ERROR MESSAGE AND TERMINATE FOR AN ERROR IN
3049
C CALLING XERMSG. THE ERROR NUMBER SHOULD BE POSITIVE,
3050
C AND THE LEVEL SHOULD BE BETWEEN 0 AND 2.
3052
IF (NERR.LT.-9999999 .OR. NERR.GT.99999999 .OR. NERR.EQ.0 .OR.
3053
* LEVEL.LT.-1 .OR. LEVEL.GT.2) THEN
3054
CALL XERPRN (' ***', -1, 'FATAL ERROR IN...$$ ' //
3055
* 'XERMSG -- INVALID ERROR NUMBER OR LEVEL$$ '//
3056
* 'JOB ABORT DUE TO FATAL ERROR.', 72)
3057
CALL XERHLT (' ***XERMSG -- INVALID INPUT')
3061
C SET DEFAULT VALUES FOR CONTROL PARAMETERS.
3066
C ANNOUNCE THE NAMES OF THE LIBRARY AND SUBROUTINE BY BUILDING A
3067
C MESSAGE IN CHARACTER VARIABLE TEMP (NOT EXCEEDING 66 CHARACTERS)
3068
C AND SENDING IT OUT VIA XERPRN. PRINT ONLY IF CONTROL FLAG
3071
IF (LKNTRL .NE. 0) THEN
3072
TEMP(1:21) = 'MESSAGE FROM ROUTINE '
3073
I = MIN(LEN(SUBROU), 16)
3074
TEMP(22:21+I) = SUBROU(1:I)
3075
TEMP(22+I:33+I) = ' IN LIBRARY '
3077
I = MIN(LEN(LIBRAR), 16)
3078
TEMP(LTEMP+1:LTEMP+I) = LIBRAR (1:I)
3079
TEMP(LTEMP+I+1:LTEMP+I+1) = '.'
3080
LTEMP = LTEMP + I + 1
3081
CALL XERPRN (' ***', -1, TEMP(1:LTEMP), 72)
3084
C IF LKNTRL IS POSITIVE, PRINT AN INTRODUCTORY LINE BEFORE
3085
C PRINTING THE MESSAGE. THE INTRODUCTORY LINE TELLS THE CHOICE
3086
C FROM EACH OF THE FOLLOWING TWO OPTIONS.
3087
C 1. LEVEL OF THE MESSAGE
3088
C 'INFORMATIVE MESSAGE'
3089
C 'POTENTIALLY RECOVERABLE ERROR'
3091
C 2. WHETHER CONTROL FLAG WILL ALLOW PROGRAM TO CONTINUE
3092
C 'PROGRAM CONTINUES'
3094
C NOTICE THAT THE LINE INCLUDING FOUR PREFIX CHARACTERS WILL NOT
3095
C EXCEED 74 CHARACTERS.
3096
C WE SKIP THE NEXT BLOCK IF THE INTRODUCTORY LINE IS NOT NEEDED.
3098
IF (LKNTRL .GT. 0) THEN
3100
C THE FIRST PART OF THE MESSAGE TELLS ABOUT THE LEVEL.
3102
IF (LEVEL .LE. 0) THEN
3103
TEMP(1:20) = 'INFORMATIVE MESSAGE,'
3105
ELSEIF (LEVEL .EQ. 1) THEN
3106
TEMP(1:30) = 'POTENTIALLY RECOVERABLE ERROR,'
3109
TEMP(1:12) = 'FATAL ERROR,'
3113
C THEN WHETHER THE PROGRAM WILL CONTINUE.
3115
IF ((MKNTRL.EQ.2 .AND. LEVEL.GE.1) .OR.
3116
* (MKNTRL.EQ.1 .AND. LEVEL.EQ.2)) THEN
3117
TEMP(LTEMP+1:LTEMP+17) = ' PROGRAM ABORTED.'
3120
TEMP(LTEMP+1:LTEMP+19) = ' PROGRAM CONTINUES.'
3124
CALL XERPRN (' ***', -1, TEMP(1:LTEMP), 72)
3127
C NOW SEND OUT THE MESSAGE.
3129
CALL XERPRN (' * ', -1, MESSG, 72)
3131
C IF LKNTRL IS POSITIVE, WRITE THE ERROR NUMBER.
3133
IF (LKNTRL .GT. 0) THEN
3134
WRITE (TEMP, '(''ERROR NUMBER = '', I8)') NERR
3136
IF (TEMP(I:I) .NE. ' ') GO TO 20
3139
20 CALL XERPRN (' * ', -1, TEMP(1:15) // TEMP(I:23), 72)
3142
C IF LKNTRL IS NOT ZERO, PRINT A BLANK LINE AND AN END OF MESSAGE.
3144
IF (LKNTRL .NE. 0) THEN
3145
CALL XERPRN (' * ', -1, ' ', 72)
3146
CALL XERPRN (' ***', -1, 'END OF MESSAGE', 72)
3147
CALL XERPRN (' ', 0, ' ', 72)
3150
C IF THE ERROR IS NOT FATAL OR THE ERROR IS RECOVERABLE AND THE
3151
C CONTROL FLAG IS SET FOR RECOVERY, THEN RETURN.
3153
30 IF (LEVEL.LE.0 .OR. (LEVEL.EQ.1 .AND. MKNTRL.LE.1)) RETURN
3155
C THE PROGRAM WILL BE STOPPED DUE TO AN UNRECOVERED ERROR OR A
3156
C FATAL ERROR. PRINT THE REASON FOR THE ABORT AND THE ERROR
3157
C SUMMARY IF THE CONTROL FLAG AND THE MAXIMUM ERROR COUNT PERMIT.
3159
IF (LKNTRL.GT.0) THEN
3160
IF (LEVEL .EQ. 1) THEN
3162
* (' ***', -1, 'JOB ABORT DUE TO UNRECOVERED ERROR.', 72)
3164
CALL XERPRN(' ***', -1, 'JOB ABORT DUE TO FATAL ERROR.', 72)
3170
SUBROUTINE XERHLT (MESSG)
3171
C***BEGIN PROLOGUE XERHLT
3173
C***PURPOSE Abort program execution and print error message.
3174
C***LIBRARY SLATEC (XERROR)
3176
C***TYPE ALL (XERHLT-A)
3177
C***KEYWORDS ERROR, XERROR
3178
C***AUTHOR JONES, R. E., (SNLA)
3182
C ***Note*** machine dependent routine
3183
C XERHLT aborts the execution of the program.
3184
C The error message causing the abort is given in the calling
3185
C sequence, in case one needs it for printing on a dayfile,
3188
C Description of Parameters
3189
C MESSG is as in XERROR.
3191
C***REFERENCES JONES R.E., KAHANER D.K., 'XERROR, THE SLATEC ERROR-
3192
C HANDLING PACKAGE', SAND82-0800, SANDIA LABORATORIES,
3194
C***ROUTINES CALLED (NONE)
3195
C***REVISION HISTORY (YYMMDD)
3196
C 790801 DATE WRITTEN as XERABT
3197
C 861211 REVISION DATE from Version 3.2
3198
C 891214 Prologue converted to Version 4.0 format. (BAB)
3199
C 900206 Routine changed from user-callable to subsidiary. (WRB)
3200
C 900510 Changed calling sequence to delete length of char string
3201
C Changed subroutine name from XERABT to XERHLT. (RWC)
3202
C***END PROLOGUE XERHLT
3204
C***FIRST EXECUTABLE STATEMENT XERHLT
3208
SUBROUTINE XERPRN (PREFIX, NPREF, MESSG, NWRAP)
3209
C***BEGIN PROLOGUE XERPRN
3211
C***PURPOSE This routine is called by XERMSG to print error messages
3215
C***KEYWORDS ERROR MESSAGES, PRINTING, XERROR
3216
C***AUTHOR FONG, KIRBY, (NMFECC AT LLNL)
3219
C This routine sends one or more lines to each of the (up to five)
3220
C logical units to which error messages are to be sent. This routine
3221
C is called several times by XERMSG, sometimes with a single line to
3222
C print and sometimes with a (potentially very long) message that may
3223
C wrap around into multiple lines.
3225
C PREFIX Input argument of type CHARACTER. This argument contains
3226
C characters to be put at the beginning of each line before
3227
C the body of the message. No more than 16 characters of
3228
C PREFIX will be used.
3230
C NPREF Input argument of type INTEGER. This argument is the number
3231
C of characters to use from PREFIX. If it is negative, the
3232
C intrinsic function LEN is used to determine its length. If
3233
C it is zero, PREFIX is not used. If it exceeds 16 or if
3234
C LEN(PREFIX) exceeds 16, only the first 16 characters will be
3235
C used. If NPREF is positive and the length of PREFIX is less
3236
C than NPREF, a copy of PREFIX extended with blanks to length
3237
C NPREF will be used.
3239
C MESSG Input argument of type CHARACTER. This is the text of a
3240
C message to be printed. If it is a long message, it will be
3241
C broken into pieces for printing on multiple lines. Each line
3242
C will start with the appropriate prefix and be followed by a
3243
C piece of the message. NWRAP is the number of characters per
3244
C piece; that is, after each NWRAP characters, we break and
3245
C start a new line. In addition the characters '$$' embedded
3246
C in MESSG are a sentinel for a new line. The counting of
3247
C characters up to NWRAP starts over for each new line. The
3248
C value of NWRAP typically used by XERMSG is 72 since many
3249
C older error messages in the SLATEC Library are laid out to
3250
C rely on wrap-around every 72 characters.
3252
C NWRAP Input argument of type INTEGER. This gives the maximum size
3253
C piece into which to break MESSG for printing on multiple
3254
C lines. An embedded '$$' ends a line, and the count restarts
3255
C at the following character. If a line break does not occur
3256
C on a blank (it would split a word) that word is moved to the
3257
C next line. Values of NWRAP less than 16 will be treated as
3258
C 16. Values of NWRAP greater than 132 will be treated as 132.
3259
C The actual line length will be NPREF + NWRAP after NPREF has
3260
C been adjusted to fall between 0 and 16 and NWRAP has been
3261
C adjusted to fall between 16 and 132.
3263
C***REFERENCES (NONE)
3264
C***ROUTINES CALLED I1MACH, XGETUA
3265
C***REVISION HISTORY (YYMMDD)
3266
C 880621 DATE WRITTEN
3267
C 880708 REVISED AFTER THE SLATEC CML SUBCOMMITTEE MEETING OF
3268
C JUNE 29 AND 30 TO CHANGE THE NAME TO XERPRN AND TO REWORK
3269
C THE HANDLING OF THE NEW LINE SENTINEL TO BEHAVE LIKE THE
3270
C SLASH CHARACTER IN FORMAT STATEMENTS.
3271
C 890706 REVISED WITH THE HELP OF FRED FRITSCH AND REG CLEMMENS TO
3272
C STREAMLINE THE CODING AND FIX A BUG THAT CAUSED EXTRA BLANK
3273
C LINES TO BE PRINTED.
3274
C 890721 REVISED TO ADD A NEW FEATURE. A NEGATIVE VALUE OF NPREF
3275
C CAUSES LEN(PREFIX) TO BE USED AS THE LENGTH.
3276
C 891013 REVISED TO CORRECT ERROR IN CALCULATING PREFIX LENGTH.
3277
C 891214 Prologue converted to Version 4.0 format. (WRB)
3278
C 900510 Added code to break messages between words. (RWC)
3279
C***END PROLOGUE XERPRN
3280
CHARACTER*(*) PREFIX, MESSG
3281
INTEGER NPREF, NWRAP
3283
INTEGER IU(5), NUNIT
3285
PARAMETER (NEWLIN = '$$')
3286
C***FIRST EXECUTABLE STATEMENT XERPRN
3287
CALL XGETUA(IU,NUNIT)
3289
C A ZERO VALUE FOR A LOGICAL UNIT NUMBER MEANS TO USE THE STANDARD
3290
C ERROR MESSAGE UNIT INSTEAD. I1MACH(4) RETRIEVES THE STANDARD
3291
C ERROR MESSAGE UNIT.
3295
IF (IU(I) .EQ. 0) IU(I) = N
3298
C LPREF IS THE LENGTH OF THE PREFIX. THE PREFIX IS PLACED AT THE
3299
C BEGINNING OF CBUFF, THE CHARACTER BUFFER, AND KEPT THERE DURING
3300
C THE REST OF THIS ROUTINE.
3302
IF ( NPREF .LT. 0 ) THEN
3307
LPREF = MIN(16, LPREF)
3308
IF (LPREF .NE. 0) CBUFF(1:LPREF) = PREFIX
3310
C LWRAP IS THE MAXIMUM NUMBER OF CHARACTERS WE WANT TO TAKE AT ONE
3311
C TIME FROM MESSG TO PRINT ON ONE LINE.
3313
LWRAP = MAX(16, MIN(132, NWRAP))
3315
C SET LENMSG TO THE LENGTH OF MESSG, IGNORE ANY TRAILING BLANKS.
3320
IF (MESSG(LENMSG:LENMSG) .NE. ' ') GO TO 30
3325
C IF THE MESSAGE IS ALL BLANKS, THEN PRINT ONE BLANK LINE.
3327
IF (LENMSG .EQ. 0) THEN
3328
CBUFF(LPREF+1:LPREF+1) = ' '
3330
WRITE(IU(I), '(A)') CBUFF(1:LPREF+1)
3335
C SET NEXTC TO THE POSITION IN MESSG WHERE THE NEXT SUBSTRING
3336
C STARTS. FROM THIS POSITION WE SCAN FOR THE NEW LINE SENTINEL.
3337
C WHEN NEXTC EXCEEDS LENMSG, THERE IS NO MORE TO PRINT.
3338
C WE LOOP BACK TO LABEL 50 UNTIL ALL PIECES HAVE BEEN PRINTED.
3340
C WE LOOK FOR THE NEXT OCCURRENCE OF THE NEW LINE SENTINEL. THE
3341
C INDEX INTRINSIC FUNCTION RETURNS ZERO IF THERE IS NO OCCURRENCE
3342
C OR IF THE LENGTH OF THE FIRST ARGUMENT IS LESS THAN THE LENGTH
3343
C OF THE SECOND ARGUMENT.
3345
C THERE ARE SEVERAL CASES WHICH SHOULD BE CHECKED FOR IN THE
3346
C FOLLOWING ORDER. WE ARE ATTEMPTING TO SET LPIECE TO THE NUMBER
3347
C OF CHARACTERS THAT SHOULD BE TAKEN FROM MESSG STARTING AT
3350
C LPIECE .EQ. 0 THE NEW LINE SENTINEL DOES NOT OCCUR IN THE
3351
C REMAINDER OF THE CHARACTER STRING. LPIECE
3352
C SHOULD BE SET TO LWRAP OR LENMSG+1-NEXTC,
3353
C WHICHEVER IS LESS.
3355
C LPIECE .EQ. 1 THE NEW LINE SENTINEL STARTS AT MESSG(NEXTC:
3356
C NEXTC). LPIECE IS EFFECTIVELY ZERO, AND WE
3357
C PRINT NOTHING TO AVOID PRODUCING UNNECESSARY
3358
C BLANK LINES. THIS TAKES CARE OF THE SITUATION
3359
C WHERE THE LIBRARY ROUTINE HAS A MESSAGE OF
3360
C EXACTLY 72 CHARACTERS FOLLOWED BY A NEW LINE
3361
C SENTINEL FOLLOWED BY MORE CHARACTERS. NEXTC
3362
C SHOULD BE INCREMENTED BY 2.
3364
C LPIECE .GT. LWRAP+1 REDUCE LPIECE TO LWRAP.
3366
C ELSE THIS LAST CASE MEANS 2 .LE. LPIECE .LE. LWRAP+1
3367
C RESET LPIECE = LPIECE-1. NOTE THAT THIS
3368
C PROPERLY HANDLES THE END CASE WHERE LPIECE .EQ.
3369
C LWRAP+1. THAT IS, THE SENTINEL FALLS EXACTLY
3370
C AT THE END OF A LINE.
3373
50 LPIECE = INDEX(MESSG(NEXTC:LENMSG), NEWLIN)
3374
IF (LPIECE .EQ. 0) THEN
3376
C THERE WAS NO NEW LINE SENTINEL FOUND.
3379
LPIECE = MIN(LWRAP, LENMSG+1-NEXTC)
3380
IF (LPIECE .LT. LENMSG+1-NEXTC) THEN
3381
DO 52 I=LPIECE+1,2,-1
3382
IF (MESSG(NEXTC+I-1:NEXTC+I-1) .EQ. ' ') THEN
3389
54 CBUFF(LPREF+1:LPREF+LPIECE) = MESSG(NEXTC:NEXTC+LPIECE-1)
3390
NEXTC = NEXTC + LPIECE + IDELTA
3391
ELSEIF (LPIECE .EQ. 1) THEN
3393
C WE HAVE A NEW LINE SENTINEL AT MESSG(NEXTC:NEXTC+1).
3394
C DON'T PRINT A BLANK LINE.
3398
ELSEIF (LPIECE .GT. LWRAP+1) THEN
3400
C LPIECE SHOULD BE SET DOWN TO LWRAP.
3404
DO 56 I=LPIECE+1,2,-1
3405
IF (MESSG(NEXTC+I-1:NEXTC+I-1) .EQ. ' ') THEN
3411
58 CBUFF(LPREF+1:LPREF+LPIECE) = MESSG(NEXTC:NEXTC+LPIECE-1)
3412
NEXTC = NEXTC + LPIECE + IDELTA
3415
C IF WE ARRIVE HERE, IT MEANS 2 .LE. LPIECE .LE. LWRAP+1.
3416
C WE SHOULD DECREMENT LPIECE BY ONE.
3419
CBUFF(LPREF+1:LPREF+LPIECE) = MESSG(NEXTC:NEXTC+LPIECE-1)
3420
NEXTC = NEXTC + LPIECE + 2
3426
WRITE(IU(I), '(A)') CBUFF(1:LPREF+LPIECE)
3429
IF (NEXTC .LE. LENMSG) GO TO 50
3433
SUBROUTINE XGETUA (IUNITA, N)
3434
C***BEGIN PROLOGUE XGETUA
3435
C***PURPOSE Return unit number(s) to which error messages are being
3437
C***LIBRARY SLATEC (XERROR)
3439
C***TYPE ALL (XGETUA-A)
3440
C***KEYWORDS ERROR, XERROR
3441
C***AUTHOR JONES, R. E., (SNLA)
3443
C FRITSCH, F. N., (LLNL)
3447
C XGETUA may be called to determine the unit number or numbers
3448
C to which error messages are being sent.
3449
C These unit numbers may have been set by a call to XSETUN,
3450
C or a call to XSETUA, or may be a default value.
3452
C Description of Parameters
3454
C IUNIT - an array of one to five unit numbers, depending
3455
C on the value of N. A value of zero refers to the
3456
C default unit, as defined by the I1MACH machine
3457
C constant routine. Only IUNIT(1),...,IUNIT(N) are
3458
C defined by XGETUA. The values of IUNIT(N+1),...,
3459
C IUNIT(5) are not defined (for N .LT. 5) or altered
3460
C in any way by XGETUA.
3461
C N - the number of units to which copies of the
3462
C error messages are being sent. N will be in the
3463
C range from 1 to 5.
3465
C CAUTION: The use of COMMON in this version is not safe for
3468
C***REFERENCES JONES R.E., KAHANER D.K., 'XERROR, THE SLATEC ERROR-
3469
C HANDLING PACKAGE', SAND82-0800, SANDIA LABORATORIES,
3471
C***ROUTINES CALLED (NONE)
3472
C***COMMON BLOCKS XERUNI
3473
C***REVISION HISTORY (YYMMDD)
3474
C 790801 DATE WRITTEN
3475
C 861211 REVISION DATE from Version 3.2
3476
C 891214 Prologue converted to Version 4.0 format. (BAB)
3477
C 901011 Rewritten to not use J4SAVE. (FNF)
3478
C 901012 Corrected initialization problem. (FNF)
3479
C***END PROLOGUE XGETUA
3481
INTEGER NUNIT, IUNIT(5)
3482
COMMON /XERUNI/ NUNIT, IUNIT
3483
C***FIRST EXECUTABLE STATEMENT XGETUA
3484
C Initialize so XERMSG will use standard error unit number if
3485
C block has not been set up by a CALL XSETUA.
3486
C CAUTION: This assumes uninitialized COMMON tests .LE.0 .
3487
IF (NUNIT.LE.0) THEN
3493
IUNITA(I) = IUNIT(I)
3498
SUBROUTINE XSETUA (IUNITA, N)
3499
C***BEGIN PROLOGUE XSETUA
3500
C***PURPOSE Set logical unit numbers (up to 5) to which error
3501
C messages are to be sent.
3502
C***LIBRARY SLATEC (XERROR)
3504
C***TYPE ALL (XSETUA-A)
3505
C***KEYWORDS ERROR, XERROR
3506
C***AUTHOR JONES, R. E., (SNLA)
3508
C FRITSCH, F. N., (LLNL)
3512
C XSETUA may be called to declare a list of up to five
3513
C logical units, each of which is to receive a copy of
3514
C each error message processed by this package.
3515
C The purpose of XSETUA is to allow simultaneous printing
3516
C of each error message on, say, a main output file,
3517
C an interactive terminal, and other files such as graphics
3518
C communication files.
3520
C Description of Parameters
3522
C IUNIT - an array of up to five unit numbers.
3523
C Normally these numbers should all be different
3524
C (but duplicates are not prohibited.)
3525
C N - the number of unit numbers provided in IUNIT
3526
C must have 1 .LE. N .LE. 5.
3528
C CAUTION: The use of COMMON in this version is not safe for
3531
C***REFERENCES JONES R.E., KAHANER D.K., 'XERROR, THE SLATEC ERROR-
3532
C HANDLING PACKAGE', SAND82-0800, SANDIA LABORATORIES,
3534
C***ROUTINES CALLED XERMSG
3535
C***COMMON BLOCKS XERUNI
3536
C***REVISION HISTORY (YYMMDD)
3537
C 790801 DATE WRITTEN
3538
C 861211 REVISION DATE from Version 3.2
3539
C 891214 Prologue converted to Version 4.0 format. (BAB)
3540
C 900510 Change call to XERRWV to XERMSG. (RWC)
3541
C 901011 Rewritten to not use J4SAVE. (FNF)
3542
C***END PROLOGUE XSETUA
3544
INTEGER NUNIT, IUNIT(5)
3545
COMMON /XERUNI/ NUNIT, IUNIT
3547
C***FIRST EXECUTABLE STATEMENT XSETUA
3549
IF (N.LT.1 .OR. N.GT.5) THEN
3550
WRITE (XERN1, '(I8)') N
3551
CALL XERMSG ('SLATEC', 'XSETUA',
3552
* 'INVALID NUMBER OF UNITS, N = ' // XERN1, 1, 2)
3557
IUNIT(I) = IUNITA(I)
added
removed
Lib/integrate/odepack/ddassl.f
