subroutine qag(f,a,b,epsabs,epsrel,key,result,abserr,neval,ier, * limit,lenw,last,iwork,work) c***begin prologue qag c***date written 800101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a1a1 c***keywords automatic integrator, general-purpose, c integrand examinator, globally adaptive, c gauss-kronrod c***author piessens,robert,appl. math. & progr. div - k.u.leuven c de doncker,elise,appl. math. & progr. div. - k.u.leuven c***purpose the routine calculates an approximation result to a given c definite integral i = integral of f over (a,b), c hopefully satisfying following claim for accuracy c abs(i-result)le.max(epsabs,epsrel*abs(i)). c***description c c computation of a definite integral c standard fortran subroutine c real version c c f - real c function subprogam 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 a - real c lower limit of integration c c b - real c upper limit of integration c c epsabs - real c absolute accuracy requested c epsrel - real c relative accuracy requested c if epsabs.le.0 c and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), c the routine will end with ier = 6. c c key - integer c key for choice of local integration rule c a gauss-kronrod pair is used with c 7 - 15 points if key.lt.2, c 10 - 21 points if key = 2, c 15 - 31 points if key = 3, c 20 - 41 points if key = 4, c 25 - 51 points if key = 5, c 30 - 61 points if key.gt.5. c c on return c result - real c approximation to the integral c c abserr - real c estimate of the modulus of the absolute error, c which should equal or exceed abs(i-result) c c neval - integer c number of integrand evaluations c c ier - integer c ier = 0 normal and reliable termination of the c routine. it is assumed that the requested c accuracy has been achieved. c ier.gt.0 abnormal termination of the routine c the estimates for result and error are c less reliable. it is assumed that the c requested accuracy has not been achieved. c error messages c ier = 1 maximum number of subdivisions allowed c has been achieved. one can allow more c subdivisions by increasing the value of c limit (and taking the according dimension c adjustments into account). however, if c this yield no improvement it is advised c to analyze the integrand in order to c determine the integration difficulaties. c if the position of a local difficulty can c be determined (i.e.singularity, c discontinuity within the interval) one c will probably gain from splitting up the c interval at this point and calling the c integrator on the subranges. if possible, c an appropriate special-purpose integrator c should be used which is designed for c handling the type of difficulty involved. c = 2 the occurrence of roundoff error is c detected, which prevents the requested c tolerance from being achieved. c = 3 extremely bad integrand behaviour occurs c at some points of the integration c interval. c = 6 the input is invalid, because c (epsabs.le.0 and c epsrel.lt.max(50*rel.mach.acc.,0.5d-28)) c or limit.lt.1 or lenw.lt.limit*4. c result, abserr, neval, last are set c to zero. c except when lenw is invalid, iwork(1), c work(limit*2+1) and work(limit*3+1) are c set to zero, work(1) is set to a and c work(limit+1) to b. c c dimensioning parameters c limit - integer c dimensioning parameter for iwork c limit determines the maximum number of subintervals c in the partition of the given integration interval c (a,b), limit.ge.1. c if limit.lt.1, the routine will end with ier = 6. c c lenw - integer c dimensioning parameter for work c lenw must be at least limit*4. c if lenw.lt.limit*4, the routine will end with c ier = 6. c c last - integer c on return, last equals the number of subintervals c produced in the subdivision process, which c determines the number of significant elements c actually in the work arrays. c c work arrays c iwork - integer c vector of dimension at least limit, the first k c elements of which contain pointers to the error c estimates over the subintervals, such that c work(limit*3+iwork(1)),... , work(limit*3+iwork(k)) c form a decreasing sequence with k = last if c last.le.(limit/2+2), and k = limit+1-last otherwise c c work - real c vector of dimension at least lenw c on return c work(1), ..., work(last) contain the left end c points of the subintervals in the partition of c (a,b), c work(limit+1), ..., work(limit+last) contain the c right end points, c work(limit*2+1), ..., work(limit*2+last) contain c the integral approximations over the subintervals, c work(limit*3+1), ..., work(limit*3+last) contain c the error estimates. c c***references (none) c***routines called qage,xerror c***end prologue qag c real a,abserr,b,epsabs,epsrel,f,result,work integer ier,iwork,key,lenw,limit,lvl,l1,l2,l3,neval c dimension iwork(limit),work(lenw) c external f c c check validity of lenw. c c***first executable statement qag ier = 6 neval = 0 last = 0 result = 0.0e+00 abserr = 0.0e+00 if(limit.lt.1.or.lenw.lt.limit*4) go to 10 c c prepare call for qage. c l1 = limit+1 l2 = limit+l1 l3 = limit+l2 c call qage(f,a,b,epsabs,epsrel,key,limit,result,abserr,neval, * ier,work(1),work(l1),work(l2),work(l3),iwork,last) c c call error handler if necessary. c lvl = 0 10 if(ier.eq.6) lvl = 1 if(ier.ne.0) call xerror(26habnormal return from qag , * 26,ier,lvl) return end subroutine qage(f,a,b,epsabs,epsrel,key,limit,result,abserr, * neval,ier,alist,blist,rlist,elist,iord,last) c***begin prologue qage c***date written 800101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a1a1 c***keywords automatic integrator, general-purpose, c integrand examinator, globally adaptive, c gauss-kronrod c***author piessens,robert,appl. math. & progr. div. - k.u.leuven c de doncker,elise,appl. math. & progr. div. - k.u.leuven c***purpose the routine calculates an approximation result to a given c definite integral i = integral of f over (a,b), c hopefully satisfying following claim for accuracy c abs(i-reslt).le.max(epsabs,epsrel*abs(i)). c***description c c computation of a definite integral 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 a - real c lower limit of integration c c b - real c upper limit of integration c c epsabs - real c absolute accuracy requested c epsrel - real c relative accuracy requested c if epsabs.le.0 c and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), c the routine will end with ier = 6. c c key - integer c key for choice of local integration rule c a gauss-kronrod pair is used with c 7 - 15 points if key.lt.2, c 10 - 21 points if key = 2, c 15 - 31 points if key = 3, c 20 - 41 points if key = 4, c 25 - 51 points if key = 5, c 30 - 61 points if key.gt.5. c c limit - integer c gives an upperbound on the number of subintervals c in the partition of (a,b), limit.ge.1. c c on return c result - real c approximation to the integral c c abserr - real c estimate of the modulus of the absolute error, c which should equal or exceed abs(i-result) c c neval - integer c number of integrand evaluations c c ier - integer c ier = 0 normal and reliable termination of the c routine. it is assumed that the requested c accuracy has been achieved. c ier.gt.0 abnormal termination of the routine c the estimates for result and error are c less reliable. it is assumed that the c requested accuracy has not been achieved. c error messages c ier = 1 maximum number of subdivisions allowed c has been achieved. one can allow more c subdivisions by increasing the value c of limit. c however, if this yields no improvement it c is rather advised to analyze the integrand c in order to determine the integration c difficulties. if the position of a local c difficulty can be determined(e.g. c singularity, discontinuity within the c interval) one will probably gain from c splitting up the interval at this point c and calling the integrator on the c subranges. if possible, an appropriate c special-purpose integrator should be used c which is designed for handling the type of c difficulty involved. c = 2 the occurrence of roundoff error is c detected, which prevents the requested c tolerance from being achieved. c = 3 extremely bad integrand behaviour occurs c at some points of the integration c interval. c = 6 the input is invalid, because c (epsabs.le.0 and c epsrel.lt.max(50*rel.mach.acc.,0.5d-28), c result, abserr, neval, last, rlist(1) , c elist(1) and iord(1) are set to zero. c alist(1) and blist(1) are set to a and b c respectively. c c alist - real c vector of dimension at least limit, the first c last elements of which are the left c end points of the subintervals in the partition c of the given integration range (a,b) c c blist - real c vector of dimension at least limit, the first c last elements of which are the right c end points of the subintervals in the partition c of the given integration range (a,b) c c rlist - real c vector of dimension at least limit, the first c last elements of which are the c integral approximations on the subintervals c c elist - real c vector of dimension at least limit, the first c last elements of which are the moduli of the c absolute error estimates on the subintervals c c iord - integer c vector of dimension at least limit, the first k c elements of which are pointers to the c error estimates over the subintervals, c such that elist(iord(1)), ..., c elist(iord(k)) form a decreasing sequence, c with k = last if last.le.(limit/2+2), and c k = limit+1-last otherwise c c last - integer c number of subintervals actually produced in the c subdivision process c c***references (none) c***routines called qk15,qk21,qk31,qk41,qk51,qk61,qpsrt,r1mach c***end prologue qage c real a,abserr,alist,area,area1,area12,area2,a1,a2,b,blist, * b1,b2,defabs,defab1,defab2,r1mach,elist,epmach, * epsabs,epsrel,errbnd,errmax,error1,error2,erro12,errsum,f, * resabs,result,rlist,uflow integer ier,iord,iroff1,iroff2,k,key,keyf,last, * limit,maxerr,neval,nrmax c dimension alist(limit),blist(limit),elist(limit),iord(limit), * rlist(limit) c external f c c list of major variables c ----------------------- c c alist - list of left end points of all subintervals c considered up to now c blist - list of right end points of all subintervals c considered up to now c rlist(i) - approximation to the integral over c (alist(i),blist(i)) c elist(i) - error estimate applying to rlist(i) c maxerr - pointer to the interval with largest c error estimate c errmax - elist(maxerr) c area - sum of the integrals over the subintervals c errsum - sum of the errors over the subintervals c errbnd - requested accuracy max(epsabs,epsrel* c abs(result)) c *****1 - variable for the left subinterval c *****2 - variable for the right subinterval c last - index for subdivision c 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 qage epmach = r1mach(4) uflow = r1mach(1) c c test on validity of parameters c ------------------------------ c ier = 0 neval = 0 last = 0 result = 0.0e+00 abserr = 0.0e+00 alist(1) = a blist(1) = b rlist(1) = 0.0e+00 elist(1) = 0.0e+00 iord(1) = 0 if(epsabs.le.0.0e+00.and. * epsrel.lt.amax1(0.5e+02*epmach,0.5e-14)) ier = 6 if(ier.eq.6) go to 999 c c first approximation to the integral c ----------------------------------- c keyf = key if(key.le.0) keyf = 1 if(key.ge.7) keyf = 6 neval = 0 if(keyf.eq.1) call qk15(f,a,b,result,abserr,defabs,resabs) if(keyf.eq.2) call qk21(f,a,b,result,abserr,defabs,resabs) if(keyf.eq.3) call qk31(f,a,b,result,abserr,defabs,resabs) if(keyf.eq.4) call qk41(f,a,b,result,abserr,defabs,resabs) if(keyf.eq.5) call qk51(f,a,b,result,abserr,defabs,resabs) if(keyf.eq.6) call qk61(f,a,b,result,abserr,defabs,resabs) last = 1 rlist(1) = result elist(1) = abserr iord(1) = 1 c c test on accuracy. c errbnd = amax1(epsabs,epsrel*abs(result)) if(abserr.le.0.5e+02*epmach*defabs.and.abserr.gt. * errbnd) ier = 2 if(limit.eq.1) ier = 1 if(ier.ne.0.or.(abserr.le.errbnd.and.abserr.ne.resabs) * .or.abserr.eq.0.0e+00) go to 60 c c initialization c -------------- c c errmax = abserr maxerr = 1 area = result errsum = abserr nrmax = 1 iroff1 = 0 iroff2 = 0 c c main do-loop c ------------ c do 30 last = 2,limit c c bisect the subinterval with the largest error estimate. c a1 = alist(maxerr) b1 = 0.5e+00*(alist(maxerr)+blist(maxerr)) a2 = b1 b2 = blist(maxerr) if(keyf.eq.1) call qk15(f,a1,b1,area1,error1,resabs,defab1) if(keyf.eq.2) call qk21(f,a1,b1,area1,error1,resabs,defab1) if(keyf.eq.3) call qk31(f,a1,b1,area1,error1,resabs,defab1) if(keyf.eq.4) call qk41(f,a1,b1,area1,error1,resabs,defab1) if(keyf.eq.5) call qk51(f,a1,b1,area1,error1,resabs,defab1) if(keyf.eq.6) call qk61(f,a1,b1,area1,error1,resabs,defab1) if(keyf.eq.1) call qk15(f,a2,b2,area2,error2,resabs,defab2) if(keyf.eq.2) call qk21(f,a2,b2,area2,error2,resabs,defab2) if(keyf.eq.3) call qk31(f,a2,b2,area2,error2,resabs,defab2) if(keyf.eq.4) call qk41(f,a2,b2,area2,error2,resabs,defab2) if(keyf.eq.5) call qk51(f,a2,b2,area2,error2,resabs,defab2) if(keyf.eq.6) call qk61(f,a2,b2,area2,error2,resabs,defab2) c c improve previous approximations to integral c and error and test for accuracy. c neval = neval+1 area12 = area1+area2 erro12 = error1+error2 errsum = errsum+erro12-errmax area = area+area12-rlist(maxerr) if(defab1.eq.error1.or.defab2.eq.error2) go to 5 if(abs(rlist(maxerr)-area12).le.0.1e-04*abs(area12) * .and.erro12.ge.0.99e+00*errmax) iroff1 = iroff1+1 if(last.gt.10.and.erro12.gt.errmax) iroff2 = iroff2+1 5 rlist(maxerr) = area1 rlist(last) = area2 errbnd = amax1(epsabs,epsrel*abs(area)) if(errsum.le.errbnd) go to 8 c c test for roundoff error and eventually c set error flag. c if(iroff1.ge.6.or.iroff2.ge.20) ier = 2 c c set error flag in the case that the number of c subintervals equals limit. c if(last.eq.limit) ier = 1 c c set error flag in the case of bad integrand behaviour c at a point of the integration range. c if(amax1(abs(a1),abs(b2)).le.(0.1e+01+0.1e+03* * epmach)*(abs(a2)+0.1e+04*uflow)) ier = 3 c c append the newly-created intervals to the list. c 8 if(error2.gt.error1) go to 10 alist(last) = a2 blist(maxerr) = b1 blist(last) = b2 elist(maxerr) = error1 elist(last) = error2 go to 20 10 alist(maxerr) = a2 alist(last) = a1 blist(last) = b1 rlist(maxerr) = area2 rlist(last) = area1 elist(maxerr) = error2 elist(last) = error1 c c call subroutine qpsrt to maintain the descending ordering c in the list of error estimates and select the c subinterval with the largest error estimate (to be c bisected next). c 20 call qpsrt(limit,last,maxerr,errmax,elist,iord,nrmax) c ***jump out of do-loop if(ier.ne.0.or.errsum.le.errbnd) go to 40 30 continue c c compute final result. c --------------------- c 40 result = 0.0e+00 do 50 k=1,last result = result+rlist(k) 50 continue abserr = errsum 60 if(keyf.ne.1) neval = (10*keyf+1)*(2*neval+1) if(keyf.eq.1) neval = 30*neval+15 999 return end subroutine qk15(f,a,b,result,abserr,resabs,resasc) c***begin prologue qk15 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 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 calling program. 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 the7-point gauss rule(resg). c c abserr - real c estimate of the modulus of the absolute error, c which should not exceed abs(i-result) c c resabs - real c approximation to the integral j c c resasc - real c approximation to the integral of abs(f-i/(b-a)) c over (a,b) c c***references (none) c***routines called r1mach c***end prologue qk15 c real a,absc,abserr,b,centr,dhlgth,epmach,f,fc,fsum,fval1,fval2, * fv1,fv2,hlgth,resabs,resasc,resg,resk,reskh,result,r1mach,uflow, * wg,wgk,xgk integer j,jtw,jtwm1 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 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.0e+00 / 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/ 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 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 qk15 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 c the integral, and estimate the absolute error. c fc = f(centr) resg = fc*wg(4) resk = fc*wgk(8) resabs = abs(resk) do 10 j=1,3 jtw = j*2 absc = hlgth*xgk(jtw) fval1 = f(centr-absc) fval2 = f(centr+absc) 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) fval1 = f(centr-absc) fval2 = f(centr+absc) 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 subroutine qk21(f,a,b,result,abserr,resabs,resasc) c***begin prologue qk21 c***date written 800101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a1a2 c***keywords 21-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 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 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 21-point c kronrod rule (resk) obtained by optimal addition c of abscissae to the 10-point gauss rule (resg). c c abserr - real c estimate of the modulus of the absolute error, c which should not exceed abs(i-result) c c resabs - real c approximation to the integral j c c resasc - real c approximation to the integral of abs(f-i/(b-a)) c over (a,b) c c***references (none) c***routines called r1mach c***end prologue qk21 c real a,absc,abserr,b,centr,dhlgth,epmach,f,fc,fsum,fval1,fval2, * fv1,fv2,hlgth,resabs,resg,resk,reskh,result,r1mach,uflow,wg,wgk, * xgk integer j,jtw,jtwm1 external f c dimension fv1(10),fv2(10),wg(5),wgk(11),xgk(11) 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 21-point kronrod rule c xgk(2), xgk(4), ... abscissae of the 10-point c gauss rule c xgk(1), xgk(3), ... abscissae which are optimally c added to the 10-point gauss rule c c wgk - weights of the 21-point kronrod rule c c wg - weights of the 10-point gauss rule c data xgk(1),xgk(2),xgk(3),xgk(4),xgk(5),xgk(6),xgk(7), * xgk(8),xgk(9),xgk(10),xgk(11)/ * 0.9956571630258081e+00, 0.9739065285171717e+00, * 0.9301574913557082e+00, 0.8650633666889845e+00, * 0.7808177265864169e+00, 0.6794095682990244e+00, * 0.5627571346686047e+00, 0.4333953941292472e+00, * 0.2943928627014602e+00, 0.1488743389816312e+00, * 0.0000000000000000e+00/ c data wgk(1),wgk(2),wgk(3),wgk(4),wgk(5),wgk(6),wgk(7), * wgk(8),wgk(9),wgk(10),wgk(11)/ * 0.1169463886737187e-01, 0.3255816230796473e-01, * 0.5475589657435200e-01, 0.7503967481091995e-01, * 0.9312545458369761e-01, 0.1093871588022976e+00, * 0.1234919762620659e+00, 0.1347092173114733e+00, * 0.1427759385770601e+00, 0.1477391049013385e+00, * 0.1494455540029169e+00/ c data wg(1),wg(2),wg(3),wg(4),wg(5)/ * 0.6667134430868814e-01, 0.1494513491505806e+00, * 0.2190863625159820e+00, 0.2692667193099964e+00, * 0.2955242247147529e+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 10-point gauss formula c resk - result of the 21-point kronrod formula c reskh - approximation to the mean value of f over (a,b), c i.e. to i/(b-a) c 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 qk21 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 21-point kronrod approximation to c the integral, and estimate the absolute error. c resg = 0.0e+00 fc = f(centr) resk = wgk(11)*fc resabs = abs(resk) do 10 j=1,5 jtw = 2*j absc = hlgth*xgk(jtw) fval1 = f(centr-absc) fval2 = f(centr+absc) 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,5 jtwm1 = 2*j-1 absc = hlgth*xgk(jtwm1) fval1 = f(centr-absc) fval2 = f(centr+absc) 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(11)*abs(fc-reskh) do 20 j=1,10 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 subroutine qk31(f,a,b,result,abserr,resabs,resasc) c***begin prologue qk31 c***date written 800101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a1a2 c***keywords 31-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 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 calling program. 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 31-point c gauss-kronrod rule (resk), obtained by optimal c addition of abscissae to the 15-point gauss c rule (resg). c c abserr - real c estimate of the modulus of the modulus, c which should not exceed abs(i-result) c c resabs - real c approximation to the integral j c c resasc - real c approximation to the integral of abs(f-i/(b-a)) c over (a,b) c c***references (none) c***routines called r1mach c***end prologue qk31 real a,absc,abserr,b,centr,dhlgth,epmach,f,fc,fsum,fval1,fval2, * fv1,fv2,hlgth,resabs,resasc,resg,resk,reskh,result,r1mach,uflow, * wg,wgk,xgk integer j,jtw,jtwm1 external f c dimension fv1(15),fv2(15),xgk(16),wgk(16),wg(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 31-point kronrod rule c xgk(2), xgk(4), ... abscissae of the 15-point c gauss rule c xgk(1), xgk(3), ... abscissae which are optimally c added to the 15-point gauss rule c c wgk - weights of the 31-point kronrod rule c c wg - weights of the 15-point gauss rule c data xgk(1),xgk(2),xgk(3),xgk(4),xgk(5),xgk(6),xgk(7),xgk(8), * xgk(9),xgk(10),xgk(11),xgk(12),xgk(13),xgk(14),xgk(15), * xgk(16)/ * 0.9980022986933971e+00, 0.9879925180204854e+00, * 0.9677390756791391e+00, 0.9372733924007059e+00, * 0.8972645323440819e+00, 0.8482065834104272e+00, * 0.7904185014424659e+00, 0.7244177313601700e+00, * 0.6509967412974170e+00, 0.5709721726085388e+00, * 0.4850818636402397e+00, 0.3941513470775634e+00, * 0.2991800071531688e+00, 0.2011940939974345e+00, * 0.1011420669187175e+00, 0.0e+00 / data wgk(1),wgk(2),wgk(3),wgk(4),wgk(5),wgk(6),wgk(7),wgk(8), * wgk(9),wgk(10),wgk(11),wgk(12),wgk(13),wgk(14),wgk(15), * wgk(16)/ * 0.5377479872923349e-02, 0.1500794732931612e-01, * 0.2546084732671532e-01, 0.3534636079137585e-01, * 0.4458975132476488e-01, 0.5348152469092809e-01, * 0.6200956780067064e-01, 0.6985412131872826e-01, * 0.7684968075772038e-01, 0.8308050282313302e-01, * 0.8856444305621177e-01, 0.9312659817082532e-01, * 0.9664272698362368e-01, 0.9917359872179196e-01, * 0.1007698455238756e+00, 0.1013300070147915e+00/ data wg(1),wg(2),wg(3),wg(4),wg(5),wg(6),wg(7),wg(8)/ * 0.3075324199611727e-01, 0.7036604748810812e-01, * 0.1071592204671719e+00, 0.1395706779261543e+00, * 0.1662692058169939e+00, 0.1861610000155622e+00, * 0.1984314853271116e+00, 0.2025782419255613e+00/ c c c list of major variables 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 15-point gauss formula c resk - result of the 31-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 epmach is the largest relative spacing. c uflow is the smallest positive magnitude. c c***first executable statement qk31 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 31-point kronrod approximation to c the integral, and estimate the absolute error. c fc = f(centr) resg = wg(8)*fc resk = wgk(16)*fc resabs = abs(resk) do 10 j=1,7 jtw = j*2 absc = hlgth*xgk(jtw) fval1 = f(centr-absc) fval2 = f(centr+absc) 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,8 jtwm1 = j*2-1 absc = hlgth*xgk(jtwm1) fval1 = f(centr-absc) fval2 = f(centr+absc) 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(16)*abs(fc-reskh) do 20 j=1,15 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 subroutine qk41(f,a,b,result,abserr,resabs,resasc) c***begin prologue qk41 c***date written 800101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a1a2 c***keywords 41-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 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 calling program. 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 41-point c gauss-kronrod rule (resk) obtained by optimal c addition of abscissae to the 20-point gauss c rule (resg). c c abserr - real c estimate of the modulus of the absolute error, c which should not exceed abs(i-result) c c resabs - real c approximation to the integral j c c resasc - real c approximation to the integal of abs(f-i/(b-a)) c over (a,b) c c***references (none) c***routines called r1mach c***end prologue qk41 c real a,absc,abserr,b,centr,dhlgth,epmach,f,fc,fsum,fval1,fval2, * fv1,fv2,hlgth,resabs, * resasc,resg,resk,reskh,result,r1mach,uflow, * wg,wgk,xgk integer j,jtw,jtwm1 external f c dimension fv1(20),fv2(20),xgk(21),wgk(21),wg(10) 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 41-point gauss-kronrod rule c xgk(2), xgk(4), ... abscissae of the 20-point c gauss rule c xgk(1), xgk(3), ... abscissae which are optimally c added to the 20-point gauss rule c c wgk - weights of the 41-point gauss-kronrod rule c c wg - weights of the 20-point gauss rule c data xgk(1),xgk(2),xgk(3),xgk(4),xgk(5),xgk(6),xgk(7),xgk(8), * xgk(9),xgk(10),xgk(11),xgk(12),xgk(13),xgk(14),xgk(15), * xgk(16),xgk(17),xgk(18),xgk(19),xgk(20),xgk(21)/ * 0.9988590315882777e+00, 0.9931285991850949e+00, * 0.9815078774502503e+00, 0.9639719272779138e+00, * 0.9408226338317548e+00, 0.9122344282513259e+00, * 0.8782768112522820e+00, 0.8391169718222188e+00, * 0.7950414288375512e+00, 0.7463319064601508e+00, * 0.6932376563347514e+00, 0.6360536807265150e+00, * 0.5751404468197103e+00, 0.5108670019508271e+00, * 0.4435931752387251e+00, 0.3737060887154196e+00, * 0.3016278681149130e+00, 0.2277858511416451e+00, * 0.1526054652409227e+00, 0.7652652113349733e-01, * 0.0e+00 / data wgk(1),wgk(2),wgk(3),wgk(4),wgk(5),wgk(6),wgk(7),wgk(8), * wgk(9),wgk(10),wgk(11),wgk(12),wgk(13),wgk(14),wgk(15),wgk(16), * wgk(17),wgk(18),wgk(19),wgk(20),wgk(21)/ * 0.3073583718520532e-02, 0.8600269855642942e-02, * 0.1462616925697125e-01, 0.2038837346126652e-01, * 0.2588213360495116e-01, 0.3128730677703280e-01, * 0.3660016975820080e-01, 0.4166887332797369e-01, * 0.4643482186749767e-01, 0.5094457392372869e-01, * 0.5519510534828599e-01, 0.5911140088063957e-01, * 0.6265323755478117e-01, 0.6583459713361842e-01, * 0.6864867292852162e-01, 0.7105442355344407e-01, * 0.7303069033278667e-01, 0.7458287540049919e-01, * 0.7570449768455667e-01, 0.7637786767208074e-01, * 0.7660071191799966e-01/ data wg(1),wg(2),wg(3),wg(4),wg(5),wg(6),wg(7),wg(8),wg(9),wg(10)/ * 0.1761400713915212e-01, 0.4060142980038694e-01, * 0.6267204833410906e-01, 0.8327674157670475e-01, * 0.1019301198172404e+00, 0.1181945319615184e+00, * 0.1316886384491766e+00, 0.1420961093183821e+00, * 0.1491729864726037e+00, 0.1527533871307259e+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 20-point gauss formula c resk - result of the 41-point kronrod formula c reskh - approximation to mean value of f over (a,b), i.e. c 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 qk41 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 41-point gauss-kronrod approximation to c the integral, and estimate the absolute error. c resg = 0.0e+00 fc = f(centr) resk = wgk(21)*fc resabs = abs(resk) do 10 j=1,10 jtw = j*2 absc = hlgth*xgk(jtw) fval1 = f(centr-absc) fval2 = f(centr+absc) 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,10 jtwm1 = j*2-1 absc = hlgth*xgk(jtwm1) fval1 = f(centr-absc) fval2 = f(centr+absc) 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(21)*abs(fc-reskh) do 20 j=1,20 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.e+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 subroutine qk51(f,a,b,result,abserr,resabs,resasc) c***begin prologue qk51 c***date written 800101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a1a2 c***keywords 51-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 real version c c parameters c on entry c f - real c function subroutine 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 - 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 51-point c kronrod rule (resk) obtained by optimal addition c of abscissae to the 25-point gauss rule (resg). c c abserr - real c estimate of the modulus of the absolute error, c which should not exceed abs(i-result) c c resabs - real c approximation to the integral j c c resasc - real c approximation to the integral of abs(f-i/(b-a)) c over (a,b) c c***references (none) c***routines called r1mach c***end prologue qk51 c real a,absc,abserr,b,centr,dhlgth,epmach,f,fc,fsum,fval1,fval2, * fv1,fv2,hlgth,resabs,resasc,resg,resk,reskh,result,r1mach,uflow, * wg,wgk,xgk integer j,jtw,jtwm1 external f c dimension fv1(25),fv2(25),xgk(26),wgk(26),wg(13) 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 51-point kronrod rule c xgk(2), xgk(4), ... abscissae of the 25-point c gauss rule c xgk(1), xgk(3), ... abscissae which are optimally c added to the 25-point gauss rule c c wgk - weights of the 51-point kronrod rule c c wg - weights of the 25-point gauss rule c data xgk(1),xgk(2),xgk(3),xgk(4),xgk(5),xgk(6),xgk(7),xgk(8), * xgk(9),xgk(10),xgk(11),xgk(12),xgk(13),xgk(14)/ * 0.9992621049926098e+00, 0.9955569697904981e+00, * 0.9880357945340772e+00, 0.9766639214595175e+00, * 0.9616149864258425e+00, 0.9429745712289743e+00, * 0.9207471152817016e+00, 0.8949919978782754e+00, * 0.8658470652932756e+00, 0.8334426287608340e+00, * 0.7978737979985001e+00, 0.7592592630373576e+00, * 0.7177664068130844e+00, 0.6735663684734684e+00/ data xgk(15),xgk(16),xgk(17),xgk(18),xgk(19),xgk(20),xgk(21), * xgk(22),xgk(23),xgk(24),xgk(25),xgk(26)/ * 0.6268100990103174e+00, 0.5776629302412230e+00, * 0.5263252843347192e+00, 0.4730027314457150e+00, * 0.4178853821930377e+00, 0.3611723058093878e+00, * 0.3030895389311078e+00, 0.2438668837209884e+00, * 0.1837189394210489e+00, 0.1228646926107104e+00, * 0.6154448300568508e-01, 0.0e+00 / data wgk(1),wgk(2),wgk(3),wgk(4),wgk(5),wgk(6),wgk(7),wgk(8), * wgk(9),wgk(10),wgk(11),wgk(12),wgk(13),wgk(14)/ * 0.1987383892330316e-02, 0.5561932135356714e-02, * 0.9473973386174152e-02, 0.1323622919557167e-01, * 0.1684781770912830e-01, 0.2043537114588284e-01, * 0.2400994560695322e-01, 0.2747531758785174e-01, * 0.3079230016738749e-01, 0.3400213027432934e-01, * 0.3711627148341554e-01, 0.4008382550403238e-01, * 0.4287284502017005e-01, 0.4550291304992179e-01/ data wgk(15),wgk(16),wgk(17),wgk(18),wgk(19),wgk(20),wgk(21) * ,wgk(22),wgk(23),wgk(24),wgk(25),wgk(26)/ * 0.4798253713883671e-01, 0.5027767908071567e-01, * 0.5236288580640748e-01, 0.5425112988854549e-01, * 0.5595081122041232e-01, 0.5743711636156783e-01, * 0.5868968002239421e-01, 0.5972034032417406e-01, * 0.6053945537604586e-01, 0.6112850971705305e-01, * 0.6147118987142532e-01, 0.6158081806783294e-01/ data wg(1),wg(2),wg(3),wg(4),wg(5),wg(6),wg(7),wg(8),wg(9), * wg(10),wg(11),wg(12),wg(13)/ * 0.1139379850102629e-01, 0.2635498661503214e-01, * 0.4093915670130631e-01, 0.5490469597583519e-01, * 0.6803833381235692e-01, 0.8014070033500102e-01, * 0.9102826198296365e-01, 0.1005359490670506e+00, * 0.1085196244742637e+00, 0.1148582591457116e+00, * 0.1194557635357848e+00, 0.1222424429903100e+00, * 0.1231760537267155e+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 25-point gauss formula c resk - result of the 51-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 qk51 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 51-point kronrod approximation to c the integral, and estimate the absolute error. c fc = f(centr) resg = wg(13)*fc resk = wgk(26)*fc resabs = abs(resk) do 10 j=1,12 jtw = j*2 absc = hlgth*xgk(jtw) fval1 = f(centr-absc) fval2 = f(centr+absc) 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,13 jtwm1 = j*2-1 absc = hlgth*xgk(jtwm1) fval1 = f(centr-absc) fval2 = f(centr+absc) 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(26)*abs(fc-reskh) do 20 j=1,25 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 subroutine qk61(f,a,b,result,abserr,resabs,resasc) c***begin prologue qk61 c***date written 800101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a1a2 c***keywords 61-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 dabs(f) over (a,b) c***description c c integration rule c standard fortran subroutine c real version c 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 calling program. 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 61-point c kronrod rule (resk) obtained by optimal addition of c abscissae to the 30-point gauss rule (resg). c c abserr - real c estimate of the modulus of the absolute error, c which should equal or exceed dabs(i-result) c c resabs - real c approximation to the integral j c c resasc - real c approximation to the integral of dabs(f-i/(b-a)) c c c***references (none) c***routines called r1mach c***end prologue qk61 c real a,absc,abserr,b,centr,dhlgth,epmach,f,fc,fsum,fval1,fval2, * fv1,fv2,hlgth,resabs,resasc,resg,resk,reskh,result,r1mach,uflow, * wg,wgk,xgk integer j,jtw,jtwm1 external f c dimension fv1(30),fv2(30),xgk(31),wgk(31),wg(15) c c the abscissae and weights are given for the c interval (-1,1). because of symmetry only the positive c abscissae and their corresponding weights are given. c c xgk - abscissae of the 61-point kronrod rule c xgk(2), xgk(4) ... abscissae of the 30-point c gauss rule c xgk(1), xgk(3) ... optimally added abscissae c to the 30-point gauss rule c c wgk - weights of the 61-point kronrod rule c c wg - weigths of the 30-point gauss rule c data xgk(1),xgk(2),xgk(3),xgk(4),xgk(5),xgk(6),xgk(7),xgk(8), * xgk(9),xgk(10)/ * 0.9994844100504906e+00, 0.9968934840746495e+00, * 0.9916309968704046e+00, 0.9836681232797472e+00, * 0.9731163225011263e+00, 0.9600218649683075e+00, * 0.9443744447485600e+00, 0.9262000474292743e+00, * 0.9055733076999078e+00, 0.8825605357920527e+00/ data xgk(11),xgk(12),xgk(13),xgk(14),xgk(15),xgk(16), * xgk(17),xgk(18),xgk(19),xgk(20)/ * 0.8572052335460611e+00, 0.8295657623827684e+00, * 0.7997278358218391e+00, 0.7677774321048262e+00, * 0.7337900624532268e+00, 0.6978504947933158e+00, * 0.6600610641266270e+00, 0.6205261829892429e+00, * 0.5793452358263617e+00, 0.5366241481420199e+00/ data xgk(21),xgk(22),xgk(23),xgk(24), * xgk(25),xgk(26),xgk(27),xgk(28),xgk(29),xgk(30),xgk(31)/ * 0.4924804678617786e+00, 0.4470337695380892e+00, * 0.4004012548303944e+00, 0.3527047255308781e+00, * 0.3040732022736251e+00, 0.2546369261678898e+00, * 0.2045251166823099e+00, 0.1538699136085835e+00, * 0.1028069379667370e+00, 0.5147184255531770e-01, * 0.0e+00 / data wgk(1),wgk(2),wgk(3),wgk(4),wgk(5),wgk(6),wgk(7),wgk(8), * wgk(9),wgk(10)/ * 0.1389013698677008e-02, 0.3890461127099884e-02, * 0.6630703915931292e-02, 0.9273279659517763e-02, * 0.1182301525349634e-01, 0.1436972950704580e-01, * 0.1692088918905327e-01, 0.1941414119394238e-01, * 0.2182803582160919e-01, 0.2419116207808060e-01/ data wgk(11),wgk(12),wgk(13),wgk(14),wgk(15),wgk(16), * wgk(17),wgk(18),wgk(19),wgk(20)/ * 0.2650995488233310e-01, 0.2875404876504129e-01, * 0.3090725756238776e-01, 0.3298144705748373e-01, * 0.3497933802806002e-01, 0.3688236465182123e-01, * 0.3867894562472759e-01, 0.4037453895153596e-01, * 0.4196981021516425e-01, 0.4345253970135607e-01/ data wgk(21),wgk(22),wgk(23),wgk(24), * wgk(25),wgk(26),wgk(27),wgk(28),wgk(29),wgk(30),wgk(31)/ * 0.4481480013316266e-01, 0.4605923827100699e-01, * 0.4718554656929915e-01, 0.4818586175708713e-01, * 0.4905543455502978e-01, 0.4979568342707421e-01, * 0.5040592140278235e-01, 0.5088179589874961e-01, * 0.5122154784925877e-01, 0.5142612853745903e-01, * 0.5149472942945157e-01/ data wg(1),wg(2),wg(3),wg(4),wg(5),wg(6),wg(7),wg(8)/ * 0.7968192496166606e-02, 0.1846646831109096e-01, * 0.2878470788332337e-01, 0.3879919256962705e-01, * 0.4840267283059405e-01, 0.5749315621761907e-01, * 0.6597422988218050e-01, 0.7375597473770521e-01/ data wg(9),wg(10),wg(11),wg(12),wg(13),wg(14),wg(15)/ * 0.8075589522942022e-01, 0.8689978720108298e-01, * 0.9212252223778613e-01, 0.9636873717464426e-01, * 0.9959342058679527e-01, 0.1017623897484055e+00, * 0.1028526528935588e+00/ 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 30-point gauss rule c resk - result of the 61-point kronrod rule c reskh - approximation to the mean value of f c over (a,b), 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 qk61 epmach = r1mach(4) uflow = r1mach(1) c centr = 0.5e+00*(b+a) hlgth = 0.5e+00*(b-a) dhlgth = abs(hlgth) c c compute the 61-point kronrod approximation to the c integral, and estimate the absolute error. c resg = 0.0e+00 fc = f(centr) resk = wgk(31)*fc resabs = abs(resk) do 10 j=1,15 jtw = j*2 absc = hlgth*xgk(jtw) fval1 = f(centr-absc) fval2 = f(centr+absc) 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,15 jtwm1 = j*2-1 absc = hlgth*xgk(jtwm1) fval1 = f(centr-absc) fval2 = f(centr+absc) 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(31)*abs(fc-reskh) do 20 j=1,30 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 subroutine qpsrt(limit,last,maxerr,ermax,elist,iord,nrmax) c***begin prologue qpsrt c***refer to qage,qagie,qagpe,qagse,qawce,qawse,qawoe c***routines called (none) c***keywords sequential sorting c***description c c 1. qpsrt c ordering routine c standard fortran subroutine c real version c c 2. purpose c this routine maintains the descending ordering c in the list of the local error estimates resulting from c the interval subdivision process. at each call two error c estimates are inserted using the sequential search c method, top-down for the largest error estimate c and bottom-up for the smallest error estimate. c c 3. calling sequence c call qpsrt(limit,last,maxerr,ermax,elist,iord,nrmax) c c parameters (meaning at output) c limit - integer c maximum number of error estimates the list c can contain c c last - integer c number of error estimates currently c in the list c c maxerr - integer c maxerr points to the nrmax-th largest error c estimate currently in the list c c ermax - real c nrmax-th largest error estimate c ermax = elist(maxerr) c c elist - real c vector of dimension last containing c the error estimates c c iord - integer c vector of dimension last, the first k c elements of which contain pointers c to the error estimates, such that c elist(iord(1)),... , elist(iord(k)) c form a decreasing sequence, with c k = last if last.le.(limit/2+2), and c k = limit+1-last otherwise c c nrmax - integer c maxerr = iord(nrmax) c c 4. no subroutines or functions needed c***end prologue qpsrt c real elist,ermax,errmax,errmin integer i,ibeg,ido,iord,isucc,j,jbnd,jupbn,k,last,limit,maxerr, * nrmax dimension elist(last),iord(last) c c check whether the list contains more than c two error estimates. c c***first executable statement qpsrt if(last.gt.2) go to 10 iord(1) = 1 iord(2) = 2 go to 90 c c this part of the routine is only executed c if, due to a difficult integrand, subdivision c increased the error estimate. in the normal case c the insert procedure should start after the c nrmax-th largest error estimate. c 10 errmax = elist(maxerr) if(nrmax.eq.1) go to 30 ido = nrmax-1 do 20 i = 1,ido isucc = iord(nrmax-1) c ***jump out of do-loop if(errmax.le.elist(isucc)) go to 30 iord(nrmax) = isucc nrmax = nrmax-1 20 continue c c compute the number of elements in the list to c be maintained in descending order. this number c depends on the number of subdivisions still c allowed. c 30 jupbn = last if(last.gt.(limit/2+2)) jupbn = limit+3-last errmin = elist(last) c c insert errmax by traversing the list top-down, c starting comparison from the element elist(iord(nrmax+1)). c jbnd = jupbn-1 ibeg = nrmax+1 if(ibeg.gt.jbnd) go to 50 do 40 i=ibeg,jbnd isucc = iord(i) c ***jump out of do-loop if(errmax.ge.elist(isucc)) go to 60 iord(i-1) = isucc 40 continue 50 iord(jbnd) = maxerr iord(jupbn) = last go to 90 c c insert errmin by traversing the list bottom-up. c 60 iord(i-1) = maxerr k = jbnd do 70 j=i,jbnd isucc = iord(k) c ***jump out of do-loop if(errmin.lt.elist(isucc)) go to 80 iord(k+1) = isucc k = k-1 70 continue iord(i) = last go to 90 80 iord(k+1) = last c c set maxerr and ermax. c 90 maxerr = iord(nrmax) ermax = elist(maxerr) return end REAL FUNCTION R1MACH(I) INTEGER I C C SINGLE-PRECISION MACHINE CONSTANTS C R1MACH(1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE. C R1MACH(2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE. C R1MACH(3) = B**(-T), THE SMALLEST RELATIVE SPACING. C R1MACH(4) = B**(1-T), THE LARGEST RELATIVE SPACING. C R1MACH(5) = LOG10(B) C INTEGER SMALL(2) INTEGER LARGE(2) INTEGER RIGHT(2) INTEGER DIVER(2) INTEGER LOG10(2) C needs to be (2) for AUTODOUBLE, HARRIS SLASH 6, ... INTEGER SC SAVE SMALL, LARGE, RIGHT, DIVER, LOG10, SC REAL RMACH(5) EQUIVALENCE (RMACH(1),SMALL(1)) EQUIVALENCE (RMACH(2),LARGE(1)) EQUIVALENCE (RMACH(3),RIGHT(1)) EQUIVALENCE (RMACH(4),DIVER(1)) EQUIVALENCE (RMACH(5),LOG10(1)) INTEGER J, K, L, T3E(3) DATA T3E(1) / 9777664 / DATA T3E(2) / 5323660 / DATA T3E(3) / 46980 / C THIS VERSION ADAPTS AUTOMATICALLY TO MOST CURRENT MACHINES, C INCLUDING AUTO-DOUBLE COMPILERS. C TO COMPILE ON OLDER MACHINES, ADD A C IN COLUMN 1 C ON THE NEXT LINE DATA SC/0/ C AND REMOVE THE C FROM COLUMN 1 IN ONE OF THE SECTIONS BELOW. C CONSTANTS FOR EVEN OLDER MACHINES CAN BE OBTAINED BY C mail netlib@research.bell-labs.com C send old1mach from blas C PLEASE SEND CORRECTIONS TO dmg OR ehg@bell-labs.com. C C MACHINE CONSTANTS FOR THE HONEYWELL DPS 8/70 SERIES. C DATA RMACH(1) / O402400000000 / C DATA RMACH(2) / O376777777777 / C DATA RMACH(3) / O714400000000 / C DATA RMACH(4) / O716400000000 / C DATA RMACH(5) / O776464202324 /, SC/987/ C C MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING C 32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C DATA SMALL(1) / 8388608 / C DATA LARGE(1) / 2147483647 / C DATA RIGHT(1) / 880803840 / C DATA DIVER(1) / 889192448 / C DATA LOG10(1) / 1067065499 /, SC/987/ C DATA RMACH(1) / O00040000000 / C DATA RMACH(2) / O17777777777 / C DATA RMACH(3) / O06440000000 / C DATA RMACH(4) / O06500000000 / C DATA RMACH(5) / O07746420233 /, SC/987/ C C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. C DATA RMACH(1) / O000400000000 / C DATA RMACH(2) / O377777777777 / C DATA RMACH(3) / O146400000000 / C DATA RMACH(4) / O147400000000 / C DATA RMACH(5) / O177464202324 /, SC/987/ C IF (SC .NE. 987) THEN * *** CHECK FOR AUTODOUBLE *** SMALL(2) = 0 RMACH(1) = 1E13 IF (SMALL(2) .NE. 0) THEN * *** AUTODOUBLED *** IF ( SMALL(1) .EQ. 1117925532 * .AND. SMALL(2) .EQ. -448790528) THEN * *** IEEE BIG ENDIAN *** SMALL(1) = 1048576 SMALL(2) = 0 LARGE(1) = 2146435071 LARGE(2) = -1 RIGHT(1) = 1017118720 RIGHT(2) = 0 DIVER(1) = 1018167296 DIVER(2) = 0 LOG10(1) = 1070810131 LOG10(2) = 1352628735 ELSE IF ( SMALL(2) .EQ. 1117925532 * .AND. SMALL(1) .EQ. -448790528) THEN * *** IEEE LITTLE ENDIAN *** SMALL(2) = 1048576 SMALL(1) = 0 LARGE(2) = 2146435071 LARGE(1) = -1 RIGHT(2) = 1017118720 RIGHT(1) = 0 DIVER(2) = 1018167296 DIVER(1) = 0 LOG10(2) = 1070810131 LOG10(1) = 1352628735 ELSE IF ( SMALL(1) .EQ. -2065213935 * .AND. SMALL(2) .EQ. 10752) THEN * *** VAX WITH D_FLOATING *** SMALL(1) = 128 SMALL(2) = 0 LARGE(1) = -32769 LARGE(2) = -1 RIGHT(1) = 9344 RIGHT(2) = 0 DIVER(1) = 9472 DIVER(2) = 0 LOG10(1) = 546979738 LOG10(2) = -805796613 ELSE IF ( SMALL(1) .EQ. 1267827943 * .AND. SMALL(2) .EQ. 704643072) THEN * *** IBM MAINFRAME *** SMALL(1) = 1048576 SMALL(2) = 0 LARGE(1) = 2147483647 LARGE(2) = -1 RIGHT(1) = 856686592 RIGHT(2) = 0 DIVER(1) = 873463808 DIVER(2) = 0 LOG10(1) = 1091781651 LOG10(2) = 1352628735 ELSE WRITE(*,9010) STOP 777 END IF ELSE RMACH(1) = 1234567. IF (SMALL(1) .EQ. 1234613304) THEN * *** IEEE *** SMALL(1) = 8388608 LARGE(1) = 2139095039 RIGHT(1) = 864026624 DIVER(1) = 872415232 LOG10(1) = 1050288283 ELSE IF (SMALL(1) .EQ. -1271379306) THEN * *** VAX *** SMALL(1) = 128 LARGE(1) = -32769 RIGHT(1) = 13440 DIVER(1) = 13568 LOG10(1) = 547045274 ELSE IF (SMALL(1) .EQ. 1175639687) THEN * *** IBM MAINFRAME *** SMALL(1) = 1048576 LARGE(1) = 2147483647 RIGHT(1) = 990904320 DIVER(1) = 1007681536 LOG10(1) = 1091781651 ELSE IF (SMALL(1) .EQ. 1251390520) THEN * *** CONVEX C-1 *** SMALL(1) = 8388608 LARGE(1) = 2147483647 RIGHT(1) = 880803840 DIVER(1) = 889192448 LOG10(1) = 1067065499 ELSE DO 10 L = 1, 3 J = SMALL(1) / 10000000 K = SMALL(1) - 10000000*J IF (K .NE. T3E(L)) GO TO 20 SMALL(1) = J 10 CONTINUE * *** CRAY T3E *** CALL I1MCRA(SMALL, K, 16, 0, 0) CALL I1MCRA(LARGE, K, 32751, 16777215, 16777215) CALL I1MCRA(RIGHT, K, 15520, 0, 0) CALL I1MCRA(DIVER, K, 15536, 0, 0) CALL I1MCRA(LOG10, K, 16339, 4461392, 10451455) GO TO 30 20 CALL I1MCRA(J, K, 16405, 9876536, 0) IF (SMALL(1) .NE. J) THEN WRITE(*,9020) STOP 777 END IF * *** CRAY 1, XMP, 2, AND 3 *** CALL I1MCRA(SMALL(1), K, 8195, 8388608, 1) CALL I1MCRA(LARGE(1), K, 24574, 16777215, 16777214) CALL I1MCRA(RIGHT(1), K, 16338, 8388608, 0) CALL I1MCRA(DIVER(1), K, 16339, 8388608, 0) CALL I1MCRA(LOG10(1), K, 16383, 10100890, 8715216) END IF END IF 30 SC = 987 END IF * SANITY CHECK IF (RMACH(4) .GE. 1.0) STOP 776 IF (I .LT. 1 .OR. I .GT. 5) THEN WRITE(*,*) 'R1MACH(I): I =',I,' is out of bounds.' STOP END IF R1MACH = RMACH(I) RETURN 9010 FORMAT(/' Adjust autodoubled R1MACH by getting data'/ *' appropriate for your machine from D1MACH.') 9020 FORMAT(/' Adjust R1MACH by uncommenting data statements'/ *' appropriate for your machine.') * /* C source for R1MACH -- remove the * in column 1 */ *#include *#include *#include *float r1mach_(long *i) *{ * switch(*i){ * case 1: return FLT_MIN; * case 2: return FLT_MAX; * case 3: return FLT_EPSILON/FLT_RADIX; * case 4: return FLT_EPSILON; * case 5: return log10((double)FLT_RADIX); * } * fprintf(stderr, "invalid argument: r1mach(%ld)\n", *i); * exit(1); return 0; /* else complaint of missing return value */ *} END SUBROUTINE I1MCRA(A, A1, B, C, D) **** SPECIAL COMPUTATION FOR CRAY MACHINES **** INTEGER A, A1, B, C, D A1 = 16777216*B + C A = 16777216*A1 + D END subroutine xerror(messg,nmessg,nerr,level) c c abstract c xerror processes a diagnostic message, in a manner c determined by the value of level and the current value c of the library error control flag, kontrl. c (see subroutine xsetf for details.) c c description of parameters c --input-- c messg - the hollerith message to be processed, containing c no more than 72 characters. c nmessg- the actual number of characters in messg. c nerr - the error number associated with this message. c nerr must not be zero. c level - error category. c =2 means this is an unconditionally fatal error. c =1 means this is a recoverable error. (i.e., it is c non-fatal if xsetf has been appropriately called.) c =0 means this is a warning message only. c =-1 means this is a warning message which is to be c printed at most once, regardless of how many c times this call is executed. c c examples c call xerror(23hsmooth -- num was zero.,23,1,2) c call xerror(43hinteg -- less than full accuracy achieved., c 43,2,1) c call xerror(65hrooter -- actual zero of f found before interval c 1 fully collapsed.,65,3,0) c call xerror(39hexp -- underflows being set to zero.,39,1,-1) c c written by ron jones, with slatec common math library subcommittee c latest revision --- 7 feb 1979 c dimension messg(nmessg) call xerrwv(messg,nmessg,nerr,level,0,0,0,0,0.,0.) return end subroutine xerrwv (msg, nmes, nerr, level, ni, i1, i2, nr, r1, r2) integer msg, nmes, nerr, level, ni, i1, i2, nr, 1 i, lun, lunit, mesflg, ncpw, nch, nwds real r1, r2 dimension msg(nmes) c----------------------------------------------------------------------- c subroutines xerrwv, xsetf, and xsetun, as given here, constitute c a simplified version of the slatec error handling package. c written by a. c. hindmarsh at llnl. version of march 30, 1987. c c all arguments are input arguments. c c msg = the message (hollerith literal or integer array). c nmes = the length of msg (number of characters). c nerr = the error number (not used). c level = the error level.. c 0 or 1 means recoverable (control returns to caller). c 2 means fatal (run is aborted--see note below). c ni = number of integers (0, 1, or 2) to be printed with message. c i1,i2 = integers to be printed, depending on ni. c nr = number of reals (0, 1, or 2) to be printed with message. c r1,r2 = reals to be printed, depending on nr. c c note.. this routine is machine-dependent and specialized for use c in limited context, in the following ways.. c 1. the number of hollerith characters stored per word, denoted c by ncpw below, is a data-loaded constant. c 2. the value of nmes is assumed to be at most 60. c (multi-line messages are generated by repeated calls.) c 3. if level = 2, control passes to the statement stop c to abort the run. this statement may be machine-dependent. c 4. r1 and r2 are assumed to be in single precision and are printed c in e21.13 format. c 5. the common block /eh0001/ below is data-loaded (a machine- c dependent feature) with default values. c this block is needed for proper retention of parameters used by c this routine which the user can reset by calling xsetf or xsetun. c the variables in this block are as follows.. c mesflg = print control flag.. c 1 means print all messages (the default). c 0 means no printing. c lunit = logical unit number for messages. c the default is 6 (machine-dependent). c----------------------------------------------------------------------- c the following are instructions for installing this routine c in different machine environments. c c to change the default output unit, change the data statement below. c c for some systems, the data statement below must be replaced c by a separate block data subprogram. c c for a different number of characters per word, change the c data statement setting ncpw below, and format 10. alternatives for c various computers are shown in comment cards. c c for a different run-abort command, change the statement following c statement 100 at the end. c----------------------------------------------------------------------- common /eh0001/ mesflg, lunit c data mesflg/1/, lunit/6/ c----------------------------------------------------------------------- c the following data-loaded value of ncpw is valid for the cdc-6600 c and cdc-7600 computers. c data ncpw/10/ c the following is valid for the cray-1 computer. data ncpw/8/ c the following is valid for the burroughs 6700 and 7800 computers. c data ncpw/6/ c the following is valid for the pdp-10 computer. c data ncpw/5/ c the following is valid for the vax computer with 4 bytes per integer, c and for the ibm-360, ibm-370, ibm-303x, and ibm-43xx computers. c data ncpw/4/ c the following is valid for the pdp-11, or vax with 2-byte integers. c data ncpw/2/ c----------------------------------------------------------------------- c if (mesflg .eq. 0) go to 100 c get logical unit number. --------------------------------------------- lun = lunit c get number of words in message. -------------------------------------- nch = min0(nmes,60) nwds = nch/ncpw if (nch .ne. nwds*ncpw) nwds = nwds + 1 c write the message. --------------------------------------------------- write (lun, 10) (msg(i),i=1,nwds) c----------------------------------------------------------------------- c the following format statement is to have the form c 10 format(1x,mmann) c where nn = ncpw and mm is the smallest integer .ge. 60/ncpw. c the following is valid for ncpw = 10. c 10 format(1x,6a10) c the following is valid for ncpw = 8. 10 format(1x,8a8) c the following is valid for ncpw = 6. c 10 format(1x,10a6) c the following is valid for ncpw = 5. c 10 format(1x,12a5) c the following is valid for ncpw = 4. c 10 format(1x,15a4) c the following is valid for ncpw = 2. c 10 format(1x,30a2) c----------------------------------------------------------------------- if (ni .eq. 1) write (lun, 20) i1 20 format(6x,23hin above message, i1 =,i10) if (ni .eq. 2) write (lun, 30) i1,i2 30 format(6x,23hin above message, i1 =,i10,3x,4hi2 =,i10) if (nr .eq. 1) write (lun, 40) r1 40 format(6x,23hin above message, r1 =,e21.13) if (nr .eq. 2) write (lun, 50) r1,r2 50 format(6x,15hin above, r1 =,e21.13,3x,4hr2 =,e21.13) c abort the run if level = 2. ------------------------------------------ 100 if (level .ne. 2) return stop c----------------------- end of subroutine xerrwv ---------------------- end