SUBROUTINE QAGI(F,BOUND,INF,EPSABS,EPSREL,RESULT,ABSERR,NEVAL, 1 IER,LIMIT,LENW,LAST,IWORK,WORK) C***BEGIN PROLOGUE QAGI C***DATE WRITTEN 800101 (YYMMDD) C***REVISION DATE 830518 (YYMMDD) C***CATEGORY NO. H2A3A1,H2A4A1 C***KEYWORDS AUTOMATIC INTEGRATOR,EXTRAPOLATION,GENERAL-PURPOSE, C GLOBALLY ADAPTIVE,INFINITE INTERVALS,TRANSFORMATION C***AUTHOR PIESSENS, ROBERT, APPLIED MATH. AND PROGR. DIV. - C K. U. LEUVEN C DE DONCKER, ELISE, APPLIED MATH. AND PROGR. DIV. - C K. U. LEUVEN C***PURPOSE The routine calculates an approximation result to a given C INTEGRAL I = Integral of F over (BOUND,+INFINITY) C OR I = Integral of F over (-INFINITY,BOUND) C OR I = Integral of F over (-INFINITY,+INFINITY) C Hopefully satisfying following claim for accuracy C ABS(I-RESULT).LE.MAX(EPSABS,EPSREL*ABS(I)). C***DESCRIPTION C C From the book "Numerical Methods and Software" C by D. Kahaner, C. Moler, S. Nash C Prentice Hall 1988 C C Integration over infinite intervals C Standard fortran subroutine C C PARAMETERS C ON ENTRY C F - Real C Function subprogram defining the integrand C function F(X). The actual name for F needs to be C declared E X T E R N A L in the driver program. C C BOUND - Real C Finite bound of integration range C (has no meaning if interval is doubly-infinite) C C INF - Integer C indicating the kind of integration range involved C INF = 1 corresponds to (BOUND,+INFINITY), C INF = -1 to (-INFINITY,BOUND), C INF = 2 to (-INFINITY,+INFINITY). C C EPSABS - Real C Absolute accuracy requested C EPSREL - Real C Relative accuracy requested C If EPSABS.LE.0 C and EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28), C the routine will end with IER = 6. C C C ON RETURN C RESULT - Real C Approximation to the integral C C ABSERR - Real C Estimate of the modulus of the absolute error, C which should equal or exceed ABS(I-RESULT) C C NEVAL - Integer C Number of integrand evaluations C C IER - Integer C IER = 0 normal and reliable termination of the C routine. It is assumed that the requested C accuracy has been achieved. C - IER.GT.0 abnormal termination of the routine. The C estimates for result and error are less C reliable. It is assumed that the requested C accuracy has not been achieved. C ERROR MESSAGES C IER = 1 Maximum number of subdivisions allowed C has been achieved. One can allow more C subdivisions by increasing the value of C LIMIT (and taking the according dimension C adjustments into account). However, if C this yields no improvement it is advised C to analyze the integrand in order to C determine the integration difficulties. If C the position of a local difficulty can be C determined (e.g. SINGULARITY, C DISCONTINUITY within the interval) one C will probably gain from splitting up the C interval at this point and calling the C integrator on the subranges. If possible, C an appropriate special-purpose integrator C should be used, which is designed for C handling the type of difficulty involved. C = 2 The occurrence of roundoff error is C detected, which prevents the requested C tolerance from being achieved. C The error may be under-estimated. C = 3 Extremely bad integrand behaviour occurs C at some points of the integration C interval. C = 4 The algorithm does not converge. C Roundoff error is detected in the C extrapolation table. C It is assumed that the requested tolerance C cannot be achieved, and that the returned C RESULT is the best which can be obtained. C = 5 The integral is probably divergent, or C slowly convergent. It must be noted that C divergence can occur with any other value C of IER. C = 6 The input is invalid, because C (EPSABS.LE.0 and C EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28)) C or LIMIT.LT.1 or LENIW.LT.LIMIT*4. C RESULT, ABSERR, NEVAL, LAST are set to C zero. Exept when LIMIT or LENIW is C invalid, IWORK(1), WORK(LIMIT*2+1) and C WORK(LIMIT*3+1) are set to ZERO, WORK(1) C is set to A and WORK(LIMIT+1) to B. C C DIMENSIONING PARAMETERS C LIMIT - Integer C Dimensioning parameter for IWORK C LIMIT determines the maximum number of subintervals C in the partition of the given integration interval C (A,B), LIMIT.GE.1. In many cases LIMIT = 100 is ok. C If LIMIT.LT.1, the routine will end with IER = 6. C C LENW - Integer C Dimensioning parameter for WORK C LENW must be at least LIMIT*4. C If LENW.LT.LIMIT*4, the routine will end C with IER = 6. C C LAST - Integer C On return, LAST equals the number of subintervals C produced in the subdivision process, which C determines the number of significant elements C actually in the WORK ARRAYS. C C WORK ARRAYS C IWORK - Integer C Vector of dimension at least LIMIT, the first C K elements of which contain pointers C to the error estimates over the subintervals, C such that WORK(LIMIT*3+IWORK(1)),... , C WORK(LIMIT*3+IWORK(K)) form a decreasing C sequence, with K = LAST if LAST.LE.(LIMIT/2+2), and C K = LIMIT+1-LAST otherwise C C WORK - Real C Vector of dimension at least LENW C on return C WORK(1), ..., WORK(LAST) contain the left C end points of the subintervals in the C partition of (A,B), C WORK(LIMIT+1), ..., WORK(LIMIT+LAST) Contain C the right end points, C WORK(LIMIT*2+1), ...,WORK(LIMIT*2+LAST) contain the C integral approximations over the subintervals, C WORK(LIMIT*3+1), ..., WORK(LIMIT*3) C contain the error estimates. C***REFERENCES (NONE) C***ROUTINES CALLED QAGIE,XERROR C***END PROLOGUE QAGI C REAL ABSERR,BOUND,EPSABS,EPSREL,F,RESULT,WORK INTEGER IER,INF,IWORK,LAST,LENW,LIMIT,LVL,L1,L2,L3,NEVAL C DIMENSION IWORK(LIMIT),WORK(LENW) C EXTERNAL F C C CHECK VALIDITY OF LIMIT AND LENW. C C***FIRST EXECUTABLE STATEMENT QAGI IER = 6 NEVAL = 0 LAST = 0 RESULT = 0.0E+00 ABSERR = 0.0E+00 IF(LIMIT.LT.1.OR.LENW.LT.LIMIT*4) GO TO 10 C C PREPARE CALL FOR QAGIE. C L1 = LIMIT+1 L2 = LIMIT+L1 L3 = LIMIT+L2 C CALL QAGIE(F,BOUND,INF,EPSABS,EPSREL,LIMIT,RESULT,ABSERR, 1 NEVAL,IER,WORK(1),WORK(L1),WORK(L2),WORK(L3),IWORK,LAST) C C CALL ERROR HANDLER IF NECESSARY. C LVL = 0 10 IF(IER.EQ.6) LVL = 1 IF(IER.NE.0) CALL XERROR( 'ABNORMAL RETURN FROM QAGI', 1 26,IER,LVL) RETURN END SUBROUTINE EA(NEWFLG,SVALUE,LIMEXP,RESULT,ABSERR,EPSTAB,IERR) C***BEGIN PROLOGUE EA C***DATE WRITTEN 800101 (YYMMDD) C***REVISION DATE 871208 (YYMMDD) C***CATEGORY NO. E5 C***KEYWORDS CONVERGENCE ACCELERATION,EPSILON ALGORITHM,EXTRAPOLATION C***AUTHOR PIESSENS, ROBERT, APPLIED MATH. AND PROGR. DIV. - C K. U. LEUVEN C DE DONCKER-KAPENGA, ELISE,WESTERN MICHIGAN UNIVERSITY C KAHANER, DAVID K., NATIONAL BUREAU OF STANDARDS C STARKENBURG, C. B., NATIONAL BUREAU OF STANDARDS C***PURPOSE Given a slowly convergent sequence, this routine attempts C to extrapolate nonlinearly to a better estimate of the C sequence's limiting value, thus improving the rate of C convergence. Routine is based on the epsilon algorithm C of P. Wynn. An estimate of the absolute error is also C given. C THIS PROLOGUE HAS BEEN REMOVED FOR REASONS OF SPACE C FOR A COMPLETE COPY OF THIS ROUTINE CONTACT THE AUTHORS C From the book "Numerical Methods and Software" C by D. Kahaner, C. Moler, S. Nash C Prentice Hall 1988 C***END PROLOGUE EA REAL ABSERR,DELTA1,DELTA2,DELTA3,EPRN,EPSTAB(*), 1 ERROR,ERR1,ERR2,ERR3,E0,E1,E2,E3,RELPR,RES,RESULT, 2 RES3LA(3),R1MACH,SS,SVALUE,TOL1,TOL2,TOL3 INTEGER I,IB,IB2,IE,IERR,IN,K1,K2,K3,LIMEXP,N,NEWELM,NUM,NRES LOGICAL NEWFLG C C***FIRST EXECUTABLE STATEMENT EA IF(LIMEXP.LT.3) THEN IERR = 1 CALL XERROR('LIMEXP IS LESS THAN 3',21,1,1) GO TO 110 ENDIF IERR = 0 RES3LA(1)=EPSTAB(LIMEXP+5) RES3LA(2)=EPSTAB(LIMEXP+6) RES3LA(3)=EPSTAB(LIMEXP+7) RESULT=SVALUE IF(NEWFLG) THEN N=1 NRES=0 NEWFLG=.FALSE. EPSTAB(N)=SVALUE ABSERR=ABS(RESULT) GO TO 100 ELSE N=INT(EPSTAB(LIMEXP+3)) NRES=INT(EPSTAB(LIMEXP+4)) IF(N.EQ.2) THEN EPSTAB(N)=SVALUE ABSERR=.6E+01*ABS(RESULT-EPSTAB(1)) GO TO 100 ENDIF ENDIF EPSTAB(N)=SVALUE RELPR=R1MACH(4) EPRN=1.0E+01*RELPR EPSTAB(N+2)=EPSTAB(N) NEWELM=(N-1)/2 NUM=N K1=N DO 40 I=1,NEWELM K2=K1-1 K3=K1-2 RES=EPSTAB(K1+2) E0=EPSTAB(K3) E1=EPSTAB(K2) E2=RES DELTA2=E2-E1 ERR2=ABS(DELTA2) TOL2=MAX(ABS(E2),ABS(E1))*RELPR DELTA3=E1-E0 ERR3=ABS(DELTA3) TOL3=MAX(ABS(E1),ABS(E0))*RELPR IF(ERR2.GT.TOL2.OR.ERR3.GT.TOL3) GO TO 10 C C IF E0, E1 AND E2 ARE EQUAL TO WITHIN MACHINE C ACCURACY, CONVERGENCE IS ASSUMED. C RESULT=E2 C ABSERR=ABS(E1-E0)+ABS(E2-E1) C RESULT=RES ABSERR=ERR2+ERR3 GO TO 50 10 IF(I.NE.1) THEN E3=EPSTAB(K1) EPSTAB(K1)=E1 DELTA1=E1-E3 ERR1=ABS(DELTA1) TOL1=MAX(ABS(E1),ABS(E3))*RELPR C C IF TWO ELEMENTS ARE VERY CLOSE TO EACH OTHER, OMIT C A PART OF THE TABLE BY ADJUSTING THE VALUE OF N C IF(ERR1.LE.TOL1.OR.ERR2.LE.TOL2.OR.ERR3.LE.TOL3) GO TO 20 SS=0.1E+01/DELTA1+0.1E+01/DELTA2-0.1E+01/DELTA3 ELSE EPSTAB(K1)=E1 IF(ERR2.LE.TOL2.OR.ERR3.LE.TOL3) GO TO 20 SS=0.1E+01/DELTA2-0.1E+01/DELTA3 ENDIF C C TEST TO DETECT IRREGULAR BEHAVIOUR IN THE TABLE, AND C EVENTUALLY OMIT A PART OF THE TABLE ADJUSTING THE VALUE C OF N C IF(ABS(SS*E1).GT.0.1E-03) GO TO 30 20 N=I+I-1 IF(NRES.EQ.0) THEN ABSERR=ERR2+ERR3 RESULT=RES ELSE IF(NRES.EQ.1) THEN RESULT=RES3LA(1) ELSE IF(NRES.EQ.2) THEN RESULT=RES3LA(2) ELSE RESULT=RES3LA(3) ENDIF GO TO 50 C C COMPUTE A NEW ELEMENT AND EVENTUALLY ADJUST C THE VALUE OF RESULT C 30 RES=E1+0.1E+01/SS EPSTAB(K1)=RES K1=K1-2 IF(NRES.EQ.0) THEN ABSERR=ERR2+ABS(RES-E2)+ERR3 RESULT=RES GO TO 40 ELSE IF(NRES.EQ.1) THEN ERROR=.6E+01*(ABS(RES-RES3LA(1))) ELSE IF(NRES.EQ.2) THEN ERROR=.2E+01*(ABS(RES-RES3LA(2))+ABS(RES-RES3LA(1))) ELSE ERROR=ABS(RES-RES3LA(3))+ABS(RES-RES3LA(2)) 1 +ABS(RES-RES3LA(1)) ENDIF IF(ERROR.GT.1.0E+01*ABSERR) GO TO 40 ABSERR=ERROR RESULT=RES 40 CONTINUE C C COMPUTE ERROR ESTIMATE C IF(NRES.EQ.1) THEN ABSERR=.6E+01*(ABS(RESULT-RES3LA(1))) ELSE IF(NRES.EQ.2) THEN ABSERR=.2E+01*ABS(RESULT-RES3LA(2))+ABS(RESULT-RES3LA(1)) ELSE IF(NRES.GT.2) THEN ABSERR=ABS(RESULT-RES3LA(3))+ABS(RESULT-RES3LA(2)) 1 +ABS(RESULT-RES3LA(1)) ENDIF C C SHIFT THE TABLE C 50 IF(N.EQ.LIMEXP) N=2*(LIMEXP/2)-1 IB=1 IF((NUM/2)*2.EQ.NUM) IB=2 IE=NEWELM+1 DO 60 I=1,IE IB2=IB+2 EPSTAB(IB)=EPSTAB(IB2) IB=IB2 60 CONTINUE IF(NUM.EQ.N) GO TO 80 IN=NUM-N+1 DO 70 I=1,N EPSTAB(I)=EPSTAB(IN) IN=IN+1 70 CONTINUE C C UPDATE RES3LA C 80 IF(NRES.EQ.0) THEN RES3LA(1)=RESULT ELSE IF(NRES.EQ.1) THEN RES3LA(2)=RESULT ELSE IF(NRES.EQ.2) THEN RES3LA(3)=RESULT ELSE RES3LA(1)=RES3LA(2) RES3LA(2)=RES3LA(3) RES3LA(3)=RESULT ENDIF 90 ABSERR=MAX(ABSERR,EPRN*ABS(RESULT)) NRES=NRES+1 100 N=N+1 * IF(N.LE.3) ABSERR = R1MACH(2) * (0.1D-03) EPSTAB(LIMEXP+3)=REAL(N) EPSTAB(LIMEXP+4)=REAL(NRES) EPSTAB(LIMEXP+5)=RES3LA(1) EPSTAB(LIMEXP+6)=RES3LA(2) EPSTAB(LIMEXP+7)=RES3LA(3) 110 RETURN END SUBROUTINE QAGIE(F,BOUND,INF,EPSABS,EPSREL,LIMIT,RESULT,ABSERR, 1 NEVAL,IER,ALIST,BLIST,RLIST,ELIST,IORD,LAST) C***BEGIN PROLOGUE QAGIE C THIS PROLOGUE HAS BEEN REMOVED FOR REASONS OF SPACE C FOR A COMPLETE COPY OF THIS ROUTINE CONTACT THE AUTHORS C From the book "Numerical Methods and Software" C by D. Kahaner, C. Moler, S. Nash C Prentice Hall 1988 C***END PROLOGUE QAGIE C REAL ABSEPS,ABSERR,ALIST,AREA,AREA1,AREA12,AREA2,A1, 1 A2,BLIST,BOUN,BOUND,B1,B2,CORREC,DEFABS,DEFAB1,DEFAB2, 2 DRES,ELIST,EPMACH,EPSABS,EPSREL,ERLARG,ERLAST,ERRBND, 3 ERRMAX,ERROR1,ERROR2,ERRO12,ERRSUM,ERTEST,F,RESABS, 4 RESEPS,RESULT,RLIST,RLIST2,R1MACH,SMALL,UFLOW INTEGER ID,IER,IERR,IERRO,INF,IORD,IROFF1,IROFF2,IROFF3,JUPBND,K, 1 KSGN,KTMIN,LAST,LIMEXP,LIMIT,MAXERR,NEVAL,NRMAX LOGICAL EXTRAP,LERR,NEWFLG,NOEXT C PARAMETER (LIMEXP = 50) C C LIMEXP IS THE SIZE OF THE EPSILON TABLE THAT CAN BE C GENERATED IN EA C DIMENSION ALIST(LIMIT),BLIST(LIMIT),ELIST(LIMIT),IORD(LIMIT), 1 RLIST(LIMIT),RLIST2(LIMEXP+7) C EXTERNAL F C C***FIRST EXECUTABLE STATEMENT QAGIE EPMACH = R1MACH(4) C C TEST ON VALIDITY OF PARAMETERS C ----------------------------- C IER = 0 NEVAL = 0 LAST = 0 RESULT = 0.0E+00 ABSERR = 0.0E+00 ALIST(1) = 0.0E+00 BLIST(1) = 0.1E+01 RLIST(1) = 0.0E+00 ELIST(1) = 0.0E+00 IORD(1) = 0 NEWFLG = .TRUE. IF(EPSABS.LE.0.0E+00.AND.EPSREL.LT.AMAX1(0.5E+02*EPMACH,0.5E-14)) 1 IER = 6 IF(IER.EQ.6) GO TO 999 C C C FIRST APPROXIMATION TO THE INTEGRAL C ----------------------------------- C C DETERMINE THE INTERVAL TO BE MAPPED ONTO (0,1). C IF INF = 2 THE INTEGRAL IS COMPUTED AS I = I1+I2, WHERE C I1 = INTEGRAL OF F OVER (-INFINITY,0), C I2 = INTEGRAL OF F OVER (0,+INFINITY). C BOUN = BOUND IF(INF.EQ.2) BOUN = 0.0E+00 CALL QK15I(F,BOUN,INF,0.0E+00,0.1E+01,RESULT,ABSERR, 1 DEFABS,RESABS) C C TEST ON ACCURACY C LAST = 1 RLIST(1) = RESULT ELIST(1) = ABSERR IORD(1) = 1 DRES = ABS(RESULT) ERRBND = AMAX1(EPSABS,EPSREL*DRES) IF(ABSERR.LE.1.0E+02*EPMACH*DEFABS.AND.ABSERR.GT. 1 ERRBND) IER = 2 IF(LIMIT.EQ.1) IER = 1 IF(IER.NE.0.OR.(ABSERR.LE.ERRBND.AND.ABSERR.NE.RESABS).OR. 1 ABSERR.EQ.0.0E+00) GO TO 130 C C INITIALIZATION C -------------- C UFLOW = R1MACH(1) * 0.2E+01 LERR = .FALSE. CALL EA(NEWFLG,RESULT,LIMEXP,RESEPS,ABSEPS,RLIST2,IERR) ERRMAX = ABSERR MAXERR = 1 AREA = RESULT ERRSUM = ABSERR NRMAX = 1 KTMIN = 0 EXTRAP = .FALSE. NOEXT = .FALSE. IERRO = 0 IROFF1 = 0 IROFF2 = 0 IROFF3 = 0 KSGN = -1 IF(DRES.GE.(0.1E+01-0.5E+02*EPMACH)*DEFABS) KSGN = 1 C C MAIN DO-LOOP C ------------ C DO 90 LAST = 2,LIMIT C C BISECT THE SUBINTERVAL WITH NRMAX-TH LARGEST C ERROR ESTIMATE. C A1 = ALIST(MAXERR) B1 = 0.5E+00*(ALIST(MAXERR)+BLIST(MAXERR)) A2 = B1 B2 = BLIST(MAXERR) ERLAST = ERRMAX CALL QK15I(F,BOUN,INF,A1,B1,AREA1,ERROR1,RESABS,DEFAB1) CALL QK15I(F,BOUN,INF,A2,B2,AREA2,ERROR2,RESABS,DEFAB2) C C IMPROVE PREVIOUS APPROXIMATIONS TO INTEGRAL C AND ERROR AND TEST FOR ACCURACY. C AREA12 = AREA1+AREA2 ERRO12 = ERROR1+ERROR2 ERRSUM = ERRSUM+ERRO12-ERRMAX AREA = AREA+AREA12-RLIST(MAXERR) IF(DEFAB1.EQ.ERROR1.OR.DEFAB2.EQ.ERROR2)GO TO 15 IF(ABS(RLIST(MAXERR)-AREA12).GT.0.1E-04*ABS(AREA12) 1 .OR.ERRO12.LT.0.99E+00*ERRMAX) GO TO 10 IF(EXTRAP) IROFF2 = IROFF2+1 IF(.NOT.EXTRAP) IROFF1 = IROFF1+1 10 IF(LAST.GT.10.AND.ERRO12.GT.ERRMAX) IROFF3 = IROFF3+1 15 RLIST(MAXERR) = AREA1 RLIST(LAST) = AREA2 ERRBND = AMAX1(EPSABS,EPSREL*ABS(AREA)) C C TEST FOR ROUNDOFF ERROR AND EVENTUALLY C SET ERROR FLAG. C IF(IROFF1+IROFF2.GE.10.OR.IROFF3.GE.20) IER = 2 IF(IROFF2.GE.5) IERRO = 3 C C SET ERROR FLAG IN THE CASE THAT THE NUMBER OF C SUBINTERVALS EQUALS LIMIT. C IF(LAST.EQ.LIMIT) IER = 1 C C SET ERROR FLAG IN THE CASE OF BAD INTEGRAND BEHAVIOUR C AT SOME POINTS OF THE INTEGRATION RANGE. C IF(AMAX1(ABS(A1),ABS(B2)).LE.(0.1E+01+0.1E+03*EPMACH)* 1 (ABS(A2)+0.1E+04*UFLOW)) IER = 4 C C APPEND THE NEWLY-CREATED INTERVALS TO THE LIST. C IF(ERROR2.GT.ERROR1) GO TO 20 ALIST(LAST) = A2 BLIST(MAXERR) = B1 BLIST(LAST) = B2 ELIST(MAXERR) = ERROR1 ELIST(LAST) = ERROR2 GO TO 30 20 ALIST(MAXERR) = A2 ALIST(LAST) = A1 BLIST(LAST) = B1 RLIST(MAXERR) = AREA2 RLIST(LAST) = AREA1 ELIST(MAXERR) = ERROR2 ELIST(LAST) = ERROR1 C C CALL SUBROUTINE QPSRT TO MAINTAIN THE DESCENDING ORDERING C IN THE LIST OF ERROR ESTIMATES AND SELECT THE C SUBINTERVAL WITH NRMAX-TH LARGEST ERROR ESTIMATE (TO BE C BISECTED NEXT). C 30 CALL QPSRT(LIMIT,LAST,MAXERR,ERRMAX,ELIST,IORD,NRMAX) IF(ERRSUM.LE.ERRBND) GO TO 115 IF(IER.NE.0) GO TO 100 IF(LAST.EQ.2) GO TO 80 IF(NOEXT) GO TO 90 ERLARG = ERLARG-ERLAST IF(ABS(B1-A1).GT.SMALL) ERLARG = ERLARG+ERRO12 IF(EXTRAP) GO TO 40 C C TEST WHETHER THE INTERVAL TO BE BISECTED NEXT IS THE C SMALLEST INTERVAL. C IF(ABS(BLIST(MAXERR)-ALIST(MAXERR)).GT.SMALL) GO TO 90 EXTRAP = .TRUE. NRMAX = 2 40 IF(IERRO.EQ.3.OR.ERLARG.LE.ERTEST) GO TO 60 C C THE SMALLEST INTERVAL HAS THE LARGEST ERROR. C BEFORE BISECTING DECREASE THE SUM OF THE ERRORS C OVER THE LARGER INTERVALS (ERLARG) AND PERFORM C EXTRAPOLATION. C ID = NRMAX JUPBND = LAST IF(LAST.GT.(2+LIMIT/2)) JUPBND = LIMIT+3-LAST DO 50 K = ID,JUPBND MAXERR = IORD(NRMAX) ERRMAX = ELIST(MAXERR) IF(ABS(BLIST(MAXERR)-ALIST(MAXERR)).GT.SMALL) GO TO 90 NRMAX = NRMAX+1 50 CONTINUE C C PERFORM EXTRAPOLATION. C 60 CALL EA(NEWFLG,AREA,LIMEXP,RESEPS,ABSEPS,RLIST2,IERR) KTMIN = KTMIN+1 IF((KTMIN.GT.5).AND.(ABSERR.LT.0.1E-02*ERRSUM).AND.(LERR)) 1 IER = 5 IF((ABSEPS.GE.ABSERR).AND.(LERR)) GO TO 70 KTMIN = 0 ABSERR = ABSEPS LERR = .TRUE. RESULT = RESEPS CORREC = ERLARG ERTEST = AMAX1(EPSABS,EPSREL*ABS(RESEPS)) IF((ABSERR.LE.ERTEST).AND.(LERR)) GO TO 100 C C PREPARE BISECTION OF THE SMALLEST INTERVAL. C 70 IF(RLIST2(LIMEXP+3).EQ.1) NOEXT = .TRUE. IF(IER.EQ.5) GO TO 100 MAXERR = IORD(1) ERRMAX = ELIST(MAXERR) NRMAX = 1 EXTRAP = .FALSE. SMALL = SMALL*0.5E+00 ERLARG = ERRSUM GO TO 90 80 SMALL = 0.375E+00 ERLARG = ERRSUM ERTEST = ERRBND CALL EA(NEWFLG,AREA,LIMEXP,RESEPS,ABSEPS,RLIST2,IERR) 90 CONTINUE C C SET FINAL RESULT AND ERROR ESTIMATE. C ------------------------------------ C 100 IF(.NOT.LERR) GO TO 115 IF((IER+IERRO).EQ.0) GO TO 110 IF(IERRO.EQ.3) ABSERR = ABSERR+CORREC IF(IER.EQ.0) IER = 3 IF(RESULT.NE.0.0E+00.AND.AREA.NE.0.0E+00)GO TO 105 IF(ABSERR.GT.ERRSUM)GO TO 115 IF(AREA.EQ.0.0E+00) GO TO 130 GO TO 110 105 IF(ABSERR/ABS(RESULT).GT.ERRSUM/ABS(AREA))GO TO 115 C C TEST ON DIVERGENCE C 110 IF(KSGN.EQ.(-1).AND.AMAX1(ABS(RESULT),ABS(AREA)).LE. 1 DEFABS*0.1E-01) GO TO 130 IF(0.1E-01.GT.(RESULT/AREA).OR.(RESULT/AREA).GT.0.1E+03. 1OR.ERRSUM.GT.ABS(AREA)) IER = 6 GO TO 130 C C COMPUTE GLOBAL INTEGRAL SUM. C 115 RESULT = 0.0E+00 DO 120 K = 1,LAST RESULT = RESULT+RLIST(K) 120 CONTINUE ABSERR = ERRSUM 130 NEVAL = 30*LAST-15 IF(INF.EQ.2) NEVAL = 2*NEVAL IF(IER.GT.2) IER=IER-1 999 RETURN END SUBROUTINE QK15I(F,BOUN,INF,A,B,RESULT,ABSERR,RESABS,RESASC) C***BEGIN PROLOGUE QK15I C THIS PROLOGUE HAS BEEN REMOVED FOR REASONS OF SPACE C FOR A COMPLETE COPY OF THIS ROUTINE CONTACT THE AUTHORS C From the book "Numerical Methods and Software" C by D. Kahaner, C. Moler, S. Nash C Prentice Hall 1988 C***END PROLOGUE QK15I C REAL A,ABSC,ABSC1,ABSC2,ABSERR,B,BOUN,CENTR, 1 DINF,R1MACH,EPMACH,F,FC,FSUM,FVAL1,FVAL2,FV1, 2 FV2,HLGTH,RESABS,RESASC,RESG,RESK,RESKH,RESULT,TABSC1,TABSC2, 3 UFLOW,WG,WGK,XGK INTEGER INF,J,MIN0 EXTERNAL F C DIMENSION FV1(7),FV2(7),XGK(8),WGK(8),WG(8) C C DATA XGK(1),XGK(2),XGK(3),XGK(4),XGK(5),XGK(6),XGK(7), 1 XGK(8)/ 2 0.9914553711208126E+00, 0.9491079123427585E+00, 3 0.8648644233597691E+00, 0.7415311855993944E+00, 4 0.5860872354676911E+00, 0.4058451513773972E+00, 5 0.2077849550078985E+00, 0.0000000000000000E+00/ C DATA WGK(1),WGK(2),WGK(3),WGK(4),WGK(5),WGK(6),WGK(7), 1 WGK(8)/ 2 0.2293532201052922E-01, 0.6309209262997855E-01, 3 0.1047900103222502E+00, 0.1406532597155259E+00, 4 0.1690047266392679E+00, 0.1903505780647854E+00, 5 0.2044329400752989E+00, 0.2094821410847278E+00/ C DATA WG(1),WG(2),WG(3),WG(4),WG(5),WG(6),WG(7),WG(8)/ 1 0.0000000000000000E+00, 0.1294849661688697E+00, 2 0.0000000000000000E+00, 0.2797053914892767E+00, 3 0.0000000000000000E+00, 0.3818300505051189E+00, 4 0.0000000000000000E+00, 0.4179591836734694E+00/ C C C***FIRST EXECUTABLE STATEMENT QK15I EPMACH = R1MACH(4) UFLOW = R1MACH(1) DINF = MIN0(1,INF) C CENTR = 0.5E+00*(A+B) HLGTH = 0.5E+00*(B-A) TABSC1 = BOUN+DINF*(0.1E+01-CENTR)/CENTR FVAL1 = F(TABSC1) IF(INF.EQ.2) FVAL1 = FVAL1+F(-TABSC1) FC = (FVAL1/CENTR)/CENTR C C COMPUTE THE 15-POINT KRONROD APPROXIMATION TO C THE INTEGRAL, AND ESTIMATE THE ERROR. C RESG = WG(8)*FC RESK = WGK(8)*FC RESABS = ABS(RESK) DO 10 J=1,7 ABSC = HLGTH*XGK(J) ABSC1 = CENTR-ABSC ABSC2 = CENTR+ABSC TABSC1 = BOUN+DINF*(0.1E+01-ABSC1)/ABSC1 TABSC2 = BOUN+DINF*(0.1E+01-ABSC2)/ABSC2 FVAL1 = F(TABSC1) FVAL2 = F(TABSC2) IF(INF.EQ.2) FVAL1 = FVAL1+F(-TABSC1) IF(INF.EQ.2) FVAL2 = FVAL2+F(-TABSC2) FVAL1 = (FVAL1/ABSC1)/ABSC1 FVAL2 = (FVAL2/ABSC2)/ABSC2 FV1(J) = FVAL1 FV2(J) = FVAL2 FSUM = FVAL1+FVAL2 RESG = RESG+WG(J)*FSUM RESK = RESK+WGK(J)*FSUM RESABS = RESABS+WGK(J)*(ABS(FVAL1)+ABS(FVAL2)) 10 CONTINUE RESKH = RESK*0.5E+00 RESASC = WGK(8)*ABS(FC-RESKH) DO 20 J=1,7 RESASC = RESASC+WGK(J)*(ABS(FV1(J)-RESKH)+ABS(FV2(J)-RESKH)) 20 CONTINUE RESULT = RESK*HLGTH RESASC = RESASC*HLGTH RESABS = RESABS*HLGTH ABSERR = ABS((RESK-RESG)*HLGTH) IF(RESASC.NE.0.0E+00.AND.ABSERR.NE.0.E0) ABSERR = RESASC* 1 AMIN1(0.1E+01,(0.2E+03*ABSERR/RESASC)**1.5E+00) IF(RESABS.GT.UFLOW/(0.5E+02*EPMACH)) ABSERR = AMAX1 1 ((EPMACH*0.5E+02)*RESABS,ABSERR) RETURN END SUBROUTINE QPSRT(LIMIT,LAST,MAXERR,ERMAX,ELIST,IORD,NRMAX) C***BEGIN PROLOGUE QPSRT C THIS PROLOGUE HAS BEEN REMOVED FOR REASONS OF SPACE C FOR A COMPLETE COPY OF THIS ROUTINE CONTACT THE AUTHORS C From the book "Numerical Methods and Software" C by D. Kahaner, C. Moler, S. Nash C Prentice Hall 1988 C***END PROLOGUE QPSRT C REAL ELIST,ERMAX,ERRMAX,ERRMIN INTEGER I,IBEG,IDO,IORD,ISUCC,J,JBND,JUPBN,K,LAST,LIMIT,MAXERR, 1 NRMAX DIMENSION ELIST(LAST),IORD(LAST) C C CHECK WHETHER THE LIST CONTAINS MORE THAN C TWO ERROR ESTIMATES. C C***FIRST EXECUTABLE STATEMENT QPSRT IF(LAST.GT.2) GO TO 10 IORD(1) = 1 IORD(2) = 2 GO TO 90 C C THIS PART OF THE ROUTINE IS ONLY EXECUTED C IF, DUE TO A DIFFICULT INTEGRAND, SUBDIVISION C INCREASED THE ERROR ESTIMATE. IN THE NORMAL CASE C THE INSERT PROCEDURE SHOULD START AFTER THE C NRMAX-TH LARGEST ERROR ESTIMATE. C 10 ERRMAX = ELIST(MAXERR) IF(NRMAX.EQ.1) GO TO 30 IDO = NRMAX-1 DO 20 I = 1,IDO ISUCC = IORD(NRMAX-1) C ***JUMP OUT OF DO-LOOP IF(ERRMAX.LE.ELIST(ISUCC)) GO TO 30 IORD(NRMAX) = ISUCC NRMAX = NRMAX-1 20 CONTINUE C C COMPUTE THE NUMBER OF ELEMENTS IN THE LIST TO C BE MAINTAINED IN DESCENDING ORDER. THIS NUMBER C DEPENDS ON THE NUMBER OF SUBDIVISIONS STILL C ALLOWED. C 30 JUPBN = LAST IF(LAST.GT.(LIMIT/2+2)) JUPBN = LIMIT+3-LAST ERRMIN = ELIST(LAST) C C INSERT ERRMAX BY TRAVERSING THE LIST TOP-DOWN, C STARTING COMPARISON FROM THE ELEMENT ELIST(IORD(NRMAX+1)). C JBND = JUPBN-1 IBEG = NRMAX+1 IF(IBEG.GT.JBND) GO TO 50 DO 40 I=IBEG,JBND ISUCC = IORD(I) C ***JUMP OUT OF DO-LOOP IF(ERRMAX.GE.ELIST(ISUCC)) GO TO 60 IORD(I-1) = ISUCC 40 CONTINUE 50 IORD(JBND) = MAXERR IORD(JUPBN) = LAST GO TO 90 C C INSERT ERRMIN BY TRAVERSING THE LIST BOTTOM-UP. C 60 IORD(I-1) = MAXERR K = JBND DO 70 J=I,JBND ISUCC = IORD(K) C ***JUMP OUT OF DO-LOOP IF(ERRMIN.LT.ELIST(ISUCC)) GO TO 80 IORD(K+1) = ISUCC K = K-1 70 CONTINUE IORD(I) = LAST GO TO 90 80 IORD(K+1) = LAST C C SET MAXERR AND ERMAX. C 90 MAXERR = IORD(NRMAX) ERMAX = ELIST(MAXERR) RETURN END