diff --git a/Examples/ChebyshevInterpolationFMM.cpp b/Examples/ChebyshevInterpolationFMM.cpp
index 70aaf4b8279170673998397ae1d5a66dbd4eb1d5..d2b7f06338737f2911d97fa9f1d57188ae917e3a 100755
--- a/Examples/ChebyshevInterpolationFMM.cpp
+++ b/Examples/ChebyshevInterpolationFMM.cpp
@@ -122,6 +122,7 @@ int main(int argc, char* argv[])
 	typedef FOctree<CellClass,ContainerClass,LeafClass>  OctreeClass;
 	//
 	typedef FInterpMatrixKernelR                                                                              MatrixKernelClass;
+  const MatrixKernelClass MatrixKernel;
 	typedef FChebSymKernel<CellClass,ContainerClass,MatrixKernelClass,ORDER>  KernelClass;
 	//
 #ifdef _OPENMP
@@ -161,7 +162,7 @@ int main(int argc, char* argv[])
 
 		time.tic();
 		//
-		KernelClass kernels(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
+		KernelClass kernels(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel);
 		//
 		FmmClass algorithm(&tree, &kernels);
 		//
diff --git a/Src/Kernels/Chebyshev/FAbstractChebKernel.hpp b/Src/Kernels/Chebyshev/FAbstractChebKernel.hpp
index 776cf0dd81d252794be91d9b39d600f0d5fe2a25..7fe03414578b05a12ff10d2f4f98584733165d0c 100755
--- a/Src/Kernels/Chebyshev/FAbstractChebKernel.hpp
+++ b/Src/Kernels/Chebyshev/FAbstractChebKernel.hpp
@@ -49,8 +49,6 @@ protected:
 
   /// Needed for P2M, M2M, L2L and L2P operators
   const FSmartPointer<InterpolatorClass,FSmartPointerMemory> Interpolator;
-  /// Needed for P2P operator
-  /*const*/ FSmartPointer<MatrixKernelClass,FSmartPointerMemory> MatrixKernel;
   /// Height of the entire oct-tree
   const unsigned int TreeHeight;
   /// Corner of oct-tree box
@@ -61,9 +59,6 @@ protected:
   const FReal BoxWidthLeaf;
   /// Extension of the box width ( same for all level! )
   const FReal BoxWidthExtension;
-   
-  /// Parameter to pass to matrix kernel (material specific or anything)
-  const FReal MatParam;
 
   /**
    * Compute center of leaf cell from its tree coordinate.
@@ -86,18 +81,15 @@ public:
   FAbstractChebKernel(const int inTreeHeight,
                       const FReal inBoxWidth,
                       const FPoint& inBoxCenter,
-                      const FReal inBoxWidthExtension = 0.0,
-                      const FReal inMatParam = 0.0)
+                      const FReal inBoxWidthExtension = 0.0)
     : Interpolator(new InterpolatorClass(inTreeHeight,
                                          inBoxWidth,
                                          inBoxWidthExtension)),
-      MatrixKernel(new MatrixKernelClass(inMatParam)),
       TreeHeight(inTreeHeight),
       BoxCorner(inBoxCenter - inBoxWidth / FReal(2.)),
       BoxWidth(inBoxWidth),
       BoxWidthLeaf(BoxWidth / FReal(FMath::pow(2, inTreeHeight - 1))),
-      BoxWidthExtension(inBoxWidthExtension),
-      MatParam(inMatParam)
+      BoxWidthExtension(inBoxWidthExtension)
   {
     /* empty */
   }
@@ -109,10 +101,6 @@ public:
   const InterpolatorClass * getPtrToInterpolator() const
   { return Interpolator.getPtr(); }
 
-   MatrixKernelClass  *  getPtrMatrixKernel()
-  {return   MatrixKernel; }
-
-
 
   virtual void P2M(CellClass* const LeafCell,
 		   const ContainerClass* const SourceParticles) = 0;
diff --git a/Src/Kernels/Chebyshev/FChebKernel.hpp b/Src/Kernels/Chebyshev/FChebKernel.hpp
index 17c7adc8b076abd75333e5bf9d764054c9c10a9f..6dbd318493e5b45cd7528c300719f3023d47a31b 100755
--- a/Src/Kernels/Chebyshev/FChebKernel.hpp
+++ b/Src/Kernels/Chebyshev/FChebKernel.hpp
@@ -51,6 +51,9 @@ class FChebKernel
     typedef FAbstractChebKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
 	AbstractBaseClass;
 
+  /// Needed for P2P and M2L operators
+  const MatrixKernelClass *const MatrixKernel;
+
 	/// Needed for M2L operator
 	FSmartPointer<  M2LHandlerClass,FSmartPointerMemory> M2LHandler;
 
@@ -60,20 +63,21 @@ public:
 	 * precomputed and compressed M2L operators from a binary file (an
 	 * runtime_error is thrown if the required file is not valid).
 	 */
-	FChebKernel(const int inTreeHeight,  const FReal inBoxWidth, const FPoint& inBoxCenter,
-		    const FReal Epsilon)
-        : FAbstractChebKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>(inTreeHeight,
-											   inBoxWidth,
-											   inBoxCenter),
-	  M2LHandler(new M2LHandlerClass(Epsilon))
+	FChebKernel(const int inTreeHeight,  const FReal inBoxWidth, const FPoint& inBoxCenter, const MatrixKernelClass *const inMatrixKernel,
+              const FReal Epsilon)
+    : FAbstractChebKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>(inTreeHeight,
+                                                                                       inBoxWidth,
+                                                                                       inBoxCenter),
+      MatrixKernel(inMatrixKernel),
+      M2LHandler(new M2LHandlerClass(MatrixKernel,Epsilon))
 	{
 		// read precomputed compressed m2l operators from binary file
 		//M2LHandler->ReadFromBinaryFileAndSet();
 		M2LHandler->ComputeAndCompressAndSet();
 	}
 
-	FChebKernel(const int inTreeHeight, const FReal inBoxWidth, const FPoint& inBoxCenter)
-		: 	FChebKernel(inTreeHeight, inBoxWidth,inBoxCenter,FMath::pow(10.0,static_cast<FReal>(-ORDER)))
+	FChebKernel(const int inTreeHeight, const FReal inBoxWidth, const FPoint& inBoxCenter, const MatrixKernelClass *const inMatrixKernel)
+		: 	FChebKernel(inTreeHeight, inBoxWidth,inBoxCenter,inMatrixKernel,FMath::pow(10.0,static_cast<FReal>(-ORDER)))
 	{
 
 	}
@@ -212,14 +216,14 @@ public:
                      ContainerClass* const NeighborSourceParticles[27],
                      const int /* size */)
     {
-        DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2P(TargetParticles,NeighborSourceParticles,AbstractBaseClass::MatrixKernel.getPtr());
+        DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2P(TargetParticles,NeighborSourceParticles,MatrixKernel);
     }
 
 
     void P2PRemote(const FTreeCoordinate& /*inPosition*/,
                    ContainerClass* const FRestrict inTargets, const ContainerClass* const FRestrict /*inSources*/,
                    ContainerClass* const inNeighbors[27], const int /*inSize*/){
-        DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2PRemote(inTargets,inNeighbors,27,AbstractBaseClass::MatrixKernel.getPtr());
+        DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2PRemote(inTargets,inNeighbors,27,MatrixKernel);
     }
 
 };
diff --git a/Src/Kernels/Chebyshev/FChebM2LHandler.hpp b/Src/Kernels/Chebyshev/FChebM2LHandler.hpp
index 782ef5963f2172387ea5c0a3f58779a3172c3e24..ff2fa37f8a5fdca5db0eb74135016abd2c9b166a 100755
--- a/Src/Kernels/Chebyshev/FChebM2LHandler.hpp
+++ b/Src/Kernels/Chebyshev/FChebM2LHandler.hpp
@@ -70,7 +70,7 @@ class FChebM2LHandler : FNoCopyable
 				nnodes = TensorTraits<ORDER>::nnodes,
 				ninteractions = 316}; // 7^3 - 3^3 (max num cells in far-field)
 
-	const MatrixKernelClass MatrixKernel;
+	const MatrixKernelClass *const MatrixKernel;
 
 	FReal *U, *C, *B;
 	const FReal epsilon; //<! accuracy which determines trucation of SVD
@@ -88,8 +88,8 @@ class FChebM2LHandler : FNoCopyable
 
 	
 public:
