MPI_LoadBalance.h

浏览该文件的文档。

#ifndef __MPI_LoadBalance_h__
#define __MPI_LoadBalance_h__

#ifdef __MPI_ULoadBalance_h__
#error "MPI_LoadBalance.h 和 MPI_ULoadBalance.h 是冲突的，不能一起用。"
#endif

#include <sys/types.h>
#include <sys/stat.h>
#include <unistd.h>
#include <cstdio>

#include <boost/program_options.hpp>
#include <boost/archive/binary_iarchive.hpp>
#include <boost/archive/binary_oarchive.hpp>
#include <boost/archive/text_iarchive.hpp>
#include <boost/archive/text_oarchive.hpp>

#include <AFEPack/Serialization.h>
#include "MPI.h"
#include "MPI_HGeometry_archive.h"
#include "MPI_Migration.h"

namespace MPI {

  template <class FOREST>
    class HLoadBalance {
  public:
    typedef FOREST forest_t;
    enum {dim = FOREST::dim, dow = FOREST::dow};
    typedef typename forest_t::matcher_t matcher_t;
  private:
  // 在计算每个几何体上的负载是使用的性质
  property_id_t<u_int> _pid_loading;

  // 这个性质的 first 为新分区中的秩， second 为对其进行存储的秩                    
  property_id_t<std::map<int, int> > _pid_rank_map;

  // 部分需要拆分和合并的几何体的全局标号
  property_id_t<unsigned long> _pid_global_idx;

  // 元素为(全局标号，对象指针)的对
  std::map<unsigned long, void *> _global_pointer_to_merge;

  // (进程，列表(全局标号，维数，对象指针)组)的映射
  typedef std::map<unsigned long, std::pair<int, void *> > tuple_map_t;
  std::map<int, tuple_map_t> _global_pointer_to_share;

  typedef BirdView<forest_t> birdview_t;
  typedef BirdViewSet<forest_t> birdview_set_t;
  typedef HLoadBalance<forest_t> this_t;

  forest_t * _forest; 
  std::vector<u_int> _cut_point; 
  std::vector<int> _new_rank; 

  public:
  HLoadBalance(forest_t& forest) : _forest(&forest) {}
  ~HLoadBalance() {}

  private:
  template <class GEO> std::map<int,int> *
  new_rank_map(GEO& geo) const {
    return geo.new_property(_pid_rank_map);
  }
  template <class GEO> unsigned long *
  new_global_idx(const GEO& geo) const {
    return geo.new_property(_pid_global_idx);
  }

  public:
  void clear() {
    free_property_id(_pid_loading);
    free_property_id(_pid_rank_map);
    free_property_id(_pid_global_idx);
    _global_pointer_to_merge.clear();
    _global_pointer_to_share.clear();
    _cut_point.clear();
    _new_rank.clear();
  }

  template <class GEO> std::map<int,int> *
  get_rank_map(GEO& geo) const {
    return geo.get_property(_pid_rank_map);
  }
  template <class GEO> unsigned long *
  get_global_idx(const GEO& geo) const {
    return geo.get_property(_pid_global_idx);
  }

  bool is_save_on_this_rank(std::map<int,int> * map, 
                            int new_rank) {
    typename std::map<int,int>::iterator
    the_pair = map->find(new_rank);
    assert (the_pair != map->end());

    return (the_pair->second == _forest->rank()); 
  }

  template <class GEO>
  void merge_global_pointer(int type,
                            unsigned long global_idx,
                            GEO *& p_geo) {
    std::map<unsigned long, void *>::iterator
    the_pair = _global_pointer_to_merge.find(global_idx);
    if (type == 2) {
      if (the_pair == _global_pointer_to_merge.end()) {
        _global_pointer_to_merge.insert(std::pair<unsigned long,void*>(global_idx, p_geo));
      } else {
        assert (the_pair->second == p_geo);
      }
    } else {
      assert (type == 4 && the_pair != _global_pointer_to_merge.end());
      p_geo = (GEO *)(the_pair->second);
    }
  }

