grid_sm.hpp

#ifndef GRID_HPP
#define GRID_HPP
 
 
#include "config.h"
#include "util/cuda_util.hpp"
#include <boost/shared_array.hpp>
#include <vector>
#include <initializer_list>
#include <array>
#include "memory/memory.hpp"
#include "Space/Shape/Box.hpp"
#include "grid_key.hpp"
#include <iostream>
#include "util/mathutil.hpp"
#include "iterators/stencil_type.hpp"
 
 
// Box need the definition of grid_key_dx_r
#define HARDWARE 1
 
 
struct No_check
{
    template<typename SparseGrid_type>
    __device__ __host__ bool check(SparseGrid_type & sggt,unsigned int dataBlockPos, unsigned int offset)
    {
        return true;
    }
};
 
template<unsigned int dim, typename T>
struct Box_check
{
    Box<dim,T> box;
 
    template<typename T2>
    explicit Box_check(Box<dim,T2> & box)
    :box(box)
    {}
 
    template<typename SparseGridGpu_type>
    __device__ __host__ bool check(SparseGridGpu_type & sggt,unsigned int dataBlockId, unsigned int offset)
    {
        auto key = sggt.getCoord(dataBlockId,offset);
 
        return box.isInsideKey(key);
    }
};
 
 
class NoCheck
{
public:
    static bool valid(size_t v_id, size_t sz)
    {
        return true;
    }
};
 
class CheckExistence
{
public:
    static bool valid(size_t v_id, size_t sz)
    {
        return v_id < sz;
    }
};
 
template<unsigned int dim>
struct ite_gpu
{
#if defined(CUDA_GPU)
 
    dim3 thr;
    dim3 wthr;
 
    grid_key_dx<dim,int> start;
    grid_key_dx<dim,int> stop;
 
    size_t nblocks() const
    {
        return wthr.x * ((wthr.y != 0)?wthr.y:1) * ((wthr.z != 0)?wthr.z:1);
    }
 
    size_t nthrs() const
    {
        return thr.x * ((thr.y != 0)?thr.y:1) * ((thr.z != 0)?thr.z:1);
    }
 
#endif
};
 
template<unsigned int dim>
bool has_work_gpu(ite_gpu<dim> & ite)
{
    size_t tot_work = 1;
 
    if (dim == 1)
    {tot_work *= ite.wthr.x * ite.thr.x;}
    else if(dim == 2)
    {
        tot_work *= ite.wthr.x * ite.thr.x;
        tot_work *= ite.wthr.y * ite.thr.y;
    }
    else
    {
        tot_work *= ite.wthr.x * ite.thr.x;
        tot_work *= ite.wthr.y * ite.thr.y;
        tot_work *= ite.wthr.z * ite.thr.z;
    }
 
    return tot_work != 0;
}
 
template<unsigned int N, typename T> class grid_sm;
 
template<unsigned int dim, typename grid_sm_type, typename T>
ite_gpu<dim> getGPUIterator_impl(const grid_sm_type & g1, const grid_key_dx<dim,T> & key1, const grid_key_dx<dim,T> & key2, size_t n_thr = default_kernel_wg_threads_);
 
template <unsigned int dim, typename linearizer> class print_warning_on_adjustment;
 
template<unsigned int dim,typename stencil=no_stencil, typename linearizer = grid_sm<dim,void>, typename warn=print_warning_on_adjustment<dim,linearizer>> class grid_key_dx_iterator_sub;
 
template<unsigned int N, typename T>
class grid_sm
{
    Box<N,size_t> box;
 
    size_t size_tot;
 
    size_t sz[N];
 
    size_t sz_s[N];
 
    __device__ __host__ inline void Initialize(const size_t sz)
    {
        sz_s[0] = sz;
        this->sz[0] = sz;
 
        // set the box
        box.setHigh(0,sz);
        box.setLow(0,0);
 
        for (size_t i = 1 ;  i < N ; i++)
        {
            /* coverity[dead_error_begin] */
            sz_s[i] = sz*sz_s[i-1];
            this->sz[i] = sz;
 
            // set the box
            box.setHigh(i,sz);
            box.setLow(i,0);
        }
    }
 
    __device__ __host__ inline void Initialize(const size_t (& sz)[N])
    {
        sz_s[0] = sz[0];
        this->sz[0] = sz[0];
 
        // set the box
        box.setHigh(0,sz[0]);
        box.setLow(0,0);
 
        for (size_t i = 1 ;  i < N ; i++)
        {
            /* coverity[dead_error_begin] */
            sz_s[i] = sz[i]*sz_s[i-1];
            this->sz[i] = sz[i];
 
            // set the box
            box.setHigh(i,sz[i]);
            box.setLow(i,0);
        }
    }
 
