src/numerical/slatec/fortran/dqk31.f

   1 *DECK DQK31
   2       SUBROUTINE DQK31 (F, A, B, RESULT, ABSERR, RESABS, RESASC)
   3 C***BEGIN PROLOGUE  DQK31
   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 (QK31-S, DQK31-D)
  10 C***KEYWORDS  31-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 subprogram 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 31-POINT
  40 C                       GAUSS-KRONROD RULE (RESK), obtained by optimal
  41 C                       addition of abscissae to the 15-POINT GAUSS
  42 C                       RULE (RESG).
  43 C
  44 C              ABSERR - Double precision
  45 C                       Estimate of the modulus of the modulus,
  46 C                       which should not exceed ABS(I-RESULT)
  47 C
  48 C              RESABS - Double precision
  49 C                       Approximation to the integral J
  50 C
  51 C              RESASC - Double precision
  52 C                       Approximation to the integral of ABS(F-I/(B-A))
  53 C                       over (A,B)
  54 C
  55 C***REFERENCES  (NONE)
  56 C***ROUTINES CALLED  D1MACH
  57 C***REVISION HISTORY  (YYMMDD)
  58 C   800101  DATE WRITTEN
  59 C   890531  Changed all specific intrinsics to generic.  (WRB)
  60 C   890531  REVISION DATE from Version 3.2
  61 C   891214  Prologue converted to Version 4.0 format.  (BAB)
  62 C***END PROLOGUE  DQK31
  63       DOUBLE PRECISION A,ABSC,ABSERR,B,CENTR,DHLGTH,
  64      1  D1MACH,EPMACH,F,FC,FSUM,FVAL1,FVAL2,FV1,FV2,HLGTH,RESABS,RESASC,
  65      2  RESG,RESK,RESKH,RESULT,UFLOW,WG,WGK,XGK
  66       INTEGER J,JTW,JTWM1
  67       EXTERNAL F
  68 C
  69       DIMENSION FV1(15),FV2(15),XGK(16),WGK(16),WG(8)
  70 C
  71 C           THE ABSCISSAE AND WEIGHTS ARE GIVEN FOR THE INTERVAL (-1,1).
  72 C           BECAUSE OF SYMMETRY ONLY THE POSITIVE ABSCISSAE AND THEIR
  73 C           CORRESPONDING WEIGHTS ARE GIVEN.
  74 C
  75 C           XGK    - ABSCISSAE OF THE 31-POINT KRONROD RULE
  76 C                    XGK(2), XGK(4), ...  ABSCISSAE OF THE 15-POINT
  77 C                    GAUSS RULE
  78 C                    XGK(1), XGK(3), ...  ABSCISSAE WHICH ARE OPTIMALLY
  79 C                    ADDED TO THE 15-POINT GAUSS RULE
  80 C
  81 C           WGK    - WEIGHTS OF THE 31-POINT KRONROD RULE
  82 C
  83 C           WG     - WEIGHTS OF THE 15-POINT GAUSS RULE
  84 C
  85 C
  86 C GAUSS QUADRATURE WEIGHTS AND KRONROD QUADRATURE ABSCISSAE AND WEIGHTS
  87 C AS EVALUATED WITH 80 DECIMAL DIGIT ARITHMETIC BY L. W. FULLERTON,
  88 C BELL LABS, NOV. 1981.
  89 C
  90       SAVE WG, XGK, WGK
  91       DATA WG  (  1) / 0.0307532419 9611726835 4628393577 204 D0 /
  92       DATA WG  (  2) / 0.0703660474 8810812470 9267416450 667 D0 /
  93       DATA WG  (  3) / 0.1071592204 6717193501 1869546685 869 D0 /
  94       DATA WG  (  4) / 0.1395706779 2615431444 7804794511 028 D0 /
  95       DATA WG  (  5) / 0.1662692058 1699393355 3200860481 209 D0 /
  96       DATA WG  (  6) / 0.1861610000 1556221102 6800561866 423 D0 /
  97       DATA WG  (  7) / 0.1984314853 2711157645 6118326443 839 D0 /
  98       DATA WG  (  8) / 0.2025782419 2556127288 0620199967 519 D0 /
  99 C
 100       DATA XGK (  1) / 0.9980022986 9339706028 5172840152 271 D0 /
 101       DATA XGK (  2) / 0.9879925180 2048542848 9565718586 613 D0 /
 102       DATA XGK (  3) / 0.9677390756 7913913425 7347978784 337 D0 /
 103       DATA XGK (  4) / 0.9372733924 0070590430 7758947710 209 D0 /
 104       DATA XGK (  5) / 0.8972645323 4408190088 2509656454 496 D0 /
 105       DATA XGK (  6) / 0.8482065834 1042721620 0648320774 217 D0 /
 106       DATA XGK (  7) / 0.7904185014 4246593296 7649294817 947 D0 /
 107       DATA XGK (  8) / 0.7244177313 6017004741 6186054613 938 D0 /
 108       DATA XGK (  9) / 0.6509967412 9741697053 3735895313 275 D0 /
 109       DATA XGK ( 10) / 0.5709721726 0853884753 7226737253 911 D0 /
 110       DATA XGK ( 11) / 0.4850818636 4023968069 3655740232 351 D0 /
 111       DATA XGK ( 12) / 0.3941513470 7756336989 7207370981 045 D0 /
 112       DATA XGK ( 13) / 0.2991800071 5316881216 6780024266 389 D0 /
 113       DATA XGK ( 14) / 0.2011940939 9743452230 0628303394 596 D0 /
 114       DATA XGK ( 15) / 0.1011420669 1871749902 7074231447 392 D0 /
 115       DATA XGK ( 16) / 0.0000000000 0000000000 0000000000 000 D0 /
 116 C
 117       DATA WGK (  1) / 0.0053774798 7292334898 7792051430 128 D0 /
 118       DATA WGK (  2) / 0.0150079473 2931612253 8374763075 807 D0 /
 119       DATA WGK (  3) / 0.0254608473 2671532018 6874001019 653 D0 /
 120       DATA WGK (  4) / 0.0353463607 9137584622 2037948478 360 D0 /
 121       DATA WGK (  5) / 0.0445897513 2476487660 8227299373 280 D0 /
 122       DATA WGK (  6) / 0.0534815246 9092808726 5343147239 430 D0 /
 123       DATA WGK (  7) / 0.0620095678 0067064028 5139230960 803 D0 /
 124       DATA WGK (  8) / 0.0698541213 1872825870 9520077099 147 D0 /
 125       DATA WGK (  9) / 0.0768496807 5772037889 4432777482 659 D0 /
 126       DATA WGK ( 10) / 0.0830805028 2313302103 8289247286 104 D0 /
 127       DATA WGK ( 11) / 0.0885644430 5621177064 7275443693 774 D0 /
 128       DATA WGK ( 12) / 0.0931265981 7082532122 5486872747 346 D0 /
 129       DATA WGK ( 13) / 0.0966427269 8362367850 5179907627 589 D0 /
 130       DATA WGK ( 14) / 0.0991735987 2179195933 2393173484 603 D0 /
 131       DATA WGK ( 15) / 0.1007698455 2387559504 4946662617 570 D0 /
 132       DATA WGK ( 16) / 0.1013300070 1479154901 7374792767 493 D0 /
 133 C
 134 C
 135 C           LIST OF MAJOR VARIABLES
 136 C           -----------------------
 137 C           CENTR  - MID POINT OF THE INTERVAL
 138 C           HLGTH  - HALF-LENGTH OF THE INTERVAL
 139 C           ABSC   - ABSCISSA
 140 C           FVAL*  - FUNCTION VALUE
 141 C           RESG   - RESULT OF THE 15-POINT GAUSS FORMULA
 142 C           RESK   - RESULT OF THE 31-POINT KRONROD FORMULA
 143 C           RESKH  - APPROXIMATION TO THE MEAN VALUE OF F OVER (A,B),
 144 C                    I.E. TO I/(B-A)
 145 C
 146 C           MACHINE DEPENDENT CONSTANTS
 147 C           ---------------------------
 148 C           EPMACH IS THE LARGEST RELATIVE SPACING.
 149 C           UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
 150 C***FIRST EXECUTABLE STATEMENT  DQK31
 151       EPMACH = D1MACH(4)
 152       UFLOW = D1MACH(1)
 153 C
 154       CENTR = 0.5D+00*(A+B)
 155       HLGTH = 0.5D+00*(B-A)
 156       DHLGTH = ABS(HLGTH)
 157 C
 158 C           COMPUTE THE 31-POINT KRONROD APPROXIMATION TO
 159 C           THE INTEGRAL, AND ESTIMATE THE ABSOLUTE ERROR.
 160 C
 161       FC = F(CENTR)
 162       RESG = WG(8)*FC
 163       RESK = WGK(16)*FC
 164       RESABS = ABS(RESK)
 165       DO 10 J=1,7
 166         JTW = J*2
 167         ABSC = HLGTH*XGK(JTW)
 168         FVAL1 = F(CENTR-ABSC)
 169         FVAL2 = F(CENTR+ABSC)
 170         FV1(JTW) = FVAL1
 171         FV2(JTW) = FVAL2
 172         FSUM = FVAL1+FVAL2
 173         RESG = RESG+WG(J)*FSUM
 174         RESK = RESK+WGK(JTW)*FSUM
 175         RESABS = RESABS+WGK(JTW)*(ABS(FVAL1)+ABS(FVAL2))
 176    10 CONTINUE
 177       DO 15 J = 1,8
 178         JTWM1 = J*2-1
 179         ABSC = HLGTH*XGK(JTWM1)
 180         FVAL1 = F(CENTR-ABSC)
 181         FVAL2 = F(CENTR+ABSC)
 182         FV1(JTWM1) = FVAL1
 183         FV2(JTWM1) = FVAL2
 184         FSUM = FVAL1+FVAL2
 185         RESK = RESK+WGK(JTWM1)*FSUM
 186         RESABS = RESABS+WGK(JTWM1)*(ABS(FVAL1)+ABS(FVAL2))
 187    15 CONTINUE
 188       RESKH = RESK*0.5D+00
 189       RESASC = WGK(16)*ABS(FC-RESKH)
 190       DO 20 J=1,15
 191         RESASC = RESASC+WGK(J)*(ABS(FV1(J)-RESKH)+ABS(FV2(J)-RESKH))
 192    20 CONTINUE
 193       RESULT = RESK*HLGTH
 194       RESABS = RESABS*DHLGTH
 195       RESASC = RESASC*DHLGTH
 196       ABSERR = ABS((RESK-RESG)*HLGTH)
 197       IF(RESASC.NE.0.0D+00.AND.ABSERR.NE.0.0D+00)
 198      1  ABSERR = RESASC*MIN(0.1D+01,(0.2D+03*ABSERR/RESASC)**1.5D+00)
 199       IF(RESABS.GT.UFLOW/(0.5D+02*EPMACH)) ABSERR = MAX
 200      1  ((EPMACH*0.5D+02)*RESABS,ABSERR)
 201       RETURN
 202       END