relibc_openlibm/slatec/qng.f

*DECK QNG
      SUBROUTINE QNG (F, A, B, EPSABS, EPSREL, RESULT, ABSERR, NEVAL,
     +   IER)
C***BEGIN PROLOGUE  QNG
C***PURPOSE  The routine calculates an approximation result to a
C            given definite integral I = integral of F over (A,B),
C            hopefully satisfying following claim for accuracy
C            ABS(I-RESULT).LE.MAX(EPSABS,EPSREL*ABS(I)).
C***LIBRARY   SLATEC (QUADPACK)
C***CATEGORY  H2A1A1
C***TYPE      SINGLE PRECISION (QNG-S, DQNG-D)
C***KEYWORDS  AUTOMATIC INTEGRATOR, GAUSS-KRONROD(PATTERSON) RULES,
C             NONADAPTIVE, QUADPACK, QUADRATURE, SMOOTH INTEGRAND
C***AUTHOR  Piessens, Robert
C             Applied Mathematics and Programming Division
C             K. U. Leuven
C           de Doncker, Elise
C             Applied Mathematics and Programming Division
C             K. U. Leuven
C***DESCRIPTION
C
C NON-ADAPTIVE INTEGRATION
C STANDARD FORTRAN SUBROUTINE
C REAL VERSION
C
C           F      - Real version
C                    Function subprogram defining the integrand function
C                    F(X). The actual name for F needs to be declared
C                    E X T E R N A L in the driver program.
C
C           A      - Real version
C                    Lower limit of integration
C
C           B      - Real version
C                    Upper limit of integration
C
C           EPSABS - Real
C                    Absolute accuracy requested
C           EPSREL - Real
C                    Relative accuracy requested
C                    If  EPSABS.LE.0
C                    And EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28),
C                    The routine will end with IER = 6.
C
C         ON RETURN
C           RESULT - Real
C                    Approximation to the integral I
C                    Result is obtained by applying the 21-POINT
C                    GAUSS-KRONROD RULE (RES21) obtained by optimal
C                    addition of abscissae to the 10-POINT GAUSS RULE
C                    (RES10), or by applying the 43-POINT RULE (RES43)
C                    obtained by optimal addition of abscissae to the
C                    21-POINT GAUSS-KRONROD RULE, or by applying the
C                    87-POINT RULE (RES87) obtained by optimal addition
C                    of abscissae to the 43-POINT RULE.
C
C           ABSERR - Real
C                    Estimate of the modulus of the absolute error,
C                    which should EQUAL or EXCEED ABS(I-RESULT)
C
C           NEVAL  - Integer
C                    Number of integrand evaluations
C
C           IER    - IER = 0 normal and reliable termination of the
C                            routine. It is assumed that the requested
C                            accuracy has been achieved.
C                    IER.GT.0 Abnormal termination of the routine. It is
C                            assumed that the requested accuracy has
C                            not been achieved.
C           ERROR MESSAGES
C                    IER = 1 The maximum number of steps has been
C                            executed. The integral is probably too
C                            difficult to be calculated by DQNG.
C                        = 6 The input is invalid, because
C                            EPSABS.LE.0 AND
C                            EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28).
C                            RESULT, ABSERR and NEVAL are set to zero.
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  R1MACH, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   800101  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C***END PROLOGUE  QNG
C
      REAL A,ABSC,ABSERR,B,CENTR,DHLGTH,EPMACH,EPSABS,EPSREL,F,FCENTR,
     1  FVAL,FVAL1,FVAL2,FV1,FV2,FV3,FV4,HLGTH,RESULT,RES10,RES21,RES43,
     2  RES87,RESABS,RESASC,RESKH,R1MACH,SAVFUN,UFLOW,W10,W21A,W43A,
     3  W43B,W87A,W87B,X1,X2,X3,X4
      INTEGER IER,IPX,K,L,NEVAL
      EXTERNAL F
C
      DIMENSION FV1(5),FV2(5),FV3(5),FV4(5),X1(5),X2(5),X3(11),X4(22),
     1  W10(5),W21A(5),W21B(6),W43A(10),W43B(12),W87A(21),W87B(23),
     2  SAVFUN(21)
