diff --git a/Src/Kernels/Uniform/FUnifCell.hpp b/Src/Kernels/Uniform/FUnifCell.hpp
index 1f538404050068623d7b406337f1b109e8325707..a56f4ceaa34f154fa0f226035f7eef5bcaa73024 100644
--- a/Src/Kernels/Uniform/FUnifCell.hpp
+++ b/Src/Kernels/Uniform/FUnifCell.hpp
@@ -4,12 +4,12 @@
 #ifndef FUNIFCELL_HPP
 #define FUNIFCELL_HPP
 
+#include "Utils/stdComplex.hpp"
 
 #include "./FUnifTensor.hpp"
 #include "../../Components/FBasicCell.hpp"
 #include "../../Extensions/FExtendCellType.hpp"
 
-#include "../../Utils/FComplex.hpp"
 
 /**
  * @author Pierre Blanchard (pierre.blanchard@inria.fr)
@@ -35,15 +35,15 @@ public:
     struct exp_impl {
         FReal exp[N * NVALS * VectorSize]; //< Multipole expansion
         // Multipole expansion in Fourier space
-        FComplex<FReal> transformed_exp[N * NVALS * TransformedVectorSize];
+        stdComplex<FReal> transformed_exp[N * NVALS * TransformedVectorSize];
 
         const FReal* get(const int inRhs) const
         { return this->exp + inRhs*VectorSize; }
         FReal* get(const int inRhs)
         { return this->exp + inRhs*VectorSize; }
-        const FComplex<FReal>* getTransformed(const int inRhs) const
+        const stdComplex<FReal>* getTransformed(const int inRhs) const
         { return this->transformed_exp + inRhs*TransformedVectorSize; }
-        FComplex<FReal>* getTransformed(const int inRhs)
+        stdComplex<FReal>* getTransformed(const int inRhs)
         { return this->transformed_exp + inRhs*TransformedVectorSize; }
 
         constexpr int getVectorSize() const {
@@ -80,7 +80,7 @@ public:
 
         FSize getSavedSize() const {
             return N * NVALS * VectorSize * (FSize) sizeof(FReal)
-                + N * NVALS * TransformedVectorSize * (FSize) sizeof(FComplex<FReal>);
+                + N * NVALS * TransformedVectorSize * (FSize) sizeof(stdComplex<FReal>);
         }
     };
 
diff --git a/Src/Kernels/Uniform/FUnifKernel.hpp b/Src/Kernels/Uniform/FUnifKernel.hpp
index 42b39a6742ed8e701fa2d729ad47dd5388b2346c..58adc254e43ffdef1b06331b0b4f7458510328a5 100644
--- a/Src/Kernels/Uniform/FUnifKernel.hpp
+++ b/Src/Kernels/Uniform/FUnifKernel.hpp
@@ -131,7 +131,7 @@ public:
         const FReal scale(MatrixKernel->getScaleFactor(CellWidth));
 
         for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
-            FComplex<FReal> *const TransformedLocalExpansion = TargetExpansion->getTransformed(idxRhs);
+            stdComplex<FReal> *const TransformedLocalExpansion = TargetExpansion->getTransformed(idxRhs);
 
             for(int idxExistingNeigh = 0 ; idxExistingNeigh < inSize ; ++idxExistingNeigh){
                 const int idxNeigh = neighborPositions[idxExistingNeigh];
diff --git a/Src/Kernels/Uniform/FUnifM2LHandler.hpp b/Src/Kernels/Uniform/FUnifM2LHandler.hpp
index 23a8956cb160c1124517c68c3e2a1f11b2e62a70..eebbdae48ed76d1446c2c5ac9bdea5f43bfaac6c 100644
--- a/Src/Kernels/Uniform/FUnifM2LHandler.hpp
+++ b/Src/Kernels/Uniform/FUnifM2LHandler.hpp
@@ -16,7 +16,7 @@
 #include "Utils/FDft.hpp"
 #include "Utils/FSmartPointer.hpp"
 
-#include "Utils/FComplex.hpp"
+#include "Utils/stdComplex.hpp"
 
 #include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
 #include "FUnifTensor.hpp"
@@ -24,7 +24,7 @@
 /*!  Precomputation of the 316 interactions by evaluation of the matrix kernel on the uniform grid and transformation into Fourier space.
 PB: Compute() does not belong to the M2LHandler like it does in the Chebyshev kernel. This allows much nicer specialization of the M2LHandler class with respect to the homogeneity of the kernel of interaction like in the ChebyshevSym kernel.*/
 template < class FReal,int ORDER, typename MatrixKernelClass>
