From 5341c8ed5d4b5686296b9544c38c19aa0e8b8b32 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Sat, 6 Oct 2012 12:49:06 -0400 Subject: [PATCH] Added general fourier/real transform (Nd euclidian, sphere to come...) --- sample/testEskow.cpp | 2 + src/fourier/base_types.hpp | 95 ++++++++++++++++ src/fourier/euclidian.hpp | 200 +++++++++++++++++++++++++++++++++ src/fourier/fft/fftw_calls.hpp | 75 +++++++++++++ 4 files changed, 372 insertions(+) create mode 100644 src/fourier/base_types.hpp create mode 100644 src/fourier/euclidian.hpp create mode 100644 src/fourier/fft/fftw_calls.hpp diff --git a/sample/testEskow.cpp b/sample/testEskow.cpp index 2f173a6..60c1821 100644 --- a/sample/testEskow.cpp +++ b/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++) diff --git a/src/fourier/base_types.hpp b/src/fourier/base_types.hpp new file mode 100644 index 0000000..669e382 --- /dev/null +++ b/src/fourier/base_types.hpp @@ -0,0 +1,95 @@ +#ifndef __BASE_FOURIER_TYPES_HPP +#define __BASE_FOURIER_TYPES_HPP + +#include +#include +#include + +namespace CosmoTool +{ + template + class FourierMap + { + protected: + FourierMap() {} + + public: + typedef Eigen::Matrix VecType; + typedef Eigen::Map MapType; + typedef Eigen::Map 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 *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 *map2) + { + assert(size() == map2->size()); + MapType m(data(), size()); + MapType m2(map2->data(), map2->size()); + + eigen() += map2->eigen(); + } + + virtual FourierMap *copy() const = 0; + virtual FourierMap *mimick() const = 0; + }; + + template + class FourierTransform + { + protected: + FourierTransform() {} + public: + virtual ~FourierTransform() { } + + virtual const FourierMap >& fourierSpace() const = 0; + virtual FourierMap >& fourierSpace() = 0; + + virtual const FourierMap& realSpace() const = 0; + virtual FourierMap& realSpace() = 0; + + virtual FourierTransform *mimick() const = 0; + + virtual void analysis() = 0; + virtual void synthesis() = 0; + virtual void analysis_conjugate() = 0; + virtual void synthesis_conjugate() = 0; + }; + + +}; + +#endif diff --git a/src/fourier/euclidian.hpp b/src/fourier/euclidian.hpp new file mode 100644 index 0000000..0f40e1a --- /dev/null +++ b/src/fourier/euclidian.hpp @@ -0,0 +1,200 @@ +#ifndef __COSMOTOOL_FOURIER_EUCLIDIAN_HPP +#define __COSMOTOOL_FOURIER_EUCLIDIAN_HPP + +#include +#include +#include "base_types.hpp" +#include "fft/fftw_calls.hpp" + +namespace CosmoTool +{ + + template + class EuclidianFourierMap: public FourierMap + { + private: + boost::shared_ptr m_data; + std::vector m_dims; + long m_size; + public: + EuclidianFourierMap(boost::shared_ptr indata, std::vector 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 *copy() const + { + FourierMap *m = this->mimick(); + m->eigen() = this->eigen(); + return m; + } + + virtual FourierMap *mimick() const + { + return new EuclidianFourierMap( + boost::shared_ptr((T *)fftw_malloc(sizeof(T)*m_size), + std::ptr_fun(fftw_free)), + m_dims); + } + }; + + + template + class EuclidianFourierTransform: public FourierTransform + { + private: + typedef FFTW_Calls calls; + EuclidianFourierMap *realMap; + EuclidianFourierMap > *fourierMap; + typename calls::plan_type m_analysis, m_synthesis; + double volume; + long N; + std::vector m_dims; + std::vector m_L; + public: + EuclidianFourierTransform(const std::vector& dims, const std::vector& 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( + boost::shared_ptr(calls::alloc_real(N), + std::ptr_fun(calls::free)), + dims); + fourierMap = new EuclidianFourierMap >( + boost::shared_ptr >((std::complex*)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 >& fourierSpace() const + { + return *fourierMap; + } + + FourierMap >& fourierSpace() + { + return *fourierMap; + } + + const FourierMap& realSpace() const + { + return *realMap; + } + + FourierMap& realSpace() + { + return *realMap; + } + + FourierTransform *mimick() const + { + return new EuclidianFourierTransform(m_dims, m_L); + } + }; + + template + class EuclidianFourierTransform_2d: public EuclidianFourierTransform + { + private: + template + static std::vector make_2d_vector(T2 a, T2 b) + { + T2 arr[2] = { a, b}; + return std::vector(&arr[0], &arr[2]); + } + public: + EuclidianFourierTransform_2d(int Nx, int Ny, double Lx, double Ly) + : EuclidianFourierTransform(make_2d_vector(Nx, Ny), make_2d_vector(Lx, Ly)) + { + } + + virtual ~EuclidianFourierTransform_2d() {} + + }; + + template + class EuclidianFourierTransform_3d: public EuclidianFourierTransform + { + private: + template + static std::vector make_3d_vector(T2 a, T2 b, T2 c) + { + T2 arr[2] = { a, b, c}; + return std::vector(&arr[0], &arr[3]); + } + + public: + EuclidianFourierTransform_3d(int Nx, int Ny, int Nz, double Lx, double Ly, double Lz) + : EuclidianFourierTransform(make_3d_vector(Nx, Ny, Nz), make_3d_vector(Lx, Ly, Lz)) + { + } + + virtual ~EuclidianFourierTransform_3d() {} + }; + + + +}; + +#endif diff --git a/src/fourier/fft/fftw_calls.hpp b/src/fourier/fft/fftw_calls.hpp new file mode 100644 index 0000000..0106055 --- /dev/null +++ b/src/fourier/fft/fftw_calls.hpp @@ -0,0 +1,75 @@ +#ifndef __FFTW_UNIFIED_CALLS_HPP +#define __FFTW_UNIFIED_CALLS_HPP + +#include + +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 class FFTW_Calls {}; + + +#define FFTW_CALLS_BASE(rtype, prefix) \ + template<> \ +class FFTW_Calls { \ +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