``````       subroutine dqk15r(f,a,b,result,abserr,resabs,resasc,i)
c***begin prologue  dqk15
c***date written   800101   (yymmdd)
c***revision date  830518   (yymmdd)
c***category no.  h2a1a2
c***keywords  15-point gauss-kronrod rules
c***author  piessens,robert,appl. math. & progr. div. - k.u.leuven
c           de doncker,elise,appl. math. & progr. div - k.u.leuven
c***purpose  to compute i = integral of f over (a,b), with error
c                           estimate
c                       j = integral of abs(f) over (a,b)
c***description
c
c           integration rules
c           standard fortran subroutine
c           double precision version
c
c           parameters
c            on entry
c              f      - double precision
c                       function subprogram defining the integrand
c                       function f(x). the actual name for f needs to be
c                       declared e x t e r n a l in the calling program.
c
c              a      - double precision
c                       lower limit of integration
c
c              b      - double precision
c                       upper limit of integration
c
c            on return
c              result - double precision
c                       approximation to the integral i
c                       result is computed by applying the 15-point
c                       kronrod rule (resk) obtained by optimal addition
c                       of abscissae to the7-point gauss rule(resg).
c
c              abserr - double precision
c                       estimate of the modulus of the absolute error,
c                       which should not exceed abs(i-result)
c
c              resabs - double precision
c                       approximation to the integral j
c
c              resasc - double precision
c                       approximation to the integral of abs(f-i/(b-a))
c                       over (a,b)
c
c***references  (none)
c***routines called  d1mach
c***end prologue  dqk15
c
double precision a,absc,abserr,b,centr,dabs,dhlgth,dmax1,dmin1,
+d1mach,epmach,f,fc,fsum,fval1,fval2,fv1,fv2,hlgth,resabs,resasc,
+resg,resk,reskh,result,uflow,wg,wgk,xgk
integer j,jtw,jtwm1,i
external f
c
dimension fv1(7),fv2(7),wg(4),wgk(8),xgk(8)
c
c           the abscissae and weights are given for the interval (-1,1).
c           because of symmetry only the positive abscissae and their
c           corresponding weights are given.
c
c           xgk    - abscissae of the 15-point kronrod rule
c                    xgk(2), xgk(4), ...  abscissae of the 7-point
c                    gauss rule
c                    xgk(1), xgk(3), ...  abscissae which are optimally
c                    added to the 7-point gauss rule
c
c           wgk    - weights of the 15-point kronrod rule
c
c           wg     - weights of the 7-point gauss rule
c
c
c as evaluated with 80 decimal digit arithmetic by l. w. fullerton,
c bell labs, nov. 1981.
c
data wg  (  1) / 0.1294849661 6886969327 0611432679 082 d0 /
data wg  (  2) / 0.2797053914 8927666790 1467771423 780 d0 /
data wg  (  3) / 0.3818300505 0511894495 0369775488 975 d0 /
data wg  (  4) / 0.4179591836 7346938775 5102040816 327 d0 /
c
data xgk (  1) / 0.9914553711 2081263920 6854697526 329 d0 /
data xgk (  2) / 0.9491079123 4275852452 6189684047 851 d0 /
data xgk (  3) / 0.8648644233 5976907278 9712788640 926 d0 /
data xgk (  4) / 0.7415311855 9939443986 3864773280 788 d0 /
data xgk (  5) / 0.5860872354 6769113029 4144838258 730 d0 /
data xgk (  6) / 0.4058451513 7739716690 6606412076 961 d0 /
data xgk (  7) / 0.2077849550 0789846760 0689403773 245 d0 /
data xgk (  8) / 0.0000000000 0000000000 0000000000 000 d0 /
c
data wgk (  1) / 0.0229353220 1052922496 3732008058 970 d0 /
data wgk (  2) / 0.0630920926 2997855329 0700663189 204 d0 /
data wgk (  3) / 0.1047900103 2225018383 9876322541 518 d0 /
data wgk (  4) / 0.1406532597 1552591874 5189590510 238 d0 /
data wgk (  5) / 0.1690047266 3926790282 6583426598 550 d0 /
data wgk (  6) / 0.1903505780 6478540991 3256402421 014 d0 /
data wgk (  7) / 0.2044329400 7529889241 4161999234 649 d0 /
data wgk (  8) / 0.2094821410 8472782801 2999174891 714 d0 /
c
c
c           list of major variables
c           -----------------------
c
c           centr  - mid point of the interval
c           hlgth  - half-length of the interval
c           absc   - abscissa
c           fval*  - function value
c           resg   - result of the 7-point gauss formula
c           resk   - result of the 15-point kronrod formula
c           reskh  - approximation to the mean value of f over (a,b),
c                    i.e. to 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  dqk15
epmach = d1mach(4)
uflow = d1mach(1)
c
centr = 0.5d+00*(a+b)
hlgth = 0.5d+00*(b-a)
dhlgth = dabs(hlgth)
c
c           compute the 15-point kronrod approximation to
c           the integral, and estimate the absolute error.
c
fc = f(centr,i)
resg = fc*wg(4)
resk = fc*wgk(8)
resabs = dabs(resk)
do 10 j=1,3
jtw = j*2
absc = hlgth*xgk(jtw)
fval1 = f(centr-absc,i)
fval2 = f(centr+absc,i)
fv1(jtw) = fval1
fv2(jtw) = fval2
fsum = fval1+fval2
resg = resg+wg(j)*fsum
resk = resk+wgk(jtw)*fsum
resabs = resabs+wgk(jtw)*(dabs(fval1)+dabs(fval2))
10 continue
do 15 j = 1,4
jtwm1 = j*2-1
absc = hlgth*xgk(jtwm1)
fval1 = f(centr-absc,i)
fval2 = f(centr+absc,i)
fv1(jtwm1) = fval1
fv2(jtwm1) = fval2
fsum = fval1+fval2
resk = resk+wgk(jtwm1)*fsum
resabs = resabs+wgk(jtwm1)*(dabs(fval1)+dabs(fval2))
15 continue
reskh = resk*0.5d+00
resasc = wgk(8)*dabs(fc-reskh)
do 20 j=1,7
resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh))
20 continue
result = resk*hlgth
resabs = resabs*dhlgth
resasc = resasc*dhlgth
abserr = dabs((resk-resg)*hlgth)
if(resasc.ne.0.0d+00.and.abserr.ne.0.0d+00)
+abserr = resasc*dmin1(0.1d+01,(0.2d+03*abserr/resasc)**1.5d+00)
if(resabs.gt.uflow/(0.5d+02*epmach)) abserr = dmax1
+((epmach*0.5d+02)*resabs,abserr)
return
end
``````