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

Examples/generateDistributions.cpp

@@ -102,7 +102,8 @@ int main(int argc, char ** argv){
                          "             -filename name: generic name for files (without extension) and save data\n"
                          "                     with following format in name.fma or name.bfma in -bin is set\n"
                          "             -visufmt  vtk, vtp, cosmo or cvs format.",
-                         FParameterDefinitions::InputFile, FParameterDefinitions::NbParticles);
+                         FParameterDefinitions::InputFile,  FParameterDefinitions::OutputBinFormat, FParameterDefinitions::NbParticles);
+
 
 
 
@@ -164,6 +165,7 @@ int main(int argc, char ** argv){
 		BoxWith = FMath::Max(A,FMath::Max(B,C) );
 		FReal halfBW = BoxWith*0.5;
 		Centre.setPosition(halfBW,halfBW,halfBW);
+		std::cout << "Cuboid "<< A << ":"<< B<<":"<<C<<std::endl;
 	}
 	else if(FParameters::existParameter(argc, argv, "-unitSphere")){
 		unifRandonPointsOnUnitSphere(NbPoints, particles) ;

Src/Containers/FParForEachOctree.hpp

+#ifndef FPARFOREACHOCTREE_HPP
+#define FPARFOREACHOCTREE_HPP
+
+#include "FOctree.hpp"
+
+#include <functional>
+
+namespace FParForEachOctree {
+
+/////////////////////////////////////////////////////////
+// Lambda function to apply to all member
+/////////////////////////////////////////////////////////
+
+/**
+ * @brief forEachLeaf iterate on the leaf and apply the function
+ * @param function
+ */
+template< template< class CellClass, class ContainerClass, class LeafClass, class CellAllocatorClass > class FOctree,
+          class CellClass, class ContainerClass, class LeafClass, class CellAllocatorClass,
+          class FunctionTemplate>
+void forEachLeaf(FOctree<CellClass,ContainerClass,LeafClass,CellAllocatorClass>* tree, FunctionTemplate function){
+    #pragma omp parallel
+    {
+        #pragma omp single
+        {
+            typename FOctree<CellClass,ContainerClass,LeafClass,CellAllocatorClass>::Iterator octreeIterator(tree);
+            octreeIterator.gotoBottomLeft();
+
+            do{
+                LeafClass* lf = octreeIterator.getCurrentLeaf();
+                #pragma omp task firstprivate(lf)
+                {
+                    function(lf);
+                }
+            } while(octreeIterator.moveRight());
+
+            #pragma omp taskwait
+        }
+    }
+}
+
+/**
+ * @brief forEachLeaf iterate on the cell and apply the function
+ * @param function
+ */
+template< template< class CellClass, class ContainerClass, class LeafClass, class CellAllocatorClass > class FOctree,
+          class CellClass, class ContainerClass, class LeafClass, class CellAllocatorClass,
+          class FunctionTemplate>
+void forEachCell(FOctree<CellClass,ContainerClass,LeafClass,CellAllocatorClass>* tree, FunctionTemplate function){
+    #pragma omp parallel
+    {
+        #pragma omp single
+        {
+            typename FOctree<CellClass,ContainerClass,LeafClass,CellAllocatorClass>::Iterator octreeIterator(tree);
+            octreeIterator.gotoBottomLeft();
+
+            typename FOctree<CellClass,ContainerClass,LeafClass,CellAllocatorClass>::Iterator avoidGoLeft(octreeIterator);
+
+            for(int idx = tree->getHeight()-1 ; idx >= 1 ; --idx ){
+                do{
+                    CellClass* cell = octreeIterator.getCurrentCell();
+                    #pragma omp task firstprivate(cell)
+                    {
+                        function(cell);
+                    }
+                } while(octreeIterator.moveRight());
+                avoidGoLeft.moveUp();
+                octreeIterator = avoidGoLeft;
+            }
+
+            #pragma omp taskwait
+        }
+    }
+}
+
+/**
+ * @brief forEachLeaf iterate on the cell and apply the function
+ * @param function
+ */
+template< template< class CellClass, class ContainerClass, class LeafClass, class CellAllocatorClass > class FOctree,
+          class CellClass, class ContainerClass, class LeafClass, class CellAllocatorClass,
+          class FunctionTemplate>
+void forEachCellWithLevel(FOctree<CellClass,ContainerClass,LeafClass,CellAllocatorClass>* tree, FunctionTemplate function){
+    #pragma omp parallel
+    {
+        #pragma omp single
+        {
+            typename FOctree<CellClass,ContainerClass,LeafClass,CellAllocatorClass>::Iterator octreeIterator(tree);
+            octreeIterator.gotoBottomLeft();
+
+            typename FOctree<CellClass,ContainerClass,LeafClass,CellAllocatorClass>::Iterator avoidGoLeft(octreeIterator);
+
+            for(int idx = tree->getHeight()-1 ; idx >= 1 ; --idx ){
+                do{
+                    CellClass* cell = octreeIterator.getCurrentCell();
+                    #pragma omp task firstprivate(cell, idx)
+                    {
+                        function(cell,idx);
+                    }
+                } while(octreeIterator.moveRight());
+                avoidGoLeft.moveUp();
+                octreeIterator = avoidGoLeft;
+            }
+
+            #pragma omp taskwait
+        }
+    }
+}
+
+/**
+ * @brief forEachLeaf iterate on the cell and apply the function
+ * @param function
+ */
+template< template< class CellClass, class ContainerClass, class LeafClass, class CellAllocatorClass > class FOctree,
+          class CellClass, class ContainerClass, class LeafClass, class CellAllocatorClass,
+          class FunctionTemplate>
+void forEachCellLeaf(FOctree<CellClass,ContainerClass,LeafClass,CellAllocatorClass>* tree, FunctionTemplate function){
+    #pragma omp parallel
+    {
+        #pragma omp single
+        {
+            typename FOctree<CellClass,ContainerClass,LeafClass,CellAllocatorClass>::Iterator octreeIterator(tree);
+            octreeIterator.gotoBottomLeft();
+
+            do{
+                CellClass* cell = octreeIterator.getCurrentCell();
+                LeafClass* lf = octreeIterator.getCurrentLeaf();
+                #pragma omp task firstprivate(cell, lf)
+                {
+                    function(cell,lf);
+                }
+            } while(octreeIterator.moveRight());
+
+            #pragma omp taskwait
+        }
+    }
+}
+
+}
+
+#endif // FPARFOREACHOCTREE_HPP

Src/Utils/FParObject.hpp

+#ifndef FPAROBJECT_HPP
+#define FPAROBJECT_HPP
+
+#include <omp.h>
+#include <functional>
+
+/**
+ * This class is a way to hide concurent access to the user.
+ * Especially with the FParForEachOctree functions.
+ * Please refer to the testOctreeParallelFuncteur to see an example.
+ */
+template <class SharedClass>
+class FParObject{
+    const int nbThreads;
+    SharedClass* objects;
+
+    FParObject(const FParObject&) = delete;
+    FParObject& operator=(const FParObject&) = delete;
+
+public:
+
+    template <class... ArgsClass>
+    explicit FParObject(ArgsClass... args) : nbThreads(omp_get_max_threads()), objects(nullptr) {
+        objects = reinterpret_cast<SharedClass*>(new unsigned char[sizeof(SharedClass)*nbThreads]);
+        for(int idxThread = 0 ; idxThread < nbThreads ; ++idxThread){
+            new (&objects[idxThread]) SharedClass(args...);
+        }
+    }
+
+    ~FParObject(){
+        for(int idxThread = 0 ; idxThread < nbThreads ; ++idxThread){
+            objects[idxThread].~SharedClass();
+        }
+        delete[] reinterpret_cast<unsigned char*>(objects);
+    }
+
+    SharedClass& getMine(){
+        return objects[omp_get_thread_num()];
+    }
+
+    const SharedClass& getMine() const {
+        return objects[omp_get_thread_num()];
+    }
+
+    void applyToAll(std::function<void(SharedClass&)> func){
+        for(int idxThread = 0 ; idxThread < nbThreads ; ++idxThread){
+            func(objects[idxThread]);
+        }
+    }
+
+    void applyToAllPar(std::function<void(SharedClass&)> func){
+        #pragma omp parallel for num_threads(nbThreads) schedule(static)
+        for(int idxThread = 0 ; idxThread < nbThreads ; ++idxThread){
+            func(objects[idxThread]);
+        }
+    }
+
+    SharedClass reduce(std::function<SharedClass(const SharedClass&, const SharedClass&)> func){
+        if(nbThreads == 0){
+            return SharedClass();
+        }
+        else if(nbThreads == 1){
+            return objects[0];
+        }
+        else{
+            SharedClass res = func(objects[0], objects[1]);
+            for(int idxThread = 2 ; idxThread < nbThreads ; ++idxThread){
+                res = func(res, objects[idxThread]);
+            }
+            return res;
+        }
+    }
+
+    void set(const SharedClass& val){
+        #pragma omp parallel for num_threads(nbThreads) schedule(static)
+        for(int idxThread = 0 ; idxThread < nbThreads ; ++idxThread){
+            objects[idxThread] = val;
+        }
+    }
+
+    FParObject& operator=(const SharedClass& val){
+        set(val);
+        return *this;
+    }
+};
+
+
+template <class SharedClass>
+class FParArray{
+    const int nbThreads;
+    SharedClass** arrays;
+    size_t arraySize;
+
+    FParArray(const FParArray&) = delete;
+    FParArray& operator=(const FParArray&) = delete;
+
+public:
+
+    explicit FParArray(const size_t inArraySize = 0)
+        : nbThreads(omp_get_max_threads()), arrays(nullptr), arraySize(inArraySize) {
+        arrays = new SharedClass*[nbThreads];
+        #pragma omp parallel for schedule(static)
+        for(int idxThread = 0 ; idxThread < nbThreads ; ++idxThread){
+            arrays[idxThread] = new SharedClass[inArraySize];
+            if(std::is_pod<SharedClass>::value){
+                memset(arrays[idxThread], 0, sizeof(SharedClass)*inArraySize);
+            }
+        }
+    }
+
+    ~FParArray(){
+        for(int idxThread = 0 ; idxThread < nbThreads ; ++idxThread){
+            delete[] arrays[idxThread];
+        }
+        delete[] arrays;
+    }
+
+    void resize(const size_t inArraySize){
+        if(inArraySize != arraySize){
+            for(int idxThread = 0 ; idxThread < nbThreads ; ++idxThread){
+                delete[] arrays[idxThread];
+            }
+
+            arraySize = inArraySize;
+
+            #pragma omp parallel for schedule(static)
+            for(int idxThread = 0 ; idxThread < nbThreads ; ++idxThread){
+                arrays[idxThread] = new SharedClass[inArraySize];
+                if(std::is_pod<SharedClass>::value){
+                    memset(arrays[idxThread], 0, sizeof(SharedClass)*inArraySize);
+                }
+            }
+        }
+    }
+
+    SharedClass* getMine(){
+        return arrays[omp_get_thread_num()];
+    }
+
+    const SharedClass* getMine() const {
+        return arrays[omp_get_thread_num()];
+    }
+
+    void applyToAll(std::function<void(SharedClass*)> func){
+        for(int idxThread = 0 ; idxThread < nbThreads ; ++idxThread){
+            func(arrays[idxThread]);
+        }
+    }
+
+    void applyToAllPar(std::function<void(SharedClass*)> func){
+        #pragma omp parallel for num_threads(nbThreads) schedule(static)
+        for(int idxThread = 0 ; idxThread < nbThreads ; ++idxThread){
+            func(arrays[idxThread]);
+        }
+    }
+
+    void reduceArray(SharedClass* resArray, std::function<SharedClass(const SharedClass&, const SharedClass&)> func){
+        if(nbThreads == 0){
+        }
+        else if(nbThreads == 1){
+            if(std::is_pod<SharedClass>::value){
+                memcpy(resArray, arrays[0], sizeof(SharedClass)*arraySize);
+            }
+            else{
+                for(int idxVal = 0 ; idxVal < arraySize ; ++idxVal){
+                    resArray = arrays[0][idxVal];
+                }
+            }
+        }
+        else{
+            for(int idxVal = 0 ; idxVal < arraySize ; ++idxVal){
+                resArray[idxVal] = func(arrays[0][idxVal], arrays[1][idxVal]);
+            }
+            for(int idxThread = 2 ; idxThread < nbThreads ; ++idxThread){
+                for(int idxVal = 0 ; idxVal < arraySize ; ++idxVal){
+                    resArray[idxVal] = func(resArray[idxVal], arrays[idxThread][idxVal]);
+                }
+            }
+        }
+    }
+
+    void copyArray(const SharedClass* array){
+        #pragma omp parallel for num_threads(nbThreads) schedule(static)
+        for(int idxThread = 0 ; idxThread < nbThreads ; ++idxThread){
+            if(std::is_pod<SharedClass>::value){
+                memcpy(arrays[idxThread], array, sizeof(SharedClass)*arraySize);
+            }
+            else{
+                for(int idxVal = 0 ; idxVal < arraySize ; ++idxVal){
+                    arrays[idxThread][idxVal] = array[idxVal];
+                }
+            }
+        }
+    }
+};
+
+#endif // FPAROBJECT_HPP

Tests/Utils/testOctreeParallelFuncteur.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 <iostream>
+
+#include <cstdio>
+#include <cstdlib>
+#include <time.h>
+
+#include "../../Src/Utils/FTic.hpp"
+#include "../../Src/Utils/FParameters.hpp"
+
+#include "../../Src/Containers/FOctree.hpp"
+#include "../../Src/Containers/FParForEachOctree.hpp"
+#include "../../Src/Containers/FVector.hpp"
+
+#include "../../Src/Utils/FAssert.hpp"
+#include "../../Src/Utils/FPoint.hpp"
+#include "../../Src/Utils/FParObject.hpp"
+
+#include "../../Src/Components/FBasicParticleContainer.hpp"
+#include "../../Src/Components/FBasicCell.hpp"
+#include "../../Src/Components/FSimpleLeaf.hpp"
+
+#include "../../Src/Files/FRandomLoader.hpp"
+
+#include "../../Src/Utils/FParameterNames.hpp"
+
+/**
+* In this file we show how to use octree's functeur
+*/
+
+int main(int argc, char ** argv){
+    FHelpDescribeAndExit(argc, argv,
+                         "Show how to use an octree functeur parallelized (only the code is interesting)",
+                         FParameterDefinitions::NbParticles);
+
+    typedef FBasicCell CellClass;
+    typedef FBasicParticleContainer<0>      ContainerClass;
+    typedef FSimpleLeaf< ContainerClass >                     LeafClass;
+    typedef FOctree<CellClass, ContainerClass , LeafClass >  OctreeClass;
+    ///////////////////////What we do/////////////////////////////
+    std::cout << ">> This executable is useless to execute.\n";
+    std::cout << ">> It is only interesting to wath the code to understand\n";
+    std::cout << ">> how to use the Octree\n";
+    //////////////////////////////////////////////////////////////
+    const int NbPart = FParameters::getValue(argc,argv,FParameterDefinitions::NbParticles.options, 2000);
+    FTic counter;
+
+    FRandomLoader loader(NbPart, 1, FPoint(0.5,0.5,0.5), 1);
+    OctreeClass tree(10, 3, loader.getBoxWidth(), loader.getCenterOfBox());
+
+    // -----------------------------------------------------
+    std::cout << "Creating and inserting " << NbPart << " particles ..." << std::endl;
+    counter.tic();
+
+    {
+        FPoint particlePosition;
+        for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
+            loader.fillParticle(&particlePosition);
+            tree.insert(particlePosition);
+        }
+    }
+
+    counter.tac();
+    std::cout << "Done  " << "(" << counter.elapsed() << ")." << std::endl;
+    // -----------------------------------------------------
+
+    // Call a function on each leaf
+    FParObject<long> nbParticles(0);
+    FParForEachOctree::forEachLeaf(&tree, [&](LeafClass* leaf){
+        nbParticles.getMine() += leaf->getSrc()->getNbParticles();
+    });
+    const long totalNbParticles = nbParticles.reduce([](long v1, long v2) -> long {return v1 + v2;});
+    std::cout << "There are " << totalNbParticles << " particles " << std::endl;
+
+    // Call a function on each cell
+    FParObject<long> nbCells;
+    nbCells = 0;
+    FParForEachOctree::forEachCell(&tree, [&nbCells](CellClass* /*cell*/){
+        nbCells.getMine() += 1;
+    });
+    std::cout << "There are " << nbCells.reduce([](long v1, long v2) -> long {return v1 + v2;}) << " cells " << std::endl;
+
+    // To get cell and particles at leaf level
+    FParForEachOctree::forEachCellLeaf(&tree, [&](CellClass* cell, LeafClass* /*leaf*/){
+        cell->resetToInitialState();
+    });
+
+    // If we need an array of long
+    FParArray<long> arrayExample(10);
+    arrayExample.resize(50);
+
+    FParForEachOctree::forEachCellLeaf(&tree, [&](CellClass* /*cell*/, LeafClass* /*leaf*/){
+        arrayExample.getMine()[0] += 10;
+    });
+
+    return 0;
+}
+
+
+
+