added ei_fftw_impl

3 files changed, 214 insertions, 14 deletions

diff --git a/unsupported/Eigen/FFT b/unsupported/Eigen/FFT
index 31d8c74c5..03f8504a4 100644
--- a/unsupported/Eigen/FFT
+++ b/unsupported/Eigen/FFT
@@ -30,9 +30,8 @@
 #define DEFAULT_FFT_IMPL ei_kissfft_impl
 
 // FFTW: faster, GPL -- incompatible with Eigen in LGPL form, bigger code size
-#ifdef FFTW_PATIENT  // definition of FFTW_PATIENT indicates the caller has included fftw3.h, we can use FFTW routines
-// TODO 
-// #include "src/FFT/ei_fftw_impl.h"
+#ifdef FFTW_ESTIMATE  // definition of FFTW_ESTIMATE indicates the caller has included fftw3.h, we can use FFTW routines
+#include "src/FFT/ei_fftw_impl.h"
 // #define DEFAULT_FFT_IMPL ei_fftw_impl
 #endif
 
diff --git a/unsupported/Eigen/src/FFT/ei_fftw_impl.h b/unsupported/Eigen/src/FFT/ei_fftw_impl.h
new file mode 100644
index 000000000..d592bbb20
--- /dev/null
+++ b/unsupported/Eigen/src/FFT/ei_fftw_impl.h
@@ -0,0 +1,198 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. 
+//
+// Copyright (C) 2009 Mark Borgerding mark a borgerding net
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+namespace Eigen {
+  // FFTW uses non-const arguments
+  // so we must use ugly const_cast calls for all the args it uses
+  //
+  // This should be safe as long as 
+  // 1. we use FFTW_ESTIMATE for all our planning
+  //       see the FFTW docs section 4.3.2 "Planner Flags"
+  // 2. fftw_complex is compatible with std::complex
+  //    This assumes std::complex<T> layout is array of size 2 with real,imag
+  template <typename T> 
+  T * ei_fftw_cast(const T* p) 
+  { 
+      return const_cast<T*>( p); 
+  }
+
+  fftw_complex * ei_fftw_cast( const std::complex<double> * p) 
+  { 
+      return const_cast<fftw_complex*>( reinterpret_cast<const fftw_complex*>(p) ); 
+  }
+
+  fftwf_complex * ei_fftw_cast( const std::complex<float> * p) 
+  { 
+      return const_cast<fftwf_complex*>( reinterpret_cast<const fftwf_complex*>(p) ); 
+  }
+
+  fftwl_complex * ei_fftw_cast( const std::complex<long double> * p) 
+  { 
+      return const_cast<fftwl_complex*>( reinterpret_cast<const fftwl_complex*>(p) ); 
+  }
+
+  template <typename T> 
+  struct ei_fftw_plan {};
+
+  template <> 
+  struct ei_fftw_plan<float>
+  {
+      typedef float scalar_type;
+      typedef fftwf_complex complex_type;
+      fftwf_plan m_plan;
+      ei_fftw_plan() :m_plan(NULL) {}
+      ~ei_fftw_plan() {if (m_plan) fftwf_destroy_plan(m_plan);}
+
+      void fwd(complex_type * dst,complex_type * src,int nfft) {
+          if (m_plan==NULL) m_plan = fftwf_plan_dft_1d(nfft,src,dst, FFTW_FORWARD, FFTW_ESTIMATE);
+          fftwf_execute_dft( m_plan, src,dst);
+      }
+      void inv(complex_type * dst,complex_type * src,int nfft) {
+          if (m_plan==NULL) m_plan = fftwf_plan_dft_1d(nfft,src,dst, FFTW_BACKWARD , FFTW_ESTIMATE);
+          fftwf_execute_dft( m_plan, src,dst);
+      }
+      void fwd(complex_type * dst,scalar_type * src,int nfft) {
+          if (m_plan==NULL) m_plan = fftwf_plan_dft_r2c_1d(nfft,src,dst,FFTW_ESTIMATE);
+          fftwf_execute_dft_r2c( m_plan,src,dst);
+      }
+      void inv(scalar_type * dst,complex_type * src,int nfft) {
+          if (m_plan==NULL)
+              m_plan = fftwf_plan_dft_c2r_1d(nfft,src,dst,FFTW_ESTIMATE);
+          fftwf_execute_dft_c2r( m_plan, src,dst);
+      }
+  };
+  template <> 
+  struct ei_fftw_plan<double>
+  {
+      typedef double scalar_type;
+      typedef fftw_complex complex_type;
+      fftw_plan m_plan;
+      ei_fftw_plan() :m_plan(NULL) {}
+      ~ei_fftw_plan() {if (m_plan) fftw_destroy_plan(m_plan);}
+
+      void fwd(complex_type * dst,complex_type * src,int nfft) {
+          if (m_plan==NULL) m_plan = fftw_plan_dft_1d(nfft,src,dst, FFTW_FORWARD, FFTW_ESTIMATE);
+          fftw_execute_dft( m_plan, src,dst);
+      }
+      void inv(complex_type * dst,complex_type * src,int nfft) {
+          if (m_plan==NULL) m_plan = fftw_plan_dft_1d(nfft,src,dst, FFTW_BACKWARD , FFTW_ESTIMATE);
+          fftw_execute_dft( m_plan, src,dst);
+      }
+      void fwd(complex_type * dst,scalar_type * src,int nfft) {
+          if (m_plan==NULL) m_plan = fftw_plan_dft_r2c_1d(nfft,src,dst,FFTW_ESTIMATE);
+          fftw_execute_dft_r2c( m_plan,src,dst);
+      }
+      void inv(scalar_type * dst,complex_type * src,int nfft) {
+          if (m_plan==NULL)
+              m_plan = fftw_plan_dft_c2r_1d(nfft,src,dst,FFTW_ESTIMATE);
+          fftw_execute_dft_c2r( m_plan, src,dst);
+      }
+  };
+  template <> 
+  struct ei_fftw_plan<long double>
+  {
+      typedef long double scalar_type;
+      typedef fftwl_complex complex_type;
+      fftwl_plan m_plan;
+      ei_fftw_plan() :m_plan(NULL) {}
+      ~ei_fftw_plan() {if (m_plan) fftwl_destroy_plan(m_plan);}
+
+      void fwd(complex_type * dst,complex_type * src,int nfft) {
+          if (m_plan==NULL) m_plan = fftwl_plan_dft_1d(nfft,src,dst, FFTW_FORWARD, FFTW_ESTIMATE);
+          fftwl_execute_dft( m_plan, src,dst);
+      }
+      void inv(complex_type * dst,complex_type * src,int nfft) {
+          if (m_plan==NULL) m_plan = fftwl_plan_dft_1d(nfft,src,dst, FFTW_BACKWARD , FFTW_ESTIMATE);
+          fftwl_execute_dft( m_plan, src,dst);
+      }
+      void fwd(complex_type * dst,scalar_type * src,int nfft) {
+          if (m_plan==NULL) m_plan = fftwl_plan_dft_r2c_1d(nfft,src,dst,FFTW_ESTIMATE);
+          fftwl_execute_dft_r2c( m_plan,src,dst);
+      }
+      void inv(scalar_type * dst,complex_type * src,int nfft) {
+          if (m_plan==NULL)
+              m_plan = fftwl_plan_dft_c2r_1d(nfft,src,dst,FFTW_ESTIMATE);
+          fftwl_execute_dft_c2r( m_plan, src,dst);
+      }
+  };
+
+  template <typename _Scalar>
+  struct ei_fftw_impl
+  {
+      typedef _Scalar Scalar;
+      typedef std::complex<Scalar> Complex;
+
+      void clear() 
+      {
+        m_plans.clear();
+      }
+
+      void fwd( Complex * dst,const Complex *src,int nfft)
+      {
+        get_plan(nfft,false,dst,src).fwd(ei_fftw_cast(dst), ei_fftw_cast(src),nfft );
+      }
+
+      // real-to-complex forward FFT
+      void fwd( Complex * dst,const Scalar * src,int nfft) 
+      {
+          get_plan(nfft,false,dst,src).fwd(ei_fftw_cast(dst), ei_fftw_cast(src) ,nfft);
+          int nhbins=(nfft>>1)+1;
+          for (int k=nhbins;k < nfft; ++k )
+              dst[k] = conj(dst[nfft-k]);
+      }
+
+      // inverse complex-to-complex
+      void inv(Complex * dst,const Complex  *src,int nfft)
+      {
+        get_plan(nfft,true,dst,src).inv(ei_fftw_cast(dst), ei_fftw_cast(src),nfft );
+        // scaling
+        Scalar s = 1./nfft;
+        for (int k=0;k<nfft;++k)
+          dst[k] *= s;
+      }
+
+      // half-complex to scalar
+      void inv( Scalar * dst,const Complex * src,int nfft) 
+      {
+        get_plan(nfft,true,dst,src).inv(ei_fftw_cast(dst), ei_fftw_cast(src),nfft );
+        Scalar s = 1./nfft;
+        for (int k=0;k<nfft;++k)
+          dst[k] *= s;
+      }
+
+  private:
+      typedef ei_fftw_plan<Scalar> PlanData;
+      typedef std::map<int,PlanData> PlanMap;
+
+      PlanMap m_plans;
+
+      PlanData & get_plan(int nfft,bool inverse,void * dst,const void * src)
+      {
+          bool inplace = (dst==src);
+          bool aligned = ( (reinterpret_cast<size_t>(src)&15) | (reinterpret_cast<size_t>(dst)&15) ) == 0;
+          int key = (nfft<<3 ) | (inverse<<2) | (inplace<<1) | aligned;
+          return m_plans[key];
+      }
+  };
+}
diff --git a/unsupported/Eigen/src/FFT/ei_kissfft_impl.h b/unsupported/Eigen/src/FFT/ei_kissfft_impl.h
index 91fa5ca18..a84ac68a0 100644
--- a/unsupported/Eigen/src/FFT/ei_kissfft_impl.h
+++ b/unsupported/Eigen/src/FFT/ei_kissfft_impl.h
@@ -39,6 +39,7 @@ namespace Eigen {
       std::vector<Complex> m_twiddles;
       std::vector<int> m_stageRadix;
       std::vector<int> m_stageRemainder;
+      std::vector<Complex> m_scratchBuf;
       bool m_inverse;
 
       void make_twiddles(int nfft,bool inverse)
@@ -75,6 +76,8 @@ namespace Eigen {
           n /= p;
           m_stageRadix.push_back(p);
           m_stageRemainder.push_back(n);
+          if ( p > 5 )
+              m_scratchBuf.resize(p); // scratchbuf will be needed in bfly_generic
         }while(n>1);
       }
 