  template <class GEO>
  void share_global_pointer(int rank, 
                            unsigned long global_idx,
                            GEO *& p_geo) {
    _global_pointer_to_share[rank][global_idx] = 
    std::pair<int,void *>(p_geo->dimension, p_geo);
  }

  private:
  void share_global_pointer() {
    std::cerr << "Exchanging shared global pointers ..." << std::endl;

    typedef std::map<int, tuple_map_t> map_t;

    int n = 0;
    std::list<int> target;
    std::list<BinaryBuffer<> > out_buf, in_buf;
    typename map_t::iterator
    the_pair = _global_pointer_to_share.begin(),
    end_pair = _global_pointer_to_share.end();
    for (;the_pair != end_pair;++ the_pair, ++ n) {
      target.push_back(the_pair->first);

      out_buf.push_back(BinaryBuffer<>());
      AFEPack::ostream<> os(out_buf.back());
      
      in_buf.push_back(BinaryBuffer<>());

      u_int n_item = the_pair->second.size();
      os << n_item;
      std::cerr << "sending " << the_pair->second.size()
                << " items from " << _forest->rank()
                << " to " << the_pair->first << std::endl;
      typename tuple_map_t::iterator
        the_tuple = the_pair->second.begin(),
        end_tuple = the_pair->second.end();
      for (;the_tuple != end_tuple;++ the_tuple) {
        os << the_tuple->first 
           << the_tuple->second.first 
           << the_tuple->second.second; 
      }
    }

    MPI_Barrier(_forest->communicator());
    sendrecv_data(_forest->communicator(), n, out_buf.begin(), 
                  in_buf.begin(), target.begin());

    typename std::list<int>::iterator 
    the_rank = target.begin(),
    end_rank = target.end();
    typename std::list<BinaryBuffer<> >::iterator 
    the_buf = in_buf.begin();
    for (;the_rank != end_rank;++ the_rank, ++ the_buf) {
      AFEPack::istream<> is(*the_buf);
      u_int n_item;
      is >> n_item;
      std::cerr << "received " << n_item
                << " items from " << *the_rank
                << " to " << _forest->rank() << std::endl;
      for (u_int i = 0;i < n_item;++ i) {
        unsigned long global_idx;
        int dimension;
        void * remote_obj;
        is >> global_idx >> dimension >> remote_obj;
        std::pair<int,void *>& the_pair = 
          _global_pointer_to_share[*the_rank][global_idx];
        assert (dimension == the_pair.first && 
                dimension >= 0 && dimension <= 3);

        void * local_obj = the_pair.second;
        if (dimension == 0) {
#define GDIM 0
#define GEO HGeometry<GDIM,dow>
#define OBJ Shared_object<GEO>
          GEO * p_geo = (GEO *)local_obj;
          OBJ * p_info = _forest->get_shared_info(*p_geo);
          if (p_info == NULL) p_info = _forest->new_shared_info(*p_geo);
          p_info->add_clone(*the_rank, (GEO *)remote_obj);
#undef OBJ
#undef GEO
#undef GDIM
        } else if (dimension == 1) {
#define GDIM 1
#define GEO HGeometry<GDIM,dow>
#define OBJ Shared_object<GEO>
          GEO * p_geo = (GEO *)local_obj;
          OBJ * p_info = _forest->get_shared_info(*p_geo);
          if (p_info == NULL) p_info = _forest->new_shared_info(*p_geo);
          p_info->add_clone(*the_rank, (GEO *)remote_obj);
#undef OBJ
#undef GEO
#undef GDIM
        } else if (dimension == 2) {
#define GDIM 2
#define GEO HGeometry<GDIM,dow>
#define OBJ Shared_object<GEO>
          GEO * p_geo = (GEO *)local_obj;
          OBJ * p_info = _forest->get_shared_info(*p_geo);
          if (p_info == NULL) p_info = _forest->new_shared_info(*p_geo);
          p_info->add_clone(*the_rank, (GEO *)remote_obj);
#undef OBJ
#undef GEO
#undef GDIM
        } else if (dimension == 3) {
#define GDIM 3
#define GEO HGeometry<GDIM,dow>
#define OBJ Shared_object<GEO>
          GEO * p_geo = (GEO *)local_obj;
          OBJ * p_info = _forest->get_shared_info(*p_geo);
          if (p_info == NULL) p_info = _forest->new_shared_info(*p_geo);
          p_info->add_clone(*the_rank, (GEO *)remote_obj);
#undef OBJ
#undef GEO
#undef GDIM
        }
      }
    }
  }

