1 SUBROUTINE DQK15I(F,BOUN,INF,A,B,RESULT,ABSERR,RESABS,RESASC) 2C***BEGIN PROLOGUE DQK15I 3C***DATE WRITTEN 800101 (YYMMDD) 4C***REVISION DATE 830518 (YYMMDD) 5C***REVISION HISTORY (YYMMDD) 6C 000601 Changed DMAX1/DMIN1/DABS to generic MAX/MIN/ABS 7C***CATEGORY NO. H2A3A2,H2A4A2 8C***KEYWORDS 15-POINT TRANSFORMED GAUSS-KRONROD RULES 9C***AUTHOR PIESSENS, ROBERT, APPLIED MATH. AND PROGR. DIV. - 10C K. U. LEUVEN 11C DE DONCKER, ELISE, APPLIED MATH. AND PROGR. DIV. - 12C K. U. LEUVEN 13C***PURPOSE The original (infinite integration range is mapped 14C onto the interval (0,1) and (A,B) is a part of (0,1). 15C it is the purpose to compute 16C I = Integral of transformed integrand over (A,B), 17C J = Integral of ABS(Transformed Integrand) over (A,B). 18C***DESCRIPTION 19C 20C Integration Rule 21C Standard Fortran subroutine 22C Double precision version 23C 24C PARAMETERS 25C ON ENTRY 26C F - Double precision 27C Fuction subprogram defining the integrand 28C FUNCTION F(X). The actual name for F needs to be 29C Declared E X T E R N A L in the calling program. 30C 31C BOUN - Double precision 32C Finite bound of original integration 33C Range (SET TO ZERO IF INF = +2) 34C 35C INF - Integer 36C If INF = -1, the original interval is 37C (-INFINITY,BOUND), 38C If INF = +1, the original interval is 39C (BOUND,+INFINITY), 40C If INF = +2, the original interval is 41C (-INFINITY,+INFINITY) AND 42C The integral is computed as the sum of two 43C integrals, one over (-INFINITY,0) and one over 44C (0,+INFINITY). 45C 46C A - Double precision 47C Lower limit for integration over subrange 48C of (0,1) 49C 50C B - Double precision 51C Upper limit for integration over subrange 52C of (0,1) 53C 54C ON RETURN 55C RESULT - Double precision 56C Approximation to the integral I 57C Result is computed by applying the 15-POINT 58C KRONROD RULE(RESK) obtained by optimal addition 59C of abscissae to the 7-POINT GAUSS RULE(RESG). 60C 61C ABSERR - Double precision 62C Estimate of the modulus of the absolute error, 63C WHICH SHOULD EQUAL or EXCEED ABS(I-RESULT) 64C 65C RESABS - Double precision 66C Approximation to the integral J 67C 68C RESASC - Double precision 69C Approximation to the integral of 70C ABS((TRANSFORMED INTEGRAND)-I/(B-A)) over (A,B) 71C***REFERENCES (NONE) 72C***ROUTINES CALLED D1MACH 73C***END PROLOGUE DQK15I 74C 75 DOUBLE PRECISION A,ABSC,ABSC1,ABSC2,ABSERR,B,BOUN,CENTR,DINF, 76 1 D1MACH,EPMACH,F,FC,FSUM,FVAL1,FVAL2,FV1,FV2,HLGTH, 77 2 RESABS,RESASC,RESG,RESK,RESKH,RESULT,TABSC1,TABSC2,UFLOW,WG,WGK, 78 3 XGK 79 INTEGER INF,J 80 EXTERNAL F 81C 82 DIMENSION FV1(7),FV2(7),XGK(8),WGK(8),WG(8) 83C 84C THE ABSCISSAE AND WEIGHTS ARE SUPPLIED FOR THE INTERVAL 85C (-1,1). BECAUSE OF SYMMETRY ONLY THE POSITIVE ABSCISSAE AND 86C THEIR CORRESPONDING WEIGHTS ARE GIVEN. 87C 88C XGK - ABSCISSAE OF THE 15-POINT KRONROD RULE 89C XGK(2), XGK(4), ... ABSCISSAE OF THE 7-POINT 90C GAUSS RULE 91C XGK(1), XGK(3), ... ABSCISSAE WHICH ARE OPTIMALLY 92C ADDED TO THE 7-POINT GAUSS RULE 93C 94C WGK - WEIGHTS OF THE 15-POINT KRONROD RULE 95C 96C WG - WEIGHTS OF THE 7-POINT GAUSS RULE, CORRESPONDING 97C TO THE ABSCISSAE XGK(2), XGK(4), ... 98C WG(1), WG(3), ... ARE SET TO ZERO. 99C 100 DATA XGK(1),XGK(2),XGK(3),XGK(4),XGK(5),XGK(6),XGK(7),XGK(8)/ 101 1 0.9914553711208126D+00, 0.9491079123427585D+00, 102 2 0.8648644233597691D+00, 0.7415311855993944D+00, 103 3 0.5860872354676911D+00, 0.4058451513773972D+00, 104 4 0.2077849550078985D+00, 0.0000000000000000D+00/ 105C 106 DATA WGK(1),WGK(2),WGK(3),WGK(4),WGK(5),WGK(6),WGK(7),WGK(8)/ 107 1 0.2293532201052922D-01, 0.6309209262997855D-01, 108 2 0.1047900103222502D+00, 0.1406532597155259D+00, 109 3 0.1690047266392679D+00, 0.1903505780647854D+00, 110 4 0.2044329400752989D+00, 0.2094821410847278D+00/ 111C 112 DATA WG(1),WG(2),WG(3),WG(4),WG(5),WG(6),WG(7),WG(8)/ 113 1 0.0000000000000000D+00, 0.1294849661688697D+00, 114 2 0.0000000000000000D+00, 0.2797053914892767D+00, 115 3 0.0000000000000000D+00, 0.3818300505051189D+00, 116 4 0.0000000000000000D+00, 0.4179591836734694D+00/ 117C 118C 119C LIST OF MAJOR VARIABLES 120C ----------------------- 121C 122C CENTR - MID POINT OF THE INTERVAL 123C HLGTH - HALF-LENGTH OF THE INTERVAL 124C ABSC* - ABSCISSA 125C TABSC* - TRANSFORMED ABSCISSA 126C FVAL* - FUNCTION VALUE 127C RESG - RESULT OF THE 7-POINT GAUSS FORMULA 128C RESK - RESULT OF THE 15-POINT KRONROD FORMULA 129C RESKH - APPROXIMATION TO THE MEAN VALUE OF THE TRANSFORMED 130C INTEGRAND OVER (A,B), I.E. TO I/(B-A) 131C 132C MACHINE DEPENDENT CONSTANTS 133C --------------------------- 134C 135C EPMACH IS THE LARGEST RELATIVE SPACING. 136C UFLOW IS THE SMALLEST POSITIVE MAGNITUDE. 137C 138C***FIRST EXECUTABLE STATEMENT DQK15I 139 EPMACH = D1MACH(4) 140 UFLOW = D1MACH(1) 141 DINF = MIN0(1,INF) 142C 143 CENTR = 0.5D+00*(A+B) 144 HLGTH = 0.5D+00*(B-A) 145 TABSC1 = BOUN+DINF*(0.1D+01-CENTR)/CENTR 146 FVAL1 = F(TABSC1) 147 IF(INF.EQ.2) FVAL1 = FVAL1+F(-TABSC1) 148 FC = (FVAL1/CENTR)/CENTR 149C 150C COMPUTE THE 15-POINT KRONROD APPROXIMATION TO 151C THE INTEGRAL, AND ESTIMATE THE ERROR. 152C 153 RESG = WG(8)*FC 154 RESK = WGK(8)*FC 155 RESABS = ABS(RESK) 156 DO 10 J=1,7 157 ABSC = HLGTH*XGK(J) 158 ABSC1 = CENTR-ABSC 159 ABSC2 = CENTR+ABSC 160 TABSC1 = BOUN+DINF*(0.1D+01-ABSC1)/ABSC1 161 TABSC2 = BOUN+DINF*(0.1D+01-ABSC2)/ABSC2 162 FVAL1 = F(TABSC1) 163 FVAL2 = F(TABSC2) 164 IF(INF.EQ.2) FVAL1 = FVAL1+F(-TABSC1) 165 IF(INF.EQ.2) FVAL2 = FVAL2+F(-TABSC2) 166 FVAL1 = (FVAL1/ABSC1)/ABSC1 167 FVAL2 = (FVAL2/ABSC2)/ABSC2 168 FV1(J) = FVAL1 169 FV2(J) = FVAL2 170 FSUM = FVAL1+FVAL2 171 RESG = RESG+WG(J)*FSUM 172 RESK = RESK+WGK(J)*FSUM 173 RESABS = RESABS+WGK(J)*(ABS(FVAL1)+ABS(FVAL2)) 174 10 CONTINUE 175 RESKH = RESK*0.5D+00 176 RESASC = WGK(8)*ABS(FC-RESKH) 177 DO 20 J=1,7 178 RESASC = RESASC+WGK(J)*(ABS(FV1(J)-RESKH)+ABS(FV2(J)-RESKH)) 179 20 CONTINUE 180 RESULT = RESK*HLGTH 181 RESASC = RESASC*HLGTH 182 RESABS = RESABS*HLGTH 183 ABSERR = ABS((RESK-RESG)*HLGTH) 184 IF(RESASC.NE.0.0D+00.AND.ABSERR.NE.0.D0) ABSERR = RESASC* 185 1 MIN(0.1D+01,(0.2D+03*ABSERR/RESASC)**1.5D+00) 186 IF(RESABS.GT.UFLOW/(0.5D+02*EPMACH)) ABSERR = MAX 187 1 ((EPMACH*0.5D+02)*RESABS,ABSERR) 188 RETURN 189 END 190