diff --git a/configure.ac b/configure.ac
index 7842383e6fd7081c82e7f147d74448ec5983b6fb..1796d15ef9d3e9cdc6b863d36b2c9b0380e5022b 100755
--- a/configure.ac
+++ b/configure.ac
@@ -122,13 +122,6 @@ if test x"$with_hdf5" = x"no"; then
     exit 207
 fi
 
-
-##########
-
-## Check for PETSC
-
-AX_LIB_PETSC()
-
 #########
 
 ###### Check for test coverage
@@ -314,6 +307,12 @@ AX_LAPACK([],[])
 
 AX_SUITESPARSE([],[])
 
+##########
+
+## Check for PETSC
+
+AX_LIB_PETSC()
+
 ###### Checking for EIGEN
 
 AX_EIGEN([],[])
diff --git a/m4/ax_petsc_lib.m4 b/m4/ax_petsc_lib.m4
index 1bd159091052c10c40e9e9317a4cca0ec579c4ec..ba7fff8942f745cb3fbd11e3374c85fbb074dbe5 100755
--- a/m4/ax_petsc_lib.m4
+++ b/m4/ax_petsc_lib.m4
@@ -105,7 +105,7 @@ AC_DEFUN([AX_LIB_PETSC], [
 			AX_OPENMP([CFLAGS="$OPENMP_CFLAGS"
 				   LDFLAGS="$OPENMP_LDFLAGS"],[])
                         CFLAGS="$CFLAGS -I$with_petsc/include $HDF5_INCLUDE $METIS_INCLUDE "
-                        LDFLAGS="$LDFLAGS -L$with_petsc/lib $HDF5_LDFLAGS  $HDF5_LIBS $METIS_LIB -lmetis "
+                        LDFLAGS="$LDFLAGS -L$with_petsc/lib $HDF5_LDFLAGS  $HDF5_LIBS $METIS_LIB -lmetis $SUITESPARSE_LIBS"
 			CC=$CXX
 
                         AC_LANG_SAVE
diff --git a/src/DMatrix/EMatrix.hpp b/src/DMatrix/EMatrix.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..278e779e413a0c60548597501be0f326180e4e87
--- /dev/null
+++ b/src/DMatrix/EMatrix.hpp
@@ -0,0 +1,431 @@
+/*
+ * EMatrix.hpp
+ *
+ *  Created on: Feb 12, 2018
+ *      Author: i-bird
+ */
+
+#ifndef OPENFPM_DATA_SRC_SPACE_EMATRIX_HPP_
+#define OPENFPM_DATA_SRC_SPACE_EMATRIX_HPP_
+
+#include "config.h"
+
+#ifdef HAVE_EIGEN
+#include <Eigen/Dense>
+#include "memory/ExtPreAlloc.hpp"
+#include "memory/HeapMemory.hpp"
+#include "Packer_Unpacker/Packer_util.hpp"
+#include "Packer_Unpacker/Packer.hpp"
+#include "Packer_Unpacker/Unpacker.hpp"
+
+template<typename _Scalar, int _Rows, int _Cols,
+         int _Options = Eigen::AutoAlign |
+#if EIGEN_GNUC_AT(3,4)
+    // workaround a bug in at least gcc 3.4.6
+    // the innermost ?: ternary operator is misparsed. We write it slightly
+    // differently and this makes gcc 3.4.6 happy, but it's ugly.
+    // The error would only show up with EIGEN_DEFAULT_TO_ROW_MAJOR is defined
+    // (when EIGEN_DEFAULT_MATRIX_STORAGE_ORDER_OPTION is RowMajor)
+                          ( (_Rows==1 && _Cols!=1) ? Eigen::RowMajor
+                          : !(_Cols==1 && _Rows!=1) ?  EIGEN_DEFAULT_MATRIX_STORAGE_ORDER_OPTION
+                          : Eigen::ColMajor ),
+#else
+                          ( (_Rows==1 && _Cols!=1) ? Eigen::RowMajor
+                          : (_Cols==1 && _Rows!=1) ? Eigen::ColMajor
+                          : EIGEN_DEFAULT_MATRIX_STORAGE_ORDER_OPTION ),
+#endif
+         int _MaxRows = _Rows,
+         int _MaxCols = _Cols
+> class EMatrix;
+
+namespace Eigen {
+
+namespace internal {
+template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
+struct traits<EMatrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols> >
+{
+private:
+  enum { size = internal::size_at_compile_time<_Rows,_Cols>::ret };
+  typedef typename find_best_packet<_Scalar,size>::type PacketScalar;
+  enum {
+      row_major_bit = _Options&RowMajor ? RowMajorBit : 0,
+      is_dynamic_size_storage = _MaxRows==Dynamic || _MaxCols==Dynamic,
+      max_size = is_dynamic_size_storage ? Dynamic : _MaxRows*_MaxCols,
+      default_alignment = compute_default_alignment<_Scalar,max_size>::value,
+      actual_alignment = ((_Options&DontAlign)==0) ? default_alignment : 0,
+      required_alignment = unpacket_traits<PacketScalar>::alignment,
+      packet_access_bit = (packet_traits<_Scalar>::Vectorizable && (EIGEN_UNALIGNED_VECTORIZE || (actual_alignment>=required_alignment))) ? PacketAccessBit : 0
+    };
+
+public:
+  typedef _Scalar Scalar;
+  typedef Dense StorageKind;
+  typedef Eigen::Index StorageIndex;
+  typedef MatrixXpr XprKind;
+  enum {
+    RowsAtCompileTime = _Rows,
+    ColsAtCompileTime = _Cols,
+    MaxRowsAtCompileTime = _MaxRows,
+    MaxColsAtCompileTime = _MaxCols,
+    Flags = compute_matrix_flags<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>::ret,
+    Options = _Options,
+    InnerStrideAtCompileTime = 1,
+    OuterStrideAtCompileTime = (Options&RowMajor) ? ColsAtCompileTime : RowsAtCompileTime,
+
+    // FIXME, the following flag in only used to define NeedsToAlign in PlainObjectBase
+    EvaluatorFlags = LinearAccessBit | DirectAccessBit | packet_access_bit | row_major_bit,
+    Alignment = actual_alignment
+  };
+};
+}
+}
+
+/*! \brief it is just a wrapper for Eigen matrix compatible for openfpm serialization
+ *
+ *
+ */
+template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
+class EMatrix : public Eigen::Matrix<_Scalar,_Rows,_Cols,_Options,_MaxRows,_MaxCols>
+{
+
+	typedef typename Eigen::Matrix<_Scalar,_Rows,_Cols,_Options,_MaxRows,_MaxCols>::Base Base;
+
+public:
+
+	//typedef typename Eigen::internal::traits<Eigen::Matrix<_Scalar,_Rows,_Cols>>::Scalar Scalar;
+
+	EIGEN_DENSE_PUBLIC_INTERFACE(EMatrix)
+
+	/*! \brief It calculate the number of byte required to serialize the object
+	 *
+	 * \tparam prp list of properties
+	 *
+	 * \param req reference to the total counter required to pack the information
+	 *
+	 */
+	template<int ... prp> inline void packRequest(size_t & req) const
+	{
+		// Memory required to serialize the Matrix
+		req += sizeof(_Scalar)*this->rows()*this->cols() + 2*sizeof(_Scalar);
+	}
+
+
+	/*! \brief pack a vector selecting the properties to pack
+	 *
+	 * \param mem preallocated memory where to pack the vector
+	 * \param sts pack-stat info
+	 *
+	 */
+	template<int ... prp> inline void pack(ExtPreAlloc<HeapMemory> & mem, Pack_stat & sts) const
+	{
+		//Pack the number of rows and colums
+		Packer<size_t, HeapMemory>::pack(mem,this->rows(),sts);
+		Packer<size_t, HeapMemory>::pack(mem,this->cols(),sts);
+
+		mem.allocate(sizeof(_Scalar)*this->rows()*this->cols());
+		void * dst = mem.getPointer();
+		void * src = (void *)this->data();
+
+		memcpy(dst,src,sizeof(_Scalar)*this->rows()*this->cols());
+		sts.incReq();
+	}
+
+	/*! \brief unpack a vector
+	 *
+	 * \param mem preallocated memory from where to unpack the vector
+	 * \param ps unpack-stat info
+	 */
+	template<int ... prp> inline void unpack(ExtPreAlloc<HeapMemory> & mem, Unpack_stat & ps)
+	{
+		size_t row;
+		size_t col;
+
+		//Pack the number of rows and colums
+		Unpacker<size_t, HeapMemory>::unpack(mem,row,ps);
+		Unpacker<size_t, HeapMemory>::unpack(mem,col,ps);
+
+		this->resize(row,col);
+
+		void * dst = this->data();
+		void * src = mem.getPointerOffset(ps.getOffset());
+
+		memcpy(dst,src,sizeof(_Scalar)*row*col);
+
+		ps.addOffset(sizeof(_Scalar)*row*col);
+	}
+
+	////// wrap all eigen operator
+
+    /** \brief Default constructor.
+      *
+      * For fixed-size matrices, does nothing.
+      *
+      * For dynamic-size matrices, creates an empty matrix of size 0. Does not allocate any array. Such a matrix
+      * is called a null matrix. This constructor is the unique way to create null matrices: resizing
+      * a matrix to 0 is not supported.
+      *
+      * \sa resize(Index,Index)
+      */
+    EIGEN_DEVICE_FUNC
+    EIGEN_STRONG_INLINE EMatrix()
+    :Eigen::Matrix<_Scalar,_Rows,_Cols,_Options,_MaxRows,_MaxCols>()
+    {}
+
+    /**
+      * \brief Assigns matrices to each other.
+      *
+      * \note This is a special case of the templated operator=. Its purpose is
+      * to prevent a default operator= from hiding the templated operator=.
+      *
+      * \callgraph
+      */
+    EIGEN_DEVICE_FUNC
+    EIGEN_STRONG_INLINE EMatrix& operator=(const Eigen::Matrix<_Scalar,_Rows,_Cols,_Options,_MaxRows,_MaxCols>& other)
+    {
+    	((Eigen::Matrix<_Scalar,_Rows,_Cols,_Options,_MaxRows,_MaxCols> *)this)->operator=(other);
+    	return *this;
+    }
+
+    /** \internal
+      * \brief Copies the value of the expression \a other into \c *this with automatic resizing.
+      *
+      * *this might be resized to match the dimensions of \a other. If *this was a null matrix (not already initialized),
+      * it will be initialized.
+      *
+      * Note that copying a row-vector into a vector (and conversely) is allowed.
+      * The resizing, if any, is then done in the appropriate way so that row-vectors
+      * remain row-vectors and vectors remain vectors.
+      */
+    template<typename OtherDerived>
+    EIGEN_DEVICE_FUNC
+    EIGEN_STRONG_INLINE EMatrix& operator=(const Eigen::DenseBase<OtherDerived>& other)
+    {
+    	((Eigen::Matrix<_Scalar,_Rows,_Cols,_Options,_MaxRows,_MaxCols> *)this)->operator=(other);
+    	return *this;
+    }
+
+    /**
+      * \brief Copies the generic expression \a other into *this.
+      * \copydetails DenseBase::operator=(const EigenBase<OtherDerived> &other)
+      */
+    template<typename OtherDerived>
+    EIGEN_DEVICE_FUNC
+    EIGEN_STRONG_INLINE EMatrix& operator=(const Eigen::EigenBase<OtherDerived> &other)
+    {
+    	((Eigen::Matrix<_Scalar,_Rows,_Cols,_Options,_MaxRows,_MaxCols> *)this)->operator=(other);
+    	return *this;
+    }
+
+    template<typename OtherDerived>
+    EIGEN_DEVICE_FUNC
+    EIGEN_STRONG_INLINE EMatrix& operator=(const Eigen::ReturnByValue<OtherDerived>& func)
+    {
+    	((Eigen::Matrix<_Scalar,_Rows,_Cols,_Options,_MaxRows,_MaxCols> *)this)->operator=(func);
+    	return *this;
+    }
+
+#if EIGEN_HAS_RVALUE_REFERENCES
+    EIGEN_DEVICE_FUNC
+    EMatrix(Eigen::Matrix<_Scalar,_Rows,_Cols,_Options,_MaxRows,_MaxCols> && other) EIGEN_NOEXCEPT_IF(std::is_nothrow_move_constructible<Scalar>::value)
+    :Eigen::Matrix<_Scalar,_Rows,_Cols,_Options,_MaxRows,_MaxCols>(other)
+    {}
+
+    EIGEN_DEVICE_FUNC
+    EMatrix& operator=(Eigen::Matrix<_Scalar,_Rows,_Cols,_Options,_MaxRows,_MaxCols> && other) EIGEN_NOEXCEPT_IF(std::is_nothrow_move_assignable<Scalar>::value)
+    {return *this;}
+
+#endif
+
+    #ifndef EIGEN_PARSED_BY_DOXYGEN
+
+    // This constructor is for both 1x1 matrices and dynamic vectors
+    template<typename T>
+    EIGEN_DEVICE_FUNC
+    EIGEN_STRONG_INLINE explicit EMatrix(const T& x)
+    :Eigen::Matrix<_Scalar,_Rows,_Cols,_Options,_MaxRows,_MaxCols>(x)
+    {}
+
+    template<typename T0, typename T1>
+    EIGEN_DEVICE_FUNC
+    EIGEN_STRONG_INLINE EMatrix(const T0& x, const T1& y)
+    :Eigen::Matrix<_Scalar,_Rows,_Cols,_Options,_MaxRows,_MaxCols>(x,y)
+    {}
+
+    #else
+    /** \brief Constructs a fixed-sized matrix initialized with coefficients starting at \a data */
+    EIGEN_DEVICE_FUNC
+    explicit EMatrix(const Scalar *data)
+    :Eigen::Matrix<_Scalar,_Rows,_Cols,_Options,_MaxRows,_MaxCols>(data)
+	 {}
+
+    /** \brief Constructs a vector or row-vector with given dimension. \only_for_vectors
+      *
+      * This is useful for dynamic-size vectors. For fixed-size vectors,
+      * it is redundant to pass these parameters, so one should use the default constructor
+      * Matrix() instead.
+      *
+      * \warning This constructor is disabled for fixed-size \c 1x1 matrices. For instance,
+      * calling Matrix<double,1,1>(1) will call the initialization constructor: Matrix(const Scalar&).
+      * For fixed-size \c 1x1 matrices it is therefore recommended to use the default
+      * constructor Matrix() instead, especially when using one of the non standard
+      * \c EIGEN_INITIALIZE_MATRICES_BY_{ZERO,\c NAN} macros (see \ref TopicPreprocessorDirectives).
+      */
+    EIGEN_STRONG_INLINE explicit EMatrix(Index dim)
+    :Eigen::Matrix<_Scalar,_Rows,_Cols,_Options,_MaxRows,_MaxCols>(dim)
+	 {}
+
+    /** \brief Constructs an initialized 1x1 matrix with the given coefficient */
+    EMatrix(const Scalar& x)
+    :Eigen::Matrix<_Scalar,_Rows,_Cols,_Options,_MaxRows,_MaxCols>(x)
+    {}
+    /** \brief Constructs an uninitialized matrix with \a rows rows and \a cols columns.
+      *
+      * This is useful for dynamic-size matrices. For fixed-size matrices,
+      * it is redundant to pass these parameters, so one should use the default constructor
+      * Matrix() instead.
+      *
+      * \warning This constructor is disabled for fixed-size \c 1x2 and \c 2x1 vectors. For instance,
+      * calling Matrix2f(2,1) will call the initialization constructor: Matrix(const Scalar& x, const Scalar& y).
+      * For fixed-size \c 1x2 or \c 2x1 vectors it is therefore recommended to use the default
+      * constructor Matrix() instead, especially when using one of the non standard
+      * \c EIGEN_INITIALIZE_MATRICES_BY_{ZERO,\c NAN} macros (see \ref TopicPreprocessorDirectives).
+      */
+    EIGEN_DEVICE_FUNC
+    EMatrix(Index rows, Index cols)
+    :Eigen::Matrix<_Scalar,_Rows,_Cols,_Options,_MaxRows,_MaxCols>(rows,cols)
+	 {}
+
+    /** \brief Constructs an initialized 2D vector with given coefficients */
+    EMatrix(const Scalar& x, const Scalar& y)
+    :Eigen::Matrix<_Scalar,_Rows,_Cols,_Options,_MaxRows,_MaxCols>(x,y)
+	 {};
+    #endif
+
+    /** \brief Constructs an initialized 3D vector with given coefficients */
+    EIGEN_DEVICE_FUNC
+    EIGEN_STRONG_INLINE EMatrix(const Scalar& x, const Scalar& y, const Scalar& z)
+    :Eigen::Matrix<_Scalar,_Rows,_Cols,_Options,_MaxRows,_MaxCols>(x,y,z)
+    {}
+
+    /** \brief Constructs an initialized 4D vector with given coefficients */
+    EIGEN_DEVICE_FUNC
+    EIGEN_STRONG_INLINE EMatrix(const Scalar& x, const Scalar& y, const Scalar& z, const Scalar& w)
+    :Eigen::Matrix<_Scalar,_Rows,_Cols,_Options,_MaxRows,_MaxCols>(x,y,z,w)
+    {}
+
+
+    /** \brief Copy constructor */
+    EIGEN_DEVICE_FUNC
+    EIGEN_STRONG_INLINE EMatrix(const Eigen::Matrix<_Scalar,_Rows,_Cols,_Options,_MaxRows,_MaxCols> & other)
+    :Eigen::Matrix<_Scalar,_Rows,_Cols,_Options,_MaxRows,_MaxCols>(other)
+    {}
+
+    /** \brief Copy constructor for generic expressions.
+      * \sa MatrixBase::operator=(const EigenBase<OtherDerived>&)
+      */
+    template<typename OtherDerived>
+    EIGEN_DEVICE_FUNC
+    EIGEN_STRONG_INLINE EMatrix(const Eigen::EigenBase<OtherDerived> &other)
+    :Eigen::Matrix<_Scalar,_Rows,_Cols,_Options,_MaxRows,_MaxCols>(other)
+    {}
+};
+
+// Standard typedef from eigen
+typedef EMatrix< std::complex< double >, 2, 2 > 	EMatrix2cd;
+typedef EMatrix< std::complex< float >, 2, 2 > 	EMatrix2cf;
+typedef EMatrix< double, 2, 2 > 	EMatrix2d;
+typedef EMatrix< float, 2, 2 > 	EMatrix2f;
+typedef EMatrix< int, 2, 2 > 	EMatrix2i;
+typedef EMatrix< std::complex< double >, 2, Eigen::Dynamic > 	EMatrix2Xcd;
+typedef EMatrix< std::complex< float >, 2, Eigen::Dynamic > 	EMatrix2Xcf;
+typedef EMatrix< double, 2, Eigen::Dynamic > 	EMatrix2Xd;
+typedef EMatrix< float, 2, Eigen::Dynamic > 	EMatrix2Xf;
+typedef EMatrix< int, 2, Eigen::Dynamic > 	EMatrix2Xi;
+typedef EMatrix< std::complex< double >, 3, 3 > 	EMatrix3cd;
+typedef EMatrix< std::complex< float >, 3, 3 > 	EMatrix3cf;
+typedef EMatrix< double, 3, 3 > 	EMatrix3d;
+typedef EMatrix< float, 3, 3 > 	EMatrix3f;
+typedef EMatrix< int, 3, 3 > 	EMatrix3i;
+typedef EMatrix< std::complex< double >, 3, Eigen::Dynamic > 	EMatrix3Xcd;
+typedef EMatrix< std::complex< float >, 3, Eigen::Dynamic > 	EMatrix3Xcf;
+typedef EMatrix< double, 3, Eigen::Dynamic > 	EMatrix3Xd;
+typedef EMatrix< float, 3, Eigen::Dynamic > 	EMatrix3Xf;
+typedef EMatrix< int, 3, Eigen::Dynamic > 	EMatrix3Xi;
+typedef EMatrix< std::complex< double >, 4, 4 > 	EMatrix4cd;
+typedef EMatrix< std::complex< float >, 4, 4 > 	EMatrix4cf;
+typedef EMatrix< double, 4, 4 > 	EMatrix4d;
+typedef EMatrix< float, 4, 4 > 	EMatrix4f;
+typedef EMatrix< int, 4, 4 > 	EMatrix4i;
+typedef EMatrix< std::complex< double >, 4, Eigen::Dynamic > 	EMatrix4Xcd;
+typedef EMatrix< std::complex< float >, 4, Eigen::Dynamic > 	EMatrix4Xcf;
+typedef EMatrix< double, 4, Eigen::Dynamic > 	EMatrix4Xd;
+typedef EMatrix< float, 4, Eigen::Dynamic > 	EMatrix4Xf;
+typedef EMatrix< int, 4, Eigen::Dynamic > 	EMatrix4Xi;
+typedef EMatrix< std::complex< double >, Eigen::Dynamic, 2 > 	EMatrixX2cd;
+typedef EMatrix< std::complex< float >, Eigen::Dynamic, 2 > 	EMatrixX2cf;
+typedef EMatrix< double, Eigen::Dynamic, 2 > 	EMatrixX2d;
+typedef EMatrix< float, Eigen::Dynamic, 2 > 	EMatrixX2f;
+typedef EMatrix< int, Eigen::Dynamic, 2 > 	EMatrixX2i;
+typedef EMatrix< std::complex< double >, Eigen::Dynamic, 3 > 	EMatrixX3cd;
+typedef EMatrix< std::complex< float >, Eigen::Dynamic, 3 > 	EMatrixX3cf;
+typedef EMatrix< double, Eigen::Dynamic, 3 > 	EMatrixX3d;
+typedef EMatrix< float, Eigen::Dynamic, 3 > 	EMatrixX3f;
+typedef EMatrix< int, Eigen::Dynamic, 3 > 	EMatrixX3i;
+typedef EMatrix< std::complex< double >, Eigen::Dynamic, 4 > 	EMatrixX4cd;
+typedef EMatrix< std::complex< float >, Eigen::Dynamic, 4 > 	EMatrixX4cf;
+typedef EMatrix< double, Eigen::Dynamic, 4 > 	EMatrixX4d;
+typedef EMatrix< float, Eigen::Dynamic, 4 > 	EMatrixX4f;
+typedef EMatrix< int, Eigen::Dynamic, 4 > 	EMatrixX4i;
+typedef EMatrix< std::complex< double >, Eigen::Dynamic, Eigen::Dynamic > 	EMatrixXcd;
+typedef EMatrix< std::complex< float >, Eigen::Dynamic, Eigen::Dynamic > 	EMatrixXcf;
+typedef EMatrix< double, Eigen::Dynamic, Eigen::Dynamic > 	EMatrixXd;
+typedef EMatrix< float, Eigen::Dynamic, Eigen::Dynamic > 	EMatrixXf;
+typedef EMatrix< int, Eigen::Dynamic, Eigen::Dynamic > 	EMatrixXi;
+typedef EMatrix< std::complex< double >, 1, 2 > 	ERowVector2cd;
+typedef EMatrix< std::complex< float >, 1, 2 > 	ERowVector2cf;
+typedef EMatrix< double, 1, 2 > 	ERowVector2d;
+typedef EMatrix< float, 1, 2 > 	ERowVector2f;
+typedef EMatrix< int, 1, 2 > 	ERowVector2i;
+typedef EMatrix< std::complex< double >, 1, 3 > 	ERowVector3cd;
+typedef EMatrix< std::complex< float >, 1, 3 > 	ERowVector3cf;
+typedef EMatrix< double, 1, 3 > 	ERowVector3d;
+typedef EMatrix< float, 1, 3 > 	ERowVector3f;
+typedef EMatrix< int, 1, 3 > 	ERowVector3i;
+typedef EMatrix< std::complex< double >, 1, 4 > 	ERowVector4cd;
+typedef EMatrix< std::complex< float >, 1, 4 > 	ERowVector4cf;
+typedef EMatrix< double, 1, 4 > 	ERowVector4d;
+typedef EMatrix< float, 1, 4 > 	ERowVector4f;
+typedef EMatrix< int, 1, 4 > 	ERowVector4i;
+typedef EMatrix< std::complex< double >, 1, Eigen::Dynamic > 	ERowVectorXcd;
+typedef EMatrix< std::complex< float >, 1, Eigen::Dynamic > 	ERowVectorXcf;
+typedef EMatrix< double, 1, Eigen::Dynamic > 	ERowVectorXd;
+typedef EMatrix< float, 1, Eigen::Dynamic > 	ERowVectorXf;
+typedef EMatrix< int, 1, Eigen::Dynamic > 	ERowVectorXi;
+typedef EMatrix< std::complex< double >, 2, 1 > 	EVector2cd;
+typedef EMatrix< std::complex< float >, 2, 1 > 	EVector2cf;
+typedef EMatrix< double, 2, 1 > 	EVector2d;
+typedef EMatrix< float, 2, 1 > 	EVector2f;
+typedef EMatrix< int, 2, 1 > 	EVector2i;
+typedef EMatrix< std::complex< double >, 3, 1 > 	EVector3cd;
+typedef EMatrix< std::complex< float >, 3, 1 > 	EVector3cf;
+typedef EMatrix< double, 3, 1 > 	EVector3d;
+typedef EMatrix< float, 3, 1 > 	EVector3f;
+typedef EMatrix< int, 3, 1 > 	EVector3i;
+typedef EMatrix< std::complex< double >, 4, 1 > 	EVector4cd;
+typedef EMatrix< std::complex< float >, 4, 1 > 	EVector4cf;
+typedef EMatrix< double, 4, 1 > 	EVector4d;
+typedef EMatrix< float, 4, 1 > 	EVector4f;
+typedef EMatrix< int, 4, 1 > 	EVector4i;
+typedef EMatrix< std::complex< double >, Eigen::Dynamic, 1 > 	EVectorXcd;
+typedef EMatrix< std::complex< float >, Eigen::Dynamic, 1 > 	EVectorXcf;
+typedef EMatrix< double, Eigen::Dynamic, 1 > 	EVectorXd;
+typedef EMatrix< float, Eigen::Dynamic, 1 > 	EVectorXf;
+typedef EMatrix< int, Eigen::Dynamic, 1 > 	EVectorXi;
+
+#else
+
+// There is not EMatrix if we do not have Eigen
+
+#endif
+
+#endif /* OPENFPM_DATA_SRC_SPACE_EMATRIX_HPP_ */
diff --git a/src/DMatrix/tests/EMatrix_unit_tests.cpp b/src/DMatrix/tests/EMatrix_unit_tests.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..741803b90b0492634d13a086325ba481e58a2c88
--- /dev/null
+++ b/src/DMatrix/tests/EMatrix_unit_tests.cpp
@@ -0,0 +1,108 @@
+/*
+ * EMatrix_unit_tests.cpp
+ *
+ *  Created on: Feb 13, 2018
+ *      Author: i-bird
+ */
+
+#define BOOST_TEST_DYN_LINK
+#include <boost/test/unit_test.hpp>
+#include "DMatrix/EMatrix.hpp"
+#include "memory/HeapMemory.hpp"
+
+#ifdef HAVE_EIGEN
+
+BOOST_AUTO_TEST_SUITE (EMatrix_test)
+
+BOOST_AUTO_TEST_CASE( EMatrix_test_use)
+{
+	{
+	EMatrixXd em;
+
+	em.resize(8,5);
+
+	for (size_t i = 0 ; i < 8 ; i++)
+	{
+		for (size_t j = 0 ; j < 5 ; j++)
+		{em(i,j) = i*8+j;}
+	}
+
+	size_t pr = 0;
+	em.packRequest(pr);
+
+	// allocate the memory
+	HeapMemory pmem;
+	pmem.allocate(pr);
+	ExtPreAlloc<HeapMemory> & mem = *(new ExtPreAlloc<HeapMemory>(pr,pmem));
+	mem.incRef();
+
+	BOOST_REQUIRE_EQUAL(pr,8*5*sizeof(double) + 2*sizeof(size_t));
+
+	Pack_stat sts;
+	em.pack(mem,sts);
+
+	// Reset to zero
+	for (size_t i = 0 ; i < 8 ; i++)
+	{
+		for (size_t j = 0 ; j < 5 ; j++)
+		{em(i,j) = 0;}
+	}
+
+	Unpack_stat ps;
+	em.unpack(mem,ps);
+
+	for (size_t i = 0 ; i < 8 ; i++)
+	{
+		for (size_t j = 0 ; j < 5 ; j++)
+		{BOOST_REQUIRE_EQUAL(em(i,j),i*8+j);}
+	}
+	}
+
+
+	{
+	EMatrix3d em;
+
+	em.resize(3,3);
+
+	for (size_t i = 0 ; i < 3 ; i++)
+	{
+		for (size_t j = 0 ; j < 3 ; j++)
+		{em(i,j) = i*8+j;}
+	}
+
+	size_t pr = 0;
+	em.packRequest(pr);
+
+	// allocate the memory
+	HeapMemory pmem;
+	pmem.allocate(pr);
+	ExtPreAlloc<HeapMemory> & mem = *(new ExtPreAlloc<HeapMemory>(pr,pmem));
+	mem.incRef();
+
+	BOOST_REQUIRE_EQUAL(pr,3*3*sizeof(double) + 2*sizeof(size_t));
+
+	Pack_stat sts;
+	em.pack(mem,sts);
+
+	// Reset to zero
+	for (size_t i = 0 ; i < 3 ; i++)
+	{
+		for (size_t j = 0 ; j < 3 ; j++)
+		{em(i,j) = 0;}
+	}
+
+	Unpack_stat ps;
+	em.unpack(mem,ps);
+
+	for (size_t i = 0 ; i < 3 ; i++)
+	{
+		for (size_t j = 0 ; j < 3 ; j++)
+		{BOOST_REQUIRE_EQUAL(em(i,j),i*8+j);}
+	}
+	}
+}
+
+BOOST_AUTO_TEST_SUITE_END()
+
+
+#endif
diff --git a/src/Draw/PointIteratorSkin.hpp b/src/Draw/PointIteratorSkin.hpp
index fc5e7fd81071107a71a05eb684ed8060af1a5b68..6fded693e06e33d29eb6e0ed87782b21fd4b4899 100644
--- a/src/Draw/PointIteratorSkin.hpp
+++ b/src/Draw/PointIteratorSkin.hpp
@@ -136,8 +136,15 @@ public:
 	 * \param sp grid spacing
 	 *
 	 */
-	PointIteratorSkin( Decomposition & dec, size_t (& sz)[dim], const Box<dim,T> & domain, const Box<dim,T> & sub_A, const Box<dim,T> & sub_B, size_t (& bc)[dim])
-	:grid_dist_id_iterator_dec_skin<Decomposition>(dec, sz, getAB(sz,domain,sub_A,sub_B,sp,RETURN_A), getAB(sz,domain,sub_A,sub_B,sp,RETURN_B), bc),domain(domain)
+	PointIteratorSkin( Decomposition & dec,
+			           size_t (& sz)[dim],
+					   const Box<dim,T> & domain,
+					   const Box<dim,T> & sub_A,
+					   const Box<dim,T> & sub_B,
+					   size_t (& bc)[dim])
+	:grid_dist_id_iterator_dec_skin<Decomposition>(dec, sz, getAB(sz,domain,sub_A,sub_B,sp,RETURN_A),
+	 getAB(sz,domain,sub_A,sub_B,sp,RETURN_B), bc),
+	 domain(domain)
 	{
 		sub_domainA.add(sub_A);
 		calculateAp();
diff --git a/src/FiniteDifference/FDScheme.hpp b/src/FiniteDifference/FDScheme.hpp
index 67db9cfde58ee7c0e667bf4ee8c3b5417488bd06..d86a4acf796adfff67508877118b0711dfe35b55 100644
--- a/src/FiniteDifference/FDScheme.hpp
+++ b/src/FiniteDifference/FDScheme.hpp
@@ -441,6 +441,46 @@ private:
 		}
 	}
 
+	/*! \brief Impose an operator
+	 *
+	 * This function the RHS no matrix is imposed. This
+	 * function is usefull if the Matrix has been constructed and only
+	 * the right hand side b must be changed
+	 *
+	 * Ax = b
+	 *
+	 * \param num right hand side of the term (b term)
+	 * \param id Equation id in the system that we are imposing
+	 * \param it_d iterator that define where you want to impose
+	 *
+	 */
+	template<typename bop, typename iterator>
+	void impose_dit_b(bop num,
+					  long int id ,
+					  const iterator & it_d)
+	{
+
+		auto it = it_d;
+		grid_sm<Sys_eqs::dims,void> gs = g_map.getGridInfoVoid();
+
+		// iterate all the grid points
+		while (it.isNext())
+		{
+			// get the position
+			auto key = it.get();
+
+			b(g_map.template get<0>(key)*Sys_eqs::nvar + id) = num.get(key);
+
+			// if SE_CLASS1 is defined check the position
+#ifdef SE_CLASS1
+//			T::position(key,gs,s_pos);
+#endif
+
+			++row_b;
+			++it;
+		}
+	}
+
 	/*! \brief Impose an operator
 	 *
 	 * This function impose an operator on a particular grid region to produce the system
@@ -456,7 +496,7 @@ private:
 	 * \param it_d iterator that define where you want to impose
 	 *
 	 */
-	template<typename T, typename bop, typename iterator> void impose_dit(const T & op ,
+	template<typename T, typename bop, typename iterator> void impose_dit(const T & op,
 			                         bop num,
 									 long int id ,
 									 const iterator & it_d)
@@ -493,9 +533,8 @@ private:
 				trpl.last().value() = it->second;
 
 				if (trpl.last().row() == trpl.last().col())
-					is_diag = true;
+				{is_diag = true;}
 
-//				std::cout << "(" << trpl.last().row() << "," << trpl.last().col() << "," << trpl.last().value() << ")" << "\n";
 			}
 
 			// If does not have a diagonal entry put it to zero
@@ -812,6 +851,44 @@ public:
 
 	}
 
