GNU Octave  4.4.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
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
d1mach
double precision function d1mach(i)
Definition: d1mach.f:20
dimension
OCTAVE_EXPORT octave_value_list etc The functions then dimension(columns)