subroutine qagie(f,bound,inf,epsabs,epsrel,limit,result,abserr,
* neval,ier,alist,blist,rlist,elist,iord,last)
c***begin prologue qagie
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 - 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 bound - real
c finite bound of integration range
c (has no meaning if interval is doubly-infinite)
c
c inf - real
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 - 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 upper bound 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. 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 - 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 transformed integration range (0,1).
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 transformed integration range (0,1).
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 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 qelg,qk15i,qpsrt,r1mach
c***end prologue qagie
c
real abseps,abserr,alist,area,area1,area12,area2,a1,
* a2,blist,boun,bound,b1,b2,correc,defabs,defab1,defab2,
* dres,r1mach,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 qelg.
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
epmach = r1mach(4)
c
c test on validity of parameters
c -----------------------------
c
c***first executable statement qagie
ier = 0
neval = 0
last = 0
result = 0.0e+00
abserr = 0.0e+00
alist(1) = 0.0e+00
blist(1) = 0.1e+01
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
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.0e+00
call qk15i(f,boun,inf,0.0e+00,0.1e+01,result,abserr,
* defabs,resabs)
c
c test on accuracy
c
last = 1
rlist(1) = result
elist(1) = abserr
iord(1) = 1
dres = abs(result)
errbnd = amax1(epsabs,epsrel*dres)
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 130
c
c initialization
c --------------
c
uflow = r1mach(1)
oflow = r1mach(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.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 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 qk15i(f,boun,inf,a1,b1,area1,error1,resabs,defab1)
call qk15i(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(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 some points 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)
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(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)
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))
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 = 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 neval = 30*last-15
if(inf.eq.2) neval = 2*neval
if(ier.gt.2) ier=ier-1
999 return
end
.