Merge branch 'master' of bitbucket.org:glavaux/cosmotool

src/fourier/fft/fftw_calls.hpp

@@ -113,7 +113,7 @@ public: \
 
 FFTW_CALLS_BASE(double, fftw);
 FFTW_CALLS_BASE(float, fftwf);
-
+#undef FFTW_CALLS_BASE
 };
 
 #endif

src/fourier/fft/fftw_calls_mpi.hpp

@@ -0,0 +1,98 @@
+#ifndef __MPI_FFTW_UNIFIED_CALLS_HPP
+#define __MPI_FFTW_UNIFIED_CALLS_HPP
+
+#include <mpi.h>
+#include <fftw3-mpi.h>
+
+namespace CosmoTool
+{
+
+static inline void init_fftw_mpi()
+{
+  fftw_mpi_init();
+}
+
+static inline void done_fftw_mpi()
+{
+  fftw_mpi_cleanup();
+}
+
+template<typename T> class FFTW_MPI_Calls {};
+
+
+#define FFTW_MPI_CALLS_BASE(rtype, prefix) \
+  template<>				\
+class FFTW_MPI_Calls<rtype> {		\
+public: \
+  typedef rtype real_type; \
+  typedef prefix ## _complex complex_type; \
+  typedef prefix ## _plan plan_type; \
+  \
+  static complex_type *alloc_complex(size_t N) { return prefix ## _alloc_complex(N); } \
+  static real_type *alloc_real(size_t N) { return prefix ## _alloc_real(N); } \
+  static void free(void *p) { fftw_free(p); } \
+\
+  static ptrdiff_t local_size_2d(ptrdiff_t N0, ptrdiff_t N1, MPI_Comm comm, \
+                                ptrdiff_t *local_n0, ptrdiff_t *local_0_start) { \
+    return prefix ## _mpi_local_size_2d(N0, N1, comm, local_n0, local_0_start); \
+  } \
+\
+  static ptrdiff_t local_size_3d(ptrdiff_t N0, ptrdiff_t N1, ptrdiff_t N2, MPI_Comm comm, \
+                                ptrdiff_t *local_n0, ptrdiff_t *local_0_start) { \
+    return prefix ## _mpi_local_size_3d(N0, N1, N2, comm, local_n0, local_0_start); \
+  } \
+\
+  static void execute(plan_type p) { prefix ## _execute(p); } \
+  static void execute_r2c(plan_type p, real_type *in, complex_type *out) { prefix ## _mpi_execute_dft_r2c(p, in, out); } \
+  static void execute_c2r(plan_type p, complex_type *in, real_type *out) { prefix ## _mpi_execute_dft_c2r(p, in, out); } \
+\
+  static plan_type plan_dft_r2c_2d(int Nx, int Ny,  \
+		       real_type *in, complex_type *out, \
+		       MPI_Comm comm, unsigned flags) \
+  { \
+    return prefix ## _mpi_plan_dft_r2c_2d(Nx, Ny, in, out, \
+				comm, flags); \
+  } \
+  static plan_type plan_dft_c2r_2d(int Nx, int Ny,  \
+                       complex_type *in, real_type *out, \
+                       MPI_Comm comm, unsigned flags) \
+  { \
+    return prefix ## _mpi_plan_dft_c2r_2d(Nx, Ny, in, out, \
+                                comm, flags); \
+  } \
+  static plan_type plan_dft_r2c_3d(int Nx, int Ny, int Nz, \
+                                   real_type *in, complex_type *out, \
+                                   MPI_Comm comm, unsigned flags) \
+  { \
+    return prefix ## _mpi_plan_dft_r2c_3d(Nx, Ny, Nz, in, out, comm, flags); \
+  } \
+  static plan_type plan_dft_c2r_3d(int Nx, int Ny, int Nz, \
+                                   complex_type *in, real_type *out, \
+                                   MPI_Comm comm, \
+                                   unsigned flags) \
+  { \
+    return prefix ## _mpi_plan_dft_c2r_3d(Nx, Ny, Nz, in, out, comm, flags); \
+  } \
+\
+  static plan_type plan_dft_r2c(int rank, const ptrdiff_t *n, real_type *in, \
+                                complex_type *out, MPI_Comm comm, unsigned flags) \
+  { \
+    return prefix ## _mpi_plan_dft_r2c(rank, n, in, out, comm, flags); \
+  } \
+  static plan_type plan_dft_c2r(int rank, const ptrdiff_t *n, complex_type *in, \
+                                real_type *out, MPI_Comm comm, unsigned flags) \
+  { \
+    return prefix ## _mpi_plan_dft_c2r(rank, n, in, out, comm, flags); \
+  } \
+  static void destroy_plan(plan_type plan) { prefix ## _destroy_plan(plan); } \
+}
+
+
+FFTW_MPI_CALLS_BASE(double, fftw);
+FFTW_MPI_CALLS_BASE(float, fftwf);
+
+#undef FFTW_MPI_CALLS_BASE
+
+};
+
+#endif

src/hdf5_array.hpp

@@ -21,7 +21,7 @@ liability.
 In this respect, the user's attention is drawn to the risks associated
 with loading,  using,  modifying and/or developing or reproducing the
 software by the user in light of its specific status of free software,
-that may mean  that it is complicated to manipulate,  and  that  also
+sthat may mean  that it is complicated to manipulate,  and  that  also
 therefore means  that it is reserved for developers  and  experienced
 professionals having in-depth computer knowledge. Users are therefore
 encouraged to load and test the software's suitability as regards their
@@ -90,26 +90,70 @@ namespace CosmoTool {
   
   #undef HDF5_TYPE
 
+  // Extent generator
+  template<std::size_t r>
+  struct hdf5_extent_gen {
+    typedef typename boost::detail::multi_array::extent_gen<r> type;
+    
+    static inline type build(hsize_t *d)
+    {
+      return (hdf5_extent_gen<r-1>::build(d))[d[r-1]];
+    }
+  };
+
+  template<>
+  struct hdf5_extent_gen<0> {    
+    static inline boost::multi_array_types::extent_gen build(hsize_t *d)
+    {
+      return boost::extents;
+    }
+  };
+  
+
 //!_______________________________________________________________________________________
 //!     
 //!     write_hdf5 multi_array
 //!         
+//!     \author Guilhem Lavaux (2014-2015)
 //!     \author leo Goodstadt (04 March 2013)
 //!     
 //!_______________________________________________________________________________________
   template<typename ArrayType,  typename hdf5_data_type>
   void hdf5_write_array(H5::CommonFG& fg, const std::string& data_set_name, 
                         const ArrayType& data, 
-                        const hdf5_data_type& datatype)
+                        const hdf5_data_type& datatype,
+                        const std::vector<hsize_t>& dimensions,
+                        bool doCreate = true,
+                        bool useBases = false)
   {
-      std::vector<hsize_t> dimensions(data.shape(), data.shape() + data.num_dimensions());
       H5::DataSpace dataspace(data.num_dimensions(), dimensions.data());
 
-      H5::DataSet dataset = fg.createDataSet(data_set_name, datatype, dataspace);
+      if (useBases) {
+        std::vector<hssize_t> offsets(data.index_bases(), data.index_bases() + data.num_dimensions());
+        dataspace.offsetSimple(offsets.data());
+      }
+
+      H5::DataSet dataset;
+      if (doCreate)
+        dataset = fg.createDataSet(data_set_name, datatype, dataspace);
+      else
+        dataset = fg.openDataSet(data_set_name);
 
       dataset.write(data.data(), datatype);
   }
 
+
+  template<typename ArrayType,  typename hdf5_data_type>
+  void hdf5_write_array(H5::CommonFG& fg, const std::string& data_set_name, 
+                        const ArrayType& data, 
+                        const hdf5_data_type& datatype,
+                        bool doCreate = true,
+                        bool useBases = false)
+  {
+      std::vector<hsize_t> dimensions(data.shape(), data.shape() + data.num_dimensions());
+      hdf5_write_array(fg, data_set_name, data, datatype, dimensions, doCreate, useBases);
+  }
+  
   /* HDF5 complex type */
   template<typename T>
   class hdf5_ComplexType
@@ -191,24 +235,11 @@ namespace CosmoTool {
   class InvalidDimensions: virtual std::exception {
   };
   
-  template<std::size_t r>
-  struct hdf5_extent_gen {
-    typedef typename boost::detail::multi_array::extent_gen<r> type;
-    
-    static inline type build(hsize_t *d)
-    {
-      return (hdf5_extent_gen<r-1>::build(d))[d[r-1]];
-    }
-  };
-
-  template<>
-  struct hdf5_extent_gen<0> {    
-    static inline boost::multi_array_types::extent_gen build(hsize_t *d)
-    {
-      return boost::extents;
-    }
-  };
   
+  // ----------------------------------------------------------------------
+  // Conditional resize support
+  // If the Array type support resize then it is called. Otherwise
+  // the dimensions are checked and lead to a failure if they are different
   
   template<typename Array> class array_has_resize {
     struct Fallback { int resize; };
@@ -258,10 +289,12 @@ namespace CosmoTool {
     hdf5_check_array(data, dims);
   }
 
+  // ----------------------------------------------------------------------
+  
   template<typename ArrayType, typename hdf5_data_type>
   void hdf5_read_array_typed(H5::CommonFG& fg, const std::string& data_set_name, 
-                       ArrayType& data, 
+                             ArrayType& data, 
-                       const hdf5_data_type& datatype, bool auto_resize = true)
+                             const hdf5_data_type& datatype, bool auto_resize = true)
   {
       H5::DataSet dataset = fg.openDataSet(data_set_name);
       H5::DataSpace dataspace = dataset.getSpace();
@@ -286,8 +319,7 @@ namespace CosmoTool {
   {
       typedef typename ArrayType::element T;
       
-      get_hdf5_data_type<T> hdf_data_type;
+      hdf5_read_array_typed(fg, data_set_name, data, get_hdf5_data_type<T>::type(), auto_resize);
-      hdf5_read_array_typed(fg, data_set_name, data, hdf_data_type.type(), auto_resize);
   }