C
C           THE FOLLOWING DATA STATEMENTS CONTAIN THE
C           ABSCISSAE AND WEIGHTS OF THE INTEGRATION RULES USED.
C
C           X1      ABSCISSAE COMMON TO THE 10-, 21-, 43-
C                   AND 87-POINT RULE
C           X2      ABSCISSAE COMMON TO THE 21-, 43- AND
C                   87-POINT RULE
C           X3      ABSCISSAE COMMON TO THE 43- AND 87-POINT
C                   RULE
C           X4      ABSCISSAE OF THE 87-POINT RULE
C           W10     WEIGHTS OF THE 10-POINT FORMULA
C           W21A    WEIGHTS OF THE 21-POINT FORMULA FOR
C                   ABSCISSAE X1
C           W21B    WEIGHTS OF THE 21-POINT FORMULA FOR
C                   ABSCISSAE X2
C           W43A    WEIGHTS OF THE 43-POINT FORMULA FOR
C                   ABSCISSAE X1, X3
C           W43B    WEIGHTS OF THE 43-POINT FORMULA FOR
C                   ABSCISSAE X3
C           W87A    WEIGHTS OF THE 87-POINT FORMULA FOR
C                   ABSCISSAE X1, X2, X3
C           W87B    WEIGHTS OF THE 87-POINT FORMULA FOR
C                   ABSCISSAE X4
C
      SAVE X1, X2, X3, X4, W10, W21A, W21B, W43A, W43B, W87A, W87B
      DATA X1(1),X1(2),X1(3),X1(4),X1(5)/
     1     0.9739065285171717E+00,     0.8650633666889845E+00,
     2     0.6794095682990244E+00,     0.4333953941292472E+00,
     3     0.1488743389816312E+00/
      DATA X2(1),X2(2),X2(3),X2(4),X2(5)/
     1     0.9956571630258081E+00,     0.9301574913557082E+00,
     2     0.7808177265864169E+00,     0.5627571346686047E+00,
     3     0.2943928627014602E+00/
      DATA X3(1),X3(2),X3(3),X3(4),X3(5),X3(6),X3(7),X3(8),
     1  X3(9),X3(10),X3(11)/
     2     0.9993333609019321E+00,     0.9874334029080889E+00,
     3     0.9548079348142663E+00,     0.9001486957483283E+00,
     4     0.8251983149831142E+00,     0.7321483889893050E+00,
     5     0.6228479705377252E+00,     0.4994795740710565E+00,
     6     0.3649016613465808E+00,     0.2222549197766013E+00,
     7     0.7465061746138332E-01/
      DATA X4(1),X4(2),X4(3),X4(4),X4(5),X4(6),X4(7),X4(8),X4(9),
     1  X4(10),X4(11),X4(12),X4(13),X4(14),X4(15),X4(16),X4(17),X4(18),
     2  X4(19),X4(20),X4(21),X4(22)/   0.9999029772627292E+00,
     3     0.9979898959866787E+00,     0.9921754978606872E+00,
     4     0.9813581635727128E+00,     0.9650576238583846E+00,
     5     0.9431676131336706E+00,     0.9158064146855072E+00,
     6     0.8832216577713165E+00,     0.8457107484624157E+00,
     7     0.8035576580352310E+00,     0.7570057306854956E+00,
     8     0.7062732097873218E+00,     0.6515894665011779E+00,
     9     0.5932233740579611E+00,     0.5314936059708319E+00,
     1     0.4667636230420228E+00,     0.3994248478592188E+00,
     2     0.3298748771061883E+00,     0.2585035592021616E+00,
     3     0.1856953965683467E+00,     0.1118422131799075E+00,
     4     0.3735212339461987E-01/
      DATA W10(1),W10(2),W10(3),W10(4),W10(5)/
     1     0.6667134430868814E-01,     0.1494513491505806E+00,
     2     0.2190863625159820E+00,     0.2692667193099964E+00,
     3     0.2955242247147529E+00/
      DATA W21A(1),W21A(2),W21A(3),W21A(4),W21A(5)/
     1     0.3255816230796473E-01,     0.7503967481091995E-01,
     2     0.1093871588022976E+00,     0.1347092173114733E+00,
     3     0.1477391049013385E+00/
      DATA W21B(1),W21B(2),W21B(3),W21B(4),W21B(5),W21B(6)/
     1     0.1169463886737187E-01,     0.5475589657435200E-01,
     2     0.9312545458369761E-01,     0.1234919762620659E+00,
     3     0.1427759385770601E+00,     0.1494455540029169E+00/
      DATA W43A(1),W43A(2),W43A(3),W43A(4),W43A(5),W43A(6),W43A(7),
     1  W43A(8),W43A(9),W43A(10)/      0.1629673428966656E-01,
     2     0.3752287612086950E-01,     0.5469490205825544E-01,
     3     0.6735541460947809E-01,     0.7387019963239395E-01,
     4     0.5768556059769796E-02,     0.2737189059324884E-01,
     5     0.4656082691042883E-01,     0.6174499520144256E-01,
     6     0.7138726726869340E-01/
      DATA W43B(1),W43B(2),W43B(3),W43B(4),W43B(5),W43B(6),
     1  W43B(7),W43B(8),W43B(9),W43B(10),W43B(11),W43B(12)/
     2     0.1844477640212414E-02,     0.1079868958589165E-01,
     3     0.2189536386779543E-01,     0.3259746397534569E-01,
     4     0.4216313793519181E-01,     0.5074193960018458E-01,
     5     0.5837939554261925E-01,     0.6474640495144589E-01,
     6     0.6956619791235648E-01,     0.7282444147183321E-01,
     7     0.7450775101417512E-01,     0.7472214751740301E-01/
      DATA W87A(1),W87A(2),W87A(3),W87A(4),W87A(5),W87A(6),
     1  W87A(7),W87A(8),W87A(9),W87A(10),W87A(11),W87A(12),
     2  W87A(13),W87A(14),W87A(15),W87A(16),W87A(17),W87A(18),
     3  W87A(19),W87A(20),W87A(21)/
     4     0.8148377384149173E-02,     0.1876143820156282E-01,
     5     0.2734745105005229E-01,     0.3367770731163793E-01,
     6     0.3693509982042791E-01,     0.2884872430211531E-02,
     7     0.1368594602271270E-01,     0.2328041350288831E-01,
     8     0.3087249761171336E-01,     0.3569363363941877E-01,
     9     0.9152833452022414E-03,     0.5399280219300471E-02,
     1     0.1094767960111893E-01,     0.1629873169678734E-01,
     2     0.2108156888920384E-01,     0.2537096976925383E-01,
     3     0.2918969775647575E-01,     0.3237320246720279E-01,
     4     0.3478309895036514E-01,     0.3641222073135179E-01,
     5     0.3725387550304771E-01/
      DATA W87B(1),W87B(2),W87B(3),W87B(4),W87B(5),W87B(6),W87B(7),
     1  W87B(8),W87B(9),W87B(10),W87B(11),W87B(12),W87B(13),W87B(14),
     2  W87B(15),W87B(16),W87B(17),W87B(18),W87B(19),W87B(20),
     3  W87B(21),W87B(22),W87B(23)/    0.2741455637620724E-03,
     4     0.1807124155057943E-02,     0.4096869282759165E-02,
     5     0.6758290051847379E-02,     0.9549957672201647E-02,
     6     0.1232944765224485E-01,     0.1501044734638895E-01,
     7     0.1754896798624319E-01,     0.1993803778644089E-01,
     8     0.2219493596101229E-01,     0.2433914712600081E-01,
     9     0.2637450541483921E-01,     0.2828691078877120E-01,
     1     0.3005258112809270E-01,     0.3164675137143993E-01,
     2     0.3305041341997850E-01,     0.3425509970422606E-01,
     3     0.3526241266015668E-01,     0.3607698962288870E-01,
     4     0.3669860449845609E-01,     0.3712054926983258E-01,
     5     0.3733422875193504E-01,     0.3736107376267902E-01/
