src/numerical/slatec/fortran/dqawce.f

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