grid_base_implementation.hpp

/*
 * grid_base_implementation.hpp
 *
 *  Created on: May 2, 2016
 *      Author: i-bird
 */
 
#ifndef OPENFPM_DATA_SRC_GRID_GRID_BASE_IMPLEMENTATION_HPP_
#define OPENFPM_DATA_SRC_GRID_GRID_BASE_IMPLEMENTATION_HPP_
 
#include "grid_base_impl_layout.hpp"
#include "util/cuda_util.hpp"
#include "cuda/cuda_grid_gpu_funcs.cuh"
#include "util/create_vmpl_sequence.hpp"
#include "util/object_si_di.hpp"
 
constexpr int DATA_ON_HOST = 32;
constexpr int DATA_ON_DEVICE = 64;
constexpr int EXACT_RESIZE = 128;
 
template<bool np,typename T>
struct skip_init
{
    static bool skip_()
    {
        return true;
    }
};
 
template<typename T>
struct skip_init<true,T>
{
    static bool skip_()
    {
        return T::noPointers();
    }
};
 
#ifdef __NVCC__
 
template<bool active>
struct copy_ndim_grid_device_active_impl
    {
    template<typename grid_type1, typename grid_type2, typename ite_gpu_type>
    static inline void copy(grid_type1 & g1, grid_type2 & g2, ite_gpu_type & ite)
    {
 
    }
 
    template<typename grid_type1, typename grid_type2, typename ite_gpu_type>
    static inline void copy_block(grid_type1 & g1, grid_type2 & g2, ite_gpu_type & ite)
    {
    }
};
 
template<>
struct copy_ndim_grid_device_active_impl<true>
{
    template<typename grid_type1, typename grid_type2, typename ite_gpu_type>
    static inline void copy(grid_type1 & g1, grid_type2 & g2, ite_gpu_type & ite)
    {
        CUDA_LAUNCH((copy_ndim_grid_device<grid_type1::dims,decltype(g1.toKernel())>),ite,g2.toKernel(),g1.toKernel());
    }
 
    template<typename grid_type1, typename grid_type2, typename ite_gpu_type>
    static inline void copy_block(grid_type1 & g1, grid_type2 & g2, ite_gpu_type & ite)
    {
        CUDA_LAUNCH((copy_ndim_grid_block_device<grid_type1::dims,decltype(g1.toKernel())>),ite,g2.toKernel(),g1.toKernel());
    }
};
 
template<typename S,typename grid_dst_type, typename grid_src_type>
void copy_grid_to_grid(grid_dst_type & gdst, const grid_src_type & gsrc, 
                       grid_key_dx<grid_dst_type::dims> & start, grid_key_dx<grid_dst_type::dims> & stop,
                       int blockSize)
{
    if (grid_dst_type::dims <= 3)
    {
        auto ite = gsrc.getGPUIterator(start,stop);
        bool has_work = has_work_gpu(ite);
 
        if (has_work == true)
        {
            if (blockSize == 1)
            {
                copy_ndim_grid_device_active_impl<S::isDeviceHostSame() == false>::copy(gdst,gsrc,ite);
            }
            else
            {
                move_work_to_blocks(ite);
 
                ite.thr.x = blockSize;
 
                copy_ndim_grid_device_active_impl<S::isDeviceHostSame() == false>::copy_block(gdst,gsrc,ite);
            }
        }
    }
    else
    {
        grid_key_dx<1> start;
        start.set_d(0,0);
        grid_key_dx<1> stop({});
        stop.set_d(0,gsrc.getGrid().size());
 
        size_t sz[1];
        sz[0]= gsrc.getGrid().size();
 
        grid_sm<1,void> g_sm_copy(sz);
 
        auto ite = getGPUIterator_impl<1>(g_sm_copy,start,stop);
 
        copy_ndim_grid_device_active_impl<S::isDeviceHostSame() == false>::copy(gdst,gsrc,ite);
    }
}
 
#endif
 
template<typename dest_type, typename src_type, unsigned int ... prp>
void copy_with_openmp_prp(const dest_type & dst, const src_type & src, ite_gpu<dest_type::dims> ite)
{
    #ifdef CUDIFY_USE_OPENMP
    auto lamb = [&dst,&src,&ite] __device__ (dim3 & blockIdx, dim3 & threadIdx) 
    {
        grid_key_dx<dest_type::dims> i;
 
        if (dest_type::dims == 1)
        {
            i.set_d(0,blockIdx.x*blockDim.x + threadIdx.x + ite.start.get(0));
            if (i.get(0) >= src.size(0))    {return;}
        }
        else if (dest_type::dims == 2)
        {
            i.set_d(0,blockIdx.x*blockDim.x + threadIdx.x + ite.start.get(0));
            i.set_d(1,blockIdx.y*blockDim.y + threadIdx.y + ite.start.get(1));
            if (i.get(0) >= src.size(0) || i.get(1) >= src.size(1)) {return;}
        }
        else if (dest_type::dims == 3)
        {
            i.set_d(0,blockIdx.x*blockDim.x + threadIdx.x + ite.start.get(0));
            i.set_d(1,blockIdx.y*blockDim.y + threadIdx.y + ite.start.get(1));
            i.set_d(2,blockIdx.z*blockDim.z + threadIdx.z + ite.start.get(2));
            if (i.get(0) >= src.size(0) || i.get(1) >= src.size(1) || i.get(2) >= src.size(2))  {return;}
        }
 
        object_si_di<decltype(src.get_o(i)),decltype(dst.get_o(i)),OBJ_ENCAP,prp ...>(src.get_o(i),dst.get_o(i));
    };
 
    CUDA_LAUNCH_LAMBDA(ite,lamb);
    #else
    std::cout << __FILE__ << ":" << __LINE__ << " error CUDA on back end is disabled" << std::endl;
    #endif
}
 
 
#ifdef CUDA_GPU
 
