subroutine qk15w(f,w,p1,p2,p3,p4,kp,a,b,result,abserr,
* resabs,resasc)
c***begin prologue qk15w
c***date written 810101 (yymmdd)
c***revision date 830518 (mmddyy)
c***category no. h2a2a2
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*w over (a,b), with error
c estimate
c j = integral of abs(f*w) over (a,b)
c***description
c
c integration rules
c standard fortran subroutine
c real version
c
c parameters
c on entry
c f - real
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 driver program.
c
c w - real
c function subprogram defining the integrand
c weight function w(x). the actual name for w
c needs to be declared e x t e r n a l in the
c calling program.
c
c p1, p2, p3, p4 - real
c parameters in the weight function
c
c kp - integer
c key for indicating the type of weight function
c
c a - real
c lower limit of integration
c
c b - real
c upper limit of integration
c
c on return
c result - real
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 the 7-point gauss rule (resg).
c
c abserr - real
c estimate of the modulus of the absolute error,
c which should equal or exceed abs(i-result)
c
c resabs - real
c approximation to the integral of abs(f)
c
c resasc - real
c approximation to the integral of abs(f-i/(b-a))
c
c***references (none)
c***routines called r1mach
c***end prologue qk15w
c
real a,absc,absc1,absc2,abserr,b,centr,dhlgth,
* r1mach,epmach,f,fc,fsum,fval1,fval2,fv1,fv2,
* hlgth,p1,p2,p3,p4,resabs,resasc,resg,resk,reskh,result,uflow,
* w,wg,wgk,xgk
integer j,jtw,jtwm1,kp
external f,w
c
dimension fv1(7),fv2(7),xgk(8),wgk(8),wg(4)
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 gauss-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 gauss-kronrod rule
c
c wg - weights of the 7-point gauss rule
c
data xgk(1),xgk(2),xgk(3),xgk(4),xgk(5),xgk(6),xgk(7),
* xgk(8)/
* 0.9914553711208126e+00, 0.9491079123427585e+00,
* 0.8648644233597691e+00, 0.7415311855993944e+00,
* 0.5860872354676911e+00, 0.4058451513773972e+00,
* 0.2077849550078985e+00, 0.0000000000000000e+00/
c
data wgk(1),wgk(2),wgk(3),wgk(4),wgk(5),wgk(6),wgk(7),
* wgk(8)/
* 0.2293532201052922e-01, 0.6309209262997855e-01,
* 0.1047900103222502e+00, 0.1406532597155259e+00,
* 0.1690047266392679e+00, 0.1903505780647854e+00,
* 0.2044329400752989e+00, 0.2094821410847278e+00/
c
data wg(1),wg(2),wg(3),wg(4)/
* 0.1294849661688697e+00, 0.2797053914892767e+00,
* 0.3818300505051189e+00, 0.4179591836734694e+00/
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*w 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 qk15w
epmach = r1mach(4)
uflow = r1mach(1)
c
centr = 0.5e+00*(a+b)
hlgth = 0.5e+00*(b-a)
dhlgth = abs(hlgth)
c
c compute the 15-point kronrod approximation to the
c integral, and estimate the error.
c
fc = f(centr)*w(centr,p1,p2,p3,p4,kp)
resg = wg(4)*fc
resk = wgk(8)*fc
resabs = abs(resk)
do 10 j=1,3
jtw = j*2
absc = hlgth*xgk(jtw)
absc1 = centr-absc
absc2 = centr+absc
fval1 = f(absc1)*w(absc1,p1,p2,p3,p4,kp)
fval2 = f(absc2)*w(absc2,p1,p2,p3,p4,kp)
fv1(jtw) = fval1
fv2(jtw) = fval2
fsum = fval1+fval2
resg = resg+wg(j)*fsum
resk = resk+wgk(jtw)*fsum
resabs = resabs+wgk(jtw)*(abs(fval1)+abs(fval2))
10 continue
do 15 j=1,4
jtwm1 = j*2-1
absc = hlgth*xgk(jtwm1)
absc1 = centr-absc
absc2 = centr+absc
fval1 = f(absc1)*w(absc1,p1,p2,p3,p4,kp)
fval2 = f(absc2)*w(absc2,p1,p2,p3,p4,kp)
fv1(jtwm1) = fval1
fv2(jtwm1) = fval2
fsum = fval1+fval2
resk = resk+wgk(jtwm1)*fsum
resabs = resabs+wgk(jtwm1)*(abs(fval1)+abs(fval2))
15 continue
reskh = resk*0.5e+00
resasc = wgk(8)*abs(fc-reskh)
do 20 j=1,7
resasc = resasc+wgk(j)*(abs(fv1(j)-reskh)+abs(fv2(j)-reskh))
20 continue
result = resk*hlgth
resabs = resabs*dhlgth
resasc = resasc*dhlgth
abserr = abs((resk-resg)*hlgth)
if(resasc.ne.0.0e+00.and.abserr.ne.0.0e+00)
* abserr = resasc*amin1(0.1e+01,
* (0.2e+03*abserr/resasc)**1.5e+00)
if(resabs.gt.uflow/(0.5e+02*epmach)) abserr = amax1((epmach*
* 0.5e+02)*resabs,abserr)
return
end
.