1*DECK DQK21 2 SUBROUTINE DQK21 (F, A, B, RESULT, ABSERR, RESABS, RESASC) 3C***BEGIN PROLOGUE DQK21 4C***PURPOSE To compute I = Integral of F over (A,B), with error 5C estimate 6C J = Integral of ABS(F) over (A,B) 7C***LIBRARY SLATEC (QUADPACK) 8C***CATEGORY H2A1A2 9C***TYPE DOUBLE PRECISION (QK21-S, DQK21-D) 10C***KEYWORDS 21-POINT GAUSS-KRONROD RULES, QUADPACK, QUADRATURE 11C***AUTHOR Piessens, Robert 12C Applied Mathematics and Programming Division 13C K. U. Leuven 14C de Doncker, Elise 15C Applied Mathematics and Programming Division 16C K. U. Leuven 17C***DESCRIPTION 18C 19C Integration rules 20C Standard fortran subroutine 21C Double precision version 22C 23C PARAMETERS 24C ON ENTRY 25C F - Double precision 26C Function subprogram defining the integrand 27C FUNCTION F(X). The actual name for F needs to be 28C Declared E X T E R N A L in the driver program. 29C 30C A - Double precision 31C Lower limit of integration 32C 33C B - Double precision 34C Upper limit of integration 35C 36C ON RETURN 37C RESULT - Double precision 38C Approximation to the integral I 39C RESULT is computed by applying the 21-POINT 40C KRONROD RULE (RESK) obtained by optimal addition 41C of abscissae to the 10-POINT GAUSS RULE (RESG). 42C 43C ABSERR - Double precision 44C Estimate of the modulus of the absolute error, 45C which should not exceed ABS(I-RESULT) 46C 47C RESABS - Double precision 48C Approximation to the integral J 49C 50C RESASC - Double precision 51C Approximation to the integral of ABS(F-I/(B-A)) 52C over (A,B) 53C 54C***REFERENCES (NONE) 55C***ROUTINES CALLED D1MACH 56C***REVISION HISTORY (YYMMDD) 57C 800101 DATE WRITTEN 58C 890531 Changed all specific intrinsics to generic. (WRB) 59C 890531 REVISION DATE from Version 3.2 60C 891214 Prologue converted to Version 4.0 format. (BAB) 61C***END PROLOGUE DQK21 62C 63 DOUBLE PRECISION A,ABSC,ABSERR,B,CENTR,DHLGTH, 64 1 D1MACH,EPMACH,F,FC,FSUM,FVAL1,FVAL2,FV1,FV2,HLGTH,RESABS,RESASC, 65 2 RESG,RESK,RESKH,RESULT,UFLOW,WG,WGK,XGK 66 INTEGER J,JTW,JTWM1 67 EXTERNAL F 68C 69 DIMENSION FV1(10),FV2(10),WG(5),WGK(11),XGK(11) 70C 71C THE ABSCISSAE AND WEIGHTS ARE GIVEN FOR THE INTERVAL (-1,1). 72C BECAUSE OF SYMMETRY ONLY THE POSITIVE ABSCISSAE AND THEIR 73C CORRESPONDING WEIGHTS ARE GIVEN. 74C 75C XGK - ABSCISSAE OF THE 21-POINT KRONROD RULE 76C XGK(2), XGK(4), ... ABSCISSAE OF THE 10-POINT 77C GAUSS RULE 78C XGK(1), XGK(3), ... ABSCISSAE WHICH ARE OPTIMALLY 79C ADDED TO THE 10-POINT GAUSS RULE 80C 81C WGK - WEIGHTS OF THE 21-POINT KRONROD RULE 82C 83C WG - WEIGHTS OF THE 10-POINT GAUSS RULE 84C 85C 86C GAUSS QUADRATURE WEIGHTS AND KRONROD QUADRATURE ABSCISSAE AND WEIGHTS 87C AS EVALUATED WITH 80 DECIMAL DIGIT ARITHMETIC BY L. W. FULLERTON, 88C BELL LABS, NOV. 1981. 89C 90 SAVE WG, XGK, WGK 91 DATA WG ( 1) / 0.0666713443 0868813759 3568809893 332 D0 / 92 DATA WG ( 2) / 0.1494513491 5058059314 5776339657 697 D0 / 93 DATA WG ( 3) / 0.2190863625 1598204399 5534934228 163 D0 / 94 DATA WG ( 4) / 0.2692667193 0999635509 1226921569 469 D0 / 95 DATA WG ( 5) / 0.2955242247 1475287017 3892994651 338 D0 / 96C 97 DATA XGK ( 1) / 0.9956571630 2580808073 5527280689 003 D0 / 98 DATA XGK ( 2) / 0.9739065285 1717172007 7964012084 452 D0 / 99 DATA XGK ( 3) / 0.9301574913 5570822600 1207180059 508 D0 / 100 DATA XGK ( 4) / 0.8650633666 8898451073 2096688423 493 D0 / 101 DATA XGK ( 5) / 0.7808177265 8641689706 3717578345 042 D0 / 102 DATA XGK ( 6) / 0.6794095682 9902440623 4327365114 874 D0 / 103 DATA XGK ( 7) / 0.5627571346 6860468333 9000099272 694 D0 / 104 DATA XGK ( 8) / 0.4333953941 2924719079 9265943165 784 D0 / 105 DATA XGK ( 9) / 0.2943928627 0146019813 1126603103 866 D0 / 106 DATA XGK ( 10) / 0.1488743389 8163121088 4826001129 720 D0 / 107 DATA XGK ( 11) / 0.0000000000 0000000000 0000000000 000 D0 / 108C 109 DATA WGK ( 1) / 0.0116946388 6737187427 8064396062 192 D0 / 110 DATA WGK ( 2) / 0.0325581623 0796472747 8818972459 390 D0 / 111 DATA WGK ( 3) / 0.0547558965 7435199603 1381300244 580 D0 / 112 DATA WGK ( 4) / 0.0750396748 1091995276 7043140916 190 D0 / 113 DATA WGK ( 5) / 0.0931254545 8369760553 5065465083 366 D0 / 114 DATA WGK ( 6) / 0.1093871588 0229764189 9210590325 805 D0 / 115 DATA WGK ( 7) / 0.1234919762 6206585107 7958109831 074 D0 / 116 DATA WGK ( 8) / 0.1347092173 1147332592 8054001771 707 D0 / 117 DATA WGK ( 9) / 0.1427759385 7706008079 7094273138 717 D0 / 118 DATA WGK ( 10) / 0.1477391049 0133849137 4841515972 068 D0 / 119 DATA WGK ( 11) / 0.1494455540 0291690566 4936468389 821 D0 / 120C 121C 122C LIST OF MAJOR VARIABLES 123C ----------------------- 124C 125C CENTR - MID POINT OF THE INTERVAL 126C HLGTH - HALF-LENGTH OF THE INTERVAL 127C ABSC - ABSCISSA 128C FVAL* - FUNCTION VALUE 129C RESG - RESULT OF THE 10-POINT GAUSS FORMULA 130C RESK - RESULT OF THE 21-POINT KRONROD FORMULA 131C RESKH - APPROXIMATION TO THE MEAN VALUE OF F OVER (A,B), 132C I.E. TO I/(B-A) 133C 134C 135C MACHINE DEPENDENT CONSTANTS 136C --------------------------- 137C 138C EPMACH IS THE LARGEST RELATIVE SPACING. 139C UFLOW IS THE SMALLEST POSITIVE MAGNITUDE. 140C 141C***FIRST EXECUTABLE STATEMENT DQK21 142 EPMACH = D1MACH(4) 143 UFLOW = D1MACH(1) 144C 145 CENTR = 0.5D+00*(A+B) 146 HLGTH = 0.5D+00*(B-A) 147 DHLGTH = ABS(HLGTH) 148C 149C COMPUTE THE 21-POINT KRONROD APPROXIMATION TO 150C THE INTEGRAL, AND ESTIMATE THE ABSOLUTE ERROR. 151C 152 RESG = 0.0D+00 153 FC = F(CENTR) 154 RESK = WGK(11)*FC 155 RESABS = ABS(RESK) 156 DO 10 J=1,5 157 JTW = 2*J 158 ABSC = HLGTH*XGK(JTW) 159 FVAL1 = F(CENTR-ABSC) 160 FVAL2 = F(CENTR+ABSC) 161 FV1(JTW) = FVAL1 162 FV2(JTW) = FVAL2 163 FSUM = FVAL1+FVAL2 164 RESG = RESG+WG(J)*FSUM 165 RESK = RESK+WGK(JTW)*FSUM 166 RESABS = RESABS+WGK(JTW)*(ABS(FVAL1)+ABS(FVAL2)) 167 10 CONTINUE 168 DO 15 J = 1,5 169 JTWM1 = 2*J-1 170 ABSC = HLGTH*XGK(JTWM1) 171 FVAL1 = F(CENTR-ABSC) 172 FVAL2 = F(CENTR+ABSC) 173 FV1(JTWM1) = FVAL1 174 FV2(JTWM1) = FVAL2 175 FSUM = FVAL1+FVAL2 176 RESK = RESK+WGK(JTWM1)*FSUM 177 RESABS = RESABS+WGK(JTWM1)*(ABS(FVAL1)+ABS(FVAL2)) 178 15 CONTINUE 179 RESKH = RESK*0.5D+00 180 RESASC = WGK(11)*ABS(FC-RESKH) 181 DO 20 J=1,10 182 RESASC = RESASC+WGK(J)*(ABS(FV1(J)-RESKH)+ABS(FV2(J)-RESKH)) 183 20 CONTINUE 184 RESULT = RESK*HLGTH 185 RESABS = RESABS*DHLGTH 186 RESASC = RESASC*DHLGTH 187 ABSERR = ABS((RESK-RESG)*HLGTH) 188 IF(RESASC.NE.0.0D+00.AND.ABSERR.NE.0.0D+00) 189 1 ABSERR = RESASC*MIN(0.1D+01,(0.2D+03*ABSERR/RESASC)**1.5D+00) 190 IF(RESABS.GT.UFLOW/(0.5D+02*EPMACH)) ABSERR = MAX 191 1 ((EPMACH*0.5D+02)*RESABS,ABSERR) 192 RETURN 193 END 194