[Scipy-svn] r4266 - branches/refactor_fft/scipy/fftpack/src/djbfft
scipy-svn at scipy.org
scipy-svn at scipy.org
Sun May 11 07:35:17 EDT 2008
Author: cdavid
Date: 2008-05-11 06:35:13 -0500 (Sun, 11 May 2008)
New Revision: 4266
Modified:
branches/refactor_fft/scipy/fftpack/src/djbfft/common.h
branches/refactor_fft/scipy/fftpack/src/djbfft/drfft.cxx
Log:
djbfft backend for drfft now uses c++ cycliccache.
Modified: branches/refactor_fft/scipy/fftpack/src/djbfft/common.h
===================================================================
--- branches/refactor_fft/scipy/fftpack/src/djbfft/common.h 2008-05-11 11:04:37 UTC (rev 4265)
+++ branches/refactor_fft/scipy/fftpack/src/djbfft/common.h 2008-05-11 11:35:13 UTC (rev 4266)
@@ -1,9 +1,15 @@
#ifndef _SCIPY_DJBFFT_COMMON_H
#define _SCIPY_DJBFFT_COMMON_H
+#include <cycliccache.h>
+
+namespace fft {
+
class DJBFFTCacheId : public CacheId {
public:
DJBFFTCacheId(int n) : CacheId(n) {};
};
+};
+
#endif
Modified: branches/refactor_fft/scipy/fftpack/src/djbfft/drfft.cxx
===================================================================
--- branches/refactor_fft/scipy/fftpack/src/djbfft/drfft.cxx 2008-05-11 11:04:37 UTC (rev 4265)
+++ branches/refactor_fft/scipy/fftpack/src/djbfft/drfft.cxx 2008-05-11 11:35:13 UTC (rev 4266)
@@ -1,15 +1,19 @@
/*
- * Last Change: Wed Aug 01 08:00 PM 2007 J
+ * Last Change: Sun May 11 08:00 PM 2008 J
*
* Original code by Pearu Peterson.
*/
-
/*
- * DJBFFT only implements size 2^N !
+ * RDJBFFT only implements size 2^N !
*
* drfft_def and drfft_def_destroy_cache are the functions used for size different
* than 2^N
*/
+#include <new>
+#include <cassert>
+
+#include "common.h"
+
#ifdef WITH_FFTW3
#define drfft_def drfft_fftw3
#define drfft_def_destroy_cache destroy_drfftw3_caches
@@ -21,6 +25,7 @@
#define drfft_def_destroy_cache destroy_drfftpack_caches
#endif
+#if 0
GEN_CACHE(drdjbfft, (int n)
, unsigned int *f;
double *ptr;,
@@ -31,101 +36,169 @@
free(caches_drdjbfft[id].f);
free(caches_drdjbfft[id].ptr);,
10)
+#endif
-/**************** ZFFT function **********************/
-static void drfft_djbfft(double * inout,
- int n, int direction, int howmany, int normalize)
+using namespace fft;
+
+class RDJBFFTCache: public Cache<DJBFFTCacheId> {
+ public:
+ RDJBFFTCache(const DJBFFTCacheId& id);
+ virtual ~RDJBFFTCache();
+
+ int compute_forward(double * inout) const;
+ int compute_backward(double * inout, int normalize) const;
+
+ protected:
+ unsigned int* m_f;
+ double* m_ptr;
+};
+
+RDJBFFTCache::RDJBFFTCache(const DJBFFTCacheId& id)
+: Cache<DJBFFTCacheId>(id)
{
- int i;
- double *ptr = inout;
- double *ptrc = NULL;
- unsigned int *f = NULL;
+ int n = id.m_n;
- switch (n) {
- case 2:;
- case 4:;
- case 8:;
- case 16:;
- case 32:;
- case 64:;
- case 128:;
- case 256:;
- case 512:;
- case 1024:;
- case 2048:;
- case 4096:;
- case 8192:
- i = get_cache_id_drdjbfft(n);
- f = caches_drdjbfft[i].f;
- ptrc = caches_drdjbfft[i].ptr;
- }
- if (f == NULL) {
- drfft_def(inout, n, direction, howmany, normalize);
- }
+ m_f = (unsigned int*)malloc(sizeof(*m_f) * n);
+ if (m_f == NULL) {
+ goto fail_f;
+ }
- switch (direction) {
- case 1:
- for (i = 0; i < howmany; ++i, ptr += n) {
- if (f != NULL) {
- COPYSTD2DJB(ptr, ptrc, n);
- switch (n) {
-#define TMPCASE(N) case N: fftr8_##N(ptrc); break
- TMPCASE(2);
- TMPCASE(4);
- TMPCASE(8);
- TMPCASE(16);
- TMPCASE(32);
- TMPCASE(64);
- TMPCASE(128);
- TMPCASE(256);
- TMPCASE(512);
- TMPCASE(1024);
- TMPCASE(2048);
- TMPCASE(4096);
- TMPCASE(8192);
+ m_ptr = (double *)malloc(sizeof(*m_ptr) * 2 * n);
+ if (m_ptr == NULL) {
+ goto clean_f;
+ }
+
+ fftfreq_rtable(m_f, id.m_n);
+ return;
+
+clean_f:
+ free(m_f);
+fail_f:
+ throw std::bad_alloc();
+}
+
+RDJBFFTCache::~RDJBFFTCache()
+{
+ free(m_ptr);
+ free(m_f);
+}
+
+int RDJBFFTCache::compute_forward(double *inout) const
+{
+ const int n = m_id.m_n;
+
+ COPYSTD2DJB(inout, m_ptr, n);
+ switch (n) {
+#define TMPCASE(N) case N: fftr8_##N(m_ptr); break
+ TMPCASE(2);
+ TMPCASE(4);
+ TMPCASE(8);
+ TMPCASE(16);
+ TMPCASE(32);
+ TMPCASE(64);
+ TMPCASE(128);
+ TMPCASE(256);
+ TMPCASE(512);
+ TMPCASE(1024);
+ TMPCASE(2048);
+ TMPCASE(4096);
+ TMPCASE(8192);
#undef TMPCASE
- }
- COPYDJB2STD(ptrc, ptr, f, n);
- }
}
- break;
+ COPYDJB2STD(m_ptr, inout, m_f, n);
- case -1:
- for (i = 0; i < howmany; ++i, ptr += n) {
- if (f != NULL) {
- COPYINVSTD2DJB(ptr, ptrc, normalize, f, n);
- switch (n) {
+ return 0;
+}
-#define TMPCASE(N)case N:if(normalize)fftr8_scale##N(ptrc);fftr8_un##N(ptrc);break
- TMPCASE(2);
- TMPCASE(4);
- TMPCASE(8);
- TMPCASE(16);
- TMPCASE(32);
- TMPCASE(64);
- TMPCASE(128);
- TMPCASE(256);
- TMPCASE(512);
- TMPCASE(1024);
- TMPCASE(2048);
- TMPCASE(4096);
- TMPCASE(8192);
+int RDJBFFTCache::compute_backward(double *inout, int normalize) const
+{
+ const int n = m_id.m_n;
+
+ COPYINVSTD2DJB(inout, m_ptr, normalize, m_f, n);
+ switch (n) {
+#define TMPCASE(N)case N:if(normalize)fftr8_scale##N(m_ptr);fftr8_un##N(m_ptr);break
+ TMPCASE(2);
+ TMPCASE(4);
+ TMPCASE(8);
+ TMPCASE(16);
+ TMPCASE(32);
+ TMPCASE(64);
+ TMPCASE(128);
+ TMPCASE(256);
+ TMPCASE(512);
+ TMPCASE(1024);
+ TMPCASE(2048);
+ TMPCASE(4096);
+ TMPCASE(8192);
#undef TMPCASE
- }
- COPYINVDJB2STD(ptrc, ptr, n);
- }
}
- break;
+ COPYINVDJB2STD(m_ptr, inout, n);
- default:
- fprintf(stderr, "drfft: invalid direction=%d\n", direction);
- }
+ return 0;
+}
- if (normalize && f != NULL && direction == 1) {
- double d = 1.0 / n;
- ptr = inout;
- for (i = n * howmany - 1; i >= 0; --i) {
- (*(ptr++)) *= d;
+static CacheManager<DJBFFTCacheId, RDJBFFTCache> rdjbfft_cmgr(10);
+
+/* Stub to make GEN_PUBLIC_API happy */
+static void destroy_drdjbfft_caches()
+{
+}
+
+/**************** ZFFT function **********************/
+static void drfft_djbfft(double * inout,
+ int n, int direction, int howmany, int normalize)
+{
+ int i;
+ double *ptr = inout;
+ RDJBFFTCache *cache;
+ unsigned int *f = NULL;
+
+ switch (n) {
+ case 2:;
+ case 4:;
+ case 8:;
+ case 16:;
+ case 32:;
+ case 64:;
+ case 128:;
+ case 256:;
+ case 512:;
+ case 1024:;
+ case 2048:;
+ case 4096:;
+ case 8192:
+ cache = rdjbfft_cmgr.get_cache(DJBFFTCacheId(n));
+ break;
+ default:
+ /* For sizes not handled by djbfft, use default
+ * implementation and returns */
+ drfft_def(inout, n, direction, howmany, normalize);
+ return;
}
- }
+
+ switch (direction) {
+ case 1:
+ for (i = 0; i < howmany; ++i, ptr += n) {
+ cache->compute_forward(ptr);
+ }
+ break;
+
+ case -1:
+ for (i = 0; i < howmany; ++i, ptr += n) {
+ cache->compute_backward(ptr, normalize);
+ }
+ break;
+
+ default:
+ fprintf(stderr, "drfft: invalid direction=%d\n",
+ direction);
+ }
+
+ if (normalize && f != NULL && direction == 1) {
+ double d = 1.0 / n;
+ ptr = inout;
+ for (i = n * howmany - 1; i >= 0; --i) {
+ (*(ptr++)) *= d;
+ }
+ }
}
More information about the Scipy-svn
mailing list