← Back to branch summary

~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

Viewing changes to scipy/integrate/quadpack/dqc25f.f

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

Lines of Context:
scipy/integrate/quadpack/dqc25f.f
 
1
      subroutine dqc25f(f,a,b,omega,integr,nrmom,maxp1,ksave,result,
 
2
     *   abserr,neval,resabs,resasc,momcom,chebmo)
 
3
c***begin prologue  dqc25f
 
4
c***date written   810101   (yymmdd)
 
5
c***revision date  830518   (yymmdd)
 
6
c***category no.  h2a2a2
 
7
c***keywords  integration rules for functions with cos or sin
 
8
c             factor, clenshaw-curtis, gauss-kronrod
 
9
c***author  piessens,robert,appl. math. & progr. div. - k.u.leuven
 
10
c           de doncker,elise,appl. math. & progr. div. - k.u.leuven
 
11
c***purpose  to compute the integral i=integral of f(x) over (a,b)
 
12
c            where w(x) = cos(omega*x) or w(x)=sin(omega*x) and to
 
13
c            compute j = integral of abs(f) over (a,b). for small value
 
14
c            of omega or small intervals (a,b) the 15-point gauss-kronro
 
15
c            rule is used. otherwise a generalized clenshaw-curtis
 
16
c            method is used.
 
17
c***description
 
18
c
 
19
c        integration rules for functions with cos or sin factor
 
20
c        standard fortran subroutine
 
21
c        double precision version
 
22
c
 
23
c        parameters
 
24
c         on entry
 
25
c           f      - double precision
 
26
c                    function subprogram defining the integrand
 
27
c                    function f(x). the actual name for f needs to
 
28
c                    be declared e x t e r n a l in the calling program.
 
29
c
 
30
c           a      - double precision
 
31
c                    lower limit of integration
 
32
c
 
33
c           b      - double precision
 
34
c                    upper limit of integration
 
35
c
 
36
c           omega  - double precision
 
37
c                    parameter in the weight function
 
38
c
 
39
c           integr - integer
 
40
c                    indicates which weight function is to be used
 
41
c                       integr = 1   w(x) = cos(omega*x)
 
42
c                       integr = 2   w(x) = sin(omega*x)
 
43
c
 
44
c           nrmom  - integer
 
45
c                    the length of interval (a,b) is equal to the length
 
46
c                    of the original integration interval divided by
 
47
c                    2**nrmom (we suppose that the routine is used in an
 
48
c                    adaptive integration process, otherwise set
 
49
c                    nrmom = 0). nrmom must be zero at the first call.
 
50
c
 
51
c           maxp1  - integer
 
52
c                    gives an upper bound on the number of chebyshev
 
53
c                    moments which can be stored, i.e. for the
 
54
c                    intervals of lengths abs(bb-aa)*2**(-l),
 
55
c                    l = 0,1,2, ..., maxp1-2.
 
56
c
 
57
c           ksave  - integer
 
58
c                    key which is one when the moments for the
 
59
c                    current interval have been computed
 
60
c
 
61
c         on return
 
62
c           result - double precision
 
63
c                    approximation to the integral i
 
64
c
 
65
c           abserr - double precision
 
66
c                    estimate of the modulus of the absolute
 
67
c                    error, which should equal or exceed abs(i-result)
 
68
c
 
69
c           neval  - integer
 
70
c                    number of integrand evaluations
 
71
c
 
72
c           resabs - double precision
 
73
c                    approximation to the integral j
 
74
c
 
75
c           resasc - double precision
 
76
c                    approximation to the integral of abs(f-i/(b-a))
 
77
c
 
78
c         on entry and return
 
79
c           momcom - integer
 
80
c                    for each interval length we need to compute the
 
81
c                    chebyshev moments. momcom counts the number of
 
82
c                    intervals for which these moments have already been
 
83
c                    computed. if nrmom.lt.momcom or ksave = 1, the
 
84
c                    chebyshev moments for the interval (a,b) have
 
85
c                    already been computed and stored, otherwise we
 
86
c                    compute them and we increase momcom.
 
87
c
 
88
c           chebmo - double precision
 
89
c                    array of dimension at least (maxp1,25) containing
 
90
c                    the modified chebyshev moments for the first momcom
 
91
c                    momcom interval lengths
 
92
c
 
93
c ......................................................................
 
94
c***references  (none)
 
95
c***routines called  d1mach,dgtsl,dqcheb,dqk15w,dqwgtf
 
96
c***end prologue  dqc25f
 
97
c
 
98
      double precision a,abserr,ac,an,an2,as,asap,ass,b,centr,chebmo,
 
