Line data Source code
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 : ! **************************************************************************************************
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 110141 : 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 1782710 : loop_data: DO i = 1, ndata1024
54 1782710 : IF (n <= idata(i)) THEN
55 110141 : n_next = idata(i)
56 110141 : RETURN
57 : END IF
58 : END DO loop_data
59 0 : WRITE (unit=*, fmt=*) "fourier_dim: ", n, " is bigger than ", idata(ndata1024)
60 0 : 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 97037 : 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 1583950 : DO i = 1, 150
202 1583950 : IF (i == 150) THEN
203 0 : PRINT *, 'VALUE OF', n, 'NOT ALLOWED FOR FFT, ALLOWED VALUES ARE:'
204 : 37 FORMAT(15(i5))
205 0 : WRITE (*, 37) (idata(1, j), j=1, 149)
206 0 : CPABORT("")
207 : END IF
208 1583950 : IF (n .EQ. idata(1, i)) THEN
209 97037 : ic = 0
210 331757 : DO j = 1, 6
211 331757 : itt = idata(1 + j, i)
212 331757 : IF (itt .GT. 1) THEN
213 234720 : ic = ic + 1
214 234720 : 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 97037 : after(1) = 1
224 97037 : before(ic) = 1
225 234720 : DO i = 2, ic
226 137683 : after(i) = after(i - 1)*now(i - 1)
227 234720 : before(ic - i + 1) = before(ic - i + 2)*now(ic - i + 2)
228 : END DO
229 :
230 97037 : twopi = 6.283185307179586_dp
231 97037 : angle = isign*twopi/n
232 97037 : IF (MOD(n, 2) .EQ. 0) THEN
233 85622 : nh = n/2
234 85622 : trig(1, 1) = 1._dp
235 85622 : trig(2, 1) = 0._dp
236 85622 : trig(1, nh + 1) = -1._dp
237 85622 : trig(2, nh + 1) = 0._dp
238 1884038 : DO 40, i = 1, nh - 1
239 1798416 : trigc = COS(i*angle)
240 1798416 : trigs = SIN(i*angle)
241 1798416 : trig(1, i + 1) = trigc
242 1798416 : trig(2, i + 1) = trigs
243 1798416 : trig(1, n - i + 1) = trigc
244 1798416 : trig(2, n - i + 1) = -trigs
245 85622 : 40 CONTINUE
246 : ELSE
247 11415 : nh = (n - 1)/2
248 11415 : trig(1, 1) = 1._dp
249 11415 : trig(2, 1) = 0._dp
250 113268 : DO 20, i = 1, nh
251 101853 : trigc = COS(i*angle)
252 101853 : trigs = SIN(i*angle)
253 101853 : trig(1, i + 1) = trigc
254 101853 : trig(2, i + 1) = trigs
255 101853 : trig(1, n - i + 1) = trigc
256 101853 : trig(2, n - i + 1) = -trigs
257 11415 : 20 CONTINUE
258 : END IF
259 :
260 97037 : 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 22651284 : 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 22651284 : atn = after*now
300 22651284 : atb = after*before
301 :
302 : ! sqrt(.5_dp)
303 22651284 : rt2i = 0.7071067811865475_dp
304 : IF (now .EQ. 2) THEN
305 0 : ia = 1
306 0 : nin1 = ia - after
307 0 : nout1 = ia - atn
308 0 : DO ib = 1, before
309 0 : nin1 = nin1 + after
310 0 : nin2 = nin1 + atb
311 0 : nout1 = nout1 + atn
312 0 : nout2 = nout1 + after
313 0 : DO j = 1, nfft
314 0 : r1 = zin(1, j, nin1)
315 0 : s1 = zin(2, j, nin1)
316 0 : r2 = zin(1, j, nin2)
317 0 : s2 = zin(2, j, nin2)
318 0 : zout(1, j, nout1) = r2 + r1
319 0 : zout(2, j, nout1) = s2 + s1
320 0 : zout(1, j, nout2) = r1 - r2
321 0 : zout(2, j, nout2) = s1 - s2
322 : END DO; END DO
323 0 : DO 2000, ia = 2, after
324 0 : ias = ia - 1
325 0 : IF (2*ias .EQ. after) THEN
326 0 : IF (isign .EQ. 1) THEN
327 0 : nin1 = ia - after
328 0 : nout1 = ia - atn
329 0 : DO ib = 1, before
330 0 : nin1 = nin1 + after
331 0 : nin2 = nin1 + atb
332 0 : nout1 = nout1 + atn
333 0 : nout2 = nout1 + after
334 0 : DO j = 1, nfft
335 0 : r1 = zin(1, j, nin1)
336 0 : s1 = zin(2, j, nin1)
337 0 : r2 = zin(2, j, nin2)
338 0 : s2 = zin(1, j, nin2)
339 0 : zout(1, j, nout1) = r1 - r2
340 0 : zout(2, j, nout1) = s2 + s1
341 0 : zout(1, j, nout2) = r2 + r1
342 0 : zout(2, j, nout2) = s1 - s2
343 : END DO; END DO
344 : ELSE
345 0 : nin1 = ia - after
346 0 : nout1 = ia - atn
347 0 : DO ib = 1, before
348 0 : nin1 = nin1 + after
349 0 : nin2 = nin1 + atb
350 0 : nout1 = nout1 + atn
351 0 : nout2 = nout1 + after
352 0 : DO j = 1, nfft
353 0 : r1 = zin(1, j, nin1)
354 0 : s1 = zin(2, j, nin1)
355 0 : r2 = zin(2, j, nin2)
356 0 : s2 = zin(1, j, nin2)
357 0 : zout(1, j, nout1) = r2 + r1
358 0 : zout(2, j, nout1) = s1 - s2
359 0 : zout(1, j, nout2) = r1 - r2
360 0 : zout(2, j, nout2) = s2 + s1
361 : END DO; END DO
362 : END IF
363 0 : ELSE IF (4*ias .EQ. after) THEN
364 0 : IF (isign .EQ. 1) THEN
365 0 : nin1 = ia - after
366 0 : nout1 = ia - atn
367 0 : DO ib = 1, before
368 0 : nin1 = nin1 + after
369 0 : nin2 = nin1 + atb
370 0 : nout1 = nout1 + atn
371 0 : nout2 = nout1 + after
372 0 : DO j = 1, nfft
373 0 : r1 = zin(1, j, nin1)
374 0 : s1 = zin(2, j, nin1)
375 0 : r = zin(1, j, nin2)
376 0 : s = zin(2, j, nin2)
377 0 : r2 = (r - s)*rt2i
378 0 : s2 = (r + s)*rt2i
379 0 : zout(1, j, nout1) = r2 + r1
380 0 : zout(2, j, nout1) = s2 + s1
381 0 : zout(1, j, nout2) = r1 - r2
382 0 : zout(2, j, nout2) = s1 - s2
383 : END DO; END DO
384 : ELSE
385 0 : nin1 = ia - after
386 0 : nout1 = ia - atn
387 0 : DO ib = 1, before
388 0 : nin1 = nin1 + after
389 0 : nin2 = nin1 + atb
390 0 : nout1 = nout1 + atn
391 0 : nout2 = nout1 + after
392 0 : DO j = 1, nfft
393 0 : r1 = zin(1, j, nin1)
394 0 : s1 = zin(2, j, nin1)
395 0 : r = zin(1, j, nin2)
396 0 : s = zin(2, j, nin2)
397 0 : r2 = (r + s)*rt2i
398 0 : s2 = (s - r)*rt2i
399 0 : zout(1, j, nout1) = r2 + r1
400 0 : zout(2, j, nout1) = s2 + s1
401 0 : zout(1, j, nout2) = r1 - r2
402 0 : zout(2, j, nout2) = s1 - s2
403 : END DO; END DO
404 : END IF
405 0 : ELSE IF (4*ias .EQ. 3*after) THEN
406 0 : IF (isign .EQ. 1) THEN
407 0 : nin1 = ia - after
408 0 : nout1 = ia - atn
409 0 : DO ib = 1, before
410 0 : nin1 = nin1 + after
411 0 : nin2 = nin1 + atb
412 0 : nout1 = nout1 + atn
413 0 : nout2 = nout1 + after
414 0 : DO j = 1, nfft
415 0 : r1 = zin(1, j, nin1)
416 0 : s1 = zin(2, j, nin1)
417 0 : r = zin(1, j, nin2)
418 0 : s = zin(2, j, nin2)
419 0 : r2 = (r + s)*rt2i
420 0 : s2 = (r - s)*rt2i
421 0 : zout(1, j, nout1) = r1 - r2
422 0 : zout(2, j, nout1) = s2 + s1
423 0 : zout(1, j, nout2) = r2 + r1
424 0 : zout(2, j, nout2) = s1 - s2
425 : END DO; END DO
426 : ELSE
427 0 : nin1 = ia - after
428 0 : nout1 = ia - atn
429 0 : DO ib = 1, before
430 0 : nin1 = nin1 + after
431 0 : nin2 = nin1 + atb
432 0 : nout1 = nout1 + atn
433 0 : nout2 = nout1 + after
434 0 : DO j = 1, nfft
435 0 : r1 = zin(1, j, nin1)
436 0 : s1 = zin(2, j, nin1)
437 0 : r = zin(1, j, nin2)
438 0 : s = zin(2, j, nin2)
439 0 : r2 = (s - r)*rt2i
440 0 : s2 = (r + s)*rt2i
441 0 : zout(1, j, nout1) = r2 + r1
442 0 : zout(2, j, nout1) = s1 - s2
443 0 : zout(1, j, nout2) = r1 - r2
444 0 : zout(2, j, nout2) = s2 + s1
445 : END DO; END DO
446 : END IF
447 : ELSE
448 0 : itrig = ias*before + 1
449 0 : cr2 = trig(1, itrig)
450 0 : ci2 = trig(2, itrig)
451 0 : nin1 = ia - after
452 0 : nout1 = ia - atn
453 0 : DO ib = 1, before
454 0 : nin1 = nin1 + after
455 0 : nin2 = nin1 + atb
456 0 : nout1 = nout1 + atn
457 0 : nout2 = nout1 + after
458 0 : DO j = 1, nfft
459 0 : r1 = zin(1, j, nin1)
460 0 : s1 = zin(2, j, nin1)
461 0 : r = zin(1, j, nin2)
462 0 : s = zin(2, j, nin2)
463 0 : r2 = r*cr2 - s*ci2
464 0 : s2 = r*ci2 + s*cr2
465 0 : zout(1, j, nout1) = r2 + r1
466 0 : zout(2, j, nout1) = s2 + s1
467 0 : zout(1, j, nout2) = r1 - r2
468 0 : zout(2, j, nout2) = s1 - s2
469 : END DO; END DO
470 : END IF
471 0 : 2000 CONTINUE
472 : ELSE IF (now .EQ. 4) THEN
473 5444247 : IF (isign .EQ. 1) THEN
474 2789345 : ia = 1
475 2789345 : nin1 = ia - after
476 2789345 : nout1 = ia - atn
477 24673614 : DO ib = 1, before
478 21884269 : nin1 = nin1 + after
479 21884269 : nin2 = nin1 + atb
480 21884269 : nin3 = nin2 + atb
481 21884269 : nin4 = nin3 + atb
482 21884269 : nout1 = nout1 + atn
483 21884269 : nout2 = nout1 + after
484 21884269 : nout3 = nout2 + after
485 21884269 : nout4 = nout3 + after
486 448864845 : DO j = 1, nfft
487 424191231 : r1 = zin(1, j, nin1)
488 424191231 : s1 = zin(2, j, nin1)
489 424191231 : r2 = zin(1, j, nin2)
490 424191231 : s2 = zin(2, j, nin2)
491 424191231 : r3 = zin(1, j, nin3)
492 424191231 : s3 = zin(2, j, nin3)
493 424191231 : r4 = zin(1, j, nin4)
494 424191231 : s4 = zin(2, j, nin4)
495 424191231 : r = r1 + r3
496 424191231 : s = r2 + r4
497 424191231 : zout(1, j, nout1) = r + s
498 424191231 : zout(1, j, nout3) = r - s
499 424191231 : r = r1 - r3
500 424191231 : s = s2 - s4
501 424191231 : zout(1, j, nout2) = r - s
502 424191231 : zout(1, j, nout4) = r + s
503 424191231 : r = s1 + s3
504 424191231 : s = s2 + s4
505 424191231 : zout(2, j, nout1) = r + s
506 424191231 : zout(2, j, nout3) = r - s
507 424191231 : r = s1 - s3
508 424191231 : s = r2 - r4
509 424191231 : zout(2, j, nout2) = r + s
510 446075500 : zout(2, j, nout4) = r - s
511 : END DO; END DO
512 28197876 : DO 4000, ia = 2, after
513 25408531 : ias = ia - 1
514 25408531 : IF (2*ias .EQ. after) THEN
515 1159813 : nin1 = ia - after
516 1159813 : nout1 = ia - atn
517 2730344 : DO ib = 1, before
518 1570531 : nin1 = nin1 + after
519 1570531 : nin2 = nin1 + atb
520 1570531 : nin3 = nin2 + atb
521 1570531 : nin4 = nin3 + atb
522 1570531 : nout1 = nout1 + atn
523 1570531 : nout2 = nout1 + after
524 1570531 : nout3 = nout2 + after
525 1570531 : nout4 = nout3 + after
526 32911685 : DO j = 1, nfft
527 30181341 : r1 = zin(1, j, nin1)
528 30181341 : s1 = zin(2, j, nin1)
529 30181341 : r = zin(1, j, nin2)
530 30181341 : s = zin(2, j, nin2)
531 30181341 : r2 = (r - s)*rt2i
532 30181341 : s2 = (r + s)*rt2i
533 30181341 : r3 = zin(2, j, nin3)
534 30181341 : s3 = zin(1, j, nin3)
535 30181341 : r = zin(1, j, nin4)
536 30181341 : s = zin(2, j, nin4)
537 30181341 : r4 = (r + s)*rt2i
538 30181341 : s4 = (r - s)*rt2i
539 30181341 : r = r1 - r3
540 30181341 : s = r2 - r4
541 30181341 : zout(1, j, nout1) = r + s
542 30181341 : zout(1, j, nout3) = r - s
543 30181341 : r = r1 + r3
544 30181341 : s = s2 - s4
545 30181341 : zout(1, j, nout2) = r - s
546 30181341 : zout(1, j, nout4) = r + s
547 30181341 : r = s1 + s3
548 30181341 : s = s2 + s4
549 30181341 : zout(2, j, nout1) = r + s
550 30181341 : zout(2, j, nout3) = r - s
551 30181341 : r = s1 - s3
552 30181341 : s = r2 + r4
553 30181341 : zout(2, j, nout2) = r + s
554 31751872 : zout(2, j, nout4) = r - s
555 : END DO; END DO
556 : ELSE
557 24248718 : itt = ias*before
558 24248718 : itrig = itt + 1
559 24248718 : cr2 = trig(1, itrig)
560 24248718 : ci2 = trig(2, itrig)
561 24248718 : itrig = itrig + itt
562 24248718 : cr3 = trig(1, itrig)
563 24248718 : ci3 = trig(2, itrig)
564 24248718 : itrig = itrig + itt
565 24248718 : cr4 = trig(1, itrig)
566 24248718 : ci4 = trig(2, itrig)
567 24248718 : nin1 = ia - after
568 24248718 : nout1 = ia - atn
569 60549024 : DO ib = 1, before
570 36300306 : nin1 = nin1 + after
571 36300306 : nin2 = nin1 + atb
572 36300306 : nin3 = nin2 + atb
573 36300306 : nin4 = nin3 + atb
574 36300306 : nout1 = nout1 + atn
575 36300306 : nout2 = nout1 + after
576 36300306 : nout3 = nout2 + after
577 36300306 : nout4 = nout3 + after
578 746324478 : DO j = 1, nfft
579 685775454 : r1 = zin(1, j, nin1)
580 685775454 : s1 = zin(2, j, nin1)
581 685775454 : r = zin(1, j, nin2)
582 685775454 : s = zin(2, j, nin2)
583 685775454 : r2 = r*cr2 - s*ci2
584 685775454 : s2 = r*ci2 + s*cr2
585 685775454 : r = zin(1, j, nin3)
586 685775454 : s = zin(2, j, nin3)
587 685775454 : r3 = r*cr3 - s*ci3
588 685775454 : s3 = r*ci3 + s*cr3
589 685775454 : r = zin(1, j, nin4)
590 685775454 : s = zin(2, j, nin4)
591 685775454 : r4 = r*cr4 - s*ci4
592 685775454 : s4 = r*ci4 + s*cr4
593 685775454 : r = r1 + r3
594 685775454 : s = r2 + r4
595 685775454 : zout(1, j, nout1) = r + s
596 685775454 : zout(1, j, nout3) = r - s
597 685775454 : r = r1 - r3
598 685775454 : s = s2 - s4
599 685775454 : zout(1, j, nout2) = r - s
600 685775454 : zout(1, j, nout4) = r + s
601 685775454 : r = s1 + s3
602 685775454 : s = s2 + s4
603 685775454 : zout(2, j, nout1) = r + s
604 685775454 : zout(2, j, nout3) = r - s
605 685775454 : r = s1 - s3
606 685775454 : s = r2 - r4
607 685775454 : zout(2, j, nout2) = r + s
608 722075760 : zout(2, j, nout4) = r - s
609 : END DO; END DO
610 : END IF
611 2789345 : 4000 CONTINUE
612 : ELSE
613 2654902 : ia = 1
614 2654902 : nin1 = ia - after
615 2654902 : nout1 = ia - atn
616 23408701 : DO ib = 1, before
617 20753799 : nin1 = nin1 + after
618 20753799 : nin2 = nin1 + atb
619 20753799 : nin3 = nin2 + atb
620 20753799 : nin4 = nin3 + atb
621 20753799 : nout1 = nout1 + atn
622 20753799 : nout2 = nout1 + after
623 20753799 : nout3 = nout2 + after
624 20753799 : nout4 = nout3 + after
625 426843633 : DO j = 1, nfft
626 403434932 : r1 = zin(1, j, nin1)
627 403434932 : s1 = zin(2, j, nin1)
628 403434932 : r2 = zin(1, j, nin2)
629 403434932 : s2 = zin(2, j, nin2)
630 403434932 : r3 = zin(1, j, nin3)
631 403434932 : s3 = zin(2, j, nin3)
632 403434932 : r4 = zin(1, j, nin4)
633 403434932 : s4 = zin(2, j, nin4)
634 403434932 : r = r1 + r3
635 403434932 : s = r2 + r4
636 403434932 : zout(1, j, nout1) = r + s
637 403434932 : zout(1, j, nout3) = r - s
638 403434932 : r = r1 - r3
639 403434932 : s = s2 - s4
640 403434932 : zout(1, j, nout2) = r + s
641 403434932 : zout(1, j, nout4) = r - s
642 403434932 : r = s1 + s3
643 403434932 : s = s2 + s4
644 403434932 : zout(2, j, nout1) = r + s
645 403434932 : zout(2, j, nout3) = r - s
646 403434932 : r = s1 - s3
647 403434932 : s = r2 - r4
648 403434932 : zout(2, j, nout2) = r - s
649 424188731 : zout(2, j, nout4) = r + s
650 : END DO; END DO
651 26582875 : DO 4100, ia = 2, after
652 23927973 : ias = ia - 1
653 23927973 : IF (2*ias .EQ. after) THEN
654 1098903 : nin1 = ia - after
655 1098903 : nout1 = ia - atn
656 2571156 : DO ib = 1, before
657 1472253 : nin1 = nin1 + after
658 1472253 : nin2 = nin1 + atb
659 1472253 : nin3 = nin2 + atb
660 1472253 : nin4 = nin3 + atb
661 1472253 : nout1 = nout1 + atn
662 1472253 : nout2 = nout1 + after
663 1472253 : nout3 = nout2 + after
664 1472253 : nout4 = nout3 + after
665 30889868 : DO j = 1, nfft
666 28318712 : r1 = zin(1, j, nin1)
667 28318712 : s1 = zin(2, j, nin1)
668 28318712 : r = zin(1, j, nin2)
669 28318712 : s = zin(2, j, nin2)
670 28318712 : r2 = (r + s)*rt2i
671 28318712 : s2 = (s - r)*rt2i
672 28318712 : r3 = zin(2, j, nin3)
673 28318712 : s3 = zin(1, j, nin3)
674 28318712 : r = zin(1, j, nin4)
675 28318712 : s = zin(2, j, nin4)
676 28318712 : r4 = (s - r)*rt2i
677 28318712 : s4 = (r + s)*rt2i
678 28318712 : r = r1 + r3
679 28318712 : s = r2 + r4
680 28318712 : zout(1, j, nout1) = r + s
681 28318712 : zout(1, j, nout3) = r - s
682 28318712 : r = r1 - r3
683 28318712 : s = s2 + s4
684 28318712 : zout(1, j, nout2) = r + s
685 28318712 : zout(1, j, nout4) = r - s
686 28318712 : r = s1 - s3
687 28318712 : s = s2 - s4
688 28318712 : zout(2, j, nout1) = r + s
689 28318712 : zout(2, j, nout3) = r - s
690 28318712 : r = s1 + s3
691 28318712 : s = r2 - r4
692 28318712 : zout(2, j, nout2) = r - s
693 29790965 : zout(2, j, nout4) = r + s
694 : END DO; END DO
695 : ELSE
696 22829070 : itt = ias*before
697 22829070 : itrig = itt + 1
698 22829070 : cr2 = trig(1, itrig)
699 22829070 : ci2 = trig(2, itrig)
700 22829070 : itrig = itrig + itt
701 22829070 : cr3 = trig(1, itrig)
702 22829070 : ci3 = trig(2, itrig)
703 22829070 : itrig = itrig + itt
704 22829070 : cr4 = trig(1, itrig)
705 22829070 : ci4 = trig(2, itrig)
706 22829070 : nin1 = ia - after
707 22829070 : nout1 = ia - atn
708 57127344 : DO ib = 1, before
709 34298274 : nin1 = nin1 + after
710 34298274 : nin2 = nin1 + atb
711 34298274 : nin3 = nin2 + atb
712 34298274 : nin4 = nin3 + atb
713 34298274 : nout1 = nout1 + atn
714 34298274 : nout2 = nout1 + after
715 34298274 : nout3 = nout2 + after
716 34298274 : nout4 = nout3 + after
717 708740256 : DO j = 1, nfft
718 651612912 : r1 = zin(1, j, nin1)
719 651612912 : s1 = zin(2, j, nin1)
720 651612912 : r = zin(1, j, nin2)
721 651612912 : s = zin(2, j, nin2)
722 651612912 : r2 = r*cr2 - s*ci2
723 651612912 : s2 = r*ci2 + s*cr2
724 651612912 : r = zin(1, j, nin3)
725 651612912 : s = zin(2, j, nin3)
726 651612912 : r3 = r*cr3 - s*ci3
727 651612912 : s3 = r*ci3 + s*cr3
728 651612912 : r = zin(1, j, nin4)
729 651612912 : s = zin(2, j, nin4)
730 651612912 : r4 = r*cr4 - s*ci4
731 651612912 : s4 = r*ci4 + s*cr4
732 651612912 : r = r1 + r3
733 651612912 : s = r2 + r4
734 651612912 : zout(1, j, nout1) = r + s
735 651612912 : zout(1, j, nout3) = r - s
736 651612912 : r = r1 - r3
737 651612912 : s = s2 - s4
738 651612912 : zout(1, j, nout2) = r + s
739 651612912 : zout(1, j, nout4) = r - s
740 651612912 : r = s1 + s3
741 651612912 : s = s2 + s4
742 651612912 : zout(2, j, nout1) = r + s
743 651612912 : zout(2, j, nout3) = r - s
744 651612912 : r = s1 - s3
745 651612912 : s = r2 - r4
746 651612912 : zout(2, j, nout2) = r - s
747 685911186 : zout(2, j, nout4) = r + s
748 : END DO; END DO
749 : END IF
750 2654902 : 4100 CONTINUE
751 : END IF
752 : ELSE IF (now .EQ. 8) THEN
753 2350230 : IF (isign .EQ. -1) THEN
754 1131882 : ia = 1
755 1131882 : nin1 = ia - after
756 1131882 : nout1 = ia - atn
757 9920474 : DO ib = 1, before
758 8788592 : nin1 = nin1 + after
759 8788592 : nin2 = nin1 + atb
760 8788592 : nin3 = nin2 + atb
761 8788592 : nin4 = nin3 + atb
762 8788592 : nin5 = nin4 + atb
763 8788592 : nin6 = nin5 + atb
764 8788592 : nin7 = nin6 + atb
765 8788592 : nin8 = nin7 + atb
766 8788592 : nout1 = nout1 + atn
767 8788592 : nout2 = nout1 + after
768 8788592 : nout3 = nout2 + after
769 8788592 : nout4 = nout3 + after
770 8788592 : nout5 = nout4 + after
771 8788592 : nout6 = nout5 + after
772 8788592 : nout7 = nout6 + after
773 8788592 : nout8 = nout7 + after
774 179285838 : DO j = 1, nfft
775 169365364 : r1 = zin(1, j, nin1)
776 169365364 : s1 = zin(2, j, nin1)
777 169365364 : r2 = zin(1, j, nin2)
778 169365364 : s2 = zin(2, j, nin2)
779 169365364 : r3 = zin(1, j, nin3)
780 169365364 : s3 = zin(2, j, nin3)
781 169365364 : r4 = zin(1, j, nin4)
782 169365364 : s4 = zin(2, j, nin4)
783 169365364 : r5 = zin(1, j, nin5)
784 169365364 : s5 = zin(2, j, nin5)
785 169365364 : r6 = zin(1, j, nin6)
786 169365364 : s6 = zin(2, j, nin6)
787 169365364 : r7 = zin(1, j, nin7)
788 169365364 : s7 = zin(2, j, nin7)
789 169365364 : r8 = zin(1, j, nin8)
790 169365364 : s8 = zin(2, j, nin8)
791 169365364 : r = r1 + r5
792 169365364 : s = r3 + r7
793 169365364 : ap = r + s
794 169365364 : am = r - s
795 169365364 : r = r2 + r6
796 169365364 : s = r4 + r8
797 169365364 : bp = r + s
798 169365364 : bm = r - s
799 169365364 : r = s1 + s5
800 169365364 : s = s3 + s7
801 169365364 : cp = r + s
802 169365364 : cm = r - s
803 169365364 : r = s2 + s6
804 169365364 : s = s4 + s8
805 169365364 : dpp = r + s
806 169365364 : dm = r - s
807 169365364 : zout(1, j, nout1) = ap + bp
808 169365364 : zout(2, j, nout1) = cp + dpp
809 169365364 : zout(1, j, nout5) = ap - bp
810 169365364 : zout(2, j, nout5) = cp - dpp
811 169365364 : zout(1, j, nout3) = am + dm
812 169365364 : zout(2, j, nout3) = cm - bm
813 169365364 : zout(1, j, nout7) = am - dm
814 169365364 : zout(2, j, nout7) = cm + bm
815 169365364 : r = r1 - r5
816 169365364 : s = s3 - s7
817 169365364 : ap = r + s
818 169365364 : am = r - s
819 169365364 : r = s1 - s5
820 169365364 : s = r3 - r7
821 169365364 : bp = r + s
822 169365364 : bm = r - s
823 169365364 : r = s4 - s8
824 169365364 : s = r2 - r6
825 169365364 : cp = r + s
826 169365364 : cm = r - s
827 169365364 : r = s2 - s6
828 169365364 : s = r4 - r8
829 169365364 : dpp = r + s
830 169365364 : dm = r - s
831 169365364 : r = (cp + dm)*rt2i
832 169365364 : s = (dm - cp)*rt2i
833 169365364 : cp = (cm + dpp)*rt2i
834 169365364 : dpp = (cm - dpp)*rt2i
835 169365364 : zout(1, j, nout2) = ap + r
836 169365364 : zout(2, j, nout2) = bm + s
837 169365364 : zout(1, j, nout6) = ap - r
838 169365364 : zout(2, j, nout6) = bm - s
839 169365364 : zout(1, j, nout4) = am + cp
840 169365364 : zout(2, j, nout4) = bp + dpp
841 169365364 : zout(1, j, nout8) = am - cp
842 178153956 : zout(2, j, nout8) = bp - dpp
843 : END DO; END DO
844 3072953 : DO 8000, ia = 2, after
845 1941071 : ias = ia - 1
846 1941071 : itt = ias*before
847 1941071 : itrig = itt + 1
848 1941071 : cr2 = trig(1, itrig)
849 1941071 : ci2 = trig(2, itrig)
850 1941071 : itrig = itrig + itt
851 1941071 : cr3 = trig(1, itrig)
852 1941071 : ci3 = trig(2, itrig)
853 1941071 : itrig = itrig + itt
854 1941071 : cr4 = trig(1, itrig)
855 1941071 : ci4 = trig(2, itrig)
856 1941071 : itrig = itrig + itt
857 1941071 : cr5 = trig(1, itrig)
858 1941071 : ci5 = trig(2, itrig)
859 1941071 : itrig = itrig + itt
860 1941071 : cr6 = trig(1, itrig)
861 1941071 : ci6 = trig(2, itrig)
862 1941071 : itrig = itrig + itt
863 1941071 : cr7 = trig(1, itrig)
864 1941071 : ci7 = trig(2, itrig)
865 1941071 : itrig = itrig + itt
866 1941071 : cr8 = trig(1, itrig)
867 1941071 : ci8 = trig(2, itrig)
868 1941071 : nin1 = ia - after
869 1941071 : nout1 = ia - atn
870 7076352 : DO ib = 1, before
871 5135281 : nin1 = nin1 + after
872 5135281 : nin2 = nin1 + atb
873 5135281 : nin3 = nin2 + atb
874 5135281 : nin4 = nin3 + atb
875 5135281 : nin5 = nin4 + atb
876 5135281 : nin6 = nin5 + atb
877 5135281 : nin7 = nin6 + atb
878 5135281 : nin8 = nin7 + atb
879 5135281 : nout1 = nout1 + atn
880 5135281 : nout2 = nout1 + after
881 5135281 : nout3 = nout2 + after
882 5135281 : nout4 = nout3 + after
883 5135281 : nout5 = nout4 + after
884 5135281 : nout6 = nout5 + after
885 5135281 : nout7 = nout6 + after
886 5135281 : nout8 = nout7 + after
887 76972424 : DO j = 1, nfft
888 69896072 : r1 = zin(1, j, nin1)
889 69896072 : s1 = zin(2, j, nin1)
890 69896072 : r = zin(1, j, nin2)
891 69896072 : s = zin(2, j, nin2)
892 69896072 : r2 = r*cr2 - s*ci2
893 69896072 : s2 = r*ci2 + s*cr2
894 69896072 : r = zin(1, j, nin3)
895 69896072 : s = zin(2, j, nin3)
896 69896072 : r3 = r*cr3 - s*ci3
897 69896072 : s3 = r*ci3 + s*cr3
898 69896072 : r = zin(1, j, nin4)
899 69896072 : s = zin(2, j, nin4)
900 69896072 : r4 = r*cr4 - s*ci4
901 69896072 : s4 = r*ci4 + s*cr4
902 69896072 : r = zin(1, j, nin5)
903 69896072 : s = zin(2, j, nin5)
904 69896072 : r5 = r*cr5 - s*ci5
905 69896072 : s5 = r*ci5 + s*cr5
906 69896072 : r = zin(1, j, nin6)
907 69896072 : s = zin(2, j, nin6)
908 69896072 : r6 = r*cr6 - s*ci6
909 69896072 : s6 = r*ci6 + s*cr6
910 69896072 : r = zin(1, j, nin7)
911 69896072 : s = zin(2, j, nin7)
912 69896072 : r7 = r*cr7 - s*ci7
913 69896072 : s7 = r*ci7 + s*cr7
914 69896072 : r = zin(1, j, nin8)
915 69896072 : s = zin(2, j, nin8)
916 69896072 : r8 = r*cr8 - s*ci8
917 69896072 : s8 = r*ci8 + s*cr8
918 69896072 : r = r1 + r5
919 69896072 : s = r3 + r7
920 69896072 : ap = r + s
921 69896072 : am = r - s
922 69896072 : r = r2 + r6
923 69896072 : s = r4 + r8
924 69896072 : bp = r + s
925 69896072 : bm = r - s
926 69896072 : r = s1 + s5
927 69896072 : s = s3 + s7
928 69896072 : cp = r + s
929 69896072 : cm = r - s
930 69896072 : r = s2 + s6
931 69896072 : s = s4 + s8
932 69896072 : dpp = r + s
933 69896072 : dm = r - s
934 69896072 : zout(1, j, nout1) = ap + bp
935 69896072 : zout(2, j, nout1) = cp + dpp
936 69896072 : zout(1, j, nout5) = ap - bp
937 69896072 : zout(2, j, nout5) = cp - dpp
938 69896072 : zout(1, j, nout3) = am + dm
939 69896072 : zout(2, j, nout3) = cm - bm
940 69896072 : zout(1, j, nout7) = am - dm
941 69896072 : zout(2, j, nout7) = cm + bm
942 69896072 : r = r1 - r5
943 69896072 : s = s3 - s7
944 69896072 : ap = r + s
945 69896072 : am = r - s
946 69896072 : r = s1 - s5
947 69896072 : s = r3 - r7
948 69896072 : bp = r + s
949 69896072 : bm = r - s
950 69896072 : r = s4 - s8
951 69896072 : s = r2 - r6
952 69896072 : cp = r + s
953 69896072 : cm = r - s
954 69896072 : r = s2 - s6
955 69896072 : s = r4 - r8
956 69896072 : dpp = r + s
957 69896072 : dm = r - s
958 69896072 : r = (cp + dm)*rt2i
959 69896072 : s = (dm - cp)*rt2i
960 69896072 : cp = (cm + dpp)*rt2i
961 69896072 : dpp = (cm - dpp)*rt2i
962 69896072 : zout(1, j, nout2) = ap + r
963 69896072 : zout(2, j, nout2) = bm + s
964 69896072 : zout(1, j, nout6) = ap - r
965 69896072 : zout(2, j, nout6) = bm - s
966 69896072 : zout(1, j, nout4) = am + cp
967 69896072 : zout(2, j, nout4) = bp + dpp
968 69896072 : zout(1, j, nout8) = am - cp
969 75031353 : zout(2, j, nout8) = bp - dpp
970 : END DO; END DO
971 1131882 : 8000 CONTINUE
972 :
973 : ELSE
974 1218348 : ia = 1
975 1218348 : nin1 = ia - after
976 1218348 : nout1 = ia - atn
977 10790006 : DO ib = 1, before
978 9571658 : nin1 = nin1 + after
979 9571658 : nin2 = nin1 + atb
980 9571658 : nin3 = nin2 + atb
981 9571658 : nin4 = nin3 + atb
982 9571658 : nin5 = nin4 + atb
983 9571658 : nin6 = nin5 + atb
984 9571658 : nin7 = nin6 + atb
985 9571658 : nin8 = nin7 + atb
986 9571658 : nout1 = nout1 + atn
987 9571658 : nout2 = nout1 + after
988 9571658 : nout3 = nout2 + after
989 9571658 : nout4 = nout3 + after
990 9571658 : nout5 = nout4 + after
991 9571658 : nout6 = nout5 + after
992 9571658 : nout7 = nout6 + after
993 9571658 : nout8 = nout7 + after
994 194474622 : DO j = 1, nfft
995 183684616 : r1 = zin(1, j, nin1)
996 183684616 : s1 = zin(2, j, nin1)
997 183684616 : r2 = zin(1, j, nin2)
998 183684616 : s2 = zin(2, j, nin2)
999 183684616 : r3 = zin(1, j, nin3)
1000 183684616 : s3 = zin(2, j, nin3)
1001 183684616 : r4 = zin(1, j, nin4)
1002 183684616 : s4 = zin(2, j, nin4)
1003 183684616 : r5 = zin(1, j, nin5)
1004 183684616 : s5 = zin(2, j, nin5)
1005 183684616 : r6 = zin(1, j, nin6)
1006 183684616 : s6 = zin(2, j, nin6)
1007 183684616 : r7 = zin(1, j, nin7)
1008 183684616 : s7 = zin(2, j, nin7)
1009 183684616 : r8 = zin(1, j, nin8)
1010 183684616 : s8 = zin(2, j, nin8)
1011 183684616 : r = r1 + r5
1012 183684616 : s = r3 + r7
1013 183684616 : ap = r + s
1014 183684616 : am = r - s
1015 183684616 : r = r2 + r6
1016 183684616 : s = r4 + r8
1017 183684616 : bp = r + s
1018 183684616 : bm = r - s
1019 183684616 : r = s1 + s5
1020 183684616 : s = s3 + s7
1021 183684616 : cp = r + s
1022 183684616 : cm = r - s
1023 183684616 : r = s2 + s6
1024 183684616 : s = s4 + s8
1025 183684616 : dpp = r + s
1026 183684616 : dm = r - s
1027 183684616 : zout(1, j, nout1) = ap + bp
1028 183684616 : zout(2, j, nout1) = cp + dpp
1029 183684616 : zout(1, j, nout5) = ap - bp
1030 183684616 : zout(2, j, nout5) = cp - dpp
1031 183684616 : zout(1, j, nout3) = am - dm
1032 183684616 : zout(2, j, nout3) = cm + bm
1033 183684616 : zout(1, j, nout7) = am + dm
1034 183684616 : zout(2, j, nout7) = cm - bm
1035 183684616 : r = r1 - r5
1036 183684616 : s = -s3 + s7
1037 183684616 : ap = r + s
1038 183684616 : am = r - s
1039 183684616 : r = s1 - s5
1040 183684616 : s = r7 - r3
1041 183684616 : bp = r + s
1042 183684616 : bm = r - s
1043 183684616 : r = -s4 + s8
1044 183684616 : s = r2 - r6
1045 183684616 : cp = r + s
1046 183684616 : cm = r - s
1047 183684616 : r = -s2 + s6
1048 183684616 : s = r4 - r8
1049 183684616 : dpp = r + s
1050 183684616 : dm = r - s
1051 183684616 : r = (cp + dm)*rt2i
1052 183684616 : s = (cp - dm)*rt2i
1053 183684616 : cp = (cm + dpp)*rt2i
1054 183684616 : dpp = (dpp - cm)*rt2i
1055 183684616 : zout(1, j, nout2) = ap + r
1056 183684616 : zout(2, j, nout2) = bm + s
1057 183684616 : zout(1, j, nout6) = ap - r
1058 183684616 : zout(2, j, nout6) = bm - s
1059 183684616 : zout(1, j, nout4) = am + cp
1060 183684616 : zout(2, j, nout4) = bp + dpp
1061 183684616 : zout(1, j, nout8) = am - cp
1062 193256274 : zout(2, j, nout8) = bp - dpp
1063 : END DO; END DO
1064 :
1065 3309794 : DO 8001, ia = 2, after
1066 2091446 : ias = ia - 1
1067 2091446 : itt = ias*before
1068 2091446 : itrig = itt + 1
1069 2091446 : cr2 = trig(1, itrig)
1070 2091446 : ci2 = trig(2, itrig)
1071 2091446 : itrig = itrig + itt
1072 2091446 : cr3 = trig(1, itrig)
1073 2091446 : ci3 = trig(2, itrig)
1074 2091446 : itrig = itrig + itt
1075 2091446 : cr4 = trig(1, itrig)
1076 2091446 : ci4 = trig(2, itrig)
1077 2091446 : itrig = itrig + itt
1078 2091446 : cr5 = trig(1, itrig)
1079 2091446 : ci5 = trig(2, itrig)
1080 2091446 : itrig = itrig + itt
1081 2091446 : cr6 = trig(1, itrig)
1082 2091446 : ci6 = trig(2, itrig)
1083 2091446 : itrig = itrig + itt
1084 2091446 : cr7 = trig(1, itrig)
1085 2091446 : ci7 = trig(2, itrig)
1086 2091446 : itrig = itrig + itt
1087 2091446 : cr8 = trig(1, itrig)
1088 2091446 : ci8 = trig(2, itrig)
1089 2091446 : nin1 = ia - after
1090 2091446 : nout1 = ia - atn
1091 7658763 : DO ib = 1, before
1092 5567317 : nin1 = nin1 + after
1093 5567317 : nin2 = nin1 + atb
1094 5567317 : nin3 = nin2 + atb
1095 5567317 : nin4 = nin3 + atb
1096 5567317 : nin5 = nin4 + atb
1097 5567317 : nin6 = nin5 + atb
1098 5567317 : nin7 = nin6 + atb
1099 5567317 : nin8 = nin7 + atb
1100 5567317 : nout1 = nout1 + atn
1101 5567317 : nout2 = nout1 + after
1102 5567317 : nout3 = nout2 + after
1103 5567317 : nout4 = nout3 + after
1104 5567317 : nout5 = nout4 + after
1105 5567317 : nout6 = nout5 + after
1106 5567317 : nout7 = nout6 + after
1107 5567317 : nout8 = nout7 + after
1108 82731440 : DO j = 1, nfft
1109 75072677 : r1 = zin(1, j, nin1)
1110 75072677 : s1 = zin(2, j, nin1)
1111 75072677 : r = zin(1, j, nin2)
1112 75072677 : s = zin(2, j, nin2)
1113 75072677 : r2 = r*cr2 - s*ci2
1114 75072677 : s2 = r*ci2 + s*cr2
1115 75072677 : r = zin(1, j, nin3)
1116 75072677 : s = zin(2, j, nin3)
1117 75072677 : r3 = r*cr3 - s*ci3
1118 75072677 : s3 = r*ci3 + s*cr3
1119 75072677 : r = zin(1, j, nin4)
1120 75072677 : s = zin(2, j, nin4)
1121 75072677 : r4 = r*cr4 - s*ci4
1122 75072677 : s4 = r*ci4 + s*cr4
1123 75072677 : r = zin(1, j, nin5)
1124 75072677 : s = zin(2, j, nin5)
1125 75072677 : r5 = r*cr5 - s*ci5
1126 75072677 : s5 = r*ci5 + s*cr5
1127 75072677 : r = zin(1, j, nin6)
1128 75072677 : s = zin(2, j, nin6)
1129 75072677 : r6 = r*cr6 - s*ci6
1130 75072677 : s6 = r*ci6 + s*cr6
1131 75072677 : r = zin(1, j, nin7)
1132 75072677 : s = zin(2, j, nin7)
1133 75072677 : r7 = r*cr7 - s*ci7
1134 75072677 : s7 = r*ci7 + s*cr7
1135 75072677 : r = zin(1, j, nin8)
1136 75072677 : s = zin(2, j, nin8)
1137 75072677 : r8 = r*cr8 - s*ci8
1138 75072677 : s8 = r*ci8 + s*cr8
1139 75072677 : r = r1 + r5
1140 75072677 : s = r3 + r7
1141 75072677 : ap = r + s
1142 75072677 : am = r - s
1143 75072677 : r = r2 + r6
1144 75072677 : s = r4 + r8
1145 75072677 : bp = r + s
1146 75072677 : bm = r - s
1147 75072677 : r = s1 + s5
1148 75072677 : s = s3 + s7
1149 75072677 : cp = r + s
1150 75072677 : cm = r - s
1151 75072677 : r = s2 + s6
1152 75072677 : s = s4 + s8
1153 75072677 : dpp = r + s
1154 75072677 : dm = r - s
1155 75072677 : zout(1, j, nout1) = ap + bp
1156 75072677 : zout(2, j, nout1) = cp + dpp
1157 75072677 : zout(1, j, nout5) = ap - bp
1158 75072677 : zout(2, j, nout5) = cp - dpp
1159 75072677 : zout(1, j, nout3) = am - dm
1160 75072677 : zout(2, j, nout3) = cm + bm
1161 75072677 : zout(1, j, nout7) = am + dm
1162 75072677 : zout(2, j, nout7) = cm - bm
1163 75072677 : r = r1 - r5
1164 75072677 : s = -s3 + s7
1165 75072677 : ap = r + s
1166 75072677 : am = r - s
1167 75072677 : r = s1 - s5
1168 75072677 : s = r7 - r3
1169 75072677 : bp = r + s
1170 75072677 : bm = r - s
1171 75072677 : r = -s4 + s8
1172 75072677 : s = r2 - r6
1173 75072677 : cp = r + s
1174 75072677 : cm = r - s
1175 75072677 : r = -s2 + s6
1176 75072677 : s = r4 - r8
1177 75072677 : dpp = r + s
1178 75072677 : dm = r - s
1179 75072677 : r = (cp + dm)*rt2i
1180 75072677 : s = (cp - dm)*rt2i
1181 75072677 : cp = (cm + dpp)*rt2i
1182 75072677 : dpp = (dpp - cm)*rt2i
1183 75072677 : zout(1, j, nout2) = ap + r
1184 75072677 : zout(2, j, nout2) = bm + s
1185 75072677 : zout(1, j, nout6) = ap - r
1186 75072677 : zout(2, j, nout6) = bm - s
1187 75072677 : zout(1, j, nout4) = am + cp
1188 75072677 : zout(2, j, nout4) = bp + dpp
1189 75072677 : zout(1, j, nout8) = am - cp
1190 80639994 : zout(2, j, nout8) = bp - dpp
1191 : END DO; END DO
1192 1218348 : 8001 CONTINUE
1193 :
1194 : END IF
1195 : ELSE IF (now .EQ. 3) THEN
1196 : ! .5_dp*sqrt(3._dp)
1197 8840676 : bb = isign*0.8660254037844387_dp
1198 8840676 : ia = 1
1199 8840676 : nin1 = ia - after
1200 8840676 : nout1 = ia - atn
1201 32184910 : DO ib = 1, before
1202 23344234 : nin1 = nin1 + after
1203 23344234 : nin2 = nin1 + atb
1204 23344234 : nin3 = nin2 + atb
1205 23344234 : nout1 = nout1 + atn
1206 23344234 : nout2 = nout1 + after
1207 23344234 : nout3 = nout2 + after
1208 486805186 : DO j = 1, nfft
1209 454620276 : r1 = zin(1, j, nin1)
1210 454620276 : s1 = zin(2, j, nin1)
1211 454620276 : r2 = zin(1, j, nin2)
1212 454620276 : s2 = zin(2, j, nin2)
1213 454620276 : r3 = zin(1, j, nin3)
1214 454620276 : s3 = zin(2, j, nin3)
1215 454620276 : r = r2 + r3
1216 454620276 : s = s2 + s3
1217 454620276 : zout(1, j, nout1) = r + r1
1218 454620276 : zout(2, j, nout1) = s + s1
1219 454620276 : r1 = r1 - .5_dp*r
1220 454620276 : s1 = s1 - .5_dp*s
1221 454620276 : r2 = bb*(r2 - r3)
1222 454620276 : s2 = bb*(s2 - s3)
1223 454620276 : zout(1, j, nout2) = r1 - s2
1224 454620276 : zout(2, j, nout2) = s1 + r2
1225 454620276 : zout(1, j, nout3) = r1 + s2
1226 477964510 : zout(2, j, nout3) = s1 - r2
1227 : END DO; END DO
1228 161073778 : DO 3000, ia = 2, after
1229 152233102 : ias = ia - 1
1230 152233102 : IF (4*ias .EQ. 3*after) THEN
1231 5740655 : IF (isign .EQ. 1) THEN
1232 2947927 : nin1 = ia - after
1233 2947927 : nout1 = ia - atn
1234 12049224 : DO ib = 1, before
1235 9101297 : nin1 = nin1 + after
1236 9101297 : nin2 = nin1 + atb
1237 9101297 : nin3 = nin2 + atb
1238 9101297 : nout1 = nout1 + atn
1239 9101297 : nout2 = nout1 + after
1240 9101297 : nout3 = nout2 + after
1241 185997352 : DO j = 1, nfft
1242 173948128 : r1 = zin(1, j, nin1)
1243 173948128 : s1 = zin(2, j, nin1)
1244 173948128 : r2 = zin(2, j, nin2)
1245 173948128 : s2 = zin(1, j, nin2)
1246 173948128 : r3 = zin(1, j, nin3)
1247 173948128 : s3 = zin(2, j, nin3)
1248 173948128 : r = r3 + r2
1249 173948128 : s = s2 - s3
1250 173948128 : zout(1, j, nout1) = r1 - r
1251 173948128 : zout(2, j, nout1) = s + s1
1252 173948128 : r1 = r1 + .5_dp*r
1253 173948128 : s1 = s1 - .5_dp*s
1254 173948128 : r2 = bb*(r2 - r3)
1255 173948128 : s2 = bb*(s2 + s3)
1256 173948128 : zout(1, j, nout2) = r1 - s2
1257 173948128 : zout(2, j, nout2) = s1 - r2
1258 173948128 : zout(1, j, nout3) = r1 + s2
1259 183049425 : zout(2, j, nout3) = s1 + r2
1260 : END DO; END DO
1261 : ELSE
1262 2792728 : nin1 = ia - after
1263 2792728 : nout1 = ia - atn
1264 11406510 : DO ib = 1, before
1265 8613782 : nin1 = nin1 + after
1266 8613782 : nin2 = nin1 + atb
1267 8613782 : nin3 = nin2 + atb
1268 8613782 : nout1 = nout1 + atn
1269 8613782 : nout2 = nout1 + after
1270 8613782 : nout3 = nout2 + after
1271 176682786 : DO j = 1, nfft
1272 165276276 : r1 = zin(1, j, nin1)
1273 165276276 : s1 = zin(2, j, nin1)
1274 165276276 : r2 = zin(2, j, nin2)
1275 165276276 : s2 = zin(1, j, nin2)
1276 165276276 : r3 = zin(1, j, nin3)
1277 165276276 : s3 = zin(2, j, nin3)
1278 165276276 : r = r2 - r3
1279 165276276 : s = s2 + s3
1280 165276276 : zout(1, j, nout1) = r + r1
1281 165276276 : zout(2, j, nout1) = s1 - s
1282 165276276 : r1 = r1 - .5_dp*r
1283 165276276 : s1 = s1 + .5_dp*s
1284 165276276 : r2 = bb*(r2 + r3)
1285 165276276 : s2 = bb*(s2 - s3)
1286 165276276 : zout(1, j, nout2) = r1 + s2
1287 165276276 : zout(2, j, nout2) = s1 + r2
1288 165276276 : zout(1, j, nout3) = r1 - s2
1289 173890058 : zout(2, j, nout3) = s1 - r2
1290 : END DO; END DO
1291 : END IF
1292 146492447 : ELSE IF (8*ias .EQ. 3*after) THEN
1293 1910407 : IF (isign .EQ. 1) THEN
1294 982133 : nin1 = ia - after
1295 982133 : nout1 = ia - atn
1296 2375112 : DO ib = 1, before
1297 1392979 : nin1 = nin1 + after
1298 1392979 : nin2 = nin1 + atb
1299 1392979 : nin3 = nin2 + atb
1300 1392979 : nout1 = nout1 + atn
1301 1392979 : nout2 = nout1 + after
1302 1392979 : nout3 = nout2 + after
1303 29906542 : DO j = 1, nfft
1304 27531430 : r1 = zin(1, j, nin1)
1305 27531430 : s1 = zin(2, j, nin1)
1306 27531430 : r = zin(1, j, nin2)
1307 27531430 : s = zin(2, j, nin2)
1308 27531430 : r2 = (r - s)*rt2i
1309 27531430 : s2 = (r + s)*rt2i
1310 27531430 : r3 = zin(2, j, nin3)
1311 27531430 : s3 = zin(1, j, nin3)
1312 27531430 : r = r2 - r3
1313 27531430 : s = s2 + s3
1314 27531430 : zout(1, j, nout1) = r + r1
1315 27531430 : zout(2, j, nout1) = s + s1
1316 27531430 : r1 = r1 - .5_dp*r
1317 27531430 : s1 = s1 - .5_dp*s
1318 27531430 : r2 = bb*(r2 + r3)
1319 27531430 : s2 = bb*(s2 - s3)
1320 27531430 : zout(1, j, nout2) = r1 - s2
1321 27531430 : zout(2, j, nout2) = s1 + r2
1322 27531430 : zout(1, j, nout3) = r1 + s2
1323 28924409 : zout(2, j, nout3) = s1 - r2
1324 : END DO; END DO
1325 : ELSE
1326 928274 : nin1 = ia - after
1327 928274 : nout1 = ia - atn
1328 2239998 : DO ib = 1, before
1329 1311724 : nin1 = nin1 + after
1330 1311724 : nin2 = nin1 + atb
1331 1311724 : nin3 = nin2 + atb
1332 1311724 : nout1 = nout1 + atn
1333 1311724 : nout2 = nout1 + after
1334 1311724 : nout3 = nout2 + after
1335 28046334 : DO j = 1, nfft
1336 25806336 : r1 = zin(1, j, nin1)
1337 25806336 : s1 = zin(2, j, nin1)
1338 25806336 : r = zin(1, j, nin2)
1339 25806336 : s = zin(2, j, nin2)
1340 25806336 : r2 = (r + s)*rt2i
1341 25806336 : s2 = (s - r)*rt2i
1342 25806336 : r3 = zin(2, j, nin3)
1343 25806336 : s3 = zin(1, j, nin3)
1344 25806336 : r = r2 + r3
1345 25806336 : s = s2 - s3
1346 25806336 : zout(1, j, nout1) = r + r1
1347 25806336 : zout(2, j, nout1) = s + s1
1348 25806336 : r1 = r1 - .5_dp*r
1349 25806336 : s1 = s1 - .5_dp*s
1350 25806336 : r2 = bb*(r2 - r3)
1351 25806336 : s2 = bb*(s2 + s3)
1352 25806336 : zout(1, j, nout2) = r1 - s2
1353 25806336 : zout(2, j, nout2) = s1 + r2
1354 25806336 : zout(1, j, nout3) = r1 + s2
1355 27118060 : zout(2, j, nout3) = s1 - r2
1356 : END DO; END DO
1357 : END IF
1358 : ELSE
1359 144582040 : itt = ias*before
1360 144582040 : itrig = itt + 1
1361 144582040 : cr2 = trig(1, itrig)
1362 144582040 : ci2 = trig(2, itrig)
1363 144582040 : itrig = itrig + itt
1364 144582040 : cr3 = trig(1, itrig)
1365 144582040 : ci3 = trig(2, itrig)
1366 144582040 : nin1 = ia - after
1367 144582040 : nout1 = ia - atn
1368 347588130 : DO ib = 1, before
1369 203006090 : nin1 = nin1 + after
1370 203006090 : nin2 = nin1 + atb
1371 203006090 : nin3 = nin2 + atb
1372 203006090 : nout1 = nout1 + atn
1373 203006090 : nout2 = nout1 + after
1374 203006090 : nout3 = nout2 + after
1375 4051267928 : DO j = 1, nfft
1376 3703679798 : r1 = zin(1, j, nin1)
1377 3703679798 : s1 = zin(2, j, nin1)
1378 3703679798 : r = zin(1, j, nin2)
1379 3703679798 : s = zin(2, j, nin2)
1380 3703679798 : r2 = r*cr2 - s*ci2
1381 3703679798 : s2 = r*ci2 + s*cr2
1382 3703679798 : r = zin(1, j, nin3)
1383 3703679798 : s = zin(2, j, nin3)
1384 3703679798 : r3 = r*cr3 - s*ci3
1385 3703679798 : s3 = r*ci3 + s*cr3
1386 3703679798 : r = r2 + r3
1387 3703679798 : s = s2 + s3
1388 3703679798 : zout(1, j, nout1) = r + r1
1389 3703679798 : zout(2, j, nout1) = s + s1
1390 3703679798 : r1 = r1 - .5_dp*r
1391 3703679798 : s1 = s1 - .5_dp*s
1392 3703679798 : r2 = bb*(r2 - r3)
1393 3703679798 : s2 = bb*(s2 - s3)
1394 3703679798 : zout(1, j, nout2) = r1 - s2
1395 3703679798 : zout(2, j, nout2) = s1 + r2
1396 3703679798 : zout(1, j, nout3) = r1 + s2
1397 3906685888 : zout(2, j, nout3) = s1 - r2
1398 : END DO; END DO
1399 : END IF
1400 8840676 : 3000 CONTINUE
1401 : ELSE IF (now .EQ. 5) THEN
1402 : ! cos(2._dp*pi/5._dp)
1403 3306964 : cos2 = 0.3090169943749474_dp
1404 : ! cos(4._dp*pi/5._dp)
1405 3306964 : cos4 = -0.8090169943749474_dp
1406 : ! sin(2._dp*pi/5._dp)
1407 3306964 : sin2 = isign*0.9510565162951536_dp
1408 : ! sin(4._dp*pi/5._dp)
1409 3306964 : sin4 = isign*0.5877852522924731_dp
1410 3306964 : ia = 1
1411 3306964 : nin1 = ia - after
1412 3306964 : nout1 = ia - atn
1413 33873470 : DO ib = 1, before
1414 30566506 : nin1 = nin1 + after
1415 30566506 : nin2 = nin1 + atb
1416 30566506 : nin3 = nin2 + atb
1417 30566506 : nin4 = nin3 + atb
1418 30566506 : nin5 = nin4 + atb
1419 30566506 : nout1 = nout1 + atn
1420 30566506 : nout2 = nout1 + after
1421 30566506 : nout3 = nout2 + after
1422 30566506 : nout4 = nout3 + after
1423 30566506 : nout5 = nout4 + after
1424 649230495 : DO j = 1, nfft
1425 615357025 : r1 = zin(1, j, nin1)
1426 615357025 : s1 = zin(2, j, nin1)
1427 615357025 : r2 = zin(1, j, nin2)
1428 615357025 : s2 = zin(2, j, nin2)
1429 615357025 : r3 = zin(1, j, nin3)
1430 615357025 : s3 = zin(2, j, nin3)
1431 615357025 : r4 = zin(1, j, nin4)
1432 615357025 : s4 = zin(2, j, nin4)
1433 615357025 : r5 = zin(1, j, nin5)
1434 615357025 : s5 = zin(2, j, nin5)
1435 615357025 : r25 = r2 + r5
1436 615357025 : r34 = r3 + r4
1437 615357025 : s25 = s2 - s5
1438 615357025 : s34 = s3 - s4
1439 615357025 : zout(1, j, nout1) = r1 + r25 + r34
1440 615357025 : r = r1 + cos2*r25 + cos4*r34
1441 615357025 : s = sin2*s25 + sin4*s34
1442 615357025 : zout(1, j, nout2) = r - s
1443 615357025 : zout(1, j, nout5) = r + s
1444 615357025 : r = r1 + cos4*r25 + cos2*r34
1445 615357025 : s = sin4*s25 - sin2*s34
1446 615357025 : zout(1, j, nout3) = r - s
1447 615357025 : zout(1, j, nout4) = r + s
1448 615357025 : r25 = r2 - r5
1449 615357025 : r34 = r3 - r4
1450 615357025 : s25 = s2 + s5
1451 615357025 : s34 = s3 + s4
1452 615357025 : zout(2, j, nout1) = s1 + s25 + s34
1453 615357025 : r = s1 + cos2*s25 + cos4*s34
1454 615357025 : s = sin2*r25 + sin4*r34
1455 615357025 : zout(2, j, nout2) = r + s
1456 615357025 : zout(2, j, nout5) = r - s
1457 615357025 : r = s1 + cos4*s25 + cos2*s34
1458 615357025 : s = sin4*r25 - sin2*r34
1459 615357025 : zout(2, j, nout3) = r + s
1460 645923531 : zout(2, j, nout4) = r - s
1461 : END DO; END DO
1462 11500572 : DO 5000, ia = 2, after
1463 8193608 : ias = ia - 1
1464 8193608 : IF (8*ias .EQ. 5*after) THEN
1465 615304 : IF (isign .EQ. 1) THEN
1466 320786 : nin1 = ia - after
1467 320786 : nout1 = ia - atn
1468 941692 : DO ib = 1, before
1469 620906 : nin1 = nin1 + after
1470 620906 : nin2 = nin1 + atb
1471 620906 : nin3 = nin2 + atb
1472 620906 : nin4 = nin3 + atb
1473 620906 : nin5 = nin4 + atb
1474 620906 : nout1 = nout1 + atn
1475 620906 : nout2 = nout1 + after
1476 620906 : nout3 = nout2 + after
1477 620906 : nout4 = nout3 + after
1478 620906 : nout5 = nout4 + after
1479 14525079 : DO j = 1, nfft
1480 13583387 : r1 = zin(1, j, nin1)
1481 13583387 : s1 = zin(2, j, nin1)
1482 13583387 : r = zin(1, j, nin2)
1483 13583387 : s = zin(2, j, nin2)
1484 13583387 : r2 = (r - s)*rt2i
1485 13583387 : s2 = (r + s)*rt2i
1486 13583387 : r3 = zin(2, j, nin3)
1487 13583387 : s3 = zin(1, j, nin3)
1488 13583387 : r = zin(1, j, nin4)
1489 13583387 : s = zin(2, j, nin4)
1490 13583387 : r4 = (r + s)*rt2i
1491 13583387 : s4 = (r - s)*rt2i
1492 13583387 : r5 = zin(1, j, nin5)
1493 13583387 : s5 = zin(2, j, nin5)
1494 13583387 : r25 = r2 - r5
1495 13583387 : r34 = r3 + r4
1496 13583387 : s25 = s2 + s5
1497 13583387 : s34 = s3 - s4
1498 13583387 : zout(1, j, nout1) = r1 + r25 - r34
1499 13583387 : r = r1 + cos2*r25 - cos4*r34
1500 13583387 : s = sin2*s25 + sin4*s34
1501 13583387 : zout(1, j, nout2) = r - s
1502 13583387 : zout(1, j, nout5) = r + s
1503 13583387 : r = r1 + cos4*r25 - cos2*r34
1504 13583387 : s = sin4*s25 - sin2*s34
1505 13583387 : zout(1, j, nout3) = r - s
1506 13583387 : zout(1, j, nout4) = r + s
1507 13583387 : r25 = r2 + r5
1508 13583387 : r34 = r4 - r3
1509 13583387 : s25 = s2 - s5
1510 13583387 : s34 = s3 + s4
1511 13583387 : zout(2, j, nout1) = s1 + s25 + s34
1512 13583387 : r = s1 + cos2*s25 + cos4*s34
1513 13583387 : s = sin2*r25 + sin4*r34
1514 13583387 : zout(2, j, nout2) = r + s
1515 13583387 : zout(2, j, nout5) = r - s
1516 13583387 : r = s1 + cos4*s25 + cos2*s34
1517 13583387 : s = sin4*r25 - sin2*r34
1518 13583387 : zout(2, j, nout3) = r + s
1519 14204293 : zout(2, j, nout4) = r - s
1520 : END DO; END DO
1521 : ELSE
1522 294518 : nin1 = ia - after
1523 294518 : nout1 = ia - atn
1524 875980 : DO ib = 1, before
1525 581462 : nin1 = nin1 + after
1526 581462 : nin2 = nin1 + atb
1527 581462 : nin3 = nin2 + atb
1528 581462 : nin4 = nin3 + atb
1529 581462 : nin5 = nin4 + atb
1530 581462 : nout1 = nout1 + atn
1531 581462 : nout2 = nout1 + after
1532 581462 : nout3 = nout2 + after
1533 581462 : nout4 = nout3 + after
1534 581462 : nout5 = nout4 + after
1535 13374060 : DO j = 1, nfft
1536 12498080 : r1 = zin(1, j, nin1)
1537 12498080 : s1 = zin(2, j, nin1)
1538 12498080 : r = zin(1, j, nin2)
1539 12498080 : s = zin(2, j, nin2)
1540 12498080 : r2 = (r + s)*rt2i
1541 12498080 : s2 = (s - r)*rt2i
1542 12498080 : r3 = zin(2, j, nin3)
1543 12498080 : s3 = zin(1, j, nin3)
1544 12498080 : r = zin(1, j, nin4)
1545 12498080 : s = zin(2, j, nin4)
1546 12498080 : r4 = (s - r)*rt2i
1547 12498080 : s4 = (r + s)*rt2i
1548 12498080 : r5 = zin(1, j, nin5)
1549 12498080 : s5 = zin(2, j, nin5)
1550 12498080 : r25 = r2 - r5
1551 12498080 : r34 = r3 + r4
1552 12498080 : s25 = s2 + s5
1553 12498080 : s34 = s4 - s3
1554 12498080 : zout(1, j, nout1) = r1 + r25 + r34
1555 12498080 : r = r1 + cos2*r25 + cos4*r34
1556 12498080 : s = sin2*s25 + sin4*s34
1557 12498080 : zout(1, j, nout2) = r - s
1558 12498080 : zout(1, j, nout5) = r + s
1559 12498080 : r = r1 + cos4*r25 + cos2*r34
1560 12498080 : s = sin4*s25 - sin2*s34
1561 12498080 : zout(1, j, nout3) = r - s
1562 12498080 : zout(1, j, nout4) = r + s
1563 12498080 : r25 = r2 + r5
1564 12498080 : r34 = r3 - r4
1565 12498080 : s25 = s2 - s5
1566 12498080 : s34 = s3 + s4
1567 12498080 : zout(2, j, nout1) = s1 + s25 - s34
1568 12498080 : r = s1 + cos2*s25 - cos4*s34
1569 12498080 : s = sin2*r25 + sin4*r34
1570 12498080 : zout(2, j, nout2) = r + s
1571 12498080 : zout(2, j, nout5) = r - s
1572 12498080 : r = s1 + cos4*s25 - cos2*s34
1573 12498080 : s = sin4*r25 - sin2*r34
1574 12498080 : zout(2, j, nout3) = r + s
1575 13079542 : zout(2, j, nout4) = r - s
1576 : END DO; END DO
1577 : END IF
1578 : ELSE
1579 7578304 : ias = ia - 1
1580 7578304 : itt = ias*before
1581 7578304 : itrig = itt + 1
1582 7578304 : cr2 = trig(1, itrig)
1583 7578304 : ci2 = trig(2, itrig)
1584 7578304 : itrig = itrig + itt
1585 7578304 : cr3 = trig(1, itrig)
1586 7578304 : ci3 = trig(2, itrig)
1587 7578304 : itrig = itrig + itt
1588 7578304 : cr4 = trig(1, itrig)
1589 7578304 : ci4 = trig(2, itrig)
1590 7578304 : itrig = itrig + itt
1591 7578304 : cr5 = trig(1, itrig)
1592 7578304 : ci5 = trig(2, itrig)
1593 7578304 : nin1 = ia - after
1594 7578304 : nout1 = ia - atn
1595 24476272 : DO ib = 1, before
1596 16897968 : nin1 = nin1 + after
1597 16897968 : nin2 = nin1 + atb
1598 16897968 : nin3 = nin2 + atb
1599 16897968 : nin4 = nin3 + atb
1600 16897968 : nin5 = nin4 + atb
1601 16897968 : nout1 = nout1 + atn
1602 16897968 : nout2 = nout1 + after
1603 16897968 : nout3 = nout2 + after
1604 16897968 : nout4 = nout3 + after
1605 16897968 : nout5 = nout4 + after
1606 354228874 : DO j = 1, nfft
1607 329752602 : r1 = zin(1, j, nin1)
1608 329752602 : s1 = zin(2, j, nin1)
1609 329752602 : r = zin(1, j, nin2)
1610 329752602 : s = zin(2, j, nin2)
1611 329752602 : r2 = r*cr2 - s*ci2
1612 329752602 : s2 = r*ci2 + s*cr2
1613 329752602 : r = zin(1, j, nin3)
1614 329752602 : s = zin(2, j, nin3)
1615 329752602 : r3 = r*cr3 - s*ci3
1616 329752602 : s3 = r*ci3 + s*cr3
1617 329752602 : r = zin(1, j, nin4)
1618 329752602 : s = zin(2, j, nin4)
1619 329752602 : r4 = r*cr4 - s*ci4
1620 329752602 : s4 = r*ci4 + s*cr4
1621 329752602 : r = zin(1, j, nin5)
1622 329752602 : s = zin(2, j, nin5)
1623 329752602 : r5 = r*cr5 - s*ci5
1624 329752602 : s5 = r*ci5 + s*cr5
1625 329752602 : r25 = r2 + r5
1626 329752602 : r34 = r3 + r4
1627 329752602 : s25 = s2 - s5
1628 329752602 : s34 = s3 - s4
1629 329752602 : zout(1, j, nout1) = r1 + r25 + r34
1630 329752602 : r = r1 + cos2*r25 + cos4*r34
1631 329752602 : s = sin2*s25 + sin4*s34
1632 329752602 : zout(1, j, nout2) = r - s
1633 329752602 : zout(1, j, nout5) = r + s
1634 329752602 : r = r1 + cos4*r25 + cos2*r34
1635 329752602 : s = sin4*s25 - sin2*s34
1636 329752602 : zout(1, j, nout3) = r - s
1637 329752602 : zout(1, j, nout4) = r + s
1638 329752602 : r25 = r2 - r5
1639 329752602 : r34 = r3 - r4
1640 329752602 : s25 = s2 + s5
1641 329752602 : s34 = s3 + s4
1642 329752602 : zout(2, j, nout1) = s1 + s25 + s34
1643 329752602 : r = s1 + cos2*s25 + cos4*s34
1644 329752602 : s = sin2*r25 + sin4*r34
1645 329752602 : zout(2, j, nout2) = r + s
1646 329752602 : zout(2, j, nout5) = r - s
1647 329752602 : r = s1 + cos4*s25 + cos2*s34
1648 329752602 : s = sin4*r25 - sin2*r34
1649 329752602 : zout(2, j, nout3) = r + s
1650 346650570 : zout(2, j, nout4) = r - s
1651 : END DO; END DO
1652 : END IF
1653 3306964 : 5000 CONTINUE
1654 : ELSE IF (now .EQ. 6) THEN
1655 : ! .5_dp*sqrt(3._dp)
1656 2709167 : bb = isign*0.8660254037844387_dp
1657 :
1658 2709167 : ia = 1
1659 2709167 : nin1 = ia - after
1660 2709167 : nout1 = ia - atn
1661 35774510 : DO ib = 1, before
1662 33065343 : nin1 = nin1 + after
1663 33065343 : nin2 = nin1 + atb
1664 33065343 : nin3 = nin2 + atb
1665 33065343 : nin4 = nin3 + atb
1666 33065343 : nin5 = nin4 + atb
1667 33065343 : nin6 = nin5 + atb
1668 33065343 : nout1 = nout1 + atn
1669 33065343 : nout2 = nout1 + after
1670 33065343 : nout3 = nout2 + after
1671 33065343 : nout4 = nout3 + after
1672 33065343 : nout5 = nout4 + after
1673 33065343 : nout6 = nout5 + after
1674 595327526 : DO j = 1, nfft
1675 559553016 : r2 = zin(1, j, nin3)
1676 559553016 : s2 = zin(2, j, nin3)
1677 559553016 : r3 = zin(1, j, nin5)
1678 559553016 : s3 = zin(2, j, nin5)
1679 559553016 : r = r2 + r3
1680 559553016 : s = s2 + s3
1681 559553016 : r1 = zin(1, j, nin1)
1682 559553016 : s1 = zin(2, j, nin1)
1683 559553016 : ur1 = r + r1
1684 559553016 : ui1 = s + s1
1685 559553016 : r1 = r1 - .5_dp*r
1686 559553016 : s1 = s1 - .5_dp*s
1687 559553016 : r = r2 - r3
1688 559553016 : s = s2 - s3
1689 559553016 : ur2 = r1 - s*bb
1690 559553016 : ui2 = s1 + r*bb
1691 559553016 : ur3 = r1 + s*bb
1692 559553016 : ui3 = s1 - r*bb
1693 :
1694 559553016 : r2 = zin(1, j, nin6)
1695 559553016 : s2 = zin(2, j, nin6)
1696 559553016 : r3 = zin(1, j, nin2)
1697 559553016 : s3 = zin(2, j, nin2)
1698 559553016 : r = r2 + r3
1699 559553016 : s = s2 + s3
1700 559553016 : r1 = zin(1, j, nin4)
1701 559553016 : s1 = zin(2, j, nin4)
1702 559553016 : vr1 = r + r1
1703 559553016 : vi1 = s + s1
1704 559553016 : r1 = r1 - .5_dp*r
1705 559553016 : s1 = s1 - .5_dp*s
1706 559553016 : r = r2 - r3
1707 559553016 : s = s2 - s3
1708 559553016 : vr2 = r1 - s*bb
1709 559553016 : vi2 = s1 + r*bb
1710 559553016 : vr3 = r1 + s*bb
1711 559553016 : vi3 = s1 - r*bb
1712 :
1713 559553016 : zout(1, j, nout1) = ur1 + vr1
1714 559553016 : zout(2, j, nout1) = ui1 + vi1
1715 559553016 : zout(1, j, nout5) = ur2 + vr2
1716 559553016 : zout(2, j, nout5) = ui2 + vi2
1717 559553016 : zout(1, j, nout3) = ur3 + vr3
1718 559553016 : zout(2, j, nout3) = ui3 + vi3
1719 559553016 : zout(1, j, nout4) = ur1 - vr1
1720 559553016 : zout(2, j, nout4) = ui1 - vi1
1721 559553016 : zout(1, j, nout2) = ur2 - vr2
1722 559553016 : zout(2, j, nout2) = ui2 - vi2
1723 559553016 : zout(1, j, nout6) = ur3 - vr3
1724 592618359 : zout(2, j, nout6) = ui3 - vi3
1725 : END DO; END DO
1726 : ELSE
1727 0 : CPABORT("error fftstp")
1728 : END IF
1729 :
1730 22651284 : END SUBROUTINE fftstp
1731 :
1732 : END MODULE ps_wavelet_fft3d
|
