(git:374b731)
Loading...
Searching...
No Matches
ps_wavelet_fft3d.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9MODULE ps_wavelet_fft3d
10
11 USE kinds, ONLY: dp
12#include "../base/base_uses.f90"
13
14 IMPLICIT NONE
15 PRIVATE
16
17 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ps_wavelet_fft3d'
18
19 ! longest fft supported, must be equal to the length of the ctrig array
20 INTEGER, PARAMETER :: ctrig_length = 8192
21
22 PUBLIC :: fourier_dim, &
23 ctrig, &
24 fftstp, ctrig_length
25
26CONTAINS
27
28! **************************************************************************************************
29!> \brief Give a number n_next > n compatible for the FFT
30!> \param n ...
31!> \param n_next ...
32! **************************************************************************************************
33 SUBROUTINE fourier_dim(n, n_next)
34 INTEGER, INTENT(in) :: n
35 INTEGER, INTENT(out) :: n_next
36
37 INTEGER, PARAMETER :: ndata = 149, ndata1024 = 149
38 INTEGER, DIMENSION(ndata), PARAMETER :: idata = (/3, 4, 5, 6, 8, 9, 12, 15, 16, 18, 20, 24, &
39 25, 27, 30, 32, 36, 40, 45, 48, 54, 60, 64, 72, 75, 80, 81, 90, 96, 100, 108, 120, 125, &
40 128, 135, 144, 150, 160, 162, 180, 192, 200, 216, 225, 240, 243, 256, 270, 288, 300, 320, &
41 324, 360, 375, 384, 400, 405, 432, 450, 480, 486, 500, 512, 540, 576, 600, 625, 640, 648, &
42 675, 720, 729, 750, 768, 800, 810, 864, 900, 960, 972, 1000, 1024, 1080, 1125, 1152, 1200,&
43 1215, 1280, 1296, 1350, 1440, 1458, 1500, 1536, 1600, 1620, 1728, 1800, 1875, 1920, 1944, &
44 2000, 2025, 2048, 2160, 2250, 2304, 2400, 2430, 2500, 2560, 2592, 2700, 2880, 3000, 3072, &
45 3125, 3200, 3240, 3375, 3456, 3600, 3750, 3840, 3888, 4000, 4050, 4096, 4320, 4500, 4608, &
46 4800, 5000, 5120, 5184, 5400, 5625, 5760, 6000, 6144, 6400, 6480, 6750, 6912, 7200, 7500, &
47 7680, 8000, ctrig_length/)
48
49 INTEGER :: i
50
51!Multiple of 2,3,5
52
53 loop_data: DO i = 1, ndata1024
54 IF (n <= idata(i)) THEN
55 n_next = idata(i)
56 RETURN
57 END IF
58 END DO loop_data
59 WRITE (unit=*, fmt=*) "fourier_dim: ", n, " is bigger than ", idata(ndata1024)
60 cpabort("")
61 END SUBROUTINE fourier_dim
62
63! Copyright (C) Stefan Goedecker, CEA Grenoble, 2002
64! This file is distributed under the terms of the
65! GNU General Public License, see http://www.gnu.org/copyleft/gpl.txt .
66
67! --------------------------------------------------------------
68! 3-dimensional complex-complex FFT routine:
69! When compared to the best vendor implementations on RISC architectures
70! it gives close to optimal performance (perhaps losing 20 percent in speed)
71! and it is significantly faster than many not so good vendor implementations
72! as well as other portable FFT's.
73! On all vector machines tested so far (Cray, NEC, Fujitsu) is
74! was significantly faster than the vendor routines
75! The theoretical background is described in :
76! 1) S. Goedecker: Rotating a three-dimensional array in optimal
77! positions for vector processing: Case study for a three-dimensional Fast
78! Fourier Transform, Comp. Phys. Commun. \underline{76}, 294 (1993)
79! Citing of this reference is greatly appreciated if the routines are used
80! for scientific work.
81
82! Presumably good compiler flags:
83! IBM, serial power 2: xlf -qarch=pwr2 -O2 -qmaxmem=-1
84! with OpenMP: IBM: xlf_r -qfree -O4 -qarch=pwr3 -qtune=pwr3 -qsmp=omp -qmaxmem=-1 ;
85! a.out
86! DEC: f90 -O3 -arch ev67 -pipeline
87! with OpenMP: DEC: f90 -O3 -arch ev67 -pipeline -omp -lelan ;
88! prun -N1 -c4 a.out
89
90!-----------------------------------------------------------
91
92! FFT PART -----------------------------------------------------------------
93
94! **************************************************************************************************
95!> \brief ...
96!> \param n ...
97!> \param trig ...
98!> \param after ...
99!> \param before ...
100!> \param now ...
101!> \param isign ...
102!> \param ic ...
103! **************************************************************************************************
104 SUBROUTINE ctrig(n, trig, after, before, now, isign, ic)
105! Copyright (C) Stefan Goedecker, Lausanne, Switzerland, August 1, 1991
106! Copyright (C) Stefan Goedecker, Cornell University, Ithaca, USA, 1994
107! Copyright (C) Stefan Goedecker, MPI Stuttgart, Germany, 1999
108! This file is distributed under the terms of the
109! GNU General Public License, see http://www.gnu.org/copyleft/gpl.txt .
110
111! Different factorizations affect the performance
112! Factoring 64 as 4*4*4 might for example be faster on some machines than 8*8.
113 INTEGER :: n
114 REAL(kind=dp) :: trig
115 INTEGER :: after, before, now, isign, ic
116
117 INTEGER :: i, itt, j, nh
118 REAL(kind=dp) :: angle, trigc, trigs, twopi
119
120 dimension now(7), after(7), before(7), trig(2, ctrig_length)
121 INTEGER, DIMENSION(7, 149) :: idata
122! The factor 6 is only allowed in the first place!
123 DATA((idata(i, j), i=1, 7), j=1, 76)/ &
124 3, 3, 1, 1, 1, 1, 1, 4, 4, 1, 1, 1, 1, 1, &
125 5, 5, 1, 1, 1, 1, 1, 6, 6, 1, 1, 1, 1, 1, &
126 8, 8, 1, 1, 1, 1, 1, 9, 3, 3, 1, 1, 1, 1, &
127 12, 4, 3, 1, 1, 1, 1, 15, 5, 3, 1, 1, 1, 1, &
128 16, 4, 4, 1, 1, 1, 1, 18, 6, 3, 1, 1, 1, 1, &
129 20, 5, 4, 1, 1, 1, 1, 24, 8, 3, 1, 1, 1, 1, &
130 25, 5, 5, 1, 1, 1, 1, 27, 3, 3, 3, 1, 1, 1, &
131 30, 6, 5, 1, 1, 1, 1, 32, 8, 4, 1, 1, 1, 1, &
132 36, 4, 3, 3, 1, 1, 1, 40, 8, 5, 1, 1, 1, 1, &
133 45, 5, 3, 3, 1, 1, 1, 48, 4, 4, 3, 1, 1, 1, &
134 54, 6, 3, 3, 1, 1, 1, 60, 5, 4, 3, 1, 1, 1, &
135 64, 8, 8, 1, 1, 1, 1, 72, 8, 3, 3, 1, 1, 1, &
136 75, 5, 5, 3, 1, 1, 1, 80, 5, 4, 4, 1, 1, 1, &
137 81, 3, 3, 3, 3, 1, 1, 90, 6, 5, 3, 1, 1, 1, &
138 96, 8, 4, 3, 1, 1, 1, 100, 5, 5, 4, 1, 1, 1, &
139 108, 4, 3, 3, 3, 1, 1, 120, 8, 5, 3, 1, 1, 1, &
140 125, 5, 5, 5, 1, 1, 1, 128, 8, 4, 4, 1, 1, 1, &
141 135, 5, 3, 3, 3, 1, 1, 144, 6, 8, 3, 1, 1, 1, &
142 150, 6, 5, 5, 1, 1, 1, 160, 8, 5, 4, 1, 1, 1, &
143 162, 6, 3, 3, 3, 1, 1, 180, 5, 4, 3, 3, 1, 1, &
144 192, 6, 8, 4, 1, 1, 1, 200, 8, 5, 5, 1, 1, 1, &
145 216, 8, 3, 3, 3, 1, 1, 225, 5, 5, 3, 3, 1, 1, &
146 240, 6, 8, 5, 1, 1, 1, 243, 3, 3, 3, 3, 3, 1, &
147 256, 8, 8, 4, 1, 1, 1, 270, 6, 5, 3, 3, 1, 1, &
148 288, 8, 4, 3, 3, 1, 1, 300, 5, 5, 4, 3, 1, 1, &
149 320, 5, 4, 4, 4, 1, 1, 324, 4, 3, 3, 3, 3, 1, &
150 360, 8, 5, 3, 3, 1, 1, 375, 5, 5, 5, 3, 1, 1, &
151 384, 8, 4, 4, 3, 1, 1, 400, 5, 5, 4, 4, 1, 1, &
152 405, 5, 3, 3, 3, 3, 1, 432, 4, 4, 3, 3, 3, 1, &
153 450, 6, 5, 5, 3, 1, 1, 480, 8, 5, 4, 3, 1, 1, &
154 486, 6, 3, 3, 3, 3, 1, 500, 5, 5, 5, 4, 1, 1, &
155 512, 8, 8, 8, 1, 1, 1, 540, 5, 4, 3, 3, 3, 1, &
156 576, 4, 4, 4, 3, 3, 1, 600, 8, 5, 5, 3, 1, 1, &
157 625, 5, 5, 5, 5, 1, 1, 640, 8, 5, 4, 4, 1, 1, &
158 648, 8, 3, 3, 3, 3, 1, 675, 5, 5, 3, 3, 3, 1, &
159 720, 5, 4, 4, 3, 3, 1, 729, 3, 3, 3, 3, 3, 3, &
160 750, 6, 5, 5, 5, 1, 1, 768, 4, 4, 4, 4, 3, 1, &
161 800, 8, 5, 5, 4, 1, 1, 810, 6, 5, 3, 3, 3, 1/
162 DATA((idata(i, j), i=1, 7), j=77, 149)/ &
163 864, 8, 4, 3, 3, 3, 1, 900, 5, 5, 4, 3, 3, 1, &
164 960, 5, 4, 4, 4, 3, 1, 972, 4, 3, 3, 3, 3, 3, &
165 1000, 8, 5, 5, 5, 1, 1, 1024, 4, 4, 4, 4, 4, 1, &
166 1080, 6, 5, 4, 3, 3, 1, 1125, 5, 5, 5, 3, 3, 1, &
167 1152, 6, 4, 4, 4, 3, 1, 1200, 6, 8, 5, 5, 1, 1, &
168 1215, 5, 3, 3, 3, 3, 3, 1280, 8, 8, 5, 4, 1, 1, &
169 1296, 6, 8, 3, 3, 3, 1, 1350, 6, 5, 5, 3, 3, 1, &
170 1440, 6, 5, 4, 4, 3, 1, 1458, 6, 3, 3, 3, 3, 3, &
171 1500, 5, 5, 5, 4, 3, 1, 1536, 6, 8, 8, 4, 1, 1, &
172 1600, 8, 8, 5, 5, 1, 1, 1620, 5, 4, 3, 3, 3, 3, &
173 1728, 6, 8, 4, 3, 3, 1, 1800, 6, 5, 5, 4, 3, 1, &
174 1875, 5, 5, 5, 5, 3, 1, 1920, 6, 5, 4, 4, 4, 1, &
175 1944, 6, 4, 3, 3, 3, 3, 2000, 5, 5, 5, 4, 4, 1, &
176 2025, 5, 5, 3, 3, 3, 3, 2048, 8, 4, 4, 4, 4, 1, &
177 2160, 6, 8, 5, 3, 3, 1, 2250, 6, 5, 5, 5, 3, 1, &
178 2304, 6, 8, 4, 4, 3, 1, 2400, 6, 5, 5, 4, 4, 1, &
179 2430, 6, 5, 3, 3, 3, 3, 2500, 5, 5, 5, 5, 4, 1, &
180 2560, 8, 5, 4, 4, 4, 1, 2592, 6, 4, 4, 3, 3, 3, &
181 2700, 5, 5, 4, 3, 3, 3, 2880, 6, 8, 5, 4, 3, 1, &
182 3000, 6, 5, 5, 5, 4, 1, 3072, 6, 8, 4, 4, 4, 1, &
183 3125, 5, 5, 5, 5, 5, 1, 3200, 8, 5, 5, 4, 4, 1, &
184 3240, 6, 5, 4, 3, 3, 3, 3375, 5, 5, 5, 3, 3, 3, &
185 3456, 6, 4, 4, 4, 3, 3, 3600, 6, 8, 5, 5, 3, 1, &
186 3750, 6, 5, 5, 5, 5, 1, 3840, 6, 8, 5, 4, 4, 1, &
187 3888, 6, 8, 3, 3, 3, 3, 4000, 8, 5, 5, 5, 4, 1, &
188 4050, 6, 5, 5, 3, 3, 3, 4096, 8, 8, 4, 4, 4, 1, &
189 4320, 6, 5, 4, 4, 3, 3, 4500, 5, 5, 5, 4, 3, 3, &
190 4608, 6, 8, 8, 4, 3, 1, 4800, 6, 8, 5, 5, 4, 1, &
191 5000, 8, 5, 5, 5, 5, 1, 5120, 8, 8, 5, 4, 4, 1, &
192 5184, 6, 8, 4, 3, 3, 3, 5400, 6, 5, 5, 4, 3, 3, &
193 5625, 5, 5, 5, 5, 3, 3, 5760, 6, 8, 8, 5, 3, 1, &
194 6000, 6, 8, 5, 5, 5, 1, 6144, 6, 8, 8, 4, 4, 1, &
195 6400, 8, 8, 5, 5, 4, 1, 6480, 6, 8, 5, 3, 3, 3, &
196 6750, 6, 5, 5, 5, 3, 3, 6912, 6, 8, 4, 4, 3, 3, &
197 7200, 6, 5, 5, 4, 4, 3, 7500, 5, 5, 5, 5, 4, 3, &
198 7680, 6, 8, 8, 5, 4, 1, 8000, 8, 8, 5, 5, 5, 1, &
199 8192, 8, 8, 8, 4, 4, 1/
200
201 DO i = 1, 150
202 IF (i == 150) THEN
203 print *, 'VALUE OF', n, 'NOT ALLOWED FOR FFT, ALLOWED VALUES ARE:'
20437 FORMAT(15(i5))
205 WRITE (*, 37) (idata(1, j), j=1, 149)
206 cpabort("")
207 END IF
208 IF (n .EQ. idata(1, i)) THEN
209 ic = 0
210 DO j = 1, 6
211 itt = idata(1 + j, i)
212 IF (itt .GT. 1) THEN
213 ic = ic + 1
214 now(j) = idata(1 + j, i)
215 ELSE
216 EXIT
217 END IF
218 END DO
219 EXIT
220 END IF
221 END DO
222
223 after(1) = 1
224 before(ic) = 1
225 DO i = 2, ic
226 after(i) = after(i - 1)*now(i - 1)
227 before(ic - i + 1) = before(ic - i + 2)*now(ic - i + 2)
228 END DO
229
230 twopi = 6.283185307179586_dp
231 angle = isign*twopi/n
232 IF (mod(n, 2) .EQ. 0) THEN
233 nh = n/2
234 trig(1, 1) = 1._dp
235 trig(2, 1) = 0._dp
236 trig(1, nh + 1) = -1._dp
237 trig(2, nh + 1) = 0._dp
238 DO 40, i = 1, nh - 1
239 trigc = cos(i*angle)
240 trigs = sin(i*angle)
241 trig(1, i + 1) = trigc
242 trig(2, i + 1) = trigs
243 trig(1, n - i + 1) = trigc
244 trig(2, n - i + 1) = -trigs
24540 CONTINUE
246 ELSE
247 nh = (n - 1)/2
248 trig(1, 1) = 1._dp
249 trig(2, 1) = 0._dp
250 DO 20, i = 1, nh
251 trigc = cos(i*angle)
252 trigs = sin(i*angle)
253 trig(1, i + 1) = trigc
254 trig(2, i + 1) = trigs
255 trig(1, n - i + 1) = trigc
256 trig(2, n - i + 1) = -trigs
25720 CONTINUE
258 END IF
259
260 END SUBROUTINE ctrig
261
262!ccccccccccccccccccccccccccccccccccccccccccccccc
263
264! **************************************************************************************************
265!> \brief ...
266!> \param mm ...
267!> \param nfft ...
268!> \param m ...
269!> \param nn ...
270!> \param n ...
271!> \param zin ...
272!> \param zout ...
273!> \param trig ...
274!> \param after ...
275!> \param now ...
276!> \param before ...
277!> \param isign ...
278! **************************************************************************************************
279 SUBROUTINE fftstp(mm, nfft, m, nn, n, zin, zout, trig, after, now, before, isign)
280! Copyright (C) Stefan Goedecker, Cornell University, Ithaca, USA, 1994
281! Copyright (C) Stefan Goedecker, MPI Stuttgart, Germany, 1995, 1999
282! This file is distributed under the terms of the
283! GNU General Public License, see http://www.gnu.org/copyleft/gpl.txt .
284
285 INTEGER :: mm, nfft, m, nn, n
286 REAL(kind=dp) :: zin, zout, trig
287 INTEGER :: after, now, before, isign
288
289 INTEGER :: atb, atn, ia, ias, ib, itrig, itt, j, &
290 nin1, nin2, nin3, nin4, nin5, nin6, &
291 nin7, nin8, nout1, nout2, nout3, &
292 nout4, nout5, nout6, nout7, nout8
293 REAL(kind=dp) :: am, ap, bb, bm, bp, ci2, ci3, ci4, ci5, ci6, ci7, ci8, cm, cos2, cos4, cp, &
294 cr2, cr3, cr4, cr5, cr6, cr7, cr8, dm, dpp, r, r1, r2, r25, r3, r34, r4, r5, r6, r7, r8, &
295 rt2i, s, s1, s2, s25, s3, s34, s4, s5, s6, s7, s8, sin2, sin4, ui1, ui2, ui3, ur1, ur2, &
296 ur3, vi1, vi2, vi3, vr1, vr2, vr3
297
298 dimension trig(2, ctrig_length), zin(2, mm, m), zout(2, nn, n)
299 atn = after*now
300 atb = after*before
301
302! sqrt(.5_dp)
303 rt2i = 0.7071067811865475_dp
304 IF (now .EQ. 2) THEN
305 ia = 1
306 nin1 = ia - after
307 nout1 = ia - atn
308 DO ib = 1, before
309 nin1 = nin1 + after
310 nin2 = nin1 + atb
311 nout1 = nout1 + atn
312 nout2 = nout1 + after
313 DO j = 1, nfft
314 r1 = zin(1, j, nin1)
315 s1 = zin(2, j, nin1)
316 r2 = zin(1, j, nin2)
317 s2 = zin(2, j, nin2)
318 zout(1, j, nout1) = r2 + r1
319 zout(2, j, nout1) = s2 + s1
320 zout(1, j, nout2) = r1 - r2
321 zout(2, j, nout2) = s1 - s2
322 END do; END DO
323 DO 2000, ia = 2, after
324 ias = ia - 1
325 IF (2*ias .EQ. after) THEN
326 IF (isign .EQ. 1) THEN
327 nin1 = ia - after
328 nout1 = ia - atn
329 DO ib = 1, before
330 nin1 = nin1 + after
331 nin2 = nin1 + atb
332 nout1 = nout1 + atn
333 nout2 = nout1 + after
334 DO j = 1, nfft
335 r1 = zin(1, j, nin1)
336 s1 = zin(2, j, nin1)
337 r2 = zin(2, j, nin2)
338 s2 = zin(1, j, nin2)
339 zout(1, j, nout1) = r1 - r2
340 zout(2, j, nout1) = s2 + s1
341 zout(1, j, nout2) = r2 + r1
342 zout(2, j, nout2) = s1 - s2
343 END do; END DO
344 ELSE
345 nin1 = ia - after
346 nout1 = ia - atn
347 DO ib = 1, before
348 nin1 = nin1 + after
349 nin2 = nin1 + atb
350 nout1 = nout1 + atn
351 nout2 = nout1 + after
352 DO j = 1, nfft
353 r1 = zin(1, j, nin1)
354 s1 = zin(2, j, nin1)
355 r2 = zin(2, j, nin2)
356 s2 = zin(1, j, nin2)
357 zout(1, j, nout1) = r2 + r1
358 zout(2, j, nout1) = s1 - s2
359 zout(1, j, nout2) = r1 - r2
360 zout(2, j, nout2) = s2 + s1
361 END do; END DO
362 END IF
363 ELSE IF (4*ias .EQ. after) THEN
364 IF (isign .EQ. 1) THEN
365 nin1 = ia - after
366 nout1 = ia - atn
367 DO ib = 1, before
368 nin1 = nin1 + after
369 nin2 = nin1 + atb
370 nout1 = nout1 + atn
371 nout2 = nout1 + after
372 DO j = 1, nfft
373 r1 = zin(1, j, nin1)
374 s1 = zin(2, j, nin1)
375 r = zin(1, j, nin2)
376 s = zin(2, j, nin2)
377 r2 = (r - s)*rt2i
378 s2 = (r + s)*rt2i
379 zout(1, j, nout1) = r2 + r1
380 zout(2, j, nout1) = s2 + s1
381 zout(1, j, nout2) = r1 - r2
382 zout(2, j, nout2) = s1 - s2
383 END do; END DO
384 ELSE
385 nin1 = ia - after
386 nout1 = ia - atn
387 DO ib = 1, before
388 nin1 = nin1 + after
389 nin2 = nin1 + atb
390 nout1 = nout1 + atn
391 nout2 = nout1 + after
392 DO j = 1, nfft
393 r1 = zin(1, j, nin1)
394 s1 = zin(2, j, nin1)
395 r = zin(1, j, nin2)
396 s = zin(2, j, nin2)
397 r2 = (r + s)*rt2i
398 s2 = (s - r)*rt2i
399 zout(1, j, nout1) = r2 + r1
400 zout(2, j, nout1) = s2 + s1
401 zout(1, j, nout2) = r1 - r2
402 zout(2, j, nout2) = s1 - s2
403 END do; END DO
404 END IF
405 ELSE IF (4*ias .EQ. 3*after) THEN
406 IF (isign .EQ. 1) THEN
407 nin1 = ia - after
408 nout1 = ia - atn
409 DO ib = 1, before
410 nin1 = nin1 + after
411 nin2 = nin1 + atb
412 nout1 = nout1 + atn
413 nout2 = nout1 + after
414 DO j = 1, nfft
415 r1 = zin(1, j, nin1)
416 s1 = zin(2, j, nin1)
417 r = zin(1, j, nin2)
418 s = zin(2, j, nin2)
419 r2 = (r + s)*rt2i
420 s2 = (r - s)*rt2i
421 zout(1, j, nout1) = r1 - r2
422 zout(2, j, nout1) = s2 + s1
423 zout(1, j, nout2) = r2 + r1
424 zout(2, j, nout2) = s1 - s2
425 END do; END DO
426 ELSE
427 nin1 = ia - after
428 nout1 = ia - atn
429 DO ib = 1, before
430 nin1 = nin1 + after
431 nin2 = nin1 + atb
432 nout1 = nout1 + atn
433 nout2 = nout1 + after
434 DO j = 1, nfft
435 r1 = zin(1, j, nin1)
436 s1 = zin(2, j, nin1)
437 r = zin(1, j, nin2)
438 s = zin(2, j, nin2)
439 r2 = (s - r)*rt2i
440 s2 = (r + s)*rt2i
441 zout(1, j, nout1) = r2 + r1
442 zout(2, j, nout1) = s1 - s2
443 zout(1, j, nout2) = r1 - r2
444 zout(2, j, nout2) = s2 + s1
445 END do; END DO
446 END IF
447 ELSE
448 itrig = ias*before + 1
449 cr2 = trig(1, itrig)
450 ci2 = trig(2, itrig)
451 nin1 = ia - after
452 nout1 = ia - atn
453 DO ib = 1, before
454 nin1 = nin1 + after
455 nin2 = nin1 + atb
456 nout1 = nout1 + atn
457 nout2 = nout1 + after
458 DO j = 1, nfft
459 r1 = zin(1, j, nin1)
460 s1 = zin(2, j, nin1)
461 r = zin(1, j, nin2)
462 s = zin(2, j, nin2)
463 r2 = r*cr2 - s*ci2
464 s2 = r*ci2 + s*cr2
465 zout(1, j, nout1) = r2 + r1
466 zout(2, j, nout1) = s2 + s1
467 zout(1, j, nout2) = r1 - r2
468 zout(2, j, nout2) = s1 - s2
469 END do; END DO
470 END IF
4712000 CONTINUE
472 ELSE IF (now .EQ. 4) THEN
473 IF (isign .EQ. 1) THEN
474 ia = 1
475 nin1 = ia - after
476 nout1 = ia - atn
477 DO ib = 1, before
478 nin1 = nin1 + after
479 nin2 = nin1 + atb
480 nin3 = nin2 + atb
481 nin4 = nin3 + atb
482 nout1 = nout1 + atn
483 nout2 = nout1 + after
484 nout3 = nout2 + after
485 nout4 = nout3 + after
486 DO j = 1, nfft
487 r1 = zin(1, j, nin1)
488 s1 = zin(2, j, nin1)
489 r2 = zin(1, j, nin2)
490 s2 = zin(2, j, nin2)
491 r3 = zin(1, j, nin3)
492 s3 = zin(2, j, nin3)
493 r4 = zin(1, j, nin4)
494 s4 = zin(2, j, nin4)
495 r = r1 + r3
496 s = r2 + r4
497 zout(1, j, nout1) = r + s
498 zout(1, j, nout3) = r - s
499 r = r1 - r3
500 s = s2 - s4
501 zout(1, j, nout2) = r - s
502 zout(1, j, nout4) = r + s
503 r = s1 + s3
504 s = s2 + s4
505 zout(2, j, nout1) = r + s
506 zout(2, j, nout3) = r - s
507 r = s1 - s3
508 s = r2 - r4
509 zout(2, j, nout2) = r + s
510 zout(2, j, nout4) = r - s
511 END do; END DO
512 DO 4000, ia = 2, after
513 ias = ia - 1
514 IF (2*ias .EQ. after) THEN
515 nin1 = ia - after
516 nout1 = ia - atn
517 DO ib = 1, before
518 nin1 = nin1 + after
519 nin2 = nin1 + atb
520 nin3 = nin2 + atb
521 nin4 = nin3 + atb
522 nout1 = nout1 + atn
523 nout2 = nout1 + after
524 nout3 = nout2 + after
525 nout4 = nout3 + after
526 DO j = 1, nfft
527 r1 = zin(1, j, nin1)
528 s1 = zin(2, j, nin1)
529 r = zin(1, j, nin2)
530 s = zin(2, j, nin2)
531 r2 = (r - s)*rt2i
532 s2 = (r + s)*rt2i
533 r3 = zin(2, j, nin3)
534 s3 = zin(1, j, nin3)
535 r = zin(1, j, nin4)
536 s = zin(2, j, nin4)
537 r4 = (r + s)*rt2i
538 s4 = (r - s)*rt2i
539 r = r1 - r3
540 s = r2 - r4
541 zout(1, j, nout1) = r + s
542 zout(1, j, nout3) = r - s
543 r = r1 + r3
544 s = s2 - s4
545 zout(1, j, nout2) = r - s
546 zout(1, j, nout4) = r + s
547 r = s1 + s3
548 s = s2 + s4
549 zout(2, j, nout1) = r + s
550 zout(2, j, nout3) = r - s
551 r = s1 - s3
552 s = r2 + r4
553 zout(2, j, nout2) = r + s
554 zout(2, j, nout4) = r - s
555 END do; END DO
556 ELSE
557 itt = ias*before
558 itrig = itt + 1
559 cr2 = trig(1, itrig)
560 ci2 = trig(2, itrig)
561 itrig = itrig + itt
562 cr3 = trig(1, itrig)
563 ci3 = trig(2, itrig)
564 itrig = itrig + itt
565 cr4 = trig(1, itrig)
566 ci4 = trig(2, itrig)
567 nin1 = ia - after
568 nout1 = ia - atn
569 DO ib = 1, before
570 nin1 = nin1 + after
571 nin2 = nin1 + atb
572 nin3 = nin2 + atb
573 nin4 = nin3 + atb
574 nout1 = nout1 + atn
575 nout2 = nout1 + after
576 nout3 = nout2 + after
577 nout4 = nout3 + after
578 DO j = 1, nfft
579 r1 = zin(1, j, nin1)
580 s1 = zin(2, j, nin1)
581 r = zin(1, j, nin2)
582 s = zin(2, j, nin2)
583 r2 = r*cr2 - s*ci2
584 s2 = r*ci2 + s*cr2
585 r = zin(1, j, nin3)
586 s = zin(2, j, nin3)
587 r3 = r*cr3 - s*ci3
588 s3 = r*ci3 + s*cr3
589 r = zin(1, j, nin4)
590 s = zin(2, j, nin4)
591 r4 = r*cr4 - s*ci4
592 s4 = r*ci4 + s*cr4
593 r = r1 + r3
594 s = r2 + r4
595 zout(1, j, nout1) = r + s
596 zout(1, j, nout3) = r - s
597 r = r1 - r3
598 s = s2 - s4
599 zout(1, j, nout2) = r - s
600 zout(1, j, nout4) = r + s
601 r = s1 + s3
602 s = s2 + s4
603 zout(2, j, nout1) = r + s
604 zout(2, j, nout3) = r - s
605 r = s1 - s3
606 s = r2 - r4
607 zout(2, j, nout2) = r + s
608 zout(2, j, nout4) = r - s
609 END do; END DO
610 END IF
6114000 CONTINUE
612 ELSE
613 ia = 1
614 nin1 = ia - after
615 nout1 = ia - atn
616 DO ib = 1, before
617 nin1 = nin1 + after
618 nin2 = nin1 + atb
619 nin3 = nin2 + atb
620 nin4 = nin3 + atb
621 nout1 = nout1 + atn
622 nout2 = nout1 + after
623 nout3 = nout2 + after
624 nout4 = nout3 + after
625 DO j = 1, nfft
626 r1 = zin(1, j, nin1)
627 s1 = zin(2, j, nin1)
628 r2 = zin(1, j, nin2)
629 s2 = zin(2, j, nin2)
630 r3 = zin(1, j, nin3)
631 s3 = zin(2, j, nin3)
632 r4 = zin(1, j, nin4)
633 s4 = zin(2, j, nin4)
634 r = r1 + r3
635 s = r2 + r4
636 zout(1, j, nout1) = r + s
637 zout(1, j, nout3) = r - s
638 r = r1 - r3
639 s = s2 - s4
640 zout(1, j, nout2) = r + s
641 zout(1, j, nout4) = r - s
642 r = s1 + s3
643 s = s2 + s4
644 zout(2, j, nout1) = r + s
645 zout(2, j, nout3) = r - s
646 r = s1 - s3
647 s = r2 - r4
648 zout(2, j, nout2) = r - s
649 zout(2, j, nout4) = r + s
650 END do; END DO
651 DO 4100, ia = 2, after
652 ias = ia - 1
653 IF (2*ias .EQ. after) THEN
654 nin1 = ia - after
655 nout1 = ia - atn
656 DO ib = 1, before
657 nin1 = nin1 + after
658 nin2 = nin1 + atb
659 nin3 = nin2 + atb
660 nin4 = nin3 + atb
661 nout1 = nout1 + atn
662 nout2 = nout1 + after
663 nout3 = nout2 + after
664 nout4 = nout3 + after
665 DO j = 1, nfft
666 r1 = zin(1, j, nin1)
667 s1 = zin(2, j, nin1)
668 r = zin(1, j, nin2)
669 s = zin(2, j, nin2)
670 r2 = (r + s)*rt2i
671 s2 = (s - r)*rt2i
672 r3 = zin(2, j, nin3)
673 s3 = zin(1, j, nin3)
674 r = zin(1, j, nin4)
675 s = zin(2, j, nin4)
676 r4 = (s - r)*rt2i
677 s4 = (r + s)*rt2i
678 r = r1 + r3
679 s = r2 + r4
680 zout(1, j, nout1) = r + s
681 zout(1, j, nout3) = r - s
682 r = r1 - r3
683 s = s2 + s4
684 zout(1, j, nout2) = r + s
685 zout(1, j, nout4) = r - s
686 r = s1 - s3
687 s = s2 - s4
688 zout(2, j, nout1) = r + s
689 zout(2, j, nout3) = r - s
690 r = s1 + s3
691 s = r2 - r4
692 zout(2, j, nout2) = r - s
693 zout(2, j, nout4) = r + s
694 END do; END DO
695 ELSE
696 itt = ias*before
697 itrig = itt + 1
698 cr2 = trig(1, itrig)
699 ci2 = trig(2, itrig)
700 itrig = itrig + itt
701 cr3 = trig(1, itrig)
702 ci3 = trig(2, itrig)
703 itrig = itrig + itt
704 cr4 = trig(1, itrig)
705 ci4 = trig(2, itrig)
706 nin1 = ia - after
707 nout1 = ia - atn
708 DO ib = 1, before
709 nin1 = nin1 + after
710 nin2 = nin1 + atb
711 nin3 = nin2 + atb
712 nin4 = nin3 + atb
713 nout1 = nout1 + atn
714 nout2 = nout1 + after
715 nout3 = nout2 + after
716 nout4 = nout3 + after
717 DO j = 1, nfft
718 r1 = zin(1, j, nin1)
719 s1 = zin(2, j, nin1)
720 r = zin(1, j, nin2)
721 s = zin(2, j, nin2)
722 r2 = r*cr2 - s*ci2
723 s2 = r*ci2 + s*cr2
724 r = zin(1, j, nin3)
725 s = zin(2, j, nin3)
726 r3 = r*cr3 - s*ci3
727 s3 = r*ci3 + s*cr3
728 r = zin(1, j, nin4)
729 s = zin(2, j, nin4)
730 r4 = r*cr4 - s*ci4
731 s4 = r*ci4 + s*cr4
732 r = r1 + r3
733 s = r2 + r4
734 zout(1, j, nout1) = r + s
735 zout(1, j, nout3) = r - s
736 r = r1 - r3
737 s = s2 - s4
738 zout(1, j, nout2) = r + s
739 zout(1, j, nout4) = r - s
740 r = s1 + s3
741 s = s2 + s4
742 zout(2, j, nout1) = r + s
743 zout(2, j, nout3) = r - s
744 r = s1 - s3
745 s = r2 - r4
746 zout(2, j, nout2) = r - s
747 zout(2, j, nout4) = r + s
748 END do; END DO
749 END IF
7504100 CONTINUE
751 END IF
752 ELSE IF (now .EQ. 8) THEN
753 IF (isign .EQ. -1) THEN
754 ia = 1
755 nin1 = ia - after
756 nout1 = ia - atn
757 DO ib = 1, before
758 nin1 = nin1 + after
759 nin2 = nin1 + atb
760 nin3 = nin2 + atb
761 nin4 = nin3 + atb
762 nin5 = nin4 + atb
763 nin6 = nin5 + atb
764 nin7 = nin6 + atb
765 nin8 = nin7 + atb
766 nout1 = nout1 + atn
767 nout2 = nout1 + after
768 nout3 = nout2 + after
769 nout4 = nout3 + after
770 nout5 = nout4 + after
771 nout6 = nout5 + after
772 nout7 = nout6 + after
773 nout8 = nout7 + after
774 DO j = 1, nfft
775 r1 = zin(1, j, nin1)
776 s1 = zin(2, j, nin1)
777 r2 = zin(1, j, nin2)
778 s2 = zin(2, j, nin2)
779 r3 = zin(1, j, nin3)
780 s3 = zin(2, j, nin3)
781 r4 = zin(1, j, nin4)
782 s4 = zin(2, j, nin4)
783 r5 = zin(1, j, nin5)
784 s5 = zin(2, j, nin5)
785 r6 = zin(1, j, nin6)
786 s6 = zin(2, j, nin6)
787 r7 = zin(1, j, nin7)
788 s7 = zin(2, j, nin7)
789 r8 = zin(1, j, nin8)
790 s8 = zin(2, j, nin8)
791 r = r1 + r5
792 s = r3 + r7
793 ap = r + s
794 am = r - s
795 r = r2 + r6
796 s = r4 + r8
797 bp = r + s
798 bm = r - s
799 r = s1 + s5
800 s = s3 + s7
801 cp = r + s
802 cm = r - s
803 r = s2 + s6
804 s = s4 + s8
805 dpp = r + s
806 dm = r - s
807 zout(1, j, nout1) = ap + bp
808 zout(2, j, nout1) = cp + dpp
809 zout(1, j, nout5) = ap - bp
810 zout(2, j, nout5) = cp - dpp
811 zout(1, j, nout3) = am + dm
812 zout(2, j, nout3) = cm - bm
813 zout(1, j, nout7) = am - dm
814 zout(2, j, nout7) = cm + bm
815 r = r1 - r5
816 s = s3 - s7
817 ap = r + s
818 am = r - s
819 r = s1 - s5
820 s = r3 - r7
821 bp = r + s
822 bm = r - s
823 r = s4 - s8
824 s = r2 - r6
825 cp = r + s
826 cm = r - s
827 r = s2 - s6
828 s = r4 - r8
829 dpp = r + s
830 dm = r - s
831 r = (cp + dm)*rt2i
832 s = (dm - cp)*rt2i
833 cp = (cm + dpp)*rt2i
834 dpp = (cm - dpp)*rt2i
835 zout(1, j, nout2) = ap + r
836 zout(2, j, nout2) = bm + s
837 zout(1, j, nout6) = ap - r
838 zout(2, j, nout6) = bm - s
839 zout(1, j, nout4) = am + cp
840 zout(2, j, nout4) = bp + dpp
841 zout(1, j, nout8) = am - cp
842 zout(2, j, nout8) = bp - dpp
843 END do; END DO
844 DO 8000, ia = 2, after
845 ias = ia - 1
846 itt = ias*before
847 itrig = itt + 1
848 cr2 = trig(1, itrig)
849 ci2 = trig(2, itrig)
850 itrig = itrig + itt
851 cr3 = trig(1, itrig)
852 ci3 = trig(2, itrig)
853 itrig = itrig + itt
854 cr4 = trig(1, itrig)
855 ci4 = trig(2, itrig)
856 itrig = itrig + itt
857 cr5 = trig(1, itrig)
858 ci5 = trig(2, itrig)
859 itrig = itrig + itt
860 cr6 = trig(1, itrig)
861 ci6 = trig(2, itrig)
862 itrig = itrig + itt
863 cr7 = trig(1, itrig)
864 ci7 = trig(2, itrig)
865 itrig = itrig + itt
866 cr8 = trig(1, itrig)
867 ci8 = trig(2, itrig)
868 nin1 = ia - after
869 nout1 = ia - atn
870 DO ib = 1, before
871 nin1 = nin1 + after
872 nin2 = nin1 + atb
873 nin3 = nin2 + atb
874 nin4 = nin3 + atb
875 nin5 = nin4 + atb
876 nin6 = nin5 + atb
877 nin7 = nin6 + atb
878 nin8 = nin7 + atb
879 nout1 = nout1 + atn
880 nout2 = nout1 + after
881 nout3 = nout2 + after
882 nout4 = nout3 + after
883 nout5 = nout4 + after
884 nout6 = nout5 + after
885 nout7 = nout6 + after
886 nout8 = nout7 + after
887 DO j = 1, nfft
888 r1 = zin(1, j, nin1)
889 s1 = zin(2, j, nin1)
890 r = zin(1, j, nin2)
891 s = zin(2, j, nin2)
892 r2 = r*cr2 - s*ci2
893 s2 = r*ci2 + s*cr2
894 r = zin(1, j, nin3)
895 s = zin(2, j, nin3)
896 r3 = r*cr3 - s*ci3
897 s3 = r*ci3 + s*cr3
898 r = zin(1, j, nin4)
899 s = zin(2, j, nin4)
900 r4 = r*cr4 - s*ci4
901 s4 = r*ci4 + s*cr4
902 r = zin(1, j, nin5)
903 s = zin(2, j, nin5)
904 r5 = r*cr5 - s*ci5
905 s5 = r*ci5 + s*cr5
906 r = zin(1, j, nin6)
907 s = zin(2, j, nin6)
908 r6 = r*cr6 - s*ci6
909 s6 = r*ci6 + s*cr6
910 r = zin(1, j, nin7)
911 s = zin(2, j, nin7)
912 r7 = r*cr7 - s*ci7
913 s7 = r*ci7 + s*cr7
914 r = zin(1, j, nin8)
915 s = zin(2, j, nin8)
916 r8 = r*cr8 - s*ci8
917 s8 = r*ci8 + s*cr8
918 r = r1 + r5
919 s = r3 + r7
920 ap = r + s
921 am = r - s
922 r = r2 + r6
923 s = r4 + r8
924 bp = r + s
925 bm = r - s
926 r = s1 + s5
927 s = s3 + s7
928 cp = r + s
929 cm = r - s
930 r = s2 + s6
931 s = s4 + s8
932 dpp = r + s
933 dm = r - s
934 zout(1, j, nout1) = ap + bp
935 zout(2, j, nout1) = cp + dpp
936 zout(1, j, nout5) = ap - bp
937 zout(2, j, nout5) = cp - dpp
938 zout(1, j, nout3) = am + dm
939 zout(2, j, nout3) = cm - bm
940 zout(1, j, nout7) = am - dm
941 zout(2, j, nout7) = cm + bm
942 r = r1 - r5
943 s = s3 - s7
944 ap = r + s
945 am = r - s
946 r = s1 - s5
947 s = r3 - r7
948 bp = r + s
949 bm = r - s
950 r = s4 - s8
951 s = r2 - r6
952 cp = r + s
953 cm = r - s
954 r = s2 - s6
955 s = r4 - r8
956 dpp = r + s
957 dm = r - s
958 r = (cp + dm)*rt2i
959 s = (dm - cp)*rt2i
960 cp = (cm + dpp)*rt2i
961 dpp = (cm - dpp)*rt2i
962 zout(1, j, nout2) = ap + r
963 zout(2, j, nout2) = bm + s
964 zout(1, j, nout6) = ap - r
965 zout(2, j, nout6) = bm - s
966 zout(1, j, nout4) = am + cp
967 zout(2, j, nout4) = bp + dpp
968 zout(1, j, nout8) = am - cp
969 zout(2, j, nout8) = bp - dpp
970 END do; END DO
9718000 CONTINUE
972
973 ELSE
974 ia = 1
975 nin1 = ia - after
976 nout1 = ia - atn
977 DO ib = 1, before
978 nin1 = nin1 + after
979 nin2 = nin1 + atb
980 nin3 = nin2 + atb
981 nin4 = nin3 + atb
982 nin5 = nin4 + atb
983 nin6 = nin5 + atb
984 nin7 = nin6 + atb
985 nin8 = nin7 + atb
986 nout1 = nout1 + atn
987 nout2 = nout1 + after
988 nout3 = nout2 + after
989 nout4 = nout3 + after
990 nout5 = nout4 + after
991 nout6 = nout5 + after
992 nout7 = nout6 + after
993 nout8 = nout7 + after
994 DO j = 1, nfft
995 r1 = zin(1, j, nin1)
996 s1 = zin(2, j, nin1)
997 r2 = zin(1, j, nin2)
998 s2 = zin(2, j, nin2)
999 r3 = zin(1, j, nin3)
1000 s3 = zin(2, j, nin3)
1001 r4 = zin(1, j, nin4)
1002 s4 = zin(2, j, nin4)
1003 r5 = zin(1, j, nin5)
1004 s5 = zin(2, j, nin5)
1005 r6 = zin(1, j, nin6)
1006 s6 = zin(2, j, nin6)
1007 r7 = zin(1, j, nin7)
1008 s7 = zin(2, j, nin7)
1009 r8 = zin(1, j, nin8)
1010 s8 = zin(2, j, nin8)
1011 r = r1 + r5
1012 s = r3 + r7
1013 ap = r + s
1014 am = r - s
1015 r = r2 + r6
1016 s = r4 + r8
1017 bp = r + s
1018 bm = r - s
1019 r = s1 + s5
1020 s = s3 + s7
1021 cp = r + s
1022 cm = r - s
1023 r = s2 + s6
1024 s = s4 + s8
1025 dpp = r + s
1026 dm = r - s
1027 zout(1, j, nout1) = ap + bp
1028 zout(2, j, nout1) = cp + dpp
1029 zout(1, j, nout5) = ap - bp
1030 zout(2, j, nout5) = cp - dpp
1031 zout(1, j, nout3) = am - dm
1032 zout(2, j, nout3) = cm + bm
1033 zout(1, j, nout7) = am + dm
1034 zout(2, j, nout7) = cm - bm
1035 r = r1 - r5
1036 s = -s3 + s7
1037 ap = r + s
1038 am = r - s
1039 r = s1 - s5
1040 s = r7 - r3
1041 bp = r + s
1042 bm = r - s
1043 r = -s4 + s8
1044 s = r2 - r6
1045 cp = r + s
1046 cm = r - s
1047 r = -s2 + s6
1048 s = r4 - r8
1049 dpp = r + s
1050 dm = r - s
1051 r = (cp + dm)*rt2i
1052 s = (cp - dm)*rt2i
1053 cp = (cm + dpp)*rt2i
1054 dpp = (dpp - cm)*rt2i
1055 zout(1, j, nout2) = ap + r
1056 zout(2, j, nout2) = bm + s
1057 zout(1, j, nout6) = ap - r
1058 zout(2, j, nout6) = bm - s
1059 zout(1, j, nout4) = am + cp
1060 zout(2, j, nout4) = bp + dpp
1061 zout(1, j, nout8) = am - cp
1062 zout(2, j, nout8) = bp - dpp
1063 END do; END DO
1064
1065 DO 8001, ia = 2, after
1066 ias = ia - 1
1067 itt = ias*before
1068 itrig = itt + 1
1069 cr2 = trig(1, itrig)
1070 ci2 = trig(2, itrig)
1071 itrig = itrig + itt
1072 cr3 = trig(1, itrig)
1073 ci3 = trig(2, itrig)
1074 itrig = itrig + itt
1075 cr4 = trig(1, itrig)
1076 ci4 = trig(2, itrig)
1077 itrig = itrig + itt
1078 cr5 = trig(1, itrig)
1079 ci5 = trig(2, itrig)
1080 itrig = itrig + itt
1081 cr6 = trig(1, itrig)
1082 ci6 = trig(2, itrig)
1083 itrig = itrig + itt
1084 cr7 = trig(1, itrig)
1085 ci7 = trig(2, itrig)
1086 itrig = itrig + itt
1087 cr8 = trig(1, itrig)
1088 ci8 = trig(2, itrig)
1089 nin1 = ia - after
1090 nout1 = ia - atn
1091 DO ib = 1, before
1092 nin1 = nin1 + after
1093 nin2 = nin1 + atb
1094 nin3 = nin2 + atb
1095 nin4 = nin3 + atb
1096 nin5 = nin4 + atb
1097 nin6 = nin5 + atb
1098 nin7 = nin6 + atb
1099 nin8 = nin7 + atb
1100 nout1 = nout1 + atn
1101 nout2 = nout1 + after
1102 nout3 = nout2 + after
1103 nout4 = nout3 + after
1104 nout5 = nout4 + after
1105 nout6 = nout5 + after
1106 nout7 = nout6 + after
1107 nout8 = nout7 + after
1108 DO j = 1, nfft
1109 r1 = zin(1, j, nin1)
1110 s1 = zin(2, j, nin1)
1111 r = zin(1, j, nin2)
1112 s = zin(2, j, nin2)
1113 r2 = r*cr2 - s*ci2
1114 s2 = r*ci2 + s*cr2
1115 r = zin(1, j, nin3)
1116 s = zin(2, j, nin3)
1117 r3 = r*cr3 - s*ci3
1118 s3 = r*ci3 + s*cr3
1119 r = zin(1, j, nin4)
1120 s = zin(2, j, nin4)
1121 r4 = r*cr4 - s*ci4
1122 s4 = r*ci4 + s*cr4
1123 r = zin(1, j, nin5)
1124 s = zin(2, j, nin5)
1125 r5 = r*cr5 - s*ci5
1126 s5 = r*ci5 + s*cr5
1127 r = zin(1, j, nin6)
1128 s = zin(2, j, nin6)
1129 r6 = r*cr6 - s*ci6
1130 s6 = r*ci6 + s*cr6
1131 r = zin(1, j, nin7)
1132 s = zin(2, j, nin7)
1133 r7 = r*cr7 - s*ci7
1134 s7 = r*ci7 + s*cr7
1135 r = zin(1, j, nin8)
1136 s = zin(2, j, nin8)
1137 r8 = r*cr8 - s*ci8
1138 s8 = r*ci8 + s*cr8
1139 r = r1 + r5
1140 s = r3 + r7
1141 ap = r + s
1142 am = r - s
1143 r = r2 + r6
1144 s = r4 + r8
1145 bp = r + s
1146 bm = r - s
1147 r = s1 + s5
1148 s = s3 + s7
1149 cp = r + s
1150 cm = r - s
1151 r = s2 + s6
1152 s = s4 + s8
1153 dpp = r + s
1154 dm = r - s
1155 zout(1, j, nout1) = ap + bp
1156 zout(2, j, nout1) = cp + dpp
1157 zout(1, j, nout5) = ap - bp
1158 zout(2, j, nout5) = cp - dpp
1159 zout(1, j, nout3) = am - dm
1160 zout(2, j, nout3) = cm + bm
1161 zout(1, j, nout7) = am + dm
1162 zout(2, j, nout7) = cm - bm
1163 r = r1 - r5
1164 s = -s3 + s7
1165 ap = r + s
1166 am = r - s
1167 r = s1 - s5
1168 s = r7 - r3
1169 bp = r + s
1170 bm = r - s
1171 r = -s4 + s8
1172 s = r2 - r6
1173 cp = r + s
1174 cm = r - s
1175 r = -s2 + s6
1176 s = r4 - r8
1177 dpp = r + s
1178 dm = r - s
1179 r = (cp + dm)*rt2i
1180 s = (cp - dm)*rt2i
1181 cp = (cm + dpp)*rt2i
1182 dpp = (dpp - cm)*rt2i
1183 zout(1, j, nout2) = ap + r
1184 zout(2, j, nout2) = bm + s
1185 zout(1, j, nout6) = ap - r
1186 zout(2, j, nout6) = bm - s
1187 zout(1, j, nout4) = am + cp
1188 zout(2, j, nout4) = bp + dpp
1189 zout(1, j, nout8) = am - cp
1190 zout(2, j, nout8) = bp - dpp
1191 END do; END DO
11928001 CONTINUE
1193
1194 END IF
1195 ELSE IF (now .EQ. 3) THEN
1196! .5_dp*sqrt(3._dp)
1197 bb = isign*0.8660254037844387_dp
1198 ia = 1
1199 nin1 = ia - after
1200 nout1 = ia - atn
1201 DO ib = 1, before
1202 nin1 = nin1 + after
1203 nin2 = nin1 + atb
1204 nin3 = nin2 + atb
1205 nout1 = nout1 + atn
1206 nout2 = nout1 + after
1207 nout3 = nout2 + after
1208 DO j = 1, nfft
1209 r1 = zin(1, j, nin1)
1210 s1 = zin(2, j, nin1)
1211 r2 = zin(1, j, nin2)
1212 s2 = zin(2, j, nin2)
1213 r3 = zin(1, j, nin3)
1214 s3 = zin(2, j, nin3)
1215 r = r2 + r3
1216 s = s2 + s3
1217 zout(1, j, nout1) = r + r1
1218 zout(2, j, nout1) = s + s1
1219 r1 = r1 - .5_dp*r
1220 s1 = s1 - .5_dp*s
1221 r2 = bb*(r2 - r3)
1222 s2 = bb*(s2 - s3)
1223 zout(1, j, nout2) = r1 - s2
1224 zout(2, j, nout2) = s1 + r2
1225 zout(1, j, nout3) = r1 + s2
1226 zout(2, j, nout3) = s1 - r2
1227 END do; END DO
1228 DO 3000, ia = 2, after
1229 ias = ia - 1
1230 IF (4*ias .EQ. 3*after) THEN
1231 IF (isign .EQ. 1) THEN
1232 nin1 = ia - after
1233 nout1 = ia - atn
1234 DO ib = 1, before
1235 nin1 = nin1 + after
1236 nin2 = nin1 + atb
1237 nin3 = nin2 + atb
1238 nout1 = nout1 + atn
1239 nout2 = nout1 + after
1240 nout3 = nout2 + after
1241 DO j = 1, nfft
1242 r1 = zin(1, j, nin1)
1243 s1 = zin(2, j, nin1)
1244 r2 = zin(2, j, nin2)
1245 s2 = zin(1, j, nin2)
1246 r3 = zin(1, j, nin3)
1247 s3 = zin(2, j, nin3)
1248 r = r3 + r2
1249 s = s2 - s3
1250 zout(1, j, nout1) = r1 - r
1251 zout(2, j, nout1) = s + s1
1252 r1 = r1 + .5_dp*r
1253 s1 = s1 - .5_dp*s
1254 r2 = bb*(r2 - r3)
1255 s2 = bb*(s2 + s3)
1256 zout(1, j, nout2) = r1 - s2
1257 zout(2, j, nout2) = s1 - r2
1258 zout(1, j, nout3) = r1 + s2
1259 zout(2, j, nout3) = s1 + r2
1260 END do; END DO
1261 ELSE
1262 nin1 = ia - after
1263 nout1 = ia - atn
1264 DO ib = 1, before
1265 nin1 = nin1 + after
1266 nin2 = nin1 + atb
1267 nin3 = nin2 + atb
1268 nout1 = nout1 + atn
1269 nout2 = nout1 + after
1270 nout3 = nout2 + after
1271 DO j = 1, nfft
1272 r1 = zin(1, j, nin1)
1273 s1 = zin(2, j, nin1)
1274 r2 = zin(2, j, nin2)
1275 s2 = zin(1, j, nin2)
1276 r3 = zin(1, j, nin3)
1277 s3 = zin(2, j, nin3)
1278 r = r2 - r3
1279 s = s2 + s3
1280 zout(1, j, nout1) = r + r1
1281 zout(2, j, nout1) = s1 - s
1282 r1 = r1 - .5_dp*r
1283 s1 = s1 + .5_dp*s
1284 r2 = bb*(r2 + r3)
1285 s2 = bb*(s2 - s3)
1286 zout(1, j, nout2) = r1 + s2
1287 zout(2, j, nout2) = s1 + r2
1288 zout(1, j, nout3) = r1 - s2
1289 zout(2, j, nout3) = s1 - r2
1290 END do; END DO
1291 END IF
1292 ELSE IF (8*ias .EQ. 3*after) THEN
1293 IF (isign .EQ. 1) THEN
1294 nin1 = ia - after
1295 nout1 = ia - atn
1296 DO ib = 1, before
1297 nin1 = nin1 + after
1298 nin2 = nin1 + atb
1299 nin3 = nin2 + atb
1300 nout1 = nout1 + atn
1301 nout2 = nout1 + after
1302 nout3 = nout2 + after
1303 DO j = 1, nfft
1304 r1 = zin(1, j, nin1)
1305 s1 = zin(2, j, nin1)
1306 r = zin(1, j, nin2)
1307 s = zin(2, j, nin2)
1308 r2 = (r - s)*rt2i
1309 s2 = (r + s)*rt2i
1310 r3 = zin(2, j, nin3)
1311 s3 = zin(1, j, nin3)
1312 r = r2 - r3
1313 s = s2 + s3
1314 zout(1, j, nout1) = r + r1
1315 zout(2, j, nout1) = s + s1
1316 r1 = r1 - .5_dp*r
1317 s1 = s1 - .5_dp*s
1318 r2 = bb*(r2 + r3)
1319 s2 = bb*(s2 - s3)
1320 zout(1, j, nout2) = r1 - s2
1321 zout(2, j, nout2) = s1 + r2
1322 zout(1, j, nout3) = r1 + s2
1323 zout(2, j, nout3) = s1 - r2
1324 END do; END DO
1325 ELSE
1326 nin1 = ia - after
1327 nout1 = ia - atn
1328 DO ib = 1, before
1329 nin1 = nin1 + after
1330 nin2 = nin1 + atb
1331 nin3 = nin2 + atb
1332 nout1 = nout1 + atn
1333 nout2 = nout1 + after
1334 nout3 = nout2 + after
1335 DO j = 1, nfft
1336 r1 = zin(1, j, nin1)
1337 s1 = zin(2, j, nin1)
1338 r = zin(1, j, nin2)
1339 s = zin(2, j, nin2)
1340 r2 = (r + s)*rt2i
1341 s2 = (s - r)*rt2i
1342 r3 = zin(2, j, nin3)
1343 s3 = zin(1, j, nin3)
1344 r = r2 + r3
1345 s = s2 - s3
1346 zout(1, j, nout1) = r + r1
1347 zout(2, j, nout1) = s + s1
1348 r1 = r1 - .5_dp*r
1349 s1 = s1 - .5_dp*s
1350 r2 = bb*(r2 - r3)
1351 s2 = bb*(s2 + s3)
1352 zout(1, j, nout2) = r1 - s2
1353 zout(2, j, nout2) = s1 + r2
1354 zout(1, j, nout3) = r1 + s2
1355 zout(2, j, nout3) = s1 - r2
1356 END do; END DO
1357 END IF
1358 ELSE
1359 itt = ias*before
1360 itrig = itt + 1
1361 cr2 = trig(1, itrig)
1362 ci2 = trig(2, itrig)
1363 itrig = itrig + itt
1364 cr3 = trig(1, itrig)
1365 ci3 = trig(2, itrig)
1366 nin1 = ia - after
1367 nout1 = ia - atn
1368 DO ib = 1, before
1369 nin1 = nin1 + after
1370 nin2 = nin1 + atb
1371 nin3 = nin2 + atb
1372 nout1 = nout1 + atn
1373 nout2 = nout1 + after
1374 nout3 = nout2 + after
1375 DO j = 1, nfft
1376 r1 = zin(1, j, nin1)
1377 s1 = zin(2, j, nin1)
1378 r = zin(1, j, nin2)
1379 s = zin(2, j, nin2)
1380 r2 = r*cr2 - s*ci2
1381 s2 = r*ci2 + s*cr2
1382 r = zin(1, j, nin3)
1383 s = zin(2, j, nin3)
1384 r3 = r*cr3 - s*ci3
1385 s3 = r*ci3 + s*cr3
1386 r = r2 + r3
1387 s = s2 + s3
1388 zout(1, j, nout1) = r + r1
1389 zout(2, j, nout1) = s + s1
1390 r1 = r1 - .5_dp*r
1391 s1 = s1 - .5_dp*s
1392 r2 = bb*(r2 - r3)
1393 s2 = bb*(s2 - s3)
1394 zout(1, j, nout2) = r1 - s2
1395 zout(2, j, nout2) = s1 + r2
1396 zout(1, j, nout3) = r1 + s2
1397 zout(2, j, nout3) = s1 - r2
1398 END do; END DO
1399 END IF
14003000 CONTINUE
1401 ELSE IF (now .EQ. 5) THEN
1402! cos(2._dp*pi/5._dp)
1403 cos2 = 0.3090169943749474_dp
1404! cos(4._dp*pi/5._dp)
1405 cos4 = -0.8090169943749474_dp
1406! sin(2._dp*pi/5._dp)
1407 sin2 = isign*0.9510565162951536_dp
1408! sin(4._dp*pi/5._dp)
1409 sin4 = isign*0.5877852522924731_dp
1410 ia = 1
1411 nin1 = ia - after
1412 nout1 = ia - atn
1413 DO ib = 1, before
1414 nin1 = nin1 + after
1415 nin2 = nin1 + atb
1416 nin3 = nin2 + atb
1417 nin4 = nin3 + atb
1418 nin5 = nin4 + atb
1419 nout1 = nout1 + atn
1420 nout2 = nout1 + after
1421 nout3 = nout2 + after
1422 nout4 = nout3 + after
1423 nout5 = nout4 + after
1424 DO j = 1, nfft
1425 r1 = zin(1, j, nin1)
1426 s1 = zin(2, j, nin1)
1427 r2 = zin(1, j, nin2)
1428 s2 = zin(2, j, nin2)
1429 r3 = zin(1, j, nin3)
1430 s3 = zin(2, j, nin3)
1431 r4 = zin(1, j, nin4)
1432 s4 = zin(2, j, nin4)
1433 r5 = zin(1, j, nin5)
1434 s5 = zin(2, j, nin5)
1435 r25 = r2 + r5
1436 r34 = r3 + r4
1437 s25 = s2 - s5
1438 s34 = s3 - s4
1439 zout(1, j, nout1) = r1 + r25 + r34
1440 r = r1 + cos2*r25 + cos4*r34
1441 s = sin2*s25 + sin4*s34
1442 zout(1, j, nout2) = r - s
1443 zout(1, j, nout5) = r + s
1444 r = r1 + cos4*r25 + cos2*r34
1445 s = sin4*s25 - sin2*s34
1446 zout(1, j, nout3) = r - s
1447 zout(1, j, nout4) = r + s
1448 r25 = r2 - r5
1449 r34 = r3 - r4
1450 s25 = s2 + s5
1451 s34 = s3 + s4
1452 zout(2, j, nout1) = s1 + s25 + s34
1453 r = s1 + cos2*s25 + cos4*s34
1454 s = sin2*r25 + sin4*r34
1455 zout(2, j, nout2) = r + s
1456 zout(2, j, nout5) = r - s
1457 r = s1 + cos4*s25 + cos2*s34
1458 s = sin4*r25 - sin2*r34
1459 zout(2, j, nout3) = r + s
1460 zout(2, j, nout4) = r - s
1461 END do; END DO
1462 DO 5000, ia = 2, after
1463 ias = ia - 1
1464 IF (8*ias .EQ. 5*after) THEN
1465 IF (isign .EQ. 1) THEN
1466 nin1 = ia - after
1467 nout1 = ia - atn
1468 DO ib = 1, before
1469 nin1 = nin1 + after
1470 nin2 = nin1 + atb
1471 nin3 = nin2 + atb
1472 nin4 = nin3 + atb
1473 nin5 = nin4 + atb
1474 nout1 = nout1 + atn
1475 nout2 = nout1 + after
1476 nout3 = nout2 + after
1477 nout4 = nout3 + after
1478 nout5 = nout4 + after
1479 DO j = 1, nfft
1480 r1 = zin(1, j, nin1)
1481 s1 = zin(2, j, nin1)
1482 r = zin(1, j, nin2)
1483 s = zin(2, j, nin2)
1484 r2 = (r - s)*rt2i
1485 s2 = (r + s)*rt2i
1486 r3 = zin(2, j, nin3)
1487 s3 = zin(1, j, nin3)
1488 r = zin(1, j, nin4)
1489 s = zin(2, j, nin4)
1490 r4 = (r + s)*rt2i
1491 s4 = (r - s)*rt2i
1492 r5 = zin(1, j, nin5)
1493 s5 = zin(2, j, nin5)
1494 r25 = r2 - r5
1495 r34 = r3 + r4
1496 s25 = s2 + s5
1497 s34 = s3 - s4
1498 zout(1, j, nout1) = r1 + r25 - r34
1499 r = r1 + cos2*r25 - cos4*r34
1500 s = sin2*s25 + sin4*s34
1501 zout(1, j, nout2) = r - s
1502 zout(1, j, nout5) = r + s
1503 r = r1 + cos4*r25 - cos2*r34
1504 s = sin4*s25 - sin2*s34
1505 zout(1, j, nout3) = r - s
1506 zout(1, j, nout4) = r + s
1507 r25 = r2 + r5
1508 r34 = r4 - r3
1509 s25 = s2 - s5
1510 s34 = s3 + s4
1511 zout(2, j, nout1) = s1 + s25 + s34
1512 r = s1 + cos2*s25 + cos4*s34
1513 s = sin2*r25 + sin4*r34
1514 zout(2, j, nout2) = r + s
1515 zout(2, j, nout5) = r - s
1516 r = s1 + cos4*s25 + cos2*s34
1517 s = sin4*r25 - sin2*r34
1518 zout(2, j, nout3) = r + s
1519 zout(2, j, nout4) = r - s
1520 END do; END DO
1521 ELSE
1522 nin1 = ia - after
1523 nout1 = ia - atn
1524 DO ib = 1, before
1525 nin1 = nin1 + after
1526 nin2 = nin1 + atb
1527 nin3 = nin2 + atb
1528 nin4 = nin3 + atb
1529 nin5 = nin4 + atb
1530 nout1 = nout1 + atn
1531 nout2 = nout1 + after
1532 nout3 = nout2 + after
1533 nout4 = nout3 + after
1534 nout5 = nout4 + after
1535 DO j = 1, nfft
1536 r1 = zin(1, j, nin1)
1537 s1 = zin(2, j, nin1)
1538 r = zin(1, j, nin2)
1539 s = zin(2, j, nin2)
1540 r2 = (r + s)*rt2i
1541 s2 = (s - r)*rt2i
1542 r3 = zin(2, j, nin3)
1543 s3 = zin(1, j, nin3)
1544 r = zin(1, j, nin4)
1545 s = zin(2, j, nin4)
1546 r4 = (s - r)*rt2i
1547 s4 = (r + s)*rt2i
1548 r5 = zin(1, j, nin5)
1549 s5 = zin(2, j, nin5)
1550 r25 = r2 - r5
1551 r34 = r3 + r4
1552 s25 = s2 + s5
1553 s34 = s4 - s3
1554 zout(1, j, nout1) = r1 + r25 + r34
1555 r = r1 + cos2*r25 + cos4*r34
1556 s = sin2*s25 + sin4*s34
1557 zout(1, j, nout2) = r - s
1558 zout(1, j, nout5) = r + s
1559 r = r1 + cos4*r25 + cos2*r34
1560 s = sin4*s25 - sin2*s34
1561 zout(1, j, nout3) = r - s
1562 zout(1, j, nout4) = r + s
1563 r25 = r2 + r5
1564 r34 = r3 - r4
1565 s25 = s2 - s5
1566 s34 = s3 + s4
1567 zout(2, j, nout1) = s1 + s25 - s34
1568 r = s1 + cos2*s25 - cos4*s34
1569 s = sin2*r25 + sin4*r34
1570 zout(2, j, nout2) = r + s
1571 zout(2, j, nout5) = r - s
1572 r = s1 + cos4*s25 - cos2*s34
1573 s = sin4*r25 - sin2*r34
1574 zout(2, j, nout3) = r + s
1575 zout(2, j, nout4) = r - s
1576 END do; END DO
1577 END IF
1578 ELSE
1579 ias = ia - 1
1580 itt = ias*before
1581 itrig = itt + 1
1582 cr2 = trig(1, itrig)
1583 ci2 = trig(2, itrig)
1584 itrig = itrig + itt
1585 cr3 = trig(1, itrig)
1586 ci3 = trig(2, itrig)
1587 itrig = itrig + itt
1588 cr4 = trig(1, itrig)
1589 ci4 = trig(2, itrig)
1590 itrig = itrig + itt
1591 cr5 = trig(1, itrig)
1592 ci5 = trig(2, itrig)
1593 nin1 = ia - after
1594 nout1 = ia - atn
1595 DO ib = 1, before
1596 nin1 = nin1 + after
1597 nin2 = nin1 + atb
1598 nin3 = nin2 + atb
1599 nin4 = nin3 + atb
1600 nin5 = nin4 + atb
1601 nout1 = nout1 + atn
1602 nout2 = nout1 + after
1603 nout3 = nout2 + after
1604 nout4 = nout3 + after
1605 nout5 = nout4 + after
1606 DO j = 1, nfft
1607 r1 = zin(1, j, nin1)
1608 s1 = zin(2, j, nin1)
1609 r = zin(1, j, nin2)
1610 s = zin(2, j, nin2)
1611 r2 = r*cr2 - s*ci2
1612 s2 = r*ci2 + s*cr2
1613 r = zin(1, j, nin3)
1614 s = zin(2, j, nin3)
1615 r3 = r*cr3 - s*ci3
1616 s3 = r*ci3 + s*cr3
1617 r = zin(1, j, nin4)
1618 s = zin(2, j, nin4)
1619 r4 = r*cr4 - s*ci4
1620 s4 = r*ci4 + s*cr4
1621 r = zin(1, j, nin5)
1622 s = zin(2, j, nin5)
1623 r5 = r*cr5 - s*ci5
1624 s5 = r*ci5 + s*cr5
1625 r25 = r2 + r5
1626 r34 = r3 + r4
1627 s25 = s2 - s5
1628 s34 = s3 - s4
1629 zout(1, j, nout1) = r1 + r25 + r34
1630 r = r1 + cos2*r25 + cos4*r34
1631 s = sin2*s25 + sin4*s34
1632 zout(1, j, nout2) = r - s
1633 zout(1, j, nout5) = r + s
1634 r = r1 + cos4*r25 + cos2*r34
1635 s = sin4*s25 - sin2*s34
1636 zout(1, j, nout3) = r - s
1637 zout(1, j, nout4) = r + s
1638 r25 = r2 - r5
1639 r34 = r3 - r4
1640 s25 = s2 + s5
1641 s34 = s3 + s4
1642 zout(2, j, nout1) = s1 + s25 + s34
1643 r = s1 + cos2*s25 + cos4*s34
1644 s = sin2*r25 + sin4*r34
1645 zout(2, j, nout2) = r + s
1646 zout(2, j, nout5) = r - s
1647 r = s1 + cos4*s25 + cos2*s34
1648 s = sin4*r25 - sin2*r34
1649 zout(2, j, nout3) = r + s
1650 zout(2, j, nout4) = r - s
1651 END do; END DO
1652 END IF
16535000 CONTINUE
1654 ELSE IF (now .EQ. 6) THEN
1655! .5_dp*sqrt(3._dp)
1656 bb = isign*0.8660254037844387_dp
1657
1658 ia = 1
1659 nin1 = ia - after
1660 nout1 = ia - atn
1661 DO ib = 1, before
1662 nin1 = nin1 + after
1663 nin2 = nin1 + atb
1664 nin3 = nin2 + atb
1665 nin4 = nin3 + atb
1666 nin5 = nin4 + atb
1667 nin6 = nin5 + atb
1668 nout1 = nout1 + atn
1669 nout2 = nout1 + after
1670 nout3 = nout2 + after
1671 nout4 = nout3 + after
1672 nout5 = nout4 + after
1673 nout6 = nout5 + after
1674 DO j = 1, nfft
1675 r2 = zin(1, j, nin3)
1676 s2 = zin(2, j, nin3)
1677 r3 = zin(1, j, nin5)
1678 s3 = zin(2, j, nin5)
1679 r = r2 + r3
1680 s = s2 + s3
1681 r1 = zin(1, j, nin1)
1682 s1 = zin(2, j, nin1)
1683 ur1 = r + r1
1684 ui1 = s + s1
1685 r1 = r1 - .5_dp*r
1686 s1 = s1 - .5_dp*s
1687 r = r2 - r3
1688 s = s2 - s3
1689 ur2 = r1 - s*bb
1690 ui2 = s1 + r*bb
1691 ur3 = r1 + s*bb
1692 ui3 = s1 - r*bb
1693
1694 r2 = zin(1, j, nin6)
1695 s2 = zin(2, j, nin6)
1696 r3 = zin(1, j, nin2)
1697 s3 = zin(2, j, nin2)
1698 r = r2 + r3
1699 s = s2 + s3
1700 r1 = zin(1, j, nin4)
1701 s1 = zin(2, j, nin4)
1702 vr1 = r + r1
1703 vi1 = s + s1
1704 r1 = r1 - .5_dp*r
1705 s1 = s1 - .5_dp*s
1706 r = r2 - r3
1707 s = s2 - s3
1708 vr2 = r1 - s*bb
1709 vi2 = s1 + r*bb
1710 vr3 = r1 + s*bb
1711 vi3 = s1 - r*bb
1712
1713 zout(1, j, nout1) = ur1 + vr1
1714 zout(2, j, nout1) = ui1 + vi1
1715 zout(1, j, nout5) = ur2 + vr2
1716 zout(2, j, nout5) = ui2 + vi2
1717 zout(1, j, nout3) = ur3 + vr3
1718 zout(2, j, nout3) = ui3 + vi3
1719 zout(1, j, nout4) = ur1 - vr1
1720 zout(2, j, nout4) = ui1 - vi1
1721 zout(1, j, nout2) = ur2 - vr2
1722 zout(2, j, nout2) = ui2 - vi2
1723 zout(1, j, nout6) = ur3 - vr3
1724 zout(2, j, nout6) = ui3 - vi3
1725 END do; END DO
1726 ELSE
1727 cpabort("error fftstp")
1728 END IF
1729
1730 END SUBROUTINE fftstp
1731
1732 END MODULE ps_wavelet_fft3d
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
ps_wavelet_fft3d
Definition ps_wavelet_fft3d.F:9
ps_wavelet_fft3d::ctrig_length
integer, parameter, public ctrig_length
Definition ps_wavelet_fft3d.F:20
ps_wavelet_fft3d::fftstp
subroutine, public fftstp(mm, nfft, m, nn, n, zin, zout, trig, after, now, before, isign)
...
Definition ps_wavelet_fft3d.F:280
ps_wavelet_fft3d::ctrig
subroutine, public ctrig(n, trig, after, before, now, isign, ic)
...
Definition ps_wavelet_fft3d.F:105
ps_wavelet_fft3d::fourier_dim
subroutine, public fourier_dim(n, n_next)
Give a number n_next > n compatible for the FFT.
Definition ps_wavelet_fft3d.F:34