[Scipy-svn] r5370 - trunk/scipy/fftpack/src
scipy-svn at scipy.org
scipy-svn at scipy.org
Wed Jan 7 13:26:40 EST 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