src/numerical/slatec/fortran/dqc25f.f

   1 *DECK DQC25F
   2       SUBROUTINE DQC25F (F, A, B, OMEGA, INTEGR, NRMOM, MAXP1, KSAVE,
   3      +   RESULT, ABSERR, NEVAL, RESABS, RESASC, MOMCOM, CHEBMO)
   4 C***BEGIN PROLOGUE  DQC25F
   5 C***PURPOSE  To compute the integral I=Integral of F(X) over (A,B)
   6 C            Where W(X) = COS(OMEGA*X) or W(X)=SIN(OMEGA*X) and to
   7 C            compute J = Integral of ABS(F) over (A,B). For small value
   8 C            of OMEGA or small intervals (A,B) the 15-point GAUSS-KRONRO
   9 C            Rule is used. Otherwise a generalized CLENSHAW-CURTIS
  10 C            method is used.
  11 C***LIBRARY   SLATEC (QUADPACK)
  12 C***CATEGORY  H2A2A2
  13 C***TYPE      DOUBLE PRECISION (QC25F-S, DQC25F-D)
  14 C***KEYWORDS  CLENSHAW-CURTIS METHOD, GAUSS-KRONROD RULES,
  15 C             INTEGRATION RULES FOR FUNCTIONS WITH COS OR SIN FACTOR,
  16 C             QUADPACK, QUADRATURE
  17 C***AUTHOR  Piessens, Robert
  18 C             Applied Mathematics and Programming Division
  19 C             K. U. Leuven
  20 C           de Doncker, Elise
  21 C             Applied Mathematics and Programming Division
  22 C             K. U. Leuven
  23 C***DESCRIPTION
  24 C
  25 C        Integration rules for functions with COS or SIN factor
  26 C        Standard fortran subroutine
  27 C        Double precision version
  28 C
  29 C        PARAMETERS
  30 C         ON ENTRY
  31 C           F      - Double precision
  32 C                    Function subprogram defining the integrand
  33 C                    function F(X). The actual name for F needs to
  34 C                    be declared E X T E R N A L in the calling program.
  35 C
  36 C           A      - Double precision
  37 C                    Lower limit of integration
  38 C
  39 C           B      - Double precision
  40 C                    Upper limit of integration
  41 C
  42 C           OMEGA  - Double precision
  43 C                    Parameter in the WEIGHT function
  44 C
  45 C           INTEGR - Integer
  46 C                    Indicates which WEIGHT function is to be used
  47 C                       INTEGR = 1   W(X) = COS(OMEGA*X)
  48 C                       INTEGR = 2   W(X) = SIN(OMEGA*X)
  49 C
  50 C           NRMOM  - Integer
  51 C                    The length of interval (A,B) is equal to the length
  52 C                    of the original integration interval divided by
  53 C                    2**NRMOM (we suppose that the routine is used in an
  54 C                    adaptive integration process, otherwise set
  55 C                    NRMOM = 0). NRMOM must be zero at the first call.
  56 C
  57 C           MAXP1  - Integer
  58 C                    Gives an upper bound on the number of Chebyshev
  59 C                    moments which can be stored, i.e. for the
  60 C                    intervals of lengths ABS(BB-AA)*2**(-L),
  61 C                    L = 0,1,2, ..., MAXP1-2.
  62 C
  63 C           KSAVE  - Integer
  64 C                    Key which is one when the moments for the
  65 C                    current interval have been computed
  66 C
  67 C         ON RETURN
  68 C           RESULT - Double precision
  69 C                    Approximation to the integral I
  70 C
  71 C           ABSERR - Double precision
  72 C                    Estimate of the modulus of the absolute
  73 C                    error, which should equal or exceed ABS(I-RESULT)
  74 C
  75 C           NEVAL  - Integer
  76 C                    Number of integrand evaluations
  77 C
  78 C           RESABS - Double precision
  79 C                    Approximation to the integral J
  80 C
  81 C           RESASC - Double precision
  82 C                    Approximation to the integral of ABS(F-I/(B-A))
  83 C
  84 C         ON ENTRY AND RETURN
  85 C           MOMCOM - Integer
  86 C                    For each interval length we need to compute the
  87 C                    Chebyshev moments. MOMCOM counts the number of
  88 C                    intervals for which these moments have already been
  89 C                    computed. If NRMOM.LT.MOMCOM or KSAVE = 1, the
  90 C                    Chebyshev moments for the interval (A,B) have
  91 C                    already been computed and stored, otherwise we
  92 C                    compute them and we increase MOMCOM.
  93 C
  94 C           CHEBMO - Double precision
  95 C                    Array of dimension at least (MAXP1,25) containing
  96 C                    the modified Chebyshev moments for the first MOMCOM
  97 C                    MOMCOM interval lengths
  98 C
  99 C ......................................................................
 100 C
 101 C***REFERENCES  (NONE)
 102 C***ROUTINES CALLED  D1MACH, DGTSL, DQCHEB, DQK15W, DQWGTF
 103 C***REVISION HISTORY  (YYMMDD)
 104 C   810101  DATE WRITTEN
 105 C   890531  Changed all specific intrinsics to generic.  (WRB)
 106 C   890531  REVISION DATE from Version 3.2
 107 C   891214  Prologue converted to Version 4.0 format.  (BAB)
 108 C***END PROLOGUE  DQC25F
 109 C
 110       DOUBLE PRECISION A,ABSERR,AC,AN,AN2,AS,ASAP,ASS,B,CENTR,CHEBMO,
 111      1  CHEB12,CHEB24,CONC,CONS,COSPAR,D,DQWGTF,D1,
 112      2  D1MACH,D2,ESTC,ESTS,F,FVAL,HLGTH,OFLOW,OMEGA,PARINT,PAR2,PAR22,
 113      3  P2,P3,P4,RESABS,RESASC,RESC12,RESC24,RESS12,RESS24,RESULT,
 114      4  SINPAR,V,X
 115       INTEGER I,IERS,INTEGR,ISYM,J,K,KSAVE,M,MOMCOM,NEVAL,MAXP1,
 116      1  NOEQU,NOEQ1,NRMOM
 117 C
 118       DIMENSION CHEBMO(MAXP1,25),CHEB12(13),CHEB24(25),D(25),D1(25),
 119      1  D2(25),FVAL(25),V(28),X(11)
 120 C
 121       EXTERNAL F, DQWGTF
 122 C
 123 C           THE VECTOR X CONTAINS THE VALUES COS(K*PI/24)
 124 C           K = 1, ...,11, TO BE USED FOR THE CHEBYSHEV EXPANSION OF F
 125 C
 126       SAVE X
 127       DATA X(1) / 0.9914448613 7381041114 4557526928 563D0 /
 128       DATA X(2) / 0.9659258262 8906828674 9743199728 897D0 /
 129       DATA X(3) / 0.9238795325 1128675612 8183189396 788D0 /
 130       DATA X(4) / 0.8660254037 8443864676 3723170752 936D0 /
 131       DATA X(5) / 0.7933533402 9123516457 9776961501 299D0 /
 132       DATA X(6) / 0.7071067811 8654752440 0844362104 849D0 /
 133       DATA X(7) / 0.6087614290 0872063941 6097542898 164D0 /
 134       DATA X(8) / 0.5000000000 0000000000 0000000000 000D0 /
 135       DATA X(9) / 0.3826834323 6508977172 8459984030 399D0 /
 136       DATA X(10) / 0.2588190451 0252076234 8898837624 048D0 /
 137       DATA X(11) / 0.1305261922 2005159154 8406227895 489D0 /
 138 C
 139 C           LIST OF MAJOR VARIABLES
 140 C           -----------------------
 141 C
 142 C           CENTR  - MID POINT OF THE INTEGRATION INTERVAL
 143 C           HLGTH  - HALF-LENGTH OF THE INTEGRATION INTERVAL
 144 C           FVAL   - VALUE OF THE FUNCTION F AT THE POINTS
 145 C                    (B-A)*0.5*COS(K*PI/12) + (B+A)*0.5, K = 0, ..., 24
 146 C           CHEB12 - COEFFICIENTS OF THE CHEBYSHEV SERIES EXPANSION
 147 C                    OF DEGREE 12, FOR THE FUNCTION F, IN THE
 148 C                    INTERVAL (A,B)
 149 C           CHEB24 - COEFFICIENTS OF THE CHEBYSHEV SERIES EXPANSION
 150 C                    OF DEGREE 24, FOR THE FUNCTION F, IN THE
 151 C                    INTERVAL (A,B)
 152 C           RESC12 - APPROXIMATION TO THE INTEGRAL OF
 153 C                    COS(0.5*(B-A)*OMEGA*X)*F(0.5*(B-A)*X+0.5*(B+A))
 154 C                    OVER (-1,+1), USING THE CHEBYSHEV SERIES
 155 C                    EXPANSION OF DEGREE 12
 156 C           RESC24 - APPROXIMATION TO THE SAME INTEGRAL, USING THE
 157 C                    CHEBYSHEV SERIES EXPANSION OF DEGREE 24
 158 C           RESS12 - THE ANALOGUE OF RESC12 FOR THE SINE
 159 C           RESS24 - THE ANALOGUE OF RESC24 FOR THE SINE
 160 C
 161 C
 162 C           MACHINE DEPENDENT CONSTANT
 163 C           --------------------------
 164 C
 165 C           OFLOW IS THE LARGEST POSITIVE MAGNITUDE.
 166 C
 167 C***FIRST EXECUTABLE STATEMENT  DQC25F
 168       OFLOW = D1MACH(2)
 169 C
 170       CENTR = 0.5D+00*(B+A)
 171       HLGTH = 0.5D+00*(B-A)
 172       PARINT = OMEGA*HLGTH
 173 C
 174 C           COMPUTE THE INTEGRAL USING THE 15-POINT GAUSS-KRONROD
 175 C           FORMULA IF THE VALUE OF THE PARAMETER IN THE INTEGRAND
 176 C           IS SMALL.
 177 C
 178       IF(ABS(PARINT).GT.0.2D+01) GO TO 10
 179       CALL DQK15W(F,DQWGTF,OMEGA,P2,P3,P4,INTEGR,A,B,RESULT,
 180      1  ABSERR,RESABS,RESASC)
 181       NEVAL = 15
 182       GO TO 170
 183 C
 184 C           COMPUTE THE INTEGRAL USING THE GENERALIZED CLENSHAW-
 185 C           CURTIS METHOD.
 186 C
 187    10 CONC = HLGTH*COS(CENTR*OMEGA)
 188       CONS = HLGTH*SIN(CENTR*OMEGA)
 189       RESASC = OFLOW
 190       NEVAL = 25
 191 C
 192 C           CHECK WHETHER THE CHEBYSHEV MOMENTS FOR THIS INTERVAL
 193 C           HAVE ALREADY BEEN COMPUTED.
 194 C
 195       IF(NRMOM.LT.MOMCOM.OR.KSAVE.EQ.1) GO TO 120
 196 C
 197 C           COMPUTE A NEW SET OF CHEBYSHEV MOMENTS.
 198 C
 199       M = MOMCOM+1
 200       PAR2 = PARINT*PARINT
 201       PAR22 = PAR2+0.2D+01
 202       SINPAR = SIN(PARINT)
 203       COSPAR = COS(PARINT)
 204 C
 205 C           COMPUTE THE CHEBYSHEV MOMENTS WITH RESPECT TO COSINE.
 206 C
 207       V(1) = 0.2D+01*SINPAR/PARINT
 208       V(2) = (0.8D+01*COSPAR+(PAR2+PAR2-0.8D+01)*SINPAR/PARINT)/PAR2
 209       V(3) = (0.32D+02*(PAR2-0.12D+02)*COSPAR+(0.2D+01*
 210      1  ((PAR2-0.80D+02)*PAR2+0.192D+03)*SINPAR)/PARINT)/(PAR2*PAR2)
 211       AC = 0.8D+01*COSPAR
 212       AS = 0.24D+02*PARINT*SINPAR
 213       IF(ABS(PARINT).GT.0.24D+02) GO TO 30
 214 C
 215 C           COMPUTE THE CHEBYSHEV MOMENTS AS THE SOLUTIONS OF A
 216 C           BOUNDARY VALUE PROBLEM WITH 1 INITIAL VALUE (V(3)) AND 1
 217 C           END VALUE (COMPUTED USING AN ASYMPTOTIC FORMULA).
 218 C
 219       NOEQU = 25
 220       NOEQ1 = NOEQU-1
 221       AN = 0.6D+01
 222       DO 20 K = 1,NOEQ1
 223         AN2 = AN*AN
 224         D(K) = -0.2D+01*(AN2-0.4D+01)*(PAR22-AN2-AN2)
 225         D2(K) = (AN-0.1D+01)*(AN-0.2D+01)*PAR2
 226         D1(K+1) = (AN+0.3D+01)*(AN+0.4D+01)*PAR2
 227         V(K+3) = AS-(AN2-0.4D+01)*AC
 228         AN = AN+0.2D+01
 229    20 CONTINUE
 230       AN2 = AN*AN
 231       D(NOEQU) = -0.2D+01*(AN2-0.4D+01)*(PAR22-AN2-AN2)
 232       V(NOEQU+3) = AS-(AN2-0.4D+01)*AC
 233       V(4) = V(4)-0.56D+02*PAR2*V(3)
 234       ASS = PARINT*SINPAR
 235       ASAP = (((((0.210D+03*PAR2-0.1D+01)*COSPAR-(0.105D+03*PAR2
 236      1  -0.63D+02)*ASS)/AN2-(0.1D+01-0.15D+02*PAR2)*COSPAR
 237      2  +0.15D+02*ASS)/AN2-COSPAR+0.3D+01*ASS)/AN2-COSPAR)/AN2
 238       V(NOEQU+3) = V(NOEQU+3)-0.2D+01*ASAP*PAR2*(AN-0.1D+01)*
 239      1   (AN-0.2D+01)
 240 C
 241 C           SOLVE THE TRIDIAGONAL SYSTEM BY MEANS OF GAUSSIAN
 242 C           ELIMINATION WITH PARTIAL PIVOTING.
 243 C
 244 C ***       CALL TO DGTSL MUST BE REPLACED BY CALL TO
 245 C ***       DOUBLE PRECISION VERSION OF LINPACK ROUTINE SGTSL
 246 C
 247       CALL DGTSL(NOEQU,D1,D,D2,V(4),IERS)
 248       GO TO 50
 249 C
 250 C           COMPUTE THE CHEBYSHEV MOMENTS BY MEANS OF FORWARD
 251 C           RECURSION.
 252 C
 253    30 AN = 0.4D+01
 254       DO 40 I = 4,13
 255         AN2 = AN*AN
 256         V(I) = ((AN2-0.4D+01)*(0.2D+01*(PAR22-AN2-AN2)*V(I-1)-AC)
 257      1  +AS-PAR2*(AN+0.1D+01)*(AN+0.2D+01)*V(I-2))/
 258      2  (PAR2*(AN-0.1D+01)*(AN-0.2D+01))
 259         AN = AN+0.2D+01
 260    40 CONTINUE
 261    50 DO 60 J = 1,13
 262         CHEBMO(M,2*J-1) = V(J)
 263    60 CONTINUE
 264 C
 265 C           COMPUTE THE CHEBYSHEV MOMENTS WITH RESPECT TO SINE.
 266 C
 267       V(1) = 0.2D+01*(SINPAR-PARINT*COSPAR)/PAR2
 268       V(2) = (0.18D+02-0.48D+02/PAR2)*SINPAR/PAR2
 269      1  +(-0.2D+01+0.48D+02/PAR2)*COSPAR/PARINT
 270       AC = -0.24D+02*PARINT*COSPAR
 271       AS = -0.8D+01*SINPAR
 272       IF(ABS(PARINT).GT.0.24D+02) GO TO 80
 273 C
 274 C           COMPUTE THE CHEBYSHEV MOMENTS AS THE SOLUTIONS OF A BOUNDARY
 275 C           VALUE PROBLEM WITH 1 INITIAL VALUE (V(2)) AND 1 END VALUE
 276 C           (COMPUTED USING AN ASYMPTOTIC FORMULA).
 277 C
 278       AN = 0.5D+01
 279       DO 70 K = 1,NOEQ1
 280         AN2 = AN*AN
 281         D(K) = -0.2D+01*(AN2-0.4D+01)*(PAR22-AN2-AN2)
 282         D2(K) = (AN-0.1D+01)*(AN-0.2D+01)*PAR2
 283         D1(K+1) = (AN+0.3D+01)*(AN+0.4D+01)*PAR2
 284         V(K+2) = AC+(AN2-0.4D+01)*AS
 285         AN = AN+0.2D+01
 286    70 CONTINUE
 287       AN2 = AN*AN
 288       D(NOEQU) = -0.2D+01*(AN2-0.4D+01)*(PAR22-AN2-AN2)
 289       V(NOEQU+2) = AC+(AN2-0.4D+01)*AS
 290       V(3) = V(3)-0.42D+02*PAR2*V(2)
 291       ASS = PARINT*COSPAR
 292       ASAP = (((((0.105D+03*PAR2-0.63D+02)*ASS+(0.210D+03*PAR2
 293      1  -0.1D+01)*SINPAR)/AN2+(0.15D+02*PAR2-0.1D+01)*SINPAR-
 294      2  0.15D+02*ASS)/AN2-0.3D+01*ASS-SINPAR)/AN2-SINPAR)/AN2
 295       V(NOEQU+2) = V(NOEQU+2)-0.2D+01*ASAP*PAR2*(AN-0.1D+01)
 296      1  *(AN-0.2D+01)
 297 C
 298 C           SOLVE THE TRIDIAGONAL SYSTEM BY MEANS OF GAUSSIAN
 299 C           ELIMINATION WITH PARTIAL PIVOTING.
 300 C
 301 C ***       CALL TO DGTSL MUST BE REPLACED BY CALL TO
 302 C ***       DOUBLE PRECISION VERSION OF LINPACK ROUTINE SGTSL
 303 C
 304       CALL DGTSL(NOEQU,D1,D,D2,V(3),IERS)
 305       GO TO 100
 306 C
 307 C           COMPUTE THE CHEBYSHEV MOMENTS BY MEANS OF FORWARD RECURSION.
 308 C
 309    80 AN = 0.3D+01
 310       DO 90 I = 3,12
 311         AN2 = AN*AN
 312         V(I) = ((AN2-0.4D+01)*(0.2D+01*(PAR22-AN2-AN2)*V(I-1)+AS)
 313      1  +AC-PAR2*(AN+0.1D+01)*(AN+0.2D+01)*V(I-2))
 314      2  /(PAR2*(AN-0.1D+01)*(AN-0.2D+01))
 315         AN = AN+0.2D+01
 316    90 CONTINUE
 317   100 DO 110 J = 1,12
 318         CHEBMO(M,2*J) = V(J)
 319   110 CONTINUE
 320   120 IF (NRMOM.LT.MOMCOM) M = NRMOM+1
 321        IF (MOMCOM.LT.(MAXP1-1).AND.NRMOM.GE.MOMCOM) MOMCOM = MOMCOM+1
 322 C
 323 C           COMPUTE THE COEFFICIENTS OF THE CHEBYSHEV EXPANSIONS
 324 C           OF DEGREES 12 AND 24 OF THE FUNCTION F.
 325 C
 326       FVAL(1) = 0.5D+00*F(CENTR+HLGTH)
 327       FVAL(13) = F(CENTR)
 328       FVAL(25) = 0.5D+00*F(CENTR-HLGTH)
 329       DO 130 I = 2,12
 330         ISYM = 26-I
 331         FVAL(I) = F(HLGTH*X(I-1)+CENTR)
 332         FVAL(ISYM) = F(CENTR-HLGTH*X(I-1))
 333   130 CONTINUE
 334       CALL DQCHEB(X,FVAL,CHEB12,CHEB24)
 335 C
 336 C           COMPUTE THE INTEGRAL AND ERROR ESTIMATES.
 337 C
 338       RESC12 = CHEB12(13)*CHEBMO(M,13)
 339       RESS12 = 0.0D+00
 340       K = 11
 341       DO 140 J = 1,6
 342         RESC12 = RESC12+CHEB12(K)*CHEBMO(M,K)
 343         RESS12 = RESS12+CHEB12(K+1)*CHEBMO(M,K+1)
 344         K = K-2
 345   140 CONTINUE
 346       RESC24 = CHEB24(25)*CHEBMO(M,25)
 347       RESS24 = 0.0D+00
 348       RESABS = ABS(CHEB24(25))
 349       K = 23
 350       DO 150 J = 1,12
 351         RESC24 = RESC24+CHEB24(K)*CHEBMO(M,K)
 352         RESS24 = RESS24+CHEB24(K+1)*CHEBMO(M,K+1)
 353         RESABS = RESABS + ABS(CHEB24(K))+ABS(CHEB24(K+1))
 354         K = K-2
 355   150 CONTINUE
 356       ESTC = ABS(RESC24-RESC12)
 357       ESTS = ABS(RESS24-RESS12)
 358       RESABS = RESABS*ABS(HLGTH)
 359       IF(INTEGR.EQ.2) GO TO 160
 360       RESULT = CONC*RESC24-CONS*RESS24
 361       ABSERR = ABS(CONC*ESTC)+ABS(CONS*ESTS)
 362       GO TO 170
 363   160 RESULT = CONC*RESS24+CONS*RESC24
 364       ABSERR = ABS(CONC*ESTS)+ABS(CONS*ESTC)
 365   170 RETURN
 366       END