(git:374b731)
Loading...
Searching...
No Matches
mltfftsg_tools.F
Go to the documentation of this file.
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! **************************************************************************************************
9MODULE 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
23CONTAINS
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 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 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 length = 2*(cache_size/4 + 1)
61
62 isig = -isign
63 tscal = (abs(scale - 1._dp) > 1.e-12_dp)
64 CALL ctrig(n, trig, after, before, now, isig, ic)
65 lot = cache_size/(4*n)
66 lot = lot - mod(lot + 1, 2)
67 lot = max(1, lot)
68
69 ! initializations for serial mode
70 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!$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 DEALLOCATE (z)
135
136 IF (transb == 'N' .OR. transb == 'n') THEN
137 b(1:ldbx, m + 1:ldby) = cmplx(0._dp, 0._dp, dp)
138 b(n + 1:ldbx, 1:m) = cmplx(0._dp, 0._dp, dp)
139 ELSE
140 b(1:ldbx, n + 1:ldby) = cmplx(0._dp, 0._dp, dp)
141 b(m + 1:ldbx, 1:n) = cmplx(0._dp, 0._dp, dp)
142 END IF
143
144 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 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 REAL(dp), DIMENSION(:, :, :), POINTER :: zin_real, zout_real
172
173 CALL c_f_pointer(c_loc(zin), zin_real, (/2, mm, m/))
174 CALL c_f_pointer(c_loc(zout), zout_real, (/2, nn, n/))
175 CALL fftstp(mm, nfft, m, nn, n, zin_real, zout_real, trig, now, after, before, isign)
176
177 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 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 REAL(dp), DIMENSION(:, :, :), POINTER :: zin_real, zout_real
204
205 CALL c_f_pointer(c_loc(zin), zin_real, (/2, mm, m/))
206 CALL c_f_pointer(c_loc(zout), zout_real, (/2, nn, n/))
207
208 CALL fftpre(mm, nfft, m, nn, n, zin_real, zout_real, trig, now, after, before, isign)
209
210 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 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 REAL(dp), DIMENSION(:, :, :), POINTER :: zin_real, zout_real
238
239 CALL c_f_pointer(c_loc(zin), zin_real, (/2, mm, m/))
240 CALL c_f_pointer(c_loc(zout), zout_real, (/2, nn, n/))
241
242 CALL fftrot(mm, nfft, m, nn, n, zin_real, zout_real, trig, now, after, before, isign)
243
244 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 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 atn = after*now
302 atb = after*before
303
304 IF (now == 4) THEN
305 IF (isign == 1) THEN
306 ia = 1
307 nin1 = ia - after
308 nout1 = ia - atn
309 DO ib = 1, before
310 nin1 = nin1 + after
311 nin2 = nin1 + atb
312 nin3 = nin2 + atb
313 nin4 = nin3 + atb
314 nout1 = nout1 + atn
315 nout2 = nout1 + after
316 nout3 = nout2 + after
317 nout4 = nout3 + after
318 DO j = 1, nfft
319 r1 = zin(1, j, nin1)
320 s1 = zin(2, j, nin1)
321 r2 = zin(1, j, nin2)
322 s2 = zin(2, j, nin2)
323 r3 = zin(1, j, nin3)
324 s3 = zin(2, j, nin3)
325 r4 = zin(1, j, nin4)
326 s4 = zin(2, j, nin4)
327 r = r1 + r3
328 s = r2 + r4
329 zout(1, nout1, j) = r + s
330 zout(1, nout3, j) = r - s
331 r = r1 - r3
332 s = s2 - s4
333 zout(1, nout2, j) = r - s
334 zout(1, nout4, j) = r + s
335 r = s1 + s3
336 s = s2 + s4
337 zout(2, nout1, j) = r + s
338 zout(2, nout3, j) = r - s
339 r = s1 - s3
340 s = r2 - r4
341 zout(2, nout2, j) = r + s
342 zout(2, nout4, j) = r - s
343 END DO
344 END DO
345 DO ia = 2, after
346 ias = ia - 1
347 IF (2*ias == after) THEN
348 nin1 = ia - after
349 nout1 = ia - atn
350 DO ib = 1, before
351 nin1 = nin1 + after
352 nin2 = nin1 + atb
353 nin3 = nin2 + atb
354 nin4 = nin3 + atb
355 nout1 = nout1 + atn
356 nout2 = nout1 + after
357 nout3 = nout2 + after
358 nout4 = nout3 + after
359 DO j = 1, nfft
360 r1 = zin(1, j, nin1)
361 s1 = zin(2, j, nin1)
362 r = zin(1, j, nin2)
363 s = zin(2, j, nin2)
364 r2 = (r - s)*rt2i
365 s2 = (r + s)*rt2i
366 r3 = -zin(2, j, nin3)
367 s3 = zin(1, j, nin3)
368 r = zin(1, j, nin4)
369 s = zin(2, j, nin4)
370 r4 = -(r + s)*rt2i
371 s4 = (r - s)*rt2i
372 r = r1 + r3
373 s = r2 + r4
374 zout(1, nout1, j) = r + s
375 zout(1, nout3, j) = r - s
376 r = r1 - r3
377 s = s2 - s4
378 zout(1, nout2, j) = r - s
379 zout(1, nout4, j) = r + s
380 r = s1 + s3
381 s = s2 + s4
382 zout(2, nout1, j) = r + s
383 zout(2, nout3, j) = r - s
384 r = s1 - s3
385 s = r2 - r4
386 zout(2, nout2, j) = r + s
387 zout(2, nout4, j) = r - s
388 END DO
389 END DO
390 ELSE
391 itt = ias*before
392 itrig = itt + 1
393 cr2 = trig(1, itrig)
394 ci2 = trig(2, itrig)
395 itrig = itrig + itt
396 cr3 = trig(1, itrig)
397 ci3 = trig(2, itrig)
398 itrig = itrig + itt
399 cr4 = trig(1, itrig)
400 ci4 = trig(2, itrig)
401 nin1 = ia - after
402 nout1 = ia - atn
403 DO ib = 1, before
404 nin1 = nin1 + after
405 nin2 = nin1 + atb
406 nin3 = nin2 + atb
407 nin4 = nin3 + atb
408 nout1 = nout1 + atn
409 nout2 = nout1 + after
410 nout3 = nout2 + after
411 nout4 = nout3 + after
412 DO j = 1, nfft
413 r1 = zin(1, j, nin1)
414 s1 = zin(2, j, nin1)
415 r = zin(1, j, nin2)
416 s = zin(2, j, nin2)
417 r2 = r*cr2 - s*ci2
418 s2 = r*ci2 + s*cr2
419 r = zin(1, j, nin3)
420 s = zin(2, j, nin3)
421 r3 = r*cr3 - s*ci3
422 s3 = r*ci3 + s*cr3
423 r = zin(1, j, nin4)
424 s = zin(2, j, nin4)
425 r4 = r*cr4 - s*ci4
426 s4 = r*ci4 + s*cr4
427 r = r1 + r3
428 s = r2 + r4
429 zout(1, nout1, j) = r + s
430 zout(1, nout3, j) = r - s
431 r = r1 - r3
432 s = s2 - s4
433 zout(1, nout2, j) = r - s
434 zout(1, nout4, j) = r + s
435 r = s1 + s3
436 s = s2 + s4
437 zout(2, nout1, j) = r + s
438 zout(2, nout3, j) = r - s
439 r = s1 - s3
440 s = r2 - r4
441 zout(2, nout2, j) = r + s
442 zout(2, nout4, j) = r - s
443 END DO
444 END DO
445 END IF
446 END DO
447 ELSE
448 ia = 1
449 nin1 = ia - after
450 nout1 = ia - atn
451 DO ib = 1, before
452 nin1 = nin1 + after
453 nin2 = nin1 + atb
454 nin3 = nin2 + atb
455 nin4 = nin3 + atb
456 nout1 = nout1 + atn
457 nout2 = nout1 + after
458 nout3 = nout2 + after
459 nout4 = nout3 + after
460 DO j = 1, nfft
461 r1 = zin(1, j, nin1)
462 s1 = zin(2, j, nin1)
463 r2 = zin(1, j, nin2)
464 s2 = zin(2, j, nin2)
465 r3 = zin(1, j, nin3)
466 s3 = zin(2, j, nin3)
467 r4 = zin(1, j, nin4)
468 s4 = zin(2, j, nin4)
469 r = r1 + r3
470 s = r2 + r4
471 zout(1, nout1, j) = r + s
472 zout(1, nout3, j) = r - s
473 r = r1 - r3
474 s = s2 - s4
475 zout(1, nout2, j) = r + s
476 zout(1, nout4, j) = r - s
477 r = s1 + s3
478 s = s2 + s4
479 zout(2, nout1, j) = r + s
480 zout(2, nout3, j) = r - s
481 r = s1 - s3
482 s = r2 - r4
483 zout(2, nout2, j) = r - s
484 zout(2, nout4, j) = r + s
485 END DO
486 END DO
487 DO ia = 2, after
488 ias = ia - 1
489 IF (2*ias == after) THEN
490 nin1 = ia - after
491 nout1 = ia - atn
492 DO ib = 1, before
493 nin1 = nin1 + after
494 nin2 = nin1 + atb
495 nin3 = nin2 + atb
496 nin4 = nin3 + atb
497 nout1 = nout1 + atn
498 nout2 = nout1 + after
499 nout3 = nout2 + after
500 nout4 = nout3 + after
501 DO j = 1, nfft
502 r1 = zin(1, j, nin1)
503 s1 = zin(2, j, nin1)
504 r = zin(1, j, nin2)
505 s = zin(2, j, nin2)
506 r2 = (r + s)*rt2i
507 s2 = (s - r)*rt2i
508 r3 = zin(2, j, nin3)
509 s3 = -zin(1, j, nin3)
510 r = zin(1, j, nin4)
511 s = zin(2, j, nin4)
512 r4 = (s - r)*rt2i
513 s4 = -(r + s)*rt2i
514 r = r1 + r3
515 s = r2 + r4
516 zout(1, nout1, j) = r + s
517 zout(1, nout3, j) = r - s
518 r = r1 - r3
519 s = s2 - s4
520 zout(1, nout2, j) = r + s
521 zout(1, nout4, j) = r - s
522 r = s1 + s3
523 s = s2 + s4
524 zout(2, nout1, j) = r + s
525 zout(2, nout3, j) = r - s
526 r = s1 - s3
527 s = r2 - r4
528 zout(2, nout2, j) = r - s
529 zout(2, nout4, j) = r + s
530 END DO
531 END DO
532 ELSE
533 itt = ias*before
534 itrig = itt + 1
535 cr2 = trig(1, itrig)
536 ci2 = trig(2, itrig)
537 itrig = itrig + itt
538 cr3 = trig(1, itrig)
539 ci3 = trig(2, itrig)
540 itrig = itrig + itt
541 cr4 = trig(1, itrig)
542 ci4 = trig(2, itrig)
543 nin1 = ia - after
544 nout1 = ia - atn
545 DO ib = 1, before
546 nin1 = nin1 + after
547 nin2 = nin1 + atb
548 nin3 = nin2 + atb
549 nin4 = nin3 + atb
550 nout1 = nout1 + atn
551 nout2 = nout1 + after
552 nout3 = nout2 + after
553 nout4 = nout3 + after
554 DO j = 1, nfft
555 r1 = zin(1, j, nin1)
556 s1 = zin(2, j, nin1)
557 r = zin(1, j, nin2)
558 s = zin(2, j, nin2)
559 r2 = r*cr2 - s*ci2
560 s2 = r*ci2 + s*cr2
561 r = zin(1, j, nin3)
562 s = zin(2, j, nin3)
563 r3 = r*cr3 - s*ci3
564 s3 = r*ci3 + s*cr3
565 r = zin(1, j, nin4)
566 s = zin(2, j, nin4)
567 r4 = r*cr4 - s*ci4
568 s4 = r*ci4 + s*cr4
569 r = r1 + r3
570 s = r2 + r4
571 zout(1, nout1, j) = r + s
572 zout(1, nout3, j) = r - s
573 r = r1 - r3
574 s = s2 - s4
575 zout(1, nout2, j) = r + s
576 zout(1, nout4, j) = r - s
577 r = s1 + s3
578 s = s2 + s4
579 zout(2, nout1, j) = r + s
580 zout(2, nout3, j) = r - s
581 r = s1 - s3
582 s = r2 - r4
583 zout(2, nout2, j) = r - s
584 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 IF (isign == -1) THEN
592 ia = 1
593 nin1 = ia - after
594 nout1 = ia - atn
595 DO ib = 1, before
596 nin1 = nin1 + after
597 nin2 = nin1 + atb
598 nin3 = nin2 + atb
599 nin4 = nin3 + atb
600 nin5 = nin4 + atb
601 nin6 = nin5 + atb
602 nin7 = nin6 + atb
603 nin8 = nin7 + atb
604 nout1 = nout1 + atn
605 nout2 = nout1 + after
606 nout3 = nout2 + after
607 nout4 = nout3 + after
608 nout5 = nout4 + after
609 nout6 = nout5 + after
610 nout7 = nout6 + after
611 nout8 = nout7 + after
612 DO j = 1, nfft
613 r1 = zin(1, j, nin1)
614 s1 = zin(2, j, nin1)
615 r2 = zin(1, j, nin2)
616 s2 = zin(2, j, nin2)
617 r3 = zin(1, j, nin3)
618 s3 = zin(2, j, nin3)
619 r4 = zin(1, j, nin4)
620 s4 = zin(2, j, nin4)
621 r5 = zin(1, j, nin5)
622 s5 = zin(2, j, nin5)
623 r6 = zin(1, j, nin6)
624 s6 = zin(2, j, nin6)
625 r7 = zin(1, j, nin7)
626 s7 = zin(2, j, nin7)
627 r8 = zin(1, j, nin8)
628 s8 = zin(2, j, nin8)
629 r = r1 + r5
630 s = r3 + r7
631 ap = r + s
632 am = r - s
633 r = r2 + r6
634 s = r4 + r8
635 bp = r + s
636 bm = r - s
637 r = s1 + s5
638 s = s3 + s7
639 cp = r + s
640 cm = r - s
641 r = s2 + s6
642 s = s4 + s8
643 dbl = r + s
644 dm = r - s
645 zout(1, nout1, j) = ap + bp
646 zout(2, nout1, j) = cp + dbl
647 zout(1, nout5, j) = ap - bp
648 zout(2, nout5, j) = cp - dbl
649 zout(1, nout3, j) = am + dm
650 zout(2, nout3, j) = cm - bm
651 zout(1, nout7, j) = am - dm
652 zout(2, nout7, j) = cm + bm
653 r = r1 - r5
654 s = s3 - s7
655 ap = r + s
656 am = r - s
657 r = s1 - s5
658 s = r3 - r7
659 bp = r + s
660 bm = r - s
661 r = s4 - s8
662 s = r2 - r6
663 cp = r + s
664 cm = r - s
665 r = s2 - s6
666 s = r4 - r8
667 dbl = r + s
668 dm = r - s
669 r = (cp + dm)*rt2i
670 s = (-cp + dm)*rt2i
671 cp = (cm + dbl)*rt2i
672 dbl = (cm - dbl)*rt2i
673 zout(1, nout2, j) = ap + r
674 zout(2, nout2, j) = bm + s
675 zout(1, nout6, j) = ap - r
676 zout(2, nout6, j) = bm - s
677 zout(1, nout4, j) = am + cp
678 zout(2, nout4, j) = bp + dbl
679 zout(1, nout8, j) = am - cp
680 zout(2, nout8, j) = bp - dbl
681 END DO
682 END DO
683 ELSE
684 ia = 1
685 nin1 = ia - after
686 nout1 = ia - atn
687 DO ib = 1, before
688 nin1 = nin1 + after
689 nin2 = nin1 + atb
690 nin3 = nin2 + atb
691 nin4 = nin3 + atb
692 nin5 = nin4 + atb
693 nin6 = nin5 + atb
694 nin7 = nin6 + atb
695 nin8 = nin7 + atb
696 nout1 = nout1 + atn
697 nout2 = nout1 + after
698 nout3 = nout2 + after
699 nout4 = nout3 + after
700 nout5 = nout4 + after
701 nout6 = nout5 + after
702 nout7 = nout6 + after
703 nout8 = nout7 + after
704 DO j = 1, nfft
705 r1 = zin(1, j, nin1)
706 s1 = zin(2, j, nin1)
707 r2 = zin(1, j, nin2)
708 s2 = zin(2, j, nin2)
709 r3 = zin(1, j, nin3)
710 s3 = zin(2, j, nin3)
711 r4 = zin(1, j, nin4)
712 s4 = zin(2, j, nin4)
713 r5 = zin(1, j, nin5)
714 s5 = zin(2, j, nin5)
715 r6 = zin(1, j, nin6)
716 s6 = zin(2, j, nin6)
717 r7 = zin(1, j, nin7)
718 s7 = zin(2, j, nin7)
719 r8 = zin(1, j, nin8)
720 s8 = zin(2, j, nin8)
721 r = r1 + r5
722 s = r3 + r7
723 ap = r + s
724 am = r - s
725 r = r2 + r6
726 s = r4 + r8
727 bp = r + s
728 bm = r - s
729 r = s1 + s5
730 s = s3 + s7
731 cp = r + s
732 cm = r - s
733 r = s2 + s6
734 s = s4 + s8
735 dbl = r + s
736 dm = r - s
737 zout(1, nout1, j) = ap + bp
738 zout(2, nout1, j) = cp + dbl
739 zout(1, nout5, j) = ap - bp
740 zout(2, nout5, j) = cp - dbl
741 zout(1, nout3, j) = am - dm
742 zout(2, nout3, j) = cm + bm
743 zout(1, nout7, j) = am + dm
744 zout(2, nout7, j) = cm - bm
745 r = r1 - r5
746 s = -s3 + s7
747 ap = r + s
748 am = r - s
749 r = s1 - s5
750 s = r7 - r3
751 bp = r + s
752 bm = r - s
753 r = -s4 + s8
754 s = r2 - r6
755 cp = r + s
756 cm = r - s
757 r = -s2 + s6
758 s = r4 - r8
759 dbl = r + s
760 dm = r - s
761 r = (cp + dm)*rt2i
762 s = (cp - dm)*rt2i
763 cp = (cm + dbl)*rt2i
764 dbl = (-cm + dbl)*rt2i
765 zout(1, nout2, j) = ap + r
766 zout(2, nout2, j) = bm + s
767 zout(1, nout6, j) = ap - r
768 zout(2, nout6, j) = bm - s
769 zout(1, nout4, j) = am + cp
770 zout(2, nout4, j) = bp + dbl
771 zout(1, nout8, j) = am - cp
772 zout(2, nout8, j) = bp - dbl
773 END DO
774 END DO
775 END IF
776 ELSE IF (now == 3) THEN
777 bbs = isign*bb
778 ia = 1
779 nin1 = ia - after
780 nout1 = ia - atn
781 DO ib = 1, before
782 nin1 = nin1 + after
783 nin2 = nin1 + atb
784 nin3 = nin2 + atb
785 nout1 = nout1 + atn
786 nout2 = nout1 + after
787 nout3 = nout2 + after
788 DO j = 1, nfft
789 r1 = zin(1, j, nin1)
790 s1 = zin(2, j, nin1)
791 r2 = zin(1, j, nin2)
792 s2 = zin(2, j, nin2)
793 r3 = zin(1, j, nin3)
794 s3 = zin(2, j, nin3)
795 r = r2 + r3
796 s = s2 + s3
797 zout(1, nout1, j) = r + r1
798 zout(2, nout1, j) = s + s1
799 r1 = r1 - 0.5_dp*r
800 s1 = s1 - 0.5_dp*s
801 r2 = bbs*(r2 - r3)
802 s2 = bbs*(s2 - s3)
803 zout(1, nout2, j) = r1 - s2
804 zout(2, nout2, j) = s1 + r2
805 zout(1, nout3, j) = r1 + s2
806 zout(2, nout3, j) = s1 - r2
807 END DO
808 END DO
809 DO ia = 2, after
810 ias = ia - 1
811 IF (4*ias == 3*after) THEN
812 IF (isign == 1) THEN
813 nin1 = ia - after
814 nout1 = ia - atn
815 DO ib = 1, before
816 nin1 = nin1 + after
817 nin2 = nin1 + atb
818 nin3 = nin2 + atb
819 nout1 = nout1 + atn
820 nout2 = nout1 + after
821 nout3 = nout2 + after
822 DO j = 1, nfft
823 r1 = zin(1, j, nin1)
824 s1 = zin(2, j, nin1)
825 r2 = -zin(2, j, nin2)
826 s2 = zin(1, j, nin2)
827 r3 = -zin(1, j, nin3)
828 s3 = -zin(2, j, nin3)
829 r = r2 + r3
830 s = s2 + s3
831 zout(1, nout1, j) = r + r1
832 zout(2, nout1, j) = s + s1
833 r1 = r1 - 0.5_dp*r
834 s1 = s1 - 0.5_dp*s
835 r2 = bbs*(r2 - r3)
836 s2 = bbs*(s2 - s3)
837 zout(1, nout2, j) = r1 - s2
838 zout(2, nout2, j) = s1 + r2
839 zout(1, nout3, j) = r1 + s2
840 zout(2, nout3, j) = s1 - r2
841 END DO
842 END DO
843 ELSE
844 nin1 = ia - after
845 nout1 = ia - atn
846 DO ib = 1, before
847 nin1 = nin1 + after
848 nin2 = nin1 + atb
849 nin3 = nin2 + atb
850 nout1 = nout1 + atn
851 nout2 = nout1 + after
852 nout3 = nout2 + after
853 DO j = 1, nfft
854 r1 = zin(1, j, nin1)
855 s1 = zin(2, j, nin1)
856 r2 = zin(2, j, nin2)
857 s2 = -zin(1, j, nin2)
858 r3 = -zin(1, j, nin3)
859 s3 = -zin(2, j, nin3)
860 r = r2 + r3
861 s = s2 + s3
862 zout(1, nout1, j) = r + r1
863 zout(2, nout1, j) = s + s1
864 r1 = r1 - 0.5_dp*r
865 s1 = s1 - 0.5_dp*s
866 r2 = bbs*(r2 - r3)
867 s2 = bbs*(s2 - s3)
868 zout(1, nout2, j) = r1 - s2
869 zout(2, nout2, j) = s1 + r2
870 zout(1, nout3, j) = r1 + s2
871 zout(2, nout3, j) = s1 - r2
872 END DO
873 END DO
874 END IF
875 ELSE IF (8*ias == 3*after) THEN
876 IF (isign == 1) THEN
877 nin1 = ia - after
878 nout1 = ia - atn
879 DO ib = 1, before
880 nin1 = nin1 + after
881 nin2 = nin1 + atb
882 nin3 = nin2 + atb
883 nout1 = nout1 + atn
884 nout2 = nout1 + after
885 nout3 = nout2 + after
886 DO j = 1, nfft
887 r1 = zin(1, j, nin1)
888 s1 = zin(2, j, nin1)
889 r = zin(1, j, nin2)
890 s = zin(2, j, nin2)
891 r2 = (r - s)*rt2i
892 s2 = (r + s)*rt2i
893 r3 = -zin(2, j, nin3)
894 s3 = zin(1, j, nin3)
895 r = r2 + r3
896 s = s2 + s3
897 zout(1, nout1, j) = r + r1
898 zout(2, nout1, j) = s + s1
899 r1 = r1 - 0.5_dp*r
900 s1 = s1 - 0.5_dp*s
901 r2 = bbs*(r2 - r3)
902 s2 = bbs*(s2 - s3)
903 zout(1, nout2, j) = r1 - s2
904 zout(2, nout2, j) = s1 + r2
905 zout(1, nout3, j) = r1 + s2
906 zout(2, nout3, j) = s1 - r2
907 END DO
908 END DO
909 ELSE
910 nin1 = ia - after
911 nout1 = ia - atn
912 DO ib = 1, before
913 nin1 = nin1 + after
914 nin2 = nin1 + atb
915 nin3 = nin2 + atb
916 nout1 = nout1 + atn
917 nout2 = nout1 + after
918 nout3 = nout2 + after
919 DO j = 1, nfft
920 r1 = zin(1, j, nin1)
921 s1 = zin(2, j, nin1)
922 r = zin(1, j, nin2)
923 s = zin(2, j, nin2)
924 r2 = (r + s)*rt2i
925 s2 = (-r + s)*rt2i
926 r3 = zin(2, j, nin3)
927 s3 = -zin(1, j, nin3)
928 r = r2 + r3
929 s = s2 + s3
930 zout(1, nout1, j) = r + r1
931 zout(2, nout1, j) = s + s1
932 r1 = r1 - 0.5_dp*r
933 s1 = s1 - 0.5_dp*s
934 r2 = bbs*(r2 - r3)
935 s2 = bbs*(s2 - s3)
936 zout(1, nout2, j) = r1 - s2
937 zout(2, nout2, j) = s1 + r2
938 zout(1, nout3, j) = r1 + s2
939 zout(2, nout3, j) = s1 - r2
940 END DO
941 END DO
942 END IF
943 ELSE
944 itt = ias*before
945 itrig = itt + 1
946 cr2 = trig(1, itrig)
947 ci2 = trig(2, itrig)
948 itrig = itrig + itt
949 cr3 = trig(1, itrig)
950 ci3 = trig(2, itrig)
951 nin1 = ia - after
952 nout1 = ia - atn
953 DO ib = 1, before
954 nin1 = nin1 + after
955 nin2 = nin1 + atb
956 nin3 = nin2 + atb
957 nout1 = nout1 + atn
958 nout2 = nout1 + after
959 nout3 = nout2 + after
960 DO j = 1, nfft
961 r1 = zin(1, j, nin1)
962 s1 = zin(2, j, nin1)
963 r = zin(1, j, nin2)
964 s = zin(2, j, nin2)
965 r2 = r*cr2 - s*ci2
966 s2 = r*ci2 + s*cr2
967 r = zin(1, j, nin3)
968 s = zin(2, j, nin3)
969 r3 = r*cr3 - s*ci3
970 s3 = r*ci3 + s*cr3
971 r = r2 + r3
972 s = s2 + s3
973 zout(1, nout1, j) = r + r1
974 zout(2, nout1, j) = s + s1
975 r1 = r1 - 0.5_dp*r
976 s1 = s1 - 0.5_dp*s
977 r2 = bbs*(r2 - r3)
978 s2 = bbs*(s2 - s3)
979 zout(1, nout2, j) = r1 - s2
980 zout(2, nout2, j) = s1 + r2
981 zout(1, nout3, j) = r1 + s2
982 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 sin2 = isign*sin2p
989 sin4 = isign*sin4p
990 ia = 1
991 nin1 = ia - after
992 nout1 = ia - atn
993 DO ib = 1, before
994 nin1 = nin1 + after
995 nin2 = nin1 + atb
996 nin3 = nin2 + atb
997 nin4 = nin3 + atb
998 nin5 = nin4 + atb
999 nout1 = nout1 + atn
1000 nout2 = nout1 + after
1001 nout3 = nout2 + after
1002 nout4 = nout3 + after
1003 nout5 = nout4 + after
1004 DO j = 1, nfft
1005 r1 = zin(1, j, nin1)
1006 s1 = zin(2, j, nin1)
1007 r2 = zin(1, j, nin2)
1008 s2 = zin(2, j, nin2)
1009 r3 = zin(1, j, nin3)
1010 s3 = zin(2, j, nin3)
1011 r4 = zin(1, j, nin4)
1012 s4 = zin(2, j, nin4)
1013 r5 = zin(1, j, nin5)
1014 s5 = zin(2, j, nin5)
1015 r25 = r2 + r5
1016 r34 = r3 + r4
1017 s25 = s2 - s5
1018 s34 = s3 - s4
1019 zout(1, nout1, j) = r1 + r25 + r34
1020 r = cos2*r25 + cos4*r34 + r1
1021 s = sin2*s25 + sin4*s34
1022 zout(1, nout2, j) = r - s
1023 zout(1, nout5, j) = r + s
1024 r = cos4*r25 + cos2*r34 + r1
1025 s = sin4*s25 - sin2*s34
1026 zout(1, nout3, j) = r - s
1027 zout(1, nout4, j) = r + s
1028 r25 = r2 - r5
1029 r34 = r3 - r4
1030 s25 = s2 + s5
1031 s34 = s3 + s4
1032 zout(2, nout1, j) = s1 + s25 + s34
1033 r = cos2*s25 + cos4*s34 + s1
1034 s = sin2*r25 + sin4*r34
1035 zout(2, nout2, j) = r + s
1036 zout(2, nout5, j) = r - s
1037 r = cos4*s25 + cos2*s34 + s1
1038 s = sin4*r25 - sin2*r34
1039 zout(2, nout3, j) = r + s
1040 zout(2, nout4, j) = r - s
1041 END DO
1042 END DO
1043 DO ia = 2, after
1044 ias = ia - 1
1045 IF (8*ias == 5*after) THEN
1046 IF (isign == 1) THEN
1047 nin1 = ia - after
1048 nout1 = ia - atn
1049 DO ib = 1, before
1050 nin1 = nin1 + after
1051 nin2 = nin1 + atb
1052 nin3 = nin2 + atb
1053 nin4 = nin3 + atb
1054 nin5 = nin4 + atb
1055 nout1 = nout1 + atn
1056 nout2 = nout1 + after
1057 nout3 = nout2 + after
1058 nout4 = nout3 + after
1059 nout5 = nout4 + after
1060 DO j = 1, nfft
1061 r1 = zin(1, j, nin1)
1062 s1 = zin(2, j, nin1)
1063 r = zin(1, j, nin2)
1064 s = zin(2, j, nin2)
1065 r2 = (r - s)*rt2i
1066 s2 = (r + s)*rt2i
1067 r3 = -zin(2, j, nin3)
1068 s3 = zin(1, j, nin3)
1069 r = zin(1, j, nin4)
1070 s = zin(2, j, nin4)
1071 r4 = -(r + s)*rt2i
1072 s4 = (r - s)*rt2i
1073 r5 = -zin(1, j, nin5)
1074 s5 = -zin(2, j, nin5)
1075 r25 = r2 + r5
1076 r34 = r3 + r4
1077 s25 = s2 - s5
1078 s34 = s3 - s4
1079 zout(1, nout1, j) = r1 + r25 + r34
1080 r = cos2*r25 + cos4*r34 + r1
1081 s = sin2*s25 + sin4*s34
1082 zout(1, nout2, j) = r - s
1083 zout(1, nout5, j) = r + s
1084 r = cos4*r25 + cos2*r34 + r1
1085 s = sin4*s25 - sin2*s34
1086 zout(1, nout3, j) = r - s
1087 zout(1, nout4, j) = r + s
1088 r25 = r2 - r5
1089 r34 = r3 - r4
1090 s25 = s2 + s5
1091 s34 = s3 + s4
1092 zout(2, nout1, j) = s1 + s25 + s34
1093 r = cos2*s25 + cos4*s34 + s1
1094 s = sin2*r25 + sin4*r34
1095 zout(2, nout2, j) = r + s
1096 zout(2, nout5, j) = r - s
1097 r = cos4*s25 + cos2*s34 + s1
1098 s = sin4*r25 - sin2*r34
1099 zout(2, nout3, j) = r + s
1100 zout(2, nout4, j) = r - s
1101 END DO
1102 END DO
1103 ELSE
1104 nin1 = ia - after
1105 nout1 = ia - atn
1106 DO ib = 1, before
1107 nin1 = nin1 + after
1108 nin2 = nin1 + atb
1109 nin3 = nin2 + atb
1110 nin4 = nin3 + atb
1111 nin5 = nin4 + atb
1112 nout1 = nout1 + atn
1113 nout2 = nout1 + after
1114 nout3 = nout2 + after
1115 nout4 = nout3 + after
1116 nout5 = nout4 + after
1117 DO j = 1, nfft
1118 r1 = zin(1, j, nin1)
1119 s1 = zin(2, j, nin1)
1120 r = zin(1, j, nin2)
1121 s = zin(2, j, nin2)
1122 r2 = (r + s)*rt2i
1123 s2 = (-r + s)*rt2i
1124 r3 = zin(2, j, nin3)
1125 s3 = -zin(1, j, nin3)
1126 r = zin(1, j, nin4)
1127 s = zin(2, j, nin4)
1128 r4 = (s - r)*rt2i
1129 s4 = -(r + s)*rt2i
1130 r5 = -zin(1, j, nin5)
1131 s5 = -zin(2, j, nin5)
1132 r25 = r2 + r5
1133 r34 = r3 + r4
1134 s25 = s2 - s5
1135 s34 = s3 - s4
1136 zout(1, nout1, j) = r1 + r25 + r34
1137 r = cos2*r25 + cos4*r34 + r1
1138 s = sin2*s25 + sin4*s34
1139 zout(1, nout2, j) = r - s
1140 zout(1, nout5, j) = r + s
1141 r = cos4*r25 + cos2*r34 + r1
1142 s = sin4*s25 - sin2*s34
1143 zout(1, nout3, j) = r - s
1144 zout(1, nout4, j) = r + s
1145 r25 = r2 - r5
1146 r34 = r3 - r4
1147 s25 = s2 + s5
1148 s34 = s3 + s4
1149 zout(2, nout1, j) = s1 + s25 + s34
1150 r = cos2*s25 + cos4*s34 + s1
1151 s = sin2*r25 + sin4*r34
1152 zout(2, nout2, j) = r + s
1153 zout(2, nout5, j) = r - s
1154 r = cos4*s25 + cos2*s34 + s1
1155 s = sin4*r25 - sin2*r34
1156 zout(2, nout3, j) = r + s
1157 zout(2, nout4, j) = r - s
1158 END DO
1159 END DO
1160 END IF
1161 ELSE
1162 ias = ia - 1
1163 itt = ias*before
1164 itrig = itt + 1
1165 cr2 = trig(1, itrig)
1166 ci2 = trig(2, itrig)
1167 itrig = itrig + itt
1168 cr3 = trig(1, itrig)
1169 ci3 = trig(2, itrig)
1170 itrig = itrig + itt
1171 cr4 = trig(1, itrig)
1172 ci4 = trig(2, itrig)
1173 itrig = itrig + itt
1174 cr5 = trig(1, itrig)
1175 ci5 = trig(2, itrig)
1176 nin1 = ia - after
1177 nout1 = ia - atn
1178 DO ib = 1, before
1179 nin1 = nin1 + after
1180 nin2 = nin1 + atb
1181 nin3 = nin2 + atb
1182 nin4 = nin3 + atb
1183 nin5 = nin4 + atb
1184 nout1 = nout1 + atn
1185 nout2 = nout1 + after
1186 nout3 = nout2 + after
1187 nout4 = nout3 + after
1188 nout5 = nout4 + after
1189 DO j = 1, nfft
1190 r1 = zin(1, j, nin1)
1191 s1 = zin(2, j, nin1)
1192 r = zin(1, j, nin2)
1193 s = zin(2, j, nin2)
1194 r2 = r*cr2 - s*ci2
1195 s2 = r*ci2 + s*cr2
1196 r = zin(1, j, nin3)
1197 s = zin(2, j, nin3)
1198 r3 = r*cr3 - s*ci3
1199 s3 = r*ci3 + s*cr3
1200 r = zin(1, j, nin4)
1201 s = zin(2, j, nin4)
1202 r4 = r*cr4 - s*ci4
1203 s4 = r*ci4 + s*cr4
1204 r = zin(1, j, nin5)
1205 s = zin(2, j, nin5)
1206 r5 = r*cr5 - s*ci5
1207 s5 = r*ci5 + s*cr5
1208 r25 = r2 + r5
1209 r34 = r3 + r4
1210 s25 = s2 - s5
1211 s34 = s3 - s4
1212 zout(1, nout1, j) = r1 + r25 + r34
1213 r = cos2*r25 + cos4*r34 + r1
1214 s = sin2*s25 + sin4*s34
1215 zout(1, nout2, j) = r - s
1216 zout(1, nout5, j) = r + s
1217 r = cos4*r25 + cos2*r34 + r1
1218 s = sin4*s25 - sin2*s34
1219 zout(1, nout3, j) = r - s
1220 zout(1, nout4, j) = r + s
1221 r25 = r2 - r5
1222 r34 = r3 - r4
1223 s25 = s2 + s5
1224 s34 = s3 + s4
1225 zout(2, nout1, j) = s1 + s25 + s34
1226 r = cos2*s25 + cos4*s34 + s1
1227 s = sin2*r25 + sin4*r34
1228 zout(2, nout2, j) = r + s
1229 zout(2, nout5, j) = r - s
1230 r = cos4*s25 + cos2*s34 + s1
1231 s = sin4*r25 - sin2*r34
1232 zout(2, nout3, j) = r + s
1233 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 bbs = isign*bb
1240 ia = 1
1241 nin1 = ia - after
1242 nout1 = ia - atn
1243 DO ib = 1, before
1244 nin1 = nin1 + after
1245 nin2 = nin1 + atb
1246 nin3 = nin2 + atb
1247 nin4 = nin3 + atb
1248 nin5 = nin4 + atb
1249 nin6 = nin5 + atb
1250 nout1 = nout1 + atn
1251 nout2 = nout1 + after
1252 nout3 = nout2 + after
1253 nout4 = nout3 + after
1254 nout5 = nout4 + after
1255 nout6 = nout5 + after
1256 DO j = 1, nfft
1257 r2 = zin(1, j, nin3)
1258 s2 = zin(2, j, nin3)
1259 r3 = zin(1, j, nin5)
1260 s3 = zin(2, j, nin5)
1261 r = r2 + r3
1262 s = s2 + s3
1263 r1 = zin(1, j, nin1)
1264 s1 = zin(2, j, nin1)
1265 ur1 = r + r1
1266 ui1 = s + s1
1267 r1 = r1 - 0.5_dp*r
1268 s1 = s1 - 0.5_dp*s
1269 r = r2 - r3
1270 s = s2 - s3
1271 ur2 = r1 - s*bbs
1272 ui2 = s1 + r*bbs
1273 ur3 = r1 + s*bbs
1274 ui3 = s1 - r*bbs
1275
1276 r2 = zin(1, j, nin6)
1277 s2 = zin(2, j, nin6)
1278 r3 = zin(1, j, nin2)
1279 s3 = zin(2, j, nin2)
1280 r = r2 + r3
1281 s = s2 + s3
1282 r1 = zin(1, j, nin4)
1283 s1 = zin(2, j, nin4)
1284 vr1 = r + r1
1285 vi1 = s + s1
1286 r1 = r1 - 0.5_dp*r
1287 s1 = s1 - 0.5_dp*s
1288 r = r2 - r3
1289 s = s2 - s3
1290 vr2 = r1 - s*bbs
1291 vi2 = s1 + r*bbs
1292 vr3 = r1 + s*bbs
1293 vi3 = s1 - r*bbs
1294
1295 zout(1, nout1, j) = ur1 + vr1
1296 zout(2, nout1, j) = ui1 + vi1
1297 zout(1, nout5, j) = ur2 + vr2
1298 zout(2, nout5, j) = ui2 + vi2
1299 zout(1, nout3, j) = ur3 + vr3
1300 zout(2, nout3, j) = ui3 + vi3
1301 zout(1, nout4, j) = ur1 - vr1
1302 zout(2, nout4, j) = ui1 - vi1
1303 zout(1, nout2, j) = ur2 - vr2
1304 zout(2, nout2, j) = ui2 - vi2
1305 zout(1, nout6, j) = ur3 - vr3
1306 zout(2, nout6, j) = ui3 - vi3
1307 END DO
1308 END DO
1309 ELSE
1310 cpabort('Error fftrot')
1311 END IF
1312
1313!-----------------------------------------------------------------------------!
1314
1315 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 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 atn = after*now
1373 atb = after*before
1374
1375 IF (now == 4) THEN
1376 IF (isign == 1) THEN
1377 ia = 1
1378 nin1 = ia - after
1379 nout1 = ia - atn
1380 DO ib = 1, before
1381 nin1 = nin1 + after
1382 nin2 = nin1 + atb
1383 nin3 = nin2 + atb
1384 nin4 = nin3 + atb
1385 nout1 = nout1 + atn
1386 nout2 = nout1 + after
1387 nout3 = nout2 + after
1388 nout4 = nout3 + after
1389 DO j = 1, nfft
1390 r1 = zin(1, nin1, j)
1391 s1 = zin(2, nin1, j)
1392 r2 = zin(1, nin2, j)
1393 s2 = zin(2, nin2, j)
1394 r3 = zin(1, nin3, j)
1395 s3 = zin(2, nin3, j)
1396 r4 = zin(1, nin4, j)
1397 s4 = zin(2, nin4, j)
1398 r = r1 + r3
1399 s = r2 + r4
1400 zout(1, j, nout1) = r + s
1401 zout(1, j, nout3) = r - s
1402 r = r1 - r3
1403 s = s2 - s4
1404 zout(1, j, nout2) = r - s
1405 zout(1, j, nout4) = r + s
1406 r = s1 + s3
1407 s = s2 + s4
1408 zout(2, j, nout1) = r + s
1409 zout(2, j, nout3) = r - s
1410 r = s1 - s3
1411 s = r2 - r4
1412 zout(2, j, nout2) = r + s
1413 zout(2, j, nout4) = r - s
1414 END DO
1415 END DO
1416 DO ia = 2, after
1417 ias = ia - 1
1418 IF (2*ias == after) THEN
1419 nin1 = ia - after
1420 nout1 = ia - atn
1421 DO ib = 1, before
1422 nin1 = nin1 + after
1423 nin2 = nin1 + atb
1424 nin3 = nin2 + atb
1425 nin4 = nin3 + atb
1426 nout1 = nout1 + atn
1427 nout2 = nout1 + after
1428 nout3 = nout2 + after
1429 nout4 = nout3 + after
1430 DO j = 1, nfft
1431 r1 = zin(1, nin1, j)
1432 s1 = zin(2, nin1, j)
1433 r = zin(1, nin2, j)
1434 s = zin(2, nin2, j)
1435 r2 = (r - s)*rt2i
1436 s2 = (r + s)*rt2i
1437 r3 = -zin(2, nin3, j)
1438 s3 = zin(1, nin3, j)
1439 r = zin(1, nin4, j)
1440 s = zin(2, nin4, j)
1441 r4 = -(r + s)*rt2i
1442 s4 = (r - s)*rt2i
1443 r = r1 + r3
1444 s = r2 + r4
1445 zout(1, j, nout1) = r + s
1446 zout(1, j, nout3) = r - s
1447 r = r1 - r3
1448 s = s2 - s4
1449 zout(1, j, nout2) = r - s
1450 zout(1, j, nout4) = r + s
1451 r = s1 + s3
1452 s = s2 + s4
1453 zout(2, j, nout1) = r + s
1454 zout(2, j, nout3) = r - s
1455 r = s1 - s3
1456 s = r2 - r4
1457 zout(2, j, nout2) = r + s
1458 zout(2, j, nout4) = r - s
1459 END DO
1460 END DO
1461 ELSE
1462 itt = ias*before
1463 itrig = itt + 1
1464 cr2 = trig(1, itrig)
1465 ci2 = trig(2, itrig)
1466 itrig = itrig + itt
1467 cr3 = trig(1, itrig)
1468 ci3 = trig(2, itrig)
1469 itrig = itrig + itt
1470 cr4 = trig(1, itrig)
1471 ci4 = trig(2, itrig)
1472 nin1 = ia - after
1473 nout1 = ia - atn
1474 DO ib = 1, before
1475 nin1 = nin1 + after
1476 nin2 = nin1 + atb
1477 nin3 = nin2 + atb
1478 nin4 = nin3 + atb
1479 nout1 = nout1 + atn
1480 nout2 = nout1 + after
1481 nout3 = nout2 + after
1482 nout4 = nout3 + after
1483 DO j = 1, nfft
1484 r1 = zin(1, nin1, j)
1485 s1 = zin(2, nin1, j)
1486 r = zin(1, nin2, j)
1487 s = zin(2, nin2, j)
1488 r2 = r*cr2 - s*ci2
1489 s2 = r*ci2 + s*cr2
1490 r = zin(1, nin3, j)
1491 s = zin(2, nin3, j)
1492 r3 = r*cr3 - s*ci3
1493 s3 = r*ci3 + s*cr3
1494 r = zin(1, nin4, j)
1495 s = zin(2, nin4, j)
1496 r4 = r*cr4 - s*ci4
1497 s4 = r*ci4 + s*cr4
1498 r = r1 + r3
1499 s = r2 + r4
1500 zout(1, j, nout1) = r + s
1501 zout(1, j, nout3) = r - s
1502 r = r1 - r3
1503 s = s2 - s4
1504 zout(1, j, nout2) = r - s
1505 zout(1, j, nout4) = r + s
1506 r = s1 + s3
1507 s = s2 + s4
1508 zout(2, j, nout1) = r + s
1509 zout(2, j, nout3) = r - s
1510 r = s1 - s3
1511 s = r2 - r4
1512 zout(2, j, nout2) = r + s
1513 zout(2, j, nout4) = r - s
1514 END DO
1515 END DO
1516 END IF
1517 END DO
1518 ELSE
1519 ia = 1
1520 nin1 = ia - after
1521 nout1 = ia - atn
1522 DO ib = 1, before
1523 nin1 = nin1 + after
1524 nin2 = nin1 + atb
1525 nin3 = nin2 + atb
1526 nin4 = nin3 + atb
1527 nout1 = nout1 + atn
1528 nout2 = nout1 + after
1529 nout3 = nout2 + after
1530 nout4 = nout3 + after
1531 DO j = 1, nfft
1532 r1 = zin(1, nin1, j)
1533 s1 = zin(2, nin1, j)
1534 r2 = zin(1, nin2, j)
1535 s2 = zin(2, nin2, j)
1536 r3 = zin(1, nin3, j)
1537 s3 = zin(2, nin3, j)
1538 r4 = zin(1, nin4, j)
1539 s4 = zin(2, nin4, j)
1540 r = r1 + r3
1541 s = r2 + r4
1542 zout(1, j, nout1) = r + s
1543 zout(1, j, nout3) = r - s
1544 r = r1 - r3
1545 s = s2 - s4
1546 zout(1, j, nout2) = r + s
1547 zout(1, j, nout4) = r - s
1548 r = s1 + s3
1549 s = s2 + s4
1550 zout(2, j, nout1) = r + s
1551 zout(2, j, nout3) = r - s
1552 r = s1 - s3
1553 s = r2 - r4
1554 zout(2, j, nout2) = r - s
1555 zout(2, j, nout4) = r + s
1556 END DO
1557 END DO
1558 DO ia = 2, after
1559 ias = ia - 1
1560 IF (2*ias == after) THEN
1561 nin1 = ia - after
1562 nout1 = ia - atn
1563 DO ib = 1, before
1564 nin1 = nin1 + after
1565 nin2 = nin1 + atb
1566 nin3 = nin2 + atb
1567 nin4 = nin3 + atb
1568 nout1 = nout1 + atn
1569 nout2 = nout1 + after
1570 nout3 = nout2 + after
1571 nout4 = nout3 + after
1572 DO j = 1, nfft
1573 r1 = zin(1, nin1, j)
1574 s1 = zin(2, nin1, j)
1575 r = zin(1, nin2, j)
1576 s = zin(2, nin2, j)
1577 r2 = (r + s)*rt2i
1578 s2 = (s - r)*rt2i
1579 r3 = zin(2, nin3, j)
1580 s3 = -zin(1, nin3, j)
1581 r = zin(1, nin4, j)
1582 s = zin(2, nin4, j)
1583 r4 = (s - r)*rt2i
1584 s4 = -(r + s)*rt2i
1585 r = r1 + r3
1586 s = r2 + r4
1587 zout(1, j, nout1) = r + s
1588 zout(1, j, nout3) = r - s
1589 r = r1 - r3
1590 s = s2 - s4
1591 zout(1, j, nout2) = r + s
1592 zout(1, j, nout4) = r - s
1593 r = s1 + s3
1594 s = s2 + s4
1595 zout(2, j, nout1) = r + s
1596 zout(2, j, nout3) = r - s
1597 r = s1 - s3
1598 s = r2 - r4
1599 zout(2, j, nout2) = r - s
1600 zout(2, j, nout4) = r + s
1601 END DO
1602 END DO
1603 ELSE
1604 itt = ias*before
1605 itrig = itt + 1
1606 cr2 = trig(1, itrig)
1607 ci2 = trig(2, itrig)
1608 itrig = itrig + itt
1609 cr3 = trig(1, itrig)
1610 ci3 = trig(2, itrig)
1611 itrig = itrig + itt
1612 cr4 = trig(1, itrig)
1613 ci4 = trig(2, itrig)
1614 nin1 = ia - after
1615 nout1 = ia - atn
1616 DO ib = 1, before
1617 nin1 = nin1 + after
1618 nin2 = nin1 + atb
1619 nin3 = nin2 + atb
1620 nin4 = nin3 + atb
1621 nout1 = nout1 + atn
1622 nout2 = nout1 + after
1623 nout3 = nout2 + after
1624 nout4 = nout3 + after
1625 DO j = 1, nfft
1626 r1 = zin(1, nin1, j)
1627 s1 = zin(2, nin1, j)
1628 r = zin(1, nin2, j)
1629 s = zin(2, nin2, j)
1630 r2 = r*cr2 - s*ci2
1631 s2 = r*ci2 + s*cr2
1632 r = zin(1, nin3, j)
1633 s = zin(2, nin3, j)
1634 r3 = r*cr3 - s*ci3
1635 s3 = r*ci3 + s*cr3
1636 r = zin(1, nin4, j)
1637 s = zin(2, nin4, j)
1638 r4 = r*cr4 - s*ci4
1639 s4 = r*ci4 + s*cr4
1640 r = r1 + r3
1641 s = r2 + r4
1642 zout(1, j, nout1) = r + s
1643 zout(1, j, nout3) = r - s
1644 r = r1 - r3
1645 s = s2 - s4
1646 zout(1, j, nout2) = r + s
1647 zout(1, j, nout4) = r - s
1648 r = s1 + s3
1649 s = s2 + s4
1650 zout(2, j, nout1) = r + s
1651 zout(2, j, nout3) = r - s
1652 r = s1 - s3
1653 s = r2 - r4
1654 zout(2, j, nout2) = r - s
1655 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 IF (isign == -1) THEN
1663 ia = 1
1664 nin1 = ia - after
1665 nout1 = ia - atn
1666 DO ib = 1, before
1667 nin1 = nin1 + after
1668 nin2 = nin1 + atb
1669 nin3 = nin2 + atb
1670 nin4 = nin3 + atb
1671 nin5 = nin4 + atb
1672 nin6 = nin5 + atb
1673 nin7 = nin6 + atb
1674 nin8 = nin7 + atb
1675 nout1 = nout1 + atn
1676 nout2 = nout1 + after
1677 nout3 = nout2 + after
1678 nout4 = nout3 + after
1679 nout5 = nout4 + after
1680 nout6 = nout5 + after
1681 nout7 = nout6 + after
1682 nout8 = nout7 + after
1683 DO j = 1, nfft
1684 r1 = zin(1, nin1, j)
1685 s1 = zin(2, nin1, j)
1686 r2 = zin(1, nin2, j)
1687 s2 = zin(2, nin2, j)
1688 r3 = zin(1, nin3, j)
1689 s3 = zin(2, nin3, j)
1690 r4 = zin(1, nin4, j)
1691 s4 = zin(2, nin4, j)
1692 r5 = zin(1, nin5, j)
1693 s5 = zin(2, nin5, j)
1694 r6 = zin(1, nin6, j)
1695 s6 = zin(2, nin6, j)
1696 r7 = zin(1, nin7, j)
1697 s7 = zin(2, nin7, j)
1698 r8 = zin(1, nin8, j)
1699 s8 = zin(2, nin8, j)
1700 r = r1 + r5
1701 s = r3 + r7
1702 ap = r + s
1703 am = r - s
1704 r = r2 + r6
1705 s = r4 + r8
1706 bp = r + s
1707 bm = r - s
1708 r = s1 + s5
1709 s = s3 + s7
1710 cp = r + s
1711 cm = r - s
1712 r = s2 + s6
1713 s = s4 + s8
1714 dbl = r + s
1715 dm = r - s
1716 zout(1, j, nout1) = ap + bp
1717 zout(2, j, nout1) = cp + dbl
1718 zout(1, j, nout5) = ap - bp
1719 zout(2, j, nout5) = cp - dbl
1720 zout(1, j, nout3) = am + dm
1721 zout(2, j, nout3) = cm - bm
1722 zout(1, j, nout7) = am - dm
1723 zout(2, j, nout7) = cm + bm
1724 r = r1 - r5
1725 s = s3 - s7
1726 ap = r + s
1727 am = r - s
1728 r = s1 - s5
1729 s = r3 - r7
1730 bp = r + s
1731 bm = r - s
1732 r = s4 - s8
1733 s = r2 - r6
1734 cp = r + s
1735 cm = r - s
1736 r = s2 - s6
1737 s = r4 - r8
1738 dbl = r + s
1739 dm = r - s
1740 r = (cp + dm)*rt2i
1741 s = (-cp + dm)*rt2i
1742 cp = (cm + dbl)*rt2i
1743 dbl = (cm - dbl)*rt2i
1744 zout(1, j, nout2) = ap + r
1745 zout(2, j, nout2) = bm + s
1746 zout(1, j, nout6) = ap - r
1747 zout(2, j, nout6) = bm - s
1748 zout(1, j, nout4) = am + cp
1749 zout(2, j, nout4) = bp + dbl
1750 zout(1, j, nout8) = am - cp
1751 zout(2, j, nout8) = bp - dbl
1752 END DO
1753 END DO
1754 ELSE
1755 ia = 1
1756 nin1 = ia - after
1757 nout1 = ia - atn
1758 DO ib = 1, before
1759 nin1 = nin1 + after
1760 nin2 = nin1 + atb
1761 nin3 = nin2 + atb
1762 nin4 = nin3 + atb
1763 nin5 = nin4 + atb
1764 nin6 = nin5 + atb
1765 nin7 = nin6 + atb
1766 nin8 = nin7 + atb
1767 nout1 = nout1 + atn
1768 nout2 = nout1 + after
1769 nout3 = nout2 + after
1770 nout4 = nout3 + after
1771 nout5 = nout4 + after
1772 nout6 = nout5 + after
1773 nout7 = nout6 + after
1774 nout8 = nout7 + after
1775 DO j = 1, nfft
1776 r1 = zin(1, nin1, j)
1777 s1 = zin(2, nin1, j)
1778 r2 = zin(1, nin2, j)
1779 s2 = zin(2, nin2, j)
1780 r3 = zin(1, nin3, j)
1781 s3 = zin(2, nin3, j)
1782 r4 = zin(1, nin4, j)
1783 s4 = zin(2, nin4, j)
1784 r5 = zin(1, nin5, j)
1785 s5 = zin(2, nin5, j)
1786 r6 = zin(1, nin6, j)
1787 s6 = zin(2, nin6, j)
1788 r7 = zin(1, nin7, j)
1789 s7 = zin(2, nin7, j)
1790 r8 = zin(1, nin8, j)
1791 s8 = zin(2, nin8, j)
1792 r = r1 + r5
1793 s = r3 + r7
1794 ap = r + s
1795 am = r - s
1796 r = r2 + r6
1797 s = r4 + r8
1798 bp = r + s
1799 bm = r - s
1800 r = s1 + s5
1801 s = s3 + s7
1802 cp = r + s
1803 cm = r - s
1804 r = s2 + s6
1805 s = s4 + s8
1806 dbl = r + s
1807 dm = r - s
1808 zout(1, j, nout1) = ap + bp
1809 zout(2, j, nout1) = cp + dbl
1810 zout(1, j, nout5) = ap - bp
1811 zout(2, j, nout5) = cp - dbl
1812 zout(1, j, nout3) = am - dm
1813 zout(2, j, nout3) = cm + bm
1814 zout(1, j, nout7) = am + dm
1815 zout(2, j, nout7) = cm - bm
1816 r = r1 - r5
1817 s = -s3 + s7
1818 ap = r + s
1819 am = r - s
1820 r = s1 - s5
1821 s = r7 - r3
1822 bp = r + s
1823 bm = r - s
1824 r = -s4 + s8
1825 s = r2 - r6
1826 cp = r + s
1827 cm = r - s
1828 r = -s2 + s6
1829 s = r4 - r8
1830 dbl = r + s
1831 dm = r - s
1832 r = (cp + dm)*rt2i
1833 s = (cp - dm)*rt2i
1834 cp = (cm + dbl)*rt2i
1835 dbl = (-cm + dbl)*rt2i
1836 zout(1, j, nout2) = ap + r
1837 zout(2, j, nout2) = bm + s
1838 zout(1, j, nout6) = ap - r
1839 zout(2, j, nout6) = bm - s
1840 zout(1, j, nout4) = am + cp
1841 zout(2, j, nout4) = bp + dbl
1842 zout(1, j, nout8) = am - cp
1843 zout(2, j, nout8) = bp - dbl
1844 END DO
1845 END DO
1846 END IF
1847 ELSE IF (now == 3) THEN
1848 ia = 1
1849 nin1 = ia - after
1850 nout1 = ia - atn
1851 bbs = isign*bb
1852 DO ib = 1, before
1853 nin1 = nin1 + after
1854 nin2 = nin1 + atb
1855 nin3 = nin2 + atb
1856 nout1 = nout1 + atn
1857 nout2 = nout1 + after
1858 nout3 = nout2 + after
1859 DO j = 1, nfft
1860 r1 = zin(1, nin1, j)
1861 s1 = zin(2, nin1, j)
1862 r2 = zin(1, nin2, j)
1863 s2 = zin(2, nin2, j)
1864 r3 = zin(1, nin3, j)
1865 s3 = zin(2, nin3, j)
1866 r = r2 + r3
1867 s = s2 + s3
1868 zout(1, j, nout1) = r + r1
1869 zout(2, j, nout1) = s + s1
1870 r1 = r1 - 0.5_dp*r
1871 s1 = s1 - 0.5_dp*s
1872 r2 = bbs*(r2 - r3)
1873 s2 = bbs*(s2 - s3)
1874 zout(1, j, nout2) = r1 - s2
1875 zout(2, j, nout2) = s1 + r2
1876 zout(1, j, nout3) = r1 + s2
1877 zout(2, j, nout3) = s1 - r2
1878 END DO
1879 END DO
1880 DO ia = 2, after
1881 ias = ia - 1
1882 IF (4*ias == 3*after) THEN
1883 IF (isign == 1) THEN
1884 nin1 = ia - after
1885 nout1 = ia - atn
1886 DO ib = 1, before
1887 nin1 = nin1 + after
1888 nin2 = nin1 + atb
1889 nin3 = nin2 + atb
1890 nout1 = nout1 + atn
1891 nout2 = nout1 + after
1892 nout3 = nout2 + after
1893 DO j = 1, nfft
1894 r1 = zin(1, nin1, j)
1895 s1 = zin(2, nin1, j)
1896 r2 = -zin(2, nin2, j)
1897 s2 = zin(1, nin2, j)
1898 r3 = -zin(1, nin3, j)
1899 s3 = -zin(2, nin3, j)
1900 r = r2 + r3
1901 s = s2 + s3
1902 zout(1, j, nout1) = r + r1
1903 zout(2, j, nout1) = s + s1
1904 r1 = r1 - 0.5_dp*r
1905 s1 = s1 - 0.5_dp*s
1906 r2 = bbs*(r2 - r3)
1907 s2 = bbs*(s2 - s3)
1908 zout(1, j, nout2) = r1 - s2
1909 zout(2, j, nout2) = s1 + r2
1910 zout(1, j, nout3) = r1 + s2
1911 zout(2, j, nout3) = s1 - r2
1912 END DO
1913 END DO
1914 ELSE
1915 nin1 = ia - after
1916 nout1 = ia - atn
1917 DO ib = 1, before
1918 nin1 = nin1 + after
1919 nin2 = nin1 + atb
1920 nin3 = nin2 + atb
1921 nout1 = nout1 + atn
1922 nout2 = nout1 + after
1923 nout3 = nout2 + after
1924 DO j = 1, nfft
1925 r1 = zin(1, nin1, j)
1926 s1 = zin(2, nin1, j)
1927 r2 = zin(2, nin2, j)
1928 s2 = -zin(1, nin2, j)
1929 r3 = -zin(1, nin3, j)
1930 s3 = -zin(2, nin3, j)
1931 r = r2 + r3
1932 s = s2 + s3
1933 zout(1, j, nout1) = r + r1
1934 zout(2, j, nout1) = s + s1
1935 r1 = r1 - 0.5_dp*r
1936 s1 = s1 - 0.5_dp*s
1937 r2 = bbs*(r2 - r3)
1938 s2 = bbs*(s2 - s3)
1939 zout(1, j, nout2) = r1 - s2
1940 zout(2, j, nout2) = s1 + r2
1941 zout(1, j, nout3) = r1 + s2
1942 zout(2, j, nout3) = s1 - r2
1943 END DO
1944 END DO
1945 END IF
1946 ELSE IF (8*ias == 3*after) THEN
1947 IF (isign == 1) THEN
1948 nin1 = ia - after
1949 nout1 = ia - atn
1950 DO ib = 1, before
1951 nin1 = nin1 + after
1952 nin2 = nin1 + atb
1953 nin3 = nin2 + atb
1954 nout1 = nout1 + atn
1955 nout2 = nout1 + after
1956 nout3 = nout2 + after
1957 DO j = 1, nfft
1958 r1 = zin(1, nin1, j)
1959 s1 = zin(2, nin1, j)
1960 r = zin(1, nin2, j)
1961 s = zin(2, nin2, j)
1962 r2 = (r - s)*rt2i
1963 s2 = (r + s)*rt2i
1964 r3 = -zin(2, nin3, j)
1965 s3 = zin(1, nin3, j)
1966 r = r2 + r3
1967 s = s2 + s3
1968 zout(1, j, nout1) = r + r1
1969 zout(2, j, nout1) = s + s1
1970 r1 = r1 - 0.5_dp*r
1971 s1 = s1 - 0.5_dp*s
1972 r2 = bbs*(r2 - r3)
1973 s2 = bbs*(s2 - s3)
1974 zout(1, j, nout2) = r1 - s2
1975 zout(2, j, nout2) = s1 + r2
1976 zout(1, j, nout3) = r1 + s2
1977 zout(2, j, nout3) = s1 - r2
1978 END DO
1979 END DO
1980 ELSE
1981 nin1 = ia - after
1982 nout1 = ia - atn
1983 DO ib = 1, before
1984 nin1 = nin1 + after
1985 nin2 = nin1 + atb
1986 nin3 = nin2 + atb
1987 nout1 = nout1 + atn
1988 nout2 = nout1 + after
1989 nout3 = nout2 + after
1990 DO j = 1, nfft
1991 r1 = zin(1, nin1, j)
1992 s1 = zin(2, nin1, j)
1993 r = zin(1, nin2, j)
1994 s = zin(2, nin2, j)
1995 r2 = (r + s)*rt2i
1996 s2 = (-r + s)*rt2i
1997 r3 = zin(2, nin3, j)
1998 s3 = -zin(1, nin3, j)
1999 r = r2 + r3
2000 s = s2 + s3
2001 zout(1, j, nout1) = r + r1
2002 zout(2, j, nout1) = s + s1
2003 r1 = r1 - 0.5_dp*r
2004 s1 = s1 - 0.5_dp*s
2005 r2 = bbs*(r2 - r3)
2006 s2 = bbs*(s2 - s3)
2007 zout(1, j, nout2) = r1 - s2
2008 zout(2, j, nout2) = s1 + r2
2009 zout(1, j, nout3) = r1 + s2
2010 zout(2, j, nout3) = s1 - r2
2011 END DO
2012 END DO
2013 END IF
2014 ELSE
2015 itt = ias*before
2016 itrig = itt + 1
2017 cr2 = trig(1, itrig)
2018 ci2 = trig(2, itrig)
2019 itrig = itrig + itt
2020 cr3 = trig(1, itrig)
2021 ci3 = trig(2, itrig)
2022 nin1 = ia - after
2023 nout1 = ia - atn
2024 DO ib = 1, before
2025 nin1 = nin1 + after
2026 nin2 = nin1 + atb
2027 nin3 = nin2 + atb
2028 nout1 = nout1 + atn
2029 nout2 = nout1 + after
2030 nout3 = nout2 + after
2031 DO j = 1, nfft
2032 r1 = zin(1, nin1, j)
2033 s1 = zin(2, nin1, j)
2034 r = zin(1, nin2, j)
2035 s = zin(2, nin2, j)
2036 r2 = r*cr2 - s*ci2
2037 s2 = r*ci2 + s*cr2
2038 r = zin(1, nin3, j)
2039 s = zin(2, nin3, j)
2040 r3 = r*cr3 - s*ci3
2041 s3 = r*ci3 + s*cr3
2042 r = r2 + r3
2043 s = s2 + s3
2044 zout(1, j, nout1) = r + r1
2045 zout(2, j, nout1) = s + s1
2046 r1 = r1 - 0.5_dp*r
2047 s1 = s1 - 0.5_dp*s
2048 r2 = bbs*(r2 - r3)
2049 s2 = bbs*(s2 - s3)
2050 zout(1, j, nout2) = r1 - s2
2051 zout(2, j, nout2) = s1 + r2
2052 zout(1, j, nout3) = r1 + s2
2053 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 sin2 = isign*sin2p
2060 sin4 = isign*sin4p
2061 ia = 1
2062 nin1 = ia - after
2063 nout1 = ia - atn
2064 DO ib = 1, before
2065 nin1 = nin1 + after
2066 nin2 = nin1 + atb
2067 nin3 = nin2 + atb
2068 nin4 = nin3 + atb
2069 nin5 = nin4 + atb
2070 nout1 = nout1 + atn
2071 nout2 = nout1 + after
2072 nout3 = nout2 + after
2073 nout4 = nout3 + after
2074 nout5 = nout4 + after
2075 DO j = 1, nfft
2076 r1 = zin(1, nin1, j)
2077 s1 = zin(2, nin1, j)
2078 r2 = zin(1, nin2, j)
2079 s2 = zin(2, nin2, j)
2080 r3 = zin(1, nin3, j)
2081 s3 = zin(2, nin3, j)
2082 r4 = zin(1, nin4, j)
2083 s4 = zin(2, nin4, j)
2084 r5 = zin(1, nin5, j)
2085 s5 = zin(2, nin5, j)
2086 r25 = r2 + r5
2087 r34 = r3 + r4
2088 s25 = s2 - s5
2089 s34 = s3 - s4
2090 zout(1, j, nout1) = r1 + r25 + r34
2091 r = cos2*r25 + cos4*r34 + r1
2092 s = sin2*s25 + sin4*s34
2093 zout(1, j, nout2) = r - s
2094 zout(1, j, nout5) = r + s
2095 r = cos4*r25 + cos2*r34 + r1
2096 s = sin4*s25 - sin2*s34
2097 zout(1, j, nout3) = r - s
2098 zout(1, j, nout4) = r + s
2099 r25 = r2 - r5
2100 r34 = r3 - r4
2101 s25 = s2 + s5
2102 s34 = s3 + s4
2103 zout(2, j, nout1) = s1 + s25 + s34
2104 r = cos2*s25 + cos4*s34 + s1
2105 s = sin2*r25 + sin4*r34
2106 zout(2, j, nout2) = r + s
2107 zout(2, j, nout5) = r - s
2108 r = cos4*s25 + cos2*s34 + s1
2109 s = sin4*r25 - sin2*r34
2110 zout(2, j, nout3) = r + s
2111 zout(2, j, nout4) = r - s
2112 END DO
2113 END DO
2114 DO ia = 2, after
2115 ias = ia - 1
2116 IF (8*ias == 5*after) THEN
2117 IF (isign == 1) THEN
2118 nin1 = ia - after
2119 nout1 = ia - atn
2120 DO ib = 1, before
2121 nin1 = nin1 + after
2122 nin2 = nin1 + atb
2123 nin3 = nin2 + atb
2124 nin4 = nin3 + atb
2125 nin5 = nin4 + atb
2126 nout1 = nout1 + atn
2127 nout2 = nout1 + after
2128 nout3 = nout2 + after
2129 nout4 = nout3 + after
2130 nout5 = nout4 + after
2131 DO j = 1, nfft
2132 r1 = zin(1, nin1, j)
2133 s1 = zin(2, nin1, j)
2134 r = zin(1, nin2, j)
2135 s = zin(2, nin2, j)
2136 r2 = (r - s)*rt2i
2137 s2 = (r + s)*rt2i
2138 r3 = -zin(2, nin3, j)
2139 s3 = zin(1, nin3, j)
2140 r = zin(1, nin4, j)
2141 s = zin(2, nin4, j)
2142 r4 = -(r + s)*rt2i
2143 s4 = (r - s)*rt2i
2144 r5 = -zin(1, nin5, j)
2145 s5 = -zin(2, nin5, j)
2146 r25 = r2 + r5
2147 r34 = r3 + r4
2148 s25 = s2 - s5
2149 s34 = s3 - s4
2150 zout(1, j, nout1) = r1 + r25 + r34
2151 r = cos2*r25 + cos4*r34 + r1
2152 s = sin2*s25 + sin4*s34
2153 zout(1, j, nout2) = r - s
2154 zout(1, j, nout5) = r + s
2155 r = cos4*r25 + cos2*r34 + r1
2156 s = sin4*s25 - sin2*s34
2157 zout(1, j, nout3) = r - s
2158 zout(1, j, nout4) = r + s
2159 r25 = r2 - r5
2160 r34 = r3 - r4
2161 s25 = s2 + s5
2162 s34 = s3 + s4
2163 zout(2, j, nout1) = s1 + s25 + s34
2164 r = cos2*s25 + cos4*s34 + s1
2165 s = sin2*r25 + sin4*r34
2166 zout(2, j, nout2) = r + s
2167 zout(2, j, nout5) = r - s
2168 r = cos4*s25 + cos2*s34 + s1
2169 s = sin4*r25 - sin2*r34
2170 zout(2, j, nout3) = r + s
2171 zout(2, j, nout4) = r - s
2172 END DO
2173 END DO
2174 ELSE
2175 nin1 = ia - after
2176 nout1 = ia - atn
2177 DO ib = 1, before
2178 nin1 = nin1 + after
2179 nin2 = nin1 + atb
2180 nin3 = nin2 + atb
2181 nin4 = nin3 + atb
2182 nin5 = nin4 + atb
2183 nout1 = nout1 + atn
2184 nout2 = nout1 + after
2185 nout3 = nout2 + after
2186 nout4 = nout3 + after
2187 nout5 = nout4 + after
2188 DO j = 1, nfft
2189 r1 = zin(1, nin1, j)
2190 s1 = zin(2, nin1, j)
2191 r = zin(1, nin2, j)
2192 s = zin(2, nin2, j)
2193 r2 = (r + s)*rt2i
2194 s2 = (-r + s)*rt2i
2195 r3 = zin(2, nin3, j)
2196 s3 = -zin(1, nin3, j)
2197 r = zin(1, nin4, j)
2198 s = zin(2, nin4, j)
2199 r4 = (s - r)*rt2i
2200 s4 = -(r + s)*rt2i
2201 r5 = -zin(1, nin5, j)
2202 s5 = -zin(2, nin5, j)
2203 r25 = r2 + r5
2204 r34 = r3 + r4
2205 s25 = s2 - s5
2206 s34 = s3 - s4
2207 zout(1, j, nout1) = r1 + r25 + r34
2208 r = cos2*r25 + cos4*r34 + r1
2209 s = sin2*s25 + sin4*s34
2210 zout(1, j, nout2) = r - s
2211 zout(1, j, nout5) = r + s
2212 r = cos4*r25 + cos2*r34 + r1
2213 s = sin4*s25 - sin2*s34
2214 zout(1, j, nout3) = r - s
2215 zout(1, j, nout4) = r + s
2216 r25 = r2 - r5
2217 r34 = r3 - r4
2218 s25 = s2 + s5
2219 s34 = s3 + s4
2220 zout(2, j, nout1) = s1 + s25 + s34
2221 r = cos2*s25 + cos4*s34 + s1
2222 s = sin2*r25 + sin4*r34
2223 zout(2, j, nout2) = r + s
2224 zout(2, j, nout5) = r - s
2225 r = cos4*s25 + cos2*s34 + s1
2226 s = sin4*r25 - sin2*r34
2227 zout(2, j, nout3) = r + s
2228 zout(2, j, nout4) = r - s
2229 END DO
2230 END DO
2231 END IF
2232 ELSE
2233 ias = ia - 1
2234 itt = ias*before
2235 itrig = itt + 1
2236 cr2 = trig(1, itrig)
2237 ci2 = trig(2, itrig)
2238 itrig = itrig + itt
2239 cr3 = trig(1, itrig)
2240 ci3 = trig(2, itrig)
2241 itrig = itrig + itt
2242 cr4 = trig(1, itrig)
2243 ci4 = trig(2, itrig)
2244 itrig = itrig + itt
2245 cr5 = trig(1, itrig)
2246 ci5 = trig(2, itrig)
2247 nin1 = ia - after
2248 nout1 = ia - atn
2249 DO ib = 1, before
2250 nin1 = nin1 + after
2251 nin2 = nin1 + atb
2252 nin3 = nin2 + atb
2253 nin4 = nin3 + atb
2254 nin5 = nin4 + atb
2255 nout1 = nout1 + atn
2256 nout2 = nout1 + after
2257 nout3 = nout2 + after
2258 nout4 = nout3 + after
2259 nout5 = nout4 + after
2260 DO j = 1, nfft
2261 r1 = zin(1, nin1, j)
2262 s1 = zin(2, nin1, j)
2263 r = zin(1, nin2, j)
2264 s = zin(2, nin2, j)
2265 r2 = r*cr2 - s*ci2
2266 s2 = r*ci2 + s*cr2
2267 r = zin(1, nin3, j)
2268 s = zin(2, nin3, j)
2269 r3 = r*cr3 - s*ci3
2270 s3 = r*ci3 + s*cr3
2271 r = zin(1, nin4, j)
2272 s = zin(2, nin4, j)
2273 r4 = r*cr4 - s*ci4
2274 s4 = r*ci4 + s*cr4
2275 r = zin(1, nin5, j)
2276 s = zin(2, nin5, j)
2277 r5 = r*cr5 - s*ci5
2278 s5 = r*ci5 + s*cr5
2279 r25 = r2 + r5
2280 r34 = r3 + r4
2281 s25 = s2 - s5
2282 s34 = s3 - s4
2283 zout(1, j, nout1) = r1 + r25 + r34
2284 r = cos2*r25 + cos4*r34 + r1
2285 s = sin2*s25 + sin4*s34
2286 zout(1, j, nout2) = r - s
2287 zout(1, j, nout5) = r + s
2288 r = cos4*r25 + cos2*r34 + r1
2289 s = sin4*s25 - sin2*s34
2290 zout(1, j, nout3) = r - s
2291 zout(1, j, nout4) = r + s
2292 r25 = r2 - r5
2293 r34 = r3 - r4
2294 s25 = s2 + s5
2295 s34 = s3 + s4
2296 zout(2, j, nout1) = s1 + s25 + s34
2297 r = cos2*s25 + cos4*s34 + s1
2298 s = sin2*r25 + sin4*r34
2299 zout(2, j, nout2) = r + s
2300 zout(2, j, nout5) = r - s
2301 r = cos4*s25 + cos2*s34 + s1
2302 s = sin4*r25 - sin2*r34
2303 zout(2, j, nout3) = r + s
2304 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 bbs = isign*bb
2311 ia = 1
2312 nin1 = ia - after
2313 nout1 = ia - atn
2314 DO ib = 1, before
2315 nin1 = nin1 + after
2316 nin2 = nin1 + atb
2317 nin3 = nin2 + atb
2318 nin4 = nin3 + atb
2319 nin5 = nin4 + atb
2320 nin6 = nin5 + atb
2321 nout1 = nout1 + atn
2322 nout2 = nout1 + after
2323 nout3 = nout2 + after
2324 nout4 = nout3 + after
2325 nout5 = nout4 + after
2326 nout6 = nout5 + after
2327 DO j = 1, nfft
2328 r2 = zin(1, nin3, j)
2329 s2 = zin(2, nin3, j)
2330 r3 = zin(1, nin5, j)
2331 s3 = zin(2, nin5, j)
2332 r = r2 + r3
2333 s = s2 + s3
2334 r1 = zin(1, nin1, j)
2335 s1 = zin(2, nin1, j)
2336 ur1 = r + r1
2337 ui1 = s + s1
2338 r1 = r1 - 0.5_dp*r
2339 s1 = s1 - 0.5_dp*s
2340 r = r2 - r3
2341 s = s2 - s3
2342 ur2 = r1 - s*bbs
2343 ui2 = s1 + r*bbs
2344 ur3 = r1 + s*bbs
2345 ui3 = s1 - r*bbs
2346
2347 r2 = zin(1, nin6, j)
2348 s2 = zin(2, nin6, j)
2349 r3 = zin(1, nin2, j)
2350 s3 = zin(2, nin2, j)
2351 r = r2 + r3
2352 s = s2 + s3
2353 r1 = zin(1, nin4, j)
2354 s1 = zin(2, nin4, j)
2355 vr1 = r + r1
2356 vi1 = s + s1
2357 r1 = r1 - 0.5_dp*r
2358 s1 = s1 - 0.5_dp*s
2359 r = r2 - r3
2360 s = s2 - s3
2361 vr2 = r1 - s*bbs
2362 vi2 = s1 + r*bbs
2363 vr3 = r1 + s*bbs
2364 vi3 = s1 - r*bbs
2365
2366 zout(1, j, nout1) = ur1 + vr1
2367 zout(2, j, nout1) = ui1 + vi1
2368 zout(1, j, nout5) = ur2 + vr2
2369 zout(2, j, nout5) = ui2 + vi2
2370 zout(1, j, nout3) = ur3 + vr3
2371 zout(2, j, nout3) = ui3 + vi3
2372 zout(1, j, nout4) = ur1 - vr1
2373 zout(2, j, nout4) = ui1 - vi1
2374 zout(1, j, nout2) = ur2 - vr2
2375 zout(2, j, nout2) = ui2 - vi2
2376 zout(1, j, nout6) = ur3 - vr3
2377 zout(2, j, nout6) = ui3 - vi3
2378 END DO
2379 END DO
2380 ELSE
2381 cpabort('Error fftpre')
2382 END IF
2383
2384!-----------------------------------------------------------------------------!
2385
2386 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 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 atn = after*now
2445 atb = after*before
2446
2447 IF (now == 4) THEN
2448 IF (isign == 1) THEN
2449 ia = 1
2450 nin1 = ia - after
2451 nout1 = ia - atn
2452 DO ib = 1, before
2453 nin1 = nin1 + after
2454 nin2 = nin1 + atb
2455 nin3 = nin2 + atb
2456 nin4 = nin3 + atb
2457 nout1 = nout1 + atn
2458 nout2 = nout1 + after
2459 nout3 = nout2 + after
2460 nout4 = nout3 + after
2461 DO j = 1, nfft
2462 r1 = zin(1, j, nin1)
2463 s1 = zin(2, j, nin1)
2464 r2 = zin(1, j, nin2)
2465 s2 = zin(2, j, nin2)
2466 r3 = zin(1, j, nin3)
2467 s3 = zin(2, j, nin3)
2468 r4 = zin(1, j, nin4)
2469 s4 = zin(2, j, nin4)
2470 r = r1 + r3
2471 s = r2 + r4
2472 zout(1, j, nout1) = r + s
2473 zout(1, j, nout3) = r - s
2474 r = r1 - r3
2475 s = s2 - s4
2476 zout(1, j, nout2) = r - s
2477 zout(1, j, nout4) = r + s
2478 r = s1 + s3
2479 s = s2 + s4
2480 zout(2, j, nout1) = r + s
2481 zout(2, j, nout3) = r - s
2482 r = s1 - s3
2483 s = r2 - r4
2484 zout(2, j, nout2) = r + s
2485 zout(2, j, nout4) = r - s
2486 END DO
2487 END DO
2488 DO ia = 2, after
2489 ias = ia - 1
2490 IF (2*ias == after) THEN
2491 nin1 = ia - after
2492 nout1 = ia - atn
2493 DO ib = 1, before
2494 nin1 = nin1 + after
2495 nin2 = nin1 + atb
2496 nin3 = nin2 + atb
2497 nin4 = nin3 + atb
2498 nout1 = nout1 + atn
2499 nout2 = nout1 + after
2500 nout3 = nout2 + after
2501 nout4 = nout3 + after
2502 DO j = 1, nfft
2503 r1 = zin(1, j, nin1)
2504 s1 = zin(2, j, nin1)
2505 r = zin(1, j, nin2)
2506 s = zin(2, j, nin2)
2507 r2 = (r - s)*rt2i
2508 s2 = (r + s)*rt2i
2509 r3 = -zin(2, j, nin3)
2510 s3 = zin(1, j, nin3)
2511 r = zin(1, j, nin4)
2512 s = zin(2, j, nin4)
2513 r4 = -(r + s)*rt2i
2514 s4 = (r - s)*rt2i
2515 r = r1 + r3
2516 s = r2 + r4
2517 zout(1, j, nout1) = r + s
2518 zout(1, j, nout3) = r - s
2519 r = r1 - r3
2520 s = s2 - s4
2521 zout(1, j, nout2) = r - s
2522 zout(1, j, nout4) = r + s
2523 r = s1 + s3
2524 s = s2 + s4
2525 zout(2, j, nout1) = r + s
2526 zout(2, j, nout3) = r - s
2527 r = s1 - s3
2528 s = r2 - r4
2529 zout(2, j, nout2) = r + s
2530 zout(2, j, nout4) = r - s
2531 END DO
2532 END DO
2533 ELSE
2534 itt = ias*before
2535 itrig = itt + 1
2536 cr2 = trig(1, itrig)
2537 ci2 = trig(2, itrig)
2538 itrig = itrig + itt
2539 cr3 = trig(1, itrig)
2540 ci3 = trig(2, itrig)
2541 itrig = itrig + itt
2542 cr4 = trig(1, itrig)
2543 ci4 = trig(2, itrig)
2544 nin1 = ia - after
2545 nout1 = ia - atn
2546 DO ib = 1, before
2547 nin1 = nin1 + after
2548 nin2 = nin1 + atb
2549 nin3 = nin2 + atb
2550 nin4 = nin3 + atb
2551 nout1 = nout1 + atn
2552 nout2 = nout1 + after
2553 nout3 = nout2 + after
2554 nout4 = nout3 + after
2555 DO j = 1, nfft
2556 r1 = zin(1, j, nin1)
2557 s1 = zin(2, j, nin1)
2558 r = zin(1, j, nin2)
2559 s = zin(2, j, nin2)
2560 r2 = r*cr2 - s*ci2
2561 s2 = r*ci2 + s*cr2
2562 r = zin(1, j, nin3)
2563 s = zin(2, j, nin3)
2564 r3 = r*cr3 - s*ci3
2565 s3 = r*ci3 + s*cr3
2566 r = zin(1, j, nin4)
2567 s = zin(2, j, nin4)
2568 r4 = r*cr4 - s*ci4
2569 s4 = r*ci4 + s*cr4
2570 r = r1 + r3
2571 s = r2 + r4
2572 zout(1, j, nout1) = r + s
2573 zout(1, j, nout3) = r - s
2574 r = r1 - r3
2575 s = s2 - s4
2576 zout(1, j, nout2) = r - s
2577 zout(1, j, nout4) = r + s
2578 r = s1 + s3
2579 s = s2 + s4
2580 zout(2, j, nout1) = r + s
2581 zout(2, j, nout3) = r - s
2582 r = s1 - s3
2583 s = r2 - r4
2584 zout(2, j, nout2) = r + s
2585 zout(2, j, nout4) = r - s
2586 END DO
2587 END DO
2588 END IF
2589 END DO
2590 ELSE
2591 ia = 1
2592 nin1 = ia - after
2593 nout1 = ia - atn
2594 DO ib = 1, before
2595 nin1 = nin1 + after
2596 nin2 = nin1 + atb
2597 nin3 = nin2 + atb
2598 nin4 = nin3 + atb
2599 nout1 = nout1 + atn
2600 nout2 = nout1 + after
2601 nout3 = nout2 + after
2602 nout4 = nout3 + after
2603 DO j = 1, nfft
2604 r1 = zin(1, j, nin1)
2605 s1 = zin(2, j, nin1)
2606 r2 = zin(1, j, nin2)
2607 s2 = zin(2, j, nin2)
2608 r3 = zin(1, j, nin3)
2609 s3 = zin(2, j, nin3)
2610 r4 = zin(1, j, nin4)
2611 s4 = zin(2, j, nin4)
2612 r = r1 + r3
2613 s = r2 + r4
2614 zout(1, j, nout1) = r + s
2615 zout(1, j, nout3) = r - s
2616 r = r1 - r3
2617 s = s2 - s4
2618 zout(1, j, nout2) = r + s
2619 zout(1, j, nout4) = r - s
2620 r = s1 + s3
2621 s = s2 + s4
2622 zout(2, j, nout1) = r + s
2623 zout(2, j, nout3) = r - s
2624 r = s1 - s3
2625 s = r2 - r4
2626 zout(2, j, nout2) = r - s
2627 zout(2, j, nout4) = r + s
2628 END DO
2629 END DO
2630 DO ia = 2, after
2631 ias = ia - 1
2632 IF (2*ias == after) THEN
2633 nin1 = ia - after
2634 nout1 = ia - atn
2635 DO ib = 1, before
2636 nin1 = nin1 + after
2637 nin2 = nin1 + atb
2638 nin3 = nin2 + atb
2639 nin4 = nin3 + atb
2640 nout1 = nout1 + atn
2641 nout2 = nout1 + after
2642 nout3 = nout2 + after
2643 nout4 = nout3 + after
2644 DO j = 1, nfft
2645 r1 = zin(1, j, nin1)
2646 s1 = zin(2, j, nin1)
2647 r = zin(1, j, nin2)
2648 s = zin(2, j, nin2)
2649 r2 = (r + s)*rt2i
2650 s2 = (s - r)*rt2i
2651 r3 = zin(2, j, nin3)
2652 s3 = -zin(1, j, nin3)
2653 r = zin(1, j, nin4)
2654 s = zin(2, j, nin4)
2655 r4 = (s - r)*rt2i
2656 s4 = -(r + s)*rt2i
2657 r = r1 + r3
2658 s = r2 + r4
2659 zout(1, j, nout1) = r + s
2660 zout(1, j, nout3) = r - s
2661 r = r1 - r3
2662 s = s2 - s4
2663 zout(1, j, nout2) = r + s
2664 zout(1, j, nout4) = r - s
2665 r = s1 + s3
2666 s = s2 + s4
2667 zout(2, j, nout1) = r + s
2668 zout(2, j, nout3) = r - s
2669 r = s1 - s3
2670 s = r2 - r4
2671 zout(2, j, nout2) = r - s
2672 zout(2, j, nout4) = r + s
2673 END DO
2674 END DO
2675 ELSE
2676 itt = ias*before
2677 itrig = itt + 1
2678 cr2 = trig(1, itrig)
2679 ci2 = trig(2, itrig)
2680 itrig = itrig + itt
2681 cr3 = trig(1, itrig)
2682 ci3 = trig(2, itrig)
2683 itrig = itrig + itt
2684 cr4 = trig(1, itrig)
2685 ci4 = trig(2, itrig)
2686 nin1 = ia - after
2687 nout1 = ia - atn
2688 DO ib = 1, before
2689 nin1 = nin1 + after
2690 nin2 = nin1 + atb
2691 nin3 = nin2 + atb
2692 nin4 = nin3 + atb
2693 nout1 = nout1 + atn
2694 nout2 = nout1 + after
2695 nout3 = nout2 + after
2696 nout4 = nout3 + after
2697 DO j = 1, nfft
2698 r1 = zin(1, j, nin1)
2699 s1 = zin(2, j, nin1)
2700 r = zin(1, j, nin2)
2701 s = zin(2, j, nin2)
2702 r2 = r*cr2 - s*ci2
2703 s2 = r*ci2 + s*cr2
2704 r = zin(1, j, nin3)
2705 s = zin(2, j, nin3)
2706 r3 = r*cr3 - s*ci3
2707 s3 = r*ci3 + s*cr3
2708 r = zin(1, j, nin4)
2709 s = zin(2, j, nin4)
2710 r4 = r*cr4 - s*ci4
2711 s4 = r*ci4 + s*cr4
2712 r = r1 + r3
2713 s = r2 + r4
2714 zout(1, j, nout1) = r + s
2715 zout(1, j, nout3) = r - s
2716 r = r1 - r3
2717 s = s2 - s4
2718 zout(1, j, nout2) = r + s
2719 zout(1, j, nout4) = r - s
2720 r = s1 + s3
2721 s = s2 + s4
2722 zout(2, j, nout1) = r + s
2723 zout(2, j, nout3) = r - s
2724 r = s1 - s3
2725 s = r2 - r4
2726 zout(2, j, nout2) = r - s
2727 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 IF (isign == -1) THEN
2735 ia = 1
2736 nin1 = ia - after
2737 nout1 = ia - atn
2738 DO ib = 1, before
2739 nin1 = nin1 + after
2740 nin2 = nin1 + atb
2741 nin3 = nin2 + atb
2742 nin4 = nin3 + atb
2743 nin5 = nin4 + atb
2744 nin6 = nin5 + atb
2745 nin7 = nin6 + atb
2746 nin8 = nin7 + atb
2747 nout1 = nout1 + atn
2748 nout2 = nout1 + after
2749 nout3 = nout2 + after
2750 nout4 = nout3 + after
2751 nout5 = nout4 + after
2752 nout6 = nout5 + after
2753 nout7 = nout6 + after
2754 nout8 = nout7 + after
2755 DO j = 1, nfft
2756 r1 = zin(1, j, nin1)
2757 s1 = zin(2, j, nin1)
2758 r2 = zin(1, j, nin2)
2759 s2 = zin(2, j, nin2)
2760 r3 = zin(1, j, nin3)
2761 s3 = zin(2, j, nin3)
2762 r4 = zin(1, j, nin4)
2763 s4 = zin(2, j, nin4)
2764 r5 = zin(1, j, nin5)
2765 s5 = zin(2, j, nin5)
2766 r6 = zin(1, j, nin6)
2767 s6 = zin(2, j, nin6)
2768 r7 = zin(1, j, nin7)
2769 s7 = zin(2, j, nin7)
2770 r8 = zin(1, j, nin8)
2771 s8 = zin(2, j, nin8)
2772 r = r1 + r5
2773 s = r3 + r7
2774 ap = r + s
2775 am = r - s
2776 r = r2 + r6
2777 s = r4 + r8
2778 bp = r + s
2779 bm = r - s
2780 r = s1 + s5
2781 s = s3 + s7
2782 cp = r + s
2783 cm = r - s
2784 r = s2 + s6
2785 s = s4 + s8
2786 dbl = r + s
2787 dm = r - s
2788 zout(1, j, nout1) = ap + bp
2789 zout(2, j, nout1) = cp + dbl
2790 zout(1, j, nout5) = ap - bp
2791 zout(2, j, nout5) = cp - dbl
2792 zout(1, j, nout3) = am + dm
2793 zout(2, j, nout3) = cm - bm
2794 zout(1, j, nout7) = am - dm
2795 zout(2, j, nout7) = cm + bm
2796 r = r1 - r5
2797 s = s3 - s7
2798 ap = r + s
2799 am = r - s
2800 r = s1 - s5
2801 s = r3 - r7
2802 bp = r + s
2803 bm = r - s
2804 r = s4 - s8
2805 s = r2 - r6
2806 cp = r + s
2807 cm = r - s
2808 r = s2 - s6
2809 s = r4 - r8
2810 dbl = r + s
2811 dm = r - s
2812 r = (cp + dm)*rt2i
2813 s = (-cp + dm)*rt2i
2814 cp = (cm + dbl)*rt2i
2815 dbl = (cm - dbl)*rt2i
2816 zout(1, j, nout2) = ap + r
2817 zout(2, j, nout2) = bm + s
2818 zout(1, j, nout6) = ap - r
2819 zout(2, j, nout6) = bm - s
2820 zout(1, j, nout4) = am + cp
2821 zout(2, j, nout4) = bp + dbl
2822 zout(1, j, nout8) = am - cp
2823 zout(2, j, nout8) = bp - dbl
2824 END DO
2825 END DO
2826 ELSE
2827 ia = 1
2828 nin1 = ia - after
2829 nout1 = ia - atn
2830 DO ib = 1, before
2831 nin1 = nin1 + after
2832 nin2 = nin1 + atb
2833 nin3 = nin2 + atb
2834 nin4 = nin3 + atb
2835 nin5 = nin4 + atb
2836 nin6 = nin5 + atb
2837 nin7 = nin6 + atb
2838 nin8 = nin7 + atb
2839 nout1 = nout1 + atn
2840 nout2 = nout1 + after
2841 nout3 = nout2 + after
2842 nout4 = nout3 + after
2843 nout5 = nout4 + after
2844 nout6 = nout5 + after
2845 nout7 = nout6 + after
2846 nout8 = nout7 + after
2847 DO j = 1, nfft
2848 r1 = zin(1, j, nin1)
2849 s1 = zin(2, j, nin1)
2850 r2 = zin(1, j, nin2)
2851 s2 = zin(2, j, nin2)
2852 r3 = zin(1, j, nin3)
2853 s3 = zin(2, j, nin3)
2854 r4 = zin(1, j, nin4)
2855 s4 = zin(2, j, nin4)
2856 r5 = zin(1, j, nin5)
2857 s5 = zin(2, j, nin5)
2858 r6 = zin(1, j, nin6)
2859 s6 = zin(2, j, nin6)
2860 r7 = zin(1, j, nin7)
2861 s7 = zin(2, j, nin7)
2862 r8 = zin(1, j, nin8)
2863 s8 = zin(2, j, nin8)
2864 r = r1 + r5
2865 s = r3 + r7
2866 ap = r + s
2867 am = r - s
2868 r = r2 + r6
2869 s = r4 + r8
2870 bp = r + s
2871 bm = r - s
2872 r = s1 + s5
2873 s = s3 + s7
2874 cp = r + s
2875 cm = r - s
2876 r = s2 + s6
2877 s = s4 + s8
2878 dbl = r + s
2879 dm = r - s
2880 zout(1, j, nout1) = ap + bp
2881 zout(2, j, nout1) = cp + dbl
2882 zout(1, j, nout5) = ap - bp
2883 zout(2, j, nout5) = cp - dbl
2884 zout(1, j, nout3) = am - dm
2885 zout(2, j, nout3) = cm + bm
2886 zout(1, j, nout7) = am + dm
2887 zout(2, j, nout7) = cm - bm
2888 r = r1 - r5
2889 s = -s3 + s7
2890 ap = r + s
2891 am = r - s
2892 r = s1 - s5
2893 s = r7 - r3
2894 bp = r + s
2895 bm = r - s
2896 r = -s4 + s8
2897 s = r2 - r6
2898 cp = r + s
2899 cm = r - s
2900 r = -s2 + s6
2901 s = r4 - r8
2902 dbl = r + s
2903 dm = r - s
2904 r = (cp + dm)*rt2i
2905 s = (cp - dm)*rt2i
2906 cp = (cm + dbl)*rt2i
2907 dbl = (-cm + dbl)*rt2i
2908 zout(1, j, nout2) = ap + r
2909 zout(2, j, nout2) = bm + s
2910 zout(1, j, nout6) = ap - r
2911 zout(2, j, nout6) = bm - s
2912 zout(1, j, nout4) = am + cp
2913 zout(2, j, nout4) = bp + dbl
2914 zout(1, j, nout8) = am - cp
2915 zout(2, j, nout8) = bp - dbl
2916 END DO
2917 END DO
2918 END IF
2919 ELSE IF (now == 3) THEN
2920 bbs = isign*bb
2921 ia = 1
2922 nin1 = ia - after
2923 nout1 = ia - atn
2924 DO ib = 1, before
2925 nin1 = nin1 + after
2926 nin2 = nin1 + atb
2927 nin3 = nin2 + atb
2928 nout1 = nout1 + atn
2929 nout2 = nout1 + after
2930 nout3 = nout2 + after
2931 DO j = 1, nfft
2932 r1 = zin(1, j, nin1)
2933 s1 = zin(2, j, nin1)
2934 r2 = zin(1, j, nin2)
2935 s2 = zin(2, j, nin2)
2936 r3 = zin(1, j, nin3)
2937 s3 = zin(2, j, nin3)
2938 r = r2 + r3
2939 s = s2 + s3
2940 zout(1, j, nout1) = r + r1
2941 zout(2, j, nout1) = s + s1
2942 r1 = r1 - 0.5_dp*r
2943 s1 = s1 - 0.5_dp*s
2944 r2 = bbs*(r2 - r3)
2945 s2 = bbs*(s2 - s3)
2946 zout(1, j, nout2) = r1 - s2
2947 zout(2, j, nout2) = s1 + r2
2948 zout(1, j, nout3) = r1 + s2
2949 zout(2, j, nout3) = s1 - r2
2950 END DO
2951 END DO
2952 DO ia = 2, after
2953 ias = ia - 1
2954 IF (4*ias == 3*after) THEN
2955 IF (isign == 1) THEN
2956 nin1 = ia - after
2957 nout1 = ia - atn
2958 DO ib = 1, before
2959 nin1 = nin1 + after
2960 nin2 = nin1 + atb
2961 nin3 = nin2 + atb
2962 nout1 = nout1 + atn
2963 nout2 = nout1 + after
2964 nout3 = nout2 + after
2965 DO j = 1, nfft
2966 r1 = zin(1, j, nin1)
2967 s1 = zin(2, j, nin1)
2968 r2 = -zin(2, j, nin2)
2969 s2 = zin(1, j, nin2)
2970 r3 = -zin(1, j, nin3)
2971 s3 = -zin(2, j, nin3)
2972 r = r2 + r3
2973 s = s2 + s3
2974 zout(1, j, nout1) = r + r1
2975 zout(2, j, nout1) = s + s1
2976 r1 = r1 - 0.5_dp*r
2977 s1 = s1 - 0.5_dp*s
2978 r2 = bbs*(r2 - r3)
2979 s2 = bbs*(s2 - s3)
2980 zout(1, j, nout2) = r1 - s2
2981 zout(2, j, nout2) = s1 + r2
2982 zout(1, j, nout3) = r1 + s2
2983 zout(2, j, nout3) = s1 - r2
2984 END DO
2985 END DO
2986 ELSE
2987 nin1 = ia - after
2988 nout1 = ia - atn
2989 DO ib = 1, before
2990 nin1 = nin1 + after
2991 nin2 = nin1 + atb
2992 nin3 = nin2 + atb
2993 nout1 = nout1 + atn
2994 nout2 = nout1 + after
2995 nout3 = nout2 + after
2996 DO j = 1, nfft
2997 r1 = zin(1, j, nin1)
2998 s1 = zin(2, j, nin1)
2999 r2 = zin(2, j, nin2)
3000 s2 = -zin(1, j, nin2)
3001 r3 = -zin(1, j, nin3)
3002 s3 = -zin(2, j, nin3)
3003 r = r2 + r3
3004 s = s2 + s3
3005 zout(1, j, nout1) = r + r1
3006 zout(2, j, nout1) = s + s1
3007 r1 = r1 - 0.5_dp*r
3008 s1 = s1 - 0.5_dp*s
3009 r2 = bbs*(r2 - r3)
3010 s2 = bbs*(s2 - s3)
3011 zout(1, j, nout2) = r1 - s2
3012 zout(2, j, nout2) = s1 + r2
3013 zout(1, j, nout3) = r1 + s2
3014 zout(2, j, nout3) = s1 - r2
3015 END DO
3016 END DO
3017 END IF
3018 ELSE IF (8*ias == 3*after) THEN
3019 IF (isign == 1) THEN
3020 nin1 = ia - after
3021 nout1 = ia - atn
3022 DO ib = 1, before
3023 nin1 = nin1 + after
3024 nin2 = nin1 + atb
3025 nin3 = nin2 + atb
3026 nout1 = nout1 + atn
3027 nout2 = nout1 + after
3028 nout3 = nout2 + after
3029 DO j = 1, nfft
3030 r1 = zin(1, j, nin1)
3031 s1 = zin(2, j, nin1)
3032 r = zin(1, j, nin2)
3033 s = zin(2, j, nin2)
3034 r2 = (r - s)*rt2i
3035 s2 = (r + s)*rt2i
3036 r3 = -zin(2, j, nin3)
3037 s3 = zin(1, j, nin3)
3038 r = r2 + r3
3039 s = s2 + s3
3040 zout(1, j, nout1) = r + r1
3041 zout(2, j, nout1) = s + s1
3042 r1 = r1 - 0.5_dp*r
3043 s1 = s1 - 0.5_dp*s
3044 r2 = bbs*(r2 - r3)
3045 s2 = bbs*(s2 - s3)
3046 zout(1, j, nout2) = r1 - s2
3047 zout(2, j, nout2) = s1 + r2
3048 zout(1, j, nout3) = r1 + s2
3049 zout(2, j, nout3) = s1 - r2
3050 END DO
3051 END DO
3052 ELSE
3053 nin1 = ia - after
3054 nout1 = ia - atn
3055 DO ib = 1, before
3056 nin1 = nin1 + after
3057 nin2 = nin1 + atb
3058 nin3 = nin2 + atb
3059 nout1 = nout1 + atn
3060 nout2 = nout1 + after
3061 nout3 = nout2 + after
3062 DO j = 1, nfft
3063 r1 = zin(1, j, nin1)
3064 s1 = zin(2, j, nin1)
3065 r = zin(1, j, nin2)
3066 s = zin(2, j, nin2)
3067 r2 = (r + s)*rt2i
3068 s2 = (-r + s)*rt2i
3069 r3 = zin(2, j, nin3)
3070 s3 = -zin(1, j, nin3)
3071 r = r2 + r3
3072 s = s2 + s3
3073 zout(1, j, nout1) = r + r1
3074 zout(2, j, nout1) = s + s1
3075 r1 = r1 - 0.5_dp*r
3076 s1 = s1 - 0.5_dp*s
3077 r2 = bbs*(r2 - r3)
3078 s2 = bbs*(s2 - s3)
3079 zout(1, j, nout2) = r1 - s2
3080 zout(2, j, nout2) = s1 + r2
3081 zout(1, j, nout3) = r1 + s2
3082 zout(2, j, nout3) = s1 - r2
3083 END DO
3084 END DO
3085 END IF
3086 ELSE
3087 itt = ias*before
3088 itrig = itt + 1
3089 cr2 = trig(1, itrig)
3090 ci2 = trig(2, itrig)
3091 itrig = itrig + itt
3092 cr3 = trig(1, itrig)
3093 ci3 = trig(2, itrig)
3094 nin1 = ia - after
3095 nout1 = ia - atn
3096 DO ib = 1, before
3097 nin1 = nin1 + after
3098 nin2 = nin1 + atb
3099 nin3 = nin2 + atb
3100 nout1 = nout1 + atn
3101 nout2 = nout1 + after
3102 nout3 = nout2 + after
3103 DO j = 1, nfft
3104 r1 = zin(1, j, nin1)
3105 s1 = zin(2, j, nin1)
3106 r = zin(1, j, nin2)
3107 s = zin(2, j, nin2)
3108 r2 = r*cr2 - s*ci2
3109 s2 = r*ci2 + s*cr2
3110 r = zin(1, j, nin3)
3111 s = zin(2, j, nin3)
3112 r3 = r*cr3 - s*ci3
3113 s3 = r*ci3 + s*cr3
3114 r = r2 + r3
3115 s = s2 + s3
3116 zout(1, j, nout1) = r + r1
3117 zout(2, j, nout1) = s + s1
3118 r1 = r1 - 0.5_dp*r
3119 s1 = s1 - 0.5_dp*s
3120 r2 = bbs*(r2 - r3)
3121 s2 = bbs*(s2 - s3)
3122 zout(1, j, nout2) = r1 - s2
3123 zout(2, j, nout2) = s1 + r2
3124 zout(1, j, nout3) = r1 + s2
3125 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 sin2 = isign*sin2p
3132 sin4 = isign*sin4p
3133 ia = 1
3134 nin1 = ia - after
3135 nout1 = ia - atn
3136 DO ib = 1, before
3137 nin1 = nin1 + after
3138 nin2 = nin1 + atb
3139 nin3 = nin2 + atb
3140 nin4 = nin3 + atb
3141 nin5 = nin4 + atb
3142 nout1 = nout1 + atn
3143 nout2 = nout1 + after
3144 nout3 = nout2 + after
3145 nout4 = nout3 + after
3146 nout5 = nout4 + after
3147 DO j = 1, nfft
3148 r1 = zin(1, j, nin1)
3149 s1 = zin(2, j, nin1)
3150 r2 = zin(1, j, nin2)
3151 s2 = zin(2, j, nin2)
3152 r3 = zin(1, j, nin3)
3153 s3 = zin(2, j, nin3)
3154 r4 = zin(1, j, nin4)
3155 s4 = zin(2, j, nin4)
3156 r5 = zin(1, j, nin5)
3157 s5 = zin(2, j, nin5)
3158 r25 = r2 + r5
3159 r34 = r3 + r4
3160 s25 = s2 - s5
3161 s34 = s3 - s4
3162 zout(1, j, nout1) = r1 + r25 + r34
3163 r = cos2*r25 + cos4*r34 + r1
3164 s = sin2*s25 + sin4*s34
3165 zout(1, j, nout2) = r - s
3166 zout(1, j, nout5) = r + s
3167 r = cos4*r25 + cos2*r34 + r1
3168 s = sin4*s25 - sin2*s34
3169 zout(1, j, nout3) = r - s
3170 zout(1, j, nout4) = r + s
3171 r25 = r2 - r5
3172 r34 = r3 - r4
3173 s25 = s2 + s5
3174 s34 = s3 + s4
3175 zout(2, j, nout1) = s1 + s25 + s34
3176 r = cos2*s25 + cos4*s34 + s1
3177 s = sin2*r25 + sin4*r34
3178 zout(2, j, nout2) = r + s
3179 zout(2, j, nout5) = r - s
3180 r = cos4*s25 + cos2*s34 + s1
3181 s = sin4*r25 - sin2*r34
3182 zout(2, j, nout3) = r + s
3183 zout(2, j, nout4) = r - s
3184 END DO
3185 END DO
3186 DO ia = 2, after
3187 ias = ia - 1
3188 IF (8*ias == 5*after) THEN
3189 IF (isign == 1) THEN
3190 nin1 = ia - after
3191 nout1 = ia - atn
3192 DO ib = 1, before
3193 nin1 = nin1 + after
3194 nin2 = nin1 + atb
3195 nin3 = nin2 + atb
3196 nin4 = nin3 + atb
3197 nin5 = nin4 + atb
3198 nout1 = nout1 + atn
3199 nout2 = nout1 + after
3200 nout3 = nout2 + after
3201 nout4 = nout3 + after
3202 nout5 = nout4 + after
3203 DO j = 1, nfft
3204 r1 = zin(1, j, nin1)
3205 s1 = zin(2, j, nin1)
3206 r = zin(1, j, nin2)
3207 s = zin(2, j, nin2)
3208 r2 = (r - s)*rt2i
3209 s2 = (r + s)*rt2i
3210 r3 = -zin(2, j, nin3)
3211 s3 = zin(1, j, nin3)
3212 r = zin(1, j, nin4)
3213 s = zin(2, j, nin4)
3214 r4 = -(r + s)*rt2i
3215 s4 = (r - s)*rt2i
3216 r5 = -zin(1, j, nin5)
3217 s5 = -zin(2, j, nin5)
3218 r25 = r2 + r5
3219 r34 = r3 + r4
3220 s25 = s2 - s5
3221 s34 = s3 - s4
3222 zout(1, j, nout1) = r1 + r25 + r34
3223 r = cos2*r25 + cos4*r34 + r1
3224 s = sin2*s25 + sin4*s34
3225 zout(1, j, nout2) = r - s
3226 zout(1, j, nout5) = r + s
3227 r = cos4*r25 + cos2*r34 + r1
3228 s = sin4*s25 - sin2*s34
3229 zout(1, j, nout3) = r - s
3230 zout(1, j, nout4) = r + s
3231 r25 = r2 - r5
3232 r34 = r3 - r4
3233 s25 = s2 + s5
3234 s34 = s3 + s4
3235 zout(2, j, nout1) = s1 + s25 + s34
3236 r = cos2*s25 + cos4*s34 + s1
3237 s = sin2*r25 + sin4*r34
3238 zout(2, j, nout2) = r + s
3239 zout(2, j, nout5) = r - s
3240 r = cos4*s25 + cos2*s34 + s1
3241 s = sin4*r25 - sin2*r34
3242 zout(2, j, nout3) = r + s
3243 zout(2, j, nout4) = r - s
3244 END DO
3245 END DO
3246 ELSE
3247 nin1 = ia - after
3248 nout1 = ia - atn
3249 DO ib = 1, before
3250 nin1 = nin1 + after
3251 nin2 = nin1 + atb
3252 nin3 = nin2 + atb
3253 nin4 = nin3 + atb
3254 nin5 = nin4 + atb
3255 nout1 = nout1 + atn
3256 nout2 = nout1 + after
3257 nout3 = nout2 + after
3258 nout4 = nout3 + after
3259 nout5 = nout4 + after
3260 DO j = 1, nfft
3261 r1 = zin(1, j, nin1)
3262 s1 = zin(2, j, nin1)
3263 r = zin(1, j, nin2)
3264 s = zin(2, j, nin2)
3265 r2 = (r + s)*rt2i
3266 s2 = (-r + s)*rt2i
3267 r3 = zin(2, j, nin3)
3268 s3 = -zin(1, j, nin3)
3269 r = zin(1, j, nin4)
3270 s = zin(2, j, nin4)
3271 r4 = (s - r)*rt2i
3272 s4 = -(r + s)*rt2i
3273 r5 = -zin(1, j, nin5)
3274 s5 = -zin(2, j, nin5)
3275 r25 = r2 + r5
3276 r34 = r3 + r4
3277 s25 = s2 - s5
3278 s34 = s3 - s4
3279 zout(1, j, nout1) = r1 + r25 + r34
3280 r = cos2*r25 + cos4*r34 + r1
3281 s = sin2*s25 + sin4*s34
3282 zout(1, j, nout2) = r - s
3283 zout(1, j, nout5) = r + s
3284 r = cos4*r25 + cos2*r34 + r1
3285 s = sin4*s25 - sin2*s34
3286 zout(1, j, nout3) = r - s
3287 zout(1, j, nout4) = r + s
3288 r25 = r2 - r5
3289 r34 = r3 - r4
3290 s25 = s2 + s5
3291 s34 = s3 + s4
3292 zout(2, j, nout1) = s1 + s25 + s34
3293 r = cos2*s25 + cos4*s34 + s1
3294 s = sin2*r25 + sin4*r34
3295 zout(2, j, nout2) = r + s
3296 zout(2, j, nout5) = r - s
3297 r = cos4*s25 + cos2*s34 + s1
3298 s = sin4*r25 - sin2*r34
3299 zout(2, j, nout3) = r + s
3300 zout(2, j, nout4) = r - s
3301 END DO
3302 END DO
3303 END IF
3304 ELSE
3305 ias = ia - 1
3306 itt = ias*before
3307 itrig = itt + 1
3308 cr2 = trig(1, itrig)
3309 ci2 = trig(2, itrig)
3310 itrig = itrig + itt
3311 cr3 = trig(1, itrig)
3312 ci3 = trig(2, itrig)
3313 itrig = itrig + itt
3314 cr4 = trig(1, itrig)
3315 ci4 = trig(2, itrig)
3316 itrig = itrig + itt
3317 cr5 = trig(1, itrig)
3318 ci5 = trig(2, itrig)
3319 nin1 = ia - after
3320 nout1 = ia - atn
3321 DO ib = 1, before
3322 nin1 = nin1 + after
3323 nin2 = nin1 + atb
3324 nin3 = nin2 + atb
3325 nin4 = nin3 + atb
3326 nin5 = nin4 + atb
3327 nout1 = nout1 + atn
3328 nout2 = nout1 + after
3329 nout3 = nout2 + after
3330 nout4 = nout3 + after
3331 nout5 = nout4 + after
3332 DO j = 1, nfft
3333 r1 = zin(1, j, nin1)
3334 s1 = zin(2, j, nin1)
3335 r = zin(1, j, nin2)
3336 s = zin(2, j, nin2)
3337 r2 = r*cr2 - s*ci2
3338 s2 = r*ci2 + s*cr2
3339 r = zin(1, j, nin3)
3340 s = zin(2, j, nin3)
3341 r3 = r*cr3 - s*ci3
3342 s3 = r*ci3 + s*cr3
3343 r = zin(1, j, nin4)
3344 s = zin(2, j, nin4)
3345 r4 = r*cr4 - s*ci4
3346 s4 = r*ci4 + s*cr4
3347 r = zin(1, j, nin5)
3348 s = zin(2, j, nin5)
3349 r5 = r*cr5 - s*ci5
3350 s5 = r*ci5 + s*cr5
3351 r25 = r2 + r5
3352 r34 = r3 + r4
3353 s25 = s2 - s5
3354 s34 = s3 - s4
3355 zout(1, j, nout1) = r1 + r25 + r34
3356 r = cos2*r25 + cos4*r34 + r1
3357 s = sin2*s25 + sin4*s34
3358 zout(1, j, nout2) = r - s
3359 zout(1, j, nout5) = r + s
3360 r = cos4*r25 + cos2*r34 + r1
3361 s = sin4*s25 - sin2*s34
3362 zout(1, j, nout3) = r - s
3363 zout(1, j, nout4) = r + s
3364 r25 = r2 - r5
3365 r34 = r3 - r4
3366 s25 = s2 + s5
3367 s34 = s3 + s4
3368 zout(2, j, nout1) = s1 + s25 + s34
3369 r = cos2*s25 + cos4*s34 + s1
3370 s = sin2*r25 + sin4*r34
3371 zout(2, j, nout2) = r + s
3372 zout(2, j, nout5) = r - s
3373 r = cos4*s25 + cos2*s34 + s1
3374 s = sin4*r25 - sin2*r34
3375 zout(2, j, nout3) = r + s
3376 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 bbs = isign*bb
3383 ia = 1
3384 nin1 = ia - after
3385 nout1 = ia - atn
3386 DO ib = 1, before
3387 nin1 = nin1 + after
3388 nin2 = nin1 + atb
3389 nin3 = nin2 + atb
3390 nin4 = nin3 + atb
3391 nin5 = nin4 + atb
3392 nin6 = nin5 + atb
3393 nout1 = nout1 + atn
3394 nout2 = nout1 + after
3395 nout3 = nout2 + after
3396 nout4 = nout3 + after
3397 nout5 = nout4 + after
3398 nout6 = nout5 + after
3399 DO j = 1, nfft
3400 r2 = zin(1, j, nin3)
3401 s2 = zin(2, j, nin3)
3402 r3 = zin(1, j, nin5)
3403 s3 = zin(2, j, nin5)
3404 r = r2 + r3
3405 s = s2 + s3
3406 r1 = zin(1, j, nin1)
3407 s1 = zin(2, j, nin1)
3408 ur1 = r + r1
3409 ui1 = s + s1
3410 r1 = r1 - 0.5_dp*r
3411 s1 = s1 - 0.5_dp*s
3412 r = r2 - r3
3413 s = s2 - s3
3414 ur2 = r1 - s*bbs
3415 ui2 = s1 + r*bbs
3416 ur3 = r1 + s*bbs
3417 ui3 = s1 - r*bbs
3418
3419 r2 = zin(1, j, nin6)
3420 s2 = zin(2, j, nin6)
3421 r3 = zin(1, j, nin2)
3422 s3 = zin(2, j, nin2)
3423 r = r2 + r3
3424 s = s2 + s3
3425 r1 = zin(1, j, nin4)
3426 s1 = zin(2, j, nin4)
3427 vr1 = r + r1
3428 vi1 = s + s1
3429 r1 = r1 - 0.5_dp*r
3430 s1 = s1 - 0.5_dp*s
3431 r = r2 - r3
3432 s = s2 - s3
3433 vr2 = r1 - s*bbs
3434 vi2 = s1 + r*bbs
3435 vr3 = r1 + s*bbs
3436 vi3 = s1 - r*bbs
3437
3438 zout(1, j, nout1) = ur1 + vr1
3439 zout(2, j, nout1) = ui1 + vi1
3440 zout(1, j, nout5) = ur2 + vr2
3441 zout(2, j, nout5) = ui2 + vi2
3442 zout(1, j, nout3) = ur3 + vr3
3443 zout(2, j, nout3) = ui3 + vi3
3444 zout(1, j, nout4) = ur1 - vr1
3445 zout(2, j, nout4) = ui1 - vi1
3446 zout(1, j, nout2) = ur2 - vr2
3447 zout(2, j, nout2) = ui2 - vi2
3448 zout(1, j, nout6) = ur3 - vr3
3449 zout(2, j, nout6) = ui3 - vi3
3450 END DO
3451 END DO
3452 ELSE
3453 cpabort('Error fftstp')
3454 END IF
3455
3456!-----------------------------------------------------------------------------!
3457
3458 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 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 mloop: DO i = 1, nt
3517 IF (n == idata(1, i)) THEN
3518 ic = 0
3519 DO j = 1, 6
3520 itt = idata(1 + j, i)
3521 IF (itt > 1) THEN
3522 ic = ic + 1
3523 now(j) = idata(1 + j, i)
3524 ELSE
3525 EXIT mloop
3526 END IF
3527 END DO
3528 EXIT mloop
3529 END IF
3530 IF (i == nt) THEN
3531 WRITE (*, '(A,i5,A)') " Value of ", n, &
3532 " not allowed for fft, allowed values are:"
3533 WRITE (*, '(15i5)') (idata(1, j), j=1, nt)
3534 cpabort('ctrig')
3535 END IF
3536 END DO mloop
3537
3538 after(1) = 1
3539 before(ic) = 1
3540 DO i = 2, ic
3541 after(i) = after(i - 1)*now(i - 1)
3542 before(ic - i + 1) = before(ic - i + 2)*now(ic - i + 2)
3543 END DO
3544
3545 twopi = 8._dp*atan(1._dp)
3546 angle = isign*twopi/real(n, dp)
3547 trig(1, 1) = 1._dp
3548 trig(2, 1) = 0._dp
3549 DO i = 1, n - 1
3550 trig(1, i + 1) = cos(real(i, dp)*angle)
3551 trig(2, i + 1) = sin(real(i, dp)*angle)
3552 END DO
3553
3554 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 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 b(1:n, 1:m) = a(1:n, 1:m)
3572 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 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 b(1:n, 1:m) = transpose(a(1:m, 1:n))
3590 END SUBROUTINE zgetmo
3591
3592! **************************************************************************************************
3593!> \brief ...
3594!> \param n ...
3595!> \param sc ...
3596!> \param a ...
3597! **************************************************************************************************
3598 SUBROUTINE scaled(n, sc, a)
3599 INTEGER :: n
3600 REAL(dp) :: sc
3601 COMPLEX(dp) :: a(n)
3602
3603 CALL dscal(n, sc, a, 1)
3604
3605 END SUBROUTINE scaled
3606
3607END MODULE mltfftsg_tools
fft_kinds
Defines the basic variable types.
Definition fft_kinds.F:13
fft_kinds::dp
integer, parameter, public dp
Definition fft_kinds.F:18
mltfftsg_tools
Definition mltfftsg_tools.F:9
mltfftsg_tools::mltfftsg
subroutine, public mltfftsg(transa, transb, a, ldax, lday, b, ldbx, ldby, n, m, isign, scale)
...
Definition mltfftsg_tools.F:41