Fixing numeric with latest ExtPreAlloc

configure.ac

@@ -143,6 +143,12 @@ AC_CHECK_LIB(rt, clock_gettime, [AC_DEFINE([HAVE_CLOCK_GETTIME],[],[Have clock t
                                  OPT_LIBS="$OPT_LIBS -lrt"
                                 ])
 
+##########
+
+## Check for LIBHILBERT
+
+AX_LIB_HILBERT([],[echo "Cannot detect libhilbert, use the --with-libhilbert option if it is not installed in the default location"
+                    exit 210])
 
 ####### include OpenFPM_numerics include path"
 

src/Equations/stoke_flow_eq_3d.hpp

@@ -12,7 +12,6 @@
 
 constexpr unsigned int v[] = {0,1,2};
 constexpr unsigned int P = 3;
-constexpr unsigned int ic = 3;
 
 typedef Field<v[x],lid_nn_3d> v_x;
 typedef Field<v[y],lid_nn_3d> v_y;

src/FiniteDifference/FDScheme.hpp

@@ -454,13 +454,32 @@ public:
 	 */
 	template<typename T> void impose(const T & op , typename Sys_eqs::stype num ,long int id ,const long int (& start)[Sys_eqs::dims], const long int (& stop)[Sys_eqs::dims], bool skip_first = false)
 	{
-		// add padding to start and stop
-		grid_key_dx<Sys_eqs::dims> start_k = grid_key_dx<Sys_eqs::dims>(start);
-		grid_key_dx<Sys_eqs::dims> stop_k = grid_key_dx<Sys_eqs::dims>(stop);
+		grid_key_dx<Sys_eqs::dims> start_k;
+		grid_key_dx<Sys_eqs::dims> stop_k;
 
-		auto it = g_map.getSubDomainIterator(start_k,stop_k);
+        bool increment = false;
+        if (skip_first == true)
+        {
+                start_k = grid_key_dx<Sys_eqs::dims>(start);
+                stop_k = grid_key_dx<Sys_eqs::dims>(start);
+
+                auto it = g_map.getSubDomainIterator(start_k,stop_k);
+
+                if (it.isNext() == true)
+                        increment++;
+        }
+
+        // add padding to start and stop
+        start_k = grid_key_dx<Sys_eqs::dims>(start);
+        stop_k = grid_key_dx<Sys_eqs::dims>(stop);
+
+        auto it = g_map.getSubDomainIterator(start_k,stop_k);
+
+        if (increment == true)
+                ++it;
+
+        impose(op,num,id,it);
 
-		impose(op,num,id,it,skip_first);
 	}
 
 	/*! \brief Impose an operator
@@ -478,10 +497,8 @@ public:
 	 * \param it_d iterator that define where you want to impose
 	 *
 	 */
-	template<typename T> void impose(const T & op , typename Sys_eqs::stype num ,long int id ,grid_dist_iterator_sub<Sys_eqs::dims,typename g_map_type::d_grid> it_d, bool skip_first = false)
+	template<typename T> void impose(const T & op , typename Sys_eqs::stype num ,long int id ,grid_dist_iterator_sub<Sys_eqs::dims,typename g_map_type::d_grid> it_d)
 	{
-		Vcluster & v_cl = create_vcluster();
-
 		openfpm::vector<triplet> & trpl = A.getMatrixTriplets();
 
 		auto it = it_d;
@@ -489,20 +506,9 @@ public:
 
 		std::unordered_map<long int,float> cols;
 
-		bool is_first = skip_first;
-
 		// iterate all the grid points
 		while (it.isNext())
 		{
-			if (is_first == true && v_cl.getProcessUnitID() == 0)
-			{
-				++it;
-				is_first = false;
-				continue;
-			}
-			else
-				is_first = false;
-
 			// get the position
 			auto key = it.get();
 

src/FiniteDifference/eq.hpp

@@ -114,4 +114,10 @@ inline size_t mat_factor(size_t nvar, size_t sz, const size_t ord)
 	return nvar;
 }
 
+#include "mul.hpp"
+#include "Average.hpp"
+#include "Derivative.hpp"
+#include "sum.hpp"
+#include "Laplacian.hpp"
+
 #endif /* OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_EQ_HPP_ */

src/Makefile.am

 
-LINKLIBS =  $(PETSC_LIB) $(SUITESPARSE_LIBS) $(LAPACK_LIBS) $(BLAS_LIBS)  $(METIS_LIB) $(PARMETIS_LIB) $(DEFAULT_LIB) $(PTHREAD_LIBS) $(OPT_LIBS) $(BOOST_LDFLAGS) $(BOOST_PROGRAM_OPTIONS_LIB) $(BOOST_IOSTREAMS_LIB) $(HDF5_LDFLAGS)  $(HDF5_LIBS) $(LIBQUADMATH)
+LINKLIBS = $(LIBHILBERT_LIB) $(PETSC_LIB) $(SUITESPARSE_LIBS) $(LAPACK_LIBS) $(BLAS_LIBS)  $(METIS_LIB) $(PARMETIS_LIB) $(DEFAULT_LIB) $(PTHREAD_LIBS) $(OPT_LIBS) $(BOOST_LDFLAGS) $(BOOST_PROGRAM_OPTIONS_LIB) $(BOOST_IOSTREAMS_LIB) $(HDF5_LDFLAGS)  $(HDF5_LIBS) $(LIBQUADMATH)
 
 noinst_PROGRAMS = numerics
 numerics_SOURCES = main.cpp ../../openfpm_vcluster/src/VCluster.cpp ../../openfpm_devices/src/memory/HeapMemory.cpp ../../openfpm_devices/src/memory/PtrMemory.cpp ../../openfpm_devices/src/Memleak_check.cpp
-numerics_CXXFLAGS = $(AM_CXXFLAGS) $(HDF5_CPPFLAGS) $(INCLUDES_PATH) $(BOOST_CPPFLAGS) $(SUITESPARSE_INCLUDE) $(METIS_INCLUDE) $(PARMETIS_INCLUDE) $(EIGEN_INCLUDE) $(PETSC_INCLUDE) -Wno-deprecated-declarations -Wno-unused-local-typedefs
+numerics_CXXFLAGS = $(LIBHILBERT_INCLUDE) $(AM_CXXFLAGS) $(HDF5_CPPFLAGS) $(INCLUDES_PATH) $(BOOST_CPPFLAGS) $(SUITESPARSE_INCLUDE) $(METIS_INCLUDE) $(PARMETIS_INCLUDE) $(EIGEN_INCLUDE) $(PETSC_INCLUDE) -Wno-deprecated-declarations -Wno-unused-local-typedefs
 numerics_CFLAGS = $(CUDA_CFLAGS)
 numerics_LDADD = $(LINKLIBS) -lparmetis -lmetis
-nobase_include_HEADERS = PSE/Kernels.hpp PSE/Kernels_test_util.hpp Operators/Vector/vector_dist_operators_extensions.hpp Operators/Vector/vector_dist_operators.hpp Operators/Vector/vector_dist_operators_apply_kernel.hpp Operators/Vector/vector_dist_operators_functions.hpp
+nobase_include_HEADERS = Matrix/SparseMatrix.hpp Matrix/SparseMatrix_Eigen.hpp Matrix/SparseMatrix_petsc.hpp \
+Vector/Vector_eigen.hpp Vector/Vector_petsc.hpp Vector/Vector_util.hpp Vector/Vector.hpp \
+Solvers/umfpack_solver.hpp Solvers/petsc_solver.hpp \
+util/petsc_util.hpp util/linalgebra_lib.hpp util/util_num.hpp util/grid_dist_testing.hpp  \
+FiniteDifference/Average.hpp FiniteDifference/Derivative.hpp FiniteDifference/eq.hpp FiniteDifference/FDScheme.hpp FiniteDifference/Laplacian.hpp FiniteDifference/mul.hpp FiniteDifference/sum.hpp FiniteDifference/util/common.hpp \
+PSE/Kernels.hpp PSE/Kernels_test_util.hpp Operators/Vector/vector_dist_operators_extensions.hpp Operators/Vector/vector_dist_operators.hpp Operators/Vector/vector_dist_operators_apply_kernel.hpp Operators/Vector/vector_dist_operators_functions.hpp Operators/Vector/vector_dist_operator_assign.hpp
 
 .cu.o :
 	$(NVCC) $(NVCCFLAGS) -o $@ -c $<

src/Matrix/SparseMatrix.hpp

@@ -8,6 +8,9 @@
 #ifndef OPENFPM_NUMERICS_SRC_MATRIX_SPARSEMATRIX_HPP_
 #define OPENFPM_NUMERICS_SRC_MATRIX_SPARSEMATRIX_HPP_
 
+#include "config/config.h"
+#include "util/linalgebra_lib.hpp"
+
 #ifdef HAVE_EIGEN
 #include <Eigen/Sparse>
 #define DEFAULT_MATRIX  = EIGEN_BASE
@@ -93,7 +96,7 @@ public:
 	openfpm::vector<triplet_type> & getMatrixTriplets() {std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use this class you must compile OpenFPM with linear algebra support" << std::endl; return stub_vt;}
 	const int & getMat() const {std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use this class you must compile OpenFPM with linear algebra support" << std::endl; return stub_i;}
 	int & getMat() {std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use this class you must compile OpenFPM with linear algebra support" << std::endl; return stub_i;}
-	void resize(size_t row, size_t col) {std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use this class you must compile OpenFPM with linear algebra support" << std::endl;}
+	void resize(size_t row, size_t col, size_t row_n, size_t col_n) {std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use this class you must compile OpenFPM with linear algebra support" << std::endl;}
 	T operator()(id_t i, id_t j) {std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use this class you must compile OpenFPM with linear algebra support" << std::endl; return stub_t;}
 	bool save(const std::string & file) const {std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use this class you must compile OpenFPM with linear algebra support" << std::endl;return true;}
 	bool load(const std::string & file) {std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use this class you must compile OpenFPM with linear algebra support" << std::endl; return false;}

src/Matrix/SparseMatrix_Eigen.hpp

@@ -212,16 +212,13 @@ public:
 	 */
 	bool save(const std::string & file) const
 	{
-		std::vector<size_t> pap_prp;
+		size_t pap_prp;
 
 		Packer<decltype(trpl),HeapMemory>::packRequest(trpl,pap_prp);
 
-		// Calculate how much preallocated memory we need to pack all the objects
-		size_t req = ExtPreAlloc<HeapMemory>::calculateMem(pap_prp);
-
 		// allocate the memory
 		HeapMemory pmem;
-		pmem.allocate(req);
+		pmem.allocate(pap_prp);
 		ExtPreAlloc<HeapMemory> mem(pap_prp,pmem);
 
 		//Packing

src/Matrix/SparseMatrix_petsc.hpp

@@ -8,7 +8,7 @@
 #ifndef OPENFPM_NUMERICS_SRC_MATRIX_SPARSEMATRIX_PETSC_HPP_
 #define OPENFPM_NUMERICS_SRC_MATRIX_SPARSEMATRIX_PETSC_HPP_
 
-
+#include "util/petsc_util.hpp"
 #include "Vector/map_vector.hpp"
 #include <boost/mpl/int.hpp>
 #include <petscmat.h>
@@ -213,7 +213,8 @@ public:
 	~SparseMatrix()
 	{
 		// Destroy the matrix
-		PETSC_SAFE_CALL(MatDestroy(&mat));
+		if (is_openfpm_init() == true)
+		{PETSC_SAFE_CALL(MatDestroy(&mat));}
 	}
 
 	/*! \brief Get the Matrix triplets bugger

src/Matrix/SparseMatrix_unit_tests.hpp

@@ -407,7 +407,7 @@ BOOST_AUTO_TEST_CASE(sparse_matrix_eigen_petsc_solve)
 	}
 
 	// Get the petsc Matrix
-	Mat & matp = sm.getMat();
+	sm.getMat();
 
 
 	petsc_solver<double> solver;

src/Operators/Vector/vector_dist_operators.hpp

@@ -11,6 +11,7 @@
 #include "Vector/vector_dist.hpp"
 
 #define PROP_POS (unsigned int)-1
+#define PROP_CUSTOM (unsigned int)-2
 
 #define VECT_SUM 1
 #define VECT_SUB 2
@@ -23,12 +24,9 @@
 #define VECT_APPLYKER_IN_GEN 10
 #define VECT_APPLYKER_OUT_GEN 11
 #define VECT_APPLYKER_REDUCE_GEN 12
-#define VECT_APPLYKER_MULTI_IN 13
-#define VECT_APPLYKER_MULTI_OUT 14
-#define VECT_APPLYKER_MULTI_REDUCE 15
-#define VECT_APPLYKER_MULTI_IN_GEN 16
-#define VECT_APPLYKER_MULTI_OUT_GEN 17
-#define VECT_APPLYKER_MULTI_REDUCE_GEN 18
+#define VECT_APPLYKER_IN_SIM 13
+#define VECT_APPLYKER_OUT_SIM 14
+#define VECT_APPLYKER_REDUCE_SIM 15
 
 #define VECT_NORM 56
 #define VECT_NORM2 57
@@ -67,6 +65,8 @@
 #define VECT_PMUL 91
 #define VECT_SUB_UNI 92
 
+#define VECT_SUM_REDUCE 93
+
 
 /*! \brief has_init check if a type has defined a
  * method called init
@@ -338,10 +338,25 @@ class vector_dist_expression
 
 public:
 
+	typedef vector vtype;
+
+	//! Property id of the point
+	static const unsigned int prop = prp;
+
 	vector_dist_expression(vector & v)
 	:v(v)
 	{}
 
+	/*! \brief Return the vector on which is acting
+	 *
+	 * It return the vector used in getVExpr, to get this object
+	 *
+	 */
+	vector & getVector()
+	{
+		return v;
+	}
+
 	/*! \brief This function must be called before value
 	 *
 	 * it initialize the expression if needed
@@ -375,7 +390,7 @@ public:
 		{
 			auto key = it.get();
 
-			v.template getProp<prp>(key) = v_exp.value(key);
+			pos_or_prop<vector,prp>::value(v,key) = v_exp.value(key);
 
 			++it;
 		}
@@ -383,7 +398,6 @@ public:
 		return v;
 	}
 
-
 	/*! \brief Fill the vector property with the evaluated expression
 	 *
 	 * \param v_exp expression to evaluate
@@ -399,7 +413,7 @@ public:
 		{
 			auto key = it.get();
 
-			v.template getProp<prp>(key) = v_exp.value(key);
+			pos_or_prop<vector,prp>::value(v,key) = v_exp.value(key);
 
 			++it;
 		}
@@ -420,7 +434,7 @@ public:
 		{
 			auto key = it.get();
 
-			v.template getProp<prp>(key) = d;
+			pos_or_prop<vector,prp>::value(v,key) = d;
 
 			++it;
 		}
@@ -454,7 +468,7 @@ class vector_dist_expression<prp,double>
 
 public:
 
-	inline vector_dist_expression(double & d)
+	inline vector_dist_expression(const double & d)
 	:d(d)
 	{}
 
@@ -477,6 +491,79 @@ public:
 	}
 };
 
+/*! \brief Main class that encapsulate a float constant
+ *
+ * \param prp no meaning
+ *
+ */
+template<unsigned int prp>
+class vector_dist_expression<prp,float>
+{
+	float d;
+
+public:
+
+	typedef float vtype;
+
+	inline vector_dist_expression(const float & d)
+	:d(d)
+	{}
+
+	/*! \brief This function must be called before value
+	 *
+	 * it initialize the expression if needed
+	 *
+	 */
+	inline void init() const
+	{}
+
+	/*! \brief Evaluate the expression
+	 *
+	 * It just return the velue set in the constructor
+	 *
+	 */
+	inline float value(const vect_dist_key_dx & k) const
+	{
+		return d;
+	}
+};
+
+
+
+/*! \brief Main class that encapsulate a float constant
+ *
+ * \param prp no meaning
+ *
+ */
+template<typename T>
+class vector_dist_expression<PROP_CUSTOM,T>
+{
+
+public:
+
+	typedef float vtype;
+
+	inline vector_dist_expression()
+	{}
+
+	/*! \brief This function must be called before value
+	 *
+	 * it initialize the expression if needed
+	 *
+	 */
+	inline void init() const
+	{}
+
+	/*! \brief Evaluate the expression
+	 *
+	 * It just return the velue set in the constructor
+	 *
+	 */
+	inline T value(const vect_dist_key_dx & k) const
+	{
+		return T(0.0);
+	}
+};
 
 /* \brief sum two distributed vector expression
  *
@@ -1010,5 +1097,6 @@ operator/(const vector_dist_expression_op<exp1,exp2,op1> & va, const vector_dist
 #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_ */

src/Operators/Vector/vector_dist_operators_apply_kernel_unit_tests.hpp

@@ -51,6 +51,9 @@ template <typename vector,typename Kernel, typename NN_type> bool check_values_a
 
 		base2 += vd.template getProp<C>(p);
 
+		if (base1 != base2)
+			std::cout << "CAZZO " << base1 << "  " << base2 << std::endl;
+
 		ret &= base1 == base2;
 
 		++it2;
@@ -176,6 +179,58 @@ template <typename vector,typename Kernel, typename NN_type> bool check_values_a
 	return ret;
 }
 
+template <typename vector,typename Kernel, typename NN_type> bool check_values_apply_kernel3(vector & vd, Kernel & ker, NN_type & NN)
+{
+	bool ret = true;
+	auto it = vd.getDomainIterator();
+
+	// ### WIKI 13 ###
+	//
+	// Check that apply kernel work
+	//
+	auto it2 = vd.getDomainIterator();
+	while (it2.isNext())
+	{
+		Point<3,float> base2 = 0.0;
+		auto p = it2.get();
+
+		Point<3,float> xp = vd.getPos(p);
+
+		Point<3,float> base1 = vd.template getProp<VA>(p);
+
+		Point<3,float> prp_x = vd.template getProp<VC>(p);
+
+		// For each neighborhood particle
+		auto Np = NN.template getNNIterator<NO_CHECK>(NN.getCell(vd.getPos(p)));
+
+		while (Np.isNext())
+		{
+			// Neighborhood particle q
+			auto q = Np.get();
+
+			if (q == p.getKey())	{++Np; continue;};
+
+			// position q
+			Point<3,float> xq = vd.getPos(q);
+
+			Point<3,float> prp_y = vd.template getProp<VC>(q);
+
+			base2 += ker.value(xp,xq,prp_x,prp_y);
+
+			++Np;
+		}
+
+		base2 += vd.template getProp<VC>(p);
+
+		ret &= base1 == base2;
+
+		++it2;
+	}
+
+	BOOST_REQUIRE_EQUAL(ret,true);
+
+	return ret;
+}
 
 template <typename vector,typename Kernel, typename NN_type> bool check_values_apply_kernel2_reduce(vector & vd, Kernel & ker, NN_type & NN)
 {
@@ -238,6 +293,51 @@ template <typename vector,typename Kernel, typename NN_type> bool check_values_a
 	return ret;
 }
 
+template <typename vector,typename Kernel, typename NN_type> bool check_values_apply_kernel3_reduce(vector & vd, Kernel & ker, NN_type & NN, const Point<2,float> & p)
+{
+	bool ret = true;
+	auto it = vd.getDomainIterator();
+
+	Point<2,float> base2 = 0.0;
+
+	auto it2 = vd.getDomainIterator();
+	while (it2.isNext())
+	{
+		auto p = it2.get();
+
+		Point<3,float> xp = vd.getPos(p);
+
+		Point<2,float> ker_accu = 0.0;
+
+		// For each neighborhood particle
+		auto Np = NN.template getNNIterator<NO_CHECK>(NN.getCell(vd.getPos(p)));
+
+		while (Np.isNext())
+		{
+			// Neighborhood particle q
+			auto q = Np.get();
+
+			if (q == p.getKey())	{++Np; continue;};
+
+			// position q
+			Point<3,float> xq = vd.getPos(q);
+
+			ker_accu += ker.value(xp,xq);
+
+			++Np;
+		}
+
+		base2 += ker_accu;
+
+		++it2;
+	}
+
+	BOOST_REQUIRE_EQUAL(p.get(0),base2.get(0));
+	BOOST_REQUIRE_EQUAL(p.get(1),base2.get(1));
+
+	return ret;
+}
+
 BOOST_AUTO_TEST_CASE( vector_dist_operators_apply_kernel_test )
 {
 	if (create_vcluster().getProcessingUnits() > 3)
@@ -257,7 +357,6 @@ BOOST_AUTO_TEST_CASE( vector_dist_operators_apply_kernel_test )
 	vector_dist<3,float,aggregate<float,float,float,VectorS<3,float>,VectorS<3,float>,VectorS<3,float>,float>> vd(512,box,bc,ghost);
 
 	auto vA = getV<A>(vd);
-	auto vB = getV<B>(vd);
 	auto vC = getV<C>(vd);
 
 	auto vVA = getV<VA>(vd);
@@ -281,23 +380,34 @@ BOOST_AUTO_TEST_CASE( vector_dist_operators_apply_kernel_test )
 	vVA = applyKernel_in(2.0*vVC + vVB ,vd,cl,ker) + vVC;
 	check_values_apply_kernel2(vd,ker,cl);
 
-	vA = applyKernel_reduce(vVC * vVB + norm(vVB),vd,cl,ker) + vC;
+	vA = rsum(applyKernel_in(vVC * vVB + norm(vVB),vd,cl,ker),vd) + vC;
 	check_values_apply_kernel_reduce(vd,ker,cl);
 
-	vVA = applyKernel_reduce(2.0*vVC + vVB ,vd,cl,ker) + vVC;
+	vVA = rsum(applyKernel_in(2.0*vVC + vVB ,vd,cl,ker),vd) + vVC;
 	check_values_apply_kernel2_reduce(vd,ker,cl);
 
-/*	vA = applyKernel_in_gen(vVC * vVB + norm(vVB),vd,cl,ker) + vC;
+	vA = applyKernel_in_gen(vVC * vVB + norm(vVB),vd,cl,ker) + vC;
 	check_values_apply_kernel(vd,ker,cl);
 
 	vVA = applyKernel_in_gen(2.0*vVC + vVB ,vd,cl,ker) + vVC;
 	check_values_apply_kernel2(vd,ker,cl);
 
-	vA = applyKernel_reduce_gen(vVC * vVB + norm(vVB),vd,cl,ker) + vC;
+	vA = rsum(applyKernel_in_gen(vVC * vVB + norm(vVB),vd,cl,ker),vd) + vC;
 	check_values_apply_kernel_reduce(vd,ker,cl);
 
-	vVA = applyKernel_reduce_gen(2.0*vVC + vVB ,vd,cl,ker) + vVC;
-	check_values_apply_kernel2_reduce(vd,ker,cl);*/
+	vVA = rsum(applyKernel_in_gen(2.0*vVC + vVB ,vd,cl,ker),vd) + vVC;
+	check_values_apply_kernel2_reduce(vd,ker,cl);
+
+	// Check it compile the code is the same
+	vVA = applyKernel_in_gen(vVC,vd,cl,ker) + vVC;
+	check_values_apply_kernel3(vd,ker,cl);
+
+	vVA = applyKernel_in    (vVC,vd,cl,ker) + vVC;
+	check_values_apply_kernel3(vd,ker,cl);
+
+	// just check it compile
+	Point<2,float> p = rsum(applyKernel_in_sim(vd,cl,ker),vd).get();
+	check_values_apply_kernel3_reduce(vd,ker,cl,p);
 }
 
 

src/Operators/Vector/vector_dist_operators_functions.hpp

@@ -212,4 +212,76 @@ fun_name(double d, const vector_dist_expression_op<exp1,exp2,op1> & va)\
 
 CREATE_VDIST_ARG2_FUNC(pmul,pmul,VECT_PMUL)
 
+
+////////// Special function reduce /////////////////////////
+
+
+template <typename exp1, typename vector_type>
+class vector_dist_expression_op<exp1,vector_type,VECT_SUM_REDUCE>
+{
+	const exp1 o1;
+
+	typedef typename apply_kernel_rtype<decltype(o1.value(vect_dist_key_dx(0)))>::rtype rtype;
+
+	// calculated value
+	mutable typename std::remove_reference<rtype>::type val;
+
+	const vector_type & vd;
+
+public:
+
+	vector_dist_expression_op(const exp1 & o1, const vector_type & vd)
+	:o1(o1), vd(vd)
+	{}
+
+	inline void init() const
+	{
+		o1.init();
+
+		val = 0.0;
+
+		auto it = vd.getDomainIterator();
+
+		while (it.isNext())
+		{
+			auto key = it.get();
+
+			val += o1.value(key);
+
+			++it;
+		}
+	}
+
+	inline typename std::remove_reference<rtype>::type get()
+	{
+		init();
+		return value(vect_dist_key_dx(0));
+	}
+
+	template<typename r_type= typename std::remove_reference<rtype>::type > inline r_type value(const vect_dist_key_dx & key) const
+	{
+		return val;
+	}
+};
+
+
+template<typename exp1, typename exp2_, unsigned int op1, typename vector_type>
+inline vector_dist_expression_op<vector_dist_expression_op<exp1,exp2_,op1>,vector_type,VECT_SUM_REDUCE>
+rsum(const vector_dist_expression_op<exp1,exp2_,op1> & va, const vector_type & vd)
+{
+	vector_dist_expression_op<vector_dist_expression_op<exp1,exp2_,op1>,vector_type,VECT_SUM_REDUCE> exp_sum(va,vd);
+
+	return exp_sum;
+}
+
+template<unsigned int prp1, typename v1, typename vector_type>
+inline vector_dist_expression_op<vector_dist_expression<prp1,v1>,vector_type,VECT_SUM_REDUCE>
+rsum(const vector_dist_expression<prp1,v1> & va, const vector_type & vd)