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