1       subroutine dqagse(f,a,b,epsabs,epsrel,limit,result,abserr,neval,
2      *   ier,alist,blist,rlist,elist,iord,last)
3 c***begin prologue  dqagse
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             (end point) singularities, 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            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)).
16 c***description
17 c
18 c        computation of a definite integral
19 c        standard fortran subroutine
20 c        double precision version
21 c
22 c        parameters
23 c         on entry
24 c            f      - double precision
25 c                     function subprogram 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            a      - double precision
30 c                     lower limit of integration
31 c
32 c            b      - double precision
33 c                     upper limit of integration
34 c
35 c            epsabs - double precision
36 c                     absolute accuracy requested
37 c            epsrel - double precision
38 c                     relative accuracy requested
39 c                     if  epsabs.le.0
40 c                     and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
41 c                     the routine will end with ier = 6.
42 c
43 c            limit  - integer
44 c                     gives an upperbound on the number of subintervals
45 c                     in the partition of (a,b)
46 c
47 c         on return
48 c            result - double precision
49 c                     approximation to the integral
50 c
51 c            abserr - double precision
52 c                     estimate of the modulus of the absolute error,
53 c                     which should equal or exceed abs(i-result)
54 c
55 c            neval  - integer
56 c                     number of integrand evaluations
57 c
58 c            ier    - integer
59 c                     ier = 0 normal and reliable termination of the
60 c                             routine. it is assumed that the requested
61 c                             accuracy has been achieved.
62 c                     ier.gt.0 abnormal termination of the routine
63 c                             the estimates for integral and error are
64 c                             less reliable. it is assumed that the
65 c                             requested accuracy has not been achieved.
66 c            error messages
67 c                         = 1 maximum number of subdivisions allowed
68 c                             has been achieved. one can allow more sub-
69 c                             divisions by increasing the value of limit
70 c                             (and taking the according dimension
71 c                             adjustments into account). however, if
72 c                             this yields no improvement it is advised
73 c                             to analyze the integrand in order to
74 c                             determine the integration difficulties. if
75 c                             the position of a local difficulty can be
76 c                             determined (e.g. singularity,
77 c                             discontinuity within the interval) one
78 c                             will probably gain from splitting up the
79 c                             interval at this point and calling the
80 c                             integrator on the subranges. if possible,
81 c                             an appropriate special-purpose integrator
82 c                             should be used, which is designed for
83 c                             handling the type of difficulty involved.
84 c                         = 2 the occurrence of roundoff error is detec-
85 c                             ted, which prevents the requested
86 c                             tolerance from being achieved.
87 c                             the error may be under-estimated.
88 c                         = 3 extremely bad integrand behaviour
89 c                             occurs at some points of the integration
90 c                             interval.
91 c                         = 4 the algorithm does not converge.
92 c                             roundoff error is detected in the
93 c                             extrapolation table.
94 c                             it is presumed that the requested
95 c                             tolerance cannot be achieved, and that the
96 c                             returned result is the best which can be
97 c                             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                             result, abserr, neval, last, rlist(1),
106 c                             iord(1) and elist(1) are set to zero.
107 c                             alist(1) and blist(1) are set to a and b
108 c                             respectively.
109 c
110 c            alist  - double precision
111 c                     vector of dimension at least limit, the first
112 c                      last  elements of which are the left end points
113 c                     of the subintervals in the partition of the
114 c                     given integration range (a,b)
115 c
116 c            blist  - double precision
117 c                     vector of dimension at least limit, the first
118 c                      last  elements of which are the right end points
119 c                     of the subintervals in the partition of the given
120 c                     integration range (a,b)
121 c
122 c            rlist  - double precision
123 c                     vector of dimension at least limit, the first
124 c                      last  elements of which are the integral
125 c                     approximations on the subintervals
126 c
127 c            elist  - double precision
128 c                     vector of dimension at least limit, the first
129 c                      last  elements of which are the moduli of the
130 c                     absolute error estimates on the subintervals
131 c
132 c            iord   - integer
133 c                     vector of dimension at least limit, the first k
134 c                     elements of which are pointers to the
135 c                     error estimates over the subintervals,
136 c                     such that elist(iord(1)), ..., elist(iord(k))
137 c                     form a decreasing sequence, with k = last
138 c                     if last.le.(limit/2+2), and k = limit+1-last
139 c                     otherwise
140 c
141 c            last   - integer
142 c                     number of subintervals actually produced in the
143 c                     subdivision process
144 c
145 c***references  (none)
146 c***routines called  d1mach,dqelg,dqk21,dqpsrt
147 c***Scilab Enterprises input
148 c        2013 - Paul Bignier
149 c        internal parameter nbisect is the number of bisections carried out.
150 c        if this parameter exceeds two, then a warning is issued to
151 c        suggest the user to reduce the integration interval if he thinks
152 c        edges have been missed.
153 c***end prologue  dqagse
154 c
155       double precision a,abseps,abserr,alist,area,area1,area12,area2,a1,
156      *  a2,b,blist,b1,b2,correc,dabs,defabs,defab1,defab2,d1mach,dmax1,
157      *  dres,elist,epmach,epsabs,epsrel,erlarg,erlast,errbnd,errmax,
158      *  error1,error2,erro12,errsum,ertest,f,oflow,resabs,reseps,result,
159      *  res3la,rlist,rlist2,small,uflow
160       integer id,ier,ierro,iord,iroff1,iroff2,iroff3,jupbnd,k,ksgn,
161      *  ktmin,last,limit,maxerr,neval,nres,nrmax,numrl2
162 c     /* Scilab Enterprises input
163       integer nbisect
164 c     */
165       logical extrap,noext
166 c
167       dimension alist(limit),blist(limit),elist(limit),iord(limit),
168      * res3la(3),rlist(limit),rlist2(52)
169 c
170       external f
171 c
172 c            the dimension of rlist2 is determined by the value of
173 c            limexp in subroutine dqelg (rlist2 should be of dimension
174 c            (limexp+2) at least).
175 c
176 c            list of major variables
177 c            -----------------------
178 c
179 c           alist     - list of left end points of all subintervals
180 c                       considered up to now
181 c           blist     - list of right end points of all subintervals
182 c                       considered up to now
183 c           rlist(i)  - approximation to the integral over
184 c                       (alist(i),blist(i))
185 c           rlist2    - array of dimension at least limexp+2 containing
186 c                       the part of the epsilon table which is still
187 c                       needed for further computations
188 c           elist(i)  - error estimate applying to rlist(i)
189 c           maxerr    - pointer to the interval with largest error
190 c                       estimate
191 c           errmax    - elist(maxerr)
192 c           erlast    - error on the interval currently subdivided
193 c                       (before that subdivision has taken place)
194 c           area      - sum of the integrals over the subintervals
195 c           errsum    - sum of the errors over the subintervals
196 c           errbnd    - requested accuracy max(epsabs,epsrel*
197 c                       abs(result))
198 c           *****1    - variable for the left interval
199 c           *****2    - variable for the right interval
200 c           last      - index for subdivision
201 c           nres      - number of calls to the extrapolation routine
202 c           numrl2    - number of elements currently in rlist2. if an
203 c                       appropriate approximation to the compounded
204 c                       integral has been obtained it is put in
205 c                       rlist2(numrl2) after numrl2 has been increased
206 c                       by one.
207 c           small     - length of the smallest interval considered up
208 c                       to now, multiplied by 1.5
209 c           erlarg    - sum of the errors over the intervals larger
210 c                       than the smallest interval considered up to now
211 c           extrap    - logical variable denoting that the routine is
212 c                       attempting to perform extrapolation i.e. before
213 c                       subdividing the smallest interval we try to
214 c                       decrease the value of erlarg.
215 c           noext     - logical variable denoting that extrapolation
216 c                       is no longer allowed (true value)
217 c
218 c            machine dependent constants
219 c            ---------------------------
220 c
221 c           epmach is the largest relative spacing.
222 c           uflow is the smallest positive magnitude.
223 c           oflow is the largest positive magnitude.
224 c
225 c***first executable statement  dqagse
226       epmach = d1mach(4)
227 c
228 c            test on validity of parameters
229 c            ------------------------------
230       ier = 0
231       neval = 0
232       last = 0
233       result = 0.0d+00
234       abserr = 0.0d+00
235       alist(1) = a
236       blist(1) = b
237       rlist(1) = 0.0d+00
238       elist(1) = 0.0d+00
239       if(epsabs.le.0.0d+00.and.epsrel.lt.dmax1(0.5d+02*epmach,0.5d-28))
240      *   ier = 6
241       if(ier.eq.6) go to 999
242 c
243 c           first approximation to the integral
244 c           -----------------------------------
245 c
246       uflow = d1mach(1)
247       oflow = d1mach(2)
248       ierro = 0
249 c     /* Scilab Enterprises input
250       nbisect = 0
251 c     */
252       call dqk21(f,a,b,result,abserr,defabs,resabs)
253 c
254 c           test on accuracy.
255 c
256       dres = dabs(result)
257       errbnd = dmax1(epsabs,epsrel*dres)
258       last = 1
259       rlist(1) = result
260       elist(1) = abserr
261       iord(1) = 1
262       if(abserr.le.1.0d+02*epmach*defabs.and.abserr.gt.errbnd) ier = 2
263       if(limit.eq.1) ier = 1
264       if(ier.ne.0.or.(abserr.le.errbnd.and.abserr.ne.resabs).or.
265      *  abserr.eq.0.0d+00) go to 140
266 c
267 c           initialization
268 c           --------------
269 c
270       rlist2(1) = result
271       errmax = abserr
272       maxerr = 1
273       area = result
274       errsum = abserr
275       abserr = oflow
276       nrmax = 1
277       nres = 0
278       numrl2 = 2
279       ktmin = 0
280       extrap = .false.
281       noext = .false.
282       iroff1 = 0
283       iroff2 = 0
284       iroff3 = 0
285       ksgn = -1
286       if(dres.ge.(0.1d+01-0.5d+02*epmach)*defabs) ksgn = 1
287 c
288 c           main do-loop
289 c           ------------
290 c
291       do 90 last = 2,limit
292 c
293 c           bisect the subinterval with the nrmax-th largest error
294 c           estimate.
295 c
296 c       /* Scilab Enterprises input
297         nbisect = nbisect + 1
298 c       */
299         a1 = alist(maxerr)
300         b1 = 0.5d+00*(alist(maxerr)+blist(maxerr))
301         a2 = b1
302         b2 = blist(maxerr)
303         erlast = errmax
304         call dqk21(f,a1,b1,area1,error1,resabs,defab1)
305         call dqk21(f,a2,b2,area2,error2,resabs,defab2)
306 c
307 c           improve previous approximations to integral
308 c           and error and test for accuracy.
309 c
310         area12 = area1+area2
311         erro12 = error1+error2
312         errsum = errsum+erro12-errmax
313         area = area+area12-rlist(maxerr)
314         if(defab1.eq.error1.or.defab2.eq.error2) go to 15
315         if(dabs(rlist(maxerr)-area12).gt.0.1d-04*dabs(area12)
316      *  .or.erro12.lt.0.99d+00*errmax) go to 10
317         if(extrap) iroff2 = iroff2+1
318         if(.not.extrap) iroff1 = iroff1+1
319    10   if(last.gt.10.and.erro12.gt.errmax) iroff3 = iroff3+1
320    15   rlist(maxerr) = area1
321         rlist(last) = area2
322         errbnd = dmax1(epsabs,epsrel*dabs(area))
323 c
324 c           test for roundoff error and eventually set error flag.
325 c
326         if(iroff1+iroff2.ge.10.or.iroff3.ge.20) ier = 2
327         if(iroff2.ge.5) ierro = 3
328 c
329 c           set error flag in the case that the number of subintervals
330 c           equals limit.
331 c
332         if(last.eq.limit) ier = 1
333 c
334 c           set error flag in the case of bad integrand behaviour
335 c           at a point of the integration range.
336 c
337         if(dmax1(dabs(a1),dabs(b2)).le.(0.1d+01+0.1d+03*epmach)*
338      *  (dabs(a2)+0.1d+04*uflow)) ier = 4
339 c
340 c           append the newly-created intervals to the list.
341 c
342         if(error2.gt.error1) go to 20
343         alist(last) = a2
344         blist(maxerr) = b1
345         blist(last) = b2
346         elist(maxerr) = error1
347         elist(last) = error2
348         go to 30
349    20   alist(maxerr) = a2
350         alist(last) = a1
351         blist(last) = b1
352         rlist(maxerr) = area2
353         rlist(last) = area1
354         elist(maxerr) = error2
355         elist(last) = error1
356 c
357 c           call subroutine dqpsrt to maintain the descending ordering
358 c           in the list of error estimates and select the subinterval
359 c           with nrmax-th largest error estimate (to be bisected next).
360 c
361    30   call dqpsrt(limit,last,maxerr,errmax,elist,iord,nrmax)
362 c ***jump out of do-loop
363         if(errsum.le.errbnd) go to 115
364 c ***jump out of do-loop
365         if(ier.ne.0) go to 100
366         if(last.eq.2) go to 80
367         if(noext) go to 90
368         erlarg = erlarg-erlast
369         if(dabs(b1-a1).gt.small) erlarg = erlarg+erro12
370         if(extrap) go to 40
371 c
372 c           test whether the interval to be bisected next is the
373 c           smallest interval.
374 c
375         if(dabs(blist(maxerr)-alist(maxerr)).gt.small) go to 90
376         extrap = .true.
377         nrmax = 2
378    40   if(ierro.eq.3.or.erlarg.le.ertest) go to 60
379 c
380 c           the smallest interval has the largest error.
381 c           before bisecting decrease the sum of the errors over the
382 c           larger intervals (erlarg) and perform extrapolation.
383 c
384         id = nrmax
385         jupbnd = last
386         if(last.gt.(2+limit/2)) jupbnd = limit+3-last
387         do 50 k = id,jupbnd
388           maxerr = iord(nrmax)
389           errmax = elist(maxerr)
390 c ***jump out of do-loop
391           if(dabs(blist(maxerr)-alist(maxerr)).gt.small) go to 90
392           nrmax = nrmax+1
393    50   continue
394 c
395 c           perform extrapolation.
396 c
397    60   numrl2 = numrl2+1
398         rlist2(numrl2) = area
399         call dqelg(numrl2,rlist2,reseps,abseps,res3la,nres)
400         ktmin = ktmin+1
401         if(ktmin.gt.5.and.abserr.lt.0.1d-02*errsum) ier = 5
402         if(abseps.ge.abserr) go to 70
403         ktmin = 0
404         abserr = abseps
405         result = reseps
406         correc = erlarg
407         ertest = dmax1(epsabs,epsrel*dabs(reseps))
408 c ***jump out of do-loop
409         if(abserr.le.ertest) go to 100
410 c
411 c           prepare bisection of the smallest interval.
412 c
413    70   if(numrl2.eq.1) noext = .true.
414         if(ier.eq.5) go to 100
415         maxerr = iord(1)
416         errmax = elist(maxerr)
417         nrmax = 1
418         extrap = .false.
419         small = small*0.5d+00
420         erlarg = errsum
421         go to 90
422    80   small = dabs(b-a)*0.375d+00
423         erlarg = errsum
424         ertest = errbnd
425         rlist2(2) = area
426    90 continue
427 c
428 c           set final result and error estimate.
429 c           ------------------------------------
430 c
431   100 if(abserr.eq.oflow) go to 115
432       if(ier+ierro.eq.0) go to 110
433       if(ierro.eq.3) abserr = abserr+correc
434       if(ier.eq.0) ier = 3
435       if(result.ne.0.0d+00.and.area.ne.0.0d+00) go to 105
436       if(abserr.gt.errsum) go to 115
437       if(area.eq.0.0d+00) go to 130
438       go to 110
439   105 if(abserr/dabs(result).gt.errsum/dabs(area)) go to 115
440 c
441 c           test on divergence.
442 c
443   110 if(ksgn.eq.(-1).and.dmax1(dabs(result),dabs(area)).le.
444      * defabs*0.1d-01) go to 130
445       if(0.1d-01.gt.(result/area).or.(result/area).gt.0.1d+03
446      * .or.errsum.gt.dabs(area)) ier = 6
447       go to 130
448 c
449 c           compute global integral sum.
450 c
451   115 result = 0.0d+00
452       do 120 k = 1,last
453          result = result+rlist(k)
454   120 continue
455       abserr = errsum
456   130 if(ier.gt.2) ier = ier-1
457   140 neval = 42*last-21
458 c     /* Scilab Enterprises input
459       if (nbisect.le.1) then
460          call msgstxt('Warning: argument function is detected to be '//
461      &      'very smooth.')
462          call msgstxt('Reduce integration interval if you think '//
463      &      'edges have been missed.')
464       endif
465 c     */
466   999 return
467       end