src/numerical/slatec/fortran/dqk41.f

   1 *DECK DQK41
   2       SUBROUTINE DQK41 (F, A, B, RESULT, ABSERR, RESABS, RESASC)
   3 C***BEGIN PROLOGUE  DQK41
   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 (QK41-S, DQK41-D)
  10 C***KEYWORDS  41-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 41-POINT
  40 C                       GAUSS-KRONROD RULE (RESK) obtained by optimal
  41 C                       addition of abscissae to the 20-POINT GAUSS
  42 C                       RULE (RESG).
  43 C
  44 C              ABSERR - Double precision
  45 C                       Estimate of the modulus of the absolute error,
  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  DQK41
  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(20),FV2(20),XGK(21),WGK(21),WG(10)
  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 41-POINT GAUSS-KRONROD RULE
  77 C                    XGK(2), XGK(4), ...  ABSCISSAE OF THE 20-POINT
  78 C                    GAUSS RULE
  79 C                    XGK(1), XGK(3), ...  ABSCISSAE WHICH ARE OPTIMALLY
  80 C                    ADDED TO THE 20-POINT GAUSS RULE
  81 C
  82 C           WGK    - WEIGHTS OF THE 41-POINT GAUSS-KRONROD RULE
  83 C
  84 C           WG     - WEIGHTS OF THE 20-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.0176140071 3915211831 1861962351 853 D0 /
  93       DATA WG  (  2) / 0.0406014298 0038694133 1039952274 932 D0 /
  94       DATA WG  (  3) / 0.0626720483 3410906356 9506535187 042 D0 /
  95       DATA WG  (  4) / 0.0832767415 7670474872 4758143222 046 D0 /
  96       DATA WG  (  5) / 0.1019301198 1724043503 6750135480 350 D0 /
  97       DATA WG  (  6) / 0.1181945319 6151841731 2377377711 382 D0 /
  98       DATA WG  (  7) / 0.1316886384 4917662689 8494499748 163 D0 /
  99       DATA WG  (  8) / 0.1420961093 1838205132 9298325067 165 D0 /
 100       DATA WG  (  9) / 0.1491729864 7260374678 7828737001 969 D0 /
 101       DATA WG  ( 10) / 0.1527533871 3072585069 8084331955 098 D0 /
 102 C
 103       DATA XGK (  1) / 0.9988590315 8827766383 8315576545 863 D0 /
 104       DATA XGK (  2) / 0.9931285991 8509492478 6122388471 320 D0 /
 105       DATA XGK (  3) / 0.9815078774 5025025919 3342994720 217 D0 /
 106       DATA XGK (  4) / 0.9639719272 7791379126 7666131197 277 D0 /
 107       DATA XGK (  5) / 0.9408226338 3175475351 9982722212 443 D0 /
 108       DATA XGK (  6) / 0.9122344282 5132590586 7752441203 298 D0 /
 109       DATA XGK (  7) / 0.8782768112 5228197607 7442995113 078 D0 /
 110       DATA XGK (  8) / 0.8391169718 2221882339 4529061701 521 D0 /
 111       DATA XGK (  9) / 0.7950414288 3755119835 0638833272 788 D0 /
 112       DATA XGK ( 10) / 0.7463319064 6015079261 4305070355 642 D0 /
 113       DATA XGK ( 11) / 0.6932376563 3475138480 5490711845 932 D0 /
 114       DATA XGK ( 12) / 0.6360536807 2651502545 2836696226 286 D0 /
 115       DATA XGK ( 13) / 0.5751404468 1971031534 2946036586 425 D0 /
 116       DATA XGK ( 14) / 0.5108670019 5082709800 4364050955 251 D0 /
 117       DATA XGK ( 15) / 0.4435931752 3872510319 9992213492 640 D0 /
 118       DATA XGK ( 16) / 0.3737060887 1541956067 2548177024 927 D0 /
 119       DATA XGK ( 17) / 0.3016278681 1491300432 0555356858 592 D0 /
 120       DATA XGK ( 18) / 0.2277858511 4164507808 0496195368 575 D0 /
 121       DATA XGK ( 19) / 0.1526054652 4092267550 5220241022 678 D0 /
 122       DATA XGK ( 20) / 0.0765265211 3349733375 4640409398 838 D0 /
 123       DATA XGK ( 21) / 0.0000000000 0000000000 0000000000 000 D0 /
 124 C
 125       DATA WGK (  1) / 0.0030735837 1852053150 1218293246 031 D0 /
 126       DATA WGK (  2) / 0.0086002698 5564294219 8661787950 102 D0 /
 127       DATA WGK (  3) / 0.0146261692 5697125298 3787960308 868 D0 /
 128       DATA WGK (  4) / 0.0203883734 6126652359 8010231432 755 D0 /
 129       DATA WGK (  5) / 0.0258821336 0495115883 4505067096 153 D0 /
 130       DATA WGK (  6) / 0.0312873067 7703279895 8543119323 801 D0 /
 131       DATA WGK (  7) / 0.0366001697 5820079803 0557240707 211 D0 /
 132       DATA WGK (  8) / 0.0416688733 2797368626 3788305936 895 D0 /
 133       DATA WGK (  9) / 0.0464348218 6749767472 0231880926 108 D0 /
 134       DATA WGK ( 10) / 0.0509445739 2372869193 2707670050 345 D0 /
 135       DATA WGK ( 11) / 0.0551951053 4828599474 4832372419 777 D0 /
 136       DATA WGK ( 12) / 0.0591114008 8063957237 4967220648 594 D0 /
 137       DATA WGK ( 13) / 0.0626532375 5478116802 5870122174 255 D0 /
 138       DATA WGK ( 14) / 0.0658345971 3361842211 1563556969 398 D0 /
 139       DATA WGK ( 15) / 0.0686486729 2852161934 5623411885 368 D0 /
 140       DATA WGK ( 16) / 0.0710544235 5344406830 5790361723 210 D0 /
 141       DATA WGK ( 17) / 0.0730306903 3278666749 5189417658 913 D0 /
 142       DATA WGK ( 18) / 0.0745828754 0049918898 6581418362 488 D0 /
 143       DATA WGK ( 19) / 0.0757044976 8455667465 9542775376 617 D0 /
 144       DATA WGK ( 20) / 0.0763778676 7208073670 5502835038 061 D0 /
 145       DATA WGK ( 21) / 0.0766007119 1799965644 5049901530 102 D0 /
 146 C
 147 C
 148 C           LIST OF MAJOR VARIABLES
 149 C           -----------------------
 150 C
 151 C           CENTR  - MID POINT OF THE INTERVAL
 152 C           HLGTH  - HALF-LENGTH OF THE INTERVAL
 153 C           ABSC   - ABSCISSA
 154 C           FVAL*  - FUNCTION VALUE
 155 C           RESG   - RESULT OF THE 20-POINT GAUSS FORMULA
 156 C           RESK   - RESULT OF THE 41-POINT KRONROD FORMULA
 157 C           RESKH  - APPROXIMATION TO MEAN VALUE OF F OVER (A,B), I.E.
 158 C                    TO I/(B-A)
 159 C
 160 C           MACHINE DEPENDENT CONSTANTS
 161 C           ---------------------------
 162 C
 163 C           EPMACH IS THE LARGEST RELATIVE SPACING.
 164 C           UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
 165 C
 166 C***FIRST EXECUTABLE STATEMENT  DQK41
 167       EPMACH = D1MACH(4)
 168       UFLOW = D1MACH(1)
 169 C
 170       CENTR = 0.5D+00*(A+B)
 171       HLGTH = 0.5D+00*(B-A)
 172       DHLGTH = ABS(HLGTH)
 173 C
 174 C           COMPUTE THE 41-POINT GAUSS-KRONROD APPROXIMATION TO
 175 C           THE INTEGRAL, AND ESTIMATE THE ABSOLUTE ERROR.
 176 C
 177       RESG = 0.0D+00
 178       FC = F(CENTR)
 179       RESK = WGK(21)*FC
 180       RESABS = ABS(RESK)
 181       DO 10 J=1,10
 182         JTW = J*2
 183         ABSC = HLGTH*XGK(JTW)
 184         FVAL1 = F(CENTR-ABSC)
 185         FVAL2 = F(CENTR+ABSC)
 186         FV1(JTW) = FVAL1
 187         FV2(JTW) = FVAL2
 188         FSUM = FVAL1+FVAL2
 189         RESG = RESG+WG(J)*FSUM
 190         RESK = RESK+WGK(JTW)*FSUM
 191         RESABS = RESABS+WGK(JTW)*(ABS(FVAL1)+ABS(FVAL2))
 192    10 CONTINUE
 193       DO 15 J = 1,10
 194         JTWM1 = J*2-1
 195         ABSC = HLGTH*XGK(JTWM1)
 196         FVAL1 = F(CENTR-ABSC)
 197         FVAL2 = F(CENTR+ABSC)
 198         FV1(JTWM1) = FVAL1
 199         FV2(JTWM1) = FVAL2
 200         FSUM = FVAL1+FVAL2
 201         RESK = RESK+WGK(JTWM1)*FSUM
 202         RESABS = RESABS+WGK(JTWM1)*(ABS(FVAL1)+ABS(FVAL2))
 203    15 CONTINUE
 204       RESKH = RESK*0.5D+00
 205       RESASC = WGK(21)*ABS(FC-RESKH)
 206       DO 20 J=1,20
 207         RESASC = RESASC+WGK(J)*(ABS(FV1(J)-RESKH)+ABS(FV2(J)-RESKH))
 208    20 CONTINUE
 209       RESULT = RESK*HLGTH
 210       RESABS = RESABS*DHLGTH
 211       RESASC = RESASC*DHLGTH
 212       ABSERR = ABS((RESK-RESG)*HLGTH)
 213       IF(RESASC.NE.0.0D+00.AND.ABSERR.NE.0.D+00)
 214      1  ABSERR = RESASC*MIN(0.1D+01,(0.2D+03*ABSERR/RESASC)**1.5D+00)
 215       IF(RESABS.GT.UFLOW/(0.5D+02*EPMACH)) ABSERR = MAX
 216      1  ((EPMACH*0.5D+02)*RESABS,ABSERR)
 217       RETURN
 218       END