# qpdoc.f

          SUBROUTINE QPDOC
    C***BEGIN PROLOGUE  QPDOC
    C***PURPOSE  Documentation for QUADPACK, a package of subprograms for
    C            automatic evaluation of one-dimensional definite integrals.
    C***LIBRARY   SLATEC (QUADPACK)
    C***CATEGORY  H2, Z
    C***TYPE      ALL (QPDOC-A)
    C***KEYWORDS  DOCUMENTATION, GUIDELINES FOR SELECTION, QUADPACK,
    C             QUADRATURE, SURVEY OF INTEGRATORS
    C***AUTHOR  Piessens, Robert
    C             Applied Mathematics and Programming Division
    C             K. U. Leuven
    C           de Doncker, Elise
    C             Applied Mathematics and Programming Division
    C             K. U. Leuven
    C           Kahaner, D. K., (NBS)
    C***DESCRIPTION
    C
    C 1. Introduction
    C    ------------
    C    QUADPACK is a FORTRAN subroutine package for the numerical
    C    computation of definite one-dimensional integrals. It originated
    C    from a joint project of R. Piessens and E. de Doncker (Appl.
    C    Math. and Progr. Div.- K.U.Leuven, Belgium), C. Ueberhuber (Inst.
    C    Fuer Math.- Techn. U. Wien, Austria), and D. Kahaner (National
    C    Bureau of Standards- Washington D.C., U.S.A.).
    C
    C    Documentation routine QPDOC describes the package in the form it
    C    was released from A.M.P.D.- Leuven, for adherence to the SLATEC
    C    library in May 1981. Apart from a survey of the integrators, some
    C    guidelines will be given in order to help the QUADPACK user with
    C    selecting an appropriate routine or a combination of several
    C    routines for handling his problem.
    C
    C    In the Long Description of QPDOC it is demonstrated how to call
    C    the integrators, by means of small example calling programs.
    C
    C    For precise guidelines involving the use of each routine in
    C    particular, we refer to the extensive introductory comments
    C    within each routine.
    C
    C 2. Survey
    C    ------
    C    The following list gives an overview of the QUADPACK integrators.
    C    The routine names for the DOUBLE PRECISION versions are preceded
    C    by the letter D.
    C
    C    - QNG  : Is a simple non-adaptive automatic integrator, based on
    C             a sequence of rules with increasing degree of algebraic
    C             precision (Patterson, 1968).
    C
    C    - QAG  : Is a simple globally adaptive integrator using the
    C             strategy of Aind (Piessens, 1973). It is possible to
    C             choose between 6 pairs of Gauss-Kronrod quadrature
    C             formulae for the rule evaluation component. The pairs
    C             of high degree of precision are suitable for handling
    C             integration difficulties due to a strongly oscillating
    C             integrand.
    C
    C    - QAGS : Is an integrator based on globally adaptive interval
    C             subdivision in connection with extrapolation (de Doncker,
    C             1978) by the Epsilon algorithm (Wynn, 1956).
    C
    C    - QAGP : Serves the same purposes as QAGS, but also allows
    C             for eventual user-supplied information, i.e. the
    C             abscissae of internal singularities, discontinuities
    C             and other difficulties of the integrand function.
    C             The algorithm is a modification of that in QAGS.
    C
    C    - QAGI : Handles integration over infinite intervals. The
    C             infinite range is mapped onto a finite interval and
    C             then the same strategy as in QAGS is applied.
    C
    C    - QAWO : Is a routine for the integration of COS(OMEGA*X)*F(X)
    C             or SIN(OMEGA*X)*F(X) over a finite interval (A,B).
    C             OMEGA is is specified by the user
    C             The rule evaluation component is based on the
    C             modified Clenshaw-Curtis technique.
    C             An adaptive subdivision scheme is used connected with
    C             an extrapolation procedure, which is a modification
    C             of that in QAGS and provides the possibility to deal
    C             even with singularities in F.
    C
    C    - QAWF : Calculates the Fourier cosine or Fourier sine
    C             transform of F(X), for user-supplied interval (A,
    C             INFINITY), OMEGA, and F. The procedure of QAWO is
    C             used on successive finite intervals, and convergence
    C             acceleration by means of the Epsilon algorithm (Wynn,
    C             1956) is applied to the series of the integral
    C             contributions.
    C
    C    - QAWS : Integrates W(X)*F(X) over (A,B) with A.LT.B finite,
    C             and   W(X) = ((X-A)**ALFA)*((B-X)**BETA)*V(X)
    C             where V(X) = 1 or LOG(X-A) or LOG(B-X)
    C                            or LOG(X-A)*LOG(B-X)
    C             and   ALFA.GT.(-1), BETA.GT.(-1).
    C             The user specifies A, B, ALFA, BETA and the type of
    C             the function V.
    C             A globally adaptive subdivision strategy is applied,
    C             with modified Clenshaw-Curtis integration on the
    C             subintervals which contain A or B.
    C
    C    - QAWC : Computes the Cauchy Principal Value of F(X)/(X-C)
    C             over a finite interval (A,B) and for
    C             user-determined C.
    C             The strategy is globally adaptive, and modified
    C             Clenshaw-Curtis integration is used on the subranges
    C             which contain the point X = C.
    C
    C  Each of the routines above also has a "more detailed" version
    C    with a name ending in E, as QAGE.  These provide more
    C    information and control than the easier versions.
    C
    C
    C   The preceding routines are all automatic.  That is, the user
    C      inputs his problem and an error tolerance.  The routine
    C      attempts to perform the integration to within the requested
    C      absolute or relative error.
    C   There are, in addition, a number of non-automatic integrators.
    C      These are most useful when the problem is such that the
    C      user knows that a fixed rule will provide the accuracy
    C      required.  Typically they return an error estimate but make
    C      no attempt to satisfy any particular input error request.
    C
    C      QK15
    C      QK21
    C      QK31
    C      QK41
    C      QK51
    C      QK61
    C           Estimate the integral on [a,b] using 15, 21,..., 61
    C           point rule and return an error estimate.
    C      QK15I 15 point rule for (semi)infinite interval.
    C      QK15W 15 point rule for special singular weight functions.
    C      QC25C 25 point rule for Cauchy Principal Values
    C      QC25F 25 point rule for sin/cos integrand.
    C      QMOMO Integrates k-th degree Chebyshev polynomial times
    C            function with various explicit singularities.
    C
    C 3. Guidelines for the use of QUADPACK
    C    ----------------------------------
    C    Here it is not our purpose to investigate the question when
    C    automatic quadrature should be used. We shall rather attempt
    C    to help the user who already made the decision to use QUADPACK,
    C    with selecting an appropriate routine or a combination of
    C    several routines for handling his problem.
    C
    C    For both quadrature over finite and over infinite intervals,
    C    one of the first questions to be answered by the user is
    C    related to the amount of computer time he wants to spend,
    C    versus his -own- time which would be needed, for example, for
    C    manual subdivision of the interval or other analytic
    C    manipulations.
    C
    C    (1) The user may not care about computer time, or not be
    C        willing to do any analysis of the problem. especially when
    C        only one or a few integrals must be calculated, this attitude
    C        can be perfectly reasonable. In this case it is clear that
    C        either the most sophisticated of the routines for finite
    C        intervals, QAGS, must be used, or its analogue for infinite
    C        intervals, GAGI. These routines are able to cope with
    C        rather difficult, even with improper integrals.
    C        This way of proceeding may be expensive. But the integrator
    C        is supposed to give you an answer in return, with additional
    C        information in the case of a failure, through its error
    C        estimate and flag. Yet it must be stressed that the programs
    C        cannot be totally reliable.
    C        ------
    C
    C    (2) The user may want to examine the integrand function.
    C        If bad local difficulties occur, such as a discontinuity, a
    C        singularity, derivative singularity or high peak at one or
    C        more points within the interval, the first advice is to
    C        split up the interval at these points. The integrand must
    C        then be examined over each of the subintervals separately,
    C        so that a suitable integrator can be selected for each of
    C        them. If this yields problems involving relative accuracies
    C        to be imposed on -finite- subintervals, one can make use of
    C        QAGP, which must be provided with the positions of the local
    C        difficulties. However, if strong singularities are present
    C        and a high accuracy is requested, application of QAGS on the
    C        subintervals may yield a better result.
    C
    C        For quadrature over finite intervals we thus dispose of QAGS
    C        and
    C        - QNG for well-behaved integrands,
    C        - QAG for functions with an oscillating behaviour of a non
    C          specific type,
    C        - QAWO for functions, eventually singular, containing a
    C          factor COS(OMEGA*X) or SIN(OMEGA*X) where OMEGA is known,
    C        - QAWS for integrands with Algebraico-Logarithmic end point
    C          singularities of known type,
    C        - QAWC for Cauchy Principal Values.
    C
    C        Remark
    C        ------
    C        On return, the work arrays in the argument lists of the
    C        adaptive integrators contain information about the interval
    C        subdivision process and hence about the integrand behaviour:
    C        the end points of the subintervals, the local integral
    C        contributions and error estimates, and eventually other
    C        characteristics. For this reason, and because of its simple
    C        globally adaptive nature, the routine QAG in particular is
    C        well-suited for integrand examination. Difficult spots can
    C        be located by investigating the error estimates on the
    C        subintervals.
    C
    C        For infinite intervals we provide only one general-purpose
    C        routine, QAGI. It is based on the QAGS algorithm applied
    C        after a transformation of the original interval into (0,1).
    C        Yet it may eventuate that another type of transformation is
    C        more appropriate, or one might prefer to break up the
    C        original interval and use QAGI only on the infinite part
    C        and so on. These kinds of actions suggest a combined use of
    C        different QUADPACK integrators. Note that, when the only
    C        difficulty is an integrand singularity at the finite
    C        integration limit, it will in general not be necessary to
    C        break up the interval, as QAGI deals with several types of
    C        singularity at the boundary point of the integration range.
    C        It also handles slowly convergent improper integrals, on
    C        the condition that the integrand does not oscillate over
    C        the entire infinite interval. If it does we would advise
    C        to sum succeeding positive and negative contributions to
    C        the integral -e.g. integrate between the zeros- with one
    C        or more of the finite-range integrators, and apply
    C        convergence acceleration eventually by means of QUADPACK
    C        subroutine QELG which implements the Epsilon algorithm.
    C        Such quadrature problems include the Fourier transform as
    C        a special case. Yet for the latter we have an automatic
    C        integrator available, QAWF.
    C
    C *Long Description:
    C
    C 4. Example Programs
    C    ----------------
    C 4.1. Calling Program for QNG
    C      -----------------------
    C
    C            REAL A,ABSERR,B,F,EPSABS,EPSREL,RESULT
    C            INTEGER IER,NEVAL
    C            EXTERNAL F
    C            A = 0.0E0
    C            B = 1.0E0
    C            EPSABS = 0.0E0
    C            EPSREL = 1.0E-3
    C            CALL QNG(F,A,B,EPSABS,EPSREL,RESULT,ABSERR,NEVAL,IER)
    C      C  INCLUDE WRITE STATEMENTS
    C            STOP
    C            END
    C      C
    C            REAL FUNCTION F(X)
    C            REAL X
    C            F = EXP(X)/(X*X+0.1E+01)
    C            RETURN
    C            END
    C
    C 4.2. Calling Program for QAG
    C      -----------------------
    C
    C            REAL A,ABSERR,B,EPSABS,EPSREL,F,RESULT,WORK
    C            INTEGER IER,IWORK,KEY,LAST,LENW,LIMIT,NEVAL
    C            DIMENSION IWORK(100),WORK(400)
    C            EXTERNAL F
    C            A = 0.0E0
    C            B = 1.0E0
    C            EPSABS = 0.0E0
    C            EPSREL = 1.0E-3
    C            KEY = 6
    C            LIMIT = 100
    C            LENW = LIMIT*4
    C            CALL QAG(F,A,B,EPSABS,EPSREL,KEY,RESULT,ABSERR,NEVAL,
    C           *  IER,LIMIT,LENW,LAST,IWORK,WORK)
    C      C  INCLUDE WRITE STATEMENTS
    C            STOP
    C            END
    C      C
    C            REAL FUNCTION F(X)
    C            REAL X
    C            F = 2.0E0/(2.0E0+SIN(31.41592653589793E0*X))
    C            RETURN
    C            END
    C
    C 4.3. Calling Program for QAGS
    C      ------------------------
    C
    C            REAL A,ABSERR,B,EPSABS,EPSREL,F,RESULT,WORK
    C            INTEGER IER,IWORK,LAST,LENW,LIMIT,NEVAL
    C            DIMENSION IWORK(100),WORK(400)
    C            EXTERNAL F
    C            A = 0.0E0
    C            B = 1.0E0
    C            EPSABS = 0.0E0
    C            EPSREL = 1.0E-3
    C            LIMIT = 100
    C            LENW = LIMIT*4
    C            CALL QAGS(F,A,B,EPSABS,EPSREL,RESULT,ABSERR,NEVAL,IER,
    C           *  LIMIT,LENW,LAST,IWORK,WORK)
    C      C  INCLUDE WRITE STATEMENTS
    C            STOP
    C            END
    C      C
    C            REAL FUNCTION F(X)
    C            REAL X
    C            F = 0.0E0
    C            IF(X.GT.0.0E0) F = 1.0E0/SQRT(X)
    C            RETURN
    C            END
    C
    C 4.4. Calling Program for QAGP
    C      ------------------------
    C
    C            REAL A,ABSERR,B,EPSABS,EPSREL,F,POINTS,RESULT,WORK
    C            INTEGER IER,IWORK,LAST,LENIW,LENW,LIMIT,NEVAL,NPTS2
    C            DIMENSION IWORK(204),POINTS(4),WORK(404)
    C            EXTERNAL F
    C            A = 0.0E0
    C            B = 1.0E0
    C            NPTS2 = 4
    C            POINTS(1) = 1.0E0/7.0E0
    C            POINTS(2) = 2.0E0/3.0E0
    C            LIMIT = 100
    C            LENIW = LIMIT*2+NPTS2
    C            LENW = LIMIT*4+NPTS2
    C            CALL QAGP(F,A,B,NPTS2,POINTS,EPSABS,EPSREL,RESULT,ABSERR,
    C           *  NEVAL,IER,LENIW,LENW,LAST,IWORK,WORK)
    C      C  INCLUDE WRITE STATEMENTS
    C            STOP
    C            END
    C      C
    C            REAL FUNCTION F(X)
    C            REAL X
    C            F = 0.0E+00
    C            IF(X.NE.1.0E0/7.0E0.AND.X.NE.2.0E0/3.0E0) F =
    C           *  ABS(X-1.0E0/7.0E0)**(-0.25E0)*
    C           *  ABS(X-2.0E0/3.0E0)**(-0.55E0)
    C            RETURN
    C            END
    C
    C 4.5. Calling Program for QAGI
    C      ------------------------
    C
    C            REAL ABSERR,BOUN,EPSABS,EPSREL,F,RESULT,WORK
    C            INTEGER IER,INF,IWORK,LAST,LENW,LIMIT,NEVAL
    C            DIMENSION IWORK(100),WORK(400)
    C            EXTERNAL F
    C            BOUN = 0.0E0
    C            INF = 1
    C            EPSABS = 0.0E0
    C            EPSREL = 1.0E-3
    C            LIMIT = 100
    C            LENW = LIMIT*4
    C            CALL QAGI(F,BOUN,INF,EPSABS,EPSREL,RESULT,ABSERR,NEVAL,
    C           *  IER,LIMIT,LENW,LAST,IWORK,WORK)
    C      C  INCLUDE WRITE STATEMENTS
    C            STOP
    C            END
    C      C
    C            REAL FUNCTION F(X)
    C            REAL X
    C            F = 0.0E0
    C            IF(X.GT.0.0E0) F = SQRT(X)*LOG(X)/
    C           *                   ((X+1.0E0)*(X+2.0E0))
    C            RETURN
    C            END
    C
    C 4.6. Calling Program for QAWO
    C      ------------------------
    C
    C            REAL A,ABSERR,B,EPSABS,EPSREL,F,RESULT,OMEGA,WORK
    C            INTEGER IER,INTEGR,IWORK,LAST,LENIW,LENW,LIMIT,MAXP1,NEVAL
    C            DIMENSION IWORK(200),WORK(925)
    C            EXTERNAL F
    C            A = 0.0E0
    C            B = 1.0E0
    C            OMEGA = 10.0E0
    C            INTEGR = 1
    C            EPSABS = 0.0E0
    C            EPSREL = 1.0E-3
    C            LIMIT = 100
    C            LENIW = LIMIT*2
    C            MAXP1 = 21
    C            LENW = LIMIT*4+MAXP1*25
    C            CALL QAWO(F,A,B,OMEGA,INTEGR,EPSABS,EPSREL,RESULT,ABSERR,
    C           *  NEVAL,IER,LENIW,MAXP1,LENW,LAST,IWORK,WORK)
    C      C  INCLUDE WRITE STATEMENTS
    C            STOP
    C            END
    C      C
    C            REAL FUNCTION F(X)
    C            REAL X
    C            F = 0.0E0
    C            IF(X.GT.0.0E0) F = EXP(-X)*LOG(X)
    C            RETURN
    C            END
    C
    C 4.7. Calling Program for QAWF
    C      ------------------------
    C
    C            REAL A,ABSERR,EPSABS,F,RESULT,OMEGA,WORK
    C            INTEGER IER,INTEGR,IWORK,LAST,LENIW,LENW,LIMIT,LIMLST,
    C           *  LST,MAXP1,NEVAL
    C            DIMENSION IWORK(250),WORK(1025)
    C            EXTERNAL F
    C            A = 0.0E0
    C            OMEGA = 8.0E0
    C            INTEGR = 2
    C            EPSABS = 1.0E-3
    C            LIMLST = 50
    C            LIMIT = 100
    C            LENIW = LIMIT*2+LIMLST
    C            MAXP1 = 21
    C            LENW = LENIW*2+MAXP1*25
    C            CALL QAWF(F,A,OMEGA,INTEGR,EPSABS,RESULT,ABSERR,NEVAL,
    C           *  IER,LIMLST,LST,LENIW,MAXP1,LENW,IWORK,WORK)
    C      C  INCLUDE WRITE STATEMENTS
    C            STOP
    C            END
    C      C
    C            REAL FUNCTION F(X)
    C            REAL X
    C            IF(X.GT.0.0E0) F = SIN(50.0E0*X)/(X*SQRT(X))
    C            RETURN
    C            END
    C
    C 4.8. Calling Program for QAWS
    C      ------------------------
    C
    C            REAL A,ABSERR,ALFA,B,BETA,EPSABS,EPSREL,F,RESULT,WORK
    C            INTEGER IER,INTEGR,IWORK,LAST,LENW,LIMIT,NEVAL
    C            DIMENSION IWORK(100),WORK(400)
    C            EXTERNAL F
    C            A = 0.0E0
    C            B = 1.0E0
    C            ALFA = -0.5E0
    C            BETA = -0.5E0
    C            INTEGR = 1
    C            EPSABS = 0.0E0
    C            EPSREL = 1.0E-3
    C            LIMIT = 100
    C            LENW = LIMIT*4
    C            CALL QAWS(F,A,B,ALFA,BETA,INTEGR,EPSABS,EPSREL,RESULT,
    C           *  ABSERR,NEVAL,IER,LIMIT,LENW,LAST,IWORK,WORK)
    C      C  INCLUDE WRITE STATEMENTS
    C            STOP
    C            END
    C      C
    C            REAL FUNCTION F(X)
    C            REAL X
    C            F = SIN(10.0E0*X)
    C            RETURN
    C            END
    C
    C 4.9. Calling Program for QAWC
    C      ------------------------
    C
    C            REAL A,ABSERR,B,C,EPSABS,EPSREL,F,RESULT,WORK
    C            INTEGER IER,IWORK,LAST,LENW,LIMIT,NEVAL
    C            DIMENSION IWORK(100),WORK(400)
    C            EXTERNAL F
    C            A = -1.0E0
    C            B = 1.0E0
    C            C = 0.5E0
    C            EPSABS = 0.0E0
    C            EPSREL = 1.0E-3
    C            LIMIT = 100
    C            LENW = LIMIT*4
    C            CALL QAWC(F,A,B,C,EPSABS,EPSREL,RESULT,ABSERR,NEVAL,
    C           *  IER,LIMIT,LENW,LAST,IWORK,WORK)
    C      C  INCLUDE WRITE STATEMENTS
    C            STOP
    C            END
    C      C
    C            REAL FUNCTION F(X)
    C            REAL X
    C            F = 1.0E0/(X*X+1.0E-4)
    C            RETURN
    C            END
    C
    C***REFERENCES  (NONE)
    C***ROUTINES CALLED  (NONE)
    C***REVISION HISTORY  (YYMMDD)
    C   810401  DATE WRITTEN
    C   890531  Changed all specific intrinsics to generic.  (WRB)
    C   890531  REVISION DATE from Version 3.2
    C   891214  Prologue converted to Version 4.0 format.  (BAB)
    C   900723  PURPOSE section revised.  (WRB)
    C***END PROLOGUE  QPDOC