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