#define GRID_ID_3_RAW(start,stop) int x[3] = {threadIdx.x + blockIdx.x * blockDim.x + start.get(0),\
                                         threadIdx.y + blockIdx.y * blockDim.y + start.get(1),\
                                         threadIdx.z + blockIdx.z * blockDim.z + start.get(2)};\
                                         \
                                         if (x[0] > stop.get(0) || x[1] > stop.get(1) || x[2] > stop.get(2))\
                                         {return;}
 
#define GRID_ID_3_TRAW(start,stop) int tx = threadIdx.x + blockIdx.x * blockDim.x + start.get(0);\
                                   int ty = threadIdx.y + blockIdx.y * blockDim.y + start.get(1);\
                                   int tz = threadIdx.z + blockIdx.z * blockDim.z + start.get(2);\
                                         \
                                         if (tx > stop.get(0) || ty > stop.get(1) || tz > stop.get(2))\
                                         {return;}
 
#define GRID_ID_3(ite_gpu) grid_key_dx<3,int> key;\
                              key.set_d(0,threadIdx.x + blockIdx.x * blockDim.x + ite_gpu.start.get(0));\
                              key.set_d(1,threadIdx.y + blockIdx.y * blockDim.y + ite_gpu.start.get(1));\
                              key.set_d(2,threadIdx.z + blockIdx.z * blockDim.z + ite_gpu.start.get(2));\
                                         \
                                         if (key.get(0) > ite_gpu.stop.get(0) || key.get(1) > ite_gpu.stop.get(1) || key.get(2) > ite_gpu.stop.get(2))\
                                         {return;}
 
#define GRID_ID_2(ite_gpu) grid_key_dx<2,int> key;\
                              key.set_d(0,threadIdx.x + blockIdx.x * blockDim.x + ite_gpu.start.get(0));\
                              key.set_d(1,threadIdx.y + blockIdx.y * blockDim.y + ite_gpu.start.get(1));\
                                         \
                                         if (key.get(0) > ite_gpu.stop.get(0) || key.get(1) > ite_gpu.stop.get(1))\
                                         {return;}
 
#ifdef __NVCC__
 
 
template<unsigned int dim, typename ids_type = int>
struct grid_p
{
    __device__ static inline grid_key_dx<dim,ids_type> get_grid_point(const grid_sm<dim,void> & g)
    {
        grid_key_dx<dim,ids_type> key;
 
        key.set_d(0,blockIdx.x * blockDim.x + threadIdx.x);
        key.set_d(1,blockIdx.y * blockDim.y + threadIdx.y);
 
        unsigned int bz = blockIdx.z * blockDim.z + threadIdx.z;
        key.set_d(2,bz % g.size(2));
 
        for (unsigned int i = 3 ; i < dim ; i++)
        {
            bz /= g.size(i);
            key.set_d(i,bz % g.size(i));
        }
 
        return key;
    }
 
    __device__ static inline grid_key_dx<dim,ids_type> get_grid_point(const openfpm::array<ids_type,dim,unsigned int> & g)
    {
        grid_key_dx<dim,ids_type> key;
 
        key.set_d(0,blockIdx.x * blockDim.x + threadIdx.x);
        key.set_d(1,blockIdx.y * blockDim.y + threadIdx.y);
 
        unsigned int bz = blockIdx.z * blockDim.z + threadIdx.z;
        key.set_d(2,bz % g[2]);
 
        for (unsigned int i = 3 ; i < dim ; i++)
        {
            bz /= g[i];
            key.set_d(i,bz % g[i]);
        }
 
        return key;
    }
};
 
template<typename ids_type>
struct grid_p<3,ids_type>
{
    __device__ static inline grid_key_dx<3,ids_type> get_grid_point(const grid_sm<3,void> & g)
    {
        grid_key_dx<3,unsigned int> key;
 
        key.set_d(0,blockIdx.x * blockDim.x + threadIdx.x);
        key.set_d(1,blockIdx.y * blockDim.y + threadIdx.y);
        key.set_d(2,blockIdx.z * blockDim.z + threadIdx.z);
 
        return key;
    }
 
    __device__ static inline grid_key_dx<3,ids_type> get_grid_point(const openfpm::array<ids_type,3,unsigned int> & g)
    {
        grid_key_dx<3,ids_type> key;
 
        key.set_d(0,blockIdx.x * blockDim.x + threadIdx.x);
        key.set_d(1,blockIdx.y * blockDim.y + threadIdx.y);
        key.set_d(2,blockIdx.z * blockDim.z + threadIdx.z);
 
        return key;
    }
};
 
template<typename ids_type>
struct grid_p<2,ids_type>
{
    __device__ static inline grid_key_dx<2,ids_type> get_grid_point(const grid_sm<2,void> & g)
    {
        grid_key_dx<2,ids_type> key;
 
        key.set_d(0,blockIdx.x * blockDim.x + threadIdx.x);
        key.set_d(1,blockIdx.y * blockDim.y + threadIdx.y);
 
        return key;
    }
 
    __device__ static inline grid_key_dx<2,ids_type> get_grid_point(const openfpm::array<ids_type,2,unsigned int> & g)
    {
        grid_key_dx<2,ids_type> key;
 
        key.set_d(0,blockIdx.x * blockDim.x + threadIdx.x);
        key.set_d(1,blockIdx.y * blockDim.y + threadIdx.y);
 
        return key;
    }
};
 
template<typename ids_type>
struct grid_p<1,ids_type>
{
    __device__ static inline grid_key_dx<1,unsigned int> get_grid_point(const grid_sm<1,void> & g)
    {
        grid_key_dx<1,unsigned int> key;
 
        key.set_d(0,blockIdx.x * blockDim.x + threadIdx.x);
 
        return key;
    }
};
 
#endif
 
 
template<unsigned int dim>
void move_work_to_blocks(ite_gpu<dim> & ite)
{
    if (dim == 1)
    {
        ite.wthr.x = ite.wthr.x * ite.thr.x;
        ite.thr.x = 1;
    }
    else if(dim == 2)
    {
        ite.wthr.x = ite.wthr.x * ite.thr.x;
        ite.wthr.y = ite.wthr.y * ite.thr.y;
        ite.thr.x = 1;
        ite.thr.y = 1;
 
    }
    else
    {
        ite.wthr.x = ite.wthr.x * ite.thr.x;
        ite.wthr.x = ite.wthr.y * ite.thr.y;
        ite.wthr.x = ite.wthr.z * ite.thr.z;
        ite.thr.x = 1;
        ite.thr.y = 1;
        ite.thr.z = 1;
    }
}
 
