1*DECK DQC25S
2      SUBROUTINE DQC25S (F, A, B, BL, BR, ALFA, BETA, RI, RJ, RG, RH,
3     +   RESULT, ABSERR, RESASC, INTEGR, NEV)
4C***BEGIN PROLOGUE  DQC25S
5C***PURPOSE  To compute I = Integral of F*W over (BL,BR), with error
6C            estimate, where the weight function W has a singular
7C            behaviour of ALGEBRAICO-LOGARITHMIC type at the points
8C            A and/or B. (BL,BR) is a part of (A,B).
9C***LIBRARY   SLATEC (QUADPACK)
10C***CATEGORY  H2A2A2
11C***TYPE      DOUBLE PRECISION (QC25S-S, DQC25S-D)
12C***KEYWORDS  25-POINT CLENSHAW-CURTIS INTEGRATION, QUADPACK, QUADRATURE
13C***AUTHOR  Piessens, Robert
14C             Applied Mathematics and Programming Division
15C             K. U. Leuven
16C           de Doncker, Elise
17C             Applied Mathematics and Programming Division
18C             K. U. Leuven
19C***DESCRIPTION
20C
21C        Integration rules for integrands having ALGEBRAICO-LOGARITHMIC
22C        end point singularities
23C        Standard fortran subroutine
24C        Double precision version
25C
26C        PARAMETERS
27C           F      - Double precision
28C                    Function subprogram defining the integrand
29C                    F(X). The actual name for F needs to be declared
30C                    E X T E R N A L  in the driver program.
31C
32C           A      - Double precision
33C                    Left end point of the original interval
34C
35C           B      - Double precision
36C                    Right end point of the original interval, B.GT.A
37C
38C           BL     - Double precision
39C                    Lower limit of integration, BL.GE.A
40C
41C           BR     - Double precision
42C                    Upper limit of integration, BR.LE.B
43C
44C           ALFA   - Double precision
45C                    PARAMETER IN THE WEIGHT FUNCTION
46C
47C           BETA   - Double precision
48C                    Parameter in the weight function
49C
50C           RI,RJ,RG,RH - Double precision
51C                    Modified CHEBYSHEV moments for the application
52C                    of the generalized CLENSHAW-CURTIS
53C                    method (computed in subroutine DQMOMO)
54C
55C           RESULT - Double precision
56C                    Approximation to the integral
57C                    RESULT is computed by using a generalized
58C                    CLENSHAW-CURTIS method if B1 = A or BR = B.
59C                    in all other cases the 15-POINT KRONROD
60C                    RULE is applied, obtained by optimal addition of
61C                    Abscissae to the 7-POINT GAUSS RULE.
62C
63C           ABSERR - Double precision
64C                    Estimate of the modulus of the absolute error,
65C                    which should equal or exceed ABS(I-RESULT)
66C
67C           RESASC - Double precision
68C                    Approximation to the integral of ABS(F*W-I/(B-A))
69C
70C           INTEGR - Integer
71C                    Which determines the weight function
72C                    = 1   W(X) = (X-A)**ALFA*(B-X)**BETA
73C                    = 2   W(X) = (X-A)**ALFA*(B-X)**BETA*LOG(X-A)
74C                    = 3   W(X) = (X-A)**ALFA*(B-X)**BETA*LOG(B-X)
75C                    = 4   W(X) = (X-A)**ALFA*(B-X)**BETA*LOG(X-A)*
76C                                 LOG(B-X)
77C
78C           NEV    - Integer
79C                    Number of integrand evaluations
80C
81C***REFERENCES  (NONE)
82C***ROUTINES CALLED  DQCHEB, DQK15W, DQWGTS
83C***REVISION HISTORY  (YYMMDD)
84C   810101  DATE WRITTEN
85C   890531  Changed all specific intrinsics to generic.  (WRB)
86C   890531  REVISION DATE from Version 3.2
87C   891214  Prologue converted to Version 4.0 format.  (BAB)
88C***END PROLOGUE  DQC25S
89C
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
94C
95      DIMENSION CHEB12(13),CHEB24(25),FVAL(25),RG(25),RH(25),RI(25),
96     1  RJ(25),X(11)
97C
98      EXTERNAL F, DQWGTS
99C
100C           THE VECTOR X CONTAINS THE VALUES COS(K*PI/24)
101C           K = 1, ..., 11, TO BE USED FOR THE COMPUTATION OF THE
102C           CHEBYSHEV SERIES EXPANSION OF F.
103C
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/
112C
113C           LIST OF MAJOR VARIABLES
114C           -----------------------
115C
116C           FVAL   - VALUE OF THE FUNCTION F AT THE POINTS
117C                    (BR-BL)*0.5*COS(K*PI/24)+(BR+BL)*0.5
118C                    K = 0, ..., 24
119C           CHEB12 - COEFFICIENTS OF THE CHEBYSHEV SERIES EXPANSION
120C                    OF DEGREE 12, FOR THE FUNCTION F, IN THE
121C                    INTERVAL (BL,BR)
122C           CHEB24 - COEFFICIENTS OF THE CHEBYSHEV SERIES EXPANSION
123C                    OF DEGREE 24, FOR THE FUNCTION F, IN THE
124C                    INTERVAL (BL,BR)
125C           RES12  - APPROXIMATION TO THE INTEGRAL OBTAINED FROM CHEB12
126C           RES24  - APPROXIMATION TO THE INTEGRAL OBTAINED FROM CHEB24
127C           DQWGTS - EXTERNAL FUNCTION SUBPROGRAM DEFINING
128C                    THE FOUR POSSIBLE WEIGHT FUNCTIONS
129C           HLGTH  - HALF-LENGTH OF THE INTERVAL (BL,BR)
130C           CENTR  - MID POINT OF THE INTERVAL (BL,BR)
131C
132C***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
138C
139C           IF A.GT.BL AND B.LT.BR, APPLY THE 15-POINT GAUSS-KRONROD
140C           SCHEME.
141C
142C
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
147C
148C           THIS PART OF THE PROGRAM IS EXECUTED ONLY IF A = BL.
149C           ----------------------------------------------------
150C
151C           COMPUTE THE CHEBYSHEV SERIES EXPANSION OF THE
152C           FOLLOWING FUNCTION
153C           F1 = (0.5*(B+B-BR-A)-0.5*(BR-A)*X)**BETA
154C                  *F(0.5*(BR-A)*X+0.5*(BR+A))
155C
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)
175C
176C           INTEGR = 1  (OR 2)
177C
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
186C
187C           INTEGR = 2
188C
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
202C
203C           COMPUTE THE CHEBYSHEV SERIES EXPANSION OF THE
204C           FOLLOWING FUNCTION
205C           F4 = F1*LOG(0.5*(B+B-BR-A)-0.5*(BR-A)*X)
206C
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)
217C
218C           INTEGR = 3  (OR 4)
219C
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
228C
229C           INTEGR = 4
230C
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
246C
247C           THIS PART OF THE PROGRAM IS EXECUTED ONLY IF B = BR.
248C           ----------------------------------------------------
249C
250C           COMPUTE THE CHEBYSHEV SERIES EXPANSION OF THE
251C           FOLLOWING FUNCTION
252C           F2 = (0.5*(B+BL-A-A)+0.5*(B-BL)*X)**ALFA
253C                *F(0.5*(B-BL)*X+0.5*(B+BL))
254C
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
273C
274C           INTEGR = 1  (OR 3)
275C
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
285C
286C           INTEGR = 3
287C
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
301C
302C           COMPUTE THE CHEBYSHEV SERIES EXPANSION OF THE
303C           FOLLOWING FUNCTION
304C           F3 = F2*LOG(0.5*(B-BL)*X+0.5*(B+BL-A-A))
305C
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)
316C
317C           INTEGR = 2  (OR 4)
318C
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
332C
333C           INTEGR = 4
334C
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
346