@@ -249,7 +252,7 @@ namespace Eigen {
         Complex * twiddles = &m_twiddles[0];
         Complex t;
         int Norig = m_twiddles.size();
-        Complex * scratchbuf = (Complex*)alloca(p*sizeof(Complex) );
+        Complex * scratchbuf = &m_scratchBuf[0];
 
         for ( u=0; u<m; ++u ) {
           k=u;
@@ -341,28 +344,28 @@ namespace Eigen {
       void inv( Scalar * dst,const Complex * src,int nfft) 
       {
         if (nfft&3) {
-          m_scratchBuf.resize(nfft);
-          inv(&m_scratchBuf[0],src,nfft);
+          m_tmpBuf.resize(nfft);
+          inv(&m_tmpBuf[0],src,nfft);
           for (int k=0;k<nfft;++k)
-            dst[k] = m_scratchBuf[k].real();
+            dst[k] = m_tmpBuf[k].real();
         }else{
           // optimized version for multiple of 4
           int ncfft = nfft>>1;
           int ncfft2 = nfft>>2;
           Complex * rtw = real_twiddles(ncfft2);
-          m_scratchBuf.resize(ncfft);
-          m_scratchBuf[0] = Complex( src[0].real() + src[ncfft].real(), src[0].real() - src[ncfft].real() );
+          m_tmpBuf.resize(ncfft);
+          m_tmpBuf[0] = Complex( src[0].real() + src[ncfft].real(), src[0].real() - src[ncfft].real() );
           for (int k = 1; k <= ncfft / 2; ++k) {
             Complex fk = src[k];
             Complex fnkc = conj(src[ncfft-k]);
             Complex fek = fk + fnkc;
             Complex tmp = fk - fnkc;
             Complex fok = tmp * conj(rtw[k-1]);
-            m_scratchBuf[k] = fek + fok;
-            m_scratchBuf[ncfft-k] = conj(fek - fok);
+            m_tmpBuf[k] = fek + fok;
+            m_tmpBuf[ncfft-k] = conj(fek - fok);
           }
-          scale(&m_scratchBuf[0], ncfft, Scalar(1)/nfft );
-          get_plan(ncfft,true).work(0, reinterpret_cast<Complex*>(dst), &m_scratchBuf[0], 1,1);
+          scale(&m_tmpBuf[0], ncfft, Scalar(1)/nfft );
+          get_plan(ncfft,true).work(0, reinterpret_cast<Complex*>(dst), &m_tmpBuf[0], 1,1);
         }
       }
 
@@ -372,7 +375,7 @@ namespace Eigen {
 
       PlanMap m_plans;
       std::map<int, std::vector<Complex> > m_realTwiddles;
-      std::vector<Complex> m_scratchBuf;
+      std::vector<Complex> m_tmpBuf;
 
       int PlanKey(int nfft,bool isinverse) const { return (nfft<<1) | isinverse; }