  public:
  template <class LOADER>
  void config(birdview_t& ir_mesh,
              LOADER& loader,
              u_int (LOADER::*value)(GeometryBM&)) {

    if (&(ir_mesh.getForest()) != _forest) {
      std::cerr << "The mesh should be on the same forest of the HLoadBalance." 
                << std::endl;
      abort();
    }

    RegularMesh<dim,dow>& mesh = ir_mesh.regularMesh();
    new_property_id(_pid_loading); 
    u_int n_ele = mesh.n_geometry(dim);
    for (u_int i = 0;i < n_ele;++ i) {
      HGeometry<dim,dow> * p_geo = mesh.template h_geometry<dim>(i);
      u_int * p_loading = p_geo->get_property(_pid_loading);
      if (p_loading == NULL) {
        p_loading = p_geo->new_property(_pid_loading);
      }
      (*p_loading) = (loader.*value)(mesh.geometry(dim,i));
    }
  }

  void config(birdview_t& ir_mesh,
              u_int (*value)(GeometryBM&) = &_default_element_loading) {
    if (&(ir_mesh.getForest()) != _forest) {
      std::cerr << "The mesh should be on the same forest of the HLoadBalance." 
                << std::endl;
      abort();
    }

    RegularMesh<dim,dow>& mesh = ir_mesh.regularMesh();
    new_property_id(_pid_loading); 
    u_int n_ele = mesh.n_geometry(dim);
    for (u_int i = 0;i < n_ele;++ i) {
      const HGeometry<dim,dow> * p_geo = mesh.template h_geometry<dim>(i);
      u_int * p_loading = p_geo->get_property(_pid_loading);
      if (p_loading == NULL) {
        p_loading = p_geo->new_property(_pid_loading);
      }
      (*p_loading) = (*value)(mesh.geometry(dim,i));
    }
  }

  static u_int _default_element_loading(GeometryBM&) { return 1; }

  public:
  void partition(u_int n_new_rank = 0) {
    int rank =  _forest->rank();
    int n_rank = _forest->n_rank();

    u_int loading = this->lump_loading();
    u_int total_loading, partial_loading;
    MPI_Comm comm = _forest->communicator();
    MPI_Scan(&loading, &partial_loading, 1, MPI_UNSIGNED, MPI_SUM, comm);

    if (rank == n_rank - 1) {
      total_loading = partial_loading;
    }
    if (n_new_rank == 0) n_new_rank = n_rank; 

    MPI_Bcast(&total_loading, 1, MPI_UNSIGNED, n_rank - 1, comm);
    u_int mean_loading = total_loading/n_new_rank;

    _cut_point.clear();
    _new_rank.clear();

    u_int idx = 0;
    _cut_point.push_back(idx); 

    loading  = partial_loading - loading;
    u_int current_rank = loading/mean_loading;
    _new_rank.push_back(current_rank);

    typename forest_t::RootIterator
    the_ele = _forest->beginRootElement(),
    end_ele = _forest->endRootElement();
    for (;the_ele != end_ele;++ the_ele) {
      idx += 1;
      loading += this->get_loading(*the_ele);
      u_int the_rank = loading/mean_loading;
      the_rank = std::min(the_rank, n_new_rank - 1); 
      if (the_rank > current_rank) {
        _cut_point.push_back(idx);
        _new_rank.push_back(++ current_rank);
      }
    }
    if (_cut_point.back() != idx) {
      _cut_point.push_back(idx);
    } else {
      _new_rank.pop_back();
    }

    free_property_id(_pid_loading); 
  }

