[Scipy-svn] r5370 - trunk/scipy/fftpack/src

scipy-svn@scip... scipy-svn@scip...
Wed Jan 7 12:26:40 CST 2009


Author: cdavid
Date: 2009-01-07 12:26:15 -0600 (Wed, 07 Jan 2009)
New Revision: 5370

Modified:
   trunk/scipy/fftpack/src/zfftnd.c
Log:
Put all nd code in zfftnd.c

Modified: trunk/scipy/fftpack/src/zfftnd.c
===================================================================
--- trunk/scipy/fftpack/src/zfftnd.c	2009-01-07 18:25:47 UTC (rev 5369)
+++ trunk/scipy/fftpack/src/zfftnd.c	2009-01-07 18:26:15 UTC (rev 5370)
@@ -5,31 +5,210 @@
  */
 #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);\
-}\
-void destroy_cfftnd_cache(void)\
-{\
-        destroy_cfftnd_##name##_caches();\
-}\
-\
-void cfftnd(complex_float * inout, int rank,\
-		           int *dims, int direction, int howmany, int normalize)\
-{\
-        cfftnd_##name(inout, rank, dims, direction, howmany, normalize);\
+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)
+
+GEN_CACHE(cfftnd_fftpack, (int n, int rank)
+	  , complex_float * ptr; int *iptr; int rank;
+	  , ((caches_cfftnd_fftpack[i].n == n)
+	     && (caches_cfftnd_fftpack[i].rank == rank))
+	  , caches_cfftnd_fftpack[id].n = n;
+	  caches_cfftnd_fftpack[id].ptr =
+	  (complex_float *) malloc(2 * sizeof(float) * n);
+	  caches_cfftnd_fftpack[id].iptr =
+	  (int *) malloc(4 * rank * sizeof(int));
+	  ,
+	  free(caches_cfftnd_fftpack[id].ptr);
+	  free(caches_cfftnd_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);
+            }
+        }
+    }
+}
 
-#include "zfftnd_fftpack.c"
-#include "cfftnd_fftpack.c"
-GEN_PUBLIC_API(fftpack)
+static
+void sflatten(complex_float * dest, complex_float * 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 cfft(complex_float * inout,
+		 int n, int direction, int howmany, int normalize);
+
+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(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(tmp, dims[axis], direction, sz / dims[axis], normalize);
+            flatten(ptr, tmp, rank, itmp[axis], dims[axis], 1, itmp);
+        }
+    }
+
+}
+
+extern void cfftnd_fftpack(complex_float * inout, int rank,
+			   int *dims, int direction, int howmany,
+			   int normalize)
+{
+    int i, sz;
+    complex_float *ptr = inout;
+    int axis;
+    complex_float *tmp;
+    int *itmp;
+    int k, j;
+
+    sz = 1;
+    for (i = 0; i < rank; ++i) {
+        sz *= dims[i];
+    }
+    cfft(ptr, dims[rank - 1], direction, howmany * sz / dims[rank - 1],
+	 normalize);
+
+    i = get_cache_id_cfftnd_fftpack(sz, rank);
+    tmp = caches_cfftnd_fftpack[i].ptr;
+    itmp = caches_cfftnd_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;
+                }
+            }
+            sflatten(tmp, ptr, rank, itmp[axis], dims[axis], 0, itmp);
+            cfft(tmp, dims[axis], direction, sz / dims[axis], normalize);
+            sflatten(ptr, tmp, rank, itmp[axis], dims[axis], 1, itmp);
+        }
+    }
+
+}



More information about the Scipy-svn mailing list