Merge branch 'master' of git+ssh://scm.gforge.inria.fr//gitroot/scalfmm/scalfmm

CMakeModules/morse/find/FindPTSCOTCH.cmake

@@ -264,10 +264,21 @@ if(PTSCOTCH_LIBRARIES)
     if(CMAKE_THREAD_LIBS_INIT)
         list(APPEND REQUIRED_LIBS "${CMAKE_THREAD_LIBS_INIT}")
     endif()
-    if(UNIX OR WIN32)
+    set(Z_LIBRARY "Z_LIBRARY-NOTFOUND")
+    find_library(Z_LIBRARY NAMES z)
+    if(Z_LIBRARY)
+        list(APPEND REQUIRED_LIBS "-lz")
+    endif()
+    set(M_LIBRARY "M_LIBRARY-NOTFOUND")
+    find_library(M_LIBRARY NAMES m)
+    if(M_LIBRARY)
         list(APPEND REQUIRED_LIBS "-lm")
     endif()
-    list(APPEND REQUIRED_LIBS "-lz -lrt")
+    set(RT_LIBRARY "RT_LIBRARY-NOTFOUND")
+    find_library(RT_LIBRARY NAMES rt)
+    if(RT_LIBRARY)
+        list(APPEND REQUIRED_LIBS "-lrt")
+    endif()
 
     # set required libraries for link
     set(CMAKE_REQUIRED_INCLUDES "${REQUIRED_INCDIRS}")

CMakeModules/morse/find/FindSCOTCH.cmake

@@ -233,10 +233,21 @@ if(SCOTCH_LIBRARIES)
     if(CMAKE_THREAD_LIBS_INIT)
         list(APPEND REQUIRED_LIBS "${CMAKE_THREAD_LIBS_INIT}")
     endif()
-    if(UNIX OR WIN32)
+    set(Z_LIBRARY "Z_LIBRARY-NOTFOUND")
+    find_library(Z_LIBRARY NAMES z)
+    if(Z_LIBRARY)
+        list(APPEND REQUIRED_LIBS "-lz")
+    endif()
+    set(M_LIBRARY "M_LIBRARY-NOTFOUND")
+    find_library(M_LIBRARY NAMES m)
+    if(M_LIBRARY)
         list(APPEND REQUIRED_LIBS "-lm")
     endif()
-    list(APPEND REQUIRED_LIBS "-lz -lrt")
+    set(RT_LIBRARY "RT_LIBRARY-NOTFOUND")
+    find_library(RT_LIBRARY NAMES rt)
+    if(RT_LIBRARY)
+        list(APPEND REQUIRED_LIBS "-lrt")
+    endif()
 
     # set required libraries for link
     set(CMAKE_REQUIRED_INCLUDES "${REQUIRED_INCDIRS}")

Src/Core/FFmmAlgorithmSectionTask.hpp

@@ -27,6 +27,7 @@
 #include "../Containers/FVector.hpp"
 
 #include "FCoreCommon.hpp"
+#include "FP2PExclusion.hpp"
 
 /**
  * @author Berenger Bramas (berenger.bramas@inria.fr)
@@ -45,7 +46,7 @@
  *
  * Upon destruction, this class does not deallocate pointers given to its constructor.
  */
