relibc_openlibm/slatec/dqk51.f

*DECK DQK51
      SUBROUTINE DQK51 (F, A, B, RESULT, ABSERR, RESABS, RESASC)
C***BEGIN PROLOGUE  DQK51
C***PURPOSE  To compute I = Integral of F over (A,B) with error
C                           estimate
C                       J = Integral of ABS(F) over (A,B)
C***LIBRARY   SLATEC (QUADPACK)
C***CATEGORY  H2A1A2
C***TYPE      DOUBLE PRECISION (QK51-S, DQK51-D)
C***KEYWORDS  51-POINT GAUSS-KRONROD RULES, QUADPACK, QUADRATURE
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           Integration rules
C           Standard fortran subroutine
C           Double precision version
C
C           PARAMETERS
C            ON ENTRY
C              F      - Double precision
C                       Function subroutine defining the integrand
C                       function F(X). The actual name for F needs to be
C                       declared E X T E R N A L in the calling program.
C
C              A      - Double precision
C                       Lower limit of integration
C
C              B      - Double precision
C                       Upper limit of integration
C
C            ON RETURN
C              RESULT - Double precision
C                       Approximation to the integral I
C                       RESULT is computed by applying the 51-point
C                       Kronrod rule (RESK) obtained by optimal addition
C                       of abscissae to the 25-point Gauss rule (RESG).
C
C              ABSERR - Double precision
C                       Estimate of the modulus of the absolute error,
C                       which should not exceed ABS(I-RESULT)
C
C              RESABS - Double precision
C                       Approximation to the integral J
C
C              RESASC - Double precision
C                       Approximation to the integral of ABS(F-I/(B-A))
C                       over (A,B)
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  D1MACH
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   910819  Added WGK(26) to code.  (WRB)
C***END PROLOGUE  DQK51
C
      DOUBLE PRECISION A,ABSC,ABSERR,B,CENTR,DHLGTH,
     1  D1MACH,EPMACH,F,FC,FSUM,FVAL1,FVAL2,FV1,FV2,HLGTH,RESABS,RESASC,
     2  RESG,RESK,RESKH,RESULT,UFLOW,WG,WGK,XGK
      INTEGER J,JTW,JTWM1
      EXTERNAL F
C
      DIMENSION FV1(25),FV2(25),XGK(26),WGK(26),WG(13)
C
C           THE ABSCISSAE AND WEIGHTS ARE GIVEN FOR THE INTERVAL (-1,1).
C           BECAUSE OF SYMMETRY ONLY THE POSITIVE ABSCISSAE AND THEIR
C           CORRESPONDING WEIGHTS ARE GIVEN.
C
C           XGK    - ABSCISSAE OF THE 51-POINT KRONROD RULE
C                    XGK(2), XGK(4), ...  ABSCISSAE OF THE 25-POINT
C                    GAUSS RULE
C                    XGK(1), XGK(3), ...  ABSCISSAE WHICH ARE OPTIMALLY
C                    ADDED TO THE 25-POINT GAUSS RULE
C
C           WGK    - WEIGHTS OF THE 51-POINT KRONROD RULE
C
C           WG     - WEIGHTS OF THE 25-POINT GAUSS RULE
C
C
C GAUSS QUADRATURE WEIGHTS AND KRONROD QUADRATURE ABSCISSAE AND WEIGHTS
C AS EVALUATED WITH 80 DECIMAL DIGIT ARITHMETIC BY L. W. FULLERTON,
C BELL LABS, NOV. 1981.
C
      SAVE WG, XGK, WGK
      DATA WG  (  1) / 0.0113937985 0102628794 7902964113 235 D0 /
      DATA WG  (  2) / 0.0263549866 1503213726 1901815295 299 D0 /
      DATA WG  (  3) / 0.0409391567 0130631265 5623487711 646 D0 /
      DATA WG  (  4) / 0.0549046959 7583519192 5936891540 473 D0 /
      DATA WG  (  5) / 0.0680383338 1235691720 7187185656 708 D0 /
      DATA WG  (  6) / 0.0801407003 3500101801 3234959669 111 D0 /
      DATA WG  (  7) / 0.0910282619 8296364981 1497220702 892 D0 /
      DATA WG  (  8) / 0.1005359490 6705064420 2206890392 686 D0 /
      DATA WG  (  9) / 0.1085196244 7426365311 6093957050 117 D0 /
      DATA WG  ( 10) / 0.1148582591 4571164833 9325545869 556 D0 /
      DATA WG  ( 11) / 0.1194557635 3578477222 8178126512 901 D0 /
      DATA WG  ( 12) / 0.1222424429 9031004168 8959518945 852 D0 /
      DATA WG  ( 13) / 0.1231760537 2671545120 3902873079 050 D0 /
