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