relibc_openlibm/slatec/dqawf.f

*DECK DQAWF
      SUBROUTINE DQAWF (F, A, OMEGA, INTEGR, EPSABS, RESULT, ABSERR,
     +   NEVAL, IER, LIMLST, LST, LENIW, MAXP1, LENW, IWORK, WORK)
C***BEGIN PROLOGUE  DQAWF
C***PURPOSE  The routine calculates an approximation result to a given
C            Fourier integral I=Integral of F(X)*W(X) over (A,INFINITY)
C            where W(X) = COS(OMEGA*X) or W(X) = SIN(OMEGA*X).
C            Hopefully satisfying following claim for accuracy
C            ABS(I-RESULT).LE.EPSABS.
C***LIBRARY   SLATEC (QUADPACK)
C***CATEGORY  H2A3A1
C***TYPE      DOUBLE PRECISION (QAWF-S, DQAWF-D)
C***KEYWORDS  AUTOMATIC INTEGRATOR, CONVERGENCE ACCELERATION,
C             FOURIER INTEGRALS, INTEGRATION BETWEEN ZEROS, QUADPACK,
C             QUADRATURE, SPECIAL-PURPOSE INTEGRAL
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***DESCRIPTION
C
C        Computation of Fourier integrals
C        Standard fortran subroutine
C        Double precision version
C
C
C        PARAMETERS
C         ON ENTRY
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            A      - Double precision
C                     Lower limit of integration
C
C            OMEGA  - Double precision
C                     Parameter in the integrand WEIGHT function
C
C            INTEGR - Integer
C                     Indicates which of the WEIGHT functions is used
C                     INTEGR = 1      W(X) = COS(OMEGA*X)
C                     INTEGR = 2      W(X) = SIN(OMEGA*X)
C                     IF INTEGR.NE.1.AND.INTEGR.NE.2, the routine
C                     will end with IER = 6.
C
C            EPSABS - Double precision
C                     Absolute accuracy requested, EPSABS.GT.0.
C                     If EPSABS.LE.0, the routine will end with IER = 6.
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.
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                    If OMEGA.NE.0
C                     IER = 1 Maximum number of cycles allowed
C                             has been achieved, i.e. of subintervals
C                             (A+(K-1)C,A+KC) where
C                             C = (2*INT(ABS(OMEGA))+1)*PI/ABS(OMEGA),
C                             FOR K = 1, 2, ..., LST.
C                             One can allow more cycles by increasing
C                             the value of LIMLST (and taking the
C                             according dimension adjustments into
C                             account). Examine the array IWORK which
C                             contains the error flags on the cycles, in
C                             order to look for eventual local
C                             integration difficulties.
C                             If the position of a local difficulty
C                             can 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
C                             appropriate integrators on the subranges.
C                         = 4 The extrapolation table constructed for
C                             convergence acceleration of the series
C                             formed by the integral contributions over
C                             the cycles, does not converge to within
C                             the requested accuracy.
C                             As in the case of IER = 1, it is advised
C                             to examine the array IWORK which contains
C                             the error flags on the cycles.
C                         = 6 The input is invalid because
C                             (INTEGR.NE.1 AND INTEGR.NE.2) or
C                              EPSABS.LE.0 or LIMLST.LT.1 or
C                              LENIW.LT.(LIMLST+2) or MAXP1.LT.1 or
C                              LENW.LT.(LENIW*2+MAXP1*25).
C                              RESULT, ABSERR, NEVAL, LST are set to
C                              zero.
C                         = 7 Bad integrand behaviour occurs within
C                             one or more of the cycles. Location and
C                             type of the difficulty involved can be
C                             determined from the first LST elements of
C                             vector IWORK.  Here LST is the number of
C                             cycles actually needed (see below).
C                             IWORK(K) = 1 The maximum number of
C                                          subdivisions (=(LENIW-LIMLST)
C                                          /2) has been achieved on the
C                                          K th cycle.
C                                      = 2 Occurrence of roundoff error
C                                          is detected and prevents the
C                                          tolerance imposed on the K th
C                                          cycle, from being achieved
C                                          on this cycle.
C                                      = 3 Extremely bad integrand
C                                          behaviour occurs at some
C                                          points of the K th cycle.
C                                      = 4 The integration procedure
C                                          over the K th cycle does
C                                          not converge (to within the
C                                          required accuracy) due to
C                                          roundoff in the extrapolation
C                                          procedure invoked on this
C                                          cycle. It is assumed that the
C                                          result on this interval is
C                                          the best which can be
C                                          obtained.
C                                      = 5 The integral over the K th
C                                          cycle is probably divergent
C                                          or slowly convergent. It must
C                                          be noted that divergence can
C                                          occur with any other value of
C                                          IWORK(K).
C                    If OMEGA = 0 and INTEGR = 1,
C                    The integral is calculated by means of DQAGIE,
C                    and IER = IWORK(1) (with meaning as described
C                    for IWORK(K),K = 1).
C
C         DIMENSIONING PARAMETERS
C            LIMLST - Integer
C                     LIMLST gives an upper bound on the number of
C                     cycles, LIMLST.GE.3.
C                     If LIMLST.LT.3, the routine will end with IER = 6.
C
C            LST    - Integer
C                     On return, LST indicates the number of cycles
C                     actually needed for the integration.
C                     If OMEGA = 0, then LST is set to 1.
C
C            LENIW  - Integer
C                     Dimensioning parameter for IWORK. On entry,
C                     (LENIW-LIMLST)/2 equals the maximum number of
C                     subintervals allowed in the partition of each
C                     cycle, LENIW.GE.(LIMLST+2).
C                     If LENIW.LT.(LIMLST+2), the routine will end with
C                     IER = 6.
C
C            MAXP1  - Integer
C                     MAXP1 gives an upper bound on the number of
C                     Chebyshev moments which can be stored, i.e. for
C                     the intervals of lengths ABS(B-A)*2**(-L),
C                     L = 0,1, ..., MAXP1-2, MAXP1.GE.1.
C                     If MAXP1.LT.1, the routine will end with IER = 6.
C            LENW   - Integer
C                     Dimensioning parameter for WORK
C                     LENW must be at least LENIW*2+MAXP1*25.
C                     If LENW.LT.(LENIW*2+MAXP1*25), the routine will
C                     end with IER = 6.
C
C         WORK ARRAYS
C            IWORK  - Integer
C                     Vector of dimension at least LENIW
C                     On return, IWORK(K) FOR K = 1, 2, ..., LST
C                     contain the error flags on the cycles.
C
C            WORK   - Double precision
C                     Vector of dimension at least
C                     On return,
C                     WORK(1), ..., WORK(LST) contain the integral
C                      approximations over the cycles,
C                     WORK(LIMLST+1), ..., WORK(LIMLST+LST) contain
C                      the error estimates over the cycles.
C                     further elements of WORK have no specific
C                     meaning for the user.
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  DQAWFE, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   800101  DATE WRITTEN
C   890831  Modified array declarations.  (WRB)
C   891009  Removed unreferenced variable.  (WRB)
C   891009  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C***END PROLOGUE  DQAWF
C
       DOUBLE PRECISION A,ABSERR,EPSABS,F,OMEGA,RESULT,WORK
       INTEGER IER,INTEGR,IWORK,LENIW,LENW,LIMIT,LIMLST,LL2,LVL,
     1  LST,L1,L2,L3,L4,L5,L6,MAXP1,NEVAL