99
     *  cheb12,cheb24,conc,cons,cospar,d,dabs,dcos,dsin,dqwgtf,d1,
 
100
     *  d1mach,d2,estc,ests,f,fval,hlgth,oflow,omega,parint,par2,par22,
 
101
     *  p2,p3,p4,resabs,resasc,resc12,resc24,ress12,ress24,result,
 
102
     *  sinpar,v,x
 
103
      integer i,iers,integr,isym,j,k,ksave,m,momcom,neval,maxp1,
 
104
     *  noequ,noeq1,nrmom
 
105
c
 
106
      dimension chebmo(maxp1,25),cheb12(13),cheb24(25),d(25),d1(25),
 
107
     *  d2(25),fval(25),v(28),x(11)
 
108
c
 
109
      external f,dqwgtf
 
110
c
 
111
c           the vector x contains the values cos(k*pi/24)
 
112
c           k = 1, ...,11, to be used for the chebyshev expansion of f
 
113
c
 
114
      data x(1) / 0.9914448613 7381041114 4557526928 563d0 /
 
115
      data x(2) / 0.9659258262 8906828674 9743199728 897d0 /
 
116
      data x(3) / 0.9238795325 1128675612 8183189396 788d0 /
 
117
      data x(4) / 0.8660254037 8443864676 3723170752 936d0 /
 
118
      data x(5) / 0.7933533402 9123516457 9776961501 299d0 /
 
119
      data x(6) / 0.7071067811 8654752440 0844362104 849d0 /
 
120
      data x(7) / 0.6087614290 0872063941 6097542898 164d0 /
 
121
      data x(8) / 0.5000000000 0000000000 0000000000 000d0 /
 
122
      data x(9) / 0.3826834323 6508977172 8459984030 399d0 /
 
123
      data x(10) / 0.2588190451 0252076234 8898837624 048d0 /
 
124
      data x(11) / 0.1305261922 2005159154 8406227895 489d0 /
 
125
c
 
126
c           list of major variables
 
127
c           -----------------------
 
128
c
 
129
c           centr  - mid point of the integration interval
 
130
c           hlgth  - half-length of the integration interval
 
131
c           fval   - value of the function f at the points
 
132
c                    (b-a)*0.5*cos(k*pi/12) + (b+a)*0.5, k = 0, ..., 24
 
133
c           cheb12 - coefficients of the chebyshev series expansion
 
134
c                    of degree 12, for the function f, in the
 
135
c                    interval (a,b)
 
136
c           cheb24 - coefficients of the chebyshev series expansion
 
137
c                    of degree 24, for the function f, in the
 
138
c                    interval (a,b)
 
139
c           resc12 - approximation to the integral of
 
140
c                    cos(0.5*(b-a)*omega*x)*f(0.5*(b-a)*x+0.5*(b+a))
 
141
c                    over (-1,+1), using the chebyshev series
 
142
c                    expansion of degree 12
 
143
c           resc24 - approximation to the same integral, using the
 
144
c                    chebyshev series expansion of degree 24
 
145
c           ress12 - the analogue of resc12 for the sine
 
146
c           ress24 - the analogue of resc24 for the sine
 
147
c
 
148
c
 
149
c           machine dependent constant
 
150
c           --------------------------
 
151
c
 
152
c           oflow is the largest positive magnitude.
 
153
c
 
154
c***first executable statement  dqc25f
 
155
      oflow = d1mach(2)
 
156
c
 
157
      centr = 0.5d+00*(b+a)
 
158
      hlgth = 0.5d+00*(b-a)
 
159
      parint = omega*hlgth
 
160
c
 
161
c           compute the integral using the 15-point gauss-kronrod
 
162
c           formula if the value of the parameter in the integrand
 
163
c           is small.
 
164
c
 
165
      if(dabs(parint).gt.0.2d+01) go to 10
 
166
      call dqk15w(f,dqwgtf,omega,p2,p3,p4,integr,a,b,result,
 
167
     *  abserr,resabs,resasc)
 
168
      neval = 15
 
169
      go to 170
 
170
c
 
171
c           compute the integral using the generalized clenshaw-
 
172
c           curtis method.
 
173
c
 
174
   10 conc = hlgth*dcos(centr*omega)
 
175
      cons = hlgth*dsin(centr*omega)
 
176
      resasc = oflow
 
177
      neval = 25
 
178
c
 
179
c           check whether the chebyshev moments for this interval
 
180
c           have already been computed.
 
181
c
 
182
      if(nrmom.lt.momcom.or.ksave.eq.1) go to 120
 
183
c
 
184
c           compute a new set of chebyshev moments.
 
185
c
 
186
      m = momcom+1
 
187
      par2 = parint*parint
 
