Commit c92b4ef6 authored by PIACIBELLO Cyrille's avatar PIACIBELLO Cyrille
Browse files

Formating the output

parent bddfa4ac
......@@ -77,153 +77,169 @@ void usage() {
// Simply create particles and try the kernels
int main(int argc, char* argv[])
{
const std::string defaultFile(/*SCALFMMDataPath+*/"../Data/test20k.fma");
const std::string filename = FParameters::getStr(argc,argv,"-f", defaultFile.c_str());
const unsigned int TreeHeight = FParameters::getValue(argc, argv, "-depth", 5);
const unsigned int SubTreeHeight = FParameters::getValue(argc, argv, "-subdepth", 2);
const unsigned int NbThreads = FParameters::getValue(argc, argv, "-t", 1);
if(FParameters::existParameter(argc, argv, "-h")||FParameters::existParameter(argc, argv, "-help")){
usage() ;
exit(-1);
}
const std::string defaultFile(/*SCALFMMDataPath+*/"../Data/test20k.fma");
const std::string filename = FParameters::getStr(argc,argv,"-f", defaultFile.c_str());
const unsigned int TreeHeight = FParameters::getValue(argc, argv, "-depth", 5);
const unsigned int SubTreeHeight = FParameters::getValue(argc, argv, "-subdepth", 2);
const unsigned int NbThreads = FParameters::getValue(argc, argv, "-t", 1);
if(FParameters::existParameter(argc, argv, "-h")||FParameters::existParameter(argc, argv, "-help")){
usage() ;
exit(-1);
}
#ifdef _OPENMP
omp_set_num_threads(NbThreads);
std::cout << "\n>> Using " << omp_get_max_threads() << " threads.\n" << std::endl;
omp_set_num_threads(NbThreads);
std::cout << "\n>> Using " << omp_get_max_threads() << " threads.\n" << std::endl;
#else
std::cout << "\n>> Sequential version.\n" << std::endl;
std::cout << "\n>> Sequential version.\n" << std::endl;
#endif
//
std::cout << "Parameters "<< std::endl
<< " Octree Depth "<< TreeHeight <<std::endl
<< " SubOctree depth "<< SubTreeHeight <<std::endl
<< " Input file name: " <<filename <<std::endl
<< " Thread number: " << NbThreads <<std::endl
<<std::endl;
//init values for MPI
FMpi app(argc,argv);
//
// init timer
FTic time;
FMpiFmaLoader loader(filename,app.global());
if(!loader.isOpen()) throw std::runtime_error("Particle file couldn't be opened!") ;
////////////////////////////////////////////////////////////////////
// begin spherical kernel
// accuracy
const unsigned int P = 22;
// typedefs
typedef FP2PParticleContainerIndexed<> ContainerClass;
typedef FSimpleLeaf< ContainerClass > LeafClass;
typedef FRotationCell<P> CellClass;
typedef FOctree<CellClass,ContainerClass,LeafClass> OctreeClass;
//
typedef FRotationKernel< CellClass, ContainerClass , P> KernelClass;
//
typedef FFmmAlgorithmThreadProc<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClassProc;
// init oct-tree
OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
{ // -----------------------------------------------------
std::cout << "Creating & Inserting " << loader.getNumberOfParticles()
<< " particles ..." << std::endl;
std::cout << "\tHeight : " << TreeHeight << " \t sub-height : " << SubTreeHeight << std::endl;
time.tic();
//
struct TestParticle{
FPoint position;
FReal physicalValue;
const FPoint& getPosition(){
return position;
}
};
TestParticle* particles = new TestParticle[loader.getNumberOfParticles()];
memset(particles, 0, (unsigned int) (sizeof(TestParticle) * loader.getNumberOfParticles()));
//
std::cout << "Parameters "<< std::endl
<< " Octree Depth "<< TreeHeight <<std::endl
<< " SubOctree depth "<< SubTreeHeight <<std::endl
<< " Input file name: " <<filename <<std::endl
<< " Thread number: " << NbThreads <<std::endl
<<std::endl;
//init values for MPI
FMpi app(argc,argv);
//
// init timer
FTic time;
FMpiFmaLoader loader(filename,app.global(),true);
if(!loader.isOpen()) throw std::runtime_error("Particle file couldn't be opened!") ;
////////////////////////////////////////////////////////////////////
// begin spherical kernel
// accuracy
const unsigned int P = 22;
// typedefs
typedef FP2PParticleContainerIndexed<> ContainerClass;
typedef FSimpleLeaf< ContainerClass > LeafClass;
typedef FRotationCell<P> CellClass;
typedef FOctree<CellClass,ContainerClass,LeafClass> OctreeClass;
//
typedef FRotationKernel< CellClass, ContainerClass , P> KernelClass;
//
typedef FFmmAlgorithmThreadProc<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClassProc;
// init oct-tree
OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
{ // -----------------------------------------------------
std::cout << "Creating & Inserting " << loader.getNumberOfParticles()
<< " particles ..." << std::endl;
std::cout << "\tHeight : " << TreeHeight << " \t sub-height : " << SubTreeHeight << std::endl;
time.tic();
//
struct TestParticle{
int indexInFile;
FPoint position;
FReal physicalValue;
const FPoint& getPosition(){
return position;
}
};
TestParticle* particles = new TestParticle[loader.getNumberOfParticles()];
memset(particles, 0, (unsigned int) (sizeof(TestParticle) * loader.getNumberOfParticles()));
for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
// Read particles from file
loader.fillParticle(&particles[idxPart].position,&particles[idxPart].physicalValue);
}
FVector<TestParticle> finalParticles;
FLeafBalance balancer;
FMpiTreeBuilder< TestParticle >::ArrayToTree(app.global(), particles, loader.getNumberOfParticles(),
//idx (in file) of the first part that will be used by this proc.
int idxStart = loader.getStart();
for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
//Storage of the index (in the original file) of each part.
particles[idxPart].indexInFile = idxPart + idxStart;
// Read particles from file
loader.fillParticle(&particles[idxPart].position,&particles[idxPart].physicalValue);
}
FVector<TestParticle> finalParticles;
FLeafBalance balancer;
FMpiTreeBuilder< TestParticle >::ArrayToTree(app.global(), particles, loader.getNumberOfParticles(),
tree.getBoxCenter(),
tree.getBoxWidth(),
tree.getHeight(), &finalParticles,&balancer);
for(int idx = 0 ; idx < finalParticles.getSize(); ++idx){
tree.insert(finalParticles[idx].position,idx,finalParticles[idx].physicalValue);
}
delete[] particles;
for(int idx = 0 ; idx < finalParticles.getSize(); ++idx){
tree.insert(finalParticles[idx].position,finalParticles[idx].indexInFile,finalParticles[idx].physicalValue);
}
delete[] particles;
time.tac();
std::cout << "Done " << "(@Creating and Inserting Particles = "
<< time.elapsed() << "s)." << std::endl;
} // -----------------------------------------------------
time.tac();
std::cout << "Done " << "(@Creating and Inserting Particles = "
<< time.elapsed() << "s)." << std::endl;
} // -----------------------------------------------------
{ // -----------------------------------------------------
std::cout << "\nRotation harmonic Spherical FMM Proc (P="<< P << ") ... " << std::endl;
{ // -----------------------------------------------------
std::cout << "\nRotation harmonic Spherical FMM Proc (P="<< P << ") ... " << std::endl;
time.tic();
//
// Here we use a pointer due to the limited size of the stack
//
KernelClass *kernels = new KernelClass(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
//
FmmClassProc algorithm(app.global(),&tree, kernels);
//
algorithm.execute(); // Here the call of the FMM algorithm
//
time.tac();
std::cout << "Done " << "(@Algorithm = " << time.elapsed() << "s)." << std::endl;
}
// -----------------------------------------------------
//
// Some output
//
//
{ // -----------------------------------------------------
long int N1=0, N2= loader.getNumberOfParticles()/4, N3= (loader.getNumberOfParticles() -1)/2; ;
FReal energy =0.0 ;
//
// Loop over all leaves
//
std::cout <<std::endl<<" &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& "<<std::endl;
std::cout << std::scientific;
std::cout.precision(10) ;
tree.forEachLeaf([&](LeafClass* leaf){
const FReal*const potentials = leaf->getTargets()->getPotentials();
const FReal*const forcesX = leaf->getTargets()->getForcesX();
const FReal*const forcesY = leaf->getTargets()->getForcesY();
const FReal*const forcesZ = leaf->getTargets()->getForcesZ();
const int nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
time.tic();
//
// Here we use a pointer due to the limited size of the stack
//
KernelClass *kernels = new KernelClass(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
//
FmmClassProc algorithm(app.global(),&tree, kernels);
//
algorithm.execute(); // Here the call of the FMM algorithm
//
time.tac();
std::cout << "Done " << "(@Algorithm = " << time.elapsed() << "s)." << std::endl;
}
// -----------------------------------------------------
//
// Some output
//
//
{ // -----------------------------------------------------
long int N1=0, N2= loader.getTotalNumberOfParticles()/2, N3= (loader.getTotalNumberOfParticles()-1); ;
FReal energy =0.0 ;
//
// Loop over all leaves
//
std::cout <<std::endl<<" &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& "<<std::endl;
std::cout << std::scientific;
std::cout.precision(10) ;
tree.forEachLeaf([&](LeafClass* leaf){
const FReal*const posX = leaf->getTargets()->getPositions()[0];
const FReal*const posY = leaf->getTargets()->getPositions()[1];
const FReal*const posZ = leaf->getTargets()->getPositions()[2];
const FReal*const potentials = leaf->getTargets()->getPotentials();
const FReal*const forcesX = leaf->getTargets()->getForcesX();
const FReal*const forcesY = leaf->getTargets()->getForcesY();
const FReal*const forcesZ = leaf->getTargets()->getForcesZ();
const int nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
const FVector<int>& indexes = leaf->getTargets()->getIndexes();
const FVector<int>& indexes = leaf->getTargets()->getIndexes();
for(int idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
const int indexPartOrig = indexes[idxPart];
if ((indexPartOrig == N1) || (indexPartOrig == N2) || (indexPartOrig == N3) ) {
std::cout << "Proc "<< app.global().processId() << " Index "<< indexPartOrig <<" potential " << potentials[idxPart]
<< " Forces: " << forcesX[idxPart] << " " << forcesY[idxPart] << " "<< forcesZ[idxPart] <<std::endl;
}
energy += potentials[idxPart]*physicalValues[idxPart] ;
}
});
std::cout <<std::endl<<"Energy: "<< energy<<std::endl;
std::cout <<std::endl<<" &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& "<<std::endl<<std::endl;
for(int idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
const int indexPartOrig = indexes[idxPart];
if ((indexPartOrig == N1) || (indexPartOrig == N2) || (indexPartOrig == N3) ) {
std::cout << "Proc "<< app.global().processId() << " Index "<< indexPartOrig << " potential " << potentials[idxPart]
<< " Pos "<<posX[idxPart]<<" "<<posY[idxPart]<<" "<<posZ[idxPart]
<< " Forces: " << forcesX[idxPart] << " " << forcesY[idxPart] << " "<< forcesZ[idxPart] <<std::endl;
}
energy += potentials[idxPart]*physicalValues[idxPart] ;
}
// -----------------------------------------------------
return 0;
});
FReal gloEnergy = app.global().reduceSum(energy);
if(0 == app.global().processId()){
std::cout <<std::endl << "Proc "<< app.global().processId() << " Energy: "<< gloEnergy <<std::endl;
std::cout <<std::endl <<" &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& "<<std::endl<<std::endl;
}
}
// -----------------------------------------------------
return 0;
}
This diff is collapsed.
#include <cstdlib>
#include <string.h>
#include <stdexcept>
#include "../../Src/Utils/FMpi.hpp"
#include "../../Src/Containers/FVector.hpp"
#include "../../Src/Files/FFmaGenericLoader.hpp"
#include "../../Src/Files/FMpiFmaLoader.hpp"
#include "../../Src/BalanceTree/FLeafBalance.hpp"
#include "../../Src/Containers/FTreeCoordinate.hpp"
#include "../../Src/Utils/FQuickSortMpi.hpp"
#include "../../Src/Files/FMpiTreeBuilder.hpp"
#include "../../Src/Core/FCoreCommon.hpp"
#include "../../Src/Utils/FPoint.hpp"
#include "../../Src/Utils/FMath.hpp"
#include "../../Src/Utils/FParameters.hpp"
/* Struct Test particle : classic particle but with index In file,
* which is the number of line where the part was found. It is a true
* constant.
*/
struct TestParticle{
MortonIndex index;
FSize indexInFile;
FPoint position;
FReal physicalValue;
FPoint& getPosition(){
return position;
}
};
/** This method has been tacken from the octree
* it computes a tree coordinate (x or y or z) from real position
*/
static int getTreeCoordinate(const FReal inRelativePosition, const FReal boxWidthAtLeafLevel) {
const FReal indexFReal = inRelativePosition / boxWidthAtLeafLevel;
const int index = int(FMath::dfloor(indexFReal));
if( index && FMath::LookEqual(inRelativePosition, boxWidthAtLeafLevel * FReal(index) ) ){
return index - 1;
}
return index;
}
/**
* This function compare IndexedParticles that are in the same leaf by
* sorting them along X,Y and then Z axis. The aim is to get the same
* sorting output for each leaf.
*/
int compareLeafI(const void *a,const void *b){
FMpiTreeBuilder<TestParticle>::IndexedParticle partA = *(reinterpret_cast<const FMpiTreeBuilder<TestParticle>::IndexedParticle*>(a));
FMpiTreeBuilder<TestParticle>::IndexedParticle partB = *(reinterpret_cast<const FMpiTreeBuilder<TestParticle>::IndexedParticle*>(b));
if(partA.particle.getPosition().getX() < partB.particle.getPosition().getX()){
return -1;
}
else{
if(partA.particle.getPosition().getX() > partB.particle.getPosition().getX()){
return 1;
}
else{
if(partA.particle.getPosition().getY() < partB.particle.getPosition().getY()){
return -1;
}
else{
if(partA.particle.getPosition().getY() > partB.particle.getPosition().getY()){
return 1;
}
else{
if(partA.particle.getPosition().getZ() < partB.particle.getPosition().getZ()){
return -1;
}
else{
if(partA.particle.getPosition().getZ() > partB.particle.getPosition().getZ()){
return 1;
}
else{
return 0;
}
}
}
}
}
}
}
/**
* Same function as above but with TestParticle instead of
* IndexedParticles (pretty much same thing).
*/
int compareLeaf(const void *a,const void *b){
TestParticle partA = *(reinterpret_cast<const TestParticle*>(a));
TestParticle partB = *(reinterpret_cast<const TestParticle*>(b));
if(partA.getPosition().getX() < partB.getPosition().getX()){
return -1;
}
else{
if(partA.getPosition().getX() > partB.getPosition().getX()){
return 1;
}
else{
if(partA.getPosition().getY() < partB.getPosition().getY()){
return -1;
}
else{
if(partA.getPosition().getY() > partB.getPosition().getY()){
return 1;
}
else{
if(partA.getPosition().getZ() < partB.getPosition().getZ()){
return -1;
}
else{
if(partA.getPosition().getZ() > partB.getPosition().getZ()){
return 1;
}
else{
return 0;
printf("Should not happen\n");
}
}
}
}
}
}
}
/**
* Compare two articles by MortonIndex, if same morton index, call to
* compareLeaf
*/
int cmpTestParticles(const void* a,const void* b){
const TestParticle partA = *(reinterpret_cast<const TestParticle*>(a));
const TestParticle partB = *(reinterpret_cast<const TestParticle*>(b));
if(partA.index < partB.index){
return -1;
}
else{
if(partA.index > partB.index){
return 1;
}
else{
return compareLeaf(a,b);
}
}
}
/** This test aimed to validate the Tree Builder algorithms
*
* Each process runs a sort (from std::qsort) and then a parallel
* sort and compare its particles with the global sort.
*
* Then, each process runs a MergeLeaves, and compare again the output
* with the one from sstd::qsort.
*
* Finally, each process runs EqualizeAndFillTree, and compare again
* the output with the one from std::sort.
*
*/
int main(int argc,char* argv[]){
//File to sort
const char* const filename = FParameters::getStr(argc,argv,"-f", "../Data/test20k.bin.fma.double");
const unsigned int TreeHeight = FParameters::getValue(argc, argv, "-depth", 5);
//First part is Sequential
//Since there is no MPI Loader with ascii datas, the file need to be in binary format.
FFmaGenericLoader loaderSeq(filename,true);
if(!loaderSeq.isOpen()) throw std::runtime_error("Particle file couldn't be opened!") ;
//Get the needed informations
FReal boxWidth = loaderSeq.getBoxWidth();
FReal boxWidthAtLeafLevel = boxWidth/FReal(1 << (TreeHeight - 1));
FPoint centerOfBox = loaderSeq.getCenterOfBox();
FPoint boxCorner = centerOfBox - boxWidth/2;
FTreeCoordinate host;
//Copy from the file to the array that will be sorted
FSize nbOfParticles = loaderSeq.getNumberOfParticles();
struct TestParticle* arrayOfParticles = new struct TestParticle[nbOfParticles];
for(FSize idxParts=0 ; idxParts<nbOfParticles ; ++idxParts){
//Fill automatically position AND physicalValue attributes
loaderSeq.fillParticle(&(arrayOfParticles[idxParts].position),&(arrayOfParticles[idxParts].physicalValue));
//We store the index in the file
arrayOfParticles[idxParts].indexInFile = idxParts;
//Build temporary TreeCoordinate
host.setX( getTreeCoordinate( arrayOfParticles[idxParts].getPosition().getX() - boxCorner.getX(), boxWidthAtLeafLevel ));
host.setY( getTreeCoordinate( arrayOfParticles[idxParts].getPosition().getY() - boxCorner.getY(), boxWidthAtLeafLevel ));
host.setZ( getTreeCoordinate( arrayOfParticles[idxParts].getPosition().getZ() - boxCorner.getZ(), boxWidthAtLeafLevel ));
//Set Morton index from Tree Coordinate
arrayOfParticles[idxParts].index = host.getMortonIndex(TreeHeight - 1);
}
//Save the original array
struct TestParticle * originalArray = new TestParticle[nbOfParticles];
memcpy(originalArray,arrayOfParticles,sizeof(struct TestParticle)*nbOfParticles);
//Sort the array
qsort((void*) arrayOfParticles,nbOfParticles,sizeof(TestParticle),cmpTestParticles);
//Start of the parallel part :
FMpi app(argc,argv);
FSize outputSize;
//Refer to ChebyshevInterpolationAlgorithmProc to know how to FMpiFmaLoader + index
FMpiFmaLoader loader(filename,app.global());
if(!loader.isOpen()) throw std::runtime_error("Particle file couldn't be opened!") ;
//Now, we sort again the particles with MPI QuickSort
int idxStart = loader.getStart();
FMpiTreeBuilder<TestParticle>::IndexedParticle * arrayToBeSorted = new FMpiTreeBuilder<TestParticle>::IndexedParticle[loader.getNumberOfParticles()];
//Copy the TestParticles into an array of indexedParticle
for(int i=0 ; i<loader.getNumberOfParticles() ; ++i){
arrayToBeSorted[i].particle = originalArray[i+idxStart];
arrayToBeSorted[i].index = arrayToBeSorted[i].particle.index;
}
FMpiTreeBuilder<TestParticle>::IndexedParticle* outputArray = 0;
FQuickSortMpi<FMpiTreeBuilder<TestParticle>::IndexedParticle,MortonIndex,FSize>::QsMpi(arrayToBeSorted,loader.getNumberOfParticles(),outputArray,outputSize,app.global());
//Sum the outputSize of every body for knowing where to start inside the sorted array
FSize starter = 0;
//We use a prefix sum
MPI_Exscan(&outputSize,&starter,1,MPI_LONG_LONG_INT,MPI_SUM,app.global().getComm());
//We sort the output array relatvely to position inside leafs
FSize inc = 0;
FMpiTreeBuilder<TestParticle>::IndexedParticle * saveForSort = outputArray;
int nbOfPartsInLeaf = 0;
while(inc < outputSize){
while(outputArray[inc].index == saveForSort->index){
inc++;
nbOfPartsInLeaf++;
}
qsort((void*)saveForSort,nbOfPartsInLeaf,sizeof(FMpiTreeBuilder<TestParticle>::IndexedParticle),compareLeafI);
nbOfPartsInLeaf = 0;
saveForSort = &outputArray[inc];
}
bool resultQsMpi = true; //passed..
//Test
for(int i=0 ; i<outputSize ; ++i){
if(outputArray[i].particle.indexInFile != arrayOfParticles[i+starter].indexInFile){
printf("QsMPI :: Proc %d i=%d Particles file : [%lld,%lld] with Mpi Morton : %lld != %lld\n",app.global().processId(),i,
arrayOfParticles[i+idxStart].indexInFile,
outputArray[i].particle.indexInFile,
outputArray[i].index,
arrayOfParticles[i+idxStart].index);
resultQsMpi = false;
}
}
if(!resultQsMpi){
printf("Proc :: %d ====> QsMpi didn't make it !!\n",app.global().processId());
}
else{
printf("Proc :: %d ====> QsMpi successfull !!\n",app.global().processId());
}
//Test MergeLeaves
bool resultMergeLeaves= true;
//inputs needed
TestParticle * leavesArray = 0;
FSize * leavesIndices = 0;
FSize leaveSize = 0;
FMpiTreeBuilder<TestParticle>::testMergeLeaves(app.global(),outputArray,&outputSize,&leavesIndices,&leavesArray,&leaveSize);
//Compare again the results with the output of std::qsort