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