grid_sm.hpp

#ifndef GRID_HPP
#define GRID_HPP


#include "config.h"
#include "util/cuda_launch.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()
    {
        return wthr.x * wthr.y * wthr.z;
    }

    size_t nthrs()
    {
        return thr.x * thr.y * thr.z;
    }

#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];

    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);
        }
    }

    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> 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];
        }
    }

    // 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

    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);
    }

    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 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 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;
    }

    ~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);
    }

    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];
    }

    inline const size_t (& getSize() const)[N]
    {
        return sz;
    }

    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)
{
    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;

    // Work to do
    ite_gpu<dim> ig;

    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;

        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