GNU Octave  9.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
qagp.f
Go to the documentation of this file.
1  subroutine qagp(f,a,b,npts2,points,epsabs,epsrel,result,abserr,
2  * neval,ier,leniw,lenw,last,iwork,work)
3 c***begin prologue qagp
4 c***date written 800101 (yymmdd)
5 c***revision date 830518 (yymmdd)
6 c***category no. h2a2a1
7 c***keywords automatic integrator, general-purpose,
8 c singularities at user specified points,
9 c extrapolation, globally adaptive
10 c***author piessens,robert,appl. math. & progr. div - k.u.leuven
11 c de doncker,elise,appl. math. & progr. div. - k.u.leuven
12 c***purpose the routine calculates an approximation result to a given
13 c definite integral i = integral of f over (a,b),
14 c hopefully satisfying following claim for accuracy
15 c break points of the integration interval, where local
16 c difficulties of the integrand may occur(e.g. singularities,
17 c discontinuities), are provided by the user.
18 c***description
19 c
20 c computation of a definite integral
21 c standard fortran subroutine
22 c real version
23 c
24 c parameters
25 c on entry
26 c f - subroutine f(x,ierr,result) defining the integrand
27 c function f(x). the actual name for f needs to be
28 c declared e x t e r n a l in the driver program.
29 c
30 c a - real
31 c lower limit of integration
32 c
33 c b - real
34 c upper limit of integration
35 c
36 c npts2 - integer
37 c number equal to two more than the number of
38 c user-supplied break points within the integration
39 c range, npts.ge.2.
40 c if npts2.lt.2, the routine will end with ier = 6.
41 c
42 c points - real
43 c vector of dimension npts2, the first (npts2-2)
44 c elements of which are the user provided break
45 c points. if these points do not constitute an
46 c ascending sequence there will be an automatic
47 c sorting.
48 c
49 c epsabs - real
50 c absolute accuracy requested
51 c epsrel - real
52 c relative accuracy requested
53 c if epsabs.le.0
54 c and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
55 c the routine will end with ier = 6.
56 c
57 c on return
58 c result - real
59 c approximation to the integral
60 c
61 c abserr - real
62 c estimate of the modulus of the absolute error,
63 c which should equal or exceed abs(i-result)
64 c
65 c neval - integer
66 c number of integrand evaluations
67 c
68 c ier - integer
69 c ier = 0 normal and reliable termination of the
70 c routine. it is assumed that the requested
71 c accuracy has been achieved.
72 c ier.gt.0 abnormal termination of the routine.
73 c the estimates for integral and error are
74 c less reliable. it is assumed that the
75 c requested accuracy has not been achieved.
76 c error messages
77 c ier = 1 maximum number of subdivisions allowed
78 c has been achieved. one can allow more
79 c subdivisions by increasing the value of
80 c limit (and taking the according dimension
81 c adjustments into account). however, if
82 c this yields no improvement it is advised
83 c to analyze the integrand in order to
84 c determine the integration difficulties. if
85 c the position of a local difficulty can be
86 c determined (i.e. singularity,
87 c discontinuity within the interval), it
88 c should be supplied to the routine as an
89 c element of the vector points. if necessary
90 c an appropriate special-purpose integrator
91 c must be used, which is designed for
92 c handling the type of difficulty involved.
93 c = 2 the occurrence of roundoff error is
94 c detected, which prevents the requested
95 c tolerance from being achieved.
96 c the error may be under-estimated.
97 c = 3 extremely bad integrand behaviour occurs
98 c at some points of the integration
99 c interval.
100 c = 4 the algorithm does not converge.
101 c roundoff error is detected in the
102 c extrapolation table.
103 c it is presumed that the requested
104 c tolerance cannot be achieved, and that
105 c the returned result is the best which
106 c can be obtained.
107 c = 5 the integral is probably divergent, or
108 c slowly convergent. it must be noted that
109 c divergence can occur with any other value
110 c of ier.gt.0.
111 c = 6 the input is invalid because
112 c npts2.lt.2 or
113 c break points are specified outside
114 c the integration range or
115 c (epsabs.le.0 and
116 c epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
117 c result, abserr, neval, last are set to
118 c zero. exept when leniw or lenw or npts2 is
119 c invalid, iwork(1), iwork(limit+1),
120 c work(limit*2+1) and work(limit*3+1)
121 c are set to zero.
122 c work(1) is set to a and work(limit+1)
123 c to b (where limit = (leniw-npts2)/2).
124 c
125 c dimensioning parameters
126 c leniw - integer
127 c dimensioning parameter for iwork
128 c leniw determines limit = (leniw-npts2)/2,
129 c which is the maximum number of subintervals in the
130 c partition of the given integration interval (a,b),
131 c leniw.ge.(3*npts2-2).
132 c if leniw.lt.(3*npts2-2), the routine will end with
133 c ier = 6.
134 c
135 c lenw - integer
136 c dimensioning parameter for work
137 c lenw must be at least leniw*2-npts2.
138 c if lenw.lt.leniw*2-npts2, the routine will end
139 c with ier = 6.
140 c
141 c last - integer
142 c on return, last equals the number of subintervals
143 c produced in the subdivision process, which
144 c determines the number of significant elements
145 c actually in the work arrays.
146 c
147 c work arrays
148 c iwork - integer
149 c vector of dimension at least leniw. on return,
150 c the first k elements of which contain
151 c pointers to the error estimates over the
152 c subintervals, such that work(limit*3+iwork(1)),...,
153 c work(limit*3+iwork(k)) form a decreasing
154 c sequence, with k = last if last.le.(limit/2+2), and
155 c k = limit+1-last otherwise
156 c iwork(limit+1), ...,iwork(limit+last) contain the
157 c subdivision levels of the subintervals, i.e.
158 c if (aa,bb) is a subinterval of (p1,p2)
159 c where p1 as well as p2 is a user-provided
160 c break point or integration limit, then (aa,bb) has
161 c level l if abs(bb-aa) = abs(p2-p1)*2**(-l),
162 c iwork(limit*2+1), ..., iwork(limit*2+npts2) have
163 c no significance for the user,
164 c note that limit = (leniw-npts2)/2.
165 c
166 c work - real
167 c vector of dimension at least lenw
168 c on return
169 c work(1), ..., work(last) contain the left
170 c end points of the subintervals in the
171 c partition of (a,b),
172 c work(limit+1), ..., work(limit+last) contain
173 c the right end points,
174 c work(limit*2+1), ..., work(limit*2+last) contain
175 c the integral approximations over the subintervals,
176 c work(limit*3+1), ..., work(limit*3+last)
177 c contain the corresponding error estimates,
178 c work(limit*4+1), ..., work(limit*4+npts2)
179 c contain the integration limits and the
180 c break points sorted in an ascending sequence.
181 c note that limit = (leniw-npts2)/2.
182 c
183 c***references (none)
184 c***routines called qagpe,xerror
185 c***end prologue qagp
186 c
187  real a,abserr,b,epsabs,epsrel,points,result,work
188  integer ier,iwork,leniw,lenw,limit,lvl,l1,l2,l3,neval,npts2
189 c
190  dimension iwork(leniw),points(npts2),work(lenw)
191 c
192  external f
193 c
194 c check validity of limit and lenw.
195 c
196 c***first executable statement qagp
197  ier = 6
198  neval = 0
199  last = 0
200  result = 0.0e+00
201  abserr = 0.0e+00
202  if(leniw.lt.(3*npts2-2).or.lenw.lt.(leniw*2-npts2).or.npts2.lt.2)
203  * go to 10
204 c
205 c prepare call for qagpe.
206 c
207  limit = (leniw-npts2)/2
208  l1 = limit+1
209  l2 = limit+l1
210  l3 = limit+l2
211  l4 = limit+l3
212 c
213  call qagpe(f,a,b,npts2,points,epsabs,epsrel,limit,result,abserr,
214  * neval,ier,work(1),work(l1),work(l2),work(l3),work(l4),
215  * iwork(1),iwork(l1),iwork(l2),last)
216 c
217 c call error handler if necessary.
218 c
219  lvl = 0
220 10 if(ier.eq.6) lvl = 1
221  if(ier.ne.0) call xerror('abnormal return from qagp',26,ier,lvl)
222  return
223  end
subroutine qagp(f, a, b, npts2, points, epsabs, epsrel, result, abserr, neval, ier, leniw, lenw, last, iwork, work)
Definition: qagp.f:3
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 xerror(MESSG, NMESSG, NERR, LEVEL)
Definition: xerror.f:2