  void reuse_partition() {
    _cut_point.clear();
    _new_rank.clear();

    _cut_point.resize(2, 0);
    _cut_point[1] = _forest->n_rootElement();
    _new_rank.resize(1, _forest->rank());
  }

  private:
  u_int lump_loading() const {
    u_int loading = 0;
    typename forest_t::RootIterator
    the_ele = _forest->beginRootElement(),
    end_ele = _forest->endRootElement();
    for (;the_ele != end_ele;++ the_ele) {
      loading += this->get_loading(*the_ele);
    }
    return loading;
  }

  template <class GEO> u_int
  get_loading(GEO& geo) const {
    u_int * p_loading = geo.get_property(_pid_loading);
    if (p_loading == NULL) {
      p_loading = geo.new_property(_pid_loading);
      for (u_int i = 0;i < geo.n_child;++ i) {
        (*p_loading) += get_loading(*geo.child[i]);
      }
    }
    return (*p_loading);
  }

  private:
  void set_new_rank() {
    new_property_id(_pid_rank_map);

    typename forest_t::RootIterator
    the_ele = _forest->beginRootElement();
    u_int n_part = _new_rank.size();
    for (u_int i = 0;i < n_part;++ i) {
      for (u_int j = _cut_point[i];j < _cut_point[i + 1];++ j, ++ the_ele) {
        geometry_set_new_rank(*the_ele, _new_rank[i]);
      }
    }

    if (_forest->n_rank() <= 1) return;
    Shared_type_filter::only<0> type_filter;
    MPI_Comm comm = _forest->communicator();
#define SYNC_DATA(D) \
    if (dim >= D) { \
      sync_data(comm, _forest->template get_shared_list<D>(), *this, \
                &this_t::template pack_set_new_rank<D>, \
                &this_t::template unpack_set_new_rank<D>, \
                type_filter); \
    }
    SYNC_DATA(0);
    SYNC_DATA(1);
    SYNC_DATA(2);
    SYNC_DATA(3);
#undef SYNC_DATA
  }

  public:
  template <int GDIM> void
  pack_set_new_rank(HGeometry<GDIM,dow> * geo,
                    int remote_rank,
                    AFEPack::ostream<>& os) {
    std::map<int,int> * p_map = get_rank_map(*geo);
    assert (p_map != NULL);

    os << p_map->size();
    typename std::map<int,int>::iterator
    the_rank = p_map->begin(),
    end_rank = p_map->end();
    for (;the_rank != end_rank;++ the_rank) {
      os << the_rank->first << the_rank->second;
    }
  }

  template <int GDIM> void
  unpack_set_new_rank(HGeometry<GDIM,dow> * geo,
                      int remote_rank,
                      AFEPack::istream<>& is) {
    std::map<int,int> * p_map = get_rank_map(*geo);
    assert (p_map != NULL);

    int new_rank, old_rank;
    std::size_t n;
    is >> n;
    for (int i = 0;i < n;++ i) {
      is >> new_rank >> old_rank;
      typename std::map<int,int>::iterator 
        the_pair = p_map->find(new_rank);
      if (the_pair == p_map->end()) {
        p_map->insert(std::pair<int,int>(new_rank, old_rank));
      } else {
        if (the_pair->second > old_rank) {
          the_pair->second = old_rank;
        }
      }
    }
  }

  private:
  template <class GEO> void
  geometry_set_new_rank(GEO& geo, int new_rank) const {
    std::map<int,int> * p_map = get_rank_map(geo);
    if (p_map == NULL) p_map = new_rank_map(geo);
    p_map->insert(std::pair<int,int>(new_rank, _forest->rank())); 
    
    for (u_int i = 0;i < geo.n_vertex;++ i) {
      geometry_set_new_rank(*geo.vertex[i], new_rank);
    }
    for (u_int i = 0;i < geo.n_boundary;++ i) {
      geometry_set_new_rank(*geo.boundary[i], new_rank);
    }
    if (geo.isRefined()) {
      for (u_int i = 0;i < geo.n_child;++ i) {
        geometry_set_new_rank(*geo.child[i], new_rank);
      }
    }
  }

