vector_dist_operators.hpp

/*
 * vector_dist_operators.hpp
 *
 *  Created on: Jun 11, 2016
 *      Author: i-bird
 */

#ifndef OPENFPM_NUMERICS_SRC_OPERATORS_VECTOR_VECTOR_DIST_OPERATORS_HPP_
#define OPENFPM_NUMERICS_SRC_OPERATORS_VECTOR_VECTOR_DIST_OPERATORS_HPP_

#include "Vector/vector_dist.hpp"
#include "lib/pdata.hpp"
#include "cuda/vector_dist_operators_cuda.cuh"

#define PROP_CUSTOM (unsigned int)-2

#define VECT_SUM 1
#define VECT_SUB 2
#define VECT_MUL 3
#define VECT_DIV 4

#define VECT_APPLYKER_IN 7
#define VECT_APPLYKER_OUT 8
#define VECT_APPLYKER_REDUCE 9
#define VECT_APPLYKER_IN_GEN 10
#define VECT_APPLYKER_OUT_GEN 11
#define VECT_APPLYKER_REDUCE_GEN 12
#define VECT_APPLYKER_IN_SIM 13
#define VECT_APPLYKER_OUT_SIM 14
#define VECT_APPLYKER_REDUCE_SIM 15

#define VECT_APPLYKER_IN_GEN_SORT 16
#define VECT_APPLYKER_IN_SORT 17
#define VECT_APPLYKER_IN_SIM_SORT 18

#define VECT_NORM 56
#define VECT_NORM2 57
#define VECT_ABS 58
#define VECT_EXP 59
#define VECT_EXP2 60
#define VECT_EXPM1 61
#define VECT_LOG 62
#define VECT_LOG10 63
#define VECT_LOG2 64
#define VECT_LOG1P 65
#define VECT_SQRT 67
#define VECT_CBRT 68
#define VECT_SIN 69
#define VECT_COS 70
#define VECT_TAN 71
#define VECT_ASIN 72
#define VECT_ACOS 73
#define VECT_ATAN 74
#define VECT_SINH 75
#define VECT_COSH 76
#define VECT_TANH 77
#define VECT_ASINH 78
#define VECT_ACOSH 79
#define VECT_ATANH 80
#define VECT_ERF 81
#define VECT_ERFC 82
#define VECT_TGAMMA 83
#define VECT_LGAMMA 84
#define VECT_CEIL 85
#define VECT_FLOOR 86
#define VECT_TRUNC 87
#define VECT_ROUND 88
#define VECT_NEARBYINT 89
#define VECT_RINT 90
#define VECT_PMUL 91
#define VECT_SUB_UNI 92
#define VECT_SUM_REDUCE 93
#define VECT_COMP 94
#define VECT_NORM_INF 95


#define VECT_DCPSE 100
#define VECT_DCPSE_V 101
#define VECT_DCPSE_V_SUM 102
#define VECT_DCPSE_V_DOT 103
#define VECT_DCPSE_V_DIV 104
#define VECT_DCPSE_V_CURL2D 105
#define VECT_COPY_1_TO_N 300
#define VECT_COPY_N_TO_N 301
#define VECT_COPY_N_TO_1 302
#define VECT_PMUL 91
#define VECT_SUB_UNI 92




template<bool cond, typename exp1, typename exp2>
struct first_or_second
{
    typedef typename exp2::vtype vtype;

    static auto getVector(const exp1 & o1, const exp2 & o2) -> decltype(o2.getVector())
    {
        return o2.getVector();
    }
};

template<typename exp1, typename exp2>
struct first_or_second<true,exp1,exp2>
{
    typedef typename exp1::vtype vtype;

    static auto getVector(const exp1 & o1, const exp2 & o2) -> decltype(o1.getVector())
    {
        return o1.getVector();
    }
};

template<typename T, typename Sfinae = void>
struct has_vtype: std::false_type {};

template<typename T>
struct has_vtype<T, typename Void<typename T::vtype>::type> : std::true_type
{};




template<typename ObjType, typename Sfinae = void>
struct has_init: std::false_type {};

template<typename ObjType>
struct has_init<ObjType, typename Void<typename ObjType::has_init>::type> : std::true_type
{};

template <typename T, bool has_init = has_init<T>::value >
struct call_init_if_needed
{
    static inline void call(T & r_exp)
    {
        r_exp.init();
    }
};

template <typename T>
struct call_init_if_needed<T,false>
{
    static inline void call(T & r_exp)
    {
    }
};



template <typename exp1, typename exp2, unsigned int op>
class vector_dist_expression_op
{

};

template<typename v1_type, typename v2_type>
struct vector_result
{
    typedef v1_type type;

    template<typename exp1, typename exp2>
    static const type & getVector(const exp1 & o1, const exp2 & o2)
    {
        return o1.getVector();
    }

    template<typename exp1>
    static const type & getVector(const exp1 & o1)
    {
        return o1.getVector();
    }
};


template<typename v2_type>
struct vector_result<void,v2_type>
{
    typedef v2_type type;

    template<typename exp1, typename exp2>
    static const type & getVector(const exp1 & o1, const exp2 & o2)
    {
        return o2.getVector();
    }

    template<typename exp2>
    static const type & getVector(exp2 & o2)
    {
        return o2.getVector();
    }
};

template<typename NN1_type, typename NN2_type>
struct nn_type_result
{
    typedef NN1_type type;

    template<typename exp1, typename exp2>
    static type * getNN(exp1 & o1, exp2 & o2)
    {
        return o1.getNN();
    }

    template<typename exp1>
    static type * getNN(exp1 & o1)
    {
        return o1.getNN();
    }
};


template<typename NN2_type>
struct nn_type_result<void,NN2_type>
{
    typedef NN2_type type;

    template<typename exp1, typename exp2>
    static type * getNN(exp1 & o1, exp2 & o2)
    {
        return o2.getNN();
    }

    template<typename exp2>
    static type * getNN(exp2 & o2)
    {
        return o2.getNN();
    }
};

template<bool s1, bool s2>
struct vector_is_sort_result
{
    typedef boost::mpl::bool_<s1 | s2> type;
};

template <typename exp1, typename exp2>
class vector_dist_expression_op<exp1,exp2,VECT_SUM>
{
    const exp1 o1;

    const exp2 o2;

public:

    typedef typename exp1::is_ker is_ker;

    typedef typename first_or_second<has_vtype<exp1>::value,exp1,exp2>::vtype vtype;

    typedef typename vector_is_sort_result<exp1::is_sort::value,exp2::is_sort::value>::type is_sort;

    typedef typename nn_type_result<typename exp1::NN_type,typename exp2::NN_type>::type NN_type;

    inline vector_dist_expression_op(const exp1 & o1, const exp2 & o2)
    :o1(o1),o2(o2)
    {}

    inline NN_type * getNN() const
    {
        return nn_type_result<typename exp1::NN_type,typename exp2::NN_type>::getNN(o1,o2);
    }

    inline void init() const
    {
        o1.init();
        o2.init();
    }

    template<typename r_type=typename std::remove_reference<decltype(o1.value(vect_dist_key_dx()) + o2.value(vect_dist_key_dx()))>::type >
    __device__ __host__ inline r_type value(const vect_dist_key_dx & key) const
    {
        return o1.value(key) + o2.value(key);
    }

