src/numerical/slatec/fortran/dqawoe.f

   1 *DECK DQAWOE
   2       SUBROUTINE DQAWOE (F, A, B, OMEGA, INTEGR, EPSABS, EPSREL, LIMIT,
   3      +   ICALL, MAXP1, RESULT, ABSERR, NEVAL, IER, LAST, ALIST, BLIST,
   4      +   RLIST, ELIST, IORD, NNLOG, MOMCOM, CHEBMO)
   5 C***BEGIN PROLOGUE  DQAWOE
   6 C***PURPOSE  Calculate an approximation to a given definite integral
   7 C            I = Integral of F(X)*W(X) over (A,B), where
   8 C                     W(X) = COS(OMEGA*X)
   9 C                 or  W(X)=SIN(OMEGA*X),
  10 C            hopefully satisfying the following claim for accuracy
  11 C                 ABS(I-RESULT).LE.MAX(EPSABS,EPSREL*ABS(I)).
  12 C***LIBRARY   SLATEC (QUADPACK)
  13 C***CATEGORY  H2A2A1
  14 C***TYPE      DOUBLE PRECISION (QAWOE-S, DQAWOE-D)
  15 C***KEYWORDS  AUTOMATIC INTEGRATOR, CLENSHAW-CURTIS METHOD,
  16 C             EXTRAPOLATION, GLOBALLY ADAPTIVE,
  17 C             INTEGRAND WITH OSCILLATORY COS OR SIN FACTOR, QUADPACK,
  18 C             QUADRATURE, SPECIAL-PURPOSE
  19 C***AUTHOR  Piessens, Robert
  20 C             Applied Mathematics and Programming Division
  21 C             K. U. Leuven
  22 C           de Doncker, Elise
  23 C             Applied Mathematics and Programming Division
  24 C             K. U. Leuven
  25 C***DESCRIPTION
  26 C
  27 C        Computation of Oscillatory integrals
  28 C        Standard fortran subroutine
  29 C        Double precision version
  30 C
  31 C        PARAMETERS
  32 C         ON ENTRY
  33 C            F      - Double precision
  34 C                     Function subprogram defining the integrand
  35 C                     function F(X). The actual name for F needs to be
  36 C                     declared E X T E R N A L in the driver program.
  37 C
  38 C            A      - Double precision
  39 C                     Lower limit of integration
  40 C
  41 C            B      - Double precision
  42 C                     Upper limit of integration
  43 C
  44 C            OMEGA  - Double precision
  45 C                     Parameter in the integrand weight function
  46 C
  47 C            INTEGR - Integer
  48 C                     Indicates which of the WEIGHT functions is to be
  49 C                     used
  50 C                     INTEGR = 1      W(X) = COS(OMEGA*X)
  51 C                     INTEGR = 2      W(X) = SIN(OMEGA*X)
  52 C                     If INTEGR.NE.1 and INTEGR.NE.2, the routine
  53 C                     will end with IER = 6.
  54 C
  55 C            EPSABS - Double precision
  56 C                     Absolute accuracy requested
  57 C            EPSREL - Double precision
  58 C                     Relative accuracy requested
  59 C                     If  EPSABS.LE.0
  60 C                     and EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28),
  61 C                     the routine will end with IER = 6.
  62 C
  63 C            LIMIT  - Integer
  64 C                     Gives an upper bound on the number of subdivisions
  65 C                     in the partition of (A,B), LIMIT.GE.1.
  66 C
  67 C            ICALL  - Integer
  68 C                     If DQAWOE is to be used only once, ICALL must
  69 C                     be set to 1.  Assume that during this call, the
  70 C                     Chebyshev moments (for CLENSHAW-CURTIS integration
  71 C                     of degree 24) have been computed for intervals of
  72 C                     lengths (ABS(B-A))*2**(-L), L=0,1,2,...MOMCOM-1.
  73 C                     If ICALL.GT.1 this means that DQAWOE has been
  74 C                     called twice or more on intervals of the same
  75 C                     length ABS(B-A). The Chebyshev moments already
  76 C                     computed are then re-used in subsequent calls.
  77 C                     If ICALL.LT.1, the routine will end with IER = 6.
  78 C
  79 C            MAXP1  - Integer
  80 C                     Gives an upper bound on the number of Chebyshev
  81 C                     moments which can be stored, i.e. for the
  82 C                     intervals of lengths ABS(B-A)*2**(-L),
  83 C                     L=0,1, ..., MAXP1-2, MAXP1.GE.1.
  84 C                     If MAXP1.LT.1, the routine will end with IER = 6.
  85 C
  86 C         ON RETURN
  87 C            RESULT - Double precision
  88 C                     Approximation to the integral
  89 C
  90 C            ABSERR - Double precision
  91 C                     Estimate of the modulus of the absolute error,
  92 C                     which should equal or exceed ABS(I-RESULT)
  93 C
  94 C            NEVAL  - Integer
  95 C                     Number of integrand evaluations
  96 C
  97 C            IER    - Integer
  98 C                     IER = 0 Normal and reliable termination of the
  99 C                             routine. It is assumed that the
 100 C                             requested accuracy has been achieved.
 101 C                   - IER.GT.0 Abnormal termination of the routine.
 102 C                             The estimates for integral and error are
 103 C                             less reliable. It is assumed that the
 104 C                             requested accuracy has not been achieved.
 105 C            ERROR MESSAGES
 106 C                     IER = 1 Maximum number of subdivisions allowed
 107 C                             has been achieved. One can allow more
 108 C                             subdivisions by increasing the value of
 109 C                             LIMIT (and taking according dimension
 110 C                             adjustments into account). However, if
 111 C                             this yields no improvement it is advised
 112 C                             to analyze the integrand, in order to
 113 C                             determine the integration difficulties.
 114 C                             If the position of a local difficulty can
 115 C                             be determined (e.g. SINGULARITY,
 116 C                             DISCONTINUITY within the interval) one
 117 C                             will probably gain from splitting up the
 118 C                             interval at this point and calling the
 119 C                             integrator on the subranges. If possible,
 120 C                             an appropriate special-purpose integrator
 121 C                             should be used which is designed for
 122 C                             handling the type of difficulty involved.
 123 C                         = 2 The occurrence of roundoff error is
 124 C                             detected, which prevents the requested
 125 C                             tolerance from being achieved.
 126 C                             The error may be under-estimated.
 127 C                         = 3 Extremely bad integrand behaviour occurs
 128 C                             at some points of the integration
 129 C                             interval.
 130 C                         = 4 The algorithm does not converge.
 131 C                             Roundoff error is detected in the
 132 C                             extrapolation table.
 133 C                             It is presumed that the requested
 134 C                             tolerance cannot be achieved due to
 135 C                             roundoff in the extrapolation table,
 136 C                             and that the returned result is the
 137 C                             best which can be obtained.
 138 C                         = 5 The integral is probably divergent, or
 139 C                             slowly convergent. It must be noted that
 140 C                             divergence can occur with any other value
 141 C                             of IER.GT.0.
 142 C                         = 6 The input is invalid, because
 143 C                             (EPSABS.LE.0 and
 144 C                              EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28))
 145 C                             or (INTEGR.NE.1 and INTEGR.NE.2) or
 146 C                             ICALL.LT.1 or MAXP1.LT.1.
 147 C                             RESULT, ABSERR, NEVAL, LAST, RLIST(1),
 148 C                             ELIST(1), IORD(1) and NNLOG(1) are set
 149 C                             to ZERO. ALIST(1) and BLIST(1) are set
 150 C                             to A and B respectively.
 151 C
 152 C            LAST  -  Integer
 153 C                     On return, LAST equals the number of
 154 C                     subintervals produces in the subdivision
 155 C                     process, which determines the number of
 156 C                     significant elements actually in the
 157 C                     WORK ARRAYS.
 158 C            ALIST  - Double precision
 159 C                     Vector of dimension at least LIMIT, the first
 160 C                      LAST  elements of which are the left
 161 C                     end points of the subintervals in the partition
 162 C                     of the given integration range (A,B)
 163 C
 164 C            BLIST  - Double precision
 165 C                     Vector of dimension at least LIMIT, the first
 166 C                      LAST  elements of which are the right
 167 C                     end points of the subintervals in the partition
 168 C                     of the given integration range (A,B)
 169 C
 170 C            RLIST  - Double precision
 171 C                     Vector of dimension at least LIMIT, the first
 172 C                      LAST  elements of which are the integral
 173 C                     approximations on the subintervals
 174 C
 175 C            ELIST  - Double precision
 176 C                     Vector of dimension at least LIMIT, the first
 177 C                      LAST  elements of which are the moduli of the
 178 C                     absolute error estimates on the subintervals
 179 C
 180 C            IORD   - Integer
 181 C                     Vector of dimension at least LIMIT, the first K
 182 C                     elements of which are pointers to the error
 183 C                     estimates over the subintervals,
 184 C                     such that ELIST(IORD(1)), ...,
 185 C                     ELIST(IORD(K)) form a decreasing sequence, with
 186 C                     K = LAST if LAST.LE.(LIMIT/2+2), and
 187 C                     K = LIMIT+1-LAST otherwise.
 188 C
 189 C            NNLOG  - Integer
 190 C                     Vector of dimension at least LIMIT, containing the
 191 C                     subdivision levels of the subintervals, i.e.
 192 C                     IWORK(I) = L means that the subinterval
 193 C                     numbered I is of length ABS(B-A)*2**(1-L)
 194 C
 195 C         ON ENTRY AND RETURN
 196 C            MOMCOM - Integer
 197 C                     Indicating that the Chebyshev moments
 198 C                     have been computed for intervals of lengths
 199 C                     (ABS(B-A))*2**(-L), L=0,1,2, ..., MOMCOM-1,
 200 C                     MOMCOM.LT.MAXP1
 201 C
 202 C            CHEBMO - Double precision
 203 C                     Array of dimension (MAXP1,25) containing the
 204 C                     Chebyshev moments
 205 C
 206 C***REFERENCES  (NONE)
 207 C***ROUTINES CALLED  D1MACH, DQC25F, DQELG, DQPSRT
 208 C***REVISION HISTORY  (YYMMDD)
 209 C   800101  DATE WRITTEN
 210 C   890531  Changed all specific intrinsics to generic.  (WRB)
 211 C   890831  Modified array declarations.  (WRB)
 212 C   890831  REVISION DATE from Version 3.2
 213 C   891214  Prologue converted to Version 4.0 format.  (BAB)
 214 C***END PROLOGUE  DQAWOE
 215 C
 216       DOUBLE PRECISION A,ABSEPS,ABSERR,ALIST,AREA,AREA1,AREA12,AREA2,A1,
 217      1  A2,B,BLIST,B1,B2,CHEBMO,CORREC,DEFAB1,DEFAB2,DEFABS,
 218      2  DOMEGA,D1MACH,DRES,ELIST,EPMACH,EPSABS,EPSREL,ERLARG,ERLAST,
 219      3  ERRBND,ERRMAX,ERROR1,ERRO12,ERROR2,ERRSUM,ERTEST,F,OFLOW,
 220      4  OMEGA,RESABS,RESEPS,RESULT,RES3LA,RLIST,RLIST2,SMALL,UFLOW,WIDTH
 221       INTEGER ICALL,ID,IER,IERRO,INTEGR,IORD,IROFF1,IROFF2,IROFF3,
 222      1  JUPBND,K,KSGN,KTMIN,LAST,LIMIT,MAXERR,MAXP1,MOMCOM,NEV,NEVAL,
 223      2  NNLOG,NRES,NRMAX,NRMOM,NUMRL2
 224       LOGICAL EXTRAP,NOEXT,EXTALL
 225 C
 226       DIMENSION ALIST(*),BLIST(*),RLIST(*),ELIST(*),
 227      1  IORD(*),RLIST2(52),RES3LA(3),CHEBMO(MAXP1,25),NNLOG(*)
 228 C
 229       EXTERNAL F
 230 C
 231 C            THE DIMENSION OF RLIST2 IS DETERMINED BY  THE VALUE OF
 232 C            LIMEXP IN SUBROUTINE DQELG (RLIST2 SHOULD BE OF
 233 C            DIMENSION (LIMEXP+2) AT LEAST).
 234 C
 235 C            LIST OF MAJOR VARIABLES
 236 C            -----------------------
 237 C
 238 C           ALIST     - LIST OF LEFT END POINTS OF ALL SUBINTERVALS
 239 C                       CONSIDERED UP TO NOW
 240 C           BLIST     - LIST OF RIGHT END POINTS OF ALL SUBINTERVALS
 241 C                       CONSIDERED UP TO NOW
 242 C           RLIST(I)  - APPROXIMATION TO THE INTEGRAL OVER
 243 C                       (ALIST(I),BLIST(I))
 244 C           RLIST2    - ARRAY OF DIMENSION AT LEAST LIMEXP+2
 245 C                       CONTAINING THE PART OF THE EPSILON TABLE
 246 C                       WHICH IS STILL NEEDED FOR FURTHER COMPUTATIONS
 247 C           ELIST(I)  - ERROR ESTIMATE APPLYING TO RLIST(I)
 248 C           MAXERR    - POINTER TO THE INTERVAL WITH LARGEST
 249 C                       ERROR ESTIMATE
 250 C           ERRMAX    - ELIST(MAXERR)
 251 C           ERLAST    - ERROR ON THE INTERVAL CURRENTLY SUBDIVIDED
 252 C           AREA      - SUM OF THE INTEGRALS OVER THE SUBINTERVALS
 253 C           ERRSUM    - SUM OF THE ERRORS OVER THE SUBINTERVALS
 254 C           ERRBND    - REQUESTED ACCURACY MAX(EPSABS,EPSREL*
 255 C                       ABS(RESULT))
 256 C           *****1    - VARIABLE FOR THE LEFT SUBINTERVAL
 257 C           *****2    - VARIABLE FOR THE RIGHT SUBINTERVAL
 258 C           LAST      - INDEX FOR SUBDIVISION
 259 C           NRES      - NUMBER OF CALLS TO THE EXTRAPOLATION ROUTINE
 260 C           NUMRL2    - NUMBER OF ELEMENTS IN RLIST2. IF AN APPROPRIATE
 261 C                       APPROXIMATION TO THE COMPOUNDED INTEGRAL HAS
 262 C                       BEEN OBTAINED IT IS PUT IN RLIST2(NUMRL2) AFTER
 263 C                       NUMRL2 HAS BEEN INCREASED BY ONE
 264 C           SMALL     - LENGTH OF THE SMALLEST INTERVAL CONSIDERED
 265 C                       UP TO NOW, MULTIPLIED BY 1.5
 266 C           ERLARG    - SUM OF THE ERRORS OVER THE INTERVALS LARGER
 267 C                       THAN THE SMALLEST INTERVAL CONSIDERED UP TO NOW
 268 C           EXTRAP    - LOGICAL VARIABLE DENOTING THAT THE ROUTINE IS
 269 C                       ATTEMPTING TO PERFORM EXTRAPOLATION, I.E. BEFORE
 270 C                       SUBDIVIDING THE SMALLEST INTERVAL WE TRY TO
 271 C                       DECREASE THE VALUE OF ERLARG
 272 C           NOEXT     - LOGICAL VARIABLE DENOTING THAT EXTRAPOLATION
 273 C                       IS NO LONGER ALLOWED (TRUE  VALUE)
 274 C
 275 C            MACHINE DEPENDENT CONSTANTS
 276 C            ---------------------------
 277 C
 278 C           EPMACH IS THE LARGEST RELATIVE SPACING.
 279 C           UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
 280 C           OFLOW IS THE LARGEST POSITIVE MAGNITUDE.
 281 C
 282 C***FIRST EXECUTABLE STATEMENT  DQAWOE
 283       EPMACH = D1MACH(4)
 284 C
 285 C         TEST ON VALIDITY OF PARAMETERS
 286 C         ------------------------------
 287 C
 288       IER = 0
 289       NEVAL = 0
 290       LAST = 0
 291       RESULT = 0.0D+00
 292       ABSERR = 0.0D+00
 293       ALIST(1) = A
 294       BLIST(1) = B
 295       RLIST(1) = 0.0D+00
 296       ELIST(1) = 0.0D+00
 297       IORD(1) = 0
 298       NNLOG(1) = 0
 299       IF((INTEGR.NE.1.AND.INTEGR.NE.2).OR.(EPSABS.LE.0.0D+00.AND.
 300      1  EPSREL.LT.MAX(0.5D+02*EPMACH,0.5D-28)).OR.ICALL.LT.1.OR.
 301      2  MAXP1.LT.1) IER = 6
 302       IF(IER.EQ.6) GO TO 999
 303 C
 304 C           FIRST APPROXIMATION TO THE INTEGRAL
 305 C           -----------------------------------
 306 C
 307       DOMEGA = ABS(OMEGA)
 308       NRMOM = 0
 309       IF (ICALL.GT.1) GO TO 5
 310       MOMCOM = 0
 311     5 CALL DQC25F(F,A,B,DOMEGA,INTEGR,NRMOM,MAXP1,0,RESULT,ABSERR,
 312      1  NEVAL,DEFABS,RESABS,MOMCOM,CHEBMO)
 313 C
 314 C           TEST ON ACCURACY.
 315 C
 316       DRES = ABS(RESULT)
 317       ERRBND = MAX(EPSABS,EPSREL*DRES)
 318       RLIST(1) = RESULT
 319       ELIST(1) = ABSERR
 320       IORD(1) = 1
 321       IF(ABSERR.LE.0.1D+03*EPMACH*DEFABS.AND.ABSERR.GT.ERRBND) IER = 2
 322       IF(LIMIT.EQ.1) IER = 1
 323       IF(IER.NE.0.OR.ABSERR.LE.ERRBND) GO TO 200
 324 C
 325 C           INITIALIZATIONS
 326 C           ---------------
 327 C
 328       UFLOW = D1MACH(1)
 329       OFLOW = D1MACH(2)
 330       ERRMAX = ABSERR
 331       MAXERR = 1
 332       AREA = RESULT
 333       ERRSUM = ABSERR
 334       ABSERR = OFLOW
 335       NRMAX = 1
 336       EXTRAP = .FALSE.
 337       NOEXT = .FALSE.
 338       IERRO = 0
 339       IROFF1 = 0
 340       IROFF2 = 0
 341       IROFF3 = 0
 342       KTMIN = 0
 343       SMALL = ABS(B-A)*0.75D+00
 344       NRES = 0
 345       NUMRL2 = 0
 346       EXTALL = .FALSE.
 347       IF(0.5D+00*ABS(B-A)*DOMEGA.GT.0.2D+01) GO TO 10
 348       NUMRL2 = 1
 349       EXTALL = .TRUE.
 350       RLIST2(1) = RESULT
 351    10 IF(0.25D+00*ABS(B-A)*DOMEGA.LE.0.2D+01) EXTALL = .TRUE.
 352       KSGN = -1
 353       IF(DRES.GE.(0.1D+01-0.5D+02*EPMACH)*DEFABS) KSGN = 1
 354 C
 355 C           MAIN DO-LOOP
 356 C           ------------
 357 C
 358       DO 140 LAST = 2,LIMIT
 359 C
 360 C           BISECT THE SUBINTERVAL WITH THE NRMAX-TH LARGEST
 361 C           ERROR ESTIMATE.
 362 C
 363         NRMOM = NNLOG(MAXERR)+1
 364         A1 = ALIST(MAXERR)
 365         B1 = 0.5D+00*(ALIST(MAXERR)+BLIST(MAXERR))
 366         A2 = B1
 367         B2 = BLIST(MAXERR)
 368         ERLAST = ERRMAX
 369         CALL DQC25F(F,A1,B1,DOMEGA,INTEGR,NRMOM,MAXP1,0,
 370      1  AREA1,ERROR1,NEV,RESABS,DEFAB1,MOMCOM,CHEBMO)
 371         NEVAL = NEVAL+NEV
 372         CALL DQC25F(F,A2,B2,DOMEGA,INTEGR,NRMOM,MAXP1,1,
 373      1  AREA2,ERROR2,NEV,RESABS,DEFAB2,MOMCOM,CHEBMO)
 374         NEVAL = NEVAL+NEV
 375 C
 376 C           IMPROVE PREVIOUS APPROXIMATIONS TO INTEGRAL
 377 C           AND ERROR AND TEST FOR ACCURACY.
 378 C
 379         AREA12 = AREA1+AREA2
 380         ERRO12 = ERROR1+ERROR2
 381         ERRSUM = ERRSUM+ERRO12-ERRMAX
 382         AREA = AREA+AREA12-RLIST(MAXERR)
 383         IF(DEFAB1.EQ.ERROR1.OR.DEFAB2.EQ.ERROR2) GO TO 25
 384         IF(ABS(RLIST(MAXERR)-AREA12).GT.0.1D-04*ABS(AREA12)
 385      1  .OR.ERRO12.LT.0.99D+00*ERRMAX) GO TO 20
 386         IF(EXTRAP) IROFF2 = IROFF2+1
 387         IF(.NOT.EXTRAP) IROFF1 = IROFF1+1
 388    20   IF(LAST.GT.10.AND.ERRO12.GT.ERRMAX) IROFF3 = IROFF3+1
 389    25   RLIST(MAXERR) = AREA1
 390         RLIST(LAST) = AREA2
 391         NNLOG(MAXERR) = NRMOM
 392         NNLOG(LAST) = NRMOM
 393         ERRBND = MAX(EPSABS,EPSREL*ABS(AREA))
 394 C
 395 C           TEST FOR ROUNDOFF ERROR AND EVENTUALLY SET ERROR FLAG.
 396 C
 397         IF(IROFF1+IROFF2.GE.10.OR.IROFF3.GE.20) IER = 2
 398         IF(IROFF2.GE.5) IERRO = 3
 399 C
 400 C           SET ERROR FLAG IN THE CASE THAT THE NUMBER OF
 401 C           SUBINTERVALS EQUALS LIMIT.
 402 C
 403         IF(LAST.EQ.LIMIT) IER = 1
 404 C
 405 C           SET ERROR FLAG IN THE CASE OF BAD INTEGRAND BEHAVIOUR
 406 C           AT A POINT OF THE INTEGRATION RANGE.
 407 C
 408         IF(MAX(ABS(A1),ABS(B2)).LE.(0.1D+01+0.1D+03*EPMACH)
 409      1  *(ABS(A2)+0.1D+04*UFLOW)) IER = 4
 410 C
 411 C           APPEND THE NEWLY-CREATED INTERVALS TO THE LIST.
 412 C
 413         IF(ERROR2.GT.ERROR1) GO TO 30
 414         ALIST(LAST) = A2
 415         BLIST(MAXERR) = B1
 416         BLIST(LAST) = B2
 417         ELIST(MAXERR) = ERROR1
 418         ELIST(LAST) = ERROR2
 419         GO TO 40
 420    30   ALIST(MAXERR) = A2
 421         ALIST(LAST) = A1
 422         BLIST(LAST) = B1
 423         RLIST(MAXERR) = AREA2
 424         RLIST(LAST) = AREA1
 425         ELIST(MAXERR) = ERROR2
 426         ELIST(LAST) = ERROR1
 427 C
 428 C           CALL SUBROUTINE DQPSRT TO MAINTAIN THE DESCENDING ORDERING
 429 C           IN THE LIST OF ERROR ESTIMATES AND SELECT THE SUBINTERVAL
 430 C           WITH NRMAX-TH LARGEST ERROR ESTIMATE (TO BISECTED NEXT).
 431 C
 432    40   CALL DQPSRT(LIMIT,LAST,MAXERR,ERRMAX,ELIST,IORD,NRMAX)
 433 C ***JUMP OUT OF DO-LOOP
 434       IF(ERRSUM.LE.ERRBND) GO TO 170
 435       IF(IER.NE.0) GO TO 150
 436         IF(LAST.EQ.2.AND.EXTALL) GO TO 120
 437         IF(NOEXT) GO TO 140
 438         IF(.NOT.EXTALL) GO TO 50
 439         ERLARG = ERLARG-ERLAST
 440         IF(ABS(B1-A1).GT.SMALL) ERLARG = ERLARG+ERRO12
 441         IF(EXTRAP) GO TO 70
 442 C
 443 C           TEST WHETHER THE INTERVAL TO BE BISECTED NEXT IS THE
 444 C           SMALLEST INTERVAL.
 445 C
 446    50   WIDTH = ABS(BLIST(MAXERR)-ALIST(MAXERR))
 447         IF(WIDTH.GT.SMALL) GO TO 140
 448         IF(EXTALL) GO TO 60
 449 C
 450 C           TEST WHETHER WE CAN START WITH THE EXTRAPOLATION PROCEDURE
 451 C           (WE DO THIS IF WE INTEGRATE OVER THE NEXT INTERVAL WITH
 452 C           USE OF A GAUSS-KRONROD RULE - SEE SUBROUTINE DQC25F).
 453 C
 454         SMALL = SMALL*0.5D+00
 455         IF(0.25D+00*WIDTH*DOMEGA.GT.0.2D+01) GO TO 140
 456         EXTALL = .TRUE.
 457         GO TO 130
 458    60   EXTRAP = .TRUE.
 459         NRMAX = 2
 460    70   IF(IERRO.EQ.3.OR.ERLARG.LE.ERTEST) GO TO 90
 461 C
 462 C           THE SMALLEST INTERVAL HAS THE LARGEST ERROR.
 463 C           BEFORE BISECTING DECREASE THE SUM OF THE ERRORS OVER
 464 C           THE LARGER INTERVALS (ERLARG) AND PERFORM EXTRAPOLATION.
 465 C
 466         JUPBND = LAST
 467         IF (LAST.GT.(LIMIT/2+2)) JUPBND = LIMIT+3-LAST
 468         ID = NRMAX
 469         DO 80 K = ID,JUPBND
 470           MAXERR = IORD(NRMAX)
 471           ERRMAX = ELIST(MAXERR)
 472           IF(ABS(BLIST(MAXERR)-ALIST(MAXERR)).GT.SMALL) GO TO 140
 473           NRMAX = NRMAX+1
 474    80   CONTINUE
 475 C
 476 C           PERFORM EXTRAPOLATION.
 477 C
 478    90   NUMRL2 = NUMRL2+1
 479         RLIST2(NUMRL2) = AREA
 480         IF(NUMRL2.LT.3) GO TO 110
 481         CALL DQELG(NUMRL2,RLIST2,RESEPS,ABSEPS,RES3LA,NRES)
 482         KTMIN = KTMIN+1
 483         IF(KTMIN.GT.5.AND.ABSERR.LT.0.1D-02*ERRSUM) IER = 5
 484         IF(ABSEPS.GE.ABSERR) GO TO 100
 485         KTMIN = 0
 486         ABSERR = ABSEPS
 487         RESULT = RESEPS
 488         CORREC = ERLARG
 489         ERTEST = MAX(EPSABS,EPSREL*ABS(RESEPS))
 490 C ***JUMP OUT OF DO-LOOP
 491         IF(ABSERR.LE.ERTEST) GO TO 150
 492 C
 493 C           PREPARE BISECTION OF THE SMALLEST INTERVAL.
 494 C
 495   100   IF(NUMRL2.EQ.1) NOEXT = .TRUE.
 496         IF(IER.EQ.5) GO TO 150
 497   110   MAXERR = IORD(1)
 498         ERRMAX = ELIST(MAXERR)
 499         NRMAX = 1
 500         EXTRAP = .FALSE.
 501         SMALL = SMALL*0.5D+00
 502         ERLARG = ERRSUM
 503         GO TO 140
 504   120   SMALL = SMALL*0.5D+00
 505         NUMRL2 = NUMRL2+1
 506         RLIST2(NUMRL2) = AREA
 507   130   ERTEST = ERRBND
 508         ERLARG = ERRSUM
 509   140 CONTINUE
 510 C
 511 C           SET THE FINAL RESULT.
 512 C           ---------------------
 513 C
 514   150 IF(ABSERR.EQ.OFLOW.OR.NRES.EQ.0) GO TO 170
 515       IF(IER+IERRO.EQ.0) GO TO 165
 516       IF(IERRO.EQ.3) ABSERR = ABSERR+CORREC
 517       IF(IER.EQ.0) IER = 3
 518       IF(RESULT.NE.0.0D+00.AND.AREA.NE.0.0D+00) GO TO 160
 519       IF(ABSERR.GT.ERRSUM) GO TO 170
 520       IF(AREA.EQ.0.0D+00) GO TO 190
 521       GO TO 165
 522   160 IF(ABSERR/ABS(RESULT).GT.ERRSUM/ABS(AREA)) GO TO 170
 523 C
 524 C           TEST ON DIVERGENCE.
 525 C
 526   165 IF(KSGN.EQ.(-1).AND.MAX(ABS(RESULT),ABS(AREA)).LE.
 527      1 DEFABS*0.1D-01) GO TO 190
 528       IF(0.1D-01.GT.(RESULT/AREA).OR.(RESULT/AREA).GT.0.1D+03
 529      1 .OR.ERRSUM.GE.ABS(AREA)) IER = 6
 530       GO TO 190
 531 C
 532 C           COMPUTE GLOBAL INTEGRAL SUM.
 533 C
 534   170 RESULT = 0.0D+00
 535       DO 180 K=1,LAST
 536         RESULT = RESULT+RLIST(K)
 537   180 CONTINUE
 538       ABSERR = ERRSUM
 539   190 IF (IER.GT.2) IER=IER-1
 540   200 IF (INTEGR.EQ.2.AND.OMEGA.LT.0.0D+00) RESULT=-RESULT
 541   999 RETURN
 542       END