From d8734b2a59901137fae3be527c9bd51554d92f4a Mon Sep 17 00:00:00 2001
From: Guilhem Lavaux <guilhem.lavaux@gmail.com>
Date: Sun, 11 Nov 2012 14:25:05 -0500
Subject: [PATCH] Fixed allocation and specification of dimensions to FFTW
 (must be swapped compared to CosmoTool convention)

---
 src/fourier/euclidian.hpp | 16 ++++++++++------
 1 file changed, 10 insertions(+), 6 deletions(-)

diff --git a/src/fourier/euclidian.hpp b/src/fourier/euclidian.hpp
index 56d6964..600ebbb 100644
--- a/src/fourier/euclidian.hpp
+++ b/src/fourier/euclidian.hpp
@@ -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()