C Copyright C.A.H. Paul 1992 C Research funded by the Science & Engineering Research Council C and the Numerical Algorithms Group (Oxford) Ltd C Subsequent funding by SERC grant number GR/H59237 C C C ******************************************************************+++++++ C C SUBROUTINE ABCNRM(INTGRS,REALS,ERRORI) C DOUBLE PRECISION PLYVAL C INTEGER INTGRS(*),LOOP DOUBLE PRECISION ERRORI(*),MAXMUM,REALS(*),VALUE C EXTERNAL PLYVAL C C Routine to calculate the maximum error in the continuous C solution at the abscissae of the Runge-Kutta method C DO 1 LOOP = 1,INTGRS(13) C C Evaluate the error polynomial at each of the abscissae C MAXMUM = PLYVAL(INTGRS,INTGRS(13),LOOP, + REALS(INTGRS(47)+INTGRS(13)),1D0/5D0) VALUE = PLYVAL(INTGRS,INTGRS(13),LOOP, + REALS(INTGRS(47)+INTGRS(13)),3D0/10D0) IF (VALUE.GT.MAXMUM) MAXMUM = VALUE VALUE = PLYVAL(INTGRS,INTGRS(13),LOOP, + REALS(INTGRS(47)+INTGRS(13)),1D0/2D0) IF (VALUE.GT.MAXMUM) MAXMUM = VALUE VALUE = PLYVAL(INTGRS,INTGRS(13),LOOP, + REALS(INTGRS(47)+INTGRS(13)),4D0/5D0) IF (VALUE.GT.MAXMUM) MAXMUM = VALUE VALUE = PLYVAL(INTGRS,INTGRS(13),LOOP, + REALS(INTGRS(47)+INTGRS(13)),8D0/9D0) IF (VALUE.GT.MAXMUM) MAXMUM = VALUE VALUE = PLYVAL(INTGRS,INTGRS(13),LOOP, + REALS(INTGRS(47)+INTGRS(13)),1D0) IF (VALUE.GT.MAXMUM) MAXMUM = VALUE C C Store the maximum value C ERRORI(LOOP) = MAXMUM*REALS(11) C 1 CONTINUE C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE ASYMP(INTGRS,REALS,DIMEN,KI,ERRORI) C INTEGER DIMEN,I,INTGRS(*) DOUBLE PRECISION REALS(*),KI(DIMEN,9),ERRORI(*) C C Routine to calculate the asymptotic embedded error estimate C of Dormand and Prince method C DO 1 I = 1,DIMEN ERRORI(I) = ((71D0/57600D0)*KI(I,1)-(71D0/16695D0)*KI(I,3)+ + (71D0/1920D0)*KI(I,4)-(17253D0/339200D0)*KI(I,5)+ + (22D0/525D0)*KI(I,6)-(1D0/40D0)*KI(I,7))*REALS(11) 1 CONTINUE C RETURN END C C C ******************************************************************+++++++ C C LOGICAL FUNCTION CHECK(INTGRS,REALS,IN,ORDER,AT) C DOUBLE PRECISION DABS,DMAX1 C INTEGER IN,INTGRS(*),ORDER,POINT,SEARCH,TEMP DOUBLE PRECISION AT,DIFF,DIST,REALS(*) C INTRINSIC DABS,DMAX1 C C Function to determine whether a lower order discontinuity C already exists in the solution component IN at the point AT C C Calculate a minimum distance between tracked discontinuities C before they are considered as being distinct, taking into C account the user's choice `REALS(16)', the position of the C discontinuity and the unit round off C DIST = DMAX1(REALS(16),3D0*DMAX1(DABS(AT),1D0)*REALS(1)) C C Get the first continuity record after the current point C SEARCH = INTGRS(61) POINT = INTGRS(39)+SEARCH*2-1 C C Search upto the last continuity record stored C 1 IF (SEARCH.LE.INTGRS(10)) THEN IF (INTGRS(POINT-1).EQ.IN) THEN C C If the solution components match, check the positions C of the discontinuities allowing for rounding error) C DIFF = DABS(REALS(INTGRS(55)+SEARCH-1)-AT) C C If the positions `match', set the order of the discontinuity C to the lowest of ORDER and the one already stored in INTGRS C and return to the calling routine C IF (DIFF.LE.DIST) THEN CHECK = .TRUE. C C Set the order of the discontinuity to the lowest C IF (ORDER.LT.INTGRS(POINT)) THEN TEMP = ORDER ORDER = INTGRS(POINT) INTGRS(POINT) = TEMP ENDIF C RETURN ENDIF ENDIF C C Increment the search variable and repeat C SEARCH = SEARCH+1 POINT = POINT+2 GOTO 1 ENDIF C C Searched all the records, so signal that that is does not exist C CHECK = .FALSE. C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE CHENRM(INTGRS,REALS,ERRORI) C DOUBLE PRECISION DCOS,PLYVAL C INTEGER INTGRS(*),LOOP DOUBLE PRECISION ERRORI(*),MAXMUM,PI,REALS(*),VALUE PARAMETER (PI = 3.14159265358979323D0) C EXTERNAL PLYVAL INTRINSIC DCOS C C Routine to calculate the maximum error in the continuous C solution at the (non-trivial) extremal values of the fifth C degree Chebyshev polynomial C DO 1 LOOP = 1,INTGRS(13) C C Evaluate the error polynomial at each of the abscissae C MAXMUM = PLYVAL(INTGRS,INTGRS(13),LOOP, + REALS(INTGRS(47)+INTGRS(13)), + (1D0/2D0)*(1D0+DCOS(PI/5D0))) VALUE = PLYVAL(INTGRS,INTGRS(13),LOOP, + REALS(INTGRS(47)+INTGRS(13)), + (1D0/2D0)*(1D0+DCOS(2D0*PI/5D0))) IF (VALUE.GT.MAXMUM) MAXMUM = VALUE VALUE = PLYVAL(INTGRS,INTGRS(13),LOOP, + REALS(INTGRS(47)+INTGRS(13)), + (1D0/2D0)*(1D0+DCOS(3D0*PI/5D0))) IF (VALUE.GT.MAXMUM) MAXMUM = VALUE VALUE = PLYVAL(INTGRS,INTGRS(13),LOOP, + REALS(INTGRS(47)+INTGRS(13)), + (1D0/2D0)*(1D0+DCOS(4D0*PI/5D0))) IF (VALUE.GT.MAXMUM) MAXMUM = VALUE VALUE = PLYVAL(INTGRS,INTGRS(13),LOOP, + REALS(INTGRS(47)+INTGRS(13)),1D0) IF (VALUE.GT.MAXMUM) MAXMUM = VALUE C C Store the maximum value C ERRORI(LOOP) = MAXMUM*REALS(11) C 1 CONTINUE C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE CONTY(INTGRS,REALS,IN,ORDER,AT) C INTEGER IN,INTGRS(*),ORDER,POS DOUBLE PRECISION AT,REALS(*) C EXTERNAL ERRHAN,MAPSET C C Routine to specify the continuity of components of the initial C function and initial point C C IN - Component of solution in which discontinuity occurs C ORDER - The order of continuity of the discontinuity C AT - The position at which the discontinuity occurs C C Only check the user-supplied continuity records (determined by C the number of derivative evaluations) C IF (INTGRS(25).EQ.0) THEN C C Ensure that IN does not exceed the dimension of the problem C IF (IN.LT.1 .OR. IN.GT.INTGRS(13)) + CALL ERRHAN(INTGRS, + 'Discontinuity occurs in a solution that cannot exist.', + 'CONTY.',-1) C C Ensure that the order of the discontinuity is non-negative C IF (ORDER.LT.0) + CALL ERRHAN(INTGRS, + 'Orders of discontinuities must be non-negative.', + 'CONTY.',-1) C C If this is the first call to either CONTY, (N)HARD or (N)SOFT set C up the pointers to start of arrays in INTGRS. After which the number C of Volterra terms and storage allocated to the continuity records C and the HARD and SOFT discontinuity networks cannot be altered. C C INTGRS(8)=0 & INTGRS(9)=0 & INTGRS(10)=0 C IF (INTGRS(8)+INTGRS(9)+INTGRS(10).EQ.0) + CALL MAPSET(INTGRS,'CONTY.') ENDIF C C Check that there is still space left for mapping C IF (INTGRS(10).GE.INTGRS(21)) THEN CALL ERRHAN(INTGRS, + 'Insufficient space for continuity record.', + 'CONTY.',2) IF (INTGRS(15).GE.2) + WRITE (*,*) 'Continuity record ignored' RETURN ENDIF C C Save the continuity record: IN, ORDER C POS = INTGRS(39)+2*INTGRS(10) INTGRS(POS) = IN INTGRS(POS+1) = ORDER POS = INTGRS(55)+INTGRS(10) REALS(POS) = AT C C Update the number of continuity records C INTGRS(10) = INTGRS(10)+1 C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE CUTOFF(INTGRS,REALS,STOPAT) C INTEGER INTGRS(*),STOPAT DOUBLE PRECISION REALS(*) C EXTERNAL ERRHAN C C Optional routine to specify after what order of continuity C discontinuities are no longer tracked C C Ensure that the order is not less than 0 (default = 5) C IF (STOPAT.LT.0) + CALL ERRHAN(INTGRS, + 'Order of discontinuity tracking too low.','CUTOFF.',-1) C C If STOPAT = 0, report only neutral discontinuity tracking is done C IF (STOPAT.EQ.0 .AND. INTGRS(15).GE.2) + CALL ERRHAN(INTGRS, + 'Tracking of order 0 neutral discontinuities only.', + 'CUTOFF.',2) C C Check if iterative refinement selected, if it is we should be C tracking discontinuities upto order 5 C IF (INTGRS(7).EQ.0 .AND. STOPAT.LT.5) CALL ERRHAN(INTGRS, + 'Refinement needs discontinuities upto order 5 to be tracked.', + 'CUTOFF.',2) C INTGRS(17) = STOPAT C RETURN END C C C ******************************************************************+++++++ C C DOUBLE PRECISION FUNCTION DELAY(INTGRS,REALS,LAGFN,SOLN) C DOUBLE PRECISION DABS,RECALL,YINIT C INTEGER INTGRS(*),LAGFN,SOLN DOUBLE PRECISION LAGVAL,REALS(*) C EXTERNAL ERRHAN,RECALL,YINIT INTRINSIC DABS C C Routine to retrieve the SOLN'th component of the solution C at the point indicated by the LAG'th lag value. The initial C solution must be given by the function YINIT C C Ensure that the given lag function and solution exist. C Use kernel number (INTGRS(58)) to determine what type of C lag function we are evaluating. Only check information ONCE C IF (INTGRS(25).EQ.0) THEN C IF (LAGFN.LT.0 .OR. (INTGRS(58).EQ.0 .AND. LAGFN.GT.INTGRS(14)) + .OR. (INTGRS(58).NE.0 .AND. LAGFN.GT.INTGRS(14)+INTGRS(44))) + CALL ERRHAN(INTGRS,'Specified lag function cannot exist.', + 'DELAY.',-1) C IF (SOLN.LT.1 .OR. SOLN.GT.INTGRS(13)) + CALL ERRHAN(INTGRS, + 'Specified solution component cannot exist.', + 'DELAY.',-1) ENDIF C C Increment the `number of back evaluations' diagnostic C INTGRS(28) = INTGRS(28)+1 C C If an ODE term is requested (the zeroth lag function), and C we are not in iterative refinement mode or tracking state C dependent discontinuities, use the appropriate component of C Y1 otherwise use the iterated continuous extension C C LAGFN=0 & INTGRS(31)=0 & INTGRS(57)=0 C IF (LAGFN+INTGRS(31)+INTGRS(57).EQ.0) THEN DELAY = REALS(INTGRS(35)+SOLN-1) RETURN ENDIF C C Calculate the lag-point C LAGVAL = REALS(INTGRS(45)+LAGFN-1) C C If we evaluate a back-solution `just before' the initial point, C set the `re-evaluate first stage' flag to ON. This is only C (really) necessary when discontinuities are NOT tracked AND C there is a jump in the solution at the initial point. C IF (INTGRS(60).NE.1 .AND. DABS(REALS(7)-LAGVAL).LT.3D0*REALS(1)) + INTGRS(59) = 1 C C Check if the initial function is to be called: C If we are evaluating the derivative at the start of the C interval (INTGRS(60)=1) or not evaluating a Runge-Kutta stage, C (INTGRS(60)=0) use a .LT. comparison C IF (LAGVAL.LE.REALS(7)) THEN DELAY = YINIT(SOLN,LAGVAL) IF (LAGVAL.LT.REALS(7) .OR. INTGRS(60).GT.1) RETURN ENDIF C C Ensure that if iterative refinement is signalled, we C at least return a valid value for the solution (the value C at the start of the current interval) C IF (INTGRS(31).EQ.1) THEN DELAY = REALS(INTGRS(34)+SOLN-1) RETURN ENDIF C C Call the actual function which calculates the solution value C passing TYPE as .TRUE. to indicate that we want to evaluate C a delay value and not a neutral value, and the position of C the storage arrays in REALS. C DELAY = RECALL(INTGRS,REALS,.TRUE.,LAGFN,LAGVAL,SOLN, + REALS(INTGRS(34)),INTGRS(13), + REALS(INTGRS(48)),REALS(INTGRS(49)), + REALS(INTGRS(50)),REALS(INTGRS(51)), + REALS(INTGRS(52)),REALS(INTGRS(53)), + REALS(INTGRS(54))) C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE DELAYS(INTGRS,REALS,NUMBER) C INTEGER INTGRS(*),NUMBER DOUBLE PRECISION REALS(*) C EXTERNAL ERRHAN C C Routine for specifying the number of delays in the equation C IF (NUMBER.LE.0) + CALL ERRHAN(INTGRS, + 'Invalid number of delays specified.','DELAYS.',-1) C C Check that the number of delay functions has not already been set C IF (INTGRS(14).GT.1) CALL ERRHAN(INTGRS, + 'The number of delay functions has already been set.', + 'DELAYS.',-1) C INTGRS(14) = NUMBER C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE DETERM(INTGRS,REALS) C DOUBLE PRECISION LAG C INTEGER COUNT,INTGRS(*),LOOP DOUBLE PRECISION REALS(*) C EXTERNAL LAG C C Routine for determining whether lag functions are state dependent C or independent by evaluating each lag in turn. To avoid requiring C solution values from the current interval, we evaluate the lag at C the initial point. A state dependent delay is detected by checking C if the number of approximant evaluations has increased C C First set the ODE lag function to the initial point C REALS(INTGRS(45)-1) = REALS(7) C C Do for each lag function C DO 1 LOOP = 1,INTGRS(14) C C Save the current number of approximant evaluations C COUNT = INTGRS(28) C C Evaluate the lag function: LAG NUMBER, AT C REALS(INTGRS(45)+LOOP-1) = LAG(INTGRS,REALS,LOOP,REALS(7)) C C Check if the number of approximant evaluations has increased C IF (INTGRS(28).EQ.COUNT) THEN C C Mark this delay as state independent C INTGRS(INTGRS(43)+LOOP-1) = 0 ELSE C C Mark this delay as state dependent C INTGRS(INTGRS(43)+LOOP-1) = 1 ENDIF C 1 CONTINUE C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE DIMEN(INTGRS,REALS,SIZES) C INTEGER INTGRS(*),SIZES DOUBLE PRECISION REALS(*) C EXTERNAL ERRHAN C C Set the dimension of the problem (optional: default is 1) C IF (SIZES.LT.1) + CALL ERRHAN(INTGRS,'Invalid dimension of problem.','DIMEN.',-1) C C Check that the dimension has not already been set C IF (INTGRS(13).GT.1) CALL ERRHAN(INTGRS, + 'The dimension of the problem has already been set.', + 'DIMEN.',-1) C INTGRS(13) = SIZES C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE DOHARD(INTGRS,REALS,MAPPNG,CHANGE) C LOGICAL CHECK C INTEGER FROM,INTAT,INTGRS(*),INTMAP,INTO INTEGER MAPPNG,ORDER,SEARCH,SMOOTH,TEMP DOUBLE PRECISION DISCTO,REALS(*) LOGICAL CHANGE,EXIST,NOFULL C EXTERNAL CHECK,CONTY C C Subroutine to update the continuity records based on the C HARD discontinuity networks mappings. MAPPNG is the number C of the HARD discontinuity mapping. SEARCH is set to the C first continuity record added with this call of NXTDSC C C Convert the HARD discontinuity mapping number into the position C in INTGRS for the mapping C INTMAP = INTGRS(40)+MAPPNG*2-2 C C Get the solution components that the discontinuities are propagated C from and into C FROM = INTGRS(INTMAP) INTO = INTGRS(INTMAP+1) C C Calculate the degree of smoothing C IF (INTO.LT.0) THEN C C For neutral terms assume no smoothing C INTO = -INTO SMOOTH = 0 ELSE C C For ODE, delay and Volterra terms assume one degree of smoothing C SMOOTH = 1 ENDIF C C Search all the continuity records added with this call of NXTDSC, C including those added by previous calls to this routine C C Get the number of the first continuity record added in this call C SEARCH = INTGRS(56) C C Calculate the position in INTGRS of this continuity record C INTAT = INTGRS(39)+SEARCH*2-2 C C Repeat for all continuity records after the first, so long as the C continuity record is not full C NOFULL = (INTGRS(10).LT.INTGRS(21)) 1 IF (NOFULL .AND. SEARCH.LE.INTGRS(10)) THEN C C If the solution components agree, get the point of the discontinuity C and its order C IF (FROM.EQ.INTGRS(INTAT)) THEN DISCTO = REALS(SEARCH+INTGRS(55)-1) ORDER = INTGRS(INTAT+1)+SMOOTH C C If the order of the propagated discontinuity is beyond C CUTOFF, skip to next continuity record C IF (ORDER.LE.INTGRS(17)) THEN C TEMP = ORDER C C Check to see if it already exists C EXIST = CHECK(INTGRS,REALS,INTO,TEMP,DISCTO) C C If TEMP has changed then set CHANGE to .TRUE. C CHANGE = (TEMP.NE.ORDER) C C If EXIST is .FALSE. then add the propagated discontinuity C to the continuity records and set CHANGE to .TRUE. C IF (.NOT.EXIST) THEN CALL CONTY(INTGRS,REALS,INTO,ORDER,DISCTO) CHANGE = .TRUE. NOFULL = (INTGRS(10).LT.INTGRS(21)) ENDIF ENDIF ENDIF C C Increment the search variable and repeat C SEARCH = SEARCH+1 INTAT = INTAT+2 GOTO 1 ENDIF C C End of propagating for this HARD discontinuity mapping C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE DOTRCK(INTGRS,REALS,POS,LIMIT) C DOUBLE PRECISION DMIN1 LOGICAL CHECK C INTEGER CODE,FROM,IN,INTAT,INTGRS(*) INTEGER LINK,NEWORD,POS,SEARCH,SMOOTH DOUBLE PRECISION DISCAT,LIMIT,REALS(*),ROOT LOGICAL EXIST,NOFULL,UNDONE C EXTERNAL CHECK,CONTY,SECANT INTRINSIC DMIN1 C C Routine to track discontinuities propagated by the SOFT C discontinuity mapping pointed to by POS. The root finder C used is the Secant method. C C Get the SOFT discontinuity mapping C FROM = INTGRS(POS) IN = INTGRS(POS+1) LINK = INTGRS(POS+2) C C Calculate the degree of smoothing C IF (IN.LT.0) THEN C C Degree of smoothing for neutral equations C IN = -IN SMOOTH = 0 ELSE C C Degree of smoothing for ODE, delay and Volterra terms C SMOOTH = 1 ENDIF C C Search for continuity records for discontinuities upto the C current point where the discontinuity occurs in the same C solution component as specified in FROM C C Initialize the pointer to the first continuity record that can C be propagated (effected by size of LBOUND) C SEARCH = INTGRS(62) C C Convert the continuity record number into a position in INTGRS C where the solution component and order of the discontinuity are C INTAT = INTGRS(39)+SEARCH*2-2 UNDONE = .TRUE. C C Repeat for all valid continuity records, so long as the continuity C record is not full C NOFULL = (INTGRS(10).LT.INTGRS(21)) 1 IF (NOFULL .AND. SEARCH.LT.INTGRS(56)) THEN C C Convert the continuity record number into the position of the C discontinuity (stored in REALS) C DISCAT = REALS(SEARCH+INTGRS(55)-1) C C Only track the discontinuity if it is within the upper bound C for the size of the lag(s) C IF (DISCAT+REALS(18).LE.REALS(9)) THEN IF (UNDONE) INTGRS(62) = SEARCH+1 GOTO 2 ENDIF UNDONE = .FALSE. C C Ensure that the solution components match C IF (INTGRS(INTAT).NE.FROM) GOTO 2 C C Calculate the new order of the propagated discontinuity C NEWORD = INTGRS(INTAT+1)+SMOOTH C C Ensure that the order of the propagated discontinuity is not C higher than CUTOFF C IF (NEWORD.GT.INTGRS(17)) GOTO 2 C C All the conditions are met, so do a Secant search to see if C the discontinuity is propagated into the interval [Tn,Tn+2h] C CALL SECANT(INTGRS,REALS,LINK,DISCAT,ROOT,CODE) C C If return CODE = -1 failure, CODE = 1 or 2 success C C Save the continuity record if successful without requiring C extrapolation for state dependent lag functions and a lower C order discontinuity does not already exist in this solution C at this point C IF (CODE.EQ.1) THEN EXIST = CHECK(INTGRS,REALS,IN,NEWORD,ROOT) IF (.NOT.EXIST) THEN CALL CONTY(INTGRS,REALS,IN,NEWORD,ROOT) NOFULL = (INTGRS(10).LT.INTGRS(21)) ENDIF ELSE C C If we were successful, but required extrapolation to `find' C the discontinuity - do not add the point to the continuity C records, but reduce LIMIT so that the next integration step C does not cross this approximate discontinuity C IF (CODE.EQ.2) LIMIT = DMIN1(LIMIT,ROOT) ENDIF C C Continuation of search loop C 2 SEARCH = SEARCH+1 INTAT = INTAT+2 GOTO 1 ENDIF C C End of search loop C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE ERRCAL(INTGRS,REALS,EST) C DOUBLE PRECISION DABS,DBLE,DMAX1,DSQRT C INTEGER INTGRS(*),LOOP,TYPE,YPOS,Y1POS DOUBLE PRECISION DENOM,EST,REALS(*) C EXTERNAL ABCNRM,ASYMP,CHENRM,ERRPLY,INFNRM,ONENRM INTRINSIC DABS,DBLE,DMAX1,DSQRT C C Selection routine for the five different error estimators C C INTGRS(6) = 1 - Asymptotic (default) C = 2 - One norm of continuous error estimate C = 3 - Infinity norm of cts error estimate C = 4 - Infinity norm of cts error estimate at abscissae points C = 5 - Infinity norm of cts error estimate at Chebyshev points C TYPE = INTGRS(6) C C Calculate the position in REALS of Y and Y1 C YPOS = INTGRS(34) Y1POS = INTGRS(35) C C Deal with the asymptotic error estimate first. The DIMEN C individual error estimates are stored in the first DIMEN elements C of the storage allocated for the calculation of the continuous C error polynomial in REALS, which start at INTGRS(46) and is of C length DIMEN for the asymptotic error estimate and DIMEN*6 for C the continuous error estimates C IF (TYPE.EQ.1) THEN CALL ASYMP(INTGRS,REALS,INTGRS(13),REALS(INTGRS(36)), + REALS(INTGRS(47))) C C Branch to where the 2-norm of all the individual relative error C estimates is calculated, to give one number representing the C overall error C GOTO 1 ENDIF C C Calculate the continuous error polynomial for all other cases. C CALL ERRPLY(INTGRS,REALS,INTGRS(13),REALS(INTGRS(36)), + REALS(INTGRS(47)+INTGRS(13))) C C Call the appropriate error estimator C IF (TYPE.EQ.2) THEN CALL ONENRM(INTGRS,REALS,INTGRS(13),REALS(INTGRS(47)), + REALS(INTGRS(47)+INTGRS(13))) ELSE IF (TYPE.EQ.3) THEN CALL INFNRM(INTGRS,REALS,INTGRS(13),REALS(INTGRS(47)), + REALS(INTGRS(47)+INTGRS(13))) ELSE IF (TYPE.EQ.4) THEN CALL ABCNRM(INTGRS,REALS,REALS(INTGRS(47))) ELSE CALL CHENRM(INTGRS,REALS,REALS(INTGRS(47))) ENDIF ENDIF ENDIF C C Calculate the 2-norm of the individual relative error estimates C 1 EST = 0D0 DO 2 LOOP = 0,INTGRS(13)-1 DENOM = DMAX1(40D0*REALS(1),DABS(REALS(YPOS+LOOP)), + DABS(REALS(Y1POS+LOOP)),2D0*REALS(1)/REALS(2)) EST = EST+REALS(INTGRS(47)+LOOP)*REALS(INTGRS(47)+LOOP)/ + (DENOM*DENOM) 2 CONTINUE C EST = DSQRT(EST/DBLE(INTGRS(13))) C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE ERRHAN(INTGRS,ERROR,WHERE,TYPE) C INTEGER INDEX C INTEGER INTGRS(*),TYPE,UPTO1,UPTO2 CHARACTER ERROR*(*),WHERE*(*) C INTRINSIC INDEX C C Routine for reporting errors and where they occur C C INTGRS(15) - level of error reporting C = 0 - no reporting and automatic error correction C = 1 - reporting of extrapolation and iterative refinement C = 2 - reporting of all errors but automatic correction C = 3 - report error and stop C TYPE states whether the error is recoverable or not C = -1 - fatal error C = 1 - extrapolation or iterative refinement used C = 2 - recoverable error C C If no errors to be printed and the error is recoverable, return C IF (INTGRS(15).EQ.0 .AND. TYPE.GT.0) RETURN IF (INTGRS(15).EQ.1 .AND. TYPE.GT.1) RETURN C C Strip out trailing spaces in error message and origin of error C UPTO1 = INDEX(ERROR,'.') UPTO2 = INDEX(WHERE,'.') C IF (UPTO1.EQ.0 .OR. UPTO2.EQ.0) THEN C C Invalid error message was passed - which means faulty code C WRITE (*,*) 'Invalid error message passed to ERRHAN.' STOP ENDIF C C Report error and its origin C WRITE (*,*) ERROR(1:UPTO1-1),' in ',WHERE(1:UPTO2) C C Handle fatal errors and level 3 reporting C IF (INTGRS(15).EQ.3 .OR. TYPE.EQ.-1) THEN IF (TYPE.EQ.-1) WRITE (*,*) 'Fatal error - unable to continue' STOP ENDIF C C The error is recoverable and automatic correction selected C so continue C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE ERROR(INTGRS,REALS,TYPE,EPS) C INTEGER INTGRS(*),TYPE DOUBLE PRECISION EPS,REALS(*) C EXTERNAL ERRHAN C C Routine for specifying type of error control C C TYPE = 0 - No error control (only selectable by fixed stepsize option) C = 1 - Asymptotic (default) C = 2 - One norm of continuous error estimate C = 3 - Infinity norm of continuous error estimate C = 4 - Infinity norm of continuous error estimate at abscissae points C = 5 - Infinity norm of continuous error estimate at Chebyshev points C C Ignore this call if fixed stepsize strategy already selected C IF (INTGRS(4).EQ.1) THEN CALL ERRHAN(INTGRS, + 'Fixed stepsize strategy already selected.','ERROR.',2) IF (INTGRS(15).GE.2) WRITE (*,*) + 'Request to use error-per-step control ignored' RETURN ENDIF C C Check that type of error control exists C IF (TYPE.LT.1 .OR. TYPE.GT.5) THEN CALL ERRHAN(INTGRS, + 'Error control specified does not exist.','ERROR.',2) IF (INTGRS(15).GE.2) + WRITE (*,*) 'Current error control retained' RETURN ENDIF C C Ensure that EPS is a valid size C IF (EPS.LE.REALS(1)) THEN CALL ERRHAN(INTGRS, + 'Error-per-step is less than unit roundoff.','ERROR.',2) IF (INTGRS(15).GE.2) WRITE (*,*) + 'Current error-per-step control ',REALS(2),' retained' ELSE REALS(2) = EPS ENDIF C INTGRS(6) = TYPE C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE ERRPLY(INTGRS,REALS,DIMEN,KI,CTSPLY) C INTEGER DIMEN,I,INTGRS(*) DOUBLE PRECISION CTSPLY(DIMEN,5),KI(DIMEN,9),REALS(*) C C Routine to calculate a continuous fourth order C error polynomial - the difference between the fifth C order Hermite interpolant and the fourth order C continuous extension (full-step fourth order as well) C DO 1 I = 1,DIMEN CTSPLY(I,1)=0D0 CTSPLY(I,2)=(-6307D0/5920D0)*KI(I,1)+ + (12248D0/1961D0)*KI(I,3)+(4599D0/592D0)*KI(I,4)- + (903231D0/313760D0)*KI(I,5)+(4422D0/13135D0)*KI(I,6)- + (3483D0/2627D0)*KI(I,7)-(40D0/37D0)*KI(I,8)- + 8D0*KI(I,9) CTSPLY(I,3)=(212723D0/53280D0)*KI(I,1)- + (413752D0/17649D0)*KI(I,3)-(52717D0/1776D0)*KI(I,4)+ + (3595671D0/313760D0)*KI(I,5)-(84326D0/39405D0)*KI(I,6)+ + (14847D0/2627D0)*KI(I,7)+(80D0/37D0)*KI(I,8)+ + 32D0*KI(I,9) CTSPLY(I,4)=(-3366191D0/710400D0)*KI(I,1)+ + (5760191D0/205905D0)*KI(I,3)+(889809D0/23680D0)*KI(I,4)- + (205438761D0/12550400D0)*KI(I,5)+ + (2293148D0/459725D0)*KI(I,6)-(877507D0/105080D0)*KI(I,7)- + (40D0/37D0)*KI(I,8)-40D0*KI(I,9) CTSPLY(I,5)=(29D0/16D0)*KI(I,1)-(4000D0/371D0)*KI(I,3)- + (125D0/8D0)*KI(I,4)+(6561D0/848D0)*KI(I,5)- + (22D0/7D0)*KI(I,6)+4D0*KI(I,7)+16D0*KI(I,9) 1 CONTINUE C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE EXTRAP(INTGRS,REALS,MODE) C INTEGER INTGRS(*),MODE DOUBLE PRECISION REALS(*) C EXTERNAL ERRHAN C C Routine for specifying that extrapolation instead of C iterative refinement is to be carried out where possible C MODE = 1 for no extrapolation (and no iterative refinement) C MODE = 2 for extrapolation in the current interval only. C MODE = 3 for extrapolation beyond the current interval. C C Check that a valid extrapolation mode has been requested. C IF (MODE.LT.1 .OR. MODE.GT.3) THEN CALL ERRHAN(INTGRS, + 'Invalid extrapolation mode requested','EXTRAP.',2) IF (INTGRS(15).GE.2) WRITE (*,*) + 'Current extrapolation mode/iterative refinement retained.' RETURN ENDIF C C Report that extrapolation is to be carried out C IF (INTGRS(15).GE.2 .AND. MODE.GT.1) + WRITE (*,*) 'Extrapolation selected' C INTGRS(7) = MODE C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE FINDRT(REALS,POLY,ROOTS) C DOUBLE PRECISION DABS,DMAX1 C INTEGER COUNT,LOOP,ROOT DOUBLE PRECISION F,G,POLY(0:3),REALS(*),ROOTS(6),SWAP,T LOGICAL FOUND C INTRINSIC DABS,DMAX1 C C Routine to find the roots of POLY in the interval [0,1] C and return them in ROOTS in ascending order C C Initialize the ROOTS array C ROOT = 2 ROOTS(1) = 0D0 DO 1 LOOP = 2,6 ROOTS(LOOP) = 1D0 1 CONTINUE C C Main loop for finding roots of POLY C C Initialize counter for the number of Newton iterations C and the starting point for the iteration C 2 COUNT = 0 T = 0D0 C C Newton iteration loop C C Ensure roots are in the interval [0,1] C 3 IF (T.GT.1D0) THEN T = 0.5D0 ELSE IF (T.LT.0D0) T = 1D0 ENDIF C C Calculate the polynomial value at T C F = POLY(0)+T*(POLY(1)+T*(POLY(2)+T*POLY(3))) C C If the current point is not a root do a Newton step C FOUND = (DABS(F).GT.5D0*DMAX1(DABS(F),1D0)*REALS(1)) IF (FOUND) THEN C C Increment the number of iterations C COUNT = COUNT+1 C C Calculate the derivative C G = POLY(1)+T*(2D0*POLY(2)+T*3D0*POLY(3)) C C Handle a zero derivative C IF (DABS(G).LT.REALS(1)) THEN C C If zero derivative, perturb the current point C T = T+1000D0*REALS(1) ELSE T = T-F/G ENDIF ENDIF C C Repeat the Newton iteration if no root found and C we have not done too many iterations C IF (FOUND .AND. COUNT.LE.10) GOTO 3 C C If a root was found, store it in ROOTS C IF (.NOT.FOUND) THEN ROOTS(ROOT) = T ROOT = ROOT+1 C C Deflate the polynomial C POLY(0) = POLY(1)+T*(POLY(2)+T*POLY(3)) POLY(1) = POLY(2)+T*POLY(3) POLY(2) = POLY(3) POLY(3) = 0D0 C C If a root was found, try another Newton iteration C IF (ROOT.LT.6) GOTO 2 ENDIF C C Found all possible roots, now order them using bubble-sort C DO 5 COUNT = 1,5 DO 4 LOOP = COUNT+1,6 IF (ROOTS(COUNT).GT.ROOTS(LOOP)) THEN SWAP = ROOTS(COUNT) ROOTS(COUNT) = ROOTS(LOOP) ROOTS(LOOP) = SWAP ENDIF 4 CONTINUE 5 CONTINUE C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE FIXSTP(INTGRS,REALS,STEP) C DOUBLE PRECISION DABS,DMAX1 C INTEGER INTGRS(*) DOUBLE PRECISION REALS(*),STEP C EXTERNAL ERRHAN INTRINSIC DABS,DMAX1 C C Routine for imposing a constant stepsize strategy C C This ignores tracking of discontinuities, error-per-step control, C the initial stepsize and the maximum allowed stepsize C C Ensure that FIXSTP is only called after the start and endpoints C have been specified C IF (INTGRS(1).EQ.0) + CALL ERRHAN(INTGRS, + 'Range of integration must be specified first.', + 'FIXSTP.',-1) C C Check for valid fixed stepsize C IF (STEP.LT.0D0) + CALL ERRHAN(INTGRS,'Fixed stepsize is negative.','FIXSTP.',-1) IF (STEP.LE.10D0*DMAX1(1D0,DABS(REALS(7)))*REALS(1)) + CALL ERRHAN(INTGRS,'Fixed stepsize is too small.','FIXSTP.',-1) IF (REALS(7)+STEP.GT.REALS(8)) + CALL ERRHAN(INTGRS,'Fixed stepsize is too large.','FIXSTP.',-1) C C Check if tracking discontinuities has already been selected C IF (INTGRS(5).GT.0) THEN CALL ERRHAN(INTGRS, + 'Fixed stepsize conflicts with tracking discontinuities.', + 'FIXSTP.',2) IF (INTGRS(15).GE.2) WRITE (*,*) + 'Request to track discontinuities ignored' INTGRS(5) = 0 ENDIF C C Check if error-per-step control has already been selected C IF (INTGRS(6).GT.0) THEN CALL ERRHAN(INTGRS, + 'Fixed stepsize conflicts with error-per-step control.', + 'FIXSTP.',2) IF (INTGRS(15).GE.2) WRITE (*,*) + 'Request to use error-per-step control ignored' INTGRS(6) = 0 ENDIF C C Ignore the initial stepsize if set C IF (INTGRS(2).EQ.1) THEN CALL ERRHAN(INTGRS, + 'Fixed stepsize conflicts with initial stepsize.', + 'FIXSTP.',2) IF (INTGRS(15).GE.2) WRITE (*,*) + 'Request to set initial stepsize ignored' INTGRS(2) = 0 ENDIF C C Ignore the maximum stepsize if set C IF (INTGRS(3).EQ.1) THEN CALL ERRHAN(INTGRS, + 'Fixed stepsize conflicts with maximum stepsize.', + 'FIXSTP.',2) IF (INTGRS(15).GE.2) WRITE (*,*) + 'Request to set maximum permitted stepsize ignored' INTGRS(3) = 0 ENDIF C C Flag that fixed stepsize strategy has been selected C INTGRS(4) = 1 C REALS(11) = STEP REALS(12) = STEP*2D0 REALS(13) = STEP C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE HARD(INTGRS,REALS,FROM,IN) C INTEGER FROM,IN,INTGRS(*),POS DOUBLE PRECISION REALS(*) C EXTERNAL ERRHAN,MAPSET C C Routine to specify a HARD delay discontinuity network link C C FROM - Component of solution in which discontinuity already exists C IN - Component of solution in which discontinuity is created C C Ensure that FROM does not exceed the dimension of the problem C IF (FROM.LT.1 .OR. FROM.GT.INTGRS(13)) + CALL ERRHAN(INTGRS, + 'Discontinuity occurs in a solution that cannot exist.', + 'HARD.',-1) C C Ensure that IN does not exceed the dimension of the problem C IF (IN.LT.1 .OR. IN.GT.INTGRS(13)) + CALL ERRHAN(INTGRS, + 'Discontinuity mapped to a solution that cannot exist.', + 'HARD.',-1) C IF (FROM.EQ.IN) THEN CALL ERRHAN(INTGRS,'Solution components must be different.', + 'HARD.',2) IF (INTGRS(15).GE.2) WRITE (*,*) + 'HARD discontinuity network mapping ignored' RETURN ENDIF C C If this is the first call to either CONTY, (N)HARD or (N)SOFT set C up the pointers to start of arrays in INTGRS. After which the number C of Volterra terms and storage allocated to the continuity records C and the HARD and SOFT discontinuity networks cannot be altered. C C INTGRS(8)=0 & INTGRS(9)=0 & INTGRS(10)=0 C IF (INTGRS(8)+INTGRS(9)+INTGRS(10).EQ.0) + CALL MAPSET(INTGRS,'HARD.') C C Check that there is still space left for mapping C IF (INTGRS(8).EQ.INTGRS(19)) THEN CALL ERRHAN(INTGRS, + 'Insufficient space for HARD discontinuity network.', + 'HARD.',2) IF (INTGRS(15).GE.2) + WRITE (*,*) 'HARD discontinuity network mapping ignored' RETURN ENDIF C C Save the HARD discontinuity mapping: IN, FROM C POS = INTGRS(40)+2*INTGRS(8) INTGRS(POS) = FROM INTGRS(POS+1) = IN C C Update the number of HARD discontinuity mappings set C INTGRS(8) = INTGRS(8)+1 C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE HSTART(INTGRS,REALS,STEP) C DOUBLE PRECISION DABS,DMAX1 C INTEGER INTGRS(*) DOUBLE PRECISION REALS(*),STEP C EXTERNAL ERRHAN INTRINSIC DABS,DMAX1 C C Routine for specifying the initial stepsize C C Ensure that HSTART is only called after the start and endpoints C have been specified C IF (INTGRS(1).EQ.0) + CALL ERRHAN(INTGRS, + 'Range of integration must be specified first.', + 'HSTART.',-1) C C Ignore this call if the fixed stepsize strategy has been selected C IF (INTGRS(4).EQ.1) THEN CALL ERRHAN(INTGRS, + 'Fixed stepsize strategy already selected.','HSTART.',2) IF (INTGRS(15).GE.2) WRITE (*,*) + 'Request to set initial stepsize ignored' RETURN ENDIF C C Ensure that the initial stepsize if positive C IF (STEP.LE.0D0) + CALL ERRHAN(INTGRS, + 'Initial stepsize must be positive.','HSTART.',-1) C C Ensure that the initial stepsize is large enough C IF (STEP.LE.10D0*DMAX1(1D0,DABS(REALS(7)))*REALS(1)) THEN CALL ERRHAN(INTGRS, + 'Initial stepsize is too small.','HSTART.',2) C C Set the initial stepsize to a `reasonable' size C STEP = 10000D0*DMAX1(1D0,DABS(REALS(7)))*REALS(1) IF (INTGRS(15).GE.2) WRITE (*,*) + 'Initial stepsize has been set to ',STEP ENDIF C C Ensure that the initial stepsize is small enough C IF (REALS(7)+STEP.GT.REALS(8)) THEN CALL ERRHAN(INTGRS, + 'Initial stepsize is too large.','HSTART.',2) C C Set the initial stepsize to the largest possible value C STEP = REALS(8)-REALS(7) C C Do not violate the maximum stepsize if set C IF (INTGRS(3).EQ.1 .AND. STEP.GT.REALS(12)) STEP = REALS(12) C IF (INTGRS(15).GE.2) WRITE (*,*) + 'The initial stepsize has been reduced to ',STEP ENDIF C C Do not violate the maximum stepsize for all initial stepsizes C IF (INTGRS(3).EQ.1 .AND. STEP.GT.REALS(12)) THEN CALL ERRHAN(INTGRS, + 'Initial stepsize exceeds maximum permitted stepsize.', + 'HSTART.',2) STEP = REALS(12) IF (INTGRS(15).GE.2) WRITE (*,*) + 'The initial stepsize has been reduced to ',STEP ENDIF C C Flag that HSTART has been specified C INTGRS(2) = 1 C REALS(11) = STEP C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE INFNRM(INTGRS,REALS,DIMEN,ERRORI,CTSPLY) C DOUBLE PRECISION DABS,DBLE C INTEGER DIMEN,INTGRS(*),LOOP,SOLN DOUBLE PRECISION MAXMUM,POLY(0:3),ROOTS(6),T,VALUE DOUBLE PRECISION CTSPLY(DIMEN,5),ERRORI(*),REALS(*) C EXTERNAL FINDRT INTRINSIC DABS,DBLE C C Routine to calculate the infinity-norm of the error polynomial C C Repeat for all components of the solution C DO 4 SOLN = 1,DIMEN C C Since the maxima and minima of the error polynomial occur C either at the start or end of the interval or at one of the C zeros of the derivative of the error polynomial, find C these zeros C C Transfer the correct component of the error polynomial C to POLY to allow for deflation when searching for roots C DO 1 LOOP = 0,3 POLY(LOOP) = CTSPLY(SOLN,LOOP+2)*DBLE(LOOP+2) 1 CONTINUE C C Find the roots in the interval C CALL FINDRT(REALS,POLY,ROOTS) C C Find the zero root in the ordered sequence C LOOP = 0 2 LOOP = LOOP+1 IF (DABS(ROOTS(LOOP)).GT.REALS(1)) GOTO 2 C C Calculate the maximum value of the error polynomial in C the interval [0,1] C MAXMUM = 0D0 C C The search loop C 3 LOOP = LOOP+1 T = ROOTS(LOOP) VALUE = DABS(T*T*(CTSPLY(SOLN,2)+T*(CTSPLY(SOLN,3)+ + T*(CTSPLY(SOLN,4)+T*CTSPLY(SOLN,5))))) IF (VALUE.GT.MAXMUM) MAXMUM = VALUE C C Repeat until the end of the continuous interval T = 1 C IF (DABS(ROOTS(LOOP)-1D0).GT.REALS(1)) GOTO 3 C C Store this component of the error estimate and repeat C for all components of the problem C ERRORI(SOLN) = MAXMUM*REALS(11) 4 CONTINUE C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE INIT(INTGRS,REALS,DIMINT,DIMREL) C INTEGER ISTART,RSTART PARAMETER (ISTART=65,RSTART=20) INTEGER DIMINT,DIMREL,INTGRS(*),LOOP DOUBLE PRECISION REALS(*),TEMP C EXTERNAL ERRHAN C C Routine to initialize default values and check storage allocation C Initialize INTGRS and REALS where necessary C C Check that the size of INTGRS is not too small to hold defaults C IF (DIMINT.LT.ISTART-1) + CALL ERRHAN(INTGRS, + 'Size of INTGRS is much too small.','INIT.',-1) C C Check that the size of REALS is not too small to hold defaults C IF (DIMREL.LT.RSTART-1) + CALL ERRHAN(INTGRS, + 'Size of REALS is much too small.','INIT.',-1) C C Comment out later initializations to zero C DO 1 LOOP = 1,ISTART-1 INTGRS(LOOP) = 0 1 CONTINUE DO 2 LOOP = 1,RSTART-1 REALS(LOOP) = 0D0 2 CONTINUE C C Set the start of the continuity records in INTGRS and REALS C INTGRS(39) = ISTART INTGRS(55) = RSTART C C Flag that the interval of integration has not been set C C INTGRS(1) = 0 C C Flag that the initial stepsize has not been specified C C INTGRS(2) = 0 C C Flag that the maximum permitted stepsize has not been specified C C INTGRS(3) = 0 C C Flag that the fixed stepsize strategy is not being used (default) C C INTGRS(4) = 0 C C Flag that tracking discontinuities is not selected (default) C C INTGRS(5) = 0 C C Flag that no error control has been selected (default) C C INTGRS(6) = 0 C C Flag that extrapolation for small delays is used (default) C INTGRS(7) = 2 C C Initialize the number of HARD discontinuity networks C C INTGRS(8) = 0 C C Initialize the number of SOFT discontinuity networks C C INTGRS(9) = 0 C C Initialize the number of discontinuity records C C INTGRS(10) = 0 C C Save the size of REALS allocated C INTGRS(11) = DIMREL C C Save the size of INTGRS allocated C INTGRS(12) = DIMINT C C Set the dimension of the delay problem (default) C INTGRS(13) = 1 C C Set the number of lag terms (default) C INTGRS(14) = 1 C C Set the level of error reporting (default) C INTGRS(15) = 3 C C Set the level of code diagnostics generated (default) C INTGRS(16) = 1 C C Set upto what order discontinuities are tracked (default) C INTGRS(17) = 5 C C Set the number of Volterra terms in functional equation (default) C INTGRS(18) = 1 C C Set the maximum number of HARD discontinuity networks (default) C INTGRS(19) = 25 C C Set the maximum number of SOFT discontinuity networks (default) C INTGRS(20) = 25 C C Set the maximum number of continuity records (default) C INTGRS(21) = 200 C C Initialize diagnostics: C number of steps C number of accepted steps C number of rejected steps C number of derivative function evaluations C number of kernel evaluations C number of extrapolations C number of approximant evaluations C C DO 3 LOOP = 22,28 C INTGRS(LOOP) = 0 C 3 CONTINUE C C Initialize the position in the cyclic storage of the first element C INTGRS(29) = 1 C C Initialize the position in the cyclic storage of the last element C C INTGRS(30) = 0 C C Set the iterative refinement flag to `off' C C INTGRS(31) = 0 C C Set the number of kernel lag functions C INTGRS(44) = 1 C C Initialize the `next discontinuity record before NXTDSC called' C pointer - also used to determine whether calls to CONTY come C from the driver program or from tracking discontinuities C C INTGRS(56) = 0 C C Set the `extrapolate during tracking discontinuities' flag to `off' C C INTGRS(57) = 0 C C Set the kernel number to 0 (see VOLTRA). This is used to determine C what type of lag function we are evaluating C C INTGRS(58) = 0 C C Set kernel smoothness flag to `smooth' C C INTGRS(64) = 0 C C Further initialization of INTGRS is postponed until any of the C default settings have been changed C C Initialize REALS, but first calculate unit roundoff using a C SCRATCH file to avoid problems with optimizing compilers C OPEN (UNIT=44,STATUS='SCRATCH') REALS(1) = 1D0 4 REALS(1) = REALS(1)/2D0 TEMP = 1D0+REALS(1)/2D0 WRITE (44,*) TEMP IF (TEMP.GT.1D0) GOTO 4 CLOSE (UNIT=44) C C Set the default error-per-step tolerance at 1000*unit roundoff C REALS(2) = 1000D0*REALS(1) C C Set the minimum distance between tracked discontinuities before C they are treated as distinct discontinuities C REALS(16) = 0.001D0 C C Set an upper bound for the size of the lag(s) (only used when tracking) C REALS(18) = 1D30 C C Further initialization of REALS is postponed until the exact C problem is specified C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE INTEG(INTGRS,REALS,DIMEN,Y,KI) C DOUBLE PRECISION DABS,DMAX1 C INTEGER DIMEN,I,INTGRS(*),POS DOUBLE PRECISION H,KI(DIMEN,9),REALS(*),T,Y(*) C EXTERNAL DERIV,LAGVAL INTRINSIC DABS,DMAX1 C C Routine for solving the Runge-Kutta equations. If iterative C refinement is being carried out, the current estimate for y C is given by the iterated approximant and not by Y1, but we C still need estimates for Y(Tn+H/2) and Y(Tn+H) to allow C construction of the iterated continuous extension C C Ensure that iterative refinement is not flagged as being C necessary, if it is exit C IF (INTGRS(31).EQ.1) RETURN C C As a shortcut put the current stepsize in H and the position C of Y1 in REALS in POS C H = REALS(11) POS = INTGRS(35)-1 C C Calculate Y(Tn+H/5) for all dimensions of problem C IF (INTGRS(31).EQ.0) THEN DO 1 I = 1,DIMEN REALS(POS+I) = Y(I)+(1D0/5D0)*KI(I,1)*H 1 CONTINUE ENDIF C C Evaluate the lag functions at Tn+H/5 and then the C derivative estimates and update diagnostic C T = REALS(9)+(1D0/5D0)*H CALL LAGVAL(INTGRS,REALS,T) IF (INTGRS(31).EQ.1) RETURN INTGRS(60) = 2 CALL DERIV(INTGRS,REALS,T,KI(1,2)) INTGRS(60) = 0 INTGRS(25) = INTGRS(25)+1 IF (INTGRS(31).EQ.1) RETURN C C Calculate Y(Tn+3H/10) for all dimensions of problem C IF (INTGRS(31).EQ.0) THEN DO 2 I = 1,DIMEN REALS(POS+I) = Y(I)+((3D0/40D0)*KI(I,1)+ + (9D0/40D0)*KI(I,2))*H 2 CONTINUE ENDIF C C Evaluate the lag functions at Tn+3H/10 and then the C derivative estimates and update diagnostic C T = REALS(9)+(3D0/10D0)*H CALL LAGVAL(INTGRS,REALS,T) IF (INTGRS(31).EQ.1) RETURN INTGRS(60) = 3 CALL DERIV(INTGRS,REALS,T,KI(1,3)) INTGRS(60) = 0 INTGRS(25) = INTGRS(25)+1 IF (INTGRS(31).EQ.1) RETURN C C Calculate Y(Tn+4H/5) for all dimensions of problem C IF (INTGRS(31).EQ.0) THEN DO 3 I = 1,DIMEN REALS(POS+I) = Y(I)+((44D0/45D0)*KI(I,1)- + (56D0/15D0)*KI(I,2)+(32D0/9D0)*KI(I,3))*H 3 CONTINUE ENDIF C C Evaluate the lag functions at Tn+4H/5 and then the C derivative estimates and update diagnostic C T = REALS(9)+(4D0/5D0)*H CALL LAGVAL(INTGRS,REALS,T) IF (INTGRS(31).EQ.1) RETURN INTGRS(60) = 4 CALL DERIV(INTGRS,REALS,T,KI(1,4)) INTGRS(60) = 0 INTGRS(25) = INTGRS(25)+1 IF (INTGRS(31).EQ.1) RETURN C C Calculate Y(Tn+8H/9) for all dimensions of problem C IF (INTGRS(31).EQ.0) THEN DO 4 I = 1,DIMEN REALS(POS+I) = Y(I)+((19372D0/6561D0)*KI(I,1)- + (25360D0/2187D0)*KI(I,2)+(64448D0/6561D0)*KI(I,3)- + (212D0/729D0)*KI(I,4))*H 4 CONTINUE ENDIF C C Evaluate the lag functions at Tn+8H/9 and then the C derivative estimates and update diagnostic C T = REALS(9)+(8D0/9D0)*H CALL LAGVAL(INTGRS,REALS,T) IF (INTGRS(31).EQ.1) RETURN INTGRS(60) = 5 CALL DERIV(INTGRS,REALS,T,KI(1,5)) INTGRS(60) = 0 INTGRS(25) = INTGRS(25)+1 IF (INTGRS(31).EQ.1) RETURN C C Calculate Y(Tn+H) for all dimensions of problem C IF (INTGRS(31).EQ.0) THEN DO 5 I = 1,DIMEN REALS(POS+I) = Y(I)+((9017D0/3168D0)*KI(I,1)- + (355D0/33D0)*KI(I,2)+(46732D0/5247D0)*KI(I,3)+ + (49D0/176D0)*KI(I,4)-(5103D0/18656D0)*KI(I,5))*H 5 CONTINUE ENDIF C C Evaluate the lag functions at Tn+H and then the C derivative estimates and update diagnostic C C Ensure that if we are solving upto the point `solve upto', C then we do not cross it (especially important if there is C a discontinuity at `solve upto'). C T = REALS(9)+H IF (DABS(T-REALS(10)).LE.2D0*DMAX1(DABS(T),1D0)*REALS(1)) + T = REALS(10)-2D0*DMAX1(DABS(T),1D0)*REALS(1) C CALL LAGVAL(INTGRS,REALS,T) IF (INTGRS(31).EQ.1) RETURN INTGRS(60) = 6 CALL DERIV(INTGRS,REALS,T,KI(1,6)) INTGRS(60) = 0 INTGRS(25) = INTGRS(25)+1 IF (INTGRS(31).EQ.1) RETURN C C Calculate Y(Tn+H) for all dimensions of problem C DO 6 I = 1,DIMEN REALS(POS+I) = Y(I)+((35D0/384D0)*KI(I,1)+ + (500D0/1113D0)*KI(I,3)+(125D0/192D0)*KI(I,4)- + (2187D0/6784D0)*KI(I,5)+(11D0/84D0)*KI(I,6))*H 6 CONTINUE C C Re-evaluate the derivative estimates at Tn+H and update diagnostic C INTGRS(60) = 7 CALL DERIV(INTGRS,REALS,T,KI(1,7)) INTGRS(60) = 0 INTGRS(25) = INTGRS(25)+1 IF (INTGRS(31).EQ.1) RETURN C C Temporarily store Y1, so that the extra stages required for C a fifth order Hermite interpolant may be calculated C DO 7 I = 1,DIMEN REALS(INTGRS(47)+I-1) = REALS(POS+I) 7 CONTINUE C C Calculate fourth order Y(Tn+H/2) for all dimensions of problem C IF (INTGRS(31).EQ.0) THEN DO 8 I = 1,DIMEN REALS(POS+I) = Y(I)+ + ((-33728713D0/104693760D0)*KI(I,1)+2D0*KI(I,2)- + (30167461D0/21674880D0)*KI(I,3)+ + (7739027D0/17448960D0)*KI(I,4)- + (19162737D0/123305984D0)*KI(I,5)- + (26949D0/363520D0)*KI(I,7))*H 8 CONTINUE ENDIF C C Evaluate the lag functions at Tn+H/2 and then the C derivative estimates and update diagnostic C T = REALS(9)+H/2D0 CALL LAGVAL(INTGRS,REALS,T) IF (INTGRS(31).EQ.1) RETURN INTGRS(60) = 8 CALL DERIV(INTGRS,REALS,T,KI(1,8)) INTGRS(60) = 0 INTGRS(25) = INTGRS(25)+1 IF (INTGRS(31).EQ.1) RETURN C C Calculate fifth order Y(Tn+H/2) for all dimensions of problem C DO 9 I = 1,DIMEN REALS(POS+I) = Y(I)+ + ((7157D0/75776D0)*KI(I,1)+(70925D0/164724D0)*KI(I,3)+ + (10825D0/113664D0)*KI(I,4)- + (220887D0/4016128D0)*KI(I,5)+ + (80069D0/3530688D0)*KI(I,6)-(107D0/5254D0)*KI(I,7)- + (5D0/74D0)*KI(I,8))*H 9 CONTINUE C C Store the fifth order solution at Tn+H/2 for use later C DO 10 I = 1,DIMEN REALS(INTGRS(46)+I-1) = REALS(POS+I) 10 CONTINUE C C Re-evaluate the derivative estimates at Tn+H/2 and update C diagnostic C INTGRS(60) = 9 CALL DERIV(INTGRS,REALS,T,KI(1,9)) INTGRS(60) = 0 INTGRS(25) = INTGRS(25)+1 C C Restore the solution at the endpoint C DO 11 I = 1,DIMEN REALS(POS+I) = REALS(INTGRS(47)+I-1) 11 CONTINUE C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE LAGVAL(INTGRS,REALS,AT) C DOUBLE PRECISION LAG C INTEGER INTGRS(*),LOOP,POS DOUBLE PRECISION AT,REALS(*) C EXTERNAL LAG C C Routine to evaluate complete set of lag functions C at the point AT. Order of evaluation is IMPORTANT C in the case of state-dependent DDEs. POS is the C position in REALS of the lag values. C POS = INTGRS(45)-1 C C Calculate the ODE zero lag function first C REALS(POS) = AT C DO 1 LOOP = 1,INTGRS(14) REALS(POS+LOOP) = LAG(INTGRS,REALS,LOOP,AT) C C Exit LAGVAL immediately if a state dependent C lag function requires a solution value from C the current interval and iterative refinement C has been specified. C IF (INTGRS(31).EQ.1) RETURN 1 CONTINUE C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE LBOUND(INTGRS,REALS,BOUND) C INTEGER INTGRS(*) DOUBLE PRECISION BOUND,REALS(*) C EXTERNAL ERRHAN C C Routine to specify an upper bound on the size of the lag(s) C C Ensure that the bound is not negative C IF (BOUND.LT.0D0) THEN CALL ERRHAN(INTGRS,'Upper bound for lag(s) is negative.', + 'LBOUND.',2) IF (INTGRS(15).GE.2) WRITE (*,*) 'Current upper bound retained' ELSE REALS(18) = BOUND ENDIF C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE LIST C C Routine to list user callable and supplied routines C WRITE (*,*) ' List of User callable routines' WRITE (*,*) ' ------------------------------' WRITE (*,*) WRITE (*,*) 'CONTY(INTGRS,REALS,solution,order,position)' WRITE (*,*) +' Routine to specify the position and order of discontinuities' WRITE (*,*) +' in the components of the initial function, the derivative' WRITE (*,*) ' function and at the initial point.' WRITE (*,*) WRITE (*,*) 'CUTOFF(INTGRS,REALS,order of discontinuity)' WRITE (*,*) +' Routine to specify the maximum order upto which' WRITE (*,*) ' discontinuities are to be tracked.' WRITE (*,*) WRITE (*,*) 'DELAY(INTGRS,REALS,lag function,solution)' WRITE (*,*) +' Function to calculate the specified solution component at' WRITE (*,*) +' the position given by the specified lag function value.' WRITE (*,*) WRITE (*,*) 'DELAYS(INTGRS,REALS,lag functions)' WRITE (*,*) +' Routine to specify the maximum number of lag functions' WRITE (*,*) ' (excluding Volterra lag functions).' WRITE (*,*) WRITE (*,*) 'DIMEN(INTGRS,REALS,dimension of problem)' WRITE (*,*) +' Routine to specify the dimension of the system of functional' WRITE (*,*) ' differential equations.' WRITE (*,*) WRITE (*,*) 'ERROR(INTGRS,REALS,error estimate type,tolerance)' WRITE (*,*) +' Routine to specify what type of error estimate to use and' WRITE (*,*) ' the error-per-step tolerance.' WRITE (*,*) WRITE (*,*) 'EXTRAP(INTGRS,REALS,mode)' WRITE (*,*) +' Routine to specify the extrapolation mode to be used. If' WRITE (*,*) ' `mode`=1, no extrapolation is allowed. If' WRITE (*,*) ' `mode`=2, only extrapolation in the current' WRITE (*,*) ' interval is allowed. If `mode`=3, extrapolation' WRITE (*,*) ' beyond the current interval allowed.' WRITE (*,*) WRITE (*,*) 'FIXSTP(INTGRS,REALS,stepsize)' WRITE (*,*) +' Routine to specify that the constant stepsize strategy is' WRITE (*,*) ' to be used and the stepsize to use.' WRITE (*,*) WRITE (*,*) 'HARD(INTGRS,REALS,from,to)' WRITE (*,*) +' Routine to specify a hard discontinuity mapping: propagating' WRITE (*,*) +' the discontinuity from the `from` solution component to' WRITE (*,*) ' the `to` component with one order of smoothing.' WRITE (*,*) WRITE (*,*) 'HSTART(INTGRS,REALS,stepsize)' WRITE (*,*) ' Routine to specify the initial stepsize.' WRITE (*,*) WRITE (*,*) 'INIT(INTGRS,REALS,intgrs size,reals size)' WRITE (*,*) ' Initialization routine: must be the first call.' WRITE (*,*) WRITE (*,*) 'LBOUND(INTGRS,REALS,upper bound)' WRITE (*,*) +' Routine to specify an upper bound on the size of the lag(s).' WRITE (*,*) WRITE (*,*) 'MAXSTP(INTGRS,REALS,stepsize)' WRITE (*,*) ' Routine to specify the maximum stepsize.' WRITE (*,*) WRITE (*,*) 'NEUTRL(INTGRS,REALS,lag function,derivative)' WRITE (*,*) +' Function to calculate the specified derivative component at' WRITE (*,*) +' the position given by the specified lag function value.' WRITE (*,*) WRITE (*,*) 'NHARD(INTGRS,REALS,from,to)' WRITE (*,*) +' Routine to specify a hard discontinuity mapping: propagating' WRITE (*,*) +' the discontinuity from the `from` solution component to the ' WRITE (*,*) ' `to` component with no smoothing.' WRITE (*,*) WRITE (*,*) 'NSOFT(INTGRS,REALS,from,to,lag function)' WRITE (*,*) +' Routine to specify a soft discontinuity mapping: propagating' WRITE (*,*) +' the discontinuity from the `from` solution component to the' WRITE (*,*) +' `to` component linked by the specified lag function, with' WRITE (*,*) ' no smoothing.' WRITE (*,*) WRITE (*,*) 'PRNTCO(INTGRS,REALS)' WRITE (*,*) ' Routine to print out all the continuity records.' WRITE (*,*) WRITE (*,*) 'PRNTDY(INTGRS,REALS,derivative,from,to,steps)' WRITE (*,*) +' Routine to print out the specified derivative from `from`' WRITE (*,*) ' to `to` in `steps` intervals.' WRITE (*,*) WRITE (*,*) 'PRNTTI(INTGRS,REALS)' WRITE (*,*) ' Routine to print out all the stored abscissae.' WRITE (*,*) WRITE (*,*) 'PRNTY(INTGRS,REALS,solution,from,to,steps)' WRITE (*,*) +' Routine to print out the specified solution from `from`' WRITE (*,*) ' to `to` in `steps` intervals.' WRITE (*,*) WRITE (*,*) 'RANGE(INTGRS,REALS,initial point,end point)' WRITE (*,*) ' Routine to specify the range of integration.' WRITE (*,*) WRITE (*,*) 'REFINE(INTGRS,REALS)' WRITE (*,*) +' Routine to specify that iterative refinement is to be used' WRITE (*,*) ' in preference to extrapolation.' WRITE (*,*) WRITE (*,*) 'REPORT(INTGRS,REALS,level)' WRITE (*,*) +' Routine to specify the level of error reporting/handling to' WRITE (*,*) ' be used.' WRITE (*,*) WRITE (*,*) 'SOFT(INTGRS,REALS,from,to,lag function)' WRITE (*,*) +' Routine to specify a soft discontinuity mapping: propagating' WRITE (*,*) +' the discontinuity from the `from` solution component to the' WRITE (*,*) +' `to` component linked by the specified lag function, with' WRITE (*,*) ' one order of smoothing.' WRITE (*,*) WRITE (*,*) 'SOLNYI(INTGRS,REALS,solution,at)' WRITE (*,*) +' Function to calculate the specified solution at `at`.' WRITE (*,*) WRITE (*,*) 'SPACE(INTGRS,REALS,distance)' WRITE (*,*) +' Routine to specify the distance between discontinuities' WRITE (*,*) +' before they are considered to be distinct.' WRITE (*,*) WRITE (*,*) 'STATS(INTGRS,REALS,level)' WRITE (*,*) +' Routine to specify the level of diagnostics to be generated' WRITE (*,*) ' during and after running.' WRITE (*,*) WRITE (*,*) 'STORCY(INTGRS,REALS,spaces)' WRITE (*,*) +' Routine to specify the amount of storage to be used for' WRITE (*,*) ' storing continuity records.' WRITE (*,*) WRITE (*,*) 'STORHD(INTGRS,REALS,spaces)' WRITE (*,*) +' Routine to specify the amount of storage to be used for' WRITE (*,*) ' storing hard discontinuity mappings.' WRITE (*,*) WRITE (*,*) 'STORST(INTGRS,REALS,spaces)' WRITE (*,*) +' Routine to specify the amount of storage to be used for' WRITE (*,*) ' storing soft discontinuity mappings.' WRITE (*,*) WRITE (*,*) 'TRACK(INTGRS,REALS,level)' WRITE (*,*) +' Routine to specify the level of discontinuity tracking.' WRITE (*,*) WRITE (*,*) 'VJUMP(INTGRS,REALS)' WRITE (*,*) +' Routine to specify that a kernel of an integral equation' WRITE (*,*) ' has poor continuity.' WRITE (*,*) WRITE (*,*) 'VLAG(INTGRS,REALS,lagfn,lag value)' WRITE (*,*) +' Routine to store `lag value` as the value of the `lagfn`th' WRITE (*,*) ' Volterra lag function.' WRITE (*,*) WRITE (*,*) 'VOLTER(INTGRS,REALS,Volterra terms,lag functions)' WRITE (*,*) +' Routine to specify the number of accelerated Volterra terms' WRITE (*,*) ' and the number of kernel lag functions.' WRITE (*,*) WRITE (*,*) 'VOLTRA(INTGRS,REALS,T,lagfn1,lagfn2,kernfn)' WRITE (*,*) +' Function to evaluate the integral with lower limit specified' WRITE (*,*) +' by `lagfn1` and upper limit specified by `lagfn2` using the' WRITE (*,*) +' `kernfn`th kernel function in KERNEL.' WRITE (*,*) WRITE (*,*) 'YSTOP(INTGRS,REALS,solution,t0,value)' WRITE (*,*) +' Function to find the point at which the specified solution' WRITE (*,*) +' has the value `value` (using starting iterate `t0`).' WRITE (*,*) WRITE (*,*) C WRITE (*,*) ' List of User supplied routines' WRITE (*,*) ' ------------------------------' WRITE (*,*) WRITE (*,*) 'DERIV(INTGRS,REALS,t,f)' WRITE (*,*) +' Routine to calculate all components of the derivative at `t`' WRITE (*,*) ' and return them in the vector `f`.' WRITE (*,*) WRITE (*,*) 'LAG(INTGRS,REALS,lagfn,t)' WRITE (*,*) +' Function to calculate the lag function `lagfn` at `t`' WRITE (*,*) WRITE (*,*) 'YINIT(soln,t)' WRITE (*,*) +' Function to calculate the `soln`th component of the initial' WRITE (*,*) ' function at `t`.' WRITE (*,*) WRITE (*,*) 'FINIT(deriv,t)' WRITE (*,*) +' Function to calculate the `soln`th component of the initial' WRITE (*,*) ' derivative at `t`.' WRITE (*,*) WRITE (*,*) 'KERNEL(INTGRS,REALS,s,t,kernfn)' WRITE (*,*) +' Function to calculate the `kernfn`th kernel of a Volterra' WRITE (*,*) ' term at the points `s` and `t`.' C STOP END C C C ******************************************************************+++++++ C C SUBROUTINE MAPSET(INTGRS,WHERE) C INTEGER INTGRS(*),LOOP CHARACTER*(*) WHERE C EXTERNAL ERRHAN C C Routine to set up the pointers to the HARD delay and neutral C discontinuity networks, and SOFT delay and neutral discontinuity C networks. Also other pointers dependent on the number of delays, C lag functions, etc. C C Pointer to start of HARD discontinuity networks in INTGRS C INTGRS(40) = INTGRS(39)+INTGRS(21)*2 C C Pointer to start of SOFT discontinuity networks in INTGRS C INTGRS(41) = INTGRS(40)+INTGRS(19)*2 C C Pointer to start of individual pointers in INTGRS into cyclic C storage for each lag function and Volterra lag function that C can exist and initialize them to all being 1. Note we allow C an extra storage position at the start of the array for the C zero lag function to correspond to ODEs C INTGRS(42) = INTGRS(41)+INTGRS(20)*3+1 DO 1 LOOP = -1,INTGRS(14)+INTGRS(44)-1 INTGRS(INTGRS(42)+LOOP) = 1 1 CONTINUE C C Pointer to start of 'types of lag function' array in INTGRS C INTGRS(43) = INTGRS(42)+INTGRS(14)+INTGRS(44) C C Check that the size of INTGRS is not too small C IF (INTGRS(43)+INTGRS(14).GT.INTGRS(12)) + CALL ERRHAN(INTGRS, + 'Size of INTGRS is too small.',WHERE,-1) C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE MAXSTP(INTGRS,REALS,STEP) C DOUBLE PRECISION DABS,DMAX1 C INTEGER INTGRS(*) DOUBLE PRECISION REALS(*),STEP C EXTERNAL ERRHAN INTRINSIC DABS,DMAX1 C C Routine for specifying the maximum permitted stepsize C C Ensure that MAXSTP is only called after the start and C endpoint have been given C IF (INTGRS(1).EQ.0) + CALL ERRHAN(INTGRS, + 'Range of integration must be specified first.', + 'MAXSTP.',-1) C C If a fixed stepsize strategy is selected, ignore this call C IF (INTGRS(4).EQ.1) THEN CALL ERRHAN(INTGRS, + 'Fixed stepsize strategy already selected.','MAXSTP.',2) IF (INTGRS(15).GE.2) WRITE (*,*) + 'Request to set maximum permitted stepsize has been ignored' RETURN ENDIF C C Ensure that the maximum stepsize is not too small C IF (STEP.LE.1000D0*DMAX1(1D0,DABS(REALS(7)))*REALS(1)) THEN CALL ERRHAN(INTGRS, + 'Maximum permitted stepsize too small.','MAXSTP.',2) C C Increase maximum permitted stepsize to a `significant' size C STEP = 10000D0*DMAX1(1D0,DABS(REALS(7)))*REALS(1) IF (INTGRS(15).GE.2) WRITE (*,*) + 'The maximum permitted stepsize has been set to ',STEP ENDIF C C Ensure that the maximum permitted stepsize is not too large C IF (REALS(7)+STEP.GT.REALS(8)) THEN CALL ERRHAN(INTGRS, + 'Maximum permitted stepsize too large.','MAXSTP.',2) C C Set maximum permitted stepsize to FINISH-START C STEP = REALS(8)-REALS(7) IF (INTGRS(15).GE.2) WRITE (*,*) + 'The maximum permitted stepsize has been set to ',STEP ENDIF C C Ensure that if the initial stepsize has been set, it does not C violate STEP C IF (INTGRS(2).EQ.1 .AND. REALS(11).GT.STEP) THEN CALL ERRHAN(INTGRS, + 'Initial stepsize exceeds maximum permitted stepsize.', + 'MAXSTP.',2) REALS(11) = STEP IF (INTGRS(15).GE.2) WRITE (*,*) + 'Initial stepsize has been reduced to ',STEP ENDIF C C Flag that maximum stepsize has been specified C INTGRS(3) = 1 C REALS(12) = STEP C RETURN END C C C ******************************************************************+++++++ C C DOUBLE PRECISION FUNCTION NEUTRL(INTGRS,REALS,LAGFN,SOLN) C DOUBLE PRECISION DABS,FINIT,RECALL C INTEGER INTGRS(*),LAGFN,SOLN DOUBLE PRECISION LAGVAL,REALS(*) C EXTERNAL ERRHAN,FINIT,RECALL INTRINSIC DABS C C Routine to retrieve the SOLN'th component of the derivative C solution at the point indicated by the LAG'th lag value. The C initial derivative must be given by the function FINIT C C Ensure that the given lag function and solution exist. C Use kernel number (INTGRS(58)) to determine what type of C lag function we are evaluating. Only check information ONCE C IF (INTGRS(25).EQ.0) THEN C IF (LAGFN.LT.0 .OR. (INTGRS(58).EQ.0 .AND. LAGFN.GT.INTGRS(14)) + .OR. (INTGRS(58).NE.0 .AND. LAGFN.GT.INTGRS(14)+INTGRS(44))) + CALL ERRHAN(INTGRS,'Specified lag function cannot exist.', + 'NEUTRL.',-1) C IF (SOLN.LT.1 .OR. SOLN.GT.INTGRS(13)) + CALL ERRHAN(INTGRS, + 'Specified solution component cannot exist.', + 'NEUTRL.',-1) ENDIF C C Increment the `number of back evaluations' diagnostic C INTGRS(28) = INTGRS(28)+1 C C Calculate the lag-point C LAGVAL = REALS(INTGRS(45)+LAGFN-1) C C If we evaluate a back-derivative `just before' the initial point, C set the `re-evaluate first stage' flag to ON. This is only C (really) necessary when discontinuities are NOT tracked AND C there is a jump in the derivative at the initial point. C IF (INTGRS(60).NE.1 .AND. DABS(REALS(7)-LAGVAL).LT.3D0*REALS(1)) + INTGRS(59) = 1 C C If an ODE term is requested (the zeroth lag function), make C a call to RECALL so that either extrapolation is done or C iterative refinement is signalled accordingly C C Check if the initial derivative is to be called C IF (LAGVAL.LT.REALS(7)) THEN NEUTRL = FINIT(SOLN,LAGVAL) RETURN ENDIF C C Ensure that if iterative refinement has been signalled, we C at least return a valid value for the derivative (the value C at the start of the current interval) C IF (INTGRS(31).EQ.1) THEN C C Do something special in the case of the first interval, since C f(t0) will not have been set up C IF (INTGRS(25).EQ.0) THEN NEUTRL = FINIT(SOLN,LAGVAL) ELSE NEUTRL = REALS(INTGRS(36)+SOLN-1) ENDIF RETURN ENDIF C C Call the actual function which calculates the derivative value C passing TYPE as .FALSE. to indicate that we want to evaluate C a neutral value and not a delay value, and the position of C the storage arrays in REALS. C NEUTRL = RECALL(INTGRS,REALS,.FALSE.,LAGFN,LAGVAL,SOLN, + REALS(INTGRS(34)),INTGRS(13), + REALS(INTGRS(48)),REALS(INTGRS(49)), + REALS(INTGRS(50)),REALS(INTGRS(51)), + REALS(INTGRS(52)),REALS(INTGRS(53)), + REALS(INTGRS(54))) C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE NHARD(INTGRS,REALS,FROM,IN) C INTEGER FROM,IN,INTGRS(*),POS DOUBLE PRECISION REALS(*) C EXTERNAL ERRHAN,MAPSET C C Routine to specify a HARD neutral discontinuity network link C C FROM - Component of solution in which discontinuity already exists C IN - Component of solution in which discontinuity is created C C Ensure that FROM does not exceed the dimension of the problem C IF (FROM.LT.1 .OR. FROM.GT.INTGRS(13)) + CALL ERRHAN(INTGRS, + 'Discontinuity occurs in a solution that cannot exist.', + 'NHARD.',-1) C C Ensure that IN does not exceed the dimension of the problem C IF (IN.LT.1 .OR. IN.GT.INTGRS(13)) + CALL ERRHAN(INTGRS, + 'Discontinuity mapped to a solution that cannot exist.', + 'NHARD.',-1) C IF (FROM.EQ.IN) THEN CALL ERRHAN(INTGRS,'Solution components must be different.', + 'NHARD.',2) IF (INTGRS(15).GE.2) WRITE (*,*) + 'HARD discontinuity network mapping ignored' RETURN ENDIF C C If this is the first call to either CONTY, (N)HARD or (N)SOFT set C up the pointers to start of arrays in INTGRS. After which the number C of Volterra terms and storage allocated to the continuity records C and the HARD and SOFT discontinuity networks cannot be altered. C C INTGRS(8)=0 & INTGRS(9)=0 & INTGRS(10)=0 C IF (INTGRS(8)+INTGRS(9)+INTGRS(10).EQ.0) + CALL MAPSET(INTGRS,'NHARD.') C C Check that there is still space left for mapping C IF (INTGRS(8).EQ.INTGRS(19)) THEN CALL ERRHAN(INTGRS, + 'Insufficient space for HARD discontinuity network.', + 'NHARD.',2) IF (INTGRS(15).GE.2) + WRITE (*,*) 'HARD discontinuity network mapping ignored' RETURN ENDIF C C Save the HARD discontinuity mapping: IN, FROM C Distinguish delay and neutral mappings by sign of IN C POS = INTGRS(40)+2*INTGRS(8) INTGRS(POS) = FROM INTGRS(POS+1) = -IN C C Update the number of HARD discontinuity mappings set C INTGRS(8) = INTGRS(8)+1 C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE NSOFT(INTGRS,REALS,FROM,IN,LINK) C INTEGER FROM,IN,INTGRS(*),LINK,POS DOUBLE PRECISION REALS(*) C EXTERNAL ERRHAN,MAPSET C C Routine to specify a SOFT neutral discontinuity network link C C FROM - Component of solution in which discontinuity already exists C IN - Component of solution in which discontinuity is created C LINK - Lag term which links IN and FROM solution components C C Ensure that FROM does not exceed the dimension of the problem C IF (FROM.LT.1 .OR. FROM.GT.INTGRS(13)) + CALL ERRHAN(INTGRS, + 'Discontinuity occurs in a solution that cannot exist.', + 'NSOFT.',-1) C C Ensure that IN does not exceed the dimension of the problem C IF (IN.LT.1 .OR. IN.GT.INTGRS(13)) + CALL ERRHAN(INTGRS, + 'Discontinuity mapped to a solution that cannot exist.', + 'NSOFT.',-1) C C Ensure that LINK corresponds to a lag function which can exist C IF (LINK.LT.1 .OR. LINK.GT.INTGRS(14)) + CALL ERRHAN(INTGRS, + 'Linked lag function cannot exist.','NSOFT.',-1) C C If this is the first call to either CONTY, (N)HARD or (N)SOFT set C up the pointers to start of arrays in INTGRS. After which the number C of Volterra terms and storage allocated to the continuity records C and the HARD and SOFT discontinuity networks cannot be altered. C C INTGRS(8)=0 & INTGRS(9)=0 & INTGRS(10)=0 C IF (INTGRS(8)+INTGRS(9)+INTGRS(10).EQ.0) + CALL MAPSET(INTGRS,'NSOFT.') C C Check that there is still space left for mapping C IF (INTGRS(9).EQ.INTGRS(20)) THEN CALL ERRHAN(INTGRS, + 'Insufficient space for SOFT discontinuity network.', + 'NSOFT.',2) IF (INTGRS(15).GE.2) + WRITE (*,*) 'SOFT discontinuity network mapping ignored' RETURN ENDIF C C Save the SOFT discontinuity mapping: IN, FROM, LINK C Distinguish delay and neutral mapping by sign of IN C POS = INTGRS(41)+3*INTGRS(9) INTGRS(POS) = FROM INTGRS(POS+1) = -IN INTGRS(POS+2) = LINK C C Update the number of SOFT discontinuity mappings set C INTGRS(9) = INTGRS(9)+1 C RETURN END C C C ******************************************************************+++++++ C C DOUBLE PRECISION FUNCTION NXTDSC(INTGRS,REALS) C DOUBLE PRECISION SCANJUMP C INTEGER INTGRS(*),LOOP,POS DOUBLE PRECISION LIMIT,REALS(*) LOGICAL CHANGE C EXTERNAL DOHARD,DOTRCK,ERRHAN,SCANJUMP C C Function to return the position of the next discontinuity C that RETARD will encounter C C If the continuity record is full, skip any further tracking C IF (INTGRS(10).GE.INTGRS(21)) THEN C C Only issue the warning ONCE! C IF (INTGRS(10).EQ.INTGRS(21)) THEN CALL ERRHAN(INTGRS,'Continuity record is full.','NXTDSC.',2) INTGRS(21) = INTGRS(21)-1 ENDIF GOTO 4 ENDIF C C Set the furthest point we can `solve upto', which may be reduced C by state dependent discontinuities which require extrapolation C (and thus are not stored as continuity records) C LIMIT = REALS(8) C C Calculate the next set of discontinuities for all SOFT C discontinuity network mappings C C Get the start of the SOFT network in INTGRS C POS = INTGRS(41) C C Repeat for all the SOFT discontinuity mappings set C DO 1 LOOP = 1,INTGRS(9) C C Obey the options set C IF (INTGRS(17).EQ.0) THEN C C Track only order 0 neutral discontinuities if CUTOFF = 0 C IF (INTGRS(POS+1).LT.0) CALL DOTRCK(INTGRS,REALS,POS,LIMIT) ELSE C C Check if tracking only state independent delays selected C IF (INTGRS(5).EQ.1) THEN C C If selected, only track state independent delays C IF (INTGRS(INTGRS(43)+INTGRS(POS+2)-1).EQ.0) + CALL DOTRCK(INTGRS,REALS,POS,LIMIT) ELSE C C Tracking all discontinuities selected C CALL DOTRCK(INTGRS,REALS,POS,LIMIT) ENDIF ENDIF C C Update the pointer to the SOFT discontinuity mappings C POS = POS+3 C 1 CONTINUE C C Track the discontinuities mapped by the HARD mappings C C Set CHANGE to indicate whether any changes are made to C the continuity records when tracking the HARD mappings C 2 CHANGE = .FALSE. C C Repeat for all HARD discontinuity mappings set C DO 3 LOOP = 1,INTGRS(8) CALL DOHARD(INTGRS,REALS,LOOP,CHANGE) 3 CONTINUE C C Repeat the tracking of HARD discontinuities until no change C IF (CHANGE) GOTO 2 C C Find the `next' discontinuity by scanning the continuity records C 4 NXTDSC = SCANJUMP(INTGRS,REALS,LIMIT) C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE ONENRM(INTGRS,REALS,DIMEN,ERRORI,CTSPLY) C DOUBLE PRECISION DABS C INTEGER DIMEN,INTGRS(*),LOOP,SOLN DOUBLE PRECISION AREA,OLDVAL,POLY(0:3),ROOTS(6),T,VALUE DOUBLE PRECISION CTSPLY(DIMEN,5),ERRORI(*),REALS(*) C EXTERNAL FINDRT INTRINSIC DABS C C Routine to calculate the one-norm of the error polynomial C C Repeat for all components of the solution C DO 4 SOLN = 1,DIMEN C C Transfer the correct component of the error polynomial C to POLY to allow for deflation when searching for roots C DO 1 LOOP = 0,3 POLY(LOOP) = CTSPLY(SOLN,LOOP+2) 1 CONTINUE C C Find the roots in the interval C CALL FINDRT(REALS,POLY,ROOTS) C C Find the zero root in the ordered sequence C LOOP = 0 2 LOOP = LOOP+1 IF (DABS(ROOTS(LOOP)).GT.REALS(1)) GOTO 2 C C Calculate the absolute area from 0 to 1 C AREA = 0D0 OLDVAL = 0D0 C C The absolute integration loop C 3 LOOP = LOOP+1 T = ROOTS(LOOP) VALUE = T*T*T*((1D0/3D0)*CTSPLY(SOLN,2)+T* + ((1D0/4D0)*CTSPLY(SOLN,3)+T*((1D0/5D0)*CTSPLY(SOLN,4)+ + T*(1D0/6D0)*CTSPLY(SOLN,5)))) AREA = AREA+DABS(VALUE-OLDVAL) OLDVAL = VALUE C C Repeat until the end of the continuous interval T = 1 C IF (DABS(ROOTS(LOOP)-1D0).GT.REALS(1)) GOTO 3 C C Store this component of the error estimate and repeat C for all components of the problem C ERRORI(SOLN) = AREA*REALS(11) 4 CONTINUE C RETURN END C C C ******************************************************************+++++++ C C DOUBLE PRECISION FUNCTION PLYVAL(INTGRS,DIMEN,COMPNT,CTSPLY,AT) C DOUBLE PRECISION DABS C INTEGER DIMEN,COMPNT,INTGRS(*) DOUBLE PRECISION AT,CTSPLY(DIMEN,5) C INTRINSIC DABS C C Function to evaluate the absolute value of the continuous error C polynomial corresponding to the COMPNTth solution at the scaled C point AT which lies in [0,1] (with no constant or linear term) C PLYVAL = DABS((((CTSPLY(COMPNT,5)*AT+CTSPLY(COMPNT,4))*AT+ + CTSPLY(COMPNT,3))*AT+CTSPLY(COMPNT,2))*AT*AT) C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE PRNTCO(INTGRS,REALS) C INTEGER INTAT,INTGRS(*),LOOP DOUBLE PRECISION REALS(*) C C Routine to print out a full list of the continuity C records calculated so far C C Calculate the position in INTGRS of the first continuity record C INTAT = INTGRS(39) C WRITE (*,*) 'Component Smoothness Position' DO 1 LOOP = 1,INTGRS(10) WRITE (*,9999) INTGRS(INTAT),INTGRS(INTAT+1),' ', + REALS(INTGRS(55)+LOOP-1) INTAT = INTAT+2 1 CONTINUE WRITE (*,*) C 9999 FORMAT (I6,I11,A4,G25.17) C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE PRNTDY(INTGRS,REALS,COMPNT,FROM,TO,STEPS) C DOUBLE PRECISION DBLE,NEUTRL C INTEGER COMPNT,INTGRS(*),LOOP,STEPS DOUBLE PRECISION FROM,OLDLAG,REALS(*),SIZES,T,TFIRST,TLAST,TO C EXTERNAL ERRHAN,NEUTRL INTRINSIC DBLE C C Routine to print out the COMPNT'th component derivative values C (still held in the cyclic storage) in STEPS intervals C C Ensure that the component specified is valid C IF (COMPNT.LT.1 .OR. COMPNT.GT.INTGRS(13)) + CALL ERRHAN(INTGRS, + 'Invalid solution component specified.', + 'PRNTDY.',-1) C C Ensure that a valid interval has been given C IF (FROM.GE.TO) + CALL ERRHAN(INTGRS, + 'An invalid interval has been specified.', + 'PRNTDY.',-1) C C Ensure the number of intervals specified is valid C IF (STEPS.LT.1) THEN CALL ERRHAN(INTGRS, + 'Invalid number of derivative intervals.', + 'PRNTDY.',2) IF (INTGRS(15).GE.2) + WRITE (*,*) 'Default number of 20 intervals assumed' STEPS = 20 ENDIF C C Calculate start and end of history queue and the size of each interval C LOOP = 0 SIZES = (TO-FROM)/DBLE(STEPS) TFIRST = REALS(INTGRS(48)+INTGRS(29)-1) TLAST = REALS(3)+DMAX1(DABS(REALS(3)),1D0)*REALS(1) C C Check if FROM is after the initial point, but before the first C interval still retained in the history queue. If it is print appropriate C message and find the first abscissae which is in the history queue C IF (FROM.GE.REALS(7) .AND. FROM.LT.TFIRST) THEN PRINT *,'First interval in the history queue starts at ',TFIRST 1 LOOP = LOOP+1 IF (FROM+SIZES*DBLE(LOOP).LT.TFIRST) GOTO 1 ENDIF C C Print out the data, using the first lag function which must be preserved C OLDLAG = REALS(INTGRS(45)) C 2 T = FROM+SIZES*DBLE(LOOP) LOOP = LOOP+1 IF (T.LT.REALS(7) .OR. (T.GE.TFIRST .AND. T.LE.TLAST)) THEN C C Print out the appropriate abscissae and solution at that point C and repeat if necessary C REALS(INTGRS(45)) = T WRITE (*,9998) T,NEUTRL(INTGRS,REALS,1,COMPNT) ELSE WRITE (*,9999) T,' Derivative is not in history queue' ENDIF IF (LOOP.LE.STEPS) GOTO 2 WRITE (*,*) C C Restore the ODE lag function C REALS(INTGRS(45)) = OLDLAG C 9998 FORMAT (2G25.17) 9999 FORMAT (G25.17,A35) C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE PRNTTI(INTGRS,REALS) C INTEGER MOD C INTEGER INTGRS(*),LOOP DOUBLE PRECISION REALS(*) C INTRINSIC MOD C C Routine to print out all abscissae still held in the C cyclic storage array C C Get the position of the first element C LOOP = INTGRS(29)-1 C C Repeat until all the data has been printed C 1 WRITE (*,9999) REALS(INTGRS(48)+LOOP) LOOP = MOD(LOOP,INTGRS(33))+1 IF (LOOP.NE.MOD(INTGRS(30),INTGRS(33))) GOTO 1 WRITE (*,*) C 9999 FORMAT (G25.17) C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE PRNTY(INTGRS,REALS,COMPNT,FROM,TO,STEPS) C DOUBLE PRECISION DELAY,DBLE C INTEGER COMPNT,INTGRS(*),LOOP,STEPS DOUBLE PRECISION FROM,OLDLAG,REALS(*),SIZES,T,TFIRST,TLAST,TO C EXTERNAL DELAY,ERRHAN INTRINSIC DBLE C C Routine to print out the COMPNT'th component solution values C (still held in the cyclic storage) in STEPS intervals C C Ensure that the component specified is valid C IF (COMPNT.LT.1 .OR. COMPNT.GT.INTGRS(13)) + CALL ERRHAN(INTGRS, + 'Invalid solution component specified.', + 'PRNTY.',-1) C C Ensure that a valid interval has been given C IF (FROM.GE.TO) + CALL ERRHAN(INTGRS, + 'An invalid interval has been specified.', + 'PRNTY.',-1) C C Ensure the number of intervals specified is valid C IF (STEPS.LT.1) THEN CALL ERRHAN(INTGRS, + 'Invalid number of solution intervals.', + 'PRNTY.',2) IF (INTGRS(15).GE.2) + WRITE (*,*) 'Default number of 20 intervals assumed' STEPS = 20 ENDIF C C Calculate start and end of history queue and the size of each interval C LOOP = 0 SIZES = (TO-FROM)/DBLE(STEPS) TFIRST = REALS(INTGRS(48)+INTGRS(29)-1) TLAST = REALS(3)+DMAX1(DABS(REALS(3)),1D0)*REALS(1) C C Check if FROM is after the initial point, but before the first C interval still retained in the history queue. If it is print appropriate C message and find the first abscissae which is in the history queue C IF (FROM.GE.REALS(7) .AND. FROM.LT.TFIRST) THEN PRINT *,'First interval in the history queue starts at ',TFIRST 1 LOOP = LOOP+1 IF (FROM+SIZES*DBLE(LOOP).LT.TFIRST) GOTO 1 ENDIF C C Print out the data, using the first lag function which must be preserved C OLDLAG = REALS(INTGRS(45)) C 2 T = FROM+SIZES*DBLE(LOOP) LOOP = LOOP+1 IF (T.LT.REALS(7) .OR. (T.GE.TFIRST .AND. T.LE.TLAST)) THEN C C Print out the appropriate abscissae and solution at that point C and repeat if necessary C REALS(INTGRS(45)) = T WRITE (*,9998) T,DELAY(INTGRS,REALS,1,COMPNT) ELSE WRITE (*,9999) T,' Solution is not in history queue' ENDIF IF (LOOP.LE.STEPS) GOTO 2 WRITE (*,*) C C Restore the ODE lag function C REALS(INTGRS(45)) = OLDLAG C 9998 FORMAT (2G25.17) 9999 FORMAT (G25.17,A33) C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE RANGE(INTGRS,REALS,START,FINISH) C DOUBLE PRECISION DABS C INTEGER INTGRS(*) DOUBLE PRECISION FINISH,REALS(*),START C EXTERNAL ERRHAN INTRINSIC DABS C C Routine for specifying the interval of integration C C Do not allow the range of integration to be reset C IF (INTGRS(1).EQ.1) THEN CALL ERRHAN(INTGRS, + 'The range of integration cannot be changed.', + 'RANGE.',2) IF (INTGRS(15).GE.2) + WRITE (*,*) 'Current range of integration retained' RETURN ENDIF C C Set last point stored, initial point and current point C REALS(3) = START REALS(7) = START REALS(9) = START C C Do not allow reverse direction of integration C IF (START.GT.FINISH) + CALL ERRHAN(INTGRS, + 'Direction of integration is incorrect.','RANGE.',-1) IF (DABS(START-FINISH).LT.REALS(1)) + CALL ERRHAN(INTGRS, + 'Range of integration is zero.','RANGE.',-1) C REALS(8) = FINISH C C Flag that the start and endpoint have been specified C INTGRS(1) = 1 C RETURN END C C C ******************************************************************+++++++ C C DOUBLE PRECISION FUNCTION RECALL(INTGRS,REALS,TYPE,LAGFN,LAGVAL, + SOLN,Y,DIMEN,TSTOR,YSTOR,DEG1,DEG2,DEG3,DEG4,DEG5) C INTEGER MOD DOUBLE PRECISION DABS,DMAX1 C INTEGER ACTUAL,DIMEN,LAGFN,LENGTH,POINTR,SOLN,INTGRS(*) DOUBLE PRECISION LAGVAL,REALS(*),SCALED,STEP,Y(*) DOUBLE PRECISION TSTOR(*),YSTOR(DIMEN,*) DOUBLE PRECISION DEG1(DIMEN,*),DEG2(DIMEN,*),DEG3(DIMEN,*) DOUBLE PRECISION DEG4(DIMEN,*),DEG5(DIMEN,*) LOGICAL TYPE C EXTERNAL ERRHAN INTRINSIC DABS,DMAX1,MOD C C Function to retrieve back solution and derivative values C TYPE = .TRUE. means retrieve a solution value C TYPE = .FALSE. means retrieve a derivative value C LAGFN indexes the lag function value in the array REALS C SOLN corresponds to the component of the solution that we require C C See if we require a value beyond the stored approximants (taking C account of rounding error), or if no approximants at all have been C stored in which case this is an initial value problem C IF (LAGVAL.GT.REALS(19) .OR. INTGRS(30).EQ.0) THEN C C If we were not in iterative refinement mode when it is C specified and we are not tracking state dependent C discontinuities, or its an initial value problem, return the C value at the start of the interval C C INTGRS(31)=0 & INTGRS(7)=0 & INTGRS(57)=0 C IF (INTGRS(31)+INTGRS(7)+INTGRS(57).EQ.0 + .OR. (INTGRS(30).EQ.0 .AND. INTGRS(7).NE.1)) THEN C C Return the solution or derivative value at the start of C the current interval as appropriate C IF (TYPE) THEN RECALL = Y(SOLN) ELSE RECALL = REALS(INTGRS(36)+SOLN-1) ENDIF C C If not in iterative refinement mode, signal it C IF (INTGRS(31).EQ.0) INTGRS(31) = 1 RETURN ENDIF C C If we are not doing the first stage of the iterative refinement C and we are not tracking state dependent discontinuities, but we C still require a function value beyond the current interval then C this is an ADVANCED equation. However do not signal an error if C extrapolation mode 3 is selected. C IF (((INTGRS(31).NE.2 .AND. INTGRS(7).EQ.0) .OR. + (INTGRS(7).EQ.1 .AND. LAGVAL.GT.REALS(3)) .OR. + (INTGRS(7).EQ.2 .AND. LAGVAL.GT.REALS(3)+REALS(11))) + .AND. INTGRS(57).EQ.0) THEN IF (INTGRS(15).GE.1) + WRITE (*,*) 'Requested ',LAGVAL,' beyond ',REALS(3) CALL ERRHAN(INTGRS,'The equation has an advanced delay.', + 'SOLVE.',-1) ENDIF ENDIF C C Check to see it the approximation data required has been C overwritten - ie the cyclic storage is FULL C IF (LAGVAL.LT.TSTOR(INTGRS(29))) + CALL ERRHAN(INTGRS,'The cyclic storage is full.', + 'SOLVE.',-1) C C Otherwise the approximation data exists or we can extrapolate C from the last available approximant so get the length of the C actual and maximum cyclic storage structure C LENGTH = INTGRS(33) ACTUAL = MOD(INTGRS(30)-1,LENGTH)+1 C C Get the LAGFN's corresponding pointer into REALS C POINTR = INTGRS(INTGRS(42)+LAGFN-1) C C If LAGVAL is after the current interval pointed to by POINTR C increase POINTR in a cyclic fashion until we reach the correct C or the current interval C 1 IF (LAGVAL.GT.TSTOR(POINTR)) THEN IF (POINTR.NE.ACTUAL) THEN POINTR = MOD(POINTR,LENGTH)+1 GOTO 1 ENDIF ENDIF C C If LAGVAL is before the current interval pointed to by POINTR C decrease POINTR in a cyclic fashion until we reach the correct C or the oldest interval C 2 IF (LAGVAL.LT.TSTOR(POINTR)) THEN IF (POINTR.NE.INTGRS(29)) THEN POINTR = MOD(POINTR+LENGTH-2,LENGTH)+1 GOTO 2 ENDIF ENDIF C C Save the value of POINTR for the next call of this lag function C INTGRS(INTGRS(42)+LAGFN-1) = POINTR C C Calculate the stepsize corresponding to this interval. If C LAGVAL lies in or beyond the last interval stored, the length C of this approximant interval is REALS(3)-TSTOR(POINTR) C IF (POINTR.EQ.ACTUAL) THEN C C Get the length of the last interval stored C STEP = REALS(3)-TSTOR(POINTR) ELSE C C Get the length of a normal stored approximant C STEP = TSTOR(MOD(POINTR,LENGTH)+1)-TSTOR(POINTR) ENDIF C C Calculate the scaled stepsize (in the range 0 to 1) C IF (STEP.LT.REALS(1)) THEN SCALED = 0D0 ELSE SCALED = (LAGVAL-TSTOR(POINTR))/STEP ENDIF C C If this is extrapolation C 3 IF (SCALED.GT.1D0) THEN C C Ensure that the extrapolation is not too `excessive', caused C by the last accepted stepsize being `small' due to the presence C of a state dependent discontinuity whose position required C extrapolation to `find' it. C IF (SCALED.GT.15D0 .AND. POINTR.NE.INTGRS(29)) THEN C C If the extrapolation is excessive (2*maximum stepsize increase) C try extrapolating from a previous approximant. Since the code C maintains an in-size stepsize when approaching discontinuities C and we do not extrapolate when solving the first interval, a C `reasonable' sized interval MUST exist. C (NOTE: We should not extrapolate over low order discontinuities) C POINTR = MOD(POINTR+LENGTH-2,LENGTH)+1 STEP = TSTOR(MOD(POINTR,LENGTH)+1)-TSTOR(POINTR) SCALED = (LAGVAL-TSTOR(POINTR))/STEP GOTO 3 ENDIF C C If we are extrapolating so that we can track state dependent C discontinuities, then change INTGRS(57) to show that we needed C extrapolation so that the position of the discontinuity may be C flagged as needing further refinement C IF (INTGRS(57).GT.0) THEN INTGRS(57) = 2 ELSE C C If extrapolation selected C IF (INTGRS(7).NE.0) THEN C C If level of error reporting = 1, report extrapolation C IF (INTGRS(15).EQ.1) WRITE (*,*) + 'Extrapolate from [',TSTOR(POINTR),',',REALS(3), + '] to ',LAGVAL C C Update the maximum amount of extrapolation `diagnostic' C INTGRS(27) = INTGRS(27)+1 IF (SCALED.GT.REALS(17)) REALS(17) = SCALED ENDIF ENDIF ENDIF C C Calculate the solution or the derivative as requested C IF (TYPE) THEN C C Calculate the back solution C RECALL = YSTOR(SOLN,POINTR)+SCALED*(DEG1(SOLN,POINTR)+ + SCALED*(DEG2(SOLN,POINTR)+SCALED*(DEG3(SOLN,POINTR)+ + SCALED*(DEG4(SOLN,POINTR)+SCALED*DEG5(SOLN,POINTR))))) ELSE C C Calculate the back derivative C RECALL = (DEG1(SOLN,POINTR)+SCALED*(2D0*DEG2(SOLN,POINTR)+ + SCALED*(3D0*DEG3(SOLN,POINTR)+SCALED*(4D0*DEG4(SOLN,POINTR)+ + 5D0*SCALED*DEG5(SOLN,POINTR)))))/STEP ENDIF C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE REFINE(INTGRS,REALS) C INTEGER INTGRS(*) DOUBLE PRECISION REALS(*) C EXTERNAL ERRHAN C C Routine for specifying that iterative refinement instead C of extrapolation is to be carried out C C Report that iterative refinement is to be carried out C IF (INTGRS(15).GE.2) WRITE (*,*) 'Iterative refinement selected' C IF (INTGRS(5).LT.2 .OR. INTGRS(17).LT.5) CALL ERRHAN(INTGRS, + 'Refinement needs all discontinuities upto order 5 tracked.', + 'REFINE.',2) C INTGRS(7) = 0 C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE REPORT(INTGRS,REALS,LEVEL) C INTEGER INTGRS(*),LEVEL DOUBLE PRECISION REALS(*) C EXTERNAL ERRHAN C C Routine for specifying level of error control C C LEVEL = 0 - No errors reports, automatic correction where possible C 1 - Reporting of extrapolation/iterative refinement only C 2 - Report all errors, automatic correction where possible C 3 - Report all errors and STOP (default) C IF (LEVEL.LT.0 .OR. LEVEL.GT.3) + CALL ERRHAN(INTGRS, + 'Invalid level of error reporting.','REPORT.',-1) C INTGRS(15) = LEVEL C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE RETARD(INTGRS,REALS) C DOUBLE PRECISION DABS,DMAX1,DMIN1,NXTDSC C INTEGER INTGRS(*),LOOP DOUBLE PRECISION ENDPNT,ERREST,HFACTR,HNEW,REALS(*),T LOGICAL NOTFIX,REJECT,TRCKNG C EXTERNAL DERIV,ERRCAL,ERRHAN,INTEG,LAGVAL,NXTDSC,SNGLAR,STORE INTRINSIC DABS,DMAX1,DMIN1 C C Routine for solving the functional equation for the interval specified C C NOTFIX indicates whether the fixed stepsize strategy is being used C NOTFIX = (INTGRS(4).EQ.0) C C REJECT indicates whether the last step was rejected C REJECT = .FALSE. C C Determine if we are tracking discontinuities C TRCKNG = (INTGRS(5).NE.0) C C Initialize the point upto which Volterra terms remain unchanged C REALS(14) = REALS(9) C C Main loop for solving upto the next discontinuity if tracking C C If we are tracking discontinuities or this is the first step, C since Dormand and Prince method uses the first-same-as-last C strategy, (re-)evaluate the derivative. Also re-evaluate if C the `re-evaluate first stage' flag is set. C 1 IF (TRCKNG .OR. INTGRS(22).EQ.0 .OR. INTGRS(59).EQ.1) THEN C C Set `re-evaluate first stage' flag to NO C INTGRS(59) = 0 C C Evaluate the lag functions at Tn, calculate the derivatives C and update the number of derivative functions diagnostic C C Ensure that if Tn is a discontinuity, then T > Tn C T = REALS(9)+2D0*DMAX1(DABS(REALS(9)),1D0)*REALS(1) CALL LAGVAL(INTGRS,REALS,T) INTGRS(60) = 1 CALL DERIV(INTGRS,REALS,REALS(9),REALS(INTGRS(36))) INTGRS(60) = 0 INTGRS(25) = INTGRS(25)+1 ENDIF C C Main loop used to integrate upto `solve upto' REALS(10) unless C the fixed stepsize strategy is specified C C Ensure that the current stepsize is not too small C 2 IF (REALS(11).LE.DMAX1(1D0,10D0*DABS(REALS(9)))*REALS(1)) + CALL ERRHAN(INTGRS,'Current stepsize is too small.', + 'SOLVE.',-1) C C If tracking, calculate the position of the next discontinuity C IF (TRCKNG) THEN C C Save the next continuity record to be set by NXTDSC C INTGRS(56) = INTGRS(10)+1 C C Calculate the position of the next discontinuity and set C the `solve upto' point to this point C REALS(10) = NXTDSC(INTGRS,REALS) ENDIF C C Shorten the current stepsize if necessary to reach `solve upto' C or to maintain an in-size stepsize if not using the fixed C stepsize strategy. But first save the unshortened stepsize in C case the shortened stepsize is due to unit roundoff error C IF (NOTFIX) THEN REALS(15) = REALS(11) IF (REALS(9)+REALS(11).GT.REALS(10)) THEN REALS(11) = REALS(10)-REALS(9) ELSE IF (REALS(9)+2.03125D0*REALS(11).GT.REALS(10)) + REALS(11) = (REALS(10)-REALS(9))/2D0 ENDIF C C If no accepted steps yet, do not save the untruncated stepsize C IF (INTGRS(23).EQ.0) REALS(15) = REALS(11) ENDIF C C Reset upto which point Volterra terms remain unchanged C REALS(14) = REALS(3) C C Reset the attempted extrapolation `diagnostic' C REALS(17) = 1D0 C C Increase the number of steps diagnostic C INTGRS(22) = INTGRS(22)+1 C C Print out current stepsize running diagnostic C IF (INTGRS(16).EQ.3) + WRITE (*,*) 'Attempting [',REALS(9),',',REALS(9)+REALS(11),')' C C Evaluate the Runge-Kutta equations C CALL INTEG(INTGRS,REALS,INTGRS(13),REALS(INTGRS(34)), + REALS(INTGRS(36))) C C Check if iterative refinement is necessary C IF (INTGRS(31).EQ.1) CALL SNGLAR(INTGRS,REALS) C C Calculate an error estimate if not using the fixed stepsize strategy C IF (NOTFIX) THEN CALL ERRCAL(INTGRS,REALS,ERREST) C C Calculate the new stepsize based on the classical stepsize control C algorithm C HFACTR = DMAX1(0.15D0,DMIN1(5D0, + (ERREST/REALS(2))**(1D0/5D0)/0.85D0)) C C Ensure the maximum permitted stepsize is obeyed if set C IF (INTGRS(3).EQ.0) THEN HNEW = REALS(11)/HFACTR ELSE HNEW = DMIN1(REALS(11)/HFACTR,REALS(12)) ENDIF ENDIF C C If error tolerance satisfied or fixed stepsize strategy C selected then calculate the approximant and store it C IF (ERREST.LE.REALS(2) .OR. .NOT.NOTFIX) THEN C C Update the number of accepted steps diagnostic C INTGRS(23) = INTGRS(23)+1 C C Calculate and store the 5th order Hermite approximant C CALL STORE(INTGRS,REALS,REALS(INTGRS(34)),REALS(INTGRS(35)), + REALS(INTGRS(46)),REALS(INTGRS(36)),INTGRS(13), + REALS(INTGRS(48)),REALS(INTGRS(49)),REALS(INTGRS(50)), + REALS(INTGRS(51)),REALS(INTGRS(52)),REALS(INTGRS(53)), + REALS(INTGRS(54))) C C Update the solution at the start of the interval and, since using C the first-same-as-last strategy, update the derivative estimates C DO 3 LOOP = 0,INTGRS(13)-1 REALS(INTGRS(34)+LOOP) = REALS(INTGRS(35)+LOOP) REALS(INTGRS(36)+LOOP) = + REALS(INTGRS(36)+6*INTGRS(13)+LOOP) 3 CONTINUE C C Reset iterative refinement flag to normal operation C INTGRS(31) = 0 C C Update the current point C REALS(9) = REALS(9)+REALS(11) C C Update the accepted stepsize and extrapolation diagnostics C IF (REALS(11).GT.REALS(4)) REALS(4) = REALS(11) IF (REALS(11).LT.REALS(5)) REALS(5) = REALS(11) IF (REALS(17).GT.REALS(6)) REALS(6) = REALS(17) C C Do not allow an increase in stepsize if the last stepsize rejected C IF (REJECT) HNEW = DMIN1(REALS(11),HNEW) REJECT = .FALSE. ELSE C C If the current was rejected C REJECT = .TRUE. C C If this was not due to the stepsize being a too large an estimate for C the initial stepsize, update the number of rejected steps diagnostic C IF (INTGRS(23).GT.0) INTGRS(24) = INTGRS(24)+1 C C If the rejected stepsize occurred after using the iterative refinement C algorithm, then reset the approximant data so that it is correct C IF (INTGRS(31).EQ.3) THEN INTGRS(30) = INTGRS(30)-1 REALS(3) = REALS(9) REALS(19) = REALS(3)+DMAX1(DABS(REALS(3)),1D0)*REALS(1) ENDIF ENDIF C C Set the next stepsize and reset iterative refinement to normal operation C REALS(11) = HNEW INTGRS(31) = 0 C C Print out the diagnostics if requested C IF (INTGRS(16).EQ.2) THEN IF (INTGRS(4).EQ.0) THEN WRITE (*,9998) INTGRS(22),INTGRS(23),INTGRS(24), + INTGRS(25),INTGRS(26),INTGRS(28),INTGRS(27) ELSE WRITE (*,9999) + INTGRS(22),INTGRS(25),INTGRS(26),INTGRS(28),INTGRS(27) ENDIF ENDIF C C Repeat main integration loop until `solve upto' reached, C checking if the `re-evaluate first stage' flag is set C ENDPNT = REALS(9)+DMAX1(DABS(REALS(9)),1D0)*REALS(1) IF (NOTFIX .AND. ENDPNT.LT.REALS(10)) THEN IF (INTGRS(59).EQ.1) GOTO 1 GOTO 2 ENDIF C C Restore the unshortened in-size stepsize C REALS(11) = REALS(15) C C If `solve upto' is not the endpoint and not using fixed stepsize C loop to find next discontinuity `solve upto' point C IF (NOTFIX .AND. ENDPNT.LT.REALS(8)) GOTO 1 C 9998 FORMAT (7I10) 9999 FORMAT (5I10) C RETURN END C C C ************************************************************************* C C DOUBLE PRECISION FUNCTION RULE1(INTGRS,REALS,H,X,T) C DOUBLE PRECISION DSQRT,KERNEL,TTRANS C INTEGER INTGRS(*) DOUBLE PRECISION ADD,H,NEWX,REALS(*),T,X C EXTERNAL KERNEL,TTRANS INTRINSIC DSQRT C C 4 point 7th order Gauss-Legendre rule C NEWX = TTRANS(H,X,-DSQRT((3D0+DSQRT(24D0/5D0))/7D0)) ADD = (1D0-DSQRT(5D0/54D0))*KERNEL(INTGRS,REALS,NEWX,T,INTGRS(58)) NEWX = TTRANS(H,X,-DSQRT((3D0-DSQRT(24D0/5D0))/7D0)) ADD = ADD+(1D0+DSQRT(5D0/54D0))* + KERNEL(INTGRS,REALS,NEWX,T,INTGRS(58)) NEWX = TTRANS(H,X,DSQRT((3D0-DSQRT(24D0/5D0))/7D0)) ADD = ADD+(1D0+DSQRT(5D0/54D0))* + KERNEL(INTGRS,REALS,NEWX,T,INTGRS(58)) NEWX = TTRANS(H,X,DSQRT((3D0+DSQRT(24D0/5D0))/7D0)) ADD = ADD+(1D0-DSQRT(5D0/54D0))* + KERNEL(INTGRS,REALS,NEWX,T,INTGRS(58)) C RULE1 = H*ADD*(1D0/4D0) C C Update the kernel evaluations diagnostic C INTGRS(26) = INTGRS(26)+4 C RETURN END C C C ******************************************************************+++++++ C C DOUBLE PRECISION FUNCTION RULE2(INTGRS,REALS,H,X,T) C DOUBLE PRECISION DSQRT,KERNEL,TTRANS C INTEGER INTGRS(*) DOUBLE PRECISION ADD,H,NEWX,REALS(*),T,X C EXTERNAL KERNEL,TTRANS INTRINSIC DSQRT C C 3 point 5th order Gauss-Legendre rule C NEWX = TTRANS(H,X,-DSQRT(3D0/5D0)) ADD = 5D0*KERNEL(INTGRS,REALS,NEWX,T,INTGRS(58)) NEWX = TTRANS(H,X,0D0) ADD = ADD+8D0*KERNEL(INTGRS,REALS,NEWX,T,INTGRS(58)) NEWX = TTRANS(H,X,DSQRT(3D0/5D0)) ADD = ADD+5D0*KERNEL(INTGRS,REALS,NEWX,T,INTGRS(58)) C RULE2 = ADD*H*(1D0/18D0) C C Update the kernel evaluations diagnostic C INTGRS(26) = INTGRS(26)+3 C RETURN END C C C ******************************************************************+++++++ C C DOUBLE PRECISION FUNCTION RULE3(INTGRS,REALS,H,X,T) C DOUBLE PRECISION DSQRT,KERNEL,TTRANS C INTEGER INTGRS(*) DOUBLE PRECISION ADD,H,NEWX,REALS(*),T,X C EXTERNAL KERNEL,TTRANS INTRINSIC DSQRT C C 2 point 3rd order Gauss-Legendre rule C NEWX = TTRANS(H,X,-DSQRT(1D0/3D0)) ADD = KERNEL(INTGRS,REALS,NEWX,T,INTGRS(58)) NEWX = TTRANS(H,X,DSQRT(1D0/3D0)) ADD = ADD+KERNEL(INTGRS,REALS,NEWX,T,INTGRS(58)) C RULE3 = ADD*H*(1D0/2D0) C C Update the kernel evaluations diagnostic C INTGRS(26) = INTGRS(26)+2 C RETURN END C C C ******************************************************************+++++++ C C DOUBLE PRECISION FUNCTION SCANJUMP(INTGRS,REALS,LIMIT) C DOUBLE PRECISION DABS,DMAX1 C INTEGER INTGRS(*),POS DOUBLE PRECISION CURRNT,DISC,LIMIT,PRTURB,REALS(*) LOGICAL PASTIT C INTRINSIC DABS,DMAX1 C C Function to scan the continuity records and return the C first discontinuity after the current point C C LIMIT is the furthest point we can `solve upto', namely the endpoint C (LIMIT may be reduced if tracking state dependent discontinuities C which require extrapolation to `guess' the lag function) C C Get the current point (allowing for some rounding error) C CURRNT = REALS(9) PRTURB = CURRNT+DMAX1(DABS(CURRNT),1D0)*REALS(1) C C Get the first continuity record after the current point C POS = INTGRS(61) C C Repeat for all continuity records set after the current point C PASTIT = .TRUE. 1 IF (POS.LE.INTGRS(10)) THEN C C Get the position of the discontinuity C DISC = REALS(INTGRS(55)+POS-1) C C Update the `first discontinuity' after current point pointer C IF (PASTIT) THEN PASTIT = (DISC.LT.CURRNT) INTGRS(61) = POS ENDIF C C Check if a nearer discontinuity has been found (allowing C for rounding error in the position of the discontinuity) C IF (DISC+DMAX1(DABS(DISC),1D0)*REALS(1).LT.LIMIT .AND. + DISC.GT.PRTURB) LIMIT = DISC C C Increment the pointer to the continuity records and repeat C POS = POS+1 GOTO 1 ENDIF C C The variable LIMIT contains the next discontinuity C SCANJUMP = LIMIT C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE SECANT(INTGRS,REALS,LAGFN,DISCAT,SECND,CODE) C DOUBLE PRECISION DABS,DMAX1,DMIN1,LAG C INTEGER CODE,COUNT,INTGRS(*),LAGFN,LOOP,POS,TYPE DOUBLE PRECISION DISCAT,FFIRST,FIRST,FSECND DOUBLE PRECISION NEXT,REALS(*),SECND,TOP C EXTERNAL LAG INTRINSIC DABS,DMAX1,DMIN1 C C Routine to implement the Secant method to track the discontinuity C at DISCAT propagated by the lag function LAGFN C C Set the Secant bracket to [Tn,min(Tn+2h,endpoint)] C FIRST = DMIN1(REALS(9)+2D0*REALS(11),REALS(8)) SECND = REALS(9) C C Set POS to point to the start of the lag values in REALS C POS = INTGRS(45)-1 C C Get the type of lag function (state in\dependent) C TYPE = INTGRS(INTGRS(43)+LAGFN-1) C C Set the top of the bracket outside of which we stop iterating, the C bottom of the bracket is the current point. C TOP = FIRST C C Set the `extrapolate during tracking discontinuities' flag to `on' C INTGRS(57) = 1 C C Set an iteration counter to ensure that the iteration terminates C and calculate the lag values for FIRST and SECND C COUNT = 0 C IF (TYPE.EQ.0) THEN C C In the case of state independent lag functions, we only C need to evaluate the one lag function C FFIRST = LAG(INTGRS,REALS,LAGFN,FIRST)-DISCAT FSECND = LAG(INTGRS,REALS,LAGFN,SECND)-DISCAT ELSE C C For state dependent delays, without knowing more information C about how the lag functions connect, we must evaluate all the C lag functions upto and including LAGFN C C First do for the point FIRST, the ODE lag function first C REALS(POS) = FIRST DO 1 LOOP = 1,LAGFN REALS(POS+LOOP) = LAG(INTGRS,REALS,LOOP,FIRST) 1 CONTINUE FFIRST = REALS(POS+LAGFN)-DISCAT C C Repeat for the point SECND, the ODE lag function first C REALS(POS) = SECND DO 2 LOOP = 1,LAGFN REALS(POS+LOOP) = LAG(INTGRS,REALS,LOOP,SECND) 2 CONTINUE FSECND = REALS(POS+LAGFN)-DISCAT ENDIF C C The Secant iteration loop, ensuring that we remain in the bracket C 3 IF (DABS(FFIRST-FSECND).LE.REALS(1)) THEN IF (DABS(FIRST-SECND).LE.REALS(1)) THEN C C If the bracket has become zero then signal a convergence failure. C This corresponds to a highly sensitive lag function C NEXT = 2D0*TOP ELSE NEXT = (FIRST+SECND)/2D0 ENDIF ELSE NEXT = SECND-FSECND*(SECND-FIRST)/(FSECND-FFIRST) ENDIF FFIRST = FSECND FIRST = SECND SECND = NEXT C C Check if termination condition met: C 1) SECND < current point or SECND > current point+2*h C IF (SECND.LT.REALS(9) .OR. SECND.GT.TOP) THEN C C Set the iteration return code CODE = -1 for failure C CODE = -1 C C Set the `extrapolate during tracking discontinuities' flag to `off' C INTGRS(57) = 0 RETURN ENDIF C C Calculate the new lag function value C IF (TYPE.EQ.0) THEN C C State independent case first C FSECND = LAG(INTGRS,REALS,LAGFN,SECND)-DISCAT ELSE C C State dependent case C C Reset the `extrapolate during tracking discontinuities' flag C INTGRS(57) = 1 C REALS(POS) = SECND DO 4 LOOP = 1,LAGFN REALS(POS+LOOP) = LAG(INTGRS,REALS,LOOP,SECND) 4 CONTINUE FSECND = REALS(POS+LAGFN)-DISCAT ENDIF C C Update the counter C COUNT = COUNT+1 C C Check if termination conditions met: C 2) COUNT > 50 C 3) iterative refinement signalled - initial interval C 4) ABS(FSECND) < scaled round off C IF (COUNT.GT.50 .OR. INTGRS(31).EQ.1) THEN C C Set the iteration return code CODE = -1 for failure C CODE = -1 ELSE IF (DABS(FSECND).LT.4D0*DMAX1(DABS(DISCAT),1D0)*REALS(1)) THEN IF (INTGRS(57).EQ.1) THEN C C Set the iteration return code CODE = 1 for success C CODE = 1 ELSE C C Set the iteration return code CODE = 2 for success but C required extrapolation to find the root C CODE = 2 ENDIF ELSE C C Repeat the iteration C GOTO 3 ENDIF ENDIF C C Set the `extrapolate during tracking discontinuities' flag to `off' C INTGRS(57) = 0 C RETURN END C C C ******************************************************************+++++++ C C DOUBLE PRECISION FUNCTION SEMIAP(INTGRS,REALS,T,FROM,TO) C DOUBLE PRECISION DABS,DBLE,DMAX1,DMIN1,RULE1,RULE2,RULE3 C INTEGER INTGRS(*) DOUBLE PRECISION FROM,REALS(*),T,TO DOUBLE PRECISION AT,ERREST,ERRORD,EXPECT,H,HIGH,HNEW,TOLEFT LOGICAL REJECT C EXTERNAL RULE1,RULE2,RULE3 INTRINSIC DABS,DBLE,DMAX1,DMIN1 C C Semi-adaptive quadrature scheme C C Set the value of the integral and `requested' tolerance C SEMIAP = 0D0 TOLEFT = DMAX1(10D0*REALS(1),REALS(2)) C C If iterative refinement signalled, abandon integration C IF (INTGRS(31).EQ.1) RETURN C C Calculate the order of the error estimate C ERRORD = 1D0/(7D0-2D0*DBLE(INTGRS(64))) C C Set last-step-rejected flag C REJECT = .FALSE. C C Choose an initial stepsize for the adaptive quadrature scheme C H = DMAX1((TO-FROM)/100D0,100D0*REALS(1)) C C Initialize the current point C AT = FROM C C Start of main integration loop C C Check to see if the integration is finished C 1 IF (AT+DMAX1(DABS(AT),1D0)*REALS(1).GE.TO) RETURN C C Calculate the high order integral estimate, choice of integration C rule depends on `kernel smoothness' flag C IF (INTGRS(64).EQ.0) THEN HIGH = RULE1(INTGRS,REALS,H,AT,T) ELSE HIGH = RULE2(INTGRS,REALS,H,AT,T) ENDIF C C If iterative refinement signalled, abandon integration C IF (INTGRS(31).EQ.1) RETURN C C Calculate a lower order integral error estimate, choice of integration C rule depends on `kernel smoothness' flag. If the size of the high order C integral estimate is less than 1 use an absolute error estimate, C otherwise use a relative error estimate. Set order of error estimate C IF (INTGRS(64).EQ.0) THEN ERREST = DABS(HIGH-RULE2(INTGRS,REALS,H,AT,T))/ + DMAX1(1D0,DABS(HIGH)) ELSE ERREST = DABS(HIGH-RULE3(INTGRS,REALS,H,AT,T))/ + DMAX1(1D0,DABS(HIGH)) ENDIF C C If iterative refinement signalled, abandon integration C IF (INTGRS(31).EQ.1) RETURN C C Calculate the expected error estimate based on the amount of C error still unused and the size of the interval being solved C EXPECT = DMAX1(4D0*REALS(1),TOLEFT*H/(TO-AT)) C C Calculate the next stepsize using a Runge-Kutta stepsize algorithm C HNEW = DMAX1(DABS(AT)*REALS(1),H/ + DMAX1(0.2D0,DMIN1(4D0,(ERREST/EXPECT)**ERRORD/0.65D0))) C C If the expected error estimate is satisfied OR the current C stepsize is approaching its level of significance, accept C the current stepsize C IF (ERREST.LE.EXPECT .OR. H.LE.DABS(AT)*REALS(1)) THEN C C Update the integral, the current point and the amount of C tolerance still remaining C AT = AT+H SEMIAP = SEMIAP+HIGH TOLEFT = TOLEFT-ERREST C C Do not allow a stepsize increase if the last step was rejected C IF (REJECT) HNEW = DMIN1(H,HNEW) REJECT = .FALSE. ELSE REJECT = .TRUE. ENDIF C C Set the next stepsize, ensuring that it does not go too far C H = DMIN1(HNEW,TO-AT) C C Repeat integration loop until whole interval is covered C GOTO 1 C END C C C ******************************************************************+++++++ C C SUBROUTINE SNGLAR(INTGRS,REALS) C INTEGER INTGRS(*) DOUBLE PRECISION REALS(*) C EXTERNAL INTEG,STORE C C Report that iterative refinement is being carried out C IF (INTGRS(15).EQ.1) + WRITE (*,*) 'Iterative refining the interval [', + REALS(9),',',REALS(9)+REALS(11),')' C C Set the iterative refinement flag so that any solution C value from the current interval is Y(Tn). C INTGRS(31) = 2 C C Solve the Runge-Kutta equations C C Number of times INTEG has been called whilst in iterative mode C INTGRS(32) = 1 C CALL INTEG(INTGRS,REALS,INTGRS(13),REALS(INTGRS(34)), + REALS(INTGRS(36))) C C Calculate the low order approximant C CALL STORE(INTGRS,REALS,REALS(INTGRS(34)),REALS(INTGRS(35)), + REALS(INTGRS(46)),REALS(INTGRS(36)),INTGRS(13), + REALS(INTGRS(48)),REALS(INTGRS(49)),REALS(INTGRS(50)), + REALS(INTGRS(51)),REALS(INTGRS(52)),REALS(INTGRS(53)), + REALS(INTGRS(54))) C C Adjust iterative refinement flag so as to use low order approximant C INTGRS(31) = 3 C C Start of iterative refinement loop C C Increase number of times INTEG has been called counter C 1 INTGRS(32) = INTGRS(32)+1 CALL INTEG(INTGRS,REALS,INTGRS(13),REALS(INTGRS(34)), + REALS(INTGRS(36))) C C If we are going to do another iteration then recalculate the C approximant, otherwise approximant is recalculated only if C the current stepsize is accepted by RETARD. Number of C iterations is determined by the order of the method. C IF (INTGRS(32).LT.5) THEN CALL STORE(INTGRS,REALS,REALS(INTGRS(34)),REALS(INTGRS(35)), + REALS(INTGRS(46)),REALS(INTGRS(36)),INTGRS(13), + REALS(INTGRS(48)),REALS(INTGRS(49)),REALS(INTGRS(50)), + REALS(INTGRS(51)),REALS(INTGRS(52)),REALS(INTGRS(53)), + REALS(INTGRS(54))) GOTO 1 ENDIF C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE SOFT(INTGRS,REALS,FROM,IN,LINK) C INTEGER FROM,IN,INTGRS(*),LINK,POS DOUBLE PRECISION REALS(*) C EXTERNAL ERRHAN,MAPSET C C Routine to specify a SOFT delay discontinuity network link C C FROM - Component of solution in which discontinuity already exists C IN - Component of solution in which discontinuity is created C LINK - Lag term which links IN and FROM solution components C C Ensure that FROM does not exceed the dimension of the problem C IF (FROM.LT.1 .OR. FROM.GT.INTGRS(13)) + CALL ERRHAN(INTGRS, + 'Discontinuity occurs in a solution that cannot exist.', + 'SOFT.',-1) C C Ensure that IN does not exceed the dimension of the problem C IF (IN.LT.1 .OR. IN.GT.INTGRS(13)) + CALL ERRHAN(INTGRS, + 'Discontinuity mapped to a solution that cannot exist.', + 'SOFT.',-1) C C Ensure that LINK corresponds to a lag function which can exist C IF (LINK.LT.0 .OR. LINK.GT.INTGRS(14)) + CALL ERRHAN(INTGRS, + 'Linked lag function cannot exist.','SOFT.',-1) C C If this is the first call to either CONTY, (N)HARD or (N)SOFT set C up the pointers to start of arrays in INTGRS. After which the number C of Volterra terms and storage allocated to the continuity records C and the HARD and SOFT discontinuity networks cannot be altered. C C INTGRS(8)=0 & INTGRS(9)=0 & INTGRS(10)=0 C IF (INTGRS(8)+INTGRS(9)+INTGRS(10).EQ.0) + CALL MAPSET(INTGRS,'SOFT.') C C Check that there is still space left for mapping C IF (INTGRS(9).EQ.INTGRS(20)) THEN CALL ERRHAN(INTGRS, + 'Insufficient space for SOFT discontinuity network.', + 'SOFT.',2) IF (INTGRS(15).GE.2) + WRITE (*,*) 'SOFT discontinuity network mapping ignored' RETURN ENDIF C C Save the SOFT discontinuity mapping: IN, FROM, LINK C POS = INTGRS(41)+3*INTGRS(9) INTGRS(POS) = FROM INTGRS(POS+1) = IN INTGRS(POS+2) = LINK C C Update the number of SOFT discontinuity mappings set C INTGRS(9) = INTGRS(9)+1 C RETURN END C C C ******************************************************************+++++++ C C DOUBLE PRECISION FUNCTION SOLNYI(INTGRS,REALS,COMPNT,AT) C DOUBLE PRECISION DELAY C INTEGER COMPNT,INTGRS(*),OLDPNT DOUBLE PRECISION AT,OLDLAG,REALS(*),TFIRST,TLAST C EXTERNAL DELAY,ERRHAN C C Function to calculate the COMPNT'th component solution value C (still held in the cyclic storage) at the point AT. C C Ensure that the component specified is valid C IF (COMPNT.LT.1 .OR. COMPNT.GT.INTGRS(13)) + CALL ERRHAN(INTGRS, + 'Invalid solution component specified.', + 'SOLNYI.',-1) C C Calculate the start and end of the history queue C TFIRST = REALS(INTGRS(48)+INTGRS(29)-1) TLAST = REALS(3)+DMAX1(DABS(REALS(3)),1D0)*REALS(1) C C Calculate the solution value using the first lag function which must C be preserved (so that we can call SOLNYI in DERIV without problems) C and preserve the original pointer in the history queue (for efficiency) C OLDLAG = REALS(INTGRS(45)) OLDPNT = INTGRS(INTGRS(42)) INTGRS(INTGRS(42)) = INTGRS(63) C IF (AT.LT.REALS(7) .OR. (AT.GE.TFIRST .AND. AT.LE.TLAST)) THEN REALS(INTGRS(45)) = AT SOLNYI = DELAY(INTGRS,REALS,1,COMPNT) ELSE PRINT *,AT,' Solution is not in history queue' SOLNYI = 0D0 ENDIF C C Restore the ODE lag function variables C INTGRS(63) = INTGRS(INTGRS(42)) INTGRS(INTGRS(42)) = OLDPNT REALS(INTGRS(45)) = OLDLAG C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE SOLVE(INTGRS,REALS,Y) C DOUBLE PRECISION DABS,DBLE,DMAX1 C INTEGER INTGRS(*),LOOP DOUBLE PRECISION ENDPNT,REALS(*),Y(*) C EXTERNAL DETERM,ERRHAN,MAPSET,RETARD INTRINSIC DABS,DBLE,DMAX1 C C Routine for solving the interval of integration in one call C of RETARD if no fixed stepsize is selected. If a fixed C stepsize is selected, then multiple calls to RETARD are made C from SOLVE, the endpoint of each interval of integration being C the next point of the fixed step. C C Ensure that an initial stepsize has been specified or that C the fixed stepsize strategy is being used C IF (INTGRS(2).EQ.0 .AND. INTGRS(4).NE.1) + CALL ERRHAN(INTGRS,'No initial stepsize specified.', + 'SOLVE.',-1) C C If tracking discontinuities, ensure at least one SOFT discontinuity C network and one continuity record have been specified C IF (INTGRS(5).GT.0 .AND. INTGRS(9).EQ.0) + CALL ERRHAN(INTGRS, + 'No discontinuity networks specified for tracking.', + 'SOLVE.',-1) IF (INTGRS(5).GT.0 .AND. INTGRS(10).EQ.0) + CALL ERRHAN(INTGRS, + 'No continuity records specified for tracking.', + 'SOLVE.',-1) C C Check that either an error-per-step control or the fixed C stepsize strategy has been selected C IF (INTGRS(4).EQ.0 .AND. INTGRS(6).EQ.0) + CALL ERRHAN(INTGRS, + 'An error control or fixed stepsize must be specified.', + 'SOLVE.',-1) C C Initialize any variables that were not set before because C of a dependence on the nature of the problem C C Initialize pointers to types of lag, etc if not already C done so by calling one of CONTY, HARD, NHARD, NSOFT, SOFT C IF (INTGRS(43).EQ.0) CALL MAPSET(INTGRS,'SOLVE.') C C Initialize the storage variable `last point stored' and a perturbed C version, used in improving the performance of RECALL C REALS(3) = REALS(7) REALS(19) = REALS(3)+DMAX1(DABS(REALS(3)),1D0)*REALS(1) C C Initialize diagnostics HMAXM, HMINM, EXTRAP C REALS(4) = 0D0 REALS(5) = 2D0*DMAX1(REALS(11),REALS(12)) REALS(6) = 1D0 C C Initialize `solve upto' to the endpoint C REALS(10) = REALS(8) C C Calculate the start of the lag values in REALS, allowing for C the ODE zero lag function at the start C INTGRS(45) = INTGRS(55)+INTGRS(21)+1 C C Calculate the start of solution at Tn+H/2 in REALS C INTGRS(46) = INTGRS(45)+INTGRS(14)+INTGRS(44) C C Calculate the start of the continuous error polynomial in REALS C INTGRS(47) = INTGRS(46)+INTGRS(13) C C Calculate the start of the solution Y in REALS C IF (INTGRS(6).EQ.1) THEN C C If asymptotic error estimate selected, no storage is required C for the continuous error polynomial but we reserve DIMEN C elements for temporary storage (of Y, when the two extra stages C of the Runge-Kutta method are calculated which are necessary C to calculate the fifth order Hermite interpolant). C INTGRS(34) = INTGRS(47)+INTGRS(13) ELSE C C If a continuous error estimate is selected, storage requirement C is dimension of problem * (degree of polynomial (5) + 1) C INTGRS(34) = INTGRS(47)+INTGRS(13)*6 ENDIF C C Calculate the start of the intermediate solution Y1 in REALS C INTGRS(35) = INTGRS(34)+INTGRS(13) C C Calculate the start of the projected derivatives KI in REALS C INTGRS(36) = INTGRS(35)+INTGRS(13) C C Calculate the start of the Volterra shortcut storage in REALS C The total number of stages in the Runge-Kutta method is 9 C INTGRS(37) = INTGRS(36)+INTGRS(13)*9 C C Calculate the start of the cyclic storage structure in REALS C The total number of stages in the Runge-Kutta method is 9 C INTGRS(38) = INTGRS(37)+INTGRS(18)*9 C C Calculate how long the cyclic storage structure in REALS can be. C The degree of the approximant is 5, so total is LENGTH*(DIMEN*6+1) C INTGRS(33) = (INTGRS(11)-INTGRS(38)+1)/(6*INTGRS(13)+1) C C Ensure that REALS allows enough space for cyclic storage array C IF (INTGRS(33).LT.10) + CALL ERRHAN(INTGRS,'Size of REALS is too small.','SOLVE.',-1) IF (INTGRS(15).GE.2) + WRITE (*,*) 'Length of cyclic storage is ',INTGRS(33) C C Calculate the offset in REALS for the 7 arrays used for storing C the approximation data C INTGRS(48) = INTGRS(38) INTGRS(49) = INTGRS(48)+INTGRS(33) DO 1 LOOP = 3,7 INTGRS(47+LOOP) = INTGRS(46+LOOP)+INTGRS(33)*INTGRS(13) 1 CONTINUE C C Transfer the initial solution to Y and Y1 in REALS C DO 2 LOOP = 1,INTGRS(13) REALS(INTGRS(34)+LOOP-1) = Y(LOOP) REALS(INTGRS(35)+LOOP-1) = Y(LOOP) 2 CONTINUE C C Determine which lags are state dependent/independent C CALL DETERM(INTGRS,REALS) C C Set the iterative refinement flag to `off' C INTGRS(31) = 0 C C Set the `re-evaluate first stage' flag to NO C INTGRS(59) = 0 C C Pointer into the continuity record for the first discontinuity C after the current point (for determining the `next' discontinuity) C INTGRS(61) = 1 C C Pointer into the continuity record to the oldest discontinuity C that can be propagated with the bound on the size of the lags C INTGRS(62) = 1 C C Pointer into the continuity record used by SOLNYI C INTGRS(63) = 1 C C Now solve the actual integration problem C C If no fixed stepsize is selected, solve the problem in one call C IF (INTGRS(4).NE.1) THEN C C If diagnostics requested print out headline C IF (INTGRS(16).EQ.2) WRITE (*,*) + ' STEPS ACCEPT REJECT FCN', + ' KERNEL YLAG EXTRAP' C CALL RETARD(INTGRS,REALS) ELSE C C Solve using a fixed stepsize C C If diagnostics requested print out headline C IF (INTGRS(16).EQ.2) WRITE (*,*) + ' STEPS FCN KERNEL YLAG EXTRAP' C C Repeat RETARD until whole interval is solved C C Calculate the current stepsize using a defect approach C based on where we should actually solve upto C 3 REALS(11) = DBLE(INTGRS(22)+1)*REALS(13)+REALS(7)-REALS(9) C C Solve the next `fixed' interval C CALL RETARD(INTGRS,REALS) C C Repeat until the whole interval has been solved C ENDPNT = REALS(9)+DMAX1(DABS(REALS(9)),1D0)*REALS(1) IF (ENDPNT.LT.REALS(8)) GOTO 3 ENDIF C C If diagnostics requested print them all out C IF (INTGRS(16).GT.0) THEN WRITE (*,*) WRITE (*,*) 'Total number of steps = ',INTGRS(22) WRITE (*,*) 'Number of accepted steps = ',INTGRS(23) WRITE (*,*) 'Number of rejected steps = ',INTGRS(24) WRITE (*,*) 'Number of derivative calls = ',INTGRS(25) WRITE (*,*) 'Number of approximant calls = ',INTGRS(28) WRITE (*,*) 'Number of kernel calls = ',INTGRS(26) WRITE (*,*) 'Number of extrapolations = ',INTGRS(27) WRITE (*,*) 'Number of continuity records = ',INTGRS(10) WRITE (*,*) WRITE (*,*) 'Maximum accepted stepsize = ',REALS(4) WRITE (*,*) 'Minimum accepted stepsize = ',REALS(5) WRITE (*,*) 'Maximum extrapolation = ',REALS(6)-1D0,'H' WRITE (*,*) WRITE (*,*) ENDIF C C Finally transfer the solution back into Y C DO 4 LOOP = 1,INTGRS(13) Y(LOOP) = REALS(INTGRS(34)+LOOP-1) 4 CONTINUE C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE SPACE(INTGRS,REALS,DIST) C INTEGER INTGRS(*) DOUBLE PRECISION DIST,REALS(*) C EXTERNAL ERRHAN C C Routine to set the minimum distance between tracked discontinuities C before they are considered as being distinct C C Ensure that the distance is not negative C IF (DIST.LT.0D0) THEN CALL ERRHAN(INTGRS,'Distance specified is negative.', + 'SPACE.',2) IF (INTGRS(15).GE.2) WRITE (*,*) 'Current distance retained' ELSE REALS(16) = DIST ENDIF C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE STATS(INTGRS,REALS,TYPE) C INTEGER INTGRS(*),TYPE DOUBLE PRECISION REALS(*) C EXTERNAL ERRHAN C C Routine for specifying the type of diagnostics to be generated C C TYPE = 0 - no diagnostics C = 1 - diagnostics to be given at end of integration C = 2 - limited diagnostics during run + all diagnostics at end C = 3 - current point and stepsize during run + all diagnostics C IF (TYPE.LT.0 .OR. TYPE.GT.3) THEN CALL ERRHAN(INTGRS, + 'Invalid level of diagnostics requested.','STATS.',2) IF (INTGRS(15).GE.2) + WRITE (*,*) 'Current diagnostic level retained' ELSE INTGRS(16) = TYPE ENDIF C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE STORCY(INTGRS,REALS,SPACE) C INTEGER INTGRS(*),SPACE DOUBLE PRECISION REALS(*) C EXTERNAL ERRHAN C C Routine to alter the storage allocated to continuity records C C C Ensure that no discontinuity networks or records have been given C IF (INTGRS(8).NE.0 .OR. INTGRS(9).NE.0 .OR. INTGRS(10).NE.0) THEN CALL ERRHAN(INTGRS, + 'Continuity record storage cannot be changed.', + 'STORCY.',2) IF (INTGRS(15).GE.2) + WRITE (*,*) 'Current continuity record storage retained' RETURN ENDIF C IF (SPACE.LT.10) THEN CALL ERRHAN(INTGRS, + 'Continuity records need at least 10 storage spaces.', + 'STORCY.',2) IF (INTGRS(15).GE.2) + WRITE (*,*) 'Current continuity record storage retained' ELSE INTGRS(21) = SPACE ENDIF C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE STORE(INTGRS,REALS,Y,Y1,Y12,KI,DIMEN,TSTOR, + YSTOR,DEG1,DEG2,DEG3,DEG4,DEG5) C INTEGER MAX0,MOD DOUBLE PRECISION DABS,DMAX1 C INTEGER ADDRSS,DIMEN,I,INTGRS(*) DOUBLE PRECISION KI(DIMEN,9),REALS(*),Y(*),Y1(*),Y12(*) DOUBLE PRECISION TSTOR(*),YSTOR(DIMEN,*) DOUBLE PRECISION DEG1(DIMEN,*),DEG2(DIMEN,*),DEG3(DIMEN,*) DOUBLE PRECISION DEG4(DIMEN,*),DEG5(DIMEN,*) C INTRINSIC DABS,DMAX1,MAX0,MOD C C Routine to calculate and store fifth order Hermite interpolant C C Update the last point stored, the perturbed version and pointer to last C stored approximant if not in iterative refinement mode C IF (INTGRS(31).NE.3) THEN REALS(3) = REALS(9)+REALS(11) REALS(19) = REALS(3)+DMAX1(DABS(REALS(3)),1D0)*REALS(1) INTGRS(30) = INTGRS(30)+1 ENDIF C C Re-calculate the position of the oldest approximant still stored C in case part of it was just overwritten C INTGRS(29) = MOD(MAX0(0,INTGRS(30)-INTGRS(33)),INTGRS(33))+1 C C Calculate where to store the approximant, since INTGRS(30) is C the total number of approximants stored C ADDRSS = MOD(INTGRS(30)-1,INTGRS(33))+1 C C Store the approximant data C TSTOR(ADDRSS) = REALS(9) C DO 1 I = 1,DIMEN YSTOR(I,ADDRSS) = Y(I) DEG1(I,ADDRSS) = KI(I,1)*REALS(11) DEG2(I,ADDRSS) = 7D0*Y1(I)-23D0*Y(I)+16D0*Y12(I)- + (6D0*KI(I,1)+KI(I,7)+8D0*KI(I,9))*REALS(11) DEG3(I,ADDRSS) = 66D0*Y(I)-34D0*Y1(I)-32D0*Y12(I)+ + (13D0*KI(I,1)+5D0*KI(I,7)+32D0*KI(I,9))*REALS(11) DEG4(I,ADDRSS) = 52D0*Y1(I)-68D0*Y(I)+16D0*Y12(I)- + (12D0*KI(I,1)+8D0*KI(I,7)+40D0*KI(I,9))*REALS(11) DEG5(I,ADDRSS) = 24D0*Y(I)-24D0*Y1(I)+ + (4D0*KI(I,1)+4D0*KI(I,7)+16D0*KI(I,9))*REALS(11) 1 CONTINUE C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE STORHD(INTGRS,REALS,SPACE) C INTEGER INTGRS(*),SPACE DOUBLE PRECISION REALS(*) C EXTERNAL ERRHAN C C Routine to alter the storage allocated to HARD network C C C Ensure that no discontinuity networks or records have been given C IF (INTGRS(8).NE.0 .OR. INTGRS(9).NE.0 .OR. INTGRS(10).NE.0) THEN CALL ERRHAN(INTGRS, + 'HARD network storage cannot be changed.','STORHD.',2) IF (INTGRS(15).GE.2) + WRITE (*,*) 'Current HARD network storage retained' RETURN ENDIF C IF (SPACE.LT.1) THEN CALL ERRHAN(INTGRS, + 'HARD network needs at least 1 storage space.', + 'STORHD.',2) IF (INTGRS(15).GE.2) + WRITE (*,*) 'Current HARD network storage retained' ELSE INTGRS(19) = SPACE ENDIF C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE STORST(INTGRS,REALS,SPACE) C INTEGER INTGRS(*),SPACE DOUBLE PRECISION REALS(*) C EXTERNAL ERRHAN C C Routine to alter the storage allocated to SOFT network C C Ensure that no discontinuity networks or records have been given C IF (INTGRS(8).NE.0 .OR. INTGRS(9).NE.0 .OR. INTGRS(10).NE.0) THEN CALL ERRHAN(INTGRS, + 'SOFT network storage cannot be changed.','STORST.',2) IF (INTGRS(15).GE.2) + WRITE (*,*) 'Current SOFT network storage retained' RETURN ENDIF C IF (SPACE.LT.1) THEN CALL ERRHAN(INTGRS, + 'SOFT network needs at least 1 storage space.', + 'STORST.',2) IF (INTGRS(15).GE.2) + WRITE (*,*) 'Current SOFT network storage retained' ELSE INTGRS(20) = SPACE ENDIF C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE TRACK(INTGRS,REALS,LEVEL) C INTEGER INTGRS(*),LEVEL DOUBLE PRECISION REALS(*) C EXTERNAL ERRHAN C C Routine for specifying that discontinuities are to be tracked C LEVEL = 0 - do not track any discontinuities (default) C LEVEL = 1 - only track state independent discontinuities C LEVEL = 2 - track state dependent discontinuities as well C C Check if a valid level of tracking has been given C IF (LEVEL.LT.0 .OR. LEVEL.GT.2) THEN CALL ERRHAN(INTGRS, + 'Invalid level of discontinuity tracking.', + 'TRACK.',2) IF (INTGRS(15).GE.2) WRITE (*,*) + 'Current level of discontinuity tracking retained' RETURN ENDIF C C Check if no tracking has been specified first, as this option C does not conflict with the fixed stepsize strategy C IF (LEVEL.EQ.0) THEN INTGRS(5) = 0 IF (INTGRS(15).GE.2) + WRITE (*,*) 'No tracking of discontinuities specified' RETURN ENDIF C C If the fixed stepsize strategy has been selected, ignore this call C IF (INTGRS(4).EQ.1) THEN CALL ERRHAN(INTGRS, + 'Fixed stepsize strategy already selected.','TRACK.',2) IF (INTGRS(15).GE.2) WRITE (*,*) + 'Request to track discontinuities ignored' RETURN ENDIF C C Issue an advisory message if error reporting level high enough C IF (INTGRS(15).GE.2) THEN IF (LEVEL.EQ.1) THEN WRITE (*,*) 'Tracking state independent discontinuities' ELSE WRITE (*,*) 'Tracking all discontinuities' ENDIF ENDIF C INTGRS(5) = LEVEL C RETURN END C C C ******************************************************************+++++++ C C DOUBLE PRECISION FUNCTION TTRANS(H,AT,POINT) C DOUBLE PRECISION AT,H,POINT C TTRANS = AT+H*(POINT+1D0)*(1D0/2D0) C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE VJUMP(INTGRS,REALS) C INTEGER INTGRS(*) DOUBLE PRECISION REALS(*) C C Routine to specify that the smoothness of a kernel function is poor C C Issue an appropriate message if required C IF (INTGRS(15).GE.2) WRITE (*,*) 'Non-smooth integrand specified' C INTGRS(64) = 1 C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE VLAG(INTGRS,REALS,LAGNO,LAG) C INTEGER INTGRS(*),LAGNO DOUBLE PRECISION LAG,REALS(*) C EXTERNAL ERRHAN C C Routine to evaluate a single Volterra kernel lag function. C C The lag number must NOT be one already specified for use C as a DELAY or NEUTRL lag function. This enables kernel C functions with state dependent delays to be evaluated (as C otherwise the call to evaluate the lag at the point S passed C to the kernel function would require a recursive call of C LAG or the interdependence of state dependent and independent C lags to be specified). This also avoids recalculating all the C lag functions for the normal DELAY and NEUTRL terms at S, C which would then need recalculating at the point T before C control is passed back to DERIV. Further not specifying the C kernel lag functions in LAG allows more exotic lag functions C to be attempted such as S-T. C C Ensure that the lag function specified is not a `normal' lag C function as specified in the subroutine LAG C IF (LAGNO.LE.INTGRS(14)) + CALL ERRHAN(INTGRS,'Kernel lag function is too low.', + 'VLAG.',-1) C C Ensure that the necessary storage is available (set in VOLTER) C IF (LAGNO.GT.INTGRS(14)+INTGRS(44)) + CALL ERRHAN(INTGRS,'Kernel lag function cannot exist.', + 'VLAG.',-1) C C The actual value of the lag function is calculated as C the argument LAG is passed, so we just have to store it C REALS(INTGRS(45)+LAGNO-1) = LAG C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE VOLTER(INTGRS,REALS,NUMBER,DELAYS) C INTEGER DELAYS,INTGRS(*),NUMBER DOUBLE PRECISION REALS(*) C EXTERNAL ERRHAN C C Optional routine to change the default number of accelerated Volterra C terms and the number of kernel lag functions. C C Ensure that no discontinuity networks or records have been given C IF (INTGRS(8).NE.0 .OR. INTGRS(9).NE.0 .OR. INTGRS(10).NE.0) THEN CALL ERRHAN(INTGRS, + 'Number of Volterra terms cannot be changed.','VOLTER.',2) IF (INTGRS(15).GE.2) + WRITE (*,*) 'Current number of Volterra terms retained' RETURN ENDIF C C Ensure that the number specified is positive C IF (NUMBER.LT.0) THEN CALL ERRHAN(INTGRS, + 'Invalid number of Volterra terms specified.','VOLTER.',2) IF (INTGRS(15).GE.2) + WRITE (*,*) 'Current number of Volterra terms retained' ELSE INTGRS(18) = NUMBER ENDIF C C Ensure that the number of kernel lag functions allowed is C non-negative C IF (DELAYS.LT.0) THEN CALL ERRHAN(INTGRS, + 'Invalid number of kernel delay terms specified.', + 'VOLTER.',2) IF (INTGRS(15).GE.2) + WRITE (*,*) 'Current number of kernel delay terms retained' ELSE INTGRS(44) = DELAYS ENDIF C RETURN END C C C ******************************************************************+++++++ C C DOUBLE PRECISION FUNCTION + VOLTRA(INTGRS,REALS,T,LAGFN1,LAGFN2,KERNNO) C DOUBLE PRECISION DABS,DMAX1,DMIN1,SEMIAP C INTEGER INTGRS(*),KERNNO,LAGFN1,LAGFN2,NFCN,VOLTNO DOUBLE PRECISION AREA,LAG1,LAG2,REALS(*),T,TEMP LOGICAL STATED C EXTERNAL ERRHAN,SEMIAP INTRINSIC DABS,DMAX1,DMIN1 C C Routine to evaluate Volterra terms from LAG1 to LAG2 using the C KERNNO'th function in KERNEL C SAVE C C Store the kernel number in INTGRS C INTGRS(58) = KERNNO C C Initialize the integral C AREA = 0D0 C C If this is the first time that VOLTRA is being called (by looking C at the number of derivative evaluations), check that both lag C functions can exist and initialize `accelerating' variables C IF (INTGRS(25).EQ.0) THEN C IF (LAGFN1.LT.0 .OR. LAGFN1.GT.INTGRS(14)) + CALL ERRHAN(INTGRS,'First lag function cannot exist.', + 'VOLTRA.',-1) C IF (LAGFN2.LT.0 .OR. LAGFN2.GT.INTGRS(14)) + CALL ERRHAN(INTGRS,'Second lag function cannot exist.', + 'VOLTRA.',-1) C C Check that the lag functions are not identical C IF (LAGFN1.EQ.LAGFN2) + CALL ERRHAN(INTGRS,'Identical lag functions specified.', + 'VOLTRA.',-1) C C Initialize `accelerating' variables - these variables are why C the SAVE statement is required C NFCN = INTGRS(25) VOLTNO = 0 ENDIF C C If iterative refinement signalled, exit immediately C IF (INTGRS(31).EQ.1) GOTO 1 C C If in the iterative refinement process, save the parts of C the integral which remain constant in the iteration C IF (INTGRS(31).GT.1) THEN C C It is possible to calculate which Volterra term is being C evaluated by how many calls to VOLTRA there have been since C the number of derivative evaluations was last increased. C IF (NFCN.EQ.INTGRS(25)) THEN C C No change in the number of derivative evaluations so this is C another Volterra term C VOLTNO = VOLTNO + 1 ELSE C C The number of derivative evaluations has changed so this must C be the first Volterra term C NFCN = INTGRS(25) VOLTNO = 1 ENDIF ENDIF C C Get the start and endpoint for the integration C LAG1 = REALS(INTGRS(45)+LAGFN1-1) LAG2 = REALS(INTGRS(45)+LAGFN2-1) C C Ensure the correct direction of integration C IF (LAG1.GT.LAG2) THEN CALL ERRHAN(INTGRS, + 'Incorrect direction of Volterra integration.', + 'VOLTRA.',2) IF (INTGRS(15).GE.2) + WRITE (*,*) 'Limits of integration swapped' TEMP = LAG1 LAG1 = LAG2 LAG2 = TEMP ENDIF C C Ensure that the endpoint is on the correct side of any C discontinuity that may occur at the endpoint C LAG1 = LAG1+2D0*DMAX1(DABS(LAG1),1D0)*REALS(1) LAG2 = LAG2-2D0*DMAX1(DABS(LAG2),1D0)*REALS(1) C C Check to see if the interval of integration is zero C to within unit-roundoff C IF (LAG1.GT.LAG2) GOTO 1 C C Determine if either of the limits of integration corresponds to a C state-dependent lag function, because STATE-DEPENDENT LIMITS OF C INTEGRATION CANNOT BE ACCELERATED UNDER ANY CIRCUMSTANCES C STATED = ((INTGRS(INTGRS(43)+LAGFN1-1).NE.0) .OR. + (INTGRS(INTGRS(43)+LAGFN2-1).NE.0)) C C Check what mode of iterative refinement we are in C IF (INTGRS(31).EQ.0 .OR. STATED) THEN C C Normal mode or a state-dependent limit of integration, so do C the integration by brute force. C C Calculate the area before the initial point C IF (LAG1.LT.REALS(7)) AREA = + SEMIAP(INTGRS,REALS,T,LAG1,DMIN1(REALS(7),LAG2)) C C Calculate the area after the initial point C IF (LAG2.GT.REALS(7)) AREA = AREA+ + SEMIAP(INTGRS,REALS,T,DMAX1(REALS(7),LAG1),LAG2) GOTO 1 ENDIF C C If iterative is signalled but not yet started exit immediately, C unless we are tracking state dependent discontinuities in which C case report an error C IF (INTGRS(31).EQ.1) THEN IF (INTGRS(57).EQ.0) GOTO 1 CALL ERRHAN(INTGRS, + 'Cannot track state dependent Volterra lag functions.', + 'SOLVE.',-1) ENDIF C C If the first iteration of iterative refinement is being carried C out, split the integral up into two parts: one part which C remains constant in the iteration and the other which varies. C IT IS IMPORTANT TO REALIZE THAT BECAUSE THE FIRST ABSCISSAE IS C ZERO, ITS DERIVATIVE VALUE REMAINS CONSTANT IN THE ITERATION C IF (INTGRS(31).EQ.2) THEN C C Calculate the area before the initial point C IF (LAG1.LT.REALS(7)) AREA = + SEMIAP(INTGRS,REALS,T,LAG1,DMIN1(REALS(7),LAG2)) C C Calculate the area after the initial point upto the last C point for which an approximant is already known C IF (LAG1.LE.REALS(14)) AREA = AREA+ + SEMIAP(INTGRS,REALS,T,DMAX1(REALS(7),LAG1), + DMIN1(REALS(14),LAG2)) C C Save the unchanging integral, to speed up the iteration if the C necessary storage is available C IF (VOLTNO.LE.INTGRS(18)) + REALS(INTGRS(37)+VOLTNO*9+INTGRS(60)-10) = AREA C C Calculate the component of the integral for which initially C no approximant is available C IF (LAG1.GT.REALS(14)) THEN AREA = SEMIAP(INTGRS,REALS,T,LAG1,LAG2) ELSE IF (LAG2.GT.REALS(14)) AREA = AREA+ + SEMIAP(INTGRS,REALS,T,REALS(14),LAG2) ENDIF GOTO 1 ENDIF C C Calculate the integral for which the approximant for the current C interval is being refined C IF (VOLTNO.GT.INTGRS(18)) THEN C C This Volterra term needs total re-evaluation because not C enough Volterra terms were specified for the problem. C C Calculate the area before the initial point C IF (LAG1.LT.REALS(7)) AREA = + SEMIAP(INTGRS,REALS,T,LAG1,DMIN1(REALS(7),LAG2)) C C Calculate the area after the initial point C IF (LAG2.GT.REALS(7)) AREA = AREA+ + SEMIAP(INTGRS,REALS,T,DMAX1(REALS(7),LAG1),LAG2) ELSE C C This is an accelerated Volterra term C C The integral is unchanged during the iteration C IF (LAG2.LE.REALS(14)) THEN AREA = REALS(INTGRS(37)+VOLTNO*9+INTGRS(60)-10) ELSE C C The lower limit of integration is not `vanishing' C IF (LAG1.LE.REALS(14)) THEN AREA = REALS(INTGRS(37)+VOLTNO*9+INTGRS(60)-10)+ + SEMIAP(INTGRS,REALS,T,REALS(14),LAG2) ELSE C C The lower limit of integration is also `vanishing' C AREA = SEMIAP(INTGRS,REALS,T,LAG1,LAG2) ENDIF ENDIF ENDIF C C Resetting INTGRS(58) allows us to determine what type of lag C function we are evaluating C 1 INTGRS(58) = 0 C VOLTRA = AREA C RETURN END C C C ******************************************************************+++++++ C C DOUBLE PRECISION FUNCTION YSTOP(INTGRS,REALS,COMPNT,START,VALUE) C DOUBLE PRECISION DABS,DELAY,NEUTRL C INTEGER COMPNT,COUNT,INTGRS(*) DOUBLE PRECISION DT,FT,REALS(*),START,T,TFIRST,TLAST,VALUE C EXTERNAL DELAY,ERRHAN,NEUTRL INTRINSIC DABS C C Function to find the point `t' at which y(t) = VALUE C using START as the starting iterate C C Ensure that the component specified is valid C IF (COMPNT.LT.1 .OR. COMPNT.GT.INTGRS(13)) + CALL ERRHAN(INTGRS, + 'Invalid solution component specified.', + 'YSTOP.',-1) C C Initialize the necessary variables C COUNT = 0 T = START TFIRST = REALS(INTGRS(48)+INTGRS(29)-1) TLAST = REALS(3)+DMAX1(DABS(REALS(3)),1D0)*REALS(1) C C Check if the root-finder has gone out of bounds C 1 IF ((T.GE.REALS(7) .AND. T.LT.TFIRST) .OR. T.GT.TLAST) THEN CALL ERRHAN(INTGRS, + 'The root-finder is out of bounds.','YSTOP.',2) YSTOP = 0D0 RETURN ENDIF C C Set the point at which the solution/derivative is evaluated C REALS(INTGRS(45)) = T FT = DELAY(INTGRS,REALS,1,COMPNT)-VALUE IF (DABS(FT).LT.5D0*REALS(1)) THEN YSTOP = T RETURN ENDIF IF (COUNT.GT.50) THEN CALL ERRHAN(INTGRS, + 'Maximum 50 iterations exceeded.','YSTOP.',2) YSTOP = 0D0 RETURN ENDIF DT = NEUTRL(INTGRS,REALS,1,COMPNT) IF (DABS(DT).LT.50D0*REALS(1)) THEN T = T+1000D0*REALS(1) ELSE T = T-FT/DT ENDIF C C Increment iteration counter and repeat iteration C COUNT = COUNT+1 GOTO 1 C END