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

scipy-svn@scip... scipy-svn@scip...
Wed Aug 8 01:40:21 CDT 2007


Author: cdavid
Date: 2007-08-08 01:39:47 -0500 (Wed, 08 Aug 2007)
New Revision: 3223

Added:
   trunk/Lib/fftpack/src/zfftnd_fftpack.c
   trunk/Lib/fftpack/src/zfftnd_fftw.c
   trunk/Lib/fftpack/src/zfftnd_fftw3.c
   trunk/Lib/fftpack/src/zfftnd_mkl.c
Modified:
   trunk/Lib/fftpack/setup.py
   trunk/Lib/fftpack/src/zfftnd.c
Log:
Cleaning fft sources for Multi dimensional fft

Modified: trunk/Lib/fftpack/setup.py
===================================================================
--- trunk/Lib/fftpack/setup.py	2007-08-07 23:55:56 UTC (rev 3222)
+++ trunk/Lib/fftpack/setup.py	2007-08-08 06:39:47 UTC (rev 3223)
@@ -33,7 +33,10 @@
         depends=['src/zfft_djbfft.c', 'src/zfft_fftpack.c', 'src/zfft_fftw.c',
             'src/zfft_fftw3.c', 'src/zfft_mkl.c', 
             'src/drfft_djbfft.c', 'src/drfft_fftpack.c',
-            'src/drfft_fftw3.c', 'src/drfft_fftw.c'],
+            'src/drfft_fftw3.c', 'src/drfft_fftw.c',
+            'src/zfftnd_fftpack.c', 'src/zfftnd_fftw.c',
+            'src/zfftnd_fftw3.c', 'src/zfftnd_mkl.c',
+            ],
     )
 
     config.add_extension('convolve',

Modified: trunk/Lib/fftpack/src/zfftnd.c
===================================================================
--- trunk/Lib/fftpack/src/zfftnd.c	2007-08-07 23:55:56 UTC (rev 3222)
+++ trunk/Lib/fftpack/src/zfftnd.c	2007-08-08 06:39:47 UTC (rev 3223)
@@ -5,6 +5,20 @@
  */
 #include "fftpack.h"
 
+/* The following macro convert private backend specific function to the public
+ * functions exported by the module  */
+#define GEN_PUBLIC_API(name) \
+void destroy_zfftnd_cache(void)\
+{\
+        destroy_zfftnd_##name##_caches();\
+}\
+\
+void zfftnd(complex_double * inout, int rank,\
+		           int *dims, int direction, int howmany, int normalize)\
+{\
+        zfftnd_##name(inout, rank, dims, direction, howmany, normalize);\
+}
+
 #if defined(WITH_FFTW) || defined(WITH_MKL)
 static
 int equal_dims(int rank,int *dims1,int *dims2) {
@@ -15,6 +29,22 @@
   return 1;
 }
 #endif
+
+#ifdef WITH_FFTW3
+    #include "zfftnd_fftw3.c"
+    GEN_PUBLIC_API(fftw3)
+#elif defined WITH_FFTW
+    #include "zfftnd_fftw.c"
+    GEN_PUBLIC_API(fftw)
+#elif defined WITH_MKL
+    #include "zfftnd_mkl.c"
+    GEN_PUBLIC_API(mkl)
+#else /* Use fftpack by default */
+    #include "zfftnd_fftpack.c"
+    GEN_PUBLIC_API(fftpack)
+#endif
+
+#if 0
 /**************** INTEL MKL **************************/
 #ifdef WITH_MKL
 long* convert_dims(int n, int *dims)
@@ -93,43 +123,6 @@
 }
 #if defined(WITH_FFTW) || defined(WITH_FFTW3) || defined(WITH_MKL)
 #else
