C ALGORITHM 703 , COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 18, NO. 2, JUNE, 1992, PP. 156-158. SUBROUTINE MEBDF(N,T0,H0,Y0,TOUT,TEND,EPS,MF,INDEX,LOUT,LWORK, + WORK,LIWORK,IWORK,FCN) C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C THIS SUBROUTINE IS FOR THE PURPOSE * C OF SPLITTING UP THE WORK ARRAYS WORK AND IWORK * C FOR USE INSIDE THE INTEGRATOR STIFF * C <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< C THIS IS THE JULY 1990 VERSION OF OVDRIV, A PACKAGE FOR C THE SOLUTION OF THE INITIAL VALUE PROBLEM FOR SYSTEMS OF C ORDINARY DIFFERENTIAL EQUATIONS C DY/DT = F(Y,T), Y=(Y(1),Y(2),Y(3), . . . ,Y(N)) C SUBROUTINE OVDRIV IS A DRIVER ROUTINE FOR THIS PACKAGE C C REFERENCES C C 1. J.R. CASH, ON THE INTEGRATION OF STIFF SYSTEMS OF O.D.E.S C USING EXTENDED BACKWARD DIFFERENTIATION FORMULAE, C NUMER. MATH. 34, 235-246 (1980) C 2. THE INTEGRATION OF STIFF INITIAL VALUE PROBLEMS IN O.D.E.S C USING MODIFIED EXTENDED BACKWARD DIFFERENTIATION FORMULAE, C COMP. AND MATHS. WITH APPLICS., 9, 645-657, (1983). C 3. J.R. CASH, STABLE RECURSIONS WITH APPLICATIONS TO THE C NUMERICAL SOLUTION OF STIFF SYSTEMS, ACADEMIC PRESS.(1979) C 4. A.C. HINDMARSH, GEAR.. ORDINARY DIFFERENTIAL EQUATION C SYSTEM SOLVER, UCID-30001 REV. 3, LAWRENCE LIVERMORE C LABORATORY, P.O. BOX 808, LIVERMORE, CA 94550, DEC. 1974. C 5. J.R. CASH AND S. CONSIDINE, AN MEBDF CODE FOR STIFF INITIAL C VALUE PROBLEMS, ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE, C 1992. C C ---------------------------------------------------------------- C OVDRIV IS TO BE CALLED ONCE FOR EACH OUTPUT VALUE OF T, AND C IN TURN MAKES REPEATED CALLS TO THE CORE INTEGRATOR STIFF. C C THE INPUT PARAMETERS ARE .. C N = THE NUMBER OF FIRST ORDER DIFFERENTIAL EQUATIONS. C T0 = THE INITIAL VALUE OF T, THE INDEPENDENT VARIABLE C (USED ONLY ON THE FIRST CALL) C H0 = THE NEXT STEP SIZE IN T (USED FOR INPUT ONLY ON THE C FIRST CALL) C Y0 = A VECTOR OF LENGTH N CONTAINING THE INITIAL VALUES OF Y C (USED FOR INPUT ONLY ON FIRST CALL) C TEND = END OF THE RANGE OF INTEGRATION. C TOUT = THE VALUE OF T AT WHICH OUTPUT IS DESIRED NEXT. C INTEGRATION WILL NORMALLY GO SLIGHTLY BEYOND TOUT C AND THE PACKAGE WILL INTERPOLATE TO T = TOUT C MUST HAVE TEND .GE. TOUT IF H0 IS POSITIVE AND C TOUT .GE. TEND OTHERWISE. C EPS = THE RELATIVE ERROR BOUND (USED ONLY ON THE FIRST C CALL ,UNLESS INDEX = -1). SINGLE STEP ERROR C ESTIMATES DIVIDED BY YMAX(I) WILL BE KEPT LESS THAN C EPS IN ROOT-MEAN-SQUARE NORM (I.E. EUCLIDEAN NORM C DIVIDED BY SQRT(N) ). THE VECTOR YMAX OF WEIGHTS IS C COMPUTED IN OVDRIV. INITIALLY YMAX(I) IS ABS(Y(I)) C WITH A DEFAULT VALUE OF 1 IF Y(I)=0 INITIALLY. C THEREAFTER, YMAX(I) IS THE LARGEST VALUE OF ABS(Y(I)) C SEEN SO FAR, OR THE INITIAL VALUE YMAX(I) IF THAT IS C LARGER C MF = THE METHOD FLAG. AT PRESENT ONLY MF=21 OR 22 IS ALLOWED C WHICH IS EXTENDED BACKWARD DIFFERENTIATION FORMULAE C USING THE CHORD METHOD WITH ANALYTIC OR NUMERICAL C JACOBIAN. THE USER NEEDS TO SPECIFY SUBROUTINE PDERV C IF MF=21 (SEE BELOW) C INDEX = INTEGER USED ON INPUT TO INDICATE TYPE OF CALL, C 1 THIS IS THE FIRST CALL FOR THE PROBLEM C INDEX HAS TO BE SET = 0 AFTER FIRST CALL. C 0 THIS IS NOT THE FIRST CALL FOR THIS PROBLEM C AND INTEGRATION IS TO CONTINUE UP TO T = TOUT C OR BEYOND. C -1 THIS IS NOT THE FIRST CALL FOR THE PROBLEM, C AND THE USER HAS RESET N, EPS, AND/OR MF. C 2 SAME AS 0 EXCEPT THAT TOUT HAS TO BE HIT C EXACTLY (NO INTERPOLATION IS DONE). C ASSUMES TOUT .GE. THE CURRENT T. C 3 ONE STEP MODE. SAME AS 0 EXCEPT CONTROL C RETURNS TO CALLING PROGRAM AFTER ONE STEP. C SINCE THE NORMAL OUTPUT VALUE OF INDEX IS 0, C IT NEED NOT BE RESET FOR NORMAL CONTINUATION. C LOUT = LOGICAL OUTPUT DEVICE. C LWORK = DIMENSION OF REAL WORK ARRAY WORK. LWORK MUST BE AT C LEAST 2*N*N + 37*N C WORK = REAL WORK ARRAY. C LIWORK= DIMENSION OF INTEGER WORK ARRAY. LIWORK MUST BE AT C LEAST N. C IWORK = INTEGER WORK ARRAY. C C AFTER THE INITIAL CALL, IF A NORMAL RETURN OCCURED AND A NORMAL C CONTINUATION IS DESIRED, SIMPLY RESET TOUT AND CALL AGAIN. C ALL OTHER PARAMETERS WILL BE READY FOR THE NEXT CALL. C A CHANGE OF PARAMETERS WITH INDEX = -1 CAN BE MADE AFTER C EITHER A SUCCESSFUL OR AN UNSUCCESSFUL RETURN. C C THE OUTPUT PARAMETERS ARE.. C T0 = THE VALUE OF T WHICH RELATES TO THE CURRENT SOLUTION C POINT Y0() C H0 = THE STEPSIZE H USED LAST, WHETHER SUCCESSFULLY OR NOT. C Y0 = THE COMPUTED VALUES OF Y AT T = TOUT C TOUT = UNCHANGED FROM ITS INPUT VALUE. C INDEX = INTEGER USED ON OUTPUT TO INDICATE RESULTS, WITH C THE FOLLOWING VALUES AND MEANINGS.. C C 3 ONE STEP MODE IS BEING USED AND A SUCCESSFUL STEP HAS C BEEN COMPLETED. THE USER SHOULD CALL MEBDF AGAIN WITH C INDEX = 3. C 1 ONE SUCCESSFUL STEP HAS BEEN COMPLETED AND THE ONE STEP C MODE IS NOT BEING USED. THE USER MUST NOW SET INDEX=0 C AND CALL MEBDF AGAIN. C 0 INTEGRATION WAS COMPLETED TO TOUT OR BEYOND. C C -1 THE INTEGRATION WAS HALTED AFTER FAILING TO PASS THE C ERROR TEST EVEN AFTER REDUCING H BY A FACTOR OF C 1.E10 FROM ITS INITIAL VALUE. C C -2 AFTER SOME INITIAL SUCCESS, THE INTEGRATION WAS C HALTED EITHER BY REPEATED ERROR TEST FAILURES OR BY C A TEST ON EPS. TOO MUCH ACCURACY HAS BEEN REQUESTED. C C -3 THE INTEGRATION WAS HALTED AFTER FAILING TO ACHIEVE C CORRECTOR CONVERGENCE EVEN AFTER REDUCING H BY A C FACTOR OF 1.E10 FROM ITS INITIAL VALUE. C C -4 IMMEDIATE HALT BECAUSE OF ILLEGAL VALUES OF INPUT C PARAMETERS. SEE PRINTED MESSAGE. C C -5 INDEX WAS -1 ON INPUT, BUT THE DESIRED CHANGES OF C PARAMETERS WERE NOT IMPLEMENTED BECAUSE TOUT C WAS NOT BEYOND T. INTERPOLATION AT T = TOUT WAS C PERFORMED AS ON A NORMAL RETURN. TO TRY AGAIN, C SIMPLY CALL AGAIN WITH INDEX = -1 AND A NEW TOUT. C C -10 INSUFFICIENT REAL WORKSPACE FOR THE INTEGRATION C C -11 INSUFFICIENT INTEGER WORKSPACE FOR THE INTEGRATION C C C IN ADDITION TO OVDRIV, THE FOLLOWING ROUTINES ARE PROVIDED C IN THE PACKAGE.. C C INTERP( - ) INTERPOLATES TO GET THE OUTPUT VALUES C AT T=TOUT FROM THE DATA IN THE Y ARRAY. C STIFF( - ) IS THE CORE INTEGRATOR ROUTINE. IT PERFORMS A C SINGLE STEP AND ASSOCIATED ERROR CONTROL. C COSET( - ) SETS COEFFICIENTS FOR BACKWARD DIFFERENTIATION C SCHEMES FOR USE IN THE CORE INTEGRATOR. C PSET( - ) COMPUTES AND PROCESSES THE JACOBIAN C MATRIX J = DF/DY C DEC( - ) PERFORMS AN LU DECOMPOSITION ON A MATRIX. C SOL( - ) SOLVES LINEAR SYSTEMS A*X = B AFTER DEC C HAS BEEN CALLED FOR THE MATRIX A C C C THE FOLLOWING ROUTINES ARE TO BE SUPPLIED BY THE USER.. C C FCN(T,Y,YDOT) COMPUTES THE FUNCTION YDOT = F(Y,T), THE C RIGHT HAND SIDE OF THE O.D.E. C HERE Y AND YDOT ARE VECTORS OF LENGTH N C PDERV(T,Y,PD) COMPUTES THE N*N JACOBIAN MATRIX OF C PARTIAL DERIVATIVES AND STORES IT IN PD C AS AN N BY N ARRAY. PD(I,J) IS TO BE SET C TO THE PARTIAL DERIVATIVE OF YDOT(I) WITH C RESPECT TO Y(J). PDERV IS CALLED ONLY IF C MITER = 1. OTHERWISE A DUMMY ROUTINE CAN C BE SUBSTITUTED. C C THE DIMENSION OF PW BELOW MUST BE AT LEAST N**2. THE DIMENSIONS C OF YMAX,ERROR,SAVE1,SAVE2,IPIV AND THE FIRST DIMENSION C OF Y SHOULD ALL BE AT LEAST N. C C COMMON BLOCKS OF INTEREST TO THE USER. C C /MEBDF2/ C C HUSED LAST STEPSIZE SUCCESSFULLY USED BY THE INTEGRATOR C NQUSED ORDER LAST SUCCESSFULLY USED C NSTEP NUMBER OF STEPS TAKEN SO FAR C NFE NUMBER OF FUNCTION EVALUATIONS USED SO FAR C NJE NUMBER OF JACOBIAN EVALUATIONS USED SO FAR C NDEC NUMBER OF LU DECOMPOSITIONS USED SO FAR C NBSOL NUMBER OF 'BACKSOLVES' USED SO FAR C NPSET NUMBER OF TIMES A NEW COEFFICIENT MATRIX HAS BEEN C FORMED SO FAR C NCOSET NUMBER OF TIMES THE ORDER OF THE METHOD USED HAS BEEN C CHANGED SO FAR C MAXORD THE MAXIMUM ORDER USED SO FAR IN THE INTEGRATION C TO ACCESS THESE THE USER SHOULD PUT THE FOLLOWING COMMON C STATEMENT IN HIS MAIN PROGRAM. C COMMON/MEBDF2/HUSED,NQUSED,NSTEP,NFE,NJE,NDEC,NBSOL, C + NPSET,NCOSET,MAXORD C C ***************************************************************** C THE USER WILL NEED TO SET THE MACHINE PRECISION CONSTANTS FOR C HIS MACHINE. THESE ARE C UROUND - THE UNIT ROUNDOFF ERROR. C EPSJAC - THE SQUARE ROOT OF UROUND. C THESE ARE SET IN A BLOCK DATA STATEMENT FOLLOWING THIS DRIVER. C ****************************************************************** C C .. SCALAR ARGUMENTS .. REAL EPS,H0,T0,TOUT INTEGER INDEX,LIWORK,LOUT,LWORK,MF,N C .. C .. ARRAY ARGUMENTS .. REAL WORK(LWORK),Y0(N) INTEGER IWORK(LIWORK) C .. C .. LOCAL SCALARS .. INTEGER I1,I2,I3,I4,I5,I6,I7,I8,I9 C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL OVDRIV,FCN C .. C .. SAVE STATEMENT .. SAVE I1,I2,I3,I4,I5,I6,I7,I8,I9 C .. IF (INDEX.EQ.1) THEN IF (N.LE.0) THEN WRITE (LOUT,9020) N INDEX = -4 ELSE I1 = N*15 + 1 I2 = I1 + N*15 I3 = I2 + N*2 I4 = I3 + N I5 = I4 + N I6 = I5 + N I7 = I6 + N I8 = I7 + N I9 = I8 + N*N IF (LWORK.LT. (I9+N*N-1)) THEN INDEX = -10 WRITE (LOUT,9000) I9 + N*N - 1 END IF IF (LIWORK.LT.N) THEN INDEX = -11 WRITE (LOUT,9010) N END IF END IF IF (INDEX.LT.0) RETURN END IF CALL OVDRIV(N,T0,H0,Y0,TOUT,TEND,EPS,MF,INDEX,LOUT,WORK,WORK(I1), + WORK(I2),WORK(I3),WORK(I4),WORK(I5),WORK(I6),WORK(I7), + WORK(I8),WORK(I9),IWORK,FCN) C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C WORK() HOUSES THE FOLLOWING ARRAYS C C Y(N,15) , YHOLD(N,15) , YNHOLD(N,2) , YMAX(N) C ERRORS(N) , SAVE1(N) , SAVE2(N) , ARH(N) , PW(N*N) C PWCOPY(N*N) C <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< RETURN 9000 FORMAT (/,/,' ***** ERROR ***** INTEGRATION HALTED IN MEBDF',/,/, + ' >>> REAL WORKSPACE IS INSUFFICIENT <<< ',/, + ' WORKSPACE MUST BE AT LEAST ',I8,' ELEMENTS LONG') 9010 FORMAT (/,/,' ***** ERROR ***** INTEGRATION HALTED IN MEBDF',/,/, + ' >>> INTEGER WORKSPACE IS INSUFFICIENT <<< ',/, + ' WORKSPACE MUST BE AT LEAST ',I6,' ELEMENTS LONG') 9020 FORMAT (/,/,' ***** ERROR ***** INTEGRATION HALTED IN MEBDF',/,/, + ' >>> ILLEGAL VALUE FOR NUMBER OF EQUATIONS <<< ',/, + ' WITH N = ',I6) END BLOCK DATA C C MACHINE DEPENDENT CONSTANTS ARE SET HERE C C UROUND THE SMALLEST +VE NUMBER SUCH THAT C C 1.0E+0 + UROUND .NE. 1.0E+0 C C EPSJAC THIS IS USED IN FORMING A NUMERICAL JACOBIAN AND IS C C EPSJAC = SQRT( UROUND ) C C .. COMMON BLOCKS .. COMMON /MACHNE/UROUND,EPSJAC REAL EPSJAC,UROUND C .. C .. DATA STATEMENTS .. DATA UROUND,EPSJAC/7.1E-15,8.660254E-8/ C .. END SUBROUTINE OVDRIV(N,T0,H0,Y0,TOUT,TEND,EPS,MF,INDEX,LOUT,Y,YHOLD, + YNHOLD,YMAX,ERRORS,SAVE1,SAVE2,ARH,PW,PWCOPY, + IPIV,FCN) C C START OF PROGRAM PROPER C C .. SCALAR ARGUMENTS .. REAL EPS,H0,T0,TOUT INTEGER INDEX,LOUT,MF,N C .. C .. ARRAY ARGUMENTS .. REAL ARH(N),ERRORS(N),PW(*),PWCOPY(*),SAVE1(N),SAVE2(N),Y(N,15), + Y0(N),YHOLD(N,15),YMAX(N),YNHOLD(N,2) INTEGER IPIV(N) C .. C .. LOCAL SCALARS .. REAL AYI,D INTEGER I,KGO,NHCUT C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL INTERP,STIFF,FCN C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS,AMAX1,REAL C .. C .. COMMON BLOCKS .. COMMON /FAILS/NFAIL1,NFAIL2,NT COMMON /MACHNE/UROUND,EPSJAC 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 /MEBDF6/CRATE1,CRATE2,TCRAT1,TCRAT2,AVNEW2,AVOLD2,AVNEWJ, + AVOLDJ COMMON /MEBDFL/CFAIL,JNEWIM,SAMPLE REAL AVNEW2,AVNEWJ,AVOLD2,AVOLDJ,CRATE1,CRATE2,EPSCPY,EPSJAC,H, + HMAX,HMIN,HUSED,T,TCRAT1,TCRAT2,TOUTCP,UROUND INTEGER 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 CFAIL,JNEWIM,SAMPLE C .. IF (INDEX.EQ.0) THEN C I.E. NORMAL CONTINUATION OF INTEGRATION T0=T HMAX = ABS(TEND-T0)*10.0E+0 IF ((T-TOUT)*H.GE.0.0E+0) THEN C 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.EQ.2) THEN C I.E. CONTINUING INTEGRATION BUT WISH TO HIT TOUT T0 = T HMAX = ABS(TEND-T0)*10.0E+0 IF (((T+H)-TOUT)*H.GT.0.0E+0) THEN C WE HAVE ALREADY OVERSHOT THE END OR WE WILL C DO SO ON THE NEXT STEP IF (((T-TOUT)*H.GE.0.0E+0) .OR. (ABS(T-TOUT).LE. + 100.0E+0*UROUND*HMAX)) THEN C HAVE OVERSHOT, SO INTERPOLATE CALL INTERP(N,JSTART,H,T,Y,TOUT,Y0) T0 = TOUT H0 = H INDEX = KFLAG RETURN ELSE C WILL PASS TOUT ON NEXT STEP WITH CURRENT STEPSIZE C SO REDUCE STEPSIZE TO HIT TOUT 'EXACTLY' H = (TOUT-T)* (1.0E+0-4.0E+0*UROUND) JSTART = -1 END IF END IF ELSE IF (INDEX.EQ.-1) THEN C NOT FIRST CALL BUT PARAMETERS RESET T0 = T IF ((T-TOUT)*H.GE.0.0E+0) THEN C 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.EQ.3) THEN T0 = T IF ((T-TOUT)*H.GE.0.0E+0) THEN C HAVE OVERSHOT,SO INTERPOLATE CALL INTERP(N,JSTART,H,T,Y,TOUT,Y0) INDEX = KFLAG T0 = TOUT H0 = H RETURN END IF ELSE C INDEX SHOULD BE 1 AND THIS IS THE FIRST CALL FOR THIS PROBLEM C CHECK THE ARGUMENTS THAT WERE PASSED FOR CORRECTNESS IF (INDEX.NE.1) THEN C VALUE OF INDEX NOT ALLOWED WRITE (LOUT,9070) INDEX INDEX = -4 END IF IF (EPS.LE.0.0E+0) THEN C ILLEGAL VALUE FOR RELATIVE ERROR TOLERENCE WRITE (LOUT,9040) INDEX = -4 END IF IF (N.LE.0) THEN C ILLEGAL VALUE FOR THE NUMBER OF EQUATIONS WRITE (LOUT,9050) INDEX = -4 END IF IF ((T0-TOUT)*H0.GE.0.0E+0) THEN C PARAMETERS FOR INTEGRATION ARE ILLEGAL WRITE (LOUT,9060) INDEX = -4 END IF IF((TEND-TOUT)*H0.LT.-100.0E+0*UROUND) THEN C TOUT IS BEYOND TEND WRITE (LOUT,9110) TOUT,TEND INDEX = -4 ENDIF IF ((MF.NE.21) .AND. (MF.NE.22)) THEN C ILLEGAL VALUE FOR METHOD FLAG WRITE (LOUT,9090) MF INDEX = -4 END IF IF (INDEX.NE.1) THEN RETURN ELSE C THE INITIAL PARAMETERS ARE O.K. SO INITIALISE EVERYTHING C ELSE NECESSARY FOR THE INTEGRATION. C IF VALUES OF YMAX OTHER THAN THOSE SET BELOW ARE DESIRED, C THEY SHOULD BE SET HERE. ALL YMAX(I) MUST BE POSITIVE. IF C VALUES FOR HMIN OR HMAX, THE BOUNDS ON ABS(H), OTHER THAN C THOSE BELOW ARE DESIRED, THEY SHOULD BE SET BELOW DO 10 I = 1,N YMAX(I) = ABS(Y0(I)) IF (YMAX(I).EQ.0.0E+0) YMAX(I) = 1.0E+0 Y(I,1) = Y0(I) 10 CONTINUE NCOPY = N T = T0 H = H0 HMIN = ABS(H0) HMAX = ABS(T0-TEND)*10.0E+0 EPSCPY = EPS MFCOPY = MF JSTART = 0 NHCUT = 0 END IF END IF C <<<<<<<<<<<<<<<<< C < TAKE A STEP > C <<<<<<<<<<<<<<<<< 20 IF ((T+H).EQ.T) THEN WRITE (LOUT,9000) NT = NT + 1 END IF CALL STIFF(EPS,H,HMAX,HMIN,JSTART,KFLAG,MF,T,TEND,Y,N,YMAX, + ERRORS,SAVE1,SAVE2,PW,PWCOPY,YHOLD,YNHOLD,ARH,IPIV,LOUT,FCN) KGO = 1 - KFLAG IF (KGO.EQ.1) THEN C NORMAL RETURN FROM STIFF GO TO 30 ELSE IF (KGO.EQ.2) THEN C COULD NOT ACHIEVE ERROR WITH HMIN C SO CHOP HMIN IF WE HAVEN'T DONE SO 10 TIMES GO TO 60 ELSE IF (KGO.EQ.3) THEN C ERROR SMALLER THAN CAN BE HANDLED FOR THIS PROBLEM WRITE (LOUT,9010) T,H GO TO 70 ELSE IF (KGO.EQ.4) THEN C COULD NOT ACHIEVE CONVERGENCE WITH HMIN WRITE (LOUT,9030) T GO TO 60 END IF 30 CONTINUE C --------------------------------------------------------------------- C NORMAL RETURN FROM THE INTEGRATOR. C C THE WEIGHTS YMAX(I) ARE UPDATED. IF DIFFERENT VALUES ARE DESIRED, C THEY SHOULD BE SET HERE. A TEST IS MADE FOR EPS BEING TOO SMALL C FOR THE MACHINE PRECISION. C C ANY OTHER TESTS OR CALCULATIONS THAT ARE REQUIRED AFTER EVERY C STEP SHOULD BE INSERTED HERE. C C IF INDEX = 3, Y0 IS SET TO THE CURRENT Y VALUES ON RETURN. C IF INDEX = 2, H IS CONTROLLED TO HIT TOUT (WITHIN ROUNDOFF C ERROR), AND THEN THE CURRENT Y VALUES ARE PUT IN Y0 ON RETURN. C FOR ANY OTHER VALUE OF INDEX, CONTROL RETURNS TO THE INTEGRATOR C UNLESS TOUT HAS BEEN REACHED. THEN INTERPOLATED VALUES OF Y ARE C COMPUTED AND STORED IN Y0 ON RETURN. C IF INTERPOLATION IS NOT DESIRED, THE CALL TO INTERP SHOULD BE C REMOVED AND CONTROL TRANSFERRED TO STATEMENT 500 INSTEAD OF 520. C -------------------------------------------------------------------- D = 0.0E+0 DO 40 I = 1,N AYI = ABS(Y(I,1)) YMAX(I) = AMAX1(YMAX(I),AYI) D = D + (AYI/YMAX(I))**2 40 CONTINUE D = D* (UROUND/EPS)**2 IF (D.GT.REAL(N)) THEN C EPS SMALLER THAN MACHINE PRECISION WRITE (LOUT,9020) T KFLAG = -2 GO TO 70 END IF IF (INDEX.EQ.3.OR.INDEX.EQ.1) GO TO 70 IF (ABS(T-TOUT).LE.ABS(15.0E+0*UROUND*TOUT)) THEN C EFFECTIVELY WE HAVE HIT TOUT INDEX = KFLAG T0 = TOUT DO 50 I = 1,N Y0(I) = Y(I,1) 50 CONTINUE H0 = H RETURN END IF IF (INDEX.EQ.2) THEN C CONTINUING INTEGRATION BUT MUST HIT TOUT EXACTLY IF (((T+H)-TOUT)*H.GT.0.0E+0) THEN C WE HAVE ALREADY OVERSHOT THE END OR WE WILL DO C SO ON THE NEXT STEP IF (((T-TOUT)*H.GE.0.0E+0) .OR. (ABS(T-TOUT).LE. + 100.0E+0*UROUND*HMAX)) THEN C HAVE OVERSHOT, SO INTERPOLATE CALL INTERP(N,JSTART,H,T,Y,TOUT,Y0) T0 = TOUT H0 = H INDEX = KFLAG RETURN ELSE C WILL PASS TOUT ON NEXT STEP WITH CURRENT STEPSIZE C SO REDUCE STEPSIZE TO HIT TOUT 'EXACTLY' H = (TOUT-T)* (1.0E+0-4.0E+0*UROUND) JSTART = -1 END IF END IF ELSE IF ((T-TOUT)*H.GE.0.0E+0) THEN C 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 C ------------------------------------------------------------------- C ON AN ERROR RETURN FROM THE INTEGRATOR, AN IMMEDIATE RETURN OCCURS C IF KFLAG = -2, AND RECOVERY ATTEMPTS ARE MADE OTHERWISE. C H AND HMIN ARE REDUCED BY A FACTOR OF .1 UP TO 10 TIMES C BEFORE GIVING UP. C -------------------------------------------------------------------- 60 CONTINUE IF (NHCUT.EQ.10) THEN C HAVE REDUCED H TEN TIMES WRITE (LOUT,9100) GO TO 70 END IF NHCUT = NHCUT + 1 HMIN = 0.1E+0*HMIN H = 0.1E+0*H IF (NSTEP.GT.0) THEN JSTART = -1 ELSE C THE INITIAL STEPSIZE WAS TOO LARGE. RESET CERTAIN VARIABLES C AND TRY AGAIN. NJE = 0 NFE = 1 CFAIL = .TRUE. NEWPAR = 0 MQ1TMP = 0 MQ2TMP = 0 MEQC1 = 0 MEQC2 = 0 TCRAT1 = 0.0E+0 TCRAT2 = 0.0E+0 CRATE1 = 1.0E+0 CRATE2 = 1.0E+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).GT.1000.0E+0*UROUND) THEN DO 80 I = 1,N Y0(I) = Y(I,1) 80 CONTINUE T0 = T ELSE C HAVE EITHER PASSED TOUT OR WE ARE EXTREMELY CLOSE TO IT C SO INTERPOLATE. CALL INTERP(N,JSTART,H,T,Y,TOUT,Y0) T0 = TOUT INDEX = KFLAG END IF H0 = H IF(KFLAG.NE.0) INDEX = KFLAG RETURN C -------------------------- 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 .GE. 0.',/,/) 9070 FORMAT (/,/,' ILLEGAL INPUT.. INDEX =',I5,/,/) 9080 FORMAT (/,/,' INDEX = -1 ON INPUT WITH (T-TOUT)*H .GE. 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 AND TEND SO THAT (TEND-TOUT)*H0 .GE.0.0') END SUBROUTINE INTERP(N,JSTART,H,T,Y,TOUT,Y0) C .. SCALAR ARGUMENTS .. REAL H,T,TOUT INTEGER JSTART,N C .. C .. ARRAY ARGUMENTS .. REAL Y(N,15),Y0(N) C .. C .. LOCAL SCALARS .. REAL S,S1 INTEGER I,J,L C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC REAL C .. DO 10 I = 1,N Y0(I) = Y(I,1) 10 CONTINUE L = JSTART + 2 S = (TOUT-T)/H S1 = 1.0E+0 DO 30 J = 2,L S1 = S1* (S+REAL(J-2))/REAL(J-1) DO 20 I = 1,N Y0(I) = Y0(I) + S1*Y(I,J) 20 CONTINUE 30 CONTINUE RETURN C -------------- END OF SUBROUTINE INTERP --------------------------- END SUBROUTINE COSET(NQ,EL,ELST,TQ,MAXDER) C -------------------------------------------------------------------- C COSET IS CALLED BY THE INTEGRATOR AND SETS THE COEFFICIENTS USED C BY THE CONVENTIONAL BACKWARD DIFFERENTIATION SCHEME. THE VECTOR C EL OF LENGTH NQ+1 DETERMINES THE BASIC METHOD. THE VECTOR ELST OF C DIMENSION NQ+2 DETERMINES THE MODIFIED EXTENDED BACKWARD C DIFFERENTIATION FORMULA COEFFICIENTS. THE VECTOR TQ OF C LENGTH 4 IS INVOLVED IN ADJUSTING THE STEPSIZE IN RELATION TO THE C TRUNCATION ERROR. ITS VALUES ARE GIVEN BY THE PERTST ARRAY. THE C VECTORS EL AND TQ BOTH DEPEND ON METH AND NQ. COSET ALSO SETS C MAXDER, THE MAXIMUM ORDER OF THE METHOD AVAILABLE. THE C COEFFICIENTS IN PERTST NEED TO BE GIVEN TO ONLY ABOUT ONE PERCENT C ACCURACY. THE ORDER IN WHICH THE GROUPS APPEAR BELOW IS C COEFFICIENTS FOR ORDER NQ-1, COEFFICIENTS FOR ORDER NQ, C COEFFICIENTS FOR ORDER NQ+1. C ------------------------------------------------------------------- C .. SCALAR ARGUMENTS .. INTEGER MAXDER,NQ C .. C .. ARRAY ARGUMENTS .. REAL EL(10),ELST(10),TQ(4) C .. C .. LOCAL SCALARS .. INTEGER K C .. C .. LOCAL ARRAYS .. REAL PERTST(8,3) C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC REAL C .. C .. COMMON BLOCKS .. COMMON /MEBDF2/HUSED,NQUSED,NSTEP,NFE,NJE,NDEC,NBSOL,NPSET,NCOSET, + MAXORD REAL HUSED INTEGER MAXORD,NBSOL,NCOSET,NDEC,NFE,NJE,NPSET,NQUSED,NSTEP C .. C .. DATA STATEMENTS .. DATA PERTST(1,1)/1./,PERTST(2,1)/2./,PERTST(3,1)/4.5/, + PERTST(4,1)/7.333/,PERTST(5,1)/10.42/,PERTST(6,1)/13.7/, + PERTST(7,1)/17.15/,PERTST(8,1)/20.74/ DATA PERTST(1,2)/2./,PERTST(2,2)/4.5/,PERTST(3,2)/7.333/, + PERTST(4,2)/10.42/,PERTST(5,2)/13.7/,PERTST(6,2)/17.15/, + PERTST(7,2)/20.74/,PERTST(8,2)/24.46/ DATA PERTST(1,3)/4.5/,PERTST(2,3)/7.333/,PERTST(3,3)/10.42/, + PERTST(4,3)/13.7/,PERTST(5,3)/17.15/,PERTST(6,3)/20.74/, + PERTST(7,3)/24.46/,PERTST(8,3)/1./ C .. C ------------------------------------------------------------------- C THE FOLLOWING COEFFICIENTS SHOULD BE DEFINED TO MACHINE ACCURACY. C THEIR DERIVATION IS GIVEN IN REFERENCE 2. C ------------------------------------------------------------------- IF (NQ.GT.MAXORD) MAXORD = NQ MAXDER = 6 NCOSET = NCOSET + 1 GO TO (10,20,30,40,50,60,70) NQ 10 EL(1) = 1.0E+0 ELST(1) = 1.5E+0 ELST(3) = -0.5E+0 GO TO 80 20 EL(1) = 6.6666666666667E-01 EL(3) = 3.3333333333333E-01 ELST(1) = 9.5652173913043E-01 ELST(3) = 2.1739130434782E-01 ELST(4) = -1.7391304347826E-01 GO TO 80 30 EL(1) = 5.4545454545455E-01 EL(3) = 4.5454545454545E-01 EL(4) = 1.8181818181818E-01 ELST(1) = 7.6142131979695E-01 ELST(3) = 3.2994923857868E-01 ELST(4) = 8.6294416243654E-02 ELST(5) = -9.1370558375634E-02 GO TO 80 40 EL(1) = 0.48E+0 EL(3) = 0.52E+0 EL(4) = 0.28E+0 EL(5) = 0.12E+0 ELST(1) = 6.5733706517393E-01 ELST(3) = 4.0023990403838E-01 ELST(4) = 1.5793682526989E-01 ELST(5) = 4.4382247101159E-02 ELST(6) = -5.7576969212315E-02 GO TO 80 50 EL(1) = 4.3795620437956E-01 EL(3) = 5.62043795620436E-01 EL(4) = 3.43065693430656E-01 EL(5) = 1.97080291970802E-01 EL(6) = 8.75912408759123E-02 ELST(1) = 5.9119243917152E-01 ELST(3) = 4.4902473356122E-01 ELST(4) = 2.1375427307460E-01 ELST(5) = 9.0421610027481503E-02 ELST(6) = 2.6409276761177E-02 ELST(7) = -4.0217172732757E-02 GO TO 80 60 EL(1) = 4.08163265306120E-01 EL(3) = 5.91836734693874E-01 EL(4) = 3.87755102040813E-01 EL(5) = 2.51700680272107E-01 EL(6) = 1.49659863945577E-01 EL(7) = 6.80272108843534E-02 ELST(1) = 5.4475876041119E-01 ELST(3) = 4.8525549636077E-01 ELST(4) = 2.5789750131312E-01 ELST(5) = 1.3133738525800E-01 ELST(6) = 5.7677396763462E-02 ELST(7) = 1.7258197643881E-02 ELST(8) = -3.0014256771967E-02 GO TO 80 70 EL(1) = 3.85674931129476E-01 EL(3) = 6.14325068870521E-01 EL(4) = 4.21487603305783E-01 EL(5) = 2.9292929292929E-01 EL(6) = 1.96510560146923E-01 EL(7) = 1.19375573921028E-01 EL(8) = 5.50964187327820E-02 ELST(1) = 5.0999746293734E-01 ELST(3) = 5.1345839935281E-01 ELST(4) = 2.9364346131937E-01 ELST(5) = 1.6664672120553E-01 ELST(6) = 8.8013735242353E-02 ELST(7) = 3.9571794884069E-02 ELST(8) = 1.2039080338722E-02 ELST(9) = -2.3455862290154E-02 80 DO 90 K = 1,3 TQ(K) = PERTST(NQ,K) 90 CONTINUE TQ(4) = 0.5E+0*TQ(2)/REAL(NQ) RETURN C --------------------- END OF SUBROUTINE COSET --------------------- END SUBROUTINE PSET(Y,N,H,T,UROUND,EPSJAC,CON,MITER,IER,NRENEW,YMAX, + SAVE2,PW,PWCOPY,WRKSPC,IPIV,FCN) C ------------------------------------------------------------------- C PSET IS CALLED BY STIFF TO COMPUTE AND PROCESS THE COEFFICIENT C MATRIX I - H*EL(1)*J WHERE J IS AN APPROXIMATION TO C THE RELEVANT JACOBIAN. THIS MATRIX IS THEN SUBJECTED TO LU C DECOMPOSITION IN PREPARATION FOR LATER SOLUTION OF LINEAR SYSTEMS C OF ALGEBRAIC EQUATIONS WITH LU AS THE COEFFICIENT MATRIX. THE C MATRIX J IS FOUND BY THE USER-SUPPLIED ROUTINE PDERV IF MITER=1 C OR BY FINITE DIFFERENCING IF MITER = 2. C IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION WITH C PSET USES THE FOLLOWING .. C EPSJAC = SQRT(UROUND), USED IN NUMERICAL JACOBIAN INCREMENTS. C ******************************************************************* C THE ARGUMENT NRENEW IS USED TO SIGNAL WHETHER OR NOT C WE REQUIRE A NEW JACOBIAN TO BE CALCULATED. C C IF NRENEW > 0 THEN WE REQUIRE A NEW J TO BE COMPUTED C = 0 THEN USE A COPY OF THE LAST J COMPUTED C ******************************************************************* C .. SCALAR ARGUMENTS .. REAL CON,EPSJAC,H,T,UROUND INTEGER IER,MITER,N,NRENEW C .. C .. ARRAY ARGUMENTS .. REAL PW(*),PWCOPY(*),SAVE2(N),WRKSPC(N),Y(N,15),YMAX(N) INTEGER IPIV(N) C .. C .. LOCAL SCALARS .. REAL D,R,R0,TEMPRY,YJ INTEGER I,J,J1,JJKK C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL DEC,FCN,PDERV C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS,AMAX1,SQRT C .. C .. COMMON BLOCKS .. COMMON /MEBDF2/HUSED,NQUSED,NSTEP,NFE,NJE,NDEC,NBSOL,NPSET,NCOSET, + MAXORD REAL HUSED INTEGER MAXORD,NBSOL,NCOSET,NDEC,NFCN1,NFE,NJE,NPSET, + NQUSED,NSTEP C .. NPSET = NPSET + 1 IF (NRENEW.EQ.0) THEN DO 10 I = 1,N*N PW(I) = PWCOPY(I)*CON 10 CONTINUE GO TO 70 END IF IF (MITER.EQ.2) GO TO 30 CALL PDERV(T,Y,PW) DO 20 I = 1,N*N PWCOPY(I) = PW(I) PW(I) = PW(I)*CON 20 CONTINUE NJE = NJE + 1 GO TO 70 30 NJE = NJE + 1 D = 0.0E+0 DO 40 I = 1,N D = D + SAVE2(I)**2 40 CONTINUE R0 = ABS(H)*SQRT(D)*1.0E+03*UROUND J1 = 0 DO 60 J = 1,N YJ = Y(J,1) R = EPSJAC*YMAX(J) R = AMAX1(R,R0) Y(J,1) = Y(J,1) + R D = CON/R CALL FCN(T,Y,WRKSPC) DO 50 I = 1,N JJKK = I + J1 TEMPRY = (WRKSPC(I)-SAVE2(I)) PWCOPY(JJKK) = TEMPRY/R PW(JJKK) = TEMPRY*D 50 CONTINUE Y(J,1) = YJ J1 = J1 + N 60 CONTINUE NFE = NFE + N 70 J = 1 DO 80 I = 1,N PW(J) = PW(J) + 1.0E+0 J = J + (N+1) 80 CONTINUE CALL DEC(N,N,PW,IPIV,IER) NDEC = NDEC + 1 RETURN C ---------------------- END OF SUBROUTINE PSET --------------------- END SUBROUTINE DEC(N,NDIM,A,IP,IER) C ------------------------------------------------------------------- C MATRIX TRIANGULARISATION BY GAUSSIAN ELIMINATION C INPUT.. C N = ORDER OF MATRIX. C NDIM = DECLARED DIMENSION OF ARRAY A. C A = MATRIX TO BE TRIANGULARISED. C OUTPUT.. C A(I,J), I.LE.J = UPPER TRIANGULAR FACTOR, U. C A(I,J), I.GT.J = MULTIPLIERS = LOWER TRIANGULAR FACTOR, I-L. C IP(K), K.LT.N = INDEX OF KTH PIVOT ROW. C IP(N) = (-1)**(NUMBER OF INTERCHANGES) OR 0. C IER = 0 IF MATRIX IS NON-SINGULAR, OR K IF FOUND TO BE SINGULAR C AT STAGE K. C USE SOL TO OBTAIN SOLUTION OF LINEAR SYSTEM. C DETERM(A) = IP(N)*A(1,1)*A(2,2)* . . . *A(N,N). C IF IP(N) = 0, A IS SINGULAR, SOL WILL DIVIDE BY ZERO. C C REFERENCE. C C.B. MOLER, ALGORITHM 423, LINEAR EQUATION SOLVER, C.A.C.M C 15 (1972), P.274. C ------------------------------------------------------------------ C .. SCALAR ARGUMENTS .. INTEGER IER,N,NDIM C .. C .. ARRAY ARGUMENTS .. REAL A(NDIM,N) INTEGER IP(N) C .. C .. LOCAL SCALARS .. REAL T INTEGER I,J,K,KP1,M,NM1 C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS C .. C .. COMMON BLOCKS .. IER = 0 IP(N) = 1 IF (N.EQ.1) GO TO 70 NM1 = N - 1 DO 60 K = 1,NM1 KP1 = K + 1 M = K DO 10 I = KP1,N IF (ABS(A(I,K)).GT.ABS(A(M,K))) M = I 10 CONTINUE IP(K) = M T = A(M,K) IF (M.EQ.K) GO TO 20 IP(N) = -IP(N) A(M,K) = A(K,K) A(K,K) = T 20 IF (T.EQ.0.0E+0) GO TO 80 T = 1.0E+0/T DO 30 I = KP1,N A(I,K) = -A(I,K)*T 30 CONTINUE DO 50 J = KP1,N T = A(M,J) A(M,J) = A(K,J) A(K,J) = T IF (T.EQ.0.0E+0) GO TO 50 DO 40 I = KP1,N A(I,J) = A(I,J) + A(I,K)*T 40 CONTINUE 50 CONTINUE 60 CONTINUE 70 K = N IF (A(N,N).EQ.0.0E+0) GO TO 80 RETURN 80 IER = K IP(N) = 0 RETURN C --------------------- END OF SUBROUTINE DEC ---------------------- END SUBROUTINE SOL(N,NDIM,A,B,IP) C .. SCALAR ARGUMENTS .. INTEGER N,NDIM C .. C .. ARRAY ARGUMENTS .. REAL A(NDIM,N),B(N) INTEGER IP(N) C .. C .. LOCAL SCALARS .. REAL T INTEGER I,K,KB,KM1,KP1,M,NM1 C .. C .. COMMON BLOCKS .. C ------------------------------------------------------------------ C SOLUTION OF LINEAR SYSTEM, A*X = B. C INPUT .. C N = ORDER OF MATRIX. C NDIM = DECLARED DIMENSION OF MATRIX A. C A = TRIANGULARISED MATRIX OBTAINED FROM DEC. C B = RIGHT HAND SIDE VECTOR. C IP = PIVOT VECTOR OBTAINED FROM DEC. C DO NOT USE IF DEC HAS SET IER .NE. 0 C OUTPUT.. C B = SOLUTION VECTOR, X. C ------------------------------------------------------------------ IF (N.EQ.1) GO TO 50 NM1 = N - 1 DO 20 K = 1,NM1 KP1 = K + 1 M = IP(K) T = B(M) B(M) = B(K) B(K) = T DO 10 I = KP1,N B(I) = B(I) + A(I,K)*T 10 CONTINUE 20 CONTINUE DO 40 KB = 1,NM1 KM1 = N - KB K = KM1 + 1 B(K) = B(K)/A(K,K) T = -B(K) DO 30 I = 1,KM1 B(I) = B(I) + A(I,K)*T 30 CONTINUE 40 CONTINUE 50 B(1) = B(1)/A(1,1) RETURN C ------------------------- END OF SUBROUTINE SOL ------------------ END SUBROUTINE ERRORS(N,TQ,EPS,EDN,E,EUP,BND) C *************************************************** C C THIS ROUTINE CALCULATES ERRORS USED IN TESTS C IN STIFF . C C *************************************************** C .. SCALAR ARGUMENTS .. REAL BND,E,EDN,EPS,EUP INTEGER N C .. C .. ARRAY ARGUMENTS .. REAL TQ(4) C .. C .. LOCAL SCALARS .. REAL SQHOL C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC REAL C .. SQHOL = REAL(N)*EPS*EPS EDN = TQ(1)*TQ(1)*SQHOL C C ** ERROR ASSOCIATED WITH LOWER ORDER METHOD C E = TQ(2)*TQ(2)*SQHOL C C ** ERROR ASSOCIATED WITH PRESENT ORDER C EUP = TQ(3)*TQ(3)*SQHOL C C ** ERROR ASSOCIATED WITH HIGHER ORDER METHOD C BND = TQ(4)*TQ(4)*SQHOL*0.3E+0 RETURN END SUBROUTINE PRDICT(T,H,Y,L,N,YPRIME,FCN) C ********************************************************************** C PREDICTS A VALUE FOR Y AT (T+H) GIVEN THE HISTORY ARRAY AT Y(T) C THEN EVALUATES THE DERIVATIVE AT THIS POINT, THE RESULT OF THIS C EVALUATION BEING STORED IN YPRIME() C ********************************************************************** C .. SCALAR ARGUMENTS .. REAL H,T INTEGER L,N C .. C .. ARRAY ARGUMENTS .. REAL Y(N,15),YPRIME(N) C .. C .. LOCAL SCALARS .. INTEGER I,J2 C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL FCN C .. C .. COMMON BLOCKS .. COMMON /MEBDF2/HUSED,NQUSED,NSTEP,NFE,NJE,NDEC,NBSOL,NPSET,NCOSET, + MAXORD REAL HUSED INTEGER MAXORD,NBSOL,NCOSET,NDEC,NFE,NJE,NPSET,NQUSED,NSTEP C .. DO 20 J2 = 2,L DO 10 I = 1,N Y(I,1) = Y(I,1) + Y(I,J2) 10 CONTINUE 20 CONTINUE T = T + H CALL FCN(T,Y,YPRIME) NFE = NFE + 1 RETURN END SUBROUTINE ITRAT2(Y,N,T,HBETA,ERRBND,ARH,CRATE,TCRATE,M,WORKED, + YMAX,ERROR,SAVE1,SAVE2,PW,IPIV,LMB,FCN) C .. SCALAR ARGUMENTS .. REAL CRATE,ERRBND,HBETA,T,TCRATE INTEGER M,N LOGICAL WORKED C .. C .. ARRAY ARGUMENTS .. REAL ARH(N),ERROR(N),PW(*),SAVE1(N),SAVE2(N),Y(N,15),YMAX(N) INTEGER IPIV(N) C .. C .. LOCAL SCALARS .. REAL D,D1,ZERO INTEGER I C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL FCN,SOL C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC AMAX1,AMIN1 C .. C .. COMMON BLOCKS .. COMMON /MEBDF2/HUSED,NQUSED,NSTEP,NFE,NJE,NDEC,NBSOL,NPSET,NCOSET, + MAXORD REAL HUSED INTEGER MAXORD,NBSOL,NCOSET,NDEC,NFE,NJE,NPSET,NQUSED,NSTEP C .. C .. DATA STATEMENTS .. DATA ZERO/0.0E+0/ C .. IF(LMB.EQ.1) GOTO 25 DO 10 I = 1,N SAVE1(I) = -SAVE1(I) + HBETA*SAVE2(I) + ARH(I) 10 CONTINUE CALL SOL(N,N,PW,SAVE1,IPIV) NBSOL = NBSOL + 1 D = ZERO DO 20 I = 1,N ERROR(I) = ERROR(I) + SAVE1(I) D = D + (SAVE1(I)/YMAX(I))**2 SAVE1(I) = Y(I,1) + ERROR(I) 20 CONTINUE TCRATE = TCRATE + CRATE D1 = D M = 1 CALL FCN(T,SAVE1,SAVE2) NFE = NFE + 1 25 CONTINUE WORKED = .TRUE. 30 DO 40 I = 1,N SAVE1(I) = -SAVE1(I) + HBETA*SAVE2(I) + ARH(I) 40 CONTINUE C C IF WE ARE HERE THEN PARTIALS ARE O.K. C CALL SOL(N,N,PW,SAVE1,IPIV) C C WE NOW CALCULATE A WEIGHTED RMS TYPE NORM C NBSOL = NBSOL + 1 D = ZERO DO 50 I = 1,N ERROR(I) = ERROR(I) + SAVE1(I) D = D + (SAVE1(I)/YMAX(I))**2 SAVE1(I) = Y(I,1) + ERROR(I) 50 CONTINUE C ------------------------------------------------------------------- C TEST FOR CONVERGENCE. IF M.GT.0 , AN ESTIMATE OF THE CONVERGENCE C RATE CONSTANT IS STORED IN CRATE, AND THIS IS USED IN THE TEST. C ------------------------------------------------------------------- IF (M.NE.0) THEN IF (D1.NE.ZERO) CRATE = AMAX1(0.9E+0*CRATE,D/D1) END IF TCRATE = TCRATE + CRATE IF ((D*AMIN1(1.0E+0,2.0E+0*CRATE)).LT.ERRBND) RETURN IF (M.NE.0) THEN IF (D.GT.D1) THEN WORKED = .FALSE. RETURN END IF END IF D1 = D IF (M.EQ.3) THEN WORKED = .FALSE. RETURN END IF M = M + 1 CALL FCN(T,SAVE1,SAVE2) NFE = NFE + 1 GO TO 30 END SUBROUTINE STIFF(EPS,H,HMAX,HMIN,JSTART,KFLAG,MF,T,TEND,Y, + N,YMAX,ERROR,SAVE1,SAVE2,PW,PWCOPY,YHOLD,YNHOLD,ARH, + IPIV,LOUT,FCN) C ------------------------------------------------------------------ C STIFF PERFORMS ONE STEP OF THE INTEGRATION OF AN INITIAL VALUE C PROBLEM FOR A SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS. C COMMUNICATION WITH STIFF IS DONE WITH THE FOLLOWING VARIABLES.. C Y AN N BY LMAX+4 ARRAY CONTAINING THE DEPENDENT VARIABLES C AND THEIR BACKWARD DIFFERENCES. LMAX-1 =MAXDER IS THE C MAXIMUM ORDER AVAILABLE. SEE SUBROUTINE COSET. C Y(I,J+1) CONTAINS THE JTH BACKWARD DIFFERENCE OF Y(I) C T THE INDEPENDENT VARIABLE. T IS UPDATED ON EACH STEP TAKEN. C H THE STEPSIZE TO BE ATTEMPTED ON THE NEXT STEP. C H IS ALTERED BY THE ERROR CONTROL ALGORITHM DURING C THE PROBLEM. H CAN BE EITHER POSITIVE OR NEGATIVE BUT C ITS SIGN MUST REMAIN CONSTANT THROUGHOUT THE PROBLEM. C HMIN THE MINIMUM AND MAXIMUM ABSOLUTE VALUE OF THE STEPSIZE C HMAX TO BE USED FOR THE STEP. THESE MAY BE CHANGED AT ANY C TIME BUT WILL NOT TAKE EFFECT UNTIL THE NEXT H CHANGE. C EPS THE RELATIVE ERROR BOUND. SEE DESCRIPTION IN MEBDF C UROUND THE UNIT ROUNDOFF OF THE MACHINE. C N THE NUMBER OF FIRST ORDER DIFFERENTIAL EQUATIONS. C MF THE METHOD FLAG. MUST BE SET TO 21 OR 22 AT PRESENT C KFLAG A COMPLETION CODE WITH THE FOLLOWING MEANINGS.. C 0 THE STEP WAS SUCCESSFUL C -1 THE REQUESTED ERROR COULD NOT BE ACHIEVED C WITH ABS(H) = HMIN. C -2 THE REQUESTED ERROR IS SMALLER THAN CAN C BE HANDLED FOR THIS PROBLEM. C -3 CORRECTOR CONVERGENCE COULD NOT BE C ACHIEVED FOR ABS(H)=HMIN. C ON A RETURN WITH KFLAG NEGATIVE, THE VALUES OF T AND C THE Y ARRAY ARE AS AT THE BEGINNING OF THE LAST C STEP ATTEMPTED, AND H IS THE LAST STEP SIZE ATTEMPTED. C JSTART AN INTEGER USED ON INPUT AND OUTPUT. C ON INPUT IT HAS THE FOLLOWING VALUES AND MEANINGS.. C 0 PERFORM THE FIRST STEP. C .GT.0 TAKE A NEW STEP CONTINUING FROM THE LAST C .LT.0 TAKE THE NEXT STEP WITH A NEW VALUE OF C H, EPS OR N C ON EXIT, JSTART IS NQUSED, THE ORDER OF THE METHOD LAST USED. C YMAX AN ARRAY OF N ELEMENTS WITH WHICH THE ESTIMATED LOCAL C ERRORS IN Y ARE COMPARED C ERROR AN ARRAY OF N ELEMENTS. C SAVE1,2 TWO ARRAYS OF WORKING SPACE BOTH OF LENGTH N. C PW A BLOCK OF LOCATIONS USED FOR PARTIAL DERIVATIVES C IPIV AN INTEGER ARRAY OF LENGTH N USED FOR PIVOT INFORMATION. C C JNEWIM IS TO INDICATE IF PRESENT ITERATION MATRIX C WAS FORMED USING A NEW J OR OLD J. C JSNOLD KEEPS TRACK OF NUMBER OF STEPS TAKEN WITH C PRESENT ITERATION MATRIX (BE IT FORMED BY C A NEW J OR NOT). C AVNEWJ STORES VALUE FOR AVERAGE CRATE WHEN ITERATION C MATRIX WAS FORMED BY A NEW J. C AVOLDJ STORES VALUE FOR AVERAGE CRATE WHEN ITERATION C MATRIX WAS FORMED BY AN OLD J. C NRENEW FLAG THAT IS USED IN COMMUNICATION WITH SUBROUTINE PSET. C IF NRENEW > 0 THEN FORM A NEW JACOBIAN BEFORE C COMPUTING THE COEFFICIENT MATRIX FOR C THE NEWTON ITERATION C = 0 FORM THE COEFFICIENT MATRIX USING A C COPY OF AN OLD JACOBIAN C NEWPAR FLAG USED IN THIS SUBROUTINE TO INDICATE IF A JACOBIAN C HAS BEEN EVALUATED FOR THE CURRENT STEP C ********************************************************************** C .. SCALAR ARGUMENTS .. REAL EPS,H,HMAX,HMIN,T INTEGER JSTART,KFLAG,LOUT,MF,N C .. C .. ARRAY ARGUMENTS .. REAL HSTPSZ(2,14) REAL ARH(N),ERROR(N),PW(*),PWCOPY(*),SAVE1(N),SAVE2(N),Y(N,15), + YHOLD(N,15),YMAX(N),YNHOLD(N,2) INTEGER IPIV(N) C .. C .. LOCAL SCALARS .. REAL BND,D,DDOWN,DSTEP2,DUP,EBDF1,EBDF2,EDN,EFAIL,ENQ1,ENQ2, + ENQ3,EPSOLD,EUP,HOLD,OLDLO,ONE,PR1,PR2,PR3,PRBDF1,PRBDF2, + PRFAIL,QQ,RHBDF1,RHBDF2,RHLMIT,TRANGE,ZERO INTEGER I,IER,IITER,IITER2,ITST,J,J1,J2,KFAIL,LL,LMAX, + M3STEP,MAXDER,METH,MFOLD,MITER,NEWQ LOGICAL FINISH,OVRIDE,WORKED C .. C .. LOCAL ARRAYS .. REAL EL(10),ELST(10),TQ(4) C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL COSET,CPYARY,ERRORS,FCN,HCHOSE,ITRAT2, + PRDICT,PSET,RSCALE,SOL C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS,AMAX1,AMIN1,INT,REAL C .. C .. COMMON BLOCKS .. COMMON/STPSZE/HSTPSZ,TOP COMMON /FAILS/NFAIL1,NFAIL2,NT COMMON /MACHNE/UROUND,EPSJAC 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 REAL AVNEW2,AVNEWJ,AVOLD2,AVOLDJ,CRATE1,CRATE2,E,EPSJAC,HUSED,QI, + RC,RH,RMAX,TCRAT1,TCRAT2,TOLD,UROUND INTEGER IDOUB,ISAMP,IWEVAL,JCHANG,JSINUP,JSNOLD,L,M1,M2,MAXORD, + MEQC1,MEQC2,MQ1TMP,MQ2TMP,NBSOL,NCOSET,NDEC,NEWPAR,NFAIL1, + NFAIL2,NFE,NJE,NPSET,NQ,NQUSED,NRENEW,NSTEP,NT LOGICAL CFAIL,JNEWIM,SAMPLE C .. C .. SAVE STATEMENT .. SAVE EPSOLD,MFOLD,HOLD,LMAX,EDN,EUP,BND,EL,TQ,ELST,E C .. C .. DATA STATEMENTS .. DATA EL(2),ELST(2),OLDLO/3*1.0E+0/ DATA ZERO,ONE/0.0E+0,1.0E+0/ C .. TOLD = T KFLAG = 0 IF (JSTART.GT.0) GO TO 60 IF (JSTART.NE.0) GO TO 30 C ------------------------------------------------------------------ C ON THE FIRST CALL, THE ORDER IS SET TO 1 AND THE INITIAL YDOT C IS CALCULATED. RMAX IS THE MAXIMUM RATIO BY WHICH H CAN BE C INCREASED IN A SINGLE STEP. RMAX IS SET EQUAL TO 1.E4 INITIALLY C TO COMPENSATE FOR THE SMALL INITIAL H, BUT THEN IS NORMALLY = 10. C IF A FAILURE OCCURS (IN CORRECTOR CONVERGENCE OR ERROR TEST), C RMAX IS SET AT 2. FOR THE NEXT INCREASE. C ------------------------------------------------------------------ CALL FCN(T,Y,SAVE1) DO 10 I = 1,N Y(I,2) = H*SAVE1(I) 10 CONTINUE METH = 2 MITER = MF - 10*METH NQ = 1 L = 2 IDOUB = 3 KFAIL = 0 RMAX = 10000.0E+0 RC = ZERO CRATE1 = ONE CRATE2 = ONE JSNOLD = 0 JNEWIM = .TRUE. TCRAT1 = ZERO TCRAT2 = ZERO TOP=0 ATOL=EPS/10.0E+0 DO 15 I=1,12 HSTPSZ(1,I)=1.0E+0 HSTPSZ(2,I)=ATOL 15 CONTINUE 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 C ************************************************** C CFAIL=.TRUE. ENSURES THAT WE CALCULATE A NEW C J ON THE FIRST CALL. C ************************************************** MEQC1 = 0 MEQC2 = 0 MQ1TMP = 0 MQ2TMP = 0 NBSOL = 0 C ----------------------------------------------------------------- C IF THE CALLER HAS CHANGED N OR EPS, THE CONSTANTS E, EDN, EUP C AND BND MUST BE RESET. E IS A COMPARISON FOR ERRORS AT THE C CURRENT ORDER NQ. EUP IS TO TEST FOR INCREASING THE ORDER, C EDN FOR DECREASING THE ORDER. BND IS USED TO TEST FOR CONVERGENCE C OF THE CORRECTOR ITERATES. IF THE CALLER HAS CHANGED H, Y MUST C BE RE-SCALED. IF H IS CHANGED, IDOUB IS SET TO L+1 TO PREVENT C FURTHER CHANGES IN H FOR THAT MANY STEPS. C ----------------------------------------------------------------- CALL COSET(NQ,EL,ELST,TQ,MAXDER) LMAX = MAXDER + 1 RC = RC*EL(1)/OLDLO OLDLO = EL(1) IWEVAL = MITER NRENEW = 1 NEWPAR = 0 C ***************************************************** C NRENEW AND NEWPAR ARE TO INSTRUCT ROUTINE THAT C WE WISH A NEW J TO BE CALCULATED FOR THIS STEP. C ***************************************************** CALL ERRORS(N,TQ,EPS,EDN,E,EUP,BND) EPSOLD = EPS DO 20 I = 1,N ARH(I) = EL(2)*Y(I,1) 20 CONTINUE CALL CPYARY(N*L,Y,YHOLD) QI = H*EL(1) QQ = ONE/QI CALL PRDICT(T,H,Y,L,N,SAVE2,FCN) GO TO 110 C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C DIFFERENT PARAMETERS ON THIS CALL < C <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 30 CALL CPYARY(N*L,YHOLD,Y) IF (MF.NE.MFOLD) THEN METH = MF/10 MITER = MF - 10*METH MFOLD = MF IWEVAL = MITER END IF IF (EPS.NE.EPSOLD) THEN CALL ERRORS(N,TQ,EPS,EDN,E,EUP,BND) EPSOLD = EPS END IF IF (H.NE.HOLD) THEN RH = H/HOLD H = HOLD GO TO 50 ELSE GO TO 60 END IF C ********************************************* C RE-SCALE Y AFTER A CHANGE OF STEPSIZE * C ********************************************* 40 RH = AMAX1(RH,HMIN/ABS(H)) 50 RH = AMIN1(RH,HMAX/ABS(H),RMAX) CALL RSCALE(N,L,RH,Y) RMAX = 10.0E+0 JCHANG = 1 H = H*RH RC = RC*RH IF (JSNOLD.GT.40) THEN CFAIL = .TRUE. NEWPAR = 0 RC = ZERO C ********************************************************************** C CFAIL=TRUE AND NEWPAR=0 SHOULD FORCE A NEW J TO BE EVALUATED C AFTER 42 STEPS WITH AN OLD J, IF WE HAVE HAD A FAILURE OF ANY C KIND ON THE FIRST, SECOND OR THIRD STAGE OF THE CURRENT STEP C ********************************************************************** END IF IDOUB = L + 1 CALL CPYARY(N*L,Y,YHOLD) 60 IF (ABS(RC-ONE).GT.0.3E+0) IWEVAL = MITER C ------------------------------------------------------------------ C THIS SECTION COMPUTES THE PREDICTED VALUES OF Y C AND THE RHS, ARH, FOR USE IN THE NEWTON ITERATION SCHEME. C RC IS THE RATIO OF THE NEW TO OLD VALUES OF THE COEFFICIENT C H*EL(1). WHEN RC DIFFERS FROM 1 BY MORE THAN 30 PERCENT, IWEVAL IS C SET TO MITER TO FORCE THE PARTIALS TO BE UPDATED. C ------------------------------------------------------------------ QI = H*EL(1) QQ = ONE/QI DO 70 I = 1,N ARH(I) = EL(2)*Y(I,1) 70 CONTINUE DO 90 J1 = 2,NQ DO 80 I = 1,N ARH(I) = ARH(I) + EL(J1+1)*Y(I,J1) 80 CONTINUE 90 CONTINUE IF (JCHANG.EQ.1) THEN C IF WE HAVE CHANGED STEPIZE THEN PREDICT A VALUE FOR Y(T+H) C AND EVALUATE THE DERIVATIVE THERE (STORED IN SAVE2()) CALL PRDICT(T,H,Y,L,N,SAVE2,FCN) ELSE C ELSE USE THE VALUES COMPUTED FOR THE SECOND BDF FROM THE LAST C STEP. Y( ,LMAX+3) HOLDS THE VALUE FOR THE DERIVATIVE AT (T+H) C AND Y( ,LMAX+4) HOLDS THE APPROXIMATION TO Y AT THIS POINT. DO 100 I = 1,N SAVE2(I) = Y(I,LMAX+3) Y(I,1) = Y(I,LMAX+4) 100 CONTINUE T = T + H END IF 110 IF (IWEVAL.LE.0) GO TO 120 C ------------------------------------------------------------------- C IF INDICATED, THE MATRIX P = I - H*EL(1)*J IS RE-EVALUATED BEFORE C STARTING THE CORRECTOR ITERATION. IWEVAL IS SET = 0 TO INDICATE C THAT THIS HAS BEEN DONE. P IS COMPUTED AND PROCESSED IN PSET. C THE PROCESSED MATRIX IS STORED IN PW C ------------------------------------------------------------------- IWEVAL = 0 RC = ONE IITER = MEQC1 - MQ1TMP IITER2 = MEQC2 - MQ2TMP IF (JNEWIM) THEN IF (JSNOLD.GE.3) THEN AVNEWJ = TCRAT1/REAL(IITER) AVNEW2 = TCRAT2/REAL(IITER2) ELSE AVNEWJ = ONE AVNEW2 = ONE END IF ELSE C C MATRIX P WAS FORMED WITH A COPY OF J C IF (JSNOLD.GE.3) THEN AVOLDJ = TCRAT1/REAL(IITER) AVOLD2 = TCRAT2/REAL(IITER2) IF (AVOLDJ.LT.AVNEWJ) THEN AVNEWJ = AVOLDJ ELSE IF (((ABS(AVOLDJ-AVNEWJ)).GT.0.3E+0) .OR. + ((AVOLDJ.GT.0.85E+0).AND. (AVOLDJ.NE.ONE))) THEN C C SINCE IN CERTAIN INSTANCES AVOLDJ WILL C BE 1.0 AND THERE WILL BE NO NEED TO C UPDATE J. C CFAIL = .TRUE. CRATE1 = ONE CRATE2 = ONE END IF ELSE CFAIL = .TRUE. CRATE1 = ONE CRATE2 = ONE C C ********************************************************* C IF WE HAVE REACHED HERE THINGS MUST HAVE GONE WRONG C ********************************************************* C END IF END IF TCRAT1 = ZERO TCRAT2 = ZERO IF (CFAIL) THEN IF (NEWPAR.EQ.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) C NOTE THAT ERROR() IS JUST BEING USED AS A WORKSPACE BY PSET IF (IER.NE.0) THEN C IF IER>0 THEN WE HAVE A SINGULARITY IN THE COEFFICIENT MATRIX IJUS=1 RED=0.5E+0 GO TO 450 END IF 120 DO 130 I = 1,N SAVE1(I) = Y(I,1) ERROR(I) = ZERO 130 CONTINUE M1 = 0 C ********************************************************************** C UP TO 4 CORRECTOR ITERATIONS ARE TAKEN. A CONVERGENCE TEST IS C MADE ON THE R.M.S. NORM OF EACH CORRECTION ,USING BND, WHICH C DEPENDS ON EPS. THE SUM OF THE CORRECTIONS IS ACCUMULATED IN THE C VECTOR ERROR(I). THE Y ARRAY IS NOT ALTERED IN THE CORRECTOR C LOOP. THE UPDATED Y VECTOR IS STORED TEMPORARILY IN SAVE1. C ********************************************************************** 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 C C NOW TEST TO SEE IF IT WAS SUCCESSFUL OR NOT C C IF (.NOT.WORKED) THEN NFAIL1 = NFAIL1 + 1 C ********************************************************************** C THE CORRECTOR ITERATION FAILED TO CONVERGE IN 4 TRIES. IF C PARTIALS ARE NOT UP TO DATE, THEY ARE RE-EVALUATED FOR THE C NEXT TRY. OTHERWISE THE Y ARRAY IS RETRACTED TO ITS VALUES C BEFORE PREDICTION AND H IS REDUCED IF POSSIBLE. IF NOT A C NON-CONVERGENCE EXIT IS TAKEN C ********************************************************************** IF (IWEVAL.EQ.-1) THEN C HAVE BEEN USING OLD PARTIALS, UPDATE THEM AND TRY AGAIN IWEVAL = MITER CFAIL = .TRUE. CALL FCN(T,Y,SAVE2) NFE = NFE + 1 GO TO 110 END IF IJUS=0 RED=0.5E+0 GO TO 450 END IF IWEVAL = -1 HUSED = H NQUSED = NQ DO 140 I = 1,N SAVE2(I) = (SAVE1(I)-ARH(I))*QQ Y(I,1) = SAVE1(I) 140 CONTINUE C C UPDATE THE DIFFERENCES AT N+1 C DO 160 J = 2,L DO 150 I = 1,N Y(I,J) = Y(I,J-1) - YHOLD(I,J-1) 150 CONTINUE 160 CONTINUE C C COMPUTE ERROR IN THE SOLUTION C D = ZERO DO 170 I = 1,N D = D + ((Y(I,L)-YHOLD(I,L))/YMAX(I))**2 170 CONTINUE C C STORING Y FROM FIRST STEP FOR USE IN THIRD STEP C DO 180 I = 1,N YNHOLD(I,1) = Y(I,1) YNHOLD(I,2) = SAVE2(I) 180 CONTINUE IF (D.GT.E) GO TO 330 IF ((M1+M2).GE.ITST) THEN M2 = 0 EBDF1 = 0.5E+0/REAL(L) PRBDF1 = ((D/ (E*0.7E+0))**EBDF1)*1.5E+0 + 1.6E-6 RHBDF1 = ONE/PRBDF1 CALL HCHOSE(RHBDF1,H,OVRIDE) RHLMIT = 1 - (5+ITST-M1-M2)*0.5E-1 IF (RHBDF1.LT.RHLMIT) THEN IJUS=1 RED=RHBDF1 GOTO 450 END IF END IF KFAIL = 0 C ---------------------------------------------------------------------- DO 190 I = 1,N ARH(I) = EL(2)*Y(I,1) 190 CONTINUE DO 210 J1 = 2,NQ DO 200 I = 1,N ARH(I) = ARH(I) + EL(J1+1)*Y(I,J1) 200 CONTINUE 210 CONTINUE CALL PRDICT(T,H,Y,L,N,SAVE2,FCN) DO 220 I = 1,N SAVE1(I) = Y(I,1) ERROR(I) = ZERO 220 CONTINUE M2 = 0 C C FOR NOW WILL ASSUME THAT WE DO NOT WISH TO SAMPLE C AT THE N+2 STEP POINT C CALL ITRAT2(Y,N,T,QI,BND,ARH,CRATE2,TCRAT2,M2,WORKED,YMAX,ERROR, + SAVE1,SAVE2,PW,IPIV,1,FCN) MEQC2 = MEQC2 + M2 + 1 C C NOW CHECK TO SEE IF IT WAS SUCCESSFUL OR NOT C IF (.NOT.WORKED) THEN NFAIL2 = NFAIL2 + 1 IJUS=0 RED=0.5E+0 GOTO 450 END IF C C IF WE ARE DOWN TO HERE THEN THINGS MUST HAVE CONVERGED C DO 230 I = 1,N Y(I,LMAX+3) = (SAVE1(I)-ARH(I))*QQ Y(I,LMAX+4) = SAVE1(I) 230 CONTINUE IF ((M2.GE.2) .AND. (CRATE2.LT.0.35E+0)) THEN DO 250 J = 2,NQ DO 240 I = 1,N SAVE1(I) = SAVE1(I) - Y(I,J) 240 CONTINUE 250 CONTINUE DSTEP2 = ZERO DO 260 I = 1,N DSTEP2 = DSTEP2 + ((SAVE1(I)-YNHOLD(I,1)-Y(I,L))/YMAX(I))**2 260 CONTINUE EBDF2 = 0.5E+0/REAL(L) PRBDF2 = ((DSTEP2/ (E*0.7E+0))**EBDF2)*1.5E+0 + 1.6E-6 RHBDF2 = ONE/PRBDF2 CALL HCHOSE(RHBDF2,H,OVRIDE) IF (RHBDF2.LT.0.9E+0) THEN IJUS=1 RED=RHBDF2 GOTO 450 END IF END IF C C WE ARE NOW COMPUTING THE THIRD STAGE C LL = L + 1 T = TOLD + H DO 280 I = 1,N ARH(I) = H* (ELST(NQ+2)*Y(I,LMAX+3)+ + (ELST(1)-EL(1))*YNHOLD(I,2)) DO 270 J1 = 1,NQ ARH(I) = ARH(I) + ELST(J1+1)*YHOLD(I,J1) 270 CONTINUE 280 CONTINUE DO 290 I = 1,N SAVE2(I) = YNHOLD(I,2) Y(I,1) = YNHOLD(I,1) 290 CONTINUE M3STEP = 0 300 DO 310 I = 1,N SAVE1(I) = -Y(I,1) + QI*SAVE2(I) + ARH(I) 310 CONTINUE CALL SOL(N,N,PW,SAVE1,IPIV) NBSOL = NBSOL + 1 D = ZERO DO 320 I = 1,N D = D + (SAVE1(I)/YMAX(I))**2 Y(I,1) = Y(I,1) + SAVE1(I) 320 CONTINUE IF ((D*AMIN1(ONE,2.0E+0*CRATE1)).LE.BND) GO TO 360 IF (M3STEP.EQ.4) THEN WRITE (LOUT,9000) IJUS=1 RED=0.5E+0 GO TO 450 END IF M3STEP = M3STEP + 1 CALL FCN(T,Y,SAVE2) NFE = NFE + 1 GO TO 300 330 KFAIL = KFAIL - 1 C ********************************************************************** C THE ERROR TEST FAILED. KFAIL KEEPS TRACK OF MULTIPLE FAILURES. C RESTORE T AND THE Y ARRAY TO THEIR PREVIOUS VALUES AND PREPARE TO C TRY THE STEP AGAIN. COMPUTE THE OPTIMAL STEP SIZE FOR THIS ORDER C AND ONE ORDER LOWER. C ********************************************************************** T = TOLD HOLD = H EFAIL = 0.5E+0/REAL(L) CALL CPYARY(N*L,YHOLD,Y) RMAX = 2.0E+0 IF (ABS(H).LE.HMIN*1.00001E+0) THEN C C REQUESTED ERROR NOT POSSIBLE WITH GIVEN HMIN C KFLAG = -1 HOLD = H RETURN END IF IF (KFAIL.LE.-3) GO TO 340 C C PREDICTING A NEW H AFTER INSUFFICIENT ACCURACY C PRFAIL = ((D/E)**EFAIL)*1.5E+0 + 1.6E-6 NEWQ = NQ RH = ONE/ (PRFAIL*REAL(-KFAIL)) CALL HCHOSE(RH,H,OVRIDE) GO TO 40 C ********************************************************************** C CONTROL REACHES THIS STAGE IF 3 OR MORE FAILURES HAVE OCCURED. C IT IS ASSUMED THAT THE DERIVATIVES THAT HAVE ACCUMULATED IN THE Y C ARRAY HAVE ERRORS OF THE WRONG ORDER. HENCE THE FIRST DERIVATIVE C IS RE-COMPUTED, AND THE ORDER IS SET TO 1. THEN H IS REDUCED BY A C FACTOR OF 10, AND THE STEP IS RETRIED. AFTER A TOTAL OF 7 C FAILURES AN EXIT IS TAKEN WITH KFLAG=-2. C ********************************************************************** 340 IF (KFAIL.EQ.-7) THEN C ERROR SMALLER THAN CAN BE HANDLED FOR PROBLEM KFLAG = -2 HOLD = H RETURN END IF C ********************************* C START FROM ORDER 1 AGAIN * C ********************************* JCHANG = 1 RH = AMAX1(HMIN/ABS(H),0.1E+0) CALL HCHOSE(RH,H,OVRIDE) H = H*RH CALL FCN(T,YHOLD,SAVE1) NFE = NFE + 1 DO 350 I = 1,N Y(I,1) = YHOLD(I,1) Y(I,2) = H*SAVE1(I) YHOLD(I,2) = Y(I,2) 350 CONTINUE IWEVAL = MITER CFAIL = .TRUE. C SINCE WE HAVE HAD PROBLEMS PROCEED WITH THIS ORDER C FOR 10 STEPS (IF WE CAN) IDOUB = 10 IF (NQ.EQ.1) GO TO 60 NQ = 1 L = 2 C 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) C NOW JUMP TO NORMAL CONTINUATION POINT GO TO 60 C ********************************************************************** C THE ITERATION FOR THE CORRECTED SOLUTION HAS CONVERGED. C UPDATE THE Y ARRAY. C ********************************************************************** 360 CONTINUE DO 380 J2 = 2,LL DO 370 I = 1,N Y(I,J2) = Y(I,J2-1) - YHOLD(I,J2-1) 370 CONTINUE 380 CONTINUE C ---------------------------------------------------------------------- C IF THE COUNTER IDOUB EQUALS 2 , STORE Y(I,LMAX+5) WHICH IS USED IN C ASSESSING THE POSSIBILITY OF INCREASING THE ORDER. IF IDOUB = 0 C CONTROL PASSES TO 480 WHERE AN ATTEMPT TO CHANGE THE STEPSIZE AND C ORDER IS MADE. C ---------------------------------------------------------------------- IF (IDOUB.EQ.2) THEN DO 390 I = 1,N Y(I,LMAX+5) = Y(I,LL) 390 CONTINUE END IF IDOUB = IDOUB - 1 TRANGE=(TEND-TOLD-H)*H IF(TRANGE.LT.0.0E+0) GOTO 440 JCHANG = 0 IF (IDOUB.EQ.0) THEN SAMPLE = .FALSE. ISAMP = ISAMP + 1 IF (ISAMP.EQ.4) THEN SAMPLE = .TRUE. ISAMP = 0 END IF C ********************************************************************** C NOW COMPUTE THE FACTORS PR1, PR2 AND PR3, BY WHICH C H COULD BE DIVIDED AT ORDER NQ-1, ORDER NQ AND ORDER NQ+1 C RESPECTIVELY. THE SMALLEST OF THESE IS DETERMINED AND THE NEW C ORDER CHOSEN ACCORDINGLY. IF THE ORDER IS TO BE INCREASED WE C MUST COMPUTE ONE MORE BACKWARD DIFFERENCE. C ********************************************************************** PR3 = 1.E+20 IF (L.NE.LMAX) THEN DUP = ZERO DO 400 I = 1,N DUP = DUP + ((Y(I,LL)-Y(I,LMAX+5))/YMAX(I))**2 400 CONTINUE ENQ3 = 0.5E+0/REAL(L+1) PR3 = ((DUP/EUP)**ENQ3)*1.7E+0 + 1.8E-6 END IF ENQ2 = 0.5E+0/REAL(L) D = ZERO DO 410 I = 1,N D = D + ((Y(I,L)-YHOLD(I,L))/YMAX(I))**2 410 CONTINUE PR2 = ((D/E)**ENQ2)*1.5E+0 + 1.6E-6 PR1 = 1.E+20 IF (NQ.GT.1) THEN DDOWN = ZERO DO 420 I = 1,N DDOWN = DDOWN + (Y(I,L)/YMAX(I))**2 420 CONTINUE ENQ1 = 0.5E+0/REAL(NQ) PR1 = ((DDOWN/EDN)**ENQ1)*1.6E+0 + 1.7E-6 END IF IF (PR2.LE.PR3) THEN IF (PR2.GT.PR1) THEN NEWQ = NQ - 1 RH = 1.0E+0/PR1 ELSE NEWQ = NQ RH = 1.0E+0/PR2 END IF ELSE IF (PR3.LT.PR1) THEN NEWQ = L RH = 1.0E+0/PR3 ELSE NEWQ = NQ - 1 RH = 1.0E+0/PR1 END IF RH = AMIN1(RH,RMAX) CALL HCHOSE(RH,H,OVRIDE) IF ((JSINUP.LE.20).AND.(KFLAG.EQ.0).AND.(RH.LT.1.1E+0)) THEN C WE HAVE RUN INTO PROBLEMS IDOUB = 10 NQ = NQUSED L = NQ + 1 GO TO 440 END IF C ********************************************************************** C IF THERE IS A CHANGE IN ORDER, RESET NQ, L AND THE C COEFFICIENTS. IN ANY CASE H IS RESET ACCORDINGLY TO RH AND THE C Y ARRAY IS RE-SCALED C ********************************************************************** IF (NEWQ.NE.NQ) THEN IF (NEWQ.GT.NQ) THEN C ADD AN EXTRA TERM TO THE HISTORY ARRAY DO 430 I = 1,N Y(I,LL) = Y(I,L) - YHOLD(I,L) 430 CONTINUE END IF NQ = NEWQ L = NQ + 1 C 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 C NOW RESCALE THE HISTORY ARRAY FOR THE NEW STEPSIZE RH = AMAX1(RH,HMIN/ABS(H)) RH = AMIN1(RH,HMAX/ABS(H),RMAX) CALL RSCALE(N,L,RH,Y) RMAX = 10.0E+0 JCHANG = 1 H = H*RH RC = RC*RH IF (JSNOLD.GT.40) THEN RC = ZERO C IF WE HAVE BEEN USING THE SAME COEFFICIENT MATRIX FOR MORE C THAN 40 STEPS FORCE A NEW ONE TO BE FORMED ON THE NEXT STEP. C CODE WILL ATTEMPT TO USE A COPY OF THE LAST JACOBIAN FORMED END IF IDOUB = L + 1 END IF 440 CONTINUE C ---------------------------------------------------------------------- C STORE THE Y ARRAY IN THE MATRIX YHOLD. STORE IN THE Y ARRAY THE C INFORMATION NECASSARY TO PERFORM AN INTERPOLATION TO FIND THE C SOLUTION AT THE SPECIFIED OUTPUT POINT IF APPROPRIATE. C ---------------------------------------------------------------------- CALL CPYARY(N*L,Y,YHOLD) NSTEP = NSTEP + 1 JSINUP = JSINUP + 1 JSNOLD = JSNOLD + 1 JSTART = NQUSED T = TOLD + HUSED HOLD = H KFAIL = 0 NEWPAR = 0 CFAIL = .FALSE. RETURN 450 CONTINUE FINISH = .FALSE. T=TOLD RMAX=2.0E+0 DO 460 J1=1,L DO 460 I=1,N Y(I,J1)=YHOLD(I,J1) 460 CONTINUE IF(ABS(H).LE.HMIN*1.00001E+0) THEN C C CORRECTOR CONVERGENCE COULD NOT BE ACHIEVED C IF(NSTEP.EQ.0) THEN KFLAG=-1 ELSE KFLAG=-3 END IF C C TO SUPPRESS ERROR MESSAGES AT START AS H MAY C HAVE BEEN TOO LARGE ON THE FIRST STEP. C HOLD=H FINISH = .TRUE. END IF RH = RED C C TRY AGAIN WITH UPDATED PARTIALS C IF(IJUS.EQ.0) CALL HCHOSE(RH,H,OVRIDE) IF(.NOT.FINISH) THEN GO TO 40 ELSE RETURN END IF C ------------------- END OF SUBROUTINE STIFF -------------------------- 9000 FORMAT (1X,' CORRECTOR HAS NOT CONVERGED') END SUBROUTINE RSCALE(N,L,RH,Y) C .. SCALAR ARGUMENTS .. REAL RH INTEGER L,N C .. C .. ARRAY ARGUMENTS .. REAL Y(N,15) C .. C .. LOCAL SCALARS .. REAL TA,TB,TC,TD,TE,TF,ZERO,ZZ INTEGER I,J,J1 C .. C .. LOCAL ARRAYS .. REAL DI(8,8) C .. C .. DATA STATEMENTS .. C SUBROUTINE IS FOR RESCALING THE HISTORY ARRAY AFTER A CHANGE IN C STEPSIZE C C N ORDER OF THE PROBLEM C L NUMBER OF TERMS IN THE HISTORY ARRAY TO BE RESCALED C RH RATIO OF THE STEPSIZE CHANGE (I.E. RH = HNEW/HOLD) C Y() THE HISTORY ARRAY C DATA ZERO/0.0E+0/ C .. DI(2,2) = RH IF (L.GT.2) THEN TA = RH*RH DI(2,3) = RH* (1.0E+0-RH)/2.0E+0 DI(3,3) = TA IF (L.GT.3) THEN TB = TA*RH DI(2,4) = RH* ((RH-3.0E+0)*RH+2.0E+0)/6.0E+0 DI(3,4) = TA* (1.0E+0-RH) DI(4,4) = TB IF (L.GT.4) THEN TC = TB*RH DI(2,5) = - (((RH-6.0E+0)*RH+11.0E+0)*RH-6.0E+0)*RH/ + 24.0E+0 DI(3,5) = TA* ((7.0E+0*RH-18.0E+0)*RH+11.0E+0)/12.0E+0 DI(4,5) = 1.5E+0*TB* (1.0E+0-RH) DI(5,5) = TC IF (L.GT.5) THEN TD = TC*RH DI(2,6) = ((((RH-10.0E+0)*RH+35.0E+0)*RH-50.0E+0) + *RH+24.0E+0)*RH/120.0E+0 DI(3,6) = - (((3.0E+0*RH-14.0E+0)*RH+21.0E+0)*RH + -10.0E+0)*TA/12.0E+0 DI(4,6) = ((5.0E+0*RH-12.0E+0)*RH+7.0E+0)*TB/4.0E+0 DI(5,6) = 2.0E+0*TC* (1.0E+0-RH) DI(6,6) = TD IF (L.GT.6) THEN TE = TD*RH DI(2,7) = -RH* (RH-1.0E+0)* (RH-2.0E+0)* + (RH-3.0E+0)*(RH-4.0E+0)*(RH-5.0E+0)/ + 720.0E+0 DI(3,7) = TA* ((((62.0E+0*RH-450.0E+0)*RH+ + 1190.0E+0)*RH-1350.0E+0)*RH+548.0E+0) + /720.0E+0 DI(4,7) = TB* (((-18.0E+0*RH+75.0E+0)*RH + -102.0E+0)*RH+45.0E+0)/24.0E+0 DI(5,7) = TC* ((13.0E+0*RH-30.0E+0)*RH+17.0E+0) + /6.0E+0 DI(6,7) = 2.5E+0*TD* (1.0E+0-RH) DI(7,7) = TE IF (L.GT.7) THEN TF = TE*RH DI(2,8) = RH*(RH-1.0E+0)*(RH-2.0E+0)*(RH + -3.0E+0)*(RH-4.0E+0)*(RH-5.0E+0) + *(RH-6.0E+0)/5040.0E+0 DI(3,8) = TA* ((((((-126.0E+0*RH)+1302.0E+0)*RH- + 5250.0E+0)*RH+10290.0E+0)*RH-9744.0E+0 + )*RH+3528.0E+0)/5040.0E+0 DI(4,8) = TB* ((((43.0E+0*RH-270.0E+0)*RH+ + 625.0E+0)*RH-630.0E+0)*RH+232.0E+0) + /120.0E+0 DI(5,8) = TC* (((-10.0E+0*RH+39.0E+0)*RH- + 50.0E+0)*RH+21.0E+0)/6.0E+0 DI(6,8) = TD*((20.0E+0*RH-45.0E+0)*RH+25.0E+0 + )/6.0E+0 DI(7,8) = 3.0E+0*TE* (1.0E+0-RH) DI(8,8) = TF END IF END IF END IF END IF END IF END IF DO 30 I = 1,N DO 20 J = 2,L ZZ = ZERO DO 10 J1 = J,L ZZ = ZZ + DI(J,J1)*Y(I,J1) 10 CONTINUE Y(I,J) = ZZ 20 CONTINUE 30 CONTINUE RETURN END SUBROUTINE CPYARY(NELEM,SOURCE,TARGET) C C COPIES THE ARRAY SOURCE() INTO THE ARRAY TARGET() C C THIS SUBROUTINE COULD BE REPLACED BY THE BLAS ROUTINE SCOPY C (AFTER CHANGING THE ARGUMENT LIST APPROPRIATELY) C C .. SCALAR ARGUMENTS .. INTEGER NELEM C .. C .. ARRAY ARGUMENTS .. REAL SOURCE(NELEM),TARGET(NELEM) C .. C .. LOCAL SCALARS .. INTEGER I C .. DO 10 I = 1,NELEM TARGET(I) = SOURCE(I) 10 CONTINUE RETURN END SUBROUTINE HCHOSE(RH,H,OVRIDE) COMMON / STPSZE / HSTPSZ(2,14),TOP LOGICAL OVRIDE C C FIRST MOVE ALL ELEMENTS DOWN ONE PLACE C IF (H.NE.HSTPSZ(2,1)) THEN DO 5 I=12,2,-1 I2=I-1 HSTPSZ(1,I)=HSTPSZ(1,I2) 5 HSTPSZ(2,I)=HSTPSZ(2,I2) C C NOW INSERT VALUE OF H USED BEFORE THIS CALL C HSTPSZ(1,2)=H/HSTPSZ(2,1) HSTPSZ(2,1)=H END IF C C NOW DECIDE ON THE NEW CHANGE C IF (RH.GT.1.0E+0) THEN OVRIDE=.FALSE. ELSE IF (HSTPSZ(1,2).LE.1.0E+0) THEN OVRIDE=.FALSE. ELSE IF ((RH*H).LE.HSTPSZ(2,2)) THEN OVRIDE=.FALSE. ELSE RH=HSTPSZ(2,2)/H OVRIDE=.TRUE. END IF HSTPSZ(1,1)=RH RETURN C C ************************************************************ C END PROGRAM RUNMEB(INPUT,OUTPUT,RESULT,TAPE5=INPUT,TAPE6=RESULT) DIMENSION Y(20) C C THE REAL WORK ARRAY NEEDS TO BE AT LEAST 2*N*N + 37*N WORDS LONG C THE INTEGER WORK ARRAY NEEDS TO BE AT LEAST N WORDS LONG C DIMENSION WORK(4680) INTEGER IWORK(40) EXTERNAL FCN COMMON/MEBDF2/HUSED,NQUSED,NSTEP,NFE,NJE,NDEC,NBSOL,NPSET, + NCOSET,MAXORD C C THE INCLUSION OF THE COMMON BLOCK MEBDF2 ALLOWS THE USER TO C ACCESS VARIOUS COUNTERS WHICH MAY BE OF INTEREST TO HIM. C THESE VARIABLES ARE EXPLAINED IN THE COMMENTS IN SUBROUTINE C MEBDF. C THE USER NEEDS TO SET THE UNIT ROUND-OFF ERROR FOR THE MACHINE C BEING USED IN THE BLOCK DATA STATEMENT FOLLOWING SUBROUTINE MEBDF C DATA LOUT / 6 / LWORK = 4680 LIWORK = 40 N=4 DO 100 NN=3,10 INDEX = 1 XEND=20.0E+0 X=0.0E+0 Y(1)=1.0E+0 Y(2)=1.0E+0 Y(3)=1.0E+0 Y(4)=1.0E+0 IF (NN.EQ. 3 .OR. NN .EQ. 4) HSTART = 1.0E-4 IF (NN.EQ. 5 .OR. NN .EQ. 6) HSTART = 1.0E-5 IF (NN.EQ. 7 .OR. NN .EQ. 8) HSTART = 1.0E-6 IF (NN.EQ. 9 .OR. NN .EQ. 10) HSTART = 1.0E-7 C C THESE ARE THE INITIAL STEPS CHOSEN BY THE DETEST PACKAGE C FOR THIS PROBLEM. C H0=HSTART TOL=10.0E+0**(-NN) XOUT=20.0E+0 10 CONTINUE CALL MEBDF(N, X, H0, Y, XOUT, XEND, TOL, 21, INDEX, LOUT, + LWORK, WORK, LIWORK, IWORK, FCN) IF ((INDEX.NE.0).AND.(INDEX.NE.3)) THEN IF(INDEX.EQ.1) THEN INDEX=0 GO TO 10 ENDIF WRITE( 6 ,15) INDEX 15 FORMAT(' ***INTEGRATION HAS FAILED*** WITH INDEX=',I3) STOP END IF C C HAVE COMPLETED ONE STEP C 60 CONTINUE C C HAVE WE FINISHED YET? C IF( INDEX.EQ.0) THEN C C THEN WE HAVE EFFECTIVELY HIT TOUT C X=XOUT ELSE INDEX=3 GOTO 10 END IF WRITE(6, 690) TOL 690 FORMAT(1X, 'REQUESTED TOLERANCE',5X, 1PE10.3) WRITE(6,700) NFE,NJE 700 FORMAT(1X,I5,2X,'FUNCTION EVALUATIONS',10X, 1 I5,2X,'JACOBIAN EVALUATIONS') ERR1 = Y(1) - 4.539992969929191E-05 ERR2 = Y(2) - 2.061153036149920E-09 ERR3 = Y(3) - 2.823006338263857E-56 ERR4 = Y(4) - 5.235792540515384E-52 C C THESE VERY SMALL CONSTANTS SHOULD BE SET TO ZERO IF THERE C ARE LIKELY TO BE DIFFICULTIES DUE TO UNDERFLOW. C WRITE(6,710) 710 FORMAT(1H ,'ERRORS AT ENDPOINT') WRITE(6,720) ERR1, ERR2, ERR3, ERR4 720 FORMAT(1H ,4HERRS, 5(1PE14.5)) 100 CONTINUE STOP END SUBROUTINE FCN(T,Y,YDOT) DIMENSION Y(4),YDOT(4) YDOT(1)=-0.5E+0*Y(1) YDOT(2) = -Y(2) YDOT(3) = -100.0E+0*Y(3) YDOT(4) = -90.0E+0*Y(4) RETURN END SUBROUTINE PDERV(T,Y,PW) DIMENSION Y(4),PW(4,4) DO 50 I = 1, 4 DO 50 J = 1, 4 PW(I,J) = 0.0E+0 50 CONTINUE PW(1,1)=-0.5E+0 PW(2,2) = -1.0E+0 PW(3,3) = -100.0E+0 PW(4,4) = -90.0E+0 RETURN END SUBROUTINE MEBDF(N,T0,H0,Y0,TOUT,TEND,EPS,MF,INDEX,LOUT,LWORK, + WORK,LIWORK,IWORK,FCN) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C THIS SUBROUTINE IS FOR THE PURPOSE * C OF SPLITTING UP THE WORK ARRAYS WORK AND IWORK * C FOR USE INSIDE THE INTEGRATOR STIFF * C <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< C THIS IS THE JULY 1990 VERSION OF OVDRIV, A PACKAGE FOR C THE SOLUTION OF THE INITIAL VALUE PROBLEM FOR SYSTEMS OF C ORDINARY DIFFERENTIAL EQUATIONS C DY/DT = F(Y,T), Y=(Y(1),Y(2),Y(3), . . . ,Y(N)) C SUBROUTINE OVDRIV IS A DRIVER ROUTINE FOR THIS PACKAGE C C REFERENCES C C 1. J.R. CASH, ON THE INTEGRATION OF STIFF SYSTEMS OF O.D.E.S C USING EXTENDED BACKWARD DIFFERENTIATION FORMULAE, C NUMER. MATH. 34, 235-246 (1980) C 2. THE INTEGRATION OF STIFF INITIAL VALUE PROBLEMS IN O.D.E.S C USING MODIFIED EXTENDED BACKWARD DIFFERENTIATION FORMULAE, C COMP. AND MATHS. WITH APPLICS., 9, 645-657, (1983). C 3. J.R. CASH, STABLE RECURSIONS WITH APPLICATIONS TO THE C NUMERICAL SOLUTION OF STIFF SYSTEMS, ACADEMIC PRESS.(1979) C 4. A.C. HINDMARSH, GEAR.. ORDINARY DIFFERENTIAL EQUATION C SYSTEM SOLVER, UCID-30001 REV. 3, LAWRENCE LIVERMORE C LABORATORY, P.O. BOX 808, LIVERMORE, CA 94550, DEC. 1974. C 5. J.R. CASH AND S. CONSIDINE, AN MEBDF CODE FOR STIFF INITIAL C VALUE PROBLEMS, ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE, C 1992. C C ---------------------------------------------------------------- C OVDRIV IS TO BE CALLED ONCE FOR EACH OUTPUT VALUE OF T, AND C IN TURN MAKES REPEATED CALLS TO THE CORE INTEGRATOR STIFF. C C THE INPUT PARAMETERS ARE .. C N = THE NUMBER OF FIRST ORDER DIFFERENTIAL EQUATIONS. C T0 = THE INITIAL VALUE OF T, THE INDEPENDENT VARIABLE C (USED ONLY ON THE FIRST CALL) C H0 = THE NEXT STEP SIZE IN T (USED FOR INPUT ONLY ON THE C FIRST CALL) C Y0 = A VECTOR OF LENGTH N CONTAINING THE INITIAL VALUES OF Y C (USED FOR INPUT ONLY ON FIRST CALL) C TEND = END OF THE RANGE OF INTEGRATION. C TOUT = THE VALUE OF T AT WHICH OUTPUT IS DESIRED NEXT. C INTEGRATION WILL NORMALLY GO SLIGHTLY BEYOND TOUT C AND THE PACKAGE WILL INTERPOLATE TO T = TOUT C MUST HAVE TEND .GE. TOUT IF H0 IS POSITIVE AND C TOUT .GE. TEND OTHERWISE. C EPS = THE RELATIVE ERROR BOUND (USED ONLY ON THE FIRST C CALL ,UNLESS INDEX = -1). SINGLE STEP ERROR C ESTIMATES DIVIDED BY YMAX(I) WILL BE KEPT LESS THAN C EPS IN ROOT-MEAN-SQUARE NORM (I.E. EUCLIDEAN NORM C DIVIDED BY SQRT(N) ). THE VECTOR YMAX OF WEIGHTS IS C COMPUTED IN OVDRIV. INITIALLY YMAX(I) IS ABS(Y(I)) C WITH A DEFAULT VALUE OF 1 IF Y(I)=0 INITIALLY. C THEREAFTER, YMAX(I) IS THE LARGEST VALUE OF ABS(Y(I)) C SEEN SO FAR, OR THE INITIAL VALUE YMAX(I) IF THAT IS C LARGER C MF = THE METHOD FLAG. AT PRESENT ONLY MF=21 OR 22 IS ALLOWED C WHICH IS EXTENDED BACKWARD DIFFERENTIATION FORMULAE C USING THE CHORD METHOD WITH ANALYTIC OR NUMERICAL C JACOBIAN. THE USER NEEDS TO SPECIFY SUBROUTINE PDERV C IF MF=21 (SEE BELOW) C INDEX = INTEGER USED ON INPUT TO INDICATE TYPE OF CALL, C 1 THIS IS THE FIRST CALL FOR THE PROBLEM C INDEX HAS TO BE SET = 0 AFTER FIRST CALL. C 0 THIS IS NOT THE FIRST CALL FOR THIS PROBLEM C AND INTEGRATION IS TO CONTINUE UP TO T=TOUT C OR BEYOND. C -1 THIS IS NOT THE FIRST CALL FOR THE PROBLEM, C AND THE USER HAS RESET N, EPS, AND/OR MF. C 2 SAME AS 0 EXCEPT THAT TOUT HAS TO BE HIT C EXACTLY (NO INTERPOLATION IS DONE). C ASSUMES TOUT .GE. THE CURRENT T. C 3 ONE STEP MODE. SAME AS 0 EXCEPT CONTROL C RETURNS TO CALLING PROGRAM AFTER ONE STEP. C SINCE THE NORMAL OUTPUT VALUE OF INDEX IS 0, C IT NEED NOT BE RESET FOR NORMAL CONTINUATION. C LOUT = LOGICAL OUTPUT DEVICE. C LWORK = DIMENSION OF REAL WORK ARRAY WORK. LWORK MUST BE AT C LEAST 2*N*N + 37*N C WORK = REAL WORK ARRAY. C LIWORK= DIMENSION OF INTEGER WORK ARRAY. LIWORK MUST BE AT C LEAST N. C IWORK = INTEGER WORK ARRAY. C C AFTER THE INITIAL CALL, IF A NORMAL RETURN OCCURED AND A NORMAL C CONTINUATION IS DESIRED, SIMPLY RESET TOUT AND CALL AGAIN. C ALL OTHER PARAMETERS WILL BE READY FOR THE NEXT CALL. C A CHANGE OF PARAMETERS WITH INDEX = -1 CAN BE MADE AFTER C EITHER A SUCCESSFUL OR AN UNSUCCESSFUL RETURN. C C THE OUTPUT PARAMETERS ARE.. C T0 = THE VALUE OF T WHICH RELATES TO THE CURRENT SOLUTION C POINT Y0() C H0 = THE STEPSIZE H USED LAST, WHETHER SUCCESSFULLY OR NOT. C Y0 = THE COMPUTED VALUES OF Y AT T = TOUT C TOUT = UNCHANGED FROM ITS INPUT VALUE. C INDEX = INTEGER USED ON OUTPUT TO INDICATE RESULTS, WITH C THE FOLLOWING VALUES AND MEANINGS.. C C 3 ONE STEP MODE IS BEING USED AND A SUCCESSFUL STEP HAS BEEN C COMPLETED. THE USER SHOULD CALL MEBDF AGAIN WITH INDEX = 3. C 1 ONE SUCCESSFUL STEP HAS BEEN COMPLETED AND THE ONE STEP MODE C IS NOT BEING USED. THE USER MUST NOW SET INDEX = 0 AND CALL C MEBDF AGAIN. C 0 INTEGRATION WAS COMPLETED TO TOUT OR BEYOND. C C -1 THE INTEGRATION WAS HALTED AFTER FAILING TO PASS THE C ERROR TEST EVEN AFTER REDUCING H BY A FACTOR OF C 1.E10 FROM ITS INITIAL VALUE. C C -2 AFTER SOME INITIAL SUCCESS, THE INTEGRATION WAS C HALTED EITHER BY REPEATED ERROR TEST FAILURES OR BY C A TEST ON EPS. TOO MUCH ACCURACY HAS BEEN REQUESTED. C C -3 THE INTEGRATION WAS HALTED AFTER FAILING TO ACHIEVE C CORRECTOR CONVERGENCE EVEN AFTER REDUCING H BY A C FACTOR OF 1.E10 FROM ITS INITIAL VALUE. C C -4 IMMEDIATE HALT BECAUSE OF ILLEGAL VALUES OF INPUT C PARAMETERS. SEE PRINTED MESSAGE. C C -5 INDEX WAS -1 ON INPUT, BUT THE DESIRED CHANGES OF C PARAMETERS WERE NOT IMPLEMENTED BECAUSE TOUT C WAS NOT BEYOND T. INTERPOLATION AT T = TOUT WAS C PERFORMED AS ON A NORMAL RETURN. TO TRY AGAIN, C SIMPLY CALL AGAIN WITH INDEX = -1 AND A NEW TOUT. C C -10 INSUFFICIENT REAL WORKSPACE FOR THE INTEGRATION C C -11 INSUFFICIENT INTEGER WORKSPACE FOR THE INTEGRATION C C C IN ADDITION TO OVDRIV, THE FOLLOWING ROUTINES ARE PROVIDED C IN THE PACKAGE.. C C INTERP( - ) INTERPOLATES TO GET THE OUTPUT VALUES C AT T=TOUT FROM THE DATA IN THE Y ARRAY. C STIFF( - ) IS THE CORE INTEGRATOR ROUTINE. IT PERFORMS A C SINGLE STEP AND ASSOCIATED ERROR CONTROL. C COSET( - ) SETS COEFFICIENTS FOR BACKWARD DIFFERENTIATION C SCHEMES FOR USE IN THE CORE INTEGRATOR. C PSET( - ) COMPUTES AND PROCESSES THE JACOBIAN C MATRIX J = DF/DY C DEC( - ) PERFORMS AN LU DECOMPOSITION ON A MATRIX. C SOL( - ) SOLVES LINEAR SYSTEMS A*X = B AFTER DEC C HAS BEEN CALLED FOR THE MATRIX A C C C THE FOLLOWING ROUTINES ARE TO BE SUPPLIED BY THE USER.. C C FCN(T,Y,YDOT) COMPUTES THE FUNCTION YDOT = F(Y,T), THE C RIGHT HAND SIDE OF THE O.D.E. C HERE Y AND YDOT ARE VECTORS OF LENGTH N C PDERV(T,Y,PD) COMPUTES THE N*N JACOBIAN MATRIX OF C PARTIAL DERIVATIVES AND STORES IT IN PD C AS AN N BY N ARRAY. PD(I,J) IS TO BE SET C TO THE PARTIAL DERIVATIVE OF YDOT(I) WITH C RESPECT TO Y(J). PDERV IS CALLED ONLY IF C MITER = 1. OTHERWISE A DUMMY ROUTINE CAN C BE SUBSTITUTED. C C THE DIMENSION OF PW BELOW MUST BE AT LEAST N**2. THE DIMENSIONS C OF YMAX,ERROR,SAVE1,SAVE2,IPIV AND THE FIRST DIMENSION C OF Y SHOULD ALL BE AT LEAST N. C C COMMON BLOCKS OF INTEREST TO THE USER. C C /MEBDF2/ C C HUSED LAST STEPSIZE SUCCESSFULLY USED BY THE INTEGRATOR C NQUSED ORDER LAST SUCCESSFULLY USED C NSTEP NUMBER OF STEPS TAKEN SO FAR C NFE NUMBER OF FUNCTION EVALUATIONS USED SO FAR C NJE NUMBER OF JACOBIAN EVALUATIONS USED SO FAR C NDEC NUMBER OF LU DECOMPOSITIONS USED SO FAR C NBSOL NUMBER OF 'BACKSOLVES' USED SO FAR C NPSET NUMBER OF TIMES A NEW COEFFICIENT MATRIX HAS BEEN C FORMED SO FAR C NCOSET NUMBER OF TIMES THE ORDER OF THE METHOD USED HAS BEEN C CHANGED SO FAR C MAXORD THE MAXIMUM ORDER USED SO FAR IN THE INTEGRATION. C TO ACCESS THESE THE USER SHOULD PUT THE FOLLOWING COMMON C STATEMENT IN HIS MAIN PROGRAM: C COMMON/MEBDF2/HUSED,NQUSED,NSTEP,NFE,NJE,NDEC,NBSOL, C + NPSET,NCOSET,MAXORD C C ***************************************************************** C THE USER WILL NEED TO SET THE MACHINE PRECISION CONSTANTS FOR C HIS MACHINE. THESE ARE C UROUND - THE UNIT ROUNDOFF ERROR. C EPSJAC - THE SQUARE ROOT OF UROUND. C THESE ARE SET IN A BLOCK DATA STATEMENT FOLLOWING THIS DRIVER. C ****************************************************************** C C .. SCALAR ARGUMENTS .. INTEGER INDEX,LIWORK,LOUT,LWORK,MF,N C .. C .. ARRAY ARGUMENTS .. DIMENSION WORK(LWORK),Y0(N) INTEGER IWORK(LIWORK) C .. C .. LOCAL SCALARS .. INTEGER I1,I2,I3,I4,I5,I6,I7,I8,I9 C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL OVDRIV,FCN C .. C .. SAVE STATEMENT .. SAVE I1,I2,I3,I4,I5,I6,I7,I8,I9 C .. IF (INDEX.EQ.1) THEN IF (N.LE.0) THEN WRITE (LOUT,9020) N INDEX = -4 ELSE I1 = N*15 + 1 I2 = I1 + N*15 I3 = I2 + N*2 I4 = I3 + N I5 = I4 + N I6 = I5 + N I7 = I6 + N I8 = I7 + N I9 = I8 + N*N IF (LWORK.LT. (I9+N*N-1)) THEN INDEX = -10 WRITE (LOUT,9000) I9 + N*N - 1 END IF IF (LIWORK.LT.N) THEN INDEX = -11 WRITE (LOUT,9010) N END IF END IF IF (INDEX.LT.0) RETURN END IF CALL OVDRIV(N,T0,H0,Y0,TOUT,TEND,EPS,MF,INDEX,LOUT,WORK,WORK(I1), + WORK(I2),WORK(I3),WORK(I4),WORK(I5),WORK(I6),WORK(I7), + WORK(I8),WORK(I9),IWORK,FCN) C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C WORK() HOUSES THE FOLLOWING ARRAYS C C Y(N,15) , YHOLD(N,15) , YNHOLD(N,2) , YMAX(N) C ERRORS(N) , SAVE1(N) , SAVE2(N) , ARH(N) , PW(N*N) C PWCOPY(N*N) C <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< RETURN 9000 FORMAT (/,/,' ***** ERROR ***** INTEGRATION HALTED IN MEBDF',/,/, + ' >>> REAL WORKSPACE IS INSUFFICIENT <<< ',/, + ' WORKSPACE MUST BE AT LEAST ',I8,' ELEMENTS LONG') 9010 FORMAT (/,/,' ***** ERROR ***** INTEGRATION HALTED IN MEBDF',/,/, + ' >>> INTEGER WORKSPACE IS INSUFFICIENT <<< ',/, + ' WORKSPACE MUST BE AT LEAST ',I6,' ELEMENTS LONG') 9020 FORMAT (/,/,' ***** ERROR ***** INTEGRATION HALTED IN MEBDF',/,/, + ' >>> ILLEGAL VALUE FOR NUMBER OF EQUATIONS <<< ',/, + ' WITH N = ',I6) END BLOCK DATA C C MACHINE DEPENDENT CONSTANTS ARE SET HERE C C UROUND THE SMALLEST +VE NUMBER SUCH THAT C C 1.0D+0 + UROUND .NE. 1.0D+0 C C EPSJAC THIS IS USED IN FORMING A NUMERICAL JACOBIAN AND IS C C EPSJAC = DSQRT( UROUND ) C C .. COMMON BLOCKS .. DOUBLE PRECISION UROUND, EPSJAC COMMON /MACHNE/UROUND,EPSJAC C .. C .. DATA STATEMENTS .. DATA UROUND,EPSJAC/2.2D-16,8.660254D-8/ C .. END SUBROUTINE OVDRIV(N,T0,H0,Y0,TOUT,TEND,EPS,MF,INDEX,LOUT,Y,YHOLD, + YNHOLD,YMAX,ERRORS,SAVE1,SAVE2,ARH,PW,PWCOPY, + IPIV,FCN) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C START OF PROGRAM PROPER C C .. SCALAR ARGUMENTS .. INTEGER INDEX,LOUT,MF,N C .. C .. ARRAY ARGUMENTS .. DIMENSION ARH(N),ERRORS(N),PW(*),PWCOPY(*),SAVE1(N),SAVE2(N) + ,Y(N,15),Y0(N),YHOLD(N,15),YMAX(N),YNHOLD(N,2) INTEGER IPIV(N) C .. C .. LOCAL SCALARS .. INTEGER I,KGO,NHCUT C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL INTERP,STIFF,FCN C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DABS,DMAX1 C .. C .. COMMON BLOCKS .. COMMON /FAILS/NFAIL1,NFAIL2,NT COMMON /MACHNE/UROUND,EPSJAC 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 /MEBDF6/CRATE1,CRATE2,TCRAT1,TCRAT2,AVNEW2,AVOLD2,AVNEWJ, + AVOLDJ COMMON /MEBDFL/CFAIL,JNEWIM,SAMPLE INTEGER 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 CFAIL,JNEWIM,SAMPLE C .. IF (INDEX.EQ.0) THEN C I.E. NORMAL CONTINUATION OF INTEGRATION T0=T HMAX = DABS(TEND-T0)*10.0D+0 IF ((T-TOUT)*H.GE.0.0D+0) THEN C 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.EQ.2) THEN C I.E. CONTINUING INTEGRATION BUT WISH TO HIT TOUT T0 = T HMAX = DABS(TEND-T0)*10.0D+0 IF (((T+H)-TOUT)*H.GT.0.0D+0) THEN C WE HAVE ALREADY OVERSHOT THE END OR WE WILL C DO SO ON THE NEXT STEP IF (((T-TOUT)*H.GE.0.0D+0) .OR. (DABS(T-TOUT).LE. + 100.0D+0*UROUND*HMAX)) THEN C HAVE OVERSHOT, SO INTERPOLATE CALL INTERP(N,JSTART,H,T,Y,TOUT,Y0) T0 = TOUT H0 = H INDEX = KFLAG RETURN ELSE C WILL PASS TOUT ON NEXT STEP WITH CURRENT STEPSIZE C 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.EQ.-1) THEN C NOT FIRST CALL BUT PARAMETERS RESET T0 = T IF ((T-TOUT)*H.GE.0.0D+0) THEN C 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.EQ.3) THEN T0 = T IF ((T-TOUT)*H.GE.0.0D+0) THEN C HAVE OVERSHOT,SO INTERPOLATE CALL INTERP(N,JSTART,H,T,Y,TOUT,Y0) INDEX = KFLAG T0 = TOUT H0 = H RETURN END IF ELSE C INDEX SHOULD BE 1 AND THIS IS THE FIRST CALL FOR THIS PROBLEM C CHECK THE ARGUMENTS THAT WERE PASSED FOR CORRECTNESS IF (INDEX.NE.1) THEN C VALUE OF INDEX NOT ALLOWED WRITE (LOUT,9070) INDEX INDEX = -4 END IF IF (EPS.LE.0.0D+0) THEN C ILLEGAL VALUE FOR RELATIVE ERROR TOLERENCE WRITE (LOUT,9040) INDEX = -4 END IF IF (N.LE.0) THEN C ILLEGAL VALUE FOR THE NUMBER OF EQUATIONS WRITE (LOUT,9050) INDEX = -4 END IF IF ((T0-TOUT)*H0.GE.0.0D+0) THEN C PARAMETERS FOR INTEGRATION ARE ILLEGAL WRITE (LOUT,9060) INDEX = -4 END IF IF((TEND-TOUT)*H0.LT.-100.0D+0*UROUND) THEN C TOUT IS BEYOND TEND. WRITE(LOUT,9110) TOUT,TEND INDEX = -4 ENDIF IF ((MF.NE.21) .AND. (MF.NE.22)) THEN C ILLEGAL VALUE FOR METHOD FLAG WRITE (LOUT,9090) MF INDEX = -4 END IF IF (INDEX.NE.1) THEN RETURN ELSE C THE INITIAL PARAMETERS ARE O.K. SO INITIALISE EVERYTHING C ELSE NECESSARY FOR THE INTEGRATION. C IF VALUES OF YMAX OTHER THAN THOSE SET BELOW ARE DESIRED, C THEY SHOULD BE SET HERE. ALL YMAX(I) MUST BE POSITIVE. IF C VALUES FOR HMIN OR HMAX, THE BOUNDS ON DABS(H), OTHER THAN C THOSE BELOW ARE DESIRED, THEY SHOULD BE SET BELOW DO 10 I = 1,N YMAX(I) = DABS(Y0(I)) IF (YMAX(I).EQ.0.0D+0) YMAX(I) = 1.0D+0 Y(I,1) = Y0(I) 10 CONTINUE NCOPY = N T = T0 H = H0 HMIN = DABS(H0) HMAX = DABS(T0-TEND)*10.0D+0 EPSCPY = EPS MFCOPY = MF JSTART = 0 NHCUT = 0 END IF END IF C <<<<<<<<<<<<<<<<< C < TAKE A STEP > C <<<<<<<<<<<<<<<<< 20 IF ((T+H).EQ.T) THEN WRITE (LOUT,9000) NT = NT + 1 END IF CALL STIFF(EPS,H,HMAX,HMIN,JSTART,KFLAG,MF,T,TEND,Y,N,YMAX, + ERRORS,SAVE1,SAVE2,PW,PWCOPY,YHOLD,YNHOLD,ARH,IPIV,LOUT,FCN) KGO = 1 - KFLAG IF (KGO.EQ.1) THEN C NORMAL RETURN FROM STIFF GO TO 30 ELSE IF (KGO.EQ.2) THEN C COULD NOT ACHIEVE ERROR WITH HMIN C SO CHOP HMIN IF WE HAVEN'T DONE SO 10 TIMES GO TO 60 ELSE IF (KGO.EQ.3) THEN C ERROR SMALLER THAN CAN BE HANDLED FOR THIS PROBLEM WRITE (LOUT,9010) T,H GO TO 70 ELSE IF (KGO.EQ.4) THEN C COULD NOT ACHIEVE CONVERGENCE WITH HMIN WRITE (LOUT,9030) T GO TO 60 END IF 30 CONTINUE C --------------------------------------------------------------------- C NORMAL RETURN FROM THE INTEGRATOR. C C THE WEIGHTS YMAX(I) ARE UPDATED. IF DIFFERENT VALUES ARE DESIRED, C THEY SHOULD BE SET HERE. A TEST IS MADE FOR EPS BEING TOO SMALL C FOR THE MACHINE PRECISION. C C ANY OTHER TESTS OR CALCULATIONS THAT ARE REQUIRED AFTER EVERY C STEP SHOULD BE INSERTED HERE. C C IF INDEX = 3, Y0 IS SET TO THE CURRENT Y VALUES ON RETURN. C IF INDEX = 2, H IS CONTROLLED TO HIT TOUT (WITHIN ROUNDOFF C ERROR), AND THEN THE CURRENT Y VALUES ARE PUT IN Y0 ON RETURN. C FOR ANY OTHER VALUE OF INDEX, CONTROL RETURNS TO THE INTEGRATOR C UNLESS TOUT HAS BEEN REACHED. THEN INTERPOLATED VALUES OF Y ARE C COMPUTED AND STORED IN Y0 ON RETURN. C IF INTERPOLATION IS NOT DESIRED, THE CALL TO INTERP SHOULD BE C REMOVED AND CONTROL TRANSFERRED TO STATEMENT 500 INSTEAD OF 520. C -------------------------------------------------------------------- D = 0.0D+0 DO 40 I = 1,N AYI = DABS(Y(I,1)) YMAX(I) = DMAX1(YMAX(I),AYI) D = D + (AYI/YMAX(I))**2 40 CONTINUE D = D* (UROUND/EPS)**2 IF (D.GT.DBLE(FLOAT(N))) THEN C EPS SMALLER THAN MACHINE PRECISION WRITE (LOUT,9020) T KFLAG = -2 GO TO 70 END IF IF (INDEX.EQ.3.OR.INDEX.EQ.1) GO TO 70 IF (DABS(T-TOUT).LE.DABS(15.0D+0*UROUND*TOUT)) THEN C EFFECTIVELY WE HAVE HIT TOUT INDEX = KFLAG T0 = TOUT DO 50 I = 1,N Y0(I) = Y(I,1) 50 CONTINUE H0 = H RETURN END IF IF (INDEX.EQ.2) THEN C CONTINUING INTEGRATION BUT MUST HIT TOUT EXACTLY IF (((T+H)-TOUT)*H.GT.0.0D+0) THEN C WE HAVE ALREADY OVERSHOT THE END OR WE WILL DO C SO ON THE NEXT STEP IF (((T-TOUT)*H.GE.0.0D+0) .OR. (DABS(T-TOUT).LE. + 100.0D+0*UROUND*HMAX)) THEN C HAVE OVERSHOT, SO INTERPOLATE CALL INTERP(N,JSTART,H,T,Y,TOUT,Y0) T0 = TOUT H0 = H INDEX = KFLAG RETURN ELSE C WILL PASS TOUT ON NEXT STEP WITH CURRENT STEPSIZE C 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.GE.0.0D+0) THEN C 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 C ------------------------------------------------------------------- C ON AN ERROR RETURN FROM THE INTEGRATOR, AN IMMEDIATE RETURN OCCURS C IF KFLAG = -2, AND RECOVERY ATTEMPTS ARE MADE OTHERWISE. C H AND HMIN ARE REDUCED BY A FACTOR OF .1 UP TO 10 TIMES C BEFORE GIVING UP. C -------------------------------------------------------------------- 60 CONTINUE IF (NHCUT.EQ.10) THEN C 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.GT.0) THEN JSTART = -1 ELSE C THE INITIAL STEPSIZE WAS TOO LARGE. RESET CERTAIN VARIABLES C 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(DABS(T-TOUT).GT.1000.0D+0*UROUND) THEN DO 80 I = 1,N Y0(I) = Y(I,1) 80 CONTINUE T0 = T ELSE C HAVE EITHER PASSED TOUT OR WE ARE EXTREMELY CLOSE TO IT C SO INTERPOLATE. CALL INTERP(N,JSTART,H,T,Y,TOUT,Y0) T0 = TOUT INDEX = KFLAG END IF H0 = H IF(KFLAG.NE.0) INDEX = KFLAG RETURN C -------------------------- 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 .GE. 0.',/,/) 9070 FORMAT (/,/,' ILLEGAL INPUT.. INDEX =',I5,/,/) 9080 FORMAT (/,/,' INDEX = -1 ON INPUT WITH (T-TOUT)*H .GE. 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 .GE.0' ) END SUBROUTINE INTERP(N,JSTART,H,T,Y,TOUT,Y0) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C .. SCALAR ARGUMENTS .. INTEGER JSTART,N C .. C .. ARRAY ARGUMENTS .. DIMENSION Y(N,15),Y0(N) C .. C .. LOCAL SCALARS .. INTEGER I,J,L C .. C .. INTRINSIC FUNCTIONS .. C .. DO 10 I = 1,N Y0(I) = Y(I,1) 10 CONTINUE L = JSTART + 2 S = (TOUT-T)/H S1 = 1.0D+0 DO 30 J = 2,L S1 = S1* (S+DBLE(FLOAT(J-2)))/DBLE(FLOAT(J-1)) DO 20 I = 1,N Y0(I) = Y0(I) + S1*Y(I,J) 20 CONTINUE 30 CONTINUE RETURN C -------------- END OF SUBROUTINE INTERP --------------------------- END SUBROUTINE COSET(NQ,EL,ELST,TQ,MAXDER) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C -------------------------------------------------------------------- C COSET IS CALLED BY THE INTEGRATOR AND SETS THE COEFFICIENTS USED C BY THE CONVENTIONAL BACKWARD DIFFERENTIATION SCHEME. THE VECTOR C EL OF LENGTH NQ+1 DETERMINES THE BASIC METHOD. THE VECTOR ELST C OF DIMENSION NQ+2 DETERMINES THE MODIFIED EXTENDED BACKWARD C DIFFERENTIATION FORMULAE COEFFICIENTS. THE VECTOR TQ OF C LENGTH 4 IS INVOLVED IN ADJUSTING THE STEPSIZE IN RELATION TO THE C TRUNCATION ERROR. ITS VALUES ARE GIVEN BY THE PERTST ARRAY. THE C VECTORS EL AND TQ BOTH DEPEND ON METH AND NQ. COSET ALSO SETS C MAXDER, THE MAXIMUM ORDER OF THE METHOD AVAILABLE. THE C COEFFICIENTS IN PERTST NEED TO BE GIVEN TO ONLY ABOUT ONE PERCENT C ACCURACY. THE ORDER IN WHICH THE GROUPS APPEAR BELOW IS C COEFFICIENTS FOR ORDER NQ-1, COEFFICIENTS FOR ORDER NQ, C COEFFICIENTS FOR ORDER NQ+1. C ------------------------------------------------------------------- C .. SCALAR ARGUMENTS .. INTEGER MAXDER,NQ C .. C .. ARRAY ARGUMENTS .. DIMENSION EL(10),ELST(10),TQ(4) C .. C .. LOCAL SCALARS .. INTEGER K C .. C .. LOCAL ARRAYS .. DIMENSION PERTST(8,3) C .. C .. INTRINSIC FUNCTIONS .. C .. C .. COMMON BLOCKS .. COMMON /MEBDF2/HUSED,NQUSED,NSTEP,NFE,NJE,NDEC,NBSOL,NPSET,NCOSET, + MAXORD INTEGER MAXORD,NBSOL,NCOSET,NDEC,NFE,NJE,NPSET,NQUSED,NSTEP C .. C .. DATA STATEMENTS .. DATA PERTST(1,1)/1./,PERTST(2,1)/2./,PERTST(3,1)/4.5/, + PERTST(4,1)/7.333/,PERTST(5,1)/10.42/,PERTST(6,1)/13.7/, + PERTST(7,1)/17.15/,PERTST(8,1)/20.74/ DATA PERTST(1,2)/2./,PERTST(2,2)/4.5/,PERTST(3,2)/7.333/, + PERTST(4,2)/10.42/,PERTST(5,2)/13.7/,PERTST(6,2)/17.15/, + PERTST(7,2)/20.74/,PERTST(8,2)/24.46/ DATA PERTST(1,3)/4.5/,PERTST(2,3)/7.333/,PERTST(3,3)/10.42/, + PERTST(4,3)/13.7/,PERTST(5,3)/17.15/,PERTST(6,3)/20.74/, + PERTST(7,3)/24.46/,PERTST(8,3)/1./ C .. C ------------------------------------------------------------------- C THE FOLLOWING COEFFICIENTS SHOULD BE DEFINED TO MACHINE ACCURACY. C THEIR DERIVATION IS GIVEN IN REFERENCE 2. C ------------------------------------------------------------------- IF (NQ.GT.MAXORD) MAXORD = NQ MAXDER = 6 NCOSET = NCOSET + 1 GO TO (10,20,30,40,50,60,70) NQ 10 EL(1) = 1.0D+0 ELST(1) = 1.5D+0 ELST(3) = -0.5D+0 GO TO 80 20 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 GO TO 80 30 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 GO TO 80 40 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 GO TO 80 50 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 GO TO 80 60 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 GO TO 80 70 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 80 DO 90 K = 1,3 TQ(K) = PERTST(NQ,K) 90 CONTINUE TQ(4) = 0.5D+0*TQ(2)/DBLE(FLOAT(NQ)) RETURN C --------------------- END OF SUBROUTINE COSET --------------------- END SUBROUTINE PSET(Y,N,H,T,UROUND,EPSJAC,CON,MITER,IER,NRENEW,YMAX, + SAVE2,PW,PWCOPY,WRKSPC,IPIV,FCN) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C ------------------------------------------------------------------- C PSET IS CALLED BY STIFF TO COMPUTE AND PROCESS THE COEFFICIENT C MATRIX I - H*EL(1)*J WHERE J IS AN APPROXIMATION TO C THE RELEVANT JACOBIAN. THIS MATRIX IS THEN SUBJECTED TO LU C DECOMPOSITION IN PREPARATION FOR LATER SOLUTION OF LINEAR SYSTEMS C OF ALGEBRAIC EQUATIONS WITH LU AS THE COEFFICIENT MATRIX. THE C MATRIX J IS FOUND BY THE USER-SUPPLIED ROUTINE PDERV IF MITER=1 C OR BY FINITE DIFFERENCING IF MITER = 2. C IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION WITH C PSET USES THE FOLLOWING .. C EPSJAC = DSQRT(UROUND), USED IN NUMERICAL JACOBIAN INCREMENTS. C ******************************************************************* C THE ARGUMENT NRENEW IS USED TO SIGNAL WHETHER OR NOT C WE REQUIRE A NEW JACOBIAN TO BE CALCULATED. C C IF NRENEW > 0 THEN WE REQUIRE A NEW J TO BE COMPUTED C = 0 THEN USE A COPY OF THE LAST J COMPUTED C ******************************************************************* C .. SCALAR ARGUMENTS .. INTEGER IER,MITER,N,NRENEW C .. C .. ARRAY ARGUMENTS .. DIMENSION PW(*),PWCOPY(*),SAVE2(N),WRKSPC(N),Y(N,15),YMAX(N) INTEGER IPIV(N) C .. C .. LOCAL SCALARS .. INTEGER I,J,J1,JJKK C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL DEC,FCN,PDERV C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DABS,DMAX1,DSQRT C .. C .. COMMON BLOCKS .. COMMON /MEBDF2/HUSED,NQUSED,NSTEP,NFE,NJE,NDEC,NBSOL,NPSET,NCOSET, + MAXORD INTEGER MAXORD,NBSOL,NCOSET,NDEC,NFCN1,NFE,NJE,NPSET, + NQUSED,NSTEP C .. NPSET = NPSET + 1 IF (NRENEW.EQ.0) THEN DO 10 I = 1,N*N PW(I) = PWCOPY(I)*CON 10 CONTINUE GO TO 70 END IF IF (MITER.EQ.2) GO TO 30 CALL PDERV(T,Y,PW) DO 20 I = 1,N*N PWCOPY(I) = PW(I) PW(I) = PW(I)*CON 20 CONTINUE NJE = NJE + 1 GO TO 70 30 NJE = NJE + 1 D = 0.0D+0 DO 40 I = 1,N D = D + SAVE2(I)**2 40 CONTINUE R0 = DABS(H)*DSQRT(D)*1.0D+03*UROUND J1 = 0 DO 60 J = 1,N YJ = Y(J,1) R = EPSJAC*YMAX(J) R = DMAX1(R,R0) Y(J,1) = Y(J,1) + R D = CON/R CALL FCN(T,Y,WRKSPC) DO 50 I = 1,N JJKK = I + J1 TEMPRY = (WRKSPC(I)-SAVE2(I)) PWCOPY(JJKK) = TEMPRY/R PW(JJKK) = TEMPRY*D 50 CONTINUE Y(J,1) = YJ J1 = J1 + N 60 CONTINUE NFE = NFE + N 70 J = 1 DO 80 I = 1,N PW(J) = PW(J) + 1.0D+0 J = J + (N+1) 80 CONTINUE CALL DEC(N,N,PW,IPIV,IER) * NDEC = NDEC + 1 RETURN C ---------------------- END OF SUBROUTINE PSET --------------------- END SUBROUTINE DEC(N,NDIM,A,IP,IER) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C ------------------------------------------------------------------- C MATRIX TRIANGULARISATION BY GAUSSIAN ELIMINATION C INPUT.. C N = ORDER OF MATRIX. C NDIM = DECLARED DIMENSION OF ARRAY A. C A = MATRIX TO BE TRIANGULARISED. C OUTPUT.. C A(I,J), I.LE.J = UPPER TRIANGULAR FACTOR, U. C A(I,J), I.GT.J = MULTIPLIERS = LOWER TRIANGULAR FACTOR, I-L. C IP(K), K.LT.N = INDEX OF KTH PIVOT ROW. C IP(N) = (-1)**(NUMBER OF INTERCHANGES) OR 0. C IER = 0 IF MATRIX IS NON-SINGULAR, OR K IF FOUND TO BE SINGULAR C AT STAGE K. C USE SOL TO OBTAIN SOLUTION OF LINEAR SYSTEM. C DETERM(A) = IP(N)*A(1,1)*A(2,2)* . . . *A(N,N). C IF IP(N) = 0, A IS SINGULAR, SOL WILL DIVIDE BY ZERO. C C REFERENCE. C C.B. MOLER, ALGORITHM 423, LINEAR EQUATION SOLVER, C.A.C.M C 15 (1972), P.274. C ------------------------------------------------------------------ C .. SCALAR ARGUMENTS .. INTEGER IER,N,NDIM C .. C .. ARRAY ARGUMENTS .. DIMENSION A(NDIM,N) INTEGER IP(N) C .. C .. LOCAL SCALARS .. INTEGER I,J,K,KP1,M,NM1 C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DABS C .. C .. COMMON BLOCKS .. IER = 0 IP(N) = 1 IF (N.EQ.1) GO TO 70 NM1 = N - 1 DO 60 K = 1,NM1 KP1 = K + 1 M = K DO 10 I = KP1,N IF (DABS(A(I,K)).GT.DABS(A(M,K))) M = I 10 CONTINUE IP(K) = M T = A(M,K) IF (M.EQ.K) GO TO 20 IP(N) = -IP(N) A(M,K) = A(K,K) A(K,K) = T 20 IF (T.EQ.0.0D+0) GO TO 80 T = 1.0D+0/T DO 30 I = KP1,N A(I,K) = -A(I,K)*T 30 CONTINUE DO 50 J = KP1,N T = A(M,J) A(M,J) = A(K,J) A(K,J) = T IF (T.EQ.0.0D+0) GO TO 50 DO 40 I = KP1,N A(I,J) = A(I,J) + A(I,K)*T 40 CONTINUE 50 CONTINUE 60 CONTINUE 70 K = N IF (A(N,N).EQ.0.0D+0) GO TO 80 RETURN 80 IER = K IP(N) = 0 RETURN C --------------------- END OF SUBROUTINE DEC ---------------------- END SUBROUTINE SOL(N,NDIM,A,B,IP) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C .. SCALAR ARGUMENTS .. INTEGER N,NDIM C .. C .. ARRAY ARGUMENTS .. DIMENSION A(NDIM,N),B(N) INTEGER IP(N) C .. C .. LOCAL SCALARS .. INTEGER I,K,KB,KM1,KP1,M,NM1 C .. C .. COMMON BLOCKS .. C ------------------------------------------------------------------ C SOLUTION OF LINEAR SYSTEM, A*X = B. C INPUT .. C N = ORDER OF MATRIX. C NDIM = DECLARED DIMENSION OF MATRIX A. C A = TRIANGULARISED MATRIX OBTAINED FROM DEC. C B = RIGHT HAND SIDE VECTOR. C IP = PIVOT VECTOR OBTAINED FROM DEC. C DO NOT USE IF DEC HAS SET IER .NE. 0 C OUTPUT.. C B = SOLUTION VECTOR, X. C ------------------------------------------------------------------ IF (N.EQ.1) GO TO 50 NM1 = N - 1 DO 20 K = 1,NM1 KP1 = K + 1 M = IP(K) T = B(M) B(M) = B(K) B(K) = T DO 10 I = KP1,N B(I) = B(I) + A(I,K)*T 10 CONTINUE 20 CONTINUE DO 40 KB = 1,NM1 KM1 = N - KB K = KM1 + 1 B(K) = B(K)/A(K,K) T = -B(K) DO 30 I = 1,KM1 B(I) = B(I) + A(I,K)*T 30 CONTINUE 40 CONTINUE 50 B(1) = B(1)/A(1,1) RETURN C ------------------------- END OF SUBROUTINE SOL ------------------ END SUBROUTINE ERRORS(N,TQ,EPS,EDN,E,EUP,BND) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C *************************************************** C C THIS ROUTINE CALCULATES ERRORS USED IN TESTS C IN STIFF . C C *************************************************** C .. SCALAR ARGUMENTS .. INTEGER N C .. C .. ARRAY ARGUMENTS .. DIMENSION TQ(4) C .. C .. LOCAL SCALARS .. C .. C .. INTRINSIC FUNCTIONS .. C .. SQHOL = DBLE(FLOAT(N))*EPS*EPS EDN = TQ(1)*TQ(1)*SQHOL C C ** ERROR ASSOCIATED WITH LOWER ORDER METHOD C E = TQ(2)*TQ(2)*SQHOL C C ** ERROR ASSOCIATED WITH PRESENT ORDER C EUP = TQ(3)*TQ(3)*SQHOL C C ** ERROR ASSOCIATED WITH HIGHER ORDER METHOD C BND = TQ(4)*TQ(4)*SQHOL*0.3D+0 RETURN END SUBROUTINE PRDICT(T,H,Y,L,N,YPRIME,FCN) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C ********************************************************************** C PREDICTS A VALUE FOR Y AT (T+H) GIVEN THE HISTORY ARRAY AT Y(T) C THEN EVALUATES THE DERIVATIVE AT THIS POINT, THE RESULT OF THIS C EVALUATION BEING STORED IN YPRIME() C ********************************************************************** C .. SCALAR ARGUMENTS .. INTEGER L,N C .. C .. ARRAY ARGUMENTS .. DIMENSION Y(N,15),YPRIME(N) C .. C .. LOCAL SCALARS .. INTEGER I,J2 C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL FCN C .. C .. COMMON BLOCKS .. COMMON /MEBDF2/HUSED,NQUSED,NSTEP,NFE,NJE,NDEC,NBSOL,NPSET,NCOSET, + MAXORD INTEGER MAXORD,NBSOL,NCOSET,NDEC,NFE,NJE,NPSET,NQUSED,NSTEP C .. DO 20 J2 = 2,L DO 10 I = 1,N Y(I,1) = Y(I,1) + Y(I,J2) 10 CONTINUE 20 CONTINUE T = T + H CALL FCN(T,Y,YPRIME) NFE = NFE + 1 RETURN END SUBROUTINE ITRAT2(Y,N,T,HBETA,ERRBND,ARH,CRATE,TCRATE,M,WORKED, + YMAX,ERROR,SAVE1,SAVE2,PW,IPIV,LMB,FCN) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C .. SCALAR ARGUMENTS .. INTEGER M,N LOGICAL WORKED C .. C .. ARRAY ARGUMENTS .. DIMENSION ARH(N),ERROR(N),PW(*),SAVE1(N),SAVE2(N),Y(N,15),YMAX(N) INTEGER IPIV(N) C .. C .. LOCAL SCALARS .. INTEGER I C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL FCN,SOL C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DMAX1,DMIN1 C .. C .. COMMON BLOCKS .. COMMON /MEBDF2/HUSED,NQUSED,NSTEP,NFE,NJE,NDEC,NBSOL,NPSET,NCOSET, + MAXORD INTEGER MAXORD,NBSOL,NCOSET,NDEC,NFE,NJE,NPSET,NQUSED,NSTEP C .. C .. DATA STATEMENTS .. DATA ZERO/0.0D+0/ C .. IF(LMB.EQ.1) GOTO 25 DO 10 I = 1,N SAVE1(I) = -SAVE1(I) + HBETA*SAVE2(I) + ARH(I) 10 CONTINUE CALL SOL(N,N,PW,SAVE1,IPIV) NBSOL = NBSOL + 1 D = ZERO DO 20 I = 1,N ERROR(I) = ERROR(I) + SAVE1(I) D = D + (SAVE1(I)/YMAX(I))**2 SAVE1(I) = Y(I,1) + ERROR(I) 20 CONTINUE TCRATE = TCRATE + CRATE D1 = D M = 1 CALL FCN(T,SAVE1,SAVE2) NFE = NFE + 1 25 CONTINUE WORKED = .TRUE. 30 DO 40 I = 1,N SAVE1(I) = -SAVE1(I) + HBETA*SAVE2(I) + ARH(I) 40 CONTINUE C C IF WE ARE HERE THEN PARTIALS ARE O.K. C CALL SOL(N,N,PW,SAVE1,IPIV) NBSOL = NBSOL + 1 C C WE NOW CALCULATE A WEIGHTED RMS TYPE NORM C D = ZERO DO 50 I = 1,N ERROR(I) = ERROR(I) + SAVE1(I) D = D + (SAVE1(I)/YMAX(I))**2 SAVE1(I) = Y(I,1) + ERROR(I) 50 CONTINUE C ------------------------------------------------------------------- C TEST FOR CONVERGENCE. IF M.GT.0 , AN ESTIMATE OF THE CONVERGENCE C RATE CONSTANT IS STORED IN CRATE, AND THIS IS USED IN THE TEST. C ------------------------------------------------------------------- IF (M.NE.0) THEN IF (D1.NE.ZERO) CRATE = DMAX1(0.9D+0*CRATE,D/D1) END IF TCRATE = TCRATE + CRATE IF ((D*DMIN1(1.0D+0,2.0D+0*CRATE)).LT.ERRBND) RETURN IF (M.NE.0) THEN IF (D.GT.D1) THEN WORKED = .FALSE. RETURN END IF END IF D1 = D IF (M.EQ.3) THEN WORKED = .FALSE. RETURN END IF M = M + 1 CALL FCN(T,SAVE1,SAVE2) NFE = NFE + 1 GO TO 30 END SUBROUTINE STIFF(EPS,H,HMAX,HMIN,JSTART,KFLAG,MF,T,TEND,Y, + N,YMAX,ERROR,SAVE1,SAVE2,PW,PWCOPY,YHOLD,YNHOLD,ARH, + IPIV,LOUT,FCN) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C ------------------------------------------------------------------ C STIFF PERFORMS ONE STEP OF THE INTEGRATION OF AN INITIAL VALUE C PROBLEM FOR A SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS. C COMMUNICATION WITH STIFF IS DONE WITH THE FOLLOWING VARIABLES.. C Y AN N BY LMAX+4 ARRAY CONTAINING THE DEPENDENT VARIABLES C AND THEIR BACKWARD DIFFERENCES. LMAX-1 =MAXDER IS THE C MAXIMUM ORDER AVAILABLE. SEE SUBROUTINE COSET. C Y(I,J+1) CONTAINS THE JTH BACKWARD DIFFERENCE OF Y(I) C T THE INDEPENDENT VARIABLE. T IS UPDATED ON EACH STEP TAKEN. C H THE STEPSIZE TO BE ATTEMPTED ON THE NEXT STEP. C H IS ALTERED BY THE ERROR CONTROL ALGORITHM DURING C THE PROBLEM. H CAN BE EITHER POSITIVE OR NEGATIVE BUT C ITS SIGN MUST REMAIN CONSTANT THROUGHOUT THE PROBLEM. C HMIN THE MINIMUM AND MAXIMUM ABSOLUTE VALUE OF THE STEPSIZE C HMAX TO BE USED FOR THE STEP. THESE MAY BE CHANGED AT ANY C TIME BUT WILL NOT TAKE EFFECT UNTIL THE NEXT H CHANGE. C EPS THE RELATIVE ERROR BOUND. SEE DESCRIPTION IN MEBDF C UROUND THE UNIT ROUNDOFF OF THE MACHINE. C N THE NUMBER OF FIRST ORDER DIFFERENTIAL EQUATIONS. C MF THE METHOD FLAG. MUST BE SET TO 21 OR 22 AT PRESENT C KFLAG A COMPLETION CODE WITH THE FOLLOWING MEANINGS.. C 0 THE STEP WAS SUCCESSFUL C -1 THE REQUESTED ERROR COULD NOT BE ACHIEVED C WITH ABS(H) = HMIN. C -2 THE REQUESTED ERROR IS SMALLER THAN CAN C BE HANDLED FOR THIS PROBLEM. C -3 CORRECTOR CONVERGENCE COULD NOT BE C ACHIEVED FOR ABS(H)=HMIN. C ON A RETURN WITH KFLAG NEGATIVE, THE VALUES OF T AND C THE Y ARRAY ARE AS AT THE BEGINNING OF THE LAST C STEP ATTEMPTED, AND H IS THE LAST STEP SIZE ATTEMPTED. C JSTART AN INTEGER USED ON INPUT AND OUTPUT. C ON INPUT IT HAS THE FOLLOWING VALUES AND MEANINGS.. C 0 PERFORM THE FIRST STEP. C .GT.0 TAKE A NEW STEP CONTINUING FROM THE LAST C .LT.0 TAKE THE NEXT STEP WITH A NEW VALUE OF C H, EPS OR N C ON EXIT, JSTART IS NQUSED, THE ORDER OF THE METHOD LAST USED. C YMAX AN ARRAY OF N ELEMENTS WITH WHICH THE ESTIMATED LOCAL C ERRORS IN Y ARE COMPARED C ERROR AN ARRAY OF N ELEMENTS. C SAVE1,2 TWO ARRAYS OF WORKING SPACE BOTH OF LENGTH N. C PW A BLOCK OF LOCATIONS USED FOR PARTIAL DERIVATIVES C IPIV AN INTEGER ARRAY OF LENGTH N USED FOR PIVOT INFORMATION. C C JNEWIM IS TO INDICATE IF PRESENT ITERATION MATRIX C WAS FORMED USING A NEW J OR OLD J. C JSNOLD KEEPS TRACK OF NUMBER OF STEPS TAKEN WITH C PRESENT ITERATION MATRIX (BE IT FORMED BY C A NEW J OR NOT). C AVNEWJ STORES VALUE FOR AVERAGE CRATE WHEN ITERATION C MATRIX WAS FORMED BY A NEW J. C AVOLDJ STORES VALUE FOR AVERAGE CRATE WHEN ITERATION C MATRIX WAS FORMED BY AN OLD J. C NRENEW FLAG THAT IS USED IN COMMUNICATION WITH SUBROUTINE PSET. C IF NRENEW > 0 THEN FORM A NEW JACOBIAN BEFORE C COMPUTING THE COEFFICIENT MATRIX FOR C THE NEWTON ITERATION C = 0 FORM THE COEFFICIENT MATRIX USING A C COPY OF AN OLD JACOBIAN C NEWPAR FLAG USED IN THIS SUBROUTINE TO INDICATE IF A JACOBIAN C HAS BEEN EVALUATED FOR THE CURRENT STEP C ********************************************************************** C .. SCALAR ARGUMENTS .. INTEGER JSTART,KFLAG,LOUT,MF,N C .. C .. ARRAY ARGUMENTS .. DIMENSION HSTPSZ(2,14) DIMENSION ARH(N),ERROR(N),PW(*),PWCOPY(*),SAVE1(N),SAVE2(N), + Y(N,15), YHOLD(N,15),YMAX(N),YNHOLD(N,2) INTEGER IPIV(N) C .. C .. LOCAL SCALARS .. INTEGER I,IER,IITER,IITER2,ITST,J,J1,J2,KFAIL,LL,LMAX, + M3STEP,MAXDER,METH,MFOLD,MITER,NEWQ LOGICAL FINISH,OVRIDE,WORKED C .. C .. LOCAL ARRAYS .. DIMENSION EL(10),ELST(10),TQ(4) C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL COSET,CPYARY,ERRORS,FCN,HCHOSE,ITRAT2, + PRDICT,PSET,RSCALE,SOL C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DABS,DMAX1,DMIN1 C .. C .. COMMON BLOCKS .. COMMON/STPSZE/HSTPSZ,TOP COMMON /FAILS/NFAIL1,NFAIL2,NT COMMON /MACHNE/UROUND,EPSJAC 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 INTEGER IDOUB,ISAMP,IWEVAL,JCHANG,JSINUP,JSNOLD,L,M1,M2,MAXORD, + MEQC1,MEQC2,MQ1TMP,MQ2TMP,NBSOL,NCOSET,NDEC,NEWPAR,NFAIL1, + NFAIL2,NFE,NJE,NPSET,NQ,NQUSED,NRENEW,NSTEP,NT LOGICAL CFAIL,JNEWIM,SAMPLE C .. C .. SAVE STATEMENT .. SAVE EPSOLD,MFOLD,HOLD,LMAX,EDN,EUP,BND,EL,TQ,ELST,E C .. C .. DATA STATEMENTS .. DATA EL(2),ELST(2),OLDLO/3*1.0D+0/ DATA ZERO,ONE/0.0D+0,1.0D+0/ C .. TOLD = T KFLAG = 0 IF (JSTART.GT.0) GO TO 60 IF (JSTART.NE.0) GO TO 30 C ------------------------------------------------------------------ C ON THE FIRST CALL, THE ORDER IS SET TO 1 AND THE INITIAL YDOT C IS CALCULATED. RMAX IS THE MAXIMUM RATIO BY WHICH H CAN BE C INCREASED IN A SINGLE STEP. RMAX IS SET EQUAL TO 1.D4 INITIALLY C TO COMPENSATE FOR THE SMALL INITIAL H, BUT THEN IS NORMALLY = 10. C IF A FAILURE OCCURS (IN CORRECTOR CONVERGENCE OR ERROR TEST), C RMAX IS SET AT 2. FOR THE NEXT INCREASE. C ------------------------------------------------------------------ CALL FCN(T,Y,SAVE1) DO 10 I = 1,N Y(I,2) = H*SAVE1(I) 10 CONTINUE 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 15 I=1,12 HSTPSZ(1,I)=1.0D+0 HSTPSZ(2,I)=ATOL 15 CONTINUE 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 C ************************************************** C CFAIL=.TRUE. ENSURES THAT WE CALCULATE A NEW C J ON THE FIRST CALL. C ************************************************** MEQC1 = 0 MEQC2 = 0 MQ1TMP = 0 MQ2TMP = 0 NBSOL = 0 C ----------------------------------------------------------------- C IF THE CALLER HAS CHANGED N OR EPS, THE CONSTANTS E, EDN, EUP C AND BND MUST BE RESET. E IS A COMPARISON FOR ERRORS AT THE C CURRENT ORDER NQ. EUP IS TO TEST FOR INCREASING THE ORDER, C EDN FOR DECREASING THE ORDER. BND IS USED TO TEST FOR CONVERGENCE C OF THE CORRECTOR ITERATES. IF THE CALLER HAS CHANGED H, Y MUST C BE RE-SCALED. IF H IS CHANGED, IDOUB IS SET TO L+1 TO PREVENT C FURTHER CHANGES IN H FOR THAT MANY STEPS. C ----------------------------------------------------------------- CALL COSET(NQ,EL,ELST,TQ,MAXDER) LMAX = MAXDER + 1 RC = RC*EL(1)/OLDLO OLDLO = EL(1) IWEVAL = MITER NRENEW = 1 NEWPAR = 0 C ***************************************************** C NRENEW AND NEWPAR ARE TO INSTRUCT ROUTINE THAT C WE WISH A NEW J TO BE CALCULATED FOR THIS STEP. C ***************************************************** CALL ERRORS(N,TQ,EPS,EDN,E,EUP,BND) EPSOLD = EPS DO 20 I = 1,N ARH(I) = EL(2)*Y(I,1) 20 CONTINUE CALL CPYARY(N*L,Y,YHOLD) QI = H*EL(1) QQ = ONE/QI CALL PRDICT(T,H,Y,L,N,SAVE2,FCN) GO TO 110 C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C DIFFERENT PARAMETERS ON THIS CALL < C <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 30 CALL CPYARY(N*L,YHOLD,Y) IF (MF.NE.MFOLD) THEN METH = MF/10 MITER = MF - 10*METH MFOLD = MF IWEVAL = MITER END IF IF (EPS.NE.EPSOLD) THEN CALL ERRORS(N,TQ,EPS,EDN,E,EUP,BND) EPSOLD = EPS END IF IF (H.NE.HOLD) THEN RH = H/HOLD H = HOLD GO TO 50 ELSE GO TO 60 END IF C ********************************************* C RE-SCALE Y AFTER A CHANGE OF STEPSIZE * C ********************************************* 40 RH = DMAX1(RH,HMIN/DABS(H)) 50 RH = DMIN1(RH,HMAX/DABS(H),RMAX) CALL RSCALE(N,L,RH,Y) RMAX = 10.0D+0 JCHANG = 1 H = H*RH RC = RC*RH IF (JSNOLD.GT.40) THEN CFAIL = .TRUE. NEWPAR = 0 RC = ZERO C ********************************************************************** C CFAIL=TRUE AND NEWPAR=0 SHOULD FORCE A NEW J TO BE EVALUATED C AFTER 42 STEPS WITH AN OLD J, IF WE HAVE HAD A FAILURE OF ANY C KIND ON THE FIRST, SECOND OR THIRD STAGE OF THE CURRENT STEP C ********************************************************************** END IF IDOUB = L + 1 CALL CPYARY(N*L,Y,YHOLD) 60 IF (DABS(RC-ONE).GT.0.3D+0) IWEVAL = MITER C ------------------------------------------------------------------ C THIS SECTION COMPUTES THE PREDICTED VALUES OF Y C AND THE RHS, ARH, FOR USE IN THE NEWTON ITERATION SCHEME. C RC IS THE RATIO OF THE NEW TO OLD VALUES OF THE COEFFICIENT C H*EL(1). WHEN RC DIFFERS FROM 1 BY MORE THAN 30 PERCENT, IWEVAL IS C SET TO MITER TO FORCE THE PARTIALS TO BE UPDATED. C ------------------------------------------------------------------ QI = H*EL(1) QQ = ONE/QI DO 70 I = 1,N ARH(I) = EL(2)*Y(I,1) 70 CONTINUE DO 90 J1 = 2,NQ DO 80 I = 1,N ARH(I) = ARH(I) + EL(J1+1)*Y(I,J1) 80 CONTINUE 90 CONTINUE IF (JCHANG.EQ.1) THEN C IF WE HAVE CHANGED STEPIZE THEN PREDICT A VALUE FOR Y(T+H) C AND EVALUATE THE DERIVATIVE THERE (STORED IN SAVE2()) CALL PRDICT(T,H,Y,L,N,SAVE2,FCN) ELSE C ELSE USE THE VALUES COMPUTED FOR THE SECOND BDF FROM THE LAST C STEP. Y( ,LMAX+3) HOLDS THE VALUE FOR THE DERIVATIVE AT (T+H) C AND Y( ,LMAX+4) HOLDS THE APPROXIMATION TO Y AT THIS POINT. DO 100 I = 1,N SAVE2(I) = Y(I,LMAX+3) Y(I,1) = Y(I,LMAX+4) 100 CONTINUE T = T + H END IF 110 IF (IWEVAL.LE.0) GO TO 120 C ------------------------------------------------------------------- C IF INDICATED, THE MATRIX P = I - H*EL(1)*J IS RE-EVALUATED BEFORE C STARTING THE CORRECTOR ITERATION. IWEVAL IS SET = 0 TO INDICATE C THAT THIS HAS BEEN DONE. P IS COMPUTED AND PROCESSED IN PSET. C THE PROCESSED MATRIX IS STORED IN PW C ------------------------------------------------------------------- IWEVAL = 0 RC = ONE IITER = MEQC1 - MQ1TMP IITER2 = MEQC2 - MQ2TMP IF (JNEWIM) THEN IF (JSNOLD.GE.3) THEN AVNEWJ = TCRAT1/DBLE(FLOAT(IITER)) AVNEW2 = TCRAT2/DBLE(FLOAT(IITER2)) ELSE AVNEWJ = ONE AVNEW2 = ONE END IF ELSE C C MATRIX P WAS FORMED WITH A COPY OF J C IF (JSNOLD.GE.3) THEN AVOLDJ = TCRAT1/DBLE(FLOAT(IITER)) AVOLD2 = TCRAT2/DBLE(FLOAT(IITER2)) IF (AVOLDJ.LT.AVNEWJ) THEN AVNEWJ = AVOLDJ ELSE IF (((DABS(AVOLDJ-AVNEWJ)).GT.0.3D+0) .OR. + ((AVOLDJ.GT.0.85D+0).AND. (AVOLDJ.NE.ONE))) THEN C C SINCE IN CERTAIN INSTANCES AVOLDJ WILL C BE 1.0 AND THERE WILL BE NO NEED TO C UPDATE J. C CFAIL = .TRUE. CRATE1 = ONE CRATE2 = ONE END IF ELSE CFAIL = .TRUE. CRATE1 = ONE CRATE2 = ONE C C ********************************************************* C IF WE HAVE REACHED HERE THINGS MUST HAVE GONE WRONG C ********************************************************* C END IF END IF TCRAT1 = ZERO TCRAT2 = ZERO IF (CFAIL) THEN IF (NEWPAR.EQ.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) C NOTE THAT ERROR() IS JUST BEING USED AS A WORKSPACE BY PSET IF (IER.NE.0) THEN C 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 130 I = 1,N SAVE1(I) = Y(I,1) ERROR(I) = ZERO 130 CONTINUE M1 = 0 C ********************************************************************** C UP TO 4 CORRECTOR ITERATIONS ARE TAKEN. A CONVERGENCE TEST IS C MADE ON THE R.M.S. NORM OF EACH CORRECTION ,USING BND, WHICH C DEPENDS ON EPS. THE SUM OF THE CORRECTIONS IS ACCUMULATED IN THE C VECTOR ERROR(I). THE Y ARRAY IS NOT ALTERED IN THE CORRECTOR C LOOP. THE UPDATED Y VECTOR IS STORED TEMPORARILY IN SAVE1. C ********************************************************************** 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 C C NOW TEST TO SEE IF IT WAS SUCCESSFUL OR NOT C C IF (.NOT.WORKED) THEN NFAIL1 = NFAIL1 + 1 C ********************************************************************** C THE CORRECTOR ITERATION FAILED TO CONVERGE IN 4 TRIES. IF C PARTIALS ARE NOT UP TO DATE, THEY ARE RE-EVALUATED FOR THE C NEXT TRY. OTHERWISE THE Y ARRAY IS RETRACTED TO ITS VALUES C BEFORE PREDICTION AND H IS REDUCED IF POSSIBLE. IF NOT A C NON-CONVERGENCE EXIT IS TAKEN C ********************************************************************** IF (IWEVAL.EQ.-1) THEN C HAVE BEEN USING OLD PARTIALS, UPDATE THEM AND TRY AGAIN IWEVAL = MITER CFAIL = .TRUE. CALL FCN(T,Y,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 140 I = 1,N SAVE2(I) = (SAVE1(I)-ARH(I))*QQ Y(I,1) = SAVE1(I) 140 CONTINUE C C UPDATE THE DIFFERENCES AT N+1 C DO 160 J = 2,L DO 150 I = 1,N Y(I,J) = Y(I,J-1) - YHOLD(I,J-1) 150 CONTINUE 160 CONTINUE C C COMPUTE ERROR IN THE SOLUTION C D = ZERO DO 170 I = 1,N D = D + ((Y(I,L)-YHOLD(I,L))/YMAX(I))**2 170 CONTINUE C C STORING Y FROM FIRST STEP FOR USE IN THIRD STEP C DO 180 I = 1,N YNHOLD(I,1) = Y(I,1) YNHOLD(I,2) = SAVE2(I) 180 CONTINUE IF (D.GT.E) GO TO 330 IF ((M1+M2).GE.ITST) THEN M2 = 0 EBDF1 = 0.5D+0/DBLE(FLOAT(L)) PRBDF1 = ((D/ (E*0.7D+0))**EBDF1)*1.5D+0 + 1.6D-6 RHBDF1 = ONE/PRBDF1 CALL HCHOSE(RHBDF1,H,OVRIDE) RHLMIT = 1 - (5+ITST-M1-M2)*0.5D-1 IF (RHBDF1.LT.RHLMIT) THEN IJUS=1 RED=RHBDF1 GOTO 450 END IF END IF KFAIL = 0 C ---------------------------------------------------------------------- DO 190 I = 1,N ARH(I) = EL(2)*Y(I,1) 190 CONTINUE DO 210 J1 = 2,NQ DO 200 I = 1,N ARH(I) = ARH(I) + EL(J1+1)*Y(I,J1) 200 CONTINUE 210 CONTINUE CALL PRDICT(T,H,Y,L,N,SAVE2,FCN) DO 220 I = 1,N SAVE1(I) = Y(I,1) ERROR(I) = ZERO 220 CONTINUE M2 = 0 C C FOR NOW WILL ASSUME THAT WE DO NOT WISH TO SAMPLE C AT THE N+2 STEP POINT C CALL ITRAT2(Y,N,T,QI,BND,ARH,CRATE2,TCRAT2,M2,WORKED,YMAX,ERROR, + SAVE1,SAVE2,PW,IPIV,1,FCN) MEQC2 = MEQC2 + M2 + 1 C C NOW CHECK TO SEE IF IT WAS SUCCESSFUL OR NOT C IF (.NOT.WORKED) THEN NFAIL2 = NFAIL2 + 1 IJUS=0 RED=0.5D+0 GOTO 450 END IF C C IF WE ARE DOWN TO HERE THEN THINGS MUST HAVE CONVERGED C DO 230 I = 1,N Y(I,LMAX+3) = (SAVE1(I)-ARH(I))*QQ Y(I,LMAX+4) = SAVE1(I) 230 CONTINUE IF ((M2.GE.2) .AND. (CRATE2.LT.0.35D+0)) THEN DO 250 J = 2,NQ DO 240 I = 1,N SAVE1(I) = SAVE1(I) - Y(I,J) 240 CONTINUE 250 CONTINUE DSTEP2 = ZERO DO 260 I = 1,N DSTEP2 = DSTEP2 + ((SAVE1(I)-YNHOLD(I,1)-Y(I,L))/YMAX(I))**2 260 CONTINUE EBDF2 = 0.5D+0/DBLE(FLOAT(L)) PRBDF2 = ((DSTEP2/ (E*0.7D+0))**EBDF2)*1.5D+0 + 1.6D-6 RHBDF2 = ONE/PRBDF2 CALL HCHOSE(RHBDF2,H,OVRIDE) IF (RHBDF2.LT.0.9D+0) THEN IJUS=1 RED=RHBDF2 GOTO 450 END IF END IF C C WE ARE NOW COMPUTING THE THIRD STAGE C LL = L + 1 T = TOLD + H DO 280 I = 1,N ARH(I) = H* (ELST(NQ+2)*Y(I,LMAX+3)+ + (ELST(1)-EL(1))*YNHOLD(I,2)) DO 270 J1 = 1,NQ ARH(I) = ARH(I) + ELST(J1+1)*YHOLD(I,J1) 270 CONTINUE 280 CONTINUE DO 290 I = 1,N SAVE2(I) = YNHOLD(I,2) Y(I,1) = YNHOLD(I,1) 290 CONTINUE M3STEP = 0 300 DO 310 I = 1,N SAVE1(I) = -Y(I,1) + QI*SAVE2(I) + ARH(I) 310 CONTINUE CALL SOL(N,N,PW,SAVE1,IPIV) NBSOL = NBSOL + 1 D = ZERO DO 320 I = 1,N D = D + (SAVE1(I)/YMAX(I))**2 Y(I,1) = Y(I,1) + SAVE1(I) 320 CONTINUE IF ((D*DMIN1(ONE,2.0D+0*CRATE1)).LE.BND) GO TO 360 IF (M3STEP.EQ.4) THEN WRITE (LOUT,9000) IJUS=1 RED=0.5D+0 GO TO 450 END IF M3STEP = M3STEP + 1 CALL FCN(T,Y,SAVE2) NFE = NFE + 1 GO TO 300 330 KFAIL = KFAIL - 1 C ********************************************************************** C THE ERROR TEST FAILED. KFAIL KEEPS TRACK OF MULTIPLE FAILURES. C RESTORE T AND THE Y ARRAY TO THEIR PREVIOUS VALUES AND PREPARE TO C TRY THE STEP AGAIN. COMPUTE THE OPTIMAL STEP SIZE FOR THIS ORDER C AND ONE ORDER LOWER. C ********************************************************************** T = TOLD HOLD = H EFAIL = 0.5D+0/DBLE(FLOAT(L)) CALL CPYARY(N*L,YHOLD,Y) RMAX = 2.0D+0 IF (DABS(H).LE.HMIN*1.00001D+0) THEN C C REQUESTED ERROR NOT POSSIBLE WITH GIVEN HMIN C KFLAG = -1 HOLD = H RETURN END IF IF (KFAIL.LE.-3) GO TO 340 C C PREDICTING A NEW H AFTER INSUFFICIENT ACCURACY C PRFAIL = ((D/E)**EFAIL)*1.5D+0 + 1.6D-6 NEWQ = NQ RH = ONE/ (PRFAIL*DBLE(FLOAT(-KFAIL))) CALL HCHOSE(RH,H,OVRIDE) GO TO 40 C ********************************************************************** C CONTROL REACHES THIS STAGE IF 3 OR MORE FAILURES HAVE OCCURED. C IT IS ASSUMED THAT THE DERIVATIVES THAT HAVE ACCUMULATED IN THE Y C ARRAY HAVE ERRORS OF THE WRONG ORDER. HENCE THE FIRST DERIVATIVE C IS RE-COMPUTED, AND THE ORDER IS SET TO 1. THEN H IS REDUCED BY A C FACTOR OF 10, AND THE STEP IS RETRIED. AFTER A TOTAL OF 7 C FAILURES AN EXIT IS TAKEN WITH KFLAG=-2. C ********************************************************************** 340 IF (KFAIL.EQ.-7) THEN C ERROR SMALLER THAN CAN BE HANDLED FOR PROBLEM KFLAG = -2 HOLD = H RETURN END IF C ********************************* C START FROM ORDER 1 AGAIN * C ********************************* JCHANG = 1 RH = DMAX1(HMIN/DABS(H),0.1D+0) CALL HCHOSE(RH,H,OVRIDE) H = H*RH CALL FCN(T,YHOLD,SAVE1) NFE = NFE + 1 DO 350 I = 1,N Y(I,1) = YHOLD(I,1) Y(I,2) = H*SAVE1(I) YHOLD(I,2) = Y(I,2) 350 CONTINUE IWEVAL = MITER CFAIL = .TRUE. C SINCE WE HAVE HAD PROBLEMS PROCEED WITH THIS ORDER C FOR 10 STEPS (IF WE CAN) IDOUB = 10 IF (NQ.EQ.1) GO TO 60 NQ = 1 L = 2 C 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) C NOW JUMP TO NORMAL CONTINUATION POINT GO TO 60 C ********************************************************************** C THE ITERATION FOR THE CORRECTED SOLUTION HAS CONVERGED. C UPDATE THE Y ARRAY. C ********************************************************************** 360 CONTINUE DO 380 J2 = 2,LL DO 370 I = 1,N Y(I,J2) = Y(I,J2-1) - YHOLD(I,J2-1) 370 CONTINUE 380 CONTINUE C ---------------------------------------------------------------------- C IF THE COUNTER IDOUB EQUALS 2 , STORE Y(I,LMAX+5) WHICH IS USED IN C ASSESSING THE POSSIBILITY OF INCREASING THE ORDER. IF IDOUB = 0 C CONTROL PASSES TO 480 WHERE AN ATTEMPT TO CHANGE THE STEPSIZE AND C ORDER IS MADE. C ---------------------------------------------------------------------- IF (IDOUB.EQ.2) THEN DO 390 I = 1,N Y(I,LMAX+5) = Y(I,LL) 390 CONTINUE END IF IDOUB = IDOUB - 1 TRANGE=(TEND-TOLD-H)*H IF(TRANGE.LT.0.0D+0) GOTO 440 JCHANG = 0 IF (IDOUB.EQ.0) THEN SAMPLE = .FALSE. ISAMP = ISAMP + 1 IF (ISAMP.EQ.4) THEN SAMPLE = .TRUE. ISAMP = 0 END IF C ********************************************************************** C NOW COMPUTE THE FACTORS PR1, PR2 AND PR3, BY WHICH C H COULD BE DIVIDED AT ORDER NQ-1, ORDER NQ AND ORDER NQ+1 C RESPECTIVELY. THE SMALLEST OF THESE IS DETERMINED AND THE NEW C ORDER CHOSEN ACCORDINGLY. IF THE ORDER IS TO BE INCREASED WE C MUST COMPUTE ONE MORE BACKWARD DIFFERENCE. C ********************************************************************** PR3 = 1.D+20 IF (L.NE.LMAX) THEN DUP = ZERO DO 400 I = 1,N DUP = DUP + ((Y(I,LL)-Y(I,LMAX+5))/YMAX(I))**2 400 CONTINUE ENQ3 = 0.5D+0/DBLE(FLOAT(L+1)) PR3 = ((DUP/EUP)**ENQ3)*1.7D+0 + 1.8D-6 END IF ENQ2 = 0.5D+0/DBLE(FLOAT(L)) D = ZERO DO 410 I = 1,N D = D + ((Y(I,L)-YHOLD(I,L))/YMAX(I))**2 410 CONTINUE PR2 = ((D/E)**ENQ2)*1.5D+0 + 1.6D-6 PR1 = 1.D+20 IF (NQ.GT.1) THEN DDOWN = ZERO DO 420 I = 1,N DDOWN = DDOWN + (Y(I,L)/YMAX(I))**2 420 CONTINUE ENQ1 = 0.5D+0/DBLE(FLOAT(NQ)) PR1 = ((DDOWN/EDN)**ENQ1)*1.6D+0 + 1.7D-6 END IF IF (PR2.LE.PR3) THEN IF (PR2.GT.PR1) THEN NEWQ = NQ - 1 RH = 1.0D+0/PR1 ELSE NEWQ = NQ RH = 1.0D+0/PR2 END IF ELSE IF (PR3.LT.PR1) THEN NEWQ = L RH = 1.0D+0/PR3 ELSE NEWQ = NQ - 1 RH = 1.0D+0/PR1 END IF RH = DMIN1(RH,RMAX) CALL HCHOSE(RH,H,OVRIDE) IF ((JSINUP.LE.20).AND.(KFLAG.EQ.0).AND.(RH.LT.1.1D+0)) THEN C WE HAVE RUN INTO PROBLEMS IDOUB = 10 NQ = NQUSED L = NQ + 1 GO TO 440 END IF C ********************************************************************** C IF THERE IS A CHANGE IN ORDER, RESET NQ, L AND THE C COEFFICIENTS. IN ANY CASE H IS RESET ACCORDINGLY TO RH AND THE C Y ARRAY IS RE-SCALED C ********************************************************************** IF (NEWQ.NE.NQ) THEN IF (NEWQ.GT.NQ) THEN C ADD AN EXTRA TERM TO THE HISTORY ARRAY DO 430 I = 1,N Y(I,LL) = Y(I,L) - YHOLD(I,L) 430 CONTINUE END IF NQ = NEWQ L = NQ + 1 C 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 C NOW RESCALE THE HISTORY ARRAY FOR THE NEW STEPSIZE RH = DMAX1(RH,HMIN/DABS(H)) RH = DMIN1(RH,HMAX/DABS(H),RMAX) CALL RSCALE(N,L,RH,Y) RMAX = 10.0D+0 JCHANG = 1 H = H*RH RC = RC*RH IF (JSNOLD.GT.40) THEN RC = ZERO C IF WE HAVE BEEN USING THE SAME COEFFICIENT MATRIX FOR MORE C THAN 40 STEPS FORCE A NEW ONE TO BE FORMED ON THE NEXT STEP. C CODE WILL ATTEMPT TO USE A COPY OF THE LAST JACOBIAN FORMED END IF IDOUB = L + 1 END IF 440 CONTINUE C ---------------------------------------------------------------------- C STORE THE Y ARRAY IN THE MATRIX YHOLD. STORE IN THE Y ARRAY THE C INFORMATION NECASSARY TO PERFORM AN INTERPOLATION TO FIND THE C SOLUTION AT THE SPECIFIED OUTPUT POINT IF APPROPRIATE. C ---------------------------------------------------------------------- CALL CPYARY(N*L,Y,YHOLD) NSTEP = NSTEP + 1 JSINUP = JSINUP + 1 JSNOLD = JSNOLD + 1 JSTART = NQUSED T = TOLD + HUSED HOLD = H KFAIL = 0 NEWPAR = 0 CFAIL = .FALSE. RETURN 450 CONTINUE FINISH = .FALSE. T=TOLD RMAX=2.0D+0 DO 460 J1=1,L DO 460 I=1,N Y(I,J1)=YHOLD(I,J1) 460 CONTINUE IF(DABS(H).LE.HMIN*1.00001D+0) THEN C C CORRECTOR CONVERGENCE COULD NOT BE ACHIEVED C IF(NSTEP.EQ.0) THEN KFLAG=-1 ELSE KFLAG=-3 END IF C C TO SUPPRESS ERROR MESSAGES AT START AS H MAY C HAVE BEEN TOO LARGE ON THE FIRST STEP. C HOLD=H FINISH = .TRUE. END IF RH = RED C C TRY AGAIN WITH UPDATED PARTIALS C IF(IJUS.EQ.0) CALL HCHOSE(RH,H,OVRIDE) IF(.NOT.FINISH) THEN GO TO 40 ELSE RETURN END IF C ------------------- END OF SUBROUTINE STIFF -------------------------- 9000 FORMAT (1X,' CORRECTOR HAS NOT CONVERGED') END SUBROUTINE RSCALE(N,L,RH,Y) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C .. SCALAR ARGUMENTS .. INTEGER L,N C .. C .. ARRAY ARGUMENTS .. DIMENSION Y(N,15) C .. C .. LOCAL SCALARS .. INTEGER I,J,J1 C .. C .. LOCAL ARRAYS .. DIMENSION DI(8,8) C .. C .. DATA STATEMENTS .. C SUBROUTINE IS FOR RESCALING THE HISTORY ARRAY AFTER A CHANGE IN C STEPSIZE C C N ORDER OF THE PROBLEM C L NUMBER OF TERMS IN THE HISTORY ARRAY TO BE RESCALED C RH RATIO OF THE STEPSIZE CHANGE (I.E. RH = HNEW/HOLD) C Y() THE HISTORY ARRAY C DATA ZERO/0.0D+0/ C .. DI(2,2) = RH IF (L.GT.2) THEN TA = RH*RH DI(2,3) = RH* (1.0D+0-RH)/2.0D+0 DI(3,3) = TA IF (L.GT.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.GT.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.GT.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.GT.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.GT.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 30 I = 1,N DO 20 J = 2,L ZZ = ZERO DO 10 J1 = J,L ZZ = ZZ + DI(J,J1)*Y(I,J1) 10 CONTINUE Y(I,J) = ZZ 20 CONTINUE 30 CONTINUE RETURN END SUBROUTINE CPYARY(NELEM,SOURCE,TARGET) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C C COPIES THE ARRAY SOURCE() INTO THE ARRAY TARGET() C C THIS SUBROUTINE COULD BE REPLACED BY THE BLAS ROUTINE SCOPY C (AFTER CHANGING THE ARGUMENT LIST APPROPRIATELY) C C .. SCALAR ARGUMENTS .. INTEGER NELEM C .. C .. ARRAY ARGUMENTS .. DIMENSION SOURCE(NELEM),TARGET(NELEM) C .. C .. LOCAL SCALARS .. INTEGER I C .. DO 10 I = 1,NELEM TARGET(I) = SOURCE(I) 10 CONTINUE RETURN END SUBROUTINE HCHOSE(RH,H,OVRIDE) IMPLICIT DOUBLE PRECISION(A-H,O-Z) COMMON / STPSZE / HSTPSZ(2,14),TOP LOGICAL OVRIDE C C FIRST MOVE ALL ELEMENTS DOWN ONE PLACE C IF (H.NE.HSTPSZ(2,1)) THEN DO 5 I=12,2,-1 I2=I-1 HSTPSZ(1,I)=HSTPSZ(1,I2) 5 HSTPSZ(2,I)=HSTPSZ(2,I2) C C NOW INSERT VALUE OF H USED BEFORE THIS CALL C HSTPSZ(1,2)=H/HSTPSZ(2,1) HSTPSZ(2,1)=H END IF C C NOW DECIDE ON THE NEW CHANGE C IF (RH.GT.1.0D+0) THEN OVRIDE=.FALSE. ELSE IF (HSTPSZ(1,2).LE.1.0D+0) THEN OVRIDE=.FALSE. ELSE IF ((RH*H).LE.HSTPSZ(2,2)) THEN OVRIDE=.FALSE. ELSE RH=HSTPSZ(2,2)/H OVRIDE=.TRUE. END IF HSTPSZ(1,1)=RH RETURN C C ************************************************************ C END PROGRAM RUNMEB(INPUT,OUTPUT,RESULT,TAPE5=INPUT,TAPE6=RESULT) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION Y(20) C C THE REAL WORK ARRAY NEEDS TO BE AT LEAST 2*N*N + 37*N WORDS LONG C THE INTEGER WORK ARRAY NEEDS TO BE AT LEAST N WORDS LONG C DIMENSION WORK(4680) INTEGER IWORK(40) EXTERNAL FCN COMMON/MEBDF2/HUSED,NQUSED,NSTEP,NFE,NJE,NDEC,NBSOL,NPSET, + NCOSET,MAXORD C C THE INCLUSION OF THE COMMON BLOCK /MEBDF2/ ALLOWS THE USER TO C ACCESS VARIOUS COUNTERS THAT MAY BE OF INTEREST TO HIM. C THESE VARIABLES ARE EXPLAINED IN THE COMMENTS IN SUBROUTINE C MEBDF. C THE USER NEEDS TO SET THE UNIT ROUND-OFF ERROR FOR THE MACHINE C BEING USED IN THE BLOCK DATA STATEMENT FOLLOWING SUBROUTINE C MEBDF. C DATA LOUT / 6 / LWORK = 4680 LIWORK = 40 N=4 DO 100 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.EQ. 3 .OR. NN .EQ. 4) HSTART = 1.0D-4 IF (NN.EQ. 5 .OR. NN .EQ. 6) HSTART = 1.0D-5 IF (NN.EQ. 7 .OR. NN .EQ. 8) HSTART = 1.0D-6 IF (NN.EQ. 9 .OR. NN .EQ. 10) HSTART = 1.0D-7 C C THESE ARE THE INITIAL STEPS CHOSEN BY THE DETEST PACKAGE C FOR THIS PROBLEM C H0=HSTART TOL=10.0D+0**(-NN) XOUT=20.0D+0 10 CONTINUE CALL MEBDF(N, X, H0, Y, XOUT, XEND, TOL, 21, INDEX, LOUT, + LWORK, WORK, LIWORK, IWORK, FCN) IF ((INDEX.NE.0).AND.(INDEX.NE.3)) THEN IF(INDEX.EQ.1) THEN INDEX=0 GO TO 10 ENDIF WRITE( 6 ,15) INDEX 15 FORMAT(' ***INTEGRATION HAS FAILED*** WITH INDEX=',I3) STOP END IF C C HAVE COMPLETED ONE STEP C 60 CONTINUE C C HAVE WE FINISHED YET? C IF( INDEX.EQ.0) THEN C C THEN WE HAVE EFFECTIVELY HIT TOUT C X=XOUT ELSE INDEX=3 GOTO 10 END IF WRITE(6, 690) TOL 690 FORMAT(1X, 'REQUESTED TOLERANCE',5X, 1PE10.3) WRITE(6,700) NFE,NJE 700 FORMAT(1X,I5,2X,'FUNCTION EVALUATIONS',10X, 1 I5,2X,'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 C C THESE VERY SMALL CONSTANTS SHOULD BE SET TO ZERO IF THERE C ARE LIKELY TO BE DIFFICULTIES DUE TO UNDERFLOW. C WRITE(6,710) 710 FORMAT(1H ,'ERRORS AT ENDPOINT') WRITE(6,720) ERR1, ERR2, ERR3, ERR4 720 FORMAT(1H ,4HERRS, 5(1PE14.5)) 100 CONTINUE STOP END SUBROUTINE FCN(T,Y,YDOT) IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION Y(4),YDOT(4) YDOT(1)=-0.5D+0*Y(1) YDOT(2) = -Y(2) YDOT(3) = -100.0D+0*Y(3) YDOT(4) = -90.0D+0*Y(4) RETURN END SUBROUTINE PDERV(T,Y,PW) IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION Y(4),PW(4,4) DO 50 I = 1, 4 DO 50 J = 1, 4 PW(I,J) = 0.0D+0 50 CONTINUE PW(1,1)=-0.5D+0 PW(2,2) = -1.0D+0 PW(3,3) = -100.0D+0 PW(4,4) = -90.0D+0 RETURN END