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