#endif
 
#include "copy_grid_fast.hpp"
 
template<typename T>
struct copy_grid_fast_caller
{
    template<typename grid_type>
    static void call(grid_type & gd, const grid_type & gs, const Box<grid_type::dims,size_t> & box_src, const Box<grid_type::dims,size_t> & box_dst)
    {
        std::cout << "Error: " << __FILE__ << ":" << __LINE__ << " copy_grid_fast_caller failure" << std::endl;
    }
};
 
template<int ... prp>
struct copy_grid_fast_caller<index_tuple_sq<prp ...>>
{
    template<typename grid_type>
    static void call(grid_type & gd, const grid_type & gs, const Box<grid_type::dims,size_t> & box_src, const Box<grid_type::dims,size_t> & box_dst)
    {
        grid_key_dx<grid_type::dims> cnt[1];
        cnt[0].zero();
 
        typedef typename std::remove_reference<decltype(gd)>::type grid_cp;
        typedef typename std::remove_reference<decltype(gd.getGrid())>::type grid_info_cp;
 
        typedef typename to_int_sequence<0,grid_type::value_type::max_prop>::type result;
 
        copy_grid_fast<!is_contiguos<prp...>::type::value || has_pack_gen<typename grid_type::value_type>::value,
                                   grid_type::dims,
                                   grid_cp,
                                   grid_info_cp>::copy(gs.getGrid(),
                                   gd.getGrid(),
                                   box_src,
                                   box_dst,
                                   gs,gd,cnt);
    }
};
 
template<unsigned int dim,
         typename T,
         typename S,
         template<typename> class layout_base,
         typename ord_type = grid_sm<dim,void> >
class grid_base_impl
{
    typedef typename layout_base<T>::type layout;
 
    typedef typename apply_transform<layout_base,T>::type T_;
 
public:
 
    typedef layout layout_type;
 
    static constexpr unsigned int dims = dim;
 
    typedef grid_key_dx<dim> access_key;
 
    typedef typename T::type T_type;
 
    typedef layout_base<T> layout_base_;
 
    typedef ord_type linearizer_type;
 
    typedef T background_type;
 
protected:
 
    layout data_;
 
    ord_type g1;
 
private:
 
    bool is_mem_init = false;
 
    bool isExternal;
 
#ifdef SE_CLASS1
 
    inline void check_init() const
    {
#ifndef __CUDA_ARCH__
        if (is_mem_init == false)
        {
            std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " you must call SetMemory before access the grid\n";
            ACTION_ON_ERROR(GRID_ERROR_OBJECT);
        }
#endif
    }
 
    inline void check_bound(const grid_key_dx<dim> & v1) const
    {
#ifndef __CUDA_ARCH__
        for (long int i = 0 ; i < dim ; i++)
        {
            if (v1.get(i) >= (long int)getGrid().size(i))
            {
                std::cerr << "Error " __FILE__ << ":" << __LINE__ <<" grid overflow " << "x=[" << i << "]=" << v1.get(i) << " >= " << getGrid().size(i) << "\n";
                ACTION_ON_ERROR(GRID_ERROR_OBJECT);
            }
            else if (v1.get(i) < 0)
            {
                std::cerr << "Error " __FILE__ << ":" << __LINE__ <<" grid overflow " << "x=[" << i << "]=" << v1.get(i) << " is negative " << "\n";
                ACTION_ON_ERROR(GRID_ERROR_OBJECT);
            }
        }
#endif
    }
 
    inline void check_bound(size_t v1) const
    {
#ifndef __CUDA_ARCH__
        if (v1 >= getGrid().size())
        {
            std::cerr << "Error " __FILE__ << ":" << __LINE__ <<" grid overflow " << v1<< " >= " << getGrid().size() << "\n";
            ACTION_ON_ERROR(GRID_ERROR_OBJECT);
        }
#endif
    }
 
    template<typename Mem> inline void check_bound(const grid_base_impl<dim,T,Mem,layout_base, ord_type> & g,const grid_key_dx<dim> & key2) const
    {
#ifndef __CUDA_ARCH__
        for (size_t i = 0 ; i < dim ; i++)
        {
            if (key2.get(i) >= (long int)g.getGrid().size(i))
            {
                std::cerr << "Error " __FILE__ << ":" << __LINE__ <<" grid overflow " << "x=[" << i << "]=" << key2.get(i) << " >= " << g.getGrid().size(i) << "\n";
                ACTION_ON_ERROR(GRID_ERROR_OBJECT);
            }
            else if (key2.get(i) < 0)
            {
                std::cerr << "Error " __FILE__ << ":" << __LINE__ <<" grid overflow " << "x=[" << i << "]=" << key2.get(i) << " is negative " << "\n";
                ACTION_ON_ERROR(GRID_ERROR_OBJECT);
            }
        }
#endif
    }
 
    template<typename Mem, template <typename> class layout_base2> 
    inline void check_bound(const grid_base_impl<dim,T,Mem,layout_base2,ord_type> & g,const grid_key_dx<dim> & key2) const
    {
#ifndef __CUDA_ARCH__
        for (size_t i = 0 ; i < dim ; i++)
        {
            if (key2.get(i) >= (long int)g.getGrid().size(i))
            {
                std::cerr << "Error " __FILE__ << ":" << __LINE__ <<" grid overflow " << "x=[" << i << "]=" << key2.get(i) << " >= " << g.getGrid().size(i) << "\n";
                ACTION_ON_ERROR(GRID_ERROR_OBJECT);
            }
            else if (key2.get(i) < 0)
            {
                std::cerr << "Error " __FILE__ << ":" << __LINE__ <<" grid overflow " << "x=[" << i << "]=" << key2.get(i) << " is negative " << "\n";
                ACTION_ON_ERROR(GRID_ERROR_OBJECT);
            }
        }
#endif
    }
 
