guilhem_lavaux/cosmotool
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Fixed allocation and specification of dimensions to FFTW (must be swapped compared to CosmoTool convention)

Browse Source
This commit is contained in:
Guilhem Lavaux 2012-11-11 14:25:05 -05:00
parent 596dd5d681
commit d8734b2a59
1 changed files with 10 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
src/fourier/euclidian.hpp
View File

@ -225,7 +225,7 @@ namespace CosmoTool
EuclidianFourierMapComplex<T> *fourierMap;
typename calls::plan_type m_analysis, m_synthesis;
double volume;
long N;
long N, Nc;
DimArray m_dims, m_dims_hc;
std::vector<double> m_L;
public:
@ -233,6 +233,7 @@ namespace CosmoTool
{
assert(L.size() == dims.size());
std::vector<double> dk(L.size());
std::vector<int> swapped_dims(dims.size());
m_dims = dims;
m_dims_hc = dims;
@ -240,12 +241,15 @@ namespace CosmoTool
m_L = L;
N = 1;
Nc = 1;
volume = 1;
for (int i = 0; i < dims.size(); i++)
{
N *= dims[i];
Nc *= m_dims_hc[i];
volume *= L[i];
dk[i] = 2*M_PI/L[i];
swapped_dims[dims.size()-1-i] = dims[i];
}
realMap = new EuclidianFourierMapReal<T>(
@ -253,15 +257,15 @@ namespace CosmoTool
std::ptr_fun(calls::free)),
m_dims);
fourierMap = new EuclidianFourierMapComplex<T>(
boost::shared_ptr<std::complex<T> >((std::complex<T>*)calls::alloc_complex(N),
boost::shared_ptr<std::complex<T> >((std::complex<T>*)calls::alloc_complex(Nc),
std::ptr_fun(calls::free)),
dims[0], m_dims_hc, dk);
m_analysis = calls::plan_dft_r2c(dims.size(), &dims[0],
m_analysis = calls::plan_dft_r2c(dims.size(), &swapped_dims[0],
realMap->data(), (typename calls::complex_type *)fourierMap->data(),
FFTW_MEASURE);
m_synthesis = calls::plan_dft_c2r(dims.size(), &dims[0],
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_MEASURE);
FFTW_DESTROY_INPUT|FFTW_MEASURE);
}
virtual ~EuclidianFourierTransform()