topic Fortran 3-D FFT real-to-complex and complex-to-real examples in Intel® oneAPI Math Kernel Library

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

sbarbot — Thu, 06 May 2010 23:56:31 GMT

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

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

Dmitry_B_Intel — Fri, 07 May 2010 03:28:07 GMT

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

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

sbarbot — Fri, 07 May 2010 18:20:29 GMT

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