C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
C C3FFT(MFFT) <---> MKL(DFT)
C IFFT = 1:INVERSE FFT
C IFFT = -1:FFT
C SUBROUTINE C3FFT_MKL(ZRB,KFT1,KFT1_1,KFT2,KFT3,IWL,IWM,IWN
SUBROUTINE C3FFT(ZRB,KFT1,KFT1_1,KFT2,KFT3,IWL,IWM,IWN
& ,IDUM1,IFFT,IDUM2,IWORK,IERR)
IMPLICIT REAL(A-H,O-Y)
IMPLICIT COMPLEX(Z)
INCLUDE 'PACVPP'
C-----ARRAYS FOR MFFT--------------------------------------------
DIMENSION IWL(8*IFX2+28),IWM(8*IFY2+28),IWN(8*IFZ2+28)
& ,IWORK(2*IFX2),IP1(3),IP2(3),IP3(3),IOPT(2)
C
C ZWK(NWK): NWK >= MAX((NX + NWKW)*NY*NZ, NX*(NY + NWKW)*NZ)
C
C PARAMETER(ISX=0,ISY=0,ISZ=0,IFX=KNX+ISX,IFY=KNY+ISY,IFZ=KNZ+ISZ)
C PARAMETER(IFX2=2*IFX-1,IFY2=2*IFY-2,IFZ2=2*IFZ-2)
C PARAMETER(KFFT=IFX2*IFY2*IFZ2)
C
PARAMETER(KFFT2=IFX2*(IFY2+1)*IFZ2)
C KFFT
DIMENSION ZTB1(IFX2-1),ZTB2(IFY2),ZTB3(IFZ2),ZWK(KFFT2)
DIMENSION ZRB(IFX2,IFY2,IFZ2)
C
KFT1 = IFX2
KFT2 = IFY2
KFT3 = IFZ2
IERR = 0
C
IF (IFFT.NE.0) THEN
C
IERR = 0
KSUM=KFT1*KFT2*KFT3
KVOL=(KFT1-1)*KFT2*KFT3
C
C SUBROUTINE C3FFT(ZRB,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN
C & ,IDUM1,IFFT,IDUM2,IWORK,IERR)
C IFFT=1:INVERSE FFT,IFFT=-1:FFT
C SUBROUTINE C3FFT_MKL(X_IN,KFT1,KFT1-1,KFT2,KFT3,IERR)
CALL C3FFT_MKL(ZRB,KFT1,KFT1-1,KFT2,KFT3,IFFT,IERR)
C
IF (IERR.NE.0) THEN
WRITE (6,*) ' C3FFT(SET UP)] IERR = ',IERR
STOP
END IF
C
END IF
C
RETURN
END
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SUBROUTINE C3FFT_MKL(ZRB,KFT1,KFT1_1,KFT2,KFT3,IFFT,IERR)
C
C IntelのMKLにあった参考プログラムを元にして、筆者使用の
C プログラムで動くように新たに作り直したルーチンである（無保証）。
C [1/7, 2009]
C
USE MKL_DFTI
C
INCLUDE 'PACVPP'
C
INTEGER M, N, K
C
INTEGER IFFT, IERR, KFT1,KFT2,KFT3,KFT1_1
C
COMPLEX(8) :: X_IN_3D(IFX2-1,IFY2,IFZ2)
COMPLEX(8) :: X_IN ((IFX2-1)*IFY2*IFZ2)
C
COMPLEX(8) :: ZRB(IFX2,IFY2,IFZ2)
C
EQUIVALENCE (X_IN, X_IN_3D)
C
TYPE(DFTI_DESCRIPTOR), POINTER :: Desc_Handle
INTEGER STATUS
REAL(8) SCALE
INTEGER LENGTHS(3)
C
INTEGER I, J, L
C
C m <-- KFT1_1
C n <-- KFT2
C k <-- KFT3
C
C IFFT = 1:INVERSE FFT
C IFFT = -1:FFT
C WRITE (6,*) "IFFT,IERR = ",IFFT,IERR
C KFT1_1 = KFT1 - 1
M = KFT1_1
N = KFT2
K = KFT3
C
DO 111 L=1,IFZ2
DO 112 J=1,IFY2
DO 113 I=1,IFX2-1
C DO 113 I=1,IFX2
X_IN_3D(I,J,L)=ZRB(I,J,L)
113 CONTINUE
112 CONTINUE
111 CONTINUE
C
LENGTHS(1) = M
LENGTHS(2) = N
LENGTHS(3) = K
C
C DftiCreateDescriptor
C
STATUS = DftiCreateDescriptor( Desc_Handle,
& DFTI_DOUBLE, DFTI_COMPLEX, 3, LENGTHS)
C
IF (IFFT .EQ. -1) THEN
C
C DftiCommitDescriptor
C
STATUS = DftiCommitDescriptor( Desc_Handle )
C
C DftiComputeForward
C
STATUS = DftiComputeForward( Desc_Handle, X_IN)
C
ELSE IF (IFFT .EQ. 1) THEN
C
C SET SCALE = 1.0 (1/6, 2008) = NUMBER FOR BACKWARD TRANSFORM
C DftiSetValue
C
C SCALE = 1.0/REAL(M*N*K, KIND=8)
SCALE = 1.0
C
STATUS = DftiSetValue(Desc_Handle, DFTI_BACKWARD_SCALE, SCALE)
C
C DftiCommitDescriptor
C
STATUS = DftiCommitDescriptor( Desc_Handle )
C
C DftiComputeBackward
C
STATUS = DftiComputeBackward( Desc_Handle, X_IN)
C
END IF
C
DO 211 L=1,IFZ2
DO 212 J=1,IFY2
DO 213 I=1,IFX2-1
C DO 113 I=1,IFX2
ZRB(I,J,L) = X_IN_3D(I,J,L)
213 CONTINUE
212 CONTINUE
211 CONTINUE
C
100 CONTINUE
STATUS = DftiFreeDescriptor(Desc_Handle)
C
101 CONTINUE
END