+	/*! \brief In case we want to impose a new b re-using FDScheme we have to call
+	 *         This function
+	 */
+	void new_b()
+	{row_b = 0;}
+
+	/*! \brief In case we want to impose a new A re-using FDScheme we have to call
+	 *         This function
+	 *
+	 */
+	void new_A()
+	{row = 0;}
+
+	/*! \brief Impose an operator
+	 *
+	 * This function impose an operator on a particular grid region to produce the system
+	 *
+	 * Ax = b
+	 *
+	 * ## Stokes equation 2D, lid driven cavity with one splipping wall
+	 * \snippet eq_unit_test.hpp Copy the solution to grid
+	 *
+	 * \param op Operator to impose (A term)
+	 * \param num right hand side of the term (b term)
+	 * \param id Equation id in the system that we are imposing
+	 * \param it_d iterator that define where you want to impose
+	 *
+	 */
+	template<unsigned int prp, typename b_term, typename iterator>
+	void impose_dit_b(b_term & b_t,
+					  const iterator & it_d,
+					  long int id = 0)
+	{
+		grid_b<b_term,prp> b(b_t);
+
+		impose_dit_b(b,id,it_d);
+	}
+
 	/*! \brief Impose an operator
 	 *
 	 * This function impose an operator on a particular grid region to produce the system
diff --git a/src/FiniteDifference/eq_unit_test_3d.cpp b/src/FiniteDifference/eq_unit_test_3d.cpp
index 1d59be60f227a065ed91e75764a3270ebce23df9..7ab1062b6367610cac7364b7e8d074cfece479f8 100644
--- a/src/FiniteDifference/eq_unit_test_3d.cpp
+++ b/src/FiniteDifference/eq_unit_test_3d.cpp
@@ -229,6 +229,11 @@ template<typename solver_type,typename lid_nn_3d> void lid_driven_cavity_3d()
         std::string file1 = std::string("test/") + s + "lid_driven_cavity_3d_p" + std::to_string(v_cl.getProcessingUnits()) + "_grid_" + std::to_string(v_cl.getProcessUnitID()) + "_test_GCC6.vtk";
         std::string file2 = s + "lid_driven_cavity_3d_p" + std::to_string(v_cl.getProcessingUnits()) + "_grid_" + std::to_string(v_cl.getProcessUnitID()) + ".vtk";
 
+    #elif __GNUC__ == 7
+
+        std::string file1 = std::string("test/") + s + "lid_driven_cavity_3d_p" + std::to_string(v_cl.getProcessingUnits()) + "_grid_" + std::to_string(v_cl.getProcessUnitID()) + "_test_GCC7.vtk";
+        std::string file2 = s + "lid_driven_cavity_3d_p" + std::to_string(v_cl.getProcessingUnits()) + "_grid_" + std::to_string(v_cl.getProcessUnitID()) + ".vtk";
+
 	#else
 
         std::string file1 = std::string("test/") + s + "lid_driven_cavity_3d_p" + std::to_string(v_cl.getProcessingUnits()) + "_grid_" + std::to_string(v_cl.getProcessUnitID()) + "_test_GCC4.vtk";
@@ -254,7 +259,7 @@ template<typename solver_type,typename lid_nn_3d> void lid_driven_cavity_3d()
 
 BOOST_AUTO_TEST_CASE(lid_driven_cavity)
 {
-#ifdef HAVE_EIGEN
+#if defined(HAVE_EIGEN) && defined(HAVE_SUITESPARSE)
 	lid_driven_cavity_3d<umfpack_solver<double>,lid_nn_3d_eigen>();
 #endif
 #ifdef HAVE_PETSC
diff --git a/src/Makefile.am b/src/Makefile.am
index daaf7d879a0e47f8ffc7411e29911abbe9a5e46d..ff728f2fa44668fd49e744ad1977f0216469ae80 100755
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -2,8 +2,13 @@
 LINKLIBS = $(HDF5_LDFLAGS)  $(HDF5_LIBS) $(OPENMP_LDFLAGS) $(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) $(LIBQUADMATH) $(OPENMP_LDFLAGS) $(LIBIFCORE)
 
 noinst_PROGRAMS = numerics
+<<<<<<< HEAD
 numerics_SOURCES = main.cpp Matrix/SparseMatrix_unit_tests.cpp interpolation/interpolation_unit_tests.cpp  Vector/Vector_unit_tests.cpp  Solvers/petsc_solver_unit_tests.cpp FiniteDifference/FDScheme_unit_tests.cpp FiniteDifference/eq_unit_test_3d.cpp FiniteDifference/eq_unit_test.cpp  Operators/Vector/vector_dist_operators_unit_tests.cpp ../../openfpm_vcluster/src/VCluster/VCluster.cpp ../../openfpm_devices/src/memory/HeapMemory.cpp ../../openfpm_devices/src/memory/PtrMemory.cpp ../../openfpm_devices/src/Memleak_check.cpp
 numerics_CXXFLAGS = -Wno-unknown-pragmas -Wno-int-in-bool-context $(HDF5_CPPFLAGS) $(OPENMP_CFLAGS) $(LIBHILBERT_INCLUDE) $(AM_CXXFLAGS) $(INCLUDES_PATH) $(BOOST_CPPFLAGS) $(SUITESPARSE_INCLUDE) $(METIS_INCLUDE) $(PARMETIS_INCLUDE) $(EIGEN_INCLUDE) $(PETSC_INCLUDE) -Wno-deprecated-declarations -Wno-unused-local-typedefs
+=======
+numerics_SOURCES = main.cpp Solvers/petsc_solver_unit_tests.cpp DMatrix/tests/EMatrix_unit_tests.cpp  ../../openfpm_vcluster/src/VCluster/VCluster.cpp ../../openfpm_devices/src/memory/HeapMemory.cpp ../../openfpm_devices/src/memory/PtrMemory.cpp ../../openfpm_devices/src/Memleak_check.cpp
+numerics_CXXFLAGS = $(HDF5_CPPFLAGS) $(OPENMP_CFLAGS) $(LIBHILBERT_INCLUDE) $(AM_CXXFLAGS) $(INCLUDES_PATH) $(BOOST_CPPFLAGS) $(SUITESPARSE_INCLUDE) $(METIS_INCLUDE) $(PARMETIS_INCLUDE) $(EIGEN_INCLUDE) $(PETSC_INCLUDE) -Wno-deprecated-declarations -Wno-unused-local-typedefs
+>>>>>>> 0b425037b139f9c2dbe989da5468eebb215bc623
 numerics_CFLAGS = $(CUDA_CFLAGS)
 numerics_LDADD = $(LINKLIBS) -lparmetis -lmetis
 nobase_include_HEADERS = Matrix/SparseMatrix.hpp Matrix/SparseMatrix_Eigen.hpp Matrix/SparseMatrix_petsc.hpp \
@@ -13,6 +18,7 @@ util/petsc_util.hpp util/linalgebra_lib.hpp util/util_num.hpp util/grid_dist_tes
 FiniteDifference/Average.hpp FiniteDifference/Derivative.hpp FiniteDifference/FD_util_include.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 \
 Draw/DrawParticles.hpp Draw/PointIterator.hpp Draw/PointIteratorSkin.hpp \
+DMatrix/EMatrix.hpp \
 interpolation/interpolation.hpp interpolation/mp4_kernel.hpp interpolation/z_spline.hpp Solvers/petsc_solver_AMG_report.hpp
 
 .cu.o :
diff --git a/src/Operators/Vector/vector_dist_operators_apply_kernel.hpp b/src/Operators/Vector/vector_dist_operators_apply_kernel.hpp
index 844360ab0ce60901e893c3606365b2116beefffd..80ed67c6c76d5a9e673c25e51efa0ede2f8d6f3d 100644
--- a/src/Operators/Vector/vector_dist_operators_apply_kernel.hpp
+++ b/src/Operators/Vector/vector_dist_operators_apply_kernel.hpp
@@ -8,6 +8,16 @@
 #ifndef OPENFPM_NUMERICS_SRC_OPERATORS_VECTOR_VECTOR_DIST_OPERATORS_APPLY_KERNEL_HPP_
 #define OPENFPM_NUMERICS_SRC_OPERATORS_VECTOR_VECTOR_DIST_OPERATORS_APPLY_KERNEL_HPP_
 
+//////// SET of small macro to make to define integral easy
+
+#define DEFINE_INTERACTION_3D(name) struct name \
+{\
+\
+	Point<3,double> value(const Point<3,double> & xp, const Point<3,double> xq)\
+	{
+
+#define END_INTERACTION }\
+						};
 
 /*! \brief is_expression check if a type is simple a type or is just an encapsulation of an expression
  *
diff --git a/src/Operators/Vector/vector_dist_operators_functions.hpp b/src/Operators/Vector/vector_dist_operators_functions.hpp
index 50e397b48a22ead8679acae99632d8c3132c8707..6f8fe79e34d3d005a0bef52d129e2ecfb87d504b 100644
--- a/src/Operators/Vector/vector_dist_operators_functions.hpp
+++ b/src/Operators/Vector/vector_dist_operators_functions.hpp
@@ -204,7 +204,6 @@ fun_name(double d, const vector_dist_expression_op<exp1,exp2,op1> & va)\
 
 CREATE_VDIST_ARG2_FUNC(pmul,pmul,VECT_PMUL)
 
-
 ////////// Special function reduce /////////////////////////
 
 
@@ -291,5 +290,17 @@ rsum(const vector_dist_expression<prp1,v1> & va, const vector_type & vd)
 	return exp_sum;
 }
 
+namespace openfpm
+{
+	/*! \brief General distance formula
+	 *
+	 *
+	 */
+	template <typename T, typename P> auto distance(T exp1, P exp2) -> decltype(norm(exp1 - exp2))
+	{
+		return norm(exp1 - exp2);
+	}
+}
+
 
 #endif /* OPENFPM_NUMERICS_SRC_OPERATORS_VECTOR_VECTOR_DIST_OPERATORS_FUNCTIONS_HPP_ */
