CarpetLib: Begin to add new load balancing mechanism

Add new class bintree<T,D,P>.
Add new routine balance<T>.

7 files changed, 1237 insertions, 1 deletions

diff --git a/Carpet/CarpetLib/src/balance.cc b/Carpet/CarpetLib/src/balance.cc
new file mode 100644
index 000000000..e9f5cad6a
--- /dev/null
+++ b/Carpet/CarpetLib/src/balance.cc
@@ -0,0 +1,495 @@
+#include <cctk.h>
+
+#include <algorithm>
+#include <cassert>
+#include <cstdlib>
+#include <limits>
+#include <list>
+#include <vector>
+
+#include "region.hh"
+
+#include "balance.hh"
+
+using namespace std;
+
+
+
+namespace CarpetLib {
+  
+  // Interface for one item
+  struct item_ifc {
+    CCTK_REAL load() const;
+    item_ifc split (CCTK_REAL ratio_new_over_old);
+  };
+  
+  // Declaration of template functions implementing an equivalen
+  // interface
+  template <typename item_t>
+  CCTK_REAL
+  item_load (item_t const& item);
+  
+  template <typename item_t>
+  item_t item_split (item_t& item, CCTK_REAL ratio_new_over_old);
+  
+  // Default definitions of these template functions using the
+  // interface above
+  template <typename item_t>
+  CCTK_REAL
+  item_load (item_t const& item)
+  {
+    return item.load();
+  }
+  
+  template <typename item_t>
+  item_t item_split (item_t& item, CCTK_REAL ratio_new_over_old)
+  {
+    return item.split (ratio_new_over_old);
+  }
+  
+  
+  
+  // A collection of items
+  template <typename item_t>
+  class items_t {
+    typedef list<item_t> coll_t;
+    coll_t items;
+    typename coll_t::iterator find_largest_item ();
+  public:
+    items_t ();
+    items_t (vector<item_t> const& items_);
+    void add_item (item_t const& item);
+    bool empty () const;
+    size_t size () const;
+    CCTK_REAL load () const;
+    item_t& get_one_item ();
+    item_t& get_largest_item ();
+    item_t get_and_remove_largest_item ();
+    void copy_out (vector<item_t>& dst) const;
+  };
+  
+  // A worker
+  template <typename item_t>
+  class worker_t: public items_t<item_t> {
+  public:
+    CCTK_REAL strength () const;
+    CCTK_REAL business () const;
+  };
+  
+  // A collection of workers
+  template <typename item_t>
+  class workers_t {
+    typedef vector<worker_t<item_t> > coll_t;
+    coll_t workers;
+    typename coll_t::iterator find_least_busy_worker ();
+    typename coll_t::iterator find_most_busy_worker ();
+  public:
+    workers_t (int nworkers);
+    bool empty () const;
+    size_t size () const;
+    CCTK_REAL load () const;
+    CCTK_REAL strength () const;
+    CCTK_REAL ideal_business () const;
+    CCTK_REAL imbalance () const;
+    worker_t<item_t>& get_least_busy_worker ();
+    worker_t<item_t>& get_most_busy_worker ();
+    void ensure_same_size ();
+    void copy_out (vector<vector<item_t> >& dst) const;
+  };
+  
+  
+  
+  template <typename item_t>
+  items_t<item_t>::items_t ()
+  {
+  }
+  
+  template <typename item_t>
+  items_t<item_t>::items_t (vector<item_t> const& items_)
+  {
+    for (typename vector<item_t>::const_iterator
+           p = items_.begin(); p != items_.end(); ++p)
+    {
+      add_item (*p);
+    }
+  }
+  
+  template <typename item_t>
+  void
+  items_t<item_t>::add_item (item_t const& item)
+  {
+    items.push_back (item);
+  }
+  
+  template <typename item_t>
+  bool
+  items_t<item_t>::empty () const
+  {
+    return items.empty();
+  }
+  
+  template <typename item_t>
+  size_t
+  items_t<item_t>::size () const
+  {
+    return items.size();
+  }
+  
+  template <typename item_t>
+  CCTK_REAL
+  items_t<item_t>::load () const
+  {
+    CCTK_REAL total_load = 0.0;
+    for (typename coll_t::const_iterator
+           p = items.begin(); p != items.end(); ++p)
+    {
+      total_load += item_load (*p);
+    }
+    return total_load;
+  }
+  
+  template <typename item_t>
+  typename items_t<item_t>::coll_t::iterator
+  items_t<item_t>::find_largest_item ()
+  {
+    typename coll_t::iterator max_item = items.end();
+    CCTK_REAL max_load = -1.0;
+    for (typename coll_t::iterator p = items.begin(); p != items.end(); ++p) {
+      if (item_load (*p) > max_load) {
+        max_item = p;
+        max_load = item_load (*max_item);
+      }
+  }
+    return max_item;
+  }
+  
+  template <typename item_t>
+  item_t&
+  items_t<item_t>::get_one_item ()
+  {
+    assert (not empty());
+    return items.front();
+  }
+  
+  template <typename item_t>
+  item_t&
+  items_t<item_t>::get_largest_item ()
+  {
+    typename coll_t::iterator const max_item = find_largest_item();
+    assert (max_item != items.end());
+    return *max_item;
+  }
+  
+  template <typename item_t>
+  item_t
+  items_t<item_t>::get_and_remove_largest_item ()
+  {
+    typename coll_t::iterator const max_item = find_largest_item();
+    assert (max_item != items.end());
+    item_t const item = *max_item;
+    items.erase (max_item);
+    return item;
+  }
+  
+  template <typename item_t>
+  void
+  items_t<item_t>::copy_out (vector<item_t>& dst) const
+  {
+    dst.resize (items.size());
+    copy (items.begin(), items.end(), dst.begin());
+  }
+  
+  
+  
+  template <typename item_t>
+  CCTK_REAL
+  worker_t<item_t>::strength () const
+  {
+    return 1.0;                 // All workers have the same strength
+  }
+  
+  template <typename item_t>
+  CCTK_REAL
+  worker_t<item_t>::business () const
+  {
+    return this->load() / strength();
+  }
+  
+  
+  
+  template <typename item_t>
+  workers_t<item_t>::workers_t (int const nworkers)
+    : workers (nworkers)
+  {
+  }
+  
+  template <typename item_t>
+  bool
+  workers_t<item_t>::empty () const
+  {
+    return workers.empty();
+  }
+  
+  template <typename item_t>
+  size_t
+  workers_t<item_t>::size () const
+  {
+    return workers.size();
+  }
+  
+  template <typename item_t>
+  CCTK_REAL
+  workers_t<item_t>::load () const
+  {
+    CCTK_REAL total_load = 0.0;
+    for (typename coll_t::const_iterator
+           w = workers.begin(); w != workers.end(); ++w)
+    {
+      total_load += w->load();
+    }
+    return total_load;
+  }
+  
+  template <typename item_t>
+  CCTK_REAL
+  workers_t<item_t>::strength () const
+  {
+    CCTK_REAL total_strength = 0.0;
+    for (typename coll_t::const_iterator
+           w = workers.begin(); w != workers.end(); ++w)
+    {
+      total_strength += w->strength();
+    }
+    return total_strength;
+  }
+  
+  template <typename item_t>
+  CCTK_REAL
+  workers_t<item_t>::ideal_business () const
+  {
+    return load() / strength();
+  }
+  
+  template <typename item_t>
+  CCTK_REAL
+  workers_t<item_t>::imbalance () const
+  {
+    assert (not empty());
+    CCTK_REAL max_load = 0.0;
+    CCTK_REAL avg_load = 0.0;
+    for (typename coll_t::const_iterator
+           w = workers.begin(); w != workers.end(); ++w)
+    {
+      max_load = max (max_load, w->business());
+      avg_load += w->business();
+    }
+    avg_load /= size();
+    return max_load - avg_load;
+  }
+  
+  template <typename item_t>
+  typename workers_t<item_t>::coll_t::iterator
+  workers_t<item_t>::find_least_busy_worker ()
+  {
+    typename coll_t::iterator min_worker = workers.end();
+    CCTK_REAL min_business = numeric_limits<CCTK_REAL>::max();
+    for (typename coll_t::iterator
+           w = workers.begin(); w != workers.end(); ++w)
+    {
+      if (w->business() < min_business) {
+        min_worker = w;
+        min_business = min_worker->business();
+      }
+    }
+    return min_worker;
+  }
+  
+  template <typename item_t>
+  typename workers_t<item_t>::coll_t::iterator
+  workers_t<item_t>::find_most_busy_worker ()
+  {
+    typename coll_t::iterator max_worker = workers.end();
+    CCTK_REAL max_business = 0.0;
+    for (typename coll_t::iterator
+           w = workers.begin(); w != workers.end(); ++w)
+    {
+      if (w->business() > max_business) {
+        max_worker = w;
+        max_business = max_worker->business();
+      }
+    }
+    return max_worker;
+  }
+  
+  template <typename item_t>
+  worker_t<item_t>&
+  workers_t<item_t>::get_least_busy_worker ()
+  {
+    typename coll_t::iterator const min_worker = find_least_busy_worker();
+    assert (min_worker != workers.end());
+    return *min_worker;
+  }
+  
+  template <typename item_t>
+  worker_t<item_t>&
+  workers_t<item_t>::get_most_busy_worker ()
+  {
+    typename coll_t::iterator const max_worker = find_most_busy_worker();
+    assert (max_worker != workers.end());
+    return *max_worker;
+  }
+  
+  template <typename item_t>
+  void
+  workers_t<item_t>::ensure_same_size ()
+  {
+    if (empty()) return;          // nothing to do
+    
+    size_t max_items = 0;
+    typename coll_t::iterator nonempty_worker = workers.end();
+    for (typename coll_t::iterator
+           w = workers.begin(); w != workers.end(); ++w)
+    {
+      if (nonempty_worker == workers.end() and not w->empty()) {
+        nonempty_worker = w;
+      }
+      max_items = max (max_items, w->size());
+    }
+    if (max_items == 0) return;   // load is already equal
+    
+    for (typename coll_t::iterator
+           w = workers.begin(); w != workers.end(); ++w)
+    {
+      // find a worker who has an item
+      typename coll_t::iterator const worker = w->empty() ? nonempty_worker : w;
+      while (w->size() < max_items) {
+        CCTK_REAL const ratio = 0.0; // create empty fill-up items
+        w->add_item (item_split (w->get_one_item(), ratio));
+      }
+    }
+    
+    for (typename coll_t::const_iterator
+           w = workers.begin(); w != workers.end(); ++w)
+    {
+      assert (w->size() == max_items);
+    }
+  }
+  
+  template <typename item_t>
+  void
+  workers_t<item_t>::copy_out (vector<vector<item_t> >& dst) const
+  {
+    dst.resize (workers.size());
+    for (size_t w = 0; w < workers.size(); ++w) {
+      workers.at(w).copy_out (dst.at(w));
+    }
+  }
+  
+  
+  
+  //////////////////////////////////////////////////////////////////////////////
+  
+  
+  
+  template <typename item_t>
+  void
+  assign_item (items_t<item_t>& items, workers_t<item_t>& workers)
+  {
+    // Assign the largest item to the least busy worker
+    item_t const item = items.get_and_remove_largest_item ();
+    worker_t<item_t>& worker = workers.get_least_busy_worker ();
+    worker.add_item (item);
+  }
+  
+  template <typename item_t>
+  void
+  split_and_distribute (workers_t<item_t>& workers)
+  {
+    // Split the largest item of the most busy worker and give the
+    // remainder to another worker
+    worker_t<item_t>& worker = workers.get_most_busy_worker ();
+    item_t&           item = worker.get_largest_item ();
+    
+    // Determine how to split the item
+    CCTK_REAL const imbalance = worker.business() - workers.ideal_business();
+    // Should we even be here?
+    assert (imbalance > 0.0);
+    CCTK_REAL const item_business = item_load (item) / worker.strength();
+    // This should be the largest item!
+    assert (item_business > 0.0);
+    // Determine how much of the item to give away
+    CCTK_REAL const ratio = imbalance / item_business;
+    // A ratio of one or more indicates that the item should be given
+    // away in its entirety -- which would mean that something went
+    // wrong before
+    assert (ratio < 1.0);
+    
+    // Split the item...
+    item_t const new_item = item_split (item, ratio);
+    // ...and give the remainder to someone else
+    worker_t<item_t>& new_worker = workers.get_least_busy_worker ();
+    // This should be someone else!
+    assert (&new_worker != &worker);
+    new_worker.add_item (new_item);
+  }
+  
+  
+  
+  template <typename item_t>
+  void
+  balance (vector<item_t> const& items_,
+           vector<vector<item_t> >& split_items,
+           int const nworkers,
+           CCTK_REAL const max_imbalance,
+           bool const ensure_same_size)
+  {
+    assert (max_imbalance >= 0.0);
+    
+    if (items_.empty()) return; // nothing to do
+    items_t<item_t> items (items_);
+    assert (nworkers > 0);
+    workers_t<item_t> workers (nworkers);
+    
+    // Assign items
+    while (not items.empty()) {
+      assign_item (items, workers);
+    }
+    
+    // TODO: To parallelise this: Group workers into groups, and
+    // assign work to these groups.  Then balance the load
+    // recursively.
+    
+    // Balance items
+    for (;;) {
+      // Measure imbalance
+      CCTK_REAL const imbalance = workers.imbalance();
+      // Are we done?
+      if (imbalance <= max_imbalance) break;
+      // Should we even be here?
+      assert (workers.size() > 1);
+      
+      // Do something
+      split_and_distribute (workers);
+      
+      // Ensure progress
+      assert (workers.imbalance() < imbalance);
+    }
+    
+    if (ensure_same_size) {
+      workers.ensure_same_size ();
+    }
+    
+    workers.copy_out (split_items);
+  }
+  
+  
+  
+  template
+  void
+  balance (vector<region_t> const& items_,
+           vector<vector<region_t> >& split_items,
+           int const nworkers,
+           CCTK_REAL const max_imbalance,
+           bool const ensure_same_size);
+  
+} // namespace CarpetLib
diff --git a/Carpet/CarpetLib/src/balance.hh b/Carpet/CarpetLib/src/balance.hh
new file mode 100644
index 000000000..ccf7dad91
--- /dev/null
+++ b/Carpet/CarpetLib/src/balance.hh
@@ -0,0 +1,29 @@
+#ifndef BALANCE_HH
+#define BALANCE_HH
+
+#include <cctk.h>
+
+#include <vector>
+
+using namespace std;
+
+
+
+namespace CarpetLib {
+  
+  // This routine splits N items over P workers.  It can split the
+  // items if necessary to ensure a maximum imbalance, and it can
+  // ensure (artificially) that each worker receives the same number
+  // of items.
+
+  template <typename T>
+  void
+  balance (vector<T> const& items_,
+           vector<vector<T> >& split_items_,
+           int nworkers,
+           CCTK_REAL max_imbalance,
+           bool ensure_same_size);
+  
+} // namespace CarpetLib
+
+#endif  // #ifndef BALANCE_HH
diff --git a/Carpet/CarpetLib/src/bintree.cc b/Carpet/CarpetLib/src/bintree.cc
new file mode 100644
index 000000000..db3b27e1c
--- /dev/null
+++ b/Carpet/CarpetLib/src/bintree.cc
@@ -0,0 +1,489 @@
+#include <algorithm>
+#include <cmath>
+#include <iostream>
+
+#include "defs.hh"
+#include "region.hh"
+
+#include "bintree.hh"
+
+
+
+// Create an empty tree
+template <typename T, int D, typename P>
+bintree<T,D,P>::bintree ()
+  : type (type_empty)
+{
+  assert (invariant());
+  // This is unused
+  assert (0);
+}
+
+
+
+// Create a tree branch from a list of bounds and subtrees
+template <typename T, int D, typename P>
+bintree<T,D,P>::bintree (int const dir_, vect <T, 3> const & bounds_,
+                         vect <bintree *, 2> const & subtrees_)
+  : type (type_branch), dir (dir_), bounds (bounds_), subtrees (subtrees_)
+{
+  assert (dir>=0 and dir<D);
+  assert (invariant());
+}
+
+
+
+// Create a tree leaf from a payload
+template <typename T, int D, typename P>
+bintree<T,D,P>::bintree (P const & p_)
+  : type (type_leaf), p (p_)
+{
+  assert (invariant());
+}
+  
+
+
+// Create a tree as copy from another tree
+template <typename T, int D, typename P>
+bintree<T,D,P>::bintree (bintree const & t)
+  : type (t.type)
+{
+  switch (type) {
+  case type_empty:
+    // do nothing
+    break;
+  case type_branch:
+    dir = t.dir;
+    bounds = t.bounds;
+    for (int i=0; i<2; ++i) {
+      subtrees[i] = new bintree (*t.subtrees[i]);
+    }
+    break;
+  case type_leaf:
+    p = t.p;
+    break;
+  default:
+    assert (0);
+  }
+  assert (invariant());
+}
+  
+
+
+// Copy a tree
+template <typename T, int D, typename P>
+bintree<T,D,P> &
+bintree<T,D,P>::operator= (bintree const & t)
+{
+  assert (invariant());
+  if (is_branch()) {
+    for (int i=0; i<2; ++i) {
+      delete subtrees[i];
+    }
+  }
+  type = t.type;
+  switch (type) {
+  case type_empty:
+    // do nothing
+    break;
+  case type_branch:
+    dir = t.dir;
+    bounds = t.bounds;
+    for (int i=0; i<2; ++i) {
+      subtrees[i] = new bintree (*t.subtrees[i]);
+    }
+    break;
+  case type_leaf:
+    p = t.p;
+    break;
+  default:
+    assert (0);
+  }
+  assert (invariant());
+  return *this;
+}
+
+
+
+// Delete a tree (and its subtrees)
+template <typename T, int D, typename P>
+bintree<T,D,P>::~bintree ()
+{
+  assert (invariant());
+  if (is_branch()) {
+    for (int i=0; i<2; ++i) {
+      delete subtrees[i];
+    }
+  }
+}
+
+
+
+// Compare trees
+template <typename T, int D, typename P>
+bool
+bintree<T,D,P>::operator== (bintree const & t) const
+{
+  assert (invariant());
+  assert (t.invariant());
+  if (type != t.type) return false;
+  switch (type) {
+  case type_empty:
+    break;
+  case type_branch:
+    if (dir != t.dir) return false;
+    if (any (bounds != t.bounds)) return false;
+    if (any (subtrees != t.subtrees)) return false;
+    break;
+  case type_leaf:
+    return p == t.p;
+  default:
+    assert (0);
+  }
+  return true;
+}
+
+
+
+// Invariant
+template <typename T, int D, typename P>
+bool
+bintree<T,D,P>::invariant () const
+{
+  return empty() + is_branch() + is_leaf() == 1;
+}
+
+
+
+// Find the leaf payload corresponding to a position
+template <typename T, int D, typename P>
+P const *
+bintree<T,D,P>::search (tvect const & ipos) const
+{
+  assert (not empty());
+  if (is_leaf()) return & p;
+  int const i = ::asearch (ipos[dir], bounds);
+  if (i<0 or i>=2) return NULL; // not found
+  return subtrees[i]->search(ipos);
+}
+
+template <typename T, int D, typename P>
+P *
+bintree<T,D,P>::search (tvect const & ipos)
+{
+  assert (not empty());
+  if (is_leaf()) return & p;
+  int const i = ::asearch (ipos[dir], bounds);
+  if (i<0 or i>=2) return NULL; // not found
+  return subtrees[i]->search(ipos);
+}
+
+
+
+// Const iterator
+template <typename T, int D, typename P>
+bintree<T,D,P>::const_iterator::const_iterator (bintree const & f_)
+  : f(f_), i(0), it(0)
+{
+  if (f.is_branch()) {
+    it = new const_iterator (* f.subtrees[i]);
+    while ((*it).done()) {
+      delete it;
+      it = 0;
+      ++ i;
+      if (done()) break;
+      // to do: use a new function "reset iterator" instead
+      it = new const_iterator (* f.subtrees[i]);
+    }
+    assert (done() or not (*it).done());
+  }
+}
+
+template <typename T, int D, typename P>
+bintree<T,D,P>::const_iterator::const_iterator (bintree const & f_, int)
+  : f(f_), it(0)
+{
+  if (f.empty()) {
+    i = 0;
+  } else if (f.is_leaf()) {
+    i = 1;
+  } else {
+    i = 2;
+  }
+}
+
+template <typename T, int D, typename P>
+bintree<T,D,P>::const_iterator::~const_iterator ()
+{
+  if (it) {
+    delete it;
+  }
+}
+
+template <typename T, int D, typename P>
+bintree<T,D,P> const & 
+bintree<T,D,P>::const_iterator::operator* () const
+{
+  assert (not done());
+  assert (not f.empty());
+  if (f.is_leaf()) {
+    return f;
+  } else {
+    assert (it);
+    return **it;
+  }
+}
+
+template <typename T, int D, typename P>
+bool
+bintree<T,D,P>::const_iterator::operator== (const_iterator const & it2) const
+{
+  // assert (f == it2.f);
+  assert (&f == &it2.f);
+  if (i != it2.i) return false;
+  if (it == 0 and it2.it == 0) return true;
+  if (it == 0 or it2.it == 0) return false;
+  return *it == *it2.it;
+}
+
+template <typename T, int D, typename P>
+typename bintree<T,D,P>::const_iterator &
+bintree<T,D,P>::const_iterator::operator++ ()
+{
+  assert (not done());
+  assert (not f.empty());
+  if (f.is_leaf()) {
+    ++ i;
+  } else {
+    ++ *it;
+    while ((*it).done()) {
+      delete it;
+      it = 0;
+      ++ i;
+      if (done()) break;
+      // to do: use a new function "reset iterator" instead
+      it = new const_iterator (* f.subtrees[i]);
+    }
+    assert (done() or not (*it).done());
+  }
+  return *this;
+}
+
+template <typename T, int D, typename P>
+bool
+bintree<T,D,P>::const_iterator::done () const
+{
+  if (f.empty()) {
+    return true;
+  } else if (f.is_leaf()) {
+    return i > 0;
+  } else {
+    return i == 2;
+  }
+}
+
+
+
+// Non-const iterator
+template <typename T, int D, typename P>
+bintree<T,D,P>::iterator::iterator (bintree & f_)
+  : f(f_), i(0), it(0)
+{
+  if (f.is_branch()) {
+    it = new iterator (* f.subtrees[i]);
+    while ((*it).done()) {
+      delete it;
+      it = 0;
+      ++ i;
+      if (done()) break;
+      // to do: use a new function "reset iterator" instead
+      it = new iterator (* f.subtrees[i]);
+    }
+    assert (done() or not (*it).done());
+  }
+}
+
+template <typename T, int D, typename P>
+bintree<T,D,P>::iterator::iterator (bintree & f_, int)
+  : f(f_), it(0)
+{
+  if (f.empty()) {
+    i = 0;
+  } else if (f.is_leaf()) {
+    i = 1;
+  } else {
+    i = 2;
+  }
+}
+
+template <typename T, int D, typename P>
+bintree<T,D,P>::iterator::~iterator ()
+{
+  if (it) {
+    delete it;
+  }
+}
+
+template <typename T, int D, typename P>
+bintree<T,D,P> & 
+bintree<T,D,P>::iterator::operator* () const
+{
+  assert (not done());
+  assert (not f.empty());
+  if (f.is_leaf()) {
+    return f;
+  } else {
+    assert (it);
+    return **it;
+  }
+}
+
+template <typename T, int D, typename P>
+bool
+bintree<T,D,P>::iterator::operator== (iterator const & it2) const
+{
+  // assert (f == it2.f);
+  assert (&f == &it2.f);
+  if (i != it2.i) return false;
+  if (it == 0 and it2.it == 0) return true;
+  if (it == 0 or it2.it == 0) return false;
+  return *it == *it2.it;
+}
+
+template <typename T, int D, typename P>
+typename bintree<T,D,P>::iterator &
+bintree<T,D,P>::iterator::operator++ ()
+{
+  assert (not done());
+  assert (not f.empty());
+  if (f.is_leaf()) {
+    ++ i;
+  } else {
+    ++ *it;
+    while ((*it).done()) {
+      delete it;
+      it = 0;
+      ++ i;
+      if (done()) break;
+      // to do: use a new function "reset iterator" instead
+      it = new iterator (* f.subtrees[i]);
+    }
+    assert (done() or not (*it).done());
+  }
+  return *this;
+}
+
+template <typename T, int D, typename P>
+bool
+bintree<T,D,P>::iterator::done () const
+{
+  if (f.empty()) {
+    return true;
+  } else if (f.is_leaf()) {
+    return i > 0;
+  } else {
+    return i == 2;
+  }
+}
+
+
+
+// Memory usage
+template <typename T, int D, typename P>
+size_t
+bintree<T,D,P>::memory () const
+{
+  size_t size = sizeof *this;
+  if (is_branch()) {
+    for (int i=0; i<2; ++i) {
+      size += memoryof(*subtrees[i]);
+    }
+  }
+  return size;
+}
+
+
+
+// Output helper
+template <typename T, int D, typename P>
+void
+bintree<T,D,P>::output (ostream & os) const
+{
+  os << "bintree{";
+  if (empty()) {
+    os << "empty";
+  } else if (is_branch()) {
+    os << "branch:"
+       << "dir=" << dir << ","
+       << "subtrees=[";
+    for (int i=0; i<2; ++i) {
+      os << bounds[i] << ":[" << i << "]=" << *subtrees[i] << ":";
+    }
+    os << bounds[2] << "]";
+  } else {
+    os << "leaf:"
+       << "payload=" << p;
+  }
+  os << "}";
+}
+
+
+
+// Generic arithmetic search
+template <typename T, int D>
+static
+int asearch (T const t, vect <T, D> const & ts)
+{
+  int imin = 0;
+  int imax = D - 1;
+  if (imax <= imin) {
+    return imin;
+  }
+  T tmin = ts[imin];
+  if (t < tmin) {
+    return -1;
+  }
+  T tmax = ts[imax];
+  if (t >= tmax) {
+    return imax;
+  }
+  int isize = imax - imin;
+  for (;;) {
+    // check loop invariants
+    assert (imin < imax);
+    assert (t >= tmin);
+    assert (t < tmax);
+    if (imax == imin+1) {
+      return imin;
+    }
+    assert (tmax > tmin);       // require that ts is strictly ordered
+    CCTK_REAL const rguess =
+      (imax - imin) * CCTK_REAL(t - tmin) / (tmax - tmin);
+    int const iguess = imin + max (1, int(floor(rguess)));
+    // handle round-off errors
+    if (iguess == imax) {
+      return imax - 1;
+    }
+    assert (iguess > imin and iguess < imax);
+    T const tguess = ts[iguess];
+    if (t < tguess) {
+      imax = iguess;
+      tmax = tguess;
+    } else {
+      imin = iguess;
+      tmin = tguess;
+    }
+    // check loop variant
+    int const newisize = imax - imin;
+    assert (newisize < isize);
+    isize = newisize;
+  }
+}
+
+
+
+template class bintree <int, dim, pseudoregion_t>;
+
+template size_t memoryof (bintree <int, dim, pseudoregion_t> const & f);
+
+template ostream & operator<< (ostream & os, bintree <int, dim, pseudoregion_t> const & f);
diff --git a/Carpet/CarpetLib/src/bintree.hh b/Carpet/CarpetLib/src/bintree.hh
new file mode 100644
index 000000000..7a7927d9a
--- /dev/null
+++ b/Carpet/CarpetLib/src/bintree.hh
@@ -0,0 +1,172 @@
+#ifndef BINTREE_HH
+#define BINTREE_HH
+
+#include <cassert>
+#include <cmath>
+#include <cstdlib>
+#include <iostream>
+#include <vector>
+
+#include <vect.hh>
+
+using namespace std;
+
+
+
+// This is a binary tree data structure, i.e. a tree data structure
+// which decomposes a cuboid domain into two non-overlapping cuboid
+// subdomains.  Each node splits a domain in exactly one direction.
+// Subdomains cannot be empty.
+
+// All intervals are closed at the lower and open at the upper
+// boundary.  This makes it possible to combine adjacent such
+// intervals, obtaining again an interval with this property.  (In
+// particular, if bboxes are stored, the upper bound of bboxes must be
+// increased by the stride to arrive at such intervals.)
+
+
+// Generic arithmetic search
+template <typename T, int D>
+static
+int asearch (T t, vect <T, D> const & ts);
+
+
+// T: index space (usually integer, or maybe real)
+// D: number of dimensions (rank)
+// P: payload (e.g. region_t)
+template <typename T, int D, typename P>
+class bintree {
+  
+public:
+  
+  // Short name for a small vector
+  typedef vect <T, D> tvect;
+  
+private:
+  
+  enum tree_type_t { type_empty, type_branch, type_leaf } type;
+  
+  // Direction in which the node is split (0 <= dir < D)
+  int dir;
+  
+  // If this is a branch:
+  // 3 bounds, splitting the domain
+  vect <T, 3> bounds;
+  // 2 pointers to subtrees
+  vect <bintree *, 2> subtrees;
+  
+  // If this is a leaf:
+  // just the payload
+  P p;
+  
+public:
+  
+  // Create an empty tree
+  bintree ();
+  
+  // Create a tree branch from a list of bounds and subtrees
+  bintree (int dir_, vect <T, 3> const & bounds_,
+           vect <bintree *, 2> const & subtrees_);
+  
+  // Create a tree leaf from a payload
+  bintree (P const & p_);
+  
+  // Create a tree as copy from another tree
+  bintree (bintree const & t);
+  
+  // Copy a tree
+  bintree & operator= (bintree const & t);
+  
+  // Delete a tree (and its subtrees)
+  ~bintree ();
+  
+  // Inquire tree properties
+  bool empty() const { return type == type_empty; }
+  bool is_branch() const { return type == type_branch; }
+  bool is_leaf() const { return type == type_leaf; }
+  
+  // Compare trees
+  bool operator== (bintree const & t) const;
+  bool operator!= (bintree const & t) const
+  { return not (*this == t); }
+  
+  // Invariant
+  bool invariant () const;
+  
+  // Access the payload
+  P const & payload () const { assert (is_leaf()); return p; }
+  P & payload () { assert (is_leaf()); return p; }
+  
+  // Find the leaf payload corresponding to a position
+  P const * search (tvect const & ipos) const;
+  P * search (tvect const & ipos);
+  
+  class const_iterator {
+    bintree const & f;
+    int i;
+    const_iterator * it;
+  public:
+    const_iterator (bintree const & f_);
+    const_iterator (bintree const & f_, int);
+    ~const_iterator ();
+  private:
+    const_iterator (const_iterator const &);
+    const_iterator & operator= (const_iterator const &);
+  public:
+    bintree const & operator* () const;
+    bool operator== (const_iterator const & it2) const;
+    bool operator!= (const_iterator const & it2) const
+    { return not (*this == it2); }
+    const_iterator & operator++ ();
+    bool done () const;
+  };
+  
+  class iterator {
+    bintree & f;
+    int i;
+    iterator * it;
+  public:
+    iterator (bintree & f_);
+    iterator (bintree & f_, int);
+    ~iterator ();
+  private:
+    iterator (iterator const &);
+    iterator & operator= (iterator const &);
+  public:
+    bintree & operator* () const;
+    bool operator== (iterator const & it2) const;
+    bool operator!= (iterator const & it2) const
+    { return not (*this == it2); }
+    iterator & operator++ ();
+    bool done () const;
+  };
+  
+  // Memory usage
+  size_t memory () const CCTK_ATTRIBUTE_PURE;
+  
+  // Output helper
+  void output (ostream & os) const;
+};
+
+
+
+// Memory usage
+template <typename T, int D, typename P>
+inline size_t memoryof (bintree<T,D,P> const & f) CCTK_ATTRIBUTE_PURE;
+template <typename T, int D, typename P>
+inline size_t memoryof (bintree<T,D,P> const & f) { return f.memory(); }
+
+
+
+// Output
+template <typename T, int D, typename P>
+ostream &
+operator<< (ostream & os, bintree<T,D,P> const & f)
+{
+  f.output (os);
+  return os;
+}
+
+
+
+#endif  // #ifndef BINTREE_HH
diff --git a/Carpet/CarpetLib/src/make.code.defn b/Carpet/CarpetLib/src/make.code.defn
index ccf896ae5..96d621172 100644
--- a/Carpet/CarpetLib/src/make.code.defn
+++ b/Carpet/CarpetLib/src/make.code.defn
@@ -1,8 +1,11 @@
 # Main make.code.defn file for thorn CarpetLib   -*-Makefile-*-
 
 # Source files in this directory
