C ALGORITHM 691, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 17, NO. 2, PP. 218-232. JUNE, 1991. This disk contains 8 files: README : this information file stest.for : which contains simple test routines (single precision) TEST (main) F squad.for : which contains quadrature routines (single precision) QXGS QXGSE QXG QXGE QXLQM QXRUL QXRRD QXCPY saux.for : which contains auxiliary routines (single precision) QPSRT QELG R1MACH (machine dependent, TO BE MODIFIED) dtest.for : which contains simple test routines (double precision) DTEST (main) F dquad.for : which contains quadrature routines (double precision) DQXGS DQXGSE DQXG DQXGE DQXLQM DQXRUL DQXRRD DQXCPY daux.for : which contains auxiliary routines (double precision) DQPSRT DQELG D1MACH (machine dependent, TO BE MODIFIED) CC+------------------------------------------------------------------+CC CC+ +CC CC+ LAWRENCE LIVERMORE NATIONAL LABORATORY +CC CC+ +CC CC+ MATHEMATICS AND STATISTICS DIVISION SOFTWARE LIBRARY +CC CC+ +CC CC+------------------------------------------------------------------+CC CC+ +CC CC+ THE SLATEC MATHEMATICAL PROGRAM LIBRARY IS ISSUED BY +CC CC+ THE FOLLOWING: +CC CC+ +CC CC+ AIR FORCE WEAPONS LABORATORY, ALBUQUERQUE +CC CC+ LAWRENCE LIVERMORE NATIONAL LABORATORY, LIVERMORE +CC CC+ LOS ALAMOS NATIONAL SCIENTIFIC LABORATORY, LOS ALAMOS +CC CC+ NATIONAL BUREAU OF STANDARDS, WASHINGTON +CC CC+ OAK RIDGE NATIONAL LABORATORY, OAK RIDGE +CC CC+ SANDIA NATIONAL LABORATORIES, ALBUQUERQUE +CC CC+ SANDIA NATIONAL LABORATORIES, LIVERMORE +CC CC+ +CC CC+ PLEASE REFER QUESTIONS REGARDING THIS SOFTWARE TO THE +CC CC+ MSD CONSULTING OFFICE, EXTENSION 3-2976. +CC CC+ +CC CC+ +CC CC+ * * * * * N O T I C E * * * * * +CC CC+ +CC CC+ THIS MATERIAL WAS PREPARED AS AN ACCOUNT OF WORK SPONSORED +CC CC+ BY THE UNITED STATES GOVERNMENT. NEITHER THE UNITED +CC CC+ STATES GOVERNMENT, NOR THE DEPARTMENT OF DEFENSE, NOR +CC CC+ ANY OF THEIR EMPLOYEES, MAKES ANY WARRANTY, EXPRESSED +CC CC+ OR IMPLIED, OR ASSUMES ANY LEGAL LIABILITY OR +CC CC+ RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR +CC CC+ USEFULNESS OF ANY INFORMATION, APPARATUS, PRODUCT, +CC CC+ OR PROCESS DISCLOSED, OR REPRESENTS ITS USE WOULD NOT +CC CC+ INFRINGE UPON PRIVATELY OWNED RIGHTS. +CC CC+ +CC CC+------------------------------------------------------------------+CC *DECK DQPSRT 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 *DECK DQELG SUBROUTINE DQELG(N,EPSTAB,RESULT,ABSERR,RES3LA,NRES) C***BEGIN PROLOGUE DQELG C***REFER TO DQAGIE,DQAGOE,DQAGPE,DQAGSE C***ROUTINES CALLED D1MACH C***REVISION DATE 830518 (YYMMDD) C***KEYWORDS CONVERGENCE ACCELERATION,EPSILON ALGORITHM,EXTRAPOLATION 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 DETERMINES THE LIMIT OF A GIVEN SEQUENCE OF C APPROXIMATIONS, BY MEANS OF THE EPSILON ALGORITHM OF C P.WYNN. AN ESTIMATE OF THE ABSOLUTE ERROR IS ALSO GIVEN. C THE CONDENSED EPSILON TABLE IS COMPUTED. ONLY THOSE C ELEMENTS NEEDED FOR THE COMPUTATION OF THE NEXT DIAGONAL C ARE PRESERVED. C***DESCRIPTION C C EPSILON ALGORITHM C STANDARD FORTRAN SUBROUTINE C DOUBLE PRECISION VERSION C C PARAMETERS C N - INTEGER C EPSTAB(N) CONTAINS THE NEW ELEMENT IN THE C FIRST COLUMN OF THE EPSILON TABLE. C C EPSTAB - DOUBLE PRECISION C VECTOR OF DIMENSION 52 CONTAINING THE ELEMENTS C OF THE TWO LOWER DIAGONALS OF THE TRIANGULAR C EPSILON TABLE. THE ELEMENTS ARE NUMBERED C STARTING AT THE RIGHT-HAND CORNER OF THE C TRIANGLE. C C RESULT - DOUBLE PRECISION C RESULTING APPROXIMATION TO THE INTEGRAL C C ABSERR - DOUBLE PRECISION C ESTIMATE OF THE ABSOLUTE ERROR COMPUTED FROM C RESULT AND THE 3 PREVIOUS RESULTS C C RES3LA - DOUBLE PRECISION C VECTOR OF DIMENSION 3 CONTAINING THE LAST 3 C RESULTS C C NRES - INTEGER C NUMBER OF CALLS TO THE ROUTINE C (SHOULD BE ZERO AT FIRST CALL) C***END PROLOGUE DQELG C DOUBLE PRECISION ABSERR,DABS,DELTA1,DELTA2,DELTA3,DMAX1,D1MACH, 1 EPMACH,EPSINF,EPSTAB,ERROR,ERR1,ERR2,ERR3,E0,E1,E1ABS,E2,E3, 2 OFLOW,RES,RESULT,RES3LA,SS,TOL1,TOL2,TOL3 INTEGER I,IB,IB2,IE,INDX,K1,K2,K3,LIMEXP,N,NEWELM,NRES,NUM DIMENSION EPSTAB(52),RES3LA(3) C C LIST OF MAJOR VARIABLES C ----------------------- C C E0 - THE 4 ELEMENTS ON WHICH THE COMPUTATION OF A NEW C E1 ELEMENT IN THE EPSILON TABLE IS BASED C E2 C E3 E0 C E3 E1 NEW C E2 C NEWELM - NUMBER OF ELEMENTS TO BE COMPUTED IN THE NEW C DIAGONAL C ERROR - ERROR = ABS(E1-E0)+ABS(E2-E1)+ABS(NEW-E2) C RESULT - THE ELEMENT IN THE NEW DIAGONAL WITH LEAST VALUE C OF ERROR C C MACHINE DEPENDENT CONSTANTS C --------------------------- C C EPMACH IS THE LARGEST RELATIVE SPACING. C OFLOW IS THE LARGEST POSITIVE MAGNITUDE. C LIMEXP IS THE MAXIMUM NUMBER OF ELEMENTS THE EPSILON C TABLE CAN CONTAIN. IF THIS NUMBER IS REACHED, THE UPPER C DIAGONAL OF THE EPSILON TABLE IS DELETED. C C***FIRST EXECUTABLE STATEMENT DQELG EPMACH = D1MACH(4) OFLOW = D1MACH(2) NRES = NRES+1 ABSERR = OFLOW RESULT = EPSTAB(N) IF(N.LT.3) GO TO 100 LIMEXP = 50 EPSTAB(N+2) = EPSTAB(N) NEWELM = (N-1)/2 EPSTAB(N) = OFLOW 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 E1ABS = DABS(E1) DELTA2 = E2-E1 ERR2 = DABS(DELTA2) TOL2 = DMAX1(DABS(E2),E1ABS)*EPMACH DELTA3 = E1-E0 ERR3 = DABS(DELTA3) TOL3 = DMAX1(E1ABS,DABS(E0))*EPMACH 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 C ***JUMP OUT OF DO-LOOP GO TO 100 10 E3 = EPSTAB(K1) EPSTAB(K1) = E1 DELTA1 = E1-E3 ERR1 = DABS(DELTA1) TOL1 = DMAX1(E1ABS,DABS(E3))*EPMACH 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 EPSINF = DABS(SS*E1) 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(EPSINF.GT.0.1D-03) GO TO 30 20 N = I+I-1 C ***JUMP OUT OF DO-LOOP 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+DABS(RES-E2)+ERR3 IF(ERROR.GT.ABSERR) GO TO 40 ABSERR = ERROR RESULT = RES 40 CONTINUE 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 INDX = NUM-N+1 DO 70 I = 1,N EPSTAB(I)= EPSTAB(INDX) INDX = INDX+1 70 CONTINUE 80 IF(NRES.GE.4) GO TO 90 RES3LA(NRES) = RESULT ABSERR = OFLOW GO TO 100 C C COMPUTE ERROR ESTIMATE C 90 ABSERR = DABS(RESULT-RES3LA(3))+DABS(RESULT-RES3LA(2)) 1 +DABS(RESULT-RES3LA(1)) RES3LA(1) = RES3LA(2) RES3LA(2) = RES3LA(3) RES3LA(3) = RESULT 100 ABSERR = DMAX1(ABSERR,0.5D+01*EPMACH*DABS(RESULT)) RETURN END *DECK D1MACH DOUBLE PRECISION FUNCTION D1MACH(I) C***BEGIN PROLOGUE D1MACH C***DATE WRITTEN 750101 (YYMMDD) C***REVISION DATE 890213 (YYMMDD) C***CATEGORY NO. R1 C***KEYWORDS LIBRARY=SLATEC,TYPE=DOUBLE PRECISION(R1MACH-S D1MACH-D), C MACHINE CONSTANTS C***AUTHOR FOX, P. A., (BELL LABS) C HALL, A. D., (BELL LABS) C SCHRYER, N. L., (BELL LABS) C***PURPOSE RETURNS DOUBLE PRECISION MACHINE DEPENDENT CONSTANTS C***DESCRIPTION C C D1MACH CAN BE USED TO OBTAIN MACHINE-DEPENDENT PARAMETERS C FOR THE LOCAL MACHINE ENVIRONMENT. IT IS A FUNCTION C SUBPROGRAM WITH ONE (INPUT) ARGUMENT, AND CAN BE CALLED C AS FOLLOWS, FOR EXAMPLE C C D = D1MACH(I) C C WHERE I=1,...,5. THE (OUTPUT) VALUE OF D ABOVE IS C DETERMINED BY THE (INPUT) VALUE OF I. THE RESULTS FOR C VARIOUS VALUES OF I ARE DISCUSSED BELOW. C C D1MACH( 1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE. C D1MACH( 2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE. C D1MACH( 3) = B**(-T), THE SMALLEST RELATIVE SPACING. C D1MACH( 4) = B**(1-T), THE LARGEST RELATIVE SPACING. C D1MACH( 5) = LOG10(B) C C ASSUME DOUBLE PRECISION NUMBERS ARE REPRESENTED IN THE T-DIGIT, C BASE-B FORM C C SIGN (B**E)*( (X(1)/B) + ... + (X(T)/B**T) ) C C WHERE 0 .LE. X(I) .LT. B FOR I=1,...,T, 0 .LT. X(1), AND C EMIN .LE. E .LE. EMAX. C C THE VALUES OF B, T, EMIN AND EMAX ARE PROVIDED IN I1MACH AS C FOLLOWS: C I1MACH(10) = B, THE BASE. C I1MACH(14) = T, THE NUMBER OF BASE-B DIGITS. C I1MACH(15) = EMIN, THE SMALLEST EXPONENT E. C I1MACH(16) = EMAX, THE LARGEST EXPONENT E. C C TO ALTER THIS FUNCTION FOR A PARTICULAR ENVIRONMENT, C THE DESIRED SET OF DATA STATEMENTS SHOULD BE ACTIVATED BY C REMOVING THE C FROM COLUMN 1. ALSO, THE VALUES OF C D1MACH(1) - D1MACH(4) SHOULD BE CHECKED FOR CONSISTENCY C WITH THE LOCAL OPERATING SYSTEM. C C***REFERENCES FOX P.A., HALL A.D., SCHRYER N.L.,*FRAMEWORK FOR A C PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHEMATICAL C SOFTWARE, VOL. 4, NO. 2, JUNE 1978, PP. 177-188. C***ROUTINES CALLED XERROR C C THE CALL TO XERROR HAS BEEN COMMENTED TO GET A STAND ALONE ROUTINE C IT CAN BE ADAPTED TO FIT THE LOCAL ERROR HANDLING PROCEDURES C NOTE ADDED BY F. ROMANI 7/11/89 C C***END PROLOGUE D1MACH C INTEGER SMALL(4) INTEGER LARGE(4) INTEGER RIGHT(4) INTEGER DIVER(4) INTEGER LOG10(4) C DOUBLE PRECISION DMACH(5) SAVE DMACH C EQUIVALENCE (DMACH(1),SMALL(1)) EQUIVALENCE (DMACH(2),LARGE(1)) EQUIVALENCE (DMACH(3),RIGHT(1)) EQUIVALENCE (DMACH(4),DIVER(1)) EQUIVALENCE (DMACH(5),LOG10(1)) C C MACHINE CONSTANTS FOR THE AMIGA C ABSOFT FORTRAN COMPILER USING THE 68020/68881 COMPILER OPTION C C DATA SMALL(1), SMALL(2) / Z'00100000', Z'00000000' / C DATA LARGE(1), LARGE(2) / Z'7FEFFFFF', Z'FFFFFFFF' / C DATA RIGHT(1), RIGHT(2) / Z'3CA00000', Z'00000000' / C DATA DIVER(1), DIVER(2) / Z'3CB00000', Z'00000000' / C DATA LOG10(1), LOG10(2) / Z'3FD34413', Z'509F79FF' / C C MACHINE CONSTANTS FOR THE AMIGA C ABSOFT FORTRAN COMPILER USING SOFTWARE FLOATING POINT C C DATA SMALL(1), SMALL(2) / Z'00100000', Z'00000000' / C DATA LARGE(1), LARGE(2) / Z'7FDFFFFF', Z'FFFFFFFF' / C DATA RIGHT(1), RIGHT(2) / Z'3CA00000', Z'00000000' / C DATA DIVER(1), DIVER(2) / Z'3CB00000', Z'00000000' / C DATA LOG10(1), LOG10(2) / Z'3FD34413', Z'509F79FF' / C C MACHINE CONSTANTS FOR THE APOLLO C C DATA SMALL(1), SMALL(2) / 16#00100000, 16#00000000 / C DATA LARGE(1), LARGE(2) / 16#7FFFFFFF, 16#FFFFFFFF / C DATA RIGHT(1), RIGHT(2) / 16#3CA00000, 16#00000000 / C DATA DIVER(1), DIVER(2) / 16#3CB00000, 16#00000000 / C DATA LOG10(1), LOG10(2) / 16#3FD34413, 16#509F79FF / C C MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM C C DATA SMALL(1) / ZC00800000 / C DATA SMALL(2) / Z000000000 / C DATA LARGE(1) / ZDFFFFFFFF / C DATA LARGE(2) / ZFFFFFFFFF / C DATA RIGHT(1) / ZCC5800000 / C DATA RIGHT(2) / Z000000000 / C DATA DIVER(1) / ZCC6800000 / C DATA DIVER(2) / Z000000000 / C DATA LOG10(1) / ZD00E730E7 / C DATA LOG10(2) / ZC77800DC0 / C C MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM C C DATA SMALL(1) / O1771000000000000 / C DATA SMALL(2) / O0000000000000000 / C DATA LARGE(1) / O0777777777777777 / C DATA LARGE(2) / O0007777777777777 / C DATA RIGHT(1) / O1461000000000000 / C DATA RIGHT(2) / O0000000000000000 / C DATA DIVER(1) / O1451000000000000 / C DATA DIVER(2) / O0000000000000000 / C DATA LOG10(1) / O1157163034761674 / C DATA LOG10(2) / O0006677466732724 / C C MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS C C DATA SMALL(1) / O1771000000000000 / C DATA SMALL(2) / O7770000000000000 / C DATA LARGE(1) / O0777777777777777 / C DATA LARGE(2) / O7777777777777777 / C DATA RIGHT(1) / O1461000000000000 / C DATA RIGHT(2) / O0000000000000000 / C DATA DIVER(1) / O1451000000000000 / C DATA DIVER(2) / O0000000000000000 / C DATA LOG10(1) / O1157163034761674 / C DATA LOG10(2) / O0006677466732724 / C C MACHINE CONSTANTS FOR THE CDC 170/180 SERIES USING NOS/VE C C DATA SMALL(1) / Z"3001800000000000" / C DATA SMALL(2) / Z"3001000000000000" / C DATA LARGE(1) / Z"4FFEFFFFFFFFFFFE" / C DATA LARGE(2) / Z"4FFE000000000000" / C DATA RIGHT(1) / Z"3FD2800000000000" / C DATA RIGHT(2) / Z"3FD2000000000000" / C DATA DIVER(1) / Z"3FD3800000000000" / C DATA DIVER(2) / Z"3FD3000000000000" / C DATA LOG10(1) / Z"3FFF9A209A84FBCF" / C DATA LOG10(2) / Z"3FFFF7988F8959AC" / C C MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES C C DATA SMALL(1) / 00564000000000000000B / C DATA SMALL(2) / 00000000000000000000B / C DATA LARGE(1) / 37757777777777777777B / C DATA LARGE(2) / 37157777777777777777B / C DATA RIGHT(1) / 15624000000000000000B / C DATA RIGHT(2) / 00000000000000000000B / C DATA DIVER(1) / 15634000000000000000B / C DATA DIVER(2) / 00000000000000000000B / C DATA LOG10(1) / 17164642023241175717B / C DATA LOG10(2) / 16367571421742254654B / C C MACHINE CONSTANTS FOR THE CELERITY C1260 C C DATA SMALL(1), SMALL(2) / Z'00100000', Z'00000000' / C DATA LARGE(1), LARGE(2) / Z'7FEFFFFF', Z'FFFFFFFF' / C DATA RIGHT(1), RIGHT(2) / Z'3CA00000', Z'00000000' / C DATA DIVER(1), DIVER(2) / Z'3CB00000', Z'00000000' / C DATA LOG10(1), LOG10(2) / Z'3FD34413', Z'509F79FF' / C C MACHINE CONSTANTS FOR THE CONVEX C-1 C C DATA SMALL(1), SMALL(2) / '00100000'X,'00000000'X / C DATA LARGE(1), LARGE(2) / '7FFFFFFF'X,'FFFFFFFF'X / C DATA RIGHT(1), RIGHT(2) / '3CC00000'X,'00000000'X / C DATA DIVER(1), DIVER(2) / '3CD00000'X,'00000000'X / C DATA LOG10(1), LOG10(2) / '3FF34413'X,'509F79FF'X / C C MACHINE CONSTANTS FOR THE CRAY-1 C C DATA SMALL(1) / 201354000000000000000B / C DATA SMALL(2) / 000000000000000000000B / C DATA LARGE(1) / 577767777777777777777B / C DATA LARGE(2) / 000007777777777777774B / C DATA RIGHT(1) / 376434000000000000000B / C DATA RIGHT(2) / 000000000000000000000B / C DATA DIVER(1) / 376444000000000000000B / C DATA DIVER(2) / 000000000000000000000B / C DATA LOG10(1) / 377774642023241175717B / C DATA LOG10(2) / 000007571421742254654B / C C MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200 C C NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING CARD - C STATIC DMACH(5) C C DATA SMALL / 20K, 3*0 / C DATA LARGE / 77777K, 3*177777K / C DATA RIGHT / 31420K, 3*0 / C DATA DIVER / 32020K, 3*0 / C DATA LOG10 / 40423K, 42023K, 50237K, 74776K / C C MACHINE CONSTANTS FOR THE ELXSI 6400 C (ASSUMING REAL*8 IS THE DEFAULT DOUBLE PRECISION) C C DATA SMALL(1), SMALL(2) / '00100000'X,'00000000'X / C DATA LARGE(1), LARGE(2) / '7FEFFFFF'X,'FFFFFFFF'X / C DATA RIGHT(1), RIGHT(2) / '3CB00000'X,'00000000'X / C DATA DIVER(1), DIVER(2) / '3CC00000'X,'00000000'X / C DATA LOG10(1), LOG10(2) / '3FD34413'X,'509F79FF'X / C C MACHINE CONSTANTS FOR THE HARRIS 220 C C DATA SMALL(1), SMALL(2) / '20000000, '00000201 / C DATA LARGE(1), LARGE(2) / '37777777, '37777577 / C DATA RIGHT(1), RIGHT(2) / '20000000, '00000333 / C DATA DIVER(1), DIVER(2) / '20000000, '00000334 / C DATA LOG10(1), LOG10(2) / '23210115, '10237777 / C C MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES C C DATA SMALL(1), SMALL(2) / O402400000000, O000000000000 / C DATA LARGE(1), LARGE(2) / O376777777777, O777777777777 / C DATA RIGHT(1), RIGHT(2) / O604400000000, O000000000000 / C DATA DIVER(1), DIVER(2) / O606400000000, O000000000000 / C DATA LOG10(1), LOG10(2) / O776464202324, O117571775714 / C C MACHINE CONSTANTS FOR THE HP 2100 C THREE WORD DOUBLE PRECISION OPTION WITH FTN4 C C DATA SMALL(1), SMALL(2), SMALL(3) / 40000B, 0, 1 / C DATA LARGE(1), LARGE(2), LARGE(3) / 77777B, 177777B, 177776B / C DATA RIGHT(1), RIGHT(2), RIGHT(3) / 40000B, 0, 265B / C DATA DIVER(1), DIVER(2), DIVER(3) / 40000B, 0, 276B / C DATA LOG10(1), LOG10(2), LOG10(3) / 46420B, 46502B, 77777B / C C MACHINE CONSTANTS FOR THE HP 2100 C FOUR WORD DOUBLE PRECISION OPTION WITH FTN4 C C DATA SMALL(1), SMALL(2) / 40000B, 0 / C DATA SMALL(3), SMALL(4) / 0, 1 / C DATA LARGE(1), LARGE(2) / 77777B, 177777B / C DATA LARGE(3), LARGE(4) / 177777B, 177776B / C DATA RIGHT(1), RIGHT(2) / 40000B, 0 / C DATA RIGHT(3), RIGHT(4) / 0, 225B / C DATA DIVER(1), DIVER(2) / 40000B, 0 / C DATA DIVER(3), DIVER(4) / 0, 227B / C DATA LOG10(1), LOG10(2) / 46420B, 46502B / C DATA LOG10(3), LOG10(4) / 76747B, 176377B / C C MACHINE CONSTANTS FOR THE HP 9000 C C DATA SMALL(1), SMALL(2) / 00040000000B, 00000000000B / C DATA LARGE(1), LARGE(2) / 17737777777B, 37777777777B / C DATA RIGHT(1), RIGHT(2) / 07454000000B, 00000000000B / C DATA DIVER(1), DIVER(2) / 07460000000B, 00000000000B / C DATA LOG10(1), LOG10(2) / 07764642023B, 12047674777B / C C MACHINE CONSTANTS FOR THE IBM 360/370 SERIES, C THE XEROX SIGMA 5/7/9, THE SEL SYSTEMS 85/86, AND C THE PERKIN ELMER (INTERDATA) 7/32. C C DATA SMALL(1), SMALL(2) / Z00100000, Z00000000 / C DATA LARGE(1), LARGE(2) / Z7FFFFFFF, ZFFFFFFFF / C DATA RIGHT(1), RIGHT(2) / Z33100000, Z00000000 / C DATA DIVER(1), DIVER(2) / Z34100000, Z00000000 / C DATA LOG10(1), LOG10(2) / Z41134413, Z509F79FF / C C MACHINE CONSTANTS FOR THE IBM PC C ASSUMES THAT ALL ARITHMETIC IS DONE IN DOUBLE PRECISION C ON 8088, I.E., NOT IN 80 BIT FORM FOR THE 8087. C C DATA SMALL(1) / 2.23D-308 / C DATA LARGE(1) / 1.79D+308 / C DATA RIGHT(1) / 1.11D-16 / C DATA DIVER(1) / 2.22D-16 / C DATA LOG10(1) / 0.301029995663981195D0 / C C MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR) C C DATA SMALL(1), SMALL(2) / "033400000000, "000000000000 / C DATA LARGE(1), LARGE(2) / "377777777777, "344777777777 / C DATA RIGHT(1), RIGHT(2) / "113400000000, "000000000000 / C DATA DIVER(1), DIVER(2) / "114400000000, "000000000000 / C DATA LOG10(1), LOG10(2) / "177464202324, "144117571776 / C C MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR) C C DATA SMALL(1), SMALL(2) / "000400000000, "000000000000 / C DATA LARGE(1), LARGE(2) / "377777777777, "377777777777 / C DATA RIGHT(1), RIGHT(2) / "103400000000, "000000000000 / C DATA DIVER(1), DIVER(2) / "104400000000, "000000000000 / C DATA LOG10(1), LOG10(2) / "177464202324, "476747767461 / C C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING C 32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C C DATA SMALL(1), SMALL(2) / 8388608, 0 / C DATA LARGE(1), LARGE(2) / 2147483647, -1 / C DATA RIGHT(1), RIGHT(2) / 612368384, 0 / C DATA DIVER(1), DIVER(2) / 620756992, 0 / C DATA LOG10(1), LOG10(2) / 1067065498, -2063872008 / C C DATA SMALL(1), SMALL(2) / O00040000000, O00000000000 / C DATA LARGE(1), LARGE(2) / O17777777777, O37777777777 / C DATA RIGHT(1), RIGHT(2) / O04440000000, O00000000000 / C DATA DIVER(1), DIVER(2) / O04500000000, O00000000000 / C DATA LOG10(1), LOG10(2) / O07746420232, O20476747770 / C C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING C 16-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C C DATA SMALL(1), SMALL(2) / 128, 0 / C DATA SMALL(3), SMALL(4) / 0, 0 / C DATA LARGE(1), LARGE(2) / 32767, -1 / C DATA LARGE(3), LARGE(4) / -1, -1 / C DATA RIGHT(1), RIGHT(2) / 9344, 0 / C DATA RIGHT(3), RIGHT(4) / 0, 0 / C DATA DIVER(1), DIVER(2) / 9472, 0 / C DATA DIVER(3), DIVER(4) / 0, 0 / C DATA LOG10(1), LOG10(2) / 16282, 8346 / C DATA LOG10(3), LOG10(4) / -31493, -12296 / C C DATA SMALL(1), SMALL(2) / O000200, O000000 / C DATA SMALL(3), SMALL(4) / O000000, O000000 / C DATA LARGE(1), LARGE(2) / O077777, O177777 / C DATA LARGE(3), LARGE(4) / O177777, O177777 / C DATA RIGHT(1), RIGHT(2) / O022200, O000000 / C DATA RIGHT(3), RIGHT(4) / O000000, O000000 / C DATA DIVER(1), DIVER(2) / O022400, O000000 / C DATA DIVER(3), DIVER(4) / O000000, O000000 / C DATA LOG10(1), LOG10(2) / O037632, O020232 / C DATA LOG10(3), LOG10(4) / O102373, O147770 / C C MACHINE CONSTANTS FOR THE SUN C C DATA SMALL(1), SMALL(2) / Z'00100000', Z'00000000' / C DATA LARGE(1), LARGE(2) / Z'7FEFFFFF', Z'FFFFFFFF' / C DATA RIGHT(1), RIGHT(2) / Z'3CA00000', Z'00000000' / C DATA DIVER(1), DIVER(2) / Z'3CB00000', Z'00000000' / C DATA LOG10(1), LOG10(2) / Z'3FD34413', Z'509F79FF' / C C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES FTN COMPILER C C DATA SMALL(1), SMALL(2) / O000040000000, O000000000000 / C DATA LARGE(1), LARGE(2) / O377777777777, O777777777777 / C DATA RIGHT(1), RIGHT(2) / O170540000000, O000000000000 / C DATA DIVER(1), DIVER(2) / O170640000000, O000000000000 / C DATA LOG10(1), LOG10(2) / O177746420232, O411757177572 / C C MACHINE CONSTANTS FOR VAX 11/780 C (EXPRESSED IN INTEGER AND HEXADECIMAL) C THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSYEMS C THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS C C DATA SMALL(1), SMALL(2) / 128, 0 / C DATA LARGE(1), LARGE(2) / -32769, -1 / C DATA RIGHT(1), RIGHT(2) / 9344, 0 / C DATA DIVER(1), DIVER(2) / 9472, 0 / C DATA LOG10(1), LOG10(2) / 546979738, -805796613 / C C DATA SMALL(1), SMALL(2) / Z00000080, Z00000000 / C DATA LARGE(1), LARGE(2) / ZFFFF7FFF, ZFFFFFFFF / C DATA RIGHT(1), RIGHT(2) / Z00002480, Z00000000 / C DATA DIVER(1), DIVER(2) / Z00002500, Z00000000 / C DATA LOG10(1), LOG10(2) / Z209A3F9A, ZCFF884FB / C C MACHINE CONSTANTS FOR VAX 11/780 (G-FLOATING) C (EXPRESSED IN INTEGER AND HEXADECIMAL) C THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSYEMS C THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS C C DATA SMALL(1), SMALL(2) / 16, 0 / C DATA LARGE(1), LARGE(2) / -32769, -1 / C DATA RIGHT(1), RIGHT(2) / 15552, 0 / C DATA DIVER(1), DIVER(2) / 15568, 0 / C DATA LOG10(1), LOG10(2) / 1142112243, 2046775455 / C C DATA SMALL(1), SMALL(2) / Z00000010, Z00000000 / C DATA LARGE(1), LARGE(2) / ZFFFF7FFF, ZFFFFFFFF / C DATA RIGHT(1), RIGHT(2) / Z00003CC0, Z00000000 / C DATA DIVER(1), DIVER(2) / Z00003CD0, Z00000000 / C DATA LOG10(1), LOG10(2) / Z44133FF3, Z79FF509F / C C C***FIRST EXECUTABLE STATEMENT D1MACH C C THE CALL TO XERROR HAS BEEN COMMENTED TO GET A STAND ALONE ROUTINE C IT CAN BE ADAPTED TO FIT THE LOCAL ERROR HANDLING PROCEDURES C NOTE ADDED BY F. ROMANI 7/11/89 C C IF (I .LT. 1 .OR. I .GT. 5) C 1 CALL XERROR ('D1MACH -- I OUT OF BOUNDS', 25, 1, 2) C D1MACH = DMACH(I) RETURN C END *DECK DQXGS SUBROUTINE DQXGS(F,A,B,EPSABS,EPSREL,RESULT,ABSERR,IER, 1 LIMIT,LENIW,LENW,LAST,IWORK,WORK) C C THE ROUTINE CALCULATES AN APPROXIMATION RESULT TO A GIVEN C DEFINITE INTEGRAL I = INTEGRAL OF F OVER (A,B), C HOPEFULLY SATISFYING FOLLOWING CLAIM FOR ACCURACY C ABS(I-RESULT).LE.MAX(EPSABS,EPSREL*ABS(I)). 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 A - DOUBLE PRECISION C LOWER LIMIT OF INTEGRATION C C B - DOUBLE PRECISION C UPPER LIMIT OF INTEGRATION 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 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 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 C THE ESTIMATES FOR INTEGRAL AND ERROR ARE C LESS RELIABLE. IT IS ASSUMED THAT THE C REQUESTED ACCURACY HAS NOT BEEN ACHIEVED. C ERROR MESSAGES C IER = 1 MAXIMUM NUMBER OF SUBDIVISIONS ALLOWED C HAS BEEN ACHIEVED. ONE CAN ALLOW MORE SUB- C DIVISIONS BY INCREASING THE VALUE OF LIMIT C (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 DETEC- C TED, WHICH PREVENTS THE REQUESTED C TOLERANCE FROM BEING ACHIEVED. C THE ERROR MAY BE UNDER-ESTIMATED. C = 3 EXTREMELY BAD INTEGRAND BEHAVIOUR C OCCURS 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. IT IS PRESUMED THAT C THE REQUESTED TOLERANCE CANNOT BE C ACHIEVED, AND THAT THE RETURNED RESULT IS C 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 LENW.LT.LIMIT*46 OR C LENIW .LT. LIMIT*3 C RESULT, ABSERR, LAST ARE SET TO C ZERO.EXCEPT WHEN LIMIT OR LENW OR LENIW C IS 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 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*46 C IF LENW.LT.LIMIT*46, THE ROUTINE WILL END C WITH IER = 6. C C LENIW - INTEGER C DIMENSIONING PARAMETER FOR IWORK C LENIW MUST BE AT LEAST LIMIT*3 C IF LENW.LT.LIMIT*3, 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, DETERMINES THE C NUMBER OF SIGNIFICANT ELEMENTS ACTUALLY IN THE WORK C ARRAYS. C C WORK ARRAYS C IWORK - INTEGER C VECTOR OF DIMENSION AT LEAST 3*LIMIT, THE FIRST K C 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), C AND 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 C THE INTEGRAL APPROXIMATIONS OVER THE SUBINTERVALS, C WORK(LIMIT*3+1), ..., WORK(LIMIT*3+LAST) C CONTAIN THE ERROR ESTIMATES. C WORK(LIMIT*4+1), ... IS THE AREA RESERVED TO STORE C FUNCTIONAL VALUES. C***ROUTINES CALLED DQXGSE,XERROR C C THE CALL TO XERROR HAS BEEN COMMENTED TO GET A STAND ALONE ROUTINE C IT CAN BE ADAPTED TO FIT THE LOCAL ERROR HANDLING PROCEDURES C DOUBLE PRECISION A,ABSERR,B,EPSABS,EPSREL,F,RESULT,WORK INTEGER IER,IWORK,LAST,LENIW,LENW,LIMIT,LVL,L1,L2,L3,L4,L5 C DIMENSION IWORK(LENIW),WORK(LENW) C EXTERNAL F C C CHECK VALIDITY OF LIMIT,LENIW AND LENW. C C***FIRST EXECUTABLE STATEMENT DQXGS IER = 6 LAST = 0 RESULT = 0.0D+00 ABSERR = 0.0D+00 IF(LIMIT.LT.1.OR.LENIW.LT.LIMIT*3.OR.LENW.LT.LIMIT*46) GO TO 10 C C PREPARE CALL FOR DQXGSE. C L1 = LIMIT+1 L2 = LIMIT+L1 L3 = LIMIT+L2 L4 = LIMIT+L3 L5 = 21*LIMIT+L3 C CALL DQXGSE(F,A,B,EPSABS,EPSREL,LIMIT,RESULT,ABSERR, * IER,WORK(1),WORK(L1),WORK(L2),WORK(L3),IWORK,LAST, * WORK(L4),WORK(L5),IWORK(L1),IWORK(L2)) C C CALL ERROR HANDLER IF NECESSARY. C LVL = 0 10 IF(IER.EQ.6) LVL = 1 C C THE CALL TO XERROR HAS BEEN COMMENTED TO GET A STAND ALONE ROUTINE C IT CAN BE ADAPTED TO FIT THE LOCAL ERROR HANDLING PROCEDURES C C IF(IER.NE.0)CALL XERROR('ABNORMAL RETURN FROM DQXGS ',28,IER,LVL) RETURN END *DECK DQXGSE SUBROUTINE DQXGSE(F,A,B,EPSABS,EPSREL,LIMIT,RESULT,ABSERR, * IER,ALIST,BLIST,RLIST,ELIST,IORD,LAST,VALP,VALN,LP,LN) C C THE ROUTINE CALCULATES AN APPROXIMATION RESULT TO A GIVEN C DEFINITE INTEGRAL I = INTEGRAL OF F OVER (A,B), C HOPEFULLY SATISFYING FOLLOWING CLAIM FOR ACCURACY C ABS(I-RESULT).LE.MAX(EPSABS,EPSREL*ABS(I)). 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 A - DOUBLE PRECISION C LOWER LIMIT OF INTEGRATION C C B - DOUBLE PRECISION C UPPER LIMIT OF INTEGRATION 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 UPPERBOUND ON THE NUMBER OF SUBINTERVALS C IN THE PARTITION OF (A,B) 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 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 C THE ESTIMATES FOR INTEGRAL AND ERROR ARE C LESS RELIABLE. IT IS ASSUMED THAT THE C REQUESTED ACCURACY HAS NOT BEEN ACHIEVED. C ERROR MESSAGES C = 1 MAXIMUM NUMBER OF SUBDIVISIONS ALLOWED C HAS BEEN ACHIEVED. ONE CAN ALLOW MORE SUB- C DIVISIONS BY INCREASING THE VALUE OF LIMIT C (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 DETEC- C TED, WHICH PREVENTS THE REQUESTED C TOLERANCE FROM BEING ACHIEVED. C THE ERROR MAY BE UNDER-ESTIMATED. C = 3 EXTREMELY BAD INTEGRAND BEHAVIOUR C OCCURS 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 PRESUMED THAT THE REQUESTED C TOLERANCE CANNOT BE ACHIEVED, AND THAT THE C RETURNED RESULT IS THE BEST WHICH CAN BE C 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, LAST, RLIST(1), C IORD(1) AND ELIST(1) ARE SET TO ZERO. C ALIST(1) AND BLIST(1) ARE SET TO A AND B C RESPECTIVELY. C C ALIST - DOUBLE PRECISION C VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST C LAST ELEMENTS OF WHICH ARE THE LEFT END POINTS C OF THE SUBINTERVALS IN THE PARTITION OF THE C GIVEN INTEGRATION RANGE (A,B) C C BLIST - DOUBLE PRECISION C VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST C LAST ELEMENTS OF WHICH ARE THE RIGHT END POINTS C OF THE SUBINTERVALS IN THE PARTITION OF THE GIVEN C INTEGRATION RANGE (A,B) 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 AT LEAST 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 IN THE C SUBDIVISION PROCESS C C VALP - DOUBLE PRECISION C VALN ARRAYS OF DIMENSION AT LEAST (21,LIMIT) USED TO C SAVE THE FUNCTIONAL VALUES C C LP - INTEGER C LN VECTORS OF DIMENSION AT LEAST LIMIT, USED TO C STORE THE ACTUAL NUMBER OF FUNCTIONAL VALUES C SAVED IN THE CORRESPONDING COLUMN C OF VALP,VALN C C***ROUTINES CALLED F,D1MACH,DQELG,DQXLQM,DQPSRT,DQXRRD,DQXCPY C DOUBLE PRECISION A,ABSEPS,ABSERR,ALIST,AREA,AREA1,AREA12,AREA2,A1, * A2,B,BLIST,B1,B2,CORREC,DABS,DEFABS,DEFAB1,DEFAB2,D1MACH,DMAX1, * DRES,ELIST,EPMACH,EPSABS,EPSREL,ERLARG,ERLAST,ERRBND,ERRMAX, * ERROR1,ERROR2,ERRO12,ERRSUM,ERTEST,F,OFLOW,RESABS,RESEPS,RESULT, * RES3LA,RLIST,RLIST2,SMALL,UFLOW, * VALP,VALN,VP1,VP2,VN1,VN2 INTEGER ID,IER,IERRO,IORD,IROFF1,IROFF2,IROFF3,JUPBND,K,KSGN, * KTMIN,LAST,LIMIT,MAXERR,NRES,NRMAX,NUMRL2, * LP,LN,LP1,LP2,LN1,LN2 LOGICAL EXTRAP,NOEXT DIMENSION ALIST(LIMIT),BLIST(LIMIT),ELIST(LIMIT),IORD(LIMIT), * RES3LA(3),RLIST(LIMIT),RLIST2(52), * VALP(21,LIMIT),VALN(21,LIMIT),LP(LIMIT),LN(LIMIT), * VP1(21),VP2(21),VN1(21),VN2(21) EXTERNAL F C C MACHINE DEPENDENT CONSTANTS C --------------------------- C C EPMACH IS THE LARGEST RELATIVE SPACING. C UFLOW IS THE SMALLEST POSITIVE MAGNITUDE. C OFLOW IS THE LARGEST POSITIVE MAGNITUDE. C C***FIRST EXECUTABLE STATEMENT DQXGSE EPMACH = D1MACH(4) C C TEST ON VALIDITY OF PARAMETERS C ------------------------------ IER = 0 LAST = 0 RESULT = 0.0D+00 ABSERR = 0.0D+00 ALIST(1) = A BLIST(1) = B RLIST(1) = 0.0D+00 ELIST(1) = 0.0D+00 IF(EPSABS.LE.0.0D+00.AND.EPSREL.LT.DMAX1(0.5D+02*EPMACH,0.5D-28)) * IER = 6 IF(IER.EQ.6) GO TO 999 C C FIRST APPROXIMATION TO THE INTEGRAL C ----------------------------------- C UFLOW = D1MACH(1) OFLOW = D1MACH(2) IERRO = 0 LP(1)=1 LN(1)=1 VALP(1,1)=F((A+B)*0.5D0) VALN(1,1)=VALP(1,1) CALL DQXLQM(F,A,B,RESULT,ABSERR,DEFABS,RESABS, * VALP(1,1),VALN(1,1),LP(1),LN(1), 2 ) C C TEST ON ACCURACY. C DRES = DABS(RESULT) ERRBND = DMAX1(EPSABS,EPSREL*DRES) LAST = 1 RLIST(1) = RESULT ELIST(1) = ABSERR IORD(1) = 1 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. * ABSERR.EQ.0.0D+00) GO TO 999 C C INITIALIZATION C -------------- C RLIST2(1) = RESULT ERRMAX = ABSERR MAXERR = 1 AREA = RESULT ERRSUM = ABSERR ABSERR = OFLOW NRMAX = 1 NRES = 0 NUMRL2 = 2 KTMIN = 0 EXTRAP = .FALSE. NOEXT = .FALSE. 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 THE NRMAX-TH LARGEST ERROR C ESTIMATE. C A1 = ALIST(MAXERR) B1 = 0.5D+00*(ALIST(MAXERR)+BLIST(MAXERR)) A2 = B1 B2 = BLIST(MAXERR) ERLAST = ERRMAX CALL DQXRRD(F,VALN(1,MAXERR),LN(MAXERR),B1,A1,VN1,VP1,LN1,LP1) CALL DQXLQM(F,A1,B1,AREA1,ERROR1,RESABS,DEFAB1,VP1,VN1,LP1,LN1, * 2) CALL DQXRRD(F,VALP(1,MAXERR),LP(MAXERR),A2,B2,VP2,VN2,LP2,LN2) CALL DQXLQM(F,A2,B2,AREA2,ERROR2,RESABS,DEFAB2,VP2,VN2,LP2,LN2, * 2) 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) * .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 SUBINTERVALS C EQUALS LIMIT. C IF(LAST.EQ.LIMIT) IER = 1 C C SET ERROR FLAG IN THE CASE OF BAD INTEGRAND BEHAVIOUR C AT A POINT OF THE INTEGRATION RANGE. C IF(DMAX1(DABS(A1),DABS(B2)).LE.(0.1D+01+0.1D+03*EPMACH)* * (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 CALL DQXCPY(VALP(1,MAXERR),VP1,LP1) LP(MAXERR)=LP1 CALL DQXCPY(VALN(1,MAXERR),VN1,LN1) LN(MAXERR)=LN1 CALL DQXCPY(VALP(1,LAST),VP2,LP2) LP(LAST)=LP2 CALL DQXCPY(VALN(1,LAST),VN2,LN2) LN(LAST)=LN2 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 CALL DQXCPY(VALP(1,MAXERR),VP2,LP2) LP(MAXERR)=LP2 CALL DQXCPY(VALN(1,MAXERR),VN2,LN2) LN(MAXERR)=LN2 CALL DQXCPY(VALP(1,LAST),VP1,LP1) LP(LAST)=LP1 CALL DQXCPY(VALN(1,LAST),VN1,LN1) LN(LAST)=LN1 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) C ***JUMP OUT OF DO-LOOP IF(ERRSUM.LE.ERRBND) GO TO 115 C ***JUMP OUT OF DO-LOOP 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 C C THE BOUND 0.3*ERTEST HAS BEEN INTRODUCED TO PERFORM A C MORE CAUTIOUS EXTRAPOLATION THAN IN THE ORIGINAL DQAGSE C ROUTINE C 40 IF(IERRO.EQ.3.OR.ERLARG.LE.0.3D0*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) C ***JUMP OUT OF DO-LOOP IF(DABS(BLIST(MAXERR)-ALIST(MAXERR)).GT.SMALL) GO TO 90 NRMAX = NRMAX+1 50 CONTINUE C C PERFORM EXTRAPOLATION. C 60 NUMRL2 = NUMRL2+1 RLIST2(NUMRL2) = AREA CALL DQELG(NUMRL2,RLIST2,RESEPS,ABSEPS,RES3LA,NRES) KTMIN = KTMIN+1 IF(KTMIN.GT.5.AND.ABSERR.LT.0.1D-02*ERRSUM) IER = 5 IF(ABSEPS.GE.ABSERR) GO TO 70 KTMIN = 0 ABSERR = ABSEPS RESULT = RESEPS CORREC = ERLARG ERTEST = DMAX1(EPSABS,EPSREL*DABS(RESEPS)) C ***JUMP OUT OF DO-LOOP IF(ABSERR.LE.ERTEST) GO TO 100 C C PREPARE BISECTION OF THE SMALLEST INTERVAL. C 70 IF(NUMRL2.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 = DABS(B-A)*0.375D+00 ERLARG = ERRSUM ERTEST = ERRBND RLIST2(2) = AREA 90 CONTINUE C C SET FINAL RESULT AND ERROR ESTIMATE. C ------------------------------------ C 100 IF(ABSERR.EQ.OFLOW) 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. * DEFABS*0.1D-01) GO TO 130 IF(0.1D-01.GT.(RESULT/AREA).OR.(RESULT/AREA).GT.0.1D+03 * .OR.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 IF(IER.GT.2) IER = IER-1 999 RETURN END *DECK DQXG SUBROUTINE DQXG(F,A,B,EPSABS,EPSREL,KEY,RESULT,ABSERR,IER, 1 LIMIT,LENIW,LENW,LAST,IWORK,WORK) C C THE ROUTINE CALCULATES AN APPROXIMATION RESULT TO A GIVEN C DEFINITE INTEGRAL I = INTEGRAL OF F OVER (A,B), C HOPEFULLY SATISFYING FOLLOWING CLAIM FOR ACCURACY C ABS(I-RESULT).LE.MAX(EPSABS,EPSREL*ABS(I)). 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 A - DOUBLE PRECISION C LOWER LIMIT OF INTEGRATION C C B - DOUBLE PRECISION C UPPER LIMIT OF INTEGRATION 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 KEY - INTEGER C KEY FOR CHOICE OF LOCAL INTEGRATION RULE C RMS FORMULAS ARE USED WITH C 13 - 19 POINTS IF KEY.LT.1, C 13 - 19 - (27) POINTS IF KEY = 1, C 13 - 19 - (27) - (41) POINTS IF KEY = 2, C 19 - 27 - (41) POINTS IF KEY = 3, C 27 - 41 POINTS IF KEY.GT.3. C C (RULES) USED IF THE FUNCTION APPEARS C ENOUGH REGULAR IN THE SUBINTERVAL 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 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 C THE ESTIMATES FOR RESULT AND ERROR ARE C LESS RELIABLE. IT IS ASSUMED THAT THE C REQUESTED ACCURACY HAS NOT BEEN ACHIEVED. C ERROR MESSAGES C IER = 1 MAXIMUM NUMBER OF SUBDIVISIONS ALLOWED C HAS BEEN ACHIEVED. ONE CAN ALLOW MORE SUB- C DIVISIONS BY INCREASING THE VALUE OF LIMIT C (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 DETEC- C TED, WHICH PREVENTS THE REQUESTED C TOLERANCE FROM BEING ACHIEVED. C THE ERROR MAY BE UNDER-ESTIMATED. C = 3 EXTREMELY BAD INTEGRAND BEHAVIOUR C OCCURS 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. IT IS PRESUMED THAT C THE REQUESTED TOLERANCE CANNOT BE C ACHIEVED, AND THAT THE RETURNED RESULT IS C 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 LENW.LT.LIMIT*46 OR C LENIW .LT. LIMIT*3 C RESULT, ABSERR, LAST ARE SET TO C ZERO.EXCEPT WHEN LIMIT OR LENW OR LENIW C IS 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 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*46 C IF LENW.LT.LIMIT*46, THE ROUTINE WILL END C WITH IER = 6. C C LENIW - INTEGER C DIMENSIONING PARAMETER FOR IWORK C LENIW MUST BE AT LEAST LIMIT*3 C IF LENW.LT.LIMIT*3, 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, DETERMINES THE C NUMBER OF SIGNIFICANT ELEMENTS ACTUALLY IN THE WORK C ARRAYS. C C WORK ARRAYS C IWORK - INTEGER C VECTOR OF DIMENSION AT LEAST 3*LIMIT, THE FIRST K C 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), C AND 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 C THE INTEGRAL APPROXIMATIONS OVER THE SUBINTERVALS, C WORK(LIMIT*3+1), ..., WORK(LIMIT*3+LAST) C CONTAIN THE ERROR ESTIMATES. C WORK(LIMIT*4+1), ... IS THE AREA RESERVED TO STORE C FUNCTIONAL VALUES. C***ROUTINES CALLED DQXGE,XERROR C C THE CALL TO XERROR HAS BEEN COMMENTED TO GET A STAND ALONE ROUTINE C IT CAN BE ADAPTED TO FIT THE LOCAL ERROR HANDLING PROCEDURES C C DOUBLE PRECISION A,ABSERR,B,EPSABS,EPSREL,F,RESULT,WORK INTEGER KEY,IER,IWORK,LAST,LENIW,LENW,LIMIT,LVL,L1,L2,L3,L4,L5 C DIMENSION IWORK(LENIW),WORK(LENW) C EXTERNAL F C C CHECK VALIDITY OF LIMIT,LENIW AND LENW. C C***FIRST EXECUTABLE STATEMENT DQXG IER = 6 LAST = 0 RESULT = 0.0D+00 ABSERR = 0.0D+00 IF(LIMIT.LT.1.OR.LENIW.LT.LIMIT*3.OR.LENW.LT.LIMIT*46) GO TO 10 C C PREPARE CALL FOR DQXGE. C L1 = LIMIT+1 L2 = LIMIT+L1 L3 = LIMIT+L2 L4 = LIMIT+L3 L5 = 21*LIMIT+L3 C CALL DQXGE(F,A,B,EPSABS,EPSREL,KEY,LIMIT,RESULT,ABSERR, * IER,WORK(1),WORK(L1),WORK(L2),WORK(L3),IWORK,LAST, * WORK(L4),WORK(L5),IWORK(L1),IWORK(L2)) C C CALL ERROR HANDLER IF NECESSARY. C LVL = 0 10 IF(IER.EQ.6) LVL = 1 C C THE CALL TO XERROR HAS BEEN COMMENTED TO GET A STAND ALONE ROUTINE C IT CAN BE ADAPTED TO FIT THE LOCAL ERROR HANDLING PROCEDURES C C IF(IER.NE.0)CALL XERROR('ABNORMAL RETURN FROM DQXG ',28,IER,LVL) RETURN END *DECK DQXGE SUBROUTINE DQXGE(F,A,B,EPSABS,EPSREL,KEY,LIMIT,RESULT,ABSERR, * IER,ALIST,BLIST,RLIST,ELIST,IORD,LAST,VALP,VALN,LP,LN) C C C THE ROUTINE CALCULATES AN APPROXIMATION RESULT TO A GIVEN C DEFINITE INTEGRAL I = INTEGRAL OF F OVER (A,B), C HOPEFULLY SATISFYING FOLLOWING CLAIM FOR ACCURACY C ABS(I-RESLT).LE.MAX(EPSABS,EPSREL*ABS(I)). 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 A - DOUBLE PRECISION C LOWER LIMIT OF INTEGRATION C C B - DOUBLE PRECISION C UPPER LIMIT OF INTEGRATION C C EPSABS - DOUBLE PRECISION C ABSOLUTE ACCURACY REQUESTED C 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 KEY - INTEGER C KEY FOR CHOICE OF LOCAL INTEGRATION RULE C RMS FORMULAS ARE USED WITH C 13 - 19 POINTS IF KEY.LT.1, C 13 - 19 - (27) POINTS IF KEY = 1, C 13 - 19 - (27) - (41) POINTS IF KEY = 2, C 19 - 27 - (41) POINTS IF KEY = 3, C 27 - 41 POINTS IF KEY.GT.3. C C (RULES) USED IF THE FUNCTION APPEARS C ENOUGH REGULAR IN THE SUBINTERVAL C C LIMIT - INTEGER C GIVES AN UPPERBOUND 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 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 C THE ESTIMATES FOR RESULT AND ERROR ARE C LESS RELIABLE. IT IS ASSUMED THAT THE C REQUESTED 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 C OF LIMIT. C HOWEVER, IF THIS YIELDS NO IMPROVEMENT IT C IS RATHER ADVISED TO ANALYZE THE INTEGRAND C IN ORDER TO DETERMINE THE INTEGRATION C DIFFICULTIES. IF THE POSITION OF A LOCAL C DIFFICULTY CAN BE DETERMINED(E.G. C SINGULARITY, DISCONTINUITY WITHIN THE C INTERVAL) ONE WILL PROBABLY GAIN FROM C SPLITTING UP THE INTERVAL AT THIS POINT C AND CALLING THE INTEGRATOR ON THE C SUBRANGES. IF POSSIBLE, AN APPROPRIATE C SPECIAL-PURPOSE INTEGRATOR SHOULD BE USED C WHICH IS DESIGNED FOR HANDLING THE TYPE OF C DIFFICULTY INVOLVED. C = 2 THE OCCURRENCE OF ROUNDOFF ERROR IS C DETECTED, WHICH PREVENTS THE REQUESTED C TOLERANCE FROM BEING ACHIEVED. C = 3 EXTREMELY BAD INTEGRAND BEHAVIOUR OCCURS C AT SOME POINTS OF THE INTEGRATION C INTERVAL. 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, LAST, RLIST(1) , C ELIST(1) AND IORD(1) ARE SET TO ZERO. C ALIST(1) AND BLIST(1) ARE SET TO A AND B C 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 GIVEN INTEGRATION RANGE (A,B) 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 GIVEN INTEGRATION RANGE (A,B) C C RLIST - DOUBLE PRECISION C VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST C LAST ELEMENTS OF WHICH ARE THE C INTEGRAL 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 AT LEAST LIMIT, THE FIRST K C ELEMENTS OF WHICH ARE POINTERS TO THE C ERROR ESTIMATES OVER THE SUBINTERVALS, C SUCH THAT ELIST(IORD(1)), ..., C ELIST(IORD(K)) FORM A DECREASING SEQUENCE, C WITH K = LAST IF LAST.LE.(LIMIT/2+2), AND C K = LIMIT+1-LAST OTHERWISE C C LAST - INTEGER C NUMBER OF SUBINTERVALS ACTUALLY PRODUCED IN THE C SUBDIVISION PROCESS C C VALP - DOUBLE PRECISION C VALN ARRAYS OF DIMENSION AT LEAST (21,LIMIT) USED TO C SAVE THE FUNCTIONAL VALUES C C LP - INTEGER C LN VECTORS OF DIMENSION AT LEAST LIMIT, USED TO C STORE THE ACTUAL NUMBER OF FUNCTIONAL VALUES C SAVED IN THE CORRESPONDING COLUMN C OF VALP,VALN C C***ROUTINES CALLED F,D1MACH,DQXLQM,DQXRRD,DQPSRT,DQXCPY C DOUBLE PRECISION A,ABSERR,ALIST,AREA,AREA1,AREA12,AREA2,A1,A2,B, * BLIST,B1,B2,DABS,DEFABS,DEFAB1,DEFAB2,DMAX1,D1MACH,ELIST,EPMACH, * EPSABS,EPSREL,ERRBND,ERRMAX,ERROR1,ERROR2,ERRO12,ERRSUM,F, * RESABS,RESULT,RLIST,UFLOW,VALP,VALN,VP1,VP2,VN1,VN2 INTEGER KEY,IER,IORD,IROFF1,IROFF2,K,LAST,LIMIT,MAXERR, * NRMAX,LP,LN,LP1,LP2,LN1,LN2 C DIMENSION ALIST(LIMIT),BLIST(LIMIT),ELIST(LIMIT),IORD(LIMIT), * RLIST(LIMIT),VALP(21,LIMIT),VALN(21,LIMIT),LP(LIMIT),LN(LIMIT), * VP1(21),VP2(21),VN1(21),VN2(21) C EXTERNAL F 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 DQXGE EPMACH = D1MACH(4) UFLOW = D1MACH(1) C C TEST ON VALIDITY OF PARAMETERS C ------------------------------ C IER = 0 LAST = 0 RESULT = 0.0D+00 ABSERR = 0.0D+00 ALIST(1) = A BLIST(1) = B RLIST(1) = 0.0D+00 ELIST(1) = 0.0D+00 IORD(1) = 0 IF(EPSABS.LE.0.0D+00.AND. * EPSREL.LT.DMAX1(0.5D+02*EPMACH,0.5D-28)) IER = 6 IF(IER.EQ.6) GO TO 999 C C FIRST APPROXIMATION TO THE INTEGRAL C ----------------------------------- C LP(1)=1 LN(1)=1 VALP(1,1)=F((A+B)*0.5D0) VALN(1,1)=VALP(1,1) CALL DQXLQM(F,A,B,RESULT,ABSERR,DEFABS,RESABS, * VALP(1,1),VALN(1,1),LP(1),LN(1),KEY) LAST = 1 RLIST(1) = RESULT ELIST(1) = ABSERR IORD(1) = 1 C C TEST ON ACCURACY. C ERRBND = DMAX1(EPSABS,EPSREL*DABS(RESULT)) IF(ABSERR.LE.0.5D+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.ABSERR.EQ.0.0D+00) GO TO 999 C C INITIALIZATION C -------------- C C ERRMAX = ABSERR MAXERR = 1 AREA = RESULT ERRSUM = ABSERR NRMAX = 1 IROFF1 = 0 IROFF2 = 0 C C MAIN DO-LOOP C ------------ C DO 30 LAST = 2,LIMIT C C BISECT THE SUBINTERVAL WITH THE LARGEST ERROR ESTIMATE. C A1 = ALIST(MAXERR) B1 = 0.5D+00*(ALIST(MAXERR)+BLIST(MAXERR)) A2 = B1 B2 = BLIST(MAXERR) CALL DQXRRD(F,VALN(1,MAXERR),LN(MAXERR),B1,A1,VN1,VP1,LN1,LP1) CALL DQXLQM(F,A1,B1,AREA1,ERROR1,RESABS,DEFAB1,VP1,VN1,LP1,LN1, * KEY) CALL DQXRRD(F,VALP(1,MAXERR),LP(MAXERR),A2,B2,VP2,VN2,LP2,LN2) CALL DQXLQM(F,A2,B2,AREA2,ERROR2,RESABS,DEFAB2,VP2,VN2,LP2,LN2, * KEY) 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 5 IF(DABS(RLIST(MAXERR)-AREA12).LE.0.1D-04*DABS(AREA12) * .AND.ERRO12.GE.0.99D+00*ERRMAX) IROFF1 = IROFF1+1 IF(LAST.GT.10.AND.ERRO12.GT.ERRMAX) IROFF2 = IROFF2+1 5 RLIST(MAXERR) = AREA1 RLIST(LAST) = AREA2 ERRBND = DMAX1(EPSABS,EPSREL*DABS(AREA)) IF(ERRSUM.LE.ERRBND) GO TO 8 C C TEST FOR ROUNDOFF ERROR AND EVENTUALLY SET ERROR FLAG. C IF(IROFF1.GE.6.OR.IROFF2.GE.20) IER = 2 C C SET ERROR FLAG IN THE CASE THAT THE NUMBER OF SUBINTERVALS C EQUALS LIMIT. C IF(LAST.EQ.LIMIT) IER = 1 C C SET ERROR FLAG IN THE CASE OF BAD INTEGRAND BEHAVIOUR C AT A POINT OF THE INTEGRATION RANGE. C IF(DMAX1(DABS(A1),DABS(B2)).LE.(0.1D+01+0.1D+03* * EPMACH)*(DABS(A2)+0.1D+04*UFLOW)) IER = 3 C C APPEND THE NEWLY-CREATED INTERVALS TO THE LIST. C 8 IF(ERROR2.GT.ERROR1) GO TO 10 ALIST(LAST) = A2 BLIST(MAXERR) = B1 BLIST(LAST) = B2 ELIST(MAXERR) = ERROR1 ELIST(LAST) = ERROR2 CALL DQXCPY(VALP(1,MAXERR),VP1,LP1) LP(MAXERR)=LP1 CALL DQXCPY(VALN(1,MAXERR),VN1,LN1) LN(MAXERR)=LN1 CALL DQXCPY(VALP(1,LAST),VP2,LP2) LP(LAST)=LP2 CALL DQXCPY(VALN(1,LAST),VN2,LN2) LN(LAST)=LN2 GO TO 20 10 ALIST(MAXERR) = A2 ALIST(LAST) = A1 BLIST(LAST) = B1 RLIST(MAXERR) = AREA2 RLIST(LAST) = AREA1 ELIST(MAXERR) = ERROR2 ELIST(LAST) = ERROR1 CALL DQXCPY(VALP(1,MAXERR),VP2,LP2) LP(MAXERR)=LP2 CALL DQXCPY(VALN(1,MAXERR),VN2,LN2) LN(MAXERR)=LN2 CALL DQXCPY(VALP(1,LAST),VP1,LP1) LP(LAST)=LP1 CALL DQXCPY(VALN(1,LAST),VN1,LN1) LN(LAST)=LN1 C C CALL SUBROUTINE DQPSRT TO MAINTAIN THE DESCENDING ORDERING C IN THE LIST OF ERROR ESTIMATES AND SELECT THE SUBINTERVAL C WITH THE LARGEST ERROR ESTIMATE (TO BE BISECTED NEXT). C 20 CALL DQPSRT(LIMIT,LAST,MAXERR,ERRMAX,ELIST,IORD,NRMAX) C ***JUMP OUT OF DO-LOOP IF(IER.NE.0.OR.ERRSUM.LE.ERRBND) GO TO 40 30 CONTINUE C C COMPUTE FINAL RESULT. C --------------------- C 40 RESULT = 0.0D+00 DO 50 K=1,LAST RESULT = RESULT+RLIST(K) 50 CONTINUE ABSERR = ERRSUM 999 RETURN END *DECK DQXLQM SUBROUTINE DQXLQM(F,A,B,RESULT,ABSERR,RESABS,RESASC,VR,VS,LR,LS, * KEY) C C TO COMPUTE I = INTEGRAL OF F OVER (A,B), WITH ERROR C ESTIMATE C J = INTEGRAL OF ABS(F) OVER (A,B) 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 A - DOUBLE PRECISION C LOWER LIMIT OF INTEGRATION C C B - DOUBLE PRECISION C UPPER LIMIT OF INTEGRATION C C VR - DOUBLE PRECISION C VECTOR OF LENGTH LR CONTAINING THE C SAVED FUNCTIONAL VALUES OF POSITIVE ABSCISSAS C C VS - DOUBLE PRECISION C VECTOR OF LENGTH LS CONTAINING THE C SAVED FUNCTIONAL VALUES OF NEGATIVE ABSCISSAS C C LR - INTEGER C LS NUMBER OF ELEMENTS IN C VR,VS RESPECTIVELY C C KEY - INTEGER C KEY FOR CHOICE OF LOCAL INTEGRATION RULE C RMS FORMULAS ARE USED WITH C 13 - 19 POINTS IF KEY.LT.1, C 13 - 19 - (27) POINTS IF KEY = 1, C 13 - 19 - (27) - (41) POINTS IF KEY = 2, C 19 - 27 - (41) POINTS IF KEY = 3, C 27 - 41 POINTS IF KEY.GT.3. C C (RULES) USED IF THE FUNCTION APPEARS C ENOUGH REGULAR C C ON RETURN C RESULT - DOUBLE PRECISION C APPROXIMATION TO THE INTEGRAL I C C ABSERR - DOUBLE PRECISION C ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR, C WHICH SHOULD NOT 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 ABS(F-I/(B-A)) C OVER (A,B) C C VR - DOUBLE PRECISION C VECTOR OF LENGTH LR CONTAINING THE C SAVED FUNCTIONAL VALUES OF POSITIVE ABSCISSAS C C VS - DOUBLE PRECISION C VECTOR OF LENGTH LS CONTAINING THE C SAVED FUNCTIONAL VALUES OF NEGATIVE ABSCISSAS C C LR - INTEGER C LS NUMBER OF ELEMENTS IN C VR,VS RESPECTIVELY C C***ROUTINES CALLED D1MACH,DQXRUL C DOUBLE PRECISION F,A,B,RESULT,ABSERR,RESABS,RESASC, * D1MACH,EPMACH,RESG,RESK,UFLOW,ERROLD,VR(21),VS(21) INTEGER K,K0,K1,K2,KEY,KEY1,LR,LS EXTERNAL F C C MACHINE DEPENDENT CONSTANTS C --------------------------- C C EPMACH IS THE LARGEST RELATIVE SPACING. C UFLOW IS THE SMALLEST POSITIVE MAGNITUDE. C ERROLD IS THE LARGEST MAGNITUDE. C C***FIRST EXECUTABLE STATEMENT DQXLQM EPMACH = D1MACH(4) UFLOW = D1MACH(1) ERROLD = D1MACH(2) C KEY1 = MAX(KEY , 0) KEY1 = MIN(KEY1, 4) K0 = MAX(KEY1-2,0) K1 = K0+1 K2 = MIN(KEY1+1,3) C CALL DQXRUL(F,A,B,RESG,RESABS,RESASC,K0,K1,VR,VS,LR,LS) DO 99 K=K1,K2 CALL DQXRUL(F,A,B,RESK,RESABS,RESASC,K,K1,VR,VS,LR,LS) RESULT=RESK ABSERR = DABS((RESK-RESG)) IF(RESASC.NE.0.0D+00.AND.ABSERR.NE.0.0D+00) * ABSERR = RESASC*DMIN1(1.D0,(200D0*ABSERR/RESASC)**1.5D+00) IF(RESABS.GT.UFLOW/(10D0*EPMACH)) ABSERR = DMAX1 * ((EPMACH*10D0)*RESABS,ABSERR) RESG=RESK IF(ABSERR.GT.ERROLD*0.16)GOTO 3000 IF(ABSERR.LT.1000*EPMACH*RESABS)GOTO 3000 ERROLD=ABSERR 99 CONTINUE 3000 CONTINUE RETURN END *DECK DQXRUL SUBROUTINE DQXRUL(F,XL,XU,Y,YA,YM,KE,K1,FV1,FV2,L1,L2) C C TO COMPUTE I = INTEGRAL OF F OVER (A,B), WITH ERROR C ESTIMATE C AND CONDITIONALLY COMPUTE C J = INTEGRAL OF ABS(F) OVER (A,B) C BY USING AN RMS RULE 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 XL - DOUBLE PRECISION C LOWER LIMIT OF INTEGRATION C C XU - DOUBLE PRECISION C UPPER LIMIT OF INTEGRATION C C C KE - INTEGER C KEY FOR CHOICE OF LOCAL INTEGRATION RULE C AN RMS RULE IS USED WITH C 13 POINTS IF KE = 2, C 19 POINTS IF KE = 3, C 27 POINTS IF KE = 4, C 42 POINTS IF KE = 5 C C K1 INTEGER C VALUE OF KEY FOR WHICH THE ADDITIONAL ESTIMATES C YA, YM ARE TO BE COMPUTED C C FV1 - DOUBLE PRECISION C VECTOR CONTAINING L1 C SAVED FUNCTIONAL VALUES OF POSITIVE ABSCISSAS C C FV2 - DOUBLE PRECISION C VECTOR CONTAINING L2 C SAVED FUNCTIONAL VALUES OF NEGATIVE ABSCISSAS C C L1 - INTEGER C L2 NUMBER OF ELEMENTS IN FV1,FV2 RESPECTIVELY C C ON RETURN C Y - DOUBLE PRECISION C APPROXIMATION TO THE INTEGRAL I C RESULT IS COMPUTED BY APPLYING THE C REQUESTED RMS RULE C C YA - DOUBLE PRECISION C IF KEY = K1 APPROXIMATION TO THE INTEGRAL J C ELSE UNCHANGED C C YM - DOUBLE PRECISION C IF KEY = K1 APPROXIMATION TO THE INTEGRAL OF C ABS(F-I/(XU-XL) OVER (XL,XU) C ELSE UNCHANGED C C FV1 - DOUBLE PRECISION C VECTOR L1 CONTAINING L1 C SAVED FUNCTIONAL VALUES OF POSITIVE ABSCISSAS C C FV2 - DOUBLE PRECISION C VECTOR CONTAINING L2 C SAVED FUNCTIONAL VALUES OF NEGATIVE ABSCISSAS C C L1 - INTEGER C L2 NUMBER OF ELEMENTS IN FV1,FV2 RESPECTIVELY C C***ROUTINES CALLED F C DOUBLE PRECISION F,XL,XU,LDL,Y,YA,YM,Y2,XX(41),WW(52), * FV1(21),FV2(21),AA,BB,C EXTERNAL F INTEGER ISTART(4),LEN(4),KE,K1,L1,L2 SAVE ISTART,LEN,XX,WW DATA ISTART/0, 7, 17, 31/,LEN/7, 10, 14, 21/ DATA XX( 1)/.0 / DATA XX( 2)/.25000000000000000000D+00/ DATA XX( 3)/.50000000000000000000D+00/ DATA XX( 4)/.75000000000000000000D+00/ DATA XX( 5)/.87500000000000000000D+00/ DATA XX( 6)/.93750000000000000000D+00/ DATA XX( 7)/.10000000000000000000D+01/ DATA XX( 8)/.37500000000000000000D+00/ DATA XX( 9)/.62500000000000000000D+00/ DATA XX( 10)/.96875000000000000000D+00/ DATA XX( 11)/.12500000000000000000D+00/ DATA XX( 12)/.68750000000000000000D+00/ DATA XX( 13)/.81250000000000000000D+00/ DATA XX( 14)/.98437500000000000000D+00/ DATA XX( 15)/.18750000000000000000D+00/ DATA XX( 16)/.31250000000000000000D+00/ DATA XX( 17)/.43750000000000000000D+00/ DATA XX( 18)/.56250000000000000000D+00/ DATA XX( 19)/.84375000000000000000D+00/ DATA XX( 20)/.90625000000000000000D+00/ DATA XX( 21)/.99218750000000000000D+00/ C NUMBER OF NODES 13 DATA WW(1)/1.303262173284849021810473057638590518409112513421D-1/ DATA WW(2)/2.390632866847646220320329836544615917290026806242D-1/ DATA WW(3)/2.630626354774670227333506083741355715758124943143D-1/ DATA WW(4)/2.186819313830574175167853094864355208948886875898D-1/ DATA WW(5)/2.757897646642836865859601197607471574336674206700D-2/ DATA WW(6)/1.055750100538458443365034879086669791305550493830D-1/ DATA WW(7)/1.571194260595182254168429283636656908546309467968D-2/ C NUMBER OF NODES 19 DATA WW(8)/1.298751627936015783241173611320651866834051160074D-1/ DATA WW(9)/2.249996826462523640447834514709508786970828213187D-1/ DATA WW(15)/5.542699233295875168406783695143646338274805359780D-2/ DATA WW(10)/1.680415725925575286319046726692683040162290325505D-1/ DATA WW(16)/9.986735247403367525720377847755415293097913496236D-2/ DATA WW(11)/1.415567675701225879892811622832845252125600939627D-1/ DATA WW(12)/1.006482260551160175038684459742336605269707889822D-1/ DATA WW(13)/2.510604860724282479058338820428989444699235030871D-2/ DATA WW(17)/4.507523056810492466415880450799432587809828791196D-2/ DATA WW(14)/9.402964360009747110031098328922608224934320397592D-3/ C NUMBER OF NODES 27 DATA WW(18)/6.300942249647773931746170540321811473310938661469D-2/ DATA WW(28)/1.239572396231834242194189674243818619042280816640D-1/ DATA WW(19)/1.261383225537664703012999637242003647020326905948D-1/ DATA WW(25)/1.235837891364555000245004813294817451524633100256D-1/ DATA WW(20)/1.273864433581028272878709981850307363453523117880D-1/ DATA WW(26)/1.148933497158144016800199601785309838604146040215D-1/ DATA WW(29)/2.501306413750310579525950767549691151739047969345D-2/ DATA WW(21)/8.576500414311820514214087864326799153427368592787D-2/ DATA WW(30)/4.915957918146130094258849161350510503556792927578D-2/ DATA WW(22)/7.102884842310253397447305465997026228407227220665D-2/ DATA WW(23)/5.026383572857942403759829860675892897279675661654D-2/ DATA WW(27)/1.252575774226122633391477702593585307254527198070D-2/ DATA WW(31)/2.259167374956474713302030584548274729936249753832D-2/ DATA WW(24)/4.683670010609093810432609684738393586390722052124D-3/ C NUMBER OF NODES 41 DATA WW(32)/6.362762978782724559269342300509058175967124446839D-2/ DATA WW(42)/1.187141856692283347609436153545356484256869129472D-1/ DATA WW(46)/1.533126874056586959338368742803997744815413565014D-2/ DATA WW(33)/9.950065827346794643193261975720606296171462239514D-2/ DATA WW(47)/3.527159369750123100455704702965541866345781113903D-2/ DATA WW(39)/8.140326425945938045967829319725797511040878579808D-2/ DATA WW(48)/5.000556431653955124212795201196389006184693561679D-2/ DATA WW(34)/7.048220002718565366098742295389607994441704889441D-2/ DATA WW(49)/5.744164831179720106340717579281831675999717767532D-2/ DATA WW(40)/6.583213447600552906273539578430361199084485578379D-2/ DATA WW(43)/5.999947605385971985589674757013565610751028128731D-2/ DATA WW(35)/6.512297339398335645872697307762912795346716454337D-2/ DATA WW(44)/5.500937980198041736910257988346101839062581489820D-2/ DATA WW(50)/1.598823797283813438301248206397233634639162043386D-2/ DATA WW(36)/3.998229150313659724790527138690215186863915308702D-2/ DATA WW(51)/2.635660410220884993472478832884065450876913559421D-2/ DATA WW(37)/3.456512257080287509832054272964315588028252136044D-2/ DATA WW(41)/2.592913726450792546064232192976262988065252032902D-2/ DATA WW(45)/5.264422421764655969760271538981443718440340270116D-3/ DATA WW(52)/1.196003937945541091670106760660561117114584656319D-2/ DATA WW(38)/2.212167975884114432760321569298651047876071264944D-3/ C C***FIRST EXECUTABLE STATEMENT DQXRUL K=KE+1 IS=ISTART(K) KS=LEN(K) LDL= XU-XL BB= LDL*0.5D0 AA= XL+BB Y =0.D0 DO 100 I=1,KS IF(I.GT.L1.OR.I.GT.L2) C=BB*XX(I) IF(I.GT.L1) FV1(I)=F(AA+C) IF(I.GT.L2) FV2(I)=F(AA-C) 100 Y=Y+(FV1(I)+FV2(I))*WW(IS+I) Y2=Y Y=Y*BB IF(L1.LT.KS)L1=KS IF(L2.LT.KS)L2=KS IF(KE.NE.K1)RETURN YA=0.D0 DO 25 I=1,KS 25 YA=YA+(DABS(FV1(I))+DABS(FV2(I)))*WW(IS+I) Y2=Y2*0.5D0 YM=0.D0 YA=YA*DABS(BB) DO 27 I=1,KS 27 YM=YM+(DABS(FV1(I)-Y2)+DABS(FV2(I)-Y2))*WW(IS+I) YM=YM*DABS(BB) RETURN END *DECK DQXRRD SUBROUTINE DQXRRD(F,Z,LZ,XL,XU,R,S,LR,LS) C C TO REORDER THE COMPUTED FUNCTIONAL VALUES BEFORE C THE BISECTION OF AN INTERVAL 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 XL - DOUBLE PRECISION C LOWER LIMIT OF INTERVAL C C XU - DOUBLE PRECISION C UPPER LIMIT OF INTERVAL C C Z - DOUBLE PRECISION C VECTOR CONTAINING LZ C SAVED FUNCTIONAL VALUES C C LZ - INTEGER C NUMBER OF ELEMENTS IN LZ C C C ON RETURN C R - DOUBLE PRECISION C S VECTORS CONTAINING LR, LS C SAVED FUNCTIONAL VALUES FOR THE NEW INTERVALS C C LR - INTEGER C LS NUMBER OF ELEMENTES IN R,S RESPECTIVELY C C***ROUTINES CALLED F C DOUBLE PRECISION F,R,S,Z,XU,XL,DLEN,CENTR INTEGER LR,LS,LZ DIMENSION R(21),S(21),Z(21) C***FIRST EXECUTABLE STATEMENT DQXRRD DLEN= (XU-XL)*0.5D0 CENTR= XL+DLEN R(1)= Z(3) R(2)= Z(9) R(3)= Z(4) R(4)= Z(5) R(5)= Z(6) R(6)= Z(10) R(7)= Z(7) S(1)= Z(3) S(2)= Z(8) S(3)= Z(2) S(7)= Z(1) IF(LZ.GT.11)GOTO 2 R(8)= F(CENTR+DLEN*0.37500000D0) R(9)= F(CENTR+DLEN*0.62500000D0) R(10)= F(CENTR+DLEN*0.96875000D0) LR= 10 IF(LZ.NE.11)S(4)= F(CENTR-DLEN*0.75000000D0) IF(LZ.EQ.11)S(4)= Z(11) S(5)= F(CENTR-DLEN*0.87500000D0) S(6)= F(CENTR-DLEN*0.93750000D0) S(8)= F(CENTR-DLEN*0.37500000D0) S(9)= F(CENTR-DLEN*0.62500000D0) S(10)= F(CENTR-DLEN*0.96875000D0) LS= 10 GOTO 3000 2 R(8)= Z(12) R(9)= Z(13) R(10)= Z(14) LR= 10 S(4)= Z(11) S(5)= F(CENTR-DLEN*0.87500000D0) S(6)= F(CENTR-DLEN*0.93750000D0) IF(LZ.GT.14)GOTO3 S(8)= F(CENTR-DLEN*0.37500000D0) S(9)= F(CENTR-DLEN*0.62500000D0) S(10)= F(CENTR-DLEN*0.96875000D0) LS= 10 GOTO 3000 3 R(11)= Z(18) R(12)= Z(19) R(13)= Z(20) R(14)= Z(21) LR= 14 S(8)= Z(16) S(9)= Z(15) S(10)= F(CENTR-DLEN*0.96875000D0) S(11)= Z(17) LS= 11 3000 RETURN END *DECK DQXCPY SUBROUTINE DQXCPY(A,B,L) C C TO COPY THE DOUBLE PRECISION VECTOR B OF LENGTH L I N T O C THE DOUBLE PRECISION VECTOR A OF LENGTH L C C***REMARK THIS ROUTINE CAN BE IMPROVED, BY CODING IT IN THE C ASSEMBLER LANGUAGE OF THE USED MACHINE C***ROUTINES CALLED (NONE) C INTEGER L DOUBLE PRECISION A(L),B(L) C***FIRST EXECUTABLE STATEMENT DQXCPY DO 1 I=1,L 1 A(I)=B(I) RETURN END *DECK DTEST *** TEST PROGRAM *** C THIS DRIVER CALLS VARIOUS INTEGRATORS WITH VARIOUS RELATIVE C TOLERACES TO INTEGRATE THE FUNCTION C F(X) = DSQRT(X/(1-X)) * DLOG(X) OVER (0, 1) C C THE EXACT INTEGRAL IS -0.606789763508705... C DOUBLE PRECISION F,A,B,ATOL,RTOL,TEE,Y,WORK(4600) INTEGER NCALLS,IER,NEVAL,LAST,IWORK(300),KKK,KEY,LIMIT,LENW,LENIW EXTERNAL F COMMON /CTEST/NCALLS DATA LIMIT/100/,LENIW/300/,LENW/4600/ WRITE(6,600) 600 FORMAT(32H1 COMPUTATION OF THE INTEGRAL OF, * /38H F(X) = DSQRT(X/(1-X)) * DLOG(X), * /13H OVER (0, 1), * /23H THE EXACT INTEGRAL IS,30X, * 22H -0.606789763508705...) A=0.0D0 B=1.0D0 WRITE(6,61) 61 FORMAT(//22H TEST OF DQXGS // * 53H RELATIVE TOLERANCE FUNCTION CALLS ERROR ESTIMATE IER, * 13H RESULT ) DO 1 I=3,13,2 ATOL=1D-20 RTOL=10D0**(-I) NCALLS=0 CALL DQXGS (F,A,B,ATOL,RTOL,Y,TEE,IER,LIMIT, * LENIW,LENW,LAST,IWORK,WORK) WRITE(6,60)I,NCALLS,TEE,IER,Y 60 FORMAT(4X,6H10**(-,I2,1H),I17,9X,D7.2,4X,I2,F20.15) 1 CONTINUE DO 2 KKK = 1,5 KEY = KKK - 1 WRITE(6,62) KEY 62 FORMAT(//23H TEST OF DQXG , KEY = ,I4 // * 53H RELATIVE TOLERANCE FUNCTION CALLS ERROR ESTIMATE IER, * 13H RESULT ) DO 2 I=3,13,2 ATOL=1D-20 RTOL=10D0**(-I) NCALLS=0 CALL DQXG (F,A,B,ATOL,RTOL,KEY,Y,TEE,IER,LIMIT, * LENIW,LENW,LAST,IWORK,WORK) WRITE(6,60)I,NCALLS,TEE,IER,Y 2 CONTINUE STOP END *DECK F DOUBLE PRECISION FUNCTION F(X) DOUBLE PRECISION X C INTEGRAND FUNCTION, IN THE INTERVAL (0,1) THE EXACT INTEGRAL C TO 15 DIGITS IS C -0.60678 97635 08705D0 C COMMON /CTEST/NCALLS NCALLS=NCALLS+1 IF (X.EQ.0D0) GO TO 10 IF (X.EQ.1D0) GO TO 10 F=DSQRT(X/(1-X))*DLOG(X) RETURN 10 F=0.0 RETURN END CC+------------------------------------------------------------------+CC CC+ +CC CC+ LAWRENCE LIVERMORE NATIONAL LABORATORY +CC CC+ +CC CC+ MATHEMATICS AND STATISTICS DIVISION SOFTWARE LIBRARY +CC CC+ +CC CC+------------------------------------------------------------------+CC CC+ +CC CC+ THE SLATEC MATHEMATICAL PROGRAM LIBRARY IS ISSUED BY +CC CC+ THE FOLLOWING: +CC CC+ +CC CC+ AIR FORCE WEAPONS LABORATORY, ALBUQUERQUE +CC CC+ LAWRENCE LIVERMORE NATIONAL LABORATORY, LIVERMORE +CC CC+ LOS ALAMOS NATIONAL SCIENTIFIC LABORATORY, LOS ALAMOS +CC CC+ NATIONAL BUREAU OF STANDARDS, WASHINGTON +CC CC+ OAK RIDGE NATIONAL LABORATORY, OAK RIDGE +CC CC+ SANDIA NATIONAL LABORATORIES, ALBUQUERQUE +CC CC+ SANDIA NATIONAL LABORATORIES, LIVERMORE +CC CC+ +CC CC+ PLEASE REFER QUESTIONS REGARDING THIS SOFTWARE TO THE +CC CC+ MSD CONSULTING OFFICE, EXTENSION 3-2976. +CC CC+ +CC CC+ +CC CC+ * * * * * N O T I C E * * * * * +CC CC+ +CC CC+ THIS MATERIAL WAS PREPARED AS AN ACCOUNT OF WORK SPONSORED +CC CC+ BY THE UNITED STATES GOVERNMENT. NEITHER THE UNITED +CC CC+ STATES GOVERNMENT, NOR THE DEPARTMENT OF DEFENSE, NOR +CC CC+ ANY OF THEIR EMPLOYEES, MAKES ANY WARRANTY, EXPRESSED +CC CC+ OR IMPLIED, OR ASSUMES ANY LEGAL LIABILITY OR +CC CC+ RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR +CC CC+ USEFULNESS OF ANY INFORMATION, APPARATUS, PRODUCT, +CC CC+ OR PROCESS DISCLOSED, OR REPRESENTS ITS USE WOULD NOT +CC CC+ INFRINGE UPON PRIVATELY OWNED RIGHTS. +CC CC+ +CC CC+------------------------------------------------------------------+CC *DECK QPSRT SUBROUTINE QPSRT(LIMIT,LAST,MAXERR,ERMAX,ELIST,IORD,NRMAX) C***BEGIN PROLOGUE QPSRT C***REFER TO QAGE,QAGIE,QAGPE,QAGSE,QAWCE,QAWOE,QAWSE C***ROUTINES CALLED (NONE) C***KEYWORDS SEQUENTIAL SORTING C***DESCRIPTION C C 1. QPSRT C ORDERING ROUTINE C STANDARD FORTRAN SUBROUTINE C REAL VERSION C C 2. PURPOSE C THIS ROUTINE MAINTAINS THE DESCENDING ORDERING C IN THE LIST OF THE LOCAL ERROR ESTIMATES RESULTING FROM C THE 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 C AND BOTTOM-UP FOR THE SMALLEST ERROR ESTIMATE. C C 3. CALLING SEQUENCE C CALL QPSRT(LIMIT,LAST,MAXERR,ERMAX,ELIST,IORD,NRMAX) 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 C 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 - REAL C NRMAX-TH LARGEST ERROR ESTIMATE C ERMAX = ELIST(MAXERR) C C ELIST - REAL C VECTOR OF DIMENSION LAST CONTAINING C THE ERROR ESTIMATES C C IORD - INTEGER C VECTOR OF DIMENSION LAST, THE FIRST K C ELEMENTS OF WHICH CONTAIN POINTERS C TO THE ERROR 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 C 4. NO SUBROUTINES OR FUNCTIONS NEEDED 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 *DECK QELG SUBROUTINE QELG(N,EPSTAB,RESULT,ABSERR,RES3LA,NRES) C***BEGIN PROLOGUE QELG C***REFER TO QAGIE,QAGOE,QAGPE,QAGSE C***ROUTINES CALLED R1MACH C***REVISION DATE 830518 (YYMMDD) C***KEYWORDS CONVERGENCE ACCELERATION,EPSILON ALGORITHM,EXTRAPOLATION 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 DETERMINES THE LIMIT OF A GIVEN SEQUENCE OF C APPROXIMATIONS, BY MEANS OF THE EPSILON ALGORITHM OF C P. WYNN. AN ESTIMATE OF THE ABSOLUTE ERROR IS ALSO GIVEN. C THE CONDENSED EPSILON TABLE IS COMPUTED. ONLY THOSE C ELEMENTS NEEDED FOR THE COMPUTATION OF THE NEXT DIAGONAL C ARE PRESERVED. C***DESCRIPTION C C EPSILON ALGORITHM C STANDARD FORTRAN SUBROUTINE C REAL VERSION C C PARAMETERS C N - INTEGER C EPSTAB(N) CONTAINS THE NEW ELEMENT IN THE C FIRST COLUMN OF THE EPSILON TABLE. C C EPSTAB - REAL C VECTOR OF DIMENSION 52 CONTAINING THE ELEMENTS C OF THE TWO LOWER DIAGONALS OF THE TRIANGULAR C EPSILON TABLE. THE ELEMENTS ARE NUMBERED C STARTING AT THE RIGHT-HAND CORNER OF THE C TRIANGLE. C C RESULT - REAL C RESULTING APPROXIMATION TO THE INTEGRAL C C ABSERR - REAL C ESTIMATE OF THE ABSOLUTE ERROR COMPUTED FROM C RESULT AND THE 3 PREVIOUS RESULTS C C RES3LA - REAL C VECTOR OF DIMENSION 3 CONTAINING THE LAST 3 C RESULTS C C NRES - INTEGER C NUMBER OF CALLS TO THE ROUTINE C (SHOULD BE ZERO AT FIRST CALL) C***END PROLOGUE QELG C REAL ABSERR,DELTA1,DELTA2,DELTA3,R1MACH, 1 EPMACH,EPSINF,EPSTAB,ERROR,ERR1,ERR2,ERR3,E0,E1,E1ABS,E2,E3, 2 OFLOW,RES,RESULT,RES3LA,SS,TOL1,TOL2,TOL3 INTEGER I,IB,IB2,IE,INDX,K1,K2,K3,LIMEXP,N,NEWELM,NRES,NUM DIMENSION EPSTAB(52),RES3LA(3) C C LIST OF MAJOR VARIABLES C ----------------------- C C E0 - THE 4 ELEMENTS ON WHICH THE C E1 COMPUTATION OF A NEW ELEMENT IN C E2 THE EPSILON TABLE IS BASED C E3 E0 C E3 E1 NEW C E2 C NEWELM - NUMBER OF ELEMENTS TO BE COMPUTED IN THE NEW C DIAGONAL C ERROR - ERROR = ABS(E1-E0)+ABS(E2-E1)+ABS(NEW-E2) C RESULT - THE ELEMENT IN THE NEW DIAGONAL WITH LEAST VALUE C OF ERROR C C MACHINE DEPENDENT CONSTANTS C --------------------------- C C EPMACH IS THE LARGEST RELATIVE SPACING. C OFLOW IS THE LARGEST POSITIVE MAGNITUDE. C LIMEXP IS THE MAXIMUM NUMBER OF ELEMENTS THE EPSILON C TABLE CAN CONTAIN. IF THIS NUMBER IS REACHED, THE UPPER C DIAGONAL OF THE EPSILON TABLE IS DELETED. C C***FIRST EXECUTABLE STATEMENT QELG EPMACH = R1MACH(4) OFLOW = R1MACH(2) NRES = NRES+1 ABSERR = OFLOW RESULT = EPSTAB(N) IF(N.LT.3) GO TO 100 LIMEXP = 50 EPSTAB(N+2) = EPSTAB(N) NEWELM = (N-1)/2 EPSTAB(N) = OFLOW 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 E1ABS = ABS(E1) DELTA2 = E2-E1 ERR2 = ABS(DELTA2) TOL2 = AMAX1(ABS(E2),E1ABS)*EPMACH DELTA3 = E1-E0 ERR3 = ABS(DELTA3) TOL3 = AMAX1(E1ABS,ABS(E0))*EPMACH 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 C ***JUMP OUT OF DO-LOOP GO TO 100 10 E3 = EPSTAB(K1) EPSTAB(K1) = E1 DELTA1 = E1-E3 ERR1 = ABS(DELTA1) TOL1 = AMAX1(E1ABS,ABS(E3))*EPMACH 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 EPSINF = ABS(SS*E1) 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(EPSINF.GT.0.1E-03) GO TO 30 20 N = I+I-1 C ***JUMP OUT OF DO-LOOP 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 ERROR = ERR2+ABS(RES-E2)+ERR3 IF(ERROR.GT.ABSERR) GO TO 40 ABSERR = ERROR RESULT = RES 40 CONTINUE 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 INDX = NUM-N+1 DO 70 I = 1,N EPSTAB(I)= EPSTAB(INDX) INDX = INDX+1 70 CONTINUE 80 IF(NRES.GE.4) GO TO 90 RES3LA(NRES) = RESULT ABSERR = OFLOW GO TO 100 C C COMPUTE ERROR ESTIMATE C 90 ABSERR = ABS(RESULT-RES3LA(3))+ABS(RESULT-RES3LA(2)) 1 +ABS(RESULT-RES3LA(1)) RES3LA(1) = RES3LA(2) RES3LA(2) = RES3LA(3) RES3LA(3) = RESULT 100 ABSERR = AMAX1(ABSERR,0.5E+01*EPMACH*ABS(RESULT)) RETURN END *DECK R1MACH REAL FUNCTION R1MACH(I) C***BEGIN PROLOGUE R1MACH C***DATE WRITTEN 790101 (YYMMDD) C***REVISION DATE 890213 (YYMMDD) C***CATEGORY NO. R1 C***KEYWORDS LIBRARY=SLATEC,TYPE=SINGLE PRECISION(R1MACH-S D1MACH-D), C MACHINE CONSTANTS C***AUTHOR FOX, P. A., (BELL LABS) C HALL, A. D., (BELL LABS) C SCHRYER, N. L., (BELL LABS) C***PURPOSE RETURNS SINGLE PRECISION MACHINE DEPENDENT CONSTANTS C***DESCRIPTION C C R1MACH CAN BE USED TO OBTAIN MACHINE-DEPENDENT PARAMETERS C FOR THE LOCAL MACHINE ENVIRONMENT. IT IS A FUNCTION C SUBROUTINE WITH ONE (INPUT) ARGUMENT, AND CAN BE CALLED C AS FOLLOWS, FOR EXAMPLE C C A = R1MACH(I) C C WHERE I=1,...,5. THE (OUTPUT) VALUE OF A ABOVE IS C DETERMINED BY THE (INPUT) VALUE OF I. THE RESULTS FOR C VARIOUS VALUES OF I ARE DISCUSSED BELOW. C C SINGLE-PRECISION MACHINE CONSTANTS C R1MACH(1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE. C R1MACH(2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE. C R1MACH(3) = B**(-T), THE SMALLEST RELATIVE SPACING. C R1MACH(4) = B**(1-T), THE LARGEST RELATIVE SPACING. C R1MACH(5) = LOG10(B) C C ASSUME SINGLE PRECISION NUMBERS ARE REPRESENTED IN THE T-DIGIT, C BASE-B FORM C C SIGN (B**E)*( (X(1)/B) + ... + (X(T)/B**T) ) C C WHERE 0 .LE. X(I) .LT. B FOR I=1,...,T, 0 .LT. X(1), AND C EMIN .LE. E .LE. EMAX. C C THE VALUES OF B, T, EMIN AND EMAX ARE PROVIDED IN I1MACH AS C FOLLOWS: C I1MACH(10) = B, THE BASE. C I1MACH(11) = T, THE NUMBER OF BASE-B DIGITS. C I1MACH(12) = EMIN, THE SMALLEST EXPONENT E. C I1MACH(13) = EMAX, THE LARGEST EXPONENT E. C C TO ALTER THIS FUNCTION FOR A PARTICULAR ENVIRONMENT, C THE DESIRED SET OF DATA STATEMENTS SHOULD BE ACTIVATED BY C REMOVING THE C FROM COLUMN 1. ALSO, THE VALUES OF C R1MACH(1) - R1MACH(4) SHOULD BE CHECKED FOR CONSISTENCY C WITH THE LOCAL OPERATING SYSTEM. C C***REFERENCES FOX, P.A., HALL, A.D., SCHRYER, N.L, *FRAMEWORK FOR C A PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHE- C MATICAL SOFTWARE, VOL. 4, NO. 2, JUNE 1978, C PP. 177-188. C***ROUTINES CALLED XERROR C C THE CALL TO XERROR HAS BEEN COMMENTED TO GET A STAND ALONE ROUTINE C IT CAN BE ADAPTED TO FIT THE LOCAL ERROR HANDLING PROCEDURES C NOTE ADDED BY F. ROMANI 7/11/89 C C***END PROLOGUE R1MACH C INTEGER SMALL(2) INTEGER LARGE(2) INTEGER RIGHT(2) INTEGER DIVER(2) INTEGER LOG10(2) C REAL RMACH(5) SAVE RMACH C EQUIVALENCE (RMACH(1),SMALL(1)) EQUIVALENCE (RMACH(2),LARGE(1)) EQUIVALENCE (RMACH(3),RIGHT(1)) EQUIVALENCE (RMACH(4),DIVER(1)) EQUIVALENCE (RMACH(5),LOG10(1)) C C MACHINE CONSTANTS FOR THE AMIGA C ABSOFT FORTRAN COMPILER USING THE 68020/68881 COMPILER OPTION C C DATA SMALL(1) / Z'00800000' / C DATA LARGE(1) / Z'7F7FFFFF' / C DATA RIGHT(1) / Z'33800000' / C DATA DIVER(1) / Z'34000000' / C DATA LOG10(1) / Z'3E9A209B' / C C MACHINE CONSTANTS FOR THE AMIGA C ABSOFT FORTRAN COMPILER USING SOFTWARE FLOATING POINT C C DATA SMALL(1) / Z'00800000' / C DATA LARGE(1) / Z'7EFFFFFF' / C DATA RIGHT(1) / Z'33800000' / C DATA DIVER(1) / Z'34000000' / C DATA LOG10(1) / Z'3E9A209B' / C C MACHINE CONSTANTS FOR THE APOLLO C C DATA SMALL(1) / 16#00800000 / C DATA LARGE(1) / 16#7FFFFFFF / C DATA RIGHT(1) / 16#33800000 / C DATA DIVER(1) / 16#34000000 / C DATA LOG10(1) / 16#3E9A209B / C C MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM C C DATA RMACH(1) / Z400800000 / C DATA RMACH(2) / Z5FFFFFFFF / C DATA RMACH(3) / Z4E9800000 / C DATA RMACH(4) / Z4EA800000 / C DATA RMACH(5) / Z500E730E8 / C C MACHINE CONSTANTS FOR THE BURROUGHS 5700/6700/7700 SYSTEMS C C DATA RMACH(1) / O1771000000000000 / C DATA RMACH(2) / O0777777777777777 / C DATA RMACH(3) / O1311000000000000 / C DATA RMACH(4) / O1301000000000000 / C DATA RMACH(5) / O1157163034761675 / C C MACHINE CONSTANTS FOR THE CDC 170/180 SERIES USING NOS/VE C C DATA RMACH(1) / Z"3001800000000000" / C DATA RMACH(2) / Z"4FFEFFFFFFFFFFFE" / C DATA RMACH(3) / Z"3FD2800000000000" / C DATA RMACH(4) / Z"3FD3800000000000" / C DATA RMACH(5) / Z"3FFF9A209A84FBCF" / C C MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES C C DATA RMACH(1) / 00564000000000000000B / C DATA RMACH(2) / 37767777777777777776B / C DATA RMACH(3) / 16414000000000000000B / C DATA RMACH(4) / 16424000000000000000B / C DATA RMACH(5) / 17164642023241175720B / C C MACHINE CONSTANTS FOR THE CELERITY C1260 C C DATA SMALL(1) / Z'00800000' / C DATA LARGE(1) / Z'7F7FFFFF' / C DATA RIGHT(1) / Z'33800000' / C DATA DIVER(1) / Z'34000000' / C DATA LOG10(1) / Z'3E9A209B' / C C MACHINE CONSTANTS FOR THE CONVEX C-1 C C DATA SMALL(1) / '00800000'X / C DATA LARGE(1) / '7FFFFFFF'X / C DATA RIGHT(1) / '34800000'X / C DATA DIVER(1) / '35000000'X / C DATA LOG10(1) / '3F9A209B'X / C C MACHINE CONSTANTS FOR THE CRAY-1 C C DATA RMACH(1) / 200034000000000000000B / C DATA RMACH(2) / 577767777777777777776B / C DATA RMACH(3) / 377224000000000000000B / C DATA RMACH(4) / 377234000000000000000B / C DATA RMACH(5) / 377774642023241175720B / C C MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200 C C NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING CARD - C STATIC RMACH(5) C C DATA SMALL / 20K, 0 / C DATA LARGE / 77777K, 177777K / C DATA RIGHT / 35420K, 0 / C DATA DIVER / 36020K, 0 / C DATA LOG10 / 40423K, 42023K / C C MACHINE CONSTANTS FOR THE HARRIS 220 C C DATA SMALL(1), SMALL(2) / '20000000, '00000201 / C DATA LARGE(1), LARGE(2) / '37777777, '00000177 / C DATA RIGHT(1), RIGHT(2) / '20000000, '00000352 / C DATA DIVER(1), DIVER(2) / '20000000, '00000353 / C DATA LOG10(1), LOG10(2) / '23210115, '00000377 / C C MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES C C DATA RMACH(1) / O402400000000 / C DATA RMACH(2) / O376777777777 / C DATA RMACH(3) / O714400000000 / C DATA RMACH(4) / O716400000000 / C DATA RMACH(5) / O776464202324 / C C MACHINE CONSTANTS FOR THE HP 2100 C 3 WORD DOUBLE PRECISION WITH FTN4 C C DATA SMALL(1), SMALL(2) / 40000B, 1 / C DATA LARGE(1), LARGE(2) / 77777B, 177776B / C DATA RIGHT(1), RIGHT(2) / 40000B, 325B / C DATA DIVER(1), DIVER(2) / 40000B, 327B / C DATA LOG10(1), LOG10(2) / 46420B, 46777B / C C MACHINE CONSTANTS FOR THE HP 2100 C 4 WORD DOUBLE PRECISION WITH FTN4 C C DATA SMALL(1), SMALL(2) / 40000B, 1 / C DATA LARGE91), LARGE(2) / 77777B, 177776B / C DATA RIGHT(1), RIGHT(2) / 40000B, 325B / C DATA DIVER(1), DIVER(2) / 40000B, 327B / C DATA LOG10(1), LOG10(2) / 46420B, 46777B / C C MACHINE CONSTANTS FOR THE HP 9000 C C DATA SMALL(1) / 00004000000B / C DATA LARGE(1) / 17677777777B / C DATA RIGHT(1) / 06340000000B / C DATA DIVER(1) / 06400000000B / C DATA LOG10(1) / 07646420233B / C C MACHINE CONSTANTS FOR THE ELXSI 6400 C ASSUMING REAL*4 IS THE DEFAULT REAL C C DATA SMALL(1) / '00800000'X / C DATA LARGE(1) / '7F7FFFFF'X / C DATA RIGHT(1) / '33800000'X / C DATA DIVER(1) / '34000000'X / C DATA LOG10(1) / '3E9A209B'X / C C MACHINE CONSTANTS FOR THE IBM 360/370 SERIES, C THE XEROX SIGMA 5/7/9, THE SEL SYSTEMS 85/86 AND C THE PERKIN ELMER (INTERDATA) 7/32. C C DATA RMACH(1) / Z00100000 / C DATA RMACH(2) / Z7FFFFFFF / C DATA RMACH(3) / Z3B100000 / C DATA RMACH(4) / Z3C100000 / C DATA RMACH(5) / Z41134413 / C C MACHINE CONSTANTS FOR THE IBM PC C C DATA SMALL(1) / 1.18E-38 / C DATA LARGE(1) / 3.40E+38 / C DATA RIGHT(1) / 0.595E-07 / C DATA DIVER(1) / 1.19E-07 / C DATA LOG10(1) / 0.30102999566 / C C MACHINE CONSTANTS FOR THE PDP-10 (KA OR KI PROCESSOR) C C DATA RMACH(1) / "000400000000 / C DATA RMACH(2) / "377777777777 / C DATA RMACH(3) / "146400000000 / C DATA RMACH(4) / "147400000000 / C DATA RMACH(5) / "177464202324 / C C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING C 32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C C DATA SMALL(1) / 8388608 / C DATA LARGE(1) / 2147483647 / C DATA RIGHT(1) / 880803840 / C DATA DIVER(1) / 889192448 / C DATA LOG10(1) / 1067065499 / C C DATA RMACH(1) / O00040000000 / C DATA RMACH(2) / O17777777777 / C DATA RMACH(3) / O06440000000 / C DATA RMACH(4) / O06500000000 / C DATA RMACH(5) / O07746420233 / C C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING C 16-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C C DATA SMALL(1), SMALL(2) / 128, 0 / C DATA LARGE(1), LARGE(2) / 32767, -1 / C DATA RIGHT(1), RIGHT(2) / 13440, 0 / C DATA DIVER(1), DIVER(2) / 13568, 0 / C DATA LOG10(1), LOG10(2) / 16282, 8347 / C C DATA SMALL(1), SMALL(2) / O000200, O000000 / C DATA LARGE(1), LARGE(2) / O077777, O177777 / C DATA RIGHT(1), RIGHT(2) / O032200, O000000 / C DATA DIVER(1), DIVER(2) / O032400, O000000 / C DATA LOG10(1), LOG10(2) / O037632, O020233 / C C MACHINE CONSTANTS FOR THE SUN C C DATA SMALL(1) / Z'00800000' / C DATA LARGE(1) / Z'7F7FFFFF' / C DATA RIGHT(1) / Z'33800000' / C DATA DIVER(1) / Z'34000000' / C DATA LOG10(1) / Z'3E9A209B' / C C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES C C DATA RMACH(1) / O000400000000 / C DATA RMACH(2) / O377777777777 / C DATA RMACH(3) / O146400000000 / C DATA RMACH(4) / O147400000000 / C DATA RMACH(5) / O177464202324 / C C MACHINE CONSTANTS FOR THE VAX 11/780 C (EXPRESSED IN INTEGER AND HEXADECIMAL) C THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSTEMS C THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS C C DATA SMALL(1) / 128 / C DATA LARGE(1) / -32769 / C DATA RIGHT(1) / 13440 / C DATA DIVER(1) / 13568 / C DATA LOG10(1) / 547045274 / C C DATA SMALL(1) / Z00000080 / C DATA LARGE(1) / ZFFFF7FFF / C DATA RIGHT(1) / Z00003480 / C DATA DIVER(1) / Z00003500 / C DATA LOG10(1) / Z209B3F9A / C C MACHINE CONSTANTS FOR THE Z80 MICROPROCESSOR C C DATA SMALL(1), SMALL(2) / 0, 256/ C DATA LARGE(1), LARGE(2) / -1, -129/ C DATA RIGHT(1), RIGHT(2) / 0, 26880/ C DATA DIVER(1), DIVER(2) / 0, 27136/ C DATA LOG10(1), LOG10(2) / 8347, 32538/ C C C***FIRST EXECUTABLE STATEMENT R1MACH C C THE CALL TO XERROR HAS BEEN COMMENTED TO GET A STAND ALONE ROUTINE C IT CAN BE ADAPTED TO FIT THE LOCAL ERROR HANDLING PROCEDURES C NOTE ADDED BY F. ROMANI 7/11/89 C C IF (I .LT. 1 .OR. I .GT. 5) C 1 CALL XERROR ('R1MACH -- I OUT OF BOUNDS', 25, 1, 2) C R1MACH = RMACH(I) RETURN C END *DECK QXGS SUBROUTINE QXGS(F,A,B,EPSABS,EPSREL,RESULT,ABSERR,IER, 1 LIMIT,LENIW,LENW,LAST,IWORK,WORK) C C THE ROUTINE CALCULATES AN APPROXIMATION RESULT TO A GIVEN C DEFINITE INTEGRAL I = INTEGRAL OF F OVER (A,B), C HOPEFULLY SATISFYING FOLLOWING CLAIM FOR ACCURACY C ABS(I-RESULT).LE.MAX(EPSABS,EPSREL*ABS(I)). 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 A - REAL C LOWER LIMIT OF INTEGRATION C C B - REAL C UPPER LIMIT OF INTEGRATION 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.5E-28), C THE ROUTINE WILL END WITH IER = 6. 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 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 C THE ESTIMATES FOR INTEGRAL AND ERROR ARE C LESS RELIABLE. IT IS ASSUMED THAT THE C REQUESTED ACCURACY HAS NOT BEEN ACHIEVED. C ERROR MESSAGES C IER = 1 MAXIMUM NUMBER OF SUBDIVISIONS ALLOWED C HAS BEEN ACHIEVED. ONE CAN ALLOW MORE SUB- C DIVISIONS BY INCREASING THE VALUE OF LIMIT C (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 DETEC- C TED, WHICH PREVENTS THE REQUESTED C TOLERANCE FROM BEING ACHIEVED. C THE ERROR MAY BE UNDER-ESTIMATED. C = 3 EXTREMELY BAD INTEGRAND BEHAVIOUR C OCCURS 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. IT IS PRESUMED THAT C THE REQUESTED TOLERANCE CANNOT BE C ACHIEVED, AND THAT THE RETURNED RESULT IS C 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.5E-28) C OR LIMIT.LT.1 OR LENW.LT.LIMIT*46 OR C LENIW .LT. LIMIT*3 C RESULT, ABSERR, LAST ARE SET TO C ZERO.EXCEPT WHEN LIMIT OR LENW OR LENIW C IS 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 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*46 C IF LENW.LT.LIMIT*46, THE ROUTINE WILL END C WITH IER = 6. C C LENIW - INTEGER C DIMENSIONING PARAMETER FOR IWORK C LENIW MUST BE AT LEAST LIMIT*3 C IF LENW.LT.LIMIT*3, 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, DETERMINES THE C NUMBER OF SIGNIFICANT ELEMENTS ACTUALLY IN THE WORK C ARRAYS. C C WORK ARRAYS C IWORK - INTEGER C VECTOR OF DIMENSION AT LEAST 3*LIMIT, THE FIRST K C 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), C AND 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 C THE INTEGRAL APPROXIMATIONS OVER THE SUBINTERVALS, C WORK(LIMIT*3+1), ..., WORK(LIMIT*3+LAST) C CONTAIN THE ERROR ESTIMATES. C WORK(LIMIT*4+1), ... IS THE AREA RESERVED TO STORE C FUNCTIONAL VALUES. C***ROUTINES CALLED QXGSE,XERROR C C THE CALL TO XERROR HAS BEEN COMMENTED TO GET A STAND ALONE ROUTINE C IT CAN BE ADAPTED TO FIT THE LOCAL ERROR HANDLING PROCEDURES C REAL A,ABSERR,B,EPSABS,EPSREL,F,RESULT,WORK INTEGER IER,IWORK,LAST,LENIW,LENW,LIMIT,LVL,L1,L2,L3,L4,L5 C DIMENSION IWORK(LENIW),WORK(LENW) C EXTERNAL F C C CHECK VALIDITY OF LIMIT,LENIW AND LENW. C C***FIRST EXECUTABLE STATEMENT QXGS IER = 6 LAST = 0 RESULT = 0.0E+00 ABSERR = 0.0E+00 IF(LIMIT.LT.1.OR.LENIW.LT.LIMIT*3.OR.LENW.LT.LIMIT*46) GO TO 10 C C PREPARE CALL FOR QXGSE. C L1 = LIMIT+1 L2 = LIMIT+L1 L3 = LIMIT+L2 L4 = LIMIT+L3 L5 = 21*LIMIT+L3 C CALL QXGSE(F,A,B,EPSABS,EPSREL,LIMIT,RESULT,ABSERR, * IER,WORK(1),WORK(L1),WORK(L2),WORK(L3),IWORK,LAST, * WORK(L4),WORK(L5),IWORK(L1),IWORK(L2)) C C CALL ERROR HANDLER IF NECESSARY. C LVL = 0 10 IF(IER.EQ.6) LVL = 1 C C THE CALL TO XERROR HAS BEEN COMMENTED TO GET A STAND ALONE ROUTINE C IT CAN BE ADAPTED TO FIT THE LOCAL ERROR HANDLING PROCEDURES C C IF(IER.NE.0)CALL XERROR('ABNORMAL RETURN FROM QXGS ',28,IER,LVL) RETURN END *DECK QXGSE SUBROUTINE QXGSE(F,A,B,EPSABS,EPSREL,LIMIT,RESULT,ABSERR, * IER,ALIST,BLIST,RLIST,ELIST,IORD,LAST,VALP,VALN,LP,LN) C C THE ROUTINE CALCULATES AN APPROXIMATION RESULT TO A GIVEN C DEFINITE INTEGRAL I = INTEGRAL OF F OVER (A,B), C HOPEFULLY SATISFYING FOLLOWING CLAIM FOR ACCURACY C ABS(I-RESULT).LE.MAX(EPSABS,EPSREL*ABS(I)). 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 A - REAL C LOWER LIMIT OF INTEGRATION C C B - REAL C UPPER LIMIT OF INTEGRATION 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.5E-28), C THE ROUTINE WILL END WITH IER = 6. C C LIMIT - INTEGER C GIVES AN UPPERBOUND ON THE NUMBER OF SUBINTERVALS C IN THE PARTITION OF (A,B) 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 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 C THE ESTIMATES FOR INTEGRAL AND ERROR ARE C LESS RELIABLE. IT IS ASSUMED THAT THE C REQUESTED ACCURACY HAS NOT BEEN ACHIEVED. C ERROR MESSAGES C = 1 MAXIMUM NUMBER OF SUBDIVISIONS ALLOWED C HAS BEEN ACHIEVED. ONE CAN ALLOW MORE SUB- C DIVISIONS BY INCREASING THE VALUE OF LIMIT C (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 DETEC- C TED, WHICH PREVENTS THE REQUESTED C TOLERANCE FROM BEING ACHIEVED. C THE ERROR MAY BE UNDER-ESTIMATED. C = 3 EXTREMELY BAD INTEGRAND BEHAVIOUR C OCCURS 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 PRESUMED THAT THE REQUESTED C TOLERANCE CANNOT BE ACHIEVED, AND THAT THE C RETURNED RESULT IS THE BEST WHICH CAN BE C 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.5E-28). C RESULT, ABSERR, LAST, RLIST(1), C IORD(1) AND ELIST(1) ARE SET TO ZERO. C ALIST(1) AND BLIST(1) ARE SET TO A AND B C RESPECTIVELY. C C ALIST - REAL C VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST C LAST ELEMENTS OF WHICH ARE THE LEFT END POINTS C OF THE SUBINTERVALS IN THE PARTITION OF THE C GIVEN INTEGRATION RANGE (A,B) C C BLIST - REAL C VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST C LAST ELEMENTS OF WHICH ARE THE RIGHT END POINTS C OF THE SUBINTERVALS IN THE PARTITION OF THE GIVEN C INTEGRATION RANGE (A,B) C C RLIST - REAL 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 - REAL 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 AT LEAST 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 IN THE C SUBDIVISION PROCESS C C VALP - REAL C VALN ARRAYS OF DIMENSION AT LEAST (21,LIMIT) USED TO C SAVE THE FUNCTIONAL VALUES C C LP - INTEGER C LN VECTORS OF DIMENSION AT LEAST LIMIT, USED TO C STORE THE ACTUAL NUMBER OF FUNCTIONAL VALUES C SAVED IN THE CORRESPONDING COLUMN C OF VALP,VALN C C***ROUTINES CALLED F,R1MACH,QELG,QXLQM,QPSRT,QXRRD,QXCPY C REAL A,ABSEPS,ABSERR,ALIST,AREA,AREA1,AREA12,AREA2,A1, * A2,B,BLIST,B1,B2,CORREC,ABS,DEFABS,DEFAB1,DEFAB2,R1MACH,AMAX1, * DRES,ELIST,EPMACH,EPSABS,EPSREL,ERLARG,ERLAST,ERRBND,ERRMAX, * ERROR1,ERROR2,ERRO12,ERRSUM,ERTEST,F,OFLOW,RESABS,RESEPS,RESULT, * RES3LA,RLIST,RLIST2,SMALL,UFLOW, * VALP,VALN,VP1,VP2,VN1,VN2 INTEGER ID,IER,IERRO,IORD,IROFF1,IROFF2,IROFF3,JUPBND,K,KSGN, * KTMIN,LAST,LIMIT,MAXERR,NRES,NRMAX,NUMRL2, * LP,LN,LP1,LP2,LN1,LN2 LOGICAL EXTRAP,NOEXT DIMENSION ALIST(LIMIT),BLIST(LIMIT),ELIST(LIMIT),IORD(LIMIT), * RES3LA(3),RLIST(LIMIT),RLIST2(52), * VALP(21,LIMIT),VALN(21,LIMIT),LP(LIMIT),LN(LIMIT), * VP1(21),VP2(21),VN1(21),VN2(21) EXTERNAL F C C MACHINE DEPENDENT CONSTANTS C --------------------------- C C EPMACH IS THE LARGEST RELATIVE SPACING. C UFLOW IS THE SMALLEST POSITIVE MAGNITUDE. C OFLOW IS THE LARGEST POSITIVE MAGNITUDE. C C***FIRST EXECUTABLE STATEMENT QXGSE EPMACH = R1MACH(4) C C TEST ON VALIDITY OF PARAMETERS C ------------------------------ IER = 0 LAST = 0 RESULT = 0.0E+00 ABSERR = 0.0E+00 ALIST(1) = A BLIST(1) = B RLIST(1) = 0.0E+00 ELIST(1) = 0.0E+00 IF(EPSABS.LE.0.0E+00.AND.EPSREL.LT.AMAX1(0.5E+02*EPMACH,0.5E-28)) * IER = 6 IF(IER.EQ.6) GO TO 999 C C FIRST APPROXIMATION TO THE INTEGRAL C ----------------------------------- C UFLOW = R1MACH(1) OFLOW = R1MACH(2) IERRO = 0 LP(1)=1 LN(1)=1 VALP(1,1)=F((A+B)*0.5E0) VALN(1,1)=VALP(1,1) CALL QXLQM(F,A,B,RESULT,ABSERR,DEFABS,RESABS, * VALP(1,1),VALN(1,1),LP(1),LN(1), 2 ) C C TEST ON ACCURACY. C DRES = ABS(RESULT) ERRBND = AMAX1(EPSABS,EPSREL*DRES) LAST = 1 RLIST(1) = RESULT ELIST(1) = ABSERR IORD(1) = 1 IF(ABSERR.LE.1.0E+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. * ABSERR.EQ.0.0E+00) GO TO 999 C C INITIALIZATION C -------------- C RLIST2(1) = RESULT ERRMAX = ABSERR MAXERR = 1 AREA = RESULT ERRSUM = ABSERR ABSERR = OFLOW NRMAX = 1 NRES = 0 NUMRL2 = 2 KTMIN = 0 EXTRAP = .FALSE. NOEXT = .FALSE. 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 THE NRMAX-TH LARGEST ERROR C ESTIMATE. C A1 = ALIST(MAXERR) B1 = 0.5E+00*(ALIST(MAXERR)+BLIST(MAXERR)) A2 = B1 B2 = BLIST(MAXERR) ERLAST = ERRMAX CALL QXRRD(F,VALN(1,MAXERR),LN(MAXERR),B1,A1,VN1,VP1,LN1,LP1) CALL QXLQM(F,A1,B1,AREA1,ERROR1,RESABS,DEFAB1,VP1,VN1,LP1,LN1, * 2) CALL QXRRD(F,VALP(1,MAXERR),LP(MAXERR),A2,B2,VP2,VN2,LP2,LN2) CALL QXLQM(F,A2,B2,AREA2,ERROR2,RESABS,DEFAB2,VP2,VN2,LP2,LN2, * 2) 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) * .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 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 SUBINTERVALS C EQUALS LIMIT. C IF(LAST.EQ.LIMIT) IER = 1 C C SET ERROR FLAG IN THE CASE OF BAD INTEGRAND BEHAVIOUR C AT A POINT OF THE INTEGRATION RANGE. C IF(AMAX1(ABS(A1),ABS(B2)).LE.(0.1E+01+0.1E+03*EPMACH)* * (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 CALL QXCPY(VALP(1,MAXERR),VP1,LP1) LP(MAXERR)=LP1 CALL QXCPY(VALN(1,MAXERR),VN1,LN1) LN(MAXERR)=LN1 CALL QXCPY(VALP(1,LAST),VP2,LP2) LP(LAST)=LP2 CALL QXCPY(VALN(1,LAST),VN2,LN2) LN(LAST)=LN2 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 CALL QXCPY(VALP(1,MAXERR),VP2,LP2) LP(MAXERR)=LP2 CALL QXCPY(VALN(1,MAXERR),VN2,LN2) LN(MAXERR)=LN2 CALL QXCPY(VALP(1,LAST),VP1,LP1) LP(LAST)=LP1 CALL QXCPY(VALN(1,LAST),VN1,LN1) LN(LAST)=LN1 C C CALL SUBROUTINE QPSRT 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 QPSRT(LIMIT,LAST,MAXERR,ERRMAX,ELIST,IORD,NRMAX) C ***JUMP OUT OF DO-LOOP IF(ERRSUM.LE.ERRBND) GO TO 115 C ***JUMP OUT OF DO-LOOP 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 C C THE BOUND 0.3*ERTEST HAS BEEN INTRODUCED TO PERFORM A C MORE CAUTIOUS EXTRAPOLATION THAN IN THE ORIGINAL DQAGSE C ROUTINE C 40 IF(IERRO.EQ.3.OR.ERLARG.LE.0.3E0*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) C ***JUMP OUT OF DO-LOOP IF(ABS(BLIST(MAXERR)-ALIST(MAXERR)).GT.SMALL) GO TO 90 NRMAX = NRMAX+1 50 CONTINUE C C PERFORM EXTRAPOLATION. C 60 NUMRL2 = NUMRL2+1 RLIST2(NUMRL2) = AREA CALL QELG(NUMRL2,RLIST2,RESEPS,ABSEPS,RES3LA,NRES) KTMIN = KTMIN+1 IF(KTMIN.GT.5.AND.ABSERR.LT.0.1E-02*ERRSUM) IER = 5 IF(ABSEPS.GE.ABSERR) GO TO 70 KTMIN = 0 ABSERR = ABSEPS RESULT = RESEPS CORREC = ERLARG ERTEST = AMAX1(EPSABS,EPSREL*ABS(RESEPS)) C ***JUMP OUT OF DO-LOOP IF(ABSERR.LE.ERTEST) GO TO 100 C C PREPARE BISECTION OF THE SMALLEST INTERVAL. C 70 IF(NUMRL2.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 = ABS(B-A)*0.375E+00 ERLARG = ERRSUM ERTEST = ERRBND RLIST2(2) = AREA 90 CONTINUE C C SET FINAL RESULT AND ERROR ESTIMATE. C ------------------------------------ C 100 IF(ABSERR.EQ.OFLOW) 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. * DEFABS*0.1E-01) GO TO 130 IF(0.1E-01.GT.(RESULT/AREA).OR.(RESULT/AREA).GT.0.1E+03 * .OR.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 IF(IER.GT.2) IER = IER-1 999 RETURN END *DECK QXG SUBROUTINE QXG(F,A,B,EPSABS,EPSREL,KEY,RESULT,ABSERR,IER, 1 LIMIT,LENIW,LENW,LAST,IWORK,WORK) C C THE ROUTINE CALCULATES AN APPROXIMATION RESULT TO A GIVEN C DEFINITE INTEGRAL I = INTEGRAL OF F OVER (A,B), C HOPEFULLY SATISFYING FOLLOWING CLAIM FOR ACCURACY C ABS(I-RESULT).LE.MAX(EPSABS,EPSREL*ABS(I)). 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 A - REAL C LOWER LIMIT OF INTEGRATION C C B - REAL C UPPER LIMIT OF INTEGRATION 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.5E-28), C THE ROUTINE WILL END WITH IER = 6. C C KEY - INTEGER C KEY FOR CHOICE OF LOCAL INTEGRATION RULE C RMS FORMULAS ARE USED WITH C 13 - 19 POINTS IF KEY.LT.1, C 13 - 19 - (27) POINTS IF KEY = 1, C 13 - 19 - (27) - (41) POINTS IF KEY = 2, C 19 - 27 - (41) POINTS IF KEY = 3, C 27 - 41 POINTS IF KEY.GT.3. C C (RULES) USED IF THE FUNCTION APPEARS C ENOUGH REGULAR IN THE SUBINTERVAL 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 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 C THE ESTIMATES FOR RESULT AND ERROR ARE C LESS RELIABLE. IT IS ASSUMED THAT THE C REQUESTED ACCURACY HAS NOT BEEN ACHIEVED. C ERROR MESSAGES C IER = 1 MAXIMUM NUMBER OF SUBDIVISIONS ALLOWED C HAS BEEN ACHIEVED. ONE CAN ALLOW MORE SUB- C DIVISIONS BY INCREASING THE VALUE OF LIMIT C (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 DETEC- C TED, WHICH PREVENTS THE REQUESTED C TOLERANCE FROM BEING ACHIEVED. C THE ERROR MAY BE UNDER-ESTIMATED. C = 3 EXTREMELY BAD INTEGRAND BEHAVIOUR C OCCURS 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. IT IS PRESUMED THAT C THE REQUESTED TOLERANCE CANNOT BE C ACHIEVED, AND THAT THE RETURNED RESULT IS C 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.5E-28) C OR LIMIT.LT.1 OR LENW.LT.LIMIT*46 OR C LENIW .LT. LIMIT*3 C RESULT, ABSERR, LAST ARE SET TO C ZERO.EXCEPT WHEN LIMIT OR LENW OR LENIW C IS 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 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*46 C IF LENW.LT.LIMIT*46, THE ROUTINE WILL END C WITH IER = 6. C C LENIW - INTEGER C DIMENSIONING PARAMETER FOR IWORK C LENIW MUST BE AT LEAST LIMIT*3 C IF LENW.LT.LIMIT*3, 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, DETERMINES THE C NUMBER OF SIGNIFICANT ELEMENTS ACTUALLY IN THE WORK C ARRAYS. C C WORK ARRAYS C IWORK - INTEGER C VECTOR OF DIMENSION AT LEAST 3*LIMIT, THE FIRST K C 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), C AND 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 C THE INTEGRAL APPROXIMATIONS OVER THE SUBINTERVALS, C WORK(LIMIT*3+1), ..., WORK(LIMIT*3+LAST) C CONTAIN THE ERROR ESTIMATES. C WORK(LIMIT*4+1), ... IS THE AREA RESERVED TO STORE C FUNCTIONAL VALUES. C***ROUTINES CALLED QXGE,XERROR C C THE CALL TO XERROR HAS BEEN COMMENTED TO GET A STAND ALONE ROUTINE C IT CAN BE ADAPTED TO FIT THE LOCAL ERROR HANDLING PROCEDURES C C REAL A,ABSERR,B,EPSABS,EPSREL,F,RESULT,WORK INTEGER KEY,IER,IWORK,LAST,LENIW,LENW,LIMIT,LVL,L1,L2,L3,L4,L5 C DIMENSION IWORK(LENIW),WORK(LENW) C EXTERNAL F C C CHECK VALIDITY OF LIMIT,LENIW AND LENW. C C***FIRST EXECUTABLE STATEMENT QXG IER = 6 LAST = 0 RESULT = 0.0E+00 ABSERR = 0.0E+00 IF(LIMIT.LT.1.OR.LENIW.LT.LIMIT*3.OR.LENW.LT.LIMIT*46) GO TO 10 C C PREPARE CALL FOR QXGE. C L1 = LIMIT+1 L2 = LIMIT+L1 L3 = LIMIT+L2 L4 = LIMIT+L3 L5 = 21*LIMIT+L3 C CALL QXGE(F,A,B,EPSABS,EPSREL,KEY,LIMIT,RESULT,ABSERR, * IER,WORK(1),WORK(L1),WORK(L2),WORK(L3),IWORK,LAST, * WORK(L4),WORK(L5),IWORK(L1),IWORK(L2)) C C CALL ERROR HANDLER IF NECESSARY. C LVL = 0 10 IF(IER.EQ.6) LVL = 1 C C THE CALL TO XERROR HAS BEEN COMMENTED TO GET A STAND ALONE ROUTINE C IT CAN BE ADAPTED TO FIT THE LOCAL ERROR HANDLING PROCEDURES C C IF(IER.NE.0)CALL XERROR('ABNORMAL RETURN FROM QXG ',28,IER,LVL) RETURN END *DECK QXGE SUBROUTINE QXGE(F,A,B,EPSABS,EPSREL,KEY,LIMIT,RESULT,ABSERR, * IER,ALIST,BLIST,RLIST,ELIST,IORD,LAST,VALP,VALN,LP,LN) C C C THE ROUTINE CALCULATES AN APPROXIMATION RESULT TO A GIVEN C DEFINITE INTEGRAL I = INTEGRAL OF F OVER (A,B), C HOPEFULLY SATISFYING FOLLOWING CLAIM FOR ACCURACY C ABS(I-RESLT).LE.MAX(EPSABS,EPSREL*ABS(I)). 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 A - REAL C LOWER LIMIT OF INTEGRATION C C B - REAL C UPPER LIMIT OF INTEGRATION C C EPSABS - REAL C ABSOLUTE ACCURACY REQUESTED C C EPSREL - REAL C RELATIVE ACCURACY REQUESTED C IF EPSABS.LE.0 C AND EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5E-28), C THE ROUTINE WILL END WITH IER = 6. C C KEY - INTEGER C KEY FOR CHOICE OF LOCAL INTEGRATION RULE C RMS FORMULAS ARE USED WITH C 13 - 19 POINTS IF KEY.LT.1, C 13 - 19 - (27) POINTS IF KEY = 1, C 13 - 19 - (27) - (41) POINTS IF KEY = 2, C 19 - 27 - (41) POINTS IF KEY = 3, C 27 - 41 POINTS IF KEY.GT.3. C C (RULES) USED IF THE FUNCTION APPEARS C ENOUGH REGULAR IN THE SUBINTERVAL C C LIMIT - INTEGER C GIVES AN UPPERBOUND ON THE NUMBER OF SUBINTERVALS C IN THE PARTITION OF (A,B), LIMIT.GE.1. 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 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 C THE ESTIMATES FOR RESULT AND ERROR ARE C LESS RELIABLE. IT IS ASSUMED THAT THE C REQUESTED 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 C OF LIMIT. C HOWEVER, IF THIS YIELDS NO IMPROVEMENT IT C IS RATHER ADVISED TO ANALYZE THE INTEGRAND C IN ORDER TO DETERMINE THE INTEGRATION C DIFFICULTIES. IF THE POSITION OF A LOCAL C DIFFICULTY CAN BE DETERMINED(E.G. C SINGULARITY, DISCONTINUITY WITHIN THE C INTERVAL) ONE WILL PROBABLY GAIN FROM C SPLITTING UP THE INTERVAL AT THIS POINT C AND CALLING THE INTEGRATOR ON THE C SUBRANGES. IF POSSIBLE, AN APPROPRIATE C SPECIAL-PURPOSE INTEGRATOR SHOULD BE USED C WHICH IS DESIGNED FOR HANDLING THE TYPE OF C DIFFICULTY INVOLVED. C = 2 THE OCCURRENCE OF ROUNDOFF ERROR IS C DETECTED, WHICH PREVENTS THE REQUESTED C TOLERANCE FROM BEING ACHIEVED. C = 3 EXTREMELY BAD INTEGRAND BEHAVIOUR OCCURS C AT SOME POINTS OF THE INTEGRATION C INTERVAL. C = 6 THE INPUT IS INVALID, BECAUSE C (EPSABS.LE.0 AND C EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5E-28), C RESULT, ABSERR, LAST, RLIST(1) , C ELIST(1) AND IORD(1) ARE SET TO ZERO. C ALIST(1) AND BLIST(1) ARE SET TO A AND B C RESPECTIVELY. C C ALIST - REAL 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 GIVEN INTEGRATION RANGE (A,B) C C BLIST - REAL 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 GIVEN INTEGRATION RANGE (A,B) C C RLIST - REAL C VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST C LAST ELEMENTS OF WHICH ARE THE C INTEGRAL APPROXIMATIONS ON THE SUBINTERVALS C C ELIST - REAL 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 AT LEAST LIMIT, THE FIRST K C ELEMENTS OF WHICH ARE POINTERS TO THE C ERROR ESTIMATES OVER THE SUBINTERVALS, C SUCH THAT ELIST(IORD(1)), ..., C ELIST(IORD(K)) FORM A DECREASING SEQUENCE, C WITH K = LAST IF LAST.LE.(LIMIT/2+2), AND C K = LIMIT+1-LAST OTHERWISE C C LAST - INTEGER C NUMBER OF SUBINTERVALS ACTUALLY PRODUCED IN THE C SUBDIVISION PROCESS C C VALP - REAL C VALN ARRAYS OF DIMENSION AT LEAST (21,LIMIT) USED TO C SAVE THE FUNCTIONAL VALUES C C LP - INTEGER C LN VECTORS OF DIMENSION AT LEAST LIMIT, USED TO C STORE THE ACTUAL NUMBER OF FUNCTIONAL VALUES C SAVED IN THE CORRESPONDING COLUMN C OF VALP,VALN C C***ROUTINES CALLED F,R1MACH,QXLQM,QXRRD,QPSRT,QXCPY C REAL A,ABSERR,ALIST,AREA,AREA1,AREA12,AREA2,A1,A2,B, * BLIST,B1,B2,ABS,DEFABS,DEFAB1,DEFAB2,AMAX1,R1MACH,ELIST,EPMACH, * EPSABS,EPSREL,ERRBND,ERRMAX,ERROR1,ERROR2,ERRO12,ERRSUM,F, * RESABS,RESULT,RLIST,UFLOW,VALP,VALN,VP1,VP2,VN1,VN2 INTEGER KEY,IER,IORD,IROFF1,IROFF2,K,LAST,LIMIT,MAXERR, * NRMAX,LP,LN,LP1,LP2,LN1,LN2 C DIMENSION ALIST(LIMIT),BLIST(LIMIT),ELIST(LIMIT),IORD(LIMIT), * RLIST(LIMIT),VALP(21,LIMIT),VALN(21,LIMIT),LP(LIMIT),LN(LIMIT), * VP1(21),VP2(21),VN1(21),VN2(21) C EXTERNAL F 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 QXGE EPMACH = R1MACH(4) UFLOW = R1MACH(1) C C TEST ON VALIDITY OF PARAMETERS C ------------------------------ C IER = 0 LAST = 0 RESULT = 0.0E+00 ABSERR = 0.0E+00 ALIST(1) = A BLIST(1) = B RLIST(1) = 0.0E+00 ELIST(1) = 0.0E+00 IORD(1) = 0 IF(EPSABS.LE.0.0E+00.AND. * EPSREL.LT.AMAX1(0.5E+02*EPMACH,0.5E-28)) IER = 6 IF(IER.EQ.6) GO TO 999 C C FIRST APPROXIMATION TO THE INTEGRAL C ----------------------------------- C LP(1)=1 LN(1)=1 VALP(1,1)=F((A+B)*0.5E0) VALN(1,1)=VALP(1,1) CALL QXLQM(F,A,B,RESULT,ABSERR,DEFABS,RESABS, * VALP(1,1),VALN(1,1),LP(1),LN(1),KEY) LAST = 1 RLIST(1) = RESULT ELIST(1) = ABSERR IORD(1) = 1 C C TEST ON ACCURACY. C ERRBND = AMAX1(EPSABS,EPSREL*ABS(RESULT)) IF(ABSERR.LE.0.5E+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.ABSERR.EQ.0.0E+00) GO TO 999 C C INITIALIZATION C -------------- C C ERRMAX = ABSERR MAXERR = 1 AREA = RESULT ERRSUM = ABSERR NRMAX = 1 IROFF1 = 0 IROFF2 = 0 C C MAIN DO-LOOP C ------------ C DO 30 LAST = 2,LIMIT C C BISECT THE SUBINTERVAL WITH THE LARGEST ERROR ESTIMATE. C A1 = ALIST(MAXERR) B1 = 0.5E+00*(ALIST(MAXERR)+BLIST(MAXERR)) A2 = B1 B2 = BLIST(MAXERR) CALL QXRRD(F,VALN(1,MAXERR),LN(MAXERR),B1,A1,VN1,VP1,LN1,LP1) CALL QXLQM(F,A1,B1,AREA1,ERROR1,RESABS,DEFAB1,VP1,VN1,LP1,LN1, * KEY) CALL QXRRD(F,VALP(1,MAXERR),LP(MAXERR),A2,B2,VP2,VN2,LP2,LN2) CALL QXLQM(F,A2,B2,AREA2,ERROR2,RESABS,DEFAB2,VP2,VN2,LP2,LN2, * KEY) 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 5 IF(ABS(RLIST(MAXERR)-AREA12).LE.0.1E-04*ABS(AREA12) * .AND.ERRO12.GE.0.99E+00*ERRMAX) IROFF1 = IROFF1+1 IF(LAST.GT.10.AND.ERRO12.GT.ERRMAX) IROFF2 = IROFF2+1 5 RLIST(MAXERR) = AREA1 RLIST(LAST) = AREA2 ERRBND = AMAX1(EPSABS,EPSREL*ABS(AREA)) IF(ERRSUM.LE.ERRBND) GO TO 8 C C TEST FOR ROUNDOFF ERROR AND EVENTUALLY SET ERROR FLAG. C IF(IROFF1.GE.6.OR.IROFF2.GE.20) IER = 2 C C SET ERROR FLAG IN THE CASE THAT THE NUMBER OF SUBINTERVALS C EQUALS LIMIT. C IF(LAST.EQ.LIMIT) IER = 1 C C SET ERROR FLAG IN THE CASE OF BAD INTEGRAND BEHAVIOUR C AT A POINT OF THE INTEGRATION RANGE. C IF(AMAX1(ABS(A1),ABS(B2)).LE.(0.1E+01+0.1E+03* * EPMACH)*(ABS(A2)+0.1E+04*UFLOW)) IER = 3 C C APPEND THE NEWLY-CREATED INTERVALS TO THE LIST. C 8 IF(ERROR2.GT.ERROR1) GO TO 10 ALIST(LAST) = A2 BLIST(MAXERR) = B1 BLIST(LAST) = B2 ELIST(MAXERR) = ERROR1 ELIST(LAST) = ERROR2 CALL QXCPY(VALP(1,MAXERR),VP1,LP1) LP(MAXERR)=LP1 CALL QXCPY(VALN(1,MAXERR),VN1,LN1) LN(MAXERR)=LN1 CALL QXCPY(VALP(1,LAST),VP2,LP2) LP(LAST)=LP2 CALL QXCPY(VALN(1,LAST),VN2,LN2) LN(LAST)=LN2 GO TO 20 10 ALIST(MAXERR) = A2 ALIST(LAST) = A1 BLIST(LAST) = B1 RLIST(MAXERR) = AREA2 RLIST(LAST) = AREA1 ELIST(MAXERR) = ERROR2 ELIST(LAST) = ERROR1 CALL QXCPY(VALP(1,MAXERR),VP2,LP2) LP(MAXERR)=LP2 CALL QXCPY(VALN(1,MAXERR),VN2,LN2) LN(MAXERR)=LN2 CALL QXCPY(VALP(1,LAST),VP1,LP1) LP(LAST)=LP1 CALL QXCPY(VALN(1,LAST),VN1,LN1) LN(LAST)=LN1 C C CALL SUBROUTINE QPSRT TO MAINTAIN THE DESCENDING ORDERING C IN THE LIST OF ERROR ESTIMATES AND SELECT THE SUBINTERVAL C WITH THE LARGEST ERROR ESTIMATE (TO BE BISECTED NEXT). C 20 CALL QPSRT(LIMIT,LAST,MAXERR,ERRMAX,ELIST,IORD,NRMAX) C ***JUMP OUT OF DO-LOOP IF(IER.NE.0.OR.ERRSUM.LE.ERRBND) GO TO 40 30 CONTINUE C C COMPUTE FINAL RESULT. C --------------------- C 40 RESULT = 0.0E+00 DO 50 K=1,LAST RESULT = RESULT+RLIST(K) 50 CONTINUE ABSERR = ERRSUM 999 RETURN END *DECK QXLQM SUBROUTINE QXLQM(F,A,B,RESULT,ABSERR,RESABS,RESASC,VR,VS,LR,LS, * KEY) C C TO COMPUTE I = INTEGRAL OF F OVER (A,B), WITH ERROR C ESTIMATE C J = INTEGRAL OF ABS(F) OVER (A,B) 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 A - REAL C LOWER LIMIT OF INTEGRATION C C B - REAL C UPPER LIMIT OF INTEGRATION C C VR - REAL C VECTOR OF LENGTH LR CONTAINING THE C SAVED FUNCTIONAL VALUES OF POSITIVE ABSCISSAS C C VS - REAL C VECTOR OF LENGTH LS CONTAINING THE C SAVED FUNCTIONAL VALUES OF NEGATIVE ABSCISSAS C C LR - INTEGER C LS NUMBER OF ELEMENTS IN C VR,VS RESPECTIVELY C C KEY - INTEGER C KEY FOR CHOICE OF LOCAL INTEGRATION RULE C RMS FORMULAS ARE USED WITH C 13 - 19 POINTS IF KEY.LT.1, C 13 - 19 - (27) POINTS IF KEY = 1, C 13 - 19 - (27) - (41) POINTS IF KEY = 2, C 19 - 27 - (41) POINTS IF KEY = 3, C 27 - 41 POINTS IF KEY.GT.3. C C (RULES) USED IF THE FUNCTION APPEARS C ENOUGH REGULAR C C ON RETURN C RESULT - REAL C APPROXIMATION TO THE INTEGRAL I C C ABSERR - REAL C ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR, C WHICH SHOULD NOT EXCEED ABS(I-RESULT) C C RESABS - REAL C APPROXIMATION TO THE INTEGRAL J C C RESASC - REAL C APPROXIMATION TO THE INTEGRAL OF ABS(F-I/(B-A)) C OVER (A,B) C C VR - REAL C VECTOR OF LENGTH LR CONTAINING THE C SAVED FUNCTIONAL VALUES OF POSITIVE ABSCISSAS C C VS - REAL C VECTOR OF LENGTH LS CONTAINING THE C SAVED FUNCTIONAL VALUES OF NEGATIVE ABSCISSAS C C LR - INTEGER C LS NUMBER OF ELEMENTS IN C VR,VS RESPECTIVELY C C***ROUTINES CALLED R1MACH,QXRUL C REAL F,A,B,RESULT,ABSERR,RESABS,RESASC, * R1MACH,EPMACH,RESG,RESK,UFLOW,ERROLD,VR(21),VS(21) INTEGER K,K0,K1,K2,KEY,KEY1,LR,LS EXTERNAL F C C MACHINE DEPENDENT CONSTANTS C --------------------------- C C EPMACH IS THE LARGEST RELATIVE SPACING. C UFLOW IS THE SMALLEST POSITIVE MAGNITUDE. C ERROLD IS THE LARGEST MAGNITUDE. C C***FIRST EXECUTABLE STATEMENT QXLQM EPMACH = R1MACH(4) UFLOW = R1MACH(1) ERROLD = R1MACH(2) C KEY1 = MAX(KEY , 0) KEY1 = MIN(KEY1, 4) K0 = MAX(KEY1-2,0) K1 = K0+1 K2 = MIN(KEY1+1,3) C CALL QXRUL(F,A,B,RESG,RESABS,RESASC,K0,K1,VR,VS,LR,LS) DO 99 K=K1,K2 CALL QXRUL(F,A,B,RESK,RESABS,RESASC,K,K1,VR,VS,LR,LS) RESULT=RESK ABSERR = ABS((RESK-RESG)) IF(RESASC.NE.0.0E+00.AND.ABSERR.NE.0.0E+00) * ABSERR = RESASC*AMIN1(1.E0,(200E0*ABSERR/RESASC)**1.5E+00) IF(RESABS.GT.UFLOW/(10E0*EPMACH)) ABSERR = AMAX1 * ((EPMACH*10E0)*RESABS,ABSERR) RESG=RESK IF(ABSERR.GT.ERROLD*0.16)GOTO 3000 IF(ABSERR.LT.1000*EPMACH*RESABS)GOTO 3000 ERROLD=ABSERR 99 CONTINUE 3000 CONTINUE RETURN END *DECK QXRUL SUBROUTINE QXRUL(F,XL,XU,Y,YA,YM,KE,K1,FV1,FV2,L1,L2) C C TO COMPUTE I = INTEGRAL OF F OVER (A,B), WITH ERROR C ESTIMATE C AND CONDITIONALLY COMPUTE C J = INTEGRAL OF ABS(F) OVER (A,B) C BY USING AN RMS RULE 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 XL - REAL C LOWER LIMIT OF INTEGRATION C C XU - REAL C UPPER LIMIT OF INTEGRATION C C C KE - INTEGER C KEY FOR CHOICE OF LOCAL INTEGRATION RULE C AN RMS RULE IS USED WITH C 13 POINTS IF KE = 2, C 19 POINTS IF KE = 3, C 27 POINTS IF KE = 4, C 42 POINTS IF KE = 5 C C K1 INTEGER C VALUE OF KEY FOR WHICH THE ADDITIONAL ESTIMATES C YA, YM ARE TO BE COMPUTED C C FV1 - REAL C VECTOR CONTAINING L1 C SAVED FUNCTIONAL VALUES OF POSITIVE ABSCISSAS C C FV2 - REAL C VECTOR CONTAINING L2 C SAVED FUNCTIONAL VALUES OF NEGATIVE ABSCISSAS C C L1 - INTEGER C L2 NUMBER OF ELEMENTS IN FV1,FV2 RESPECTIVELY C C ON RETURN C Y - REAL C APPROXIMATION TO THE INTEGRAL I C RESULT IS COMPUTED BY APPLYING THE C REQUESTED RMS RULE C C YA - REAL C IF KEY = K1 APPROXIMATION TO THE INTEGRAL J C ELSE UNCHANGED C C YM - REAL C IF KEY = K1 APPROXIMATION TO THE INTEGRAL OF C ABS(F-I/(XU-XL) OVER (XL,XU) C ELSE UNCHANGED C C FV1 - REAL C VECTOR L1 CONTAINING L1 C SAVED FUNCTIONAL VALUES OF POSITIVE ABSCISSAS C C FV2 - REAL C VECTOR CONTAINING L2 C SAVED FUNCTIONAL VALUES OF NEGATIVE ABSCISSAS C C L1 - INTEGER C L2 NUMBER OF ELEMENTS IN FV1,FV2 RESPECTIVELY C C***ROUTINES CALLED F C REAL F,XL,XU,LDL,Y,YA,YM,Y2,XX(41),WW(52), * FV1(21),FV2(21),AA,BB,C EXTERNAL F INTEGER ISTART(4),LEN(4),KE,K1,L1,L2 SAVE ISTART,LEN,XX,WW DATA ISTART/0, 7, 17, 31/,LEN/7, 10, 14, 21/ DATA XX( 1)/.0 / DATA XX( 2)/.25000000000000000000E+00/ DATA XX( 3)/.50000000000000000000E+00/ DATA XX( 4)/.75000000000000000000E+00/ DATA XX( 5)/.87500000000000000000E+00/ DATA XX( 6)/.93750000000000000000E+00/ DATA XX( 7)/.10000000000000000000E+01/ DATA XX( 8)/.37500000000000000000E+00/ DATA XX( 9)/.62500000000000000000E+00/ DATA XX( 10)/.96875000000000000000E+00/ DATA XX( 11)/.12500000000000000000E+00/ DATA XX( 12)/.68750000000000000000E+00/ DATA XX( 13)/.81250000000000000000E+00/ DATA XX( 14)/.98437500000000000000E+00/ DATA XX( 15)/.18750000000000000000E+00/ DATA XX( 16)/.31250000000000000000E+00/ DATA XX( 17)/.43750000000000000000E+00/ DATA XX( 18)/.56250000000000000000E+00/ DATA XX( 19)/.84375000000000000000E+00/ DATA XX( 20)/.90625000000000000000E+00/ DATA XX( 21)/.99218750000000000000E+00/ C NUMBER OF NODES 13 DATA WW(1)/1.303262173284849021810473057638590518409112513421E-1/ DATA WW(2)/2.390632866847646220320329836544615917290026806242E-1/ DATA WW(3)/2.630626354774670227333506083741355715758124943143E-1/ DATA WW(4)/2.186819313830574175167853094864355208948886875898E-1/ DATA WW(5)/2.757897646642836865859601197607471574336674206700E-2/ DATA WW(6)/1.055750100538458443365034879086669791305550493830E-1/ DATA WW(7)/1.571194260595182254168429283636656908546309467968E-2/ C NUMBER OF NODES 19 DATA WW(8)/1.298751627936015783241173611320651866834051160074E-1/ DATA WW(9)/2.249996826462523640447834514709508786970828213187E-1/ DATA WW(15)/5.542699233295875168406783695143646338274805359780E-2/ DATA WW(10)/1.680415725925575286319046726692683040162290325505E-1/ DATA WW(16)/9.986735247403367525720377847755415293097913496236E-2/ DATA WW(11)/1.415567675701225879892811622832845252125600939627E-1/ DATA WW(12)/1.006482260551160175038684459742336605269707889822E-1/ DATA WW(13)/2.510604860724282479058338820428989444699235030871E-2/ DATA WW(17)/4.507523056810492466415880450799432587809828791196E-2/ DATA WW(14)/9.402964360009747110031098328922608224934320397592E-3/ C NUMBER OF NODES 27 DATA WW(18)/6.300942249647773931746170540321811473310938661469E-2/ DATA WW(28)/1.239572396231834242194189674243818619042280816640E-1/ DATA WW(19)/1.261383225537664703012999637242003647020326905948E-1/ DATA WW(25)/1.235837891364555000245004813294817451524633100256E-1/ DATA WW(20)/1.273864433581028272878709981850307363453523117880E-1/ DATA WW(26)/1.148933497158144016800199601785309838604146040215E-1/ DATA WW(29)/2.501306413750310579525950767549691151739047969345D-2/ DATA WW(21)/8.576500414311820514214087864326799153427368592787E-2/ DATA WW(30)/4.915957918146130094258849161350510503556792927578E-2/ DATA WW(22)/7.102884842310253397447305465997026228407227220665E-2/ DATA WW(23)/5.026383572857942403759829860675892897279675661654E-2/ DATA WW(27)/1.252575774226122633391477702593585307254527198070E-2/ DATA WW(31)/2.259167374956474713302030584548274729936249753832E-2/ DATA WW(24)/4.683670010609093810432609684738393586390722052124E-3/ C NUMBER OF NODES 41 DATA WW(32)/6.362762978782724559269342300509058175967124446839E-2/ DATA WW(42)/1.187141856692283347609436153545356484256869129472E-1/ DATA WW(46)/1.533126874056586959338368742803997744815413565014E-2/ DATA WW(33)/9.950065827346794643193261975720606296171462239514E-2/ DATA WW(47)/3.527159369750123100455704702965541866345781113903E-2/ DATA WW(39)/8.140326425945938045967829319725797511040878579808E-2/ DATA WW(48)/5.000556431653955124212795201196389006184693561679E-2/ DATA WW(34)/7.048220002718565366098742295389607994441704889441E-2/ DATA WW(49)/5.744164831179720106340717579281831675999717767532E-2/ DATA WW(40)/6.583213447600552906273539578430361199084485578379E-2/ DATA WW(43)/5.999947605385971985589674757013565610751028128731E-2/ DATA WW(35)/6.512297339398335645872697307762912795346716454337E-2/ DATA WW(44)/5.500937980198041736910257988346101839062581489820E-2/ DATA WW(50)/1.598823797283813438301248206397233634639162043386E-2/ DATA WW(36)/3.998229150313659724790527138690215186863915308702E-2/ DATA WW(51)/2.635660410220884993472478832884065450876913559421E-2/ DATA WW(37)/3.456512257080287509832054272964315588028252136044E-2/ DATA WW(41)/2.592913726450792546064232192976262988065252032902E-2/ DATA WW(45)/5.264422421764655969760271538981443718440340270116E-3/ DATA WW(52)/1.196003937945541091670106760660561117114584656319E-2/ DATA WW(38)/2.212167975884114432760321569298651047876071264944E-3/ C C***FIRST EXECUTABLE STATEMENT QXRUL K=KE+1 IS=ISTART(K) KS=LEN(K) LDL= XU-XL BB= LDL*0.5E0 AA= XL+BB Y =0.E0 DO 100 I=1,KS IF(I.GT.L1.OR.I.GT.L2) C=BB*XX(I) IF(I.GT.L1) FV1(I)=F(AA+C) IF(I.GT.L2) FV2(I)=F(AA-C) 100 Y=Y+(FV1(I)+FV2(I))*WW(IS+I) Y2=Y Y=Y*BB IF(L1.LT.KS)L1=KS IF(L2.LT.KS)L2=KS IF(KE.NE.K1)RETURN YA=0.E0 DO 25 I=1,KS 25 YA=YA+(ABS(FV1(I))+ABS(FV2(I)))*WW(IS+I) Y2=Y2*0.5E0 YM=0.E0 YA=YA*ABS(BB) DO 27 I=1,KS 27 YM=YM+(ABS(FV1(I)-Y2)+ABS(FV2(I)-Y2))*WW(IS+I) YM=YM*ABS(BB) RETURN END *DECK QXRRD SUBROUTINE QXRRD(F,Z,LZ,XL,XU,R,S,LR,LS) C C TO REORDER THE COMPUTED FUNCTIONAL VALUES BEFORE C THE BISECTION OF AN INTERVAL 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 XL - REAL C LOWER LIMIT OF INTERVAL C C XU - REAL C UPPER LIMIT OF INTERVAL C C Z - REAL C VECTOR CONTAINING LZ C SAVED FUNCTIONAL VALUES C C LZ - INTEGER C NUMBER OF ELEMENTS IN LZ C C C ON RETURN C R - REAL C S VECTORS CONTAINING LR, LS C SAVED FUNCTIONAL VALUES FOR THE NEW INTERVALS C C LR - INTEGER C LS NUMBER OF ELEMENTES IN R,S RESPECTIVELY C C***ROUTINES CALLED F C REAL F,R,S,Z,XU,XL,DLEN,CENTR INTEGER LR,LS,LZ DIMENSION R(21),S(21),Z(21) C***FIRST EXECUTABLE STATEMENT QXRRD DLEN= (XU-XL)*0.5E0 CENTR= XL+DLEN R(1)= Z(3) R(2)= Z(9) R(3)= Z(4) R(4)= Z(5) R(5)= Z(6) R(6)= Z(10) R(7)= Z(7) S(1)= Z(3) S(2)= Z(8) S(3)= Z(2) S(7)= Z(1) IF(LZ.GT.11)GOTO 2 R(8)= F(CENTR+DLEN*0.37500000E0) R(9)= F(CENTR+DLEN*0.62500000E0) R(10)= F(CENTR+DLEN*0.96875000E0) LR= 10 IF(LZ.NE.11)S(4)= F(CENTR-DLEN*0.75000000E0) IF(LZ.EQ.11)S(4)= Z(11) S(5)= F(CENTR-DLEN*0.87500000E0) S(6)= F(CENTR-DLEN*0.93750000E0) S(8)= F(CENTR-DLEN*0.37500000E0) S(9)= F(CENTR-DLEN*0.62500000E0) S(10)= F(CENTR-DLEN*0.96875000E0) LS= 10 GOTO 3000 2 R(8)= Z(12) R(9)= Z(13) R(10)= Z(14) LR= 10 S(4)= Z(11) S(5)= F(CENTR-DLEN*0.87500000E0) S(6)= F(CENTR-DLEN*0.93750000E0) IF(LZ.GT.14)GOTO3 S(8)= F(CENTR-DLEN*0.37500000E0) S(9)= F(CENTR-DLEN*0.62500000E0) S(10)= F(CENTR-DLEN*0.96875000E0) LS= 10 GOTO 3000 3 R(11)= Z(18) R(12)= Z(19) R(13)= Z(20) R(14)= Z(21) LR= 14 S(8)= Z(16) S(9)= Z(15) S(10)= F(CENTR-DLEN*0.96875000E0) S(11)= Z(17) LS= 11 3000 RETURN END *DECK QXCPY SUBROUTINE QXCPY(A,B,L) C C TO COPY THE REAL VECTOR B OF LENGTH L I N T O C THE REAL VECTOR A OF LENGTH L C C***REMARK THIS ROUTINE CAN BE IMPROVED, BY CODING IT IN THE C ASSEMBLER LANGUAGE OF THE USED MACHINE C***ROUTINES CALLED (NONE) C INTEGER L REAL A(L),B(L) C***FIRST EXECUTABLE STATEMENT QXCPY DO 1 I=1,L 1 A(I)=B(I) RETURN END *DECK MAIN *** TEST PROGRAM *** C THIS DRIVER CALLS VARIOUS INTEGRATORS WITH VARIOUS RELATIVE C TOLERACES TO INTEGRATE THE FUNCTION C F(X) = SQRT(X/(1-X)) * LOG(X) OVER (0, 1) C C THE EXACT INTEGRAL IS -0.606789763508705... C REAL F,A,B,ATOL,RTOL,TEE,Y,WORK(4600) INTEGER NCALLS,IER,NEVAL,LAST,IWORK(300),KKK,KEY,LIMIT,LENW,LENIW EXTERNAL F COMMON /CTEST/NCALLS DATA LIMIT/100/,LENIW/300/,LENW/4600/ WRITE(6,600) 600 FORMAT(32H1 COMPUTATION OF THE INTEGRAL OF, * /38H F(X) = SQRT(X/(1-X)) * LOG(X), * /13H OVER (0, 1), * /23H THE EXACT INTEGRAL IS,30X, * 22H -0.606789763508705...) A=0.0 B=1.0 WRITE(6,61) 61 FORMAT(//22H TEST OF QXGS // * 53H RELATIVE TOLERANCE FUNCTION CALLS ERROR ESTIMATE IER, * 13H RESULT ) DO 1 I=1,4 ATOL=1E-8 RTOL=10E0**(-I) NCALLS=0 CALL QXGS (F,A,B,ATOL,RTOL,Y,TEE,IER,LIMIT, * LENIW,LENW,LAST,IWORK,WORK) WRITE(6,60)I,NCALLS,TEE,IER,Y 60 FORMAT(4X,6H10**(-,I2,1H),I17,9X,E7.2,4X,I2,F20.8) 1 CONTINUE DO 2 KKK = 1,5 KEY = KKK - 1 WRITE(6,62) KEY 62 FORMAT(//23H TEST OF QXG , KEY = ,I4 // * 53H RELATIVE TOLERANCE FUNCTION CALLS ERROR ESTIMATE IER, * 13H RESULT ) DO 2 I=1,4 ATOL=1E-8 RTOL=10E0**(-I) NCALLS=0 CALL QXG (F,A,B,ATOL,RTOL,KEY,Y,TEE,IER,LIMIT, * LENIW,LENW,LAST,IWORK,WORK) WRITE(6,60)I,NCALLS,TEE,IER,Y 2 CONTINUE STOP END *DECK F REAL FUNCTION F(X) REAL X C INTEGRAND FUNCTION, IN THE INTERVAL (0,1) THE EXACT INTEGRAL C TO 15 DIGITS IS C -0.60678 97635 08705D0 C COMMON /CTEST/NCALLS NCALLS=NCALLS+1 IF (X.EQ.0.0) GO TO 10 IF (X.EQ.1.0) GO TO 10 F=SQRT(X/(1-X))*LOG(X) RETURN 10 F=0.0 RETURN END