SUBROUTINE DQAGI(F,BOUND,INF,EPSABS,EPSREL,RESULT,ABSERR,NEVAL, 1 IER,LIMIT,LENW,LAST,IWORK,WORK) C***BEGIN PROLOGUE DQAGI 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 Integration over infinite intervals C Standard fortran subroutine C C PARAMETERS C ON ENTRY C F - Double precision 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 - Double precision 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 - Double precision C Absolute accuracy requested C EPSREL - Double precision 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 - Double precision C Approximation to the integral C C ABSERR - Double precision 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. 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 - Double precision 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 DQAGIE,XERROR C***END PROLOGUE DQAGI C DOUBLE PRECISION 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 DQAGI IER = 6 NEVAL = 0 LAST = 0 RESULT = 0.0D+00 ABSERR = 0.0D+00 IF(LIMIT.LT.1.OR.LENW.LT.LIMIT*4) GO TO 10 C C PREPARE CALL FOR DQAGIE. C L1 = LIMIT+1 L2 = LIMIT+L1 L3 = LIMIT+L2 C CALL DQAGIE(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 DQAGI',26,IER,LVL) RETURN END SUBROUTINE DEA(NEWFLG,SVALUE,LIMEXP,RESULT,ABSERR,EPSTAB,IERR) C***BEGIN PROLOGUE DEA C***DATE WRITTEN 800101 (YYMMDD) C***REVISION DATE 880908 (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 DONKER-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***DESCRIPTION C C Epsilon algorithm. Standard fortran subroutine. C Double precision version. C C A R G U M E N T S I N T H E C A L L S E Q U E N C E C C NEWFLG - LOGICAL (INPUT and OUTPUT) C On the first call to DEA set NEWFLG to .TRUE. C (indicating a new sequence). DEA will set NEWFLG C to .FALSE. C C SVALUE - DOUBLE PRECISION (INPUT) C On the first call to DEA set SVALUE to the first C term in the sequence. On subsequent calls set C SVALUE to the subsequent sequence value. C C LIMEXP - INTEGER (INPUT) C An integer equal to or greater than the total C number of sequence terms to be evaluated. Do not C change the value of LIMEXP until a new sequence C is evaluated (NEWFLG=.TRUE.). LIMEXP .GE. 3 C C RESULT - DOUBLE PRECISION (OUTPUT) C Best approximation to the sequence's limit. C C ABSERR - DOUBLE PRECISION (OUTPUT) C Estimate of the absolute error. C C EPSTAB - DOUBLE PRECISION (OUTPUT) C Workvector of DIMENSION at least (LIMEXP+7). C C IERR - INTEGER (OUTPUT) C IERR=0 Normal termination of the routine. C IERR=1 The input is invalid because LIMEXP.LT.3. C C T Y P I C A L P R O B L E M S E T U P C C This sample problem uses the trapezoidal rule to evaluate the C integral of the sin function from 0.0 to 0.5*PI (value = 1.0). The C program implements the trapezoidal rule 8 times creating an C increasingly accurate sequence of approximations to the integral. C Each time the trapezoidal rule is used, it uses twice as many C panels as the time before. DEA is called to obtain even more C accurate estimates. C C PROGRAM SAMPLE C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C DOUBLE PRECISION EPSTAB(57) CC [57 = LIMEXP + 7] C LOGICAL NEWFLG C EXTERNAL F C DATA LIMEXP/50/ C WRITE(*,*) ' NO. PANELS TRAP. APPROX' C * ,' APPROX W/EA ABSERR' C WRITE(*,*) C HALFPI = DASIN(1.0D+00) CC [UPPER INTEGRATION LIMIT = PI/2] C NEWFLG = .TRUE. CC [SET FLAG - 1ST DEA CALL] C DO 10 I = 0,7 C NPARTS = 2 ** I C WIDTH = HALFPI/NPARTS C APPROX = 0.5D+00 * WIDTH * (F(0.0D+00) + F(HALFPI)) C DO 11 J = 1,NPARTS-1 C APPROX = APPROX + F(J * WIDTH) * WIDTH C 11 CONTINUE CC [END TRAPEZOIDAL RULE APPROX] C SVALUE = APPROX CC [SVALUE = NEW SEQUENCE VALUE] C CALL DEA(NEWFLG,SVALUE,LIMEXP,RESULT,ABSERR,EPSTAB,IERR) CC [CALL DEA FOR BETTER ESTIMATE] C WRITE(*,12) NPARTS,APPROX,RESULT,ABSERR C 12 FORMAT(' ',I4,T20,F16.13,T40,F16.13,T60,D11.4) C 10 CONTINUE C STOP C END C C DOUBLE PRECISION FUNCTION F(X) C DOUBLE PRECISION X C F = DSIN(X) CC [INTEGRAND] C RETURN C END C C Output from the above program will be: C C NO. PANELS TRAP. APPROX APPROX W/EA ABSERR C C 1 .7853981633974 .7853981633974 .7854D+00 C 2 .9480594489685 .9480594489685 .9760D+00 C 4 .9871158009728 .9994567212570 .2141D+00 C 8 .9967851718862 .9999667417647 .5100D-03 C 16 .9991966804851 .9999998781041 .5763D-03 C 32 .9997991943200 .9999999981026 .5767D-03 C 64 .9999498000921 .9999999999982 .3338D-04 C 128 .9999874501175 1.0000000000000 .1238D-06 C C----------------------------------------------------------------------- C***REFERENCES "Acceleration de la convergence en analyse numerique", C C. Brezinski, "Lecture Notes in Math.", vol. 584, C Springer-Verlag, New York, 1977. C***ROUTINES CALLED D1MACH,XERROR C***END PROLOGUE DEA DOUBLE PRECISION ABSERR,DELTA1,DELTA2,DELTA3,DRELPR,DEPRN, 1 D1MACH,EPSTAB(*),ERROR,ERR1,ERR2,ERR3,E0,E1,E2,E3,RES, 2 RESULT,RES3LA(3),SS,SVALUE,TOL1,TOL2,TOL3 INTEGER I,IB,IB2,IE,IERR,IN,K1,K2,K3,LIMEXP,N,NEWELM,NUM,NRES LOGICAL NEWFLG C C C LIMEXP is the maximum number of elements the C epsilon table data can contain. The epsilon table C is stored in the first (LIMEXP+2) entries of EPSTAB. C C C LIST OF MAJOR VARIABLES C ----------------------- C E0,E1,E2,E3 - DOUBLE PRECISION C The 4 elements on which the computation of C a new element in the epsilon table is based. C NRES - INTEGER C Number of extrapolation results actually C generated by the epsilon algorithm in prior C calls to the routine. C NEWELM - INTEGER C Number of elements to be computed in the C new diagonal of the epsilon table. The C condensed epsilon table is computed. Only C those elements needed for the computation of C the next diagonal are preserved. C RES - DOUBLE PRECISION C New element in the new diagonal of the C epsilon table. C ERROR - DOUBLE PRECISION C An estimate of the absolute error of RES. C Routine decides whether RESULT=RES or C RESULT=SVALUE by comparing ERROR with C ABSERR from the previous call. C RES3LA - DOUBLE PRECISION C Vector of DIMENSION 3 containing at most C the last 3 results. C C C MACHINE DEPENDENT CONSTANTS C --------------------------- C DRELPR is the largest relative spacing. C C***FIRST EXECUTABLE STATEMENT DEA 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=.6D+01*ABS(RESULT-EPSTAB(1)) GO TO 100 ENDIF ENDIF EPSTAB(N)=SVALUE DRELPR=D1MACH(4) DEPRN=1.0D+01*DRELPR 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))*DRELPR DELTA3=E1-E0 ERR3=ABS(DELTA3) TOL3=MAX(ABS(E1),ABS(E0))*DRELPR 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))*DRELPR 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.1D+01/DELTA1+0.1D+01/DELTA2-0.1D+01/DELTA3 ELSE EPSTAB(K1)=E1 IF(ERR2.LE.TOL2.OR.ERR3.LE.TOL3) GO TO 20 SS=0.1D+01/DELTA2-0.1D+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.1D-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.1D+01/SS EPSTAB(K1)=RES K1=K1-2 ERROR=ERR2+ABS(RES-E2)+ERR3 IF(NRES.EQ.0) GO TO 35 IF(ERROR.GT.0.1E3*ABSERR) GO TO 40 35 ABSERR=ERROR RESULT=RES 40 CONTINUE C C COMPUTE ERROR ESTIMATE C IF(NRES.EQ.1) THEN ABSERR=ABS(RESULT-RES3LA(1)) ELSE IF(NRES.EQ.2) THEN ABSERR=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,DEPRN*ABS(RESULT)) NRES=NRES+1 100 N=N+1 EPSTAB(LIMEXP+3)=DBLE(N) EPSTAB(LIMEXP+4)=DBLE(NRES) EPSTAB(LIMEXP+5)=RES3LA(1) EPSTAB(LIMEXP+6)=RES3LA(2) EPSTAB(LIMEXP+7)=RES3LA(3) 110 RETURN END SUBROUTINE DQPSRT(LIMIT,LAST,MAXERR,ERMAX,ELIST,IORD,NRMAX) C***BEGIN PROLOGUE DQPSRT C***REFER TO DQAGE,DQAGIE,DQAGPE,DQAWSE C***ROUTINES CALLED (NONE) C***REVISION DATE 810101 (YYMMDD) C***KEYWORDS SEQUENTIAL SORTING 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 This routine maintains the descending ordering in the C list of the local error estimated resulting from the C interval subdivision process. At each call two error C estimates are inserted using the sequential search C method, top-down for the largest error estimate and C bottom-up for the smallest error estimate. C***DESCRIPTION C C Ordering routine C Standard fortran subroutine C Double precision version C C PARAMETERS (MEANING AT OUTPUT) C LIMIT - Integer C Maximum number of error estimates the list C can contain C C LAST - Integer C Number of error estimates currently in the list C C MAXERR - Integer C Maxerr points to the NRMAX-th largest error C estimate currently in the list C C ERMAX - Double precision C NRMAX-th largest error estimate C ERMAX = ELIST(MAXERR) C C ELIST - Double precision C Vector of dimension LAST containing C the error estimates C C IORD - Integer C Vector of dimension LAST, the first K elements C of which contain pointers to the error C estimates, such that C ELIST(IORD(1)),..., ELIST(IORD(K)) C form a decreasing sequence, with C K = LAST if LAST.LE.(LIMIT/2+2), and C K = LIMIT+1-LAST otherwise C C NRMAX - Integer C MAXERR = IORD(NRMAX) C***END PROLOGUE DQPSRT C DOUBLE PRECISION 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 DQPSRT 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 IF, DUE TO A C DIFFICULT INTEGRAND, SUBDIVISION INCREASED THE ERROR C ESTIMATE. IN THE NORMAL CASE THE INSERT PROCEDURE SHOULD C START AFTER THE 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 BE MAINTAINED C IN DESCENDING ORDER. THIS NUMBER DEPENDS ON THE NUMBER OF C SUBDIVISIONS STILL 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 SUBROUTINE DQAGIE(F,BOUND,INF,EPSABS,EPSREL,LIMIT,RESULT,ABSERR, 1 NEVAL,IER,ALIST,BLIST,RLIST,ELIST,IORD,LAST) C***BEGIN PROLOGUE DQAGIE 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 Integration over infinite intervals C Standard fortran subroutine C C F - Double precision 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 - Double precision C Finite bound of integration range C (has no meaning if interval is doubly-infinite) C C INF - Double precision 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 - Double precision C Absolute accuracy requested C EPSREL - Double precision 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 LIMIT - Integer C Gives an upper bound on the number of subintervals C in the partition of (A,B), LIMIT.GE.1 C C ON RETURN C RESULT - Double precision C Approximation to the integral C C ABSERR - Double precision 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. C If the position of a local difficulty can C be 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 RESULT, ABSERR, NEVAL, LAST, RLIST(1), C ELIST(1) and IORD(1) are set to zero. C ALIST(1) and BLIST(1) are set to 0 C and 1 respectively. C C ALIST - Double precision C Vector of dimension at least LIMIT, the first C LAST elements of which are the left C end points of the subintervals in the partition C of the transformed integration range (0,1). C C BLIST - Double precision C Vector of dimension at least LIMIT, the first C LAST elements of which are the right C end points of the subintervals in the partition C of the transformed integration range (0,1). C C RLIST - Double precision C Vector of dimension at least LIMIT, the first C LAST elements of which are the integral C approximations on the subintervals C C ELIST - Double precision C Vector of dimension at least LIMIT, the first C LAST elements of which are the moduli of the C absolute error estimates on the subintervals C C IORD - Integer C Vector of dimension LIMIT, the first K C elements of which are pointers to the C error estimates over the subintervals, C such that ELIST(IORD(1)), ..., ELIST(IORD(K)) C form a decreasing sequence, with K = LAST C If LAST.LE.(LIMIT/2+2), and K = LIMIT+1-LAST C otherwise C C LAST - Integer C Number of subintervals actually produced C in the subdivision process C***REFERENCES (NONE) C***ROUTINES CALLED DEA,DQK15I,DQPSRT,D1MACH C***END PROLOGUE DQAGIE DOUBLE PRECISION ABSEPS,ABSERR,ALIST,AREA,AREA1,AREA12,AREA2,A1, 1 A2,BLIST,BOUN,BOUND,B1,B2,CORREC,DABS,DEFABS,DEFAB1,DEFAB2, 2 DMAX1,DRES,D1MACH,ELIST,EPMACH,EPSABS,EPSREL,ERLARG,ERLAST, 3 ERRBND,ERRMAX,ERROR1,ERROR2,ERRO12,ERRSUM,ERTEST,F,RESABS, 4 RESEPS,RESULT,RLIST,RLIST2,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 LIST OF MAJOR VARIABLES C ----------------------- C C ALIST - LIST OF LEFT END POINTS OF ALL SUBINTERVALS C CONSIDERED UP TO NOW C BLIST - LIST OF RIGHT END POINTS OF ALL SUBINTERVALS C CONSIDERED UP TO NOW C RLIST(I) - APPROXIMATION TO THE INTEGRAL OVER C (ALIST(I),BLIST(I)) C RLIST2 - ARRAY OF DIMENSION AT LEAST (LIMEXP+7), C CONTAINING THE PART OF THE EPSILON TABLE C WICH IS STILL NEEDED FOR FURTHER COMPUTATIONS C ELIST(I) - ERROR ESTIMATE APPLYING TO RLIST(I) C MAXERR - POINTER TO THE INTERVAL WITH LARGEST ERROR C ESTIMATE C ERRMAX - ELIST(MAXERR) C ERLAST - ERROR ON THE INTERVAL CURRENTLY SUBDIVIDED C (BEFORE THAT SUBDIVISION HAS TAKEN PLACE) C AREA - SUM OF THE INTEGRALS OVER THE SUBINTERVALS C ERRSUM - SUM OF THE ERRORS OVER THE SUBINTERVALS C ERRBND - REQUESTED ACCURACY MAX(EPSABS,EPSREL* C ABS(RESULT)) C *****1 - VARIABLE FOR THE LEFT SUBINTERVAL C *****2 - VARIABLE FOR THE RIGHT SUBINTERVAL C LAST - INDEX FOR SUBDIVISION C SMALL - LENGTH OF THE SMALLEST INTERVAL CONSIDERED UP C TO NOW, MULTIPLIED BY 1.5 C ERLARG - SUM OF THE ERRORS OVER THE INTERVALS LARGER C THAN THE SMALLEST INTERVAL CONSIDERED UP TO NOW C EXTRAP - LOGICAL VARIABLE DENOTING THAT THE ROUTINE C IS ATTEMPTING TO PERFORM EXTRAPOLATION. I.E. C BEFORE SUBDIVIDING THE SMALLEST INTERVAL WE C TRY TO DECREASE THE VALUE OF ERLARG. C LERR - LOGICAL VARIABLE INDICATING WHETHER ABSERR HAS C BEEN CHANGED IN MAIN DO-LOOP. C NOEXT - LOGICAL VARIABLE DENOTING THAT EXTRAPOLATION C IS NO LONGER ALLOWED (TRUE-VALUE) C C MACHINE DEPENDENT CONSTANTS C --------------------------- C C EPMACH IS THE LARGEST RELATIVE SPACING. C UFLOW IS THE SMALLEST POSITIVE MAGNITUDE. C C***FIRST EXECUTABLE STATEMENT DQAGIE EPMACH = D1MACH(4) C C TEST ON VALIDITY OF PARAMETERS C ----------------------------- C IER = 0 NEVAL = 0 LAST = 0 RESULT = 0.0D+00 ABSERR = 0.0D+00 ALIST(1) = 0.0D+00 BLIST(1) = 0.1D+01 RLIST(1) = 0.0D+00 ELIST(1) = 0.0D+00 IORD(1) = 0 NEWFLG = .TRUE. IF(EPSABS.LE.0.0D+00.AND.EPSREL.LT.DMAX1(0.5D+02*EPMACH,0.5D-28)) 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.0D+00 CALL DQK15I(F,BOUN,INF,0.0D+00,0.1D+01,RESULT,ABSERR, 1 DEFABS,RESABS) C C TEST ON ACCURACY C LAST = 1 RLIST(1) = RESULT ELIST(1) = ABSERR IORD(1) = 1 DRES = DABS(RESULT) ERRBND = DMAX1(EPSABS,EPSREL*DRES) IF(ABSERR.LE.1.0D+02*EPMACH*DEFABS.AND.ABSERR.GT.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.0D+00) GO TO 130 C C INITIALIZATION C -------------- C UFLOW = D1MACH(1) * 0.2E+01 LERR = .FALSE. CALL DEA(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.1D+01-0.5D+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 ERROR ESTIMATE. C A1 = ALIST(MAXERR) B1 = 0.5D+00*(ALIST(MAXERR)+BLIST(MAXERR)) A2 = B1 B2 = BLIST(MAXERR) ERLAST = ERRMAX CALL DQK15I(F,BOUN,INF,A1,B1,AREA1,ERROR1,RESABS,DEFAB1) CALL DQK15I(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(DABS(RLIST(MAXERR)-AREA12).GT.0.1D-04*DABS(AREA12) 1 .OR.ERRO12.LT.0.99D+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 = DMAX1(EPSABS,EPSREL*DABS(AREA)) C C TEST FOR ROUNDOFF ERROR AND EVENTUALLY 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(DMAX1(DABS(A1),DABS(B2)).LE.(0.1D+01+0.1D+03*EPMACH)* 1 (DABS(A2)+0.1D+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 DQPSRT TO MAINTAIN THE DESCENDING ORDERING C IN THE LIST OF ERROR ESTIMATES AND SELECT THE SUBINTERVAL C WITH NRMAX-TH LARGEST ERROR ESTIMATE (TO BE BISECTED NEXT). C 30 CALL DQPSRT(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(DABS(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(DABS(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 OVER THE C LARGER INTERVALS (ERLARG) AND PERFORM 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(DABS(BLIST(MAXERR)-ALIST(MAXERR)).GT.SMALL) GO TO 90 NRMAX = NRMAX+1 50 CONTINUE C C PERFORM EXTRAPOLATION. C 60 CALL DEA(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 = DMAX1(EPSABS,EPSREL*DABS(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.5D+00 ERLARG = ERRSUM GO TO 90 80 SMALL = 0.375D+00 ERLARG = ERRSUM ERTEST = ERRBND CALL DEA(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.0D+00.AND.AREA.NE.0.0D+00)GO TO 105 IF(ABSERR.GT.ERRSUM)GO TO 115 IF(AREA.EQ.0.0D+00) GO TO 130 GO TO 110 105 IF(ABSERR/DABS(RESULT).GT.ERRSUM/DABS(AREA))GO TO 115 C C TEST ON DIVERGENCE C 110 IF(KSGN.EQ.(-1).AND.DMAX1(DABS(RESULT),DABS(AREA)).LE. 1 DEFABS*0.1D-01) GO TO 130 IF(0.1D-01.GT.(RESULT/AREA).OR.(RESULT/AREA).GT.0.1D+03. 1OR.ERRSUM.GT.DABS(AREA)) IER = 6 GO TO 130 C C COMPUTE GLOBAL INTEGRAL SUM. C 115 RESULT = 0.0D+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 DQK15I(F,BOUN,INF,A,B,RESULT,ABSERR,RESABS,RESASC) C***BEGIN PROLOGUE DQK15I C***DATE WRITTEN 800101 (YYMMDD) C***REVISION DATE 830518 (YYMMDD) C***CATEGORY NO. H2A3A2,H2A4A2 C***KEYWORDS 15-POINT TRANSFORMED GAUSS-KRONROD RULES 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 original (infinite integration range is mapped C onto the interval (0,1) and (A,B) is a part of (0,1). C it is the purpose to compute C I = Integral of transformed integrand over (A,B), C J = Integral of ABS(Transformed Integrand) over (A,B). C***DESCRIPTION C C Integration Rule C Standard Fortran subroutine C Double precision version C C PARAMETERS C ON ENTRY C F - Double precision C Fuction 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 calling program. C C BOUN - Double precision C Finite bound of original integration C Range (SET TO ZERO IF INF = +2) C C INF - Integer C If INF = -1, the original interval is C (-INFINITY,BOUND), C If INF = +1, the original interval is C (BOUND,+INFINITY), C If INF = +2, the original interval is C (-INFINITY,+INFINITY) AND C The integral is computed as the sum of two C integrals, one over (-INFINITY,0) and one over C (0,+INFINITY). C C A - Double precision C Lower limit for integration over subrange C of (0,1) C C B - Double precision C Upper limit for integration over subrange C of (0,1) C C ON RETURN C RESULT - Double precision C Approximation to the integral I C Result is computed by applying the 15-POINT C KRONROD RULE(RESK) obtained by optimal addition C of abscissae to the 7-POINT GAUSS RULE(RESG). C C ABSERR - Double precision C Estimate of the modulus of the absolute error, C WHICH SHOULD EQUAL or EXCEED ABS(I-RESULT) C C RESABS - Double precision C Approximation to the integral J C C RESASC - Double precision C Approximation to the integral of C ABS((TRANSFORMED INTEGRAND)-I/(B-A)) over (A,B) C***REFERENCES (NONE) C***ROUTINES CALLED D1MACH C***END PROLOGUE DQK15I C DOUBLE PRECISION A,ABSC,ABSC1,ABSC2,ABSERR,B,BOUN,CENTR,DABS,DINF, 1 DMAX1,DMIN1,D1MACH,EPMACH,F,FC,FSUM,FVAL1,FVAL2,FV1,FV2,HLGTH, 2 RESABS,RESASC,RESG,RESK,RESKH,RESULT,TABSC1,TABSC2,UFLOW,WG,WGK, 3 XGK INTEGER INF,J EXTERNAL F C DIMENSION FV1(7),FV2(7),XGK(8),WGK(8),WG(8) C C THE ABSCISSAE AND WEIGHTS ARE SUPPLIED FOR THE INTERVAL C (-1,1). BECAUSE OF SYMMETRY ONLY THE POSITIVE ABSCISSAE AND C THEIR CORRESPONDING WEIGHTS ARE GIVEN. C C XGK - ABSCISSAE OF THE 15-POINT KRONROD RULE C XGK(2), XGK(4), ... ABSCISSAE OF THE 7-POINT C GAUSS RULE C XGK(1), XGK(3), ... ABSCISSAE WHICH ARE OPTIMALLY C ADDED TO THE 7-POINT GAUSS RULE C C WGK - WEIGHTS OF THE 15-POINT KRONROD RULE C C WG - WEIGHTS OF THE 7-POINT GAUSS RULE, CORRESPONDING C TO THE ABSCISSAE XGK(2), XGK(4), ... C WG(1), WG(3), ... ARE SET TO ZERO. C DATA XGK(1),XGK(2),XGK(3),XGK(4),XGK(5),XGK(6),XGK(7),XGK(8)/ 1 0.9914553711208126D+00, 0.9491079123427585D+00, 2 0.8648644233597691D+00, 0.7415311855993944D+00, 3 0.5860872354676911D+00, 0.4058451513773972D+00, 4 0.2077849550078985D+00, 0.0000000000000000D+00/ C DATA WGK(1),WGK(2),WGK(3),WGK(4),WGK(5),WGK(6),WGK(7),WGK(8)/ 1 0.2293532201052922D-01, 0.6309209262997855D-01, 2 0.1047900103222502D+00, 0.1406532597155259D+00, 3 0.1690047266392679D+00, 0.1903505780647854D+00, 4 0.2044329400752989D+00, 0.2094821410847278D+00/ C DATA WG(1),WG(2),WG(3),WG(4),WG(5),WG(6),WG(7),WG(8)/ 1 0.0000000000000000D+00, 0.1294849661688697D+00, 2 0.0000000000000000D+00, 0.2797053914892767D+00, 3 0.0000000000000000D+00, 0.3818300505051189D+00, 4 0.0000000000000000D+00, 0.4179591836734694D+00/ C C C LIST OF MAJOR VARIABLES C ----------------------- C C CENTR - MID POINT OF THE INTERVAL C HLGTH - HALF-LENGTH OF THE INTERVAL C ABSC* - ABSCISSA C TABSC* - TRANSFORMED ABSCISSA C FVAL* - FUNCTION VALUE C RESG - RESULT OF THE 7-POINT GAUSS FORMULA C RESK - RESULT OF THE 15-POINT KRONROD FORMULA C RESKH - APPROXIMATION TO THE MEAN VALUE OF THE TRANSFORMED C INTEGRAND OVER (A,B), I.E. TO I/(B-A) C C MACHINE DEPENDENT CONSTANTS C --------------------------- C C EPMACH IS THE LARGEST RELATIVE SPACING. C UFLOW IS THE SMALLEST POSITIVE MAGNITUDE. C C***FIRST EXECUTABLE STATEMENT DQK15I EPMACH = D1MACH(4) UFLOW = D1MACH(1) DINF = MIN0(1,INF) C CENTR = 0.5D+00*(A+B) HLGTH = 0.5D+00*(B-A) TABSC1 = BOUN+DINF*(0.1D+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 = DABS(RESK) DO 10 J=1,7 ABSC = HLGTH*XGK(J) ABSC1 = CENTR-ABSC ABSC2 = CENTR+ABSC TABSC1 = BOUN+DINF*(0.1D+01-ABSC1)/ABSC1 TABSC2 = BOUN+DINF*(0.1D+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)*(DABS(FVAL1)+DABS(FVAL2)) 10 CONTINUE RESKH = RESK*0.5D+00 RESASC = WGK(8)*DABS(FC-RESKH) DO 20 J=1,7 RESASC = RESASC+WGK(J)*(DABS(FV1(J)-RESKH)+DABS(FV2(J)-RESKH)) 20 CONTINUE RESULT = RESK*HLGTH RESASC = RESASC*HLGTH RESABS = RESABS*HLGTH ABSERR = DABS((RESK-RESG)*HLGTH) IF(RESASC.NE.0.0D+00.AND.ABSERR.NE.0.D0) ABSERR = RESASC* 1 DMIN1(0.1D+01,(0.2D+03*ABSERR/RESASC)**1.5D+00) IF(RESABS.GT.UFLOW/(0.5D+02*EPMACH)) ABSERR = DMAX1 1 ((EPMACH*0.5D+02)*RESABS,ABSERR) RETURN END