diff --git a/src/Graph/dist_map_graph.hpp b/src/Graph/dist_map_graph.hpp
index 208ea06936b5323e44e7103892610a9f52f6bcd4..2a181a8c40428926d6955c2558c8eaf2af32d3ca 100644
--- a/src/Graph/dist_map_graph.hpp
+++ b/src/Graph/dist_map_graph.hpp
@@ -572,43 +572,42 @@ class DistGraph_CSR
 			size_t pc = i;
 
 			size_t vp_size = sgp.get(pc).send_v.size();
-			std::vector<size_t> pap_prp;
+
+			size_t req = 0;
 
 			// prepare slot for number of vertices
-			Packer<size_t, HeapMemory>::packRequest(pap_prp);
+			Packer<size_t, HeapMemory>::packRequest(req);
 
 			for (size_t j = 0; j < vp_size; j++)
 			{
 				// prepare slot for vertex
-				Packer<V, HeapMemory>::packRequest(pap_prp);
+				Packer<V, HeapMemory>::packRequest(req);
 
 				// prepare slot info for vertex
-				Packer<v_info, HeapMemory>::packRequest(pap_prp);
+				Packer<v_info, HeapMemory>::packRequest(req);
 
 				// prepare slot for the number of children
-				Packer<size_t, HeapMemory>::packRequest(pap_prp);
+				Packer<size_t, HeapMemory>::packRequest(req);
 
 				// prepare slots for the children
 				for (size_t k = 0; k < sgp.get(pc).send_es.get(j); k++)
 				{
 					// prepare slot for edge
-					Packer<E, HeapMemory>::packRequest(pap_prp);
+					Packer<E, HeapMemory>::packRequest(req);
 
 					// prepare slot for edge info
-					Packer<e_info, HeapMemory>::packRequest(pap_prp);
+					Packer<e_info, HeapMemory>::packRequest(req);
 
 					// prepare slot for edge target id
-					Packer<size_t, HeapMemory>::packRequest(pap_prp);
+					Packer<size_t, HeapMemory>::packRequest(req);
 				}
 			}
 
-			// Calculate how much preallocated memory we need to pack all the objects for each vector
-			ExtPreAlloc<HeapMemory>::calculateMem(pap_prp);
-
 			// allocate the memory
 			HeapMemory & pmem = *(new HeapMemory());
 //			pmem.allocate(req);
-			ExtPreAlloc<HeapMemory> & mem = *(new ExtPreAlloc<HeapMemory>(pap_prp, pmem));
+			ExtPreAlloc<HeapMemory> & mem = *(new ExtPreAlloc<HeapMemory>(req, pmem));
+
 			mem.incRef();
 
 			Pack_stat sts;
@@ -646,7 +645,7 @@ class DistGraph_CSR
 
 			prc.add(i);
 			size.add(pmem.size());
-			ptr.add(mem.getPointer(0));
+			ptr.add(mem.getPointerBase());
 		}
 
 		// Exchange informations through processors
diff --git a/src/Grid/grid_dist_id.hpp b/src/Grid/grid_dist_id.hpp
index 593766d29f81eb69386cee266170a65ee422c5d2..54d187aae12ae267a4537ab0b06095707d08800f 100644
--- a/src/Grid/grid_dist_id.hpp
+++ b/src/Grid/grid_dist_id.hpp
@@ -1198,8 +1198,8 @@ public:
 		// Convert the local external ghost boxes into grid unit boxes
 		create_local_eg_box();
 
-		// total number of sending vector
-		std::vector<size_t> pap_prp;
+
+		size_t req = 0;
 
 		// Create a packing request vector
 		for ( size_t i = 0 ; i < ig_box.size() ; i++ )
@@ -1218,19 +1218,20 @@ public:
 				g_ig_box -= gdb_ext.get(sub_id).origin.template convertPoint<size_t>();
 
 				// Pack a size_t for the internal ghost id
