1      SUBROUTINE DQAGPE(F,A,B,NPTS2,POINTS,EPSABS,EPSREL,LIMIT,RESULT,
2     1   ABSERR,NEVAL,IER,ALIST,BLIST,RLIST,ELIST,PTS,IORD,LEVEL,NDIN,
3     2   LAST)
4C***BEGIN PROLOGUE  DQAGPE
5C***DATE WRITTEN   800101   (YYMMDD)
6C***REVISION DATE  840425   (YYMMDD)
7C***REVISION HISTORY (YYMMDD)
8C   000601   Changed DMAX1/DMIN1/DABS to generic MAX/MIN/ABS
9C***CATEGORY NO.  H2A2A1
10C***KEYWORDS  AUTOMATIC INTEGRATOR,EXTRAPOLATION,GENERAL-PURPOSE,
11C             GLOBALLY ADAPTIVE.,SINGULARITIES AT USER SPECIFIED POINTS
12C***AUTHOR  PIESSENS, ROBERT, APPLIED MATH. AND PROGR. DIV. -
13C             K. U. LEUVEN
14C           DE DONCKER, ELISE, APPLIED MATH. AND PROGR. DIV. -
15C             K. U. LEUVEN
16C***PURPOSE  The routine calculates an approximation result to a given
17C            definite integral I = Integral of F over (A,B), hopefully
18C            satisfying following claim for accuracy ABS(I-RESULT).LE.
19C            MAX(EPSABS,EPSREL*ABS(I)). Break points of the integration
20C            interval, where local difficulties of the integrand may
21C            occur(e.g. singularities,discontinuities),provided by user.
22C***DESCRIPTION
23C
24C        Computation of a definite integral
25C        Standard fortran subroutine
26C        Double precision version
27C
28C        PARAMETERS
29C         ON ENTRY
30C            F      - Double precision
31C                     Function subprogram defining the integrand
32C                     function F(X). The actual name for F needs to be
33C                     declared E X T E R N A L in the driver program.
34C
35C            A      - Double precision
36C                     Lower limit of integration
37C
38C            B      - Double precision
39C                     Upper limit of integration
40C
41C            NPTS2  - Integer
42C                     Number equal to two more than the number of
43C                     user-supplied break points within the integration
44C                     range, NPTS2.GE.2.
45C                     If NPTS2.LT.2, the routine will end with IER = 6.
46C
47C            POINTS - Double precision
48C                     Vector of dimension NPTS2, the first (NPTS2-2)
49C                     elements of which are the user provided break
50C                     POINTS. If these POINTS do not constitute an
51C                     ascending sequence there will be an automatic
52C                     sorting.
53C
54C            EPSABS - Double precision
55C                     Absolute accuracy requested
56C            EPSREL - Double precision
57C                     Relative accuracy requested
58C                     If  EPSABS.LE.0
59C                     and EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28),
60C                     the routine will end with IER = 6.
61C
62C            LIMIT  - Integer
63C                     Gives an upper bound on the number of subintervals
64C                     in the partition of (A,B), LIMIT.GE.NPTS2
65C                     If LIMIT.LT.NPTS2, the routine will end with
66C                     IER = 6.
67C
68C         ON RETURN
69C            RESULT - Double precision
70C                     Approximation to the integral
71C
72C            ABSERR - Double precision
73C                     Estimate of the modulus of the absolute error,
74C                     which should equal or exceed ABS(I-RESULT)
75C
76C            NEVAL  - Integer
77C                     Number of integrand evaluations
78C
79C            IER    - Integer
80C                     IER = 0 Normal and reliable termination of the
81C                             routine. It is assumed that the requested
82C                             accuracy has been achieved.
83C                     IER.GT.0 Abnormal termination of the routine.
84C                             The estimates for integral and error are
85C                             less reliable. It is assumed that the
86C                             requested accuracy has not been achieved.
87C            ERROR MESSAGES
88C                     IER = 1 Maximum number of subdivisions allowed
89C                             has been achieved. One can allow more
90C                             subdivisions by increasing the value of
91C                             LIMIT (and taking the according dimension
92C                             adjustments into account). However, if
93C                             this yields no improvement it is advised
94C                             to analyze the integrand in order to
95C                             determine the integration difficulties. If
96C                             the position of a local difficulty can be
97C                             determined (i.e. SINGULARITY,
98C                             DISCONTINUITY within the interval), it
99C                             should be supplied to the routine as an
100C                             element of the vector points. If necessary
101C                             an appropriate special-purpose integrator
102C                             must be used, which is designed for
103C                             handling the type of difficulty involved.
104C                         = 2 The occurrence of roundoff error is
105C                             detected, which prevents the requested
106C                             tolerance from being achieved.
107C                             The error may be under-estimated.
108C                         = 3 Extremely bad integrand behaviour occurs
109C                             At some points of the integration
110C                             interval.
111C                         = 4 The algorithm does not converge.
112C                             Roundoff error is detected in the
113C                             extrapolation table. It is presumed that
114C                             the requested tolerance cannot be
115C                             achieved, and that the returned result is
116C                             the best which can be obtained.
117C                         = 5 The integral is probably divergent, or
118C                             slowly convergent. It must be noted that
119C                             divergence can occur with any other value
120C                             of IER.GT.0.
121C                         = 6 The input is invalid because
122C                             NPTS2.LT.2 or
123C                             Break points are specified outside
124C                             the integration range or
125C                             (EPSABS.LE.0 and
126C                              EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28))
127C                             or LIMIT.LT.NPTS2.
128C                             RESULT, ABSERR, NEVAL, LAST, RLIST(1),
129C                             and ELIST(1) are set to zero. ALIST(1) and
130C                             BLIST(1) are set to A and B respectively.
131C
132C            ALIST  - Double precision
133C                     Vector of dimension at least LIMIT, the first
134C                      LAST  elements of which are the left end points
135C                     of the subintervals in the partition of the given
136C                     integration range (A,B)
137C
138C            BLIST  - Double precision
139C                     Vector of dimension at least LIMIT, the first
140C                      LAST  elements of which are the right end points
141C                     of the subintervals in the partition of the given
142C                     integration range (A,B)
143C
144C            RLIST  - Double precision
145C                     Vector of dimension at least LIMIT, the first
146C                      LAST  elements of which are the integral
147C                     approximations on the subintervals
148C
149C            ELIST  - Double precision
150C                     Vector of dimension at least LIMIT, the first
151C                      LAST  elements of which are the moduli of the
152C                     absolute error estimates on the subintervals
153C
154C            PTS    - Double precision
155C                     Vector of dimension at least NPTS2, containing the
156C                     integration limits and the break points of the
157C                     interval in ascending sequence.
158C
159C            LEVEL  - Integer
160C                     Vector of dimension at least LIMIT, containing the
161C                     subdivision levels of the subinterval, i.e. if
162C                     (AA,BB) is a subinterval of (P1,P2) where P1 as
163C                     well as P2 is a user-provided break point or
164C                     integration limit, then (AA,BB) has level L if
165C                     ABS(BB-AA) = ABS(P2-P1)*2**(-L).
166C
167C            NDIN   - Integer
168C                     Vector of dimension at least NPTS2, after first
169C                     integration over the intervals (PTS(I)),PTS(I+1),
170C                     I = 0,1, ..., NPTS2-2, the error estimates over
171C                     some of the intervals may have been increased
172C                     artificially, in order to put their subdivision
173C                     forward. If this happens for the subinterval
174C                     numbered K, NDIN(K) is put to 1, otherwise
175C                     NDIN(K) = 0.
176C
177C            IORD   - Integer
178C                     Vector of dimension at least LIMIT, the first K
179C                     elements of which are pointers to the
180C                     error estimates over the subintervals,
181C                     such that ELIST(IORD(1)), ..., ELIST(IORD(K))
182C                     form a decreasing sequence, with K = LAST
183C                     If LAST.LE.(LIMIT/2+2), and K = LIMIT+1-LAST
184C                     otherwise
185C
186C            LAST   - Integer
187C                     Number of subintervals actually produced in the
188C                     subdivisions process
189C***REFERENCES  (NONE)
190C***ROUTINES CALLED  D1MACH,DQELG,DQK21,DQPSRT
191C***END PROLOGUE  DQAGPE
192      DOUBLE PRECISION A,ABSEPS,ABSERR,ALIST,AREA,AREA1,AREA12,AREA2,A1,
193     1  A2,B,BLIST,B1,B2,CORREC,DEFABS,DEFAB1,DEFAB2,
194     2  DRES,D1MACH,ELIST,EPMACH,EPSABS,EPSREL,ERLARG,ERLAST,ERRBND,
195     3  ERRMAX,ERROR1,ERRO12,ERROR2,ERRSUM,ERTEST,F,OFLOW,POINTS,PTS,
196     4  RESA,RESABS,RESEPS,RESULT,RES3LA,RLIST,RLIST2,SIGN,TEMP,UFLOW
197      INTEGER I,ID,IER,IERRO,IND1,IND2,IORD,IP1,IROFF1,IROFF2,IROFF3,J,
198     1  JLOW,JUPBND,K,KSGN,KTMIN,LAST,LEVCUR,LEVEL,LEVMAX,LIMIT,MAXERR,
199     2  NDIN,NEVAL,NINT,NINTP1,NPTS,NPTS2,NRES,NRMAX,NUMRL2
200      LOGICAL EXTRAP,NOEXT
201C
202C
203      DIMENSION ALIST(LIMIT),BLIST(LIMIT),ELIST(LIMIT),IORD(LIMIT),
204     1  LEVEL(LIMIT),NDIN(NPTS2),POINTS(NPTS2),PTS(NPTS2),RES3LA(3),
205     2  RLIST(LIMIT),RLIST2(52)
206C
207      EXTERNAL F
208C
209C            THE DIMENSION OF RLIST2 IS DETERMINED BY THE VALUE OF
210C            LIMEXP IN SUBROUTINE EPSALG (RLIST2 SHOULD BE OF DIMENSION
211C            (LIMEXP+2) AT LEAST).
212C
213C
214C            LIST OF MAJOR VARIABLES
215C            -----------------------
216C
217C           ALIST     - LIST OF LEFT END POINTS OF ALL SUBINTERVALS
218C                       CONSIDERED UP TO NOW
219C           BLIST     - LIST OF RIGHT END POINTS OF ALL SUBINTERVALS
220C                       CONSIDERED UP TO NOW
221C           RLIST(I)  - APPROXIMATION TO THE INTEGRAL OVER
222C                       (ALIST(I),BLIST(I))
223C           RLIST2    - ARRAY OF DIMENSION AT LEAST LIMEXP+2
224C                       CONTAINING THE PART OF THE EPSILON TABLE WHICH
225C                       IS STILL NEEDED FOR FURTHER COMPUTATIONS
226C           ELIST(I)  - ERROR ESTIMATE APPLYING TO RLIST(I)
227C           MAXERR    - POINTER TO THE INTERVAL WITH LARGEST ERROR
228C                       ESTIMATE
229C           ERRMAX    - ELIST(MAXERR)
230C           ERLAST    - ERROR ON THE INTERVAL CURRENTLY SUBDIVIDED
231C                       (BEFORE THAT SUBDIVISION HAS TAKEN PLACE)
232C           AREA      - SUM OF THE INTEGRALS OVER THE SUBINTERVALS
233C           ERRSUM    - SUM OF THE ERRORS OVER THE SUBINTERVALS
234C           ERRBND    - REQUESTED ACCURACY MAX(EPSABS,EPSREL*
235C                       ABS(RESULT))
236C           *****1    - VARIABLE FOR THE LEFT SUBINTERVAL
237C           *****2    - VARIABLE FOR THE RIGHT SUBINTERVAL
238C           LAST      - INDEX FOR SUBDIVISION
239C           NRES      - NUMBER OF CALLS TO THE EXTRAPOLATION ROUTINE
240C           NUMRL2    - NUMBER OF ELEMENTS IN RLIST2. IF AN APPROPRIATE
241C                       APPROXIMATION TO THE COMPOUNDED INTEGRAL HAS
242C                       BEEN OBTAINED, IT IS PUT IN RLIST2(NUMRL2) AFTER
243C                       NUMRL2 HAS BEEN INCREASED BY ONE.
244C           ERLARG    - SUM OF THE ERRORS OVER THE INTERVALS LARGER
245C                       THAN THE SMALLEST INTERVAL CONSIDERED UP TO NOW
246C           EXTRAP    - LOGICAL VARIABLE DENOTING THAT THE ROUTINE
247C                       IS ATTEMPTING TO PERFORM EXTRAPOLATION. I.E.
248C                       BEFORE SUBDIVIDING THE SMALLEST INTERVAL WE
249C                       TRY TO DECREASE THE VALUE OF ERLARG.
250C           NOEXT     - LOGICAL VARIABLE DENOTING THAT EXTRAPOLATION IS
251C                       NO LONGER ALLOWED (TRUE-VALUE)
252C
253C            MACHINE DEPENDENT CONSTANTS
254C            ---------------------------
255C
256C           EPMACH IS THE LARGEST RELATIVE SPACING.
257C           UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
258C           OFLOW IS THE LARGEST POSITIVE MAGNITUDE.
259C
260C***FIRST EXECUTABLE STATEMENT  DQAGPE
261      EPMACH = D1MACH(4)
262C
263C            TEST ON VALIDITY OF PARAMETERS
264C            -----------------------------
265C
266      IER = 0
267      NEVAL = 0
268      LAST = 0
269      RESULT = 0.0D+00
270      ABSERR = 0.0D+00
271      ALIST(1) = A
272      BLIST(1) = B
273      RLIST(1) = 0.0D+00
274      ELIST(1) = 0.0D+00
275      IORD(1) = 0
276      LEVEL(1) = 0
277      NPTS = NPTS2-2
278      IF(NPTS2.LT.2.OR.LIMIT.LE.NPTS.OR.(EPSABS.LE.0.0D+00.AND.
279     1  EPSREL.LT.MAX(0.5D+02*EPMACH,0.5D-28))) IER = 6
280      IF(IER.EQ.6) GO TO 999
281C
282C            IF ANY BREAK POINTS ARE PROVIDED, SORT THEM INTO AN
283C            ASCENDING SEQUENCE.
284C
285      SIGN = 1.0D+00
286      IF(A.GT.B) SIGN = -1.0D+00
287      PTS(1) = MIN(A,B)
288      IF(NPTS.EQ.0) GO TO 15
289      DO 10 I = 1,NPTS
290        PTS(I+1) = POINTS(I)
291   10 CONTINUE
292   15 PTS(NPTS+2) = MAX(A,B)
293      NINT = NPTS+1
294      A1 = PTS(1)
295      IF(NPTS.EQ.0) GO TO 40
296      NINTP1 = NINT+1
297      DO 20 I = 1,NINT
298        IP1 = I+1
299        DO 20 J = IP1,NINTP1
300          IF(PTS(I).LE.PTS(J)) GO TO 20
301          TEMP = PTS(I)
302          PTS(I) = PTS(J)
303          PTS(J) = TEMP
304   20 CONTINUE
305      IF(PTS(1).NE.MIN(A,B).OR.PTS(NINTP1).NE.MAX(A,B)) IER = 6
306      IF(IER.EQ.6) GO TO 999
307C
308C            COMPUTE FIRST INTEGRAL AND ERROR APPROXIMATIONS.
309C            ------------------------------------------------
310C
311   40 RESABS = 0.0D+00
312      DO 50 I = 1,NINT
313        B1 = PTS(I+1)
314        CALL DQK21(F,A1,B1,AREA1,ERROR1,DEFABS,RESA)
315        ABSERR = ABSERR+ERROR1
316        RESULT = RESULT+AREA1
317        NDIN(I) = 0
318        IF(ERROR1.EQ.RESA.AND.ERROR1.NE.0.0D+00) NDIN(I) = 1
319        RESABS = RESABS+DEFABS
320        LEVEL(I) = 0
321        ELIST(I) = ERROR1
322        ALIST(I) = A1
323        BLIST(I) = B1
324        RLIST(I) = AREA1
325        IORD(I) = I
326        A1 = B1
327   50 CONTINUE
328      ERRSUM = 0.0D+00
329      DO 55 I = 1,NINT
330        IF(NDIN(I).EQ.1) ELIST(I) = ABSERR
331        ERRSUM = ERRSUM+ELIST(I)
332   55 CONTINUE
333C
334C           TEST ON ACCURACY.
335C
336      LAST = NINT
337      NEVAL = 21*NINT
338      DRES = ABS(RESULT)
339      ERRBND = MAX(EPSABS,EPSREL*DRES)
340      IF(ABSERR.LE.0.1D+03*EPMACH*RESABS.AND.ABSERR.GT.ERRBND) IER = 2
341      IF(NINT.EQ.1) GO TO 80
342      DO 70 I = 1,NPTS
343        JLOW = I+1
344        IND1 = IORD(I)
345        DO 60 J = JLOW,NINT
346          IND2 = IORD(J)
347          IF(ELIST(IND1).GT.ELIST(IND2)) GO TO 60
348          IND1 = IND2
349          K = J
350   60   CONTINUE
351        IF(IND1.EQ.IORD(I)) GO TO 70
352        IORD(K) = IORD(I)
353        IORD(I) = IND1
354   70 CONTINUE
355      IF(LIMIT.LT.NPTS2) IER = 1
356   80 IF(IER.NE.0.OR.ABSERR.LE.ERRBND) GO TO 210
357C
358C           INITIALIZATION
359C           --------------
360C
361      RLIST2(1) = RESULT
362      MAXERR = IORD(1)
363      ERRMAX = ELIST(MAXERR)
364      AREA = RESULT
365      NRMAX = 1
366      NRES = 0
367      NUMRL2 = 1
368      KTMIN = 0
369      EXTRAP = .FALSE.
370      NOEXT = .FALSE.
371      ERLARG = ERRSUM
372      ERTEST = ERRBND
373      LEVMAX = 1
374      IROFF1 = 0
375      IROFF2 = 0
376      IROFF3 = 0
377      IERRO = 0
378      UFLOW = D1MACH(1)
379      OFLOW = D1MACH(2)
380      ABSERR = OFLOW
381      KSGN = -1
382      IF(DRES.GE.(0.1D+01-0.5D+02*EPMACH)*RESABS) KSGN = 1
383C
384C           MAIN DO-LOOP
385C           ------------
386C
387      DO 160 LAST = NPTS2,LIMIT
388C
389C           BISECT THE SUBINTERVAL WITH THE NRMAX-TH LARGEST ERROR
390C           ESTIMATE.
391C
392        LEVCUR = LEVEL(MAXERR)+1
393        A1 = ALIST(MAXERR)
394        B1 = 0.5D+00*(ALIST(MAXERR)+BLIST(MAXERR))
395        A2 = B1
396        B2 = BLIST(MAXERR)
397        ERLAST = ERRMAX
398        CALL DQK21(F,A1,B1,AREA1,ERROR1,RESA,DEFAB1)
399        CALL DQK21(F,A2,B2,AREA2,ERROR2,RESA,DEFAB2)
400C
401C           IMPROVE PREVIOUS APPROXIMATIONS TO INTEGRAL
402C           AND ERROR AND TEST FOR ACCURACY.
403C
404        NEVAL = NEVAL+42
405        AREA12 = AREA1+AREA2
406        ERRO12 = ERROR1+ERROR2
407        ERRSUM = ERRSUM+ERRO12-ERRMAX
408        AREA = AREA+AREA12-RLIST(MAXERR)
409        IF(DEFAB1.EQ.ERROR1.OR.DEFAB2.EQ.ERROR2) GO TO 95
410        IF(ABS(RLIST(MAXERR)-AREA12).GT.0.1D-04*ABS(AREA12)
411     1  .OR.ERRO12.LT.0.99D+00*ERRMAX) GO TO 90
412        IF(EXTRAP) IROFF2 = IROFF2+1
413        IF(.NOT.EXTRAP) IROFF1 = IROFF1+1
414   90   IF(LAST.GT.10.AND.ERRO12.GT.ERRMAX) IROFF3 = IROFF3+1
415   95   LEVEL(MAXERR) = LEVCUR
416        LEVEL(LAST) = LEVCUR
417        RLIST(MAXERR) = AREA1
418        RLIST(LAST) = AREA2
419        ERRBND = MAX(EPSABS,EPSREL*ABS(AREA))
420C
421C           TEST FOR ROUNDOFF ERROR AND EVENTUALLY SET ERROR FLAG.
422C
423        IF(IROFF1+IROFF2.GE.10.OR.IROFF3.GE.20) IER = 2
424        IF(IROFF2.GE.5) IERRO = 3
425C
426C           SET ERROR FLAG IN THE CASE THAT THE NUMBER OF
427C           SUBINTERVALS EQUALS LIMIT.
428C
429        IF(LAST.EQ.LIMIT) IER = 1
430C
431C           SET ERROR FLAG IN THE CASE OF BAD INTEGRAND BEHAVIOUR
432C           AT A POINT OF THE INTEGRATION RANGE
433C
434        IF(MAX(ABS(A1),ABS(B2)).LE.(0.1D+01+0.1D+03*EPMACH)*
435     1  (ABS(A2)+0.1D+04*UFLOW)) IER = 4
436C
437C           APPEND THE NEWLY-CREATED INTERVALS TO THE LIST.
438C
439        IF(ERROR2.GT.ERROR1) GO TO 100
440        ALIST(LAST) = A2
441        BLIST(MAXERR) = B1
442        BLIST(LAST) = B2
443        ELIST(MAXERR) = ERROR1
444        ELIST(LAST) = ERROR2
445        GO TO 110
446  100   ALIST(MAXERR) = A2
447        ALIST(LAST) = A1
448        BLIST(LAST) = B1
449        RLIST(MAXERR) = AREA2
450        RLIST(LAST) = AREA1
451        ELIST(MAXERR) = ERROR2
452        ELIST(LAST) = ERROR1
453C
454C           CALL SUBROUTINE DQPSRT TO MAINTAIN THE DESCENDING ORDERING
455C           IN THE LIST OF ERROR ESTIMATES AND SELECT THE SUBINTERVAL
456C           WITH NRMAX-TH LARGEST ERROR ESTIMATE (TO BE BISECTED NEXT).
457C
458  110   CALL DQPSRT(LIMIT,LAST,MAXERR,ERRMAX,ELIST,IORD,NRMAX)
459C ***JUMP OUT OF DO-LOOP
460        IF(ERRSUM.LE.ERRBND) GO TO 190
461C ***JUMP OUT OF DO-LOOP
462        IF(IER.NE.0) GO TO 170
463        IF(NOEXT) GO TO 160
464        ERLARG = ERLARG-ERLAST
465        IF(LEVCUR+1.LE.LEVMAX) ERLARG = ERLARG+ERRO12
466        IF(EXTRAP) GO TO 120
467C
468C           TEST WHETHER THE INTERVAL TO BE BISECTED NEXT IS THE
469C           SMALLEST INTERVAL.
470C
471        IF(LEVEL(MAXERR)+1.LE.LEVMAX) GO TO 160
472        EXTRAP = .TRUE.
473        NRMAX = 2
474  120   IF(IERRO.EQ.3.OR.ERLARG.LE.ERTEST) GO TO 140
475C
476C           THE SMALLEST INTERVAL HAS THE LARGEST ERROR.
477C           BEFORE BISECTING DECREASE THE SUM OF THE ERRORS OVER
478C           THE LARGER INTERVALS (ERLARG) AND PERFORM EXTRAPOLATION.
479C
480        ID = NRMAX
481        JUPBND = LAST
482        IF(LAST.GT.(2+LIMIT/2)) JUPBND = LIMIT+3-LAST
483        DO 130 K = ID,JUPBND
484          MAXERR = IORD(NRMAX)
485          ERRMAX = ELIST(MAXERR)
486C ***JUMP OUT OF DO-LOOP
487          IF(LEVEL(MAXERR)+1.LE.LEVMAX) GO TO 160
488          NRMAX = NRMAX+1
489  130   CONTINUE
490C
491C           PERFORM EXTRAPOLATION.
492C
493  140   NUMRL2 = NUMRL2+1
494        RLIST2(NUMRL2) = AREA
495        IF(NUMRL2.LE.2) GO TO 155
496        CALL DQELG(NUMRL2,RLIST2,RESEPS,ABSEPS,RES3LA,NRES)
497        KTMIN = KTMIN+1
498        IF(KTMIN.GT.5.AND.ABSERR.LT.0.1D-02*ERRSUM) IER = 5
499        IF(ABSEPS.GE.ABSERR) GO TO 150
500        KTMIN = 0
501        ABSERR = ABSEPS
502        RESULT = RESEPS
503        CORREC = ERLARG
504        ERTEST = MAX(EPSABS,EPSREL*ABS(RESEPS))
505C ***JUMP OUT OF DO-LOOP
506        IF(ABSERR.LT.ERTEST) GO TO 170
507C
508C           PREPARE BISECTION OF THE SMALLEST INTERVAL.
509C
510  150   IF(NUMRL2.EQ.1) NOEXT = .TRUE.
511        IF(IER.GE.5) GO TO 170
512  155   MAXERR = IORD(1)
513        ERRMAX = ELIST(MAXERR)
514        NRMAX = 1
515        EXTRAP = .FALSE.
516        LEVMAX = LEVMAX+1
517        ERLARG = ERRSUM
518  160 CONTINUE
519C
520C           SET THE FINAL RESULT.
521C           ---------------------
522C
523C
524  170 IF(ABSERR.EQ.OFLOW) GO TO 190
525      IF((IER+IERRO).EQ.0) GO TO 180
526      IF(IERRO.EQ.3) ABSERR = ABSERR+CORREC
527      IF(IER.EQ.0) IER = 3
528      IF(RESULT.NE.0.0D+00.AND.AREA.NE.0.0D+00)GO TO 175
529      IF(ABSERR.GT.ERRSUM)GO TO 190
530      IF(AREA.EQ.0.0D+00) GO TO 210
531      GO TO 180
532  175 IF(ABSERR/ABS(RESULT).GT.ERRSUM/ABS(AREA))GO TO 190
533C
534C           TEST ON DIVERGENCE.
535C
536  180 IF(KSGN.EQ.(-1).AND.MAX(ABS(RESULT),ABS(AREA)).LE.
537     1  RESABS*0.1D-01) GO TO 210
538      IF(0.1D-01.GT.(RESULT/AREA).OR.(RESULT/AREA).GT.0.1D+03.OR.
539     1  ERRSUM.GT.ABS(AREA)) IER = 6
540      GO TO 210
541C
542C           COMPUTE GLOBAL INTEGRAL SUM.
543C
544  190 RESULT = 0.0D+00
545      DO 200 K = 1,LAST
546        RESULT = RESULT+RLIST(K)
547  200 CONTINUE
548      ABSERR = ERRSUM
549  210 IF(IER.GT.2) IER = IER-1
550      RESULT = RESULT*SIGN
551  999 RETURN
552      END
553