-	FChebM2LHandler(const FReal _epsilon)
-		: MatrixKernel(), U(nullptr), C(nullptr), B(nullptr), epsilon(_epsilon), rank(0)
+	FChebM2LHandler(const MatrixKernelClass *const inMatrixKernel, const FReal _epsilon)
+		: MatrixKernel(inMatrixKernel), U(nullptr), C(nullptr), B(nullptr), epsilon(_epsilon), rank(0)
 	{}
 
 	~FChebM2LHandler()
@@ -108,7 +108,7 @@ public:
 		FTic time; time.tic();
 		// check if aready set
 		if (U||C||B) throw std::runtime_error("Compressed M2L operator already set");
-		rank = ComputeAndCompress(epsilon, U, C, B);
+		rank = ComputeAndCompress(MatrixKernel, epsilon, U, C, B);
 
     unsigned long sizeM2L = 343*rank*rank*sizeof(FReal);
 
@@ -135,13 +135,13 @@ public:
 	 * @param[out] C matrix of size \f$r\times 316 r\f$ storing \f$[C_1,\dots,C_{316}]\f$
 	 * @param[out] B matrix of size \f$\ell^3\times r\f$
 	 */
-	static unsigned int ComputeAndCompress(const FReal epsilon, FReal* &U, FReal* &C, FReal* &B);
+	static unsigned int ComputeAndCompress(const MatrixKernelClass *const MatrixKernel, const FReal epsilon, FReal* &U, FReal* &C, FReal* &B);
 
 	/**
 	 * Computes, compresses and stores the matrices \f$Y, C_t, B\f$ in a binary
 	 * file
 	 */
-	static void ComputeAndCompressAndStoreInBinaryFile(const FReal epsilon);
+	static void ComputeAndCompressAndStoreInBinaryFile(const MatrixKernelClass *const MatrixKernel, const FReal epsilon);
 
 	/**
 	 * Reads the matrices \f$Y, C_t, B\f$ from the respective binary file
@@ -182,19 +182,19 @@ public:
   {
 		const unsigned int idx
 			= (transfer[0]+3)*7*7 + (transfer[1]+3)*7 + (transfer[2]+3);
-		const FReal scale(MatrixKernel.getScaleFactor(CellWidth));
+		const FReal scale(MatrixKernel->getScaleFactor(CellWidth));
     FBlas::gemva(rank, rank, scale, C + idx*rank*rank, const_cast<FReal*>(Y), X);
   }
   void applyC(const unsigned int idx, FReal CellWidth,
 							const FReal *const Y, FReal *const X) const
   {
-		const FReal scale(MatrixKernel.getScaleFactor(CellWidth));
+		const FReal scale(MatrixKernel->getScaleFactor(CellWidth));
     FBlas::gemva(rank, rank, scale, C + idx*rank*rank, const_cast<FReal*>(Y), X);
   }
   void applyC(FReal CellWidth,
 							const FReal *const Y, FReal *const X) const
   {
-		const FReal scale(MatrixKernel.getScaleFactor(CellWidth));
+		const FReal scale(MatrixKernel->getScaleFactor(CellWidth));
     FBlas::gemva(rank, rank * 343, scale, C, const_cast<FReal*>(Y), X);
   }
 
@@ -229,7 +229,8 @@ public:
 
 template <int ORDER, class MatrixKernelClass>
 unsigned int
-FChebM2LHandler<ORDER, MatrixKernelClass>::ComputeAndCompress(const FReal epsilon,
+FChebM2LHandler<ORDER, MatrixKernelClass>::ComputeAndCompress(const MatrixKernelClass *const MatrixKernel, 
+                                                              const FReal epsilon,
 																															FReal* &U,
 																															FReal* &C,
 																															FReal* &B)
@@ -241,8 +242,6 @@ FChebM2LHandler<ORDER, MatrixKernelClass>::ComputeAndCompress(const FReal epsilo
 	FPoint X[nnodes], Y[nnodes];
 	// set roots of target cell (X)
 	FChebTensor<order>::setRoots(FPoint(0.,0.,0.), FReal(2.), X);
-	// init matrix kernel
-	const MatrixKernelClass MatrixKernel;
 
 	// allocate memory and compute 316 m2l operators
 	FReal *_U, *_C, *_B;
@@ -260,7 +259,7 @@ FChebM2LHandler<ORDER, MatrixKernelClass>::ComputeAndCompress(const FReal epsilo
 					for (unsigned int n=0; n<nnodes; ++n)
 						for (unsigned int m=0; m<nnodes; ++m)
 							_C[counter*nnodes*nnodes + n*nnodes + m]
-								= MatrixKernel.evaluate(X[m], Y[n]);
+								= MatrixKernel->evaluate(X[m], Y[n]);
 					// increment interaction counter
 					counter++;
 				}
@@ -330,14 +329,14 @@ FChebM2LHandler<ORDER, MatrixKernelClass>::ComputeAndCompress(const FReal epsilo
 
 template <int ORDER, class MatrixKernelClass>
 void
-FChebM2LHandler<ORDER, MatrixKernelClass>::ComputeAndCompressAndStoreInBinaryFile(const FReal epsilon)
+FChebM2LHandler<ORDER, MatrixKernelClass>::ComputeAndCompressAndStoreInBinaryFile(const MatrixKernelClass *const MatrixKernel, const FReal epsilon)
 {
 	// measure time
 	FTic time; time.tic();
 	// start computing process
 	FReal *U, *C, *B;
 	U = C = B = nullptr;
-	const unsigned int rank = ComputeAndCompress(epsilon, U, C, B);
+	const unsigned int rank = ComputeAndCompress(MatrixKernel, epsilon, U, C, B);
 	// store into binary file
 	const std::string filename(getFileName(epsilon));
 	std::ofstream stream(filename.c_str(),
diff --git a/Src/Kernels/Chebyshev/FChebSymKernel.hpp b/Src/Kernels/Chebyshev/FChebSymKernel.hpp
index 93ed3e065006787de4e9bd638f1aafa257f9b15c..df3394853f2ce41bbf2a52ef17da9051771cfe73 100755
--- a/Src/Kernels/Chebyshev/FChebSymKernel.hpp
+++ b/Src/Kernels/Chebyshev/FChebSymKernel.hpp
@@ -56,6 +56,9 @@ class FChebSymKernel
 	typedef SymmetryHandler<ORDER, MatrixKernelClass::Type> SymmetryHandlerClass;
   enum {nnodes = AbstractBaseClass::nnodes};
 
+  /// Needed for P2P and M2L operators
+  const MatrixKernelClass *const MatrixKernel;
+
 	/// Needed for handling all symmetries
 	const FSmartPointer<SymmetryHandlerClass,FSmartPointerMemory> SymHandler;
 
@@ -105,11 +108,13 @@ public:
 	 * runtime_error is thrown if the required file is not valid).
 	 */
 	FChebSymKernel(const int inTreeHeight,
-		       const FReal inBoxWidth,
-		       const FPoint& inBoxCenter,
-		       const FReal Epsilon)
+                 const FReal inBoxWidth,
+                 const FPoint& inBoxCenter,
+                 const MatrixKernelClass *const inMatrixKernel,
+                 const FReal Epsilon)
 	  : AbstractBaseClass(inTreeHeight, inBoxWidth, inBoxCenter),
-			SymHandler(new SymmetryHandlerClass(AbstractBaseClass::MatrixKernel.getPtr(), Epsilon, inBoxWidth, inTreeHeight)),
+      MatrixKernel(inMatrixKernel),
+			SymHandler(new SymmetryHandlerClass(MatrixKernel, Epsilon, inBoxWidth, inTreeHeight)),
 			Loc(nullptr), Mul(nullptr), countExp(nullptr)
 	{
 		this->allocateMemoryForPermutedExpansions();
@@ -126,8 +131,9 @@ public:
 	 * runtime_error is thrown if the required file is not valid).
 	 */
 	FChebSymKernel(const int inTreeHeight,
-		       const FReal inBoxWidth,
-		       const FPoint& inBoxCenter) :FChebSymKernel(inTreeHeight,inBoxWidth, inBoxCenter,FMath::pow(10.0,static_cast<FReal>(-ORDER)))
+                 const FReal inBoxWidth,
+                 const FPoint& inBoxCenter,
+                 const MatrixKernelClass *const inMatrixKernel) :FChebSymKernel(inTreeHeight,inBoxWidth, inBoxCenter,inMatrixKernel,FMath::pow(10.0,static_cast<FReal>(-ORDER)))
 		       {}
 
 	
