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