← Back to branch summary

~ubuntu-branches/ubuntu/wily/julia/wily

« back to all changes in this revision

Viewing changes to deps/openlibm/slatec/qc25f.f

  • Committer: Package Import Robot
  • Author(s): Sébastien Villemot
  • Date: 2013-01-16 12:29:42 UTC
  • mfrom: (1.1.2)
  • Revision ID: package-import@ubuntu.com-20130116122942-x86e42akjq31repw
Tags: 0.0.0+20130107.gitd9656f41-1
* New upstream snashot
* No longer try to rebuild helpdb.jl.
   + debian/rules: remove helpdb.jl from build-arch rule
   + debian/control: move back python-sphinx to Build-Depends-Indep
* debian/copyright: reflect upstream changes
* Add Build-Conflicts on libatlas3-base (makes linalg tests fail)
* debian/rules: replace obsolete USE_DEBIAN makeflag by a list of
  USE_SYSTEM_* flags
* debian/rules: on non-x86 systems, use libm instead of openlibm
* dpkg-buildflags.patch: remove patch, applied upstream
* Refreshed other patches

Show diffs side-by-side

added added

removed removed

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