GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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
Generated on Mon Dec 30 2013 03:04:47 for GNU Octave by   doxygen 1.8.1.2