src/numerical/slatec/fortran/dqk51.f

   1 *DECK DQK51
   2       SUBROUTINE DQK51 (F, A, B, RESULT, ABSERR, RESABS, RESASC)
   3 C***BEGIN PROLOGUE  DQK51
   4 C***PURPOSE  To compute I = Integral of F over (A,B) with error
   5 C                           estimate
   6 C                       J = Integral of ABS(F) over (A,B)
   7 C***LIBRARY   SLATEC (QUADPACK)
   8 C***CATEGORY  H2A1A2
   9 C***TYPE      DOUBLE PRECISION (QK51-S, DQK51-D)
  10 C***KEYWORDS  51-POINT GAUSS-KRONROD RULES, QUADPACK, QUADRATURE
  11 C***AUTHOR  Piessens, Robert
  12 C             Applied Mathematics and Programming Division
  13 C             K. U. Leuven
  14 C           de Doncker, Elise
  15 C             Applied Mathematics and Programming Division
  16 C             K. U. Leuven
  17 C***DESCRIPTION
  18 C
  19 C           Integration rules
  20 C           Standard fortran subroutine
  21 C           Double precision version
  22 C
  23 C           PARAMETERS
  24 C            ON ENTRY
  25 C              F      - Double precision
  26 C                       Function subroutine defining the integrand
  27 C                       function F(X). The actual name for F needs to be
  28 C                       declared E X T E R N A L in the calling program.
  29 C
  30 C              A      - Double precision
  31 C                       Lower limit of integration
  32 C
  33 C              B      - Double precision
  34 C                       Upper limit of integration
  35 C
  36 C            ON RETURN
  37 C              RESULT - Double precision
  38 C                       Approximation to the integral I
  39 C                       RESULT is computed by applying the 51-point
  40 C                       Kronrod rule (RESK) obtained by optimal addition
  41 C                       of abscissae to the 25-point Gauss rule (RESG).
  42 C
  43 C              ABSERR - Double precision
  44 C                       Estimate of the modulus of the absolute error,
  45 C                       which should not exceed ABS(I-RESULT)
  46 C
  47 C              RESABS - Double precision
  48 C                       Approximation to the integral J
  49 C
  50 C              RESASC - Double precision
  51 C                       Approximation to the integral of ABS(F-I/(B-A))
  52 C                       over (A,B)
  53 C
  54 C***REFERENCES  (NONE)
  55 C***ROUTINES CALLED  D1MACH
  56 C***REVISION HISTORY  (YYMMDD)
  57 C   800101  DATE WRITTEN
  58 C   890531  Changed all specific intrinsics to generic.  (WRB)
  59 C   890531  REVISION DATE from Version 3.2
  60 C   891214  Prologue converted to Version 4.0 format.  (BAB)
  61 C   910819  Added WGK(26) to code.  (WRB)
  62 C***END PROLOGUE  DQK51
  63 C
  64       DOUBLE PRECISION A,ABSC,ABSERR,B,CENTR,DHLGTH,
  65      1  D1MACH,EPMACH,F,FC,FSUM,FVAL1,FVAL2,FV1,FV2,HLGTH,RESABS,RESASC,
  66      2  RESG,RESK,RESKH,RESULT,UFLOW,WG,WGK,XGK
  67       INTEGER J,JTW,JTWM1
  68       EXTERNAL F
  69 C
  70       DIMENSION FV1(25),FV2(25),XGK(26),WGK(26),WG(13)
  71 C
  72 C           THE ABSCISSAE AND WEIGHTS ARE GIVEN FOR THE INTERVAL (-1,1).
  73 C           BECAUSE OF SYMMETRY ONLY THE POSITIVE ABSCISSAE AND THEIR
  74 C           CORRESPONDING WEIGHTS ARE GIVEN.
  75 C
  76 C           XGK    - ABSCISSAE OF THE 51-POINT KRONROD RULE
  77 C                    XGK(2), XGK(4), ...  ABSCISSAE OF THE 25-POINT
  78 C                    GAUSS RULE
  79 C                    XGK(1), XGK(3), ...  ABSCISSAE WHICH ARE OPTIMALLY
  80 C                    ADDED TO THE 25-POINT GAUSS RULE
  81 C
  82 C           WGK    - WEIGHTS OF THE 51-POINT KRONROD RULE
  83 C
  84 C           WG     - WEIGHTS OF THE 25-POINT GAUSS RULE
  85 C
  86 C
  87 C GAUSS QUADRATURE WEIGHTS AND KRONROD QUADRATURE ABSCISSAE AND WEIGHTS
  88 C AS EVALUATED WITH 80 DECIMAL DIGIT ARITHMETIC BY L. W. FULLERTON,
  89 C BELL LABS, NOV. 1981.
  90 C
  91       SAVE WG, XGK, WGK
  92       DATA WG  (  1) / 0.0113937985 0102628794 7902964113 235 D0 /
  93       DATA WG  (  2) / 0.0263549866 1503213726 1901815295 299 D0 /
  94       DATA WG  (  3) / 0.0409391567 0130631265 5623487711 646 D0 /
  95       DATA WG  (  4) / 0.0549046959 7583519192 5936891540 473 D0 /
  96       DATA WG  (  5) / 0.0680383338 1235691720 7187185656 708 D0 /
  97       DATA WG  (  6) / 0.0801407003 3500101801 3234959669 111 D0 /
  98       DATA WG  (  7) / 0.0910282619 8296364981 1497220702 892 D0 /
  99       DATA WG  (  8) / 0.1005359490 6705064420 2206890392 686 D0 /
 100       DATA WG  (  9) / 0.1085196244 7426365311 6093957050 117 D0 /
 101       DATA WG  ( 10) / 0.1148582591 4571164833 9325545869 556 D0 /
 102       DATA WG  ( 11) / 0.1194557635 3578477222 8178126512 901 D0 /
 103       DATA WG  ( 12) / 0.1222424429 9031004168 8959518945 852 D0 /
 104       DATA WG  ( 13) / 0.1231760537 2671545120 3902873079 050 D0 /
 105 C
 106       DATA XGK (  1) / 0.9992621049 9260983419 3457486540 341 D0 /
 107       DATA XGK (  2) / 0.9955569697 9049809790 8784946893 902 D0 /
 108       DATA XGK (  3) / 0.9880357945 3407724763 7331014577 406 D0 /
 109       DATA XGK (  4) / 0.9766639214 5951751149 8315386479 594 D0 /
 110       DATA XGK (  5) / 0.9616149864 2584251241 8130033660 167 D0 /
 111       DATA XGK (  6) / 0.9429745712 2897433941 4011169658 471 D0 /
 112       DATA XGK (  7) / 0.9207471152 8170156174 6346084546 331 D0 /
 113       DATA XGK (  8) / 0.8949919978 7827536885 1042006782 805 D0 /
 114       DATA XGK (  9) / 0.8658470652 9327559544 8996969588 340 D0 /
 115       DATA XGK ( 10) / 0.8334426287 6083400142 1021108693 570 D0 /
 116       DATA XGK ( 11) / 0.7978737979 9850005941 0410904994 307 D0 /
 117       DATA XGK ( 12) / 0.7592592630 3735763057 7282865204 361 D0 /
 118       DATA XGK ( 13) / 0.7177664068 1308438818 6654079773 298 D0 /
 119       DATA XGK ( 14) / 0.6735663684 7346836448 5120633247 622 D0 /
 120       DATA XGK ( 15) / 0.6268100990 1031741278 8122681624 518 D0 /
 121       DATA XGK ( 16) / 0.5776629302 4122296772 3689841612 654 D0 /
 122       DATA XGK ( 17) / 0.5263252843 3471918259 9623778158 010 D0 /
 123       DATA XGK ( 18) / 0.4730027314 4571496052 2182115009 192 D0 /
 124       DATA XGK ( 19) / 0.4178853821 9303774885 1814394594 572 D0 /
 125       DATA XGK ( 20) / 0.3611723058 0938783773 5821730127 641 D0 /
 126       DATA XGK ( 21) / 0.3030895389 3110783016 7478909980 339 D0 /
 127       DATA XGK ( 22) / 0.2438668837 2098843204 5190362797 452 D0 /
 128       DATA XGK ( 23) / 0.1837189394 2104889201 5969888759 528 D0 /
 129       DATA XGK ( 24) / 0.1228646926 1071039638 7359818808 037 D0 /
 130       DATA XGK ( 25) / 0.0615444830 0568507888 6546392366 797 D0 /
 131       DATA XGK ( 26) / 0.0000000000 0000000000 0000000000 000 D0 /
 132 C
 133       DATA WGK (  1) / 0.0019873838 9233031592 6507851882 843 D0 /
 134       DATA WGK (  2) / 0.0055619321 3535671375 8040236901 066 D0 /
 135       DATA WGK (  3) / 0.0094739733 8617415160 7207710523 655 D0 /
 136       DATA WGK (  4) / 0.0132362291 9557167481 3656405846 976 D0 /
 137       DATA WGK (  5) / 0.0168478177 0912829823 1516667536 336 D0 /
 138       DATA WGK (  6) / 0.0204353711 4588283545 6568292235 939 D0 /
 139       DATA WGK (  7) / 0.0240099456 0695321622 0092489164 881 D0 /
 140       DATA WGK (  8) / 0.0274753175 8785173780 2948455517 811 D0 /
 141       DATA WGK (  9) / 0.0307923001 6738748889 1109020215 229 D0 /
 142       DATA WGK ( 10) / 0.0340021302 7432933783 6748795229 551 D0 /
 143       DATA WGK ( 11) / 0.0371162714 8341554356 0330625367 620 D0 /
 144       DATA WGK ( 12) / 0.0400838255 0403238207 4839284467 076 D0 /
 145       DATA WGK ( 13) / 0.0428728450 2017004947 6895792439 495 D0 /
 146       DATA WGK ( 14) / 0.0455029130 4992178890 9870584752 660 D0 /
 147       DATA WGK ( 15) / 0.0479825371 3883671390 6392255756 915 D0 /
 148       DATA WGK ( 16) / 0.0502776790 8071567196 3325259433 440 D0 /
 149       DATA WGK ( 17) / 0.0523628858 0640747586 4366712137 873 D0 /
 150       DATA WGK ( 18) / 0.0542511298 8854549014 4543370459 876 D0 /
 151       DATA WGK ( 19) / 0.0559508112 2041231730 8240686382 747 D0 /
 152       DATA WGK ( 20) / 0.0574371163 6156783285 3582693939 506 D0 /
 153       DATA WGK ( 21) / 0.0586896800 2239420796 1974175856 788 D0 /
 154       DATA WGK ( 22) / 0.0597203403 2417405997 9099291932 562 D0 /
 155       DATA WGK ( 23) / 0.0605394553 7604586294 5360267517 565 D0 /
 156       DATA WGK ( 24) / 0.0611285097 1705304830 5859030416 293 D0 /
 157       DATA WGK ( 25) / 0.0614711898 7142531666 1544131965 264 D0 /
 158       DATA WGK ( 26) / 0.0615808180 6783293507 8759824240 055 D0 /
 159 C
 160 C
 161 C           LIST OF MAJOR VARIABLES
 162 C           -----------------------
 163 C
 164 C           CENTR  - MID POINT OF THE INTERVAL
 165 C           HLGTH  - HALF-LENGTH OF THE INTERVAL
 166 C           ABSC   - ABSCISSA
 167 C           FVAL*  - FUNCTION VALUE
 168 C           RESG   - RESULT OF THE 25-POINT GAUSS FORMULA
 169 C           RESK   - RESULT OF THE 51-POINT KRONROD FORMULA
 170 C           RESKH  - APPROXIMATION TO THE MEAN VALUE OF F OVER (A,B),
 171 C                    I.E. TO I/(B-A)
 172 C
 173 C           MACHINE DEPENDENT CONSTANTS
 174 C           ---------------------------
 175 C
 176 C           EPMACH IS THE LARGEST RELATIVE SPACING.
 177 C           UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
 178 C
 179 C***FIRST EXECUTABLE STATEMENT  DQK51
 180       EPMACH = D1MACH(4)
 181       UFLOW = D1MACH(1)
 182 C
 183       CENTR = 0.5D+00*(A+B)
 184       HLGTH = 0.5D+00*(B-A)
 185       DHLGTH = ABS(HLGTH)
 186 C
 187 C           COMPUTE THE 51-POINT KRONROD APPROXIMATION TO
 188 C           THE INTEGRAL, AND ESTIMATE THE ABSOLUTE ERROR.
 189 C
 190       FC = F(CENTR)
 191       RESG = WG(13)*FC
 192       RESK = WGK(26)*FC
 193       RESABS = ABS(RESK)
 194       DO 10 J=1,12
 195         JTW = J*2
 196         ABSC = HLGTH*XGK(JTW)
 197         FVAL1 = F(CENTR-ABSC)
 198         FVAL2 = F(CENTR+ABSC)
 199         FV1(JTW) = FVAL1
 200         FV2(JTW) = FVAL2
 201         FSUM = FVAL1+FVAL2
 202         RESG = RESG+WG(J)*FSUM
 203         RESK = RESK+WGK(JTW)*FSUM
 204         RESABS = RESABS+WGK(JTW)*(ABS(FVAL1)+ABS(FVAL2))
 205    10 CONTINUE
 206       DO 15 J = 1,13
 207         JTWM1 = J*2-1
 208         ABSC = HLGTH*XGK(JTWM1)
 209         FVAL1 = F(CENTR-ABSC)
 210         FVAL2 = F(CENTR+ABSC)
 211         FV1(JTWM1) = FVAL1
 212         FV2(JTWM1) = FVAL2
 213         FSUM = FVAL1+FVAL2
 214         RESK = RESK+WGK(JTWM1)*FSUM
 215         RESABS = RESABS+WGK(JTWM1)*(ABS(FVAL1)+ABS(FVAL2))
 216    15 CONTINUE
 217       RESKH = RESK*0.5D+00
 218       RESASC = WGK(26)*ABS(FC-RESKH)
 219       DO 20 J=1,25
 220         RESASC = RESASC+WGK(J)*(ABS(FV1(J)-RESKH)+ABS(FV2(J)-RESKH))
 221    20 CONTINUE
 222       RESULT = RESK*HLGTH
 223       RESABS = RESABS*DHLGTH
 224       RESASC = RESASC*DHLGTH
 225       ABSERR = ABS((RESK-RESG)*HLGTH)
 226       IF(RESASC.NE.0.0D+00.AND.ABSERR.NE.0.0D+00)
 227      1  ABSERR = RESASC*MIN(0.1D+01,(0.2D+03*ABSERR/RESASC)**1.5D+00)
 228       IF(RESABS.GT.UFLOW/(0.5D+02*EPMACH)) ABSERR = MAX
 229      1  ((EPMACH*0.5D+02)*RESABS,ABSERR)
 230       RETURN
 231       END