diff --git a/src/fourier/details/euclidian_maps.hpp b/src/fourier/details/euclidian_maps.hpp index 5b808b2..1e82a90 100644 --- a/src/fourier/details/euclidian_maps.hpp +++ b/src/fourier/details/euclidian_maps.hpp @@ -37,10 +37,15 @@ knowledge of the CeCILL license and that you accept its terms. #define __DETAILS_EUCLIDIAN_MAPS #include +#include namespace CosmoTool { + namespace details { + static void no_free_euclidian_map(void *) {} + } + template class EuclidianFourierMapBase: public FourierMap { @@ -58,9 +63,22 @@ namespace CosmoTool m_dims = indims; m_size = 1; for (int i = 0; i < m_dims.size(); i++) - m_size *= m_dims[i]; + m_size *= m_dims[i]; } + template + EuclidianFourierMapBase(ArrayType& indata) + { + m_data = boost::shared_ptr( + indata.origin(), + std::ptr_fun(details::no_free_euclidian_map)); + m_dims = DimArray(indata.num_dimensions()); + m_size = indata.num_elements(); + for (int i = 0; i < m_dims.size(); i++) + m_dims[i] = indata.shape()[i]; + } + + virtual ~EuclidianFourierMapBase() { } @@ -71,6 +89,14 @@ namespace CosmoTool virtual T *data() { return m_data.get(); } virtual long size() const { return m_size; } + boost::multi_array_ref& array() { + return boost::multi_array_ref(m_data.get(), boost::extents[size]); + } + + boost::const_multi_array_ref& array() const { + return boost::const_multi_array_ref(m_data.get(), boost::extents[size]); + } + virtual FourierMap *copy() const { FourierMap *m = this->mimick(); @@ -90,6 +116,12 @@ namespace CosmoTool : EuclidianFourierMapBase(indata, indims) {} + template + EuclidianFourierMapReal(ArrayType& indata) + : EuclidianFourierMapBase(indata) + {} + + virtual FourierMap *mimick() const { return new EuclidianFourierMapReal( @@ -125,16 +157,32 @@ namespace CosmoTool int dim0, const DimArray& indims, const std::vector& dk) - : EuclidianFourierMapBase >(indata, indims), delta_k(dk), m_dim0(dim0), even0((dim0 % 2)==0) + : EuclidianFourierMapBase >(indata, indims), + delta_k(dk), m_dim0(dim0), even0((dim0 % 2)==0) { assert(dk.size() == indims.size()); plane_size = 1; alleven = true; - for (int q = 1; q < indims.size(); q++) - { - plane_size *= indims[q]; - alleven = alleven && ((indims[q]%2)==0); - } + for (int q = 1; q < indims.size(); q++) { + plane_size *= indims[q]; + alleven = alleven && ((indims[q]%2)==0); + } + } + + template + EuclidianFourierMapComplex(ArrayType& indata, + int dim0, + const std::vector& dk) + : EuclidianFourierMapBase >(indata), + delta_k(dk), m_dim0(dim0), even0((dim0 % 2)==0) + { + assert(dk.size() == indata.num_dimensions()); + plane_size = 1; + alleven = true; + for (int q = 1; q < indims.size(); q++) { + plane_size *= m_dims[q]; + alleven = alleven && ((indims[q]%2)==0); + } } virtual FourierMap > *mimick() const @@ -162,14 +210,13 @@ namespace CosmoTool assert(kvec.size() == dims.size()); kvec[0] = ik[0] * delta_k[0]; - for (int q = 1; q < ik.size(); q++) - { - int dk = ik[q]; - if (dk > dims[q]/2) - dk = dk - dims[q]; + for (int q = 1; q < ik.size(); q++) { + int dk = ik[q]; + if (dk > dims[q]/2) + dk = dk - dims[q]; - kvec[q] = dk*delta_k[q]; - } + kvec[q] = dk*delta_k[q]; + } } template @@ -201,15 +248,14 @@ namespace CosmoTool double k2 = 0; k2 += CosmoTool::square(ik[0]*delta_k[0]); - for (int q = 1; q < ik.size(); q++) - { - int dk = ik[q]; + for (int q = 1; q < ik.size(); q++) { + int dk = ik[q]; - if (dk > dims[q]/2) - dk = dk - dims[q]; - - k2 += CosmoTool::square(delta_k[q]*dk); - } + if (dk > dims[q]/2) + dk = dk - dims[q]; + + k2 += CosmoTool::square(delta_k[q]*dk); + } return std::sqrt(k2); } @@ -231,7 +277,8 @@ namespace CosmoTool virtual std::complex dot_product(const FourierMap >& other) const throw(std::bad_cast) { - const EuclidianFourierMapComplex& m2 = dynamic_cast&>(other); + const EuclidianFourierMapComplex& m2 = + dynamic_cast&>(other); if (this->size() != m2.size()) throw std::bad_cast(); @@ -241,24 +288,21 @@ namespace CosmoTool int N0 = dims[0] + (even0 ? 0 : 1); std::complex result = 0; - for (long q0 = 1; q0 < N0-1; q0++) - { - for (long p = 0; p < plane_size; p++) - { - long idx = q0+dims[0]*p; - assert(idx < this->size()); - result += T(2)*(std::conj(d1[idx]) * d2[idx]).real(); - } - } - if (even0) - { - for (long p = 0; p < plane_size; p++) - { - long q0 = N0*p, q1 = (p+1)*N0-1; - result += T(2)*std::conj(d1[q0]) * d2[q0]; - result += T(2)*std::conj(d1[q1]) * d2[q1]; - } - } + for (long q0 = 1; q0 < N0-1; q0++) { + for (long p = 0; p < plane_size; p++) { + long idx = q0+dims[0]*p; + assert(idx < this->size()); + result += T(2)*(std::conj(d1[idx]) * d2[idx]).real(); + } + } + if (even0) { + for (long p = 0; p < plane_size; p++) + { + long q0 = N0*p, q1 = (p+1)*N0-1; + result += T(2)*std::conj(d1[q0]) * d2[q0]; + result += T(2)*std::conj(d1[q1]) * d2[q1]; + } + } return result; } diff --git a/src/fourier/details/euclidian_transform.hpp b/src/fourier/details/euclidian_transform.hpp index ac16fa9..4754361 100644 --- a/src/fourier/details/euclidian_transform.hpp +++ b/src/fourier/details/euclidian_transform.hpp @@ -96,16 +96,20 @@ namespace CosmoTool std::ptr_fun(calls::free)), m_dims); fourierMap = new EuclidianFourierMapComplex( - boost::shared_ptr >((std::complex*)calls::alloc_complex(Nc), - std::ptr_fun(calls::free)), + boost::shared_ptr >( + (std::complex*)calls::alloc_complex(Nc), + std::ptr_fun(calls::free)), dims[0], m_dims_hc, dk); { - m_analysis = calls::plan_dft_r2c(dims.size(), &swapped_dims[0], - realMap->data(), (typename calls::complex_type *)fourierMap->data(), - FFTW_DESTROY_INPUT|FFTW_MEASURE); - m_synthesis = calls::plan_dft_c2r(dims.size(), &swapped_dims[0], - (typename calls::complex_type *)fourierMap->data(), realMap->data(), - FFTW_DESTROY_INPUT|FFTW_MEASURE); + m_analysis = calls::plan_dft_r2c( + dims.size(), &swapped_dims[0], + realMap->data(), + (typename calls::complex_type *)fourierMap->data(), + FFTW_DESTROY_INPUT|FFTW_MEASURE); + m_synthesis = calls::plan_dft_c2r( + dims.size(), &swapped_dims[0], + (typename calls::complex_type *)fourierMap->data(), realMap->data(), + FFTW_DESTROY_INPUT|FFTW_MEASURE); } } diff --git a/src/fourier/fft/fftw_calls.hpp b/src/fourier/fft/fftw_calls.hpp index 7d1138f..b5f37db 100644 --- a/src/fourier/fft/fftw_calls.hpp +++ b/src/fourier/fft/fftw_calls.hpp @@ -68,6 +68,8 @@ public: \ static void free(void *p) { fftw_free(p); } \ \ static void execute(plan_type p) { prefix ## _execute(p); } \ + static void execute_r2c(plan_type p, real_type *in, complex_type *out) { prefix ## _execute_dft_r2c(p, in, out); } \ + static void execute_c2r(plan_type p, complex_type *in, real_type *out) { prefix ## _execute_dft_c2r(p, in, out); } \ static plan_type plan_dft_r2c_2d(int Nx, int Ny, \ real_type *in, complex_type *out, \ unsigned flags) \