testInsert.cpp 10.6 KB
Newer Older
1
// ===================================================================================
2 3 4 5
// Copyright ScalFmm 2016 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.
6
//
7
// This software is governed by the CeCILL-C and LGPL licenses and
8
// abiding by the rules of distribution of free software.
9 10 11
// An extension to the license is given to allow static linking of scalfmm
// inside a proprietary application (no matter its license).
// See the main license file for more details.
12 13 14 15
//
// 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
16 17 18
// GNU General Public and CeCILL-C Licenses for more details.
// "http://www.cecill.info".
// "http://www.gnu.org/licenses".
19 20
// ===================================================================================

21 22 23
// Keep in private GIT
// @SCALFMM_PRIVATE

24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55

#include <iostream>

#include <cstdlib>
#include <string>
#include <stdexcept>
#include <algorithm>
#include <vector>

#include <cstdio>

#include "ScalFmmConfig.h"

#include "Files/FFmaGenericLoader.hpp"

#include "Components/FSimpleLeaf.hpp"
#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"

#include "../../Src/Kernels/P2P/FP2PParticleContainer.hpp"

#include "../../Src/Kernels/Rotation/FRotationKernel.hpp"
#include "../../Src/Kernels/Rotation/FRotationCell.hpp"

#include "../Src/BalanceTree/FLeafBalance.hpp"
#include "../Src/Containers/FTreeCoordinate.hpp"

#include "Utils/FTic.hpp"
#include "Utils/FQuickSort.hpp"
#include "Utils/FParameters.hpp"
#include "../../Src/Utils/FParameterNames.hpp"

#include "Containers/FOctree.hpp"
56
#include "Containers/FCoordinateComputer.hpp"
57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73

#ifdef _OPENMP
#include "Core/FFmmAlgorithmThread.hpp"
#else
#include "Core/FFmmAlgorithm.hpp"
#endif

#include "Utils/FTemplate.hpp"


/**
 * This program build a tree and insert the parts inside.
 * Time needed for the insert is outputed
 *
 */
