diff --git a/sample/test_fft_calls.cpp b/sample/test_fft_calls.cpp index 35eeb12..d633c14 100644 --- a/sample/test_fft_calls.cpp +++ b/sample/test_fft_calls.cpp @@ -1,12 +1,17 @@ +#include #include "fourier/euclidian.hpp" using namespace CosmoTool; +using namespace std; int main() { EuclidianFourierTransform_2d dft(128,128,1.0,1.0); + double volume = 128*128; dft.realSpace().eigen().setRandom(); dft.analysis(); + cout << "Map dot-product = " << dft.realSpace().dot_product(dft.realSpace()) << endl; + cout << "Fourier dot-product = " << dft.fourierSpace().dot_product(dft.fourierSpace()).real()*volume << endl; return 0; } diff --git a/src/fourier/euclidian.hpp b/src/fourier/euclidian.hpp index a29c623..fac27f6 100644 --- a/src/fourier/euclidian.hpp +++ b/src/fourier/euclidian.hpp @@ -179,7 +179,9 @@ namespace CosmoTool { for (long p = 0; p < plane_size; p++) { - result += 2*(conj(d1[q0+N0*p]) * d2[q0+N0*p]).real(); + long idx = q0+N0*p; + assert(idx < this->size()); + result += 2*(conj(d1[idx]) * d2[idx]).real(); } } for (long p = 0; p < plane_size; p++) @@ -205,7 +207,7 @@ namespace CosmoTool typename calls::plan_type m_analysis, m_synthesis; double volume; long N; - DimArray m_dims; + DimArray m_dims, m_dims_hc; std::vector m_L; public: EuclidianFourierTransform(const DimArray& dims, const std::vector& L) @@ -214,6 +216,8 @@ namespace CosmoTool std::vector dk(L.size()); m_dims = dims; + m_dims_hc = dims; + m_dims_hc[0] = dims[0]/2+1; m_L = L; N = 1; @@ -228,11 +232,11 @@ namespace CosmoTool realMap = new EuclidianFourierMapReal( boost::shared_ptr(calls::alloc_real(N), std::ptr_fun(calls::free)), - dims); + m_dims); fourierMap = new EuclidianFourierMapComplex( boost::shared_ptr >((std::complex*)calls::alloc_complex(N), std::ptr_fun(calls::free)), - dims, dk); + m_dims_hc, dk); m_analysis = calls::plan_dft_r2c(dims.size(), &dims[0], realMap->data(), (typename calls::complex_type *)fourierMap->data(), FFTW_MEASURE);