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