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