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

« back to all changes in this revision

Viewing changes to Lib/fftpack/dfftpack/zfftb1.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
* 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:
1
 
      SUBROUTINE ZFFTB1 (N,C,CH,WA,IFAC)
2
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3
 
      DIMENSION       CH(*)      ,C(*)       ,WA(*)      ,IFAC(*)
4
 
      NF = IFAC(2)
5
 
      NA = 0
6
 
      L1 = 1
7
 
      IW = 1
8
 
      DO 116 K1=1,NF
9
 
         IP = IFAC(K1+2)
10
 
         L2 = IP*L1
11
 
         IDO = N/L2
12
 
         IDOT = IDO+IDO
13
 
         IDL1 = IDOT*L1
14
 
         IF (IP .NE. 4) GO TO 103
15
 
         IX2 = IW+IDOT
16
 
         IX3 = IX2+IDOT
17
 
         IF (NA .NE. 0) GO TO 101
18
 
         CALL DPASSB4 (IDOT,L1,C,CH,WA(IW),WA(IX2),WA(IX3))
19
 
         GO TO 102
20
 
  101    CALL DPASSB4 (IDOT,L1,CH,C,WA(IW),WA(IX2),WA(IX3))
21
 
  102    NA = 1-NA
22
 
         GO TO 115
23
 
  103    IF (IP .NE. 2) GO TO 106
24
 
         IF (NA .NE. 0) GO TO 104
25
 
         CALL DPASSB2 (IDOT,L1,C,CH,WA(IW))
26
 
         GO TO 105
27
 
  104    CALL DPASSB2 (IDOT,L1,CH,C,WA(IW))
28
 
  105    NA = 1-NA
29
 
         GO TO 115
30
 
  106    IF (IP .NE. 3) GO TO 109
31
 
         IX2 = IW+IDOT
32
 
         IF (NA .NE. 0) GO TO 107
33
 
         CALL DPASSB3 (IDOT,L1,C,CH,WA(IW),WA(IX2))
34
 
         GO TO 108
35
 
  107    CALL DPASSB3 (IDOT,L1,CH,C,WA(IW),WA(IX2))
36
 
  108    NA = 1-NA
37
 
         GO TO 115
38
 
  109    IF (IP .NE. 5) GO TO 112
39
 
         IX2 = IW+IDOT
40
 
         IX3 = IX2+IDOT
41
 
         IX4 = IX3+IDOT
42
 
         IF (NA .NE. 0) GO TO 110
43
 
         CALL DPASSB5 (IDOT,L1,C,CH,WA(IW),WA(IX2),WA(IX3),WA(IX4))
44
 
         GO TO 111
45
 
  110    CALL DPASSB5 (IDOT,L1,CH,C,WA(IW),WA(IX2),WA(IX3),WA(IX4))
46
 
  111    NA = 1-NA
47
 
         GO TO 115
48
 
  112    IF (NA .NE. 0) GO TO 113
49
 
         CALL DPASSB (NAC,IDOT,IP,L1,IDL1,C,C,C,CH,CH,WA(IW))
50
 
         GO TO 114
51
 
  113    CALL DPASSB (NAC,IDOT,IP,L1,IDL1,CH,CH,CH,C,C,WA(IW))
52
 
  114    IF (NAC .NE. 0) NA = 1-NA
53
 
  115    L1 = L2
54
 
         IW = IW+(IP-1)*IDOT
55
 
  116 CONTINUE
56
 
      IF (NA .EQ. 0) RETURN
57
 
      N2 = N+N
58
 
      DO 117 I=1,N2
59
 
         C(I) = CH(I)
60
 
  117 CONTINUE
61
 
      RETURN
62
 
      END
63
 
 
64
 
      SUBROUTINE DPASSB (NAC,IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA)
65
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
66
 
      DIMENSION       CH(IDO,L1,IP)          ,CC(IDO,IP,L1)          ,
67
 
     1                C1(IDO,L1,IP)          ,WA(1)      ,C2(IDL1,IP),
68
 
     2                CH2(IDL1,IP)
69
 
      IDOT = IDO/2
70
 
      NT = IP*IDL1
71
 
      IPP2 = IP+2
72
 
      IPPH = (IP+1)/2
73
 
      IDP = IP*IDO
74
 
C
75
 
      IF (IDO .LT. L1) GO TO 106
76
 
      DO 103 J=2,IPPH
77
 
         JC = IPP2-J
78
 
         DO 102 K=1,L1
79
 
            DO 101 I=1,IDO
80
 
               CH(I,K,J) = CC(I,J,K)+CC(I,JC,K)