C
C           LIST OF MAJOR VARIABLES
C           -----------------------
C
C           CENTR  - MID POINT OF THE INTEGRATION INTERVAL
C           HLGTH  - HALF-LENGTH OF THE INTEGRATION INTERVAL
C           FCENTR - FUNCTION VALUE AT MID POINT
C           ABSC   - ABSCISSA
C           FVAL   - FUNCTION VALUE
C           SAVFUN - ARRAY OF FUNCTION VALUES WHICH
C                    HAVE ALREADY BEEN COMPUTED
C           RES10  - 10-POINT GAUSS RESULT
C           RES21  - 21-POINT KRONROD RESULT
C           RES43  - 43-POINT RESULT
C           RES87  - 87-POINT RESULT
C           RESABS - APPROXIMATION TO THE INTEGRAL OF ABS(F)
C           RESASC - APPROXIMATION TO THE INTEGRAL OF ABS(F-I/(B-A))
C
C           MACHINE DEPENDENT CONSTANTS
C           ---------------------------
C
C           EPMACH IS THE LARGEST RELATIVE SPACING.
C           UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
C
C***FIRST EXECUTABLE STATEMENT  QNG
      EPMACH = R1MACH(4)
      UFLOW = R1MACH(1)
C
C           TEST ON VALIDITY OF PARAMETERS
C           ------------------------------
C
      RESULT = 0.0E+00
      ABSERR = 0.0E+00
      NEVAL = 0
      IER = 6
      IF(EPSABS.LE.0.0E+00.AND.EPSREL.LT.MAX(0.5E-14,0.5E+02*EPMACH))
     1  GO TO 80
      HLGTH = 0.5E+00*(B-A)
      DHLGTH = ABS(HLGTH)
      CENTR = 0.5E+00*(B+A)
      FCENTR = F(CENTR)
      NEVAL = 21
      IER = 1
