TOMS703, algorithm to solve systems of stiff ODEs. ! ALGORITHM 703 , COLLECTED ALGORITHMS FROM ACM. ! THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, ! VOL. 18, NO. 2, JUNE, 1992, PP. 156-158. ! This ELF90-compatible version by Alan Miller ! Alan.Miller @ vic.cmis.csiro.au ! http://www.ozemail.com.au/~milleraj ! Latest revision - 10 August 1998 MODULE common703 IMPLICIT NONE ! Contains the COMMON declarations from the Fortran 77 version, ! and the precision declaration. INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(14, 60) ! COMMON BLOCKS OF INTEREST TO THE USER. ! /MEBDF2/ ! HUSED LAST STEPSIZE SUCCESSFULLY USED BY THE INTEGRATOR ! NQUSED ORDER LAST SUCCESSFULLY USED ! NSTEP NUMBER OF STEPS TAKEN SO FAR ! NFE NUMBER OF FUNCTION EVALUATIONS USED SO FAR ! NJE NUMBER OF JACOBIAN EVALUATIONS USED SO FAR ! NDEC NUMBER OF LU DECOMPOSITIONS USED SO FAR ! NBSOL NUMBER OF 'BACKSOLVES' USED SO FAR ! NPSET NUMBER OF TIMES A NEW COEFFICIENT MATRIX HAS BEEN FORMED SO FAR ! NCOSET NUMBER OF TIMES THE ORDER OF THE METHOD USED HAS BEEN ! CHANGED SO FAR ! MAXORD THE MAXIMUM ORDER USED SO FAR IN THE INTEGRATION. ! PUT THE FOLLOWING STATEMENT IN THE MAIN PROGRAM: ! USE common703 ! BLOCK DATA ! MACHINE DEPENDENT CONSTANTS ARE SET HERE ! UROUND THE SMALLEST +VE NUMBER SUCH THAT 1.0 + UROUND > 1.0 ! EPSJAC THIS IS USED IN FORMING A NUMERICAL JACOBIAN AND IS ! EPSJAC = SQRT( UROUND ) ! COMMON /machne/uround,epsjac REAL (dp), PARAMETER :: uround = EPSILON(1.0_dp) REAL (dp), SAVE :: epsjac ! COMMON /fails/nfail1,nfail2,nt ! COMMON /mebdf1/t,h,hmin,hmax,epscpy,ncopy,mfcopy,kflag,jstart,toutcp ! COMMON /mebdf2/hused,nqused,nstep,nfe,nje,ndec,nbsol,npset,ncoset,maxord ! COMMON /mebdf3/iweval,newpar,nrenew,jsinup,jsnold,idoub,jchang,l,nq ! COMMON /mebdf4/m1,m2,meqc1,meqc2,mq1tmp,mq2tmp,isamp ! COMMON /mebdf5/qi,rc,rh,rmax,told ! COMMON /mebdf6/crate1,crate2,tcrat1,tcrat2,avnew2,avold2,avnewj,avoldj ! COMMON /mebdfl/cfail,jnewim,sample ! COMMON /stpsze/hstpsz,top INTEGER, SAVE :: idoub,isamp,iweval,jchang,jsinup,jsnold,jstart,kflag,l, & m1,m2,maxord,meqc1,meqc2,mfcopy,mq1tmp,mq2tmp,nbsol, & ncopy,ncoset,ndec,newpar,nfail1,nfail2,nfe,nje,npset,nq, & nqused,nrenew,nstep,nt LOGICAL, SAVE :: cfail,jnewim,sample REAL (dp), SAVE :: avnew2,avold2,avnewj,avoldj,epscpy,crate1,crate2,h,hmin, & hmax,hstpsz(2,14),hused,qi,rc,rh,rmax,t,tcrat1,tcrat2, & told,top,toutcp END MODULE common703 MODULE toms703 USE common703 IMPLICIT NONE INTERFACE SUBROUTINE pderv(t,y,pw) USE common703, ONLY: dp IMPLICIT NONE REAL (dp), INTENT(IN) :: t REAL (dp), INTENT(IN) :: y(:,:) REAL (dp), INTENT(OUT) :: pw(:,:) ! pw(4,4) END SUBROUTINE pderv END INTERFACE CONTAINS SUBROUTINE mebdf(n, t0, h0, y0, tout, tend, eps, mf, INDEX, lout, fcn) IMPLICIT NONE INTEGER, INTENT(IN) :: n REAL (dp), INTENT(IN OUT) :: t0 REAL (dp), INTENT(IN OUT) :: h0 REAL (dp), INTENT(IN OUT) :: y0(:) REAL (dp), INTENT(IN OUT) :: tout REAL (dp), INTENT(IN OUT) :: tend REAL (dp), INTENT(IN OUT) :: eps INTEGER, INTENT(IN) :: mf INTEGER, INTENT(OUT) :: INDEX INTEGER, INTENT(IN) :: lout INTERFACE SUBROUTINE fcn(t,y,ydot) USE common703, ONLY: dp IMPLICIT NONE REAL (dp), INTENT(IN) :: t REAL (dp), INTENT(IN) :: y(:) REAL (dp), INTENT(OUT) :: ydot(:) END SUBROUTINE fcn END INTERFACE ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< ! THIS IS THE JULY 1990 VERSION OF OVDRIV, A PACKAGE FOR ! THE SOLUTION OF THE INITIAL VALUE PROBLEM FOR SYSTEMS OF ! ORDINARY DIFFERENTIAL EQUATIONS ! DY/DT = F(Y,T), Y=(Y(1),Y(2),Y(3), . . . ,Y(N)) ! SUBROUTINE OVDRIV IS A DRIVER ROUTINE FOR THIS PACKAGE ! REFERENCES ! 1. J.R. CASH, ON THE INTEGRATION OF STIFF SYSTEMS OF O.D.E.S ! USING EXTENDED BACKWARD DIFFERENTIATION FORMULAE, ! NUMER. MATH. 34, 235-246 (1980) ! 2. THE INTEGRATION OF STIFF INITIAL VALUE PROBLEMS IN O.D.E.S ! USING MODIFIED EXTENDED BACKWARD DIFFERENTIATION FORMULAE, ! COMP. AND MATHS. WITH APPLICS., 9, 645-657, (1983). ! 3. J.R. CASH, STABLE RECURSIONS WITH APPLICATIONS TO THE ! NUMERICAL SOLUTION OF STIFF SYSTEMS, ACADEMIC PRESS.(1979) ! 4. A.C. HINDMARSH, GEAR.. ORDINARY DIFFERENTIAL EQUATION ! SYSTEM SOLVER, UCID-30001 REV. 3, LAWRENCE LIVERMORE ! LABORATORY, P.O. BOX 808, LIVERMORE, CA 94550, DEC. 1974. ! 5. J.R. CASH AND S. CONSIDINE, AN MEBDF CODE FOR STIFF INITIAL ! VALUE PROBLEMS, ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE, 1992. ! ---------------------------------------------------------------- ! OVDRIV IS TO BE CALLED ONCE FOR EACH OUTPUT VALUE OF T, AND ! IN TURN MAKES REPEATED CALLS TO THE CORE INTEGRATOR STIFF. ! THE INPUT PARAMETERS ARE .. ! N = THE NUMBER OF FIRST ORDER DIFFERENTIAL EQUATIONS. ! T0 = THE INITIAL VALUE OF T, THE INDEPENDENT VARIABLE ! (USED ONLY ON THE FIRST CALL) ! H0 = THE NEXT STEP SIZE IN T (USED FOR INPUT ONLY ON THE FIRST CALL) ! Y0 = A VECTOR OF LENGTH N CONTAINING THE INITIAL VALUES OF Y ! (USED FOR INPUT ONLY ON FIRST CALL) ! TEND = END OF THE RANGE OF INTEGRATION. ! TOUT = THE VALUE OF T AT WHICH OUTPUT IS DESIRED NEXT. ! INTEGRATION WILL NORMALLY GO SLIGHTLY BEYOND TOUT AND THE ! PACKAGE WILL INTERPOLATE TO T = TOUT MUST HAVE TEND >= TOUT IF ! H0 IS POSITIVE AND TOUT >= TEND OTHERWISE. ! EPS = THE RELATIVE ERROR BOUND (USED ONLY ON THE FIRST CALL, ! UNLESS INDEX = -1). SINGLE STEP ERROR ESTIMATES DIVIDED BY ! YMAX(I) WILL BE KEPT LESS THAN EPS IN ROOT-MEAN-SQUARE NORM ! (I.E. EUCLIDEAN NORM DIVIDED BY SQRT(N) ). ! THE VECTOR YMAX OF WEIGHTS IS COMPUTED IN OVDRIV. ! INITIALLY YMAX(I) IS ABS(Y(I)) WITH A DEFAULT VALUE OF 1 IF ! Y(I)=0 INITIALLY. ! THEREAFTER, YMAX(I) IS THE LARGEST VALUE OF ABS(Y(I)) ! SEEN SO FAR, OR THE INITIAL VALUE YMAX(I) IF THAT IS LARGER ! MF = THE METHOD FLAG. AT PRESENT ONLY MF=21 OR 22 IS ALLOWED ! WHICH IS EXTENDED BACKWARD DIFFERENTIATION FORMULAE USING THE ! CHORD METHOD WITH ANALYTIC OR NUMERICAL JACOBIAN. ! THE USER NEEDS TO SPECIFY SUBROUTINE PDERV IF MF=21 (SEE BELOW) ! INDEX = INTEGER USED ON INPUT TO INDICATE TYPE OF CALL, ! 1 THIS IS THE FIRST CALL FOR THE PROBLEM ! INDEX HAS TO BE SET = 0 AFTER FIRST CALL. ! 0 THIS IS NOT THE FIRST CALL FOR THIS PROBLEM ! AND INTEGRATION IS TO CONTINUE UP TO T=TOUT OR BEYOND. ! -1 THIS IS NOT THE FIRST CALL FOR THE PROBLEM, ! AND THE USER HAS RESET N, EPS, AND/OR MF. ! 2 SAME AS 0 EXCEPT THAT TOUT HAS TO BE HIT ! EXACTLY (NO INTERPOLATION IS DONE). ! ASSUMES TOUT >= THE CURRENT T. ! 3 ONE STEP MODE. SAME AS 0 EXCEPT CONTROL ! RETURNS TO CALLING PROGRAM AFTER ONE STEP. ! SINCE THE NORMAL OUTPUT VALUE OF INDEX IS 0, ! IT NEED NOT BE RESET FOR NORMAL CONTINUATION. ! LOUT = LOGICAL OUTPUT DEVICE. ! AFTER THE INITIAL CALL, IF A NORMAL RETURN OCCURED AND A NORMAL ! CONTINUATION IS DESIRED, SIMPLY RESET TOUT AND CALL AGAIN. ! ALL OTHER PARAMETERS WILL BE READY FOR THE NEXT CALL. ! A CHANGE OF PARAMETERS WITH INDEX = -1 CAN BE MADE AFTER ! EITHER A SUCCESSFUL OR AN UNSUCCESSFUL RETURN. ! THE OUTPUT PARAMETERS ARE.. ! T0 = THE VALUE OF T WHICH RELATES TO THE CURRENT SOLUTION POINT Y0() ! H0 = THE STEPSIZE H USED LAST, WHETHER SUCCESSFULLY OR NOT. ! Y0 = THE COMPUTED VALUES OF Y AT T = TOUT ! TOUT = UNCHANGED FROM ITS INPUT VALUE. ! INDEX = INTEGER USED ON OUTPUT TO INDICATE RESULTS, WITH ! THE FOLLOWING VALUES AND MEANINGS.. ! 3 ONE STEP MODE IS BEING USED AND A SUCCESSFUL STEP HAS BEEN ! COMPLETED. THE USER SHOULD CALL MEBDF AGAIN WITH INDEX = 3. ! 1 ONE SUCCESSFUL STEP HAS BEEN COMPLETED AND THE ONE STEP MODE ! IS NOT BEING USED. THE USER MUST NOW SET INDEX = 0 AND CALL ! MEBDF AGAIN. ! 0 INTEGRATION WAS COMPLETED TO TOUT OR BEYOND. ! -1 THE INTEGRATION WAS HALTED AFTER FAILING TO PASS THE ERROR TEST ! EVEN AFTER REDUCING H BY A FACTOR OF 1.E10 FROM ITS INITIAL VALUE. ! -2 AFTER SOME INITIAL SUCCESS, THE INTEGRATION WAS HALTED EITHER BY ! REPEATED ERROR TEST FAILURES OR BY A TEST ON EPS. ! TOO MUCH ACCURACY HAS BEEN REQUESTED. ! -3 THE INTEGRATION WAS HALTED AFTER FAILING TO ACHIEVE CORRECTOR ! CONVERGENCE EVEN AFTER REDUCING H BY A FACTOR OF 1.E10 FROM ITS ! INITIAL VALUE. ! -4 IMMEDIATE HALT BECAUSE OF ILLEGAL VALUES OF INPUT PARAMETERS. ! SEE PRINTED MESSAGE. ! -5 INDEX WAS -1 ON INPUT, BUT THE DESIRED CHANGES OF PARAMETERS WERE ! NOT IMPLEMENTED BECAUSE TOUT WAS NOT BEYOND T. ! INTERPOLATION AT T = TOUT WAS PERFORMED AS ON A NORMAL RETURN. ! TO TRY AGAIN, SIMPLY CALL AGAIN WITH INDEX = -1 AND A NEW TOUT. ! IN ADDITION TO OVDRIV, THE FOLLOWING ROUTINES ARE PROVIDED IN THE PACKAGE.. ! INTERP( - ) INTERPOLATES TO GET THE OUTPUT VALUES ! AT T=TOUT FROM THE DATA IN THE Y ARRAY. ! STIFF( - ) IS THE CORE INTEGRATOR ROUTINE. IT PERFORMS A ! SINGLE STEP AND ASSOCIATED ERROR CONTROL. ! COSET( - ) SETS COEFFICIENTS FOR BACKWARD DIFFERENTIATION ! SCHEMES FOR USE IN THE CORE INTEGRATOR. ! PSET( - ) COMPUTES AND PROCESSES THE JACOBIAN MATRIX J = DF/DY ! DEC( - ) PERFORMS AN LU DECOMPOSITION ON A MATRIX. ! SOL( - ) SOLVES LINEAR SYSTEMS A*X = B AFTER DEC ! HAS BEEN CALLED FOR THE MATRIX A ! THE FOLLOWING ROUTINES ARE TO BE SUPPLIED BY THE USER.. ! FCN(T,Y,YDOT) COMPUTES THE FUNCTION YDOT = F(Y,T), THE ! RIGHT HAND SIDE OF THE O.D.E. ! HERE Y AND YDOT ARE VECTORS OF LENGTH N ! PDERV(T,Y,PD) COMPUTES THE N*N JACOBIAN MATRIX OF PARTIAL DERIVATIVES ! AND STORES IT IN PD AS AN N BY N ARRAY. ! PD(I,J) IS TO BE SET TO THE PARTIAL DERIVATIVE OF ! YDOT(I) WITH RESPECT TO Y(J). ! PDERV IS CALLED ONLY IF MITER = 1. ! OTHERWISE A DUMMY ROUTINE CAN BE SUBSTITUTED. ! THE DIMENSIONS OF PW BELOW MUST BE AT LEAST N x N. THE DIMENSIONS OF ! YMAX, ERROR, SAVE1, SAVE2, IPIV AND THE FIRST DIMENSION OF Y SHOULD ALL ! BE AT LEAST N. ! ***************************************************************** epsjac = SQRT(uround) IF (INDEX == 1) THEN IF (n <= 0) THEN WRITE (lout,9020) n INDEX = -4 END IF IF (INDEX < 0) RETURN END IF CALL ovdriv(n, t0, h0, y0, tout, tend, eps, mf, INDEX, lout, fcn) RETURN 9020 FORMAT (/,/,' ***** ERROR ***** INTEGRATION HALTED IN MEBDF',/,/, & ' >>> ILLEGAL VALUE FOR NUMBER OF EQUATIONS <<< ',/, & ' WITH N = ',i6) END SUBROUTINE mebdf SUBROUTINE ovdriv(n, t0, h0, y0, tout, tend, eps, mf, INDEX, lout, fcn) IMPLICIT NONE INTEGER, INTENT(IN) :: n REAL (dp), INTENT(OUT) :: t0 REAL (dp), INTENT(OUT) :: h0 REAL (dp), INTENT(IN OUT) :: y0(:) REAL (dp), INTENT(IN) :: tout REAL (dp), INTENT(IN OUT) :: tend REAL (dp), INTENT(IN) :: eps INTEGER, INTENT(IN) :: mf INTEGER, INTENT(OUT) :: INDEX INTEGER, INTENT(IN) :: lout INTERFACE SUBROUTINE fcn(t,y,ydot) USE common703, ONLY: dp IMPLICIT NONE REAL (dp), INTENT(IN) :: t REAL (dp), INTENT(IN) :: y(:) REAL (dp), INTENT(OUT) :: ydot(:) END SUBROUTINE fcn END INTERFACE ! START OF PROGRAM PROPER ! .. ! .. LOCAL SCALARS .. INTEGER :: i,kgo,nhcut REAL (dp) :: ayi, d ! .. ! Local arrays REAL (dp) :: y(n,15), yhold(n,15), ymax(n), pw(n,n), pwcopy(n*n) INTEGER :: ipiv(n) ! .. IF (INDEX == 0) THEN ! I.E. NORMAL CONTINUATION OF INTEGRATION t0=t hmax = ABS(tend-t0)*10.0D+0 IF ((t-tout)*h >= 0.0D+0) THEN ! HAVE OVERSHOT THE ENDPOINT, SO INTERPOLATE CALL interp(n,jstart,h,t,y,tout,y0) INDEX = kflag t0 = tout h0 = h RETURN END IF ELSE IF (INDEX == 2) THEN ! I.E. CONTINUING INTEGRATION BUT WISH TO HIT TOUT t0 = t hmax = ABS(tend-t0)*10.0D+0 IF (((t+h)-tout)*h > 0.0D+0) THEN ! WE HAVE ALREADY OVERSHOT THE END OR WE WILL ! DO SO ON THE NEXT STEP IF (((t-tout)*h >= 0.0D+0) .OR. (ABS(t-tout) <= & 100.0D+0*uround*hmax)) THEN ! HAVE OVERSHOT, SO INTERPOLATE CALL interp(n,jstart,h,t,y,tout,y0) t0 = tout h0 = h INDEX = kflag RETURN ELSE ! WILL PASS TOUT ON NEXT STEP WITH CURRENT STEPSIZE ! SO REDUCE STEPSIZE TO HIT TOUT 'EXACTLY' h = (tout-t)* (1.0D+0-4.0D+0*uround) jstart = -1 END IF END IF ELSE IF (INDEX == -1) THEN ! NOT FIRST CALL BUT PARAMETERS RESET t0 = t IF ((t-tout)*h >= 0.0D+0) THEN ! HAVE OVERSHOT TOUT WRITE (lout,9080) t,tout,h CALL interp(n,jstart,h,t,y,tout,y0) h0 = h t0 = tout INDEX = -5 RETURN ELSE jstart = -1 ncopy = n epscpy = eps mfcopy = mf END IF ELSE IF (INDEX == 3) THEN t0 = t IF ((t-tout)*h >= 0.0D+0) THEN ! HAVE OVERSHOT,SO INTERPOLATE CALL interp(n,jstart,h,t,y,tout,y0) INDEX = kflag t0 = tout h0 = h RETURN END IF ELSE ! INDEX SHOULD BE 1 AND THIS IS THE FIRST CALL FOR THIS PROBLEM ! CHECK THE ARGUMENTS THAT WERE PASSED FOR CORRECTNESS IF (INDEX /= 1) THEN ! VALUE OF INDEX NOT ALLOWED WRITE (lout,9070) INDEX INDEX = -4 END IF IF (eps <= 0.0D+0) THEN ! ILLEGAL VALUE FOR RELATIVE ERROR TOLERENCE WRITE (lout,9040) INDEX = -4 END IF IF (n <= 0) THEN ! ILLEGAL VALUE FOR THE NUMBER OF EQUATIONS WRITE (lout,9050) INDEX = -4 END IF IF ((t0-tout)*h0 >= 0.0D+0) THEN ! PARAMETERS FOR INTEGRATION ARE ILLEGAL WRITE (lout,9060) INDEX = -4 END IF IF((tend-tout)*h0 < -100.0D+0*uround) THEN ! TOUT IS BEYOND TEND. WRITE(lout,9110) tout,tend INDEX = -4 END IF IF ((mf /= 21) .AND. (mf /= 22)) THEN ! ILLEGAL VALUE FOR METHOD FLAG WRITE (lout,9090) mf INDEX = -4 END IF IF (INDEX /= 1) THEN RETURN ELSE ! THE INITIAL PARAMETERS ARE O.K. SO INITIALISE EVERYTHING ! ELSE NECESSARY FOR THE INTEGRATION. ! IF VALUES OF YMAX OTHER THAN THOSE SET BELOW ARE DESIRED, ! THEY SHOULD BE SET HERE. ALL YMAX(I) MUST BE POSITIVE. IF ! VALUES FOR HMIN OR HMAX, THE BOUNDS ON ABS(H), OTHER THAN ! THOSE BELOW ARE DESIRED, THEY SHOULD BE SET BELOW DO i = 1,n ymax(i) = ABS(y0(i)) IF (ymax(i) == 0.0D+0) ymax(i) = 1.0D+0 y(i,1) = y0(i) END DO ncopy = n t = t0 h = h0 hmin = ABS(h0) hmax = ABS(t0-tend)*10.0D+0 epscpy = eps mfcopy = mf jstart = 0 nhcut = 0 END IF END IF ! <<<<<<<<<<<<<<<<< ! < TAKE A STEP > ! <<<<<<<<<<<<<<<<< 20 IF ((t+h) == t) THEN WRITE (lout,9000) nt = nt + 1 END IF CALL stiff(eps, h, hmax, hmin, jstart, kflag, mf, t, tend, y, n, ymax, & pw, pwcopy, yhold, ipiv, lout, fcn) kgo = 1 - kflag IF (kgo == 1) THEN ! NORMAL RETURN FROM STIFF GO TO 30 ELSE IF (kgo == 2) THEN ! COULD NOT ACHIEVE ERROR WITH HMIN ! SO CHOP HMIN IF WE HAVEN'T DONE SO 10 TIMES GO TO 60 ELSE IF (kgo == 3) THEN ! ERROR SMALLER THAN CAN BE HANDLED FOR THIS PROBLEM WRITE (lout, 9010) t, h GO TO 70 ELSE IF (kgo == 4) THEN ! COULD NOT ACHIEVE CONVERGENCE WITH HMIN WRITE (lout, 9030) t GO TO 60 END IF ! --------------------------------------------------------------------- ! NORMAL RETURN FROM THE INTEGRATOR. ! THE WEIGHTS YMAX(I) ARE UPDATED. IF DIFFERENT VALUES ARE DESIRED, ! THEY SHOULD BE SET HERE. A TEST IS MADE FOR EPS BEING TOO SMALL ! FOR THE MACHINE PRECISION. ! ANY OTHER TESTS OR CALCULATIONS THAT ARE REQUIRED AFTER EVERY ! STEP SHOULD BE INSERTED HERE. ! IF INDEX = 3, Y0 IS SET TO THE CURRENT Y VALUES ON RETURN. ! IF INDEX = 2, H IS CONTROLLED TO HIT TOUT (WITHIN ROUNDOFF ! ERROR), AND THEN THE CURRENT Y VALUES ARE PUT IN Y0 ON RETURN. ! FOR ANY OTHER VALUE OF INDEX, CONTROL RETURNS TO THE INTEGRATOR ! UNLESS TOUT HAS BEEN REACHED. THEN INTERPOLATED VALUES OF Y ARE ! COMPUTED AND STORED IN Y0 ON RETURN. ! IF INTERPOLATION IS NOT DESIRED, THE CALL TO INTERP SHOULD BE ! REMOVED AND CONTROL TRANSFERRED TO STATEMENT 500 INSTEAD OF 520. ! -------------------------------------------------------------------- 30 d = 0.0D+0 DO i = 1,n ayi = ABS(y(i,1)) ymax(i) = MAX(ymax(i),ayi) d = d + (ayi/ymax(i))**2 END DO d = d* (uround/eps)**2 IF (d > DBLE(n)) THEN ! EPS SMALLER THAN MACHINE PRECISION WRITE (lout,9020) t kflag = -2 GO TO 70 END IF IF (INDEX == 3.OR.INDEX == 1) GO TO 70 IF (ABS(t-tout) <= ABS(15.0D+0*uround*tout)) THEN ! EFFECTIVELY WE HAVE HIT TOUT INDEX = kflag t0 = tout DO i = 1,n y0(i) = y(i,1) END DO h0 = h RETURN END IF IF (INDEX == 2) THEN ! CONTINUING INTEGRATION BUT MUST HIT TOUT EXACTLY IF (((t+h)-tout)*h > 0.0D+0) THEN ! WE HAVE ALREADY OVERSHOT THE END OR WE WILL DO ! SO ON THE NEXT STEP IF (((t-tout)*h >= 0.0D+0) .OR. (ABS(t-tout) <= & 100.0D+0*uround*hmax)) THEN ! HAVE OVERSHOT, SO INTERPOLATE CALL interp(n,jstart,h,t,y,tout,y0) t0 = tout h0 = h INDEX = kflag RETURN ELSE ! WILL PASS TOUT ON NEXT STEP WITH CURRENT STEPSIZE ! SO REDUCE STEPSIZE TO HIT TOUT 'EXACTLY' h = (tout-t)* (1.0D+0-4.0D+0*uround) jstart = -1 END IF END IF ELSE IF ((t-tout)*h >= 0.0D+0) THEN ! HAVE OVERSHOT, SO INTERPOLATE CALL interp(n,jstart,h,t,y,tout,y0) INDEX = kflag h0 = h t0 = tout RETURN END IF GO TO 20 ! ------------------------------------------------------------------- ! ON AN ERROR RETURN FROM THE INTEGRATOR, AN IMMEDIATE RETURN OCCURS ! IF KFLAG = -2, AND RECOVERY ATTEMPTS ARE MADE OTHERWISE. ! H AND HMIN ARE REDUCED BY A FACTOR OF .1 UP TO 10 TIMES ! BEFORE GIVING UP. ! -------------------------------------------------------------------- 60 IF (nhcut == 10) THEN ! HAVE REDUCED H TEN TIMES WRITE (lout,9100) GO TO 70 END IF nhcut = nhcut + 1 hmin = 0.1D+0*hmin h = 0.1D+0*h IF (nstep > 0) THEN jstart = -1 ELSE ! THE INITIAL STEPSIZE WAS TOO LARGE. RESET CERTAIN VARIABLES ! AND TRY AGAIN. nje = 0 nfe = 1 cfail = .true. newpar = 0 mq1tmp = 0 mq2tmp = 0 meqc1 = 0 meqc2 = 0 tcrat1 = 0.0D+0 tcrat2 = 0.0D+0 crate1 = 1.0D+0 crate2 = 1.0D+0 nt = 0 nstep = 0 nbsol = 0 npset = 0 ncoset = 0 ndec = 0 jstart = -1 END IF GO TO 20 70 IF(ABS(t-tout) > 1000.0D+0*uround) THEN DO i = 1,n y0(i) = y(i,1) END DO t0 = t ELSE ! HAVE EITHER PASSED TOUT OR WE ARE EXTREMELY CLOSE TO IT ! SO INTERPOLATE. CALL interp(n,jstart,h,t,y,tout,y0) t0 = tout INDEX = kflag END IF h0 = h IF(kflag /= 0) INDEX = kflag RETURN ! -------------------------- END OF SUBROUTINE OVDRIV ----------------- 9000 FORMAT (' WARNING.. T + H = T ON NEXT STEP.') 9010 FORMAT (/,/,' KFLAG = -2 FROM INTEGRATOR AT T = ',e16.8,' H =', & e16.8,/, ' THE REQUESTED ERROR IS SMALLER THAN CAN BE HANDLED',/,/) 9020 FORMAT (/,/,' INTEGRATION HALTED BY MEBDF AT T = ',e16.8,/, & ' EPS TOO SMALL TO BE ATTAINED FOR THE MACHINE PRECISION', /) 9030 FORMAT (/,/,' KFLAG = -3 FROM INTEGRATOR AT T = ',e16.8,/, & ' CORRECTOR CONVERGENCE COULD NOT BE ACHIEVED',/) 9040 FORMAT (/,/,' ILLEGAL INPUT.. EPS .LE. 0.',/,/) 9050 FORMAT (/,/,' ILLEGAL INPUT.. N .LE. 0',/,/) 9060 FORMAT (/,/,' ILLEGAL INPUT.. (T0-TOUT)*H >= 0.',/,/) 9070 FORMAT (/,/,' ILLEGAL INPUT.. INDEX =',i5,/,/) 9080 FORMAT (/,/,' INDEX = -1 ON INPUT WITH (T-TOUT)*H >= 0.',/, & ' T =',e16.8,' TOUT =',e16.8,' H =',e16.8,/, & ' INTERPOLATION WAS DONE AS ON NORMAL RETURN.',/, & ' DESIRED PARAMETER CHANGES WERE NOT MADE.') 9090 FORMAT (/,/,' ILLEGAL INPUT.. METHOD FLAG, MF, = ',i6,/, & ' ALLOWED VALUES ARE 21 OR 22',/) 9100 FORMAT (/,/,' PROBLEM APPEARS UNSOLVABLE WITH GIVEN INPUT',/, & ' HMIN REDUCED BY A FACTOR OF 1.0E10',/,/) 9110 FORMAT(/,/,' TOUT IS BEYOND TEND' ,/, & 'TOUT =' ,e16.8, ' TEND = ', e16.8,/, & 'RESET TOUT,TEND SO THAT (TEND-TOUT)*H0 >=0' ) END SUBROUTINE ovdriv SUBROUTINE interp(n,jstart,h,t,y,tout,y0) IMPLICIT NONE INTEGER, INTENT(IN) :: n INTEGER, INTENT(IN) :: jstart REAL (dp), INTENT(IN) :: h REAL (dp), INTENT(IN) :: t REAL (dp), INTENT(IN) :: y(:,:) ! y(n,15) REAL (dp), INTENT(IN) :: tout REAL (dp), INTENT(OUT) :: y0(n) ! .. ! .. LOCAL SCALARS .. INTEGER :: i,j,l REAL (dp) :: s, s1 ! .. DO i = 1,n y0(i) = y(i,1) END DO l = jstart + 2 s = (tout-t)/h s1 = 1.0D+0 DO j = 2,l s1 = s1* (s + DBLE(j-2))/DBLE(j-1) DO i = 1,n y0(i) = y0(i) + s1*y(i,j) END DO END DO RETURN ! -------------- END OF SUBROUTINE INTERP --------------------------- END SUBROUTINE interp SUBROUTINE coset(nq,el,elst,tq,maxder) IMPLICIT NONE INTEGER, INTENT(IN) :: nq REAL (dp), INTENT(OUT) :: el(:) REAL (dp), INTENT(IN OUT) :: elst(:) REAL (dp), INTENT(OUT) :: tq(:) INTEGER, INTENT(OUT) :: maxder ! ------------------------------------------------------------------------- ! COSET IS CALLED BY THE INTEGRATOR AND SETS THE COEFFICIENTS USED BY THE ! CONVENTIONAL BACKWARD DIFFERENTIATION SCHEME. THE VECTOR EL OF LENGTH ! NQ+1 DETERMINES THE BASIC METHOD. THE VECTOR ELST OF DIMENSION NQ+2 ! DETERMINES THE MODIFIED EXTENDED BACKWARD DIFFERENTIATION FORMULAE ! COEFFICIENTS. THE VECTOR TQ OF LENGTH 4 IS INVOLVED IN ADJUSTING THE ! STEPSIZE IN RELATION TO THE TRUNCATION ERROR. ITS VALUES ARE GIVEN BY ! THE PERTST ARRAY. THE VECTORS EL AND TQ BOTH DEPEND ON METH AND NQ. ! COSET ALSO SETS MAXDER, THE MAXIMUM ORDER OF THE METHOD AVAILABLE. THE ! COEFFICIENTS IN PERTST NEED TO BE GIVEN TO ONLY ABOUT ONE PERCENT ! ACCURACY. THE ORDER IN WHICH THE GROUPS APPEAR BELOW IS COEFFICIENTS ! FOR ORDER NQ-1, COEFFICIENTS FOR ORDER NQ, COEFFICIENTS FOR ORDER NQ+1. ! ------------------------------------------------------------------------ ! .. ! .. LOCAL SCALARS .. INTEGER :: k ! .. ! .. LOCAL ARRAYS .. REAL (dp), SAVE :: pertst(8,3) = RESHAPE( (/ 1.0_dp, 2.0_dp, 4.5_dp, & 7.333_dp, 10.42_dp, 13.7_dp, 17.15_dp, & 20.74_dp, 2.0_dp, 4.5_dp, 7.333_dp, & 10.42_dp, 13.7_dp, 17.15_dp, 20.74_dp, & 24.46_dp, 4.5_dp, 7.333_dp, 10.42_dp, & 13.7_dp, 17.15_dp, 20.74_dp, 24.46_dp, & 1.0_dp /), (/ 8, 3 /) ) ! .. ! ------------------------------------------------------------------- ! THE FOLLOWING COEFFICIENTS SHOULD BE DEFINED TO MACHINE ACCURACY. ! THEIR DERIVATION IS GIVEN IN REFERENCE 2. ! ------------------------------------------------------------------- IF (nq > maxord) maxord = nq maxder = 6 ncoset = ncoset + 1 SELECT CASE (nq) CASE (1) el(1) = 1.0D+0 elst(1) = 1.5D+0 elst(3) = -0.5D+0 CASE (2) el(1) = 6.6666666666667D-01 el(3) = 3.3333333333333D-01 elst(1) = 9.5652173913043D-01 elst(3) = 2.1739130434782D-01 elst(4) = -1.7391304347826D-01 CASE (3) el(1) = 5.4545454545455D-01 el(3) = 4.5454545454545D-01 el(4) = 1.8181818181818D-01 elst(1) = 7.6142131979695D-01 elst(3) = 3.2994923857868D-01 elst(4) = 8.6294416243654D-02 elst(5) = -9.1370558375634D-02 CASE (4) el(1) = 0.48D+0 el(3) = 0.52D+0 el(4) = 0.28D+0 el(5) = 0.12D+0 elst(1) = 6.5733706517393D-01 elst(3) = 4.0023990403838D-01 elst(4) = 1.5793682526989D-01 elst(5) = 4.4382247101159D-02 elst(6) = -5.7576969212315D-02 CASE (5) el(1) = 4.3795620437956D-01 el(3) = 5.62043795620436D-01 el(4) = 3.43065693430656D-01 el(5) = 1.97080291970802D-01 el(6) = 8.75912408759123D-02 elst(1) = 5.9119243917152D-01 elst(3) = 4.4902473356122D-01 elst(4) = 2.1375427307460D-01 elst(5) = 9.0421610027481503D-02 elst(6) = 2.6409276761177D-02 elst(7) = -4.0217172732757D-02 CASE (6) el(1) = 4.08163265306120D-01 el(3) = 5.91836734693874D-01 el(4) = 3.87755102040813D-01 el(5) = 2.51700680272107D-01 el(6) = 1.49659863945577D-01 el(7) = 6.80272108843534D-02 elst(1) = 5.4475876041119D-01 elst(3) = 4.8525549636077D-01 elst(4) = 2.5789750131312D-01 elst(5) = 1.3133738525800D-01 elst(6) = 5.7677396763462D-02 elst(7) = 1.7258197643881D-02 elst(8) = -3.0014256771967D-02 CASE (7) el(1) = 3.85674931129476D-01 el(3) = 6.14325068870521D-01 el(4) = 4.21487603305783D-01 el(5) = 2.9292929292929D-01 el(6) = 1.96510560146923D-01 el(7) = 1.19375573921028D-01 el(8) = 5.50964187327820D-02 elst(1) = 5.0999746293734D-01 elst(3) = 5.1345839935281D-01 elst(4) = 2.9364346131937D-01 elst(5) = 1.6664672120553D-01 elst(6) = 8.8013735242353D-02 elst(7) = 3.9571794884069D-02 elst(8) = 1.2039080338722D-02 elst(9) = -2.3455862290154D-02 END SELECT DO k = 1,3 tq(k) = pertst(nq,k) END DO tq(4) = 0.5D+0*tq(2) / nq RETURN ! --------------------- END OF SUBROUTINE COSET --------------------- END SUBROUTINE coset SUBROUTINE pset(y,n,h,t,uround,epsjac,con,miter,ier,nrenew,ymax, & save2,pw,pwcopy,wrkspc,ipiv,fcn) IMPLICIT NONE REAL (dp), INTENT(IN OUT) :: y(:,:) ! y(n,15) INTEGER, INTENT(IN) :: n REAL (dp), INTENT(IN OUT) :: h REAL (dp), INTENT(IN OUT) :: t REAL (dp), INTENT(IN) :: uround REAL (dp), INTENT(IN) :: epsjac REAL (dp), INTENT(IN) :: con INTEGER, INTENT(IN OUT) :: miter INTEGER, INTENT(OUT) :: ier INTEGER, INTENT(IN OUT) :: nrenew REAL (dp), INTENT(IN) :: ymax(:) REAL (dp), INTENT(IN) :: save2(:) REAL (dp), INTENT(OUT) :: pw(:,:) REAL (dp), INTENT(IN OUT) :: pwcopy(:) REAL (dp), INTENT(IN OUT) :: wrkspc(:) INTEGER, INTENT(IN OUT) :: ipiv(:) INTERFACE SUBROUTINE fcn(t,y,ydot) USE common703, ONLY: dp IMPLICIT NONE REAL (dp), INTENT(IN) :: t REAL (dp), INTENT(IN) :: y(:) REAL (dp), INTENT(OUT) :: ydot(:) END SUBROUTINE fcn END INTERFACE ! ------------------------------------------------------------------- ! PSET IS CALLED BY STIFF TO COMPUTE AND PROCESS THE COEFFICIENT ! MATRIX I - H*EL(1)*J WHERE J IS AN APPROXIMATION TO ! THE RELEVANT JACOBIAN. THIS MATRIX IS THEN SUBJECTED TO LU ! DECOMPOSITION IN PREPARATION FOR LATER SOLUTION OF LINEAR SYSTEMS ! OF ALGEBRAIC EQUATIONS WITH LU AS THE COEFFICIENT MATRIX. THE ! MATRIX J IS FOUND BY THE USER-SUPPLIED ROUTINE PDERV IF MITER=1 ! OR BY FINITE DIFFERENCING IF MITER = 2. ! IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION WITH ! PSET USES THE FOLLOWING .. ! EPSJAC = DSQRT(UROUND), USED IN NUMERICAL JACOBIAN INCREMENTS. ! ******************************************************************* ! THE ARGUMENT NRENEW IS USED TO SIGNAL WHETHER OR NOT ! WE REQUIRE A NEW JACOBIAN TO BE CALCULATED. ! IF NRENEW > 0 THEN WE REQUIRE A NEW J TO BE COMPUTED ! = 0 THEN USE A COPY OF THE LAST J COMPUTED ! ******************************************************************* ! .. ! .. LOCAL SCALARS .. INTEGER :: i, j, j1, jjkk REAL (dp) :: d, r, r0, tempry, yj npset = npset + 1 IF (nrenew == 0) THEN pw = RESHAPE( pwcopy(1:n*n) * con, (/ n, n /) ) GO TO 70 END IF IF (miter == 2) GO TO 30 CALL pderv(t,y,pw) pwcopy(1:n*n) = RESHAPE( pw(1:n,1:n), (/ n*n /) ) pw = con * pw nje = nje + 1 GO TO 70 30 nje = nje + 1 d = 0.0D+0 DO i = 1,n d = d + save2(i)**2 END DO r0 = ABS(h)*SQRT(d)*1.0D+03*uround j1 = 0 DO j = 1,n yj = y(j,1) r = epsjac*ymax(j) r = MAX(r,r0) y(j,1) = y(j,1) + r d = con/r CALL fcn(t, y(:,1), wrkspc) DO i = 1,n jjkk = i + j1 tempry = (wrkspc(i) - save2(i)) pwcopy(jjkk) = tempry/r pw(i,j) = tempry*d END DO y(j,1) = yj j1 = j1 + n END DO nfe = nfe + n 70 DO i = 1,n pw(i,i) = pw(i,i) + 1.0D+0 END DO CALL dec(n, pw, ipiv, ier) ! NDEC = NDEC + 1 RETURN ! ---------------------- END OF SUBROUTINE PSET --------------------- END SUBROUTINE pset SUBROUTINE dec(n, a, ip, ier) IMPLICIT NONE INTEGER, INTENT(IN) :: n REAL (dp), INTENT(IN OUT) :: a(:,:) ! a(ndim,n) INTEGER, INTENT(OUT) :: ip(:) INTEGER, INTENT(OUT) :: ier ! ------------------------------------------------------------------- ! MATRIX TRIANGULARISATION BY GAUSSIAN ELIMINATION ! INPUT.. ! N = ORDER OF MATRIX. ! NDIM = DECLARED DIMENSION OF ARRAY A. ! A = MATRIX TO BE TRIANGULARISED. ! OUTPUT.. ! A(I,J), I.LE.J = UPPER TRIANGULAR FACTOR, U. ! A(I,J), I.GT.J = MULTIPLIERS = LOWER TRIANGULAR FACTOR, I-L. ! IP(K), K.LT.N = INDEX OF KTH PIVOT ROW. ! IP(N) = (-1)**(NUMBER OF INTERCHANGES) OR 0. ! IER = 0 IF MATRIX IS NON-SINGULAR, OR K IF FOUND TO BE SINGULAR ! AT STAGE K. ! USE SOL TO OBTAIN SOLUTION OF LINEAR SYSTEM. ! DETERM(A) = IP(N)*A(1,1)*A(2,2)* . . . *A(N,N). ! IF IP(N) = 0, A IS SINGULAR, SOL WILL DIVIDE BY ZERO. ! REFERENCE. ! C.B. MOLER, ALGORITHM 423, LINEAR EQUATION SOLVER, C.A.C.M ! 15 (1972), P.274. ! ------------------------------------------------------------------ ! .. ! .. LOCAL SCALARS .. INTEGER :: i,j,k,kp1,m,nm1 ! .. ier = 0 ip(n) = 1 IF (n == 1) GO TO 70 nm1 = n - 1 DO k = 1,nm1 kp1 = k + 1 m = k DO i = kp1,n IF (ABS(a(i,k)) > ABS(a(m,k))) m = i END DO ip(k) = m t = a(m,k) IF (m == k) GO TO 20 ip(n) = -ip(n) a(m,k) = a(k,k) a(k,k) = t 20 IF (t == 0.0D+0) GO TO 80 t = 1.0D+0/t DO i = kp1,n a(i,k) = -a(i,k)*t END DO DO j = kp1,n t = a(m,j) a(m,j) = a(k,j) a(k,j) = t IF (t == 0.0D+0) CYCLE DO i = kp1,n a(i,j) = a(i,j) + a(i,k)*t END DO END DO END DO 70 k = n IF (a(n,n) == 0.0D+0) GO TO 80 RETURN 80 ier = k ip(n) = 0 RETURN ! --------------------- END OF SUBROUTINE DEC ---------------------- END SUBROUTINE dec SUBROUTINE sol(n, a, b, ip) IMPLICIT NONE INTEGER, INTENT(IN) :: n REAL (dp), INTENT(IN) :: a(:,:) ! a(ndim,n) REAL (dp), INTENT(IN OUT) :: b(:) INTEGER, INTENT(IN) :: ip(:) ! .. ! .. LOCAL SCALARS .. INTEGER :: i,k,kb,km1,kp1,m,nm1 ! .. ! ------------------------------------------------------------------ ! SOLUTION OF LINEAR SYSTEM, A*X = B. ! INPUT .. ! N = ORDER OF MATRIX. ! NDIM = DECLARED DIMENSION OF MATRIX A. ! A = TRIANGULARISED MATRIX OBTAINED FROM DEC. ! B = RIGHT HAND SIDE VECTOR. ! IP = PIVOT VECTOR OBTAINED FROM DEC. ! DO NOT USE IF DEC HAS SET IER .NE. 0 ! OUTPUT.. ! B = SOLUTION VECTOR, X. ! ------------------------------------------------------------------ IF (n == 1) GO TO 50 nm1 = n - 1 DO k = 1,nm1 kp1 = k + 1 m = ip(k) t = b(m) b(m) = b(k) b(k) = t DO i = kp1,n b(i) = b(i) + a(i,k)*t END DO END DO DO kb = 1,nm1 km1 = n - kb k = km1 + 1 b(k) = b(k)/a(k,k) t = -b(k) DO i = 1,km1 b(i) = b(i) + a(i,k)*t END DO END DO 50 b(1) = b(1)/a(1,1) RETURN ! ------------------------- END OF SUBROUTINE SOL ------------------ END SUBROUTINE sol SUBROUTINE errors(n, tq, eps, edn, e, eup, bnd) INTEGER, INTENT(IN) :: n REAL (dp), INTENT(IN) :: tq(:) REAL (dp), INTENT(IN) :: eps REAL (dp), INTENT(OUT) :: edn REAL (dp), INTENT(OUT) :: e REAL (dp), INTENT(OUT) :: eup REAL (dp), INTENT(OUT) :: bnd ! *************************************************** ! THIS ROUTINE CALCULATES ERRORS USED IN TESTS IN STIFF . ! *************************************************** ! .. ! .. LOCAL SCALARS .. ! .. REAL (dp) :: sqhol sqhol = n*eps*eps edn = tq(1)*tq(1)*sqhol ! ** ERROR ASSOCIATED WITH LOWER ORDER METHOD e = tq(2)*tq(2)*sqhol ! ** ERROR ASSOCIATED WITH PRESENT ORDER eup = tq(3)*tq(3)*sqhol ! ** ERROR ASSOCIATED WITH HIGHER ORDER METHOD bnd = tq(4)*tq(4)*sqhol*0.3D+0 RETURN END SUBROUTINE errors SUBROUTINE prdict(t, h, y, l, n, yprime, fcn) IMPLICIT NONE REAL (dp), INTENT(IN OUT) :: t REAL (dp), INTENT(IN) :: h REAL (dp), INTENT(OUT) :: y(:,:) ! y(n,15) INTEGER, INTENT(IN) :: l INTEGER, INTENT(IN) :: n REAL (dp), INTENT(IN OUT) :: yprime(:) INTERFACE SUBROUTINE fcn(t,y,ydot) USE common703, ONLY: dp IMPLICIT NONE REAL (dp), INTENT(IN) :: t REAL (dp), INTENT(IN) :: y(:) REAL (dp), INTENT(OUT) :: ydot(:) END SUBROUTINE fcn END INTERFACE ! ********************************************************************** ! PREDICTS A VALUE FOR Y AT (T+H) GIVEN THE HISTORY ARRAY AT Y(T) ! THEN EVALUATES THE DERIVATIVE AT THIS POINT, THE RESULT OF THIS ! EVALUATION BEING STORED IN YPRIME() ! ********************************************************************** ! .. ! .. LOCAL SCALARS .. INTEGER :: i,j2 ! .. DO j2 = 2,l DO i = 1,n y(i,1) = y(i,1) + y(i,j2) END DO END DO t = t + h CALL fcn(t, y(:,1), yprime) nfe = nfe + 1 RETURN END SUBROUTINE prdict SUBROUTINE itrat2(y,n,t,hbeta,errbnd,arh,crate,tcrate,m,worked, & ymax,error,save1,save2,pw,ipiv,lmb,fcn) IMPLICIT NONE REAL (dp), INTENT(IN) :: y(:,:) ! y(n,15) INTEGER, INTENT(IN) :: n REAL (dp), INTENT(IN OUT) :: t REAL (dp), INTENT(IN) :: hbeta REAL (dp), INTENT(IN OUT) :: errbnd REAL (dp), INTENT(IN) :: arh(:) REAL (dp), INTENT(IN OUT) :: crate REAL (dp), INTENT(OUT) :: tcrate INTEGER, INTENT(OUT) :: m LOGICAL, INTENT(OUT) :: worked REAL (dp), INTENT(IN) :: ymax(:) REAL (dp), INTENT(OUT) :: error(:) REAL (dp), INTENT(IN OUT) :: save1(:) REAL (dp), INTENT(IN OUT) :: save2(:) REAL (dp), INTENT(IN OUT) :: pw(:,:) INTEGER, INTENT(IN OUT) :: ipiv(:) INTEGER, INTENT(IN) :: lmb INTERFACE SUBROUTINE fcn(t,y,ydot) USE common703, ONLY: dp IMPLICIT NONE REAL (dp), INTENT(IN) :: t REAL (dp), INTENT(IN) :: y(:) REAL (dp), INTENT(OUT) :: ydot(:) END SUBROUTINE fcn END INTERFACE ! .. ! .. LOCAL SCALARS .. INTEGER :: i REAL (dp) :: d, d1 ! .. ! .. DATA STATEMENTS .. REAL (dp), PARAMETER :: zero = 0.0_dp ! .. IF(lmb == 1) GO TO 25 DO i = 1,n save1(i) = -save1(i) + hbeta*save2(i) + arh(i) END DO CALL sol(n, pw, save1, ipiv) nbsol = nbsol + 1 d = zero DO i = 1,n error(i) = error(i) + save1(i) d = d + (save1(i)/ymax(i))**2 save1(i) = y(i,1) + error(i) END DO tcrate = tcrate + crate d1 = d m = 1 CALL fcn(t,save1,save2) nfe = nfe + 1 25 worked = .true. 30 DO i = 1,n save1(i) = -save1(i) + hbeta*save2(i) + arh(i) END DO ! IF WE ARE HERE THEN PARTIALS ARE O.K. CALL sol(n, pw, save1, ipiv) nbsol = nbsol + 1 ! WE NOW CALCULATE A WEIGHTED RMS TYPE NORM d = zero DO i = 1,n error(i) = error(i) + save1(i) d = d + (save1(i)/ymax(i))**2 save1(i) = y(i,1) + error(i) END DO ! ------------------------------------------------------------------- ! TEST FOR CONVERGENCE. IF M.GT.0 , AN ESTIMATE OF THE CONVERGENCE ! RATE CONSTANT IS STORED IN CRATE, AND THIS IS USED IN THE TEST. ! ------------------------------------------------------------------- IF (m /= 0) THEN IF (d1 /= zero) crate = MAX(0.9D+0*crate,d/d1) END IF tcrate = tcrate + crate IF ((d*MIN(1.0D+0, 2.0D+0*crate)) < errbnd) RETURN IF (m /= 0) THEN IF (d > d1) THEN worked = .false. RETURN END IF END IF d1 = d IF (m == 3) THEN worked = .false. RETURN END IF m = m + 1 CALL fcn(t,save1,save2) nfe = nfe + 1 GO TO 30 END SUBROUTINE itrat2 SUBROUTINE stiff(eps, h, hmax, hmin, jstart, kflag, mf, t, tend, y, n, ymax, & pw, pwcopy, yhold, ipiv, lout, fcn) IMPLICIT NONE REAL (dp), INTENT(IN) :: eps REAL (dp), INTENT(IN OUT) :: h REAL (dp), INTENT(IN) :: hmax REAL (dp), INTENT(IN) :: hmin INTEGER, INTENT(OUT) :: jstart INTEGER, INTENT(OUT) :: kflag INTEGER, INTENT(IN) :: mf REAL (dp), INTENT(IN OUT) :: t REAL (dp), INTENT(IN) :: tend REAL (dp), INTENT(OUT) :: y(:,:) ! y(n,15) INTEGER, INTENT(IN) :: n REAL (dp), INTENT(IN) :: ymax(:) REAL (dp), INTENT(IN OUT) :: pw(:,:) REAL (dp), INTENT(IN OUT) :: pwcopy(:) REAL (dp), INTENT(IN OUT) :: yhold(:,:) ! yhold(n,15) INTEGER, INTENT(IN OUT) :: ipiv(:) INTEGER, INTENT(IN) :: lout INTERFACE SUBROUTINE fcn(t,y,ydot) USE common703, ONLY: dp IMPLICIT NONE REAL (dp), INTENT(IN) :: t REAL (dp), INTENT(IN) :: y(:) REAL (dp), INTENT(OUT) :: ydot(:) END SUBROUTINE fcn END INTERFACE ! ------------------------------------------------------------------ ! STIFF PERFORMS ONE STEP OF THE INTEGRATION OF AN INITIAL VALUE ! PROBLEM FOR A SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS. ! COMMUNICATION WITH STIFF IS DONE WITH THE FOLLOWING VARIABLES.. ! Y AN N BY LMAX+4 ARRAY CONTAINING THE DEPENDENT VARIABLES ! AND THEIR BACKWARD DIFFERENCES. LMAX-1 =MAXDER IS THE ! MAXIMUM ORDER AVAILABLE. SEE SUBROUTINE COSET. ! Y(I,J+1) CONTAINS THE JTH BACKWARD DIFFERENCE OF Y(I) ! T THE INDEPENDENT VARIABLE. T IS UPDATED ON EACH STEP TAKEN. ! H THE STEPSIZE TO BE ATTEMPTED ON THE NEXT STEP. ! H IS ALTERED BY THE ERROR CONTROL ALGORITHM DURING ! THE PROBLEM. H CAN BE EITHER POSITIVE OR NEGATIVE BUT ! ITS SIGN MUST REMAIN CONSTANT THROUGHOUT THE PROBLEM. ! HMIN THE MINIMUM AND MAXIMUM ABSOLUTE VALUE OF THE STEPSIZE ! HMAX TO BE USED FOR THE STEP. THESE MAY BE CHANGED AT ANY ! TIME BUT WILL NOT TAKE EFFECT UNTIL THE NEXT H CHANGE. ! EPS THE RELATIVE ERROR BOUND. SEE DESCRIPTION IN MEBDF ! UROUND THE UNIT ROUNDOFF OF THE MACHINE. ! N THE NUMBER OF FIRST ORDER DIFFERENTIAL EQUATIONS. ! MF THE METHOD FLAG. MUST BE SET TO 21 OR 22 AT PRESENT ! KFLAG A COMPLETION CODE WITH THE FOLLOWING MEANINGS.. ! 0 THE STEP WAS SUCCESSFUL ! -1 THE REQUESTED ERROR COULD NOT BE ACHIEVED ! WITH ABS(H) = HMIN. ! -2 THE REQUESTED ERROR IS SMALLER THAN CAN ! BE HANDLED FOR THIS PROBLEM. ! -3 CORRECTOR CONVERGENCE COULD NOT BE ! ACHIEVED FOR ABS(H)=HMIN. ! ON A RETURN WITH KFLAG NEGATIVE, THE VALUES OF T AND THE Y ARRAY ARE ! AS AT THE BEGINNING OF THE LAST STEP ATTEMPTED, AND H IS THE LAST ! STEP SIZE ATTEMPTED. ! JSTART AN INTEGER USED ON INPUT AND OUTPUT. ! ON INPUT IT HAS THE FOLLOWING VALUES AND MEANINGS.. ! 0 PERFORM THE FIRST STEP. ! .GT.0 TAKE A NEW STEP CONTINUING FROM THE LAST ! .LT.0 TAKE THE NEXT STEP WITH A NEW VALUE OF ! H, EPS OR N ! ON EXIT, JSTART IS NQUSED, THE ORDER OF THE METHOD LAST USED. ! YMAX AN ARRAY OF N ELEMENTS WITH WHICH THE ESTIMATED LOCAL ! ERRORS IN Y ARE COMPARED ! PW A BLOCK OF LOCATIONS USED FOR PARTIAL DERIVATIVES ! IPIV AN INTEGER ARRAY OF LENGTH N USED FOR PIVOT INFORMATION. ! JNEWIM IS TO INDICATE IF PRESENT ITERATION MATRIX ! WAS FORMED USING A NEW J OR OLD J. ! JSNOLD KEEPS TRACK OF NUMBER OF STEPS TAKEN WITH ! PRESENT ITERATION MATRIX (BE IT FORMED BY ! A NEW J OR NOT). ! AVNEWJ STORES VALUE FOR AVERAGE CRATE WHEN ITERATION ! MATRIX WAS FORMED BY A NEW J. ! AVOLDJ STORES VALUE FOR AVERAGE CRATE WHEN ITERATION ! MATRIX WAS FORMED BY AN OLD J. ! NRENEW FLAG THAT IS USED IN COMMUNICATION WITH SUBROUTINE PSET. ! IF NRENEW > 0 THEN FORM A NEW JACOBIAN BEFORE ! COMPUTING THE COEFFICIENT MATRIX FOR ! THE NEWTON ITERATION ! = 0 FORM THE COEFFICIENT MATRIX USING A ! COPY OF AN OLD JACOBIAN ! NEWPAR FLAG USED IN THIS SUBROUTINE TO INDICATE IF A JACOBIAN ! HAS BEEN EVALUATED FOR THE CURRENT STEP ! ********************************************************************** ! .. ! .. LOCAL SCALARS .. INTEGER :: i, ier, iiter, iiter2, ijus, itst, j, j1, j2, kfail, ll, & m3step, maxder, meth, miter, newq LOGICAL :: finish, worked REAL (dp) :: atol, d, ddown, dstep2, dup, ebdf1, ebdf2, efail, enq1, enq2, & enq3, prbdf1, prbdf2, prfail, pr1, pr2, pr3, qq, red, rhbdf1, & rhbdf2, rhlmit, trange ! .. ! .. LOCAL ARRAYS .. REAL (dp), SAVE :: el(10), elst(10), tq(4) REAL (dp) :: arh(n), error(n), save1(n), save2(n), ynhold(n,2) ! .. ! .. SAVE STATEMENT .. REAL (dp), SAVE :: epsold, hold, edn, eup, bnd, e INTEGER, SAVE :: mfold, lmax ! .. ! .. DATA STATEMENTS .. REAL (dp), SAVE :: oldlo = 1.0_dp REAL (dp), PARAMETER :: zero = 0.0_dp, one = 1.0_dp ! .. el(2) = one elst(2) = one told = t kflag = 0 IF (jstart > 0) GO TO 60 IF (jstart /= 0) GO TO 30 ! ------------------------------------------------------------------ ! ON THE FIRST CALL, THE ORDER IS SET TO 1 AND THE INITIAL YDOT ! IS CALCULATED. RMAX IS THE MAXIMUM RATIO BY WHICH H CAN BE ! INCREASED IN A SINGLE STEP. RMAX IS SET EQUAL TO 1.D4 INITIALLY ! TO COMPENSATE FOR THE SMALL INITIAL H, BUT THEN IS NORMALLY = 10. ! IF A FAILURE OCCURS (IN CORRECTOR CONVERGENCE OR ERROR TEST), ! RMAX IS SET AT 2. FOR THE NEXT INCREASE. ! ------------------------------------------------------------------ CALL fcn(t, y(:,1), save1) DO i = 1,n y(i,2) = h*save1(i) END DO meth = 2 miter = mf - 10*meth nq = 1 l = 2 idoub = 3 kfail = 0 rmax = 10000.0D+0 rc = zero crate1 = one crate2 = one jsnold = 0 jnewim = .true. tcrat1 = zero tcrat2 = zero top=0 atol=eps/10.0D+0 DO i=1,12 hstpsz(1,i)=1.0D+0 hstpsz(2,i)=atol END DO hold = h mfold = mf nstep = 0 nfe = 1 nje = 0 ndec = 0 npset = 0 ncoset = 0 maxord = 1 nfail1 = 0 nfail2 = 0 cfail = .true. avnewj = zero avoldj = zero avnew2 = zero avold2 = zero sample = .false. isamp = 0 ! ************************************************** ! CFAIL=.TRUE. ENSURES THAT WE CALCULATE A NEW ! J ON THE FIRST CALL. ! ************************************************** meqc1 = 0 meqc2 = 0 mq1tmp = 0 mq2tmp = 0 nbsol = 0 ! ----------------------------------------------------------------- ! IF THE CALLER HAS CHANGED N OR EPS, THE CONSTANTS E, EDN, EUP ! AND BND MUST BE RESET. E IS A COMPARISON FOR ERRORS AT THE ! CURRENT ORDER NQ. EUP IS TO TEST FOR INCREASING THE ORDER, ! EDN FOR DECREASING THE ORDER. BND IS USED TO TEST FOR CONVERGENCE ! OF THE CORRECTOR ITERATES. IF THE CALLER HAS CHANGED H, Y MUST ! BE RE-SCALED. IF H IS CHANGED, IDOUB IS SET TO L+1 TO PREVENT ! FURTHER CHANGES IN H FOR THAT MANY STEPS. ! ----------------------------------------------------------------- CALL coset(nq,el,elst,tq,maxder) lmax = maxder + 1 rc = rc*el(1)/oldlo oldlo = el(1) iweval = miter nrenew = 1 newpar = 0 ! ***************************************************** ! NRENEW AND NEWPAR ARE TO INSTRUCT ROUTINE THAT ! WE WISH A NEW J TO BE CALCULATED FOR THIS STEP. ! ***************************************************** CALL errors(n,tq,eps,edn,e,eup,bnd) epsold = eps DO i = 1,n arh(i) = el(2)*y(i,1) END DO yhold(1:n,1:l) = y(1:n,1:l) qi = h*el(1) qq = one/qi CALL prdict(t,h,y,l,n,save2,fcn) GO TO 110 ! >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> ! DIFFERENT PARAMETERS ON THIS CALL < ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 30 y(1:n,1:l) = yhold(1:n,1:l) IF (mf /= mfold) THEN meth = mf/10 miter = mf - 10*meth mfold = mf iweval = miter END IF IF (eps /= epsold) THEN CALL errors(n,tq,eps,edn,e,eup,bnd) epsold = eps END IF IF (h /= hold) THEN rh = h/hold h = hold GO TO 50 ELSE GO TO 60 END IF ! ********************************************* ! RE-SCALE Y AFTER A CHANGE OF STEPSIZE * ! ********************************************* 40 rh = MAX(rh,hmin/ABS(h)) 50 rh = MIN(rh,hmax/ABS(h),rmax) CALL rscale(n,l,rh,y) rmax = 10.0D+0 jchang = 1 h = h*rh rc = rc*rh IF (jsnold > 40) THEN cfail = .true. newpar = 0 rc = zero ! ********************************************************************** ! CFAIL=TRUE AND NEWPAR=0 SHOULD FORCE A NEW J TO BE EVALUATED ! AFTER 42 STEPS WITH AN OLD J, IF WE HAVE HAD A FAILURE OF ANY ! KIND ON THE FIRST, SECOND OR THIRD STAGE OF THE CURRENT STEP ! ********************************************************************** END IF idoub = l + 1 yhold(1:n,1:l) = y(1:n,1:l) 60 IF (ABS(rc-one) > 0.3D+0) iweval = miter ! ------------------------------------------------------------------ ! THIS SECTION COMPUTES THE PREDICTED VALUES OF Y ! AND THE RHS, ARH, FOR USE IN THE NEWTON ITERATION SCHEME. ! RC IS THE RATIO OF THE NEW TO OLD VALUES OF THE COEFFICIENT ! H*EL(1). WHEN RC DIFFERS FROM 1 BY MORE THAN 30 PERCENT, IWEVAL IS ! SET TO MITER TO FORCE THE PARTIALS TO BE UPDATED. ! ------------------------------------------------------------------ qi = h*el(1) qq = one/qi DO i = 1,n arh(i) = el(2)*y(i,1) END DO DO j1 = 2,nq DO i = 1,n arh(i) = arh(i) + el(j1+1)*y(i,j1) END DO END DO IF (jchang == 1) THEN ! IF WE HAVE CHANGED STEPIZE THEN PREDICT A VALUE FOR Y(T+H) ! AND EVALUATE THE DERIVATIVE THERE (STORED IN SAVE2()) CALL prdict(t,h,y,l,n,save2,fcn) ELSE ! ELSE USE THE VALUES COMPUTED FOR THE SECOND BDF FROM THE LAST ! STEP. Y( ,LMAX+3) HOLDS THE VALUE FOR THE DERIVATIVE AT (T+H) ! AND Y( ,LMAX+4) HOLDS THE APPROXIMATION TO Y AT THIS POINT. DO i = 1,n save2(i) = y(i,lmax+3) y(i,1) = y(i,lmax+4) END DO t = t + h END IF 110 IF (iweval <= 0) GO TO 120 ! ------------------------------------------------------------------- ! IF INDICATED, THE MATRIX P = I - H*EL(1)*J IS RE-EVALUATED BEFORE ! STARTING THE CORRECTOR ITERATION. IWEVAL IS SET = 0 TO INDICATE ! THAT THIS HAS BEEN DONE. P IS COMPUTED AND PROCESSED IN PSET. ! THE PROCESSED MATRIX IS STORED IN PW ! ------------------------------------------------------------------- iweval = 0 rc = one iiter = meqc1 - mq1tmp iiter2 = meqc2 - mq2tmp IF (jnewim) THEN IF (jsnold >= 3) THEN avnewj = tcrat1/DBLE(iiter) avnew2 = tcrat2/DBLE(iiter2) ELSE avnewj = one avnew2 = one END IF ELSE ! MATRIX P WAS FORMED WITH A COPY OF J IF (jsnold >= 3) THEN avoldj = tcrat1/DBLE(iiter) avold2 = tcrat2/DBLE(iiter2) IF (avoldj < avnewj) THEN avnewj = avoldj ELSE IF (((ABS(avoldj-avnewj)) > 0.3D+0) .OR. & ((avoldj > 0.85D+0).AND. (avoldj /= one))) THEN ! SINCE IN CERTAIN INSTANCES AVOLDJ WILL ! BE 1.0 AND THERE WILL BE NO NEED TO UPDATE J. cfail = .true. crate1 = one crate2 = one END IF ELSE cfail = .true. crate1 = one crate2 = one ! ********************************************************* ! IF WE HAVE REACHED HERE THINGS MUST HAVE GONE WRONG ! ********************************************************* END IF END IF tcrat1 = zero tcrat2 = zero IF (cfail) THEN IF (newpar == 1) THEN nrenew = 0 jnewim = .true. ELSE nrenew = 1 newpar = 1 jsinup = -1 jnewim = .true. END IF ELSE nrenew = 0 jnewim = .false. END IF cfail = .false. jsnold = 0 mq1tmp = meqc1 mq2tmp = meqc2 CALL pset(y,n,h,t,uround,epsjac,-qi,miter,ier,nrenew,ymax,save2, & pw,pwcopy,error,ipiv,fcn) ! NOTE THAT ERROR() IS JUST BEING USED AS A WORKSPACE BY PSET IF (ier /= 0) THEN ! IF IER>0 THEN WE HAVE A SINGULARITY IN THE COEFFICIENT MATRIX ijus=1 red=0.5D+0 GO TO 450 END IF 120 DO i = 1,n save1(i) = y(i,1) error(i) = zero END DO m1 = 0 ! ********************************************************************** ! UP TO 4 CORRECTOR ITERATIONS ARE TAKEN. A CONVERGENCE TEST IS ! MADE ON THE R.M.S. NORM OF EACH CORRECTION ,USING BND, WHICH ! DEPENDS ON EPS. THE SUM OF THE CORRECTIONS IS ACCUMULATED IN THE ! VECTOR ERROR(I). THE Y ARRAY IS NOT ALTERED IN THE CORRECTOR ! LOOP. THE UPDATED Y VECTOR IS STORED TEMPORARILY IN SAVE1. ! ********************************************************************** IF (.NOT.sample) THEN CALL itrat2(y,n,t,qi,bnd,arh,crate1,tcrat1,m1,worked,ymax, & error,save1,save2,pw,ipiv,1,fcn) itst = 2 ELSE CALL itrat2(y,n,t,qi,bnd,arh,crate1,tcrat1,m1,worked,ymax, & error,save1,save2,pw,ipiv,0,fcn) itst = 3 END IF meqc1 = meqc1 + m1 + 1 ! NOW TEST TO SEE IF IT WAS SUCCESSFUL OR NOT IF (.NOT.worked) THEN nfail1 = nfail1 + 1 ! ********************************************************************** ! THE CORRECTOR ITERATION FAILED TO CONVERGE IN 4 TRIES. IF ! PARTIALS ARE NOT UP TO DATE, THEY ARE RE-EVALUATED FOR THE ! NEXT TRY. OTHERWISE THE Y ARRAY IS RETRACTED TO ITS VALUES ! BEFORE PREDICTION AND H IS REDUCED IF POSSIBLE. IF NOT A ! NON-CONVERGENCE EXIT IS TAKEN ! ********************************************************************** IF (iweval == -1) THEN ! HAVE BEEN USING OLD PARTIALS, UPDATE THEM AND TRY AGAIN iweval = miter cfail = .true. CALL fcn(t, y(:,1), save2) nfe = nfe + 1 GO TO 110 END IF ijus=0 red=0.5D+0 GO TO 450 END IF iweval = -1 hused = h nqused = nq DO i = 1,n save2(i) = (save1(i)-arh(i))*qq y(i,1) = save1(i) END DO ! UPDATE THE DIFFERENCES AT N+1 DO j = 2,l DO i = 1,n y(i,j) = y(i,j-1) - yhold(i,j-1) END DO END DO ! COMPUTE ERROR IN THE SOLUTION d = zero DO i = 1,n d = d + ((y(i,l)-yhold(i,l))/ymax(i))**2 END DO ! STORING Y FROM FIRST STEP FOR USE IN THIRD STEP DO i = 1,n ynhold(i,1) = y(i,1) ynhold(i,2) = save2(i) END DO IF (d > e) GO TO 330 IF ((m1+m2) >= itst) THEN m2 = 0 ebdf1 = 0.5D+0 / l prbdf1 = ((d/ (e*0.7D+0))**ebdf1)*1.5D+0 + 1.6D-6 rhbdf1 = one/prbdf1 CALL hchose(rhbdf1, h) rhlmit = 1 - (5+itst-m1-m2)*0.5D-1 IF (rhbdf1 < rhlmit) THEN ijus=1 red=rhbdf1 GO TO 450 END IF END IF kfail = 0 ! ---------------------------------------------------------------------- DO i = 1,n arh(i) = el(2)*y(i,1) END DO DO j1 = 2,nq DO i = 1,n arh(i) = arh(i) + el(j1+1)*y(i,j1) END DO END DO CALL prdict(t,h,y,l,n,save2,fcn) DO i = 1,n save1(i) = y(i,1) error(i) = zero END DO m2 = 0 ! FOR NOW WILL ASSUME THAT WE DO NOT WISH TO SAMPLE ! AT THE N+2 STEP POINT CALL itrat2(y,n,t,qi,bnd,arh,crate2,tcrat2,m2,worked,ymax,error, & save1,save2,pw,ipiv,1,fcn) meqc2 = meqc2 + m2 + 1 ! NOW CHECK TO SEE IF IT WAS SUCCESSFUL OR NOT IF (.NOT.worked) THEN nfail2 = nfail2 + 1 ijus=0 red=0.5D+0 GO TO 450 END IF ! IF WE ARE DOWN TO HERE THEN THINGS MUST HAVE CONVERGED DO i = 1,n y(i,lmax+3) = (save1(i)-arh(i))*qq y(i,lmax+4) = save1(i) END DO IF ((m2 >= 2) .AND. (crate2 < 0.35D+0)) THEN DO j = 2,nq DO i = 1,n save1(i) = save1(i) - y(i,j) END DO END DO dstep2 = zero DO i = 1,n dstep2 = dstep2 + ((save1(i)-ynhold(i,1)-y(i,l))/ymax(i))**2 END DO ebdf2 = 0.5D+0 / l prbdf2 = ((dstep2/ (e*0.7D+0))**ebdf2)*1.5D+0 + 1.6D-6 rhbdf2 = one/prbdf2 CALL hchose(rhbdf2, h) IF (rhbdf2 < 0.9D+0) THEN ijus=1 red=rhbdf2 GO TO 450 END IF END IF ! WE ARE NOW COMPUTING THE THIRD STAGE ll = l + 1 t = told + h DO i = 1,n arh(i) = h* (elst(nq+2)*y(i,lmax+3) + (elst(1)-el(1))*ynhold(i,2)) DO j1 = 1,nq arh(i) = arh(i) + elst(j1+1)*yhold(i,j1) END DO END DO DO i = 1,n save2(i) = ynhold(i,2) y(i,1) = ynhold(i,1) END DO m3step = 0 300 DO i = 1,n save1(i) = -y(i,1) + qi*save2(i) + arh(i) END DO CALL sol(n, pw, save1, ipiv) nbsol = nbsol + 1 d = zero DO i = 1,n d = d + (save1(i)/ymax(i))**2 y(i,1) = y(i,1) + save1(i) END DO IF ((d*MIN(one, 2.0D+0*crate1)) <= bnd) GO TO 360 IF (m3step == 4) THEN WRITE (lout,9000) ijus=1 red=0.5D+0 GO TO 450 END IF m3step = m3step + 1 CALL fcn(t, y(:,1), save2) nfe = nfe + 1 GO TO 300 330 kfail = kfail - 1 ! ************************************************************************** ! THE ERROR TEST FAILED. KFAIL KEEPS TRACK OF MULTIPLE FAILURES. RESTORE T ! AND THE Y ARRAY TO THEIR PREVIOUS VALUES AND PREPARE TO TRY THE STEP ! AGAIN. COMPUTE THE OPTIMAL STEP SIZE FOR THIS ORDER AND ONE ORDER LOWER. ! ************************************************************************** t = told hold = h efail = 0.5D+0 / l y(1:n,1:l) = yhold(1:n,1:l) rmax = 2.0D+0 IF (ABS(h) <= hmin*1.00001D+0) THEN ! REQUESTED ERROR NOT POSSIBLE WITH GIVEN HMIN kflag = -1 hold = h RETURN END IF IF (kfail <= -3) GO TO 340 ! PREDICTING A NEW H AFTER INSUFFICIENT ACCURACY prfail = ((d/e)**efail)*1.5D+0 + 1.6D-6 newq = nq rh = one / (prfail*DBLE(-kfail)) CALL hchose(rh, h) GO TO 40 ! ********************************************************************** ! CONTROL REACHES THIS STAGE IF 3 OR MORE FAILURES HAVE OCCURED. ! IT IS ASSUMED THAT THE DERIVATIVES THAT HAVE ACCUMULATED IN THE Y ! ARRAY HAVE ERRORS OF THE WRONG ORDER. HENCE THE FIRST DERIVATIVE ! IS RE-COMPUTED, AND THE ORDER IS SET TO 1. THEN H IS REDUCED BY A ! FACTOR OF 10, AND THE STEP IS RETRIED. AFTER A TOTAL OF 7 ! FAILURES AN EXIT IS TAKEN WITH KFLAG=-2. ! ********************************************************************** 340 IF (kfail == -7) THEN ! ERROR SMALLER THAN CAN BE HANDLED FOR PROBLEM kflag = -2 hold = h RETURN END IF ! ********************************* ! START FROM ORDER 1 AGAIN * ! ********************************* jchang = 1 rh = MAX(hmin/ABS(h),0.1D+0) CALL hchose(rh, h) h = h*rh CALL fcn(t, yhold(:,1), save1) nfe = nfe + 1 DO i = 1,n y(i,1) = yhold(i,1) y(i,2) = h*save1(i) yhold(i,2) = y(i,2) END DO iweval = miter cfail = .true. ! SINCE WE HAVE HAD PROBLEMS PROCEED WITH THIS ORDER ! FOR 10 STEPS (IF WE CAN) idoub = 10 IF (nq == 1) GO TO 60 nq = 1 l = 2 ! RESET ORDER, RECALCULATE ERROR BOUNDS CALL coset(nq,el,elst,tq,maxder) lmax = maxder + 1 rc = rc*el(1)/oldlo oldlo = el(1) CALL errors(n,tq,eps,edn,e,eup,bnd) ! NOW JUMP TO NORMAL CONTINUATION POINT GO TO 60 ! ********************************************************************** ! THE ITERATION FOR THE CORRECTED SOLUTION HAS CONVERGED. ! UPDATE THE Y ARRAY. ! ********************************************************************** 360 DO j2 = 2,ll DO i = 1,n y(i,j2) = y(i,j2-1) - yhold(i,j2-1) END DO END DO ! ---------------------------------------------------------------------- ! IF THE COUNTER IDOUB EQUALS 2 , STORE Y(I,LMAX+5) WHICH IS USED IN ! ASSESSING THE POSSIBILITY OF INCREASING THE ORDER. IF IDOUB = 0 ! CONTROL PASSES TO 480 WHERE AN ATTEMPT TO CHANGE THE STEPSIZE AND ! ORDER IS MADE. ! ---------------------------------------------------------------------- IF (idoub == 2) THEN DO i = 1,n y(i,lmax+5) = y(i,ll) END DO END IF idoub = idoub - 1 trange=(tend-told-h)*h IF(trange < 0.0D+0) GO TO 440 jchang = 0 IF (idoub == 0) THEN sample = .false. isamp = isamp + 1 IF (isamp == 4) THEN sample = .true. isamp = 0 END IF ! ********************************************************************** ! NOW COMPUTE THE FACTORS PR1, PR2 AND PR3, BY WHICH ! H COULD BE DIVIDED AT ORDER NQ-1, ORDER NQ AND ORDER NQ+1 ! RESPECTIVELY. THE SMALLEST OF THESE IS DETERMINED AND THE NEW ! ORDER CHOSEN ACCORDINGLY. IF THE ORDER IS TO BE INCREASED WE ! MUST COMPUTE ONE MORE BACKWARD DIFFERENCE. ! ********************************************************************** pr3 = 1.d+20 IF (l /= lmax) THEN dup = zero DO i = 1,n dup = dup + ((y(i,ll) - y(i,lmax+5)) / ymax(i))**2 END DO enq3 = 0.5D+0 / (l+1) pr3 = ((dup/eup)**enq3)*1.7D+0 + 1.8D-6 END IF enq2 = 0.5D+0 / l d = zero DO i = 1,n d = d + ((y(i,l) - yhold(i,l)) / ymax(i))**2 END DO pr2 = ((d/e)**enq2)*1.5D+0 + 1.6D-6 pr1 = 1.d+20 IF (nq > 1) THEN ddown = zero DO i = 1,n ddown = ddown + (y(i,l)/ymax(i))**2 END DO enq1 = 0.5D+0 / nq pr1 = ((ddown/edn)**enq1)*1.6D+0 + 1.7D-6 END IF IF (pr2 <= pr3) THEN IF (pr2 > pr1) THEN newq = nq - 1 rh = 1.0D+0/pr1 ELSE newq = nq rh = 1.0D+0/pr2 END IF ELSE IF (pr3 < pr1) THEN newq = l rh = 1.0D+0/pr3 ELSE newq = nq - 1 rh = 1.0D+0/pr1 END IF rh = MIN(rh,rmax) CALL hchose(rh, h) IF ((jsinup <= 20).AND.(kflag == 0).AND.(rh < 1.1D+0)) THEN ! WE HAVE RUN INTO PROBLEMS idoub = 10 nq = nqused l = nq + 1 GO TO 440 END IF ! ********************************************************************** ! IF THERE IS A CHANGE IN ORDER, RESET NQ, L AND THE ! COEFFICIENTS. IN ANY CASE H IS RESET ACCORDINGLY TO RH AND THE ! Y ARRAY IS RE-SCALED ! ********************************************************************** IF (newq /= nq) THEN IF (newq > nq) THEN ! ADD AN EXTRA TERM TO THE HISTORY ARRAY DO i = 1,n y(i,ll) = y(i,l) - yhold(i,l) END DO END IF nq = newq l = nq + 1 ! RESET ORDER,RECALCULATE ERROR BOUNDS CALL coset(nq,el,elst,tq,maxder) lmax = maxder + 1 rc = rc*el(1)/oldlo oldlo = el(1) CALL errors(n,tq,eps,edn,e,eup,bnd) END IF ! NOW RESCALE THE HISTORY ARRAY FOR THE NEW STEPSIZE rh = MAX(rh,hmin/ABS(h)) rh = MIN(rh,hmax/ABS(h),rmax) CALL rscale(n,l,rh,y) rmax = 10.0D+0 jchang = 1 h = h*rh rc = rc*rh IF (jsnold > 40) THEN rc = zero ! IF WE HAVE BEEN USING THE SAME COEFFICIENT MATRIX FOR MORE ! THAN 40 STEPS FORCE A NEW ONE TO BE FORMED ON THE NEXT STEP. ! CODE WILL ATTEMPT TO USE A COPY OF THE LAST JACOBIAN FORMED END IF idoub = l + 1 END IF ! ---------------------------------------------------------------------- ! STORE THE Y ARRAY IN THE MATRIX YHOLD. STORE IN THE Y ARRAY THE ! INFORMATION NECESSARY TO PERFORM AN INTERPOLATION TO FIND THE ! SOLUTION AT THE SPECIFIED OUTPUT POINT IF APPROPRIATE. ! ---------------------------------------------------------------------- 440 yhold(1:n,1:l) = y(1:n,1:l) nstep = nstep + 1 jsinup = jsinup + 1 jsnold = jsnold + 1 jstart = nqused t = told + hused hold = h kfail = 0 newpar = 0 cfail = .false. RETURN 450 finish = .false. t=told rmax=2.0D+0 DO j1=1,l DO i=1,n y(i,j1)=yhold(i,j1) END DO END DO IF(ABS(h) <= hmin*1.00001D+0) THEN ! CORRECTOR CONVERGENCE COULD NOT BE ACHIEVED IF(nstep == 0) THEN kflag=-1 ELSE kflag=-3 END IF ! TO SUPPRESS ERROR MESSAGES AT START AS H MAY ! HAVE BEEN TOO LARGE ON THE FIRST STEP. hold=h finish = .true. END IF rh = red ! TRY AGAIN WITH UPDATED PARTIALS IF(ijus == 0) CALL hchose(rh, h) IF(.NOT.finish) THEN GO TO 40 END IF RETURN ! ------------------- END OF SUBROUTINE STIFF -------------------------- 9000 FORMAT (' CORRECTOR HAS NOT CONVERGED') END SUBROUTINE stiff SUBROUTINE rscale(n,l,rh,y) IMPLICIT NONE INTEGER, INTENT(IN) :: n INTEGER, INTENT(IN) :: l REAL (dp), INTENT(IN) :: rh REAL (dp), INTENT(IN OUT) :: y(:,:) ! y(n,15) ! .. SCALAR ARGUMENTS .. ! .. ! .. ARRAY ARGUMENTS .. ! .. ! .. LOCAL SCALARS .. INTEGER :: i, j, j1 REAL (dp) :: ta, tb, tc, td, te, tf, zz ! .. ! .. LOCAL ARRAYS .. REAL (dp) :: di(8,8) ! SUBROUTINE IS FOR RESCALING THE HISTORY ARRAY AFTER A CHANGE IN ! STEPSIZE ! N ORDER OF THE PROBLEM ! L NUMBER OF TERMS IN THE HISTORY ARRAY TO BE RESCALED ! RH RATIO OF THE STEPSIZE CHANGE (I.E. RH = HNEW/HOLD) ! Y() THE HISTORY ARRAY REAL (dp), PARAMETER :: zero = 0.0_dp ! .. di(2,2) = rh IF (l > 2) THEN ta = rh*rh di(2,3) = rh* (1.0D+0-rh)/2.0D+0 di(3,3) = ta IF (l > 3) THEN tb = ta*rh di(2,4) = rh* ((rh-3.0D+0)*rh+2.0D+0)/6.0D+0 di(3,4) = ta* (1.0D+0-rh) di(4,4) = tb IF (l > 4) THEN tc = tb*rh di(2,5) = - (((rh-6.0D+0)*rh+11.0D+0)*rh-6.0D+0)*rh/ 24.0D+0 di(3,5) = ta* ((7.0D+0*rh-18.0D+0)*rh+11.0D+0)/12.0D+0 di(4,5) = 1.5D+0*tb* (1.0D+0-rh) di(5,5) = tc IF (l > 5) THEN td = tc*rh di(2,6) = ((((rh-10.0D+0)*rh+35.0D+0)*rh-50.0D+0) & *rh+24.0D+0)*rh / 120.0D+0 di(3,6) = - (((3.0D+0*rh-14.0D+0)*rh+21.0D+0)*rh -10.0D+0)*ta/12.0D+0 di(4,6) = ((5.0D+0*rh-12.0D+0)*rh+7.0D+0)*tb/4.0D+0 di(5,6) = 2.0D+0*tc* (1.0D+0-rh) di(6,6) = td IF (l > 6) THEN te = td*rh di(2,7) = -rh* (rh-1.0D+0)* (rh-2.0D+0)* & (rh-3.0D+0)*(rh-4.0D+0)*(rh-5.0D+0)/ 720.0D+0 di(3,7) = ta* ((((62.0D+0*rh-450.0D+0)*rh+ & 1190.0D+0)*rh-1350.0D+0)*rh+548.0D+0) /720.0D+0 di(4,7) = tb* (((-18.0D+0*rh+75.0D+0)*rh & -102.0D+0)*rh+45.0D+0)/24.0D+0 di(5,7) = tc* ((13.0D+0*rh-30.0D+0)*rh+17.0D+0) /6.0D+0 di(6,7) = 2.5D+0*td* (1.0D+0-rh) di(7,7) = te IF (l > 7) THEN tf = te*rh di(2,8) = rh*(rh-1.0D+0)*(rh-2.0D+0)*(rh & -3.0D+0)*(rh-4.0D+0)*(rh-5.0D+0) *(rh-6.0D+0)/5040.0D+0 di(3,8) = ta* ((((((-126.0D+0*rh)+1302.0D+0)*rh- & 5250.0D+0)*rh+10290.0D+0)*rh-9744.0D+0 )*rh+3528.0D+0)/5040.0D+0 di(4,8) = tb* ((((43.0D+0*rh-270.0D+0)*rh+ & 625.0D+0)*rh-630.0D+0)*rh+232.0D+0) /120.0D+0 di(5,8) = tc* (((-10.0D+0*rh+39.0D+0)*rh- & 50.0D+0)*rh+21.0D+0)/6.0D+0 di(6,8) = td* ((20.0D+0*rh-45.0D+0)*rh+25.0D+0 )/6.0D+0 di(7,8) = 3.0D+0*te* (1.0D+0-rh) di(8,8) = tf END IF END IF END IF END IF END IF END IF DO i = 1,n DO j = 2,l zz = zero DO j1 = j,l zz = zz + di(j,j1)*y(i,j1) END DO y(i,j) = zz END DO END DO RETURN END SUBROUTINE rscale SUBROUTINE hchose(rh, h) IMPLICIT NONE REAL (dp), INTENT(IN OUT) :: rh REAL (dp), INTENT(IN) :: h ! Local variables INTEGER :: i, i2 ! FIRST MOVE ALL ELEMENTS DOWN ONE PLACE IF (h /= hstpsz(2,1)) THEN DO i=12,2,-1 i2=i-1 hstpsz(1,i)=hstpsz(1,i2) hstpsz(2,i)=hstpsz(2,i2) END DO ! NOW INSERT VALUE OF H USED BEFORE THIS CALL hstpsz(1,2)=h/hstpsz(2,1) hstpsz(2,1)=h END IF ! NOW DECIDE ON THE NEW CHANGE IF (rh > 1.0D+0) THEN ELSE IF (hstpsz(1,2) <= 1.0D+0) THEN ELSE IF ((rh*h) <= hstpsz(2,2)) THEN ELSE rh=hstpsz(2,2)/h END IF hstpsz(1,1)=rh RETURN END SUBROUTINE hchose END MODULE toms703 PROGRAM runmeb USE common703 USE toms703 IMPLICIT NONE INTERFACE SUBROUTINE fcn(t,y,ydot) USE common703, ONLY: dp IMPLICIT NONE REAL (dp), INTENT(IN) :: t REAL (dp), INTENT(IN) :: y(:) REAL (dp), INTENT(OUT) :: ydot(:) END SUBROUTINE fcn END INTERFACE REAL (dp) :: y(20) ! THE INCLUSION OF THE COMMON BLOCK /MEBDF2/ ALLOWS THE USER TO ! ACCESS VARIOUS COUNTERS THAT MAY BE OF INTEREST TO HIM. ! THESE VARIABLES ARE EXPLAINED IN THE COMMENTS IN SUBROUTINE MEBDF. INTEGER, PARAMETER :: lout = 6 REAL (dp) :: err1, err2, err3, err4, hstart, h0, tol, x, xend, xout INTEGER :: index, n, nn n=4 DO nn=3,10 INDEX = 1 xend=20.0D+0 x=0.0D+0 y(1)=1.0D+0 y(2)=1.0D+0 y(3)=1.0D+0 y(4)=1.0D+0 IF (nn == 3 .OR. nn == 4) hstart = 1.0D-4 IF (nn == 5 .OR. nn == 6) hstart = 1.0D-5 IF (nn == 7 .OR. nn == 8) hstart = 1.0D-6 IF (nn == 9 .OR. nn == 10) hstart = 1.0D-7 ! THESE ARE THE INITIAL STEPS CHOSEN BY THE DETEST PACKAGE FOR THIS PROBLEM h0=hstart tol=10.0D+0**(-nn) xout=20.0D+0 10 CALL mebdf(n, x, h0, y, xout, xend, tol, 21, INDEX, lout, fcn) IF ((INDEX /= 0) .AND. (INDEX /= 3)) THEN IF(INDEX == 1) THEN INDEX=0 GO TO 10 END IF WRITE( 6 ,15) INDEX 15 FORMAT(' ***INTEGRATION HAS FAILED*** WITH INDEX=', i3) STOP END IF ! HAVE COMPLETED ONE STEP ! HAVE WE FINISHED YET? IF( INDEX == 0) THEN ! THEN WE HAVE EFFECTIVELY HIT TOUT x=xout ELSE INDEX=3 GO TO 10 END IF WRITE(6, 690) tol 690 FORMAT(' REQUESTED TOLERANCE ', G10.3) WRITE(6,700) nfe,nje 700 FORMAT(' ', i5, ' FUNCTION EVALUATIONS ', & i5, ' JACOBIAN EVALUATIONS') err1 = y(1) - 4.539992969929191D-05 err2 = y(2) - 2.061153036149920D-09 err3 = y(3) - 2.823006338263857D-56 err4 = y(4) - 5.235792540515384D-52 ! THESE VERY SMALL CONSTANTS SHOULD BE SET TO ZERO IF THERE ! ARE LIKELY TO BE DIFFICULTIES DUE TO UNDERFLOW. WRITE(6,710) 710 FORMAT(' ERRORS AT ENDPOINT') WRITE(6,720) err1, err2, err3, err4 720 FORMAT(' ERRS', 5G14.5) END DO STOP END PROGRAM runmeb SUBROUTINE fcn(t,y,ydot) USE common703, ONLY: dp IMPLICIT NONE REAL (dp), INTENT(IN) :: t REAL (dp), INTENT(IN) :: y(:) REAL (dp), INTENT(OUT) :: ydot(:) ydot(1) = -0.5_dp*y(1) ydot(2) = -y(2) ydot(3) = -100.0_dp*y(3) ydot(4) = -90.0_dp*y(4) ! N.B. This application does not use argument t, but for ELF90: IF (t < -9.9E99_dp) STOP RETURN END SUBROUTINE fcn SUBROUTINE pderv(t,y,pw) USE common703, ONLY: dp IMPLICIT NONE REAL (dp), INTENT(IN) :: t REAL (dp), INTENT(IN) :: y(:,:) REAL (dp), INTENT(OUT) :: pw(:,:) ! pw(4,4) pw(1:4,1:4) = 0.0_dp pw(1,1) = -0.5_dp pw(2,2) = -1.0_dp pw(3,3) = -100.0_dp pw(4,4) = -90.0_dp ! N.B. This application does not use arguments t & y, but for ELF90: IF (t < -9.9E99_dp) STOP IF (y(1,1) < -9.9E99_dp) STOP RETURN END SUBROUTINE pderv