81
 
               CH(I,K,JC) = CC(I,J,K)-CC(I,JC,K)
82
 
  101       CONTINUE
83
 
  102    CONTINUE
84
 
  103 CONTINUE
85
 
      DO 105 K=1,L1
86
 
         DO 104 I=1,IDO
87
 
            CH(I,K,1) = CC(I,1,K)
88
 
  104    CONTINUE
89
 
  105 CONTINUE
90
 
      GO TO 112
91
 
  106 DO 109 J=2,IPPH
92
 
         JC = IPP2-J
93
 
         DO 108 I=1,IDO
94
 
            DO 107 K=1,L1
95
 
               CH(I,K,J) = CC(I,J,K)+CC(I,JC,K)
96
 
               CH(I,K,JC) = CC(I,J,K)-CC(I,JC,K)
97
 
  107       CONTINUE
98
 
  108    CONTINUE
99
 
  109 CONTINUE
100
 
      DO 111 I=1,IDO
101
 
         DO 110 K=1,L1
102
 
            CH(I,K,1) = CC(I,1,K)
103
 
  110    CONTINUE
104
 
  111 CONTINUE
105
 
  112 IDL = 2-IDO
106
 
      INC = 0
107
 
      DO 116 L=2,IPPH
108
 
         LC = IPP2-L
109
 
         IDL = IDL+IDO
110
 
         DO 113 IK=1,IDL1
111
 
            C2(IK,L) = CH2(IK,1)+WA(IDL-1)*CH2(IK,2)
112
 
            C2(IK,LC) = WA(IDL)*CH2(IK,IP)
113
 
  113    CONTINUE
114
 
         IDLJ = IDL
115
 
         INC = INC+IDO
116
 
         DO 115 J=3,IPPH
117
 
            JC = IPP2-J
118
 
            IDLJ = IDLJ+INC
119
 
            IF (IDLJ .GT. IDP) IDLJ = IDLJ-IDP
120
 
            WAR = WA(IDLJ-1)
121
 
            WAI = WA(IDLJ)
122
 
            DO 114 IK=1,IDL1
123
 
               C2(IK,L) = C2(IK,L)+WAR*CH2(IK,J)
124
 
               C2(IK,LC) = C2(IK,LC)+WAI*CH2(IK,JC)
125
 
  114       CONTINUE
126
 
  115    CONTINUE
127
 
  116 CONTINUE
128
 
      DO 118 J=2,IPPH
129
 
         DO 117 IK=1,IDL1
130
 
            CH2(IK,1) = CH2(IK,1)+CH2(IK,J)
131
 
  117    CONTINUE
132
 
  118 CONTINUE
133
 
      DO 120 J=2,IPPH
134
 
         JC = IPP2-J
135
 
         DO 119 IK=2,IDL1,2
136
 
            CH2(IK-1,J) = C2(IK-1,J)-C2(IK,JC)
137
 
            CH2(IK-1,JC) = C2(IK-1,J)+C2(IK,JC)
138
 
            CH2(IK,J) = C2(IK,J)+C2(IK-1,JC)
139
 
            CH2(IK,JC) = C2(IK,J)-C2(IK-1,JC)
140
 
  119    CONTINUE
141
 
  120 CONTINUE
142
 
      NAC = 1
143
 
      IF (IDO .EQ. 2) RETURN
144
 
      NAC = 0
145
 
      DO 121 IK=1,IDL1
146
 
         C2(IK,1) = CH2(IK,1)
147
 
  121 CONTINUE
148
 
      DO 123 J=2,IP
149
 
         DO 122 K=1,L1
150
 
            C1(1,K,J) = CH(1,K,J)
151
 
            C1(2,K,J) = CH(2,K,J)
152
 
  122    CONTINUE
153
 
  123 CONTINUE
154
 
      IF (IDOT .GT. L1) GO TO 127
155
 
      IDIJ = 0
156
 
      DO 126 J=2,IP
157
 
         IDIJ = IDIJ+2
158
 
         DO 125 I=4,IDO,2
159
 
            IDIJ = IDIJ+2
160
 
            DO 124 K=1,L1
161
 
               C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)-WA(IDIJ)*CH(I,K,J)
162
 
               C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)+WA(IDIJ)*CH(I-1,K,J)
163
 
  124       CONTINUE
164
 
  125    CONTINUE
165
 
  126 CONTINUE
166
 
      RETURN
167
 
  127 IDJ = 2-IDO
168
 
      DO 130 J=2,IP
169
 
         IDJ = IDJ+IDO
170
 
         DO 129 K=1,L1
171
 
            IDIJ = IDJ
172
 
            DO 128 I=4,IDO,2
