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