src/numerical/slatec/fortran/dqk15i.f

   1 *DECK DQK15I
   2       SUBROUTINE DQK15I (F, BOUN, INF, A, B, RESULT, ABSERR, RESABS,
   3      +   RESASC)
   4 C***BEGIN PROLOGUE  DQK15I
   5 C***PURPOSE  The original (infinite integration range is mapped
   6 C            onto the interval (0,1) and (A,B) is a part of (0,1).
   7 C            it is the purpose to compute
   8 C            I = Integral of transformed integrand over (A,B),
   9 C            J = Integral of ABS(Transformed Integrand) over (A,B).
  10 C***LIBRARY   SLATEC (QUADPACK)
  11 C***CATEGORY  H2A3A2, H2A4A2
  12 C***TYPE      DOUBLE PRECISION (QK15I-S, DQK15I-D)
  13 C***KEYWORDS  15-POINT GAUSS-KRONROD RULES, QUADPACK, QUADRATURE
  14 C***AUTHOR  Piessens, Robert
  15 C             Applied Mathematics and Programming Division
  16 C             K. U. Leuven
  17 C           de Doncker, Elise
  18 C             Applied Mathematics and Programming Division
  19 C             K. U. Leuven
  20 C***DESCRIPTION
  21 C
  22 C           Integration Rule
  23 C           Standard Fortran subroutine
  24 C           Double precision version
  25 C
  26 C           PARAMETERS
  27 C            ON ENTRY
  28 C              F      - Double precision
  29 C                       Function subprogram defining the integrand
  30 C                       FUNCTION F(X). The actual name for F needs to be
  31 C                       Declared E X T E R N A L in the calling program.
  32 C
  33 C              BOUN   - Double precision
  34 C                       Finite bound of original integration
  35 C                       Range (SET TO ZERO IF INF = +2)
  36 C
  37 C              INF    - Integer
  38 C                       If INF = -1, the original interval is
  39 C                                   (-INFINITY,BOUND),
  40 C                       If INF = +1, the original interval is
  41 C                                   (BOUND,+INFINITY),
  42 C                       If INF = +2, the original interval is
  43 C                                   (-INFINITY,+INFINITY) AND
  44 C                       The integral is computed as the sum of two
  45 C                       integrals, one over (-INFINITY,0) and one over
  46 C                       (0,+INFINITY).
  47 C
  48 C              A      - Double precision
  49 C                       Lower limit for integration over subrange
  50 C                       of (0,1)
  51 C
  52 C              B      - Double precision
  53 C                       Upper limit for integration over subrange
  54 C                       of (0,1)
  55 C
  56 C            ON RETURN
  57 C              RESULT - Double precision
  58 C                       Approximation to the integral I
  59 C                       Result is computed by applying the 15-POINT
  60 C                       KRONROD RULE(RESK) obtained by optimal addition
  61 C                       of abscissae to the 7-POINT GAUSS RULE(RESG).
  62 C
  63 C              ABSERR - Double precision
  64 C                       Estimate of the modulus of the absolute error,
  65 C                       WHICH SHOULD EQUAL or EXCEED ABS(I-RESULT)
  66 C
  67 C              RESABS - Double precision
  68 C                       Approximation to the integral J
  69 C
  70 C              RESASC - Double precision
  71 C                       Approximation to the integral of
  72 C                       ABS((TRANSFORMED INTEGRAND)-I/(B-A)) over (A,B)
  73 C
  74 C***REFERENCES  (NONE)
  75 C***ROUTINES CALLED  D1MACH
  76 C***REVISION HISTORY  (YYMMDD)
  77 C   800101  DATE WRITTEN
  78 C   890531  Changed all specific intrinsics to generic.  (WRB)
  79 C   890531  REVISION DATE from Version 3.2
  80 C   891214  Prologue converted to Version 4.0 format.  (BAB)
  81 C***END PROLOGUE  DQK15I
  82 C
  83       DOUBLE PRECISION A,ABSC,ABSC1,ABSC2,ABSERR,B,BOUN,CENTR,DINF,
  84      1  D1MACH,EPMACH,F,FC,FSUM,FVAL1,FVAL2,FV1,FV2,HLGTH,
  85      2  RESABS,RESASC,RESG,RESK,RESKH,RESULT,TABSC1,TABSC2,UFLOW,WG,WGK,
  86      3  XGK
  87       INTEGER INF,J
  88       EXTERNAL F
  89 C
  90       DIMENSION FV1(7),FV2(7),XGK(8),WGK(8),WG(8)
  91 C
  92 C           THE ABSCISSAE AND WEIGHTS ARE SUPPLIED FOR THE INTERVAL
  93 C           (-1,1).  BECAUSE OF SYMMETRY ONLY THE POSITIVE ABSCISSAE AND
  94 C           THEIR CORRESPONDING WEIGHTS ARE GIVEN.
  95 C
  96 C           XGK    - ABSCISSAE OF THE 15-POINT KRONROD RULE
  97 C                    XGK(2), XGK(4), ... ABSCISSAE OF THE 7-POINT
  98 C                    GAUSS RULE
  99 C                    XGK(1), XGK(3), ...  ABSCISSAE WHICH ARE OPTIMALLY
 100 C                    ADDED TO THE 7-POINT GAUSS RULE
 101 C
 102 C           WGK    - WEIGHTS OF THE 15-POINT KRONROD RULE
 103 C
 104 C           WG     - WEIGHTS OF THE 7-POINT GAUSS RULE, CORRESPONDING
 105 C                    TO THE ABSCISSAE XGK(2), XGK(4), ...
 106 C                    WG(1), WG(3), ... ARE SET TO ZERO.
 107 C
 108       SAVE XGK, WGK, WG
 109       DATA WG(1) / 0.0D0 /
 110       DATA WG(2) / 0.1294849661 6886969327 0611432679 082D0 /
 111       DATA WG(3) / 0.0D0 /
 112       DATA WG(4) / 0.2797053914 8927666790 1467771423 780D0 /
 113       DATA WG(5) / 0.0D0 /
 114       DATA WG(6) / 0.3818300505 0511894495 0369775488 975D0 /
 115       DATA WG(7) / 0.0D0 /
 116       DATA WG(8) / 0.4179591836 7346938775 5102040816 327D0 /
 117 C
 118       DATA XGK(1) / 0.9914553711 2081263920 6854697526 329D0 /
 119       DATA XGK(2) / 0.9491079123 4275852452 6189684047 851D0 /
 120       DATA XGK(3) / 0.8648644233 5976907278 9712788640 926D0 /
 121       DATA XGK(4) / 0.7415311855 9939443986 3864773280 788D0 /
 122       DATA XGK(5) / 0.5860872354 6769113029 4144838258 730D0 /
 123       DATA XGK(6) / 0.4058451513 7739716690 6606412076 961D0 /
 124       DATA XGK(7) / 0.2077849550 0789846760 0689403773 245D0 /
 125       DATA XGK(8) / 0.0000000000 0000000000 0000000000 000D0 /
 126 C
 127       DATA WGK(1) / 0.0229353220 1052922496 3732008058 970D0 /
 128       DATA WGK(2) / 0.0630920926 2997855329 0700663189 204D0 /
 129       DATA WGK(3) / 0.1047900103 2225018383 9876322541 518D0 /
 130       DATA WGK(4) / 0.1406532597 1552591874 5189590510 238D0 /
 131       DATA WGK(5) / 0.1690047266 3926790282 6583426598 550D0 /
 132       DATA WGK(6) / 0.1903505780 6478540991 3256402421 014D0 /
 133       DATA WGK(7) / 0.2044329400 7529889241 4161999234 649D0 /
 134       DATA WGK(8) / 0.2094821410 8472782801 2999174891 714D0 /
 135 C
 136 C
 137 C           LIST OF MAJOR VARIABLES
 138 C           -----------------------
 139 C
 140 C           CENTR  - MID POINT OF THE INTERVAL
 141 C           HLGTH  - HALF-LENGTH OF THE INTERVAL
 142 C           ABSC*  - ABSCISSA
 143 C           TABSC* - TRANSFORMED ABSCISSA
 144 C           FVAL*  - FUNCTION VALUE
 145 C           RESG   - RESULT OF THE 7-POINT GAUSS FORMULA
 146 C           RESK   - RESULT OF THE 15-POINT KRONROD FORMULA
 147 C           RESKH  - APPROXIMATION TO THE MEAN VALUE OF THE TRANSFORMED
 148 C                    INTEGRAND OVER (A,B), I.E. TO I/(B-A)
 149 C
 150 C           MACHINE DEPENDENT CONSTANTS
 151 C           ---------------------------
 152 C
 153 C           EPMACH IS THE LARGEST RELATIVE SPACING.
 154 C           UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
 155 C
 156 C***FIRST EXECUTABLE STATEMENT  DQK15I
 157       EPMACH = D1MACH(4)
 158       UFLOW = D1MACH(1)
 159       DINF = MIN(1,INF)
 160 C
 161       CENTR = 0.5D+00*(A+B)
 162       HLGTH = 0.5D+00*(B-A)
 163       TABSC1 = BOUN+DINF*(0.1D+01-CENTR)/CENTR
 164       FVAL1 = F(TABSC1)
 165       IF(INF.EQ.2) FVAL1 = FVAL1+F(-TABSC1)
 166       FC = (FVAL1/CENTR)/CENTR
 167 C
 168 C           COMPUTE THE 15-POINT KRONROD APPROXIMATION TO
 169 C           THE INTEGRAL, AND ESTIMATE THE ERROR.
 170 C
 171       RESG = WG(8)*FC
 172       RESK = WGK(8)*FC
 173       RESABS = ABS(RESK)
 174       DO 10 J=1,7
 175         ABSC = HLGTH*XGK(J)
 176         ABSC1 = CENTR-ABSC
 177         ABSC2 = CENTR+ABSC
 178         TABSC1 = BOUN+DINF*(0.1D+01-ABSC1)/ABSC1
 179         TABSC2 = BOUN+DINF*(0.1D+01-ABSC2)/ABSC2
 180         FVAL1 = F(TABSC1)
 181         FVAL2 = F(TABSC2)
 182         IF(INF.EQ.2) FVAL1 = FVAL1+F(-TABSC1)
 183         IF(INF.EQ.2) FVAL2 = FVAL2+F(-TABSC2)
 184         FVAL1 = (FVAL1/ABSC1)/ABSC1
 185         FVAL2 = (FVAL2/ABSC2)/ABSC2
 186         FV1(J) = FVAL1
 187         FV2(J) = FVAL2
 188         FSUM = FVAL1+FVAL2
 189         RESG = RESG+WG(J)*FSUM
 190         RESK = RESK+WGK(J)*FSUM
 191         RESABS = RESABS+WGK(J)*(ABS(FVAL1)+ABS(FVAL2))
 192    10 CONTINUE
 193       RESKH = RESK*0.5D+00
 194       RESASC = WGK(8)*ABS(FC-RESKH)
 195       DO 20 J=1,7
 196         RESASC = RESASC+WGK(J)*(ABS(FV1(J)-RESKH)+ABS(FV2(J)-RESKH))
 197    20 CONTINUE
 198       RESULT = RESK*HLGTH
 199       RESASC = RESASC*HLGTH
 200       RESABS = RESABS*HLGTH
 201       ABSERR = ABS((RESK-RESG)*HLGTH)
 202       IF(RESASC.NE.0.0D+00.AND.ABSERR.NE.0.D0) ABSERR = RESASC*
 203      1 MIN(0.1D+01,(0.2D+03*ABSERR/RESASC)**1.5D+00)
 204       IF(RESABS.GT.UFLOW/(0.5D+02*EPMACH)) ABSERR = MAX
 205      1 ((EPMACH*0.5D+02)*RESABS,ABSERR)
 206       RETURN
 207       END