-static
-/*inline : disabled because MSVC6.0 fails to compile it. */
-int next_comb(int *ia,int *da,int m) {
-  while (m>=0 && ia[m]==da[m]) ia[m--] = 0;
-  if (m<0) return 0;
-  ia[m]++;
-  return 1;
-}
-static
-void flatten(complex_double *dest,complex_double *src,
-	     int rank,int strides_axis,int dims_axis,int unflat,int *tmp) {
-  int *new_strides = tmp+rank;
-  int *new_dims = tmp+2*rank;
-  int *ia = tmp+3*rank;
-  int rm1=rank-1,rm2=rank-2;
-  int i,j,k;
-  for (i=0;i<rm2;++i) ia[i]=0;ia[rm2] = -1;
-  j = 0;
-  if (unflat)
-    while (next_comb(ia,new_dims,rm2)) {
-      k = 0;
-      for(i=0;i<rm1;++i)
-	k += ia[i] * new_strides[i];
-      for(i=0;i<dims_axis;++i)
-	*(dest+k+i*strides_axis) = *(src+j++);
-    }
-  else
-    while (next_comb(ia,new_dims,rm2)) {
-      k = 0;
-      for(i=0;i<rm1;++i)
-	k += ia[i] * new_strides[i];
-      for(i=0;i<dims_axis;++i)
-	*(dest+j++) = *(src+k+i*strides_axis);
-    }
-}
-extern void zfft(complex_double *inout,
-		 int n,int direction,int howmany,int normalize);
 #endif
 /**************** ZFFTND function **********************/
 extern void zfftnd(complex_double *inout,int rank,
@@ -223,3 +216,4 @@
   }
 #endif
 }
+#endif

Added: trunk/Lib/fftpack/src/zfftnd_fftpack.c
===================================================================
--- trunk/Lib/fftpack/src/zfftnd_fftpack.c	2007-08-07 23:55:56 UTC (rev 3222)
+++ trunk/Lib/fftpack/src/zfftnd_fftpack.c	2007-08-08 06:39:47 UTC (rev 3223)
@@ -0,0 +1,121 @@
+/*
+ * fftpack backend for multi dimensional fft
+ *
+ * Original code by Pearu Peaterson
+ *
+ * Last Change: Wed Aug 08 02:00 PM 2007 J
+ */
+
+GEN_CACHE(zfftnd_fftpack, (int n, int rank)
+	  , complex_double * ptr; int *iptr; int rank;
+	  , ((caches_zfftnd_fftpack[i].n == n)
+	     && (caches_zfftnd_fftpack[i].rank == rank))
+	  , caches_zfftnd_fftpack[id].n = n;
+	  caches_zfftnd_fftpack[id].ptr =
+	  (complex_double *) malloc(2 * sizeof(double) * n);
+	  caches_zfftnd_fftpack[id].iptr =
+	  (int *) malloc(4 * rank * sizeof(int));
+	  ,
+	  free(caches_zfftnd_fftpack[id].ptr);
+	  free(caches_zfftnd_fftpack[id].iptr);
+	  , 10)
+
+static
+/*inline : disabled because MSVC6.0 fails to compile it. */
+int next_comb(int *ia, int *da, int m)
+{
+    while (m >= 0 && ia[m] == da[m]) {
+        ia[m--] = 0;
+    }
+    if (m < 0) {
+        return 0;
+    }
+    ia[m]++;
+    return 1;
+}
+
+static
+void flatten(complex_double * dest, complex_double * src,
+	     int rank, int strides_axis, int dims_axis, int unflat,
+	     int *tmp)
+{
+    int *new_strides = tmp + rank;
+    int *new_dims = tmp + 2 * rank;
+    int *ia = tmp + 3 * rank;
+    int rm1 = rank - 1, rm2 = rank - 2;
+    int i, j, k;
+    for (i = 0; i < rm2; ++i)
+	ia[i] = 0;
+    ia[rm2] = -1;
+    j = 0;
+    if (unflat) {
+        while (next_comb(ia, new_dims, rm2)) {
+            k = 0;
+            for (i = 0; i < rm1; ++i) {
+                k += ia[i] * new_strides[i];
+            }
+            for (i = 0; i < dims_axis; ++i) {
+                *(dest + k + i * strides_axis) = *(src + j++);
+            }
+        }
+    } else {
+        while (next_comb(ia, new_dims, rm2)) {
+            k = 0;
+            for (i = 0; i < rm1; ++i) {
+                k += ia[i] * new_strides[i];
+            }
+            for (i = 0; i < dims_axis; ++i) {
+                *(dest + j++) = *(src + k + i * strides_axis);
+            }
+        }
+    }
+}
+
+extern void zfft(complex_double * inout,
+		 int n, int direction, int howmany, int normalize);
+
+extern void zfftnd_fftpack(complex_double * inout, int rank,
+			   int *dims, int direction, int howmany,
+			   int normalize)
+{
+    int i, sz;
+    complex_double *ptr = inout;
+    int axis;
+    complex_double *tmp;
+    int *itmp;
+    int k, j;
+
+    sz = 1;
+    for (i = 0; i < rank; ++i) {
+        sz *= dims[i];
+    }
+    //zfft_fftpack(ptr, dims[rank - 1], direction, howmany * sz / dims[rank - 1],
+    //     normalize);
+    zfft(ptr, dims[rank - 1], direction, howmany * sz / dims[rank - 1],
+	 normalize);
+
+    i = get_cache_id_zfftnd_fftpack(sz, rank);
+    tmp = caches_zfftnd_fftpack[i].ptr;
+    itmp = caches_zfftnd_fftpack[i].iptr;
+
+    itmp[rank - 1] = 1;
+    for (i = 2; i <= rank; ++i) {
+        itmp[rank - i] = itmp[rank - i + 1] * dims[rank - i + 1];
+    }
+
+    for (i = 0; i < howmany; ++i, ptr += sz) {
+        for (axis = 0; axis < rank - 1; ++axis) {
+            for (k = j = 0; k < rank; ++k) {
+                if (k != axis) {
+                    *(itmp + rank + j) = itmp[k];
+                    *(itmp + 2 * rank + j++) = dims[k] - 1;
+                }
+            }
+            flatten(tmp, ptr, rank, itmp[axis], dims[axis], 0, itmp);
+            //zfft_fftpack(tmp, dims[axis], direction, sz / dims[axis], normalize);
+            zfft(tmp, dims[axis], direction, sz / dims[axis], normalize);
+            flatten(ptr, tmp, rank, itmp[axis], dims[axis], 1, itmp);
+        }
+    }
+
+}

