Home Explore Help
Sign In
DragonOS-Community
/
relibc_openlibm
mirror of https://gitlab.redox-os.org/redox-os/openlibm.git
2
0
Fork 0
Files Issues 0 Wiki
Tree: 50ce8eff11
Branches Tags
relibc_openlibm
/
slatec
/
dqk41.f

dqk41.f 9.0 KB
History Raw

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218 
  1. *DECK DQK41
    2. 

  2.       SUBROUTINE DQK41 (F, A, B, RESULT, ABSERR, RESABS, RESASC)
    3. 

  3. C***BEGIN PROLOGUE  DQK41
    4. 

  4. C***PURPOSE  To compute I = Integral of F over (A,B), with error
    5. 

  5. C                           estimate
    6. 

  6. C                       J = Integral of ABS(F) over (A,B)
    7. 

  7. C***LIBRARY   SLATEC (QUADPACK)
    8. 

  8. C***CATEGORY  H2A1A2
    9. 

  9. C***TYPE      DOUBLE PRECISION (QK41-S, DQK41-D)
    10. 

  10. C***KEYWORDS  41-POINT GAUSS-KRONROD RULES, QUADPACK, QUADRATURE
    11. 

  11. C***AUTHOR  Piessens, Robert
    12. 

  12. C             Applied Mathematics and Programming Division
    13. 

  13. C             K. U. Leuven
    14. 

  14. C           de Doncker, Elise
    15. 

  15. C             Applied Mathematics and Programming Division
    16. 

  16. C             K. U. Leuven
    17. 

  17. C***DESCRIPTION
    18. 

  18. C
    19. 

  19. C           Integration rules
    20. 

  20. C           Standard fortran subroutine
    21. 

  21. C           Double precision version
    22. 

  22. C
    23. 

  23. C           PARAMETERS
    24. 

  24. C            ON ENTRY
    25. 

  25. C              F      - Double precision
    26. 

  26. C                       Function subprogram defining the integrand
    27. 

  27. C                       FUNCTION F(X). The actual name for F needs to be
    28. 

  28. C                       declared E X T E R N A L in the calling program.
    29. 

  29. C
    30. 

  30. C              A      - Double precision
    31. 

  31. C                       Lower limit of integration
    32. 

  32. C
    33. 

  33. C              B      - Double precision
    34. 

  34. C                       Upper limit of integration
    35. 

  35. C
    36. 

  36. C            ON RETURN
    37. 

  37. C              RESULT - Double precision
    38. 

  38. C                       Approximation to the integral I
    39. 

  39. C                       RESULT is computed by applying the 41-POINT
    40. 

  40. C                       GAUSS-KRONROD RULE (RESK) obtained by optimal
    41. 

  41. C                       addition of abscissae to the 20-POINT GAUSS
    42. 

  42. C                       RULE (RESG).
    43. 

  43. C
    44. 

  44. C              ABSERR - Double precision
    45. 

  45. C                       Estimate of the modulus of the absolute error,
    46. 

  46. C                       which should not exceed ABS(I-RESULT)
    47. 

  47. C
    48. 

  48. C              RESABS - Double precision
    49. 

  49. C                       Approximation to the integral J
    50. 

  50. C
    51. 

  51. C              RESASC - Double precision
    52. 

  52. C                       Approximation to the integral of ABS(F-I/(B-A))
    53. 

  53. C                       over (A,B)
    54. 

  54. C
    55. 

  55. C***REFERENCES  (NONE)
    56. 

  56. C***ROUTINES CALLED  D1MACH
    57. 

  57. C***REVISION HISTORY  (YYMMDD)
    58. 

  58. C   800101  DATE WRITTEN
    59. 

  59. C   890531  Changed all specific intrinsics to generic.  (WRB)
    60. 

  60. C   890531  REVISION DATE from Version 3.2
    61. 

  61. C   891214  Prologue converted to Version 4.0 format.  (BAB)
    62. 

  62. C***END PROLOGUE  DQK41
    63. 

  63. C
    64. 

  64.       DOUBLE PRECISION A,ABSC,ABSERR,B,CENTR,DHLGTH,
    65. 

  65.      1  D1MACH,EPMACH,F,FC,FSUM,FVAL1,FVAL2,FV1,FV2,HLGTH,RESABS,RESASC,
    66. 

  66.      2  RESG,RESK,RESKH,RESULT,UFLOW,WG,WGK,XGK
    67. 

  67.       INTEGER J,JTW,JTWM1
    68. 

  68.       EXTERNAL F
    69. 

  69. C
    70. 

  70.       DIMENSION FV1(20),FV2(20),XGK(21),WGK(21),WG(10)
    71. 

  71. C
    72. 

  72. C           THE ABSCISSAE AND WEIGHTS ARE GIVEN FOR THE INTERVAL (-1,1).
    73. 

  73. C           BECAUSE OF SYMMETRY ONLY THE POSITIVE ABSCISSAE AND THEIR
    74. 

  74. C           CORRESPONDING WEIGHTS ARE GIVEN.
    75. 

  75. C
    76. 

  76. C           XGK    - ABSCISSAE OF THE 41-POINT GAUSS-KRONROD RULE
    77. 

  77. C                    XGK(2), XGK(4), ...  ABSCISSAE OF THE 20-POINT
    78. 

  78. C                    GAUSS RULE
    79. 

  79. C                    XGK(1), XGK(3), ...  ABSCISSAE WHICH ARE OPTIMALLY
    80. 

  80. C                    ADDED TO THE 20-POINT GAUSS RULE
    81. 

  81. C
    82. 

  82. C           WGK    - WEIGHTS OF THE 41-POINT GAUSS-KRONROD RULE
    83. 

  83. C
    84. 

  84. C           WG     - WEIGHTS OF THE 20-POINT GAUSS RULE
    85. 

  85. C
    86. 

  86. C
    87. 

  87. C GAUSS QUADRATURE WEIGHTS AND KRONROD QUADRATURE ABSCISSAE AND WEIGHTS
    88. 

  88. C AS EVALUATED WITH 80 DECIMAL DIGIT ARITHMETIC BY L. W. FULLERTON,
    89. 

  89. C BELL LABS, NOV. 1981.
    90. 

  90. C
    91. 

  91.       SAVE WG, XGK, WGK
    92. 

  92.       DATA WG  (  1) / 0.0176140071 3915211831 1861962351 853 D0 /
    93. 

  93.       DATA WG  (  2) / 0.0406014298 0038694133 1039952274 932 D0 /
    94. 

  94.       DATA WG  (  3) / 0.0626720483 3410906356 9506535187 042 D0 /
    95. 

  95.       DATA WG  (  4) / 0.0832767415 7670474872 4758143222 046 D0 /
    96. 

  96.       DATA WG  (  5) / 0.1019301198 1724043503 6750135480 350 D0 /
    97. 

  97.       DATA WG  (  6) / 0.1181945319 6151841731 2377377711 382 D0 /
    98. 

  98.       DATA WG  (  7) / 0.1316886384 4917662689 8494499748 163 D0 /
    99. 

  99.       DATA WG  (  8) / 0.1420961093 1838205132 9298325067 165 D0 /
    100. 

  100.       DATA WG  (  9) / 0.1491729864 7260374678 7828737001 969 D0 /
    101. 

  101.       DATA WG  ( 10) / 0.1527533871 3072585069 8084331955 098 D0 /
    102. 

  102. C
    103. 

  103.       DATA XGK (  1) / 0.9988590315 8827766383 8315576545 863 D0 /
    104. 

  104.       DATA XGK (  2) / 0.9931285991 8509492478 6122388471 320 D0 /
    105. 

  105.       DATA XGK (  3) / 0.9815078774 5025025919 3342994720 217 D0 /
    106. 

  106.       DATA XGK (  4) / 0.9639719272 7791379126 7666131197 277 D0 /
    107. 

  107.       DATA XGK (  5) / 0.9408226338 3175475351 9982722212 443 D0 /
    108. 

  108.       DATA XGK (  6) / 0.9122344282 5132590586 7752441203 298 D0 /
    109. 

  109.       DATA XGK (  7) / 0.8782768112 5228197607 7442995113 078 D0 /
    110. 

  110.       DATA XGK (  8) / 0.8391169718 2221882339 4529061701 521 D0 /
    111. 

  111.       DATA XGK (  9) / 0.7950414288 3755119835 0638833272 788 D0 /
    112. 

  112.       DATA XGK ( 10) / 0.7463319064 6015079261 4305070355 642 D0 /
    113. 

  113.       DATA XGK ( 11) / 0.6932376563 3475138480 5490711845 932 D0 /
    114. 

  114.       DATA XGK ( 12) / 0.6360536807 2651502545 2836696226 286 D0 /
    115. 

  115.       DATA XGK ( 13) / 0.5751404468 1971031534 2946036586 425 D0 /
    116. 

  116.       DATA XGK ( 14) / 0.5108670019 5082709800 4364050955 251 D0 /
    117. 

  117.       DATA XGK ( 15) / 0.4435931752 3872510319 9992213492 640 D0 /
    118. 

  118.       DATA XGK ( 16) / 0.3737060887 1541956067 2548177024 927 D0 /
    119. 

  119.       DATA XGK ( 17) / 0.3016278681 1491300432 0555356858 592 D0 /
    120. 

  120.       DATA XGK ( 18) / 0.2277858511 4164507808 0496195368 575 D0 /
    121. 

  121.       DATA XGK ( 19) / 0.1526054652 4092267550 5220241022 678 D0 /
    122. 

  122.       DATA XGK ( 20) / 0.0765265211 3349733375 4640409398 838 D0 /
    123. 

  123.       DATA XGK ( 21) / 0.0000000000 0000000000 0000000000 000 D0 /
    124. 

  124. C
    125. 

  125.       DATA WGK (  1) / 0.0030735837 1852053150 1218293246 031 D0 /
    126. 

  126.       DATA WGK (  2) / 0.0086002698 5564294219 8661787950 102 D0 /
    127. 

  127.       DATA WGK (  3) / 0.0146261692 5697125298 3787960308 868 D0 /
    128. 

  128.       DATA WGK (  4) / 0.0203883734 6126652359 8010231432 755 D0 /
    129. 

  129.       DATA WGK (  5) / 0.0258821336 0495115883 4505067096 153 D0 /
    130. 

  130.       DATA WGK (  6) / 0.0312873067 7703279895 8543119323 801 D0 /
    131. 

  131.       DATA WGK (  7) / 0.0366001697 5820079803 0557240707 211 D0 /
    132. 

  132.       DATA WGK (  8) / 0.0416688733 2797368626 3788305936 895 D0 /
    133. 

  133.       DATA WGK (  9) / 0.0464348218 6749767472 0231880926 108 D0 /
    134. 

  134.       DATA WGK ( 10) / 0.0509445739 2372869193 2707670050 345 D0 /
    135. 

  135.       DATA WGK ( 11) / 0.0551951053 4828599474 4832372419 777 D0 /
    136. 

  136.       DATA WGK ( 12) / 0.0591114008 8063957237 4967220648 594 D0 /
    137. 

  137.       DATA WGK ( 13) / 0.0626532375 5478116802 5870122174 255 D0 /
    138. 

  138.       DATA WGK ( 14) / 0.0658345971 3361842211 1563556969 398 D0 /
    139. 

  139.       DATA WGK ( 15) / 0.0686486729 2852161934 5623411885 368 D0 /
    140. 

  140.       DATA WGK ( 16) / 0.0710544235 5344406830 5790361723 210 D0 /
    141. 

  141.       DATA WGK ( 17) / 0.0730306903 3278666749 5189417658 913 D0 /
    142. 

  142.       DATA WGK ( 18) / 0.0745828754 0049918898 6581418362 488 D0 /
    143. 

  143.       DATA WGK ( 19) / 0.0757044976 8455667465 9542775376 617 D0 /
    144. 

  144.       DATA WGK ( 20) / 0.0763778676 7208073670 5502835038 061 D0 /
    145. 

  145.       DATA WGK ( 21) / 0.0766007119 1799965644 5049901530 102 D0 /
    146. 

  146. C
    147. 

  147. C
    148. 

  148. C           LIST OF MAJOR VARIABLES
    149. 

  149. C           -----------------------
    150. 

  150. C
    151. 

  151. C           CENTR  - MID POINT OF THE INTERVAL
    152. 

  152. C           HLGTH  - HALF-LENGTH OF THE INTERVAL
    153. 

  153. C           ABSC   - ABSCISSA
    154. 

  154. C           FVAL*  - FUNCTION VALUE
    155. 

  155. C           RESG   - RESULT OF THE 20-POINT GAUSS FORMULA
    156. 

  156. C           RESK   - RESULT OF THE 41-POINT KRONROD FORMULA
    157. 

  157. C           RESKH  - APPROXIMATION TO MEAN VALUE OF F OVER (A,B), I.E.
    158. 

  158. C                    TO I/(B-A)
    159. 

  159. C
    160. 

  160. C           MACHINE DEPENDENT CONSTANTS
    161. 

  161. C           ---------------------------
    162. 

  162. C
    163. 

  163. C           EPMACH IS THE LARGEST RELATIVE SPACING.
    164. 

  164. C           UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
    165. 

  165. C
    166. 

  166. C***FIRST EXECUTABLE STATEMENT  DQK41
    167. 

  167.       EPMACH = D1MACH(4)
    168. 

  168.       UFLOW = D1MACH(1)
    169. 

  169. C
    170. 

  170.       CENTR = 0.5D+00*(A+B)
    171. 

  171.       HLGTH = 0.5D+00*(B-A)
    172. 

  172.       DHLGTH = ABS(HLGTH)
    173. 

  173. C
    174. 

  174. C           COMPUTE THE 41-POINT GAUSS-KRONROD APPROXIMATION TO
    175. 

  175. C           THE INTEGRAL, AND ESTIMATE THE ABSOLUTE ERROR.
    176. 

  176. C
    177. 

  177.       RESG = 0.0D+00
    178. 

  178.       FC = F(CENTR)
    179. 

  179.       RESK = WGK(21)*FC
    180. 

  180.       RESABS = ABS(RESK)
    181. 

  181.       DO 10 J=1,10
    182. 

  182.         JTW = J*2
    183. 

  183.         ABSC = HLGTH*XGK(JTW)
    184. 

  184.         FVAL1 = F(CENTR-ABSC)
    185. 

  185.         FVAL2 = F(CENTR+ABSC)
    186. 

  186.         FV1(JTW) = FVAL1
    187. 

  187.         FV2(JTW) = FVAL2
    188. 

  188.         FSUM = FVAL1+FVAL2
    189. 

  189.         RESG = RESG+WG(J)*FSUM
    190. 

  190.         RESK = RESK+WGK(JTW)*FSUM
    191. 

  191.         RESABS = RESABS+WGK(JTW)*(ABS(FVAL1)+ABS(FVAL2))
    192. 

  192.    10 CONTINUE
    193. 

  193.       DO 15 J = 1,10
    194. 

  194.         JTWM1 = J*2-1
    195. 

  195.         ABSC = HLGTH*XGK(JTWM1)
    196. 

  196.         FVAL1 = F(CENTR-ABSC)
    197. 

  197.         FVAL2 = F(CENTR+ABSC)
    198. 

  198.         FV1(JTWM1) = FVAL1
    199. 

  199.         FV2(JTWM1) = FVAL2
    200. 

  200.         FSUM = FVAL1+FVAL2
    201. 

  201.         RESK = RESK+WGK(JTWM1)*FSUM
    202. 

  202.         RESABS = RESABS+WGK(JTWM1)*(ABS(FVAL1)+ABS(FVAL2))
    203. 

  203.    15 CONTINUE
    204. 

  204.       RESKH = RESK*0.5D+00
    205. 

  205.       RESASC = WGK(21)*ABS(FC-RESKH)
    206. 

  206.       DO 20 J=1,20
    207. 

  207.         RESASC = RESASC+WGK(J)*(ABS(FV1(J)-RESKH)+ABS(FV2(J)-RESKH))
    208. 

  208.    20 CONTINUE
    209. 

  209.       RESULT = RESK*HLGTH
    210. 

  210.       RESABS = RESABS*DHLGTH
    211. 

  211.       RESASC = RESASC*DHLGTH
    212. 

  212.       ABSERR = ABS((RESK-RESG)*HLGTH)
    213. 

  213.       IF(RESASC.NE.0.0D+00.AND.ABSERR.NE.0.D+00)
    214. 

  214.      1  ABSERR = RESASC*MIN(0.1D+01,(0.2D+03*ABSERR/RESASC)**1.5D+00)
    215. 

  215.       IF(RESABS.GT.UFLOW/(0.5D+02*EPMACH)) ABSERR = MAX
    216. 

  216.      1  ((EPMACH*0.5D+02)*RESABS,ABSERR)
    217. 

  217.       RETURN
    218. 

  218.       END
    219.