Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.
6956 Discussions

Fortran 3-D FFT real-to-complex and complex-to-real examples

sbarbot
Beginner
412 Views
Hi,
I have tried to put together a real-to-complex, complex-to-real Fourier transform with ifort and the MKL library. I found that often the documentation was spotty: for example, there's no example of a back Fourier transform in the documentation or in the forum. I essentially want to use the CCE storage format: there is an option to specify CCE storage, but it can be forced with the strides system.
I provide an example of 3-D real-to-complex, and complex-to-real Fourier transforms in Fortran. I obtain similar results with FFTW3. The example is functional, but some updates of the descriptor might be redundant:
with the definitions
#define WRITE_MKL_DEBUG_INFO(i) IF (i .NE. 0) THEN; IF (.NOT. DftiErrorClass(i,DFTI_NO_ERROR)) THEN; WRITE_DEBUG_INFO; WRITE (0,*) DftiErrorMessage(i); STOP 1; END IF; END IF
#define FFT_FORWARD -1
#define FFT_INVERSE=1
the FFT subroutine is
SUBROUTINE fft3(data,sx1,sx2,sx3,dx1,dx2,dx3,direction)
REAL*4, DIMENSION(0:*), INTENT(INOUT) :: data
REAL*8, INTENT(IN) :: dx1,dx2,dx3
INTEGER, INTENT(IN) :: sx1,sx2,sx3,direction
INTEGER :: iret,size(3),rstrides(4),cstrides(4)
TYPE(DFTI_DESCRIPTOR), POINTER :: desc
REAL*4 :: scale
rstrides=(/ 0,1,sx1+2,(sx1+2)*sx2 /)
cstrides=(/ 0,1,sx1/2+1,(sx1/2+1)*sx2 /)
size=(/ sx1,sx2,sx3 /)
iret=DftiCreateDescriptor(desc,DFTI_SINGLE,DFTI_REAL,3,size)
iret=DftiSetValue(desc,DFTI_FORWARD_DOMAIN,DFTI_REAL)
iret=DftiSetValue(desc,DFTI_CONJUGATE_EVEN_STORAGE,DFTI_COMPLEX_COMPLEX)
iret=DftiSetValue(desc,DFTI_PACKED_FORMAT,DFTI_CCE_FORMAT)
WRITE_MKL_DEBUG_INFO(iret)
IF (FFT_FORWARD == direction) THEN
scale=dx1*dx2*dx3
iret=DftiSetValue(desc,DFTI_FORWARD_SCALE,scale)
iret=DftiSetValue(desc,DFTI_INPUT_STRIDES,rstrides);
iret=DftiSetValue(desc,DFTI_OUTPUT_STRIDES,cstrides);
iret=DftiCommitDescriptor(desc)
iret=DftiComputeForward(desc,data)
ELSE
scale=1._4/(sx1*dx1*sx2*dx2*sx3*dx3)
iret=DftiSetValue(desc,DFTI_BACKWARD_SCALE,scale)
iret=DftiSetValue(desc,DFTI_INPUT_STRIDES,cstrides);
iret=DftiSetValue(desc,DFTI_OUTPUT_STRIDES,rstrides);
iret=DftiCommitDescriptor(desc)
iret=DftiComputeBackward(desc,data)
END IF
iret=DftiFreeDescriptor(desc)
WRITE_MKL_DEBUG_INFO(iret)
END SUBROUTINE fft3
and an exert of the makefile:
LIBPATH=-L/sw/lib -L/opt/intel/Compiler/11.1/084/Frameworks/mkl/lib/em64t/ -L/opt/intel/Compiler/11.1/084/lib/
INCPATH=-I/sw/include -I/opt/intel/Compiler/11.1/084/Frameworks/mkl/include

LIBS=-lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -openmp -lpthread

FC=ifort -openmp

OBJ = *.o
%.o : %.f90
$(FC) $(FFLAGS) -c $^

program: $(OBJ)
$(FC) $(FFLAGS) -o $@ $^ $(LIBPATH) $(LIBS)
Let me know if the code can be trimmed a bit.
Thanks,
Sylvain Barbot
0 Kudos
1 Solution
Dmitry_B_Intel
Employee
412 Views
Hi Sylvain,