diff --git a/src/Operators/Vector/vector_dist_operators_unit_tests.cpp b/src/Operators/Vector/vector_dist_operators_unit_tests.cpp
index 2d6f2a958f620f1ff8f2723f9629f29730d903ab..4215733d60b81ecf1193e5c4ac98af05e8f7dfa6 100644
--- a/src/Operators/Vector/vector_dist_operators_unit_tests.cpp
+++ b/src/Operators/Vector/vector_dist_operators_unit_tests.cpp
@@ -940,7 +940,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_operators_test )
 
 	// normalization function
 
-	vA = vVB * vVC + norm(vVC + vVB) + distance(vVC,vVB);
+	vA = vVB * vVC + norm(vVC + vVB) + openfpm::distance(vVC,vVB);
 	check_values_scal_norm_dist(vd);
 
 	Point<3,float> p0({2.0,2.0,2.0});
diff --git a/src/PSE/Kernels_test_util.hpp b/src/PSE/Kernels_test_util.hpp
index f74682b718bcf17a3a279bd37ee362edc748be46..c179f6b4411d0e4c0aeaf0b1fb9341ebfbf4c4fa 100644
--- a/src/PSE/Kernels_test_util.hpp
+++ b/src/PSE/Kernels_test_util.hpp
@@ -86,7 +86,7 @@ template<typename T, typename Kernel> void PSE_test(size_t Npart, size_t overlap
     size_t bc[1]={NON_PERIODIC};
 	Ghost<1,T> g(20*eps);
 
-	vector_dist<1,T, aggregate<T>, CartDecomposition<1,T> > vd(Npart,box,bc,g);
+	vector_dist<1,T, aggregate<T> > vd(Npart,box,bc,g);
 
 	auto it2 = vd.getIterator();
 
diff --git a/src/Solvers/umfpack_solver.hpp b/src/Solvers/umfpack_solver.hpp
index 2a8287b6047e2df4526c4f2c4aa52e3bf3ff188f..b38b152a74cd41c13d8311bb28acc9a05a13d77d 100644
--- a/src/Solvers/umfpack_solver.hpp
+++ b/src/Solvers/umfpack_solver.hpp
@@ -14,7 +14,7 @@
 #define SOLVER_PRINT_RESIDUAL_NORM_INFINITY 1
 #define SOLVER_PRINT_DETERMINANT 2
 
-#ifdef HAVE_EIGEN
+#if defined(HAVE_EIGEN) && defined(HAVE_SUITESPARSE)
 
 /////// Compiled with EIGEN support
 
diff --git a/src/Vector/Vector_eigen.hpp b/src/Vector/Vector_eigen.hpp
index b20d5162ee6969e61135d474e386aed60d8cfe97..3af417f45fd232389d1249e82bf215430841ce9e 100644
--- a/src/Vector/Vector_eigen.hpp
+++ b/src/Vector/Vector_eigen.hpp
@@ -75,6 +75,16 @@ public:
 		row() = i;
 		value() = val;
 	}