    inline void Initialize()
    {
        sz_s[0] = 0;
        this->sz[0] = 0;
 
        // set the box
        box.setHigh(0,0);
        box.setLow(0,0);
 
        for (size_t i = 1 ;  i < N ; i++)
        {
            /* coverity[dead_error_begin] */
            sz_s[i] = sz[i]*sz_s[i-1];
 
            // set the box
            box.setHigh(i,sz[i]);
            box.setLow(i,0);
        }
    }
 
    template<typename a, typename ...lT>
    __device__ __host__ inline mem_id Lin_impl(a v,lT...t) const
    {
        return v*sz_s[sizeof...(t)-1] + Lin_impl(t...);
    }
 
    template<typename a> __device__ __host__ inline mem_id Lin_impl(a v) const
    {
        return v;
    }
 
public:
 
 
    inline Box<N,size_t> getBox() const
    {
        return box;
    }
 
    inline const Box<N,size_t> getBoxKey() const
    {
        Box<N,size_t> bx;
 
        for (size_t i = 0 ; i < N ; i++)
        {
            bx.setLow(i,box.getLow(i));
            bx.setHigh(i,box.getHigh(i) - 1);
        }
 
        return bx;
    }
 
    inline void setDimensions(const size_t  (& dims)[N])
    {
        Initialize(dims);
        size_tot = totalSize(dims);
    }
 
    inline grid_sm()
    :size_tot(0)
    {
        // Initialize sz
        for (size_t i = 0 ; i < N ; i++)
        {sz[i] = 0;}
 
        Initialize();
    }
 
    template<typename S>
    __device__ __host__ inline grid_sm(const grid_sm<N,S> & g)
    {
        size_tot = g.size_tot;
 
        for (size_t i = 0 ; i < N ; i++)
        {
            sz[i] = g.sz[i];
            sz_s[i] = g.sz_s[i];
        }
    }
 
    __device__ __host__ inline grid_sm(const grid_sm<N,T> & g)
    {
        size_tot = g.size_tot;
 
        for (size_t i = 0 ; i < N ; i++)
        {
            sz[i] = g.sz[i];
            sz_s[i] = g.sz_s[i];
        }
    }
 
    // Static element to calculate total size
 
    inline size_t totalSize(const size_t sz)
    {
        size_t tSz = 1;
 
        for (size_t i = 0 ;  i < N ; i++)
        {
            tSz *= sz;
        }
 
        return tSz;
    }
 
    // Static element to calculate total size
 
    __device__ __host__ inline size_t totalSize(const size_t (& sz)[N])
    {
        size_t tSz = 1;
 
        for (size_t i = 0 ;  i < N ; i++)
        {
            tSz *= sz[i];
        }
 
        return tSz;
    }
 
 
    inline grid_sm(const size_t & sz)
    : size_tot(totalSize(sz))
    {
        Initialize(sz);
    }
 
    __device__ __host__ inline grid_sm(const size_t (& sz)[N])
    : size_tot(totalSize(sz))
    {
        Initialize(sz);
    }
 
    template<typename check=NoCheck, typename ids_type>
    inline mem_id LinId(const grid_key_dx<N,ids_type> & gk, const signed char sum_id[N]) const
    {
        mem_id lid;
 
        // Check the sum produce a valid key
 
        if (check::valid(gk.get(0) + sum_id[0],sz[0]) == false)
            return -1;
 
        lid = gk.get(0) + sum_id[0];
 
 
        for (mem_id i = 1 ; i < N ; i++)
        {
            // Check the sum produce a valid key
 
            if (check::valid(gk.get(i) + sum_id[i],sz[i]) == false)
                return -1;
 
            lid += (gk.get(i) + sum_id[i]) * sz_s[i-1];
        }
 
        return lid;
    }
 
