FDistributedGroupTreeBuilder.hpp 55.1 KB
Newer Older
1
/**
2 3 4 5 6
 * This file contain function to manage the FGroupLinearTree and build a
 * GroupTree with LET
 * The LET is the Local Essential Tree
 * The LET is the symbolic information of leaf for P2P and M2L operation
 *
7 8 9
 * @author benjamin.dufoyer@inria.fr
 */

10

11 12 13 14
#ifndef _FDISTRIBUTED_GROUPTREE_BUILDER_HPP_
#define _FDISTRIBUTED_GROUPTREE_BUILDER_HPP_

#include "FGroupTree.hpp"
15
#include "FOutOfBlockInteraction.hpp"
16
#include <cmath>
17 18 19 20 21
#include <algorithm>

#include <stdint.h>
#include <limits.h>

22
// Define a MPI type for std::size_t
23 24 25 26 27 28 29 30 31 32 33 34 35 36
#if SIZE_MAX == UCHAR_MAX
   #define my_MPI_SIZE_T MPI_UNSIGNED_CHAR
#elif SIZE_MAX == USHRT_MAX
   #define my_MPI_SIZE_T MPI_UNSIGNED_SHORT
#elif SIZE_MAX == UINT_MAX
   #define my_MPI_SIZE_T MPI_UNSIGNED
#elif SIZE_MAX == ULONG_MAX
   #define my_MPI_SIZE_T MPI_UNSIGNED_LONG
#elif SIZE_MAX == ULLONG_MAX
   #define my_MPI_SIZE_T MPI_UNSIGNED_LONG_LONG
#else
   #error "what is happening here?"
#endif

37
#define MAX_SIZE_MPI_MESSAGE 4000000
38

39

