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