-SRCS =	bbox.cc					\
+SRCS =	balance.cc				\
+	bbox.cc					\
 	bboxset.cc				\
+	bboxtree.cc				\
+	bintree.cc				\
 	commstate.cc				\
 	data.cc					\
 	defs.cc					\
diff --git a/Carpet/CarpetLib/src/region.cc b/Carpet/CarpetLib/src/region.cc
index 00dabc410..eaf04d490 100644
--- a/Carpet/CarpetLib/src/region.cc
+++ b/Carpet/CarpetLib/src/region.cc
@@ -88,6 +88,50 @@ operator== (region_t const & a, region_t const & b)
 
 
 
+// Assign a load to a region
+CCTK_REAL
+region_t::load () const
+{
+  return extent.size();
+}
+
+
+
+// Split a region into two
+region_t
+region_t::split (CCTK_REAL const ratio_new_over_old)
+{
+  assert (ratio_new_over_old >= 0 and ratio_new_over_old <= 1);
+  if (extent.empty()) {
+    // Don't do anything for empty regions
+    return *this;
+  }
+  // Choose a direction (prefer the z direction)
+  int const idir = maxloc1 (extent.shape());
+  int const np = extent.shape()[idir];
+  // Keep the lower part, and split off the upper part
+  int const new_np = floor (np * ratio_new_over_old + 0.5);
+  int const keep_np = np - new_np;
+  // Calculate new region extents
+  ivect const lo = extent.lower();
+  ivect const up = extent.upper();
+  ivect const str = extent.stride();
+  ivect const cut = lo + str * ivect::dir(idir) * keep_np;
+  
+  // Copy the region
+  region_t newreg = *this;
+  // Set new extents
+  extent = ibbox (lo, cut-str, str);
+  newreg.extent = ibbox (cut, up, str);
+  // Mark cutting boundary as not outer boundary
+  outer_boundaries[idir][1] = false;
+  newreg.outer_boundaries[idir][0] = false;
+  
+  return newreg;
+}
+
+
+
 // Combine a collection of regions.  Regions can be combined if they
 // abutt on boundaries which are not outer boundaries, ignoring the
 // processor distribution.  This should lead to a canonical
diff --git a/Carpet/CarpetLib/src/region.hh b/Carpet/CarpetLib/src/region.hh
index ed14b024d..f47c38617 100644
--- a/Carpet/CarpetLib/src/region.hh
+++ b/Carpet/CarpetLib/src/region.hh
@@ -26,6 +26,10 @@ struct region_t {
   ~region_t ();
   
   bool invariant () const CCTK_ATTRIBUTE_PURE;
+  
+  // For regridding
+  CCTK_REAL load () const;
+  region_t split (CCTK_REAL ratio_new_over_old);
 };