-				Packer<size_t,HeapMemory>::packRequest(pap_prp);
+				Packer<size_t,HeapMemory>::packRequest(req);
 				// Create a sub grid iterator spanning the internal ghost layer
 				grid_key_dx_iterator_sub<dim> sub_it(loc_grid.get(sub_id).getGrid(),g_ig_box.getKP1(),g_ig_box.getKP2());
 				// and pack the internal ghost grid
-				Packer<device_grid,HeapMemory>::template packRequest<prp...>(loc_grid.get(sub_id),sub_it,pap_prp);
+				Packer<device_grid,HeapMemory>::template packRequest<prp...>(loc_grid.get(sub_id),sub_it,req);
 			}
 		}
 
 		// resize the property buffer memory
-		g_send_prp_mem.resize(ExtPreAlloc<Memory>::calculateMem(pap_prp));
+		g_send_prp_mem.resize(req);
 
 		// Create an object of preallocated memory for properties
-		ExtPreAlloc<Memory> & prAlloc_prp = *(new ExtPreAlloc<Memory>(pap_prp,g_send_prp_mem));
+		ExtPreAlloc<Memory> & prAlloc_prp = *(new ExtPreAlloc<Memory>(req,g_send_prp_mem));
+
 		prAlloc_prp.incRef();
 
 		// Pack information
@@ -1239,7 +1240,9 @@ public:
 		// Pack the information for each processor and send it
 		for ( size_t i = 0 ; i < ig_box.size() ; i++ )
 		{
+
 			sts.mark();
+			void * pointer = prAlloc_prp.getPointerEnd();
 
 			// for each ghost box
 			for (size_t j = 0 ; j < ig_box.get(i).bid.size() ; j++)
@@ -1264,7 +1267,12 @@ public:
 				Packer<device_grid,HeapMemory>::template pack<prp...>(prAlloc_prp,loc_grid.get(sub_id),sub_it,sts);
 			}
 			// send the request
-			v_cl.send(ig_box.get(i).prc,0,sts.getMarkPointer(prAlloc_prp),sts.getMarkSize(prAlloc_prp));
+
+			void * pointer2 = prAlloc_prp.getPointerEnd();
+
+			v_cl.send(ig_box.get(i).prc,0,pointer/*sts.getMarkPointer(prAlloc_prp)*/,(char *)pointer2 - (char *)pointer /*sts.getMarkSize(prAlloc_prp)*/);
+
+//			pointer = prAlloc_prp.getPointerEnd();
 		}
 
 		// Calculate the total information to receive from each processors
@@ -1288,7 +1296,8 @@ public:
 		g_recv_prp_mem.resize(ExtPreAlloc<Memory>::calculateMem(prp_recv));
 
 		// Create an object of preallocated memory for properties
-		ExtPreAlloc<Memory> & prRecv_prp = *(new ExtPreAlloc<Memory>(prp_recv,g_recv_prp_mem));
+		ExtPreAlloc<Memory> & prRecv_prp = *(new ExtPreAlloc<Memory>(req,g_recv_prp_mem));
+
 		prRecv_prp.incRef();
 
 		// queue the receives
diff --git a/src/Makefile.am b/src/Makefile.am
index e6809078ee7bd1fc2cf631db4be93907b16a3513..1928568150c0a7178b5457bcf69b1ac79cb0e3b8 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -4,8 +4,8 @@ noinst_PROGRAMS = pdata
 pdata_SOURCES = main.cpp Grid/grid_dist_id_unit_test.cpp  lib/pdata.cpp test_multiple_o.cpp ../openfpm_devices/src/memory/HeapMemory.cpp ../openfpm_devices/src/memory/PtrMemory.cpp ../openfpm_vcluster/src/VCluster.cpp ../openfpm_devices/src/Memleak_check.cpp
 pdata_CXXFLAGS = $(PETSC_INCLUDE) $(HDF5_CPPFLAGS) $(CUDA_CFLAGS) $(INCLUDES_PATH) $(PARMETIS_INCLUDE) $(METIS_INCLUDE) $(BOOST_CPPFLAGS) $(H5PART_INCLUDE) -DPARALLEL_IO  -Wno-unused-local-typedefs
 pdata_CFLAGS = $(CUDA_CFLAGS)
