subroutine qagse(f,a,b,epsabs,epsrel,limit,result,abserr,neval, * ier,alist,blist,rlist,elist,iord,last) c***begin prologue qagse c***date written 800101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a1a1 c***keywords automatic integrator, general-purpose, c (end point) singularities, extrapolation, c globally adaptive 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 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 limit - integer c gives an upperbound on the number of subintervals c in the partition of (a,b) 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 integral and error are c less reliable. it is assumed that the c requested accuracy has not been achieved. c error messages c = 1 maximum number of subdivisions allowed c has been achieved. one can allow more sub- c divisions by increasing the value of limit c (and taking the according dimension c adjustments into account). however, if c this yields no improvement it is advised c to analyze the integrand in order to c determine the integration difficulties. if c the position of a local difficulty can be c determined (e.g. 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 detec- c ted, which prevents the requested c tolerance from being achieved. c the error may be under-estimated. c = 3 extremely bad integrand behaviour c occurs at some points of the integration c interval. c = 4 the algorithm does not converge. c roundoff error is detected in the c extrapolation table. c it is presumed that the requested c tolerance cannot be achieved, and that the c returned result is the best which can be c obtained. c = 5 the integral is probably divergent, or c slowly convergent. it must be noted that c divergence can occur with any other value c of ier. 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 iord(1) and elist(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 end points c of the subintervals in the partition of the c 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 end points c of the subintervals in the partition of the given c integration range (a,b) c c rlist - real c vector of dimension at least limit, the first c last elements of which are the integral c 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)), ..., elist(iord(k)) c form a decreasing sequence, with k = last c if last.le.(limit/2+2), and k = limit+1-last c otherwise c c last - integer c number of subintervals actually produced in the c subdivision process c c***references (none) c***routines called qelg,qk21,qpsrt,r1mach c***end prologue qagse c real a,abseps,abserr,alist,area,area1,area12,area2,a1, * a2,b,blist,b1,b2,correc,defabs,defab1,defab2,r1mach, * dres,elist,epmach,epsabs,epsrel,erlarg,erlast,errbnd, * errmax,error1,error2,erro12,errsum,ertest,f,oflow,resabs, * reseps,result,res3la,rlist,rlist2,small,uflow integer id,ier,ierro,iord,iroff1,iroff2,iroff3,jupbnd,k,ksgn, * ktmin,last,limit,maxerr,neval,nres,nrmax,numrl2 logical extrap,noext c dimension alist(limit),blist(limit),elist(limit),iord(limit), * res3la(3),rlist(limit),rlist2(52) c external f c c the dimension of rlist2 is determined by the value of c limexp in subroutine qelg (rlist2 should be of dimension c (limexp+2) at least). 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 rlist2 - array of dimension at least limexp+2 c containing the part of the epsilon table c which is still needed for further computations c elist(i) - error estimate applying to rlist(i) c maxerr - pointer to the interval with largest error c estimate c errmax - elist(maxerr) c erlast - error on the interval currently subdivided c (before that subdivision has taken place) 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 interval c *****2 - variable for the right interval c last - index for subdivision c nres - number of calls to the extrapolation routine c numrl2 - number of elements currently in rlist2. if an c appropriate approximation to the compounded c integral has been obtained it is put in c rlist2(numrl2) after numrl2 has been increased c by one. c small - length of the smallest interval considered c up to now, multiplied by 1.5 c erlarg - sum of the errors over the intervals larger c than the smallest interval considered up to now c extrap - logical variable denoting that the routine c is attempting to perform extrapolation c i.e. before subdividing the smallest interval c we try to decrease the value of erlarg. c noext - logical variable denoting that extrapolation c is no longer allowed (true value) c c machine dependent constants c --------------------------- c c epmach is the largest relative spacing. c uflow is the smallest positive magnitude. c oflow is the largest positive magnitude. c c***first executable statement qagse epmach = r1mach(4) c c test on validity of parameters 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 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 uflow = r1mach(1) oflow = r1mach(2) ierro = 0 call qk21(f,a,b,result,abserr,defabs,resabs) c c test on accuracy. c dres = abs(result) errbnd = amax1(epsabs,epsrel*dres) last = 1 rlist(1) = result elist(1) = abserr iord(1) = 1 if(abserr.le.1.0e+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 140 c c initialization c -------------- c rlist2(1) = result errmax = abserr maxerr = 1 area = result errsum = abserr abserr = oflow nrmax = 1 nres = 0 numrl2 = 2 ktmin = 0 extrap = .false. noext = .false. iroff1 = 0 iroff2 = 0 iroff3 = 0 ksgn = -1 if(dres.ge.(0.1e+01-0.5e+02*epmach)*defabs) ksgn = 1 c c main do-loop c ------------ c do 90 last = 2,limit c c bisect the subinterval with the nrmax-th largest c error estimate. c a1 = alist(maxerr) b1 = 0.5e+00*(alist(maxerr)+blist(maxerr)) a2 = b1 b2 = blist(maxerr) erlast = errmax call qk21(f,a1,b1,area1,error1,resabs,defab1) call qk21(f,a2,b2,area2,error2,resabs,defab2) c c improve previous approximations to integral c and error and test for accuracy. c 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 15 if(abs(rlist(maxerr)-area12).gt.0.1e-04*abs(area12) * .or.erro12.lt.0.99e+00*errmax) go to 10 if(extrap) iroff2 = iroff2+1 if(.not.extrap) iroff1 = iroff1+1 10 if(last.gt.10.and.erro12.gt.errmax) iroff3 = iroff3+1 15 rlist(maxerr) = area1 rlist(last) = area2 errbnd = amax1(epsabs,epsrel*abs(area)) c c test for roundoff error and eventually c set error flag. c if(iroff1+iroff2.ge.10.or.iroff3.ge.20) ier = 2 if(iroff2.ge.5) ierro = 3 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 = 4 c c append the newly-created intervals to the list. c if(error2.gt.error1) go to 20 alist(last) = a2 blist(maxerr) = b1 blist(last) = b2 elist(maxerr) = error1 elist(last) = error2 go to 30 20 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 nrmax-th largest error estimate (to be c bisected next). c 30 call qpsrt(limit,last,maxerr,errmax,elist,iord,nrmax) c ***jump out of do-loop if(errsum.le.errbnd) go to 115 c ***jump out of do-loop if(ier.ne.0) go to 100 if(last.eq.2) go to 80 if(noext) go to 90 erlarg = erlarg-erlast if(abs(b1-a1).gt.small) erlarg = erlarg+erro12 if(extrap) go to 40 c c test whether the interval to be bisected next is the c smallest interval. c if(abs(blist(maxerr)-alist(maxerr)).gt.small) go to 90 extrap = .true. nrmax = 2 40 if(ierro.eq.3.or.erlarg.le.ertest) go to 60 c c the smallest interval has the largest error. c before bisecting decrease the sum of the errors c over the larger intervals (erlarg) and perform c extrapolation. c id = nrmax jupbnd = last if(last.gt.(2+limit/2)) jupbnd = limit+3-last do 50 k = id,jupbnd maxerr = iord(nrmax) errmax = elist(maxerr) c ***jump out of do-loop if(abs(blist(maxerr)-alist(maxerr)).gt.small) go to 90 nrmax = nrmax+1 50 continue c c perform extrapolation. c 60 numrl2 = numrl2+1 rlist2(numrl2) = area call qelg(numrl2,rlist2,reseps,abseps,res3la,nres) ktmin = ktmin+1 if(ktmin.gt.5.and.abserr.lt.0.1e-02*errsum) ier = 5 if(abseps.ge.abserr) go to 70 ktmin = 0 abserr = abseps result = reseps correc = erlarg ertest = amax1(epsabs,epsrel*abs(reseps)) c ***jump out of do-loop if(abserr.le.ertest) go to 100 c c prepare bisection of the smallest interval. c 70 if(numrl2.eq.1) noext = .true. if(ier.eq.5) go to 100 maxerr = iord(1) errmax = elist(maxerr) nrmax = 1 extrap = .false. small = small*0.5e+00 erlarg = errsum go to 90 80 small = abs(b-a)*0.375e+00 erlarg = errsum ertest = errbnd rlist2(2) = area 90 continue c c set final result and error estimate. c ------------------------------------ c 100 if(abserr.eq.oflow) go to 115 if(ier+ierro.eq.0) go to 110 if(ierro.eq.3) abserr = abserr+correc if(ier.eq.0) ier = 3 if(result.ne.0.0e+00.and.area.ne.0.0e+00) go to 105 if(abserr.gt.errsum) go to 115 if(area.eq.0.0e+00) go to 130 go to 110 105 if(abserr/abs(result).gt.errsum/abs(area)) go to 115 c c test on divergence. c 110 if(ksgn.eq.(-1).and.amax1(abs(result),abs(area)).le. * defabs*0.1e-01) go to 130 if(0.1e-01.gt.(result/area).or.(result/area).gt.0.1e+03 * .or.errsum.gt.abs(area)) ier = 6 go to 130 c c compute global integral sum. c 115 result = 0.0e+00 do 120 k = 1,last result = result+rlist(k) 120 continue abserr = errsum 130 if(ier.gt.2) ier = ier-1 140 neval = 42*last-21 999 return end