1 subroutine qk15i(f,boun,inf,a,b,result,abserr,resabs,resasc,ierr) 2c***begin prologue qk15i 3c***date written 800101 (yymmdd) 4c***revision date 830518 (yymmdd) 5c***category no. h2a3a2,h2a4a2 6c***keywords 15-point transformed 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 the original (infinite integration range is mapped 10c onto the interval (0,1) and (a,b) is a part of (0,1). 11c it is the purpose to compute 12c i = integral of transformed integrand over (a,b), 13c j = integral of abs(transformed integrand) over (a,b). 14c***description 15c 16c integration rule 17c standard fortran subroutine 18c real version 19c 20c parameters 21c on entry 22c f - subroutine f(x,ierr,result) defining the integrand 23c function f(x). the actual name for f needs to be 24c declared e x t e r n a l in the calling program. 25c 26c boun - real 27c finite bound of original integration 28c range (set to zero if inf = +2) 29c 30c inf - integer 31c if inf = -1, the original interval is 32c (-infinity,bound), 33c if inf = +1, the original interval is 34c (bound,+infinity), 35c if inf = +2, the original interval is 36c (-infinity,+infinity) and 37c the integral is computed as the sum of two 38c integrals, one over (-infinity,0) and one over 39c (0,+infinity). 40c 41c a - real 42c lower limit for integration over subrange 43c of (0,1) 44c 45c b - real 46c upper limit for integration over subrange 47c of (0,1) 48c 49c on return 50c result - real 51c approximation to the integral i 52c result is computed by applying the 15-point 53c kronrod rule(resk) obtained by optimal addition 54c of abscissae to the 7-point gauss rule(resg). 55c 56c abserr - real 57c estimate of the modulus of the absolute error, 58c which should equal or exceed abs(i-result) 59c 60c resabs - real 61c approximation to the integral j 62c 63c resasc - real 64c approximation to the integral of 65c abs((transformed integrand)-i/(b-a)) over (a,b) 66c 67c***references (none) 68c***routines called r1mach 69c***end prologue qk15i 70c 71 real a,absc,absc1,absc2,abserr,b,boun,centr, 72 * dinf,r1mach,epmach,fc,fsum,fval1,fval2,fvalt,fv1, 73 * fv2,hlgth,resabs,resasc,resg,resk,reskh,result,tabsc1,tabsc2, 74 * uflow,wg,wgk,xgk 75 integer inf,j,min0 76 external f 77c 78 dimension fv1(7),fv2(7),xgk(8),wgk(8),wg(8) 79c 80c the abscissae and weights are supplied for the interval 81c (-1,1). because of symmetry only the positive abscissae and 82c their corresponding weights are given. 83c 84c xgk - abscissae of the 15-point kronrod rule 85c xgk(2), xgk(4), ... abscissae of the 7-point 86c gauss rule 87c xgk(1), xgk(3), ... abscissae which are optimally 88c added to the 7-point gauss rule 89c 90c wgk - weights of the 15-point kronrod rule 91c 92c wg - weights of the 7-point gauss rule, corresponding 93c to the abscissae xgk(2), xgk(4), ... 94c wg(1), wg(3), ... are set to zero. 95c 96 data xgk(1),xgk(2),xgk(3),xgk(4),xgk(5),xgk(6),xgk(7), 97 * xgk(8)/ 98 * 0.9914553711208126e+00, 0.9491079123427585e+00, 99 * 0.8648644233597691e+00, 0.7415311855993944e+00, 100 * 0.5860872354676911e+00, 0.4058451513773972e+00, 101 * 0.2077849550078985e+00, 0.0000000000000000e+00/ 102c 103 data wgk(1),wgk(2),wgk(3),wgk(4),wgk(5),wgk(6),wgk(7), 104 * wgk(8)/ 105 * 0.2293532201052922e-01, 0.6309209262997855e-01, 106 * 0.1047900103222502e+00, 0.1406532597155259e+00, 107 * 0.1690047266392679e+00, 0.1903505780647854e+00, 108 * 0.2044329400752989e+00, 0.2094821410847278e+00/ 109c 110 data wg(1),wg(2),wg(3),wg(4),wg(5),wg(6),wg(7),wg(8)/ 111 * 0.0000000000000000e+00, 0.1294849661688697e+00, 112 * 0.0000000000000000e+00, 0.2797053914892767e+00, 113 * 0.0000000000000000e+00, 0.3818300505051189e+00, 114 * 0.0000000000000000e+00, 0.4179591836734694e+00/ 115c 116c 117c list of major variables 118c ----------------------- 119c 120c centr - mid point of the interval 121c hlgth - half-length of the interval 122c absc* - abscissa 123c tabsc* - transformed abscissa 124c fval* - function value 125c resg - result of the 7-point gauss formula 126c resk - result of the 15-point kronrod formula 127c reskh - approximation to the mean value of the transformed 128c integrand over (a,b), i.e. to i/(b-a) 129c 130c machine dependent constants 131c --------------------------- 132c 133c epmach is the largest relative spacing. 134c uflow is the smallest positive magnitude. 135c 136c***first executable statement qk15i 137 epmach = r1mach(4) 138 uflow = r1mach(1) 139 dinf = min0(1,inf) 140c 141 centr = 0.5e+00*(a+b) 142 hlgth = 0.5e+00*(b-a) 143 tabsc1 = boun+dinf*(0.1e+01-centr)/centr 144 call f(tabsc1, ierr, fval1) 145 if (ierr.lt.0) return 146 if(inf.eq.2) then 147 call f(-tabsc1, ierr, fval1) 148 if (ierr.lt.0) return 149 fval1 = fval1 + fvalt 150 endif 151 fc = (fval1/centr)/centr 152c 153c compute the 15-point kronrod approximation to 154c the integral, and estimate the error. 155c 156 resg = wg(8)*fc 157 resk = wgk(8)*fc 158 resabs = abs(resk) 159 do 10 j=1,7 160 absc = hlgth*xgk(j) 161 absc1 = centr-absc 162 absc2 = centr+absc 163 tabsc1 = boun+dinf*(0.1e+01-absc1)/absc1 164 tabsc2 = boun+dinf*(0.1e+01-absc2)/absc2 165 call f(tabsc1, ierr, fval1) 166 if (ierr.lt.0) return 167 call f(tabsc2, ierr, fval2) 168 if (ierr.lt.0) return 169 if(inf.eq.2) then 170 call f(-tabsc1,ierr,fvalt) 171 if (ierr.lt.0) return 172 fval1 = fval1 + fvalt 173 endif 174 if(inf.eq.2) then 175 call f(-tabsc2,ierr,fvalt) 176 if (ierr.lt.0) return 177 fval2 = fval2 + fvalt 178 endif 179 fval1 = (fval1/absc1)/absc1 180 fval2 = (fval2/absc2)/absc2 181 fv1(j) = fval1 182 fv2(j) = fval2 183 fsum = fval1+fval2 184 resg = resg+wg(j)*fsum 185 resk = resk+wgk(j)*fsum 186 resabs = resabs+wgk(j)*(abs(fval1)+abs(fval2)) 187 10 continue 188 reskh = resk*0.5e+00 189 resasc = wgk(8)*abs(fc-reskh) 190 do 20 j=1,7 191 resasc = resasc+wgk(j)*(abs(fv1(j)-reskh)+abs(fv2(j)-reskh)) 192 20 continue 193 result = resk*hlgth 194 resasc = resasc*hlgth 195 resabs = resabs*hlgth 196 abserr = abs((resk-resg)*hlgth) 197 if(resasc.ne.0.0e+00.and.abserr.ne.0.e0) abserr = resasc* 198 * amin1(0.1e+01,(0.2e+03*abserr/resasc)**1.5e+00) 199 if(resabs.gt.uflow/(0.5e+02*epmach)) abserr = amax1 200 * ((epmach*0.5e+02)*resabs,abserr) 201 return 202 end 203