Adding boost multi_array strip out

src/Grid/Encap.hpp

@@ -227,7 +227,7 @@ struct compare_cpu_encap_encap
  *
  */
 
-template<unsigned int p,typename Mem>
+/*template<unsigned int p,typename Mem>
 struct type_gpu_prop
 {
 	//! return a boost::fusion::vector<memory_c<....>....>
@@ -238,7 +238,7 @@ struct type_gpu_prop
 	typedef typename boost::remove_reference<rtype>::type mtype;
 	//! get the base type that the buffer is storing
 	typedef typename mtype::type type;
-};
+};*/
 
 /*! Stub encap
  *
@@ -506,7 +506,7 @@ public:
 	 * \return The reference of the data
 	 *
 	 */
-	template <unsigned int p> typename type_gpu_prop<p,typename memory_traits_inte<T>::type>::type::reference get()
+	template <unsigned int p> auto get() -> decltype(boost::fusion::at_c<p>(data).mem_r.operator[](k))
 	{
 		return boost::fusion::at_c<p>(data).mem_r.operator[](k);
 	}
@@ -518,7 +518,7 @@ public:
 	 * \return The reference of the data
 	 *
 	 */
-	template <unsigned int p> const typename type_gpu_prop<p,typename memory_traits_inte<T>::type>::type::reference get() const
+	template <unsigned int p> auto get() const -> decltype(boost::fusion::at_c<p>(data).mem_r.operator[](k))
 	{
 		return boost::fusion::at_c<p>(data).mem_r.operator[](k);
 	}

src/Grid/gpu_test/cuda_grid_unit_tests.cu

 #include "config.h"
 #include <Grid/map_grid.hpp>
 #include "Point_test.hpp"
+#include <stdio.h>
 
-__global__ void compute_stencil_grid(grid_gpu_ker<3,Point_test<float>> g1, grid_gpu_ker<3,Point_test<float>> g2, ite_gpu & ite_gpu)
+/*__global__ void grid_fill_vector(grid_gpu_ker<3,Point_test<float>> g1,  ite_gpu ite_gpu)
 {
+	GRID_ID_3(ite_gpu);
+
+	g1.template get<4>(key)[0] = 1.0;
+}*/
 
+__global__ void compute_stencil_grid(grid_gpu_ker<3,Point_test<float>> g1, grid_gpu_ker<3,Point_test<float>> g2, ite_gpu ite_gpu)
+{
 	GRID_ID_3(ite_gpu);
 
-/*	g2.template get<0>(key) = g1.template get<0>(key.move(0,1)) + g1.template get<0>(key.move(0,-1)) +
+	g2.template get<0>(key) = g1.template get<0>(key.move(0,1)) + g1.template get<0>(key.move(0,-1)) +
 			                  g1.template get<0>(key.move(1,1)) + g1.template get<0>(key.move(1,-1)) +
 							  g1.template get<0>(key.move(2,1)) + g1.template get<0>(key.move(2,-1)) -
-							  6.0*g1.template get<0>(key);*/
-
+							  6.0*g1.template get<0>(key);
 }
 
 __global__ void compute_stencil(float * prp_0, float * prp_1, int sz, grid_key_dx<3> start, grid_key_dx<3> stop)
@@ -89,7 +95,15 @@ void gpu_grid_3D_compute_grid_stencil(grid_gpu<3,Point_test<float>> & g1, grid_g
 	auto gpu_it = g2.getGPUIterator(start,stop);
 
 	auto g1k = g1.toGPU();
+	auto g2k = g2.toGPU();
 
-	compute_stencil_grid<<< gpu_it.thr, gpu_it.wthr >>>(g1k,g1k,gpu_it);
+	compute_stencil_grid<<< gpu_it.thr, gpu_it.wthr >>>(g1k,g2k,gpu_it);
 }
 
+/*void gpu_grid_fill_vector(grid_gpu<3,Point_test<float>> & g1, grid_key_dx<3> & start, grid_key_dx<3> & stop)
+{
+	auto gpu_it = g1.getGPUIterator(start,stop);
+
+	grid_fill_vector<<< gpu_it.thr, gpu_it.wthr >>>(g1.toGPU(),gpu_it);
+}*/
+

