[Scipy-svn] r2210 - in trunk/Lib/fftpack: . src

scipy-svn at scipy.org scipy-svn at scipy.org
Thu Sep 21 17:19:59 CDT 2006


Author: cookedm
Date: 2006-09-21 17:19:54 -0500 (Thu, 21 Sep 2006)
New Revision: 2210

Modified:
   trunk/Lib/fftpack/setup.py
   trunk/Lib/fftpack/src/fftpack.h
   trunk/Lib/fftpack/src/zfft.c
   trunk/Lib/fftpack/src/zfftnd.c
Log:
#268: patch from jtravs for Intel MKL support for fft.

Also, replace using fft_opt_info in setup.py with explicitly checking
for each of the FFT packages that we support.


Modified: trunk/Lib/fftpack/setup.py
===================================================================
--- trunk/Lib/fftpack/setup.py	2006-09-21 16:21:58 UTC (rev 2209)
+++ trunk/Lib/fftpack/setup.py	2006-09-21 22:19:54 UTC (rev 2210)
@@ -8,8 +8,17 @@
     from numpy.distutils.misc_util import Configuration
     from numpy.distutils.system_info import get_info
     config = Configuration('fftpack',parent_package, top_path)
-    fft_opt_info = get_info('fft_opt')
 
+    djbfft_info = {}
+    mkl_info = get_info('mkl')
+    if mkl_info:
+        mkl_info.setdefault('define_macros', []).append(('SCIPY_MKL_H', None))
+        fft_opt_info = mkl_info
+    else:
+        fft_opt_info = get_info('fftw3') or get_info('fftw2') \
+                        or get_info('dfftw')
+        djbfft_info = get_info('djbfft')
+
     config.add_data_dir('tests')
 
     config.add_library('dfftpack',
@@ -27,7 +36,7 @@
     config.add_extension('convolve',
                          sources = ['convolve.pyf','src/convolve.c'],
                          libraries = ['dfftpack'],
-                         extra_info = fft_opt_info
+                         extra_info = [fft_opt_info, djbfft_info],
                          )
     return config
 

Modified: trunk/Lib/fftpack/src/fftpack.h
===================================================================
--- trunk/Lib/fftpack/src/fftpack.h	2006-09-21 16:21:58 UTC (rev 2209)
+++ trunk/Lib/fftpack/src/fftpack.h	2006-09-21 22:19:54 UTC (rev 2210)
@@ -39,6 +39,11 @@
 #include <fftr8.h>
 #endif
 
+#ifdef SCIPY_MKL_H
+#define WITH_MKL
+#include <mkl_dfti.h>
+#endif
+
 #ifdef SCIPY_FFTW3_H
 #define WITH_FFTW3
 #include <fftw3.h>

Modified: trunk/Lib/fftpack/src/zfft.c
===================================================================
--- trunk/Lib/fftpack/src/zfft.c	2006-09-21 16:21:58 UTC (rev 2209)
+++ trunk/Lib/fftpack/src/zfft.c	2006-09-21 22:19:54 UTC (rev 2210)
@@ -19,6 +19,7 @@
 #endif
 
 /**************** DJBFFT *****************************/
+#ifndef WITH_MKL
 #ifdef WITH_DJBFFT
 GEN_CACHE(zdjbfft,(int n)
 	  ,unsigned int* f;
@@ -33,8 +34,20 @@
 	   free(caches_zdjbfft[id].ptr);
 	  ,10)
 #endif
+#endif
 
-#if defined WITH_FFTW3
+/**************** INTEL MKL **************************/
+#ifdef WITH_MKL
+GEN_CACHE(zmklfft,(int n)
+	  ,DFTI_DESCRIPTOR_HANDLE desc_handle;
+	  ,(caches_zmklfft[i].n==n)
+      ,DftiCreateDescriptor(&caches_zmklfft[id].desc_handle, DFTI_DOUBLE, DFTI_COMPLEX, 1, (long)n); 
+       DftiCommitDescriptor(caches_zmklfft[id].desc_handle);
+	  ,DftiFreeDescriptor(&caches_zmklfft[id].desc_handle);
+	  ,10)
+
+/**************** FFTW3 *****************************/
+#elif defined WITH_FFTW3
 /*
  *don't cache anything
  */
@@ -69,10 +82,14 @@
 #ifdef WITH_FFTWORK
   destroy_zfftwork_caches();
 #endif
+#ifndef WITH_MKL
 #ifdef WITH_DJBFFT
   destroy_zdjbfft_caches();
 #endif
-#ifdef WITH_FFTW3
+#endif
+#ifdef WITH_MKL
+  destroy_zmklfft_caches();
+#elif defined WITH_FFTW3
 #elif defined WITH_FFTW
   destroy_zfftw_caches();
 #else
@@ -85,28 +102,36 @@
 		 int n,int direction,int howmany,int normalize) {
   int i;
   complex_double *ptr = inout;
+#ifndef WITH_MKL
 #ifdef WITH_FFTW3
   fftw_complex *ptrm = NULL;
 #endif
 #if defined(WITH_FFTW) || defined(WITH_FFTW3)
   fftw_plan plan = NULL;
+#endif
+#endif
+#if defined WITH_MKL
+  DFTI_DESCRIPTOR_HANDLE desc_handle;
 #else
   double* wsave = NULL;
 #endif
 #ifdef WITH_FFTWORK
   coef_dbl* coef = NULL;
 #endif
+#ifndef WITH_MKL
 #ifdef WITH_DJBFFT
   int j;
   complex_double *ptrc = NULL;
   unsigned int *f = NULL;
 #endif
+#endif
 #ifdef WITH_FFTWORK
   if (ispow2le2e30(n)) {
     i = get_cache_id_zfftwork(n);
     coef = caches_zfftwork[i].coef;
   } else
 #endif
+#ifndef WITH_MKL
 #ifdef WITH_DJBFFT
   switch (n) {
   case 2:;case 4:;case 8:;case 16:;case 32:;case 64:;case 128:;case 256:;
@@ -117,7 +142,10 @@
   }
   if (f==0)
 #endif
-#ifdef WITH_FFTW3
+#endif
+#ifdef WITH_MKL
+  desc_handle = caches_zmklfft[get_cache_id_zmklfft(n)].desc_handle;
+#elif defined WITH_FFTW3
 #elif defined WITH_FFTW
     plan = caches_zfftw[get_cache_id_zfftw(n,direction)].plan;
 #else
@@ -133,6 +161,7 @@
         fft_for_cplx_flt((cplx_dbl*)ptr,coef,n);
       } else
 #endif
