home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-src.tgz / tar.out / fsf / octave / libcruft / quadpack / dqk15i.f < prev    next >
Text File  |  1996-09-28  |  8KB  |  213 lines
  1.       SUBROUTINE DQK15I(F,BOUN,INF,A,B,RESULT,ABSERR,RESABS,RESASC,
  2.      1   IERR)
  3. C***BEGIN PROLOGUE  DQK15I
  4. C***DATE WRITTEN   800101   (YYMMDD)
  5. C***REVISION DATE  830518   (YYMMDD)
  6. C***CATEGORY NO.  H2A3A2,H2A4A2
  7. C***KEYWORDS  15-POINT TRANSFORMED GAUSS-KRONROD RULES
  8. C***AUTHOR  PIESSENS,ROBERT,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN
  9. C           DE DONCKER,ELISE,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN
  10. C***PURPOSE  THE ORIGINAL (INFINITE INTEGRATION RANGE IS MAPPED
  11. C            ONTO THE INTERVAL (0,1) AND (A,B) IS A PART OF (0,1).
  12. C            IT IS THE PURPOSE TO COMPUTE
  13. C            I = INTEGRAL OF TRANSFORMED INTEGRAND OVER (A,B),
  14. C            J = INTEGRAL OF ABS(TRANSFORMED INTEGRAND) OVER (A,B).
  15. C***DESCRIPTION
  16. C
  17. C           INTEGRATION RULE
  18. C           STANDARD FORTRAN SUBROUTINE
  19. C           DOUBLE PRECISION VERSION
  20. C
  21. C           PARAMETERS
  22. C            ON ENTRY
  23. C              F      - DOUBLE PRECISION
  24. C                       FUCTION SUBPROGRAM DEFINING THE INTEGRAND
  25. C                       FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
  26. C                       DECLARED E X T E R N A L IN THE CALLING PROGRAM.
  27. C
  28. C              BOUN   - DOUBLE PRECISION
  29. C                       FINITE BOUND OF ORIGINAL INTEGRATION
  30. C                       RANGE (SET TO ZERO IF INF = +2)
  31. C
  32. C              INF    - INTEGER
  33. C                       IF INF = -1, THE ORIGINAL INTERVAL IS
  34. C                                   (-INFINITY,BOUND),
  35. C                       IF INF = +1, THE ORIGINAL INTERVAL IS
  36. C                                   (BOUND,+INFINITY),
  37. C                       IF INF = +2, THE ORIGINAL INTERVAL IS
  38. C                                   (-INFINITY,+INFINITY) AND
  39. C                       THE INTEGRAL IS COMPUTED AS THE SUM OF TWO
  40. C                       INTEGRALS, ONE OVER (-INFINITY,0) AND ONE OVER
  41. C                       (0,+INFINITY).
  42. C
  43. C              A      - DOUBLE PRECISION
  44. C                       LOWER LIMIT FOR INTEGRATION OVER SUBRANGE
  45. C                       OF (0,1)
  46. C
  47. C              B      - DOUBLE PRECISION
  48. C                       UPPER LIMIT FOR INTEGRATION OVER SUBRANGE
  49. C                       OF (0,1)
  50. C
  51. C            ON RETURN
  52. C              RESULT - DOUBLE PRECISION
  53. C                       APPROXIMATION TO THE INTEGRAL I
  54. C                       RESULT IS COMPUTED BY APPLYING THE 15-POINT
  55. C                       KRONROD RULE(RESK) OBTAINED BY OPTIMAL ADDITION
  56. C                       OF ABSCISSAE TO THE 7-POINT GAUSS RULE(RESG).
  57. C
  58. C              ABSERR - DOUBLE PRECISION
  59. C                       ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR,
  60. C                       WHICH SHOULD EQUAL OR EXCEED ABS(I-RESULT)
  61. C
  62. C              RESABS - DOUBLE PRECISION
  63. C                       APPROXIMATION TO THE INTEGRAL J
  64. C
  65. C              RESASC - DOUBLE PRECISION
  66. C                       APPROXIMATION TO THE INTEGRAL OF
  67. C                       ABS((TRANSFORMED INTEGRAND)-I/(B-A)) OVER (A,B)
  68. C
  69. C***REFERENCES  (NONE)
  70. C***ROUTINES CALLED  D1MACH
  71. C***END PROLOGUE  DQK15I
  72. C
  73.       DOUBLE PRECISION A,ABSC,ABSC1,ABSC2,ABSERR,B,BOUN,CENTR,DABS,DINF,
  74.      *  DMAX1,DMIN1,D1MACH,EPMACH,F,FC,FSUM,FVAL1,FVAL2,FV1,FV2,HLGTH,
  75.      *  RESABS,RESASC,RESG,RESK,RESKH,RESULT,TABSC1,TABSC2,UFLOW,WG,WGK,
  76.      *  XGK,FVALT
  77.       INTEGER INF,J
  78.       EXTERNAL F
  79. C
  80.       DIMENSION FV1(7),FV2(7),XGK(8),WGK(8),WG(8)
  81. C
  82. C           THE ABSCISSAE AND WEIGHTS ARE SUPPLIED FOR THE INTERVAL
  83. C           (-1,1).  BECAUSE OF SYMMETRY ONLY THE POSITIVE ABSCISSAE AND
  84. C           THEIR CORRESPONDING WEIGHTS ARE GIVEN.
  85. C
  86. C           XGK    - ABSCISSAE OF THE 15-POINT KRONROD RULE
  87. C                    XGK(2), XGK(4), ... ABSCISSAE OF THE 7-POINT
  88. C                    GAUSS RULE
  89. C                    XGK(1), XGK(3), ...  ABSCISSAE WHICH ARE OPTIMALLY
  90. C                    ADDED TO THE 7-POINT GAUSS RULE
  91. C
  92. C           WGK    - WEIGHTS OF THE 15-POINT KRONROD RULE
  93. C
  94. C           WG     - WEIGHTS OF THE 7-POINT GAUSS RULE, CORRESPONDING
  95. C                    TO THE ABSCISSAE XGK(2), XGK(4), ...
  96. C                    WG(1), WG(3), ... ARE SET TO ZERO.
  97. C
  98.       DATA WG(1) / 0.0D0 /
  99.       DATA WG(2) / 0.1294849661 6886969327 0611432679 082D0 /
  100.       DATA WG(3) / 0.0D0 /
  101.       DATA WG(4) / 0.2797053914 8927666790 1467771423 780D0 /
  102.       DATA WG(5) / 0.0D0 /
  103.       DATA WG(6) / 0.3818300505 0511894495 0369775488 975D0 /
  104.       DATA WG(7) / 0.0D0 /
  105.       DATA WG(8) / 0.4179591836 7346938775 5102040816 327D0 /
  106. C
  107.       DATA XGK(1) / 0.9914553711 2081263920 6854697526 329D0 /
  108.       DATA XGK(2) / 0.9491079123 4275852452 6189684047 851D0 /
  109.       DATA XGK(3) / 0.8648644233 5976907278 9712788640 926D0 /
  110.       DATA XGK(4) / 0.7415311855 9939443986 3864773280 788D0 /
  111.       DATA XGK(5) / 0.5860872354 6769113029 4144838258 730D0 /
  112.       DATA XGK(6) / 0.4058451513 7739716690 6606412076 961D0 /
  113.       DATA XGK(7) / 0.2077849550 0789846760 0689403773 245D0 /
  114.       DATA XGK(8) / 0.0000000000 0000000000 0000000000 000D0 /
  115. C
  116.       DATA WGK(1) / 0.0229353220 1052922496 3732008058 970D0 /
  117.       DATA WGK(2) / 0.0630920926 2997855329 0700663189 204D0 /
  118.       DATA WGK(3) / 0.1047900103 2225018383 9876322541 518D0 /
  119.       DATA WGK(4) / 0.1406532597 1552591874 5189590510 238D0 /
  120.       DATA WGK(5) / 0.1690047266 3926790282 6583426598 550D0 /
  121.       DATA WGK(6) / 0.1903505780 6478540991 3256402421 014D0 /
  122.       DATA WGK(7) / 0.2044329400 7529889241 4161999234 649D0 /
  123.       DATA WGK(8) / 0.2094821410 8472782801 2999174891 714D0 /
  124. C
  125. C
  126. C           LIST OF MAJOR VARIABLES
  127. C           -----------------------
  128. C
  129. C           CENTR  - MID POINT OF THE INTERVAL
  130. C           HLGTH  - HALF-LENGTH OF THE INTERVAL
  131. C           ABSC*  - ABSCISSA
  132. C           TABSC* - TRANSFORMED ABSCISSA
  133. C           FVAL*  - FUNCTION VALUE
  134. C           RESG   - RESULT OF THE 7-POINT GAUSS FORMULA
  135. C           RESK   - RESULT OF THE 15-POINT KRONROD FORMULA
  136. C           RESKH  - APPROXIMATION TO THE MEAN VALUE OF THE TRANSFORMED
  137. C                    INTEGRAND OVER (A,B), I.E. TO I/(B-A)
  138. C
  139. C           MACHINE DEPENDENT CONSTANTS
  140. C           ---------------------------
  141. C
  142. C           EPMACH IS THE LARGEST RELATIVE SPACING.
  143. C           UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
  144. C
  145. C***FIRST EXECUTABLE STATEMENT  DQK15I
  146.       EPMACH = D1MACH(4)
  147.       UFLOW = D1MACH(1)
  148.       DINF = MIN0(1,INF)
  149. C
  150.       CENTR = 0.5D+00*(A+B)
  151.       HLGTH = 0.5D+00*(B-A)
  152.       TABSC1 = BOUN+DINF*(0.1D+01-CENTR)/CENTR
  153.       IERR = 0
  154.       FVAL1 = F(TABSC1,IERR)
  155.       IF (IERR .LT. 0) RETURN
  156.       IF(INF.EQ.2) THEN
  157.         FVALT = F(-TABSC1,IERR)
  158.         IF (IERR .LT. 0) RETURN
  159.         FVAL1 = FVAL1+FVALT
  160.       ENDIF
  161.       FC = (FVAL1/CENTR)/CENTR
  162. C
  163. C           COMPUTE THE 15-POINT KRONROD APPROXIMATION TO
  164. C           THE INTEGRAL, AND ESTIMATE THE ERROR.
  165. C
  166.       RESG = WG(8)*FC
  167.       RESK = WGK(8)*FC
  168.       RESABS = DABS(RESK)
  169.       DO 10 J=1,7
  170.         ABSC = HLGTH*XGK(J)
  171.         ABSC1 = CENTR-ABSC
  172.         ABSC2 = CENTR+ABSC
  173.         TABSC1 = BOUN+DINF*(0.1D+01-ABSC1)/ABSC1
  174.         TABSC2 = BOUN+DINF*(0.1D+01-ABSC2)/ABSC2
  175.         FVAL1 = F(TABSC1,IERR)
  176.         IF (IERR .LT. 0) RETURN
  177.         FVAL2 = F(TABSC2,IERR)
  178.         IF (IERR .LT. 0) RETURN
  179.         IF(INF.EQ.2) THEN
  180.           FVALT = F(-TABSC1,IERR)
  181.           IF (IERR .LT. 0) RETURN
  182.           FVAL1 = FVAL1+FVALT
  183.         ENDIF
  184.         IF(INF.EQ.2) THEN
  185.           FVALT = F(-TABSC2,IERR)
  186.           IF (IERR .LT. 0) RETURN
  187.           FVAL2 = FVAL2+FVALT
  188.         ENDIF
  189.         FVAL1 = (FVAL1/ABSC1)/ABSC1
  190.         FVAL2 = (FVAL2/ABSC2)/ABSC2
  191.         FV1(J) = FVAL1
  192.         FV2(J) = FVAL2
  193.         FSUM = FVAL1+FVAL2
  194.         RESG = RESG+WG(J)*FSUM
  195.         RESK = RESK+WGK(J)*FSUM
  196.         RESABS = RESABS+WGK(J)*(DABS(FVAL1)+DABS(FVAL2))
  197.    10 CONTINUE
  198.       RESKH = RESK*0.5D+00
  199.       RESASC = WGK(8)*DABS(FC-RESKH)
  200.       DO 20 J=1,7
  201.         RESASC = RESASC+WGK(J)*(DABS(FV1(J)-RESKH)+DABS(FV2(J)-RESKH))
  202.    20 CONTINUE
  203.       RESULT = RESK*HLGTH
  204.       RESASC = RESASC*HLGTH
  205.       RESABS = RESABS*HLGTH
  206.       ABSERR = DABS((RESK-RESG)*HLGTH)
  207.       IF(RESASC.NE.0.0D+00.AND.ABSERR.NE.0.D0) ABSERR = RESASC*
  208.      * DMIN1(0.1D+01,(0.2D+03*ABSERR/RESASC)**1.5D+00)
  209.       IF(RESABS.GT.UFLOW/(0.5D+02*EPMACH)) ABSERR = DMAX1
  210.      * ((EPMACH*0.5D+02)*RESABS,ABSERR)
  211.       RETURN
  212.       END
  213.