  private:
  void global_index() {
    new_property_id(_pid_global_idx);
    unsigned long idx = 0;
    typename forest_t::RootIterator
    the_ele = _forest->beginRootElement(),
    end_ele = _forest->endRootElement();
    for (;the_ele != end_ele;++ the_ele) {
      geometry_global_index(*the_ele, idx);
    }
    free_property_id(_pid_global_idx); 

    MPI_Comm comm = _forest->communicator();
    unsigned long global_idx = idx;
    MPI_Scan(&idx, &global_idx, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm);
    idx = global_idx - idx;

    new_property_id(_pid_global_idx);
    the_ele = _forest->beginRootElement();
    for (;the_ele != end_ele;++ the_ele) {
      geometry_global_index(*the_ele, idx);
    }
    
    if (_forest->n_rank() <= 1) return;
    Shared_type_filter::only<0> type_filter;
#define SYNC_DATA(D) \
    if (dim >= D) { \
      sync_data(comm, _forest->template get_shared_list<D>(), *this, \
                &this_t::template pack_global_index<D>, \
                &this_t::template unpack_global_index<D>, \
                type_filter); \
    }
    SYNC_DATA(0);
    SYNC_DATA(1);
    SYNC_DATA(2);
    SYNC_DATA(3);
#undef SYNC_DATA
  }

  public:
  template <int GDIM> void
  pack_global_index(HGeometry<GDIM,dow> * geo,
                    int remote_rank,
                    AFEPack::ostream<>& os) {
    unsigned long * p_idx = get_global_idx(*geo);
    assert (p_idx != NULL);

    os << *p_idx; 
  }

  template <int GDIM> void
  unpack_global_index(HGeometry<GDIM,dow> * geo,
                      int remote_rank,
                      AFEPack::istream<>& is) {
    unsigned long * p_idx = get_global_idx(*geo);
    assert (p_idx != NULL);

    unsigned long remote_idx;
    is >> remote_idx;
    if (*p_idx > remote_idx) {
      *p_idx = remote_idx; 
    }
  }

  private:
  template <class GEO> void
  geometry_global_index(GEO& geo, 
                        unsigned long& idx) const {
    unsigned long * p_idx = get_global_idx(geo);
    if (p_idx != NULL) return; 

    std::map<int,int> * p_map = get_rank_map(geo);
    if (_forest->get_shared_info(geo) != NULL || 
        p_map->size() > 1) { 
      p_idx = new_global_idx(geo);
      *p_idx = idx ++;
    }

    for (u_int i = 0;i < geo.n_vertex;++ i) {
      geometry_global_index(*geo.vertex[i], idx);
    }
    for (u_int i = 0;i < geo.n_boundary;++ i) {
      geometry_global_index(*geo.boundary[i], idx);
    }
    if (geo.isRefined()) {
      for (u_int i = 0;i < geo.n_child;++ i) {
        geometry_global_index(*geo.child[i], idx);
      }
    }
  }

  public:
  void write_config_file(const std::string& dirname) {
    char filename[1024];
    sprintf(filename, "%s/.config", dirname.c_str());
    std::ofstream os(filename);
    os << "old-rank=" << _forest->n_rank() << std::endl;
    os.close();
    Migration::save_config(dirname);
  }

  private:
  void export_birdview(birdview_t& ir_mesh, 
                       Migration::data_id_t data_id, 
                       u_int idx) {
    typename birdview_t::ActiveIterator
    the_ele = ir_mesh.beginActiveElement(),
    end_ele = ir_mesh.endActiveElement();
    for (;the_ele != end_ele;++ the_ele) {
      HGeometry<dim,dow> * p_geo = the_ele->h_element;
      AFEPack::ostream<> os(p_geo->buffer[data_id]);
      os << idx;
    }
  }

