src/numerical/slatec/fortran/dqage.f

   1 *DECK DQAGE
   2       SUBROUTINE DQAGE (F, A, B, EPSABS, EPSREL, KEY, LIMIT, RESULT,
   3      +   ABSERR, NEVAL, IER, ALIST, BLIST, RLIST, ELIST, IORD, LAST)
   4 C***BEGIN PROLOGUE  DQAGE
   5 C***PURPOSE  The routine calculates an approximation result to a given
   6 C            definite integral   I = Integral of F over (A,B),
   7 C            hopefully satisfying following claim for accuracy
   8 C            ABS(I-RESLT).LE.MAX(EPSABS,EPSREL*ABS(I)).
   9 C***LIBRARY   SLATEC (QUADPACK)
  10 C***CATEGORY  H2A1A1
  11 C***TYPE      DOUBLE PRECISION (QAGE-S, DQAGE-D)
  12 C***KEYWORDS  AUTOMATIC INTEGRATOR, GAUSS-KRONROD RULES,
  13 C             GENERAL-PURPOSE, GLOBALLY ADAPTIVE, INTEGRAND EXAMINATOR,
  14 C             QUADPACK, QUADRATURE
  15 C***AUTHOR  Piessens, Robert
  16 C             Applied Mathematics and Programming Division
  17 C             K. U. Leuven
  18 C           de Doncker, Elise
  19 C             Applied Mathematics and Programming Division
  20 C             K. U. Leuven
  21 C***DESCRIPTION
  22 C
  23 C        Computation of a definite integral
  24 C        Standard fortran subroutine
  25 C        Double precision version
  26 C
  27 C        PARAMETERS
  28 C         ON ENTRY
  29 C            F      - Double precision
  30 C                     Function subprogram defining the integrand
  31 C                     function F(X). The actual name for F needs to be
  32 C                     declared E X T E R N A L in the driver program.
  33 C
  34 C            A      - Double precision
  35 C                     Lower limit of integration
  36 C
  37 C            B      - Double precision
  38 C                     Upper limit of integration
  39 C
  40 C            EPSABS - Double precision
  41 C                     Absolute accuracy requested
  42 C            EPSREL - Double precision
  43 C                     Relative accuracy requested
  44 C                     If  EPSABS.LE.0
  45 C                     and EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28),
  46 C                     the routine will end with IER = 6.
  47 C
  48 C            KEY    - Integer
  49 C                     Key for choice of local integration rule
  50 C                     A Gauss-Kronrod pair is used with
  51 C                          7 - 15 points if KEY.LT.2,
  52 C                         10 - 21 points if KEY = 2,
  53 C                         15 - 31 points if KEY = 3,
  54 C                         20 - 41 points if KEY = 4,
  55 C                         25 - 51 points if KEY = 5,
  56 C                         30 - 61 points if KEY.GT.5.
  57 C
  58 C            LIMIT  - Integer
  59 C                     Gives an upper bound on the number of subintervals
  60 C                     in the partition of (A,B), LIMIT.GE.1.
  61 C
  62 C         ON RETURN
  63 C            RESULT - Double precision
  64 C                     Approximation to the integral
  65 C
  66 C            ABSERR - Double precision
  67 C                     Estimate of the modulus of the absolute error,
  68 C                     which should equal or exceed ABS(I-RESULT)
  69 C
  70 C            NEVAL  - Integer
  71 C                     Number of integrand evaluations
  72 C
  73 C            IER    - Integer
  74 C                     IER = 0 Normal and reliable termination of the
  75 C                             routine. It is assumed that the requested
  76 C                             accuracy has been achieved.
  77 C                     IER.GT.0 Abnormal termination of the routine
  78 C                             The estimates for result and error are
  79 C                             less reliable. It is assumed that the
  80 C                             requested accuracy has not been achieved.
  81 C            ERROR MESSAGES
  82 C                     IER = 1 Maximum number of subdivisions allowed
  83 C                             has been achieved. One can allow more
  84 C                             subdivisions by increasing the value
  85 C                             of LIMIT.
  86 C                             However, if this yields no improvement it
  87 C                             is rather advised to analyze the integrand
  88 C                             in 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                         = 3 Extremely bad integrand behaviour occurs
 103 C                             at some points of the integration
 104 C                             interval.
 105 C                         = 6 The input is invalid, because
 106 C                             (EPSABS.LE.0 and
 107 C                              EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28),
 108 C                             RESULT, ABSERR, NEVAL, LAST, RLIST(1) ,
 109 C                             ELIST(1) and IORD(1) are set to zero.
 110 C                             ALIST(1) and BLIST(1) are set to A and B
 111 C                             respectively.
 112 C
 113 C            ALIST   - Double precision
 114 C                      Vector of dimension at least LIMIT, the first
 115 C                       LAST  elements of which are the left
 116 C                      end points of the subintervals in the partition
 117 C                      of the given integration range (A,B)
 118 C
 119 C            BLIST   - Double precision
 120 C                      Vector of dimension at least LIMIT, the first
 121 C                       LAST  elements of which are the right
 122 C                      end points of the subintervals in the partition
 123 C                      of the given integration range (A,B)
 124 C
 125 C            RLIST   - Double precision
 126 C                      Vector of dimension at least LIMIT, the first
 127 C                       LAST  elements of which are the
 128 C                      integral approximations on the subintervals
 129 C
 130 C            ELIST   - Double precision
 131 C                      Vector of dimension at least LIMIT, the first
 132 C                       LAST  elements of which are the moduli of the
 133 C                      absolute error estimates on the subintervals
 134 C
 135 C            IORD    - Integer
 136 C                      Vector of dimension at least LIMIT, the first K
 137 C                      elements of which are pointers to the
 138 C                      error estimates over the subintervals,
 139 C                      such that ELIST(IORD(1)), ...,
 140 C                      ELIST(IORD(K)) form a decreasing sequence,
 141 C                      with K = LAST if LAST.LE.(LIMIT/2+2), and
 142 C                      K = LIMIT+1-LAST otherwise
 143 C
 144 C            LAST    - Integer
 145 C                      Number of subintervals actually produced in the
 146 C                      subdivision process
 147 C
 148 C***REFERENCES  (NONE)
 149 C***ROUTINES CALLED  D1MACH, DQK15, DQK21, DQK31, DQK41, DQK51, DQK61,
 150 C                    DQPSRT
 151 C***REVISION HISTORY  (YYMMDD)
 152 C   800101  DATE WRITTEN
 153 C   890531  Changed all specific intrinsics to generic.  (WRB)
 154 C   890831  Modified array declarations.  (WRB)
 155 C   890831  REVISION DATE from Version 3.2
 156 C   891214  Prologue converted to Version 4.0 format.  (BAB)
 157 C***END PROLOGUE  DQAGE
 158 C
 159       DOUBLE PRECISION A,ABSERR,ALIST,AREA,AREA1,AREA12,AREA2,A1,A2,B,
 160      1  BLIST,B1,B2,DEFABS,DEFAB1,DEFAB2,D1MACH,ELIST,EPMACH,
 161      2  EPSABS,EPSREL,ERRBND,ERRMAX,ERROR1,ERROR2,ERRO12,ERRSUM,F,
 162      3  RESABS,RESULT,RLIST,UFLOW
 163       INTEGER IER,IORD,IROFF1,IROFF2,K,KEY,KEYF,LAST,LIMIT,MAXERR,NEVAL,
 164      1  NRMAX
 165 C
 166       DIMENSION ALIST(*),BLIST(*),ELIST(*),IORD(*),
 167      1  RLIST(*)
 168 C
 169       EXTERNAL F
 170 C
 171 C            LIST OF MAJOR VARIABLES
 172 C            -----------------------
 173 C
 174 C           ALIST     - LIST OF LEFT END POINTS OF ALL SUBINTERVALS
 175 C                       CONSIDERED UP TO NOW
 176 C           BLIST     - LIST OF RIGHT END POINTS OF ALL SUBINTERVALS
 177 C                       CONSIDERED UP TO NOW
 178 C           RLIST(I)  - APPROXIMATION TO THE INTEGRAL OVER
 179 C                      (ALIST(I),BLIST(I))
 180 C           ELIST(I)  - ERROR ESTIMATE APPLYING TO RLIST(I)
 181 C           MAXERR    - POINTER TO THE INTERVAL WITH LARGEST
 182 C                       ERROR ESTIMATE
 183 C           ERRMAX    - ELIST(MAXERR)
 184 C           AREA      - SUM OF THE INTEGRALS OVER THE SUBINTERVALS
 185 C           ERRSUM    - SUM OF THE ERRORS OVER THE SUBINTERVALS
 186 C           ERRBND    - REQUESTED ACCURACY MAX(EPSABS,EPSREL*
 187 C                       ABS(RESULT))
 188 C           *****1    - VARIABLE FOR THE LEFT SUBINTERVAL
 189 C           *****2    - VARIABLE FOR THE RIGHT SUBINTERVAL
 190 C           LAST      - INDEX FOR SUBDIVISION
 191 C
 192 C
 193 C           MACHINE DEPENDENT CONSTANTS
 194 C           ---------------------------
 195 C
 196 C           EPMACH  IS THE LARGEST RELATIVE SPACING.
 197 C           UFLOW  IS THE SMALLEST POSITIVE MAGNITUDE.
 198 C
 199 C***FIRST EXECUTABLE STATEMENT  DQAGE
 200       EPMACH = D1MACH(4)
 201       UFLOW = D1MACH(1)
 202 C
 203 C           TEST ON VALIDITY OF PARAMETERS
 204 C           ------------------------------
 205 C
 206       IER = 0
 207       NEVAL = 0
 208       LAST = 0
 209       RESULT = 0.0D+00
 210       ABSERR = 0.0D+00
 211       ALIST(1) = A
 212       BLIST(1) = B
 213       RLIST(1) = 0.0D+00
 214       ELIST(1) = 0.0D+00
 215       IORD(1) = 0
 216       IF(EPSABS.LE.0.0D+00.AND.
 217      1  EPSREL.LT.MAX(0.5D+02*EPMACH,0.5D-28)) IER = 6
 218       IF(IER.EQ.6) GO TO 999
 219 C
 220 C           FIRST APPROXIMATION TO THE INTEGRAL
 221 C           -----------------------------------
 222 C
 223       KEYF = KEY
 224       IF(KEY.LE.0) KEYF = 1
 225       IF(KEY.GE.7) KEYF = 6
 226       NEVAL = 0
 227       IF(KEYF.EQ.1) CALL DQK15(F,A,B,RESULT,ABSERR,DEFABS,RESABS)
 228       IF(KEYF.EQ.2) CALL DQK21(F,A,B,RESULT,ABSERR,DEFABS,RESABS)
 229       IF(KEYF.EQ.3) CALL DQK31(F,A,B,RESULT,ABSERR,DEFABS,RESABS)
 230       IF(KEYF.EQ.4) CALL DQK41(F,A,B,RESULT,ABSERR,DEFABS,RESABS)
 231       IF(KEYF.EQ.5) CALL DQK51(F,A,B,RESULT,ABSERR,DEFABS,RESABS)
 232       IF(KEYF.EQ.6) CALL DQK61(F,A,B,RESULT,ABSERR,DEFABS,RESABS)
 233       LAST = 1
 234       RLIST(1) = RESULT
 235       ELIST(1) = ABSERR
 236       IORD(1) = 1
 237 C
 238 C           TEST ON ACCURACY.
 239 C
 240       ERRBND = MAX(EPSABS,EPSREL*ABS(RESULT))
 241       IF(ABSERR.LE.0.5D+02*EPMACH*DEFABS.AND.ABSERR.GT.ERRBND) IER = 2
 242       IF(LIMIT.EQ.1) IER = 1
 243       IF(IER.NE.0.OR.(ABSERR.LE.ERRBND.AND.ABSERR.NE.RESABS)
 244      1  .OR.ABSERR.EQ.0.0D+00) GO TO 60
 245 C
 246 C           INITIALIZATION
 247 C           --------------
 248 C
 249 C
 250       ERRMAX = ABSERR
 251       MAXERR = 1
 252       AREA = RESULT
 253       ERRSUM = ABSERR
 254       NRMAX = 1
 255       IROFF1 = 0
 256       IROFF2 = 0
 257 C
 258 C           MAIN DO-LOOP
 259 C           ------------
 260 C
 261       DO 30 LAST = 2,LIMIT
 262 C
 263 C           BISECT THE SUBINTERVAL WITH THE LARGEST ERROR ESTIMATE.
 264 C
 265         A1 = ALIST(MAXERR)
 266         B1 = 0.5D+00*(ALIST(MAXERR)+BLIST(MAXERR))
 267         A2 = B1
 268         B2 = BLIST(MAXERR)
 269         IF(KEYF.EQ.1) CALL DQK15(F,A1,B1,AREA1,ERROR1,RESABS,DEFAB1)
 270         IF(KEYF.EQ.2) CALL DQK21(F,A1,B1,AREA1,ERROR1,RESABS,DEFAB1)
 271         IF(KEYF.EQ.3) CALL DQK31(F,A1,B1,AREA1,ERROR1,RESABS,DEFAB1)
 272         IF(KEYF.EQ.4) CALL DQK41(F,A1,B1,AREA1,ERROR1,RESABS,DEFAB1)
 273         IF(KEYF.EQ.5) CALL DQK51(F,A1,B1,AREA1,ERROR1,RESABS,DEFAB1)
 274         IF(KEYF.EQ.6) CALL DQK61(F,A1,B1,AREA1,ERROR1,RESABS,DEFAB1)
 275         IF(KEYF.EQ.1) CALL DQK15(F,A2,B2,AREA2,ERROR2,RESABS,DEFAB2)
 276         IF(KEYF.EQ.2) CALL DQK21(F,A2,B2,AREA2,ERROR2,RESABS,DEFAB2)
 277         IF(KEYF.EQ.3) CALL DQK31(F,A2,B2,AREA2,ERROR2,RESABS,DEFAB2)
 278         IF(KEYF.EQ.4) CALL DQK41(F,A2,B2,AREA2,ERROR2,RESABS,DEFAB2)
 279         IF(KEYF.EQ.5) CALL DQK51(F,A2,B2,AREA2,ERROR2,RESABS,DEFAB2)
 280         IF(KEYF.EQ.6) CALL DQK61(F,A2,B2,AREA2,ERROR2,RESABS,DEFAB2)
 281 C
 282 C           IMPROVE PREVIOUS APPROXIMATIONS TO INTEGRAL
 283 C           AND ERROR AND TEST FOR ACCURACY.
 284 C
 285         NEVAL = NEVAL+1
 286         AREA12 = AREA1+AREA2
 287         ERRO12 = ERROR1+ERROR2
 288         ERRSUM = ERRSUM+ERRO12-ERRMAX
 289         AREA = AREA+AREA12-RLIST(MAXERR)
 290         IF(DEFAB1.EQ.ERROR1.OR.DEFAB2.EQ.ERROR2) GO TO 5
 291         IF(ABS(RLIST(MAXERR)-AREA12).LE.0.1D-04*ABS(AREA12)
 292      1  .AND.ERRO12.GE.0.99D+00*ERRMAX) IROFF1 = IROFF1+1
 293         IF(LAST.GT.10.AND.ERRO12.GT.ERRMAX) IROFF2 = IROFF2+1
 294     5   RLIST(MAXERR) = AREA1
 295         RLIST(LAST) = AREA2
 296         ERRBND = MAX(EPSABS,EPSREL*ABS(AREA))
 297         IF(ERRSUM.LE.ERRBND) GO TO 8
 298 C
 299 C           TEST FOR ROUNDOFF ERROR AND EVENTUALLY SET ERROR FLAG.
 300 C
 301         IF(IROFF1.GE.6.OR.IROFF2.GE.20) IER = 2
 302 C
 303 C           SET ERROR FLAG IN THE CASE THAT THE NUMBER OF SUBINTERVALS
 304 C           EQUALS LIMIT.
 305 C
 306         IF(LAST.EQ.LIMIT) IER = 1
 307 C
 308 C           SET ERROR FLAG IN THE CASE OF BAD INTEGRAND BEHAVIOUR
 309 C           AT A POINT OF THE INTEGRATION RANGE.
 310 C
 311         IF(MAX(ABS(A1),ABS(B2)).LE.(0.1D+01+0.1D+03*
 312      1  EPMACH)*(ABS(A2)+0.1D+04*UFLOW)) IER = 3
 313 C
 314 C           APPEND THE NEWLY-CREATED INTERVALS TO THE LIST.
 315 C
 316     8   IF(ERROR2.GT.ERROR1) GO TO 10
 317         ALIST(LAST) = A2
 318         BLIST(MAXERR) = B1
 319         BLIST(LAST) = B2
 320         ELIST(MAXERR) = ERROR1
 321         ELIST(LAST) = ERROR2
 322         GO TO 20
 323    10   ALIST(MAXERR) = A2
 324         ALIST(LAST) = A1
 325         BLIST(LAST) = B1
 326         RLIST(MAXERR) = AREA2
 327         RLIST(LAST) = AREA1
 328         ELIST(MAXERR) = ERROR2
 329         ELIST(LAST) = ERROR1
 330 C
 331 C           CALL SUBROUTINE DQPSRT TO MAINTAIN THE DESCENDING ORDERING
 332 C           IN THE LIST OF ERROR ESTIMATES AND SELECT THE SUBINTERVAL
 333 C           WITH THE LARGEST ERROR ESTIMATE (TO BE BISECTED NEXT).
 334 C
 335    20   CALL DQPSRT(LIMIT,LAST,MAXERR,ERRMAX,ELIST,IORD,NRMAX)
 336 C ***JUMP OUT OF DO-LOOP
 337         IF(IER.NE.0.OR.ERRSUM.LE.ERRBND) GO TO 40
 338    30 CONTINUE
 339 C
 340 C           COMPUTE FINAL RESULT.
 341 C           ---------------------
 342 C
 343    40 RESULT = 0.0D+00
 344       DO 50 K=1,LAST
 345         RESULT = RESULT+RLIST(K)
 346    50 CONTINUE
 347       ABSERR = ERRSUM
 348    60 IF(KEYF.NE.1) NEVAL = (10*KEYF+1)*(2*NEVAL+1)
 349       IF(KEYF.EQ.1) NEVAL = 30*NEVAL+15
 350   999 RETURN
 351       END