diff --git a/Src/Kernels/Uniform/FUnifM2LHandler.hpp b/Src/Kernels/Uniform/FUnifM2LHandler.hpp
index b550a1c9fbba28086af46429031b4ed49fda8ddb..e2cf8063cfed7c672e626a99bf12768594ac2b7b 100644
--- a/Src/Kernels/Uniform/FUnifM2LHandler.hpp
+++ b/Src/Kernels/Uniform/FUnifM2LHandler.hpp
@@ -66,8 +66,10 @@ static void Compute(const MatrixKernelClass *const MatrixKernel, const FReal Cel
   TensorType::setNodeIdsPairs(node_ids_pairs);
 
   // init Discrete Fourier Transformator
+  const int dimfft = 1; // unidim FFT since fully circulant embedding
+  const int steps[dimfft] = {rc};
 //	FDft Dft(rc);
-	FFft<1> Dft(rc);
+	FFft<dimfft> Dft(steps);
 
   // get first column of K via permutation
   unsigned int perm[rc];
@@ -177,15 +179,17 @@ class FUnifM2LHandler<ORDER,HOMOGENEOUS> : FNoCopyable
 
 	FComplexe *FC;
 
-  // for real valued kernel only n/2+1 complex values are stored 
-  // after performing the DFT (the rest is deduced by conjugation)
-  unsigned int opt_rc; 
-
   typedef FUnifTensor<ORDER> TensorType;
   unsigned int node_diff[nnodes*nnodes];
 
-  //  FDft Dft; // Direct Discrete Fourier Transformator
-  FFft<1> Dft; // Fast Discrete Fourier Transformator
+  /// DFT specific
+  static const int dimfft = 1; // unidim FFT since fully circulant embedding
+//  FDft Dft; // Direct Discrete Fourier Transformator
+  typedef FFft<dimfft> DftClass; // Fast Discrete Fourier Transformator
+  FSmartPointer<DftClass,FSmartPointerMemory> Dft;
+
+  const unsigned int opt_rc; // specific to real valued kernel
+
 
 	static const std::string getFileName()
 	{
@@ -201,9 +205,12 @@ public:
 	template <typename MatrixKernelClass>
 	FUnifM2LHandler(const MatrixKernelClass *const MatrixKernel, const unsigned int, const FReal)
 		: FC(NULL),
-      opt_rc(rc/2+1),
-      Dft(rc) // initialize Discrete Fourier Transformator
+      opt_rc(rc/2+1)
 	{    
+    // init DFT
+    const int steps[dimfft] = {rc};
+    Dft = new DftClass(steps);
+
     // initialize root node ids
     TensorType::setNodeIdsDiff(node_diff);
 
@@ -253,7 +260,7 @@ public:
     FReal Px[rc];
     FBlas::setzero(rc,Px);
     // Apply forward Discrete Fourier Transform
-    Dft.applyIDFT(FX,Px);
+    Dft->applyIDFT(FX,Px);
 
     // Unapply Zero Padding
     for (unsigned int j=0; j<nnodes; ++j)
@@ -319,7 +326,7 @@ public:
       Py[node_diff[i*nnodes]]=y[i];
 
     // Apply forward Discrete Fourier Transform
-    Dft.applyDFT(Py,FY);
+    Dft->applyDFT(Py,FY);
 
   }
 
@@ -342,15 +349,17 @@ class FUnifM2LHandler<ORDER,NON_HOMOGENEOUS> : FNoCopyable
   const unsigned int TreeHeight;
   const FReal RootCellWidth;
 
-  // for real valued kernel only n/2+1 complex values are stored 
-  // after performing the DFT (the rest is deduced by conjugation)
-  unsigned int opt_rc; 
-
   typedef FUnifTensor<ORDER> TensorType;
   unsigned int node_diff[nnodes*nnodes];
 
-  //  FDft Dft; // Direct Discrete Fourier Transformator
-  FFft<1> Dft; // Fast Discrete Fourier Transformator
+  // DFT specific
+  static const int dimfft = 1; // unidim FFT since fully circulant embedding
+//  FDft Dft; // Direct Discrete Fourier Transformator
+  typedef FFft<dimfft> DftClass; // Fast Discrete Fourier Transformator
+  FSmartPointer<DftClass,FSmartPointerMemory> Dft;
+
+  const unsigned int opt_rc; // specific to real valued kernel
+
 
 	static const std::string getFileName()
 	{
@@ -367,9 +376,12 @@ public:
 	FUnifM2LHandler(const MatrixKernelClass *const MatrixKernel, const unsigned int inTreeHeight, const FReal inRootCellWidth)
 		: TreeHeight(inTreeHeight),
       RootCellWidth(inRootCellWidth),
-      opt_rc(rc/2+1),
-      Dft(rc) // initialize Discrete Fourier Transformator
-	{    
+      opt_rc(rc/2+1)
+	{
+    // init DFT
+    const int steps[dimfft] = {rc};
+    Dft = new DftClass(steps);
+
     // initialize root node ids
     TensorType::setNodeIdsDiff(node_diff);
     
@@ -431,7 +443,7 @@ public:
     FReal Px[rc];
     FBlas::setzero(rc,Px);
     // Apply forward Discrete Fourier Transform
-    Dft.applyIDFT(FX,Px);
+    Dft->applyIDFT(FX,Px);
 
     // Unapply Zero Padding
     for (unsigned int j=0; j<nnodes; ++j)
@@ -496,7 +508,7 @@ public:
       Py[node_diff[i*nnodes]]=y[i];
 
     // Apply forward Discrete Fourier Transform
-    Dft.applyDFT(Py,FY);
+    Dft->applyDFT(Py,FY);
 
   }
 
diff --git a/Src/Kernels/Uniform/FUnifSymM2LHandler.hpp b/Src/Kernels/Uniform/FUnifSymM2LHandler.hpp
index a7e4843b651eadebe7ae35815815213ffee0f982..4fa1f77011420e34eb29a64679c63ec99d1cca6c 100755
--- a/Src/Kernels/Uniform/FUnifSymM2LHandler.hpp
+++ b/Src/Kernels/Uniform/FUnifSymM2LHandler.hpp
@@ -60,8 +60,10 @@ static void precompute(const MatrixKernelClass *const MatrixKernel, const FReal
 	_FC = new FComplexe [rc]; // TODO: do it in the non-sym version!!!
 
   // init Discrete Fourier Transformator
+  const int dimfft = 1; // unidim FFT since fully circulant embedding
+  const int steps[dimfft] = {rc};
 //	FDft Dft(rc);
-	FFft Dft(rc);
+	FFft<dimfft> Dft(steps);
 
   // reduce storage if real valued kernel
   const unsigned int opt_rc = rc/2+1;
diff --git a/Src/Kernels/Uniform/FUnifTensorialKernel.hpp b/Src/Kernels/Uniform/FUnifTensorialKernel.hpp
index 74f1cc760cb78ea1e3ec7c2b3fd58b5d95fea000..5734d70a0dec1b2a4d5b7472251b5536c5484b7b 100644
--- a/Src/Kernels/Uniform/FUnifTensorialKernel.hpp
+++ b/Src/Kernels/Uniform/FUnifTensorialKernel.hpp
@@ -48,9 +48,9 @@ class FTreeCoordinate;
  * In fact, in the ChebyshevSym variant the matrix kernel needs to be 
  * evaluated compo-by-compo since we currently use a scalar ACA.
  *
- * 3) We currently use multiple 1D FFT instead of multidim FFT. 
- * TODO switch to multidim if relevant in considered range of size 
- * (see testFFTW and testFFTWMultidim).
+ * 3) We currently use multiple 1D FFT instead of multidim FFT since embedding
+ * is circulant. Multidim FFT could be used if embedding were block circulant.
+ * TODO investigate possibility of block circulant embedding
  *
  * @tparam CellClass Type of cell
  * @tparam ContainerClass Type of container to store particles
diff --git a/Src/Kernels/Uniform/FUnifTensorialM2LHandler.hpp b/Src/Kernels/Uniform/FUnifTensorialM2LHandler.hpp
index a4c1d62db8ca62ab5873fc2de13ba452ef99fe2b..84c9dfa68a0bfba21d2c9b3256fc3d4b5b252345 100644
--- a/Src/Kernels/Uniform/FUnifTensorialM2LHandler.hpp
+++ b/Src/Kernels/Uniform/FUnifTensorialM2LHandler.hpp
@@ -75,8 +75,10 @@ static void Compute(const MatrixKernelClass *const MatrixKernel,
   TensorType::setNodeIdsPairs(node_ids_pairs);
 
   // init Discrete Fourier Transformator
+  const int dimfft = 1; // unidim FFT since fully circulant embedding
+  const int steps[dimfft] = {rc};
 //	FDft Dft(rc);
-	FFft<1> Dft(rc);
+	FFft<dimfft> Dft(steps);
 
   // get first column of K via permutation
   unsigned int perm[rc];
@@ -198,15 +200,17 @@ class FUnifTensorialM2LHandler<ORDER,MatrixKernelClass,HOMOGENEOUS> : FNoCopyabl
   // Tensorial MatrixKernel specific
 	FComplexe** FC;
 
-  // for real valued kernel only n/2+1 complex values are stored 
-  // after performing the DFT (the rest is deduced by conjugation)
-  unsigned int opt_rc; 
-
   typedef FUnifTensor<ORDER> TensorType;
   unsigned int node_diff[nnodes*nnodes];
 
-  //  FDft Dft; // Direct Discrete Fourier Transformator
-  FFft<1> Dft; // Fast Discrete Fourier Transformator
+  // DFT specific
+  static const int dimfft = 1; // unidim FFT since fully circulant embedding
+//  FDft Dft; // Direct Discrete Fourier Transformator
+  typedef FFft<dimfft> DftClass; // Fast Discrete Fourier Transformator
+  FSmartPointer<DftClass,FSmartPointerMemory> Dft;
+
+  const unsigned int opt_rc; // specific to real valued kernel
+
 
 	static const std::string getFileName()
 	{
@@ -220,9 +224,12 @@ class FUnifTensorialM2LHandler<ORDER,MatrixKernelClass,HOMOGENEOUS> : FNoCopyabl
 	
 public:
 	FUnifTensorialM2LHandler(const MatrixKernelClass *const MatrixKernel, const unsigned int, const FReal)
-		: opt_rc(rc/2+1), 
-      Dft(rc) // initialize Discrete Fourier Transformator
+		: opt_rc(rc/2+1)
 	{
+    // init DFT
+    const int steps[dimfft] = {rc};
+    Dft = new DftClass(steps);
+
     // allocate FC
     FC = new FComplexe*[dim];
     for (unsigned int d=0; d<dim; ++d)
@@ -277,7 +284,7 @@ public:
     FReal Px[rc];
     FBlas::setzero(rc,Px);
     // Apply forward Discrete Fourier Transform
-    Dft.applyIDFT(FX,Px);
+    Dft->applyIDFT(FX,Px);
 
     // Unapply Zero Padding
     for (unsigned int j=0; j<nnodes; ++j)
@@ -344,7 +351,7 @@ public:
       Py[node_diff[i*nnodes]]=y[i];
 
     // Apply forward Discrete Fourier Transform
-    Dft.applyDFT(Py,FY);
+    Dft->applyDFT(Py,FY);
 
   }
 
@@ -369,15 +376,17 @@ class FUnifTensorialM2LHandler<ORDER,MatrixKernelClass,NON_HOMOGENEOUS> : FNoCop
   const unsigned int TreeHeight;
   const FReal RootCellWidth;
 
-  // for real valued kernel only n/2+1 complex values are stored 
-  // after performing the DFT (the rest is deduced by conjugation)
-  unsigned int opt_rc; 
-
   typedef FUnifTensor<ORDER> TensorType;
   unsigned int node_diff[nnodes*nnodes];
 
-  //  FDft Dft; // Direct Discrete Fourier Transformator
-  FFft<1> Dft; // Fast Discrete Fourier Transformator
+  // DFT specific
+  static const int dimfft = 1; // unidim FFT since fully circulant embedding
+//  FDft Dft; // Direct Discrete Fourier Transformator
+  typedef FFft<dimfft> DftClass; // Fast Discrete Fourier Transformator
+  FSmartPointer<DftClass,FSmartPointerMemory> Dft;
+
+  const unsigned int opt_rc; // specific to real valued kernel
+
 
 	static const std::string getFileName()
 	{
@@ -393,9 +402,12 @@ public:
 	FUnifTensorialM2LHandler(const MatrixKernelClass *const MatrixKernel, const unsigned int inTreeHeight, const FReal inRootCellWidth)
 		: TreeHeight(inTreeHeight),
       RootCellWidth(inRootCellWidth),
-      opt_rc(rc/2+1), 
-      Dft(rc) // initialize Discrete Fourier Transformator
+      opt_rc(rc/2+1)
 	{
+    // init DFT
+    const int steps[dimfft] = {rc};
+    Dft = new DftClass(steps);
+
     // allocate FC
     FC = new FComplexe**[TreeHeight];
 		for (unsigned int l=0; l<TreeHeight; ++l){ 
@@ -458,7 +470,7 @@ public:
     FReal Px[rc];
     FBlas::setzero(rc,Px);
     // Apply forward Discrete Fourier Transform
-    Dft.applyIDFT(FX,Px);
+    Dft->applyIDFT(FX,Px);
 
     // Unapply Zero Padding
     for (unsigned int j=0; j<nnodes; ++j)
@@ -524,7 +536,7 @@ public:
       Py[node_diff[i*nnodes]]=y[i];
 
     // Apply forward Discrete Fourier Transform
-    Dft.applyDFT(Py,FY);
+    Dft->applyDFT(Py,FY);
 
   }
 
diff --git a/Src/Utils/FBlas.hpp b/Src/Utils/FBlas.hpp
index 01fea5dd01a4e61e9ce93e145a29955ab127eae2..fe1731a8958ab8d7f5fbcc6e22c8cac605915f2c 100755
--- a/Src/Utils/FBlas.hpp
+++ b/Src/Utils/FBlas.hpp
@@ -70,6 +70,7 @@ extern "C"
 							 double*, double*, const unsigned*, int*);
 	void dorgqr_(const unsigned*, const unsigned*, const unsigned*,
 							 double*, const unsigned*, double*, double*, const unsigned*, int*);
+  void dpotrf_(const char*, const unsigned*, double*, const unsigned*, int*);
 
 #ifdef ScalFMM_USE_MKL_AS_BLAS
   // mkl: hadamard product is not implemented in mkl_blas
@@ -103,6 +104,7 @@ extern "C"
 							 float*, float*, const unsigned*, int*);
 	void sorgqr_(const unsigned*, const unsigned*, const unsigned*,
 							 float*, const unsigned*, float*, float*, const unsigned*, int*);
+  void spotrf_(const char*, const unsigned*, float*, const unsigned*, int*);
 
 	// double complex //////////////////////////////////////////////////
 	// blas 1
@@ -119,7 +121,7 @@ extern "C"
 							double*, const unsigned*, const double*, double*, const unsigned*);
 	void zgeqrf_(const unsigned*, const unsigned*, double*, const unsigned*,
 							 double*, double*, const unsigned*, int*);
-
+  void zpotrf_(const char*, const unsigned*, double*, const unsigned*, int*);
 
 	// single complex //////////////////////////////////////////////////
 	// blas 1
@@ -136,7 +138,7 @@ extern "C"
 							float*, const unsigned*, const float*, float*, const unsigned*);
 	void cgeqrf_(const unsigned*, const unsigned*, float*, const unsigned*,
 							 float*, float*, const unsigned*, int*);
-
+  void cpotrf_(const char*, const unsigned*, float*, const unsigned*, int*);
 
 
 }
@@ -506,6 +508,19 @@ namespace FBlas {
 
 
 
+  // Cholesky decomposition: A=LL^T (if A is symmetric definite positive)
+  inline int potrf(const unsigned m, double* A, const unsigned n)
+  { 
+		int INF;  
+    dpotrf_("L", &m, A, &n, &INF); 
+    return INF;
+  }
+  inline int potrf(const unsigned m, float* A, const unsigned n)
+  { 
+		int INF;  
+    spotrf_("L", &m, A, &n, &INF); 
+    return INF;
+  }
 
 } // end namespace FCBlas
 
diff --git a/Src/Utils/FDft.hpp b/Src/Utils/FDft.hpp
index ce9c0b86b73febacafca80e8a8a74f74abb42e38..7d81afed087c1f2d0977d20161c0acf32d07f0cb 100644
--- a/Src/Utils/FDft.hpp
+++ b/Src/Utils/FDft.hpp
@@ -36,7 +36,7 @@
 
 /**
  * @author Pierre Blanchard (pierre.blanchard@inria.fr)
- * @class FDft and @class FFft
+ * @class FDft, @class FFft and @class FCFft
  * Please read the license
  *
  * These classes handle the forward and backward Discete Fourier Transform 
@@ -45,6 +45,10 @@
  * Fourier Transform (FFT). The FFT algorithm can either be provided by the 
  * FFTW(3) library itself or a version that is wrapped in Intel MKL.
  *
+ * The direct DFT is templatized with the input value type (FReal or FComplexe),
+ * while 2 distinct classes resp. @class FFft and @class FCFft are defined 
+ * resp. for the real and complex valued FFT.
+ *
  * The aim of writing these specific classes is to allow further customization 
  * of the DFT such as implementing a tensorial variant, a weighted variant 
  * or any other feature.
@@ -59,25 +63,25 @@
  *
  * @tparam nsteps number of sampled values \f$N\f$
  */
-
+template<typename ValueType>
 class FDft
 {
 
-  FReal* fftR_;
-  FComplexe* fftC_;
+  ValueType* data_;  //< data in physical space
+  FComplexe* dataF_; //< data in Fourier space
 
   FReal *cosRS_, *sinRS_; 
 
 private:
-  unsigned int nsteps_;
+  unsigned int nsteps_; //< number of steps
 
 public:
 
   FDft(const unsigned int nsteps)
   : nsteps_(nsteps)
   {
-    fftR_ = new FReal[nsteps_];
-    fftC_ = new FComplexe[nsteps_];
+    data_  = new ValueType[nsteps_];
+    dataF_ = new FComplexe[nsteps_];
 
     // Beware this is extremely HEAVY to store!!! => Use FDft only for debug!
     cosRS_ = new FReal[nsteps_*nsteps_];
@@ -93,67 +97,104 @@ public:
 
   virtual ~FDft()
   {
-    delete [] fftR_;
-    delete [] fftC_;
+    delete [] data_;
+    delete [] dataF_;
     delete [] cosRS_;
     delete [] sinRS_;
   }
 
+  /// Forward DFT
+  // Real valued DFT
   void applyDFT(const FReal* sampledData,
                 FComplexe* transformedData) const
   {
-//    FTic time;
-
     // read sampled data
-//    std::cout<< "copy(";
-//    time.tic();
-    FBlas::c_setzero(nsteps_,reinterpret_cast<FReal*>(fftC_));
-    FBlas::copy(nsteps_, sampledData,fftR_);
-//    std::cout << time.tacAndElapsed() << ")";
-
-//    std::cout<< " - exe(";
-//    time.tic();
+    FBlas::c_setzero(nsteps_,reinterpret_cast<FReal*>(dataF_));
+    FBlas::copy(nsteps_, sampledData,data_);
+
+    // perform direct forward transformation
     for(unsigned int r=0; r<nsteps_; ++r)
       for(unsigned int s=0; s<nsteps_; ++s){
-        fftC_[r] += FComplexe(fftR_[s]*cosRS_[r*nsteps_+s], 
-                              -fftR_[s]*sinRS_[r*nsteps_+s]);
+        dataF_[r] += FComplexe(data_[s]*cosRS_[r*nsteps_+s], 
+                              -data_[s]*sinRS_[r*nsteps_+s]);
       }
-//    std::cout << time.tacAndElapsed() << ")";
 
     // write transformed data
-//    std::cout<< " - copy(";
-//    time.tic();
-    FBlas::c_copy(nsteps_,reinterpret_cast<FReal*>(fftC_),
+    FBlas::c_copy(nsteps_,reinterpret_cast<FReal*>(dataF_),
                   reinterpret_cast<FReal*>(transformedData));
-//    std::cout << time.tacAndElapsed() << ") ";
-
   }
+  // Complexe valued DFT
+  void applyDFT(const FComplexe* sampledData,
+                FComplexe* transformedData) const
+  {
+    // read sampled data
+    FBlas::c_setzero(nsteps_,reinterpret_cast<FReal*>(dataF_));
+    FBlas::c_copy(nsteps_,reinterpret_cast<const FReal*>(sampledData),
+                  reinterpret_cast<FReal*>(data_));
+
+    // perform direct forward transformation
+    for(unsigned int r=0; r<nsteps_; ++r)
+      for(unsigned int s=0; s<nsteps_; ++s){
+        dataF_[r] += FComplexe(data_[s].getReal()*cosRS_[r*nsteps_+s]
+                               + data_[s].getImag()*sinRS_[r*nsteps_+s], 
+                               data_[s].getImag()*cosRS_[r*nsteps_+s]
+                               - data_[s].getReal()*sinRS_[r*nsteps_+s]);
+      }
 
+    // write transformed data
+    FBlas::c_copy(nsteps_,reinterpret_cast<FReal*>(dataF_),
+                  reinterpret_cast<FReal*>(transformedData));
+  }
 
+  /// Backward DFT
+  // Real valued IDFT
   void applyIDFT(const FComplexe* transformedData,
                  FReal* sampledData) const
   {
     // read transformed data
-    FBlas::setzero(nsteps_,fftR_);
+    FBlas::setzero(nsteps_,data_);
     FBlas::c_copy(nsteps_,reinterpret_cast<const FReal*>(transformedData),
-                  reinterpret_cast<FReal*>(fftC_));
+                  reinterpret_cast<FReal*>(dataF_));
 
+    // perform direct backward transformation
     for(unsigned int r=0; r<nsteps_; ++r){
       for(unsigned int s=0; s<nsteps_; ++s){
-        fftR_[r] += fftC_[s].getReal()*cosRS_[r*nsteps_+s]
-          -fftC_[s].getImag()*sinRS_[r*nsteps_+s];
+        data_[r] += dataF_[s].getReal()*cosRS_[r*nsteps_+s]
+          + dataF_[s].getImag()*sinRS_[r*nsteps_+s];
       }
-      fftR_[r]*=1./nsteps_;
+      data_[r]*=1./nsteps_;
     }  
 
     // write sampled data
-    FBlas::copy(nsteps_,fftR_,sampledData);
-
+    FBlas::copy(nsteps_,data_,sampledData);
   }
 
+  // Complexe valued IDFT
+  void applyIDFT(const FComplexe* transformedData,
+                 FComplexe* sampledData) const
+  {
+    // read transformed data
+    FBlas::c_setzero(nsteps_,reinterpret_cast<FReal*>(data_));
+    FBlas::c_copy(nsteps_,reinterpret_cast<const FReal*>(transformedData),
+                  reinterpret_cast<FReal*>(dataF_));
 
-};
+    // perform direct backward transformation
+    for(unsigned int r=0; r<nsteps_; ++r){
+      for(unsigned int s=0; s<nsteps_; ++s){
+        data_[r] += FComplexe(dataF_[s].getReal()*cosRS_[r*nsteps_+s]
+                              - dataF_[s].getImag()*sinRS_[r*nsteps_+s],
+                              dataF_[s].getImag()*cosRS_[r*nsteps_+s]
+                              + dataF_[s].getReal()*sinRS_[r*nsteps_+s]);
+      }
+      data_[r]*=1./nsteps_;
+    }  
 
+    // write sampled data
+    FBlas::c_copy(nsteps_,reinterpret_cast<FReal*>(data_),
+                  reinterpret_cast<FReal*>(sampledData));
+  }
+
+};
 
 
 /**
@@ -167,50 +208,56 @@ class FFft
   enum{dim = DIM};
 
   // arrays
-  FReal* fftR_;
-  FComplexe* fftC_;
+  FReal* data_;      //< data in physical space
+  FComplexe* dataF_; //< data in Fourier space
 
   // plans
-  fftw_plan plan_c2r_; // backward FFT plan
-  fftw_plan plan_r2c_; // forward FFT plan 
+  fftw_plan plan_c2r_; //< backward FFT plan
+  fftw_plan plan_r2c_; //< forward FFT plan 
 
 private:
-  /*unsigned*/ int nsteps_; // all dimensions have same nb of steps
-  /*unsigned*/ int nsteps_opt_;
+  int steps_[dim]; //< number of steps per dimension
+  int nsteps_;     //< total number of steps
+  int nsteps_opt_; //< reduced number of steps for real valued FFT
 
 public:
 
-  FFft(const unsigned int nsteps)
-    : nsteps_(nsteps),
-      nsteps_opt_(nsteps/2+1) // SPECIFIC TO REAL VALUED FTT
+  FFft(const int steps[dim])
   {
+    // init number of steps
+    nsteps_=1; nsteps_opt_=1;
+    for(int d = 0; d<dim;++d){
+      steps_[d]=steps[d];
+      nsteps_*=steps[d];
+      nsteps_opt_*=steps[d]/2+1; // real valued DFT specific
+    }
+
     // allocate arrays
-    fftR_ = (FReal*) fftw_malloc(dim * sizeof(FReal) * dim * nsteps_);
-    fftC_ = (FComplexe*) fftw_malloc(dim * sizeof(FComplexe) * dim * nsteps_opt_);
+    data_ = (FReal*) fftw_malloc(sizeof(FReal) * nsteps_);
+    dataF_ = (FComplexe*) fftw_malloc(sizeof(FComplexe) * nsteps_opt_);
 
-//    // fftw plans
+//    // unidim fftw plans
 //    plan_c2r_ =
 //      fftw_plan_dft_c2r_1d(nsteps_, 
-//                           reinterpret_cast<fftw_complex*>(fftC_),
-//                           fftR_, 
+//                           reinterpret_cast<fftw_complex*>(dataF_),
+//                           data_, 
 //                           FFTW_MEASURE);// TODO: test FFTW_ESTIMATE
 //    plan_r2c_ =
 //      fftw_plan_dft_r2c_1d(nsteps_, 
-//                           fftR_, 
-//                           reinterpret_cast<fftw_complex*>(fftC_), 
+//                           data_, 
+//                           reinterpret_cast<fftw_complex*>(dataF_), 
 //                           FFTW_MEASURE);
 
     // multidim fftw plans
-    const int steps[2] = {dim, nsteps_};
     plan_c2r_ =
-      fftw_plan_dft_c2r(2, steps, 
-                        reinterpret_cast<fftw_complex*>(fftC_),
-                        fftR_, 
+      fftw_plan_dft_c2r(dim, steps_, 
+                        reinterpret_cast<fftw_complex*>(dataF_),
+                        data_, 
                         FFTW_MEASURE);// TODO: test FFTW_ESTIMATE
     plan_r2c_ =
-      fftw_plan_dft_r2c(2, steps, 
-                        fftR_, 
-                        reinterpret_cast<fftw_complex*>(fftC_), 
+      fftw_plan_dft_r2c(dim, steps_, 
+                        data_, 
+                        reinterpret_cast<fftw_complex*>(dataF_), 
                         FFTW_MEASURE);
 
   }
@@ -219,38 +266,22 @@ public:
   {
     fftw_destroy_plan(plan_c2r_);
     fftw_destroy_plan(plan_r2c_);
-    fftw_free(fftR_);
-    fftw_free(fftC_);
+    fftw_free(data_);
+    fftw_free(dataF_);
   }
 
   void applyDFT(const FReal* sampledData,
                 FComplexe* transformedData) const
   {
-//    FTic time;
-
-    // read sampled data
-//    std::cout<< "copy(";
-//    time.tic();
-    FBlas::c_setzero(dim*nsteps_opt_,reinterpret_cast<FReal*>(fftC_));
-    FBlas::copy(dim*nsteps_, sampledData,fftR_);
-//    std::cout << time.tacAndElapsed() << ")";
+    FBlas::c_setzero(nsteps_opt_,reinterpret_cast<FReal*>(dataF_));
+    FBlas::copy(nsteps_, sampledData,data_);
 
     // perform fft
-//    std::cout<< " - exe(";
-//    time.tic();
     fftw_execute( plan_r2c_ );
-//    std::cout << time.tacAndElapsed() << ")";
 
     // write transformed data
-//    std::cout<< " - copy(";
-//    time.tic();
-    FBlas::c_copy(dim*nsteps_opt_,reinterpret_cast<FReal*>(fftC_),
+    FBlas::c_copy(nsteps_opt_,reinterpret_cast<FReal*>(dataF_),
                   reinterpret_cast<FReal*>(transformedData));
-//    for(unsigned int s=0; s<nsteps_opt_; ++s)
-//      transformedData[s]=fftC_[s];
-
-//    std::cout << time.tacAndElapsed() << ") ";
-
   }
 
 
@@ -258,18 +289,118 @@ public:
                  FReal* sampledData) const
   {
     // read transformed data
-    FBlas::setzero(dim*nsteps_,fftR_);
-    FBlas::c_copy(dim*nsteps_opt_,reinterpret_cast<const FReal*>(transformedData),
-                  reinterpret_cast<FReal*>(fftC_));
+    FBlas::setzero(nsteps_,data_);
+    FBlas::c_copy(nsteps_opt_,reinterpret_cast<const FReal*>(transformedData),
+                  reinterpret_cast<FReal*>(dataF_));
 
     // perform ifft
     fftw_execute( plan_c2r_ );
 
     for(unsigned int s=0; s<dim*nsteps_; ++s)
-      fftR_[s]/=(dim*nsteps_); // the fft from fftw is not scaled !!!!!!!!
+      data_[s]/=nsteps_; // the fft from fftw is not scaled !!!!!!!!
+
+    // write sampled data
+    FBlas::copy(nsteps_,data_,sampledData);
+
+  }
+
+
+};
+
+
+/**
+ * @class FCFft
+ *
+ * @tparam nsteps number of sampled values \f$N\f$
+ */
+template<int DIM = 1>
+class FCFft
+{
+  enum{dim = DIM};
+
+  // arrays
+  FComplexe* data_;  //< data in physical space
+  FComplexe* dataF_; //< data in Fourier space
+
+  // plans
+  fftw_plan plan_b_; // backward FFT plan
+  fftw_plan plan_f_; // forward FFT plan 
+
+private:
+  int steps_[dim]; //< number of steps per dimension
+  int nsteps_;     //< total number of steps 
+public:
+
+  FCFft(const int steps[dim])
+  {
+    // init number of steps
+    nsteps_=1;
+    for(int d = 0; d<dim;++d){
+      steps_[d]=steps[d];
+      nsteps_*=steps[d];
+    }
+
+    // allocate arrays
+    data_  = (FComplexe*) fftw_malloc(sizeof(FComplexe) * nsteps_);
+    dataF_ = (FComplexe*) fftw_malloc(sizeof(FComplexe) * nsteps_);
+
+    // multidim fftw plans
+    plan_b_ =
+      fftw_plan_dft(dim, steps_, 
+                    reinterpret_cast<fftw_complex*>(dataF_),
+                    reinterpret_cast<fftw_complex*>(data_), 
+                    FFTW_BACKWARD,
+                    FFTW_MEASURE);// TODO: test FFTW_ESTIMATE
+    plan_f_ =
+      fftw_plan_dft(dim, steps_, 
+                    reinterpret_cast<fftw_complex*>(data_), 
+                    reinterpret_cast<fftw_complex*>(dataF_), 
+                    FFTW_FORWARD,
+                    FFTW_MEASURE);
+
+  }
+
+  virtual ~FCFft()
+  {
+    fftw_destroy_plan(plan_b_);
+    fftw_destroy_plan(plan_f_);
+    fftw_free(data_);
+    fftw_free(dataF_);
+  }
+
+  void applyDFT(const FComplexe* sampledData,
+                FComplexe* transformedData) const
+  {
+    FBlas::c_setzero(nsteps_,reinterpret_cast<FReal*>(dataF_));
+    FBlas::c_copy(nsteps_, reinterpret_cast<const FReal*>(sampledData),
+                  reinterpret_cast<FReal*>(data_));
+
+    // perform fft
+    fftw_execute( plan_f_ );
+
+    // write transformed data
+    FBlas::c_copy(nsteps_,reinterpret_cast<FReal*>(dataF_),
+                  reinterpret_cast<FReal*>(transformedData));
+  }
+
+
+  void applyIDFT(const FComplexe* transformedData,
+                 FComplexe* sampledData) const
+  {
+    // read transformed data
+    FBlas::c_setzero(nsteps_,reinterpret_cast<FReal*>(data_));
+    FBlas::c_copy(nsteps_,reinterpret_cast<const FReal*>(transformedData),
+                  reinterpret_cast<FReal*>(dataF_));
+
+    // perform ifft
+    fftw_execute( plan_b_ );
+
+    for(unsigned int s=0; s<nsteps_; ++s)
+      data_[s]*=1./(nsteps_); // the fft from fftw is not scaled !!!!!!!!
 
     // write sampled data
-    FBlas::copy(dim*nsteps_,fftR_,sampledData);
+    FBlas::c_copy(nsteps_,reinterpret_cast<FReal*>(data_),
+                  reinterpret_cast<FReal*>(sampledData));
 
   }
 
diff --git a/Tests/Kernels/testChebTensorialAlgorithm.cpp b/Tests/Kernels/testChebTensorialAlgorithm.cpp
index 2d522bbac4608331d40121344884815003b81462..876f365bc5e1bc0aeaf966d3deeaca5d7712b4f5 100644
--- a/Tests/Kernels/testChebTensorialAlgorithm.cpp
+++ b/Tests/Kernels/testChebTensorialAlgorithm.cpp
@@ -15,7 +15,7 @@
 // ===================================================================================
 
 // ==== CMAKE =====
-// @FUSE_FFT
+// @FUSE_BLAS
 // ================
 
 #include <iostream>
diff --git a/Tests/Utils/testFFTWMultidim.cpp b/Tests/Utils/testFFTWMultidim.cpp
index da668758ce1a2a4a498d6c08788f2c3a856a8129..44290c78be973e4d7f6d7fd9988085e4614fcd99 100644
--- a/Tests/Utils/testFFTWMultidim.cpp
+++ b/Tests/Utils/testFFTWMultidim.cpp
@@ -44,23 +44,24 @@ int main()
 
   //////////////////////////////////////////////////////////////////////////////
   // INITIALIZATION
-
   // size (pick a power of 2 for better performance of the FFT algorithm)
-  unsigned int rank = 2;
-  unsigned int nsteps_ = 500; 
-  unsigned int dim = 10;
-  const int steps[2]={static_cast<int>(dim),
-                      static_cast<int>(nsteps_)};
-  unsigned int size = dim*nsteps_;
-
+  const int dim=2;
+  const int pow_nsteps_=8;
+  int steps_[dim];
+  int nsteps_=1;
+  for(int d=0; d<dim; ++d) {
+    steps_[d]=FMath::pow(2,pow_nsteps_);
+    nsteps_*=steps_[d];
+  }
   //////////////////////////////////////////////////////////////////////////////
   // Multidimensionnal FFT PLANS
+  std::cout<< "Test "<< dim <<"D FFT."<<std::endl;
   // fftw arrays
   FReal* fftR_;
   FComplexe* fftC_;
 
-  fftR_ = (FReal*) fftw_malloc(sizeof(FReal) * size );
-  fftC_ = (FComplexe*) fftw_malloc(sizeof(FComplexe) * size );
+  fftR_ = (FReal*) fftw_malloc(sizeof(FReal) * nsteps_ );
+  fftC_ = (FComplexe*) fftw_malloc(sizeof(FComplexe) * nsteps_ );
 
   // fftw plans
   // use routine defined in file:
@@ -72,12 +73,12 @@ int main()
   time.tic();
 
   plan_c2r_ =
-    fftw_plan_dft_c2r(rank, steps, 
+    fftw_plan_dft_c2r(dim, steps_, 
                       reinterpret_cast<fftw_complex*>(fftC_),
                       fftR_, 
                       FFTW_MEASURE);
   plan_r2c_ =
-    fftw_plan_dft_r2c(rank, steps, 
+    fftw_plan_dft_r2c(dim, steps_, 
                       fftR_, 
                       reinterpret_cast<fftw_complex*>(fftC_), 
                       FFTW_MEASURE);
@@ -87,20 +88,20 @@ int main()
   //////////////////////////////////////////////////////////////////////////////
   // EXECUTION
   // generate random physical data
-  for(unsigned int s=0; s<size; ++s)
+  for(int s=0; s<nsteps_; ++s)
     fftR_[s] = FReal(rand())/FRandMax; 
 
 // // display data in  physical space
 // std::cout<< "Physical data: "<<std::endl;
-// for(unsigned int d=0; d<dim; ++d){
-//   for(unsigned int s=0; s<nsteps_; ++s)
-//     std::cout<< fftR_[s+d*nsteps_] << ", ";
+// for( int r=0; r<steps_[0]; ++r) {
+//   for( int s=0; s<steps_[1]; ++s)
+//     std::cout<< fftR_[r*steps_[1]+s] << ", ";
 //   std::cout<<std::endl;
 // }
 // std::cout<<std::endl;
 
   // perform fft
-  std::cout<< "Perform Forward FFT: ";
+ std::cout<< "Perform Forward FFT: ";
   time.tic();
   fftw_execute( plan_r2c_ );
   std::cout << "took " << time.tacAndElapsed() << "sec." << std::endl;
@@ -108,20 +109,23 @@ int main()
 //  // display transform in Fourier space
 //  // beware the real data FFT stores only N/2+1 complex output values
 //  std::cout<< "Transformed data : "<<std::endl;
-//  for(unsigned int d=0; d<dim; ++d){
-//    for(unsigned int s=0; s<nsteps_/2+1; ++s)
-//      std::cout<< fftC_[s+d*(nsteps_/2+1)] << ", ";
+//  for( int r=0; r<steps_[0]/2+1; ++r) {
+//    for( int s=0; s<steps_[1]/2+1; ++s)
+//      std::cout<< fftC_[r*(steps_[1]/2+1)+s] << ", ";
 //    std::cout<<std::endl;
 //  }
 //  std::cout<<std::endl;
 
-//  for(unsigned int s=0; s<nsteps_/2+1; ++s){
-//    fftC_[nsteps_-s]=FComplexe(fftC_[s].getReal(),-fftC_[s].getImag());
-//  }
-//
+////  for( int s=0; s<steps_[1]/2+1; ++s){
+////    fftC_[nsteps_-s]=FComplexe(fftC_[s].getReal(),-fftC_[s].getImag());
+////  }
+
 //  std::cout<< "Full Transformed data : "<<std::endl;
-//  for(unsigned int s=0; s<nsteps_; ++s)
-//    std::cout<< fftC_[s] << ", ";
+//  for( int r=0; r<steps_[0]; ++r){
+//    for( int s=0; s<steps_[1]; ++s)
+//      std::cout<< fftC_[r*steps_[1]+s] << ", ";
+//    std::cout<<std::endl;
+//  }
 //  std::cout<<std::endl;
 
   // perform ifft of tranformed data (in order to get physical values back)
@@ -132,9 +136,9 @@ int main()
 
 //  // display data in physical space
 //  std::cout<< "Physical data (from 1/N*IFFT(FFT(Physical data))): "<<std::endl;
-//  for(unsigned int d=0; d<dim; ++d){
-//    for(unsigned int s=0; s<nsteps_; ++s)
-//      std::cout<< fftR_[s+d*nsteps_]/(nsteps_*dim) << ", ";
+//  for( int r=0; r<steps_[0]; ++r) {
+//    for( int s=0; s<steps_[1]; ++s)
+//      std::cout<< fftR_[r*steps_[1]+s]/(nsteps_) << ", ";
 //    std::cout<<std::endl;
 //  }
 //  std::cout<<std::endl;
diff --git a/Tests/Utils/testFastDiscreteConvolution.cpp b/Tests/Utils/testFastDiscreteConvolution.cpp
index b9d109db85fa61ac84a171eb7c613ac0aed53e1b..107d53552d579c1f750a319f388ffb8c936656d5 100644
--- a/Tests/Utils/testFastDiscreteConvolution.cpp
+++ b/Tests/Utils/testFastDiscreteConvolution.cpp
@@ -106,8 +106,10 @@ int main(int, char **){
   // Init DFTor
   std::cout<< "Set DFT: ";
   time.tic();
+  const int dim = 1;
+  const int steps[dim] = {N};
   //FDft Dft(N);// direct version
-  FFft<1> Dft(N);// fast version
+  FFft<dim> Dft(steps);// fast version
   std::cout << "took " << time.tacAndElapsed() << "sec." << std::endl;
 
   // Initialize manually
@@ -135,7 +137,8 @@ int main(int, char **){
 
   // Transform first column of K
   FReal tK[N];
-  for(unsigned int i=0; i<N; ++i) tK[i]=K[i*N];
+  for(unsigned int i=0; i<N; ++i) tK[i]=K[i*N]; // first column
+//  for(unsigned int i=0; i<N; ++i) tK[i]=K[i]; // first row
   std::cout<< "Transform tK->FK: ";
   time.tic();
   Dft.applyDFT(tK,FK);
diff --git a/Tests/Utils/testLapack.cpp b/Tests/Utils/testLapack.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ecf3f327e5ebf0b835a23933767fb6d2eb8f1f63
--- /dev/null
+++ b/Tests/Utils/testLapack.cpp
@@ -0,0 +1,121 @@
+// ===================================================================================
+// Ce LOGICIEL "ScalFmm" est couvert par le copyright Inria 20xx-2012.
+// Inria détient tous les droits de propriété sur le LOGICIEL, et souhaite que
+// la communauté scientifique l'utilise afin de le tester et de l'évaluer.
+// Inria donne gracieusement le droit d'utiliser ce LOGICIEL. Toute utilisation
+// dans un but lucratif ou à des fins commerciales est interdite sauf autorisation
+// expresse et préalable d'Inria.
+// Toute utilisation hors des limites précisées ci-dessus et réalisée sans l'accord
+// expresse préalable d'Inria constituerait donc le délit de contrefaçon.
+// Le LOGICIEL étant un produit en cours de développement, Inria ne saurait assurer
+// aucune responsabilité et notamment en aucune manière et en aucun cas, être tenu
+// de répondre d'éventuels dommages directs ou indirects subits par l'utilisateur.
+// Tout utilisateur du LOGICIEL s'engage à communiquer à Inria ses remarques
+// relatives à l'usage du LOGICIEL
+// ===================================================================================
+
+// ==== CMAKE =====
+// @FUSE_BLAS
+// ================
+
+#include <iostream>
+#include <stdlib.h>
+
+#include "../../Src/Utils/FBlas.hpp"
+#include "../../Src/Utils/FTic.hpp"
+
+/**
+ * Test functionality of C - interfaced LAPACK functions
+ */
+
+int main()
+{
+  FTic time;
+
+  /*
+   * List of tested functions:
+   * Cholesky decomposition: FBlas::potrf()
+   * TODO SVD: FBlas::gesvd()
+   * TODO QR decomposition: FBlas::geqrf()
+   */
+
+	const unsigned int m = 4, n = 4;
+	FReal* A = new FReal [m * n]; // matrix: column major ordering
+
+  // A= LL^T ////////////////////////////////////
+  // define symmetric definite positive matrix A
+  A[0]=5; A[10]=4; A[15]=7;
+  A[1]=A[3]=A[4]=A[12]=2; 
+  A[6]=A[7]=A[9]=A[13]=1; 
+  A[2]=A[5]=A[8]=3;
+  A[11]=A[14]=-1; 
+
+  // copy A in C
+	FReal* C = new FReal [m * n]; // matrix: column major ordering
+	for (unsigned int ii=0; ii<m; ++ii) 
+		for (unsigned int jj=0; jj<n; ++jj)
+      C[ii*m + jj]=A[ii*m + jj];
+
+  std::cout<<"\nA=["<<std::endl;
+  for (unsigned int i=0; i<m; ++i) {
+    for (unsigned int j=0; j<n; ++j)
+      std::cout << A[i*n+j] << " ";
+    std::cout<< std::endl;
+  }
+  std::cout<<"]"<<std::endl;
+
+  // perform Cholesky decomposition
+  std::cout<<"\nCholesky decomposition ";
+  int INF = FBlas::potrf(m, A, n);
+  if(INF==0) {std::cout<<"succeeded!"<<std::endl;}
+  else {std::cout<<"failed!"<<std::endl;}
+
+  std::cout<<"\nA_out=["<<std::endl;
+  for (unsigned int i=0; i<m; ++i) {
+    for (unsigned int j=0; j<n; ++j)
+      std::cout << A[i*n+j] << " ";
+    std::cout<<std::endl;
+  }
+  std::cout<<"]"<<std::endl;
+
+  // build lower matrix
+  FReal* L = new FReal [m * n]; // matrix: column major ordering
+	for (unsigned int ii=0; ii<m; ++ii) 
+		for (unsigned int jj=0; jj<n; ++jj){ 
+      if(ii<=jj)
+        L[ii*m + jj]=A[ii*m + jj];
+      else
+        L[ii*m + jj]=0.;
+    }
+
+  std::cout<<"\nL=["<<std::endl;
+  for (unsigned int i=0; i<m; ++i) {
+    for (unsigned int j=0; j<n; ++j)
+      std::cout << L[i*n+j] << " ";
+    std::cout<< std::endl;
+  }
+  std::cout<<"]"<<std::endl;
+
+  // verify result by computing B=LL^T
+	FReal* B = new FReal [m * n]; // matrix: column major ordering
+	for (unsigned int ii=0; ii<m; ++ii)
+		for (unsigned int jj=0; jj<n; ++jj){ 
+      B[ii*m + jj]=0.;
+      for (unsigned int j=0; j<n; ++j)
+        B[ii*m + jj]+=L[j*m + ii]*L[j*m + jj];
+    }
+
+  std::cout<<"\nA-LL^T=["<<std::endl;
+  for (unsigned int i=0; i<m; ++i) {
+    for (unsigned int j=0; j<n; ++j)
+      std::cout << B[i*n+j]-C[i*n+j] << " ";
+    std::cout<< std::endl;
+  }
+  std::cout<<"]"<<std::endl;
+
+	delete [] A;
+	delete [] B;
+	delete [] C;
+
+	return 0;
+}
diff --git a/Tests/Utils/testUnifInterpolator.cpp b/Tests/Utils/testUnifInterpolator.cpp
index ac1a22838bf63b39f5592931f1d17f9d5bfeaff8..ffc0cb427e04a4ce5e7778e658b075d86b0d6451 100755
--- a/Tests/Utils/testUnifInterpolator.cpp
+++ b/Tests/Utils/testUnifInterpolator.cpp
@@ -211,14 +211,14 @@ int main(int, char **){
   // In order to actually embed K into a circulant matrix C one just
   // needs to insert (ORDER-1) extra lines/columns (to each block).
 
-//  std::cout<< "K=" <<std::endl;
-//  for (unsigned int i=0; i<nnodes; ++i){
-//    for (unsigned int j=0; j<nnodes; ++j){
-//      std::cout<< K[i*nnodes+j]<<", ";
-//    }
-//    std::cout<<std::endl;
-//  }
-//  std::cout<<std::endl;
+  std::cout<< "K=" <<std::endl;
+  for (unsigned int i=0; i<nnodes; ++i){
+    for (unsigned int j=0; j<nnodes; ++j){
+      std::cout<< K[i*nnodes+j]<<", ";
+    }
+    std::cout<<std::endl;
+  }
+  std::cout<<std::endl;
 
   // Check multi-index
   std::cout<< "node_ids=" <<std::endl;
@@ -324,9 +324,10 @@ int main(int, char **){
   FBlas::setzero(rc,PLocalExp);
 
   // Init DFT
-  // NB: only one FFTor is defined a since scalar problems involve scalar r/lhs and scalar matrix kernel. All dimensions are 1. 
+  const int dimfft = 1;
+  const int steps[dimfft] = {rc};
   //FDft Dft(rc); // direct version
-  FFft<1> Dft(rc); // fast version
+  FFft<dimfft> Dft(steps); // fast version
 
   // Get first COLUMN of K and Store in T
   FReal T[rc];
diff --git a/Tests/Utils/testUnifTensorialInterpolator.cpp b/Tests/Utils/testUnifTensorialInterpolator.cpp
index 1602b0052ae7423f18d272aa2fbcb06108f40067..5429e743b1ecae8b8d6395eb50c2e07ce38dc7ab 100755
--- a/Tests/Utils/testUnifTensorialInterpolator.cpp
+++ b/Tests/Utils/testUnifTensorialInterpolator.cpp
@@ -356,10 +356,10 @@ int main(int, char **){
   // Efficient application of the Toeplitz system in FOURIER SPACE
 
   // Init DFT
+  const int dimfft = 1;
+  const int steps[dimfft] = {rc};
   //FDft Dft(rc); // direct version
-  FFft<1/*(TODO: fix MultidimFFT) nK*/> DftK(rc); // fast version
-  FFft<1/*nrhs*/> DftRhs(rc);
-  FFft<1/*nlhs*/> DftLhs(rc);
+  FFft<dimfft> Dft(steps); // fast version
 
   // Get first COLUMN of K and Store in T
   FReal T[dim*rc];
@@ -393,9 +393,9 @@ int main(int, char **){
 
   // if first COLUMN (T) of C is used
   for (unsigned int d=0; d<dim; ++d)
-    DftK.applyDFT(T+d*rc,FT+d*rc);
+    Dft.applyDFT(T+d*rc,FT+d*rc);
 //  // if first ROW of C is used
-//  DftK.applyDFT(C,FT);
+//  Dft.applyDFT(C,FT);
 
   FComplexe FPMultExp[nrhs*rc];
   FComplexe FPLocalExp[nlhs*rc];
@@ -428,7 +428,7 @@ int main(int, char **){
 
   // Transform PaddedMultExp
   for (unsigned int idxRhs=0; idxRhs<nrhs; ++idxRhs) // apply nrhs 1 dimensionnal FFT
-    DftRhs.applyDFT(PaddedMultExp+idxRhs*rc,FPMultExp+idxRhs*rc);
+    Dft.applyDFT(PaddedMultExp+idxRhs*rc,FPMultExp+idxRhs*rc);
 
   std::cout<< "Apply M2L in Fourier space: ";
   time.tic();
@@ -463,7 +463,7 @@ int main(int, char **){
 //    std::cout<<std::endl;
 
   for (unsigned int idxLhs=0; idxLhs<nlhs; ++idxLhs) // apply nrhs 1 dimensionnal FFT
-    DftLhs.applyIDFT(FPLocalExp+idxLhs*rc,PLocalExp+idxLhs*rc);
+    Dft.applyIDFT(FPLocalExp+idxLhs*rc,PLocalExp+idxLhs*rc);
 
 //  std::cout<< "Padded LocalExp: "<<std::endl;
 //  for (unsigned int p=0; p<rc; ++p)