    template<typename Mem> inline void check_bound(const grid_base_impl<dim,T,Mem,layout_base> & g,const size_t & key2) const
    {
#ifndef __CUDA_ARCH__
        if (key2 >= g.getGrid().size())
        {
            std::cerr << "Error " __FILE__ << ":" << __LINE__ <<" grid overflow " << key2 << " >= " << getGrid().size() << "\n";
            ACTION_ON_ERROR(GRID_ERROR_OBJECT);
        }
#endif
    }
 
#endif
 
    void resize_impl_device(const size_t (& sz)[dim],grid_base_impl<dim,T,S,layout_base,ord_type> & grid_new, unsigned int blockSize = 1)
    {
#if defined(CUDA_GPU) && defined(__NVCC__)
 
            // Compile time-cheking that make sense to call a GPU kernel to copy.
            
 
            grid_key_dx<dim> start;
            grid_key_dx<dim> stop;
 
            for (size_t i = 0 ; i < dim ; i++)
            {
                start.set_d(i,0);
 
                // take the smallest
                if (grid_new.g1.size(i) < sz[i])
                {stop.set_d(i,grid_new.g1.size(i)-1);}
                else
                {stop.set_d(i,sz[i]-1);}
            }
 
            copy_grid_to_grid<S>(grid_new,*this,start,stop,blockSize);
 
#else
 
            std::cout << __FILE__ << ":" << __LINE__ << " error: the function resize require the launch of a kernel, but it seem that this" <<
                                                        " file (grid_base_implementation.hpp) has not been compiled with NVCC  " << std::endl;
 
#endif
    }
 
    void resize_impl_host(const size_t (& sz)[dim], grid_base_impl<dim,T,S,layout_base,ord_type> & grid_new)
    {
        size_t sz_c[dim];
        for (size_t i = 0 ; i < dim ; i++)
        {sz_c[i] = (g1.size(i) < sz[i])?g1.size(i):sz[i];}
 
        grid_sm<dim,void> g1_c(sz_c);
 
        grid_key_dx_iterator<dim> it(g1_c);
 
        while(it.isNext())
        {
            // get the grid key
            grid_key_dx<dim> key = it.get();
 
            // create a copy element
 
            grid_new.get_o(key) = this->get_o(key);
 
            ++it;
        }
    }
 
    void resize_impl_memset(grid_base_impl<dim,T,S,layout_base,ord_type> & grid_new)
    {
        if (isExternal == true)
        {
            mem_setext<typename std::remove_reference<decltype(grid_new)>::type,S,layout_base<T>,decltype(data_)>::set(grid_new,*this,this->data_);
        }
        else
            grid_new.setMemory();
 
#if defined(CUDIFY_USE_SEQUENTIAL) || defined(CUDIFY_USE_OPENMP)
 
        base_gpu = grid_toKernelImpl<is_layout_inte<layout_base<T_>>::value,grid_gpu_ker<dim,T_,layout_base,linearizer_type>,dim,T_,S>::toKernel(*this);
 
#endif
    }
 
public:
 
    // Implementation of packer and unpacker for grid
    #include "grid_pack_unpack.ipp"
 
    typedef int yes_i_am_grid;
 
    typedef typename memory_traits_lin<typename T::type>::type memory_lin;
 
    // you can access all the properties of T
    typedef encapc<dim,T,layout> container;
 
    typedef T value_type;
 
    grid_base_impl() THROW
    :g1(0),isExternal(false)
    {
        // Add this pointer
    }
 
    grid_base_impl(const grid_base_impl & g) THROW
    :isExternal(false)
    {
        this->operator=(g);
    }
 
    grid_base_impl(const size_t & sz) THROW
    :g1(sz),isExternal(false)
    {
        // Add this pointer
    }
 
    grid_base_impl(const size_t (& sz)[dim]) THROW
    :g1(sz),isExternal(false)
    {
        // Add this pointer
    }
 
    ~grid_base_impl() THROW
    {
        // delete this pointer
    }
 
    grid_base_impl<dim,T,S,layout_base> & operator=(const grid_base_impl<dim,T,S,layout_base> & g)
    {
        swap(g.duplicate());
 
#if defined(CUDIFY_USE_SEQUENTIAL) || defined(CUDIFY_USE_OPENMP)
 
        base_gpu = grid_toKernelImpl<is_layout_inte<layout_base<T_>>::value,grid_gpu_ker<dim,T_,layout_base,linearizer_type>,dim,T_,S>::toKernel(*this);
        g.base_gpu = grid_toKernelImpl<is_layout_inte<layout_base<T_>>::value,grid_gpu_ker<dim,T_,layout_base,linearizer_type>,dim,T_,S>::toKernel(g);
 
#endif
 
        return *this;
    }
 
    grid_base_impl<dim,T,S,layout_base> & operator=(grid_base_impl<dim,T,S,layout_base> && g)
    {
        swap(g);
 
#if defined(CUDIFY_USE_SEQUENTIAL) || defined(CUDIFY_USE_OPENMP)
 
        base_gpu = grid_toKernelImpl<is_layout_inte<layout_base<T_>>::value,grid_gpu_ker<dim,T,layout_base,linearizer_type>,dim,T_,S>::toKernel(*this);
 
#endif
 
        return *this;
    }
 
    bool operator==(const grid_base_impl<dim,T,S,layout_base> & g)
    {
        // check if the have the same size
        if (g1 != g.g1)
            return false;
 
        auto it = getIterator();
 
        while (it.isNext())
        {
            auto key = it.get();
 
            if (this->get_o(key) != g.get_o(key))
                return false;
 
            ++it;
        }
 
        return true;
    }
 
    grid_base_impl<dim,T,S,layout_base> duplicate() const THROW
    {
 
        grid_base_impl<dim,T,S,layout_base> grid_new(g1.getSize());
 
        grid_new.setMemory();
 
        // We know that, if it is 1D we can safely copy the memory
//      if (dim == 1)
//      {
//          grid_new.data_.mem->copy(*data_.mem);
//      }
//      else
//      {
 
            grid_key_dx_iterator<dim> it(g1);
 
            while(it.isNext())
            {
                grid_new.set(it.get(),*this,it.get());
 
                ++it;
            }
//      }
 
        // copy grid_new to the base
 
        return grid_new;
    }
 
#ifdef CUDA_GPU
 