    template<typename check=NoCheck,typename ids_type>
    inline mem_id LinId(const grid_key_dx<N,ids_type> & gk, const signed char sum_id[N], const size_t (&bc)[N]) const
    {
        mem_id lid;
 
        // Check the sum produce a valid key
 
        if (bc[0] == NON_PERIODIC)
        {
            if (check::valid(gk.get(0) + sum_id[0],sz[0]) == false)
                return -1;
 
            lid = gk.get(0) + sum_id[0];
        }
        else
        {
            lid = openfpm::math::positive_modulo(gk.get(0) + sum_id[0],sz[0]);
        }
 
        for (mem_id i = 1 ; i < N ; i++)
        {
            // Check the sum produce a valid key
 
            /* coverity[dead_error_line] */
            if (bc[i] == NON_PERIODIC)
            {
                if (check::valid(gk.get(i) + sum_id[i],sz[i]) == false)
                    return -1;
 
                lid += (gk.get(i) + sum_id[i]) * sz_s[i-1];
            }
            else
            {
                lid += (openfpm::math::positive_modulo(gk.get(i) + sum_id[i],sz[i])) * sz_s[i-1];
            }
        }
 
        return lid;
    }
 
    inline mem_id LinIdPtr(size_t * k) const
    {
        mem_id lid = k[0];
        for (mem_id i = 1 ; i < N ; i++)
        {
            lid += k[i] * sz_s[i-1];
        }
 
        return lid;
    }
 
    __device__ __host__ inline mem_id LinId(const size_t (& k)[N]) const
    {
        mem_id lid = k[0];
        for (mem_id i = 1 ; i < N ; i++)
        {
            /* coverity[dead_error_line] */
            lid += k[i] * sz_s[i-1];
        }
 
        return lid;
    }
 
    template<typename ids_type> __device__ __host__ inline mem_id LinId(const grid_key_dx<N,ids_type> & gk) const
    {
        mem_id lid = gk.get(0);
        for (mem_id i = 1 ; i < N ; i++)
        {
            /* coverity[dead_error_line] */
            lid += gk.get(i) * sz_s[i-1];
        }
 
        return lid;
    }
 
    template<typename a, typename ...lT, typename enabler = typename std::enable_if<sizeof...(lT) == N-1>::type >
    __device__ __host__ inline mem_id Lin(a v,lT...t) const
    {
        return v*sz_s[sizeof...(t)-1] + Lin_impl(t...);
    }
 
 
    __device__ __host__ inline grid_key_dx<N> InvLinId(mem_id id) const
    {
        // Inversion of linearize
 
        grid_key_dx<N> gk;
 
        for (mem_id i = 0 ; i < N ; i++)
        {
            gk.set_d(i,id % sz[i]);
            id /= sz[i];
        }
 
        return gk;
    }
 
 
    //#pragma openfpm layout(get)
    inline mem_id LinId(mem_id * id) const
    {
        mem_id lid = 0;
        lid += id[0];
        for (mem_id i = 1 ; i < N ; i++)
        {
            lid += id[i] * sz_s[i-1];
        }
 
        return lid;
    }
 
    __device__ __host__ ~grid_sm() {};
 
    __device__ __host__ inline size_t size() const
    {
        return size_tot;
    };
 
    __device__ __host__ inline grid_sm<N,T> & operator=(const grid_sm<N,T> & g)
    {
        size_tot = g.size_tot;
 
        for (size_t i = 0 ; i < N ; i++)
        {
            sz[i] = g.sz[i];
            sz_s[i] = g.sz_s[i];
        }
 
        box = g.box;
 
        return *this;
    }
 
    inline bool operator==(const grid_sm<N,T> & g)
    {
        for (size_t i = 0 ; i < N ; i++)
        {
            if (sz[i] != g.sz[i])
                return false;
        }
 
#ifdef SE_CLASS1
 
        if (size_tot != g.size_tot)
            return false;
 
        for (size_t i = 0 ; i < N ; i++)
        {
            if (sz_s[i] != g.sz_s[i])
                return false;
        }
 
#endif
        return true;
    }
 
    inline bool operator!=(const grid_sm<N,T> & g)
    {
        return ! this->operator==(g);
    }
 
    __device__ __host__ inline size_t size_s(unsigned int i) const
    {
        return sz_s[i];
    }
 
    __device__ __host__ inline size_t size(unsigned int i) const
    {
        return sz[i];
    }
 
    __device__ __host__ inline const size_t (& getSize() const)[N]
    {
        return sz;
    }
 
    __device__ __host__ inline void getSize(size_t (&size) [N]) const
    {
        for (size_t i = 0; i < N; ++i) size[i] = sz[i];
    }
 
    inline grid_key_dx_iterator_sub<N> getSubIterator(const grid_key_dx<N> & start, const grid_key_dx<N> & stop) const
    {
        return grid_key_dx_iterator_sub<N>(*this,start,stop);
    }
 
#if defined(CUDA_GPU)
 
