1 SUBROUTINE DQAGPE(F,A,B,NPTS2,POINTS,EPSABS,EPSREL,LIMIT,RESULT, 2 1 ABSERR,NEVAL,IER,ALIST,BLIST,RLIST,ELIST,PTS,IORD,LEVEL,NDIN, 3 2 LAST) 4C***BEGIN PROLOGUE DQAGPE 5C***DATE WRITTEN 800101 (YYMMDD) 6C***REVISION DATE 840425 (YYMMDD) 7C***REVISION HISTORY (YYMMDD) 8C 000601 Changed DMAX1/DMIN1/DABS to generic MAX/MIN/ABS 9C***CATEGORY NO. H2A2A1 10C***KEYWORDS AUTOMATIC INTEGRATOR,EXTRAPOLATION,GENERAL-PURPOSE, 11C GLOBALLY ADAPTIVE.,SINGULARITIES AT USER SPECIFIED POINTS 12C***AUTHOR PIESSENS, ROBERT, APPLIED MATH. AND PROGR. DIV. - 13C K. U. LEUVEN 14C DE DONCKER, ELISE, APPLIED MATH. AND PROGR. DIV. - 15C K. U. LEUVEN 16C***PURPOSE The routine calculates an approximation result to a given 17C definite integral I = Integral of F over (A,B), hopefully 18C satisfying following claim for accuracy ABS(I-RESULT).LE. 19C MAX(EPSABS,EPSREL*ABS(I)). Break points of the integration 20C interval, where local difficulties of the integrand may 21C occur(e.g. singularities,discontinuities),provided by user. 22C***DESCRIPTION 23C 24C Computation of a definite integral 25C Standard fortran subroutine 26C Double precision version 27C 28C PARAMETERS 29C ON ENTRY 30C F - Double precision 31C Function subprogram defining the integrand 32C function F(X). The actual name for F needs to be 33C declared E X T E R N A L in the driver program. 34C 35C A - Double precision 36C Lower limit of integration 37C 38C B - Double precision 39C Upper limit of integration 40C 41C NPTS2 - Integer 42C Number equal to two more than the number of 43C user-supplied break points within the integration 44C range, NPTS2.GE.2. 45C If NPTS2.LT.2, the routine will end with IER = 6. 46C 47C POINTS - Double precision 48C Vector of dimension NPTS2, the first (NPTS2-2) 49C elements of which are the user provided break 50C POINTS. If these POINTS do not constitute an 51C ascending sequence there will be an automatic 52C sorting. 53C 54C EPSABS - Double precision 55C Absolute accuracy requested 56C EPSREL - Double precision 57C Relative accuracy requested 58C If EPSABS.LE.0 59C and EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28), 60C the routine will end with IER = 6. 61C 62C LIMIT - Integer 63C Gives an upper bound on the number of subintervals 64C in the partition of (A,B), LIMIT.GE.NPTS2 65C If LIMIT.LT.NPTS2, the routine will end with 66C IER = 6. 67C 68C ON RETURN 69C RESULT - Double precision 70C Approximation to the integral 71C 72C ABSERR - Double precision 73C Estimate of the modulus of the absolute error, 74C which should equal or exceed ABS(I-RESULT) 75C 76C NEVAL - Integer 77C Number of integrand evaluations 78C 79C IER - Integer 80C IER = 0 Normal and reliable termination of the 81C routine. It is assumed that the requested 82C accuracy has been achieved. 83C IER.GT.0 Abnormal termination of the routine. 84C The estimates for integral and error are 85C less reliable. It is assumed that the 86C requested accuracy has not been achieved. 87C ERROR MESSAGES 88C IER = 1 Maximum number of subdivisions allowed 89C has been achieved. One can allow more 90C subdivisions by increasing the value of 91C LIMIT (and taking the according dimension 92C adjustments into account). However, if 93C this yields no improvement it is advised 94C to analyze the integrand in order to 95C determine the integration difficulties. If 96C the position of a local difficulty can be 97C determined (i.e. SINGULARITY, 98C DISCONTINUITY within the interval), it 99C should be supplied to the routine as an 100C element of the vector points. If necessary 101C an appropriate special-purpose integrator 102C must be used, which is designed for 103C handling the type of difficulty involved. 104C = 2 The occurrence of roundoff error is 105C detected, which prevents the requested 106C tolerance from being achieved. 107C The error may be under-estimated. 108C = 3 Extremely bad integrand behaviour occurs 109C At some points of the integration 110C interval. 111C = 4 The algorithm does not converge. 112C Roundoff error is detected in the 113C extrapolation table. It is presumed that 114C the requested tolerance cannot be 115C achieved, and that the returned result is 116C the best which can be obtained. 117C = 5 The integral is probably divergent, or 118C slowly convergent. It must be noted that 119C divergence can occur with any other value 120C of IER.GT.0. 121C = 6 The input is invalid because 122C NPTS2.LT.2 or 123C Break points are specified outside 124C the integration range or 125C (EPSABS.LE.0 and 126C EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28)) 127C or LIMIT.LT.NPTS2. 128C RESULT, ABSERR, NEVAL, LAST, RLIST(1), 129C and ELIST(1) are set to zero. ALIST(1) and 130C BLIST(1) are set to A and B respectively. 131C 132C ALIST - Double precision 133C Vector of dimension at least LIMIT, the first 134C LAST elements of which are the left end points 135C of the subintervals in the partition of the given 136C integration range (A,B) 137C 138C BLIST - Double precision 139C Vector of dimension at least LIMIT, the first 140C LAST elements of which are the right end points 141C of the subintervals in the partition of the given 142C integration range (A,B) 143C 144C RLIST - Double precision 145C Vector of dimension at least LIMIT, the first 146C LAST elements of which are the integral 147C approximations on the subintervals 148C 149C ELIST - Double precision 150C Vector of dimension at least LIMIT, the first 151C LAST elements of which are the moduli of the 152C absolute error estimates on the subintervals 153C 154C PTS - Double precision 155C Vector of dimension at least NPTS2, containing the 156C integration limits and the break points of the 157C interval in ascending sequence. 158C 159C LEVEL - Integer 160C Vector of dimension at least LIMIT, containing the 161C subdivision levels of the subinterval, i.e. if 162C (AA,BB) is a subinterval of (P1,P2) where P1 as 163C well as P2 is a user-provided break point or 164C integration limit, then (AA,BB) has level L if 165C ABS(BB-AA) = ABS(P2-P1)*2**(-L). 166C 167C NDIN - Integer 168C Vector of dimension at least NPTS2, after first 169C integration over the intervals (PTS(I)),PTS(I+1), 170C I = 0,1, ..., NPTS2-2, the error estimates over 171C some of the intervals may have been increased 172C artificially, in order to put their subdivision 173C forward. If this happens for the subinterval 174C numbered K, NDIN(K) is put to 1, otherwise 175C NDIN(K) = 0. 176C 177C IORD - Integer 178C Vector of dimension at least LIMIT, the first K 179C elements of which are pointers to the 180C error estimates over the subintervals, 181C such that ELIST(IORD(1)), ..., ELIST(IORD(K)) 182C form a decreasing sequence, with K = LAST 183C If LAST.LE.(LIMIT/2+2), and K = LIMIT+1-LAST 184C otherwise 185C 186C LAST - Integer 187C Number of subintervals actually produced in the 188C subdivisions process 189C***REFERENCES (NONE) 190C***ROUTINES CALLED D1MACH,DQELG,DQK21,DQPSRT 191C***END PROLOGUE DQAGPE 192 DOUBLE PRECISION A,ABSEPS,ABSERR,ALIST,AREA,AREA1,AREA12,AREA2,A1, 193 1 A2,B,BLIST,B1,B2,CORREC,DEFABS,DEFAB1,DEFAB2, 194 2 DRES,D1MACH,ELIST,EPMACH,EPSABS,EPSREL,ERLARG,ERLAST,ERRBND, 195 3 ERRMAX,ERROR1,ERRO12,ERROR2,ERRSUM,ERTEST,F,OFLOW,POINTS,PTS, 196 4 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 1 JLOW,JUPBND,K,KSGN,KTMIN,LAST,LEVCUR,LEVEL,LEVMAX,LIMIT,MAXERR, 199 2 NDIN,NEVAL,NINT,NINTP1,NPTS,NPTS2,NRES,NRMAX,NUMRL2 200 LOGICAL EXTRAP,NOEXT 201C 202C 203 DIMENSION ALIST(LIMIT),BLIST(LIMIT),ELIST(LIMIT),IORD(LIMIT), 204 1 LEVEL(LIMIT),NDIN(NPTS2),POINTS(NPTS2),PTS(NPTS2),RES3LA(3), 205 2 RLIST(LIMIT),RLIST2(52) 206C 207 EXTERNAL F 208C 209C THE DIMENSION OF RLIST2 IS DETERMINED BY THE VALUE OF 210C LIMEXP IN SUBROUTINE EPSALG (RLIST2 SHOULD BE OF DIMENSION 211C (LIMEXP+2) AT LEAST). 212C 213C 214C LIST OF MAJOR VARIABLES 215C ----------------------- 216C 217C ALIST - LIST OF LEFT END POINTS OF ALL SUBINTERVALS 218C CONSIDERED UP TO NOW 219C BLIST - LIST OF RIGHT END POINTS OF ALL SUBINTERVALS 220C CONSIDERED UP TO NOW 221C RLIST(I) - APPROXIMATION TO THE INTEGRAL OVER 222C (ALIST(I),BLIST(I)) 223C RLIST2 - ARRAY OF DIMENSION AT LEAST LIMEXP+2 224C CONTAINING THE PART OF THE EPSILON TABLE WHICH 225C IS STILL NEEDED FOR FURTHER COMPUTATIONS 226C ELIST(I) - ERROR ESTIMATE APPLYING TO RLIST(I) 227C MAXERR - POINTER TO THE INTERVAL WITH LARGEST ERROR 228C ESTIMATE 229C ERRMAX - ELIST(MAXERR) 230C ERLAST - ERROR ON THE INTERVAL CURRENTLY SUBDIVIDED 231C (BEFORE THAT SUBDIVISION HAS TAKEN PLACE) 232C AREA - SUM OF THE INTEGRALS OVER THE SUBINTERVALS 233C ERRSUM - SUM OF THE ERRORS OVER THE SUBINTERVALS 234C ERRBND - REQUESTED ACCURACY MAX(EPSABS,EPSREL* 235C ABS(RESULT)) 236C *****1 - VARIABLE FOR THE LEFT SUBINTERVAL 237C *****2 - VARIABLE FOR THE RIGHT SUBINTERVAL 238C LAST - INDEX FOR SUBDIVISION 239C NRES - NUMBER OF CALLS TO THE EXTRAPOLATION ROUTINE 240C NUMRL2 - NUMBER OF ELEMENTS IN RLIST2. IF AN APPROPRIATE 241C APPROXIMATION TO THE COMPOUNDED INTEGRAL HAS 242C BEEN OBTAINED, IT IS PUT IN RLIST2(NUMRL2) AFTER 243C NUMRL2 HAS BEEN INCREASED BY ONE. 244C ERLARG - SUM OF THE ERRORS OVER THE INTERVALS LARGER 245C THAN THE SMALLEST INTERVAL CONSIDERED UP TO NOW 246C EXTRAP - LOGICAL VARIABLE DENOTING THAT THE ROUTINE 247C IS ATTEMPTING TO PERFORM EXTRAPOLATION. I.E. 248C BEFORE SUBDIVIDING THE SMALLEST INTERVAL WE 249C TRY TO DECREASE THE VALUE OF ERLARG. 250C NOEXT - LOGICAL VARIABLE DENOTING THAT EXTRAPOLATION IS 251C NO LONGER ALLOWED (TRUE-VALUE) 252C 253C MACHINE DEPENDENT CONSTANTS 254C --------------------------- 255C 256C EPMACH IS THE LARGEST RELATIVE SPACING. 257C UFLOW IS THE SMALLEST POSITIVE MAGNITUDE. 258C OFLOW IS THE LARGEST POSITIVE MAGNITUDE. 259C 260C***FIRST EXECUTABLE STATEMENT DQAGPE 261 EPMACH = D1MACH(4) 262C 263C TEST ON VALIDITY OF PARAMETERS 264C ----------------------------- 265C 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 1 EPSREL.LT.MAX(0.5D+02*EPMACH,0.5D-28))) IER = 6 280 IF(IER.EQ.6) GO TO 999 281C 282C IF ANY BREAK POINTS ARE PROVIDED, SORT THEM INTO AN 283C ASCENDING SEQUENCE. 284C 285 SIGN = 1.0D+00 286 IF(A.GT.B) SIGN = -1.0D+00 287 PTS(1) = MIN(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) = MAX(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.MIN(A,B).OR.PTS(NINTP1).NE.MAX(A,B)) IER = 6 306 IF(IER.EQ.6) GO TO 999 307C 308C COMPUTE FIRST INTEGRAL AND ERROR APPROXIMATIONS. 309C ------------------------------------------------ 310C 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) 315 ABSERR = ABSERR+ERROR1 316 RESULT = RESULT+AREA1 317 NDIN(I) = 0 318 IF(ERROR1.EQ.RESA.AND.ERROR1.NE.0.0D+00) NDIN(I) = 1 319 RESABS = RESABS+DEFABS 320 LEVEL(I) = 0 321 ELIST(I) = ERROR1 322 ALIST(I) = A1 323 BLIST(I) = B1 324 RLIST(I) = AREA1 325 IORD(I) = I 326 A1 = B1 327 50 CONTINUE 328 ERRSUM = 0.0D+00 329 DO 55 I = 1,NINT 330 IF(NDIN(I).EQ.1) ELIST(I) = ABSERR 331 ERRSUM = ERRSUM+ELIST(I) 332 55 CONTINUE 333C 334C TEST ON ACCURACY. 335C 336 LAST = NINT 337 NEVAL = 21*NINT 338 DRES = ABS(RESULT) 339 ERRBND = MAX(EPSABS,EPSREL*DRES) 340 IF(ABSERR.LE.0.1D+03*EPMACH*RESABS.AND.ABSERR.GT.ERRBND) IER = 2 341 IF(NINT.EQ.1) GO TO 80 342 DO 70 I = 1,NPTS 343 JLOW = I+1 344 IND1 = IORD(I) 345 DO 60 J = JLOW,NINT 346 IND2 = IORD(J) 347 IF(ELIST(IND1).GT.ELIST(IND2)) GO TO 60 348 IND1 = IND2 349 K = J 350 60 CONTINUE 351 IF(IND1.EQ.IORD(I)) GO TO 70 352 IORD(K) = IORD(I) 353 IORD(I) = IND1 354 70 CONTINUE 355 IF(LIMIT.LT.NPTS2) IER = 1 356 80 IF(IER.NE.0.OR.ABSERR.LE.ERRBND) GO TO 210 357C 358C INITIALIZATION 359C -------------- 360C 361 RLIST2(1) = RESULT 362 MAXERR = IORD(1) 363 ERRMAX = ELIST(MAXERR) 364 AREA = RESULT 365 NRMAX = 1 366 NRES = 0 367 NUMRL2 = 1 368 KTMIN = 0 369 EXTRAP = .FALSE. 370 NOEXT = .FALSE. 371 ERLARG = ERRSUM 372 ERTEST = ERRBND 373 LEVMAX = 1 374 IROFF1 = 0 375 IROFF2 = 0 376 IROFF3 = 0 377 IERRO = 0 378 UFLOW = D1MACH(1) 379 OFLOW = D1MACH(2) 380 ABSERR = OFLOW 381 KSGN = -1 382 IF(DRES.GE.(0.1D+01-0.5D+02*EPMACH)*RESABS) KSGN = 1 383C 384C MAIN DO-LOOP 385C ------------ 386C 387 DO 160 LAST = NPTS2,LIMIT 388C 389C BISECT THE SUBINTERVAL WITH THE NRMAX-TH LARGEST ERROR 390C ESTIMATE. 391C 392 LEVCUR = LEVEL(MAXERR)+1 393 A1 = ALIST(MAXERR) 394 B1 = 0.5D+00*(ALIST(MAXERR)+BLIST(MAXERR)) 395 A2 = B1 396 B2 = BLIST(MAXERR) 397 ERLAST = ERRMAX 398 CALL DQK21(F,A1,B1,AREA1,ERROR1,RESA,DEFAB1) 399 CALL DQK21(F,A2,B2,AREA2,ERROR2,RESA,DEFAB2) 400C 401C IMPROVE PREVIOUS APPROXIMATIONS TO INTEGRAL 402C AND ERROR AND TEST FOR ACCURACY. 403C 404 NEVAL = NEVAL+42 405 AREA12 = AREA1+AREA2 406 ERRO12 = ERROR1+ERROR2 407 ERRSUM = ERRSUM+ERRO12-ERRMAX 408 AREA = AREA+AREA12-RLIST(MAXERR) 409 IF(DEFAB1.EQ.ERROR1.OR.DEFAB2.EQ.ERROR2) GO TO 95 410 IF(ABS(RLIST(MAXERR)-AREA12).GT.0.1D-04*ABS(AREA12) 411 1 .OR.ERRO12.LT.0.99D+00*ERRMAX) GO TO 90 412 IF(EXTRAP) IROFF2 = IROFF2+1 413 IF(.NOT.EXTRAP) IROFF1 = IROFF1+1 414 90 IF(LAST.GT.10.AND.ERRO12.GT.ERRMAX) IROFF3 = IROFF3+1 415 95 LEVEL(MAXERR) = LEVCUR 416 LEVEL(LAST) = LEVCUR 417 RLIST(MAXERR) = AREA1 418 RLIST(LAST) = AREA2 419 ERRBND = MAX(EPSABS,EPSREL*ABS(AREA)) 420C 421C TEST FOR ROUNDOFF ERROR AND EVENTUALLY SET ERROR FLAG. 422C 423 IF(IROFF1+IROFF2.GE.10.OR.IROFF3.GE.20) IER = 2 424 IF(IROFF2.GE.5) IERRO = 3 425C 426C SET ERROR FLAG IN THE CASE THAT THE NUMBER OF 427C SUBINTERVALS EQUALS LIMIT. 428C 429 IF(LAST.EQ.LIMIT) IER = 1 430C 431C SET ERROR FLAG IN THE CASE OF BAD INTEGRAND BEHAVIOUR 432C AT A POINT OF THE INTEGRATION RANGE 433C 434 IF(MAX(ABS(A1),ABS(B2)).LE.(0.1D+01+0.1D+03*EPMACH)* 435 1 (ABS(A2)+0.1D+04*UFLOW)) IER = 4 436C 437C APPEND THE NEWLY-CREATED INTERVALS TO THE LIST. 438C 439 IF(ERROR2.GT.ERROR1) GO TO 100 440 ALIST(LAST) = A2 441 BLIST(MAXERR) = B1 442 BLIST(LAST) = B2 443 ELIST(MAXERR) = ERROR1 444 ELIST(LAST) = ERROR2 445 GO TO 110 446 100 ALIST(MAXERR) = A2 447 ALIST(LAST) = A1 448 BLIST(LAST) = B1 449 RLIST(MAXERR) = AREA2 450 RLIST(LAST) = AREA1 451 ELIST(MAXERR) = ERROR2 452 ELIST(LAST) = ERROR1 453C 454C CALL SUBROUTINE DQPSRT TO MAINTAIN THE DESCENDING ORDERING 455C IN THE LIST OF ERROR ESTIMATES AND SELECT THE SUBINTERVAL 456C WITH NRMAX-TH LARGEST ERROR ESTIMATE (TO BE BISECTED NEXT). 457C 458 110 CALL DQPSRT(LIMIT,LAST,MAXERR,ERRMAX,ELIST,IORD,NRMAX) 459C ***JUMP OUT OF DO-LOOP 460 IF(ERRSUM.LE.ERRBND) GO TO 190 461C ***JUMP OUT OF DO-LOOP 462 IF(IER.NE.0) GO TO 170 463 IF(NOEXT) GO TO 160 464 ERLARG = ERLARG-ERLAST 465 IF(LEVCUR+1.LE.LEVMAX) ERLARG = ERLARG+ERRO12 466 IF(EXTRAP) GO TO 120 467C 468C TEST WHETHER THE INTERVAL TO BE BISECTED NEXT IS THE 469C SMALLEST INTERVAL. 470C 471 IF(LEVEL(MAXERR)+1.LE.LEVMAX) GO TO 160 472 EXTRAP = .TRUE. 473 NRMAX = 2 474 120 IF(IERRO.EQ.3.OR.ERLARG.LE.ERTEST) GO TO 140 475C 476C THE SMALLEST INTERVAL HAS THE LARGEST ERROR. 477C BEFORE BISECTING DECREASE THE SUM OF THE ERRORS OVER 478C THE LARGER INTERVALS (ERLARG) AND PERFORM EXTRAPOLATION. 479C 480 ID = NRMAX 481 JUPBND = LAST 482 IF(LAST.GT.(2+LIMIT/2)) JUPBND = LIMIT+3-LAST 483 DO 130 K = ID,JUPBND 484 MAXERR = IORD(NRMAX) 485 ERRMAX = ELIST(MAXERR) 486C ***JUMP OUT OF DO-LOOP 487 IF(LEVEL(MAXERR)+1.LE.LEVMAX) GO TO 160 488 NRMAX = NRMAX+1 489 130 CONTINUE 490C 491C PERFORM EXTRAPOLATION. 492C 493 140 NUMRL2 = NUMRL2+1 494 RLIST2(NUMRL2) = AREA 495 IF(NUMRL2.LE.2) GO TO 155 496 CALL DQELG(NUMRL2,RLIST2,RESEPS,ABSEPS,RES3LA,NRES) 497 KTMIN = KTMIN+1 498 IF(KTMIN.GT.5.AND.ABSERR.LT.0.1D-02*ERRSUM) IER = 5 499 IF(ABSEPS.GE.ABSERR) GO TO 150 500 KTMIN = 0 501 ABSERR = ABSEPS 502 RESULT = RESEPS 503 CORREC = ERLARG 504 ERTEST = MAX(EPSABS,EPSREL*ABS(RESEPS)) 505C ***JUMP OUT OF DO-LOOP 506 IF(ABSERR.LT.ERTEST) GO TO 170 507C 508C PREPARE BISECTION OF THE SMALLEST INTERVAL. 509C 510 150 IF(NUMRL2.EQ.1) NOEXT = .TRUE. 511 IF(IER.GE.5) GO TO 170 512 155 MAXERR = IORD(1) 513 ERRMAX = ELIST(MAXERR) 514 NRMAX = 1 515 EXTRAP = .FALSE. 516 LEVMAX = LEVMAX+1 517 ERLARG = ERRSUM 518 160 CONTINUE 519C 520C SET THE FINAL RESULT. 521C --------------------- 522C 523C 524 170 IF(ABSERR.EQ.OFLOW) GO TO 190 525 IF((IER+IERRO).EQ.0) GO TO 180 526 IF(IERRO.EQ.3) ABSERR = ABSERR+CORREC 527 IF(IER.EQ.0) IER = 3 528 IF(RESULT.NE.0.0D+00.AND.AREA.NE.0.0D+00)GO TO 175 529 IF(ABSERR.GT.ERRSUM)GO TO 190 530 IF(AREA.EQ.0.0D+00) GO TO 210 531 GO TO 180 532 175 IF(ABSERR/ABS(RESULT).GT.ERRSUM/ABS(AREA))GO TO 190 533C 534C TEST ON DIVERGENCE. 535C 536 180 IF(KSGN.EQ.(-1).AND.MAX(ABS(RESULT),ABS(AREA)).LE. 537 1 RESABS*0.1D-01) GO TO 210 538 IF(0.1D-01.GT.(RESULT/AREA).OR.(RESULT/AREA).GT.0.1D+03.OR. 539 1 ERRSUM.GT.ABS(AREA)) IER = 6 540 GO TO 210 541C 542C COMPUTE GLOBAL INTEGRAL SUM. 543C 544 190 RESULT = 0.0D+00 545 DO 200 K = 1,LAST 546 RESULT = RESULT+RLIST(K) 547 200 CONTINUE 548 ABSERR = ERRSUM 549 210 IF(IER.GT.2) IER = IER-1 550 RESULT = RESULT*SIGN 551 999 RETURN 552 END 553