    struct ite_gpu<dim> getGPUIterator(const grid_key_dx<dim,long int> & key1, const grid_key_dx<dim,long int> & key2, size_t n_thr = default_kernel_wg_threads_) const
    {
        return getGPUIterator_impl<dim>(g1,key1,key2,n_thr);
    }
 
#endif
 
    int size(int i) const
    {
        return g1.size(i);
    }
 
    const ord_type & getGrid() const
    {
        return g1;
    }
 
    void setMemory()
    {
        mem_setm<S,layout_base<T>,decltype(this->data_),decltype(this->g1)>::setMemory(data_,g1,is_mem_init);
 
#if defined(CUDIFY_USE_SEQUENTIAL) || defined(CUDIFY_USE_OPENMP)
 
        base_gpu = grid_toKernelImpl<is_layout_inte<layout_base<T_>>::value,grid_gpu_ker<dim,T_,layout_base,linearizer_type>,dim,T_,S>::toKernel(*this);
 
#endif
 
    }
 
    template<unsigned int p>
    auto getMemory() -> decltype(boost::fusion::at_c<p>(data_).getMemory())
    {
        return boost::fusion::at_c<p>(data_).getMemory();
    }
 
    template<unsigned int p = 0> void setMemory(S & m)
    {
        isExternal = true;
 
        bool skip_ini = skip_init<has_noPointers<T>::value,T>::skip_();
 
        mem_setmemory<decltype(data_),S,layout_base<T>>::template setMemory<p>(data_,m,g1.size(),skip_ini);
 
        is_mem_init = true;
 
#if defined(CUDIFY_USE_SEQUENTIAL) || defined(CUDIFY_USE_OPENMP)
 
        base_gpu = grid_toKernelImpl<is_layout_inte<layout_base<T_>>::value,grid_gpu_ker<dim,T_,layout_base,linearizer_type>,dim,T_,S>::toKernel(*this);
 
#endif
    }
 
    void setMemoryArray(S * m)
    {
        isExternal = true;
 
        bool skip_ini = skip_init<has_noPointers<T>::value,T>::skip_();
 
        mem_setmemory<decltype(data_),S,layout_base<T>>::setMemoryArray(*this,m,g1.size(),skip_ini);
 
        is_mem_init = true;
 
#if defined(CUDIFY_USE_SEQUENTIAL) || defined(CUDIFY_USE_OPENMP)
 
        base_gpu = grid_toKernelImpl<is_layout_inte<layout_base<T_>>::value,grid_gpu_ker<dim,T_,layout_base,linearizer_type>,dim,T_,S>::toKernel(*this);
 
#endif
    }
 
    template<unsigned int p = 0> void * getPointer()
    {
        return mem_getpointer<decltype(data_),layout_base_>::template getPointer<p>(data_);
    }
 
    template<unsigned int p = 0> const void * getPointer() const
    {
        return mem_getpointer<decltype(data_),layout_base_>::template getPointer<p>(data_);
    }
 
 
    template <unsigned int p, typename r_type=decltype(layout_base<T>::template get<p>(data_,g1,grid_key_dx<dim>()))>
    inline r_type insert(const grid_key_dx<dim> & v1)
    {
#ifdef SE_CLASS1
        check_init();
        check_bound(v1);
#endif
        return this->get<p>(v1);
    }
 
 
    template <unsigned int p, typename r_type=decltype(layout_base<T>::template get<p>(data_,g1,grid_key_dx<dim>()))>
    __device__ __host__ inline r_type get_usafe(const grid_key_dx<dim> & v1)
    {
#ifdef SE_CLASS1
        check_init();
#endif
        return layout_base<T>::template get<p>(data_,g1,v1);
    }
 
    int getBlockEdgeSize()
    {
        return 1;
    }
 
    void setGPUInsertBuffer(unsigned int nb, unsigned int nt)
    {}
 
    template <unsigned int p, typename r_type=decltype(layout_base<T>::template get_c<p>(data_,g1,grid_key_dx<dim>()))>
    __device__ __host__ inline r_type get_unsafe(const grid_key_dx<dim> & v1) const
    {
#ifdef SE_CLASS1
        check_init();
#endif
        return layout_base<T>::template get_c<p>(data_,g1,v1);
    }
 
    template <unsigned int p, typename r_type=decltype(layout_base<T>::template get<p>(data_,g1,grid_key_dx<dim>()))>
    __device__ __host__ inline r_type get(const grid_key_dx<dim> & v1)
    {
#ifdef SE_CLASS1
        check_init();
        check_bound(v1);
#endif
        return layout_base<T>::template get<p>(data_,g1,v1);
    }
 
    __device__ __host__ inline unsigned char getFlag(const grid_key_dx<dim> & v1) const
    {
        return 0;
    }
 
    template <unsigned int p, typename r_type=decltype(layout_base<T>::template get_c<p>(data_,g1,grid_key_dx<dim>()))>
    __device__ __host__ inline r_type get(const grid_key_dx<dim> & v1) const
    {
#ifdef SE_CLASS1
        check_init();
        check_bound(v1);
#endif
        return layout_base<T>::template get_c<p>(data_,g1,v1);
    }
 
    template <unsigned int p, typename r_type=decltype(layout_base<T>::template get_lin<p>(data_,g1,0))>
    __device__ __host__ inline r_type get(const size_t lin_id)
    {
#ifdef SE_CLASS1
        check_init();
        check_bound(lin_id);
#endif
        return layout_base<T>::template get_lin<p>(data_,g1,lin_id);
    }
 
