- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
1 Solution
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
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
Link Copied
2 Replies
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
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
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
Reply
Topic Options
- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page