188
      par22 = par2+0.2d+01
 
189
      sinpar = dsin(parint)
 
190
      cospar = dcos(parint)
 
191
c
 
192
c           compute the chebyshev moments with respect to cosine.
 
193
c
 
194
      v(1) = 0.2d+01*sinpar/parint
 
195
      v(2) = (0.8d+01*cospar+(par2+par2-0.8d+01)*sinpar/parint)/par2
 
196
      v(3) = (0.32d+02*(par2-0.12d+02)*cospar+(0.2d+01*
 
197
     *  ((par2-0.80d+02)*par2+0.192d+03)*sinpar)/parint)/(par2*par2)
 
198
      ac = 0.8d+01*cospar
 
199
      as = 0.24d+02*parint*sinpar
 
200
      if(dabs(parint).gt.0.24d+02) go to 30
 
201
c
 
202
c           compute the chebyshev moments as the solutions of a
 
203
c           boundary value problem with 1 initial value (v(3)) and 1
 
204
c           end value (computed using an asymptotic formula).
 
205
c
 
206
      noequ = 25
 
207
      noeq1 = noequ-1
 
208
      an = 0.6d+01
 
209
      do 20 k = 1,noeq1
 
210
        an2 = an*an
 
211
        d(k) = -0.2d+01*(an2-0.4d+01)*(par22-an2-an2)
 
212
        d2(k) = (an-0.1d+01)*(an-0.2d+01)*par2
 
213
        d1(k+1) = (an+0.3d+01)*(an+0.4d+01)*par2
 
214
        v(k+3) = as-(an2-0.4d+01)*ac
 
215
        an = an+0.2d+01
 
216
   20 continue
 
217
      an2 = an*an
 
218
      d(noequ) = -0.2d+01*(an2-0.4d+01)*(par22-an2-an2)
 
219
      v(noequ+3) = as-(an2-0.4d+01)*ac
 
220
      v(4) = v(4)-0.56d+02*par2*v(3)
 
221
      ass = parint*sinpar
 
222
      asap = (((((0.210d+03*par2-0.1d+01)*cospar-(0.105d+03*par2
 
223
     *  -0.63d+02)*ass)/an2-(0.1d+01-0.15d+02*par2)*cospar
 
224
     *  +0.15d+02*ass)/an2-cospar+0.3d+01*ass)/an2-cospar)/an2
 
225
      v(noequ+3) = v(noequ+3)-0.2d+01*asap*par2*(an-0.1d+01)*
 
226
     *   (an-0.2d+01)
 
227
c
 
228
c           solve the tridiagonal system by means of gaussian
 
229
c           elimination with partial pivoting.
 
230
c
 
231
c***        call to dgtsl must be replaced by call to
 
232
c***        double precision version of linpack routine sgtsl
 
233
c
 
234
      call dgtsl(noequ,d1,d,d2,v(4),iers)
 
235
      go to 50
 
236
c
 
237
c           compute the chebyshev moments by means of forward
 
238
c           recursion.
 
239
c
 
240
   30 an = 0.4d+01
 
241
      do 40 i = 4,13
 
242
        an2 = an*an
 
243
        v(i) = ((an2-0.4d+01)*(0.2d+01*(par22-an2-an2)*v(i-1)-ac)
 
244
     *  +as-par2*(an+0.1d+01)*(an+0.2d+01)*v(i-2))/
 
245
     *  (par2*(an-0.1d+01)*(an-0.2d+01))
 
246
        an = an+0.2d+01
 
247
   40 continue
 
248
   50 do 60 j = 1,13
 
249
        chebmo(m,2*j-1) = v(j)
 
250
   60 continue
 
251
c
 
252
c           compute the chebyshev moments with respect to sine.
 
253
c
 
254
      v(1) = 0.2d+01*(sinpar-parint*cospar)/par2
 
255
      v(2) = (0.18d+02-0.48d+02/par2)*sinpar/par2
 
256
     *  +(-0.2d+01+0.48d+02/par2)*cospar/parint
 
257
      ac = -0.24d+02*parint*cospar
 
258
      as = -0.8d+01*sinpar
 
259
      if(dabs(parint).gt.0.24d+02) go to 80
 
260
c
 
261
c           compute the chebyshev moments as the solutions of a boundary
 
262
c           value problem with 1 initial value (v(2)) and 1 end value
 
263
c           (computed using an asymptotic formula).
 
264
c
 
265
      an = 0.5d+01
 
266
      do 70 k = 1,noeq1
 
267
        an2 = an*an
 
268
        d(k) = -0.2d+01*(an2-0.4d+01)*(par22-an2-an2)
 
