[Scipy-svn] r4360 - branches/refactor_fft/scipy/fftpack/backends/fftw
scipy-svn at scipy.org
scipy-svn at scipy.org
Fri May 16 02:03:51 EDT 2008
Author: cdavid
Date: 2008-05-16 01:03:43 -0500 (Fri, 16 May 2008)
New Revision: 4360
Added:
branches/refactor_fft/scipy/fftpack/backends/fftw/__init__.py
branches/refactor_fft/scipy/fftpack/backends/fftw/fftw.pyf
branches/refactor_fft/scipy/fftpack/backends/fftw/setup.py
Removed:
branches/refactor_fft/scipy/fftpack/backends/fftw/api.h
branches/refactor_fft/scipy/fftpack/backends/fftw/common.h
branches/refactor_fft/scipy/fftpack/backends/fftw/convolve.cxx
branches/refactor_fft/scipy/fftpack/backends/fftw/drfft.cxx
branches/refactor_fft/scipy/fftpack/backends/fftw/zfft.cxx
branches/refactor_fft/scipy/fftpack/backends/fftw/zfftnd.cxx
Log:
fftw backend now builds
Added: branches/refactor_fft/scipy/fftpack/backends/fftw/__init__.py
===================================================================
--- branches/refactor_fft/scipy/fftpack/backends/fftw/__init__.py 2008-05-16 06:02:54 UTC (rev 4359)
+++ branches/refactor_fft/scipy/fftpack/backends/fftw/__init__.py 2008-05-16 06:03:43 UTC (rev 4360)
@@ -0,0 +1,3 @@
+from _fftw import zfft_fftw as zfft, \
+ drfft_fftw as drfft, \
+ zfftnd_fftw as zfftnd
Deleted: branches/refactor_fft/scipy/fftpack/backends/fftw/api.h
===================================================================
--- branches/refactor_fft/scipy/fftpack/backends/fftw/api.h 2008-05-16 06:02:54 UTC (rev 4359)
+++ branches/refactor_fft/scipy/fftpack/backends/fftw/api.h 2008-05-16 06:03:43 UTC (rev 4360)
@@ -1,30 +0,0 @@
-#ifndef _SCIPY_FFTPACK_FFTW_API_H_
-#define _SCIPY_FFTPACK_FFTW_API_H_
-
-#include "misc.h"
-
-/*
- * straight FFT api
- */
-void drfft_fftw(double * inout, int n, int direction, int howmany,
- int normalize);
-
-void zfft_fftw(complex_double * inout,
- int n, int direction, int howmany, int normalize);
-
-void zfftnd_fftw(complex_double * inout, int rank,
- int *dims, int direction, int howmany,
- int normalize);
-
-/*
- * Convolution api
- */
-void convolve_fftw(int n, double *inout, double *omega, int swap_real_imag);
-void convolve_z_fftw(int n, double *inout, double *omega_real,
- double* omega_imag);
-
-void init_convolution_kernel_fftw(int n, double *omega, int d,
- double (*kernel_func) (int),
- int zero_nyquist);
-
-#endif
Deleted: branches/refactor_fft/scipy/fftpack/backends/fftw/common.h
===================================================================
--- branches/refactor_fft/scipy/fftpack/backends/fftw/common.h 2008-05-16 06:02:54 UTC (rev 4359)
+++ branches/refactor_fft/scipy/fftpack/backends/fftw/common.h 2008-05-16 06:03:43 UTC (rev 4360)
@@ -1,31 +0,0 @@
-#ifndef _SCIPY_FFTW_COMMON_H
-#define _SCIPY_FFTW_COMMON_H
-
-#define COPYRFFTW2STD(SRC,DEST,N) { \
- int j,n2=(N)/2; \
- *(DEST) = *(SRC); \
- for (j=1;j<n2;++j) { \
- *(DEST+2*j-1) = *(SRC+j); \
- *(DEST+2*j) = *(SRC+(N)-j); \
- } \
- if (N>1) { \
- *(DEST+2*n2-1) = *(SRC+n2); \
- if ((N)%2) \
- *(DEST+2*n2) = *(SRC+(N)-n2); \
- } \
-}
-#define COPYINVRFFTW2STD(SRC,DEST,N) { \
- int j,n2=(N)/2; \
- *(DEST) = *(SRC); \
- for (j=1;j<n2;++j) { \
- *(DEST+j) = *(SRC+2*j-1); \
- *(DEST+(N)-j) = *(SRC+2*j); \
- } \
- if (N>1) {\
- *(DEST+n2) = *(SRC+2*n2-1); \
- if ((N)%2) \
- *(DEST+(N)-n2) = *(SRC+2*n2); \
- } \
-}
-
-#endif
Deleted: branches/refactor_fft/scipy/fftpack/backends/fftw/convolve.cxx
===================================================================
--- branches/refactor_fft/scipy/fftpack/backends/fftw/convolve.cxx 2008-05-16 06:02:54 UTC (rev 4359)
+++ branches/refactor_fft/scipy/fftpack/backends/fftw/convolve.cxx 2008-05-16 06:03:43 UTC (rev 4360)
@@ -1,183 +0,0 @@
-#include <new>
-
-#include <rfftw.h>
-#include <fftw.h>
-
-#include "api.h"
-
-#include "cycliccache.h"
-
-using namespace fft;
-
-class DRFFTWCacheId : public CacheId {
- public:
- DRFFTWCacheId(int n);
-
-};
-
-DRFFTWCacheId::DRFFTWCacheId(int n):
- CacheId(n)
-{
-}
-
-class DRFFTWCache : public Cache<DRFFTWCacheId> {
- public:
- DRFFTWCache(const DRFFTWCacheId& id);
- virtual ~DRFFTWCache();
-
- int convolve(double* inout, double* omega,
- int swap_real_imag) const;
- int convolve_z(double* inout, double* omega_real,
- double* omega_imag) const;
-
- protected:
- rfftw_plan m_plan1;
- rfftw_plan m_plan2;
-};
-
-DRFFTWCache::DRFFTWCache(const DRFFTWCacheId& id)
-: Cache<DRFFTWCacheId>(id)
-{
- int flags = FFTW_ESTIMATE | FFTW_IN_PLACE;
-
- m_plan1 = rfftw_create_plan(id.m_n, FFTW_REAL_TO_COMPLEX, flags);
- if (m_plan1 == NULL) {
- goto fail;
- }
-
- m_plan2 = rfftw_create_plan(id.m_n, FFTW_COMPLEX_TO_REAL, flags);
- if (m_plan2 == NULL) {
- goto clean_plan1;
- }
-
- return;
-
-clean_plan1:
- rfftw_destroy_plan(m_plan1);
-fail:
- throw std::bad_alloc();
-}
-
-DRFFTWCache::~DRFFTWCache()
-{
- rfftw_destroy_plan(m_plan2);
- rfftw_destroy_plan(m_plan1);
-}
-
-int DRFFTWCache::convolve(double* inout, double* omega, int swap_real_imag)
- const
-{
- int n = m_id.m_n;
- int l = (n-1)/2+1;
- int i;
- double c;
-
- rfftw_one(m_plan1, (fftw_real *)inout, NULL);
- if (swap_real_imag) {
- inout[0] *= omega[0];
- if (!(n%2)) {
- inout[n/2] *= omega[n/2];
- }
- for(i=1;i<l;++i) {
- c = inout[i] * omega[i];
- inout[i] = omega[n-i] * inout[n-i];
- inout[n-i] = c;
- }
- } else {
- for(i=0;i<n;++i) {
- inout[i] *= omega[i];
- }
- }
- rfftw_one(m_plan2, (fftw_real *)inout, NULL);
-
- return 0;
-}
-
-int DRFFTWCache::convolve_z(double* inout, double* omega_real,
- double* omega_imag) const
-{
- int n = m_id.m_n;
- int i;
- int l = (n-1)/2+1;
- double c;
-
- rfftw_one(m_plan1, (fftw_real *)inout, NULL);
- inout[0] *= (omega_real[0]+omega_imag[0]);
-
- if (!(n%2)) {
- inout[n/2] *= (omega_real[n/2]+omega_imag[n/2]);
- }
- for(i=1;i<l;++i) {
- c = inout[i] * omega_imag[i];
- inout[i] *= omega_real[i];
- inout[i] += omega_imag[n-i] * inout[n-i];
- inout[n-i] *= omega_real[n-i];
- inout[n-i] += c;
- }
- rfftw_one(m_plan2, (fftw_real *)inout, NULL);
-
- return 0;
-}
-
-CacheManager<DRFFTWCacheId, DRFFTWCache> drfftw_cmgr(20);
-
-/**************** convolve **********************/
-void convolve_fftw(int n,double* inout,double* omega,int swap_real_imag)
-{
- DRFFTWCache *cache;
-
- cache = drfftw_cmgr.get_cache(DRFFTWCacheId(n));
- cache->convolve(inout, omega, swap_real_imag);
-}
-
-/**************** convolve **********************/
-void convolve_z_fftw(int n,double* inout,double* omega_real,double* omega_imag)
-{
- DRFFTWCache *cache;
-
- cache = drfftw_cmgr.get_cache(DRFFTWCacheId(n));
- cache->convolve_z(inout, omega_real, omega_imag);
-}
-
-void init_convolution_kernel_fftw(int n,double* omega, int d,
- double (*kernel_func)(int),
- int zero_nyquist)
-{
- /*
- * omega[k] = pow(sqrt(-1),d) * kernel_func(k)
- * omega[0] = kernel_func(0)
- * conjugate(omega[-k]) == omega[k]
- */
- int k,l=(n-1)/2+1;
- omega[0] = (*kernel_func)(0)/n;;
- switch (d%4) {
- case 0:
- for (k=1;k<l;++k)
- omega[k] = omega[n-k] = (*kernel_func)(k)/n;
- if (!(n%2))
- omega[n/2] = (zero_nyquist?0.0:(*kernel_func)(n/2)/n);
- break;
- case 1:;case -3:
- for (k=1;k<l;++k) {
- omega[k] = (*kernel_func)(k)/n;
- omega[n-k] = -omega[k];
- }
- if (!(n%2))
- omega[n/2] = (zero_nyquist?0.0:(*kernel_func)(n/2)/n);
- break;
- case 2:;case -2:
- for (k=1;k<l;++k)
- omega[k] = omega[n-k] = -(*kernel_func)(k)/n;
- if (!(n%2))
- omega[n/2] = (zero_nyquist?0.0:-(*kernel_func)(n/2)/n);
- break;
- case 3:;case -1:
- for (k=1;k<l;++k) {
- omega[k] = -(*kernel_func)(k)/n;
- omega[n-k] = -omega[k];
- }
- if (!(n%2))
- omega[n/2] = (zero_nyquist?0.0:-(*kernel_func)(n/2)/n);
- break;
- }
-}
Deleted: branches/refactor_fft/scipy/fftpack/backends/fftw/drfft.cxx
===================================================================
--- branches/refactor_fft/scipy/fftpack/backends/fftw/drfft.cxx 2008-05-16 06:02:54 UTC (rev 4359)
+++ branches/refactor_fft/scipy/fftpack/backends/fftw/drfft.cxx 2008-05-16 06:03:43 UTC (rev 4360)
@@ -1,135 +0,0 @@
-/*
- * Last Change: Tue May 13 06:00 PM 2008 J
- *
- * FFTW2 implementation
- *
- * Original code by Pearu Peterson.
- */
-#include <new>
-
-#include <rfftw.h>
-#include <fftw.h>
-
-#include "cycliccache.h"
-#include "api.h"
-#include "common.h"
-
-using namespace fft;
-
-class RFFTWCacheId : public CacheId {
- public:
- RFFTWCacheId(int n, int dir, int flags);
-
- virtual bool operator==(const RFFTWCacheId& other) const
- {
- return is_equal(other);
- }
-
- virtual bool is_equal(const RFFTWCacheId& other) const
- {
- const CacheId *ot = &other;
- const CacheId *th = this;
-
- return m_dir == other.m_dir &&
- m_flags == other.m_flags &&
- th->is_equal(*ot);
- }
-
- public:
- int m_dir;
- int m_flags;
-};
-
-RFFTWCacheId::RFFTWCacheId(int n, int dir, int flags):
- CacheId(n),
- m_dir(dir),
- m_flags(flags)
-{
-}
-
-class RFFTWCache : public Cache<RFFTWCacheId> {
- public:
- RFFTWCache(const RFFTWCacheId& id);
- virtual ~RFFTWCache();
-
- int compute(double* inout) const;
-
- protected:
- rfftw_plan m_plan;
- double *m_wrk;
-};
-
-RFFTWCache::RFFTWCache(const RFFTWCacheId& id)
-: Cache<RFFTWCacheId>(id)
-{
- m_wrk = (double *) malloc(sizeof(double) * id.m_n);
- if(m_wrk == NULL) {
- goto fail;
- }
-
- m_plan = rfftw_create_plan(id.m_n,
- (id.m_dir > 0 ? FFTW_REAL_TO_COMPLEX :
- FFTW_COMPLEX_TO_REAL),
- id.m_flags);
-
- if (m_plan == NULL) {
- goto clean_wrk;
- }
-
- return;
-
-clean_wrk:
- free(m_wrk);
-fail:
- throw std::bad_alloc();
-}
-
-RFFTWCache::~RFFTWCache()
-{
- free(m_wrk);
- rfftw_destroy_plan(m_plan);
-}
-
-int RFFTWCache::compute(double* inout) const
-{
- if(m_id.m_dir == 1) {
- memcpy(m_wrk, inout, sizeof(double) * m_id.m_n);
- rfftw(m_plan, 1, (fftw_real *) m_wrk, 1, 1, NULL, 1, 1);
- COPYRFFTW2STD(m_wrk, inout, m_id.m_n);
- } else {
- COPYINVRFFTW2STD(inout, m_wrk, m_id.m_n);
- rfftw(m_plan, 1, (fftw_real *) m_wrk, 1, 1, NULL, 1, 1);
- memcpy(inout, m_wrk, sizeof(double) * m_id.m_n);
- }
- return 0;
-};
-
-CacheManager<RFFTWCacheId, RFFTWCache> rfftw_cmgr(10);
-
-void drfft_fftw(double *inout, int n, int dir,
- int howmany, int normalize)
-{
- int i;
- double *ptr = inout;
- RFFTWCache *cache;
-
- cache = rfftw_cmgr.get_cache(
- RFFTWCacheId(n, dir, FFTW_IN_PLACE | FFTW_ESTIMATE));
-
- if (dir != -1 && dir != 1) {
- fprintf(stderr, "drfft: invalid dir=%d\n", dir);
- return;
- } else {
- for (i = 0; i < howmany; ++i, ptr += n) {
- cache->compute(ptr);
- }
- }
-
- if (normalize) {
- double d = 1.0 / n;
- ptr = inout;
- for (i = n * howmany - 1; i >= 0; --i) {
- (*(ptr++)) *= d;
- }
- }
-}
Added: branches/refactor_fft/scipy/fftpack/backends/fftw/fftw.pyf
===================================================================
--- branches/refactor_fft/scipy/fftpack/backends/fftw/fftw.pyf 2008-05-16 06:02:54 UTC (rev 4359)
+++ branches/refactor_fft/scipy/fftpack/backends/fftw/fftw.pyf 2008-05-16 06:03:43 UTC (rev 4360)
@@ -0,0 +1,77 @@
+!%f90 -*- f90 -*-
+! Author: Pearu Peterson, August 2002
+
+python module _fftw
+ interface
+
+ subroutine zfft_fftw(x,n,direction,howmany,normalize)
+ ! y = fft(x[,n,direction,normalize,overwrite_x])
+ intent(c) zfft_fftw
+ complex*16 intent(c,in,out,copy,out=y) :: x(*)
+ integer optional,depend(x),intent(c,in) :: n=size(x)
+ check(n>0) n
+ integer depend(x,n),intent(c,hide) :: howmany = size(x)/n
+ check(n*howmany==size(x)) howmany
+ integer optional,intent(c,in) :: direction = 1
+ integer optional,intent(c,in),depend(direction) &
+ :: normalize = (direction<0)
+ end subroutine zfft_fftw
+
+ subroutine drfft_fftw(x,n,direction,howmany,normalize)
+ ! y = drfft(x[,n,direction,normalize,overwrite_x])
+ intent(c) drfft_fftw
+ real*8 intent(c,in,out,copy,out=y) :: x(*)
+ integer optional,depend(x),intent(c,in) :: n=size(x)
+ check(n>0&&n<=size(x)) n
+ integer depend(x,n),intent(c,hide) :: howmany = size(x)/n
+ check(n*howmany==size(x)) howmany
+ integer optional,intent(c,in) :: direction = 1
+ integer optional,intent(c,in),depend(direction) &
+ :: normalize = (direction<0)
+ end subroutine drfft_fftw
+
+! subroutine zrfft_fftw(x,n,direction,howmany,normalize)
+! ! y = zrfft(x[,n,direction,normalize,overwrite_x])
+! intent(c) zrfft_fftw
+! complex*16 intent(c,in,out,overwrite,out=y) :: x(*)
+! integer optional,depend(x),intent(c,in) :: n=size(x)
+! check(n>0&&n<=size(x)) n
+! integer depend(x,n),intent(c,hide) :: howmany = size(x)/n
+! check(n*howmany==size(x)) howmany
+! integer optional,intent(c,in) :: direction = 1
+! integer optional,intent(c,in),depend(direction) &
+! :: normalize = (direction<0)
+! end subroutine zrfft_fftw
+
+ subroutine zfftnd_fftw(x,r,s,direction,howmany,normalize,j)
+ ! y = zfftnd(x[,s,direction,normalize,overwrite_x])
+ intent(c) zfftnd_fftw
+ complex*16 intent(c,in,out,copy,out=y) :: x(*)
+ integer intent(c,hide),depend(x) :: r=old_rank(x)
+ integer intent(c,hide) :: j=0
+ integer optional,depend(r),dimension(r),intent(c,in) &
+ :: s=old_shape(x,j++)
+ check(r>=len(s)) s
+ integer intent(c,hide) :: howmany = 1
+ integer optional,intent(c,in) :: direction = 1
+ integer optional,intent(c,in),depend(direction) :: &
+ normalize = (direction<0)
+ callprotoargument complex_double*,int,int*,int,int,int
+ callstatement {&
+ int i,sz=1,xsz=size(x); &
+ for (i=0;i<r;++i) sz *= s[i]; &
+ howmany = xsz/sz; &
+ if (sz*howmany==xsz) &
+ (*f2py_func)(x,r,s,direction,howmany,normalize); &
+ else {&
+ f2py_success = 0; &
+ PyErr_SetString(_fftw_error, &
+ "inconsistency in x.shape and s argument"); &
+ } &
+ }
+ end subroutine zfftnd_fftw
+
+ end interface
+end python module _fftw
+
+! See http://cens.ioc.ee/projects/f2py2e/
Added: branches/refactor_fft/scipy/fftpack/backends/fftw/setup.py
===================================================================
--- branches/refactor_fft/scipy/fftpack/backends/fftw/setup.py 2008-05-16 06:02:54 UTC (rev 4359)
+++ branches/refactor_fft/scipy/fftpack/backends/fftw/setup.py 2008-05-16 06:03:43 UTC (rev 4360)
@@ -0,0 +1,28 @@
+#!/usr/bin/env python
+# Created by Pearu Peterson, August 2002
+
+from os.path import join
+
+def build_backends(config):
+ from numpy.distutils.system_info import get_info
+ src = [join('src', i) for i in
+ ['zfft.cxx', 'drfft.cxx', 'zfftnd.cxx']]
+ info = get_info("fftw2")
+ if info:
+ config.add_library("fftw_backend",
+ sources = src,
+ include_dirs = ["../common",
+ info['include_dirs']])
+ config.add_extension("_fftw", sources = ["fftw.pyf"], extra_info = info, libraries = ["fftw_backend"])
+
+def configuration(parent_package='',top_path=None):
+ from numpy.distutils.misc_util import Configuration
+
+ config = Configuration('fftw', parent_package, top_path)
+ build_backends(config)
+
+ return config
+
+if __name__ == '__main__':
+ from numpy.distutils.core import setup
+ setup(**configuration(top_path='').todict())
Property changes on: branches/refactor_fft/scipy/fftpack/backends/fftw/setup.py
___________________________________________________________________
Name: svn:executable
+ *
Deleted: branches/refactor_fft/scipy/fftpack/backends/fftw/zfft.cxx
===================================================================
--- branches/refactor_fft/scipy/fftpack/backends/fftw/zfft.cxx 2008-05-16 06:02:54 UTC (rev 4359)
+++ branches/refactor_fft/scipy/fftpack/backends/fftw/zfft.cxx 2008-05-16 06:03:43 UTC (rev 4360)
@@ -1,95 +0,0 @@
-#include <new>
-
-#include <rfftw.h>
-#include <fftw.h>
-
-#include "cycliccache.h"
-#include "api.h"
-
-using namespace fft;
-
-class FFTWCacheId : public CacheId {
- public:
- FFTWCacheId(int n, int dir) : CacheId(n), m_dir(dir) {};
-
- virtual bool operator==(const FFTWCacheId& other) const
- {
- return is_equal(other);
- }
-
- virtual bool is_equal(const FFTWCacheId& other) const
- {
- const CacheId *ot = &other;
- const CacheId *th = this;
-
- return m_dir == other.m_dir && th->is_equal(*ot);
- }
-
- public:
- int m_dir;
-};
-
-class FFTWCache : public Cache<FFTWCacheId> {
- public:
- FFTWCache(const FFTWCacheId& id);
- virtual ~FFTWCache();
-
- int compute(fftw_complex* inout)
- {
- fftw_one(m_plan, inout, NULL);
- return 0;
- };
-
- protected:
- fftw_plan m_plan;
-};
-
-FFTWCache::FFTWCache(const FFTWCacheId& id)
-: Cache<FFTWCacheId>(id)
-{
- m_plan = fftw_create_plan(id.m_n,
- (id.m_dir > 0 ? FFTW_FORWARD:FFTW_BACKWARD),
- FFTW_ESTIMATE | FFTW_IN_PLACE);
-
- if (m_plan == NULL) {
- goto fail;
- }
-
- return ;
-
-fail:
- throw std::bad_alloc();
-}
-
-FFTWCache::~FFTWCache()
-{
- fftw_destroy_plan(m_plan);
-}
-
-CacheManager<FFTWCacheId, FFTWCache> fftw_cmgr(10);
-
-void zfft_fftw(complex_double * inout, int n,
- int dir, int howmany, int normalize)
-{
- int i;
- complex_double *ptr = inout;
- FFTWCache* cache;
-
- cache = fftw_cmgr.get_cache(FFTWCacheId(n, dir));
-
- if (dir != -1 && dir != 1) {
- fprintf(stderr, "zfft: invalid dir=%d\n", dir);
- } else {
- for (i = 0; i < howmany; ++i, ptr += n) {
- cache->compute((fftw_complex *) ptr);
- }
- }
-
- if (normalize) {
- ptr = inout;
- for (i = n * howmany - 1; i >= 0; --i) {
- *((double *) (ptr)) /= n;
- *((double *) (ptr++) + 1) /= n;
- }
- }
-}
Deleted: branches/refactor_fft/scipy/fftpack/backends/fftw/zfftnd.cxx
===================================================================
--- branches/refactor_fft/scipy/fftpack/backends/fftw/zfftnd.cxx 2008-05-16 06:02:54 UTC (rev 4359)
+++ branches/refactor_fft/scipy/fftpack/backends/fftw/zfftnd.cxx 2008-05-16 06:03:43 UTC (rev 4360)
@@ -1,192 +0,0 @@
-/*
- * fftw2 backend for multi dimensional fft
- *
- * Original code by Pearu Peaterson
- *
- * Last Change: Tue May 13 02:00 PM 2008 J
- */
-#include <new>
-#include <cassert>
-
-#include <fftw.h>
-#include <rfftw.h>
-
-#include <cycliccache.h>
-#include "api.h"
-
-using namespace fft;
-
-static int equal_dims(int rank,int *dims1,int *dims2)
-{
- int i;
- for (i = 0; i < rank; ++i) {
- if (dims1[i] != dims2[i]) {
- return 0;
- }
- }
- return 1;
-}
-
-class NDFFTWCacheId {
- public:
- NDFFTWCacheId(int rank, int *dims, int dir, int flags);
- virtual ~NDFFTWCacheId();
-
- NDFFTWCacheId(const NDFFTWCacheId &);
-
- virtual bool operator==(const NDFFTWCacheId & other) const {
- return is_equal(other);
- };
-
- virtual bool is_equal(const NDFFTWCacheId & other) const;
-
- public:
- int m_dir;
- int m_rank;
- int *m_dims;
- int m_flags;
-
- private:
- int init(int rank, int *dims);
-};
-
-int NDFFTWCacheId::init(int rank, int *dims)
-{
- m_dims = (int *) malloc(sizeof(int) * rank);
- if (m_dims == NULL) {
- return -1;
- }
- memcpy(m_dims, dims, rank * sizeof(*m_dims));
-
- return 0;
-
-}
-
-NDFFTWCacheId::NDFFTWCacheId(int rank, int *dims, int dir, int flags) :
- m_dir(dir),
- m_rank(rank),
- m_flags(flags)
-{
- if (init(rank, dims)) {
- goto fail;
- }
-
-fail:
- std::bad_alloc();
-}
-
-NDFFTWCacheId::NDFFTWCacheId(const NDFFTWCacheId & copy) :
- m_dir(copy.m_dir),
- m_rank(copy.m_rank),
- m_flags(copy.m_flags)
-{
- if (init(copy.m_rank, copy.m_dims)) {
- goto fail;
- }
-
-fail:
- std::bad_alloc();
-}
-
-NDFFTWCacheId::~NDFFTWCacheId()
-{
- free(m_dims);
-}
-
-bool NDFFTWCacheId::is_equal(const NDFFTWCacheId & other) const
-{
- bool res;
-
- res = (other.m_dir == m_dir);
- res = res && (other.m_flags == m_flags);
-
- if (m_rank == other.m_rank) {
- res = res && equal_dims(m_rank, m_dims, other.m_dims);
- } else {
- return false;
- }
-
- return res;
-}
-
-class NDFFTWCache : public Cache < NDFFTWCacheId > {
- public:
- NDFFTWCache(const NDFFTWCacheId & id);
- virtual ~ NDFFTWCache();
-
- int compute(fftw_complex * inout) const
- {
- fftwnd_one(m_plan, inout, NULL);
- return 0;
- };
-
- protected:
- fftwnd_plan m_plan;
-};
-
-NDFFTWCache::NDFFTWCache(const NDFFTWCacheId & id)
-: Cache < NDFFTWCacheId > (id)
-{
- int flags = FFTW_ESTIMATE | FFTW_IN_PLACE;
-#if 0
- int sz;
- int i;
-
- sz = 1;
- for (i = 0; i < m_id.m_rank; ++i) {
- sz *= m_id.m_dims[i];
- }
-#endif
-
- m_plan = fftwnd_create_plan(m_id.m_rank,
- m_id.m_dims,
- (id.m_dir > 0 ?
- FFTW_FORWARD : FFTW_BACKWARD),
- flags);
-
- if (m_plan == NULL) {
- goto fail;
- }
-
- return;
-
-fail:
- throw std::bad_alloc();
-}
-
-NDFFTWCache::~NDFFTWCache()
-{
- fftwnd_destroy_plan(m_plan);
-}
-
-static CacheManager < NDFFTWCacheId, NDFFTWCache > fftwnd_cmgr(10);
-
-void zfftnd_fftw(complex_double * inout, int rank,
- int *dims, int direction, int howmany,
- int normalize)
-{
- int i, sz;
- complex_double *ptr = inout;
- NDFFTWCache *cache;
-
- sz = 1;
- for (i = 0; i < rank; ++i) {
- sz *= dims[i];
- }
-
- cache = fftwnd_cmgr.get_cache(
- NDFFTWCacheId(rank, dims, direction,
- FFTW_IN_PLACE | FFTW_ESTIMATE));
-
- for (i = 0; i < howmany; ++i, ptr += sz) {
- cache->compute((fftw_complex*)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