From 261804394dd2ac70c47204c0e1751667c1a4b32b Mon Sep 17 00:00:00 2001
From: Abhinav <absingh@mpi-cbg.de>
Date: Wed, 1 Mar 2023 02:23:05 +0100
Subject: [PATCH] Odeint GPU Support under construction.

---
 src/OdeIntegrators/OdeIntegrators.hpp         |  29 +
 .../tests/OdeIntegratores_base_tests.cpp      |   4 +-
 src/OdeIntegrators/vector_algebra_ofp_gpu.hpp | 942 ++++++++++++++++++
 .../Vector/vector_dist_operators.hpp          | 404 +++++---
 4 files changed, 1221 insertions(+), 158 deletions(-)
 create mode 100644 src/OdeIntegrators/vector_algebra_ofp_gpu.hpp

diff --git a/src/OdeIntegrators/OdeIntegrators.hpp b/src/OdeIntegrators/OdeIntegrators.hpp
index 63de81be..a485b42f 100644
--- a/src/OdeIntegrators/OdeIntegrators.hpp
+++ b/src/OdeIntegrators/OdeIntegrators.hpp
@@ -63,6 +63,29 @@ struct state_type_1d_ofp{
     }
 };
 
+/*! \brief A 1d Odeint and Openfpm compatible structure.
+ *
+ *  Use the method this.data.get<d>() to refer to property of all the particles in the dimension d.
+ *
+ * d starts with 0.
+ *
+ */
+struct state_type_1d_ofp_gpu{
+    state_type_1d_ofp_gpu(){
+    }
+    typedef size_t size_type;
+    typedef int is_state_vector;
+    aggregate<texp_v_gpu<double>> data;
+
+    size_t size() const
+    { return data.get<0>().size(); }
+
+    void resize(size_t n)
+    {
+        data.get<0>().resize(n);
+    }
+};
+
 /*! \brief A 2d Odeint and Openfpm compatible structure.
  *
  *  Use the method this.data.get<d>() to refer to property of all the particles in the dimension d.
@@ -175,6 +198,12 @@ namespace boost {
             static const bool value = type::value;
             };
 
+            template<>
+            struct is_resizeable<state_type_1d_ofp_gpu> {
+                typedef boost::true_type type;
+                static const bool value = type::value;
+            };
+
             template<>
             struct is_resizeable<state_type_2d_ofp> {
                 typedef boost::true_type type;
diff --git a/src/OdeIntegrators/tests/OdeIntegratores_base_tests.cpp b/src/OdeIntegrators/tests/OdeIntegratores_base_tests.cpp
index f0ee5530..c1455d08 100644
--- a/src/OdeIntegrators/tests/OdeIntegratores_base_tests.cpp
+++ b/src/OdeIntegrators/tests/OdeIntegratores_base_tests.cpp
@@ -15,7 +15,9 @@
 #include "Vector/vector_dist_subset.hpp"
 #include "Decomposition/Distribution/SpaceDistribution.hpp"
 #include "OdeIntegrators/OdeIntegrators.hpp"
+#ifdef HAVE_EIGEN
 #include "DCPSE/DCPSE_op/DCPSE_op.hpp"
+#endif
 #include "OdeIntegrators/vector_algebra_ofp.hpp"
 
 typedef texp_v<double> state_type;
@@ -466,7 +468,6 @@ BOOST_AUTO_TEST_CASE(odeint_base_test3)
 
 
 #ifdef HAVE_EIGEN
-
 BOOST_AUTO_TEST_CASE(dcpse_op_react_diff_test) {
         size_t edgeSemiSize = 5;
         const size_t sz[2] = {2 * edgeSemiSize+1, 2 * edgeSemiSize+1};
@@ -568,5 +569,4 @@ BOOST_AUTO_TEST_CASE(dcpse_op_react_diff_test) {
         }
 }
 #endif
-
 BOOST_AUTO_TEST_SUITE_END()
diff --git a/src/OdeIntegrators/vector_algebra_ofp_gpu.hpp b/src/OdeIntegrators/vector_algebra_ofp_gpu.hpp
new file mode 100644
index 00000000..4be1d2a0
--- /dev/null
+++ b/src/OdeIntegrators/vector_algebra_ofp_gpu.hpp
@@ -0,0 +1,942 @@
+//
+// Created by Abhinav Singh on 18.02.21.
+//
+
+#ifndef OPENFPM_PDATA_VECTOR_ALGEBRA_OFP_HPP
+#define OPENFPM_PDATA_VECTOR_ALGEBRA_OFP_HPP
+
+namespace boost {
+    namespace numeric {
+        namespace odeint {
+
+/*
+ * This class template has to be overload in order to call vector_space_algebra::norm_inf
+ */
+ //           template< class State, class Enabler = void > struct vector_space_norm_inf;
+
+/*
+ * Example: instantiation for sole doubles and complex
+ */
+/*            template<>
+            struct vector_space_norm_inf< double >
+            {
+                typedef double result_type;
+                double operator()( double x ) const
+                {
+                    using std::abs;
+                    return abs(x);
+                }
+            };
+
+            template<>
+            struct vector_space_norm_inf< float >
+            {
+                typedef float result_type;
+                result_type operator()( float x ) const
+                {
+                    using std::abs;
+                    return abs(x);
+                }
+            };
+
+            template< typename T >
+            struct vector_space_norm_inf< std::complex<T> >
+        {
+            typedef T result_type;
+            result_type operator()( std::complex<T> x ) const
+            {
+                using std::abs;
+                return abs( x );
+            }
+        };*/
+
+        template<typename S1,typename S2>
+        struct for_each_prop_resize{
+            S1 &v1;
+            S2 &v2;
+            /*! \brief constructor
+             *
+             *
+             * \param src source encapsulated object
+             * \param dst destination encapsulated object
+             *
+             */
+            inline for_each_prop_resize(S1 &v1,S2 &v2)
+            :v1(v1),v2(v2)
+            {};
+            //! It call the copy function for each property
+            template<typename T>
+            inline void operator()(T& t) const
+            {
+                v1.data.template get<T::value>().getVector().resize(v2.data.template get<T::value>().getVector().size());
+            }
+
+        };
+
+        /* It copy one element of the chunk for each property
+         *
+         */
+        template<typename vector_type,typename index_type,typename op_type>
+        struct for_each_prop1
+        {
+
+            vector_type &v;
+            index_type &p;
+            op_type &op;
+            /*! \brief constructor
+             *
+             *
+             * \param src source encapsulated object
+             * \param dst destination encapsulated object
+             *
+             */
+            inline for_each_prop1(vector_type &v,index_type &p,op_type &op)
+            :v(v),p(p),op(op)
+            {};
+            //! It call the copy function for each property
+            template<typename T>
+            inline void operator()(T& t) const
+            {
+
+                op(v.data.template get<T::value>().getVector().template get<0>(p));
+            }
+        };
+
+        struct vector_space_algebra_ofp_gpu
+        {
+            template< class S1 , class Op >
+            static void for_each1( S1 &s1 , Op op )
+            {
+
+                // ToDo : build checks, that the +-*/ operators are well defined
+                auto it=s1.data.template get<0>().getVector().getDomainIteratorGPU();
+                while(){
+                    auto p=it.get();
+                    //converting to boost vector ids.
+                    for_each_prop1<S1,size_t,Op> cp(s1,p,op);
+                    //creating an iterator on v_ids[0] [1] [2]
+                    boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
+
+                    ++it;
+                }
+            }
+
+
+            template<typename S1,typename S2,typename index_type,typename op_type>
+            struct for_each_prop2
+            {
+
+                S1 &v1;
+                S2 &v2;
+                index_type &p;
+                op_type &op;
+                /*! \brief constructor
+                 *
+                 *
+                 * \param src source encapsulated object
+                 * \param dst destination encapsulated object
+                 *
+                 */
+                inline for_each_prop2(S1 &v1,S2 &v2,index_type &p,op_type &op)
+                :v1(v1),v2(v2),p(p),op(op)
+                {};
+                //! It call the copy function for each property
+                template<typename T>
+                inline void operator()(T& t) const
+                {
+                    op(v1.data.template get<T::value>().getVector().template get<0>(p),v2.data.template get<T::value>().getVector().template get<0>(p));
+                }
+            };
+            template< class S1 , class S2 , class Op >
+            static void for_each2( S1 &s1 , S2 &s2 , Op op )
+            {
+                for_each_prop_resize<S1,S2> the_resize(s1,s2);
+                boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
+                //s1.data.template get<0>().getVector().resize(s2.data.template get<0>().getVector().size());
+                // ToDo : build checks, that the +-*/ operators are well defined
+                auto it=s1.data.template get<0>().getVector().getIterator();
+                while(it.isNext()){
+                    auto p=it.get();
+                    //converting to boost vector ids.
+                    for_each_prop2<S1,S2,size_t,Op> cp(s1,s2,p,op);
+                    //creating an iterator on v_ids[0] [1] [2]
+                    boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
+
+                    ++it;
+                }
+            }
+
+
+
+            template<typename S1,typename S2,typename S3,typename index_type,typename op_type>
+            struct for_each_prop3
+            {
+
+                S1 &v1;
+                S2 &v2;
+                S3 &v3;
+                index_type &p;
+                op_type &op;
+                /*! \brief constructor
+                 *
+                 *
+                 * \param src source encapsulated object
+                 * \param dst destination encapsulated object
+                 *
+                 */
+                vector_algebra_ofp.hpp                inline for_each_prop3(S1 &v1,S2 &v2,S3 &v3,index_type &p,op_type &op)
+                :v1(v1),v2(v2),v3(v3),p(p),op(op)
+                {};
+                //! It call the copy function for each property
+                template<typename T>
+                inline void operator()(T& t) const
+                {
+                    //std::cout<<v1.data.template get<T::value>().getVector().size()<<":"<<v2.data.template get<T::value>().getVector().size()<<":"<<v3.data.template get<T::value>().getVector().size()<<std::endl;
+                    op(v1.data.template get<T::value>().getVector().template get<0>(p),v2.data.template get<T::value>().getVector().template get<0>(p),v3.data.template get<T::value>().getVector().template get<0>(p));
+                }
+            };
+
+
+            template< class S1 , class S2 , class S3 , class Op >
+            static void for_each3( S1 &s1 , S2 &s2 , S3 &s3 , Op op )
+            {
+
+//
+                for_each_prop_resize<S1,S2> the_resize(s1,s2);
+                boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(the_resize);
+                // ToDo : build checks, that the +-*/ operators are well defined
+                auto it=s1.data.template get<0>().getVector().getIterator();
+                while(it.isNext()){
+                    auto p=it.get();
+                    //converting to boost vector ids.
+                    for_each_prop3<S1,S2,S3,size_t,Op> cp(s1,s2,s3,p,op);
+                    //creating an iterator on v_ids[0] [1] [2]
+                    boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
+
+                    ++it;
+                }
+            }
+
+
+            template<typename S1,typename S2,typename S3,typename S4,typename index_type,typename op_type>
+            struct for_each_prop4
+            {
+
+                S1 &v1;
+                S2 &v2;
+                S3 &v3;
+                S4 &v4;
+
+                index_type &p;
+                op_type &op;
+                /*! \brief constructor
+                 *
+                 *
+                 * \param src source encapsulated object
+                 * \param dst destination encapsulated object
+                 *
+                 */
+                inline for_each_prop4(S1 &v1,S2 &v2,S3 &v3,S4 &v4,index_type &p,op_type &op)
+                :v1(v1),v2(v2),v3(v3),v4(v4),p(p),op(op)
+                {};
+                //! It call the copy function for each property
+                template<typename T>
+                inline void operator()(T& t) const
+                {
+                    op(v1.data.template get<T::value>().getVector().template get<0>(p),v2.data.template get<T::value>().getVector().template get<0>(p),v3.data.template get<T::value>().getVector().template get<0>(p),v4.data.template get<T::value>().getVector().template get<0>(p));
+                }
+            };
+
+            template< class S1 , class S2 , class S3 , class S4 , class Op >
+            static void for_each4( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4 , Op op )
+            {
+                for_each_prop_resize<S1,S2> the_resize(s1,s2);
+                boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
+                // ToDo : build checks, that the +-*/ operators are well defined
+                auto it=s1.data.template get<0>().getVector().getIterator();
+                while(it.isNext()){
+                    auto p=it.get();
+                    //converting to boost vector ids.
+                    for_each_prop4<S1,S2,S3,S4,size_t,Op> cp(s1,s2,s3,s4,p,op);
+                    //creating an iterator on v_ids[0] [1] [2]
+                    boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
+
+                    ++it;
+                }
+            }
+
+
+            template<typename S1,typename S2,typename S3,typename S4,typename S5,typename index_type,typename op_type>
+            struct for_each_prop5
+            {
+
+                S1 &v1;
+                S2 &v2;
+                S3 &v3;
+                S4 &v4;
+                S5 &v5;
+
+                index_type &p;
+                op_type &op;
+                /*! \brief constructor
+                 *
+                 *
+                 * \param src source encapsulated object
+                 * \param dst destination encapsulated object
+                 *
+                 */
+                inline for_each_prop5(S1 &v1,S2 &v2,S3 &v3,S4 &v4,S5 &v5,index_type &p,op_type &op)
+                        :v1(v1),v2(v2),v3(v3),v4(v4),v5(v5),p(p),op(op)
+                {};
+                //! It call the copy function for each property
+                template<typename T>
+                inline void operator()(T& t) const
+                {
+                    op(v1.data.template get<T::value>().getVector().template get<0>(p),v2.data.template get<T::value>().getVector().template get<0>(p),v3.data.template get<T::value>().getVector().template get<0>(p),v4.data.template get<T::value>().getVector().template get<0>(p),v5.data.template get<T::value>().getVector().template get<0>(p));
+                }
+            };
+
+            template< class S1 , class S2 , class S3 , class S4,class S5 , class Op >
+            static void for_each5( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4,S5 &s5 , Op op )
+            {
+                for_each_prop_resize<S1,S2> the_resize(s1,s2);
+                boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
+                // ToDo : build checks, that the +-*/ operators are well defined
+                auto it=s1.data.template get<0>().getVector().getIterator();
+                while(it.isNext()){
+                    auto p=it.get();
+                    //converting to boost vector ids.
+                    for_each_prop5<S1,S2,S3,S4,S5,size_t,Op> cp(s1,s2,s3,s4,s5,p,op);
+                    //creating an iterator on v_ids[0] [1] [2]
+                    boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
+
+                    ++it;
+                }
+            }
+
+
+            template<typename S1,typename S2,typename S3,typename S4,typename S5,typename S6,typename index_type,typename op_type>
+            struct for_each_prop6
+            {
+
+                S1 &v1;
+                S2 &v2;
+                S3 &v3;
+                S4 &v4;
+                S5 &v5;
+                S6 &v6;
+
+
+                index_type &p;
+                op_type &op;
+                /*! \brief constructor
+                 *
+                 *
+                 * \param src source encapsulated object
+                 * \param dst destination encapsulated object
+                 *
+                 */
+                inline for_each_prop6(S1 &v1,S2 &v2,S3 &v3,S4 &v4,S5 &v5,S6 &v6,index_type &p,op_type &op)
+                        :v1(v1),v2(v2),v3(v3),v4(v4),v5(v5),v6(v6),p(p),op(op)
+                {};
+                //! It call the copy function for each property
+                template<typename T>
+                inline void operator()(T& t) const
+                {
+                    op(v1.data.template get<T::value>().getVector().template get<0>(p),v2.data.template get<T::value>().getVector().template get<0>(p),v3.data.template get<T::value>().getVector().template get<0>(p),v4.data.template get<T::value>().getVector().template get<0>(p),v5.data.template get<T::value>().getVector().template get<0>(p),v6.data.template get<T::value>().getVector().template get<0>(p));
+                }
+            };
+            template< class S1 , class S2 , class S3 , class S4,class S5,class S6 , class Op >
+            static void for_each6( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4,S5 &s5,S6 &s6 , Op op )
+            {
+                for_each_prop_resize<S1,S2> the_resize(s1,s2);
+                boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
+                // ToDo : build checks, that the +-*/ operators are well defined
+                auto it=s1.data.template get<0>().getVector().getIterator();
+                while(it.isNext()){
+                    auto p=it.get();
+                    //converting to boost vector ids.
+                    for_each_prop6<S1,S2,S3,S4,S5,S6,size_t,Op> cp(s1,s2,s3,s4,s5,s6,p,op);
+                    //creating an iterator on v_ids[0] [1] [2]
+                    boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
+
+                    ++it;
+                }
+            }
+
+
+
+            template<typename S1,typename S2,typename S3,typename S4,typename S5,typename S6,typename S7,typename index_type,typename op_type>
+            struct for_each_prop7
+            {
+
+                S1 &v1;
+                S2 &v2;
+                S3 &v3;
+                S4 &v4;
+                S5 &v5;
+                S6 &v6;
+                S7 &v7;
+
+
+                index_type &p;
+                op_type &op;
+                /*! \brief constructor
+                 *
+                 *
+                 * \param src source encapsulated object
+                 * \param dst destination encapsulated object
+                 *
+                 */
+                inline for_each_prop7(S1 &v1,S2 &v2,S3 &v3,S4 &v4,S5 &v5,S6 &v6,S7 &v7,index_type &p,op_type &op)
+                        :v1(v1),v2(v2),v3(v3),v4(v4),v5(v5),v6(v6),v7(v7),p(p),op(op)
+                {};
+                //! It call the copy function for each property
+                template<typename T>
+                inline void operator()(T& t) const
+                {
+                    op(v1.data.template get<T::value>().getVector().template get<0>(p),v2.data.template get<T::value>().getVector().template get<0>(p),v3.data.template get<T::value>().getVector().template get<0>(p),v4.data.template get<T::value>().getVector().template get<0>(p),v5.data.template get<T::value>().getVector().template get<0>(p),v6.data.template get<T::value>().getVector().template get<0>(p),v7.data.template get<T::value>().getVector().template get<0>(p));
+                }
+            };
+
+
+            template< class S1 , class S2 , class S3 , class S4,class S5,class S6 ,class S7, class Op >
+            static void for_each7( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4,S5 &s5,S6 &s6,S7 &s7 , Op op )
+            {
+                for_each_prop_resize<S1,S2> the_resize(s1,s2);
+                boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
+                // ToDo : build checks, that the +-*/ operators are well defined
+                auto it=s1.data.template get<0>().getVector().getIterator();
+                while(it.isNext()){
+                    auto p=it.get();
+                    //converting to boost vector ids.
+                    for_each_prop7<S1,S2,S3,S4,S5,S6,S7,size_t,Op> cp(s1,s2,s3,s4,s5,s6,s7,p,op);
+                    //creating an iterator on v_ids[0] [1] [2]
+                    boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
+
+                    ++it;
+                }
+            }
+
+            template<typename S1,typename S2,typename S3,typename S4,typename S5,typename S6,typename S7,typename S8,typename index_type,typename op_type>
+            struct for_each_prop8
+            {
+
+                S1 &v1;
+                S2 &v2;
+                S3 &v3;
+                S4 &v4;
+                S5 &v5;
+                S6 &v6;
+                S7 &v7;
+                S8 &v8;
+
+
+                index_type &p;
+                op_type &op;
+                /*! \brief constructor
+                 *
+                 *
+                 * \param src source encapsulated object
+                 * \param dst destination encapsulated object
+                 *
+                 */
+                inline for_each_prop8(S1 &v1,S2 &v2,S3 &v3,S4 &v4,S5 &v5,S6 &v6,S7 &v7,S8 &v8,index_type &p,op_type &op)
+                        :v1(v1),v2(v2),v3(v3),v4(v4),v5(v5),v6(v6),v7(v7),v8(v8),p(p),op(op)
+                {};
+                //! It call the copy function for each property
+                template<typename T>
+                inline void operator()(T& t) const
+                {
+                    op(v1.data.template get<T::value>().getVector().template get<0>(p),v2.data.template get<T::value>().getVector().template get<0>(p),v3.data.template get<T::value>().getVector().template get<0>(p),v4.data.template get<T::value>().getVector().template get<0>(p),v5.data.template get<T::value>().getVector().template get<0>(p),v6.data.template get<T::value>().getVector().template get<0>(p),v7.data.template get<T::value>().getVector().template get<0>(p),v8.data.template get<T::value>().getVector().template get<0>(p));
+                }
+            };
+
+
+            template< class S1 , class S2 , class S3 , class S4,class S5,class S6 ,class S7,class S8, class Op >
+            static void for_each8( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4,S5 &s5,S6 &s6,S7 &s7,S8 &s8 , Op op )
+            {
+                for_each_prop_resize<S1,S2> the_resize(s1,s2);
+                boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
+                // ToDo : build checks, that the +-*/ operators are well defined
+                auto it=s1.data.template get<0>().getVector().getIterator();
+                while(it.isNext()){
+                    auto p=it.get();
+                    //converting to boost vector ids.
+                    for_each_prop8<S1,S2,S3,S4,S5,S6,S7,S8,size_t,Op> cp(s1,s2,s3,s4,s5,s6,s7,s8,p,op);
+                    //creating an iterator on v_ids[0] [1] [2]
+                    boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
+
+                    ++it;
+                }
+            }
+
+            template<typename S1,typename S2,typename S3,typename S4,typename S5,typename S6,typename S7,typename S8, typename S9,typename index_type,typename op_type>
+            struct for_each_prop9
+            {
+
+                S1 &v1;
+                S2 &v2;
+                S3 &v3;
+                S4 &v4;
+                S5 &v5;
+                S6 &v6;
+                S7 &v7;
+                S8 &v8;
+                S9 &v9;
+
+
+                index_type &p;
+                op_type &op;
+                /*! \brief constructor
+                 *
+                 *
+                 * \param src source encapsulated object
+                 * \param dst destination encapsulated object
+                 *
+                 */
+                inline for_each_prop9(S1 &v1,S2 &v2,S3 &v3,S4 &v4,S5 &v5,S6 &v6,S7 &v7,S8 &v8,S9 &v9,index_type &p,op_type &op)
+                        :v1(v1),v2(v2),v3(v3),v4(v4),v5(v5),v6(v6),v7(v7),v8(v8),v9(v9),p(p),op(op)
+                {};
+                //! It call the copy function for each property
+                template<typename T>
+                inline void operator()(T& t) const
+                {
+                    op(v1.data.template get<T::value>().getVector().template get<0>(p),v2.data.template get<T::value>().getVector().template get<0>(p),v3.data.template get<T::value>().getVector().template get<0>(p),v4.data.template get<T::value>().getVector().template get<0>(p),v5.data.template get<T::value>().getVector().template get<0>(p),v6.data.template get<T::value>().getVector().template get<0>(p),v7.data.template get<T::value>().getVector().template get<0>(p),v8.data.template get<T::value>().getVector().template get<0>(p),v9.data.template get<T::value>().getVector().template get<0>(p));
+                }
+            };
+
+
+            template< class S1 , class S2 , class S3 , class S4,class S5,class S6 ,class S7,class S8, class S9, class Op >
+            static void for_each9( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4,S5 &s5,S6 &s6,S7 &s7,S8 &s8, S9 &s9 , Op op )
+            {
+                for_each_prop_resize<S1,S2> the_resize(s1,s2);
+                boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
+                // ToDo : build checks, that the +-*/ operators are well defined
+                auto it=s1.data.template get<0>().getVector().getIterator();
+                while(it.isNext()){
+                    auto p=it.get();
+                    //converting to boost vector ids.
+                    for_each_prop9<S1,S2,S3,S4,S5,S6,S7,S8,S9,size_t,Op> cp(s1,s2,s3,s4,s5,s6,s7,s8,s9,p,op);
+                    //creating an iterator on v_ids[0] [1] [2]
+                    boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
+
+                    ++it;
+                }
+            }
+
+            template<typename S1,typename S2,typename S3,typename S4,typename S5,typename S6,typename S7,typename S8, typename S9, typename S10,typename index_type,typename op_type>
+            struct for_each_prop10
+            {
+
+                S1 &v1;
+                S2 &v2;
+                S3 &v3;
+                S4 &v4;
+                S5 &v5;
+                S6 &v6;
+                S7 &v7;
+                S8 &v8;
+                S9 &v9;
+                S10 &v10;
+
+
+                index_type &p;
+                op_type &op;
+                /*! \brief constructor
+                 *
+                 *
+                 * \param src source encapsulated object
+                 * \param dst destination encapsulated object
+                 *
+                 */
+                inline for_each_prop10(S1 &v1,S2 &v2,S3 &v3,S4 &v4,S5 &v5,S6 &v6,S7 &v7,S8 &v8,S9 &v9,S10 &v10,index_type &p,op_type &op)
+                        :v1(v1),v2(v2),v3(v3),v4(v4),v5(v5),v6(v6),v7(v7),v8(v8),v9(v9),v10(v10),p(p),op(op)
+                {};
+                //! It call the copy function for each property
+                template<typename T>
+                inline void operator()(T& t) const
+                {
+                    op(v1.data.template get<T::value>().getVector().template get<0>(p),v2.data.template get<T::value>().getVector().template get<0>(p),v3.data.template get<T::value>().getVector().template get<0>(p),v4.data.template get<T::value>().getVector().template get<0>(p),v5.data.template get<T::value>().getVector().template get<0>(p),v6.data.template get<T::value>().getVector().template get<0>(p),v7.data.template get<T::value>().getVector().template get<0>(p),v8.data.template get<T::value>().getVector().template get<0>(p),v9.data.template get<T::value>().getVector().template get<0>(p),v10.data.template get<T::value>().getVector().template get<0>(p));
+                }
+            };
+
+
+            template< class S1 , class S2 , class S3 , class S4,class S5,class S6 ,class S7,class S8, class S9, class S10, class Op >
+            static void for_each10( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4,S5 &s5,S6 &s6,S7 &s7,S8 &s8, S9 &s9 , S10 &s10, Op op )
+            {
+                for_each_prop_resize<S1,S2> the_resize(s1,s2);
+                boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
+                // ToDo : build checks, that the +-*/ operators are well defined
+                auto it=s1.data.template get<0>().getVector().getIterator();
+                while(it.isNext()){
+                    auto p=it.get();
+                    //converting to boost vector ids.
+                    for_each_prop10<S1,S2,S3,S4,S5,S6,S7,S8,S9,S10,size_t,Op> cp(s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,p,op);
+                    //creating an iterator on v_ids[0] [1] [2]
+                    boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
+
+                    ++it;
+                }
+            }
+
+            template<typename S1,typename S2,typename S3,typename S4,typename S5,typename S6,typename S7,typename S8, typename S9, typename S10, typename S11,typename index_type,typename op_type>
+            struct for_each_prop11
+            {
+
+                S1 &v1;
+                S2 &v2;
+                S3 &v3;
+                S4 &v4;
+                S5 &v5;
+                S6 &v6;
+                S7 &v7;
+                S8 &v8;
+                S9 &v9;
+                S10 &v10;
+                S11 &v11;
+
+
+                index_type &p;
+                op_type &op;
+                /*! \brief constructor
+                 *
+                 *
+                 * \param src source encapsulated object
+                 * \param dst destination encapsulated object
+                 *
+                 */
+                inline for_each_prop11(S1 &v1,S2 &v2,S3 &v3,S4 &v4,S5 &v5,S6 &v6,S7 &v7,S8 &v8,S9 &v9,S10 &v10, S11 &v11,index_type &p,op_type &op)
+                        :v1(v1),v2(v2),v3(v3),v4(v4),v5(v5),v6(v6),v7(v7),v8(v8),v9(v9),v10(v10),v11(v11),p(p),op(op)
+                {};
+                //! It call the copy function for each property
+                template<typename T>
+                inline void operator()(T& t) const
+                {
+                    op(v1.data.template get<T::value>().getVector().template get<0>(p),v2.data.template get<T::value>().getVector().template get<0>(p),v3.data.template get<T::value>().getVector().template get<0>(p),v4.data.template get<T::value>().getVector().template get<0>(p),v5.data.template get<T::value>().getVector().template get<0>(p),v6.data.template get<T::value>().getVector().template get<0>(p),v7.data.template get<T::value>().getVector().template get<0>(p),v8.data.template get<T::value>().getVector().template get<0>(p),v9.data.template get<T::value>().getVector().template get<0>(p),v10.data.template get<T::value>().getVector().template get<0>(p),v11.data.template get<T::value>().getVector().template get<0>(p));
+                }
+            };
+
+
+            template< class S1 , class S2 , class S3 , class S4,class S5,class S6 ,class S7,class S8, class S9, class S10, class S11, class Op >
+            static void for_each11( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4,S5 &s5,S6 &s6,S7 &s7,S8 &s8, S9 &s9 , S10 &s10,S11 &s11, Op op )
+            {
+                for_each_prop_resize<S1,S2> the_resize(s1,s2);
+                boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
+                // ToDo : build checks, that the +-*/ operators are well defined
+                auto it=s1.data.template get<0>().getVector().getIterator();
+                while(it.isNext()){
+                    auto p=it.get();
+                    //converting to boost vector ids.
+                    for_each_prop11<S1,S2,S3,S4,S5,S6,S7,S8,S9,S10,S11,size_t,Op> cp(s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,s11,p,op);
+                    //creating an iterator on v_ids[0] [1] [2]
+                    boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
+
+                    ++it;
+                }
+            }
+
+            template<typename S1,typename S2,typename S3,typename S4,typename S5,typename S6,typename S7,typename S8, typename S9, typename S10, typename S11, typename S12,typename index_type,typename op_type>
+            struct for_each_prop12
+            {
+
+                S1 &v1;
+                S2 &v2;
+                S3 &v3;
+                S4 &v4;
+                S5 &v5;
+                S6 &v6;
+                S7 &v7;
+                S8 &v8;
+                S9 &v9;
+                S10 &v10;
+                S11 &v11;
+                S12 &v12;
+
+
+                index_type &p;
+                op_type &op;
+                /*! \brief constructor
+                 *
+                 *
+                 * \param src source encapsulated object
+                 * \param dst destination encapsulated object
+                 *
+                 */
+                inline for_each_prop12(S1 &v1,S2 &v2,S3 &v3,S4 &v4,S5 &v5,S6 &v6,S7 &v7,S8 &v8,S9 &v9,S10 &v10, S11 &v11,S12 &v12,index_type &p,op_type &op)
+                        :v1(v1),v2(v2),v3(v3),v4(v4),v5(v5),v6(v6),v7(v7),v8(v8),v9(v9),v10(v10),v11(v11), v12(v12),p(p),op(op)
+                {};
+                //! It call the copy function for each property
+                template<typename T>
+                inline void operator()(T& t) const
+                {
+                    op(v1.data.template get<T::value>().getVector().template get<0>(p),v2.data.template get<T::value>().getVector().template get<0>(p),v3.data.template get<T::value>().getVector().template get<0>(p),v4.data.template get<T::value>().getVector().template get<0>(p),v5.data.template get<T::value>().getVector().template get<0>(p),v6.data.template get<T::value>().getVector().template get<0>(p),v7.data.template get<T::value>().getVector().template get<0>(p),v8.data.template get<T::value>().getVector().template get<0>(p),v9.data.template get<T::value>().getVector().template get<0>(p),v10.data.template get<T::value>().getVector().template get<0>(p),v11.data.template get<T::value>().getVector().template get<0>(p),v12.data.template get<T::value>().getVector().template get<0>(p));
+                }
+            };
+
+
+            template< class S1 , class S2 , class S3 , class S4,class S5,class S6 ,class S7,class S8, class S9, class S10, class S11, class S12, class Op >
+            static void for_each12( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4,S5 &s5,S6 &s6,S7 &s7,S8 &s8, S9 &s9 , S10 &s10,S11 &s11,S12 &s12, Op op )
+            {
+                for_each_prop_resize<S1,S2> the_resize(s1,s2);
+                boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
+                // ToDo : build checks, that the +-*/ operators are well defined
+                auto it=s1.data.template get<0>().getVector().getIterator();
+                while(it.isNext()){
+                    auto p=it.get();
+                    //converting to boost vector ids.
+                    for_each_prop12<S1,S2,S3,S4,S5,S6,S7,S8,S9,S10,S11,S12,size_t,Op> cp(s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,s11,s12,p,op);
+                    //creating an iterator on v_ids[0] [1] [2]
+                    boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
+
+                    ++it;
+                }
+            }
+
+            template<typename S1,typename S2,typename S3,typename S4,typename S5,typename S6,typename S7,typename S8, typename S9, typename S10, typename S11, typename S12, typename S13,typename index_type,typename op_type>
+            struct for_each_prop13
+            {
+
+                S1 &v1;
+                S2 &v2;
+                S3 &v3;
+                S4 &v4;
+                S5 &v5;
+                S6 &v6;
+                S7 &v7;
+                S8 &v8;
+                S9 &v9;
+                S10 &v10;
+                S11 &v11;
+                S12 &v12;
+                S13 &v13;
+
+
+                index_type &p;
+                op_type &op;
+                /*! \brief constructor
+                 *
+                 *
+                 * \param src source encapsulated object
+                 * \param dst destination encapsulated object
+                 *
+                 */
+                inline for_each_prop13(S1 &v1,S2 &v2,S3 &v3,S4 &v4,S5 &v5,S6 &v6,S7 &v7,S8 &v8,S9 &v9,S10 &v10, S11 &v11,S12 &v12,S13 &v13,index_type &p,op_type &op)
+                        :v1(v1),v2(v2),v3(v3),v4(v4),v5(v5),v6(v6),v7(v7),v8(v8),v9(v9),v10(v10),v11(v11), v12(v12),v13(v13),p(p),op(op)
+                {};
+                //! It call the copy function for each property
+                template<typename T>
+                inline void operator()(T& t) const
+                {
+                    op(v1.data.template get<T::value>().getVector().template get<0>(p),v2.data.template get<T::value>().getVector().template get<0>(p),v3.data.template get<T::value>().getVector().template get<0>(p),v4.data.template get<T::value>().getVector().template get<0>(p),v5.data.template get<T::value>().getVector().template get<0>(p),v6.data.template get<T::value>().getVector().template get<0>(p),v7.data.template get<T::value>().getVector().template get<0>(p),v8.data.template get<T::value>().getVector().template get<0>(p),v9.data.template get<T::value>().getVector().template get<0>(p),v10.data.template get<T::value>().getVector().template get<0>(p),v11.data.template get<T::value>().getVector().template get<0>(p),v12.data.template get<T::value>().getVector().template get<0>(p),v13.data.template get<T::value>().getVector().template get<0>(p));
+                }
+            };
+
+
+            template< class S1 , class S2 , class S3 , class S4,class S5,class S6 ,class S7,class S8, class S9, class S10, class S11, class S12, class S13, class Op >
+            static void for_each13( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4,S5 &s5,S6 &s6,S7 &s7,S8 &s8, S9 &s9 , S10 &s10,S11 &s11,S12 &s12,S13 &s13, Op op )
+            {
+                for_each_prop_resize<S1,S2> the_resize(s1,s2);
+                boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
+                // ToDo : build checks, that the +-*/ operators are well defined
+                auto it=s1.data.template get<0>().getVector().getIterator();
+                while(it.isNext()){
+                    auto p=it.get();
+                    //converting to boost vector ids.
+                    for_each_prop13<S1,S2,S3,S4,S5,S6,S7,S8,S9,S10,S11,S12,S13,size_t,Op> cp(s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,s11,s12,s13,p,op);
+                    //creating an iterator on v_ids[0] [1] [2]
+                    boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
+
+                    ++it;
+                }
+            }
+
+            template<typename S1,typename S2,typename S3,typename S4,typename S5,typename S6,typename S7,typename S8, typename S9, typename S10, typename S11, typename S12, typename S13, typename S14,typename index_type,typename op_type>
+            struct for_each_prop14
+            {
+
+                S1 &v1;
+                S2 &v2;
+                S3 &v3;
+                S4 &v4;
+                S5 &v5;
+                S6 &v6;
+                S7 &v7;
+                S8 &v8;
+                S9 &v9;
+                S10 &v10;
+                S11 &v11;
+                S12 &v12;
+                S13 &v13;
+                S14 &v14;
+
+
+                index_type &p;
+                op_type &op;
+                /*! \brief constructor
+                 *
+                 *
+                 * \param src source encapsulated object
+                 * \param dst destination encapsulated object
+                 *
+                 */
+                inline for_each_prop14(S1 &v1,S2 &v2,S3 &v3,S4 &v4,S5 &v5,S6 &v6,S7 &v7,S8 &v8,S9 &v9,S10 &v10, S11 &v11,S12 &v12,S13 &v13,S14 &v14,index_type &p,op_type &op)
+                        :v1(v1),v2(v2),v3(v3),v4(v4),v5(v5),v6(v6),v7(v7),v8(v8),v9(v9),v10(v10),v11(v11), v12(v12),v13(v13),v14(v14),p(p),op(op)
+                {};
+                //! It call the copy function for each property
+                template<typename T>
+                inline void operator()(T& t) const
+                {
+                    op(v1.data.template get<T::value>().getVector().template get<0>(p),v2.data.template get<T::value>().getVector().template get<0>(p),v3.data.template get<T::value>().getVector().template get<0>(p),v4.data.template get<T::value>().getVector().template get<0>(p),v5.data.template get<T::value>().getVector().template get<0>(p),v6.data.template get<T::value>().getVector().template get<0>(p),v7.data.template get<T::value>().getVector().template get<0>(p),v8.data.template get<T::value>().getVector().template get<0>(p),v9.data.template get<T::value>().getVector().template get<0>(p),v10.data.template get<T::value>().getVector().template get<0>(p),v11.data.template get<T::value>().getVector().template get<0>(p),v12.data.template get<T::value>().getVector().template get<0>(p),v13.data.template get<T::value>().getVector().template get<0>(p),v14.data.template get<T::value>().getVector().template get<0>(p));
+                }
+            };
+
+
+            template< class S1 , class S2 , class S3 , class S4,class S5,class S6 ,class S7,class S8, class S9, class S10, class S11, class S12, class S13, class S14, class Op >
+            static void for_each14( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4,S5 &s5,S6 &s6,S7 &s7,S8 &s8, S9 &s9, S10 &s10,S11 &s11,S12 &s12,S13 &s13,S14 &s14, Op op )
+            {
+                for_each_prop_resize<S1,S2> the_resize(s1,s2);
+                boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
+                // ToDo : build checks, that the +-*/ operators are well defined
+                auto it=s1.data.template get<0>().getVector().getIterator();
+                while(it.isNext()){
+                    auto p=it.get();
+                    //converting to boost vector ids.
+                    for_each_prop14<S1,S2,S3,S4,S5,S6,S7,S8,S9,S10,S11,S12,S13,S14,size_t,Op> cp(s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,s11,s12,s13,s14,p,op);
+                    //creating an iterator on v_ids[0] [1] [2]
+                    boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
+
+                    ++it;
+                }
+            }
+
+            template<typename S1,typename S2,typename S3,typename S4,typename S5,typename S6,typename S7,typename S8, typename S9, typename S10, typename S11, typename S12, typename S13, typename S14, typename S15,typename index_type,typename op_type>
+            struct for_each_prop15
+            {
+
+                S1 &v1;
+                S2 &v2;
+                S3 &v3;
+                S4 &v4;
+                S5 &v5;
+                S6 &v6;
+                S7 &v7;
+                S8 &v8;
+                S9 &v9;
+                S10 &v10;
+                S11 &v11;
+                S12 &v12;
+                S13 &v13;
+                S14 &v14;
+                S15 &v15;
+
+
+                index_type &p;
+                op_type &op;
+                /*! \brief constructor
+                 *
+                 *
+                 * \param src source encapsulated object
+                 * \param dst destination encapsulated object
+                 *
+                 */
+                inline for_each_prop15(S1 &v1,S2 &v2,S3 &v3,S4 &v4,S5 &v5,S6 &v6,S7 &v7,S8 &v8,S9 &v9,S10 &v10, S11 &v11,S12 &v12,S13 &v13,S14 &v14,S15 &v15,index_type &p,op_type &op)
+                        :v1(v1),v2(v2),v3(v3),v4(v4),v5(v5),v6(v6),v7(v7),v8(v8),v9(v9),v10(v10),v11(v11), v12(v12),v13(v13),v14(v14),v15(v15),p(p),op(op)
+                {};
+                //! It call the copy function for each property
+                template<typename T>
+                inline void operator()(T& t) const
+                {
+                    op(v1.data.template get<T::value>().getVector().template get<0>(p),v2.data.template get<T::value>().getVector().template get<0>(p),v3.data.template get<T::value>().getVector().template get<0>(p),v4.data.template get<T::value>().getVector().template get<0>(p),v5.data.template get<T::value>().getVector().template get<0>(p),v6.data.template get<T::value>().getVector().template get<0>(p),v7.data.template get<T::value>().getVector().template get<0>(p),v8.data.template get<T::value>().getVector().template get<0>(p),v9.data.template get<T::value>().getVector().template get<0>(p),v10.data.template get<T::value>().getVector().template get<0>(p),v11.data.template get<T::value>().getVector().template get<0>(p),v12.data.template get<T::value>().getVector().template get<0>(p),v13.data.template get<T::value>().getVector().template get<0>(p),v14.data.template get<T::value>().getVector().template get<0>(p),v15.data.template get<T::value>().getVector().template get<0>(p));
+                }
+            };
+
+
+            template< class S1 , class S2 , class S3 , class S4,class S5,class S6 ,class S7,class S8, class S9, class S10, class S11, class S12, class S13, class S14, class S15, class Op >
+            static void for_each15( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4,S5 &s5,S6 &s6,S7 &s7,S8 &s8, S9 &s9, S10 &s10,S11 &s11,S12 &s12,S13 &s13,S14 &s14,S15 &s15, Op op )
+            {
+                for_each_prop_resize<S1,S2> the_resize(s1,s2);
+                boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
+                // ToDo : build checks, that the +-*/ operators are well defined
+                auto it=s1.data.template get<0>().getVector().getIterator();
+                while(it.isNext()){
+                    auto p=it.get();
+                    //converting to boost vector ids.
+                    for_each_prop15<S1,S2,S3,S4,S5,S6,S7,S8,S9,S10,S11,S12,S13,S14,S15,size_t,Op> cp(s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,s11,s12,s13,s14,s15,p,op);
+                    //creating an iterator on v_ids[0] [1] [2]
+                    boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
+
+                    ++it;
+                }
+            }
+
+
+
+
+
+
+           template<typename vector_type,typename index_type,typename norm_result_type>
+           struct for_each_norm
+           {
+               const vector_type &v;
+               index_type &p;
+               norm_result_type &n;
+               /*! \brief constructor
+                *
+                *
+                * \param src source encapsulated object
+                * \param dst destination encapsulated object
+                *
+                */
+               inline for_each_norm(const vector_type &v,index_type &p,norm_result_type &n)
+               :v(v),p(p),n(n)
+               {};
+               //! It call the copy function for each property
+               template<typename T>
+               inline void operator()(T& t) const
+               {
+                    if(fabs(v.data.template get<T::value>().getVector().template get<0>(p)) > n)
+                    {
+                        n=fabs(v.data.template get<T::value>().getVector().template get<0>(p));
+                    }
+
+               }
+           };
+
+            template< class S >
+            static typename boost::numeric::odeint::vector_space_norm_inf< S >::result_type norm_inf( const S &s )
+            {
+                typename boost::numeric::odeint::vector_space_norm_inf< S >::result_type n=0;
+                auto it=s.data.template get<0>().getVector().getIterator();
+                while(it.isNext()){
+                    auto p=it.get();
+                    //converting to boost vector ids.
+                    for_each_norm<S,size_t,typename boost::numeric::odeint::vector_space_norm_inf< S >::result_type> cp(s,p,n);
+                    //creating an iterator on v_ids[0] [1] [2]
+                    boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s.data)::max_prop>>(cp);
+
+                    ++it;
+                }
+                auto &v_cl = create_vcluster();
+                v_cl.max(n);
+                v_cl.execute();
+                //std::max();
+                //std::cout<<n<<std::endl;
+                return n;
+            }
+        };
+
+
+
+    } // odeint
+} // numeric
+} // boost
+
+
+
+
+
+
+
+
+#endif //OPENFPM_PDATA_VECTOR_ALGEBRA_OFP_HPP
diff --git a/src/Operators/Vector/vector_dist_operators.hpp b/src/Operators/Vector/vector_dist_operators.hpp
index c96188a5..9e148984 100644
--- a/src/Operators/Vector/vector_dist_operators.hpp
+++ b/src/Operators/Vector/vector_dist_operators.hpp
@@ -781,7 +781,7 @@ public:
 
 
 	//! return the result of the expression
