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
dqagpe.f
Go to the documentation of this file.
1  SUBROUTINE dqagpe(F,A,B,NPTS2,POINTS,EPSABS,EPSREL,LIMIT,RESULT,
2  * abserr,neval,ier,alist,blist,rlist,elist,pts,iord,level,ndin,
3  * last)
4 C***BEGIN PROLOGUE DQAGPE
5 C***DATE WRITTEN 800101 (YYMMDD)
6 C***REVISION DATE 830518 (YYMMDD)
7 C***CATEGORY NO. H2A2A1
8 C***KEYWORDS AUTOMATIC INTEGRATOR, GENERAL-PURPOSE,
9 C SINGULARITIES AT USER SPECIFIED POINTS,
10 C EXTRAPOLATION, GLOBALLY ADAPTIVE.
11 C***AUTHOR PIESSENS,ROBERT ,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN
12 C DE DONCKER,ELISE,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN
13 C***PURPOSE THE ROUTINE CALCULATES AN APPROXIMATION RESULT TO A GIVEN
14 C DEFINITE INTEGRAL I = INTEGRAL OF F OVER (A,B), HOPEFULLY
15 C SATISFYING FOLLOWING CLAIM FOR ACCURACY ABS(I-RESULT).LE.
16 C MAX(EPSABS,EPSREL*ABS(I)). BREAK POINTS OF THE INTEGRATION
17 C INTERVAL, WHERE LOCAL DIFFICULTIES OF THE INTEGRAND MAY
18 C OCCUR(E.G. SINGULARITIES,DISCONTINUITIES),PROVIDED BY USER.
19 C***DESCRIPTION
20 C
21 C COMPUTATION OF A DEFINITE INTEGRAL
22 C STANDARD FORTRAN SUBROUTINE
23 C DOUBLE PRECISION VERSION
24 C
25 C PARAMETERS
26 C ON ENTRY
27 C F - SUBROUTINE F(X,IERR,RESULT) DEFINING THE INTEGRAND
28 C FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
29 C DECLARED E X T E R N A L IN THE DRIVER PROGRAM.
30 C
31 C A - DOUBLE PRECISION
32 C LOWER LIMIT OF INTEGRATION
33 C
34 C B - DOUBLE PRECISION
35 C UPPER LIMIT OF INTEGRATION
36 C
37 C NPTS2 - INTEGER
38 C NUMBER EQUAL TO TWO MORE THAN THE NUMBER OF
39 C USER-SUPPLIED BREAK POINTS WITHIN THE INTEGRATION
40 C RANGE, NPTS2.GE.2.
41 C IF NPTS2.LT.2, THE ROUTINE WILL END WITH IER = 6.
42 C
43 C POINTS - DOUBLE PRECISION
44 C VECTOR OF DIMENSION NPTS2, THE FIRST (NPTS2-2)
45 C ELEMENTS OF WHICH ARE THE USER PROVIDED BREAK
46 C POINTS. IF THESE POINTS DO NOT CONSTITUTE AN
47 C ASCENDING SEQUENCE THERE WILL BE AN AUTOMATIC
48 C SORTING.
49 C
50 C EPSABS - DOUBLE PRECISION
51 C ABSOLUTE ACCURACY REQUESTED
52 C EPSREL - DOUBLE PRECISION
53 C RELATIVE ACCURACY REQUESTED
54 C IF EPSABS.LE.0
55 C AND EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28),
56 C THE ROUTINE WILL END WITH IER = 6.
57 C
58 C LIMIT - INTEGER
59 C GIVES AN UPPER BOUND ON THE NUMBER OF SUBINTERVALS
60 C IN THE PARTITION OF (A,B), LIMIT.GE.NPTS2
61 C IF LIMIT.LT.NPTS2, THE ROUTINE WILL END WITH
62 C IER = 6.
63 C
64 C ON RETURN
65 C RESULT - DOUBLE PRECISION
66 C APPROXIMATION TO THE INTEGRAL
67 C
68 C ABSERR - DOUBLE PRECISION
69 C ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR,
70 C WHICH SHOULD EQUAL OR EXCEED ABS(I-RESULT)
71 C
72 C NEVAL - INTEGER
73 C NUMBER OF INTEGRAND EVALUATIONS
74 C
75 C IER - INTEGER
76 C IER = 0 NORMAL AND RELIABLE TERMINATION OF THE
77 C ROUTINE. IT IS ASSUMED THAT THE REQUESTED
78 C ACCURACY HAS BEEN ACHIEVED.
79 C IER.GT.0 ABNORMAL TERMINATION OF THE ROUTINE.
80 C THE ESTIMATES FOR INTEGRAL AND ERROR ARE
81 C LESS RELIABLE. IT IS ASSUMED THAT THE
82 C REQUESTED ACCURACY HAS NOT BEEN ACHIEVED.
83 C IER.LT.0 EXIT REQUESTED FROM USER-SUPPLIED
84 C FUNCTION.
85 C
86 C ERROR MESSAGES
87 C IER = 1 MAXIMUM NUMBER OF SUBDIVISIONS ALLOWED
88 C HAS BEEN ACHIEVED. ONE CAN ALLOW MORE
89 C SUBDIVISIONS BY INCREASING THE VALUE OF
90 C LIMIT (AND TAKING THE ACCORDING DIMENSION
91 C ADJUSTMENTS INTO ACCOUNT). HOWEVER, IF
92 C THIS YIELDS NO IMPROVEMENT IT IS ADVISED
93 C TO ANALYZE THE INTEGRAND IN ORDER TO
94 C DETERMINE THE INTEGRATION DIFFICULTIES. IF
95 C THE POSITION OF A LOCAL DIFFICULTY CAN BE
96 C DETERMINED (I.E. SINGULARITY,
97 C DISCONTINUITY WITHIN THE INTERVAL), IT
98 C SHOULD BE SUPPLIED TO THE ROUTINE AS AN
99 C ELEMENT OF THE VECTOR POINTS. IF NECESSARY
100 C AN APPROPRIATE SPECIAL-PURPOSE INTEGRATOR
101 C MUST BE USED, WHICH IS DESIGNED FOR
102 C HANDLING THE TYPE OF DIFFICULTY INVOLVED.
103 C = 2 THE OCCURRENCE OF ROUNDOFF ERROR IS
104 C DETECTED, WHICH PREVENTS THE REQUESTED
105 C TOLERANCE FROM BEING ACHIEVED.
106 C THE ERROR MAY BE UNDER-ESTIMATED.
107 C = 3 EXTREMELY BAD INTEGRAND BEHAVIOUR OCCURS
108 C AT SOME POINTS OF THE INTEGRATION
109 C INTERVAL.
110 C = 4 THE ALGORITHM DOES NOT CONVERGE.
111 C ROUNDOFF ERROR IS DETECTED IN THE
112 C EXTRAPOLATION TABLE. IT IS PRESUMED THAT
113 C THE REQUESTED TOLERANCE CANNOT BE
114 C ACHIEVED, AND THAT THE RETURNED RESULT IS
115 C THE BEST WHICH CAN BE OBTAINED.
116 C = 5 THE INTEGRAL IS PROBABLY DIVERGENT, OR
117 C SLOWLY CONVERGENT. IT MUST BE NOTED THAT
118 C DIVERGENCE CAN OCCUR WITH ANY OTHER VALUE
119 C OF IER.GT.0.
120 C = 6 THE INPUT IS INVALID BECAUSE
121 C NPTS2.LT.2 OR
122 C BREAK POINTS ARE SPECIFIED OUTSIDE
123 C THE INTEGRATION RANGE OR
124 C (EPSABS.LE.0 AND
125 C EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28))
126 C OR LIMIT.LT.NPTS2.
127 C RESULT, ABSERR, NEVAL, LAST, RLIST(1),
128 C AND ELIST(1) ARE SET TO ZERO. ALIST(1) AND
129 C BLIST(1) ARE SET TO A AND B RESPECTIVELY.
130 C
131 C ALIST - DOUBLE PRECISION
132 C VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
133 C LAST ELEMENTS OF WHICH ARE THE LEFT END POINTS
134 C OF THE SUBINTERVALS IN THE PARTITION OF THE GIVEN
135 C INTEGRATION RANGE (A,B)
136 C
137 C BLIST - DOUBLE PRECISION
138 C VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
139 C LAST ELEMENTS OF WHICH ARE THE RIGHT END POINTS
140 C OF THE SUBINTERVALS IN THE PARTITION OF THE GIVEN
141 C INTEGRATION RANGE (A,B)
142 C
143 C RLIST - DOUBLE PRECISION
144 C VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
145 C LAST ELEMENTS OF WHICH ARE THE INTEGRAL
146 C APPROXIMATIONS ON THE SUBINTERVALS
147 C
148 C ELIST - DOUBLE PRECISION
149 C VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
150 C LAST ELEMENTS OF WHICH ARE THE MODULI OF THE
151 C ABSOLUTE ERROR ESTIMATES ON THE SUBINTERVALS
152 C
153 C PTS - DOUBLE PRECISION
154 C VECTOR OF DIMENSION AT LEAST NPTS2, CONTAINING THE
155 C INTEGRATION LIMITS AND THE BREAK POINTS OF THE
156 C INTERVAL IN ASCENDING SEQUENCE.
157 C
158 C LEVEL - INTEGER
159 C VECTOR OF DIMENSION AT LEAST LIMIT, CONTAINING THE
160 C SUBDIVISION LEVELS OF THE SUBINTERVAL, I.E. IF
161 C (AA,BB) IS A SUBINTERVAL OF (P1,P2) WHERE P1 AS
162 C WELL AS P2 IS A USER-PROVIDED BREAK POINT OR
163 C INTEGRATION LIMIT, THEN (AA,BB) HAS LEVEL L IF
164 C ABS(BB-AA) = ABS(P2-P1)*2**(-L).
165 C
166 C NDIN - INTEGER
167 C VECTOR OF DIMENSION AT LEAST NPTS2, AFTER FIRST
168 C INTEGRATION OVER THE INTERVALS (PTS(I)),PTS(I+1),
169 C I = 0,1, ..., NPTS2-2, THE ERROR ESTIMATES OVER
170 C SOME OF THE INTERVALS MAY HAVE BEEN INCREASED
171 C ARTIFICIALLY, IN ORDER TO PUT THEIR SUBDIVISION
172 C FORWARD. IF THIS HAPPENS FOR THE SUBINTERVAL
173 C NUMBERED K, NDIN(K) IS PUT TO 1, OTHERWISE
174 C NDIN(K) = 0.
175 C
176 C IORD - INTEGER
177 C VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST K
178 C ELEMENTS OF WHICH ARE POINTERS TO THE
179 C ERROR ESTIMATES OVER THE SUBINTERVALS,
180 C SUCH THAT ELIST(IORD(1)), ..., ELIST(IORD(K))
181 C FORM A DECREASING SEQUENCE, WITH K = LAST
182 C IF LAST.LE.(LIMIT/2+2), AND K = LIMIT+1-LAST
183 C OTHERWISE
184 C
185 C LAST - INTEGER
186 C NUMBER OF SUBINTERVALS ACTUALLY PRODUCED IN THE
187 C SUBDIVISIONS PROCESS
188 C
189 C***REFERENCES (NONE)
190 C***ROUTINES CALLED D1MACH,DQELG,DQK21,DQPSRT
191 C***END PROLOGUE DQAGPE
192  DOUBLE PRECISION a,abseps,abserr,alist,area,area1,area12,area2,a1,
193  * a2,b,blist,b1,b2,correc,dabs,defabs,defab1,defab2,dmax1,dmin1,
194  * dres,d1mach,elist,epmach,epsabs,epsrel,erlarg,erlast,errbnd,
195  * errmax,error1,erro12,error2,errsum,ertest,oflow,points,pts,
196  * 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  * jlow,jupbnd,k,ksgn,ktmin,last,levcur,level,levmax,limit,maxerr,
199  * ndin,neval,nint,nintp1,npts,npts2,nres,nrmax,numrl2
200  LOGICAL extrap,noext
201 C
202 C
203  dimension alist(limit),blist(limit),elist(limit),iord(limit),
204  * level(limit),ndin(npts2),points(npts2),pts(npts2),res3la(3),
205  * rlist(limit),rlist2(52)
206 C
207  EXTERNAL f
208 C
209 C THE DIMENSION OF RLIST2 IS DETERMINED BY THE VALUE OF
210 C LIMEXP IN SUBROUTINE EPSALG (RLIST2 SHOULD BE OF DIMENSION
211 C (LIMEXP+2) AT LEAST).
212 C
213 C
214 C LIST OF MAJOR VARIABLES
215 C -----------------------
216 C
217 C ALIST - LIST OF LEFT END POINTS OF ALL SUBINTERVALS
218 C CONSIDERED UP TO NOW
219 C BLIST - LIST OF RIGHT END POINTS OF ALL SUBINTERVALS
220 C CONSIDERED UP TO NOW
221 C RLIST(I) - APPROXIMATION TO THE INTEGRAL OVER
222 C (ALIST(I),BLIST(I))
223 C RLIST2 - ARRAY OF DIMENSION AT LEAST LIMEXP+2
224 C CONTAINING THE PART OF THE EPSILON TABLE WHICH
225 C IS STILL NEEDED FOR FURTHER COMPUTATIONS
226 C ELIST(I) - ERROR ESTIMATE APPLYING TO RLIST(I)
227 C MAXERR - POINTER TO THE INTERVAL WITH LARGEST ERROR
228 C ESTIMATE
229 C ERRMAX - ELIST(MAXERR)
230 C ERLAST - ERROR ON THE INTERVAL CURRENTLY SUBDIVIDED
231 C (BEFORE THAT SUBDIVISION HAS TAKEN PLACE)
232 C AREA - SUM OF THE INTEGRALS OVER THE SUBINTERVALS
233 C ERRSUM - SUM OF THE ERRORS OVER THE SUBINTERVALS
234 C ERRBND - REQUESTED ACCURACY MAX(EPSABS,EPSREL*
235 C ABS(RESULT))
236 C *****1 - VARIABLE FOR THE LEFT SUBINTERVAL
237 C *****2 - VARIABLE FOR THE RIGHT SUBINTERVAL
238 C LAST - INDEX FOR SUBDIVISION
239 C NRES - NUMBER OF CALLS TO THE EXTRAPOLATION ROUTINE
240 C NUMRL2 - NUMBER OF ELEMENTS IN RLIST2. IF AN APPROPRIATE
241 C APPROXIMATION TO THE COMPOUNDED INTEGRAL HAS
242 C BEEN OBTAINED, IT IS PUT IN RLIST2(NUMRL2) AFTER
243 C NUMRL2 HAS BEEN INCREASED BY ONE.
244 C ERLARG - SUM OF THE ERRORS OVER THE INTERVALS LARGER
245 C THAN THE SMALLEST INTERVAL CONSIDERED UP TO NOW
246 C EXTRAP - LOGICAL VARIABLE DENOTING THAT THE ROUTINE
247 C IS ATTEMPTING TO PERFORM EXTRAPOLATION. I.E.
248 C BEFORE SUBDIVIDING THE SMALLEST INTERVAL WE
249 C TRY TO DECREASE THE VALUE OF ERLARG.
250 C NOEXT - LOGICAL VARIABLE DENOTING THAT EXTRAPOLATION IS
251 C NO LONGER ALLOWED (TRUE-VALUE)
252 C
253 C MACHINE DEPENDENT CONSTANTS
254 C ---------------------------
255 C
256 C EPMACH IS THE LARGEST RELATIVE SPACING.
257 C UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
258 C OFLOW IS THE LARGEST POSITIVE MAGNITUDE.
259 C
260 C***FIRST EXECUTABLE STATEMENT DQAGPE
261  epmach = d1mach(4)
262 C
263 C TEST ON VALIDITY OF PARAMETERS
264 C -----------------------------
265 C
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  * epsrel.LT.dmax1(0.5d+02*epmach,0.5d-28))) ier = 6
280  IF(ier.EQ.6) go to 999
281 C
282 C IF ANY BREAK POINTS ARE PROVIDED, SORT THEM INTO AN
283 C ASCENDING SEQUENCE.
284 C
285  sign = 1.0d+00
286  IF(a.GT.b) sign = -1.0d+00
287  pts(1) = dmin1(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) = dmax1(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.dmin1(a,b).OR.pts(nintp1).NE.dmax1(a,b)) ier = 6
306  IF(ier.EQ.6) go to 999
307 C
308 C COMPUTE FIRST INTEGRAL AND ERROR APPROXIMATIONS.
309 C ------------------------------------------------
310 C
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,ier)
315  IF (ier .LT. 0) RETURN
316  abserr = abserr+error1
317  result = result+area1
318  ndin(i) = 0
319  IF(error1.EQ.resa.AND.error1.NE.0.0d+00) ndin(i) = 1
320  resabs = resabs+defabs
321  level(i) = 0
322  elist(i) = error1
323  alist(i) = a1
324  blist(i) = b1
325  rlist(i) = area1
326  iord(i) = i
327  a1 = b1
328  50 CONTINUE
329  errsum = 0.0d+00
330  DO 55 i = 1,nint
331  IF(ndin(i).EQ.1) elist(i) = abserr
332  errsum = errsum+elist(i)
333  55 CONTINUE
334 C
335 C TEST ON ACCURACY.
336 C
337  last = nint
338  neval = 21*nint
339  dres = dabs(result)
340  errbnd = dmax1(epsabs,epsrel*dres)
341  IF(abserr.LE.0.1d+03*epmach*resabs.AND.abserr.GT.errbnd) ier = 2
342  IF(nint.EQ.1) go to 80
343  DO 70 i = 1,npts
344  jlow = i+1
345  ind1 = iord(i)
346  DO 60 j = jlow,nint
347  ind2 = iord(j)
348  IF(elist(ind1).GT.elist(ind2)) go to 60
349  ind1 = ind2
350  k = j
351  60 CONTINUE
352  IF(ind1.EQ.iord(i)) go to 70
353  iord(k) = iord(i)
354  iord(i) = ind1
355  70 CONTINUE
356  IF(limit.LT.npts2) ier = 1
357  80 IF(ier.NE.0.OR.abserr.LE.errbnd) go to 210
358 C
359 C INITIALIZATION
360 C --------------
361 C
362  rlist2(1) = result
363  maxerr = iord(1)
364  errmax = elist(maxerr)
365  area = result
366  nrmax = 1
367  nres = 0
368  numrl2 = 1
369  ktmin = 0
370  extrap = .false.
371  noext = .false.
372  erlarg = errsum
373  ertest = errbnd
374  levmax = 1
375  iroff1 = 0
376  iroff2 = 0
377  iroff3 = 0
378  ierro = 0
379  uflow = d1mach(1)
380  oflow = d1mach(2)
381  abserr = oflow
382  ksgn = -1
383  IF(dres.GE.(0.1d+01-0.5d+02*epmach)*resabs) ksgn = 1
384 C
385 C MAIN DO-LOOP
386 C ------------
387 C
388  DO 160 last = npts2,limit
389 C
390 C BISECT THE SUBINTERVAL WITH THE NRMAX-TH LARGEST ERROR
391 C ESTIMATE.
392 C
393  levcur = level(maxerr)+1
394  a1 = alist(maxerr)
395  b1 = 0.5d+00*(alist(maxerr)+blist(maxerr))
396  a2 = b1
397  b2 = blist(maxerr)
398  erlast = errmax
399  CALL dqk21(f,a1,b1,area1,error1,resa,defab1,ier)
400  IF (ier .LT. 0) RETURN
401  CALL dqk21(f,a2,b2,area2,error2,resa,defab2,ier)
402  IF (ier .LT. 0) RETURN
403 C
404 C IMPROVE PREVIOUS APPROXIMATIONS TO INTEGRAL
405 C AND ERROR AND TEST FOR ACCURACY.
406 C
407  neval = neval+42
408  area12 = area1+area2
409  erro12 = error1+error2
410  errsum = errsum+erro12-errmax
411  area = area+area12-rlist(maxerr)
412  IF(defab1.EQ.error1.OR.defab2.EQ.error2) go to 95
413  IF(dabs(rlist(maxerr)-area12).GT.0.1d-04*dabs(area12)
414  * .OR.erro12.LT.0.99d+00*errmax) go to 90
415  IF(extrap) iroff2 = iroff2+1
416  IF(.NOT.extrap) iroff1 = iroff1+1
417  90 IF(last.GT.10.AND.erro12.GT.errmax) iroff3 = iroff3+1
418  95 level(maxerr) = levcur
419  level(last) = levcur
420  rlist(maxerr) = area1
421  rlist(last) = area2
422  errbnd = dmax1(epsabs,epsrel*dabs(area))
423 C
424 C TEST FOR ROUNDOFF ERROR AND EVENTUALLY SET ERROR FLAG.
425 C
426  IF(iroff1+iroff2.GE.10.OR.iroff3.GE.20) ier = 2
427  IF(iroff2.GE.5) ierro = 3
428 C
429 C SET ERROR FLAG IN THE CASE THAT THE NUMBER OF
430 C SUBINTERVALS EQUALS LIMIT.
431 C
432  IF(last.EQ.limit) ier = 1
433 C
434 C SET ERROR FLAG IN THE CASE OF BAD INTEGRAND BEHAVIOUR
435 C AT A POINT OF THE INTEGRATION RANGE
436 C
437  IF(dmax1(dabs(a1),dabs(b2)).LE.(0.1d+01+0.1d+03*epmach)*
438  * (dabs(a2)+0.1d+04*uflow)) ier = 4
439 C
440 C APPEND THE NEWLY-CREATED INTERVALS TO THE LIST.
441 C
442  IF(error2.GT.error1) go to 100
443  alist(last) = a2
444  blist(maxerr) = b1
445  blist(last) = b2
446  elist(maxerr) = error1
447  elist(last) = error2
448  go to 110
449  100 alist(maxerr) = a2
450  alist(last) = a1
451  blist(last) = b1
452  rlist(maxerr) = area2
453  rlist(last) = area1
454  elist(maxerr) = error2
455  elist(last) = error1
456 C
457 C CALL SUBROUTINE DQPSRT TO MAINTAIN THE DESCENDING ORDERING
458 C IN THE LIST OF ERROR ESTIMATES AND SELECT THE SUBINTERVAL
459 C WITH NRMAX-TH LARGEST ERROR ESTIMATE (TO BE BISECTED NEXT).
460 C
461  110 CALL dqpsrt(limit,last,maxerr,errmax,elist,iord,nrmax)
462 C ***JUMP OUT OF DO-LOOP
463  IF(errsum.LE.errbnd) go to 190
464 C ***JUMP OUT OF DO-LOOP
465  IF(ier.NE.0) go to 170
466  IF(noext) go to 160
467  erlarg = erlarg-erlast
468  IF(levcur+1.LE.levmax) erlarg = erlarg+erro12
469  IF(extrap) go to 120
470 C
471 C TEST WHETHER THE INTERVAL TO BE BISECTED NEXT IS THE
472 C SMALLEST INTERVAL.
473 C
474  IF(level(maxerr)+1.LE.levmax) go to 160
475  extrap = .true.
476  nrmax = 2
477  120 IF(ierro.EQ.3.OR.erlarg.LE.ertest) go to 140
478 C
479 C THE SMALLEST INTERVAL HAS THE LARGEST ERROR.
480 C BEFORE BISECTING DECREASE THE SUM OF THE ERRORS OVER
481 C THE LARGER INTERVALS (ERLARG) AND PERFORM EXTRAPOLATION.
482 C
483  id = nrmax
484  jupbnd = last
485  IF(last.GT.(2+limit/2)) jupbnd = limit+3-last
486  DO 130 k = id,jupbnd
487  maxerr = iord(nrmax)
488  errmax = elist(maxerr)
489 C ***JUMP OUT OF DO-LOOP
490  IF(level(maxerr)+1.LE.levmax) go to 160
491  nrmax = nrmax+1
492  130 CONTINUE
493 C
494 C PERFORM EXTRAPOLATION.
495 C
496  140 numrl2 = numrl2+1
497  rlist2(numrl2) = area
498  IF(numrl2.LE.2) go to 155
499  CALL dqelg(numrl2,rlist2,reseps,abseps,res3la,nres)
500  ktmin = ktmin+1
501  IF(ktmin.GT.5.AND.abserr.LT.0.1d-02*errsum) ier = 5
502  IF(abseps.GE.abserr) go to 150
503  ktmin = 0
504  abserr = abseps
505  result = reseps
506  correc = erlarg
507  ertest = dmax1(epsabs,epsrel*dabs(reseps))
508 C ***JUMP OUT OF DO-LOOP
509  IF(abserr.LT.ertest) go to 170
510 C
511 C PREPARE BISECTION OF THE SMALLEST INTERVAL.
512 C
513  150 IF(numrl2.EQ.1) noext = .true.
514  IF(ier.GE.5) go to 170
515  155 maxerr = iord(1)
516  errmax = elist(maxerr)
517  nrmax = 1
518  extrap = .false.
519  levmax = levmax+1
520  erlarg = errsum
521  160 CONTINUE
522 C
523 C SET THE FINAL RESULT.
524 C ---------------------
525 C
526 C
527  170 IF(abserr.EQ.oflow) go to 190
528  IF((ier+ierro).EQ.0) go to 180
529  IF(ierro.EQ.3) abserr = abserr+correc
530  IF(ier.EQ.0) ier = 3
531  IF(result.NE.0.0d+00.AND.area.NE.0.0d+00)go to 175
532  IF(abserr.GT.errsum)go to 190
533  IF(area.EQ.0.0d+00) go to 210
534  go to 180
535  175 IF(abserr/dabs(result).GT.errsum/dabs(area))go to 190
536 C
537 C TEST ON DIVERGENCE.
538 C
539  180 IF(ksgn.EQ.(-1).AND.dmax1(dabs(result),dabs(area)).LE.
540  * resabs*0.1d-01) go to 210
541  IF(0.1d-01.GT.(result/area).OR.(result/area).GT.0.1d+03.OR.
542  * errsum.GT.dabs(area)) ier = 6
543  go to 210
544 C
545 C COMPUTE GLOBAL INTEGRAL SUM.
546 C
547  190 result = 0.0d+00
548  DO 200 k = 1,last
549  result = result+rlist(k)
550  200 CONTINUE
551  abserr = errsum
552  210 IF(ier.GT.2) ier = ier-1
553  result = result*sign
554  999 RETURN
555  END
Generated on Mon Dec 30 2013 03:04:47 for GNU Octave by   doxygen 1.8.1.2