-pdata_LDADD = $(LINKLIBS) -lmetis -lparmetis
-nobase_include_HEADERS = Decomposition/CartDecomposition.hpp Decomposition/CartDecomposition_ext.hpp Decomposition/common.hpp Decomposition/Decomposition.hpp  Decomposition/ie_ghost.hpp \
+pdata_LDADD = $(LINKLIBS) -lparmetis -lmetis -L/usr/local/libhilbert/lib -lhilbert
+nobase_include_HEADERS = Decomposition/CartDecomposition.hpp Decomposition/common.hpp Decomposition/Decomposition.hpp  Decomposition/ie_ghost.hpp \
          Decomposition/nn_processor.hpp Decomposition/ie_loc_ghost.hpp Decomposition/ORB.hpp \
          Graph/CartesianGraphFactory.hpp \
          Grid/grid_dist_id.hpp Grid/grid_dist_id_iterator_dec.hpp Grid/grid_dist_util.hpp  Grid/grid_dist_id_iterator_sub.hpp Grid/grid_dist_id_iterator.hpp Grid/grid_dist_key.hpp Grid/staggered_dist_grid.hpp Grid/staggered_dist_grid_util.hpp Grid/staggered_dist_grid_copy.hpp \
diff --git a/src/Vector/vector_dist.hpp b/src/Vector/vector_dist.hpp
index 9bd6989c2b971c1c1e06404380ab9487a0eab993..14e76723819960ed1e7de05b4dffef3e624b4d69 100644
--- a/src/Vector/vector_dist.hpp
+++ b/src/Vector/vector_dist.hpp
@@ -17,6 +17,7 @@
 #include "memory/PreAllocHeapMemory.hpp"
 #include "memory/PtrMemory.hpp"
 #include "NN/CellList/CellList.hpp"
+#include "NN/CellList/CellListFast_hilb.hpp"
 #include "util/common.hpp"
 #include "util/object_util.hpp"
 #include "memory/ExtPreAlloc.hpp"
@@ -24,6 +25,7 @@
 #include "VTKWriter/VTKWriter.hpp"
 #include "Decomposition/common.hpp"
 #include "Grid/grid_dist_id_iterator_dec.hpp"
+#include "Grid/grid_key_dx_iterator_hilbert.hpp"
 #include "Vector/vector_dist_ofb.hpp"
 #include "Decomposition/CartDecomposition.hpp"
 #include "data_type/aggregate.hpp"
@@ -773,6 +775,35 @@ private:
 		}
 	}
 
+	/*! \brief Calculate parameters for the cell list
+	 *
+	 * \param div Division array
+	 * \param r_cut interation radius or size of each cell
+	 * \param enlarge In case of padding particles the cell list must be enlarged, like a ghost. This parameter says how much must be enlarged
+	 *
+	 * \return the processor bounding box
+	 */
+	inline Box<dim, St> cl_param_calculate(size_t (&div)[dim], St r_cut, const Ghost<dim, St> & enlarge)
+
+	{
+		// calculate the parameters of the cell list
+
+		// get the processor bounding box
+		Box<dim, St> pbox = dec.getProcessorBounds();
+
+		// extend by the ghost
+		pbox.enlarge(enlarge);
+
+		// Calculate the division array and the cell box
+		for (size_t i = 0; i < dim; i++)
+		{
+			div[i] = static_cast<size_t>((pbox.getP2().get(i) - pbox.getP1().get(i)) / r_cut);
+			div[i]++;
+			pbox.setHigh(i,pbox.getLow(i) + div[i]*r_cut);
+		}
+		return pbox;
+	}
+
 public:
 
 	//! space type
@@ -1120,7 +1151,11 @@ public:
 
 		// Create and fill send buffer for particle properties
 