    template<typename T2>
    struct ite_gpu<N> getGPUIterator(const grid_key_dx<N,T2> & key1, const grid_key_dx<N,T2> & key2, size_t n_thr = default_kernel_wg_threads_) const
    {
        return getGPUIterator_impl<N>(*this,key1,key2,n_thr);
    }
 
    struct ite_gpu<N> getGPUIterator(size_t n_thr = default_kernel_wg_threads_) const
    {
        grid_key_dx<N> k1;
        grid_key_dx<N> k2;
 
        for (size_t i = 0 ; i < N ; i++)
        {
            k1.set_d(i,0);
            k2.set_d(i,size(i));
        }
 
        return getGPUIterator_impl<N>(*this,k1,k2,n_thr);
    }
 
 
#endif
 
    inline void swap(grid_sm<N,T> & g)
    {
        size_t tmp = size_tot;
        size_tot = g.size_tot;
        g.size_tot = tmp;
 
        for (size_t i = 0 ; i < N ; i++)
        {
            tmp = sz[i];
            sz[i] = g.sz[i];
            g.sz[i] = tmp;
 
            tmp = sz_s[i];
            sz_s[i] = g.sz_s[i];
            g.sz_s[i] = tmp;
        }
    }
 
    std::string toString() const
    {
        std::stringstream str;
 
        for (size_t i = 0 ; i < N ; i++)
            str << "sz[" << i << "]=" << size(i) << " ";
 
        return str.str();
    }
 
    template <unsigned int,typename> friend class grid_sm;
};
 
 
template<unsigned int dim, typename grid_sm_type, typename T>
ite_gpu<dim> getGPUIterator_impl(const grid_sm_type & g1, const grid_key_dx<dim,T> & key1, const grid_key_dx<dim,T> & key2, const size_t n_thr)
{
    // Work to do
    ite_gpu<dim> ig;
 
    if (n_thr == static_cast<size_t>(-1)) {
        ig.thr.x = 0xFFFFFFFF;
        ig.thr.y = 0xFFFFFFFF;
        ig.thr.z = 0xFFFFFFFF;
 
        if (dim >= 1)
            ig.wthr.x = key2.get(0) - key1.get(0) + 1;
        else
            ig.wthr.x = 0;
 
        if (dim >= 2)
            ig.wthr.y = key2.get(1) - key1.get(1) + 1;
        else
            ig.wthr.y = 0;
 
        if (dim >= 3)
        {
            ig.wthr.z = 1;
            for (size_t i = 2 ; i < dim ; i++)
            {ig.wthr.z *= (key2.get(i) - key1.get(i) + 1);}
        }
        else
            ig.wthr.z = 0;
 
        for (size_t i = 0 ; i < dim ; i++)
        {
            ig.start.set_d(i,key1.get(i));
            ig.stop.set_d(i,key2.get(i));
        }
 
        return ig;
    
    }
 
    size_t tot_work = 1;
    for (size_t i = 0 ; i < dim ; i++)
    {tot_work *= key2.get(i) - key1.get(i) + 1;}
 
    size_t n = (tot_work <= n_thr)?openfpm::math::round_big_2(tot_work):n_thr;
 
    if (tot_work == 0)
    {
        ig.thr.x = 0;
        ig.thr.y = 0;
        ig.thr.z = 0;
 
        ig.wthr.x = 0;
        ig.wthr.y = 0;
        ig.wthr.z = 0;
 
        for (size_t i = 0 ; i < dim ; i++)
        {
            ig.start.set_d(i,key1.get(i));
            ig.stop.set_d(i,key2.get(i));
        }
 
        return ig;
    }
 
    ig.thr.x = 1;
    ig.thr.y = 1;
    ig.thr.z = 1;
 
    int dir = 0;
 
    while (n != 1)
    {
        if (dir % 3 == 0)
        {ig.thr.x = ig.thr.x << 1;}
        else if (dir % 3 == 1)
        {ig.thr.y = ig.thr.y << 1;}
        else if (dir % 3 == 2)
        {ig.thr.z = ig.thr.z << 1;}
 
        n = n >> 1;
 
        dir++;
        dir %= dim;
    }
 
    if (dim >= 1)
    {ig.wthr.x = (key2.get(0) - key1.get(0) + 1) / ig.thr.x + (((key2.get(0) - key1.get(0) + 1)%ig.thr.x != 0)?1:0);}
 
    if (dim >= 2)
    {ig.wthr.y = (key2.get(1) - key1.get(1) + 1) / ig.thr.y + (((key2.get(1) - key1.get(1) + 1)%ig.thr.y != 0)?1:0);}
    else
    {ig.wthr.y = 1;}
 
    if (dim >= 3)
    {
        // Roll the other dimensions on z
        ig.wthr.z = 1;
        for (size_t i = 2 ; i < dim ; i++)
        {ig.wthr.z *= (key2.get(i) - key1.get(i) + 1) / ig.thr.z + (((key2.get(i) - key1.get(i) + 1)%ig.thr.z != 0)?1:0);}
    }
    else
    {ig.wthr.z = 1;}
 
    // crop if wthr == 1
 
    if (dim >= 1 && ig.wthr.x == 1)
    {ig.thr.x = (key2.get(0) - key1.get(0) + 1);}
 
    if (dim >= 2 && ig.wthr.y == 1)
    {ig.thr.y = key2.get(1) - key1.get(1) + 1;}
 
    if (dim == 3 && ig.wthr.z == 1)
    {ig.thr.z = key2.get(2) - key1.get(2) + 1;}
 
    for (size_t i = 0 ; i < dim ; i++)
    {
        ig.start.set_d(i,key1.get(i));
        ig.stop.set_d(i,key2.get(i));
    }
 
    return ig;
}
 
 
class grid_key_dx_r
{
    size_t dim;
 
public:
 
