balance.h

Go to the documentation of this file.

/** @file balance.h
 */
// Dynamic load-balancing
 
typedef struct {
  short leaf, prolongation;
  int pid;
} NewPid;
 
#define NEWPID() ((NewPid *)&val(newpid,0,0,0))
 
@if TRASH
@ define is_newpid() (!isnan(val(newpid,0,0,0)) && NEWPID()->pid > 0)
@else
@ define is_newpid() (NEWPID()->pid > 0)
@endif
 
Array * linear_tree (size_t size, scalar newpid)
{
  const unsigned short sent = 1 << user, next = 1 << (user + 1);  
  Array * a = array_new();
 
  foreach_cell_post_all (true)
    if (level > 0 && (cell.flags & (sent|next)))
      aparent(0).flags |= next;
 
  bool empty = true;
  foreach_cell_all() {
    if (cell.flags & sent) {
      array_append (a, &cell, size);
      cell.flags &= ~sent;
      empty = false;
    }
    else {
      if (cell.pid >= 0 && NEWPID()->leaf)
  assert (is_leaf(cell));
      if (is_refined_check()) {
  /* check that we are not too refined i.e. that none of our
     children are prolongations */
  bool prolo = false;
  for (int _c = 0; _c < 4; _c++) /* foreach_child */
    if (NEWPID()->prolongation)
      prolo = true;
  if (prolo) {
    /* if we are too refined, we unrefine before sending */
    cell.flags |= leaf;
    array_append (a, &cell, sizeof(Cell));
    cell.flags &= ~leaf;
  }
  else
    array_append (a, &cell, sizeof(Cell));
      }
      else
  array_append (a, &cell, sizeof(Cell));
    }
    if (cell.flags & next)
      cell.flags &= ~next;
    else
      continue;
  }
 
  if (empty)
    a->len = 0;
  return a;
}
 
macro2 foreach_tree (Array * t, size_t size, scalar * list, scalar newpid = newpid)
{
  {
    const unsigned short _sent = 1 << user, _next = 1 << (user + 1);
    scalar * _list = list;
    char * _i = (char *) (t)->p;
    foreach_cell_all() {
      Cell * c = (Cell *) _i;
      if (c->flags & _sent) {
  _i += size;
  {...}
      }
      else
  _i += sizeof(Cell);
      if (c->flags & _next) {
  assert (c->neighbors);
  if (!(c->flags & leaf) && is_leaf(cell) &&
      (!is_newpid() || !NEWPID()->leaf))
    /* refined */
    refine_cell (point, _list, 0, NULL);
  else if (!cell.neighbors)
    /* prolongation */
    alloc_children (point);
      }
      else
  continue;
    }
  }
}
 
Array * neighborhood (scalar newpid, int nextpid, FILE * fp)
{
  const unsigned short sent = 1 << user;
  for (int _i = 0; _i < _N; _i++) /* foreach_cell */ {
    // root cells
    bool root = false;
    if ((!is_local(cell) || NEWPID()->pid - 1 != nextpid) && is_refined(cell)) {
      for (int _c = 0; _c < 4; _c++) /* foreach_child */
  if (is_local(cell) && NEWPID()->pid - 1 == nextpid) {
    root = true; break;
  }
      if (root && cell.pid != nextpid) {
  for (int _n = 0; _n < ; _n++)
    if (cell.pid != nextpid && is_newpid()) {
      if (fp)
        fprintf (fp, "%g %g %g %d %d root\n",
           x, y, z, NEWPID()->pid - 1, cell.pid);
      cell.flags |= sent;
    }
      }
    }
    // children
    if ((is_local(cell) && NEWPID()->pid - 1 == nextpid) || root) {
      for (int _n = 0; _n < 1; _n++)
  if (cell.neighbors && cell.pid != nextpid)
    for (int _c = 0; _c < 4; _c++) /* foreach_child */
      if (cell.pid != nextpid && is_newpid()) {
        if (fp)
    fprintf (fp, "%g %g %g %d %d nextpid\n",
       x, y, z, NEWPID()->pid - 1, cell.pid);
        cell.flags |= sent;
      }
    }
    if (is_leaf(cell))
      continue;
  }
 
  return linear_tree (sizeof(Cell) + datasize, newpid);
}
 
static void send_tree (Array * a, int to, MPI_Request * r)
{
  MPI_Isend (&a->len, 1, MPI_LONG, to, MOVED_TAG(), MPI_COMM_WORLD, &r[0]);
  if (a->len > 0) {
    MPI_Isend (a->p, a->len, MPI_BYTE, to, MOVED_TAG(), MPI_COMM_WORLD, &r[1]);
    tree->dirty = true;
  }
}
 