    template<typename Sys_eqs, typename pmap_type, typename unordered_map_type, typename coeff_type>
    inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff, unsigned int comp) const
    {
        o1.template value_nz<Sys_eqs>(p_map,key,cols,coeff, comp);
        o2.template value_nz<Sys_eqs>(p_map,key,cols,coeff, comp);
    }

    const vtype & getVector()
    {
        return first_or_second<has_vtype<exp1>::value,exp1,exp2>::getVector(o1,o2);
    }

    const vtype & getVector() const
    {
            return first_or_second<has_vtype<exp1>::value,exp1,exp2>::getVector(o1,o2);
    }

    template<typename r_type=typename std::remove_reference<decltype(o1.value(0) + o2.value(0))>::type >
    __device__ __host__ inline r_type value(const unsigned int & key) const
    {
        return o1.value(key) + o2.value(key);
    }

};

template <typename exp1, typename exp2>
class vector_dist_expression_op<exp1,exp2,VECT_SUB>
{
    const exp1 o1;

    const exp2 o2;

public:

    typedef typename first_or_second<has_vtype<exp1>::value,exp1,exp2>::vtype vtype;

    typedef typename exp1::is_ker is_ker;

    typedef typename vector_is_sort_result<exp1::is_sort::value,exp2::is_sort::value>::type is_sort;

    typedef typename nn_type_result<typename exp1::NN_type,typename exp2::NN_type>::type NN_type;

    inline vector_dist_expression_op(const exp1 & o1, const exp2 & o2)
    :o1(o1),o2(o2)
    {}

    inline NN_type * getNN() const
    {
        return nn_type_result<typename exp1::NN_type,typename exp2::NN_type>::getNN(o1,o2);
    }


    inline void init() const
    {
        o1.init();
        o2.init();
    }

    template<typename r_type=typename std::remove_reference<decltype(o1.value(vect_dist_key_dx()) - o2.value(vect_dist_key_dx()))>::type >
    inline r_type value(const vect_dist_key_dx & key) const
    {
        return o1.value(key) - o2.value(key);
    }

    template<typename r_type=typename std::remove_reference<decltype(o1.value(vect_dist_key_dx()) - o2.value(vect_dist_key_dx()))>::type >
    __device__ __host__ inline r_type value(const unsigned int & key) const
    {
        return o1.value(key) - o2.value(key);
    }

    const vtype & getVector()
    {
        return first_or_second<has_vtype<exp1>::value,exp1,exp2>::getVector(o1,o2);
    }

    const vtype & getVector() const
    {
        return first_or_second<has_vtype<exp1>::value,exp1,exp2>::getVector(o1,o2);
    }

    template<typename Sys_eqs,typename pmap_type, typename unordered_map_type, typename coeff_type>
    inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff, unsigned int comp) const
    {
        o1.template value_nz<Sys_eqs>(p_map,key,cols,coeff,comp);
        coeff_type tmp = -coeff;
        o2.template value_nz<Sys_eqs>(p_map,key,cols,tmp,comp);
    }
};

template <typename exp1, typename exp2>
class vector_dist_expression_op<exp1,exp2,VECT_MUL>
{
    const exp1 o1;

    const exp2 o2;

public:

    typedef typename first_or_second<has_vtype<exp1>::value,exp1,exp2>::vtype vtype;

    typedef typename exp1::is_ker is_ker;

    typedef typename vector_is_sort_result<exp1::is_sort::value,exp2::is_sort::value>::type is_sort;

    typedef typename nn_type_result<typename exp1::NN_type,typename exp2::NN_type>::type NN_type;

    vector_dist_expression_op(const exp1 & o1, const exp2 & o2)
    :o1(o1),o2(o2)
    {}

    inline void init() const
    {
        o1.init();
        o2.init();
    }

    template<typename r_type=typename std::remove_reference<decltype(o1.value(vect_dist_key_dx()) * o2.value(vect_dist_key_dx()))>::type >
    __device__ __host__ inline r_type value(const vect_dist_key_dx & key) const
    {
        return o1.value(key) * o2.value(key);
    }

    template<typename Sys_eqs,typename pmap_type, typename unordered_map_type, typename coeff_type>
    inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff, unsigned int comp) const
    {
            //o1.template value_nz<Sys_eqs>(p_map,key,cols,coeff,comp);
            auto coeff_tmp = o1.value(key) * coeff;
            o2.template value_nz<Sys_eqs>(p_map,key,cols,coeff_tmp,comp);
        }

    const vtype & getVector()
    {
            return first_or_second<has_vtype<exp1>::value,exp1,exp2>::getVector(o1,o2);
    }

    const vtype & getVector() const
    {
            return first_or_second<has_vtype<exp1>::value,exp1,exp2>::getVector(o1,o2);
    }

    template<typename r_type=typename std::remove_reference<decltype(o1.value(vect_dist_key_dx()) * o2.value(vect_dist_key_dx()))>::type >
    __device__ __host__ inline r_type value(const unsigned int & key) const
    {
        return o1.value(key) * o2.value(key);
    }
};

template <typename exp1, typename exp2>
class vector_dist_expression_op<exp1,exp2,VECT_DIV>
{
    const exp1 o1;

    const exp2 o2;

public:

    typedef typename first_or_second<has_vtype<exp1>::value,exp1,exp2>::vtype vtype;

    typedef typename exp1::is_ker is_ker;

    typedef typename vector_is_sort_result<exp1::is_sort::value,exp2::is_sort::value>::type is_sort;

    typedef typename nn_type_result<typename exp1::NN_type,typename exp2::NN_type>::type NN_type;

    vector_dist_expression_op(const exp1 & o1, const exp2 & o2)
    :o1(o1),o2(o2)
    {}

    inline NN_type * getNN() const
    {
        return nn_type_result<typename exp1::NN_type,typename exp2::NN_type>::getNN(o1,o2);
    }


    inline void init() const
    {
        o1.init();
        o2.init();
    }

    template<typename r_type=typename std::remove_reference<decltype(o1.value(vect_dist_key_dx()) / o2.value(vect_dist_key_dx()))>::type > inline r_type value(const vect_dist_key_dx & key) const
    {
        return o1.value(key) / o2.value(key);
    }

    template<typename Sys_eqs,typename pmap_type, typename unordered_map_type, typename coeff_type>
    inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff, unsigned int comp) const
    {
        std::cout << __FILE__ << ":" << __LINE__ << " You are trying to divide by an operator,  this is not possible " << std::endl;
    }

    const vtype & getVector()
    {
        return first_or_second<has_vtype<exp1>::value,exp1,exp2>::getVector(o1,o2);
    }

    const vtype & getVector() const
    {
        return first_or_second<has_vtype<exp1>::value,exp1,exp2>::getVector(o1,o2);
    }

        template<typename r_type=typename std::remove_reference<decltype(o1.value(vect_dist_key_dx()) / o2.value(vect_dist_key_dx()))>::type >
        __device__ __host__ inline r_type value(const unsigned int & key) const
        {
                return o1.value(key) / o2.value(key);
        }
};


template <typename exp1>
class vector_dist_expression_op<exp1,void,VECT_SUB_UNI>
{
    const exp1 o1;

public:


    typedef typename exp1::is_ker is_ker;

    typedef typename vector_result<typename exp1::vtype,void>::type vtype;

    typedef typename vector_is_sort_result<exp1::is_sort::value,false>::type is_sort;

    typedef typename nn_type_result<typename exp1::NN_type,void>::type NN_type;

    vector_dist_expression_op(const exp1 & o1)
    :o1(o1)
    {}

