Permalink
Fetching contributors…
Cannot retrieve contributors at this time
3606 lines (3518 sloc) 132 KB
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 - 2018 CP2K developers group !
!--------------------------------------------------------------------------------------------------!
! **************************************************************************************************
MODULE mltfftsg_tools
USE ISO_C_BINDING, ONLY: C_F_POINTER,&
C_LOC
USE fft_kinds, ONLY: dp
!$ USE OMP_LIB, ONLY: omp_get_num_threads, omp_get_thread_num
#include "../../base/base_uses.f90"
PRIVATE
INTEGER, PARAMETER :: ctrig_length = 1024
INTEGER, PARAMETER :: cache_size = 2048
PUBLIC :: mltfftsg
CONTAINS
! **************************************************************************************************
!> \brief ...
!> \param transa ...
!> \param transb ...
!> \param a ...
!> \param ldax ...
!> \param lday ...
!> \param b ...
!> \param ldbx ...
!> \param ldby ...
!> \param n ...
!> \param m ...
!> \param isign ...
!> \param scale ...
! **************************************************************************************************
SUBROUTINE mltfftsg(transa, transb, a, ldax, lday, b, ldbx, ldby, n, m, isign, scale)
CHARACTER(LEN=1), INTENT(IN) :: transa, transb
INTEGER, INTENT(IN) :: ldax, lday
COMPLEX(dp), INTENT(INOUT) :: a(ldax, lday)
INTEGER, INTENT(IN) :: ldbx, ldby
COMPLEX(dp), INTENT(INOUT) :: b(ldbx, ldby)
INTEGER, INTENT(IN) :: n, m, isign
REAL(dp), INTENT(IN) :: scale
COMPLEX(dp), ALLOCATABLE, DIMENSION(:, :, :) :: z
INTEGER :: after(20), before(20), chunk, i, ic, id, &
iend, inzee, isig, istart, iterations, &
itr, length, lot, nfft, now(20), &
num_threads
LOGICAL :: tscal
REAL(dp) :: trig(2, 1024)
! Variables
LENGTH = 2*(cache_size/4+1)
ISIG = -ISIGN
TSCAL = (ABS(SCALE-1._dp) > 1.e-12_dp)
CALL ctrig(N, TRIG, AFTER, BEFORE, NOW, ISIG, IC)
LOT = cache_size/(4*N)
LOT = LOT-MOD(LOT+1, 2)
LOT = MAX(1, LOT)
! initializations for serial mode
id = 0; num_threads = 1
!$OMP PARALLEL &
!$OMP PRIVATE ( id, istart, iend, nfft, i, inzee, itr) DEFAULT(NONE) &
!$OMP SHARED (NUM_THREADS,z,iterations,chunk,LOT,length,m,transa,isig, &
!$OMP before,after,now,trig,A,n,ldax,tscal,scale,ic,transb,ldbx,b)
!$OMP SINGLE
!$ num_threads = omp_get_num_threads()
ALLOCATE (Z(LENGTH, 2, 0:num_threads-1))
iterations = (M+LOT-1)/LOT
chunk = LOT*((iterations+num_threads-1)/num_threads)
!$OMP END SINGLE
!$OMP BARRIER
!$ id = omp_get_thread_num()
istart = id*chunk+1
iend = MIN((id+1)*chunk, M)
DO ITR = istart, iend, LOT
NFFT = MIN(M-ITR+1, LOT)
IF (TRANSA == 'N' .OR. TRANSA == 'n') THEN
CALL fftpre_cmplx(NFFT, NFFT, LDAX, LOT, N, A(1, ITR), Z(1, 1, id), &
TRIG, NOW(1), AFTER(1), BEFORE(1), ISIG)
ELSE
CALL fftstp_cmplx(LDAX, NFFT, N, LOT, N, A(ITR, 1), Z(1, 1, id), &
TRIG, NOW(1), AFTER(1), BEFORE(1), ISIG)
ENDIF
IF (TSCAL) THEN
IF (LOT == NFFT) THEN
CALL scaled(2*LOT*N, SCALE, Z(1, 1, id))
ELSE
DO I = 1, N
CALL scaled(2*NFFT, SCALE, Z(LOT*(I-1)+1, 1, id))
END DO
END IF
END IF
IF (IC .EQ. 1) THEN
IF (TRANSB == 'N' .OR. TRANSB == 'n') THEN
CALL zgetmo(Z(1, 1, id), LOT, NFFT, N, B(1, ITR), LDBX)
ELSE
CALL matmov(NFFT, N, Z(1, 1, id), LOT, B(ITR, 1), LDBX)
ENDIF
ELSE
INZEE = 1
DO I = 2, IC-1
CALL fftstp_cmplx(LOT, NFFT, N, LOT, N, Z(1, INZEE, id), &
Z(1, 3-INZEE, id), TRIG, NOW(I), AFTER(I), &
BEFORE(I), ISIG)
INZEE = 3-INZEE
ENDDO
IF (TRANSB == 'N' .OR. TRANSB == 'n') THEN
CALL fftrot_cmplx(LOT, NFFT, N, NFFT, LDBX, Z(1, INZEE, id), &
B(1, ITR), TRIG, NOW(IC), AFTER(IC), BEFORE(IC), ISIG)
ELSE
CALL fftstp_cmplx(LOT, NFFT, N, LDBX, N, Z(1, INZEE, id), &
B(ITR, 1), TRIG, NOW(IC), AFTER(IC), BEFORE(IC), ISIG)
ENDIF
ENDIF
ENDDO
!$OMP END PARALLEL
DEALLOCATE (Z)
IF (TRANSB == 'N' .OR. TRANSB == 'n') THEN
B(1:LDBX, M+1:LDBY) = CMPLX(0._dp, 0._dp, dp)
B(N+1:LDBX, 1:M) = CMPLX(0._dp, 0._dp, dp)
ELSE
B(1:LDBX, N+1:LDBY) = CMPLX(0._dp, 0._dp, dp)
B(M+1:LDBX, 1:N) = CMPLX(0._dp, 0._dp, dp)
ENDIF
END SUBROUTINE mltfftsg
! this formalizes what we have been assuming before, call with a complex(*) array, and passing to a real(2,*)
! **************************************************************************************************
!> \brief ...
!> \param mm ...
!> \param nfft ...
!> \param m ...
!> \param nn ...
!> \param n ...
!> \param zin ...
!> \param zout ...
!> \param trig ...
!> \param now ...
!> \param after ...
!> \param before ...
!> \param isign ...
! **************************************************************************************************
SUBROUTINE fftstp_cmplx(mm, nfft, m, nn, n, zin, zout, trig, now, after, before, isign)
INTEGER, INTENT(IN) :: mm, nfft, m, nn, n
COMPLEX(dp), DIMENSION(mm, m), INTENT(IN), TARGET :: zin
COMPLEX(dp), DIMENSION(nn, n), INTENT(INOUT), &
TARGET :: zout
REAL(dp), DIMENSION(2, ctrig_length), INTENT(IN) :: trig
INTEGER, INTENT(IN) :: now, after, before, isign
REAL(dp), DIMENSION(:, :, :), POINTER :: zin_real, zout_real
CALL C_F_POINTER(C_LOC(zin), zin_real, (/2, mm, m/))
CALL C_F_POINTER(C_LOC(zout), zout_real, (/2, nn, n/))
CALL fftstp(mm, nfft, m, nn, n, zin_real, zout_real, trig, now, after, before, isign)
END SUBROUTINE
! **************************************************************************************************
!> \brief ...
!> \param mm ...
!> \param nfft ...
!> \param m ...
!> \param nn ...
!> \param n ...
!> \param zin ...
!> \param zout ...
!> \param trig ...
!> \param now ...
!> \param after ...
!> \param before ...
!> \param isign ...
! **************************************************************************************************
SUBROUTINE fftpre_cmplx(mm, nfft, m, nn, n, zin, zout, trig, now, after, before, isign)
INTEGER, INTENT(IN) :: mm, nfft, m, nn, n
COMPLEX(dp), DIMENSION(m, mm), INTENT(IN), TARGET :: zin
COMPLEX(dp), DIMENSION(nn, n), INTENT(INOUT), &
TARGET :: zout
REAL(dp), DIMENSION(2, ctrig_length), INTENT(IN) :: trig
INTEGER, INTENT(IN) :: now, after, before, isign
REAL(dp), DIMENSION(:, :, :), POINTER :: zin_real, zout_real
CALL C_F_POINTER(C_LOC(zin), zin_real, (/2, mm, m/))
CALL C_F_POINTER(C_LOC(zout), zout_real, (/2, nn, n/))
CALL fftpre(mm, nfft, m, nn, n, zin_real, zout_real, trig, now, after, before, isign)
END SUBROUTINE
! **************************************************************************************************
!> \brief ...
!> \param mm ...
!> \param nfft ...
!> \param m ...
!> \param nn ...
!> \param n ...
!> \param zin ...
!> \param zout ...
!> \param trig ...
!> \param now ...
!> \param after ...
!> \param before ...
!> \param isign ...
! **************************************************************************************************
SUBROUTINE fftrot_cmplx(mm, nfft, m, nn, n, zin, zout, trig, now, after, before, isign)
USE fft_kinds, ONLY: dp
INTEGER, INTENT(IN) :: mm, nfft, m, nn, n
COMPLEX(dp), DIMENSION(mm, m), INTENT(IN), TARGET :: zin
COMPLEX(dp), DIMENSION(n, nn), INTENT(INOUT), &
TARGET :: zout
REAL(dp), DIMENSION(2, ctrig_length), INTENT(IN) :: trig
INTEGER, INTENT(IN) :: now, after, before, isign
REAL(dp), DIMENSION(:, :, :), POINTER :: zin_real, zout_real
CALL C_F_POINTER(C_LOC(zin), zin_real, (/2, mm, m/))
CALL C_F_POINTER(C_LOC(zout), zout_real, (/2, nn, n/))
CALL fftrot(mm, nfft, m, nn, n, zin_real, zout_real, trig, now, after, before, isign)
END SUBROUTINE
!-----------------------------------------------------------------------------!
! Copyright (C) Stefan Goedecker, Lausanne, Switzerland, August 1, 1991
! Copyright (C) Stefan Goedecker, Cornell University, Ithaca, USA, 1994
! Copyright (C) Stefan Goedecker, MPI Stuttgart, Germany, 1999
! This file is distributed under the terms of the
! GNU General Public License version 2 (or later),
! see http://www.gnu.org/copyleft/gpl.txt .
!-----------------------------------------------------------------------------!
! S. Goedecker: Rotating a three-dimensional array in optimal
! positions for vector processing: Case study for a three-dimensional Fast
! Fourier Transform, Comp. Phys. Commun. 76, 294 (1993)
! **************************************************************************************************
!> \brief ...
!> \param mm ...
!> \param nfft ...
!> \param m ...
!> \param nn ...
!> \param n ...
!> \param zin ...
!> \param zout ...
!> \param trig ...
!> \param now ...
!> \param after ...
!> \param before ...
!> \param isign ...
! **************************************************************************************************
SUBROUTINE fftrot(mm, nfft, m, nn, n, zin, zout, trig, now, after, before, isign)
USE fft_kinds, ONLY: dp
INTEGER, INTENT(IN) :: mm, nfft, m, nn, n
REAL(dp), DIMENSION(2, mm, m), INTENT(IN) :: zin
REAL(dp), DIMENSION(2, n, nn), INTENT(INOUT) :: zout
REAL(dp), DIMENSION(2, ctrig_length), INTENT(IN) :: trig
INTEGER, INTENT(IN) :: now, after, before, isign
REAL(dp), PARAMETER :: bb = 0.8660254037844387_dp, cos2 = 0.3090169943749474_dp, &
cos4 = -0.8090169943749474_dp, rt2i = 0.7071067811865475_dp, &
sin2p = 0.9510565162951536_dp, sin4p = 0.5877852522924731_dp
INTEGER :: atb, atn, ia, ias, ib, itrig, itt, j, &
nin1, nin2, nin3, nin4, nin5, nin6, &
nin7, nin8, nout1, nout2, nout3, &
nout4, nout5, nout6, nout7, nout8
REAL(dp) :: am, ap, bbs, bm, bp, ci2, ci3, ci4, ci5, cm, cp, cr2, cr3, cr4, cr5, dbl, dm, r, &
r1, r2, r25, r3, r34, r4, r5, r6, r7, r8, s, s1, s2, s25, s3, s34, s4, s5, s6, s7, s8, &
sin2, sin4, ui1, ui2, ui3, ur1, ur2, ur3, vi1, vi2, vi3, vr1, vr2, vr3
! sqrt(0.5)
! sqrt(3)/2
! cos(2*pi/5)
! cos(4*pi/5)
! sin(2*pi/5)
! sin(4*pi/5)
!-----------------------------------------------------------------------------!
atn = after*now
atb = after*before
IF (now == 4) THEN
IF (isign == 1) THEN
ia = 1
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nin4 = nin3+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
nout4 = nout3+after
DO j = 1, nfft
r1 = zin(1, j, nin1)
s1 = zin(2, j, nin1)
r2 = zin(1, j, nin2)
s2 = zin(2, j, nin2)
r3 = zin(1, j, nin3)
s3 = zin(2, j, nin3)
r4 = zin(1, j, nin4)
s4 = zin(2, j, nin4)
r = r1+r3
s = r2+r4
zout(1, nout1, j) = r+s
zout(1, nout3, j) = r-s
r = r1-r3
s = s2-s4
zout(1, nout2, j) = r-s
zout(1, nout4, j) = r+s
r = s1+s3
s = s2+s4
zout(2, nout1, j) = r+s
zout(2, nout3, j) = r-s
r = s1-s3
s = r2-r4
zout(2, nout2, j) = r+s
zout(2, nout4, j) = r-s
END DO
END DO
DO ia = 2, after
ias = ia-1
IF (2*ias == after) THEN
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nin4 = nin3+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
nout4 = nout3+after
DO j = 1, nfft
r1 = zin(1, j, nin1)
s1 = zin(2, j, nin1)
r = zin(1, j, nin2)
s = zin(2, j, nin2)
r2 = (r-s)*rt2i
s2 = (r+s)*rt2i
r3 = -zin(2, j, nin3)
s3 = zin(1, j, nin3)
r = zin(1, j, nin4)
s = zin(2, j, nin4)
r4 = -(r+s)*rt2i
s4 = (r-s)*rt2i
r = r1+r3
s = r2+r4
zout(1, nout1, j) = r+s
zout(1, nout3, j) = r-s
r = r1-r3
s = s2-s4
zout(1, nout2, j) = r-s
zout(1, nout4, j) = r+s
r = s1+s3
s = s2+s4
zout(2, nout1, j) = r+s
zout(2, nout3, j) = r-s
r = s1-s3
s = r2-r4
zout(2, nout2, j) = r+s
zout(2, nout4, j) = r-s
END DO
END DO
ELSE
itt = ias*before
itrig = itt+1
cr2 = trig(1, itrig)
ci2 = trig(2, itrig)
itrig = itrig+itt
cr3 = trig(1, itrig)
ci3 = trig(2, itrig)
itrig = itrig+itt
cr4 = trig(1, itrig)
ci4 = trig(2, itrig)
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nin4 = nin3+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
nout4 = nout3+after
DO j = 1, nfft
r1 = zin(1, j, nin1)
s1 = zin(2, j, nin1)
r = zin(1, j, nin2)
s = zin(2, j, nin2)
r2 = r*cr2-s*ci2
s2 = r*ci2+s*cr2
r = zin(1, j, nin3)
s = zin(2, j, nin3)
r3 = r*cr3-s*ci3
s3 = r*ci3+s*cr3
r = zin(1, j, nin4)
s = zin(2, j, nin4)
r4 = r*cr4-s*ci4
s4 = r*ci4+s*cr4
r = r1+r3
s = r2+r4
zout(1, nout1, j) = r+s
zout(1, nout3, j) = r-s
r = r1-r3
s = s2-s4
zout(1, nout2, j) = r-s
zout(1, nout4, j) = r+s
r = s1+s3
s = s2+s4
zout(2, nout1, j) = r+s
zout(2, nout3, j) = r-s
r = s1-s3
s = r2-r4
zout(2, nout2, j) = r+s
zout(2, nout4, j) = r-s
END DO
END DO
END IF
END DO
ELSE
ia = 1
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nin4 = nin3+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
nout4 = nout3+after
DO j = 1, nfft
r1 = zin(1, j, nin1)
s1 = zin(2, j, nin1)
r2 = zin(1, j, nin2)
s2 = zin(2, j, nin2)
r3 = zin(1, j, nin3)
s3 = zin(2, j, nin3)
r4 = zin(1, j, nin4)
s4 = zin(2, j, nin4)
r = r1+r3
s = r2+r4
zout(1, nout1, j) = r+s
zout(1, nout3, j) = r-s
r = r1-r3
s = s2-s4
zout(1, nout2, j) = r+s
zout(1, nout4, j) = r-s
r = s1+s3
s = s2+s4
zout(2, nout1, j) = r+s
zout(2, nout3, j) = r-s
r = s1-s3
s = r2-r4
zout(2, nout2, j) = r-s
zout(2, nout4, j) = r+s
END DO
END DO
DO ia = 2, after
ias = ia-1
IF (2*ias == after) THEN
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nin4 = nin3+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
nout4 = nout3+after
DO j = 1, nfft
r1 = zin(1, j, nin1)
s1 = zin(2, j, nin1)
r = zin(1, j, nin2)
s = zin(2, j, nin2)
r2 = (r+s)*rt2i
s2 = (s-r)*rt2i
r3 = zin(2, j, nin3)
s3 = -zin(1, j, nin3)
r = zin(1, j, nin4)
s = zin(2, j, nin4)
r4 = (s-r)*rt2i
s4 = -(r+s)*rt2i
r = r1+r3
s = r2+r4
zout(1, nout1, j) = r+s
zout(1, nout3, j) = r-s
r = r1-r3
s = s2-s4
zout(1, nout2, j) = r+s
zout(1, nout4, j) = r-s
r = s1+s3
s = s2+s4
zout(2, nout1, j) = r+s
zout(2, nout3, j) = r-s
r = s1-s3
s = r2-r4
zout(2, nout2, j) = r-s
zout(2, nout4, j) = r+s
END DO
END DO
ELSE
itt = ias*before
itrig = itt+1
cr2 = trig(1, itrig)
ci2 = trig(2, itrig)
itrig = itrig+itt
cr3 = trig(1, itrig)
ci3 = trig(2, itrig)
itrig = itrig+itt
cr4 = trig(1, itrig)
ci4 = trig(2, itrig)
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nin4 = nin3+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
nout4 = nout3+after
DO j = 1, nfft
r1 = zin(1, j, nin1)
s1 = zin(2, j, nin1)
r = zin(1, j, nin2)
s = zin(2, j, nin2)
r2 = r*cr2-s*ci2
s2 = r*ci2+s*cr2
r = zin(1, j, nin3)
s = zin(2, j, nin3)
r3 = r*cr3-s*ci3
s3 = r*ci3+s*cr3
r = zin(1, j, nin4)
s = zin(2, j, nin4)
r4 = r*cr4-s*ci4
s4 = r*ci4+s*cr4
r = r1+r3
s = r2+r4
zout(1, nout1, j) = r+s
zout(1, nout3, j) = r-s
r = r1-r3
s = s2-s4
zout(1, nout2, j) = r+s
zout(1, nout4, j) = r-s
r = s1+s3
s = s2+s4
zout(2, nout1, j) = r+s
zout(2, nout3, j) = r-s
r = s1-s3
s = r2-r4
zout(2, nout2, j) = r-s
zout(2, nout4, j) = r+s
END DO
END DO
END IF
END DO
END IF
ELSE IF (now == 8) THEN
IF (isign == -1) THEN
ia = 1
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nin4 = nin3+atb
nin5 = nin4+atb
nin6 = nin5+atb
nin7 = nin6+atb
nin8 = nin7+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
nout4 = nout3+after
nout5 = nout4+after
nout6 = nout5+after
nout7 = nout6+after
nout8 = nout7+after
DO j = 1, nfft
r1 = zin(1, j, nin1)
s1 = zin(2, j, nin1)
r2 = zin(1, j, nin2)
s2 = zin(2, j, nin2)
r3 = zin(1, j, nin3)
s3 = zin(2, j, nin3)
r4 = zin(1, j, nin4)
s4 = zin(2, j, nin4)
r5 = zin(1, j, nin5)
s5 = zin(2, j, nin5)
r6 = zin(1, j, nin6)
s6 = zin(2, j, nin6)
r7 = zin(1, j, nin7)
s7 = zin(2, j, nin7)
r8 = zin(1, j, nin8)
s8 = zin(2, j, nin8)
r = r1+r5
s = r3+r7
ap = r+s
am = r-s
r = r2+r6
s = r4+r8
bp = r+s
bm = r-s
r = s1+s5
s = s3+s7
cp = r+s
cm = r-s
r = s2+s6
s = s4+s8
dbl = r+s
dm = r-s
zout(1, nout1, j) = ap+bp
zout(2, nout1, j) = cp+dbl
zout(1, nout5, j) = ap-bp
zout(2, nout5, j) = cp-dbl
zout(1, nout3, j) = am+dm
zout(2, nout3, j) = cm-bm
zout(1, nout7, j) = am-dm
zout(2, nout7, j) = cm+bm
r = r1-r5
s = s3-s7
ap = r+s
am = r-s
r = s1-s5
s = r3-r7
bp = r+s
bm = r-s
r = s4-s8
s = r2-r6
cp = r+s
cm = r-s
r = s2-s6
s = r4-r8
dbl = r+s
dm = r-s
r = (cp+dm)*rt2i
s = (-cp+dm)*rt2i
cp = (cm+dbl)*rt2i
dbl = (cm-dbl)*rt2i
zout(1, nout2, j) = ap+r
zout(2, nout2, j) = bm+s
zout(1, nout6, j) = ap-r
zout(2, nout6, j) = bm-s
zout(1, nout4, j) = am+cp
zout(2, nout4, j) = bp+dbl
zout(1, nout8, j) = am-cp
zout(2, nout8, j) = bp-dbl
END DO
END DO
ELSE
ia = 1
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nin4 = nin3+atb
nin5 = nin4+atb
nin6 = nin5+atb
nin7 = nin6+atb
nin8 = nin7+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
nout4 = nout3+after
nout5 = nout4+after
nout6 = nout5+after
nout7 = nout6+after
nout8 = nout7+after
DO j = 1, nfft
r1 = zin(1, j, nin1)
s1 = zin(2, j, nin1)
r2 = zin(1, j, nin2)
s2 = zin(2, j, nin2)
r3 = zin(1, j, nin3)
s3 = zin(2, j, nin3)
r4 = zin(1, j, nin4)
s4 = zin(2, j, nin4)
r5 = zin(1, j, nin5)
s5 = zin(2, j, nin5)
r6 = zin(1, j, nin6)
s6 = zin(2, j, nin6)
r7 = zin(1, j, nin7)
s7 = zin(2, j, nin7)
r8 = zin(1, j, nin8)
s8 = zin(2, j, nin8)
r = r1+r5
s = r3+r7
ap = r+s
am = r-s
r = r2+r6
s = r4+r8
bp = r+s
bm = r-s
r = s1+s5
s = s3+s7
cp = r+s
cm = r-s
r = s2+s6
s = s4+s8
dbl = r+s
dm = r-s
zout(1, nout1, j) = ap+bp
zout(2, nout1, j) = cp+dbl
zout(1, nout5, j) = ap-bp
zout(2, nout5, j) = cp-dbl
zout(1, nout3, j) = am-dm
zout(2, nout3, j) = cm+bm
zout(1, nout7, j) = am+dm
zout(2, nout7, j) = cm-bm
r = r1-r5
s = -s3+s7
ap = r+s
am = r-s
r = s1-s5
s = r7-r3
bp = r+s
bm = r-s
r = -s4+s8
s = r2-r6
cp = r+s
cm = r-s
r = -s2+s6
s = r4-r8
dbl = r+s
dm = r-s
r = (cp+dm)*rt2i
s = (cp-dm)*rt2i
cp = (cm+dbl)*rt2i
dbl = (-cm+dbl)*rt2i
zout(1, nout2, j) = ap+r
zout(2, nout2, j) = bm+s
zout(1, nout6, j) = ap-r
zout(2, nout6, j) = bm-s
zout(1, nout4, j) = am+cp
zout(2, nout4, j) = bp+dbl
zout(1, nout8, j) = am-cp
zout(2, nout8, j) = bp-dbl
END DO
END DO
END IF
ELSE IF (now == 3) THEN
bbs = isign*bb
ia = 1
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
DO j = 1, nfft
r1 = zin(1, j, nin1)
s1 = zin(2, j, nin1)
r2 = zin(1, j, nin2)
s2 = zin(2, j, nin2)
r3 = zin(1, j, nin3)
s3 = zin(2, j, nin3)
r = r2+r3
s = s2+s3
zout(1, nout1, j) = r+r1
zout(2, nout1, j) = s+s1
r1 = r1-0.5_dp*r
s1 = s1-0.5_dp*s
r2 = bbs*(r2-r3)
s2 = bbs*(s2-s3)
zout(1, nout2, j) = r1-s2
zout(2, nout2, j) = s1+r2
zout(1, nout3, j) = r1+s2
zout(2, nout3, j) = s1-r2
END DO
END DO
DO ia = 2, after
ias = ia-1
IF (4*ias == 3*after) THEN
IF (isign == 1) THEN
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
DO j = 1, nfft
r1 = zin(1, j, nin1)
s1 = zin(2, j, nin1)
r2 = -zin(2, j, nin2)
s2 = zin(1, j, nin2)
r3 = -zin(1, j, nin3)
s3 = -zin(2, j, nin3)
r = r2+r3
s = s2+s3
zout(1, nout1, j) = r+r1
zout(2, nout1, j) = s+s1
r1 = r1-0.5_dp*r
s1 = s1-0.5_dp*s
r2 = bbs*(r2-r3)
s2 = bbs*(s2-s3)
zout(1, nout2, j) = r1-s2
zout(2, nout2, j) = s1+r2
zout(1, nout3, j) = r1+s2
zout(2, nout3, j) = s1-r2
END DO
END DO
ELSE
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
DO j = 1, nfft
r1 = zin(1, j, nin1)
s1 = zin(2, j, nin1)
r2 = zin(2, j, nin2)
s2 = -zin(1, j, nin2)
r3 = -zin(1, j, nin3)
s3 = -zin(2, j, nin3)
r = r2+r3
s = s2+s3
zout(1, nout1, j) = r+r1
zout(2, nout1, j) = s+s1
r1 = r1-0.5_dp*r
s1 = s1-0.5_dp*s
r2 = bbs*(r2-r3)
s2 = bbs*(s2-s3)
zout(1, nout2, j) = r1-s2
zout(2, nout2, j) = s1+r2
zout(1, nout3, j) = r1+s2
zout(2, nout3, j) = s1-r2
END DO
END DO
END IF
ELSE IF (8*ias == 3*after) THEN
IF (isign == 1) THEN
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
DO j = 1, nfft
r1 = zin(1, j, nin1)
s1 = zin(2, j, nin1)
r = zin(1, j, nin2)
s = zin(2, j, nin2)
r2 = (r-s)*rt2i
s2 = (r+s)*rt2i
r3 = -zin(2, j, nin3)
s3 = zin(1, j, nin3)
r = r2+r3
s = s2+s3
zout(1, nout1, j) = r+r1
zout(2, nout1, j) = s+s1
r1 = r1-0.5_dp*r
s1 = s1-0.5_dp*s
r2 = bbs*(r2-r3)
s2 = bbs*(s2-s3)
zout(1, nout2, j) = r1-s2
zout(2, nout2, j) = s1+r2
zout(1, nout3, j) = r1+s2
zout(2, nout3, j) = s1-r2
END DO
END DO
ELSE
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
DO j = 1, nfft
r1 = zin(1, j, nin1)
s1 = zin(2, j, nin1)
r = zin(1, j, nin2)
s = zin(2, j, nin2)
r2 = (r+s)*rt2i
s2 = (-r+s)*rt2i
r3 = zin(2, j, nin3)
s3 = -zin(1, j, nin3)
r = r2+r3
s = s2+s3
zout(1, nout1, j) = r+r1
zout(2, nout1, j) = s+s1
r1 = r1-0.5_dp*r
s1 = s1-0.5_dp*s
r2 = bbs*(r2-r3)
s2 = bbs*(s2-s3)
zout(1, nout2, j) = r1-s2
zout(2, nout2, j) = s1+r2
zout(1, nout3, j) = r1+s2
zout(2, nout3, j) = s1-r2
END DO
END DO
END IF
ELSE
itt = ias*before
itrig = itt+1
cr2 = trig(1, itrig)
ci2 = trig(2, itrig)
itrig = itrig+itt
cr3 = trig(1, itrig)
ci3 = trig(2, itrig)
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
DO j = 1, nfft
r1 = zin(1, j, nin1)
s1 = zin(2, j, nin1)
r = zin(1, j, nin2)
s = zin(2, j, nin2)
r2 = r*cr2-s*ci2
s2 = r*ci2+s*cr2
r = zin(1, j, nin3)
s = zin(2, j, nin3)
r3 = r*cr3-s*ci3
s3 = r*ci3+s*cr3
r = r2+r3
s = s2+s3
zout(1, nout1, j) = r+r1
zout(2, nout1, j) = s+s1
r1 = r1-0.5_dp*r
s1 = s1-0.5_dp*s
r2 = bbs*(r2-r3)
s2 = bbs*(s2-s3)
zout(1, nout2, j) = r1-s2
zout(2, nout2, j) = s1+r2
zout(1, nout3, j) = r1+s2
zout(2, nout3, j) = s1-r2
END DO
END DO
END IF
END DO
ELSE IF (now == 5) THEN
sin2 = isign*sin2p
sin4 = isign*sin4p
ia = 1
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nin4 = nin3+atb
nin5 = nin4+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
nout4 = nout3+after
nout5 = nout4+after
DO j = 1, nfft
r1 = zin(1, j, nin1)
s1 = zin(2, j, nin1)
r2 = zin(1, j, nin2)
s2 = zin(2, j, nin2)
r3 = zin(1, j, nin3)
s3 = zin(2, j, nin3)
r4 = zin(1, j, nin4)
s4 = zin(2, j, nin4)
r5 = zin(1, j, nin5)
s5 = zin(2, j, nin5)
r25 = r2+r5
r34 = r3+r4
s25 = s2-s5
s34 = s3-s4
zout(1, nout1, j) = r1+r25+r34
r = cos2*r25+cos4*r34+r1
s = sin2*s25+sin4*s34
zout(1, nout2, j) = r-s
zout(1, nout5, j) = r+s
r = cos4*r25+cos2*r34+r1
s = sin4*s25-sin2*s34
zout(1, nout3, j) = r-s
zout(1, nout4, j) = r+s
r25 = r2-r5
r34 = r3-r4
s25 = s2+s5
s34 = s3+s4
zout(2, nout1, j) = s1+s25+s34
r = cos2*s25+cos4*s34+s1
s = sin2*r25+sin4*r34
zout(2, nout2, j) = r+s
zout(2, nout5, j) = r-s
r = cos4*s25+cos2*s34+s1
s = sin4*r25-sin2*r34
zout(2, nout3, j) = r+s
zout(2, nout4, j) = r-s
END DO
END DO
DO ia = 2, after
ias = ia-1
IF (8*ias == 5*after) THEN
IF (isign == 1) THEN
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nin4 = nin3+atb
nin5 = nin4+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
nout4 = nout3+after
nout5 = nout4+after
DO j = 1, nfft
r1 = zin(1, j, nin1)
s1 = zin(2, j, nin1)
r = zin(1, j, nin2)
s = zin(2, j, nin2)
r2 = (r-s)*rt2i
s2 = (r+s)*rt2i
r3 = -zin(2, j, nin3)
s3 = zin(1, j, nin3)
r = zin(1, j, nin4)
s = zin(2, j, nin4)
r4 = -(r+s)*rt2i
s4 = (r-s)*rt2i
r5 = -zin(1, j, nin5)
s5 = -zin(2, j, nin5)
r25 = r2+r5
r34 = r3+r4
s25 = s2-s5
s34 = s3-s4
zout(1, nout1, j) = r1+r25+r34
r = cos2*r25+cos4*r34+r1
s = sin2*s25+sin4*s34
zout(1, nout2, j) = r-s
zout(1, nout5, j) = r+s
r = cos4*r25+cos2*r34+r1
s = sin4*s25-sin2*s34
zout(1, nout3, j) = r-s
zout(1, nout4, j) = r+s
r25 = r2-r5
r34 = r3-r4
s25 = s2+s5
s34 = s3+s4
zout(2, nout1, j) = s1+s25+s34
r = cos2*s25+cos4*s34+s1
s = sin2*r25+sin4*r34
zout(2, nout2, j) = r+s
zout(2, nout5, j) = r-s
r = cos4*s25+cos2*s34+s1
s = sin4*r25-sin2*r34
zout(2, nout3, j) = r+s
zout(2, nout4, j) = r-s
END DO
END DO
ELSE
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nin4 = nin3+atb
nin5 = nin4+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
nout4 = nout3+after
nout5 = nout4+after
DO j = 1, nfft
r1 = zin(1, j, nin1)
s1 = zin(2, j, nin1)
r = zin(1, j, nin2)
s = zin(2, j, nin2)
r2 = (r+s)*rt2i
s2 = (-r+s)*rt2i
r3 = zin(2, j, nin3)
s3 = -zin(1, j, nin3)
r = zin(1, j, nin4)
s = zin(2, j, nin4)
r4 = (s-r)*rt2i
s4 = -(r+s)*rt2i
r5 = -zin(1, j, nin5)
s5 = -zin(2, j, nin5)
r25 = r2+r5
r34 = r3+r4
s25 = s2-s5
s34 = s3-s4
zout(1, nout1, j) = r1+r25+r34
r = cos2*r25+cos4*r34+r1
s = sin2*s25+sin4*s34
zout(1, nout2, j) = r-s
zout(1, nout5, j) = r+s
r = cos4*r25+cos2*r34+r1
s = sin4*s25-sin2*s34
zout(1, nout3, j) = r-s
zout(1, nout4, j) = r+s
r25 = r2-r5
r34 = r3-r4
s25 = s2+s5
s34 = s3+s4
zout(2, nout1, j) = s1+s25+s34
r = cos2*s25+cos4*s34+s1
s = sin2*r25+sin4*r34
zout(2, nout2, j) = r+s
zout(2, nout5, j) = r-s
r = cos4*s25+cos2*s34+s1
s = sin4*r25-sin2*r34
zout(2, nout3, j) = r+s
zout(2, nout4, j) = r-s
END DO
END DO
END IF
ELSE
ias = ia-1
itt = ias*before
itrig = itt+1
cr2 = trig(1, itrig)
ci2 = trig(2, itrig)
itrig = itrig+itt
cr3 = trig(1, itrig)
ci3 = trig(2, itrig)
itrig = itrig+itt
cr4 = trig(1, itrig)
ci4 = trig(2, itrig)
itrig = itrig+itt
cr5 = trig(1, itrig)
ci5 = trig(2, itrig)
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nin4 = nin3+atb
nin5 = nin4+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
nout4 = nout3+after
nout5 = nout4+after
DO j = 1, nfft
r1 = zin(1, j, nin1)
s1 = zin(2, j, nin1)
r = zin(1, j, nin2)
s = zin(2, j, nin2)
r2 = r*cr2-s*ci2
s2 = r*ci2+s*cr2
r = zin(1, j, nin3)
s = zin(2, j, nin3)
r3 = r*cr3-s*ci3
s3 = r*ci3+s*cr3
r = zin(1, j, nin4)
s = zin(2, j, nin4)
r4 = r*cr4-s*ci4
s4 = r*ci4+s*cr4
r = zin(1, j, nin5)
s = zin(2, j, nin5)
r5 = r*cr5-s*ci5
s5 = r*ci5+s*cr5
r25 = r2+r5
r34 = r3+r4
s25 = s2-s5
s34 = s3-s4
zout(1, nout1, j) = r1+r25+r34
r = cos2*r25+cos4*r34+r1
s = sin2*s25+sin4*s34
zout(1, nout2, j) = r-s
zout(1, nout5, j) = r+s
r = cos4*r25+cos2*r34+r1
s = sin4*s25-sin2*s34
zout(1, nout3, j) = r-s
zout(1, nout4, j) = r+s
r25 = r2-r5
r34 = r3-r4
s25 = s2+s5
s34 = s3+s4
zout(2, nout1, j) = s1+s25+s34
r = cos2*s25+cos4*s34+s1
s = sin2*r25+sin4*r34
zout(2, nout2, j) = r+s
zout(2, nout5, j) = r-s
r = cos4*s25+cos2*s34+s1
s = sin4*r25-sin2*r34
zout(2, nout3, j) = r+s
zout(2, nout4, j) = r-s
END DO
END DO
END IF
END DO
ELSE IF (now == 6) THEN
bbs = isign*bb
ia = 1
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nin4 = nin3+atb
nin5 = nin4+atb
nin6 = nin5+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
nout4 = nout3+after
nout5 = nout4+after
nout6 = nout5+after
DO j = 1, nfft
r2 = zin(1, j, nin3)
s2 = zin(2, j, nin3)
r3 = zin(1, j, nin5)
s3 = zin(2, j, nin5)
r = r2+r3
s = s2+s3
r1 = zin(1, j, nin1)
s1 = zin(2, j, nin1)
ur1 = r+r1
ui1 = s+s1
r1 = r1-0.5_dp*r
s1 = s1-0.5_dp*s
r = r2-r3
s = s2-s3
ur2 = r1-s*bbs
ui2 = s1+r*bbs
ur3 = r1+s*bbs
ui3 = s1-r*bbs
r2 = zin(1, j, nin6)
s2 = zin(2, j, nin6)
r3 = zin(1, j, nin2)
s3 = zin(2, j, nin2)
r = r2+r3
s = s2+s3
r1 = zin(1, j, nin4)
s1 = zin(2, j, nin4)
vr1 = r+r1
vi1 = s+s1
r1 = r1-0.5_dp*r
s1 = s1-0.5_dp*s
r = r2-r3
s = s2-s3
vr2 = r1-s*bbs
vi2 = s1+r*bbs
vr3 = r1+s*bbs
vi3 = s1-r*bbs
zout(1, nout1, j) = ur1+vr1
zout(2, nout1, j) = ui1+vi1
zout(1, nout5, j) = ur2+vr2
zout(2, nout5, j) = ui2+vi2
zout(1, nout3, j) = ur3+vr3
zout(2, nout3, j) = ui3+vi3
zout(1, nout4, j) = ur1-vr1
zout(2, nout4, j) = ui1-vi1
zout(1, nout2, j) = ur2-vr2
zout(2, nout2, j) = ui2-vi2
zout(1, nout6, j) = ur3-vr3
zout(2, nout6, j) = ui3-vi3
END DO
END DO
ELSE
CPABORT('Error fftrot')
END IF
!-----------------------------------------------------------------------------!
END SUBROUTINE fftrot
!-----------------------------------------------------------------------------!
!-----------------------------------------------------------------------------!
! Copyright (C) Stefan Goedecker, Lausanne, Switzerland, August 1, 1991
! Copyright (C) Stefan Goedecker, Cornell University, Ithaca, USA, 1994
! Copyright (C) Stefan Goedecker, MPI Stuttgart, Germany, 1999
! This file is distributed under the terms of the
! GNU General Public License version 2 (or later),
! see http://www.gnu.org/copyleft/gpl.txt .
!-----------------------------------------------------------------------------!
! S. Goedecker: Rotating a three-dimensional array in optimal
! positions for vector processing: Case study for a three-dimensional Fast
! Fourier Transform, Comp. Phys. Commun. 76, 294 (1993)
! **************************************************************************************************
!> \brief ...
!> \param mm ...
!> \param nfft ...
!> \param m ...
!> \param nn ...
!> \param n ...
!> \param zin ...
!> \param zout ...
!> \param trig ...
!> \param now ...
!> \param after ...
!> \param before ...
!> \param isign ...
! **************************************************************************************************
SUBROUTINE fftpre(mm, nfft, m, nn, n, zin, zout, trig, now, after, before, isign)
INTEGER, INTENT(IN) :: mm, nfft, m, nn, n
REAL(dp), DIMENSION(2, m, mm), INTENT(IN) :: zin
REAL(dp), DIMENSION(2, nn, n), INTENT(INOUT) :: zout
REAL(dp), DIMENSION(2, ctrig_length), INTENT(IN) :: trig
INTEGER, INTENT(IN) :: now, after, before, isign
REAL(dp), PARAMETER :: bb = 0.8660254037844387_dp, cos2 = 0.3090169943749474_dp, &
cos4 = -0.8090169943749474_dp, rt2i = 0.7071067811865475_dp, &
sin2p = 0.9510565162951536_dp, sin4p = 0.5877852522924731_dp
INTEGER :: atb, atn, ia, ias, ib, itrig, itt, j, &
nin1, nin2, nin3, nin4, nin5, nin6, &
nin7, nin8, nout1, nout2, nout3, &
nout4, nout5, nout6, nout7, nout8
REAL(dp) :: am, ap, bbs, bm, bp, ci2, ci3, ci4, ci5, cm, cp, cr2, cr3, cr4, cr5, dbl, dm, r, &
r1, r2, r25, r3, r34, r4, r5, r6, r7, r8, s, s1, s2, s25, s3, s34, s4, s5, s6, s7, s8, &
sin2, sin4, ui1, ui2, ui3, ur1, ur2, ur3, vi1, vi2, vi3, vr1, vr2, vr3
! sqrt(0.5)
! sqrt(3)/2
! cos(2*pi/5)
! cos(4*pi/5)
! sin(2*pi/5)
! sin(4*pi/5)
!-----------------------------------------------------------------------------!
atn = after*now
atb = after*before
IF (now == 4) THEN
IF (isign == 1) THEN
ia = 1
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nin4 = nin3+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
nout4 = nout3+after
DO j = 1, nfft
r1 = zin(1, nin1, j)
s1 = zin(2, nin1, j)
r2 = zin(1, nin2, j)
s2 = zin(2, nin2, j)
r3 = zin(1, nin3, j)
s3 = zin(2, nin3, j)
r4 = zin(1, nin4, j)
s4 = zin(2, nin4, j)
r = r1+r3
s = r2+r4
zout(1, j, nout1) = r+s
zout(1, j, nout3) = r-s
r = r1-r3
s = s2-s4
zout(1, j, nout2) = r-s
zout(1, j, nout4) = r+s
r = s1+s3
s = s2+s4
zout(2, j, nout1) = r+s
zout(2, j, nout3) = r-s
r = s1-s3
s = r2-r4
zout(2, j, nout2) = r+s
zout(2, j, nout4) = r-s
END DO
END DO
DO ia = 2, after
ias = ia-1
IF (2*ias == after) THEN
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nin4 = nin3+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
nout4 = nout3+after
DO j = 1, nfft
r1 = zin(1, nin1, j)
s1 = zin(2, nin1, j)
r = zin(1, nin2, j)
s = zin(2, nin2, j)
r2 = (r-s)*rt2i
s2 = (r+s)*rt2i
r3 = -zin(2, nin3, j)
s3 = zin(1, nin3, j)
r = zin(1, nin4, j)
s = zin(2, nin4, j)
r4 = -(r+s)*rt2i
s4 = (r-s)*rt2i
r = r1+r3
s = r2+r4
zout(1, j, nout1) = r+s
zout(1, j, nout3) = r-s
r = r1-r3
s = s2-s4
zout(1, j, nout2) = r-s
zout(1, j, nout4) = r+s
r = s1+s3
s = s2+s4
zout(2, j, nout1) = r+s
zout(2, j, nout3) = r-s
r = s1-s3
s = r2-r4
zout(2, j, nout2) = r+s
zout(2, j, nout4) = r-s
END DO
END DO
ELSE
itt = ias*before
itrig = itt+1
cr2 = trig(1, itrig)
ci2 = trig(2, itrig)
itrig = itrig+itt
cr3 = trig(1, itrig)
ci3 = trig(2, itrig)
itrig = itrig+itt
cr4 = trig(1, itrig)
ci4 = trig(2, itrig)
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nin4 = nin3+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
nout4 = nout3+after
DO j = 1, nfft
r1 = zin(1, nin1, j)
s1 = zin(2, nin1, j)
r = zin(1, nin2, j)
s = zin(2, nin2, j)
r2 = r*cr2-s*ci2
s2 = r*ci2+s*cr2
r = zin(1, nin3, j)
s = zin(2, nin3, j)
r3 = r*cr3-s*ci3
s3 = r*ci3+s*cr3
r = zin(1, nin4, j)
s = zin(2, nin4, j)
r4 = r*cr4-s*ci4
s4 = r*ci4+s*cr4
r = r1+r3
s = r2+r4
zout(1, j, nout1) = r+s
zout(1, j, nout3) = r-s
r = r1-r3
s = s2-s4
zout(1, j, nout2) = r-s
zout(1, j, nout4) = r+s
r = s1+s3
s = s2+s4
zout(2, j, nout1) = r+s
zout(2, j, nout3) = r-s
r = s1-s3
s = r2-r4
zout(2, j, nout2) = r+s
zout(2, j, nout4) = r-s
END DO
END DO
END IF
END DO
ELSE
ia = 1
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nin4 = nin3+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
nout4 = nout3+after
DO j = 1, nfft
r1 = zin(1, nin1, j)
s1 = zin(2, nin1, j)
r2 = zin(1, nin2, j)
s2 = zin(2, nin2, j)
r3 = zin(1, nin3, j)
s3 = zin(2, nin3, j)
r4 = zin(1, nin4, j)
s4 = zin(2, nin4, j)
r = r1+r3
s = r2+r4
zout(1, j, nout1) = r+s
zout(1, j, nout3) = r-s
r = r1-r3
s = s2-s4
zout(1, j, nout2) = r+s
zout(1, j, nout4) = r-s
r = s1+s3
s = s2+s4
zout(2, j, nout1) = r+s
zout(2, j, nout3) = r-s
r = s1-s3
s = r2-r4
zout(2, j, nout2) = r-s
zout(2, j, nout4) = r+s
END DO
END DO
DO ia = 2, after
ias = ia-1
IF (2*ias == after) THEN
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nin4 = nin3+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
nout4 = nout3+after
DO j = 1, nfft
r1 = zin(1, nin1, j)
s1 = zin(2, nin1, j)
r = zin(1, nin2, j)
s = zin(2, nin2, j)
r2 = (r+s)*rt2i
s2 = (s-r)*rt2i
r3 = zin(2, nin3, j)
s3 = -zin(1, nin3, j)
r = zin(1, nin4, j)
s = zin(2, nin4, j)
r4 = (s-r)*rt2i
s4 = -(r+s)*rt2i
r = r1+r3
s = r2+r4
zout(1, j, nout1) = r+s
zout(1, j, nout3) = r-s
r = r1-r3
s = s2-s4
zout(1, j, nout2) = r+s
zout(1, j, nout4) = r-s
r = s1+s3
s = s2+s4
zout(2, j, nout1) = r+s
zout(2, j, nout3) = r-s
r = s1-s3
s = r2-r4
zout(2, j, nout2) = r-s
zout(2, j, nout4) = r+s
END DO
END DO
ELSE
itt = ias*before
itrig = itt+1
cr2 = trig(1, itrig)
ci2 = trig(2, itrig)
itrig = itrig+itt
cr3 = trig(1, itrig)
ci3 = trig(2, itrig)
itrig = itrig+itt
cr4 = trig(1, itrig)
ci4 = trig(2, itrig)
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nin4 = nin3+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
nout4 = nout3+after
DO j = 1, nfft
r1 = zin(1, nin1, j)
s1 = zin(2, nin1, j)
r = zin(1, nin2, j)
s = zin(2, nin2, j)
r2 = r*cr2-s*ci2
s2 = r*ci2+s*cr2
r = zin(1, nin3, j)
s = zin(2, nin3, j)
r3 = r*cr3-s*ci3
s3 = r*ci3+s*cr3
r = zin(1, nin4, j)
s = zin(2, nin4, j)
r4 = r*cr4-s*ci4
s4 = r*ci4+s*cr4
r = r1+r3
s = r2+r4
zout(1, j, nout1) = r+s
zout(1, j, nout3) = r-s
r = r1-r3
s = s2-s4
zout(1, j, nout2) = r+s
zout(1, j, nout4) = r-s
r = s1+s3
s = s2+s4
zout(2, j, nout1) = r+s
zout(2, j, nout3) = r-s
r = s1-s3
s = r2-r4
zout(2, j, nout2) = r-s
zout(2, j, nout4) = r+s
END DO
END DO
END IF
END DO
END IF
ELSE IF (now == 8) THEN
IF (isign == -1) THEN
ia = 1
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nin4 = nin3+atb
nin5 = nin4+atb
nin6 = nin5+atb
nin7 = nin6+atb
nin8 = nin7+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
nout4 = nout3+after
nout5 = nout4+after
nout6 = nout5+after
nout7 = nout6+after
nout8 = nout7+after
DO j = 1, nfft
r1 = zin(1, nin1, j)
s1 = zin(2, nin1, j)
r2 = zin(1, nin2, j)
s2 = zin(2, nin2, j)
r3 = zin(1, nin3, j)
s3 = zin(2, nin3, j)
r4 = zin(1, nin4, j)
s4 = zin(2, nin4, j)
r5 = zin(1, nin5, j)
s5 = zin(2, nin5, j)
r6 = zin(1, nin6, j)
s6 = zin(2, nin6, j)
r7 = zin(1, nin7, j)
s7 = zin(2, nin7, j)
r8 = zin(1, nin8, j)
s8 = zin(2, nin8, j)
r = r1+r5
s = r3+r7
ap = r+s
am = r-s
r = r2+r6
s = r4+r8
bp = r+s
bm = r-s
r = s1+s5
s = s3+s7
cp = r+s
cm = r-s
r = s2+s6
s = s4+s8
dbl = r+s
dm = r-s
zout(1, j, nout1) = ap+bp
zout(2, j, nout1) = cp+dbl
zout(1, j, nout5) = ap-bp
zout(2, j, nout5) = cp-dbl
zout(1, j, nout3) = am+dm
zout(2, j, nout3) = cm-bm
zout(1, j, nout7) = am-dm
zout(2, j, nout7) = cm+bm
r = r1-r5
s = s3-s7
ap = r+s
am = r-s
r = s1-s5
s = r3-r7
bp = r+s
bm = r-s
r = s4-s8
s = r2-r6
cp = r+s
cm = r-s
r = s2-s6
s = r4-r8
dbl = r+s
dm = r-s
r = (cp+dm)*rt2i
s = (-cp+dm)*rt2i
cp = (cm+dbl)*rt2i
dbl = (cm-dbl)*rt2i
zout(1, j, nout2) = ap+r
zout(2, j, nout2) = bm+s
zout(1, j, nout6) = ap-r
zout(2, j, nout6) = bm-s
zout(1, j, nout4) = am+cp
zout(2, j, nout4) = bp+dbl
zout(1, j, nout8) = am-cp
zout(2, j, nout8) = bp-dbl
END DO
END DO
ELSE
ia = 1
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nin4 = nin3+atb
nin5 = nin4+atb
nin6 = nin5+atb
nin7 = nin6+atb
nin8 = nin7+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
nout4 = nout3+after
nout5 = nout4+after
nout6 = nout5+after
nout7 = nout6+after
nout8 = nout7+after
DO j = 1, nfft
r1 = zin(1, nin1, j)
s1 = zin(2, nin1, j)
r2 = zin(1, nin2, j)
s2 = zin(2, nin2, j)
r3 = zin(1, nin3, j)
s3 = zin(2, nin3, j)
r4 = zin(1, nin4, j)
s4 = zin(2, nin4, j)
r5 = zin(1, nin5, j)
s5 = zin(2, nin5, j)
r6 = zin(1, nin6, j)
s6 = zin(2, nin6, j)
r7 = zin(1, nin7, j)
s7 = zin(2, nin7, j)
r8 = zin(1, nin8, j)
s8 = zin(2, nin8, j)
r = r1+r5
s = r3+r7
ap = r+s
am = r-s
r = r2+r6
s = r4+r8
bp = r+s
bm = r-s
r = s1+s5
s = s3+s7
cp = r+s
cm = r-s
r = s2+s6
s = s4+s8
dbl = r+s
dm = r-s
zout(1, j, nout1) = ap+bp
zout(2, j, nout1) = cp+dbl
zout(1, j, nout5) = ap-bp
zout(2, j, nout5) = cp-dbl
zout(1, j, nout3) = am-dm
zout(2, j, nout3) = cm+bm
zout(1, j, nout7) = am+dm
zout(2, j, nout7) = cm-bm
r = r1-r5
s = -s3+s7
ap = r+s
am = r-s
r = s1-s5
s = r7-r3
bp = r+s
bm = r-s
r = -s4+s8
s = r2-r6
cp = r+s
cm = r-s
r = -s2+s6
s = r4-r8
dbl = r+s
dm = r-s
r = (cp+dm)*rt2i
s = (cp-dm)*rt2i
cp = (cm+dbl)*rt2i
dbl = (-cm+dbl)*rt2i
zout(1, j, nout2) = ap+r
zout(2, j, nout2) = bm+s
zout(1, j, nout6) = ap-r
zout(2, j, nout6) = bm-s
zout(1, j, nout4) = am+cp
zout(2, j, nout4) = bp+dbl
zout(1, j, nout8) = am-cp
zout(2, j, nout8) = bp-dbl
END DO
END DO
END IF
ELSE IF (now == 3) THEN
ia = 1
nin1 = ia-after
nout1 = ia-atn
bbs = isign*bb
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
DO j = 1, nfft
r1 = zin(1, nin1, j)
s1 = zin(2, nin1, j)
r2 = zin(1, nin2, j)
s2 = zin(2, nin2, j)
r3 = zin(1, nin3, j)
s3 = zin(2, nin3, j)
r = r2+r3
s = s2+s3
zout(1, j, nout1) = r+r1
zout(2, j, nout1) = s+s1
r1 = r1-0.5_dp*r
s1 = s1-0.5_dp*s
r2 = bbs*(r2-r3)
s2 = bbs*(s2-s3)
zout(1, j, nout2) = r1-s2
zout(2, j, nout2) = s1+r2
zout(1, j, nout3) = r1+s2
zout(2, j, nout3) = s1-r2
END DO
END DO
DO ia = 2, after
ias = ia-1
IF (4*ias == 3*after) THEN
IF (isign == 1) THEN
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
DO j = 1, nfft
r1 = zin(1, nin1, j)
s1 = zin(2, nin1, j)
r2 = -zin(2, nin2, j)
s2 = zin(1, nin2, j)
r3 = -zin(1, nin3, j)
s3 = -zin(2, nin3, j)
r = r2+r3
s = s2+s3
zout(1, j, nout1) = r+r1
zout(2, j, nout1) = s+s1
r1 = r1-0.5_dp*r
s1 = s1-0.5_dp*s
r2 = bbs*(r2-r3)
s2 = bbs*(s2-s3)
zout(1, j, nout2) = r1-s2
zout(2, j, nout2) = s1+r2
zout(1, j, nout3) = r1+s2
zout(2, j, nout3) = s1-r2
END DO
END DO
ELSE
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
DO j = 1, nfft
r1 = zin(1, nin1, j)
s1 = zin(2, nin1, j)
r2 = zin(2, nin2, j)
s2 = -zin(1, nin2, j)
r3 = -zin(1, nin3, j)
s3 = -zin(2, nin3, j)
r = r2+r3
s = s2+s3
zout(1, j, nout1) = r+r1
zout(2, j, nout1) = s+s1
r1 = r1-0.5_dp*r
s1 = s1-0.5_dp*s
r2 = bbs*(r2-r3)
s2 = bbs*(s2-s3)
zout(1, j, nout2) = r1-s2
zout(2, j, nout2) = s1+r2
zout(1, j, nout3) = r1+s2
zout(2, j, nout3) = s1-r2
END DO
END DO
END IF
ELSE IF (8*ias == 3*after) THEN
IF (isign == 1) THEN
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
DO j = 1, nfft
r1 = zin(1, nin1, j)
s1 = zin(2, nin1, j)
r = zin(1, nin2, j)
s = zin(2, nin2, j)
r2 = (r-s)*rt2i
s2 = (r+s)*rt2i
r3 = -zin(2, nin3, j)
s3 = zin(1, nin3, j)
r = r2+r3
s = s2+s3
zout(1, j, nout1) = r+r1
zout(2, j, nout1) = s+s1
r1 = r1-0.5_dp*r
s1 = s1-0.5_dp*s
r2 = bbs*(r2-r3)
s2 = bbs*(s2-s3)
zout(1, j, nout2) = r1-s2
zout(2, j, nout2) = s1+r2
zout(1, j, nout3) = r1+s2
zout(2, j, nout3) = s1-r2
END DO
END DO
ELSE
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
DO j = 1, nfft
r1 = zin(1, nin1, j)
s1 = zin(2, nin1, j)
r = zin(1, nin2, j)
s = zin(2, nin2, j)
r2 = (r+s)*rt2i
s2 = (-r+s)*rt2i
r3 = zin(2, nin3, j)
s3 = -zin(1, nin3, j)
r = r2+r3
s = s2+s3
zout(1, j, nout1) = r+r1
zout(2, j, nout1) = s+s1
r1 = r1-0.5_dp*r
s1 = s1-0.5_dp*s
r2 = bbs*(r2-r3)
s2 = bbs*(s2-s3)
zout(1, j, nout2) = r1-s2
zout(2, j, nout2) = s1+r2
zout(1, j, nout3) = r1+s2
zout(2, j, nout3) = s1-r2
END DO
END DO
END IF
ELSE
itt = ias*before
itrig = itt+1
cr2 = trig(1, itrig)
ci2 = trig(2, itrig)
itrig = itrig+itt
cr3 = trig(1, itrig)
ci3 = trig(2, itrig)
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
DO j = 1, nfft
r1 = zin(1, nin1, j)
s1 = zin(2, nin1, j)
r = zin(1, nin2, j)
s = zin(2, nin2, j)
r2 = r*cr2-s*ci2
s2 = r*ci2+s*cr2
r = zin(1, nin3, j)
s = zin(2, nin3, j)
r3 = r*cr3-s*ci3
s3 = r*ci3+s*cr3
r = r2+r3
s = s2+s3
zout(1, j, nout1) = r+r1
zout(2, j, nout1) = s+s1
r1 = r1-0.5_dp*r
s1 = s1-0.5_dp*s
r2 = bbs*(r2-r3)
s2 = bbs*(s2-s3)
zout(1, j, nout2) = r1-s2
zout(2, j, nout2) = s1+r2
zout(1, j, nout3) = r1+s2
zout(2, j, nout3) = s1-r2
END DO
END DO
END IF
END DO
ELSE IF (now == 5) THEN
sin2 = isign*sin2p
sin4 = isign*sin4p
ia = 1
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nin4 = nin3+atb
nin5 = nin4+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
nout4 = nout3+after
nout5 = nout4+after
DO j = 1, nfft
r1 = zin(1, nin1, j)
s1 = zin(2, nin1, j)
r2 = zin(1, nin2, j)
s2 = zin(2, nin2, j)
r3 = zin(1, nin3, j)
s3 = zin(2, nin3, j)
r4 = zin(1, nin4, j)
s4 = zin(2, nin4, j)
r5 = zin(1, nin5, j)
s5 = zin(2, nin5, j)
r25 = r2+r5
r34 = r3+r4
s25 = s2-s5
s34 = s3-s4
zout(1, j, nout1) = r1+r25+r34
r = cos2*r25+cos4*r34+r1
s = sin2*s25+sin4*s34
zout(1, j, nout2) = r-s
zout(1, j, nout5) = r+s
r = cos4*r25+cos2*r34+r1
s = sin4*s25-sin2*s34
zout(1, j, nout3) = r-s
zout(1, j, nout4) = r+s
r25 = r2-r5
r34 = r3-r4
s25 = s2+s5
s34 = s3+s4
zout(2, j, nout1) = s1+s25+s34
r = cos2*s25+cos4*s34+s1
s = sin2*r25+sin4*r34
zout(2, j, nout2) = r+s
zout(2, j, nout5) = r-s
r = cos4*s25+cos2*s34+s1
s = sin4*r25-sin2*r34
zout(2, j, nout3) = r+s
zout(2, j, nout4) = r-s
END DO
END DO
DO ia = 2, after
ias = ia-1
IF (8*ias == 5*after) THEN
IF (isign == 1) THEN
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nin4 = nin3+atb
nin5 = nin4+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
nout4 = nout3+after
nout5 = nout4+after
DO j = 1, nfft
r1 = zin(1, nin1, j)
s1 = zin(2, nin1, j)
r = zin(1, nin2, j)
s = zin(2, nin2, j)
r2 = (r-s)*rt2i
s2 = (r+s)*rt2i
r3 = -zin(2, nin3, j)
s3 = zin(1, nin3, j)
r = zin(1, nin4, j)
s = zin(2, nin4, j)
r4 = -(r+s)*rt2i
s4 = (r-s)*rt2i
r5 = -zin(1, nin5, j)
s5 = -zin(2, nin5, j)
r25 = r2+r5
r34 = r3+r4
s25 = s2-s5
s34 = s3-s4
zout(1, j, nout1) = r1+r25+r34
r = cos2*r25+cos4*r34+r1
s = sin2*s25+sin4*s34
zout(1, j, nout2) = r-s
zout(1, j, nout5) = r+s
r = cos4*r25+cos2*r34+r1
s = sin4*s25-sin2*s34
zout(1, j, nout3) = r-s
zout(1, j, nout4) = r+s
r25 = r2-r5
r34 = r3-r4
s25 = s2+s5
s34 = s3+s4
zout(2, j, nout1) = s1+s25+s34
r = cos2*s25+cos4*s34+s1
s = sin2*r25+sin4*r34
zout(2, j, nout2) = r+s
zout(2, j, nout5) = r-s
r = cos4*s25+cos2*s34+s1
s = sin4*r25-sin2*r34
zout(2, j, nout3) = r+s
zout(2, j, nout4) = r-s
END DO
END DO
ELSE
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nin4 = nin3+atb
nin5 = nin4+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
nout4 = nout3+after
nout5 = nout4+after
DO j = 1, nfft
r1 = zin(1, nin1, j)
s1 = zin(2, nin1, j)
r = zin(1, nin2, j)
s = zin(2, nin2, j)
r2 = (r+s)*rt2i
s2 = (-r+s)*rt2i
r3 = zin(2, nin3, j)
s3 = -zin(1, nin3, j)
r = zin(1, nin4, j)
s = zin(2, nin4, j)
r4 = (s-r)*rt2i
s4 = -(r+s)*rt2i
r5 = -zin(1, nin5, j)
s5 = -zin(2, nin5, j)
r25 = r2+r5
r34 = r3+r4
s25 = s2-s5
s34 = s3-s4
zout(1, j, nout1) = r1+r25+r34
r = cos2*r25+cos4*r34+r1
s = sin2*s25+sin4*s34
zout(1, j, nout2) = r-s
zout(1, j, nout5) = r+s
r = cos4*r25+cos2*r34+r1
s = sin4*s25-sin2*s34
zout(1, j, nout3) = r-s
zout(1, j, nout4) = r+s
r25 = r2-r5
r34 = r3-r4
s25 = s2+s5
s34 = s3+s4
zout(2, j, nout1) = s1+s25+s34
r = cos2*s25+cos4*s34+s1
s = sin2*r25+sin4*r34
zout(2, j, nout2) = r+s
zout(2, j, nout5) = r-s
r = cos4*s25+cos2*s34+s1
s = sin4*r25-sin2*r34
zout(2, j, nout3) = r+s
zout(2, j, nout4) = r-s
END DO
END DO
END IF
ELSE
ias = ia-1
itt = ias*before
itrig = itt+1
cr2 = trig(1, itrig)
ci2 = trig(2, itrig)
itrig = itrig+itt
cr3 = trig(1, itrig)
ci3 = trig(2, itrig)
itrig = itrig+itt
cr4 = trig(1, itrig)
ci4 = trig(2, itrig)
itrig = itrig+itt
cr5 = trig(1, itrig)
ci5 = trig(2, itrig)
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nin4 = nin3+atb
nin5 = nin4+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
nout4 = nout3+after
nout5 = nout4+after
DO j = 1, nfft
r1 = zin(1, nin1, j)
s1 = zin(2, nin1, j)
r = zin(1, nin2, j)
s = zin(2, nin2, j)
r2 = r*cr2-s*ci2
s2 = r*ci2+s*cr2
r = zin(1, nin3, j)
s = zin(2, nin3, j)
r3 = r*cr3-s*ci3
s3 = r*ci3+s*cr3
r = zin(1, nin4, j)
s = zin(2, nin4, j)
r4 = r*cr4-s*ci4
s4 = r*ci4+s*cr4
r = zin(1, nin5, j)
s = zin(2, nin5, j)
r5 = r*cr5-s*ci5
s5 = r*ci5+s*cr5
r25 = r2+r5
r34 = r3+r4
s25 = s2-s5
s34 = s3-s4
zout(1, j, nout1) = r1+r25+r34
r = cos2*r25+cos4*r34+r1
s = sin2*s25+sin4*s34
zout(1, j, nout2) = r-s
zout(1, j, nout5) = r+s
r = cos4*r25+cos2*r34+r1
s = sin4*s25-sin2*s34
zout(1, j, nout3) = r-s
zout(1, j, nout4) = r+s
r25 = r2-r5
r34 = r3-r4
s25 = s2+s5
s34 = s3+s4
zout(2, j, nout1) = s1+s25+s34
r = cos2*s25+cos4*s34+s1
s = sin2*r25+sin4*r34
zout(2, j, nout2) = r+s
zout(2, j, nout5) = r-s
r = cos4*s25+cos2*s34+s1
s = sin4*r25-sin2*r34
zout(2, j, nout3) = r+s
zout(2, j, nout4) = r-s
END DO
END DO
END IF
END DO
ELSE IF (now == 6) THEN
bbs = isign*bb
ia = 1
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nin4 = nin3+atb
nin5 = nin4+atb
nin6 = nin5+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
nout4 = nout3+after
nout5 = nout4+after
nout6 = nout5+after
DO j = 1, nfft
r2 = zin(1, nin3, j)
s2 = zin(2, nin3, j)
r3 = zin(1, nin5, j)
s3 = zin(2, nin5, j)
r = r2+r3
s = s2+s3
r1 = zin(1, nin1, j)
s1 = zin(2, nin1, j)
ur1 = r+r1
ui1 = s+s1
r1 = r1-0.5_dp*r
s1 = s1-0.5_dp*s
r = r2-r3
s = s2-s3
ur2 = r1-s*bbs
ui2 = s1+r*bbs
ur3 = r1+s*bbs
ui3 = s1-r*bbs
r2 = zin(1, nin6, j)
s2 = zin(2, nin6, j)
r3 = zin(1, nin2, j)
s3 = zin(2, nin2, j)
r = r2+r3
s = s2+s3
r1 = zin(1, nin4, j)
s1 = zin(2, nin4, j)
vr1 = r+r1
vi1 = s+s1
r1 = r1-0.5_dp*r
s1 = s1-0.5_dp*s
r = r2-r3
s = s2-s3
vr2 = r1-s*bbs
vi2 = s1+r*bbs
vr3 = r1+s*bbs
vi3 = s1-r*bbs
zout(1, j, nout1) = ur1+vr1
zout(2, j, nout1) = ui1+vi1
zout(1, j, nout5) = ur2+vr2
zout(2, j, nout5) = ui2+vi2
zout(1, j, nout3) = ur3+vr3
zout(2, j, nout3) = ui3+vi3
zout(1, j, nout4) = ur1-vr1
zout(2, j, nout4) = ui1-vi1
zout(1, j, nout2) = ur2-vr2
zout(2, j, nout2) = ui2-vi2
zout(1, j, nout6) = ur3-vr3
zout(2, j, nout6) = ui3-vi3
END DO
END DO
ELSE
CPABORT('Error fftpre')
END IF
!-----------------------------------------------------------------------------!
END SUBROUTINE fftpre
!-----------------------------------------------------------------------------!
!-----------------------------------------------------------------------------!
! Copyright (C) Stefan Goedecker, Lausanne, Switzerland, August 1, 1991
! Copyright (C) Stefan Goedecker, Cornell University, Ithaca, USA, 1994
! Copyright (C) Stefan Goedecker, MPI Stuttgart, Germany, 1999
! This file is distributed under the terms of the
! GNU General Public License version 2 (or later),
! see http://www.gnu.org/copyleft/gpl.txt .
!-----------------------------------------------------------------------------!
! S. Goedecker: Rotating a three-dimensional array in optimal
! positions for vector processing: Case study for a three-dimensional Fast
! Fourier Transform, Comp. Phys. Commun. 76, 294 (1993)
! **************************************************************************************************
!> \brief ...
!> \param mm ...
!> \param nfft ...
!> \param m ...
!> \param nn ...
!> \param n ...
!> \param zin ...
!> \param zout ...
!> \param trig ...
!> \param now ...
!> \param after ...
!> \param before ...
!> \param isign ...
! **************************************************************************************************
SUBROUTINE fftstp(mm, nfft, m, nn, n, zin, zout, trig, now, after, before, isign)
INTEGER, INTENT(IN) :: mm, nfft, m, nn, n
REAL(dp), DIMENSION(2, mm, m), INTENT(IN) :: zin
REAL(dp), DIMENSION(2, nn, n), INTENT(INOUT) :: zout
REAL(dp), DIMENSION(2, ctrig_length), INTENT(IN) :: trig
INTEGER, INTENT(IN) :: now, after, before, isign
REAL(dp), PARAMETER :: bb = 0.8660254037844387_dp, cos2 = 0.3090169943749474_dp, &
cos4 = -0.8090169943749474_dp, rt2i = 0.7071067811865475_dp, &
sin2p = 0.9510565162951536_dp, sin4p = 0.5877852522924731_dp
INTEGER :: atb, atn, ia, ias, ib, itrig, itt, j, &
nin1, nin2, nin3, nin4, nin5, nin6, &
nin7, nin8, nout1, nout2, nout3, &
nout4, nout5, nout6, nout7, nout8
REAL(dp) :: am, ap, bbs, bm, bp, ci2, ci3, ci4, ci5, cm, cp, cr2, cr3, cr4, cr5, dbl, dm, r, &
r1, r2, r25, r3, r34, r4, r5, r6, r7, r8, s, s1, s2, s25, s3, s34, s4, s5, s6, s7, s8, &
sin2, sin4, ui1, ui2, ui3, ur1, ur2, ur3, vi1, vi2, vi3, vr1, vr2, vr3
! sqrt(0.5)
! sqrt(3)/2
! cos(2*pi/5)
! cos(4*pi/5)
! sin(2*pi/5)
! sin(4*pi/5)
!-----------------------------------------------------------------------------!
atn = after*now
atb = after*before
IF (now == 4) THEN
IF (isign == 1) THEN
ia = 1
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nin4 = nin3+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
nout4 = nout3+after
DO j = 1, nfft
r1 = zin(1, j, nin1)
s1 = zin(2, j, nin1)
r2 = zin(1, j, nin2)
s2 = zin(2, j, nin2)
r3 = zin(1, j, nin3)
s3 = zin(2, j, nin3)
r4 = zin(1, j, nin4)
s4 = zin(2, j, nin4)
r = r1+r3
s = r2+r4
zout(1, j, nout1) = r+s
zout(1, j, nout3) = r-s
r = r1-r3
s = s2-s4
zout(1, j, nout2) = r-s
zout(1, j, nout4) = r+s
r = s1+s3
s = s2+s4
zout(2, j, nout1) = r+s
zout(2, j, nout3) = r-s
r = s1-s3
s = r2-r4
zout(2, j, nout2) = r+s
zout(2, j, nout4) = r-s
END DO
END DO
DO ia = 2, after
ias = ia-1
IF (2*ias == after) THEN
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nin4 = nin3+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
nout4 = nout3+after
DO j = 1, nfft
r1 = zin(1, j, nin1)
s1 = zin(2, j, nin1)
r = zin(1, j, nin2)
s = zin(2, j, nin2)
r2 = (r-s)*rt2i
s2 = (r+s)*rt2i
r3 = -zin(2, j, nin3)
s3 = zin(1, j, nin3)
r = zin(1, j, nin4)
s = zin(2, j, nin4)
r4 = -(r+s)*rt2i
s4 = (r-s)*rt2i
r = r1+r3
s = r2+r4
zout(1, j, nout1) = r+s
zout(1, j, nout3) = r-s
r = r1-r3
s = s2-s4
zout(1, j, nout2) = r-s
zout(1, j, nout4) = r+s
r = s1+s3
s = s2+s4
zout(2, j, nout1) = r+s
zout(2, j, nout3) = r-s
r = s1-s3
s = r2-r4
zout(2, j, nout2) = r+s
zout(2, j, nout4) = r-s
END DO
END DO
ELSE
itt = ias*before
itrig = itt+1
cr2 = trig(1, itrig)
ci2 = trig(2, itrig)
itrig = itrig+itt
cr3 = trig(1, itrig)
ci3 = trig(2, itrig)
itrig = itrig+itt
cr4 = trig(1, itrig)
ci4 = trig(2, itrig)
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nin4 = nin3+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
nout4 = nout3+after
DO j = 1, nfft
r1 = zin(1, j, nin1)
s1 = zin(2, j, nin1)
r = zin(1, j, nin2)
s = zin(2, j, nin2)
r2 = r*cr2-s*ci2
s2 = r*ci2+s*cr2
r = zin(1, j, nin3)
s = zin(2, j, nin3)
r3 = r*cr3-s*ci3
s3 = r*ci3+s*cr3
r = zin(1, j, nin4)
s = zin(2, j, nin4)
r4 = r*cr4-s*ci4
s4 = r*ci4+s*cr4
r = r1+r3
s = r2+r4
zout(1, j, nout1) = r+s
zout(1, j, nout3) = r-s
r = r1-r3
s = s2-s4
zout(1, j, nout2) = r-s
zout(1, j, nout4) = r+s
r = s1+s3
s = s2+s4
zout(2, j, nout1) = r+s
zout(2, j, nout3) = r-s
r = s1-s3
s = r2-r4
zout(2, j, nout2) = r+s
zout(2, j, nout4) = r-s
END DO
END DO
END IF
END DO
ELSE
ia = 1
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nin4 = nin3+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
nout4 = nout3+after
DO j = 1, nfft
r1 = zin(1, j, nin1)
s1 = zin(2, j, nin1)
r2 = zin(1, j, nin2)
s2 = zin(2, j, nin2)
r3 = zin(1, j, nin3)
s3 = zin(2, j, nin3)
r4 = zin(1, j, nin4)
s4 = zin(2, j, nin4)
r = r1+r3
s = r2+r4
zout(1, j, nout1) = r+s
zout(1, j, nout3) = r-s
r = r1-r3
s = s2-s4
zout(1, j, nout2) = r+s
zout(1, j, nout4) = r-s
r = s1+s3
s = s2+s4
zout(2, j, nout1) = r+s
zout(2, j, nout3) = r-s
r = s1-s3
s = r2-r4
zout(2, j, nout2) = r-s
zout(2, j, nout4) = r+s
END DO
END DO
DO ia = 2, after
ias = ia-1
IF (2*ias == after) THEN
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nin4 = nin3+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
nout4 = nout3+after
DO j = 1, nfft
r1 = zin(1, j, nin1)
s1 = zin(2, j, nin1)
r = zin(1, j, nin2)
s = zin(2, j, nin2)
r2 = (r+s)*rt2i
s2 = (s-r)*rt2i
r3 = zin(2, j, nin3)
s3 = -zin(1, j, nin3)
r = zin(1, j, nin4)
s = zin(2, j, nin4)
r4 = (s-r)*rt2i
s4 = -(r+s)*rt2i
r = r1+r3
s = r2+r4
zout(1, j, nout1) = r+s
zout(1, j, nout3) = r-s
r = r1-r3
s = s2-s4
zout(1, j, nout2) = r+s
zout(1, j, nout4) = r-s
r = s1+s3
s = s2+s4
zout(2, j, nout1) = r+s
zout(2, j, nout3) = r-s
r = s1-s3
s = r2-r4
zout(2, j, nout2) = r-s
zout(2, j, nout4) = r+s
END DO
END DO
ELSE
itt = ias*before
itrig = itt+1
cr2 = trig(1, itrig)
ci2 = trig(2, itrig)
itrig = itrig+itt
cr3 = trig(1, itrig)
ci3 = trig(2, itrig)
itrig = itrig+itt
cr4 = trig(1, itrig)
ci4 = trig(2, itrig)
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nin4 = nin3+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
nout4 = nout3+after
DO j = 1, nfft
r1 = zin(1, j, nin1)
s1 = zin(2, j, nin1)
r = zin(1, j, nin2)
s = zin(2, j, nin2)
r2 = r*cr2-s*ci2
s2 = r*ci2+s*cr2
r = zin(1, j, nin3)
s = zin(2, j, nin3)
r3 = r*cr3-s*ci3
s3 = r*ci3+s*cr3
r = zin(1, j, nin4)
s = zin(2, j, nin4)
r4 = r*cr4-s*ci4
s4 = r*ci4+s*cr4
r = r1+r3
s = r2+r4
zout(1, j, nout1) = r+s
zout(1, j, nout3) = r-s
r = r1-r3
s = s2-s4
zout(1, j, nout2) = r+s
zout(1, j, nout4) = r-s
r = s1+s3
s = s2+s4
zout(2, j, nout1) = r+s
zout(2, j, nout3) = r-s
r = s1-s3
s = r2-r4
zout(2, j, nout2) = r-s
zout(2, j, nout4) = r+s
END DO
END DO
END IF
END DO
END IF
ELSE IF (now == 8) THEN
IF (isign == -1) THEN
ia = 1
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nin4 = nin3+atb
nin5 = nin4+atb
nin6 = nin5+atb
nin7 = nin6+atb
nin8 = nin7+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
nout4 = nout3+after
nout5 = nout4+after
nout6 = nout5+after
nout7 = nout6+after
nout8 = nout7+after
DO j = 1, nfft
r1 = zin(1, j, nin1)
s1 = zin(2, j, nin1)
r2 = zin(1, j, nin2)
s2 = zin(2, j, nin2)
r3 = zin(1, j, nin3)
s3 = zin(2, j, nin3)
r4 = zin(1, j, nin4)
s4 = zin(2, j, nin4)
r5 = zin(1, j, nin5)
s5 = zin(2, j, nin5)
r6 = zin(1, j, nin6)
s6 = zin(2, j, nin6)
r7 = zin(1, j, nin7)
s7 = zin(2, j, nin7)
r8 = zin(1, j, nin8)
s8 = zin(2, j, nin8)
r = r1+r5
s = r3+r7
ap = r+s
am = r-s
r = r2+r6
s = r4+r8
bp = r+s
bm = r-s
r = s1+s5
s = s3+s7
cp = r+s
cm = r-s
r = s2+s6
s = s4+s8
dbl = r+s
dm = r-s
zout(1, j, nout1) = ap+bp
zout(2, j, nout1) = cp+dbl
zout(1, j, nout5) = ap-bp
zout(2, j, nout5) = cp-dbl
zout(1, j, nout3) = am+dm
zout(2, j, nout3) = cm-bm
zout(1, j, nout7) = am-dm
zout(2, j, nout7) = cm+bm
r = r1-r5
s = s3-s7
ap = r+s
am = r-s
r = s1-s5
s = r3-r7
bp = r+s
bm = r-s
r = s4-s8
s = r2-r6
cp = r+s
cm = r-s
r = s2-s6
s = r4-r8
dbl = r+s
dm = r-s
r = (cp+dm)*rt2i
s = (-cp+dm)*rt2i
cp = (cm+dbl)*rt2i
dbl = (cm-dbl)*rt2i
zout(1, j, nout2) = ap+r
zout(2, j, nout2) = bm+s
zout(1, j, nout6) = ap-r
zout(2, j, nout6) = bm-s
zout(1, j, nout4) = am+cp
zout(2, j, nout4) = bp+dbl
zout(1, j, nout8) = am-cp
zout(2, j, nout8) = bp-dbl
END DO
END DO
ELSE
ia = 1
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nin4 = nin3+atb
nin5 = nin4+atb
nin6 = nin5+atb
nin7 = nin6+atb
nin8 = nin7+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
nout4 = nout3+after
nout5 = nout4+after
nout6 = nout5+after
nout7 = nout6+after
nout8 = nout7+after
DO j = 1, nfft
r1 = zin(1, j, nin1)
s1 = zin(2, j, nin1)
r2 = zin(1, j, nin2)
s2 = zin(2, j, nin2)
r3 = zin(1, j, nin3)
s3 = zin(2, j, nin3)
r4 = zin(1, j, nin4)
s4 = zin(2, j, nin4)
r5 = zin(1, j, nin5)
s5 = zin(2, j, nin5)
r6 = zin(1, j, nin6)
s6 = zin(2, j, nin6)
r7 = zin(1, j, nin7)
s7 = zin(2, j, nin7)
r8 = zin(1, j, nin8)
s8 = zin(2, j, nin8)
r = r1+r5
s = r3+r7
ap = r+s
am = r-s
r = r2+r6
s = r4+r8
bp = r+s
bm = r-s
r = s1+s5
s = s3+s7
cp = r+s
cm = r-s
r = s2+s6
s = s4+s8
dbl = r+s
dm = r-s
zout(1, j, nout1) = ap+bp
zout(2, j, nout1) = cp+dbl
zout(1, j, nout5) = ap-bp
zout(2, j, nout5) = cp-dbl
zout(1, j, nout3) = am-dm
zout(2, j, nout3) = cm+bm
zout(1, j, nout7) = am+dm
zout(2, j, nout7) = cm-bm
r = r1-r5
s = -s3+s7
ap = r+s
am = r-s
r = s1-s5
s = r7-r3
bp = r+s
bm = r-s
r = -s4+s8
s = r2-r6
cp = r+s
cm = r-s
r = -s2+s6
s = r4-r8
dbl = r+s
dm = r-s
r = (cp+dm)*rt2i
s = (cp-dm)*rt2i
cp = (cm+dbl)*rt2i
dbl = (-cm+dbl)*rt2i
zout(1, j, nout2) = ap+r
zout(2, j, nout2) = bm+s
zout(1, j, nout6) = ap-r
zout(2, j, nout6) = bm-s
zout(1, j, nout4) = am+cp
zout(2, j, nout4) = bp+dbl
zout(1, j, nout8) = am-cp
zout(2, j, nout8) = bp-dbl
END DO
END DO
END IF
ELSE IF (now == 3) THEN
bbs = isign*bb
ia = 1
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
DO j = 1, nfft
r1 = zin(1, j, nin1)
s1 = zin(2, j, nin1)
r2 = zin(1, j, nin2)
s2 = zin(2, j, nin2)
r3 = zin(1, j, nin3)
s3 = zin(2, j, nin3)
r = r2+r3
s = s2+s3
zout(1, j, nout1) = r+r1
zout(2, j, nout1) = s+s1
r1 = r1-0.5_dp*r
s1 = s1-0.5_dp*s
r2 = bbs*(r2-r3)
s2 = bbs*(s2-s3)
zout(1, j, nout2) = r1-s2
zout(2, j, nout2) = s1+r2
zout(1, j, nout3) = r1+s2
zout(2, j, nout3) = s1-r2
END DO
END DO
DO ia = 2, after
ias = ia-1
IF (4*ias == 3*after) THEN
IF (isign == 1) THEN
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
DO j = 1, nfft
r1 = zin(1, j, nin1)
s1 = zin(2, j, nin1)
r2 = -zin(2, j, nin2)
s2 = zin(1, j, nin2)
r3 = -zin(1, j, nin3)
s3 = -zin(2, j, nin3)
r = r2+r3
s = s2+s3
zout(1, j, nout1) = r+r1
zout(2, j, nout1) = s+s1
r1 = r1-0.5_dp*r
s1 = s1-0.5_dp*s
r2 = bbs*(r2-r3)
s2 = bbs*(s2-s3)
zout(1, j, nout2) = r1-s2
zout(2, j, nout2) = s1+r2
zout(1, j, nout3) = r1+s2
zout(2, j, nout3) = s1-r2
END DO
END DO
ELSE
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
DO j = 1, nfft
r1 = zin(1, j, nin1)
s1 = zin(2, j, nin1)
r2 = zin(2, j, nin2)
s2 = -zin(1, j, nin2)
r3 = -zin(1, j, nin3)
s3 = -zin(2, j, nin3)
r = r2+r3
s = s2+s3
zout(1, j, nout1) = r+r1
zout(2, j, nout1) = s+s1
r1 = r1-0.5_dp*r
s1 = s1-0.5_dp*s
r2 = bbs*(r2-r3)
s2 = bbs*(s2-s3)
zout(1, j, nout2) = r1-s2
zout(2, j, nout2) = s1+r2
zout(1, j, nout3) = r1+s2
zout(2, j, nout3) = s1-r2
END DO
END DO
END IF
ELSE IF (8*ias == 3*after) THEN
IF (isign == 1) THEN
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
DO j = 1, nfft
r1 = zin(1, j, nin1)
s1 = zin(2, j, nin1)
r = zin(1, j, nin2)
s = zin(2, j, nin2)
r2 = (r-s)*rt2i
s2 = (r+s)*rt2i
r3 = -zin(2, j, nin3)
s3 = zin(1, j, nin3)
r = r2+r3
s = s2+s3
zout(1, j, nout1) = r+r1
zout(2, j, nout1) = s+s1
r1 = r1-0.5_dp*r
s1 = s1-0.5_dp*s
r2 = bbs*(r2-r3)
s2 = bbs*(s2-s3)
zout(1, j, nout2) = r1-s2
zout(2, j, nout2) = s1+r2
zout(1, j, nout3) = r1+s2
zout(2, j, nout3) = s1-r2
END DO
END DO
ELSE
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
DO j = 1, nfft
r1 = zin(1, j, nin1)
s1 = zin(2, j, nin1)
r = zin(1, j, nin2)
s = zin(2, j, nin2)
r2 = (r+s)*rt2i
s2 = (-r+s)*rt2i
r3 = zin(2, j, nin3)
s3 = -zin(1, j, nin3)
r = r2+r3
s = s2+s3
zout(1, j, nout1) = r+r1
zout(2, j, nout1) = s+s1
r1 = r1-0.5_dp*r
s1 = s1-0.5_dp*s
r2 = bbs*(r2-r3)
s2 = bbs*(s2-s3)
zout(1, j, nout2) = r1-s2
zout(2, j, nout2) = s1+r2
zout(1, j, nout3) = r1+s2
zout(2, j, nout3) = s1-r2
END DO
END DO
END IF
ELSE
itt = ias*before
itrig = itt+1
cr2 = trig(1, itrig)
ci2 = trig(2, itrig)
itrig = itrig+itt
cr3 = trig(1, itrig)
ci3 = trig(2, itrig)
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
DO j = 1, nfft
r1 = zin(1, j, nin1)
s1 = zin(2, j, nin1)
r = zin(1, j, nin2)
s = zin(2, j, nin2)
r2 = r*cr2-s*ci2
s2 = r*ci2+s*cr2
r = zin(1, j, nin3)
s = zin(2, j, nin3)
r3 = r*cr3-s*ci3
s3 = r*ci3+s*cr3
r = r2+r3
s = s2+s3
zout(1, j, nout1) = r+r1
zout(2, j, nout1) = s+s1
r1 = r1-0.5_dp*r
s1 = s1-0.5_dp*s
r2 = bbs*(r2-r3)
s2 = bbs*(s2-s3)
zout(1, j, nout2) = r1-s2
zout(2, j, nout2) = s1+r2
zout(1, j, nout3) = r1+s2
zout(2, j, nout3) = s1-r2
END DO
END DO
END IF
END DO
ELSE IF (now == 5) THEN
sin2 = isign*sin2p
sin4 = isign*sin4p
ia = 1
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nin4 = nin3+atb
nin5 = nin4+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
nout4 = nout3+after
nout5 = nout4+after
DO j = 1, nfft
r1 = zin(1, j, nin1)
s1 = zin(2, j, nin1)
r2 = zin(1, j, nin2)
s2 = zin(2, j, nin2)
r3 = zin(1, j, nin3)
s3 = zin(2, j, nin3)
r4 = zin(1, j, nin4)
s4 = zin(2, j, nin4)
r5 = zin(1, j, nin5)
s5 = zin(2, j, nin5)
r25 = r2+r5
r34 = r3+r4
s25 = s2-s5
s34 = s3-s4
zout(1, j, nout1) = r1+r25+r34
r = cos2*r25+cos4*r34+r1
s = sin2*s25+sin4*s34
zout(1, j, nout2) = r-s
zout(1, j, nout5) = r+s
r = cos4*r25+cos2*r34+r1
s = sin4*s25-sin2*s34
zout(1, j, nout3) = r-s
zout(1, j, nout4) = r+s
r25 = r2-r5
r34 = r3-r4
s25 = s2+s5
s34 = s3+s4
zout(2, j, nout1) = s1+s25+s34
r = cos2*s25+cos4*s34+s1
s = sin2*r25+sin4*r34
zout(2, j, nout2) = r+s
zout(2, j, nout5) = r-s
r = cos4*s25+cos2*s34+s1
s = sin4*r25-sin2*r34
zout(2, j, nout3) = r+s
zout(2, j, nout4) = r-s
END DO
END DO
DO ia = 2, after
ias = ia-1
IF (8*ias == 5*after) THEN
IF (isign == 1) THEN
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nin4 = nin3+atb
nin5 = nin4+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
nout4 = nout3+after
nout5 = nout4+after
DO j = 1, nfft
r1 = zin(1, j, nin1)
s1 = zin(2, j, nin1)
r = zin(1, j, nin2)
s = zin(2, j, nin2)
r2 = (r-s)*rt2i
s2 = (r+s)*rt2i
r3 = -zin(2, j, nin3)
s3 = zin(1, j, nin3)
r = zin(1, j, nin4)
s = zin(2, j, nin4)
r4 = -(r+s)*rt2i
s4 = (r-s)*rt2i
r5 = -zin(1, j, nin5)
s5 = -zin(2, j, nin5)
r25 = r2+r5
r34 = r3+r4
s25 = s2-s5
s34 = s3-s4
zout(1, j, nout1) = r1+r25+r34
r = cos2*r25+cos4*r34+r1
s = sin2*s25+sin4*s34
zout(1, j, nout2) = r-s
zout(1, j, nout5) = r+s
r = cos4*r25+cos2*r34+r1
s = sin4*s25-sin2*s34
zout(1, j, nout3) = r-s
zout(1, j, nout4) = r+s
r25 = r2-r5
r34 = r3-r4
s25 = s2+s5
s34 = s3+s4
zout(2, j, nout1) = s1+s25+s34
r = cos2*s25+cos4*s34+s1
s = sin2*r25+sin4*r34
zout(2, j, nout2) = r+s
zout(2, j, nout5) = r-s
r = cos4*s25+cos2*s34+s1
s = sin4*r25-sin2*r34
zout(2, j, nout3) = r+s
zout(2, j, nout4) = r-s
END DO
END DO
ELSE
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nin4 = nin3+atb
nin5 = nin4+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
nout4 = nout3+after
nout5 = nout4+after
DO j = 1, nfft
r1 = zin(1, j, nin1)
s1 = zin(2, j, nin1)
r = zin(1, j, nin2)
s = zin(2, j, nin2)
r2 = (r+s)*rt2i
s2 = (-r+s)*rt2i
r3 = zin(2, j, nin3)
s3 = -zin(1, j, nin3)
r = zin(1, j, nin4)
s = zin(2, j, nin4)
r4 = (s-r)*rt2i
s4 = -(r+s)*rt2i
r5 = -zin(1, j, nin5)
s5 = -zin(2, j, nin5)
r25 = r2+r5
r34 = r3+r4
s25 = s2-s5
s34 = s3-s4
zout(1, j, nout1) = r1+r25+r34
r = cos2*r25+cos4*r34+r1
s = sin2*s25+sin4*s34
zout(1, j, nout2) = r-s
zout(1, j, nout5) = r+s
r = cos4*r25+cos2*r34+r1
s = sin4*s25-sin2*s34
zout(1, j, nout3) = r-s
zout(1, j, nout4) = r+s
r25 = r2-r5
r34 = r3-r4
s25 = s2+s5
s34 = s3+s4
zout(2, j, nout1) = s1+s25+s34
r = cos2*s25+cos4*s34+s1
s = sin2*r25+sin4*r34
zout(2, j, nout2) = r+s
zout(2, j, nout5) = r-s
r = cos4*s25+cos2*s34+s1
s = sin4*r25-sin2*r34
zout(2, j, nout3) = r+s
zout(2, j, nout4) = r-s
END DO
END DO
END IF
ELSE
ias = ia-1
itt = ias*before
itrig = itt+1
cr2 = trig(1, itrig)
ci2 = trig(2, itrig)
itrig = itrig+itt
cr3 = trig(1, itrig)
ci3 = trig(2, itrig)
itrig = itrig+itt
cr4 = trig(1, itrig)
ci4 = trig(2, itrig)
itrig = itrig+itt
cr5 = trig(1, itrig)
ci5 = trig(2, itrig)
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nin4 = nin3+atb
nin5 = nin4+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
nout4 = nout3+after
nout5 = nout4+after
DO j = 1, nfft
r1 = zin(1, j, nin1)
s1 = zin(2, j, nin1)
r = zin(1, j, nin2)
s = zin(2, j, nin2)
r2 = r*cr2-s*ci2
s2 = r*ci2+s*cr2
r = zin(1, j, nin3)
s = zin(2, j, nin3)
r3 = r*cr3-s*ci3
s3 = r*ci3+s*cr3
r = zin(1, j, nin4)
s = zin(2, j, nin4)
r4 = r*cr4-s*ci4
s4 = r*ci4+s*cr4
r = zin(1, j, nin5)
s = zin(2, j, nin5)
r5 = r*cr5-s*ci5
s5 = r*ci5+s*cr5
r25 = r2+r5
r34 = r3+r4
s25 = s2-s5
s34 = s3-s4
zout(1, j, nout1) = r1+r25+r34
r = cos2*r25+cos4*r34+r1
s = sin2*s25+sin4*s34
zout(1, j, nout2) = r-s
zout(1, j, nout5) = r+s
r = cos4*r25+cos2*r34+r1
s = sin4*s25-sin2*s34
zout(1, j, nout3) = r-s
zout(1, j, nout4) = r+s
r25 = r2-r5
r34 = r3-r4
s25 = s2+s5
s34 = s3+s4
zout(2, j, nout1) = s1+s25+s34
r = cos2*s25+cos4*s34+s1
s = sin2*r25+sin4*r34
zout(2, j, nout2) = r+s
zout(2, j, nout5) = r-s
r = cos4*s25+cos2*s34+s1
s = sin4*r25-sin2*r34
zout(2, j, nout3) = r+s
zout(2, j, nout4) = r-s
END DO
END DO
END IF
END DO
ELSE IF (now == 6) THEN
bbs = isign*bb
ia = 1
nin1 = ia-after
nout1 = ia-atn
DO ib = 1, before
nin1 = nin1+after
nin2 = nin1+atb
nin3 = nin2+atb
nin4 = nin3+atb
nin5 = nin4+atb
nin6 = nin5+atb
nout1 = nout1+atn
nout2 = nout1+after
nout3 = nout2+after
nout4 = nout3+after
nout5 = nout4+after
nout6 = nout5+after
DO j = 1, nfft
r2 = zin(1, j, nin3)
s2 = zin(2, j, nin3)
r3 = zin(1, j, nin5)
s3 = zin(2, j, nin5)
r = r2+r3
s = s2+s3
r1 = zin(1, j, nin1)
s1 = zin(2, j, nin1)
ur1 = r+r1
ui1 = s+s1
r1 = r1-0.5_dp*r
s1 = s1-0.5_dp*s
r = r2-r3
s = s2-s3
ur2 = r1-s*bbs
ui2 = s1+r*bbs
ur3 = r1+s*bbs
ui3 = s1-r*bbs
r2 = zin(1, j, nin6)
s2 = zin(2, j, nin6)
r3 = zin(1, j, nin2)
s3 = zin(2, j, nin2)
r = r2+r3
s = s2+s3
r1 = zin(1, j, nin4)
s1 = zin(2, j, nin4)
vr1 = r+r1
vi1 = s+s1
r1 = r1-0.5_dp*r
s1 = s1-0.5_dp*s
r = r2-r3
s = s2-s3
vr2 = r1-s*bbs
vi2 = s1+r*bbs
vr3 = r1+s*bbs
vi3 = s1-r*bbs
zout(1, j, nout1) = ur1+vr1
zout(2, j, nout1) = ui1+vi1
zout(1, j, nout5) = ur2+vr2
zout(2, j, nout5) = ui2+vi2
zout(1, j, nout3) = ur3+vr3
zout(2, j, nout3) = ui3+vi3
zout(1, j, nout4) = ur1-vr1
zout(2, j, nout4) = ui1-vi1
zout(1, j, nout2) = ur2-vr2
zout(2, j, nout2) = ui2-vi2
zout(1, j, nout6) = ur3-vr3
zout(2, j, nout6) = ui3-vi3
END DO
END DO
ELSE
CPABORT('Error fftstp')
END IF
!-----------------------------------------------------------------------------!
END SUBROUTINE fftstp
!-----------------------------------------------------------------------------!
!-----------------------------------------------------------------------------!
! Copyright (C) Stefan Goedecker, Lausanne, Switzerland, August 1, 1991
! Copyright (C) Stefan Goedecker, Cornell University, Ithaca, USA, 1994
! Copyright (C) Stefan Goedecker, MPI Stuttgart, Germany, 1999
! This file is distributed under the terms of the
! GNU General Public License version 2 (or later),
! see http://www.gnu.org/copyleft/gpl.txt .
!-----------------------------------------------------------------------------!
! S. Goedecker: Rotating a three-dimensional array in optimal
! positions for vector processing: Case study for a three-dimensional Fast
! Fourier Transform, Comp. Phys. Commun. 76, 294 (1993)
! **************************************************************************************************
!> \brief ...
!> \param n ...
!> \param trig ...
!> \param after ...
!> \param before ...
!> \param now ...
!> \param isign ...
!> \param ic ...
! **************************************************************************************************
SUBROUTINE ctrig(n, trig, after, before, now, isign, ic)
INTEGER, INTENT(IN) :: n
REAL(dp), DIMENSION(2, ctrig_length), INTENT(OUT) :: trig
INTEGER, DIMENSION(7), INTENT(OUT) :: after, before, now
INTEGER, INTENT(IN) :: isign
INTEGER, INTENT(OUT) :: ic
INTEGER, PARAMETER :: nt = 82
INTEGER, DIMENSION(7, nt), PARAMETER :: idata = RESHAPE((/3, 3, 1, 1, 1, 1, 1, 4, 4, 1, 1, 1,&
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, &
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&
, 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, &
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 &
, 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, &
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 &
, 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&
, 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 &
, 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,&
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&
, 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, &
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&
, 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,&
1, 384, 8, 4, 4, 3, 1, 1, 400