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