C
       DIMENSION IWORK(*),WORK(*)
C
       EXTERNAL F
C
C         CHECK VALIDITY OF LIMLST, LENIW, MAXP1 AND LENW.
C
C***FIRST EXECUTABLE STATEMENT  DQAWF
      IER = 6
      NEVAL = 0
      RESULT = 0.0D+00
      ABSERR = 0.0D+00
      IF(LIMLST.LT.3.OR.LENIW.LT.(LIMLST+2).OR.MAXP1.LT.1.OR.LENW.LT.
     1   (LENIW*2+MAXP1*25)) GO TO 10
C
C         PREPARE CALL FOR DQAWFE
C
      LIMIT = (LENIW-LIMLST)/2
      L1 = LIMLST+1
      L2 = LIMLST+L1
      L3 = LIMIT+L2
      L4 = LIMIT+L3
      L5 = LIMIT+L4
      L6 = LIMIT+L5
      LL2 = LIMIT+L1
      CALL DQAWFE(F,A,OMEGA,INTEGR,EPSABS,LIMLST,LIMIT,MAXP1,RESULT,
     1  ABSERR,NEVAL,IER,WORK(1),WORK(L1),IWORK(1),LST,WORK(L2),
     2  WORK(L3),WORK(L4),WORK(L5),IWORK(L1),IWORK(LL2),WORK(L6))
C
C         CALL ERROR HANDLER IF NECESSARY
C
      LVL = 0
10    IF(IER.EQ.6) LVL = 1
      IF (IER .NE. 0) CALL XERMSG ('SLATEC', 'DQAWF',
     +   'ABNORMAL RETURN', IER, LVL)
      RETURN
      END