Had clean the Blas version

+
The blas fonction (only one yet) used is in a global file
+
A simple test to compare Targets and Sources model (Tsm)
and no-Tsm has been added to the Tests directory.
It simply checks that all cells have the same values (and it's right)
But it does not test the L2P and P2P.

git-svn-id: svn+ssh://scm.gforge.inria.fr/svn/scalfmm/scalfmm/trunk@42 2616d619-271b-44dc-8df4-d4a8f33a7222

Sources/Fmb/FAbstractFmbKernels.hpp

@@ -84,10 +84,6 @@ protected:
     // p_Y_theta_derivated
     FComplexe current_thread_Y_theta_derivated[FMB_Info_exp_size];
 
-    // [OK?] p_precomputed_cos_and_sin_array
-    FComplexe cosSin[FMB_Info_M2L_P + 1];
-    // [OK?] p_associated_Legendre_function_Array
-    FReal legendre[FMB_Info_M2L_exp_size];
 
     // pow_of_I_array
     static const FReal PiArrayInner[4];
@@ -282,19 +278,24 @@ protected:
     //2.7 these
     void harmonicInner(const Spherical& inSphere, FComplexe* const outResults){
 
-        FComplexe* ptrCosSin = this->cosSin;
-        for(int idxl = 0 , idxlMod4 = 0; idxl <= FMB_Info_P ; ++idxl, ++idxlMod4, ++ptrCosSin){
+        // p_precomputed_cos_and_sin_array
+        FComplexe cosSin[FMB_Info_M2L_P + 1];
+
+        for(int idxl = 0 , idxlMod4 = 0; idxl <= FMB_Info_P ; ++idxl, ++idxlMod4){
             if(idxlMod4 == 4) idxlMod4 = 0;
             const FReal angleinter = FReal(idxl) * inSphere.phi;
             const FReal angle = angleinter + this->PiArrayInner[idxlMod4];
 
-            ptrCosSin->setReal( FMath::Sin(angle + FMath::FPiDiv2) );
-            ptrCosSin->setImag( FMath::Sin(angle) );
+            cosSin[idxl].setReal( FMath::Sin(angle + FMath::FPiDiv2) );
+            cosSin[idxl].setImag( FMath::Sin(angle) );
 
-            //printf("%d=%f/%f (%d/%f/%f)\n",idxl,ptrCosSin->getReal(),ptrCosSin->getImag(),idxl,inSphere.phi,this->PiArrayInner[idxlMod4]);
+            //printf("%d=%f/%f (%d/%f/%f)\n",idxl,cosSin[idxl].getReal(),cosSin[idxl].getImag(),idxl,inSphere.phi,this->PiArrayInner[idxlMod4]);
         }
 
-        legendreFunction(FMB_Info_P,inSphere.cosTheta, inSphere.sinTheta, this->legendre);
+        // p_associated_Legendre_function_Array
+        FReal legendre[FMB_Info_M2L_exp_size];
+
+        legendreFunction(FMB_Info_P,inSphere.cosTheta, inSphere.sinTheta, legendre);
         /*printf("FMB_Info_M2L_exp_size=%d\n",FMB_Info_M2L_exp_size);
         for(int temp = 0 ; temp < FMB_Info_M2L_exp_size ; ++temp){
             printf("%f\n",this->legendre[temp]);
@@ -309,8 +310,8 @@ protected:
         for(int idxl = 0; idxl <= FMB_Info_P ; ++idxl, idxRl *= inSphere.r){
             for(int idxm = 0 ; idxm <= idxl ; ++idxm, ++currentResult, ++idxSphereHarmoCoef, ++idxLegendre){
                 const FReal magnitude = this->sphereHarmoInnerCoef[idxSphereHarmoCoef] * idxRl * legendre[idxLegendre];
-                currentResult->setReal( magnitude * this->cosSin[idxm].getReal() );
-                currentResult->setImag( magnitude * this->cosSin[idxm].getImag() );
+                currentResult->setReal( magnitude * cosSin[idxm].getReal() );
+                currentResult->setImag( magnitude * cosSin[idxm].getImag() );
 
                 //printf("l = %d m = %d\n",idxl,idxm);
                 //printf("magnitude=%f idxRl=%f sphereHarmoInnerCoef=%f real=%f imag=%f\n",magnitude,idxRl,this->sphereHarmoInnerCoef[idxSphereHarmoCoef],currentResult->getReal(),currentResult->getImag());
@@ -320,20 +321,23 @@ protected:
     }
     // spherical_harmonic_Outer
     void harmonicOuter(const Spherical& inSphere, FComplexe* const outResults){
-
-        FComplexe* ptrCosSin = this->cosSin;
-        for(int idxl = 0, idxlMod4 = 0; idxl <= FMB_Info_M2L_P ; ++idxl, ++idxlMod4, ++ptrCosSin){
+        // p_precomputed_cos_and_sin_array
+        FComplexe cosSin[FMB_Info_M2L_P + 1];
+        for(int idxl = 0, idxlMod4 = 0; idxl <= FMB_Info_M2L_P ; ++idxl, ++idxlMod4){
             if(idxlMod4 == 4) idxlMod4 = 0;
             const FReal angle = idxl * inSphere.phi + this->PiArrayOuter[idxlMod4];
 
-            ptrCosSin->setReal( FMath::Sin(angle + FMath::FPiDiv2) );
-            ptrCosSin->setImag( FMath::Sin(angle) );
+            cosSin[idxl].setReal( FMath::Sin(angle + FMath::FPiDiv2) );
+            cosSin[idxl].setImag( FMath::Sin(angle) );
 
             //printf("l=%d \t inSphere.phi=%f \t this->PiArrayOuter[idxlMod4]=%f \t angle=%f \t FMath::Sin(angle + FMath::FPiDiv2)=%f \t FMath::Sin(angle)=%f\n",
             //        idxl, inSphere.phi, this->PiArrayOuter[idxlMod4], angle, FMath::Sin(angle + FMath::FPiDiv2) , FMath::Sin(angle));
         }
 
-        legendreFunction(FMB_Info_M2L_P,inSphere.cosTheta, inSphere.sinTheta, this->legendre);
+        // p_associated_Legendre_function_Array
+        FReal legendre[FMB_Info_M2L_exp_size];
+
+        legendreFunction(FMB_Info_M2L_P,inSphere.cosTheta, inSphere.sinTheta, legendre);
 
         int idxLegendre = 0;
         FComplexe* currentResult = outResults;
@@ -342,12 +346,12 @@ protected:
         FReal idxRl1 = invR;
         for(int idxl = 0 ; idxl <= FMB_Info_M2L_P ; ++idxl, idxRl1 *= invR){
             for(int idxm = 0 ; idxm <= idxl ; ++idxm, ++currentResult, ++idxLegendre){
-                const FReal magnitude = this->sphereHarmoOuterCoef[idxl-idxm] * idxRl1 * this->legendre[idxLegendre];
-                currentResult->setReal( magnitude * this->cosSin[idxm].getReal() );
-                currentResult->setImag( magnitude * this->cosSin[idxm].getImag() );
+                const FReal magnitude = this->sphereHarmoOuterCoef[idxl-idxm] * idxRl1 * legendre[idxLegendre];
+                currentResult->setReal( magnitude * cosSin[idxm].getReal() );
+                currentResult->setImag( magnitude * cosSin[idxm].getImag() );
                 //printf("l=%d\t m=%d\t idxRl1=%f\t magnitude=%f\n",idxl,idxm,idxRl1,magnitude);
-                //printf("l=%d\t m=%d\t this->cosSin[idxm].getReal()=%f\t this->cosSin[idxm].getImag()=%f\n",
-                //       idxl,idxm,this->cosSin[idxm].getReal(),this->cosSin[idxm].getImag());
+                //printf("l=%d\t m=%d\t cosSin[idxm].getReal()=%f\t cosSin[idxm].getImag()=%f\n",
+                //       idxl,idxm,cosSin[idxm].getReal(),cosSin[idxm].getImag());
             }
         }
     }
@@ -387,12 +391,8 @@ protected:
                          FComplexe * results_array,
                          FComplexe * theta_derivated_results_array
                          ){
-        FComplexe *p_term = results_array;
-        FComplexe *p_theta_derivated_term = theta_derivated_results_array;
-        FReal *p_spherical_harmonic_Inner_coefficients_array = sphereHarmoInnerCoef;
-        FReal *ptr_associated_Legendre_function_Array = legendre;
-        FReal *start_ptr_associated_Legendre_function_Array = ptr_associated_Legendre_function_Array;
-
+        // p_precomputed_cos_and_sin_array
+        FComplexe cosSin[FMB_Info_M2L_P + 1];
 
         //printf("HarmoInnerTheta \t lmax = %d \t r = %f \t cos_theta = %f \t sin_theta = %f \t phi = %f\n",
         //       FMB_Info_P,inSphere.r,inSphere.cosTheta,inSphere.sinTheta,inSphere.phi);
@@ -401,15 +401,26 @@ protected:
         for(int idxm = 0 , idxmMod4 = 0; idxm <= FMB_Info_P ; ++idxm, ++idxmMod4){
             if(idxmMod4 == 4) idxmMod4 = 0;
             const FReal angle = idxm*inSphere.phi + PiArrayInner[idxmMod4];
-            (cosSin + idxm)->setReal(FMath::Sin(angle + FMath::FPiDiv2));
-            (cosSin + idxm)->setImag(FMath::Sin(angle));
+            cosSin[idxm].setReal(FMath::Sin(angle + FMath::FPiDiv2));
+            cosSin[idxm].setImag(FMath::Sin(angle));
 
             //printf("l=%d \t inSphere.phi=%f \t this->PiArrayOuter[idxlMod4]=%f \t angle=%f \t FMath::Sin(angle + FMath::FPiDiv2)=%f \t FMath::Sin(angle)=%f\n",
             //        idxm, inSphere.phi, this->PiArrayInner[idxmMod4], angle, FMath::Sin(angle + FMath::FPiDiv2) , FMath::Sin(angle));
         }
 
+
+        // p_associated_Legendre_function_Array
+        FReal legendre[FMB_Info_M2L_exp_size];
+
         // Initialization of associated_Legendre_function_Array:
-        legendreFunction(FMB_Info_P, inSphere.cosTheta, inSphere.sinTheta, this->legendre);
+        legendreFunction(FMB_Info_P, inSphere.cosTheta, inSphere.sinTheta, legendre);
+
+
+        FComplexe *p_term = results_array;
+        FComplexe *p_theta_derivated_term = theta_derivated_results_array;
+        FReal *p_spherical_harmonic_Inner_coefficients_array = sphereHarmoInnerCoef;
+        FReal *ptr_associated_Legendre_function_Array = legendre;
+        FReal *start_ptr_associated_Legendre_function_Array = ptr_associated_Legendre_function_Array;
 
         // r^l
         FReal r_l = 1.0;
@@ -421,14 +432,14 @@ protected:
                 magnitude = (*p_spherical_harmonic_Inner_coefficients_array) * r_l * (*ptr_associated_Legendre_function_Array);
 
                 // Computation of Inner_l^m(r, theta, phi):
-                p_term->setReal( magnitude * (cosSin + m)->getReal());
-                p_term->setImag( magnitude * (cosSin + m)->getImag());
+                p_term->setReal( magnitude * cosSin[m].getReal());
+                p_term->setImag( magnitude * cosSin[m].getImag());
 
                 /*printf("%d/%d - magnitude=%e ptr_precomputed_cos_and_sin_array real=%e imag=%e p_term real=%e imag=%e\n",
                        l,m,
                        magnitude,
-                       (cosSin + m)->getReal(),
-                       (cosSin + m)->getImag(),
+                       cosSin[m].getReal(),
+                       cosSin[m].getImag(),
                        p_term->getReal(),
                        p_term->getImag());*/
                 /*printf("\t p_spherical_harmonic_Inner_coefficients_array = %e \t ptr_associated_Legendre_function_Array = %e \t r_l = %e\n",
@@ -440,8 +451,8 @@ protected:
                 // Computation of {\partial Inner_l^m(r, theta, phi)}/{\partial theta}:
                 magnitude = (*p_spherical_harmonic_Inner_coefficients_array) * r_l * ((l*inSphere.cosTheta*(*ptr_associated_Legendre_function_Array)
                                         - (l+m)*(*(start_ptr_associated_Legendre_function_Array + expansion_Redirection_array_for_j[l-1] + m))) / inSphere.sinTheta);
-                p_theta_derivated_term->setReal(magnitude * (cosSin + m)->getReal());
-                p_theta_derivated_term->setImag(magnitude * (cosSin + m)->getImag());
+                p_theta_derivated_term->setReal(magnitude * cosSin[m].getReal());
+                p_theta_derivated_term->setImag(magnitude * cosSin[m].getImag());
 
                 //printf("magnitude=%f r_l=%f p_spherical_harmonic_Inner_coefficients_array=%f real=%f imag=%f\n",
                 //       magnitude,r_l,*p_spherical_harmonic_Inner_coefficients_array,p_theta_derivated_term->getReal(),p_theta_derivated_term->getImag());
@@ -450,14 +461,14 @@ protected:
             // m=l:
             // Computation of Inner_m^m(r, theta, phi):
             magnitude = (*p_spherical_harmonic_Inner_coefficients_array) * r_l * (*ptr_associated_Legendre_function_Array);
-            p_term->setReal(magnitude * (cosSin + m)->getReal());
-            p_term->setImag(magnitude * (cosSin + m)->getImag());
+            p_term->setReal(magnitude * cosSin[m].getReal());
+            p_term->setImag(magnitude * cosSin[m].getImag());
 
             /*printf("%d - magnitude=%e ptr_precomputed_cos_and_sin_array real=%e imag=%e p_term real=%e imag=%e\n",
                    l,
                    magnitude,
-                   (cosSin + m)->getReal(),
-                   (cosSin + m)->getImag(),
+                   cosSin[m].getReal(),
+                   cosSin[m].getImag(),
                    p_term->getReal(),
                    p_term->getImag());*/
             /*printf("\t p_spherical_harmonic_Inner_coefficients_array = %e \t ptr_associated_Legendre_function_Array = %e \t r_l = %e\n",
@@ -468,8 +479,8 @@ protected:
 
             // Computation of {\partial Inner_m^m(r, theta, phi)}/{\partial theta}:
             magnitude = (*p_spherical_harmonic_Inner_coefficients_array) * r_l * (m * inSphere.cosTheta * (*ptr_associated_Legendre_function_Array) / inSphere.sinTheta);
-            p_theta_derivated_term->setReal(magnitude * (cosSin + m)->getReal());
-            p_theta_derivated_term->setImag(magnitude * (cosSin + m)->getImag());
+            p_theta_derivated_term->setReal(magnitude * cosSin[m].getReal());
+            p_theta_derivated_term->setImag(magnitude * cosSin[m].getImag());
 
             //printf("magnitude=%f r_l=%f p_spherical_harmonic_Inner_coefficients_array=%f real=%f imag=%f\n",
             //       magnitude,r_l,*p_spherical_harmonic_Inner_coefficients_array,p_theta_derivated_term->getReal(),p_theta_derivated_term->getImag());
@@ -662,7 +673,7 @@ public:
 
             harmonicInner(positionTsmphere(iterParticule.value()->getPosition() - inPole->getPosition()),current_thread_Y);
 
-            //ok printf("\tr=%f\tcos_theta=%f\tsin_theta=%f\tphi=%f\n",spherical.r,spherical.cosTheta,spherical.sinTheta,spherical.phi);
+            //printf("\tr=%f\tcos_theta=%f\tsin_theta=%f\tphi=%f\n",spherical.r,spherical.cosTheta,spherical.sinTheta,spherical.phi);
 
             FComplexe* p_exp_term = inPole->getMultipole();
             FComplexe* p_Y_term = current_thread_Y;

Sources/Fmb/FExtendFmbCell.hpp

@@ -17,10 +17,11 @@
 * It is needed by the Fmb Kernels.
 */
 class FExtendFmbCell {
-protected:
+public:
     // FMB_Info_P is declared in FAbstractFmbKernels
     static const int MultipoleSize = int(((FMB_Info_P)+1) * ((FMB_Info_P)+2) * 0.5); //< The size of the multipole
 
+protected:
     FComplexe multipole_exp[MultipoleSize]; //< For multipole extenssion
     FComplexe local_exp[MultipoleSize];     //< For local extenssion
 

Sources/Fmb/FFmbKernelsBlas.hpp

@@ -3,6 +3,7 @@
 // /!\ Please, you must read the license at the bottom of this page
 
 #include "../Utils/FGlobal.hpp"
+#include "../Utils/FBlas.hpp"
 #include "../Components/FAbstractKernels.hpp"
 
 #include "../Containers/FTreeCoordinate.hpp"
@@ -14,36 +15,6 @@
 #include "../Utils/FTrace.hpp"
 
 #include <iostream>
-#include <cblas.h>
-
-
-template <typename T>
-void cblas_gemv(const enum CBLAS_ORDER order ,
-                   const enum CBLAS_TRANSPOSE TransA , const int M , const int N ,
-                   const void *alpha , const void *A , const int lda ,
-                   const void *X , const int incX , const void *beta ,
-                   void *Y , const int incY){
-    T t;
-    t.you_cannot_use_this_function_with_this_type();
-}
-
-template <>
-void cblas_gemv<double>(const enum CBLAS_ORDER order ,
-        const enum CBLAS_TRANSPOSE TransA , const int M , const int N ,
-        const void *alpha , const void *A , const int lda ,
-        const void *X , const int incX , const void *beta ,
-        void *Y , const int incY){
-    cblas_zgemv(order,TransA,M,N,alpha,A,lda,X,incX,beta,Y,incY);
-}
-
-template <>
-void cblas_gemv<float>(const enum CBLAS_ORDER order ,
-        const enum CBLAS_TRANSPOSE TransA , const int M , const int N ,
-        const void *alpha , const void *A , const int lda ,
-        const void *X , const int incX , const void *beta ,
-        void *Y , const int incY){
-    cblas_cgemv(order,TransA,M,N,alpha,A,lda,X,incX,beta,Y,incY);
-}
 
 
 /**
@@ -321,11 +292,12 @@ protected:
             cosSin[idxl].setReal( FMath::Sin(angle + FMath::FPiDiv2) );
             cosSin[idxl].setImag( FMath::Sin(angle) );
 
-            //printf("%d=%.15e/%.15e (angle %.15e=%d * %.15e + %.15e // (%.15e))\n",idxl,cosSin[idxl].getReal(),cosSin[idxl].getImag(),
-            //       angle,idxl,inSphere.phi,this->PiArrayInner[idxlMod4],angleinter);
+            /*printf("%d=%.15e/%.15e (angle %.15e=%d * %.15e + %.15e // (%.15e))\n",
+                   idxl,cosSin[idxl].getReal(),cosSin[idxl].getImag(),
+                   angle,idxl,inSphere.phi,this->PiArrayInner[idxlMod4],angleinter);*/
 
-            //printf("sin(%.15e)  = %.15e\n",
-            //       angle + FMath::FPiDiv2,FMath::Sin(angle + FMath::FPiDiv2));
+            /*printf("sin(%.15e)  = %.15e\n",
+                   angle + FMath::FPiDiv2,FMath::Sin(angle + FMath::FPiDiv2));*/
         }
 
         // p_associated_Legendre_function_Array

Sources/Utils/FBlas.hpp

+#ifndef FBLAS_HPP
+#define FBLAS_HPP
+
+#include <cblas.h>
+//#include <mkl.h>
+
+///////////////////////////////////////////////////////
+// GEMV
+///////////////////////////////////////////////////////
+
+template <typename T>
+void cblas_gemv(const CBLAS_ORDER order ,
+                   const CBLAS_TRANSPOSE TransA , const int M , const int N ,
+                   const void *alpha , const void *A , const int lda ,
+                   const void *X , const int incX , const void *beta ,
+                   void *Y , const int incY){
+    T t;
+    t.you_cannot_use_this_function_with_this_type();
+}
+
+template <>
+void cblas_gemv<double>(const CBLAS_ORDER order ,
+        const CBLAS_TRANSPOSE TransA , const int M , const int N ,
+        const void *alpha , const void *A , const int lda ,
+        const void *X , const int incX , const void *beta ,
+        void *Y , const int incY){
+    cblas_zgemv(order,TransA,M,N,alpha,A,lda,X,incX,beta,Y,incY);
+}
+
+template <>
+void cblas_gemv<float>(const CBLAS_ORDER order ,
+        const CBLAS_TRANSPOSE TransA , const int M , const int N ,
+        const void *alpha , const void *A , const int lda ,
+        const void *X , const int incX , const void *beta ,
+        void *Y , const int incY){
+    cblas_cgemv(order,TransA,M,N,alpha,A,lda,X,incX,beta,Y,incY);
+}
+
+#endif //FBLAS_HPP
+

Sources/Utils/FComplexe.hpp

@@ -10,6 +10,9 @@
 * Please read the license
 *
 * Propose basic complexe class.
+* Do not modify the attributes of this class.
+* It can be passed to blas fonction and has to be
+* 2 x real size only.
 */
 class FComplexe {
     FReal real;    //< Real

Tests/CMakeLists.txt

@@ -12,6 +12,7 @@ endif()
 
 #add blass
 set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -lcblas")
+#ADD_DEFINITIONS(-Wl,--start-group $(MKLROOT)/lib/ia32/libmkl_intel.a $(MKLROOT)/lib/ia32/libmkl_sequential.a $(MKLROOT)/lib/ia32/libmkl_core.a -Wl,--end-group -lpthread)
 
 #test
 # ADD_DEFINITIONS(-fexpensive-optimizations -funroll-loops -frerun-cse-after-loop -frerun-loop-opt)

Tests/testFmbBlasAlgorithm.cpp

@@ -29,9 +29,10 @@
 
 #include "../Sources/Files/FFMALoader.hpp"
 
-// With openmp : g++ testFmbBlasAlgorithm.cpp ../Sources/Utils/FAssertable.cpp ../Sources/Utils/FDebug.cpp ../Sources/Utils/FTrace.cpp -lgomp -fopenmp -lblas -O2 -o testFmbBlasAlgorithm.exe
+// With openmp : g++ testFmbBlasAlgorithm.cpp ../Sources/Utils/FAssertable.cpp ../Sources/Utils/FDebug.cpp ../Sources/Utils/FTrace.cpp -lgomp -fopenmp -lcblas -O2 -o testFmbBlasAlgorithm.exe
 // icpc -openmp -openmp-lib=compat testFmbAlgorithm.cpp ../Sources/Utils/FAssertable.cpp ../Sources/Utils/FDebug.cpp -O2 -o testFmbAlgorithm.exe
 
+// g++ testFmbBlasAlgorithm.cpp ../Sources/Utils/FAssertable.cpp ../Sources/Utils/FDebug.cpp ../Sources/Utils/FTrace.cpp -lgomp -fopenmp -lmkl_ia32 -O2 -o testFmbBlasAlgorithm.exe
 /** This program show an example of use of
   * the fmm basic algo
   * it also check that eachh particules is little or longer
@@ -116,9 +117,9 @@ int main(int argc, char ** argv){
     counter.tic();
 
     //FFmbKernelsBlas FFmbKernelsPotentialForces
-    FFmbKernelsPotentialForces<FmbParticule, FmbCell, NbLevels> kernels(loader.getBoxWidth());
+    FFmbKernelsBlas<FmbParticule, FmbCell, NbLevels> kernels(loader.getBoxWidth());
     //FFMMAlgorithm FFMMAlgorithmThreaded FFMMAlgorithmArray FFMMAlgorithmTask
-    FFmmAlgorithm<FFmbKernelsPotentialForces, FmbParticule, FmbCell, FSimpleLeaf, NbLevels, SizeSubLevels> algo(&tree,&kernels);
+    FFmmAlgorithm<FFmbKernelsBlas, FmbParticule, FmbCell, FSimpleLeaf, NbLevels, SizeSubLevels> algo(&tree,&kernels);
     algo.execute();
 
     counter.tac();

Tests/testFmbTsmNoTsm.cpp

+// /!\ Please, you must read the license at the bottom of this page
+
+#include <iostream>
+
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "../Sources/Utils/FTic.hpp"
+
+#include "../Sources/Containers/FOctree.hpp"
+#include "../Sources/Containers/FList.hpp"
+
+#include "../Sources/Components/FFmaParticule.hpp"
+#include "../Sources/Extenssions/FExtendForces.hpp"
+#include "../Sources/Extenssions/FExtendPotential.hpp"
+
+#include "../Sources/Extenssions/FExtendParticuleType.hpp"
+#include "../Sources/Extenssions/FExtendCellType.hpp"
+
+
+#include "../Sources/Components/FBasicCell.hpp"
+#include "../Sources/Fmb/FExtendFmbCell.hpp"
+
+#include "../Sources/Core/FFmmAlgorithm.hpp"
+#include "../Sources/Core/FFmmAlgorithmTsm.hpp"
+#include "../Sources/Core/FFmmAlgorithmArray.hpp"
+#include "../Sources/Core/FFmmAlgorithmArrayTsm.hpp"
+
+#include "../Sources/Components/FSimpleLeaf.hpp"
+#include "../Sources/Components/FTypedLeaf.hpp"
+
+#include "../Sources/Fmb/FFmbKernelsPotentialForces.hpp"
+#include "../Sources/Fmb/FFmbKernelsForces.hpp"
+#include "../Sources/Fmb/FFmbKernelsPotential.hpp"
+
+
+// With openmp : g++ testFmbTsmAlgorithm.cpp ../Sources/Utils/FAssertable.cpp ../Sources/Utils/FDebug.cpp ../Sources/Utils/FTrace.cpp -lgomp -fopenmp -O2 -o testFmbTsmAlgorithm.exe
+// icpc -openmp -openmp-lib=compat testFmbTsmAlgorithm.cpp ../Sources/Utils/FAssertable.cpp ../Sources/Utils/FDebug.cpp -O2 -o testFmbTsmAlgorithm.exe
+
+/** This program show an example of use of
+  * the fmm basic algo
+  * it also check that eachh particules is little or longer
+  * related that each other
+  */
+
+
+/** Fmb class has to extend {FExtendForces,FExtendPotential,FExtendPhysicalValue}
+  * Because we use fma loader it needs {FFmaParticule}
+  */
+class FmbParticule : public FFmaParticule, public FExtendForces, public FExtendPotential {
+public:
+};
+
+class FmbParticuleTyped : public FmbParticule, public FExtendParticuleType {
+public:
+};
+
+/** Custom cell
+  *
+  */
+class FmbCell : public FBasicCell, public FExtendFmbCell {
+public:
+};
+
+class FmbCellTyped : public FmbCell, public FExtendCellType {
+public:
+};
+
+// Simply create particules and try the kernels
+int main(int argc, char ** argv){
+    ///////////////////////What we do/////////////////////////////
+    std::cout << ">> This executable has to be used to test Fmb on a Tsm system.\n";
+    std::cout << ">> It compares the results between Tms and no Tms (except P2P & L2P).\n";
+    //////////////////////////////////////////////////////////////
+
+    const int NbLevels = 9;//10;
+    const int SizeSubLevels = 3;//3
+    FTic counter;
+    const long NbPart = 200000;//2000000
+    const double BoxWidth = 1.0;
+    const F3DPosition CenterOfBox(0.5,0.5,0.5);
+
+    // -----------------------------------------------------
+
+    std::cout << "Creating " << NbPart << " particules ..." << std::endl;
+    counter.tic();
+
+    FmbParticule* particules = new FmbParticule[NbPart];
+    FmbParticuleTyped* particulesTyped = new FmbParticuleTyped[NbPart*2];
+
+    for(int idxPart = 0 ; idxPart < NbPart ; ++idxPart){
+        const double x = FReal(rand())/RAND_MAX;
+        const double y = FReal(rand())/RAND_MAX;
+        const double z = FReal(rand())/RAND_MAX;
+        // Particule for standart model
+        particules[idxPart].setPosition(x,y,z);
+        particules[idxPart].setPhysicalValue(1);
+
+        // Create a clone for typed (Tsm) version
+        particulesTyped[idxPart*2].setPosition(x,y,z);
+        particulesTyped[idxPart*2+1].setPosition(x,y,z);
+
+        particulesTyped[idxPart*2].setPhysicalValue(1);
+        particulesTyped[idxPart*2+1].setPhysicalValue(1);
+
+        particulesTyped[idxPart*2].setAsSource();
+        particulesTyped[idxPart*2+1].setAsTarget();
+    }
+
+    counter.tac();
+    std::cout << "Done  " << "(" << counter.elapsed() << "s)." << std::endl;
+
+    // -----------------------------------------------------
+    FOctree<FmbParticule, FmbCell, FSimpleLeaf, NbLevels, SizeSubLevels> tree(BoxWidth,CenterOfBox);
+    FOctree<FmbParticuleTyped, FmbCellTyped, FTypedLeaf, NbLevels, SizeSubLevels> treeTyped(BoxWidth,CenterOfBox);
+
+
+    std::cout << "Inserting particules ..." << std::endl;
+    counter.tic();
+    for(long idxPart = 0 ; idxPart < NbPart ; ++idxPart){
+        tree.insert(&particules[idxPart]);
+        treeTyped.insert(&particulesTyped[idxPart*2]);
+        treeTyped.insert(&particulesTyped[idxPart*2+1]);
+    }
+    counter.tac();
+    std::cout << "Done  " << "(" << counter.elapsed() << "s)." << std::endl;
+
+    // -----------------------------------------------------
+
+    std::cout << "Working on particules ..." << std::endl;
+    counter.tic();
+
+    //FFmbKernelsPotentialForces FFmbKernelsForces FFmbKernelsPotential
+    FFmbKernelsPotentialForces<FmbParticule, FmbCell, NbLevels> kernels(BoxWidth);
+    FFmbKernelsPotentialForces<FmbParticuleTyped, FmbCellTyped, NbLevels> kernelsTyped(BoxWidth);
+    //FFmmAlgorithm FFmmAlgorithmArray
+    FFmmAlgorithmArray<FFmbKernelsPotentialForces, FmbParticule, FmbCell, FSimpleLeaf, NbLevels, SizeSubLevels> algo(&tree,&kernels);
+    FFmmAlgorithmArrayTsm<FFmbKernelsPotentialForces, FmbParticuleTyped, FmbCellTyped, FTypedLeaf, NbLevels, SizeSubLevels> algoTyped(&treeTyped,&kernelsTyped);
+    algo.execute();
+    algoTyped.execute();
+
+    counter.tac();
+    std::cout << "Done  " << "(" << counter.elapsed() << "s)." << std::endl;
+
+    // -----------------------------------------------------
+
+    // Here we compare the cells of the trees that must contains the same values
+
+    std::cout << "Start checking ..." << std::endl;
+    {
+        FOctree<FmbParticule, FmbCell, FSimpleLeaf, NbLevels, SizeSubLevels>::Iterator octreeIterator(&tree);
+        octreeIterator.gotoBottomLeft();
+
+        FOctree<FmbParticuleTyped, FmbCellTyped, FTypedLeaf, NbLevels, SizeSubLevels>::Iterator octreeIteratorTyped(&treeTyped);
+        octreeIteratorTyped.gotoBottomLeft();
+
+        for(int idxLevel = NbLevels - 1 ; idxLevel > 1 ; --idxLevel ){
+            std::cout << "\t test level " << idxLevel << "\n";
+
+            do{
+                bool poleDiff = false;
+                bool localDiff = false;
+                for(int idxValues = 0 ; idxValues < FExtendFmbCell::MultipoleSize && !(poleDiff && localDiff); ++idxValues){
+                    const FComplexe pole = octreeIterator.getCurrentCell()->getMultipole()[idxValues];
+                    const FComplexe poleTyped = octreeIteratorTyped.getCurrentCell()->getMultipole()[idxValues];
+                    if(!FMath::LookEqual(pole.getImag(),poleTyped.getImag()) || !FMath::LookEqual(pole.getReal(),poleTyped.getReal())){
+                        poleDiff = true;
+                        printf("Pole diff imag( %.15e , %.15e ) real( %.15e , %.15e)\n",
+                               pole.getImag(),poleTyped.getImag(),pole.getReal(),poleTyped.getReal());
+                    }
+                    const FComplexe local = octreeIterator.getCurrentCell()->getLocal()[idxValues];
+                    const FComplexe localTyped = octreeIteratorTyped.getCurrentCell()->getLocal()[idxValues];
+                    if(!FMath::LookEqual(local.getImag(),localTyped.getImag()) || !FMath::LookEqual(local.getReal(),localTyped.getReal())){
+                        localDiff = true;
+                        printf("Pole diff imag( %.15e , %.15e ) real( %.15e , %.15e)\n",
+                               local.getImag(),localTyped.getImag(),local.getReal(),localTyped.getReal());
+                    }
+                }
+                if(poleDiff){
+                    std::cout << "Multipole error at level " << idxLevel << "\n";
+                }
+                if(localDiff){
+                    std::cout << "Locale error at level " << idxLevel << "\n";
+                }
+            } while(octreeIterator.moveRight() && octreeIteratorTyped.moveRight());
+
+            octreeIterator.moveUp();
+            octreeIterator.gotoLeft();
+
+            octreeIteratorTyped.moveUp();
+            octreeIteratorTyped.gotoLeft();
+        }
+    }
+
+    std::cout << "Done ..." << std::endl;
+
+    delete [] particules;
+    delete [] particulesTyped;
+
+    return 0;
+}
+
+
+// [--LICENSE--]