Added: trunk/Lib/fftpack/src/zfftnd_fftw.c
===================================================================
--- trunk/Lib/fftpack/src/zfftnd_fftw.c	2007-08-07 23:55:56 UTC (rev 3222)
+++ trunk/Lib/fftpack/src/zfftnd_fftw.c	2007-08-08 06:39:47 UTC (rev 3223)
@@ -0,0 +1,53 @@
+/*
+ * fftw2 backend for multi dimensional fft
+ *
+ * Original code by Pearu Peaterson
+ *
+ * Last Change: Wed Aug 08 03:00 PM 2007 J
+ */
+
+GEN_CACHE(zfftnd_fftw, (int n, int *dims, int d, int flags)
+	  , int direction;
+	  int *dims;
+	  fftwnd_plan plan;, ((caches_zfftnd_fftw[i].n == n) &&
+			      (caches_zfftnd_fftw[i].direction == d) &&
+			      (equal_dims
+			       (n, caches_zfftnd_fftw[i].dims, dims)))
+	  , caches_zfftnd_fftw[id].direction = d;
+	  caches_zfftnd_fftw[id].n = n;
+	  caches_zfftnd_fftw[id].dims = (int *) malloc(sizeof(int) * n);
+	  memcpy(caches_zfftnd_fftw[id].dims, dims, sizeof(int) * n);
+	  caches_zfftnd_fftw[id].plan =
+	  fftwnd_create_plan(n, dims,
+			     (d > 0 ? FFTW_FORWARD : FFTW_BACKWARD),
+			     flags);,
+	  fftwnd_destroy_plan(caches_zfftnd_fftw[id].plan);
+	  free(caches_zfftnd_fftw[id].dims);, 10)
+
+
+extern void zfftnd_mkl(complex_double * inout, int rank,
+		       int *dims, int direction, int howmany,
+		       int normalize)
+{
+    int i, sz;
+    complex_double *ptr = inout;
+    fftwnd_plan plan = NULL;
+
+    sz = 1;
+    for (i = 0; i < rank; ++i) {
+	sz *= dims[i];
+    }
+    i = get_cache_id_zfftnd_fftw(rank, dims, direction,
+				 FFTW_IN_PLACE | FFTW_ESTIMATE);
+    plan = caches_zfftnd_fftw[i].plan;
+    for (i = 0; i < howmany; ++i, ptr += sz) {
+	fftwnd_one(plan, (fftw_complex *) ptr, NULL);
+    }
+    if (normalize) {
+	ptr = inout;
+	for (i = sz * howmany - 1; i >= 0; --i) {
+	    *((double *) (ptr)) /= sz;
+	    *((double *) (ptr++) + 1) /= sz;
+	}
+    }
+}