40 41
namespace dstr_grp_tree_builder{

42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 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 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210

/**
 * Return the number of MPI Message needed to send this buffer
 * @author benjamin.dufoyer@inria.fr
 * @param  size_buffer  size of the buffer
 * @return              the number of message
 */
unsigned get_nb_mpi_msg(long unsigned size_buffer){
    unsigned nb_message = 1;
    while(size_buffer > MAX_SIZE_MPI_MESSAGE){
        size_buffer -= size_buffer - MAX_SIZE_MPI_MESSAGE;
        nb_message++;
    }
    return nb_message;
}
unsigned get_nb_mpi_msg(int size_obj, int nb_obj_to_send){
    return get_nb_mpi_msg(size_obj*nb_obj_to_send);
}
template<class send_type>
unsigned get_nb_mpi_msg(send_type obj_to_send, int nb_obj_to_send){
    return get_nb_mpi_msg(sizeof(send_type)*nb_obj_to_send);
}

/**
 * This function return the number of element to put into a message
 * @author benjamin.dufoyer@inria.fr
 * @param  obj_to_send      a object to send
 * @param  nb_obj_to_send   number of object to send
 * @return                  the number of element per message
 */
unsigned get_nb_elt_interval(unsigned long int size_buffer){
    return (unsigned)MAX_SIZE_MPI_MESSAGE/(unsigned)size_buffer;
}
template<class send_type>
unsigned get_nb_elt_interval(send_type obj_to_send, int nb_obj_to_send){
    return get_nb_elt_interval(size_of(obj_to_send)*nb_obj_to_send);
}
unsigned get_nb_elt_interval(int size_obj, int nb_obj_to_send){
    return get_nb_elt_interval(size_obj*nb_obj_to_send);
}

/**
 * This function split MPI message if the buffer is too high
 * this funciton run with the irecv_splited
 * She send more 1 message if the buffer have a good size and send more than 1
 * message if the size of the buffer is bigger than MAX_SIZE_MPI_MESSAGE define
 * at the front of this documents
 * @author benjamin.dufoyer@inria.fr
 * @param  conf                 MPI Conf
 * @param  addr_send            Vector address of data to send
 * @param  idx_to_send          index where the data to send start
 * @param  nb_element_to_send   number of element to send
 * @param  destination          number of the destination proc
 * @param  tag                  MPI tag
 */
template<class class_sended>
void isend_splited(const inria::mpi_config& conf,
                   std::vector<class_sended>* addr_send,
                   unsigned* idx_to_send,
                   std::size_t nb_element_to_send,
                   int& destination,
                   int tag = 1 )
{
    // getting usefull variable
    unsigned  size_buffer = (unsigned)sizeof(class_sended)*(unsigned)nb_element_to_send;
    unsigned nb_message = get_nb_mpi_msg(size_buffer);
    // Check the number of message
    if(nb_message > 1){
        unsigned nb_elt_interval = get_nb_elt_interval(size_buffer);
        // Send all messages
        unsigned nb_elt;
        for(unsigned i = 0; i < nb_message ; ++i ){
            if(nb_element_to_send > nb_elt_interval){
                nb_elt = nb_elt_interval;
                nb_element_to_send -= nb_elt_interval;
            } else {
                // last message
                nb_elt = (unsigned)nb_element_to_send;
            }
            conf.comm.isend(
                &addr_send->data()[*idx_to_send],
                (int)sizeof(class_sended)*(int)nb_elt,
                MPI_CHAR,
                destination,tag
            );
            *idx_to_send += nb_elt;
        }
    } else {
        // send 1 message if the buffer is not too big
        conf.comm.isend(
            &addr_send->data()[*idx_to_send],
            (int)sizeof(class_sended)*(int)nb_element_to_send,
            MPI_CHAR,
            destination,tag
        );
        *idx_to_send += (unsigned)nb_element_to_send;
    }
}

/**
 * This function post 1 or more MPI Irecv. She check if the buffer is too big
 * She modify dynamicly the vector of request for the waitAll, she realloc
 * every time when it's needly
 * @author benjamin.dufoyer@inria.fr
 * @param  conf                 MPI conf
 * @param  vector_request       Adress of the vector with MPI status
 * @param  idx_request          Index of the current MPI status
 * @param  addr_recev           Address of the vector where data will be stock
 * @param  idx_reception        Index of the vector where data will be stock
 * @param  nb_element_to_recv   Number of element to recv
 * @param  destination          number of the destination proc
 * @param  tag                  tag of the communication
 */
template<class class_recv>
void irecv_splited(const inria::mpi_config& conf,
                    std::vector<inria::mpi::request>* vector_request,
                    int* idx_request,
                    std::vector<class_recv>* addr_recev,
                    unsigned* idx_reception,
                    std::size_t  nb_element_to_recv,
                    int& destination,
                    int tag = 1 )
{
    // getting usefull variable
    unsigned long int size_buffer = sizeof(class_recv)*nb_element_to_recv;
    unsigned nb_message = get_nb_mpi_msg(size_buffer);
    // check if this function is call at good time
    if( nb_message > 1){
        unsigned nb_elt_interval = get_nb_elt_interval(size_buffer);
        // resize the vector of request
        {
            // we do -1 because, we don't count the message already allocate
            unsigned current_nb_msg = (unsigned)vector_request->size()-1;
            vector_request->resize(current_nb_msg+nb_message);
        }
        // send the good number of message
        unsigned nb_elt = 0;
        for(unsigned i = 0; i < nb_message ; ++i ){
            // compute the number of element recev
            if(nb_element_to_recv > nb_elt_interval){
                nb_elt = nb_elt_interval;
                nb_element_to_recv -= nb_elt_interval;
            } else {
                // last message
                nb_elt = (unsigned)nb_element_to_recv;
            }
            vector_request->data()[*idx_request] =
                conf.comm.irecv(
                    &addr_recev->data()[*idx_reception],
                    (int)sizeof(class_recv)*(int)nb_elt,
                    MPI_CHAR,
                    destination,tag
                );
                *idx_reception += nb_elt;
                *idx_request+=1;
        }
    } else {
        vector_request->data()[*idx_request] =
            conf.comm.irecv(
                &addr_recev->data()[*idx_reception],
                (int)sizeof(class_recv)*(int)nb_element_to_recv,
                MPI_CHAR,
                destination,tag
            );
        *idx_request += 1;
        *idx_reception =+ (unsigned)nb_element_to_recv;
    }
}

211 212 213
/**
 * fill_new_linear_tree this function fill the new linear tree with the value of the current linear tree who are conserved.
 * @author benjamin.dufoyer@inria.fr
214 215 216 217 218 219
 * @param  source               old linear tree
 * @param  destination          new linear tree
 * @param  nb_leaf_recev_left   number of leaf recev from left
 * @param  nb_leaf_recev_right  number of leaf recev from right
 * @param  nb_leaf_send_left    number of leaf send to left
 * @param  nb_leaf_send_right   number of leaf send to right
220 221 222 223 224 225 226 227 228 229 230 231 232 233 234
 */
template<class node_t>
void fill_new_linear_tree(
    std::vector<node_t>* source,
    std::vector<node_t>* destination,
    unsigned nb_leaf_recev_left,
    unsigned nb_leaf_recev_right,
    unsigned nb_leaf_send_left,
    unsigned nb_leaf_send_right
){
    unsigned min_copy = nb_leaf_send_left;
    unsigned max_copy = (unsigned)source->size() - nb_leaf_send_right;
    unsigned min_destination = nb_leaf_recev_left;
    unsigned max_destination = (unsigned)destination->size() - nb_leaf_recev_right;

235 236 237 238 239 240 241
    // find the smallest interval
    unsigned destination_interval = max_destination-min_destination;
    unsigned source_interval = max_copy-min_copy;
    if(source_interval < destination_interval){
        memcpy(&destination->data()[min_destination],&source->data()[min_copy],sizeof(node_t)*source_interval);
    } else {
        memcpy(&destination->data()[min_destination],&source->data()[min_copy],sizeof(node_t)*destination_interval);
242 243 244 245
    }
}

/**
246 247 248 249
 * This function redistribute leaf between proc according to groupSize
 * [RESTRICTION] : the number of leaf is very important, we cannot be in case
 *                  where a proc havn't leaf
 * the leaf on proc according to the groupsize
250 251 252 253 254 255
 * on a distributed linear tree
 * @author benjamin.dufoyer@inria.fr
 * @param  conf         MPI conf
 * @param  linear_tree  current distributed linear tree
 * @param  group_size   size of the group of leaf
 */
256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399
 template<class node_t>
 void parrallel_build_block(const inria::mpi_config& conf, std::vector<node_t>* linear_tree, const int& group_size){
     int  nb_local_leaf  = (int)linear_tree->size();
     const int  nb_proc        = conf.comm.size();
     int* array_global_nb_leaf  = (int *)malloc(sizeof(int) * nb_proc); //nb leaf
     const int  my_rank        = conf.comm.rank();
     // Check if i have leaf on my proc
     FAssert(nb_local_leaf > 0);
     // Distribute the local number of leaf to every process
     conf.comm.allgather(&nb_local_leaf,
                     1,
                     MPI_INT,
                     array_global_nb_leaf,
                     1,
                     MPI_INT);

     int nb_global_leaf = 0;
     for(int i = 0 ; i < nb_proc ; i++)
         nb_global_leaf += array_global_nb_leaf[i];

     int nb_global_group = nb_global_leaf / group_size;
     int nb_local_group  = nb_global_group / nb_proc;
     int nb_leaf_needed  = nb_local_group * group_size;
     // Check if we habe enought leafs for every proc
     if( (nb_leaf_needed*(nb_proc-1)) > nb_global_leaf ){
         std::cout << " nb_leaf_needed : " << nb_leaf_needed << std::endl;
         std::cout << " nb_global_leaf : " << nb_global_leaf << std::endl;
         std::cout << " res :  " << (nb_leaf_needed*(nb_proc-1)) << std::endl;
     }
     FAssert( (nb_leaf_needed*(nb_proc-1)) < nb_global_leaf );

     struct message_info{
         int process_rank;
         int nb_leaf;
     };
     // We stock the future interaction in 2 vector
     std::vector<message_info> interaction_send;
     std::vector<message_info> interaction_recev;

     // The number of leaf send and revev from left
     // it's used to fill the new linear_tree
     int nb_leaf_recev_left = 0;
     int nb_leaf_recev_right = 0;

     int nb_leaf_send_right  = 0;
     int nb_leaf_send_left   = 0;
     // COMPUTE FOR LEFT PROCESS
     // Check to know if the current proc need to send leaf
     // The compute is from left to right because it's the right process
     // who don't have a fix number of particle
     if(!my_rank == 0){ //The first process don't have computation on his left
         for(int i = 1 ; i < my_rank ; i++ ){
             array_global_nb_leaf[i] += array_global_nb_leaf[i-1];
         }
         // Check if on left process need leaf or have too many leaf
         // The idea is, if the total of leaf is a multiple of nb_leaf_needed, no communication is needed
         int nb_leaf_to_send = array_global_nb_leaf[my_rank-1] - (nb_leaf_needed * my_rank);
         if(nb_leaf_to_send != 0 ){
             message_info mess;
             mess.process_rank = my_rank-1;
             //Check if the left process have too much leaf
             if(nb_leaf_to_send > 0){
                 mess.nb_leaf = nb_leaf_to_send;
                 interaction_recev.push_back(mess);
                 //Update the array global with future value for the current proc
                 array_global_nb_leaf[my_rank] += nb_leaf_to_send;
                 nb_leaf_recev_left += nb_leaf_to_send;
             // Else the left process don't have leaf to create block_size
             } else {
                 nb_leaf_to_send = -nb_leaf_to_send;
                 mess.nb_leaf = nb_leaf_to_send;
                 interaction_send.push_back(mess);
                 //Update the array global with future value for the current proc
                 array_global_nb_leaf[my_rank] -= nb_leaf_to_send;
                 nb_leaf_send_left += nb_leaf_to_send;
             }
         }
     }

     // COMPUTE FOR RIGHT PROCESS
     // every proc compute his number of leaf after the first communication
     int nb_leaf_to_send = array_global_nb_leaf[my_rank] - nb_leaf_needed;
     // The last proc don't have proc on his right
     if( (my_rank+1) != nb_proc) {
         if(nb_leaf_to_send != 0){
             message_info mess;
             mess.process_rank = my_rank+1;
             //Check if the current process DON'T have too much leaf
             if(nb_leaf_to_send < 0){
                 nb_leaf_to_send = -nb_leaf_to_send;
                 mess.nb_leaf = nb_leaf_to_send;
                 interaction_recev.push_back(mess);
                 // Else the left process don't have leaf to form block_size
                 //Update the array global with future value for the current proc
                 array_global_nb_leaf[my_rank] += nb_leaf_to_send;
                 nb_leaf_recev_right += nb_leaf_to_send;
             } else {
                 mess.nb_leaf = nb_leaf_to_send;
                 interaction_send.push_back(mess);
                 //Update the array global with future value for the current proc
                 array_global_nb_leaf[my_rank] -= nb_leaf_to_send;
                 nb_leaf_send_right += nb_leaf_to_send;
             }
         }
     }

     // Now we have 2 vector with all interaction with other process
     // in the first we will post every recev message
     // in a second time we post every send message

     // We declare the new linear_tree who have the new number of cell
     std::vector<node_t>* new_linear_tree = new std::vector<node_t>(array_global_nb_leaf[my_rank]);
     //Posting receiv message
     //For every interaction
     // Array of request to know the status
     inria::mpi::request tab_mpi_status[interaction_recev.size()];
     for(unsigned i = 0 ; i < interaction_recev.size(); i++ ){
         // Size of buffer
         int size_recev = (int) (sizeof(node_t)*interaction_recev[i].nb_leaf);
         // Compute the pointer to write result
         unsigned start = 0;
         if(my_rank < interaction_recev[i].process_rank){
             start = (unsigned)new_linear_tree->size() - interaction_recev[i].nb_leaf;
         }
         // Sending request
         tab_mpi_status[i] =
             conf.comm.irecv(&new_linear_tree->data()[start],
                         size_recev,
                         MPI_CHAR,
                         interaction_recev[i].process_rank,1);
     }

     ////Posting sending message
     for(unsigned i = 0 ; i < (unsigned)interaction_send.size(); i++ ){
         int size_send = (int)sizeof(node_t)*interaction_send[i].nb_leaf;
         // Compute the pointer to send cell
         unsigned start = 0;
         if(my_rank < interaction_send[i].process_rank){
             start = (unsigned)linear_tree->size() - interaction_send[i].nb_leaf;
         }

         //sending leaf
         conf.comm.isend(&linear_tree->data()[start],
                        size_send,
400
                        MPI_CHAR,
401 402
                        interaction_send[i].process_rank,1);
     }
403

404 405 406 407 408 409 410 411 412 413
     // Filling vector with the old linear_tree
     // The function need to know how many cell the current process send on
     // the right and on the left because the MPI request write on the same
     //  pointer
     fill_new_linear_tree(linear_tree,
                          new_linear_tree,
                          nb_leaf_recev_left,
                          nb_leaf_recev_right,
                          nb_leaf_send_left,
                          nb_leaf_send_right);
414

415

416 417


418 419
     // waiting for the end of MPI request
     inria::mpi::request::waitall(interaction_recev.size(),tab_mpi_status);
420

421 422 423
     free(array_global_nb_leaf);
     // swaping linear_tree pointer
     std::swap(*linear_tree,*new_linear_tree);
424

425
 }
426

427 428 429 430 431 432 433
/**
 * The idea of this function is to know the repartition of particle
 * morton index on every proc
 * Every proc send he first and his last particle morton index with a MPI
 * allgather
 * @author benjamin.dufoyer@inria.fr
 * @param   conf MPI
434 435
 * @param   particle : vector where particles a stock, they will be sorted by
 *                     morton index BEFORE calling function
436 437 438
 * @param   particle_repartition : a empty vector with the size of the
 *                                 number  of process
 */
DUFOYER Benjamin's avatar
DUFOYER Benjamin committed
439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456

