C
C     FFTE: A FAST FOURIER TRANSFORM PACKAGE
C
C     (C) COPYRIGHT SOFTWARE, 2000-2004, 2008-2011, 2020
C         ALL RIGHTS RESERVED
C                BY
C         DAISUKE TAKAHASHI
C         CENTER FOR COMPUTATIONAL SCIENCES
C         UNIVERSITY OF TSUKUBA
C         1-1-1 TENNODAI, TSUKUBA, IBARAKI 305-8577, JAPAN
C         E-MAIL: daisuke@cs.tsukuba.ac.jp
C
C
C     3-D COMPLEX FFT ROUTINE
C
C     FORTRAN77 SOURCE PROGRAM
C
C     CALL ZFFT3D(A,NX,NY,NZ,IOPT)
C
C     A(NX,NY,NZ) IS COMPLEX INPUT/OUTPUT VECTOR (COMPLEX*16)
C     NX IS THE LENGTH OF THE TRANSFORMS IN THE X-DIRECTION (INTEGER*4)
C     NY IS THE LENGTH OF THE TRANSFORMS IN THE Y-DIRECTION (INTEGER*4)
C     NZ IS THE LENGTH OF THE TRANSFORMS IN THE Z-DIRECTION (INTEGER*4)
C       ------------------------------------
C         NX = (2**IP) * (3**IQ) * (5**IR)
C         NY = (2**JP) * (3**JQ) * (5**JR)
C         NZ = (2**KP) * (3**KQ) * (5**KR)
C       ------------------------------------
C     IOPT = 0 FOR INITIALIZING THE COEFFICIENTS (INTEGER*4)
C          = -1 FOR FORWARD TRANSFORM
C          = +1 FOR INVERSE TRANSFORM
C
C     WRITTEN BY DAISUKE TAKAHASHI
C
      SUBROUTINE ZFFT3D(A,NX,NY,NZ,IOPT)
      IMPLICIT REAL*8 (A-H,O-Z)
      INCLUDE 'param.h'
      COMPLEX*16 A(*)
      COMPLEX*16 B((NDA3+NP)*NBLK),C(NDA3)
      COMPLEX*16 WX(NDA3),WY(NDA3),WZ(NDA3)
      DIMENSION LNX(3),LNY(3),LNZ(3)
      SAVE WX,WY,WZ
C
      IF (IOPT .EQ. 0) THEN
        CALL SETTBL(WX,NX)
        CALL SETTBL(WY,NY)
        CALL SETTBL(WZ,NZ)
        RETURN
      END IF
C
      IF (IOPT .EQ. 1) THEN
!$OMP PARALLEL DO
!DIR$ VECTOR ALIGNED
        DO 10 I=1,NX*NY*NZ
          A(I)=DCONJG(A(I))
   10   CONTINUE
      END IF
C
      CALL FACTOR(NX,LNX)
      CALL FACTOR(NY,LNY)
      CALL FACTOR(NZ,LNZ)
C
!$OMP PARALLEL PRIVATE(B,C)
      CALL ZFFT3D0(A,B,B,C,WX,WY,WZ,NX,NY,NZ,LNX,LNY,LNZ)
!$OMP END PARALLEL
C
      IF (IOPT .EQ. 1) THEN
        DN=1.0D0/(DBLE(NX)*DBLE(NY)*DBLE(NZ))
!$OMP PARALLEL DO
!DIR$ VECTOR ALIGNED
        DO 20 I=1,NX*NY*NZ
          A(I)=DCONJG(A(I))*DN
   20   CONTINUE
      END IF
      RETURN
      END
      SUBROUTINE ZFFT3D0(A,BY,BZ,C,WX,WY,WZ,NX,NY,NZ,LNX,LNY,LNZ)
      IMPLICIT REAL*8 (A-H,O-Z)
      INCLUDE 'param.h'
      COMPLEX*16 A(NX,NY,*)
      COMPLEX*16 BY(NY+NP,*),BZ(NZ+NP,*),C(*)
      COMPLEX*16 WX(*),WY(*),WZ(*)
      DIMENSION LNX(*),LNY(*),LNZ(*)
C
!$OMP DO PRIVATE(I,II,K,KK)
      DO 80 J=1,NY
        DO 70 II=1,NX,NBLK
          DO 30 KK=1,NZ,NBLK
            DO 20 I=II,MIN0(II+NBLK-1,NX)
!DIR$ VECTOR ALIGNED
              DO 10 K=KK,MIN0(KK+NBLK-1,NZ)
                BZ(K,I-II+1)=A(I,J,K)
   10         CONTINUE
   20       CONTINUE
   30     CONTINUE
          DO 40 I=II,MIN0(II+NBLK-1,NX)
            CALL FFT235(BZ(1,I-II+1),C,WZ,NZ,LNZ)
   40     CONTINUE
          DO 60 K=1,NZ
!DIR$ VECTOR ALIGNED
            DO 50 I=II,MIN0(II+NBLK-1,NX)
              A(I,J,K)=BZ(K,I-II+1)
   50       CONTINUE
   60     CONTINUE
   70   CONTINUE
   80 CONTINUE
!$OMP DO PRIVATE(I,II,J,JJ)
      DO 170 K=1,NZ
        DO 150 II=1,NX,NBLK
          DO 110 JJ=1,NY,NBLK
            DO 100 I=II,MIN0(II+NBLK-1,NX)
!DIR$ VECTOR ALIGNED
              DO 90 J=JJ,MIN0(JJ+NBLK-1,NY)
                BY(J,I-II+1)=A(I,J,K)
   90         CONTINUE
  100       CONTINUE
  110     CONTINUE
          DO 120 I=II,MIN0(II+NBLK-1,NX)
            CALL FFT235(BY(1,I-II+1),C,WY,NY,LNY)
  120     CONTINUE
          DO 140 J=1,NY
!DIR$ VECTOR ALIGNED
            DO 130 I=II,MIN0(II+NBLK-1,NX)
              A(I,J,K)=BY(J,I-II+1)
  130       CONTINUE
  140     CONTINUE
  150   CONTINUE
        DO 160 J=1,NY
          CALL FFT235(A(1,J,K),C,WX,NX,LNX)
  160   CONTINUE
  170 CONTINUE
      RETURN
      END