    inline NN_type * getNN() const
    {
        return nn_type_result<typename exp1::NN_type,void>::getNN(o1);
    }

    const vtype & getVector()
    {
        return vector_result<typename exp1::vtype,void>::getVector(o1);
    }

    inline void init() const
    {
        o1.init();
    }

    const vtype & getVector() const
    {
        return o1.getVector();
    }


    template<typename r_type=typename std::remove_reference<decltype(-(o1.value(vect_dist_key_dx(0))))>::type > 
    __device__ __host__  inline r_type value(const vect_dist_key_dx & key) const
    {
        return -(o1.value(key));
    }

    template<typename Sys_eqs, typename pmap_type, typename unordered_map_type, typename coeff_type>
    inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff, unsigned int comp) const
    {
        coeff_type coeff_tmp = -coeff;
        o1.template value_nz<Sys_eqs>(p_map,key,cols,coeff_tmp, comp);
    }
    template<typename r_type=typename std::remove_reference<decltype(-(o1.value(vect_dist_key_dx(0))))>::type >
    __device__ __host__ inline r_type value(const unsigned int & key) const
    {
        return -(o1.value(key));
    }
};

template<int impl, bool vect_ker>
struct vector_dist_expression_comp_sel
{
    typedef boost::mpl::int_<impl> type;
};

template<>
struct vector_dist_expression_comp_sel<comp_host,true>
{
    typedef boost::mpl::int_<-1> type;
};

template<>
struct vector_dist_expression_comp_sel<comp_dev,false>
{
    typedef boost::mpl::int_<-1> type;
};

template<typename vector, bool is_ker = has_vector_kernel<vector>::type::value>
struct vector_expression_transform
{
    typedef vector& type;
};

template<typename vector>
struct vector_expression_transform<vector,true>
{
    typedef vector type;
};

template<typename vector>
struct v_mem_mutable
{
    vector v;

    v_mem_mutable(vector & v)
    :v(v)
    {}
};



template<unsigned int prp, typename vector>
class vector_dist_expression
{
    mutable v_mem_mutable<typename vector_expression_transform<vector>::type> v;

    vector_dist_ker_list<vector> * vdl;

public:

    typedef typename has_vector_kernel<vector>::type is_ker;

    typedef vector vtype;

    typedef boost::mpl::bool_<false> is_sort;

    typedef void NN_type;

    static const unsigned int prop = prp;

    int var_id = 0;

    void setVarId(int var_id)
    {
        this->var_id = var_id;
    }

    vector_dist_expression(vector & v)
    :v(v),vdl(NULL)
    {}

    ~vector_dist_expression()
    {
        if (vdl != NULL)
        {vdl->remove(v.v);}
    }

    inline void * getNN() const
    {
        return NULL;
    }

    __device__ __host__ const vector & getVector() const
    {
        return v.v;
    }

    __device__ __host__ vector & getVector()
    {
        return v.v;
    }

    void set_vector_dist_ker_list(vector_dist_ker_list<vector> & vdkl, bool is_sort)
    {
        vdkl.add(v.v,is_sort);
        vdl = &vdkl;
    }

    inline void init() const
    {}

    __host__ __device__ inline auto value(const vect_dist_key_dx & k) const -> decltype(pos_or_propR<vector,prp>::value(v.v,k))
    {
        return pos_or_propR<vector,prp>::value(v.v,k);
    }

    __host__ __device__ inline auto value(const vect_dist_key_dx & k) -> decltype(pos_or_propR<vector,prp>::value(v.v,k))
    {
        return pos_or_propR<vector,prp>::value(v.v,k);
    }

    __device__ __host__ inline auto value(const unsigned int & k) const -> decltype(pos_or_propR<vector,prp>::value(v.v,k))
    {
        return pos_or_propR<vector,prp>::value(v.v,k);
    }

    __device__ __host__ inline auto value(const unsigned int & k) -> decltype(pos_or_propR<vector,prp>::value(v.v,k))
    {
        return pos_or_propR<vector,prp>::value(v.v,k);
    }

    inline auto get(const vect_dist_key_dx & key) const -> decltype(value(key))
    {
        return this->value(key);
    }

    template<unsigned int prp2> vector & operator=(const vector_dist_expression<prp2,vector> & v_exp)
    {
        if (v_exp.getVector().isSubset() == true)
        {
                        std::cout << __FILE__ << ":" << __LINE__ << " error on the right hand side of the expression you have to use non-subset properties" << std::endl;
                        return v.v;
        }

        if (has_vector_kernel<vector>::type::value == false)
        {
            vector_dist_op_compute_op<prp,false,vector_dist_expression_comp_sel<comp_host,
                                                                          has_vector_kernel<vector>::type::value>::type::value>
            ::compute_expr(v.v,v_exp);
        }
        else
        {
            vector_dist_op_compute_op<prp,false,vector_dist_expression_comp_sel<comp_dev,
                                                                          has_vector_kernel<vector>::type::value>::type::value>
            ::compute_expr(v.v,v_exp);
        }

        return v.v;
    }

    template<typename T> vector & operator=(const vector_dist_expression<0,openfpm::vector<aggregate<T>>> & v_exp)
    {
        vector_dist_op_compute_op<prp,false,vector_dist_expression_comp_sel<comp_host,
                                                                          has_vector_kernel<vector>::type::value>::type::value>
        ::compute_expr(v.v,v_exp);

        return v.v;
    }

    template<typename exp1, typename exp2, unsigned int op>
    vector & operator=(const vector_dist_expression_op<exp1,exp2,op> & v_exp)
    {
        if (v_exp.getVector().isSubset() == true)
        {
            std::cout << __FILE__ << ":" << __LINE__ << " error on the right hand side of the expression you have to use non-subset properties" << std::endl;
            return v.v;
        }

        if (has_vector_kernel<vector>::type::value == false)
        {
            vector_dist_op_compute_op<prp,
                                      vector_dist_expression_op<exp1,exp2,op>::is_sort::value,
                                      vector_dist_expression_comp_sel<comp_host,
                                                                      has_vector_kernel<vector>::type::value>::type::value>
            ::compute_expr(v.v,v_exp);
        }
        else
        {
            vector_dist_op_compute_op<prp,
                                      vector_dist_expression_op<exp1,exp2,op>::is_sort::value,
                                      vector_dist_expression_comp_sel<comp_dev,
                                                                      has_vector_kernel<vector>::type::value>::type::value>
            ::compute_expr(v.v,v_exp);
        }

        return v.v;
    }

    vector & operator=(double d)
    {
        if (has_vector_kernel<vector>::type::value == false)
        {
            vector_dist_op_compute_op<prp,
                                      false,
                                      vector_dist_expression_comp_sel<comp_host,
                                                                      has_vector_kernel<vector>::type::value>::type::value>
            ::compute_const(v.v,d);
        }
        else
        {
            vector_dist_op_compute_op<prp,
                                      false,
                                      vector_dist_expression_comp_sel<comp_dev,
                                                                      has_vector_kernel<vector>::type::value>::type::value>
            ::compute_const(v.v,d);
        }

        return v.v;
    }


    template<typename Sys_eqs, typename pmap_type, typename unordered_map_type, typename coeff_type>
    inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff, unsigned int comp) const
    {
            cols[p_map. template getProp<0>(key)*Sys_eqs::nvar + var_id + comp] += coeff;
    }

    inline vector_dist_expression_op<vector_dist_expression<prp,vector>,boost::mpl::int_<1>,VECT_COMP> operator[](int comp)
    {
        int comp_n[1];

        comp_n[0] = comp;

        vector_dist_expression_op<vector_dist_expression<prp,vector>,boost::mpl::int_<1>,VECT_COMP> v_exp(*this,comp_n,var_id);

        return v_exp;
    }
};

