src/numerical/slatec/fortran/dqawfe.f

   1 *DECK DQAWFE
   2       SUBROUTINE DQAWFE (F, A, OMEGA, INTEGR, EPSABS, LIMLST, LIMIT,
   3      +   MAXP1, RESULT, ABSERR, NEVAL, IER, RSLST, ERLST, IERLST, LST,
   4      +   ALIST, BLIST, RLIST, ELIST, IORD, NNLOG, CHEBMO)
   5 C***BEGIN PROLOGUE  DQAWFE
   6 C***PURPOSE  The routine calculates an approximation result to a
   7 C            given Fourier integral
   8 C            I = Integral of F(X)*W(X) over (A,INFINITY)
   9 C            where W(X)=COS(OMEGA*X) or W(X)=SIN(OMEGA*X),
  10 C            hopefully satisfying following claim for accuracy
  11 C            ABS(I-RESULT).LE.EPSABS.
  12 C***LIBRARY   SLATEC (QUADPACK)
  13 C***CATEGORY  H2A3A1
  14 C***TYPE      DOUBLE PRECISION (QAWFE-S, DQAWFE-D)
  15 C***KEYWORDS  AUTOMATIC INTEGRATOR, CONVERGENCE ACCELERATION,
  16 C             FOURIER INTEGRALS, INTEGRATION BETWEEN ZEROS, QUADPACK,
  17 C             QUADRATURE, SPECIAL-PURPOSE INTEGRAL
  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 Fourier 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 integrand
  34 C                     Function F(X). The actual name for F needs to
  35 C                     be 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            OMEGA  - Double precision
  41 C                     Parameter in the WEIGHT function
  42 C
  43 C            INTEGR - Integer
  44 C                     Indicates which WEIGHT function is used
  45 C                     INTEGR = 1      W(X) = COS(OMEGA*X)
  46 C                     INTEGR = 2      W(X) = SIN(OMEGA*X)
  47 C                     If INTEGR.NE.1.AND.INTEGR.NE.2, the routine will
  48 C                     end with IER = 6.
  49 C
  50 C            EPSABS - Double precision
  51 C                     absolute accuracy requested, EPSABS.GT.0
  52 C                     If EPSABS.LE.0, the routine will end with IER = 6.
  53 C
  54 C            LIMLST - Integer
  55 C                     LIMLST gives an upper bound on the number of
  56 C                     cycles, LIMLST.GE.1.
  57 C                     If LIMLST.LT.3, the routine will end with IER = 6.
  58 C
  59 C            LIMIT  - Integer
  60 C                     Gives an upper bound on the number of subintervals
  61 C                     allowed in the partition of each cycle, LIMIT.GE.1
  62 C                     each cycle, LIMIT.GE.1.
  63 C
  64 C            MAXP1  - Integer
  65 C                     Gives an upper bound on the number of
  66 C                     Chebyshev moments which can be stored, I.E.
  67 C                     for the intervals of lengths ABS(B-A)*2**(-L),
  68 C                     L=0,1, ..., MAXP1-2, MAXP1.GE.1
  69 C
  70 C         ON RETURN
  71 C            RESULT - Double precision
  72 C                     Approximation to the integral X
  73 C
  74 C            ABSERR - Double precision
  75 C                     Estimate of the modulus of the absolute error,
  76 C                     which should equal or exceed ABS(I-RESULT)
  77 C
  78 C            NEVAL  - Integer
  79 C                     Number of integrand evaluations
  80 C
  81 C            IER    - IER = 0 Normal and reliable termination of
  82 C                             the routine. It is assumed that the
  83 C                             requested accuracy has been achieved.
  84 C                     IER.GT.0 Abnormal termination of the routine. The
  85 C                             estimates for integral and error are less
  86 C                             reliable. It is assumed that the requested
  87 C                             accuracy has not been achieved.
  88 C            ERROR MESSAGES
  89 C                    If OMEGA.NE.0
  90 C                     IER = 1 Maximum number of  cycles  allowed
  91 C                             Has been achieved., i.e. of subintervals
  92 C                             (A+(K-1)C,A+KC) where
  93 C                             C = (2*INT(ABS(OMEGA))+1)*PI/ABS(OMEGA),
  94 C                             for K = 1, 2, ..., LST.
  95 C                             One can allow more cycles by increasing
  96 C                             the value of LIMLST (and taking the
  97 C                             according dimension adjustments into
  98 C                             account).
  99 C                             Examine the array IWORK which contains
 100 C                             the error flags on the cycles, in order to
 101 C                             look for eventual local integration
 102 C                             difficulties. If the position of a local
 103 C                             difficulty can be determined (e.g.
 104 C                             SINGULARITY, DISCONTINUITY within the
 105 C                             interval) one will probably gain from
 106 C                             splitting up the interval at this point
 107 C                             and calling appropriate integrators on
 108 C                             the subranges.
 109 C                         = 4 The extrapolation table constructed for
 110 C                             convergence acceleration of the series
 111 C                             formed by the integral contributions over
 112 C                             the cycles, does not converge to within
 113 C                             the requested accuracy. As in the case of
 114 C                             IER = 1, it is advised to examine the
 115 C                             array IWORK which contains the error
 116 C                             flags on the cycles.
 117 C                         = 6 The input is invalid because
 118 C                             (INTEGR.NE.1 AND INTEGR.NE.2) or
 119 C                              EPSABS.LE.0 or LIMLST.LT.3.
 120 C                              RESULT, ABSERR, NEVAL, LST are set
 121 C                              to zero.
 122 C                         = 7 Bad integrand behaviour occurs within one
 123 C                             or more of the cycles. Location and type
 124 C                             of the difficulty involved can be
 125 C                             determined from the vector IERLST. Here
 126 C                             LST is the number of cycles actually
 127 C                             needed (see below).
 128 C                             IERLST(K) = 1 The maximum number of
 129 C                                           subdivisions (= LIMIT) has
 130 C                                           been achieved on the K th
 131 C                                           cycle.
 132 C                                       = 2 Occurrence of roundoff error
 133 C                                           is detected and prevents the
 134 C                                           tolerance imposed on the
 135 C                                           K th cycle, from being
 136 C                                           achieved.
 137 C                                       = 3 Extremely bad integrand
 138 C                                           behaviour occurs at some
 139 C                                           points of the K th cycle.
 140 C                                       = 4 The integration procedure
 141 C                                           over the K th cycle does
 142 C                                           not converge (to within the
 143 C                                           required accuracy) due to
 144 C                                           roundoff in the
 145 C                                           extrapolation procedure
 146 C                                           invoked on this cycle. It
 147 C                                           is assumed that the result
 148 C                                           on this interval is the
 149 C                                           best which can be obtained.
 150 C                                       = 5 The integral over the K th
 151 C                                           cycle is probably divergent
 152 C                                           or slowly convergent. It
 153 C                                           must be noted that
 154 C                                           divergence can occur with
 155 C                                           any other value of
 156 C                                           IERLST(K).
 157 C                    If OMEGA = 0 and INTEGR = 1,
 158 C                    The integral is calculated by means of DQAGIE
 159 C                    and IER = IERLST(1) (with meaning as described
 160 C                    for IERLST(K), K = 1).
 161 C
 162 C            RSLST  - Double precision
 163 C                     Vector of dimension at least LIMLST
 164 C                     RSLST(K) contains the integral contribution
 165 C                     over the interval (A+(K-1)C,A+KC) where
 166 C                     C = (2*INT(ABS(OMEGA))+1)*PI/ABS(OMEGA),
 167 C                     K = 1, 2, ..., LST.
 168 C                     Note that, if OMEGA = 0, RSLST(1) contains
 169 C                     the value of the integral over (A,INFINITY).
 170 C
 171 C            ERLST  - Double precision
 172 C                     Vector of dimension at least LIMLST
 173 C                     ERLST(K) contains the error estimate corresponding
 174 C                     with RSLST(K).
 175 C
 176 C            IERLST - Integer
 177 C                     Vector of dimension at least LIMLST
 178 C                     IERLST(K) contains the error flag corresponding
 179 C                     with RSLST(K). For the meaning of the local error
 180 C                     flags see description of output parameter IER.
 181 C
 182 C            LST    - Integer
 183 C                     Number of subintervals needed for the integration
 184 C                     If OMEGA = 0 then LST is set to 1.
 185 C
 186 C            ALIST, BLIST, RLIST, ELIST - Double precision
 187 C                     vector of dimension at least LIMIT,
 188 C
 189 C            IORD, NNLOG - Integer
 190 C                     Vector of dimension at least LIMIT, providing
 191 C                     space for the quantities needed in the subdivision
 192 C                     process of each cycle
 193 C
 194 C            CHEBMO - Double precision
 195 C                     Array of dimension at least (MAXP1,25), providing
 196 C                     space for the Chebyshev moments needed within the
 197 C                     cycles
 198 C
 199 C***REFERENCES  (NONE)
 200 C***ROUTINES CALLED  D1MACH, DQAGIE, DQAWOE, DQELG
 201 C***REVISION HISTORY  (YYMMDD)
 202 C   800101  DATE WRITTEN
 203 C   890531  Changed all specific intrinsics to generic.  (WRB)
 204 C   890831  Modified array declarations.  (WRB)
 205 C   891009  Removed unreferenced variable.  (WRB)
 206 C   891009  REVISION DATE from Version 3.2
 207 C   891214  Prologue converted to Version 4.0 format.  (BAB)
 208 C***END PROLOGUE  DQAWFE
 209 C
 210       DOUBLE PRECISION A,ABSEPS,ABSERR,ALIST,BLIST,CHEBMO,CORREC,CYCLE,
 211      1  C1,C2,DL,DRL,D1MACH,ELIST,ERLST,EP,EPS,EPSA,
 212      2  EPSABS,ERRSUM,F,FACT,OMEGA,P,PI,P1,PSUM,RESEPS,RESULT,RES3LA,
 213      3  RLIST,RSLST,UFLOW
 214       INTEGER IER,IERLST,INTEGR,IORD,KTMIN,L,LAST,LST,LIMIT,LIMLST,LL,
 215      1    MAXP1,MOMCOM,NEV,NEVAL,NNLOG,NRES,NUMRL2
 216 C
 217       DIMENSION ALIST(*),BLIST(*),CHEBMO(MAXP1,25),ELIST(*),
 218      1  ERLST(*),IERLST(*),IORD(*),NNLOG(*),PSUM(52),
 219      2  RES3LA(3),RLIST(*),RSLST(*)
 220 C
 221       EXTERNAL F
 222 C
 223 C
 224 C            THE DIMENSION OF  PSUM  IS DETERMINED BY THE VALUE OF
 225 C            LIMEXP IN SUBROUTINE DQELG (PSUM MUST BE OF DIMENSION
 226 C            (LIMEXP+2) AT LEAST).
 227 C
 228 C           LIST OF MAJOR VARIABLES
 229 C           -----------------------
 230 C
 231 C           C1, C2    - END POINTS OF SUBINTERVAL (OF LENGTH CYCLE)
 232 C           CYCLE     - (2*INT(ABS(OMEGA))+1)*PI/ABS(OMEGA)
 233 C           PSUM      - VECTOR OF DIMENSION AT LEAST (LIMEXP+2)
 234 C                       (SEE ROUTINE DQELG)
 235 C                       PSUM CONTAINS THE PART OF THE EPSILON TABLE
 236 C                       WHICH IS STILL NEEDED FOR FURTHER COMPUTATIONS.
 237 C                       EACH ELEMENT OF PSUM IS A PARTIAL SUM OF THE
 238 C                       SERIES WHICH SHOULD SUM TO THE VALUE OF THE
 239 C                       INTEGRAL.
 240 C           ERRSUM    - SUM OF ERROR ESTIMATES OVER THE SUBINTERVALS,
 241 C                       CALCULATED CUMULATIVELY
 242 C           EPSA      - ABSOLUTE TOLERANCE REQUESTED OVER CURRENT
 243 C                       SUBINTERVAL
 244 C           CHEBMO    - ARRAY CONTAINING THE MODIFIED CHEBYSHEV
 245 C                       MOMENTS (SEE ALSO ROUTINE DQC25F)
 246 C
 247       SAVE P, PI
 248       DATA P/0.9D+00/
 249       DATA PI / 3.1415926535 8979323846 2643383279 50 D0 /
 250 C
 251 C           TEST ON VALIDITY OF PARAMETERS
 252 C           ------------------------------
 253 C
 254 C***FIRST EXECUTABLE STATEMENT  DQAWFE
 255       RESULT = 0.0D+00
 256       ABSERR = 0.0D+00
 257       NEVAL = 0
 258       LST = 0
 259       IER = 0
 260       IF((INTEGR.NE.1.AND.INTEGR.NE.2).OR.EPSABS.LE.0.0D+00.OR.
 261      1  LIMLST.LT.3) IER = 6
 262       IF(IER.EQ.6) GO TO 999
 263       IF(OMEGA.NE.0.0D+00) GO TO 10
 264 C
 265 C           INTEGRATION BY DQAGIE IF OMEGA IS ZERO
 266 C           --------------------------------------
 267 C
 268       IF(INTEGR.EQ.1) CALL DQAGIE(F,A,1,EPSABS,0.0D+00,LIMIT,
 269      1  RESULT,ABSERR,NEVAL,IER,ALIST,BLIST,RLIST,ELIST,IORD,LAST)
 270       RSLST(1) = RESULT
 271       ERLST(1) = ABSERR
 272       IERLST(1) = IER
 273       LST = 1
 274       GO TO 999
 275 C
 276 C           INITIALIZATIONS
 277 C           ---------------
 278 C
 279    10 L = ABS(OMEGA)
 280       DL = 2*L+1
 281       CYCLE = DL*PI/ABS(OMEGA)
 282       IER = 0
 283       KTMIN = 0
 284       NEVAL = 0
 285       NUMRL2 = 0
 286       NRES = 0
 287       C1 = A
 288       C2 = CYCLE+A
 289       P1 = 0.1D+01-P
 290       UFLOW = D1MACH(1)
 291       EPS = EPSABS
 292       IF(EPSABS.GT.UFLOW/P1) EPS = EPSABS*P1
 293       EP = EPS
 294       FACT = 0.1D+01
 295       CORREC = 0.0D+00
 296       ABSERR = 0.0D+00
 297       ERRSUM = 0.0D+00
 298 C
 299 C           MAIN DO-LOOP
 300 C           ------------
 301 C
 302       DO 50 LST = 1,LIMLST
 303 C
 304 C           INTEGRATE OVER CURRENT SUBINTERVAL.
 305 C
 306         EPSA = EPS*FACT
 307         CALL DQAWOE(F,C1,C2,OMEGA,INTEGR,EPSA,0.0D+00,LIMIT,LST,MAXP1,
 308      1  RSLST(LST),ERLST(LST),NEV,IERLST(LST),LAST,ALIST,BLIST,RLIST,
 309      2  ELIST,IORD,NNLOG,MOMCOM,CHEBMO)
 310         NEVAL = NEVAL+NEV
 311         FACT = FACT*P
 312         ERRSUM = ERRSUM+ERLST(LST)
 313         DRL = 0.5D+02*ABS(RSLST(LST))
 314 C
 315 C           TEST ON ACCURACY WITH PARTIAL SUM
 316 C
 317         IF((ERRSUM+DRL).LE.EPSABS.AND.LST.GE.6) GO TO 80
 318         CORREC = MAX(CORREC,ERLST(LST))
 319         IF(IERLST(LST).NE.0) EPS = MAX(EP,CORREC*P1)
 320         IF(IERLST(LST).NE.0) IER = 7
 321         IF(IER.EQ.7.AND.(ERRSUM+DRL).LE.CORREC*0.1D+02.AND.
 322      1  LST.GT.5) GO TO 80
 323         NUMRL2 = NUMRL2+1
 324         IF(LST.GT.1) GO TO 20
 325         PSUM(1) = RSLST(1)
 326         GO TO 40
 327    20   PSUM(NUMRL2) = PSUM(LL)+RSLST(LST)
 328         IF(LST.EQ.2) GO TO 40
 329 C
 330 C           TEST ON MAXIMUM NUMBER OF SUBINTERVALS
 331 C
 332         IF(LST.EQ.LIMLST) IER = 1
 333 C
 334 C           PERFORM NEW EXTRAPOLATION
 335 C
 336         CALL DQELG(NUMRL2,PSUM,RESEPS,ABSEPS,RES3LA,NRES)
 337 C
 338 C           TEST WHETHER EXTRAPOLATED RESULT IS INFLUENCED BY ROUNDOFF
 339 C
 340         KTMIN = KTMIN+1
 341         IF(KTMIN.GE.15.AND.ABSERR.LE.0.1D-02*(ERRSUM+DRL)) IER = 4
 342         IF(ABSEPS.GT.ABSERR.AND.LST.NE.3) GO TO 30
 343         ABSERR = ABSEPS
 344         RESULT = RESEPS
 345         KTMIN = 0
 346 C
 347 C           IF IER IS NOT 0, CHECK WHETHER DIRECT RESULT (PARTIAL SUM)
 348 C           OR EXTRAPOLATED RESULT YIELDS THE BEST INTEGRAL
 349 C           APPROXIMATION
 350 C
 351         IF((ABSERR+0.1D+02*CORREC).LE.EPSABS.OR.
 352      1  (ABSERR.LE.EPSABS.AND.0.1D+02*CORREC.GE.EPSABS)) GO TO 60
 353    30   IF(IER.NE.0.AND.IER.NE.7) GO TO 60
 354    40   LL = NUMRL2
 355         C1 = C2
 356         C2 = C2+CYCLE
 357    50 CONTINUE
 358 C
 359 C         SET FINAL RESULT AND ERROR ESTIMATE
 360 C         -----------------------------------
 361 C
 362    60 ABSERR = ABSERR+0.1D+02*CORREC
 363       IF(IER.EQ.0) GO TO 999
 364       IF(RESULT.NE.0.0D+00.AND.PSUM(NUMRL2).NE.0.0D+00) GO TO 70
 365       IF(ABSERR.GT.ERRSUM) GO TO 80
 366       IF(PSUM(NUMRL2).EQ.0.0D+00) GO TO 999
 367    70 IF(ABSERR/ABS(RESULT).GT.(ERRSUM+DRL)/ABS(PSUM(NUMRL2)))
 368      1  GO TO 80
 369       IF(IER.GE.1.AND.IER.NE.7) ABSERR = ABSERR+DRL
 370       GO TO 999
 371    80 RESULT = PSUM(NUMRL2)
 372       ABSERR = ERRSUM+DRL
 373   999 RETURN
 374       END