173
 
               IDIJ = IDIJ+2
174
 
               C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)-WA(IDIJ)*CH(I,K,J)
175
 
               C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)+WA(IDIJ)*CH(I-1,K,J)
176
 
  128       CONTINUE
177
 
  129    CONTINUE
178
 
  130 CONTINUE
179
 
      RETURN
180
 
      END
181
 
 
182
 
      SUBROUTINE DPASSB2 (IDO,L1,CC,CH,WA1)
183
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
184
 
      DIMENSION       CC(IDO,2,L1)           ,CH(IDO,L1,2)           ,
185
 
     1                WA1(1)
186
 
      IF (IDO .GT. 2) GO TO 102
187
 
      DO 101 K=1,L1
188
 
         CH(1,K,1) = CC(1,1,K)+CC(1,2,K)
189
 
         CH(1,K,2) = CC(1,1,K)-CC(1,2,K)
190
 
         CH(2,K,1) = CC(2,1,K)+CC(2,2,K)
191
 
         CH(2,K,2) = CC(2,1,K)-CC(2,2,K)
192
 
  101 CONTINUE
193
 
      RETURN
194
 
  102 DO 104 K=1,L1
195
 
         DO 103 I=2,IDO,2
196
 
            CH(I-1,K,1) = CC(I-1,1,K)+CC(I-1,2,K)
197
 
            TR2 = CC(I-1,1,K)-CC(I-1,2,K)
198
 
            CH(I,K,1) = CC(I,1,K)+CC(I,2,K)
199
 
            TI2 = CC(I,1,K)-CC(I,2,K)
200
 
            CH(I,K,2) = WA1(I-1)*TI2+WA1(I)*TR2
201
 
            CH(I-1,K,2) = WA1(I-1)*TR2-WA1(I)*TI2
202
 
  103    CONTINUE
203
 
  104 CONTINUE
204
 
      RETURN
205
 
      END
206
 
 
207
 
      SUBROUTINE DPASSB3 (IDO,L1,CC,CH,WA1,WA2)
208
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
209
 
      DIMENSION       CC(IDO,3,L1)           ,CH(IDO,L1,3)           ,
210
 
     1                WA1(1)     ,WA2(1)
211
 
C     *** TAUI IS SQRT(3)/2 *** 
212
 
      DATA TAUR,TAUI /-0.5D0,0.86602540378443864676D0/
213
 
      IF (IDO .NE. 2) GO TO 102
214
 
      DO 101 K=1,L1
215
 
         TR2 = CC(1,2,K)+CC(1,3,K)
216
 
         CR2 = CC(1,1,K)+TAUR*TR2
217
 
         CH(1,K,1) = CC(1,1,K)+TR2
218
 
         TI2 = CC(2,2,K)+CC(2,3,K)
219
 
         CI2 = CC(2,1,K)+TAUR*TI2
220
 
         CH(2,K,1) = CC(2,1,K)+TI2
221
 
         CR3 = TAUI*(CC(1,2,K)-CC(1,3,K))
222
 
         CI3 = TAUI*(CC(2,2,K)-CC(2,3,K))
223
 
         CH(1,K,2) = CR2-CI3
224
 
         CH(1,K,3) = CR2+CI3
225
 
         CH(2,K,2) = CI2+CR3
226
 
         CH(2,K,3) = CI2-CR3
227
 
  101 CONTINUE
228
 
      RETURN
229
 
  102 DO 104 K=1,L1
230
 
         DO 103 I=2,IDO,2
231
 
            TR2 = CC(I-1,2,K)+CC(I-1,3,K)
232
 
            CR2 = CC(I-1,1,K)+TAUR*TR2
233
 
            CH(I-1,K,1) = CC(I-1,1,K)+TR2
234
 
            TI2 = CC(I,2,K)+CC(I,3,K)
235
 
            CI2 = CC(I,1,K)+TAUR*TI2
236
 
            CH(I,K,1) = CC(I,1,K)+TI2
237
 
            CR3 = TAUI*(CC(I-1,2,K)-CC(I-1,3,K))
238
 
            CI3 = TAUI*(CC(I,2,K)-CC(I,3,K))
239
 
            DR2 = CR2-CI3
240
 
            DR3 = CR2+CI3
241
 
            DI2 = CI2+CR3
242
 
            DI3 = CI2-CR3
243
 
            CH(I,K,2) = WA1(I-1)*DI2+WA1(I)*DR2
244
 
            CH(I-1,K,2) = WA1(I-1)*DR2-WA1(I)*DI2
245
 
            CH(I,K,3) = WA2(I-1)*DI3+WA2(I)*DR3