-	template<typename r_type=typename std::remove_reference<decltype(-(o1.value(vect_dist_key_dx(0))))>::type > 
+	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));
@@ -823,6 +823,27 @@ struct vector_dist_expression_comp_sel<comp_dev,false>
 	typedef boost::mpl::int_<-1> type;
 };
 
+/*! \brief Expression implementation computation selector
+ *
+ */
+template<bool cond>
+struct vector_dist_expression_comp_proxy_sel
+{
+    template<bool cond_, typename v_type, typename exp_type>
+    static void compute(v_type &v,exp_type &v_exp)
+    { vector_dist_op_compute_op<0,false,vector_dist_expression_comp_sel<comp_dev,cond_>::type::value>
+        ::compute_expr(v,v_exp);}
+};
+template<>
+struct vector_dist_expression_comp_proxy_sel<false>
+{
+    template<bool cond, typename v_type, typename exp_type>
+    static void compute(v_type &v, exp_type &v_exp)
+    {   auto v_ker=v.toKernel();
+        vector_dist_op_compute_op<0,false,vector_dist_expression_comp_sel<comp_dev,cond>::type::value>
+        ::compute_expr(v_ker,v_exp);}
+};
+
 template<typename vector, bool is_ker = has_vector_kernel<vector>::type::value>
 struct vector_expression_transform
 {
@@ -1040,15 +1061,34 @@ public:
 	 * \return itself
 	 *
 	 */
-	template<typename T> vector & operator=(const vector_dist_expression<0,openfpm::vector<aggregate<T>>> & v_exp)
+	template<typename T,typename memory,template <typename> class layout_base > vector & operator=(const vector_dist_expression<0,openfpm::vector<aggregate<T>, memory, layout_base>> & 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);
+		//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);
+            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;
 	}
 