C
C          COMPUTE THE INTEGRAL USING THE 10- AND 21-POINT FORMULA.
C
      DO 70 L = 1,3
      GO TO (5,25,45),L
    5 RES10 = 0.0E+00
      RES21 = W21B(6)*FCENTR
      RESABS = W21B(6)*ABS(FCENTR)
      DO 10 K=1,5
        ABSC = HLGTH*X1(K)
        FVAL1 = F(CENTR+ABSC)
        FVAL2 = F(CENTR-ABSC)
        FVAL = FVAL1+FVAL2
        RES10 = RES10+W10(K)*FVAL
        RES21 = RES21+W21A(K)*FVAL
        RESABS = RESABS+W21A(K)*(ABS(FVAL1)+ABS(FVAL2))
        SAVFUN(K) = FVAL
        FV1(K) = FVAL1
        FV2(K) = FVAL2
   10 CONTINUE
      IPX = 5
      DO 15 K=1,5
        IPX = IPX+1
        ABSC = HLGTH*X2(K)
        FVAL1 = F(CENTR+ABSC)
        FVAL2 = F(CENTR-ABSC)
        FVAL = FVAL1+FVAL2
        RES21 = RES21+W21B(K)*FVAL
        RESABS = RESABS+W21B(K)*(ABS(FVAL1)+ABS(FVAL2))
        SAVFUN(IPX) = FVAL
        FV3(K) = FVAL1
        FV4(K) = FVAL2
   15 CONTINUE
C
C          TEST FOR CONVERGENCE.
C
      RESULT = RES21*HLGTH
      RESABS = RESABS*DHLGTH
      RESKH = 0.5E+00*RES21
      RESASC = W21B(6)*ABS(FCENTR-RESKH)
      DO 20 K = 1,5
        RESASC = RESASC+W21A(K)*(ABS(FV1(K)-RESKH)+ABS(FV2(K)-RESKH))
     1                  +W21B(K)*(ABS(FV3(K)-RESKH)+ABS(FV4(K)-RESKH))
   20 CONTINUE
      ABSERR = ABS((RES21-RES10)*HLGTH)
      RESASC = RESASC*DHLGTH
      GO TO 65
C
C          COMPUTE THE INTEGRAL USING THE 43-POINT FORMULA.
C
   25 RES43 = W43B(12)*FCENTR
      NEVAL = 43
      DO 30 K=1,10
        RES43 = RES43+SAVFUN(K)*W43A(K)
   30 CONTINUE
      DO 40 K=1,11
        IPX = IPX+1
        ABSC = HLGTH*X3(K)
        FVAL = F(ABSC+CENTR)+F(CENTR-ABSC)
        RES43 = RES43+FVAL*W43B(K)
        SAVFUN(IPX) = FVAL
   40 CONTINUE
C
C          TEST FOR CONVERGENCE.
C
      RESULT = RES43*HLGTH
      ABSERR = ABS((RES43-RES21)*HLGTH)
      GO TO 65
C
C          COMPUTE THE INTEGRAL USING THE 87-POINT FORMULA.
C
   45 RES87 = W87B(23)*FCENTR
      NEVAL = 87
      DO 50 K=1,21
        RES87 = RES87+SAVFUN(K)*W87A(K)
   50 CONTINUE
      DO 60 K=1,22
        ABSC = HLGTH*X4(K)
        RES87 = RES87+W87B(K)*(F(ABSC+CENTR)+F(CENTR-ABSC))
   60 CONTINUE
      RESULT = RES87*HLGTH
      ABSERR = ABS((RES87-RES43)*HLGTH)
   65 IF(RESASC.NE.0.0E+00.AND.ABSERR.NE.0.0E+00)
     1  ABSERR = RESASC*MIN(0.1E+01,
     2  (0.2E+03*ABSERR/RESASC)**1.5E+00)
      IF (RESABS.GT.UFLOW/(0.5E+02*EPMACH)) ABSERR = MAX
     1  ((EPMACH*0.5E+02)*RESABS,ABSERR)
      IF (ABSERR.LE.MAX(EPSABS,EPSREL*ABS(RESULT))) IER = 0
C ***JUMP OUT OF DO-LOOP
      IF (IER.EQ.0) GO TO 999
   70 CONTINUE
   80 CALL XERMSG ('SLATEC', 'QNG', 'ABNORMAL RETURN', IER, 0)
  999 RETURN
      END