    template <unsigned int p, typename r_type=decltype(layout_base<T>::template get_lin<p>(data_,g1,0))>
    __device__ __host__ inline const r_type get(size_t lin_id) const
    {
#ifdef SE_CLASS1
        check_init();
        check_bound(lin_id);
#endif
        return layout_base<T>::get_lin_const(data_,g1,lin_id);
    }
 
 
    inline encapc<dim,T,layout> get_o(const grid_key_dx<dim> & v1)
    {
#ifdef SE_CLASS1
        check_init();
        check_bound(v1);
#endif
        return mem_geto<dim,T,layout_base<T>,decltype(this->data_),decltype(this->g1),decltype(v1)>::get(data_,g1,v1);
    }
 
    inline const encapc<dim,T,layout> get_o(const grid_key_dx<dim> & v1) const
    {
#ifdef SE_CLASS1
        check_init();
        check_bound(v1);
#endif
        return mem_geto<dim,T,layout_base<T>,decltype(this->data_),decltype(this->g1),decltype(v1)>
               ::get(const_cast<typename std::add_lvalue_reference<decltype(this->data_)>::type>(data_),
                     g1,v1);
    }
 
 
    inline encapc<dim,T,layout> insert_o(const grid_key_dx<dim> & v1)
    {
#ifdef SE_CLASS1
        check_init();
        check_bound(v1);
#endif
        return mem_geto<dim,T,layout_base<T>,decltype(this->data_),decltype(this->g1),decltype(v1)>::get(data_,g1,v1);
    }
 
    inline encapc<dim,T,layout> get_o(size_t v1)
    {
#ifdef SE_CLASS1
        check_init();
        check_bound(v1);
#endif
        return mem_geto<dim,T,layout_base<T>,decltype(this->data_),decltype(this->g1),decltype(v1)>::get_lin(data_,v1);
    }
 
    inline const encapc<dim,T,layout> get_o(size_t v1) const
    {
#ifdef SE_CLASS1
        check_init();
        check_bound(v1);
#endif
        return mem_geto<dim,T,layout_base<T>,decltype(this->data_),decltype(this->g1),decltype(v1)>
               ::get_lin(const_cast<typename std::add_lvalue_reference<decltype(this->data_)>::type>(data_),v1);
    }
 
    template<int prp>
    void fill(unsigned char fl)
    {
        if (prp != 0 || is_layout_mlin<layout_base<T>>::type::value == false)
        {
            std::cout << "Error: " << __FILE__ << ":" << __LINE__ << " unsupported fill operation " << std::endl;
        }
 
        memset(getPointer(),fl,size() * sizeof(T));
    }
 
    void remove(Box<dim,long int> & section_to_delete)
    {}
 
    void copyRemoveReset()
    {
    }
 
    bool isSkipLabellingPossible()
    {
        return false;
    }
 
    void copy_to(const grid_base_impl<dim,T,S,layout_base> & grid_src,
                 const Box<dim,long int> & box_src,
                 const Box<dim,long int> & box_dst)
    {
        // fix box_dst
 
         Box<dim,size_t> box_src_;
         Box<dim,size_t> box_dst_;
 
        for (size_t i = 0 ; i < dim ; i++)
        {
            if (box_dst.getHigh(i) >= (long int)g1.size(i))
            {
                long int shift = box_dst.getHigh(i) - g1.size(i) + 1;
                box_dst_.setHigh(i,box_dst.getHigh(i) - shift);
                box_src_.setHigh(i,box_src.getHigh(i) - shift);
            }
            else
            {
                box_dst_.setHigh(i,box_dst.getHigh(i));
                box_src_.setHigh(i,box_src.getHigh(i));
            }
 
            if (box_dst.getLow(i) < 0)
            {
                long int shift = -box_dst.getLow(i);
                box_dst_.setLow(i,box_dst.getLow(i) - shift);
                box_src_.setLow(i,box_src.getLow(i) - shift);
            }
            else
            {
                box_dst_.setLow(i,box_dst.getLow(i));
                box_src_.setLow(i,box_src.getLow(i));
            }
        }
 
        typedef typename to_int_sequence<0,T::max_prop>::type result;
 
        copy_grid_fast_caller<result>::call(*this,grid_src,box_src_,box_dst_);
 
/*        copy_grid_fast<!is_contiguos<prp...>::type::value || has_pack_gen<typename device_grid::value_type>::value,
                                   dim,
                                   grid_cp,
                                   grid_info_cp>::copy(grid_src.getGrid(),
                                   this->getGrid(),
                                   box_src,
                                   box_dst,
                                   grid_src,*this,cnt);*/
 
    }
 
    template<unsigned int ... prp>
    void copy_to_prp(const grid_base_impl<dim,T,S,layout_base> & grid_src,
                 const Box<dim,size_t> & box_src,
                 const Box<dim,size_t> & box_dst)
    {
        typedef typename std::remove_reference<decltype(grid_src)>::type grid_cp;
        typedef typename std::remove_reference<decltype(grid_src.getGrid())>::type grid_info_cp;
 
        grid_key_dx<dim> cnt[1];
        cnt[0].zero();
 
        copy_grid_fast<!is_contiguos<prp...>::type::value || has_pack_gen<T>::value,
                                   dim,
                                   grid_cp,
                                   grid_info_cp>::copy(grid_src.getGrid(),
                                   this->getGrid(),
                                   box_src,
                                   box_dst,
                                   grid_src,*this,cnt);
    }
 
    void clear()
    {}
 
    template<template<typename,typename> class op, unsigned int ... prp>
    void copy_to_op(const grid_base_impl<dim,T,S,layout_base> & gs,
                 const Box<dim,size_t> & bx_src,
                 const Box<dim,size_t> & bx_dst)
    {
        grid_key_dx_iterator_sub<dim> sub_src(gs.getGrid(),bx_src.getKP1(),bx_src.getKP2());
        grid_key_dx_iterator_sub<dim> sub_dst(this->getGrid(),bx_dst.getKP1(),bx_dst.getKP2());
 
        while (sub_src.isNext())
        {
            // write the object in the last element
            object_si_di_op<op,decltype(gs.get_o(sub_src.get())),decltype(this->get_o(sub_dst.get())),OBJ_ENCAP,prp...>(gs.get_o(sub_src.get()),this->get_o(sub_dst.get()));
 
            ++sub_src;
            ++sub_dst;
        }
    }
 