-static void Compute(const MatrixKernelClass *const MatrixKernel, const FReal CellWidth, FComplex<FReal>* &FC, const int SeparationCriterion = 1)
+static void Compute(const MatrixKernelClass *const MatrixKernel, const FReal CellWidth, stdComplex<FReal>* &FC, const int SeparationCriterion = 1)
 {
     // allocate memory and store compressed M2L operators
     if (FC) throw std::runtime_error("M2L operators are already set");
@@ -41,19 +41,19 @@ static void Compute(const MatrixKernelClass *const MatrixKernel, const FReal Cel
 
     // allocate memory and compute 316 m2l operators (342 if separation equals 0, 343 if separation equals -1)
     FReal           *_C;
-    FComplex<FReal> *_FC;
+    stdComplex<FReal> *_FC;
 
     // reduce storage from nnodes^2=order^6 to (2order-1)^3
     const unsigned int rc = (2*order-1)*(2*order-1)*(2*order-1);
     _C = new FReal [rc];
-    _FC = new FComplex<FReal> [rc * ninteractions];
+    _FC = new stdComplex<FReal> [rc * ninteractions];
 
     // initialize root node ids pairs
     unsigned int node_ids_pairs[rc][2];
     TensorType::setNodeIdsPairs(node_ids_pairs);
     // init Discrete Fourier Transformator
     const int dimfft = 1; // unidim FFT since fully circulant embedding
-    FFftw<FReal,FComplex<FReal>,dimfft> Dft(rc);
+    FFftw<FReal,stdComplex<FReal>,dimfft> Dft(rc);
     // get first column of K via permutation
     unsigned int perm[rc];
     TensorType::setStoragePermutation(perm);
@@ -104,7 +104,7 @@ static void Compute(const MatrixKernelClass *const MatrixKernel, const FReal Cel
     // reduce storage if real valued kernel
     const unsigned int opt_rc = rc/2+1;
     // allocate M2L
-    FC = new FComplex<FReal>[343 * opt_rc];
+    FC = new stdComplex<FReal>[343 * opt_rc];
 
     for (int i=-3; i<=3; ++i)
         for (int j=-3; j<=3; ++j)
@@ -159,10 +159,10 @@ public:
 
     // DFT specific
     static const int dimfft = 1; // unidim FFT since fully circulant embedding
-    typedef FFftw<FReal,FComplex<FReal>,dimfft> DftClass; // Fast Discrete Fourier Transformator
+    typedef FFftw<FReal,stdComplex<FReal>,dimfft> DftClass; // Fast Discrete Fourier Transformator
 private:
     /// M2L Operators (stored in Fourier space)
-    FSmartPointer< FComplex<FReal>,FSmartArrayMemory> FC;
+    FSmartPointer< stdComplex<FReal>,FSmartArrayMemory> FC;
 
     /// Utils
     unsigned int node_diff[nnodes*nnodes];
@@ -223,12 +223,12 @@ public:
         if (FC) throw std::runtime_error("M2L operator already set");
         // Compute matrix of interactions
         const FReal ReferenceCellWidth = FReal(2.);
-        FComplex<FReal>* pFC = NULL;
+        stdComplex<FReal>* pFC = NULL;
         Compute<FReal,order>(MatrixKernel,ReferenceCellWidth,pFC,LeafLevelSeparationCriterion);
         FC.assign(pFC);
 
         // Compute memory usage
-        unsigned long sizeM2L = 343*opt_rc*sizeof(FComplex<FReal>);
+        unsigned long sizeM2L = 343*opt_rc*sizeof(stdComplex<FReal>);
 
 
         // write info
@@ -237,7 +237,7 @@ public:
     }
 
     unsigned long long getMemory() const {
-        return 343*opt_rc*sizeof(FComplex<FReal>);
+        return 343*opt_rc*sizeof(stdComplex<FReal>);
     }
 
     /**
@@ -247,7 +247,7 @@ public:
     * @param[in] X transformed local expansion of size \f$r\f$
     * @param[out] x local expansion of size \f$\ell^3\f$
     */
-    void unapplyZeroPaddingAndDFT(const FComplex<FReal> *const FX, FReal *const x) const
+    void unapplyZeroPaddingAndDFT(const stdComplex<FReal> *const FX, FReal *const x) const
     {
         FReal Px[rc];
         FBlas::setzero(rc,Px);
@@ -273,13 +273,16 @@ public:
      * @param[in] CellWidth needed for the scaling of the compressed M2L operators which are based on a homogeneous matrix kernel computed for the reference cell width \f$w=2\f$, ie in \f$[-1,1]^3\f$.
      */
     void applyFC(const unsigned int idx, const unsigned int, const FReal scale,
-                 const FComplex<FReal> *const FY, FComplex<FReal> *const FX) const
+                 const stdComplex<FReal> *const FY, stdComplex<FReal> *const FX) const
     {
         // Perform entrywise product manually
         for (unsigned int j=0; j<opt_rc; ++j){
-            FX[j].addMul(FComplex<FReal>(scale*FC[idx*opt_rc + j].getReal(),
-                                  scale*FC[idx*opt_rc + j].getImag()),
-                         FY[j]);
+//            FX[j].addMul(newFComplex<FReal>(scale*FC[idx*opt_rc + j].real(),
+//                                  scale*FC[idx*opt_rc + j].imag()),
+//                         FY[j]);
+            FX[j] += stdComplex<FReal>(scale*FC[idx*opt_rc + j].real(),
+                                  scale*FC[idx*opt_rc + j].imag())
+                         *FY[j];
         }
     }
 
@@ -291,7 +294,7 @@ public:
      * @param[in] y multipole expansion of size \f$\ell^3\f$
      * @param[out] Y transformed multipole expansion of size \f$r\f$
      */
-    void applyZeroPaddingAndDFT(FReal *const y, FComplex<FReal> *const FY) const
+    void applyZeroPaddingAndDFT(FReal *const y, stdComplex<FReal> *const FY) const
     {
         FReal Py[rc];
         FBlas::setzero(rc,Py);
@@ -303,7 +306,7 @@ public:
     }
 
 
-    const FComplex<FReal>& getFc(const int i, const int j) const{
+    const stdComplex<FReal>& getFc(const int i, const int j) const{
         return FC[i*opt_rc + j];
     }
 };
@@ -319,7 +322,7 @@ class FUnifM2LHandler<FReal,ORDER,NON_HOMOGENEOUS>
           rc = (2*ORDER-1)*(2*ORDER-1)*(2*ORDER-1)};
 
     /// M2L Operators (stored in Fourier space for each level)
-    FSmartPointer< FComplex<FReal>*,FSmartArrayMemory> FC;
+    FSmartPointer< stdComplex<FReal>*,FSmartArrayMemory> FC;
     /// Homogeneity specific variables
     const unsigned int TreeHeight;
     const FReal RootCellWidth;
@@ -328,7 +331,7 @@ class FUnifM2LHandler<FReal,ORDER,NON_HOMOGENEOUS>
     unsigned int node_diff[nnodes*nnodes];
     /// DFT specific
     static const int dimfft = 1; // unidim FFT since fully circulant embedding
-    typedef FFftw<FReal,FComplex<FReal>,dimfft> DftClass; // Fast real-valued Discrete Fourier Transformator
+    typedef FFftw<FReal,stdComplex<FReal>,dimfft> DftClass; // Fast real-valued Discrete Fourier Transformator
     DftClass Dft;
     const unsigned int opt_rc; // specific to real valued kernel
 
@@ -357,7 +360,7 @@ public:
         TensorType::setNodeIdsDiff(node_diff);
 
         // init M2L operators
-        FC = new FComplex<FReal>*[TreeHeight];
+        FC = new stdComplex<FReal>*[TreeHeight];
         FC[0]       = NULL; FC[1]       = NULL;
         for (unsigned int l=2; l<TreeHeight; ++l) FC[l] = NULL;
 
@@ -412,7 +415,7 @@ public:
         }
 
         // Compute memory usage
-        unsigned long sizeM2L = (TreeHeight-2)*343*opt_rc*sizeof(FComplex<FReal>);
+        unsigned long sizeM2L = (TreeHeight-2)*343*opt_rc*sizeof(stdComplex<FReal>);
 
         // write info
         std::cout << "Compute and set M2L operators ("<< long(sizeM2L/**1e-6*/) <<" B) in "
@@ -420,7 +423,7 @@ public:
     }
 
     unsigned long long getMemory() const {
-        return (TreeHeight-2)*343*opt_rc*sizeof(FComplex<FReal>);
+        return (TreeHeight-2)*343*opt_rc*sizeof(stdComplex<FReal>);
     }
 
     /**
@@ -430,7 +433,7 @@ public:
     * @param[in] X transformed local expansion of size \f$r\f$
     * @param[out] x local expansion of size \f$\ell^3\f$
     */
