Made fourier library more flexible. Compatible with boost::multi_array now.

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 <cmath>
+#include <boost/multi_array.hpp>
 
 namespace CosmoTool
 {
   
+  namespace details {
+    static void no_free_euclidian_map(void *) {}
+  }
+  
   template<typename T>
   class EuclidianFourierMapBase: public FourierMap<T>
   {
@@ -61,6 +66,19 @@ namespace CosmoTool
         m_size *= m_dims[i];
     }
 
+    template<typename ArrayType>
+    EuclidianFourierMapBase(ArrayType& indata)
+    {
+      m_data = boost::shared_ptr<T>(
+            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<T, 1>& array() { 
+      return boost::multi_array_ref<T, 1>(m_data.get(), boost::extents[size]);
+    }
+
+    boost::const_multi_array_ref<T, 1>& array() const { 
+      return boost::const_multi_array_ref<T, 1>(m_data.get(), boost::extents[size]);
+    }
+    
     virtual FourierMap<T> *copy() const
     {
       FourierMap<T> *m = this->mimick();
@@ -90,6 +116,12 @@ namespace CosmoTool
       : EuclidianFourierMapBase<T>(indata, indims) 
     {}
 
+    template<typename ArrayType>
+    EuclidianFourierMapReal(ArrayType& indata) 
+      : EuclidianFourierMapBase<T>(indata) 
+    {}
+
+
     virtual FourierMap<T> *mimick() const
     {
       return new EuclidianFourierMapReal<T>(
@@ -125,18 +157,34 @@ namespace CosmoTool
 			       int dim0,
 			       const DimArray& indims, 
 			       const std::vector<double>& dk)
-      : EuclidianFourierMapBase<std::complex<T> >(indata, indims), delta_k(dk), m_dim0(dim0), even0((dim0 % 2)==0)
+      : EuclidianFourierMapBase<std::complex<T> >(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++)
+      for (int q = 1; q < indims.size(); q++) {
-	{
         plane_size *= indims[q];
         alleven = alleven && ((indims[q]%2)==0);
       }
     }
 
+    template<typename ArrayType>
+    EuclidianFourierMapComplex(ArrayType& indata, 
+                      	       int dim0,
+                               const std::vector<double>& dk) 
+      : EuclidianFourierMapBase<std::complex<T> >(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<std::complex<T> > *mimick() const
     {
       return 
@@ -162,8 +210,7 @@ namespace CosmoTool
       assert(kvec.size() == dims.size());
 
       kvec[0] = ik[0] * delta_k[0];
-      for (int q = 1; q < ik.size(); q++)
+      for (int q = 1; q < ik.size(); q++) {
-	{
         int dk = ik[q];
         if (dk > dims[q]/2)
           dk = dk - dims[q];
@@ -201,8 +248,7 @@ namespace CosmoTool
       double k2 = 0;
       k2 += CosmoTool::square(ik[0]*delta_k[0]);
 
-      for (int q = 1; q < ik.size(); q++)
+      for (int q = 1; q < ik.size(); q++) {
-	{
         int dk = ik[q];
 
         if (dk > dims[q]/2)
@@ -231,7 +277,8 @@ namespace CosmoTool
     virtual std::complex<T> dot_product(const FourierMap<std::complex<T> >& other) const
       throw(std::bad_cast)
     {
-      const EuclidianFourierMapComplex<T>& m2 = dynamic_cast<const EuclidianFourierMapComplex<T>&>(other);
+      const EuclidianFourierMapComplex<T>& m2 = 
+        dynamic_cast<const EuclidianFourierMapComplex<T>&>(other);
       if (this->size() != m2.size())
         throw std::bad_cast();
       
@@ -241,17 +288,14 @@ namespace CosmoTool
       int N0 = dims[0] + (even0 ? 0 : 1);
       std::complex<T> result = 0;
 
-      for (long q0 = 1; q0 < N0-1; q0++)
+      for (long q0 = 1; q0 < N0-1; q0++) {
-	{
+        for (long p = 0; p < plane_size; p++) {
-	  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)
+      if (even0) {
-	{
         for (long p = 0; p < plane_size; p++)
           {
             long q0 = N0*p, q1 = (p+1)*N0-1;

src/fourier/details/euclidian_transform.hpp

@@ -96,14 +96,18 @@ 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(Nc),
+                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(), &swapped_dims[0],
+        m_analysis = calls::plan_dft_r2c(
-                                        realMap->data(), (typename calls::complex_type *)fourierMap->data(),
+              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],
+        m_synthesis = calls::plan_dft_c2r(
+              dims.size(), &swapped_dims[0],
               (typename calls::complex_type *)fourierMap->data(), realMap->data(),
               FFTW_DESTROY_INPUT|FFTW_MEASURE);
       }

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) \