inlining,all namespace declaration moved to FFT, removed preprocessor definitions,

unsupported/Eigen/FFT

@@ -25,29 +25,39 @@
 #ifndef EIGEN_FFT_H
 #define EIGEN_FFT_H
 
-// ei_kissfft_impl:  small, free, reasonably efficient default, derived from kissfft
+#include <complex>
-#include "src/FFT/ei_kissfft_impl.h"
+#include <vector>
-#define DEFAULT_FFT_IMPL ei_kissfft_impl
+#include <map>
 
+#ifdef EIGEN_FFTW_DEFAULT
 // FFTW: faster, GPL -- incompatible with Eigen in LGPL form, bigger code size
-#ifdef FFTW_ESTIMATE  // definition of FFTW_ESTIMATE indicates the caller has included fftw3.h, we can use FFTW routines
+#  include <fftw3.h>
-#include "src/FFT/ei_fftw_impl.h"
+   namespace Eigen {
-#undef DEFAULT_FFT_IMPL
+#    include "src/FFT/ei_fftw_impl.h"
-#define DEFAULT_FFT_IMPL ei_fftw_impl
+     //template <typename T> typedef struct ei_fftw_impl  default_fft_impl; this does not work
-#endif
+     template <typename T> struct default_fft_impl : public ei_fftw_impl<T> {};
-
+   }
-// intel Math Kernel Library: fastest, commercial -- incompatible with Eigen in GPL form
+#elif defined EIGEN_MKL_DEFAULT
-#ifdef _MKL_DFTI_H_ // mkl_dfti.h has been included, we can use MKL FFT routines
 // TODO 
-// #include "src/FFT/ei_imkl_impl.h"
+// intel Math Kernel Library: fastest, commercial -- may be incompatible with Eigen in GPL form
-// #define DEFAULT_FFT_IMPL ei_imkl_impl
+   namespace Eigen {
+#    include "src/FFT/ei_imklfft_impl.h"
+     template <typename T> struct default_fft_impl : public ei_imklfft_impl {};
+   }
+#else
+// ei_kissfft_impl:  small, free, reasonably efficient default, derived from kissfft
+//
+  namespace Eigen {
+#   include "src/FFT/ei_kissfft_impl.h"
+     template <typename T> 
+       struct default_fft_impl : public ei_kissfft_impl<T> {};
+  }
 #endif
 
 namespace Eigen {
 
 template <typename _Scalar,
-         typename _Impl=DEFAULT_FFT_IMPL<_Scalar> 
+         typename _Impl=default_fft_impl<_Scalar> >
-         >
 class FFT
 {
   public:
@@ -120,7 +130,6 @@ class FFT
   private:
     impl_type m_impl;
 };
-#undef DEFAULT_FFT_IMPL
 }
 #endif
 /* vim: set filetype=cpp et sw=2 ts=2 ai: */

unsupported/Eigen/src/FFT/ei_fftw_impl.h

@@ -22,7 +22,8 @@
 // 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
   //
@@ -32,21 +33,25 @@ namespace Eigen {
   // 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> 
+  inline 
   T * ei_fftw_cast(const T* p) 
   { 
       return const_cast<T*>( p); 
   }
 
+  inline 
   fftw_complex * ei_fftw_cast( const std::complex<double> * p) 
   {
       return const_cast<fftw_complex*>( reinterpret_cast<const fftw_complex*>(p) ); 
   }
 
+  inline 
   fftwf_complex * ei_fftw_cast( const std::complex<float> * p) 
   { 
       return const_cast<fftwf_complex*>( reinterpret_cast<const fftwf_complex*>(p) ); 
   }
 
+  inline 
   fftwl_complex * ei_fftw_cast( const std::complex<long double> * p) 
   { 
       return const_cast<fftwl_complex*>( reinterpret_cast<const fftwl_complex*>(p) ); 
@@ -64,18 +69,22 @@ namespace Eigen {
       ei_fftw_plan() :m_plan(NULL) {}
       ~ei_fftw_plan() {if (m_plan) fftwf_destroy_plan(m_plan);}
 
+      inline
       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);
       }
+      inline
       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);
       }
+      inline
       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);
       }
+      inline
       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);
@@ -91,18 +100,22 @@ namespace Eigen {
       ei_fftw_plan() :m_plan(NULL) {}
       ~ei_fftw_plan() {if (m_plan) fftw_destroy_plan(m_plan);}
 
+      inline
       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);
       }
+      inline
       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);
       }
+      inline
       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);
       }
+      inline
       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);
@@ -118,18 +131,22 @@ namespace Eigen {
       ei_fftw_plan() :m_plan(NULL) {}
       ~ei_fftw_plan() {if (m_plan) fftwl_destroy_plan(m_plan);}
 
+      inline
       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);
       }
+      inline
       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);
       }
+      inline
       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);
       }
+      inline
       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);
@@ -143,17 +160,20 @@ namespace Eigen {
       typedef _Scalar Scalar;
       typedef std::complex<Scalar> Complex;
 
+      inline
       void clear() 
       {
         m_plans.clear();
       }
 
+      inline
       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
+      inline
       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);
@@ -163,30 +183,37 @@ namespace Eigen {
       }
 
       // inverse complex-to-complex
+      inline
       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 );
+
+        //TODO move scaling to Eigen::FFT
         // scaling
-        Scalar s = 1./nfft;
+        Scalar s = Scalar(1.)/nfft;
         for (int k=0;k<nfft;++k)
           dst[k] *= s;
       }
 
       // half-complex to scalar