static void receive_tree (int from, scalar newpid, FILE * fp)
{
  Array a;
  mpi_recv_check (&a.len, 1, MPI_LONG, from, MOVED_TAG(),
      MPI_COMM_WORLD, MPI_STATUS_IGNORE, "receive_tree (len)");
  if (a.len > 0) {
    a.p = malloc (a.len);
    if (fp)
      fprintf (fp, "receiving %ld from %d\n", a.len, from);
    mpi_recv_check (a.p, a.len, MPI_BYTE, from, MOVED_TAG(),
        MPI_COMM_WORLD, MPI_STATUS_IGNORE, "receive_tree (p)");
    //    const unsigned short next = 1 << (user + 1);
    foreach_tree (&a, sizeof(Cell) + datasize, NULL) {
      memcpy (((char *)&cell) + sizeof(Cell), ((char *)c) + sizeof(Cell),
        datasize);
      assert (NEWPID()->pid > 0);
      if (fp)
  fprintf (fp, "%g %g %g %d %d %d %d %d %d recv\n",
     x, y, z, NEWPID()->pid - 1, cell.pid,
     c->flags & leaf,
     cell.flags & leaf, from, NEWPID()->leaf);
    }
    free (a.p);
    tree->dirty = true;
  }
}
 
static void wait_tree (Array * a, MPI_Request * r)
{
  MPI_Wait (&r[0], MPI_STATUS_IGNORE);
  if (a->len > 0)
    MPI_Wait (&r[1], MPI_STATUS_IGNORE);
}
 
static void check_flags()
{
#if DEBUG_MPI
  for (int _i = 0; _i < _N; _i++) /* foreach_cell */
    for (int _n = 0; _n < ; _n++)
      if (allocated(0))
  for (int i = user; i <= user + 7; i++)
    assert (!(cell.flags & (1 << i)));
#endif
}
 
struct {
  int  min;    // minimum number of points per process
  bool leaves; // balance leaves only
  
  int npe; // number of active processes
} mpi = {
  1,
  true
};
 
