1 subroutine dqawoe (f,a,b,omega,integr,epsabs,epsrel,limit,icall, 2 * maxp1,result,abserr,neval,ier,last,alist,blist,rlist,elist,iord, 3 * nnlog,momcom,chebmo) 4c***begin prologue dqawoe 5c***date written 800101 (yymmdd) 6c***revision date 830518 (yymmdd) 7c***category no. h2a2a1 8c***keywords automatic integrator, special-purpose, 9c integrand with oscillatory cos or sin factor, 10c clenshaw-curtis method, (end point) singularities, 11c extrapolation, globally adaptive 12c***author piessens,robert,appl. math. & progr. div. - k.u.leuven 13c de doncker,elise,appl. math. & progr. div. - k.u.leuven 14c***purpose the routine calculates an approximation result to a given 15c definite integral 16c i = integral of f(x)*w(x) over (a,b) 17c where w(x) = cos(omega*x) or w(x)=sin(omega*x), 18c hopefully satisfying following claim for accuracy 19c abs(i-result).le.max(epsabs,epsrel*abs(i)). 20c***description 21c 22c computation of oscillatory integrals 23c standard fortran subroutine 24c double precision version 25c 26c parameters 27c on entry 28c f - double precision 29c function subprogram defining the integrand 30c function f(x). the actual name for f needs to be 31c declared e x t e r n a l in the driver program. 32c 33c a - double precision 34c lower limit of integration 35c 36c b - double precision 37c upper limit of integration 38c 39c omega - double precision 40c parameter in the integrand weight function 41c 42c integr - integer 43c indicates which of the weight functions is to be 44c used 45c integr = 1 w(x) = cos(omega*x) 46c integr = 2 w(x) = sin(omega*x) 47c if integr.ne.1 and integr.ne.2, the routine 48c will end with ier = 6. 49c 50c epsabs - double precision 51c absolute accuracy requested 52c epsrel - double precision 53c relative accuracy requested 54c if epsabs.le.0 55c and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), 56c the routine will end with ier = 6. 57c 58c limit - integer 59c gives an upper bound on the number of subdivisions 60c in the partition of (a,b), limit.ge.1. 61c 62c icall - integer 63c if dqawoe is to be used only once, icall must 64c be set to 1. assume that during this call, the 65c chebyshev moments (for clenshaw-curtis integration 66c of degree 24) have been computed for intervals of 67c lengths (abs(b-a))*2**(-l), l=0,1,2,...momcom-1. 68c if icall.gt.1 this means that dqawoe has been 69c called twice or more on intervals of the same 70c length abs(b-a). the chebyshev moments already 71c computed are then re-used in subsequent calls. 72c if icall.lt.1, the routine will end with ier = 6. 73c 74c maxp1 - integer 75c gives an upper bound on the number of chebyshev 76c moments which can be stored, i.e. for the 77c intervals of lengths abs(b-a)*2**(-l), 78c l=0,1, ..., maxp1-2, maxp1.ge.1. 79c if maxp1.lt.1, the routine will end with ier = 6. 80c 81c on return 82c result - double precision 83c approximation to the integral 84c 85c abserr - double precision 86c estimate of the modulus of the absolute error, 87c which should equal or exceed abs(i-result) 88c 89c neval - integer 90c number of integrand evaluations 91c 92c ier - integer 93c ier = 0 normal and reliable termination of the 94c routine. it is assumed that the 95c requested accuracy has been achieved. 96c - ier.gt.0 abnormal termination of the routine. 97c the estimates for integral and error are 98c less reliable. it is assumed that the 99c requested accuracy has not been achieved. 100c error messages 101c ier = 1 maximum number of subdivisions allowed 102c has been achieved. one can allow more 103c subdivisions by increasing the value of 104c limit (and taking according dimension 105c adjustments into account). however, if 106c this yields no improvement it is advised 107c to analyze the integrand, in order to 108c determine the integration difficulties. 109c if the position of a local difficulty can 110c be determined (e.g. singularity, 111c discontinuity within the interval) one 112c will probably gain from splitting up the 113c interval at this point and calling the 114c integrator on the subranges. if possible, 115c an appropriate special-purpose integrator 116c should be used which is designed for 117c handling the type of difficulty involved. 118c = 2 the occurrence of roundoff error is 119c detected, which prevents the requested 120c tolerance from being achieved. 121c the error may be under-estimated. 122c = 3 extremely bad integrand behaviour occurs 123c at some points of the integration 124c interval. 125c = 4 the algorithm does not converge. 126c roundoff error is detected in the 127c extrapolation table. 128c it is presumed that the requested 129c tolerance cannot be achieved due to 130c roundoff in the extrapolation table, 131c and that the returned result is the 132c best which can be obtained. 133c = 5 the integral is probably divergent, or 134c slowly convergent. it must be noted that 135c divergence can occur with any other value 136c of ier.gt.0. 137c = 6 the input is invalid, because 138c (epsabs.le.0 and 139c epsrel.lt.max(50*rel.mach.acc.,0.5d-28)) 140c or (integr.ne.1 and integr.ne.2) or 141c icall.lt.1 or maxp1.lt.1. 142c result, abserr, neval, last, rlist(1), 143c elist(1), iord(1) and nnlog(1) are set 144c to zero. alist(1) and blist(1) are set 145c to a and b respectively. 146c 147c last - integer 148c on return, last equals the number of 149c subintervals produces in the subdivision 150c process, which determines the number of 151c significant elements actually in the 152c work arrays. 153c alist - double precision 154c vector of dimension at least limit, the first 155c last elements of which are the left 156c end points of the subintervals in the partition 157c of the given integration range (a,b) 158c 159c blist - double precision 160c vector of dimension at least limit, the first 161c last elements of which are the right 162c end points of the subintervals in the partition 163c of the given integration range (a,b) 164c 165c rlist - double precision 166c vector of dimension at least limit, the first 167c last elements of which are the integral 168c approximations on the subintervals 169c 170c elist - double precision 171c vector of dimension at least limit, the first 172c last elements of which are the moduli of the 173c absolute error estimates on the subintervals 174c 175c iord - integer 176c vector of dimension at least limit, the first k 177c elements of which are pointers to the error 178c estimates over the subintervals, 179c such that elist(iord(1)), ..., 180c elist(iord(k)) form a decreasing sequence, with 181c k = last if last.le.(limit/2+2), and 182c k = limit+1-last otherwise. 183c 184c nnlog - integer 185c vector of dimension at least limit, containing the 186c subdivision levels of the subintervals, i.e. 187c iwork(i) = l means that the subinterval 188c numbered i is of length abs(b-a)*2**(1-l) 189c 190c on entry and return 191c momcom - integer 192c indicating that the chebyshev moments 193c have been computed for intervals of lengths 194c (abs(b-a))*2**(-l), l=0,1,2, ..., momcom-1, 195c momcom.lt.maxp1 196c 197c chebmo - double precision 198c array of dimension (maxp1,25) containing the 199c chebyshev moments 200c 201c***references (none) 202c***routines called d1mach,dqc25f,dqelg,dqpsrt 203c***end prologue dqawoe 204c 205 double precision a,abseps,abserr,alist,area,area1,area12,area2,a1, 206 * a2,b,blist,b1,b2,chebmo,correc,dabs,defab1,defab2,defabs,dmax1, 207 * domega,d1mach,dres,elist,epmach,epsabs,epsrel,erlarg,erlast, 208 * errbnd,errmax,error1,erro12,error2,errsum,ertest,f,oflow, 209 * omega,resabs,reseps,result,res3la,rlist,rlist2,small,uflow,width 210 integer icall,id,ier,ierro,integr,iord,iroff1,iroff2,iroff3, 211 * jupbnd,k,ksgn,ktmin,last,limit,maxerr,maxp1,momcom,nev,neval, 212 * nnlog,nres,nrmax,nrmom,numrl2 213 logical extrap,noext,extall 214c 215 dimension alist(limit),blist(limit),rlist(limit),elist(limit), 216 * iord(limit),rlist2(52),res3la(3),chebmo(maxp1,25),nnlog(limit) 217c 218 external f 219c 220c the dimension of rlist2 is determined by the value of 221c limexp in subroutine dqelg (rlist2 should be of 222c dimension (limexp+2) at least). 223c 224c list of major variables 225c ----------------------- 226c 227c alist - list of left end points of all subintervals 228c considered up to now 229c blist - list of right end points of all subintervals 230c considered up to now 231c rlist(i) - approximation to the integral over 232c (alist(i),blist(i)) 233c rlist2 - array of dimension at least limexp+2 234c containing the part of the epsilon table 235c which is still needed for further computations 236c elist(i) - error estimate applying to rlist(i) 237c maxerr - pointer to the interval with largest 238c error estimate 239c errmax - elist(maxerr) 240c erlast - error on the interval currently subdivided 241c area - sum of the integrals over the subintervals 242c errsum - sum of the errors over the subintervals 243c errbnd - requested accuracy max(epsabs,epsrel* 244c abs(result)) 245c *****1 - variable for the left subinterval 246c *****2 - variable for the right subinterval 247c last - index for subdivision 248c nres - number of calls to the extrapolation routine 249c numrl2 - number of elements in rlist2. if an appropriate 250c approximation to the compounded integral has 251c been obtained it is put in rlist2(numrl2) after 252c numrl2 has been increased by one 253c small - length of the smallest interval considered 254c up to now, multiplied by 1.5 255c erlarg - sum of the errors over the intervals larger 256c than the smallest interval considered up to now 257c extrap - logical variable denoting that the routine is 258c attempting to perform extrapolation, i.e. before 259c subdividing the smallest interval we try to 260c decrease the value of erlarg 261c noext - logical variable denoting that extrapolation 262c is no longer allowed (true value) 263c 264c machine dependent constants 265c --------------------------- 266c 267c epmach is the largest relative spacing. 268c uflow is the smallest positive magnitude. 269c oflow is the largest positive magnitude. 270c 271c***first executable statement dqawoe 272 epmach = d1mach(4) 273c 274c test on validity of parameters 275c ------------------------------ 276c 277 ier = 0 278 neval = 0 279 last = 0 280 result = 0.0d+00 281 abserr = 0.0d+00 282 alist(1) = a 283 blist(1) = b 284 rlist(1) = 0.0d+00 285 elist(1) = 0.0d+00 286 iord(1) = 0 287 nnlog(1) = 0 288 if((integr.ne.1.and.integr.ne.2).or.(epsabs.le.0.0d+00.and. 289 * epsrel.lt.dmax1(0.5d+02*epmach,0.5d-28)).or.icall.lt.1.or. 290 * maxp1.lt.1) ier = 6 291 if(ier.eq.6) go to 999 292c 293c first approximation to the integral 294c ----------------------------------- 295c 296 domega = dabs(omega) 297 nrmom = 0 298 if (icall.gt.1) go to 5 299 momcom = 0 300 5 call dqc25f(f,a,b,domega,integr,nrmom,maxp1,0,result,abserr, 301 * neval,defabs,resabs,momcom,chebmo) 302c 303c test on accuracy. 304c 305 dres = dabs(result) 306 errbnd = dmax1(epsabs,epsrel*dres) 307 rlist(1) = result 308 elist(1) = abserr 309 iord(1) = 1 310 if(abserr.le.0.1d+03*epmach*defabs.and.abserr.gt.errbnd) ier = 2 311 if(limit.eq.1) ier = 1 312 if(ier.ne.0.or.abserr.le.errbnd) go to 200 313c 314c initializations 315c --------------- 316c 317 uflow = d1mach(1) 318 oflow = d1mach(2) 319 errmax = abserr 320 maxerr = 1 321 area = result 322 errsum = abserr 323 abserr = oflow 324 nrmax = 1 325 extrap = .false. 326 noext = .false. 327 ierro = 0 328 iroff1 = 0 329 iroff2 = 0 330 iroff3 = 0 331 ktmin = 0 332 small = dabs(b-a)*0.75d+00 333 nres = 0 334 numrl2 = 0 335 extall = .false. 336 if(0.5d+00*dabs(b-a)*domega.gt.0.2d+01) go to 10 337 numrl2 = 1 338 extall = .true. 339 rlist2(1) = result 340 10 if(0.25d+00*dabs(b-a)*domega.le.0.2d+01) extall = .true. 341 ksgn = -1 342 if(dres.ge.(0.1d+01-0.5d+02*epmach)*defabs) ksgn = 1 343c 344c main do-loop 345c ------------ 346c 347 do 140 last = 2,limit 348c 349c bisect the subinterval with the nrmax-th largest 350c error estimate. 351c 352 nrmom = nnlog(maxerr)+1 353 a1 = alist(maxerr) 354 b1 = 0.5d+00*(alist(maxerr)+blist(maxerr)) 355 a2 = b1 356 b2 = blist(maxerr) 357 erlast = errmax 358 call dqc25f(f,a1,b1,domega,integr,nrmom,maxp1,0, 359 * area1,error1,nev,resabs,defab1,momcom,chebmo) 360 neval = neval+nev 361 call dqc25f(f,a2,b2,domega,integr,nrmom,maxp1,1, 362 * area2,error2,nev,resabs,defab2,momcom,chebmo) 363 neval = neval+nev 364c 365c improve previous approximations to integral 366c and error and test for accuracy. 367c 368 area12 = area1+area2 369 erro12 = error1+error2 370 errsum = errsum+erro12-errmax 371 area = area+area12-rlist(maxerr) 372 if(defab1.eq.error1.or.defab2.eq.error2) go to 25 373 if(dabs(rlist(maxerr)-area12).gt.0.1d-04*dabs(area12) 374 * .or.erro12.lt.0.99d+00*errmax) go to 20 375 if(extrap) iroff2 = iroff2+1 376 if(.not.extrap) iroff1 = iroff1+1 377 20 if(last.gt.10.and.erro12.gt.errmax) iroff3 = iroff3+1 378 25 rlist(maxerr) = area1 379 rlist(last) = area2 380 nnlog(maxerr) = nrmom 381 nnlog(last) = nrmom 382 errbnd = dmax1(epsabs,epsrel*dabs(area)) 383c 384c test for roundoff error and eventually set error flag. 385c 386 if(iroff1+iroff2.ge.10.or.iroff3.ge.20) ier = 2 387 if(iroff2.ge.5) ierro = 3 388c 389c set error flag in the case that the number of 390c subintervals equals limit. 391c 392 if(last.eq.limit) ier = 1 393c 394c set error flag in the case of bad integrand behaviour 395c at a point of the integration range. 396c 397 if(dmax1(dabs(a1),dabs(b2)).le.(0.1d+01+0.1d+03*epmach) 398 * *(dabs(a2)+0.1d+04*uflow)) ier = 4 399c 400c append the newly-created intervals to the list. 401c 402 if(error2.gt.error1) go to 30 403 alist(last) = a2 404 blist(maxerr) = b1 405 blist(last) = b2 406 elist(maxerr) = error1 407 elist(last) = error2 408 go to 40 409 30 alist(maxerr) = a2 410 alist(last) = a1 411 blist(last) = b1 412 rlist(maxerr) = area2 413 rlist(last) = area1 414 elist(maxerr) = error2 415 elist(last) = error1 416c 417c call subroutine dqpsrt to maintain the descending ordering 418c in the list of error estimates and select the subinterval 419c with nrmax-th largest error estimate (to bisected next). 420c 421 40 call dqpsrt(limit,last,maxerr,errmax,elist,iord,nrmax) 422c ***jump out of do-loop 423 if(errsum.le.errbnd) go to 170 424 if(ier.ne.0) go to 150 425 if(last.eq.2.and.extall) go to 120 426 if(noext) go to 140 427 if(.not.extall) go to 50 428 erlarg = erlarg-erlast 429 if(dabs(b1-a1).gt.small) erlarg = erlarg+erro12 430 if(extrap) go to 70 431c 432c test whether the interval to be bisected next is the 433c smallest interval. 434c 435 50 width = dabs(blist(maxerr)-alist(maxerr)) 436 if(width.gt.small) go to 140 437 if(extall) go to 60 438c 439c test whether we can start with the extrapolation procedure 440c (we do this if we integrate over the next interval with 441c use of a gauss-kronrod rule - see subroutine dqc25f). 442c 443 small = small*0.5d+00 444 if(0.25d+00*width*domega.gt.0.2d+01) go to 140 445 extall = .true. 446 go to 130 447 60 extrap = .true. 448 nrmax = 2 449 70 if(ierro.eq.3.or.erlarg.le.ertest) go to 90 450c 451c the smallest interval has the largest error. 452c before bisecting decrease the sum of the errors over 453c the larger intervals (erlarg) and perform extrapolation. 454c 455 jupbnd = last 456 if (last.gt.(limit/2+2)) jupbnd = limit+3-last 457 id = nrmax 458 do 80 k = id,jupbnd 459 maxerr = iord(nrmax) 460 errmax = elist(maxerr) 461 if(dabs(blist(maxerr)-alist(maxerr)).gt.small) go to 140 462 nrmax = nrmax+1 463 80 continue 464c 465c perform extrapolation. 466c 467 90 numrl2 = numrl2+1 468 rlist2(numrl2) = area 469 if(numrl2.lt.3) go to 110 470 call dqelg(numrl2,rlist2,reseps,abseps,res3la,nres) 471 ktmin = ktmin+1 472 if(ktmin.gt.5.and.abserr.lt.0.1d-02*errsum) ier = 5 473 if(abseps.ge.abserr) go to 100 474 ktmin = 0 475 abserr = abseps 476 result = reseps 477 correc = erlarg 478 ertest = dmax1(epsabs,epsrel*dabs(reseps)) 479c ***jump out of do-loop 480 if(abserr.le.ertest) go to 150 481c 482c prepare bisection of the smallest interval. 483c 484 100 if(numrl2.eq.1) noext = .true. 485 if(ier.eq.5) go to 150 486 110 maxerr = iord(1) 487 errmax = elist(maxerr) 488 nrmax = 1 489 extrap = .false. 490 small = small*0.5d+00 491 erlarg = errsum 492 go to 140 493 120 small = small*0.5d+00 494 numrl2 = numrl2+1 495 rlist2(numrl2) = area 496 130 ertest = errbnd 497 erlarg = errsum 498 140 continue 499c 500c set the final result. 501c --------------------- 502c 503 150 if(abserr.eq.oflow.or.nres.eq.0) go to 170 504 if(ier+ierro.eq.0) go to 165 505 if(ierro.eq.3) abserr = abserr+correc 506 if(ier.eq.0) ier = 3 507 if(result.ne.0.0d+00.and.area.ne.0.0d+00) go to 160 508 if(abserr.gt.errsum) go to 170 509 if(area.eq.0.0d+00) go to 190 510 go to 165 511 160 if(abserr/dabs(result).gt.errsum/dabs(area)) go to 170 512c 513c test on divergence. 514c 515 165 if(ksgn.eq.(-1).and.dmax1(dabs(result),dabs(area)).le. 516 * defabs*0.1d-01) go to 190 517 if(0.1d-01.gt.(result/area).or.(result/area).gt.0.1d+03 518 * .or.errsum.ge.dabs(area)) ier = 6 519 go to 190 520c 521c compute global integral sum. 522c 523 170 result = 0.0d+00 524 do 180 k=1,last 525 result = result+rlist(k) 526 180 continue 527 abserr = errsum 528 190 if (ier.gt.2) ier=ier-1 529 200 if (integr.eq.2.and.omega.lt.0.0d+00) result=-result 530 999 return 531 end 532