1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
|
C Copyright 1981-2007 ECMWF
C
C Licensed under the GNU Lesser General Public License which
C incorporates the terms and conditions of version 3 of the GNU
C General Public License.
C See LICENSE and gpl-3.0.txt for details.
C
SUBROUTINE FFT99(A,WORK,TRIGS,IFAX,INC,JUMP,N,LOT,ISIGN)
C
C---->
C**** FFT99
C
C PURPOSE
C _______
C
C Multiple fast real periodic transform.
C
C
C INTERFACE
C _________
C
C CALL FFT99(A,WORK,TRIGS,IFAX,INC,JUMP,N,LOT,ISIGN)
C
C
C Input parameters
C ________________
C
C A - the array containing input data
C WORK - an area of size (N+1)*MIN(LOT,64)
C TRIGS - a previously prepared list of trig function values
C IFAX - a previously prepared list of factors of N
C INC - the increment within each data 'vector'
C (e.g. INC=1 for consecutively stored data)
C JUMP - the increment between the start of each data vector
C N - the length of the data vectors
C LOT - the number of data vectors
C ISIGN - +1 for transform from spectral to gridpoint
C - -1 for transform from gridpoint to spectral
C
C Output parameters
C ________________
C
C A - the array containing output data
C
C
C Method
C ______
C
C Ordering of coefficients:
C A(0),B(0),A(1),B(1),A(2),B(2),...,A(N/2),B(N/2)
C where B(0) = B(N/2)=0; (N+2) locations required
C
C Ordering of data:
C X(N-1),X(0),X(1),X(2),...,X(N-1),X(0)
C i.e. explicit cyclic continuity; (N+2) locations required
C
C Vectorization is achieved by doing the transforms in parallel
C
C N must be composed of factors 2,3 & 5 but does not have to be even
C
C
C Real transform of length N performed by removing redundant
C operations from complex transform of length N
C Definition of transforms:
C
C ISIGN = +1:
C X(J) = SUM(K=0,...,N-1)(C(K)*EXP(2*I*J*K*PI/N))
C where C(K) = A(K)+I*B(K) and C(N-K) = A(K)-I*B(K)
C
C ISIGN = -1:
C A(K) = (1/N)*SUM(J=0,...,N-1)(X(J)*COS(2*J*K*PI/N))
C B(K) = -(1/N)*SUM(J=0,...,N-1)(X(J)*SIN(2*J*K*PI/N))
C
C Externals
C _________
C
C RPASSM - Performs one pass through data as part of multiple real
C FFT (fourier analysis) routine.
C QPASSM - Performs one pass through data as part of multiple real
C FFT (fourier synthesis) routine.
C
C
C Reference
C _________
C
C None.
C
C
C Comments
C ________
C
C Tidy up of code in older version of same routine.
C
C
C AUTHOR
C ______
C
C J.D.Chambers *ECMWF* Nov 1996
C
C
C MODIFICATIONS
C _____________
C
C None.
C
C----<
C _______________________________________________________
C
C* Section 0. Definition of variables.
C _______________________________________________________
C
IMPLICIT NONE
C
C Subroutine arguments
REAL A, WORK, TRIGS
INTEGER IFAX, INC, JUMP, N, LOT, ISIGN
DIMENSION A(N),WORK(N),TRIGS(N),IFAX(10)
C
C Local variables
INTEGER NFAX, NX, NBLOX, NVEX, ISTART, NB, IA, I, J
INTEGER LA, IGO, K, IFAC, IERR, IBASE, JBASE, JJ, II
INTEGER IX, IZ
C
C ------------------------------------------------------------------
C Section 1. Initialise
C ------------------------------------------------------------------
C
100 CONTINUE
C
C Ensure factorization of N has been done.
IF(IFAX(10).NE.N) CALL SET99(TRIGS,IFAX,N)
C
NFAX = IFAX(1)
IF (MOD(N,2).EQ.1) THEN
NX = N
ELSE
NX = N + 1
ENDIF
C
C Calculate number of blocks of 64 vectors and number of vectors
C 'left over'. This remainder is transformed first.
C
NBLOX = 1 + (LOT-1)/64
NVEX = LOT-(NBLOX-1)*64
C
IF (ISIGN.EQ.1) THEN
C
C ------------------------------------------------------------------
C Section 2. Spectral to gridpoint transform.
C ------------------------------------------------------------------
C
200 CONTINUE
C
C Loop through the blocks of vectors
ISTART = 1
DO 270 NB = 1,NBLOX
IA = ISTART
I = ISTART
DO 210 J = 1,NVEX
A(I+INC) = 0.5*A(I)
I = I + JUMP
210 CONTINUE
C
IF (MOD(N,2).NE.1) THEN
I = ISTART + N*INC
DO 220 J = 1,NVEX
A(I) = 0.5*A(I)
I = I + JUMP
220 CONTINUE
ENDIF
C
IA = ISTART + INC
LA = 1
IGO = 1
C
C Work through the factors
DO 230 K = 1,NFAX
IFAC = IFAX(K+1)
IERR = -1
IF (IGO.NE.-1) THEN
CALL RPASSM(A(IA),A(IA+LA*INC),WORK(1),WORK(IFAC*LA+1),
X TRIGS,INC,1,JUMP,NX,NVEX,N,IFAC,LA,IERR)
ELSE
CALL RPASSM(WORK(1),WORK(LA+1),A(IA),A(IA+IFAC*LA*INC),
X TRIGS,1,INC,NX,JUMP,NVEX,N,IFAC,LA,IERR)
ENDIF
C
IF (IERR.NE.0) GO TO 950
C
LA = IFAC*LA
IGO = -IGO
IA = ISTART + INC
230 CONTINUE
C
C If necessary, copy results back to A
C
IF (MOD(NFAX,2).NE.0) THEN
IBASE = 1
JBASE = IA
DO 250 JJ = 1,NVEX
I = IBASE
J = JBASE
DO 240 II = 1,N
A(J) = WORK(I)
I = I + 1
J = J + INC
240 CONTINUE
IBASE = IBASE + NX
JBASE = JBASE + JUMP
250 CONTINUE
ENDIF
C
C Fill in cyclic boundary values (ie repeat the data vector
C end points at opposite end of the vector)
C
IX = ISTART
IZ = ISTART + N*INC
CDIR$ IVDEP
DO 260 J = 1,NVEX
A(IX) = A(IZ)
A(IZ+INC) = A(IX+INC)
IX = IX + JUMP
IZ = IZ + JUMP
260 CONTINUE
C
C Adjust pointers for next block
ISTART = ISTART + NVEX*JUMP
NVEX = 64
270 CONTINUE
C
ELSE
C
C ------------------------------------------------------------------
C Section 3. Gridpoint to spectral transform.
C ------------------------------------------------------------------
C
300 CONTINUE
C
C Loop through the blocks of vectors
ISTART = 1
DO 390 NB = 1,NBLOX
IA = ISTART + INC
LA = N
IGO = 1
C
DO 310 K = 1,NFAX
IFAC = IFAX(NFAX+2-K)
LA = LA/IFAC
IERR = -1
IF (IGO.NE.-1) THEN
CALL QPASSM(A(IA),A(IA+IFAC*LA*INC),WORK(1),WORK(LA+1),
X TRIGS,INC,1,JUMP,NX,NVEX,N,IFAC,LA,IERR)
ELSE
CALL QPASSM(WORK(1),WORK(IFAC*LA+1),A(IA),A(IA+LA*INC),
X TRIGS,1,INC,NX,JUMP,NVEX,N,IFAC,LA,IERR)
ENDIF
IF (IERR.NE.0) GO TO 950
IGO = -IGO
IA = ISTART + INC
310 CONTINUE
C
C If necessary, copy results back to A
C
IF (MOD(NFAX,2).NE.0) THEN
IBASE = 1
JBASE = IA
DO 330 JJ = 1,NVEX
I = IBASE
J = JBASE
DO 320 II = 1,N
A(J) = WORK(I)
I = I + 1
J = J + INC
320 CONTINUE
IBASE = IBASE + NX
JBASE = JBASE + JUMP
330 CONTINUE
ENDIF
C
C Shift A(0) and fill in zero imaginary parts
C
IX = ISTART
DO 340 J = 1,NVEX
A(IX) = A(IX+INC)
A(IX+INC) = 0.0
IX = IX + JUMP
340 CONTINUE
C
IF (MOD(N,2).NE.1) THEN
IZ = ISTART + (N+1)*INC
DO 350 J = 1,NVEX
A(IZ) = 0.0
IZ = IZ + JUMP
350 CONTINUE
C
ENDIF
C
C Adjust pointers for next block
ISTART = ISTART + NVEX*JUMP
NVEX = 64
390 CONTINUE
C
ENDIF
C
C ------------------------------------------------------------------
C Section 9. Closedown.
C ------------------------------------------------------------------
C
900 CONTINUE
C
RETURN
C
C Error messages
C
950 CONTINUE
GO TO (960,970,980) IERR
C
960 CONTINUE
WRITE(*,*) 'Vector length greater than 64, = ', NVEX
GO TO 900
970 CONTINUE
WRITE(*,*) 'Factor not handled, =', IFAC
GO TO 900
980 CONTINUE
WRITE(*,*) 'Factor only handled if LA*IFAC=N. Factor = ', IFAC
C
END
|
