00001 subroutine qagpe(f,a,b,npts2,points,epsabs,epsrel,limit,result, 00002 * abserr,neval,ier,alist,blist,rlist,elist,pts,iord,level,ndin, 00003 * last) 00004 c***begin prologue qagpe 00005 c***date written 800101 (yymmdd) 00006 c***revision date 830518 (yymmdd) 00007 c***category no. h2a2a1 00008 c***keywords automatic integrator, general-purpose, 00009 c singularities at user specified points, 00010 c extrapolation, globally adaptive. 00011 c***author piessens,robert ,appl. math. & progr. div. - k.u.leuven 00012 c de doncker,elise,appl. math. & progr. div. - k.u.leuven 00013 c***purpose the routine calculates an approximation result to a given 00014 c definite integral i = integral of f over (a,b),hopefully 00015 c satisfying following claim for accuracy abs(i-result).le. 00016 c max(epsabs,epsrel*abs(i)). break points of the integration 00017 c interval, where local difficulties of the integrand may 00018 c occur(e.g. singularities,discontinuities),provided by user. 00019 c***description 00020 c 00021 c computation of a definite integral 00022 c standard fortran subroutine 00023 c real version 00024 c 00025 c parameters 00026 c on entry 00027 c f - subroutine f(x,ierr,result) defining the integrand 00028 c function f(x). the actual name for f needs to be 00029 c declared e x t e r n a l in the driver program. 00030 c 00031 c a - real 00032 c lower limit of integration 00033 c 00034 c b - real 00035 c upper limit of integration 00036 c 00037 c npts2 - integer 00038 c number equal to two more than the number of 00039 c user-supplied break points within the integration 00040 c range, npts2.ge.2. 00041 c if npts2.lt.2, the routine will end with ier = 6. 00042 c 00043 c points - real 00044 c vector of dimension npts2, the first (npts2-2) 00045 c elements of which are the user provided break 00046 c points. if these points do not constitute an 00047 c ascending sequence there will be an automatic 00048 c sorting. 00049 c 00050 c epsabs - real 00051 c absolute accuracy requested 00052 c epsrel - real 00053 c relative accuracy requested 00054 c if epsabs.le.0 00055 c and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), 00056 c the routine will end with ier = 6. 00057 c 00058 c limit - integer 00059 c gives an upper bound on the number of subintervals 00060 c in the partition of (a,b), limit.ge.npts2 00061 c if limit.lt.npts2, the routine will end with 00062 c ier = 6. 00063 c 00064 c on return 00065 c result - real 00066 c approximation to the integral 00067 c 00068 c abserr - real 00069 c estimate of the modulus of the absolute error, 00070 c which should equal or exceed abs(i-result) 00071 c 00072 c neval - integer 00073 c number of integrand evaluations 00074 c 00075 c ier - integer 00076 c ier = 0 normal and reliable termination of the 00077 c routine. it is assumed that the requested 00078 c accuracy has been achieved. 00079 c ier.gt.0 abnormal termination of the routine. 00080 c the estimates for integral and error are 00081 c less reliable. it is assumed that the 00082 c requested accuracy has not been achieved. 00083 c error messages 00084 c ier = 1 maximum number of subdivisions allowed 00085 c has been achieved. one can allow more 00086 c subdivisions by increasing the value of 00087 c limit (and taking the according dimension 00088 c adjustments into account). however, if 00089 c this yields no improvement it is advised 00090 c to analyze the integrand in order to 00091 c determine the integration difficulties. if 00092 c the position of a local difficulty can be 00093 c determined (i.e. singularity, 00094 c discontinuity within the interval), it 00095 c should be supplied to the routine as an 00096 c element of the vector points. if necessary 00097 c an appropriate special-purpose integrator 00098 c must be used, which is designed for 00099 c handling the type of difficulty involved. 00100 c = 2 the occurrence of roundoff error is 00101 c detected, which prevents the requested 00102 c tolerance from being achieved. 00103 c the error may be under-estimated. 00104 c = 3 extremely bad integrand behaviour occurs 00105 c at some points of the integration 00106 c interval. 00107 c = 4 the algorithm does not converge. 00108 c roundoff error is detected in the 00109 c extrapolation table. it is presumed that 00110 c the requested tolerance cannot be 00111 c achieved, and that the returned result is 00112 c the best which can be obtained. 00113 c = 5 the integral is probably divergent, or 00114 c slowly convergent. it must be noted that 00115 c divergence can occur with any other value 00116 c of ier.gt.0. 00117 c = 6 the input is invalid because 00118 c npts2.lt.2 or 00119 c break points are specified outside 00120 c the integration range or 00121 c (epsabs.le.0 and 00122 c epsrel.lt.max(50*rel.mach.acc.,0.5d-28)) 00123 c or limit.lt.npts2. 00124 c result, abserr, neval, last, rlist(1), 00125 c and elist(1) are set to zero. alist(1) and 00126 c blist(1) are set to a and b respectively. 00127 c 00128 c alist - real 00129 c vector of dimension at least limit, the first 00130 c last elements of which are the left end points 00131 c of the subintervals in the partition of the given 00132 c integration range (a,b) 00133 c 00134 c blist - real 00135 c vector of dimension at least limit, the first 00136 c last elements of which are the right end points 00137 c of the subintervals in the partition of the given 00138 c integration range (a,b) 00139 c 00140 c rlist - real 00141 c vector of dimension at least limit, the first 00142 c last elements of which are the integral 00143 c approximations on the subintervals 00144 c 00145 c elist - real 00146 c vector of dimension at least limit, the first 00147 c last elements of which are the moduli of the 00148 c absolute error estimates on the subintervals 00149 c 00150 c pts - real 00151 c vector of dimension at least npts2, containing the 00152 c integration limits and the break points of the 00153 c interval in ascending sequence. 00154 c 00155 c level - integer 00156 c vector of dimension at least limit, containing the 00157 c subdivision levels of the subinterval, i.e. if 00158 c (aa,bb) is a subinterval of (p1,p2) where p1 as 00159 c well as p2 is a user-provided break point or 00160 c integration limit, then (aa,bb) has level l if 00161 c abs(bb-aa) = abs(p2-p1)*2**(-l). 00162 c 00163 c ndin - integer 00164 c vector of dimension at least npts2, after first 00165 c integration over the intervals (pts(i)),pts(i+1), 00166 c i = 0,1, ..., npts2-2, the error estimates over 00167 c some of the intervals may have been increased 00168 c artificially, in order to put their subdivision 00169 c forward. if this happens for the subinterval 00170 c numbered k, ndin(k) is put to 1, otherwise 00171 c ndin(k) = 0. 00172 c 00173 c iord - integer 00174 c vector of dimension at least limit, the first k 00175 c elements of which are pointers to the 00176 c error estimates over the subintervals, 00177 c such that elist(iord(1)), ..., elist(iord(k)) 00178 c form a decreasing sequence, with k = last 00179 c if last.le.(limit/2+2), and k = limit+1-last 00180 c otherwise 00181 c 00182 c last - integer 00183 c number of subintervals actually produced in the 00184 c subdivisions process 00185 c 00186 c***references (none) 00187 c***routines called qelg,qk21,qpsrt,r1mach 00188 c***end prologue qagpe 00189 real a,abseps,abserr,alist,area,area1,area12,area2,a1, 00190 * a2,b,blist,b1,b2,correc,defabs,defab1,defab2, 00191 * dres,r1mach,elist,epmach,epsabs,epsrel,erlarg,erlast,errbnd, 00192 * errmax,error1,erro12,error2,errsum,ertest,oflow,points,pts, 00193 * resa,resabs,reseps,result,res3la,rlist,rlist2,sign,temp, 00194 * uflow 00195 integer i,id,ier,ierro,ind1,ind2,iord,ip1,iroff1,iroff2, 00196 * iroff3,j,jlow,jupbnd,k,ksgn,ktmin,last,levcur,level,levmax, 00197 * limit,maxerr,ndin,neval,nint,nintp1,npts,npts2,nres, 00198 * nrmax,numrl2 00199 logical extrap,noext 00200 c 00201 c 00202 dimension alist(limit),blist(limit),elist(limit),iord(limit), 00203 * level(limit),ndin(npts2),points(npts2),pts(npts2),res3la(3), 00204 * rlist(limit),rlist2(52) 00205 c 00206 external f 00207 c 00208 c the dimension of rlist2 is determined by the value of 00209 c limexp in subroutine epsalg (rlist2 should be of dimension 00210 c (limexp+2) at least). 00211 c 00212 c 00213 c list of major variables 00214 c ----------------------- 00215 c 00216 c alist - list of left end points of all subintervals 00217 c considered up to now 00218 c blist - list of right end points of all subintervals 00219 c considered up to now 00220 c rlist(i) - approximation to the integral over 00221 c (alist(i),blist(i)) 00222 c rlist2 - array of dimension at least limexp+2 00223 c containing the part of the epsilon table which 00224 c is still needed for further computations 00225 c elist(i) - error estimate applying to rlist(i) 00226 c maxerr - pointer to the interval with largest error 00227 c estimate 00228 c errmax - elist(maxerr) 00229 c erlast - error on the interval currently subdivided 00230 c (before that subdivision has taken place) 00231 c area - sum of the integrals over the subintervals 00232 c errsum - sum of the errors over the subintervals 00233 c errbnd - requested accuracy max(epsabs,epsrel* 00234 c abs(result)) 00235 c *****1 - variable for the left subinterval 00236 c *****2 - variable for the right subinterval 00237 c last - index for subdivision 00238 c nres - number of calls to the extrapolation routine 00239 c numrl2 - number of elements in rlist2. if an 00240 c appropriate approximation to the compounded 00241 c integral has been obtained, it is put in 00242 c rlist2(numrl2) after numrl2 has been increased 00243 c by one. 00244 c erlarg - sum of the errors over the intervals larger 00245 c than the smallest interval considered up to now 00246 c extrap - logical variable denoting that the routine 00247 c is attempting to perform extrapolation. i.e. 00248 c before subdividing the smallest interval we 00249 c try to decrease the value of erlarg. 00250 c noext - logical variable denoting that extrapolation is 00251 c no longer allowed (true-value) 00252 c 00253 c machine dependent constants 00254 c --------------------------- 00255 c 00256 c epmach is the largest relative spacing. 00257 c uflow is the smallest positive magnitude. 00258 c oflow is the largest positive magnitude. 00259 c 00260 c***first executable statement qagpe 00261 epmach = r1mach(4) 00262 c 00263 c test on validity of parameters 00264 c ----------------------------- 00265 c 00266 ier = 0 00267 neval = 0 00268 last = 0 00269 result = 0.0e+00 00270 abserr = 0.0e+00 00271 alist(1) = a 00272 blist(1) = b 00273 rlist(1) = 0.0e+00 00274 elist(1) = 0.0e+00 00275 iord(1) = 0 00276 level(1) = 0 00277 npts = npts2-2 00278 if(npts2.lt.2.or.limit.le.npts.or.(epsabs.le.0.0e+00.and. 00279 * epsrel.lt.amax1(0.5e+02*epmach,0.5e-14))) ier = 6 00280 if(ier.eq.6) go to 210 00281 c 00282 c if any break points are provided, sort them into an 00283 c ascending sequence. 00284 c 00285 sign = 1.0e+00 00286 if(a.gt.b) sign = -1.0e+00 00287 pts(1) = amin1(a,b) 00288 if(npts.eq.0) go to 15 00289 do 10 i = 1,npts 00290 pts(i+1) = points(i) 00291 10 continue 00292 15 pts(npts+2) = amax1(a,b) 00293 nint = npts+1 00294 a1 = pts(1) 00295 if(npts.eq.0) go to 40 00296 nintp1 = nint+1 00297 do 20 i = 1,nint 00298 ip1 = i+1 00299 do 20 j = ip1,nintp1 00300 if(pts(i).le.pts(j)) go to 20 00301 temp = pts(i) 00302 pts(i) = pts(j) 00303 pts(j) = temp 00304 20 continue 00305 if(pts(1).ne.amin1(a,b).or.pts(nintp1).ne. 00306 * amax1(a,b)) ier = 6 00307 if(ier.eq.6) go to 999 00308 c 00309 c compute first integral and error approximations. 00310 c ------------------------------------------------ 00311 c 00312 40 resabs = 0.0e+00 00313 do 50 i = 1,nint 00314 b1 = pts(i+1) 00315 call qk21(f,a1,b1,area1,error1,defabs,resa,ier) 00316 if (ier.lt.0) return 00317 abserr = abserr+error1 00318 result = result+area1 00319 ndin(i) = 0 00320 if(error1.eq.resa.and.error1.ne.0.0e+00) ndin(i) = 1 00321 resabs = resabs+defabs 00322 level(i) = 0 00323 elist(i) = error1 00324 alist(i) = a1 00325 blist(i) = b1 00326 rlist(i) = area1 00327 iord(i) = i 00328 a1 = b1 00329 50 continue 00330 errsum = 0.0e+00 00331 do 55 i = 1,nint 00332 if(ndin(i).eq.1) elist(i) = abserr 00333 errsum = errsum+elist(i) 00334 55 continue 00335 c 00336 c test on accuracy. 00337 c 00338 last = nint 00339 neval = 21*nint 00340 dres = abs(result) 00341 errbnd = amax1(epsabs,epsrel*dres) 00342 if(abserr.le.0.1e+03*epmach*resabs.and.abserr.gt. 00343 * errbnd) ier = 2 00344 if(nint.eq.1) go to 80 00345 do 70 i = 1,npts 00346 jlow = i+1 00347 ind1 = iord(i) 00348 do 60 j = jlow,nint 00349 ind2 = iord(j) 00350 if(elist(ind1).gt.elist(ind2)) go to 60 00351 ind1 = ind2 00352 k = j 00353 60 continue 00354 if(ind1.eq.iord(i)) go to 70 00355 iord(k) = iord(i) 00356 iord(i) = ind1 00357 70 continue 00358 if(limit.lt.npts2) ier = 1 00359 80 if(ier.ne.0.or.abserr.le.errbnd) go to 999 00360 c 00361 c initialization 00362 c -------------- 00363 c 00364 rlist2(1) = result 00365 maxerr = iord(1) 00366 errmax = elist(maxerr) 00367 area = result 00368 nrmax = 1 00369 nres = 0 00370 numrl2 = 1 00371 ktmin = 0 00372 extrap = .false. 00373 noext = .false. 00374 erlarg = errsum 00375 ertest = errbnd 00376 levmax = 1 00377 iroff1 = 0 00378 iroff2 = 0 00379 iroff3 = 0 00380 ierro = 0 00381 uflow = r1mach(1) 00382 oflow = r1mach(2) 00383 abserr = oflow 00384 ksgn = -1 00385 if(dres.ge.(0.1e+01-0.5e+02*epmach)*resabs) ksgn = 1 00386 c 00387 c main do-loop 00388 c ------------ 00389 c 00390 do 160 last = npts2,limit 00391 c 00392 c bisect the subinterval with the nrmax-th largest 00393 c error estimate. 00394 c 00395 levcur = level(maxerr)+1 00396 a1 = alist(maxerr) 00397 b1 = 0.5e+00*(alist(maxerr)+blist(maxerr)) 00398 a2 = b1 00399 b2 = blist(maxerr) 00400 erlast = errmax 00401 call qk21(f,a1,b1,area1,error1,resa,defab1,ier) 00402 if (ier.lt.0) return 00403 call qk21(f,a2,b2,area2,error2,resa,defab2,ier) 00404 if (ier.lt.0) return 00405 c 00406 c improve previous approximations to integral 00407 c and error and test for accuracy. 00408 c 00409 neval = neval+42 00410 area12 = area1+area2 00411 erro12 = error1+error2 00412 errsum = errsum+erro12-errmax 00413 area = area+area12-rlist(maxerr) 00414 if(defab1.eq.error1.or.defab2.eq.error2) go to 95 00415 if(abs(rlist(maxerr)-area12).gt.0.1e-04*abs(area12) 00416 * .or.erro12.lt.0.99e+00*errmax) go to 90 00417 if(extrap) iroff2 = iroff2+1 00418 if(.not.extrap) iroff1 = iroff1+1 00419 90 if(last.gt.10.and.erro12.gt.errmax) iroff3 = iroff3+1 00420 95 level(maxerr) = levcur 00421 level(last) = levcur 00422 rlist(maxerr) = area1 00423 rlist(last) = area2 00424 errbnd = amax1(epsabs,epsrel*abs(area)) 00425 c 00426 c test for roundoff error and eventually 00427 c set error flag. 00428 c 00429 if(iroff1+iroff2.ge.10.or.iroff3.ge.20) ier = 2 00430 if(iroff2.ge.5) ierro = 3 00431 c 00432 c set error flag in the case that the number of 00433 c subintervals equals limit. 00434 c 00435 if(last.eq.limit) ier = 1 00436 c 00437 c set error flag in the case of bad integrand behaviour 00438 c at a point of the integration range 00439 c 00440 if(amax1(abs(a1),abs(b2)).le.(0.1e+01+0.1e+03*epmach)* 00441 * (abs(a2)+0.1e+04*uflow)) ier = 4 00442 c 00443 c append the newly-created intervals to the list. 00444 c 00445 if(error2.gt.error1) go to 100 00446 alist(last) = a2 00447 blist(maxerr) = b1 00448 blist(last) = b2 00449 elist(maxerr) = error1 00450 elist(last) = error2 00451 go to 110 00452 100 alist(maxerr) = a2 00453 alist(last) = a1 00454 blist(last) = b1 00455 rlist(maxerr) = area2 00456 rlist(last) = area1 00457 elist(maxerr) = error2 00458 elist(last) = error1 00459 c 00460 c call subroutine qpsrt to maintain the descending ordering 00461 c in the list of error estimates and select the 00462 c subinterval with nrmax-th largest error estimate (to be 00463 c bisected next). 00464 c 00465 110 call qpsrt(limit,last,maxerr,errmax,elist,iord,nrmax) 00466 c ***jump out of do-loop 00467 if(errsum.le.errbnd) go to 190 00468 c ***jump out of do-loop 00469 if(ier.ne.0) go to 170 00470 if(noext) go to 160 00471 erlarg = erlarg-erlast 00472 if(levcur+1.le.levmax) erlarg = erlarg+erro12 00473 if(extrap) go to 120 00474 c 00475 c test whether the interval to be bisected next is the 00476 c smallest interval. 00477 c 00478 if(level(maxerr)+1.le.levmax) go to 160 00479 extrap = .true. 00480 nrmax = 2 00481 120 if(ierro.eq.3.or.erlarg.le.ertest) go to 140 00482 c 00483 c the smallest interval has the largest error. 00484 c before bisecting decrease the sum of the errors 00485 c over the larger intervals (erlarg) and perform 00486 c extrapolation. 00487 c 00488 id = nrmax 00489 jupbnd = last 00490 if(last.gt.(2+limit/2)) jupbnd = limit+3-last 00491 do 130 k = id,jupbnd 00492 maxerr = iord(nrmax) 00493 errmax = elist(maxerr) 00494 c ***jump out of do-loop 00495 if(level(maxerr)+1.le.levmax) go to 160 00496 nrmax = nrmax+1 00497 130 continue 00498 c 00499 c perform extrapolation. 00500 c 00501 140 numrl2 = numrl2+1 00502 rlist2(numrl2) = area 00503 if(numrl2.le.2) go to 155 00504 call qelg(numrl2,rlist2,reseps,abseps,res3la,nres) 00505 ktmin = ktmin+1 00506 if(ktmin.gt.5.and.abserr.lt.0.1e-02*errsum) ier = 5 00507 if(abseps.ge.abserr) go to 150 00508 ktmin = 0 00509 abserr = abseps 00510 result = reseps 00511 correc = erlarg 00512 ertest = amax1(epsabs,epsrel*abs(reseps)) 00513 c ***jump out of do-loop 00514 if(abserr.lt.ertest) go to 170 00515 c 00516 c prepare bisection of the smallest interval. 00517 c 00518 150 if(numrl2.eq.1) noext = .true. 00519 if(ier.ge.5) go to 170 00520 155 maxerr = iord(1) 00521 errmax = elist(maxerr) 00522 nrmax = 1 00523 extrap = .false. 00524 levmax = levmax+1 00525 erlarg = errsum 00526 160 continue 00527 c 00528 c set the final result. 00529 c --------------------- 00530 c 00531 c 00532 170 if(abserr.eq.oflow) go to 190 00533 if((ier+ierro).eq.0) go to 180 00534 if(ierro.eq.3) abserr = abserr+correc 00535 if(ier.eq.0) ier = 3 00536 if(result.ne.0.0e+00.and.area.ne.0.0e+00)go to 175 00537 if(abserr.gt.errsum)go to 190 00538 if(area.eq.0.0e+00) go to 210 00539 go to 180 00540 175 if(abserr/abs(result).gt.errsum/abs(area))go to 190 00541 c 00542 c test on divergence. 00543 c 00544 180 if(ksgn.eq.(-1).and.amax1(abs(result),abs(area)).le. 00545 * resabs*0.1e-01) go to 210 00546 if(0.1e-01.gt.(result/area).or.(result/area).gt.0.1e+03.or. 00547 * errsum.gt.abs(area)) ier = 6 00548 go to 210 00549 c 00550 c compute global integral sum. 00551 c 00552 190 result = 0.0e+00 00553 do 200 k = 1,last 00554 result = result+rlist(k) 00555 200 continue 00556 abserr = errsum 00557 210 if(ier.gt.2) ier = ier - 1 00558 result = result*sign 00559 999 return 00560 end
1.7.1
