src_stdpar/Communication.cpp
// Apply the Unbalanced Recursive Bisection (URB) algorithm to establish the boxes
// (parts of the computational domain) that are assigned to MPI processes. The result
// is stored within dis.ulist data structure
Apply the Unbalanced Recursive Bisection (URB) algorithm to establish the boxes
(parts of the computational domain) that are assigned to MPI processes. The result
is stored within dis.ulist data structure
Apply the Unbalanced Recursive Bisection ( URB ) algorithm to establish the boxes ( parts of the computational domain ) that are assigned to MPI processes . The result is stored within dis . ulist data structure
void Distrib::createDecomposition(Config *conf) {
  Domain &dom = conf->dom_;
  ulist_.resize(nbp_);
  for(int id=0; id<nbp_; id++)
    ulist_[id].pid_ = -1;
  ulist_[0].pid_ = 0;
  ulist_[0].nbp_ = nbp_;
  for(int j=0; j<DIMENSION; j++)
    ulist_[0].xmin_[j] = 0;
  for(int j=0; j<DIMENSION; j++) {
    ulist_[0].xmax_[j] = dom.nxmax_[j] - 1;
    nxmax_[j] = dom.nxmax_[j];
  }
  for(int count = 1; count < nbp_; count++) {
    const float64 cutbound = 10000000000.;
    float64 cutbest = cutbound;
    int cutpart = -1;
    int cutid = -1;
    int cutwidth = -1;
    int cutdir = -1;
    for(int direction = 0; direction < NB_DIMS_DECOMPOSITION; direction++) {
      for(int id = 0; id < count; id++) {
        int nbp = ulist_[id].nbp_;
        if(nbp > 1) {
          int center = floor(nbp * 0.5);
          for(int part = center; part < nbp; part++) {
            const int width = ulist_[id].xmax_[direction] - ulist_[id].xmin_[direction] + 1;
            const float64 fract = static_cast<float64>(part) / static_cast<float64>(nbp);
            const float64 splitscore = 1. - width / static_cast<float64>(dom.nxmax_[direction]);
            const int tmpwidth = round(fract * width);
            const float64 ratioscore = fabs(fract * width - tmpwidth);
            const float64 cutscore = ratioscore + splitscore + cutbound * (tmpwidth < MMAX + 2) + cutbound * ((width - tmpwidth) < (MMAX + 2));
            if(cutscore < cutbest) {
              cutid = id;
              cutbest = cutscore;
              cutpart = part;
              cutwidth = tmpwidth;
              cutdir = direction;
            }
          }
        }
      }
    }
    if(cutbest == cutbound) {
      printf("No cut found to create domain decomposition, reduce nb of processes\n");
    }
    for(int j=0; j<DIMENSION; j++) {
      ulist_[count].xmax_[j] = ulist_[cutid].xmax_[j];
      ulist_[count].xmin_[j] = ulist_[cutid].xmin_[j];
    }
    ulist_[cutid].nbp_          -= cutpart;
    ulist_[cutid].xmin_[cutdir] += cutwidth;
    ulist_[count].nbp_          = cutpart;
    ulist_[count].xmax_[cutdir] = ulist_[cutid].xmin_[cutdir] - 1;
  }
  uint64 msum = 0;
  for(int id = 0; id < nbp_; id++) {
    const uint64 mcontrib = (ulist_[id].xmax_[3] - ulist_[id].xmin_[3] + 1) * (ulist_[id].xmax_[2] - ulist_[id].xmin_[2] + 1)
                          * (ulist_[id].xmax_[1] - ulist_[id].xmin_[1] + 1) * (ulist_[id].xmax_[0] - ulist_[id].xmin_[0] + 1);
    if(master())
      printf("[%d] local domain [%4u:%4u,%4u:%4u,%4u:%4u,%4u:%4u] cost %lu\n", id, ulist_[id].xmin_[0], ulist_[id].xmax_[0], ulist_[id].xmin_[1],
             ulist_[id].xmax_[1], ulist_[id].xmin_[2], ulist_[id].xmax_[2], ulist_[id].xmin_[3], ulist_[id].xmax_[3], mcontrib);
    msum += mcontrib;
  }
  uint64 mref = static_cast<uint64>(dom.nxmax_[3]) * static_cast<uint64>(dom.nxmax_[2]) * static_cast<uint64>(dom.nxmax_[1]) * static_cast<uint64>(dom.nxmax_[0]);
  if(mref != msum) {
    printf("Problem: sum check mref %lf mcalc %lu \n", mref, msum);
    MPI_Barrier(MPI_COMM_WORLD);
    MPI_Finalize();
    exit(1);
  }
  node_ = &(ulist_[pid_]);
  for(int j=0; j<DIMENSION; j++) {
    dom.local_nx_[j]    = node_->xmax_[j] - node_->xmin_[j] + 1;
    dom.local_nxmin_[j] = node_->xmin_[j];
    dom.local_nxmax_[j] = node_->xmax_[j];
  }
}
void Distrib :: createDecomposition ( Config * conf ) { Domain & dom = conf -> dom_ ; ulist_ . resize ( nbp_ ) ; for ( int id = 0 ; id < nbp_ ; id ++ ) ulist_ [ id ] . pid_ = -1 ; ulist_ [ 0 ] . pid_ = 0 ; ulist_ [ 0 ] . nbp_ = nbp_ ; for ( int j = 0 ; j < DIMENSION ; j ++ ) ulist_ [ 0 ] . xmin_ [ j ] = 0 ; for ( int j = 0 ; j < DIMENSION ; j ++ ) { ulist_ [ 0 ] . xmax_ [ j ] = dom . nxmax_ [ j ] - 1 ; nxmax_ [ j ] = dom . nxmax_ [ j ] ; } for ( int count = 1 ; count < nbp_ ; count ++ ) { const float64 cutbound = 10000000000. ; float64 cutbest = cutbound ; int cutpart = -1 ; int cutid = -1 ; int cutwidth = -1 ; int cutdir = -1 ; for ( int direction = 0 ; direction < NB_DIMS_DECOMPOSITION ; direction ++ ) { for ( int id = 0 ; id < count ; id ++ ) { int nbp = ulist_ [ id ] . nbp_ ; if ( nbp > 1 ) { int center = floor ( nbp * 0.5 ) ; for ( int part = center ; part < nbp ; part ++ ) { const int width = ulist_ [ id ] . xmax_ [ direction ] - ulist_ [ id ] . xmin_ [ direction ] + 1 ; const float64 fract = static_cast < float64 > ( part ) / static_cast < float64 > ( nbp ) ; const float64 splitscore = 1. - width / static_cast < float64 > ( dom . nxmax_ [ direction ] ) ; const int tmpwidth = round ( fract * width ) ; const float64 ratioscore = fabs ( fract * width - tmpwidth ) ; const float64 cutscore = ratioscore + splitscore + cutbound * ( tmpwidth < MMAX + 2 ) + cutbound * ( ( width - tmpwidth ) < ( MMAX + 2 ) ) ; if ( cutscore < cutbest ) { cutid = id ; cutbest = cutscore ; cutpart = part ; cutwidth = tmpwidth ; cutdir = direction ; } } } } } if ( cutbest == cutbound ) { printf ( " \n " ) ; } for ( int j = 0 ; j < DIMENSION ; j ++ ) { ulist_ [ count ] . xmax_ [ j ] = ulist_ [ cutid ] . xmax_ [ j ] ; ulist_ [ count ] . xmin_ [ j ] = ulist_ [ cutid ] . xmin_ [ j ] ; } ulist_ [ cutid ] . nbp_ -= cutpart ; ulist_ [ cutid ] . xmin_ [ cutdir ] += cutwidth ; ulist_ [ count ] . nbp_ = cutpart ; ulist_ [ count ] . xmax_ [ cutdir ] = ulist_ [ cutid ] . xmin_ [ cutdir ] - 1 ; } uint64 msum = 0 ; for ( int id = 0 ; id < nbp_ ; id ++ ) { const uint64 mcontrib = ( ulist_ [ id ] . xmax_ [ 3 ] - ulist_ [ id ] . xmin_ [ 3 ] + 1 ) * ( ulist_ [ id ] . xmax_ [ 2 ] - ulist_ [ id ] . xmin_ [ 2 ] + 1 ) * ( ulist_ [ id ] . xmax_ [ 1 ] - ulist_ [ id ] . xmin_ [ 1 ] + 1 ) * ( ulist_ [ id ] . xmax_ [ 0 ] - ulist_ [ id ] . xmin_ [ 0 ] + 1 ) ; if ( master ( ) ) printf ( " \n " , id , ulist_ [ id ] . xmin_ [ 0 ] , ulist_ [ id ] . xmax_ [ 0 ] , ulist_ [ id ] . xmin_ [ 1 ] , ulist_ [ id ] . xmax_ [ 1 ] , ulist_ [ id ] . xmin_ [ 2 ] , ulist_ [ id ] . xmax_ [ 2 ] , ulist_ [ id ] . xmin_ [ 3 ] , ulist_ [ id ] . xmax_ [ 3 ] , mcontrib ) ; msum += mcontrib ; } uint64 mref = static_cast < uint64 > ( dom . nxmax_ [ 3 ] ) * static_cast < uint64 > ( dom . nxmax_ [ 2 ] ) * static_cast < uint64 > ( dom . nxmax_ [ 1 ] ) * static_cast < uint64 > ( dom . nxmax_ [ 0 ] ) ; if ( mref != msum ) { printf ( " \n " , mref , msum ) ; MPI_Barrier ( MPI_COMM_WORLD ) ; MPI_Finalize ( ) ; exit ( 1 ) ; } node_ = & ( ulist_ [ pid_ ] ) ; for ( int j = 0 ; j < DIMENSION ; j ++ ) { dom . local_nx_ [ j ] = node_ -> xmax_ [ j ] - node_ -> xmin_ [ j ] + 1 ; dom . local_nxmin_ [ j ] = node_ -> xmin_ [ j ] ; dom . local_nxmax_ [ j ] = node_ -> xmax_ [ j ] ; } }
Apply the Unbalanced Recursive Bisection (URB) algorithm to establish the boxes
(parts of the computational domain) that are assigned to MPI processes.
Apply the Unbalanced Recursive Bisection ( URB ) algorithm to establish the boxes ( parts of the computational domain ) that are assigned to MPI processes .
// Unbalanced Recursive Bisection // to create the domain decomposition among nbp processes // Best score to choose a cut worse >= 1. / best = 0. // Number of processes chosen for the cut // Id of the process to be cut // Number of points to put into the newly created process // Direction of the cut // Consider cutting only in the first NB_DIMS_DECOMPOSITION dimensions. // Valid values are in-between 1 and 4. // select the best score to cut the former domain //for(int part = center; part < nbp; part++) //if(nbp > 1) //for(int id = 0; id < count; id++) //for(int direction = 0; direction < NB_DIMS_DECOMPOSITION; direction++) //for(int count = 1; count < nbp_; count++) // Update local domain size