    static bool is_unpack_header_supported()
    {return false;}
 
    void resize(const size_t (& sz)[dim], size_t opt = DATA_ON_HOST | DATA_ON_DEVICE, unsigned int blockSize = 1)
    {
 
        grid_base_impl<dim,T,S,layout_base,ord_type> grid_new(sz);
 
        resize_impl_memset(grid_new);
 
 
        if (opt & DATA_ON_HOST)
        {resize_impl_host(sz,grid_new);}
 
        if (opt & DATA_ON_DEVICE && S::isDeviceHostSame() == false)
        {resize_impl_device(sz,grid_new,blockSize);}
 
//      }
 
        // copy grid_new to the base
 
        this->swap(grid_new);
    }
 
    void resize_no_device(const size_t (& sz)[dim])
    {
 
        grid_base_impl<dim,T,S,layout_base> grid_new(sz);
 
        resize_impl_memset(grid_new);
        resize_impl_host(sz,grid_new);
 
        this->swap(grid_new);
    }
 
    void remove(size_t key)
    {
        if (dim != 1)
        {
#ifdef SE_CLASS1
            std::cerr << "Error: " << __FILE__ << " " << __LINE__ << " remove work only on dimension == 1 " << "\n";
#endif
            return;
        }
 
        // It is safe to do a memory copy
 
        data_.move(&this->template get<0>());
    }
 
    void swap_nomode(grid_base_impl<dim,T,S,layout_base> & grid)
    {
        mem_swap<T,layout_base<T>,decltype(data_),decltype(grid)>::template swap_nomode<S>(data_,grid.data_);
 
        // exchange the grid info
        g1.swap(grid.g1);
 
        // exchange the init status
        bool exg = is_mem_init;
        is_mem_init = grid.is_mem_init;
        grid.is_mem_init = exg;
 
#if defined(CUDIFY_USE_SEQUENTIAL) || defined(CUDIFY_USE_OPENMP)
 
        base_gpu = grid_toKernelImpl<is_layout_inte<layout_base<T_>>::value,grid_gpu_ker<dim,T_,layout_base,linearizer_type>,dim,T_,S>::toKernel(*this);
 
#endif
    }
 
    void swap(grid_base_impl<dim,T,S,layout_base,ord_type> & grid)
    {
        mem_swap<T,layout_base<T>,decltype(data_),decltype(grid)>::swap(data_,grid.data_);
 
        // exchange the grid info
        g1.swap(grid.g1);
 
        // exchange the init status
        bool exg = is_mem_init;
        is_mem_init = grid.is_mem_init;
        grid.is_mem_init = exg;
 
        // exchange the is external status
        exg = isExternal;
        isExternal = grid.isExternal;
        grid.isExternal = exg;
 
#if defined(CUDIFY_USE_SEQUENTIAL) || defined(CUDIFY_USE_OPENMP)
 
        base_gpu = grid_toKernelImpl<is_layout_inte<layout_base<T_>>::value,grid_gpu_ker<dim,T_,layout_base,linearizer_type>,dim,T_,S>::toKernel(*this);
        grid.base_gpu = grid_toKernelImpl<is_layout_inte<layout_base<T_>>::value,grid_gpu_ker<dim,T_,layout_base,linearizer_type>,dim,T_,S>::toKernel(grid);
 
#endif
    }
 
    void swap(grid_base_impl<dim,T,S,layout_base> && grid)
    {
        swap(grid);
    }
 
    template<unsigned int ... prp>
    __device__ __host__ inline void set(const grid_key_dx<dim> & key1,const grid_base_impl & g, const grid_key_dx<dim> & key2)
    {
        auto edest = this->get_o(key1);
        auto esrc = g.get_o(key2);
 
        copy_cpu_encap_encap_prp<decltype(g.get_o(key2)),decltype(this->get_o(key1)),prp...> ec(esrc,edest);
 
        boost::mpl::for_each_ref<boost::mpl::range_c<int,0,sizeof...(prp)>>(ec);
    }
 
    template<typename Memory> inline void set(grid_key_dx<dim> dx, const encapc<1,T,Memory> & obj)
    {
#ifdef SE_CLASS1
        check_init();
        check_bound(dx);
#endif
 
        // create the object to copy the properties
        copy_cpu_encap<dim,grid_base_impl<dim,T,S,layout_base>,layout> cp(dx,*this,obj);
 
        // copy each property
        boost::mpl::for_each_ref< boost::mpl::range_c<int,0,T::max_prop> >(cp);
    }
 
    inline void set(grid_key_dx<dim> dx, const T & obj)
    {
#ifdef SE_CLASS1
        check_init();
        check_bound(dx);
#endif
 
        this->get_o(dx) = obj;
    }
 
 
    inline void set(const grid_key_dx<dim> & key1,
                    const grid_base_impl<dim,T,S,layout_base> & g,
                    const grid_key_dx<dim> & key2)
    {
#ifdef SE_CLASS1
        check_init();
        check_bound(key1);
        check_bound(g,key2);
#endif
 
        this->get_o(key1) = g.get_o(key2);
    }
 
    inline void set(const size_t key1,
                    const grid_base_impl<dim,T,S,layout_base> & g,
                    const size_t key2)
    {
#ifdef SE_CLASS1
        check_init();
        check_bound(key1);
        check_bound(g,key2);
#endif
 
        this->get_o(key1) = g.get_o(key2);
    }
 
    template<typename Mem> inline void set(const grid_key_dx<dim> & key1,const grid_base_impl<dim,T,Mem,layout_base> & g, const grid_key_dx<dim> & key2)
    {
#ifdef SE_CLASS1
        check_init();
        check_bound(key1);
        check_bound(g,key2);
#endif
 
        this->get_o(key1) = g.get_o(key2);
    }
 
    template<typename Mem, template <typename> class layout_base2> inline void set_general(const grid_key_dx<dim> & key1,
                                                   const grid_base_impl<dim,T,Mem,layout_base2> & g,
                                                   const grid_key_dx<dim> & key2)
    {
#ifdef SE_CLASS1
        check_init();
        check_bound(key1);
        check_bound(g,key2);
#endif
 
        this->get_o(key1) = g.get_o(key2);
    }
 