+
+	/*! \brief Indicate that the structure has no pointer
+	 *
+	 * \return true
+	 *
+	 */
+	static inline bool noPointers()
+	{
+		return true;
+	}
 };
 
 template<typename T>
diff --git a/src/Vector/Vector_petsc.hpp b/src/Vector/Vector_petsc.hpp
index 23a5e4919c4fe353a7108f703c15a00ce2089d5a..607e205bf712f37ad25d6adc89e764e77967926d 100644
--- a/src/Vector/Vector_petsc.hpp
+++ b/src/Vector/Vector_petsc.hpp
@@ -79,6 +79,16 @@ public:
 		rw() = i;
 		val() = val;
 	}
+
+	/*! \brief Indicate that the structure has no pointer
+	 *
+	 * \return true
+	 *
+	 */
+	static inline bool noPointers()
+	{
+		return true;
+	}
 };
 
 constexpr unsigned int row_id = 0;
diff --git a/src/interpolation/interpolation.hpp b/src/interpolation/interpolation.hpp
index 8fd0ca06827ae68409fa40374c3a055cb5355518..28d14ad98626d0130b83597dbb7b69ea04bad7fc 100644
--- a/src/interpolation/interpolation.hpp
+++ b/src/interpolation/interpolation.hpp
@@ -11,6 +11,7 @@
 #include "NN/Mem_type/MemFast.hpp"
 #include "NN/CellList/CellList.hpp"
 #include "Grid/grid_dist_key.hpp"