trace
bool balance()
{
  if (npe() == 1)
    return false;
 
  assert (sizeof(NewPid) == sizeof(double));
 
  check_flags();
 
  long nl = 0, nt = 0;
  for (int _i = 0; _i < _N; _i++) /* foreach_cell */ {
    if (is_local(cell)) {
      nt++;
      if (is_leaf(cell))
  nl++;
    }
    if (is_leaf(cell))
      continue;
  }
 
  grid->n = grid->tn = nl;
  grid->maxdepth = depth();
  long nmin = nl, nmax = nl;
  // fixme: do all reductions in one go
  mpi_all_reduce (nmax, MPI_LONG, MPI_MAX);
  mpi_all_reduce (nmin, MPI_LONG, MPI_MIN);
  mpi_all_reduce (grid->tn, MPI_LONG, MPI_SUM);
  mpi_all_reduce (grid->maxdepth, MPI_INT, MPI_MAX);
  if (mpi.leaves)
    nt = grid->tn;
  else
    mpi_all_reduce (nt, MPI_LONG, MPI_SUM);
    
  long ne = max(1, nt/npe());
 
  if (ne < mpi.min) {
    mpi.npe = max(1, nt/mpi.min);
    ne = max(1, nt/mpi.npe);
  }
  else
    mpi.npe = npe();
 
  if (nmax - nmin <= 1)
    return false;
  
  scalar newpid[];
  double zn = z_indexing (newpid, mpi.leaves);
  if (pid() == 0)
    assert (zn + 1 == nt);
  
  FILE * fp = NULL;
#ifdef DEBUGCOND
  extern double t;
  if (DEBUGCOND)
    fp = lfopen ("bal", "w");
#elif DEBUG_MPI
  fp = lfopen ("bal", "w");
#endif
 
  // compute new pid, stored in newpid[]
  bool next = false, prev = false;
  foreach_cell_all() {
    if (is_local(cell)) {
      int pid = balanced_pid (newpid[], nt, mpi.npe);
      pid = clamp (pid, cell.pid - 1, cell.pid + 1);
      if (pid == pid() + 1)
  next = true;
      else if (pid == pid() - 1)
  prev = true;
      NEWPID()->pid = pid + 1;
      NEWPID()->leaf = is_leaf(cell);
      NEWPID()->prolongation = is_prolongation(cell);
      if (fp)
  fprintf (fp, "%g %g %d %d newpid\n", x, y, NEWPID()->pid - 1, cell.pid);
    }
    else
      newpid[] = 0;
  }
  for (int l = 0; l <= depth(); l++)
    boundary_iterate (level, {newpid}, l);
 
#ifdef DEBUGCOND
  extern double t;
  NOT_UNUSED(t);
  if (DEBUGCOND) {
    char name[80];
    sprintf (name, "colls-before-%d", pid());
    FILE * fp = fopen (name, "w");
    output_cells (fp);
    fclose (fp);
 
    sprintf (name, "pid-before-%d", pid());
    fp = fopen (name, "w");
    for (int _i = 0; _i < _N; _i++) /* foreach_cell */ {
      fprintf (fp, "%g %g %g %d %d %d %d\n",
         x, y, z, cell.pid, NEWPID()->pid - 1,
         NEWPID()->leaf, is_leaf(cell));
      if (NEWPID()->leaf)
  assert (is_leaf(cell));
    }
    fclose (fp);  
  }
#endif // DEBUGCOND
  
  Array * anext = next ? neighborhood (newpid, pid() + 1, fp) : array_new();
  Array * aprev = prev ? neighborhood (newpid, pid() - 1, fp) : array_new();
 
  if (fp)
    fflush (fp);
  
  check_flags();
  
  // send mesh to previous/next process
  MPI_Request rprev[2], rnext[2];
  if (pid() > 0)
    send_tree (aprev, pid() - 1, rprev);
  if (pid() < npe() - 1)
    send_tree (anext, pid() + 1, rnext);
 
  // receive mesh from next/previous process
  if (pid() < npe() - 1)
    receive_tree (pid() + 1, newpid, fp);
  if (pid() > 0)
    receive_tree (pid() - 1, newpid, fp);
 
  /* check that mesh was received OK and free send buffers */
  if (pid() > 0)
    wait_tree (aprev, rprev);
  array_free (aprev);
  if (pid() < npe() - 1)
    wait_tree (anext, rnext);
  array_free (anext);
 
  if (fp)
    fflush (fp);
  
  // set new pids
  int pid_changed = false;
  foreach_cell_all() {
    if (cell.pid >= 0) {
      if (is_newpid()) {
  if (fp)
    fprintf (fp, "%g %g %g %d %d %d %d %d new\n",
       x, y, z, NEWPID()->pid - 1, cell.pid,
       is_leaf(cell), cell.neighbors, NEWPID()->leaf);
  if (cell.pid != NEWPID()->pid - 1) {
    cell.pid = NEWPID()->pid - 1;
    cell.flags &= ~(active|border);
    if (is_local(cell))
      cell.flags |= active;
    pid_changed = true;
  }
  if (NEWPID()->leaf && !is_leaf(cell) && cell.neighbors)
    coarsen_cell_recursive (point, NULL);
      }
      else if (level > 0 && ((NewPid *)&coarse(newpid))->leaf)
  cell.pid = aparent(0).pid;
    }
    // cleanup unused prolongations
    if (!cell.neighbors && allocated_child(0)) {
      if (fp)
  fprintf (fp, "%g %g %g %d %d freechildren\n",
     x, y, z, NEWPID()->pid - 1, cell.pid);
      free_children (point);
    }
  }
 
  if (tree->dirty || pid_changed) {
#if 1
    // update active cells: fixme: can this be done above
    foreach_cell_post (!is_leaf (cell))
      if (!is_leaf(cell) && !is_local(cell)) {
  unsigned short flags = cell.flags & ~active;
  for (int _c = 0; _c < 4; _c++) /* foreach_child */
    if (is_active(cell)) {
      flags |= active; break;
    }
  cell.flags = flags;
      }
#endif
    flag_border_cells(); // fixme: can this be done above?
    pid_changed = true;
  }
  
  if (fp)
    fclose (fp);
 
  mpi_all_reduce (pid_changed, MPI_INT, MPI_MAX);
  if (pid_changed)
    mpi_boundary_update_buffers();
  
  return pid_changed;
}
 
void mpi_boundary_update (scalar * list)
{
  mpi_boundary_update_buffers();
  for (int _s = 0; _s < 1; _s++) /* scalar in list */
    set_dirty_stencil (s);
  grid->tn = 0; // so that tree is not "full" for the call below
  boundary (list);
  while (balance());
}