-		ExtPreAlloc<Memory> * prAlloc_prp = new ExtPreAlloc<Memory>(pap_prp, g_prp_mem);
+////////////////////////////////////////////////
+		size_t req = ExtPreAlloc<Memory>::calculateMem(pap_prp);
+
+		ExtPreAlloc<Memory> * prAlloc_prp = new ExtPreAlloc<Memory>(req, g_prp_mem);
+/////////////////////////////////////////////////
 		openfpm::vector<send_vector> g_send_prp;
 		fill_send_ghost_prp_buf<send_vector, prp_object, prp...>(g_send_prp, prAlloc_prp);
 
@@ -1130,7 +1165,11 @@ public:
 		openfpm::vector<send_pos_vector> g_pos_send;
 		if (opt != NO_POSITION)
 		{
-			prAlloc_pos = new ExtPreAlloc<Memory>(pap_pos, g_pos_mem);
+////////////////////////////////////////////////
+			size_t req1 = ExtPreAlloc<Memory>::calculateMem(pap_pos);
+
+			prAlloc_pos = new ExtPreAlloc<Memory>(req1, g_pos_mem);
+////////////////////////////////////////////////////////
 			fill_send_ghost_pos_buf(g_pos_send, prAlloc_pos);
 		}
 
@@ -1264,6 +1303,24 @@ public:
 		return getCellList(r_cut, g);
 	}
 
+	/*! \brief Construct an hilbert cell list starting from the stored particles
+	 *
+	 * \tparam CellL CellList type to construct
+	 *
+	 * \param r_cut interation radius, or size of each cell
+	 *
+	 * \return the Cell list
+	 *
+	 */
+	template<typename CellL = CellList_hilb<dim, St, FAST, shift<dim, St> > > CellL getCellList_hilb(St r_cut)
+	{
+		// Get ghost and anlarge by 1%
+		Ghost<dim,St> g = dec.getGhost();
+		g.magnify(1.01);
+
+		return getCellList_hilb(r_cut, g);
+	}
+
 	/*! \brief Update a cell list using the stored particles
 	 *
 	 * \tparam CellL CellList type to construct
@@ -1307,24 +1364,43 @@ public:
 	{
 		CellL cell_list;
 
-		// calculate the parameters of the cell list
+		// Division array
+		size_t div[dim];
 
-		// get the processor bounding box
-		Box<dim, St> pbox = dec.getProcessorBounds();
-		// extend by the ghost
-		pbox.enlarge(enlarge);
+		// Processor bounding box
+		auto pbox = cl_param_calculate(div, r_cut, enlarge);
 
+		cell_list.Initialize(pbox, div);
+
+		updateCellList(cell_list);
+
+		return cell_list;
+	}
+
+	/*! \brief Construct an hilbert cell list starting from the stored particles
+	 *
+	 * It differ from the get getCellList for an additional parameter, in case the
+	 * domain + ghost is not big enough to contain additional padding particles, a Cell list
+	 * with bigger space can be created
+	 * (padding particles in general are particles added by the user out of the domains)
+	 *
+	 * \tparam CellL CellList type to construct
+	 *
+	 * \param r_cut interation radius, or size of each cell
+	 * \param enlarge In case of padding particles the cell list must be enlarged, like a ghost this parameter say how much must be enlarged
+	 *
+	 */
+	template<typename CellL = CellList_hilb<dim, St, FAST, shift<dim, St> > > CellL getCellList_hilb(St r_cut, const Ghost<dim, St> & enlarge)
+	{
+		CellL cell_list;
+
+		// Division array
 		size_t div[dim];
 
-		// Calculate the division array and the cell box
-		for (size_t i = 0; i < dim; i++)
-		{
-			div[i] = static_cast<size_t>((pbox.getP2().get(i) - pbox.getP1().get(i)) / r_cut);
-			div[i]++;
-			pbox.setHigh(i,pbox.getLow(i) + div[i]*r_cut);
-		}
+		// Processor bounding box
+		auto pbox = cl_param_calculate(div, r_cut, enlarge);
 
-		cell_list.Initialize(pbox, div);
+		cell_list.Initialize(pbox, div, g_m);
 
 		updateCellList(cell_list);
 
@@ -1388,6 +1464,116 @@ public:
 		}
 	}
 