C
      DATA XGK (  1) / 0.9992621049 9260983419 3457486540 341 D0 /
      DATA XGK (  2) / 0.9955569697 9049809790 8784946893 902 D0 /
      DATA XGK (  3) / 0.9880357945 3407724763 7331014577 406 D0 /
      DATA XGK (  4) / 0.9766639214 5951751149 8315386479 594 D0 /
      DATA XGK (  5) / 0.9616149864 2584251241 8130033660 167 D0 /
      DATA XGK (  6) / 0.9429745712 2897433941 4011169658 471 D0 /
      DATA XGK (  7) / 0.9207471152 8170156174 6346084546 331 D0 /
      DATA XGK (  8) / 0.8949919978 7827536885 1042006782 805 D0 /
      DATA XGK (  9) / 0.8658470652 9327559544 8996969588 340 D0 /
      DATA XGK ( 10) / 0.8334426287 6083400142 1021108693 570 D0 /
      DATA XGK ( 11) / 0.7978737979 9850005941 0410904994 307 D0 /
      DATA XGK ( 12) / 0.7592592630 3735763057 7282865204 361 D0 /
      DATA XGK ( 13) / 0.7177664068 1308438818 6654079773 298 D0 /
      DATA XGK ( 14) / 0.6735663684 7346836448 5120633247 622 D0 /
      DATA XGK ( 15) / 0.6268100990 1031741278 8122681624 518 D0 /
      DATA XGK ( 16) / 0.5776629302 4122296772 3689841612 654 D0 /
      DATA XGK ( 17) / 0.5263252843 3471918259 9623778158 010 D0 /
      DATA XGK ( 18) / 0.4730027314 4571496052 2182115009 192 D0 /
      DATA XGK ( 19) / 0.4178853821 9303774885 1814394594 572 D0 /
      DATA XGK ( 20) / 0.3611723058 0938783773 5821730127 641 D0 /
      DATA XGK ( 21) / 0.3030895389 3110783016 7478909980 339 D0 /
      DATA XGK ( 22) / 0.2438668837 2098843204 5190362797 452 D0 /
      DATA XGK ( 23) / 0.1837189394 2104889201 5969888759 528 D0 /
      DATA XGK ( 24) / 0.1228646926 1071039638 7359818808 037 D0 /
      DATA XGK ( 25) / 0.0615444830 0568507888 6546392366 797 D0 /
      DATA XGK ( 26) / 0.0000000000 0000000000 0000000000 000 D0 /
C
      DATA WGK (  1) / 0.0019873838 9233031592 6507851882 843 D0 /
      DATA WGK (  2) / 0.0055619321 3535671375 8040236901 066 D0 /
      DATA WGK (  3) / 0.0094739733 8617415160 7207710523 655 D0 /
      DATA WGK (  4) / 0.0132362291 9557167481 3656405846 976 D0 /
      DATA WGK (  5) / 0.0168478177 0912829823 1516667536 336 D0 /
      DATA WGK (  6) / 0.0204353711 4588283545 6568292235 939 D0 /
      DATA WGK (  7) / 0.0240099456 0695321622 0092489164 881 D0 /
      DATA WGK (  8) / 0.0274753175 8785173780 2948455517 811 D0 /
      DATA WGK (  9) / 0.0307923001 6738748889 1109020215 229 D0 /
      DATA WGK ( 10) / 0.0340021302 7432933783 6748795229 551 D0 /
      DATA WGK ( 11) / 0.0371162714 8341554356 0330625367 620 D0 /
      DATA WGK ( 12) / 0.0400838255 0403238207 4839284467 076 D0 /
      DATA WGK ( 13) / 0.0428728450 2017004947 6895792439 495 D0 /
      DATA WGK ( 14) / 0.0455029130 4992178890 9870584752 660 D0 /
      DATA WGK ( 15) / 0.0479825371 3883671390 6392255756 915 D0 /
      DATA WGK ( 16) / 0.0502776790 8071567196 3325259433 440 D0 /
      DATA WGK ( 17) / 0.0523628858 0640747586 4366712137 873 D0 /
      DATA WGK ( 18) / 0.0542511298 8854549014 4543370459 876 D0 /
      DATA WGK ( 19) / 0.0559508112 2041231730 8240686382 747 D0 /
      DATA WGK ( 20) / 0.0574371163 6156783285 3582693939 506 D0 /
      DATA WGK ( 21) / 0.0586896800 2239420796 1974175856 788 D0 /
      DATA WGK ( 22) / 0.0597203403 2417405997 9099291932 562 D0 /
      DATA WGK ( 23) / 0.0605394553 7604586294 5360267517 565 D0 /
      DATA WGK ( 24) / 0.0611285097 1705304830 5859030416 293 D0 /
      DATA WGK ( 25) / 0.0614711898 7142531666 1544131965 264 D0 /
      DATA WGK ( 26) / 0.0615808180 6783293507 8759824240 055 D0 /
