Math-Lsoda
view release on metacpan or search on metacpan
C of the type of chord method used, or the Jacobian structure.
C Communication with DSTODE is done with the following variables:
C
C NEQ = integer array containing problem size in NEQ(1), and
C passed as the NEQ argument in all calls to F and JAC.
C Y = an array of length .ge. N used as the Y argument in
C all calls to F and JAC.
C YH = an NYH by LMAX array containing the dependent variables
C and their approximate scaled derivatives, where
C LMAX = MAXORD + 1. YH(i,j+1) contains the approximate
C j-th derivative of y(i), scaled by h**j/factorial(j)
C (j = 0,1,...,NQ). on entry for the first step, the first
C two columns of YH must be set from the initial values.
C NYH = a constant integer .ge. N, the first dimension of YH.
C YH1 = a one-dimensional array occupying the same space as YH.
C EWT = an array of length N containing multiplicative weights
C for local error measurements. Local errors in Y(i) are
C compared to 1.0/EWT(i) in various error tests.
C SAVF = an array of working storage, of length N.
C Also used for input of YH(*,MAXORD+2) when JSTART = -1
C and MAXORD .lt. the current order NQ.
C ACOR = a work array of length N, used for the accumulated
C corrections. On a successful return, ACOR(i) contains
C the estimated one-step local error in Y(i).
C WM,IWM = real and integer work arrays associated with matrix
C operations in chord iteration (MITER .ne. 0).
C PJAC = name of routine to evaluate and preprocess Jacobian matrix
C and P = I - h*el0*JAC, if a chord method is being used.
C SLVS = name of routine to solve linear system in chord iteration.
C CCMAX = maximum relative change in h*el0 before PJAC is called.
C H = the step size to be attempted on the next step.
C H is altered by the error control algorithm during the
C problem. H can be either positive or negative, but its
C sign must remain constant throughout the problem.
C HMIN = the minimum absolute value of the step size h to be used.
C HMXI = inverse of the maximum absolute value of h to be used.
C HMXI = 0.0 is allowed and corresponds to an infinite hmax.
C HMIN and HMXI may be changed at any time, but will not
C take effect until the next change of h is considered.
C TN = the independent variable. TN is updated on each step taken.
C JSTART = an integer used for input only, with the following
C values and meanings:
C 0 perform the first step.
C .gt.0 take a new step continuing from the last.
C -1 take the next step with a new value of H, MAXORD,
C N, METH, MITER, and/or matrix parameters.
C -2 take the next step with a new value of H,
C but with other inputs unchanged.
C On return, JSTART is set to 1 to facilitate continuation.
C KFLAG = a completion code with the following meanings:
C 0 the step was succesful.
C -1 the requested error could not be achieved.
C -2 corrector convergence could not be achieved.
C -3 fatal error in PJAC or SLVS.
C A return with KFLAG = -1 or -2 means either
C abs(H) = HMIN or 10 consecutive failures occurred.
C On a return with KFLAG negative, the values of TN and
C the YH array are as of the beginning of the last
C step, and H is the last step size attempted.
C MAXORD = the maximum order of integration method to be allowed.
C MAXCOR = the maximum number of corrector iterations allowed.
C MSBP = maximum number of steps between PJAC calls (MITER .gt. 0).
C MXNCF = maximum number of convergence failures allowed.
C METH/MITER = the method flags. See description in driver.
C N = the number of first-order differential equations.
C The values of CCMAX, H, HMIN, HMXI, TN, JSTART, KFLAG, MAXORD,
C MAXCOR, MSBP, MXNCF, METH, MITER, and N are communicated via COMMON.
C
C***SEE ALSO DLSODE
C***ROUTINES CALLED DCFODE, DVNORM
C***COMMON BLOCKS DLS001
C***REVISION HISTORY (YYMMDD)
C 791129 DATE WRITTEN
C 890501 Modified prologue to SLATEC/LDOC format. (FNF)
C 890503 Minor cosmetic changes. (FNF)
C 930809 Renamed to allow single/double precision versions. (ACH)
C 010418 Reduced size of Common block /DLS001/. (ACH)
C 031105 Restored 'own' variables to Common block /DLS001/, to
C enable interrupt/restart feature. (ACH)
C***END PROLOGUE DSTODE
C**End
EXTERNAL F, JAC, PJAC, SLVS
INTEGER NEQ, NYH, IWM
DOUBLE PRECISION Y, YH, YH1, EWT, SAVF, ACOR, WM
DIMENSION NEQ(*), Y(*), YH(NYH,*), YH1(*), EWT(*), SAVF(*),
1 ACOR(*), WM(*), IWM(*)
INTEGER IOWND, IALTH, IPUP, LMAX, MEO, NQNYH, NSLP,
1 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L,
2 LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER,
3 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
INTEGER I, I1, IREDO, IRET, J, JB, M, NCF, NEWQ
DOUBLE PRECISION CONIT, CRATE, EL, ELCO, HOLD, RMAX, TESCO,
2 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
DOUBLE PRECISION DCON, DDN, DEL, DELP, DSM, DUP, EXDN, EXSM, EXUP,
1 R, RH, RHDN, RHSM, RHUP, TOLD, DVNORM
COMMON /DLS001/ CONIT, CRATE, EL(13), ELCO(13,12),
1 HOLD, RMAX, TESCO(3,12),
2 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND,
3 IOWND(6), IALTH, IPUP, LMAX, MEO, NQNYH, NSLP,
3 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L,
4 LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER,
5 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
C
C***FIRST EXECUTABLE STATEMENT DSTODE
KFLAG = 0
TOLD = TN
NCF = 0
IERPJ = 0
IERSL = 0
JCUR = 0
ICF = 0
DELP = 0.0D0
IF (JSTART .GT. 0) GO TO 200
IF (JSTART .EQ. -1) GO TO 100
IF (JSTART .EQ. -2) GO TO 160
C-----------------------------------------------------------------------
C On the first call, the order is set to 1, and other variables are
C initialized. RMAX is the maximum ratio by which H can be increased
C in a single step. It is initially 1.E4 to compensate for the small
C initial H, but then is normally equal to 10. If a failure
C occurs (in corrector convergence or error test), RMAX is set to 2
IREDO = 3
IF (H .EQ. HOLD) GO TO 170
RH = MIN(RH,ABS(H/HOLD))
H = HOLD
GO TO 175
C-----------------------------------------------------------------------
C DCFODE is called to get all the integration coefficients for the
C current METH. Then the EL vector and related constants are reset
C whenever the order NQ is changed, or at the start of the problem.
C-----------------------------------------------------------------------
140 CALL DCFODE (METH, ELCO, TESCO)
150 DO 155 I = 1,L
155 EL(I) = ELCO(I,NQ)
NQNYH = NQ*NYH
RC = RC*EL(1)/EL0
EL0 = EL(1)
CONIT = 0.5D0/(NQ+2)
GO TO (160, 170, 200), IRET
C-----------------------------------------------------------------------
C If H is being changed, the H ratio RH is checked against
C RMAX, HMIN, and HMXI, and the YH array rescaled. IALTH is set to
C L = NQ + 1 to prevent a change of H for that many steps, unless
C forced by a convergence or error test failure.
C-----------------------------------------------------------------------
160 IF (H .EQ. HOLD) GO TO 200
RH = H/HOLD
H = HOLD
IREDO = 3
GO TO 175
170 RH = MAX(RH,HMIN/ABS(H))
175 RH = MIN(RH,RMAX)
RH = RH/MAX(1.0D0,ABS(H)*HMXI*RH)
R = 1.0D0
DO 180 J = 2,L
R = R*RH
DO 180 I = 1,N
180 YH(I,J) = YH(I,J)*R
H = H*RH
RC = RC*RH
IALTH = L
IF (IREDO .EQ. 0) GO TO 690
C-----------------------------------------------------------------------
C This section computes the predicted values by effectively
C multiplying the YH array by the Pascal Triangle matrix.
C RC is the ratio of new to old values of the coefficient H*EL(1).
C When RC differs from 1 by more than CCMAX, IPUP is set to MITER
C to force PJAC to be called, if a Jacobian is involved.
C In any case, PJAC is called at least every MSBP steps.
C-----------------------------------------------------------------------
200 IF (ABS(RC-1.0D0) .GT. CCMAX) IPUP = MITER
IF (NST .GE. NSLP+MSBP) IPUP = MITER
TN = TN + H
I1 = NQNYH + 1
DO 215 JB = 1,NQ
I1 = I1 - NYH
Cdir$ ivdep
DO 210 I = I1,NQNYH
210 YH1(I) = YH1(I) + YH1(I+NYH)
215 CONTINUE
C-----------------------------------------------------------------------
C Up to MAXCOR corrector iterations are taken. A convergence test is
C made on the R.M.S. norm of each correction, weighted by the error
C weight vector EWT. The sum of the corrections is accumulated in the
C vector ACOR(i). The YH array is not altered in the corrector loop.
C-----------------------------------------------------------------------
220 M = 0
DO 230 I = 1,N
230 Y(I) = YH(I,1)
CALL F (NEQ, TN, Y, SAVF)
NFE = NFE + 1
IF (IPUP .LE. 0) GO TO 250
C-----------------------------------------------------------------------
C If indicated, the matrix P = I - h*el(1)*J is reevaluated and
C preprocessed before starting the corrector iteration. IPUP is set
C to 0 as an indicator that this has been done.
C-----------------------------------------------------------------------
CALL PJAC (NEQ, Y, YH, NYH, EWT, ACOR, SAVF, WM, IWM, F, JAC)
IPUP = 0
RC = 1.0D0
NSLP = NST
CRATE = 0.7D0
IF (IERPJ .NE. 0) GO TO 430
250 DO 260 I = 1,N
260 ACOR(I) = 0.0D0
270 IF (MITER .NE. 0) GO TO 350
C-----------------------------------------------------------------------
C In the case of functional iteration, update Y directly from
C the result of the last function evaluation.
C-----------------------------------------------------------------------
DO 290 I = 1,N
SAVF(I) = H*SAVF(I) - YH(I,2)
290 Y(I) = SAVF(I) - ACOR(I)
DEL = DVNORM (N, Y, EWT)
DO 300 I = 1,N
Y(I) = YH(I,1) + EL(1)*SAVF(I)
300 ACOR(I) = SAVF(I)
GO TO 400
C-----------------------------------------------------------------------
C In the case of the chord method, compute the corrector error,
C and solve the linear system with that as right-hand side and
C P as coefficient matrix.
C-----------------------------------------------------------------------
350 DO 360 I = 1,N
360 Y(I) = H*SAVF(I) - (YH(I,2) + ACOR(I))
CALL SLVS (WM, IWM, Y, SAVF)
IF (IERSL .LT. 0) GO TO 430
IF (IERSL .GT. 0) GO TO 410
DEL = DVNORM (N, Y, EWT)
DO 380 I = 1,N
ACOR(I) = ACOR(I) + Y(I)
380 Y(I) = YH(I,1) + EL(1)*ACOR(I)
C-----------------------------------------------------------------------
C Test for convergence. If M.gt.0, an estimate of the convergence
C rate constant is stored in CRATE, and this is used in the test.
C-----------------------------------------------------------------------
400 IF (M .NE. 0) CRATE = MAX(0.2D0*CRATE,DEL/DELP)
DCON = DEL*MIN(1.0D0,1.5D0*CRATE)/(TESCO(2,NQ)*CONIT)
IF (DCON .LE. 1.0D0) GO TO 450
M = M + 1
IF (M .EQ. MAXCOR) GO TO 410
IF (M .GE. 2 .AND. DEL .GT. 2.0D0*DELP) GO TO 410
C indicator MITER, when this is .ne. 0, and hence is independent
C of the type of chord method used, or the Jacobian structure.
C Communication with DSTODA is done with the following variables:
C
C Y = an array of length .ge. N used as the Y argument in
C all calls to F and JAC.
C NEQ = integer array containing problem size in NEQ(1), and
C passed as the NEQ argument in all calls to F and JAC.
C YH = an NYH by LMAX array containing the dependent variables
C and their approximate scaled derivatives, where
C LMAX = MAXORD + 1. YH(i,j+1) contains the approximate
C j-th derivative of y(i), scaled by H**j/factorial(j)
C (j = 0,1,...,NQ). On entry for the first step, the first
C two columns of YH must be set from the initial values.
C NYH = a constant integer .ge. N, the first dimension of YH.
C YH1 = a one-dimensional array occupying the same space as YH.
C EWT = an array of length N containing multiplicative weights
C for local error measurements. Local errors in y(i) are
C compared to 1.0/EWT(i) in various error tests.
C SAVF = an array of working storage, of length N.
C ACOR = a work array of length N, used for the accumulated
C corrections. On a successful return, ACOR(i) contains
C the estimated one-step local error in y(i).
C WM,IWM = real and integer work arrays associated with matrix
C operations in chord iteration (MITER .ne. 0).
C PJAC = name of routine to evaluate and preprocess Jacobian matrix
C and P = I - H*EL0*Jac, if a chord method is being used.
C It also returns an estimate of norm(Jac) in PDNORM.
C SLVS = name of routine to solve linear system in chord iteration.
C CCMAX = maximum relative change in H*EL0 before PJAC is called.
C H = the step size to be attempted on the next step.
C H is altered by the error control algorithm during the
C problem. H can be either positive or negative, but its
C sign must remain constant throughout the problem.
C HMIN = the minimum absolute value of the step size H to be used.
C HMXI = inverse of the maximum absolute value of H to be used.
C HMXI = 0.0 is allowed and corresponds to an infinite HMAX.
C HMIN and HMXI may be changed at any time, but will not
C take effect until the next change of H is considered.
C TN = the independent variable. TN is updated on each step taken.
C JSTART = an integer used for input only, with the following
C values and meanings:
C 0 perform the first step.
C .gt.0 take a new step continuing from the last.
C -1 take the next step with a new value of H,
C N, METH, MITER, and/or matrix parameters.
C -2 take the next step with a new value of H,
C but with other inputs unchanged.
C On return, JSTART is set to 1 to facilitate continuation.
C KFLAG = a completion code with the following meanings:
C 0 the step was succesful.
C -1 the requested error could not be achieved.
C -2 corrector convergence could not be achieved.
C -3 fatal error in PJAC or SLVS.
C A return with KFLAG = -1 or -2 means either
C ABS(H) = HMIN or 10 consecutive failures occurred.
C On a return with KFLAG negative, the values of TN and
C the YH array are as of the beginning of the last
C step, and H is the last step size attempted.
C MAXORD = the maximum order of integration method to be allowed.
C MAXCOR = the maximum number of corrector iterations allowed.
C MSBP = maximum number of steps between PJAC calls (MITER .gt. 0).
C MXNCF = maximum number of convergence failures allowed.
C METH = current method.
C METH = 1 means Adams method (nonstiff)
C METH = 2 means BDF method (stiff)
C METH may be reset by DSTODA.
C MITER = corrector iteration method.
C MITER = 0 means functional iteration.
C MITER = JT .gt. 0 means a chord iteration corresponding
C to Jacobian type JT. (The DLSODA/DLSODAR argument JT is
C communicated here as JTYP, but is not used in DSTODA
C except to load MITER following a method switch.)
C MITER may be reset by DSTODA.
C N = the number of first-order differential equations.
C-----------------------------------------------------------------------
KFLAG = 0
TOLD = TN
NCF = 0
IERPJ = 0
IERSL = 0
JCUR = 0
ICF = 0
DELP = 0.0D0
IF (JSTART .GT. 0) GO TO 200
IF (JSTART .EQ. -1) GO TO 100
IF (JSTART .EQ. -2) GO TO 160
C-----------------------------------------------------------------------
C On the first call, the order is set to 1, and other variables are
C initialized. RMAX is the maximum ratio by which H can be increased
C in a single step. It is initially 1.E4 to compensate for the small
C initial H, but then is normally equal to 10. If a failure
C occurs (in corrector convergence or error test), RMAX is set at 2
C for the next increase.
C DCFODE is called to get the needed coefficients for both methods.
C-----------------------------------------------------------------------
LMAX = MAXORD + 1
NQ = 1
L = 2
IALTH = 2
RMAX = 10000.0D0
RC = 0.0D0
EL0 = 1.0D0
CRATE = 0.7D0
HOLD = H
NSLP = 0
IPUP = MITER
IRET = 3
C Initialize switching parameters. METH = 1 is assumed initially. -----
ICOUNT = 20
IRFLAG = 0
PDEST = 0.0D0
PDLAST = 0.0D0
RATIO = 5.0D0
CALL DCFODE (2, ELCO, TESCO)
DO 10 I = 1,5
10 CM2(I) = TESCO(2,I)*ELCO(I+1,I)
CALL DCFODE (1, ELCO, TESCO)
DO 20 I = 1,12
20 CM1(I) = TESCO(2,I)*ELCO(I+1,I)
GO TO 150
NQNYH = NQ*NYH
RC = RC*EL(1)/EL0
EL0 = EL(1)
CONIT = 0.5D0/(NQ+2)
GO TO (160, 170, 200), IRET
C-----------------------------------------------------------------------
C If H is being changed, the H ratio RH is checked against
C RMAX, HMIN, and HMXI, and the YH array rescaled. IALTH is set to
C L = NQ + 1 to prevent a change of H for that many steps, unless
C forced by a convergence or error test failure.
C-----------------------------------------------------------------------
160 IF (H .EQ. HOLD) GO TO 200
RH = H/HOLD
H = HOLD
IREDO = 3
GO TO 175
170 RH = MAX(RH,HMIN/ABS(H))
175 RH = MIN(RH,RMAX)
RH = RH/MAX(1.0D0,ABS(H)*HMXI*RH)
C-----------------------------------------------------------------------
C If METH = 1, also restrict the new step size by the stability region.
C If this reduces H, set IRFLAG to 1 so that if there are roundoff
C problems later, we can assume that is the cause of the trouble.
C-----------------------------------------------------------------------
IF (METH .EQ. 2) GO TO 178
IRFLAG = 0
PDH = MAX(ABS(H)*PDLAST,0.000001D0)
IF (RH*PDH*1.00001D0 .LT. SM1(NQ)) GO TO 178
RH = SM1(NQ)/PDH
IRFLAG = 1
178 CONTINUE
R = 1.0D0
DO 180 J = 2,L
R = R*RH
DO 180 I = 1,N
180 YH(I,J) = YH(I,J)*R
H = H*RH
RC = RC*RH
IALTH = L
IF (IREDO .EQ. 0) GO TO 690
C-----------------------------------------------------------------------
C This section computes the predicted values by effectively
C multiplying the YH array by the Pascal triangle matrix.
C RC is the ratio of new to old values of the coefficient H*EL(1).
C When RC differs from 1 by more than CCMAX, IPUP is set to MITER
C to force PJAC to be called, if a Jacobian is involved.
C In any case, PJAC is called at least every MSBP steps.
C-----------------------------------------------------------------------
200 IF (ABS(RC-1.0D0) .GT. CCMAX) IPUP = MITER
IF (NST .GE. NSLP+MSBP) IPUP = MITER
TN = TN + H
I1 = NQNYH + 1
DO 215 JB = 1,NQ
I1 = I1 - NYH
CDIR$ IVDEP
DO 210 I = I1,NQNYH
210 YH1(I) = YH1(I) + YH1(I+NYH)
215 CONTINUE
PNORM = DMNORM (N, YH1, EWT)
C-----------------------------------------------------------------------
C Up to MAXCOR corrector iterations are taken. A convergence test is
C made on the RMS-norm of each correction, weighted by the error
C weight vector EWT. The sum of the corrections is accumulated in the
C vector ACOR(i). The YH array is not altered in the corrector loop.
C-----------------------------------------------------------------------
220 M = 0
RATE = 0.0D0
DEL = 0.0D0
DO 230 I = 1,N
230 Y(I) = YH(I,1)
CALL F (NEQ, TN, Y, SAVF)
NFE = NFE + 1
IF (IPUP .LE. 0) GO TO 250
C-----------------------------------------------------------------------
C If indicated, the matrix P = I - H*EL(1)*J is reevaluated and
C preprocessed before starting the corrector iteration. IPUP is set
C to 0 as an indicator that this has been done.
C-----------------------------------------------------------------------
CALL PJAC (NEQ, Y, YH, NYH, EWT, ACOR, SAVF, WM, IWM, F, JAC)
IPUP = 0
RC = 1.0D0
NSLP = NST
CRATE = 0.7D0
IF (IERPJ .NE. 0) GO TO 430
250 DO 260 I = 1,N
260 ACOR(I) = 0.0D0
270 IF (MITER .NE. 0) GO TO 350
C-----------------------------------------------------------------------
C In the case of functional iteration, update Y directly from
C the result of the last function evaluation.
C-----------------------------------------------------------------------
DO 290 I = 1,N
SAVF(I) = H*SAVF(I) - YH(I,2)
290 Y(I) = SAVF(I) - ACOR(I)
DEL = DMNORM (N, Y, EWT)
DO 300 I = 1,N
Y(I) = YH(I,1) + EL(1)*SAVF(I)
300 ACOR(I) = SAVF(I)
GO TO 400
C-----------------------------------------------------------------------
C In the case of the chord method, compute the corrector error,
C and solve the linear system with that as right-hand side and
C P as coefficient matrix.
C-----------------------------------------------------------------------
350 DO 360 I = 1,N
360 Y(I) = H*SAVF(I) - (YH(I,2) + ACOR(I))
CALL SLVS (WM, IWM, Y, SAVF)
IF (IERSL .LT. 0) GO TO 430
IF (IERSL .GT. 0) GO TO 410
DEL = DMNORM (N, Y, EWT)
DO 380 I = 1,N
ACOR(I) = ACOR(I) + Y(I)
380 Y(I) = YH(I,1) + EL(1)*ACOR(I)
C-----------------------------------------------------------------------
C Test for convergence. If M .gt. 0, an estimate of the convergence
C rate constant is stored in CRATE, and this is used in the test.
C
C We first check for a change of iterates that is the size of
C roundoff error. If this occurs, the iteration has converged, and a
C new rate estimate is not formed.
C In all other cases, force at least two iterations to estimate a
C local Lipschitz constant estimate for Adams methods.
C On convergence, form PDEST = local maximum Lipschitz constant
C estimate. PDLAST is the most recent nonzero estimate.
C-----------------------------------------------------------------------
400 CONTINUE
IF (DEL .LE. 100.0D0*PNORM*UROUND) GO TO 450
IF (M .EQ. 0 .AND. METH .EQ. 1) GO TO 405
IF (M .EQ. 0) GO TO 402
RM = 1024.0D0
IF (DEL .LE. 1024.0D0*DELP) RM = DEL/DELP
RATE = MAX(RATE,RM)
CRATE = MAX(0.2D0*CRATE,RM)
402 DCON = DEL*MIN(1.0D0,1.5D0*CRATE)/(TESCO(2,NQ)*CONIT)
IF (DCON .GT. 1.0D0) GO TO 405
PDEST = MAX(PDEST,RATE/ABS(H*EL(1)))
IF (PDEST .NE. 0.0D0) PDLAST = PDEST
GO TO 450
405 CONTINUE
M = M + 1
IF (M .EQ. MAXCOR) GO TO 410
IF (M .GE. 2 .AND. DEL .GT. 2.0D0*DELP) GO TO 410
DELP = DEL
CALL F (NEQ, TN, Y, SAVF)
NFE = NFE + 1
GO TO 270
C-----------------------------------------------------------------------
C The corrector iteration failed to converge.
C If MITER .ne. 0 and the Jacobian is out of date, PJAC is called for
C the next try. Otherwise the YH array is retracted to its values
C before prediction, and H is reduced, if possible. If H cannot be
C reduced or MXNCF failures have occurred, exit with KFLAG = -2.
C-----------------------------------------------------------------------
410 IF (MITER .EQ. 0 .OR. JCUR .EQ. 1) GO TO 430
ICF = 1
IPUP = MITER
GO TO 220
430 ICF = 2
NCF = NCF + 1
RMAX = 2.0D0
TN = TOLD
I1 = NQNYH + 1
DO 445 JB = 1,NQ
I1 = I1 - NYH
CDIR$ IVDEP
DO 440 I = I1,NQNYH
440 YH1(I) = YH1(I) - YH1(I+NYH)
445 CONTINUE
IF (IERPJ .LT. 0 .OR. IERSL .LT. 0) GO TO 680
IF (ABS(H) .LE. HMIN*1.00001D0) GO TO 670
IF (NCF .EQ. MXNCF) GO TO 670
RH = 0.25D0
IPUP = MITER
IREDO = 1
GO TO 170
C-----------------------------------------------------------------------
C The corrector has converged. JCUR is set to 0
C to signal that the Jacobian involved may need updating later.
C The local error test is made and control passes to statement 500
C if it fails.
C-----------------------------------------------------------------------
C Note: DSTODPK is independent of the value of the iteration method
C indicator MITER, when this is .ne. 0, and hence is independent
C of the type of chord method used, or the Jacobian structure.
C Communication with DSTODPK is done with the following variables:
C
C NEQ = integer array containing problem size in NEQ(1), and
C passed as the NEQ argument in all calls to F and JAC.
C Y = an array of length .ge. N used as the Y argument in
C all calls to F and JAC.
C YH = an NYH by LMAX array containing the dependent variables
C and their approximate scaled derivatives, where
C LMAX = MAXORD + 1. YH(i,j+1) contains the approximate
C j-th derivative of y(i), scaled by H**j/factorial(j)
C (j = 0,1,...,NQ). On entry for the first step, the first
C two columns of YH must be set from the initial values.
C NYH = a constant integer .ge. N, the first dimension of YH.
C YH1 = a one-dimensional array occupying the same space as YH.
C EWT = an array of length N containing multiplicative weights
C for local error measurements. Local errors in y(i) are
C compared to 1.0/EWT(i) in various error tests.
C SAVF = an array of working storage, of length N.
C Also used for input of YH(*,MAXORD+2) when JSTART = -1
C and MAXORD .lt. the current order NQ.
C SAVX = an array of working storage, of length N.
C ACOR = a work array of length N, used for the accumulated
C corrections. On a successful return, ACOR(i) contains
C the estimated one-step local error in y(i).
C WM,IWM = real and integer work arrays associated with matrix
C operations in chord iteration (MITER .ne. 0).
C CCMAX = maximum relative change in H*EL0 before DPKSET is called.
C H = the step size to be attempted on the next step.
C H is altered by the error control algorithm during the
C problem. H can be either positive or negative, but its
C sign must remain constant throughout the problem.
C HMIN = the minimum absolute value of the step size H to be used.
C HMXI = inverse of the maximum absolute value of H to be used.
C HMXI = 0.0 is allowed and corresponds to an infinite HMAX.
C HMIN and HMXI may be changed at any time, but will not
C take effect until the next change of H is considered.
C TN = the independent variable. TN is updated on each step taken.
C JSTART = an integer used for input only, with the following
C values and meanings:
C 0 perform the first step.
C .gt.0 take a new step continuing from the last.
C -1 take the next step with a new value of H, MAXORD,
C N, METH, MITER, and/or matrix parameters.
C -2 take the next step with a new value of H,
C but with other inputs unchanged.
C On return, JSTART is set to 1 to facilitate continuation.
C KFLAG = a completion code with the following meanings:
C 0 the step was succesful.
C -1 the requested error could not be achieved.
C -2 corrector convergence could not be achieved.
C -3 fatal error in DPKSET or DSOLPK.
C A return with KFLAG = -1 or -2 means either
C ABS(H) = HMIN or 10 consecutive failures occurred.
C On a return with KFLAG negative, the values of TN and
C the YH array are as of the beginning of the last
C step, and H is the last step size attempted.
C MAXORD = the maximum order of integration method to be allowed.
C MAXCOR = the maximum number of corrector iterations allowed.
C MSBP = maximum number of steps between DPKSET calls (MITER .gt. 0).
C MXNCF = maximum number of convergence failures allowed.
C METH/MITER = the method flags. See description in driver.
C N = the number of first-order differential equations.
C-----------------------------------------------------------------------
INTEGER I, I1, IREDO, IRET, J, JB, M, NCF, NEWQ
DOUBLE PRECISION DCON, DDN, DEL, DELP, DSM, DUP, EXDN, EXSM, EXUP,
1 R, RH, RHDN, RHSM, RHUP, TOLD, DVNORM
C
KFLAG = 0
TOLD = TN
NCF = 0
IERPJ = 0
IERSL = 0
JCUR = 0
ICF = 0
DELP = 0.0D0
IF (JSTART .GT. 0) GO TO 200
IF (JSTART .EQ. -1) GO TO 100
IF (JSTART .EQ. -2) GO TO 160
C-----------------------------------------------------------------------
C On the first call, the order is set to 1, and other variables are
C initialized. RMAX is the maximum ratio by which H can be increased
C in a single step. It is initially 1.E4 to compensate for the small
C initial H, but then is normally equal to 10. If a failure
C occurs (in corrector convergence or error test), RMAX is set at 2
C for the next increase.
C-----------------------------------------------------------------------
LMAX = MAXORD + 1
NQ = 1
L = 2
IALTH = 2
RMAX = 10000.0D0
RC = 0.0D0
EL0 = 1.0D0
CRATE = 0.7D0
HOLD = H
MEO = METH
NSLP = 0
IPUP = MITER
IRET = 3
GO TO 140
C-----------------------------------------------------------------------
C The following block handles preliminaries needed when JSTART = -1.
C IPUP is set to MITER to force a matrix update.
C If an order increase is about to be considered (IALTH = 1),
C IALTH is reset to 2 to postpone consideration one more step.
C If the caller has changed METH, DCFODE is called to reset
C the coefficients of the method.
C If the caller has changed MAXORD to a value less than the current
C order NQ, NQ is reduced to MAXORD, and a new H chosen accordingly.
C If H is to be changed, YH must be rescaled.
C If H or METH is being changed, IALTH is reset to L = NQ + 1
C to prevent further changes in H for that many steps.
C-----------------------------------------------------------------------
100 IPUP = MITER
LMAX = MAXORD + 1
IF (IALTH .EQ. 1) IALTH = 2
IF (METH .EQ. MEO) GO TO 110
CALL DCFODE (METH, ELCO, TESCO)
C DCFODE is called to get all the integration coefficients for the
C current METH. Then the EL vector and related constants are reset
C whenever the order NQ is changed, or at the start of the problem.
C-----------------------------------------------------------------------
140 CALL DCFODE (METH, ELCO, TESCO)
150 DO 155 I = 1,L
155 EL(I) = ELCO(I,NQ)
NQNYH = NQ*NYH
RC = RC*EL(1)/EL0
EL0 = EL(1)
CONIT = 0.5D0/(NQ+2)
EPCON = CONIT*TESCO(2,NQ)
GO TO (160, 170, 200), IRET
C-----------------------------------------------------------------------
C If H is being changed, the H ratio RH is checked against
C RMAX, HMIN, and HMXI, and the YH array rescaled. IALTH is set to
C L = NQ + 1 to prevent a change of H for that many steps, unless
C forced by a convergence or error test failure.
C-----------------------------------------------------------------------
160 IF (H .EQ. HOLD) GO TO 200
RH = H/HOLD
H = HOLD
IREDO = 3
GO TO 175
170 RH = MAX(RH,HMIN/ABS(H))
175 RH = MIN(RH,RMAX)
RH = RH/MAX(1.0D0,ABS(H)*HMXI*RH)
R = 1.0D0
DO 180 J = 2,L
R = R*RH
DO 180 I = 1,N
180 YH(I,J) = YH(I,J)*R
H = H*RH
RC = RC*RH
IALTH = L
IF (IREDO .EQ. 0) GO TO 690
C-----------------------------------------------------------------------
C This section computes the predicted values by effectively
C multiplying the YH array by the Pascal triangle matrix.
C The flag IPUP is set according to whether matrix data is involved
C (JACFLG .ne. 0) or not (JACFLG = 0), to trigger a call to DPKSET.
C IPUP is set to MITER when RC differs from 1 by more than CCMAX,
C and at least every MSBP steps, when JACFLG = 1.
C RC is the ratio of new to old values of the coefficient H*EL(1).
C-----------------------------------------------------------------------
200 IF (JACFLG .NE. 0) GO TO 202
IPUP = 0
CRATE = 0.7D0
GO TO 205
202 IF (ABS(RC-1.0D0) .GT. CCMAX) IPUP = MITER
IF (NST .GE. NSLP+MSBP) IPUP = MITER
205 TN = TN + H
I1 = NQNYH + 1
DO 215 JB = 1,NQ
I1 = I1 - NYH
CDIR$ IVDEP
DO 210 I = I1,NQNYH
210 YH1(I) = YH1(I) + YH1(I+NYH)
215 CONTINUE
C-----------------------------------------------------------------------
C Up to MAXCOR corrector iterations are taken. A convergence test is
C made on the RMS-norm of each correction, weighted by the error
C weight vector EWT. The sum of the corrections is accumulated in the
C vector ACOR(i). The YH array is not altered in the corrector loop.
C-----------------------------------------------------------------------
220 M = 0
MNEWT = 0
DO 230 I = 1,N
230 Y(I) = YH(I,1)
CALL F (NEQ, TN, Y, SAVF)
NFE = NFE + 1
IF (IPUP .LE. 0) GO TO 250
C-----------------------------------------------------------------------
C If indicated, DPKSET is called to update any matrix data needed,
C before starting the corrector iteration.
C IPUP is set to 0 as an indicator that this has been done.
C-----------------------------------------------------------------------
CALL DPKSET (NEQ, Y, YH1, EWT, ACOR, SAVF, WM, IWM, F, JAC)
IPUP = 0
RC = 1.0D0
NSLP = NST
CRATE = 0.7D0
IF (IERPJ .NE. 0) GO TO 430
250 DO 260 I = 1,N
260 ACOR(I) = 0.0D0
270 IF (MITER .NE. 0) GO TO 350
C-----------------------------------------------------------------------
C In the case of functional iteration, update Y directly from
C the result of the last function evaluation.
C-----------------------------------------------------------------------
DO 290 I = 1,N
SAVF(I) = H*SAVF(I) - YH(I,2)
290 Y(I) = SAVF(I) - ACOR(I)
DEL = DVNORM (N, Y, EWT)
DO 300 I = 1,N
Y(I) = YH(I,1) + EL(1)*SAVF(I)
300 ACOR(I) = SAVF(I)
GO TO 400
C-----------------------------------------------------------------------
C In the case of the chord method, compute the corrector error,
C and solve the linear system with that as right-hand side and
C P as coefficient matrix.
C-----------------------------------------------------------------------
350 DO 360 I = 1,N
360 SAVX(I) = H*SAVF(I) - (YH(I,2) + ACOR(I))
CALL DSOLPK (NEQ, Y, SAVF, SAVX, EWT, WM, IWM, F, PSOL)
IF (IERSL .LT. 0) GO TO 430
IF (IERSL .GT. 0) GO TO 410
DEL = DVNORM (N, SAVX, EWT)
DO 380 I = 1,N
ACOR(I) = ACOR(I) + SAVX(I)
380 Y(I) = YH(I,1) + EL(1)*ACOR(I)
C-----------------------------------------------------------------------
C Test for convergence. If M .gt. 0, an estimate of the convergence
C rate constant is stored in CRATE, and this is used in the test.
C-----------------------------------------------------------------------
400 IF (M .NE. 0) CRATE = MAX(0.2D0*CRATE,DEL/DELP)
DCON = DEL*MIN(1.0D0,1.5D0*CRATE)/EPCON
IF (DCON .LE. 1.0D0) GO TO 450
M = M + 1
IF (M .EQ. MAXCOR) GO TO 410
C preconditioned version of the Incomplete Orthogonalization Method.
C An initial guess of x = 0 is assumed.
C-----------------------------------------------------------------------
C
C On entry
C
C NEQ = problem size, passed to F and PSOL (NEQ(1) = N).
C
C TN = current value of t.
C
C Y = array containing current dependent variable vector.
C
C SAVF = array containing current value of f(t,y).
C
C B = the right hand side of the system A*x = b.
C B is also used as work space when computing the
C final approximation.
C (B is the same as V(*,MAXL+1) in the call to DSPIOM.)
C
C WGHT = array of length N containing scale factors.
C 1/WGHT(i) are the diagonal elements of the diagonal
C scaling matrix D.
C
C N = the order of the matrix A, and the lengths
C of the vectors Y, SAVF, B, WGHT, and X.
C
C MAXL = the maximum allowable order of the matrix HES.
C
C KMP = the number of previous vectors the new vector VNEW
C must be made orthogonal to. KMP .le. MAXL.
C
C DELTA = tolerance on residuals b - A*x in weighted RMS-norm.
C
C HL0 = current value of (step size h) * (coefficient l0).
C
C JPRE = preconditioner type flag.
C
C MNEWT = Newton iteration counter (.ge. 0).
C
C WK = real work array of length N used by DATV and PSOL.
C
C WP = real work array used by preconditioner PSOL.
C
C IWP = integer work array used by preconditioner PSOL.
C
C On return
C
C X = the final computed approximation to the solution
C of the system A*x = b.
C
C V = the N by (LIOM+1) array containing the LIOM
C orthogonal vectors V(*,1) to V(*,LIOM).
C
C HES = the LU factorization of the LIOM by LIOM upper
C Hessenberg matrix whose entries are the
C scaled inner products of A*V(*,k) and V(*,i).
C
C IPVT = an integer array containg pivoting information.
C It is loaded in DHEFA and used in DHESL.
C
C LIOM = the number of iterations performed, and current
C order of the upper Hessenberg matrix HES.
C
C NPSL = the number of calls to PSOL.
C
C IFLAG = integer error flag:
C 0 means convergence in LIOM iterations, LIOM.le.MAXL.
C 1 means the convergence test did not pass in MAXL
C iterations, but the residual norm is .lt. 1,
C or .lt. norm(b) if MNEWT = 0, and so X is computed.
C 2 means the convergence test did not pass in MAXL
C iterations, residual .gt. 1, and X is undefined.
C 3 means there was a recoverable error in PSOL
C caused by the preconditioner being out of date.
C -1 means there was a nonrecoverable error in PSOL.
C
C-----------------------------------------------------------------------
INTEGER I, IER, INFO, J, K, LL, LM1
DOUBLE PRECISION BNRM, BNRM0, PROD, RHO, SNORMW, DNRM2, TEM
C
IFLAG = 0
LIOM = 0
NPSL = 0
C-----------------------------------------------------------------------
C The initial residual is the vector b. Apply scaling to b, and test
C for an immediate return with X = 0 or X = b.
C-----------------------------------------------------------------------
DO 10 I = 1,N
10 V(I,1) = B(I)*WGHT(I)
BNRM0 = DNRM2 (N, V, 1)
BNRM = BNRM0
IF (BNRM0 .GT. DELTA) GO TO 30
IF (MNEWT .GT. 0) GO TO 20
CALL DCOPY (N, B, 1, X, 1)
RETURN
20 DO 25 I = 1,N
25 X(I) = 0.0D0
RETURN
30 CONTINUE
C Apply inverse of left preconditioner to vector b. --------------------
IER = 0
IF (JPRE .EQ. 0 .OR. JPRE .EQ. 2) GO TO 55
CALL PSOL (NEQ, TN, Y, SAVF, WK, HL0, WP, IWP, B, 1, IER)
NPSL = 1
IF (IER .NE. 0) GO TO 300
C Calculate norm of scaled vector V(*,1) and normalize it. -------------
DO 50 I = 1,N
50 V(I,1) = B(I)*WGHT(I)
BNRM = DNRM2(N, V, 1)
DELTA = DELTA*(BNRM/BNRM0)
55 TEM = 1.0D0/BNRM
CALL DSCAL (N, TEM, V(1,1), 1)
C Zero out the HES array. ----------------------------------------------
DO 65 J = 1,MAXL
DO 60 I = 1,MAXL
60 HES(I,J) = 0.0D0
65 CONTINUE
C-----------------------------------------------------------------------
C Main loop on LL = l to compute the vectors V(*,2) to V(*,MAXL).
C The running product PROD is needed for the convergence test.
C-----------------------------------------------------------------------
PROD = 1.0D0
DO 90 LL = 1,MAXL
LIOM = LL
C-----------------------------------------------------------------------
C Call routine DATV to compute VNEW = Abar*v(l), where Abar is
C the matrix A with scaling and inverse preconditioner factors applied.
C Call routine DORTHOG to orthogonalize the new vector vnew = V(*,l+1).
C Call routine DHEFA to update the factors of HES.
C-----------------------------------------------------------------------
CALL DATV (NEQ, Y, SAVF, V(1,LL), WGHT, X, F, PSOL, V(1,LL+1),
1 WK, WP, IWP, HL0, JPRE, IER, NPSL)
DOUBLE PRECISION TN,Y,SAVF,B,WGHT,DELTA,HL0,X,V,HES,Q,WP,WK,DL
DIMENSION NEQ(*), Y(*), SAVF(*), B(*), WGHT(*), X(*), V(N,*),
1 HES(MAXLP1,*), Q(*), WP(*), IWP(*), WK(*), DL(*)
C-----------------------------------------------------------------------
C This routine solves the linear system A * x = b using a scaled
C preconditioned version of the Generalized Minimal Residual method.
C An initial guess of x = 0 is assumed.
C-----------------------------------------------------------------------
C
C On entry
C
C NEQ = problem size, passed to F and PSOL (NEQ(1) = N).
C
C TN = current value of t.
C
C Y = array containing current dependent variable vector.
C
C SAVF = array containing current value of f(t,y).
C
C B = the right hand side of the system A*x = b.
C B is also used as work space when computing
C the final approximation.
C (B is the same as V(*,MAXL+1) in the call to DSPIGMR.)
C
C WGHT = the vector of length N containing the nonzero
C elements of the diagonal scaling matrix.
C
C N = the order of the matrix A, and the lengths
C of the vectors WGHT, B and X.
C
C MAXL = the maximum allowable order of the matrix HES.
C
C MAXLP1 = MAXL + 1, used for dynamic dimensioning of HES.
C
C KMP = the number of previous vectors the new vector VNEW
C must be made orthogonal to. KMP .le. MAXL.
C
C DELTA = tolerance on residuals b - A*x in weighted RMS-norm.
C
C HL0 = current value of (step size h) * (coefficient l0).
C
C JPRE = preconditioner type flag.
C
C MNEWT = Newton iteration counter (.ge. 0).
C
C WK = real work array used by routine DATV and PSOL.
C
C DL = real work array used for calculation of the residual
C norm RHO when the method is incomplete (KMP .lt. MAXL).
C Not needed or referenced in complete case (KMP = MAXL).
C
C WP = real work array used by preconditioner PSOL.
C
C IWP = integer work array used by preconditioner PSOL.
C
C On return
C
C X = the final computed approximation to the solution
C of the system A*x = b.
C
C LGMR = the number of iterations performed and
C the current order of the upper Hessenberg
C matrix HES.
C
C NPSL = the number of calls to PSOL.
C
C V = the N by (LGMR+1) array containing the LGMR
C orthogonal vectors V(*,1) to V(*,LGMR).
C
C HES = the upper triangular factor of the QR decomposition
C of the (LGMR+1) by lgmr upper Hessenberg matrix whose
C entries are the scaled inner-products of A*V(*,i)
C and V(*,k).
C
C Q = real array of length 2*MAXL containing the components
C of the Givens rotations used in the QR decomposition
C of HES. It is loaded in DHEQR and used in DHELS.
C
C IFLAG = integer error flag:
C 0 means convergence in LGMR iterations, LGMR .le. MAXL.
C 1 means the convergence test did not pass in MAXL
C iterations, but the residual norm is .lt. 1,
C or .lt. norm(b) if MNEWT = 0, and so x is computed.
C 2 means the convergence test did not pass in MAXL
C iterations, residual .gt. 1, and X is undefined.
C 3 means there was a recoverable error in PSOL
C caused by the preconditioner being out of date.
C -1 means there was a nonrecoverable error in PSOL.
C
C-----------------------------------------------------------------------
INTEGER I, IER, INFO, IP1, I2, J, K, LL, LLP1
DOUBLE PRECISION BNRM,BNRM0,C,DLNRM,PROD,RHO,S,SNORMW,DNRM2,TEM
C
IFLAG = 0
LGMR = 0
NPSL = 0
C-----------------------------------------------------------------------
C The initial residual is the vector b. Apply scaling to b, and test
C for an immediate return with X = 0 or X = b.
C-----------------------------------------------------------------------
DO 10 I = 1,N
10 V(I,1) = B(I)*WGHT(I)
BNRM0 = DNRM2 (N, V, 1)
BNRM = BNRM0
IF (BNRM0 .GT. DELTA) GO TO 30
IF (MNEWT .GT. 0) GO TO 20
CALL DCOPY (N, B, 1, X, 1)
RETURN
20 DO 25 I = 1,N
25 X(I) = 0.0D0
RETURN
30 CONTINUE
C Apply inverse of left preconditioner to vector b. --------------------
IER = 0
IF (JPRE .EQ. 0 .OR. JPRE .EQ. 2) GO TO 55
CALL PSOL (NEQ, TN, Y, SAVF, WK, HL0, WP, IWP, B, 1, IER)
NPSL = 1
IF (IER .NE. 0) GO TO 300
C Calculate norm of scaled vector V(*,1) and normalize it. -------------
DO 50 I = 1,N
50 V(I,1) = B(I)*WGHT(I)
BNRM = DNRM2 (N, V, 1)
DELTA = DELTA*(BNRM/BNRM0)
55 TEM = 1.0D0/BNRM
CALL DSCAL (N, TEM, V(1,1), 1)
C Zero out the HES array. ----------------------------------------------
DO 65 J = 1,MAXL
DO 60 I = 1,MAXLP1
60 HES(I,J) = 0.0D0
65 CONTINUE
C-----------------------------------------------------------------------
C Main loop to compute the vectors V(*,2) to V(*,MAXL).
C The running product PROD is needed for the convergence test.
C-----------------------------------------------------------------------
PROD = 1.0D0
DO 90 LL = 1,MAXL
LGMR = LL
C-----------------------------------------------------------------------
C Call routine DATV to compute VNEW = Abar*v(ll), where Abar is
C the matrix A with scaling and inverse preconditioner factors applied.
C Call routine DORTHOG to orthogonalize the new vector VNEW = V(*,LL+1).
C Call routine DHEQR to update the factors of HES.
C-----------------------------------------------------------------------
CALL DATV (NEQ, Y, SAVF, V(1,LL), WGHT, X, F, PSOL, V(1,LL+1),
1 WK, WP, IWP, HL0, JPRE, IER, NPSL)
IF (IER .GT. 0) IFLAG = 3
C
RETURN
C----------------------- End of Subroutine DSPIGMR ---------------------
END
*DECK DPCG
SUBROUTINE DPCG (NEQ, TN, Y, SAVF, R, WGHT, N, MAXL, DELTA, HL0,
1 JPRE, MNEWT, F, PSOL, NPSL, X, P, W, Z, LPCG, WP, IWP, WK, IFLAG)
EXTERNAL F, PSOL
INTEGER NEQ, N, MAXL, JPRE, MNEWT, NPSL, LPCG, IWP, IFLAG
DOUBLE PRECISION TN,Y,SAVF,R,WGHT,DELTA,HL0,X,P,W,Z,WP,WK
DIMENSION NEQ(*), Y(*), SAVF(*), R(*), WGHT(*), X(*), P(*), W(*),
1 Z(*), WP(*), IWP(*), WK(*)
C-----------------------------------------------------------------------
C This routine computes the solution to the system A*x = b using a
C preconditioned version of the Conjugate Gradient algorithm.
C It is assumed here that the matrix A and the preconditioner
C matrix M are symmetric positive definite or nearly so.
C-----------------------------------------------------------------------
C
C On entry
C
C NEQ = problem size, passed to F and PSOL (NEQ(1) = N).
C
C TN = current value of t.
C
C Y = array containing current dependent variable vector.
C
C SAVF = array containing current value of f(t,y).
C
C R = the right hand side of the system A*x = b.
C
C WGHT = array of length N containing scale factors.
C 1/WGHT(i) are the diagonal elements of the diagonal
C scaling matrix D.
C
C N = the order of the matrix A, and the lengths
C of the vectors Y, SAVF, R, WGHT, P, W, Z, WK, and X.
C
C MAXL = the maximum allowable number of iterates.
C
C DELTA = tolerance on residuals b - A*x in weighted RMS-norm.
C
C HL0 = current value of (step size h) * (coefficient l0).
C
C JPRE = preconditioner type flag.
C
C MNEWT = Newton iteration counter (.ge. 0).
C
C WK = real work array used by routine DATP.
C
C WP = real work array used by preconditioner PSOL.
C
C IWP = integer work array used by preconditioner PSOL.
C
C On return
C
C X = the final computed approximation to the solution
C of the system A*x = b.
C
C LPCG = the number of iterations performed, and current
C order of the upper Hessenberg matrix HES.
C
C NPSL = the number of calls to PSOL.
C
C IFLAG = integer error flag:
C 0 means convergence in LPCG iterations, LPCG .le. MAXL.
C 1 means the convergence test did not pass in MAXL
C iterations, but the residual norm is .lt. 1,
C or .lt. norm(b) if MNEWT = 0, and so X is computed.
C 2 means the convergence test did not pass in MAXL
C iterations, residual .gt. 1, and X is undefined.
C 3 means there was a recoverable error in PSOL
C caused by the preconditioner being out of date.
C 4 means there was a zero denominator in the algorithm.
C The system matrix or preconditioner matrix is not
C sufficiently close to being symmetric pos. definite.
C -1 means there was a nonrecoverable error in PSOL.
C
C-----------------------------------------------------------------------
INTEGER I, IER
DOUBLE PRECISION ALPHA,BETA,BNRM,PTW,RNRM,DDOT,DVNORM,ZTR,ZTR0
C
IFLAG = 0
NPSL = 0
LPCG = 0
DO 10 I = 1,N
10 X(I) = 0.0D0
BNRM = DVNORM (N, R, WGHT)
C Test for immediate return with X = 0 or X = b. -----------------------
IF (BNRM .GT. DELTA) GO TO 20
IF (MNEWT .GT. 0) RETURN
CALL DCOPY (N, R, 1, X, 1)
RETURN
C
20 ZTR = 0.0D0
C Loop point for PCG iterations. ---------------------------------------
30 CONTINUE
LPCG = LPCG + 1
CALL DCOPY (N, R, 1, Z, 1)
IER = 0
IF (JPRE .EQ. 0) GO TO 40
CALL PSOL (NEQ, TN, Y, SAVF, WK, HL0, WP, IWP, Z, 3, IER)
NPSL = NPSL + 1
IF (IER .NE. 0) GO TO 100
40 CONTINUE
ZTR0 = ZTR
ZTR = DDOT (N, Z, 1, R, 1)
IF (LPCG .NE. 1) GO TO 50
CALL DCOPY (N, Z, 1, P, 1)
GO TO 70
50 CONTINUE
IF (ZTR0 .EQ. 0.0D0) GO TO 200
BETA = ZTR/ZTR0
DO 60 I = 1,N
60 P(I) = Z(I) + BETA*P(I)
70 CONTINUE
C-----------------------------------------------------------------------
C Call DATP to compute A*p and return the answer in W.
C-----------------------------------------------------------------------
CALL DATP (NEQ, Y, SAVF, P, WGHT, HL0, WK, F, W)
C
PTW = DDOT (N, P, 1, W, 1)
IF (PTW .EQ. 0.0D0) GO TO 200
ALPHA = ZTR/PTW
CALL DAXPY (N, ALPHA, P, 1, X, 1)
ALPHA = -ALPHA
CALL DAXPY (N, ALPHA, W, 1, R, 1)
RNRM = DVNORM (N, R, WGHT)
IF (RNRM .LE. DELTA) RETURN
IF (LPCG .LT. MAXL) GO TO 30
IFLAG = 2
IF (RNRM .LE. 1.0D0) IFLAG = 1
IF (RNRM .LE. BNRM .AND. MNEWT .EQ. 0) IFLAG = 1
RETURN
C-----------------------------------------------------------------------
C This block handles error returns from PSOL.
C-----------------------------------------------------------------------
100 CONTINUE
IF (IER .LT. 0) IFLAG = -1
IF (IER .GT. 0) IFLAG = 3
RETURN
C-----------------------------------------------------------------------
C This block handles division by zero errors.
C-----------------------------------------------------------------------
200 CONTINUE
IFLAG = 4
RETURN
C----------------------- End of Subroutine DPCG ------------------------
END
*DECK DPCGS
SUBROUTINE DPCGS (NEQ, TN, Y, SAVF, R, WGHT, N, MAXL, DELTA, HL0,
1 JPRE, MNEWT, F, PSOL, NPSL, X, P, W, Z, LPCG, WP, IWP, WK, IFLAG)
EXTERNAL F, PSOL
INTEGER NEQ, N, MAXL, JPRE, MNEWT, NPSL, LPCG, IWP, IFLAG
DOUBLE PRECISION TN,Y,SAVF,R,WGHT,DELTA,HL0,X,P,W,Z,WP,WK
DIMENSION NEQ(*), Y(*), SAVF(*), R(*), WGHT(*), X(*), P(*), W(*),
1 Z(*), WP(*), IWP(*), WK(*)
C-----------------------------------------------------------------------
C This routine computes the solution to the system A*x = b using a
C scaled preconditioned version of the Conjugate Gradient algorithm.
C It is assumed here that the scaled matrix D**-1 * A * D and the
C scaled preconditioner D**-1 * M * D are close to being
C symmetric positive definite.
C-----------------------------------------------------------------------
C
C On entry
C
C NEQ = problem size, passed to F and PSOL (NEQ(1) = N).
C
C TN = current value of t.
C
C Y = array containing current dependent variable vector.
C
C SAVF = array containing current value of f(t,y).
C
C R = the right hand side of the system A*x = b.
C
C WGHT = array of length N containing scale factors.
C 1/WGHT(i) are the diagonal elements of the diagonal
C scaling matrix D.
C
C N = the order of the matrix A, and the lengths
C of the vectors Y, SAVF, R, WGHT, P, W, Z, WK, and X.
C
C MAXL = the maximum allowable number of iterates.
C
C DELTA = tolerance on residuals b - A*x in weighted RMS-norm.
C
C HL0 = current value of (step size h) * (coefficient l0).
C
C JPRE = preconditioner type flag.
C
C MNEWT = Newton iteration counter (.ge. 0).
C
C WK = real work array used by routine DATP.
C
C WP = real work array used by preconditioner PSOL.
C
C IWP = integer work array used by preconditioner PSOL.
C
C On return
C
C X = the final computed approximation to the solution
C of the system A*x = b.
C
C LPCG = the number of iterations performed, and current
C order of the upper Hessenberg matrix HES.
C
C NPSL = the number of calls to PSOL.
C
C IFLAG = integer error flag:
C 0 means convergence in LPCG iterations, LPCG .le. MAXL.
C 1 means the convergence test did not pass in MAXL
C iterations, but the residual norm is .lt. 1,
C or .lt. norm(b) if MNEWT = 0, and so X is computed.
C 2 means the convergence test did not pass in MAXL
C iterations, residual .gt. 1, and X is undefined.
C 3 means there was a recoverable error in PSOL
C caused by the preconditioner being out of date.
C 4 means there was a zero denominator in the algorithm.
C the scaled matrix or scaled preconditioner is not
C sufficiently close to being symmetric pos. definite.
C -1 means there was a nonrecoverable error in PSOL.
C
C-----------------------------------------------------------------------
INTEGER I, IER
DOUBLE PRECISION ALPHA, BETA, BNRM, PTW, RNRM, DVNORM, ZTR, ZTR0
C
IFLAG = 0
NPSL = 0
LPCG = 0
DO 10 I = 1,N
10 X(I) = 0.0D0
BNRM = DVNORM (N, R, WGHT)
C Test for immediate return with X = 0 or X = b. -----------------------
IF (BNRM .GT. DELTA) GO TO 20
IF (MNEWT .GT. 0) RETURN
CALL DCOPY (N, R, 1, X, 1)
RETURN
C
20 ZTR = 0.0D0
C Loop point for PCG iterations. ---------------------------------------
30 CONTINUE
LPCG = LPCG + 1
CALL DCOPY (N, R, 1, Z, 1)
IER = 0
IF (JPRE .EQ. 0) GO TO 40
CALL PSOL (NEQ, TN, Y, SAVF, WK, HL0, WP, IWP, Z, 3, IER)
NPSL = NPSL + 1
IF (IER .NE. 0) GO TO 100
40 CONTINUE
ZTR0 = ZTR
ZTR = 0.0D0
DO 45 I = 1,N
45 ZTR = ZTR + Z(I)*R(I)*WGHT(I)**2
IF (LPCG .NE. 1) GO TO 50
CALL DCOPY (N, Z, 1, P, 1)
GO TO 70
50 CONTINUE
IF (ZTR0 .EQ. 0.0D0) GO TO 200
BETA = ZTR/ZTR0
DO 60 I = 1,N
60 P(I) = Z(I) + BETA*P(I)
70 CONTINUE
C-----------------------------------------------------------------------
C Call DATP to compute A*p and return the answer in W.
C-----------------------------------------------------------------------
CALL DATP (NEQ, Y, SAVF, P, WGHT, HL0, WK, F, W)
C
PTW = 0.0D0
DO 80 I = 1,N
80 PTW = PTW + P(I)*W(I)*WGHT(I)**2
IF (PTW .EQ. 0.0D0) GO TO 200
ALPHA = ZTR/PTW
CALL DAXPY (N, ALPHA, P, 1, X, 1)
ALPHA = -ALPHA
CALL DAXPY (N, ALPHA, W, 1, R, 1)
RNRM = DVNORM (N, R, WGHT)
IF (RNRM .LE. DELTA) RETURN
IF (LPCG .LT. MAXL) GO TO 30
IFLAG = 2
IF (RNRM .LE. 1.0D0) IFLAG = 1
IF (RNRM .LE. BNRM .AND. MNEWT .EQ. 0) IFLAG = 1
RETURN
C-----------------------------------------------------------------------
C This block handles error returns from PSOL.
C-----------------------------------------------------------------------
100 CONTINUE
IF (IER .LT. 0) IFLAG = -1
IF (IER .GT. 0) IFLAG = 3
RETURN
C-----------------------------------------------------------------------
C This block handles division by zero errors.
C-----------------------------------------------------------------------
200 CONTINUE
IFLAG = 4
RETURN
C----------------------- End of Subroutine DPCGS -----------------------
END
*DECK DATP
SUBROUTINE DATP (NEQ, Y, SAVF, P, WGHT, HL0, WK, F, W)
EXTERNAL F
IQ = 2*(K-1) + 1
C = Q(IQ)
S = Q(IQ+1)
T1 = B(K)
T2 = B(KP1)
B(K) = C*T1 - S*T2
B(KP1) = S*T1 + C*T2
20 CONTINUE
C
C Now solve R*x = Q*b.
C
DO 40 KB = 1, N
K = N + 1 - KB
B(K) = B(K)/A(K,K)
T = -B(K)
CALL DAXPY (K-1, T, A(1,K), 1, B(1), 1)
40 CONTINUE
RETURN
C----------------------- End of Subroutine DHELS -----------------------
END
*DECK DLHIN
SUBROUTINE DLHIN (NEQ, N, T0, Y0, YDOT, F, TOUT, UROUND,
1 EWT, ITOL, ATOL, Y, TEMP, H0, NITER, IER)
EXTERNAL F
DOUBLE PRECISION T0, Y0, YDOT, TOUT, UROUND, EWT, ATOL, Y,
1 TEMP, H0
INTEGER NEQ, N, ITOL, NITER, IER
DIMENSION NEQ(*), Y0(*), YDOT(*), EWT(*), ATOL(*), Y(*), TEMP(*)
C-----------------------------------------------------------------------
C Call sequence input -- NEQ, N, T0, Y0, YDOT, F, TOUT, UROUND,
C EWT, ITOL, ATOL, Y, TEMP
C Call sequence output -- H0, NITER, IER
C Common block variables accessed -- None
C
C Subroutines called by DLHIN: F, DCOPY
C Function routines called by DLHIN: DVNORM
C-----------------------------------------------------------------------
C This routine computes the step size, H0, to be attempted on the
C first step, when the user has not supplied a value for this.
C
C First we check that TOUT - T0 differs significantly from zero. Then
C an iteration is done to approximate the initial second derivative
C and this is used to define H from WRMS-norm(H**2 * yddot / 2) = 1.
C A bias factor of 1/2 is applied to the resulting h.
C The sign of H0 is inferred from the initial values of TOUT and T0.
C
C Communication with DLHIN is done with the following variables:
C
C NEQ = NEQ array of solver, passed to F.
C N = size of ODE system, input.
C T0 = initial value of independent variable, input.
C Y0 = vector of initial conditions, input.
C YDOT = vector of initial first derivatives, input.
C F = name of subroutine for right-hand side f(t,y), input.
C TOUT = first output value of independent variable
C UROUND = machine unit roundoff
C EWT, ITOL, ATOL = error weights and tolerance parameters
C as described in the driver routine, input.
C Y, TEMP = work arrays of length N.
C H0 = step size to be attempted, output.
C NITER = number of iterations (and of f evaluations) to compute H0,
C output.
C IER = the error flag, returned with the value
C IER = 0 if no trouble occurred, or
C IER = -1 if TOUT and t0 are considered too close to proceed.
C-----------------------------------------------------------------------
C
C Type declarations for local variables --------------------------------
C
DOUBLE PRECISION AFI, ATOLI, DELYI, HALF, HG, HLB, HNEW, HRAT,
1 HUB, HUN, PT1, T1, TDIST, TROUND, TWO, DVNORM, YDDNRM
INTEGER I, ITER
C-----------------------------------------------------------------------
C The following Fortran-77 declaration is to cause the values of the
C listed (local) variables to be saved between calls to this integrator.
C-----------------------------------------------------------------------
SAVE HALF, HUN, PT1, TWO
DATA HALF /0.5D0/, HUN /100.0D0/, PT1 /0.1D0/, TWO /2.0D0/
C
NITER = 0
TDIST = ABS(TOUT - T0)
TROUND = UROUND*MAX(ABS(T0),ABS(TOUT))
IF (TDIST .LT. TWO*TROUND) GO TO 100
C
C Set a lower bound on H based on the roundoff level in T0 and TOUT. ---
HLB = HUN*TROUND
C Set an upper bound on H based on TOUT-T0 and the initial Y and YDOT. -
HUB = PT1*TDIST
ATOLI = ATOL(1)
DO 10 I = 1,N
IF (ITOL .EQ. 2 .OR. ITOL .EQ. 4) ATOLI = ATOL(I)
DELYI = PT1*ABS(Y0(I)) + ATOLI
AFI = ABS(YDOT(I))
IF (AFI*HUB .GT. DELYI) HUB = DELYI/AFI
10 CONTINUE
C
C Set initial guess for H as geometric mean of upper and lower bounds. -
ITER = 0
HG = SQRT(HLB*HUB)
C If the bounds have crossed, exit with the mean value. ----------------
IF (HUB .LT. HLB) THEN
H0 = HG
GO TO 90
ENDIF
C
C Looping point for iteration. -----------------------------------------
50 CONTINUE
C Estimate the second derivative as a difference quotient in f. --------
T1 = T0 + HG
DO 60 I = 1,N
60 Y(I) = Y0(I) + HG*YDOT(I)
CALL F (NEQ, T1, Y, TEMP)
DO 70 I = 1,N
70 TEMP(I) = (TEMP(I) - YDOT(I))/HG
YDDNRM = DVNORM (N, TEMP, EWT)
C Get the corresponding new value of H. --------------------------------
IF (YDDNRM*HUB*HUB .GT. TWO) THEN
HNEW = SQRT(TWO/YDDNRM)
ELSE
HNEW = SQRT(HG*HUB)
ENDIF
ITER = ITER + 1
C-----------------------------------------------------------------------
C Test the stopping conditions.
C Stop if the new and previous H values differ by a factor of .lt. 2.
C Stop if four iterations have been done. Also, stop with previous H
C if hnew/hg .gt. 2 after first iteration, as this probably means that
C the second derivative value is bad because of cancellation error.
C-----------------------------------------------------------------------
IF (ITER .GE. 4) GO TO 80
HRAT = HNEW/HG
IF ( (HRAT .GT. HALF) .AND. (HRAT .LT. TWO) ) GO TO 80
IF ( (ITER .GE. 2) .AND. (HNEW .GT. TWO*HG) ) THEN
HNEW = HG
GO TO 80
ENDIF
HG = HNEW
GO TO 50
C
C Iteration done. Apply bounds, bias factor, and sign. ----------------
80 H0 = HNEW*HALF
IF (H0 .LT. HLB) H0 = HLB
IF (H0 .GT. HUB) H0 = HUB
90 H0 = SIGN(H0, TOUT - T0)
C Restore Y array from Y0, then exit. ----------------------------------
CALL DCOPY (N, Y0, 1, Y, 1)
NITER = ITER
IER = 0
RETURN
C Error return for TOUT - T0 too small. --------------------------------
100 IER = -1
RETURN
C----------------------- End of Subroutine DLHIN -----------------------
END
*DECK DSTOKA
SUBROUTINE DSTOKA (NEQ, Y, YH, NYH, YH1, EWT, SAVF, SAVX, ACOR,
1 WM, IWM, F, JAC, PSOL)
EXTERNAL F, JAC, PSOL
INTEGER NEQ, NYH, IWM
DOUBLE PRECISION Y, YH, YH1, EWT, SAVF, SAVX, ACOR, WM
DIMENSION NEQ(*), Y(*), YH(NYH,*), YH1(*), EWT(*), SAVF(*),
1 SAVX(*), ACOR(*), WM(*), IWM(*)
INTEGER IOWND, IALTH, IPUP, LMAX, MEO, NQNYH, NSLP,
1 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L,
2 LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER,
3 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
INTEGER NEWT, NSFI, NSLJ, NJEV
INTEGER JPRE, JACFLG, LOCWP, LOCIWP, LSAVX, KMP, MAXL, MNEWT,
1 NNI, NLI, NPS, NCFN, NCFL
DOUBLE PRECISION CONIT, CRATE, EL, ELCO, HOLD, RMAX, TESCO,
2 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
DOUBLE PRECISION STIFR
DOUBLE PRECISION DELT, EPCON, SQRTN, RSQRTN
COMMON /DLS001/ CONIT, CRATE, EL(13), ELCO(13,12),
1 HOLD, RMAX, TESCO(3,12),
2 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND,
3 IOWND(6), IALTH, IPUP, LMAX, MEO, NQNYH, NSLP,
3 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L,
4 LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER,
5 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
COMMON /DLS002/ STIFR, NEWT, NSFI, NSLJ, NJEV
COMMON /DLPK01/ DELT, EPCON, SQRTN, RSQRTN,
1 JPRE, JACFLG, LOCWP, LOCIWP, LSAVX, KMP, MAXL, MNEWT,
2 NNI, NLI, NPS, NCFN, NCFL
C-----------------------------------------------------------------------
C DSTOKA performs one step of the integration of an initial value
C Note: DSTOKA is independent of the value of the iteration method
C indicator MITER, when this is .ne. 0, and hence is independent
C of the type of chord method used, or the Jacobian structure.
C Communication with DSTOKA is done with the following variables:
C
C NEQ = integer array containing problem size in NEQ(1), and
C passed as the NEQ argument in all calls to F and JAC.
C Y = an array of length .ge. N used as the Y argument in
C all calls to F and JAC.
C YH = an NYH by LMAX array containing the dependent variables
C and their approximate scaled derivatives, where
C LMAX = MAXORD + 1. YH(i,j+1) contains the approximate
C j-th derivative of y(i), scaled by H**j/factorial(j)
C (j = 0,1,...,NQ). On entry for the first step, the first
C two columns of YH must be set from the initial values.
C NYH = a constant integer .ge. N, the first dimension of YH.
C YH1 = a one-dimensional array occupying the same space as YH.
C EWT = an array of length N containing multiplicative weights
C for local error measurements. Local errors in y(i) are
C compared to 1.0/EWT(i) in various error tests.
C SAVF = an array of working storage, of length N.
C Also used for input of YH(*,MAXORD+2) when JSTART = -1
C and MAXORD .lt. the current order NQ.
C SAVX = an array of working storage, of length N.
C ACOR = a work array of length N, used for the accumulated
C corrections. On a successful return, ACOR(i) contains
C the estimated one-step local error in y(i).
C WM,IWM = real and integer work arrays associated with matrix
C operations in chord iteration (MITER .ne. 0).
C CCMAX = maximum relative change in H*EL0 before DSETPK is called.
C H = the step size to be attempted on the next step.
C H is altered by the error control algorithm during the
C problem. H can be either positive or negative, but its
C sign must remain constant throughout the problem.
C HMIN = the minimum absolute value of the step size H to be used.
C HMXI = inverse of the maximum absolute value of H to be used.
C HMXI = 0.0 is allowed and corresponds to an infinite HMAX.
C HMIN and HMXI may be changed at any time, but will not
C take effect until the next change of H is considered.
C TN = the independent variable. TN is updated on each step taken.
C JSTART = an integer used for input only, with the following
C values and meanings:
C 0 perform the first step.
C .gt.0 take a new step continuing from the last.
C -1 take the next step with a new value of H, MAXORD,
C N, METH, MITER, and/or matrix parameters.
C -2 take the next step with a new value of H,
C but with other inputs unchanged.
C On return, JSTART is set to 1 to facilitate continuation.
C KFLAG = a completion code with the following meanings:
C 0 the step was succesful.
C -1 the requested error could not be achieved.
C -2 corrector convergence could not be achieved.
C -3 fatal error in DSETPK or DSOLPK.
C A return with KFLAG = -1 or -2 means either
C ABS(H) = HMIN or 10 consecutive failures occurred.
C On a return with KFLAG negative, the values of TN and
C the YH array are as of the beginning of the last
C step, and H is the last step size attempted.
C MAXORD = the maximum order of integration method to be allowed.
C MAXCOR = the maximum number of corrector iterations allowed.
C MSBP = maximum number of steps between DSETPK calls (MITER .gt. 0).
C MXNCF = maximum number of convergence failures allowed.
C METH/MITER = the method flags. See description in driver.
C N = the number of first-order differential equations.
C-----------------------------------------------------------------------
INTEGER I, I1, IREDO, IRET, J, JB, JOK, M, NCF, NEWQ, NSLOW
DOUBLE PRECISION DCON, DDN, DEL, DELP, DRC, DSM, DUP, EXDN, EXSM,
1 EXUP, DFNORM, R, RH, RHDN, RHSM, RHUP, ROC, STIFF, TOLD, DVNORM
C
KFLAG = 0
TOLD = TN
NCF = 0
IERPJ = 0
IERSL = 0
JCUR = 0
ICF = 0
DELP = 0.0D0
IF (JSTART .GT. 0) GO TO 200
IF (JSTART .EQ. -1) GO TO 100
IF (JSTART .EQ. -2) GO TO 160
C-----------------------------------------------------------------------
C On the first call, the order is set to 1, and other variables are
C initialized. RMAX is the maximum ratio by which H can be increased
C in a single step. It is initially 1.E4 to compensate for the small
C initial H, but then is normally equal to 10. If a failure
C occurs (in corrector convergence or error test), RMAX is set at 2
C for the next increase.
C-----------------------------------------------------------------------
LMAX = MAXORD + 1
NQ = 1
L = 2
IALTH = 2
RMAX = 10000.0D0
RC = 0.0D0
EL0 = 1.0D0
CRATE = 0.7D0
HOLD = H
MEO = METH
NSLP = 0
NSLJ = 0
IPUP = 0
IRET = 3
NEWT = 0
STIFR = 0.0D0
GO TO 140
C-----------------------------------------------------------------------
C The following block handles preliminaries needed when JSTART = -1.
C IPUP is set to MITER to force a matrix update.
C If an order increase is about to be considered (IALTH = 1),
C IALTH is reset to 2 to postpone consideration one more step.
C If the caller has changed METH, DCFODE is called to reset
C the coefficients of the method.
C If the caller has changed MAXORD to a value less than the current
C order NQ, NQ is reduced to MAXORD, and a new H chosen accordingly.
C If H is to be changed, YH must be rescaled.
C If H or METH is being changed, IALTH is reset to L = NQ + 1
C to prevent further changes in H for that many steps.
C-----------------------------------------------------------------------
100 IPUP = MITER
LMAX = MAXORD + 1
C-----------------------------------------------------------------------
140 CALL DCFODE (METH, ELCO, TESCO)
150 DO 155 I = 1,L
155 EL(I) = ELCO(I,NQ)
NQNYH = NQ*NYH
RC = RC*EL(1)/EL0
EL0 = EL(1)
CONIT = 0.5D0/(NQ+2)
EPCON = CONIT*TESCO(2,NQ)
GO TO (160, 170, 200), IRET
C-----------------------------------------------------------------------
C If H is being changed, the H ratio RH is checked against
C RMAX, HMIN, and HMXI, and the YH array rescaled. IALTH is set to
C L = NQ + 1 to prevent a change of H for that many steps, unless
C forced by a convergence or error test failure.
C-----------------------------------------------------------------------
160 IF (H .EQ. HOLD) GO TO 200
RH = H/HOLD
H = HOLD
IREDO = 3
GO TO 175
170 RH = MAX(RH,HMIN/ABS(H))
175 RH = MIN(RH,RMAX)
RH = RH/MAX(1.0D0,ABS(H)*HMXI*RH)
R = 1.0D0
DO 180 J = 2,L
R = R*RH
DO 180 I = 1,N
180 YH(I,J) = YH(I,J)*R
H = H*RH
RC = RC*RH
IALTH = L
IF (IREDO .EQ. 0) GO TO 690
C-----------------------------------------------------------------------
C This section computes the predicted values by effectively
C multiplying the YH array by the Pascal triangle matrix.
C The flag IPUP is set according to whether matrix data is involved
C (NEWT .gt. 0 .and. JACFLG .ne. 0) or not, to trigger a call to DSETPK.
C IPUP is set to MITER when RC differs from 1 by more than CCMAX,
C and at least every MSBP steps, when JACFLG = 1.
C RC is the ratio of new to old values of the coefficient H*EL(1).
C-----------------------------------------------------------------------
200 IF (NEWT .EQ. 0 .OR. JACFLG .EQ. 0) THEN
DRC = 0.0D0
IPUP = 0
CRATE = 0.7D0
ELSE
DRC = ABS(RC - 1.0D0)
IF (DRC .GT. CCMAX) IPUP = MITER
IF (NST .GE. NSLP+MSBP) IPUP = MITER
ENDIF
TN = TN + H
I1 = NQNYH + 1
DO 215 JB = 1,NQ
I1 = I1 - NYH
CDIR$ IVDEP
DO 210 I = I1,NQNYH
210 YH1(I) = YH1(I) + YH1(I+NYH)
215 CONTINUE
C-----------------------------------------------------------------------
C Up to MAXCOR corrector iterations are taken. A convergence test is
C made on the RMS-norm of each correction, weighted by the error
C weight vector EWT. The sum of the corrections is accumulated in the
C vector ACOR(i). The YH array is not altered in the corrector loop.
C Within the corrector loop, an estimated rate of convergence (ROC)
C and a stiffness ratio estimate (STIFF) are kept. Corresponding
C global estimates are kept as CRATE and stifr.
C-----------------------------------------------------------------------
220 M = 0
MNEWT = 0
STIFF = 0.0D0
ROC = 0.05D0
NSLOW = 0
DO 230 I = 1,N
230 Y(I) = YH(I,1)
CALL F (NEQ, TN, Y, SAVF)
NFE = NFE + 1
IF (NEWT .EQ. 0 .OR. IPUP .LE. 0) GO TO 250
C-----------------------------------------------------------------------
C If indicated, DSETPK is called to update any matrix data needed,
C before starting the corrector iteration.
C JOK is set to indicate if the matrix data need not be recomputed.
C IPUP is set to 0 as an indicator that the matrix data is up to date.
C-----------------------------------------------------------------------
JOK = 1
IF (NST .EQ. 0 .OR. NST .GT. NSLJ+50) JOK = -1
IF (ICF .EQ. 1 .AND. DRC .LT. 0.2D0) JOK = -1
IF (ICF .EQ. 2) JOK = -1
IF (JOK .EQ. -1) THEN
NSLJ = NST
NJEV = NJEV + 1
ENDIF
CALL DSETPK (NEQ, Y, YH1, EWT, ACOR, SAVF, JOK, WM, IWM, F, JAC)
IPUP = 0
RC = 1.0D0
DRC = 0.0D0
NSLP = NST
CRATE = 0.7D0
IF (IERPJ .NE. 0) GO TO 430
250 DO 260 I = 1,N
260 ACOR(I) = 0.0D0
270 IF (NEWT .NE. 0) GO TO 350
C-----------------------------------------------------------------------
C In the case of functional iteration, update Y directly from
C the result of the last function evaluation, and STIFF is set to 1.0.
C-----------------------------------------------------------------------
DO 290 I = 1,N
SAVF(I) = H*SAVF(I) - YH(I,2)
290 Y(I) = SAVF(I) - ACOR(I)
DEL = DVNORM (N, Y, EWT)
DO 300 I = 1,N
Y(I) = YH(I,1) + EL(1)*SAVF(I)
300 ACOR(I) = SAVF(I)
STIFF = 1.0D0
GO TO 400
C-----------------------------------------------------------------------
C In the case of the chord method, compute the corrector error,
C and solve the linear system with that as right-hand side and
C P as coefficient matrix. STIFF is set to the ratio of the norms
C of the residual and the correction vector.
C-----------------------------------------------------------------------
C all calls to RES, JAC, and ADDA.
C NEQ = integer array containing problem size in NEQ(1), and
C passed as the NEQ argument in all calls tO RES, G, ADDA,
C and JAC.
C YH = an NYH by LMAX array containing the dependent variables
C and their approximate scaled derivatives, where
C LMAX = MAXORD + 1. YH(i,j+1) contains the approximate
C j-th derivative of y(i), scaled by H**j/factorial(j)
C (j = 0,1,...,NQ). On entry for the first step, the first
C two columns of YH must be set from the initial values.
C NYH = a constant integer .ge. N, the first dimension of YH.
C YH1 = a one-dimensional array occupying the same space as YH.
C EWT = an array of length N containing multiplicative weights
C for local error measurements. Local errors in y(i) are
C compared to 1.0/EWT(i) in various error tests.
C SAVF = an array of working storage, of length N. also used for
C input of YH(*,MAXORD+2) when JSTART = -1 and MAXORD is less
C than the current order NQ.
C Same as YDOTI in the driver.
C SAVR = an array of working storage, of length N.
C ACOR = a work array of length N used for the accumulated
C corrections. On a succesful return, ACOR(i) contains
C the estimated one-step local error in y(i).
C WM,IWM = real and integer work arrays associated with matrix
C operations in chord iteration.
C PJAC = name of routine to evaluate and preprocess Jacobian matrix.
C SLVS = name of routine to solve linear system in chord iteration.
C CCMAX = maximum relative change in H*EL0 before PJAC is called.
C H = the step size to be attempted on the next step.
C H is altered by the error control algorithm during the
C problem. H can be either positive or negative, but its
C sign must remain constant throughout the problem.
C HMIN = the minimum absolute value of the step size H to be used.
C HMXI = inverse of the maximum absolute value of H to be used.
C HMXI = 0.0 is allowed and corresponds to an infinite HMAX.
C HMIN and HMXI may be changed at any time, but will not
C take effect until the next change of H is considered.
C TN = the independent variable. TN is updated on each step taken.
C JSTART = an integer used for input only, with the following
C values and meanings:
C 0 perform the first step.
C .gt.0 take a new step continuing from the last.
C -1 take the next step with a new value of H, MAXORD,
C N, METH, MITER, and/or matrix parameters.
C -2 take the next step with a new value of H,
C but with other inputs unchanged.
C On return, JSTART is set to 1 to facilitate continuation.
C KFLAG = a completion code with the following meanings:
C 0 the step was succesful.
C -1 the requested error could not be achieved.
C -2 corrector convergence could not be achieved.
C -3 RES ordered immediate return.
C -4 error condition from RES could not be avoided.
C -5 fatal error in PJAC or SLVS.
C A return with KFLAG = -1, -2, or -4 means either
C ABS(H) = HMIN or 10 consecutive failures occurred.
C On a return with KFLAG negative, the values of TN and
C the YH array are as of the beginning of the last
C step, and H is the last step size attempted.
C MAXORD = the maximum order of integration method to be allowed.
C MAXCOR = the maximum number of corrector iterations allowed.
C MSBP = maximum number of steps between PJAC calls.
C MXNCF = maximum number of convergence failures allowed.
C METH/MITER = the method flags. See description in driver.
C N = the number of first-order differential equations.
C-----------------------------------------------------------------------
KFLAG = 0
TOLD = TN
NCF = 0
IERPJ = 0
IERSL = 0
JCUR = 0
ICF = 0
DELP = 0.0D0
IF (JSTART .GT. 0) GO TO 200
IF (JSTART .EQ. -1) GO TO 100
IF (JSTART .EQ. -2) GO TO 160
C-----------------------------------------------------------------------
C On the first call, the order is set to 1, and other variables are
C initialized. RMAX is the maximum ratio by which H can be increased
C in a single step. It is initially 1.E4 to compensate for the small
C initial H, but then is normally equal to 10. If a failure
C occurs (in corrector convergence or error test), RMAX is set at 2
C for the next increase.
C-----------------------------------------------------------------------
LMAX = MAXORD + 1
NQ = 1
L = 2
IALTH = 2
RMAX = 10000.0D0
RC = 0.0D0
EL0 = 1.0D0
CRATE = 0.7D0
HOLD = H
MEO = METH
NSLP = 0
IPUP = MITER
IRET = 3
GO TO 140
C-----------------------------------------------------------------------
C The following block handles preliminaries needed when JSTART = -1.
C IPUP is set to MITER to force a matrix update.
C If an order increase is about to be considered (IALTH = 1),
C IALTH is reset to 2 to postpone consideration one more step.
C If the caller has changed METH, DCFODE is called to reset
C the coefficients of the method.
C If the caller has changed MAXORD to a value less than the current
C order NQ, NQ is reduced to MAXORD, and a new H chosen accordingly.
C If H is to be changed, YH must be rescaled.
C If H or METH is being changed, IALTH is reset to L = NQ + 1
C to prevent further changes in H for that many steps.
C-----------------------------------------------------------------------
100 IPUP = MITER
LMAX = MAXORD + 1
IF (IALTH .EQ. 1) IALTH = 2
IF (METH .EQ. MEO) GO TO 110
CALL DCFODE (METH, ELCO, TESCO)
MEO = METH
IF (NQ .GT. MAXORD) GO TO 120
IALTH = L
IRET = 1
IREDO = 3
IF (H .EQ. HOLD) GO TO 170
RH = MIN(RH,ABS(H/HOLD))
H = HOLD
GO TO 175
C-----------------------------------------------------------------------
C DCFODE is called to get all the integration coefficients for the
C current METH. Then the EL vector and related constants are reset
C whenever the order NQ is changed, or at the start of the problem.
C-----------------------------------------------------------------------
140 CALL DCFODE (METH, ELCO, TESCO)
150 DO 155 I = 1,L
155 EL(I) = ELCO(I,NQ)
NQNYH = NQ*NYH
RC = RC*EL(1)/EL0
EL0 = EL(1)
CONIT = 0.5D0/(NQ+2)
GO TO (160, 170, 200), IRET
C-----------------------------------------------------------------------
C If H is being changed, the H ratio RH is checked against
C RMAX, HMIN, and HMXI, and the YH array rescaled. IALTH is set to
C L = NQ + 1 to prevent a change of H for that many steps, unless
C forced by a convergence or error test failure.
C-----------------------------------------------------------------------
160 IF (H .EQ. HOLD) GO TO 200
RH = H/HOLD
H = HOLD
IREDO = 3
GO TO 175
170 RH = MAX(RH,HMIN/ABS(H))
175 RH = MIN(RH,RMAX)
RH = RH/MAX(1.0D0,ABS(H)*HMXI*RH)
R = 1.0D0
DO 180 J = 2,L
R = R*RH
DO 180 I = 1,N
180 YH(I,J) = YH(I,J)*R
H = H*RH
RC = RC*RH
IALTH = L
IF (IREDO .EQ. 0) GO TO 690
C-----------------------------------------------------------------------
C This section computes the predicted values by effectively
C multiplying the YH array by the Pascal triangle matrix.
C RC is the ratio of new to old values of the coefficient H*EL(1).
C When RC differs from 1 by more than CCMAX, IPUP is set to MITER
C to force PJAC to be called.
C In any case, PJAC is called at least every MSBP steps.
C-----------------------------------------------------------------------
200 IF (ABS(RC-1.0D0) .GT. CCMAX) IPUP = MITER
IF (NST .GE. NSLP+MSBP) IPUP = MITER
TN = TN + H
I1 = NQNYH + 1
DO 215 JB = 1,NQ
I1 = I1 - NYH
CDIR$ IVDEP
DO 210 I = I1,NQNYH
210 YH1(I) = YH1(I) + YH1(I+NYH)
215 CONTINUE
C-----------------------------------------------------------------------
C Up to MAXCOR corrector iterations are taken. A convergence test is
C made on the RMS-norm of each correction, weighted by H and the
C error weight vector EWT. The sum of the corrections is accumulated
C in ACOR(i). The YH array is not altered in the corrector loop.
C-----------------------------------------------------------------------
220 M = 0
DO 230 I = 1,N
SAVF(I) = YH(I,2) / H
230 Y(I) = YH(I,1)
IF (IPUP .LE. 0) GO TO 240
C-----------------------------------------------------------------------
C If indicated, the matrix P = A - H*EL(1)*dr/dy is reevaluated and
C preprocessed before starting the corrector iteration. IPUP is set
C to 0 as an indicator that this has been done.
C-----------------------------------------------------------------------
CALL PJAC (NEQ, Y, YH, NYH, EWT, ACOR, SAVR, SAVF, WM, IWM,
1 RES, JAC, ADDA )
IPUP = 0
RC = 1.0D0
NSLP = NST
CRATE = 0.7D0
IF (IERPJ .EQ. 0) GO TO 250
IF (IERPJ .LT. 0) GO TO 435
IRES = IERPJ
GO TO (430, 435, 430), IRES
C Get residual at predicted values, if not already done in PJAC. -------
240 IRES = 1
CALL RES ( NEQ, TN, Y, SAVF, SAVR, IRES )
NFE = NFE + 1
KGO = ABS(IRES)
GO TO ( 250, 435, 430 ) , KGO
250 DO 260 I = 1,N
260 ACOR(I) = 0.0D0
C-----------------------------------------------------------------------
C Solve the linear system with the current residual as
C right-hand side and P as coefficient matrix.
C-----------------------------------------------------------------------
270 CONTINUE
CALL SLVS (WM, IWM, SAVR, SAVF)
IF (IERSL .LT. 0) GO TO 430
IF (IERSL .GT. 0) GO TO 410
EL1H = EL(1) * H
DEL = DVNORM (N, SAVR, EWT) * ABS(H)
DO 380 I = 1,N
ACOR(I) = ACOR(I) + SAVR(I)
SAVF(I) = ACOR(I) + YH(I,2)/H
380 Y(I) = YH(I,1) + EL1H*ACOR(I)
C-----------------------------------------------------------------------
C Test for convergence. If M .gt. 0, an estimate of the convergence
C rate constant is stored in CRATE, and this is used in the test.
C-----------------------------------------------------------------------
IF (M .NE. 0) CRATE = MAX(0.2D0*CRATE,DEL/DELP)
DCON = DEL*MIN(1.0D0,1.5D0*CRATE)/(TESCO(2,NQ)*CONIT)
IF (DCON .LE. 1.0D0) GO TO 460
M = M + 1
IF (M .EQ. MAXCOR) GO TO 410
IF (M .GE. 2 .AND. DEL .GT. 2.0D0*DELP) GO TO 410
DELP = DEL
IRES = 1
CALL RES ( NEQ, TN, Y, SAVF, SAVR, IRES )
NFE = NFE + 1
( run in 0.539 second using v1.01-cache-2.11-cpan-71847e10f99 )