template<typename T>
class vector_dist_expression<0,openfpm::vector<aggregate<T>> >
{
    typedef openfpm::vector<aggregate<T>> vector;

    mutable vector v;

public:

    typedef T * iterator;
    typedef const  T * const_iterator;

    typedef typename has_vector_kernel<vector>::type is_ker;

    typedef vector vtype;

    typedef T value_type;

    typedef boost::mpl::bool_<false> is_sort;

    typedef void NN_type;

    static const unsigned int prop = 0;

    int var_id = 0;

    void setVarId(int var_id)
    {
        this->var_id = var_id;
    }

    iterator begin()
    { return &v.template get<0>(0); }

    const_iterator begin() const
    { return &v.template get<0>(0); }

    iterator end()
    { return &v.template get<0>(v.size()-1)+1; }

    const_iterator end() const
    { return &v.template get<0>(v.size()-1)+1; }

    size_t size() const
    { return v.size(); }

    void resize(size_t n)
    {
        // Here

        v.resize(n);
    }

/*  T * begin() {
        return &v.template get<0>(0);
    }

    T * end() {
        return &v.template get<0>(v.size()-1);
    }*/

    // ... [ implement container interface ]
//]
    //const double & operator[]( const size_t n ) const
    //{ return m_v[n]; }

    //double & operator[]( const size_t n )
    //{ return m_v[n]; }



    vector_dist_expression()
    {}

    template<typename exp1, typename exp2, unsigned int op>
    vector_dist_expression(const vector_dist_expression_op<exp1,exp2,op> & v_exp)
    {
        this->operator=(v_exp);
    }

    inline void * getNN() const
    {
        return NULL;
    }

    __device__ __host__ const vector & getVector() const
    {
        return v;
    }

    __device__ __host__ vector & getVector()
    {
        return v;
    }

    inline void init() const
    {}

    __host__ inline auto value(const vect_dist_key_dx & k) const -> decltype(v.template get<0>(k.getKey()))
    {
        return v.template get<0>(k.getKey());
    }


    template<unsigned int prp2, typename vector2> vector & operator=(const vector_dist_expression<prp2,vector2> & v_exp)
    {
        if (v_exp.getVector().isSubset() == true)
        {
                        std::cout << __FILE__ << ":" << __LINE__ << " error on the right hand side of the expression you have to use non-subset properties" << std::endl;
                        return v;
        }

        v.resize(v_exp.getVector().size_local());

        if (has_vector_kernel<vector>::type::value == false)
        {
            vector_dist_op_compute_op<0,false,vector_dist_expression_comp_sel<comp_host,
                                                                          has_vector_kernel<vector>::type::value>::type::value>
            ::compute_expr(v,v_exp);
        }
        else
        {
            vector_dist_op_compute_op<0,false,vector_dist_expression_comp_sel<comp_dev,
                                                                          has_vector_kernel<vector>::type::value>::type::value>
            ::compute_expr(v,v_exp);
        }

        return v;
    }


    template<typename exp1, typename exp2, unsigned int op>
    vector & operator=(const vector_dist_expression_op<exp1,exp2,op> & v_exp)
    {
        if (v_exp.getVector().isSubset() == true)
        {
            std::cout << __FILE__ << ":" << __LINE__ << " error on the right hand side of the expression you have to use non-subset properties" << std::endl;
            return v;
        }

        v.resize(v_exp.getVector().size_local());

        if (has_vector_kernel<vector>::type::value == false)
        {
            vector_dist_op_compute_op<0,
                                      vector_dist_expression_op<exp1,exp2,op>::is_sort::value,
                                      vector_dist_expression_comp_sel<comp_host,
                                                                      has_vector_kernel<vector>::type::value>::type::value>
            ::compute_expr(v,v_exp);
        }
        else
        {
            vector_dist_op_compute_op<0,
                                      vector_dist_expression_op<exp1,exp2,op>::is_sort::value,
                                      vector_dist_expression_comp_sel<comp_dev,
                                                                      has_vector_kernel<vector>::type::value>::type::value>
            ::compute_expr(v,v_exp);
        }

        return v;
    }

    vector & operator=(double d)
    {
        std::cout << __FILE__ << ":" << __LINE__ << " Error: temporal with constants is unsupported" << std::endl;
    }


    template<typename Sys_eqs, typename pmap_type, typename unordered_map_type, typename coeff_type>
    inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff, unsigned int comp) const
    {
        std::cout << __FILE__ << ":" << __LINE__ << " Error: use of temporal is not supported to construct equations";
    }

    inline vector_dist_expression_op<vector_dist_expression<0,vector>,boost::mpl::int_<1>,VECT_COMP> operator[](int comp)
    {
        int comp_n[1];

        comp_n[0] = comp;

        vector_dist_expression_op<vector_dist_expression<0,vector>,boost::mpl::int_<1>,VECT_COMP> v_exp(*this,comp_n,var_id);

        return v_exp;
    }
};


template<typename T> using texp_v = vector_dist_expression<0,openfpm::vector<aggregate<T>> >;

template<typename vector, unsigned int impl>
struct switcher_get_v
{
    typedef vector type;

    static vector & get(vector & v) {return v;};

    template<typename exp_type, typename vector_klist>
    static void register_vector(exp_type & exp_v, vector_klist & v,bool is_sort)
    {
    }
};

template<typename vector>
struct switcher_get_v<vector,comp_dev>
{
    typedef decltype(std::declval<vector>().toKernel()) type;

    static type get(vector & v) {return v.toKernel();};

    template<typename exp_type, typename vector_klist>
    static void register_vector(exp_type & exp_v, vector_klist & v,bool is_sort)
    {
        exp_v.set_vector_dist_ker_list(v.private_get_vector_dist_ker_list(),is_sort);
    }
};

template <typename exp1,int n>
class vector_dist_expression_op<exp1,boost::mpl::int_<n>,VECT_COMP>
{
    exp1 o1;

    int comp[n];

    int var_id = 0;
    void setVarId(int var_id)
    {
        this->var_id = var_id;
    }

    typedef vector_dist_expression_op<exp1,boost::mpl::int_<n>,VECT_COMP> myself;

public:

        typedef std::false_type is_ker;

        typedef boost::mpl::bool_<false> is_sort;

        typedef boost::mpl::bool_<false> NN_type;

    typedef typename exp1::vtype vtype;


    vector_dist_expression_op(const exp1 & o1, int (& comp)[n], int var_id)
    :o1(o1),var_id(var_id)
    {
        for (int i = 0 ; i < n ; i++)
        {this->comp[i] = comp[i];}
    }

    __device__ __host__ const vtype & getVector() const
    {
        return o1.getVector();
    }

    __device__ __host__ vtype & getVector()
    {
        return o1.getVector();
    }

    inline void init() const
    {
        o1.init();
    }

    //typedef typename boost::mpl::at<typename vtype::value_type::type,boost::mpl::int_<exp1::prop>>::type property_act;
    typedef typename pos_or_propL<vtype,exp1::prop>::property_act property_act;

    __host__ __device__ inline auto value(const vect_dist_key_dx & key) const -> decltype(get_vector_dist_expression_op<n,n == rank_gen<property_act>::type::value>::get(o1,vect_dist_key_dx(0),comp))
    {
        return get_vector_dist_expression_op<n,n == rank_gen<property_act>::type::value>::get(o1,key,comp);
    }

