1
subroutine dqag(f,a,b,epsabs,epsrel,key,result,abserr,neval,ier,
2
* limit,lenw,last,iwork,work)
3
c***begin prologue dqag
4
c***date written 800101 (yymmdd)
5
c***revision date 830518 (yymmdd)
6
c***category no. h2a1a1
7
c***keywords automatic integrator, general-purpose,
8
c integrand examinator, 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 abs(i-result)le.max(epsabs,epsrel*abs(i)).
18
c computation of a definite integral
19
c standard fortran subroutine
20
c double precision version
22
c f - double precision
23
c function subprogam defining the integrand
24
c function f(x). the actual name for f needs to be
25
c declared e x t e r n a l in the driver program.
27
c a - double precision
28
c lower limit of integration
30
c b - double precision
31
c upper limit of integration
33
c epsabs - double precision
34
c absolute accoracy requested
35
c epsrel - double precision
36
c relative accuracy requested
38
c and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
39
c the routine will end with ier = 6.
42
c key for choice of local integration rule
43
c a gauss-kronrod pair is used with
44
c 7 - 15 points if key.lt.2,
45
c 10 - 21 points if key = 2,
46
c 15 - 31 points if key = 3,
47
c 20 - 41 points if key = 4,
48
c 25 - 51 points if key = 5,
49
c 30 - 61 points if key.gt.5.
52
c result - double precision
53
c approximation to the integral
55
c abserr - double precision
56
c estimate of the modulus of the absolute error,
57
c which should equal or exceed abs(i-result)
60
c number of integrand evaluations
63
c ier = 0 normal and reliable termination of the
64
c routine. it is assumed that the requested
65
c accuracy has been achieved.
66
c ier.gt.0 abnormal termination of the routine
67
c the estimates for result and error are
68
c less reliable. it is assumed that the
69
c requested accuracy has not been achieved.
71
c ier = 1 maximum number of subdivisions allowed
72
c has been achieved. one can allow more
73
c subdivisions by increasing the value of
74
c limit (and taking the according dimension
75
c adjustments into account). however, if
76
c this yield no improvement it is advised
77
c to analyze the integrand in order to
78
c determine the integration difficulaties.
79
c if the position of a local difficulty can
80
c be determined (i.e.singularity,
81
c discontinuity within the interval) one
82
c will probably gain from splitting up the
83
c interval at this point and calling the
84
c integrator on the subranges. if possible,
85
c an appropriate special-purpose integrator
86
c should be used which is designed for
87
c handling the type of difficulty involved.
88
c = 2 the occurrence of roundoff error is
89
c detected, which prevents the requested
90
c tolerance from being achieved.
91
c = 3 extremely bad integrand behaviour occurs
92
c at some points of the integration
94
c = 6 the input is invalid, because
96
c epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
97
c or limit.lt.1 or lenw.lt.limit*4.
98
c result, abserr, neval, last are set
100
c except when lenw is invalid, iwork(1),
101
c work(limit*2+1) and work(limit*3+1) are
102
c set to zero, work(1) is set to a and
103
c work(limit+1) to b.
105
c dimensioning parameters
107
c dimensioning parameter for iwork
108
c limit determines the maximum number of subintervals
109
c in the partition of the given integration interval
111
c if limit.lt.1, the routine will end with ier = 6.
114
c dimensioning parameter for work
115
c lenw must be at least limit*4.
116
c if lenw.lt.limit*4, the routine will end with
120
c on return, last equals the number of subintervals
121
c produced in the subdiviosion process, which
122
c determines the number of significant elements
123
c actually in the work arrays.
127
c vector of dimension at least limit, the first k
128
c elements of which contain pointers to the error
129
c estimates over the subintervals, such that
130
c work(limit*3+iwork(1)),... , work(limit*3+iwork(k))
131
c form a decreasing sequence with k = last if
132
c last.le.(limit/2+2), and k = limit+1-last otherwise
134
c work - double precision
135
c vector of dimension at least lenw
137
c work(1), ..., work(last) contain the left end
138
c points of the subintervals in the partition of
140
c work(limit+1), ..., work(limit+last) contain the
142
c work(limit*2+1), ..., work(limit*2+last) contain
143
c the integral approximations over the subintervals,
144
c work(limit*3+1), ..., work(limit*3+last) contain
145
c the error estimates.
147
c***references (none)
148
c***routines called dqage,xerror
149
c***end prologue dqag
150
double precision a,abserr,b,epsabs,epsrel,f,result,work
151
integer ier,iwork,key,last,lenw,limit,lvl,l1,l2,l3,neval
153
dimension iwork(limit),work(lenw)
157
c check validity of lenw.
159
c***first executable statement dqag
165
if(limit.lt.1.or.lenw.lt.limit*4) go to 10
167
c prepare call for dqage.
173
call dqage(f,a,b,epsabs,epsrel,key,limit,result,abserr,neval,
174
* ier,work(1),work(l1),work(l2),work(l3),iwork,last)
176
c call error handler if necessary.
179
10 if(ier.eq.6) lvl = 1
180
if(ier.ne.0) call xerror('abnormal return from dqag' ,26,ier,lvl)