+#include "Vector/vector_dist_key.hpp"
 
 #define INTERPOLATION_ERROR_OBJECT std::runtime_error("Runtime interpolation error");
 
@@ -407,8 +408,8 @@ struct inte_calc_impl
 	 *        kernel stencil for each local grid (sub-domain)
 	 *
 	 */
-	template<unsigned int prp_g, unsigned int prp_v, unsigned int m2p_or_p2m, unsigned int np_a_int, typename iterator, typename grid>
-																	 static inline void inte_calc(iterator & it,
+	template<unsigned int prp_g, unsigned int prp_v, unsigned int m2p_or_p2m, unsigned int np_a_int, typename grid>
+																	 static inline void inte_calc(const vect_dist_key_dx & key_p,
 																	 vector & vd,
 																	 const Box<vector::dims,typename vector::stype> & domain,
 																	 int (& ip)[vector::dims][kernel::np],
@@ -422,8 +423,6 @@ struct inte_calc_impl
 																	 const CellList<vector::dims,typename vector::stype,Mem_fast<>,shift<vector::dims,typename vector::stype>> & geo_cell,
 																	 openfpm::vector<agg_arr<openfpm::math::pow(kernel::np,vector::dims)>> & offsets)
 	{
-		auto key_p = it.get();
-
 		Point<vector::dims,typename vector::stype> p = vd.getPos(key_p);
 
 		// On which sub-domain we interpolate the particle
@@ -434,24 +433,24 @@ struct inte_calc_impl
 		// calculate the position of the particle in the grid
 		// coordinated
 		for (size_t i = 0 ; i < vector::dims ; i++)
-			x0[i] = (p.get(i)-domain.getLow(i))*dx[i];
+		{x0[i] = (p.get(i)-domain.getLow(i))*dx[i];}
 
 		// convert into integer
 		for (size_t i = 0 ; i < vector::dims ; i++)
-			ip[i][0] = (int)x0[i];
+		{ip[i][0] = (int)x0[i];}
 
 		// convert the global grid position into local grid position
 		grid_key_dx<vector::dims> base;
 
 		for (size_t i = 0 ; i < vector::dims ; i++)
-			base.set_d(i,ip[i][0] - gd.getLocalGridsInfo().get(sub).origin.get(i) - (long int)kernel::np/2 + 1);
+		{base.set_d(i,ip[i][0] - gd.getLocalGridsInfo().get(sub).origin.get(i) - (long int)kernel::np/2 + 1);}
 
 		// convenient grid of points
 
 		for (size_t j = 0 ; j < kernel::np-1 ; j++)
 		{
 			for (size_t i = 0 ; i < vector::dims ; i++)
-				ip[i][j+1] = (int)ip[i][j]+1;
+			{ip[i][j+1] = (int)ip[i][j]+1;}
 		}
 
 		for (size_t i = 0 ; i < vector::dims ; i++)
@@ -460,13 +459,13 @@ struct inte_calc_impl
 		for (long int j = 0 ; j < kernel::np ; j++)
 		{
 			for (size_t i = 0 ; i < vector::dims ; i++)
-				x[i][j] = - xp[i] + typename vector::stype((long int)j - (long int)kernel::np/2 + 1);
+			{x[i][j] = - xp[i] + typename vector::stype((long int)j - (long int)kernel::np/2 + 1);}
 		}
 
 		for (size_t j = 0 ; j < kernel::np ; j++)
 		{
 			for (size_t i = 0 ; i < vector::dims ; i++)
-				a[i][j] = kernel::value(x[i][j],j);
+			{a[i][j] = kernel::value(x[i][j],j);}
 		}
 
 		calculate_aint<vector::dims,vector,kernel::np>::value(sz,a_int,a);
@@ -541,6 +540,18 @@ class interpolate
 	//! Type of the calculations
 	typedef typename vector::stype arr_type;
 
+	//! the offset for each sub-sub-domain
+	openfpm::vector<agg_arr<openfpm::math::pow(kernel::np,vector::dims)>> offsets;
+
+	//! kernel size
+	size_t sz[vector::dims];
+
+	//! grid spacing
+	typename vector::stype dx[vector::dims];
+
+	//! Simulation domain
+	Box<vector::dims,typename vector::stype> domain;
+
 	/*! \brief It calculate the interpolation stencil offsets
 	 *
 	 * \param offsets array where to store the linearized offset of the
@@ -624,6 +635,18 @@ public:
 				++g_sub;
 			}
 		}
+
+		for (size_t i = 0 ; i < vector::dims ; i++)
+		{sz[i] = kernel::np;}
+
+		calculate_the_offsets(offsets,sz);
+
+		// calculate spacing
+		for (size_t i = 0 ; i < vector::dims ; i++)
+		{dx[i] = 1.0/gd.spacing(i);}
+
+		// simulation domain
+		domain = vd.getDecomposition().getDomain();
 	};
 
 	/*! \brief Interpolate particles to mesh
@@ -651,23 +674,51 @@ public:
 
 #endif
 
-		Box<vector::dims,typename vector::stype> domain = vd.getDecomposition().getDomain();
+		// point position
+		typename vector::stype xp[vector::dims];
 
-		// grid spacing
-		typename vector::stype dx[vector::dims];
+		int ip[vector::dims][kernel::np];
+		typename vector::stype x[vector::dims][kernel::np];
+		typename vector::stype a[vector::dims][kernel::np];
 
-		for (size_t i = 0 ; i < vector::dims ; i++)
-			dx[i] = 1.0/gd.spacing(i);
+		typename vector::stype a_int[openfpm::math::pow(kernel::np,vector::dims)];
 
-		size_t sz[vector::dims];
+		auto it = vd.getDomainIterator();
 
-		for (size_t i = 0 ; i < vector::dims ; i++)
-			sz[i] = kernel::np;
+		while (it.isNext() == true)
+		{
+			auto key_p = it.get();
 
-		// Precalculate the offset for each sub-sub-domain
-		openfpm::vector<agg_arr<openfpm::math::pow(kernel::np,vector::dims)>> offsets;
+			inte_calc_impl<vector,kernel>::template inte_calc<prp_g,prp_v,inte_p2m,openfpm::math::pow(kernel::np,vector::dims)>(key_p,vd,domain,ip,gd,dx,xp,a_int,a,x,sz,geo_cell,offsets);
 
-		calculate_the_offsets(offsets,sz);
+			++it;
+		}
+	}
+
+	/*! \brief Interpolate mesh to particle
+	 *
+	 * Most of the time the particle set and the mesh are the same
+	 * as the one used in the constructor. They also can be different
+	 * as soon as they have the same decomposition
+	 *
+	 * \param gd grid or mesh
+	 * \param vd particle set
+	 *
+	 */
+	template<unsigned int prp_g, unsigned int prp_v> void m2p(grid & gd, vector & vd)
+	{
+#ifdef SE_CLASS1
+
+		if (!vd.getDecomposition().is_equal_ng(gd.getDecomposition()) )
+		{
+			std::cerr << __FILE__ << ":" << __LINE__ << " Error: the distribution of the vector of particles" <<
+					" and the grid is different. In order to interpolate the two data structure must have the" <<
+					" same decomposition" << std::endl;
+
+			ACTION_ON_ERROR(INTERPOLATION_ERROR_OBJECT)
+		}
+
+#endif
 
 		// point position
 		typename vector::stype xp[vector::dims];
@@ -676,19 +727,24 @@ public:
 		typename vector::stype x[vector::dims][kernel::np];
 		typename vector::stype a[vector::dims][kernel::np];
 
+		//		grid_cpu<vector::dims,aggregate<typename vector::stype>> a_int(sz);
+		//		a_int.setMemory();
 		typename vector::stype a_int[openfpm::math::pow(kernel::np,vector::dims)];
 
 		auto it = vd.getDomainIterator();
 
 		while (it.isNext() == true)
 		{
-			inte_calc_impl<vector,kernel>::template inte_calc<prp_g,prp_v,inte_p2m,openfpm::math::pow(kernel::np,vector::dims)>(it,vd,domain,ip,gd,dx,xp,a_int,a,x,sz,geo_cell,offsets);
+			auto key_p = it.get();
+
+			inte_calc_impl<vector,kernel>::template inte_calc<prp_g,prp_v,inte_m2p,openfpm::math::pow(kernel::np,vector::dims)>(key_p,vd,domain,ip,gd,dx,xp,a_int,a,x,sz,geo_cell,offsets);
 
 			++it;
 		}
 	}
 
-	/*! \brief Interpolate mesh to particle
+
+	/*! \brief Interpolate particles to mesh
 	 *
 	 * Most of the time the particle set and the mesh are the same
 	 * as the one used in the constructor. They also can be different
@@ -696,9 +752,10 @@ public:
 	 *
 	 * \param gd grid or mesh
 	 * \param vd particle set
+	 * \param p particle
 	 *
 	 */
-	template<unsigned int prp_g, unsigned int prp_v> void m2p(grid & gd, vector & vd)
+	template<unsigned int prp_v, unsigned int prp_g> inline void p2m(vector & vd, grid & gd,const vect_dist_key_dx & p)
 	{
 #ifdef SE_CLASS1
 
@@ -713,14 +770,6 @@ public:
 
 #endif
 
-		Box<vector::dims,typename vector::stype> domain = vd.getDecomposition().getDomain();
-
-		// grid spacing
-		typename vector::stype dx[vector::dims];
-
-		for (size_t i = 0 ; i < vector::dims ; i++)
-			dx[i] = 1.0/gd.spacing(i);
-
 		// point position
 		typename vector::stype xp[vector::dims];
 
@@ -728,28 +777,49 @@ public:
 		typename vector::stype x[vector::dims][kernel::np];
 		typename vector::stype a[vector::dims][kernel::np];
 
-		size_t sz[vector::dims];
+		typename vector::stype a_int[openfpm::math::pow(kernel::np,vector::dims)];
 
-		for (size_t i = 0 ; i < vector::dims ; i++)
-			sz[i] = kernel::np;
+		inte_calc_impl<vector,kernel>::template inte_calc<prp_g,prp_v,inte_p2m,openfpm::math::pow(kernel::np,vector::dims)>(p,vd,domain,ip,gd,dx,xp,a_int,a,x,sz,geo_cell,offsets);
+	}
 
-		// Precalculate the offset for each sub-sub-domain
-		openfpm::vector<agg_arr<openfpm::math::pow(kernel::np,vector::dims)>> offsets;
+	/*! \brief Interpolate mesh to particle
+	 *
+	 * Most of the time the particle set and the mesh are the same
+	 * as the one used in the constructor. They also can be different
+	 * as soon as they have the same decomposition
+	 *
+	 * \param gd grid or mesh
+	 * \param vd particle set
+	 * \param p particle
+	 *
+	 */
+	template<unsigned int prp_g, unsigned int prp_v> inline void m2p(grid & gd, vector & vd, const vect_dist_key_dx & p)
+	{
+#ifdef SE_CLASS1
 
-		calculate_the_offsets(offsets,sz);
+		if (!vd.getDecomposition().is_equal_ng(gd.getDecomposition()) )
+		{
+			std::cerr << __FILE__ << ":" << __LINE__ << " Error: the distribution of the vector of particles" <<
+					" and the grid is different. In order to interpolate the two data structure must have the" <<
+					" same decomposition" << std::endl;
+
+			ACTION_ON_ERROR(INTERPOLATION_ERROR_OBJECT)
+		}
+
+#endif
+
+		// point position
+		typename vector::stype xp[vector::dims];
+
+		int ip[vector::dims][kernel::np];
+		typename vector::stype x[vector::dims][kernel::np];
+		typename vector::stype a[vector::dims][kernel::np];
 
 		//		grid_cpu<vector::dims,aggregate<typename vector::stype>> a_int(sz);
 		//		a_int.setMemory();
 		typename vector::stype a_int[openfpm::math::pow(kernel::np,vector::dims)];
 
-		auto it = vd.getDomainIterator();
-
-		while (it.isNext() == true)
-		{
-			inte_calc_impl<vector,kernel>::template inte_calc<prp_g,prp_v,inte_m2p,openfpm::math::pow(kernel::np,vector::dims)>(it,vd,domain,ip,gd,dx,xp,a_int,a,x,sz,geo_cell,offsets);
-
-			++it;
-		}
+		inte_calc_impl<vector,kernel>::template inte_calc<prp_g,prp_v,inte_m2p,openfpm::math::pow(kernel::np,vector::dims)>(p,vd,domain,ip,gd,dx,xp,a_int,a,x,sz,geo_cell,offsets);
 	}
 
 };
diff --git a/src/interpolation/interpolation_unit_tests.cpp b/src/interpolation/interpolation_unit_tests.cpp
index 9512a73fe31ea2539e418a789002d5c6654989f3..ed795c1d0f794d016b39b418c255c532ee36b043 100644
--- a/src/interpolation/interpolation_unit_tests.cpp
+++ b/src/interpolation/interpolation_unit_tests.cpp
@@ -88,8 +88,63 @@ template<typename vector, unsigned int mom_p> void momenta_vector(vector & vd,ty
 	}
 }
 
