include "mkl_dfti.f90" program main use omp_lib use mkl_dfti implicit none !fundamental constant integer, parameter :: DP = kind(0.0) real(DP), parameter :: PI = acos(-1.0) complex(DP), parameter :: USO = (0.0, 1.0) !contants integer, parameter :: ThreadNumber = 36, Nx = 512, ID_CPUTime = 40, ID_debug = 90 !Variables real(DP), save, allocatable, dimension(:, :, :, :) :: R_Velo complex(DP), save, dimension(:), allocatable :: Tmp_K_Field real(DP), save, dimension(:), allocatable :: Tmp_R_Field character, save :: FileName*48 !For cpu time integer, save, dimension(8) :: jikoku real(DP), save :: time_start, time_end !For FFT integer :: stat TYPE(DFTI_DESCRIPTOR), pointer, save :: Handle1, Handle2 integer, dimension(3) :: Array_D = [Nx/2+1, Nx, Nx] integer, dimension(4) :: Array_Stride_R = [0, 1, Nx, Nx*Nx], & & Array_Stride_K = [0, 1, Nx/2+1, (Nx/2+1)*Nx] !For loops integer :: ix, iy, iz, iDire, NTime !Initialize FFT stat = DftiCreateDescriptor(Handle1, DFTI_SINGLE, DFTI_REAL, 3, Array_D) stat = DftiSetValue(Handle1, DFTI_PLACEMENT, DFTI_NOT_INPLACE) stat = DftiSetValue(Handle1, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX) stat = DftiSetValue(Handle1, DFTI_INPUT_STRIDES, Array_Stride_R) stat = DftiSetValue(Handle1, DFTI_OUTPUT_STRIDES, Array_Stride_K) stat = DftiSetValue(Handle1, DFTI_FORWARD_SCALE, 1.0/real(Nx*Nx*Nx, kind=DP)) stat = DftiSetValue(Handle1, DFTI_BACKWARD_SCALE, 1.0) stat = DftiCommitDescriptor(Handle1) stat = DftiCreateDescriptor(Handle2, DFTI_SINGLE, DFTI_REAL, 3, Array_D) stat = DftiSetValue(Handle2, DFTI_PLACEMENT, DFTI_NOT_INPLACE) stat = DftiSetValue(Handle2, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX) stat = DftiSetValue(Handle2, DFTI_INPUT_STRIDES, Array_Stride_K) stat = DftiSetValue(Handle2, DFTI_OUTPUT_STRIDES, Array_Stride_R) stat = DftiSetValue(Handle2, DFTI_FORWARD_SCALE, 1.0/real(Nx*Nx*Nx, kind=DP)) stat = DftiSetValue(Handle2, DFTI_BACKWARD_SCALE, 1.0) stat = DftiCommitDescriptor(Handle2) !Set Thread Number call OMP_set_num_threads( ThreadNumber ) !Save Starting time write(FileName, "('CPUTime.dat')") call cpu_time(time_start) call date_and_time(values=jikoku) open(ID_CPUTime, file = FileName) write(ID_CPUTime, *) 'Start' write(ID_CPUTime, *) jikoku(3), 'th', jikoku(5), 'h', jikoku(6), 'min', jikoku(7), 'sec' !Initialize varibales allocate(Tmp_K_Field((Nx/2 + 1) * Nx * Nx)) allocate(R_Velo(0:Nx-1, 0:Nx-1, 0:Nx-1, 3)) allocate(Tmp_R_Field(Nx * Nx * Nx)) Tmp_K_Field(:) = 0.0 !Loop using FFT do Ntime = 1, 10 do iDire = 1, 3 stat = DftiComputeBackward(Handle2, Tmp_K_Field, Tmp_R_Field) do iz = 0, Nx-1; do iy = 0, Nx-1; do ix = 0, Nx-1 R_Velo(ix, iy, iz, iDire) = Tmp_R_Field(ix + Nx*(iy + iz * Nx) + 1) enddo; enddo; enddo enddo stat = DftiComputeForward(Handle1, Tmp_R_Field, Tmp_K_Field) enddo !Save ending time call date_and_time(values=jikoku) write(ID_CPUTime, *) 'End' write(ID_CPUTime, *) jikoku(3), 'th', jikoku(5), 'h', jikoku(6), 'min', jikoku(7), 'sec' call cpu_time(time_end) write(ID_CPUTime, *) 'CPU Time' write(ID_CPUTime, *) time_end - time_start close(ID_CPUTime) deallocate(Tmp_R_Field) deallocate(Tmp_K_Field) deallocate(R_Velo) end program main