+	/*! \brief Construct a cell list starting from the stored particles and reorder a vector according to the Hilberts curve
+	 *
+	 * \tparam CellL CellList type to construct
+	 *
+	 * \param m an order of a hilbert curve
+	 *
+	 *
+	 *
+	 */
+	template<typename CellL=CellList<dim,St,FAST,shift<dim,St> > > void reorder (int32_t m)
+	{
+		reorder(m,dec.getGhost());
+	}
+
+
+	/*! \brief Construct a cell list starting from the stored particles and reorder a vector according to the Hilberts curve
+	 *
+	 *
+	 *It differs from the reorder(m) for an additional parameter, in case the
+	 * domain + ghost is not big enough to contain additional padding particles, a Cell list
+	 * with bigger space can be created
+	 * (padding particles in general are particles added by the user out of the domains)
+	 *
+	 * \param m order of a curve
+	 * \param enlarge In case of padding particles the cell list must be enlarged, like a ghost this parameter say how much must be enlarged
+	 *
+	 */
+	template<typename CellL=CellList<dim,St,FAST,shift<dim,St> > > void reorder(int32_t m, const Ghost<dim,St> & enlarge)
+	{
+		// reset the ghost part
+		v_pos.resize(g_m);
+		v_prp.resize(g_m);
+
+
+		CellL cell_list;
+
+		// calculate the parameters of the cell list
+
+		// get the processor bounding box
+		Box<dim,St> pbox = dec.getProcessorBounds();
+		// extend by the ghost
+		pbox.enlarge(enlarge);
+
+		size_t div[dim];
+
+		// Calculate the division array and the cell box
+		for (size_t i = 0 ; i < dim ; i++)
+		{
+			div[i] = 1 << m;
+		}
+
+		cell_list.Initialize(pbox,div);
+
+		// for each particle add the particle to the cell list
+
+		auto it = getIterator();
+
+		while (it.isNext())
+		{
+			auto key = it.get();
+
+			cell_list.add(this->getPos(key),key.getKey());
+
+			++it;
+		}
+
+		// Use cell_list to reorder v_pos
+
+		//destination vector
+		openfpm::vector<Point<dim,St>> v_pos_dest;
+		openfpm::vector<prop> v_prp_dest;
+
+		v_pos_dest.resize(v_pos.size());
+		v_prp_dest.resize(v_prp.size());
+
+		//hilberts curve iterator
+		grid_key_dx_iterator_hilbert<dim> h_it(m);
+
+		//Index for v_pos_dest
+		size_t count = 0;
+
+		grid_key_dx<dim> ksum;
+
+		for (size_t i = 0; i < dim ; i++)
+			ksum.set_d(i,cell_list.getPadding(i));
+
+		while (h_it.isNext())
+		{
+		  auto key = h_it.get();
+		  key += ksum;
+
+		  size_t lin = cell_list.getGrid().LinId(key);
+
+		  // for each particle in the Cell "lin"
+		  for (size_t i = 0; i < cell_list.getNelements(lin); i++)
+		  {
+			  //reorder
+			  auto v = cell_list.get(lin,i);
+			  v_pos_dest.get(count) = v_pos.get(v);
+			  v_prp_dest.get(count) = v_prp.get(v);
+
+			  count++;
+		  }
+		  ++h_it;
+		}
+
+		v_pos.swap(v_pos_dest);
+		v_prp.swap(v_prp_dest);
+	}
+
 	/*! \brief It return the number of particles contained by the previous processors
 	 *
 	 * \Warning It only work with the initial decomposition
diff --git a/src/Vector/vector_dist_iterator.hpp b/src/Vector/vector_dist_iterator.hpp
index 795e4290b9c2d24067599d891cefe8f3842007e1..4b8ee861f90e1fa5955065df98e53598d462813e 100644
--- a/src/Vector/vector_dist_iterator.hpp
+++ b/src/Vector/vector_dist_iterator.hpp
@@ -74,6 +74,15 @@ class vector_dist_iterator
 	{
 		return vect_dist_key_dx(v_it);
 	}
+
+	/*! \brief Reset the iterator
+	 *
+	 *
+	 */
+	void reset()
+	{
+		v_it = 0;
+	}
 };
 
 