+template<unsigned int dim, typename T, typename grid, typename vector>
+void interp_test(grid & gd, vector & vd, bool single_particle)
+{
+	// Reset the grid
 
-BOOST_AUTO_TEST_CASE( interpolation_full_test_2D )
+	auto it2 = gd.getDomainGhostIterator();
+
+	while (it2.isNext())
+	{
+		auto key = it2.get();
+
+		gd.template get<0>(key) = 0.0;
+
+		++it2;
+	}
+
+	interpolate<vector,grid,mp4_kernel<float>> inte(vd,gd);
+
+	if (single_particle == false)
+	{inte.template p2m<0,0>(vd,gd);}
+	else
+	{
+		auto it = vd.getDomainIterator();
+
+		while (it.isNext())
+		{
+			auto p = it.get();
+
+			inte.template p2m<0,0>(vd,gd,p);
+
+			++it;
+		}
+	}
+
+	T mg[dim];
+	T mv[dim];
+
+	momenta_grid<grid,0>(gd,mg);
+	momenta_vector<vector,0>(vd,mv);
+
+	for (size_t i = 0 ; i < dim ; i++)
+	{BOOST_REQUIRE_CLOSE(mg[i],mv[i],0.001);}
+
+	momenta_grid<grid,1>(gd,mg);
+	momenta_vector<vector,1>(vd,mv);
+
+	for (size_t i = 0 ; i < dim ; i++)
+	{BOOST_REQUIRE_CLOSE(mg[i],mv[i],0.001);}
+
+	momenta_grid<grid,2>(gd,mg);
+	momenta_vector<vector,2>(vd,mv);
+
+	for (size_t i = 0 ; i < dim ; i++)
+	{BOOST_REQUIRE_CLOSE(mg[i],mv[i],0.001);}
+}
+
+BOOST_AUTO_TEST_CASE( interpolation_full_single_test_2D )
 {
 	Box<2,float> domain({0.0,0.0},{1.0,1.0});
 	size_t sz[2] = {64,64};
@@ -99,8 +154,7 @@ BOOST_AUTO_TEST_CASE( interpolation_full_test_2D )
 
 	size_t bc_v[2] = {PERIODIC,PERIODIC};
 
-	{
-	vector_dist<2,float,aggregate<float>> vd(4096,domain,bc_v,gv);
+	vector_dist<2,float,aggregate<float>> vd(65536,domain,bc_v,gv);
 	grid_dist_id<2,float,aggregate<float>> gd(vd.getDecomposition(),sz,gg);
 
 	// set one particle on vd
@@ -121,46 +175,52 @@ BOOST_AUTO_TEST_CASE( interpolation_full_test_2D )
 
 	vd.map();
 
-	// Reset the grid
+	interp_test<2,float>(gd,vd,true);
+}
 
-	auto it2 = gd.getDomainGhostIterator();
 
-	while (it2.isNext())
-	{
-		auto key = it2.get();
+BOOST_AUTO_TEST_CASE( interpolation_full_test_2D )
+{
+	Box<2,float> domain({0.0,0.0},{1.0,1.0});
+	size_t sz[2] = {64,64};
 
-		gd.template get<0>(key) = 0.0;
+	Ghost<2,long int> gg(2);
+	Ghost<2,float> gv(0.01);
 
-		++it2;
-	}
+	size_t bc_v[2] = {PERIODIC,PERIODIC};
 
-	interpolate<decltype(vd),decltype(gd),mp4_kernel<float>> inte(vd,gd);
+	{
+	vector_dist<2,float,aggregate<float>> vd(4096,domain,bc_v,gv);
+	grid_dist_id<2,float,aggregate<float>> gd(vd.getDecomposition(),sz,gg);
 
-	inte.p2m<0,0>(vd,gd);
+	// set one particle on vd
 
-	float mg[2];
-	float mv[2];
+	auto it = vd.getDomainIterator();
 
-	momenta_grid<decltype(gd),0>(gd,mg);
-	momenta_vector<decltype(vd),0>(vd,mv);
+	while (it.isNext())
+	{
+		auto p = it.get();
 
-	BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
-	BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
+		vd.getPos(p)[0] = (double)rand()/RAND_MAX;
+		vd.getPos(p)[1] = (double)rand()/RAND_MAX;
 
-	momenta_grid<decltype(gd),1>(gd,mg);
-	momenta_vector<decltype(vd),1>(vd,mv);
+		vd.getProp<0>(p) = 5.0;
 
-	BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
-	BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
+		++it;
+	}
 
-	momenta_grid<decltype(gd),2>(gd,mg);
-	momenta_vector<decltype(vd),2>(vd,mv);
+	vd.map();
 
-	BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
-	BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
+	// Reset the grid
+	interp_test<2,float>(gd,vd,false);
+
+	float mg[2];
+	float mv[2];
 
 	auto & v_cl = create_vcluster();
 
+	interpolate<decltype(vd),decltype(gd),mp4_kernel<float>> inte(vd,gd);
+
 	// We have to do a ghost get before interpolating m2p
 	// Before doing mesh to particle particle must be arranged
 	// into a grid like
@@ -254,8 +314,7 @@ BOOST_AUTO_TEST_CASE( interpolation_full_test_2D )
 	}
 }
 
-
-BOOST_AUTO_TEST_CASE( interpolation_full_test_3D )
+BOOST_AUTO_TEST_CASE( interpolation_full_single_test_3D )
 {
 	Box<3,double> domain({0.0,0.0,0.0},{1.0,1.0,1.0});
 	size_t sz[3] = {64,64,64};
@@ -265,7 +324,6 @@ BOOST_AUTO_TEST_CASE( interpolation_full_test_3D )
 
 	size_t bc_v[3] = {PERIODIC,PERIODIC,PERIODIC};
 
-	{
 	vector_dist<3,double,aggregate<double>> vd(65536,domain,bc_v,gv);
 	grid_dist_id<3,double,aggregate<double>> gd(vd.getDecomposition(),sz,gg);
 
@@ -289,47 +347,51 @@ BOOST_AUTO_TEST_CASE( interpolation_full_test_3D )
 	vd.map();
 
 	// Reset the grid
+	interp_test<3,double>(gd,vd,true);
+}
 
-	auto it2 = gd.getDomainGhostIterator();
+BOOST_AUTO_TEST_CASE( interpolation_full_test_3D )
+{
+	Box<3,double> domain({0.0,0.0,0.0},{1.0,1.0,1.0});
+	size_t sz[3] = {64,64,64};
 
-	while (it2.isNext())
-	{
-		auto key = it2.get();
+	Ghost<3,long int> gg(2);
+	Ghost<3,double> gv(0.01);
 
-		gd.template get<0>(key) = 0.0;
+	size_t bc_v[3] = {PERIODIC,PERIODIC,PERIODIC};
 
-		++it2;
-	}
+	{
+	vector_dist<3,double,aggregate<double>> vd(65536,domain,bc_v,gv);
+	grid_dist_id<3,double,aggregate<double>> gd(vd.getDecomposition(),sz,gg);
 
-	interpolate<decltype(vd),decltype(gd),mp4_kernel<double>> inte(vd,gd);
+	// set one particle on vd
 
-	inte.p2m<0,0>(vd,gd);
+	auto it = vd.getDomainIterator();
 
-	double mg[3];
-	double mv[3];
+	while (it.isNext())
+	{
+		auto p = it.get();
 
-	momenta_grid<decltype(gd),0>(gd,mg);
-	momenta_vector<decltype(vd),0>(vd,mv);
+		vd.getPos(p)[0] = (double)rand()/RAND_MAX;
+		vd.getPos(p)[1] = (double)rand()/RAND_MAX;
+		vd.getPos(p)[2] = (double)rand()/RAND_MAX;
 
-	BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
-	BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
-	BOOST_REQUIRE_CLOSE(mg[2],mv[2],0.001);
+		vd.getProp<0>(p) = 5.0;
 
-	momenta_grid<decltype(gd),1>(gd,mg);
-	momenta_vector<decltype(vd),1>(vd,mv);
+		++it;
+	}
 
-	BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
-	BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
-	BOOST_REQUIRE_CLOSE(mg[2],mv[2],0.001);
+	vd.map();
 
-	momenta_grid<decltype(gd),2>(gd,mg);
-	momenta_vector<decltype(vd),2>(vd,mv);
+	// Reset the grid
 
-	BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
-	BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
-	BOOST_REQUIRE_CLOSE(mg[2],mv[2],0.001);
+	// Reset the grid
+	interp_test<3,double>(gd,vd,false);
 
 	auto & v_cl = create_vcluster();
+	double mg[3];
+	double mv[3];
+	interpolate<decltype(vd),decltype(gd),mp4_kernel<double>> inte(vd,gd);
 
 	// We have to do a ghost get before interpolating m2p
 	// Before doing mesh to particle particle must be arranged