+#ifndef WITH_MKL
 #ifdef WITH_DJBFFT
       if (f!=NULL) {
 	memcpy(ptrc,ptr,2*n*sizeof(double));
@@ -146,7 +175,10 @@
 	for (j=0;j<n;++j) *(ptr+f[j]) = *(ptrc+j);
       } else
 #endif
-#ifdef WITH_FFTW3
+#endif
+#ifdef WITH_MKL
+    DftiComputeForward(desc_handle, (double *)ptr);
+#elif defined WITH_FFTW3
 	plan = fftw_plan_dft_1d(n, (fftw_complex*)ptr, (fftw_complex*)ptr,
 				(direction>0?FFTW_FORWARD:FFTW_BACKWARD),
 				FFTW_ESTIMATE);
@@ -167,6 +199,7 @@
 	fft_bak_cplx_flt((cplx_dbl*)ptr,coef,n);
       } else
 #endif
+#ifndef WITH_MKL
 #ifdef WITH_DJBFFT
       if (f!=NULL) {
 	for (j=0;j<n;++j) *(ptrc+j) = *(ptr+f[j]);
@@ -180,7 +213,10 @@
 	memcpy(ptr,ptrc,2*n*sizeof(double));
       } else
 #endif
-#ifdef WITH_FFTW3
+#endif
+#ifdef WITH_MKL
+    DftiComputeBackward(desc_handle, (double *)ptr);
+#elif defined WITH_FFTW3
 	plan = fftw_plan_dft_1d(n, (fftw_complex*)ptr, (fftw_complex*)ptr,
 				(direction>0?FFTW_FORWARD:FFTW_BACKWARD),
 				FFTW_ESTIMATE);
@@ -197,9 +233,10 @@
   default:
     fprintf(stderr,"zfft: invalid direction=%d\n",direction);
   }
