*DECK DDSTP
SUBROUTINE DDSTP (EPS, F, FA, HMAX, IMPL, IERROR, JACOBN, MATDIM,
8 MAXORD, MINT, MITER, ML, MU, N, NDE, YWT, UROUND, USERS, AVGH,
8 AVGORD, H, HUSED, JTASK, MNTOLD, MTROLD, NFE, NJE, NQUSED,
8 NSTEP, T, Y, YH, A, CONVRG, DFDY, EL, FAC, HOLD, IPVT, JSTATE,
8 JSTEPL, NQ, NWAIT, RC, RMAX, SAVE1, SAVE2, TQ, TREND, ISWFLG,
8 MTRSV, MXRDSV)
C***BEGIN PROLOGUE DDSTP
C***SUBSIDIARY
C***PURPOSE DDSTP performs one step of the integration of an initial
C value problem for a system of ordinary differential
C equations.
C***LIBRARY SLATEC (SDRIVE)
C***TYPE DOUBLE PRECISION (SDSTP-S, DDSTP-D, CDSTP-C)
C***AUTHOR Kahaner, D. K., (NIST)
C National Institute of Standards and Technology
C Gaithersburg, MD 20899
C Sutherland, C. D., (LANL)
C Mail Stop D466
C Los Alamos National Laboratory
C Los Alamos, NM 87545
C***DESCRIPTION
C
C Communication with DDSTP is done with the following variables:
C
C YH An N by MAXORD+1 array containing the dependent variables
C and their scaled derivatives. MAXORD, the maximum order
C used, is currently 12 for the Adams methods and 5 for the
C Gear methods. YH(I,J+1) contains the J-th derivative of
C Y(I), scaled by H**J/factorial(J). Only Y(I),
C 1 .LE. I .LE. N, need be set by the calling program on
C the first entry. The YH array should not be altered by
C the calling program. When referencing YH as a
C 2-dimensional array, use a column length of N, as this is
C the value used in DDSTP.
C DFDY A block of locations used for partial derivatives if MITER
C is not 0. If MITER is 1 or 2 its length must be at least
C N*N. If MITER is 4 or 5 its length must be at least
C (2*ML+MU+1)*N.
C YWT An array of N locations used in convergence and error tests
C SAVE1
C SAVE2 Arrays of length N used for temporary storage.
C IPVT An integer array of length N used by the linear system
C solvers for the storage of row interchange information.
C A A block of locations used to store the matrix A, when using
C the implicit method. If IMPL is 1, A is a MATDIM by N
C array. If MITER is 1 or 2 MATDIM is N, and if MITER is 4
C or 5 MATDIM is 2*ML+MU+1. If IMPL is 2 its length is N.
C If IMPL is 3, A is a MATDIM by NDE array.
C JTASK An integer used on input.
C It has the following values and meanings:
C .EQ. 0 Perform the first step. This value enables
C the subroutine to initialize itself.
C .GT. 0 Take a new step continuing from the last.
C Assumes the last step was successful and
C user has not changed any parameters.
C .LT. 0 Take a new step with a new value of H and/or
C MINT and/or MITER.
C JSTATE A completion code with the following meanings:
C 1 The step was successful.
C 2 A solution could not be obtained with H .NE. 0.
C 3 A solution was not obtained in MXTRY attempts.
C 4 For IMPL .NE. 0, the matrix A is singular.
C On a return with JSTATE .GT. 1, the values of T and
C the YH array are as of the beginning of the last
C step, and H is the last step size attempted.
C
C***ROUTINES CALLED DDCOR, DDCST, DDNTL, DDPSC, DDPST, DDSCL, DNRM2
C***REVISION HISTORY (YYMMDD)
C 790601 DATE WRITTEN
C 900329 Initial submission to SLATEC.
C***END PROLOGUE DDSTP
EXTERNAL F, JACOBN, FA, USERS
INTEGER I, IERROR, IMPL, IPVT(*), ISWFLG, ITER, J, JSTATE, JSTEPL,
8 JTASK, MATDIM, MAXORD, MINT, MITER, ML, MNTOLD, MTROLD,
8 MTRSV, MU, MXFAIL, MXITER, MXRDSV, MXTRY, N, NDE, NDJSTP,
8 NFAIL, NFE, NJE, NQ, NQUSED, NSTEP, NSV, NTRY, NWAIT
DOUBLE PRECISION A(MATDIM,*), AVGH, AVGORD, BIAS1, BIAS2, BIAS3,
8 BND, CTEST, D, DENOM, DFDY(MATDIM,*), D1, EL(13,12), EPS,
8 ERDN, ERUP, ETEST, FAC(*), H, HMAX, HN, HOLD, HS, HUSED,
8 NUMER, RC, RCTEST, RH, RH1, RH2, RH3, RMAX, RMFAIL, RMNORM,
8 SAVE1(*), SAVE2(*), DNRM2, T, TOLD, TQ(3,12), TREND, TRSHLD,
8 UROUND, Y(*), YH(N,*), YWT(*), Y0NRM
LOGICAL CONVRG, EVALFA, EVALJC, IER, SWITCH
PARAMETER(BIAS1 = 1.3D0, BIAS2 = 1.2D0, BIAS3 = 1.4D0, MXFAIL = 3,
8 MXITER = 3, MXTRY = 50, RCTEST = .3D0, RMFAIL = 2.D0,
8 RMNORM = 10.D0, TRSHLD = 1.D0)
PARAMETER (NDJSTP = 10)
DATA IER /.FALSE./
C***FIRST EXECUTABLE STATEMENT DDSTP
NSV = N
BND = 0.D0
SWITCH = .FALSE.
NTRY = 0
TOLD = T
NFAIL = 0
IF (JTASK .LE. 0) THEN
CALL DDNTL (EPS, F, FA, HMAX, HOLD, IMPL, JTASK, MATDIM,
8 MAXORD, MINT, MITER, ML, MU, N, NDE, SAVE1, T,
8 UROUND, USERS, Y, YWT, H, MNTOLD, MTROLD, NFE, RC,
8 YH, A, CONVRG, EL, FAC, IER, IPVT, NQ, NWAIT, RH,
8 RMAX, SAVE2, TQ, TREND, ISWFLG, JSTATE)
IF (N .EQ. 0) GO TO 440
IF (H .EQ. 0.D0) GO TO 400
IF (IER) GO TO 420
END IF
100 NTRY = NTRY + 1
IF (NTRY .GT. MXTRY) GO TO 410
T = T + H
CALL DDPSC (1, N, NQ, YH)
EVALJC = (((ABS(RC - 1.D0) .GT. RCTEST) .OR.
8 (NSTEP .GE. JSTEPL + NDJSTP)) .AND. (MITER .NE. 0))
EVALFA = .NOT. EVALJC
C
110 ITER = 0
DO 115 I = 1,N
115 Y(I) = YH(I,1)
CALL F (N, T, Y, SAVE2)
IF (N .EQ. 0) THEN
JSTATE = 6
GO TO 430
END IF
NFE = NFE + 1
IF (EVALJC .OR. IER) THEN
CALL DDPST (EL, F, FA, H, IMPL, JACOBN, MATDIM, MITER, ML,
8 MU, N, NDE, NQ, SAVE2, T, USERS, Y, YH, YWT, UROUND,
8 NFE, NJE, A, DFDY, FAC, IER, IPVT, SAVE1, ISWFLG,
8 BND, JSTATE)
IF (N .EQ. 0) GO TO 430
IF (IER) GO TO 160
CONVRG = .FALSE.
RC = 1.D0
JSTEPL = NSTEP
END IF
DO 125 I = 1,N
125 SAVE1(I) = 0.D0
C Up to MXITER corrector iterations are taken.
C Convergence is tested by requiring the r.m.s.
C norm of changes to be less than EPS. The sum of
C the corrections is accumulated in the vector
C SAVE1(I). It is approximately equal to the L-th
C derivative of Y multiplied by
C H**L/(factorial(L-1)*EL(L,NQ)), and is thus
C proportional to the actual errors to the lowest
C power of H present (H**L). The YH array is not
C altered in the correction loop. The norm of the
C iterate difference is stored in D. If
C ITER .GT. 0, an estimate of the convergence rate
C constant is stored in TREND, and this is used in
C the convergence test.
C
130 CALL DDCOR (DFDY, EL, FA, H, IERROR, IMPL, IPVT, MATDIM, MITER,
8 ML, MU, N, NDE, NQ, T, USERS, Y, YH, YWT, EVALFA,
8 SAVE1, SAVE2, A, D, JSTATE)
IF (N .EQ. 0) GO TO 430
IF (ISWFLG .EQ. 3 .AND. MINT .EQ. 1) THEN
IF (ITER .EQ. 0) THEN
NUMER = DNRM2(N, SAVE1, 1)
DO 132 I = 1,N
132 DFDY(1,I) = SAVE1(I)
Y0NRM = DNRM2(N, YH, 1)
ELSE
DENOM = NUMER
DO 134 I = 1,N
134 DFDY(1,I) = SAVE1(I) - DFDY(1,I)
NUMER = DNRM2(N, DFDY, MATDIM)
IF (EL(1,NQ)*NUMER .LE. 100.D0*UROUND*Y0NRM) THEN
IF (RMAX .EQ. RMFAIL) THEN
SWITCH = .TRUE.
GO TO 170
END IF
END IF
DO 136 I = 1,N
136 DFDY(1,I) = SAVE1(I)
IF (DENOM .NE. 0.D0)
8 BND = MAX(BND, NUMER/(DENOM*ABS(H)*EL(1,NQ)))
END IF
END IF
IF (ITER .GT. 0) TREND = MAX(.9D0*TREND, D/D1)
D1 = D
CTEST = MIN(2.D0*TREND, 1.D0)*D
IF (CTEST .LE. EPS) GO TO 170
ITER = ITER + 1
IF (ITER .LT. MXITER) THEN
DO 140 I = 1,N
140 Y(I) = YH(I,1) + EL(1,NQ)*SAVE1(I)
CALL F (N, T, Y, SAVE2)
IF (N .EQ. 0) THEN
JSTATE = 6
GO TO 430
END IF
NFE = NFE + 1
GO TO 130
END IF
C The corrector iteration failed to converge in
C MXITER tries. If partials are involved but are
C not up to date, they are reevaluated for the next
C try. Otherwise the YH array is retracted to its
C values before prediction, and H is reduced, if
C possible. If not, a no-convergence exit is taken.
IF (CONVRG) THEN
EVALJC = .TRUE.
EVALFA = .FALSE.
GO TO 110
END IF
160 T = TOLD
CALL DDPSC (-1, N, NQ, YH)
NWAIT = NQ + 2
IF (JTASK .NE. 0 .AND. JTASK .NE. 2) RMAX = RMFAIL
IF (ITER .EQ. 0) THEN
RH = .3D0
ELSE
RH = .9D0*(EPS/CTEST)**(.2D0)
END IF
IF (RH*H .EQ. 0.D0) GO TO 400
CALL DDSCL (HMAX, N, NQ, RMAX, H, RC, RH, YH)
GO TO 100
C The corrector has converged. CONVRG is set
C to .TRUE. if partial derivatives were used,
C to indicate that they may need updating on
C subsequent steps. The error test is made.
170 CONVRG = (MITER .NE. 0)
IF (IERROR .EQ. 1 .OR. IERROR .EQ. 5) THEN
DO 180 I = 1,NDE
180 SAVE2(I) = SAVE1(I)/YWT(I)
ELSE
DO 185 I = 1,NDE
185 SAVE2(I) = SAVE1(I)/MAX(ABS(Y(I)), YWT(I))
END IF
ETEST = DNRM2(NDE, SAVE2, 1)/(TQ(2,NQ)*SQRT(DBLE(NDE)))
C
C The error test failed. NFAIL keeps track of
C multiple failures. Restore T and the YH
C array to their previous values, and prepare
C to try the step again. Compute the optimum
C step size for this or one lower order.
IF (ETEST .GT. EPS) THEN
T = TOLD
CALL DDPSC (-1, N, NQ, YH)
NFAIL = NFAIL + 1
IF (NFAIL .LT. MXFAIL .OR. NQ .EQ. 1) THEN
IF (JTASK .NE. 0 .AND. JTASK .NE. 2) RMAX = RMFAIL
RH2 = 1.D0/(BIAS2*(ETEST/EPS)**(1.D0/(NQ+1)))
IF (NQ .GT. 1) THEN
IF (IERROR .EQ. 1 .OR. IERROR .EQ. 5) THEN
DO 190 I = 1,NDE
190 SAVE2(I) = YH(I,NQ+1)/YWT(I)
ELSE
DO 195 I = 1,NDE
195 SAVE2(I) = YH(I,NQ+1)/MAX(ABS(Y(I)), YWT(I))
END IF
ERDN = DNRM2(NDE, SAVE2, 1)/(TQ(1,NQ)*SQRT(DBLE(NDE)))
RH1 = 1.D0/MAX(1.D0, BIAS1*(ERDN/EPS)**(1.D0/NQ))
IF (RH2 .LT. RH1) THEN
NQ = NQ - 1
RC = RC*EL(1,NQ)/EL(1,NQ+1)
RH = RH1
ELSE
RH = RH2
END IF
ELSE
RH = RH2
END IF
NWAIT = NQ + 2
IF (RH*H .EQ. 0.D0) GO TO 400
CALL DDSCL (HMAX, N, NQ, RMAX, H, RC, RH, YH)
GO TO 100
END IF
C Control reaches this section if the error test has
C failed MXFAIL or more times. It is assumed that the
C derivatives that have accumulated in the YH array have
C errors of the wrong order. Hence the first derivative
C is recomputed, the order is set to 1, and the step is
C retried.
NFAIL = 0
JTASK = 2
DO 215 I = 1,N
215 Y(I) = YH(I,1)
CALL DDNTL (EPS, F, FA, HMAX, HOLD, IMPL, JTASK, MATDIM,
8 MAXORD, MINT, MITER, ML, MU, N, NDE, SAVE1, T,
8 UROUND, USERS, Y, YWT, H, MNTOLD, MTROLD, NFE, RC,
8 YH, A, CONVRG, EL, FAC, IER, IPVT, NQ, NWAIT, RH,
8 RMAX, SAVE2, TQ, TREND, ISWFLG, JSTATE)
RMAX = RMNORM
IF (N .EQ. 0) GO TO 440
IF (H .EQ. 0.D0) GO TO 400
IF (IER) GO TO 420
GO TO 100
END IF
C After a successful step, update the YH array.
NSTEP = NSTEP + 1
HUSED = H
NQUSED = NQ
AVGH = ((NSTEP-1)*AVGH + H)/NSTEP
AVGORD = ((NSTEP-1)*AVGORD + NQ)/NSTEP
DO 230 J = 1,NQ+1
DO 230 I = 1,N
230 YH(I,J) = YH(I,J) + EL(J,NQ)*SAVE1(I)
DO 235 I = 1,N
235 Y(I) = YH(I,1)
C If ISWFLG is 3, consider
C changing integration methods.
IF (ISWFLG .EQ. 3) THEN
IF (BND .NE. 0.D0) THEN
IF (MINT .EQ. 1 .AND. NQ .LE. 5) THEN
HN = ABS(H)/MAX(UROUND, (ETEST/EPS)**(1.D0/(NQ+1)))
HN = MIN(HN, 1.D0/(2.D0*EL(1,NQ)*BND))
HS = ABS(H)/MAX(UROUND,
8 (ETEST/(EPS*EL(NQ+1,1)))**(1.D0/(NQ+1)))
IF (HS .GT. 1.2D0*HN) THEN
MINT = 2
MNTOLD = MINT
MITER = MTRSV
MTROLD = MITER
MAXORD = MIN(MXRDSV, 5)
RC = 0.D0
RMAX = RMNORM
TREND = 1.D0
CALL DDCST (MAXORD, MINT, ISWFLG, EL, TQ)
NWAIT = NQ + 2
END IF
ELSE IF (MINT .EQ. 2) THEN
HS = ABS(H)/MAX(UROUND, (ETEST/EPS)**(1.D0/(NQ+1)))
HN = ABS(H)/MAX(UROUND,
8 (ETEST*EL(NQ+1,1)/EPS)**(1.D0/(NQ+1)))
HN = MIN(HN, 1.D0/(2.D0*EL(1,NQ)*BND))
IF (HN .GE. HS) THEN
MINT = 1
MNTOLD = MINT
MITER = 0
MTROLD = MITER
MAXORD = MIN(MXRDSV, 12)
RMAX = RMNORM
TREND = 1.D0
CONVRG = .FALSE.
CALL DDCST (MAXORD, MINT, ISWFLG, EL, TQ)
NWAIT = NQ + 2
END IF
END IF
END IF
END IF
IF (SWITCH) THEN
MINT = 2
MNTOLD = MINT
MITER = MTRSV
MTROLD = MITER
MAXORD = MIN(MXRDSV, 5)
NQ = MIN(NQ, MAXORD)
RC = 0.D0
RMAX = RMNORM
TREND = 1.D0
CALL DDCST (MAXORD, MINT, ISWFLG, EL, TQ)
NWAIT = NQ + 2
END IF
C Consider changing H if NWAIT = 1. Otherwise
C decrease NWAIT by 1. If NWAIT is then 1 and
C NQ.LT.MAXORD, then SAVE1 is saved for use in
C a possible order increase on the next step.
C
IF (JTASK .EQ. 0 .OR. JTASK .EQ. 2) THEN
RH = 1.D0/MAX(UROUND, BIAS2*(ETEST/EPS)**(1.D0/(NQ+1)))
IF (RH.GT.TRSHLD) CALL DDSCL (HMAX, N, NQ, RMAX, H, RC, RH, YH)
ELSE IF (NWAIT .GT. 1) THEN
NWAIT = NWAIT - 1
IF (NWAIT .EQ. 1 .AND. NQ .LT. MAXORD) THEN
DO 250 I = 1,NDE
250 YH(I,MAXORD+1) = SAVE1(I)
END IF
C If a change in H is considered, an increase or decrease in
C order by one is considered also. A change in H is made
C only if it is by a factor of at least TRSHLD. Factors
C RH1, RH2, and RH3 are computed, by which H could be
C multiplied at order NQ - 1, order NQ, or order NQ + 1,
C respectively. The largest of these is determined and the
C new order chosen accordingly. If the order is to be
C increased, we compute one additional scaled derivative.
C If there is a change of order, reset NQ and the
C coefficients. In any case H is reset according to RH and
C the YH array is rescaled.
ELSE
IF (NQ .EQ. 1) THEN
RH1 = 0.D0
ELSE
IF (IERROR .EQ. 1 .OR. IERROR .EQ. 5) THEN
DO 270 I = 1,NDE
270 SAVE2(I) = YH(I,NQ+1)/YWT(I)
ELSE
DO 275 I = 1,NDE
275 SAVE2(I) = YH(I,NQ+1)/MAX(ABS(Y(I)), YWT(I))
END IF
ERDN = DNRM2(NDE, SAVE2, 1)/(TQ(1,NQ)*SQRT(DBLE(NDE)))
RH1 = 1.D0/MAX(UROUND, BIAS1*(ERDN/EPS)**(1.D0/NQ))
END IF
RH2 = 1.D0/MAX(UROUND, BIAS2*(ETEST/EPS)**(1.D0/(NQ+1)))
IF (NQ .EQ. MAXORD) THEN
RH3 = 0.D0
ELSE
IF (IERROR .EQ. 1 .OR. IERROR .EQ. 5) THEN
DO 290 I = 1,NDE
290 SAVE2(I) = (SAVE1(I) - YH(I,MAXORD+1))/YWT(I)
ELSE
DO 295 I = 1,NDE
SAVE2(I) = (SAVE1(I) - YH(I,MAXORD+1))/
8 MAX(ABS(Y(I)), YWT(I))
295 CONTINUE
END IF
ERUP = DNRM2(NDE, SAVE2, 1)/(TQ(3,NQ)*SQRT(DBLE(NDE)))
RH3 = 1.D0/MAX(UROUND, BIAS3*(ERUP/EPS)**(1.D0/(NQ+2)))
END IF
IF (RH1 .GT. RH2 .AND. RH1 .GE. RH3) THEN
RH = RH1
IF (RH .LE. TRSHLD) GO TO 380
NQ = NQ - 1
RC = RC*EL(1,NQ)/EL(1,NQ+1)
ELSE IF (RH2 .GE. RH1 .AND. RH2 .GE. RH3) THEN
RH = RH2
IF (RH .LE. TRSHLD) GO TO 380
ELSE
RH = RH3
IF (RH .LE. TRSHLD) GO TO 380
DO 360 I = 1,N
360 YH(I,NQ+2) = SAVE1(I)*EL(NQ+1,NQ)/(NQ+1)
NQ = NQ + 1
RC = RC*EL(1,NQ)/EL(1,NQ-1)
END IF
IF (ISWFLG .EQ. 3 .AND. MINT .EQ. 1) THEN
IF (BND.NE.0.D0) RH = MIN(RH, 1.D0/(2.D0*EL(1,NQ)*BND*ABS(H)))
END IF
CALL DDSCL (HMAX, N, NQ, RMAX, H, RC, RH, YH)
RMAX = RMNORM
380 NWAIT = NQ + 2
END IF
C All returns are made through this section. H is saved
C in HOLD to allow the caller to change H on the next step
JSTATE = 1
HOLD = H
RETURN
C
400 JSTATE = 2
HOLD = H
DO 405 I = 1,N
405 Y(I) = YH(I,1)
RETURN
C
410 JSTATE = 3
HOLD = H
RETURN
C
420 JSTATE = 4
HOLD = H
RETURN
C
430 T = TOLD
CALL DDPSC (-1, NSV, NQ, YH)
DO 435 I = 1,NSV
435 Y(I) = YH(I,1)
440 HOLD = H
RETURN
END