It looks likeodd sx1 will be mishandled.The following changewould do processing them:

- rstrides=(/ 0,1,sx1+2,(sx1+2)*sx2 /)
+ rstrides=(/ 0,1,(sx1/2+1)*2,((sx1/2+1)*2)*sx2 /)

The line setting DFTI_FORWARD_DOMAIN is redundant, and in fact is an error because this parameter is read-only (it is already set by DftiCreateDescriptor).

The line setting DFTI_PACKED_FORMAT is redundant after DFTI_CONJUGATE_EVEN_STORAGE is set to DFTI_COMPLEX_COMPLEX.

Thanks
Dima

View solution in original post

0 Kudos
2 Replies
Dmitry_B_Intel
Employee
413 Views
Hi Sylvain,

It looks likeodd sx1 will be mishandled.The following changewould do processing them:

- rstrides=(/ 0,1,sx1+2,(sx1+2)*sx2 /)
+ rstrides=(/ 0,1,(sx1/2+1)*2,((sx1/2+1)*2)*sx2 /)

The line setting DFTI_FORWARD_DOMAIN is redundant, and in fact is an error because this parameter is read-only (it is already set by DftiCreateDescriptor).

The line setting DFTI_PACKED_FORMAT is redundant after DFTI_CONJUGATE_EVEN_STORAGE is set to DFTI_COMPLEX_COMPLEX.

Thanks
Dima
0 Kudos
sbarbot
Beginner
412 Views
Hi Dima,
Thanks for your answer. I updated the function as you suggested. Everything works well.
For completeness, I want to mention that one needs to set the environment variableDYLD_LIBRARY_PATH to
export DYLD_LIBRARY_PATH="/opt/intel/Compiler/11.1/084/lib:$DYLD_LIBRARY_PATH"
in order to load the libraries successfully at runtime.I am running the same program under two implementations of the FFT. One is FFTW3 with 16 threads and another is Intel MKL FFT with 16 threads.
I find a factor of x4 improvement of the Intel FFT over the FFTW3.
Pretty impressive.
I show the updated code for reference. Thank you again. Sylvain
SUBROUTINE fft3(data,sx1,sx2,sx3,dx1,dx2,dx3,direction)
REAL*4, DIMENSION(0:*), INTENT(INOUT) :: data
REAL*8, INTENT(IN) :: dx1,dx2,dx3
INTEGER, INTENT(IN) :: sx1,sx2,sx3,direction
INTEGER :: iret,size(3),rstrides(4),cstrides(4)
TYPE(DFTI_DESCRIPTOR), POINTER :: desc
REAL*4 :: scale
rstrides=(/ 0,1,(sx1/2+1)*2,(sx1/2+1)*2*sx2 /)
cstrides=(/ 0,1,sx1/2+1,(sx1/2+1)*sx2 /)
size=(/ sx1,sx2,sx3 /)
iret=DftiCreateDescriptor(desc,DFTI_SINGLE,DFTI_REAL,3,size)
iret=DftiSetValue(desc,DFTI_CONJUGATE_EVEN_STORAGE,DFTI_COMPLEX_COMPLEX)
WRITE_MKL_DEBUG_INFO(iret)
IF (FFT_FORWARD == direction) THEN
scale=dx1*dx2*dx3
iret=DftiSetValue(desc,DFTI_FORWARD_SCALE,scale)
iret=DftiSetValue(desc,DFTI_INPUT_STRIDES,rstrides);
iret=DftiSetValue(desc,DFTI_OUTPUT_STRIDES,cstrides);
iret=DftiCommitDescriptor(desc)
iret=DftiComputeForward(desc,data)
ELSE
scale=1._4/(sx1*dx1*sx2*dx2*sx3*dx3)
iret=DftiSetValue(desc,DFTI_BACKWARD_SCALE,scale)
iret=DftiSetValue(desc,DFTI_INPUT_STRIDES,cstrides);
iret=DftiSetValue(desc,DFTI_OUTPUT_STRIDES,rstrides);
iret=DftiCommitDescriptor(desc)
iret=DftiComputeBackward(desc,data)
END IF
iret=DftiFreeDescriptor(desc)
WRITE_MKL_DEBUG_INFO(iret)
END SUBROUTINE fft3
0 Kudos
Reply