Merge pull request #342 from dunkpa/mkl_fft

Add MKL DFTI FFT bindings.

dlib/CMakeLists.txt

@@ -84,6 +84,8 @@ if (NOT TARGET dlib)
    set (DLIB_LINK_WITH_SQLITE3_STR
    "Disable this if you don't want to link against sqlite3" )
    #set (DLIB_USE_FFTW_STR "Disable this if you don't want to link against fftw" )
+   set (DLIB_USE_MKL_FFT_STR
+   "Disable this is you don't want to use the MKL DFTI FFT implementation" )
 
    option(DLIB_ISO_CPP_ONLY ${DLIB_ISO_CPP_ONLY_STR} OFF)
    toggle_preprocessor_switch(DLIB_ISO_CPP_ONLY)
@@ -125,6 +127,7 @@ if (NOT TARGET dlib)
       option(DLIB_PNG_SUPPORT ${DLIB_PNG_SUPPORT_STR} OFF)
       option(DLIB_GIF_SUPPORT ${DLIB_GIF_SUPPORT_STR} OFF)
       #option(DLIB_USE_FFTW ${DLIB_USE_FFTW_STR} OFF)
+      option(DLIB_USE_MKL_FFT ${DLIB_USE_MKL_FFT_STR} OFF)
    else()
       option(DLIB_JPEG_SUPPORT ${DLIB_JPEG_SUPPORT_STR} ON)
       option(DLIB_LINK_WITH_SQLITE3 ${DLIB_LINK_WITH_SQLITE3_STR} ON)
@@ -134,6 +137,7 @@ if (NOT TARGET dlib)
       option(DLIB_PNG_SUPPORT ${DLIB_PNG_SUPPORT_STR} ON)
       option(DLIB_GIF_SUPPORT ${DLIB_GIF_SUPPORT_STR} ON)
       #option(DLIB_USE_FFTW ${DLIB_USE_FFTW_STR} ON)
+      option(DLIB_USE_MKL_FFT ${DLIB_USE_MKL_FFT_STR} ON)
    endif()
    toggle_preprocessor_switch(DLIB_JPEG_SUPPORT)
    toggle_preprocessor_switch(DLIB_USE_BLAS)
@@ -142,6 +146,7 @@ if (NOT TARGET dlib)
    toggle_preprocessor_switch(DLIB_PNG_SUPPORT) 
    toggle_preprocessor_switch(DLIB_GIF_SUPPORT) 
    #toggle_preprocessor_switch(DLIB_USE_FFTW)
