1 subroutine dqawse(f,a,b,alfa,beta,integr,epsabs,epsrel,limit, 2 * result,abserr,neval,ier,alist,blist,rlist,elist,iord,last) 3c***begin prologue dqawse 4c***date written 800101 (yymmdd) 5c***revision date 830518 (yymmdd) 6c***category no. h2a2a1 7c***keywords automatic integrator, special-purpose, 8c algebraico-logarithmic end point singularities, 9c clenshaw-curtis method 10c***author piessens,robert,appl. math. & progr. div. - k.u.leuven 11c de doncker,elise,appl. math. & progr. div. - k.u.leuven 12c***purpose the routine calculates an approximation result to a given 13c definite integral i = integral of f*w over (a,b), 14c (where w shows a singular behaviour at the end points, 15c see parameter integr). 16c hopefully satisfying following claim for accuracy 17c abs(i-result).le.max(epsabs,epsrel*abs(i)). 18c***description 19c 20c integration of functions having algebraico-logarithmic 21c end point singularities 22c standard fortran subroutine 23c double precision version 24c 25c parameters 26c on entry 27c f - double precision 28c function subprogram defining the integrand 29c function f(x). the actual name for f needs to be 30c declared e x t e r n a l in the driver program. 31c 32c a - double precision 33c lower limit of integration 34c 35c b - double precision 36c upper limit of integration, b.gt.a 37c if b.le.a, the routine will end with ier = 6. 38c 39c alfa - double precision 40c parameter in the weight function, alfa.gt.(-1) 41c if alfa.le.(-1), the routine will end with 42c ier = 6. 43c 44c beta - double precision 45c parameter in the weight function, beta.gt.(-1) 46c if beta.le.(-1), the routine will end with 47c ier = 6. 48c 49c integr - integer 50c indicates which weight function is to be used 51c = 1 (x-a)**alfa*(b-x)**beta 52c = 2 (x-a)**alfa*(b-x)**beta*log(x-a) 53c = 3 (x-a)**alfa*(b-x)**beta*log(b-x) 54c = 4 (x-a)**alfa*(b-x)**beta*log(x-a)*log(b-x) 55c if integr.lt.1 or integr.gt.4, the routine 56c will end with ier = 6. 57c 58c epsabs - double precision 59c absolute accuracy requested 60c epsrel - double precision 61c relative accuracy requested 62c if epsabs.le.0 63c and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), 64c the routine will end with ier = 6. 65c 66c limit - integer 67c gives an upper bound on the number of subintervals 68c in the partition of (a,b), limit.ge.2 69c if limit.lt.2, the routine will end with ier = 6. 70c 71c on return 72c result - double precision 73c approximation to the integral 74c 75c abserr - double precision 76c estimate of the modulus of the absolute error, 77c which should equal or exceed abs(i-result) 78c 79c neval - integer 80c number of integrand evaluations 81c 82c ier - integer 83c ier = 0 normal and reliable termination of the 84c routine. it is assumed that the requested 85c accuracy has been achieved. 86c ier.gt.0 abnormal termination of the routine 87c the estimates for the integral and error 88c are less reliable. it is assumed that the 89c requested accuracy has not been achieved. 90c error messages 91c = 1 maximum number of subdivisions allowed 92c has been achieved. one can allow more 93c subdivisions by increasing the value of 94c limit. however, if this yields no 95c improvement, it is advised to analyze the 96c integrand in order to determine the 97c integration difficulties which prevent the 98c requested tolerance from being achieved. 99c in case of a jump discontinuity or a local 100c singularity of algebraico-logarithmic type 101c at one or more interior points of the 102c integration range, one should proceed by 103c splitting up the interval at these 104c points and calling the integrator on the 105c subranges. 106c = 2 the occurrence of roundoff error is 107c detected, which prevents the requested 108c tolerance from being achieved. 109c = 3 extremely bad integrand behaviour occurs 110c at some points of the integration 111c interval. 112c = 6 the input is invalid, because 113c b.le.a or alfa.le.(-1) or beta.le.(-1), or 114c integr.lt.1 or integr.gt.4, or 115c (epsabs.le.0 and 116c epsrel.lt.max(50*rel.mach.acc.,0.5d-28), 117c or limit.lt.2. 118c result, abserr, neval, rlist(1), elist(1), 119c iord(1) and last are set to zero. alist(1) 120c and blist(1) are set to a and b 121c respectively. 122c 123c alist - double precision 124c vector of dimension at least limit, the first 125c last elements of which are the left 126c end points of the subintervals in the partition 127c of the given integration range (a,b) 128c 129c blist - double precision 130c vector of dimension at least limit, the first 131c last elements of which are the right 132c end points of the subintervals in the partition 133c of the given integration range (a,b) 134c 135c rlist - double precision 136c vector of dimension at least limit,the first 137c last elements of which are the integral 138c approximations on the subintervals 139c 140c elist - double precision 141c vector of dimension at least limit, the first 142c last elements of which are the moduli of the 143c absolute error estimates on the subintervals 144c 145c iord - integer 146c vector of dimension at least limit, the first k 147c of which are pointers to the error 148c estimates over the subintervals, so that 149c elist(iord(1)), ..., elist(iord(k)) with k = last 150c if last.le.(limit/2+2), and k = limit+1-last 151c otherwise form a decreasing sequence 152c 153c last - integer 154c number of subintervals actually produced in 155c the subdivision process 156c 157c***references (none) 158c***routines called d1mach,dqc25s,dqmomo,dqpsrt 159c***end prologue dqawse 160c 161 double precision a,abserr,alfa,alist,area,area1,area12,area2,a1, 162 * a2,b,beta,blist,b1,b2,centre,dabs,dmax1,d1mach,elist,epmach, 163 * epsabs,epsrel,errbnd,errmax,error1,erro12,error2,errsum,f, 164 * resas1,resas2,result,rg,rh,ri,rj,rlist,uflow 165 integer ier,integr,iord,iroff1,iroff2,k,last,limit,maxerr,nev, 166 * neval,nrmax 167c 168 external f 169c 170 dimension alist(limit),blist(limit),rlist(limit),elist(limit), 171 * iord(limit),ri(25),rj(25),rh(25),rg(25) 172c 173c list of major variables 174c ----------------------- 175c 176c alist - list of left end points of all subintervals 177c considered up to now 178c blist - list of right end points of all subintervals 179c considered up to now 180c rlist(i) - approximation to the integral over 181c (alist(i),blist(i)) 182c elist(i) - error estimate applying to rlist(i) 183c maxerr - pointer to the interval with largest 184c error estimate 185c errmax - elist(maxerr) 186c area - sum of the integrals over the subintervals 187c errsum - sum of the errors over the subintervals 188c errbnd - requested accuracy max(epsabs,epsrel* 189c abs(result)) 190c *****1 - variable for the left subinterval 191c *****2 - variable for the right subinterval 192c last - index for subdivision 193c 194c 195c machine dependent constants 196c --------------------------- 197c 198c epmach is the largest relative spacing. 199c uflow is the smallest positive magnitude. 200c 201c***first executable statement dqawse 202 epmach = d1mach(4) 203 uflow = d1mach(1) 204c 205c test on validity of parameters 206c ------------------------------ 207c 208 ier = 6 209 neval = 0 210 last = 0 211 rlist(1) = 0.0d+00 212 elist(1) = 0.0d+00 213 iord(1) = 0 214 result = 0.0d+00 215 abserr = 0.0d+00 216 if(b.le.a.or.(epsabs.eq.0.0d+00.and. 217 * epsrel.lt.dmax1(0.5d+02*epmach,0.5d-28)).or.alfa.le.(-0.1d+01). 218 * or.beta.le.(-0.1d+01).or.integr.lt.1.or.integr.gt.4.or. 219 * limit.lt.2) go to 999 220 ier = 0 221c 222c compute the modified chebyshev moments. 223c 224 call dqmomo(alfa,beta,ri,rj,rg,rh,integr) 225c 226c integrate over the intervals (a,(a+b)/2) and ((a+b)/2,b). 227c 228 centre = 0.5d+00*(b+a) 229 call dqc25s(f,a,b,a,centre,alfa,beta,ri,rj,rg,rh,area1, 230 * error1,resas1,integr,nev) 231 neval = nev 232 call dqc25s(f,a,b,centre,b,alfa,beta,ri,rj,rg,rh,area2, 233 * error2,resas2,integr,nev) 234 last = 2 235 neval = neval+nev 236 result = area1+area2 237 abserr = error1+error2 238c 239c test on accuracy. 240c 241 errbnd = dmax1(epsabs,epsrel*dabs(result)) 242c 243c initialization 244c -------------- 245c 246 if(error2.gt.error1) go to 10 247 alist(1) = a 248 alist(2) = centre 249 blist(1) = centre 250 blist(2) = b 251 rlist(1) = area1 252 rlist(2) = area2 253 elist(1) = error1 254 elist(2) = error2 255 go to 20 256 10 alist(1) = centre 257 alist(2) = a 258 blist(1) = b 259 blist(2) = centre 260 rlist(1) = area2 261 rlist(2) = area1 262 elist(1) = error2 263 elist(2) = error1 264 20 iord(1) = 1 265 iord(2) = 2 266 if(limit.eq.2) ier = 1 267 if(abserr.le.errbnd.or.ier.eq.1) go to 999 268 errmax = elist(1) 269 maxerr = 1 270 nrmax = 1 271 area = result 272 errsum = abserr 273 iroff1 = 0 274 iroff2 = 0 275c 276c main do-loop 277c ------------ 278c 279 do 60 last = 3,limit 280c 281c bisect the subinterval with largest error estimate. 282c 283 a1 = alist(maxerr) 284 b1 = 0.5d+00*(alist(maxerr)+blist(maxerr)) 285 a2 = b1 286 b2 = blist(maxerr) 287c 288 call dqc25s(f,a,b,a1,b1,alfa,beta,ri,rj,rg,rh,area1, 289 * error1,resas1,integr,nev) 290 neval = neval+nev 291 call dqc25s(f,a,b,a2,b2,alfa,beta,ri,rj,rg,rh,area2, 292 * error2,resas2,integr,nev) 293 neval = neval+nev 294c 295c improve previous approximations integral and error 296c and test for accuracy. 297c 298 area12 = area1+area2 299 erro12 = error1+error2 300 errsum = errsum+erro12-errmax 301 area = area+area12-rlist(maxerr) 302 if(a.eq.a1.or.b.eq.b2) go to 30 303 if(resas1.eq.error1.or.resas2.eq.error2) go to 30 304c 305c test for roundoff error. 306c 307 if(dabs(rlist(maxerr)-area12).lt.0.1d-04*dabs(area12) 308 * .and.erro12.ge.0.99d+00*errmax) iroff1 = iroff1+1 309 if(last.gt.10.and.erro12.gt.errmax) iroff2 = iroff2+1 310 30 rlist(maxerr) = area1 311 rlist(last) = area2 312c 313c test on accuracy. 314c 315 errbnd = dmax1(epsabs,epsrel*dabs(area)) 316 if(errsum.le.errbnd) go to 35 317c 318c set error flag in the case that the number of interval 319c bisections exceeds limit. 320c 321 if(last.eq.limit) ier = 1 322c 323c 324c set error flag in the case of roundoff error. 325c 326 if(iroff1.ge.6.or.iroff2.ge.20) ier = 2 327c 328c set error flag in the case of bad integrand behaviour 329c at interior points of integration range. 330c 331 if(dmax1(dabs(a1),dabs(b2)).le.(0.1d+01+0.1d+03*epmach)* 332 * (dabs(a2)+0.1d+04*uflow)) ier = 3 333c 334c append the newly-created intervals to the list. 335c 336 35 if(error2.gt.error1) go to 40 337 alist(last) = a2 338 blist(maxerr) = b1 339 blist(last) = b2 340 elist(maxerr) = error1 341 elist(last) = error2 342 go to 50 343 40 alist(maxerr) = a2 344 alist(last) = a1 345 blist(last) = b1 346 rlist(maxerr) = area2 347 rlist(last) = area1 348 elist(maxerr) = error2 349 elist(last) = error1 350c 351c call subroutine dqpsrt to maintain the descending ordering 352c in the list of error estimates and select the subinterval 353c with largest error estimate (to be bisected next). 354c 355 50 call dqpsrt(limit,last,maxerr,errmax,elist,iord,nrmax) 356c ***jump out of do-loop 357 if (ier.ne.0.or.errsum.le.errbnd) go to 70 358 60 continue 359c 360c compute final result. 361c --------------------- 362c 363 70 result = 0.0d+00 364 do 80 k=1,last 365 result = result+rlist(k) 366 80 continue 367 abserr = errsum 368 999 return 369 end 370