Added general fourier/real transform (Nd euclidian, sphere to come...)

sample/testEskow.cpp

@@ -50,6 +50,8 @@ int main(int argc, char **argv)
 
   chol.cholesky(M, M.N, norm_E);
 
+  cout << "norm_E = " << norm_E << endl;
+
   for (int i = 0; i < M.N; i++)
     {
       for (int j = 0; j < M.N; j++)

src/fourier/base_types.hpp

@@ -0,0 +1,95 @@
+#ifndef __BASE_FOURIER_TYPES_HPP
+#define __BASE_FOURIER_TYPES_HPP
+
+#include <string>
+#include <Eigen/Dense>
+#include <complex>
+
+namespace CosmoTool
+{
+  template<typename T>
+  class FourierMap
+  {
+  protected:
+    FourierMap() {}
+    
+  public:
+    typedef Eigen::Matrix<T, 1, Eigen::Dynamic> VecType;
+    typedef Eigen::Map<VecType, Eigen::Aligned> MapType;
+    typedef Eigen::Map<VecType const, Eigen::Aligned> ConstMapType;
+
+    virtual ~FourierMap() {}
+
+    virtual const T* data() const = 0;
+    virtual T* data()  = 0;
+
+    virtual long size() const = 0;
+
+    MapType eigen()
+    {
+      return MapType(data(), size());
+    }
+
+    ConstMapType eigen() const
+    {
+      return ConstMapType(data(), size());
+    }
+ 
+    void scale(const T& factor)
+    {
+      MapType m(data(), size());
+      m *= factor;
+    }
+
+    void scale(const FourierMap<T> *map2)
+    {
+      assert(size() == map2->size());
+      MapType m(data(), size());
+      MapType m2(map2->data(), map2->size());
+      m *= m2;
+    }
+
+    void add(const T&  factor)
+    {
+      eigen() += factor;
+    }
+
+    void add(const FourierMap<T> *map2)
+    {
+      assert(size() == map2->size());
+      MapType m(data(), size());
+      MapType m2(map2->data(), map2->size());
+
+      eigen() += map2->eigen();
+    }
+
+    virtual FourierMap<T> *copy() const = 0;
+    virtual FourierMap<T> *mimick() const = 0;
+  };
+
+  template<typename T>
+  class FourierTransform
+  {
+  protected:
+    FourierTransform() {}
+  public:
+    virtual ~FourierTransform() { }
+    
+    virtual const FourierMap<std::complex<T> >& fourierSpace() const = 0;
+    virtual FourierMap<std::complex<T> >& fourierSpace() = 0;
+
+    virtual const FourierMap<T>& realSpace() const = 0;
+    virtual FourierMap<T>& realSpace() = 0;
+
+    virtual FourierTransform<T> *mimick() const = 0;
+
+    virtual void analysis() = 0;
+    virtual void synthesis() = 0;
+    virtual void analysis_conjugate() = 0;
+    virtual void synthesis_conjugate() = 0;
+  };
+
+
+};
+
+#endif

src/fourier/euclidian.hpp

@@ -0,0 +1,200 @@
+#ifndef __COSMOTOOL_FOURIER_EUCLIDIAN_HPP
+#define __COSMOTOOL_FOURIER_EUCLIDIAN_HPP
+
+#include <vector>
+#include <boost/shared_ptr.hpp>
+#include "base_types.hpp"
+#include "fft/fftw_calls.hpp"
+
+namespace CosmoTool
+{
+  
+  template<typename T>
+  class EuclidianFourierMap: public FourierMap<T>
+  {
+  private:
+    boost::shared_ptr<T> m_data;
+    std::vector<int> m_dims;
+    long m_size;
+  public:
+    EuclidianFourierMap(boost::shared_ptr<T> indata, std::vector<int> indims)
+    {
+      m_data = indata;
+      m_dims = indims;
+      m_size = 1;
+      for (int i = 0; i < m_dims.size(); i++)
+	m_size *= m_dims[i];
+    }
+
+    virtual ~EuclidianFourierMap()
+    {
+    }
+
+    virtual const T *data() const { return m_data.get(); }
+    virtual T *data() { return m_data.get(); }
+    virtual long size() const  { return m_size; } 
+    
+    virtual FourierMap<T> *copy() const
+    {
+      FourierMap<T> *m = this->mimick();
+      m->eigen() = this->eigen();
+      return m;
+    }
+
+    virtual FourierMap<T> *mimick() const
+    {
+      return new EuclidianFourierMap<T>(
+		      boost::shared_ptr<T>((T *)fftw_malloc(sizeof(T)*m_size), 
+					   std::ptr_fun(fftw_free)),
+		       m_dims);
+    }
+  };
+
+  
+  template<typename T>
+  class EuclidianFourierTransform: public FourierTransform<T>
+  {
+  private:
+    typedef FFTW_Calls<T> calls;
+    EuclidianFourierMap<T> *realMap;
+    EuclidianFourierMap<std::complex<T> > *fourierMap;
+    typename calls::plan_type m_analysis, m_synthesis;
+    double volume;
+    long N;
+    std::vector<int> m_dims;
+    std::vector<double> m_L;
+  public:
+    EuclidianFourierTransform(const std::vector<int>& dims, const std::vector<double>& L)
+    {
+      assert(L.size() == dims.size());
+
+      m_dims = dims;
+      m_L = L;
+
+      N = 1;
+      volume = 1;
+      for (int i = 0; i < dims.size(); i++)
+	{
+	  N *= dims[i];
+	  volume *= L[i];
+	}
+
+      realMap = new EuclidianFourierMap<T>(
+		   boost::shared_ptr<T>(calls::alloc_real(N),
+					std::ptr_fun(calls::free)),
+		   dims);
+      fourierMap = new EuclidianFourierMap<std::complex<T> >(
+		   boost::shared_ptr<std::complex<T> >((std::complex<T>*)calls::alloc_complex(N),
+						       std::ptr_fun(calls::free)), 
+		   dims);
+      m_analysis = calls::plan_dft_r2c(dims.size(), &dims[0],
+				       realMap->data(), (typename calls::complex_type *)fourierMap->data(),
+				       FFTW_MEASURE);
+      m_synthesis = calls::plan_dft_c2r(dims.size(), &dims[0],
+					(typename calls::complex_type *)fourierMap->data(), realMap->data(),
+					FFTW_MEASURE);
+    }
+    
+    virtual ~EuclidianFourierTransform()
+    {
+      delete realMap;
+      delete fourierMap;
+      calls::destroy_plan(m_synthesis);
+      calls::destroy_plan(m_analysis);
+    }
+
+    void synthesis()
+    {
+      calls::execute(m_synthesis);
+      realMap->scale(1/volume);
+    }
+    
+    void analysis()
+    {
+      calls::execute(m_analysis);
+      fourierMap->scale(volume/N);
+    }
+    
+    void synthesis_conjugate()
+    {
+      calls::execute(m_analysis);
+      fourierMap->scale(1/volume);
+    }
+    
+    void analysis_conjugate()
+    {
+      calls::execute(m_synthesis);
+      realMap->scale(volume/N);
+    }
+
+    const FourierMap<std::complex<T> >& fourierSpace() const
+    {
+      return *fourierMap;
+    }
+
+    FourierMap<std::complex<T> >& fourierSpace()
+    {
+      return *fourierMap;
+    }
+
+    const FourierMap<T>& realSpace() const
+    {
+      return *realMap;
+    }
+
+    FourierMap<T>& realSpace()
+    {
+      return *realMap;
+    }
+
+    FourierTransform<T> *mimick() const
+    {
+      return new EuclidianFourierTransform(m_dims, m_L);
+    }
+  };
+
+  template<typename T>
+  class EuclidianFourierTransform_2d: public EuclidianFourierTransform<T>
+  {
+  private:
+    template<typename T2>
+    static std::vector<T2> make_2d_vector(T2 a, T2 b)
+    {
+      T2 arr[2] = { a, b};
+      return std::vector<T2>(&arr[0], &arr[2]);
+    }
+  public:
+    EuclidianFourierTransform_2d(int Nx, int Ny, double Lx, double Ly)
+      : EuclidianFourierTransform<T>(make_2d_vector<int>(Nx, Ny), make_2d_vector<double>(Lx, Ly))
+    {
+    }
+
+    virtual ~EuclidianFourierTransform_2d() {}
+
+  };
+
+  template<typename T>
+  class EuclidianFourierTransform_3d: public EuclidianFourierTransform<T>
+  {
+  private:
+    template<typename T2>
+    static std::vector<T2> make_3d_vector(T2 a, T2 b, T2 c)
+    {
+      T2 arr[2] = { a, b, c};
+      return std::vector<T2>(&arr[0], &arr[3]);
+    }
+
+  public:
+    EuclidianFourierTransform_3d(int Nx, int Ny, int Nz, double Lx, double Ly, double Lz)
+      : EuclidianFourierTransform<T>(make_3d_vector<int>(Nx, Ny, Nz), make_3d_vector<double>(Lx, Ly, Lz))
+    {
+    }
+
+    virtual ~EuclidianFourierTransform_3d() {}
+  };
+
+
+
+};
+
+#endif

src/fourier/fft/fftw_calls.hpp

@@ -0,0 +1,75 @@
+#ifndef __FFTW_UNIFIED_CALLS_HPP
+#define __FFTW_UNIFIED_CALLS_HPP
+
+#include <fftw3.h>
+
+namespace CosmoTool
+{
+
+static inline void init_fftw_wisdom()
+{
+  fftw_import_system_wisdom();
+  fftw_import_wisdom_from_filename("fft_wisdom");
+}
+
+static inline void save_fftw_wisdom()
+{
+  fftw_export_wisdom_to_filename("fft_wisdom");
+}
+
+template<typename T> class FFTW_Calls {};
+
+
+#define FFTW_CALLS_BASE(rtype, prefix) \
+  template<>				\
+class FFTW_Calls<rtype> {		\
+public: \
+  typedef rtype real_type; \
+  typedef prefix ## _complex complex_type; \
+  typedef prefix ## _plan plan_type; \
+  \
+  static complex_type *alloc_complex(int N) { return prefix ## _alloc_complex(N); } \
+  static real_type *alloc_real(int N) { return prefix ## _alloc_real(N); } \
+  static void free(void *p) { fftw_free(p); } \
+\
+  static void execute(plan_type p) { prefix ## _execute(p); } \
+  static plan_type plan_dft_r2c_2d(int Nx, int Ny,  \
+		       real_type *in, complex_type *out, \
+		       unsigned flags) \
+  { \
+    return prefix ## _plan_dft_r2c_2d(Nx, Ny, in, out, \
+				flags); \
+  } \
+  static plan_type plan_dft_c2r_2d(int Nx, int Ny,  \
+                       complex_type *in, real_type *out, \
+                       unsigned flags) \
+  { \
+    return prefix ## _plan_dft_c2r_2d(Nx, Ny, in, out, \
+                                flags); \
+  } \
+  static plan_type plan_dft_r2c_3d(int Nx, int Ny, int Nz, \
+                                   real_type *in, complex_type *out, \
+                                   unsigned flags) \
+  { \
+    return prefix ## _plan_dft_r2c_3d(Nx, Ny, Nz, in, out, flags); \
+  } \
+  static plan_type plan_dft_r2c(int rank, const int *n, real_type *in, \
+                                complex_type *out, unsigned flags) \
+  { \
+    return prefix ## _plan_dft_r2c(rank, n, in, out, flags); \
+  } \
+  static plan_type plan_dft_c2r(int rank, const int *n, complex_type *in, \
+                                real_type *out, unsigned flags) \
+  { \
+    return prefix ## _plan_dft_c2r(rank, n, in, out, flags); \
+  } \
+  static void destroy_plan(plan_type plan) { prefix ## _destroy_plan(plan); } \
+}
+
+
+FFTW_CALLS_BASE(double, fftw);
+FFTW_CALLS_BASE(float, fftwf);
+
+};
+
+#endif