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