    size_t getDim()
    {
        return dim;
    }
 
    grid_key_dx_r(grid_key_dx_r & key)
    :dim(key.dim)
    {
        // Allocate the key
        k = new mem_id[dim];
 
        // Copy the key
        for(unsigned int i = 0 ; i < dim ; i++)
        {
            k[i] = key.k[i];
        }
    }
 
    grid_key_dx_r(size_t dim)
    :dim(dim)
    {
        // Allocate the key
        k = new mem_id[dim];
    }
 
    ~grid_key_dx_r()
    {
        delete [] k;
    }
 
    template<typename a, typename ...T>void set(a v, T...t)
    {
        k[dim-1] = v;
        invert_assign(t...);
    }
 
    mem_id get(size_t i)
    {
        return k[i];
    }
 
    void set_d(size_t i, mem_id id)
    {
        k[i] = id;
    }
 
    mem_id * k;
 
private:
 
    template<typename a, typename ...T>void invert_assign(a v,T...t)
    {
        k[sizeof...(T)] = v;
        invert_assign(t...);
    }
 
    template<typename a, typename ...T>void invert_assign(a v)
    {
        k[0] = v;
    }
 
    void invert_assign()
    {
    }
 
};
 
 
class Iterator_g_const
{
    size_t dim;
 
    size_t sz;
 
    // Actual grid_key position
    grid_key_dx_r gk;
 
public:
 
    size_t getDim()
    {
        return dim;
    }
 
    Iterator_g_const(size_t n, size_t sz)
    :dim(n),sz(sz),gk(n)
    {
        // fill gk with the first grid element that satisfied the constrain: 0,1,2,3... dim
 
        for (size_t i = 0 ; i < dim ; i++)
        {
            gk.set_d(i,dim-i-1);
        }
    }
 
    Iterator_g_const & operator++()
            {
 
        gk.set_d(0,gk.get(0)+1);
 
 
        unsigned int i = 0;
        for ( ; i < dim-1 ; i++)
        {
            size_t id = gk.get(i);
            if (id >= sz)
            {
                // ! overflow, increment the next index
 
                id = gk.get(i+1);
                if (id+i+2 >= sz)
                {
                    // there is no-way to produce a valid key
                    // there is not need to check the previous index
                    // overflow i+1
 
                    gk.set_d(i+1,sz);
                }
                else
                {
 
                    // reinitialize the previous index
 
                    for (unsigned int s = 0 ; s <= i+1 ; s++)
                    {
                        gk.set_d(i+1-s,id+1+s);
                    }
                }
            }
            else
            {
                break;
            }
        }
 
        return *this;
            }
 
    bool isNext()
    {
        // If dimensionless return immediately
        if (dim == 0)
            return false;
 
        if (gk.get(dim-1) < static_cast<mem_id>(sz-dim+1))
        {
 
            return true;
        }
 
        return false;
    }
 
    grid_key_dx_r & get()
    {
        return gk;
    }
};
 
#endif