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 SUBROUTINE BOUND(INTGRS,REALS,PARAM,LOWER,UPPER) C INTEGER INTGRS(*),PARAM DOUBLE PRECISION LOWER,REALS(*),UPPER 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(73) 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 SUBROUTINE DATUM(INTGRS,REALS,NAME) C INTEGER COUNT,INTGRS(*),LOOP,LOOP2,OFFSET,TPOINT,YPOINT CHARACTER*(*) NAME DOUBLE PRECISION REALS(*),TDATA,YDATA C EXTERNAL ERRHAN C C Routine to specify the file containing the data to be fitted and C to read it in C C Check that the data has not already been read in C IF (INTGRS(64).NE.0) THEN CALL ERRHAN(INTGRS,'Data to be fitted can only be read once.', + 'DATUM.',2) IF (INTGRS(15).GE.2) + WRITE (*,*) 'Previous data to be fitted retained.' RETURN ENDIF C C Check that the file really exists by specifying STATUS C OPEN (UNIT=44,STATUS='OLD',FILE=NAME) C C Read in the number of data points for each component of the C solution. Use space in INTGRS to store this data. NOTE: Dimension C of the problem must already have been specified. C INTGRS(12) = INTGRS(12)-INTGRS(13) C C Use of variable OFFSET overcomes a bug in the DEC Alpha (OSF1) C Fortran compiler when optimization is used C COUNT = 0 OFFSET = INTGRS(12) READ (44,*) (INTGRS(OFFSET+LOOP), LOOP = 1,INTGRS(13)) C C Check that all components are OK C DO 1 LOOP = 1,INTGRS(13) IF (INTGRS(INTGRS(12)+LOOP).LT.0) CALL ERRHAN(INTGRS, + 'Number of data points must be non-negative.','DATUM.',-1) COUNT = COUNT+INTGRS(INTGRS(12)+LOOP) 1 CONTINUE IF (COUNT.EQ.0) CALL ERRHAN(INTGRS, + 'No data points were specified.','DATUM.',-1) C INTGRS(64) = COUNT C C Be very sneaky and grab dataspace from the end of REALS, by C modifying its apparent length. C INTGRS(11) = INTGRS(11)-2*COUNT C C Calculate pointers into REALS for the TDATA and YDATA C TPOINT = INTGRS(11) YPOINT = TPOINT+COUNT INTGRS(66) = TPOINT C C Ensure there is enough space available, and read in the data C for each solution component C IF (INTGRS(11).LT.INTGRS(55)) + CALL ERRHAN(INTGRS,'Size of REALS is too small.','DATUM.',-1) C COUNT = 1 DO 3 LOOP = 1,INTGRS(13) C C Check to see if there is any data for this solution component C IF (INTGRS(INTGRS(12)+LOOP).EQ.0) GOTO 3 C C There is some data to read in, so read it into REALS C DO 2 LOOP2 = 1,INTGRS(INTGRS(12)+LOOP) READ (44,*) TDATA,YDATA REALS(TPOINT+COUNT) = TDATA REALS(YPOINT+COUNT) = YDATA COUNT = COUNT+1 2 CONTINUE C C Read the blank line which separates the data corresponding to C separate solution C IF (COUNT.LE.INTGRS(64)) READ (44,*) C 3 CONTINUE C CLOSE (UNIT=44) 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(REALS(INTGRS(62)),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,REALS(INTGRS(62)),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 C Do not allow the dimension of the problem to be specified after C the parameter estimation `data to be fitted' has been read in C IF (INTGRS(64).NE.0) CALL ERRHAN(INTGRS, + 'Dimension of problem cannot be specified after fitting data.', + 'DIMEN.',-1) C C Check that a valid dimension has been specified 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 Get the solution component from which discontinuities are propagated C FROM = INTGRS(INTMAP) 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(74) 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(21).LE.REALS(9)) THEN IF (UNDONE) INTGRS(74) = 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 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 C If parameter tolerance already set, issue message if error-per-step C tolerance is greater than parameter tolerance C IF (REALS(19).GT.0D0 .AND. EPS.GT.REALS(19)) CALL ERRHAN(INTGRS, + 'Error-per-step is greater than parameter tolerance.', + 'ERROR.',2) 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 Initialise 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 Initialise 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 bubblesort 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=77,RSTART=24) INTEGER DIMINT,DIMREL,INTGRS(*),LOOP DOUBLE PRECISION REALS(*),TEMP C EXTERNAL ERRHAN C C Routine to initialise default values and check storage allocation C Initialise 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 Initialise the number of HARD discontinuity networks C C INTGRS(8) = 0 C C Initialise the number of SOFT discontinuity networks C C INTGRS(9) = 0 C C Initialise 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 Initialise 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 Initialise the position in the cyclic storage of the first element C INTGRS(29) = 1 C C Initialise 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 Initialise 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 the number of parameters to be estimated C C INTGRS(61) = 0 C C Set the number of initial solution <-> parameter mappings C C INTGRS(65) = 0 C C Initialise number of bounds on parameters, non-linear constraints C and the pointer to them in REALS C C INTGRS(66) = 0 C INTGRS(67) = 0 C INTGRS(68) = 0 C C Set the `initial point is a parameter' flag/mapping to `no' C C INTGRS(70) = 0 C C Set kernel smoothness flag to `smooth' C C INTGRS(76) = 0 C C Further initialisation of INTGRS is postponed until any of the C default settings have been changed C C Initialise 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 C C Calculate the smallest positive number C REALS(22) = REALS(1) 5 REALS(22) = REALS(22)/2D0 TEMP = REALS(22)/2D0 WRITE (44,*) TEMP IF (TEMP.GT.0D0) GOTO 5 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 Ensure `initial' parameter tolerance is zero C C REALS(19) = 0D0 C C Set an upper bound for the size of the lag(s) (only used when tracking) C REALS(21) = 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,REALS(INTGRS(62)),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,REALS(INTGRS(62)),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,REALS(INTGRS(62)),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,REALS(INTGRS(62)),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,REALS(INTGRS(62)),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,REALS(INTGRS(62)),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,REALS(INTGRS(62)),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,REALS(INTGRS(62)),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,REALS(INTGRS(62)),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(21) = BOUND ENDIF C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE LIMITS(INTGRS,REALS,NONLIN) C INTEGER INTGRS(*),NONLIN DOUBLE PRECISION REALS(*) 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 (*,*) 'DATUM(INTGRS,REALS,filename)' WRITE (*,*) +' Routine to specify the file containing data to be fitted.' 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 (*,*) 'PARMAP(INTGRS,REALS,parameter,solution)' WRITE (*,*) +' Routine to specify initial solution value as a parameter.' WRITE (*,*) WRITE (*,*) 'PART0(INTGRS,REALS,parameter)' WRITE (*,*) +' Routine to specify that the initial point is a parameter.' WRITE (*,*) WRITE (*,*) 'PARTOL(INTGRS,REALS,tol)' WRITE (*,*) +' Routine to specify the tolerance for parameter estimation.' WRITE (*,*) WRITE (*,*) 'PDIMEN(INTGRS,REALS,parameters)' WRITE (*,*) +' Routine to specify number of parameters for parameter' WRITE (*,*) ' estimation.' 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 (*,*) 'PRNTRX(INTGRS,REALS)' WRITE (*,*) +' Routine to print out the residual vector in fitting data' WRITE (*,*) ' after parameter estimation.' 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 (*,*) 'RESID(INTGRS,REALS)' WRITE (*,*) +' Function to calculate the 2-norm of the residual.' 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,param,t,f)' WRITE (*,*) +' Routine to calculate all components of the derivative at `t`' WRITE (*,*) ' using the parameter estimation parameters in' WRITE (*,*) ' `param`, and return them in the vector `f`.' WRITE (*,*) WRITE (*,*) 'LAG(INTGRS,REALS,param,lagfn,t)' WRITE (*,*) +' Function to calculate the lag function `lagfn` at `t` -' WRITE (*,*) ' parameter estimation parameters in `param`.' WRITE (*,*) WRITE (*,*) 'YINIT(param,soln,t)' WRITE (*,*) +' Function to calculate the `soln`th component of the initial' WRITE (*,*) ' function at `t` using the parameter estimation' WRITE (*,*) ' parameters in `param`.' WRITE (*,*) WRITE (*,*) 'FINIT(param,deriv,t)' WRITE (*,*) +' Function to calculate the `soln`th component of the initial' WRITE (*,*) ' derivative at `t` using the parameter estimation' WRITE (*,*) ' parameters in `param`.' WRITE (*,*) WRITE (*,*) 'KERNEL(INTGRS,REALS,param,s,t,kernfn)' WRITE (*,*) +' Function to calculate the `kernfn`th kernel of a Volterra' WRITE (*,*) ' term at the points `s` and `t` using the' WRITE (*,*) ' parameter estimation parameters in `param`.' 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 initialise 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 SUBROUTINE MINIM(INTGRS,REALS,PARAM,Y) C INTEGER DUMMYI,FPOINT,INFO,INTGRS(*),LENGTH,LOOP,YPOINT DOUBLE PRECISION PARAM(*),REALS(*),Y(*) C EXTERNAL ERRHAN,LMDIF1,OPTIM C C Header routine to interface LMDIF1 and ARCHI C C Check that the number of parameters and parameter tolerance C have both been set C IF (INTGRS(61).LT.1) CALL ERRHAN(INTGRS, + 'Number of parameters has not been set.','MINIM.',-1) IF (INTGRS(64).LT.1) CALL ERRHAN(INTGRS, + 'Data to be fitted has not been given..','MINIM.',-1) IF (REALS(19).LT.REALS(1)) CALL ERRHAN(INTGRS, + 'Requested parameter tolerance not set.','MINIM.',-1) C C Calculate length of REAL work array required for minimization C LENGTH = (INTGRS(61)+1)*INTGRS(64)+5*INTGRS(61) C C Be sneaky again and use storage space at the end of INTGRS and C REALS to save declaring two sets of work arrays. Calculate pointers C to start of minimization work arrays, FVEC and backup Y in INTGRS C and REALS C INTGRS(11) = INTGRS(11)-LENGTH-INTGRS(64)-INTGRS(13)*2 INTGRS(12) = INTGRS(12)-INTGRS(61) C C Check if we still have enough storage space in INTGRS and REALS C IF (INTGRS(11).LT.INTGRS(55)-INTGRS(13)) + CALL ERRHAN(INTGRS,'Size of REALS is too small.','MINIM.',-1) IF (INTGRS(12).LT.INTGRS(39)) + CALL ERRHAN(INTGRS,'Size of INTGRS is too small.','MINIM.',-1) C C Store initial Y0 solution vector next to TDATA at end of REALS C YPOINT = INTGRS(11)+LENGTH+INTGRS(64) DO 1 LOOP = 1,INTGRS(13) REALS(YPOINT+LOOP) = Y(LOOP) 1 CONTINUE C C Calculate pointer to FVEC workspace next to Y0 at end of REALS C FPOINT = INTGRS(11)+LENGTH INTGRS(69) = FPOINT C C Do the minimization C CALL LMDIF1(OPTIM,INTGRS(64),INTGRS(61),PARAM,REALS(FPOINT+1), + REALS(19),INFO,INTGRS(INTGRS(12)+1), + REALS(INTGRS(11)+1),LENGTH,INTGRS,REALS) C C Re-solve the problem with the best-fit parameter values so that C the DDE solution is correct (if appropriate) C IF (INFO.NE.0) CALL OPTIM(INTGRS(64),INTGRS(61),PARAM, + REALS(FPOINT+1),DUMMYI,INTGRS,REALS) C C Restore final number of continuity records C INTGRS(10) = INTGRS(71) C C Only issue information about minimization results if requested. C Issue the appropriate message depending on INFO C GOTO (2,3,4,5,6,7,8,9) INFO+1 C 2 CALL ERRHAN(INTGRS,'Improper input parameters.','MINIM.',-1) C 3 IF (INTGRS(15).NE.0) THEN WRITE (*,*) 'Algorithm estimates the relative error in the sum' WRITE (*,*) 'of squares is at most ',REALS(19) ENDIF RETURN C 4 IF (INTGRS(15).NE.0) THEN WRITE (*,*) 'Algorithm estimates the relative error between x' WRITE (*,*) 'and the solution is at most ',REALS(19) ENDIF RETURN C 5 IF (INTGRS(15).NE.0) THEN WRITE (*,*) 'Algorithm estimates the relative error in the sum' WRITE (*,*) 'of squares and between x and the solution is at' WRITE (*,*) 'most ',REALS(19) ENDIF RETURN C 6 CALL ERRHAN(INTGRS, + 'Derivative vector is orthogonal to jacobian.','MINIM.',-1) C 7 CALL ERRHAN(INTGRS, + 'Number of Archi calls exceeded 200*(n+1).','MINIM.',-1) C 8 CALL ERRHAN(INTGRS, + 'No further reduction in the sum of squares is possible.', + 'MINIM.',2) RETURN C 9 CALL ERRHAN(INTGRS, + 'No further improvement in the solution is possible.', + 'MINIM.',2) C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE NBOUND(INTGRS,REALS,CONSTR,LOWER,UPPER) C INTEGER CONSTR,INTGRS(*) DOUBLE PRECISION LOWER,REALS(*),UPPER 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(REALS(INTGRS(62)),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(REALS(INTGRS(62)),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 SUBROUTINE OPTIM(DATUM,PARAMS,PARAM,FVEC,IFLAG,INTGRS,REALS) C DOUBLE PRECISION DABS,DELAY,DMAX1,DSQRT,YINIT C INTEGER COUNT,DATUM,IFLAG,LOOP,LOOP2,NUMBER,PARAMS INTEGER CPOINT,PPOINT,TDATA,YDATA,YPOINT INTEGER INTGRS(*) DOUBLE PRECISION FVEC(DATUM),PARAM(*),REALS(*) DOUBLE PRECISION DIST,NORM2,PRTURB C EXTERNAL DELAY,SOLVE,YINIT INTRINSIC DABS,DMAX1,DSQRT C C Calculate the position of the backup initial value in REALS C YPOINT = INTGRS(11)+(INTGRS(61)+2)*INTGRS(64)+INTGRS(61)*5 C C Copy the initial value into a new pseudo-Y vector which exists C in REALS. The pseudo-Y vector is changed after returning from C SOLVE. The original Y vector cannot be used because it cannot C be passed through LMDIF (or the NAG routine). C DO 1 LOOP = 1,INTGRS(13) REALS(YPOINT+INTGRS(13)+LOOP) = REALS(YPOINT+LOOP) 1 CONTINUE C C If the initial point is a parameter then modify the start point, C the first stored interval, the current point BUT FIRST modify C the continuity record. C PRTURB = 5D0*DMAX1(1D0,DABS(REALS(7)))*REALS(1) IF (INTGRS(70).NE.0) THEN DO 2 LOOP = 0,INTGRS(10)-1 DIST = DABS(REALS(INTGRS(55)+LOOP)-REALS(7)) IF (DIST.LE.PRTURB) + REALS(INTGRS(55)+LOOP) = PARAM(INTGRS(70)) 2 CONTINUE REALS(3) = PARAM(INTGRS(70)) REALS(7) = PARAM(INTGRS(70)) REALS(9) = PARAM(INTGRS(70)) ENDIF C C Check the user-supplied continuity record and if the C continuity at the initial point is not zero, then map the C appropriate value of the initial function into the initial C solution (overwriting any user-specified initial solution) C CPOINT = INTGRS(39) DO 3 LOOP = 0,INTGRS(10)-1 IF (INTGRS(CPOINT+1).GT.0) THEN C C Calculate distance of initial point from discontinuity C DIST = DABS(REALS(INTGRS(55)+LOOP)-REALS(7)) IF (DIST.LE.PRTURB) + REALS(YPOINT+INTGRS(13)+INTGRS(CPOINT)) = + YINIT(PARAM,INTGRS(CPOINT),REALS(7)) ENDIF CPOINT = CPOINT+2 3 CONTINUE C C Now modify the initial solution if any initial solution <-> C parameter mappings have been specified C IF (INTGRS(65).NE.0) THEN PPOINT = INTGRS(12)+INTGRS(61)+INTGRS(13)+1 DO 4 LOOP = 1,INTGRS(65) REALS(YPOINT+INTGRS(13)+INTGRS(PPOINT)) = + PARAM(INTGRS(PPOINT+1)) PPOINT = PPOINT+2 4 CONTINUE ENDIF C C (re)set the pointers to first/last interval stored and then C solve the DDE C INTGRS(29) = 1 INTGRS(30) = 0 CALL SOLVE(INTGRS,REALS,PARAM,REALS(YPOINT+INTGRS(13)+1)) C C Calculate the position of the TDATA and YDATA in REALS C TDATA = INTGRS(66) YDATA = TDATA+INTGRS(64) C C calculate the error in fitting the data for each component of C the solution and the 2-norm of the residual vector C COUNT = 1 NORM2 = 0D0 DO 6 LOOP = 1,INTGRS(13) C C Check (and read) if any data to fitted for this component C NUMBER = INTGRS(INTGRS(12)+INTGRS(61)+LOOP) IF (NUMBER.EQ.0) GOTO 6 C DO 5 LOOP2 = 1,NUMBER REALS(INTGRS(45)) = REALS(TDATA+COUNT) FVEC(COUNT) = REALS(YDATA+COUNT)-DELAY(INTGRS,REALS,1,LOOP) NORM2 = NORM2+FVEC(COUNT)*FVEC(COUNT) COUNT = COUNT+1 5 CONTINUE 6 CONTINUE C C Save the 2-norm where it can be readily accessed by the function RESID C REALS(20) = DSQRT(NORM2) C C restore the renewable Archi variables C 1) Initial number of continuity records set (but keep a copy) C 2) Total number of steps C 3) Number of accepted steps C 4) Initial point C 5) Initial stepsize C INTGRS(71) = INTGRS(10) INTGRS(10) = INTGRS(63) INTGRS(22) = 0 INTGRS(23) = 0 C REALS(9) = REALS(7) REALS(11) = REALS(18) C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE PARMAP(INTGRS,REALS,PARAM,SOLN) C INTEGER INTGRS(*),LOOP,PARAM,SOLN DOUBLE PRECISION REALS(*) C EXTERNAL ERRHAN C C Routine to specify which components of the initial solution are C actually parameters C C Ensure that the dimension of the solution and number of parameters C have both been specified. Also ensure that the data to be fitted C has not yet been specified C IF (INTGRS(61).EQ.0) CALL ERRHAN(INTGRS, + 'Number of parameters has not been specified.','PARMAP.',-1) IF (INTGRS(64).NE.0) CALL ERRHAN(INTGRS, + 'Parameter mappings cannot be specified after fitting data.', + 'PARMAP.',-1) C C Check that the solution component and parameter are both valid C IF (SOLN.LT.1 .OR. SOLN.GT.INTGRS(13)) CALL ERRHAN(INTGRS, + 'Invalid solution component specified.','PARMAP.',-1) IF (PARAM.LT.1 .OR. PARAM.GT.INTGRS(61)) CALL ERRHAN(INTGRS, + 'Invalid parameter number specified.','PARMAP.',-1) C C Check that the solution component has not already been mapped C IF (INTGRS(65).NE.0) THEN DO 1 LOOP = 0,INTGRS(65)-1 IF (INTGRS(INTGRS(12)+LOOP*2+1).EQ.SOLN) THEN CALL ERRHAN(INTGRS,'Solution component already mapped.', + 'PARMAP.',2) IF (INTGRS(15).GE.2) WRITE (*,*) + 'Request to map initial solution component to parameter ignored' RETURN ENDIF 1 CONTINUE ENDIF C C All OK, so sneakily add the mapping (at the end of INTGRS array) C and modify the apparent length of INTGRS accordingly. Order of C data is SOLN then PARAM C INTGRS(12) = INTGRS(12)-2 INTGRS(65) = INTGRS(65)+1 INTGRS(INTGRS(12)+1) = SOLN INTGRS(INTGRS(12)+2) = PARAM C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE PART0(INTGRS,REALS,MAPTO) C INTEGER INTGRS(*),MAPTO DOUBLE PRECISION REALS(*) C EXTERNAL ERRHAN C C Routine to specify that the initial point is a parameter, and to C map the appropriate parameter C C Ensure that the dimension of the solution and number of parameters C have both been specified. C IF (INTGRS(61).EQ.0) CALL ERRHAN(INTGRS, + 'Number of parameters has not been specified.','PART0.',-1) C C Check that PART0 has not already been called C IF (INTGRS(70).NE.0) CALL ERRHAN(INTGRS, + 'Parameter mapping to initial point cannot be changed.', + 'PART0.',-1) C C Check that the parameter is valid C IF (MAPTO.LT.1 .OR. MAPTO.GT.INTGRS(61)) CALL ERRHAN(INTGRS, + 'Invalid parameter number specified.','PART0.',-1) C INTGRS(70) = MAPTO C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE PARTOL(INTGRS,REALS,TOL) C INTEGER INTGRS(*) DOUBLE PRECISION REALS(*),TOL C EXTERNAL ERRHAN C C Routine to specify the tolerance to which parameters are to found C IF (TOL.LT.100D0*REALS(1)) CALL ERRHAN(INTGRS, + 'Parameter tolerance is too small.','PARTOL.',-1) C C Check if parameter tolerance already set C IF (REALS(19).GT.0D0) CALL ERRHAN(INTGRS, + 'Parameter tolerance already set.','PARTOL.',2) C C If DDE tolerance set, issue message if parameter tolerance is C less than the error-per-step tolerance C IF (INTGRS(6).GT.0 .AND. TOL.LT.REALS(2) .AND. INTGRS(15).GT.1) + CALL ERRHAN(INTGRS, + 'Parameter tolerance is less than error-per-step tolerance.', + 'PARTOL.',2) C REALS(19) = TOL C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE PDIMEN(INTGRS,REALS,NUMBER) C INTEGER INTGRS(*),NUMBER DOUBLE PRECISION REALS(*) C EXTERNAL ERRHAN C C Routine to specify the number of parameters to be estimated C IF (NUMBER.LT.1) CALL ERRHAN(INTGRS, + 'Number of parameters must be positive.','PDIMEN.',-1) IF (INTGRS(61).NE.0) CALL ERRHAN(INTGRS, + 'Number of parameters already set.','PDIMEN.',2) C INTGRS(61) = NUMBER 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 PRNTRX(INTGRS,REALS) C INTEGER INTGRS(*),LOOP DOUBLE PRECISION REALS(*) C EXTERNAL ERRHAN C C Routine to print out the residual vector in fitting data for C parameter estimation C C Check that LMDIF/E04UPF has been called at least once, by checking C the number of intervals solved C IF (INTGRS(25).EQ.0) THEN CALL ERRHAN(INTGRS, + 'The residual is not available.','PRNTRX.',2) RETURN ENDIF C WRITE (*,*) 'Residual vector =' DO 1 LOOP = 1,INTGRS(64) WRITE (*,9999) REALS(INTGRS(66)+LOOP),REALS(INTGRS(69)+LOOP) 1 CONTINUE C 9999 FORMAT (2G25.17) 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 DBLE,DELAY 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 indexs 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(23) .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.', + 'MINIM.',-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.', + 'MINIM.',-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 information/warnings 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 DOUBLE PRECISION FUNCTION RESID(INTGRS,REALS) C INTEGER INTGRS(*) DOUBLE PRECISION REALS(*) C EXTERNAL ERRHAN C C Function to return the (current) value of the 2-norm of the C residual vector C C Check that LMDIF/E04UPF has been called at least once, by checking C the number of intervals solved C IF (INTGRS(25).EQ.0) THEN CALL ERRHAN(INTGRS, + 'The residual is not available.','RESID.',2) RESID = 0D0 RETURN ENDIF C RESID = REALS(20) 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 Initialise 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(INTGRS(62)),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.', + 'MINIM.',-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 a 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(23) = 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,REALS(INTGRS(62)),NEWX,T,INTGRS(58)) NEWX = TTRANS(H,X,-DSQRT((3D0-DSQRT(24D0/5D0))/7D0)) ADD = ADD+(1D0+DSQRT(5D0/54D0))* + KERNEL(INTGRS,REALS,REALS(INTGRS(62)),NEWX,T,INTGRS(58)) NEWX = TTRANS(H,X,DSQRT((3D0-DSQRT(24D0/5D0))/7D0)) ADD = ADD+(1D0+DSQRT(5D0/54D0))* + KERNEL(INTGRS,REALS,REALS(INTGRS(62)),NEWX,T,INTGRS(58)) NEWX = TTRANS(H,X,DSQRT((3D0+DSQRT(24D0/5D0))/7D0)) ADD = ADD+(1D0-DSQRT(5D0/54D0))* + KERNEL(INTGRS,REALS,REALS(INTGRS(62)),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,REALS(INTGRS(62)),NEWX,T,INTGRS(58)) NEWX = TTRANS(H,X,0D0) ADD = ADD+8D0* + KERNEL(INTGRS,REALS,REALS(INTGRS(62)),NEWX,T,INTGRS(58)) NEWX = TTRANS(H,X,DSQRT(3D0/5D0)) ADD = ADD+5D0* + KERNEL(INTGRS,REALS,REALS(INTGRS(62)),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,REALS(INTGRS(62)),NEWX,T,INTGRS(58)) NEWX = TTRANS(H,X,DSQRT(1D0/3D0)) ADD = ADD+ + KERNEL(INTGRS,REALS,REALS(INTGRS(62)),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(73) 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(73) = 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,REALS(INTGRS(62)),LAGFN,FIRST)-DISCAT FSECND = LAG(INTGRS,REALS,REALS(INTGRS(62)),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,REALS(INTGRS(62)),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,REALS(INTGRS(62)),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,REALS(INTGRS(62)),LAGFN,SECND)-DISCAT ELSE C C State dependent case CC 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,REALS(INTGRS(62)),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(76))) 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(76).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(76).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(75) 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(75) = INTGRS(INTGRS(42)) INTGRS(INTGRS(42)) = OLDPNT REALS(INTGRS(45)) = OLDLAG C RETURN END C C C ******************************************************************+++++++ C C SUBROUTINE SOLVE(INTGRS,REALS,PARAM,Y) C DOUBLE PRECISION DABS,DBLE,DMAX1 C INTEGER INTGRS(*),LOOP DOUBLE PRECISION ENDPNT,REALS(*),PARAM(*),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.', + 'MINIM.',-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.', + 'MINIM.',-1) IF (INTGRS(5).GT.0 .AND. INTGRS(10).EQ.0) + CALL ERRHAN(INTGRS, + 'No continuity records specified for tracking.', + 'MINIM.',-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.', + 'MINIM.',-1) C C Initialise any variables that were not set before because C of a dependence on the nature of the problem C C Initialise 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,'MINIM.') C C Initialise the storage variable `last point stored' and a perturbed C version, used in improving the performance of RECALL C REALS(3) = REALS(7) REALS(23) = REALS(3)+DMAX1(DABS(REALS(3)),1D0)*REALS(1) C C Initialise diagnostics HMAXM, HMINM, EXTRAP C REALS(4) = 0D0 REALS(5) = 2D0*DMAX1(REALS(11),REALS(12)) REALS(6) = 1D0 C C Initialise `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 parameter estimation parameters in REALS C The total number of stages in the Runge-Kutta method is 9 C INTGRS(62) = INTGRS(37)+INTGRS(18)*9 C C Calculate the start of the cyclic storage structure in REALS C INTGRS(38) = INTGRS(62)+INTGRS(61) 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 Only print length of cyclic storage once for minimization problems C IF (INTGRS(33).LT.10) + CALL ERRHAN(INTGRS,'Size of REALS is too small.','MINIM.',-1) IF (INTGRS(15).GE.2 .AND. INTGRS(71).EQ.0) + 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 Copy the parameter estimates into REALS, and store the initial C stepsize and initial number of continuity records set - for use C in parameter estimation problems C DO 2 LOOP = 1,INTGRS(61) REALS(INTGRS(62)+LOOP-1) = PARAM(LOOP) 2 CONTINUE C INTGRS(63) = INTGRS(10) REALS(18) = REALS(11) C C Transfer the initial solution to Y and Y1 in REALS C DO 3 LOOP = 1,INTGRS(13) REALS(INTGRS(34)+LOOP-1) = Y(LOOP) REALS(INTGRS(35)+LOOP-1) = Y(LOOP) 3 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(73) = 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(74) = 1 C C Pointer into the continuity record used by SOLNYI C INTGRS(75) = 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 4 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 4 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 5 LOOP = 1,INTGRS(13) Y(LOOP) = REALS(INTGRS(34)+LOOP-1) 5 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(23) = 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.1) 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(76) = 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.', + 'MINIM.',-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 REALISE 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 Initialise 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 C C C ******************************************************************+++++++ C ******************************************************************+++++++ C ******************************************************************+++++++ C C subroutine lmdif1(fcn,m,n,x,fvec,tol,info,iwa,wa,lwa,INTGRS,REALS) integer m,n,info,lwa integer iwa(n) double precision tol double precision x(n),fvec(m),wa(lwa) INTEGER INTGRS(*) DOUBLE PRECISION REALS(*) external fcn c ********** c c subroutine lmdif1 c c the purpose of lmdif1 is to minimize the sum of the squares of c m nonlinear functions in n variables by a modification of the c levenberg-marquardt algorithm. this is done by using the more c general least-squares solver lmdif. the user must provide a c subroutine which calculates the functions. the jacobian is c then calculated by a forward-difference approximation. c c the subroutine statement is c c subroutine lmdif1(fcn,m,n,x,fvec,tol,info,iwa,wa,lwa,INTGRS,REALS) c c where c c fcn is the name of the user-supplied subroutine which c calculates the functions. fcn must be declared c in an external statement in the user calling c program, and should be written as follows. c c subroutine fcn(m,n,x,fvec,iflag,INTGRS,REALS) c integer m,n,iflag c double precision x(n),fvec(m) c ---------- c calculate the functions at x and c return this vector in fvec. c ---------- c return c end c c the value of iflag should not be changed by fcn unless c the user wants to terminate execution of lmdif1. c in this case set iflag to a negative integer. c c m is a positive integer input variable set to the number c of functions. c c n is a positive integer input variable set to the number c of variables. n must not exceed m. c c x is an array of length n. on input x must contain c an initial estimate of the solution vector. on output x c contains the final estimate of the solution vector. c c fvec is an output array of length m which contains c the functions evaluated at the output x. c c tol is a nonnegative input variable. termination occurs c when the algorithm estimates either that the relative c error in the sum of squares is at most tol or that c the relative error between x and the solution is at c most tol. c c info is an integer output variable. if the user has c terminated execution, info is set to the (negative) c value of iflag. see description of fcn. otherwise, c info is set as follows. c c info = 0 improper input parameters. c c info = 1 algorithm estimates that the relative error c in the sum of squares is at most tol. c c info = 2 algorithm estimates that the relative error c between x and the solution is at most tol. c c info = 3 conditions for info = 1 and info = 2 both hold. c c info = 4 fvec is orthogonal to the columns of the c jacobian to machine precision. c c info = 5 number of calls to fcn has reached or c exceeded 200*(n+1). c c info = 6 tol is too small. no further reduction in c the sum of squares is possible. c c info = 7 tol is too small. no further improvement in c the approximate solution x is possible. c c iwa is an integer work array of length n. c c wa is a work array of length lwa. c c lwa is a positive integer input variable not less than c m*n+5*n+m. c c subprograms called c c user-supplied ...... fcn c c minpack-supplied ... lmdif c c argonne national laboratory. minpack project. march 1980. c burton s. garbow, kenneth e. hillstrom, jorge j. more c c ********** integer maxfev,mode,mp5n,nfev,nprint double precision epsfcn,factor,ftol,gtol,xtol,zero data factor,zero /10d0,0.0d0/ info = 0 c c check the input parameters for errors. c if (n .le. 0 .or. m .lt. n .or. tol .lt. zero * .or. lwa .lt. m*n + 5*n + m) go to 10 c c call lmdif. c maxfev = 200*(n + 1) ftol = tol xtol = tol gtol = zero epsfcn = zero mode = 1 nprint = 0 mp5n = m + 5*n call lmdif(fcn,m,n,x,fvec,ftol,xtol,gtol,maxfev,epsfcn,wa(1), * mode,factor,nprint,info,nfev,wa(mp5n+1),m,iwa, * wa(n+1),wa(2*n+1),wa(3*n+1),wa(4*n+1),wa(5*n+1), + INTGRS,REALS) if (info .eq. 8) info = 4 10 continue return c c last card of subroutine lmdif1. c end subroutine lmdif(fcn,m,n,x,fvec,ftol,xtol,gtol,maxfev,epsfcn, * diag,mode,factor,nprint,info,nfev,fjac,ldfjac, * ipvt,qtf,wa1,wa2,wa3,wa4,INTGRS,REALS) integer m,n,maxfev,mode,nprint,info,nfev,ldfjac integer ipvt(n) double precision ftol,xtol,gtol,epsfcn,factor double precision x(n),fvec(m),diag(n),fjac(ldfjac,n),qtf(n), * wa1(n),wa2(n),wa3(n),wa4(m) INTEGER INTGRS(*) DOUBLE PRECISION REALS(*) external fcn c ********** c c subroutine lmdif c c the purpose of lmdif is to minimize the sum of the squares of c m nonlinear functions in n variables by a modification of c the levenberg-marquardt algorithm. the user must provide a c subroutine which calculates the functions. the jacobian is c then calculated by a forward-difference approximation. c c the subroutine statement is c c subroutine lmdif(fcn,m,n,x,fvec,ftol,xtol,gtol,maxfev,epsfcn, c diag,mode,factor,nprint,info,nfev,fjac, c ldfjac,ipvt,qtf,wa1,wa2,wa3,wa4,INTGRS,REALS) c c where c c fcn is the name of the user-supplied subroutine which c calculates the functions. fcn must be declared c in an external statement in the user calling c program, and should be written as follows. c c subroutine fcn(m,n,x,fvec,iflag,INTGRS,REALS) c integer m,n,iflag c double precision x(n),fvec(m) c ---------- c calculate the functions at x and c return this vector in fvec. c ---------- c return c end c c the value of iflag should not be changed by fcn unless c the user wants to terminate execution of lmdif. c in this case set iflag to a negative integer. c c m is a positive integer input variable set to the number c of functions. c c n is a positive integer input variable set to the number c of variables. n must not exceed m. c c x is an array of length n. on input x must contain c an initial estimate of the solution vector. on output x c contains the final estimate of the solution vector. c c fvec is an output array of length m which contains c the functions evaluated at the output x. c c ftol is a nonnegative input variable. termination c occurs when both the actual and predicted relative c reductions in the sum of squares are at most ftol. c therefore, ftol measures the relative error desired c in the sum of squares. c c xtol is a nonnegative input variable. termination c occurs when the relative error between two consecutive c iterates is at most xtol. therefore, xtol measures the c relative error desired in the approximate solution. c c gtol is a nonnegative input variable. termination c occurs when the cosine of the angle between fvec and c any column of the jacobian is at most gtol in absolute c value. therefore, gtol measures the orthogonality c desired between the function vector and the columns c of the jacobian. c c maxfev is a positive integer input variable. termination c occurs when the number of calls to fcn is at least c maxfev by the end of an iteration. c c epsfcn is an input variable used in determining a suitable c step length for the forward-difference approximation. this c approximation assumes that the relative errors in the c functions are of the order of epsfcn. if epsfcn is less c than the machine precision, it is assumed that the relative c errors in the functions are of the order of the machine c precision. c c diag is an array of length n. if mode = 1 (see c below), diag is internally set. if mode = 2, diag c must contain positive entries that serve as c multiplicative scale factors for the variables. c c mode is an integer input variable. if mode = 1, the c variables will be scaled internally. if mode = 2, c the scaling is specified by the input diag. other c values of mode are equivalent to mode = 1. c c factor is a positive input variable used in determining the c initial step bound. this bound is set to the product of c factor and the euclidean norm of diag*x if nonzero, or else c to factor itself. in most cases factor should lie in the c interval (.1,100.). 100. is a generally recommended value. c c nprint is an integer input variable that enables controlled c printing of iterates if it is positive. in this case, c fcn is called with iflag = 0 at the beginning of the first c iteration and every nprint iterations thereafter and c immediately prior to return, with x and fvec available c for printing. if nprint is not positive, no special calls c of fcn with iflag = 0 are made. c c info is an integer output variable. if the user has c terminated execution, info is set to the (negative) c value of iflag. see description of fcn. otherwise, c info is set as follows. c c info = 0 improper input parameters. c c info = 1 both actual and predicted relative reductions c in the sum of squares are at most ftol. c c info = 2 relative error between two consecutive iterates c is at most xtol. c c info = 3 conditions for info = 1 and info = 2 both hold. c c info = 4 the cosine of the angle between fvec and any c column of the jacobian is at most gtol in c absolute value. c c info = 5 number of calls to fcn has reached or c exceeded maxfev. c c info = 6 ftol is too small. no further reduction in c the sum of squares is possible. c c info = 7 xtol is too small. no further improvement in c the approximate solution x is possible. c c info = 8 gtol is too small. fvec is orthogonal to the c columns of the jacobian to machine precision. c c nfev is an integer output variable set to the number of c calls to fcn. c c fjac is an output m by n array. the upper n by n submatrix c of fjac contains an upper triangular matrix r with c diagonal elements of nonincreasing magnitude such that c c t t t c p *(jac *jac)*p = r *r, c c where p is a permutation matrix and jac is the final c calculated jacobian. column j of p is column ipvt(j) c (see below) of the identity matrix. the lower trapezoidal c part of fjac contains information generated during c the computation of r. c c ldfjac is a positive integer input variable not less than m c which specifies the leading dimension of the array fjac. c c ipvt is an integer output array of length n. ipvt c defines a permutation matrix p such that jac*p = q*r, c where jac is the final calculated jacobian, q is c orthogonal (not stored), and r is upper triangular c with diagonal elements of nonincreasing magnitude. c column j of p is column ipvt(j) of the identity matrix. c c qtf is an output array of length n which contains c the first n elements of the vector (q transpose)*fvec. c c wa1, wa2, and wa3 are work arrays of length n. c c wa4 is a work array of length m. c c subprograms called c c user-supplied ...... fcn c c minpack-supplied ... enorm,fdjac2,lmpar,qrfac c c fortran-supplied ... dabs,dmax1,dmin1,dsqrt,mod c c argonne national laboratory. minpack project. march 1980. c burton s. garbow, kenneth e. hillstrom, jorge j. more c c ********** integer i,iflag,iter,j,l double precision actred,delta,dirder,epsmch,fnorm,fnorm1,gnorm, * one,par,pnorm,prered,p1,p5,p25,p75,p0001,ratio, * sum,temp,temp1,temp2,xnorm,zero double precision enorm data one,p1,p5,p25,p75,p0001,zero * /1.0d0,1.0d-1,5.0d-1,2.5d-1,7.5d-1,1.0d-4,0.0d0/ c c epsmch is the machine precision. c epsmch = REALS(1) c info = 0 iflag = 0 nfev = 0 c c check the input parameters for errors. c if (n .le. 0 .or. m .lt. n .or. ldfjac .lt. m * .or. ftol .lt. zero .or. xtol .lt. zero .or. gtol .lt. zero * .or. maxfev .le. 0 .or. factor .le. zero) go to 300 if (mode .ne. 2) go to 20 do 10 j = 1, n if (diag(j) .le. zero) go to 300 10 continue 20 continue c c evaluate the function at the starting point c and calculate its norm. c iflag = 1 call fcn(m,n,x,fvec,iflag,INTGRS,REALS) nfev = 1 if (iflag .lt. 0) go to 300 fnorm = enorm(m,fvec) c c initialize levenberg-marquardt parameter and iteration counter. c par = zero iter = 1 c c beginning of the outer loop. c 30 continue c c calculate the jacobian matrix. c iflag = 2 call fdjac2(fcn,m,n,x,fvec,fjac,ldfjac,iflag,epsfcn,wa4, + INTGRS,REALS) nfev = nfev + n if (iflag .lt. 0) go to 300 c c if requested, call fcn to enable printing of iterates. c if (nprint .le. 0) go to 40 iflag = 0 if (mod(iter-1,nprint) .eq. 0) call fcn(m,n,x,fvec,iflag, + INTGRS,REALS) if (iflag .lt. 0) go to 300 40 continue c c compute the qr factorization of the jacobian. c call qrfac(m,n,fjac,ldfjac,.true.,ipvt,n,wa1,wa2,wa3, + INTGRS,REALS) c c on the first iteration and if mode is 1, scale according c to the norms of the columns of the initial jacobian. c if (iter .ne. 1) go to 80 if (mode .eq. 2) go to 60 do 50 j = 1, n diag(j) = wa2(j) if (wa2(j) .eq. zero) diag(j) = one 50 continue 60 continue c c on the first iteration, calculate the norm of the scaled x c and initialize the step bound delta. c do 70 j = 1, n wa3(j) = diag(j)*x(j) 70 continue xnorm = enorm(n,wa3) delta = factor*xnorm if (delta .eq. zero) delta = factor 80 continue c c form (q transpose)*fvec and store the first n components in c qtf. c do 90 i = 1, m wa4(i) = fvec(i) 90 continue do 130 j = 1, n if (fjac(j,j) .eq. zero) go to 120 sum = zero do 100 i = j, m sum = sum + fjac(i,j)*wa4(i) 100 continue temp = -sum/fjac(j,j) do 110 i = j, m wa4(i) = wa4(i) + fjac(i,j)*temp 110 continue 120 continue fjac(j,j) = wa1(j) qtf(j) = wa4(j) 130 continue c c compute the norm of the scaled gradient. c gnorm = zero if (fnorm .eq. zero) go to 170 do 160 j = 1, n l = ipvt(j) if (wa2(l) .eq. zero) go to 150 sum = zero do 140 i = 1, j sum = sum + fjac(i,j)*(qtf(i)/fnorm) 140 continue gnorm = dmax1(gnorm,dabs(sum/wa2(l))) 150 continue 160 continue 170 continue c c test for convergence of the gradient norm. c if (gnorm .le. gtol) info = 4 if (info .ne. 0) go to 300 c c rescale if necessary. c if (mode .eq. 2) go to 190 do 180 j = 1, n diag(j) = dmax1(diag(j),wa2(j)) 180 continue 190 continue c c beginning of the inner loop. c 200 continue c c determine the levenberg-marquardt parameter. c call lmpar(n,fjac,ldfjac,ipvt,diag,qtf,delta,par,wa1,wa2, * wa3,wa4,INTGRS,REALS) c c store the direction p and x + p. calculate the norm of p. c do 210 j = 1, n wa1(j) = -wa1(j) wa2(j) = x(j) + wa1(j) wa3(j) = diag(j)*wa1(j) 210 continue pnorm = enorm(n,wa3) c c on the first iteration, adjust the initial step bound. c if (iter .eq. 1) delta = dmin1(delta,pnorm) c c evaluate the function at x + p and calculate its norm. c iflag = 1 call fcn(m,n,wa2,wa4,iflag,INTGRS,REALS) nfev = nfev + 1 if (iflag .lt. 0) go to 300 fnorm1 = enorm(m,wa4) c c compute the scaled actual reduction. c actred = -one if (p1*fnorm1 .lt. fnorm) actred = one - (fnorm1/fnorm)**2 c c compute the scaled predicted reduction and c the scaled directional derivative. c do 230 j = 1, n wa3(j) = zero l = ipvt(j) temp = wa1(l) do 220 i = 1, j wa3(i) = wa3(i) + fjac(i,j)*temp 220 continue 230 continue temp1 = enorm(n,wa3)/fnorm temp2 = (dsqrt(par)*pnorm)/fnorm prered = temp1**2 + temp2**2/p5 dirder = -(temp1**2 + temp2**2) c c compute the ratio of the actual to the predicted c reduction. c ratio = zero if (prered .ne. zero) ratio = actred/prered c c update the step bound. c if (ratio .gt. p25) go to 240 if (actred .ge. zero) temp = p5 if (actred .lt. zero) * temp = p5*dirder/(dirder + p5*actred) if (p1*fnorm1 .ge. fnorm .or. temp .lt. p1) temp = p1 delta = temp*dmin1(delta,pnorm/p1) par = par/temp go to 260 240 continue if (par .ne. zero .and. ratio .lt. p75) go to 250 delta = pnorm/p5 par = p5*par 250 continue 260 continue c c test for successful iteration. c if (ratio .lt. p0001) go to 290 c c successful iteration. update x, fvec, and their norms. c do 270 j = 1, n x(j) = wa2(j) wa2(j) = diag(j)*x(j) 270 continue do 280 i = 1, m fvec(i) = wa4(i) 280 continue xnorm = enorm(n,wa2) fnorm = fnorm1 iter = iter + 1 290 continue c c tests for convergence. c if (dabs(actred) .le. ftol .and. prered .le. ftol * .and. p5*ratio .le. one) info = 1 if (delta .le. xtol*xnorm) info = 2 if (dabs(actred) .le. ftol .and. prered .le. ftol * .and. p5*ratio .le. one .and. info .eq. 2) info = 3 if (info .ne. 0) go to 300 c c tests for termination and stringent tolerances. c if (nfev .ge. maxfev) info = 5 if (dabs(actred) .le. epsmch .and. prered .le. epsmch * .and. p5*ratio .le. one) info = 6 if (delta .le. epsmch*xnorm) info = 7 if (gnorm .le. epsmch) info = 8 if (info .ne. 0) go to 300 c c end of the inner loop. repeat if iteration unsuccessful. c if (ratio .lt. p0001) go to 200 c c end of the outer loop. c go to 30 300 continue c c termination, either normal or user imposed. c if (iflag .lt. 0) info = iflag iflag = 0 if (nprint .gt. 0) call fcn(m,n,x,fvec,iflag,INTGRS,REALS) C C Save the number of iterations C INTGRS(72) = iter C return c c last card of subroutine lmdif. c end subroutine lmpar(n,r,ldr,ipvt,diag,qtb,delta,par,x,sdiag,wa1, * wa2,INTGRS,REALS) integer n,ldr integer ipvt(n) double precision delta,par double precision r(ldr,n),diag(n),qtb(n),x(n),sdiag(n),wa1(n), * wa2(n) INTEGER INTGRS(*) DOUBLE PRECISION REALS(*) c ********** c c subroutine lmpar c c given an m by n matrix a, an n by n nonsingular diagonal c matrix d, an m-vector b, and a positive number delta, c the problem is to determine a value for the parameter c par such that if x solves the system c c a*x = b , sqrt(par)*d*x = 0 , c c in the least squares sense, and dxnorm is the euclidean c norm of d*x, then either par is zero and c c (dxnorm-delta) .le. 0.1*delta , c c or par is positive and c c abs(dxnorm-delta) .le. 0.1*delta . c c this subroutine completes the solution of the problem c if it is provided with the necessary information from the c qr factorization, with column pivoting, of a. that is, if c a*p = q*r, where p is a permutation matrix, q has orthogonal c columns, and r is an upper triangular matrix with diagonal c elements of nonincreasing magnitude, then lmpar expects c the full upper triangle of r, the permutation matrix p, c and the first n components of (q transpose)*b. on output c lmpar also provides an upper triangular matrix s such that c c t t t c p *(a *a + par*d*d)*p = s *s . c c s is employed within lmpar and may be of separate interest. c c only a few iterations are generally needed for convergence c of the algorithm. if, however, the limit of 10 iterations c is reached, then the output par will contain the best c value obtained so far. c c the subroutine statement is c c subroutine lmpar(n,r,ldr,ipvt,diag,qtb,delta,par,x,sdiag, c wa1,wa2,INTGRS,REALS) c c where c c n is a positive integer input variable set to the order of r. c c r is an n by n array. on input the full upper triangle c must contain the full upper triangle of the matrix r. c on output the full upper triangle is unaltered, and the c strict lower triangle contains the strict upper triangle c (transposed) of the upper triangular matrix s. c c ldr is a positive integer input variable not less than n c which specifies the leading dimension of the array r. c c ipvt is an integer input array of length n which defines the c permutation matrix p such that a*p = q*r. column j of p c is column ipvt(j) of the identity matrix. c c diag is an input array of length n which must contain the c diagonal elements of the matrix d. c c qtb is an input array of length n which must contain the first c n elements of the vector (q transpose)*b. c c delta is a positive input variable which specifies an upper c bound on the euclidean norm of d*x. c c par is a nonnegative variable. on input par contains an c initial estimate of the levenberg-marquardt parameter. c on output par contains the final estimate. c c x is an output array of length n which contains the least c squares solution of the system a*x = b, sqrt(par)*d*x = 0, c for the output par. c c sdiag is an output array of length n which contains the c diagonal elements of the upper triangular matrix s. c c wa1 and wa2 are work arrays of length n. c c subprograms called c c minpack-supplied ... enorm,qrsolv c c fortran-supplied ... dabs,dmax1,dmin1,dsqrt c c argonne national laboratory. minpack project. march 1980. c burton s. garbow, kenneth e. hillstrom, jorge j. more c c ********** integer i,iter,j,jm1,jp1,k,l,nsing double precision dxnorm,dwarf,fp,gnorm,parc,parl,paru,p1,p001, * sum,temp,zero double precision enorm data p1,p001,zero /1.0d-1,1.0d-3,0.0d0/ c c dwarf is the smallest positive magnitude. c dwarf = REALS(22) c c compute and store in x the gauss-newton direction. if the c jacobian is rank-deficient, obtain a least squares solution. c nsing = n do 10 j = 1, n wa1(j) = qtb(j) if (r(j,j) .eq. zero .and. nsing .eq. n) nsing = j - 1 if (nsing .lt. n) wa1(j) = zero 10 continue if (nsing .lt. 1) go to 50 do 40 k = 1, nsing j = nsing - k + 1 wa1(j) = wa1(j)/r(j,j) temp = wa1(j) jm1 = j - 1 if (jm1 .lt. 1) go to 30 do 20 i = 1, jm1 wa1(i) = wa1(i) - r(i,j)*temp 20 continue 30 continue 40 continue 50 continue do 60 j = 1, n l = ipvt(j) x(l) = wa1(j) 60 continue c c initialize the iteration counter. c evaluate the function at the origin, and test c for acceptance of the gauss-newton direction. c iter = 0 do 70 j = 1, n wa2(j) = diag(j)*x(j) 70 continue dxnorm = enorm(n,wa2) fp = dxnorm - delta if (fp .le. p1*delta) go to 220 c c if the jacobian is not rank deficient, the newton c step provides a lower bound, parl, for the zero of c the function. otherwise set this bound to zero. c parl = zero if (nsing .lt. n) go to 120 do 80 j = 1, n l = ipvt(j) wa1(j) = diag(l)*(wa2(l)/dxnorm) 80 continue do 110 j = 1, n sum = zero jm1 = j - 1 if (jm1 .lt. 1) go to 100 do 90 i = 1, jm1 sum = sum + r(i,j)*wa1(i) 90 continue 100 continue wa1(j) = (wa1(j) - sum)/r(j,j) 110 continue temp = enorm(n,wa1) parl = ((fp/delta)/temp)/temp 120 continue c c calculate an upper bound, paru, for the zero of the function. c do 140 j = 1, n sum = zero do 130 i = 1, j sum = sum + r(i,j)*qtb(i) 130 continue l = ipvt(j) wa1(j) = sum/diag(l) 140 continue gnorm = enorm(n,wa1) paru = gnorm/delta if (paru .eq. zero) paru = dwarf/dmin1(delta,p1) c c if the input par lies outside of the interval (parl,paru), c set par to the closer endpoint. c par = dmax1(par,parl) par = dmin1(par,paru) if (par .eq. zero) par = gnorm/dxnorm c c beginning of an iteration. c 150 continue iter = iter + 1 c c evaluate the function at the current value of par. c if (par .eq. zero) par = dmax1(dwarf,p001*paru) temp = dsqrt(par) do 160 j = 1, n wa1(j) = temp*diag(j) 160 continue call qrsolv(n,r,ldr,ipvt,wa1,qtb,x,sdiag,wa2) do 170 j = 1, n wa2(j) = diag(j)*x(j) 170 continue dxnorm = enorm(n,wa2) temp = fp fp = dxnorm - delta c c if the function is small enough, accept the current value c of par. also test for the exceptional cases where parl c is zero or the number of iterations has reached 10. c if (dabs(fp) .le. p1*delta * .or. parl .eq. zero .and. fp .le. temp * .and. temp .lt. zero .or. iter .eq. 10) go to 220 c c compute the newton correction. c do 180 j = 1, n l = ipvt(j) wa1(j) = diag(l)*(wa2(l)/dxnorm) 180 continue do 210 j = 1, n wa1(j) = wa1(j)/sdiag(j) temp = wa1(j) jp1 = j + 1 if (n .lt. jp1) go to 200 do 190 i = jp1, n wa1(i) = wa1(i) - r(i,j)*temp 190 continue 200 continue 210 continue temp = enorm(n,wa1) parc = ((fp/delta)/temp)/temp c c depending on the sign of the function, update parl or paru. c if (fp .gt. zero) parl = dmax1(parl,par) if (fp .lt. zero) paru = dmin1(paru,par) c c compute an improved estimate for par. c par = dmax1(parl,par+parc) c c end of an iteration. c go to 150 220 continue c c termination. c if (iter .eq. 0) par = zero return c c last card of subroutine lmpar. c end subroutine qrsolv(n,r,ldr,ipvt,diag,qtb,x,sdiag,wa) integer n,ldr integer ipvt(n) double precision r(ldr,n),diag(n),qtb(n),x(n),sdiag(n),wa(n) c ********** c c subroutine qrsolv c c given an m by n matrix a, an n by n diagonal matrix d, c and an m-vector b, the problem is to determine an x which c solves the system c c a*x = b , d*x = 0 , c c in the least squares sense. c c this subroutine completes the solution of the problem c if it is provided with the necessary information from the c qr factorization, with column pivoting, of a. that is, if c a*p = q*r, where p is a permutation matrix, q has orthogonal c columns, and r is an upper triangular matrix with diagonal c elements of nonincreasing magnitude, then qrsolv expects c the full upper triangle of r, the permutation matrix p, c and the first n components of (q transpose)*b. the system c a*x = b, d*x = 0, is then equivalent to c c t t c r*z = q *b , p *d*p*z = 0 , c c where x = p*z. if this system does not have full rank, c then a least squares solution is obtained. on output qrsolv c also provides an upper triangular matrix s such that c c t t t c p *(a *a + d*d)*p = s *s . c c s is computed within qrsolv and may be of separate interest. c c the subroutine statement is c c subroutine qrsolv(n,r,ldr,ipvt,diag,qtb,x,sdiag,wa) c c where c c n is a positive integer input variable set to the order of r. c c r is an n by n array. on input the full upper triangle c must contain the full upper triangle of the matrix r. c on output the full upper triangle is unaltered, and the c strict lower triangle contains the strict upper triangle c (transposed) of the upper triangular matrix s. c c ldr is a positive integer input variable not less than n c which specifies the leading dimension of the array r. c c ipvt is an integer input array of length n which defines the c permutation matrix p such that a*p = q*r. column j of p c is column ipvt(j) of the identity matrix. c c diag is an input array of length n which must contain the c diagonal elements of the matrix d. c c qtb is an input array of length n which must contain the first c n elements of the vector (q transpose)*b. c c x is an output array of length n which contains the least c squares solution of the system a*x = b, d*x = 0. c c sdiag is an output array of length n which contains the c diagonal elements of the upper triangular matrix s. c c wa is a work array of length n. c c subprograms called c c fortran-supplied ... dabs,dsqrt c c argonne national laboratory. minpack project. march 1980. c burton s. garbow, kenneth e. hillstrom, jorge j. more c c ********** integer i,j,jp1,k,kp1,l,nsing double precision cos,cotan,p5,p25,qtbpj,sin,sum,tan,temp,zero data p5,p25,zero /5.0d-1,2.5d-1,0.0d0/ c c copy r and (q transpose)*b to preserve input and initialize s. c in particular, save the diagonal elements of r in x. c do 20 j = 1, n do 10 i = j, n r(i,j) = r(j,i) 10 continue x(j) = r(j,j) wa(j) = qtb(j) 20 continue c c eliminate the diagonal matrix d using a givens rotation. c do 100 j = 1, n c c prepare the row of d to be eliminated, locating the c diagonal element using p from the qr factorization. c l = ipvt(j) if (diag(l) .eq. zero) go to 90 do 30 k = j, n sdiag(k) = zero 30 continue sdiag(j) = diag(l) c c the transformations to eliminate the row of d c modify only a single element of (q transpose)*b c beyond the first n, which is initially zero. c qtbpj = zero do 80 k = j, n c c determine a givens rotation which eliminates the c appropriate element in the current row of d. c if (sdiag(k) .eq. zero) go to 70 if (dabs(r(k,k)) .ge. dabs(sdiag(k))) go to 40 cotan = r(k,k)/sdiag(k) sin = p5/dsqrt(p25+p25*cotan**2) cos = sin*cotan go to 50 40 continue tan = sdiag(k)/r(k,k) cos = p5/dsqrt(p25+p25*tan**2) sin = cos*tan 50 continue c c compute the modified diagonal element of r and c the modified element of ((q transpose)*b,0). c r(k,k) = cos*r(k,k) + sin*sdiag(k) temp = cos*wa(k) + sin*qtbpj qtbpj = -sin*wa(k) + cos*qtbpj wa(k) = temp c c accumulate the tranformation in the row of s. c kp1 = k + 1 if (n .lt. kp1) go to 70 do 60 i = kp1, n temp = cos*r(i,k) + sin*sdiag(i) sdiag(i) = -sin*r(i,k) + cos*sdiag(i) r(i,k) = temp 60 continue 70 continue 80 continue 90 continue c c store the diagonal element of s and restore c the corresponding diagonal element of r. c sdiag(j) = r(j,j) r(j,j) = x(j) 100 continue c c solve the triangular system for z. if the system is c singular, then obtain a least squares solution. c nsing = n do 110 j = 1, n if (sdiag(j) .eq. zero .and. nsing .eq. n) nsing = j - 1 if (nsing .lt. n) wa(j) = zero 110 continue if (nsing .lt. 1) go to 150 do 140 k = 1, nsing j = nsing - k + 1 sum = zero jp1 = j + 1 if (nsing .lt. jp1) go to 130 do 120 i = jp1, nsing sum = sum + r(i,j)*wa(i) 120 continue 130 continue wa(j) = (wa(j) - sum)/sdiag(j) 140 continue 150 continue c c permute the components of z back to components of x. c do 160 j = 1, n l = ipvt(j) x(l) = wa(j) 160 continue return c c last card of subroutine qrsolv. c end subroutine qrfac(m,n,a,lda,pivot,ipvt,lipvt,rdiag,acnorm,wa, + INTGRS,REALS) integer m,n,lda,lipvt integer ipvt(lipvt) logical pivot double precision a(lda,n),rdiag(n),acnorm(n),wa(n) INTEGER INTGRS(*) DOUBLE PRECISION REALS(*) c ********** c c subroutine qrfac c c this subroutine uses householder transformations with column c pivoting (optional) to compute a qr factorization of the c m by n matrix a. that is, qrfac determines an orthogonal c matrix q, a permutation matrix p, and an upper trapezoidal c matrix r with diagonal elements of nonincreasing magnitude, c such that a*p = q*r. the householder transformation for c column k, k = 1,2,...,min(m,n), is of the form c c t c i - (1/u(k))*u*u c c where u has zeros in the first k-1 positions. the form of c this transformation and the method of pivoting first c appeared in the corresponding linpack subroutine. c c the subroutine statement is c c subroutine qrfac(m,n,a,lda,pivot,ipvt,lipvt,rdiag,acnorm,wa, c INTGRS,REALS) c c where c c m is a positive integer input variable set to the number c of rows of a. c c n is a positive integer input variable set to the number c of columns of a. c c a is an m by n array. on input a contains the matrix for c which the qr factorization is to be computed. on output c the strict upper trapezoidal part of a contains the strict c upper trapezoidal part of r, and the lower trapezoidal c part of a contains a factored form of q (the non-trivial c elements of the u vectors described above). c c lda is a positive integer input variable not less than m c which specifies the leading dimension of the array a. c c pivot is a logical input variable. if pivot is set true, c then column pivoting is enforced. if pivot is set false, c then no column pivoting is done. c c ipvt is an integer output array of length lipvt. ipvt c defines the permutation matrix p such that a*p = q*r. c column j of p is column ipvt(j) of the identity matrix. c if pivot is false, ipvt is not referenced. c c lipvt is a positive integer input variable. if pivot is false, c then lipvt may be as small as 1. if pivot is true, then c lipvt must be at least n. c c rdiag is an output array of length n which contains the c diagonal elements of r. c c acnorm is an output array of length n which contains the c norms of the corresponding columns of the input matrix a. c if this information is not needed, then acnorm can coincide c with rdiag. c c wa is a work array of length n. if pivot is false, then wa c can coincide with rdiag. c c subprograms called c c minpack-supplied ... enorm c c fortran-supplied ... dmax1,dsqrt,min0 c c argonne national laboratory. minpack project. march 1980. c burton s. garbow, kenneth e. hillstrom, jorge j. more c c ********** integer i,j,jp1,k,kmax,minmn double precision ajnorm,epsmch,one,p05,sum,temp,zero double precision enorm data one,p05,zero /1.0d0,5.0d-2,0.0d0/ c c epsmch is the machine precision. c epsmch = REALS(1) c c compute the initial column norms and initialize several arrays. c do 10 j = 1, n acnorm(j) = enorm(m,a(1,j)) rdiag(j) = acnorm(j) wa(j) = rdiag(j) if (pivot) ipvt(j) = j 10 continue c c reduce a to r with householder transformations. c minmn = min0(m,n) do 110 j = 1, minmn if (.not.pivot) go to 40 c c bring the column of largest norm into the pivot position. c kmax = j do 20 k = j, n if (rdiag(k) .gt. rdiag(kmax)) kmax = k 20 continue if (kmax .eq. j) go to 40 do 30 i = 1, m temp = a(i,j) a(i,j) = a(i,kmax) a(i,kmax) = temp 30 continue rdiag(kmax) = rdiag(j) wa(kmax) = wa(j) k = ipvt(j) ipvt(j) = ipvt(kmax) ipvt(kmax) = k 40 continue c c compute the householder transformation to reduce the c j-th column of a to a multiple of the j-th unit vector. c ajnorm = enorm(m-j+1,a(j,j)) if (ajnorm .eq. zero) go to 100 if (a(j,j) .lt. zero) ajnorm = -ajnorm do 50 i = j, m a(i,j) = a(i,j)/ajnorm 50 continue a(j,j) = a(j,j) + one c c apply the transformation to the remaining columns c and update the norms. c jp1 = j + 1 if (n .lt. jp1) go to 100 do 90 k = jp1, n sum = zero do 60 i = j, m sum = sum + a(i,j)*a(i,k) 60 continue temp = sum/a(j,j) do 70 i = j, m a(i,k) = a(i,k) - temp*a(i,j) 70 continue if (.not.pivot .or. rdiag(k) .eq. zero) go to 80 temp = a(j,k)/rdiag(k) rdiag(k) = rdiag(k)*dsqrt(dmax1(zero,one-temp**2)) if (p05*(rdiag(k)/wa(k))**2 .gt. epsmch) go to 80 rdiag(k) = enorm(m-j,a(jp1,k)) wa(k) = rdiag(k) 80 continue 90 continue 100 continue rdiag(j) = -ajnorm 110 continue return c c last card of subroutine qrfac. c end double precision function enorm(n,x) integer n double precision x(n) c ********** c c function enorm c c given an n-vector x, this function calculates the c euclidean norm of x. c c the euclidean norm is computed by accumulating the sum of c squares in three different sums. the sums of squares for the c small and large components are scaled so that no overflows c occur. non-destructive underflows are permitted. underflows c and overflows do not occur in the computation of the unscaled c sum of squares for the intermediate components. c the definitions of small, intermediate and large components c depend on two constants, rdwarf and rgiant. the main c restrictions on these constants are that rdwarf**2 not c underflow and rgiant**2 not overflow. the constants c given here are suitable for every known computer. c c the function statement is c c double precision function enorm(n,x) c c where c c n is a positive integer input variable. c c x is an input array of length n. c c subprograms called c c fortran-supplied ... dabs,dsqrt c c argonne national laboratory. minpack project. march 1980. c burton s. garbow, kenneth e. hillstrom, jorge j. more c c ********** integer i double precision agiant,floatn,one,rdwarf,rgiant,s1,s2,s3,xabs, * x1max,x3max,zero data one,zero,rdwarf,rgiant /1.0d0,0.0d0,3.834d-20,1.304d19/ s1 = zero s2 = zero s3 = zero x1max = zero x3max = zero floatn = n agiant = rgiant/floatn do 90 i = 1, n xabs = dabs(x(i)) if (xabs .gt. rdwarf .and. xabs .lt. agiant) go to 70 if (xabs .le. rdwarf) go to 30 c c sum for large components. c if (xabs .le. x1max) go to 10 s1 = one + s1*(x1max/xabs)**2 x1max = xabs go to 20 10 continue s1 = s1 + (xabs/x1max)**2 20 continue go to 60 30 continue c c sum for small components. c if (xabs .le. x3max) go to 40 s3 = one + s3*(x3max/xabs)**2 x3max = xabs go to 50 40 continue if (xabs .ne. zero) s3 = s3 + (xabs/x3max)**2 50 continue 60 continue go to 80 70 continue c c sum for intermediate components. c s2 = s2 + xabs**2 80 continue 90 continue c c calculation of norm. c if (s1 .eq. zero) go to 100 enorm = x1max*dsqrt(s1+(s2/x1max)/x1max) go to 130 100 continue if (s2 .eq. zero) go to 110 if (s2 .ge. x3max) * enorm = dsqrt(s2*(one+(x3max/s2)*(x3max*s3))) if (s2 .lt. x3max) * enorm = dsqrt(x3max*((s2/x3max)+(x3max*s3))) go to 120 110 continue enorm = x3max*dsqrt(s3) 120 continue 130 continue return c c last card of function enorm. c end subroutine fdjac2(fcn,m,n,x,fvec,fjac,ldfjac,iflag,epsfcn,wa, + INTGRS,REALS) integer m,n,ldfjac,iflag double precision epsfcn double precision x(n),fvec(m),fjac(ldfjac,n),wa(m) INTEGER INTGRS(*) DOUBLE PRECISION REALS(*) c ********** c c subroutine fdjac2 c c this subroutine computes a forward-difference approximation c to the m by n jacobian matrix associated with a specified c problem of m functions in n variables. c c the subroutine statement is c c subroutine fdjac2(fcn,m,n,x,fvec,fjac,ldfjac,iflag,epsfcn,wa, c INTGRS,REALS) c c where c c fcn is the name of the user-supplied subroutine which c calculates the functions. fcn must be declared c in an external statement in the user calling c program, and should be written as follows. c c subroutine fcn(m,n,x,fvec,iflag,INTGRS,REALS) c integer m,n,iflag c double precision x(n),fvec(m) c ---------- c calculate the functions at x and c return this vector in fvec. c ---------- c return c end c c the value of iflag should not be changed by fcn unless c the user wants to terminate execution of fdjac2. c in this case set iflag to a negative integer. c c m is a positive integer input variable set to the number c of functions. c c n is a positive integer input variable set to the number c of variables. n must not exceed m. c c x is an input array of length n. c c fvec is an input array of length m which must contain the c functions evaluated at x. c c fjac is an output m by n array which contains the c approximation to the jacobian matrix evaluated at x. c c ldfjac is a positive integer input variable not less than m c which specifies the leading dimension of the array fjac. c c iflag is an integer variable which can be used to terminate c the execution of fdjac2. see description of fcn. c c epsfcn is an input variable used in determining a suitable c step length for the forward-difference approximation. this c approximation assumes that the relative errors in the c functions are of the order of epsfcn. if epsfcn is less c than the machine precision, it is assumed that the relative c errors in the functions are of the order of the machine c precision. c c wa is a work array of length m. c c subprograms called c c user-supplied ...... fcn c c fortran-supplied ... dabs,dmax1,dsqrt c c argonne national laboratory. minpack project. march 1980. c burton s. garbow, kenneth e. hillstrom, jorge j. more c c ********** integer i,j double precision eps,epsmch,h,temp,zero external fcn data zero /0.0d0/ c c epsmch is the machine precision. c epsmch = REALS(1) c eps = dsqrt(dmax1(epsfcn,epsmch)) do 20 j = 1, n temp = x(j) h = eps*dabs(temp) if (h .eq. zero) h = eps x(j) = temp + h call fcn(m,n,x,wa,iflag,INTGRS,REALS) if (iflag .lt. 0) go to 30 x(j) = temp do 10 i = 1, m fjac(i,j) = (wa(i) - fvec(i))/h 10 continue 20 continue 30 continue return c c last card of subroutine fdjac2. c end