246
 
            CH(I-1,K,3) = WA2(I-1)*DR3-WA2(I)*DI3
247
 
  103    CONTINUE
248
 
  104 CONTINUE
249
 
      RETURN
250
 
      END
251
 
 
252
 
 
253
 
      SUBROUTINE DPASSB4 (IDO,L1,CC,CH,WA1,WA2,WA3)
254
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z) 
255
 
      DIMENSION       CC(IDO,4,L1)           ,CH(IDO,L1,4)           ,
256
 
     1                WA1(1)     ,WA2(1)     ,WA3(1)
257
 
      IF (IDO .NE. 2) GO TO 102
258
 
      DO 101 K=1,L1
259
 
         TI1 = CC(2,1,K)-CC(2,3,K)
260
 
         TI2 = CC(2,1,K)+CC(2,3,K)
261
 
         TR4 = CC(2,4,K)-CC(2,2,K)
262
 
         TI3 = CC(2,2,K)+CC(2,4,K)
263
 
         TR1 = CC(1,1,K)-CC(1,3,K)
264
 
         TR2 = CC(1,1,K)+CC(1,3,K)
265
 
         TI4 = CC(1,2,K)-CC(1,4,K)
266
 
         TR3 = CC(1,2,K)+CC(1,4,K)
267
 
         CH(1,K,1) = TR2+TR3
268
 
         CH(1,K,3) = TR2-TR3
269
 
         CH(2,K,1) = TI2+TI3
270
 
         CH(2,K,3) = TI2-TI3
271
 
         CH(1,K,2) = TR1+TR4
272
 
         CH(1,K,4) = TR1-TR4
273
 
         CH(2,K,2) = TI1+TI4
274
 
         CH(2,K,4) = TI1-TI4
275
 
  101 CONTINUE
276
 
      RETURN
277
 
  102 DO 104 K=1,L1
278
 
         DO 103 I=2,IDO,2
279
 
            TI1 = CC(I,1,K)-CC(I,3,K)
280
 
            TI2 = CC(I,1,K)+CC(I,3,K)
281
 
            TI3 = CC(I,2,K)+CC(I,4,K)
282
 
            TR4 = CC(I,4,K)-CC(I,2,K)
283
 
            TR1 = CC(I-1,1,K)-CC(I-1,3,K)
284
 
            TR2 = CC(I-1,1,K)+CC(I-1,3,K)
285
 
            TI4 = CC(I-1,2,K)-CC(I-1,4,K)
286
 
            TR3 = CC(I-1,2,K)+CC(I-1,4,K)
287
 
            CH(I-1,K,1) = TR2+TR3
288
 
            CR3 = TR2-TR3
289
 
            CH(I,K,1) = TI2+TI3
290
 
            CI3 = TI2-TI3
291
 
            CR2 = TR1+TR4
292
 
            CR4 = TR1-TR4
293
 
            CI2 = TI1+TI4
294
 
            CI4 = TI1-TI4
295
 
            CH(I-1,K,2) = WA1(I-1)*CR2-WA1(I)*CI2
296
 
            CH(I,K,2) = WA1(I-1)*CI2+WA1(I)*CR2
297
 
            CH(I-1,K,3) = WA2(I-1)*CR3-WA2(I)*CI3
298
 
            CH(I,K,3) = WA2(I-1)*CI3+WA2(I)*CR3
299
 
            CH(I-1,K,4) = WA3(I-1)*CR4-WA3(I)*CI4
300
 
            CH(I,K,4) = WA3(I-1)*CI4+WA3(I)*CR4
301
 
  103    CONTINUE
302
 
  104 CONTINUE
303
 
      RETURN
304
 
      END
305
 
 
306
 
      SUBROUTINE DPASSB5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4)
307
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
308
 
      DIMENSION       CC(IDO,5,L1)           ,CH(IDO,L1,5)           ,
309
 
     1                WA1(1)     ,WA2(1)     ,WA3(1)     ,WA4(1)
310
 
C     *** TR11=COS(2*PI/5), TI11=SIN(2*PI/5)
311
 
C     *** TR12=COS(4*PI/5), TI12=SIN(4*PI/5)      
312
 
      DATA TR11,TI11,TR12,TI12 /0.3090169943749474241D0,
313
 
     +     0.95105651629515357212D0,
314
 
     +     -0.8090169943749474241D0,0.58778525229247312917D0/
315
 
      IF (IDO .NE. 2) GO TO 102
316
 
      DO 101 K=1,L1
317
 
         TI5 = CC(2,2,K)-CC(2,5,K)
318
 
         TI2 = CC(2,2,K)+CC(2,5,K)
