Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : MODULE mltfftsg_tools
10 : USE ISO_C_BINDING, ONLY: C_F_POINTER,&
11 : C_LOC
12 : USE fft_kinds, ONLY: dp
13 :
14 : !$ USE OMP_LIB, ONLY: omp_get_num_threads, omp_get_thread_num
15 :
16 : #include "../../base/base_uses.f90"
17 :
18 : PRIVATE
19 : INTEGER, PARAMETER :: ctrig_length = 1024
20 : INTEGER, PARAMETER :: cache_size = 2048
21 : PUBLIC :: mltfftsg
22 :
23 : CONTAINS
24 :
25 : ! **************************************************************************************************
26 : !> \brief ...
27 : !> \param transa ...
28 : !> \param transb ...
29 : !> \param a ...
30 : !> \param ldax ...
31 : !> \param lday ...
32 : !> \param b ...
33 : !> \param ldbx ...
34 : !> \param ldby ...
35 : !> \param n ...
36 : !> \param m ...
37 : !> \param isign ...
38 : !> \param scale ...
39 : ! **************************************************************************************************
40 17820 : SUBROUTINE mltfftsg(transa, transb, a, ldax, lday, b, ldbx, ldby, n, m, isign, scale)
41 :
42 : CHARACTER(LEN=1), INTENT(IN) :: transa, transb
43 : INTEGER, INTENT(IN) :: ldax, lday
44 : COMPLEX(dp), INTENT(INOUT) :: a(ldax, lday)
45 : INTEGER, INTENT(IN) :: ldbx, ldby
46 : COMPLEX(dp), INTENT(INOUT) :: b(ldbx, ldby)
47 : INTEGER, INTENT(IN) :: n, m, isign
48 : REAL(dp), INTENT(IN) :: scale
49 :
50 17820 : COMPLEX(dp), ALLOCATABLE, DIMENSION(:, :, :) :: z
51 : INTEGER :: after(20), before(20), chunk, i, ic, id, &
52 : iend, inzee, isig, istart, iterations, &
53 : itr, length, lot, nfft, now(20), &
54 : num_threads
55 : LOGICAL :: tscal
56 : REAL(dp) :: trig(2, 1024)
57 :
58 : ! Variables
59 :
60 17820 : LENGTH = 2*(cache_size/4 + 1)
61 :
62 17820 : ISIG = -ISIGN
63 17820 : TSCAL = (ABS(SCALE - 1._dp) > 1.e-12_dp)
64 17820 : CALL ctrig(N, TRIG, AFTER, BEFORE, NOW, ISIG, IC)
65 17820 : LOT = cache_size/(4*N)
66 17820 : LOT = LOT - MOD(LOT + 1, 2)
67 17820 : LOT = MAX(1, LOT)
68 :
69 : ! initializations for serial mode
70 17820 : id = 0; num_threads = 1
71 :
72 : !$OMP PARALLEL &
73 : !$OMP PRIVATE ( id, istart, iend, nfft, i, inzee, itr) DEFAULT(NONE) &
74 : !$OMP SHARED (NUM_THREADS,z,iterations,chunk,LOT,length,m,transa,isig, &
75 17820 : !$OMP before,after,now,trig,A,n,ldax,tscal,scale,ic,transb,ldbx,b)
76 :
77 : !$OMP SINGLE
78 : !$ num_threads = omp_get_num_threads()
79 : ALLOCATE (Z(LENGTH, 2, 0:num_threads - 1))
80 : iterations = (M + LOT - 1)/LOT
81 : chunk = LOT*((iterations + num_threads - 1)/num_threads)
82 : !$OMP END SINGLE
83 : !$OMP BARRIER
84 :
85 : !$ id = omp_get_thread_num()
86 : istart = id*chunk + 1
87 : iend = MIN((id + 1)*chunk, M)
88 :
89 : DO ITR = istart, iend, LOT
90 :
91 : NFFT = MIN(M - ITR + 1, LOT)
92 : IF (TRANSA == 'N' .OR. TRANSA == 'n') THEN
93 : CALL fftpre_cmplx(NFFT, NFFT, LDAX, LOT, N, A(1, ITR), Z(1, 1, id), &
94 : TRIG, NOW(1), AFTER(1), BEFORE(1), ISIG)
95 : ELSE
96 : CALL fftstp_cmplx(LDAX, NFFT, N, LOT, N, A(ITR, 1), Z(1, 1, id), &
97 : TRIG, NOW(1), AFTER(1), BEFORE(1), ISIG)
98 : END IF
99 : IF (TSCAL) THEN
100 : IF (LOT == NFFT) THEN
101 : CALL scaled(2*LOT*N, SCALE, Z(1, 1, id))
102 : ELSE
103 : DO I = 1, N
104 : CALL scaled(2*NFFT, SCALE, Z(LOT*(I - 1) + 1, 1, id))
105 : END DO
106 : END IF
107 : END IF
108 : IF (IC .EQ. 1) THEN
109 : IF (TRANSB == 'N' .OR. TRANSB == 'n') THEN
110 : CALL zgetmo(Z(1, 1, id), LOT, NFFT, N, B(1, ITR), LDBX)
111 : ELSE
112 : CALL matmov(NFFT, N, Z(1, 1, id), LOT, B(ITR, 1), LDBX)
113 : END IF
114 : ELSE
115 : INZEE = 1
116 : DO I = 2, IC - 1
117 : CALL fftstp_cmplx(LOT, NFFT, N, LOT, N, Z(1, INZEE, id), &
118 : Z(1, 3 - INZEE, id), TRIG, NOW(I), AFTER(I), &
119 : BEFORE(I), ISIG)
120 : INZEE = 3 - INZEE
121 : END DO
122 : IF (TRANSB == 'N' .OR. TRANSB == 'n') THEN
123 : CALL fftrot_cmplx(LOT, NFFT, N, NFFT, LDBX, Z(1, INZEE, id), &
124 : B(1, ITR), TRIG, NOW(IC), AFTER(IC), BEFORE(IC), ISIG)
125 : ELSE
126 : CALL fftstp_cmplx(LOT, NFFT, N, LDBX, N, Z(1, INZEE, id), &
127 : B(ITR, 1), TRIG, NOW(IC), AFTER(IC), BEFORE(IC), ISIG)
128 : END IF
129 : END IF
130 : END DO
131 :
132 : !$OMP END PARALLEL
133 :
134 17820 : DEALLOCATE (Z)
135 :
136 17820 : IF (TRANSB == 'N' .OR. TRANSB == 'n') THEN
137 7452 : B(1:LDBX, M + 1:LDBY) = CMPLX(0._dp, 0._dp, dp)
138 2749124 : B(N + 1:LDBX, 1:M) = CMPLX(0._dp, 0._dp, dp)
139 : ELSE
140 10368 : B(1:LDBX, N + 1:LDBY) = CMPLX(0._dp, 0._dp, dp)
141 234750 : B(M + 1:LDBX, 1:N) = CMPLX(0._dp, 0._dp, dp)
142 : END IF
143 :
144 17820 : END SUBROUTINE mltfftsg
145 :
146 : ! this formalizes what we have been assuming before, call with a complex(*) array, and passing to a real(2,*)
147 : ! **************************************************************************************************
148 : !> \brief ...
149 : !> \param mm ...
150 : !> \param nfft ...
151 : !> \param m ...
152 : !> \param nn ...
153 : !> \param n ...
154 : !> \param zin ...
155 : !> \param zout ...
156 : !> \param trig ...
157 : !> \param now ...
158 : !> \param after ...
159 : !> \param before ...
160 : !> \param isign ...
161 : ! **************************************************************************************************
162 1244146 : SUBROUTINE fftstp_cmplx(mm, nfft, m, nn, n, zin, zout, trig, now, after, before, isign)
163 :
164 : INTEGER, INTENT(IN) :: mm, nfft, m, nn, n
165 : COMPLEX(dp), DIMENSION(mm, m), INTENT(IN), TARGET :: zin
166 : COMPLEX(dp), DIMENSION(nn, n), INTENT(INOUT), &
167 : TARGET :: zout
168 : REAL(dp), DIMENSION(2, ctrig_length), INTENT(IN) :: trig
169 : INTEGER, INTENT(IN) :: now, after, before, isign
170 :
171 1244146 : REAL(dp), DIMENSION(:, :, :), POINTER :: zin_real, zout_real
172 :
173 4976584 : CALL C_F_POINTER(C_LOC(zin), zin_real, (/2, mm, m/))
174 4976584 : CALL C_F_POINTER(C_LOC(zout), zout_real, (/2, nn, n/))
175 1244146 : CALL fftstp(mm, nfft, m, nn, n, zin_real, zout_real, trig, now, after, before, isign)
176 :
177 1244146 : END SUBROUTINE
178 :
179 : ! **************************************************************************************************
180 : !> \brief ...
181 : !> \param mm ...
182 : !> \param nfft ...
183 : !> \param m ...
184 : !> \param nn ...
185 : !> \param n ...
186 : !> \param zin ...
187 : !> \param zout ...
188 : !> \param trig ...
189 : !> \param now ...
190 : !> \param after ...
191 : !> \param before ...
192 : !> \param isign ...
193 : ! **************************************************************************************************
194 424773 : SUBROUTINE fftpre_cmplx(mm, nfft, m, nn, n, zin, zout, trig, now, after, before, isign)
195 :
196 : INTEGER, INTENT(IN) :: mm, nfft, m, nn, n
197 : COMPLEX(dp), DIMENSION(m, mm), INTENT(IN), TARGET :: zin
198 : COMPLEX(dp), DIMENSION(nn, n), INTENT(INOUT), &
199 : TARGET :: zout
200 : REAL(dp), DIMENSION(2, ctrig_length), INTENT(IN) :: trig
201 : INTEGER, INTENT(IN) :: now, after, before, isign
202 :
203 424773 : REAL(dp), DIMENSION(:, :, :), POINTER :: zin_real, zout_real
204 :
205 1699092 : CALL C_F_POINTER(C_LOC(zin), zin_real, (/2, mm, m/))
206 1699092 : CALL C_F_POINTER(C_LOC(zout), zout_real, (/2, nn, n/))
207 :
208 424773 : CALL fftpre(mm, nfft, m, nn, n, zin_real, zout_real, trig, now, after, before, isign)
209 :
210 424773 : END SUBROUTINE
211 :
212 : ! **************************************************************************************************
213 : !> \brief ...
214 : !> \param mm ...
215 : !> \param nfft ...
216 : !> \param m ...
217 : !> \param nn ...
218 : !> \param n ...
219 : !> \param zin ...
220 : !> \param zout ...
221 : !> \param trig ...
222 : !> \param now ...
223 : !> \param after ...
224 : !> \param before ...
225 : !> \param isign ...
226 : ! **************************************************************************************************
227 298145 : SUBROUTINE fftrot_cmplx(mm, nfft, m, nn, n, zin, zout, trig, now, after, before, isign)
228 :
229 : USE fft_kinds, ONLY: dp
230 : INTEGER, INTENT(IN) :: mm, nfft, m, nn, n
231 : COMPLEX(dp), DIMENSION(mm, m), INTENT(IN), TARGET :: zin
232 : COMPLEX(dp), DIMENSION(n, nn), INTENT(INOUT), &
233 : TARGET :: zout
234 : REAL(dp), DIMENSION(2, ctrig_length), INTENT(IN) :: trig
235 : INTEGER, INTENT(IN) :: now, after, before, isign
236 :
237 298145 : REAL(dp), DIMENSION(:, :, :), POINTER :: zin_real, zout_real
238 :
239 1192580 : CALL C_F_POINTER(C_LOC(zin), zin_real, (/2, mm, m/))
240 1192580 : CALL C_F_POINTER(C_LOC(zout), zout_real, (/2, nn, n/))
241 :
242 298145 : CALL fftrot(mm, nfft, m, nn, n, zin_real, zout_real, trig, now, after, before, isign)
243 :
244 298145 : END SUBROUTINE
245 :
246 : !-----------------------------------------------------------------------------!
247 : ! Copyright (C) Stefan Goedecker, Lausanne, Switzerland, August 1, 1991
248 : ! Copyright (C) Stefan Goedecker, Cornell University, Ithaca, USA, 1994
249 : ! Copyright (C) Stefan Goedecker, MPI Stuttgart, Germany, 1999
250 : ! This file is distributed under the terms of the
251 : ! GNU General Public License version 2 (or later),
252 : ! see http://www.gnu.org/copyleft/gpl.txt .
253 : !-----------------------------------------------------------------------------!
254 : ! S. Goedecker: Rotating a three-dimensional array in optimal
255 : ! positions for vector processing: Case study for a three-dimensional Fast
256 : ! Fourier Transform, Comp. Phys. Commun. 76, 294 (1993)
257 : ! **************************************************************************************************
258 : !> \brief ...
259 : !> \param mm ...
260 : !> \param nfft ...
261 : !> \param m ...
262 : !> \param nn ...
263 : !> \param n ...
264 : !> \param zin ...
265 : !> \param zout ...
266 : !> \param trig ...
267 : !> \param now ...
268 : !> \param after ...
269 : !> \param before ...
270 : !> \param isign ...
271 : ! **************************************************************************************************
272 298145 : SUBROUTINE fftrot(mm, nfft, m, nn, n, zin, zout, trig, now, after, before, isign)
273 :
274 : USE fft_kinds, ONLY: dp
275 : INTEGER, INTENT(IN) :: mm, nfft, m, nn, n
276 : REAL(dp), DIMENSION(2, mm, m), INTENT(IN) :: zin
277 : REAL(dp), DIMENSION(2, n, nn), INTENT(INOUT) :: zout
278 : REAL(dp), DIMENSION(2, ctrig_length), INTENT(IN) :: trig
279 : INTEGER, INTENT(IN) :: now, after, before, isign
280 :
281 : REAL(dp), PARAMETER :: bb = 0.8660254037844387_dp, cos2 = 0.3090169943749474_dp, &
282 : cos4 = -0.8090169943749474_dp, rt2i = 0.7071067811865475_dp, &
283 : sin2p = 0.9510565162951536_dp, sin4p = 0.5877852522924731_dp
284 :
285 : INTEGER :: atb, atn, ia, ias, ib, itrig, itt, j, &
286 : nin1, nin2, nin3, nin4, nin5, nin6, &
287 : nin7, nin8, nout1, nout2, nout3, &
288 : nout4, nout5, nout6, nout7, nout8
289 : REAL(dp) :: am, ap, bbs, bm, bp, ci2, ci3, ci4, ci5, cm, cp, cr2, cr3, cr4, cr5, dbl, dm, r, &
290 : r1, r2, r25, r3, r34, r4, r5, r6, r7, r8, s, s1, s2, s25, s3, s34, s4, s5, s6, s7, s8, &
291 : sin2, sin4, ui1, ui2, ui3, ur1, ur2, ur3, vi1, vi2, vi3, vr1, vr2, vr3
292 :
293 : ! sqrt(0.5)
294 : ! sqrt(3)/2
295 : ! cos(2*pi/5)
296 : ! cos(4*pi/5)
297 : ! sin(2*pi/5)
298 : ! sin(4*pi/5)
299 : !-----------------------------------------------------------------------------!
300 :
301 298145 : atn = after*now
302 298145 : atb = after*before
303 :
304 298145 : IF (now == 4) THEN
305 49350 : IF (isign == 1) 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 : nin3 = nin2 + atb
313 0 : nin4 = nin3 + atb
314 0 : nout1 = nout1 + atn
315 0 : nout2 = nout1 + after
316 0 : nout3 = nout2 + after
317 0 : nout4 = nout3 + after
318 0 : DO j = 1, nfft
319 0 : r1 = zin(1, j, nin1)
320 0 : s1 = zin(2, j, nin1)
321 0 : r2 = zin(1, j, nin2)
322 0 : s2 = zin(2, j, nin2)
323 0 : r3 = zin(1, j, nin3)
324 0 : s3 = zin(2, j, nin3)
325 0 : r4 = zin(1, j, nin4)
326 0 : s4 = zin(2, j, nin4)
327 0 : r = r1 + r3
328 0 : s = r2 + r4
329 0 : zout(1, nout1, j) = r + s
330 0 : zout(1, nout3, j) = r - s
331 0 : r = r1 - r3
332 0 : s = s2 - s4
333 0 : zout(1, nout2, j) = r - s
334 0 : zout(1, nout4, j) = r + s
335 0 : r = s1 + s3
336 0 : s = s2 + s4
337 0 : zout(2, nout1, j) = r + s
338 0 : zout(2, nout3, j) = r - s
339 0 : r = s1 - s3
340 0 : s = r2 - r4
341 0 : zout(2, nout2, j) = r + s
342 0 : zout(2, nout4, j) = r - s
343 : END DO
344 : END DO
345 0 : DO ia = 2, after
346 0 : ias = ia - 1
347 0 : IF (2*ias == after) THEN
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 : nin3 = nin2 + atb
354 0 : nin4 = nin3 + atb
355 0 : nout1 = nout1 + atn
356 0 : nout2 = nout1 + after
357 0 : nout3 = nout2 + after
358 0 : nout4 = nout3 + after
359 0 : DO j = 1, nfft
360 0 : r1 = zin(1, j, nin1)
361 0 : s1 = zin(2, j, nin1)
362 0 : r = zin(1, j, nin2)
363 0 : s = zin(2, j, nin2)
364 0 : r2 = (r - s)*rt2i
365 0 : s2 = (r + s)*rt2i
366 0 : r3 = -zin(2, j, nin3)
367 0 : s3 = zin(1, j, nin3)
368 0 : r = zin(1, j, nin4)
369 0 : s = zin(2, j, nin4)
370 0 : r4 = -(r + s)*rt2i
371 0 : s4 = (r - s)*rt2i
372 0 : r = r1 + r3
373 0 : s = r2 + r4
374 0 : zout(1, nout1, j) = r + s
375 0 : zout(1, nout3, j) = r - s
376 0 : r = r1 - r3
377 0 : s = s2 - s4
378 0 : zout(1, nout2, j) = r - s
379 0 : zout(1, nout4, j) = r + s
380 0 : r = s1 + s3
381 0 : s = s2 + s4
382 0 : zout(2, nout1, j) = r + s
383 0 : zout(2, nout3, j) = r - s
384 0 : r = s1 - s3
385 0 : s = r2 - r4
386 0 : zout(2, nout2, j) = r + s
387 0 : zout(2, nout4, j) = r - s
388 : END DO
389 : END DO
390 : ELSE
391 0 : itt = ias*before
392 0 : itrig = itt + 1
393 0 : cr2 = trig(1, itrig)
394 0 : ci2 = trig(2, itrig)
395 0 : itrig = itrig + itt
396 0 : cr3 = trig(1, itrig)
397 0 : ci3 = trig(2, itrig)
398 0 : itrig = itrig + itt
399 0 : cr4 = trig(1, itrig)
400 0 : ci4 = trig(2, itrig)
401 0 : nin1 = ia - after
402 0 : nout1 = ia - atn
403 0 : DO ib = 1, before
404 0 : nin1 = nin1 + after
405 0 : nin2 = nin1 + atb
406 0 : nin3 = nin2 + atb
407 0 : nin4 = nin3 + atb
408 0 : nout1 = nout1 + atn
409 0 : nout2 = nout1 + after
410 0 : nout3 = nout2 + after
411 0 : nout4 = nout3 + after
412 0 : DO j = 1, nfft
413 0 : r1 = zin(1, j, nin1)
414 0 : s1 = zin(2, j, nin1)
415 0 : r = zin(1, j, nin2)
416 0 : s = zin(2, j, nin2)
417 0 : r2 = r*cr2 - s*ci2
418 0 : s2 = r*ci2 + s*cr2
419 0 : r = zin(1, j, nin3)
420 0 : s = zin(2, j, nin3)
421 0 : r3 = r*cr3 - s*ci3
422 0 : s3 = r*ci3 + s*cr3
423 0 : r = zin(1, j, nin4)
424 0 : s = zin(2, j, nin4)
425 0 : r4 = r*cr4 - s*ci4
426 0 : s4 = r*ci4 + s*cr4
427 0 : r = r1 + r3
428 0 : s = r2 + r4
429 0 : zout(1, nout1, j) = r + s
430 0 : zout(1, nout3, j) = r - s
431 0 : r = r1 - r3
432 0 : s = s2 - s4
433 0 : zout(1, nout2, j) = r - s
434 0 : zout(1, nout4, j) = r + s
435 0 : r = s1 + s3
436 0 : s = s2 + s4
437 0 : zout(2, nout1, j) = r + s
438 0 : zout(2, nout3, j) = r - s
439 0 : r = s1 - s3
440 0 : s = r2 - r4
441 0 : zout(2, nout2, j) = r + s
442 0 : zout(2, nout4, j) = r - s
443 : END DO
444 : END DO
445 : END IF
446 : END DO
447 : ELSE
448 49350 : ia = 1
449 49350 : nin1 = ia - after
450 49350 : nout1 = ia - atn
451 98700 : DO ib = 1, before
452 49350 : nin1 = nin1 + after
453 49350 : nin2 = nin1 + atb
454 49350 : nin3 = nin2 + atb
455 49350 : nin4 = nin3 + atb
456 49350 : nout1 = nout1 + atn
457 49350 : nout2 = nout1 + after
458 49350 : nout3 = nout2 + after
459 49350 : nout4 = nout3 + after
460 820620 : DO j = 1, nfft
461 721920 : r1 = zin(1, j, nin1)
462 721920 : s1 = zin(2, j, nin1)
463 721920 : r2 = zin(1, j, nin2)
464 721920 : s2 = zin(2, j, nin2)
465 721920 : r3 = zin(1, j, nin3)
466 721920 : s3 = zin(2, j, nin3)
467 721920 : r4 = zin(1, j, nin4)
468 721920 : s4 = zin(2, j, nin4)
469 721920 : r = r1 + r3
470 721920 : s = r2 + r4
471 721920 : zout(1, nout1, j) = r + s
472 721920 : zout(1, nout3, j) = r - s
473 721920 : r = r1 - r3
474 721920 : s = s2 - s4
475 721920 : zout(1, nout2, j) = r + s
476 721920 : zout(1, nout4, j) = r - s
477 721920 : r = s1 + s3
478 721920 : s = s2 + s4
479 721920 : zout(2, nout1, j) = r + s
480 721920 : zout(2, nout3, j) = r - s
481 721920 : r = s1 - s3
482 721920 : s = r2 - r4
483 721920 : zout(2, nout2, j) = r - s
484 771270 : zout(2, nout4, j) = r + s
485 : END DO
486 : END DO
487 394800 : DO ia = 2, after
488 345450 : ias = ia - 1
489 394800 : IF (2*ias == after) THEN
490 49350 : nin1 = ia - after
491 49350 : nout1 = ia - atn
492 98700 : DO ib = 1, before
493 49350 : nin1 = nin1 + after
494 49350 : nin2 = nin1 + atb
495 49350 : nin3 = nin2 + atb
496 49350 : nin4 = nin3 + atb
497 49350 : nout1 = nout1 + atn
498 49350 : nout2 = nout1 + after
499 49350 : nout3 = nout2 + after
500 49350 : nout4 = nout3 + after
501 820620 : DO j = 1, nfft
502 721920 : r1 = zin(1, j, nin1)
503 721920 : s1 = zin(2, j, nin1)
504 721920 : r = zin(1, j, nin2)
505 721920 : s = zin(2, j, nin2)
506 721920 : r2 = (r + s)*rt2i
507 721920 : s2 = (s - r)*rt2i
508 721920 : r3 = zin(2, j, nin3)
509 721920 : s3 = -zin(1, j, nin3)
510 721920 : r = zin(1, j, nin4)
511 721920 : s = zin(2, j, nin4)
512 721920 : r4 = (s - r)*rt2i
513 721920 : s4 = -(r + s)*rt2i
514 721920 : r = r1 + r3
515 721920 : s = r2 + r4
516 721920 : zout(1, nout1, j) = r + s
517 721920 : zout(1, nout3, j) = r - s
518 721920 : r = r1 - r3
519 721920 : s = s2 - s4
520 721920 : zout(1, nout2, j) = r + s
521 721920 : zout(1, nout4, j) = r - s
522 721920 : r = s1 + s3
523 721920 : s = s2 + s4
524 721920 : zout(2, nout1, j) = r + s
525 721920 : zout(2, nout3, j) = r - s
526 721920 : r = s1 - s3
527 721920 : s = r2 - r4
528 721920 : zout(2, nout2, j) = r - s
529 771270 : zout(2, nout4, j) = r + s
530 : END DO
531 : END DO
532 : ELSE
533 296100 : itt = ias*before
534 296100 : itrig = itt + 1
535 296100 : cr2 = trig(1, itrig)
536 296100 : ci2 = trig(2, itrig)
537 296100 : itrig = itrig + itt
538 296100 : cr3 = trig(1, itrig)
539 296100 : ci3 = trig(2, itrig)
540 296100 : itrig = itrig + itt
541 296100 : cr4 = trig(1, itrig)
542 296100 : ci4 = trig(2, itrig)
543 296100 : nin1 = ia - after
544 296100 : nout1 = ia - atn
545 592200 : DO ib = 1, before
546 296100 : nin1 = nin1 + after
547 296100 : nin2 = nin1 + atb
548 296100 : nin3 = nin2 + atb
549 296100 : nin4 = nin3 + atb
550 296100 : nout1 = nout1 + atn
551 296100 : nout2 = nout1 + after
552 296100 : nout3 = nout2 + after
553 296100 : nout4 = nout3 + after
554 4923720 : DO j = 1, nfft
555 4331520 : r1 = zin(1, j, nin1)
556 4331520 : s1 = zin(2, j, nin1)
557 4331520 : r = zin(1, j, nin2)
558 4331520 : s = zin(2, j, nin2)
559 4331520 : r2 = r*cr2 - s*ci2
560 4331520 : s2 = r*ci2 + s*cr2
561 4331520 : r = zin(1, j, nin3)
562 4331520 : s = zin(2, j, nin3)
563 4331520 : r3 = r*cr3 - s*ci3
564 4331520 : s3 = r*ci3 + s*cr3
565 4331520 : r = zin(1, j, nin4)
566 4331520 : s = zin(2, j, nin4)
567 4331520 : r4 = r*cr4 - s*ci4
568 4331520 : s4 = r*ci4 + s*cr4
569 4331520 : r = r1 + r3
570 4331520 : s = r2 + r4
571 4331520 : zout(1, nout1, j) = r + s
572 4331520 : zout(1, nout3, j) = r - s
573 4331520 : r = r1 - r3
574 4331520 : s = s2 - s4
575 4331520 : zout(1, nout2, j) = r + s
576 4331520 : zout(1, nout4, j) = r - s
577 4331520 : r = s1 + s3
578 4331520 : s = s2 + s4
579 4331520 : zout(2, nout1, j) = r + s
580 4331520 : zout(2, nout3, j) = r - s
581 4331520 : r = s1 - s3
582 4331520 : s = r2 - r4
583 4331520 : zout(2, nout2, j) = r - s
584 4627620 : zout(2, nout4, j) = r + s
585 : END DO
586 : END DO
587 : END IF
588 : END DO
589 : END IF
590 248795 : ELSE IF (now == 8) THEN
591 0 : IF (isign == -1) THEN
592 0 : ia = 1
593 0 : nin1 = ia - after
594 0 : nout1 = ia - atn
595 0 : DO ib = 1, before
596 0 : nin1 = nin1 + after
597 0 : nin2 = nin1 + atb
598 0 : nin3 = nin2 + atb
599 0 : nin4 = nin3 + atb
600 0 : nin5 = nin4 + atb
601 0 : nin6 = nin5 + atb
602 0 : nin7 = nin6 + atb
603 0 : nin8 = nin7 + atb
604 0 : nout1 = nout1 + atn
605 0 : nout2 = nout1 + after
606 0 : nout3 = nout2 + after
607 0 : nout4 = nout3 + after
608 0 : nout5 = nout4 + after
609 0 : nout6 = nout5 + after
610 0 : nout7 = nout6 + after
611 0 : nout8 = nout7 + after
612 0 : DO j = 1, nfft
613 0 : r1 = zin(1, j, nin1)
614 0 : s1 = zin(2, j, nin1)
615 0 : r2 = zin(1, j, nin2)
616 0 : s2 = zin(2, j, nin2)
617 0 : r3 = zin(1, j, nin3)
618 0 : s3 = zin(2, j, nin3)
619 0 : r4 = zin(1, j, nin4)
620 0 : s4 = zin(2, j, nin4)
621 0 : r5 = zin(1, j, nin5)
622 0 : s5 = zin(2, j, nin5)
623 0 : r6 = zin(1, j, nin6)
624 0 : s6 = zin(2, j, nin6)
625 0 : r7 = zin(1, j, nin7)
626 0 : s7 = zin(2, j, nin7)
627 0 : r8 = zin(1, j, nin8)
628 0 : s8 = zin(2, j, nin8)
629 0 : r = r1 + r5
630 0 : s = r3 + r7
631 0 : ap = r + s
632 0 : am = r - s
633 0 : r = r2 + r6
634 0 : s = r4 + r8
635 0 : bp = r + s
636 0 : bm = r - s
637 0 : r = s1 + s5
638 0 : s = s3 + s7
639 0 : cp = r + s
640 0 : cm = r - s
641 0 : r = s2 + s6
642 0 : s = s4 + s8
643 0 : dbl = r + s
644 0 : dm = r - s
645 0 : zout(1, nout1, j) = ap + bp
646 0 : zout(2, nout1, j) = cp + dbl
647 0 : zout(1, nout5, j) = ap - bp
648 0 : zout(2, nout5, j) = cp - dbl
649 0 : zout(1, nout3, j) = am + dm
650 0 : zout(2, nout3, j) = cm - bm
651 0 : zout(1, nout7, j) = am - dm
652 0 : zout(2, nout7, j) = cm + bm
653 0 : r = r1 - r5
654 0 : s = s3 - s7
655 0 : ap = r + s
656 0 : am = r - s
657 0 : r = s1 - s5
658 0 : s = r3 - r7
659 0 : bp = r + s
660 0 : bm = r - s
661 0 : r = s4 - s8
662 0 : s = r2 - r6
663 0 : cp = r + s
664 0 : cm = r - s
665 0 : r = s2 - s6
666 0 : s = r4 - r8
667 0 : dbl = r + s
668 0 : dm = r - s
669 0 : r = (cp + dm)*rt2i
670 0 : s = (-cp + dm)*rt2i
671 0 : cp = (cm + dbl)*rt2i
672 0 : dbl = (cm - dbl)*rt2i
673 0 : zout(1, nout2, j) = ap + r
674 0 : zout(2, nout2, j) = bm + s
675 0 : zout(1, nout6, j) = ap - r
676 0 : zout(2, nout6, j) = bm - s
677 0 : zout(1, nout4, j) = am + cp
678 0 : zout(2, nout4, j) = bp + dbl
679 0 : zout(1, nout8, j) = am - cp
680 0 : zout(2, nout8, j) = bp - dbl
681 : END DO
682 : END DO
683 : ELSE
684 0 : ia = 1
685 0 : nin1 = ia - after
686 0 : nout1 = ia - atn
687 0 : DO ib = 1, before
688 0 : nin1 = nin1 + after
689 0 : nin2 = nin1 + atb
690 0 : nin3 = nin2 + atb
691 0 : nin4 = nin3 + atb
692 0 : nin5 = nin4 + atb
693 0 : nin6 = nin5 + atb
694 0 : nin7 = nin6 + atb
695 0 : nin8 = nin7 + atb
696 0 : nout1 = nout1 + atn
697 0 : nout2 = nout1 + after
698 0 : nout3 = nout2 + after
699 0 : nout4 = nout3 + after
700 0 : nout5 = nout4 + after
701 0 : nout6 = nout5 + after
702 0 : nout7 = nout6 + after
703 0 : nout8 = nout7 + after
704 0 : DO j = 1, nfft
705 0 : r1 = zin(1, j, nin1)
706 0 : s1 = zin(2, j, nin1)
707 0 : r2 = zin(1, j, nin2)
708 0 : s2 = zin(2, j, nin2)
709 0 : r3 = zin(1, j, nin3)
710 0 : s3 = zin(2, j, nin3)
711 0 : r4 = zin(1, j, nin4)
712 0 : s4 = zin(2, j, nin4)
713 0 : r5 = zin(1, j, nin5)
714 0 : s5 = zin(2, j, nin5)
715 0 : r6 = zin(1, j, nin6)
716 0 : s6 = zin(2, j, nin6)
717 0 : r7 = zin(1, j, nin7)
718 0 : s7 = zin(2, j, nin7)
719 0 : r8 = zin(1, j, nin8)
720 0 : s8 = zin(2, j, nin8)
721 0 : r = r1 + r5
722 0 : s = r3 + r7
723 0 : ap = r + s
724 0 : am = r - s
725 0 : r = r2 + r6
726 0 : s = r4 + r8
727 0 : bp = r + s
728 0 : bm = r - s
729 0 : r = s1 + s5
730 0 : s = s3 + s7
731 0 : cp = r + s
732 0 : cm = r - s
733 0 : r = s2 + s6
734 0 : s = s4 + s8
735 0 : dbl = r + s
736 0 : dm = r - s
737 0 : zout(1, nout1, j) = ap + bp
738 0 : zout(2, nout1, j) = cp + dbl
739 0 : zout(1, nout5, j) = ap - bp
740 0 : zout(2, nout5, j) = cp - dbl
741 0 : zout(1, nout3, j) = am - dm
742 0 : zout(2, nout3, j) = cm + bm
743 0 : zout(1, nout7, j) = am + dm
744 0 : zout(2, nout7, j) = cm - bm
745 0 : r = r1 - r5
746 0 : s = -s3 + s7
747 0 : ap = r + s
748 0 : am = r - s
749 0 : r = s1 - s5
750 0 : s = r7 - r3
751 0 : bp = r + s
752 0 : bm = r - s
753 0 : r = -s4 + s8
754 0 : s = r2 - r6
755 0 : cp = r + s
756 0 : cm = r - s
757 0 : r = -s2 + s6
758 0 : s = r4 - r8
759 0 : dbl = r + s
760 0 : dm = r - s
761 0 : r = (cp + dm)*rt2i
762 0 : s = (cp - dm)*rt2i
763 0 : cp = (cm + dbl)*rt2i
764 0 : dbl = (-cm + dbl)*rt2i
765 0 : zout(1, nout2, j) = ap + r
766 0 : zout(2, nout2, j) = bm + s
767 0 : zout(1, nout6, j) = ap - r
768 0 : zout(2, nout6, j) = bm - s
769 0 : zout(1, nout4, j) = am + cp
770 0 : zout(2, nout4, j) = bp + dbl
771 0 : zout(1, nout8, j) = am - cp
772 0 : zout(2, nout8, j) = bp - dbl
773 : END DO
774 : END DO
775 : END IF
776 248795 : ELSE IF (now == 3) THEN
777 216635 : bbs = isign*bb
778 216635 : ia = 1
779 216635 : nin1 = ia - after
780 216635 : nout1 = ia - atn
781 433270 : DO ib = 1, before
782 216635 : nin1 = nin1 + after
783 216635 : nin2 = nin1 + atb
784 216635 : nin3 = nin2 + atb
785 216635 : nout1 = nout1 + atn
786 216635 : nout2 = nout1 + after
787 216635 : nout3 = nout2 + after
788 1843497 : DO j = 1, nfft
789 1410227 : r1 = zin(1, j, nin1)
790 1410227 : s1 = zin(2, j, nin1)
791 1410227 : r2 = zin(1, j, nin2)
792 1410227 : s2 = zin(2, j, nin2)
793 1410227 : r3 = zin(1, j, nin3)
794 1410227 : s3 = zin(2, j, nin3)
795 1410227 : r = r2 + r3
796 1410227 : s = s2 + s3
797 1410227 : zout(1, nout1, j) = r + r1
798 1410227 : zout(2, nout1, j) = s + s1
799 1410227 : r1 = r1 - 0.5_dp*r
800 1410227 : s1 = s1 - 0.5_dp*s
801 1410227 : r2 = bbs*(r2 - r3)
802 1410227 : s2 = bbs*(s2 - s3)
803 1410227 : zout(1, nout2, j) = r1 - s2
804 1410227 : zout(2, nout2, j) = s1 + r2
805 1410227 : zout(1, nout3, j) = r1 + s2
806 1626862 : zout(2, nout3, j) = s1 - r2
807 : END DO
808 : END DO
809 5901760 : DO ia = 2, after
810 5685125 : ias = ia - 1
811 5901760 : IF (4*ias == 3*after) THEN
812 16775 : IF (isign == 1) THEN
813 0 : nin1 = ia - after
814 0 : nout1 = ia - atn
815 0 : DO ib = 1, before
816 0 : nin1 = nin1 + after
817 0 : nin2 = nin1 + atb
818 0 : nin3 = nin2 + atb
819 0 : nout1 = nout1 + atn
820 0 : nout2 = nout1 + after
821 0 : nout3 = nout2 + after
822 0 : DO j = 1, nfft
823 0 : r1 = zin(1, j, nin1)
824 0 : s1 = zin(2, j, nin1)
825 0 : r2 = -zin(2, j, nin2)
826 0 : s2 = zin(1, j, nin2)
827 0 : r3 = -zin(1, j, nin3)
828 0 : s3 = -zin(2, j, nin3)
829 0 : r = r2 + r3
830 0 : s = s2 + s3
831 0 : zout(1, nout1, j) = r + r1
832 0 : zout(2, nout1, j) = s + s1
833 0 : r1 = r1 - 0.5_dp*r
834 0 : s1 = s1 - 0.5_dp*s
835 0 : r2 = bbs*(r2 - r3)
836 0 : s2 = bbs*(s2 - s3)
837 0 : zout(1, nout2, j) = r1 - s2
838 0 : zout(2, nout2, j) = s1 + r2
839 0 : zout(1, nout3, j) = r1 + s2
840 0 : zout(2, nout3, j) = s1 - r2
841 : END DO
842 : END DO
843 : ELSE
844 16775 : nin1 = ia - after
845 16775 : nout1 = ia - atn
846 33550 : DO ib = 1, before
847 16775 : nin1 = nin1 + after
848 16775 : nin2 = nin1 + atb
849 16775 : nin3 = nin2 + atb
850 16775 : nout1 = nout1 + atn
851 16775 : nout2 = nout1 + after
852 16775 : nout3 = nout2 + after
853 188925 : DO j = 1, nfft
854 155375 : r1 = zin(1, j, nin1)
855 155375 : s1 = zin(2, j, nin1)
856 155375 : r2 = zin(2, j, nin2)
857 155375 : s2 = -zin(1, j, nin2)
858 155375 : r3 = -zin(1, j, nin3)
859 155375 : s3 = -zin(2, j, nin3)
860 155375 : r = r2 + r3
861 155375 : s = s2 + s3
862 155375 : zout(1, nout1, j) = r + r1
863 155375 : zout(2, nout1, j) = s + s1
864 155375 : r1 = r1 - 0.5_dp*r
865 155375 : s1 = s1 - 0.5_dp*s
866 155375 : r2 = bbs*(r2 - r3)
867 155375 : s2 = bbs*(s2 - s3)
868 155375 : zout(1, nout2, j) = r1 - s2
869 155375 : zout(2, nout2, j) = s1 + r2
870 155375 : zout(1, nout3, j) = r1 + s2
871 172150 : zout(2, nout3, j) = s1 - r2
872 : END DO
873 : END DO
874 : END IF
875 5668350 : ELSE IF (8*ias == 3*after) THEN
876 0 : IF (isign == 1) THEN
877 0 : nin1 = ia - after
878 0 : nout1 = ia - atn
879 0 : DO ib = 1, before
880 0 : nin1 = nin1 + after
881 0 : nin2 = nin1 + atb
882 0 : nin3 = nin2 + atb
883 0 : nout1 = nout1 + atn
884 0 : nout2 = nout1 + after
885 0 : nout3 = nout2 + after
886 0 : DO j = 1, nfft
887 0 : r1 = zin(1, j, nin1)
888 0 : s1 = zin(2, j, nin1)
889 0 : r = zin(1, j, nin2)
890 0 : s = zin(2, j, nin2)
891 0 : r2 = (r - s)*rt2i
892 0 : s2 = (r + s)*rt2i
893 0 : r3 = -zin(2, j, nin3)
894 0 : s3 = zin(1, j, nin3)
895 0 : r = r2 + r3
896 0 : s = s2 + s3
897 0 : zout(1, nout1, j) = r + r1
898 0 : zout(2, nout1, j) = s + s1
899 0 : r1 = r1 - 0.5_dp*r
900 0 : s1 = s1 - 0.5_dp*s
901 0 : r2 = bbs*(r2 - r3)
902 0 : s2 = bbs*(s2 - s3)
903 0 : zout(1, nout2, j) = r1 - s2
904 0 : zout(2, nout2, j) = s1 + r2
905 0 : zout(1, nout3, j) = r1 + s2
906 0 : zout(2, nout3, j) = s1 - r2
907 : END DO
908 : END DO
909 : ELSE
910 0 : nin1 = ia - after
911 0 : nout1 = ia - atn
912 0 : DO ib = 1, before
913 0 : nin1 = nin1 + after
914 0 : nin2 = nin1 + atb
915 0 : nin3 = nin2 + atb
916 0 : nout1 = nout1 + atn
917 0 : nout2 = nout1 + after
918 0 : nout3 = nout2 + after
919 0 : DO j = 1, nfft
920 0 : r1 = zin(1, j, nin1)
921 0 : s1 = zin(2, j, nin1)
922 0 : r = zin(1, j, nin2)
923 0 : s = zin(2, j, nin2)
924 0 : r2 = (r + s)*rt2i
925 0 : s2 = (-r + s)*rt2i
926 0 : r3 = zin(2, j, nin3)
927 0 : s3 = -zin(1, j, nin3)
928 0 : r = r2 + r3
929 0 : s = s2 + s3
930 0 : zout(1, nout1, j) = r + r1
931 0 : zout(2, nout1, j) = s + s1
932 0 : r1 = r1 - 0.5_dp*r
933 0 : s1 = s1 - 0.5_dp*s
934 0 : r2 = bbs*(r2 - r3)
935 0 : s2 = bbs*(s2 - s3)
936 0 : zout(1, nout2, j) = r1 - s2
937 0 : zout(2, nout2, j) = s1 + r2
938 0 : zout(1, nout3, j) = r1 + s2
939 0 : zout(2, nout3, j) = s1 - r2
940 : END DO
941 : END DO
942 : END IF
943 : ELSE
944 5668350 : itt = ias*before
945 5668350 : itrig = itt + 1
946 5668350 : cr2 = trig(1, itrig)
947 5668350 : ci2 = trig(2, itrig)
948 5668350 : itrig = itrig + itt
949 5668350 : cr3 = trig(1, itrig)
950 5668350 : ci3 = trig(2, itrig)
951 5668350 : nin1 = ia - after
952 5668350 : nout1 = ia - atn
953 11336700 : DO ib = 1, before
954 5668350 : nin1 = nin1 + after
955 5668350 : nin2 = nin1 + atb
956 5668350 : nin3 = nin2 + atb
957 5668350 : nout1 = nout1 + atn
958 5668350 : nout2 = nout1 + after
959 5668350 : nout3 = nout2 + after
960 42300822 : DO j = 1, nfft
961 30964122 : r1 = zin(1, j, nin1)
962 30964122 : s1 = zin(2, j, nin1)
963 30964122 : r = zin(1, j, nin2)
964 30964122 : s = zin(2, j, nin2)
965 30964122 : r2 = r*cr2 - s*ci2
966 30964122 : s2 = r*ci2 + s*cr2
967 30964122 : r = zin(1, j, nin3)
968 30964122 : s = zin(2, j, nin3)
969 30964122 : r3 = r*cr3 - s*ci3
970 30964122 : s3 = r*ci3 + s*cr3
971 30964122 : r = r2 + r3
972 30964122 : s = s2 + s3
973 30964122 : zout(1, nout1, j) = r + r1
974 30964122 : zout(2, nout1, j) = s + s1
975 30964122 : r1 = r1 - 0.5_dp*r
976 30964122 : s1 = s1 - 0.5_dp*s
977 30964122 : r2 = bbs*(r2 - r3)
978 30964122 : s2 = bbs*(s2 - s3)
979 30964122 : zout(1, nout2, j) = r1 - s2
980 30964122 : zout(2, nout2, j) = s1 + r2
981 30964122 : zout(1, nout3, j) = r1 + s2
982 36632472 : zout(2, nout3, j) = s1 - r2
983 : END DO
984 : END DO
985 : END IF
986 : END DO
987 32160 : ELSE IF (now == 5) THEN
988 32160 : sin2 = isign*sin2p
989 32160 : sin4 = isign*sin4p
990 32160 : ia = 1
991 32160 : nin1 = ia - after
992 32160 : nout1 = ia - atn
993 64320 : DO ib = 1, before
994 32160 : nin1 = nin1 + after
995 32160 : nin2 = nin1 + atb
996 32160 : nin3 = nin2 + atb
997 32160 : nin4 = nin3 + atb
998 32160 : nin5 = nin4 + atb
999 32160 : nout1 = nout1 + atn
1000 32160 : nout2 = nout1 + after
1001 32160 : nout3 = nout2 + after
1002 32160 : nout4 = nout3 + after
1003 32160 : nout5 = nout4 + after
1004 652245 : DO j = 1, nfft
1005 587925 : r1 = zin(1, j, nin1)
1006 587925 : s1 = zin(2, j, nin1)
1007 587925 : r2 = zin(1, j, nin2)
1008 587925 : s2 = zin(2, j, nin2)
1009 587925 : r3 = zin(1, j, nin3)
1010 587925 : s3 = zin(2, j, nin3)
1011 587925 : r4 = zin(1, j, nin4)
1012 587925 : s4 = zin(2, j, nin4)
1013 587925 : r5 = zin(1, j, nin5)
1014 587925 : s5 = zin(2, j, nin5)
1015 587925 : r25 = r2 + r5
1016 587925 : r34 = r3 + r4
1017 587925 : s25 = s2 - s5
1018 587925 : s34 = s3 - s4
1019 587925 : zout(1, nout1, j) = r1 + r25 + r34
1020 587925 : r = cos2*r25 + cos4*r34 + r1
1021 587925 : s = sin2*s25 + sin4*s34
1022 587925 : zout(1, nout2, j) = r - s
1023 587925 : zout(1, nout5, j) = r + s
1024 587925 : r = cos4*r25 + cos2*r34 + r1
1025 587925 : s = sin4*s25 - sin2*s34
1026 587925 : zout(1, nout3, j) = r - s
1027 587925 : zout(1, nout4, j) = r + s
1028 587925 : r25 = r2 - r5
1029 587925 : r34 = r3 - r4
1030 587925 : s25 = s2 + s5
1031 587925 : s34 = s3 + s4
1032 587925 : zout(2, nout1, j) = s1 + s25 + s34
1033 587925 : r = cos2*s25 + cos4*s34 + s1
1034 587925 : s = sin2*r25 + sin4*r34
1035 587925 : zout(2, nout2, j) = r + s
1036 587925 : zout(2, nout5, j) = r - s
1037 587925 : r = cos4*s25 + cos2*s34 + s1
1038 587925 : s = sin4*r25 - sin2*r34
1039 587925 : zout(2, nout3, j) = r + s
1040 620085 : zout(2, nout4, j) = r - s
1041 : END DO
1042 : END DO
1043 162114 : DO ia = 2, after
1044 129954 : ias = ia - 1
1045 162114 : IF (8*ias == 5*after) THEN
1046 438 : IF (isign == 1) THEN
1047 0 : nin1 = ia - after
1048 0 : nout1 = ia - atn
1049 0 : DO ib = 1, before
1050 0 : nin1 = nin1 + after
1051 0 : nin2 = nin1 + atb
1052 0 : nin3 = nin2 + atb
1053 0 : nin4 = nin3 + atb
1054 0 : nin5 = nin4 + atb
1055 0 : nout1 = nout1 + atn
1056 0 : nout2 = nout1 + after
1057 0 : nout3 = nout2 + after
1058 0 : nout4 = nout3 + after
1059 0 : nout5 = nout4 + after
1060 0 : DO j = 1, nfft
1061 0 : r1 = zin(1, j, nin1)
1062 0 : s1 = zin(2, j, nin1)
1063 0 : r = zin(1, j, nin2)
1064 0 : s = zin(2, j, nin2)
1065 0 : r2 = (r - s)*rt2i
1066 0 : s2 = (r + s)*rt2i
1067 0 : r3 = -zin(2, j, nin3)
1068 0 : s3 = zin(1, j, nin3)
1069 0 : r = zin(1, j, nin4)
1070 0 : s = zin(2, j, nin4)
1071 0 : r4 = -(r + s)*rt2i
1072 0 : s4 = (r - s)*rt2i
1073 0 : r5 = -zin(1, j, nin5)
1074 0 : s5 = -zin(2, j, nin5)
1075 0 : r25 = r2 + r5
1076 0 : r34 = r3 + r4
1077 0 : s25 = s2 - s5
1078 0 : s34 = s3 - s4
1079 0 : zout(1, nout1, j) = r1 + r25 + r34
1080 0 : r = cos2*r25 + cos4*r34 + r1
1081 0 : s = sin2*s25 + sin4*s34
1082 0 : zout(1, nout2, j) = r - s
1083 0 : zout(1, nout5, j) = r + s
1084 0 : r = cos4*r25 + cos2*r34 + r1
1085 0 : s = sin4*s25 - sin2*s34
1086 0 : zout(1, nout3, j) = r - s
1087 0 : zout(1, nout4, j) = r + s
1088 0 : r25 = r2 - r5
1089 0 : r34 = r3 - r4
1090 0 : s25 = s2 + s5
1091 0 : s34 = s3 + s4
1092 0 : zout(2, nout1, j) = s1 + s25 + s34
1093 0 : r = cos2*s25 + cos4*s34 + s1
1094 0 : s = sin2*r25 + sin4*r34
1095 0 : zout(2, nout2, j) = r + s
1096 0 : zout(2, nout5, j) = r - s
1097 0 : r = cos4*s25 + cos2*s34 + s1
1098 0 : s = sin4*r25 - sin2*r34
1099 0 : zout(2, nout3, j) = r + s
1100 0 : zout(2, nout4, j) = r - s
1101 : END DO
1102 : END DO
1103 : ELSE
1104 438 : nin1 = ia - after
1105 438 : nout1 = ia - atn
1106 876 : DO ib = 1, before
1107 438 : nin1 = nin1 + after
1108 438 : nin2 = nin1 + atb
1109 438 : nin3 = nin2 + atb
1110 438 : nin4 = nin3 + atb
1111 438 : nin5 = nin4 + atb
1112 438 : nout1 = nout1 + atn
1113 438 : nout2 = nout1 + after
1114 438 : nout3 = nout2 + after
1115 438 : nout4 = nout3 + after
1116 438 : nout5 = nout4 + after
1117 5676 : DO j = 1, nfft
1118 4800 : r1 = zin(1, j, nin1)
1119 4800 : s1 = zin(2, j, nin1)
1120 4800 : r = zin(1, j, nin2)
1121 4800 : s = zin(2, j, nin2)
1122 4800 : r2 = (r + s)*rt2i
1123 4800 : s2 = (-r + s)*rt2i
1124 4800 : r3 = zin(2, j, nin3)
1125 4800 : s3 = -zin(1, j, nin3)
1126 4800 : r = zin(1, j, nin4)
1127 4800 : s = zin(2, j, nin4)
1128 4800 : r4 = (s - r)*rt2i
1129 4800 : s4 = -(r + s)*rt2i
1130 4800 : r5 = -zin(1, j, nin5)
1131 4800 : s5 = -zin(2, j, nin5)
1132 4800 : r25 = r2 + r5
1133 4800 : r34 = r3 + r4
1134 4800 : s25 = s2 - s5
1135 4800 : s34 = s3 - s4
1136 4800 : zout(1, nout1, j) = r1 + r25 + r34
1137 4800 : r = cos2*r25 + cos4*r34 + r1
1138 4800 : s = sin2*s25 + sin4*s34
1139 4800 : zout(1, nout2, j) = r - s
1140 4800 : zout(1, nout5, j) = r + s
1141 4800 : r = cos4*r25 + cos2*r34 + r1
1142 4800 : s = sin4*s25 - sin2*s34
1143 4800 : zout(1, nout3, j) = r - s
1144 4800 : zout(1, nout4, j) = r + s
1145 4800 : r25 = r2 - r5
1146 4800 : r34 = r3 - r4
1147 4800 : s25 = s2 + s5
1148 4800 : s34 = s3 + s4
1149 4800 : zout(2, nout1, j) = s1 + s25 + s34
1150 4800 : r = cos2*s25 + cos4*s34 + s1
1151 4800 : s = sin2*r25 + sin4*r34
1152 4800 : zout(2, nout2, j) = r + s
1153 4800 : zout(2, nout5, j) = r - s
1154 4800 : r = cos4*s25 + cos2*s34 + s1
1155 4800 : s = sin4*r25 - sin2*r34
1156 4800 : zout(2, nout3, j) = r + s
1157 5238 : zout(2, nout4, j) = r - s
1158 : END DO
1159 : END DO
1160 : END IF
1161 : ELSE
1162 129516 : ias = ia - 1
1163 129516 : itt = ias*before
1164 129516 : itrig = itt + 1
1165 129516 : cr2 = trig(1, itrig)
1166 129516 : ci2 = trig(2, itrig)
1167 129516 : itrig = itrig + itt
1168 129516 : cr3 = trig(1, itrig)
1169 129516 : ci3 = trig(2, itrig)
1170 129516 : itrig = itrig + itt
1171 129516 : cr4 = trig(1, itrig)
1172 129516 : ci4 = trig(2, itrig)
1173 129516 : itrig = itrig + itt
1174 129516 : cr5 = trig(1, itrig)
1175 129516 : ci5 = trig(2, itrig)
1176 129516 : nin1 = ia - after
1177 129516 : nout1 = ia - atn
1178 259032 : DO ib = 1, before
1179 129516 : nin1 = nin1 + after
1180 129516 : nin2 = nin1 + atb
1181 129516 : nin3 = nin2 + atb
1182 129516 : nin4 = nin3 + atb
1183 129516 : nin5 = nin4 + atb
1184 129516 : nout1 = nout1 + atn
1185 129516 : nout2 = nout1 + after
1186 129516 : nout3 = nout2 + after
1187 129516 : nout4 = nout3 + after
1188 129516 : nout5 = nout4 + after
1189 2620332 : DO j = 1, nfft
1190 2361300 : r1 = zin(1, j, nin1)
1191 2361300 : s1 = zin(2, j, nin1)
1192 2361300 : r = zin(1, j, nin2)
1193 2361300 : s = zin(2, j, nin2)
1194 2361300 : r2 = r*cr2 - s*ci2
1195 2361300 : s2 = r*ci2 + s*cr2
1196 2361300 : r = zin(1, j, nin3)
1197 2361300 : s = zin(2, j, nin3)
1198 2361300 : r3 = r*cr3 - s*ci3
1199 2361300 : s3 = r*ci3 + s*cr3
1200 2361300 : r = zin(1, j, nin4)
1201 2361300 : s = zin(2, j, nin4)
1202 2361300 : r4 = r*cr4 - s*ci4
1203 2361300 : s4 = r*ci4 + s*cr4
1204 2361300 : r = zin(1, j, nin5)
1205 2361300 : s = zin(2, j, nin5)
1206 2361300 : r5 = r*cr5 - s*ci5
1207 2361300 : s5 = r*ci5 + s*cr5
1208 2361300 : r25 = r2 + r5
1209 2361300 : r34 = r3 + r4
1210 2361300 : s25 = s2 - s5
1211 2361300 : s34 = s3 - s4
1212 2361300 : zout(1, nout1, j) = r1 + r25 + r34
1213 2361300 : r = cos2*r25 + cos4*r34 + r1
1214 2361300 : s = sin2*s25 + sin4*s34
1215 2361300 : zout(1, nout2, j) = r - s
1216 2361300 : zout(1, nout5, j) = r + s
1217 2361300 : r = cos4*r25 + cos2*r34 + r1
1218 2361300 : s = sin4*s25 - sin2*s34
1219 2361300 : zout(1, nout3, j) = r - s
1220 2361300 : zout(1, nout4, j) = r + s
1221 2361300 : r25 = r2 - r5
1222 2361300 : r34 = r3 - r4
1223 2361300 : s25 = s2 + s5
1224 2361300 : s34 = s3 + s4
1225 2361300 : zout(2, nout1, j) = s1 + s25 + s34
1226 2361300 : r = cos2*s25 + cos4*s34 + s1
1227 2361300 : s = sin2*r25 + sin4*r34
1228 2361300 : zout(2, nout2, j) = r + s
1229 2361300 : zout(2, nout5, j) = r - s
1230 2361300 : r = cos4*s25 + cos2*s34 + s1
1231 2361300 : s = sin4*r25 - sin2*r34
1232 2361300 : zout(2, nout3, j) = r + s
1233 2490816 : zout(2, nout4, j) = r - s
1234 : END DO
1235 : END DO
1236 : END IF
1237 : END DO
1238 0 : ELSE IF (now == 6) THEN
1239 0 : bbs = isign*bb
1240 0 : ia = 1
1241 0 : nin1 = ia - after
1242 0 : nout1 = ia - atn
1243 0 : DO ib = 1, before
1244 0 : nin1 = nin1 + after
1245 0 : nin2 = nin1 + atb
1246 0 : nin3 = nin2 + atb
1247 0 : nin4 = nin3 + atb
1248 0 : nin5 = nin4 + atb
1249 0 : nin6 = nin5 + atb
1250 0 : nout1 = nout1 + atn
1251 0 : nout2 = nout1 + after
1252 0 : nout3 = nout2 + after
1253 0 : nout4 = nout3 + after
1254 0 : nout5 = nout4 + after
1255 0 : nout6 = nout5 + after
1256 0 : DO j = 1, nfft
1257 0 : r2 = zin(1, j, nin3)
1258 0 : s2 = zin(2, j, nin3)
1259 0 : r3 = zin(1, j, nin5)
1260 0 : s3 = zin(2, j, nin5)
1261 0 : r = r2 + r3
1262 0 : s = s2 + s3
1263 0 : r1 = zin(1, j, nin1)
1264 0 : s1 = zin(2, j, nin1)
1265 0 : ur1 = r + r1
1266 0 : ui1 = s + s1
1267 0 : r1 = r1 - 0.5_dp*r
1268 0 : s1 = s1 - 0.5_dp*s
1269 0 : r = r2 - r3
1270 0 : s = s2 - s3
1271 0 : ur2 = r1 - s*bbs
1272 0 : ui2 = s1 + r*bbs
1273 0 : ur3 = r1 + s*bbs
1274 0 : ui3 = s1 - r*bbs
1275 :
1276 0 : r2 = zin(1, j, nin6)
1277 0 : s2 = zin(2, j, nin6)
1278 0 : r3 = zin(1, j, nin2)
1279 0 : s3 = zin(2, j, nin2)
1280 0 : r = r2 + r3
1281 0 : s = s2 + s3
1282 0 : r1 = zin(1, j, nin4)
1283 0 : s1 = zin(2, j, nin4)
1284 0 : vr1 = r + r1
1285 0 : vi1 = s + s1
1286 0 : r1 = r1 - 0.5_dp*r
1287 0 : s1 = s1 - 0.5_dp*s
1288 0 : r = r2 - r3
1289 0 : s = s2 - s3
1290 0 : vr2 = r1 - s*bbs
1291 0 : vi2 = s1 + r*bbs
1292 0 : vr3 = r1 + s*bbs
1293 0 : vi3 = s1 - r*bbs
1294 :
1295 0 : zout(1, nout1, j) = ur1 + vr1
1296 0 : zout(2, nout1, j) = ui1 + vi1
1297 0 : zout(1, nout5, j) = ur2 + vr2
1298 0 : zout(2, nout5, j) = ui2 + vi2
1299 0 : zout(1, nout3, j) = ur3 + vr3
1300 0 : zout(2, nout3, j) = ui3 + vi3
1301 0 : zout(1, nout4, j) = ur1 - vr1
1302 0 : zout(2, nout4, j) = ui1 - vi1
1303 0 : zout(1, nout2, j) = ur2 - vr2
1304 0 : zout(2, nout2, j) = ui2 - vi2
1305 0 : zout(1, nout6, j) = ur3 - vr3
1306 0 : zout(2, nout6, j) = ui3 - vi3
1307 : END DO
1308 : END DO
1309 : ELSE
1310 0 : CPABORT('Error fftrot')
1311 : END IF
1312 :
1313 : !-----------------------------------------------------------------------------!
1314 :
1315 298145 : END SUBROUTINE fftrot
1316 :
1317 : !-----------------------------------------------------------------------------!
1318 : !-----------------------------------------------------------------------------!
1319 : ! Copyright (C) Stefan Goedecker, Lausanne, Switzerland, August 1, 1991
1320 : ! Copyright (C) Stefan Goedecker, Cornell University, Ithaca, USA, 1994
1321 : ! Copyright (C) Stefan Goedecker, MPI Stuttgart, Germany, 1999
1322 : ! This file is distributed under the terms of the
1323 : ! GNU General Public License version 2 (or later),
1324 : ! see http://www.gnu.org/copyleft/gpl.txt .
1325 : !-----------------------------------------------------------------------------!
1326 : ! S. Goedecker: Rotating a three-dimensional array in optimal
1327 : ! positions for vector processing: Case study for a three-dimensional Fast
1328 : ! Fourier Transform, Comp. Phys. Commun. 76, 294 (1993)
1329 : ! **************************************************************************************************
1330 : !> \brief ...
1331 : !> \param mm ...
1332 : !> \param nfft ...
1333 : !> \param m ...
1334 : !> \param nn ...
1335 : !> \param n ...
1336 : !> \param zin ...
1337 : !> \param zout ...
1338 : !> \param trig ...
1339 : !> \param now ...
1340 : !> \param after ...
1341 : !> \param before ...
1342 : !> \param isign ...
1343 : ! **************************************************************************************************
1344 424773 : SUBROUTINE fftpre(mm, nfft, m, nn, n, zin, zout, trig, now, after, before, isign)
1345 :
1346 : INTEGER, INTENT(IN) :: mm, nfft, m, nn, n
1347 : REAL(dp), DIMENSION(2, m, mm), INTENT(IN) :: zin
1348 : REAL(dp), DIMENSION(2, nn, n), INTENT(INOUT) :: zout
1349 : REAL(dp), DIMENSION(2, ctrig_length), INTENT(IN) :: trig
1350 : INTEGER, INTENT(IN) :: now, after, before, isign
1351 :
1352 : REAL(dp), PARAMETER :: bb = 0.8660254037844387_dp, cos2 = 0.3090169943749474_dp, &
1353 : cos4 = -0.8090169943749474_dp, rt2i = 0.7071067811865475_dp, &
1354 : sin2p = 0.9510565162951536_dp, sin4p = 0.5877852522924731_dp
1355 :
1356 : INTEGER :: atb, atn, ia, ias, ib, itrig, itt, j, &
1357 : nin1, nin2, nin3, nin4, nin5, nin6, &
1358 : nin7, nin8, nout1, nout2, nout3, &
1359 : nout4, nout5, nout6, nout7, nout8
1360 : REAL(dp) :: am, ap, bbs, bm, bp, ci2, ci3, ci4, ci5, cm, cp, cr2, cr3, cr4, cr5, dbl, dm, r, &
1361 : r1, r2, r25, r3, r34, r4, r5, r6, r7, r8, s, s1, s2, s25, s3, s34, s4, s5, s6, s7, s8, &
1362 : sin2, sin4, ui1, ui2, ui3, ur1, ur2, ur3, vi1, vi2, vi3, vr1, vr2, vr3
1363 :
1364 : ! sqrt(0.5)
1365 : ! sqrt(3)/2
1366 : ! cos(2*pi/5)
1367 : ! cos(4*pi/5)
1368 : ! sin(2*pi/5)
1369 : ! sin(4*pi/5)
1370 : !-----------------------------------------------------------------------------!
1371 :
1372 424773 : atn = after*now
1373 424773 : atb = after*before
1374 :
1375 424773 : IF (now == 4) THEN
1376 60102 : IF (isign == 1) THEN
1377 30696 : ia = 1
1378 30696 : nin1 = ia - after
1379 30696 : nout1 = ia - atn
1380 343968 : DO ib = 1, before
1381 313272 : nin1 = nin1 + after
1382 313272 : nin2 = nin1 + atb
1383 313272 : nin3 = nin2 + atb
1384 313272 : nin4 = nin3 + atb
1385 313272 : nout1 = nout1 + atn
1386 313272 : nout2 = nout1 + after
1387 313272 : nout3 = nout2 + after
1388 313272 : nout4 = nout3 + after
1389 3659136 : DO j = 1, nfft
1390 3315168 : r1 = zin(1, nin1, j)
1391 3315168 : s1 = zin(2, nin1, j)
1392 3315168 : r2 = zin(1, nin2, j)
1393 3315168 : s2 = zin(2, nin2, j)
1394 3315168 : r3 = zin(1, nin3, j)
1395 3315168 : s3 = zin(2, nin3, j)
1396 3315168 : r4 = zin(1, nin4, j)
1397 3315168 : s4 = zin(2, nin4, j)
1398 3315168 : r = r1 + r3
1399 3315168 : s = r2 + r4
1400 3315168 : zout(1, j, nout1) = r + s
1401 3315168 : zout(1, j, nout3) = r - s
1402 3315168 : r = r1 - r3
1403 3315168 : s = s2 - s4
1404 3315168 : zout(1, j, nout2) = r - s
1405 3315168 : zout(1, j, nout4) = r + s
1406 3315168 : r = s1 + s3
1407 3315168 : s = s2 + s4
1408 3315168 : zout(2, j, nout1) = r + s
1409 3315168 : zout(2, j, nout3) = r - s
1410 3315168 : r = s1 - s3
1411 3315168 : s = r2 - r4
1412 3315168 : zout(2, j, nout2) = r + s
1413 3628440 : zout(2, j, nout4) = r - s
1414 : END DO
1415 : END DO
1416 30696 : DO ia = 2, after
1417 0 : ias = ia - 1
1418 30696 : IF (2*ias == after) THEN
1419 0 : nin1 = ia - after
1420 0 : nout1 = ia - atn
1421 0 : DO ib = 1, before
1422 0 : nin1 = nin1 + after
1423 0 : nin2 = nin1 + atb
1424 0 : nin3 = nin2 + atb
1425 0 : nin4 = nin3 + atb
1426 0 : nout1 = nout1 + atn
1427 0 : nout2 = nout1 + after
1428 0 : nout3 = nout2 + after
1429 0 : nout4 = nout3 + after
1430 0 : DO j = 1, nfft
1431 0 : r1 = zin(1, nin1, j)
1432 0 : s1 = zin(2, nin1, j)
1433 0 : r = zin(1, nin2, j)
1434 0 : s = zin(2, nin2, j)
1435 0 : r2 = (r - s)*rt2i
1436 0 : s2 = (r + s)*rt2i
1437 0 : r3 = -zin(2, nin3, j)
1438 0 : s3 = zin(1, nin3, j)
1439 0 : r = zin(1, nin4, j)
1440 0 : s = zin(2, nin4, j)
1441 0 : r4 = -(r + s)*rt2i
1442 0 : s4 = (r - s)*rt2i
1443 0 : r = r1 + r3
1444 0 : s = r2 + r4
1445 0 : zout(1, j, nout1) = r + s
1446 0 : zout(1, j, nout3) = r - s
1447 0 : r = r1 - r3
1448 0 : s = s2 - s4
1449 0 : zout(1, j, nout2) = r - s
1450 0 : zout(1, j, nout4) = r + s
1451 0 : r = s1 + s3
1452 0 : s = s2 + s4
1453 0 : zout(2, j, nout1) = r + s
1454 0 : zout(2, j, nout3) = r - s
1455 0 : r = s1 - s3
1456 0 : s = r2 - r4
1457 0 : zout(2, j, nout2) = r + s
1458 0 : zout(2, j, nout4) = r - s
1459 : END DO
1460 : END DO
1461 : ELSE
1462 0 : itt = ias*before
1463 0 : itrig = itt + 1
1464 0 : cr2 = trig(1, itrig)
1465 0 : ci2 = trig(2, itrig)
1466 0 : itrig = itrig + itt
1467 0 : cr3 = trig(1, itrig)
1468 0 : ci3 = trig(2, itrig)
1469 0 : itrig = itrig + itt
1470 0 : cr4 = trig(1, itrig)
1471 0 : ci4 = trig(2, itrig)
1472 0 : nin1 = ia - after
1473 0 : nout1 = ia - atn
1474 0 : DO ib = 1, before
1475 0 : nin1 = nin1 + after
1476 0 : nin2 = nin1 + atb
1477 0 : nin3 = nin2 + atb
1478 0 : nin4 = nin3 + atb
1479 0 : nout1 = nout1 + atn
1480 0 : nout2 = nout1 + after
1481 0 : nout3 = nout2 + after
1482 0 : nout4 = nout3 + after
1483 0 : DO j = 1, nfft
1484 0 : r1 = zin(1, nin1, j)
1485 0 : s1 = zin(2, nin1, j)
1486 0 : r = zin(1, nin2, j)
1487 0 : s = zin(2, nin2, j)
1488 0 : r2 = r*cr2 - s*ci2
1489 0 : s2 = r*ci2 + s*cr2
1490 0 : r = zin(1, nin3, j)
1491 0 : s = zin(2, nin3, j)
1492 0 : r3 = r*cr3 - s*ci3
1493 0 : s3 = r*ci3 + s*cr3
1494 0 : r = zin(1, nin4, j)
1495 0 : s = zin(2, nin4, j)
1496 0 : r4 = r*cr4 - s*ci4
1497 0 : s4 = r*ci4 + s*cr4
1498 0 : r = r1 + r3
1499 0 : s = r2 + r4
1500 0 : zout(1, j, nout1) = r + s
1501 0 : zout(1, j, nout3) = r - s
1502 0 : r = r1 - r3
1503 0 : s = s2 - s4
1504 0 : zout(1, j, nout2) = r - s
1505 0 : zout(1, j, nout4) = r + s
1506 0 : r = s1 + s3
1507 0 : s = s2 + s4
1508 0 : zout(2, j, nout1) = r + s
1509 0 : zout(2, j, nout3) = r - s
1510 0 : r = s1 - s3
1511 0 : s = r2 - r4
1512 0 : zout(2, j, nout2) = r + s
1513 0 : zout(2, j, nout4) = r - s
1514 : END DO
1515 : END DO
1516 : END IF
1517 : END DO
1518 : ELSE
1519 29406 : ia = 1
1520 29406 : nin1 = ia - after
1521 29406 : nout1 = ia - atn
1522 338748 : DO ib = 1, before
1523 309342 : nin1 = nin1 + after
1524 309342 : nin2 = nin1 + atb
1525 309342 : nin3 = nin2 + atb
1526 309342 : nin4 = nin3 + atb
1527 309342 : nout1 = nout1 + atn
1528 309342 : nout2 = nout1 + after
1529 309342 : nout3 = nout2 + after
1530 309342 : nout4 = nout3 + after
1531 3511836 : DO j = 1, nfft
1532 3173088 : r1 = zin(1, nin1, j)
1533 3173088 : s1 = zin(2, nin1, j)
1534 3173088 : r2 = zin(1, nin2, j)
1535 3173088 : s2 = zin(2, nin2, j)
1536 3173088 : r3 = zin(1, nin3, j)
1537 3173088 : s3 = zin(2, nin3, j)
1538 3173088 : r4 = zin(1, nin4, j)
1539 3173088 : s4 = zin(2, nin4, j)
1540 3173088 : r = r1 + r3
1541 3173088 : s = r2 + r4
1542 3173088 : zout(1, j, nout1) = r + s
1543 3173088 : zout(1, j, nout3) = r - s
1544 3173088 : r = r1 - r3
1545 3173088 : s = s2 - s4
1546 3173088 : zout(1, j, nout2) = r + s
1547 3173088 : zout(1, j, nout4) = r - s
1548 3173088 : r = s1 + s3
1549 3173088 : s = s2 + s4
1550 3173088 : zout(2, j, nout1) = r + s
1551 3173088 : zout(2, j, nout3) = r - s
1552 3173088 : r = s1 - s3
1553 3173088 : s = r2 - r4
1554 3173088 : zout(2, j, nout2) = r - s
1555 3482430 : zout(2, j, nout4) = r + s
1556 : END DO
1557 : END DO
1558 29406 : DO ia = 2, after
1559 0 : ias = ia - 1
1560 29406 : IF (2*ias == after) THEN
1561 0 : nin1 = ia - after
1562 0 : nout1 = ia - atn
1563 0 : DO ib = 1, before
1564 0 : nin1 = nin1 + after
1565 0 : nin2 = nin1 + atb
1566 0 : nin3 = nin2 + atb
1567 0 : nin4 = nin3 + atb
1568 0 : nout1 = nout1 + atn
1569 0 : nout2 = nout1 + after
1570 0 : nout3 = nout2 + after
1571 0 : nout4 = nout3 + after
1572 0 : DO j = 1, nfft
1573 0 : r1 = zin(1, nin1, j)
1574 0 : s1 = zin(2, nin1, j)
1575 0 : r = zin(1, nin2, j)
1576 0 : s = zin(2, nin2, j)
1577 0 : r2 = (r + s)*rt2i
1578 0 : s2 = (s - r)*rt2i
1579 0 : r3 = zin(2, nin3, j)
1580 0 : s3 = -zin(1, nin3, j)
1581 0 : r = zin(1, nin4, j)
1582 0 : s = zin(2, nin4, j)
1583 0 : r4 = (s - r)*rt2i
1584 0 : s4 = -(r + s)*rt2i
1585 0 : r = r1 + r3
1586 0 : s = r2 + r4
1587 0 : zout(1, j, nout1) = r + s
1588 0 : zout(1, j, nout3) = r - s
1589 0 : r = r1 - r3
1590 0 : s = s2 - s4
1591 0 : zout(1, j, nout2) = r + s
1592 0 : zout(1, j, nout4) = r - s
1593 0 : r = s1 + s3
1594 0 : s = s2 + s4
1595 0 : zout(2, j, nout1) = r + s
1596 0 : zout(2, j, nout3) = r - s
1597 0 : r = s1 - s3
1598 0 : s = r2 - r4
1599 0 : zout(2, j, nout2) = r - s
1600 0 : zout(2, j, nout4) = r + s
1601 : END DO
1602 : END DO
1603 : ELSE
1604 0 : itt = ias*before
1605 0 : itrig = itt + 1
1606 0 : cr2 = trig(1, itrig)
1607 0 : ci2 = trig(2, itrig)
1608 0 : itrig = itrig + itt
1609 0 : cr3 = trig(1, itrig)
1610 0 : ci3 = trig(2, itrig)
1611 0 : itrig = itrig + itt
1612 0 : cr4 = trig(1, itrig)
1613 0 : ci4 = trig(2, itrig)
1614 0 : nin1 = ia - after
1615 0 : nout1 = ia - atn
1616 0 : DO ib = 1, before
1617 0 : nin1 = nin1 + after
1618 0 : nin2 = nin1 + atb
1619 0 : nin3 = nin2 + atb
1620 0 : nin4 = nin3 + atb
1621 0 : nout1 = nout1 + atn
1622 0 : nout2 = nout1 + after
1623 0 : nout3 = nout2 + after
1624 0 : nout4 = nout3 + after
1625 0 : DO j = 1, nfft
1626 0 : r1 = zin(1, nin1, j)
1627 0 : s1 = zin(2, nin1, j)
1628 0 : r = zin(1, nin2, j)
1629 0 : s = zin(2, nin2, j)
1630 0 : r2 = r*cr2 - s*ci2
1631 0 : s2 = r*ci2 + s*cr2
1632 0 : r = zin(1, nin3, j)
1633 0 : s = zin(2, nin3, j)
1634 0 : r3 = r*cr3 - s*ci3
1635 0 : s3 = r*ci3 + s*cr3
1636 0 : r = zin(1, nin4, j)
1637 0 : s = zin(2, nin4, j)
1638 0 : r4 = r*cr4 - s*ci4
1639 0 : s4 = r*ci4 + s*cr4
1640 0 : r = r1 + r3
1641 0 : s = r2 + r4
1642 0 : zout(1, j, nout1) = r + s
1643 0 : zout(1, j, nout3) = r - s
1644 0 : r = r1 - r3
1645 0 : s = s2 - s4
1646 0 : zout(1, j, nout2) = r + s
1647 0 : zout(1, j, nout4) = r - s
1648 0 : r = s1 + s3
1649 0 : s = s2 + s4
1650 0 : zout(2, j, nout1) = r + s
1651 0 : zout(2, j, nout3) = r - s
1652 0 : r = s1 - s3
1653 0 : s = r2 - r4
1654 0 : zout(2, j, nout2) = r - s
1655 0 : zout(2, j, nout4) = r + s
1656 : END DO
1657 : END DO
1658 : END IF
1659 : END DO
1660 : END IF
1661 364671 : ELSE IF (now == 8) THEN
1662 90792 : IF (isign == -1) THEN
1663 22296 : ia = 1
1664 22296 : nin1 = ia - after
1665 22296 : nout1 = ia - atn
1666 288336 : DO ib = 1, before
1667 266040 : nin1 = nin1 + after
1668 266040 : nin2 = nin1 + atb
1669 266040 : nin3 = nin2 + atb
1670 266040 : nin4 = nin3 + atb
1671 266040 : nin5 = nin4 + atb
1672 266040 : nin6 = nin5 + atb
1673 266040 : nin7 = nin6 + atb
1674 266040 : nin8 = nin7 + atb
1675 266040 : nout1 = nout1 + atn
1676 266040 : nout2 = nout1 + after
1677 266040 : nout3 = nout2 + after
1678 266040 : nout4 = nout3 + after
1679 266040 : nout5 = nout4 + after
1680 266040 : nout6 = nout5 + after
1681 266040 : nout7 = nout6 + after
1682 266040 : nout8 = nout7 + after
1683 1625808 : DO j = 1, nfft
1684 1337472 : r1 = zin(1, nin1, j)
1685 1337472 : s1 = zin(2, nin1, j)
1686 1337472 : r2 = zin(1, nin2, j)
1687 1337472 : s2 = zin(2, nin2, j)
1688 1337472 : r3 = zin(1, nin3, j)
1689 1337472 : s3 = zin(2, nin3, j)
1690 1337472 : r4 = zin(1, nin4, j)
1691 1337472 : s4 = zin(2, nin4, j)
1692 1337472 : r5 = zin(1, nin5, j)
1693 1337472 : s5 = zin(2, nin5, j)
1694 1337472 : r6 = zin(1, nin6, j)
1695 1337472 : s6 = zin(2, nin6, j)
1696 1337472 : r7 = zin(1, nin7, j)
1697 1337472 : s7 = zin(2, nin7, j)
1698 1337472 : r8 = zin(1, nin8, j)
1699 1337472 : s8 = zin(2, nin8, j)
1700 1337472 : r = r1 + r5
1701 1337472 : s = r3 + r7
1702 1337472 : ap = r + s
1703 1337472 : am = r - s
1704 1337472 : r = r2 + r6
1705 1337472 : s = r4 + r8
1706 1337472 : bp = r + s
1707 1337472 : bm = r - s
1708 1337472 : r = s1 + s5
1709 1337472 : s = s3 + s7
1710 1337472 : cp = r + s
1711 1337472 : cm = r - s
1712 1337472 : r = s2 + s6
1713 1337472 : s = s4 + s8
1714 1337472 : dbl = r + s
1715 1337472 : dm = r - s
1716 1337472 : zout(1, j, nout1) = ap + bp
1717 1337472 : zout(2, j, nout1) = cp + dbl
1718 1337472 : zout(1, j, nout5) = ap - bp
1719 1337472 : zout(2, j, nout5) = cp - dbl
1720 1337472 : zout(1, j, nout3) = am + dm
1721 1337472 : zout(2, j, nout3) = cm - bm
1722 1337472 : zout(1, j, nout7) = am - dm
1723 1337472 : zout(2, j, nout7) = cm + bm
1724 1337472 : r = r1 - r5
1725 1337472 : s = s3 - s7
1726 1337472 : ap = r + s
1727 1337472 : am = r - s
1728 1337472 : r = s1 - s5
1729 1337472 : s = r3 - r7
1730 1337472 : bp = r + s
1731 1337472 : bm = r - s
1732 1337472 : r = s4 - s8
1733 1337472 : s = r2 - r6
1734 1337472 : cp = r + s
1735 1337472 : cm = r - s
1736 1337472 : r = s2 - s6
1737 1337472 : s = r4 - r8
1738 1337472 : dbl = r + s
1739 1337472 : dm = r - s
1740 1337472 : r = (cp + dm)*rt2i
1741 1337472 : s = (-cp + dm)*rt2i
1742 1337472 : cp = (cm + dbl)*rt2i
1743 1337472 : dbl = (cm - dbl)*rt2i
1744 1337472 : zout(1, j, nout2) = ap + r
1745 1337472 : zout(2, j, nout2) = bm + s
1746 1337472 : zout(1, j, nout6) = ap - r
1747 1337472 : zout(2, j, nout6) = bm - s
1748 1337472 : zout(1, j, nout4) = am + cp
1749 1337472 : zout(2, j, nout4) = bp + dbl
1750 1337472 : zout(1, j, nout8) = am - cp
1751 1603512 : zout(2, j, nout8) = bp - dbl
1752 : END DO
1753 : END DO
1754 : ELSE
1755 68496 : ia = 1
1756 68496 : nin1 = ia - after
1757 68496 : nout1 = ia - atn
1758 519336 : DO ib = 1, before
1759 450840 : nin1 = nin1 + after
1760 450840 : nin2 = nin1 + atb
1761 450840 : nin3 = nin2 + atb
1762 450840 : nin4 = nin3 + atb
1763 450840 : nin5 = nin4 + atb
1764 450840 : nin6 = nin5 + atb
1765 450840 : nin7 = nin6 + atb
1766 450840 : nin8 = nin7 + atb
1767 450840 : nout1 = nout1 + atn
1768 450840 : nout2 = nout1 + after
1769 450840 : nout3 = nout2 + after
1770 450840 : nout4 = nout3 + after
1771 450840 : nout5 = nout4 + after
1772 450840 : nout6 = nout5 + after
1773 450840 : nout7 = nout6 + after
1774 450840 : nout8 = nout7 + after
1775 4560168 : DO j = 1, nfft
1776 4040832 : r1 = zin(1, nin1, j)
1777 4040832 : s1 = zin(2, nin1, j)
1778 4040832 : r2 = zin(1, nin2, j)
1779 4040832 : s2 = zin(2, nin2, j)
1780 4040832 : r3 = zin(1, nin3, j)
1781 4040832 : s3 = zin(2, nin3, j)
1782 4040832 : r4 = zin(1, nin4, j)
1783 4040832 : s4 = zin(2, nin4, j)
1784 4040832 : r5 = zin(1, nin5, j)
1785 4040832 : s5 = zin(2, nin5, j)
1786 4040832 : r6 = zin(1, nin6, j)
1787 4040832 : s6 = zin(2, nin6, j)
1788 4040832 : r7 = zin(1, nin7, j)
1789 4040832 : s7 = zin(2, nin7, j)
1790 4040832 : r8 = zin(1, nin8, j)
1791 4040832 : s8 = zin(2, nin8, j)
1792 4040832 : r = r1 + r5
1793 4040832 : s = r3 + r7
1794 4040832 : ap = r + s
1795 4040832 : am = r - s
1796 4040832 : r = r2 + r6
1797 4040832 : s = r4 + r8
1798 4040832 : bp = r + s
1799 4040832 : bm = r - s
1800 4040832 : r = s1 + s5
1801 4040832 : s = s3 + s7
1802 4040832 : cp = r + s
1803 4040832 : cm = r - s
1804 4040832 : r = s2 + s6
1805 4040832 : s = s4 + s8
1806 4040832 : dbl = r + s
1807 4040832 : dm = r - s
1808 4040832 : zout(1, j, nout1) = ap + bp
1809 4040832 : zout(2, j, nout1) = cp + dbl
1810 4040832 : zout(1, j, nout5) = ap - bp
1811 4040832 : zout(2, j, nout5) = cp - dbl
1812 4040832 : zout(1, j, nout3) = am - dm
1813 4040832 : zout(2, j, nout3) = cm + bm
1814 4040832 : zout(1, j, nout7) = am + dm
1815 4040832 : zout(2, j, nout7) = cm - bm
1816 4040832 : r = r1 - r5
1817 4040832 : s = -s3 + s7
1818 4040832 : ap = r + s
1819 4040832 : am = r - s
1820 4040832 : r = s1 - s5
1821 4040832 : s = r7 - r3
1822 4040832 : bp = r + s
1823 4040832 : bm = r - s
1824 4040832 : r = -s4 + s8
1825 4040832 : s = r2 - r6
1826 4040832 : cp = r + s
1827 4040832 : cm = r - s
1828 4040832 : r = -s2 + s6
1829 4040832 : s = r4 - r8
1830 4040832 : dbl = r + s
1831 4040832 : dm = r - s
1832 4040832 : r = (cp + dm)*rt2i
1833 4040832 : s = (cp - dm)*rt2i
1834 4040832 : cp = (cm + dbl)*rt2i
1835 4040832 : dbl = (-cm + dbl)*rt2i
1836 4040832 : zout(1, j, nout2) = ap + r
1837 4040832 : zout(2, j, nout2) = bm + s
1838 4040832 : zout(1, j, nout6) = ap - r
1839 4040832 : zout(2, j, nout6) = bm - s
1840 4040832 : zout(1, j, nout4) = am + cp
1841 4040832 : zout(2, j, nout4) = bp + dbl
1842 4040832 : zout(1, j, nout8) = am - cp
1843 4491672 : zout(2, j, nout8) = bp - dbl
1844 : END DO
1845 : END DO
1846 : END IF
1847 273879 : ELSE IF (now == 3) THEN
1848 2430 : ia = 1
1849 2430 : nin1 = ia - after
1850 2430 : nout1 = ia - atn
1851 2430 : bbs = isign*bb
1852 22392 : DO ib = 1, before
1853 19962 : nin1 = nin1 + after
1854 19962 : nin2 = nin1 + atb
1855 19962 : nin3 = nin2 + atb
1856 19962 : nout1 = nout1 + atn
1857 19962 : nout2 = nout1 + after
1858 19962 : nout3 = nout2 + after
1859 375957 : DO j = 1, nfft
1860 353565 : r1 = zin(1, nin1, j)
1861 353565 : s1 = zin(2, nin1, j)
1862 353565 : r2 = zin(1, nin2, j)
1863 353565 : s2 = zin(2, nin2, j)
1864 353565 : r3 = zin(1, nin3, j)
1865 353565 : s3 = zin(2, nin3, j)
1866 353565 : r = r2 + r3
1867 353565 : s = s2 + s3
1868 353565 : zout(1, j, nout1) = r + r1
1869 353565 : zout(2, j, nout1) = s + s1
1870 353565 : r1 = r1 - 0.5_dp*r
1871 353565 : s1 = s1 - 0.5_dp*s
1872 353565 : r2 = bbs*(r2 - r3)
1873 353565 : s2 = bbs*(s2 - s3)
1874 353565 : zout(1, j, nout2) = r1 - s2
1875 353565 : zout(2, j, nout2) = s1 + r2
1876 353565 : zout(1, j, nout3) = r1 + s2
1877 373527 : zout(2, j, nout3) = s1 - r2
1878 : END DO
1879 : END DO
1880 2430 : DO ia = 2, after
1881 0 : ias = ia - 1
1882 2430 : IF (4*ias == 3*after) THEN
1883 0 : IF (isign == 1) THEN
1884 0 : nin1 = ia - after
1885 0 : nout1 = ia - atn
1886 0 : DO ib = 1, before
1887 0 : nin1 = nin1 + after
1888 0 : nin2 = nin1 + atb
1889 0 : nin3 = nin2 + atb
1890 0 : nout1 = nout1 + atn
1891 0 : nout2 = nout1 + after
1892 0 : nout3 = nout2 + after
1893 0 : DO j = 1, nfft
1894 0 : r1 = zin(1, nin1, j)
1895 0 : s1 = zin(2, nin1, j)
1896 0 : r2 = -zin(2, nin2, j)
1897 0 : s2 = zin(1, nin2, j)
1898 0 : r3 = -zin(1, nin3, j)
1899 0 : s3 = -zin(2, nin3, j)
1900 0 : r = r2 + r3
1901 0 : s = s2 + s3
1902 0 : zout(1, j, nout1) = r + r1
1903 0 : zout(2, j, nout1) = s + s1
1904 0 : r1 = r1 - 0.5_dp*r
1905 0 : s1 = s1 - 0.5_dp*s
1906 0 : r2 = bbs*(r2 - r3)
1907 0 : s2 = bbs*(s2 - s3)
1908 0 : zout(1, j, nout2) = r1 - s2
1909 0 : zout(2, j, nout2) = s1 + r2
1910 0 : zout(1, j, nout3) = r1 + s2
1911 0 : zout(2, j, nout3) = s1 - r2
1912 : END DO
1913 : END DO
1914 : ELSE
1915 0 : nin1 = ia - after
1916 0 : nout1 = ia - atn
1917 0 : DO ib = 1, before
1918 0 : nin1 = nin1 + after
1919 0 : nin2 = nin1 + atb
1920 0 : nin3 = nin2 + atb
1921 0 : nout1 = nout1 + atn
1922 0 : nout2 = nout1 + after
1923 0 : nout3 = nout2 + after
1924 0 : DO j = 1, nfft
1925 0 : r1 = zin(1, nin1, j)
1926 0 : s1 = zin(2, nin1, j)
1927 0 : r2 = zin(2, nin2, j)
1928 0 : s2 = -zin(1, nin2, j)
1929 0 : r3 = -zin(1, nin3, j)
1930 0 : s3 = -zin(2, nin3, j)
1931 0 : r = r2 + r3
1932 0 : s = s2 + s3
1933 0 : zout(1, j, nout1) = r + r1
1934 0 : zout(2, j, nout1) = s + s1
1935 0 : r1 = r1 - 0.5_dp*r
1936 0 : s1 = s1 - 0.5_dp*s
1937 0 : r2 = bbs*(r2 - r3)
1938 0 : s2 = bbs*(s2 - s3)
1939 0 : zout(1, j, nout2) = r1 - s2
1940 0 : zout(2, j, nout2) = s1 + r2
1941 0 : zout(1, j, nout3) = r1 + s2
1942 0 : zout(2, j, nout3) = s1 - r2
1943 : END DO
1944 : END DO
1945 : END IF
1946 0 : ELSE IF (8*ias == 3*after) THEN
1947 0 : IF (isign == 1) THEN
1948 0 : nin1 = ia - after
1949 0 : nout1 = ia - atn
1950 0 : DO ib = 1, before
1951 0 : nin1 = nin1 + after
1952 0 : nin2 = nin1 + atb
1953 0 : nin3 = nin2 + atb
1954 0 : nout1 = nout1 + atn
1955 0 : nout2 = nout1 + after
1956 0 : nout3 = nout2 + after
1957 0 : DO j = 1, nfft
1958 0 : r1 = zin(1, nin1, j)
1959 0 : s1 = zin(2, nin1, j)
1960 0 : r = zin(1, nin2, j)
1961 0 : s = zin(2, nin2, j)
1962 0 : r2 = (r - s)*rt2i
1963 0 : s2 = (r + s)*rt2i
1964 0 : r3 = -zin(2, nin3, j)
1965 0 : s3 = zin(1, nin3, j)
1966 0 : r = r2 + r3
1967 0 : s = s2 + s3
1968 0 : zout(1, j, nout1) = r + r1
1969 0 : zout(2, j, nout1) = s + s1
1970 0 : r1 = r1 - 0.5_dp*r
1971 0 : s1 = s1 - 0.5_dp*s
1972 0 : r2 = bbs*(r2 - r3)
1973 0 : s2 = bbs*(s2 - s3)
1974 0 : zout(1, j, nout2) = r1 - s2
1975 0 : zout(2, j, nout2) = s1 + r2
1976 0 : zout(1, j, nout3) = r1 + s2
1977 0 : zout(2, j, nout3) = s1 - r2
1978 : END DO
1979 : END DO
1980 : ELSE
1981 0 : nin1 = ia - after
1982 0 : nout1 = ia - atn
1983 0 : DO ib = 1, before
1984 0 : nin1 = nin1 + after
1985 0 : nin2 = nin1 + atb
1986 0 : nin3 = nin2 + atb
1987 0 : nout1 = nout1 + atn
1988 0 : nout2 = nout1 + after
1989 0 : nout3 = nout2 + after
1990 0 : DO j = 1, nfft
1991 0 : r1 = zin(1, nin1, j)
1992 0 : s1 = zin(2, nin1, j)
1993 0 : r = zin(1, nin2, j)
1994 0 : s = zin(2, nin2, j)
1995 0 : r2 = (r + s)*rt2i
1996 0 : s2 = (-r + s)*rt2i
1997 0 : r3 = zin(2, nin3, j)
1998 0 : s3 = -zin(1, nin3, j)
1999 0 : r = r2 + r3
2000 0 : s = s2 + s3
2001 0 : zout(1, j, nout1) = r + r1
2002 0 : zout(2, j, nout1) = s + s1
2003 0 : r1 = r1 - 0.5_dp*r
2004 0 : s1 = s1 - 0.5_dp*s
2005 0 : r2 = bbs*(r2 - r3)
2006 0 : s2 = bbs*(s2 - s3)
2007 0 : zout(1, j, nout2) = r1 - s2
2008 0 : zout(2, j, nout2) = s1 + r2
2009 0 : zout(1, j, nout3) = r1 + s2
2010 0 : zout(2, j, nout3) = s1 - r2
2011 : END DO
2012 : END DO
2013 : END IF
2014 : ELSE
2015 0 : itt = ias*before
2016 0 : itrig = itt + 1
2017 0 : cr2 = trig(1, itrig)
2018 0 : ci2 = trig(2, itrig)
2019 0 : itrig = itrig + itt
2020 0 : cr3 = trig(1, itrig)
2021 0 : ci3 = trig(2, itrig)
2022 0 : nin1 = ia - after
2023 0 : nout1 = ia - atn
2024 0 : DO ib = 1, before
2025 0 : nin1 = nin1 + after
2026 0 : nin2 = nin1 + atb
2027 0 : nin3 = nin2 + atb
2028 0 : nout1 = nout1 + atn
2029 0 : nout2 = nout1 + after
2030 0 : nout3 = nout2 + after
2031 0 : DO j = 1, nfft
2032 0 : r1 = zin(1, nin1, j)
2033 0 : s1 = zin(2, nin1, j)
2034 0 : r = zin(1, nin2, j)
2035 0 : s = zin(2, nin2, j)
2036 0 : r2 = r*cr2 - s*ci2
2037 0 : s2 = r*ci2 + s*cr2
2038 0 : r = zin(1, nin3, j)
2039 0 : s = zin(2, nin3, j)
2040 0 : r3 = r*cr3 - s*ci3
2041 0 : s3 = r*ci3 + s*cr3
2042 0 : r = r2 + r3
2043 0 : s = s2 + s3
2044 0 : zout(1, j, nout1) = r + r1
2045 0 : zout(2, j, nout1) = s + s1
2046 0 : r1 = r1 - 0.5_dp*r
2047 0 : s1 = s1 - 0.5_dp*s
2048 0 : r2 = bbs*(r2 - r3)
2049 0 : s2 = bbs*(s2 - s3)
2050 0 : zout(1, j, nout2) = r1 - s2
2051 0 : zout(2, j, nout2) = s1 + r2
2052 0 : zout(1, j, nout3) = r1 + s2
2053 0 : zout(2, j, nout3) = s1 - r2
2054 : END DO
2055 : END DO
2056 : END IF
2057 : END DO
2058 271449 : ELSE IF (now == 5) THEN
2059 94863 : sin2 = isign*sin2p
2060 94863 : sin4 = isign*sin4p
2061 94863 : ia = 1
2062 94863 : nin1 = ia - after
2063 94863 : nout1 = ia - atn
2064 891189 : DO ib = 1, before
2065 796326 : nin1 = nin1 + after
2066 796326 : nin2 = nin1 + atb
2067 796326 : nin3 = nin2 + atb
2068 796326 : nin4 = nin3 + atb
2069 796326 : nin5 = nin4 + atb
2070 796326 : nout1 = nout1 + atn
2071 796326 : nout2 = nout1 + after
2072 796326 : nout3 = nout2 + after
2073 796326 : nout4 = nout3 + after
2074 796326 : nout5 = nout4 + after
2075 9221919 : DO j = 1, nfft
2076 8330730 : r1 = zin(1, nin1, j)
2077 8330730 : s1 = zin(2, nin1, j)
2078 8330730 : r2 = zin(1, nin2, j)
2079 8330730 : s2 = zin(2, nin2, j)
2080 8330730 : r3 = zin(1, nin3, j)
2081 8330730 : s3 = zin(2, nin3, j)
2082 8330730 : r4 = zin(1, nin4, j)
2083 8330730 : s4 = zin(2, nin4, j)
2084 8330730 : r5 = zin(1, nin5, j)
2085 8330730 : s5 = zin(2, nin5, j)
2086 8330730 : r25 = r2 + r5
2087 8330730 : r34 = r3 + r4
2088 8330730 : s25 = s2 - s5
2089 8330730 : s34 = s3 - s4
2090 8330730 : zout(1, j, nout1) = r1 + r25 + r34
2091 8330730 : r = cos2*r25 + cos4*r34 + r1
2092 8330730 : s = sin2*s25 + sin4*s34
2093 8330730 : zout(1, j, nout2) = r - s
2094 8330730 : zout(1, j, nout5) = r + s
2095 8330730 : r = cos4*r25 + cos2*r34 + r1
2096 8330730 : s = sin4*s25 - sin2*s34
2097 8330730 : zout(1, j, nout3) = r - s
2098 8330730 : zout(1, j, nout4) = r + s
2099 8330730 : r25 = r2 - r5
2100 8330730 : r34 = r3 - r4
2101 8330730 : s25 = s2 + s5
2102 8330730 : s34 = s3 + s4
2103 8330730 : zout(2, j, nout1) = s1 + s25 + s34
2104 8330730 : r = cos2*s25 + cos4*s34 + s1
2105 8330730 : s = sin2*r25 + sin4*r34
2106 8330730 : zout(2, j, nout2) = r + s
2107 8330730 : zout(2, j, nout5) = r - s
2108 8330730 : r = cos4*s25 + cos2*s34 + s1
2109 8330730 : s = sin4*r25 - sin2*r34
2110 8330730 : zout(2, j, nout3) = r + s
2111 9127056 : zout(2, j, nout4) = r - s
2112 : END DO
2113 : END DO
2114 94863 : DO ia = 2, after
2115 0 : ias = ia - 1
2116 94863 : IF (8*ias == 5*after) THEN
2117 0 : IF (isign == 1) THEN
2118 0 : nin1 = ia - after
2119 0 : nout1 = ia - atn
2120 0 : DO ib = 1, before
2121 0 : nin1 = nin1 + after
2122 0 : nin2 = nin1 + atb
2123 0 : nin3 = nin2 + atb
2124 0 : nin4 = nin3 + atb
2125 0 : nin5 = nin4 + atb
2126 0 : nout1 = nout1 + atn
2127 0 : nout2 = nout1 + after
2128 0 : nout3 = nout2 + after
2129 0 : nout4 = nout3 + after
2130 0 : nout5 = nout4 + after
2131 0 : DO j = 1, nfft
2132 0 : r1 = zin(1, nin1, j)
2133 0 : s1 = zin(2, nin1, j)
2134 0 : r = zin(1, nin2, j)
2135 0 : s = zin(2, nin2, j)
2136 0 : r2 = (r - s)*rt2i
2137 0 : s2 = (r + s)*rt2i
2138 0 : r3 = -zin(2, nin3, j)
2139 0 : s3 = zin(1, nin3, j)
2140 0 : r = zin(1, nin4, j)
2141 0 : s = zin(2, nin4, j)
2142 0 : r4 = -(r + s)*rt2i
2143 0 : s4 = (r - s)*rt2i
2144 0 : r5 = -zin(1, nin5, j)
2145 0 : s5 = -zin(2, nin5, j)
2146 0 : r25 = r2 + r5
2147 0 : r34 = r3 + r4
2148 0 : s25 = s2 - s5
2149 0 : s34 = s3 - s4
2150 0 : zout(1, j, nout1) = r1 + r25 + r34
2151 0 : r = cos2*r25 + cos4*r34 + r1
2152 0 : s = sin2*s25 + sin4*s34
2153 0 : zout(1, j, nout2) = r - s
2154 0 : zout(1, j, nout5) = r + s
2155 0 : r = cos4*r25 + cos2*r34 + r1
2156 0 : s = sin4*s25 - sin2*s34
2157 0 : zout(1, j, nout3) = r - s
2158 0 : zout(1, j, nout4) = r + s
2159 0 : r25 = r2 - r5
2160 0 : r34 = r3 - r4
2161 0 : s25 = s2 + s5
2162 0 : s34 = s3 + s4
2163 0 : zout(2, j, nout1) = s1 + s25 + s34
2164 0 : r = cos2*s25 + cos4*s34 + s1
2165 0 : s = sin2*r25 + sin4*r34
2166 0 : zout(2, j, nout2) = r + s
2167 0 : zout(2, j, nout5) = r - s
2168 0 : r = cos4*s25 + cos2*s34 + s1
2169 0 : s = sin4*r25 - sin2*r34
2170 0 : zout(2, j, nout3) = r + s
2171 0 : zout(2, j, nout4) = r - s
2172 : END DO
2173 : END DO
2174 : ELSE
2175 0 : nin1 = ia - after
2176 0 : nout1 = ia - atn
2177 0 : DO ib = 1, before
2178 0 : nin1 = nin1 + after
2179 0 : nin2 = nin1 + atb
2180 0 : nin3 = nin2 + atb
2181 0 : nin4 = nin3 + atb
2182 0 : nin5 = nin4 + atb
2183 0 : nout1 = nout1 + atn
2184 0 : nout2 = nout1 + after
2185 0 : nout3 = nout2 + after
2186 0 : nout4 = nout3 + after
2187 0 : nout5 = nout4 + after
2188 0 : DO j = 1, nfft
2189 0 : r1 = zin(1, nin1, j)
2190 0 : s1 = zin(2, nin1, j)
2191 0 : r = zin(1, nin2, j)
2192 0 : s = zin(2, nin2, j)
2193 0 : r2 = (r + s)*rt2i
2194 0 : s2 = (-r + s)*rt2i
2195 0 : r3 = zin(2, nin3, j)
2196 0 : s3 = -zin(1, nin3, j)
2197 0 : r = zin(1, nin4, j)
2198 0 : s = zin(2, nin4, j)
2199 0 : r4 = (s - r)*rt2i
2200 0 : s4 = -(r + s)*rt2i
2201 0 : r5 = -zin(1, nin5, j)
2202 0 : s5 = -zin(2, nin5, j)
2203 0 : r25 = r2 + r5
2204 0 : r34 = r3 + r4
2205 0 : s25 = s2 - s5
2206 0 : s34 = s3 - s4
2207 0 : zout(1, j, nout1) = r1 + r25 + r34
2208 0 : r = cos2*r25 + cos4*r34 + r1
2209 0 : s = sin2*s25 + sin4*s34
2210 0 : zout(1, j, nout2) = r - s
2211 0 : zout(1, j, nout5) = r + s
2212 0 : r = cos4*r25 + cos2*r34 + r1
2213 0 : s = sin4*s25 - sin2*s34
2214 0 : zout(1, j, nout3) = r - s
2215 0 : zout(1, j, nout4) = r + s
2216 0 : r25 = r2 - r5
2217 0 : r34 = r3 - r4
2218 0 : s25 = s2 + s5
2219 0 : s34 = s3 + s4
2220 0 : zout(2, j, nout1) = s1 + s25 + s34
2221 0 : r = cos2*s25 + cos4*s34 + s1
2222 0 : s = sin2*r25 + sin4*r34
2223 0 : zout(2, j, nout2) = r + s
2224 0 : zout(2, j, nout5) = r - s
2225 0 : r = cos4*s25 + cos2*s34 + s1
2226 0 : s = sin4*r25 - sin2*r34
2227 0 : zout(2, j, nout3) = r + s
2228 0 : zout(2, j, nout4) = r - s
2229 : END DO
2230 : END DO
2231 : END IF
2232 : ELSE
2233 0 : ias = ia - 1
2234 0 : itt = ias*before
2235 0 : itrig = itt + 1
2236 0 : cr2 = trig(1, itrig)
2237 0 : ci2 = trig(2, itrig)
2238 0 : itrig = itrig + itt
2239 0 : cr3 = trig(1, itrig)
2240 0 : ci3 = trig(2, itrig)
2241 0 : itrig = itrig + itt
2242 0 : cr4 = trig(1, itrig)
2243 0 : ci4 = trig(2, itrig)
2244 0 : itrig = itrig + itt
2245 0 : cr5 = trig(1, itrig)
2246 0 : ci5 = trig(2, itrig)
2247 0 : nin1 = ia - after
2248 0 : nout1 = ia - atn
2249 0 : DO ib = 1, before
2250 0 : nin1 = nin1 + after
2251 0 : nin2 = nin1 + atb
2252 0 : nin3 = nin2 + atb
2253 0 : nin4 = nin3 + atb
2254 0 : nin5 = nin4 + atb
2255 0 : nout1 = nout1 + atn
2256 0 : nout2 = nout1 + after
2257 0 : nout3 = nout2 + after
2258 0 : nout4 = nout3 + after
2259 0 : nout5 = nout4 + after
2260 0 : DO j = 1, nfft
2261 0 : r1 = zin(1, nin1, j)
2262 0 : s1 = zin(2, nin1, j)
2263 0 : r = zin(1, nin2, j)
2264 0 : s = zin(2, nin2, j)
2265 0 : r2 = r*cr2 - s*ci2
2266 0 : s2 = r*ci2 + s*cr2
2267 0 : r = zin(1, nin3, j)
2268 0 : s = zin(2, nin3, j)
2269 0 : r3 = r*cr3 - s*ci3
2270 0 : s3 = r*ci3 + s*cr3
2271 0 : r = zin(1, nin4, j)
2272 0 : s = zin(2, nin4, j)
2273 0 : r4 = r*cr4 - s*ci4
2274 0 : s4 = r*ci4 + s*cr4
2275 0 : r = zin(1, nin5, j)
2276 0 : s = zin(2, nin5, j)
2277 0 : r5 = r*cr5 - s*ci5
2278 0 : s5 = r*ci5 + s*cr5
2279 0 : r25 = r2 + r5
2280 0 : r34 = r3 + r4
2281 0 : s25 = s2 - s5
2282 0 : s34 = s3 - s4
2283 0 : zout(1, j, nout1) = r1 + r25 + r34
2284 0 : r = cos2*r25 + cos4*r34 + r1
2285 0 : s = sin2*s25 + sin4*s34
2286 0 : zout(1, j, nout2) = r - s
2287 0 : zout(1, j, nout5) = r + s
2288 0 : r = cos4*r25 + cos2*r34 + r1
2289 0 : s = sin4*s25 - sin2*s34
2290 0 : zout(1, j, nout3) = r - s
2291 0 : zout(1, j, nout4) = r + s
2292 0 : r25 = r2 - r5
2293 0 : r34 = r3 - r4
2294 0 : s25 = s2 + s5
2295 0 : s34 = s3 + s4
2296 0 : zout(2, j, nout1) = s1 + s25 + s34
2297 0 : r = cos2*s25 + cos4*s34 + s1
2298 0 : s = sin2*r25 + sin4*r34
2299 0 : zout(2, j, nout2) = r + s
2300 0 : zout(2, j, nout5) = r - s
2301 0 : r = cos4*s25 + cos2*s34 + s1
2302 0 : s = sin4*r25 - sin2*r34
2303 0 : zout(2, j, nout3) = r + s
2304 0 : zout(2, j, nout4) = r - s
2305 : END DO
2306 : END DO
2307 : END IF
2308 : END DO
2309 176586 : ELSE IF (now == 6) THEN
2310 176586 : bbs = isign*bb
2311 176586 : ia = 1
2312 176586 : nin1 = ia - after
2313 176586 : nout1 = ia - atn
2314 2663940 : DO ib = 1, before
2315 2487354 : nin1 = nin1 + after
2316 2487354 : nin2 = nin1 + atb
2317 2487354 : nin3 = nin2 + atb
2318 2487354 : nin4 = nin3 + atb
2319 2487354 : nin5 = nin4 + atb
2320 2487354 : nin6 = nin5 + atb
2321 2487354 : nout1 = nout1 + atn
2322 2487354 : nout2 = nout1 + after
2323 2487354 : nout3 = nout2 + after
2324 2487354 : nout4 = nout3 + after
2325 2487354 : nout5 = nout4 + after
2326 2487354 : nout6 = nout5 + after
2327 15972672 : DO j = 1, nfft
2328 13308732 : r2 = zin(1, nin3, j)
2329 13308732 : s2 = zin(2, nin3, j)
2330 13308732 : r3 = zin(1, nin5, j)
2331 13308732 : s3 = zin(2, nin5, j)
2332 13308732 : r = r2 + r3
2333 13308732 : s = s2 + s3
2334 13308732 : r1 = zin(1, nin1, j)
2335 13308732 : s1 = zin(2, nin1, j)
2336 13308732 : ur1 = r + r1
2337 13308732 : ui1 = s + s1
2338 13308732 : r1 = r1 - 0.5_dp*r
2339 13308732 : s1 = s1 - 0.5_dp*s
2340 13308732 : r = r2 - r3
2341 13308732 : s = s2 - s3
2342 13308732 : ur2 = r1 - s*bbs
2343 13308732 : ui2 = s1 + r*bbs
2344 13308732 : ur3 = r1 + s*bbs
2345 13308732 : ui3 = s1 - r*bbs
2346 :
2347 13308732 : r2 = zin(1, nin6, j)
2348 13308732 : s2 = zin(2, nin6, j)
2349 13308732 : r3 = zin(1, nin2, j)
2350 13308732 : s3 = zin(2, nin2, j)
2351 13308732 : r = r2 + r3
2352 13308732 : s = s2 + s3
2353 13308732 : r1 = zin(1, nin4, j)
2354 13308732 : s1 = zin(2, nin4, j)
2355 13308732 : vr1 = r + r1
2356 13308732 : vi1 = s + s1
2357 13308732 : r1 = r1 - 0.5_dp*r
2358 13308732 : s1 = s1 - 0.5_dp*s
2359 13308732 : r = r2 - r3
2360 13308732 : s = s2 - s3
2361 13308732 : vr2 = r1 - s*bbs
2362 13308732 : vi2 = s1 + r*bbs
2363 13308732 : vr3 = r1 + s*bbs
2364 13308732 : vi3 = s1 - r*bbs
2365 :
2366 13308732 : zout(1, j, nout1) = ur1 + vr1
2367 13308732 : zout(2, j, nout1) = ui1 + vi1
2368 13308732 : zout(1, j, nout5) = ur2 + vr2
2369 13308732 : zout(2, j, nout5) = ui2 + vi2
2370 13308732 : zout(1, j, nout3) = ur3 + vr3
2371 13308732 : zout(2, j, nout3) = ui3 + vi3
2372 13308732 : zout(1, j, nout4) = ur1 - vr1
2373 13308732 : zout(2, j, nout4) = ui1 - vi1
2374 13308732 : zout(1, j, nout2) = ur2 - vr2
2375 13308732 : zout(2, j, nout2) = ui2 - vi2
2376 13308732 : zout(1, j, nout6) = ur3 - vr3
2377 15796086 : zout(2, j, nout6) = ui3 - vi3
2378 : END DO
2379 : END DO
2380 : ELSE
2381 0 : CPABORT('Error fftpre')
2382 : END IF
2383 :
2384 : !-----------------------------------------------------------------------------!
2385 :
2386 424773 : END SUBROUTINE fftpre
2387 :
2388 : !-----------------------------------------------------------------------------!
2389 :
2390 : !-----------------------------------------------------------------------------!
2391 : ! Copyright (C) Stefan Goedecker, Lausanne, Switzerland, August 1, 1991
2392 : ! Copyright (C) Stefan Goedecker, Cornell University, Ithaca, USA, 1994
2393 : ! Copyright (C) Stefan Goedecker, MPI Stuttgart, Germany, 1999
2394 : ! This file is distributed under the terms of the
2395 : ! GNU General Public License version 2 (or later),
2396 : ! see http://www.gnu.org/copyleft/gpl.txt .
2397 : !-----------------------------------------------------------------------------!
2398 : ! S. Goedecker: Rotating a three-dimensional array in optimal
2399 : ! positions for vector processing: Case study for a three-dimensional Fast
2400 : ! Fourier Transform, Comp. Phys. Commun. 76, 294 (1993)
2401 : ! **************************************************************************************************
2402 : !> \brief ...
2403 : !> \param mm ...
2404 : !> \param nfft ...
2405 : !> \param m ...
2406 : !> \param nn ...
2407 : !> \param n ...
2408 : !> \param zin ...
2409 : !> \param zout ...
2410 : !> \param trig ...
2411 : !> \param now ...
2412 : !> \param after ...
2413 : !> \param before ...
2414 : !> \param isign ...
2415 : ! **************************************************************************************************
2416 1244146 : SUBROUTINE fftstp(mm, nfft, m, nn, n, zin, zout, trig, now, after, before, isign)
2417 :
2418 : INTEGER, INTENT(IN) :: mm, nfft, m, nn, n
2419 : REAL(dp), DIMENSION(2, mm, m), INTENT(IN) :: zin
2420 : REAL(dp), DIMENSION(2, nn, n), INTENT(INOUT) :: zout
2421 : REAL(dp), DIMENSION(2, ctrig_length), INTENT(IN) :: trig
2422 : INTEGER, INTENT(IN) :: now, after, before, isign
2423 :
2424 : REAL(dp), PARAMETER :: bb = 0.8660254037844387_dp, cos2 = 0.3090169943749474_dp, &
2425 : cos4 = -0.8090169943749474_dp, rt2i = 0.7071067811865475_dp, &
2426 : sin2p = 0.9510565162951536_dp, sin4p = 0.5877852522924731_dp
2427 :
2428 : INTEGER :: atb, atn, ia, ias, ib, itrig, itt, j, &
2429 : nin1, nin2, nin3, nin4, nin5, nin6, &
2430 : nin7, nin8, nout1, nout2, nout3, &
2431 : nout4, nout5, nout6, nout7, nout8
2432 : REAL(dp) :: am, ap, bbs, bm, bp, ci2, ci3, ci4, ci5, cm, cp, cr2, cr3, cr4, cr5, dbl, dm, r, &
2433 : r1, r2, r25, r3, r34, r4, r5, r6, r7, r8, s, s1, s2, s25, s3, s34, s4, s5, s6, s7, s8, &
2434 : sin2, sin4, ui1, ui2, ui3, ur1, ur2, ur3, vi1, vi2, vi3, vr1, vr2, vr3
2435 :
2436 : ! sqrt(0.5)
2437 : ! sqrt(3)/2
2438 : ! cos(2*pi/5)
2439 : ! cos(4*pi/5)
2440 : ! sin(2*pi/5)
2441 : ! sin(4*pi/5)
2442 : !-----------------------------------------------------------------------------!
2443 :
2444 1244146 : atn = after*now
2445 1244146 : atb = after*before
2446 :
2447 1244146 : IF (now == 4) THEN
2448 202748 : IF (isign == 1) THEN
2449 139269 : ia = 1
2450 139269 : nin1 = ia - after
2451 139269 : nout1 = ia - atn
2452 464676 : DO ib = 1, before
2453 325407 : nin1 = nin1 + after
2454 325407 : nin2 = nin1 + atb
2455 325407 : nin3 = nin2 + atb
2456 325407 : nin4 = nin3 + atb
2457 325407 : nout1 = nout1 + atn
2458 325407 : nout2 = nout1 + after
2459 325407 : nout3 = nout2 + after
2460 325407 : nout4 = nout3 + after
2461 3106539 : DO j = 1, nfft
2462 2641863 : r1 = zin(1, j, nin1)
2463 2641863 : s1 = zin(2, j, nin1)
2464 2641863 : r2 = zin(1, j, nin2)
2465 2641863 : s2 = zin(2, j, nin2)
2466 2641863 : r3 = zin(1, j, nin3)
2467 2641863 : s3 = zin(2, j, nin3)
2468 2641863 : r4 = zin(1, j, nin4)
2469 2641863 : s4 = zin(2, j, nin4)
2470 2641863 : r = r1 + r3
2471 2641863 : s = r2 + r4
2472 2641863 : zout(1, j, nout1) = r + s
2473 2641863 : zout(1, j, nout3) = r - s
2474 2641863 : r = r1 - r3
2475 2641863 : s = s2 - s4
2476 2641863 : zout(1, j, nout2) = r - s
2477 2641863 : zout(1, j, nout4) = r + s
2478 2641863 : r = s1 + s3
2479 2641863 : s = s2 + s4
2480 2641863 : zout(2, j, nout1) = r + s
2481 2641863 : zout(2, j, nout3) = r - s
2482 2641863 : r = s1 - s3
2483 2641863 : s = r2 - r4
2484 2641863 : zout(2, j, nout2) = r + s
2485 2967270 : zout(2, j, nout4) = r - s
2486 : END DO
2487 : END DO
2488 876753 : DO ia = 2, after
2489 737484 : ias = ia - 1
2490 876753 : IF (2*ias == after) THEN
2491 92904 : nin1 = ia - after
2492 92904 : nout1 = ia - atn
2493 279216 : DO ib = 1, before
2494 186312 : nin1 = nin1 + after
2495 186312 : nin2 = nin1 + atb
2496 186312 : nin3 = nin2 + atb
2497 186312 : nin4 = nin3 + atb
2498 186312 : nout1 = nout1 + atn
2499 186312 : nout2 = nout1 + after
2500 186312 : nout3 = nout2 + after
2501 186312 : nout4 = nout3 + after
2502 1950384 : DO j = 1, nfft
2503 1671168 : r1 = zin(1, j, nin1)
2504 1671168 : s1 = zin(2, j, nin1)
2505 1671168 : r = zin(1, j, nin2)
2506 1671168 : s = zin(2, j, nin2)
2507 1671168 : r2 = (r - s)*rt2i
2508 1671168 : s2 = (r + s)*rt2i
2509 1671168 : r3 = -zin(2, j, nin3)
2510 1671168 : s3 = zin(1, j, nin3)
2511 1671168 : r = zin(1, j, nin4)
2512 1671168 : s = zin(2, j, nin4)
2513 1671168 : r4 = -(r + s)*rt2i
2514 1671168 : s4 = (r - s)*rt2i
2515 1671168 : r = r1 + r3
2516 1671168 : s = r2 + r4
2517 1671168 : zout(1, j, nout1) = r + s
2518 1671168 : zout(1, j, nout3) = r - s
2519 1671168 : r = r1 - r3
2520 1671168 : s = s2 - s4
2521 1671168 : zout(1, j, nout2) = r - s
2522 1671168 : zout(1, j, nout4) = r + s
2523 1671168 : r = s1 + s3
2524 1671168 : s = s2 + s4
2525 1671168 : zout(2, j, nout1) = r + s
2526 1671168 : zout(2, j, nout3) = r - s
2527 1671168 : r = s1 - s3
2528 1671168 : s = r2 - r4
2529 1671168 : zout(2, j, nout2) = r + s
2530 1857480 : zout(2, j, nout4) = r - s
2531 : END DO
2532 : END DO
2533 : ELSE
2534 644580 : itt = ias*before
2535 644580 : itrig = itt + 1
2536 644580 : cr2 = trig(1, itrig)
2537 644580 : ci2 = trig(2, itrig)
2538 644580 : itrig = itrig + itt
2539 644580 : cr3 = trig(1, itrig)
2540 644580 : ci3 = trig(2, itrig)
2541 644580 : itrig = itrig + itt
2542 644580 : cr4 = trig(1, itrig)
2543 644580 : ci4 = trig(2, itrig)
2544 644580 : nin1 = ia - after
2545 644580 : nout1 = ia - atn
2546 2023920 : DO ib = 1, before
2547 1379340 : nin1 = nin1 + after
2548 1379340 : nin2 = nin1 + atb
2549 1379340 : nin3 = nin2 + atb
2550 1379340 : nin4 = nin3 + atb
2551 1379340 : nout1 = nout1 + atn
2552 1379340 : nout2 = nout1 + after
2553 1379340 : nout3 = nout2 + after
2554 1379340 : nout4 = nout3 + after
2555 13279500 : DO j = 1, nfft
2556 11255580 : r1 = zin(1, j, nin1)
2557 11255580 : s1 = zin(2, j, nin1)
2558 11255580 : r = zin(1, j, nin2)
2559 11255580 : s = zin(2, j, nin2)
2560 11255580 : r2 = r*cr2 - s*ci2
2561 11255580 : s2 = r*ci2 + s*cr2
2562 11255580 : r = zin(1, j, nin3)
2563 11255580 : s = zin(2, j, nin3)
2564 11255580 : r3 = r*cr3 - s*ci3
2565 11255580 : s3 = r*ci3 + s*cr3
2566 11255580 : r = zin(1, j, nin4)
2567 11255580 : s = zin(2, j, nin4)
2568 11255580 : r4 = r*cr4 - s*ci4
2569 11255580 : s4 = r*ci4 + s*cr4
2570 11255580 : r = r1 + r3
2571 11255580 : s = r2 + r4
2572 11255580 : zout(1, j, nout1) = r + s
2573 11255580 : zout(1, j, nout3) = r - s
2574 11255580 : r = r1 - r3
2575 11255580 : s = s2 - s4
2576 11255580 : zout(1, j, nout2) = r - s
2577 11255580 : zout(1, j, nout4) = r + s
2578 11255580 : r = s1 + s3
2579 11255580 : s = s2 + s4
2580 11255580 : zout(2, j, nout1) = r + s
2581 11255580 : zout(2, j, nout3) = r - s
2582 11255580 : r = s1 - s3
2583 11255580 : s = r2 - r4
2584 11255580 : zout(2, j, nout2) = r + s
2585 12634920 : zout(2, j, nout4) = r - s
2586 : END DO
2587 : END DO
2588 : END IF
2589 : END DO
2590 : ELSE
2591 63479 : ia = 1
2592 63479 : nin1 = ia - after
2593 63479 : nout1 = ia - atn
2594 253916 : DO ib = 1, before
2595 190437 : nin1 = nin1 + after
2596 190437 : nin2 = nin1 + atb
2597 190437 : nin3 = nin2 + atb
2598 190437 : nin4 = nin3 + atb
2599 190437 : nout1 = nout1 + atn
2600 190437 : nout2 = nout1 + after
2601 190437 : nout3 = nout2 + after
2602 190437 : nout4 = nout3 + after
2603 1715369 : DO j = 1, nfft
2604 1461453 : r1 = zin(1, j, nin1)
2605 1461453 : s1 = zin(2, j, nin1)
2606 1461453 : r2 = zin(1, j, nin2)
2607 1461453 : s2 = zin(2, j, nin2)
2608 1461453 : r3 = zin(1, j, nin3)
2609 1461453 : s3 = zin(2, j, nin3)
2610 1461453 : r4 = zin(1, j, nin4)
2611 1461453 : s4 = zin(2, j, nin4)
2612 1461453 : r = r1 + r3
2613 1461453 : s = r2 + r4
2614 1461453 : zout(1, j, nout1) = r + s
2615 1461453 : zout(1, j, nout3) = r - s
2616 1461453 : r = r1 - r3
2617 1461453 : s = s2 - s4
2618 1461453 : zout(1, j, nout2) = r + s
2619 1461453 : zout(1, j, nout4) = r - s
2620 1461453 : r = s1 + s3
2621 1461453 : s = s2 + s4
2622 1461453 : zout(2, j, nout1) = r + s
2623 1461453 : zout(2, j, nout3) = r - s
2624 1461453 : r = s1 - s3
2625 1461453 : s = r2 - r4
2626 1461453 : zout(2, j, nout2) = r - s
2627 1651890 : zout(2, j, nout4) = r + s
2628 : END DO
2629 : END DO
2630 353923 : DO ia = 2, after
2631 290444 : ias = ia - 1
2632 353923 : IF (2*ias == after) THEN
2633 46704 : nin1 = ia - after
2634 46704 : nout1 = ia - atn
2635 186816 : DO ib = 1, before
2636 140112 : nin1 = nin1 + after
2637 140112 : nin2 = nin1 + atb
2638 140112 : nin3 = nin2 + atb
2639 140112 : nin4 = nin3 + atb
2640 140112 : nout1 = nout1 + atn
2641 140112 : nout2 = nout1 + after
2642 140112 : nout3 = nout2 + after
2643 140112 : nout4 = nout3 + after
2644 1182144 : DO j = 1, nfft
2645 995328 : r1 = zin(1, j, nin1)
2646 995328 : s1 = zin(2, j, nin1)
2647 995328 : r = zin(1, j, nin2)
2648 995328 : s = zin(2, j, nin2)
2649 995328 : r2 = (r + s)*rt2i
2650 995328 : s2 = (s - r)*rt2i
2651 995328 : r3 = zin(2, j, nin3)
2652 995328 : s3 = -zin(1, j, nin3)
2653 995328 : r = zin(1, j, nin4)
2654 995328 : s = zin(2, j, nin4)
2655 995328 : r4 = (s - r)*rt2i
2656 995328 : s4 = -(r + s)*rt2i
2657 995328 : r = r1 + r3
2658 995328 : s = r2 + r4
2659 995328 : zout(1, j, nout1) = r + s
2660 995328 : zout(1, j, nout3) = r - s
2661 995328 : r = r1 - r3
2662 995328 : s = s2 - s4
2663 995328 : zout(1, j, nout2) = r + s
2664 995328 : zout(1, j, nout4) = r - s
2665 995328 : r = s1 + s3
2666 995328 : s = s2 + s4
2667 995328 : zout(2, j, nout1) = r + s
2668 995328 : zout(2, j, nout3) = r - s
2669 995328 : r = s1 - s3
2670 995328 : s = r2 - r4
2671 995328 : zout(2, j, nout2) = r - s
2672 1135440 : zout(2, j, nout4) = r + s
2673 : END DO
2674 : END DO
2675 : ELSE
2676 243740 : itt = ias*before
2677 243740 : itrig = itt + 1
2678 243740 : cr2 = trig(1, itrig)
2679 243740 : ci2 = trig(2, itrig)
2680 243740 : itrig = itrig + itt
2681 243740 : cr3 = trig(1, itrig)
2682 243740 : ci3 = trig(2, itrig)
2683 243740 : itrig = itrig + itt
2684 243740 : cr4 = trig(1, itrig)
2685 243740 : ci4 = trig(2, itrig)
2686 243740 : nin1 = ia - after
2687 243740 : nout1 = ia - atn
2688 974960 : DO ib = 1, before
2689 731220 : nin1 = nin1 + after
2690 731220 : nin2 = nin1 + atb
2691 731220 : nin3 = nin2 + atb
2692 731220 : nin4 = nin3 + atb
2693 731220 : nout1 = nout1 + atn
2694 731220 : nout2 = nout1 + after
2695 731220 : nout3 = nout2 + after
2696 731220 : nout4 = nout3 + after
2697 5586980 : DO j = 1, nfft
2698 4612020 : r1 = zin(1, j, nin1)
2699 4612020 : s1 = zin(2, j, nin1)
2700 4612020 : r = zin(1, j, nin2)
2701 4612020 : s = zin(2, j, nin2)
2702 4612020 : r2 = r*cr2 - s*ci2
2703 4612020 : s2 = r*ci2 + s*cr2
2704 4612020 : r = zin(1, j, nin3)
2705 4612020 : s = zin(2, j, nin3)
2706 4612020 : r3 = r*cr3 - s*ci3
2707 4612020 : s3 = r*ci3 + s*cr3
2708 4612020 : r = zin(1, j, nin4)
2709 4612020 : s = zin(2, j, nin4)
2710 4612020 : r4 = r*cr4 - s*ci4
2711 4612020 : s4 = r*ci4 + s*cr4
2712 4612020 : r = r1 + r3
2713 4612020 : s = r2 + r4
2714 4612020 : zout(1, j, nout1) = r + s
2715 4612020 : zout(1, j, nout3) = r - s
2716 4612020 : r = r1 - r3
2717 4612020 : s = s2 - s4
2718 4612020 : zout(1, j, nout2) = r + s
2719 4612020 : zout(1, j, nout4) = r - s
2720 4612020 : r = s1 + s3
2721 4612020 : s = s2 + s4
2722 4612020 : zout(2, j, nout1) = r + s
2723 4612020 : zout(2, j, nout3) = r - s
2724 4612020 : r = s1 - s3
2725 4612020 : s = r2 - r4
2726 4612020 : zout(2, j, nout2) = r - s
2727 5343240 : zout(2, j, nout4) = r + s
2728 : END DO
2729 : END DO
2730 : END IF
2731 : END DO
2732 : END IF
2733 1041398 : ELSE IF (now == 8) THEN
2734 49788 : IF (isign == -1) THEN
2735 49788 : ia = 1
2736 49788 : nin1 = ia - after
2737 49788 : nout1 = ia - atn
2738 249378 : DO ib = 1, before
2739 199590 : nin1 = nin1 + after
2740 199590 : nin2 = nin1 + atb
2741 199590 : nin3 = nin2 + atb
2742 199590 : nin4 = nin3 + atb
2743 199590 : nin5 = nin4 + atb
2744 199590 : nin6 = nin5 + atb
2745 199590 : nin7 = nin6 + atb
2746 199590 : nin8 = nin7 + atb
2747 199590 : nout1 = nout1 + atn
2748 199590 : nout2 = nout1 + after
2749 199590 : nout3 = nout2 + after
2750 199590 : nout4 = nout3 + after
2751 199590 : nout5 = nout4 + after
2752 199590 : nout6 = nout5 + after
2753 199590 : nout7 = nout6 + after
2754 199590 : nout8 = nout7 + after
2755 3161058 : DO j = 1, nfft
2756 2911680 : r1 = zin(1, j, nin1)
2757 2911680 : s1 = zin(2, j, nin1)
2758 2911680 : r2 = zin(1, j, nin2)
2759 2911680 : s2 = zin(2, j, nin2)
2760 2911680 : r3 = zin(1, j, nin3)
2761 2911680 : s3 = zin(2, j, nin3)
2762 2911680 : r4 = zin(1, j, nin4)
2763 2911680 : s4 = zin(2, j, nin4)
2764 2911680 : r5 = zin(1, j, nin5)
2765 2911680 : s5 = zin(2, j, nin5)
2766 2911680 : r6 = zin(1, j, nin6)
2767 2911680 : s6 = zin(2, j, nin6)
2768 2911680 : r7 = zin(1, j, nin7)
2769 2911680 : s7 = zin(2, j, nin7)
2770 2911680 : r8 = zin(1, j, nin8)
2771 2911680 : s8 = zin(2, j, nin8)
2772 2911680 : r = r1 + r5
2773 2911680 : s = r3 + r7
2774 2911680 : ap = r + s
2775 2911680 : am = r - s
2776 2911680 : r = r2 + r6
2777 2911680 : s = r4 + r8
2778 2911680 : bp = r + s
2779 2911680 : bm = r - s
2780 2911680 : r = s1 + s5
2781 2911680 : s = s3 + s7
2782 2911680 : cp = r + s
2783 2911680 : cm = r - s
2784 2911680 : r = s2 + s6
2785 2911680 : s = s4 + s8
2786 2911680 : dbl = r + s
2787 2911680 : dm = r - s
2788 2911680 : zout(1, j, nout1) = ap + bp
2789 2911680 : zout(2, j, nout1) = cp + dbl
2790 2911680 : zout(1, j, nout5) = ap - bp
2791 2911680 : zout(2, j, nout5) = cp - dbl
2792 2911680 : zout(1, j, nout3) = am + dm
2793 2911680 : zout(2, j, nout3) = cm - bm
2794 2911680 : zout(1, j, nout7) = am - dm
2795 2911680 : zout(2, j, nout7) = cm + bm
2796 2911680 : r = r1 - r5
2797 2911680 : s = s3 - s7
2798 2911680 : ap = r + s
2799 2911680 : am = r - s
2800 2911680 : r = s1 - s5
2801 2911680 : s = r3 - r7
2802 2911680 : bp = r + s
2803 2911680 : bm = r - s
2804 2911680 : r = s4 - s8
2805 2911680 : s = r2 - r6
2806 2911680 : cp = r + s
2807 2911680 : cm = r - s
2808 2911680 : r = s2 - s6
2809 2911680 : s = r4 - r8
2810 2911680 : dbl = r + s
2811 2911680 : dm = r - s
2812 2911680 : r = (cp + dm)*rt2i
2813 2911680 : s = (-cp + dm)*rt2i
2814 2911680 : cp = (cm + dbl)*rt2i
2815 2911680 : dbl = (cm - dbl)*rt2i
2816 2911680 : zout(1, j, nout2) = ap + r
2817 2911680 : zout(2, j, nout2) = bm + s
2818 2911680 : zout(1, j, nout6) = ap - r
2819 2911680 : zout(2, j, nout6) = bm - s
2820 2911680 : zout(1, j, nout4) = am + cp
2821 2911680 : zout(2, j, nout4) = bp + dbl
2822 2911680 : zout(1, j, nout8) = am - cp
2823 3111270 : zout(2, j, nout8) = bp - dbl
2824 : END DO
2825 : END DO
2826 : ELSE
2827 0 : ia = 1
2828 0 : nin1 = ia - after
2829 0 : nout1 = ia - atn
2830 0 : DO ib = 1, before
2831 0 : nin1 = nin1 + after
2832 0 : nin2 = nin1 + atb
2833 0 : nin3 = nin2 + atb
2834 0 : nin4 = nin3 + atb
2835 0 : nin5 = nin4 + atb
2836 0 : nin6 = nin5 + atb
2837 0 : nin7 = nin6 + atb
2838 0 : nin8 = nin7 + atb
2839 0 : nout1 = nout1 + atn
2840 0 : nout2 = nout1 + after
2841 0 : nout3 = nout2 + after
2842 0 : nout4 = nout3 + after
2843 0 : nout5 = nout4 + after
2844 0 : nout6 = nout5 + after
2845 0 : nout7 = nout6 + after
2846 0 : nout8 = nout7 + after
2847 0 : DO j = 1, nfft
2848 0 : r1 = zin(1, j, nin1)
2849 0 : s1 = zin(2, j, nin1)
2850 0 : r2 = zin(1, j, nin2)
2851 0 : s2 = zin(2, j, nin2)
2852 0 : r3 = zin(1, j, nin3)
2853 0 : s3 = zin(2, j, nin3)
2854 0 : r4 = zin(1, j, nin4)
2855 0 : s4 = zin(2, j, nin4)
2856 0 : r5 = zin(1, j, nin5)
2857 0 : s5 = zin(2, j, nin5)
2858 0 : r6 = zin(1, j, nin6)
2859 0 : s6 = zin(2, j, nin6)
2860 0 : r7 = zin(1, j, nin7)
2861 0 : s7 = zin(2, j, nin7)
2862 0 : r8 = zin(1, j, nin8)
2863 0 : s8 = zin(2, j, nin8)
2864 0 : r = r1 + r5
2865 0 : s = r3 + r7
2866 0 : ap = r + s
2867 0 : am = r - s
2868 0 : r = r2 + r6
2869 0 : s = r4 + r8
2870 0 : bp = r + s
2871 0 : bm = r - s
2872 0 : r = s1 + s5
2873 0 : s = s3 + s7
2874 0 : cp = r + s
2875 0 : cm = r - s
2876 0 : r = s2 + s6
2877 0 : s = s4 + s8
2878 0 : dbl = r + s
2879 0 : dm = r - s
2880 0 : zout(1, j, nout1) = ap + bp
2881 0 : zout(2, j, nout1) = cp + dbl
2882 0 : zout(1, j, nout5) = ap - bp
2883 0 : zout(2, j, nout5) = cp - dbl
2884 0 : zout(1, j, nout3) = am - dm
2885 0 : zout(2, j, nout3) = cm + bm
2886 0 : zout(1, j, nout7) = am + dm
2887 0 : zout(2, j, nout7) = cm - bm
2888 0 : r = r1 - r5
2889 0 : s = -s3 + s7
2890 0 : ap = r + s
2891 0 : am = r - s
2892 0 : r = s1 - s5
2893 0 : s = r7 - r3
2894 0 : bp = r + s
2895 0 : bm = r - s
2896 0 : r = -s4 + s8
2897 0 : s = r2 - r6
2898 0 : cp = r + s
2899 0 : cm = r - s
2900 0 : r = -s2 + s6
2901 0 : s = r4 - r8
2902 0 : dbl = r + s
2903 0 : dm = r - s
2904 0 : r = (cp + dm)*rt2i
2905 0 : s = (cp - dm)*rt2i
2906 0 : cp = (cm + dbl)*rt2i
2907 0 : dbl = (-cm + dbl)*rt2i
2908 0 : zout(1, j, nout2) = ap + r
2909 0 : zout(2, j, nout2) = bm + s
2910 0 : zout(1, j, nout6) = ap - r
2911 0 : zout(2, j, nout6) = bm - s
2912 0 : zout(1, j, nout4) = am + cp
2913 0 : zout(2, j, nout4) = bp + dbl
2914 0 : zout(1, j, nout8) = am - cp
2915 0 : zout(2, j, nout8) = bp - dbl
2916 : END DO
2917 : END DO
2918 : END IF
2919 991610 : ELSE IF (now == 3) THEN
2920 368319 : bbs = isign*bb
2921 368319 : ia = 1
2922 368319 : nin1 = ia - after
2923 368319 : nout1 = ia - atn
2924 825270 : DO ib = 1, before
2925 456951 : nin1 = nin1 + after
2926 456951 : nin2 = nin1 + atb
2927 456951 : nin3 = nin2 + atb
2928 456951 : nout1 = nout1 + atn
2929 456951 : nout2 = nout1 + after
2930 456951 : nout3 = nout2 + after
2931 4766403 : DO j = 1, nfft
2932 3941133 : r1 = zin(1, j, nin1)
2933 3941133 : s1 = zin(2, j, nin1)
2934 3941133 : r2 = zin(1, j, nin2)
2935 3941133 : s2 = zin(2, j, nin2)
2936 3941133 : r3 = zin(1, j, nin3)
2937 3941133 : s3 = zin(2, j, nin3)
2938 3941133 : r = r2 + r3
2939 3941133 : s = s2 + s3
2940 3941133 : zout(1, j, nout1) = r + r1
2941 3941133 : zout(2, j, nout1) = s + s1
2942 3941133 : r1 = r1 - 0.5_dp*r
2943 3941133 : s1 = s1 - 0.5_dp*s
2944 3941133 : r2 = bbs*(r2 - r3)
2945 3941133 : s2 = bbs*(s2 - s3)
2946 3941133 : zout(1, j, nout2) = r1 - s2
2947 3941133 : zout(2, j, nout2) = s1 + r2
2948 3941133 : zout(1, j, nout3) = r1 + s2
2949 4398084 : zout(2, j, nout3) = s1 - r2
2950 : END DO
2951 : END DO
2952 8377314 : DO ia = 2, after
2953 8008995 : ias = ia - 1
2954 8377314 : IF (4*ias == 3*after) THEN
2955 151029 : IF (isign == 1) THEN
2956 99357 : nin1 = ia - after
2957 99357 : nout1 = ia - atn
2958 198714 : DO ib = 1, before
2959 99357 : nin1 = nin1 + after
2960 99357 : nin2 = nin1 + atb
2961 99357 : nin3 = nin2 + atb
2962 99357 : nout1 = nout1 + atn
2963 99357 : nout2 = nout1 + after
2964 99357 : nout3 = nout2 + after
2965 1077831 : DO j = 1, nfft
2966 879117 : r1 = zin(1, j, nin1)
2967 879117 : s1 = zin(2, j, nin1)
2968 879117 : r2 = -zin(2, j, nin2)
2969 879117 : s2 = zin(1, j, nin2)
2970 879117 : r3 = -zin(1, j, nin3)
2971 879117 : s3 = -zin(2, j, nin3)
2972 879117 : r = r2 + r3
2973 879117 : s = s2 + s3
2974 879117 : zout(1, j, nout1) = r + r1
2975 879117 : zout(2, j, nout1) = s + s1
2976 879117 : r1 = r1 - 0.5_dp*r
2977 879117 : s1 = s1 - 0.5_dp*s
2978 879117 : r2 = bbs*(r2 - r3)
2979 879117 : s2 = bbs*(s2 - s3)
2980 879117 : zout(1, j, nout2) = r1 - s2
2981 879117 : zout(2, j, nout2) = s1 + r2
2982 879117 : zout(1, j, nout3) = r1 + s2
2983 978474 : zout(2, j, nout3) = s1 - r2
2984 : END DO
2985 : END DO
2986 : ELSE
2987 51672 : nin1 = ia - after
2988 51672 : nout1 = ia - atn
2989 103344 : DO ib = 1, before
2990 51672 : nin1 = nin1 + after
2991 51672 : nin2 = nin1 + atb
2992 51672 : nin3 = nin2 + atb
2993 51672 : nout1 = nout1 + atn
2994 51672 : nout2 = nout1 + after
2995 51672 : nout3 = nout2 + after
2996 611376 : DO j = 1, nfft
2997 508032 : r1 = zin(1, j, nin1)
2998 508032 : s1 = zin(2, j, nin1)
2999 508032 : r2 = zin(2, j, nin2)
3000 508032 : s2 = -zin(1, j, nin2)
3001 508032 : r3 = -zin(1, j, nin3)
3002 508032 : s3 = -zin(2, j, nin3)
3003 508032 : r = r2 + r3
3004 508032 : s = s2 + s3
3005 508032 : zout(1, j, nout1) = r + r1
3006 508032 : zout(2, j, nout1) = s + s1
3007 508032 : r1 = r1 - 0.5_dp*r
3008 508032 : s1 = s1 - 0.5_dp*s
3009 508032 : r2 = bbs*(r2 - r3)
3010 508032 : s2 = bbs*(s2 - s3)
3011 508032 : zout(1, j, nout2) = r1 - s2
3012 508032 : zout(2, j, nout2) = s1 + r2
3013 508032 : zout(1, j, nout3) = r1 + s2
3014 559704 : zout(2, j, nout3) = s1 - r2
3015 : END DO
3016 : END DO
3017 : END IF
3018 7857966 : ELSE IF (8*ias == 3*after) THEN
3019 93744 : IF (isign == 1) THEN
3020 46872 : nin1 = ia - after
3021 46872 : nout1 = ia - atn
3022 93744 : DO ib = 1, before
3023 46872 : nin1 = nin1 + after
3024 46872 : nin2 = nin1 + atb
3025 46872 : nin3 = nin2 + atb
3026 46872 : nout1 = nout1 + atn
3027 46872 : nout2 = nout1 + after
3028 46872 : nout3 = nout2 + after
3029 428976 : DO j = 1, nfft
3030 335232 : r1 = zin(1, j, nin1)
3031 335232 : s1 = zin(2, j, nin1)
3032 335232 : r = zin(1, j, nin2)
3033 335232 : s = zin(2, j, nin2)
3034 335232 : r2 = (r - s)*rt2i
3035 335232 : s2 = (r + s)*rt2i
3036 335232 : r3 = -zin(2, j, nin3)
3037 335232 : s3 = zin(1, j, nin3)
3038 335232 : r = r2 + r3
3039 335232 : s = s2 + s3
3040 335232 : zout(1, j, nout1) = r + r1
3041 335232 : zout(2, j, nout1) = s + s1
3042 335232 : r1 = r1 - 0.5_dp*r
3043 335232 : s1 = s1 - 0.5_dp*s
3044 335232 : r2 = bbs*(r2 - r3)
3045 335232 : s2 = bbs*(s2 - s3)
3046 335232 : zout(1, j, nout2) = r1 - s2
3047 335232 : zout(2, j, nout2) = s1 + r2
3048 335232 : zout(1, j, nout3) = r1 + s2
3049 382104 : zout(2, j, nout3) = s1 - r2
3050 : END DO
3051 : END DO
3052 : ELSE
3053 46872 : nin1 = ia - after
3054 46872 : nout1 = ia - atn
3055 93744 : DO ib = 1, before
3056 46872 : nin1 = nin1 + after
3057 46872 : nin2 = nin1 + atb
3058 46872 : nin3 = nin2 + atb
3059 46872 : nout1 = nout1 + atn
3060 46872 : nout2 = nout1 + after
3061 46872 : nout3 = nout2 + after
3062 428976 : DO j = 1, nfft
3063 335232 : r1 = zin(1, j, nin1)
3064 335232 : s1 = zin(2, j, nin1)
3065 335232 : r = zin(1, j, nin2)
3066 335232 : s = zin(2, j, nin2)
3067 335232 : r2 = (r + s)*rt2i
3068 335232 : s2 = (-r + s)*rt2i
3069 335232 : r3 = zin(2, j, nin3)
3070 335232 : s3 = -zin(1, j, nin3)
3071 335232 : r = r2 + r3
3072 335232 : s = s2 + s3
3073 335232 : zout(1, j, nout1) = r + r1
3074 335232 : zout(2, j, nout1) = s + s1
3075 335232 : r1 = r1 - 0.5_dp*r
3076 335232 : s1 = s1 - 0.5_dp*s
3077 335232 : r2 = bbs*(r2 - r3)
3078 335232 : s2 = bbs*(s2 - s3)
3079 335232 : zout(1, j, nout2) = r1 - s2
3080 335232 : zout(2, j, nout2) = s1 + r2
3081 335232 : zout(1, j, nout3) = r1 + s2
3082 382104 : zout(2, j, nout3) = s1 - r2
3083 : END DO
3084 : END DO
3085 : END IF
3086 : ELSE
3087 7764222 : itt = ias*before
3088 7764222 : itrig = itt + 1
3089 7764222 : cr2 = trig(1, itrig)
3090 7764222 : ci2 = trig(2, itrig)
3091 7764222 : itrig = itrig + itt
3092 7764222 : cr3 = trig(1, itrig)
3093 7764222 : ci3 = trig(2, itrig)
3094 7764222 : nin1 = ia - after
3095 7764222 : nout1 = ia - atn
3096 15856380 : DO ib = 1, before
3097 8092158 : nin1 = nin1 + after
3098 8092158 : nin2 = nin1 + atb
3099 8092158 : nin3 = nin2 + atb
3100 8092158 : nout1 = nout1 + atn
3101 8092158 : nout2 = nout1 + after
3102 8092158 : nout3 = nout2 + after
3103 65276424 : DO j = 1, nfft
3104 49420044 : r1 = zin(1, j, nin1)
3105 49420044 : s1 = zin(2, j, nin1)
3106 49420044 : r = zin(1, j, nin2)
3107 49420044 : s = zin(2, j, nin2)
3108 49420044 : r2 = r*cr2 - s*ci2
3109 49420044 : s2 = r*ci2 + s*cr2
3110 49420044 : r = zin(1, j, nin3)
3111 49420044 : s = zin(2, j, nin3)
3112 49420044 : r3 = r*cr3 - s*ci3
3113 49420044 : s3 = r*ci3 + s*cr3
3114 49420044 : r = r2 + r3
3115 49420044 : s = s2 + s3
3116 49420044 : zout(1, j, nout1) = r + r1
3117 49420044 : zout(2, j, nout1) = s + s1
3118 49420044 : r1 = r1 - 0.5_dp*r
3119 49420044 : s1 = s1 - 0.5_dp*s
3120 49420044 : r2 = bbs*(r2 - r3)
3121 49420044 : s2 = bbs*(s2 - s3)
3122 49420044 : zout(1, j, nout2) = r1 - s2
3123 49420044 : zout(2, j, nout2) = s1 + r2
3124 49420044 : zout(1, j, nout3) = r1 + s2
3125 57512202 : zout(2, j, nout3) = s1 - r2
3126 : END DO
3127 : END DO
3128 : END IF
3129 : END DO
3130 623291 : ELSE IF (now == 5) THEN
3131 427043 : sin2 = isign*sin2p
3132 427043 : sin4 = isign*sin4p
3133 427043 : ia = 1
3134 427043 : nin1 = ia - after
3135 427043 : nout1 = ia - atn
3136 1816259 : DO ib = 1, before
3137 1389216 : nin1 = nin1 + after
3138 1389216 : nin2 = nin1 + atb
3139 1389216 : nin3 = nin2 + atb
3140 1389216 : nin4 = nin3 + atb
3141 1389216 : nin5 = nin4 + atb
3142 1389216 : nout1 = nout1 + atn
3143 1389216 : nout2 = nout1 + after
3144 1389216 : nout3 = nout2 + after
3145 1389216 : nout4 = nout3 + after
3146 1389216 : nout5 = nout4 + after
3147 12033719 : DO j = 1, nfft
3148 10217460 : r1 = zin(1, j, nin1)
3149 10217460 : s1 = zin(2, j, nin1)
3150 10217460 : r2 = zin(1, j, nin2)
3151 10217460 : s2 = zin(2, j, nin2)
3152 10217460 : r3 = zin(1, j, nin3)
3153 10217460 : s3 = zin(2, j, nin3)
3154 10217460 : r4 = zin(1, j, nin4)
3155 10217460 : s4 = zin(2, j, nin4)
3156 10217460 : r5 = zin(1, j, nin5)
3157 10217460 : s5 = zin(2, j, nin5)
3158 10217460 : r25 = r2 + r5
3159 10217460 : r34 = r3 + r4
3160 10217460 : s25 = s2 - s5
3161 10217460 : s34 = s3 - s4
3162 10217460 : zout(1, j, nout1) = r1 + r25 + r34
3163 10217460 : r = cos2*r25 + cos4*r34 + r1
3164 10217460 : s = sin2*s25 + sin4*s34
3165 10217460 : zout(1, j, nout2) = r - s
3166 10217460 : zout(1, j, nout5) = r + s
3167 10217460 : r = cos4*r25 + cos2*r34 + r1
3168 10217460 : s = sin4*s25 - sin2*s34
3169 10217460 : zout(1, j, nout3) = r - s
3170 10217460 : zout(1, j, nout4) = r + s
3171 10217460 : r25 = r2 - r5
3172 10217460 : r34 = r3 - r4
3173 10217460 : s25 = s2 + s5
3174 10217460 : s34 = s3 + s4
3175 10217460 : zout(2, j, nout1) = s1 + s25 + s34
3176 10217460 : r = cos2*s25 + cos4*s34 + s1
3177 10217460 : s = sin2*r25 + sin4*r34
3178 10217460 : zout(2, j, nout2) = r + s
3179 10217460 : zout(2, j, nout5) = r - s
3180 10217460 : r = cos4*s25 + cos2*s34 + s1
3181 10217460 : s = sin4*r25 - sin2*r34
3182 10217460 : zout(2, j, nout3) = r + s
3183 11606676 : zout(2, j, nout4) = r - s
3184 : END DO
3185 : END DO
3186 2268347 : DO ia = 2, after
3187 1841304 : ias = ia - 1
3188 2268347 : IF (8*ias == 5*after) THEN
3189 0 : IF (isign == 1) THEN
3190 0 : nin1 = ia - after
3191 0 : nout1 = ia - atn
3192 0 : DO ib = 1, before
3193 0 : nin1 = nin1 + after
3194 0 : nin2 = nin1 + atb
3195 0 : nin3 = nin2 + atb
3196 0 : nin4 = nin3 + atb
3197 0 : nin5 = nin4 + atb
3198 0 : nout1 = nout1 + atn
3199 0 : nout2 = nout1 + after
3200 0 : nout3 = nout2 + after
3201 0 : nout4 = nout3 + after
3202 0 : nout5 = nout4 + after
3203 0 : DO j = 1, nfft
3204 0 : r1 = zin(1, j, nin1)
3205 0 : s1 = zin(2, j, nin1)
3206 0 : r = zin(1, j, nin2)
3207 0 : s = zin(2, j, nin2)
3208 0 : r2 = (r - s)*rt2i
3209 0 : s2 = (r + s)*rt2i
3210 0 : r3 = -zin(2, j, nin3)
3211 0 : s3 = zin(1, j, nin3)
3212 0 : r = zin(1, j, nin4)
3213 0 : s = zin(2, j, nin4)
3214 0 : r4 = -(r + s)*rt2i
3215 0 : s4 = (r - s)*rt2i
3216 0 : r5 = -zin(1, j, nin5)
3217 0 : s5 = -zin(2, j, nin5)
3218 0 : r25 = r2 + r5
3219 0 : r34 = r3 + r4
3220 0 : s25 = s2 - s5
3221 0 : s34 = s3 - s4
3222 0 : zout(1, j, nout1) = r1 + r25 + r34
3223 0 : r = cos2*r25 + cos4*r34 + r1
3224 0 : s = sin2*s25 + sin4*s34
3225 0 : zout(1, j, nout2) = r - s
3226 0 : zout(1, j, nout5) = r + s
3227 0 : r = cos4*r25 + cos2*r34 + r1
3228 0 : s = sin4*s25 - sin2*s34
3229 0 : zout(1, j, nout3) = r - s
3230 0 : zout(1, j, nout4) = r + s
3231 0 : r25 = r2 - r5
3232 0 : r34 = r3 - r4
3233 0 : s25 = s2 + s5
3234 0 : s34 = s3 + s4
3235 0 : zout(2, j, nout1) = s1 + s25 + s34
3236 0 : r = cos2*s25 + cos4*s34 + s1
3237 0 : s = sin2*r25 + sin4*r34
3238 0 : zout(2, j, nout2) = r + s
3239 0 : zout(2, j, nout5) = r - s
3240 0 : r = cos4*s25 + cos2*s34 + s1
3241 0 : s = sin4*r25 - sin2*r34
3242 0 : zout(2, j, nout3) = r + s
3243 0 : zout(2, j, nout4) = r - s
3244 : END DO
3245 : END DO
3246 : ELSE
3247 0 : nin1 = ia - after
3248 0 : nout1 = ia - atn
3249 0 : DO ib = 1, before
3250 0 : nin1 = nin1 + after
3251 0 : nin2 = nin1 + atb
3252 0 : nin3 = nin2 + atb
3253 0 : nin4 = nin3 + atb
3254 0 : nin5 = nin4 + atb
3255 0 : nout1 = nout1 + atn
3256 0 : nout2 = nout1 + after
3257 0 : nout3 = nout2 + after
3258 0 : nout4 = nout3 + after
3259 0 : nout5 = nout4 + after
3260 0 : DO j = 1, nfft
3261 0 : r1 = zin(1, j, nin1)
3262 0 : s1 = zin(2, j, nin1)
3263 0 : r = zin(1, j, nin2)
3264 0 : s = zin(2, j, nin2)
3265 0 : r2 = (r + s)*rt2i
3266 0 : s2 = (-r + s)*rt2i
3267 0 : r3 = zin(2, j, nin3)
3268 0 : s3 = -zin(1, j, nin3)
3269 0 : r = zin(1, j, nin4)
3270 0 : s = zin(2, j, nin4)
3271 0 : r4 = (s - r)*rt2i
3272 0 : s4 = -(r + s)*rt2i
3273 0 : r5 = -zin(1, j, nin5)
3274 0 : s5 = -zin(2, j, nin5)
3275 0 : r25 = r2 + r5
3276 0 : r34 = r3 + r4
3277 0 : s25 = s2 - s5
3278 0 : s34 = s3 - s4
3279 0 : zout(1, j, nout1) = r1 + r25 + r34
3280 0 : r = cos2*r25 + cos4*r34 + r1
3281 0 : s = sin2*s25 + sin4*s34
3282 0 : zout(1, j, nout2) = r - s
3283 0 : zout(1, j, nout5) = r + s
3284 0 : r = cos4*r25 + cos2*r34 + r1
3285 0 : s = sin4*s25 - sin2*s34
3286 0 : zout(1, j, nout3) = r - s
3287 0 : zout(1, j, nout4) = r + s
3288 0 : r25 = r2 - r5
3289 0 : r34 = r3 - r4
3290 0 : s25 = s2 + s5
3291 0 : s34 = s3 + s4
3292 0 : zout(2, j, nout1) = s1 + s25 + s34
3293 0 : r = cos2*s25 + cos4*s34 + s1
3294 0 : s = sin2*r25 + sin4*r34
3295 0 : zout(2, j, nout2) = r + s
3296 0 : zout(2, j, nout5) = r - s
3297 0 : r = cos4*s25 + cos2*s34 + s1
3298 0 : s = sin4*r25 - sin2*r34
3299 0 : zout(2, j, nout3) = r + s
3300 0 : zout(2, j, nout4) = r - s
3301 : END DO
3302 : END DO
3303 : END IF
3304 : ELSE
3305 1841304 : ias = ia - 1
3306 1841304 : itt = ias*before
3307 1841304 : itrig = itt + 1
3308 1841304 : cr2 = trig(1, itrig)
3309 1841304 : ci2 = trig(2, itrig)
3310 1841304 : itrig = itrig + itt
3311 1841304 : cr3 = trig(1, itrig)
3312 1841304 : ci3 = trig(2, itrig)
3313 1841304 : itrig = itrig + itt
3314 1841304 : cr4 = trig(1, itrig)
3315 1841304 : ci4 = trig(2, itrig)
3316 1841304 : itrig = itrig + itt
3317 1841304 : cr5 = trig(1, itrig)
3318 1841304 : ci5 = trig(2, itrig)
3319 1841304 : nin1 = ia - after
3320 1841304 : nout1 = ia - atn
3321 6987408 : DO ib = 1, before
3322 5146104 : nin1 = nin1 + after
3323 5146104 : nin2 = nin1 + atb
3324 5146104 : nin3 = nin2 + atb
3325 5146104 : nin4 = nin3 + atb
3326 5146104 : nin5 = nin4 + atb
3327 5146104 : nout1 = nout1 + atn
3328 5146104 : nout2 = nout1 + after
3329 5146104 : nout3 = nout2 + after
3330 5146104 : nout4 = nout3 + after
3331 5146104 : nout5 = nout4 + after
3332 35245908 : DO j = 1, nfft
3333 28258500 : r1 = zin(1, j, nin1)
3334 28258500 : s1 = zin(2, j, nin1)
3335 28258500 : r = zin(1, j, nin2)
3336 28258500 : s = zin(2, j, nin2)
3337 28258500 : r2 = r*cr2 - s*ci2
3338 28258500 : s2 = r*ci2 + s*cr2
3339 28258500 : r = zin(1, j, nin3)
3340 28258500 : s = zin(2, j, nin3)
3341 28258500 : r3 = r*cr3 - s*ci3
3342 28258500 : s3 = r*ci3 + s*cr3
3343 28258500 : r = zin(1, j, nin4)
3344 28258500 : s = zin(2, j, nin4)
3345 28258500 : r4 = r*cr4 - s*ci4
3346 28258500 : s4 = r*ci4 + s*cr4
3347 28258500 : r = zin(1, j, nin5)
3348 28258500 : s = zin(2, j, nin5)
3349 28258500 : r5 = r*cr5 - s*ci5
3350 28258500 : s5 = r*ci5 + s*cr5
3351 28258500 : r25 = r2 + r5
3352 28258500 : r34 = r3 + r4
3353 28258500 : s25 = s2 - s5
3354 28258500 : s34 = s3 - s4
3355 28258500 : zout(1, j, nout1) = r1 + r25 + r34
3356 28258500 : r = cos2*r25 + cos4*r34 + r1
3357 28258500 : s = sin2*s25 + sin4*s34
3358 28258500 : zout(1, j, nout2) = r - s
3359 28258500 : zout(1, j, nout5) = r + s
3360 28258500 : r = cos4*r25 + cos2*r34 + r1
3361 28258500 : s = sin4*s25 - sin2*s34
3362 28258500 : zout(1, j, nout3) = r - s
3363 28258500 : zout(1, j, nout4) = r + s
3364 28258500 : r25 = r2 - r5
3365 28258500 : r34 = r3 - r4
3366 28258500 : s25 = s2 + s5
3367 28258500 : s34 = s3 + s4
3368 28258500 : zout(2, j, nout1) = s1 + s25 + s34
3369 28258500 : r = cos2*s25 + cos4*s34 + s1
3370 28258500 : s = sin2*r25 + sin4*r34
3371 28258500 : zout(2, j, nout2) = r + s
3372 28258500 : zout(2, j, nout5) = r - s
3373 28258500 : r = cos4*s25 + cos2*s34 + s1
3374 28258500 : s = sin4*r25 - sin2*r34
3375 28258500 : zout(2, j, nout3) = r + s
3376 33404604 : zout(2, j, nout4) = r - s
3377 : END DO
3378 : END DO
3379 : END IF
3380 : END DO
3381 196248 : ELSE IF (now == 6) THEN
3382 196248 : bbs = isign*bb
3383 196248 : ia = 1
3384 196248 : nin1 = ia - after
3385 196248 : nout1 = ia - atn
3386 2975424 : DO ib = 1, before
3387 2779176 : nin1 = nin1 + after
3388 2779176 : nin2 = nin1 + atb
3389 2779176 : nin3 = nin2 + atb
3390 2779176 : nin4 = nin3 + atb
3391 2779176 : nin5 = nin4 + atb
3392 2779176 : nin6 = nin5 + atb
3393 2779176 : nout1 = nout1 + atn
3394 2779176 : nout2 = nout1 + after
3395 2779176 : nout3 = nout2 + after
3396 2779176 : nout4 = nout3 + after
3397 2779176 : nout5 = nout4 + after
3398 2779176 : nout6 = nout5 + after
3399 17746152 : DO j = 1, nfft
3400 14770728 : r2 = zin(1, j, nin3)
3401 14770728 : s2 = zin(2, j, nin3)
3402 14770728 : r3 = zin(1, j, nin5)
3403 14770728 : s3 = zin(2, j, nin5)
3404 14770728 : r = r2 + r3
3405 14770728 : s = s2 + s3
3406 14770728 : r1 = zin(1, j, nin1)
3407 14770728 : s1 = zin(2, j, nin1)
3408 14770728 : ur1 = r + r1
3409 14770728 : ui1 = s + s1
3410 14770728 : r1 = r1 - 0.5_dp*r
3411 14770728 : s1 = s1 - 0.5_dp*s
3412 14770728 : r = r2 - r3
3413 14770728 : s = s2 - s3
3414 14770728 : ur2 = r1 - s*bbs
3415 14770728 : ui2 = s1 + r*bbs
3416 14770728 : ur3 = r1 + s*bbs
3417 14770728 : ui3 = s1 - r*bbs
3418 :
3419 14770728 : r2 = zin(1, j, nin6)
3420 14770728 : s2 = zin(2, j, nin6)
3421 14770728 : r3 = zin(1, j, nin2)
3422 14770728 : s3 = zin(2, j, nin2)
3423 14770728 : r = r2 + r3
3424 14770728 : s = s2 + s3
3425 14770728 : r1 = zin(1, j, nin4)
3426 14770728 : s1 = zin(2, j, nin4)
3427 14770728 : vr1 = r + r1
3428 14770728 : vi1 = s + s1
3429 14770728 : r1 = r1 - 0.5_dp*r
3430 14770728 : s1 = s1 - 0.5_dp*s
3431 14770728 : r = r2 - r3
3432 14770728 : s = s2 - s3
3433 14770728 : vr2 = r1 - s*bbs
3434 14770728 : vi2 = s1 + r*bbs
3435 14770728 : vr3 = r1 + s*bbs
3436 14770728 : vi3 = s1 - r*bbs
3437 :
3438 14770728 : zout(1, j, nout1) = ur1 + vr1
3439 14770728 : zout(2, j, nout1) = ui1 + vi1
3440 14770728 : zout(1, j, nout5) = ur2 + vr2
3441 14770728 : zout(2, j, nout5) = ui2 + vi2
3442 14770728 : zout(1, j, nout3) = ur3 + vr3
3443 14770728 : zout(2, j, nout3) = ui3 + vi3
3444 14770728 : zout(1, j, nout4) = ur1 - vr1
3445 14770728 : zout(2, j, nout4) = ui1 - vi1
3446 14770728 : zout(1, j, nout2) = ur2 - vr2
3447 14770728 : zout(2, j, nout2) = ui2 - vi2
3448 14770728 : zout(1, j, nout6) = ur3 - vr3
3449 17549904 : zout(2, j, nout6) = ui3 - vi3
3450 : END DO
3451 : END DO
3452 : ELSE
3453 0 : CPABORT('Error fftstp')
3454 : END IF
3455 :
3456 : !-----------------------------------------------------------------------------!
3457 :
3458 1244146 : END SUBROUTINE fftstp
3459 :
3460 : !-----------------------------------------------------------------------------!
3461 : !-----------------------------------------------------------------------------!
3462 : ! Copyright (C) Stefan Goedecker, Lausanne, Switzerland, August 1, 1991
3463 : ! Copyright (C) Stefan Goedecker, Cornell University, Ithaca, USA, 1994
3464 : ! Copyright (C) Stefan Goedecker, MPI Stuttgart, Germany, 1999
3465 : ! This file is distributed under the terms of the
3466 : ! GNU General Public License version 2 (or later),
3467 : ! see http://www.gnu.org/copyleft/gpl.txt .
3468 : !-----------------------------------------------------------------------------!
3469 : ! S. Goedecker: Rotating a three-dimensional array in optimal
3470 : ! positions for vector processing: Case study for a three-dimensional Fast
3471 : ! Fourier Transform, Comp. Phys. Commun. 76, 294 (1993)
3472 : ! **************************************************************************************************
3473 : !> \brief ...
3474 : !> \param n ...
3475 : !> \param trig ...
3476 : !> \param after ...
3477 : !> \param before ...
3478 : !> \param now ...
3479 : !> \param isign ...
3480 : !> \param ic ...
3481 : ! **************************************************************************************************
3482 17820 : SUBROUTINE ctrig(n, trig, after, before, now, isign, ic)
3483 : INTEGER, INTENT(IN) :: n
3484 : REAL(dp), DIMENSION(2, ctrig_length), INTENT(OUT) :: trig
3485 : INTEGER, DIMENSION(7), INTENT(OUT) :: after, before, now
3486 : INTEGER, INTENT(IN) :: isign
3487 : INTEGER, INTENT(OUT) :: ic
3488 :
3489 : INTEGER, PARAMETER :: nt = 82
3490 : INTEGER, DIMENSION(7, nt), PARAMETER :: idata = RESHAPE((/3, 3, 1, 1, 1, 1, 1, 4, 4, 1, 1, 1,&
3491 : 1, 1, 5, 5, 1, 1, 1, 1, 1, 6, 6, 1, 1, 1, 1, 1, 8, 8, 1, 1, 1, 1, 1, 9, 3, 3, 1, 1, 1, 1, &
3492 : 12, 4, 3, 1, 1, 1, 1, 15, 5, 3, 1, 1, 1, 1, 16, 4, 4, 1, 1, 1, 1, 18, 6, 3, 1, 1, 1, 1, 20&
3493 : , 5, 4, 1, 1, 1, 1, 24, 8, 3, 1, 1, 1, 1, 25, 5, 5, 1, 1, 1, 1, 27, 3, 3, 3, 1, 1, 1, 30, &
3494 : 6, 5, 1, 1, 1, 1, 32, 8, 4, 1, 1, 1, 1, 36, 4, 3, 3, 1, 1, 1, 40, 8, 5, 1, 1, 1, 1, 45, 5 &
3495 : , 3, 3, 1, 1, 1, 48, 4, 4, 3, 1, 1, 1, 54, 6, 3, 3, 1, 1, 1, 60, 5, 4, 3, 1, 1, 1, 64, 4, &
3496 : 4, 4, 1, 1, 1, 72, 8, 3, 3, 1, 1, 1, 75, 5, 5, 3, 1, 1, 1, 80, 5, 4, 4, 1, 1, 1, 81, 3, 3 &
3497 : , 3, 3, 1, 1, 90, 6, 5, 3, 1, 1, 1, 96, 8, 4, 3, 1, 1, 1, 100, 5, 5, 4, 1, 1, 1, 108, 4, 3&
3498 : , 3, 3, 1, 1, 120, 8, 5, 3, 1, 1, 1, 125, 5, 5, 5, 1, 1, 1, 128, 8, 4, 4, 1, 1, 1, 135, 5 &
3499 : , 3, 3, 3, 1, 1, 144, 4, 4, 3, 3, 1, 1, 150, 6, 5, 5, 1, 1, 1, 160, 8, 5, 4, 1, 1, 1, 162,&
3500 : 6, 3, 3, 3, 1, 1, 180, 5, 4, 3, 3, 1, 1, 192, 4, 4, 4, 3, 1, 1, 200, 8, 5, 5, 1, 1, 1, 216&
3501 : , 8, 3, 3, 3, 1, 1, 225, 5, 5, 3, 3, 1, 1, 240, 5, 4, 4, 3, 1, 1, 243, 3, 3, 3, 3, 3, 1, &
3502 : 256, 4, 4, 4, 4, 1, 1, 270, 6, 5, 3, 3, 1, 1, 288, 8, 4, 3, 3, 1, 1, 300, 5, 5, 4, 3, 1, 1&
3503 : , 320, 5, 4, 4, 4, 1, 1, 324, 4, 3, 3, 3, 3, 1, 360, 8, 5, 3, 3, 1, 1, 375, 5, 5, 5, 3, 1,&
3504 : 1, 384, 8, 4, 4, 3, 1, 1, 400, 5, 5, 4, 4, 1, 1, 405, 5, 3, 3, 3, 3, 1, 432, 4, 4, 3, 3, 3&
3505 : , 1, 450, 6, 5, 5, 3, 1, 1, 480, 8, 5, 4, 3, 1, 1, 486, 6, 3, 3, 3, 3, 1, 500, 5, 5, 5, 4,&
3506 : 1, 1, 512, 8, 4, 4, 4, 1, 1, 540, 5, 4, 3, 3, 3, 1, 576, 4, 4, 4, 3, 3, 1, 600, 8, 5, 5, 3&
3507 : , 1, 1, 625, 5, 5, 5, 5, 1, 1, 640, 8, 5, 4, 4, 1, 1, 648, 8, 3, 3, 3, 3, 1, 675, 5, 5, 3,&
3508 : 3, 3, 1, 720, 5, 4, 4, 3, 3, 1, 729, 3, 3, 3, 3, 3, 3, 750, 6, 5, 5, 5, 1, 1, 768, 4, 4, 4&
3509 : , 4, 3, 1, 800, 8, 5, 5, 4, 1, 1, 810, 6, 5, 3, 3, 3, 1, 864, 8, 4, 3, 3, 3, 1, 900, 5, 5,&
3510 : 4, 3, 3, 1, 960, 5, 4, 4, 4, 3, 1, 972, 4, 3, 3, 3, 3, 3, 1000, 8, 5, 5, 5, 1, 1, &
3511 : ctrig_length, 4, 4, 4, 4, 4, 1/), (/7, nt/))
3512 :
3513 : INTEGER :: i, itt, j
3514 : REAL(dp) :: angle, twopi
3515 :
3516 194736 : mloop: DO i = 1, nt
3517 194736 : IF (n == idata(1, i)) THEN
3518 17820 : ic = 0
3519 52524 : DO j = 1, 6
3520 52524 : itt = idata(1 + j, i)
3521 52524 : IF (itt > 1) THEN
3522 34704 : ic = ic + 1
3523 34704 : now(j) = idata(1 + j, i)
3524 : ELSE
3525 : EXIT mloop
3526 : END IF
3527 : END DO
3528 : EXIT mloop
3529 : END IF
3530 194736 : IF (i == nt) THEN
3531 0 : WRITE (*, '(A,i5,A)') " Value of ", n, &
3532 0 : " not allowed for fft, allowed values are:"
3533 0 : WRITE (*, '(15i5)') (idata(1, j), j=1, nt)
3534 0 : CPABORT('ctrig')
3535 : END IF
3536 : END DO mloop
3537 :
3538 17820 : after(1) = 1
3539 17820 : before(ic) = 1
3540 34704 : DO i = 2, ic
3541 16884 : after(i) = after(i - 1)*now(i - 1)
3542 34704 : before(ic - i + 1) = before(ic - i + 2)*now(ic - i + 2)
3543 : END DO
3544 :
3545 17820 : twopi = 8._dp*ATAN(1._dp)
3546 17820 : angle = isign*twopi/REAL(n, dp)
3547 17820 : trig(1, 1) = 1._dp
3548 17820 : trig(2, 1) = 0._dp
3549 389256 : DO i = 1, n - 1
3550 371436 : trig(1, i + 1) = COS(REAL(i, dp)*angle)
3551 389256 : trig(2, i + 1) = SIN(REAL(i, dp)*angle)
3552 : END DO
3553 :
3554 17820 : END SUBROUTINE ctrig
3555 :
3556 : ! **************************************************************************************************
3557 : !> \brief ...
3558 : !> \param n ...
3559 : !> \param m ...
3560 : !> \param a ...
3561 : !> \param lda ...
3562 : !> \param b ...
3563 : !> \param ldb ...
3564 : ! **************************************************************************************************
3565 1008 : SUBROUTINE matmov(n, m, a, lda, b, ldb)
3566 : INTEGER :: n, m, lda
3567 : COMPLEX(dp) :: a(lda, *)
3568 : INTEGER :: ldb
3569 : COMPLEX(dp) :: b(ldb, *)
3570 :
3571 114540 : b(1:n, 1:m) = a(1:n, 1:m)
3572 1008 : END SUBROUTINE matmov
3573 :
3574 : ! **************************************************************************************************
3575 : !> \brief ...
3576 : !> \param a ...
3577 : !> \param lda ...
3578 : !> \param m ...
3579 : !> \param n ...
3580 : !> \param b ...
3581 : !> \param ldb ...
3582 : ! **************************************************************************************************
3583 1200 : SUBROUTINE zgetmo(a, lda, m, n, b, ldb)
3584 : INTEGER :: lda, m, n
3585 : COMPLEX(dp) :: a(lda, n)
3586 : INTEGER :: ldb
3587 : COMPLEX(dp) :: b(ldb, m)
3588 :
3589 152400 : b(1:n, 1:m) = TRANSPOSE(a(1:m, 1:n))
3590 1200 : END SUBROUTINE zgetmo
3591 :
3592 : ! **************************************************************************************************
3593 : !> \brief ...
3594 : !> \param n ...
3595 : !> \param sc ...
3596 : !> \param a ...
3597 : ! **************************************************************************************************
3598 160118 : SUBROUTINE scaled(n, sc, a)
3599 : INTEGER :: n
3600 : REAL(dp) :: sc
3601 : COMPLEX(dp) :: a(n)
3602 :
3603 160118 : CALL dscal(n, sc, a, 1)
3604 :
3605 160118 : END SUBROUTINE scaled
3606 :
3607 : END MODULE mltfftsg_tools
|
