2
SUBROUTINE DSTOD (NEQ, Y, YH, NYH, YH1, EWT, SAVF, ACOR, WM, IWM,
3
+ DF, DJAC, RPAR, IPAR)
4
C***BEGIN PROLOGUE DSTOD
6
C***PURPOSE Subsidiary to DDEBDF
8
C***TYPE DOUBLE PRECISION (STOD-S, DSTOD-D)
9
C***AUTHOR Watts, H. A., (SNLA)
12
C DSTOD integrates a system of first order odes over one step in the
13
C integrator package DDEBDF.
14
C ----------------------------------------------------------------------
15
C DSTOD performs one step of the integration of an initial value
16
C problem for a system of ordinary differential equations.
17
C Note.. DSTOD is independent of the value of the iteration method
18
C indicator MITER, when this is .NE. 0, and hence is independent
19
C of the type of chord method used, or the Jacobian structure.
20
C Communication with DSTOD is done with the following variables..
22
C Y = An array of length .GE. N used as the Y argument in
23
C all calls to DF and DJAC.
24
C NEQ = Integer array containing problem size in NEQ(1), and
25
C passed as the NEQ argument in all calls to DF and DJAC.
26
C YH = An NYH by LMAX array containing the dependent variables
27
C and their approximate scaled derivatives, where
28
C LMAX = MAXORD + 1. YH(I,J+1) contains the approximate
29
C J-th derivative of Y(I), scaled by H**J/FACTORIAL(J)
30
C (J = 0,1,...,NQ). On entry for the first step, the first
31
C two columns of YH must be set from the initial values.
32
C NYH = A constant integer .GE. N, the first dimension of YH.
33
C YH1 = A one-dimensional array occupying the same space as YH.
34
C EWT = An array of N elements with which the estimated local
35
C errors in YH are compared.
36
C SAVF = An array of working storage, of length N.
37
C ACOR = A work array of length N, used for the accumulated
38
C corrections. On a successful return, ACOR(I) contains
39
C the estimated one-step local error in Y(I).
40
C WM,IWM = DOUBLE PRECISION and INTEGER work arrays associated with
41
C matrix operations in chord iteration (MITER .NE. 0).
42
C DPJAC = Name of routine to evaluate and preprocess Jacobian matrix
43
C if a chord method is being used.
44
C DSLVS = Name of routine to solve linear system in chord iteration.
45
C H = The step size to be attempted on the next step.
46
C H is altered by the error control algorithm during the
47
C problem. H can be either positive or negative, but its
48
C sign must remain constant throughout the problem.
49
C HMIN = The minimum absolute value of the step size H to be used.
50
C HMXI = Inverse of the maximum absolute value of H to be used.
51
C HMXI = 0.0 is allowed and corresponds to an infinite HMAX.
52
C HMIN and HMXI may be changed at any time, but will not
53
C take effect until the next change of H is considered.
54
C TN = The independent variable. TN is updated on each step taken.
55
C JSTART = An integer used for input only, with the following
56
C values and meanings..
57
C 0 Perform the first step.
58
C .GT.0 Take a new step continuing from the last.
59
C -1 Take the next step with a new value of H, MAXORD,
60
C N, METH, MITER, and/or matrix parameters.
61
C -2 Take the next step with a new value of H,
62
C but with other inputs unchanged.
63
C On return, JSTART is set to 1 to facilitate continuation.
64
C KFLAG = a completion code with the following meanings..
65
C 0 The step was successful.
66
C -1 The requested error could not be achieved.
67
C -2 Corrector convergence could not be achieved.
68
C A return with KFLAG = -1 or -2 means either
69
C ABS(H) = HMIN or 10 consecutive failures occurred.
70
C On a return with KFLAG negative, the values of TN and
71
C the YH array are as of the beginning of the last
72
C step, and H is the last step size attempted.
73
C MAXORD = The maximum order of integration method to be allowed.
74
C METH/MITER = The method flags. See description in driver.
75
C N = The number of first-order differential equations.
76
C ----------------------------------------------------------------------
79
C***ROUTINES CALLED DCFOD, DPJAC, DSLVS, DVNRMS
80
C***COMMON BLOCKS DDEBD1
81
C***REVISION HISTORY (YYMMDD)
83
C 890531 Changed all specific intrinsics to generic. (WRB)
84
C 890911 Removed unnecessary intrinsics. (WRB)
85
C 891214 Prologue converted to Version 4.0 format. (BAB)
86
C 900328 Added TYPE section. (WRB)
87
C 910722 Updated AUTHOR section. (ALS)
88
C 920422 Changed DIMENSION statement. (WRB)
89
C***END PROLOGUE DSTOD
91
INTEGER I, I1, IALTH, IER, IOD, IOWND, IPAR, IPUP, IREDO, IRET,
92
1 IWM, J, JB, JSTART, KFLAG, KSTEPS, L, LMAX, M, MAXORD,
93
2 MEO, METH, MITER, N, NCF, NEQ, NEWQ, NFE, NJE, NQ, NQNYH,
94
3 NQU, NST, NSTEPJ, NYH
95
DOUBLE PRECISION ACOR, CONIT, CRATE, DCON, DDN,
96
1 DEL, DELP, DSM, DUP, DVNRMS, EL, EL0, ELCO,
97
2 EWT, EXDN, EXSM, EXUP, H, HMIN, HMXI, HOLD, HU, R, RC,
98
3 RH, RHDN, RHSM, RHUP, RMAX, ROWND, RPAR, SAVF, TESCO,
99
4 TN, TOLD, UROUND, WM, Y, YH, YH1
102
DIMENSION Y(*),YH(NYH,*),YH1(*),EWT(*),SAVF(*),ACOR(*),WM(*),
103
1 IWM(*),RPAR(*),IPAR(*)
104
COMMON /DDEBD1/ ROWND,CONIT,CRATE,EL(13),ELCO(13,12),HOLD,RC,RMAX,
105
1 TESCO(3,12),EL0,H,HMIN,HMXI,HU,TN,UROUND,IOWND(7),
106
2 KSTEPS,IOD(6),IALTH,IPUP,LMAX,MEO,NQNYH,NSTEPJ,
107
3 IER,JSTART,KFLAG,L,METH,MITER,MAXORD,N,NQ,NST,NFE,
111
C BEGIN BLOCK PERMITTING ...EXITS TO 690
112
C BEGIN BLOCK PERMITTING ...EXITS TO 60
113
C***FIRST EXECUTABLE STATEMENT DSTOD
117
IF (JSTART .GT. 0) GO TO 160
118
IF (JSTART .EQ. -1) GO TO 10
119
IF (JSTART .EQ. -2) GO TO 90
120
C ---------------------------------------------------------
121
C ON THE FIRST CALL, THE ORDER IS SET TO 1, AND OTHER
122
C VARIABLES ARE INITIALIZED. RMAX IS THE MAXIMUM RATIO BY
123
C WHICH H CAN BE INCREASED IN A SINGLE STEP. IT IS
124
C INITIALLY 1.E4 TO COMPENSATE FOR THE SMALL INITIAL H,
125
C BUT THEN IS NORMALLY EQUAL TO 10. IF A FAILURE OCCURS
126
C (IN CORRECTOR CONVERGENCE OR ERROR TEST), RMAX IS SET AT
127
C 2 FOR THE NEXT INCREASE.
128
C ---------------------------------------------------------
144
C BEGIN BLOCK PERMITTING ...EXITS TO 30
145
C ------------------------------------------------------
146
C THE FOLLOWING BLOCK HANDLES PRELIMINARIES NEEDED WHEN
147
C JSTART = -1. IPUP IS SET TO MITER TO FORCE A MATRIX
148
C UPDATE. IF AN ORDER INCREASE IS ABOUT TO BE
149
C CONSIDERED (IALTH = 1), IALTH IS RESET TO 2 TO
150
C POSTPONE CONSIDERATION ONE MORE STEP. IF THE CALLER
151
C HAS CHANGED METH, DCFOD IS CALLED TO RESET THE
152
C COEFFICIENTS OF THE METHOD. IF THE CALLER HAS
153
C CHANGED MAXORD TO A VALUE LESS THAN THE CURRENT
154
C ORDER NQ, NQ IS REDUCED TO MAXORD, AND A NEW H CHOSEN
155
C ACCORDINGLY. IF H IS TO BE CHANGED, YH MUST BE
156
C RESCALED. IF H OR METH IS BEING CHANGED, IALTH IS
157
C RESET TO L = NQ + 1 TO PREVENT FURTHER CHANGES IN H
158
C FOR THAT MANY STEPS.
159
C ------------------------------------------------------
162
IF (IALTH .EQ. 1) IALTH = 2
163
IF (METH .EQ. MEO) GO TO 20
164
CALL DCFOD(METH,ELCO,TESCO)
167
IF (NQ .GT. MAXORD) GO TO 30
173
IF (NQ .LE. MAXORD) GO TO 90
184
DDN = DVNRMS(N,SAVF,EWT)/TESCO(1,L)
186
RHDN = 1.0D0/(1.3D0*DDN**EXDN + 0.0000013D0)
189
IF (H .EQ. HOLD) GO TO 660
190
RH = MIN(RH,ABS(H/HOLD))
194
C ------------------------------------------------------------
195
C DCFOD IS CALLED TO GET ALL THE INTEGRATION COEFFICIENTS
196
C FOR THE CURRENT METH. THEN THE EL VECTOR AND RELATED
197
C CONSTANTS ARE RESET WHENEVER THE ORDER NQ IS CHANGED, OR AT
198
C THE START OF THE PROBLEM.
199
C ------------------------------------------------------------
200
CALL DCFOD(METH,ELCO,TESCO)
203
C BEGIN BLOCK PERMITTING ...EXITS TO 680
211
GO TO (90,660,160), IRET
212
C ---------------------------------------------------------
213
C IF H IS BEING CHANGED, THE H RATIO RH IS CHECKED AGAINST
214
C RMAX, HMIN, AND HMXI, AND THE YH ARRAY RESCALED. IALTH
215
C IS SET TO L = NQ + 1 TO PREVENT A CHANGE OF H FOR THAT
216
C MANY STEPS, UNLESS FORCED BY A CONVERGENCE OR ERROR TEST
218
C ---------------------------------------------------------
220
IF (H .EQ. HOLD) GO TO 160
227
RH = RH/MAX(1.0D0,ABS(H)*HMXI*RH)
238
IF (IREDO .NE. 0) GO TO 150
240
R = 1.0D0/TESCO(2,NQU)
244
C ...............EXIT
247
C ------------------------------------------------------
248
C THIS SECTION COMPUTES THE PREDICTED VALUES BY
249
C EFFECTIVELY MULTIPLYING THE YH ARRAY BY THE PASCAL
250
C TRIANGLE MATRIX. RC IS THE RATIO OF NEW TO OLD
251
C VALUES OF THE COEFFICIENT H*EL(1). WHEN RC DIFFERS
252
C FROM 1 BY MORE THAN 30 PERCENT, IPUP IS SET TO MITER
253
C TO FORCE DPJAC TO BE CALLED, IF A JACOBIAN IS
254
C INVOLVED. IN ANY CASE, DPJAC IS CALLED AT LEAST
256
C ------------------------------------------------------
259
C BEGIN BLOCK PERMITTING ...EXITS TO 610
260
C BEGIN BLOCK PERMITTING ...EXITS TO 490
261
IF (ABS(RC-1.0D0) .GT. 0.3D0) IPUP = MITER
262
IF (NST .GE. NSTEPJ + 20) IPUP = MITER
268
YH1(I) = YH1(I) + YH1(I+NYH)
272
C ---------------------------------------------
273
C UP TO 3 CORRECTOR ITERATIONS ARE TAKEN. A
274
C CONVERGENCE TEST IS MADE ON THE R.M.S. NORM
275
C OF EACH CORRECTION, WEIGHTED BY THE ERROR
276
C WEIGHT VECTOR EWT. THE SUM OF THE
277
C CORRECTIONS IS ACCUMULATED IN THE VECTOR
278
C ACOR(I). THE YH ARRAY IS NOT ALTERED IN THE
280
C ---------------------------------------------
286
CALL DF(TN,Y,SAVF,RPAR,IPAR)
288
IF (IPUP .LE. 0) GO TO 220
289
C ---------------------------------------
290
C IF INDICATED, THE MATRIX P = I -
291
C H*EL(1)*J IS REEVALUATED AND
292
C PREPROCESSED BEFORE STARTING THE
293
C CORRECTOR ITERATION. IPUP IS SET TO 0
294
C AS AN INDICATOR THAT THIS HAS BEEN
296
C ---------------------------------------
301
CALL DPJAC(NEQ,Y,YH,NYH,EWT,ACOR,SAVF,
302
1 WM,IWM,DF,DJAC,RPAR,IPAR)
304
IF (IER .NE. 0) GO TO 440
310
IF (MITER .NE. 0) GO TO 270
311
C ------------------------------------
312
C IN THE CASE OF FUNCTIONAL
313
C ITERATION, UPDATE Y DIRECTLY FROM
314
C THE RESULT OF THE LAST FUNCTION
316
C ------------------------------------
318
SAVF(I) = H*SAVF(I) - YH(I,2)
319
Y(I) = SAVF(I) - ACOR(I)
321
DEL = DVNRMS(N,Y,EWT)
323
Y(I) = YH(I,1) + EL(1)*SAVF(I)
328
C ------------------------------------
329
C IN THE CASE OF THE CHORD METHOD,
330
C COMPUTE THE CORRECTOR ERROR, AND
331
C SOLVE THE LINEAR SYSTEM WITH THAT
332
C AS RIGHT-HAND SIDE AND P AS
333
C COEFFICIENT MATRIX.
334
C ------------------------------------
337
1 - (YH(I,2) + ACOR(I))
339
CALL DSLVS(WM,IWM,Y,SAVF)
341
IF (IER .NE. 0) GO TO 430
342
DEL = DVNRMS(N,Y,EWT)
344
ACOR(I) = ACOR(I) + Y(I)
345
Y(I) = YH(I,1) + EL(1)*ACOR(I)
348
C ---------------------------------------
349
C TEST FOR CONVERGENCE. IF M.GT.0, AN
350
C ESTIMATE OF THE CONVERGENCE RATE
351
C CONSTANT IS STORED IN CRATE, AND THIS
352
C IS USED IN THE TEST.
353
C ---------------------------------------
355
1 CRATE = MAX(0.2D0*CRATE,DEL/DELP)
356
DCON = DEL*MIN(1.0D0,1.5D0*CRATE)
357
1 /(TESCO(2,NQ)*CONIT)
358
IF (DCON .GT. 1.0D0) GO TO 420
359
C ------------------------------------
360
C THE CORRECTOR HAS CONVERGED. IPUP
361
C IS SET TO -1 IF MITER .NE. 0, TO
362
C SIGNAL THAT THE JACOBIAN INVOLVED
363
C MAY NEED UPDATING LATER. THE LOCAL
364
C ERROR TEST IS MADE AND CONTROL
365
C PASSES TO STATEMENT 500 IF IT
367
C ------------------------------------
368
IF (MITER .NE. 0) IPUP = -1
369
IF (M .EQ. 0) DSM = DEL/TESCO(2,NQ)
371
1 DSM = DVNRMS(N,ACOR,EWT)
373
IF (DSM .GT. 1.0D0) GO TO 380
375
C PERMITTING ...EXITS TO 360
376
C ------------------------------
377
C AFTER A SUCCESSFUL STEP,
378
C UPDATE THE YH ARRAY.
379
C CONSIDER CHANGING H IF IALTH
380
C = 1. OTHERWISE DECREASE
381
C IALTH BY 1. IF IALTH IS THEN
382
C 1 AND NQ .LT. MAXORD, THEN
383
C ACOR IS SAVED FOR USE IN A
384
C POSSIBLE ORDER INCREASE ON
385
C THE NEXT STEP. IF A CHANGE
386
C IN H IS CONSIDERED, AN
387
C INCREASE OR DECREASE IN ORDER
388
C BY ONE IS CONSIDERED ALSO. A
389
C CHANGE IN H IS MADE ONLY IF
390
C IT IS BY A FACTOR OF AT LEAST
391
C 1.1. IF NOT, IALTH IS SET TO
392
C 3 TO PREVENT TESTING FOR THAT
394
C ------------------------------
408
IF (IALTH .NE. 0) GO TO 340
409
C ---------------------------
410
C REGARDLESS OF THE SUCCESS
411
C OR FAILURE OF THE STEP,
412
C FACTORS RHDN, RHSM, AND
413
C RHUP ARE COMPUTED, BY
415
C MULTIPLIED AT ORDER NQ -
416
C 1, ORDER NQ, OR ORDER NQ +
417
C 1, RESPECTIVELY. IN THE
418
C CASE OF FAILURE, RHUP =
419
C 0.0 TO AVOID AN ORDER
420
C INCREASE. THE LARGEST OF
421
C THESE IS DETERMINED AND
422
C THE NEW ORDER CHOSEN
423
C ACCORDINGLY. IF THE ORDER
424
C IS TO BE INCREASED, WE
425
C COMPUTE ONE ADDITIONAL
427
C ---------------------------
429
C .....................EXIT
430
IF (L .EQ. LMAX) GO TO 490
435
DUP = DVNRMS(N,SAVF,EWT)
441
C .....................EXIT
445
IF (IALTH .GT. 1) GO TO 360
447
IF (L .EQ. LMAX) GO TO 360
452
R = 1.0D0/TESCO(2,NQU)
456
C .................................EXIT
459
C ------------------------------------
460
C THE ERROR TEST FAILED. KFLAG KEEPS
461
C TRACK OF MULTIPLE FAILURES.
462
C RESTORE TN AND THE YH ARRAY TO
463
C THEIR PREVIOUS VALUES, AND PREPARE
464
C TO TRY THE STEP AGAIN. COMPUTE THE
465
C OPTIMUM STEP SIZE FOR THIS OR ONE
466
C LOWER ORDER. AFTER 2 OR MORE
467
C FAILURES, H IS FORCED TO DECREASE
468
C BY A FACTOR OF 0.2 OR LESS.
469
C ------------------------------------
476
YH1(I) = YH1(I) - YH1(I+NYH)
480
IF (ABS(H) .GT. HMIN*1.00001D0)
482
C ---------------------------------
483
C ALL RETURNS ARE MADE THROUGH
484
C THIS SECTION. H IS SAVED IN
485
C HOLD TO ALLOW THE CALLER TO
486
C CHANGE H ON THE NEXT STEP.
487
C ---------------------------------
489
C .................................EXIT
492
C ...............EXIT
493
IF (KFLAG .LE. -3) GO TO 610
501
IF (M .EQ. 3) GO TO 430
503
IF (M .GE. 2 .AND. DEL .GT. 2.0D0*DELP)
506
CALL DF(TN,Y,SAVF,RPAR,IPAR)
510
C ------------------------------------------
511
C THE CORRECTOR ITERATION FAILED TO
512
C CONVERGE IN 3 TRIES. IF MITER .NE. 0 AND
513
C THE JACOBIAN IS OUT OF DATE, DPJAC IS
514
C CALLED FOR THE NEXT TRY. OTHERWISE THE
515
C YH ARRAY IS RETRACTED TO ITS VALUES
516
C BEFORE PREDICTION, AND H IS REDUCED, IF
517
C POSSIBLE. IF H CANNOT BE REDUCED OR 10
518
C FAILURES HAVE OCCURRED, EXIT WITH KFLAG =
520
C ------------------------------------------
522
IF (IPUP .EQ. 0) GO TO 440
533
YH1(I) = YH1(I) - YH1(I+NYH)
536
IF (ABS(H) .GT. HMIN*1.00001D0) GO TO 470
538
C ........................EXIT
541
IF (NCF .NE. 10) GO TO 480
543
C ........................EXIT
553
RHSM = 1.0D0/(1.2D0*DSM**EXSM + 0.0000012D0)
555
IF (NQ .EQ. 1) GO TO 500
556
DDN = DVNRMS(N,YH(1,L),EWT)/TESCO(1,NQ)
558
RHDN = 1.0D0/(1.3D0*DDN**EXDN + 0.0000013D0)
560
IF (RHSM .GE. RHUP) GO TO 550
561
IF (RHUP .LE. RHDN) GO TO 540
564
IF (RH .GE. 1.1D0) GO TO 520
566
R = 1.0D0/TESCO(2,NQU)
570
C ...........................EXIT
575
YH(I,NEWQ+1) = ACOR(I)*R
580
C ..................EXIT
585
IF (RHSM .LT. RHDN) GO TO 580
588
IF (KFLAG .EQ. 0 .AND. RH .LT. 1.1D0)
590
IF (KFLAG .LE. -2) RH = MIN(RH,0.2D0)
591
C ------------------------------------------
592
C IF THERE IS A CHANGE OF ORDER, RESET NQ,
593
C L, AND THE COEFFICIENTS. IN ANY CASE H
594
C IS RESET ACCORDING TO RH AND THE YH ARRAY
595
C IS RESCALED. THEN EXIT FROM 680 IF THE
596
C STEP WAS OK, OR REDO THE STEP OTHERWISE.
597
C ------------------------------------------
599
IF (NEWQ .EQ. NQ) GO TO 650
603
C ..................EXIT
607
R = 1.0D0/TESCO(2,NQU)
611
C .....................EXIT
616
IF (KFLAG .LT. 0 .AND. RH .GT. 1.0D0) RH = 1.0D0
617
IF (KFLAG .EQ. 0 .AND. RH .LT. 1.1D0) GO TO 590
618
IF (KFLAG .LE. -2) RH = MIN(RH,0.2D0)
619
C ---------------------------------------------
620
C IF THERE IS A CHANGE OF ORDER, RESET NQ, L,
621
C AND THE COEFFICIENTS. IN ANY CASE H IS
622
C RESET ACCORDING TO RH AND THE YH ARRAY IS
623
C RESCALED. THEN EXIT FROM 680 IF THE STEP
624
C WAS OK, OR REDO THE STEP OTHERWISE.
625
C ---------------------------------------------
627
IF (NEWQ .EQ. NQ) GO TO 650
631
C ...............EXIT
635
R = 1.0D0/TESCO(2,NQU)
639
C ..................EXIT
642
C ---------------------------------------------------
643
C CONTROL REACHES THIS SECTION IF 3 OR MORE FAILURES
644
C HAVE OCCURRED. IF 10 FAILURES HAVE OCCURRED, EXIT
645
C WITH KFLAG = -1. IT IS ASSUMED THAT THE
646
C DERIVATIVES THAT HAVE ACCUMULATED IN THE YH ARRAY
647
C HAVE ERRORS OF THE WRONG ORDER. HENCE THE FIRST
648
C DERIVATIVE IS RECOMPUTED, AND THE ORDER IS SET TO
649
C 1. THEN H IS REDUCED BY A FACTOR OF 10, AND THE
650
C STEP IS RETRIED, UNTIL IT SUCCEEDS OR H REACHES
652
C ---------------------------------------------------
653
IF (KFLAG .NE. -10) GO TO 620
654
C ------------------------------------------------
655
C ALL RETURNS ARE MADE THROUGH THIS SECTION. H
656
C IS SAVED IN HOLD TO ALLOW THE CALLER TO CHANGE
657
C H ON THE NEXT STEP.
658
C ------------------------------------------------
660
C ..................EXIT
664
RH = MAX(HMIN/ABS(H),RH)
669
CALL DF(TN,Y,SAVF,RPAR,IPAR)
677
IF (NQ .NE. 1) GO TO 670
681
RH = MAX(RH,HMIN/ABS(H))
693
C ----------------------- END OF SUBROUTINE DSTOD
694
C -----------------------
added
removed
deps/openlibm/slatec/dstod.f