src/Grid/gpu_test/cuda_grid_unit_tests.cuh

@@ -16,4 +16,6 @@ void gpu_grid_3D_one(grid_gpu<3,Point_test<float>> & g);
 void gpu_grid_3D_compute_grid_stencil(grid_gpu<3,Point_test<float>> & g1, grid_gpu<3,Point_test<float>> & g2,
 		 	 	 	 	 	 	 	 grid_key_dx<3> & start, grid_key_dx<3> & stop);
 
+//void gpu_grid_fill_vector(grid_gpu<3,Point_test<float>> & g1, grid_key_dx<3> & start, grid_key_dx<3> & stop);
+
 #endif /* OPENFPM_DATA_SRC_GRID_CUDA_GPU_COMPUTE_CUH_ */

src/Grid/grid_key.hpp

@@ -39,7 +39,7 @@ public:
 	 * \param p1 initializer list
 	 *
 	 */
-	inline grid_key_dx(std::initializer_list<long int> p1)
+	__device__ __host__ inline grid_key_dx(std::initializer_list<long int> p1)
 	{
 		size_t i = 0;
 		for(long int x : p1)

src/Grid/grid_sm.hpp

@@ -439,7 +439,7 @@ public:
 	 */
 	template<typename a, typename ...lT> inline mem_id Lin(a v,lT...t) const
 	{
-#ifdef DEBUG
+#ifdef SE_CLASS1
 		if (sizeof...(t)+1 > N)
 		{
 			std::cerr << "Error incorrect grid cannot linearize more index than its dimensionality" << "\n";
@@ -532,7 +532,7 @@ public:
 	 * \return the size of the grid
 	 *
 	 */
-	inline size_t size() const
+	__device__ __host__ inline size_t size() const
 	{
 		return size_tot;
 	};
@@ -623,7 +623,7 @@ public:
 	 *
 	 */
 
-	inline size_t size(unsigned int i) const
+	__device__ __host__ inline size_t size(unsigned int i) const
 	{
 		return sz[i];
 	}

src/Grid/grid_unit_tests.hpp

@@ -878,7 +878,9 @@ BOOST_AUTO_TEST_CASE (gpu_computation_grid_stencil)
 	grid_gpu<3, Point_test<float> > c3(sz);
 	grid_gpu<3, Point_test<float> > c2(sz);
 	grid_key_dx<3> key1({1,1,1});
+	grid_key_dx<3> zero({0,0,0});
 	grid_key_dx<3> key2({62,62,62});
+	grid_key_dx<3> keyl({63,63,63});
 
 
 	c3.setMemory();
@@ -936,7 +938,6 @@ BOOST_AUTO_TEST_CASE (gpu_computation_grid_stencil)
 	while(it.isNext())
 	{
 		auto key = it.get();
-
 		good &= c2.get<0>(key) == 0;
 
 		++it;
@@ -944,6 +945,10 @@ BOOST_AUTO_TEST_CASE (gpu_computation_grid_stencil)
 
 	BOOST_REQUIRE_EQUAL(good,true);
 
+	// We also try to fill a vectorial quantity
+
+//	gpu_grid_fill_vector(c3,zero,keyl);
+
 	}
 
 	#endif

src/Grid/map_grid.hpp

@@ -219,6 +219,8 @@ struct copy_switch_memory_c
 	inline void operator()(T& t) const
 	{
 		boost::fusion::at_c<T::value>(dst).mem = boost::fusion::at_c<T::value>(src).mem;
+		// Increment the reference of mem
+		boost::fusion::at_c<T::value>(dst).mem->incRef();
 		boost::fusion::at_c<T::value>(dst).mem_r.bind_ref(boost::fusion::at_c<T::value>(src).mem_r);
 		boost::fusion::at_c<T::value>(dst).switchToDevicePtr();
 	}
@@ -245,6 +247,10 @@ struct grid_gpu_ker
 	grid_gpu_ker()
 	{}
 
+	grid_gpu_ker(const grid_sm<dim,void> & g1)
+	:g1(g1)
+	{}
+
 	grid_gpu_ker(const grid_gpu_ker & cpy)
 	:g1(cpy.g1)
 	{
@@ -253,6 +259,18 @@ struct grid_gpu_ker
 		boost::mpl::for_each_ref< boost::mpl::range_c<int,0,T::max_prop> >(bp_mc);
 	}
 
+	/*! \brief Return the internal grid information
+	 *
+	 * Return the internal grid information
+	 *
+	 * \return the internal grid
+	 *
+	 */
+	__device__ __host__ const grid_sm<dim,void> & getGrid() const
+	{
+		return g1;
+	}
+
 	/*! \brief Get the reference of the selected element
 	 *
 	 * \param v1 grid_key that identify the element in the grid
@@ -261,7 +279,7 @@ struct grid_gpu_ker
 	 *
 	 */
 	template <unsigned int p, typename r_type=decltype(mem_get<p,memory_traits_inte<T>,layout,grid_sm<dim,T>,grid_key_dx<dim>>::get(data_,g1,grid_key_dx<dim>()))>
-	__device__ inline r_type get(const grid_key_dx<dim> & v1)
+	__device__ __host__ inline r_type get(const grid_key_dx<dim> & v1)
 	{
 		return mem_get<p,memory_traits_inte<T>,decltype(this->data_),decltype(this->g1),decltype(v1)>::get(data_,g1,v1);
 	}
@@ -274,7 +292,7 @@ struct grid_gpu_ker
 	 *
 	 */
 	template <unsigned int p, typename r_type=decltype(mem_get<p,memory_traits_inte<T>,layout,grid_sm<dim,T>,grid_key_dx<dim>>::get(data_,g1,grid_key_dx<dim>()))>
-	__device__ inline const r_type get(const grid_key_dx<dim> & v1) const
+	__device__ __host__ inline const r_type get(const grid_key_dx<dim> & v1) const
 	{
 		return mem_get<p,memory_traits_inte<T>,decltype(this->data_),decltype(this->g1),decltype(v1)>::get(data_,g1,v1);
 	}
@@ -287,7 +305,7 @@ struct grid_gpu_ker
 	 *
 	 */
 	template <unsigned int p, typename r_type=decltype(mem_get<p,memory_traits_inte<T>,layout,grid_sm<dim,T>,grid_key_dx<dim>>::get_lin(data_,g1,0))>
-	__device__ inline r_type get(const size_t lin_id)
+	__device__ __host__ inline r_type get(const size_t lin_id)
 	{
 		return mem_get<p,memory_traits_inte<T>,decltype(this->data_),decltype(this->g1),grid_key_dx<dim>>::get_lin(data_,g1,lin_id);
 	}
@@ -300,7 +318,7 @@ struct grid_gpu_ker
 	 *
 	 */
 	template <unsigned int p, typename r_type=decltype(mem_get<p,memory_traits_inte<T>,layout,grid_sm<dim,T>,grid_key_dx<dim>>::get_lin(data_,g1,0))>
-	__device__ inline const r_type get(size_t lin_id) const
+	__device__ __host__ inline const r_type get(size_t lin_id) const
 	{
 		return mem_get<p,memory_traits_inte<T>,decltype(this->data_),decltype(this->g1),grid_key_dx<dim>>::get_lin(data_,g1,lin_id);
 	}
@@ -420,7 +438,7 @@ public:
 	 */
 	grid_gpu_ker<dim,T> toGPU()
 	{
-		grid_gpu_ker<dim,T> g;
+		grid_gpu_ker<dim,T> g(this->g1);
 		copy_switch_memory_c<T> cp_mc(this->data_,g.data_);
 
 		boost::mpl::for_each_ref< boost::mpl::range_c<int,0,T::max_prop> >(cp_mc);

src/memory_ly/memory_c.hpp

@@ -10,10 +10,11 @@
 #include <boost/fusion/include/mpl.hpp>
 #include <boost/mpl/vector.hpp>
 #include <array>
+
+#include "util/boost_multi_array_openfpm.hpp"
 #include "util/ct_array.hpp"
 #include "memory_array.hpp"
 #include "memory/memory.hpp"
-#include "util/boost_multi_array_openfpm.hpp"
 
 #ifndef MEMORY_C_HPP_
 #define MEMORY_C_HPP_

src/util/boost/boost_multi_array_base_openfpm.hpp

+/*
+ * boost_multi_array_base_openfpm.hpp
+ *
+ *  Created on: Jun 6, 2018
+ *      Author: i-bird
+ */
+
+#ifndef OPENFPM_DATA_SRC_UTIL_BOOST_BOOST_MULTI_ARRAY_BASE_OPENFPM_HPP_
+#define OPENFPM_DATA_SRC_UTIL_BOOST_BOOST_MULTI_ARRAY_BASE_OPENFPM_HPP_
+
+
+//
+// base.hpp - some implementation base classes for from which
+// functionality is acquired
+//
+
+#include "boost/multi_array/extent_range.hpp"
+#include "boost/multi_array/extent_gen.hpp"
+#include "boost/multi_array/index_range.hpp"
+#include "boost/multi_array/index_gen.hpp"
+#include "boost/multi_array/storage_order.hpp"
+#include "boost/multi_array/types.hpp"
+#include "boost/config.hpp"
+#include "boost/multi_array/concept_checks.hpp" //for ignore_unused_...
+#include "boost/mpl/eval_if.hpp"
+#include "boost/mpl/if.hpp"
+#include "boost/mpl/size_t.hpp"
+#include "boost/iterator/reverse_iterator.hpp"
+#include "boost/static_assert.hpp"
+#include "boost/type.hpp"
+#include "boost/assert.hpp"
+#include <cstddef>
+#include <memory>
+
+namespace boost {
+
+/////////////////////////////////////////////////////////////////////////
+// class declarations
+/////////////////////////////////////////////////////////////////////////
+
+template<typename T, std::size_t NumDims,
+  typename Allocator = std::allocator<T> >
+class multi_array_openfpm;
+
+// This is a public interface for use by end users!
+namespace multi_array_types {
+  typedef boost::detail::multi_array::size_type size_type;
+  typedef std::ptrdiff_t difference_type;
+  typedef boost::detail::multi_array::index index;
+  typedef detail::multi_array::index_range<index,size_type> index_range;
+  typedef detail::multi_array::extent_range<index,size_type> extent_range;
+  typedef detail::multi_array::index_gen<0,0> index_gen;
+  typedef detail::multi_array::extent_gen<0> extent_gen;
+}
+
+
+// boost::extents and boost::indices are now a part of the public
+// interface.  That way users don't necessarily have to create their
+// own objects.  On the other hand, one may not want the overhead of
+// object creation in small-memory environments.  Thus, the objects
+// can be left undefined by defining BOOST_MULTI_ARRAY_NO_GENERATORS
+// before loading multi_array.hpp.
+//#ifndef BOOST_MULTI_ARRAY_NO_GENERATORS
+//namespace {
+//  multi_array_types::extent_gen extents;
+//  multi_array_types::index_gen indices;
+//}
+//#endif // BOOST_MULTI_ARRAY_NO_GENERATORS
+
+namespace detail {
+namespace multi_array {
+
+template <typename T, std::size_t NumDims>
+class sub_array_openfpm;
+
+template <typename T, std::size_t NumDims, typename TPtr = const T*>
+class const_sub_array_openfpm;
+
+  template <typename T, typename TPtr, typename NumDims, typename Reference,
+            typename IteratorCategory>
+class array_iterator_openfpm;
+
+template <typename T, std::size_t NumDims, typename TPtr = const T*>
+class const_multi_array_view_openfpm;
+
+template <typename T, std::size_t NumDims>
+class multi_array_view_openfpm;
+
+/////////////////////////////////////////////////////////////////////////
+// class interfaces
+/////////////////////////////////////////////////////////////////////////
+
+class multi_array_base_openfpm {
+public:
+  typedef multi_array_types::size_type size_type;
+  typedef multi_array_types::difference_type difference_type;
+  typedef multi_array_types::index index;
+  typedef multi_array_types::index_range index_range;
+  typedef multi_array_types::extent_range extent_range;
+  typedef multi_array_types::index_gen index_gen;
+  typedef multi_array_types::extent_gen extent_gen;
+};
+
+//
+// value_accessor_n
+//  contains the routines for accessing elements from
+//  N-dimensional views.
+//
+template<typename T, std::size_t NumDims>
+class value_accessor_n_openfpm : public multi_array_base_openfpm {
+  typedef multi_array_base_openfpm super_type;
+public:
+  typedef typename super_type::index index;
+
+  //
+  // public typedefs used by classes that inherit from this base
+  //
+  typedef T element;
+  typedef boost::multi_array_openfpm<T,NumDims-1> value_type;
+  typedef sub_array_openfpm<T,NumDims-1> reference;
+  typedef const_sub_array_openfpm<T,NumDims-1> const_reference;
+
+protected:
+  // used by array operator[] and iterators to get reference types.
+  template <typename Reference, typename TPtr>
+  Reference access(boost::type<Reference>,index idx,TPtr base,
+                   const size_type* extents,
+                   const index* strides,
+                   const index* index_bases) const {
+
+    BOOST_ASSERT(idx - index_bases[0] >= 0);
+    BOOST_ASSERT(size_type(idx - index_bases[0]) < extents[0]);
+    // return a sub_array<T,NDims-1> proxy object
+    TPtr newbase = base + idx * strides[0];
+    return Reference(newbase,extents+1,strides+1,index_bases+1);
+
+  }
+
+  value_accessor_n_openfpm() { }
+  ~value_accessor_n_openfpm() { }
+};
+
+
+
+//
+// value_accessor_one
+//  contains the routines for accessing reference elements from
+//  1-dimensional views.
+//
+template<typename T>
+class value_accessor_one_openfpm : public multi_array_base_openfpm {
+  typedef multi_array_base_openfpm super_type;
+public:
+  typedef typename super_type::index index;
+  //
+  // public typedefs for use by classes that inherit it.
+  //
+  typedef T element;
+  typedef T value_type;
+  typedef T& reference;
+  typedef T const& const_reference;
+
+protected:
+  // used by array operator[] and iterators to get reference types.
+  template <typename Reference, typename TPtr>
+  Reference access(boost::type<Reference>,index idx,TPtr base,
+                   const size_type* extents,
+                   const index* strides,
+                   const index* index_bases) const {
+
+    ignore_unused_variable_warning(index_bases);
+    ignore_unused_variable_warning(extents);
+    BOOST_ASSERT(idx - index_bases[0] >= 0);
+    BOOST_ASSERT(size_type(idx - index_bases[0]) < extents[0]);
+    return *(base + idx * strides[0]);
+  }
+
+  value_accessor_one_openfpm() { }
+  ~value_accessor_one_openfpm() { }
+};
+
+
+/////////////////////////////////////////////////////////////////////////
+// choose value accessor begins
+//
+
+template <typename T, std::size_t NumDims>
+struct choose_value_accessor_n_openfpm {
+  typedef value_accessor_n_openfpm<T,NumDims> type;
+};
+
+template <typename T>
+struct choose_value_accessor_one_openfpm {
+  typedef value_accessor_one_openfpm<T> type;
+};
+
+template <typename T, typename NumDims>
+struct value_accessor_generator_openfpm {
+    BOOST_STATIC_CONSTANT(std::size_t, dimensionality = NumDims::value);
+
+  typedef typename
+  mpl::eval_if_c<(dimensionality == 1),
+                  choose_value_accessor_one_openfpm<T>,
+                  choose_value_accessor_n_openfpm<T,dimensionality>
+  >::type type;
+};
+
+
+template <class T, class NumDims>
+struct associated_types_openfpm
+  : value_accessor_generator_openfpm<T,NumDims>::type
+{};
+
+//
+// choose value accessor ends
+/////////////////////////////////////////////////////////////////////////
+
+// Due to some imprecision in the C++ Standard,
+// MSVC 2010 is broken in debug mode: it requires
+// that an Output Iterator have output_iterator_tag in its iterator_category if
+// that iterator is not bidirectional_iterator or random_access_iterator.
+#if BOOST_WORKAROUND(BOOST_MSVC, >= 1600)
+struct mutable_iterator_tag
+ : boost::random_access_traversal_tag, std::input_iterator_tag
+{
+  operator std::output_iterator_tag() const {
+    return std::output_iterator_tag();
+  }
+};
+#endif
+
+////////////////////////////////////////////////////////////////////////
+// multi_array_base
+////////////////////////////////////////////////////////////////////////
+template <typename T, std::size_t NumDims>
+class multi_array_impl_base_openfpm
+  :
+      public value_accessor_generator_openfpm<T,mpl::size_t<NumDims> >::type
+{
+  typedef associated_types_openfpm<T,mpl::size_t<NumDims> > types;
+public:
+  typedef typename types::index index;
+  typedef typename types::size_type size_type;
+  typedef typename types::element element;
+  typedef typename types::index_range index_range;
+  typedef typename types::value_type value_type;
+  typedef typename types::reference reference;
+  typedef typename types::const_reference const_reference;
+
+  template <std::size_t NDims>
+  struct subarray {
+    typedef boost::detail::multi_array::sub_array_openfpm<T,NDims> type;
+  };
+
+  template <std::size_t NDims>
+  struct const_subarray {
+    typedef boost::detail::multi_array::const_sub_array_openfpm<T,NDims> type;
+  };
+
+  template <std::size_t NDims>
+  struct array_view_openfpm {
+    typedef boost::detail::multi_array::multi_array_view_openfpm<T,NDims> type;
+  };
+
+  template <std::size_t NDims>
+  struct const_array_view_openfpm {
+  public:
+    typedef boost::detail::multi_array::const_multi_array_view_openfpm<T,NDims> type;
+  };
+
+  //
+  // iterator support
+  //
+#if BOOST_WORKAROUND(BOOST_MSVC, >= 1600)
+  // Deal with VC 2010 output_iterator_tag requirement
+  typedef array_iterator_openfpm<T,T*,mpl::size_t<NumDims>,reference,
+                         mutable_iterator_tag> iterator;
+#else
+  typedef array_iterator_openfpm<T,T*,mpl::size_t<NumDims>,reference,
+                         boost::random_access_traversal_tag> iterator;
+#endif
+  typedef array_iterator_openfpm<T,T const*,mpl::size_t<NumDims>,const_reference,
+                         boost::random_access_traversal_tag> const_iterator;
+
+  typedef ::boost::reverse_iterator<iterator> reverse_iterator;
+  typedef ::boost::reverse_iterator<const_iterator> const_reverse_iterator;
+
+  BOOST_STATIC_CONSTANT(std::size_t, dimensionality = NumDims);
+protected:
+
+  multi_array_impl_base_openfpm() { }
+  ~multi_array_impl_base_openfpm() { }
+
+  // Used by operator() in our array classes
+  template <typename Reference, typename IndexList, typename TPtr>
+  Reference access_element(boost::type<Reference>,
+                           const IndexList& indices,
+                           TPtr base,
+                           const size_type* extents,
+                           const index* strides,
+                           const index* index_bases) const {
+    boost::function_requires<
+      CollectionConcept<IndexList> >();
+    ignore_unused_variable_warning(index_bases);
+    ignore_unused_variable_warning(extents);
+#if !defined(NDEBUG) && !defined(BOOST_DISABLE_ASSERTS)
+    for (size_type i = 0; i != NumDims; ++i) {
+      BOOST_ASSERT(indices[i] - index_bases[i] >= 0);
+      BOOST_ASSERT(size_type(indices[i] - index_bases[i]) < extents[i]);
+    }
+#endif
+
+    index offset = 0;
+    {
+      typename IndexList::const_iterator i = indices.begin();
+      size_type n = 0;
+      while (n != NumDims) {
+        offset += (*i) * strides[n];
+        ++n;
+        ++i;
+      }
+    }
+    return base[offset];
+  }
+
+  template <typename StrideList, typename ExtentList>
+  void compute_strides(StrideList& stride_list, ExtentList& extent_list,
+                       const general_storage_order<NumDims>& storage)
+  {
+    // invariant: stride = the stride for dimension n
+    index stride = 1;
+    for (size_type n = 0; n != NumDims; ++n) {
+      index stride_sign = +1;
+
+      if (!storage.ascending(storage.ordering(n)))
+        stride_sign = -1;
+
+      // The stride for this dimension is the product of the
+      // lengths of the ranks minor to it.
+      stride_list[storage.ordering(n)] = stride * stride_sign;
+
+      stride *= extent_list[storage.ordering(n)];
+    }
+  }
+
+  // This calculates the offset to the array base pointer due to:
+  // 1. dimensions stored in descending order
+  // 2. non-zero dimension index bases