|
2329
|
1 SUBROUTINE DQAGIE(F,BOUND,INF,EPSABS,EPSREL,LIMIT,RESULT,ABSERR, |
|
|
2 * NEVAL,IER,ALIST,BLIST,RLIST,ELIST,IORD,LAST) |
|
|
3 C***BEGIN PROLOGUE DQAGIE |
|
|
4 C***DATE WRITTEN 800101 (YYMMDD) |
|
|
5 C***REVISION DATE 830518 (YYMMDD) |
|
|
6 C***CATEGORY NO. H2A3A1,H2A4A1 |
|
|
7 C***KEYWORDS AUTOMATIC INTEGRATOR, INFINITE INTERVALS, |
|
|
8 C GENERAL-PURPOSE, TRANSFORMATION, 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 INTEGRAL I = INTEGRAL OF F OVER (BOUND,+INFINITY) |
|
|
14 C OR I = INTEGRAL OF F OVER (-INFINITY,BOUND) |
|
|
15 C OR I = INTEGRAL OF F OVER (-INFINITY,+INFINITY), |
|
|
16 C HOPEFULLY SATISFYING FOLLOWING CLAIM FOR ACCURACY |
|
|
17 C ABS(I-RESULT).LE.MAX(EPSABS,EPSREL*ABS(I)) |
|
|
18 C***DESCRIPTION |
|
|
19 C |
|
|
20 C INTEGRATION OVER INFINITE INTERVALS |
|
|
21 C STANDARD FORTRAN SUBROUTINE |
|
|
22 C |
|
3135
|
23 C F - SUBROUTINE F(X,IERR,RESULT) DEFINING THE INTEGRAND |
|
2329
|
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. |
|
|
26 C |
|
|
27 C BOUND - DOUBLE PRECISION |
|
|
28 C FINITE BOUND OF INTEGRATION RANGE |
|
|
29 C (HAS NO MEANING IF INTERVAL IS DOUBLY-INFINITE) |
|
|
30 C |
|
|
31 C INF - DOUBLE PRECISION |
|
|
32 C INDICATING THE KIND OF INTEGRATION RANGE INVOLVED |
|
|
33 C INF = 1 CORRESPONDS TO (BOUND,+INFINITY), |
|
|
34 C INF = -1 TO (-INFINITY,BOUND), |
|
|
35 C INF = 2 TO (-INFINITY,+INFINITY). |
|
|
36 C |
|
|
37 C EPSABS - DOUBLE PRECISION |
|
|
38 C ABSOLUTE ACCURACY REQUESTED |
|
|
39 C EPSREL - DOUBLE PRECISION |
|
|
40 C RELATIVE ACCURACY REQUESTED |
|
|
41 C IF EPSABS.LE.0 |
|
|
42 C AND EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28), |
|
|
43 C THE ROUTINE WILL END WITH IER = 6. |
|
|
44 C |
|
|
45 C LIMIT - INTEGER |
|
|
46 C GIVES AN UPPER BOUND ON THE NUMBER OF SUBINTERVALS |
|
|
47 C IN THE PARTITION OF (A,B), LIMIT.GE.1 |
|
|
48 C |
|
|
49 C ON RETURN |
|
|
50 C RESULT - DOUBLE PRECISION |
|
|
51 C APPROXIMATION TO THE INTEGRAL |
|
|
52 C |
|
|
53 C ABSERR - DOUBLE PRECISION |
|
|
54 C ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR, |
|
|
55 C WHICH SHOULD EQUAL OR EXCEED ABS(I-RESULT) |
|
|
56 C |
|
|
57 C NEVAL - INTEGER |
|
|
58 C NUMBER OF INTEGRAND EVALUATIONS |
|
|
59 C |
|
|
60 C IER - INTEGER |
|
|
61 C IER = 0 NORMAL AND RELIABLE TERMINATION OF THE |
|
|
62 C ROUTINE. IT IS ASSUMED THAT THE REQUESTED |
|
|
63 C ACCURACY HAS BEEN ACHIEVED. |
|
|
64 C IER.GT.0 ABNORMAL TERMINATION OF THE ROUTINE. THE |
|
|
65 C ESTIMATES FOR RESULT AND ERROR ARE LESS |
|
|
66 C RELIABLE. IT IS ASSUMED THAT THE REQUESTED |
|
|
67 C ACCURACY HAS NOT BEEN ACHIEVED. |
|
|
68 C IER.LT.0 EXIT REQUESTED FROM USER-SUPPLIED |
|
|
69 C FUNCTION. |
|
|
70 C |
|
|
71 C ERROR MESSAGES |
|
|
72 C IER = 1 MAXIMUM NUMBER OF SUBDIVISIONS ALLOWED |
|
|
73 C HAS BEEN ACHIEVED. ONE CAN ALLOW MORE |
|
|
74 C SUBDIVISIONS BY INCREASING THE VALUE OF |
|
|
75 C LIMIT (AND TAKING THE ACCORDING DIMENSION |
|
|
76 C ADJUSTMENTS INTO ACCOUNT). HOWEVER,IF |
|
|
77 C THIS YIELDS NO IMPROVEMENT IT IS ADVISED |
|
|
78 C TO ANALYZE THE INTEGRAND IN ORDER TO |
|
|
79 C DETERMINE THE INTEGRATION DIFFICULTIES. |
|
|
80 C IF THE POSITION OF A LOCAL DIFFICULTY CAN |
|
|
81 C BE DETERMINED (E.G. SINGULARITY, |
|
|
82 C DISCONTINUITY WITHIN THE INTERVAL) ONE |
|
|
83 C WILL PROBABLY GAIN FROM SPLITTING UP THE |
|
|
84 C INTERVAL AT THIS POINT AND CALLING THE |
|
|
85 C INTEGRATOR ON THE SUBRANGES. IF POSSIBLE, |
|
|
86 C AN APPROPRIATE SPECIAL-PURPOSE INTEGRATOR |
|
|
87 C SHOULD BE USED, WHICH IS DESIGNED FOR |
|
|
88 C HANDLING THE TYPE OF DIFFICULTY INVOLVED. |
|
|
89 C = 2 THE OCCURRENCE OF ROUNDOFF ERROR IS |
|
|
90 C DETECTED, WHICH PREVENTS THE REQUESTED |
|
|
91 C TOLERANCE FROM BEING ACHIEVED. |
|
|
92 C THE ERROR MAY BE UNDER-ESTIMATED. |
|
|
93 C = 3 EXTREMELY BAD INTEGRAND BEHAVIOUR OCCURS |
|
|
94 C AT SOME POINTS OF THE INTEGRATION |
|
|
95 C INTERVAL. |
|
|
96 C = 4 THE ALGORITHM DOES NOT CONVERGE. |
|
|
97 C ROUNDOFF ERROR IS DETECTED IN THE |
|
|
98 C EXTRAPOLATION TABLE. |
|
|
99 C IT IS ASSUMED THAT THE REQUESTED TOLERANCE |
|
|
100 C CANNOT BE ACHIEVED, AND THAT THE RETURNED |
|
|
101 C RESULT IS THE BEST WHICH CAN BE OBTAINED. |
|
|
102 C = 5 THE INTEGRAL IS PROBABLY DIVERGENT, OR |
|
|
103 C SLOWLY CONVERGENT. IT MUST BE NOTED THAT |
|
|
104 C DIVERGENCE CAN OCCUR WITH ANY OTHER VALUE |
|
|
105 C OF IER. |
|
|
106 C = 6 THE INPUT IS INVALID, BECAUSE |
|
|
107 C (EPSABS.LE.0 AND |
|
|
108 C EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28), |
|
|
109 C RESULT, ABSERR, NEVAL, LAST, RLIST(1), |
|
|
110 C ELIST(1) AND IORD(1) ARE SET TO ZERO. |
|
|
111 C ALIST(1) AND BLIST(1) ARE SET TO 0 |
|
|
112 C AND 1 RESPECTIVELY. |
|
|
113 C |
|
|
114 C ALIST - DOUBLE PRECISION |
|
|
115 C VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST |
|
|
116 C LAST ELEMENTS OF WHICH ARE THE LEFT |
|
|
117 C END POINTS OF THE SUBINTERVALS IN THE PARTITION |
|
|
118 C OF THE TRANSFORMED INTEGRATION RANGE (0,1). |
|
|
119 C |
|
|
120 C BLIST - DOUBLE PRECISION |
|
|
121 C VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST |
|
|
122 C LAST ELEMENTS OF WHICH ARE THE RIGHT |
|
|
123 C END POINTS OF THE SUBINTERVALS IN THE PARTITION |
|
|
124 C OF THE TRANSFORMED INTEGRATION RANGE (0,1). |
|
|
125 C |
|
|
126 C RLIST - DOUBLE PRECISION |
|
|
127 C VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST |
|
|
128 C LAST ELEMENTS OF WHICH ARE THE INTEGRAL |
|
|
129 C APPROXIMATIONS ON THE SUBINTERVALS |
|
|
130 C |
|
|
131 C ELIST - DOUBLE PRECISION |
|
|
132 C VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST |
|
|
133 C LAST ELEMENTS OF WHICH ARE THE MODULI OF THE |
|
|
134 C ABSOLUTE ERROR ESTIMATES ON THE SUBINTERVALS |
|
|
135 C |
|
|
136 C IORD - INTEGER |
|
|
137 C VECTOR OF DIMENSION LIMIT, THE FIRST K |
|
|
138 C ELEMENTS OF WHICH ARE POINTERS TO THE |
|
|
139 C ERROR ESTIMATES OVER THE SUBINTERVALS, |
|
|
140 C SUCH THAT ELIST(IORD(1)), ..., ELIST(IORD(K)) |
|
|
141 C FORM A DECREASING SEQUENCE, WITH K = LAST |
|
|
142 C IF LAST.LE.(LIMIT/2+2), AND K = LIMIT+1-LAST |
|
|
143 C OTHERWISE |
|
|
144 C |
|
|
145 C LAST - INTEGER |
|
|
146 C NUMBER OF SUBINTERVALS ACTUALLY PRODUCED |
|
|
147 C IN THE SUBDIVISION PROCESS |
|
|
148 C |
|
|
149 C***REFERENCES (NONE) |
|
|
150 C***ROUTINES CALLED D1MACH,DQELG,DQK15I,DQPSRT |
|
|
151 C***END PROLOGUE DQAGIE |
|
|
152 DOUBLE PRECISION ABSEPS,ABSERR,ALIST,AREA,AREA1,AREA12,AREA2,A1, |
|
|
153 * A2,BLIST,BOUN,BOUND,B1,B2,CORREC,DABS,DEFABS,DEFAB1,DEFAB2, |
|
|
154 * DMAX1,DRES,D1MACH,ELIST,EPMACH,EPSABS,EPSREL,ERLARG,ERLAST, |
|
3135
|
155 * ERRBND,ERRMAX,ERROR1,ERROR2,ERRO12,ERRSUM,ERTEST,OFLOW,RESABS, |
|
2329
|
156 * RESEPS,RESULT,RES3LA,RLIST,RLIST2,SMALL,UFLOW |
|
|
157 INTEGER ID,IER,IERRO,INF,IORD,IROFF1,IROFF2,IROFF3,JUPBND,K,KSGN, |
|
|
158 * KTMIN,LAST,LIMIT,MAXERR,NEVAL,NRES,NRMAX,NUMRL2 |
|
|
159 LOGICAL EXTRAP,NOEXT |
|
|
160 C |
|
|
161 DIMENSION ALIST(LIMIT),BLIST(LIMIT),ELIST(LIMIT),IORD(LIMIT), |
|
|
162 * RES3LA(3),RLIST(LIMIT),RLIST2(52) |
|
|
163 C |
|
|
164 EXTERNAL F |
|
|
165 C |
|
|
166 C THE DIMENSION OF RLIST2 IS DETERMINED BY THE VALUE OF |
|
|
167 C LIMEXP IN SUBROUTINE DQELG. |
|
|
168 C |
|
|
169 C |
|
|
170 C LIST OF MAJOR VARIABLES |
|
|
171 C ----------------------- |
|
|
172 C |
|
|
173 C ALIST - LIST OF LEFT END POINTS OF ALL SUBINTERVALS |
|
|
174 C CONSIDERED UP TO NOW |
|
|
175 C BLIST - LIST OF RIGHT END POINTS OF ALL SUBINTERVALS |
|
|
176 C CONSIDERED UP TO NOW |
|
|
177 C RLIST(I) - APPROXIMATION TO THE INTEGRAL OVER |
|
|
178 C (ALIST(I),BLIST(I)) |
|
|
179 C RLIST2 - ARRAY OF DIMENSION AT LEAST (LIMEXP+2), |
|
|
180 C CONTAINING THE PART OF THE EPSILON TABLE |
|
|
181 C WICH IS STILL NEEDED FOR FURTHER COMPUTATIONS |
|
|
182 C ELIST(I) - ERROR ESTIMATE APPLYING TO RLIST(I) |
|
|
183 C MAXERR - POINTER TO THE INTERVAL WITH LARGEST ERROR |
|
|
184 C ESTIMATE |
|
|
185 C ERRMAX - ELIST(MAXERR) |
|
|
186 C ERLAST - ERROR ON THE INTERVAL CURRENTLY SUBDIVIDED |
|
|
187 C (BEFORE THAT SUBDIVISION HAS TAKEN PLACE) |
|
|
188 C AREA - SUM OF THE INTEGRALS OVER THE SUBINTERVALS |
|
|
189 C ERRSUM - SUM OF THE ERRORS OVER THE SUBINTERVALS |
|
|
190 C ERRBND - REQUESTED ACCURACY MAX(EPSABS,EPSREL* |
|
|
191 C ABS(RESULT)) |
|
|
192 C *****1 - VARIABLE FOR THE LEFT SUBINTERVAL |
|
|
193 C *****2 - VARIABLE FOR THE RIGHT SUBINTERVAL |
|
|
194 C LAST - INDEX FOR SUBDIVISION |
|
|
195 C NRES - NUMBER OF CALLS TO THE EXTRAPOLATION ROUTINE |
|
|
196 C NUMRL2 - NUMBER OF ELEMENTS CURRENTLY IN RLIST2. IF AN |
|
|
197 C APPROPRIATE APPROXIMATION TO THE COMPOUNDED |
|
|
198 C INTEGRAL HAS BEEN OBTAINED, IT IS PUT IN |
|
|
199 C RLIST2(NUMRL2) AFTER NUMRL2 HAS BEEN INCREASED |
|
|
200 C BY ONE. |
|
|
201 C SMALL - LENGTH OF THE SMALLEST INTERVAL CONSIDERED UP |
|
|
202 C TO NOW, MULTIPLIED BY 1.5 |
|
|
203 C ERLARG - SUM OF THE ERRORS OVER THE INTERVALS LARGER |
|
|
204 C THAN THE SMALLEST INTERVAL CONSIDERED UP TO NOW |
|
|
205 C EXTRAP - LOGICAL VARIABLE DENOTING THAT THE ROUTINE |
|
|
206 C IS ATTEMPTING TO PERFORM EXTRAPOLATION. I.E. |
|
|
207 C BEFORE SUBDIVIDING THE SMALLEST INTERVAL WE |
|
|
208 C TRY TO DECREASE THE VALUE OF ERLARG. |
|
|
209 C NOEXT - LOGICAL VARIABLE DENOTING THAT EXTRAPOLATION |
|
|
210 C IS NO LONGER ALLOWED (TRUE-VALUE) |
|
|
211 C |
|
|
212 C MACHINE DEPENDENT CONSTANTS |
|
|
213 C --------------------------- |
|
|
214 C |
|
|
215 C EPMACH IS THE LARGEST RELATIVE SPACING. |
|
|
216 C UFLOW IS THE SMALLEST POSITIVE MAGNITUDE. |
|
|
217 C OFLOW IS THE LARGEST POSITIVE MAGNITUDE. |
|
|
218 C |
|
|
219 C***FIRST EXECUTABLE STATEMENT DQAGIE |
|
|
220 EPMACH = D1MACH(4) |
|
|
221 C |
|
|
222 C TEST ON VALIDITY OF PARAMETERS |
|
|
223 C ----------------------------- |
|
|
224 C |
|
|
225 IER = 0 |
|
|
226 NEVAL = 0 |
|
|
227 LAST = 0 |
|
|
228 RESULT = 0.0D+00 |
|
|
229 ABSERR = 0.0D+00 |
|
|
230 ALIST(1) = 0.0D+00 |
|
|
231 BLIST(1) = 0.1D+01 |
|
|
232 RLIST(1) = 0.0D+00 |
|
|
233 ELIST(1) = 0.0D+00 |
|
|
234 IORD(1) = 0 |
|
|
235 IF(EPSABS.LE.0.0D+00.AND.EPSREL.LT.DMAX1(0.5D+02*EPMACH,0.5D-28)) |
|
|
236 * IER = 6 |
|
|
237 IF(IER.EQ.6) GO TO 999 |
|
|
238 C |
|
|
239 C |
|
|
240 C FIRST APPROXIMATION TO THE INTEGRAL |
|
|
241 C ----------------------------------- |
|
|
242 C |
|
|
243 C DETERMINE THE INTERVAL TO BE MAPPED ONTO (0,1). |
|
|
244 C IF INF = 2 THE INTEGRAL IS COMPUTED AS I = I1+I2, WHERE |
|
|
245 C I1 = INTEGRAL OF F OVER (-INFINITY,0), |
|
|
246 C I2 = INTEGRAL OF F OVER (0,+INFINITY). |
|
|
247 C |
|
|
248 BOUN = BOUND |
|
|
249 IF(INF.EQ.2) BOUN = 0.0D+00 |
|
|
250 CALL DQK15I(F,BOUN,INF,0.0D+00,0.1D+01,RESULT,ABSERR, |
|
|
251 * DEFABS,RESABS,IER) |
|
|
252 IF (IER .LT. 0) RETURN |
|
|
253 C |
|
|
254 C TEST ON ACCURACY |
|
|
255 C |
|
|
256 LAST = 1 |
|
|
257 RLIST(1) = RESULT |
|
|
258 ELIST(1) = ABSERR |
|
|
259 IORD(1) = 1 |
|
|
260 DRES = DABS(RESULT) |
|
|
261 ERRBND = DMAX1(EPSABS,EPSREL*DRES) |
|
|
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 130 |
|
|
266 C |
|
|
267 C INITIALIZATION |
|
|
268 C -------------- |
|
|
269 C |
|
|
270 UFLOW = D1MACH(1) |
|
|
271 OFLOW = D1MACH(2) |
|
|
272 RLIST2(1) = RESULT |
|
|
273 ERRMAX = ABSERR |
|
|
274 MAXERR = 1 |
|
|
275 AREA = RESULT |
|
|
276 ERRSUM = ABSERR |
|
|
277 ABSERR = OFLOW |
|
|
278 NRMAX = 1 |
|
|
279 NRES = 0 |
|
|
280 KTMIN = 0 |
|
|
281 NUMRL2 = 2 |
|
|
282 EXTRAP = .FALSE. |
|
|
283 NOEXT = .FALSE. |
|
|
284 IERRO = 0 |
|
|
285 IROFF1 = 0 |
|
|
286 IROFF2 = 0 |
|
|
287 IROFF3 = 0 |
|
|
288 KSGN = -1 |
|
|
289 IF(DRES.GE.(0.1D+01-0.5D+02*EPMACH)*DEFABS) KSGN = 1 |
|
|
290 C |
|
|
291 C MAIN DO-LOOP |
|
|
292 C ------------ |
|
|
293 C |
|
|
294 DO 90 LAST = 2,LIMIT |
|
|
295 C |
|
|
296 C BISECT THE SUBINTERVAL WITH NRMAX-TH LARGEST ERROR ESTIMATE. |
|
|
297 C |
|
|
298 A1 = ALIST(MAXERR) |
|
|
299 B1 = 0.5D+00*(ALIST(MAXERR)+BLIST(MAXERR)) |
|
|
300 A2 = B1 |
|
|
301 B2 = BLIST(MAXERR) |
|
|
302 ERLAST = ERRMAX |
|
|
303 CALL DQK15I(F,BOUN,INF,A1,B1,AREA1,ERROR1,RESABS,DEFAB1,IER) |
|
|
304 IF (IER .LT. 0) RETURN |
|
|
305 CALL DQK15I(F,BOUN,INF,A2,B2,AREA2,ERROR2,RESABS,DEFAB2,IER) |
|
|
306 IF (IER .LT. 0) RETURN |
|
|
307 C |
|
|
308 C IMPROVE PREVIOUS APPROXIMATIONS TO INTEGRAL |
|
|
309 C AND ERROR AND TEST FOR ACCURACY. |
|
|
310 C |
|
|
311 AREA12 = AREA1+AREA2 |
|
|
312 ERRO12 = ERROR1+ERROR2 |
|
|
313 ERRSUM = ERRSUM+ERRO12-ERRMAX |
|
|
314 AREA = AREA+AREA12-RLIST(MAXERR) |
|
|
315 IF(DEFAB1.EQ.ERROR1.OR.DEFAB2.EQ.ERROR2)GO TO 15 |
|
|
316 IF(DABS(RLIST(MAXERR)-AREA12).GT.0.1D-04*DABS(AREA12) |
|
|
317 * .OR.ERRO12.LT.0.99D+00*ERRMAX) GO TO 10 |
|
|
318 IF(EXTRAP) IROFF2 = IROFF2+1 |
|
|
319 IF(.NOT.EXTRAP) IROFF1 = IROFF1+1 |
|
|
320 10 IF(LAST.GT.10.AND.ERRO12.GT.ERRMAX) IROFF3 = IROFF3+1 |
|
|
321 15 RLIST(MAXERR) = AREA1 |
|
|
322 RLIST(LAST) = AREA2 |
|
|
323 ERRBND = DMAX1(EPSABS,EPSREL*DABS(AREA)) |
|
|
324 C |
|
|
325 C TEST FOR ROUNDOFF ERROR AND EVENTUALLY SET ERROR FLAG. |
|
|
326 C |
|
|
327 IF(IROFF1+IROFF2.GE.10.OR.IROFF3.GE.20) IER = 2 |
|
|
328 IF(IROFF2.GE.5) IERRO = 3 |
|
|
329 C |
|
|
330 C SET ERROR FLAG IN THE CASE THAT THE NUMBER OF |
|
|
331 C SUBINTERVALS EQUALS LIMIT. |
|
|
332 C |
|
|
333 IF(LAST.EQ.LIMIT) IER = 1 |
|
|
334 C |
|
|
335 C SET ERROR FLAG IN THE CASE OF BAD INTEGRAND BEHAVIOUR |
|
|
336 C AT SOME POINTS OF THE INTEGRATION RANGE. |
|
|
337 C |
|
|
338 IF(DMAX1(DABS(A1),DABS(B2)).LE.(0.1D+01+0.1D+03*EPMACH)* |
|
|
339 * (DABS(A2)+0.1D+04*UFLOW)) IER = 4 |
|
|
340 C |
|
|
341 C APPEND THE NEWLY-CREATED INTERVALS TO THE LIST. |
|
|
342 C |
|
|
343 IF(ERROR2.GT.ERROR1) GO TO 20 |
|
|
344 ALIST(LAST) = A2 |
|
|
345 BLIST(MAXERR) = B1 |
|
|
346 BLIST(LAST) = B2 |
|
|
347 ELIST(MAXERR) = ERROR1 |
|
|
348 ELIST(LAST) = ERROR2 |
|
|
349 GO TO 30 |
|
|
350 20 ALIST(MAXERR) = A2 |
|
|
351 ALIST(LAST) = A1 |
|
|
352 BLIST(LAST) = B1 |
|
|
353 RLIST(MAXERR) = AREA2 |
|
|
354 RLIST(LAST) = AREA1 |
|
|
355 ELIST(MAXERR) = ERROR2 |
|
|
356 ELIST(LAST) = ERROR1 |
|
|
357 C |
|
|
358 C CALL SUBROUTINE DQPSRT TO MAINTAIN THE DESCENDING ORDERING |
|
|
359 C IN THE LIST OF ERROR ESTIMATES AND SELECT THE SUBINTERVAL |
|
|
360 C WITH NRMAX-TH LARGEST ERROR ESTIMATE (TO BE BISECTED NEXT). |
|
|
361 C |
|
|
362 30 CALL DQPSRT(LIMIT,LAST,MAXERR,ERRMAX,ELIST,IORD,NRMAX) |
|
|
363 IF(ERRSUM.LE.ERRBND) GO TO 115 |
|
|
364 IF(IER.NE.0) GO TO 100 |
|
|
365 IF(LAST.EQ.2) GO TO 80 |
|
|
366 IF(NOEXT) GO TO 90 |
|
|
367 ERLARG = ERLARG-ERLAST |
|
|
368 IF(DABS(B1-A1).GT.SMALL) ERLARG = ERLARG+ERRO12 |
|
|
369 IF(EXTRAP) GO TO 40 |
|
|
370 C |
|
|
371 C TEST WHETHER THE INTERVAL TO BE BISECTED NEXT IS THE |
|
|
372 C SMALLEST INTERVAL. |
|
|
373 C |
|
|
374 IF(DABS(BLIST(MAXERR)-ALIST(MAXERR)).GT.SMALL) GO TO 90 |
|
|
375 EXTRAP = .TRUE. |
|
|
376 NRMAX = 2 |
|
|
377 40 IF(IERRO.EQ.3.OR.ERLARG.LE.ERTEST) GO TO 60 |
|
|
378 C |
|
|
379 C THE SMALLEST INTERVAL HAS THE LARGEST ERROR. |
|
|
380 C BEFORE BISECTING DECREASE THE SUM OF THE ERRORS OVER THE |
|
|
381 C LARGER INTERVALS (ERLARG) AND PERFORM EXTRAPOLATION. |
|
|
382 C |
|
|
383 ID = NRMAX |
|
|
384 JUPBND = LAST |
|
|
385 IF(LAST.GT.(2+LIMIT/2)) JUPBND = LIMIT+3-LAST |
|
|
386 DO 50 K = ID,JUPBND |
|
|
387 MAXERR = IORD(NRMAX) |
|
|
388 ERRMAX = ELIST(MAXERR) |
|
|
389 IF(DABS(BLIST(MAXERR)-ALIST(MAXERR)).GT.SMALL) GO TO 90 |
|
|
390 NRMAX = NRMAX+1 |
|
|
391 50 CONTINUE |
|
|
392 C |
|
|
393 C PERFORM EXTRAPOLATION. |
|
|
394 C |
|
|
395 60 NUMRL2 = NUMRL2+1 |
|
|
396 RLIST2(NUMRL2) = AREA |
|
|
397 CALL DQELG(NUMRL2,RLIST2,RESEPS,ABSEPS,RES3LA,NRES) |
|
|
398 KTMIN = KTMIN+1 |
|
|
399 IF(KTMIN.GT.5.AND.ABSERR.LT.0.1D-02*ERRSUM) IER = 5 |
|
|
400 IF(ABSEPS.GE.ABSERR) GO TO 70 |
|
|
401 KTMIN = 0 |
|
|
402 ABSERR = ABSEPS |
|
|
403 RESULT = RESEPS |
|
|
404 CORREC = ERLARG |
|
|
405 ERTEST = DMAX1(EPSABS,EPSREL*DABS(RESEPS)) |
|
|
406 IF(ABSERR.LE.ERTEST) GO TO 100 |
|
|
407 C |
|
|
408 C PREPARE BISECTION OF THE SMALLEST INTERVAL. |
|
|
409 C |
|
|
410 70 IF(NUMRL2.EQ.1) NOEXT = .TRUE. |
|
|
411 IF(IER.EQ.5) GO TO 100 |
|
|
412 MAXERR = IORD(1) |
|
|
413 ERRMAX = ELIST(MAXERR) |
|
|
414 NRMAX = 1 |
|
|
415 EXTRAP = .FALSE. |
|
|
416 SMALL = SMALL*0.5D+00 |
|
|
417 ERLARG = ERRSUM |
|
|
418 GO TO 90 |
|
|
419 80 SMALL = 0.375D+00 |
|
|
420 ERLARG = ERRSUM |
|
|
421 ERTEST = ERRBND |
|
|
422 RLIST2(2) = AREA |
|
|
423 90 CONTINUE |
|
|
424 C |
|
|
425 C SET FINAL RESULT AND ERROR ESTIMATE. |
|
|
426 C ------------------------------------ |
|
|
427 C |
|
|
428 100 IF(ABSERR.EQ.OFLOW) GO TO 115 |
|
|
429 IF((IER+IERRO).EQ.0) GO TO 110 |
|
|
430 IF(IERRO.EQ.3) ABSERR = ABSERR+CORREC |
|
|
431 IF(IER.EQ.0) IER = 3 |
|
|
432 IF(RESULT.NE.0.0D+00.AND.AREA.NE.0.0D+00)GO TO 105 |
|
|
433 IF(ABSERR.GT.ERRSUM)GO TO 115 |
|
|
434 IF(AREA.EQ.0.0D+00) GO TO 130 |
|
|
435 GO TO 110 |
|
|
436 105 IF(ABSERR/DABS(RESULT).GT.ERRSUM/DABS(AREA))GO TO 115 |
|
|
437 C |
|
|
438 C TEST ON DIVERGENCE |
|
|
439 C |
|
|
440 110 IF(KSGN.EQ.(-1).AND.DMAX1(DABS(RESULT),DABS(AREA)).LE. |
|
|
441 * DEFABS*0.1D-01) GO TO 130 |
|
|
442 IF(0.1D-01.GT.(RESULT/AREA).OR.(RESULT/AREA).GT.0.1D+03. |
|
|
443 *OR.ERRSUM.GT.DABS(AREA)) IER = 6 |
|
|
444 GO TO 130 |
|
|
445 C |
|
|
446 C COMPUTE GLOBAL INTEGRAL SUM. |
|
|
447 C |
|
|
448 115 RESULT = 0.0D+00 |
|
|
449 DO 120 K = 1,LAST |
|
|
450 RESULT = RESULT+RLIST(K) |
|
|
451 120 CONTINUE |
|
|
452 ABSERR = ERRSUM |
|
|
453 130 NEVAL = 30*LAST-15 |
|
|
454 IF(INF.EQ.2) NEVAL = 2*NEVAL |
|
|
455 IF(IER.GT.2) IER=IER-1 |
|
|
456 999 RETURN |
|
|
457 END |