-
+  
   if (normalize) {
     ptr = inout;
+#ifndef WITH_MKL
 #ifdef WITH_DJBFFT
     if (f!=NULL) {
       for (i=0;i<howmany;++i,ptr+=n) {
@@ -213,6 +250,7 @@
       }
     } else
 #endif
+#endif
       for (i=n*howmany-1;i>=0;--i) {
 	*((double*)(ptr)) /= n;
 	*((double*)(ptr++)+1) /= n;

Modified: trunk/Lib/fftpack/src/zfftnd.c
===================================================================
--- trunk/Lib/fftpack/src/zfftnd.c	2006-09-21 16:21:58 UTC (rev 2209)
+++ trunk/Lib/fftpack/src/zfftnd.c	2006-09-21 22:19:54 UTC (rev 2210)
@@ -5,13 +5,7 @@
  */
 #include "fftpack.h"
 
-#ifdef WITH_FFTW3
-/* Don't worry about caching for fftw3 - plans take specific arrays and
- * keeping around a lot of memory for such a small speed up isn't
- * worth it.
- */
-#elif defined WITH_FFTW
-/**************** FFTW *****************************/
+#if defined(WITH_FFTW) || defined(WITH_MKL)
 static
 int equal_dims(int rank,int *dims1,int *dims2) {
   int i;
@@ -20,6 +14,43 @@
       return 0;
   return 1;
 }
+#endif
+/**************** INTEL MKL **************************/
+#ifdef WITH_MKL
+long* convert_dims(int n, int *dims)
+{
+    long * ndim;
+    int i;
+    ndim = (long*)malloc(sizeof(long)*n);
+    for(i=0;i<n;i++){ ndim[i] = (long)dims[i];}
+    return ndim;
+}
+
+GEN_CACHE(zmklfftnd,(int n,int *dims)
+	  ,DFTI_DESCRIPTOR_HANDLE desc_handle;
+	   int *dims;
+       long *ndims;
+      ,((caches_zmklfftnd[i].n==n) &&
+	   (equal_dims(n,caches_zmklfftnd[i].dims,dims)))
+	  ,caches_zmklfftnd[id].ndims = convert_dims(n, dims);
+       caches_zmklfftnd[id].n = n;
+       caches_zmklfftnd[id].dims = (int*)malloc(sizeof(int)*n);
+       memcpy(caches_zmklfftnd[id].dims,dims,sizeof(int)*n);
+       DftiCreateDescriptor(&caches_zmklfftnd[id].desc_handle, DFTI_DOUBLE, DFTI_COMPLEX, (long)n, caches_zmklfftnd[id].ndims); 
+       DftiCommitDescriptor(caches_zmklfftnd[id].desc_handle);
+	  ,DftiFreeDescriptor(&caches_zmklfftnd[id].desc_handle);
+	   free(caches_zmklfftnd[id].dims);
+       free(caches_zmklfftnd[id].ndims);
+	  ,10)
+
+/**************** FFTW3 *****************************/
+#elif defined WITH_FFTW3
+/* Don't worry about caching for fftw3 - plans take specific arrays and
+ * keeping around a lot of memory for such a small speed up isn't
+ * worth it.
+ */
+/**************** FFTW2 *****************************/
+#elif defined WITH_FFTW
 GEN_CACHE(zfftwnd,(int n,int *dims,int d,int flags)
 	  ,int direction;
 	   int *dims;
@@ -51,14 +82,16 @@
 #endif
 
 extern void destroy_zfftnd_cache(void) {
-#ifdef WITH_FFTW3
+#ifdef WITH_MKL
+  destroy_zmklfftnd_caches();
+#elif defined WITH_FFTW3
 #elif defined WITH_FFTW
   destroy_zfftwnd_caches();
 #else
   destroy_zfftnd_caches();
 #endif
 }
-#if defined WITH_FFTW || defined WITH_FFTW3
+#if defined(WITH_FFTW) || defined(WITH_FFTW3) || defined(WITH_MKL)
 #else
 static
 /*inline : disabled because MSVC6.0 fails to compile it. */
@@ -103,7 +136,9 @@
 		   int *dims,int direction,int howmany,int normalize) {
   int i,sz;
   complex_double *ptr = inout;
-#ifdef WITH_FFTW3
+#if defined WITH_MKL
+  DFTI_DESCRIPTOR_HANDLE desc_handle;
+#elif defined WITH_FFTW3
   fftw_plan plan = NULL;
 #elif defined WITH_FFTW
   fftwnd_plan plan = NULL;
@@ -116,7 +151,23 @@
   sz = 1;
   for(i=0;i<rank;++i)
     sz *= dims[i];
-#ifdef WITH_FFTW3
+#ifdef WITH_MKL
+  desc_handle = caches_zmklfftnd[get_cache_id_zmklfftnd(rank, dims)].desc_handle;
+  for (i=0;i<howmany;++i,ptr+=sz){
+    if (direction == 1){
+      DftiComputeForward(desc_handle, (double *)ptr);
+    }else if (direction == -1){
+      DftiComputeBackward(desc_handle, (double *)ptr);
+    }
+  }
+  if (normalize) {
+    ptr = inout;
+    for (i=sz*howmany-1;i>=0;--i) {
+      *((double*)(ptr)) /= sz;
+      *((double*)(ptr++)+1) /= sz;
+    }
+  }
+#elif defined WITH_FFTW3
   plan = fftw_plan_many_dft(rank,dims,howmany,
 			    (fftw_complex*)ptr,NULL,1,sz,
 			    (fftw_complex*)ptr,NULL,1,sz,



More information about the Scipy-svn mailing list