int main(int argc, char** argv){

74
    typedef double FReal;
75
    struct TestParticle{
76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103
        MortonIndex index;
        FSize indexInFile;
        FPoint<FReal> position;
        FReal physicalValue;
        const FPoint<FReal>& getPosition()const{
            return position;
        }
        TestParticle& operator=(const TestParticle& other){
            index=other.index;
            indexInFile=other.indexInFile;
            position=other.position;
            physicalValue=other.physicalValue;
            return *this;
        }
        bool operator<=(const TestParticle& rhs)const{
            if(rhs.index < this->index){return false;}
            else{
                if(rhs.index > this->index){return true;}
                else{
                    if(rhs.indexInFile == this->indexInFile){
                        return true;
                    }
                    else {
                        return rhs.indexInFile> this->indexInFile ;
                    }
                }
            }
        }
104 105 106
    };


107
    typedef double FReal;
108 109
    static const int P = 9;

110
    typedef FRotationCell<FReal,P>               CellClass;
111
    typedef FP2PParticleContainer<FReal>          ContainerClass;
112

113 114
    typedef FSimpleLeaf<FReal, ContainerClass >                     LeafClass;
    typedef FOctree<FReal, CellClass, ContainerClass , LeafClass >  OctreeClass;
115 116 117 118 119 120 121 122 123 124 125 126

    const int NbLevels = FParameters::getValue(argc,argv,FParameterDefinitions::OctreeHeight.options, 5);
    const int SizeSubLevels = FParameters::getValue(argc,argv,FParameterDefinitions::OctreeSubHeight.options, 3);

    const char* const filename = FParameters::getStr(argc,argv,FParameterDefinitions::InputFile.options, "../Data/test20k.fma");
    const unsigned int NbThreads = FParameters::getValue(argc, argv, FParameterDefinitions::NbThreads.options, omp_get_max_threads());

    omp_set_num_threads(NbThreads);
    std::cout << "Using " << omp_get_max_threads() <<" threads" << std::endl;
    std::cout << "Opening : " << filename << "\n";


127
    FFmaGenericLoader<FReal> loaderRef(filename);
128
    if(!loaderRef.isOpen()){
129 130
        std::cout << "LoaderRef Error, " << filename << " is missing\n";
        return 1;
131 132 133 134
    }
    FTic regInsert;
    // -----------------------------------------------------
    {
135
        OctreeClass treeRef(NbLevels, SizeSubLevels, loaderRef.getBoxWidth(), loaderRef.getCenterOfBox());
136 137


138 139
        std::cout << "Creating & Inserting " << loaderRef.getNumberOfParticles() << " particles ..." << std::endl;
        std::cout << "\tHeight : " << NbLevels << " \t sub-height : " << SizeSubLevels << std::endl;
140

141
        regInsert.tic();
142

143 144 145 146 147 148
        for(FSize idxPart = 0 ; idxPart < loaderRef.getNumberOfParticles() ; ++idxPart){
            FPoint<FReal> particlePosition;
            FReal physicalValue;
            loaderRef.fillParticle(&particlePosition,&physicalValue);
            treeRef.insert(particlePosition, physicalValue );
        }
149

150 151
        regInsert.tac();
        std::cout << "Time needed for regular insert : " << regInsert.elapsed() << " secondes" << std::endl;
152 153 154 155 156 157 158 159 160 161 162 163
    }
    //Second solution, parts must be sorted for that

    //Timers :
    FTic sortTimer;
    FTic insertTimer;
    FTic fillTimer;
    FTic enumTimer;
    FTic counter;
    FTic leavesOffset;
    FTic leavesPtr;

164
    FFmaGenericLoader<FReal> loader(filename);
165
    if(!loader.isOpen()){
166 167
        std::cout << "Loader Error, " << filename << " is missing\n";
        return 1;
168 169 170 171 172
    }

    //Get the needed informations
    FReal boxWidth = loader.getBoxWidth();
    FReal boxWidthAtLeafLevel = boxWidth/FReal(1 << (NbLevels - 1));
173 174
    FPoint<FReal> centerOfBox = loader.getCenterOfBox();
    FPoint<FReal> boxCorner   = centerOfBox - boxWidth/2;
175 176 177 178 179 180 181 182 183
    FSize nbOfParticles = loader.getNumberOfParticles();

    //Temporary TreeCoordinate
    FTreeCoordinate host;
    TestParticle * arrayOfParts = new TestParticle[nbOfParticles];
    memset(arrayOfParts,0,sizeof(TestParticle)*nbOfParticles);

    fillTimer.tic();

184 185 186 187 188 189
    for(FSize idxPart = 0 ; idxPart < nbOfParticles ; ++idxPart){
        loader.fillParticle(&arrayOfParts[idxPart].position,&arrayOfParts[idxPart].physicalValue);
        //Build temporary TreeCoordinate
        host.setX( FCoordinateComputer::GetTreeCoordinate<FReal>( arrayOfParts[idxPart].getPosition().getX() - boxCorner.getX(), boxWidth, boxWidthAtLeafLevel, NbLevels ));
        host.setY( FCoordinateComputer::GetTreeCoordinate<FReal>( arrayOfParts[idxPart].getPosition().getY() - boxCorner.getY(), boxWidth, boxWidthAtLeafLevel, NbLevels ));
        host.setZ( FCoordinateComputer::GetTreeCoordinate<FReal>( arrayOfParts[idxPart].getPosition().getZ() - boxCorner.getZ(), boxWidth, boxWidthAtLeafLevel, NbLevels ));
190

191
        //Set Morton index from Tree Coordinate
192
        arrayOfParts[idxPart].index = host.getMortonIndex();
193
        arrayOfParts[idxPart].indexInFile = idxPart;
194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212

    }


    fillTimer.tac();
    std::cout << "Time needed for filling the array : "<< fillTimer.elapsed() << " secondes !" << std::endl;

    sortTimer.tic();
    //std::sort(arrayOfParts,&arrayOfParts[nbOfParticles-1]);
    FQuickSort<TestParticle,MortonIndex>::QsOmp(arrayOfParts,nbOfParticles);
    sortTimer.tac();

    std::cout << "Time needed for sorting the array : "<< sortTimer.elapsed() << " secondes !" << std::endl;


    //Enumerate the different leaves
    unsigned int numberOfLeaves = 0;
    enumTimer.tic();

213 214 215 216
    for(FSize idxParts = 1 ; idxParts < nbOfParticles ; ++idxParts){
        if(arrayOfParts[idxParts].index != arrayOfParts[idxParts-1].index){
            numberOfLeaves++;
        }
217 218 219 220 221 222
    }
    enumTimer.tac();
    std::cout << "Time needed for enumerate the leaves : "<< enumTimer.elapsed() << " secondes !" << std::endl;
    std::cout << "Found " << numberOfLeaves << " leaves differents." << std::endl;

    //Store the size of each subOctree
223 224
    FSize * arrayOfSizeNbLeaves = new FSize[numberOfLeaves];
    memset(arrayOfSizeNbLeaves,0,sizeof(FSize)*(numberOfLeaves));
225 226 227 228 229
    //Init
    int indexInLeafArray = 0;
    arrayOfSizeNbLeaves[0] = 1;

    leavesOffset.tic();
230 231 232 233 234 235 236
    for(FSize idxParts = 1 ; idxParts < nbOfParticles ; ++idxParts){
        if(arrayOfParts[idxParts].index == arrayOfParts[idxParts-1].index){
            arrayOfSizeNbLeaves[indexInLeafArray]++;
        }
        else{
            indexInLeafArray++;
        }
237 238 239 240 241 242 243 244 245 246
    }
    leavesOffset.tac();

    std::cout << "Time needed for setting the offset of each leaves : "<< leavesOffset.elapsed() << " secondes !" << std::endl;

    OctreeClass tree(NbLevels, SizeSubLevels, loader.getBoxWidth(), loader.getCenterOfBox());

    //Then, we create the leaves inside the tree

    //Idx of the first part in this leaf in the array of part
247
    FSize idxOfFirstPartInLeaf =	0;
248 249

    struct LeafToFill{
250 251
        LeafClass * leaf;
        FSize idxOfLeafInPartArray;
252 253 254 255 256 257 258 259 260 261
    };

    leavesPtr.tic();

    LeafToFill * leavesToFill = new LeafToFill[numberOfLeaves];
    memset(leavesToFill,0,sizeof(struct LeafToFill)*numberOfLeaves);
    insertTimer.tic();


    for(FSize idxLeaf = 0; idxLeaf < numberOfLeaves ; ++idxLeaf){
262 263 264
        leavesToFill[idxLeaf].leaf = tree.createLeaf(arrayOfParts[idxOfFirstPartInLeaf].index);
        idxOfFirstPartInLeaf += arrayOfSizeNbLeaves[idxLeaf];
        leavesToFill[idxLeaf].idxOfLeafInPartArray += idxOfFirstPartInLeaf;
265 266 267 268 269 270
    }

    leavesPtr.tac();
    std::cout << "Time needed for creating empty leaves : "<< leavesPtr.elapsed() << " secondes !" << std::endl;


271
    #pragma omp parallel for schedule(auto)
272
    for(FSize idxLeaf=0 ; idxLeaf<numberOfLeaves ; ++idxLeaf ){
273 274 275 276
        //Task consists in copy the parts inside the leaf
        leavesToFill[idxLeaf].leaf->pushArray(&arrayOfParts[leavesToFill[idxLeaf].idxOfLeafInPartArray].position,
                arrayOfSizeNbLeaves[idxLeaf],
                &arrayOfParts[leavesToFill[idxLeaf].idxOfLeafInPartArray].physicalValue);
277 278 279 280 281 282 283 284

    }


    insertTimer.tac();

    std::cout << "Time needed for inserting the parts " << insertTimer.elapsed() << " secondes !" << std::endl;
    double totalTimeForThatStrat = insertTimer.elapsed() + sortTimer.elapsed()
285 286
            + fillTimer.elapsed() + enumTimer.elapsed()
            + leavesPtr.elapsed() + leavesOffset.elapsed();
287 288 289 290 291 292 293 294 295 296
    std::cout << "Total time for that strat : " << totalTimeForThatStrat << " secondes !" << std::endl;
    std::cout << "---------Total------------> Difference : " << regInsert.elapsed() - totalTimeForThatStrat << " secondes !" << std::endl;
    std::cout << "-------On Insert----------> Difference : " << regInsert.elapsed() - insertTimer.elapsed() << " secondes !" << std::endl;


    delete [] arrayOfSizeNbLeaves;
    delete [] arrayOfParts;

    return 0;
}