GNU Octave  4.0.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
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 xerror(MESSG, NMESSG, NERR, LEVEL)
Definition: xerror.f:1
std::string dimension(void) const
subroutine qagp(f, a, b, npts2, points, epsabs, epsrel, result, abserr, neval, ier, leniw, lenw, last, iwork, work)
Definition: qagp.f:1
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:1