    __host__ __device__ inline auto value(const vect_dist_key_dx & key) -> decltype(get_vector_dist_expression_op<n,n == rank_gen<property_act>::type::value>::get(o1,vect_dist_key_dx(0),comp))
    {
        return get_vector_dist_expression_op<n,n == rank_gen<property_act>::type::value>::get(o1,key,comp);
    }

    inline auto get(const vect_dist_key_dx & key) const -> decltype(value(key))
    {
        return this->value(key);
    }

    template<typename Sys_eqs, typename pmap_type, typename unordered_map_type, typename coeff_type>
    inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff, unsigned int comp_) const
    {
#ifdef SE_CLASS1

        if (n != 1)
        {
            std::cout << __FILE__ << ":" << __LINE__ << " Error it only work for tensore of rank 1 ... like vectors " << std::endl;
        }

#endif

        o1.template value_nz<Sys_eqs>(p_map,key,cols,coeff,comp_ + var_id + comp[0]);
    }

    inline vector_dist_expression_op<exp1,boost::mpl::int_<2>,VECT_COMP> operator[](int comp_)
    {
        int comp_n[n+1];

        for (int i = 0 ; i < n ; i++)
        {comp_n[i] = comp[i];}
        comp_n[n] = comp_;

        vector_dist_expression_op<exp1,boost::mpl::int_<2>,VECT_COMP> v_exp(o1,comp_n,var_id);

        return v_exp;
    }

    template<typename T> vtype & operator=(const vector_dist_expression<0,openfpm::vector<aggregate<T>>> & v_exp)
    {
        v_exp.init();

        auto & v = getVector();
/*#ifdef SE_CLASS1
        auto &v2=v_exp.getVector();

        SubsetSelector_impl<std::remove_reference<decltype(v)>::type::is_it_a_subset::value>::check(v2,v);
#endif*/

        auto it = v.getDomainIterator();

        while (it.isNext())
        {
            auto key = it.get();
            auto key_orig = v.getOriginKey(key);

            get_vector_dist_expression_op<n,n == rank_gen<property_act>::type::value>::template assign<exp1::prop>(v_exp,v,key,key_orig,comp);

            ++it;
        }

        return v;
    }

    template<unsigned int prp2> vtype & operator=(const vector_dist_expression<prp2,vtype> & v_exp)
    {
        v_exp.init();

        auto & v = getVector();
#ifdef SE_CLASS1
        auto &v2=v_exp.getVector();

        SubsetSelector_impl<std::remove_reference<decltype(v)>::type::is_it_a_subset::value>::check(v2,v);
#endif
        auto it = v.getDomainIterator();

        while (it.isNext())
        {
            auto key = it.get();
            auto key_orig = v.getOriginKey(key);

            get_vector_dist_expression_op<n,n == rank_gen<property_act>::type::value>::template assign<exp1::prop>(v_exp,v,key,key_orig,comp);

            ++it;
        }

        return v;
    }

    template<typename exp1_, typename exp2_, unsigned int op> vtype & operator=(const vector_dist_expression_op<exp1_,exp2_,op> & v_exp)
    {
        if (v_exp.getVector().isSubset() == true)
        {
            std::cout << __FILE__ << ":" << __LINE__ << " error on the right hand side of the expression you have to use non-subset properties" << std::endl;
            return this->getVector();
        }

        if (has_vector_kernel<vtype>::type::value == false)
        {
            vector_dist_op_compute_op<exp1::prop,false,vector_dist_expression_comp_sel<comp_host,
                                                                          has_vector_kernel<vtype>::type::value>::type::value>
            ::compute_expr_slice(o1.getVector(),v_exp,comp);
        }
        else
        {
            vector_dist_op_compute_op<exp1::prop,false,vector_dist_expression_comp_sel<comp_dev,
                                                                          has_vector_kernel<vtype>::type::value>::type::value>
            ::compute_expr_slice(o1.getVector(),v_exp,comp);
        }

        return this->getVector();
    }

    vtype & operator=(double d)
    {
        auto & v = getVector();

        auto it = v.getDomainIterator();

        while (it.isNext())
        {
            auto key = it.get();

            //pos_or_propL<vtype,exp1::prp>::value(v,key) = d;
            get_vector_dist_expression_op<n,n == rank_gen<property_act>::type::value>::template assign_double<exp1::prop>(d,v,key,comp);


            ++it;
        }

        return v;
    }
};

template <unsigned int prp,unsigned int impl = comp_host, typename vector>
inline vector_dist_expression<prp,typename switcher_get_v<vector,impl>::type > getV(vector & v)
{
    decltype(switcher_get_v<vector,impl>::get(v)) vk = switcher_get_v<vector,impl>::get(v);
    vector_dist_expression<prp,typename switcher_get_v<vector,impl>::type > exp_v(vk);

    switcher_get_v<vector,impl>::register_vector(exp_v,v,false);

    return exp_v;
}


template <unsigned int prp,typename vector>
inline vector_dist_expression<prp, typename switcher_get_v<vector,comp_dev>::type > getV_sort(vector & v)
{
    auto vk = v.toKernel_sorted();
    vector_dist_expression<prp,typename switcher_get_v<vector, comp_dev>::type > exp_v(vk);

    exp_v.set_vector_dist_ker_list(v.private_get_vector_dist_ker_list(),true);

    return exp_v;
}

template<unsigned int prp>
class vector_dist_expression<prp,double>
{
    double d;

public:

    typedef std::false_type is_ker;

    typedef boost::mpl::bool_<false> is_sort;

    typedef void NN_type;

    inline vector_dist_expression(const double & d)
    :d(d)
    {}

    inline void init() const
    {}

    __device__ __host__ inline double value(const vect_dist_key_dx & k) const
    {
        return d;
    }


    template<typename Sys_eqs, typename pmap_type, typename unordered_map_type, typename coeff_type>
    inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff, unsigned int comp) const
    {
        cols[p_map. template getProp<0>(key)*Sys_eqs::nvar + comp] += coeff;
    }
    __host__ __device__ inline double value(const unsigned int & k) const
    {
        return d;
    }
};


template<unsigned int prp>
class vector_dist_expression<prp,float>
{
    float d;

public:

    typedef std::false_type is_ker;


    typedef boost::mpl::bool_<false> is_sort;

    typedef void NN_type;

    inline vector_dist_expression(const float & d)
    :d(d)
    {}

    inline void init() const
    {}

    inline float value(const vect_dist_key_dx & k) const
    {
        return d;
    }

    __device__ __host__ inline float value(const unsigned int & k) const
    {
        return d;
    }
};

