Clean the new kernel, add a dynamic-mem cell class,

try to test the new kenerl.

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

Src/Core/FFmmAlgorithm.hpp

@@ -222,7 +222,7 @@ private:
         // for each leafs
         do{
             FDEBUG(computationCounterL2P.tic());
-            ///kernels->L2P(octreeIterator.getCurrentCell(), octreeIterator.getCurrentListTargets());
+            kernels->L2P(octreeIterator.getCurrentCell(), octreeIterator.getCurrentListTargets());
             FDEBUG(computationCounterL2P.tac());
             // need the current particles and neighbors particles
             const int counter = tree->getLeafsNeighborsWithIndex(neighbors, neighborsIndex, octreeIterator.getCurrentGlobalIndex(),heightMinusOne);

Src/Kernels/FComputeCell.hpp

+#ifndef FCOMPUTECELL_HPP
+#define FCOMPUTECELL_HPP
+// [--License--]
+
+
+#include "../Utils/FComplexe.hpp"
+#include "../Utils/FMemUtils.hpp"
+
+#include "../Components/FBasicCell.hpp"
+
+/**
+* @author Berenger Bramas (berenger.bramas@inria.fr)
+* @class FComputeCell
+* Please read the license.
+*
+*/
+class FComputeCell : public FBasicCell {
+protected:
+    static int DevP;
+    static int ExpP;
+
+    FComplexe* multipole_exp; //< For multipole extenssion
+    FComplexe* local_exp;     //< For local extenssion
+
+public:
+    static void Init(const int inDevP){
+        DevP = inDevP;
+        ExpP = int((DevP+1) * (DevP+2) * 0.5);
+    }
+
+
+    /** Default constructor */
+    FComputeCell()
+        : multipole_exp(0), local_exp(0){
+        multipole_exp = new FComplexe[ExpP];
+        local_exp = new FComplexe[ExpP];
+    }
+
+    /** Constructor */
+    FComputeCell(const FComputeCell& other)
+        : multipole_exp(0), local_exp(0){
+        multipole_exp = new FComplexe[ExpP];
+        local_exp = new FComplexe[ExpP];
+        (*this) = other;
+    }
+
+    /** Default destructor */
+    virtual ~FComputeCell(){
+        delete[] multipole_exp;
+        delete[] local_exp;
+    }
+
+    /** Copy constructor */
+    FComputeCell& operator=(const FComputeCell& other) {
+        FMemUtils::copyall(multipole_exp, other.multipole_exp, ExpP);
+        FMemUtils::copyall(local_exp, other.local_exp, ExpP);
+        return *this;
+    }
+
+    /** Get Multipole */
+    const FComplexe* getMultipole() const {
+        return multipole_exp;
+    }
+    /** Get Local */
+    const FComplexe* getLocal() const {
+        return local_exp;
+    }
+
+    /** Get Multipole */
+    FComplexe* getMultipole() {
+        return multipole_exp;
+    }
+    /** Get Local */
+    FComplexe* getLocal() {
+        return local_exp;
+    }
+};
+
+
+int FComputeCell::DevP(-1);
+int FComputeCell::ExpP(-1);
+
+
+#endif //FCOMPUTECELL_HPP
+
+// [--END--]

Src/Kernels/FElecForcesKernels.hpp

@@ -517,7 +517,7 @@ private:
 
                     for(/*l = -n-1 or l = -k-1 */; l >= -n ; --l){ // we have -k-l>=0 and l<=0
                         const FComplexe M_n_l = multipole_exp_src[index_n - l];
-                        const FComplexe O_n_j__k_l = M2L_Outer_transfer[index_n_j + k - l + k];
+                        const FComplexe O_n_j__k_l = M2L_Outer_transfer[index_n_j + (k - l)];
 
                         L_j_k.incReal( pow_of_minus_1_for_l *
                                                     ((M_n_l.getReal() * O_n_j__k_l.getReal()) +
@@ -638,9 +638,9 @@ private:
         const FSpherical spherical(particle->getPosition() - local_position);
         harmonic.computeInnerTheta( spherical );
 
-        int index_j_k = 1;
+        int index_j_k = 0;
 
-        for (int j = 1 ; j <= devP ; ++j ){
+        for (int j = 0 ; j <= devP ; ++j ){
             {
                 // k=0:
                 // F_r:
@@ -754,7 +754,6 @@ private:
 
         //--------------- Potential ----------------//
 
-        // TODO !! check here!!!!
         FReal result = 0.0;
         index_j_k = 0;
 
@@ -771,11 +770,11 @@ private:
             }
         }
 
-        particle->incPotential(result);
+        particle->incPotential(result * physicalValue);
 
     }
 
-
+public:
     /** P2P mutual interaction
       * F = q * q' / r²
       */
@@ -796,7 +795,7 @@ private:
         dz *= inv_square_distance;
 
         target->incForces( dx, dy, dz);
-        target->incPotential( inv_distance  * source->getPhysicalValue() );
+        target->incPotential( inv_distance * source->getPhysicalValue() );
 
         source->incForces( (-dx), (-dy), (-dz));
         source->incPotential( inv_distance * target->getPhysicalValue() );

Src/Kernels/FHarmonic.hpp

@@ -36,6 +36,7 @@ public:
 
 
 class FHarmonic {
+public: //TODO delete
     const int devP;     //< P
     const int expSize;  //< Exponen Size
 
@@ -61,7 +62,7 @@ class FHarmonic {
         FReal factOuter = 1.0;
         for(int idxP = 0 ; idxP <= devP; ++idxP ){
             sphereHarmoOuterCoef[idxP] = factOuter;
-            factOuter *= FReal(idxP);
+            factOuter *= FReal(idxP+1);
         }
 
         // Pre compute coef
@@ -73,11 +74,12 @@ class FHarmonic {
             FReal fact_l_m = factInner;
             for(int m = 0 ; m <= l ; ++m){
                 sphereHarmoInnerCoef[index_l_m] = minus_1_pow_l / fact_l_m;
-                fact_l_m *= FReal(l + m);
+
+                fact_l_m *= FReal(l + m + 1);
                 ++index_l_m;
             }
             minus_1_pow_l = -minus_1_pow_l;
-            factInner *= FReal(l);
+            factInner *= FReal(l+1);
         }
 
         // Smart redirection
@@ -248,7 +250,7 @@ public:
 
         for(int l = 0 ; l <= devP ; ++l){
             for(int m = 0 ; m <= l ; ++m){
-                const FReal magnitude = this->sphereHarmoOuterCoef[l-m] * invR_pow_l1 * legendre[index_l_m];
+                const FReal magnitude = sphereHarmoOuterCoef[l-m] * invR_pow_l1 * legendre[index_l_m];
                 harmonic[index_l_m].setReal( magnitude * cosSin[m].getReal() );
                 harmonic[index_l_m].setImag( magnitude * cosSin[m].getImag() );
 

Src/Utils/FComplexe.hpp

@@ -137,6 +137,11 @@ public:
         this->imag = (tempReal * other.imag) + (this->imag * other.real);
         return *this;
     }
+
+    /** Test if a complex is not a number */
+    bool isNan() const {
+        return FMath::IsNan(imag) || FMath::IsNan(real);
+    }
 };
 
 /** Global operator Mul a complexe by another "c=c1*c2" */

Src/Utils/FMath.hpp

@@ -2,7 +2,9 @@
 #define FMATH_HPP
 // [--License--]
 
-#include <math.h>
+#include <cmath>
+#include <climits>
+
 #include "FGlobal.hpp"
 
 /**
@@ -121,6 +123,21 @@ struct FMath{
     static double Fmod(const double inValue1,const double inValue2){
         return fmod(inValue1,inValue2);
     }
+
+    /** To know if a variable is nan, based on the C++0x */
+    template <class TestClass>
+    static bool IsNan(const TestClass& value){
+        //volatile const TestClass* const pvalue = &value;
+        //return (*pvalue) != value;
+        return std::isnan(value);
+    }
+
+    /** To know if a variable is not inf, based on the C++0x */
+    template <class TestClass>
+    static bool IsFinite(const TestClass& value){
+        // return !(value <= std::numeric_limits<T>::min()) && !(std::numeric_limits<T>::max() <= value);
+        return std::isfinite(value);
+    }
 };
 
 

Tests/testFmbAlgorithm.cpp

@@ -20,8 +20,6 @@
 #include "../Src/Files/FFmaLoader.hpp"
 
 
-#include "../Src/Kernels/FElecForcesKernels.hpp"
-
 // With openmp : g++ testFmbAlgorithm.cpp ../Src/Utils/FDebug.cpp ../Src/Utils/FTrace.cpp -lgomp -fopenmp -O2 -o testFmbAlgorithm.exe
 // icpc -openmp -openmp-lib=compat testFmbAlgorithm.cpp ../Src/Utils/FAssertable.cpp ../Src/Utils/FDebug.cpp -O2 -o testFmbAlgorithm.exe
 
@@ -40,8 +38,7 @@ int main(int argc, char ** argv){
 
     typedef FSimpleLeaf<ParticleClass, ContainerClass >                     LeafClass;
     typedef FOctree<ParticleClass, CellClass, ContainerClass , LeafClass >  OctreeClass;
-    //typedef FFmbKernels<ParticleClass, CellClass, ContainerClass >          KernelClass;
-    typedef FElecForcesKernels<ParticleClass, CellClass, ContainerClass >          KernelClass;
+    typedef FFmbKernels<ParticleClass, CellClass, ContainerClass >          KernelClass;
 
     typedef FFmmAlgorithm<OctreeClass, ParticleClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;
     ///////////////////////What we do/////////////////////////////
@@ -96,8 +93,7 @@ int main(int argc, char ** argv){
     std::cout << "Working on particles ..." << std::endl;
     counter.tic();
 
-    //KernelClass kernels(NbLevels,loader.getBoxWidth());
-    KernelClass kernels(6,NbLevels,loader.getBoxWidth());
+    KernelClass kernels(NbLevels,loader.getBoxWidth());
     FmmClass algo(&tree,&kernels);
     algo.execute();
 

UTests/utestFmbDirect.cpp

@@ -7,6 +7,10 @@
 #include "../Src/Fmb/FFmbKernels.hpp"
 #include "../Src/Fmb/FFmbComponents.hpp"
 
+#include "../Src/Kernels/FComputeCell.hpp"
+
+#include "../Src/Kernels/FElecForcesKernels.hpp"
+
 #include "../Src/Files/FFmaBinLoader.hpp"
 #include "../Src/Files/FTreeIO.hpp"
 
@@ -32,10 +36,13 @@ public:
 
 class TestFmbDirect : public FUTester<TestFmbDirect> {
     typedef IndexedParticle         ParticleClass;
-    typedef FmbCell                 CellClass;
+    //typedef FmbCell            CellClass;
+    typedef FComputeCell            CellClass;
     typedef FVector<ParticleClass>  ContainerClass;
 
-    typedef FFmbKernels<ParticleClass, CellClass, ContainerClass >          KernelClass;
+    //typedef FFmbKernels<ParticleClass, CellClass, ContainerClass >          KernelClass;
+    typedef FElecForcesKernels<ParticleClass, CellClass, ContainerClass >          KernelClass;
+
     typedef FSimpleLeaf<ParticleClass, ContainerClass >                     LeafClass;
     typedef FOctree<ParticleClass, CellClass, ContainerClass , LeafClass >  OctreeClass;
 
@@ -43,9 +50,6 @@ class TestFmbDirect : public FUTester<TestFmbDirect> {
 
 
     void TestDirect(){
-        const int NbLevels      = 4;
-        const int SizeSubLevels = 2;
-
         FFmaBinLoader<ParticleClass> loader("../Data/utestFmbDirect.bin.fma");
         if(!loader.isOpen()){
             Print("Cannot open particles file.");
@@ -54,6 +58,10 @@ class TestFmbDirect : public FUTester<TestFmbDirect> {
         Print("Number of particles:");
         Print(loader.getNumberOfParticles());
 
+        const int NbLevels      = 4;
+        const int SizeSubLevels = 2;
+        const int DevP = 5;
+        FComputeCell::Init(DevP);
         OctreeClass tree(NbLevels, SizeSubLevels, loader.getBoxWidth(), loader.getCenterOfBox());
 
 
@@ -66,7 +74,8 @@ class TestFmbDirect : public FUTester<TestFmbDirect> {
 
 
         Print("Fmm...");
-        KernelClass kernels(NbLevels,loader.getBoxWidth());
+        //KernelClass kernels(NbLevels,loader.getBoxWidth());
+        KernelClass kernels(DevP,NbLevels,loader.getBoxWidth());
         FmmClass algo(&tree,&kernels);
         algo.execute();
 
@@ -74,7 +83,8 @@ class TestFmbDirect : public FUTester<TestFmbDirect> {
         Print("Direct...");
         for(int idxTarget = 0 ; idxTarget < loader.getNumberOfParticles() ; ++idxTarget){
             for(int idxOther = idxTarget + 1 ; idxOther < loader.getNumberOfParticles() ; ++idxOther){
-                kernels.DIRECT_COMPUTATION_MUTUAL_SOFT(particles[idxTarget], particles[idxOther]);
+                //kernels.DIRECT_COMPUTATION_MUTUAL_SOFT(particles[idxTarget], particles[idxOther]);
+                kernels.directInteractionMutual(&particles[idxTarget], &particles[idxOther]);
             }
         }
 
@@ -90,22 +100,25 @@ class TestFmbDirect : public FUTester<TestFmbDirect> {
                 typename ContainerClass::BasicIterator leafIter(*octreeIterator.getCurrentListTargets());
 
                 while( leafIter.hasNotFinished() ){
-                    const FReal currentPotentialDiff = FMath::Abs(particles[leafIter.data().getIndex()].getPotential() - leafIter.data().getPotential());
+                    const ParticleClass& currentParticle = leafIter.data();//TODO delete
+                    const ParticleClass& other = particles[leafIter.data().getIndex()];//TODO delete
+
+                    const FReal currentPotentialDiff = FMath::Abs(other.getPotential() - leafIter.data().getPotential());
                     if( potentialDiff < currentPotentialDiff ){
                         potentialDiff = currentPotentialDiff;
                     }
 
-                    const FReal currentFx = FMath::Abs(particles[leafIter.data().getIndex()].getPosition().getX() - leafIter.data().getPosition().getX());
+                    const FReal currentFx = FMath::Abs(other.getForces().getX() - leafIter.data().getForces().getX());
                     if( fx < currentFx ){
                         fx = currentFx;
                     }
 
-                    const FReal currentFy = FMath::Abs(particles[leafIter.data().getIndex()].getPosition().getY() - leafIter.data().getPosition().getY());
+                    const FReal currentFy = FMath::Abs(other.getForces().getY() - leafIter.data().getForces().getY());
                     if( fy < currentFy ){
                         fy = currentFy;
                     }
 
-                    const FReal currentFz = FMath::Abs(particles[leafIter.data().getIndex()].getPosition().getZ() - leafIter.data().getPosition().getZ());
+                    const FReal currentFz = FMath::Abs(other.getForces().getZ() - leafIter.data().getForces().getZ());
                     if( fz < currentFz ){
                         fz = currentFz;
                     }