 template<class type1_t,
          class type2_t>
 void share_particle_division(
     const inria::mpi_config& conf,
     std::pair<type1_t,type2_t> my_pair,
     std::vector<std::pair<type1_t,type2_t>>& particle_repartition
 ){
     conf.comm.allgather(
         &my_pair,
         sizeof(my_pair),
         MPI_CHAR,
         particle_repartition.data(),
         sizeof(my_pair),
         MPI_CHAR);
 }


457 458 459
template<class particle_t,
         class type1_t,
         class type2_t>
460
void share_particle_division(
461 462 463 464 465
    const inria::mpi_config& conf,
    std::vector<particle_t>& particle,
    std::vector<std::pair<type1_t,type2_t>>& particle_repartition)
{
    FAssert(particle_repartition.size() == (unsigned)conf.comm.size());
466 467
    FAssert(particle.size() > 0);

468 469 470 471
    std::pair<type1_t,type2_t> my_idx;
    my_idx.first = particle.front().morton_index;
    my_idx.second = particle.back().morton_index;
    // Distribute the local min max morton_index to every process
DUFOYER Benjamin's avatar
DUFOYER Benjamin committed
472
    share_particle_division(conf,my_idx,particle_repartition);
473 474 475 476 477
}



/**
478
 * IDEA Factoriser la fin avec la fonction pour le M2L IDEA
479 480 481 482 483 484 485 486 487 488 489
 * This function compute the morton index of every leaf needed for the P2P
 * First we compute every morton index needed for every leaf
 * We sort the result
 * And to finish we remove multiple iteration
 * @author benjamin.dufoyer@inria.fr
 * @param   GroupOctreeClass Actual group tree
 * @param global_min_m_idx min morton index of the simulation
 * @param global_max_m_idx max morton index of the simulation
 * @return  std::vector with morton_index of cell who is needed to compute P2P
 */
 template<class GroupOctreeClass>
490
std::vector<MortonIndex> get_leaf_P2P_interaction(
491 492 493 494 495 496 497
    GroupOctreeClass& tree,
    const MortonIndex& global_min_m_idx,
    const MortonIndex& global_max_m_idx,
    const MortonIndex& local_min_m_idx,
    const MortonIndex& local_max_m_idx
){
    // 26 is for every interaction
498
    std::vector<std::size_t> externalInteractionsLeafLevel(tree.getTotalNbLeaf()*26,0);
499 500
    // Reset interactions
    // idx to know where we are in the vector
501
    std::size_t  idx_vector= 0;
502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541
    // First leaf level
    {
        // We iterate on every particle group
        for(int idxGroup = 0 ; idxGroup < tree.getNbParticleGroup() ; ++idxGroup){
            // get the particle group
            // it's ugly but, if i use template, it's not convenient
            auto* containers = tree.getParticleGroup(idxGroup);
            {
                // Iterate on every leaf
                for(int leafIdx = 0;
                    leafIdx < containers->getNumberOfLeavesInBlock();
                    ++leafIdx){
                    // Getting the morton index of the leaf
                    const MortonIndex mindex = containers->getLeafMortonIndex(leafIdx);

                    // Define array to receive computed information
                    MortonIndex interactionsIndexes[26];
                    int interactionsPosition[26];
                    //Getting coordinate of the current leaf
                    FTreeCoordinate coord(mindex);
                    // Getting neigbors
                    int counter = coord.getNeighborsIndexes(tree.getHeight(),interactionsIndexes,interactionsPosition);

                    // Iterate on every neighbors
                    for(int idxInter = 0 ; idxInter < counter ; ++idxInter){
                        // Check if the current proc already have the leaf
                        if(interactionsIndexes[idxInter] < local_min_m_idx  || interactionsIndexes[idxInter] > local_max_m_idx ){
                            // Check if the leaf exist
                            if(interactionsIndexes[idxInter] >= global_min_m_idx && interactionsIndexes[idxInter] <= global_max_m_idx ){
                                externalInteractionsLeafLevel[idx_vector] = interactionsIndexes[idxInter];
                                ++idx_vector;
                            }
                        }
                    }
                }
            }
        }
    }
    if(idx_vector != 0) {
        // Sort every morton index
542
        FQuickSort<std::size_t>::QsSequential(externalInteractionsLeafLevel.data(),idx_vector);
543
        // Compute the number of different morton index
544 545
        std::size_t nb_leaf = 1;
        std::size_t last_leaf = externalInteractionsLeafLevel[0];
546 547 548 549 550 551 552
        for(unsigned i = 1 ; i < idx_vector ; ++i){
            if(last_leaf != externalInteractionsLeafLevel[i]){
                last_leaf = externalInteractionsLeafLevel[i];
                nb_leaf++;
            }
        }
        // Alloc the returned vector
553
        std::vector<MortonIndex> leaf_needed(nb_leaf,0);
554
        // Fill the returned vector
555
        MortonIndex current_idx = 1;
556 557 558 559 560 561 562 563
        last_leaf = externalInteractionsLeafLevel[0];
        leaf_needed[0] = externalInteractionsLeafLevel[0];
        for(unsigned i = 1 ; i < idx_vector ; ++i){
            if(last_leaf != externalInteractionsLeafLevel[i]){
                last_leaf = externalInteractionsLeafLevel[i];
                leaf_needed[current_idx] = externalInteractionsLeafLevel[i];
                ++current_idx;
            }
564
            if(current_idx == (int)nb_leaf)
565 566 567 568
                break;
        }
        return leaf_needed;
    } else {
569
        std::vector<MortonIndex> leaf_needed(0,0);
570 571
        return leaf_needed;
    }
572
}
573

574 575

/**
DUFOYER Benjamin's avatar
DUFOYER Benjamin committed
576
 * IDEA on peut factoriser le post traitement de du P2P avec celui qui est fait
577
 *  ici IDEA
DUFOYER Benjamin's avatar
DUFOYER Benjamin committed
578
 *
579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604
 * This function compute the leaf needed for the M2L operation
 * We take every leaf of the tree, get her parent, get the neigbors of
 * the parents and take every child of the parent's neighbors.
 * After we sort the result and remove duplicate data
 * @author benjamin.dufoyer@inria.fr
 * @param global_min_m_idx  global min morton index
 * @param global_max_m_idx  global max morton index
 * @param local_min_m_idx   local min morton index
 * @param local_max_m_idx   local max morton index
 * @param tree              local group tree
 * @param dim               optionnal parameter to compute in other dimension
 * @return vector with leaf needed for the M2M operation
 */
template<class GroupOctreeClass>
std::vector<MortonIndex> get_leaf_M2L_interaction_at_level(
        const MortonIndex& global_min_m_idx,
        const MortonIndex& global_max_m_idx,
        const MortonIndex& local_min_m_idx,
        const MortonIndex& local_max_m_idx,
        int level,
        GroupOctreeClass& tree,
        int dim = 3)
{
    // idx to fill the vector
    std::size_t idx_vector = 0;
    // All External leaf
DUFOYER Benjamin's avatar
DUFOYER Benjamin committed
605
    std::vector<MortonIndex> external_leaf(tree.getNbCellGroupAtLevel(level)*tree.getNbElementsPerBlock()*216,0);
606 607 608 609 610 611 612 613 614 615 616 617
    // iterate on the group
    for(int idxGroup = 0 ; idxGroup < tree.getNbCellGroupAtLevel(level) ; ++idxGroup){
        auto* containers = tree.getCellGroup(level,idxGroup);
        MortonIndex curr_m_idx;
        // +1 is to pass the test at the first try
        for(int leafIdx = 0;
                leafIdx < containers->getNumberOfCellsInBlock();
                ++leafIdx){
            // Getting the current morton index
            curr_m_idx  = containers->getCellMortonIndex(leafIdx);
            // Compute the morton index of the father
            // If it's a new father
618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633
            // Compute coordinate
            MortonIndex interactionsIndexes[216];
            int interactionsPosition[216];
            FTreeCoordinate coord(curr_m_idx);
            // Getting neigbors of the father
            int counter = coord.getInteractionNeighbors(level,interactionsIndexes,interactionsPosition);
            for(int idxNeighbor = 0 ; idxNeighbor < counter ; ++idxNeighbor){
                MortonIndex tmp = interactionsIndexes[idxNeighbor];
                if( tmp    >= global_min_m_idx
                    && tmp <= global_max_m_idx)
                {
                    if(tmp < local_min_m_idx ||
                       tmp > local_max_m_idx){
                        //Stock the leaf
                        if(idx_vector > external_leaf.size()){
                            std::cout << "ERROR " << std::endl;
DUFOYER Benjamin's avatar
DUFOYER Benjamin committed
634
                        }
635 636
                        external_leaf[idx_vector] = tmp;
                        ++idx_vector;
DUFOYER Benjamin's avatar
DUFOYER Benjamin committed
637
                    }
638 639
                }
            } // end for neigbors
640 641 642 643 644 645 646
        } // end for leaf
    } // end for group
    // if we have leaf in the vector
    if(idx_vector != 0 ){
        //sort the result
        FQuickSort<MortonIndex>::QsSequential(external_leaf.data(),idx_vector);
        // compute the number of leaf
647

648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680
        std::size_t nb_leaf = 1;
        MortonIndex last_leaf = external_leaf[0];
        for(unsigned i = 1 ; i < idx_vector ; ++i){
            if(last_leaf != external_leaf[i]){
                last_leaf = external_leaf[i];
                nb_leaf++;
            }
        }
        // Alloc the returned vector
        std::vector<MortonIndex> leaf_needed(nb_leaf,0);
        // Fill the returned vector
        MortonIndex current_idx = 1;
        last_leaf = external_leaf[0];
        leaf_needed[0] = external_leaf[0];
        for(unsigned i = 1 ; i < idx_vector ; ++i){
            if(last_leaf != external_leaf[i]){
                last_leaf = external_leaf[i];
                leaf_needed[current_idx] = external_leaf[i];
                ++current_idx;
            }
            if(current_idx == (long)nb_leaf)
                break;
        }
        return leaf_needed;
    } else {
        std::vector<MortonIndex> leaf_needed(0,0);
        return leaf_needed;
    }


}


681 682 683 684 685 686 687 688 689 690
/**
 * The goal of this function is to concat the two vector
 * and erase duplicate data
 * first we compute the number of leaf
 * second we fill the vector with erase duplicate data
 * @author benjamin.dufoyer@inria.fr
 * @param  leaf_P2P vector of P2P leaf
 * @param  leaf_M2L vector of M2L leaf
 * @return vector with morton index of leaf needed
 */
691 692 693
std::vector<MortonIndex> concat_M2L_P2P(
    std::vector<MortonIndex>& leaf_P2P,
    std::vector<MortonIndex>& leaf_M2L)
694 695 696 697 698 699 700 701 702
{
    // Checking if one of vector is empty
    if(leaf_P2P.size() == 0)
        return leaf_M2L;
    if(leaf_M2L.size() == 0)
        return leaf_P2P;
    // Compute the total number of leaf needed
    unsigned idx_P2P = 0;
    unsigned idx_M2L = 0;
703
    std::size_t nb_leaf = 0;
704
    while(idx_P2P != leaf_P2P.size() && idx_M2L != leaf_M2L.size()){
705 706
        MortonIndex curr_P2P = leaf_P2P[idx_P2P];
        MortonIndex curr_M2M = leaf_M2L[idx_M2L];
707 708 709 710 711 712 713 714 715 716 717 718 719
        if(curr_P2P != curr_M2M){
            if(curr_P2P < curr_M2M){
                ++idx_P2P;
            } else {
                ++idx_M2L;
            }
        } else {
            ++idx_P2P;
            ++idx_M2L;
        }
        ++nb_leaf;
    }
    if(idx_P2P == leaf_P2P.size()){
720
        nb_leaf += (leaf_M2L.size()) - idx_M2L;
721
    } else {
722
        nb_leaf += (leaf_P2P.size()) - idx_P2P;
723 724
    }
    // Allocate the vector
725
    std::vector<MortonIndex> leaf_needed(nb_leaf,0);
726 727
    idx_P2P = 0;
    idx_M2L = 0;
728
    std::size_t idx_leaf = 0;
729 730
    // fill the vector
    while(idx_P2P != leaf_P2P.size() && idx_M2L != leaf_M2L.size()){
731 732
        MortonIndex curr_P2P = leaf_P2P[idx_P2P];
        MortonIndex curr_M2M = leaf_M2L[idx_M2L];
733 734 735 736 737 738 739 740
        if(curr_P2P != curr_M2M){
            if(curr_P2P < curr_M2M){
                leaf_needed[idx_leaf] = leaf_P2P[idx_P2P];
                ++idx_P2P;
            } else {
                leaf_needed[idx_leaf] = leaf_M2L[idx_M2L];
                ++idx_M2L;
            }
741
            //++nb_leaf;
742
        } else {
743
            leaf_needed.at(idx_leaf) = leaf_M2L[idx_M2L];
744 745 746 747 748 749
            ++idx_P2P;
            ++idx_M2L;
        }
        idx_leaf++;
    }
    // Copy the rest of leaf with a memcpy
750 751 752 753 754 755 756 757 758 759 760 761
    if(idx_leaf < nb_leaf){
        void* destination = &leaf_needed.data()[idx_leaf];
        void* source;
        std::size_t num = 0;
        if(idx_P2P == leaf_P2P.size()){
            source = &leaf_M2L[idx_M2L];
            num = sizeof(MortonIndex)* ((leaf_M2L.size()-1) - idx_M2L);
        } else {
            source = &leaf_P2P[idx_P2P];
            num = sizeof(MortonIndex)* ((leaf_P2P.size()-1) - idx_P2P);
        }
        memcpy(destination,source,num);
762 763 764 765 766
    }
    return leaf_needed;

}

767 768


769 770 771 772 773 774 775 776 777 778
/**
 * This function compute the number of leaf needed on every proc,
 * she check if leaf exit, if a leaf doesn't exist this function
 * remove here at the end of the function
 * @author benjamin.dufoyer@inria.fr
 * @param  needed_leaf vector with the morton index of every leaf
 * needed for my proc
 * @param  particle_repartition vector with the repartation of particle on
 * every proc
 * @param  conf conf MPI
779
 * @return vector<std::size_t> with a size nb_proc*nb_proc
780 781
 * it's a interaction matrix
 */
782
std::vector<std::vector<std::size_t>> get_matrix_interaction(
783 784
    std::vector<MortonIndex>& needed_leaf,
    std::vector<std::pair<MortonIndex,MortonIndex>>& particle_distribution,
785 786 787 788 789
    const inria::mpi_config& conf)
{
    // Getting MPI Info
    const int  nb_proc        = conf.comm.size();
    // Alloc interaction matrix
790 791
    std::vector<std::vector<std::size_t>> matrix_interaction(2,std::vector<std::size_t>(nb_proc,0));
    std::vector<std::size_t> global_matrix_interaction(nb_proc,0);
792 793
    // Initialise idx on particle_distribution
    std::size_t idx_part = 0;
794 795
    // Interate on every leaf to know where she is
    for(unsigned idx_leaf = 0; idx_leaf < needed_leaf.size(); ++idx_leaf){
796
        MortonIndex current_leaf = needed_leaf[idx_leaf];
797
        // if she is on the current proc
798 799
        if(current_leaf >= particle_distribution[idx_part].first
        && current_leaf <= particle_distribution[idx_part].second){
800
            matrix_interaction[0][idx_part] += 1;
801 802
        } else {
            // While the current leaf is not on the good interval
803
            while(particle_distribution[idx_part].second < current_leaf){
804 805
                idx_part += 1;
            }
806
            if(particle_distribution[idx_part].first > current_leaf){
807 808 809 810 811
                // in this case the leaf is not in interval, so she doesn't exist
                needed_leaf[idx_leaf] = 0;
            } else {
                // In the case it's a normal case, we juste increment the
                // number of leaf send at the proc idx_part
812
                matrix_interaction[0][idx_part] += 1;
813 814 815 816
            }
        }
    }
    // now we have the number of leaf to send at every proc
817 818 819 820 821 822 823
    // we proceed a AlltoAll to share this information at every proc
    conf.comm.alltoall(matrix_interaction[0].data(),
                       1,
                       my_MPI_SIZE_T,
                       matrix_interaction[1].data(),
                       1,
                       my_MPI_SIZE_T);
824 825 826
    // removing bad leaf
    needed_leaf.erase(std::remove(needed_leaf.begin(),needed_leaf.end(),0),needed_leaf.end());

827
    return {begin(matrix_interaction),end(matrix_interaction)};
828 829
}

830

831 832 833 834 835 836 837 838 839 840 841
/**
* This function compute the number of block needed to send all leaf
* stock in leaf_needed.
* This function return a vector with all idx of block needed by the proc
* @author benjamin.dufoyer@inria.fr
* @param  tree         GroupTree
* @param  leaf_needed  Vector where leaf are stock
* @return Vector with all block idx
*/
template<class GroupOctreeClass>
std::vector<MortonIndex> get_nb_block_from_leaf(GroupOctreeClass& tree,
DUFOYER Benjamin's avatar
DUFOYER Benjamin committed
842 843
                  MortonIndex* leaf_needed,
                  std::size_t nb_leaf)
844 845
{
    std::vector<MortonIndex> block_to_send(tree.getNbParticleGroup(),0);
846 847
    if(nb_leaf == 0)
        return {block_to_send.begin(),block_to_send.begin()};
848 849 850 851 852
    // declaration of idx varaibles
    unsigned idx_vector = 0;
    unsigned idx_leaf = 0;
    // iterate on every group
    for(int idx_group = 0 ; idx_group < tree.getNbParticleGroup() ; ++idx_group){
853
        if(idx_leaf >= nb_leaf)
854 855 856 857
            break;
        // get the current block
        auto* container = tree.getParticleGroup(idx_group);
        // get first leaf in this interval
858
        while( idx_leaf < nb_leaf && container->getStartingIndex() > leaf_needed[idx_leaf]){
859 860
            ++idx_leaf;
        }
861
        if(idx_leaf >= nb_leaf)
862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879
            break;
        while( container->getEndingIndex() < leaf_needed[idx_leaf] &&
            idx_leaf < nb_leaf){
            // if the leaf exist, keep the leaf
            if(container->exists(leaf_needed[idx_leaf])){
                block_to_send[idx_vector] = idx_group;
                ++idx_vector;
                ++idx_leaf;
                break;
            }
            ++idx_leaf;
        }
        if(idx_leaf == nb_leaf)
            break;
    }
    return {block_to_send.begin(),block_to_send.begin()+idx_vector};
}
/*
DUFOYER Benjamin's avatar
DUFOYER Benjamin committed
880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904
 template<class GroupOctreeClass>
 std::vector<MortonIndex> get_nb_block_from_node(GroupOctreeClass& tree,
              MortonIndex* node_needed,
              std::size_t nb_node,
              int level,
              std::vector<bool>* block_already_send)
 {
     std::vector<int> block_to_send(tree.getNbCellGroupAtLevel(level),0);
     int idx_vect = 0 ;
     // iterate of every node
     for(unsigned i = 0 ; i < nb_node; ++i){
         // iteracte of every block
         for(unsigned idxGroup = 0 ; idxGroup < (unsigned)tree.getNbCellGroupAtLevel(level) ; ++idxGroup){
             // If the block is not already send
             if(block_already_send->at(idxGroup) == false){
                 auto* containers = tree.getCellGroup(level,idxGroup);
                 if(containers->isInside(node_needed[i])){
                     block_to_send[idx_vect] = idxGroup;
                     ++idx_vect;
                     block_already_send->at(idxGroup) = true;
                 }
             }
         }
     }
     return {block_to_send.begin(),block_to_send.begin()+idx_vect};
905
 }*/
906

DUFOYER Benjamin's avatar
DUFOYER Benjamin committed
907

908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928
 template<class GroupOctreeClass>
 std::vector<MortonIndex> get_nb_block_from_node(GroupOctreeClass& tree,
              MortonIndex* node_needed,
              std::size_t nb_node,
              int level,
              std::vector<bool>* block_already_send)
 {
     int idx_vect = 0 ;
     std::vector<int> block_to_send(tree.getNbCellGroupAtLevel(level),0);

     unsigned idx_node = 0;
     // iterate on every group
     for(unsigned idx_group = 0; idx_group < (unsigned)tree.getNbCellGroupAtLevel(level) ;++idx_group){
         // if the current block hasnt been already send
        if(!block_already_send->at(idx_group)){
            auto* containers = tree.getCellGroup(level,idx_group);
            // check if we have check every node
            if(idx_node == nb_node){
                break;
            }
            // while the morton index of the current node is not high
929
            while(idx_node < nb_node && node_needed[idx_node] < containers->getStartingIndex()){
930 931
                ++idx_node;
            }
932
            while(idx_node < nb_node && node_needed[idx_node] < containers->getEndingIndex()){
933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961
                 if(containers->isInside(node_needed[idx_node])){
                     block_to_send[idx_vect] = idx_group;
                     ++idx_vect;
                     ++idx_node;
                     break;
                 }
                 ++idx_node;
            }
            if(idx_node == nb_node){
                break;
            }
        }
     }
     return {block_to_send.begin(),block_to_send.begin()+idx_vect};
 }

