1
subroutine dqagi(f,bound,inf,epsabs,epsrel,result,abserr,neval,
2
* ier,limit,lenw,last,iwork,work)
3
c***begin prologue dqagi
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,
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)).
20
c integration over infinite intervals
21
c standard fortran subroutine
25
c f - double precision
26
c function subprogram 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.
30
c bound - double precision
31
c finite bound of integration range
32
c (has no meaning if interval is doubly-infinite)
35
c indicating the kind of integration range involved
36
c inf = 1 corresponds to (bound,+infinity),
37
c inf = -1 to (-infinity,bound),
38
c inf = 2 to (-infinity,+infinity).
40
c epsabs - double precision
41
c absolute accuracy requested
42
c epsrel - double precision
43
c relative accuracy requested
45
c and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
46
c the routine will end with ier = 6.
50
c result - double precision
51
c approximation to the integral
53
c abserr - double precision
54
c estimate of the modulus of the absolute error,
55
c which should equal or exceed abs(i-result)
58
c number of integrand evaluations
61
c ier = 0 normal and reliable termination of the
62
c routine. it is assumed that the requested
63
c accuracy has been achieved.
64
c - ier.gt.0 abnormal termination of the routine. the
65
c estimates for result and error are less
66
c reliable. it is assumed that the requested
67
c accuracy has not been achieved.
69
c ier = 1 maximum number of subdivisions allowed
70
c has been achieved. one can allow more
71
c subdivisions by increasing the value of
72
c limit (and taking the according dimension
73
c adjustments into account). however, if
74
c this yields no improvement it is advised
75
c to analyze the integrand in order to
76
c determine the integration difficulties. if
77
c the position of a local difficulty can be
78
c determined (e.g. singularity,
79
c discontinuity within the interval) one
80
c will probably gain from splitting up the
81
c interval at this point and calling the
82
c integrator on the subranges. if possible,
83
c an appropriate special-purpose integrator
84
c should be used, which is designed for
85
c handling the type of difficulty involved.
86
c = 2 the occurrence of roundoff error is
87
c detected, which prevents the requested
88
c tolerance from being achieved.
89
c the error may be under-estimated.
90
c = 3 extremely bad integrand behaviour occurs
91
c at some points of the integration
93
c = 4 the algorithm does not converge.
94
c roundoff error is detected in the
95
c extrapolation table.
96
c it is assumed that the requested tolerance
97
c cannot be achieved, and that the returned
98
c result is the best which can be obtained.
99
c = 5 the integral is probably divergent, or
100
c slowly convergent. it must be noted that
101
c divergence can occur with any other value
103
c = 6 the input is invalid, because
105
c epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
106
c or limit.lt.1 or leniw.lt.limit*4.
107
c result, abserr, neval, last are set to
108
c zero. exept when limit or leniw is
109
c invalid, iwork(1), work(limit*2+1) and
110
c work(limit*3+1) are set to zero, work(1)
111
c is set to a and work(limit+1) to b.
113
c dimensioning parameters
115
c dimensioning parameter for iwork
116
c limit determines the maximum number of subintervals
117
c in the partition of the given integration interval
119
c if limit.lt.1, the routine will end with ier = 6.
122
c dimensioning parameter for work
123
c lenw must be at least limit*4.
124
c if lenw.lt.limit*4, the routine will end
128
c on return, last equals the number of subintervals
129
c produced in the subdivision process, which
130
c determines the number of significant elements
131
c actually in the work arrays.
135
c vector of dimension at least limit, the first
136
c k elements of which contain pointers
137
c to the error estimates over the subintervals,
138
c such that work(limit*3+iwork(1)),... ,
139
c work(limit*3+iwork(k)) form a decreasing
140
c sequence, with k = last if last.le.(limit/2+2), and
141
c k = limit+1-last otherwise
143
c work - double precision
144
c vector of dimension at least lenw
146
c work(1), ..., work(last) contain the left
147
c end points of the subintervals in the
148
c partition of (a,b),
149
c work(limit+1), ..., work(limit+last) contain
150
c the right end points,
151
c work(limit*2+1), ...,work(limit*2+last) contain the
152
c integral approximations over the subintervals,
153
c work(limit*3+1), ..., work(limit*3+last)
154
c contain the error estimates.
155
c***references (none)
156
c***routines called dqagie,xerror
157
c***end prologue dqagi
159
double precision abserr,bound,epsabs,epsrel,f,result,work
160
integer ier,inf,iwork,last,lenw,limit,lvl,l1,l2,l3,neval
162
dimension iwork(limit),work(lenw)
166
c check validity of limit and lenw.
168
c***first executable statement dqagi
174
if(limit.lt.1.or.lenw.lt.limit*4) go to 10
176
c prepare call for dqagie.
182
call dqagie(f,bound,inf,epsabs,epsrel,limit,result,abserr,
183
* neval,ier,work(1),work(l1),work(l2),work(l3),iwork,last)
185
c call error handler if necessary.
188
10 if(ier.eq.6) lvl = 1
189
if(ier.ne.0) call xerror(26habnormal return from dqagi,26,ier,lvl)
added
removed
Lib/integrate/quadpack/dqagi.f