-template<class OctreeClass, class CellClass, class ContainerClass, class KernelClass, class LeafClass>
+template<class OctreeClass, class CellClass, class ContainerClass, class KernelClass, class LeafClass, class P2PExclusionClass = FP2PMiddleExclusion>
 class FFmmAlgorithmSectionTask : public FAbstractAlgorithm, public FAlgorithmTimers {
     
     OctreeClass* const tree;  ///< The octree to work on
@@ -328,7 +329,7 @@ protected:
         // There is a maximum of 26 neighbors
         ContainerClass* neighbors[27];
 
-        const int SizeShape = 3*3*3;
+        const int SizeShape = P2PExclusionClass::SizeShape;
         FVector<typename OctreeClass::Iterator> shapes[SizeShape];
 
         typename OctreeClass::Iterator octreeIterator(tree);
@@ -338,7 +339,7 @@ protected:
         // Coloring all the cells
         do{
             const FTreeCoordinate& coord = octreeIterator.getCurrentGlobalCoordinate();
-            const int shapePosition = (coord.getX()%3)*9 + (coord.getY()%3)*3 + (coord.getZ()%3);
+            const int shapePosition = P2PExclusionClass::GetShapeIdx(coord);
 
             shapes[shapePosition].push(octreeIterator);
 

Src/Core/FFmmAlgorithmTask.hpp

@@ -27,6 +27,7 @@
 #include "../Containers/FVector.hpp"
 
 #include "FCoreCommon.hpp"
+#include "FP2PExclusion.hpp"
 
 /**
  * @author Berenger Bramas (berenger.bramas@inria.fr)
@@ -39,7 +40,7 @@
  *
  * Of course this class does not deallocate pointer given in arguements.
  */
-template<class OctreeClass, class CellClass, class ContainerClass, class KernelClass, class LeafClass>
+template<class OctreeClass, class CellClass, class ContainerClass, class KernelClass, class LeafClass, class P2PExclusionClass = FP2PMiddleExclusion>
 class FFmmAlgorithmTask : public FAbstractAlgorithm, public FAlgorithmTimers {
 
 	OctreeClass* const tree;       //< The octree to work on
@@ -389,7 +390,7 @@ protected:
 				// There is a maximum of 26 neighbors
 				ContainerClass* neighbors[27];
 
-				const int SizeShape = 3*3*3;
+                const int SizeShape = P2PExclusionClass::SizeShape;
 				FVector<typename OctreeClass::Iterator> shapes[SizeShape];
 
 				typename OctreeClass::Iterator octreeIterator(tree);
@@ -398,7 +399,7 @@ protected:
 				// for each leafs
 				do{
 					const FTreeCoordinate& coord = octreeIterator.getCurrentGlobalCoordinate();
-					const int shapePosition = (coord.getX()%3)*9 + (coord.getY()%3)*3 + (coord.getZ()%3);
+                    const int shapePosition = P2PExclusionClass::GetShapeIdx(coord);
 
 					shapes[shapePosition].push(octreeIterator);
 

Src/Core/FFmmAlgorithmThread.hpp

@@ -27,6 +27,7 @@
 #include "../Containers/FOctree.hpp"
 
 #include "FCoreCommon.hpp"
+#include "FP2PExclusion.hpp"
 
 #include <omp.h>
 
@@ -45,7 +46,7 @@
 *
 * This class does not deallocate pointers given to its constructor.
 */
-template<class OctreeClass, class CellClass, class ContainerClass, class KernelClass, class LeafClass>
+template<class OctreeClass, class CellClass, class ContainerClass, class KernelClass, class LeafClass, class P2PExclusionClass = FP2PMiddleExclusion>
 class FFmmAlgorithmThread : public FAbstractAlgorithm, public FAlgorithmTimers{
     OctreeClass* const tree;                  ///< The octree to work on.
     KernelClass** kernels;                    ///< The kernels.
@@ -53,7 +54,7 @@ class FFmmAlgorithmThread : public FAbstractAlgorithm, public FAlgorithmTimers{
     typename OctreeClass::Iterator* iterArray;
     int leafsNumber;
 
-    static const int SizeShape = 3*3*3;
+    static const int SizeShape = P2PExclusionClass::SizeShape;
     int shapeLeaf[SizeShape];
 
     const int MaxThreads;                     ///< The maximum number of threads.
@@ -140,7 +141,7 @@ protected:
         do{
             ++leafsNumber;
             const FTreeCoordinate& coord = octreeIterator.getCurrentCell()->getCoordinate();
-            ++this->shapeLeaf[(coord.getX()%3)*9 + (coord.getY()%3)*3 + (coord.getZ()%3)];
+            ++this->shapeLeaf[P2PExclusionClass::GetShapeIdx(coord)];
 
         } while(octreeIterator.moveRight());
         iterArray = new typename OctreeClass::Iterator[leafsNumber];
@@ -441,7 +442,7 @@ protected:
                 //iterArray[leafs] = octreeIterator;
                 //++leafs;
                 const FTreeCoordinate& coord = octreeIterator.getCurrentGlobalCoordinate();
-                const int shapePosition = (coord.getX()%3)*9 + (coord.getY()%3)*3 + (coord.getZ()%3);
+                const int shapePosition = P2PExclusionClass::GetShapeIdx(coord);
 
                 omp_set_lock(&lockShape[shapePosition]);
                 const int positionToWork = startPosAtShape[shapePosition]++;

Src/Core/FFmmAlgorithmThreadBalance.hpp

@@ -11,6 +11,7 @@
 #include "../Containers/FOctree.hpp"
 
 #include "FCoreCommon.hpp"
+#include "FP2PExclusion.hpp"
 
 #include <omp.h>
 #include <vector>
@@ -29,12 +30,12 @@
 *
 * This class does not deallocate pointers given to its constructor.
 */
-template<class OctreeClass, class CellClass, class ContainerClass, class KernelClass, class LeafClass>
+template<class OctreeClass, class CellClass, class ContainerClass, class KernelClass, class LeafClass, class P2PExclusionClass = FP2PMiddleExclusion>
 class FFmmAlgorithmThreadBalance : public FAbstractAlgorithm, public FAlgorithmTimers{
     OctreeClass* const tree;                  ///< The octree to work on.
     KernelClass** kernels;                    ///< The kernels.
 
-    static const int SizeShape = 3*3*3;
+    static const int SizeShape = P2PExclusionClass::SizeShape;
 
     const int MaxThreads;                     ///< The maximum number of threads.
 
@@ -206,7 +207,7 @@ protected:
             do{
                 ++leafsNumber;
                 const FTreeCoordinate& coord = octreeIterator.getCurrentCell()->getCoordinate();
-                ++shapeLeaves[(coord.getX()%3)*9 + (coord.getY()%3)*3 + (coord.getZ()%3)];
+                ++shapeLeaves[P2PExclusionClass::GetShapeIdx(coord)];
             } while(octreeIterator.moveRight());
         }
 
@@ -367,7 +368,7 @@ protected:
                     // for each leafs
                     for(int idxLeaf = 0 ; idxLeaf < leafsNumber ; ++idxLeaf){
                         const FTreeCoordinate& coord = octreeIterator.getCurrentGlobalCoordinate();
-                        const int shapePosition = (coord.getX()%3)*9 + (coord.getY()%3)*3 + (coord.getZ()%3);
+                        const int shapePosition = P2PExclusionClass::GetShapeIdx(coord);
 
                         const int positionToWork = startPosAtShape[shapePosition]++;
 

Src/Core/FFmmAlgorithmThreadProc.hpp

@@ -40,6 +40,7 @@
 #include <sys/time.h>
 
 #include "FCoreCommon.hpp"
+#include "FP2PExclusion.hpp"
 
 #include <memory>
 
@@ -63,7 +64,7 @@
  *        --tool=memcheck --leak-check=yes --show-reachable=yes --num-callers=20
  *        --track-fds=yes ./Tests/testFmmAlgorithmProc ../Data/testLoaderSmall.fma.tmp
  */
-template<class OctreeClass, class CellClass, class ContainerClass, class KernelClass, class LeafClass>
+template<class OctreeClass, class CellClass, class ContainerClass, class KernelClass, class LeafClass, class P2PExclusionClass = FP2PMiddleExclusion>
 class FFmmAlgorithmThreadProc : public FAbstractAlgorithm, public FAlgorithmTimers {
 private:
     OctreeClass* const tree;     ///< The octree to work on
@@ -1233,7 +1234,7 @@ protected:
 
         // init
         const int LeafIndex = OctreeHeight - 1;
-        const int SizeShape = 3*3*3;
+        const int SizeShape = P2PExclusionClass::SizeShape;
 
         int shapeLeaf[SizeShape];
         memset(shapeLeaf,0,SizeShape*sizeof(int));
@@ -1394,7 +1395,7 @@ protected:
                     myLeafs[idxLeaf] = octreeIterator;
 
                     const FTreeCoordinate& coord = octreeIterator.getCurrentCell()->getCoordinate();
-                    const int shape = (coord.getX()%3)*9 + (coord.getY()%3)*3 + (coord.getZ()%3);
+                    const int shape = P2PExclusionClass::GetShapeIdx(coord);
                     shapeType[idxLeaf] = shape;
 
                     ++shapeLeaf[shape];

Src/Core/FFmmAlgorithmThreadProcPeriodic.hpp

@@ -38,6 +38,7 @@
 #include <omp.h>
 
 #include "FCoreCommon.hpp"
+#include "FP2PExclusion.hpp"
 
 #include <memory>
 
@@ -61,7 +62,7 @@
  * --tool=memcheck --leak-check=yes --show-reachable=yes --num-callers=20 --track-fds=yes
  * ./Tests/testFmmAlgorithmProc ../Data/testLoaderSmall.fma.tmp
  */
-template<class FReal, class OctreeClass, class CellClass, class ContainerClass, class KernelClass, class LeafClass>
+template<class FReal, class OctreeClass, class CellClass, class ContainerClass, class KernelClass, class LeafClass, class P2PExclusionClass = FP2PMiddleExclusion>
 class FFmmAlgorithmThreadProcPeriodic : public FAbstractAlgorithm {
     OctreeClass* const tree;                 //< The octree to work on
     KernelClass** kernels;                   //< The kernels
@@ -1348,7 +1349,7 @@ protected:
 
         // init
         const int LeafIndex = OctreeHeight - 1;
-        const int SizeShape = 3*3*3;
+        const int SizeShape = P2PExclusionClass::SizeShape;
 
         int shapeLeaf[SizeShape];
         memset(shapeLeaf,0,SizeShape*sizeof(int));
@@ -1511,7 +1512,7 @@ protected:
                     myLeafs[idxLeaf] = octreeIterator;
 
                     const FTreeCoordinate& coord = octreeIterator.getCurrentCell()->getCoordinate();
-                    const int shape = (coord.getX()%3)*9 + (coord.getY()%3)*3 + (coord.getZ()%3);
+                    const int shape = P2PExclusionClass::GetShapeIdx(coord);
                     shapeType[idxLeaf] = shape;
 
                     ++shapeLeaf[shape];

Src/Core/FP2PExclusion.hpp

+#ifndef FP2PEXCLUSION_HPP
+#define FP2PEXCLUSION_HPP
+
+#include "../Containers/FTreeCoordinate.hpp"
+
+/**
+ * This class gives is responsible of the separation of the leaves
+ * using the coloring algorithm.
+ * In case of classic P2P and mutual interaction the BoxSeparations = 2 should be used.
+ * For our current mutual P2P is is a little more complicated because we need
+ * 2 boxes of separation but only in some directions.
+ */
+template <int BoxSeparations = 2>
+class FP2PExclusion{
+public:
+    static const int BoxesPerDim = (BoxSeparations+1);
+    static const int SizeShape = BoxesPerDim*BoxesPerDim*BoxesPerDim;
+
+    static int GetShapeIdx(const int inX, const int inY, const int inZ){
+        return (inX%BoxesPerDim)*(BoxesPerDim*BoxesPerDim) + (inY%BoxesPerDim)*BoxesPerDim + (inZ%BoxesPerDim);
+    }
+
+    static int GetShapeIdx(const FTreeCoordinate& coord){
+        return GetShapeIdx(coord.getX(), coord.getY(), coord.getZ());
+    }
+};
+
+/**
+ * Here the formula is related to the octree construction of neighbors list:
+ * const int index = (((idxX + 1) * 3) + (idxY +1)) * 3 + idxZ + 1;
+ * If go from 0 to 27,
+ * if we loop from 0 to 14, then we need "x" in [0;2[
+ * "y" "z" in [0;3[
+ */
+class FP2PMiddleExclusion{
+public:
+    static const int SizeShape = 3*3*2;
+
+    static int GetShapeIdx(const int inX, const int inY, const int inZ){
+        return (inX%2)*9 + (inY%3)*3 + (inZ%3);
+    }
+
+    static int GetShapeIdx(const FTreeCoordinate& coord){
+        return GetShapeIdx(coord.getX(), coord.getY(), coord.getZ());
+    }
+};
+
+
+#endif // FP2PEXCLUSION_HPP
+

UTests/utestP2PExclusion.cpp

+
+// ===================================================================================
+// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Bérenger Bramas, Matthias Messner
+// 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".
+// ===================================================================================
+#include "FUTester.hpp"
+
+#include "Core/FP2PExclusion.hpp"
+#include "Utils/FMath.hpp"
+
+#include <memory>
+
+/**
+* This file is a unit test for the FNeigborIndexes classes
+*/
+
+
+/** this class test the list container */
+class TestExclusion : public FUTester<TestExclusion> {
+    const int Size = 100;
+
+    void Exclusion2(){
+        const int Width = 2;
+        std::unique_ptr<int[]> grid(new int[Size*Size*Size]);
+        for(int idxShape = 0 ; idxShape < FP2PExclusion<Width>::SizeShape ; ++idxShape){
+            memset(grid.get(), 0, sizeof(int)*Size*Size*Size);
+
+            for(int idxX = 0 ; idxX < Size ; ++idxX){
+                for(int idxY = 0 ; idxY < Size ; ++idxY){
+                    for(int idxZ = 0 ; idxZ < Size ; ++idxZ){
+                        if(FP2PExclusion<Width>::GetShapeIdx(idxX,idxY,idxZ) == idxShape){
+                            for(int idxX_neig = FMath::Max(0,idxX-1) ; idxX_neig < FMath::Min(Size,idxX+1) ; ++idxX_neig){
+                                for(int idxY_neig = FMath::Max(0,idxY-1) ; idxY_neig < FMath::Min(Size,idxY+1) ; ++idxY_neig){
+                                    for(int idxZ_neig = FMath::Max(0,idxZ-1) ; idxZ_neig < FMath::Min(Size,idxZ+1) ; ++idxZ_neig){
+                                        uassert(grid[(idxX_neig*Size + idxY_neig)*Size + idxZ_neig] == 0);
+                                        grid[grid[(idxX_neig*Size + idxY_neig)*Size + idxZ_neig]] = 1;
+                                    }
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    void Exclusion1(){
+        const int Width = 1;
+        std::unique_ptr<int[]> grid(new int[Size*Size*Size]);
+        for(int idxShape = 0 ; idxShape < FP2PExclusion<Width>::SizeShape ; ++idxShape){
+            memset(grid.get(), 0, sizeof(int)*Size*Size*Size);
+
+            for(int idxX = 0 ; idxX < Size ; ++idxX){
+                for(int idxY = 0 ; idxY < Size ; ++idxY){
+                    for(int idxZ = 0 ; idxZ < Size ; ++idxZ){
+                        if(FP2PExclusion<Width>::GetShapeIdx(idxX,idxY,idxZ) == idxShape){
+                            for(int idxX_neig = FMath::Max(0,idxX-1) ; idxX_neig < idxX ; ++idxX_neig){
+                                for(int idxY_neig = FMath::Max(0,idxY-1) ; idxY_neig < idxY ; ++idxY_neig){
+                                    for(int idxZ_neig = FMath::Max(0,idxZ-1) ; idxZ_neig < idxZ ; ++idxZ_neig){
+                                        uassert(grid[(idxX_neig*Size + idxY_neig)*Size + idxZ_neig] == 0);
+                                        grid[grid[(idxX_neig*Size + idxY_neig)*Size + idxZ_neig]] = 1;
+                                    }
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    void Middle(){
+        std::unique_ptr<int[]> grid(new int[Size*Size*Size]);
+        for(int idxShape = 0 ; idxShape < FP2PMiddleExclusion::SizeShape ; ++idxShape){
+            memset(grid.get(), 0, sizeof(int)*Size*Size*Size);
+
+            for(int idxX = 0 ; idxX < Size ; ++idxX){
+                for(int idxY = 0 ; idxY < Size ; ++idxY){
+                    for(int idxZ = 0 ; idxZ < Size ; ++idxZ){
+                        if(FP2PMiddleExclusion::GetShapeIdx(idxX,idxY,idxZ) == idxShape){
+                            for(int idxX_neig = FMath::Max(0,idxX-1) ; idxX_neig < FMath::Min(Size,idxX+1) ; ++idxX_neig){
+                                for(int idxY_neig = FMath::Max(0,idxY-1) ; idxY_neig < FMath::Min(Size,idxY+1) ; ++idxY_neig){
+                                    for(int idxZ_neig = FMath::Max(0,idxZ-1) ; idxZ_neig < FMath::Min(Size,idxZ+1) ; ++idxZ_neig){
+                                        const int diffx = idxX_neig-idxX;
+                                        const int diffy = idxY_neig-idxY;
+                                        const int diffz = idxZ_neig-idxZ;
+                                        const int idx = (diffx+1)*9 + (diffy+1)*3 + (diffz+1);
+                                        if(idx < 14){
+                                            uassert(grid[(idxX_neig*Size + idxY_neig)*Size + idxZ_neig] == 0);
+                                            grid[grid[(idxX_neig*Size + idxY_neig)*Size + idxZ_neig]] = 1;
+                                        }
+                                    }
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+
+    // set test
+    void SetTests(){
+        AddTest(&TestExclusion::Exclusion2,"Test 2 exclustion");
+        AddTest(&TestExclusion::Exclusion1,"Test 1 exclustion");
+        AddTest(&TestExclusion::Middle,"Test middle exclustion");
+    }
+};
+
+// You must do this
+TestClass(TestExclusion)
+
+
+