*DECK DQNG SUBROUTINE DQNG (F, A, B, EPSABS, EPSREL, RESULT, ABSERR, NEVAL, + IER) C***BEGIN PROLOGUE DQNG C***PURPOSE The routine calculates an approximation result to a C given 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***LIBRARY SLATEC (QUADPACK) C***CATEGORY H2A1A1 C***TYPE DOUBLE PRECISION (QNG-S, DQNG-D) C***KEYWORDS AUTOMATIC INTEGRATOR, GAUSS-KRONROD(PATTERSON) RULES, C NONADAPTIVE, QUADPACK, QUADRATURE, SMOOTH INTEGRAND C***AUTHOR Piessens, Robert C Applied Mathematics and Programming Division C K. U. Leuven C de Doncker, Elise C Applied Mathematics and Programming Division C K. U. Leuven C***DESCRIPTION C C NON-ADAPTIVE INTEGRATION C STANDARD FORTRAN SUBROUTINE C DOUBLE PRECISION VERSION C C F - Double precision C Function subprogram defining the integrand function C F(X). The actual name for F needs to be declared C 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 I C Result is obtained by applying the 21-POINT C GAUSS-KRONROD RULE (RES21) obtained by optimal C addition of abscissae to the 10-POINT GAUSS RULE C (RES10), or by applying the 43-POINT RULE (RES43) C obtained by optimal addition of abscissae to the C 21-POINT GAUSS-KRONROD RULE, or by applying the C 87-POINT RULE (RES87) obtained by optimal addition C of abscissae to the 43-POINT RULE. C C ABSERR - Double precision C Estimate of the modulus of the absolute error, C which should EQUAL or EXCEED ABS(I-RESULT) C C NEVAL - Integer C Number of integrand evaluations C C IER - 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. It is C assumed that the requested accuracy has C not been achieved. C ERROR MESSAGES C IER = 1 The maximum number of steps has been C executed. The integral is probably too C difficult to be calculated by DQNG. 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 and NEVAL are set to zero. C C***REFERENCES (NONE) C***ROUTINES CALLED D1MACH, XERMSG C***REVISION HISTORY (YYMMDD) C 800101 DATE WRITTEN C 890531 Changed all specific intrinsics to generic. (WRB) C 890531 REVISION DATE from Version 3.2 C 891214 Prologue converted to Version 4.0 format. (BAB) C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) C***END PROLOGUE DQNG C DOUBLE PRECISION A,ABSC,ABSERR,B,CENTR,DHLGTH, 1 D1MACH,EPMACH,EPSABS,EPSREL,F,FCENTR,FVAL,FVAL1,FVAL2,FV1,FV2, 2 FV3,FV4,HLGTH,RESULT,RES10,RES21,RES43,RES87,RESABS,RESASC, 3 RESKH,SAVFUN,UFLOW,W10,W21A,W21B,W43A,W43B,W87A,W87B,X1,X2,X3,X4 INTEGER IER,IPX,K,L,NEVAL EXTERNAL F C DIMENSION FV1(5),FV2(5),FV3(5),FV4(5),X1(5),X2(5),X3(11),X4(22), 1 W10(5),W21A(5),W21B(6),W43A(10),W43B(12),W87A(21),W87B(23), 2 SAVFUN(21) C C THE FOLLOWING DATA STATEMENTS CONTAIN THE C ABSCISSAE AND WEIGHTS OF THE INTEGRATION RULES USED. C C X1 ABSCISSAE COMMON TO THE 10-, 21-, 43- AND 87- C POINT RULE C X2 ABSCISSAE COMMON TO THE 21-, 43- AND 87-POINT RULE C X3 ABSCISSAE COMMON TO THE 43- AND 87-POINT RULE C X4 ABSCISSAE OF THE 87-POINT RULE C W10 WEIGHTS OF THE 10-POINT FORMULA C W21A WEIGHTS OF THE 21-POINT FORMULA FOR ABSCISSAE X1 C W21B WEIGHTS OF THE 21-POINT FORMULA FOR ABSCISSAE X2 C W43A WEIGHTS OF THE 43-POINT FORMULA FOR ABSCISSAE X1, X3 C W43B WEIGHTS OF THE 43-POINT FORMULA FOR ABSCISSAE X3 C W87A WEIGHTS OF THE 87-POINT FORMULA FOR ABSCISSAE X1, C X2, X3 C W87B WEIGHTS OF THE 87-POINT FORMULA FOR ABSCISSAE X4 C C C GAUSS-KRONROD-PATTERSON QUADRATURE COEFFICIENTS FOR USE IN C QUADPACK ROUTINE QNG. THESE COEFFICIENTS WERE CALCULATED WITH C 101 DECIMAL DIGIT ARITHMETIC BY L. W. FULLERTON, BELL LABS, NOV 1981. C SAVE X1, W10, X2, W21A, W21B, X3, W43A, W43B, X4, W87A, W87B DATA X1 ( 1) / 0.9739065285 1717172007 7964012084 452 D0 / DATA X1 ( 2) / 0.8650633666 8898451073 2096688423 493 D0 / DATA X1 ( 3) / 0.6794095682 9902440623 4327365114 874 D0 / DATA X1 ( 4) / 0.4333953941 2924719079 9265943165 784 D0 / DATA X1 ( 5) / 0.1488743389 8163121088 4826001129 720 D0 / DATA W10 ( 1) / 0.0666713443 0868813759 3568809893 332 D0 / DATA W10 ( 2) / 0.1494513491 5058059314 5776339657 697 D0 / DATA W10 ( 3) / 0.2190863625 1598204399 5534934228 163 D0 / DATA W10 ( 4) / 0.2692667193 0999635509 1226921569 469 D0 / DATA W10 ( 5) / 0.2955242247 1475287017 3892994651 338 D0 / C DATA X2 ( 1) / 0.9956571630 2580808073 5527280689 003 D0 / DATA X2 ( 2) / 0.9301574913 5570822600 1207180059 508 D0 / DATA X2 ( 3) / 0.7808177265 8641689706 3717578345 042 D0 / DATA X2 ( 4) / 0.5627571346 6860468333 9000099272 694 D0 / DATA X2 ( 5) / 0.2943928627 0146019813 1126603103 866 D0 / DATA W21A ( 1) / 0.0325581623 0796472747 8818972459 390 D0 / DATA W21A ( 2) / 0.0750396748 1091995276 7043140916 190 D0 / DATA W21A ( 3) / 0.1093871588 0229764189 9210590325 805 D0 / DATA W21A ( 4) / 0.1347092173 1147332592 8054001771 707 D0 / DATA W21A ( 5) / 0.1477391049 0133849137 4841515972 068 D0 / DATA W21B ( 1) / 0.0116946388 6737187427 8064396062 192 D0 / DATA W21B ( 2) / 0.0547558965 7435199603 1381300244 580 D0 / DATA W21B ( 3) / 0.0931254545 8369760553 5065465083 366 D0 / DATA W21B ( 4) / 0.1234919762 6206585107 7958109831 074 D0 / DATA W21B ( 5) / 0.1427759385 7706008079 7094273138 717 D0 / DATA W21B ( 6) / 0.1494455540 0291690566 4936468389 821 D0 / C DATA X3 ( 1) / 0.9993333609 0193208139 4099323919 911 D0 / DATA X3 ( 2) / 0.9874334029 0808886979 5961478381 209 D0 / DATA X3 ( 3) / 0.9548079348 1426629925 7919200290 473 D0 / DATA X3 ( 4) / 0.9001486957 4832829362 5099494069 092 D0 / DATA X3 ( 5) / 0.8251983149 8311415084 7066732588 520 D0 / DATA X3 ( 6) / 0.7321483889 8930498261 2354848755 461 D0 / DATA X3 ( 7) / 0.6228479705 3772523864 1159120344 323 D0 / DATA X3 ( 8) / 0.4994795740 7105649995 2214885499 755 D0 / DATA X3 ( 9) / 0.3649016613 4658076804 3989548502 644 D0 / DATA X3 ( 10) / 0.2222549197 7660129649 8260928066 212 D0 / DATA X3 ( 11) / 0.0746506174 6138332204 3914435796 506 D0 / DATA W43A ( 1) / 0.0162967342 8966656492 4281974617 663 D0 / DATA W43A ( 2) / 0.0375228761 2086950146 1613795898 115 D0 / DATA W43A ( 3) / 0.0546949020 5825544214 7212685465 005 D0 / DATA W43A ( 4) / 0.0673554146 0947808607 5553166302 174 D0 / DATA W43A ( 5) / 0.0738701996 3239395343 2140695251 367 D0 / DATA W43A ( 6) / 0.0057685560 5976979618 4184327908 655 D0 / DATA W43A ( 7) / 0.0273718905 9324884208 1276069289 151 D0 / DATA W43A ( 8) / 0.0465608269 1042883074 3339154433 824 D0 / DATA W43A ( 9) / 0.0617449952 0144256449 6240336030 883 D0 / DATA W43A ( 10) / 0.0713872672 6869339776 8559114425 516 D0 / DATA W43B ( 1) / 0.0018444776 4021241410 0389106552 965 D0 / DATA W43B ( 2) / 0.0107986895 8589165174 0465406741 293 D0 / DATA W43B ( 3) / 0.0218953638 6779542810 2523123075 149 D0 / DATA W43B ( 4) / 0.0325974639 7534568944 3882222526 137 D0 / DATA W43B ( 5) / 0.0421631379 3519181184 7627924327 955 D0 / DATA W43B ( 6) / 0.0507419396 0018457778 0189020092 084 D0 / DATA W43B ( 7) / 0.0583793955 4261924837 5475369330 206 D0 / DATA W43B ( 8) / 0.0647464049 5144588554 4689259517 511 D0 / DATA W43B ( 9) / 0.0695661979 1235648452 8633315038 405 D0 / DATA W43B ( 10) / 0.0728244414 7183320815 0939535192 842 D0 / DATA W43B ( 11) / 0.0745077510 1417511827 3571813842 889 D0 / DATA W43B ( 12) / 0.0747221475 1740300559 4425168280 423 D0 / C DATA X4 ( 1) / 0.9999029772 6272923449 0529830591 582 D0 / DATA X4 ( 2) / 0.9979898959 8667874542 7496322365 960 D0 / DATA X4 ( 3) / 0.9921754978 6068722280 8523352251 425 D0 / DATA X4 ( 4) / 0.9813581635 7271277357 1916941623 894 D0 / DATA X4 ( 5) / 0.9650576238 5838461912 8284110607 926 D0 / DATA X4 ( 6) / 0.9431676131 3367059681 6416634507 426 D0 / DATA X4 ( 7) / 0.9158064146 8550720959 1826430720 050 D0 / DATA X4 ( 8) / 0.8832216577 7131650137 2117548744 163 D0 / DATA X4 ( 9) / 0.8457107484 6241566660 5902011504 855 D0 / DATA X4 ( 10) / 0.8035576580 3523098278 8739474980 964 D0 / DATA X4 ( 11) / 0.7570057306 8549555832 8942793432 020 D0 / DATA X4 ( 12) / 0.7062732097 8732181982 4094274740 840 D0 / DATA X4 ( 13) / 0.6515894665 0117792253 4422205016 736 D0 / DATA X4 ( 14) / 0.5932233740 5796108887 5273770349 144 D0 / DATA X4 ( 15) / 0.5314936059 7083193228 5268948562 671 D0 / DATA X4 ( 16) / 0.4667636230 4202284487 1966781659 270 D0 / DATA X4 ( 17) / 0.3994248478 5921880473 2101665817 923 D0 / DATA X4 ( 18) / 0.3298748771 0618828826 5053371824 597 D0 / DATA X4 ( 19) / 0.2585035592 0216155180 2280975429 025 D0 / DATA X4 ( 20) / 0.1856953965 6834665201 5917141167 606 D0 / DATA X4 ( 21) / 0.1118422131 7990746817 2398359241 362 D0 / DATA X4 ( 22) / 0.0373521233 9461987081 4998165437 704 D0 / DATA W87A ( 1) / 0.0081483773 8414917290 0002878448 190 D0 / DATA W87A ( 2) / 0.0187614382 0156282224 3935059003 794 D0 / DATA W87A ( 3) / 0.0273474510 5005228616 1582829741 283 D0 / DATA W87A ( 4) / 0.0336777073 1163793004 6581056957 588 D0 / DATA W87A ( 5) / 0.0369350998 2042790761 4589586742 499 D0 / DATA W87A ( 6) / 0.0028848724 3021153050 1334156248 695 D0 / DATA W87A ( 7) / 0.0136859460 2271270188 8950035273 128 D0 / DATA W87A ( 8) / 0.0232804135 0288831112 3409291030 404 D0 / DATA W87A ( 9) / 0.0308724976 1171335867 5466394126 442 D0 / DATA W87A ( 10) / 0.0356936336 3941877071 9351355457 044 D0 / DATA W87A ( 11) / 0.0009152833 4520224136 0843392549 948 D0 / DATA W87A ( 12) / 0.0053992802 1930047136 7738743391 053 D0 / DATA W87A ( 13) / 0.0109476796 0111893113 4327826856 808 D0 / DATA W87A ( 14) / 0.0162987316 9678733526 2665703223 280 D0 / DATA W87A ( 15) / 0.0210815688 8920383511 2433060188 190 D0 / DATA W87A ( 16) / 0.0253709697 6925382724 3467999831 710 D0 / DATA W87A ( 17) / 0.0291896977 5647575250 1446154084 920 D0 / DATA W87A ( 18) / 0.0323732024 6720278968 5788194889 595 D0 / DATA W87A ( 19) / 0.0347830989 5036514275 0781997949 596 D0 / DATA W87A ( 20) / 0.0364122207 3135178756 2801163687 577 D0 / DATA W87A ( 21) / 0.0372538755 0304770853 9592001191 226 D0 / DATA W87B ( 1) / 0.0002741455 6376207235 0016527092 881 D0 / DATA W87B ( 2) / 0.0018071241 5505794294 8341311753 254 D0 / DATA W87B ( 3) / 0.0040968692 8275916486 4458070683 480 D0 / DATA W87B ( 4) / 0.0067582900 5184737869 9816577897 424 D0 / DATA W87B ( 5) / 0.0095499576 7220164653 6053581325 377 D0 / DATA W87B ( 6) / 0.0123294476 5224485369 4626639963 780 D0 / DATA W87B ( 7) / 0.0150104473 4638895237 6697286041 943 D0 / DATA W87B ( 8) / 0.0175489679 8624319109 9665352925 900 D0 / DATA W87B ( 9) / 0.0199380377 8644088820 2278192730 714 D0 / DATA W87B ( 10) / 0.0221949359 6101228679 6332102959 499 D0 / DATA W87B ( 11) / 0.0243391471 2600080547 0360647041 454 D0 / DATA W87B ( 12) / 0.0263745054 1483920724 1503786552 615 D0 / DATA W87B ( 13) / 0.0282869107 8877120065 9968002987 960 D0 / DATA W87B ( 14) / 0.0300525811 2809269532 2521110347 341 D0 / DATA W87B ( 15) / 0.0316467513 7143992940 4586051078 883 D0 / DATA W87B ( 16) / 0.0330504134 1997850329 0785944862 689 D0 / DATA W87B ( 17) / 0.0342550997 0422606178 7082821046 821 D0 / DATA W87B ( 18) / 0.0352624126 6015668103 3782717998 428 D0 / DATA W87B ( 19) / 0.0360769896 2288870118 5500318003 895 D0 / DATA W87B ( 20) / 0.0366986044 9845609449 8018047441 094 D0 / DATA W87B ( 21) / 0.0371205492 6983257611 4119958413 599 D0 / DATA W87B ( 22) / 0.0373342287 5193504032 1235449094 698 D0 / DATA W87B ( 23) / 0.0373610737 6267902341 0321241766 599 D0 / C C LIST OF MAJOR VARIABLES C ----------------------- C C CENTR - MID POINT OF THE INTEGRATION INTERVAL C HLGTH - HALF-LENGTH OF THE INTEGRATION INTERVAL C FCENTR - FUNCTION VALUE AT MID POINT C ABSC - ABSCISSA C FVAL - FUNCTION VALUE C SAVFUN - ARRAY OF FUNCTION VALUES WHICH HAVE ALREADY BEEN C COMPUTED C RES10 - 10-POINT GAUSS RESULT C RES21 - 21-POINT KRONROD RESULT C RES43 - 43-POINT RESULT C RES87 - 87-POINT RESULT C RESABS - APPROXIMATION TO THE INTEGRAL OF ABS(F) C RESASC - APPROXIMATION TO THE INTEGRAL OF ABS(F-I/(B-A)) C C MACHINE DEPENDENT CONSTANTS C --------------------------- C C EPMACH IS THE LARGEST RELATIVE SPACING. C UFLOW IS THE SMALLEST POSITIVE MAGNITUDE. C C***FIRST EXECUTABLE STATEMENT DQNG EPMACH = D1MACH(4) UFLOW = D1MACH(1) C C TEST ON VALIDITY OF PARAMETERS C ------------------------------ C RESULT = 0.0D+00 ABSERR = 0.0D+00 NEVAL = 0 IER = 6 IF(EPSABS.LE.0.0D+00.AND.EPSREL.LT.MAX(0.5D+02*EPMACH,0.5D-28)) 1 GO TO 80 HLGTH = 0.5D+00*(B-A) DHLGTH = ABS(HLGTH) CENTR = 0.5D+00*(B+A) FCENTR = F(CENTR) NEVAL = 21 IER = 1 C C COMPUTE THE INTEGRAL USING THE 10- AND 21-POINT FORMULA. C DO 70 L = 1,3 GO TO (5,25,45),L 5 RES10 = 0.0D+00 RES21 = W21B(6)*FCENTR RESABS = W21B(6)*ABS(FCENTR) DO 10 K=1,5 ABSC = HLGTH*X1(K) FVAL1 = F(CENTR+ABSC) FVAL2 = F(CENTR-ABSC) FVAL = FVAL1+FVAL2 RES10 = RES10+W10(K)*FVAL RES21 = RES21+W21A(K)*FVAL RESABS = RESABS+W21A(K)*(ABS(FVAL1)+ABS(FVAL2)) SAVFUN(K) = FVAL FV1(K) = FVAL1 FV2(K) = FVAL2 10 CONTINUE IPX = 5 DO 15 K=1,5 IPX = IPX+1 ABSC = HLGTH*X2(K) FVAL1 = F(CENTR+ABSC) FVAL2 = F(CENTR-ABSC) FVAL = FVAL1+FVAL2 RES21 = RES21+W21B(K)*FVAL RESABS = RESABS+W21B(K)*(ABS(FVAL1)+ABS(FVAL2)) SAVFUN(IPX) = FVAL FV3(K) = FVAL1 FV4(K) = FVAL2 15 CONTINUE C C TEST FOR CONVERGENCE. C RESULT = RES21*HLGTH RESABS = RESABS*DHLGTH RESKH = 0.5D+00*RES21 RESASC = W21B(6)*ABS(FCENTR-RESKH) DO 20 K = 1,5 RESASC = RESASC+W21A(K)*(ABS(FV1(K)-RESKH)+ABS(FV2(K)-RESKH)) 1 +W21B(K)*(ABS(FV3(K)-RESKH)+ABS(FV4(K)-RESKH)) 20 CONTINUE ABSERR = ABS((RES21-RES10)*HLGTH) RESASC = RESASC*DHLGTH GO TO 65 C C COMPUTE THE INTEGRAL USING THE 43-POINT FORMULA. C 25 RES43 = W43B(12)*FCENTR NEVAL = 43 DO 30 K=1,10 RES43 = RES43+SAVFUN(K)*W43A(K) 30 CONTINUE DO 40 K=1,11 IPX = IPX+1 ABSC = HLGTH*X3(K) FVAL = F(ABSC+CENTR)+F(CENTR-ABSC) RES43 = RES43+FVAL*W43B(K) SAVFUN(IPX) = FVAL 40 CONTINUE C C TEST FOR CONVERGENCE. C RESULT = RES43*HLGTH ABSERR = ABS((RES43-RES21)*HLGTH) GO TO 65 C C COMPUTE THE INTEGRAL USING THE 87-POINT FORMULA. C 45 RES87 = W87B(23)*FCENTR NEVAL = 87 DO 50 K=1,21 RES87 = RES87+SAVFUN(K)*W87A(K) 50 CONTINUE DO 60 K=1,22 ABSC = HLGTH*X4(K) RES87 = RES87+W87B(K)*(F(ABSC+CENTR)+F(CENTR-ABSC)) 60 CONTINUE RESULT = RES87*HLGTH ABSERR = ABS((RES87-RES43)*HLGTH) 65 IF(RESASC.NE.0.0D+00.AND.ABSERR.NE.0.0D+00) 1 ABSERR = RESASC*MIN(0.1D+01,(0.2D+03*ABSERR/RESASC)**1.5D+00) IF (RESABS.GT.UFLOW/(0.5D+02*EPMACH)) ABSERR = MAX 1 ((EPMACH*0.5D+02)*RESABS,ABSERR) IF (ABSERR.LE.MAX(EPSABS,EPSREL*ABS(RESULT))) IER = 0 C ***JUMP OUT OF DO-LOOP IF (IER.EQ.0) GO TO 999 70 CONTINUE 80 CALL XERMSG ('SLATEC', 'DQNG', 'ABNORMAL RETURN', IER, 0) 999 RETURN END