Commit 72822426 authored by COULAUD Olivier's avatar COULAUD Olivier

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

# By bramas (9) and others
# Via aetcheve
* 'master' of git+ssh://scm.gforge.inria.fr//gitroot/scalfmm/scalfmm:
  fix some warnings
  Add a exec file to print the stats of the adaptive execution
  simple update of the adaptative test
  Remove old custom openmp barrier
  Manage the case when a node has all the same and so pivot equal max
  update the FMM api to make it compile from the latest changes in ScalFmm
  Add an fuse to protect from disabled fft
  Implemented adaptive routines for AdaptChebSym FMM kernel, Switched base kernel from dense to classic (fft accelerated) Unif FMM kernel in AdaptUnif FMM kernel.
  Debug the mpi quick sort in case the pivot is the max value
  changes and fix in utestMpiTreeBuilder
  Add an assert in the rotation matrix operator
  Add a division method
  Fix for the error occuring in M2M, still bugged
parents 59df215d 58613d08
......@@ -5,6 +5,7 @@ if(insource)
endif(insource)
project(Addons_fmmapi_scalfmm CXX)
ADD_DEFINITIONS( ${ScaLFMM_CXX_FLAGS})
# Active language
# -----------------------
......@@ -29,7 +30,11 @@ if(SCALFMM_ADDON_FMMAPI)
target_link_libraries( scalfmmapi scalfmm)
# Adding the entire project dir as an include dir
INCLUDE_DIRECTORIES( ${CMAKE_BINARY_DIR}/Src )
INCLUDE_DIRECTORIES(
${CMAKE_BINARY_DIR}/Src
${CMAKE_SOURCE_DIR}/Src
${SCALFMM_INCLUDES}
)
# Install lib
install( TARGETS scalfmmapi ARCHIVE DESTINATION lib )
......@@ -80,6 +85,7 @@ if(SCALFMM_ADDON_FMMAPI)
scalfmmapi
${BLAS_LIBRARIES}
${LAPACK_LIBRARIES}
${SCALFMM_LIBRARIES}
)
endif()
endforeach(exec)
......
......@@ -107,6 +107,7 @@ int FmmKernel_L2P(void *fmmCore, void* boxId);
int FmmKernel_M2M(void *fmmCore, void *boxIdFather, void *boxIdSon);
int FmmKernel_L2L(void *fmmCore, void *boxIdFather, void *boxIdSon);
int FmmKernel_M2L(void *fmmCore, void *boxIdSrc, void *boxIdDest);
int FmmKernel_P2P_inner(void *fmmCore, void *boxIdSrcDest);
int FmmKernel_P2P(void *fmmCore, void *boxIdSrc, void *boxIdDest); /* pas mutuel, i.e. on fait seulement dans 1 sens. */
......
......@@ -12,7 +12,6 @@
#include "../../Src/Containers/FVector.hpp"
#include "../../Src/Components/FSimpleLeaf.hpp"
#include "../../Src/Components/FBasicParticle.hpp"
#include "../../Src/Components/FBasicCell.hpp"
#include "../../Src/Utils/FPoint.hpp"
......@@ -23,27 +22,14 @@
#include "../../Src/Components/FBasicKernels.hpp"
#include "../../Src/Components/FBasicParticleContainer.hpp"
#ifdef ScalFMM_USE_MPI
#include "../../Src/Utils/FMpi.hpp"
#endif
#include "FmmApi.h"
template <class ParticleClass>
class CoreParticle : public ParticleClass {
int index;
public:
CoreParticle() : index(0) {
}
void setIndex(const int inIndex){
index = inIndex;
}
int getIndex() const {
return index;
}
};
template <class ContainerClass>
class CoreCell : public FBasicCell {
char* multipole;
......@@ -53,7 +39,7 @@ class CoreCell : public FBasicCell {
mutable ContainerClass* container;
public:
CoreCell() : multipole(0), local(0), level(0), container(0) {
CoreCell() : multipole(nullptr), local(nullptr), level(0), container(nullptr) {
}
void createArrays(const int multipoleSize, const int localSize){
multipole = new char[multipoleSize];
......@@ -112,8 +98,8 @@ public:
template< class ParticleClass, class CellClass, class ContainerClass>
class CoreKernel : public FAbstractKernels<ParticleClass,CellClass,ContainerClass> {
template< class CellClass, class ContainerClass>
class CoreKernel : public FAbstractKernels<CellClass,ContainerClass> {
void* fmmCore;
public:
......@@ -169,7 +155,7 @@ public:
ContainerClass* const FRestrict targets, const ContainerClass* const FRestrict /*sources*/,
ContainerClass* const neighbors[27], const int ){
CellClass* cell = (CellClass*)targets->getParentCell();
FmmKernel_P2P(fmmCore, cell, cell);
FmmKernel_P2P_inner(fmmCore, cell);
for(int idx = 0 ; idx < 27 ; ++idx){
if( neighbors[idx] ){
......@@ -188,7 +174,6 @@ public:
};
template <class ObjectClass>
class CoreVector {
mutable void* parentCell;
void* fields;
......@@ -198,8 +183,8 @@ class CoreVector {
FVector<FReal> positions;
FVector<int> indexes;
public:
CoreVector() : parentCell(0), fields(0), sizeOfField(0),
potentials(0), sizeOfPential(0){
CoreVector() : parentCell(nullptr), fields(nullptr), sizeOfField(0),
potentials(nullptr), sizeOfPential(0){
}
~CoreVector(){
......@@ -257,24 +242,24 @@ public:
return indexes.getSize();
}
void push(const ObjectClass& particle){
positions.push(particle.getPosition().getX());
positions.push(particle.getPosition().getY());
positions.push(particle.getPosition().getZ());
indexes.push(particle.getIndex());
void push(const FPoint& partPosition, const int partIndex){
positions.push(partPosition.getX());
positions.push(partPosition.getY());
positions.push(partPosition.getZ());
indexes.push(partIndex);
}
};
typedef CoreParticle<FBasicParticle> CoreParticleClass;
typedef CoreVector<CoreParticleClass> CoreContainerClass;
typedef CoreVector CoreContainerClass;
typedef CoreCell<CoreContainerClass> CoreCellClass;
typedef FSimpleLeaf<CoreParticleClass, CoreContainerClass > LeafClass;
typedef FOctree<CoreParticleClass, CoreCellClass, CoreContainerClass , LeafClass > OctreeClass;
typedef CoreKernel<CoreParticleClass, CoreCellClass, CoreContainerClass> CoreKernelClass;
typedef FSimpleLeaf<CoreContainerClass > LeafClass;
typedef FOctree<CoreCellClass, CoreContainerClass , LeafClass > OctreeClass;
typedef CoreKernel<CoreCellClass, CoreContainerClass> CoreKernelClass;
typedef FFmmAlgorithm<OctreeClass, CoreParticleClass, CoreCellClass, CoreContainerClass, CoreKernelClass, LeafClass > FmmClass;
typedef FFmmAlgorithmThread<OctreeClass, CoreParticleClass, CoreCellClass, CoreContainerClass, CoreKernelClass, LeafClass > FmmClassThread;
typedef FFmmAlgorithm<OctreeClass, CoreCellClass, CoreContainerClass, CoreKernelClass, LeafClass > FmmClass;
typedef FFmmAlgorithmThread<OctreeClass, CoreCellClass, CoreContainerClass, CoreKernelClass, LeafClass > FmmClassThread;
struct ScalFmmCoreHandle {
......@@ -554,12 +539,9 @@ int FmmCore_setPositions(void *fmmCore, int *nb, FReal *position) {
omp_set_num_threads(corehandle->config.nbThreads);
}
CoreParticleClass part;
for(int idxPart = 0 ; idxPart < (*nb) ; ++idxPart){
FReal* pos = &position[idxPart * 3];
part.setPosition(pos[0], pos[1], pos[2]);
part.setIndex(idxPart);
corehandle->octree->insert(part);
const FReal* pos = &position[idxPart * 3];
corehandle->octree->insert(FPoint(pos[0], pos[1], pos[2]),idxPart);
}
return FMMAPI_NO_ERROR;
......@@ -606,7 +588,7 @@ int FmmCore_doComputation(void *fmmCore) {
typename OctreeClass::Iterator octreeIterator(corehandle->octree);
octreeIterator.gotoBottomLeft();
for(int idxLevel = corehandle->config.treeHeight - 1 ; idxLevel > 1 ; --idxLevel ){
for(int idxLevel = corehandle->config.treeHeight - 1 ; idxLevel >= 1 ; --idxLevel ){
int multipoleSize;
int localSize;
......
This diff is collapsed.
......@@ -18,8 +18,6 @@
#include "../../Src/Files/FRandomLoader.hpp"
#include "../../Src/Components/FBasicParticle.hpp"
#include "../Src/FmmApi.h"
/** This program show an example of use of the fmm api
......@@ -39,7 +37,7 @@ int main(int argc, char ** argv){
//////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////
FRandomLoader<FBasicParticle> loader(NbPart, 1, FPoint(0.5,0.5,0.5), 1);
FRandomLoader loader(NbPart, 1, FPoint(0.5,0.5,0.5), 1);
void* FmmCoreHandle;
FmmCore_init(&FmmCoreHandle);
......@@ -66,16 +64,16 @@ int main(int argc, char ** argv){
FReal* positions = new FReal[nbPart*3];
{
FBasicParticle part;
FPoint part;
const FReal physicalValue = 0.1;
for(int idx = 0 ; idx < nbPart ; ++idx){
loader.fillParticle(part);
loader.fillParticle(&part);
potentials[idx] = physicalValue;
positions[3*idx] = part.getPosition().getX();
positions[3*idx+1] = part.getPosition().getY();
positions[3*idx+2] = part.getPosition().getZ();
positions[3*idx] = part.getX();
positions[3*idx+1] = part.getY();
positions[3*idx+2] = part.getZ();
}
FmmCore_setPositions(FmmCoreHandle, &nbPart, positions);
......@@ -135,15 +133,16 @@ int main(int argc, char ** argv){
dy *= inv_square_distance;
dz *= inv_square_distance;
fieldsdirect[4*idxTarget] += dx;
fieldsdirect[4*idxTarget+1] += dy;
fieldsdirect[4*idxTarget+2] += dz;
fieldsdirect[4*idxTarget+3] += inv_distance * potentials[idxSource];
fieldsdirect[4*idxTarget+1] += dx;
fieldsdirect[4*idxTarget+2] += dy;
fieldsdirect[4*idxTarget+3] += dz;
fieldsdirect[4*idxTarget+0] += inv_distance * potentials[idxSource];
fieldsdirect[4*idxSource+1] += -dx;
fieldsdirect[4*idxSource+2] += -dy;
fieldsdirect[4*idxSource+3] += -dz;
fieldsdirect[4*idxSource+0] += inv_distance * potentials[idxTarget];
fieldsdirect[4*idxSource] += -dx;
fieldsdirect[4*idxSource+1] += -dy;
fieldsdirect[4*idxSource+2] += -dz;
fieldsdirect[4*idxSource+3] += inv_distance * potentials[idxTarget];
}
}
......
This diff is collapsed.
......@@ -22,8 +22,7 @@
#include "Adaptative/FAdaptiveCell.hpp"
#include "Adaptative/FAdaptiveKernelWrapper.hpp"
#include "Adaptative/FAbstractAdaptiveKernel.hpp"
//#include "Kernels/Uniform/FUnifKernel.hpp"
#include "Kernels/Uniform/FUnifDenseKernel.hpp"
#include "Kernels/Uniform/FUnifKernel.hpp"
#include "Kernels/Uniform/FUnifM2LHandler.hpp"
class FTreeCoordinate;
......@@ -55,15 +54,21 @@ class FTreeCoordinate;
*/
template< class CellClass, class ContainerClass, class MatrixKernelClass, int ORDER, int NVALS = 1>
class FAdaptiveUnifKernel : FUnifDenseKernel<CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
class FAdaptiveUnifKernel : FUnifKernel<CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
, public FAbstractAdaptiveKernel<CellClass, ContainerClass> {
//
typedef FUnifDenseKernel<CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS> KernelBaseClass;
typedef FUnifKernel<CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS> KernelBaseClass;
enum {order = ORDER,
nnodes = KernelBaseClass::nnodes};
nnodes = TensorTraits<ORDER>::nnodes};
/// Needed for M2L operator
typedef FUnifM2LHandler<ORDER,MatrixKernelClass::Type> M2LHandlerClass;
const M2LHandlerClass M2LHandler;
// // If we choose to pre-assemble adaptive M2L operators
// // then we need to provide an adaptive M2L handler
// // and the transfer in Fourier space is straightforward.
// typedef FUnifM2LHandler<ORDER,MatrixKernelClass::Type> M2LHandlerClass;
// const M2LHandlerClass M2LHandler;
const MatrixKernelClass *const MatrixKernel;
public:
......@@ -82,11 +87,11 @@ public:
// * runtime_error is thrown if the required file is not valid).
// */
FAdaptiveUnifKernel(const int inTreeHeight, const FReal inBoxWidth,
const FPoint& inBoxCenter, const MatrixKernelClass *const inMatrixKernel) : KernelBaseClass(inTreeHeight, inBoxWidth, inBoxCenter, inMatrixKernel), M2LHandler(inMatrixKernel, inTreeHeight, inBoxWidth), MatrixKernel(inMatrixKernel)
const FPoint& inBoxCenter, const MatrixKernelClass *const inMatrixKernel) : KernelBaseClass(inTreeHeight, inBoxWidth, inBoxCenter, inMatrixKernel)/*, M2LHandler(inMatrixKernel, inTreeHeight, inBoxWidth)*/, MatrixKernel(inMatrixKernel)
{}
// /** Copy constructor */
FAdaptiveUnifKernel(const FAdaptiveUnifKernel& other)
: KernelBaseClass(other), M2LHandler(other.M2LHandler), MatrixKernel(other.MatrixKernel)
: KernelBaseClass(other)/*, M2LHandler(other.M2LHandler)*/, MatrixKernel(other.MatrixKernel)
{ }
//
......@@ -435,7 +440,6 @@ public:
// for(int idxPart = 0 ; idxPart < target->getNbParticles() ; ++idxPart){
// particlesAttributes[idxPart] += sources->getNbParticles();
// }
std::cout << "AdapP2P" << std::endl;
}
bool preferP2M(const ContainerClass* const particles) override {
......
......@@ -135,6 +135,7 @@ protected:
// Manage border limits
if(arrayIndex == this->leftLeafIndex && arrayIndex == this->rightLeafIndex){
this->rightLeafIndex = -1;
// only one cells, return true
return true;
}
......
This diff is collapsed.
......@@ -72,6 +72,27 @@ protected:
BoxCorner.getZ() + (FReal(Coordinate.getZ()) + FReal(.5)) * BoxWidthLeaf);
}
/**
* @brief Return the position of the center of a cell from its tree
* coordinate
* @param FTreeCoordinate
* @param inLevel the current level of Cell
*/
FPoint getCellCenter(const FTreeCoordinate coordinate, int inLevel)
{
//Set the boxes width needed
FReal widthAtCurrentLevel = BoxWidthLeaf*FReal(1 << (TreeHeight-(inLevel+1)));
FReal widthAtCurrentLevelDiv2 = widthAtCurrentLevel/FReal(2.);
//Set the center real coordinates from box corner and widths.
FReal X = BoxCorner.getX() + FReal(coordinate.getX())*widthAtCurrentLevel + widthAtCurrentLevelDiv2;
FReal Y = BoxCorner.getY() + FReal(coordinate.getY())*widthAtCurrentLevel + widthAtCurrentLevelDiv2;
FReal Z = BoxCorner.getZ() + FReal(coordinate.getZ())*widthAtCurrentLevel + widthAtCurrentLevelDiv2;
return FPoint(X,Y,Z);
}
public:
/**
* The constructor initializes all constant attributes and it reads the
......
......@@ -357,7 +357,7 @@ template <int ORDER, class MatrixKernelClass>
unsigned int ComputeAndCompress(const MatrixKernelClass *const MatrixKernel,
const FReal CellWidth,
const FReal CellWidthExtension,
const FReal epsilon,
const FReal /*epsilon*/,
FReal* &U,
FReal** &C,
FReal* &B)
......
......@@ -83,7 +83,7 @@ struct FInterpMatrixKernelR : FInterpAbstractMatrixKernel
FInterpMatrixKernelR() {}
// copy ctor
FInterpMatrixKernelR(const FInterpMatrixKernelR& other)
FInterpMatrixKernelR(const FInterpMatrixKernelR& /*other*/)
{}
// returns position in reduced storage
......@@ -151,7 +151,7 @@ struct FInterpMatrixKernelRH :FInterpMatrixKernelR{
// copy ctor
FInterpMatrixKernelRH(const FInterpMatrixKernelRH& other)
: LX(other.LX), LY(other.LY), LZ(other.LZ)
:FInterpMatrixKernelR(other), LX(other.LX), LY(other.LY), LZ(other.LZ)
{}
// evaluate interaction
......@@ -218,7 +218,7 @@ struct FInterpMatrixKernelRR : FInterpAbstractMatrixKernel
FInterpMatrixKernelRR() {}
// copy ctor
FInterpMatrixKernelRR(const FInterpMatrixKernelRR& other)
FInterpMatrixKernelRR(const FInterpMatrixKernelRR& /*other*/)
{}
// returns position in reduced storage
......@@ -287,7 +287,7 @@ struct FInterpMatrixKernelLJ : FInterpAbstractMatrixKernel
FInterpMatrixKernelLJ() {}
// copy ctor
FInterpMatrixKernelLJ(const FInterpMatrixKernelLJ& other)
FInterpMatrixKernelLJ(const FInterpMatrixKernelLJ& /*other*/)
{}
// returns position in reduced storage
......
......@@ -653,6 +653,7 @@ class FRotationKernel : public FAbstractKernels<CellClass,ContainerClass> {
*/
// theta should be between [0;pi] as the inclinaison angle
void DlmkBuild(FReal dlmk[P+1][P2+1][P2+1], const FReal inTheta) const {
FAssertLF(0 <= inTheta && inTheta < FMath::FTwoPi);
// To have constants for very used values
const FReal F0 = FReal(0.0);
const FReal F1 = FReal(1.0);
......
......@@ -98,7 +98,7 @@ public:
{
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// 1) apply Sy
FBlas::scal(AbstractBaseClass::nnodes, FReal(0.), ParentCell->getMultipole(idxRhs));
//FBlas::scal(AbstractBaseClass::nnodes, FReal(0.), ParentCell->getMultipole(idxRhs));
for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){
if (ChildCells[ChildIndex]){
AbstractBaseClass::Interpolator->applyM2M(ChildIndex, ChildCells[ChildIndex]->getMultipole(idxRhs),
......
......@@ -179,6 +179,13 @@ public:
return *this;
}
/** Mul a complexe by another "c*=c2" */
FComplexe& operator/=(const FReal& real){
this->complex[0] /= real;
this->complex[1] /= real;
return *this;
}
/** Test if a complex is not a number */
bool isNan() const {
return FMath::IsNan(complex[1]) || FMath::IsNan(complex[0]);
......
......@@ -262,7 +262,7 @@ struct FMath{
FReal getmax() const{
return max;
}
FReal getNbElements() const{
int getNbElements() const{
return nbElements;
}
void setNbElements(const int & n) {
......
......@@ -263,31 +263,31 @@ public:
// Mpi Types meta function
////////////////////////////////////////////////////////////
static const MPI_Datatype GetType(const long long&){
static MPI_Datatype GetType(const long long&){
return MPI_LONG_LONG;
}
static const MPI_Datatype GetType(const long int&){
static MPI_Datatype GetType(const long int&){
return MPI_LONG;
}
static const MPI_Datatype GetType(const double&){
static MPI_Datatype GetType(const double&){
return MPI_DOUBLE;
}
static const MPI_Datatype GetType(const float&){
static MPI_Datatype GetType(const float&){
return MPI_FLOAT;
}
static const MPI_Datatype GetType(const int&){
static MPI_Datatype GetType(const int&){
return MPI_INT;
}
static const MPI_Datatype GetType(const char&){
static MPI_Datatype GetType(const char&){
return MPI_CHAR;
}
static const MPI_Datatype GetType(const FComplexe& a){
static MPI_Datatype GetType(const FComplexe& a){
MPI_Datatype FMpiComplexe;
MPI_Type_contiguous(2, GetType(a.getReal()) , &FMpiComplexe);
return FMpiComplexe;
......
// ===================================================================================
// 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 FOMPBARRIER_HPP
#define FOMPBARRIER_HPP
#include <omp.h>
#include <climits>
/** This function is a custom omp barrier
* Because openmp give only a global barrier we need
* to be ablo to peform a barrier operation between a group
* of threads only.
*/
class FOmpBarrier {
private:
int nbThreads; //<The number of threads for this barrier
int currentNbThread; //<The current number of threads waiting
bool sense; //<Direct barrier feedback protection
omp_lock_t mutex; //<To have an atomic int
FOmpBarrier(FOmpBarrier&){}
FOmpBarrier& operator=(FOmpBarrier&){return *this;}
public:
/** Constructor with the number of threads */
FOmpBarrier(const int inNbThreads = INT_MAX)
: nbThreads(inNbThreads), currentNbThread(0), sense(false) {
omp_init_lock( &mutex );
}
/** Destructor, release the omp lock */
~FOmpBarrier(){
omp_destroy_lock( &mutex );
}
/** Perform a barrier */
void wait(){
const bool mySense = sense;
omp_set_lock( &mutex );
const int nbThreadsArrived = (++currentNbThread);
omp_unset_lock( &mutex );
if(nbThreadsArrived == nbThreads) {
currentNbThread = 0;
sense = !sense;
#pragma omp flush
}
else {
volatile const bool* const ptSense = &sense;
while( (*ptSense) == mySense){
}
}
}
/** Change the number of threads */
void setNbThreads(const int inNbThread){
omp_set_lock( &mutex );
nbThreads = inNbThread;
omp_unset_lock( &mutex );
}
};
#endif // FOMPBARRIER_HPP
......@@ -27,8 +27,6 @@
#include "FMemUtils.hpp"
#include "FTrace.hpp"
#include "FOmpBarrier.hpp"
/** This class is parallel quick sort
* It hold a mpi version
* + 2 openmp versions (one on tasks and the other like mpi)
......
......@@ -177,6 +177,7 @@ class FQuickSortMpi : public FQuickSort< SortType, CompareType, IndexType> {
}
// Wait to complete
FMpi::Assert( MPI_Waitall(whatToRecvFromWho.size(), requests, MPI_STATUSES_IGNORE), __LINE__ );
////FLOG( FLog::Controller << currentComm.processId() << "] Recv Done \n"; )
// Copy to ouput variables
(*inPartRecv) = recvBuffer;
(*inNbElementsRecv) = totalToRecv;
......@@ -200,6 +201,7 @@ class FQuickSortMpi : public FQuickSort< SortType, CompareType, IndexType> {
}
// Wait to complete
FMpi::Assert( MPI_Waitall(whatToSendToWho.size(), requests, MPI_STATUSES_IGNORE), __LINE__ );
////FLOG( FLog::Controller << currentComm.processId() << "] Send Done \n"; )
}
static CompareType SelectPivot(const SortType workingArray[], const IndexType currentSize, const FMpi::FComm& currentComm, bool* shouldStop){
......@@ -208,12 +210,16 @@ class FQuickSortMpi : public FQuickSort< SortType, CompareType, IndexType> {
NO_VALUES,
AVERAGE_2
};
// We need to know the max value to ensure that the pivot will be different
CompareType maxFoundValue = CompareType(workingArray[0]);
// Check if all the same
bool allTheSame = true;
for(int idx = 1 ; idx < currentSize && allTheSame; ++idx){
if(workingArray[0] != workingArray[idx]){
allTheSame = false;
}
// Keep the max
maxFoundValue = FMath::Max(maxFoundValue , CompareType(workingArray[idx]));
}
// Check if empty
const bool noValues = (currentSize == 0);
......@@ -221,6 +227,11 @@ class FQuickSortMpi : public FQuickSort< SortType, CompareType, IndexType> {
CompareType localPivot = CompareType(0);
if(!noValues){
localPivot = (CompareType(workingArray[currentSize/3])+CompareType(workingArray[(2*currentSize)/3]))/2;
// The pivot must be different (to ensure that the partition will return two parts)
if( localPivot == maxFoundValue && !allTheSame){
////FLOG( FLog::Controller << currentComm.processId() << "] Pivot " << localPivot << " is equal max and allTheSame equal " << allTheSame << "\n"; )
localPivot -= 1;
}
}
////FLOG( FLog::Controller << currentComm.processId() << "] localPivot = " << localPivot << "\n" );
......@@ -350,7 +361,9 @@ public:
workingArray = fullLowerPart;
currentSize = fullNbLowerElementsRecv;
// Reduce working group
////FLOG( FLog::Controller << currentComm.processId() << "] Reduce group to " << 0 << " / " << procInTheMiddle << "\n"; )
currentComm.groupReduce( 0, procInTheMiddle);
////FLOG( FLog::Controller << currentComm.processId() << "] Done\n" );
}
else {
// I am in the group of the greater elements
......@@ -370,10 +383,13 @@ public:
workingArray = fullGreaterPart;
currentSize = fullNbGreaterElementsRecv;
// Reduce working group
////FLOG( FLog::Controller << currentComm.processId() << "] Reduce group to " << procInTheMiddle + 1 << " / " << currentNbProcs - 1 << "\n"; )
currentComm.groupReduce( procInTheMiddle + 1, currentNbProcs - 1);
////FLOG( FLog::Controller << currentComm.processId() << "] Done\n"; )
}
}
////FLOG( FLog::Controller << currentComm.processId() << "] Sequential sort\n"; )
// Finish by a local sort
FQuickSort< SortType, CompareType, IndexType>::QsOmp(workingArray, currentSize);
(*outputSize) = currentSize;
......
......@@ -98,7 +98,7 @@ int main(int argc, char ** argv){
//////////////////////////////////////////////////////////////