/* \brief sum two distributed vector expression
 *
 * \param va vector expression one
 * \param vb vector expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<unsigned int p1, unsigned int p2, typename v1, typename v2>
inline vector_dist_expression_op<vector_dist_expression<p1,v1>,vector_dist_expression<p2,v2>,VECT_SUM>
operator+(const vector_dist_expression<p1,v1> & va, const vector_dist_expression<p2,v2> & vb)
{
    vector_dist_expression_op<vector_dist_expression<p1,v1>,vector_dist_expression<p2,v2>,VECT_SUM> exp_sum(va,vb);

    return exp_sum;
}

/* \brief sum two distributed vector expression
 *
 * \param va vector expression one
 * \param vb vector expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename exp1 , typename exp2, unsigned int op1, unsigned int prp1, typename v1>
inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<prp1,v1>,VECT_SUM>
operator+(const vector_dist_expression_op<exp1,exp2,op1> & va, const vector_dist_expression<prp1,v1> & vb)
{
    vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<prp1,v1>,VECT_SUM> exp_sum(va,vb);

    return exp_sum;
}

/* \brief sum two distributed vector expression
 *
 * \param va vector expression one
 * \param vb vector expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename exp1 , typename exp2, unsigned int op1, unsigned int prp1, typename v1>
inline vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression_op<exp1,exp2,op1>,VECT_SUM>
operator+(const vector_dist_expression<prp1,v1> & va, const vector_dist_expression_op<exp1,exp2,op1> & vb)
{
    vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression_op<exp1,exp2,op1>,VECT_SUM> exp_sum(va,vb);

    return exp_sum;
}

/* \brief sum two distributed vector expression
 *
 * \param va vector expression one
 * \param vb vector expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename exp1 , typename exp2, unsigned int op1, typename exp3 , typename exp4, unsigned int op2>
inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression_op<exp3,exp4,op2>,VECT_SUM>
operator+(const vector_dist_expression_op<exp1,exp2,op1> & va, const vector_dist_expression_op<exp3,exp4,op2> & vb)
{
    vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression_op<exp3,exp4,op2>,VECT_SUM> exp_sum(va,vb);

    return exp_sum;
}

/* \brief sum two distributed vector expression
 *
 * \param va vector expression one
 * \param vb vector expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<unsigned int prp1 , typename v1>
inline vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression<0,double>,VECT_SUM>
operator+(const vector_dist_expression<prp1,v1> & va, double d)
{
    vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression<0,double>,VECT_SUM> exp_sum(va,vector_dist_expression<0,double>(d));

    return exp_sum;
}

/* \brief sum two distributed vector expression
 *
 * \param va vector expression one
 * \param vb vector expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<unsigned int prp1 , typename v1>
inline vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression<0,float>,VECT_SUM>
operator+(const vector_dist_expression<prp1,v1> & va, float d)
{
    vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression<0,float>,VECT_SUM> exp_sum(va,vector_dist_expression<0,float>(d));

    return exp_sum;
}

/* \brief sum two distributed vector expression
 *
 * \param va vector expression one
 * \param vb vector expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<unsigned int prp1 , typename v1>
inline vector_dist_expression_op<vector_dist_expression<0,double>,vector_dist_expression<prp1,v1>,VECT_SUM>
operator+(double d, const vector_dist_expression<prp1,v1> & vb)
{
    vector_dist_expression_op<vector_dist_expression<0,double>,vector_dist_expression<prp1,v1>,VECT_SUM> exp_sum(vector_dist_expression<0,double>(d),vb);

    return exp_sum;
}

/* \brief sum two distributed vector expression
 *
 * \param va vector expression one
 * \param vb vector expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<unsigned int prp1 , typename v1>
inline vector_dist_expression_op<vector_dist_expression<0,float>,vector_dist_expression<prp1,v1>,VECT_SUM>
operator+(float d, const vector_dist_expression<prp1,v1> & vb)
{
    vector_dist_expression_op<vector_dist_expression<0,float>,vector_dist_expression<prp1,v1>,VECT_SUM> exp_sum(vector_dist_expression<0,float>(d),vb);

    return exp_sum;
}

/* \brief sum two distributed vector expression
 *
 * \param va vector expression one
 * \param vb vector expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename exp1 , typename exp2, unsigned int op1>
inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<0,double>,VECT_SUM>
operator+(const vector_dist_expression_op<exp1,exp2,op1> & va, double d)
{
    vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<0,double>,VECT_SUM> exp_sum(va,vector_dist_expression<0,double>(d));

    return exp_sum;
}

/* \brief sum two distributed vector expression
 *
 * \param va vector expression one
 * \param vb vector expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename exp1 , typename exp2, unsigned int op1>
inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<0,float>,VECT_SUM>
operator+(const vector_dist_expression_op<exp1,exp2,op1> & va, float d)
{
    vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<0,float>,VECT_SUM> exp_sum(va,vector_dist_expression<0,float>(d));

    return exp_sum;
}

/* \brief subtract two distributed vector expression
 *
 * \param va vector expression one
 * \param vb vector expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<unsigned int p1, unsigned int p2, typename v1, typename v2>
inline vector_dist_expression_op<vector_dist_expression<p1,v1>,vector_dist_expression<p2,v2>,VECT_SUB>
operator-(const vector_dist_expression<p1,v1> & va, const vector_dist_expression<p2,v2> & vb)
{
    vector_dist_expression_op<vector_dist_expression<p1,v1>,vector_dist_expression<p2,v2>,VECT_SUB> exp_sum(va,vb);

    return exp_sum;
}


/* \brief subtract two distributed vector expression
 *
 * \param va vector expression one
 * \param vb vector expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename exp1, typename exp2, unsigned int op1, unsigned int p2, typename v2>
inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<p2,v2>,VECT_SUB>
operator-(const vector_dist_expression_op<exp1,exp2,op1> & va, const vector_dist_expression<p2,v2> & vb)
{
    vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<p2,v2>,VECT_SUB> exp_sum(va,vb);

    return exp_sum;
}

/* \brief minus of a distributed vector expression operator
 *
 * \param va vector expression one
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename exp1, typename exp2_, unsigned int op1>
inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2_,op1>,void,VECT_SUB_UNI>
operator-(const vector_dist_expression_op<exp1,exp2_,op1> & va)
{
    vector_dist_expression_op<vector_dist_expression_op<exp1,exp2_,op1>,void,VECT_SUB_UNI> exp_sum(va);

    return exp_sum;
}

/* \brief minus of a distributed vector expression
 *
 * \param va vector expression one
 *
 * \return an object that encapsulate the expression
 *
 */
template<unsigned int p1, typename v1>
inline vector_dist_expression_op<vector_dist_expression<p1,v1>,void,VECT_SUB_UNI>
operator-(const vector_dist_expression<p1,v1> & va)
{
    vector_dist_expression_op<vector_dist_expression<p1,v1>,void,VECT_SUB_UNI> exp_sum(va);

    return exp_sum;
}