-    void unapplyZeroPaddingAndDFT(const FComplex<FReal> *const FX, FReal *const x) const
+    void unapplyZeroPaddingAndDFT(const stdComplex<FReal> *const FX, FReal *const x) const
     {
         FReal Px[rc];
         FBlas::setzero(rc,Px);
@@ -456,11 +459,12 @@ public:
      * @param[in] CellWidth needed for the scaling of the compressed M2L operators which are based on a homogeneous matrix kernel computed for the reference cell width \f$w=2\f$, ie in \f$[-1,1]^3\f$.
      */
     void applyFC(const unsigned int idx, const unsigned int TreeLevel, const FReal,
-                 const FComplex<FReal> *const FY, FComplex<FReal> *const FX) const
+                 const stdComplex<FReal> *const FY, stdComplex<FReal> *const FX) const
     {
         // Perform entrywise product manually
         for (unsigned int j=0; j<opt_rc; ++j){
-            FX[j].addMul(FC[TreeLevel][idx*opt_rc + j],FY[j]);
+  //          FX[j].addMul(FC[TreeLevel][idx*opt_rc + j],FY[j]);
+            FX[j] += (FC[TreeLevel][idx*opt_rc + j]*FY[j]);
         }
     }
 
@@ -472,7 +476,7 @@ public:
      * @param[in] y multipole expansion of size \f$\ell^3\f$
      * @param[out] Y transformed multipole expansion of size \f$r\f$
      */
-    void applyZeroPaddingAndDFT(FReal *const y, FComplex<FReal> *const FY) const
+    void applyZeroPaddingAndDFT(FReal *const y, stdComplex<FReal> *const FY) const
     {
         FReal Py[rc];
         FBlas::setzero(rc,Py);
@@ -483,7 +487,7 @@ public:
         Dft.applyDFT(Py,FY);
     }
 
