1 SUBROUTINE DQK21(F,A,B,RESULT,ABSERR,RESABS,RESASC,IERR) 2C***BEGIN PROLOGUE DQK21 3C***DATE WRITTEN 800101 (YYMMDD) 4C***REVISION DATE 830518 (YYMMDD) 5C***CATEGORY NO. H2A1A2 6C***KEYWORDS 21-POINT GAUSS-KRONROD RULES 7C***AUTHOR PIESSENS,ROBERT,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN 8C DE DONCKER,ELISE,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN 9C***PURPOSE TO COMPUTE I = INTEGRAL OF F OVER (A,B), WITH ERROR 10C ESTIMATE 11C J = INTEGRAL OF ABS(F) OVER (A,B) 12C***DESCRIPTION 13C 14C INTEGRATION RULES 15C STANDARD FORTRAN SUBROUTINE 16C DOUBLE PRECISION VERSION 17C 18C PARAMETERS 19C ON ENTRY 20C F - SUBROUTINE F(X,IERR,RESULT) DEFINING THE INTEGRAND 21C FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE 22C DECLARED E X T E R N A L IN THE DRIVER PROGRAM. 23C 24C A - DOUBLE PRECISION 25C LOWER LIMIT OF INTEGRATION 26C 27C B - DOUBLE PRECISION 28C UPPER LIMIT OF INTEGRATION 29C 30C ON RETURN 31C RESULT - DOUBLE PRECISION 32C APPROXIMATION TO THE INTEGRAL I 33C RESULT IS COMPUTED BY APPLYING THE 21-POINT 34C KRONROD RULE (RESK) OBTAINED BY OPTIMAL ADDITION 35C OF ABSCISSAE TO THE 10-POINT GAUSS RULE (RESG). 36C 37C ABSERR - DOUBLE PRECISION 38C ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR, 39C WHICH SHOULD NOT EXCEED ABS(I-RESULT) 40C 41C RESABS - DOUBLE PRECISION 42C APPROXIMATION TO THE INTEGRAL J 43C 44C RESASC - DOUBLE PRECISION 45C APPROXIMATION TO THE INTEGRAL OF ABS(F-I/(B-A)) 46C OVER (A,B) 47C 48C***REFERENCES (NONE) 49C***ROUTINES CALLED D1MACH 50C***END PROLOGUE DQK21 51C 52 DOUBLE PRECISION A,ABSC,ABSERR,B,CENTR,DABS,DHLGTH,DMAX1,DMIN1, 53 * D1MACH,EPMACH,FC,FSUM,FVAL1,FVAL2,FV1,FV2,HLGTH,RESABS,RESASC, 54 * RESG,RESK,RESKH,RESULT,UFLOW,WG,WGK,XGK 55 INTEGER J,JTW,JTWM1 56 EXTERNAL F 57C 58 DIMENSION FV1(10),FV2(10),WG(5),WGK(11),XGK(11) 59C 60C THE ABSCISSAE AND WEIGHTS ARE GIVEN FOR THE INTERVAL (-1,1). 61C BECAUSE OF SYMMETRY ONLY THE POSITIVE ABSCISSAE AND THEIR 62C CORRESPONDING WEIGHTS ARE GIVEN. 63C 64C XGK - ABSCISSAE OF THE 21-POINT KRONROD RULE 65C XGK(2), XGK(4), ... ABSCISSAE OF THE 10-POINT 66C GAUSS RULE 67C XGK(1), XGK(3), ... ABSCISSAE WHICH ARE OPTIMALLY 68C ADDED TO THE 10-POINT GAUSS RULE 69C 70C WGK - WEIGHTS OF THE 21-POINT KRONROD RULE 71C 72C WG - WEIGHTS OF THE 10-POINT GAUSS RULE 73C 74C 75C GAUSS QUADRATURE WEIGHTS AND KRONRON QUADRATURE ABSCISSAE AND WEIGHTS 76C AS EVALUATED WITH 80 DECIMAL DIGIT ARITHMETIC BY L. W. FULLERTON, 77C BELL LABS, NOV. 1981. 78C 79 DATA WG ( 1) / 0.0666713443 0868813759 3568809893 332 D0 / 80 DATA WG ( 2) / 0.1494513491 5058059314 5776339657 697 D0 / 81 DATA WG ( 3) / 0.2190863625 1598204399 5534934228 163 D0 / 82 DATA WG ( 4) / 0.2692667193 0999635509 1226921569 469 D0 / 83 DATA WG ( 5) / 0.2955242247 1475287017 3892994651 338 D0 / 84C 85 DATA XGK ( 1) / 0.9956571630 2580808073 5527280689 003 D0 / 86 DATA XGK ( 2) / 0.9739065285 1717172007 7964012084 452 D0 / 87 DATA XGK ( 3) / 0.9301574913 5570822600 1207180059 508 D0 / 88 DATA XGK ( 4) / 0.8650633666 8898451073 2096688423 493 D0 / 89 DATA XGK ( 5) / 0.7808177265 8641689706 3717578345 042 D0 / 90 DATA XGK ( 6) / 0.6794095682 9902440623 4327365114 874 D0 / 91 DATA XGK ( 7) / 0.5627571346 6860468333 9000099272 694 D0 / 92 DATA XGK ( 8) / 0.4333953941 2924719079 9265943165 784 D0 / 93 DATA XGK ( 9) / 0.2943928627 0146019813 1126603103 866 D0 / 94 DATA XGK ( 10) / 0.1488743389 8163121088 4826001129 720 D0 / 95 DATA XGK ( 11) / 0.0000000000 0000000000 0000000000 000 D0 / 96C 97 DATA WGK ( 1) / 0.0116946388 6737187427 8064396062 192 D0 / 98 DATA WGK ( 2) / 0.0325581623 0796472747 8818972459 390 D0 / 99 DATA WGK ( 3) / 0.0547558965 7435199603 1381300244 580 D0 / 100 DATA WGK ( 4) / 0.0750396748 1091995276 7043140916 190 D0 / 101 DATA WGK ( 5) / 0.0931254545 8369760553 5065465083 366 D0 / 102 DATA WGK ( 6) / 0.1093871588 0229764189 9210590325 805 D0 / 103 DATA WGK ( 7) / 0.1234919762 6206585107 7958109831 074 D0 / 104 DATA WGK ( 8) / 0.1347092173 1147332592 8054001771 707 D0 / 105 DATA WGK ( 9) / 0.1427759385 7706008079 7094273138 717 D0 / 106 DATA WGK ( 10) / 0.1477391049 0133849137 4841515972 068 D0 / 107 DATA WGK ( 11) / 0.1494455540 0291690566 4936468389 821 D0 / 108C 109C 110C LIST OF MAJOR VARIABLES 111C ----------------------- 112C 113C CENTR - MID POINT OF THE INTERVAL 114C HLGTH - HALF-LENGTH OF THE INTERVAL 115C ABSC - ABSCISSA 116C FVAL* - FUNCTION VALUE 117C RESG - RESULT OF THE 10-POINT GAUSS FORMULA 118C RESK - RESULT OF THE 21-POINT KRONROD FORMULA 119C RESKH - APPROXIMATION TO THE MEAN VALUE OF F OVER (A,B), 120C I.E. TO I/(B-A) 121C 122C 123C MACHINE DEPENDENT CONSTANTS 124C --------------------------- 125C 126C EPMACH IS THE LARGEST RELATIVE SPACING. 127C UFLOW IS THE SMALLEST POSITIVE MAGNITUDE. 128C 129C***FIRST EXECUTABLE STATEMENT DQK21 130 EPMACH = D1MACH(4) 131 UFLOW = D1MACH(1) 132C 133 CENTR = 0.5D+00*(A+B) 134 HLGTH = 0.5D+00*(B-A) 135 DHLGTH = DABS(HLGTH) 136C 137C COMPUTE THE 21-POINT KRONROD APPROXIMATION TO 138C THE INTEGRAL, AND ESTIMATE THE ABSOLUTE ERROR. 139C 140 RESG = 0.0D+00 141 IERR = 0 142 CALL F(CENTR,IERR,FC) 143 IF (IERR .LT. 0) RETURN 144 RESK = WGK(11)*FC 145 RESABS = DABS(RESK) 146 DO 10 J=1,5 147 JTW = 2*J 148 ABSC = HLGTH*XGK(JTW) 149 CALL F(CENTR-ABSC,IERR,FVAL1) 150 IF (IERR .LT. 0) RETURN 151 CALL F(CENTR+ABSC,IERR,FVAL2) 152 IF (IERR .LT. 0) RETURN 153 FV1(JTW) = FVAL1 154 FV2(JTW) = FVAL2 155 FSUM = FVAL1+FVAL2 156 RESG = RESG+WG(J)*FSUM 157 RESK = RESK+WGK(JTW)*FSUM 158 RESABS = RESABS+WGK(JTW)*(DABS(FVAL1)+DABS(FVAL2)) 159 10 CONTINUE 160 DO 15 J = 1,5 161 JTWM1 = 2*J-1 162 ABSC = HLGTH*XGK(JTWM1) 163 CALL F(CENTR-ABSC,IERR,FVAL1) 164 IF (IERR .LT. 0) RETURN 165 CALL F(CENTR+ABSC,IERR,FVAL2) 166 IF (IERR .LT. 0) RETURN 167 FV1(JTWM1) = FVAL1 168 FV2(JTWM1) = FVAL2 169 FSUM = FVAL1+FVAL2 170 RESK = RESK+WGK(JTWM1)*FSUM 171 RESABS = RESABS+WGK(JTWM1)*(DABS(FVAL1)+DABS(FVAL2)) 172 15 CONTINUE 173 RESKH = RESK*0.5D+00 174 RESASC = WGK(11)*DABS(FC-RESKH) 175 DO 20 J=1,10 176 RESASC = RESASC+WGK(J)*(DABS(FV1(J)-RESKH)+DABS(FV2(J)-RESKH)) 177 20 CONTINUE 178 RESULT = RESK*HLGTH 179 RESABS = RESABS*DHLGTH 180 RESASC = RESASC*DHLGTH 181 ABSERR = DABS((RESK-RESG)*HLGTH) 182 IF(RESASC.NE.0.0D+00.AND.ABSERR.NE.0.0D+00) 183 * ABSERR = RESASC*DMIN1(0.1D+01,(0.2D+03*ABSERR/RESASC)**1.5D+00) 184 IF(RESABS.GT.UFLOW/(0.5D+02*EPMACH)) ABSERR = DMAX1 185 * ((EPMACH*0.5D+02)*RESABS,ABSERR) 186 RETURN 187 END 188