2
SUBROUTINE DSTOKA (NEQ, Y, YH, NYH, YH1, EWT, SAVF, SAVX, ACOR,
3
1 WM, IWM, F, JAC, PSOL)
6
DOUBLE PRECISION Y, YH, YH1, EWT, SAVF, SAVX, ACOR, WM
7
DIMENSION NEQ(*), Y(*), YH(NYH,*), YH1(*), EWT(*), SAVF(*),
8
1 SAVX(*), ACOR(*), WM(*), IWM(*)
9
INTEGER IOWND, IALTH, IPUP, LMAX, MEO, NQNYH, NSLP,
10
1 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L,
11
2 LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER,
12
3 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
13
INTEGER NEWT, NSFI, NSLJ, NJEV
14
INTEGER JPRE, JACFLG, LOCWP, LOCIWP, LSAVX, KMP, MAXL, MNEWT,
15
1 NNI, NLI, NPS, NCFN, NCFL
16
DOUBLE PRECISION CONIT, CRATE, EL, ELCO, HOLD, RMAX, TESCO,
17
2 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
18
DOUBLE PRECISION STIFR
19
DOUBLE PRECISION DELT, EPCON, SQRTN, RSQRTN
20
COMMON /DLS001/ CONIT, CRATE, EL(13), ELCO(13,12),
21
1 HOLD, RMAX, TESCO(3,12),
22
2 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND,
23
3 IOWND(6), IALTH, IPUP, LMAX, MEO, NQNYH, NSLP,
24
3 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L,
25
4 LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER,
26
5 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
27
COMMON /DLS002/ STIFR, NEWT, NSFI, NSLJ, NJEV
28
COMMON /DLPK01/ DELT, EPCON, SQRTN, RSQRTN,
29
1 JPRE, JACFLG, LOCWP, LOCIWP, LSAVX, KMP, MAXL, MNEWT,
30
2 NNI, NLI, NPS, NCFN, NCFL
31
C-----------------------------------------------------------------------
32
C DSTOKA performs one step of the integration of an initial value
33
C problem for a system of Ordinary Differential Equations.
35
C This routine was derived from Subroutine DSTODPK in the DLSODPK
36
C package by the addition of automatic functional/Newton iteration
37
C switching and logic for re-use of Jacobian data.
38
C-----------------------------------------------------------------------
39
C Note: DSTOKA is independent of the value of the iteration method
40
C indicator MITER, when this is .ne. 0, and hence is independent
41
C of the type of chord method used, or the Jacobian structure.
42
C Communication with DSTOKA is done with the following variables:
44
C NEQ = integer array containing problem size in NEQ(1), and
45
C passed as the NEQ argument in all calls to F and JAC.
46
C Y = an array of length .ge. N used as the Y argument in
47
C all calls to F and JAC.
48
C YH = an NYH by LMAX array containing the dependent variables
49
C and their approximate scaled derivatives, where
50
C LMAX = MAXORD + 1. YH(i,j+1) contains the approximate
51
C j-th derivative of y(i), scaled by H**j/factorial(j)
52
C (j = 0,1,...,NQ). On entry for the first step, the first
53
C two columns of YH must be set from the initial values.
54
C NYH = a constant integer .ge. N, the first dimension of YH.
55
C YH1 = a one-dimensional array occupying the same space as YH.
56
C EWT = an array of length N containing multiplicative weights
57
C for local error measurements. Local errors in y(i) are
58
C compared to 1.0/EWT(i) in various error tests.
59
C SAVF = an array of working storage, of length N.
60
C Also used for input of YH(*,MAXORD+2) when JSTART = -1
61
C and MAXORD .lt. the current order NQ.
62
C SAVX = an array of working storage, of length N.
63
C ACOR = a work array of length N, used for the accumulated
64
C corrections. On a successful return, ACOR(i) contains
65
C the estimated one-step local error in y(i).
66
C WM,IWM = real and integer work arrays associated with matrix
67
C operations in chord iteration (MITER .ne. 0).
68
C CCMAX = maximum relative change in H*EL0 before DSETPK is called.
69
C H = the step size to be attempted on the next step.
70
C H is altered by the error control algorithm during the
71
C problem. H can be either positive or negative, but its
72
C sign must remain constant throughout the problem.
73
C HMIN = the minimum absolute value of the step size H to be used.
74
C HMXI = inverse of the maximum absolute value of H to be used.
75
C HMXI = 0.0 is allowed and corresponds to an infinite HMAX.
76
C HMIN and HMXI may be changed at any time, but will not
77
C take effect until the next change of H is considered.
78
C TN = the independent variable. TN is updated on each step taken.
79
C JSTART = an integer used for input only, with the following
80
C values and meanings:
81
C 0 perform the first step.
82
C .gt.0 take a new step continuing from the last.
83
C -1 take the next step with a new value of H, MAXORD,
84
C N, METH, MITER, and/or matrix parameters.
85
C -2 take the next step with a new value of H,
86
C but with other inputs unchanged.
87
C On return, JSTART is set to 1 to facilitate continuation.
88
C KFLAG = a completion code with the following meanings:
89
C 0 the step was succesful.
90
C -1 the requested error could not be achieved.
91
C -2 corrector convergence could not be achieved.
92
C -3 fatal error in DSETPK or DSOLPK.
93
C A return with KFLAG = -1 or -2 means either
94
C ABS(H) = HMIN or 10 consecutive failures occurred.
95
C On a return with KFLAG negative, the values of TN and
96
C the YH array are as of the beginning of the last
97
C step, and H is the last step size attempted.
98
C MAXORD = the maximum order of integration method to be allowed.
99
C MAXCOR = the maximum number of corrector iterations allowed.
100
C MSBP = maximum number of steps between DSETPK calls (MITER .gt. 0).
101
C MXNCF = maximum number of convergence failures allowed.
102
C METH/MITER = the method flags. See description in driver.
103
C N = the number of first-order differential equations.
104
C-----------------------------------------------------------------------
105
INTEGER I, I1, IREDO, IRET, J, JB, JOK, M, NCF, NEWQ, NSLOW
106
DOUBLE PRECISION DCON, DDN, DEL, DELP, DRC, DSM, DUP, EXDN, EXSM,
107
1 EXUP, DFNORM, R, RH, RHDN, RHSM, RHUP, ROC, STIFF, TOLD, DVNORM
117
IF (JSTART .GT. 0) GO TO 200
118
IF (JSTART .EQ. -1) GO TO 100
119
IF (JSTART .EQ. -2) GO TO 160
120
C-----------------------------------------------------------------------
121
C On the first call, the order is set to 1, and other variables are
122
C initialized. RMAX is the maximum ratio by which H can be increased
123
C in a single step. It is initially 1.E4 to compensate for the small
124
C initial H, but then is normally equal to 10. If a failure
125
C occurs (in corrector convergence or error test), RMAX is set at 2
126
C for the next increase.
127
C-----------------------------------------------------------------------
145
C-----------------------------------------------------------------------
146
C The following block handles preliminaries needed when JSTART = -1.
147
C IPUP is set to MITER to force a matrix update.
148
C If an order increase is about to be considered (IALTH = 1),
149
C IALTH is reset to 2 to postpone consideration one more step.
150
C If the caller has changed METH, DCFODE is called to reset
151
C the coefficients of the method.
152
C If the caller has changed MAXORD to a value less than the current
153
C order NQ, NQ is reduced to MAXORD, and a new H chosen accordingly.
154
C If H is to be changed, YH must be rescaled.
155
C If H or METH is being changed, IALTH is reset to L = NQ + 1
156
C to prevent further changes in H for that many steps.
157
C-----------------------------------------------------------------------
160
IF (IALTH .EQ. 1) IALTH = 2
161
IF (METH .EQ. MEO) GO TO 110
162
CALL DCFODE (METH, ELCO, TESCO)
164
IF (NQ .GT. MAXORD) GO TO 120
168
110 IF (NQ .LE. MAXORD) GO TO 160
172
125 EL(I) = ELCO(I,NQ)
177
EPCON = CONIT*TESCO(2,NQ)
178
DDN = DVNORM (N, SAVF, EWT)/TESCO(1,L)
180
RHDN = 1.0D0/(1.3D0*DDN**EXDN + 0.0000013D0)
183
IF (H .EQ. HOLD) GO TO 170
184
RH = MIN(RH,ABS(H/HOLD))
187
C-----------------------------------------------------------------------
188
C DCFODE is called to get all the integration coefficients for the
189
C current METH. Then the EL vector and related constants are reset
190
C whenever the order NQ is changed, or at the start of the problem.
191
C-----------------------------------------------------------------------
192
140 CALL DCFODE (METH, ELCO, TESCO)
194
155 EL(I) = ELCO(I,NQ)
199
EPCON = CONIT*TESCO(2,NQ)
200
GO TO (160, 170, 200), IRET
201
C-----------------------------------------------------------------------
202
C If H is being changed, the H ratio RH is checked against
203
C RMAX, HMIN, and HMXI, and the YH array rescaled. IALTH is set to
204
C L = NQ + 1 to prevent a change of H for that many steps, unless
205
C forced by a convergence or error test failure.
206
C-----------------------------------------------------------------------
207
160 IF (H .EQ. HOLD) GO TO 200
212
170 RH = MAX(RH,HMIN/ABS(H))
213
175 RH = MIN(RH,RMAX)
214
RH = RH/MAX(1.0D0,ABS(H)*HMXI*RH)
219
180 YH(I,J) = YH(I,J)*R
223
IF (IREDO .EQ. 0) GO TO 690
224
C-----------------------------------------------------------------------
225
C This section computes the predicted values by effectively
226
C multiplying the YH array by the Pascal triangle matrix.
227
C The flag IPUP is set according to whether matrix data is involved
228
C (NEWT .gt. 0 .and. JACFLG .ne. 0) or not, to trigger a call to DSETPK.
229
C IPUP is set to MITER when RC differs from 1 by more than CCMAX,
230
C and at least every MSBP steps, when JACFLG = 1.
231
C RC is the ratio of new to old values of the coefficient H*EL(1).
232
C-----------------------------------------------------------------------
233
200 IF (NEWT .EQ. 0 .OR. JACFLG .EQ. 0) THEN
238
DRC = ABS(RC - 1.0D0)
239
IF (DRC .GT. CCMAX) IPUP = MITER
240
IF (NST .GE. NSLP+MSBP) IPUP = MITER
248
210 YH1(I) = YH1(I) + YH1(I+NYH)
250
C-----------------------------------------------------------------------
251
C Up to MAXCOR corrector iterations are taken. A convergence test is
252
C made on the RMS-norm of each correction, weighted by the error
253
C weight vector EWT. The sum of the corrections is accumulated in the
254
C vector ACOR(i). The YH array is not altered in the corrector loop.
255
C Within the corrector loop, an estimated rate of convergence (ROC)
256
C and a stiffness ratio estimate (STIFF) are kept. Corresponding
257
C global estimates are kept as CRATE and stifr.
258
C-----------------------------------------------------------------------
266
CALL F (NEQ, TN, Y, SAVF)
268
IF (NEWT .EQ. 0 .OR. IPUP .LE. 0) GO TO 250
269
C-----------------------------------------------------------------------
270
C If indicated, DSETPK is called to update any matrix data needed,
271
C before starting the corrector iteration.
272
C JOK is set to indicate if the matrix data need not be recomputed.
273
C IPUP is set to 0 as an indicator that the matrix data is up to date.
274
C-----------------------------------------------------------------------
276
IF (NST .EQ. 0 .OR. NST .GT. NSLJ+50) JOK = -1
277
IF (ICF .EQ. 1 .AND. DRC .LT. 0.2D0) JOK = -1
278
IF (ICF .EQ. 2) JOK = -1
279
IF (JOK .EQ. -1) THEN
283
CALL DSETPK (NEQ, Y, YH1, EWT, ACOR, SAVF, JOK, WM, IWM, F, JAC)
289
IF (IERPJ .NE. 0) GO TO 430
292
270 IF (NEWT .NE. 0) GO TO 350
293
C-----------------------------------------------------------------------
294
C In the case of functional iteration, update Y directly from
295
C the result of the last function evaluation, and STIFF is set to 1.0.
296
C-----------------------------------------------------------------------
298
SAVF(I) = H*SAVF(I) - YH(I,2)
299
290 Y(I) = SAVF(I) - ACOR(I)
300
DEL = DVNORM (N, Y, EWT)
302
Y(I) = YH(I,1) + EL(1)*SAVF(I)
303
300 ACOR(I) = SAVF(I)
306
C-----------------------------------------------------------------------
307
C In the case of the chord method, compute the corrector error,
308
C and solve the linear system with that as right-hand side and
309
C P as coefficient matrix. STIFF is set to the ratio of the norms
310
C of the residual and the correction vector.
311
C-----------------------------------------------------------------------
313
360 SAVX(I) = H*SAVF(I) - (YH(I,2) + ACOR(I))
314
DFNORM = DVNORM (N, SAVX, EWT)
315
CALL DSOLPK (NEQ, Y, SAVF, SAVX, EWT, WM, IWM, F, PSOL)
316
IF (IERSL .LT. 0) GO TO 430
317
IF (IERSL .GT. 0) GO TO 410
318
DEL = DVNORM (N, SAVX, EWT)
319
IF (DEL .GT. 1.0D-8) STIFF = MAX(STIFF, DFNORM/DEL)
321
ACOR(I) = ACOR(I) + SAVX(I)
322
380 Y(I) = YH(I,1) + EL(1)*ACOR(I)
323
C-----------------------------------------------------------------------
324
C Test for convergence. If M .gt. 0, an estimate of the convergence
325
C rate constant is made for the iteration switch, and is also used
326
C in the convergence test. If the iteration seems to be diverging or
327
C converging at a slow rate (.gt. 0.8 more than once), it is stopped.
328
C-----------------------------------------------------------------------
329
400 IF (M .NE. 0) THEN
330
ROC = MAX(0.05D0, DEL/DELP)
331
CRATE = MAX(0.2D0*CRATE,ROC)
333
DCON = DEL*MIN(1.0D0,1.5D0*CRATE)/EPCON
334
IF (DCON .LE. 1.0D0) GO TO 450
336
IF (M .EQ. MAXCOR) GO TO 410
337
IF (M .GE. 2 .AND. DEL .GT. 2.0D0*DELP) GO TO 410
338
IF (ROC .GT. 10.0D0) GO TO 410
339
IF (ROC .GT. 0.8D0) NSLOW = NSLOW + 1
340
IF (NSLOW .GE. 2) GO TO 410
343
CALL F (NEQ, TN, Y, SAVF)
346
C-----------------------------------------------------------------------
347
C The corrector iteration failed to converge.
348
C If functional iteration is being done (NEWT = 0) and MITER .gt. 0
349
C (and this is not the first step), then switch to Newton
350
C (NEWT = MITER), and retry the step. (Setting STIFR = 1023 insures
351
C that a switch back will not occur for 10 step attempts.)
352
C If Newton iteration is being done, but using a preconditioner that
353
C is out of date (JACFLG .ne. 0 .and. JCUR = 0), then signal for a
354
C re-evalutation of the preconditioner, and retry the step.
355
C In all other cases, the YH array is retracted to its values
356
C before prediction, and H is reduced, if possible. If H cannot be
357
C reduced or MXNCF failures have occurred, exit with KFLAG = -2.
358
C-----------------------------------------------------------------------
360
IF (NEWT .EQ. 0) THEN
361
IF (NST .EQ. 0) GO TO 430
362
IF (MITER .EQ. 0) GO TO 430
368
IF (JCUR.EQ.1 .OR. JACFLG.EQ.0) GO TO 430
381
440 YH1(I) = YH1(I) - YH1(I+NYH)
383
IF (IERPJ .LT. 0 .OR. IERSL .LT. 0) GO TO 680
384
IF (ABS(H) .LE. HMIN*1.00001D0) GO TO 670
385
IF (NCF .EQ. MXNCF) GO TO 670
390
C-----------------------------------------------------------------------
391
C The corrector has converged. JCUR is set to 0 to signal that the
392
C preconditioner involved may need updating later.
393
C The stiffness ratio STIFR is updated using the latest STIFF value.
394
C The local error test is made and control passes to statement 500
396
C-----------------------------------------------------------------------
398
IF (NEWT .GT. 0) STIFR = 0.5D0*(STIFR + STIFF)
399
IF (M .EQ. 0) DSM = DEL/TESCO(2,NQ)
400
IF (M .GT. 0) DSM = DVNORM (N, ACOR, EWT)/TESCO(2,NQ)
401
IF (DSM .GT. 1.0D0) GO TO 500
402
C-----------------------------------------------------------------------
403
C After a successful step, update the YH array.
404
C If Newton iteration is being done and STIFR is less than 1.5,
405
C then switch to functional iteration.
406
C Consider changing H if IALTH = 1. Otherwise decrease IALTH by 1.
407
C If IALTH is then 1 and NQ .lt. MAXORD, then ACOR is saved for
408
C use in a possible order increase on the next step.
409
C If a change in H is considered, an increase or decrease in order
410
C by one is considered also. A change in H is made only if it is by a
411
C factor of at least 1.1. If not, IALTH is set to 3 to prevent
412
C testing for that many steps.
413
C-----------------------------------------------------------------------
417
IF (NEWT .EQ. 0) NSFI = NSFI + 1
418
IF (NEWT .GT. 0 .AND. STIFR .LT. 1.5D0) NEWT = 0
423
470 YH(I,J) = YH(I,J) + EL(J)*ACOR(I)
425
IF (IALTH .EQ. 0) GO TO 520
426
IF (IALTH .GT. 1) GO TO 700
427
IF (L .EQ. LMAX) GO TO 700
429
490 YH(I,LMAX) = ACOR(I)
431
C-----------------------------------------------------------------------
432
C The error test failed. KFLAG keeps track of multiple failures.
433
C Restore TN and the YH array to their previous values, and prepare
434
C to try the step again. Compute the optimum step size for this or
435
C one lower order. After 2 or more failures, H is forced to decrease
436
C by a factor of 0.2 or less.
437
C-----------------------------------------------------------------------
438
500 KFLAG = KFLAG - 1
445
510 YH1(I) = YH1(I) - YH1(I+NYH)
448
IF (ABS(H) .LE. HMIN*1.00001D0) GO TO 660
449
IF (KFLAG .LE. -3) GO TO 640
453
C-----------------------------------------------------------------------
454
C Regardless of the success or failure of the step, factors
455
C RHDN, RHSM, and RHUP are computed, by which H could be multiplied
456
C at order NQ - 1, order NQ, or order NQ + 1, respectively.
457
C in the case of failure, RHUP = 0.0 to avoid an order increase.
458
C the largest of these is determined and the new order chosen
459
C accordingly. If the order is to be increased, we compute one
460
C additional scaled derivative.
461
C-----------------------------------------------------------------------
463
IF (L .EQ. LMAX) GO TO 540
465
530 SAVF(I) = ACOR(I) - YH(I,LMAX)
466
DUP = DVNORM (N, SAVF, EWT)/TESCO(3,NQ)
468
RHUP = 1.0D0/(1.4D0*DUP**EXUP + 0.0000014D0)
470
RHSM = 1.0D0/(1.2D0*DSM**EXSM + 0.0000012D0)
472
IF (NQ .EQ. 1) GO TO 560
473
DDN = DVNORM (N, YH(1,L), EWT)/TESCO(1,NQ)
475
RHDN = 1.0D0/(1.3D0*DDN**EXDN + 0.0000013D0)
476
560 IF (RHSM .GE. RHUP) GO TO 570
477
IF (RHUP .GT. RHDN) GO TO 590
479
570 IF (RHSM .LT. RHDN) GO TO 580
485
IF (KFLAG .LT. 0 .AND. RH .GT. 1.0D0) RH = 1.0D0
489
IF (RH .LT. 1.1D0) GO TO 610
492
600 YH(I,NEWQ+1) = ACOR(I)*R
496
620 IF ((KFLAG .EQ. 0) .AND. (RH .LT. 1.1D0)) GO TO 610
497
IF (KFLAG .LE. -2) RH = MIN(RH,0.2D0)
498
C-----------------------------------------------------------------------
499
C If there is a change of order, reset NQ, L, and the coefficients.
500
C In any case H is reset according to RH and the YH array is rescaled.
501
C Then exit from 690 if the step was OK, or redo the step otherwise.
502
C-----------------------------------------------------------------------
503
IF (NEWQ .EQ. NQ) GO TO 170
508
C-----------------------------------------------------------------------
509
C Control reaches this section if 3 or more failures have occured.
510
C If 10 failures have occurred, exit with KFLAG = -1.
511
C It is assumed that the derivatives that have accumulated in the
512
C YH array have errors of the wrong order. Hence the first
513
C derivative is recomputed, and the order is set to 1. Then
514
C H is reduced by a factor of 10, and the step is retried,
515
C until it succeeds or H reaches HMIN.
516
C-----------------------------------------------------------------------
517
640 IF (KFLAG .EQ. -10) GO TO 660
519
RH = MAX(HMIN/ABS(H),RH)
523
CALL F (NEQ, TN, Y, SAVF)
526
650 YH(I,2) = H*SAVF(I)
529
IF (NQ .EQ. 1) GO TO 200
534
C-----------------------------------------------------------------------
535
C All returns are made through this section. H is saved in HOLD
536
C to allow the caller to change H on the next step.
537
C-----------------------------------------------------------------------
545
700 R = 1.0D0/TESCO(2,NQU)
547
710 ACOR(I) = ACOR(I)*R
551
C----------------------- End of Subroutine DSTOKA ----------------------
added
removed
packages/odepack/dstoka.f