+    /*! \brief Fill the vector property with the evaluated expression
+ *
+ * \param v_exp expression to evaluate
+ *
+ * \return itself
+ *
+ */
+    template<typename T> vector & operator=(const vector_dist_expression<0,openfpm::vector_gpu<aggregate<T>>> & v_exp)
+    {
+            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;
+    }
+
 	/*! \brief Fill the vector property with the evaluated expression
 	 *
 	 * \param v_exp expression to evaluate
@@ -1139,45 +1179,45 @@ public:
  * \tparam vector involved
  *
  */
-template<typename T>
-class vector_dist_expression<0,openfpm::vector<aggregate<T>> >
+template<typename vector_type>
+class vector_dist_expression_impl
 {
-	//! Internal vector
-	typedef openfpm::vector<aggregate<T>> vector;
-
-	//! The temporal vector
-	mutable vector v;
+    //! Internal vector
+    typedef vector_type vector;
+    typedef typename boost::mpl::at<typename vector_type::value_type::type,boost::mpl::int_<0>>::type T;
+    //! The temporal vector
+    mutable vector v;
 
 public:
 
     typedef T * iterator;
     typedef const  T * const_iterator;
 
-	typedef typename has_vector_kernel<vector>::type is_ker;
+    typedef typename has_vector_kernel<vector>::type is_ker;
 
-	//! The type of the internal vector
-	typedef vector vtype;
+    //! The type of the internal vector
+    typedef vector vtype;
 
     //! The type of the internal value
     typedef T value_type;
 
     //! result for is sort
-	typedef boost::mpl::bool_<false> is_sort;
+    typedef boost::mpl::bool_<false> is_sort;
 
-	//! NN_type
-	typedef void NN_type;
+    //! NN_type
+    typedef void NN_type;
 
-	//! Property id of the point
-	static const unsigned int prop = 0;
+    //! Property id of the point
+    static const unsigned int prop = 0;
 
-	int var_id = 0;
+    int var_id = 0;
 
-	void setVarId(int var_id)
-	{
-		this->var_id = var_id;
-	}
+    void setVarId(int var_id)
+    {
+        this->var_id = var_id;
+    }
 
-	///////// BOOST ODEINT interface
+    ///////// BOOST ODEINT interface
     iterator begin()
     { return &v.template get<0>(0); }
 
@@ -1193,11 +1233,11 @@ public:
     size_t size() const
     { return v.size(); }
 
-	void resize(size_t n)
+    void resize(size_t n)
     {
-	    // Here
+        // Here
 
-	    v.resize(n);
+        v.resize(n);
     }
 
 /*	T * begin() {
@@ -1217,161 +1257,165 @@ public:
     //{ 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);
-	}
+    vector_dist_expression_impl()
+    {}
 
-	/*! \brief get the NN object
-	 *
-	 * \return the NN object
-	 *
-	 */
-	inline void * getNN() const
-	{
-		return NULL;
-	}
+    template<typename exp1, typename exp2, unsigned int op>
+    vector_dist_expression_impl(const vector_dist_expression_op<exp1,exp2,op> & v_exp)
+    {
+        this->operator=(v_exp);
+    }
 
-	/*! \brief Return the vector on which is acting
-	 *
-	 * It return the vector used in getVExpr, to get this object
-	 *
-	 * \return the vector
-	 *
-	 */
-	__device__ __host__ const vector & getVector() const
-	{
-		return v;
-	}
+    /*! \brief get the NN object
+     *
+     * \return the NN object
+     *
+     */
+    inline void * getNN() const
+    {
+        return NULL;
+    }
 
-	/*! \brief Return the vector on which is acting
-	 *
-	 * It return the vector used in getVExpr, to get this object
-	 *
-	 * \return the vector
-	 *
-	 */
-	__device__ __host__ vector & getVector()
-	{
-		return v;
-	}
+    /*! \brief Return the vector on which is acting
+     *
+     * It return the vector used in getVExpr, to get this object
+     *
+     * \return the vector
+     *
+     */
+    __device__ __host__ const vector & getVector() const
+    {
+        return v;
+    }
 
-	/*! \brief This function must be called before value
-	 *
-	 * it initialize the expression if needed
-	 *
-	 */
-	inline void init() const
-	{}
+    /*! \brief Return the vector on which is acting
+     *
+     * It return the vector used in getVExpr, to get this object
+     *
+     * \return the vector
+     *
+     */
+    __device__ __host__ vector & getVector()
+    {
+        return v;
+    }
 
-	/*! \brief Evaluate the expression
-	 *
-	 * \param k where to evaluate the expression
-	 *
-	 * \return the result of the expression
-	 *
-	 */
-	__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());
-	}
+    /*! \brief This function must be called before value
+     *
+     * it initialize the expression if needed
+     *
+     */
+    inline void init() const
+    {}
+
+    /*! \brief Evaluate the expression
+     *
+     * \param k where to evaluate the expression
+     *
+     * \return the result of the expression
+     *
+     */
+    __host__ __device__ inline auto value(const vect_dist_key_dx & k) const -> decltype(v.template get<0>(k.getKey()))
+    {
+        return v.template get<0>(k.getKey());
+    }
 
 
-	/*! \brief Fill the vector property with the evaluated expression
-	 *
-	 * \param v_exp expression to evaluate
-	 *
-	 * \return itself
-	 *
-	 */
-	template<unsigned int prp2, typename vector2> vector & operator=(const vector_dist_expression<prp2,vector2> & v_exp)
-	{
+    /*! \brief Fill the vector property with the evaluated expression
+     *
+     * \param v_exp expression to evaluate
+     *
+     * \return itself
+     *
+     */
+    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;
+            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());
+        constexpr bool cond=has_vector_kernel<vector>::type::value || std::is_same<vector,openfpm::vector<aggregate<T>,CudaMemory,memory_traits_inte>>::value;
+        //std::cout<<cond<<std::endl;
+        //std::cout<< (vector_dist_expression_comp_sel<comp_host,has_vector_kernel<vector>::type::value>::type::value || std::is_same<vector,openfpm::vector<aggregate<T>,CudaMemory,memory_traits_inte>>::value)<<std::endl;
+        //std::cout<<(vector_dist_expression_comp_sel<2,
+           //     has_vector_kernel<vector>::type::value>::type::value || std::is_same<vector,openfpm::vector<aggregate<T>,CudaMemory,memory_traits_inte>>::value)<<std::endl;
+        //std::cout<<has_vector_kernel<vector>::type::value<<std::endl;
+        //std::cout<<vector_dist_expression_comp_sel<2,false>::type::value<<std::endl;
+        //std::cout<<!std::is_same<vector,openfpm::vector<aggregate<T>,CudaMemory,memory_traits_inte>>::value<<std::endl;
+        if (has_vector_kernel<vector>::type::value == false && !std::is_same<vector,openfpm::vector<aggregate<T>,CudaMemory,memory_traits_inte>>::value)
+        {
+            vector_dist_op_compute_op<0,false,vector_dist_expression_comp_sel<comp_host,cond>::type::value>
+            ::compute_expr(v,v_exp);
+        }
+        else
+        {
+            vector_dist_expression_comp_proxy_sel<!std::is_same<vector,openfpm::vector<aggregate<T>,CudaMemory,memory_traits_inte>>::value>::template compute<cond>(v,v_exp);
+        }
 
