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