SUBROUTINE SDRIV2 (N,T,Y,F,TOUT,MSTATE,NROOT,EPS,EWT,MINT,WORK,
8 LENW,IWORK,LENIW,G)
C***BEGIN PROLOGUE SDRIV2
C***DATE WRITTEN 790601 (YYMMDD)
C***REVISION DATE 871105 (YYMMDD)
C***CATEGORY NO. I1A2,I1A1B
C***KEYWORDS ODE,STIFF,ORDINARY DIFFERENTIAL EQUATIONS,
C INITIAL VALUE PROBLEMS,GEAR'S METHOD,
C SINGLE PRECISION
C***AUTHOR KAHANER, D. K., NATIONAL BUREAU OF STANDARDS,
C SUTHERLAND, C. D., LOS ALAMOS NATIONAL LABORATORY
C***PURPOSE The function of SDRIV2 is to solve N ordinary differential
C equations of the form dY(I)/dT = F(Y(I),T), given the
C initial conditions Y(I) = YI. The program has options to
C allow the solution of both stiff and non-stiff differential
C equations. SDRIV2 uses single precision arithmetic.
C***DESCRIPTION
C From the book "Numerical Methods and Software"
C by D. Kahaner, C. Moler, S. Nash
C Prentice Hall 1988
C
C I. ABSTRACT .......................................................
C
C The function of SDRIV2 is to solve N ordinary differential
C equations of the form dY(I)/dT = F(Y(I),T), given the initial
C conditions Y(I) = YI. The program has options to allow the
C solution of both stiff and non-stiff differential equations.
C SDRIV2 is to be called once for each output point of T.
C
C II. PARAMETERS ....................................................
C
C The user should use parameter names in the call sequence of SDRIV2
C for those quantities whose value may be altered by SDRIV2. The
C parameters in the call sequence are:
C
C N = (Input) The number of differential equations.
C
C T = The independent variable. On input for the first call, T
C is the initial point. On output, T is the point at which
C the solution is given.
C
C Y = The vector of dependent variables. Y is used as input on
C the first call, to set the initial values. On output, Y
C is the computed solution vector. This array Y is passed
C in the call sequence of the user-provided routines F and
C G. Thus parameters required by F and G can be stored in
C this array in components N+1 and above. (Note: Changes
C by the user to the first N components of this array will
C take effect only after a restart, i.e., after setting
C MSTATE to +1(-1).)
C
C F = A subroutine supplied by the user. The name must be
C declared EXTERNAL in the user's calling program. This
C subroutine is of the form:
C SUBROUTINE F (N, T, Y, YDOT)
C REAL Y(*), YDOT(*)
C .
C .
C YDOT(1) = ...
C .
C .
C YDOT(N) = ...
C END (Sample)
C This computes YDOT = F(Y,T), the right hand side of the
C differential equations. Here Y is a vector of length at
C least N. The actual length of Y is determined by the
C user's declaration in the program which calls SDRIV2.
C Thus the dimensioning of Y in F, while required by FORTRAN
C convention, does not actually allocate any storage. When
C this subroutine is called, the first N components of Y are
C intermediate approximations to the solution components.
C The user should not alter these values. Here YDOT is a
C vector of length N. The user should only compute YDOT(I)
C for I from 1 to N. Normally a return from F passes
C control back to SDRIV2. However, if the user would like
C to abort the calculation, i.e., return control to the
C program which calls SDRIV2, he should set N to zero.
C SDRIV2 will signal this by returning a value of MSTATE
C equal to +6(-6). Altering the value of N in F has no
C effect on the value of N in the call sequence of SDRIV2.
C
C TOUT = (Input) The point at which the solution is desired.
C
C MSTATE = An integer describing the status of integration. The user
C must initialize MSTATE to +1 or -1. If MSTATE is
C positive, the routine will integrate past TOUT and
C interpolate the solution. This is the most efficient
C mode. If MSTATE is negative, the routine will adjust its
C internal step to reach TOUT exactly (useful if a
C singularity exists beyond TOUT.) The meaning of the
C magnitude of MSTATE:
C 1 (Input) Means the first call to the routine. This
C value must be set by the user. On all subsequent
C calls the value of MSTATE should be tested by the
C user. Unless SDRIV2 is to be reinitialized, only the
C sign of MSTATE may be changed by the user. (As a
C convenience to the user who may wish to put out the
C initial conditions, SDRIV2 can be called with
C MSTATE=+1(-1), and TOUT=T. In this case the program
C will return with MSTATE unchanged, i.e.,
C MSTATE=+1(-1).)
C 2 (Output) Means a successful integration. If a normal
C continuation is desired (i.e., a further integration
C in the same direction), simply advance TOUT and call
C again. All other parameters are automatically set.
C 3 (Output)(Unsuccessful) Means the integrator has taken
C 1000 steps without reaching TOUT. The user can
C continue the integration by simply calling SDRIV2
C again. Other than an error in problem setup, the
C most likely cause for this condition is trying to
C integrate a stiff set of equations with the non-stiff
C integrator option. (See description of MINT below.)
C 4 (Output)(Unsuccessful) Means too much accuracy has
C been requested. EPS has been increased to a value
C the program estimates is appropriate. The user can
C continue the integration by simply calling SDRIV2
C again.
C 5 (Output) A root was found at a point less than TOUT.
C The user can continue the integration toward TOUT by
C simply calling SDRIV2 again.
C 6 (Output)(Unsuccessful) N has been set to zero in
C SUBROUTINE F.
C 7 (Output)(Unsuccessful) N has been set to zero in
C FUNCTION G. See description of G below.
C
C NROOT = (Input) The number of equations whose roots are desired.
C If NROOT is zero, the root search is not active. This
C option is useful for obtaining output at points which are
C not known in advance, but depend upon the solution, e.g.,
C when some solution component takes on a specified value.
C The root search is carried out using the user-written
C function G (see description of G below.) SDRIV2 attempts
C to find the value of T at which one of the equations
C changes sign. SDRIV2 can find at most one root per
C equation per internal integration step, and will then
C return the solution either at TOUT or at a root, whichever
C occurs first in the direction of integration. The index
C of the equation whose root is being reported is stored in
C the sixth element of IWORK.
C NOTE: NROOT is never altered by this program.
C
C EPS = On input, the requested relative accuracy in all solution
C components. EPS = 0 is allowed. On output, the adjusted
C relative accuracy if the input value was too small. The
C value of EPS should be set as large as is reasonable,
C because the amount of work done by SDRIV2 increases as
C EPS decreases.
C
C EWT = (Input) Problem zero, i.e., the smallest physically
C meaningful value for the solution. This is used inter-
C nally to compute an array YWT(I) = MAX(ABS(Y(I)), EWT).
C One step error estimates divided by YWT(I) are kept less
C than EPS. Setting EWT to zero provides pure relative
C error control. However, setting EWT smaller than
C necessary can adversely affect the running time.
C
C MINT = (Input) The integration method flag.
C MINT = 1 Means the Adams methods, and is used for
C non-stiff problems.
C MINT = 2 Means the stiff methods of Gear (i.e., the
C backward differentiation formulas), and is
C used for stiff problems.
C MINT = 3 Means the program dynamically selects the
C Adams methods when the problem is non-stiff
C and the Gear methods when the problem is
C stiff.
C MINT may not be changed without restarting, i.e., setting
C the magnitude of MSTATE to 1.
C
C WORK
C LENW = (Input)
C WORK is an array of LENW real words used
C internally for temporary storage. The user must allocate
C space for this array in the calling program by a statement
C such as
C REAL WORK(...)
C The length of WORK should be at least
C 16*N + 2*NROOT + 204 if MINT is 1, or
C N*N + 10*N + 2*NROOT + 204 if MINT is 2, or
C N*N + 17*N + 2*NROOT + 204 if MINT is 3,
C and LENW should be set to the value used. The contents of
C WORK should not be disturbed between calls to SDRIV2.
C
C IWORK
C LENIW = (Input)
C IWORK is an integer array of length LENIW used internally
C for temporary storage. The user must allocate space for
C this array in the calling program by a statement such as
C INTEGER IWORK(...)
C The length of IWORK should be at least
C 21 if MINT is 1, or
C N+21 if MINT is 2 or 3,
C and LENIW should be set to the value used. The contents
C of IWORK should not be disturbed between calls to SDRIV2.
C
C G = A real FORTRAN function supplied by the user
C if NROOT is not 0. In this case, the name must be
C declared EXTERNAL in the user's calling program. G is
C repeatedly called with different values of IROOT to
C obtain the value of each of the NROOT equations for which
C a root is desired. G is of the form:
C REAL FUNCTION G (N, T, Y, IROOT)
C REAL Y(*)
C GO TO (10, ...), IROOT
C 10 G = ...
C .
C .
C END (Sample)
C Here, Y is a vector of length at least N, whose first N
C components are the solution components at the point T.
C The user should not alter these values. The actual length
C of Y is determined by the user's declaration in the
C program which calls SDRIV2. Thus the dimensioning of Y in
C G, while required by FORTRAN convention, does not actually
C allocate any storage. Normally a return from G passes
C control back to SDRIV2. However, if the user would like
C to abort the calculation, i.e., return control to the
C program which calls SDRIV2, he should set N to zero.
C SDRIV2 will signal this by returning a value of MSTATE
C equal to +7(-7). In this case, the index of the equation
C being evaluated is stored in the sixth element of IWORK.
C Altering the value of N in G has no effect on the value of
C N in the call sequence of SDRIV2.
C
C***LONG DESCRIPTION
C
C III. OTHER COMMUNICATION TO THE USER ..............................
C
C A. The solver communicates to the user through the parameters
C above. In addition it writes diagnostic messages through the
C standard error handling program XERROR. That program will
C terminate the user's run if it detects a probable problem setup
C error, e.g., insufficient storage allocated by the user for the
C WORK array. Messages are written on the standard error message
C file. At installations which have this error handling package
C the user should determine the standard error handling file from
C the local documentation. Otherwise the short but serviceable
C routine, XERROR, available with this package, can be used. That
C program writes on logical unit 6 to transmit messages. A
C complete description of XERROR is given in the Sandia
C Laboratories report SAND78-1189 by R. E. Jones.
C
C B. The first three elements of WORK and the first five elements of
C IWORK will contain the following statistical data:
C AVGH The average step size used.
C HUSED The step size last used (successfully).
C AVGORD The average order used.
C IMXERR The index of the element of the solution vector that
C contributed most to the last error test.
C NQUSED The order last used (successfully).
C NSTEP The number of steps taken since last initialization.
C NFE The number of evaluations of the right hand side.
C NJE The number of evaluations of the Jacobian matrix.
C
C IV. REMARKS .......................................................
C
C A. On any return from SDRIV2 all information necessary to continue
C the calculation is contained in the call sequence parameters,
C including the work arrays. Thus it is possible to suspend one
C problem, integrate another, and then return to the first.
C
C B. If this package is to be used in an overlay situation, the user
C must declare in the primary overlay the variables in the call
C sequence to SDRIV2.
C
C C. When the routine G is not required, difficulties associated with
C an unsatisfied external can be avoided by using the name of the
C routine which calculates the right hand side of the differential
C equations in place of G in the call sequence of SDRIV2.
C
C V. USAGE ..........................................................
C
C PROGRAM SAMPLE
C EXTERNAL F
C PARAMETER(MINT = 1, NROOT = 0, N = ...,
C 8 LENW = 16*N + 2*NROOT + 204, LENIW = 21)
C N is the number of equations
C REAL EPS, EWT, T, TOUT, WORK(LENW), Y(N)
C INTEGER IWORK(LENIW)
C OPEN(FILE='TAPE6', UNIT=6, STATUS='NEW')
C T = 0. Initial point
C DO 10 I = 1,N
C 10 Y(I) = ... Set initial conditions
C TOUT = T
C EWT = ...
C MSTATE = 1
C EPS = ...
C 20 CALL SDRIV2 (N, T, Y, F, TOUT, MSTATE, NROOT, EPS, EWT,
C 8 MINT, WORK, LENW, IWORK, LENIW, F)
C Last argument is not the same
C as F if rootfinding is used.
C IF (MSTATE .GT. 2) STOP
C WRITE(6, 100) TOUT, (Y(I), I=1,N)
C TOUT = TOUT + 1.
C IF (TOUT .LE. 10.) GO TO 20
C 100 FORMAT(...)
C END (Sample)
C
C***REFERENCES GEAR, C. W., "NUMERICAL INITIAL VALUE PROBLEMS IN
C ORDINARY DIFFERENTIAL EQUATIONS", PRENTICE-HALL, 1971.
C***ROUTINES CALLED SDRIV3,XERROR
C***END PROLOGUE SDRIV2
EXTERNAL F, G
REAL EPS, EWT, EWTCOM(1), G, HMAX, T, TOUT,
8 WORK(*), Y(*)
INTEGER IWORK(*)
CHARACTER MSG*81
PARAMETER(IMPL = 0, MXSTEP = 1000)
C***FIRST EXECUTABLE STATEMENT SDRIV2
IF (MINT .LT. 1 .OR. MINT .GT. 3) THEN
WRITE(MSG, '(''SDRIV21FE Illegal input. Improper value for '',
8 ''the integration method flag,'', I8)') MINT
CALL XERROR(MSG(1:81), 81, 21, 2)
RETURN
END IF
IF (MSTATE .GE. 0) THEN
NSTATE = MSTATE
NTASK = 1
ELSE
NSTATE = - MSTATE
NTASK = 3
END IF
EWTCOM(1) = EWT
IF (EWT .NE. 0.E0) THEN
IERROR = 3
ELSE
IERROR = 2
END IF
IF (MINT .EQ. 1) THEN
MITER = 0
MXORD = 12
ELSE IF (MINT .EQ. 2) THEN
MITER = 2
MXORD = 5
ELSE IF (MINT .EQ. 3) THEN
MITER = 2
MXORD = 12
END IF
HMAX = 2.E0*ABS(TOUT - T)
CALL SDRIV3 (N, T, Y, F, NSTATE, TOUT, NTASK, NROOT, EPS, EWTCOM,
8 IERROR, MINT, MITER, IMPL, ML, MU, MXORD, HMAX, WORK,
8 LENW, IWORK, LENIW, F, F, NDE, MXSTEP, G, F)
IF (MSTATE .GE. 0) THEN
MSTATE = NSTATE
ELSE
MSTATE = - NSTATE
END IF
END
SUBROUTINE SDCOR (DFDY,EL,FA,H,IMPL,IPVT,MATDIM,MITER,ML,MU,N,
8 NDE,NQ,T,USERS,Y,YH,YWT,EVALFA,SAVE1,SAVE2,A,D,JSTATE)
C***BEGIN PROLOGUE SDCOR
C***REFER TO SDRIV3
C Subroutine SDCOR is called to compute corrections to the Y array.
C In the case of functional iteration, update Y directly from the
C result of the last call to F.
C In the case of the chord method, compute the corrector error and
C solve the linear system with that as right hand side and DFDY as
C coefficient matrix, using the LU decomposition if MITER is 1, 2, 4,
C or 5.
C***ROUTINES CALLED SGESL,SGBSL,SNRM2
C***DATE WRITTEN 790601 (YYMMDD)
C***REVISION DATE 870401 (YYMMDD)
C***CATEGORY NO. I1A2,I1A1B
C***AUTHOR KAHANER, D. K., NATIONAL BUREAU OF STANDARDS,
C SUTHERLAND, C. D., LOS ALAMOS NATIONAL LABORATORY
C***END PROLOGUE SDCOR
REAL A(MATDIM,*), D, DFDY(MATDIM,*), EL(13,12), H,
8 SAVE1(*), SAVE2(*), SNRM2, T, Y(*), YH(N,*), YWT(*)
INTEGER IPVT(*)
LOGICAL EVALFA
C***FIRST EXECUTABLE STATEMENT SDCOR
IF (MITER .EQ. 0) THEN
DO 100 I = 1,N
100 SAVE1(I) = (H*SAVE2(I) - YH(I,2) - SAVE1(I))/YWT(I)
D = SNRM2(N, SAVE1, 1)/SQRT(REAL(N))
DO 105 I = 1,N
105 SAVE1(I) = H*SAVE2(I) - YH(I,2)
ELSE IF (MITER .EQ. 1 .OR. MITER .EQ. 2) THEN
IF (IMPL .EQ. 0) THEN
DO 130 I = 1,N
130 SAVE2(I) = H*SAVE2(I) - YH(I,2) - SAVE1(I)
ELSE IF (IMPL .EQ. 1) THEN
IF (EVALFA) THEN
CALL FA (N, T, Y, A, MATDIM, ML, MU, NDE)
IF (N .EQ. 0) THEN
JSTATE = 9
RETURN
END IF
ELSE
EVALFA = .TRUE.
END IF
DO 150 I = 1,N
150 SAVE2(I) = H*SAVE2(I)
DO 160 J = 1,N
DO 160 I = 1,N
160 SAVE2(I) = SAVE2(I) - A(I,J)*(YH(J,2) + SAVE1(J))
ELSE IF (IMPL .EQ. 2) THEN
IF (EVALFA) THEN
CALL FA (N, T, Y, A, MATDIM, ML, MU, NDE)
IF (N .EQ. 0) THEN
JSTATE = 9
RETURN
END IF
ELSE
EVALFA = .TRUE.
END IF
DO 180 I = 1,N
180 SAVE2(I) = H*SAVE2(I) - A(I,1)*(YH(I,2) + SAVE1(I))
END IF
CALL SGESL (DFDY, MATDIM, N, IPVT, SAVE2, 0)
DO 200 I = 1,N
SAVE1(I) = SAVE1(I) + SAVE2(I)
200 SAVE2(I) = SAVE2(I)/YWT(I)
D = SNRM2(N, SAVE2, 1)/SQRT(REAL(N))
ELSE IF (MITER .EQ. 4 .OR. MITER .EQ. 5) THEN
IF (IMPL .EQ. 0) THEN
DO 230 I = 1,N
230 SAVE2(I) = H*SAVE2(I) - YH(I,2) - SAVE1(I)
ELSE IF (IMPL .EQ. 1) THEN
IF (EVALFA) THEN
CALL FA (N, T, Y, A(ML+1,1), MATDIM, ML, MU, NDE)
IF (N .EQ. 0) THEN
JSTATE = 9
RETURN
END IF
ELSE
EVALFA = .TRUE.
END IF
DO 250 I = 1,N
250 SAVE2(I) = H*SAVE2(I)
MW = ML + 1 + MU
DO 260 J = 1,N
I1 = MAX(ML+1, MW+1-J)
I2 = MIN(MW+N-J, MW+ML)
DO 260 I = I1,I2
I3 = I + J - MW
260 SAVE2(I3) = SAVE2(I3) - A(I,J)*(YH(J,2) + SAVE1(J))
ELSE IF (IMPL .EQ. 2) THEN
IF (EVALFA) THEN
CALL FA (N, T, Y, A, MATDIM, ML, MU, NDE)
IF (N .EQ. 0) THEN
JSTATE = 9
RETURN
END IF
ELSE
EVALFA = .TRUE.
END IF
DO 280 I = 1,N
280 SAVE2(I) = H*SAVE2(I) - A(I,1)*(YH(I,2) + SAVE1(I))
END IF
CALL SGBSL (DFDY, MATDIM, N, ML, MU, IPVT, SAVE2, 0)
DO 300 I = 1,N
SAVE1(I) = SAVE1(I) + SAVE2(I)
300 SAVE2(I) = SAVE2(I)/YWT(I)
D = SNRM2(N, SAVE2, 1)/SQRT(REAL(N))
ELSE IF (MITER .EQ. 3) THEN
IFLAG = 2
CALL USERS (Y, YH(1,2), YWT, SAVE1, SAVE2, T, H, EL(1,NQ), IMPL,
8 N, NDE, IFLAG)
IF (N .EQ. 0) THEN
JSTATE = 10
RETURN
END IF
DO 320 I = 1,N
SAVE1(I) = SAVE1(I) + SAVE2(I)
320 SAVE2(I) = SAVE2(I)/YWT(I)
D = SNRM2(N, SAVE2, 1)/SQRT(REAL(N))
END IF
END
SUBROUTINE SDCST (MAXORD,MINT,ISWFLG,EL,TQ)
C***BEGIN PROLOGUE SDCST
C***REFER TO SDRIV3
C SDCST is called by SDNTL and sets coefficients used by the core
C integrator SDSTP. The array EL determines the basic method.
C The array TQ is involved in adjusting the step size in relation
C to truncation error. EL and TQ depend upon MINT, and are calculated
C for orders 1 to MAXORD(.LE. 12). For each order NQ, the coefficients
C EL are calculated from the generating polynomial:
C L(T) = EL(1,NQ) + EL(2,NQ)*T + ... + EL(NQ+1,NQ)*T**NQ.
C For the implicit Adams methods, L(T) is given by
C dL/dT = (1+T)*(2+T)* ... *(NQ-1+T)/K, L(-1) = 0,
C where K = factorial(NQ-1).
C For the Gear methods,
C L(T) = (1+T)*(2+T)* ... *(NQ+T)/K,
C where K = factorial(NQ)*(1 + 1/2 + ... + 1/NQ).
C For each order NQ, there are three components of TQ.
C***ROUTINES CALLED (NONE)
C***DATE WRITTEN 790601 (YYMMDD)
C***REVISION DATE 870216 (YYMMDD)
C***CATEGORY NO. I1A2,I1A1B
C***AUTHOR KAHANER, D. K., NATIONAL BUREAU OF STANDARDS,
C SUTHERLAND, C. D., LOS ALAMOS NATIONAL LABORATORY
C***END PROLOGUE SDCST
REAL EL(13,12), FACTRL(12), GAMMA(14), SUM, TQ(3,12)
C***FIRST EXECUTABLE STATEMENT SDCST
FACTRL(1) = 1.E0
DO 10 I = 2,MAXORD
10 FACTRL(I) = REAL(I)*FACTRL(I-1)
C COMPUTE ADAMS COEFFICIENTS
IF (MINT .EQ. 1) THEN
GAMMA(1) = 1.E0
DO 40 I = 1,MAXORD+1
SUM = 0.E0
DO 30 J = 1,I
30 SUM = SUM - GAMMA(J)/REAL(I-J+2)
40 GAMMA(I+1) = SUM
EL(1,1) = 1.E0
EL(2,1) = 1.E0
EL(2,2) = 1.E0
EL(3,2) = 1.E0
DO 60 J = 3,MAXORD
EL(2,J) = FACTRL(J-1)
DO 50 I = 3,J
50 EL(I,J) = REAL(J-1)*EL(I,J-1) + EL(I-1,J-1)
60 EL(J+1,J) = 1.E0
DO 80 J = 2,MAXORD
EL(1,J) = EL(1,J-1) + GAMMA(J)
EL(2,J) = 1.E0
DO 80 I = 3,J+1
80 EL(I,J) = EL(I,J)/(REAL(I-1)*FACTRL(J-1))
DO 100 J = 1,MAXORD
TQ(1,J) = -1.E0/(FACTRL(J)*GAMMA(J))
TQ(2,J) = -1.E0/GAMMA(J+1)
100 TQ(3,J) = -1.E0/GAMMA(J+2)
C COMPUTE GEAR COEFFICIENTS
ELSE IF (MINT .EQ. 2) THEN
EL(1,1) = 1.E0
EL(2,1) = 1.E0
DO 130 J = 2,MAXORD
EL(1,J) = FACTRL(J)
DO 120 I = 2,J
120 EL(I,J) = REAL(J)*EL(I,J-1) + EL(I-1,J-1)
130 EL(J+1,J) = 1.E0
SUM = 1.E0
DO 150 J = 2,MAXORD
SUM = SUM + 1.E0/REAL(J)
DO 150 I = 1,J+1
150 EL(I,J) = EL(I,J)/(FACTRL(J)*SUM)
DO 170 J = 1,MAXORD
IF (J .GT. 1) TQ(1,J) = 1.E0/FACTRL(J-1)
TQ(2,J) = REAL(J+1)/EL(1,J)
170 TQ(3,J) = REAL(J+2)/EL(1,J)
END IF
C Compute constants used in the stiffness test.
C These are the ratio of TQ(2,NQ) for the Gear
C methods to those for the Adams methods.
IF (ISWFLG .EQ. 3) THEN
MXRD = MIN(MAXORD, 5)
IF (MINT .EQ. 2) THEN
GAMMA(1) = 1.E0
DO 190 I = 1,MXRD
SUM = 0.E0
DO 180 J = 1,I
180 SUM = SUM - GAMMA(J)/REAL(I-J+2)
190 GAMMA(I+1) = SUM
END IF
SUM = 1.E0
DO 200 I = 2,MXRD
SUM = SUM + 1.E0/REAL(I)
200 EL(1+I,1) = -REAL(I+1)*SUM*GAMMA(I+1)
END IF
END
SUBROUTINE SDNTL (EPS,F,FA,HMAX,HOLD,IMPL,JTASK,MATDIM,MAXORD,
8 MINT,MITER,ML,MU,N,NDE,SAVE1,T,UROUND,USERS,Y,YWT,H,MNTOLD,
8 MTROLD,NFE,RC,YH,A,CONVRG,EL,FAC,IER,IPVT,NQ,NWAIT,RH,RMAX,
8 SAVE2,TQ,TREND,ISWFLG,JSTATE)
C***BEGIN PROLOGUE SDNTL
C***REFER TO SDRIV3
C Subroutine SDNTL is called to set parameters on the first call
C to SDSTP, on an internal restart, or when the user has altered
C MINT, MITER, and/or H.
C On the first call, the order is set to 1 and the initial derivatives
C are calculated. RMAX is the maximum ratio by which H can be
C increased in one step. It is initially RMINIT to compensate
C for the small initial H, but then is normally equal to RMNORM.
C If a failure occurs (in corrector convergence or error test), RMAX
C is set at RMFAIL for the next increase.
C If the caller has changed MINT, or if JTASK = 0, SDCST is called
C to set the coefficients of the method. If the caller has changed H,
C YH must be rescaled. If H or MINT has been changed, NWAIT is
C reset to NQ + 2 to prevent further increases in H for that many
C steps. Also, RC is reset. RC is the ratio of new to old values of
C the coefficient L(0)*H. If the caller has changed MITER, RC is
C set to 0 to force the partials to be updated, if partials are used.
C***ROUTINES CALLED SDCST,SDSCL,SGEFA,SGESL,SGBFA,SGBSL,SNRM2
C***DATE WRITTEN 790601 (YYMMDD)
C***REVISION DATE 870810 (YYMMDD)
C***CATEGORY NO. I1A2,I1A1B
C***AUTHOR KAHANER, D. K., NATIONAL BUREAU OF STANDARDS,
C SUTHERLAND, C. D., LOS ALAMOS NATIONAL LABORATORY
C***END PROLOGUE SDNTL
REAL A(MATDIM,*), EL(13,12), EPS, FAC(*), H, HMAX,
8 HOLD, OLDL0, RC, RH, RMAX, RMINIT, SAVE1(*), SAVE2(*), SMAX,
8 SMIN, SNRM2, SUM, SUM0, T, TQ(3,12), TREND, UROUND, Y(*),
8 YH(N,*), YWT(*)
INTEGER IPVT(*)
LOGICAL CONVRG, IER
PARAMETER(RMINIT = 10000.E0)
C***FIRST EXECUTABLE STATEMENT SDNTL
IER = .FALSE.
IF (JTASK .GE. 0) THEN
IF (JTASK .EQ. 0) THEN
CALL SDCST (MAXORD, MINT, ISWFLG, EL, TQ)
RMAX = RMINIT
END IF
RC = 0.E0
CONVRG = .FALSE.
TREND = 1.E0
NQ = 1
NWAIT = 3
CALL F (N, T, Y, SAVE2)
IF (N .EQ. 0) THEN
JSTATE = 6
RETURN
END IF
NFE = NFE + 1
IF (IMPL .NE. 0) THEN
IF (MITER .EQ. 3) THEN
IFLAG = 0
CALL USERS (Y, YH, YWT, SAVE1, SAVE2, T, H, EL, IMPL, N,
8 NDE, IFLAG)
IF (N .EQ. 0) THEN
JSTATE = 10
RETURN
END IF
ELSE IF (IMPL .EQ. 1) THEN
IF (MITER .EQ. 1 .OR. MITER .EQ. 2) THEN
CALL FA (N, T, Y, A, MATDIM, ML, MU, NDE)
IF (N .EQ. 0) THEN
JSTATE = 9
RETURN
END IF
CALL SGEFA (A, MATDIM, N, IPVT, INFO)
IF (INFO .NE. 0) THEN
IER = .TRUE.
RETURN
END IF
CALL SGESL (A, MATDIM, N, IPVT, SAVE2, 0)
ELSE IF (MITER .EQ. 4 .OR. MITER .EQ. 5) THEN
CALL FA (N, T, Y, A(ML+1,1), MATDIM, ML, MU, NDE)
IF (N .EQ. 0) THEN
JSTATE = 9
RETURN
END IF
CALL SGBFA (A, MATDIM, N, ML, MU, IPVT, INFO)
IF (INFO .NE. 0) THEN
IER = .TRUE.
RETURN
END IF
CALL SGBSL (A, MATDIM, N, ML, MU, IPVT, SAVE2, 0)
END IF
ELSE IF (IMPL .EQ. 2) THEN
CALL FA (N, T, Y, A, MATDIM, ML, MU, NDE)
IF (N .EQ. 0) THEN
JSTATE = 9
RETURN
END IF
DO 150 I = 1,NDE
IF (A(I,1) .EQ. 0.E0) THEN
IER = .TRUE.
RETURN
ELSE
SAVE2(I) = SAVE2(I)/A(I,1)
END IF
150 CONTINUE
DO 155 I = NDE+1,N
155 A(I,1) = 0.E0
END IF
END IF
DO 170 I = 1,NDE
170 SAVE1(I) = SAVE2(I)/YWT(I)
SUM = SNRM2(NDE, SAVE1, 1)
SUM0 = 1.E0/MAX(1.E0, ABS(T))
SMAX = MAX(SUM0, SUM)
SMIN = MIN(SUM0, SUM)
SUM = SMAX*SQRT(1.E0 + (SMIN/SMAX)**2)/SQRT(REAL(NDE))
H = SIGN(MIN(2.E0*EPS/SUM, ABS(H)), H)
DO 180 I = 1,N
180 YH(I,2) = H*SAVE2(I)
IF (MITER .EQ. 2 .OR. MITER .EQ. 5 .OR. ISWFLG .EQ. 3) THEN
DO 20 I = 1,N
20 FAC(I) = SQRT(UROUND)
END IF
ELSE
IF (MITER .NE. MTROLD) THEN
MTROLD = MITER
RC = 0.E0
CONVRG = .FALSE.
END IF
IF (MINT .NE. MNTOLD) THEN
MNTOLD = MINT
OLDL0 = EL(1,NQ)
CALL SDCST (MAXORD, MINT, ISWFLG, EL, TQ)
RC = RC*EL(1,NQ)/OLDL0
NWAIT = NQ + 2
END IF
IF (H .NE. HOLD) THEN
NWAIT = NQ + 2
RH = H/HOLD
CALL SDSCL (HMAX, N, NQ, RMAX, HOLD, RC, RH, YH)
END IF
END IF
END
SUBROUTINE SDNTP (H,K,N,NQ,T,TOUT,YH,Y)
C***BEGIN PROLOGUE SDNTP
C***REFER TO SDRIV3
C Subroutine SDNTP interpolates the K-th derivative of Y at TOUT,
C using the data in the YH array. If K has a value greater than NQ,
C the NQ-th derivative is calculated.
C***ROUTINES CALLED (NONE)
C***DATE WRITTEN 790601 (YYMMDD)
C***REVISION DATE 870216 (YYMMDD)
C***CATEGORY NO. I1A2,I1A1B
C***AUTHOR KAHANER, D. K., NATIONAL BUREAU OF STANDARDS,
C SUTHERLAND, C. D., LOS ALAMOS NATIONAL LABORATORY
C***END PROLOGUE SDNTP
REAL FACTOR, H, R, T, TOUT, Y(*), YH(N,*)
C***FIRST EXECUTABLE STATEMENT SDNTP
IF (K .EQ. 0) THEN
DO 10 I = 1,N
10 Y(I) = YH(I,NQ+1)
R = ((TOUT - T)/H)
DO 20 JJ = 1,NQ
J = NQ + 1 - JJ
DO 20 I = 1,N
20 Y(I) = YH(I,J) + R*Y(I)
ELSE
KUSED = MIN(K, NQ)
FACTOR = 1.E0
DO 40 KK = 1,KUSED
40 FACTOR = FACTOR*REAL(NQ+1-KK)
DO 50 I = 1,N
50 Y(I) = FACTOR*YH(I,NQ+1)
DO 80 JJ = KUSED+1,NQ
J = K + 1 + NQ - JJ
FACTOR = 1.E0
DO 60 KK = 1,KUSED
60 FACTOR = FACTOR*REAL(J-KK)
DO 70 I = 1,N
70 Y(I) = FACTOR*YH(I,J) + R*Y(I)
80 CONTINUE
DO 100 I = 1,N
100 Y(I) = Y(I)*H**(-KUSED)
END IF
END
SUBROUTINE SDPSC (KSGN,N,NQ,YH)
C***BEGIN PROLOGUE SDPSC
C***REFER TO SDRIV3
C This subroutine computes the predicted YH values by effectively
C multiplying the YH array by the Pascal triangle matrix when KSGN
C is +1, and performs the inverse function when KSGN is -1.
C***ROUTINES CALLED (NONE)
C***DATE WRITTEN 790601 (YYMMDD)
C***REVISION DATE 841119 (YYMMDD)
C***CATEGORY NO. I1A2,I1A1B
C***AUTHOR KAHANER, D. K., NATIONAL BUREAU OF STANDARDS,
C SUTHERLAND, C. D., LOS ALAMOS NATIONAL LABORATORY
C***END PROLOGUE SDPSC
REAL YH(N,*)
C***FIRST EXECUTABLE STATEMENT SDPSC
IF (KSGN .GT. 0) THEN
DO 10 J1 = 1,NQ
DO 10 J2 = J1,NQ
J = NQ - J2 + J1
DO 10 I = 1,N
10 YH(I,J) = YH(I,J) + YH(I,J+1)
ELSE
DO 30 J1 = 1,NQ
DO 30 J2 = J1,NQ
J = NQ - J2 + J1
DO 30 I = 1,N
30 YH(I,J) = YH(I,J) - YH(I,J+1)
END IF
END
SUBROUTINE SDPST (EL,F,FA,H,IMPL,JACOBN,MATDIM,MITER,ML,MU,N,NDE,
8 NQ,SAVE2,T,USERS,Y,YH,YWT,UROUND,NFE,NJE,A,DFDY,FAC,IER,IPVT,
8 SAVE1,ISWFLG,BND,JSTATE)
C***BEGIN PROLOGUE SDPST
C***REFER TO SDRIV3
C Subroutine SDPST is called to reevaluate the partials.
C If MITER is 1, 2, 4, or 5, the matrix
C P = I - L(0)*H*Jacobian is stored in DFDY and subjected to LU
C decomposition, with the results also stored in DFDY.
C***ROUTINES CALLED SGEFA,SGBFA,SNRM2
C***DATE WRITTEN 790601 (YYMMDD)
C***REVISION DATE 870401 (YYMMDD)
C***CATEGORY NO. I1A2,I1A1B
C***AUTHOR KAHANER, D. K., NATIONAL BUREAU OF STANDARDS,
C SUTHERLAND, C. D., LOS ALAMOS NATIONAL LABORATORY
C***END PROLOGUE SDPST
REAL A(MATDIM,*), BL, BND, BP, BR, BU, DFDY(MATDIM,*),
8 DFDYMX, DIFF, DY, EL(13,12), FAC(*), FACMAX, FACMIN, FACTOR,
8 H, SAVE1(*), SAVE2(*), SCALE, SNRM2, T, UROUND, Y(*),
8 YH(N,*), YJ, YS, YWT(*)
INTEGER IPVT(*)
LOGICAL IER
PARAMETER(FACMAX = .5E0)
C***FIRST EXECUTABLE STATEMENT SDPST
NJE = NJE + 1
IER = .FALSE.
IF (MITER .EQ. 1 .OR. MITER .EQ. 2) THEN
IF (MITER .EQ. 1) THEN
CALL JACOBN (N, T, Y, DFDY, MATDIM, ML, MU)
IF (N .EQ. 0) THEN
JSTATE = 8
RETURN
END IF
IF (ISWFLG .EQ. 3) BND = SNRM2(N*N, DFDY, 1)
FACTOR = -EL(1,NQ)*H
DO 110 J = 1,N
DO 110 I = 1,N
110 DFDY(I,J) = FACTOR*DFDY(I,J)
ELSE IF (MITER .EQ. 2) THEN
BR = UROUND**(.875E0)
BL = UROUND**(.75E0)
BU = UROUND**(.25E0)
BP = UROUND**(-.15E0)
FACMIN = UROUND**(.78E0)
DO 170 J = 1,N
YS = MAX(ABS(YWT(J)), ABS(Y(J)))
120 DY = FAC(J)*YS
IF (DY .EQ. 0.E0) THEN
IF (FAC(J) .LT. FACMAX) THEN
FAC(J) = MIN(100.E0*FAC(J), FACMAX)
GO TO 120
ELSE
DY = YS
END IF
END IF
IF (NQ .EQ. 1) THEN
DY = SIGN(DY, SAVE2(J))
ELSE
DY = SIGN(DY, YH(J,3))
END IF
DY = (Y(J) + DY) - Y(J)
YJ = Y(J)
Y(J) = Y(J) + DY
CALL F (N, T, Y, SAVE1)
IF (N .EQ. 0) THEN
JSTATE = 6
RETURN
END IF
Y(J) = YJ
FACTOR = -EL(1,NQ)*H/DY
DO 140 I = 1,N
140 DFDY(I,J) = (SAVE1(I) - SAVE2(I))*FACTOR
C Step 1
DIFF = ABS(SAVE2(1) - SAVE1(1))
IMAX = 1
DO 150 I = 2,N
IF (ABS(SAVE2(I) - SAVE1(I)) .GT. DIFF) THEN
IMAX = I
DIFF = ABS(SAVE2(I) - SAVE1(I))
END IF
150 CONTINUE
C Step 2
IF (MIN(ABS(SAVE2(IMAX)), ABS(SAVE1(IMAX))) .GT. 0.E0) THEN
SCALE = MAX(ABS(SAVE2(IMAX)), ABS(SAVE1(IMAX)))
C Step 3
IF (DIFF .GT. BU*SCALE) THEN
FAC(J) = MAX(FACMIN, FAC(J)*.1E0)
ELSE IF (BR*SCALE .LE. DIFF .AND. DIFF .LE. BL*SCALE) THEN
FAC(J) = MIN(FAC(J)*10.E0, FACMAX)
C Step 4
ELSE IF (DIFF .LT. BR*SCALE) THEN
FAC(J) = MIN(BP*FAC(J), FACMAX)
END IF
END IF
170 CONTINUE
IF (ISWFLG .EQ. 3) BND = SNRM2(N*N, DFDY, 1)/(-EL(1,NQ)*H)
NFE = NFE + N
END IF
IF (IMPL .EQ. 0) THEN
DO 190 I = 1,N
190 DFDY(I,I) = DFDY(I,I) + 1.E0
ELSE IF (IMPL .EQ. 1) THEN
CALL FA (N, T, Y, A, MATDIM, ML, MU, NDE)
IF (N .EQ. 0) THEN
JSTATE = 9
RETURN
END IF
DO 210 J = 1,N
DO 210 I = 1,N
210 DFDY(I,J) = DFDY(I,J) + A(I,J)
ELSE IF (IMPL .EQ. 2) THEN
CALL FA (N, T, Y, A, MATDIM, ML, MU, NDE)
IF (N .EQ. 0) THEN
JSTATE = 9
RETURN
END IF
DO 230 I = 1,NDE
230 DFDY(I,I) = DFDY(I,I) + A(I,1)
END IF
CALL SGEFA (DFDY, MATDIM, N, IPVT, INFO)
IF (INFO .NE. 0) IER = .TRUE.
ELSE IF (MITER .EQ. 4 .OR. MITER .EQ. 5) THEN
IF (MITER .EQ. 4) THEN
CALL JACOBN (N, T, Y, DFDY(ML+1,1), MATDIM, ML, MU)
IF (N .EQ. 0) THEN
JSTATE = 8
RETURN
END IF
FACTOR = -EL(1,NQ)*H
MW = ML + MU + 1
DO 260 J = 1,N
I1 = MAX(ML+1, MW+1-J)
I2 = MIN(MW+N-J, MW+ML)
DO 260 I = I1,I2
260 DFDY(I,J) = FACTOR*DFDY(I,J)
ELSE IF (MITER .EQ. 5) THEN
BR = UROUND**(.875E0)
BL = UROUND**(.75E0)
BU = UROUND**(.25E0)
BP = UROUND**(-.15E0)
FACMIN = UROUND**(.78E0)
MW = ML + MU + 1
J2 = MIN(MW, N)
DO 340 J = 1,J2
DO 290 K = J,N,MW
YS = MAX(ABS(YWT(K)), ABS(Y(K)))
280 DY = FAC(K)*YS
IF (DY .EQ. 0.E0) THEN
IF (FAC(K) .LT. FACMAX) THEN
FAC(K) = MIN(100.E0*FAC(K), FACMAX)
GO TO 280
ELSE
DY = YS
END IF
END IF
IF (NQ .EQ. 1) THEN
DY = SIGN(DY, SAVE2(K))
ELSE
DY = SIGN(DY, YH(K,3))
END IF
DY = (Y(K) + DY) - Y(K)
DFDY(MW,K) = Y(K)
290 Y(K) = Y(K) + DY
CALL F (N, T, Y, SAVE1)
IF (N .EQ. 0) THEN
JSTATE = 6
RETURN
END IF
DO 330 K = J,N,MW
Y(K) = DFDY(MW,K)
YS = MAX(ABS(YWT(K)), ABS(Y(K)))
DY = FAC(K)*YS
IF (DY .EQ. 0.E0) DY = YS
IF (NQ .EQ. 1) THEN
DY = SIGN(DY, SAVE2(K))
ELSE
DY = SIGN(DY, YH(K,3))
END IF
DY = (Y(K) + DY) - Y(K)
FACTOR = -EL(1,NQ)*H/DY
I1 = MAX(ML+1, MW+1-K)
I2 = MIN(MW+N-K, MW+ML)
DO 300 I = I1,I2
I3 = K + I - MW
300 DFDY(I,K) = FACTOR*(SAVE1(I3) - SAVE2(I3))
C Step 1
IMAX = MAX(1, K - MU)
DIFF = ABS(SAVE2(IMAX) - SAVE1(IMAX))
I1 = IMAX
I2 = MIN(K + ML, N)
DO 310 I = I1+1,I2
IF (ABS(SAVE2(I) - SAVE1(I)) .GT. DIFF) THEN
IMAX = I
DIFF = ABS(SAVE2(I) - SAVE1(I))
END IF
310 CONTINUE
C Step 2
IF (MIN(ABS(SAVE2(IMAX)), ABS(SAVE1(IMAX))) .GT.0.E0) THEN
SCALE = MAX(ABS(SAVE2(IMAX)), ABS(SAVE1(IMAX)))
C Step 3
IF (DIFF .GT. BU*SCALE) THEN
FAC(K) = MAX(FACMIN, FAC(K)*.1E0)
ELSE IF (BR*SCALE .LE.DIFF .AND. DIFF .LE.BL*SCALE) THEN
FAC(K) = MIN(FAC(K)*10.E0, FACMAX)
C Step 4
ELSE IF (DIFF .LT. BR*SCALE) THEN
FAC(K) = MIN(BP*FAC(K), FACMAX)
END IF
END IF
330 CONTINUE
340 CONTINUE
NFE = NFE + J2
END IF
IF (ISWFLG .EQ. 3) THEN
DFDYMX = 0.E0
DO 345 J = 1,N
I1 = MAX(ML+1, MW+1-J)
I2 = MIN(MW+N-J, MW+ML)
DO 345 I = I1,I2
345 DFDYMX = MAX(DFDYMX, ABS(DFDY(I,J)))
BND = 0.E0
IF (DFDYMX .NE. 0.E0) THEN
DO 350 J = 1,N
I1 = MAX(ML+1, MW+1-J)
I2 = MIN(MW+N-J, MW+ML)
DO 350 I = I1,I2
350 BND = BND + (DFDY(I,J)/DFDYMX)**2
BND = DFDYMX*SQRT(BND)/(-EL(1,NQ)*H)
END IF
END IF
IF (IMPL .EQ. 0) THEN
DO 360 J = 1,N
360 DFDY(MW,J) = DFDY(MW,J) + 1.E0
ELSE IF (IMPL .EQ. 1) THEN
CALL FA (N, T, Y, A(ML+1,1), MATDIM, ML, MU, NDE)
IF (N .EQ. 0) THEN
JSTATE = 9
RETURN
END IF
DO 380 J = 1,N
I1 = MAX(ML+1, MW+1-J)
I2 = MIN(MW+N-J, MW+ML)
DO 380 I = I1,I2
380 DFDY(I,J) = DFDY(I,J) + A(I,J)
ELSE IF (IMPL .EQ. 2) THEN
CALL FA (N, T, Y, A, MATDIM, ML, MU, NDE)
IF (N .EQ. 0) THEN
JSTATE = 9
RETURN
END IF
DO 400 J = 1,NDE
400 DFDY(MW,J) = DFDY(MW,J) + A(J,1)
END IF
CALL SGBFA (DFDY, MATDIM, N, ML, MU, IPVT, INFO)
IF (INFO .NE. 0) IER = .TRUE.
ELSE IF (MITER .EQ. 3) THEN
IFLAG = 1
CALL USERS (Y, YH(1,2), YWT, SAVE1, SAVE2, T, H, EL(1,NQ), IMPL,
8 N, NDE, IFLAG)
IF (N .EQ. 0) THEN
JSTATE = 10
RETURN
END IF
END IF
END
SUBROUTINE SDSCL (HMAX,N,NQ,RMAX,H,RC,RH,YH)
C***BEGIN PROLOGUE SDSCL
C***REFER TO SDRIV3
C This subroutine rescales the YH array whenever the step size
C is changed.
C***ROUTINES CALLED (NONE)
C***DATE WRITTEN 790601 (YYMMDD)
C***REVISION DATE 850319 (YYMMDD)
C***CATEGORY NO. I1A2,I1A1B
C***AUTHOR KAHANER, D. K., NATIONAL BUREAU OF STANDARDS,
C SUTHERLAND, C. D., LOS ALAMOS NATIONAL LABORATORY
C***END PROLOGUE SDSCL
REAL H, HMAX, RC, RH, RMAX, R1, YH(N,*)
C***FIRST EXECUTABLE STATEMENT SDSCL
IF (H .LT. 1.E0) THEN
RH = MIN(ABS(H)*RH, ABS(H)*RMAX, HMAX)/ABS(H)
ELSE
RH = MIN(RH, RMAX, HMAX/ABS(H))
END IF
R1 = 1.E0
DO 10 J = 1,NQ
R1 = R1*RH
DO 10 I = 1,N
10 YH(I,J+1) = YH(I,J+1)*R1
H = H*RH
RC = RC*RH
END
SUBROUTINE SDSTP (EPS,F,FA,HMAX,IMPL,JACOBN,MATDIM,MAXORD,MINT,
8 MITER,ML,MU,N,NDE,YWT,UROUND,USERS,AVGH,AVGORD,H,HUSED,JTASK,
8 MNTOLD,MTROLD,NFE,NJE,NQUSED,NSTEP,T,Y,YH,A,CONVRG,DFDY,EL,FAC,
8 HOLD,IPVT,JSTATE,NQ,NWAIT,RC,RMAX,SAVE1,SAVE2,TQ,TREND,ISWFLG,
8 MTRSV,MXRDSV)
C***BEGIN PROLOGUE SDSTP
C***REFER TO SDRIV3
C SDSTP performs one step of the integration of an initial value
C problem for a system of ordinary differential equations.
C Communication with SDSTP 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 SDSTP.
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 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***ROUTINES CALLED SDNTL,SDPST,SDCOR,SDPSC,SDSCL,SNRM2
C***DATE WRITTEN 790601 (YYMMDD)
C***REVISION DATE 870810 (YYMMDD)
C***CATEGORY NO. I1A2,I1A1B
C***AUTHOR KAHANER, D. K., NATIONAL BUREAU OF STANDARDS,
C SUTHERLAND, C. D., LOS ALAMOS NATIONAL LABORATORY
C***END PROLOGUE SDSTP
EXTERNAL F, JACOBN, FA, USERS
REAL 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(*), SNRM2, T, TOLD, TQ(3,12), TREND, TRSHLD,
8 UROUND, Y(*), YH(N,*), YWT(*), Y0NRM
INTEGER IPVT(*)
LOGICAL CONVRG, EVALFA, EVALJC, IER, SWITCH
PARAMETER(BIAS1 = 1.3E0, BIAS2 = 1.2E0, BIAS3 = 1.4E0, MXFAIL = 3,
8 MXITER = 3, MXTRY = 50, RCTEST = .3E0, RMFAIL = 2.E0,
8 RMNORM = 10.E0, TRSHLD = 1.E0)
DATA IER /.FALSE./
C***FIRST EXECUTABLE STATEMENT SDSTP
NSV = N
BND = 0.E0
SWITCH = .FALSE.
NTRY = 0
TOLD = T
NFAIL = 0
IF (JTASK .LE. 0) THEN
CALL SDNTL (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.E0) 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 SDPSC (1, N, NQ, YH)
EVALJC = ((ABS(RC - 1.E0) .GT. RCTEST) .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 SDPST (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.E0
END IF
DO 125 I = 1,N
125 SAVE1(I) = 0.E0
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 SDCOR (DFDY, EL, FA, H, IMPL, IPVT, MATDIM, MITER, ML,
8 MU, N, NDE, NQ, T, USERS, Y, YH, YWT, EVALFA, SAVE1,
8 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 = SNRM2(N, SAVE1, 1)
DO 132 I = 1,N
132 DFDY(1,I) = SAVE1(I)
Y0NRM = SNRM2(N, YH, 1)
ELSE
DENOM = NUMER
DO 134 I = 1,N
134 DFDY(1,I) = SAVE1(I) - DFDY(1,I)
NUMER = SNRM2(N, DFDY, MATDIM)
IF (EL(1,NQ)*NUMER .LE. 100.E0*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.E0)
8 BND = MAX(BND, NUMER/(DENOM*ABS(H)*EL(1,NQ)))
END IF
END IF
IF (ITER .GT. 0) TREND = MAX(.9E0*TREND, D/D1)
D1 = D
CTEST = MIN(2.E0*TREND, 1.E0)*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 SDPSC (-1, N, NQ, YH)
NWAIT = NQ + 2
IF (JTASK .NE. 0 .AND. JTASK .NE. 2) RMAX = RMFAIL
IF (ITER .EQ. 0) THEN
RH = .3E0
ELSE
RH = .9E0*(EPS/CTEST)**(.2E0)
END IF
IF (RH*H .EQ. 0.E0) GO TO 400
CALL SDSCL (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)
DO 180 I = 1,NDE
180 SAVE2(I) = SAVE1(I)/YWT(I)
ETEST = SNRM2(NDE, SAVE2, 1)/(TQ(2,NQ)*SQRT(REAL(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 SDPSC (-1, N, NQ, YH)
NFAIL = NFAIL + 1
IF (NFAIL .LT. MXFAIL) THEN
IF (JTASK .NE. 0 .AND. JTASK .NE. 2) RMAX = RMFAIL
RH2 = 1.E0/(BIAS2*(ETEST/EPS)**(1.E0/REAL(NQ+1)))
IF (NQ .GT. 1) THEN
DO 190 I = 1,NDE
190 SAVE2(I) = YH(I,NQ+1)/YWT(I)
ERDN = SNRM2(NDE, SAVE2, 1)/(TQ(1,NQ)*SQRT(REAL(NDE)))
RH1 = 1.E0/MAX(1.E0, BIAS1*(ERDN/EPS)**(1.E0/REAL(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.E0) GO TO 400
CALL SDSCL (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 SDNTL (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.E0) 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 = (REAL(NSTEP-1)*AVGH + H)/REAL(NSTEP)
AVGORD = (REAL(NSTEP-1)*AVGORD + REAL(NQ))/REAL(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.E0) THEN
IF (MINT .EQ. 1 .AND. NQ .LE. 5) THEN
HN = ABS(H)/MAX(UROUND, (ETEST/EPS)**(1.E0/REAL(NQ+1)))
HN = MIN(HN, 1.E0/(2.E0*EL(1,NQ)*BND))
HS = ABS(H)/MAX(UROUND,
8 (ETEST/(EPS*EL(NQ+1,1)))**(1.E0/REAL(NQ+1)))
IF (HS .GT. 1.2E0*HN) THEN
MINT = 2
MNTOLD = MINT
MITER = MTRSV
MTROLD = MITER
MAXORD = MIN(MXRDSV, 5)
RC = 0.E0
RMAX = RMNORM
TREND = 1.E0
CALL SDCST (MAXORD, MINT, ISWFLG, EL, TQ)
NWAIT = NQ + 2
END IF
ELSE IF (MINT .EQ. 2) THEN
HS = ABS(H)/MAX(UROUND, (ETEST/EPS)**(1.E0/REAL(NQ+1)))
HN = ABS(H)/MAX(UROUND,
8 (ETEST*EL(NQ+1,1)/EPS)**(1.E0/REAL(NQ+1)))
HN = MIN(HN, 1.E0/(2.E0*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.E0
CONVRG = .FALSE.
CALL SDCST (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.E0
RMAX = RMNORM
TREND = 1.E0
CALL SDCST (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.E0/MAX(UROUND, BIAS2*(ETEST/EPS)**(1.E0/REAL(NQ+1)))
IF (RH.GT.TRSHLD) CALL SDSCL (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.E0
ELSE
DO 270 I = 1,NDE
270 SAVE2(I) = YH(I,NQ+1)/YWT(I)
ERDN = SNRM2(NDE, SAVE2, 1)/(TQ(1,NQ)*SQRT(REAL(NDE)))
RH1 = 1.E0/MAX(UROUND, BIAS1*(ERDN/EPS)**(1.E0/REAL(NQ)))
END IF
RH2 = 1.E0/MAX(UROUND, BIAS2*(ETEST/EPS)**(1.E0/REAL(NQ+1)))
IF (NQ .EQ. MAXORD) THEN
RH3 = 0.E0
ELSE
DO 290 I = 1,NDE
290 SAVE2(I) = (SAVE1(I) - YH(I,MAXORD+1))/YWT(I)
ERUP = SNRM2(NDE, SAVE2, 1)/(TQ(3,NQ)*SQRT(REAL(NDE)))
RH3 = 1.E0/MAX(UROUND, BIAS3*(ERUP/EPS)**(1.E0/REAL(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)/REAL(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.E0) RH = MIN(RH, 1.E0/(2.E0*EL(1,NQ)*BND*ABS(H)))
END IF
CALL SDSCL (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 SDPSC (-1, NSV, NQ, YH)
DO 435 I = 1,NSV
435 Y(I) = YH(I,1)
440 HOLD = H
RETURN
END
SUBROUTINE SDZRO (AE,F,H,N,NQ,IROOT,RE,T,YH,UROUND,B,C,FB,FC,Y)
C***BEGIN PROLOGUE SDZRO
C***REFER TO SDRIV3
C This is a special purpose version of ZEROIN, modified for use with
C the SDRIV1 package.
C
C Sandia Mathematical Program Library
C Mathematical Computing Services Division 5422
C Sandia Laboratories
C P. O. Box 5800
C Albuquerque, New Mexico 87115
C Control Data 6600 Version 4.5, 1 November 1971
C
C ABSTRACT
C ZEROIN searches for a zero of a function F(N, T, Y, IROOT)
C between the given values B and C until the width of the
C interval (B, C) has collapsed to within a tolerance specified
C by the stopping criterion, ABS(B - C) .LE. 2.*(RW*ABS(B) + AE).
C
C Description of parameters
C F - Name of the external function, which returns a
C real result. This name must be in an
C EXTERNAL statement in the calling program.
C B - One end of the interval (B, C). The value returned for
C B usually is the better approximation to a zero of F.
C C - The other end of the interval (B, C).
C RE - Relative error used for RW in the stopping criterion.
C If the requested RE is less than machine precision,
C then RW is set to approximately machine precision.
C AE - Absolute error used in the stopping criterion. If the
C given interval (B, C) contains the origin, then a
C nonzero value should be chosen for AE.
C
C REFERENCES
C 1. L F Shampine and H A Watts, ZEROIN, A Root-Solving Routine,
C SC-TM-70-631, Sept 1970.
C 2. T J Dekker, Finding a Zero by Means of Successive Linear
C Interpolation, "Constructive Aspects of the Fundamental
C Theorem of Algebra", edited by B Dejon and P Henrici, 1969.
C***ROUTINES CALLED SDNTP
C***DATE WRITTEN 790601 (YYMMDD)
C***REVISION DATE 870511 (YYMMDD)
C***CATEGORY NO. I1A2,I1A1B
C***AUTHOR KAHANER, D. K., NATIONAL BUREAU OF STANDARDS,
C SUTHERLAND, C. D., LOS ALAMOS NATIONAL LABORATORY
C***END PROLOGUE SDZRO
REAL A, ACBS, ACMB, AE, B, C, CMB, ER, F, FA, FB, FC,
8 H, P, Q, RE, RW, T, TOL, UROUND, Y(*), YH(N,*)
C***FIRST EXECUTABLE STATEMENT SDZRO
ER = 4.E0*UROUND
RW = MAX(RE, ER)
IC = 0
ACBS = ABS(B - C)
A = C
FA = FC
KOUNT = 0
C Perform interchange
10 IF (ABS(FC) .LT. ABS(FB)) THEN
A = B
FA = FB
B = C
FB = FC
C = A
FC = FA
END IF
CMB = 0.5E0*(C - B)
ACMB = ABS(CMB)
TOL = RW*ABS(B) + AE
C Test stopping criterion
IF (ACMB .LE. TOL) RETURN
IF (KOUNT .GT. 50) RETURN
C Calculate new iterate implicitly as
C B + P/Q, where we arrange P .GE. 0.
C The implicit form is used to prevent overflow.
P = (B - A)*FB
Q = FA - FB
IF (P .LT. 0.E0) THEN
P = -P
Q = -Q
END IF
C Update A and check for satisfactory reduction
C in the size of our bounding interval.
A = B
FA = FB
IC = IC + 1
IF (IC .GE. 4) THEN
IF (8.E0*ACMB .GE. ACBS) THEN
C Bisect
B = 0.5E0*(C + B)
GO TO 20
END IF
IC = 0
END IF
ACBS = ACMB
C Test for too small a change
IF (P .LE. ABS(Q)*TOL) THEN
C Increment by tolerance
B = B + SIGN(TOL, CMB)
C Root ought to be between
C B and (C + B)/2.
ELSE IF (P .LT. CMB*Q) THEN
C Interpolate
B = B + P/Q
ELSE
C Bisect
B = 0.5E0*(C + B)
END IF
C Have completed computation
C for new iterate B.
20 CALL SDNTP (H, 0, N, NQ, T, B, YH, Y)
FB = F(N, B, Y, IROOT)
IF (N .EQ. 0) RETURN
IF (FB .EQ. 0.E0) RETURN
KOUNT = KOUNT + 1
C
C Decide whether next step is interpolation or extrapolation
C
IF (SIGN(1.0E0, FB) .EQ. SIGN(1.0E0, FC)) THEN
C = A
FC = FA
END IF
GO TO 10
END
SUBROUTINE SDRIV3 (N,T,Y,F,NSTATE,TOUT,NTASK,NROOT,EPS,EWT,IERROR,
8 MINT,MITER,IMPL,ML,MU,MXORD,HMAX,WORK,LENW,IWORK,LENIW,JACOBN,
8 FA,NDE,MXSTEP,G,USERS)
C***BEGIN PROLOGUE SDRIV3
C THIS PROLOGUE HAS BEEN REMOVED FOR REASONS OF SPACE
C FOR A COMPLETE COPY OF THIS ROUTINE CONTACT THE AUTHORS
C
C From the book "Numerical Methods and Software"
C by D. Kahaner, C. Moler, S. Nash
C Prentice Hall 1988
C
C***END PROLOGUE SDRIV3
EXTERNAL F, JACOBN, FA, G, USERS
REAL AE, BIG, EPS, EWT(*), G, GLAST, H, HMAX, HSIGN,
8 NROUND, RE, R1MACH, SIZE, SNRM2, SUM, T, TLAST, TOUT, TROOT,
8 UROUND, WORK(*), Y(*)
INTEGER IWORK(*)
LOGICAL CONVRG
CHARACTER MSG*205
PARAMETER(NROUND = 20.E0)
PARAMETER(IAVGH = 1, IHUSED = 2, IAVGRD = 3,
8 IEL = 4, IH = 160, IHMAX = 161, IHOLD = 162,
8 IHSIGN = 163, IRC = 164, IRMAX = 165, IT = 166,
8 ITOUT = 167, ITQ = 168, ITREND = 204, IYH = 205,
8 INDMXR = 1, INQUSD = 2, INSTEP = 3, INFE = 4, INJE = 5,
8 INROOT = 6, ICNVRG = 7, IJROOT = 8, IJTASK = 9,
8 IMNTLD = 10, IMTRLD = 11, INQ = 12, INRTLD = 13,
8 INDTRT = 14, INWAIT = 15, IMNT = 16, IMTRSV = 17,
8 IMTR = 18, IMXRDS = 19, IMXORD = 20, INDPRT = 21,
8 INDPVT = 22)
C***FIRST EXECUTABLE STATEMENT SDRIV3
NPAR = N
UROUND = R1MACH (4)
IF (NROOT .NE. 0) THEN
AE = R1MACH(1)
RE = UROUND
END IF
IF (EPS .LT. 0.E0) THEN
WRITE(MSG, '(''SDRIV36FE Illegal input. EPS,'', E16.8,
8 '', is negative.'')') EPS
CALL XERROR(MSG(1:60), 60, 6, 2)
RETURN
END IF
IF (N .LE. 0) THEN
WRITE(MSG, '(''SDRIV37FE Illegal input. Number of equations,'',
8 I8, '', is not positive.'')') N
CALL XERROR(MSG(1:72), 72, 7, 2)
RETURN
END IF
IF (MXORD .LE. 0) THEN
WRITE(MSG, '(''SDRIV314FE Illegal input. Maximum order,'', I8,
8 '', is not positive.'')') MXORD
CALL XERROR(MSG(1:67), 67, 14, 2)
RETURN
END IF
IF ((MINT .LT. 1 .OR. MINT .GT. 3) .OR. (MINT .EQ. 3 .AND.
8 (MITER .EQ. 0 .OR. MITER .EQ. 3 .OR. IMPL .NE. 0))
8 .OR. (MITER .LT. 0 .OR. MITER .GT. 5) .OR.
8 (IMPL .NE. 0 .AND. IMPL .NE. 1 .AND. IMPL .NE. 2) .OR.
8 ((IMPL .EQ. 1 .OR. IMPL .EQ. 2) .AND. MITER .EQ. 0) .OR.
8 (IMPL .EQ. 2 .AND. MINT .EQ. 1) .OR.
8 (NSTATE .LT. 1 .OR. NSTATE .GT. 10)) THEN
WRITE(MSG, '(''SDRIV39FE Illegal input. Improper value for '',
8 ''NSTATE(MSTATE), MINT, MITER or IMPL.'')')
CALL XERROR(MSG(1:81), 81, 9, 2)
RETURN
END IF
IF (MITER .EQ. 0 .OR. MITER .EQ. 3) THEN
LIWCHK = INDPVT - 1
ELSE IF (MITER .EQ. 1 .OR. MITER .EQ. 2 .OR. MITER .EQ. 4 .OR.
8 MITER .EQ. 5) THEN
LIWCHK = INDPVT + N - 1
END IF
IF (LENIW .LT. LIWCHK) THEN
WRITE(MSG, '(''SDRIV310FE Illegal input. Insufficient '',
8 ''storage allocated for the IWORK array. Based on the '')')
WRITE(MSG(94:), '(''value of the input parameters involved, '',
8 ''the required storage is'', I8)') LIWCHK
CALL XERROR(MSG(1:164), 164, 10, 2)
RETURN
END IF
C Allocate the WORK array
C IYH is the index of YH in WORK
IF (MINT .EQ. 1 .OR. MINT .EQ. 3) THEN
MAXORD = MIN(MXORD, 12)
ELSE IF (MINT .EQ. 2) THEN
MAXORD = MIN(MXORD, 5)
END IF
IDFDY = IYH + (MAXORD + 1)*N
C IDFDY is the index of DFDY
C
IF (MITER .EQ. 0 .OR. MITER .EQ. 3) THEN
IYWT = IDFDY
ELSE IF (MITER .EQ. 1 .OR. MITER .EQ. 2) THEN
IYWT = IDFDY + N*N
ELSE IF (MITER .EQ. 4 .OR. MITER .EQ. 5) THEN
IYWT = IDFDY + (2*ML + MU + 1)*N
END IF
C IYWT is the index of YWT
ISAVE1 = IYWT + N
C ISAVE1 is the index of SAVE1
ISAVE2 = ISAVE1 + N
C ISAVE2 is the index of SAVE2
IGNOW = ISAVE2 + N
C IGNOW is the index of GNOW
ITROOT = IGNOW + NROOT
C ITROOT is the index of TROOT
IFAC = ITROOT + NROOT
C IFAC is the index of FAC
IF (MITER .EQ. 2 .OR. MITER .EQ. 5 .OR. MINT .EQ. 3) THEN
IA = IFAC + N
ELSE
IA = IFAC
END IF
C IA is the index of A
IF (IMPL .EQ. 0 .OR. MITER .EQ. 3) THEN
LENCHK = IA - 1
ELSE IF (IMPL .EQ. 1 .AND. (MITER .EQ. 1 .OR. MITER .EQ. 2)) THEN
LENCHK = IA - 1 + N*N
ELSE IF (IMPL .EQ. 1 .AND. (MITER .EQ. 4 .OR. MITER .EQ. 5)) THEN
LENCHK = IA - 1 + (2*ML + MU + 1)*N
ELSE IF (IMPL .EQ. 2 .AND. MITER .NE. 3) THEN
LENCHK = IA - 1 + N
END IF
IF (LENW .LT. LENCHK) THEN
WRITE(MSG, '(''SDRIV38FE Illegal input. Insufficient '',
8 ''storage allocated for the WORK array. Based on the '')')
WRITE(MSG(92:), '(''value of the input parameters involved, '',
8 ''the required storage is'', I8)') LENCHK
CALL XERROR(MSG(1:162), 162, 8, 2)
RETURN
END IF
IF (MITER .EQ. 0 .OR. MITER .EQ. 3) THEN
MATDIM = 1
ELSE IF (MITER .EQ. 1 .OR. MITER .EQ. 2) THEN
MATDIM = N
ELSE IF (MITER .EQ. 4 .OR. MITER .EQ. 5) THEN
MATDIM = 2*ML + MU + 1
END IF
IF (IMPL .EQ. 0 .OR. IMPL .EQ. 1) THEN
NDECOM = N
ELSE IF (IMPL .EQ. 2) THEN
NDECOM = NDE
END IF
IF (NSTATE .EQ. 1) THEN
C Initialize parameters
IF (MINT .EQ. 1 .OR. MINT .EQ. 3) THEN
IWORK(IMXORD) = MIN(MXORD, 12)
ELSE IF (MINT .EQ. 2) THEN
IWORK(IMXORD) = MIN(MXORD, 5)
END IF
IWORK(IMXRDS) = MXORD
IF (MINT .EQ. 1 .OR. MINT .EQ. 2) THEN
IWORK(IMNT) = MINT
IWORK(IMTR) = MITER
IWORK(IMNTLD) = MINT
IWORK(IMTRLD) = MITER
ELSE IF (MINT .EQ. 3) THEN
IWORK(IMNT) = 1
IWORK(IMTR) = 0
IWORK(IMNTLD) = IWORK(IMNT)
IWORK(IMTRLD) = IWORK(IMTR)
IWORK(IMTRSV) = MITER
END IF
WORK(IHMAX) = HMAX
H = (TOUT - T)*(1.E0 - 4.E0*UROUND)
H = SIGN(MIN(ABS(H), HMAX), H)
WORK(IH) = H
HSIGN = SIGN(1.E0, H)
WORK(IHSIGN) = HSIGN
IWORK(IJTASK) = 0
WORK(IAVGH) = 0.E0
WORK(IHUSED) = 0.E0
WORK(IAVGRD) = 0.E0
IWORK(INDMXR) = 0
IWORK(INQUSD) = 0
IWORK(INSTEP) = 0
IWORK(INFE) = 0
IWORK(INJE) = 0
IWORK(INROOT) = 0
WORK(IT) = T
IWORK(ICNVRG) = 0
IWORK(INDPRT) = 0
C Set initial conditions
DO 30 I = 1,N
JYH = I + IYH - 1
30 WORK(JYH) = Y(I)
IF (T .EQ. TOUT) RETURN
GO TO 180
END IF
C On a continuation, check
C that output points have
C been or will be overtaken.
IF (IWORK(ICNVRG) .EQ. 1) THEN
CONVRG = .TRUE.
ELSE
CONVRG = .FALSE.
END IF
T = WORK(IT)
H = WORK(IH)
HSIGN = WORK(IHSIGN)
IF (IWORK(IJTASK) .EQ. 0) GO TO 180
C
C IWORK(IJROOT) flags unreported
C roots, and is set to the value of
C NTASK when a root was last selected.
C It is set to zero when all roots
C have been reported. IWORK(INROOT)
C contains the index and WORK(ITOUT)
C contains the value of the root last
C selected to be reported.
C IWORK(INRTLD) contains the value of
C NROOT and IWORK(INDTRT) contains
C the value of ITROOT when the array
C of roots was last calculated.
IF (NROOT .NE. 0) THEN
JROOT = IWORK(IJROOT)
IF (JROOT .GT. 0) THEN
C TOUT has just been reported.
C If TROOT .LE. TOUT, report TROOT.
IF (NSTATE .NE. 5) THEN
IF (TOUT*HSIGN .GE. WORK(ITOUT)*HSIGN) THEN
TROOT = WORK(ITOUT)
CALL SDNTP(H, 0, N, IWORK(INQ), T, TROOT, WORK(IYH), Y)
T = TROOT
NSTATE = 5
GO TO 580
END IF
C A root has just been reported.
C Select the next root.
ELSE
TROOT = T
IROOT = 0
DO 50 I = 1,IWORK(INRTLD)
JTROOT = IWORK(INDTRT) + I - 1
IF (WORK(JTROOT)*HSIGN .LE. TROOT*HSIGN) THEN
C
C Check for multiple roots.
C
IF (WORK(JTROOT) .EQ. WORK(ITOUT) .AND.
8 I .GT. IWORK(INROOT)) THEN
IROOT = I
TROOT = WORK(JTROOT)
GO TO 60
END IF
IF (WORK(JTROOT)*HSIGN .GT. WORK(ITOUT)*HSIGN) THEN
IROOT = I
TROOT = WORK(JTROOT)
END IF
END IF
50 CONTINUE
60 IWORK(INROOT) = IROOT
WORK(ITOUT) = TROOT
IWORK(IJROOT) = NTASK
IF (NTASK .EQ. 1) THEN
IF (IROOT .EQ. 0) THEN
IWORK(IJROOT) = 0
ELSE
IF (TOUT*HSIGN .GE. TROOT*HSIGN) THEN
CALL SDNTP(H, 0, N, IWORK(INQ), T, TROOT,WORK(IYH),Y)
NSTATE = 5
T = TROOT
GO TO 580
END IF
END IF
ELSE IF (NTASK .EQ. 2 .OR. NTASK .EQ. 3) THEN
C
C If there are no more roots, or the
C user has altered TOUT to be less
C than a root, set IJROOT to zero.
C
IF (IROOT .EQ. 0 .OR. (TOUT*HSIGN .LT. TROOT*HSIGN)) THEN
IWORK(IJROOT) = 0
ELSE
CALL SDNTP(H, 0, N, IWORK(INQ), T, TROOT, WORK(IYH), Y)
NSTATE = 5
T = TROOT
GO TO 580
END IF
END IF
END IF
END IF
END IF
C
IF (NTASK .EQ. 1) THEN
NSTATE = 2
IF (T*HSIGN .GE. TOUT*HSIGN) THEN
CALL SDNTP (H, 0, N, IWORK(INQ), T, TOUT, WORK(IYH), Y)
T = TOUT
GO TO 580
END IF
ELSE IF (NTASK .EQ. 2) THEN
C Check if TOUT has
C been reset .LT. T
IF (T*HSIGN .GT. TOUT*HSIGN) THEN
WRITE(MSG, '(''SDRIV32WRN With NTASK='', I1, '' on input, '',
8 ''T,'', E16.8, '', was beyond TOUT,'', E16.8, ''. Solution'',
8 '' obtained by interpolation.'')') NTASK, T, TOUT
CALL XERROR(MSG(1:124), 124, 2, 0)
CALL SDNTP (H, 0, N, IWORK(INQ), T, TOUT, WORK(IYH), Y)
T = TOUT
NSTATE = 2
GO TO 580
END IF
C Determine if TOUT has been overtaken
C
IF (ABS(TOUT - T).LE.NROUND*UROUND*MAX(ABS(T), ABS(TOUT))) THEN
T = TOUT
NSTATE = 2
GO TO 560
END IF
C If there are no more roots
C to report, report T.
IF (NSTATE .EQ. 5) THEN
NSTATE = 2
GO TO 560
END IF
NSTATE = 2
C See if TOUT will
C be overtaken.
IF ((T + H)*HSIGN .GT. TOUT*HSIGN) THEN
H = TOUT - T
IF ((T + H)*HSIGN .GT. TOUT*HSIGN) H = H*(1.E0 - 4.E0*UROUND)
WORK(IH) = H
IF (H .EQ. 0.E0) GO TO 670
IWORK(IJTASK) = -1
END IF
ELSE IF (NTASK .EQ. 3) THEN
NSTATE = 2
IF (T*HSIGN .GT. TOUT*HSIGN) THEN
WRITE(MSG, '(''SDRIV32WRN With NTASK='', I1, '' on input, '',
8 ''T,'', E16.8, '', was beyond TOUT,'', E16.8, ''. Solution'',
8 '' obtained by interpolation.'')') NTASK, T, TOUT
CALL XERROR(MSG(1:124), 124, 2, 0)
CALL SDNTP (H, 0, N, IWORK(INQ), T, TOUT, WORK(IYH), Y)
T = TOUT
GO TO 580
END IF
IF (ABS(TOUT - T).LE.NROUND*UROUND*MAX(ABS(T), ABS(TOUT))) THEN
T = TOUT
GO TO 560
END IF
IF ((T + H)*HSIGN .GT. TOUT*HSIGN) THEN
H = TOUT - T
IF ((T + H)*HSIGN .GT. TOUT*HSIGN) H = H*(1.E0 - 4.E0*UROUND)
WORK(IH) = H
IF (H .EQ. 0.E0) GO TO 670
IWORK(IJTASK) = -1
END IF
END IF
C Implement changes in MINT, MITER, and/or HMAX.
C
IF ((MINT .NE. IWORK(IMNTLD) .OR. MITER .NE. IWORK(IMTRLD)) .AND.
8 MINT .NE. 3 .AND. IWORK(IMNTLD) .NE. 3) IWORK(IJTASK) = -1
IF (HMAX .NE. WORK(IHMAX)) THEN
H = SIGN(MIN(ABS(H), HMAX), H)
IF (H .NE. WORK(IH)) THEN
IWORK(IJTASK) = -1
WORK(IH) = H
END IF
WORK(IHMAX) = HMAX
END IF
C
180 NSTEPL = IWORK(INSTEP)
DO 190 I = 1,N
JYH = IYH + I - 1
190 Y(I) = WORK(JYH)
IF (NROOT .NE. 0) THEN
DO 200 I = 1,NROOT
JGNOW = IGNOW + I - 1
WORK(JGNOW) = G (NPAR, T, Y, I)
IF (NPAR .EQ. 0) THEN
IWORK(INROOT) = I
NSTATE = 7
RETURN
END IF
200 CONTINUE
END IF
IF (IERROR .EQ. 1) THEN
DO 230 I = 1,N
JYWT = I + IYWT - 1
230 WORK(JYWT) = 1.E0
GO TO 410
ELSE IF (IERROR .EQ. 5) THEN
DO 250 I = 1,N
JYWT = I + IYWT - 1
250 WORK(JYWT) = EWT(I)
GO TO 410
END IF
C Reset YWT array. Looping point.
260 IF (IERROR .EQ. 2) THEN
DO 280 I = 1,N
IF (Y(I) .EQ. 0.E0) GO TO 290
JYWT = I + IYWT - 1
280 WORK(JYWT) = ABS(Y(I))
GO TO 410
290 IF (IWORK(IJTASK) .EQ. 0) THEN
CALL F (NPAR, T, Y, WORK(ISAVE2))
IF (NPAR .EQ. 0) THEN
NSTATE = 6
RETURN
END IF
IWORK(INFE) = IWORK(INFE) + 1
IF (MITER .EQ. 3 .AND. IMPL .NE. 0) THEN
IFLAG = 0
CALL USERS(Y, WORK(IYH), WORK(IYWT), WORK(ISAVE1),
8 WORK(ISAVE2), T, H, WORK(IEL), IMPL, NPAR,
8 NDECOM, IFLAG)
IF (NPAR .EQ. 0) THEN
NSTATE = 10
RETURN
END IF
ELSE IF (IMPL .EQ. 1) THEN
IF (MITER .EQ. 1 .OR. MITER .EQ. 2) THEN
CALL FA (NPAR, T, Y, WORK(IA), MATDIM, ML, MU, NDECOM)
IF (NPAR .EQ. 0) THEN
NSTATE = 9
RETURN
END IF
CALL SGEFA (WORK(IA), MATDIM, N, IWORK(INDPVT), INFO)
IF (INFO .NE. 0) GO TO 690
CALL SGESL(WORK(IA),MATDIM,N,IWORK(INDPVT),WORK(ISAVE2),0)
ELSE IF (MITER .EQ. 4 .OR. MITER .EQ. 5) THEN
JAML = IA + ML
CALL FA (NPAR, T, Y, WORK(JAML), MATDIM, ML, MU, NDECOM)
IF (NPAR .EQ. 0) THEN
NSTATE = 9
RETURN
END IF
CALL SGBFA (WORK(IA),MATDIM,N,ML,MU,IWORK(INDPVT),INFO)
IF (INFO .NE. 0) GO TO 690
CALL SGBSL (WORK(IA), MATDIM, N, ML, MU, IWORK(INDPVT),
8 WORK(ISAVE2), 0)
END IF
ELSE IF (IMPL .EQ. 2) THEN
CALL FA (NPAR, T, Y, WORK(IA), MATDIM, ML, MU, NDECOM)
IF (NPAR .EQ. 0) THEN
NSTATE = 9
RETURN
END IF
DO 340 I = 1,NDECOM
JA = I + IA - 1
JSAVE2 = I + ISAVE2 - 1
IF (WORK(JA) .EQ. 0.E0) GO TO 690
340 WORK(JSAVE2) = WORK(JSAVE2)/WORK(JA)
END IF
END IF
DO 360 J = I,N
JYWT = J + IYWT - 1
IF (Y(J) .NE. 0.E0) THEN
WORK(JYWT) = ABS(Y(J))
ELSE
IF (IWORK(IJTASK) .EQ. 0) THEN
JSAVE2 = J + ISAVE2 - 1
WORK(JYWT) = ABS(H*WORK(JSAVE2))
ELSE
JHYP = J + IYH + N - 1
WORK(JYWT) = ABS(WORK(JHYP))
END IF
END IF
IF (WORK(JYWT) .EQ. 0.E0) WORK(JYWT) = UROUND
360 CONTINUE
ELSE IF (IERROR .EQ. 3) THEN
DO 380 I = 1,N
JYWT = I + IYWT - 1
380 WORK(JYWT) = MAX(EWT(1), ABS(Y(I)))
ELSE IF (IERROR .EQ. 4) THEN
DO 400 I = 1,N
JYWT = I + IYWT - 1
400 WORK(JYWT) = MAX(EWT(I), ABS(Y(I)))
END IF
C
410 DO 420 I = 1,N
JYWT = I + IYWT - 1
JSAVE2 = I + ISAVE2 - 1
420 WORK(JSAVE2) = Y(I)/WORK(JYWT)
SUM = SNRM2(N, WORK(ISAVE2), 1)/SQRT(REAL(N))
IF (EPS .LT. SUM*UROUND) THEN
EPS = SUM*UROUND*(1.E0 + 10.E0*UROUND)
WRITE(MSG, '(''SDRIV34REC At T,'', E16.8, '', the requested '',
8 ''accuracy, EPS, was not obtainable with the machine '',
8 ''precision. EPS has been increased to'')') T
WRITE(MSG(137:), '(E16.8)') EPS
CALL XERROR(MSG(1:152), 152, 4, 1)
NSTATE = 4
GO TO 560
END IF
IF (ABS(H) .GE. UROUND*ABS(T)) THEN
IWORK(INDPRT) = 0
ELSE IF (IWORK(INDPRT) .EQ. 0) THEN
WRITE(MSG, '(''SDRIV35WRN At T,'', E16.8, '', the step size,'',
8 E16.8, '', is smaller than the roundoff level of T. '')') T, H
WRITE(MSG(109:), '(''This may occur if there is an abrupt '',
8 ''change in the right hand side of the differential '',
8 ''equations.'')')
CALL XERROR(MSG(1:205), 205, 5, 0)
IWORK(INDPRT) = 1
END IF
IF (NTASK.NE.2) THEN
IF ((IWORK(INSTEP)-NSTEPL) .GT. MXSTEP) THEN
WRITE(MSG, '(''SDRIV33WRN At T,'', E16.8, '', '', I8,
8 '' steps have been taken without reaching TOUT,'', E16.8)')
8 T, MXSTEP, TOUT
CALL XERROR(MSG(1:103), 103, 3, 0)
NSTATE = 3
GO TO 560
END IF
END IF
C
C CALL SDSTP (EPS, F, FA, HMAX, IMPL, JACOBN, MATDIM, MAXORD,
C 8 MINT, MITER, ML, MU, N, NDE, YWT, UROUND, USERS,
C 8 AVGH, AVGORD, H, HUSED, JTASK, MNTOLD, MTROLD,
C 8 NFE, NJE, NQUSED, NSTEP, T, Y, YH, A, CONVRG,
C 8 DFDY, EL, FAC, HOLD, IPVT, JSTATE, NQ, NWAIT, RC,
C 8 RMAX, SAVE1, SAVE2, TQ, TREND, ISWFLG, MTRSV, MXRDSV)
C
CALL SDSTP (EPS, F, FA, WORK(IHMAX), IMPL, JACOBN, MATDIM,
8 IWORK(IMXORD), IWORK(IMNT), IWORK(IMTR), ML, MU, NPAR,
8 NDECOM, WORK(IYWT), UROUND, USERS, WORK(IAVGH),
8 WORK(IAVGRD), WORK(IH), WORK(IHUSED), IWORK(IJTASK),
8 IWORK(IMNTLD), IWORK(IMTRLD), IWORK(INFE), IWORK(INJE),
8 IWORK(INQUSD), IWORK(INSTEP), WORK(IT), Y, WORK(IYH),
8 WORK(IA), CONVRG, WORK(IDFDY), WORK(IEL), WORK(IFAC),
8 WORK(IHOLD), IWORK(INDPVT), JSTATE, IWORK(INQ),
8 IWORK(INWAIT), WORK(IRC), WORK(IRMAX), WORK(ISAVE1),
8 WORK(ISAVE2), WORK(ITQ), WORK(ITREND), MINT,
8 IWORK(IMTRSV), IWORK(IMXRDS))
T = WORK(IT)
H = WORK(IH)
GO TO (470, 670, 680, 690, 690, 660, 660, 660, 660, 660), JSTATE
470 IWORK(IJTASK) = 1
C Determine if a root has been overtaken
IF (NROOT .NE. 0) THEN
IROOT = 0
DO 500 I = 1,NROOT
JTROOT = ITROOT + I - 1
JGNOW = IGNOW + I - 1
GLAST = WORK(JGNOW)
WORK(JGNOW) = G (NPAR, T, Y, I)
IF (NPAR .EQ. 0) THEN
IWORK(INROOT) = I
NSTATE = 7
RETURN
END IF
IF (GLAST*WORK(JGNOW) .GT. 0.E0) THEN
WORK(JTROOT) = T + H
ELSE
IF (WORK(JGNOW) .EQ. 0.E0) THEN
WORK(JTROOT) = T
IROOT = I
ELSE
IF (GLAST .EQ. 0.E0) THEN
WORK(JTROOT) = T + H
ELSE
IF (ABS(WORK(IHUSED)) .GE. UROUND*ABS(T)) THEN
TLAST = T - WORK(IHUSED)
IROOT = I
TROOT = T
CALL SDZRO (AE, G, H, NPAR, IWORK(INQ), IROOT, RE, T,
8 WORK(IYH), UROUND, TROOT, TLAST,
8 WORK(JGNOW), GLAST, Y)
DO 480 J = 1,N
480 Y(J) = WORK(IYH + J -1)
IF (NPAR .EQ. 0) THEN
IWORK(INROOT) = I
NSTATE = 7
RETURN
END IF
WORK(JTROOT) = TROOT
ELSE
WORK(JTROOT) = T
IROOT = I
END IF
END IF
END IF
END IF
500 CONTINUE
IF (IROOT .EQ. 0) THEN
IWORK(IJROOT) = 0
C Select the first root
ELSE
IWORK(IJROOT) = NTASK
IWORK(INRTLD) = NROOT
IWORK(INDTRT) = ITROOT
TROOT = T + H
DO 510 I = 1,NROOT
JTROOT = ITROOT + I - 1
IF (WORK(JTROOT)*HSIGN .LT. TROOT*HSIGN) THEN
TROOT = WORK(JTROOT)
IROOT = I
END IF
510 CONTINUE
IWORK(INROOT) = IROOT
WORK(ITOUT) = TROOT
IF (TROOT*HSIGN .LE. TOUT*HSIGN) THEN
CALL SDNTP (H, 0, N, IWORK(INQ), T, TROOT, WORK(IYH), Y)
NSTATE = 5
T = TROOT
GO TO 580
END IF
END IF
END IF
C Test for NTASK condition to be satisfied
NSTATE = 2
IF (NTASK .EQ. 1) THEN
IF (T*HSIGN .LT. TOUT*HSIGN) GO TO 260
CALL SDNTP (H, 0, N, IWORK(INQ), T, TOUT, WORK(IYH), Y)
T = TOUT
GO TO 580
C TOUT is assumed to have been attained
C exactly if T is within twenty roundoff
C units of TOUT, relative to max(TOUT, T).
ELSE IF (NTASK .EQ. 2) THEN
IF (ABS(TOUT - T).LE.NROUND*UROUND*MAX(ABS(T), ABS(TOUT))) THEN
T = TOUT
ELSE
IF ((T + H)*HSIGN .GT. TOUT*HSIGN) THEN
H = TOUT - T
IF ((T + H)*HSIGN.GT.TOUT*HSIGN) H = H*(1.E0 - 4.E0*UROUND)
WORK(IH) = H
IF (H .EQ. 0.E0) GO TO 670
IWORK(IJTASK) = -1
END IF
END IF
ELSE IF (NTASK .EQ. 3) THEN
IF (ABS(TOUT - T).LE.NROUND*UROUND*MAX(ABS(T), ABS(TOUT))) THEN
T = TOUT
ELSE
IF ((T + H)*HSIGN .GT. TOUT*HSIGN) THEN
H = TOUT - T
IF ((T + H)*HSIGN.GT.TOUT*HSIGN) H = H*(1.E0 - 4.E0*UROUND)
WORK(IH) = H
IF (H .EQ. 0.E0) GO TO 670
IWORK(IJTASK) = -1
END IF
GO TO 260
END IF
END IF
C All returns are made through this
C section. IMXERR is determined.
560 DO 570 I = 1,N
JYH = I + IYH - 1
570 Y(I) = WORK(JYH)
580 IF (CONVRG) THEN
IWORK(ICNVRG) = 1
ELSE
IWORK(ICNVRG) = 0
END IF
IF (IWORK(IJTASK) .EQ. 0) RETURN
BIG = 0.E0
IMXERR = 1
IWORK(INDMXR) = IMXERR
DO 590 I = 1,N
C SIZE = ABS(ERROR(I)/YWT(I))
JYWT = I + IYWT - 1
JERROR = I + ISAVE1 - 1
SIZE = ABS(WORK(JERROR)/WORK(JYWT))
IF (BIG .LT. SIZE) THEN
BIG = SIZE
IMXERR = I
IWORK(INDMXR) = IMXERR
END IF
590 CONTINUE
RETURN
C
660 NSTATE = JSTATE
RETURN
C Fatal errors are processed here
C
670 WRITE(MSG, '(''SDRIV311FE At T,'', E16.8, '', the attempted '',
8 ''step size has gone to zero. Often this occurs if the '',
8 ''problem setup is incorrect.'')') T
CALL XERROR(MSG(1:129), 129, 11, 2)
RETURN
C
680 WRITE(MSG, '(''SDRIV312FE At T,'', E16.8, '', the step size has'',
8 '' been reduced about 50 times without advancing the '')') T
WRITE(MSG(103:), '(''solution. Often this occurs if the '',
8 ''problem setup is incorrect.'')')
CALL XERROR(MSG(1:165), 165, 12, 2)
RETURN
C
690 WRITE(MSG, '(''SDRIV313FE At T,'', E16.8, '', while solving'',
8 '' A*YDOT = F, A is singular.'')') T
CALL XERROR(MSG(1:74), 74, 13, 2)
RETURN
END
SUBROUTINE SGBFA(ABD,LDA,N,ML,MU,IPVT,INFO)
C***BEGIN PROLOGUE SGBFA
C THIS PROLOGUE HAS BEEN REMOVED FOR REASONS OF SPACE
C FOR A COMPLETE COPY OF THIS ROUTINE CONTACT THE AUTHORS
C
C From the book "Numerical Methods and Software"
C by D. Kahaner, C. Moler, S. Nash
C Prentice Hall 1988
C
C***END PROLOGUE SGBFA
INTEGER LDA,N,ML,MU,IPVT(*),INFO
REAL ABD(LDA,*)
C
REAL T
INTEGER I,ISAMAX,I0,J,JU,JZ,J0,J1,K,KP1,L,LM,M,MM,NM1
C
C***FIRST EXECUTABLE STATEMENT SGBFA
M = ML + MU + 1
INFO = 0
C
C ZERO INITIAL FILL-IN COLUMNS
C
J0 = MU + 2
J1 = MIN0(N,M) - 1
IF (J1 .LT. J0) GO TO 30
DO 20 JZ = J0, J1
I0 = M + 1 - JZ
DO 10 I = I0, ML
ABD(I,JZ) = 0.0E0
10 CONTINUE
20 CONTINUE
30 CONTINUE
JZ = J1
JU = 0
C
C GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING
C
NM1 = N - 1
IF (NM1 .LT. 1) GO TO 130
DO 120 K = 1, NM1
KP1 = K + 1
C
C ZERO NEXT FILL-IN COLUMN
C
JZ = JZ + 1
IF (JZ .GT. N) GO TO 50
IF (ML .LT. 1) GO TO 50
DO 40 I = 1, ML
ABD(I,JZ) = 0.0E0
40 CONTINUE
50 CONTINUE
C
C FIND L = PIVOT INDEX
C
LM = MIN0(ML,N-K)
L = ISAMAX(LM+1,ABD(M,K),1) + M - 1
IPVT(K) = L + K - M
C
C ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED
C
IF (ABD(L,K) .EQ. 0.0E0) GO TO 100
C
C INTERCHANGE IF NECESSARY
C
IF (L .EQ. M) GO TO 60
T = ABD(L,K)
ABD(L,K) = ABD(M,K)
ABD(M,K) = T
60 CONTINUE
C
C COMPUTE MULTIPLIERS
C
T = -1.0E0/ABD(M,K)
CALL SSCAL(LM,T,ABD(M+1,K),1)
C
C ROW ELIMINATION WITH COLUMN INDEXING
C
JU = MIN0(MAX0(JU,MU+IPVT(K)),N)
MM = M
IF (JU .LT. KP1) GO TO 90
DO 80 J = KP1, JU
L = L - 1
MM = MM - 1
T = ABD(L,J)
IF (L .EQ. MM) GO TO 70
ABD(L,J) = ABD(MM,J)
ABD(MM,J) = T
70 CONTINUE
CALL SAXPY(LM,T,ABD(M+1,K),1,ABD(MM+1,J),1)
80 CONTINUE
90 CONTINUE
GO TO 110
100 CONTINUE
INFO = K
110 CONTINUE
120 CONTINUE
130 CONTINUE
IPVT(N) = N
IF (ABD(M,N) .EQ. 0.0E0) INFO = N
RETURN
END
SUBROUTINE SGBSL(ABD,LDA,N,ML,MU,IPVT,B,JOB)
C***BEGIN PROLOGUE SGBSL
C THIS PROLOGUE HAS BEEN REMOVED FOR REASONS OF SPACE
C FOR A COMPLETE COPY OF THIS ROUTINE CONTACT THE AUTHORS
C
C From the book "Numerical Methods and Software"
C by D. Kahaner, C. Moler, S. Nash
C Prentice Hall 1988
C
C***END PROLOGUE SGBSL
INTEGER LDA,N,ML,MU,IPVT(*),JOB
REAL ABD(LDA,*),B(*)
C
REAL SDOT,T
INTEGER K,KB,L,LA,LB,LM,M,NM1
C***FIRST EXECUTABLE STATEMENT SGBSL
M = MU + ML + 1
NM1 = N - 1
IF (JOB .NE. 0) GO TO 50
C
C JOB = 0 , SOLVE A * X = B
C FIRST SOLVE L*Y = B
C
IF (ML .EQ. 0) GO TO 30
IF (NM1 .LT. 1) GO TO 30
DO 20 K = 1, NM1
LM = MIN0(ML,N-K)
L = IPVT(K)
T = B(L)
IF (L .EQ. K) GO TO 10
B(L) = B(K)
B(K) = T
10 CONTINUE
CALL SAXPY(LM,T,ABD(M+1,K),1,B(K+1),1)
20 CONTINUE
30 CONTINUE
C
C NOW SOLVE U*X = Y
C
DO 40 KB = 1, N
K = N + 1 - KB
B(K) = B(K)/ABD(M,K)
LM = MIN0(K,M) - 1
LA = M - LM
LB = K - LM
T = -B(K)
CALL SAXPY(LM,T,ABD(LA,K),1,B(LB),1)
40 CONTINUE
GO TO 100
50 CONTINUE
C
C JOB = NONZERO, SOLVE TRANS(A) * X = B
C FIRST SOLVE TRANS(U)*Y = B
C
DO 60 K = 1, N
LM = MIN0(K,M) - 1
LA = M - LM
LB = K - LM
T = SDOT(LM,ABD(LA,K),1,B(LB),1)
B(K) = (B(K) - T)/ABD(M,K)
60 CONTINUE
C
C NOW SOLVE TRANS(L)*X = Y
C
IF (ML .EQ. 0) GO TO 90
IF (NM1 .LT. 1) GO TO 90
DO 80 KB = 1, NM1
K = N - KB
LM = MIN0(ML,N-K)
B(K) = B(K) + SDOT(LM,ABD(M+1,K),1,B(K+1),1)
L = IPVT(K)
IF (L .EQ. K) GO TO 70
T = B(L)
B(L) = B(K)
B(K) = T
70 CONTINUE
80 CONTINUE
90 CONTINUE
100 CONTINUE
RETURN
END