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