-    const FComplex<FReal>& getFc(const int i, const int j) const{
+    const stdComplex<FReal>& getFc(const int i, const int j) const{
         return FC[i*opt_rc + j];
     }
 
diff --git a/Src/Kernels/Uniform/FUnifTensorialKernel.hpp b/Src/Kernels/Uniform/FUnifTensorialKernel.hpp
index 0e2d070335d8e753106b4ea5321e199571e044f9..8d05c6d290c29b70842d28bc7f988b951ec2e8e5 100644
--- a/Src/Kernels/Uniform/FUnifTensorialKernel.hpp
+++ b/Src/Kernels/Uniform/FUnifTensorialKernel.hpp
@@ -172,7 +172,7 @@ public:
                 const int idxLoc = idxV*nLhs + idxLhs;
 
                 // load transformed local expansion
-                FComplex<FReal> *const TransformedLocalExpansion = TargetExpansion->getTransformed(idxLoc);
+                stdComplex<FReal> *const TransformedLocalExpansion = TargetExpansion->getTransformed(idxLoc);
 
                 // update idxRhs
                 const int idxRhs = idxLhs % nPV;
diff --git a/Src/Kernels/Uniform/FUnifTensorialM2LHandler.hpp b/Src/Kernels/Uniform/FUnifTensorialM2LHandler.hpp
index 90041e5b12450c7b65d02c82a41c4f587733ee9e..64fbd11ebd54bb0654f24a6c3873c90d224cff1c 100644
--- a/Src/Kernels/Uniform/FUnifTensorialM2LHandler.hpp
+++ b/Src/Kernels/Uniform/FUnifTensorialM2LHandler.hpp
@@ -15,7 +15,7 @@
 #include "Utils/FTic.hpp"
 #include "Utils/FDft.hpp"
 
-#include "Utils/FComplex.hpp"
+#include "Utils/stdComplex.hpp"
 
 
 #include "FUnifTensor.hpp"
@@ -27,7 +27,7 @@ template < class FReal, int ORDER, class MatrixKernelClass>
 static void Compute(const MatrixKernelClass *const MatrixKernel, 
                     const FReal CellWidth, 
                     const FReal CellWidthExtension, 
-                    FComplex<FReal>** &FC, 
+                    stdComplex<FReal>** &FC,
                     const int SeparationCriterion = 1)
 {
     // dimensions of operators
@@ -51,16 +51,16 @@ static void Compute(const MatrixKernelClass *const MatrixKernel,
 
     // allocate memory and compute 316 m2l operators
     FReal** _C;
-    FComplex<FReal>** _FC;
+    stdComplex<FReal>** _FC;
 
     // reduce storage from nnodes^2=order^6 to (2order-1)^3
     const unsigned int rc = (2*order-1)*(2*order-1)*(2*order-1);
     _C  = new FReal* [ncmp];
-    _FC = new FComplex<FReal>* [ncmp];
+    _FC = new stdComplex<FReal>* [ncmp];
     for (unsigned int d=0; d<ncmp; ++d) 
     _C[d]  = new FReal[rc];
     for (unsigned int d=0; d<ncmp; ++d) 
-    _FC[d] = new FComplex<FReal>[rc * ninteractions];
+    _FC[d] = new stdComplex<FReal>[rc * ninteractions];
 
     // initialize root node ids pairs
     unsigned int node_ids_pairs[rc][2];
@@ -68,7 +68,7 @@ static void Compute(const MatrixKernelClass *const MatrixKernel,
 
     // init Discrete Fourier Transformator
     const int dimfft = 1; // unidim FFT since fully circulant embedding
-    FFftw<FReal,FComplex<FReal>,dimfft> Dft(rc);
+    FFftw<FReal,stdComplex<FReal>,dimfft> Dft(rc);
 
     // get first column of K via permutation
     unsigned int perm[rc];
@@ -110,7 +110,7 @@ static void Compute(const MatrixKernelClass *const MatrixKernel,
                         Dft.applyDFT(_C[d],_FC[d]+counter*rc);
 
                     // increment interaction counter
-                    counter++;
+                    ++counter;
                 }
             }
         }
@@ -131,7 +131,7 @@ static void Compute(const MatrixKernelClass *const MatrixKernel,
     const unsigned int opt_rc = rc/2+1;
     // allocate M2L
     for (unsigned int d=0; d<ncmp; ++d) 
-    FC[d] = new FComplex<FReal>[343 * opt_rc];
+    FC[d] = new stdComplex<FReal>[343 * opt_rc];
 
     for (int i=-3; i<=3; ++i)
         for (int j=-3; j<=3; ++j)
@@ -141,7 +141,7 @@ static void Compute(const MatrixKernelClass *const MatrixKernel,
                     for (unsigned int d=0; d<ncmp; ++d) 
                         FBlas::c_copy(opt_rc, reinterpret_cast<FReal*>(_FC[d] + counter*rc), 
                                       reinterpret_cast<FReal*>(FC[d] + idx*opt_rc));
-                    counter++;
+                    ++counter;
                 } else { 
                     for (unsigned int d=0; d<ncmp; ++d) 
                         FBlas::c_setzero(opt_rc, reinterpret_cast<FReal*>(FC[d] + idx*opt_rc));
@@ -181,7 +181,7 @@ class FUnifTensorialM2LHandler<FReal, ORDER,MatrixKernelClass,HOMOGENEOUS>
           ncmp = MatrixKernelClass::NCMP};
 
     /// M2L Operators (stored in Fourier space for each component of tensor)
-    FSmartPointer< FComplex<FReal>*,FSmartArrayMemory> FC;
+    FSmartPointer< stdComplex<FReal>*,FSmartArrayMemory> FC;
 
     const FReal CellWidthExtension; //<! extension of cells width
 
@@ -191,7 +191,7 @@ class FUnifTensorialM2LHandler<FReal, ORDER,MatrixKernelClass,HOMOGENEOUS>
 
     /// DFT specific
     static const int dimfft = 1; // unidim FFT since fully circulant embedding
-    typedef FFftw<FReal,FComplex<FReal>,dimfft> DftClass; // Fast Discrete Fourier Transformator
+    typedef FFftw<FReal,stdComplex<FReal>,dimfft> DftClass; // Fast Discrete Fourier Transformator
     DftClass Dft;
     const unsigned int opt_rc; // specific to real valued kernel
 
@@ -215,7 +215,7 @@ public:
     {
 
         // allocate FC
-        FC = new FComplex<FReal>*[ncmp];
+        FC = new stdComplex<FReal>*[ncmp];
         for (unsigned int d=0; d<ncmp; ++d)
             FC[d]=nullptr;
 
@@ -264,7 +264,7 @@ public:
         Compute<FReal,order>(MatrixKernel,ReferenceCellWidth, 0., FC, LeafLevelSeparationCriterion);
         
         // Compute memory usage
-        unsigned long sizeM2L = 343*ncmp*opt_rc*sizeof(FComplex<FReal>);
+        unsigned long sizeM2L = 343*ncmp*opt_rc*sizeof(stdComplex<FReal>);
 
         // write info
         std::cout << "Compute and set M2L operators ("<< long(sizeM2L/**1e-6*/) <<" B) in "
@@ -272,7 +272,7 @@ public:
     }
 
     unsigned long long getMemory() const {
-        return 343*ncmp*opt_rc*sizeof(FComplex<FReal>);
+        return 343*ncmp*opt_rc*sizeof(stdComplex<FReal>);
     }        
 
     /**
@@ -282,7 +282,7 @@ public:
      * @param[in] X transformed local expansion of size \f$r\f$
      * @param[out] x local expansion of size \f$\ell^3\f$
      */
-    void unapplyZeroPaddingAndDFT(const FComplex<FReal> *const FX, FReal *const x) const
+    void unapplyZeroPaddingAndDFT(const stdComplex<FReal> *const FX, FReal *const x) const
     { 
         FReal Px[rc];
         FBlas::setzero(rc,Px);
@@ -309,11 +309,11 @@ public:
      */
     void applyFC(const unsigned int idx, const unsigned int , 
                  const FReal scale, const unsigned int d,
-                 const FComplex<FReal> *const FY, FComplex<FReal> *const FX) const
+                 const stdComplex<FReal> *const FY, stdComplex<FReal> *const FX) const
     {
         // Perform entrywise product manually
         for (unsigned int j=0; j<opt_rc; ++j){
-            FX[j].addMul(FComplex<FReal>(scale*FC[d][idx*opt_rc + j].getReal(),
+            FX[j].addMul(stdComplex<FReal>(scale*FC[d][idx*opt_rc + j].getReal(),
                                   scale*FC[d][idx*opt_rc + j].getImag()), 
                          FY[j]);
         }
@@ -328,7 +328,7 @@ public:
      * @param[in] y multipole expansion of size \f$\ell^3\f$
      * @param[out] Y transformed multipole expansion of size \f$r\f$
      */
-    void applyZeroPaddingAndDFT(FReal *const y, FComplex<FReal> *const FY) const
+    void applyZeroPaddingAndDFT(FReal *const y, stdComplex<FReal> *const FY) const
     {
         FReal Py[rc];
         FBlas::setzero(rc,Py);
@@ -355,7 +355,7 @@ class FUnifTensorialM2LHandler<FReal,ORDER,MatrixKernelClass,NON_HOMOGENEOUS>
 
 
     /// M2L Operators (stored in Fourier space for each component and each level)
-    FSmartPointer< FComplex<FReal>**,FSmartArrayMemory> FC;
+    FSmartPointer< stdComplex<FReal>**,FSmartArrayMemory> FC;
 
     /// Homogeneity specific variables
     const unsigned int TreeHeight;
@@ -368,7 +368,7 @@ class FUnifTensorialM2LHandler<FReal,ORDER,MatrixKernelClass,NON_HOMOGENEOUS>
 
     /// DFT specific
     static const int dimfft = 1; // unidim FFT since fully circulant embedding
-    typedef FFftw<FReal,FComplex<FReal>,dimfft> DftClass; // Fast Discrete Fourier Transformator
+    typedef FFftw<FReal,stdComplex<FReal>,dimfft> DftClass; // Fast Discrete Fourier Transformator
     DftClass Dft;
     const unsigned int opt_rc; // specific to real valued kernel
 
@@ -394,9 +394,9 @@ public:
     {
 
         // allocate FC
-        FC = new FComplex<FReal>**[TreeHeight];
+        FC = new stdComplex<FReal>**[TreeHeight];
         for (unsigned int l=0; l<TreeHeight; ++l){ 
-            FC[l] = new FComplex<FReal>*[ncmp];
+            FC[l] = new stdComplex<FReal>*[ncmp];
             for (unsigned int d=0; d<ncmp; ++d)
                 FC[l][d]=nullptr;
         }
@@ -453,7 +453,7 @@ public:
         }
 
         // Compute memory usage
-        unsigned long sizeM2L = (TreeHeight-2)*343*ncmp*opt_rc*sizeof(FComplex<FReal>);
+        unsigned long sizeM2L = (TreeHeight-2)*343*ncmp*opt_rc*sizeof(stdComplex<FReal>);
 
         // write info
         std::cout << "Compute and set M2L operators ("<< long(sizeM2L/**1e-6*/) <<" B) in "
@@ -461,7 +461,7 @@ public:
     }
 
     unsigned long long getMemory() const {
-        return (TreeHeight-2)*343*ncmp*opt_rc*sizeof(FComplex<FReal>);
+        return (TreeHeight-2)*343*ncmp*opt_rc*sizeof(stdComplex<FReal>);
     }        
 
     /**
@@ -471,7 +471,7 @@ public:
      * @param[in] X transformed local expansion of size \f$r\f$
      * @param[out] x local expansion of size \f$\ell^3\f$
      */
-    void unapplyZeroPaddingAndDFT(const FComplex<FReal> *const FX, FReal *const x) const
+    void unapplyZeroPaddingAndDFT(const stdComplex<FReal> *const FX, FReal *const x) const
     { 
         FReal Px[rc];
         FBlas::setzero(rc,Px);
@@ -498,11 +498,12 @@ public:
      */
     void applyFC(const unsigned int idx, const unsigned int TreeLevel, 
                  const FReal, const unsigned int d,
-                 const FComplex<FReal> *const FY, FComplex<FReal> *const FX) const
+                 const stdComplex<FReal> *const FY, stdComplex<FReal> *const FX) const
     {
         // Perform entrywise product manually
         for (unsigned int j=0; j<opt_rc; ++j){
-            FX[j].addMul(FC[TreeLevel][d][idx*opt_rc + j], FY[j]);
+//           FX[j].addMul(FC[TreeLevel][d][idx*opt_rc + j], FY[j]);
+            FX[j] += (FC[TreeLevel][d][idx*opt_rc + j]*FY[j]);
         }
 
     }
@@ -515,7 +516,7 @@ public:
      * @param[in] y multipole expansion of size \f$\ell^3\f$
      * @param[out] Y transformed multipole expansion of size \f$r\f$
      */
-    void applyZeroPaddingAndDFT(FReal *const y, FComplex<FReal> *const FY) const
+    void applyZeroPaddingAndDFT(FReal *const y, stdComplex<FReal> *const FY) const
     {
         FReal Py[rc];
         FBlas::setzero(rc,Py);
diff --git a/Src/Utils/FComplex.hpp b/Src/Utils/FComplex.hpp
index 587a783feb31906f9dd66a0c85669e0ebda54549..06b5bc6b9f1eae043b5ee3590471790dd319f856 100644
--- a/Src/Utils/FComplex.hpp
+++ b/Src/Utils/FComplex.hpp
@@ -74,11 +74,19 @@ public:
         return !(*this == other);
     }
 
+    /** Get imaginary */
+    FReal imag() const{
+        return this->complex[1];
+    }
     /** Get imaginary */
     FReal getImag() const{
         return this->complex[1];
     }
 
+    /** Get complex[0] */
+    FReal real() const{
+        return this->complex[0];
+    }
     /** Get complex[0] */
     FReal getReal() const{
         return this->complex[0];
diff --git a/Src/Utils/FDft.hpp b/Src/Utils/FDft.hpp
index 4854e38951cdf6a40c8b4e808470d7299a0c8cb1..967630bc6ebd8d49d3d1998461007e2e9a70cbd6 100644
--- a/Src/Utils/FDft.hpp
+++ b/Src/Utils/FDft.hpp
@@ -10,7 +10,10 @@
 
 #ifdef SCALFMM_USE_FFT
 #include <iostream>
-#include <stdlib.h>
+#include <cstdlib>
+#include <complex>
+#include "Utils/stdComplex.hpp"
+
 
 #ifdef SCALFMM_USE_ESSL_AS_FFTW
 #include <fftw3_essl.h>
@@ -21,7 +24,7 @@
 // for @class FDft only
 #include "FBlas.hpp"
 
-#include "FComplex.hpp"
+
 
 
 /**
@@ -35,7 +38,7 @@
  * 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 @class FDft is templatized with the input value type (FReal or FComplex<FReal>),
+ * The @class FDft is templatized with the input value type (FReal or stdComplex<FReal>),
  * while @class FFftw is templatized with input and output value types and the dimension.
  *
  * The aim of writing these specific classes is to allow further customization
@@ -50,12 +53,12 @@
  *
  * @tparam size_ number of sampled values \f$N\f$
  */
-template< class FReal, typename ValueType = FReal>
+template< class FReal, typename ValueType = FReal >
 class FDft
 {
 
-    ValueType* data_;  //< data in physical space
-    FComplex<FReal>* dataF_; //< data in Fourier space
+    ValueType* data_;        //< data in physical space
+    stdComplex<FReal>* dataF_; //< data in Fourier space
 
     FReal *cosRS_, *sinRS_;
 
@@ -66,7 +69,7 @@ private:
     {
         // allocate arrays
         data_  = new ValueType[size_];
-        dataF_ = new FComplex<FReal>[size_];
+        dataF_ = new stdComplex<FReal>[size_];
 
         // Beware this is extremely HEAVY to store!!! => Use FDft only for debug!
         cosRS_ = new FReal[size_*size_];
@@ -107,7 +110,7 @@ public:
     /// Forward DFT
     // Real valued DFT
     void applyDFT(const FReal* sampledData,
-                  FComplex<FReal>* transformedData) const
+                  stdComplex<FReal>* transformedData) const
     {
         // read sampled data
         FBlas::c_setzero(size_,reinterpret_cast<FReal*>(dataF_));
@@ -116,7 +119,7 @@ public:
         // perform direct forward transformation
         for(unsigned int r=0; r<size_; ++r)
             for(unsigned int s=0; s<size_; ++s){
-                dataF_[r] += FComplex<FReal>(data_[s]*cosRS_[r*size_+s],
+                dataF_[r] += stdComplex<FReal>(data_[s]*cosRS_[r*size_+s],
                         -data_[s]*sinRS_[r*size_+s]);
             }
 
@@ -125,8 +128,8 @@ public:
                       reinterpret_cast<FReal*>(transformedData));
     }
     // Complexe valued DFT
-    void applyDFT(const FComplex<FReal>* sampledData,
-                  FComplex<FReal>* transformedData) const
+    void applyDFT(const stdComplex<FReal>* sampledData,
+                  stdComplex<FReal>* transformedData) const
     {
         // read sampled data
         FBlas::c_setzero(size_,reinterpret_cast<FReal*>(dataF_));
@@ -136,7 +139,7 @@ public:
         // perform direct forward transformation
         for(unsigned int r=0; r<size_; ++r)
             for(unsigned int s=0; s<size_; ++s){
-                dataF_[r] += FComplex<FReal>(data_[s].getReal()*cosRS_[r*size_+s]
+                dataF_[r] += stdComplex<FReal>(data_[s].getReal()*cosRS_[r*size_+s]
                         + data_[s].getImag()*sinRS_[r*size_+s],
                         data_[s].getImag()*cosRS_[r*size_+s]
                         - data_[s].getReal()*sinRS_[r*size_+s]);
@@ -149,7 +152,7 @@ public:
 
     /// Backward DFT
     // Real valued IDFT
-    void applyIDFT(const FComplex<FReal>* transformedData,
+    void applyIDFT(const stdComplex<FReal>* transformedData,
                    FReal* sampledData) const
     {
         // read transformed data
@@ -171,8 +174,8 @@ public:
     }
 
     // Complexe valued IDFT
-    void applyIDFT(const FComplex<FReal>* transformedData,
-                   FComplex<FReal>* sampledData) const
+    void applyIDFT(const stdComplex<FReal>* transformedData,
+                   stdComplex<FReal>* sampledData) const
     {
         // read transformed data
         FBlas::c_setzero(size_,reinterpret_cast<FReal*>(data_));
@@ -182,7 +185,7 @@ public:
         // perform direct backward transformation
         for(unsigned int r=0; r<size_; ++r){
             for(unsigned int s=0; s<size_; ++s){
-                data_[r] += FComplex<FReal>(dataF_[s].getReal()*cosRS_[r*size_+s]
+                data_[r] += stdComplex<FReal>(dataF_[s].getReal()*cosRS_[r*size_+s]
                         - dataF_[s].getImag()*sinRS_[r*size_+s],
                         dataF_[s].getImag()*cosRS_[r*size_+s]
                         + dataF_[s].getReal()*sinRS_[r*size_+s]);
@@ -225,10 +228,10 @@ protected:
             dest[idxVal] += src[idxVal];
         }
     }
-    static void PlusEqual(FComplex<double> dest[], const fftw_complex src[], const int nbElements){
+    static void PlusEqual(stdComplex<double> dest[], const fftw_complex src[], const int nbElements){
         const double* ptrSrc = reinterpret_cast<const double*>(src);
         for(int idxVal = 0 ; idxVal < nbElements ; ++idxVal){
-            dest[idxVal] += FComplex<double>(ptrSrc[2*idxVal], ptrSrc[2*idxVal+1]);
+            dest[idxVal] += stdComplex<double>(ptrSrc[2*idxVal], ptrSrc[2*idxVal+1]);
         }
     }
     static void Equal(double dest[], const double src[], const int nbElements){
@@ -236,10 +239,10 @@ protected:
             dest[idxVal] = src[idxVal];
         }
     }
-    static void Equal(FComplex<double> dest[], const fftw_complex src[], const int nbElements){
+    static void Equal(stdComplex<double> dest[], const fftw_complex src[], const int nbElements){
         const double* ptrSrc = reinterpret_cast<const double*>(src);
         for(int idxVal = 0 ; idxVal < nbElements ; ++idxVal){
-            dest[idxVal] = FComplex<double>(ptrSrc[2*idxVal], ptrSrc[2*idxVal+1]);
+            dest[idxVal] = stdComplex<double>(ptrSrc[2*idxVal], ptrSrc[2*idxVal+1]);
         }
     }
     static void Bind_fftw_alloc(double** ptr, const int size){
@@ -291,10 +294,10 @@ protected:
             dest[idxVal] += src[idxVal];
         }
     }
-    static void PlusEqual(FComplex<float> dest[], const fftwf_complex src[], const int nbElements){
+    static void PlusEqual(stdComplex<float> dest[], const fftwf_complex src[], const int nbElements){
         const float* ptrSrc = reinterpret_cast<const float*>(src);
         for(int idxVal = 0 ; idxVal < nbElements ; ++idxVal){
-            dest[idxVal] += FComplex<float>(ptrSrc[2*idxVal], ptrSrc[2*idxVal+1]);
+            dest[idxVal] += stdComplex<float>(ptrSrc[2*idxVal], ptrSrc[2*idxVal+1]);
         }
     }
     static void Equal(float dest[], const float src[], const int nbElements){
@@ -302,10 +305,10 @@ protected:
             dest[idxVal] = src[idxVal];
         }
     }
-    static void Equal(FComplex<float> dest[], const fftwf_complex src[], const int nbElements){
+    static void Equal(stdComplex<float> dest[], const fftwf_complex src[], const int nbElements){
         const float* ptrSrc = reinterpret_cast<const float*>(src);
         for(int idxVal = 0 ; idxVal < nbElements ; ++idxVal){
-            dest[idxVal] = FComplex<float>(ptrSrc[2*idxVal], ptrSrc[2*idxVal+1]);
+            dest[idxVal] = stdComplex<float>(ptrSrc[2*idxVal], ptrSrc[2*idxVal+1]);
         }
     }
     static void Bind_fftw_alloc(float** ptr, const int size){
@@ -537,7 +540,7 @@ public:
             resultSignal[idxVal] /= realNbPoints;
         }
     }
-    void normalize(FComplex<double>* resultSignal) const {
+    void normalize(stdComplex<double>* resultSignal) const {
         const double realNbPoints = static_cast<double>(nbPoints);
         for(int idxVal = 0 ; idxVal < nbPoints ; ++idxVal){
             resultSignal[idxVal] /= realNbPoints;
@@ -549,7 +552,7 @@ public:
             resultSignal[idxVal] /= realNbPoints;
         }
     }
-    void normalize(FComplex<float>* resultSignal) const {
+    void normalize(stdComplex<float>* resultSignal) const {
         const float realNbPoints = static_cast<float>(nbPoints);
         for(int idxVal = 0 ; idxVal < nbPoints ; ++idxVal){
             resultSignal[idxVal] /= realNbPoints;
@@ -572,14 +575,14 @@ class FFftw;
 /// Double precision
 //////////////////////////////////////////////////////////////////////////////
 template <int DIM>
-class FFftw <double, FComplex<double>, DIM> : public FFftwCore<double, FComplex<double>, double, fftw_complex, fftw_plan, DIM>{
-    typedef FFftwCore<double, FComplex<double>, double, fftw_complex, fftw_plan, DIM> ParentClass;
+class FFftw <double, stdComplex<double>, DIM> : public FFftwCore<double, stdComplex<double>, double, fftw_complex, fftw_plan, DIM>{
+    typedef FFftwCore<double, stdComplex<double>, double, fftw_complex, fftw_plan, DIM> ParentClass;
 public:
     using ParentClass::ParentClass;
 };
 template <int DIM>
-class FFftw <FComplex<double>, FComplex<double>, DIM> : public FFftwCore<FComplex<double>, FComplex<double>, fftw_complex, fftw_complex, fftw_plan, DIM>{
-    typedef FFftwCore<FComplex<double>, FComplex<double>, fftw_complex, fftw_complex, fftw_plan, DIM> ParentClass;
+class FFftw <stdComplex<double>, stdComplex<double>, DIM> : public FFftwCore<stdComplex<double>, stdComplex<double>, fftw_complex, fftw_complex, fftw_plan, DIM>{
+    typedef FFftwCore<stdComplex<double>, stdComplex<double>, fftw_complex, fftw_complex, fftw_plan, DIM> ParentClass;
 public:
     using ParentClass::ParentClass;
 };
@@ -587,14 +590,14 @@ public:
 /// Single precision
 //////////////////////////////////////////////////////////////////////////////
 template <int DIM>
-class FFftw <float, FComplex<float>, DIM> : public FFftwCore<float, FComplex<float>, float, fftwf_complex, fftwf_plan, DIM>{
-    typedef FFftwCore<float, FComplex<float>, float, fftwf_complex, fftwf_plan, DIM> ParentClass;
+class FFftw <float, stdComplex<float>, DIM> : public FFftwCore<float, stdComplex<float>, float, fftwf_complex, fftwf_plan, DIM>{
+    typedef FFftwCore<float, stdComplex<float>, float, fftwf_complex, fftwf_plan, DIM> ParentClass;
 public:
     using ParentClass::ParentClass;
 };
 template <int DIM>
-class FFftw <FComplex<float>, FComplex<float>, DIM> : public FFftwCore<FComplex<float>, FComplex<float>, fftwf_complex, fftwf_complex, fftwf_plan, DIM>{
-    typedef FFftwCore<FComplex<float>, FComplex<float>, fftwf_complex, fftwf_complex, fftwf_plan, DIM> ParentClass;
+class FFftw <stdComplex<float>, stdComplex<float>, DIM> : public FFftwCore<stdComplex<float>, stdComplex<float>, fftwf_complex, fftwf_complex, fftwf_plan, DIM>{
+    typedef FFftwCore<stdComplex<float>, stdComplex<float>, fftwf_complex, fftwf_complex, fftwf_plan, DIM> ParentClass;
 public:
     using ParentClass::ParentClass;
 };
diff --git a/Src/Utils/stdComplex.hpp b/Src/Utils/stdComplex.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..4811beb463783c3904ad36fa2645c72180614c5b
--- /dev/null
+++ b/Src/Utils/stdComplex.hpp
@@ -0,0 +1,26 @@
+// ===================================================================================
+// olivier.coulaud@inria.fr, berenger.bramas@inria.fr
+// This software is a computer program whose purpose is to compute the FMM.
+//
+// This software is governed by the CeCILL-C and LGPL licenses and
+// abiding by the rules of distribution of free software.  
+// 
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public and CeCILL-C Licenses for more details.
+// "http://www.cecill.info". 
+// "http://www.gnu.org/licenses".
+// ===================================================================================
+#ifndef STDCOMPLEXE_HPP
+#define STDCOMPLEXE_HPP
+
+#include <complex>
+template<typename T>
+using stdComplex = std::complex<T> ;
+
+
+
+#endif //STDCOMPLEXE_HPP
+
+
diff --git a/Tests/Utils/testFFT.cpp b/Tests/Utils/testFFT.cpp
index 1795f9fbc640378d566aeb83d6bd22d38b902d26..81517cee8ff1d366150bf455758bd595086e1b9a 100644
--- a/Tests/Utils/testFFT.cpp
+++ b/Tests/Utils/testFFT.cpp
@@ -28,7 +28,8 @@
 #endif
 
 #include "Utils/FGlobal.hpp"
-#include "Utils/FComplex.hpp"
+#include "Utils/FMath.hpp"
+#include "Utils/stdComplex.hpp"
 
 #include "Utils/FTic.hpp"
 
@@ -59,9 +60,9 @@ int main(int argc, char** argv)
 
     // input/output types
     typedef FReal           FftR2CInputType;
-    typedef FComplex<FReal> FftR2COutputType;
-    typedef FComplex<FReal> FftC2CInputType;
-    typedef FComplex<FReal> FftC2COutputType;
+    typedef stdComplex<FReal> FftR2COutputType;
+    typedef stdComplex<FReal> FftC2CInputType;
+    typedef stdComplex<FReal> FftC2COutputType;
 
     // fftw arrays
     FftR2CInputType*  FftR2CInput  = new FftR2CInputType[total_size];
@@ -74,8 +75,8 @@ int main(int argc, char** argv)
 
 
     // fftw wrappers
-    typedef FFftw<FReal          ,FComplex<FReal>,dim> FftwR2CClass;
-    typedef FFftw<FComplex<FReal>,FComplex<FReal>,dim> FftwC2CClass;
+    typedef FFftw<FReal          ,stdComplex<FReal>,dim> FftwR2CClass;
+    typedef FFftw<stdComplex<FReal>,stdComplex<FReal>,dim> FftwC2CClass;
 
     std::cout<< "Init FFTW wrappers: ";
     time.tic();
@@ -160,8 +161,8 @@ int main(int argc, char** argv)
     FReal* FftC2CInvOutReal = new FReal[total_size]; FReal* FftC2CInputReal = new FReal[total_size];
     FReal* FftC2CInvOutImag = new FReal[total_size]; FReal* FftC2CInputImag = new FReal[total_size];
     for( int s=0; s<total_size; ++s){
-        FftC2CInvOutReal[s] = FftC2CInvOut[s].getReal(); FftC2CInputReal[s] = FftC2CInput[s].getReal();
-        FftC2CInvOutImag[s] = FftC2CInvOut[s].getImag(); FftC2CInputImag[s] = FftC2CInput[s].getImag();
+        FftC2CInvOutReal[s] = FftC2CInvOut[s].real(); FftC2CInputReal[s] = FftC2CInput[s].real();
+        FftC2CInvOutImag[s] = FftC2CInvOut[s].imag(); FftC2CInputImag[s] = FftC2CInput[s].imag();
     }
 
     if(FMath::FAccurater<FReal>(FftR2CInput,     FftR2CInvOut,     total_size).getRelativeInfNorm() > threshold) isWorking=false;
diff --git a/Tests/Utils/testFastDiscreteConvolution.cpp b/Tests/Utils/testFastDiscreteConvolution.cpp
index 8d09cb011ba0ae7bff9a2dd2695c1a4375b6d600..150ffda447131dec9bd99e0570af000ce5f82a6e 100644
--- a/Tests/Utils/testFastDiscreteConvolution.cpp
+++ b/Tests/Utils/testFastDiscreteConvolution.cpp
@@ -33,7 +33,7 @@
 #include "Utils/FBlas.hpp" 
 
 // 
-#include "Utils/FComplex.hpp"
+#include "Utils/stdComplex.hpp"
 #include "Utils/FDft.hpp"
 
 #include "Utils/FParameterNames.hpp"
@@ -101,9 +101,9 @@ int main(int argc, char ** argv){
     std::cout<<std::endl;
 
     // now compute via DFT and use convolution theorem
-    FComplex<FReal> FK[N];
-    FComplex<FReal> FY[N];
-    FComplex<FReal> FX[N];
+    stdComplex<FReal> FK[N];
+    stdComplex<FReal> FY[N];
+    stdComplex<FReal> FX[N];
     FReal iFX[N];
 
     // Init DFTor
@@ -111,12 +111,12 @@ int main(int argc, char ** argv){
     time.tic();
     const int dim = 1;
     //FDft<FReal> Dft(N);// direct version (Beware! Ordering of output differs from REAL valued-FFT)
-    FFftw<FReal,FComplex<FReal>,dim> Dft(N);// fast version
+    FFftw<FReal,stdComplex<FReal>,dim> Dft(N);// fast version
     std::cout << "took " << time.tacAndElapsed() << "sec." << std::endl;
 
     // Initialize manually
     //  for(unsigned int s=0; s<N; ++s){
-    //    FX[s] = FY[s] = FK[s] =FComplex<FReal>(0.0,0.0); // init
+    //    FX[s] = FY[s] = FK[s] =stdComplex<FReal>(0.0,0.0); // init
     //    iFX[s]=0.0;
     //  }
     // ... or using Blas routines
diff --git a/Tests/Utils/testUnifInterpolator.cpp b/Tests/Utils/testUnifInterpolator.cpp
index acf6294b09329fc1d31aaf68f1691690862f97be..4ab0a637220a516ee0b7acfafd430d37a3c307ba 100644
--- a/Tests/Utils/testUnifInterpolator.cpp
+++ b/Tests/Utils/testUnifInterpolator.cpp
@@ -42,7 +42,7 @@
 #include "Kernels/Uniform/FUnifTensor.hpp"
 
 // Check DFT
-#include "Utils/FComplex.hpp"
+#include "Utils/stdComplex.hpp"
 #include "Utils/FDft.hpp"
 
 
@@ -333,11 +333,11 @@ int main(int argc, char ** argv){
 
     /////////////////////////////////////////////////////////////////////////////////////
     // Efficient application of the Toeplitz system in FOURIER SPACE
-    FComplex<FReal> FPMultExp[rc];
-    FComplex<FReal> FPLocalExp[rc];
+    stdComplex<FReal> FPMultExp[rc];
+    stdComplex<FReal> FPLocalExp[rc];
     FReal PLocalExp[rc];
 
-    //for (unsigned int n=0; n<rc; ++n) FPLocalExp[n]=FComplex<FReal>(0.0,0.0);
+    //for (unsigned int n=0; n<rc; ++n) FPLocalExp[n]=stdComplex<FReal>(0.0,0.0);
     FBlas::c_setzero(rc,reinterpret_cast<FReal*>(FPLocalExp));
 
     FBlas::setzero(rc,PLocalExp);
@@ -345,7 +345,7 @@ int main(int argc, char ** argv){
     // Init DFT
     const int dimfft = 1;
     //FDft Dft(rc); // direct version
-    FFftw<FReal,FComplex<FReal>,dimfft> Dft(rc); // fast version
+    FFftw<FReal,stdComplex<FReal>,dimfft> Dft(rc); // fast version
 
     // Get first COLUMN of K and Store in T
     FReal T[rc];
@@ -372,8 +372,8 @@ int main(int argc, char ** argv){
     //  std::cout<<std::endl;
 
     // Apply DFT to T
-    FComplex<FReal> FT[rc];
-    //  for (unsigned int n=0; n<rc; ++n) FT[n]=FComplex<FReal>(0.0,0.0);
+    stdComplex<FReal> FT[rc];
+    //  for (unsigned int n=0; n<rc; ++n) FT[n]=stdComplex<FReal>(0.0,0.0);
     FBlas::c_setzero(rc,reinterpret_cast<FReal*>(FT));
 
     // if first COLUMN (T) of C is used
@@ -397,7 +397,7 @@ int main(int argc, char ** argv){
 
 
     // Set transformed MultExp to 0
-    //  for (unsigned int n=0; n<rc; ++n) FPMultExp[n]=FComplex<FReal>(0.0,0.0);
+    //  for (unsigned int n=0; n<rc; ++n) FPMultExp[n]=stdComplex<FReal>(0.0,0.0);
     FBlas::c_setzero(rc,reinterpret_cast<FReal*>(FPMultExp));
 
     // Transform PaddedMultExp
diff --git a/Tests/Utils/testUnifTensorialInterpolator.cpp b/Tests/Utils/testUnifTensorialInterpolator.cpp
index 84249309317b2a70d448568192a4641fbbd1124f..f6ba6612e63bd256ad59bcb32fe31ecaaa308fc5 100644
--- a/Tests/Utils/testUnifTensorialInterpolator.cpp
+++ b/Tests/Utils/testUnifTensorialInterpolator.cpp
@@ -42,7 +42,7 @@
 #include "Kernels/Uniform/FUnifTensor.hpp"
 
 // Check DFT
-#include "Utils/FComplex.hpp"
+#include "Utils/stdComplex.hpp"
 #include "Utils/FDft.hpp"
 
 
@@ -362,7 +362,7 @@ int main(int argc, char ** argv){
     // Init DFT
     const int dimfft = 1;
     //FDft Dft(rc); // direct version
-    FFftw<FReal,FComplex<FReal>,dimfft> Dft(rc); // fast version
+    FFftw<FReal,stdComplex<FReal>,dimfft> Dft(rc); // fast version
 
     // Get first COLUMN of K and Store in T
     FReal T[ncmp*rc];
@@ -390,8 +390,8 @@ int main(int argc, char ** argv){
     //  std::cout<<std::endl;
 
     // Apply DFT to T
-    FComplex<FReal> FT[ncmp*rc];
-    //  for (unsigned int n=0; n<rc; ++n) FT[n]=FComplex<FReal>(0.0,0.0);
+    stdComplex<FReal> FT[ncmp*rc];
+    //  for (unsigned int n=0; n<rc; ++n) FT[n]=stdComplex<FReal>(0.0,0.0);
     FBlas::c_setzero(ncmp*rc,reinterpret_cast<FReal*>(FT));
 
     // if first COLUMN (T) of C is used
@@ -400,11 +400,11 @@ int main(int argc, char ** argv){
     //  // if first ROW of C is used
     //  Dft.applyDFT(C,FT);
 
-    FComplex<FReal> FPMultExp[nrhs*rc];
-    FComplex<FReal> FPLocalExp[nlhs*rc];
+    stdComplex<FReal> FPMultExp[nrhs*rc];
+    stdComplex<FReal> FPLocalExp[nlhs*rc];
     FReal PLocalExp[nlhs*rc];
 
-    //for (unsigned int n=0; n<rc; ++n) FPLocalExp[n]=FComplex<FReal>(0.0,0.0);
+    //for (unsigned int n=0; n<rc; ++n) FPLocalExp[n]=stdComplex<FReal>(0.0,0.0);
     FBlas::c_setzero(nlhs*rc,reinterpret_cast<FReal*>(FPLocalExp));
 
     FBlas::setzero(nlhs*rc,PLocalExp);
@@ -426,7 +426,7 @@ int main(int argc, char ** argv){
 
 
     // Set transformed MultExp to 0
-    //  for (unsigned int n=0; n<rc; ++n) FPMultExp[n]=FComplex<FReal>(0.0,0.0);
+    //  for (unsigned int n=0; n<rc; ++n) FPMultExp[n]=stdComplex<FReal>(0.0,0.0);
     FBlas::c_setzero(nrhs*rc,reinterpret_cast<FReal*>(FPMultExp));
 
     // Transform PaddedMultExp
@@ -437,7 +437,7 @@ int main(int argc, char ** argv){
     time.tic();
 
     // Application of M2L in FOURIER SPACE
-    FComplex<FReal> tmpFX;
+    stdComplex<FReal> tmpFX;
     for (unsigned int idxLhs=0; idxLhs<nlhs; ++idxLhs){
         unsigned int idxRhs = idxLhs % npot;
         unsigned int d = MatrixKernel.getPosition(idxLhs);