diff --git a/src/Vector/vector_dist_unit_test.hpp b/src/Vector/vector_dist_unit_test.hpp
index 7b09bbac89dd93425b26f619dbe7aa1665bb5a3e..1439ea1e2c3c8bc5a043b16ef5c393c1eaa47ac0 100644
--- a/src/Vector/vector_dist_unit_test.hpp
+++ b/src/Vector/vector_dist_unit_test.hpp
@@ -10,6 +10,8 @@
 
 #include <random>
 #include "Vector/vector_dist.hpp"
+#include "data_type/aggregate.hpp"
+#include "Vector/vector_dist_performance_util.hpp"
 
 /*! \brief Count the total number of particles
  *
@@ -1035,6 +1037,8 @@ BOOST_AUTO_TEST_CASE( vector_dist_periodic_test_interacting_particles )
 			size_t cnt = total_n_part_lc(vd,bc);
 
 			BOOST_REQUIRE_EQUAL((size_t)k,cnt);
+
+
 		}
 	}
 }
@@ -1159,6 +1163,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_cell_verlet_test )
 	}
 }
 
+
 BOOST_AUTO_TEST_CASE( vector_dist_periodic_map_list )
 {
 	typedef Point<3,float> s;
@@ -1284,6 +1289,274 @@ BOOST_AUTO_TEST_CASE( vector_dist_periodic_map_list )
 	}
 }
 
+///////////////////////// test hilb ///////////////////////////////
+
+BOOST_AUTO_TEST_CASE( vector_dist_reorder_2d_test )
+{
+	typedef Point<2,float> s;
+
+	Vcluster & v_cl = create_vcluster();
+
+    // set the seed
+	// create the random generator engine
+	std::srand(v_cl.getProcessUnitID());
+    std::default_random_engine eg;
+    std::uniform_real_distribution<float> ud(0.0f, 1.0f);
+
+    long int k = 524288 * v_cl.getProcessingUnits();
+
+	long int big_step = k / 4;
+	big_step = (big_step == 0)?1:big_step;
+
+	print_test_v( "Testing 2D vector with hilbert curve reordering k<=",k);
+
+	// 2D test
+	for ( ; k >= 2 ; k-= decrement(k,big_step) )
+	{
+		BOOST_TEST_CHECKPOINT( "Testing 2D vector with hilbert curve reordering k=" << k );
+
+		Box<2,float> box({0.0,0.0},{1.0,1.0});
+
+		// Boundary conditions
+		size_t bc[2]={NON_PERIODIC,NON_PERIODIC};
+
+		vector_dist<2,float, Point_test<float>, CartDecomposition<2,float> > vd(k,box,bc,Ghost<2,float>(0.0));
+
+		auto it = vd.getIterator();
+
+		while (it.isNext())
+		{
+			auto key = it.get();
+
+			vd.getPos(key)[0] = ud(eg);
+			vd.getPos(key)[1] = ud(eg);
+
+			++it;
+		}
+
+		vd.map();
+
+		// Create first cell list
+
+		auto NN1 = vd.getCellList(0.01);
+
+		//An order of a curve
+		int32_t m = 6;
+
+		//Reorder a vector
+		vd.reorder(m);
+
+		// Create second cell list
+		auto NN2 = vd.getCellList(0.01);
+
+		//Check equality of cell sizes
+		for (size_t i = 0 ; i < NN1.getGrid().size() ; i++)
+		{
+			size_t n1 = NN1.getNelements(i);
+			size_t n2 = NN2.getNelements(i);
+
+			BOOST_REQUIRE_EQUAL(n1,n2);
+		}
+	}
+}
+
+BOOST_AUTO_TEST_CASE( vector_dist_cl_random_vs_hilb_forces_test )
+{
+	///////////////////// INPUT DATA //////////////////////
+
+	// Dimensionality of the space
+	const size_t dim = 3;
+	// Cut-off radiuses. Can be put different number of values
+	openfpm::vector<float> cl_r_cutoff {0.05};
+	// The starting amount of particles (remember that this number is multiplied by number of processors you use for testing)
+	size_t cl_k_start = 10000;
+	// The lower threshold for number of particles
+	size_t cl_k_min = 1000;
+	// Ghost part of distributed vector
+	double ghost_part = 0.05;
+
+	///////////////////////////////////////////////////////
+
+	//For different r_cut
+	for (size_t r = 0; r < cl_r_cutoff.size(); r++ )
+	{
+		Vcluster & v_cl = create_vcluster();
+
+		//Cut-off radius
+		float r_cut = cl_r_cutoff.get(r);
+
+		//Number of particles
+		size_t k = cl_k_start * v_cl.getProcessingUnits();
+
+		std::string str("Testing " + std::to_string(dim) + "D vector's forces (random vs hilb celllist) k<=");
+
+		vector_dist_test::print_test_v(str,k);
+
+		//For different number of particles
+		for (size_t k_int = k ; k_int >= cl_k_min ; k_int/=2 )
+		{
+			BOOST_TEST_CHECKPOINT( "Testing " << dim << "D vector's forces (random vs hilb celllist) k<=" << k_int );
+
+			Box<dim,float> box;
+
+			for (size_t i = 0; i < dim; i++)
+			{
+				box.setLow(i,0.0);
+				box.setHigh(i,1.0);
+			}
+
+			// Boundary conditions
+			size_t bc[dim];
+
+			for (size_t i = 0; i < dim; i++)
+				bc[i] = PERIODIC;
+
+			vector_dist<dim,float, aggregate<float[dim]>, CartDecomposition<dim,float> > vd(k_int,box,bc,Ghost<dim,float>(ghost_part));
+
+			vector_dist<dim,float, aggregate<float[dim]>, CartDecomposition<dim,float> > vd2(k_int,box,bc,Ghost<dim,float>(ghost_part));
+
+			// Initialize dist vectors
+			vd_initialize_double<dim>(vd, vd2, v_cl, k_int);
+
+			vd.template ghost_get<0>();
+			vd2.template ghost_get<0>();
+
+			//Get a cell list
+
+			auto NN = vd.getCellList(r_cut);
+
+			//Calculate forces
+
+			calc_forces<dim>(NN,vd,r_cut);
+
+			//Get a cell list hilb
+
+			auto NN_hilb = vd2.getCellList_hilb(r_cut);
+
+			//Calculate forces
+			calc_forces_hilb<dim>(NN_hilb,vd2,r_cut);
+
+			auto it_v = vd.getIterator();
+
+			while (it_v.isNext())
+			{
+				//key
+				vect_dist_key_dx key = it_v.get();
+
+				for (size_t i = 0; i < dim; i++)
+				{
+					auto a1 = vd.template getProp<0>(key)[i];
+					auto a2 = vd2.template getProp<0>(key)[i];
+					//Check that the forces are equal
+					BOOST_REQUIRE_CLOSE(a1,a2,1);
+				}
+
+				++it_v;
+			}
+		}
+	}
+}
+
+BOOST_AUTO_TEST_CASE( vector_dist_cl_random_vs_reorder_forces_test )
+{
+	///////////////////// INPUT DATA //////////////////////
+
+	// Dimensionality of the space
+	const size_t dim = 3;
+	// Cut-off radiuses. Can be put different number of values
+	openfpm::vector<float> cl_r_cutoff {0.01};
+	// The starting amount of particles (remember that this number is multiplied by number of processors you use for testing)
+	size_t cl_k_start = 10000;
+	// The lower threshold for number of particles
+	size_t cl_k_min = 1000;
+	// Ghost part of distributed vector
+	double ghost_part = 0.01;
+
+	///////////////////////////////////////////////////////
+
+	//For different r_cut
+	for (size_t r = 0; r < cl_r_cutoff.size(); r++ )
+	{
+		Vcluster & v_cl = create_vcluster();
+
+		//Cut-off radius
+		float r_cut = cl_r_cutoff.get(r);
+
+		//Number of particles
+		size_t k = cl_k_start * v_cl.getProcessingUnits();
+
+		std::string str("Testing " + std::to_string(dim) + "D vector's forces (random vs reorder) k<=");
+
+		vector_dist_test::print_test_v(str,k);
+
+		//For different number of particles
+		for (size_t k_int = k ; k_int >= cl_k_min ; k_int/=2 )
+		{
+			BOOST_TEST_CHECKPOINT( "Testing " << dim << "D vector's forces (random vs reorder) k<=" << k_int );
+
+			Box<dim,float> box;
+
+			for (size_t i = 0; i < dim; i++)
+			{
+				box.setLow(i,0.0);
+				box.setHigh(i,1.0);
+			}
+
+			// Boundary conditions
+			size_t bc[dim];
+
+			for (size_t i = 0; i < dim; i++)
+				bc[i] = PERIODIC;
+
+			vector_dist<dim,float, aggregate<float[dim], float[dim]>, CartDecomposition<dim,float> > vd(k_int,box,bc,Ghost<dim,float>(ghost_part));
+
+			// Initialize vd
+			vd_initialize<dim,decltype(vd)>(vd, v_cl, k_int);
+
+			vd.template ghost_get<0>();
+
+			//Get a cell list
+
+			auto NN1 = vd.getCellList(r_cut);
+
+			//Calculate forces '0'
+
+			calc_forces<dim>(NN1,vd,r_cut);
+
+			//Reorder and get a cell list again
+
+			vd.reorder(4);
+
+			vd.template ghost_get<0>();
+
+			auto NN2 = vd.getCellList(r_cut);
+
+			//Calculate forces '1'
+			calc_forces<dim,1>(NN2,vd,r_cut);
+
+			//Test for equality of forces
+			auto it_v = vd.getDomainIterator();
+
+			while (it_v.isNext())
+			{
+				//key
+				vect_dist_key_dx key = it_v.get();
+
+				for (size_t i = 0; i < dim; i++)
+				{
+					auto a1 = vd.template getProp<0>(key)[i];
+					auto a2 = vd.template getProp<1>(key)[i];
+					//Check that the forces are (almost) equal
+					BOOST_REQUIRE_CLOSE(a1,a2,1);
+				}
+
+				++it_v;
+			}
+		}
+	}
+}
+
+
 BOOST_AUTO_TEST_SUITE_END()
 
 #endif /* VECTOR_DIST_UNIT_TEST_HPP_ */
diff --git a/src/main.cpp b/src/main.cpp
index 3da9319f0a056e40d032dc90e54177035c9fa4fa..103a485104501343d28a7ebbc503d12baac747e0 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -37,6 +37,9 @@ int main(int argc, char* argv[])
 #include "Decomposition/Distribution/metis_util_unit_test.hpp"
 #include "dec_optimizer_unit_test.hpp"
 #include "Vector/vector_dist_unit_test.hpp"
+#ifdef PERFORMANCE_TEST
+#include "pdata_performance.hpp"
+#endif
 #include "Decomposition/Distribution/Distribution_unit_tests.hpp"
 //#include "DLB/DLB_unit_test.hpp"
 #include "Graph/dist_map_graph_unit_test.hpp"