Commit 794015bd authored by Berenger Bramas's avatar Berenger Bramas

debug mpi from tinker branch feedback (it does not include the tinker part)

parent f2533aae
......@@ -164,6 +164,8 @@ int main(int argc, char* argv[])
tree.getBoxWidth(),
tree.getHeight(), &finalParticles,&balancer);
std::cout << "Local nb particles after sort "
<< finalParticles.getSize() << std::endl;
for(FSize idx = 0 ; idx < finalParticles.getSize(); ++idx){
tree.insert(finalParticles[idx].position,finalParticles[idx].indexInFile,finalParticles[idx].physicalValue);
}
......
......@@ -4,6 +4,9 @@ the more restrictive has to be used.
ScalFmm est régi par la licence CeCILL-C & LGPL, en cas de conflit
la plus restrictive prime.
Folders under Addons might have separate Licence, in such case
one can find a dedicated Licence file where appropriate.
//////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////
......
......@@ -81,7 +81,10 @@ public:
pack.elementTo = Min(allObjectives[idxProc].second , myCurrentInterval.second) - myCurrentInterval.first;
// Next time give from the previous end
currentElement = pack.elementTo;
packToSend.push_back(pack);
if(pack.elementTo - pack.elementFrom != 0){
packToSend.push_back(pack);
}
// Progress
idxProc += 1;
}
......
// ===================================================================================
// 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".
// ===================================================================================
#ifndef FPARTITIONSMAPPING_HPP
#define FPARTITIONSMAPPING_HPP
#include "Utils/FGlobal.hpp"
#include "Utils/FMpi.hpp"
#include "Containers/FVector.hpp"
#include "FLeafBalance.hpp"
#include "Files/FMpiTreeBuilder.hpp"
template <class FReal>
class FPartitionsMapping {
protected:
FMpi::FComm comm;
//! The number of particles from the initial decomposition
FSize nbParticlesInitial;
//! The number of particles from the scalfmm decomposition
FSize nbParticlesWorking;
std::unique_ptr<FSize[]> nbParticlesSentToOthers;
std::unique_ptr<FSize[]> offsetNbParticlesSentToOthers;
std::unique_ptr<FSize[]> nbParticlesRecvFromOthers;
std::unique_ptr<FSize[]> offsetNbParticlesRecvFromOthers;
std::unique_ptr<FSize[]> mappingToOthers;
public:
FPartitionsMapping(const FMpi::FComm& inComm)
: comm(inComm), nbParticlesInitial(0), nbParticlesWorking(0) {
}
void setComm(const FMpi::FComm& inComm){
comm = inComm;
}
template< int NbPhysicalValuesPerPart>
struct TestParticle{
FPoint<FReal> position;
std::array<FReal, NbPhysicalValuesPerPart> physicalValues;
FSize localIndex;
int initialProcOwner;
const FPoint<FReal>& getPosition() const {
return position;
}
};
template< int NbPhysicalValuesPerPart, class FillerClass>
FVector<TestParticle<NbPhysicalValuesPerPart>> distributeParticles(const FSize inNbParticles,
const FPoint<FReal>& centerOfBox, const FReal boxWidth,
const int TreeHeight, FillerClass filler){
nbParticlesInitial = inNbParticles;
////////////////////////////////////////////////////////
std::unique_ptr<TestParticle<NbPhysicalValuesPerPart>[]> initialParticles(new TestParticle<NbPhysicalValuesPerPart>[inNbParticles]);
// Create the array to distribute
for(int idxPart = 0 ; idxPart < nbParticlesInitial ; ++idxPart){
filler(idxPart, &initialParticles[idxPart].position, &initialParticles[idxPart].physicalValues);
initialParticles[idxPart].localIndex = idxPart;
initialParticles[idxPart].initialProcOwner = comm.processId();
}
FVector<TestParticle<NbPhysicalValuesPerPart>> finalParticles;
FLeafBalance balancer;
FMpiTreeBuilder< FReal,TestParticle<NbPhysicalValuesPerPart> >::DistributeArrayToContainer(comm,initialParticles.get(),
nbParticlesInitial,
centerOfBox,
boxWidth,
TreeHeight,
&finalParticles, &balancer);
FQuickSort<TestParticle<NbPhysicalValuesPerPart>,FSize>::QsOmp(finalParticles.data(), finalParticles.getSize(),
[](const TestParticle<NbPhysicalValuesPerPart>& p1,
const TestParticle<NbPhysicalValuesPerPart>& p2){
return p1.initialProcOwner < p2.initialProcOwner
|| (p1.initialProcOwner == p2.initialProcOwner
&& p1.localIndex < p2.localIndex);
});
////////////////////////////////////////////////////////
nbParticlesWorking = finalParticles.getSize();
nbParticlesRecvFromOthers.reset(new FSize[comm.processCount()]);
memset(nbParticlesRecvFromOthers.get(), 0 , sizeof(FSize)*comm.processCount());
for(FSize idxPart = 0 ; idxPart < finalParticles.getSize() ; ++idxPart){
// Count the particles received from each proc
nbParticlesRecvFromOthers[finalParticles[idxPart].initialProcOwner] += 1;
}
offsetNbParticlesRecvFromOthers.reset(new FSize[comm.processCount()+1]);
offsetNbParticlesRecvFromOthers[0] = 0;
for(int idxProc = 0 ; idxProc < comm.processCount() ; ++idxProc){
offsetNbParticlesRecvFromOthers[idxProc+1] = offsetNbParticlesRecvFromOthers[idxProc]
+ nbParticlesRecvFromOthers[idxProc];
}
////////////////////////////////////////////////////////
std::unique_ptr<FSize[]> nbParticlesRecvFromOthersAllAll(new FSize[comm.processCount()*comm.processCount()]);
// Exchange how many each proc receive from aother
FMpi::MpiAssert( MPI_Allgather( nbParticlesRecvFromOthers.get(), comm.processCount(), FMpi::GetType(FSize()),
nbParticlesRecvFromOthersAllAll.get(), comm.processCount(),
FMpi::GetType(FSize()), comm.getComm()), __LINE__ );
////////////////////////////////////////////////////////
nbParticlesSentToOthers.reset(new FSize[comm.processCount()]);
FSize checkerSent = 0;
offsetNbParticlesSentToOthers.reset(new FSize[comm.processCount()+1]);
offsetNbParticlesSentToOthers[0] = 0;
for(int idxProc = 0 ; idxProc < comm.processCount() ; ++idxProc){
nbParticlesSentToOthers[idxProc] = nbParticlesRecvFromOthersAllAll[comm.processCount()*idxProc + comm.processId()];
checkerSent += nbParticlesSentToOthers[idxProc];
offsetNbParticlesSentToOthers[idxProc+1] = offsetNbParticlesSentToOthers[idxProc]
+ nbParticlesSentToOthers[idxProc];
}
// I must have send what I was owning at the beginning
FAssertLF(checkerSent == nbParticlesInitial);
////////////////////////////////////////////////////////
std::unique_ptr<FSize[]> localIdsRecvOrdered(new FSize[nbParticlesWorking]);
// We list the local id in order of the different proc
for(FSize idxPart = 0 ; idxPart < finalParticles.getSize() ; ++idxPart){
const int procOwner = finalParticles[idxPart].initialProcOwner;
localIdsRecvOrdered[idxPart] = finalParticles[idxPart].localIndex;
FAssertLF(offsetNbParticlesRecvFromOthers[procOwner] <= idxPart
&& idxPart < offsetNbParticlesRecvFromOthers[procOwner+1]);
}
////////////////////////////////////////////////////////
std::unique_ptr<FSize[]> localIdsSendOrdered(new FSize[nbParticlesInitial]);
std::unique_ptr<MPI_Request[]> requests(new MPI_Request[comm.processCount()*2]);
int iterRequest = 0;
for(int idxProc = 0 ; idxProc < comm.processCount() ; ++idxProc){
if(idxProc == comm.processId()){
FAssertLF(nbParticlesRecvFromOthers[idxProc] == nbParticlesSentToOthers[idxProc]);
memcpy(&localIdsSendOrdered[offsetNbParticlesSentToOthers[idxProc]],
&localIdsRecvOrdered[offsetNbParticlesRecvFromOthers[idxProc]],
sizeof(FSize)*nbParticlesRecvFromOthers[idxProc]);
}
else{
const FSize nbRecvFromProc = nbParticlesRecvFromOthers[idxProc];
if(nbRecvFromProc){
FMpi::MpiAssert( MPI_Isend(&localIdsRecvOrdered[offsetNbParticlesRecvFromOthers[idxProc]],
int(nbRecvFromProc),
FMpi::GetType(FSize()), idxProc,
99, comm.getComm(), &requests[iterRequest++]), __LINE__ );
}
const FSize nbSendToProc = nbParticlesSentToOthers[idxProc];
if(nbSendToProc){
FMpi::MpiAssert( MPI_Irecv(&localIdsSendOrdered[offsetNbParticlesSentToOthers[idxProc]],
int(nbSendToProc),
FMpi::GetType(FSize()), idxProc,
99, comm.getComm(), &requests[iterRequest++]), __LINE__ );
}
}
}
FMpi::MpiAssert( MPI_Waitall( iterRequest, requests.get(), MPI_STATUSES_IGNORE), __LINE__ );
////////////////////////////////////////////////////////
mappingToOthers.reset(new FSize[nbParticlesInitial]);
for(FSize idxPart = 0; idxPart < nbParticlesInitial ; ++idxPart){
mappingToOthers[localIdsSendOrdered[idxPart]] = idxPart;
}
return std::move(finalParticles);
}
////////////////////////////////////////////////////////
// physicalValues must be of size nbParticlesInitial
template< int NbPhysicalValuesPerPart>
std::unique_ptr<std::array<FReal, NbPhysicalValuesPerPart>[]> distributeData(
const std::array<FReal, NbPhysicalValuesPerPart> physicalValues[]){
std::unique_ptr<std::array<FReal, NbPhysicalValuesPerPart>[]> physicalValuesRorder(
new std::array<FReal, NbPhysicalValuesPerPart>[nbParticlesInitial]);
for(FSize idxPart = 0; idxPart < nbParticlesInitial ; ++idxPart){
physicalValuesRorder[mappingToOthers[idxPart]] = physicalValues[idxPart];
}
// Allocate the array to store the physical values of my working interval
std::unique_ptr<std::array<FReal, NbPhysicalValuesPerPart>[]> recvPhysicalValues(new std::array<FReal, NbPhysicalValuesPerPart>[nbParticlesWorking]);
std::unique_ptr<MPI_Request[]> requests(new MPI_Request[comm.processCount()*2]);
int iterRequest = 0;
for(int idxProc = 0 ; idxProc < comm.processCount() ; ++idxProc){
if(idxProc == comm.processId()){
FAssertLF(nbParticlesRecvFromOthers[idxProc] == nbParticlesSentToOthers[idxProc]);
memcpy(&recvPhysicalValues[offsetNbParticlesRecvFromOthers[idxProc]],
&physicalValuesRorder[offsetNbParticlesSentToOthers[idxProc]],
sizeof(std::array<FReal, NbPhysicalValuesPerPart>)*nbParticlesRecvFromOthers[idxProc]);
}
else{
const FSize nbSendToProc = nbParticlesSentToOthers[idxProc];
if(nbSendToProc){
FMpi::MpiAssert( MPI_Isend(
const_cast<std::array<FReal, NbPhysicalValuesPerPart>*>(&physicalValuesRorder[offsetNbParticlesSentToOthers[idxProc]]),
int(nbSendToProc*sizeof(std::array<FReal, NbPhysicalValuesPerPart>)),
MPI_BYTE, idxProc,
2222, comm.getComm(), &requests[iterRequest++]), __LINE__ );;
}
const FSize nbRecvFromProc = nbParticlesRecvFromOthers[idxProc];
if(nbRecvFromProc){
FMpi::MpiAssert( MPI_Irecv(
(void*)&recvPhysicalValues[offsetNbParticlesRecvFromOthers[idxProc]],
int(nbRecvFromProc*sizeof(std::array<FReal, NbPhysicalValuesPerPart>)),
MPI_INT, idxProc,
2222, comm.getComm(), &requests[iterRequest++]), __LINE__ );;
}
}
}
FMpi::MpiAssert( MPI_Waitall( iterRequest, requests.get(), MPI_STATUSES_IGNORE), __LINE__ );
return std::move(recvPhysicalValues);
}
////////////////////////////////////////////////////////
// resValues must be of size nbParticlesWorking
template< int NbResValuesPerPart>
std::unique_ptr<std::array<FReal, NbResValuesPerPart>[]> getResultingData(
const std::array<FReal, NbResValuesPerPart> resValues[]){
// First allocate the array to store the result
std::unique_ptr<std::array<FReal, NbResValuesPerPart>[]> recvPhysicalValues(
new std::array<FReal, NbResValuesPerPart>[nbParticlesInitial]);
std::unique_ptr<MPI_Request[]> requests(new MPI_Request[comm.processCount()*2]);
int iterRequest = 0;
for(int idxProc = 0 ; idxProc < comm.processCount() ; ++idxProc){
if(idxProc == comm.processId()){
FAssertLF(nbParticlesRecvFromOthers[idxProc] == nbParticlesSentToOthers[idxProc]);
memcpy(&recvPhysicalValues[offsetNbParticlesSentToOthers[idxProc]],
&resValues[offsetNbParticlesRecvFromOthers[idxProc]],
sizeof(std::array<FReal, NbResValuesPerPart>)*nbParticlesRecvFromOthers[idxProc]);
}
else{
// I originaly receive nbRecvFromProc, so I should
// send nbRecvFromProc back to the real owner
const FSize nbRecvFromProc = nbParticlesRecvFromOthers[idxProc];
if(nbRecvFromProc){
FMpi::MpiAssert( MPI_Isend(
const_cast<std::array<FReal, NbResValuesPerPart>*>(&resValues[offsetNbParticlesRecvFromOthers[idxProc]]),
int(nbRecvFromProc*sizeof(std::array<FReal, NbResValuesPerPart>)),
MPI_BYTE, idxProc,
1111, comm.getComm(), &requests[iterRequest++]), __LINE__ );;
}
// I sent nbSendToProc to idxProc,
// so I should receive nbSendToProc in my interval
const FSize nbSendToProc = nbParticlesSentToOthers[idxProc];
if(nbSendToProc){
FMpi::MpiAssert( MPI_Irecv(
&recvPhysicalValues[offsetNbParticlesSentToOthers[idxProc]],
int(nbSendToProc*sizeof(std::array<FReal, NbResValuesPerPart>)),
MPI_BYTE, idxProc,
1111, comm.getComm(), &requests[iterRequest++]), __LINE__ );;
}
}
}
FMpi::MpiAssert( MPI_Waitall( iterRequest, requests.get(), MPI_STATUSES_IGNORE), __LINE__ );
std::unique_ptr<std::array<FReal, NbResValuesPerPart>[]> recvPhysicalValuesOrder(
new std::array<FReal, NbResValuesPerPart>[nbParticlesInitial]);
for(FSize idxPart = 0; idxPart < nbParticlesInitial ; ++idxPart){
recvPhysicalValuesOrder[idxPart] = recvPhysicalValues[mappingToOthers[idxPart]];
}
return std::move(recvPhysicalValuesOrder);
}
////////////////////////////////////////////////////////
FSize getNbParticlesWorking() const{
return nbParticlesWorking;
}
FSize getMappingResultToLocal(const FSize inIdx) const{
return mappingToOthers[inIdx];
}
};
#endif
......@@ -1685,6 +1685,10 @@ public:
* @param function
*/
void forEachLeaf(std::function<void(LeafClass*)> function){
if(isEmpty()){
return;
}
Iterator octreeIterator(this);
octreeIterator.gotoBottomLeft();
......@@ -1698,6 +1702,10 @@ public:
* @param function
*/
void forEachCell(std::function<void(CellClass*)> function){
if(isEmpty()){
return;
}
Iterator octreeIterator(this);
octreeIterator.gotoBottomLeft();
......@@ -1717,6 +1725,10 @@ public:
* @param function
*/
void forEachCellWithLevel(std::function<void(CellClass*,const int)> function){
if(isEmpty()){
return;
}
Iterator octreeIterator(this);
octreeIterator.gotoBottomLeft();
......@@ -1736,6 +1748,10 @@ public:
* @param function
*/
void forEachCellLeaf(std::function<void(CellClass*,LeafClass*)> function){
if(isEmpty()){
return;
}
Iterator octreeIterator(this);
octreeIterator.gotoBottomLeft();
......
......@@ -72,17 +72,20 @@ private:
OctreeClass* const tree; ///< The octree to work on
KernelClass** kernels; ///< The kernels
const FMpi::FComm& comm; ///< MPI comm
const FMpi::FComm comm; ///< MPI comm
FMpi::FComm fcomCompute;
/// Used to store pointers to cells/leafs to work with
typename OctreeClass::Iterator* iterArray;
typename OctreeClass::Iterator* iterArray;
/// Used to store pointers to cells/leafs to send/rcv
typename OctreeClass::Iterator* iterArrayComm;
int numberOfLeafs; ///< To store the size at the previous level
const int MaxThreads; ///< Max number of thread allowed by openmp
const int nbProcess; ///< Process count
const int idProcess; ///< Current process id
const int nbProcessOrig; ///< Process count
const int idProcessOrig; ///< Current process id
int nbProcess; ///< Process count
int idProcess; ///< Current process id
const int OctreeHeight; ///< Tree height
const int userChunkSize;
......@@ -149,12 +152,15 @@ public:
tree(inTree),
kernels(nullptr),
comm(inComm),
fcomCompute(inComm),
iterArray(nullptr),
iterArrayComm(nullptr),
numberOfLeafs(0),
MaxThreads(FEnv::GetValue("SCALFMM_ALGO_NUM_THREADS",omp_get_max_threads())),
nbProcess(inComm.processCount()),
idProcess(inComm.processId()),
nbProcessOrig(inComm.processCount()),
idProcessOrig(inComm.processId()),
nbProcess(0),
idProcess(0),
OctreeHeight(tree->getHeight()),
userChunkSize(inUserChunkSize),
leafLevelSeparationCriteria(inLeafLevelSeperationCriteria),
......@@ -164,9 +170,9 @@ public:
FAssertLF(leafLevelSeparationCriteria < 3, "Separation criteria should be < 3");
this->kernels = new KernelClass*[MaxThreads];
#pragma omp parallel num_threads(MaxThreads)
#pragma omp parallel num_threads(MaxThreads)
{
#pragma omp critical (InitFFmmAlgorithmThreadProc)
#pragma omp critical (InitFFmmAlgorithmThreadProc)
{
this->kernels[omp_get_thread_num()] = new KernelClass(*inKernels);
}
......@@ -175,7 +181,7 @@ public:
FAbstractAlgorithm::setNbLevelsInTree(tree->getHeight());
FLOG(FLog::Controller << "FFmmAlgorithmThreadProc\n");
FLOG(FLog::Controller << "Max threads = " << MaxThreads << ", Procs = " << nbProcess << ", I am " << idProcess << ".\n");
FLOG(FLog::Controller << "Max threads = " << MaxThreads << ", Procs = " << nbProcessOrig << ", I am " << idProcessOrig << ".\n");
FLOG(FLog::Controller << "Chunck Size = " << userChunkSize << "\n");
}
......@@ -196,136 +202,158 @@ protected:
* Call this function to run the complete algorithm
*/
void executeCore(const unsigned operationsToProceed) override {
// Count leaf
// We are not involve if the tree is empty
const int iHaveParticles = (!tree->isEmpty());
std::unique_ptr<int[]> hasParticles(new int[comm.processCount()]);
FMpi::Assert( MPI_Allgather(const_cast<int*>(&iHaveParticles), 1,MPI_INT,
hasParticles.get(), 1, MPI_INT,
comm.getComm()), __LINE__);
fcomCompute = FMpi::FComm(comm);
fcomCompute.groupReduce(hasParticles.get());
if(iHaveParticles){
nbProcess = fcomCompute.processCount();
idProcess = fcomCompute.processId();
FLOG(FLog::Controller << "Max threads = " << MaxThreads << ", Procs = " << nbProcess << ", I am " << idProcess << ".\n");
// Count leaf
#ifdef SCALFMM_TRACE_ALGO
eztrace_resume();
eztrace_resume();
#endif
this->numberOfLeafs = 0;
{
Interval myFullInterval;
{//Building the interval with the first and last leaves (and count the number of leaves)
typename OctreeClass::Iterator octreeIterator(tree);
octreeIterator.gotoBottomLeft();
myFullInterval.leftIndex = octreeIterator.getCurrentGlobalIndex();
do{
++this->numberOfLeafs;
} while(octreeIterator.moveRight());
myFullInterval.rightIndex = octreeIterator.getCurrentGlobalIndex();
}
// Allocate a number to store the pointer of the cells at a level
iterArray = new typename OctreeClass::Iterator[numberOfLeafs];
iterArrayComm = new typename OctreeClass::Iterator[numberOfLeafs];
FAssertLF(iterArray, "iterArray bad alloc");
FAssertLF(iterArrayComm, "iterArrayComm bad alloc");
// We get the leftIndex/rightIndex indexes from each procs
FMpi::MpiAssert( MPI_Allgather( &myFullInterval, sizeof(Interval), MPI_BYTE, intervals, sizeof(Interval), MPI_BYTE, comm.getComm()), __LINE__ );
// Build my intervals for all levels
std::unique_ptr<Interval[]> myIntervals(new Interval[OctreeHeight]);
// At leaf level we know it is the full interval
myIntervals[OctreeHeight - 1] = myFullInterval;
// We can estimate the interval for each level by using the parent/child relation
for(int idxLevel = OctreeHeight - 2 ; idxLevel >= 0 ; --idxLevel){
myIntervals[idxLevel].leftIndex = myIntervals[idxLevel+1].leftIndex >> 3;
myIntervals[idxLevel].rightIndex = myIntervals[idxLevel+1].rightIndex >> 3;
}
this->numberOfLeafs = 0;
{
Interval myFullInterval;
{//Building the interval with the first and last leaves (and count the number of leaves)
typename OctreeClass::Iterator octreeIterator(tree);
octreeIterator.gotoBottomLeft();
myFullInterval.leftIndex = octreeIterator.getCurrentGlobalIndex();
do{
++this->numberOfLeafs;
} while(octreeIterator.moveRight());
myFullInterval.rightIndex = octreeIterator.getCurrentGlobalIndex();
}
// Allocate a number to store the pointer of the cells at a level
iterArray = new typename OctreeClass::Iterator[numberOfLeafs];
iterArrayComm = new typename OctreeClass::Iterator[numberOfLeafs];
FAssertLF(iterArray, "iterArray bad alloc");
FAssertLF(iterArrayComm, "iterArrayComm bad alloc");
// We get the leftIndex/rightIndex indexes from each procs
FMpi::MpiAssert( MPI_Allgather( &myFullInterval, sizeof(Interval), MPI_BYTE, intervals, sizeof(Interval), MPI_BYTE, fcomCompute.getComm()), __LINE__ );
// Build my intervals for all levels
std::unique_ptr<Interval[]> myIntervals(new Interval[OctreeHeight]);
// At leaf level we know it is the full interval
myIntervals[OctreeHeight - 1] = myFullInterval;
// We can estimate the interval for each level by using the parent/child relation
for(int idxLevel = OctreeHeight - 2 ; idxLevel >= 0 ; --idxLevel){
myIntervals[idxLevel].leftIndex = myIntervals[idxLevel+1].leftIndex >> 3;
myIntervals[idxLevel].rightIndex = myIntervals[idxLevel+1].rightIndex >> 3;
}
// Process 0 uses the estimates as real intervals, but other processes
// should remove cells that belong to others
if(idProcess != 0){
//We test for each level if process on left (idProcess-1) own cell I thought I owned
typename OctreeClass::Iterator octreeIterator(tree);
octreeIterator.gotoBottomLeft();
octreeIterator.moveUp();
// Process 0 uses the estimates as real intervals, but other processes
// should remove cells that belong to others
if(idProcess != 0){
//We test for each level if process on left (idProcess-1) own cell I thought I owned
typename OctreeClass::Iterator octreeIterator(tree);
octreeIterator.gotoBottomLeft();
octreeIterator.moveUp();
// At h-1 the working limit is the parent of the right cell of the proc on the left
MortonIndex workingLimitAtLevel = intervals[idProcess-1].rightIndex >> 3;
// At h-1 the working limit is the parent of the right cell of the proc on the left
MortonIndex workingLimitAtLevel = intervals[idProcess-1].rightIndex >> 3;
// We check if we have no more work to do
int nullIntervalFromLevel = 0;
// We check if we have no more work to do
int nullIntervalFromLevel = 0;
for(int idxLevel = OctreeHeight - 2 ; idxLevel >= 1 && nullIntervalFromLevel == 0 ; --idxLevel){
while(octreeIterator.getCurrentGlobalIndex() <= workingLimitAtLevel){
if( !octreeIterator.moveRight() ){
// We cannot move right we are not owner of any more cell
nullIntervalFromLevel = idxLevel;
break;
for(int idxLevel = OctreeHeight - 2 ; idxLevel >= 1 && nullIntervalFromLevel == 0 ; --idxLevel){
while(octreeIterator.getCurrentGlobalIndex() <= workingLimitAtLevel){
if( !octreeIterator.moveRight() ){
// We cannot move right we are not owner of any more cell
nullIntervalFromLevel = idxLevel;
break;
}
}
// If we are responsible for some cells at this level keep the first index
if(nullIntervalFromLevel == 0){
myIntervals[idxLevel].leftIndex = octreeIterator.getCurrentGlobalIndex();
octreeIterator.moveUp();
workingLimitAtLevel >>= 3;
}
}
// If we are responsible for some cells at this level keep the first index
if(nullIntervalFromLevel == 0){
myIntervals[idxLevel].leftIndex = octreeIterator.getCurrentGlobalIndex();
octreeIterator.moveUp();
workingLimitAtLevel >>= 3;
// In case we are not responsible for any cells we put the leftIndex = rightIndex+1
for(int idxLevel = nullIntervalFromLevel ; idxLevel >= 1 ; --idxLevel){
myIntervals[idxLevel].leftIndex = myIntervals[idxLevel].rightIndex + 1;
}
}
// In case we are not responsible for any cells we put the leftIndex = rightIndex+1
for(int idxLevel = nullIntervalFromLevel ; idxLevel >= 1 ; --idxLevel){
myIntervals[idxLevel].leftIndex = myIntervals[idxLevel].rightIndex + 1;
}
}
// We get the leftIndex/rightIndex indexes from each procs
FMpi::MpiAssert( MPI_Allgather( myIntervals.get(), int(sizeof(Interval)) * OctreeHeight, MPI_BYTE,
workingIntervalsPerLevel, int(sizeof(Interval)) * OctreeHeight, MPI_BYTE, comm.getComm()), __LINE__ );
}
// We get the leftIndex/rightIndex indexes from each procs
FMpi::MpiAssert( MPI_Allgather( myIntervals.get(), int(sizeof(Interval)) * OctreeHeight, MPI_BYTE,
workingIntervalsPerLevel, int(sizeof(Interval)) * OctreeHeight, MPI_BYTE, fcomCompute.getComm()), __LINE__ );