    inline size_t size() const
    {
        return g1.size();
    }
 
    inline grid_key_dx_iterator_sub<dim> getSubIterator(const grid_key_dx<dim> & start, const grid_key_dx<dim> & stop) const
    {
        return g1.getSubIterator(start,stop);
    }
 
    inline grid_key_dx_iterator_sub<dim> getSubIterator(size_t m)
    {
        return grid_key_dx_iterator_sub<dim>(g1,m);
    }
 
    inline grid_key_dx_iterator<dim> getIterator() const
    {
        size_t sz[dim];
 
        for (int i = 0 ; i < dim ; i++)
        {sz[i] = g1.size(i);}
 
        grid_sm<dim,void> gvoid(sz);
 
        return grid_key_dx_iterator<dim>(gvoid);
    }
 
    template<unsigned int ... prp> void deviceToHost()
    {
        layout_base<T>::template deviceToHost<decltype(data_), prp ...>(data_,0,this->getGrid().size() - 1);
    }
 
    template<unsigned int ... prp> void deviceToHost(size_t start, size_t stop)
    {
        layout_base<T>::template deviceToHost<decltype(data_), prp ...>(data_,start,stop);
    }
 
    template<unsigned int ... prp> void hostToDeviceNUMA(size_t start, size_t stop)
    {
        #ifdef CUDIFY_USE_OPENMP
        grid_key_dx<dim> start_;
        grid_key_dx<dim> stop_;
 
        for (size_t i = 0 ; i < dim ; i++)
        {
            start_.set_d(i,start);
            stop_.set_d(i,stop);
        }
        
        auto ite = this->getGPUIterator(start_,stop_);
 
        // We have to carefull with openmp, numtiple thread can end up in numa
        copy_with_openmp_prp<decltype(this->toKernel()),typename std::remove_reference<decltype(*this)>::type,prp ...>(this->toKernel(),*this,ite);
        #else
        this->template hostToDevice<prp ...>(start,stop);
        #endif
    }
 
    template<unsigned int ... prp> void hostToDeviceNUMA()
    {
        #ifdef CUDIFY_USE_OPENMP
 
        grid_key_dx<dim> start_;
        grid_key_dx<dim> stop_;
 
        for (size_t i = 0 ; i < dim ; i++)
        {
            start_.set_d(i,0);
            stop_.set_d(i,this->g1.size() - 1);
        }
        
        auto ite = this->getGPUIterator(start_,stop_);
 
        // We have to carefull with openmp, numtiple thread can end up in numa
        copy_with_openmp_prp<decltype(this->toKernel()),typename std::remove_reference<decltype(*this)>::type,prp ...>(this->toKernel(),*this,ite);
        #else
        this->template hostToDevice<prp ...>();
        #endif
    }
 
    template<unsigned int ... prp> void hostToDevice(size_t start, size_t stop)
    {
        layout_base<T>::template hostToDevice<S, decltype(data_),prp ...>(data_,start,stop);
    }
 
 
    template<unsigned int ... prp> void hostToDevice()
    {
        layout_base<T>::template hostToDevice<S,decltype(data_),prp ...>(data_,0,this->getGrid().size() - 1);
    }
 
#if defined(CUDIFY_USE_SEQUENTIAL) || defined(CUDIFY_USE_OPENMP)
 
    mutable grid_gpu_ker<dim,T_,layout_base,linearizer_type> base_gpu;
 
    grid_gpu_ker_ref<dim,T_,layout_base,linearizer_type> toKernel()
    {
        return grid_gpu_ker_ref<dim,T_,layout_base,linearizer_type>(base_gpu);
    }
 
    const grid_gpu_ker_ref<dim,T_,layout_base,linearizer_type> toKernel() const
    {
        return grid_gpu_ker_ref<dim,T_,layout_base,linearizer_type>(base_gpu);
    }
 
#else 
 
    grid_gpu_ker<dim,T_,layout_base,linearizer_type> toKernel()
    {
        return grid_toKernelImpl<is_layout_inte<layout_base<T_>>::value,grid_gpu_ker<dim,T_,layout_base,linearizer_type>,dim,T_,S>::toKernel(*this);
    }
 
    const grid_gpu_ker<dim,T_,layout_base,linearizer_type> toKernel() const
    {
        return grid_toKernelImpl<is_layout_inte<layout_base<T_>>::value,grid_gpu_ker<dim,T_,layout_base,linearizer_type>,dim,T_,S>::toKernel(*this);
    }
 
#endif
 
    template<unsigned int Np>
    inline grid_key_dx_iterator<dim,stencil_offset_compute<dim,Np>>
    getIteratorStencil(const grid_key_dx<dim> (& stencil_pnt)[Np]) const
    {
        return grid_key_dx_iterator<dim,stencil_offset_compute<dim,Np>>(g1,stencil_pnt);
    }
 
    inline grid_key_dx_iterator_sub<dim> getIterator(const grid_key_dx<dim> & start, const grid_key_dx<dim> & stop, bool to_init = false) const
    {
        size_t sz[dim];
 
        for (int i = 0 ; i < dim ; i++)
        {sz[i] = g1.size(i);}
 
        grid_sm<dim,void> gvoid(sz);
 
        // get the starting point and the end point of the real domain
        return grid_key_dx_iterator_sub<dim>(gvoid,start,stop);
    }
 
    layout & get_internal_data_()
    {
        return data_;
    }
 
    const layout & get_internal_data_() const
    {
        return data_;
    }
 
    void removeAddUnpackReset()
    {}
 
    template<unsigned int ... prp, typename context_type>
    void removeAddUnpackFinalize(const context_type& gpuContext, int opt)
    {}
 
    template<unsigned int ... prp, typename context_type>
    void removeCopyToFinalize(const context_type& gpuContext, int opt)
    {}
 
    void resetFlush()
    {}
};
 
 
#endif /* OPENFPM_DATA_SRC_GRID_GRID_BASE_IMPLEMENTATION_HPP_ */