subroutine qagp(f,a,b,npts2,points,epsabs,epsrel,result,abserr,
     *   neval,ier,leniw,lenw,last,iwork,work)
c***begin prologue  qagp
c***date written   800101   (yymmdd)
c***revision date  830518   (yymmdd)
c***category no.  h2a2a1
c***keywords  automatic integrator, general-purpose,
c             singularities at user specified points,
c             extrapolation, 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            break points of the integration interval, where local
c            difficulties of the integrand may occur(e.g. singularities,
c            discontinuities), are provided by the user.
c***description
c
c        computation of a definite integral
c        standard fortran subroutine
c        real version
c
c        parameters
c         on entry
c            f      - subroutine f(x,ierr,result) 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            npts2  - integer
c                     number equal to two more than the number of
c                     user-supplied break points within the integration
c                     range, npts.ge.2.
c                     if npts2.lt.2, the routine will end with ier = 6.
c
c            points - real
c                     vector of dimension npts2, the first (npts2-2)
c                     elements of which are the user provided break
c                     points. if these points do not constitute an
c                     ascending sequence there will be an automatic
c                     sorting.
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         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                     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. if
c                             the position of a local difficulty can be
c                             determined (i.e. singularity,
c                             discontinuity within the interval), it
c                             should be supplied to the routine as an
c                             element of the vector points. if necessary
c                             an appropriate special-purpose integrator
c                             must 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 presumed that the requested
c                             tolerance cannot be achieved, and that
c                             the returned result is the best which
c                             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.gt.0.
c                         = 6 the input is invalid because
c                             npts2.lt.2 or
c                             break points are specified outside
c                             the integration range or
c                             (epsabs.le.0 and
c                              epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
c                             result, abserr, neval, last are set to
c                             zero. exept when leniw or lenw or npts2 is
c                             invalid, iwork(1), iwork(limit+1),
c                             work(limit*2+1) and work(limit*3+1)
c                             are set to zero.
c                             work(1) is set to a and work(limit+1)
c                             to b (where limit = (leniw-npts2)/2).
c
c         dimensioning parameters
c            leniw - integer
c                    dimensioning parameter for iwork
c                    leniw determines limit = (leniw-npts2)/2,
c                    which is the maximum number of subintervals in the
c                    partition of the given integration interval (a,b),
c                    leniw.ge.(3*npts2-2).
c                    if leniw.lt.(3*npts2-2), the routine will end with
c                    ier = 6.
c
c            lenw  - integer
c                    dimensioning parameter for work
c                    lenw must be at least leniw*2-npts2.
c                    if lenw.lt.leniw*2-npts2, the routine will end
c                    with 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 leniw. on return,
c                    the first k elements of which contain
c                    pointers to the error estimates over the
c                    subintervals, such that work(limit*3+iwork(1)),...,
c                    work(limit*3+iwork(k)) form a decreasing
c                    sequence, with k = last if last.le.(limit/2+2), and
c                    k = limit+1-last otherwise
c                    iwork(limit+1), ...,iwork(limit+last) contain the
c                     subdivision levels of the subintervals, i.e.
c                     if (aa,bb) is a subinterval of (p1,p2)
c                     where p1 as well as p2 is a user-provided
c                     break point or integration limit, then (aa,bb) has
c                     level l if abs(bb-aa) = abs(p2-p1)*2**(-l),
c                    iwork(limit*2+1), ..., iwork(limit*2+npts2) have
c                     no significance for the user,
c                    note that limit = (leniw-npts2)/2.
c
c            work  - real
c                    vector of dimension at least lenw
c                    on return
c                    work(1), ..., work(last) contain the left
c                     end points of the subintervals in the
c                     partition of (a,b),
c                    work(limit+1), ..., work(limit+last) contain
c                     the 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)
c                     contain the corresponding error estimates,
c                    work(limit*4+1), ..., work(limit*4+npts2)
c                     contain the integration limits and the
c                     break points sorted in an ascending sequence.
c                    note that limit = (leniw-npts2)/2.
c
c***references  (none)
c***routines called  qagpe,xerror
c***end prologue  qagp
c
      real a,abserr,b,epsabs,epsrel,points,result,work
      integer ier,iwork,leniw,lenw,limit,lvl,l1,l2,l3,neval,npts2
c
      dimension iwork(leniw),points(npts2),work(lenw)
c
      external f
c
c         check validity of limit and lenw.
c
c***first executable statement  qagp
      ier = 6
      neval = 0
      last = 0
      result = 0.0e+00
      abserr = 0.0e+00
      if(leniw.lt.(3*npts2-2).or.lenw.lt.(leniw*2-npts2).or.npts2.lt.2)
     *  go to 10
c
c         prepare call for qagpe.
c
      limit = (leniw-npts2)/2
      l1 = limit+1
      l2 = limit+l1
      l3 = limit+l2
      l4 = limit+l3
c
      call qagpe(f,a,b,npts2,points,epsabs,epsrel,limit,result,abserr,
     *  neval,ier,work(1),work(l1),work(l2),work(l3),work(l4),
     *  iwork(1),iwork(l1),iwork(l2),last)
c
c         call error handler if necessary.
c
      lvl = 0
10    if(ier.eq.6) lvl = 1
      if(ier.ne.0) call xerror('abnormal return from  qagp',26,ier,lvl)
      return
      end