/* \brief subtract two distributed vector expression
 *
 * \param va vector expression one
 * \param vb vector expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename exp1, typename exp2, unsigned int op1, unsigned int p2, typename v2>
inline vector_dist_expression_op<vector_dist_expression<p2,v2>,vector_dist_expression_op<exp1,exp2,op1>,VECT_SUB>
operator-(const vector_dist_expression<p2,v2> & va, const vector_dist_expression_op<exp1,exp2,op1> & vb)
{
    vector_dist_expression_op<vector_dist_expression<p2,v2>, vector_dist_expression_op<exp1,exp2,op1>,VECT_SUB> exp_sum(va,vb);

    return exp_sum;
}

/* \brief subtract two distributed vector expression
 *
 * \param va vector expression one
 * \param vb vector expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename exp1, typename exp2, unsigned int op1, typename exp3, typename exp4, unsigned int op2>
inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression_op<exp3,exp4,op2>,VECT_SUB>
operator-(const vector_dist_expression_op<exp1,exp2,op1> & va, const vector_dist_expression_op<exp3,exp4,op2> & vb)
{
    vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression_op<exp3,exp4,op2>,VECT_SUB> exp_sum(va,vb);

    return exp_sum;
}

/* \brief subtract two distributed vector expression
 *
 * \param va vector expression one
 * \param vb vector expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<unsigned int prp1, typename v1>
inline vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression<0,double>,VECT_SUB>
operator-(const vector_dist_expression<prp1,v1> & va, double d)
{
    vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression<0,double>,VECT_SUB> exp_sum(va,vector_dist_expression<0,double>(d));

    return exp_sum;
}

/* \brief subtract two distributed vector expression
 *
 * \param va vector expression one
 * \param vb vector expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<unsigned int prp1, typename v1>
inline vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression<0,float>,VECT_SUB>
operator-(const vector_dist_expression<prp1,v1> & va, float d)
{
    vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression<0,float>,VECT_SUB> exp_sum(va,vector_dist_expression<0,float>(d));

    return exp_sum;
}

/* \brief subtract two distributed vector expression
 *
 * \param va vector expression one
 * \param vb vector expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<unsigned int prp1, typename v1>
inline vector_dist_expression_op<vector_dist_expression<0,double>,vector_dist_expression<prp1,v1>,VECT_SUB>
operator-(double d, const vector_dist_expression<prp1,v1> & vb)
{
    vector_dist_expression_op<vector_dist_expression<0,double>,vector_dist_expression<prp1,v1>,VECT_SUB> exp_sum(vector_dist_expression<0,double>(d),vb);

    return exp_sum;
}

/* \brief subtract two distributed vector expression
 *
 * \param va vector expression one
 * \param vb vector expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<unsigned int prp1, typename v1>
inline vector_dist_expression_op<vector_dist_expression<0,float>,vector_dist_expression<prp1,v1>,VECT_SUB>
operator-(float d, const vector_dist_expression<prp1,v1> & vb)
{
    vector_dist_expression_op<vector_dist_expression<0,float>,vector_dist_expression<prp1,v1>,VECT_SUB> exp_sum(vector_dist_expression<0,float>(d),vb);

    return exp_sum;
}

/* \brief Multiply two distributed vector expression
 *
 * \param va vector expression one
 * \param vb vector expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<unsigned int p2, typename v2>
inline vector_dist_expression_op<vector_dist_expression<0,double>,vector_dist_expression<p2,v2>,VECT_MUL>
operator*(double d, const vector_dist_expression<p2,v2> & vb)
{
    vector_dist_expression_op<vector_dist_expression<0,double>,vector_dist_expression<p2,v2>,VECT_MUL> exp_sum(vector_dist_expression<0,double>(d),vb);

    return exp_sum;
}

/* \brief Multiply two distributed vector expression
 *
 * \param va vector expression one
 * \param vb vector expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<unsigned int p2, typename v2>
inline vector_dist_expression_op<vector_dist_expression<0,float>,vector_dist_expression<p2,v2>,VECT_MUL>
operator*(float d, const vector_dist_expression<p2,v2> & vb)
{
    vector_dist_expression_op<vector_dist_expression<0,float>,vector_dist_expression<p2,v2>,VECT_MUL> exp_sum(vector_dist_expression<0,float>(d),vb);

    return exp_sum;
}

/* \brief Multiply two distributed vector expression
 *
 * \param va vector expression one
 * \param vb vector expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<unsigned int p2, typename v2>
inline vector_dist_expression_op<vector_dist_expression<p2,v2>,vector_dist_expression<0,double>,VECT_MUL>
operator*(const vector_dist_expression<p2,v2> & va, double d)
{
    vector_dist_expression_op<vector_dist_expression<p2,v2>,vector_dist_expression<0,double>,VECT_MUL> exp_sum(va,vector_dist_expression<0,double>(d));

    return exp_sum;
}

/* \brief Multiply two distributed vector expression
 *
 * \param va vector expression one
 * \param vb vector expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<unsigned int p2, typename v2>
inline vector_dist_expression_op<vector_dist_expression<p2,v2>,vector_dist_expression<0,float>,VECT_MUL>
operator*(const vector_dist_expression<p2,v2> & va, float d)
{
    vector_dist_expression_op<vector_dist_expression<p2,v2>,vector_dist_expression<0,float>,VECT_MUL> exp_sum(va,vector_dist_expression<0,float>(d));

    return exp_sum;
}

/* \brief Multiply two distributed vector expression
 *
 * \param va vector expression one
 * \param vb vector expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<unsigned int p1, typename v1,unsigned int p2, typename v2>
inline vector_dist_expression_op<vector_dist_expression<p1,v1>,vector_dist_expression<p2,v2>,VECT_MUL>
operator*(const vector_dist_expression<p1,v1> & va, const vector_dist_expression<p2,v2> & vb)
{
    vector_dist_expression_op<vector_dist_expression<p1,v1>,vector_dist_expression<p2,v2>,VECT_MUL> exp_sum(va,vb);

    return exp_sum;
}

/* \brief Multiply two distributed vector expression
 *
 * \param va vector expression one
 * \param vb vector expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<unsigned int p1, typename v1, typename exp1, typename exp2, unsigned int op1>
inline vector_dist_expression_op<vector_dist_expression<p1,v1>,vector_dist_expression_op<exp1,exp2,op1>,VECT_MUL>
operator*(const vector_dist_expression<p1,v1> & va, const vector_dist_expression_op<exp1,exp2,op1> & vb)
{
    vector_dist_expression_op<vector_dist_expression<p1,v1>,vector_dist_expression_op<exp1,exp2,op1>,VECT_MUL> exp_sum(va,vb);

    return exp_sum;
}

/* \brief Multiply two distributed vector expression
 *
 * \param va vector expression one
 * \param vb vector expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<unsigned int p1, typename v1, typename exp1, typename exp2, unsigned int op1>
inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<p1,v1>,VECT_MUL>
operator*(const vector_dist_expression_op<exp1,exp2,op1> & va, const vector_dist_expression<p1,v1> & vb)
{
    vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<p1,v1>,VECT_MUL> exp_sum(va,vb);

    return exp_sum;
}

/* \brief Multiply two distributed vector expression
 *
 * \param va vector expression one
 * \param vb vector expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename exp1, typename exp2, unsigned int op1, typename exp3 , typename exp4, unsigned int op2>
inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression_op<exp3,exp4,op2>,VECT_MUL>
operator*(const vector_dist_expression_op<exp1,exp2,op1> & va, const vector_dist_expression_op<exp3,exp4,op2> & vb)
{
    vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression_op<exp3,exp4,op2>,VECT_MUL> exp_sum(va,vb);

    return exp_sum;
}

/* \brief Multiply a distributed vector expression by a number
 *
 * \param va vector expression
 * \param d number
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename exp1 , typename exp2, unsigned int op1>
inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<0,double>,VECT_MUL>
operator*(const vector_dist_expression_op<exp1,exp2,op1> & va, double d)
{
    vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<0,double>,VECT_MUL> exp_sum(va,vector_dist_expression<0,double>(d));

    return exp_sum;
}

/* \brief Multiply a distributed vector expression by a number
 *
 * \param va vector expression
 * \param d number
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename exp1 , typename exp2, unsigned int op1>
inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<0,float>,VECT_MUL>
operator*(const vector_dist_expression_op<exp1,exp2,op1> & va, float d)
{
    vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<0,float>,VECT_MUL> exp_sum(va,vector_dist_expression<0,float>(d));

    return exp_sum;
}

/* \brief Multiply a distributed vector expression by a number
 *
 * \param d number
 * \param vb vector expression
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename exp1 , typename exp2, unsigned int op1>
inline vector_dist_expression_op<vector_dist_expression<0,double>,vector_dist_expression_op<exp1,exp2,op1>,VECT_MUL>
operator*(double d, const vector_dist_expression_op<exp1,exp2,op1> & vb)
{
    vector_dist_expression_op<vector_dist_expression<0,double>,vector_dist_expression_op<exp1,exp2,op1>,VECT_MUL> exp_sum(vector_dist_expression<0,double>(d),vb);

    return exp_sum;
}

/* \brief Multiply a distributed vector expression by a number
 *
 * \param d number
 * \param vb vector expression
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename exp1 , typename exp2, unsigned int op1>
inline vector_dist_expression_op<vector_dist_expression<0,float>,vector_dist_expression_op<exp1,exp2,op1>,VECT_MUL>
operator*(float d, const vector_dist_expression_op<exp1,exp2,op1> & vb)
{
    vector_dist_expression_op<vector_dist_expression<0,float>,vector_dist_expression_op<exp1,exp2,op1>,VECT_MUL> exp_sum(vector_dist_expression<0,float>(d),vb);

    return exp_sum;
}

/* \brief Divide two distributed vector expression
 *
 * \param va vector expression one
 * \param vb vector expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename exp1, typename exp2, unsigned int op1>
inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<0,double>,VECT_DIV>
operator/(const vector_dist_expression_op<exp1,exp2,op1> & va, double d)
{
    vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<0,double>,VECT_DIV> exp_sum(va,vector_dist_expression<0,double>(d));

    return exp_sum;
}

/* \brief Divide two distributed vector expression
 *
 * \param va vector expression one
 * \param vb vector expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename exp1, typename exp2, unsigned int op1>
inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<0,float>,VECT_DIV>
operator/(const vector_dist_expression_op<exp1,exp2,op1> & va, float d)
{
    vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<0,float>,VECT_DIV> exp_sum(va,vector_dist_expression<0,float>(d));

    return exp_sum;
}

/* \brief Divide two distributed vector expression
 *
 * \param va vector expression one
 * \param vb vector expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename exp1, typename exp2, unsigned int op1>
inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<0,double>,VECT_DIV>
operator/(double d, const vector_dist_expression_op<exp1,exp2,op1> & va)
{
    vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<0,double>,VECT_DIV> exp_sum(vector_dist_expression<0,double>(d),va);

    return exp_sum;
}

/* \brief Divide two distributed vector expression
 *
 * \param va vector expression one
 * \param vb vector expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename exp1, typename exp2, unsigned int op1>
inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<0,float>,VECT_DIV>
operator/(float d, const vector_dist_expression_op<exp1,exp2,op1> & va)
{
    vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<0,float>,VECT_DIV> exp_sum(vector_dist_expression<0,float>(d),va);

    return exp_sum;
}

/* \brief Divide two distributed vector expression
 *
 * \param va vector expression one
 * \param vb vector expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<unsigned int prp1, typename v1>
inline vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression<0,double>,VECT_DIV>
operator/(const vector_dist_expression<prp1,v1> & va, double d)
{
    vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression<0,double>,VECT_DIV> exp_sum(va,vector_dist_expression<0,double>(d));

    return exp_sum;
}

/* \brief Divide two distributed vector expression
 *
 * \param va vector expression one
 * \param vb vector expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<unsigned int prp1, typename v1>
inline vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression<0,float>,VECT_DIV>
operator/(const vector_dist_expression<prp1,v1> & va, float d)
{
    vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression<0,float>,VECT_DIV> exp_sum(va,vector_dist_expression<0,float>(d));

    return exp_sum;
}

/* \brief Divide two distributed vector expression
 *
 * \param va vector expression one
 * \param vb vector expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<unsigned int prp1, typename v1>
inline vector_dist_expression_op<vector_dist_expression<0,double>,vector_dist_expression<prp1,v1>,VECT_DIV>
operator/(double d, const vector_dist_expression<prp1,v1> & va)
{
    vector_dist_expression_op<vector_dist_expression<0,double>,vector_dist_expression<prp1,v1>,VECT_DIV> exp_sum(vector_dist_expression<0,double>(d),va);

    return exp_sum;
}

/* \brief Divide two distributed vector expression
 *
 * \param va vector expression one
 * \param vb vector expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<unsigned int prp1, typename v1>
inline vector_dist_expression_op<vector_dist_expression<0,float>,vector_dist_expression<prp1,v1>,VECT_DIV>
operator/(float d, const vector_dist_expression<prp1,v1> & va)
{
    vector_dist_expression_op<vector_dist_expression<0,float>,vector_dist_expression<prp1,v1>,VECT_DIV> exp_sum(vector_dist_expression<0,float>(d),va);

    return exp_sum;
}

/* \brief Divide two distributed vector expression
 *
 * \param va vector expression one
 * \param vb vector expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<unsigned int prp1, typename v1, unsigned int prp2, typename v2>
inline vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression<prp2,v2>,VECT_DIV>
operator/(const vector_dist_expression<prp1,v1> & va, const vector_dist_expression<prp2,v2> & vb)
{
    vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression<prp2,v2>,VECT_DIV> exp_sum(va,vb);

    return exp_sum;
}

/* \brief Divide two distributed vector expression
 *
 * \param va vector expression one
 * \param vb vector expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<unsigned int prp1, typename v1, typename exp1,typename exp2, unsigned int op1>
inline vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression_op<exp1,exp2,op1>,VECT_DIV>
operator/(const vector_dist_expression<prp1,v1> & va, const vector_dist_expression_op<exp1,exp2,op1> & vb)
{
    vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_dist_expression_op<exp1,exp2,op1>,VECT_DIV> exp_sum(va,vb);

    return exp_sum;
}

/* \brief Divide two distributed vector expression
 *
 * \param va vector expression one
 * \param vb vector expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<unsigned int prp1, typename v1, typename exp1,typename exp2, unsigned int op1>
inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<prp1,v1>,VECT_DIV>
operator/(const vector_dist_expression_op<exp1,exp2,op1> & va, const vector_dist_expression<prp1,v1> & vb)
{
    vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression<prp1,v1>,VECT_DIV> exp_sum(va,vb);

    return exp_sum;
}

/* \brief Divide two distributed vector expression
 *
 * \param va vector expression one
 * \param vb vector expression two
 *
 * \return an object that encapsulate the expression
 *
 */
template<typename exp1,typename exp2, unsigned int op1, typename exp3, typename exp4, unsigned int op2>
inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression_op<exp3,exp4,op2>,VECT_DIV>
operator/(const vector_dist_expression_op<exp1,exp2,op1> & va, const vector_dist_expression_op<exp3,exp4,op2> & vb)
{
    vector_dist_expression_op<vector_dist_expression_op<exp1,exp2,op1>,vector_dist_expression_op<exp3,exp4,op2>,VECT_DIV> exp_sum(va,vb);

    return exp_sum;
}

#include "vector_dist_operators_apply_kernel.hpp"
#include "vector_dist_operators_functions.hpp"
#include "vector_dist_operators_extensions.hpp"
#include "Operators/Vector/vector_dist_operator_assign.hpp"


#endif /* OPENFPM_NUMERICS_SRC_OPERATORS_VECTOR_VECTOR_DIST_OPERATORS_HPP_ */