269
        d2(k) = (an-0.1d+01)*(an-0.2d+01)*par2
 
270
        d1(k+1) = (an+0.3d+01)*(an+0.4d+01)*par2
 
271
        v(k+2) = ac+(an2-0.4d+01)*as
 
272
        an = an+0.2d+01
 
273
   70 continue
 
274
      an2 = an*an
 
275
      d(noequ) = -0.2d+01*(an2-0.4d+01)*(par22-an2-an2)
 
276
      v(noequ+2) = ac+(an2-0.4d+01)*as
 
277
      v(3) = v(3)-0.42d+02*par2*v(2)
 
278
      ass = parint*cospar
 
279
      asap = (((((0.105d+03*par2-0.63d+02)*ass+(0.210d+03*par2
 
280
     *  -0.1d+01)*sinpar)/an2+(0.15d+02*par2-0.1d+01)*sinpar-
 
281
     *  0.15d+02*ass)/an2-0.3d+01*ass-sinpar)/an2-sinpar)/an2
 
282
      v(noequ+2) = v(noequ+2)-0.2d+01*asap*par2*(an-0.1d+01)
 
283
     *  *(an-0.2d+01)
 
284
c
 
285
c           solve the tridiagonal system by means of gaussian
 
286
c           elimination with partial pivoting.
 
287
c
 
288
c***        call to dgtsl must be replaced by call to
 
289
c***        double precision version of linpack routine sgtsl
 
290
c
 
291
      call dgtsl(noequ,d1,d,d2,v(3),iers)
 
292
      go to 100
 
293
c
 
294
c           compute the chebyshev moments by means of forward recursion.
 
295
c
 
296
   80 an = 0.3d+01
 
297
      do 90 i = 3,12
 
298
        an2 = an*an
 
299
        v(i) = ((an2-0.4d+01)*(0.2d+01*(par22-an2-an2)*v(i-1)+as)
 
300
     *  +ac-par2*(an+0.1d+01)*(an+0.2d+01)*v(i-2))
 
301
     *  /(par2*(an-0.1d+01)*(an-0.2d+01))
 
302
        an = an+0.2d+01
 
303
   90 continue
 
304
  100 do 110 j = 1,12
 
305
        chebmo(m,2*j) = v(j)
 
306
  110 continue
 
307
  120 if (nrmom.lt.momcom) m = nrmom+1
 
308
       if (momcom.lt.(maxp1-1).and.nrmom.ge.momcom) momcom = momcom+1
 
309
c
 
310
c           compute the coefficients of the chebyshev expansions
 
311
c           of degrees 12 and 24 of the function f.
 
312
c
 
313
      fval(1) = 0.5d+00*f(centr+hlgth)
 
314
      fval(13) = f(centr)
 
315
      fval(25) = 0.5d+00*f(centr-hlgth)
 
316
      do 130 i = 2,12
 
317
        isym = 26-i
 
318
        fval(i) = f(hlgth*x(i-1)+centr)
 
319
        fval(isym) = f(centr-hlgth*x(i-1))
 
320
  130 continue
 
321
      call dqcheb(x,fval,cheb12,cheb24)
 
322
c
 
323
c           compute the integral and error estimates.
 
324
c
 
325
      resc12 = cheb12(13)*chebmo(m,13)
 
326
      ress12 = 0.0d+00
 
327
      k = 11
 
328
      do 140 j = 1,6
 
329
        resc12 = resc12+cheb12(k)*chebmo(m,k)
 
330
        ress12 = ress12+cheb12(k+1)*chebmo(m,k+1)
 
331
        k = k-2
 
332
  140 continue
 
333
      resc24 = cheb24(25)*chebmo(m,25)
 
334
      ress24 = 0.0d+00
 
335
      resabs = dabs(cheb24(25))
 
336
      k = 23
 
337
      do 150 j = 1,12
 
338
        resc24 = resc24+cheb24(k)*chebmo(m,k)
 
339
        ress24 = ress24+cheb24(k+1)*chebmo(m,k+1)
 
340
        resabs = dabs(cheb24(k))+dabs(cheb24(k+1))
 
341
        k = k-2
 
342
  150 continue
 
343
      estc = dabs(resc24-resc12)
 
344
      ests = dabs(ress24-ress12)
 
345
      resabs = resabs*dabs(hlgth)
 
346
      if(integr.eq.2) go to 160
 
347
      result = conc*resc24-cons*ress24
 
348
      abserr = dabs(conc*estc)+dabs(cons*ests)
 
349
      go to 170
 
350
  160 result = conc*ress24+cons*resc24
 
351
      abserr = dabs(conc*ests)+dabs(cons*estc)
 
352
  170 return
 
353
      end