@@ -135,6 +141,7 @@ public:
 	/** Copy constructor */
 	FChebSymKernel(const FChebSymKernel& other)
 		: AbstractBaseClass(other),
+      MatrixKernel(other.MatrixKernel),
 			SymHandler(other.SymHandler),
 			Loc(nullptr), Mul(nullptr), countExp(nullptr)
 	{
@@ -264,7 +271,7 @@ public:
 
             // multiply (mat-mat-mul)
             FReal Compressed [nnodes * 24];
-            const FReal scale = AbstractBaseClass::MatrixKernel->getScaleFactor(AbstractBaseClass::BoxWidth, TreeLevel);
+            const FReal scale = MatrixKernel->getScaleFactor(AbstractBaseClass::BoxWidth, TreeLevel);
             for (unsigned int pidx=0; pidx<343; ++pidx) {
                 const unsigned int count = countExp[pidx];
                 if (count) {
@@ -462,14 +469,14 @@ public:
                      ContainerClass* const NeighborSourceParticles[27],
                      const int /* size */)
     {
-        DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2P(TargetParticles,NeighborSourceParticles,AbstractBaseClass::MatrixKernel.getPtr());
+        DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2P(TargetParticles,NeighborSourceParticles,MatrixKernel);
     }
 
 
     void P2PRemote(const FTreeCoordinate& /*inPosition*/,
                    ContainerClass* const FRestrict inTargets, const ContainerClass* const FRestrict /*inSources*/,
                    ContainerClass* const inNeighbors[27], const int /*inSize*/){
-        DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2PRemote(inTargets,inNeighbors,27,AbstractBaseClass::MatrixKernel.getPtr());
+        DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2PRemote(inTargets,inNeighbors,27,MatrixKernel);
     }
 
 };
diff --git a/Src/Kernels/Chebyshev/FChebSymM2LHandler.hpp b/Src/Kernels/Chebyshev/FChebSymM2LHandler.hpp
index b24225822225362fd34939b9ff20dfc4e4266333..28e6e324c7b94a78212a29b4309d5e1c164b8326 100755
--- a/Src/Kernels/Chebyshev/FChebSymM2LHandler.hpp
+++ b/Src/Kernels/Chebyshev/FChebSymM2LHandler.hpp
@@ -303,7 +303,7 @@ static void precompute(const MatrixKernelClass *const MatrixKernel, const FReal
 				FChebTensor<ORDER>::setRootOfWeights(weights);
 
 				// now the entry-computer is responsible for weighting the matrix entries
-				EntryComputer<MatrixKernelClass> Computer(nnodes, X, nnodes, Y, weights);
+				EntryComputer<MatrixKernelClass> Computer(MatrixKernel, nnodes, X, nnodes, Y, weights);
 
 				// start timer
 				time.tic();
diff --git a/Src/Kernels/Chebyshev/FChebTensorialKernel.hpp b/Src/Kernels/Chebyshev/FChebTensorialKernel.hpp
index 4966e08429b31f83ddee3ae2180cf5e4b3d85388..dcf24d479badfe3dbe434b997d8710649cfa1e66 100755
--- a/Src/Kernels/Chebyshev/FChebTensorialKernel.hpp
+++ b/Src/Kernels/Chebyshev/FChebTensorialKernel.hpp
@@ -59,6 +59,9 @@ protected://PB: for OptiDis
     typedef FAbstractChebKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
 	AbstractBaseClass;
 
+  /// Needed for P2P and M2L operators
+  const MatrixKernelClass *const MatrixKernel;
+
 	/// Needed for M2L operator
 	FSmartPointer<  M2LHandlerClass,FSmartPointerMemory> M2LHandler;
 
@@ -71,11 +74,12 @@ public:
 	FChebTensorialKernel(const int inTreeHeight,
                        const FReal inBoxWidth,
                        const FPoint& inBoxCenter,
+                       const MatrixKernelClass *const inMatrixKernel,
                        const FReal inBoxWidthExtension, 
-                       const FReal Epsilon,
-                       const double inMatParam = 0.0)
-    : FAbstractChebKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>(inTreeHeight,inBoxWidth,inBoxCenter,inBoxWidthExtension,inMatParam),
-	  M2LHandler(new M2LHandlerClass(AbstractBaseClass::MatrixKernel.getPtr(),
+                       const FReal Epsilon)
+    : FAbstractChebKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>(inTreeHeight,inBoxWidth,inBoxCenter,inBoxWidthExtension),
+      MatrixKernel(inMatrixKernel),
+	  M2LHandler(new M2LHandlerClass(MatrixKernel,
                                    inTreeHeight,
                                    inBoxWidth,
                                    inBoxWidthExtension,
@@ -143,7 +147,7 @@ public:
     // scale factor for homogeneous case
     const FReal CellWidth(AbstractBaseClass::BoxWidth / FReal(FMath::pow(2, TreeLevel)));
     const FReal ExtendedCellWidth(CellWidth + AbstractBaseClass::BoxWidthExtension);
-    const FReal scale(AbstractBaseClass::MatrixKernel.getPtr()->getScaleFactor(ExtendedCellWidth));
+    const FReal scale(MatrixKernel->getScaleFactor(ExtendedCellWidth));
 
     for(int idxV = 0 ; idxV < NVALS ; ++idxV){
       for (int idxLhs=0; idxLhs < nLhs; ++idxLhs){
@@ -158,8 +162,7 @@ public:
           const int idxMul = idxV*nRhs + idxRhs;
 
           // get index in matrix kernel
-          const unsigned int d 
-            = AbstractBaseClass::MatrixKernel.getPtr()->getPosition(idxLhs);
+          const unsigned int d = MatrixKernel->getPosition(idxLhs);
 
           for (int idx=0; idx<343; ++idx){
             if (SourceCells[idx]){
@@ -237,14 +240,14 @@ public:
                      ContainerClass* const NeighborSourceParticles[27],
                      const int /* size */)
     {
-        DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2P(TargetParticles,NeighborSourceParticles,AbstractBaseClass::MatrixKernel.getPtr());
+        DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2P(TargetParticles,NeighborSourceParticles,MatrixKernel);
     }
 
 
     void P2PRemote(const FTreeCoordinate& /*inPosition*/,
                    ContainerClass* const FRestrict inTargets, const ContainerClass* const FRestrict /*inSources*/,
                    ContainerClass* const inNeighbors[27], const int /*inSize*/){
-        DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2PRemote(inTargets,inNeighbors,27,AbstractBaseClass::MatrixKernel.getPtr());
+        DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2PRemote(inTargets,inNeighbors,27,MatrixKernel);
     }
 
 };
diff --git a/Src/Kernels/Chebyshev/FChebTensorialM2LHandler.hpp b/Src/Kernels/Chebyshev/FChebTensorialM2LHandler.hpp
index 40614bfd368fa07e91a2008da2192c0309213ccc..d65a7c594fe29426bc5567266a52983ee209a8da 100755
--- a/Src/Kernels/Chebyshev/FChebTensorialM2LHandler.hpp
+++ b/Src/Kernels/Chebyshev/FChebTensorialM2LHandler.hpp
@@ -82,8 +82,6 @@ class FChebTensorialM2LHandler<ORDER,MatrixKernelClass,HOMOGENEOUS> : FNoCopyabl
 				ninteractions = 316,// 7^3 - 3^3 (max num cells in far-field)
         ncmp = MatrixKernelClass::NCMP};
 
-//	const MatrixKernelClass MatrixKernel;
-
 //	FReal *U, *C, *B;
 	FReal *U, *B;
 	FReal** C;
diff --git a/Src/Kernels/Interpolation/FInterpMatrixKernel.hpp b/Src/Kernels/Interpolation/FInterpMatrixKernel.hpp
index 16a15dee82899ec836a6345a21cf16e8f60d0ad1..52a19d7e4f85eb6096015821960559db680b72ef 100755
--- a/Src/Kernels/Interpolation/FInterpMatrixKernel.hpp
+++ b/Src/Kernels/Interpolation/FInterpMatrixKernel.hpp
@@ -80,7 +80,11 @@ struct FInterpMatrixKernelR : FInterpAbstractMatrixKernel
   static const unsigned int NRHS = 1; //< dim of mult exp
   static const unsigned int NLHS = 1; //< dim of loc exp
 
-  FInterpMatrixKernelR(const FReal = 0.0, const unsigned int = 0) {}
+  FInterpMatrixKernelR() {}
+
+  // copy ctor
+  FInterpMatrixKernelR(const FInterpMatrixKernelR& other) 
+  {}
 
   // returns position in reduced storage
   int getPosition(const unsigned int) const
@@ -142,9 +146,14 @@ struct FInterpMatrixKernelRH :FInterpMatrixKernelR{
 	  static const unsigned int NLHS = 1; //< dim of loc exp
 	  FReal LX,LY,LZ ;
 
-	FInterpMatrixKernelRH(const FReal = 0.0, const unsigned int = 0) : LX(1.0),LY(1.0),LZ(1.0)
+	FInterpMatrixKernelRH() : LX(1.0),LY(1.0),LZ(1.0)
 	{	 }
 
+  // copy ctor
+  FInterpMatrixKernelRH(const FInterpMatrixKernelRH& other) 
+  : LX(other.LX), LY(other.LY), LZ(other.LZ) 
+  {}
+
   // evaluate interaction
 	  FReal evaluate(const FPoint& x, const FPoint& y) const
 	  {
@@ -207,7 +216,11 @@ struct FInterpMatrixKernelRR : FInterpAbstractMatrixKernel
   static const unsigned int NRHS = 1; //< dim of mult exp
   static const unsigned int NLHS = 1; //< dim of loc exp
 
-  FInterpMatrixKernelRR(const FReal = 0.0, const unsigned int = 0) {}
+  FInterpMatrixKernelRR() {}
+
+  // copy ctor
+  FInterpMatrixKernelRR(const FInterpMatrixKernelRR& other) 
+  {}
 
   // returns position in reduced storage
   int getPosition(const unsigned int) const
@@ -272,7 +285,11 @@ struct FInterpMatrixKernelLJ : FInterpAbstractMatrixKernel
   static const unsigned int NRHS = 1; //< dim of mult exp
   static const unsigned int NLHS = 1; //< dim of loc exp
 
-  FInterpMatrixKernelLJ(const FReal = 0.0, const unsigned int = 0) {}
+  FInterpMatrixKernelLJ() {}
+
+  // copy ctor
+  FInterpMatrixKernelLJ(const FInterpMatrixKernelLJ& other) 
+  {}
 
   // returns position in reduced storage
   int getPosition(const unsigned int) const
@@ -388,6 +405,11 @@ struct FInterpMatrixKernel_R_IJ : FInterpAbstractMatrixKernel
     : _i(indexTab[d]), _j(indexTab[d+NCMP]), _CoreWidth2(CoreWidth*CoreWidth)
   {}
 
+  // copy ctor
+  FInterpMatrixKernel_R_IJ(const FInterpMatrixKernel_R_IJ& other) 
+  : _i(other._i), _j(other._j), _CoreWidth2(other._CoreWidth2) 
+  {}
+
   // returns position in reduced storage from position in full 3x3 matrix
  unsigned  int getPosition(const unsigned int n) const
   {return applyTab[n];} 
@@ -487,6 +509,11 @@ struct FInterpMatrixKernel_R_IJK : FInterpAbstractMatrixKernel
   : _i(indexTab[d]), _j(indexTab[d+NCMP]), _k(indexTab[d+2*NCMP]), _CoreWidth2(CoreWidth*CoreWidth)
   {}
 
+  // copy ctor
+  FInterpMatrixKernel_R_IJK(const FInterpMatrixKernel_R_IJK& other) 
+  : _i(other._i), _j(other._j), _k(other._k), _CoreWidth2(other._CoreWidth2) 
+  {}
+
   // returns position in reduced storage from position in full 3x3x3 matrix
   unsigned int getPosition(const unsigned int n) const
   {return applyTab[n];} 
@@ -598,10 +625,10 @@ struct FInterpMatrixKernel_R_IJK : FInterpAbstractMatrixKernel
   number of rows and cols and on the coordinates x and y and the type of the
   generating matrix-kernel function.
 */
-template <typename MatrixKernel>
+template <typename MatrixKernelClass>
 class EntryComputer
 {
-  const MatrixKernel Kernel;
+  const MatrixKernelClass *const MatrixKernel;
 
   const unsigned int nx, ny;
   const FPoint *const px, *const py;
@@ -609,12 +636,11 @@ class EntryComputer
   const FReal *const weights;
 
 public:
-  explicit EntryComputer(const unsigned int _nx, const FPoint *const _px,
+  explicit EntryComputer(const MatrixKernelClass *const inMatrixKernel,
+                         const unsigned int _nx, const FPoint *const _px,
                          const unsigned int _ny, const FPoint *const _py,
-                         const FReal *const _weights = NULL,
-                         const unsigned int idxK = 0,
-                         const FReal matparam = 0.0)
-    : Kernel(matparam,idxK),	nx(_nx), ny(_ny), px(_px), py(_py), weights(_weights) {}
+                         const FReal *const _weights = NULL)
+    : MatrixKernel(inMatrixKernel),	nx(_nx), ny(_ny), px(_px), py(_py), weights(_weights) {}
 
   void operator()(const unsigned int xbeg, const unsigned int xend,
                   const unsigned int ybeg, const unsigned int yend,
@@ -624,11 +650,11 @@ public:
     if (weights) {
       for (unsigned int j=ybeg; j<yend; ++j)
         for (unsigned int i=xbeg; i<xend; ++i)
-          data[idx++] = weights[i] * weights[j] * Kernel.evaluate(px[i], py[j]);
+          data[idx++] = weights[i] * weights[j] * MatrixKernel->evaluate(px[i], py[j]);
     } else {
       for (unsigned int j=ybeg; j<yend; ++j)
         for (unsigned int i=xbeg; i<xend; ++i)
-          data[idx++] = Kernel.evaluate(px[i], py[j]);
+          data[idx++] = MatrixKernel->evaluate(px[i], py[j]);
     }
 
     /*
diff --git a/Src/Kernels/Uniform/FAbstractUnifKernel.hpp b/Src/Kernels/Uniform/FAbstractUnifKernel.hpp
index 90c847754a9d93b3e65cc43d8c636151e9ef7c63..d9529a326b9a3d6cbbfc4555378e764e4a302a9f 100644
--- a/Src/Kernels/Uniform/FAbstractUnifKernel.hpp
+++ b/Src/Kernels/Uniform/FAbstractUnifKernel.hpp
@@ -49,8 +49,6 @@ protected:
 
   /// Needed for P2M, M2M, L2L and L2P operators
   const FSmartPointer<InterpolatorClass,FSmartPointerMemory> Interpolator;
-  /// Needed for P2P operator
-  const FSmartPointer<MatrixKernelClass,FSmartPointerMemory> MatrixKernel;
   /// Height of the entire oct-tree
   const unsigned int TreeHeight;
   /// Corner of oct-tree box
@@ -62,9 +60,6 @@ protected:
   /// Extension of the box width ( same for all level! )
   const FReal BoxWidthExtension;
 
-  /// Parameter to pass to matrix kernel (material specific or anything)
-  const double MatParam;
-
   /**
    * Compute center of leaf cell from its tree coordinate.
    * @param[in] Coordinate tree coordinate
@@ -86,18 +81,15 @@ public:
   FAbstractUnifKernel(const int inTreeHeight,
                       const FReal inBoxWidth,
                       const FPoint& inBoxCenter,
-                      const FReal inBoxWidthExtension = 0.0,
-                      const double inMatParam = 0.0)
+                      const FReal inBoxWidthExtension = 0.0)
     : Interpolator(new InterpolatorClass(inTreeHeight,
                                          inBoxWidth,
                                          inBoxWidthExtension)),
-      MatrixKernel(new MatrixKernelClass(inMatParam)),
       TreeHeight(inTreeHeight),
       BoxCorner(inBoxCenter - inBoxWidth / FReal(2.)),
       BoxWidth(inBoxWidth),
       BoxWidthLeaf(BoxWidth / FReal(FMath::pow(2, inTreeHeight - 1))),
-      BoxWidthExtension(inBoxWidthExtension),
-      MatParam(inMatParam)
+      BoxWidthExtension(inBoxWidthExtension)
   {
     /* empty */
   }
diff --git a/Src/Kernels/Uniform/FUnifKernel.hpp b/Src/Kernels/Uniform/FUnifKernel.hpp
index 1fcfad13265cd1f217dedbe6164c57e9751ea40a..6c2324af486c05d4a6c4aad819b2b561d7302a74 100644
--- a/Src/Kernels/Uniform/FUnifKernel.hpp
+++ b/Src/Kernels/Uniform/FUnifKernel.hpp
@@ -51,6 +51,9 @@ class FUnifKernel
   typedef FAbstractUnifKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
   AbstractBaseClass;
 
+  /// Needed for P2P and M2L operators
+  const MatrixKernelClass *const MatrixKernel;
+
   /// Needed for M2L operator
   const M2LHandlerClass M2LHandler;
 
@@ -63,9 +66,11 @@ public:
    */
   FUnifKernel(const int inTreeHeight,
               const FReal inBoxWidth,
-              const FPoint& inBoxCenter)
+              const FPoint& inBoxCenter,
+              const MatrixKernelClass *const inMatrixKernel)
     : FAbstractUnifKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>(inTreeHeight,inBoxWidth,inBoxCenter),
-      M2LHandler(AbstractBaseClass::MatrixKernel.getPtr(),
+      MatrixKernel(inMatrixKernel),
+      M2LHandler(MatrixKernel,
                  inTreeHeight,
                  inBoxWidth) 
   { }
@@ -131,7 +136,7 @@ public:
            const int TreeLevel)
   {
     const FReal CellWidth(AbstractBaseClass::BoxWidth / FReal(FMath::pow(2, TreeLevel)));
-    const FReal scale(AbstractBaseClass::MatrixKernel.getPtr()->getScaleFactor(CellWidth));
+    const FReal scale(MatrixKernel->getScaleFactor(CellWidth));
 
     for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
       FComplexe *const TransformedLocalExpansion = TargetCell->getTransformedLocal(idxRhs);
@@ -208,7 +213,7 @@ public:
            ContainerClass* const NeighborSourceParticles[27],
            const int /* size */)
   {
-    DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2P(TargetParticles,NeighborSourceParticles,AbstractBaseClass::MatrixKernel.getPtr());
+    DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2P(TargetParticles,NeighborSourceParticles,MatrixKernel);
   }
 
 
@@ -216,7 +221,7 @@ public:
                  ContainerClass* const FRestrict inTargets, const ContainerClass* const FRestrict /*inSources*/,
                  ContainerClass* const inNeighbors[27], const int /*inSize*/)
   {
-    DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2PRemote(inTargets,inNeighbors,27,AbstractBaseClass::MatrixKernel.getPtr());
+    DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2PRemote(inTargets,inNeighbors,27,MatrixKernel);
   }
 
 };
diff --git a/Src/Kernels/Uniform/FUnifTensorialKernel.hpp b/Src/Kernels/Uniform/FUnifTensorialKernel.hpp
index 52edf1825e4930a2e3d67cb46e43c35d232855c9..05a5c289cf9f35f0fbba2b7421f47020a3165287 100644
--- a/Src/Kernels/Uniform/FUnifTensorialKernel.hpp
+++ b/Src/Kernels/Uniform/FUnifTensorialKernel.hpp
@@ -75,6 +75,9 @@ protected://PB: for OptiDis
   typedef FAbstractUnifKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
   AbstractBaseClass;
 
+  /// Needed for P2P and M2L operators
+  const MatrixKernelClass *const MatrixKernel;
+
   /// Needed for M2L operator
   const M2LHandlerClass M2LHandler;
 
@@ -87,10 +90,11 @@ public:
   FUnifTensorialKernel(const int inTreeHeight,
                        const FReal inBoxWidth,
                        const FPoint& inBoxCenter,
-                       const FReal inBoxWidthExtension, 
-                       const double inMatParam = 0.0)
-    : FAbstractUnifKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>(inTreeHeight,inBoxWidth,inBoxCenter,inBoxWidthExtension,inMatParam),
-      M2LHandler(AbstractBaseClass::MatrixKernel.getPtr(),
+                       const MatrixKernelClass *const inMatrixKernel,
+                       const FReal inBoxWidthExtension)
+    : FAbstractUnifKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>(inTreeHeight,inBoxWidth,inBoxCenter,inBoxWidthExtension),
+      MatrixKernel(inMatrixKernel),
+      M2LHandler(MatrixKernel,
                  inTreeHeight,
                  inBoxWidth,
                  inBoxWidthExtension) 
@@ -157,7 +161,7 @@ public:
   {
     const FReal CellWidth(AbstractBaseClass::BoxWidth / FReal(FMath::pow(2, TreeLevel)));
     const FReal ExtendedCellWidth(CellWidth + AbstractBaseClass::BoxWidthExtension);
-    const FReal scale(AbstractBaseClass::MatrixKernel.getPtr()->getScaleFactor(ExtendedCellWidth));
+    const FReal scale(MatrixKernel->getScaleFactor(ExtendedCellWidth));
 
     for(int idxV = 0 ; idxV < NVALS ; ++idxV){
       for (int idxLhs=0; idxLhs < nLhs; ++idxLhs){
@@ -175,8 +179,7 @@ public:
           const int idxMul = idxV*nRhs + idxRhs;
 
           // get index in matrix kernel
-          const unsigned int d 
-            = AbstractBaseClass::MatrixKernel.getPtr()->getPosition(idxLhs);
+          const unsigned int d = MatrixKernel->getPosition(idxLhs);
 
           for (int idx=0; idx<343; ++idx){
             if (SourceCells[idx]){
@@ -248,14 +251,14 @@ public:
            ContainerClass* const NeighborSourceParticles[27],
            const int /* size */)
   {
-    DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2P(TargetParticles,NeighborSourceParticles,AbstractBaseClass::MatrixKernel.getPtr());
+    DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2P(TargetParticles,NeighborSourceParticles,MatrixKernel);
   }
 
 
   void P2PRemote(const FTreeCoordinate& /*inPosition*/,
                  ContainerClass* const FRestrict inTargets, const ContainerClass* const FRestrict /*inSources*/,
                  ContainerClass* const inNeighbors[27], const int /*inSize*/){
-    DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2PRemote(inTargets,inNeighbors,27,AbstractBaseClass::MatrixKernel.getPtr());
+    DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2PRemote(inTargets,inNeighbors,27,MatrixKernel);
   }
 
 };
diff --git a/Tests/Kernels/testChebTensorialAlgorithm.cpp b/Tests/Kernels/testChebTensorialAlgorithm.cpp
index 0c928dd3d5d630d0c727198721854faf7c95fbc3..01f762c3c0c14229b1ba518ecff2e8a88a79b15c 100644
--- a/Tests/Kernels/testChebTensorialAlgorithm.cpp
+++ b/Tests/Kernels/testChebTensorialAlgorithm.cpp
@@ -78,7 +78,7 @@ int main(int argc, char* argv[])
   const unsigned int NLHS = MatrixKernelClass::NLHS;
 
   const double CoreWidth = 0.1;
-  const MatrixKernelClass DirectMatrixKernel(CoreWidth);
+  const MatrixKernelClass MatrixKernel(CoreWidth);
   std::cout<< "Interaction kernel: ";
   if(MK_ID == ONE_OVER_R) std::cout<< "1/R";
   else if(MK_ID == R_IJ)
@@ -143,16 +143,16 @@ int main(int argc, char* argv[])
                                      particles[idxOther].position.getZ(), particles[idxOther].physicalValue,
                                      particles[idxOther].forces[0], particles[idxOther].forces[1],
                                      particles[idxOther].forces[2], particles[idxOther].potential,
-                                     &DirectMatrixKernel);
-          else if(MK_ID == ONE_OVER_R)
-            FP2P::MutualParticles(particles[idxTarget].position.getX(), particles[idxTarget].position.getY(),
-                                  particles[idxTarget].position.getZ(), particles[idxTarget].physicalValue[0],
-                                  particles[idxTarget].forces[0], particles[idxTarget].forces[1],
-                                  particles[idxTarget].forces[2], particles[idxTarget].potential,
-                                  particles[idxOther].position.getX(), particles[idxOther].position.getY(),
-                                  particles[idxOther].position.getZ(), particles[idxOther].physicalValue[0],
-                                  particles[idxOther].forces[0], particles[idxOther].forces[1],
-                                  particles[idxOther].forces[2], particles[idxOther].potential);
+                                     &MatrixKernel);
+//          else if(MK_ID == ONE_OVER_R)
+//            FP2P::MutualParticles(particles[idxTarget].position.getX(), particles[idxTarget].position.getY(),
+//                                  particles[idxTarget].position.getZ(), particles[idxTarget].physicalValue[0],
+//                                  particles[idxTarget].forces[0], particles[idxTarget].forces[1],
+//                                  particles[idxTarget].forces[2], particles[idxTarget].potential,
+//                                  particles[idxOther].position.getX(), particles[idxOther].position.getY(),
+//                                  particles[idxOther].position.getZ(), particles[idxOther].physicalValue[0],
+//                                  particles[idxOther].forces[0], particles[idxOther].forces[1],
+//                                  particles[idxOther].forces[2], particles[idxOther].potential,&MatrixKernel);
           else if(MK_ID == R_IJK)
             FP2P::MutualParticlesRIJK(particles[idxTarget].position.getX(), particles[idxTarget].position.getY(),
                                       particles[idxTarget].position.getZ(), particles[idxTarget].physicalValue,
@@ -162,7 +162,7 @@ int main(int argc, char* argv[])
                                       particles[idxOther].position.getZ(), particles[idxOther].physicalValue,
                                       particles[idxOther].forces[0], particles[idxOther].forces[1],
                                       particles[idxOther].forces[2], particles[idxOther].potential,
-                                      &DirectMatrixKernel);
+                                      &MatrixKernel);
           else 
             std::runtime_error("Provide a valid matrix kernel!");
         }
@@ -240,7 +240,7 @@ int main(int argc, char* argv[])
     { // -----------------------------------------------------
       std::cout << "\nChebyshev FMM (ORDER="<< ORDER << ") ... " << std::endl;
       time.tic();
-      KernelClass kernels(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox(), BoxWidthExtension, epsilon, CoreWidth);
+      KernelClass kernels(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox(), &MatrixKernel, BoxWidthExtension, epsilon);
       FmmClass algorithm(&tree, &kernels);
       algorithm.execute();
       time.tac();
diff --git a/Tests/Kernels/testUnifAlgorithm.cpp b/Tests/Kernels/testUnifAlgorithm.cpp
index 599e0ef7665b9165a585d145effb16bf10c14fc1..6dd142667238ce489c4e5ff1b75acad5cba78b5e 100644
--- a/Tests/Kernels/testUnifAlgorithm.cpp
+++ b/Tests/Kernels/testUnifAlgorithm.cpp
@@ -166,7 +166,7 @@ int main(int argc, char* argv[])
     { // -----------------------------------------------------
       std::cout << "\nLagrange/Uniform grid FMM (ORDER="<< ORDER << ") ... " << std::endl;
       time.tic();
-      KernelClass kernels(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
+      KernelClass kernels(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel);
       FmmClass algorithm(&tree, &kernels);
       algorithm.execute();
       time.tac();
diff --git a/Tests/Kernels/testUnifAlgorithmProc.cpp b/Tests/Kernels/testUnifAlgorithmProc.cpp
index aa4451ad18b726e2eeb29fa86d3e41f4de7e19eb..7525f934c39c5b8d10404ad100d8eb9f9366976f 100644
--- a/Tests/Kernels/testUnifAlgorithmProc.cpp
+++ b/Tests/Kernels/testUnifAlgorithmProc.cpp
@@ -87,9 +87,12 @@ int main(int argc, char* argv[])
 #endif
     
   std::cout << "Opening : " <<filename << "\n" << std::endl;
-    // init timer
-    FTic time;
+  // init timer
+  FTic time;
   
+  // Create Matrix Kernel
+  const MatrixKernelClass MatrixKernel;
+
   // init particles position and physical value
   struct TestParticle{
     FPoint position;
@@ -140,7 +143,7 @@ int main(int argc, char* argv[])
   { // -----------------------------------------------------
     std::cout << "\nChebyshev FMM (ORDER="<< ORDER << ",EPS="<< epsilon <<") ... " << std::endl;
     time.tic();
-    KernelClass kernels(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
+    KernelClass kernels(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel);
     FmmClass algorithm(app.global(),&tree, &kernels);
     algorithm.execute();
     time.tac();
diff --git a/Tests/Kernels/testUnifTensorialAlgorithm.cpp b/Tests/Kernels/testUnifTensorialAlgorithm.cpp
index 4ee40d376a3cdbc74db2c6e5cc126852e74f5474..6bb1ce672887cb90d4c9767dc4c0fcbb2854634a 100644
--- a/Tests/Kernels/testUnifTensorialAlgorithm.cpp
+++ b/Tests/Kernels/testUnifTensorialAlgorithm.cpp
@@ -78,8 +78,8 @@ int main(int argc, char* argv[])
   const unsigned int NRHS = MatrixKernelClass::NRHS;
   const unsigned int NLHS = MatrixKernelClass::NLHS;
 
-  const double CoreWidth = 0.1;
-  const MatrixKernelClass DirectMatrixKernel(CoreWidth);
+  const FReal CoreWidth = 0.1;
+  const MatrixKernelClass MatrixKernel(CoreWidth);
   std::cout<< "Interaction kernel: ";
   if(MK_ID == ONE_OVER_R) std::cout<< "1/R";
   else if(MK_ID == R_IJ)
@@ -143,7 +143,7 @@ int main(int argc, char* argv[])
                                      particles[idxOther].position.getZ(), particles[idxOther].physicalValue,
                                      particles[idxOther].forces[0], particles[idxOther].forces[1],
                                      particles[idxOther].forces[2], particles[idxOther].potential,
-                                     &DirectMatrixKernel);
+                                     &MatrixKernel);
           else if(MK_ID == R_IJK)
             FP2P::MutualParticlesRIJK(particles[idxTarget].position.getX(), particles[idxTarget].position.getY(),
                                       particles[idxTarget].position.getZ(), particles[idxTarget].physicalValue,
@@ -153,7 +153,7 @@ int main(int argc, char* argv[])
                                       particles[idxOther].position.getZ(), particles[idxOther].physicalValue,
                                       particles[idxOther].forces[0], particles[idxOther].forces[1],
                                       particles[idxOther].forces[2], particles[idxOther].potential,
-                                      &DirectMatrixKernel);
+                                      &MatrixKernel);
           else 
             std::runtime_error("Provide a valid matrix kernel!");
         }
@@ -234,7 +234,7 @@ int main(int argc, char* argv[])
     { // -----------------------------------------------------
       std::cout << "\nLagrange/Uniform grid FMM (ORDER="<< ORDER << ") ... " << std::endl;
       time.tic();
-      KernelClass kernels(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox(),BoxWidthExtension,CoreWidth);
+      KernelClass kernels(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel,BoxWidthExtension);
       FmmClass algorithm(&tree, &kernels);
       algorithm.execute();
       time.tac();
diff --git a/Tests/noDist/testChebAlgorithm.cpp b/Tests/noDist/testChebAlgorithm.cpp
index 3bc15576c52fb16b3678eafce430eb81e746c9ae..1576d148cfb9335599f19ac4671a32ed7d38d138 100755
--- a/Tests/noDist/testChebAlgorithm.cpp
+++ b/Tests/noDist/testChebAlgorithm.cpp
@@ -159,7 +159,7 @@ int main(int argc, char* argv[])
     { // -----------------------------------------------------
       std::cout << "\nChebyshev FMM (ORDER="<< ORDER << ",EPS="<< epsilon <<") ... " << std::endl;
       time.tic();
-      KernelClass kernels(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
+      KernelClass kernels(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel);
       FmmClass algorithm(&tree, &kernels);
       algorithm.execute();
       time.tac();
diff --git a/Tests/noDist/testChebAlgorithmProc.cpp b/Tests/noDist/testChebAlgorithmProc.cpp
index 979d66c97335673703415f87fc64d038b0fd2360..bdb6c4f18aa3286a567decca25ca209f4b38918f 100644
--- a/Tests/noDist/testChebAlgorithmProc.cpp
+++ b/Tests/noDist/testChebAlgorithmProc.cpp
@@ -87,9 +87,12 @@ int main(int argc, char* argv[])
 #endif
     
   std::cout << "Opening : " <<filename << "\n" << std::endl;
-    // init timer
-    FTic time;
+  // init timer
+  FTic time;
   
+  // Create Matrix Kernel
+  const MatrixKernelClass MatrixKernel;
+
   // init particles position and physical value
   struct TestParticle{
     FPoint position;
@@ -141,7 +144,7 @@ int main(int argc, char* argv[])
   { // -----------------------------------------------------
     std::cout << "\nChebyshev FMM (ORDER="<< ORDER << ",EPS="<< epsilon <<") ... " << std::endl;
     time.tic();
-    KernelClass kernels(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox(), epsilon);
+    KernelClass kernels(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox(), &MatrixKernel, epsilon);
     FmmClass algorithm(app.global(),&tree, &kernels);
     algorithm.execute();
     time.tac();
diff --git a/Tests/noDist/testCompareKernels.cpp b/Tests/noDist/testCompareKernels.cpp
index 5751f087e0226ac2376644ad9aeed17c3d0fc59b..dad4f39a5e2b4ef9426cf85bd75b651df2b216d3 100755
--- a/Tests/noDist/testCompareKernels.cpp
+++ b/Tests/noDist/testCompareKernels.cpp
@@ -26,7 +26,7 @@
 #include "../../Src/Utils/FTic.hpp"
 #include "../../Src/Utils/FParameters.hpp"
 
-#include "../../Src/Files/FFmaScanfLoader.hpp"
+#include "../../Src/Files/FFmaGenericLoader.hpp"
 
 #include "../../Src/Containers/FOctree.hpp"
 #include "../../Src/Containers/FVector.hpp"
@@ -84,7 +84,7 @@ int main(int argc, char* argv[])
         FReal potential;
     };
     // open particle file
-    FFmaScanfLoader loader(filename);
+    FFmaGenericLoader loader(filename);
     if(!loader.isOpen()) throw std::runtime_error("Particle file couldn't be opened!");
 
     TestParticle* const particles = new TestParticle[loader.getNumberOfParticles()];
@@ -156,7 +156,7 @@ int main(int argc, char* argv[])
         { // -----------------------------------------------------
             std::cout << "\nChebyshev FMM ... " << std::endl;
             time.tic();
-            KernelClass kernels(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox(), epsilon);
+            KernelClass kernels(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox(), &MatrixKernel, epsilon);
             FmmClass algorithm(&tree, &kernels);
             algorithm.execute();
             time.tac();
diff --git a/UTests/FUKernelTester.hpp b/UTests/FUKernelTester.hpp
index 9bb0779dc8c62a76225f49fbcd964b2cbf626c48..1c1e0015622b16df47140946b063ba07f8d6e929 100644
--- a/UTests/FUKernelTester.hpp
+++ b/UTests/FUKernelTester.hpp
@@ -48,9 +48,9 @@ public:
     using FUTester<TestClass>::uassert;
 
     // The run function is performing the test for the given configuration
-    template <class CellClass, class ContainerClass, class KernelClass,
+    template <class CellClass, class ContainerClass, class KernelClass, class MatrixKernelClass,
               class LeafClass, class OctreeClass, class FmmClass>
-    void RunTest(std::function<std::unique_ptr<KernelClass>(int NbLevels, FReal boxWidth, FPoint centerOfBox)> GetKernelFunc)	{
+    void RunTest(std::function<std::unique_ptr<KernelClass>(int NbLevels, FReal boxWidth, FPoint centerOfBox, const MatrixKernelClass *const MatrixKernel)> GetKernelFunc)	{
         //
         // Load particles
         //
@@ -85,11 +85,14 @@ public:
             tree.insert(particles[idxPart].getPosition() , idxPart, particles[idxPart].getPhysicalValue() );
         }
         //
+        // Create Matrix Kernel
+        const MatrixKernelClass MatrixKernel; // FUKernelTester is only designed to work with 1/R, i.e. matrix kernel ctor takes no argument.
+        //
         /////////////////////////////////////////////////////////////////////////////////////////////////
         // Run FMM computation
         /////////////////////////////////////////////////////////////////////////////////////////////////
         Print("Fmm...");
-        std::unique_ptr<KernelClass> kernels(GetKernelFunc(NbLevels, loader.getBoxWidth(), loader.getCenterOfBox()));
+        std::unique_ptr<KernelClass> kernels(GetKernelFunc(NbLevels, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel));
         FmmClass algo(&tree,kernels.get());
         algo.execute();
         //
diff --git a/UTests/utestChebyshev.cpp b/UTests/utestChebyshev.cpp
index 1fd76ade868d346edc4dcec3609a92405392ff0b..7a5bcbea18698a82e33f9cbce90a29fb633a1b90 100755
--- a/UTests/utestChebyshev.cpp
+++ b/UTests/utestChebyshev.cpp
@@ -53,9 +53,9 @@ class TestChebyshevDirect : public FUKernelTester<TestChebyshevDirect> {
 		typedef FChebKernel<CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
 		typedef FFmmAlgorithm<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
 		// run test
-        RunTest<CellClass,ContainerClass,KernelClass,LeafClass,OctreeClass,FmmClass>(
-                    [&](int NbLevels, FReal boxWidth, FPoint centerOfBox){
-            return std::unique_ptr<KernelClass>(new KernelClass(NbLevels, boxWidth, centerOfBox));
+    RunTest<CellClass,ContainerClass,KernelClass,MatrixKernelClass,LeafClass,OctreeClass,FmmClass>(
+                                                                                     [&](int NbLevels, FReal boxWidth, FPoint centerOfBox, const MatrixKernelClass *const MatrixKernel){
+                                                                                       return std::unique_ptr<KernelClass>(new KernelClass(NbLevels, boxWidth, centerOfBox, MatrixKernel));
         });
 	}
 
@@ -70,9 +70,9 @@ class TestChebyshevDirect : public FUKernelTester<TestChebyshevDirect> {
 		typedef FChebSymKernel<CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
 		typedef FFmmAlgorithm<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
 		// run test
-        RunTest<CellClass,ContainerClass,KernelClass,LeafClass,OctreeClass,FmmClass>(
-                    [&](int NbLevels, FReal boxWidth, FPoint centerOfBox){
-            return std::unique_ptr<KernelClass>(new KernelClass(NbLevels, boxWidth, centerOfBox));
+    RunTest<CellClass,ContainerClass,KernelClass,MatrixKernelClass,LeafClass,OctreeClass,FmmClass>(
+                    [&](int NbLevels, FReal boxWidth, FPoint centerOfBox, const MatrixKernelClass *const MatrixKernel){
+            return std::unique_ptr<KernelClass>(new KernelClass(NbLevels, boxWidth, centerOfBox, MatrixKernel));
         });
 	}
 
diff --git a/UTests/utestChebyshevDirectPeriodic.cpp b/UTests/utestChebyshevDirectPeriodic.cpp
index 170fa84ceb42942bace770357fb92b27ffd0ed58..f7b461cc4a6e89cc3d9a77319ae65e9c96ccedaa 100755
--- a/UTests/utestChebyshevDirectPeriodic.cpp
+++ b/UTests/utestChebyshevDirectPeriodic.cpp
@@ -104,7 +104,7 @@ class TestChebyshevDirect : public FUTester<TestChebyshevDirect> {
 		/////////////////////////////////////////////////////////////////////////////////////////////////
 		Print("Fmm...");
 		FmmClass algo(&tree,PeriodicDeep );
-		KernelClass kernels(algo.extendedTreeHeight(), algo.extendedBoxWidth(), algo.extendedBoxCenter());
+		KernelClass kernels(algo.extendedTreeHeight(), algo.extendedBoxWidth(), algo.extendedBoxCenter(),&MatrixKernel);
 		algo.setKernel(&kernels);
 		algo.execute();
 		/////////////////////////////////////////////////////////////////////////////////////////////////
diff --git a/UTests/utestChebyshevMpi.cpp b/UTests/utestChebyshevMpi.cpp
index 45735a5bdb263a88866c75045e94f7ac410afb09..ec8ceed4e7a4b56f65c18878ad9c872ab219b6ad 100644
--- a/UTests/utestChebyshevMpi.cpp
+++ b/UTests/utestChebyshevMpi.cpp
@@ -67,6 +67,9 @@ class TestChebyshevMpiDirect : public FUTesterMpi<TestChebyshevMpiDirect>{
     const int nbLevels = 4; 
     const int sizeOfSubLevel = 2;
 
+   // Create Matrix Kernel
+    const MatrixKernelClass MatrixKernel; // FUKernelTester is only designed to work with 1/R, i.e. matrix kernel ctor takes no argument.
+
     // Create octree
     
     struct TestParticle : public FmaRWParticle<8,8>{
@@ -107,7 +110,7 @@ class TestChebyshevMpiDirect : public FUTesterMpi<TestChebyshevMpiDirect>{
     }
     
     
-    KernelClass* kernels= new KernelClass(nbLevels, loader.getBoxWidth(), loader.getCenterOfBox());
+    KernelClass* kernels= new KernelClass(nbLevels, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel);
     FmmClassProc algorithm(app.global(),&tree, kernels);
     algorithm.execute();
     
diff --git a/UTests/utestChebyshevMultiRhs.cpp b/UTests/utestChebyshevMultiRhs.cpp
index d2f52657b67c6e52cac8e77fd4734e0c9b10f9a6..731e0dc9f41f5244c5cda1c48b2d285b575cd83a 100644
--- a/UTests/utestChebyshevMultiRhs.cpp
+++ b/UTests/utestChebyshevMultiRhs.cpp
@@ -83,6 +83,8 @@ class TestChebyshevDirect : public FUTester<TestChebyshevDirect> {
         const int NbLevels        = 4;
         const int SizeSubLevels = 2;
 
+        // Create Matrix Kernel
+        const MatrixKernelClass MatrixKernel; // FUKernelTester is only designed to work with 1/R, i.e. matrix kernel ctor takes no argument.
         //
 		FSize nbParticles = loader.getNumberOfParticles() ;
 		FmaRWParticle<8,8>* const particles = new FmaRWParticle<8,8>[nbParticles];
@@ -116,7 +118,7 @@ class TestChebyshevDirect : public FUTester<TestChebyshevDirect> {
 
         // Run FMM
         Print("Fmm...");
-        KernelClass kernels(NbLevels, loader.getBoxWidth(), loader.getCenterOfBox());
+        KernelClass kernels(NbLevels, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel);
         FmmClass algo(&tree,&kernels);
         algo.execute();
 //
diff --git a/UTests/utestChebyshevThread.cpp b/UTests/utestChebyshevThread.cpp
index 1c929b71dd709d40741d399413e4e542a959ada0..124965c7e0f44a6fdb567dece5b5ce00c17beb36 100755
--- a/UTests/utestChebyshevThread.cpp
+++ b/UTests/utestChebyshevThread.cpp
@@ -81,16 +81,20 @@ class TestChebyshevDirect : public FUTester<TestChebyshevDirect> {
 		const int NbLevels        = 4;
 		const int SizeSubLevels = 2;
 
-		// Create octree
+    // Create Matrix Kernel
+    const MatrixKernelClass MatrixKernel; // FUKernelTester is only designed to work with 1/R, i.e. matrix kernel ctor takes no argument.
 
+    // Load particles
 		FSize nbParticles = loader.getNumberOfParticles() ;
 		FmaRWParticle<8,8>* const particles = new FmaRWParticle<8,8>[nbParticles];
 
 		loader.fillParticle(particles,nbParticles);
-         //
+
+		// Create octree
+		OctreeClass tree(NbLevels, SizeSubLevels, loader.getBoxWidth(), loader.getCenterOfBox());
+    //
 		//   Insert particle in the tree
 		//
-		OctreeClass tree(NbLevels, SizeSubLevels, loader.getBoxWidth(), loader.getCenterOfBox());
 		for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
 		    tree.insert(particles[idxPart].getPosition() , idxPart, particles[idxPart].getPhysicalValue() );
 		}
@@ -99,7 +103,7 @@ class TestChebyshevDirect : public FUTester<TestChebyshevDirect> {
 		// Run FMM computation
 		/////////////////////////////////////////////////////////////////////////////////////////////////
 		Print("Fmm...");
-		KernelClass kernels(NbLevels, loader.getBoxWidth(), loader.getCenterOfBox());
+		KernelClass kernels(NbLevels, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel);
 		FmmClass algo(&tree,&kernels);
 		algo.execute();
 		//0
diff --git a/UTests/utestLagrange.cpp b/UTests/utestLagrange.cpp
index 93d42880947fae07eb572a876ace8f19d7cc81bc..d83d95133020c9bf0c1c36cc00bedea36b7f71ce 100755
--- a/UTests/utestLagrange.cpp
+++ b/UTests/utestLagrange.cpp
@@ -73,16 +73,20 @@ class TestLagrange : public FUTester<TestLagrange> {
 		const int NbLevels        = 4;
 		const int SizeSubLevels = 2;
 
-		// Create octree
+    // Create Matrix Kernel
+    const MatrixKernelClass MatrixKernel; // FUKernelTester is only designed to work with 1/R, i.e. matrix kernel ctor takes no argument.
 
+    // Load particles
 		FSize nbParticles = loader.getNumberOfParticles() ;
 		FmaRWParticle<8,8>* const particles = new FmaRWParticle<8,8>[nbParticles];
 
 		loader.fillParticle(particles,nbParticles);
-         //
+
+		// Create octree
+		OctreeClass tree(NbLevels, SizeSubLevels, loader.getBoxWidth(), loader.getCenterOfBox());
+    //
 		//   Insert particle in the tree
 		//
-		OctreeClass tree(NbLevels, SizeSubLevels, loader.getBoxWidth(), loader.getCenterOfBox());
 		for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
 		    tree.insert(particles[idxPart].getPosition(), idxPart, particles[idxPart].getPhysicalValue() );
 		}
@@ -91,7 +95,7 @@ class TestLagrange : public FUTester<TestLagrange> {
 		// Run FMM computation
 		/////////////////////////////////////////////////////////////////////////////////////////////////
 		Print("Fmm...");
-		KernelClass kernels(NbLevels, loader.getBoxWidth(), loader.getCenterOfBox());
+		KernelClass kernels(NbLevels, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel);
 		FmmClass algo(&tree,&kernels);
 		algo.execute();
 		//0
diff --git a/UTests/utestLagrangeMpi.cpp b/UTests/utestLagrangeMpi.cpp
index 569f237e8658380f07ca25d653dacec10b63353d..931412bbfbed408ffe9cb7b7ef6958244b7df785 100644
--- a/UTests/utestLagrangeMpi.cpp
+++ b/UTests/utestLagrangeMpi.cpp
@@ -66,6 +66,9 @@ class TestLagrangeMpiDirect : public FUTesterMpi<TestLagrangeMpiDirect>{
     const int nbLevels = 4; 
     const int sizeOfSubLevel = 2;
 
+    // Create Matrix Kernel
+    const MatrixKernelClass MatrixKernel; // FUKernelTester is only designed to work with 1/R, i.e. matrix kernel ctor takes no argument.
+
     // Create octree
     
     struct TestParticle : public FmaRWParticle<8,8>{
@@ -106,7 +109,7 @@ class TestLagrangeMpiDirect : public FUTesterMpi<TestLagrangeMpiDirect>{
     }
     
     
-    KernelClass* kernels= new KernelClass(nbLevels, loader.getBoxWidth(), loader.getCenterOfBox());
+    KernelClass* kernels= new KernelClass(nbLevels, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel);
     FmmClassProc algorithm(app.global(),&tree, kernels);
     algorithm.execute();
     
diff --git a/UTests/utestLagrangeThread.cpp b/UTests/utestLagrangeThread.cpp
index 45fe8bb484bbc33a88a64cf3602f6ef0f2deeb40..c5ee30dad5013d591d500902708684f724845a29 100755
--- a/UTests/utestLagrangeThread.cpp
+++ b/UTests/utestLagrangeThread.cpp
@@ -81,16 +81,19 @@ class TestLagrange : public FUTester<TestLagrange> {
 		const int NbLevels        = 4;
 		const int SizeSubLevels = 2;
 
-		// Create octree
+    // Create Matrix Kernel
+    const MatrixKernelClass MatrixKernel; // FUKernelTester is only designed to work with 1/R, i.e. matrix kernel ctor takes no argument.
 
 		FSize nbParticles = loader.getNumberOfParticles() ;
 		FmaRWParticle<8,8>* const particles = new FmaRWParticle<8,8>[nbParticles];
 
 		loader.fillParticle(particles,nbParticles);
-         //
+
+		// Create octree
+		OctreeClass tree(NbLevels, SizeSubLevels, loader.getBoxWidth(), loader.getCenterOfBox());
+    //
 		//   Insert particle in the tree
 		//
-		OctreeClass tree(NbLevels, SizeSubLevels, loader.getBoxWidth(), loader.getCenterOfBox());
 		for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
 		    tree.insert(particles[idxPart].getPosition() , idxPart, particles[idxPart].getPhysicalValue() );
 		}
@@ -99,7 +102,7 @@ class TestLagrange : public FUTester<TestLagrange> {
 		// Run FMM computation
 		/////////////////////////////////////////////////////////////////////////////////////////////////
 		Print("Fmm...");
-		KernelClass kernels(NbLevels, loader.getBoxWidth(), loader.getCenterOfBox());
+		KernelClass kernels(NbLevels, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel);
 		FmmClass algo(&tree,&kernels);
 		algo.execute();
 		//0