+   toggle_preprocessor_switch(DLIB_USE_MKL_FFT)
 
 
    set(source_files 
@@ -446,8 +451,8 @@ if (NOT TARGET dlib)
       endif()
 
 
-      if (DLIB_USE_BLAS OR DLIB_USE_LAPACK)
-         # Try to find BLAS and LAPACK 
+      if (DLIB_USE_BLAS OR DLIB_USE_LAPACK OR DLIB_USE_MKL_FFT)
+          # Try to find BLAS, LAPACK and MKL
          include(cmake_utils/cmake_find_blas.txt)
 
          if (DLIB_USE_BLAS)
@@ -467,6 +472,16 @@ if (NOT TARGET dlib)
                toggle_preprocessor_switch(DLIB_USE_LAPACK)
             endif()
          endif()
+
+         if (DLIB_USE_MKL_FFT)
+            if (found_intel_mkl AND found_intel_mkl_headers)
+               set (dlib_needed_libraries ${dlib_needed_libraries} ${mkl_libraries})
+               include_directories(${mkl_include_dir})
+            else()
+               set(DLIB_USE_MKL_FFT OFF CACHE STRING ${DLIB_USE_MKL_FFT_STR} FORCE )
+               toggle_preprocessor_switch(DLIB_USE_MKL_FFT)
+            endif()
+         endif()
       endif()
 
 

dlib/cmake_utils/cmake_find_blas.txt

@@ -10,14 +10,20 @@
 #
 #  blas_found               - True if BLAS is available
 #  lapack_found             - True if LAPACK is available
+#  found_intel_mkl          - True if the Intel MKL library is available
+#  found_intel_mkl_headers  - True if Intel MKL headers are available
 #  blas_libraries           - link against these to use BLAS library 
 #  lapack_libraries         - link against these to use LAPACK library 
+#  mkl_libraries            - link against these to use the MKL library
+#  mkl_include_dir          - add to the include path to use the MKL library
 
 # setting this makes CMake allow normal looking if else statements
 SET(CMAKE_ALLOW_LOOSE_LOOP_CONSTRUCTS true)
 
 SET(blas_found 0)
 SET(lapack_found 0)
+SET(found_intel_mkl 0)
+SET(found_intel_mkl_headers 0)
 
 
 if (UNIX)
@@ -79,6 +85,12 @@ if (UNIX)
 
    include(CheckLibraryExists)
 
+    # Get mkl_include_dir
+    set(mkl_include_search_path
+        /opt/intel/mkl/include
+        /opt/intel/include
+    )
+    find_path(mkl_include_dir mkl_version.h ${mkl_include_search_path})
 
     # Search for the needed libraries from the MKL.  We will try to link against the mkl_rt
     # file first since this way avoids linking bugs in some cases.
@@ -86,6 +98,7 @@ if (UNIX)
     mark_as_advanced(  mkl_rt )
     # if we found the MKL 
     if ( mkl_rt)
+        set(mkl_libraries  ${mkl_rt} )
         set(blas_libraries  ${mkl_rt} )
         set(lapack_libraries  ${mkl_rt} )
         set(blas_found 1)
@@ -94,7 +107,6 @@ if (UNIX)
         message(STATUS "Found Intel MKL BLAS/LAPACK library")
     endif()
 
-
     if (NOT found_intel_mkl)
         # Search for the needed libraries from the MKL.  This time try looking for a different
         # set of MKL files and try to link against those.
@@ -106,6 +118,7 @@ if (UNIX)
         mark_as_advanced( mkl_intel mkl_core mkl_thread mkl_iomp mkl_pthread)
         # If we found the MKL 
         if (mkl_intel AND mkl_core AND mkl_thread AND mkl_iomp AND mkl_pthread)
+            set(mkl_libraries ${mkl_intel} ${mkl_core} ${mkl_thread} ${mkl_iomp} ${mkl_pthread})
             set(blas_libraries ${mkl_intel} ${mkl_core} ${mkl_thread} ${mkl_iomp} ${mkl_pthread})
             set(lapack_libraries ${mkl_intel} ${mkl_core} ${mkl_thread} ${mkl_iomp} ${mkl_pthread})
             set(blas_found 1)
@@ -115,6 +128,9 @@ if (UNIX)
         endif()
     endif()
 
+    if (found_intel_mkl AND mkl_include_dir)
+        set(found_intel_mkl_headers 1)
+    endif()
 
    # try to find some other LAPACK libraries if we didn't find the MKL
    set(extra_paths

dlib/config.h.in

@@ -23,3 +23,4 @@
 #cmakedefine DLIB_USE_BLAS
 #cmakedefine DLIB_USE_LAPACK
 #cmakedefine DLIB_USE_CUDA
+#cmakedefine DLIB_USE_MKL_FFT

dlib/matrix/matrix_fft.h

@@ -8,6 +8,9 @@
 #include "../hash.h"
 #include "../algs.h"
 
+#ifdef DLIB_USE_MKL_FFT
+#include <mkl_dfti.h>
+#endif
 
 // No using FFTW until it becomes thread safe!
 #if 0
@@ -614,6 +617,206 @@ namespace dlib
 
 // ----------------------------------------------------------------------------------------
 
+#ifdef DLIB_USE_MKL_FFT
+
+#define DLIB_DFTI_CHECK_STATUS(s) \
+    if((s) != 0 && !DftiErrorClass((s), DFTI_NO_ERROR)) \
+    { \
+        throw dlib::error(DftiErrorMessage((s))); \
+    }
+
+    template < long NR, long NC, typename MM, typename L >
+    matrix<std::complex<double>,NR,NC,MM,L> call_mkl_fft(
+        const matrix<std::complex<double>,NR,NC,MM,L>& data,
+        bool do_backward_fft)
+    {
+        // make sure requires clause is not broken
+        DLIB_CASSERT(is_power_of_two(data.nr()) && is_power_of_two(data.nc()),
+            "\t matrix fft(data)"
+            << "\n\t The number of rows and columns must be powers of two."
+            << "\n\t data.nr(): "<< data.nr()
+            << "\n\t data.nc(): "<< data.nc()
+            << "\n\t is_power_of_two(data.nr()): " << is_power_of_two(data.nr())
+            << "\n\t is_power_of_two(data.nc()): " << is_power_of_two(data.nc())
+            );
+
+        if (data.size() == 0)
+            return data;
+
+        DFTI_DESCRIPTOR_HANDLE h;
+        MKL_LONG status;
+
+        if (data.nr() == 1 || data.nc() == 1)
+        {
+            status = DftiCreateDescriptor(&h, DFTI_DOUBLE, DFTI_COMPLEX, 1, data.size());
+            DLIB_DFTI_CHECK_STATUS(status);
+        }
+        else
+        {
+            MKL_LONG size[2];
+            size[0] = data.nr();
+            size[1] = data.nc();
+
+            status = DftiCreateDescriptor(&h, DFTI_DOUBLE, DFTI_COMPLEX, 2, size);
+            DLIB_DFTI_CHECK_STATUS(status);
+
+            MKL_LONG strides[3];
+            strides[0] = 0;
+            strides[1] = size[1];
+            strides[2] = 1;
+
+            status = DftiSetValue(h, DFTI_INPUT_STRIDES, strides);
+            DLIB_DFTI_CHECK_STATUS(status);
+            status = DftiSetValue(h, DFTI_OUTPUT_STRIDES, strides);
+            DLIB_DFTI_CHECK_STATUS(status);
+        }
+
+        status = DftiSetValue(h, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
+        DLIB_DFTI_CHECK_STATUS(status);
+
+        // Unless we use sequential mode, the fft results are not correct.
+        status = DftiSetValue(h, DFTI_THREAD_LIMIT, 1);
+        DLIB_DFTI_CHECK_STATUS(status);
+
+        status = DftiCommitDescriptor(h);
+        DLIB_DFTI_CHECK_STATUS(status);
+
+        matrix<std::complex<double>,NR,NC,MM,L> out(data.nr(), data.nc());
+
+        if (do_backward_fft)
+            status = DftiComputeBackward(h, (void *)(&data(0, 0)), &out(0,0));
+        else
+            status = DftiComputeForward(h, (void *)(&data(0, 0)), &out(0,0));
+        DLIB_DFTI_CHECK_STATUS(status);
+
+        status = DftiFreeDescriptor(&h);
+        DLIB_DFTI_CHECK_STATUS(status);
+
+        return out;
+    }
+
+    template < long NR, long NC, typename MM, typename L >
+    void call_mkl_fft_inplace(
+        matrix<std::complex<double>,NR,NC,MM,L>& data,
+        bool do_backward_fft
+    )
+    {
+        // make sure requires clause is not broken
+        DLIB_CASSERT(is_power_of_two(data.nr()) && is_power_of_two(data.nc()),
+            "\t void ifft_inplace(data)"
+            << "\n\t The number of rows and columns must be powers of two."
+            << "\n\t data.nr(): "<< data.nr()
+            << "\n\t data.nc(): "<< data.nc()
+            << "\n\t is_power_of_two(data.nr()): " << is_power_of_two(data.nr())
+            << "\n\t is_power_of_two(data.nc()): " << is_power_of_two(data.nc())
+            );
+
+        if (data.size() == 0)
+            return;
+
+        DFTI_DESCRIPTOR_HANDLE h;
+        MKL_LONG status;
+
+        if (data.nr() == 1 || data.nc() == 1)
+        {
+            status = DftiCreateDescriptor(&h, DFTI_DOUBLE, DFTI_COMPLEX, 1, data.size());
+            DLIB_DFTI_CHECK_STATUS(status);
+        }
+        else
+        {
+            MKL_LONG size[2];
+            size[0] = data.nr();
+            size[1] = data.nc();
+
+            status = DftiCreateDescriptor(&h, DFTI_DOUBLE, DFTI_COMPLEX, 2, size);
+            DLIB_DFTI_CHECK_STATUS(status);
+
+            MKL_LONG strides[3];
+            strides[0] = 0;
+            strides[1] = size[1];
+            strides[2] = 1;
+
+            status = DftiSetValue(h, DFTI_INPUT_STRIDES, strides);
+            DLIB_DFTI_CHECK_STATUS(status);
+        }
+
+        // Unless we use sequential mode, the fft results are not correct.
+        status = DftiSetValue(h, DFTI_THREAD_LIMIT, 1);
+        DLIB_DFTI_CHECK_STATUS(status);
+
+        status = DftiCommitDescriptor(h);
+        DLIB_DFTI_CHECK_STATUS(status);
+
+        if (do_backward_fft)
+            status = DftiComputeBackward(h, &data(0, 0));
+        else
+            status = DftiComputeForward(h, &data(0, 0));
+        DLIB_DFTI_CHECK_STATUS(status);
+
+        status = DftiFreeDescriptor(&h);
+        DLIB_DFTI_CHECK_STATUS(status);
+
+        return;
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    // Call the MKL DFTI implementation in these cases 
+
+    inline matrix<std::complex<double>,0,1> fft (const matrix<std::complex<double>,0,1>& data)
+    {
+        return call_mkl_fft(data, false);
+    }
+    inline matrix<std::complex<double>,0,1> ifft(const matrix<std::complex<double>,0,1>& data)
+    {
+        return call_mkl_fft(data, true) / data.size();
+    }
+    inline matrix<std::complex<double>,1,0> fft (const matrix<std::complex<double>,1,0>& data)
+    {
+        return call_mkl_fft(data, false);
+    }
+    inline matrix<std::complex<double>,1,0> ifft(const matrix<std::complex<double>,1,0>& data)
+    {
+        return call_mkl_fft(data, true) / data.size();
+    }
+    inline matrix<std::complex<double> > fft (const matrix<std::complex<double> >& data)
+    {
+        return call_mkl_fft(data, false);
+    }
+    inline matrix<std::complex<double> > ifft(const matrix<std::complex<double> >& data)
+    {
+        return call_mkl_fft(data, true) / data.size();
+    }
+
+    inline void fft_inplace (matrix<std::complex<double>,0,1>& data)
+    {
+        call_mkl_fft_inplace(data, false);
+    }
+    inline void ifft_inplace(matrix<std::complex<double>,0,1>& data)
+    {
+        call_mkl_fft_inplace(data, true);
+    }
+    inline void fft_inplace (matrix<std::complex<double>,1,0>& data)
+    {
+        call_mkl_fft_inplace(data, false);
+    }
+    inline void ifft_inplace(matrix<std::complex<double>,1,0>& data)
+    {
+        call_mkl_fft_inplace(data, true);
+    }
+
+    inline void fft_inplace (matrix<std::complex<double> >& data)
+    {
+        call_mkl_fft_inplace(data, false);
+    }
+    inline void ifft_inplace(matrix<std::complex<double> >& data)
+    {
+        call_mkl_fft_inplace(data, true);
+    }
+
+#endif // DLIB_USE_MKL_FFT
+
+// ----------------------------------------------------------------------------------------
 }
 
 #endif // DLIB_FFt_Hh_