00001 SUBROUTINE DQAGIE(F,BOUND,INF,EPSABS,EPSREL,LIMIT,RESULT,ABSERR, 00002 * NEVAL,IER,ALIST,BLIST,RLIST,ELIST,IORD,LAST) 00003 C***BEGIN PROLOGUE DQAGIE 00004 C***DATE WRITTEN 800101 (YYMMDD) 00005 C***REVISION DATE 830518 (YYMMDD) 00006 C***CATEGORY NO. H2A3A1,H2A4A1 00007 C***KEYWORDS AUTOMATIC INTEGRATOR, INFINITE INTERVALS, 00008 C GENERAL-PURPOSE, TRANSFORMATION, EXTRAPOLATION, 00009 C GLOBALLY ADAPTIVE 00010 C***AUTHOR PIESSENS,ROBERT,APPL. MATH & PROGR. DIV - K.U.LEUVEN 00011 C DE DONCKER,ELISE,APPL. MATH & PROGR. DIV - K.U.LEUVEN 00012 C***PURPOSE THE ROUTINE CALCULATES AN APPROXIMATION RESULT TO A GIVEN 00013 C INTEGRAL I = INTEGRAL OF F OVER (BOUND,+INFINITY) 00014 C OR I = INTEGRAL OF F OVER (-INFINITY,BOUND) 00015 C OR I = INTEGRAL OF F OVER (-INFINITY,+INFINITY), 00016 C HOPEFULLY SATISFYING FOLLOWING CLAIM FOR ACCURACY 00017 C ABS(I-RESULT).LE.MAX(EPSABS,EPSREL*ABS(I)) 00018 C***DESCRIPTION 00019 C 00020 C INTEGRATION OVER INFINITE INTERVALS 00021 C STANDARD FORTRAN SUBROUTINE 00022 C 00023 C F - SUBROUTINE F(X,IERR,RESULT) DEFINING THE INTEGRAND 00024 C FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE 00025 C DECLARED E X T E R N A L IN THE DRIVER PROGRAM. 00026 C 00027 C BOUND - DOUBLE PRECISION 00028 C FINITE BOUND OF INTEGRATION RANGE 00029 C (HAS NO MEANING IF INTERVAL IS DOUBLY-INFINITE) 00030 C 00031 C INF - DOUBLE PRECISION 00032 C INDICATING THE KIND OF INTEGRATION RANGE INVOLVED 00033 C INF = 1 CORRESPONDS TO (BOUND,+INFINITY), 00034 C INF = -1 TO (-INFINITY,BOUND), 00035 C INF = 2 TO (-INFINITY,+INFINITY). 00036 C 00037 C EPSABS - DOUBLE PRECISION 00038 C ABSOLUTE ACCURACY REQUESTED 00039 C EPSREL - DOUBLE PRECISION 00040 C RELATIVE ACCURACY REQUESTED 00041 C IF EPSABS.LE.0 00042 C AND EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28), 00043 C THE ROUTINE WILL END WITH IER = 6. 00044 C 00045 C LIMIT - INTEGER 00046 C GIVES AN UPPER BOUND ON THE NUMBER OF SUBINTERVALS 00047 C IN THE PARTITION OF (A,B), LIMIT.GE.1 00048 C 00049 C ON RETURN 00050 C RESULT - DOUBLE PRECISION 00051 C APPROXIMATION TO THE INTEGRAL 00052 C 00053 C ABSERR - DOUBLE PRECISION 00054 C ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR, 00055 C WHICH SHOULD EQUAL OR EXCEED ABS(I-RESULT) 00056 C 00057 C NEVAL - INTEGER 00058 C NUMBER OF INTEGRAND EVALUATIONS 00059 C 00060 C IER - INTEGER 00061 C IER = 0 NORMAL AND RELIABLE TERMINATION OF THE 00062 C ROUTINE. IT IS ASSUMED THAT THE REQUESTED 00063 C ACCURACY HAS BEEN ACHIEVED. 00064 C IER.GT.0 ABNORMAL TERMINATION OF THE ROUTINE. THE 00065 C ESTIMATES FOR RESULT AND ERROR ARE LESS 00066 C RELIABLE. IT IS ASSUMED THAT THE REQUESTED 00067 C ACCURACY HAS NOT BEEN ACHIEVED. 00068 C IER.LT.0 EXIT REQUESTED FROM USER-SUPPLIED 00069 C FUNCTION. 00070 C 00071 C ERROR MESSAGES 00072 C IER = 1 MAXIMUM NUMBER OF SUBDIVISIONS ALLOWED 00073 C HAS BEEN ACHIEVED. ONE CAN ALLOW MORE 00074 C SUBDIVISIONS BY INCREASING THE VALUE OF 00075 C LIMIT (AND TAKING THE ACCORDING DIMENSION 00076 C ADJUSTMENTS INTO ACCOUNT). HOWEVER,IF 00077 C THIS YIELDS NO IMPROVEMENT IT IS ADVISED 00078 C TO ANALYZE THE INTEGRAND IN ORDER TO 00079 C DETERMINE THE INTEGRATION DIFFICULTIES. 00080 C IF THE POSITION OF A LOCAL DIFFICULTY CAN 00081 C BE DETERMINED (E.G. SINGULARITY, 00082 C DISCONTINUITY WITHIN THE INTERVAL) ONE 00083 C WILL PROBABLY GAIN FROM SPLITTING UP THE 00084 C INTERVAL AT THIS POINT AND CALLING THE 00085 C INTEGRATOR ON THE SUBRANGES. IF POSSIBLE, 00086 C AN APPROPRIATE SPECIAL-PURPOSE INTEGRATOR 00087 C SHOULD BE USED, WHICH IS DESIGNED FOR 00088 C HANDLING THE TYPE OF DIFFICULTY INVOLVED. 00089 C = 2 THE OCCURRENCE OF ROUNDOFF ERROR IS 00090 C DETECTED, WHICH PREVENTS THE REQUESTED 00091 C TOLERANCE FROM BEING ACHIEVED. 00092 C THE ERROR MAY BE UNDER-ESTIMATED. 00093 C = 3 EXTREMELY BAD INTEGRAND BEHAVIOUR OCCURS 00094 C AT SOME POINTS OF THE INTEGRATION 00095 C INTERVAL. 00096 C = 4 THE ALGORITHM DOES NOT CONVERGE. 00097 C ROUNDOFF ERROR IS DETECTED IN THE 00098 C EXTRAPOLATION TABLE. 00099 C IT IS ASSUMED THAT THE REQUESTED TOLERANCE 00100 C CANNOT BE ACHIEVED, AND THAT THE RETURNED 00101 C RESULT IS THE BEST WHICH CAN BE OBTAINED. 00102 C = 5 THE INTEGRAL IS PROBABLY DIVERGENT, OR 00103 C SLOWLY CONVERGENT. IT MUST BE NOTED THAT 00104 C DIVERGENCE CAN OCCUR WITH ANY OTHER VALUE 00105 C OF IER. 00106 C = 6 THE INPUT IS INVALID, BECAUSE 00107 C (EPSABS.LE.0 AND 00108 C EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28), 00109 C RESULT, ABSERR, NEVAL, LAST, RLIST(1), 00110 C ELIST(1) AND IORD(1) ARE SET TO ZERO. 00111 C ALIST(1) AND BLIST(1) ARE SET TO 0 00112 C AND 1 RESPECTIVELY. 00113 C 00114 C ALIST - DOUBLE PRECISION 00115 C VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST 00116 C LAST ELEMENTS OF WHICH ARE THE LEFT 00117 C END POINTS OF THE SUBINTERVALS IN THE PARTITION 00118 C OF THE TRANSFORMED INTEGRATION RANGE (0,1). 00119 C 00120 C BLIST - DOUBLE PRECISION 00121 C VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST 00122 C LAST ELEMENTS OF WHICH ARE THE RIGHT 00123 C END POINTS OF THE SUBINTERVALS IN THE PARTITION 00124 C OF THE TRANSFORMED INTEGRATION RANGE (0,1). 00125 C 00126 C RLIST - DOUBLE PRECISION 00127 C VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST 00128 C LAST ELEMENTS OF WHICH ARE THE INTEGRAL 00129 C APPROXIMATIONS ON THE SUBINTERVALS 00130 C 00131 C ELIST - DOUBLE PRECISION 00132 C VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST 00133 C LAST ELEMENTS OF WHICH ARE THE MODULI OF THE 00134 C ABSOLUTE ERROR ESTIMATES ON THE SUBINTERVALS 00135 C 00136 C IORD - INTEGER 00137 C VECTOR OF DIMENSION LIMIT, THE FIRST K 00138 C ELEMENTS OF WHICH ARE POINTERS TO THE 00139 C ERROR ESTIMATES OVER THE SUBINTERVALS, 00140 C SUCH THAT ELIST(IORD(1)), ..., ELIST(IORD(K)) 00141 C FORM A DECREASING SEQUENCE, WITH K = LAST 00142 C IF LAST.LE.(LIMIT/2+2), AND K = LIMIT+1-LAST 00143 C OTHERWISE 00144 C 00145 C LAST - INTEGER 00146 C NUMBER OF SUBINTERVALS ACTUALLY PRODUCED 00147 C IN THE SUBDIVISION PROCESS 00148 C 00149 C***REFERENCES (NONE) 00150 C***ROUTINES CALLED D1MACH,DQELG,DQK15I,DQPSRT 00151 C***END PROLOGUE DQAGIE 00152 DOUBLE PRECISION ABSEPS,ABSERR,ALIST,AREA,AREA1,AREA12,AREA2,A1, 00153 * A2,BLIST,BOUN,BOUND,B1,B2,CORREC,DABS,DEFABS,DEFAB1,DEFAB2, 00154 * DMAX1,DRES,D1MACH,ELIST,EPMACH,EPSABS,EPSREL,ERLARG,ERLAST, 00155 * ERRBND,ERRMAX,ERROR1,ERROR2,ERRO12,ERRSUM,ERTEST,OFLOW,RESABS, 00156 * RESEPS,RESULT,RES3LA,RLIST,RLIST2,SMALL,UFLOW 00157 INTEGER ID,IER,IERRO,INF,IORD,IROFF1,IROFF2,IROFF3,JUPBND,K,KSGN, 00158 * KTMIN,LAST,LIMIT,MAXERR,NEVAL,NRES,NRMAX,NUMRL2 00159 LOGICAL EXTRAP,NOEXT 00160 C 00161 DIMENSION ALIST(LIMIT),BLIST(LIMIT),ELIST(LIMIT),IORD(LIMIT), 00162 * RES3LA(3),RLIST(LIMIT),RLIST2(52) 00163 C 00164 EXTERNAL F 00165 C 00166 C THE DIMENSION OF RLIST2 IS DETERMINED BY THE VALUE OF 00167 C LIMEXP IN SUBROUTINE DQELG. 00168 C 00169 C 00170 C LIST OF MAJOR VARIABLES 00171 C ----------------------- 00172 C 00173 C ALIST - LIST OF LEFT END POINTS OF ALL SUBINTERVALS 00174 C CONSIDERED UP TO NOW 00175 C BLIST - LIST OF RIGHT END POINTS OF ALL SUBINTERVALS 00176 C CONSIDERED UP TO NOW 00177 C RLIST(I) - APPROXIMATION TO THE INTEGRAL OVER 00178 C (ALIST(I),BLIST(I)) 00179 C RLIST2 - ARRAY OF DIMENSION AT LEAST (LIMEXP+2), 00180 C CONTAINING THE PART OF THE EPSILON TABLE 00181 C WICH IS STILL NEEDED FOR FURTHER COMPUTATIONS 00182 C ELIST(I) - ERROR ESTIMATE APPLYING TO RLIST(I) 00183 C MAXERR - POINTER TO THE INTERVAL WITH LARGEST ERROR 00184 C ESTIMATE 00185 C ERRMAX - ELIST(MAXERR) 00186 C ERLAST - ERROR ON THE INTERVAL CURRENTLY SUBDIVIDED 00187 C (BEFORE THAT SUBDIVISION HAS TAKEN PLACE) 00188 C AREA - SUM OF THE INTEGRALS OVER THE SUBINTERVALS 00189 C ERRSUM - SUM OF THE ERRORS OVER THE SUBINTERVALS 00190 C ERRBND - REQUESTED ACCURACY MAX(EPSABS,EPSREL* 00191 C ABS(RESULT)) 00192 C *****1 - VARIABLE FOR THE LEFT SUBINTERVAL 00193 C *****2 - VARIABLE FOR THE RIGHT SUBINTERVAL 00194 C LAST - INDEX FOR SUBDIVISION 00195 C NRES - NUMBER OF CALLS TO THE EXTRAPOLATION ROUTINE 00196 C NUMRL2 - NUMBER OF ELEMENTS CURRENTLY IN RLIST2. IF AN 00197 C APPROPRIATE APPROXIMATION TO THE COMPOUNDED 00198 C INTEGRAL HAS BEEN OBTAINED, IT IS PUT IN 00199 C RLIST2(NUMRL2) AFTER NUMRL2 HAS BEEN INCREASED 00200 C BY ONE. 00201 C SMALL - LENGTH OF THE SMALLEST INTERVAL CONSIDERED UP 00202 C TO NOW, MULTIPLIED BY 1.5 00203 C ERLARG - SUM OF THE ERRORS OVER THE INTERVALS LARGER 00204 C THAN THE SMALLEST INTERVAL CONSIDERED UP TO NOW 00205 C EXTRAP - LOGICAL VARIABLE DENOTING THAT THE ROUTINE 00206 C IS ATTEMPTING TO PERFORM EXTRAPOLATION. I.E. 00207 C BEFORE SUBDIVIDING THE SMALLEST INTERVAL WE 00208 C TRY TO DECREASE THE VALUE OF ERLARG. 00209 C NOEXT - LOGICAL VARIABLE DENOTING THAT EXTRAPOLATION 00210 C IS NO LONGER ALLOWED (TRUE-VALUE) 00211 C 00212 C MACHINE DEPENDENT CONSTANTS 00213 C --------------------------- 00214 C 00215 C EPMACH IS THE LARGEST RELATIVE SPACING. 00216 C UFLOW IS THE SMALLEST POSITIVE MAGNITUDE. 00217 C OFLOW IS THE LARGEST POSITIVE MAGNITUDE. 00218 C 00219 C***FIRST EXECUTABLE STATEMENT DQAGIE 00220 EPMACH = D1MACH(4) 00221 C 00222 C TEST ON VALIDITY OF PARAMETERS 00223 C ----------------------------- 00224 C 00225 IER = 0 00226 NEVAL = 0 00227 LAST = 0 00228 RESULT = 0.0D+00 00229 ABSERR = 0.0D+00 00230 ALIST(1) = 0.0D+00 00231 BLIST(1) = 0.1D+01 00232 RLIST(1) = 0.0D+00 00233 ELIST(1) = 0.0D+00 00234 IORD(1) = 0 00235 IF(EPSABS.LE.0.0D+00.AND.EPSREL.LT.DMAX1(0.5D+02*EPMACH,0.5D-28)) 00236 * IER = 6 00237 IF(IER.EQ.6) GO TO 999 00238 C 00239 C 00240 C FIRST APPROXIMATION TO THE INTEGRAL 00241 C ----------------------------------- 00242 C 00243 C DETERMINE THE INTERVAL TO BE MAPPED ONTO (0,1). 00244 C IF INF = 2 THE INTEGRAL IS COMPUTED AS I = I1+I2, WHERE 00245 C I1 = INTEGRAL OF F OVER (-INFINITY,0), 00246 C I2 = INTEGRAL OF F OVER (0,+INFINITY). 00247 C 00248 BOUN = BOUND 00249 IF(INF.EQ.2) BOUN = 0.0D+00 00250 CALL DQK15I(F,BOUN,INF,0.0D+00,0.1D+01,RESULT,ABSERR, 00251 * DEFABS,RESABS,IER) 00252 IF (IER .LT. 0) RETURN 00253 C 00254 C TEST ON ACCURACY 00255 C 00256 LAST = 1 00257 RLIST(1) = RESULT 00258 ELIST(1) = ABSERR 00259 IORD(1) = 1 00260 DRES = DABS(RESULT) 00261 ERRBND = DMAX1(EPSABS,EPSREL*DRES) 00262 IF(ABSERR.LE.1.0D+02*EPMACH*DEFABS.AND.ABSERR.GT.ERRBND) IER = 2 00263 IF(LIMIT.EQ.1) IER = 1 00264 IF(IER.NE.0.OR.(ABSERR.LE.ERRBND.AND.ABSERR.NE.RESABS).OR. 00265 * ABSERR.EQ.0.0D+00) GO TO 130 00266 C 00267 C INITIALIZATION 00268 C -------------- 00269 C 00270 UFLOW = D1MACH(1) 00271 OFLOW = D1MACH(2) 00272 RLIST2(1) = RESULT 00273 ERRMAX = ABSERR 00274 MAXERR = 1 00275 AREA = RESULT 00276 ERRSUM = ABSERR 00277 ABSERR = OFLOW 00278 NRMAX = 1 00279 NRES = 0 00280 KTMIN = 0 00281 NUMRL2 = 2 00282 EXTRAP = .FALSE. 00283 NOEXT = .FALSE. 00284 IERRO = 0 00285 IROFF1 = 0 00286 IROFF2 = 0 00287 IROFF3 = 0 00288 KSGN = -1 00289 IF(DRES.GE.(0.1D+01-0.5D+02*EPMACH)*DEFABS) KSGN = 1 00290 C 00291 C MAIN DO-LOOP 00292 C ------------ 00293 C 00294 DO 90 LAST = 2,LIMIT 00295 C 00296 C BISECT THE SUBINTERVAL WITH NRMAX-TH LARGEST ERROR ESTIMATE. 00297 C 00298 A1 = ALIST(MAXERR) 00299 B1 = 0.5D+00*(ALIST(MAXERR)+BLIST(MAXERR)) 00300 A2 = B1 00301 B2 = BLIST(MAXERR) 00302 ERLAST = ERRMAX 00303 CALL DQK15I(F,BOUN,INF,A1,B1,AREA1,ERROR1,RESABS,DEFAB1,IER) 00304 IF (IER .LT. 0) RETURN 00305 CALL DQK15I(F,BOUN,INF,A2,B2,AREA2,ERROR2,RESABS,DEFAB2,IER) 00306 IF (IER .LT. 0) RETURN 00307 C 00308 C IMPROVE PREVIOUS APPROXIMATIONS TO INTEGRAL 00309 C AND ERROR AND TEST FOR ACCURACY. 00310 C 00311 AREA12 = AREA1+AREA2 00312 ERRO12 = ERROR1+ERROR2 00313 ERRSUM = ERRSUM+ERRO12-ERRMAX 00314 AREA = AREA+AREA12-RLIST(MAXERR) 00315 IF(DEFAB1.EQ.ERROR1.OR.DEFAB2.EQ.ERROR2)GO TO 15 00316 IF(DABS(RLIST(MAXERR)-AREA12).GT.0.1D-04*DABS(AREA12) 00317 * .OR.ERRO12.LT.0.99D+00*ERRMAX) GO TO 10 00318 IF(EXTRAP) IROFF2 = IROFF2+1 00319 IF(.NOT.EXTRAP) IROFF1 = IROFF1+1 00320 10 IF(LAST.GT.10.AND.ERRO12.GT.ERRMAX) IROFF3 = IROFF3+1 00321 15 RLIST(MAXERR) = AREA1 00322 RLIST(LAST) = AREA2 00323 ERRBND = DMAX1(EPSABS,EPSREL*DABS(AREA)) 00324 C 00325 C TEST FOR ROUNDOFF ERROR AND EVENTUALLY SET ERROR FLAG. 00326 C 00327 IF(IROFF1+IROFF2.GE.10.OR.IROFF3.GE.20) IER = 2 00328 IF(IROFF2.GE.5) IERRO = 3 00329 C 00330 C SET ERROR FLAG IN THE CASE THAT THE NUMBER OF 00331 C SUBINTERVALS EQUALS LIMIT. 00332 C 00333 IF(LAST.EQ.LIMIT) IER = 1 00334 C 00335 C SET ERROR FLAG IN THE CASE OF BAD INTEGRAND BEHAVIOUR 00336 C AT SOME POINTS OF THE INTEGRATION RANGE. 00337 C 00338 IF(DMAX1(DABS(A1),DABS(B2)).LE.(0.1D+01+0.1D+03*EPMACH)* 00339 * (DABS(A2)+0.1D+04*UFLOW)) IER = 4 00340 C 00341 C APPEND THE NEWLY-CREATED INTERVALS TO THE LIST. 00342 C 00343 IF(ERROR2.GT.ERROR1) GO TO 20 00344 ALIST(LAST) = A2 00345 BLIST(MAXERR) = B1 00346 BLIST(LAST) = B2 00347 ELIST(MAXERR) = ERROR1 00348 ELIST(LAST) = ERROR2 00349 GO TO 30 00350 20 ALIST(MAXERR) = A2 00351 ALIST(LAST) = A1 00352 BLIST(LAST) = B1 00353 RLIST(MAXERR) = AREA2 00354 RLIST(LAST) = AREA1 00355 ELIST(MAXERR) = ERROR2 00356 ELIST(LAST) = ERROR1 00357 C 00358 C CALL SUBROUTINE DQPSRT TO MAINTAIN THE DESCENDING ORDERING 00359 C IN THE LIST OF ERROR ESTIMATES AND SELECT THE SUBINTERVAL 00360 C WITH NRMAX-TH LARGEST ERROR ESTIMATE (TO BE BISECTED NEXT). 00361 C 00362 30 CALL DQPSRT(LIMIT,LAST,MAXERR,ERRMAX,ELIST,IORD,NRMAX) 00363 IF(ERRSUM.LE.ERRBND) GO TO 115 00364 IF(IER.NE.0) GO TO 100 00365 IF(LAST.EQ.2) GO TO 80 00366 IF(NOEXT) GO TO 90 00367 ERLARG = ERLARG-ERLAST 00368 IF(DABS(B1-A1).GT.SMALL) ERLARG = ERLARG+ERRO12 00369 IF(EXTRAP) GO TO 40 00370 C 00371 C TEST WHETHER THE INTERVAL TO BE BISECTED NEXT IS THE 00372 C SMALLEST INTERVAL. 00373 C 00374 IF(DABS(BLIST(MAXERR)-ALIST(MAXERR)).GT.SMALL) GO TO 90 00375 EXTRAP = .TRUE. 00376 NRMAX = 2 00377 40 IF(IERRO.EQ.3.OR.ERLARG.LE.ERTEST) GO TO 60 00378 C 00379 C THE SMALLEST INTERVAL HAS THE LARGEST ERROR. 00380 C BEFORE BISECTING DECREASE THE SUM OF THE ERRORS OVER THE 00381 C LARGER INTERVALS (ERLARG) AND PERFORM EXTRAPOLATION. 00382 C 00383 ID = NRMAX 00384 JUPBND = LAST 00385 IF(LAST.GT.(2+LIMIT/2)) JUPBND = LIMIT+3-LAST 00386 DO 50 K = ID,JUPBND 00387 MAXERR = IORD(NRMAX) 00388 ERRMAX = ELIST(MAXERR) 00389 IF(DABS(BLIST(MAXERR)-ALIST(MAXERR)).GT.SMALL) GO TO 90 00390 NRMAX = NRMAX+1 00391 50 CONTINUE 00392 C 00393 C PERFORM EXTRAPOLATION. 00394 C 00395 60 NUMRL2 = NUMRL2+1 00396 RLIST2(NUMRL2) = AREA 00397 CALL DQELG(NUMRL2,RLIST2,RESEPS,ABSEPS,RES3LA,NRES) 00398 KTMIN = KTMIN+1 00399 IF(KTMIN.GT.5.AND.ABSERR.LT.0.1D-02*ERRSUM) IER = 5 00400 IF(ABSEPS.GE.ABSERR) GO TO 70 00401 KTMIN = 0 00402 ABSERR = ABSEPS 00403 RESULT = RESEPS 00404 CORREC = ERLARG 00405 ERTEST = DMAX1(EPSABS,EPSREL*DABS(RESEPS)) 00406 IF(ABSERR.LE.ERTEST) GO TO 100 00407 C 00408 C PREPARE BISECTION OF THE SMALLEST INTERVAL. 00409 C 00410 70 IF(NUMRL2.EQ.1) NOEXT = .TRUE. 00411 IF(IER.EQ.5) GO TO 100 00412 MAXERR = IORD(1) 00413 ERRMAX = ELIST(MAXERR) 00414 NRMAX = 1 00415 EXTRAP = .FALSE. 00416 SMALL = SMALL*0.5D+00 00417 ERLARG = ERRSUM 00418 GO TO 90 00419 80 SMALL = 0.375D+00 00420 ERLARG = ERRSUM 00421 ERTEST = ERRBND 00422 RLIST2(2) = AREA 00423 90 CONTINUE 00424 C 00425 C SET FINAL RESULT AND ERROR ESTIMATE. 00426 C ------------------------------------ 00427 C 00428 100 IF(ABSERR.EQ.OFLOW) GO TO 115 00429 IF((IER+IERRO).EQ.0) GO TO 110 00430 IF(IERRO.EQ.3) ABSERR = ABSERR+CORREC 00431 IF(IER.EQ.0) IER = 3 00432 IF(RESULT.NE.0.0D+00.AND.AREA.NE.0.0D+00)GO TO 105 00433 IF(ABSERR.GT.ERRSUM)GO TO 115 00434 IF(AREA.EQ.0.0D+00) GO TO 130 00435 GO TO 110 00436 105 IF(ABSERR/DABS(RESULT).GT.ERRSUM/DABS(AREA))GO TO 115 00437 C 00438 C TEST ON DIVERGENCE 00439 C 00440 110 IF(KSGN.EQ.(-1).AND.DMAX1(DABS(RESULT),DABS(AREA)).LE. 00441 * DEFABS*0.1D-01) GO TO 130 00442 IF(0.1D-01.GT.(RESULT/AREA).OR.(RESULT/AREA).GT.0.1D+03. 00443 *OR.ERRSUM.GT.DABS(AREA)) IER = 6 00444 GO TO 130 00445 C 00446 C COMPUTE GLOBAL INTEGRAL SUM. 00447 C 00448 115 RESULT = 0.0D+00 00449 DO 120 K = 1,LAST 00450 RESULT = RESULT+RLIST(K) 00451 120 CONTINUE 00452 ABSERR = ERRSUM 00453 130 NEVAL = 30*LAST-15 00454 IF(INF.EQ.2) NEVAL = 2*NEVAL 00455 IF(IER.GT.2) IER=IER-1 00456 999 RETURN 00457 END
1.7.1
