src/numerical/slatec/fortran/dqc25s.f

   1 *DECK DQC25S
   2       SUBROUTINE DQC25S (F, A, B, BL, BR, ALFA, BETA, RI, RJ, RG, RH,
   3      +   RESULT, ABSERR, RESASC, INTEGR, NEV)
   4 C***BEGIN PROLOGUE  DQC25S
   5 C***PURPOSE  To compute I = Integral of F*W over (BL,BR), with error
   6 C            estimate, where the weight function W has a singular
   7 C            behaviour of ALGEBRAICO-LOGARITHMIC type at the points
   8 C            A and/or B. (BL,BR) is a part of (A,B).
   9 C***LIBRARY   SLATEC (QUADPACK)
  10 C***CATEGORY  H2A2A2
  11 C***TYPE      DOUBLE PRECISION (QC25S-S, DQC25S-D)
  12 C***KEYWORDS  25-POINT CLENSHAW-CURTIS INTEGRATION, QUADPACK, QUADRATURE
  13 C***AUTHOR  Piessens, Robert
  14 C             Applied Mathematics and Programming Division
  15 C             K. U. Leuven
  16 C           de Doncker, Elise
  17 C             Applied Mathematics and Programming Division
  18 C             K. U. Leuven
  19 C***DESCRIPTION
  20 C
  21 C        Integration rules for integrands having ALGEBRAICO-LOGARITHMIC
  22 C        end point singularities
  23 C        Standard fortran subroutine
  24 C        Double precision version
  25 C
  26 C        PARAMETERS
  27 C           F      - Double precision
  28 C                    Function subprogram defining the integrand
  29 C                    F(X). The actual name for F needs to be declared
  30 C                    E X T E R N A L  in the driver program.
  31 C
  32 C           A      - Double precision
  33 C                    Left end point of the original interval
  34 C
  35 C           B      - Double precision
  36 C                    Right end point of the original interval, B.GT.A
  37 C
  38 C           BL     - Double precision
  39 C                    Lower limit of integration, BL.GE.A
  40 C
  41 C           BR     - Double precision
  42 C                    Upper limit of integration, BR.LE.B
  43 C
  44 C           ALFA   - Double precision
  45 C                    PARAMETER IN THE WEIGHT FUNCTION
  46 C
  47 C           BETA   - Double precision
  48 C                    Parameter in the weight function
  49 C
  50 C           RI,RJ,RG,RH - Double precision
  51 C                    Modified CHEBYSHEV moments for the application
  52 C                    of the generalized CLENSHAW-CURTIS
  53 C                    method (computed in subroutine DQMOMO)
  54 C
  55 C           RESULT - Double precision
  56 C                    Approximation to the integral
  57 C                    RESULT is computed by using a generalized
  58 C                    CLENSHAW-CURTIS method if B1 = A or BR = B.
  59 C                    in all other cases the 15-POINT KRONROD
  60 C                    RULE is applied, obtained by optimal addition of
  61 C                    Abscissae to the 7-POINT GAUSS RULE.
  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           RESASC - Double precision
  68 C                    Approximation to the integral of ABS(F*W-I/(B-A))
  69 C
  70 C           INTEGR - Integer
  71 C                    Which determines the weight function
  72 C                    = 1   W(X) = (X-A)**ALFA*(B-X)**BETA
  73 C                    = 2   W(X) = (X-A)**ALFA*(B-X)**BETA*LOG(X-A)
  74 C                    = 3   W(X) = (X-A)**ALFA*(B-X)**BETA*LOG(B-X)
  75 C                    = 4   W(X) = (X-A)**ALFA*(B-X)**BETA*LOG(X-A)*
  76 C                                 LOG(B-X)
  77 C
  78 C           NEV    - Integer
  79 C                    Number of integrand evaluations
  80 C
  81 C***REFERENCES  (NONE)
  82 C***ROUTINES CALLED  DQCHEB, DQK15W, DQWGTS
  83 C***REVISION HISTORY  (YYMMDD)
  84 C   810101  DATE WRITTEN
  85 C   890531  Changed all specific intrinsics to generic.  (WRB)
  86 C   890531  REVISION DATE from Version 3.2
  87 C   891214  Prologue converted to Version 4.0 format.  (BAB)
  88 C***END PROLOGUE  DQC25S
  89 C
  90       DOUBLE PRECISION A,ABSERR,ALFA,B,BETA,BL,BR,CENTR,CHEB12,CHEB24,
  91      1  DC,F,FACTOR,FIX,FVAL,HLGTH,RESABS,RESASC,RESULT,RES12,
  92      2  RES24,RG,RH,RI,RJ,U,DQWGTS,X
  93       INTEGER I,INTEGR,ISYM,NEV
  94 C
  95       DIMENSION CHEB12(13),CHEB24(25),FVAL(25),RG(25),RH(25),RI(25),
  96      1  RJ(25),X(11)
  97 C
  98       EXTERNAL F, DQWGTS
  99 C
 100 C           THE VECTOR X CONTAINS THE VALUES COS(K*PI/24)
 101 C           K = 1, ..., 11, TO BE USED FOR THE COMPUTATION OF THE
 102 C           CHEBYSHEV SERIES EXPANSION OF F.
 103 C
 104       SAVE X
 105        DATA X(1),X(2),X(3),X(4),X(5),X(6),X(7),X(8),X(9),X(10),X(11)/
 106      1     0.9914448613738104D+00,     0.9659258262890683D+00,
 107      2     0.9238795325112868D+00,     0.8660254037844386D+00,
 108      3     0.7933533402912352D+00,     0.7071067811865475D+00,
 109      4     0.6087614290087206D+00,     0.5000000000000000D+00,
 110      5     0.3826834323650898D+00,     0.2588190451025208D+00,
 111      6     0.1305261922200516D+00/
 112 C
 113 C           LIST OF MAJOR VARIABLES
 114 C           -----------------------
 115 C
 116 C           FVAL   - VALUE OF THE FUNCTION F AT THE POINTS
 117 C                    (BR-BL)*0.5*COS(K*PI/24)+(BR+BL)*0.5
 118 C                    K = 0, ..., 24
 119 C           CHEB12 - COEFFICIENTS OF THE CHEBYSHEV SERIES EXPANSION
 120 C                    OF DEGREE 12, FOR THE FUNCTION F, IN THE
 121 C                    INTERVAL (BL,BR)
 122 C           CHEB24 - COEFFICIENTS OF THE CHEBYSHEV SERIES EXPANSION
 123 C                    OF DEGREE 24, FOR THE FUNCTION F, IN THE
 124 C                    INTERVAL (BL,BR)
 125 C           RES12  - APPROXIMATION TO THE INTEGRAL OBTAINED FROM CHEB12
 126 C           RES24  - APPROXIMATION TO THE INTEGRAL OBTAINED FROM CHEB24
 127 C           DQWGTS - EXTERNAL FUNCTION SUBPROGRAM DEFINING
 128 C                    THE FOUR POSSIBLE WEIGHT FUNCTIONS
 129 C           HLGTH  - HALF-LENGTH OF THE INTERVAL (BL,BR)
 130 C           CENTR  - MID POINT OF THE INTERVAL (BL,BR)
 131 C
 132 C***FIRST EXECUTABLE STATEMENT  DQC25S
 133       NEV = 25
 134       IF(BL.EQ.A.AND.(ALFA.NE.0.0D+00.OR.INTEGR.EQ.2.OR.INTEGR.EQ.4))
 135      1 GO TO 10
 136       IF(BR.EQ.B.AND.(BETA.NE.0.0D+00.OR.INTEGR.EQ.3.OR.INTEGR.EQ.4))
 137      1 GO TO 140
 138 C
 139 C           IF A.GT.BL AND B.LT.BR, APPLY THE 15-POINT GAUSS-KRONROD
 140 C           SCHEME.
 141 C
 142 C
 143       CALL DQK15W(F,DQWGTS,A,B,ALFA,BETA,INTEGR,BL,BR,
 144      1    RESULT,ABSERR,RESABS,RESASC)
 145       NEV = 15
 146       GO TO 270
 147 C
 148 C           THIS PART OF THE PROGRAM IS EXECUTED ONLY IF A = BL.
 149 C           ----------------------------------------------------
 150 C
 151 C           COMPUTE THE CHEBYSHEV SERIES EXPANSION OF THE
 152 C           FOLLOWING FUNCTION
 153 C           F1 = (0.5*(B+B-BR-A)-0.5*(BR-A)*X)**BETA
 154 C                  *F(0.5*(BR-A)*X+0.5*(BR+A))
 155 C
 156    10 HLGTH = 0.5D+00*(BR-BL)
 157       CENTR = 0.5D+00*(BR+BL)
 158       FIX = B-CENTR
 159       FVAL(1) = 0.5D+00*F(HLGTH+CENTR)*(FIX-HLGTH)**BETA
 160       FVAL(13) = F(CENTR)*(FIX**BETA)
 161       FVAL(25) = 0.5D+00*F(CENTR-HLGTH)*(FIX+HLGTH)**BETA
 162       DO 20 I=2,12
 163         U = HLGTH*X(I-1)
 164         ISYM = 26-I
 165         FVAL(I) = F(U+CENTR)*(FIX-U)**BETA
 166         FVAL(ISYM) = F(CENTR-U)*(FIX+U)**BETA
 167    20 CONTINUE
 168       FACTOR = HLGTH**(ALFA+0.1D+01)
 169       RESULT = 0.0D+00
 170       ABSERR = 0.0D+00
 171       RES12 = 0.0D+00
 172       RES24 = 0.0D+00
 173       IF(INTEGR.GT.2) GO TO 70
 174       CALL DQCHEB(X,FVAL,CHEB12,CHEB24)
 175 C
 176 C           INTEGR = 1  (OR 2)
 177 C
 178       DO 30 I=1,13
 179         RES12 = RES12+CHEB12(I)*RI(I)
 180         RES24 = RES24+CHEB24(I)*RI(I)
 181    30 CONTINUE
 182       DO 40 I=14,25
 183         RES24 = RES24+CHEB24(I)*RI(I)
 184    40 CONTINUE
 185       IF(INTEGR.EQ.1) GO TO 130
 186 C
 187 C           INTEGR = 2
 188 C
 189       DC = LOG(BR-BL)
 190       RESULT = RES24*DC
 191       ABSERR = ABS((RES24-RES12)*DC)
 192       RES12 = 0.0D+00
 193       RES24 = 0.0D+00
 194       DO 50 I=1,13
 195         RES12 = RES12+CHEB12(I)*RG(I)
 196         RES24 = RES24+CHEB24(I)*RG(I)
 197    50 CONTINUE
 198       DO 60 I=14,25
 199         RES24 = RES24+CHEB24(I)*RG(I)
 200    60 CONTINUE
 201       GO TO 130
 202 C
 203 C           COMPUTE THE CHEBYSHEV SERIES EXPANSION OF THE
 204 C           FOLLOWING FUNCTION
 205 C           F4 = F1*LOG(0.5*(B+B-BR-A)-0.5*(BR-A)*X)
 206 C
 207    70 FVAL(1) = FVAL(1)*LOG(FIX-HLGTH)
 208       FVAL(13) = FVAL(13)*LOG(FIX)
 209       FVAL(25) = FVAL(25)*LOG(FIX+HLGTH)
 210       DO 80 I=2,12
 211         U = HLGTH*X(I-1)
 212         ISYM = 26-I
 213         FVAL(I) = FVAL(I)*LOG(FIX-U)
 214         FVAL(ISYM) = FVAL(ISYM)*LOG(FIX+U)
 215    80 CONTINUE
 216       CALL DQCHEB(X,FVAL,CHEB12,CHEB24)
 217 C
 218 C           INTEGR = 3  (OR 4)
 219 C
 220       DO 90 I=1,13
 221         RES12 = RES12+CHEB12(I)*RI(I)
 222         RES24 = RES24+CHEB24(I)*RI(I)
 223    90 CONTINUE
 224       DO 100 I=14,25
 225         RES24 = RES24+CHEB24(I)*RI(I)
 226   100 CONTINUE
 227       IF(INTEGR.EQ.3) GO TO 130
 228 C
 229 C           INTEGR = 4
 230 C
 231       DC = LOG(BR-BL)
 232       RESULT = RES24*DC
 233       ABSERR = ABS((RES24-RES12)*DC)
 234       RES12 = 0.0D+00
 235       RES24 = 0.0D+00
 236       DO 110 I=1,13
 237         RES12 = RES12+CHEB12(I)*RG(I)
 238         RES24 = RES24+CHEB24(I)*RG(I)
 239   110 CONTINUE
 240       DO 120 I=14,25
 241         RES24 = RES24+CHEB24(I)*RG(I)
 242   120 CONTINUE
 243   130 RESULT = (RESULT+RES24)*FACTOR
 244       ABSERR = (ABSERR+ABS(RES24-RES12))*FACTOR
 245       GO TO 270
 246 C
 247 C           THIS PART OF THE PROGRAM IS EXECUTED ONLY IF B = BR.
 248 C           ----------------------------------------------------
 249 C
 250 C           COMPUTE THE CHEBYSHEV SERIES EXPANSION OF THE
 251 C           FOLLOWING FUNCTION
 252 C           F2 = (0.5*(B+BL-A-A)+0.5*(B-BL)*X)**ALFA
 253 C                *F(0.5*(B-BL)*X+0.5*(B+BL))
 254 C
 255   140 HLGTH = 0.5D+00*(BR-BL)
 256       CENTR = 0.5D+00*(BR+BL)
 257       FIX = CENTR-A
 258       FVAL(1) = 0.5D+00*F(HLGTH+CENTR)*(FIX+HLGTH)**ALFA
 259       FVAL(13) = F(CENTR)*(FIX**ALFA)
 260       FVAL(25) = 0.5D+00*F(CENTR-HLGTH)*(FIX-HLGTH)**ALFA
 261       DO 150 I=2,12
 262         U = HLGTH*X(I-1)
 263         ISYM = 26-I
 264         FVAL(I) = F(U+CENTR)*(FIX+U)**ALFA
 265         FVAL(ISYM) = F(CENTR-U)*(FIX-U)**ALFA
 266   150 CONTINUE
 267       FACTOR = HLGTH**(BETA+0.1D+01)
 268       RESULT = 0.0D+00
 269       ABSERR = 0.0D+00
 270       RES12 = 0.0D+00
 271       RES24 = 0.0D+00
 272       IF(INTEGR.EQ.2.OR.INTEGR.EQ.4) GO TO 200
 273 C
 274 C           INTEGR = 1  (OR 3)
 275 C
 276       CALL DQCHEB(X,FVAL,CHEB12,CHEB24)
 277       DO 160 I=1,13
 278         RES12 = RES12+CHEB12(I)*RJ(I)
 279         RES24 = RES24+CHEB24(I)*RJ(I)
 280   160 CONTINUE
 281       DO 170 I=14,25
 282         RES24 = RES24+CHEB24(I)*RJ(I)
 283   170 CONTINUE
 284       IF(INTEGR.EQ.1) GO TO 260
 285 C
 286 C           INTEGR = 3
 287 C
 288       DC = LOG(BR-BL)
 289       RESULT = RES24*DC
 290       ABSERR = ABS((RES24-RES12)*DC)
 291       RES12 = 0.0D+00
 292       RES24 = 0.0D+00
 293       DO 180 I=1,13
 294         RES12 = RES12+CHEB12(I)*RH(I)
 295         RES24 = RES24+CHEB24(I)*RH(I)
 296   180 CONTINUE
 297       DO 190 I=14,25
 298         RES24 = RES24+CHEB24(I)*RH(I)
 299   190 CONTINUE
 300       GO TO 260
 301 C
 302 C           COMPUTE THE CHEBYSHEV SERIES EXPANSION OF THE
 303 C           FOLLOWING FUNCTION
 304 C           F3 = F2*LOG(0.5*(B-BL)*X+0.5*(B+BL-A-A))
 305 C
 306   200 FVAL(1) = FVAL(1)*LOG(FIX+HLGTH)
 307       FVAL(13) = FVAL(13)*LOG(FIX)
 308       FVAL(25) = FVAL(25)*LOG(FIX-HLGTH)
 309       DO 210 I=2,12
 310         U = HLGTH*X(I-1)
 311         ISYM = 26-I
 312         FVAL(I) = FVAL(I)*LOG(U+FIX)
 313         FVAL(ISYM) = FVAL(ISYM)*LOG(FIX-U)
 314   210 CONTINUE
 315       CALL DQCHEB(X,FVAL,CHEB12,CHEB24)
 316 C
 317 C           INTEGR = 2  (OR 4)
 318 C
 319       DO 220 I=1,13
 320         RES12 = RES12+CHEB12(I)*RJ(I)
 321         RES24 = RES24+CHEB24(I)*RJ(I)
 322   220 CONTINUE
 323       DO 230 I=14,25
 324         RES24 = RES24+CHEB24(I)*RJ(I)
 325   230 CONTINUE
 326       IF(INTEGR.EQ.2) GO TO 260
 327       DC = LOG(BR-BL)
 328       RESULT = RES24*DC
 329       ABSERR = ABS((RES24-RES12)*DC)
 330       RES12 = 0.0D+00
 331       RES24 = 0.0D+00
 332 C
 333 C           INTEGR = 4
 334 C
 335       DO 240 I=1,13
 336         RES12 = RES12+CHEB12(I)*RH(I)
 337         RES24 = RES24+CHEB24(I)*RH(I)
 338   240 CONTINUE
 339       DO 250 I=14,25
 340         RES24 = RES24+CHEB24(I)*RH(I)
 341   250 CONTINUE
 342   260 RESULT = (RESULT+RES24)*FACTOR
 343       ABSERR = (ABSERR+ABS(RES24-RES12))*FACTOR
 344   270 RETURN
 345       END