src/numerical/slatec/fortran/dqawo.f

   1 *DECK DQAWO
   2       SUBROUTINE DQAWO (F, A, B, OMEGA, INTEGR, EPSABS, EPSREL, RESULT,
   3      +   ABSERR, NEVAL, IER, LENIW, MAXP1, LENW, LAST, IWORK, WORK)
   4 C***BEGIN PROLOGUE  DQAWO
   5 C***PURPOSE  Calculate an approximation to a given definite integral
   6 C            I= Integral of F(X)*W(X) over (A,B), where
   7 C                   W(X) = COS(OMEGA*X)
   8 C               or  W(X) = SIN(OMEGA*X),
   9 C            hopefully satisfying the following claim for accuracy
  10 C                ABS(I-RESULT).LE.MAX(EPSABS,EPSREL*ABS(I)).
  11 C***LIBRARY   SLATEC (QUADPACK)
  12 C***CATEGORY  H2A2A1
  13 C***TYPE      DOUBLE PRECISION (QAWO-S, DQAWO-D)
  14 C***KEYWORDS  AUTOMATIC INTEGRATOR, CLENSHAW-CURTIS METHOD,
  15 C             EXTRAPOLATION, GLOBALLY ADAPTIVE,
  16 C             INTEGRAND WITH OSCILLATORY COS OR SIN FACTOR, QUADPACK,
  17 C             QUADRATURE, SPECIAL-PURPOSE
  18 C***AUTHOR  Piessens, Robert
  19 C             Applied Mathematics and Programming Division
  20 C             K. U. Leuven
  21 C           de Doncker, Elise
  22 C             Applied Mathematics and Programming Division
  23 C             K. U. Leuven
  24 C***DESCRIPTION
  25 C
  26 C        Computation of oscillatory integrals
  27 C        Standard fortran subroutine
  28 C        Double precision version
  29 C
  30 C        PARAMETERS
  31 C         ON ENTRY
  32 C            F      - Double precision
  33 C                     Function subprogram defining the function
  34 C                     F(X).  The actual name for F needs to be
  35 C                     declared E X T E R N A L in the driver program.
  36 C
  37 C            A      - Double precision
  38 C                     Lower limit of integration
  39 C
  40 C            B      - Double precision
  41 C                     Upper limit of integration
  42 C
  43 C            OMEGA  - Double precision
  44 C                     Parameter in the integrand weight function
  45 C
  46 C            INTEGR - Integer
  47 C                     Indicates which of the weight functions is used
  48 C                     INTEGR = 1      W(X) = COS(OMEGA*X)
  49 C                     INTEGR = 2      W(X) = SIN(OMEGA*X)
  50 C                     If INTEGR.NE.1.AND.INTEGR.NE.2, the routine will
  51 C                     end with IER = 6.
  52 C
  53 C            EPSABS - Double precision
  54 C                     Absolute accuracy requested
  55 C            EPSREL - Double precision
  56 C                     Relative accuracy requested
  57 C                     If EPSABS.LE.0 and
  58 C                     EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28),
  59 C                     the routine will end with IER = 6.
  60 C
  61 C         ON RETURN
  62 C            RESULT - Double precision
  63 C                     Approximation to the integral
  64 C
  65 C            ABSERR - Double precision
  66 C                     Estimate of the modulus of the absolute error,
  67 C                     which should equal or exceed ABS(I-RESULT)
  68 C
  69 C            NEVAL  - Integer
  70 C                     Number of integrand evaluations
  71 C
  72 C            IER    - Integer
  73 C                     IER = 0 Normal and reliable termination of the
  74 C                             routine. It is assumed that the requested
  75 C                             accuracy has been achieved.
  76 C                   - IER.GT.0 Abnormal termination of the routine.
  77 C                             The estimates for integral and error are
  78 C                             less reliable. It is assumed that the
  79 C                             requested accuracy has not been achieved.
  80 C            ERROR MESSAGES
  81 C                     IER = 1 Maximum number of subdivisions allowed
  82 C                             has been achieved (= LENIW/2). One can
  83 C                             allow more subdivisions by increasing the
  84 C                             value of LENIW (and taking the according
  85 C                             dimension adjustments into account).
  86 C                             However, if this yields no improvement it
  87 C                             is advised to analyze the integrand in
  88 C                             order to determine the integration
  89 C                             difficulties. If the position of a local
  90 C                             difficulty can be determined (e.g.
  91 C                             SINGULARITY, DISCONTINUITY within the
  92 C                             interval) one will probably gain from
  93 C                             splitting up the interval at this point
  94 C                             and calling the integrator on the
  95 C                             subranges. If possible, an appropriate
  96 C                             special-purpose integrator should be used
  97 C                             which is designed for handling the type of
  98 C                             difficulty involved.
  99 C                         = 2 The occurrence of roundoff error is
 100 C                             detected, which prevents the requested
 101 C                             tolerance from being achieved.
 102 C                             The error may be under-estimated.
 103 C                         = 3 Extremely bad integrand behaviour occurs
 104 C                             at some interior points of the
 105 C                             integration interval.
 106 C                         = 4 The algorithm does not converge.
 107 C                             Roundoff error is detected in the
 108 C                             extrapolation table. It is presumed that
 109 C                             the requested tolerance cannot be achieved
 110 C                             due to roundoff in the extrapolation
 111 C                             table, and that the returned result is
 112 C                             the best which can be obtained.
 113 C                         = 5 The integral is probably divergent, or
 114 C                             slowly convergent. It must be noted that
 115 C                             divergence can occur with any other value
 116 C                             of IER.
 117 C                         = 6 The input is invalid, because
 118 C                             (EPSABS.LE.0 and
 119 C                              EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28))
 120 C                             or (INTEGR.NE.1 AND INTEGR.NE.2),
 121 C                             or LENIW.LT.2 OR MAXP1.LT.1 or
 122 C                             LENW.LT.LENIW*2+MAXP1*25.
 123 C                             RESULT, ABSERR, NEVAL, LAST are set to
 124 C                             zero. Except when LENIW, MAXP1 or LENW are
 125 C                             invalid, WORK(LIMIT*2+1), WORK(LIMIT*3+1),
 126 C                             IWORK(1), IWORK(LIMIT+1) are set to zero,
 127 C                             WORK(1) is set to A and WORK(LIMIT+1) to
 128 C                             B.
 129 C
 130 C         DIMENSIONING PARAMETERS
 131 C            LENIW  - Integer
 132 C                     Dimensioning parameter for IWORK.
 133 C                     LENIW/2 equals the maximum number of subintervals
 134 C                     allowed in the partition of the given integration
 135 C                     interval (A,B), LENIW.GE.2.
 136 C                     If LENIW.LT.2, the routine will end with IER = 6.
 137 C
 138 C            MAXP1  - Integer
 139 C                     Gives an upper bound on the number of Chebyshev
 140 C                     moments which can be stored, i.e. for the
 141 C                     intervals of lengths ABS(B-A)*2**(-L),
 142 C                     L=0,1, ..., MAXP1-2, MAXP1.GE.1
 143 C                     If MAXP1.LT.1, the routine will end with IER = 6.
 144 C
 145 C            LENW   - Integer
 146 C                     Dimensioning parameter for WORK
 147 C                     LENW must be at least LENIW*2+MAXP1*25.
 148 C                     If LENW.LT.(LENIW*2+MAXP1*25), the routine will
 149 C                     end with IER = 6.
 150 C
 151 C            LAST   - Integer
 152 C                     On return, LAST equals the number of subintervals
 153 C                     produced in the subdivision process, which
 154 C                     determines the number of significant elements
 155 C                     actually in the WORK ARRAYS.
 156 C
 157 C         WORK ARRAYS
 158 C            IWORK  - Integer
 159 C                     Vector of dimension at least LENIW
 160 C                     on return, the first K elements of which contain
 161 C                     pointers to the error estimates over the
 162 C                     subintervals, such that WORK(LIMIT*3+IWORK(1)), ..
 163 C                     WORK(LIMIT*3+IWORK(K)) form a decreasing
 164 C                     sequence, with LIMIT = LENW/2 , and K = LAST
 165 C                     if LAST.LE.(LIMIT/2+2), and K = LIMIT+1-LAST
 166 C                     otherwise.
 167 C                     Furthermore, IWORK(LIMIT+1), ..., IWORK(LIMIT+
 168 C                     LAST) indicate the subdivision levels of the
 169 C                     subintervals, such that IWORK(LIMIT+I) = L means
 170 C                     that the subinterval numbered I is of length
 171 C                     ABS(B-A)*2**(1-L).
 172 C
 173 C            WORK   - Double precision
 174 C                     Vector of dimension at least LENW
 175 C                     On return
 176 C                     WORK(1), ..., WORK(LAST) contain the left
 177 C                      end points of the subintervals in the
 178 C                      partition of (A,B),
 179 C                     WORK(LIMIT+1), ..., WORK(LIMIT+LAST) contain
 180 C                      the right end points,
 181 C                     WORK(LIMIT*2+1), ..., WORK(LIMIT*2+LAST) contain
 182 C                      the integral approximations over the
 183 C                      subintervals,
 184 C                     WORK(LIMIT*3+1), ..., WORK(LIMIT*3+LAST)
 185 C                      contain the error estimates.
 186 C                     WORK(LIMIT*4+1), ..., WORK(LIMIT*4+MAXP1*25)
 187 C                      Provide space for storing the Chebyshev moments.
 188 C                     Note that LIMIT = LENW/2.
 189 C
 190 C***REFERENCES  (NONE)
 191 C***ROUTINES CALLED  DQAWOE, XERMSG
 192 C***REVISION HISTORY  (YYMMDD)
 193 C   800101  DATE WRITTEN
 194 C   890831  Modified array declarations.  (WRB)
 195 C   890831  REVISION DATE from Version 3.2
 196 C   891214  Prologue converted to Version 4.0 format.  (BAB)
 197 C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
 198 C***END PROLOGUE  DQAWO
 199 C
 200        DOUBLE PRECISION A,ABSERR,B,EPSABS,EPSREL,F,OMEGA,RESULT,WORK
 201        INTEGER IER,INTEGR,IWORK,LAST,LIMIT,LENW,LENIW,LVL,L1,L2,L3,L4,
 202      1  MAXP1,MOMCOM,NEVAL
 203 C
 204        DIMENSION IWORK(*),WORK(*)
 205 C
 206        EXTERNAL F
 207 C
 208 C         CHECK VALIDITY OF LENIW, MAXP1 AND LENW.
 209 C
 210 C***FIRST EXECUTABLE STATEMENT  DQAWO
 211       IER = 6
 212       NEVAL = 0
 213       LAST = 0
 214       RESULT = 0.0D+00
 215       ABSERR = 0.0D+00
 216       IF(LENIW.LT.2.OR.MAXP1.LT.1.OR.LENW.LT.(LENIW*2+MAXP1*25))
 217      1  GO TO 10
 218 C
 219 C         PREPARE CALL FOR DQAWOE
 220 C
 221       LIMIT = LENIW/2
 222       L1 = LIMIT+1
 223       L2 = LIMIT+L1
 224       L3 = LIMIT+L2
 225       L4 = LIMIT+L3
 226       CALL DQAWOE(F,A,B,OMEGA,INTEGR,EPSABS,EPSREL,LIMIT,1,MAXP1,RESULT,
 227      1   ABSERR,NEVAL,IER,LAST,WORK(1),WORK(L1),WORK(L2),WORK(L3),
 228      2   IWORK(1),IWORK(L1),MOMCOM,WORK(L4))
 229 C
 230 C         CALL ERROR HANDLER IF NECESSARY
 231 C
 232       LVL = 0
 233 10    IF(IER.EQ.6) LVL = 0
 234       IF (IER .NE. 0) CALL XERMSG ('SLATEC', 'DQAWO',
 235      +   'ABNORMAL RETURN', IER, LVL)
 236       RETURN
 237       END