src/numerical/slatec/fortran/dqagie.f

   1 *DECK DQAGIE
   2       SUBROUTINE DQAGIE (F, BOUND, INF, EPSABS, EPSREL, LIMIT, RESULT,
   3      +   ABSERR, NEVAL, IER, ALIST, BLIST, RLIST, ELIST, IORD, LAST)
   4 C***BEGIN PROLOGUE  DQAGIE
   5 C***PURPOSE  The routine calculates an approximation result to a given
   6 C            integral   I = Integral of F over (BOUND,+INFINITY)
   7 C            or I = Integral of F over (-INFINITY,BOUND)
   8 C            or I = Integral of F over (-INFINITY,+INFINITY),
   9 C            hopefully satisfying following claim for accuracy
  10 C            ABS(I-RESULT).LE.MAX(EPSABS,EPSREL*ABS(I))
  11 C***LIBRARY   SLATEC (QUADPACK)
  12 C***CATEGORY  H2A3A1, H2A4A1
  13 C***TYPE      DOUBLE PRECISION (QAGIE-S, DQAGIE-D)
  14 C***KEYWORDS  AUTOMATIC INTEGRATOR, EXTRAPOLATION, GENERAL-PURPOSE,
  15 C             GLOBALLY ADAPTIVE, INFINITE INTERVALS, QUADPACK,
  16 C             QUADRATURE, TRANSFORMATION
  17 C***AUTHOR  Piessens, Robert
  18 C             Applied Mathematics and Programming Division
  19 C             K. U. Leuven
  20 C           de Doncker, Elise
  21 C             Applied Mathematics and Programming Division
  22 C             K. U. Leuven
  23 C***DESCRIPTION
  24 C
  25 C Integration over infinite intervals
  26 C Standard fortran subroutine
  27 C
  28 C            F      - Double precision
  29 C                     Function subprogram defining the integrand
  30 C                     function F(X). The actual name for F needs to be
  31 C                     declared E X T E R N A L in the driver program.
  32 C
  33 C            BOUND  - Double precision
  34 C                     Finite bound of integration range
  35 C                     (has no meaning if interval is doubly-infinite)
  36 C
  37 C            INF    - Double precision
  38 C                     Indicating the kind of integration range involved
  39 C                     INF = 1 corresponds to  (BOUND,+INFINITY),
  40 C                     INF = -1            to  (-INFINITY,BOUND),
  41 C                     INF = 2             to (-INFINITY,+INFINITY).
  42 C
  43 C            EPSABS - Double precision
  44 C                     Absolute accuracy requested
  45 C            EPSREL - Double precision
  46 C                     Relative accuracy requested
  47 C                     If  EPSABS.LE.0
  48 C                     and EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28),
  49 C                     the routine will end with IER = 6.
  50 C
  51 C            LIMIT  - Integer
  52 C                     Gives an upper bound on the number of subintervals
  53 C                     in the partition of (A,B), LIMIT.GE.1
  54 C
  55 C         ON RETURN
  56 C            RESULT - Double precision
  57 C                     Approximation to the integral
  58 C
  59 C            ABSERR - Double precision
  60 C                     Estimate of the modulus of the absolute error,
  61 C                     which should equal or exceed ABS(I-RESULT)
  62 C
  63 C            NEVAL  - Integer
  64 C                     Number of integrand evaluations
  65 C
  66 C            IER    - Integer
  67 C                     IER = 0 Normal and reliable termination of the
  68 C                             routine. It is assumed that the requested
  69 C                             accuracy has been achieved.
  70 C                   - IER.GT.0 Abnormal termination of the routine. The
  71 C                             estimates for result and error are less
  72 C                             reliable. It is assumed that the requested
  73 C                             accuracy has not been achieved.
  74 C            ERROR MESSAGES
  75 C                     IER = 1 Maximum number of subdivisions allowed
  76 C                             has been achieved. One can allow more
  77 C                             subdivisions by increasing the value of
  78 C                             LIMIT (and taking the according dimension
  79 C                             adjustments into account).  However, if
  80 C                             this yields no improvement it is advised
  81 C                             to analyze the integrand in order to
  82 C                             determine the integration difficulties.
  83 C                             If the position of a local difficulty can
  84 C                             be determined (e.g. SINGULARITY,
  85 C                             DISCONTINUITY within the interval) one
  86 C                             will probably gain from splitting up the
  87 C                             interval at this point and calling the
  88 C                             integrator on the subranges. If possible,
  89 C                             an appropriate special-purpose integrator
  90 C                             should be used, which is designed for
  91 C                             handling the type of difficulty involved.
  92 C                         = 2 The occurrence of roundoff error is
  93 C                             detected, which prevents the requested
  94 C                             tolerance from being achieved.
  95 C                             The error may be under-estimated.
  96 C                         = 3 Extremely bad integrand behaviour occurs
  97 C                             at some points of the integration
  98 C                             interval.
  99 C                         = 4 The algorithm does not converge.
 100 C                             Roundoff error is detected in the
 101 C                             extrapolation table.
 102 C                             It is assumed that the requested tolerance
 103 C                             cannot be achieved, and that the returned
 104 C                             result is the best which can be obtained.
 105 C                         = 5 The integral is probably divergent, or
 106 C                             slowly convergent. It must be noted that
 107 C                             divergence can occur with any other value
 108 C                             of IER.
 109 C                         = 6 The input is invalid, because
 110 C                             (EPSABS.LE.0 and
 111 C                              EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28),
 112 C                             RESULT, ABSERR, NEVAL, LAST, RLIST(1),
 113 C                             ELIST(1) and IORD(1) are set to zero.
 114 C                             ALIST(1) and BLIST(1) are set to 0
 115 C                             and 1 respectively.
 116 C
 117 C            ALIST  - Double precision
 118 C                     Vector of dimension at least LIMIT, the first
 119 C                      LAST  elements of which are the left
 120 C                     end points of the subintervals in the partition
 121 C                     of the transformed integration range (0,1).
 122 C
 123 C            BLIST  - Double precision
 124 C                     Vector of dimension at least LIMIT, the first
 125 C                      LAST  elements of which are the right
 126 C                     end points of the subintervals in the partition
 127 C                     of the transformed integration range (0,1).
 128 C
 129 C            RLIST  - Double precision
 130 C                     Vector of dimension at least LIMIT, the first
 131 C                      LAST  elements of which are the integral
 132 C                     approximations on the subintervals
 133 C
 134 C            ELIST  - Double precision
 135 C                     Vector of dimension at least LIMIT,  the first
 136 C                     LAST elements of which are the moduli of the
 137 C                     absolute error estimates on the subintervals
 138 C
 139 C            IORD   - Integer
 140 C                     Vector of dimension LIMIT, the first K
 141 C                     elements of which are pointers to the
 142 C                     error estimates over the subintervals,
 143 C                     such that ELIST(IORD(1)), ..., ELIST(IORD(K))
 144 C                     form a decreasing sequence, with K = LAST
 145 C                     If LAST.LE.(LIMIT/2+2), and K = LIMIT+1-LAST
 146 C                     otherwise
 147 C
 148 C            LAST   - Integer
 149 C                     Number of subintervals actually produced
 150 C                     in the subdivision process
 151 C
 152 C***REFERENCES  (NONE)
 153 C***ROUTINES CALLED  D1MACH, DQELG, DQK15I, DQPSRT
 154 C***REVISION HISTORY  (YYMMDD)
 155 C   800101  DATE WRITTEN
 156 C   890531  Changed all specific intrinsics to generic.  (WRB)
 157 C   890831  Modified array declarations.  (WRB)
 158 C   890831  REVISION DATE from Version 3.2
 159 C   891214  Prologue converted to Version 4.0 format.  (BAB)
 160 C***END PROLOGUE  DQAGIE
 161       DOUBLE PRECISION ABSEPS,ABSERR,ALIST,AREA,AREA1,AREA12,AREA2,A1,
 162      1  A2,BLIST,BOUN,BOUND,B1,B2,CORREC,DEFABS,DEFAB1,DEFAB2,
 163      2  DRES,D1MACH,ELIST,EPMACH,EPSABS,EPSREL,ERLARG,ERLAST,
 164      3  ERRBND,ERRMAX,ERROR1,ERROR2,ERRO12,ERRSUM,ERTEST,F,OFLOW,RESABS,
 165      4  RESEPS,RESULT,RES3LA,RLIST,RLIST2,SMALL,UFLOW
 166       INTEGER ID,IER,IERRO,INF,IORD,IROFF1,IROFF2,IROFF3,JUPBND,K,KSGN,
 167      1  KTMIN,LAST,LIMIT,MAXERR,NEVAL,NRES,NRMAX,NUMRL2
 168       LOGICAL EXTRAP,NOEXT
 169 C
 170       DIMENSION ALIST(*),BLIST(*),ELIST(*),IORD(*),
 171      1  RES3LA(3),RLIST(*),RLIST2(52)
 172 C
 173       EXTERNAL F
 174 C
 175 C            THE DIMENSION OF RLIST2 IS DETERMINED BY THE VALUE OF
 176 C            LIMEXP IN SUBROUTINE DQELG.
 177 C
 178 C
 179 C            LIST OF MAJOR VARIABLES
 180 C            -----------------------
 181 C
 182 C           ALIST     - LIST OF LEFT END POINTS OF ALL SUBINTERVALS
 183 C                       CONSIDERED UP TO NOW
 184 C           BLIST     - LIST OF RIGHT END POINTS OF ALL SUBINTERVALS
 185 C                       CONSIDERED UP TO NOW
 186 C           RLIST(I)  - APPROXIMATION TO THE INTEGRAL OVER
 187 C                       (ALIST(I),BLIST(I))
 188 C           RLIST2    - ARRAY OF DIMENSION AT LEAST (LIMEXP+2),
 189 C                       CONTAINING THE PART OF THE EPSILON TABLE
 190 C                       WHICH IS STILL NEEDED FOR FURTHER COMPUTATIONS
 191 C           ELIST(I)  - ERROR ESTIMATE APPLYING TO RLIST(I)
 192 C           MAXERR    - POINTER TO THE INTERVAL WITH LARGEST ERROR
 193 C                       ESTIMATE
 194 C           ERRMAX    - ELIST(MAXERR)
 195 C           ERLAST    - ERROR ON THE INTERVAL CURRENTLY SUBDIVIDED
 196 C                       (BEFORE THAT SUBDIVISION HAS TAKEN PLACE)
 197 C           AREA      - SUM OF THE INTEGRALS OVER THE SUBINTERVALS
 198 C           ERRSUM    - SUM OF THE ERRORS OVER THE SUBINTERVALS
 199 C           ERRBND    - REQUESTED ACCURACY MAX(EPSABS,EPSREL*
 200 C                       ABS(RESULT))
 201 C           *****1    - VARIABLE FOR THE LEFT SUBINTERVAL
 202 C           *****2    - VARIABLE FOR THE RIGHT SUBINTERVAL
 203 C           LAST      - INDEX FOR SUBDIVISION
 204 C           NRES      - NUMBER OF CALLS TO THE EXTRAPOLATION ROUTINE
 205 C           NUMRL2    - NUMBER OF ELEMENTS CURRENTLY IN RLIST2. IF AN
 206 C                       APPROPRIATE APPROXIMATION TO THE COMPOUNDED
 207 C                       INTEGRAL HAS BEEN OBTAINED, IT IS PUT IN
 208 C                       RLIST2(NUMRL2) AFTER NUMRL2 HAS BEEN INCREASED
 209 C                       BY ONE.
 210 C           SMALL     - LENGTH OF THE SMALLEST INTERVAL CONSIDERED UP
 211 C                       TO NOW, MULTIPLIED BY 1.5
 212 C           ERLARG    - SUM OF THE ERRORS OVER THE INTERVALS LARGER
 213 C                       THAN THE SMALLEST INTERVAL CONSIDERED UP TO NOW
 214 C           EXTRAP    - LOGICAL VARIABLE DENOTING THAT THE ROUTINE
 215 C                       IS ATTEMPTING TO PERFORM EXTRAPOLATION. I.E.
 216 C                       BEFORE SUBDIVIDING THE SMALLEST INTERVAL WE
 217 C                       TRY TO DECREASE THE VALUE OF ERLARG.
 218 C           NOEXT     - LOGICAL VARIABLE DENOTING THAT EXTRAPOLATION
 219 C                       IS NO LONGER ALLOWED (TRUE-VALUE)
 220 C
 221 C            MACHINE DEPENDENT CONSTANTS
 222 C            ---------------------------
 223 C
 224 C           EPMACH IS THE LARGEST RELATIVE SPACING.
 225 C           UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
 226 C           OFLOW IS THE LARGEST POSITIVE MAGNITUDE.
 227 C
 228 C***FIRST EXECUTABLE STATEMENT  DQAGIE
 229       EPMACH = D1MACH(4)
 230 C
 231 C           TEST ON VALIDITY OF PARAMETERS
 232 C           -----------------------------
 233 C
 234       IER = 0
 235       NEVAL = 0
 236       LAST = 0
 237       RESULT = 0.0D+00
 238       ABSERR = 0.0D+00
 239       ALIST(1) = 0.0D+00
 240       BLIST(1) = 0.1D+01
 241       RLIST(1) = 0.0D+00
 242       ELIST(1) = 0.0D+00
 243       IORD(1) = 0
 244       IF(EPSABS.LE.0.0D+00.AND.EPSREL.LT.MAX(0.5D+02*EPMACH,0.5D-28))
 245      1  IER = 6
 246        IF(IER.EQ.6) GO TO 999
 247 C
 248 C
 249 C           FIRST APPROXIMATION TO THE INTEGRAL
 250 C           -----------------------------------
 251 C
 252 C           DETERMINE THE INTERVAL TO BE MAPPED ONTO (0,1).
 253 C           IF INF = 2 THE INTEGRAL IS COMPUTED AS I = I1+I2, WHERE
 254 C           I1 = INTEGRAL OF F OVER (-INFINITY,0),
 255 C           I2 = INTEGRAL OF F OVER (0,+INFINITY).
 256 C
 257       BOUN = BOUND
 258       IF(INF.EQ.2) BOUN = 0.0D+00
 259       CALL DQK15I(F,BOUN,INF,0.0D+00,0.1D+01,RESULT,ABSERR,
 260      1  DEFABS,RESABS)
 261 C
 262 C           TEST ON ACCURACY
 263 C
 264       LAST = 1
 265       RLIST(1) = RESULT
 266       ELIST(1) = ABSERR
 267       IORD(1) = 1
 268       DRES = ABS(RESULT)
 269       ERRBND = MAX(EPSABS,EPSREL*DRES)
 270       IF(ABSERR.LE.1.0D+02*EPMACH*DEFABS.AND.ABSERR.GT.ERRBND) IER = 2
 271       IF(LIMIT.EQ.1) IER = 1
 272       IF(IER.NE.0.OR.(ABSERR.LE.ERRBND.AND.ABSERR.NE.RESABS).OR.
 273      1  ABSERR.EQ.0.0D+00) GO TO 130
 274 C
 275 C           INITIALIZATION
 276 C           --------------
 277 C
 278       UFLOW = D1MACH(1)
 279       OFLOW = D1MACH(2)
 280       RLIST2(1) = RESULT
 281       ERRMAX = ABSERR
 282       MAXERR = 1
 283       AREA = RESULT
 284       ERRSUM = ABSERR
 285       ABSERR = OFLOW
 286       NRMAX = 1
 287       NRES = 0
 288       KTMIN = 0
 289       NUMRL2 = 2
 290       EXTRAP = .FALSE.
 291       NOEXT = .FALSE.
 292       IERRO = 0
 293       IROFF1 = 0
 294       IROFF2 = 0
 295       IROFF3 = 0
 296       KSGN = -1
 297       IF(DRES.GE.(0.1D+01-0.5D+02*EPMACH)*DEFABS) KSGN = 1
 298 C
 299 C           MAIN DO-LOOP
 300 C           ------------
 301 C
 302       DO 90 LAST = 2,LIMIT
 303 C
 304 C           BISECT THE SUBINTERVAL WITH NRMAX-TH LARGEST ERROR ESTIMATE.
 305 C
 306         A1 = ALIST(MAXERR)
 307         B1 = 0.5D+00*(ALIST(MAXERR)+BLIST(MAXERR))
 308         A2 = B1
 309         B2 = BLIST(MAXERR)
 310         ERLAST = ERRMAX
 311         CALL DQK15I(F,BOUN,INF,A1,B1,AREA1,ERROR1,RESABS,DEFAB1)
 312         CALL DQK15I(F,BOUN,INF,A2,B2,AREA2,ERROR2,RESABS,DEFAB2)
 313 C
 314 C           IMPROVE PREVIOUS APPROXIMATIONS TO INTEGRAL
 315 C           AND ERROR AND TEST FOR ACCURACY.
 316 C
 317         AREA12 = AREA1+AREA2
 318         ERRO12 = ERROR1+ERROR2
 319         ERRSUM = ERRSUM+ERRO12-ERRMAX
 320         AREA = AREA+AREA12-RLIST(MAXERR)
 321         IF(DEFAB1.EQ.ERROR1.OR.DEFAB2.EQ.ERROR2)GO TO 15
 322         IF(ABS(RLIST(MAXERR)-AREA12).GT.0.1D-04*ABS(AREA12)
 323      1  .OR.ERRO12.LT.0.99D+00*ERRMAX) GO TO 10
 324         IF(EXTRAP) IROFF2 = IROFF2+1
 325         IF(.NOT.EXTRAP) IROFF1 = IROFF1+1
 326    10   IF(LAST.GT.10.AND.ERRO12.GT.ERRMAX) IROFF3 = IROFF3+1
 327    15   RLIST(MAXERR) = AREA1
 328         RLIST(LAST) = AREA2
 329         ERRBND = MAX(EPSABS,EPSREL*ABS(AREA))
 330 C
 331 C           TEST FOR ROUNDOFF ERROR AND EVENTUALLY SET ERROR FLAG.
 332 C
 333         IF(IROFF1+IROFF2.GE.10.OR.IROFF3.GE.20) IER = 2
 334         IF(IROFF2.GE.5) IERRO = 3
 335 C
 336 C           SET ERROR FLAG IN THE CASE THAT THE NUMBER OF
 337 C           SUBINTERVALS EQUALS LIMIT.
 338 C
 339         IF(LAST.EQ.LIMIT) IER = 1
 340 C
 341 C           SET ERROR FLAG IN THE CASE OF BAD INTEGRAND BEHAVIOUR
 342 C           AT SOME POINTS OF THE INTEGRATION RANGE.
 343 C
 344         IF(MAX(ABS(A1),ABS(B2)).LE.(0.1D+01+0.1D+03*EPMACH)*
 345      1  (ABS(A2)+0.1D+04*UFLOW)) IER = 4
 346 C
 347 C           APPEND THE NEWLY-CREATED INTERVALS TO THE LIST.
 348 C
 349         IF(ERROR2.GT.ERROR1) GO TO 20
 350         ALIST(LAST) = A2
 351         BLIST(MAXERR) = B1
 352         BLIST(LAST) = B2
 353         ELIST(MAXERR) = ERROR1
 354         ELIST(LAST) = ERROR2
 355         GO TO 30
 356    20   ALIST(MAXERR) = A2
 357         ALIST(LAST) = A1
 358         BLIST(LAST) = B1
 359         RLIST(MAXERR) = AREA2
 360         RLIST(LAST) = AREA1
 361         ELIST(MAXERR) = ERROR2
 362         ELIST(LAST) = ERROR1
 363 C
 364 C           CALL SUBROUTINE DQPSRT TO MAINTAIN THE DESCENDING ORDERING
 365 C           IN THE LIST OF ERROR ESTIMATES AND SELECT THE SUBINTERVAL
 366 C           WITH NRMAX-TH LARGEST ERROR ESTIMATE (TO BE BISECTED NEXT).
 367 C
 368    30   CALL DQPSRT(LIMIT,LAST,MAXERR,ERRMAX,ELIST,IORD,NRMAX)
 369         IF(ERRSUM.LE.ERRBND) GO TO 115
 370         IF(IER.NE.0) GO TO 100
 371         IF(LAST.EQ.2) GO TO 80
 372         IF(NOEXT) GO TO 90
 373         ERLARG = ERLARG-ERLAST
 374         IF(ABS(B1-A1).GT.SMALL) ERLARG = ERLARG+ERRO12
 375         IF(EXTRAP) GO TO 40
 376 C
 377 C           TEST WHETHER THE INTERVAL TO BE BISECTED NEXT IS THE
 378 C           SMALLEST INTERVAL.
 379 C
 380         IF(ABS(BLIST(MAXERR)-ALIST(MAXERR)).GT.SMALL) GO TO 90
 381         EXTRAP = .TRUE.
 382         NRMAX = 2
 383    40   IF(IERRO.EQ.3.OR.ERLARG.LE.ERTEST) GO TO 60
 384 C
 385 C           THE SMALLEST INTERVAL HAS THE LARGEST ERROR.
 386 C           BEFORE BISECTING DECREASE THE SUM OF THE ERRORS OVER THE
 387 C           LARGER INTERVALS (ERLARG) AND PERFORM EXTRAPOLATION.
 388 C
 389         ID = NRMAX
 390         JUPBND = LAST
 391         IF(LAST.GT.(2+LIMIT/2)) JUPBND = LIMIT+3-LAST
 392         DO 50 K = ID,JUPBND
 393           MAXERR = IORD(NRMAX)
 394           ERRMAX = ELIST(MAXERR)
 395           IF(ABS(BLIST(MAXERR)-ALIST(MAXERR)).GT.SMALL) GO TO 90
 396           NRMAX = NRMAX+1
 397    50   CONTINUE
 398 C
 399 C           PERFORM EXTRAPOLATION.
 400 C
 401    60   NUMRL2 = NUMRL2+1
 402         RLIST2(NUMRL2) = AREA
 403         CALL DQELG(NUMRL2,RLIST2,RESEPS,ABSEPS,RES3LA,NRES)
 404         KTMIN = KTMIN+1
 405         IF(KTMIN.GT.5.AND.ABSERR.LT.0.1D-02*ERRSUM) IER = 5
 406         IF(ABSEPS.GE.ABSERR) GO TO 70
 407         KTMIN = 0
 408         ABSERR = ABSEPS
 409         RESULT = RESEPS
 410         CORREC = ERLARG
 411         ERTEST = MAX(EPSABS,EPSREL*ABS(RESEPS))
 412         IF(ABSERR.LE.ERTEST) GO TO 100
 413 C
 414 C            PREPARE BISECTION OF THE SMALLEST INTERVAL.
 415 C
 416    70   IF(NUMRL2.EQ.1) NOEXT = .TRUE.
 417         IF(IER.EQ.5) GO TO 100
 418         MAXERR = IORD(1)
 419         ERRMAX = ELIST(MAXERR)
 420         NRMAX = 1
 421         EXTRAP = .FALSE.
 422         SMALL = SMALL*0.5D+00
 423         ERLARG = ERRSUM
 424         GO TO 90
 425    80   SMALL = 0.375D+00
 426         ERLARG = ERRSUM
 427         ERTEST = ERRBND
 428         RLIST2(2) = AREA
 429    90 CONTINUE
 430 C
 431 C           SET FINAL RESULT AND ERROR ESTIMATE.
 432 C           ------------------------------------
 433 C
 434   100 IF(ABSERR.EQ.OFLOW) GO TO 115
 435       IF((IER+IERRO).EQ.0) GO TO 110
 436       IF(IERRO.EQ.3) ABSERR = ABSERR+CORREC
 437       IF(IER.EQ.0) IER = 3
 438       IF(RESULT.NE.0.0D+00.AND.AREA.NE.0.0D+00)GO TO 105
 439       IF(ABSERR.GT.ERRSUM)GO TO 115
 440       IF(AREA.EQ.0.0D+00) GO TO 130
 441       GO TO 110
 442   105 IF(ABSERR/ABS(RESULT).GT.ERRSUM/ABS(AREA))GO TO 115
 443 C
 444 C           TEST ON DIVERGENCE
 445 C
 446   110 IF(KSGN.EQ.(-1).AND.MAX(ABS(RESULT),ABS(AREA)).LE.
 447      1 DEFABS*0.1D-01) GO TO 130
 448       IF (0.1D-01.GT.(RESULT/AREA).OR.(RESULT/AREA).GT.0.1D+03
 449      1    .OR.ERRSUM.GT.ABS(AREA)) IER = 6
 450       GO TO 130
 451 C
 452 C           COMPUTE GLOBAL INTEGRAL SUM.
 453 C
 454   115 RESULT = 0.0D+00
 455       DO 120 K = 1,LAST
 456         RESULT = RESULT+RLIST(K)
 457   120 CONTINUE
 458       ABSERR = ERRSUM
 459   130 NEVAL = 30*LAST-15
 460       IF(INF.EQ.2) NEVAL = 2*NEVAL
 461       IF(IER.GT.2) IER=IER-1
 462   999 RETURN
 463       END