319
 
         TI4 = CC(2,3,K)-CC(2,4,K)
320
 
         TI3 = CC(2,3,K)+CC(2,4,K)
321
 
         TR5 = CC(1,2,K)-CC(1,5,K)
322
 
         TR2 = CC(1,2,K)+CC(1,5,K)
323
 
         TR4 = CC(1,3,K)-CC(1,4,K)
324
 
         TR3 = CC(1,3,K)+CC(1,4,K)
325
 
         CH(1,K,1) = CC(1,1,K)+TR2+TR3
326
 
         CH(2,K,1) = CC(2,1,K)+TI2+TI3
327
 
         CR2 = CC(1,1,K)+TR11*TR2+TR12*TR3
328
 
         CI2 = CC(2,1,K)+TR11*TI2+TR12*TI3
329
 
         CR3 = CC(1,1,K)+TR12*TR2+TR11*TR3
330
 
         CI3 = CC(2,1,K)+TR12*TI2+TR11*TI3
331
 
         CR5 = TI11*TR5+TI12*TR4
332
 
         CI5 = TI11*TI5+TI12*TI4
333
 
         CR4 = TI12*TR5-TI11*TR4
334
 
         CI4 = TI12*TI5-TI11*TI4
335
 
         CH(1,K,2) = CR2-CI5
336
 
         CH(1,K,5) = CR2+CI5
337
 
         CH(2,K,2) = CI2+CR5
338
 
         CH(2,K,3) = CI3+CR4
339
 
         CH(1,K,3) = CR3-CI4
340
 
         CH(1,K,4) = CR3+CI4
341
 
         CH(2,K,4) = CI3-CR4
342
 
         CH(2,K,5) = CI2-CR5
343
 
  101 CONTINUE
344
 
      RETURN
345
 
  102 DO 104 K=1,L1
346
 
         DO 103 I=2,IDO,2
347
 
            TI5 = CC(I,2,K)-CC(I,5,K)
348
 
            TI2 = CC(I,2,K)+CC(I,5,K)
349
 
            TI4 = CC(I,3,K)-CC(I,4,K)
350
 
            TI3 = CC(I,3,K)+CC(I,4,K)
351
 
            TR5 = CC(I-1,2,K)-CC(I-1,5,K)
352
 
            TR2 = CC(I-1,2,K)+CC(I-1,5,K)
353
 
            TR4 = CC(I-1,3,K)-CC(I-1,4,K)
354
 
            TR3 = CC(I-1,3,K)+CC(I-1,4,K)
355
 
            CH(I-1,K,1) = CC(I-1,1,K)+TR2+TR3
356
 
            CH(I,K,1) = CC(I,1,K)+TI2+TI3
357
 
            CR2 = CC(I-1,1,K)+TR11*TR2+TR12*TR3
358
 
            CI2 = CC(I,1,K)+TR11*TI2+TR12*TI3
359
 
            CR3 = CC(I-1,1,K)+TR12*TR2+TR11*TR3
360
 
            CI3 = CC(I,1,K)+TR12*TI2+TR11*TI3
361
 
            CR5 = TI11*TR5+TI12*TR4
362
 
            CI5 = TI11*TI5+TI12*TI4
363
 
            CR4 = TI12*TR5-TI11*TR4
364
 
            CI4 = TI12*TI5-TI11*TI4
365
 
            DR3 = CR3-CI4
366
 
            DR4 = CR3+CI4
367
 
            DI3 = CI3+CR4
368
 
            DI4 = CI3-CR4
369
 
            DR5 = CR2+CI5
370
 
            DR2 = CR2-CI5
371
 
            DI5 = CI2-CR5
372
 
            DI2 = CI2+CR5
373
 
            CH(I-1,K,2) = WA1(I-1)*DR2-WA1(I)*DI2
374
 
            CH(I,K,2) = WA1(I-1)*DI2+WA1(I)*DR2
375
 
            CH(I-1,K,3) = WA2(I-1)*DR3-WA2(I)*DI3
376
 
            CH(I,K,3) = WA2(I-1)*DI3+WA2(I)*DR3
377
 
            CH(I-1,K,4) = WA3(I-1)*DR4-WA3(I)*DI4
378
 
            CH(I,K,4) = WA3(I-1)*DI4+WA3(I)*DR4
379
 
            CH(I-1,K,5) = WA4(I-1)*DR5-WA4(I)*DI5
380
 
            CH(I,K,5) = WA4(I-1)*DI5+WA4(I)*DR5
381
 
  103    CONTINUE
382
 
  104 CONTINUE
383
 
      RETURN
384
 
      END