 /*
  * In this function know the block needed by every proc
  * So we need need to compute the parents of every block and send the
  * number of parent's block to every proc so we have
  *
  * 1st Step : Compute the parents block
  * 2nd Step : Send and recev the number of parents block per level
  * @author benjamin.dufoyer@inria.fr
  * @param  nb_block_to_receiv number of block to receiv by every proc per level
  * @param  list_of_block_to_send list of block to send for every proc per level
  * @param  tree         local group octree
 * @param  conf          MPI CONF
  */
DUFOYER Benjamin's avatar
DUFOYER Benjamin committed
962 963 964
 template<class GroupOctreeClass>
 void send_get_number_of_block_node_level(
     std::vector<MortonIndex>& vect_recv,
965 966
     std::vector<std::vector<std::size_t>> global_matrix_interaction,
     //std::vector<std::size_t>& global_matrix_interaction,
DUFOYER Benjamin's avatar
DUFOYER Benjamin committed
967 968 969 970 971 972 973 974 975 976 977 978 979 980 981
     std::size_t& nb_msg_recv,
     GroupOctreeClass& tree,
     std::vector<std::pair<int,int>>& nb_block_to_receiv,
     std::vector<std::pair<int,std::vector<MortonIndex>>>& leaf_to_send,
     int level,
     const inria::mpi_config& conf
 )
 {
     int idx_status = 0;
     int idx_proc = 0;
     inria::mpi::request tab_mpi_status[nb_msg_recv];
     bool leaf_level = (tree.getHeight()-1 == level);
     std::vector<bool> block_already_send(tree.getNbCellGroupAtLevel(level),false);

     // Post the number reception of the number of block
982
     for(unsigned i = 0; i < global_matrix_interaction[0].size() ; ++i)
DUFOYER Benjamin's avatar
DUFOYER Benjamin committed
983 984
     {
         // If we have interaction with this proc
985
         if(global_matrix_interaction[0][i] != 0){
DUFOYER Benjamin's avatar
DUFOYER Benjamin committed
986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000
             nb_block_to_receiv[idx_status].first = idx_proc;
             tab_mpi_status[idx_status] = conf.comm.irecv(
                 &nb_block_to_receiv[idx_status].second,
                     1,
                     MPI_INT,
                     idx_proc,1
                 );
                 idx_status += 1;
         }
         idx_proc += 1;
     }
     idx_proc = 0;
     std::size_t idx_vect = 0;
     idx_status = 0;
     // Posting every send message
1001
     for(unsigned i = 0; i < global_matrix_interaction[1].size() ; ++i){
DUFOYER Benjamin's avatar
DUFOYER Benjamin committed
1002
         // If we have interaction with this proc
1003
         if(global_matrix_interaction[1][i] != 0){
DUFOYER Benjamin's avatar
DUFOYER Benjamin committed
1004 1005 1006 1007 1008
             // Compute the number of leaf
             if(leaf_level){
                 leaf_to_send[idx_status].second = get_nb_block_from_leaf(
                                     tree,
                                     &vect_recv.data()[idx_vect],
1009
                                     global_matrix_interaction[1][i]);
DUFOYER Benjamin's avatar
DUFOYER Benjamin committed
1010 1011 1012 1013
            } else {
                 leaf_to_send[idx_status].second = get_nb_block_from_node(
                                     tree,
                                 &vect_recv.data()[idx_vect],
1014
                                     global_matrix_interaction[1][i],
DUFOYER Benjamin's avatar
DUFOYER Benjamin committed
1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026
                                     level,
                                     &block_already_send);
             }
             int nb_block = (int)leaf_to_send[idx_status].second.size();
             leaf_to_send[idx_status].first = idx_proc;
             // send the number of leaf
             conf.comm.isend(
                 &nb_block,
                 1,
                 MPI_INT,
                 idx_proc,1
             );
1027
             idx_vect += global_matrix_interaction[1][i];
DUFOYER Benjamin's avatar
DUFOYER Benjamin committed
1028 1029 1030 1031 1032 1033 1034 1035 1036
             idx_status += 1;
         }
         idx_proc += 1;
     }
     // Waiting for all request
     if(nb_msg_recv != 0 ){
         inria::mpi::request::waitall(nb_msg_recv,tab_mpi_status);
     }
 }
1037

1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050

/**
 * This function compute the leaf needed by every proc
 * @author benjamin.dufoyer@inria.fr
 * @param  needed_leaf      List of leaf needed by me
 * @param  global_matrix_interaction
 * @param  nb_msg_recv      number of recev message
 * @param  nb_leaf_recv     number of recev leaf
 * @param  conf             MPI conf
 * @return                  Vector with all leaf needed by other proc
 */
std::vector<MortonIndex> send_get_leaf_morton(
    std::vector<MortonIndex>&     needed_leaf,
1051 1052
    std::vector<std::vector<std::size_t>>&     global_matrix_interaction,
    //std::vector<std::size_t>&     global_matrix_interaction,
1053 1054
    std::size_t&                  nb_msg_recv,
    std::size_t&                  nb_leaf_recv,
1055 1056 1057
    const inria::mpi_config& conf)
{
    // allocate tab of mpi status for synchronisazion
1058 1059
    std::vector<inria::mpi::request> tab_mpi_status(nb_msg_recv);

1060
    // Initialiser the reception vector
1061 1062
    std::vector<std::size_t> vect_recv(nb_leaf_recv,0);

1063 1064
    // Initialize every variable
    int    idx_status = 0;
1065
    unsigned idx_vect   = 0;
1066 1067
    int    idx_proc   = 0;
    // Posting every recv message
1068 1069
    for(unsigned i = 0; i < global_matrix_interaction[1].size() ; ++i ){
        if(global_matrix_interaction[1][i] != 0){
1070 1071 1072 1073 1074 1075
            irecv_splited(
                conf,
                &tab_mpi_status,
                &idx_status,
                &vect_recv,
                &idx_vect,
1076
                global_matrix_interaction[1][i],
1077 1078 1079 1080 1081 1082 1083 1084 1085
                idx_proc,1
            );
        }
        idx_proc += 1;
    }

    // Posting every send message
    idx_proc = 0;
    idx_vect = 0;
1086 1087
    for(unsigned i = 0; i < global_matrix_interaction[0].size() ; ++i){
        if(global_matrix_interaction[0][i] != 0){
1088 1089 1090 1091
            isend_splited(
                conf,
                &needed_leaf,
                &idx_vect,
1092
                global_matrix_interaction[0][i],
1093 1094 1095 1096 1097 1098
                idx_proc,1
            );
        }
        idx_proc += 1;
    }
    if(nb_msg_recv != 0 ){
1099
        inria::mpi::request::waitall(tab_mpi_status.size(),tab_mpi_status.data());
1100 1101 1102 1103 1104 1105 1106
    }
    return{begin(vect_recv),end(vect_recv)};

}



1107 1108 1109 1110 1111 1112 1113 1114
/**
 * This struct is used to stock information who wille be send to other proc
 */
struct block_t{
    std::size_t n_block;
    MortonIndex start_index;
    MortonIndex end_index;
    int nb_leaf_in_block;
1115 1116 1117
    // used to show the block
    friend
    std::ostream& operator<<(std::ostream& os, const block_t& n) {
1118
        return os << "--> n_block : " << n.n_block << " start : " << n.start_index << " end : " << n.end_index << " nb_leaf " << n.nb_leaf_in_block  <<  "<--";
1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133 1134 1135 1136 1137 1138 1139
    }

};


/**
 * This function echange symbolic information of block of a GroupTree
 *
 * @author benjamin.dufoyer@inria.fr
 * @param  nb_block_to_receiv It's a vector of pair of int,int the first int
 *                            is for the number of the proc who send me block
 *                            the second int is the number of block he send to
 *                            me
 * @param  block_to_send      it's a vector of pair, the first is the number of
 *                            the proc who i send the block, the second is the
 *                            list of the block needed by the proc
 * @param  nb_msg_recv        it's the number of message i will receiv
 * @param  tree              it's the GroupTree where block are stock
 * @param  conf                it's the MPI conf
 */
template<class GroupOctreeClass>
DUFOYER Benjamin's avatar
DUFOYER Benjamin committed
1140
std::vector<std::vector<block_t>> exchange_block(
1141
    std::vector<std::pair<int,int>> nb_block_to_receiv,
1142
    std::vector<std::pair<int,std::vector<MortonIndex>>> block_to_send,
1143
    GroupOctreeClass& tree,
DUFOYER Benjamin's avatar
DUFOYER Benjamin committed
1144
    int level,
1145 1146 1147 1148
    const inria::mpi_config& conf
)
{
    // declaration of the array of MPI status for synchro
1149 1150 1151 1152 1153 1154 1155 1156
    unsigned nb_message_recv = 0;
    for(unsigned i = 0 ; i <  nb_block_to_receiv.size() ;++i ){
        if(nb_block_to_receiv[i].second != 0 ){
            ++nb_message_recv;
        }
    }
    std::vector<inria::mpi::request> tab_mpi_status(nb_message_recv);

1157
    // Compute the total number of block i will send
1158
    std::size_t total_size = 0;
1159 1160 1161 1162 1163 1164
    for(unsigned i = 0 ; i < block_to_send.size(); ++i){
        total_size += block_to_send[i].second.size();
    }

    // Declaration of the buffer of block
    std::vector<block_t> data_to_send(total_size);
1165
    std::size_t idx_vect_to_send = 0;
1166 1167 1168
    // Filling the buffer of block
    for(unsigned i = 0 ; i < block_to_send.size(); ++i){
        for(unsigned j = 0 ; j < block_to_send[i].second.size() ; j++){
DUFOYER Benjamin's avatar
DUFOYER Benjamin committed
1169
            auto* container = tree.getCellGroup(level ,((int)block_to_send[i].second[j]));
1170
            block_t block_to_add{
1171 1172 1173 1174
                (size_t)block_to_send[i].second[j],
                container->getStartingIndex(),
                container->getEndingIndex(),
                container->getNumberOfCellsInBlock()
1175 1176 1177 1178 1179 1180 1181 1182 1183 1184 1185 1186 1187 1188 1189 1190
            };
            data_to_send[idx_vect_to_send] = block_to_add;
            ++idx_vect_to_send;
        }
    }
    // Posting recv
    std::vector<std::vector<block_t>> block_t_recv(nb_block_to_receiv.size());
    int idx_status = 0;
    for(unsigned i = 0; i < nb_block_to_receiv.size(); ++i)
    {
        // Setting parameter
        int source   = nb_block_to_receiv[i].first;
        int nb_block = nb_block_to_receiv[i].second;
        if(nb_block != 0){
            block_t_recv[i].resize(nb_block);
            // Posting reveiv message
1191 1192 1193 1194 1195 1196 1197 1198
            unsigned idx = 0;
            irecv_splited(
                conf,
                &tab_mpi_status,
                &idx_status,
                &block_t_recv.data()[i],
                &idx,
                nb_block,
1199 1200 1201 1202 1203 1204
                source,1
            );
        }
    }

    // post sending message
1205
    unsigned offset_block = 0;
1206 1207 1208
    for(unsigned i = 0 ; i < block_to_send.size(); ++i){
        // Setting parameters
        int destination = block_to_send[i].first;
1209
        size_t nb_block    = (int)block_to_send[i].second.size();
1210 1211
        // Posting send message
        if(nb_block != 0){
1212 1213 1214 1215 1216
            isend_splited(
                conf,
                &data_to_send,
                &offset_block,
                nb_block,
1217 1218 1219 1220 1221
                destination,1
            );
        }
    }
    // Waiting for all request
1222 1223
    if(nb_message_recv != 0){
        inria::mpi::request::waitall(tab_mpi_status.size(),tab_mpi_status.data());
1224 1225 1226 1227 1228
    }
    return{begin(block_t_recv),end(block_t_recv)};
}


1229 1230 1231 1232 1233 1234 1235 1236 1237 1238 1239 1240 1241 1242 1243 1244 1245 1246 1247 1248 1249 1250 1251 1252 1253 1254 1255 1256 1257 1258 1259 1260 1261 1262 1263 1264 1265 1266 1267 1268 1269 1270 1271 1272 1273 1274 1275 1276 1277 1278 1279 1280 1281 1282 1283 1284 1285 1286 1287 1288 1289 1290 1291 1292 1293 1294 1295 1296 1297 1298 1299 1300 1301 1302 1303 1304 1305 1306 1307 1308 1309 1310 1311 1312 1313 1314 1315 1316 1317 1318 1319 1320 1321 1322 1323 1324 1325 1326 1327 1328 1329 1330 1331 1332 1333 1334 1335 1336 1337 1338 1339 1340

/**
 * Compute the block needed by other proc, this function stock the idx of
 * block into the current_level variable
 * This function have 2 step, first we compute the number of parents of
 * the block needed a the level -1
 * after we allocate the vector and fill me
 * @author benjamin.dufoyer@inria.fr
 * @param  under_level      block for the under level
 * @param  current_level    block for the current level
 * @param  level            level of the tree
 * @param  tree             Local GroupTree
 */
template<class GroupOctreeClass>
void compute_block_node_level(
    std::vector<std::pair<int,std::vector<std::size_t>>>& under_level,
    std::vector<std::pair<int,std::vector<std::size_t>>>& current_level,
    int level,
    GroupOctreeClass& tree
){
    FAssert(under_level.size() == current_level.size() );
    // Iterate on every interaction of the under level
    for(unsigned i = 0 ; i < under_level.size() ; ++i){
        // Init variables for the search
        std::size_t start_m_idx = tree.getCellGroup(level+1,int(tree.getNbCellGroupAtLevel(level+1)-1))->getEndingIndex() ;
        std::size_t last_m_idx  = 0;
        std::size_t current_min_m_idx;
        std::size_t current_max_m_idx;
        int nb_block = 0;
        int last_block = -1;
        // Iterate on every block to find his parents
        for(unsigned j = 0 ; j < under_level[i].second.size() ; j++){
            // get the morton index at level
            current_min_m_idx = (tree.getCellGroup(level+1,
                                    (int)under_level[i].second[j])->getStartingIndex() >> 3);
            current_max_m_idx = (tree.getCellGroup(level+1,
                                    (int)under_level[i].second[j])->getEndingIndex() >> 3);
            // If The parents of this block is unknow, it's a new block to add
            if(current_min_m_idx < start_m_idx || current_max_m_idx > last_m_idx){
                // iterate on all group of the current level to find the
                // good block
                for(int idx_block = 0 ;  idx_block < tree.getNbCellGroupAtLevel(level) ; ++idx_block){
                    // get the current block
                    auto* container = tree.getCellGroup(level,idx_block);
                    // if the morton index is inside
                    if(container->isInside(current_min_m_idx) ||
                       container->isInside(current_max_m_idx) ){
                        // update the number of block
                        // if it's not the last block
                        if(idx_block != last_block){
                            nb_block += 1;
                            last_block = idx_block;
                            // update idx
                            if((std::size_t)container->getStartingIndex() < start_m_idx){
                                start_m_idx = container->getStartingIndex();
                            }
                            if((std::size_t)container->getEndingIndex() > last_m_idx){
                                last_m_idx = container->getEndingIndex();
                            }
                        }
                        // stop seeking for this idx
                        break;
                    }
                }
            }
        }
        // Modify current_level
        current_level[i].first = under_level[i].first;
        current_level[i].second.resize(nb_block);
        // reinit variables
        unsigned idx_vect = 0;
        start_m_idx = tree.getCellGroup(level+1,int(tree.getNbCellGroupAtLevel(level+1)-1))->getEndingIndex();
        last_m_idx  = 0;
        last_block = -1;
        // fill the current level with the same technic
        for(unsigned j = 0 ; j < under_level[i].second.size() ; j++){
            // get the morton index at level
            current_min_m_idx = (tree.getCellGroup(level+1,
                                (int)under_level[i].second[j])->getStartingIndex() >> 3);
            current_max_m_idx = (tree.getCellGroup(level+1,
                                (int)under_level[i].second[j])->getEndingIndex() >> 3);
            // The parents of this block is unknow
            if(current_min_m_idx < start_m_idx || current_max_m_idx > last_m_idx){
                // iterate on all group of the current level
                for(int idx_block = 0 ;  idx_block < tree.getNbCellGroupAtLevel(level) ; ++idx_block){
                    // get the current block
                    auto* container = tree.getCellGroup(level,idx_block);
                    if(container->isInside(current_min_m_idx)  ||
                       container->isInside(current_max_m_idx)){
                        if(idx_block != last_block){
                            // update the number of block
                            current_level[i].second[idx_vect] = idx_block;
                            idx_vect += 1;
                            last_block = idx_block;
                            // update idx
                            if((std::size_t)container->getStartingIndex() < start_m_idx){
                                start_m_idx = container->getStartingIndex();
                            }
                            if((std::size_t)container->getEndingIndex() > last_m_idx){
                                last_m_idx = container->getEndingIndex();
                            }
                        }
                        break;
                    }
                }
            }
        }
    }
}



1341 1342 1343 1344 1345 1346 1347 1348 1349 1350
/**
 * this function return a vector with the symbolic information of the needed
 * block to build a LET Group tree
 *
 * @author benjamin.dufoyer@inria.fr
 * @param  needed_leaf list of the leaf needed
 * @param  global_matrix_interaction matrix of interaction
 * @param  tree local group tree
 * @param  conf MPI conf
 */
1351
template<class GroupOctreeClass>
DUFOYER Benjamin's avatar
DUFOYER Benjamin committed
1352
std::vector<std::vector<block_t>> send_get_symbolic_block_at_level(