! MODULE RKSUITE
! CONTAINS
SUBROUTINE SETUP(NEQ,TSTART,YSTART,TEND,TOL,THRES,METHOD,TASK,
& ERRASS,HSTART,WORK,LENWRK,MESAGE)
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
C If you are not familiar with the code SETUP and how it is used in
C conjunction with UT or CT to solve initial value problems, you should study
C the document file rksuite.doc carefully before attempting to use the code.
C The following "Brief Reminder" is intended only to remind you of the
C meaning, type, and size requirements of the arguments.
C
C The environmental parameters OUTCH, MCHEPS, and DWARF are used in the
C following description. To find out their values
C
C CALL ENVIRN(OUTCH,MCHEPS,DWARF)
C
C INPUT VARIABLES
C
C NEQ - INTEGER
C The number of differential equations in the system.
C Constraint: NEQ >= 1
C TSTART - DOUBLE PRECISION
C The initial value of the independent variable.
C YSTART(*) - DOUBLE PRECISION array of length NEQ
C The vector of initial values of the solution components.
C TEND - DOUBLE PRECISION
C The integration proceeds from TSTART in the direction of
C TEND. You cannot go past TEND.
C Constraint: TEND must be clearly distinguishable from TSTART
C in the precision available.
C TOL - DOUBLE PRECISION
C The relative error tolerance.
C Constraint: 0.01D0 >= TOL >= 10*MCHEPS
C THRES(*) - DOUBLE PRECISION array of length NEQ
C THRES(L) is the threshold for the Ith solution component.
C Constraint: THRES(L) >= SQRT(DWARF)
C METHOD - INTEGER
C Specifies which Runge-Kutta pair is to be used.
C = 1 - use the (2,3) pair
C = 2 - use the (4,5) pair
C = 3 - use the (7,8) pair
C TASK - CHARACTER*(*)
C Only the first character of TASK is significant.
C TASK(1:1) = `U' or `u' - UT is to be used
C = `C' or `c' - CT is to be used
C Constraint: TASK(1:1) = `U'or `u' or`C' or `c'
C ERRASS - LOGICAL
C = .FALSE. - do not attempt to assess the true error.
C = .TRUE. - assess the true error. Costs roughly twice
C as much as the integration with METHODs 2 and
C 3, and three times with METHOD = 1.
C HSTART - DOUBLE PRECISION
C 0.0D0 - select automatically the first step size.
C non-zero - try HSTART for the first step.
C
C WORKSPACE
C
C WORK(*) - DOUBLE PRECISION array of length LENWRK
C Do not alter the contents of this array after calling SETUP.
C
C INPUT VARIABLES
C
C LENWRK - INTEGER
C Length of WORK(*): How big LENWRK must be depends
C on the task and how it is to be solved.
C
C LENWRK = 32*NEQ is sufficient for all cases.
C
C If storage is a problem, the least storage possible
C in the various cases is:
C
C If TASK = `U' or `u', then
C if ERRASS = .FALSE. and
C METHOD = 1, LENWRK must be at least 10*NEQ
C = 2 20*NEQ
C = 3 16*NEQ
C if ERRASS = .TRUE. and
C METHOD = 1, LENWRK must be at least 15*NEQ
C = 2 32*NEQ
C = 3 21*NEQ
C
C If TASK = `C' or `c', then
C if ERRASS = .FALSE. and
C METHOD = 1, LENWRK must be at least 10*NEQ
C = 2 14*NEQ
C = 3 16*NEQ
C if ERRASS = .TRUE. and
C METHOD = 1, LENWRK must be at least 15*NEQ
C = 2 26*NEQ
C = 3 21*NEQ
C
C Warning: To exploit the interpolation capability
C of METHODs 1 and 2, you have to call INTRP. This
C subroutine requires working storage in addition to
C that specified here.
C
C MESAGE - LOGICAL
C Specifies whether you want informative messages written to
C the standard output channel OUTCH.
C = .TRUE. - provide messages
C = .FALSE. - do not provide messages
C
C In the event of a "catastrophic" failure to call SETUP correctly, the
C nature of the catastrophe is reported on the standard output channel,
C regardless of the value of MESAGE. Unless special provision was made
C in advance (see rksuite.doc), the computation then comes to a STOP.
C
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C .. Scalar Arguments ..
DOUBLE PRECISION HSTART, TEND, TOL, TSTART
INTEGER LENWRK, METHOD, NEQ
LOGICAL ERRASS, MESAGE
CHARACTER*(*) TASK
C .. Array Arguments ..
DOUBLE PRECISION THRES(*), WORK(*), YSTART(*)
C .. Common Block for Problem Definition ..
DOUBLE PRECISION TSTRT, TND, DIR, HSTRT, TOLR
INTEGER NEQN
COMMON /RKCOM1/ TSTRT, TND, DIR, HSTRT, TOLR, NEQN
SAVE /RKCOM1/
C .. Common Block to hold Problem Status ..
DOUBLE PRECISION T, H, TOLD, HOLD
INTEGER NFCN, SVNFCN, OKSTP, FLSTP
LOGICAL FIRST, LAST
COMMON /RKCOM2/ T, H, TOLD, HOLD, NFCN, SVNFCN, OKSTP, FLSTP,
& FIRST, LAST
SAVE /RKCOM2/
C .. Common Block for General Workspace Pointers ..
INTEGER PRTHRS, PRERST, PRWT, PRYOLD, PRSCR, PRY, PRYP,
& PRSTGS, PRINTP, LNINTP
COMMON /RKCOM3/ PRTHRS, PRERST, PRWT, PRYOLD, PRSCR, PRY, PRYP,
& PRSTGS, PRINTP, LNINTP
SAVE /RKCOM3/
C .. Common Block to hold Formula Definitions ..
DOUBLE PRECISION A(13,13), B(13), C(13), BHAT(13), R(11,6),
& E(7)
INTEGER PTR(13), NSTAGE, METHD, MINTP
LOGICAL INTP
COMMON /RKCOM4/ A, B, C, BHAT, R, E, PTR, NSTAGE, METHD,
& MINTP, INTP
SAVE /RKCOM4/
C .. Common Block to hold Formula Characterisitcs ..
DOUBLE PRECISION TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
& RS, RS1, RS2, RS3, RS4
INTEGER ORDER, LSTSTG, MAXTRY, NSEC
LOGICAL FSAL
COMMON /RKCOM5/ TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
& RS, RS1, RS2, RS3, RS4, ORDER, LSTSTG, MAXTRY,
& NSEC, FSAL
SAVE /RKCOM5/
C .. Common Block for Global Error Assessment ..
DOUBLE PRECISION MAXERR, LOCMAX
INTEGER GNFCN, PRZSTG, PRZY, PRZYP, PRZERS, PRZERR,
& PRZYNU
LOGICAL ERASON, ERASFL
COMMON /RKCOM6/ MAXERR, LOCMAX, GNFCN, PRZSTG, PRZY, PRZYP,
& PRZERS, PRZERR, PRZYNU, ERASON, ERASFL
SAVE /RKCOM6/
C .. Common Block for Environment Parameters ..
DOUBLE PRECISION MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY
INTEGER OUTCH
COMMON /RKCOM7/ MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY,
& OUTCH
SAVE /RKCOM7/
C .. Common Block for Integrator Options ..
LOGICAL MSG, UTASK
COMMON /RKCOM8/ MSG, UTASK
SAVE /RKCOM8/
C .. Common Block for Error Message ..
CHARACTER*80 REC(10)
COMMON /RKCOM9/ REC
SAVE /RKCOM9/
C .. Parameters ..
CHARACTER*6 SRNAME
PARAMETER (SRNAME='SETUP')
INTEGER MINUS1
LOGICAL TELL
PARAMETER (MINUS1=-1,TELL=.FALSE.)
DOUBLE PRECISION ONE, ZERO, PT01
PARAMETER (ONE=1.0D+0,ZERO=0.0D+0,PT01=0.01D+0)
C .. Local Scalars ..
DOUBLE PRECISION HMIN
INTEGER FLAG, FREEPR, IER, L, LINTPL, LREQ, NREC, VECSTG
LOGICAL LEGALT, REQSTG
CHARACTER TASK1
C .. External Subroutines ..
EXTERNAL CONST, MCONST, RKMSG, RKSIT
C .. Intrinsic Functions ..
INTRINSIC ABS, MAX, SIGN
C .. Executable Statements ..
C
C Clear previous flag values of subprograms in the suite.
C
IER = MINUS1
CALL RKSIT(TELL,SRNAME,IER)
C
IER = 1
NREC = 0
C
C Fetch output channel and machine constants; initialise common
C block /RKCOM7/
C
CALL MCONST(METHOD)
C
C Check for valid input of trivial arguments
TASK1 = TASK(1:1)
LEGALT = TASK1 .EQ. 'U' .OR. TASK1 .EQ. 'u' .OR.
& TASK1 .EQ. 'C' .OR. TASK1 .EQ. 'c'
IF (.NOT.LEGALT) THEN
IER = 911
NREC = 3
WRITE (REC,'(A/A,A,A/A)')
&' ** You have set the first character of ',
&' ** TASK to be ''',TASK1,'''. It must be one of ',
&' ** ''U'',''u'',''C'' or ''c''.'
ELSE IF (NEQ.LT.1) THEN
IER = 911
NREC = 1
WRITE (REC,'(A,I6,A)')
&' ** You have set NEQ = ',NEQ,' which is less than 1.'
ELSE IF (METHOD.LT.1 .OR. METHOD.GT.3) THEN
IER = 911
NREC = 1
WRITE (REC,'(A,I6,A)')
&' ** You have set METHOD = ',METHOD,' which is not 1, 2, or 3.'
ELSE IF (TSTART.EQ.TEND) THEN
IER = 911
NREC = 1
WRITE (REC,'(A,D13.5,A)')
&' ** You have set TSTART = TEND = ',TSTART,'.'
ELSE IF ((TOL.GT.PT01) .OR. (TOL.LT.RNDOFF)) THEN
IER = 911
NREC = 2
WRITE (REC,'(A,D13.5,A/A,D13.5,A)')
&' ** You have set TOL = ',TOL,' which is not permitted. The',
&' ** range of permitted values is (',RNDOFF,',0.01D0).'
ELSE
L = 1
20 CONTINUE
IF (THRES(L).LT.TINY) THEN
IER = 911
NREC = 2
WRITE (REC,'(A,I6,A,D13.5,A/A,D13.5,A)')
&' ** You have set THRES(',L,') to be ',THRES(L),' which is ',
&' ** less than the permitted minimum,',TINY,'.'
END IF
L = L + 1
IF (IER.NE.911 .AND. L.LE.NEQ) GO TO 20
END IF
C
C Return if error detected
C
IF (IER.NE.1) GO TO 80
C
C Set formula definitions and characteristics by means of arguments
C in the call list and COMMON blocks /RKCOM4/ and /RKCOM5/
C
CALL CONST(METHOD,VECSTG,REQSTG,LINTPL)
C
C Set options in /RKCOM8/
UTASK = TASK1 .EQ. 'U' .OR. TASK1 .EQ. 'u'
MSG = MESAGE
C
C Initialise problem status in /RKCOM1/ and /RKCOM2/
NEQN = NEQ
TSTRT = TSTART
TND = TEND
T = TSTART
TOLD = TSTART
DIR = SIGN(ONE,TEND-TSTART)
C
C In CT the first step taken will have magnitude H. If HSTRT = ABS(HSTART)
C is not equal to zero, H = HSTRT. If HSTRT is equal to zero, the code is
C to find an on-scale initial step size H. To start this process, H is set
C here to an upper bound on the first step size that reflects the scale of
C the independent variable. UT has some additional information, namely the
C first output point, that is used to refine this bound in UT when UTASK
C is .TRUE.. If HSTRT is not zero, but it is either too big or too small,
C the input HSTART is ignored and HSTRT is set to zero to activate the
C automatic determination of an on-scale initial step size.
C
HSTRT = ABS(HSTART)
HMIN = MAX(TINY,TOOSML*MAX(ABS(TSTART),ABS(TEND)))
IF (HSTRT.GT.ABS(TEND-TSTART) .OR. HSTRT.LT.HMIN) HSTRT = ZERO
IF (HSTRT.EQ.ZERO) THEN
H = MAX(ABS(TEND-TSTART)/RS3,HMIN)
ELSE
H = HSTRT
END IF
HOLD = ZERO
TOLR = TOL
NFCN = 0
SVNFCN = 0
OKSTP = 0
FLSTP = 0
FIRST = .TRUE.
LAST = .FALSE.
C
C WORK(*) is partioned into a number of arrays using pointers. These
C pointers are set in /RKCOM3/.
PRTHRS = 1
C the threshold values
PRERST = PRTHRS + NEQ
C the error estimates
PRWT = PRERST + NEQ
C the weights used in the local error test
PRYOLD = PRWT + NEQ
C the previous value of the solution
PRSCR = PRYOLD + NEQ
C scratch array used for the higher order
C approximate solution and for the previous
C value of the derivative of the solution
PRY = PRSCR + NEQ
C the dependent variables
PRYP = PRY + NEQ
C the derivatives
PRSTGS = PRYP + NEQ
C intermediate stages held in an internal
C array STAGES(NEQ,VECSTG)
C
FREEPR = PRSTGS + VECSTG*NEQ
C
C Allocate storage for interpolation if the TASK = `U' or `u' was
C specified. INTP and LINTPL returned by CONST indicate whether there
C is an interpolation scheme associated with the pair and how much
C storage is required.
C
PRINTP = 1
LNINTP = 1
IF (UTASK) THEN
IF (INTP) THEN
LNINTP = LINTPL*NEQ
IF (REQSTG) THEN
PRINTP = FREEPR
FREEPR = PRINTP + LNINTP
ELSE
PRINTP = PRSTGS
FREEPR = MAX(PRINTP+VECSTG*NEQ,PRINTP+LNINTP)
END IF
END IF
END IF
C
C Initialise state and allocate storage for global error assessment
C using /RKCOM6/
GNFCN = 0
MAXERR = ZERO
LOCMAX = TSTART
ERASON = ERRASS
ERASFL = .FALSE.
IF (ERRASS) THEN
C
C Storage is required for the stages of a secondary integration. The
C stages of the primary intergration can only be overwritten in the
C cases where there is no interpolant or the interpolant does not
C require information about the stages (e.g. METHOD 3 and METHOD 1,
C respectively).
IF (.NOT.REQSTG) THEN
PRZSTG = PRSTGS
ELSE
PRZSTG = FREEPR
FREEPR = PRZSTG + VECSTG*NEQ
END IF
PRZY = FREEPR
PRZYP = PRZY + NEQ
PRZERS = PRZYP + NEQ
PRZERR = PRZERS + NEQ
PRZYNU = PRZERR + NEQ
FREEPR = PRZYNU + NEQ
ELSE
PRZSTG = 1
PRZY = 1
PRZYP = 1
PRZERS = 1
PRZERR = 1
PRZYNU = 1
END IF
C
LREQ = FREEPR - 1
C
C Check for enough workspace and suitable range of integration
C
IF (LENWRK.LT.LREQ) THEN
IER = 911
NREC = 2
WRITE (REC,'(A/A,I6,A,I6,A)')
&' ** You have not supplied enough workspace. You gave LENWRK ',
&' ** as', LENWRK, ', but it must be at least ',LREQ,'.'
ELSE
HMIN = MAX(TINY,TOOSML*MAX(ABS(TSTART),ABS(TEND)))
IF (ABS(TEND-TSTART).LT.HMIN) THEN
IER = 911
NREC = 4
WRITE (REC,'(A/A/A,D13.5/A,D13.5,A)')
&' ** You have set values for TEND and TSTART that are not ',
&' ** clearly distinguishable for the method and the precision ',
&' ** of the computer being used. ABS(TEND-TSTART) is ',
&ABS(TEND-TSTART),
&' ** but should be at least ',HMIN,'.'
END IF
END IF
C
C Return if error detected
C
IF (IER.NE.1) GO TO 80
C
C Initialize elements of the workspace
DO 40 L = 1, NEQ
WORK(PRTHRS-1+L) = THRES(L)
WORK(PRY-1+L) = YSTART(L)
40 CONTINUE
C
C Initialize the global error to zero when ERRASS = .TRUE.
IF (ERRASS) THEN
DO 60 L = 1, NEQ
WORK(PRZERR-1+L) = ZERO
60 CONTINUE
END IF
C
80 CONTINUE
C
CALL RKMSG(IER,SRNAME,NREC,FLAG)
C
RETURN
END SUBROUTINE SETUP
C
SUBROUTINE UT(F,TWANT,TGOT,YGOT,YPGOT,YMAX,WORK,UFLAG)
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
C If you are not familiar with the code UT and how it is used in
C conjunction with SETUP to solve initial value problems, you should study
C the document file rksuite.doc carefully before proceeding further. The
C following "Brief Reminder" is intended only to remind you of the meaning,
C type, and size requirements of the arguments.
C
C NAME DECLARED IN AN EXTERNAL STATEMENT IN THE CALLING PROGRAM:
C
C F - name of the subroutine for evaluating the differential
C equations.
C
C The subroutine F must have the form
C
C SUBROUTINE F(T,Y,YP)
C DOUBLE PRECISION T,Y(*),YP(*)
C Given input values of the independent variable T and the solution
C components Y(*), for each L = 1,2,...,NEQ evaluate the differential
C equation for the derivative of the Ith solution component and place the
C value in YP(L). Do not alter the input values of T and Y(*).
C RETURN
C END
C
C INPUT VARIABLE
C
C TWANT - DOUBLE PRECISION
C The next value of the independent variable where a
C solution is desired.
C
C Constraints: TWANT must lie between the previous value
C of TGOT (TSTART on the first call) and TEND. TWANT can be
C equal to TEND, but it must be clearly distinguishable from
C the previous value of TGOT (TSTART on the first call) in
C the precision available.
C
C OUTPUT VARIABLES
C
C TGOT - DOUBLE PRECISION
C A solution has been computed at this value of the
C independent variable.
C YGOT(*) - DOUBLE PRECISION array of length NEQ
C Approximation to the true solution at TGOT. Do not alter
C the contents of this array
C YPGOT(*) - DOUBLE PRECISION array of length NEQ
C Approximation to the first derivative of the true
C solution at TGOT.
C YMAX(*) - DOUBLE PRECISION array of length NEQ
C YMAX(L) is the largest magnitude of YGOT(L) computed at any
C time in the integration from TSTART to TGOT. Do not alter
C the contents of this array.
C
C WORKSPACE
C
C WORK(*) - DOUBLE PRECISION array as used in SETUP
C Do not alter the contents of this array.
C
C OUTPUT VARIABLE
C
C UFLAG - INTEGER
C
C SUCCESS. TGOT = TWANT.
C = 1 - Complete success.
C
C "SOFT" FAILURES
C = 2 - Warning: You are using METHOD = 3 inefficiently
C by computing answers at many values of TWANT. If
C you really need answers at so many specific points,
C it would be more efficient to compute them with
C METHOD = 2. To do this you would need to restart
C from TGOT, YGOT(*) by a call to SETUP. If you wish
C to continue as you are, you may.
C = 3 - Warning: A considerable amount of work has been
C expended. If you wish to continue on to TWANT, just
C call UT again.
C = 4 - Warning: It appears that this problem is "stiff".
C You really should change to another code that is
C intended for such problems, but if you insist, you can
C continue with UT by calling it again.
C
C "HARD" FAILURES
C = 5 - You are asking for too much accuracy. You cannot
C continue integrating this problem.
C = 6 - The global error assessment may not be reliable beyond
C the current point in the integration. You cannot
C continue integrating this problem.
C
C "CATASTROPHIC" FAILURES
C = 911 - The nature of the catastrophe is reported on
C the standard output channel. Unless special
C provision was made in advance (see rksuite.doc),
C the computation then comes to a STOP.
C
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C .. Scalar Arguments ..
DOUBLE PRECISION TGOT, TWANT
INTEGER UFLAG
C .. Array Arguments ..
DOUBLE PRECISION WORK(*), YGOT(*), YMAX(*), YPGOT(*)
C .. Subroutine Arguments ..
EXTERNAL F
C .. Common Block for Problem Definition ..
DOUBLE PRECISION TSTRT, TND, DIR, HSTRT, TOLR
INTEGER NEQN
COMMON /RKCOM1/ TSTRT, TND, DIR, HSTRT, TOLR, NEQN
SAVE /RKCOM1/
C .. Common Block to hold Problem Status ..
DOUBLE PRECISION T, H, TOLD, HOLD
INTEGER NFCN, SVNFCN, OKSTP, FLSTP
LOGICAL FIRST, LAST
COMMON /RKCOM2/ T, H, TOLD, HOLD, NFCN, SVNFCN, OKSTP, FLSTP,
& FIRST, LAST
SAVE /RKCOM2/
C .. Common Block for General Workspace Pointers ..
INTEGER PRTHRS, PRERST, PRWT, PRYOLD, PRSCR, PRY, PRYP,
& PRSTGS, PRINTP, LNINTP
COMMON /RKCOM3/ PRTHRS, PRERST, PRWT, PRYOLD, PRSCR, PRY, PRYP,
& PRSTGS, PRINTP, LNINTP
SAVE /RKCOM3/
C .. Common Block to hold Formula Definitions ..
DOUBLE PRECISION A(13,13), B(13), C(13), BHAT(13), R(11,6),
& E(7)
INTEGER PTR(13), NSTAGE, METHD, MINTP
LOGICAL INTP
COMMON /RKCOM4/ A, B, C, BHAT, R, E, PTR, NSTAGE, METHD,
& MINTP, INTP
SAVE /RKCOM4/
C .. Common Block to hold Formula Characterisitcs ..
DOUBLE PRECISION TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
& RS, RS1, RS2, RS3, RS4
INTEGER ORDER, LSTSTG, MAXTRY, NSEC
LOGICAL FSAL
COMMON /RKCOM5/ TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
& RS, RS1, RS2, RS3, RS4, ORDER, LSTSTG, MAXTRY,
& NSEC, FSAL
SAVE /RKCOM5/
C .. Common Block for Environment Parameters ..
DOUBLE PRECISION MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY
INTEGER OUTCH
COMMON /RKCOM7/ MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY,
& OUTCH
SAVE /RKCOM7/
C .. Common Block for Integrator Options ..
LOGICAL MSG, UTASK
COMMON /RKCOM8/ MSG, UTASK
SAVE /RKCOM8/
C .. Common Block for Error Message ..
CHARACTER*80 REC(10)
COMMON /RKCOM9/ REC
SAVE /RKCOM9/
C .. Parameters ..
CHARACTER*6 SRNAME
PARAMETER (SRNAME='UT')
LOGICAL ASK, TELL
PARAMETER (ASK=.TRUE.,TELL=.FALSE.)
INTEGER MINUS1, MINUS2
PARAMETER (MINUS1=-1,MINUS2=-2)
DOUBLE PRECISION ZERO
PARAMETER (ZERO=0.0D+0)
C .. Local Scalars ..
DOUBLE PRECISION HMIN, TLAST, TNOW, UTEND
INTEGER CFLAG, IER, L, NREC, STATE
LOGICAL BADERR, GOBACK
C .. External Subroutines ..
EXTERNAL CHKFL, CT, INTRP, RESET, RKMSG, RKSIT
C .. Intrinsic Functions ..
INTRINSIC ABS, MAX, MIN
C .. Save statement ..
SAVE UTEND, TLAST
C .. Executable Statements ..
IER = 1
NREC = 0
GOBACK = .FALSE.
BADERR = .FALSE.
C
C Is it permissible to call UT?
C
CALL RKSIT(ASK,'SETUP',STATE)
IF (STATE.EQ.911) THEN
IER = 912
NREC = 1
WRITE (REC,'(A)')
&' ** A catastrophic error has already been detected elsewhere.'
GO TO 100
END IF
IF (STATE.EQ.MINUS1) THEN
IER = 911
NREC = 1
WRITE (REC,'(A)')
&' ** You have not called SETUP, so you cannot use UT.'
GO TO 100
END IF
IF (.NOT.UTASK) THEN
IER = 911
NREC = 2
WRITE (REC,'(A/A)')
&' ** You have called UT after you specified in SETUP that ',
&' ** you were going to use CT. This is not permitted.'
GO TO 100
END IF
CALL RKSIT(ASK,SRNAME,STATE)
IF (STATE.EQ.5 .OR. STATE.EQ.6) THEN
IER = 911
NREC = 1
WRITE (REC,'(A/A)')
&' ** This routine has already returned with a hard failure.',
&' ** You must call SETUP to start another problem.'
GO TO 100
END IF
STATE = MINUS2
CALL RKSIT(TELL,SRNAME,STATE)
C
IF (FIRST) THEN
C
C First call.
C
C A value of TND is specified in SETUP. When INTP = .FALSE., as with
C METHD = 3, output is obtained at the specified TWANT by resetting TND
C to TWANT. At this point, before the integration gets started, this can
C be done with a simple assignment. Later it is done with a call to RESET.
C The original TND is SAVEd as a local variable UTEND.
C
UTEND = TND
IF (.NOT.INTP) TND = TWANT
C
C The last TGOT returned is SAVEd in the variable TLAST. T (a variable
C passed through the common block RKCOM2) records how far the integration
C has advanced towards the specified TND. When output is obtained by
C interpolation, the integration goes past the TGOT returned (T is closer
C to the specified TND than TGOT). Initialize these variables and YMAX(*).
TLAST = TSTRT
TGOT = TSTRT
DO 20 L = 1, NEQN
YMAX(L) = ABS(WORK(PRY-1+L))
20 CONTINUE
C
C If the code is to find an on-scale initial step size H, a bound was placed
C on H in SETUP. Here the first output point is used to refine this bound.
IF (HSTRT.EQ.ZERO) THEN
H = MIN(ABS(H),ABS(TWANT-TSTRT))
HMIN = MAX(TINY,TOOSML*MAX(ABS(TSTRT),ABS(TND)))
H = MAX(H,HMIN)
END IF
C
ELSE
C
C Subsequent call.
C
IF (TLAST.EQ.UTEND) THEN
IER = 911
NREC = 3
WRITE (REC,'(A/A/A)')
&' ** You have called UT after reaching TEND. (Your last ',
&' ** call to UT resulted in TGOT = TEND.) To start a new ',
&' ** problem, you will need to call SETUP.'
GO TO 100
END IF
C
END IF
C
C Check for valid TWANT.
C
IF (DIR*(TWANT-TLAST).LE.ZERO) THEN
IER = 911
NREC = 4
WRITE (REC,'(A/A/A/A)')
&' ** You have made a call to UT with a TWANT that does ',
&' ** not lie between the previous value of TGOT (TSTART ',
&' ** on the first call) and TEND. This is not permitted. ',
&' ** Check your program carefully.'
GO TO 100
END IF
IF (DIR*(TWANT-UTEND).GT.ZERO) THEN
HMIN = MAX(TINY,TOOSML*MAX(ABS(TWANT),ABS(UTEND)))
IF (ABS(TWANT-UTEND).LT.HMIN) THEN
IER = 911
NREC = 5
WRITE (REC,'(A/A/A/A)')
&' ** You have made a call to UT with a TWANT that does ',
&' ** not lie between the previous value of TGOT (TSTART on ',
&' ** the first call) and TEND. This is not permitted. TWANT ',
&' ** is very close to TEND, so you may have meant to set ',
&' ** it to be TEND exactly. Check your program carefully. '
ELSE
IER = 911
NREC = 4
WRITE (REC,'(A/A/A/A)')
&' ** You have made a call to UT with a TWANT that does ',
&' ** not lie between the previous value of TGOT (TSTART ',
&' ** on the first call) and TEND. This is not permitted. ',
&' ** Check your program carefully.'
END IF
GO TO 100
END IF
IF (.NOT.INTP) THEN
HMIN = MAX(TINY,TOOSML*MAX(ABS(TLAST),ABS(TWANT)))
IF (ABS(TWANT-TLAST).LT.HMIN) THEN
IER = 911
NREC = 4
WRITE (REC,'(A/A/A/A,D13.5,A)')
&' ** You have made a call to UT with a TWANT that is not ',
&' ** sufficiently different from the last value of TGOT ',
&' ** (TSTART on the first call). When using METHOD = 3, ',
&' ** it must differ by at least ',HMIN,'.'
GO TO 100
END IF
C
C We have a valid TWANT. There is no interpolation with this METHD and
C therefore we step to TWANT exactly by resetting TND with a call to RESET.
C On the first step this matter is handled differently as explained above.
C
IF (.NOT.FIRST) THEN
CALL RESET(TWANT)
CALL CHKFL(ASK,BADERR)
IF (BADERR) GO TO 100
END IF
END IF
C
C Process output, decide whether to take another step.
C
40 CONTINUE
C
IF (INTP) THEN
C
C Interpolation is possible with this METHD. The integration has
C already reached T. If this is past TWANT, GOBACK is set .TRUE. and
C the answers are obtained by interpolation.
C
GOBACK = DIR*(T-TWANT) .GE. ZERO
IF (GOBACK) THEN
CALL INTRP(TWANT,'Both solution and derivative',NEQN,YGOT,
& YPGOT,F,WORK,WORK(PRINTP),LNINTP)
CALL CHKFL(ASK,BADERR)
IF (BADERR) GO TO 100
TGOT = TWANT
END IF
ELSE
C
C Interpolation is not possible with this METHD, so output is obtained
C by integrating to TWANT = TND. Both YGOT(*) and YPGOT(*) are then
C already loaded with the solution at TWANT by CT.
C
GOBACK = T .EQ. TWANT
IF (GOBACK) TGOT = TWANT
END IF
C
C Updating of YMAX(*) is done here to account for the fact that when
C interpolation is done, the integration goes past TGOT. Note that YGOT(*)
C is not defined until CT is called. YMAX(*) was initialized at TSTRT
C from values stored in WORK(*), so only needs to be updated for T
C different from TSTRT.
IF (T.NE.TSTRT) THEN
DO 60 L = 1, NEQN
YMAX(L) = MAX(YMAX(L),ABS(YGOT(L)))
60 CONTINUE
END IF
C
C If done, go to the exit point.
IF (GOBACK) GO TO 100
C
C Take a step with CT in the direction of TND. On exit, the solution is
C advanced to TNOW. The way CT is written, the approximate solution at
C TNOW is available in both YGOT(*) and in WORK(*). If output is obtained by
C stepping to the end (TNOW = TWANT = TND), YGOT(*) can be returned directly.
C If output is obtained by interpolation, the subroutine INTRP that does this
C uses the values in WORK(*) for its computations and places the approximate
C solution at TWANT in the array YGOT(*) for return to the calling program.
C The approximate derivative is handled in the same way. TNOW is output from
C CT and is actually a copy of T declared above in a common block.
C
CALL CT(F,TNOW,YGOT,YPGOT,WORK,CFLAG)
IER = CFLAG
C
C A successful step by CT is indicated by CFLAG = 1 or = 2.
IF (CFLAG.EQ.1) THEN
GO TO 40
ELSE IF (CFLAG.EQ.2) THEN
C
C Supplement the warning message written in CT.
NREC = 3
WRITE (REC,'(A/A/A)')
&' ** The last message was produced on a call to CT from UT. ',
&' ** In UT the appropriate action is to change to METHOD = 2,',
&' ** or, if insufficient memory is available, to METHOD = 1. '
ELSE IF (CFLAG.LE.6) THEN
NREC = 1
WRITE (REC,'(A)')
&' ** The last message was produced on a call to CT from UT.'
ELSE
BADERR = .TRUE.
END IF
TGOT = T
C
C Update YMAX(*) before the return.
DO 80 L = 1, NEQN
YMAX(L) = MAX(YMAX(L),ABS(YGOT(L)))
80 CONTINUE
C
C Exit point for UT.
C
100 CONTINUE
C
IF (BADERR) THEN
IER = 911
NREC = 4
WRITE (REC,'(A/A/A/A)')
&' ** An internal call by UT to a subroutine resulted in an ',
&' ** error that should not happen. Check your program ',
&' ** carefully for array sizes, correct number of arguments,',
&' ** type mismatches ... .'
END IF
C
TLAST = TGOT
C
C All exits are done here after a call to RKMSG to report
C what happened and set UFLAG.
C
CALL RKMSG(IER,SRNAME,NREC,UFLAG)
C
RETURN
END SUBROUTINE UT
C
SUBROUTINE STAT(TOTFCN,STPCST,WASTE,STPSOK,HNEXT)
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
C If you are not familiar with the code STAT and how it is used in
C conjunction with the integrators CT and UT, you should study the
C document file rksuite.doc carefully before attempting to use the code.
C The following "Brief Reminder" is intended only to remind you of the
C meaning, type, and size requirements of the arguments.
C
C STAT is called to obtain some details about the integration.
C
C OUTPUT VARIABLES
C
C TOTFCN - INTEGER
C Total number of calls to F in the integration so far --
C a measure of the cost of the integration.
C STPCST - INTEGER
C Cost of a typical step with this METHOD measured in
C calls to F.
C WASTE - DOUBLE PRECISION
C The number of attempted steps that failed to meet the
C local error requirement divided by the total number of
C steps attempted so far in the integration.
C STPSOK - INTEGER
C The number of accepted steps.
C HNEXT - DOUBLE PRECISION
C The step size the integrator plans to use for the next step.
C
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C .. Scalar Arguments ..
DOUBLE PRECISION HNEXT, WASTE
INTEGER STPCST, STPSOK, TOTFCN
C .. Common Block to hold Problem Status ..
DOUBLE PRECISION T, H, TOLD, HOLD
INTEGER NFCN, SVNFCN, OKSTP, FLSTP
LOGICAL FIRST, LAST
COMMON /RKCOM2/ T, H, TOLD, HOLD, NFCN, SVNFCN, OKSTP, FLSTP,
& FIRST, LAST
SAVE /RKCOM2/
C .. Common Block to hold Formula Characterisitcs ..
DOUBLE PRECISION TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
& RS, RS1, RS2, RS3, RS4
INTEGER ORDER, LSTSTG, MAXTRY, NSEC
LOGICAL FSAL
COMMON /RKCOM5/ TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
& RS, RS1, RS2, RS3, RS4, ORDER, LSTSTG, MAXTRY,
& NSEC, FSAL
SAVE /RKCOM5/
C .. Common Block for Global Error Assessment ..
DOUBLE PRECISION MAXERR, LOCMAX
INTEGER GNFCN, PRZSTG, PRZY, PRZYP, PRZERS, PRZERR,
& PRZYNU
LOGICAL ERASON, ERASFL
COMMON /RKCOM6/ MAXERR, LOCMAX, GNFCN, PRZSTG, PRZY, PRZYP,
& PRZERS, PRZERR, PRZYNU, ERASON, ERASFL
SAVE /RKCOM6/
C .. Common Block for Integrator Options ..
LOGICAL MSG, UTASK
COMMON /RKCOM8/ MSG, UTASK
SAVE /RKCOM8/
C .. Common Block for Error Message ..
CHARACTER*80 REC(10)
COMMON /RKCOM9/ REC
SAVE /RKCOM9/
C .. Parameters ..
CHARACTER*6 SRNAME
PARAMETER (SRNAME='STAT')
LOGICAL ASK
INTEGER MINUS1, MINUS2
PARAMETER (ASK=.TRUE.,MINUS1=-1,MINUS2=-2)
DOUBLE PRECISION ZERO
PARAMETER (ZERO=0.0D+0)
C .. Local Scalars ..
INTEGER FLAG, IER, NREC, STATE
C .. External Subroutines ..
EXTERNAL RKMSG, RKSIT
C .. Intrinsic Functions ..
INTRINSIC DBLE
C .. Executable Statements ..
C
IER = 1
NREC = 0
C
C Is it permissible to call STAT?
C
CALL RKSIT(ASK,SRNAME,STATE)
IF (STATE.EQ.911) THEN
IER = 912
NREC = 1
WRITE (REC,'(A)')
&' ** A catastrophic error has already been detected elsewhere.'
GO TO 20
END IF
IF (STATE.EQ.MINUS2) THEN
IER = 911
NREC = 3
WRITE (REC,'(A/A/A)')
&' ** You have already made a call to STAT after a hard ',
&' ** failure was reported from the integrator. You cannot',
&' ** call STAT again.'
GO TO 20
END IF
CALL RKSIT(ASK,'CT',STATE)
IF (STATE.EQ.MINUS1) THEN
IER = 911
NREC = 1
IF (UTASK) THEN
WRITE (REC,'(A)')
&' ** You have not called UT, so you cannot use STAT.'
ELSE
WRITE (REC,'(A)')
&' ** You have not called CT, so you cannot use STAT.'
END IF
GO TO 20
END IF
C
C Set flag so that the routine can only be called once after a hard
C failure from the integrator.
IF (STATE.EQ.5 .OR. STATE.EQ.6) IER = MINUS2
C
TOTFCN = SVNFCN + NFCN
IF (ERASON) TOTFCN = TOTFCN + GNFCN
STPCST = COST
STPSOK = OKSTP
IF (OKSTP.LE.1) THEN
WASTE = ZERO
ELSE
WASTE = DBLE(FLSTP)/DBLE(FLSTP+OKSTP)
END IF
HNEXT = H
C
20 CONTINUE
C
CALL RKMSG(IER,SRNAME,NREC,FLAG)
C
RETURN
END SUBROUTINE STAT
C
SUBROUTINE GLBERR(RMSERR,ERRMAX,TERRMX,WORK)
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
C If you are not familiar with the code GLBERR and how it is used in
C conjunction with UT and CT to solve initial value problems, you should
C study the document file rksuite.doc carefully before attempting to use
C the code. The following "Brief Reminder" is intended only to remind you
C of the meaning, type, and size requirements of the arguments.
C
C If ERRASS was set .TRUE. in the call to SETUP, then after any call to UT
C or CT to advance the integration to TNOW or TWANT, the subroutine GLBERR
C may be called to obtain an assessment of the true error of the integration.
C At each step and for each solution component Y(L), a more accurate "true"
C solution YT(L), an average magnitude "size(L)" of its size, and its error
C abs(Y(L) - YT(L))/max("size(L)",THRES(L))
C are computed. The assessment returned in RMSERR(L) is the RMS (root-mean-
C square) average of the error in the Lth solution component over all steps
C of the integration from TSTART through TNOW.
C
C OUTPUT VARIABLES
C
C RMSERR(*) - DOUBLE PRECISION array of length NEQ
C RMSERR(L) approximates the RMS average of the true error
C of the numerical solution for the Ith solution component,
C L = 1,2,...,NEQ. The average is taken over all steps from
C TSTART to TNOW.
C ERRMAX - DOUBLE PRECISION
C The maximum (approximate) true error taken over all
C solution components and all steps from TSTART to TNOW.
C TERRMX - DOUBLE PRECISION
C First value of the independent variable where the
C (approximate) true error attains the maximum value ERRMAX.
C
C WORKSPACE
C
C WORK(*) - DOUBLE PRECISION array as used in SETUP and UT or CT
C Do not alter the contents of this array.
C
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C .. Scalar Arguments ..
DOUBLE PRECISION ERRMAX, TERRMX
C .. Array Arguments ..
DOUBLE PRECISION RMSERR(*), WORK(*)
C .. Common Block for Problem Definition ..
DOUBLE PRECISION TSTRT, TND, DIR, HSTRT, TOLR
INTEGER NEQN
COMMON /RKCOM1/ TSTRT, TND, DIR, HSTRT, TOLR, NEQN
SAVE /RKCOM1/
C .. Common Block to hold Problem Status ..
DOUBLE PRECISION T, H, TOLD, HOLD
INTEGER NFCN, SVNFCN, OKSTP, FLSTP
LOGICAL FIRST, LAST
COMMON /RKCOM2/ T, H, TOLD, HOLD, NFCN, SVNFCN, OKSTP, FLSTP,
& FIRST, LAST
SAVE /RKCOM2/
C .. Common Block for Global Error Assessment ..
DOUBLE PRECISION MAXERR, LOCMAX
INTEGER GNFCN, PRZSTG, PRZY, PRZYP, PRZERS, PRZERR,
& PRZYNU
LOGICAL ERASON, ERASFL
COMMON /RKCOM6/ MAXERR, LOCMAX, GNFCN, PRZSTG, PRZY, PRZYP,
& PRZERS, PRZERR, PRZYNU, ERASON, ERASFL
SAVE /RKCOM6/
C .. Common Block for Integrator Options ..
LOGICAL MSG, UTASK
COMMON /RKCOM8/ MSG, UTASK
SAVE /RKCOM8/
C .. Common Block for Error Message ..
CHARACTER*80 REC(10)
COMMON /RKCOM9/ REC
SAVE /RKCOM9/
C .. Parameters ..
CHARACTER*6 SRNAME
PARAMETER (SRNAME='GLBERR')
LOGICAL ASK
PARAMETER (ASK=.TRUE.)
INTEGER MINUS1, MINUS2
PARAMETER (MINUS1=-1,MINUS2=-2)
C .. Local Scalars ..
INTEGER FLAG, IER, L, NREC, STATE
C .. External Subroutines ..
EXTERNAL RKMSG, RKSIT
C .. Intrinsic Functions ..
INTRINSIC SQRT
C .. Executable Statements ..
C
IER = 1
NREC = 0
C
C Is it permissible to call GLBERR?
C
CALL RKSIT(ASK,SRNAME,STATE)
IF (STATE.EQ.911) THEN
IER = 912
NREC = 1
WRITE (REC,'(A)')
&' ** A catastrophic error has already been detected elsewhere.'
GO TO 40
END IF
IF (STATE.EQ.MINUS2) THEN
IER = 911
NREC = 3
WRITE (REC,'(A/A/A)')
&' ** You have already made a call to GLBERR after a hard ',
&' ** failure was reported from the integrator. You cannot',
&' ** call GLBERR again.'
GO TO 40
END IF
CALL RKSIT(ASK,'CT',STATE)
IF (STATE.EQ.MINUS1) THEN
IER = 911
NREC = 1
IF (UTASK) THEN
WRITE (REC,'(A)')
&' ** You have not yet called UT, so you cannot call GLBERR.'
ELSE
WRITE (REC,'(A)')
&' ** You have not yet called CT, so you cannot call GLBERR.'
END IF
GO TO 40
END IF
C
C Set flag so that the routine can only be called once after a hard
C failure from the integrator.
IF (STATE.EQ.5 .OR. STATE.EQ.6) IER = MINUS2
C
C Check that ERRASS was set properly for error assessment in SETUP.
C
IF (.NOT.ERASON) THEN
IER = 911
NREC = 3
WRITE (REC,'(A/A/A)')
&' ** No error assessment is available since you did not ',
&' ** ask for it in your call to the routine SETUP.',
&' ** Check your program carefully.'
GO TO 40
END IF
C
C Check to see if the integrator has not actually taken a step.
C
IF (OKSTP.EQ.0) THEN
IER = 911
NREC = 3
WRITE (REC,'(A/A/A)')
&' ** The integrator has not actually taken any successful ',
&' ** steps. You cannot call GLBERR in this circumstance. ',
&' ** Check your program carefully.'
GO TO 40
END IF
C
C Compute RMS error and set output variables.
C
ERRMAX = MAXERR
TERRMX = LOCMAX
DO 20 L = 1, NEQN
RMSERR(L) = SQRT(WORK(PRZERR-1+L)/OKSTP)
20 CONTINUE
C
40 CONTINUE
C
CALL RKMSG(IER,SRNAME,NREC,FLAG)
C
RETURN
END SUBROUTINE GLBERR
C
SUBROUTINE CT(F,TNOW,YNOW,YPNOW,WORK,CFLAG)
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
C If you are not familiar with the code CT and how it is used in
C conjunction with SETUP to solve initial value problems, you should study
C the document file rksuite.doc carefully before attempting to use the code.
C The following "Brief Reminder" is intended only to remind you of the
C meaning, type, and size requirements of the arguments.
C
C NAME DECLARED IN AN EXTERNAL STATEMENT IN THE CALLING PROGRAM:
C
C F - name of the subroutine for evaluating the differential
C equations.
C
C The subroutine F must have the form
C
C SUBROUTINE F(T,Y,YP)
C DOUBLE PRECISION T,Y(*),YP(*)
C Using the input values of the independent variable T and the solution
C components Y(*), for each L = 1,2,...,NEQ evaluate the differential
C equation for the derivative of the Lth solution component and place the
C value in YP(L). Do not alter the input values of T and Y(*).
C RETURN
C END
C
C OUTPUT VARIABLES
C
C TNOW - DOUBLE PRECISION
C Current value of the independent variable.
C YNOW(*) - DOUBLE PRECISION array of length NEQ
C Approximation to the true solution at TNOW.
C YPNOW(*) - DOUBLE PRECISION array of length NEQ
C Approximation to the first derivative of the
C true solution at TNOW.
C
C WORKSPACE
C
C WORK(*) - DOUBLE PRECISION array as used in SETUP
C Do not alter the contents of this array.
C
C OUTPUT VARIABLE
C
C CFLAG - INTEGER
C
C SUCCESS. A STEP WAS TAKEN TO TNOW.
C = 1 - Complete success.
C
C "SOFT" FAILURES
C = 2 - Warning: You have obtained an answer by integrating
C to TEND (TNOW = TEND). You have done this at least
C 100 times, and monitoring of the computation reveals
C that this way of getting output has degraded the
C efficiency of the code. If you really need answers at
C so many specific points, it would be more efficient to
C get them with INTRP. (If METHOD = 3, you would need
C to change METHOD and restart from TNOW, YNOW(*) by a
C call to SETUP.) If you wish to continue as you are,
C you may.
C = 3 - Warning: A considerable amount of work has been
C expended. To continue the integration, just call
C CT again.
C = 4 - Warning: It appears that this problem is "stiff".
C You really should change to another code that is
C intended for such problems, but if you insist, you
C can continue with CT by calling it again.
C
C "HARD" FAILURES
C = 5 - You are asking for too much accuracy. You cannot
C continue integrating this problem.
C = 6 - The global error assessment may not be reliable beyond
C the current point in the integration. You cannot
C continue integrating this problem.
C
C "CATASTROPHIC" FAILURES
C = 911 - The nature of the catastrophe is reported on
C the standard output channel. Unless special
C provision was made in advance (see rksuite.doc),
C the computation then comes to a STOP.
C
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C .. Scalar Arguments ..
DOUBLE PRECISION TNOW
INTEGER CFLAG
C .. Array Arguments ..
DOUBLE PRECISION WORK(*), YNOW(*), YPNOW(*)
C .. Subroutine Arguments ..
EXTERNAL F
C .. Common Block for Problem Definition ..
DOUBLE PRECISION TSTRT, TND, DIR, HSTRT, TOLR
INTEGER NEQN
COMMON /RKCOM1/ TSTRT, TND, DIR, HSTRT, TOLR, NEQN
SAVE /RKCOM1/
C .. Common Block to hold Problem Status ..
DOUBLE PRECISION T, H, TOLD, HOLD
INTEGER NFCN, SVNFCN, OKSTP, FLSTP
LOGICAL FIRST, LAST
COMMON /RKCOM2/ T, H, TOLD, HOLD, NFCN, SVNFCN, OKSTP, FLSTP,
& FIRST, LAST
SAVE /RKCOM2/
C .. Common Block for General Workspace Pointers ..
INTEGER PRTHRS, PRERST, PRWT, PRYOLD, PRSCR, PRY, PRYP,
& PRSTGS, PRINTP, LNINTP
COMMON /RKCOM3/ PRTHRS, PRERST, PRWT, PRYOLD, PRSCR, PRY, PRYP,
& PRSTGS, PRINTP, LNINTP
SAVE /RKCOM3/
C .. Common Block to hold Formula Characterisitcs ..
DOUBLE PRECISION TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
& RS, RS1, RS2, RS3, RS4
INTEGER ORDER, LSTSTG, MAXTRY, NSEC
LOGICAL FSAL
COMMON /RKCOM5/ TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
& RS, RS1, RS2, RS3, RS4, ORDER, LSTSTG, MAXTRY,
& NSEC, FSAL
SAVE /RKCOM5/
C .. Common Block for Global Error Assessment ..
DOUBLE PRECISION MAXERR, LOCMAX
INTEGER GNFCN, PRZSTG, PRZY, PRZYP, PRZERS, PRZERR,
& PRZYNU
LOGICAL ERASON, ERASFL
COMMON /RKCOM6/ MAXERR, LOCMAX, GNFCN, PRZSTG, PRZY, PRZYP,
& PRZERS, PRZERR, PRZYNU, ERASON, ERASFL
SAVE /RKCOM6/
C .. Common Block for Environment Parameters ..
DOUBLE PRECISION MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY
INTEGER OUTCH
COMMON /RKCOM7/ MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY,
& OUTCH
SAVE /RKCOM7/
C .. Common Block for Integrator Options ..
LOGICAL MSG, UTASK
COMMON /RKCOM8/ MSG, UTASK
SAVE /RKCOM8/
C .. Common Block for Error Message ..
CHARACTER*80 REC(10)
COMMON /RKCOM9/ REC
SAVE /RKCOM9/
C .. Parameters ..
CHARACTER*6 SRNAME
PARAMETER (SRNAME='CT')
LOGICAL ASK, TELL
PARAMETER (ASK=.TRUE.,TELL=.FALSE.)
INTEGER MINUS1, MINUS2
PARAMETER (MINUS1=-1,MINUS2=-2)
INTEGER MAXFCN
PARAMETER (MAXFCN=5000)
DOUBLE PRECISION ZERO, PT1, PT9, ONE, TWO, HUNDRD
PARAMETER (ZERO=0.0D+0,PT1=0.1D+0,PT9=0.9D+0,ONE=1.0D+0,
& TWO=2.0D+0,HUNDRD=100.0D+0)
C .. Local Scalars ..
DOUBLE PRECISION ALPHA, BETA, ERR, ERROLD, HAVG, HMIN, HTRY, TAU,
& TEMP1, TEMP2, YPNORM
INTEGER IER, JFLSTP, L, NREC, NTEND, POINT, STATE, YNEW,
& YPOLD
LOGICAL CHKEFF, FAILED, MAIN, PHASE1, PHASE2, PHASE3,
& TOOMCH
C .. External Subroutines ..
EXTERNAL RKMSG, RKSIT, STEP, STIFF, TRUERR
C .. Intrinsic Functions ..
INTRINSIC ABS, MAX, MIN, SIGN
C .. Save statement ..
SAVE JFLSTP, NTEND, ERROLD, HAVG, PHASE2, YNEW,
& YPOLD, CHKEFF
C .. Executable Statements ..
C
IER = 1
NREC = 0
C
C Is it permissible to call CT?
C
CALL RKSIT(ASK,'SETUP',STATE)
IF (STATE.EQ.911) THEN
IER = 912
NREC = 1
WRITE (REC,'(A)')
&' ** A catastrophic error has already been detected elsewhere.'
GO TO 180
END IF
IF (STATE.EQ.MINUS1) THEN
IER = 911
NREC = 1
WRITE (REC,'(A)')
&' ** You have not called SETUP, so you cannot use CT.'
GO TO 180
END IF
IF (UTASK) THEN
CALL RKSIT(ASK,'UT',STATE)
IF (STATE.NE.MINUS2) THEN
IER = 911
NREC = 2
WRITE (REC,'(A/A)')
&' ** You have called CT after you specified in SETUP that ',
&' ** you were going to use UT. This is not permitted.'
UTASK = .FALSE.
GO TO 180
END IF
END IF
CALL RKSIT(ASK,SRNAME,STATE)
IF (STATE.EQ.5 .OR. STATE.EQ.6) THEN
IER = 911
NREC = 3
WRITE (REC,'(A/A)')
&' ** CT has already returned with a flag value of 5 or 6.',
&' ** You cannot continue integrating this problem. You must ',
&' ** call SETUP to start another problem.'
GO TO 180
END IF
C
IF (FIRST) THEN
C
C First call in an integration -- initialize everything.
C
CHKEFF = .FALSE.
NTEND = 0
JFLSTP = 0
C
C A scratch area of WORK(*) starting at PRSCR is used to hold two
C arrays in this subroutine: the higher order approximate solution at
C the end of a step and the approximate derivative of the solution at
C the end of the last step. To make this clear, local pointers YNEW and
C YPOLD are used.
YNEW = PRSCR
YPOLD = PRSCR
C
C For this first step T was initialized to TSTRT in SETUP and the
C starting values YSTART(*) were loaded into the area of WORK(*) reserved
C for the current solution approximation starting at location PRY. The
C derivative is now computed and stored in WORK(*) starting at PRYP.
C Subsequently these arrays are copied to the output vectors YNOW(*)
C and YPNOW(*).
CALL F(T,WORK(PRY),WORK(PRYP))
NFCN = NFCN + 1
DO 20 L = 1, NEQN
YNOW(L) = WORK(PRY-1+L)
YPNOW(L) = WORK(PRYP-1+L)
20 CONTINUE
C
C Set dependent variables for error assessment.
IF (ERASON) THEN
DO 40 L = 1, NEQN
WORK(PRZY-1+L) = YNOW(L)
WORK(PRZYP-1+L) = YPNOW(L)
40 CONTINUE
END IF
C
C The weights for the control of the error depend on the size of the
C solution at the beginning and at the end of the step. On the first
C step we do not have all this information. Whilst determining the
C initial step size we initialize the weight vector to the larger of
C abs(Y(L)) and the threshold for this component.
DO 60 L = 1, NEQN
WORK(PRWT-1+L) = MAX(ABS(YNOW(L)),WORK(PRTHRS-1+L))
60 CONTINUE
C
C If HSTRT is equal to zero, the code is to find an on-scale initial step
C size H. CT has an elaborate scheme of three phases for finding such an H,
C and some preparations are made earlier. In SETUP an upper bound is placed
C on H that reflects the scale of the independent variable. When UTASK is
C .TRUE., UT refines this bound using the first output point. Here in CT
C PHASE1 applies a rule of thumb based on the error control, the order of the
C the formula, and the size of the initial slope to get a crude approximation
C to an on-scale H. PHASE2 may reduce H in the course of taking the first
C step. PHASE3 repeatedly adjusts H and retakes the first step until H is
C on scale.
C
C A guess for the magnitude of the first step size H can be provided to SETUP
C as HSTART. If it is too big or too small, it is ignored and the automatic
C determination of an on-scale initial step size is activated. If it is
C acceptable, H is set to HSTART in SETUP. Even when H is supplied to CT,
C PHASE3 of the scheme for finding an on-scale initial step size is made
C active so that the code can deal with a bad guess.
C
PHASE1 = HSTRT .EQ. ZERO
PHASE2 = PHASE1
PHASE3 = .TRUE.
IF (PHASE1) THEN
H = ABS(H)
YPNORM = ZERO
DO 80 L = 1, NEQN
IF (ABS(YNOW(L)).NE.ZERO) THEN
YPNORM = MAX(YPNORM,ABS(YPNOW(L))/WORK(PRWT-1+L))
END IF
80 CONTINUE
TAU = TOLR**EXPON
IF (H*YPNORM.GT.TAU) H = TAU/YPNORM
HMIN = MAX(TINY,TOOSML*MAX(ABS(TSTRT),ABS(TND)))
H = DIR*MAX(H,HMIN)
PHASE1 = .FALSE.
END IF
C
ELSE
C
C Continuation call
C
IF (LAST) THEN
IER = 911
NREC = 3
WRITE (REC,'(A,D13.5,A/A/A)')
&' ** You have already reached TEND ( = ', TND, ').',
&' ** To integrate further with the same problem you must ',
&' ** call the routine RESET with a new value of TEND.'
GO TO 180
END IF
END IF
C
C Begin computation of a step here.
C
FAILED = .FALSE.
C
100 CONTINUE
H = SIGN(ABS(H),DIR)
C
C Reduce the step size if necessary so that the code will not step
C past TND. "Look ahead" to prevent unnecessarily small step sizes.
LAST = DIR*((T+H)-TND) .GE. ZERO
IF (LAST) THEN
H = TND - T
ELSE IF (DIR*((T+TWO*H)-TND).GE.ZERO) THEN
H = (TND-T)/TWO
END IF
C
C When the integrator is at T and attempts a step of H, the function
C defining the differential equations will be evaluated at a number of
C arguments between T and T+H. If H is too small, these arguments cannot
C be clearly distinguished in the precision available.
C
HMIN = MAX(TINY,TOOSML*MAX(ABS(T),ABS(T+H)))
IF (ABS(H).LT.HMIN) THEN
IER = 5
NREC = 3
WRITE (REC,'(A/A,D13.5,A,D13.5/A)')
&' ** In order to satisfy your error requirements CT would ',
&' ** have to use a step size of ',H,' at TNOW = ',T,
&' ** This is too small for the machine precision.'
GO TO 180
END IF
C
C Monitor the impact of output on the efficiency of the integration.
C
IF (CHKEFF) THEN
NTEND = NTEND + 1
IF (NTEND.GE.100 .AND. NTEND.GE.OKSTP/3) THEN
IER = 2
NREC = 6
WRITE (REC,'(A/A/A/A/A/A)')
&' ** More than 100 output points have been obtained by ',
&' ** integrating to TEND. They have been sufficiently close ',
&' ** to one another that the efficiency of the integration has ',
&' ** been degraded. It would probably be (much) more efficient ',
&' ** to obtain output by interpolating with INTRP (after ',
&' ** changing to METHOD=2 if you are using METHOD = 3).'
NTEND = 0
GO TO 180
END IF
END IF
C
C Check for stiffness and for too much work. Stiffness can be
C checked only after a successful step.
C
IF (.NOT.FAILED) THEN
C
C Check for too much work.
TOOMCH = NFCN .GT. MAXFCN
IF (TOOMCH) THEN
IER = 3
NREC = 3
WRITE (REC,'(A,I6,A/A/A)')
&' ** Approximately ',MAXFCN,' function evaluations have been ',
&' ** used to compute the solution since the integration ',
&' ** started or since this message was last printed.'
C
C After this warning message, NFCN is reset to permit the integration
C to continue. The total number of function evaluations in the primary
C integration is SVNFCN + NFCN.
SVNFCN = SVNFCN + NFCN
NFCN = 0
END IF
C
C Check for stiffness. NREC is passed on to STIFF because when
C TOOMCH = .TRUE. and stiffness is diagnosed, the message about too
C much work is augmented inside STIFF to explain that it is due to
C stiffness.
CALL STIFF(F,HAVG,JFLSTP,TOOMCH,MAXFCN,WORK,IER,NREC)
C
IF (IER.NE.1) GO TO 180
END IF
C
C Take a step. Whilst finding an on-scale H (PHASE2 = .TRUE.), the input
C value of H might be reduced (repeatedly), but it will not be reduced
C below HMIN. The local error is estimated, a weight vector is formed,
C and a weighted maximum norm, ERR, of the local error is returned.
C The variable MAIN is input as .TRUE. to tell STEP that this is the
C primary, or "main", integration.
C
C H resides in the common block /RKCOM2/ which is used by both CT and STEP;
C since it may be changed inside STEP, a local copy is made to ensure
C portability of the code.
C
MAIN = .TRUE.
HTRY = H
CALL STEP(F,NEQN,T,WORK(PRY),WORK(PRYP),WORK(PRSTGS),TOLR,HTRY,
& WORK(PRWT),WORK(YNEW),WORK(PRERST),ERR,MAIN,HMIN,
& WORK(PRTHRS),PHASE2)
H = HTRY
C
C Compare the norm of the local error to the tolerance.
C
IF (ERR.GT.TOLR) THEN
C
C Failed step. Reduce the step size and try again.
C
C First step: Terminate PHASE3 of the search for an on-scale step size.
C The step size is not on scale, so ERR may not be accurate;
C reduce H by a fixed factor. Failed attempts to take the
C first step are not counted.
C Later step: Use ERR to compute an "optimal" reduction of H. More than
C one failure indicates a difficulty with the problem and an
C ERR that may not be accurate, so reduce H by a fixed factor.
C
IF (FIRST) THEN
PHASE3 = .FALSE.
ALPHA = RS1
ELSE
FLSTP = FLSTP + 1
JFLSTP = JFLSTP + 1
IF (FAILED) THEN
ALPHA = RS1
ELSE
ALPHA = SAFETY*(TOLR/ERR)**EXPON
ALPHA = MAX(ALPHA,RS1)
END IF
END IF
H = ALPHA*H
FAILED = .TRUE.
GO TO 100
END IF
C
C Successful step.
C
C Predict a step size appropriate for the next step. After the first
C step the prediction can be refined using an idea of H.A. Watts that
C takes account of how well the prediction worked on the previous step.
BETA = (ERR/TOLR)**EXPON
IF (.NOT.FIRST) THEN
TEMP1 = (ERR**EXPON)/H
TEMP2 = (ERROLD**EXPON)/HOLD
IF (TEMP1.LT.TEMP2*HUNDRD .AND. TEMP2.LT.TEMP1*HUNDRD) THEN
BETA = BETA*(TEMP1/TEMP2)
END IF
END IF
ALPHA = RS3
IF (SAFETY.LT.BETA*ALPHA) ALPHA = SAFETY/BETA
C
C On the first step a search is made for an on-scale step size. PHASE2
C of the scheme comes to an end here because a step size has been found
C that is both successful and has a credible local error estimate. Except
C in the special case that the first step is also the last, the step is
C repeated in PHASE3 as long as an increase greater than RS2 appears
C possible. An increase as big as RS3 is permitted. A step failure
C terminates PHASE3.
C
IF (FIRST) THEN
PHASE2 = .FALSE.
PHASE3 = PHASE3 .AND. .NOT. LAST .AND. (ALPHA.GT.RS2)
IF (PHASE3) THEN
H = ALPHA*H
GO TO 100
END IF
END IF
C
C After getting on scale, step size changes are more restricted.
ALPHA = MIN(ALPHA,RS)
IF (FAILED) ALPHA = MIN(ALPHA,ONE)
ALPHA = MAX(ALPHA,RS1)
HOLD = H
H = ALPHA*H
C
C For the diagnosis of stiffness, an average accepted step size, HAVG,
C must be computed and SAVEd.
IF (FIRST) THEN
HAVG = HOLD
ELSE
HAVG = PT9*HAVG + PT1*HOLD
END IF
C
FIRST = .FALSE.
ERROLD = ERR
TOLD = T
C
C Take care that T is set to precisely TND when the end of the
C integration is reached.
IF (LAST) THEN
T = TND
ELSE
T = T + HOLD
END IF
C
C Increment counter on accepted steps. Note that successful steps
C that are repeated whilst getting on scale are not counted.
OKSTP = OKSTP + 1
C
C Advance the current solution and its derivative. (Stored in WORK(*)
C with the first location being PRY and PRYP, respectively.) Update the
C previous solution and its derivative. (Stored in WORK(*) with the first
C location being PRYOLD and YPOLD, respectively.) Note that the previous
C derivative will overwrite YNEW(*).
C
DO 120 L = 1, NEQN
WORK(PRYOLD-1+L) = WORK(PRY-1+L)
WORK(PRY-1+L) = WORK(YNEW-1+L)
WORK(YPOLD-1+L) = WORK(PRYP-1+L)
120 CONTINUE
C
IF (FSAL) THEN
C
C When FSAL = .TRUE., YP(*) is the last stage of the step.
POINT = PRSTGS + (LSTSTG-1)*NEQN
DO 140 L = 1, NEQN
WORK(PRYP-1+L) = WORK(POINT-1+L)
140 CONTINUE
ELSE
C
C Call F to evaluate YP(*).
CALL F(T,WORK(PRY),WORK(PRYP))
NFCN = NFCN + 1
END IF
C
C If global error assessment is desired, advance the secondary
C integration from TOLD to T.
C
IF (ERASON) THEN
CALL TRUERR(F,NEQN,WORK(PRY),TOLR,WORK(PRWT),WORK(PRZY),
& WORK(PRZYP),WORK(PRZERR),WORK(PRZYNU),WORK(PRZERS),
& WORK(PRZSTG),IER)
IF (IER.EQ.6) THEN
C
C The global error estimating procedure has broken down. Treat it as a
C failed step. The solution and derivative are reset to their values at
C the beginning of the step since the last valid error assessment refers
C to them.
OKSTP = OKSTP - 1
ERASFL = .TRUE.
LAST = .FALSE.
T = TOLD
H = HOLD
DO 160 L = 1, NEQN
WORK(PRY-1+L) = WORK(PRYOLD-1+L)
WORK(PRYP-1+L) = WORK(YPOLD-1+L)
160 CONTINUE
IF (OKSTP.GT.1) THEN
NREC = 2
WRITE (REC,'(A/A,D13.5,A)')
&' ** The global error assessment may not be reliable for T past ',
&' ** TNOW = ',T,'. The integration is being terminated.'
ELSE
NREC = 2
WRITE (REC,'(A/A)')
&' ** The global error assessment algorithm failed at the start',
&' ** the integration. The integration is being terminated.'
END IF
GO TO 180
END IF
END IF
C
C
C Exit point for CT
C
180 CONTINUE
C
C Set the output variables and flag that interpolation is permitted
C
IF (IER.LT.911) THEN
TNOW = T
LAST = TNOW .EQ. TND
CHKEFF = LAST
DO 200 L = 1, NEQN
YNOW(L) = WORK(PRY-1+L)
YPNOW(L) = WORK(PRYP-1+L)
200 CONTINUE
IF (IER.EQ.1) THEN
STATE = MINUS2
CALL RKSIT(TELL,'INTRP',STATE)
END IF
END IF
C
C Call RKMSG to report what happened and set CFLAG.
C
CALL RKMSG(IER,SRNAME,NREC,CFLAG)
C
RETURN
END SUBROUTINE CT
C
SUBROUTINE INTRP(TWANT,REQEST,NWANT,YWANT,YPWANT,F,WORK,WRKINT,
& LENINT)
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
C If you are not familiar with the code INTRP and how it is used in
C conjunction with CT to solve initial value problems, you should study the
C document file rksuite.doc carefully before attempting to use the code. The
C following "Brief Reminder" is intended only to remind you of the meaning,
C type, and size requirements of the arguments.
C
C When integrating with METHOD = 1 or 2, answers may be obtained inexpensively
C between steps by interpolation. INTRP is called after a step by CT from a
C previous value of T, called TOLD below, to the current value of T to get
C an answer at TWANT. You can specify any value of TWANT you wish, but
C specifying a value outside the interval [TOLD,T] is likely to yield
C answers with unsatisfactory accuracy.
C
C INPUT VARIABLE
C
C TWANT - DOUBLE PRECISION
C The value of the independent variable where a solution
C is desired.
C
C The interpolant is to be evaluated at TWANT to approximate the solution
C and/or its first derivative there. There are three cases:
C
C INPUT VARIABLE
C
C REQEST - CHARACTER*(*)
C Only the first character of REQEST is significant.
C REQEST(1:1) = `S' or `s'- compute approximate `S'olution
C only.
C = `D' or `d'- compute approximate first
C `D'erivative of the solution only.
C = `B' or `b'- compute `B'oth approximate solution
C and first derivative.
C Constraint: REQEST(1:1) must be `S',`s',`D',`d',`B' or `b'.
C
C If you intend to interpolate at many points, you should arrange for the
C the interesting components to be the first ones because the code
C approximates only the first NWANT components.
C
C INPUT VARIABLE
C
C NWANT - INTEGER
C Only the first NWANT components of the answer are to be
C computed.
C Constraint: NEQ >= NWANT >= 1
C
C OUTPUT VARIABLES
C
C YWANT(*) - DOUBLE PRECISION array of length NWANT
C Approximation to the first NWANT components of the true
C solution at TWANT when REQESTed.
C YPWANT(*) - DOUBLE PRECISION array of length NWANT
C Approximation to the first NWANT components of the first
C derivative of the true solution at TWANT when REQESTed.
C
C NAME DECLARED IN AN EXTERNAL STATEMENT IN THE PROGRAM CALLING INTRP:
C
C F - name of the subroutine for evaluating the differential
C equations as provided to CT.
C
C WORKSPACE
C
C WORK(*) - DOUBLE PRECISION array as used in SETUP and CT
C Do not alter the contents of this array.
C
C WRKINT(*) - DOUBLE PRECISION array of length LENINT
C Do not alter the contents of this array.
C
C LENINT - INTEGER
C Length of WRKINT. If
C METHOD = 1, LENINT must be at least 1
C = 2, LENINT must be at least NEQ+MAX(NEQ,5*NWANT)
C = 3--CANNOT BE USED WITH THIS SUBROUTINE
C
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C .. Scalar Arguments ..
DOUBLE PRECISION TWANT
INTEGER LENINT, NWANT
CHARACTER*(*) REQEST
C .. Array Arguments ..
DOUBLE PRECISION WORK(*), WRKINT(*), YPWANT(*), YWANT(*)
C .. Subroutine Arguments ..
EXTERNAL F
C .. Common Block for Problem Definition ..
DOUBLE PRECISION TSTRT, TND, DIR, HSTRT, TOLR
INTEGER NEQN
COMMON /RKCOM1/ TSTRT, TND, DIR, HSTRT, TOLR, NEQN
SAVE /RKCOM1/
C .. Common Block for General Workspace Pointers ..
INTEGER PRTHRS, PRERST, PRWT, PRYOLD, PRSCR, PRY, PRYP,
& PRSTGS, PRINTP, LNINTP
COMMON /RKCOM3/ PRTHRS, PRERST, PRWT, PRYOLD, PRSCR, PRY, PRYP,
& PRSTGS, PRINTP, LNINTP
SAVE /RKCOM3/
C .. Common Block to hold Formula Definitions ..
DOUBLE PRECISION A(13,13), B(13), C(13), BHAT(13), R(11,6),
& E(7)
INTEGER PTR(13), NSTAGE, METHD, MINTP
LOGICAL INTP
COMMON /RKCOM4/ A, B, C, BHAT, R, E, PTR, NSTAGE, METHD,
& MINTP, INTP
SAVE /RKCOM4/
C .. Common Block for Integrator Options ..
LOGICAL MSG, UTASK
COMMON /RKCOM8/ MSG, UTASK
SAVE /RKCOM8/
C .. Common Block for Error Message ..
CHARACTER*80 REC(10)
COMMON /RKCOM9/ REC
SAVE /RKCOM9/
C .. Parameters ..
CHARACTER*6 SRNAME
PARAMETER (SRNAME='INTRP')
LOGICAL ASK
INTEGER PLUS1, MINUS1, MINUS2
PARAMETER (ASK=.TRUE.,PLUS1=1,MINUS1=-1,MINUS2=-2)
C .. Local Scalars ..
INTEGER FLAG, ICHK, IER, NREC, NWNTSV, STARTP, STATE,
& STATE1
LOGICAL ININTP, LEGALR
CHARACTER REQST1
C .. External Subroutines ..
EXTERNAL EVALI, FORMI, RKMSG, RKSIT
C .. Intrinsic Functions ..
INTRINSIC MAX
C .. Save statement ..
SAVE NWNTSV, ININTP, STARTP
C .. Data statements ..
DATA NWNTSV/MINUS1/
C .. Executable Statements ..
C
IER = 1
NREC = 0
C
C Is it permissible to call INTRP?
C
CALL RKSIT(ASK,'CT',STATE)
IF (STATE.EQ.911) THEN
IER = 912
NREC = 1
WRITE (REC,'(A)')
&' ** A catastrophic error has already been detected elsewhere.'
GO TO 20
END IF
IF (UTASK) THEN
CALL RKSIT(ASK,'UT',STATE1)
IF (STATE1.NE.MINUS2) THEN
IER = 911
NREC = 2
WRITE (REC,'(A/A)')
&' ** You have called INTRP after you specified to SETUP ',
&' ** that you were going to use UT. This is not permitted.'
GO TO 20
END IF
END IF
IF (STATE.EQ.MINUS1) THEN
IER = 911
NREC = 1
WRITE (REC,'(A)')
&' ** You have not called CT, so you cannot use INTRP.'
GO TO 20
END IF
IF (STATE.GT.PLUS1) THEN
IER = 911
NREC = 2
WRITE (REC,'(A/A)')
&' ** CT has returned with a flag value greater than 1.',
&' ** You cannot call INTRP in this circumstance.'
GO TO 20
END IF
C
C Check input
C
REQST1 = REQEST(1:1)
LEGALR = REQST1 .EQ. 'S' .OR. REQST1 .EQ. 's' .OR.
& REQST1 .EQ. 'D' .OR. REQST1 .EQ. 'd' .OR.
& REQST1 .EQ. 'B' .OR. REQST1 .EQ. 'b'
IF (.NOT.LEGALR) THEN
IER = 911
NREC = 3
WRITE (REC,'(A/A,A,A/A)')
&' ** You have set the first character of ',
&' ** REQEST to be ''',REQST1,'''. It must be one of ',
&' ** ''S'',''s'',''D'',''d'',''B'' or ''b''.'
GO TO 20
END IF
C
IF (NWANT.GT.NEQN) THEN
IER = 911
NREC = 3
WRITE (REC,'(A,I6,A/A,I6,A/A)')
&' ** You have specified the value of NWANT to be ',NWANT,'. This',
&' ** is greater than ',NEQN,', which is the number of equations ',
&' ** in the system being integrated.'
GO TO 20
ELSE IF (NWANT.LT.1) THEN
IER = 911
NREC = 3
WRITE (REC,'(A,I6,A/A/A)')
&' ** You have specified the value of NWANT to be ',NWANT,', but ',
&' ** this is less than 1. You cannot interpolate a zero or ',
&' ** negative number of components.'
GO TO 20
END IF
C
IF (METHD.EQ.1) THEN
IF (LENINT.LT.1) THEN
IER = 911
NREC = 2
WRITE (REC,'(A,I6,A/A)')
&' ** You have specified LENINT to be ',LENINT,'.',
&' ** This is too small. LENINT must be at least 1.'
GO TO 20
END IF
STARTP = 1
ELSE IF (METHD.EQ.2) THEN
ICHK = NEQN + MAX(NEQN,5*NWANT)
IF (LENINT.LT.ICHK) THEN
IER = 911
NREC = 3
WRITE (REC,'(A,I6,A/A/A,I6,A)')
&' ** You have specified LENINT to be ',LENINT,'. This is too',
&' ** small. NINT must be at least NEQ + MAX(NEQ, 5*NWANT) ',
&' ** which is ', ICHK,'.'
GO TO 20
END IF
STARTP = NEQN + 1
ELSE IF (METHD.EQ.3) THEN
IER = 911
NREC = 5
WRITE (REC,'(A/A/A/A/A)')
&' ** You have been using CT with METHOD = 3 to integrate your ',
&' ** equations. You have just called INTRP, but interpolation ',
&' ** is not available for this METHOD. Either use METHOD = 2, ',
&' ** for which interpolation is available, or use RESET to make',
&' ** CT step exactly to the points where you want output.'
GO TO 20
END IF
C
C Has the interpolant been initialised for this step?
C
CALL RKSIT(ASK,SRNAME,STATE)
ININTP = STATE .NE. MINUS2
C
C Some initialization must be done before interpolation is possible.
C To reduce the overhead, the interpolating polynomial is formed for
C the first NWANT components. In the unusual circumstance that NWANT
C is changed while still interpolating within the span of the current
C step, the scheme must be reinitialized to accomodate the additional
C components.
C
IF (.NOT.ININTP .OR. NWANT.NE.NWNTSV) THEN
C
C At present the derivative of the solution at the previous step, YPOLD(*),
C is stored in the scratch array area starting at PRSCR. In the case of
C METHD = 1 we can overwrite the stages.
C
IF (METHD.EQ.1) THEN
CALL FORMI(F,NEQN,NWANT,WORK(PRY),WORK(PRYP),WORK(PRYOLD),
& WORK(PRSCR),WORK(PRSTGS),.NOT.ININTP,
& WORK(PRSTGS),WORK(PRSTGS))
ELSE
CALL FORMI(F,NEQN,NWANT,WORK(PRY),WORK(PRYP),WORK(PRYOLD),
& WORK(PRSCR),WORK(PRSTGS),.NOT.ININTP,WRKINT,
& WRKINT(STARTP))
END IF
C
C Set markers to show that interpolation has been initialized for
C NWANT components.
NWNTSV = NWANT
ININTP = .TRUE.
END IF
C
C The actual evaluation of the interpolating polynomial and/or its first
C derivative is done in EVALI.
C
IF (METHD.EQ.1) THEN
CALL EVALI(WORK(PRY),WORK(PRYP),WORK(PRSTGS),TWANT,REQEST,
& NWANT,YWANT,YPWANT)
ELSE
CALL EVALI(WORK(PRY),WORK(PRYP),WRKINT(STARTP),TWANT,REQEST,
& NWANT,YWANT,YPWANT)
END IF
C
20 CONTINUE
C
CALL RKMSG(IER,SRNAME,NREC,FLAG)
C
RETURN
END SUBROUTINE INTRP
C
SUBROUTINE RESET(TENDNU)
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
C If you are not familiar with the code RESET and how it is used in
C conjunction with CT to solve initial value problems, you should study the
C document file rksuite.doc carefully before attempting to use the code. The
C following "Brief Reminder" is intended only to remind you of the meaning,
C type, and size requirements of the arguments.
C
C The integration using CT proceeds from TSTART in the direction of TEND, and
C is now at TNOW. To reset TEND to a new value TENDNU, just call RESET with
C TENDNU as the argument. You must continue integrating in the same
C direction, so the sign of (TENDNU - TNOW) must be the same as the sign of
C (TEND - TSTART). To change direction you must restart by a call to SETUP.
C
C INPUT VARIABLE
C
C TENDNU - DOUBLE PRECISION
C The new value of TEND.
C Constraint: TENDNU and TNOW must be clearly distinguishable
C in the precision used. The sign of (TENDNU - TNOW) must be
C the same as the sign of (TEND - TSTART).
C
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C .. Scalar Arguments ..
DOUBLE PRECISION TENDNU
C .. Common Block for Problem Definition ..
DOUBLE PRECISION TSTRT, TND, DIR, HSTRT, TOLR
INTEGER NEQN
COMMON /RKCOM1/ TSTRT, TND, DIR, HSTRT, TOLR, NEQN
SAVE /RKCOM1/
C .. Common Block to hold Problem Status ..
DOUBLE PRECISION T, H, TOLD, HOLD
INTEGER NFCN, SVNFCN, OKSTP, FLSTP
LOGICAL FIRST, LAST
COMMON /RKCOM2/ T, H, TOLD, HOLD, NFCN, SVNFCN, OKSTP, FLSTP,
& FIRST, LAST
SAVE /RKCOM2/
C .. Common Block to hold Formula Characterisitcs ..
DOUBLE PRECISION TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
& RS, RS1, RS2, RS3, RS4
INTEGER ORDER, LSTSTG, MAXTRY, NSEC
LOGICAL FSAL
COMMON /RKCOM5/ TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
& RS, RS1, RS2, RS3, RS4, ORDER, LSTSTG, MAXTRY,
& NSEC, FSAL
SAVE /RKCOM5/
C .. Common Block for Environment Parameters ..
DOUBLE PRECISION MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY
INTEGER OUTCH
COMMON /RKCOM7/ MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY,
& OUTCH
SAVE /RKCOM7/
C .. Common Block for Integrator Options ..
LOGICAL MSG, UTASK
COMMON /RKCOM8/ MSG, UTASK
SAVE /RKCOM8/
C .. Common Block for Error Message ..
CHARACTER*80 REC(10)
COMMON /RKCOM9/ REC
SAVE /RKCOM9/
C .. Parameters ..
CHARACTER*6 SRNAME
PARAMETER (SRNAME='RESET')
LOGICAL ASK
INTEGER MINUS1, MINUS2
PARAMETER (ASK=.TRUE.,MINUS1=-1,MINUS2=-2)
DOUBLE PRECISION ZERO
PARAMETER (ZERO=0.0D+0)
C .. Local Scalars ..
DOUBLE PRECISION HMIN, TDIFF
INTEGER FLAG, IER, NREC, STATE, STATE1
C .. External Subroutines ..
EXTERNAL RKMSG, RKSIT
C .. Intrinsic Functions ..
INTRINSIC ABS, MAX
C .. Executable Statements ..
IER = 1
NREC = 0
C
C Is it permissible to call RESET?
C
CALL RKSIT(ASK,'CT',STATE)
IF (STATE.EQ.911) THEN
IER = 912
NREC = 1
WRITE (REC,'(A)')
&' ** A catastrophic error has already been detected elsewhere.'
GO TO 20
END IF
IF (UTASK) THEN
CALL RKSIT(ASK,'UT',STATE1)
IF (STATE1.NE.MINUS2) THEN
IER = 911
NREC = 2
WRITE (REC,'(A/A)')
&' ** You have called RESET after you specified to SETUP that ',
&' ** you were going to use UT. This is not permitted.'
GO TO 20
END IF
END IF
IF (STATE.EQ.MINUS1) THEN
IER = 911
NREC = 1
WRITE (REC,'(A)')
&' ** You have not called CT, so you cannot use RESET.'
GO TO 20
END IF
IF (STATE.EQ.5 .OR. STATE.EQ.6) THEN
IER = 911
NREC = 2
WRITE (REC,'(A,I1,A/A)')
&' ** CT has returned with CFLAG = ',STATE,'.',
&' ** You cannot call RESET in this circumstance.'
GO TO 20
END IF
C
C Check value of TENDNU
C
IF (DIR.GT.ZERO .AND. TENDNU.LE.T) THEN
IER = 911
NREC = 4
WRITE (REC,'(A/A,D13.5/A,D13.5,A/A)')
&' ** Integration is proceeding in the positive direction. The ',
&' ** current value for the independent variable is ',T,
&' ** and you have set TENDNU = ',TENDNU,'. TENDNU must be ',
&' ** greater than T.'
ELSE IF (DIR.LT.ZERO .AND. TENDNU.GE.T) THEN
IER = 911
NREC = 4
WRITE (REC,'(A/A,D13.5/A,D13.5,A/A)')
&' ** Integration is proceeding in the negative direction. The ',
&' ** current value for the independent variable is ',T,
&' ** and you have set TENDNU = ',TENDNU,'. TENDNU must be ',
&' ** less than T.'
ELSE
HMIN = MAX(TINY,TOOSML*MAX(ABS(T),ABS(TENDNU)))
TDIFF = ABS(TENDNU-T)
IF (TDIFF.LT.HMIN) THEN
IER = 911
NREC = 4
WRITE (REC,'(A,D13.5,A/A,D13.5,A/A/A,D13.5,A)')
&' ** The current value of the independent variable T is ',T,'.',
&' ** The TENDNU you supplied has ABS(TENDNU-T) = ',TDIFF,'.',
&' ** For the METHOD and the precision of the computer being ',
&' ** used, this difference must be at least ',HMIN,'.'
END IF
END IF
IF (IER.EQ.911) GO TO 20
C
TND = TENDNU
LAST = .FALSE.
C
20 CONTINUE
C
CALL RKMSG(IER,SRNAME,NREC,FLAG)
C
RETURN
END SUBROUTINE RESET
C
SUBROUTINE MCONST(METHOD)
C************************************************
C**** NOT A DESIGNATED USER-CALLABLE ROUTINE ****
C************************************************
C
C Purpose: Sets machine-dependent global quantities
C
C Common: Initializes: /RKCOM7/ OUTCH, MCHEPS, DWARF, RNDOFF,
C SQRRMC, CUBRMC, TINY
C Reads: none
C Alters: none
C
C Comments:
C =========
C OUTCH, MCHEPS, DWARF are pure environmental parameters with values
C obtained from a call to ENVIRN. The other quantities set depend on
C the environmental parameters, the implementation, and, possibly,
C METHOD. At present the METHODs implemented in the RK suite do not
C influence the values of these quantities.
C OUTCH - Standard output channel
C MCHEPS - Largest positive number such that 1.0D0 + MCHEPS = 1.0D0
C DWARF - Smallest positive number
C RNDOFF - 10 times MCHEPS
C SQRRMC - Square root of MCHEPS
C CUBRMC - Cube root of MCHEPS
C TINY - Square root of DWARF
C
C .. Scalar Arguments ..
INTEGER METHOD
C .. Common Block for Environment Parameters ..
DOUBLE PRECISION MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY
INTEGER OUTCH
COMMON /RKCOM7/ MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY,
& OUTCH
SAVE /RKCOM7/
C .. Parameters ..
DOUBLE PRECISION TEN, THIRD
PARAMETER (TEN=10.0D+0,THIRD=1.0D+0/3.0D+0)
C .. External Subroutines ..
EXTERNAL ENVIRN
C .. Intrinsic Functions ..
INTRINSIC SQRT
C .. Executable Statements ..
C
CALL ENVIRN(OUTCH,MCHEPS,DWARF)
C
RNDOFF = TEN*MCHEPS
SQRRMC = SQRT(MCHEPS)
CUBRMC = MCHEPS**THIRD
TINY = SQRT(DWARF)
C
RETURN
END SUBROUTINE MCONST
C
SUBROUTINE ENVIRN(OUTCH,MCHEPS,DWARF)
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
C The RK suite requires some environmental parameters that are provided by
C this subroutine. The values provided with the distribution codes are those
C appropriate to the IEEE standard. They must be altered, if necessary, to
C those appropriate to the computing system you are using before calling the
C codes of the suite.
C
C ================================================================
C ================================================================
C TO MAKE SURE THAT THESE MACHINE AND INSTALLATION DEPENDENT
C QUANTITIES ARE SPECIFIED PROPERLY, THE DISTRIBUTION VERSION
C WRITES A MESSAGE ABOUT THE MATTER TO THE STANDARD OUTPUT CHANNEL
C AND TERMINATES THE RUN. THE VALUES PROVIDED IN THE DISTRIBUTION
C VERSION SHOULD BE ALTERED, IF NECESSARY, AND THE "WRITE" AND
C "STOP" STATEMENTS COMMENTED OUT.
C ================================================================
C ================================================================
C
C OUTPUT VARIABLES
C
C OUTCH - INTEGER
C Standard output channel
C MCHEPS - DOUBLE PRECISION
C MCHEPS is the largest positive number such that
C 1.0D0 + MCHEPS = 1.0D0.
C DWARF - DOUBLE PRECISION
C DWARF is the smallest positive number.
C
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
C .. Scalar Arguments ..
INTEGER OUTCH
DOUBLE PRECISION DWARF, MCHEPS
C .. Executable Statements ..
C
C The following six statements are to be Commented out after verification that
C the machine and installation dependent quantities are specified correctly.
C If you pass copies of RKSUITE on to others, please give them the whole
C distribution version of RKSUITE, and in particular, give them a version
C of ENVIRN that does not have the following six statements Commented out.
C WRITE(*,*) ' Before using RKSUITE, you must verify that the '
C WRITE(*,*) ' machine- and installation-dependent quantities '
C WRITE(*,*) ' specified in the subroutine ENVIRN are correct, '
C WRITE(*,*) ' and then Comment these WRITE statements and the '
C WRITE(*,*) ' STOP statement out of ENVIRN. '
C STOP
C
C The following values are appropriate to IEEE arithmetic with the typical
C standard output channel.
C
OUTCH = 6
MCHEPS = 1.11D-16
DWARF = 2.23D-308
C
C------------------------------------------------------------------------------
C If you have the routines D1MACH and I1MACH on your system, you could
C replace the preceding statements by the following ones to obtain the
C appropriate machine dependent numbers. The routines D1MACH and I1MACH
C are public domain software. They are available from NETLIB.
C .. Scalar Arguments ..
C INTEGER OUTCH
C DOUBLE PRECISION DWARF, MCHEPS
C .. External Functions ..
C INTEGER I1MACH
C DOUBLE PRECISION D1MACH
C .. Executable Statements ..
C
C OUTCH = I1MACH(2)
C MCHEPS = D1MACH(3)
C DWARF = D1MACH(1)
C
C If you have the NAG Fortran Library available on your system, you could
C replace the preceding statements by the following ones to obtain the
C appropriate machine dependent numbers.
C
C .. Scalar Arguments ..
C INTEGER OUTCH
C DOUBLE PRECISION DWARF, MCHEPS
C .. External Functions ..
C DOUBLE PRECISION X02AJF, X02AMF
C .. Executable Statements ..
C
C CALL X04AAF(0,OUTCH)
C MCHEPS = X02AJF()
C DWARF = X02AMF()
C
C If you have the IMSL MATH/LIBRARY available on your system, you could
C replace the preceding statements by the following ones to obtain the
C appropriate machine dependent numbers.
C
C .. Scalar Arguments ..
C INTEGER OUTCH
C DOUBLE PRECISION DWARF, MCHEPS
C .. External Functions ..
C DOUBLE PRECISION DMACH
C .. Executable Statements ..
C
C CALL UMACH(2,OUTCH)
C MCHEPS = DMACH(4)
C DWARF = DMACH(1)
C------------------------------------------------------------------------------
C
RETURN
END SUBROUTINE ENVIRN
C
SUBROUTINE CONST(METHOD,VECSTG,REQSTG,LINTPL)
C************************************************
C**** NOT A DESIGNATED USER-CALLABLE ROUTINE ****
*************************************************
C
C Purpose: Set formula definitions and formula characteristics for
C selected method. Return storage requirements for the
C selected method.
C
C Input: METHOD
C Output: VECSTG, REQSTG, LINTPL
C
C Common: Initializes: /RKCOM4/ A(*,*), B(*), C(*), BHAT(*), R(*),
C E(*), PTR(*), NSTAGE, METHD, INTP, MINTP
C /RKCOM5/ TOOSML, COST, SAFETY, EXPON, STBRAD,
C TANANG, RS, RS1, RS2, RS3, RS4, ORDER,
C LSTSTG, MAXTRY, NSEC, FSAL
C Reads: /RKCOM7/ RNDOFF
C Alters: none
C
C Comments:
C =========
C Runge-Kutta formula pairs are described by a set of coefficients
C and by the setting of a number of parameters that describe the
C characteristics of the pair. The input variable METHD indicates
C which of the three pairs available is to be set. In the case of
C METHD = 2 additional coefficients are defined that make interpolation
C of the results possible and provide an additional error estimator.
C VECSTG is the number of columns of workspace required to compute the
C stages of a METHD. For interpolation purposes the routine returns via
C COMMON the logical variable INTP indicating whether interpolation is
C possible and via the call list:
C REQSTG - whether the stages are required to form the
C interpolant (set .FALSE. if INTP=.FALSE.)
C LINTPL - the number of extra columns of storage required for use
C with UT (set 0 if INTP=.FALSE.)
C
C Quantities set in common blocks:
C METHD - copy of METHOD
C A, B, C, BHAT - coefficients of the selected method
C R - extra coefficents for interpolation with METHD = 2
C E - extra coefficients for additional local error estimate
C with METHD = 2
C PTR - vector of pointers indicating how individual stages are to
C be stored. With it zero coefficients of the formulas can
C be exploited to reduce the storage required
C NSTAGE - number of stages for the specified METHD
C INTP - indicates whether there is an associated interpolant
C (depending on the method, the user may have to supply
C extra workspace)
C MINTP - the degree of the interpolating polynomial, if one exists
C FSAL - indicates whether the last stage of a step can be used as
C the first stage of the following step
C LSTSTG - pointer to location of last stage for use with FSAL=.TRUE.
C ORDER - the lower order of the pair of Runge-Kutta formulas that
C constitute a METHD
C TANANG,
C STBRAD - the stability region of the formula used to advance
C the integration is approximated by a sector in the left half
C complex plane. TANANG is the tangent of the interior angle
C of the sector and STBRAD is the radius of the sector.
C COST - cost of a successful step in function evaluations
C MAXTRY - limit on the number of iterations in the stiffness check. As
C set, no more than 24 function evaluations are made in the check.
C NSEC - each step of size H in the primary integration corresponds to
C NSEC steps of size H/NSEC in the secondary integration when
C global error assessment is done.
C EXPON - used to adjust the step size; this code implements an error
C per step control for which EXPON = 1/(ORDER + 1).
C SAFETY - quantity used in selecting the step size
C TOOSML - quantity used to determine when a step size is too small for
C the precision available
C RS, RS1,
C RS2, RS3,
C RS4 - quantities used in determining the maximum and minimum change
C change in step size (set independently of METHD)
C
C Further comments on SAFETY:
C ============================
C The code estimates the largest step size that will yield the specified
C accuracy. About half the time this step size would result in a local
C error that is a little too large, and the step would be rejected. For
C this reason a SAFETY factor is used to reduce the "optimal" value to one
C more likely to succeed. Unfortunately, there is some art in choosing this
C value. The more expensive a failed step, the smaller one is inclined to
C take this factor. However, a small factor means that if the prediction were
C good, more accuracy than desired would be obtained and the behavior of the
C error would then be irregular. The more stringent the tolerance, the better
C the prediction, except near limiting precision. Thus the general range of
C tolerances expected influences the choice of SAFETY.
C
C .. Scalar Arguments ..
INTEGER LINTPL, METHOD, VECSTG
LOGICAL REQSTG
C .. Common Block to hold Formula Definitions ..
DOUBLE PRECISION A(13,13), B(13), C(13), BHAT(13), R(11,6),
& E(7)
INTEGER PTR(13), NSTAGE, METHD, MINTP
LOGICAL INTP
COMMON /RKCOM4/ A, B, C, BHAT, R, E, PTR, NSTAGE, METHD,
& MINTP, INTP
SAVE /RKCOM4/
C .. Common Block to hold Formula Characterisitcs ..
DOUBLE PRECISION TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
& RS, RS1, RS2, RS3, RS4
INTEGER ORDER, LSTSTG, MAXTRY, NSEC
LOGICAL FSAL
COMMON /RKCOM5/ TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
& RS, RS1, RS2, RS3, RS4, ORDER, LSTSTG, MAXTRY,
& NSEC, FSAL
SAVE /RKCOM5/
C .. Common Block for Environment Parameters ..
DOUBLE PRECISION MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY
INTEGER OUTCH
COMMON /RKCOM7/ MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY,
& OUTCH
SAVE /RKCOM7/
C .. Parameters ..
DOUBLE PRECISION ONE, ZERO, TWO, FIFTY, FIVEPC
PARAMETER (ONE=1.0D+0,ZERO=0.0D+0,TWO=2.0D+0,FIFTY=50.D+0,
& FIVEPC=0.05D+0)
C .. Local Scalars ..
DOUBLE PRECISION CDIFF, DIFF
INTEGER I, J
C .. Intrinsic Functions ..
INTRINSIC ABS, DBLE, INT, MAX, MIN
C .. Executable Statements ..
C
METHD = METHOD
C
GO TO (20,40,100) METHD
C
C METHD = 1.
C This pair is from "A 3(2) Pair of Runge-Kutta Formulas" by P. Bogacki
C and L.F. Shampine, Appl. Math. Lett., 2, pp. 321-325, 1989. The authors
C are grateful to P. Bogacki for his assistance in implementing the pair.
C
20 CONTINUE
NSTAGE = 4
FSAL = .TRUE.
ORDER = 2
TANANG = 8.9D0
STBRAD = 2.3D0
SAFETY = 0.8D0
INTP = .TRUE.
MINTP = 3
REQSTG = .FALSE.
LINTPL = 2
NSEC = 3
C
PTR(1) = 0
PTR(2) = 1
PTR(3) = 2
PTR(4) = 3
C
A(2,1) = 1.D0/2.D0
A(3,1) = 0.D0
A(3,2) = 3.D0/4.D0
A(4,1) = 2.D0/9.D0
A(4,2) = 1.D0/3.D0
A(4,3) = 4.D0/9.D0
C
C The coefficients BHAT(*) refer to the formula used to advance the
C integration, here the one of order 3. The coefficients B(*) refer
C to the other formula, here the one of order 2. For this pair, BHAT(*)
C is not needed since FSAL = .TRUE.
C
B(1) = 7.D0/24.D0
B(2) = 1.D0/4.D0
B(3) = 1.D0/3.D0
B(4) = 1.D0/8.D0
C
C(1) = 0.D0
C(2) = 1.D0/2.D0
C(3) = 3.D0/4.D0
C(4) = 1.D0
C
GO TO 120
C
C METHD = 2
C This pair is from "An Efficient Runge-Kutta (4,5) Pair" by P. Bogacki
C and L.F. Shampine, Rept. 89-20, Math. Dept., Southern Methodist
C University, Dallas, Texas, USA, 1989. The authors are grateful to
C P. Bogacki for his assistance in implementing the pair. Shampine and
C Bogacki subsequently modified the formula to enhance the reliability of
C the pair. The original fourth order formula is used in an estimate of
C the local error. If the step fails, the computation is broken off. If
C the step is acceptable, the first evaluation of the next step is done,
C i.e., the pair is implemented as FSAL and the local error of the step
C is again estimated with a fourth order formula using the additional data.
C The step must succeed with both estimators to be accepted. When the
C second estimate is formed, it is used for the subsequent adjustment of
C the step size because it is of higher quality. The two fourth order
C formulas are well matched to leading order, and only exceptionally do
C the estimators disagree -- problems with discontinuous coefficients are
C handled more reliably by using two estimators as is global error
C estimation.
C
40 CONTINUE
NSTAGE = 8
FSAL = .TRUE.
ORDER = 4
TANANG = 5.2D0
STBRAD = 3.9D0
SAFETY = 0.8D0
INTP = .TRUE.
REQSTG = .TRUE.
MINTP = 6
LINTPL = 6
NSEC = 2
C
PTR(1) = 0
PTR(2) = 1
PTR(3) = 2
PTR(4) = 3
PTR(5) = 4
PTR(6) = 5
PTR(7) = 6
PTR(8) = 7
C
A(2,1) = 1.D0/6.D0
A(3,1) = 2.D0/27.D0
A(3,2) = 4.D0/27.D0
A(4,1) = 183.D0/1372.D0
A(4,2) = -162.D0/343.D0
A(4,3) = 1053.D0/1372.D0
A(5,1) = 68.D0/297.D0
A(5,2) = -4.D0/11.D0
A(5,3) = 42.D0/143.D0
A(5,4) = 1960.D0/3861.D0
A(6,1) = 597.D0/22528.D0
A(6,2) = 81.D0/352.D0
A(6,3) = 63099.D0/585728.D0
A(6,4) = 58653.D0/366080.D0
A(6,5) = 4617.D0/20480.D0
A(7,1) = 174197.D0/959244.D0
A(7,2) = -30942.D0/79937.D0
A(7,3) = 8152137.D0/19744439.D0
A(7,4) = 666106.D0/1039181.D0
A(7,5) = -29421.D0/29068.D0
A(7,6) = 482048.D0/414219.D0
A(8,1) = 587.D0/8064.D0
A(8,2) = 0.D0
A(8,3) = 4440339.D0/15491840.D0
A(8,4) = 24353.D0/124800.D0
A(8,5) = 387.D0/44800.D0
A(8,6) = 2152.D0/5985.D0
A(8,7) = 7267.D0/94080.D0
C
C The coefficients B(*) refer to the formula of order 4.
C
B(1) = 2479.D0/34992.D0
B(2) = 0.D0
B(3) = 123.D0/416.D0
B(4) = 612941.D0/3411720.D0
B(5) = 43.D0/1440.D0
B(6) = 2272.D0/6561.D0
B(7) = 79937.D0/1113912.D0
B(8) = 3293.D0/556956.D0
C
C The coefficients E(*) refer to an estimate of the local error based on
C the first formula of order 4. It is the difference of the fifth order
C result, here located in A(8,*), and the fourth order result. By
C construction both E(2) and E(7) are zero.
C
E(1) = -3.D0/1280.D0
E(2) = 0.D0
E(3) = 6561.D0/632320.D0
E(4) = -343.D0/20800.D0
E(5) = 243.D0/12800.D0
E(6) = -1.D0/95.D0
E(7) = 0.D0
C
C(1) = 0.D0
C(2) = 1.D0/6.D0
C(3) = 2.D0/9.D0
C(4) = 3.D0/7.D0
C(5) = 2.D0/3.D0
C(6) = 3.D0/4.D0
C(7) = 1.D0
C(8) = 1.D0
C
C To do interpolation with this pair, some extra stages have to be computed.
C The following additional A(*,*) and C(*) coefficients are for this purpose.
C In addition there is an array R(*,*) that plays a role for interpolation
C analogous to that of BHAT(*) for the basic step.
C
C(9) = 1.D0/2.D0
C(10) = 5.D0/6.D0
C(11) = 1.D0/9.D0
C
A(9,1) = 455.D0/6144.D0
A(10,1) = -837888343715.D0/13176988637184.D0
A(11,1) = 98719073263.D0/1551965184000.D0
A(9,2) = 0.D0
A(10,2) = 30409415.D0/52955362.D0
A(11,2) = 1307.D0/123552.D0
A(9,3) = 10256301.D0/35409920.D0
A(10,3) = -48321525963.D0/759168069632.D0
A(11,3) = 4632066559387.D0/70181753241600.D0
A(9,4) = 2307361.D0/17971200.D0
A(10,4) = 8530738453321.D0/197654829557760.D0
A(11,4) = 7828594302389.D0/382182512025600.D0
A(9,5) = -387.D0/102400.D0
A(10,5) = 1361640523001.D0/1626788720640.D0
A(11,5) = 40763687.D0/11070259200.D0
A(9,6) = 73.D0/5130.D0
A(10,6) = -13143060689.D0/38604458898.D0
A(11,6) = 34872732407.D0/224610586200.D0
A(9,7) = -7267.D0/215040.D0
A(10,7) = 18700221969.D0/379584034816.D0
A(11,7) = -2561897.D0/30105600.D0
A(9,8) = 1.D0/32.D0
A(10,8) = -5831595.D0/847285792.D0
A(11,8) = 1.D0/10.D0
A(10,9) = -5183640.D0/26477681.D0
A(11,9) = -1.D0/10.D0
A(11,10) = -1403317093.D0/11371610250.D0
C
DO 60 I = 1, 11
R(I,1) = 0.D0
60 CONTINUE
DO 80 I = 1, 6
R(2,I) = 0.D0
80 CONTINUE
R(1,6) = -12134338393.D0/1050809760.D0
R(1,5) = -1620741229.D0/50038560.D0
R(1,4) = -2048058893.D0/59875200.D0
R(1,3) = -87098480009.D0/5254048800.D0
R(1,2) = -11513270273.D0/3502699200.D0
C
R(3,6) = -33197340367.D0/1218433216.D0
R(3,5) = -539868024987.D0/6092166080.D0
R(3,4) = -39991188681.D0/374902528.D0
R(3,3) = -69509738227.D0/1218433216.D0
R(3,2) = -29327744613.D0/2436866432.D0
C
R(4,6) = -284800997201.D0/19905339168.D0
R(4,5) = -7896875450471.D0/165877826400.D0
R(4,4) = -333945812879.D0/5671036800.D0
R(4,3) = -16209923456237.D0/497633479200.D0
R(4,2) = -2382590741699.D0/331755652800.D0
C
R(5,6) = -540919.D0/741312.D0
R(5,5) = -103626067.D0/43243200.D0
R(5,4) = -633779.D0/211200.D0
R(5,3) = -32406787.D0/18532800.D0
R(5,2) = -36591193.D0/86486400.D0
C
R(6,6) = 7157998304.D0/374350977.D0
R(6,5) = 30405842464.D0/623918295.D0
R(6,4) = 183022264.D0/5332635.D0
R(6,3) = -3357024032.D0/1871754885.D0
R(6,2) = -611586736.D0/89131185.D0
C
R(7,6) = -138073.D0/9408.D0
R(7,5) = -719433.D0/15680.D0
R(7,4) = -1620541.D0/31360.D0
R(7,3) = -385151.D0/15680.D0
R(7,2) = -65403.D0/15680.D0
C
R(8,6) = 1245.D0/64.D0
R(8,5) = 3991.D0/64.D0
R(8,4) = 4715.D0/64.D0
R(8,3) = 2501.D0/64.D0
R(8,2) = 149.D0/16.D0
R(8,1) = 1.D0
C
R(9,6) = 55.D0/3.D0
R(9,5) = 71.D0
R(9,4) = 103.D0
R(9,3) = 199.D0/3.D0
R(9,2) = 16.0D0
C
R(10,6) = -1774004627.D0/75810735.D0
R(10,5) = -1774004627.D0/25270245.D0
R(10,4) = -26477681.D0/359975.D0
R(10,3) = -11411880511.D0/379053675.D0
R(10,2) = -423642896.D0/126351225.D0
C
R(11,6) = 35.D0
R(11,5) = 105.D0
R(11,4) = 117.D0
R(11,3) = 59.D0
R(11,2) = 12.D0
C
GO TO 120
C
C METHD = 3
C This pair is from "High Order Embedded Runge-Kutta Formulae" by P.J.
C Prince and J.R. Dormand, J. Comp. Appl. Math.,7, pp. 67-75, 1981. The
C authors are grateful to P. Prince and J. Dormand for their assistance in
C implementing the pair.
C
100 CONTINUE
NSTAGE = 13
FSAL = .FALSE.
ORDER = 7
TANANG = 11.0D0
STBRAD = 5.2D0
SAFETY = 0.8D0
INTP = .FALSE.
REQSTG = .FALSE.
MINTP = 0
LINTPL = 0
NSEC = 2
C
PTR(1) = 0
PTR(2) = 1
PTR(3) = 2
PTR(4) = 1
PTR(5) = 3
PTR(6) = 2
PTR(7) = 4
PTR(8) = 5
PTR(9) = 6
PTR(10) = 7
PTR(11) = 8
PTR(12) = 9
PTR(13) = 1
C
A(2,1) = 5.55555555555555555555555555556D-2
A(3,1) = 2.08333333333333333333333333333D-2
A(3,2) = 6.25D-2
A(4,1) = 3.125D-2
A(4,2) = 0.D0
A(4,3) = 9.375D-2
A(5,1) = 3.125D-1
A(5,2) = 0.D0
A(5,3) = -1.171875D0
A(5,4) = 1.171875D0
A(6,1) = 3.75D-2
A(6,2) = 0.D0
A(6,3) = 0.D0
A(6,4) = 1.875D-1
A(6,5) = 1.5D-1
A(7,1) = 4.79101371111111111111111111111D-2
A(7,2) = 0.D0
A(7,3) = 0.0D0
A(7,4) = 1.12248712777777777777777777778D-1
A(7,5) = -2.55056737777777777777777777778D-2
A(7,6) = 1.28468238888888888888888888889D-2
A(8,1) = 1.6917989787292281181431107136D-2
A(8,2) = 0.D0
A(8,3) = 0.D0
A(8,4) = 3.87848278486043169526545744159D-1
A(8,5) = 3.59773698515003278967008896348D-2
A(8,6) = 1.96970214215666060156715256072D-1
A(8,7) = -1.72713852340501838761392997002D-1
A(9,1) = 6.90957533591923006485645489846D-2
A(9,2) = 0.D0
A(9,3) = 0.D0
A(9,4) = -6.34247976728854151882807874972D-1
A(9,5) = -1.61197575224604080366876923982D-1
A(9,6) = 1.38650309458825255419866950133D-1
A(9,7) = 9.4092861403575626972423968413D-1
A(9,8) = 2.11636326481943981855372117132D-1
A(10,1) = 1.83556996839045385489806023537D-1
A(10,2) = 0.D0
A(10,3) = 0.D0
A(10,4) = -2.46876808431559245274431575997D0
A(10,5) = -2.91286887816300456388002572804D-1
A(10,6) = -2.6473020233117375688439799466D-2
A(10,7) = 2.84783876419280044916451825422D0
A(10,8) = 2.81387331469849792539403641827D-1
A(10,9) = 1.23744899863314657627030212664D-1
A(11,1) = -1.21542481739588805916051052503D0
A(11,2) = 0.D0
A(11,3) = 0.D0
A(11,4) = 1.66726086659457724322804132886D1
A(11,5) = 9.15741828416817960595718650451D-1
A(11,6) = -6.05660580435747094755450554309D0
A(11,7) = -1.60035735941561781118417064101D1
A(11,8) = 1.4849303086297662557545391898D1
A(11,9) = -1.33715757352898493182930413962D1
A(11,10) = 5.13418264817963793317325361166D0
A(12,1) = 2.58860916438264283815730932232D-1
A(12,2) = 0.D0
A(12,3) = 0.D0
A(12,4) = -4.77448578548920511231011750971D0
A(12,5) = -4.3509301377703250944070041181D-1
A(12,6) = -3.04948333207224150956051286631D0
A(12,7) = 5.57792003993609911742367663447D0
A(12,8) = 6.15583158986104009733868912669D0
A(12,9) = -5.06210458673693837007740643391D0
A(12,10) = 2.19392617318067906127491429047D0
A(12,11) = 1.34627998659334941535726237887D-1
A(13,1) = 8.22427599626507477963168204773D-1
A(13,2) = 0.D0
A(13,3) = 0.D0
A(13,4) = -1.16586732572776642839765530355D1
A(13,5) = -7.57622116690936195881116154088D-1
A(13,6) = 7.13973588159581527978269282765D-1
A(13,7) = 1.20757749868900567395661704486D1
A(13,8) = -2.12765911392040265639082085897D0
A(13,9) = 1.99016620704895541832807169835D0
A(13,10) = -2.34286471544040292660294691857D-1
A(13,11) = 1.7589857770794226507310510589D-1
A(13,12) = 0.D0
C
C The coefficients BHAT(*) refer to the formula used to advance the
C integration, here the one of order 8. The coefficients B(*) refer
C to the other formula, here the one of order 7.
C
BHAT(1) = 4.17474911415302462220859284685D-2
BHAT(2) = 0.D0
BHAT(3) = 0.D0
BHAT(4) = 0.D0
BHAT(5) = 0.D0
BHAT(6) = -5.54523286112393089615218946547D-2
BHAT(7) = 2.39312807201180097046747354249D-1
BHAT(8) = 7.0351066940344302305804641089D-1
BHAT(9) = -7.59759613814460929884487677085D-1
BHAT(10) = 6.60563030922286341461378594838D-1
BHAT(11) = 1.58187482510123335529614838601D-1
BHAT(12) = -2.38109538752862804471863555306D-1
BHAT(13) = 2.5D-1
C
B(1) = 2.9553213676353496981964883112D-2
B(2) = 0.D0
B(3) = 0.D0
B(4) = 0.D0
B(5) = 0.D0
B(6) = -8.28606276487797039766805612689D-1
B(7) = 3.11240900051118327929913751627D-1
B(8) = 2.46734519059988698196468570407D0
B(9) = -2.54694165184190873912738007542D0
B(10) = 1.44354858367677524030187495069D0
B(11) = 7.94155958811272872713019541622D-2
B(12) = 4.44444444444444444444444444445D-2
B(13) = 0.D0
C
C(1) = 0.D0
C(2) = 5.55555555555555555555555555556D-2
C(3) = 8.33333333333333333333333333334D-2
C(4) = 1.25D-1
C(5) = 3.125D-1
C(6) = 3.75D-1
C(7) = 1.475D-1
C(8) = 4.65D-1
C(9) = 5.64865451382259575398358501426D-1
C(10) = 6.5D-1
C(11) = 9.24656277640504446745013574318D-1
C(12) = 1.D0
C(13) = C(12)
C
GO TO 120
C
C The definitions of all pairs come here for the calculation of
C LSTSTG, RS1, RS2, RS3, RS4, COST, MAXTRY, EXPON, TOOSML, and VECSTG.
C
120 CONTINUE
LSTSTG = PTR(NSTAGE)
IF (FSAL) THEN
COST = DBLE(NSTAGE-1)
ELSE
COST = DBLE(NSTAGE)
END IF
C
C MAXTRY - limit on the number of iterations of a computation made in
C diagnosing stiffness. There are at most Q = 3 function calls per
C iteration. MAXTRY is determined so that Q*MAXTRY <= 5% of the cost of
C 50 steps and 1 <= MAXTRY <= 8. This limits the number of calls to FCN
C in each diagnosis of stiffness to 24 calls.
C
MAXTRY = MIN(8,MAX(1,INT(FIVEPC*COST*FIFTY)))
C
EXPON = ONE/(ORDER+ONE)
C
C In calculating CDIFF it is assumed that there will be a non-zero
C difference |C(I) - C(J)| less than one. If C(I) = C(J) for any I not
C equal to J, they should be made precisely equal by assignment.
C
CDIFF = ONE
DO 160 I = 1, NSTAGE - 1
DO 140 J = I + 1, NSTAGE
DIFF = ABS(C(I)-C(J))
IF (DIFF.NE.ZERO) CDIFF = MIN(CDIFF,DIFF)
140 CONTINUE
160 CONTINUE
TOOSML = RNDOFF/CDIFF
C
C Determine the number of columns needed in STAGES(1:NEQ,*) (this will be
C at most NSTAGE-1 since the first stage is held in a separate array).
C The PTR array contains the column positions of the stages.
C
VECSTG = 0
DO 180 I = 2, NSTAGE
VECSTG = MAX(PTR(I),VECSTG)
180 CONTINUE
C
RS = TWO
RS1 = ONE/RS
RS2 = RS**2
RS3 = RS**3
RS4 = ONE/RS3
C
RETURN
END SUBROUTINE CONST
C
SUBROUTINE FORMI(F,NEQ,NWANT,Y,YP,YOLD,YPOLD,STAGES,CALSTG,
& XSTAGE,P)
C************************************************
C**** NOT A DESIGNATED USER-CALLABLE ROUTINE ****
C************************************************
C
C Purpose: Forms an interpolating polynomial for use with
C METHDs 1 or 2.
C
C Input: NEQ, NWANT, T, Y(*), YP(*), HOLD, YOLD(*), YPOLD(*),
C STAGES(NEQ,*), CALSTG
C Output: P(*), XSTAGE(NEQ)
C External: F
C
C Common: Initializes: none
C Reads: /RKCOM4/ A(*,*), C(*), R(*), METHD, MINTP
C /RKCOM2/ T, TOLD, HOLD
C Alters: /RKCOM2/ NFCN
C
C Comments:
C =========
C The integration has reached T with a step HOLD from TOLD = T-HOLD.
C Y(*),YP(*) and YOLD(*),YPOLD(*) approximate the solution and its
C derivative at T and TOLD respectively. STAGES(NEQ,*) holds the stages
C computed in taking this step. In the case of METHD = 2 it is necessary
C to compute some more stages in this subroutine. CALSTG indicates whether
C or not the extra stages need to be computed. A(*,*) and C(*) are used in
C computing these stages. The extra stages are stored in STAGES(NEQ,*) and
C XSTAGE(*). The coefficients of the interpolating polynomials for the first
C NWANT components of the solution are returned in the array P(*). The
C polynomial is of degree MINTP = 3 for METHD = 1 and of degree MINTP = 6
C for METHD = 2. The vector R(*) is used for workspace when METHD = 2.
C
C .. Scalar Arguments ..
INTEGER NEQ, NWANT
LOGICAL CALSTG
C .. Array Arguments ..
DOUBLE PRECISION P(*), STAGES(NEQ,*), XSTAGE(*), Y(*), YOLD(*),
& YP(*), YPOLD(*)
C .. Subroutine Arguments ..
EXTERNAL F
C .. Common Block to hold Problem Status ..
DOUBLE PRECISION T, H, TOLD, HOLD
INTEGER NFCN, SVNFCN, OKSTP, FLSTP
LOGICAL FIRST, LAST
COMMON /RKCOM2/ T, H, TOLD, HOLD, NFCN, SVNFCN, OKSTP, FLSTP,
& FIRST, LAST
SAVE /RKCOM2/
C .. Common Block to hold Formula Definitions ..
DOUBLE PRECISION A(13,13), B(13), C(13), BHAT(13), R(11,6),
& E(7)
INTEGER PTR(13), NSTAGE, METHD, MINTP
LOGICAL INTP
COMMON /RKCOM4/ A, B, C, BHAT, R, E, PTR, NSTAGE, METHD,
& MINTP, INTP
SAVE /RKCOM4/
C .. Local Scalars ..
DOUBLE PRECISION D1, D2, D3, D4, HYP, HYPOLD
INTEGER I, J, K, L
C .. Executable Statements ..
C
IF (METHD.EQ.1) THEN
C
C METHD = 1. Use the cubic Hermite interpolant that is is fully
C specified by the values and slopes at the two ends of the step.
C
DO 20 L = 1, NWANT
D1 = Y(L) - YOLD(L)
HYP = HOLD*YP(L)
HYPOLD = HOLD*YPOLD(L)
D2 = HYP - D1
D3 = D1 - HYPOLD
D4 = D2 - D3
P(L) = D2 + D4
P(NWANT+L) = D4
20 CONTINUE
C
ELSE
C
C METHD = 2.
C
IF (CALSTG) THEN
C
C Compute the extra stages needed for interpolation using the facts that
C 1. Stage 1 is YPOLD(*).
C 2. Stage i (i>1) is stored in STAGES(1:NEQ,i).
C 3. This pair is FSAL, i.e. STAGES(1:NEQ,7)=YP(1:NEQ), which frees
C up STAGES(1:NEQ,7) for use by stage 9.
C 4. XSTAGE(1:NEQ) is used for stage 10.
C 5. The coefficient of stage 2 in the interpolant is always 0, so
C STAGES(1:NEQ,1) is used for stage 11.
C The vector P(1:NEQ) is used as workspace for computing the stages.
C
DO 180 I = 9, 11
DO 40 L = 1, NEQ
P(L) = A(I,1)*YPOLD(L)
40 CONTINUE
DO 140 J = 2, I - 1
IF (J.LE.7) THEN
DO 60 L = 1, NEQ
P(L) = P(L) + A(I,J)*STAGES(L,J-1)
60 CONTINUE
ELSE IF (J.EQ.8) THEN
DO 80 L = 1, NEQ
P(L) = P(L) + A(I,J)*YP(L)
80 CONTINUE
ELSE IF (J.EQ.9) THEN
DO 100 L = 1, NEQ
P(L) = P(L) + A(I,J)*STAGES(L,7)
100 CONTINUE
ELSE IF (J.EQ.10) THEN
DO 120 L = 1, NEQ
P(L) = P(L) + A(I,J)*XSTAGE(L)
120 CONTINUE
END IF
140 CONTINUE
DO 160 L = 1, NEQ
P(L) = YOLD(L) + HOLD*P(L)
160 CONTINUE
IF (I.EQ.9) THEN
CALL F(TOLD+C(I)*HOLD,P,STAGES(1,7))
NFCN = NFCN + 1
ELSE IF (I.EQ.10) THEN
CALL F(TOLD+C(I)*HOLD,P,XSTAGE)
NFCN = NFCN + 1
ELSE
CALL F(TOLD+C(I)*HOLD,P,STAGES(1,1))
NFCN = NFCN + 1
END IF
180 CONTINUE
END IF
C
C Form the coefficients of the interpolating polynomial in its shifted
C and scaled form. The transformation from the form in which the
C polynomial is derived can be somewhat ill-conditioned. The terms
C are grouped so as to minimize the errors of the transformation.
C
C Coefficient of SIGMA**6
K = 4*NWANT
DO 200 L = 1, NWANT
P(K+L) = R(5,6)*STAGES(L,4) +
& ((R(10,6)*XSTAGE(L)+R(8,6)*YP(L))+
& (R(7,6)*STAGES(L,6)+R(6,6)*STAGES(L,5))) +
& ((R(4,6)*STAGES(L,3)+R(9,6)*STAGES(L,7))+
& (R(3,6)*STAGES(L,2)+R(11,6)*STAGES(L,1))+
& R(1,6)*YPOLD(L))
200 CONTINUE
C
C Coefficient of SIGMA**5
K = 3*NWANT
DO 220 L = 1, NWANT
P(K+L) = (R(10,5)*XSTAGE(L)+R(9,5)*STAGES(L,7)) +
& ((R(7,5)*STAGES(L,6)+R(6,5)*STAGES(L,5))+
& R(5,5)*STAGES(L,4)) + ((R(4,5)*STAGES(L,3)+
& R(8,5)*YP(L))+(R(3,5)*STAGES(L,2)+R(11,5)*
& STAGES(L,1))+R(1,5)*YPOLD(L))
220 CONTINUE
C
C Coefficient of SIGMA**4
K = 2*NWANT
DO 240 L = 1, NWANT
P(K+L) = ((R(4,4)*STAGES(L,3)+R(8,4)*YP(L))+
& (R(7,4)*STAGES(L,6)+R(6,4)*STAGES(L,5))+
& R(5,4)*STAGES(L,4)) + ((R(10,4)*XSTAGE(L)+
& R(9,4)*STAGES(L,7))+(R(3,4)*STAGES(L,2)+
& R(11,4)*STAGES(L,1))+R(1,4)*YPOLD(L))
240 CONTINUE
C
C Coefficient of SIGMA**3
K = NWANT
DO 260 L = 1, NWANT
P(K+L) = R(5,3)*STAGES(L,4) + R(6,3)*STAGES(L,5) +
& ((R(3,3)*STAGES(L,2)+R(9,3)*STAGES(L,7))+
& (R(10,3)*XSTAGE(L)+R(8,3)*YP(L))+R(1,3)*
& YPOLD(L))+((R(4,3)*STAGES(L,3)+R(11,3)*
& STAGES(L,1))+R(7,3)*STAGES(L,6))
260 CONTINUE
C
C Coefficient of SIGMA**2
C
DO 280 L = 1, NWANT
P(L) = R(5,2)*STAGES(L,4) + ((R(6,2)*STAGES(L,5)+
& R(8,2)*YP(L))+R(1,2)*YPOLD(L)) +
& ((R(3,2)*STAGES(L,2)+R(9,2)*STAGES(L,7))+
& R(10,2)*XSTAGE(L)) + ((R(4,2)*STAGES(L,3)+
& R(11,2)*STAGES(L,1))+R(7,2)*STAGES(L,6))
280 CONTINUE
C
C Scale all the coefficients by the step size.
C
DO 300 L = 1, NWANT*(MINTP-1)
P(L) = HOLD*P(L)
300 CONTINUE
C
END IF
C
RETURN
END SUBROUTINE FORMI
C
SUBROUTINE EVALI(Y,YP,P,TWANT,REQEST,NWANT,YWANT,YPWANT)
C************************************************
C**** NOT A DESIGNATED USER-CALLABLE ROUTINE ****
*************************************************
C
C Purpose: Evaluation of an interpolating polynomial and/or its
C first derivative.
C
C Input: Y(*), YP(*), P(NWANT,*), TWANT, REQEST, NWANT
C Output: YWANT(*), YPWANT(*)
C
C Common: Initializes: none
C Reads: /RKCOM2/ HOLD, T
C /RKCOM4/ MINTP
C Alters: none
C
C Comments:
C =========
C The interpolant is evaluated at TWANT to approximate the solution,
C YWANT, and/or its first derivative there, YPWANT. Only the first
C NWANT components of the answer are computed. There are three cases
C that are indicated by the first character of REQEST:
C REQEST(1:1) = `S' or `s'- compute approximate `S'olution only.
C = `D' or `d'- compute approximate first `D'erivative
C of the solution only.
C = `B' or `b'- compute `B'oth approximate solution and
C first derivative.
C The coefficents of the polynomial are contained in Y(*), YP(*) and
C P(NWANT,*).
C
C .. Scalar Arguments ..
DOUBLE PRECISION TWANT
INTEGER NWANT
CHARACTER*(*) REQEST
C .. Array Arguments ..
DOUBLE PRECISION P(NWANT,*), Y(*), YP(*), YPWANT(*), YWANT(*)
C .. Common Block to hold Problem Status ..
DOUBLE PRECISION T, H, TOLD, HOLD
INTEGER NFCN, SVNFCN, OKSTP, FLSTP
LOGICAL FIRST, LAST
COMMON /RKCOM2/ T, H, TOLD, HOLD, NFCN, SVNFCN, OKSTP, FLSTP,
& FIRST, LAST
SAVE /RKCOM2/
C .. Common Block to hold Formula Definitions ..
DOUBLE PRECISION A(13,13), B(13), C(13), BHAT(13), R(11,6),
& E(7)
INTEGER PTR(13), NSTAGE, METHD, MINTP
LOGICAL INTP
COMMON /RKCOM4/ A, B, C, BHAT, R, E, PTR, NSTAGE, METHD,
& MINTP, INTP
SAVE /RKCOM4/
C .. Local Scalars ..
DOUBLE PRECISION SIGMA
INTEGER K, L
CHARACTER REQST1
C .. Executable Statements ..
C
C Evaluate the interpolating polynomial of degree MINTP in terms of the
C shifted and scaled independent variable SIGMA.
C
SIGMA = (TWANT-T)/HOLD
C
REQST1 = REQEST(1:1)
IF (REQST1.EQ.'S' .OR. REQST1.EQ.'s' .OR.
& REQST1.EQ.'B' .OR. REQST1.EQ.'b') THEN
C
DO 20 L = 1, NWANT
YWANT(L) = P(L,MINTP-1)*SIGMA
20 CONTINUE
DO 60 K = MINTP - 2, 1, -1
DO 40 L = 1, NWANT
YWANT(L) = (YWANT(L)+P(L,K))*SIGMA
40 CONTINUE
60 CONTINUE
DO 80 L = 1, NWANT
YWANT(L) = (YWANT(L)+HOLD*YP(L))*SIGMA + Y(L)
80 CONTINUE
END IF
C
C Evaluate the derivative of the interpolating polynomial.
C
IF (REQST1.EQ.'D' .OR. REQST1.EQ.'d' .OR.
& REQST1.EQ.'B' .OR. REQST1.EQ.'b') THEN
C
C The derivative of the interpolating polynomial with respect to TWANT
C is the derivative with respect to S divided by HOLD.
C
DO 100 L = 1, NWANT
YPWANT(L) = MINTP*P(L,MINTP-1)*SIGMA
100 CONTINUE
DO 140 K = MINTP - 1, 2, -1
DO 120 L = 1, NWANT
YPWANT(L) = (YPWANT(L)+K*P(L,K-1))*SIGMA
120 CONTINUE
140 CONTINUE
DO 160 L = 1, NWANT
YPWANT(L) = (YPWANT(L)+HOLD*YP(L))/HOLD
160 CONTINUE
END IF
C
RETURN
END SUBROUTINE EVALI
C
SUBROUTINE RKMSG(IER,SRNAME,NREC,FLAG)
C************************************************
C**** NOT A DESIGNATED USER-CALLABLE ROUTINE ****
C************************************************
C
C Purpose: To process error messages and terminate the program
C in the event of a "catastrophic" failure.
C
C Input: IER, SRNAME, NREC
C Output: FLAG
C
C Common: Initializes: none
C Reads: /RKCOM7/ OUTCH
C /RKCOM8/ MSG, UTASK
C /RKCOM9/ REC
C Alters: none
C
C Comments:
C =========
C The output variable FLAG is assigned the value of the input variable IER.
C
C IER = -2 reports a successful call of the subroutine SRNAME and
C indicates that special action is to be taken elsewhere
C in the suite. FLAG is set and a return is effected.
C
C IER = 1 reports a successful call of the subroutine SRNAME. FLAG
C is set and a return is effected.
C
C 1 < IER < 911 and MSG = .TRUE.: a message of NREC records contained in
C the array REC(*) is written to the standard output channel,
C OUTCH. FLAG is set and a return is effected.
C
C IER = 911 reports a "catastrophic" error was detected in SRNAME. A
C message is written to OUTCH regardless of the value of MSG and
C normally the execution of the program is terminated. The
C execution is not terminated if the error is the result of an
C indirect call to CT, RESET, or INTRP through UT (UTASK = .TRUE.).
C Termination can be prevented by using the subroutine SOFTFL.
C
C IER = 912 reports that a "catastrophic" error was detected earlier and
C termination was prevented, but the user has failed to take
C appropriate remedial action. Execution is terminated.
C
C .. Scalar Arguments ..
INTEGER FLAG, IER, NREC
CHARACTER*(*) SRNAME
C .. Common Block for Environment Parameters ..
DOUBLE PRECISION MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY
INTEGER OUTCH
COMMON /RKCOM7/ MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY,
& OUTCH
SAVE /RKCOM7/
C .. Common Block for Integrator Options ..
LOGICAL MSG, UTASK
COMMON /RKCOM8/ MSG, UTASK
SAVE /RKCOM8/
C .. Common Block for Error Message ..
CHARACTER*80 REC(10)
COMMON /RKCOM9/ REC
SAVE /RKCOM9/
C .. Parameters ..
INTEGER PLUS1
LOGICAL ASK, TELL
PARAMETER (PLUS1=1,ASK=.TRUE.,TELL=.FALSE.)
C .. Local Scalars ..
INTEGER I
LOGICAL BADERR, OK, ON, UTCALL
C .. External Subroutines ..
EXTERNAL CHKFL, RKSIT, SOFTFL
C .. Executable Statements ..
C
C Check where the call came from - if it is an indirect call from UT,
C the run is not STOPped.
UTCALL = (SRNAME.EQ.'RESET' .OR. SRNAME.EQ.'CT' .OR.
& SRNAME.EQ.'INTRP') .AND. UTASK
C
C Check if can continue with integrator.
OK = (SRNAME.EQ.'CT' .OR. SRNAME.EQ.'UT') .AND.
& (IER.EQ.2 .OR. IER.EQ.3 .OR. IER.EQ.4)
C
C Check if program termination has been overridden.
CALL SOFTFL(ASK,ON)
C
IF ((MSG.AND.IER.GT.PLUS1) .OR. IER.GE.911) THEN
WRITE (OUTCH,'(/A)') ' **'
WRITE (OUTCH,'(A)') (REC(I),I=1,NREC)
IF (IER.GE.911) THEN
WRITE (OUTCH,'(A/A,A,A/A/)')
&' **',
&' ** Catastrophic error detected in ', SRNAME, '.',
&' **'
IF ((.NOT.(UTCALL.OR.ON).AND.IER.EQ.911) .OR.
& IER.EQ.912) THEN
WRITE (OUTCH,'(A/A/A)')
&' **',
&' ** Execution of your program is being terminated.',
&' **'
STOP
END IF
ELSE IF (OK) THEN
WRITE (OUTCH,'(A/A,A,A,I2,A/A/A)')
&' **',
&' ** Warning from routine ', SRNAME, ' with flag set ',IER, '.',
&' ** You can continue integrating this problem.',
&' **'
ELSE
WRITE (OUTCH,'(A/A,A,A,I2,A/A/A)')
&' **',
&' ** Warning from routine ', SRNAME, ' with flag set ',IER, '.',
&' ** You cannot continue integrating this problem.',
&' **'
END IF
END IF
DO 20 I = NREC + 1, 10
REC(I) = ' '
20 CONTINUE
FLAG = IER
C
C TELL RKSIT the status of the routine associated with SRNAME
CALL RKSIT(TELL,SRNAME,FLAG)
C
C Indicate that a catastrophic error has been detected
BADERR = FLAG .GE. 911
CALL CHKFL(TELL,BADERR)
C
RETURN
C
END SUBROUTINE RKMSG
C
SUBROUTINE RKSIT(ASK,SRNAME,STATE)
C************************************************
C**** NOT A DESIGNATED USER-CALLABLE ROUTINE ****
C************************************************
C
C Purpose: To save or enquire about the status of each
C subprogram in the suite.
C
C Input: ASK, SRNAME
C Input/output: STATE
C
C
C Comments:
C =========
C SRNAME indicates which routine is of interest in the call to RKSIT.
C
C If ASK=.FALSE., then the value of STATE (which as used is the error
C flag) for the routine SRNAME is saved internally. This value of STATE
C is usually positive. There are two exceptions:
C 1. SRNAME='SETUP' and STATE=-1 indicates a completely new problem,
C so all SAVEd states are cleared.
C 2. STATE=-2 is used by some routines in the suite to indicate
C circumstances which require special action.
C
C If ASK=.TRUE., then RKSIT first checks to see if there were any
C catastrophic errors, that is, a SAVEd state has value 911. This should
C happen only when the user has overridden program termination in the event
C of catastrophic failures from routines in the package but has failed to
C take appropriate action. If this is the case, then RKSIT returns a value
C of STATE = 911 which forces a termination of execution inside RKMSG. If
C no catastrophic errors are flagged, then STATE returns the saved state
C value for the routine specified by SRNAME.
C
C .. Scalar Arguments ..
INTEGER STATE
LOGICAL ASK
CHARACTER*(*) SRNAME
C .. Parameters ..
INTEGER STATES, MINUS1
PARAMETER (STATES=7,MINUS1=-1)
C .. Local Scalars ..
INTEGER I, NAME
C .. Local Arrays ..
INTEGER SVSTA(STATES)
C .. Save statement ..
SAVE SVSTA
C .. Data statements ..
DATA SVSTA/STATES*MINUS1/
C .. Executable Statements ..
C
IF (SRNAME.EQ.'SETUP') THEN
NAME = 1
ELSE IF (SRNAME.EQ.'UT') THEN
NAME = 2
ELSE IF (SRNAME.EQ.'STAT') THEN
NAME = 3
ELSE IF (SRNAME.EQ.'GLBERR') THEN
NAME = 4
ELSE IF (SRNAME.EQ.'CT') THEN
NAME = 5
ELSE IF (SRNAME.EQ.'INTRP') THEN
NAME = 6
ELSE IF (SRNAME.EQ.'RESET') THEN
NAME = 7
ELSE
NAME = 0
END IF
C
C (Re)initialize if SETUP is telling RKSIT to do so.
IF (.NOT.ASK .AND. NAME.EQ.1 .AND. STATE.EQ.MINUS1) THEN
DO 20 I = 1, STATES
SVSTA(I) = MINUS1
20 CONTINUE
GO TO 60
END IF
C
C Check for 911 on exit from a previous call.
IF (ASK) THEN
DO 40 I = 1, STATES
IF (SVSTA(I).EQ.911) THEN
STATE = 911
GO TO 60
END IF
40 CONTINUE
END IF
C
IF (ASK) THEN
STATE = SVSTA(NAME)
ELSE
SVSTA(NAME) = STATE
END IF
C
60 CONTINUE
C
RETURN
END SUBROUTINE RKSIT
C
SUBROUTINE TRUERR(F,NEQ,Y,TOL,WEIGHT,ZY,ZYP,ZERROR,ZYNEW,ZERRES,
& ZSTAGE,IER)
C************************************************
C**** NOT A DESIGNATED USER-CALLABLE ROUTINE ****
C************************************************
C
C Purpose: Compute a running RMS measure of the true (global) error
C for a general Runge-Kutta pair.
C
C
C Input: NEQ, Y(*), TOL, WEIGHT(*),
C Input/output: ZY(*), ZYP(*), ZERROR(*)
C Workspace: ZYNEW(*), ZERRES(*), ZSTAGE(NEQ,*)
C Output: IER
C External: F
C
C Common: Initializes: none
C Reads: /RKCOM2/ T, HOLD
C /RKCOM5/ TOOSML, ORDER, NSEC
C /RKCOM7/ TINY
C Alters: /RKCOM6/ MAXERR, LOCMAX, GNFCN
C
C Comments:
C =========
C A secondary integration is performed using a fraction of the step size
C of the primary integration. ZY(*) and ZYP(*) are the approximate solution
C and first derivative of this secondary integration. ZERRES(*) contains the
C error estimates for the secondary integration. ZYNEW(*) and ZSTAGE(*,*) are
C workspace for taking a step. The error assessment is computed using the
C difference of the primary and secondary solutions at the primary
C integration points as an estimate of the true error there. The weights
C used are those of the error test of the primary integration. This error
C assessment is maintained in the vector ZERROR(*). MAXERR and LOCMAX
C contain the maximum contribution to the assessment and its location,
C respectively. The number of calls to F is counted by GNFCN.
C
C .. Scalar Arguments ..
DOUBLE PRECISION TOL
INTEGER IER, NEQ
C .. Array Arguments ..
DOUBLE PRECISION WEIGHT(*), Y(*), ZERRES(*), ZERROR(*),
& ZSTAGE(NEQ,*), ZY(*), ZYNEW(*), ZYP(*)
C .. Subroutine Arguments ..
EXTERNAL F
C .. Common Block to hold Problem Status ..
DOUBLE PRECISION T, H, TOLD, HOLD
INTEGER NFCN, SVNFCN, OKSTP, FLSTP
LOGICAL FIRST, LAST
COMMON /RKCOM2/ T, H, TOLD, HOLD, NFCN, SVNFCN, OKSTP, FLSTP,
& FIRST, LAST
SAVE /RKCOM2/
C .. Common Block to hold Formula Characterisitcs ..
DOUBLE PRECISION TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
& RS, RS1, RS2, RS3, RS4
INTEGER ORDER, LSTSTG, MAXTRY, NSEC
LOGICAL FSAL
COMMON /RKCOM5/ TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
& RS, RS1, RS2, RS3, RS4, ORDER, LSTSTG, MAXTRY,
& NSEC, FSAL
SAVE /RKCOM5/
C .. Common Block for Global Error Assessment ..
DOUBLE PRECISION MAXERR, LOCMAX
INTEGER GNFCN, PRZSTG, PRZY, PRZYP, PRZERS, PRZERR,
& PRZYNU
LOGICAL ERASON, ERASFL
COMMON /RKCOM6/ MAXERR, LOCMAX, GNFCN, PRZSTG, PRZY, PRZYP,
& PRZERS, PRZERR, PRZYNU, ERASON, ERASFL
SAVE /RKCOM6/
C .. Common Block for Environment Parameters ..
DOUBLE PRECISION MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY
INTEGER OUTCH
COMMON /RKCOM7/ MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY,
& OUTCH
SAVE /RKCOM7/
C .. Parameters ..
DOUBLE PRECISION PT1, TEN, DUMMY
PARAMETER (PT1=0.1D0,TEN=10.0D0,DUMMY=1.0D0)
C .. Local Scalars ..
DOUBLE PRECISION DIFF, ERRMAX, HMIN, HSEC, MXERLC, TSEC, ZLERR,
& ZTEST1, ZTEST2
INTEGER ISTEP, L, LEVEL
LOGICAL LDUMMY, MAIN
C .. Local Arrays ..
DOUBLE PRECISION DUMARR(1)
C .. External Subroutines ..
EXTERNAL STEP
C .. Intrinsic Functions ..
INTRINSIC ABS, DBLE, MAX
C .. Executable Statements ..
TSEC = T - HOLD
HSEC = HOLD/DBLE(NSEC)
HMIN = MAX(TINY,TOOSML*MAX(ABS(TSEC),ABS(T)))
IF (ABS(HSEC).LT.HMIN) THEN
IER = 6
GO TO 120
END IF
ZTEST1 = TOL/DBLE(NSEC)
ZTEST2 = TOL/TEN
LEVEL = 0
C
C The subroutine STEP is used to take a step. In its use in the primary
C integration provision is made for getting on scale in the first step.
C In this situation only the subroutine might reduce the step size. By
C setting MAIN = .FALSE., the subroutine will take a step of the size input.
C In this use of the subroutine, all items of the call list appearing after
C MAIN are dummy variables.
C
C Perform secondary integration.
MAIN = .FALSE.
LDUMMY = .FALSE.
DO 60 ISTEP = 1, NSEC
C
C Take a step.
CALL STEP(F,NEQ,TSEC,ZY,ZYP,ZSTAGE,ZTEST1,HSEC,WEIGHT,ZYNEW,
& ZERRES,ZLERR,MAIN,DUMMY,DUMARR,LDUMMY)
C
C The primary integration is using a step size of HUSED and the secondary
C integration is using the smaller step size HSEC = HUSED/NSEC. If steps
C of this size were taken from the same starting point and the asymptotic
C behavior were evident, the smaller step size would result in a local error
C that is considerably smaller, namely by a factor of 1/(NSEC**(ORDER+1)).
C If the two approximate solutions are close and TOLR is neither too large nor
C too small, this should be approximately true. The step size is chosen in
C the primary integration so that the local error ERR is no larger than TOLR.
C The local error, ZLERR, of the secondary integration is compared to TOLR in
C an attempt to diagnose a secondary integration that is not rather more
C accurate than the primary integration.
C
IF (ZLERR.GE.ZTEST1) THEN
LEVEL = 2
ELSE IF (ZLERR.GT.ZTEST2) THEN
LEVEL = LEVEL + 1
END IF
IF (LEVEL.GE.2) THEN
IER = 6
GO TO 120
END IF
C
C Advance TSEC and the dependent variables ZY(*) and ZYP(*).
TSEC = T - DBLE(NSEC-ISTEP)*HSEC
DO 20 L = 1, NEQ
ZY(L) = ZYNEW(L)
20 CONTINUE
C
IF (FSAL) THEN
C
C When FSAL = .TRUE., the derivative ZYP(*) is the last stage of the step.
DO 40 L = 1, NEQ
ZYP(L) = ZSTAGE(L,LSTSTG)
40 CONTINUE
ELSE
C
C Call F to evaluate ZYP(*).
CALL F(TSEC,ZY,ZYP)
GNFCN = GNFCN + 1
END IF
C
60 CONTINUE
C
C Update the maximum error seen, MAXERR, and its location, LOCMAX.
C Use local variables ERRMAX and MXERLC.
C
ERRMAX = MAXERR
MXERLC = LOCMAX
DO 80 L = 1, NEQ
DIFF = ABS(ZY(L)-Y(L))/WEIGHT(L)
IF (DIFF.GT.ERRMAX) THEN
ERRMAX = DIFF
MXERLC = T
END IF
80 CONTINUE
C
C If the global error is greater than 0.1D0, the solutions have diverged so
C far that comparing them may not provide a reliable estimate of the global
C error. The test is made before ZERROR(*) and MAXERR, LCMXER are updated so
C that on a failure, they refer to the last reliable results.
C
IF (ERRMAX.GT.PT1) THEN
IER = 6
GO TO 120
ELSE
MAXERR = ERRMAX
LOCMAX = MXERLC
DO 100 L = 1, NEQ
DIFF = ABS(ZY(L)-Y(L))/WEIGHT(L)
ZERROR(L) = ZERROR(L) + DIFF**2
100 CONTINUE
IER = 1
END IF
C
C Exit point for TRUERR
120 CONTINUE
C
RETURN
END SUBROUTINE TRUERR
C
SUBROUTINE STEP(F,NEQ,TNOW,Y,YP,STAGES,TOL,HTRY,WEIGHT,YNEW,
& ERREST,ERR,MAIN,HMIN,THRES,PHASE2)
C************************************************
C**** NOT A DESIGNATED USER-CALLABLE ROUTINE ****
C************************************************
C
C Purpose: To compute a step of an explicit Runge-Kutta
C method and estimate the local error of the step.
C
C Input: NEQ, TNOW, Y(*), YP(*), TOL, MAIN, HMIN, THRES(*)
C Input/output: HTRY, PHASE2, LAST, WEIGHT(*)
C Output: STAGES(NEQ,*), YNEW(*), ERREST(*), ERR
C
C Common: Initializes: none
C Reads: /RKCOM1/ TND
C /RKCOM2/ LAST
C /RKCOM4/ A, B, C, BHAT, PTR, NSTAGE, METHD
C /RKCOM5/ FSAL
C Alters: /RKCOM2/ NFCN, LAST
C /RKCOM6/ GNFCN
C
C Comments:
C =========
C From an approximate solution Y(*) at TNOW and first derivative there,
C YP(*) = F(TNOW,Y,YP), a step is taken to get an approximation YNEW(*)
C at TNOW + HTRY. The Runge-Kutta method and how it is used are defined
C by A, B, C, BHAT, PTR, NSTAGE, METHD and FSAL. Intermediate stages
C of the method are stored in the array STAGES(NEQ,*). The error in
C each solution component is estimated and returned in ERREST(*). A
C weighted maximum norm of the local error, ERR, is formed. For some
C methods an intermediate error estimate can be computed before completion
C of the step (see routine STEPB); if the estimate is greater than the
C specified tolerance TOL, the computation of the step is terminated.
C
C When global error estimation is desired, two integrations are done.
C The usual integration is referred to as the "primary", or "main",
C integration (MAIN=.TRUE.). For global error estimation another,
C "secondary" integration (MAIN=.FALSE.) is carried out with a smaller
C step size. The weight vector WEIGHT(*) used in computing ERR is
C determined by the main integration. Thus this argument is output when
C MAIN = .TRUE. and input when MAIN = .FALSE..
C
C When taking the first step in an integration, the logical variable
C PHASE2 may be input as .TRUE. and if the first step is the whole of
C the range of integration, then LAST will be .TRUE.. When PHASE2=.TRUE.,
C the first three stages are monitored to help assure that the step
C size H is small enough for the integration to be stable and for the
C estimate of the error of the step to be credible. Calls are made to
C the subroutine STEPA for this purpose. If necessary, H will be
C reduced in STEPA (and LAST altered accordingly) and the step retried
C in STEP until an acceptable value is found.
C
C In the primary integration the number of calls to F is counted by
C NFCN, and in the secondary integration, by GNFCN.
C
C .. Scalar Arguments ..
DOUBLE PRECISION ERR, HMIN, HTRY, TNOW, TOL
INTEGER NEQ
LOGICAL MAIN, PHASE2
C .. Array Arguments ..
DOUBLE PRECISION ERREST(*), STAGES(NEQ,*), THRES(*), WEIGHT(*),
& Y(*), YNEW(*), YP(*)
C .. Subroutine Arguments ..
EXTERNAL F
C .. Common Block for Problem Definition ..
DOUBLE PRECISION TSTRT, TND, DIR, HSTRT, TOLR
INTEGER NEQN
COMMON /RKCOM1/ TSTRT, TND, DIR, HSTRT, TOLR, NEQN
SAVE /RKCOM1/
C .. Common Block to hold Problem Status ..
DOUBLE PRECISION T, H, TOLD, HOLD
INTEGER NFCN, SVNFCN, OKSTP, FLSTP
LOGICAL FIRST, LAST
COMMON /RKCOM2/ T, H, TOLD, HOLD, NFCN, SVNFCN, OKSTP, FLSTP,
& FIRST, LAST
SAVE /RKCOM2/
C .. Common Block to hold Formula Definitions ..
DOUBLE PRECISION A(13,13), B(13), C(13), BHAT(13), R(11,6),
& E(7)
INTEGER PTR(13), NSTAGE, METHD, MINTP
LOGICAL INTP
COMMON /RKCOM4/ A, B, C, BHAT, R, E, PTR, NSTAGE, METHD,
& MINTP, INTP
SAVE /RKCOM4/
C .. Common Block to hold Formula Characterisitcs ..
DOUBLE PRECISION TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
& RS, RS1, RS2, RS3, RS4
INTEGER ORDER, LSTSTG, MAXTRY, NSEC
LOGICAL FSAL
COMMON /RKCOM5/ TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
& RS, RS1, RS2, RS3, RS4, ORDER, LSTSTG, MAXTRY,
& NSEC, FSAL
SAVE /RKCOM5/
C .. Common Block for Global Error Assessment ..
DOUBLE PRECISION MAXERR, LOCMAX
INTEGER GNFCN, PRZSTG, PRZY, PRZYP, PRZERS, PRZERR,
& PRZYNU
LOGICAL ERASON, ERASFL
COMMON /RKCOM6/ MAXERR, LOCMAX, GNFCN, PRZSTG, PRZY, PRZYP,
& PRZERS, PRZERR, PRZYNU, ERASON, ERASFL
SAVE /RKCOM6/
C .. Parameters ..
DOUBLE PRECISION ZERO, HALF, ONE
PARAMETER (ZERO=0.0D0,HALF=0.5D0,ONE=1.0D0)
C .. Local Scalars ..
DOUBLE PRECISION AVGY, TSTG
INTEGER I, J, L
LOGICAL CUTBAK
C .. External Subroutines ..
EXTERNAL STEPA, STEPB
C .. Intrinsic Functions ..
INTRINSIC ABS, MAX, SIGN
C .. Executable Statements ..
C
C Many of the following loops over L = 1, NEQ have constant array values
C inside. The code is written with clarity in mind. Any optimizing
C compiler will identify these occurrences and take appropriate action.
C A check for zero multipliers has been included so as to prevent
C needless computation resulting from the storing of zero coefficients
C in the arrays for the sake of clarity. The array ERREST(*) is used
C for working storage in this computation.
C
20 CONTINUE
IF (MAIN) THEN
IF (PHASE2) THEN
C
C Initialize weights for measuring the local error.
DO 40 L = 1, NEQ
WEIGHT(L) = MAX(THRES(L),ABS(Y(L)))
40 CONTINUE
END IF
END IF
C
DO 140 I = 2, NSTAGE
DO 100 J = 1, I - 1
IF (J.EQ.1) THEN
DO 60 L = 1, NEQ
ERREST(L) = A(I,1)*YP(L)
60 CONTINUE
ELSE
IF (A(I,J).NE.ZERO) THEN
DO 80 L = 1, NEQ
ERREST(L) = ERREST(L) + A(I,J)*STAGES(L,PTR(J))
80 CONTINUE
END IF
END IF
100 CONTINUE
DO 120 L = 1, NEQ
YNEW(L) = Y(L) + HTRY*ERREST(L)
120 CONTINUE
C
C METHD = 2 is special in that an estimate of the local error can be
C formed before the step is completed. If the step is a failure,
C return immediately. Otherwise, complete the step and compute a more
C accurate error estimate.
IF (METHD.EQ.2 .AND. I.EQ.7) THEN
CALL STEPB(NEQ,Y,YP,HTRY,YNEW,STAGES,THRES,ERR,MAIN,WEIGHT)
IF (ERR.GT.TOL) RETURN
END IF
C
TSTG = TNOW + C(I)*HTRY
IF (MAIN .AND. LAST .AND. C(I).EQ.ONE) TSTG = TND
CALL F(TSTG,YNEW,STAGES(1,PTR(I)))
C
C Increment the counter for the number of function evaluations
C depending on whether the primary or secondary integration is taking
C place.
IF (MAIN) THEN
NFCN = NFCN + 1
ELSE
GNFCN = GNFCN + 1
END IF
C
C----------------------------------------------------------------------
C When PHASE2 is .TRUE. we are in the second phase of the automatic
C selection of the initial step size. The results of the first three
C stages are monitored in the subroutine STEPA for evidence that H is
C too large -- instability and/or an unreliable estimate of the error
C of the step is then possible. When the subroutine believes H to be
C too large, it returns CUTBAK = .TRUE. and a suitably reduced H for
C another try.
C
IF (MAIN) THEN
IF (PHASE2) THEN
IF (I.LE.3 .AND. ABS(HTRY).GT.HMIN) THEN
CALL STEPA(TNOW,Y,YP,TSTG,YNEW,STAGES(1,PTR(I)),
& HTRY,WEIGHT,CUTBAK)
IF (CUTBAK) THEN
LAST = .FALSE.
C
C Make sure that STEPA does not reduce the step size below the
C minimum. If it does, reset H to HMIN and deactivate PHASE2.
IF (ABS(HTRY).LE.HMIN) THEN
HTRY = SIGN(HMIN,HTRY)
PHASE2 = .FALSE.
END IF
GO TO 20
END IF
END IF
END IF
END IF
C----------------------------------------------------------------------
C
140 CONTINUE
C
C Some formulas are constructed so that the last stage represents
C the result of the step (FSAL=.TRUE.), hence if the step is acceptable,
C it will be the first stage for the next step. When FSAL=.FALSE., we
C have to complete the computation of the step.
C
IF (.NOT.FSAL) THEN
DO 200 I = 1, NSTAGE
IF (I.EQ.1) THEN
DO 160 L = 1, NEQ
ERREST(L) = BHAT(1)*YP(L)
160 CONTINUE
ELSE
IF (BHAT(I).NE.ZERO) THEN
DO 180 L = 1, NEQ
ERREST(L) = ERREST(L) + BHAT(I)*STAGES(L,PTR(I))
180 CONTINUE
END IF
END IF
200 CONTINUE
DO 220 L = 1, NEQ
YNEW(L) = Y(L) + HTRY*ERREST(L)
220 CONTINUE
END IF
C
C Form an estimate of the error in the lower order formula by comparing
C it to the higher order formula of the pair. ERREST(*) has been used
C as working storage above. The higher order approximation has been
C formed as YNEW(*) = Y(*) + HTRY*ERREST(*) where ERREST(*) is a linear
C combination of the stages of the formula. The lower order result also
C has the form Y(*) plus HTRY times a different linear combination of
C the stages. Hence, this different linear combination of stages for
C the lower order formula can just be subtracted from the combination
C stored in ERREST(*) to produce the errors. The result is then
C multiplied by HTRY to obtain the error estimate.
C
DO 280 I = 1, NSTAGE
IF (I.EQ.1 .AND. B(1).NE.ZERO) THEN
DO 240 L = 1, NEQ
ERREST(L) = ERREST(L) - B(1)*YP(L)
240 CONTINUE
ELSE
IF (B(I).NE.ZERO) THEN
DO 260 L = 1, NEQ
ERREST(L) = ERREST(L) - B(I)*STAGES(L,PTR(I))
260 CONTINUE
END IF
END IF
280 CONTINUE
DO 300 L = 1, NEQ
ERREST(L) = HTRY*ERREST(L)
300 CONTINUE
C
C The error in a solution component is measured relative to a weight
C that is the larger of a threshold and the size of the solution over
C the step. Using the magnitude of a solution component at both ends
C of the step in the definition of "size" increases the robustness of
C the test. When global error estimation is specified, the weight
C vector WEIGHT(*) is defined by the primary integration and is then
C used in the secondary integration.
C
IF (MAIN) THEN
DO 320 L = 1, NEQ
AVGY = HALF*(ABS(Y(L))+ABS(YNEW(L)))
WEIGHT(L) = MAX(AVGY,THRES(L))
320 CONTINUE
END IF
C
ERR = ZERO
DO 340 L = 1, NEQ
ERR = MAX(ERR,ABS(ERREST(L)/WEIGHT(L)))
340 CONTINUE
C
RETURN
END SUBROUTINE STEP
C
SUBROUTINE STEPA(TNOW,Y,YP,TSTG,YSTG,YPSTG,HTRY,WEIGHT,CUTBAK)
C************************************************
C**** NOT A DESIGNATED USER-CALLABLE ROUTINE ****
C************************************************
C
C Purpose: To calculate an "on-scale" step size for phase 2 of
C the initial step size computation.
C
C Input: TNOW, Y(*), YP(*), TSTG, YSTG(*), YPSTG(*)
C Input/output: HTRY, WEIGHT
C Output: CUTBAK
C
C Common: Initializes: none
C Reads: /RKCOM1/ TND, NEQ
C /RKCOM5/ STBRAD, RS1, RS4
C /RKCOM7/ RNDOFF
C Alters: none
C
C Comments:
C =========
C This subroutine is used during the first three stages of the first step.
C A Lipschitz constant L for the differential equation in autonomous form
C is approximated, and the product abs(HTRY)*L is compared to an approximate
C radius, STBRAD, of the stability region of the method. The step size is
C reduced as necessary, within a range specified by the step size control
C parameters RS1 and RS4, to assure stability and give some confidence in
C the error estimator. If HTRY is reduced, CUTBAK is set .TRUE..
C
C Y(*) and YP(*) contain the solution and its derivative at TNOW and
C similarly YSTG(*) and YPSTG(*) contain approximations at TSTG.
C
C Normally the weights used in the control of the error depend on the
C size of the solution at the beginning and at the end of the step, but
C at this time we do not have a solution at the end of the step. Each
C stage YSTG(*) of the Runge - Kutta process represents a low order
C approximation to the solution at TSTG. Because the initial value of
C WEIGHT(*) provided in the first phase of the scheme is based only on
C the solution at T and THRES(*), it is continually updated in STEPA to
C account for the size of the solution throughout the step as revealed
C by the intermediate stages YSTG(*). Inside this subroutine only, the
C differential equation is converted to autonomous form. After the
C conversion, the end of the interval of integration, TND, is used
C to define a suitable weight for the independent variable.
C
C .. Scalar Arguments ..
DOUBLE PRECISION HTRY, TNOW, TSTG
LOGICAL CUTBAK
C .. Array Arguments ..
DOUBLE PRECISION WEIGHT(*), Y(*), YP(*), YPSTG(*), YSTG(*)
C .. Common Block for Problem Definition ..
DOUBLE PRECISION TSTRT, TND, DIR, HSTRT, TOLR
INTEGER NEQN
COMMON /RKCOM1/ TSTRT, TND, DIR, HSTRT, TOLR, NEQN
SAVE /RKCOM1/
C .. Common Block to hold Formula Characterisitcs ..
DOUBLE PRECISION TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
& RS, RS1, RS2, RS3, RS4
INTEGER ORDER, LSTSTG, MAXTRY, NSEC
LOGICAL FSAL
COMMON /RKCOM5/ TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
& RS, RS1, RS2, RS3, RS4, ORDER, LSTSTG, MAXTRY,
& NSEC, FSAL
SAVE /RKCOM5/
C .. Common Block for Environment Parameters ..
DOUBLE PRECISION MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY
INTEGER OUTCH
COMMON /RKCOM7/ MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY,
& OUTCH
SAVE /RKCOM7/
C .. Parameters ..
DOUBLE PRECISION ZERO
PARAMETER (ZERO=0.0D0)
C .. Local Scalars ..
DOUBLE PRECISION ARGDIF, FDIFF, SCL, TDIFF, TWT, WT, YNRM, YSTGNM
INTEGER L
C .. Intrinsic Functions ..
INTRINSIC ABS, MAX, MIN
C .. Executable Statements ..
C
C Update the weights to account for the current intermediate solution
C approximation YSTG(*). Compute the sizes of Y(*) and YSTG(*) in the
C new norm. The size of the Lipschitz constant is assessed by a difference
C in the arguments Y(*), YSTG(*) and a difference in the function evaluated
C at these arguments.
C
YNRM = ZERO
YSTGNM = ZERO
ARGDIF = ZERO
FDIFF = ZERO
DO 20 L = 1, NEQN
WT = MAX(WEIGHT(L),ABS(YSTG(L)))
WEIGHT(L) = WT
YNRM = MAX(YNRM,ABS(Y(L))/WT)
YSTGNM = MAX(YSTGNM,ABS(YSTG(L))/WT)
ARGDIF = MAX(ARGDIF,ABS(YSTG(L)-Y(L))/WT)
FDIFF = MAX(FDIFF,ABS(YPSTG(L)-YP(L))/WT)
20 CONTINUE
C
C The transformation of the equation to autonomous form is done
C implicitly. The difference of the arguments must take into account
C the difference between the values of the independent variable T and
C TSTG. The difference of the corresponding component of the function
C is zero because of the way the standard transformation is done.
C
TDIFF = TSTG - TNOW
TWT = ABS(TND-TNOW)
YNRM = MAX(YNRM,ABS(TNOW)/TWT)
YSTGNM = MAX(YSTGNM,ABS(TSTG)/TWT)
ARGDIF = MAX(ARGDIF,ABS(TDIFF)/TWT)
C
C The ratio FDIFF/ARGDIF is a lower bound for, and an approximation to, a
C Lipschitz constant L for the differential equation written in autonomous
C form. First we must ask if the difference ARGDIF is significant in the
C precision available. If it appears to be, we insist that abs(HTRY)*L be
C less than an approximate radius, STBRAD, of the stability region of the
C method. This is more stringent than necessary for stability, possibly a
C lot more stringent, but the aim is to get an HTRY small enough that the
C error estimate for the step is credible. The reduction is required to be
C at least as much as the step control parameter RS1. It is necessary to
C limit the reduction of HTRY at any one time because we may be misled in
C the size of the reduction that is appropriate due to nonlinearity of the
C differential equation and to inaccurate weights caused by HTRY much too
C large. The reduction is not permitted to be more than the step control
C parameter RS4.
C
CUTBAK = .FALSE.
IF (ARGDIF.GT.RNDOFF*MAX(YNRM,YSTGNM)) THEN
IF ((ABS(HTRY)*FDIFF).GT.(STBRAD*ARGDIF)) THEN
SCL = (STBRAD*ARGDIF)/(ABS(HTRY)*FDIFF)
SCL = MIN(SCL,RS1)
SCL = MAX(SCL,RS4)
HTRY = SCL*HTRY
CUTBAK = .TRUE.
END IF
END IF
C
RETURN
END SUBROUTINE STEPA
C
SUBROUTINE STEPB(NEQ,Y,YP,H,YNEW,STAGES,THRES,ERR,MAIN,WEIGHT)
C************************************************
C**** NOT A DESIGNATED USER-CALLABLE ROUTINE ****
C************************************************
C
C Purpose: To compute an error estimate for METHD = 2 prior
C to completing the step.
C
C Input: NEQ, Y(*), YP(*), H, STAGES(NEQ,*), THRES(*), MAIN,
C WEIGHT(*)
C Output: ERR
C
C Common: Initializes: none
C Reads: /RKCOM4/ E, PTR
C Alters: none
C
C Comments:
C =========
C If global error assessment is taking place, then MAIN = .FALSE. and
C the weight vector generated by the primary integration is used. The
C error estimate is a linear combination (with coefficients in E(*))
C of the stages stored in STAGES(*,*) (located by PTR(*)).
C
C .. Scalar Arguments ..
DOUBLE PRECISION ERR, H
INTEGER NEQ
LOGICAL MAIN
C .. Array Arguments ..
DOUBLE PRECISION STAGES(NEQ,*), THRES(*), WEIGHT(*), Y(*),
& YNEW(*), YP(*)
C .. Common Block to hold Formula Definitions ..
DOUBLE PRECISION A(13,13), B(13), C(13), BHAT(13), R(11,6),
& E(7)
INTEGER PTR(13), NSTAGE, METHD, MINTP
LOGICAL INTP
COMMON /RKCOM4/ A, B, C, BHAT, R, E, PTR, NSTAGE, METHD,
& MINTP, INTP
SAVE /RKCOM4/
C .. Parameters ..
DOUBLE PRECISION ZERO, HALF
PARAMETER (ZERO=0.0D0,HALF=0.5D0)
C .. Local Scalars ..
DOUBLE PRECISION AVGY, SUM, WT
INTEGER INDEX, L
C .. Intrinsic Functions ..
INTRINSIC ABS, MAX
C .. Executable Statements ..
C
ERR = ZERO
DO 40 L = 1, NEQ
C
C Estimate the local error of component L. The coding makes use of
C E(2) = 0.0D0 and E(7) = 0.0D0.
C
SUM = E(1)*YP(L)
DO 20 INDEX = 3, 6
SUM = SUM + E(INDEX)*STAGES(L,PTR(INDEX))
20 CONTINUE
C
C The local error is H*SUM. A weighted maximum norm of SUM is formed
C and then the factor of H is taken into account.
C
IF (MAIN) THEN
AVGY = HALF*(ABS(Y(L))+ABS(YNEW(L)))
WT = MAX(AVGY,THRES(L))
ELSE
WT = WEIGHT(L)
END IF
C
ERR = MAX(ERR,ABS(SUM/WT))
40 CONTINUE
ERR = ABS(H)*ERR
C
RETURN
END SUBROUTINE STEPB
C
SUBROUTINE STIFF(F,HAVG,JFLSTP,TOOMCH,MAXFCN,WORK,IER,NREC)
C************************************************
C**** NOT A DESIGNATED USER-CALLABLE ROUTINE ****
C************************************************
C
C Purpose: Diagnose stiffness. This depends on two things: whether
C the step size is being restricted on grounds of stability
C and whether the integration to TND can be completed in no
C more than MAXFCN function evaluations.
C
C Input: HAVG, TOOMCH, MAXFCN, WORK(*)
C Input/output: JFLSTP
C Output: IER, NREC
C Workspace: WORK(*)
C External: F
C
C Common: Initializes: /RKCOM9/ REC
C Reads: /RKCOM1/ TND, NEQN
C /RKCOM2/ T, H, NFCN, SVNFCN, OKSTP
C /RKCOM3/ PRY, PRYP, PRTHRS, PRWT, PRSCR,
C PRSTGS, PRYOLD
C /RKCOM5/ COST
C Alters: /RKCOM2/ NFCN
C /RKCOM9/ REC
C
C .. Scalar Arguments ..
DOUBLE PRECISION HAVG
INTEGER IER, JFLSTP, MAXFCN, NREC
LOGICAL TOOMCH
C .. Array Arguments ..
DOUBLE PRECISION WORK(*)
C .. Subroutine Arguments ..
EXTERNAL F
C .. Common Block for Problem Definition ..
DOUBLE PRECISION TSTRT, TND, DIR, HSTRT, TOLR
INTEGER NEQN
COMMON /RKCOM1/ TSTRT, TND, DIR, HSTRT, TOLR, NEQN
SAVE /RKCOM1/
C .. Common Block to hold Problem Status ..
DOUBLE PRECISION T, H, TOLD, HOLD
INTEGER NFCN, SVNFCN, OKSTP, FLSTP
LOGICAL FIRST, LAST
COMMON /RKCOM2/ T, H, TOLD, HOLD, NFCN, SVNFCN, OKSTP, FLSTP,
& FIRST, LAST
SAVE /RKCOM2/
C .. Common Block for General Workspace Pointers ..
INTEGER PRTHRS, PRERST, PRWT, PRYOLD, PRSCR, PRY, PRYP,
& PRSTGS, PRINTP, LNINTP
COMMON /RKCOM3/ PRTHRS, PRERST, PRWT, PRYOLD, PRSCR, PRY, PRYP,
& PRSTGS, PRINTP, LNINTP
SAVE /RKCOM3/
C .. Common Block to hold Formula Characterisitcs ..
DOUBLE PRECISION TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
& RS, RS1, RS2, RS3, RS4
INTEGER ORDER, LSTSTG, MAXTRY, NSEC
LOGICAL FSAL
COMMON /RKCOM5/ TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
& RS, RS1, RS2, RS3, RS4, ORDER, LSTSTG, MAXTRY,
& NSEC, FSAL
SAVE /RKCOM5/
C .. Common Block for Error Message ..
CHARACTER*80 REC(10)
COMMON /RKCOM9/ REC
SAVE /RKCOM9/
C .. Parameters ..
DOUBLE PRECISION HALF
PARAMETER (HALF=0.5D0)
C .. Local Scalars ..
DOUBLE PRECISION AVGY, XTRAWK
INTEGER L
LOGICAL LOTSFL, STIF, UNSURE
C .. External Subroutines ..
EXTERNAL STIFFA
C .. Intrinsic Functions ..
INTRINSIC ABS, DBLE, MAX, MOD
C .. Executable Statements ..
C
IF (MOD(OKSTP-10,40).EQ.0) THEN
LOTSFL = JFLSTP .GE. 10
JFLSTP = 0
ELSE
LOTSFL = .FALSE.
END IF
C
C If either too much work has been done or there are lots of failed steps,
C test for stiffness.
C
IF (TOOMCH .OR. LOTSFL) THEN
C
C Regenerate weight vector
DO 20 L = 1, NEQN
AVGY = HALF*(ABS(WORK(PRY-1+L))+ABS(WORK(PRYOLD-1+L)))
WORK(PRWT-1+L) = MAX(AVGY,WORK(PRTHRS-1+L))
20 CONTINUE
C
C STIFFA determines whether the problem is STIFF. In some circumstances it
C is UNSURE. The decision depends on two things: whether the step size is
C being restricted on grounds of stability and whether the integration to
C TND can be completed in no more than MAXFCN function evaluations. The
C last four arguments of STIFFA are vectors of length NEQN used for working
C storage. Some storage in WORK(*) reserved for the stages (there are a
C minimum of three such vectors reserved for the METHDs implemented) and
C the scratch vector starting at PRSCR are used for this purpose.
C
CALL STIFFA(F,T,WORK(PRY),H,HAVG,TND,MAXFCN,WORK(PRWT),
& WORK(PRYP),WORK(PRERST),UNSURE,STIF,WORK(PRSTGS),
& WORK(PRSTGS+NEQN),WORK(PRSTGS+2*NEQN),WORK(PRSCR))
IF (.NOT.UNSURE) THEN
IF (STIF) THEN
C
C Predict how much eXTRA WorK will be needed to reach TND.
XTRAWK = (COST*ABS((TND-T)/HAVG))/DBLE(SVNFCN+NFCN)
IER = 4
WRITE (REC(NREC+1),'(A)')
&' ** Your problem has been diagnosed as stiff. If the '
WRITE (REC(NREC+2),'(A,D13.5)')
&' ** situation persists, it will cost roughly ', XTRAWK
WRITE (REC(NREC+3),'(A)')
&' ** times as much to reach TEND as it has cost to reach TNOW.'
WRITE (REC(NREC+4),'(A)')
&' ** You should probably change to a code intended for '
WRITE (REC(NREC+5),'(A)')
&' ** stiff problems. '
NREC = NREC + 5
END IF
END IF
END IF
C
RETURN
END SUBROUTINE STIFF
C
SUBROUTINE STIFFA(F,X,Y,HNOW,HAVG,XEND,MAXFCN,WT,FXY,V0,UNSURE,
& STIF,V1,V2,V3,VTEMP)
C************************************************
C**** NOT A DESIGNATED USER-CALLABLE ROUTINE ****
C************************************************
C
C External: F
C Input: X, Y(*), HNOW, HAVG, XEND, MAXFCN, WT(*), FXY(*)
C Input/Output V0(*)
C Output: UNSURE, STIF
C Workspace: V1(*), V2(*), V3(*), VTEMP(*)
C
C Common: Initializes: none
C Reads: /RKCOM1/ TND, NEQN
C /RKCOM5/ COST, STBRAD, TANANG
C /RKCOM7/ SQRRMC, CUBRMC
C Alters: none
C
C STIFFA diagnoses stiffness for an explicit Runge-Kutta code. When it
C is called, either many step failures have been observed, or a lot of
C work has been done.
C
C The NEQ equations of the problem are defined by the subroutine F(X,Y,YP).
C When STIFFA is called, the integration has reached X where the approximate
C solution is Y(*). The vector FXY(*) is defined by a call of F(X,Y,FXY).
C It is an input argument because it is usually available from the integrator.
C
C The last successful step was of size HNOW, and an average step size is
C HAVG. A weighted norm is used to measure the local error with the error
C in solution component L divided by the positive weight WT(L) provided in
C the vector WT(*).
C
C Explicit Runge - Kutta codes estimate the local error of Y(*) by
C forming the difference of two approximate solutions. This difference
C must be provided in the vector V0(*). When this difference is too
C small to be significant, STIFFA will replace it with a "random" vector.
C
C STIF is set .TRUE. when the average step size appears to be restricted
C on grounds of stability. In certain cases the variable UNSURE is set
C .TRUE.; the value of STIF is then not defined.
C
C The stability region of the explicit Runge-Kutta formula is described
C by quantities TANANG and STBRAD that are communicated by the setup routine
C via COMMON. Stability regions often change sharply near the imaginary
C axis so that it is difficult to classify the stiffness of a problem with
C eigenvalues of a local Jacobian that are "near" the imaginary axis. For
C this reason, we consider only points Z in the upper left half complex
C plane for which TAN( IMAG(Z)/( - RE(Z))) <= TANANG. Eigenvalues outside
C this region are one reason for the code being UNSURE. The stability
C region is approximated by the intersection of a disk with this sector.
C The radius of this disk is called STBRAD.
C
C Working storage must be provided via the four vectors V1(*),V2(*),
C V3(*),VTEMP(*). These vectors must be of length at least NEQ.
C
C .. Scalar Arguments ..
DOUBLE PRECISION HAVG, HNOW, X, XEND
INTEGER MAXFCN
LOGICAL STIF, UNSURE
C .. Array Arguments ..
DOUBLE PRECISION FXY(*), V0(*), V1(*), V2(*), V3(*), VTEMP(*),
& WT(*), Y(*)
C .. Subroutine Arguments ..
EXTERNAL F
C .. Common Block for Problem Definition ..
DOUBLE PRECISION TSTRT, TND, DIR, HSTRT, TOLR
INTEGER NEQN
COMMON /RKCOM1/ TSTRT, TND, DIR, HSTRT, TOLR, NEQN
SAVE /RKCOM1/
C .. Common Block to hold Formula Characterisitcs ..
DOUBLE PRECISION TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
& RS, RS1, RS2, RS3, RS4
INTEGER ORDER, LSTSTG, MAXTRY, NSEC
LOGICAL FSAL
COMMON /RKCOM5/ TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG,
& RS, RS1, RS2, RS3, RS4, ORDER, LSTSTG, MAXTRY,
& NSEC, FSAL
SAVE /RKCOM5/
C .. Common Block for Environment Parameters ..
DOUBLE PRECISION MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY
INTEGER OUTCH
COMMON /RKCOM7/ MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY,
& OUTCH
SAVE /RKCOM7/
C .. Parameters ..
DOUBLE PRECISION LARGE
PARAMETER (LARGE=1.0D+10)
DOUBLE PRECISION ZERO, P001, P9, ONE, TWO, FIVE, FIFTH
PARAMETER (ZERO=0.0D+0,P001=0.001D+0,P9=0.9D+0,ONE=1.0D+0,
& TWO=2.0D+0,FIVE=5.0D+0,FIFTH=0.2D+0)
C .. Local Scalars ..
DOUBLE PRECISION ALPHA1, ALPHA2, BETA1, BETA2, D1, D2, DET1,
& DET2, DIST, RES2, RHO, RHO2, ROLD, SCALE, V0NRM,
& V0V0, V0V1, V0V2, V1V1, V1V2, V1V3, V2V2, V2V3,
& V3NRM, V3V3, XTRFCN, YNRM
INTEGER L, NTRY
LOGICAL ROOTRE
C .. Local Arrays ..
DOUBLE PRECISION R1(2), R2(2), ROOT1(2), ROOT2(2)
C .. External Functions ..
DOUBLE PRECISION DOTPRD
EXTERNAL DOTPRD
C .. External Subroutines ..
EXTERNAL STIFFB, STIFFC, STIFFD
C .. Intrinsic Functions ..
INTRINSIC ABS, MIN, SQRT
C .. Executable Statements ..
C
C If the current step size differs substantially from the average,
C the problem is not stiff.
C
IF (ABS(HNOW/HAVG).GT.FIVE .OR. ABS(HNOW/HAVG).LT.FIFTH) THEN
STIF = .FALSE.
UNSURE = .FALSE.
RETURN
ELSE
UNSURE = .TRUE.
END IF
C
C The average step size is used to predict the cost in function evaluations
C of finishing the integration to XEND. If this cost is no more than MAXFCN,
C the problem is declared not stiff: If the step size is being restricted on
C grounds of stability, it will stay close to HAVG. The prediction will
C then be good, but the cost is too low to consider the problem stiff. If
C the step size is not close to HAVG, the problem is not stiff. Either way
C there is no point to testing for a step size restriction due to stability.
C
XTRFCN = COST*ABS((XEND-X)/HAVG)
IF (XTRFCN.LE.MAXFCN) THEN
STIF = .FALSE.
UNSURE = .FALSE.
RETURN
ELSE
UNSURE = .TRUE.
END IF
C
C There have been many step failures or a lot of work has been done. Now
C we must determine if this is due to the stability characteristics of the
C formula. This is done by calculating the dominant eigenvalues of the
C local Jacobian and then testing whether HAVG corresponds to being on the
C boundary of the stability region.
C
C The size of Y(*) provides scale information needed to approximate
C the Jacobian by differences.
C
YNRM = SQRT(DOTPRD(Y,Y,WT,NEQN))
SCALE = YNRM*SQRRMC
IF (SCALE.EQ.ZERO) THEN
C
C Degenerate case. Y(*) is (almost) the zero vector so the scale is not
C defined. The input vector V0(*) is the difference between Y(*) and a
C lower order approximation to the solution that is within the error
C tolerance. When Y(*) vanishes, V0(*) is itself an acceptable approximate
C solution, so we take SCALE from it, if this is possible.
C
YNRM = SQRT(DOTPRD(V0,V0,WT,NEQN))
SCALE = YNRM*SQRRMC
IF (SCALE.EQ.ZERO) THEN
UNSURE = .TRUE.
RETURN
END IF
END IF
C
V0V0 = DOTPRD(V0,V0,WT,NEQN)
IF (V0V0.EQ.ZERO) THEN
C
C Degenerate case. V0(*) is (almost) the zero vector so cannot
C be used to define a direction for an increment to Y(*). Try a
C "random" direction.
C
DO 20 L = 1, NEQN
V0(L) = ONE
20 CONTINUE
V0V0 = DOTPRD(V0,V0,WT,NEQN)
END IF
V0NRM = SQRT(V0V0)
DO 40 L = 1, NEQN
V0(L) = V0(L)/V0NRM
40 CONTINUE
V0V0 = ONE
C
C Use a nonlinear power method to estimate the two dominant eigenvalues.
C V0(*) is often very rich in the two associated eigenvectors. For this
C reason the computation is organized with the expectation that a minimal
C number of iterations will suffice. Indeed, it is necessary to recognize
C a kind of degeneracy when there is a dominant real eigenvalue. The
C subroutine STIFFB does this. In the first try, NTRY = 1, a Rayleigh
C quotient for such an eigenvalue is initialized as ROLD. After each
C iteration, REROOT computes a new Rayleigh quotient and tests whether the
C two approximations agree to one tenth of one per cent and the eigenvalue,
C eigenvector pair satisfy a stringent test on the residual. ROOTRE = .TRUE.
C signals that a single dominant real root has been found.
C
NTRY = 1
60 CONTINUE
C
CALL STIFFD(V0,HAVG,X,Y,F,FXY,WT,SCALE,V0V0,V1,V1V1,VTEMP)
C
C The quantity SQRT(V1V1/V0V0) is a lower bound for the product of HAVG
C and a Lipschitz constant. If it should be LARGE, stiffness is not
C restricting the step size to the stability region. The principle is
C clear enough, but the real reason for this test is to recognize an
C extremely inaccurate computation of V1V1 due to finite precision
C arithmetic in certain degenerate circumstances.
C
IF (SQRT(V1V1).GT.LARGE*SQRT(V0V0)) THEN
UNSURE = .TRUE.
RETURN
END IF
C
V0V1 = DOTPRD(V0,V1,WT,NEQN)
IF (NTRY.EQ.1) THEN
ROLD = V0V1/V0V0
C
C This is the first Rayleigh quotient approximating the product of HAVG
C and a dominant real eigenvalue. If it should be very small, the
C problem is not stiff. It is important to test for this possibility so
C as to prevent underflow and degeneracies in the subsequent iteration.
C
IF (ABS(ROLD).LT.CUBRMC) THEN
UNSURE = .FALSE.
STIF = .FALSE.
RETURN
END IF
ELSE
CALL STIFFB(V1V1,V0V1,V0V0,ROLD,RHO,ROOT1,ROOT2,ROOTRE)
IF (ROOTRE) GO TO 100
END IF
CALL STIFFD(V1,HAVG,X,Y,F,FXY,WT,SCALE,V1V1,V2,V2V2,VTEMP)
V0V2 = DOTPRD(V0,V2,WT,NEQN)
V1V2 = DOTPRD(V1,V2,WT,NEQN)
CALL STIFFB(V2V2,V1V2,V1V1,ROLD,RHO,ROOT1,ROOT2,ROOTRE)
IF (ROOTRE) GO TO 100
C
C Fit a quadratic in the eigenvalue to the three successive iterates
C V0(*),V1(*),V2(*) of the power method to get a first approximation to
C a pair of eigenvalues. A test made earlier in STIFFB implies that
C the quantity DET1 here will not be too small.
C
DET1 = V0V0*V1V1 - V0V1**2
ALPHA1 = (-V0V0*V1V2+V0V1*V0V2)/DET1
BETA1 = (V0V1*V1V2-V1V1*V0V2)/DET1
C
C Iterate again to get V3, test again for degeneracy, and then fit a
C quadratic to V1(*),V2(*),V3(*) to get a second approximation to a pair
C of eigenvalues.
C
CALL STIFFD(V2,HAVG,X,Y,F,FXY,WT,SCALE,V2V2,V3,V3V3,VTEMP)
V1V3 = DOTPRD(V1,V3,WT,NEQN)
V2V3 = DOTPRD(V2,V3,WT,NEQN)
CALL STIFFB(V3V3,V2V3,V2V2,ROLD,RHO,ROOT1,ROOT2,ROOTRE)
IF (ROOTRE) GO TO 100
DET2 = V1V1*V2V2 - V1V2**2
ALPHA2 = (-V1V1*V2V3+V1V2*V1V3)/DET2
BETA2 = (V1V2*V2V3-V2V2*V1V3)/DET2
C
C First test the residual of the quadratic fit to see if we might
C have determined a pair of eigenvalues.
C
RES2 = ABS(V3V3+V2V2*ALPHA2**2+V1V1*BETA2**2+TWO*V2V3*ALPHA2+
& TWO*V1V3*BETA2+TWO*V1V2*ALPHA2*BETA2)
IF (RES2.LE.V3V3*P001**2) THEN
C
C Calculate the two approximate pairs of eigenvalues.
C
CALL STIFFC(ALPHA1,BETA1,R1,R2)
CALL STIFFC(ALPHA2,BETA2,ROOT1,ROOT2)
C
C The test for convergence is done on the larger root of the second
C approximation. It is complicated by the fact that one pair of roots
C might be real and the other complex. First calculate the spectral
C radius RHO of HAVG*J as the magnitude of ROOT1. Then see if one of
C the roots R1,R2 is within one per cent of ROOT1. A subdominant root
C may be very poorly approximated if its magnitude is much smaller than
C RHO -- this does not matter in our use of these eigenvalues.
C
RHO = SQRT(ROOT1(1)**2+ROOT1(2)**2)
D1 = (ROOT1(1)-R1(1))**2 + (ROOT1(2)-R1(2))**2
D2 = (ROOT1(1)-R2(1))**2 + (ROOT1(2)-R2(2))**2
DIST = SQRT(MIN(D1,D2))
IF (DIST.LE.P001*RHO) GO TO 100
END IF
C
C Do not have convergence yet. Because the iterations are cheap, and
C because the convergence criterion is stringent, we are willing to try
C a few iterations.
C
IF (NTRY.LT.MAXTRY) THEN
NTRY = NTRY + 1
V3NRM = SQRT(V3V3)
DO 80 L = 1, NEQN
V0(L) = V3(L)/V3NRM
80 CONTINUE
V0V0 = ONE
GO TO 60
ELSE
UNSURE = .TRUE.
RETURN
END IF
C
C **************
C
C We now have the dominant eigenvalues. Decide if the average step
C size is being restricted on grounds of stability. Check the real
C parts of the eigenvalues. First see if the dominant eigenvalue is
C in the left half plane -- there won't be a stability restriction
C unless it is. If there is another eigenvalue of comparable magnitude
C with a positive real part, the problem is not stiff. If the dominant
C eigenvalue is too close to the imaginary axis, we cannot diagnose
C stiffness.
C
100 CONTINUE
IF (ROOT1(1).GT.ZERO) THEN
STIF = .FALSE.
UNSURE = .FALSE.
RETURN
END IF
RHO2 = SQRT(ROOT2(1)**2+ROOT2(2)**2)
IF (RHO2.GE.P9*RHO .AND. ROOT2(1).GT.ZERO) THEN
STIF = .FALSE.
UNSURE = .FALSE.
RETURN
END IF
IF (ABS(ROOT1(2)).GT.ABS(ROOT1(1))*TANANG) THEN
UNSURE = .TRUE.
RETURN
END IF
C
C If the average step size corresponds to being well within the
C stability region, the step size is not being restricted because
C of stability.
C
STIF = RHO .GE. P9*STBRAD
UNSURE = .FALSE.
RETURN
END SUBROUTINE STIFFA
C
SUBROUTINE STIFFB(V1V1,V0V1,V0V0,ROLD,RHO,ROOT1,ROOT2,ROOTRE)
C************************************************
C**** NOT A DESIGNATED USER-CALLABLE ROUTINE ****
C************************************************
C
C Input: V1V1, V0V1, V0V0
C Input/output: ROLD
C Output: RHO, ROOT1(*),ROOT2(*),ROOTRE
C
C Decide if the iteration has degenerated because of a strongly
C dominant real eigenvalue. Have just computed the latest iterate.
C V1V1 is its dot product with itself, V0V1 is the dot product
C of the previous iterate with the current one, and V0V0 is the
C dot product of the previous iterate with itself. ROLD is a
C previous Rayleigh quotient approximating a dominant real
C eigenvalue. It must be computed directly the first time the
C subroutine is called. It is updated each call to STIFFB, hence
C is available for subsequent calls.
C
C If there is a strongly dominant real eigenvalue, ROOTRE is set
C .TRUE., ROOT1(*) returns the eigenvalue, RHO returns the magnitude
C of the eigenvalue, and ROOT2(*) is set to zero.
C
C .. Scalar Arguments ..
DOUBLE PRECISION RHO, ROLD, V0V0, V0V1, V1V1
LOGICAL ROOTRE
C .. Array Arguments ..
DOUBLE PRECISION ROOT1(2), ROOT2(2)
C .. Parameters ..
DOUBLE PRECISION ZERO, P001
PARAMETER (ZERO=0.0D+0,P001=0.001D+0)
C .. Local Scalars ..
DOUBLE PRECISION DET, R, RES
C .. Intrinsic Functions ..
INTRINSIC ABS
C .. Executable Statements ..
C
R = V0V1/V0V0
RHO = ABS(R)
DET = V0V0*V1V1 - V0V1**2
RES = ABS(DET/V0V0)
ROOTRE = DET .EQ. ZERO .OR. (RES.LE.V1V1*P001**2 .AND.
& ABS(R-ROLD).LE.P001*RHO)
IF (ROOTRE) THEN
ROOT1(1) = R
ROOT1(2) = ZERO
ROOT2(1) = ZERO
ROOT2(2) = ZERO
END IF
ROLD = R
C
RETURN
END SUBROUTINE STIFFB
C
SUBROUTINE STIFFC(ALPHA,BETA,R1,R2)
C************************************************
C**** NOT A DESIGNATED USER-CALLABLE ROUTINE ****
C************************************************
C
C Input: ALPHA, BETA
C Output: R1(*), R2(*)
C
C This subroutine computes the two complex roots R1 and R2 of
C the quadratic equation X**2 + ALPHA*X + BETA = 0. The magnitude
C of R1 is greater than or equal to the magnitude of R2. R1 and R2 are
C returned as vectors of two components with the first being the real
C part of the complex number and the second being the imaginary part.
C
C .. Scalar Arguments ..
DOUBLE PRECISION ALPHA, BETA
C .. Array Arguments ..
DOUBLE PRECISION R1(2), R2(2)
C .. Parameters ..
DOUBLE PRECISION ZERO, TWO
PARAMETER (ZERO=0.0D+0,TWO=2.0D+0)
C .. Local Scalars ..
DOUBLE PRECISION DISC, SQDISC, TEMP
C .. Intrinsic Functions ..
INTRINSIC ABS, SQRT
C .. Executable Statements ..
TEMP = ALPHA/TWO
DISC = TEMP**2 - BETA
IF (DISC.EQ.ZERO) THEN
C
C Double root.
C
R1(1) = -TEMP
R1(2) = ZERO
R2(1) = R1(1)
R2(2) = R1(2)
RETURN
END IF
C
SQDISC = SQRT(ABS(DISC))
IF (DISC.LT.ZERO) THEN
C
C Complex conjugate roots.
C
R1(1) = -TEMP
R1(2) = SQDISC
R2(1) = R1(1)
R2(2) = -R1(2)
ELSE
C
C Real pair of roots. Calculate the bigger one in R1(1).
C
IF (TEMP.GT.ZERO) THEN
R1(1) = -TEMP - SQDISC
ELSE
R1(1) = -TEMP + SQDISC
END IF
R1(2) = ZERO
R2(1) = BETA/R1(1)
R2(2) = ZERO
END IF
C
RETURN
END SUBROUTINE STIFFC
C
SUBROUTINE STIFFD(V,HAVG,X,Y,F,FXY,WT,SCALE,VDOTV,Z,ZDOTZ,VTEMP)
C************************************************
C**** NOT A DESIGNATED USER-CALLABLE ROUTINE ****
C************************************************
C
C External: F
C Input: V(*), HAVG, X, Y(*), FXY(*), WT(*), SCALE, VDOTV,
C Output: Z(*), ZDOTZ
C Workspace: VTEMP(*)
C
C For an input vector V(*) of length NEQ, this subroutine computes a vector
C Z(*) that approximates the product HAVG*J*V where HAVG is an input scalar
C and J is the Jacobian matrix of a function F evaluated at the input
C arguments (X,Y(*)). This function is defined by a subroutine of the form
C F(T,U,F) that when given T and U(*), returns the value of the function in
C F(*). The input vector FXY(*) is defined by F(X,Y,FXY). Scaling is a
C delicate matter. A weighted Euclidean norm is used with the (positive)
C weights provided in WT(*). The input scalar SCALE is the square root of
C the unit roundoff times the norm of Y(*). The square of the norm of the
C input vector V(*) is input as VDOTV. The routine outputs the square of
C the norm of the output vector Z(*) as ZDOTZ. The subroutine calls the
C DOUBLE PRECISION FUNCTION DOTPRD(U,V,WT,NEQ) to compute the dot (inner)
C product. The vector VTEMP(*) is used for working storage.
C
C .. Scalar Arguments ..
DOUBLE PRECISION HAVG, SCALE, VDOTV, X, ZDOTZ
C .. Array Arguments ..
DOUBLE PRECISION FXY(*), V(*), VTEMP(*), WT(*), Y(*), Z(*)
C .. Subroutine Arguments ..
EXTERNAL F
C .. Common Block for Problem Definition ..
DOUBLE PRECISION TSTRT, TND, DIR, HSTRT, TOLR
INTEGER NEQN
COMMON /RKCOM1/ TSTRT, TND, DIR, HSTRT, TOLR, NEQN
SAVE /RKCOM1/
C .. Common Block to hold Problem Status ..
DOUBLE PRECISION T, H, TOLD, HOLD
INTEGER NFCN, SVNFCN, OKSTP, FLSTP
LOGICAL FIRST, LAST
COMMON /RKCOM2/ T, H, TOLD, HOLD, NFCN, SVNFCN, OKSTP, FLSTP,
& FIRST, LAST
SAVE /RKCOM2/
C .. Local Scalars ..
DOUBLE PRECISION TEMP1, TEMP2
INTEGER L
C .. External Functions ..
DOUBLE PRECISION DOTPRD
EXTERNAL DOTPRD
C .. Intrinsic Functions ..
INTRINSIC SQRT
C .. Executable Statements ..
C
C Scale V(*) so that it can be used as an increment to Y(*)
C for an accurate difference approximation to the Jacobian.
C
TEMP1 = SCALE/SQRT(VDOTV)
DO 20 L = 1, NEQN
VTEMP(L) = Y(L) + TEMP1*V(L)
20 CONTINUE
C
CALL F(X,VTEMP,Z)
NFCN = NFCN + 1
C
C Form the difference approximation. At the same time undo
C the scaling of V(*) and introduce the factor of HAVG.
C
TEMP2 = HAVG/TEMP1
DO 40 L = 1, NEQN
Z(L) = TEMP2*(Z(L)-FXY(L))
40 CONTINUE
C
ZDOTZ = DOTPRD(Z,Z,WT,NEQN)
C
RETURN
END SUBROUTINE STIFFD
C
FUNCTION DOTPRD(U,V,WT,NEQ)
C************************************************
C**** NOT A DESIGNATED USER-CALLABLE ROUTINE ****
C************************************************
C
C Purpose: To compute a weighted Euclidean dot (inner) product of
C two vectors.
C
C Input: U(*), V(*), WT(*), NEQ
C Output: the result DOTPRD is returned via the subprogram name
C
C Comments:
C =========
C The vectors U(*), V(*), and WT(*) are of length NEQ. The components
C of WT(*) are weights that must be non-zero.
C
C .. Scalar Arguments ..
DOUBLE PRECISION DOTPRD
INTEGER NEQ
C .. Array Arguments ..
DOUBLE PRECISION U(*), V(*), WT(*)
C .. Parameters ..
DOUBLE PRECISION ZERO
PARAMETER (ZERO=0.0D0)
C .. Local Scalars ..
DOUBLE PRECISION SUM
INTEGER L
C .. Executable Statements ..
C
SUM = ZERO
DO 20 L = 1, NEQ
SUM = SUM + (U(L)/WT(L))*(V(L)/WT(L))
20 CONTINUE
C
DOTPRD = SUM
C
RETURN
END FUNCTION DOTPRD
C
SUBROUTINE SOFTFL(ASK,ON)
C
C Purpose: To prevent a program STOP after a "catastrophic"
C failure when using a routine from RKSUITE.
C
C Input: ASK
C Input/output: ON
C
C Comments:
C =========
C When a "catastrophic" failure is detected, the default action of
C RKSUITE is to write an explanation to the standard output channel,
C OUTCH, and STOP. This subroutine can be used to prevent the STOP and
C so allow the main program to continue. To do this, you call SOFTFL with
C ASK = .FALSE. and ON = .TRUE. You must then call the subroutine CHKFL
C after every call to a user-callable routine in RKSUITE to check whether
C a catastrophic error occurred and take appropriate action if it did. Of
C course, you may call SETUP at any time to start a new problem, but calling
C any other user-callable routine in RKSUITE after a catastrophic error will
C lead to a STOP (even when "soft failure" has been set "on").
C
C When ON is set by a call to SOFTFL with ASK = .FALSE., the value of ON
C is SAVEd. The subroutine RKMSG in RKSUITE calls SOFTFL with ASK = .TRUE.
C to find out the SAVEd value of ON.
C
C .. Scalar Arguments ..
LOGICAL ASK, ON
C .. Local Scalars ..
LOGICAL SOFT
C .. Save statement ..
SAVE SOFT
C .. Data statements ..
DATA SOFT/.FALSE./
C .. Executable Statements ..
C
IF (ASK) THEN
ON = SOFT
ELSE
SOFT = ON
END IF
C
RETURN
END SUBROUTINE SOFTFL
C
SUBROUTINE CHKFL(ASK,ERROR)
C
C Purpose: Enquiry routine used in conjunction with SOFTFL.
C Reports whether a "catastrophic" error was detected.
C
C Input: ASK
C Input/output: ERROR
C
C Comments:
C =========
C When a "catastrophic" failure is detected, the default action of
C RKSUITE is to write an explanation to the standard output channel,
C OUTCH, and STOP. SOFTFL can be used to prevent the STOP and so
C allow the main program to continue. It is then necessary to call
C CHKFL with ASK = .TRUE. after every call to a user-callable routine
C in RKSUITE to check whether a catastrophic error occurred and take
C appropriate action if it did. If there was a catastrophic error,
C ERROR is returned .TRUE. Of course, you may call SETUP at any time
C to start a new problem, but calling any other user-callable routine
C in RKSUITE after a catastrophic error will lead to a STOP (even when
C "soft failure" has been set "on").
C
C When a catastrophic failure (IER = 911) is detected in one of
C the routines in RKSUITE, it calls CHKFL with ASK = .FALSE. and
C ERROR = .TRUE. This value of ERROR is SAVEd.
C
C .. Scalar Arguments ..
LOGICAL ASK, ERROR
C .. Local Scalars ..
LOGICAL SAVERR
C .. Save statement ..
SAVE SAVERR
C .. Data statements ..
DATA SAVERR/.FALSE./
C .. Executable Statements ..
C
IF (ASK) THEN
ERROR = SAVERR
ELSE
SAVERR = ERROR
END IF
C
RETURN
END SUBROUTINE CHKFL
C
! END MODULE RKSUITE