  public:
  void save_data(const std::string& dirname, 
                 birdview_set_t& bvs) {
    set_new_rank();
    global_index();

    char command[1024], *filename = command;

    Migration::data_id_t id = Migration::name_to_id("__internal_birdviewflag__"); 
    if (! Migration::is_valid(id)) {
      id = Migration::register_data_name("__internal_birdviewflag__");
    }
    typename birdview_set_t::iterator
    the_bv = bvs.begin(), end_bv = bvs.end();
    for (u_int idx = 0;the_bv != end_bv;++ the_bv) {
      export_birdview(*the_bv, id, idx ++);
    }

    typedef boost::archive::HGeometry_oarchive<this_t> archive_t;
    typename forest_t::RootIterator
    the_ele = _forest->beginRootElement();
    u_int n_part = _new_rank.size();
    for (u_int i = 0;i < n_part;++ i) {
      sprintf(command, "mkdir -p %s && mkdir -p %s/%d", 
              dirname.c_str(), dirname.c_str(), _new_rank[i]);
      system(command);
      sprintf(filename, "%s/%d/%d.dat", dirname.c_str(), 
              _new_rank[i], _forest->rank());
      std::ofstream os(filename, std::ios::binary);

      {
        archive_t oa(*this, _new_rank[i], os);

        u_int n_ele = _cut_point[i + 1] - _cut_point[i];
        oa & boost::serialization::make_binary_object(&n_ele, sizeof(u_int));
        std::cerr << "saving data in " << filename 
                  << ", # macro element = " << n_ele << std::endl;

        for (u_int j = _cut_point[i];j < _cut_point[i + 1];++ j) {
          HGeometry<dim,dow> * p_ele = *(the_ele ++);
          oa & p_ele;
        }
      }

      os.close();
    }

    if (_forest->rank() == 0) write_config_file(dirname);
    MPI_Barrier(_forest->communicator()); 
  }

  public:
  int load_config_file(const std::string& dirname) {
    namespace po = boost::program_options;
    char filename[1024];
    sprintf(filename, "%s/.config", dirname.c_str());
    Migration::load_config(dirname);

    po::options_description desc("Hidden options");
    desc.add_options()("old-rank", po::value<int>(), "old number of rank");
    po::variables_map vm;
    std::ifstream is(filename);
    po::store(po::parse_config_file(is, desc), vm);
    po::notify(vm);
    return vm["old-rank"].as<int>();
  }

  private:
  void reconstruct_birdview(HElement<dim,dow>& ele,
                            Migration::data_id_t data_id,
                            u_int idx) const {
    HGeometry<dim,dow> * p_geo = ele.h_element;
    typename Migration::buffer_t::iterator it = p_geo->buffer.find(data_id);
    bool is_active = true;
    ele.value = 0;
    if (it == p_geo->buffer.end()) {
      is_active = false;
    } else {
      BinaryBuffer<>& buf = it->second;
      int n = buf.size()/sizeof(u_int), i;
      u_int idx1;
      AFEPack::istream<> is(buf);
      for (i = 0;i < n;++ i) {
        is >> idx1;
        if (idx1 == idx) break;
      }
      if (i == n) is_active = false;
    }
    if (! is_active) {
      ele.value = 1;
      ele.refine();
      for (int i = 0;i < ele.n_child;++ i) {
        reconstruct_birdview(*ele.child[i], data_id, idx);
      }
    }
  }

  void reconstruct_birdview(birdview_t& ir_mesh,
                            Migration::data_id_t data_id,
                            u_int idx) const {
    typename birdview_t::RootIterator
    the_ele = ir_mesh.beginRootElement(),
    end_ele = ir_mesh.endRootElement();
    for (;the_ele != end_ele;++ the_ele) {
      reconstruct_birdview(*the_ele, data_id, idx);
    }
  }

  template <class GEO>
    void set_shared_info_sent(GEO& geo) const {
    if ((_forest->get_shared_info(geo) != NULL) &
        (! _forest->is_shared_info_sent(geo))) {
      _forest->set_shared_info_sent(geo);
    }

    for (u_int i = 0;i < geo.n_vertex;++ i) {
      if (! _forest->is_shared_info_sent(*geo.vertex[i])) {
        _forest->set_shared_info_sent(*geo.vertex[i]);
      }
    }
    for (u_int i = 0;i < geo.n_boundary;++ i) {
      this->set_shared_info_sent(*geo.boundary[i]);
    }
    if (geo.isRefined()) {
      for (u_int i = 0;i < geo.n_child;++ i) {
        this->set_shared_info_sent(*geo.child[i]);
      }
    }
  }

