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
qagi.f
Go to the documentation of this file.
1  subroutine qagi(f,bound,inf,epsabs,epsrel,result,abserr,neval,
2  * ier,limit,lenw,last,iwork,work)
3 c***begin prologue qagi
4 c***date written 800101 (yymmdd)
5 c***revision date 830518 (yymmdd)
6 c***category no. h2a3a1,h2a4a1
7 c***keywords automatic integrator, infinite intervals,
8 c general-purpose, transformation, extrapolation,
9 c 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 integral i = integral of f over (bound,+infinity)
14 c or i = integral of f over (-infinity,bound)
15 c or i = integral of f over (-infinity,+infinity)
16 c hopefully satisfying following claim for accuracy
17 c abs(i-result).le.max(epsabs,epsrel*abs(i)).
18 c***description
19 c
20 c integration over infinite intervals
21 c standard fortran subroutine
22 c
23 c parameters
24 c on entry
25 c f - subroutine f(x,result) defining the integrand
26 c function f(x). the actual name for f needs to be
27 c declared e x t e r n a l in the driver program.
28 c
29 c bound - real
30 c finite bound of integration range
31 c (has no meaning if interval is doubly-infinite)
32 c
33 c inf - integer
34 c indicating the kind of integration range involved
35 c inf = 1 corresponds to (bound,+infinity),
36 c inf = -1 to (-infinity,bound),
37 c inf = 2 to (-infinity,+infinity).
38 c
39 c epsabs - real
40 c absolute accuracy requested
41 c epsrel - real
42 c relative accuracy requested
43 c if epsabs.le.0
44 c and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
45 c the routine will end with ier = 6.
46 c
47 c
48 c on return
49 c result - real
50 c approximation to the integral
51 c
52 c abserr - real
53 c estimate of the modulus of the absolute error,
54 c which should equal or exceed abs(i-result)
55 c
56 c neval - integer
57 c number of integrand evaluations
58 c
59 c ier - integer
60 c ier = 0 normal and reliable termination of the
61 c routine. it is assumed that the requested
62 c accuracy has been achieved.
63 c - ier.gt.0 abnormal termination of the routine. the
64 c estimates for result and error are less
65 c reliable. it is assumed that the requested
66 c accuracy has not been achieved.
67 c error messages
68 c ier = 1 maximum number of subdivisions allowed
69 c has been achieved. one can allow more
70 c subdivisions by increasing the value of
71 c limit (and taking the according dimension
72 c adjustments into account). however, if
73 c this yields no improvement it is advised
74 c to analyze the integrand in order to
75 c determine the integration difficulties. if
76 c the position of a local difficulty can be
77 c determined (e.g. singularity,
78 c discontinuity within the interval) one
79 c will probably gain from splitting up the
80 c interval at this point and calling the
81 c integrator on the subranges. if possible,
82 c an appropriate special-purpose integrator
83 c should be used, which is designed for
84 c handling the type of difficulty involved.
85 c = 2 the occurrence of roundoff error is
86 c detected, which prevents the requested
87 c tolerance from being achieved.
88 c the error may be under-estimated.
89 c = 3 extremely bad integrand behaviour occurs
90 c at some points of the integration
91 c interval.
92 c = 4 the algorithm does not converge.
93 c roundoff error is detected in the
94 c extrapolation table.
95 c it is assumed that the requested tolerance
96 c cannot be achieved, and that the returned
97 c result is the best which can be obtained.
98 c = 5 the integral is probably divergent, or
99 c slowly convergent. it must be noted that
100 c divergence can occur with any other value
101 c of ier.
102 c = 6 the input is invalid, because
103 c (epsabs.le.0 and
104 c epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
105 c or limit.lt.1 or leniw.lt.limit*4.
106 c result, abserr, neval, last are set to
107 c zero. exept when limit or leniw is
108 c invalid, iwork(1), work(limit*2+1) and
109 c work(limit*3+1) are set to zero, work(1)
110 c is set to a and work(limit+1) to b.
111 c
112 c dimensioning parameters
113 c limit - integer
114 c dimensioning parameter for iwork
115 c limit determines the maximum number of subintervals
116 c in the partition of the given integration interval
117 c (a,b), limit.ge.1.
118 c if limit.lt.1, the routine will end with ier = 6.
119 c
120 c lenw - integer
121 c dimensioning parameter for work
122 c lenw must be at least limit*4.
123 c if lenw.lt.limit*4, the routine will end
124 c with ier = 6.
125 c
126 c last - integer
127 c on return, last equals the number of subintervals
128 c produced in the subdivision process, which
129 c determines the number of significant elements
130 c actually in the work arrays.
131 c
132 c work arrays
133 c iwork - integer
134 c vector of dimension at least limit, the first
135 c k elements of which contain pointers
136 c to the error estimates over the subintervals,
137 c such that work(limit*3+iwork(1)),... ,
138 c work(limit*3+iwork(k)) form a decreasing
139 c sequence, with k = last if last.le.(limit/2+2), and
140 c k = limit+1-last otherwise
141 c
142 c work - real
143 c vector of dimension at least lenw
144 c on return
145 c work(1), ..., work(last) contain the left
146 c end points of the subintervals in the
147 c partition of (a,b),
148 c work(limit+1), ..., work(limit+last) contain
149 c the right end points,
150 c work(limit*2+1), ...,work(limit*2+last) contain the
151 c integral approximations over the subintervals,
152 c work(limit*3+1), ..., work(limit*3)
153 c contain the error estimates.
154 c***references (none)
155 c***routines called qagie,xerror
156 c***end prologue qagi
157 c
158  real abserr, epsabs,epsrel,result,work
159  integer ier,iwork, lenw,limit,lvl,l1,l2,l3,neval
160 c
161  dimension iwork(limit),work(lenw)
162 c
163  external f
164 c
165 c check validity of limit and lenw.
166 c
167 c***first executable statement qagi
168  ier = 6
169  neval = 0
170  last = 0
171  result = 0.0e+00
172  abserr = 0.0e+00
173  if(limit.lt.1.or.lenw.lt.limit*4) go to 10
174 c
175 c prepare call for qagie.
176 c
177  l1 = limit+1
178  l2 = limit+l1
179  l3 = limit+l2
180 c
181  call qagie(f,bound,inf,epsabs,epsrel,limit,result,abserr,
182  * neval,ier,work(1),work(l1),work(l2),work(l3),iwork,last)
183 c
184 c call error handler if necessary.
185 c
186  lvl = 0
187 10 if(ier.eq.6) lvl = 1
188  if(ier.ne.0) call xerror('abnormal return from qagi',26,ier,lvl)
189  return
190  end
Generated on Mon Dec 30 2013 03:04:47 for GNU Octave by   doxygen 1.8.1.2