subroutine dqagie(f,bound,inf,epsabs,epsrel,limit,result,abserr, * neval,ier,alist,blist,rlist,elist,iord,last) c***begin prologue dqagie c***date written 800101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a3a1,h2a4a1 c***keywords automatic integrator, infinite intervals, c general-purpose, transformation, 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 integral i = integral of f over (bound,+infinity) c or i = integral of f over (-infinity,bound) c or i = integral of f over (-infinity,+infinity), c hopefully satisfying following claim for accuracy c abs(i-result).le.max(epsabs,epsrel*abs(i)) c***description c c integration over infinite intervals c standard fortran subroutine c 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 driver program. c c bound - double precision c finite bound of integration range c (has no meaning if interval is doubly-infinite) c c inf - double precision c indicating the kind of integration range involved c inf = 1 corresponds to (bound,+infinity), c inf = -1 to (-infinity,bound), c inf = 2 to (-infinity,+infinity). c c epsabs - double precision c absolute accuracy requested c epsrel - double precision 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 upper bound on the number of subintervals c in the partition of (a,b), limit.ge.1 c c on return c result - double precision c approximation to the integral c c abserr - double precision 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. the c estimates for result and error are less c reliable. it is assumed that the requested c 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 yields no improvement it is advised c to analyze the integrand in order to c determine the integration difficulties. c if the position of a local difficulty can c be 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 c detected, which prevents the requested c tolerance from being achieved. c the error may be under-estimated. c = 3 extremely bad integrand behaviour occurs c 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 assumed that the requested tolerance c cannot be achieved, and that the returned c result is the best which can be 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 elist(1) and iord(1) are set to zero. c alist(1) and blist(1) are set to 0 c and 1 respectively. c c alist - double precision 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 transformed integration range (0,1). c c blist - double precision 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 transformed integration range (0,1). c c rlist - double precision 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 - double precision 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 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 c in the subdivision process c c***references (none) c***routines called d1mach,dqelg,dqk15i,dqpsrt c***end prologue dqagie double precision abseps,abserr,alist,area,area1,area12,area2,a1, * a2,blist,boun,bound,b1,b2,correc,dabs,defabs,defab1,defab2, * dmax1,dres,d1mach,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,inf,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 dqelg. c 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 wich 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 subinterval c *****2 - variable for the right subinterval 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 up c 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. i.e. c before subdividing the smallest interval we c 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 dqagie epmach = d1mach(4) c c test on validity of parameters c ----------------------------- c ier = 0 neval = 0 last = 0 result = 0.0d+00 abserr = 0.0d+00 alist(1) = 0.0d+00 blist(1) = 0.1d+01 rlist(1) = 0.0d+00 elist(1) = 0.0d+00 iord(1) = 0 if(epsabs.le.0.0d+00.and.epsrel.lt.dmax1(0.5d+02*epmach,0.5d-28)) * ier = 6 if(ier.eq.6) go to 999 c c c first approximation to the integral c ----------------------------------- c c determine the interval to be mapped onto (0,1). c if inf = 2 the integral is computed as i = i1+i2, where c i1 = integral of f over (-infinity,0), c i2 = integral of f over (0,+infinity). c boun = bound if(inf.eq.2) boun = 0.0d+00 call dqk15i(f,boun,inf,0.0d+00,0.1d+01,result,abserr, * defabs,resabs) c c test on accuracy c last = 1 rlist(1) = result elist(1) = abserr iord(1) = 1 dres = dabs(result) errbnd = dmax1(epsabs,epsrel*dres) if(abserr.le.1.0d+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.0d+00) go to 130 c c initialization c -------------- c uflow = d1mach(1) oflow = d1mach(2) rlist2(1) = result errmax = abserr maxerr = 1 area = result errsum = abserr abserr = oflow nrmax = 1 nres = 0 ktmin = 0 numrl2 = 2 extrap = .false. noext = .false. ierro = 0 iroff1 = 0 iroff2 = 0 iroff3 = 0 ksgn = -1 if(dres.ge.(0.1d+01-0.5d+02*epmach)*defabs) ksgn = 1 c c main do-loop c ------------ c do 90 last = 2,limit c c bisect the subinterval with nrmax-th largest error estimate. c a1 = alist(maxerr) b1 = 0.5d+00*(alist(maxerr)+blist(maxerr)) a2 = b1 b2 = blist(maxerr) erlast = errmax call dqk15i(f,boun,inf,a1,b1,area1,error1,resabs,defab1) call dqk15i(f,boun,inf,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(dabs(rlist(maxerr)-area12).gt.0.1d-04*dabs(area12) * .or.erro12.lt.0.99d+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 = dmax1(epsabs,epsrel*dabs(area)) c c test for roundoff error and eventually 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 some points of the integration range. c if(dmax1(dabs(a1),dabs(b2)).le.(0.1d+01+0.1d+03*epmach)* * (dabs(a2)+0.1d+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 dqpsrt to maintain the descending ordering c in the list of error estimates and select the subinterval c with nrmax-th largest error estimate (to be bisected next). c 30 call dqpsrt(limit,last,maxerr,errmax,elist,iord,nrmax) if(errsum.le.errbnd) go to 115 if(ier.ne.0) go to 100 if(last.eq.2) go to 80 if(noext) go to 90 erlarg = erlarg-erlast if(dabs(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(dabs(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 over the c larger intervals (erlarg) and perform 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) if(dabs(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 dqelg(numrl2,rlist2,reseps,abseps,res3la,nres) ktmin = ktmin+1 if(ktmin.gt.5.and.abserr.lt.0.1d-02*errsum) ier = 5 if(abseps.ge.abserr) go to 70 ktmin = 0 abserr = abseps result = reseps correc = erlarg ertest = dmax1(epsabs,epsrel*dabs(reseps)) 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.5d+00 erlarg = errsum go to 90 80 small = 0.375d+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.0d+00.and.area.ne.0.0d+00)go to 105 if(abserr.gt.errsum)go to 115 if(area.eq.0.0d+00) go to 130 go to 110 105 if(abserr/dabs(result).gt.errsum/dabs(area))go to 115 c c test on divergence c 110 if(ksgn.eq.(-1).and.dmax1(dabs(result),dabs(area)).le. * defabs*0.1d-01) go to 130 if(0.1d-01.gt.(result/area).or.(result/area).gt.0.1d+03. *or.errsum.gt.dabs(area)) ier = 6 go to 130 c c compute global integral sum. c 115 result = 0.0d+00 do 120 k = 1,last result = result+rlist(k) 120 continue abserr = errsum 130 neval = 30*last-15 if(inf.eq.2) neval = 2*neval if(ier.gt.2) ier=ier-1 999 return end