C
C
C           LIST OF MAJOR VARIABLES
C           -----------------------
C
C           CENTR  - MID POINT OF THE INTERVAL
C           HLGTH  - HALF-LENGTH OF THE INTERVAL
C           ABSC   - ABSCISSA
C           FVAL*  - FUNCTION VALUE
C           RESG   - RESULT OF THE 25-POINT GAUSS FORMULA
C           RESK   - RESULT OF THE 51-POINT KRONROD FORMULA
C           RESKH  - APPROXIMATION TO THE MEAN VALUE OF F OVER (A,B),
C                    I.E. TO 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  DQK51
      EPMACH = D1MACH(4)
      UFLOW = D1MACH(1)
C
      CENTR = 0.5D+00*(A+B)
      HLGTH = 0.5D+00*(B-A)
      DHLGTH = ABS(HLGTH)
C
C           COMPUTE THE 51-POINT KRONROD APPROXIMATION TO
C           THE INTEGRAL, AND ESTIMATE THE ABSOLUTE ERROR.
C
      FC = F(CENTR)
      RESG = WG(13)*FC
      RESK = WGK(26)*FC
      RESABS = ABS(RESK)
      DO 10 J=1,12
        JTW = J*2
        ABSC = HLGTH*XGK(JTW)
        FVAL1 = F(CENTR-ABSC)
        FVAL2 = F(CENTR+ABSC)
        FV1(JTW) = FVAL1
        FV2(JTW) = FVAL2
        FSUM = FVAL1+FVAL2
        RESG = RESG+WG(J)*FSUM
        RESK = RESK+WGK(JTW)*FSUM
        RESABS = RESABS+WGK(JTW)*(ABS(FVAL1)+ABS(FVAL2))
   10 CONTINUE
      DO 15 J = 1,13
        JTWM1 = J*2-1
        ABSC = HLGTH*XGK(JTWM1)
        FVAL1 = F(CENTR-ABSC)
        FVAL2 = F(CENTR+ABSC)
        FV1(JTWM1) = FVAL1
        FV2(JTWM1) = FVAL2
        FSUM = FVAL1+FVAL2
        RESK = RESK+WGK(JTWM1)*FSUM
        RESABS = RESABS+WGK(JTWM1)*(ABS(FVAL1)+ABS(FVAL2))
   15 CONTINUE
      RESKH = RESK*0.5D+00
      RESASC = WGK(26)*ABS(FC-RESKH)
      DO 20 J=1,25
        RESASC = RESASC+WGK(J)*(ABS(FV1(J)-RESKH)+ABS(FV2(J)-RESKH))
   20 CONTINUE
      RESULT = RESK*HLGTH
      RESABS = RESABS*DHLGTH
      RESASC = RESASC*DHLGTH
      ABSERR = ABS((RESK-RESG)*HLGTH)
      IF(RESASC.NE.0.0D+00.AND.ABSERR.NE.0.0D+00)
     1  ABSERR = RESASC*MIN(0.1D+01,(0.2D+03*ABSERR/RESASC)**1.5D+00)
      IF(RESABS.GT.UFLOW/(0.5D+02*EPMACH)) ABSERR = MAX
     1  ((EPMACH*0.5D+02)*RESABS,ABSERR)
      RETURN
      END