(git:34ef472)
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 ! **************************************************************************************************
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, &
25 
26 CONTAINS
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:'
204 37 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
245 40 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
257 20 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
471 2000 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
611 4000 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
750 4100 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
971 8000 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
1192 8001 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
1400 3000 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
1653 5000 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
pure real(kind=dp) function angle(a, b)
Calculation of the angle between the vectors a and b. The angle is returned in radians.
Definition: dumpdcd.F:1008
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
integer, parameter, public ctrig_length
subroutine, public fftstp(mm, nfft, m, nn, n, zin, zout, trig, after, now, before, isign)
...
subroutine, public ctrig(n, trig, after, before, now, isign, ic)
...
subroutine, public fourier_dim(n, n_next)
Give a number n_next > n compatible for the FFT.