-		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;
-	}
+        return v;
+    }
 
 
-	/*! \brief Fill the vector property with the evaluated expression
-	 *
-	 * \param v_exp expression to evaluate
-	 *
-	 * \return itself
-	 *
-	 */
-	template<typename exp1, typename exp2, unsigned int op>
-	vector & operator=(const vector_dist_expression_op<exp1,exp2,op> & v_exp)
-	{
+    /*! \brief Fill the vector property with the evaluated expression
+     *
+     * \param v_exp expression to evaluate
+     *
+     * \return itself
+     *
+     */
+    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;
+            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);
-		}
+        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;
-	}
+        return v;
+    }
 
-	/*! \brief Fill the vector property with the double
-	 *
-	 * \param d value to fill
-	 *
-	 * \return the internal vector
-	 *
-	 */
-	vector & operator=(double d)
-	{
-		std::cout << __FILE__ << ":" << __LINE__ << " Error: temporal with constants is unsupported" << std::endl;
-	}
+    /*! \brief Fill the vector property with the double
+     *
+     * \param d value to fill
+     *
+     * \return the internal vector
+     *
+     */
+    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";
+        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)
@@ -1386,8 +1430,56 @@ public:
     }
 };
 
+/*! \brief Sub class that encapsulate a vector properties operand to be used for expressions construction
+ *  Temporal Expressions
+ * \tparam prp property involved
+ * \tparam vector involved
+ *
+ */
+template<typename T, typename memory,template <typename> class layout_base >
+class vector_dist_expression<0,openfpm::vector<aggregate<T>,memory, layout_base> > : public vector_dist_expression_impl<openfpm::vector<aggregate<T>,memory, layout_base>>
+{
+    typedef openfpm::vector<aggregate<T>,memory, layout_base> vector;
+    typedef vector_dist_expression_impl<vector> base;
+public:
+    template<unsigned int prp2, typename vector2> vector & operator=(const vector_dist_expression<prp2,vector2> & v_exp)
+    {
+        return base::operator=(v_exp);
+    }
+    template<typename exp1, typename exp2, unsigned int op>
+    vector & operator=(const vector_dist_expression_op<exp1,exp2,op> & v_exp)
+    {
+        return base::operator=(v_exp);
+    }
+};
+
+/*! \brief Sub class that encapsulate a GPU vector properties operand to be used for expressions construction
+ *  Temporal Expressions
+ * \tparam prp property involved
+ * \tparam vector involved
+ *
+ */
+template<typename T>
+class vector_dist_expression<0,openfpm::vector_gpu<aggregate<T>>> : public vector_dist_expression_impl<openfpm::vector_gpu<aggregate<T>>>
+{
+    typedef openfpm::vector_gpu<aggregate<T>> vector;
+    typedef vector_dist_expression_impl<vector> base;
+public:
+    template<unsigned int prp2, typename vector2> vector & operator=(const vector_dist_expression<prp2,vector2> & v_exp)
+    {
+        return base::operator=(v_exp);
+    }
+    template<typename exp1, typename exp2, unsigned int op>
+    vector & operator=(const vector_dist_expression_op<exp1,exp2,op> & v_exp)
+    {
+        return base::operator=(v_exp);
+    }
+
+};
+
+template<typename T> using texp_v = vector_dist_expression<0,openfpm::vector<aggregate<T>>>;
+template<typename T> using texp_v_gpu = vector_dist_expression<0,openfpm::vector_gpu<aggregate<T>>>;
 
-template<typename T> using texp_v = vector_dist_expression<0,openfpm::vector<aggregate<T>> >;
 
 template<typename vector, unsigned int impl>
 struct switcher_get_v
@@ -1572,7 +1664,7 @@ public:
 	 * \return itself
 	 *
 	 */
-    template<typename T> vtype & operator=(const vector_dist_expression<0,openfpm::vector<aggregate<T>>> & v_exp)
+    template<typename T, typename memory> vtype & operator=(const vector_dist_expression<0,openfpm::vector<aggregate<T>,memory>> & v_exp)
     {
         v_exp.init();
 
-- 
GitLab