Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 110819 : 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 1801626 : loop_data: DO i = 1, ndata1024
54 1801626 : IF (n <= idata(i)) THEN
55 110819 : n_next = idata(i)
56 110819 : 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 97697 : 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 1602174 : DO i = 1, 150
202 1602174 : 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 1602174 : IF (n .EQ. idata(1, i)) THEN
209 97697 : ic = 0
210 334211 : DO j = 1, 6
211 334211 : itt = idata(1 + j, i)
212 334211 : IF (itt .GT. 1) THEN
213 236514 : ic = ic + 1
214 236514 : 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 97697 : after(1) = 1
224 97697 : before(ic) = 1
225 236514 : DO i = 2, ic
226 138817 : after(i) = after(i - 1)*now(i - 1)
227 236514 : before(ic - i + 1) = before(ic - i + 2)*now(ic - i + 2)
228 : END DO
229 :
230 97697 : twopi = 6.283185307179586_dp
231 97697 : angle = isign*twopi/n
232 97697 : IF (MOD(n, 2) .EQ. 0) THEN
233 86282 : nh = n/2
234 86282 : trig(1, 1) = 1._dp
235 86282 : trig(2, 1) = 0._dp
236 86282 : trig(1, nh + 1) = -1._dp
237 86282 : trig(2, nh + 1) = 0._dp
238 1914884 : DO 40, i = 1, nh - 1
239 1828602 : trigc = COS(i*angle)
240 1828602 : trigs = SIN(i*angle)
241 1828602 : trig(1, i + 1) = trigc
242 1828602 : trig(2, i + 1) = trigs
243 1828602 : trig(1, n - i + 1) = trigc
244 1828602 : trig(2, n - i + 1) = -trigs
245 86282 : 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 97697 : 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 23211676 : 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 23211676 : atn = after*now
300 23211676 : atb = after*before
301 :
302 : ! sqrt(.5_dp)
303 23211676 : 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 5491527 : IF (isign .EQ. 1) THEN
474 2814215 : ia = 1
475 2814215 : nin1 = ia - after
476 2814215 : nout1 = ia - atn
477 24585474 : DO ib = 1, before
478 21771259 : nin1 = nin1 + after
479 21771259 : nin2 = nin1 + atb
480 21771259 : nin3 = nin2 + atb
481 21771259 : nin4 = nin3 + atb
482 21771259 : nout1 = nout1 + atn
483 21771259 : nout2 = nout1 + after
484 21771259 : nout3 = nout2 + after
485 21771259 : nout4 = nout3 + after
486 446887100 : DO j = 1, nfft
487 422301626 : r1 = zin(1, j, nin1)
488 422301626 : s1 = zin(2, j, nin1)
489 422301626 : r2 = zin(1, j, nin2)
490 422301626 : s2 = zin(2, j, nin2)
491 422301626 : r3 = zin(1, j, nin3)
492 422301626 : s3 = zin(2, j, nin3)
493 422301626 : r4 = zin(1, j, nin4)
494 422301626 : s4 = zin(2, j, nin4)
495 422301626 : r = r1 + r3
496 422301626 : s = r2 + r4
497 422301626 : zout(1, j, nout1) = r + s
498 422301626 : zout(1, j, nout3) = r - s
499 422301626 : r = r1 - r3
500 422301626 : s = s2 - s4
501 422301626 : zout(1, j, nout2) = r - s
502 422301626 : zout(1, j, nout4) = r + s
503 422301626 : r = s1 + s3
504 422301626 : s = s2 + s4
505 422301626 : zout(2, j, nout1) = r + s
506 422301626 : zout(2, j, nout3) = r - s
507 422301626 : r = s1 - s3
508 422301626 : s = r2 - r4
509 422301626 : zout(2, j, nout2) = r + s
510 444072885 : zout(2, j, nout4) = r - s
511 : END DO; END DO
512 28587546 : DO 4000, ia = 2, after
513 25773331 : ias = ia - 1
514 25773331 : IF (2*ias .EQ. after) THEN
515 1174933 : nin1 = ia - after
516 1174933 : nout1 = ia - atn
517 2757164 : DO ib = 1, before
518 1582231 : nin1 = nin1 + after
519 1582231 : nin2 = nin1 + atb
520 1582231 : nin3 = nin2 + atb
521 1582231 : nin4 = nin3 + atb
522 1582231 : nout1 = nout1 + atn
523 1582231 : nout2 = nout1 + after
524 1582231 : nout3 = nout2 + after
525 1582231 : nout4 = nout3 + after
526 33184260 : DO j = 1, nfft
527 30427096 : r1 = zin(1, j, nin1)
528 30427096 : s1 = zin(2, j, nin1)
529 30427096 : r = zin(1, j, nin2)
530 30427096 : s = zin(2, j, nin2)
531 30427096 : r2 = (r - s)*rt2i
532 30427096 : s2 = (r + s)*rt2i
533 30427096 : r3 = zin(2, j, nin3)
534 30427096 : s3 = zin(1, j, nin3)
535 30427096 : r = zin(1, j, nin4)
536 30427096 : s = zin(2, j, nin4)
537 30427096 : r4 = (r + s)*rt2i
538 30427096 : s4 = (r - s)*rt2i
539 30427096 : r = r1 - r3
540 30427096 : s = r2 - r4
541 30427096 : zout(1, j, nout1) = r + s
542 30427096 : zout(1, j, nout3) = r - s
543 30427096 : r = r1 + r3
544 30427096 : s = s2 - s4
545 30427096 : zout(1, j, nout2) = r - s
546 30427096 : zout(1, j, nout4) = r + s
547 30427096 : r = s1 + s3
548 30427096 : s = s2 + s4
549 30427096 : zout(2, j, nout1) = r + s
550 30427096 : zout(2, j, nout3) = r - s
551 30427096 : r = s1 - s3
552 30427096 : s = r2 + r4
553 30427096 : zout(2, j, nout2) = r + s
554 32009327 : zout(2, j, nout4) = r - s
555 : END DO; END DO
556 : ELSE
557 24598398 : itt = ias*before
558 24598398 : itrig = itt + 1
559 24598398 : cr2 = trig(1, itrig)
560 24598398 : ci2 = trig(2, itrig)
561 24598398 : itrig = itrig + itt
562 24598398 : cr3 = trig(1, itrig)
563 24598398 : ci3 = trig(2, itrig)
564 24598398 : itrig = itrig + itt
565 24598398 : cr4 = trig(1, itrig)
566 24598398 : ci4 = trig(2, itrig)
567 24598398 : nin1 = ia - after
568 24598398 : nout1 = ia - atn
569 61438824 : DO ib = 1, before
570 36840426 : nin1 = nin1 + after
571 36840426 : nin2 = nin1 + atb
572 36840426 : nin3 = nin2 + atb
573 36840426 : nin4 = nin3 + atb
574 36840426 : nout1 = nout1 + atn
575 36840426 : nout2 = nout1 + after
576 36840426 : nout3 = nout2 + after
577 36840426 : nout4 = nout3 + after
578 758184028 : DO j = 1, nfft
579 696745204 : r1 = zin(1, j, nin1)
580 696745204 : s1 = zin(2, j, nin1)
581 696745204 : r = zin(1, j, nin2)
582 696745204 : s = zin(2, j, nin2)
583 696745204 : r2 = r*cr2 - s*ci2
584 696745204 : s2 = r*ci2 + s*cr2
585 696745204 : r = zin(1, j, nin3)
586 696745204 : s = zin(2, j, nin3)
587 696745204 : r3 = r*cr3 - s*ci3
588 696745204 : s3 = r*ci3 + s*cr3
589 696745204 : r = zin(1, j, nin4)
590 696745204 : s = zin(2, j, nin4)
591 696745204 : r4 = r*cr4 - s*ci4
592 696745204 : s4 = r*ci4 + s*cr4
593 696745204 : r = r1 + r3
594 696745204 : s = r2 + r4
595 696745204 : zout(1, j, nout1) = r + s
596 696745204 : zout(1, j, nout3) = r - s
597 696745204 : r = r1 - r3
598 696745204 : s = s2 - s4
599 696745204 : zout(1, j, nout2) = r - s
600 696745204 : zout(1, j, nout4) = r + s
601 696745204 : r = s1 + s3
602 696745204 : s = s2 + s4
603 696745204 : zout(2, j, nout1) = r + s
604 696745204 : zout(2, j, nout3) = r - s
605 696745204 : r = s1 - s3
606 696745204 : s = r2 - r4
607 696745204 : zout(2, j, nout2) = r + s
608 733585630 : zout(2, j, nout4) = r - s
609 : END DO; END DO
610 : END IF
611 2814215 : 4000 CONTINUE
612 : ELSE
613 2677312 : ia = 1
614 2677312 : nin1 = ia - after
615 2677312 : nout1 = ia - atn
616 23311951 : DO ib = 1, before
617 20634639 : nin1 = nin1 + after
618 20634639 : nin2 = nin1 + atb
619 20634639 : nin3 = nin2 + atb
620 20634639 : nin4 = nin3 + atb
621 20634639 : nout1 = nout1 + atn
622 20634639 : nout2 = nout1 + after
623 20634639 : nout3 = nout2 + after
624 20634639 : nout4 = nout3 + after
625 424733253 : DO j = 1, nfft
626 401421302 : r1 = zin(1, j, nin1)
627 401421302 : s1 = zin(2, j, nin1)
628 401421302 : r2 = zin(1, j, nin2)
629 401421302 : s2 = zin(2, j, nin2)
630 401421302 : r3 = zin(1, j, nin3)
631 401421302 : s3 = zin(2, j, nin3)
632 401421302 : r4 = zin(1, j, nin4)
633 401421302 : s4 = zin(2, j, nin4)
634 401421302 : r = r1 + r3
635 401421302 : s = r2 + r4
636 401421302 : zout(1, j, nout1) = r + s
637 401421302 : zout(1, j, nout3) = r - s
638 401421302 : r = r1 - r3
639 401421302 : s = s2 - s4
640 401421302 : zout(1, j, nout2) = r + s
641 401421302 : zout(1, j, nout4) = r - s
642 401421302 : r = s1 + s3
643 401421302 : s = s2 + s4
644 401421302 : zout(2, j, nout1) = r + s
645 401421302 : zout(2, j, nout3) = r - s
646 401421302 : r = s1 - s3
647 401421302 : s = r2 - r4
648 401421302 : zout(2, j, nout2) = r - s
649 422055941 : zout(2, j, nout4) = r + s
650 : END DO; END DO
651 26941795 : DO 4100, ia = 2, after
652 24264483 : ias = ia - 1
653 24264483 : IF (2*ias .EQ. after) THEN
654 1112793 : nin1 = ia - after
655 1112793 : nout1 = ia - atn
656 2595516 : DO ib = 1, before
657 1482723 : nin1 = nin1 + after
658 1482723 : nin2 = nin1 + atb
659 1482723 : nin3 = nin2 + atb
660 1482723 : nin4 = nin3 + atb
661 1482723 : nout1 = nout1 + atn
662 1482723 : nout2 = nout1 + after
663 1482723 : nout3 = nout2 + after
664 1482723 : nout4 = nout3 + after
665 31135178 : DO j = 1, nfft
666 28539662 : r1 = zin(1, j, nin1)
667 28539662 : s1 = zin(2, j, nin1)
668 28539662 : r = zin(1, j, nin2)
669 28539662 : s = zin(2, j, nin2)
670 28539662 : r2 = (r + s)*rt2i
671 28539662 : s2 = (s - r)*rt2i
672 28539662 : r3 = zin(2, j, nin3)
673 28539662 : s3 = zin(1, j, nin3)
674 28539662 : r = zin(1, j, nin4)
675 28539662 : s = zin(2, j, nin4)
676 28539662 : r4 = (s - r)*rt2i
677 28539662 : s4 = (r + s)*rt2i
678 28539662 : r = r1 + r3
679 28539662 : s = r2 + r4
680 28539662 : zout(1, j, nout1) = r + s
681 28539662 : zout(1, j, nout3) = r - s
682 28539662 : r = r1 - r3
683 28539662 : s = s2 + s4
684 28539662 : zout(1, j, nout2) = r + s
685 28539662 : zout(1, j, nout4) = r - s
686 28539662 : r = s1 - s3
687 28539662 : s = s2 - s4
688 28539662 : zout(2, j, nout1) = r + s
689 28539662 : zout(2, j, nout3) = r - s
690 28539662 : r = s1 + s3
691 28539662 : s = r2 - r4
692 28539662 : zout(2, j, nout2) = r - s
693 30022385 : zout(2, j, nout4) = r + s
694 : END DO; END DO
695 : ELSE
696 23151690 : itt = ias*before
697 23151690 : itrig = itt + 1
698 23151690 : cr2 = trig(1, itrig)
699 23151690 : ci2 = trig(2, itrig)
700 23151690 : itrig = itrig + itt
701 23151690 : cr3 = trig(1, itrig)
702 23151690 : ci3 = trig(2, itrig)
703 23151690 : itrig = itrig + itt
704 23151690 : cr4 = trig(1, itrig)
705 23151690 : ci4 = trig(2, itrig)
706 23151690 : nin1 = ia - after
707 23151690 : nout1 = ia - atn
708 57948264 : DO ib = 1, before
709 34796574 : nin1 = nin1 + after
710 34796574 : nin2 = nin1 + atb
711 34796574 : nin3 = nin2 + atb
712 34796574 : nin4 = nin3 + atb
713 34796574 : nout1 = nout1 + atn
714 34796574 : nout2 = nout1 + after
715 34796574 : nout3 = nout2 + after
716 34796574 : nout4 = nout3 + after
717 719687556 : DO j = 1, nfft
718 661739292 : r1 = zin(1, j, nin1)
719 661739292 : s1 = zin(2, j, nin1)
720 661739292 : r = zin(1, j, nin2)
721 661739292 : s = zin(2, j, nin2)
722 661739292 : r2 = r*cr2 - s*ci2
723 661739292 : s2 = r*ci2 + s*cr2
724 661739292 : r = zin(1, j, nin3)
725 661739292 : s = zin(2, j, nin3)
726 661739292 : r3 = r*cr3 - s*ci3
727 661739292 : s3 = r*ci3 + s*cr3
728 661739292 : r = zin(1, j, nin4)
729 661739292 : s = zin(2, j, nin4)
730 661739292 : r4 = r*cr4 - s*ci4
731 661739292 : s4 = r*ci4 + s*cr4
732 661739292 : r = r1 + r3
733 661739292 : s = r2 + r4
734 661739292 : zout(1, j, nout1) = r + s
735 661739292 : zout(1, j, nout3) = r - s
736 661739292 : r = r1 - r3
737 661739292 : s = s2 - s4
738 661739292 : zout(1, j, nout2) = r + s
739 661739292 : zout(1, j, nout4) = r - s
740 661739292 : r = s1 + s3
741 661739292 : s = s2 + s4
742 661739292 : zout(2, j, nout1) = r + s
743 661739292 : zout(2, j, nout3) = r - s
744 661739292 : r = s1 - s3
745 661739292 : s = r2 - r4
746 661739292 : zout(2, j, nout2) = r - s
747 696535866 : zout(2, j, nout4) = r + s
748 : END DO; END DO
749 : END IF
750 2677312 : 4100 CONTINUE
751 : END IF
752 : ELSE IF (now .EQ. 8) THEN
753 2533208 : IF (isign .EQ. -1) THEN
754 1220866 : ia = 1
755 1220866 : nin1 = ia - after
756 1220866 : nout1 = ia - atn
757 10348846 : DO ib = 1, before
758 9127980 : nin1 = nin1 + after
759 9127980 : nin2 = nin1 + atb
760 9127980 : nin3 = nin2 + atb
761 9127980 : nin4 = nin3 + atb
762 9127980 : nin5 = nin4 + atb
763 9127980 : nin6 = nin5 + atb
764 9127980 : nin7 = nin6 + atb
765 9127980 : nin8 = nin7 + atb
766 9127980 : nout1 = nout1 + atn
767 9127980 : nout2 = nout1 + after
768 9127980 : nout3 = nout2 + after
769 9127980 : nout4 = nout3 + after
770 9127980 : nout5 = nout4 + after
771 9127980 : nout6 = nout5 + after
772 9127980 : nout7 = nout6 + after
773 9127980 : nout8 = nout7 + after
774 185635266 : DO j = 1, nfft
775 175286420 : r1 = zin(1, j, nin1)
776 175286420 : s1 = zin(2, j, nin1)
777 175286420 : r2 = zin(1, j, nin2)
778 175286420 : s2 = zin(2, j, nin2)
779 175286420 : r3 = zin(1, j, nin3)
780 175286420 : s3 = zin(2, j, nin3)
781 175286420 : r4 = zin(1, j, nin4)
782 175286420 : s4 = zin(2, j, nin4)
783 175286420 : r5 = zin(1, j, nin5)
784 175286420 : s5 = zin(2, j, nin5)
785 175286420 : r6 = zin(1, j, nin6)
786 175286420 : s6 = zin(2, j, nin6)
787 175286420 : r7 = zin(1, j, nin7)
788 175286420 : s7 = zin(2, j, nin7)
789 175286420 : r8 = zin(1, j, nin8)
790 175286420 : s8 = zin(2, j, nin8)
791 175286420 : r = r1 + r5
792 175286420 : s = r3 + r7
793 175286420 : ap = r + s
794 175286420 : am = r - s
795 175286420 : r = r2 + r6
796 175286420 : s = r4 + r8
797 175286420 : bp = r + s
798 175286420 : bm = r - s
799 175286420 : r = s1 + s5
800 175286420 : s = s3 + s7
801 175286420 : cp = r + s
802 175286420 : cm = r - s
803 175286420 : r = s2 + s6
804 175286420 : s = s4 + s8
805 175286420 : dpp = r + s
806 175286420 : dm = r - s
807 175286420 : zout(1, j, nout1) = ap + bp
808 175286420 : zout(2, j, nout1) = cp + dpp
809 175286420 : zout(1, j, nout5) = ap - bp
810 175286420 : zout(2, j, nout5) = cp - dpp
811 175286420 : zout(1, j, nout3) = am + dm
812 175286420 : zout(2, j, nout3) = cm - bm
813 175286420 : zout(1, j, nout7) = am - dm
814 175286420 : zout(2, j, nout7) = cm + bm
815 175286420 : r = r1 - r5
816 175286420 : s = s3 - s7
817 175286420 : ap = r + s
818 175286420 : am = r - s
819 175286420 : r = s1 - s5
820 175286420 : s = r3 - r7
821 175286420 : bp = r + s
822 175286420 : bm = r - s
823 175286420 : r = s4 - s8
824 175286420 : s = r2 - r6
825 175286420 : cp = r + s
826 175286420 : cm = r - s
827 175286420 : r = s2 - s6
828 175286420 : s = r4 - r8
829 175286420 : dpp = r + s
830 175286420 : dm = r - s
831 175286420 : r = (cp + dm)*rt2i
832 175286420 : s = (dm - cp)*rt2i
833 175286420 : cp = (cm + dpp)*rt2i
834 175286420 : dpp = (cm - dpp)*rt2i
835 175286420 : zout(1, j, nout2) = ap + r
836 175286420 : zout(2, j, nout2) = bm + s
837 175286420 : zout(1, j, nout6) = ap - r
838 175286420 : zout(2, j, nout6) = bm - s
839 175286420 : zout(1, j, nout4) = am + cp
840 175286420 : zout(2, j, nout4) = bp + dpp
841 175286420 : zout(1, j, nout8) = am - cp
842 184414400 : zout(2, j, nout8) = bp - dpp
843 : END DO; END DO
844 3533341 : DO 8000, ia = 2, after
845 2312475 : ias = ia - 1
846 2312475 : itt = ias*before
847 2312475 : itrig = itt + 1
848 2312475 : cr2 = trig(1, itrig)
849 2312475 : ci2 = trig(2, itrig)
850 2312475 : itrig = itrig + itt
851 2312475 : cr3 = trig(1, itrig)
852 2312475 : ci3 = trig(2, itrig)
853 2312475 : itrig = itrig + itt
854 2312475 : cr4 = trig(1, itrig)
855 2312475 : ci4 = trig(2, itrig)
856 2312475 : itrig = itrig + itt
857 2312475 : cr5 = trig(1, itrig)
858 2312475 : ci5 = trig(2, itrig)
859 2312475 : itrig = itrig + itt
860 2312475 : cr6 = trig(1, itrig)
861 2312475 : ci6 = trig(2, itrig)
862 2312475 : itrig = itrig + itt
863 2312475 : cr7 = trig(1, itrig)
864 2312475 : ci7 = trig(2, itrig)
865 2312475 : itrig = itrig + itt
866 2312475 : cr8 = trig(1, itrig)
867 2312475 : ci8 = trig(2, itrig)
868 2312475 : nin1 = ia - after
869 2312475 : nout1 = ia - atn
870 8563760 : DO ib = 1, before
871 6251285 : nin1 = nin1 + after
872 6251285 : nin2 = nin1 + atb
873 6251285 : nin3 = nin2 + atb
874 6251285 : nin4 = nin3 + atb
875 6251285 : nin5 = nin4 + atb
876 6251285 : nin6 = nin5 + atb
877 6251285 : nin7 = nin6 + atb
878 6251285 : nin8 = nin7 + atb
879 6251285 : nout1 = nout1 + atn
880 6251285 : nout2 = nout1 + after
881 6251285 : nout3 = nout2 + after
882 6251285 : nout4 = nout3 + after
883 6251285 : nout5 = nout4 + after
884 6251285 : nout6 = nout5 + after
885 6251285 : nout7 = nout6 + after
886 6251285 : nout8 = nout7 + after
887 92622360 : DO j = 1, nfft
888 84058600 : r1 = zin(1, j, nin1)
889 84058600 : s1 = zin(2, j, nin1)
890 84058600 : r = zin(1, j, nin2)
891 84058600 : s = zin(2, j, nin2)
892 84058600 : r2 = r*cr2 - s*ci2
893 84058600 : s2 = r*ci2 + s*cr2
894 84058600 : r = zin(1, j, nin3)
895 84058600 : s = zin(2, j, nin3)
896 84058600 : r3 = r*cr3 - s*ci3
897 84058600 : s3 = r*ci3 + s*cr3
898 84058600 : r = zin(1, j, nin4)
899 84058600 : s = zin(2, j, nin4)
900 84058600 : r4 = r*cr4 - s*ci4
901 84058600 : s4 = r*ci4 + s*cr4
902 84058600 : r = zin(1, j, nin5)
903 84058600 : s = zin(2, j, nin5)
904 84058600 : r5 = r*cr5 - s*ci5
905 84058600 : s5 = r*ci5 + s*cr5
906 84058600 : r = zin(1, j, nin6)
907 84058600 : s = zin(2, j, nin6)
908 84058600 : r6 = r*cr6 - s*ci6
909 84058600 : s6 = r*ci6 + s*cr6
910 84058600 : r = zin(1, j, nin7)
911 84058600 : s = zin(2, j, nin7)
912 84058600 : r7 = r*cr7 - s*ci7
913 84058600 : s7 = r*ci7 + s*cr7
914 84058600 : r = zin(1, j, nin8)
915 84058600 : s = zin(2, j, nin8)
916 84058600 : r8 = r*cr8 - s*ci8
917 84058600 : s8 = r*ci8 + s*cr8
918 84058600 : r = r1 + r5
919 84058600 : s = r3 + r7
920 84058600 : ap = r + s
921 84058600 : am = r - s
922 84058600 : r = r2 + r6
923 84058600 : s = r4 + r8
924 84058600 : bp = r + s
925 84058600 : bm = r - s
926 84058600 : r = s1 + s5
927 84058600 : s = s3 + s7
928 84058600 : cp = r + s
929 84058600 : cm = r - s
930 84058600 : r = s2 + s6
931 84058600 : s = s4 + s8
932 84058600 : dpp = r + s
933 84058600 : dm = r - s
934 84058600 : zout(1, j, nout1) = ap + bp
935 84058600 : zout(2, j, nout1) = cp + dpp
936 84058600 : zout(1, j, nout5) = ap - bp
937 84058600 : zout(2, j, nout5) = cp - dpp
938 84058600 : zout(1, j, nout3) = am + dm
939 84058600 : zout(2, j, nout3) = cm - bm
940 84058600 : zout(1, j, nout7) = am - dm
941 84058600 : zout(2, j, nout7) = cm + bm
942 84058600 : r = r1 - r5
943 84058600 : s = s3 - s7
944 84058600 : ap = r + s
945 84058600 : am = r - s
946 84058600 : r = s1 - s5
947 84058600 : s = r3 - r7
948 84058600 : bp = r + s
949 84058600 : bm = r - s
950 84058600 : r = s4 - s8
951 84058600 : s = r2 - r6
952 84058600 : cp = r + s
953 84058600 : cm = r - s
954 84058600 : r = s2 - s6
955 84058600 : s = r4 - r8
956 84058600 : dpp = r + s
957 84058600 : dm = r - s
958 84058600 : r = (cp + dm)*rt2i
959 84058600 : s = (dm - cp)*rt2i
960 84058600 : cp = (cm + dpp)*rt2i
961 84058600 : dpp = (cm - dpp)*rt2i
962 84058600 : zout(1, j, nout2) = ap + r
963 84058600 : zout(2, j, nout2) = bm + s
964 84058600 : zout(1, j, nout6) = ap - r
965 84058600 : zout(2, j, nout6) = bm - s
966 84058600 : zout(1, j, nout4) = am + cp
967 84058600 : zout(2, j, nout4) = bp + dpp
968 84058600 : zout(1, j, nout8) = am - cp
969 90309885 : zout(2, j, nout8) = bp - dpp
970 : END DO; END DO
971 1220866 : 8000 CONTINUE
972 :
973 : ELSE
974 1312342 : ia = 1
975 1312342 : nin1 = ia - after
976 1312342 : nout1 = ia - atn
977 11250386 : DO ib = 1, before
978 9938044 : nin1 = nin1 + after
979 9938044 : nin2 = nin1 + atb
980 9938044 : nin3 = nin2 + atb
981 9938044 : nin4 = nin3 + atb
982 9938044 : nin5 = nin4 + atb
983 9938044 : nin6 = nin5 + atb
984 9938044 : nin7 = nin6 + atb
985 9938044 : nin8 = nin7 + atb
986 9938044 : nout1 = nout1 + atn
987 9938044 : nout2 = nout1 + after
988 9938044 : nout3 = nout2 + after
989 9938044 : nout4 = nout3 + after
990 9938044 : nout5 = nout4 + after
991 9938044 : nout6 = nout5 + after
992 9938044 : nout7 = nout6 + after
993 9938044 : nout8 = nout7 + after
994 201484352 : DO j = 1, nfft
995 190233966 : r1 = zin(1, j, nin1)
996 190233966 : s1 = zin(2, j, nin1)
997 190233966 : r2 = zin(1, j, nin2)
998 190233966 : s2 = zin(2, j, nin2)
999 190233966 : r3 = zin(1, j, nin3)
1000 190233966 : s3 = zin(2, j, nin3)
1001 190233966 : r4 = zin(1, j, nin4)
1002 190233966 : s4 = zin(2, j, nin4)
1003 190233966 : r5 = zin(1, j, nin5)
1004 190233966 : s5 = zin(2, j, nin5)
1005 190233966 : r6 = zin(1, j, nin6)
1006 190233966 : s6 = zin(2, j, nin6)
1007 190233966 : r7 = zin(1, j, nin7)
1008 190233966 : s7 = zin(2, j, nin7)
1009 190233966 : r8 = zin(1, j, nin8)
1010 190233966 : s8 = zin(2, j, nin8)
1011 190233966 : r = r1 + r5
1012 190233966 : s = r3 + r7
1013 190233966 : ap = r + s
1014 190233966 : am = r - s
1015 190233966 : r = r2 + r6
1016 190233966 : s = r4 + r8
1017 190233966 : bp = r + s
1018 190233966 : bm = r - s
1019 190233966 : r = s1 + s5
1020 190233966 : s = s3 + s7
1021 190233966 : cp = r + s
1022 190233966 : cm = r - s
1023 190233966 : r = s2 + s6
1024 190233966 : s = s4 + s8
1025 190233966 : dpp = r + s
1026 190233966 : dm = r - s
1027 190233966 : zout(1, j, nout1) = ap + bp
1028 190233966 : zout(2, j, nout1) = cp + dpp
1029 190233966 : zout(1, j, nout5) = ap - bp
1030 190233966 : zout(2, j, nout5) = cp - dpp
1031 190233966 : zout(1, j, nout3) = am - dm
1032 190233966 : zout(2, j, nout3) = cm + bm
1033 190233966 : zout(1, j, nout7) = am + dm
1034 190233966 : zout(2, j, nout7) = cm - bm
1035 190233966 : r = r1 - r5
1036 190233966 : s = -s3 + s7
1037 190233966 : ap = r + s
1038 190233966 : am = r - s
1039 190233966 : r = s1 - s5
1040 190233966 : s = r7 - r3
1041 190233966 : bp = r + s
1042 190233966 : bm = r - s
1043 190233966 : r = -s4 + s8
1044 190233966 : s = r2 - r6
1045 190233966 : cp = r + s
1046 190233966 : cm = r - s
1047 190233966 : r = -s2 + s6
1048 190233966 : s = r4 - r8
1049 190233966 : dpp = r + s
1050 190233966 : dm = r - s
1051 190233966 : r = (cp + dm)*rt2i
1052 190233966 : s = (cp - dm)*rt2i
1053 190233966 : cp = (cm + dpp)*rt2i
1054 190233966 : dpp = (dpp - cm)*rt2i
1055 190233966 : zout(1, j, nout2) = ap + r
1056 190233966 : zout(2, j, nout2) = bm + s
1057 190233966 : zout(1, j, nout6) = ap - r
1058 190233966 : zout(2, j, nout6) = bm - s
1059 190233966 : zout(1, j, nout4) = am + cp
1060 190233966 : zout(2, j, nout4) = bp + dpp
1061 190233966 : zout(1, j, nout8) = am - cp
1062 200172010 : zout(2, j, nout8) = bp - dpp
1063 : END DO; END DO
1064 :
1065 3787602 : DO 8001, ia = 2, after
1066 2475260 : ias = ia - 1
1067 2475260 : itt = ias*before
1068 2475260 : itrig = itt + 1
1069 2475260 : cr2 = trig(1, itrig)
1070 2475260 : ci2 = trig(2, itrig)
1071 2475260 : itrig = itrig + itt
1072 2475260 : cr3 = trig(1, itrig)
1073 2475260 : ci3 = trig(2, itrig)
1074 2475260 : itrig = itrig + itt
1075 2475260 : cr4 = trig(1, itrig)
1076 2475260 : ci4 = trig(2, itrig)
1077 2475260 : itrig = itrig + itt
1078 2475260 : cr5 = trig(1, itrig)
1079 2475260 : ci5 = trig(2, itrig)
1080 2475260 : itrig = itrig + itt
1081 2475260 : cr6 = trig(1, itrig)
1082 2475260 : ci6 = trig(2, itrig)
1083 2475260 : itrig = itrig + itt
1084 2475260 : cr7 = trig(1, itrig)
1085 2475260 : ci7 = trig(2, itrig)
1086 2475260 : itrig = itrig + itt
1087 2475260 : cr8 = trig(1, itrig)
1088 2475260 : ci8 = trig(2, itrig)
1089 2475260 : nin1 = ia - after
1090 2475260 : nout1 = ia - atn
1091 9195811 : DO ib = 1, before
1092 6720551 : nin1 = nin1 + after
1093 6720551 : nin2 = nin1 + atb
1094 6720551 : nin3 = nin2 + atb
1095 6720551 : nin4 = nin3 + atb
1096 6720551 : nin5 = nin4 + atb
1097 6720551 : nin6 = nin5 + atb
1098 6720551 : nin7 = nin6 + atb
1099 6720551 : nin8 = nin7 + atb
1100 6720551 : nout1 = nout1 + atn
1101 6720551 : nout2 = nout1 + after
1102 6720551 : nout3 = nout2 + after
1103 6720551 : nout4 = nout3 + after
1104 6720551 : nout5 = nout4 + after
1105 6720551 : nout6 = nout5 + after
1106 6720551 : nout7 = nout6 + after
1107 6720551 : nout8 = nout7 + after
1108 98906246 : DO j = 1, nfft
1109 89710435 : r1 = zin(1, j, nin1)
1110 89710435 : s1 = zin(2, j, nin1)
1111 89710435 : r = zin(1, j, nin2)
1112 89710435 : s = zin(2, j, nin2)
1113 89710435 : r2 = r*cr2 - s*ci2
1114 89710435 : s2 = r*ci2 + s*cr2
1115 89710435 : r = zin(1, j, nin3)
1116 89710435 : s = zin(2, j, nin3)
1117 89710435 : r3 = r*cr3 - s*ci3
1118 89710435 : s3 = r*ci3 + s*cr3
1119 89710435 : r = zin(1, j, nin4)
1120 89710435 : s = zin(2, j, nin4)
1121 89710435 : r4 = r*cr4 - s*ci4
1122 89710435 : s4 = r*ci4 + s*cr4
1123 89710435 : r = zin(1, j, nin5)
1124 89710435 : s = zin(2, j, nin5)
1125 89710435 : r5 = r*cr5 - s*ci5
1126 89710435 : s5 = r*ci5 + s*cr5
1127 89710435 : r = zin(1, j, nin6)
1128 89710435 : s = zin(2, j, nin6)
1129 89710435 : r6 = r*cr6 - s*ci6
1130 89710435 : s6 = r*ci6 + s*cr6
1131 89710435 : r = zin(1, j, nin7)
1132 89710435 : s = zin(2, j, nin7)
1133 89710435 : r7 = r*cr7 - s*ci7
1134 89710435 : s7 = r*ci7 + s*cr7
1135 89710435 : r = zin(1, j, nin8)
1136 89710435 : s = zin(2, j, nin8)
1137 89710435 : r8 = r*cr8 - s*ci8
1138 89710435 : s8 = r*ci8 + s*cr8
1139 89710435 : r = r1 + r5
1140 89710435 : s = r3 + r7
1141 89710435 : ap = r + s
1142 89710435 : am = r - s
1143 89710435 : r = r2 + r6
1144 89710435 : s = r4 + r8
1145 89710435 : bp = r + s
1146 89710435 : bm = r - s
1147 89710435 : r = s1 + s5
1148 89710435 : s = s3 + s7
1149 89710435 : cp = r + s
1150 89710435 : cm = r - s
1151 89710435 : r = s2 + s6
1152 89710435 : s = s4 + s8
1153 89710435 : dpp = r + s
1154 89710435 : dm = r - s
1155 89710435 : zout(1, j, nout1) = ap + bp
1156 89710435 : zout(2, j, nout1) = cp + dpp
1157 89710435 : zout(1, j, nout5) = ap - bp
1158 89710435 : zout(2, j, nout5) = cp - dpp
1159 89710435 : zout(1, j, nout3) = am - dm
1160 89710435 : zout(2, j, nout3) = cm + bm
1161 89710435 : zout(1, j, nout7) = am + dm
1162 89710435 : zout(2, j, nout7) = cm - bm
1163 89710435 : r = r1 - r5
1164 89710435 : s = -s3 + s7
1165 89710435 : ap = r + s
1166 89710435 : am = r - s
1167 89710435 : r = s1 - s5
1168 89710435 : s = r7 - r3
1169 89710435 : bp = r + s
1170 89710435 : bm = r - s
1171 89710435 : r = -s4 + s8
1172 89710435 : s = r2 - r6
1173 89710435 : cp = r + s
1174 89710435 : cm = r - s
1175 89710435 : r = -s2 + s6
1176 89710435 : s = r4 - r8
1177 89710435 : dpp = r + s
1178 89710435 : dm = r - s
1179 89710435 : r = (cp + dm)*rt2i
1180 89710435 : s = (cp - dm)*rt2i
1181 89710435 : cp = (cm + dpp)*rt2i
1182 89710435 : dpp = (dpp - cm)*rt2i
1183 89710435 : zout(1, j, nout2) = ap + r
1184 89710435 : zout(2, j, nout2) = bm + s
1185 89710435 : zout(1, j, nout6) = ap - r
1186 89710435 : zout(2, j, nout6) = bm - s
1187 89710435 : zout(1, j, nout4) = am + cp
1188 89710435 : zout(2, j, nout4) = bp + dpp
1189 89710435 : zout(1, j, nout8) = am - cp
1190 96430986 : zout(2, j, nout8) = bp - dpp
1191 : END DO; END DO
1192 1312342 : 8001 CONTINUE
1193 :
1194 : END IF
1195 : ELSE IF (now .EQ. 3) THEN
1196 : ! .5_dp*sqrt(3._dp)
1197 8993218 : bb = isign*0.8660254037844387_dp
1198 8993218 : ia = 1
1199 8993218 : nin1 = ia - after
1200 8993218 : nout1 = ia - atn
1201 32399898 : DO ib = 1, before
1202 23406680 : nin1 = nin1 + after
1203 23406680 : nin2 = nin1 + atb
1204 23406680 : nin3 = nin2 + atb
1205 23406680 : nout1 = nout1 + atn
1206 23406680 : nout2 = nout1 + after
1207 23406680 : nout3 = nout2 + after
1208 487884536 : DO j = 1, nfft
1209 455484638 : r1 = zin(1, j, nin1)
1210 455484638 : s1 = zin(2, j, nin1)
1211 455484638 : r2 = zin(1, j, nin2)
1212 455484638 : s2 = zin(2, j, nin2)
1213 455484638 : r3 = zin(1, j, nin3)
1214 455484638 : s3 = zin(2, j, nin3)
1215 455484638 : r = r2 + r3
1216 455484638 : s = s2 + s3
1217 455484638 : zout(1, j, nout1) = r + r1
1218 455484638 : zout(2, j, nout1) = s + s1
1219 455484638 : r1 = r1 - .5_dp*r
1220 455484638 : s1 = s1 - .5_dp*s
1221 455484638 : r2 = bb*(r2 - r3)
1222 455484638 : s2 = bb*(s2 - s3)
1223 455484638 : zout(1, j, nout2) = r1 - s2
1224 455484638 : zout(2, j, nout2) = s1 + r2
1225 455484638 : zout(1, j, nout3) = r1 + s2
1226 478891318 : zout(2, j, nout3) = s1 - r2
1227 : END DO; END DO
1228 168265186 : DO 3000, ia = 2, after
1229 159271968 : ias = ia - 1
1230 159271968 : IF (4*ias .EQ. 3*after) THEN
1231 5904693 : IF (isign .EQ. 1) THEN
1232 3032915 : nin1 = ia - after
1233 3032915 : nout1 = ia - atn
1234 12179276 : DO ib = 1, before
1235 9146361 : nin1 = nin1 + after
1236 9146361 : nin2 = nin1 + atb
1237 9146361 : nin3 = nin2 + atb
1238 9146361 : nout1 = nout1 + atn
1239 9146361 : nout2 = nout1 + after
1240 9146361 : nout3 = nout2 + after
1241 186881444 : DO j = 1, nfft
1242 174702168 : r1 = zin(1, j, nin1)
1243 174702168 : s1 = zin(2, j, nin1)
1244 174702168 : r2 = zin(2, j, nin2)
1245 174702168 : s2 = zin(1, j, nin2)
1246 174702168 : r3 = zin(1, j, nin3)
1247 174702168 : s3 = zin(2, j, nin3)
1248 174702168 : r = r3 + r2
1249 174702168 : s = s2 - s3
1250 174702168 : zout(1, j, nout1) = r1 - r
1251 174702168 : zout(2, j, nout1) = s + s1
1252 174702168 : r1 = r1 + .5_dp*r
1253 174702168 : s1 = s1 - .5_dp*s
1254 174702168 : r2 = bb*(r2 - r3)
1255 174702168 : s2 = bb*(s2 + s3)
1256 174702168 : zout(1, j, nout2) = r1 - s2
1257 174702168 : zout(2, j, nout2) = s1 - r2
1258 174702168 : zout(1, j, nout3) = r1 + s2
1259 183848529 : zout(2, j, nout3) = s1 + r2
1260 : END DO; END DO
1261 : ELSE
1262 2871778 : nin1 = ia - after
1263 2871778 : nout1 = ia - atn
1264 11521230 : DO ib = 1, before
1265 8649452 : nin1 = nin1 + after
1266 8649452 : nin2 = nin1 + atb
1267 8649452 : nin3 = nin2 + atb
1268 8649452 : nout1 = nout1 + atn
1269 8649452 : nout2 = nout1 + after
1270 8649452 : nout3 = nout2 + after
1271 177353976 : DO j = 1, nfft
1272 165832746 : r1 = zin(1, j, nin1)
1273 165832746 : s1 = zin(2, j, nin1)
1274 165832746 : r2 = zin(2, j, nin2)
1275 165832746 : s2 = zin(1, j, nin2)
1276 165832746 : r3 = zin(1, j, nin3)
1277 165832746 : s3 = zin(2, j, nin3)
1278 165832746 : r = r2 - r3
1279 165832746 : s = s2 + s3
1280 165832746 : zout(1, j, nout1) = r + r1
1281 165832746 : zout(2, j, nout1) = s1 - s
1282 165832746 : r1 = r1 - .5_dp*r
1283 165832746 : s1 = s1 + .5_dp*s
1284 165832746 : r2 = bb*(r2 + r3)
1285 165832746 : s2 = bb*(s2 - s3)
1286 165832746 : zout(1, j, nout2) = r1 + s2
1287 165832746 : zout(2, j, nout2) = s1 + r2
1288 165832746 : zout(1, j, nout3) = r1 - s2
1289 174482198 : zout(2, j, nout3) = s1 - r2
1290 : END DO; END DO
1291 : END IF
1292 153367275 : ELSE IF (8*ias .EQ. 3*after) THEN
1293 2116025 : IF (isign .EQ. 1) THEN
1294 1087911 : nin1 = ia - after
1295 1087911 : nout1 = ia - atn
1296 2616044 : DO ib = 1, before
1297 1528133 : nin1 = nin1 + after
1298 1528133 : nin2 = nin1 + atb
1299 1528133 : nin3 = nin2 + atb
1300 1528133 : nout1 = nout1 + atn
1301 1528133 : nout2 = nout1 + after
1302 1528133 : nout3 = nout2 + after
1303 32523134 : DO j = 1, nfft
1304 29907090 : r1 = zin(1, j, nin1)
1305 29907090 : s1 = zin(2, j, nin1)
1306 29907090 : r = zin(1, j, nin2)
1307 29907090 : s = zin(2, j, nin2)
1308 29907090 : r2 = (r - s)*rt2i
1309 29907090 : s2 = (r + s)*rt2i
1310 29907090 : r3 = zin(2, j, nin3)
1311 29907090 : s3 = zin(1, j, nin3)
1312 29907090 : r = r2 - r3
1313 29907090 : s = s2 + s3
1314 29907090 : zout(1, j, nout1) = r + r1
1315 29907090 : zout(2, j, nout1) = s + s1
1316 29907090 : r1 = r1 - .5_dp*r
1317 29907090 : s1 = s1 - .5_dp*s
1318 29907090 : r2 = bb*(r2 + r3)
1319 29907090 : s2 = bb*(s2 - s3)
1320 29907090 : zout(1, j, nout2) = r1 - s2
1321 29907090 : zout(2, j, nout2) = s1 + r2
1322 29907090 : zout(1, j, nout3) = r1 + s2
1323 31435223 : zout(2, j, nout3) = s1 - r2
1324 : END DO; END DO
1325 : ELSE
1326 1028114 : nin1 = ia - after
1327 1028114 : nout1 = ia - atn
1328 2465598 : DO ib = 1, before
1329 1437484 : nin1 = nin1 + after
1330 1437484 : nin2 = nin1 + atb
1331 1437484 : nin3 = nin2 + atb
1332 1437484 : nout1 = nout1 + atn
1333 1437484 : nout2 = nout1 + after
1334 1437484 : nout3 = nout2 + after
1335 30450024 : DO j = 1, nfft
1336 27984426 : r1 = zin(1, j, nin1)
1337 27984426 : s1 = zin(2, j, nin1)
1338 27984426 : r = zin(1, j, nin2)
1339 27984426 : s = zin(2, j, nin2)
1340 27984426 : r2 = (r + s)*rt2i
1341 27984426 : s2 = (s - r)*rt2i
1342 27984426 : r3 = zin(2, j, nin3)
1343 27984426 : s3 = zin(1, j, nin3)
1344 27984426 : r = r2 + r3
1345 27984426 : s = s2 - s3
1346 27984426 : zout(1, j, nout1) = r + r1
1347 27984426 : zout(2, j, nout1) = s + s1
1348 27984426 : r1 = r1 - .5_dp*r
1349 27984426 : s1 = s1 - .5_dp*s
1350 27984426 : r2 = bb*(r2 - r3)
1351 27984426 : s2 = bb*(s2 + s3)
1352 27984426 : zout(1, j, nout2) = r1 - s2
1353 27984426 : zout(2, j, nout2) = s1 + r2
1354 27984426 : zout(1, j, nout3) = r1 + s2
1355 29421910 : zout(2, j, nout3) = s1 - r2
1356 : END DO; END DO
1357 : END IF
1358 : ELSE
1359 151251250 : itt = ias*before
1360 151251250 : itrig = itt + 1
1361 151251250 : cr2 = trig(1, itrig)
1362 151251250 : ci2 = trig(2, itrig)
1363 151251250 : itrig = itrig + itt
1364 151251250 : cr3 = trig(1, itrig)
1365 151251250 : ci3 = trig(2, itrig)
1366 151251250 : nin1 = ia - after
1367 151251250 : nout1 = ia - atn
1368 360670110 : DO ib = 1, before
1369 209418860 : nin1 = nin1 + after
1370 209418860 : nin2 = nin1 + atb
1371 209418860 : nin3 = nin2 + atb
1372 209418860 : nout1 = nout1 + atn
1373 209418860 : nout2 = nout1 + after
1374 209418860 : nout3 = nout2 + after
1375 4148801246 : DO j = 1, nfft
1376 3788131136 : r1 = zin(1, j, nin1)
1377 3788131136 : s1 = zin(2, j, nin1)
1378 3788131136 : r = zin(1, j, nin2)
1379 3788131136 : s = zin(2, j, nin2)
1380 3788131136 : r2 = r*cr2 - s*ci2
1381 3788131136 : s2 = r*ci2 + s*cr2
1382 3788131136 : r = zin(1, j, nin3)
1383 3788131136 : s = zin(2, j, nin3)
1384 3788131136 : r3 = r*cr3 - s*ci3
1385 3788131136 : s3 = r*ci3 + s*cr3
1386 3788131136 : r = r2 + r3
1387 3788131136 : s = s2 + s3
1388 3788131136 : zout(1, j, nout1) = r + r1
1389 3788131136 : zout(2, j, nout1) = s + s1
1390 3788131136 : r1 = r1 - .5_dp*r
1391 3788131136 : s1 = s1 - .5_dp*s
1392 3788131136 : r2 = bb*(r2 - r3)
1393 3788131136 : s2 = bb*(s2 - s3)
1394 3788131136 : zout(1, j, nout2) = r1 - s2
1395 3788131136 : zout(2, j, nout2) = s1 + r2
1396 3788131136 : zout(1, j, nout3) = r1 + s2
1397 3997549996 : zout(2, j, nout3) = s1 - r2
1398 : END DO; END DO
1399 : END IF
1400 8993218 : 3000 CONTINUE
1401 : ELSE IF (now .EQ. 5) THEN
1402 : ! cos(2._dp*pi/5._dp)
1403 3341362 : cos2 = 0.3090169943749474_dp
1404 : ! cos(4._dp*pi/5._dp)
1405 3341362 : cos4 = -0.8090169943749474_dp
1406 : ! sin(2._dp*pi/5._dp)
1407 3341362 : sin2 = isign*0.9510565162951536_dp
1408 : ! sin(4._dp*pi/5._dp)
1409 3341362 : sin4 = isign*0.5877852522924731_dp
1410 3341362 : ia = 1
1411 3341362 : nin1 = ia - after
1412 3341362 : nout1 = ia - atn
1413 34431008 : DO ib = 1, before
1414 31089646 : nin1 = nin1 + after
1415 31089646 : nin2 = nin1 + atb
1416 31089646 : nin3 = nin2 + atb
1417 31089646 : nin4 = nin3 + atb
1418 31089646 : nin5 = nin4 + atb
1419 31089646 : nout1 = nout1 + atn
1420 31089646 : nout2 = nout1 + after
1421 31089646 : nout3 = nout2 + after
1422 31089646 : nout4 = nout3 + after
1423 31089646 : nout5 = nout4 + after
1424 660417457 : DO j = 1, nfft
1425 625986449 : r1 = zin(1, j, nin1)
1426 625986449 : s1 = zin(2, j, nin1)
1427 625986449 : r2 = zin(1, j, nin2)
1428 625986449 : s2 = zin(2, j, nin2)
1429 625986449 : r3 = zin(1, j, nin3)
1430 625986449 : s3 = zin(2, j, nin3)
1431 625986449 : r4 = zin(1, j, nin4)
1432 625986449 : s4 = zin(2, j, nin4)
1433 625986449 : r5 = zin(1, j, nin5)
1434 625986449 : s5 = zin(2, j, nin5)
1435 625986449 : r25 = r2 + r5
1436 625986449 : r34 = r3 + r4
1437 625986449 : s25 = s2 - s5
1438 625986449 : s34 = s3 - s4
1439 625986449 : zout(1, j, nout1) = r1 + r25 + r34
1440 625986449 : r = r1 + cos2*r25 + cos4*r34
1441 625986449 : s = sin2*s25 + sin4*s34
1442 625986449 : zout(1, j, nout2) = r - s
1443 625986449 : zout(1, j, nout5) = r + s
1444 625986449 : r = r1 + cos4*r25 + cos2*r34
1445 625986449 : s = sin4*s25 - sin2*s34
1446 625986449 : zout(1, j, nout3) = r - s
1447 625986449 : zout(1, j, nout4) = r + s
1448 625986449 : r25 = r2 - r5
1449 625986449 : r34 = r3 - r4
1450 625986449 : s25 = s2 + s5
1451 625986449 : s34 = s3 + s4
1452 625986449 : zout(2, j, nout1) = s1 + s25 + s34
1453 625986449 : r = s1 + cos2*s25 + cos4*s34
1454 625986449 : s = sin2*r25 + sin4*r34
1455 625986449 : zout(2, j, nout2) = r + s
1456 625986449 : zout(2, j, nout5) = r - s
1457 625986449 : r = s1 + cos4*s25 + cos2*s34
1458 625986449 : s = sin4*r25 - sin2*r34
1459 625986449 : zout(2, j, nout3) = r + s
1460 657076095 : zout(2, j, nout4) = r - s
1461 : END DO; END DO
1462 11552910 : DO 5000, ia = 2, after
1463 8211548 : ias = ia - 1
1464 8211548 : IF (8*ias .EQ. 5*after) THEN
1465 621304 : IF (isign .EQ. 1) THEN
1466 324186 : nin1 = ia - after
1467 324186 : nout1 = ia - atn
1468 948492 : DO ib = 1, before
1469 624306 : nin1 = nin1 + after
1470 624306 : nin2 = nin1 + atb
1471 624306 : nin3 = nin2 + atb
1472 624306 : nin4 = nin3 + atb
1473 624306 : nin5 = nin4 + atb
1474 624306 : nout1 = nout1 + atn
1475 624306 : nout2 = nout1 + after
1476 624306 : nout3 = nout2 + after
1477 624306 : nout4 = nout3 + after
1478 624306 : nout5 = nout4 + after
1479 14667879 : DO j = 1, nfft
1480 13719387 : r1 = zin(1, j, nin1)
1481 13719387 : s1 = zin(2, j, nin1)
1482 13719387 : r = zin(1, j, nin2)
1483 13719387 : s = zin(2, j, nin2)
1484 13719387 : r2 = (r - s)*rt2i
1485 13719387 : s2 = (r + s)*rt2i
1486 13719387 : r3 = zin(2, j, nin3)
1487 13719387 : s3 = zin(1, j, nin3)
1488 13719387 : r = zin(1, j, nin4)
1489 13719387 : s = zin(2, j, nin4)
1490 13719387 : r4 = (r + s)*rt2i
1491 13719387 : s4 = (r - s)*rt2i
1492 13719387 : r5 = zin(1, j, nin5)
1493 13719387 : s5 = zin(2, j, nin5)
1494 13719387 : r25 = r2 - r5
1495 13719387 : r34 = r3 + r4
1496 13719387 : s25 = s2 + s5
1497 13719387 : s34 = s3 - s4
1498 13719387 : zout(1, j, nout1) = r1 + r25 - r34
1499 13719387 : r = r1 + cos2*r25 - cos4*r34
1500 13719387 : s = sin2*s25 + sin4*s34
1501 13719387 : zout(1, j, nout2) = r - s
1502 13719387 : zout(1, j, nout5) = r + s
1503 13719387 : r = r1 + cos4*r25 - cos2*r34
1504 13719387 : s = sin4*s25 - sin2*s34
1505 13719387 : zout(1, j, nout3) = r - s
1506 13719387 : zout(1, j, nout4) = r + s
1507 13719387 : r25 = r2 + r5
1508 13719387 : r34 = r4 - r3
1509 13719387 : s25 = s2 - s5
1510 13719387 : s34 = s3 + s4
1511 13719387 : zout(2, j, nout1) = s1 + s25 + s34
1512 13719387 : r = s1 + cos2*s25 + cos4*s34
1513 13719387 : s = sin2*r25 + sin4*r34
1514 13719387 : zout(2, j, nout2) = r + s
1515 13719387 : zout(2, j, nout5) = r - s
1516 13719387 : r = s1 + cos4*s25 + cos2*s34
1517 13719387 : s = sin4*r25 - sin2*r34
1518 13719387 : zout(2, j, nout3) = r + s
1519 14343693 : zout(2, j, nout4) = r - s
1520 : END DO; END DO
1521 : ELSE
1522 297118 : nin1 = ia - after
1523 297118 : nout1 = ia - atn
1524 881180 : DO ib = 1, before
1525 584062 : nin1 = nin1 + after
1526 584062 : nin2 = nin1 + atb
1527 584062 : nin3 = nin2 + atb
1528 584062 : nin4 = nin3 + atb
1529 584062 : nin5 = nin4 + atb
1530 584062 : nout1 = nout1 + atn
1531 584062 : nout2 = nout1 + after
1532 584062 : nout3 = nout2 + after
1533 584062 : nout4 = nout3 + after
1534 584062 : nout5 = nout4 + after
1535 13483260 : DO j = 1, nfft
1536 12602080 : r1 = zin(1, j, nin1)
1537 12602080 : s1 = zin(2, j, nin1)
1538 12602080 : r = zin(1, j, nin2)
1539 12602080 : s = zin(2, j, nin2)
1540 12602080 : r2 = (r + s)*rt2i
1541 12602080 : s2 = (s - r)*rt2i
1542 12602080 : r3 = zin(2, j, nin3)
1543 12602080 : s3 = zin(1, j, nin3)
1544 12602080 : r = zin(1, j, nin4)
1545 12602080 : s = zin(2, j, nin4)
1546 12602080 : r4 = (s - r)*rt2i
1547 12602080 : s4 = (r + s)*rt2i
1548 12602080 : r5 = zin(1, j, nin5)
1549 12602080 : s5 = zin(2, j, nin5)
1550 12602080 : r25 = r2 - r5
1551 12602080 : r34 = r3 + r4
1552 12602080 : s25 = s2 + s5
1553 12602080 : s34 = s4 - s3
1554 12602080 : zout(1, j, nout1) = r1 + r25 + r34
1555 12602080 : r = r1 + cos2*r25 + cos4*r34
1556 12602080 : s = sin2*s25 + sin4*s34
1557 12602080 : zout(1, j, nout2) = r - s
1558 12602080 : zout(1, j, nout5) = r + s
1559 12602080 : r = r1 + cos4*r25 + cos2*r34
1560 12602080 : s = sin4*s25 - sin2*s34
1561 12602080 : zout(1, j, nout3) = r - s
1562 12602080 : zout(1, j, nout4) = r + s
1563 12602080 : r25 = r2 + r5
1564 12602080 : r34 = r3 - r4
1565 12602080 : s25 = s2 - s5
1566 12602080 : s34 = s3 + s4
1567 12602080 : zout(2, j, nout1) = s1 + s25 - s34
1568 12602080 : r = s1 + cos2*s25 - cos4*s34
1569 12602080 : s = sin2*r25 + sin4*r34
1570 12602080 : zout(2, j, nout2) = r + s
1571 12602080 : zout(2, j, nout5) = r - s
1572 12602080 : r = s1 + cos4*s25 - cos2*s34
1573 12602080 : s = sin4*r25 - sin2*r34
1574 12602080 : zout(2, j, nout3) = r + s
1575 13186142 : zout(2, j, nout4) = r - s
1576 : END DO; END DO
1577 : END IF
1578 : ELSE
1579 7590244 : ias = ia - 1
1580 7590244 : itt = ias*before
1581 7590244 : itrig = itt + 1
1582 7590244 : cr2 = trig(1, itrig)
1583 7590244 : ci2 = trig(2, itrig)
1584 7590244 : itrig = itrig + itt
1585 7590244 : cr3 = trig(1, itrig)
1586 7590244 : ci3 = trig(2, itrig)
1587 7590244 : itrig = itrig + itt
1588 7590244 : cr4 = trig(1, itrig)
1589 7590244 : ci4 = trig(2, itrig)
1590 7590244 : itrig = itrig + itt
1591 7590244 : cr5 = trig(1, itrig)
1592 7590244 : ci5 = trig(2, itrig)
1593 7590244 : nin1 = ia - after
1594 7590244 : nout1 = ia - atn
1595 24453112 : DO ib = 1, before
1596 16862868 : nin1 = nin1 + after
1597 16862868 : nin2 = nin1 + atb
1598 16862868 : nin3 = nin2 + atb
1599 16862868 : nin4 = nin3 + atb
1600 16862868 : nin5 = nin4 + atb
1601 16862868 : nout1 = nout1 + atn
1602 16862868 : nout2 = nout1 + after
1603 16862868 : nout3 = nout2 + after
1604 16862868 : nout4 = nout3 + after
1605 16862868 : nout5 = nout4 + after
1606 354440434 : DO j = 1, nfft
1607 329987322 : r1 = zin(1, j, nin1)
1608 329987322 : s1 = zin(2, j, nin1)
1609 329987322 : r = zin(1, j, nin2)
1610 329987322 : s = zin(2, j, nin2)
1611 329987322 : r2 = r*cr2 - s*ci2
1612 329987322 : s2 = r*ci2 + s*cr2
1613 329987322 : r = zin(1, j, nin3)
1614 329987322 : s = zin(2, j, nin3)
1615 329987322 : r3 = r*cr3 - s*ci3
1616 329987322 : s3 = r*ci3 + s*cr3
1617 329987322 : r = zin(1, j, nin4)
1618 329987322 : s = zin(2, j, nin4)
1619 329987322 : r4 = r*cr4 - s*ci4
1620 329987322 : s4 = r*ci4 + s*cr4
1621 329987322 : r = zin(1, j, nin5)
1622 329987322 : s = zin(2, j, nin5)
1623 329987322 : r5 = r*cr5 - s*ci5
1624 329987322 : s5 = r*ci5 + s*cr5
1625 329987322 : r25 = r2 + r5
1626 329987322 : r34 = r3 + r4
1627 329987322 : s25 = s2 - s5
1628 329987322 : s34 = s3 - s4
1629 329987322 : zout(1, j, nout1) = r1 + r25 + r34
1630 329987322 : r = r1 + cos2*r25 + cos4*r34
1631 329987322 : s = sin2*s25 + sin4*s34
1632 329987322 : zout(1, j, nout2) = r - s
1633 329987322 : zout(1, j, nout5) = r + s
1634 329987322 : r = r1 + cos4*r25 + cos2*r34
1635 329987322 : s = sin4*s25 - sin2*s34
1636 329987322 : zout(1, j, nout3) = r - s
1637 329987322 : zout(1, j, nout4) = r + s
1638 329987322 : r25 = r2 - r5
1639 329987322 : r34 = r3 - r4
1640 329987322 : s25 = s2 + s5
1641 329987322 : s34 = s3 + s4
1642 329987322 : zout(2, j, nout1) = s1 + s25 + s34
1643 329987322 : r = s1 + cos2*s25 + cos4*s34
1644 329987322 : s = sin2*r25 + sin4*r34
1645 329987322 : zout(2, j, nout2) = r + s
1646 329987322 : zout(2, j, nout5) = r - s
1647 329987322 : r = s1 + cos4*s25 + cos2*s34
1648 329987322 : s = sin4*r25 - sin2*r34
1649 329987322 : zout(2, j, nout3) = r + s
1650 346850190 : zout(2, j, nout4) = r - s
1651 : END DO; END DO
1652 : END IF
1653 3341362 : 5000 CONTINUE
1654 : ELSE IF (now .EQ. 6) THEN
1655 : ! .5_dp*sqrt(3._dp)
1656 2852361 : bb = isign*0.8660254037844387_dp
1657 :
1658 2852361 : ia = 1
1659 2852361 : nin1 = ia - after
1660 2852361 : nout1 = ia - atn
1661 39449688 : DO ib = 1, before
1662 36597327 : nin1 = nin1 + after
1663 36597327 : nin2 = nin1 + atb
1664 36597327 : nin3 = nin2 + atb
1665 36597327 : nin4 = nin3 + atb
1666 36597327 : nin5 = nin4 + atb
1667 36597327 : nin6 = nin5 + atb
1668 36597327 : nout1 = nout1 + atn
1669 36597327 : nout2 = nout1 + after
1670 36597327 : nout3 = nout2 + after
1671 36597327 : nout4 = nout3 + after
1672 36597327 : nout5 = nout4 + after
1673 36597327 : nout6 = nout5 + after
1674 643144404 : DO j = 1, nfft
1675 603694716 : r2 = zin(1, j, nin3)
1676 603694716 : s2 = zin(2, j, nin3)
1677 603694716 : r3 = zin(1, j, nin5)
1678 603694716 : s3 = zin(2, j, nin5)
1679 603694716 : r = r2 + r3
1680 603694716 : s = s2 + s3
1681 603694716 : r1 = zin(1, j, nin1)
1682 603694716 : s1 = zin(2, j, nin1)
1683 603694716 : ur1 = r + r1
1684 603694716 : ui1 = s + s1
1685 603694716 : r1 = r1 - .5_dp*r
1686 603694716 : s1 = s1 - .5_dp*s
1687 603694716 : r = r2 - r3
1688 603694716 : s = s2 - s3
1689 603694716 : ur2 = r1 - s*bb
1690 603694716 : ui2 = s1 + r*bb
1691 603694716 : ur3 = r1 + s*bb
1692 603694716 : ui3 = s1 - r*bb
1693 :
1694 603694716 : r2 = zin(1, j, nin6)
1695 603694716 : s2 = zin(2, j, nin6)
1696 603694716 : r3 = zin(1, j, nin2)
1697 603694716 : s3 = zin(2, j, nin2)
1698 603694716 : r = r2 + r3
1699 603694716 : s = s2 + s3
1700 603694716 : r1 = zin(1, j, nin4)
1701 603694716 : s1 = zin(2, j, nin4)
1702 603694716 : vr1 = r + r1
1703 603694716 : vi1 = s + s1
1704 603694716 : r1 = r1 - .5_dp*r
1705 603694716 : s1 = s1 - .5_dp*s
1706 603694716 : r = r2 - r3
1707 603694716 : s = s2 - s3
1708 603694716 : vr2 = r1 - s*bb
1709 603694716 : vi2 = s1 + r*bb
1710 603694716 : vr3 = r1 + s*bb
1711 603694716 : vi3 = s1 - r*bb
1712 :
1713 603694716 : zout(1, j, nout1) = ur1 + vr1
1714 603694716 : zout(2, j, nout1) = ui1 + vi1
1715 603694716 : zout(1, j, nout5) = ur2 + vr2
1716 603694716 : zout(2, j, nout5) = ui2 + vi2
1717 603694716 : zout(1, j, nout3) = ur3 + vr3
1718 603694716 : zout(2, j, nout3) = ui3 + vi3
1719 603694716 : zout(1, j, nout4) = ur1 - vr1
1720 603694716 : zout(2, j, nout4) = ui1 - vi1
1721 603694716 : zout(1, j, nout2) = ur2 - vr2
1722 603694716 : zout(2, j, nout2) = ui2 - vi2
1723 603694716 : zout(1, j, nout6) = ur3 - vr3
1724 640292043 : zout(2, j, nout6) = ui3 - vi3
1725 : END DO; END DO
1726 : ELSE
1727 0 : CPABORT("error fftstp")
1728 : END IF
1729 :
1730 23211676 : END SUBROUTINE fftstp
1731 :
1732 : END MODULE ps_wavelet_fft3d
|