  void set_shared_info_sent() const {
    typename forest_t::RootIterator
      the_ele = _forest->beginRootElement(),
      end_ele = _forest->endRootElement();
    for (;the_ele != end_ele;++ the_ele) {
      this->set_shared_info_sent(*the_ele);
    }
  }

  public:
  void load_data(const std::string& dirname,
                 birdview_set_t& bvs) {
    clear();

    int n_old_rank = load_config_file(dirname);

    MPI_Barrier(_forest->communicator()); 

    char filename[1024];

    typedef boost::archive::HGeometry_iarchive<this_t> archive_t;

    int rank = 0;
    do {
      std::ifstream is;
      do {
        sprintf(filename, "%s/%d/%d.dat", dirname.c_str(), 
                _forest->rank(), rank ++);
        if (rank > n_old_rank) {
          share_global_pointer();
          set_shared_info_sent();

          Migration::data_id_t id = Migration::name_to_id("__internal_birdviewflag__"); 
          typename birdview_set_t::iterator
            the_bv = bvs.begin(), end_bv = bvs.end();
          for (u_int idx = 0;the_bv != end_bv;++ the_bv) {
            the_bv->reinit(*_forest);
            reconstruct_birdview(*the_bv, id, idx ++);
          }
          return;
        }
        is.open(filename, std::ios::binary); 
        if (is.good()) break; 
      } while (1);

      {
        u_int n_ele;
        std::cerr << "loading data from " << filename;
        archive_t ia(*this, _forest->rank(), is); 
        ia & boost::serialization::make_binary_object(&n_ele, sizeof(u_int));
        std::cerr << ", # macro element = " << n_ele << std::endl;

        for (u_int i = 0;i < n_ele;++ i) {
          HGeometry<dim,dow> * p_ele;
          ia & p_ele;
          _forest->rootElement().push_back(p_ele);
        }
      }

      is.close();
    } while (1);

  }

  };

  template <class FOREST>
    void load_forest(const std::string& dirname,
                     FOREST& forest) {
    HLoadBalance<FOREST> hlb(forest);
    BirdViewSet<FOREST> bvs;
    hlb.load_data(dirname, bvs);
  }

  template <class FOREST>
    void load_mesh(const std::string& dirname,
                   FOREST& forest,
                   BirdView<FOREST>& mesh) {
    HLoadBalance<FOREST> hlb(forest);
    BirdViewSet<FOREST> bvs;
    bvs.add(mesh);
    hlb.load_data(dirname, bvs);
  }

  template <class FOREST>
    void load_mesh(const std::string& dirname,
                   FOREST& forest,
                   u_int n_mesh, ...) {
    HLoadBalance<FOREST> hlb(forest);
    typedef BirdView<FOREST> * mesh_ptr_t;

    BirdViewSet<FOREST> bvs;
    va_list ap;
    va_start(ap, n_mesh);
    for (u_int i = 0;i < n_mesh;++ i) {
      mesh_ptr_t p_mesh = va_arg(ap, mesh_ptr_t);
      bvs.add(*p_mesh);
    }
    va_end(ap);

    hlb.load_data(dirname, bvs);
  }

  template <class FOREST>
    void load_mesh_set(const std::string& dirname,
                       FOREST& forest,
                       BirdViewSet<FOREST>& bvs) {
    HLoadBalance<FOREST> hlb(forest);
    hlb.load_data(dirname, bvs);
  }

  void lb_collect_local_data_dir(MPI_Comm comm,
                                 const std::string& src_dir,
                                 const std::string& dst_dir);
  void lb_sync_local_data_dir(MPI_Comm comm,
                              const std::string& src_dir,
                              const std::string& dst_dir);
}

#endif // __MPI_LoadBalance_h__