+      inline
       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;
+
+        //TODO move scaling to Eigen::FFT
+        Scalar s = Scalar(1.)/nfft;
         for (int k=0;k<nfft;++k)
           dst[k] *= s;
       }
 
-  private:
+  protected:
       typedef ei_fftw_plan<Scalar> PlanData;
       typedef std::map<int,PlanData> PlanMap;
 
       PlanMap m_plans;
 
+      inline
       PlanData & get_plan(int nfft,bool inverse,void * dst,const void * src)
       {
           bool inplace = (dst==src);
@@ -195,4 +222,3 @@ namespace Eigen {
           return m_plans[key];
       }
   };
-}

unsupported/Eigen/src/FFT/ei_kissfft_impl.h

@@ -22,11 +22,7 @@
 // License and a copy of the GNU General Public License along with
 // Eigen. If not, see <http://www.gnu.org/licenses/>.
 
-#include <complex>
-#include <vector>
-#include <map>
 
-namespace Eigen {
 
   // This FFT implementation was derived from kissfft http:sourceforge.net/projects/kissfft
   // Copyright 2003-2009 Mark Borgerding
@@ -51,13 +47,6 @@ namespace Eigen {
           m_twiddles[i] = exp( Complex(0,i*phinc) );
       }
 
-      void conjugate()
-      {
-        m_inverse = !m_inverse;
-        for ( size_t i=0;i<m_twiddles.size() ;++i)
-          m_twiddles[i] = conj( m_twiddles[i] );
-      }
-
       void factorize(int nfft)
       {
         //start factoring out 4's, then 2's, then 3,5,7,9,...
@@ -116,6 +105,7 @@ namespace Eigen {
           }
         }
 
+      inline
       void bfly2( Complex * Fout, const size_t fstride, int m)
       {
         for (int k=0;k<m;++k) {
@@ -125,6 +115,7 @@ namespace Eigen {
         }
       }
 
+      inline
       void bfly4( Complex * Fout, const size_t fstride, const size_t m)
       {
         Complex scratch[6];
@@ -147,6 +138,7 @@ namespace Eigen {
         }
       }
 
+      inline
       void bfly3( Complex * Fout, const size_t fstride, const size_t m)
       {
         size_t k=m;
@@ -175,6 +167,7 @@ namespace Eigen {
         }while(--k);
       }
 
+      inline
       void bfly5( Complex * Fout, const size_t fstride, const size_t m)
       {
         Complex *Fout0,*Fout1,*Fout2,*Fout3,*Fout4;
@@ -241,6 +234,7 @@ namespace Eigen {
       }
 
       /* perform the butterfly for one stage of a mixed radix FFT */
+      inline
       void bfly_generic(
           Complex * Fout,
           const size_t fstride,
@@ -290,6 +284,7 @@ namespace Eigen {
       }
 
       template <typename _Src>
+      inline
         void fwd( Complex * dst,const _Src *src,int nfft)
         {
           get_plan(nfft,false).work(0, dst, src, 1,1);
@@ -299,6 +294,7 @@ namespace Eigen {
       // perform two FFTs of src even and src odd
       // then twiddle to recombine them into the half-spectrum format
       // then fill in the conjugate symmetric half
+      inline
       void fwd( Complex * dst,const Scalar * src,int nfft) 
       {
         if ( nfft&3  ) {
@@ -334,6 +330,7 @@ namespace Eigen {
       }
 
       // inverse complex-to-complex
+      inline
       void inv(Complex * dst,const Complex  *src,int nfft)
       {
         get_plan(nfft,true).work(0, dst, src, 1,1);
@@ -341,6 +338,7 @@ namespace Eigen {
       }
 
       // half-complex to scalar
+      inline
       void inv( Scalar * dst,const Complex * src,int nfft) 
       {
         if (nfft&3) {
@@ -369,7 +367,7 @@ namespace Eigen {
         }
       }
 
-      private:
+      protected:
       typedef ei_kiss_cpx_fft<Scalar> PlanData;
       typedef std::map<int,PlanData> PlanMap;
 
@@ -377,8 +375,10 @@ namespace Eigen {
       std::map<int, std::vector<Complex> > m_realTwiddles;
       std::vector<Complex> m_tmpBuf;
 
+      inline
       int PlanKey(int nfft,bool isinverse) const { return (nfft<<1) | isinverse; }
 
+      inline
       PlanData & get_plan(int nfft,bool inverse)
       {
         // TODO look for PlanKey(nfft, ! inverse) and conjugate the twiddles
@@ -390,6 +390,7 @@ namespace Eigen {
         return pd;
       }
 
+      inline
       Complex * real_twiddles(int ncfft2)
       {
         std::vector<Complex> & twidref = m_realTwiddles[ncfft2];// creates new if not there
@@ -403,10 +404,11 @@ namespace Eigen {
         return &twidref[0];
       }
 
+      // TODO move scaling up into Eigen::FFT
+      inline
       void scale(Complex *dst,int n,Scalar s) 
       {
         for (int k=0;k<n;++k)
           dst[k] *= s;
       }
     };
-}

unsupported/test/CMakeLists.txt

@@ -23,6 +23,6 @@ ei_add_test(FFT)
 
 find_package(FFTW)
 if(FFTW_FOUND)
-  ei_add_test(FFTW  " " "-lfftw3 -lfftw3f -lfftw3l" )
+  ei_add_test(FFTW  "-DEIGEN_FFTW_DEFAULT " "-lfftw3 -lfftw3f -lfftw3l" )
 endif(FFTW_FOUND)