Added: trunk/Lib/fftpack/src/zfftnd_fftw3.c
===================================================================
--- trunk/Lib/fftpack/src/zfftnd_fftw3.c	2007-08-07 23:55:56 UTC (rev 3222)
+++ trunk/Lib/fftpack/src/zfftnd_fftw3.c	2007-08-08 06:39:47 UTC (rev 3223)
@@ -0,0 +1,40 @@
+/*
+ * fftw3 backend for multi dimensional fft
+ *
+ * Original code by Pearu Peaterson
+ *
+ * Last Change: Wed Aug 08 02:00 PM 2007 J
+ */
+
+/* stub because fftw3 has no cache mechanism (yet) */
+static void destroy_zfftnd_fftw3_caches(void) {}
+
+extern void zfftnd_fftw3(complex_double * inout, int rank,
+			   int *dims, int direction, int howmany,
+			   int normalize)
+{
+    int i, sz;
+    complex_double *ptr = inout;
+
+    fftw_plan plan = NULL;
+    sz = 1;
+    for (i = 0; i < rank; ++i) {
+        sz *= dims[i];
+    }
+    plan = fftw_plan_many_dft(rank, dims, howmany,
+			      (fftw_complex *) ptr, NULL, 1, sz,
+			      (fftw_complex *) ptr, NULL, 1, sz,
+			      (direction >
+			       0 ? FFTW_FORWARD : FFTW_BACKWARD),
+			      FFTW_ESTIMATE);
+    fftw_execute(plan);
+    fftw_destroy_plan(plan);
+
+    if (normalize) {
+        ptr = inout;
+        for (i = sz * howmany - 1; i >= 0; --i) {
+            *((double *) (ptr)) /= sz;
+            *((double *) (ptr++) + 1) /= sz;
+        }
+    }
+}

Added: trunk/Lib/fftpack/src/zfftnd_mkl.c
===================================================================
--- trunk/Lib/fftpack/src/zfftnd_mkl.c	2007-08-07 23:55:56 UTC (rev 3222)
+++ trunk/Lib/fftpack/src/zfftnd_mkl.c	2007-08-08 06:39:47 UTC (rev 3223)
@@ -0,0 +1,66 @@
+/*
+ * MKL backend for multi dimensional fft
+ *
+ * Original code by David M. Cooke
+ *
+ * Last Change: Wed Aug 08 03:00 PM 2007 J
+ */
+
+GEN_CACHE(zfftnd_mkl, (int n, int *dims)
+	  , DFTI_DESCRIPTOR_HANDLE desc_handle;
+	  int *dims;
+	  long *ndims;, ((caches_zfftnd_mkl[i].n == n) &&
+			 (equal_dims(n, caches_zfftnd_mkl[i].dims, dims)))
+	  , caches_zfftnd_mkl[id].ndims = convert_dims(n, dims);
+	  caches_zfftnd_mkl[id].n = n;
+	  caches_zfftnd_mkl[id].dims = (int *) malloc(sizeof(int) * n);
+	  memcpy(caches_zfftnd_mkl[id].dims, dims, sizeof(int) * n);
+	  DftiCreateDescriptor(&caches_zfftnd_mkl[id].desc_handle,
+			       DFTI_DOUBLE, DFTI_COMPLEX, (long) n,
+			       caches_zfftnd_mkl[id].ndims);
+	  DftiCommitDescriptor(caches_zfftnd_mkl[id].desc_handle);,
+	  DftiFreeDescriptor(&caches_zfftnd_mkl[id].desc_handle);
+	  free(caches_zfftnd_mkl[id].dims);
+	  free(caches_zfftnd_mkl[id].ndims);, 10)
+
+static 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;
+}
+
+extern void zfftnd_mkl(complex_double * inout, int rank,
+		       int *dims, int direction, int howmany,
+		       int normalize)
+{
+    int i, sz;
+    complex_double *ptr = inout;
+
+    DFTI_DESCRIPTOR_HANDLE desc_handle;
+    sz = 1;
+    for (i = 0; i < rank; ++i) {
+        sz *= dims[i];
+    }
+
+    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;
+        }
+    }
+}



More information about the Scipy-svn mailing list