MPI_ULoadBalance.h

浏览该文件的文档。

#ifndef __MPI_ULoadBalance_h__
#define __MPI_ULoadBalance_h__

#ifdef __MPI_LoadBalance_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_UGeometry_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;
    
    struct rank_map_val {
      int old_rank; 
      bool is_dummy; 
    };
    typedef std::map<int,rank_map_val> rank_map_t;
  private:
  // 在计算每个几何体上的负载是使用的性质
  property_id_t<u_int> _pid_loading;

  property_id_t<rank_map_t> _pid_rank_map;

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

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

  // (进程，列表(全局标号，维数，对象指针)组)的映射
  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; 

  typename forest_t::container_t _nre; 

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

  private:
  template <class GEO> rank_map_t *
  new_rank_map(const 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_merge_later.clear();
    _global_pointer_to_share.clear();
    _cut_point.clear();
    _new_rank.clear();
  }

  template <class GEO> rank_map_t *
  get_rank_map(const 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);
  }

  template <class GEO> bool 
    is_on_this_new_rank(const GEO& geo,
                        int new_rank) const {
    rank_map_t * p_map = this->get_rank_map(geo);
    if (p_map == NULL) return false;
    typename rank_map_t::iterator
      the_pair = p_map->find(new_rank);
    return (the_pair != p_map->end());
  }

  template <class GEO> bool
    is_save_on_this_rank(const GEO& geo,
                         int new_rank) const {
    rank_map_t * p_map = this->get_rank_map(geo);
    if (p_map == NULL) return false;
    typename rank_map_t::iterator
      the_pair = p_map->find(new_rank);
    if (the_pair == p_map->end()) return false;

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

  template <class GEO> bool 
    is_dummy_on_this_new_rank(const GEO& geo,
                                int new_rank) const {
    rank_map_t * p_map = this->get_rank_map(geo);
    if (p_map == NULL) return false;
    typename rank_map_t::iterator
      the_pair = p_map->find(new_rank);
    if (the_pair == p_map->end()) return false;
    return (the_pair->second.is_dummy);
  }

  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));

        std::map<unsigned long, std::list<void *> >::iterator
          the_entry = _global_pointer_to_merge_later.find(global_idx);
        if (the_entry != _global_pointer_to_merge_later.end()) {
          std::list<void *>::iterator
            the_ptr = the_entry->second.begin(),
            end_ptr = the_entry->second.end();
          for (;the_ptr != end_ptr;++ the_ptr) {
            GEO ** geo_ptr_ptr = (GEO **)(*the_ptr);
            (*geo_ptr_ptr) = p_geo;
          }
          _global_pointer_to_merge_later.erase(the_entry);
        }
      } else {
        assert (the_pair->second == p_geo);
      }
    } else {
      assert (type == 4);
      if (the_pair != _global_pointer_to_merge.end()) {
        p_geo = (GEO *)(the_pair->second);
      } else {
        _global_pointer_to_merge_later[global_idx].push_back((void *)(&p_geo));
      }
    }
  }

  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, (void *)(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) {
        const unsigned long& global_idx = the_tuple->first;
        int& dimension = the_tuple->second.first;
        void *& p_obj = the_tuple->second.second;
        os << global_idx << dimension << p_obj;
      }
    }

    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));

        if (dimension == 0) {
#define GDIM 0
#define GEO HGeometry<GDIM,dow>
#define OBJ Shared_object<GEO>
          GEO * p_geo = (GEO *)(the_pair.second);
          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 *)(the_pair.second);
          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 *)(the_pair.second);
          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 *)(the_pair.second);
          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
        }
      }
    }

#if 1 

    the_pair = _global_pointer_to_share.begin();
    for (;the_pair != end_pair;++ the_pair) {
      typename tuple_map_t::iterator
        the_tuple = the_pair->second.begin(),
        end_tuple = the_pair->second.end();
      for (;the_tuple != end_tuple;++ the_tuple) {
        void *& p_obj = the_tuple->second.second;
        if (dimension == 0) {
#define GDIM 0
#define GEO HGeometry<GDIM,dow>
          GEO * p_geo = (GEO *)p_obj;
          assert (_forest->get_shared_info(*p_geo) != NULL);
#undef GEO
#undef GDIM
        } else if (dimension == 1) {
#define GDIM 1
#define GEO HGeometry<GDIM,dow>
          GEO * p_geo = (GEO *)p_obj;
          assert (_forest->get_shared_info(*p_geo) != NULL);
#undef GEO
#undef GDIM
        } else if (dimension == 2) {
#define GDIM 2
#define GEO HGeometry<GDIM,dow>
          GEO * p_geo = (GEO *)p_obj;
          assert (_forest->get_shared_info(*p_geo) != NULL);
#undef GEO
#undef GDIM
        } else if (dimension == 3) {
#define GDIM 3
#define GEO HGeometry<GDIM,dow>
          GEO * p_geo = (GEO *)p_obj;
          assert (_forest->get_shared_info(*p_geo) != NULL);
#undef GEO
#undef GDIM
        }
      }
    }
#endif
  }

  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,
                 double tol_percent = 0.05,
                 bool is_renumerate_root = false,
                 void (*f)(double,double,double,
                           double&,double&,double&) = NULL) {
    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;

    this->refine_root_element(mean_loading*tol_percent, is_renumerate_root, f);

    _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); 
  }

  private:
  u_int lump_loading() {
    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);
    }

    MPI_Comm comm = _forest->communicator();
    sync_data(comm, _forest->template get_shared_list<dim>(), 
              *this, &this_t::pack_loading, &this_t::unpack_loading);
    return loading;
  }

  public:
  void pack_loading(HGeometry<dim,dow> * geo,
                    int remote_rank,
                    AFEPack::ostream<>& os) {
    u_int * p_loading = geo->get_property(_pid_loading);
    if (p_loading != NULL) {
      os << *p_loading;
    } else {
      u_int loading = 0;
      os << loading;
    }
  }
  void unpack_loading(HGeometry<dim,dow> * geo,
                      int remote_rank,
                      AFEPack::istream<>& is) {
    assert ((_forest->get_shared_info(*geo) != NULL));
    u_int remote_loading;
    is >> remote_loading;
    u_int * p_loading = geo->get_property(_pid_loading);
    if (p_loading != NULL) {
      (*p_loading) += remote_loading;
    }
  }

  private:

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


  void refine_root_element(double tol, bool is_renum_root, 
                           void (*f)(double,double,double,
                                     double&,double&,double&)) {
    _nre.clear();

    typename forest_t::RootIterator
    the_ele = _forest->beginRootElement(),
    end_ele = _forest->endRootElement();
    for (;the_ele != end_ele;++ the_ele) {
      this->refine_root_element(tol, *the_ele, _nre);
    }
    _forest->rootElement().swap(_nre);
    if (is_renum_root) {
      _forest->renumerateRootElement(f);
    }
  }

  void refine_root_element(double tol,
                           HGeometry<dim,dow>& geo,
                           typename forest_t::container_t& nre) const {
    u_int * p_loading = geo.get_property(_pid_loading);
    if (*p_loading > tol) {
      for (u_int i = 0;i < geo.n_child;++ i) {
        this->refine_root_element(tol, *geo.child[i], nre);
      }
    } else {
      nre.push_back(&geo);
    }
  }
            
  void coarse_root_element() {
    property_id_t<bool> pid_is_found;
    new_property_id(pid_is_found);

    _nre.clear();
    typename forest_t::RootIterator
    the_ele = _forest->beginRootElement(),
    end_ele = _forest->endRootElement();
    for (;the_ele != end_ele;++ the_ele) {
      this->coarse_root_element(pid_is_found, *the_ele, _nre);
    }
    _forest->rootElement().swap(_nre);
  }

  void coarse_root_element(const property_id_t<bool>& pid,
                           HGeometry<dim,dow>& geo,
                           typename forest_t::container_t& nre) const {
    if (geo.get_property(pid) != NULL) return;
    if ((geo.parent != NULL) &&
        (_forest->get_shared_info(*geo.parent) == NULL)) {
      this->coarse_root_element(pid, *geo.parent, nre);
    } else {
      nre.push_back(&geo);
      geo.new_property(pid);
    }
  }


  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]);
      }
    }
    this->set_new_rank_up();

    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) {
    rank_map_t * p_map = get_rank_map(*geo);
    if (p_map == NULL) p_map = new_rank_map(*geo);

    os << p_map->size();
    typename rank_map_t::iterator
    the_pair = p_map->begin(),
    end_pair = p_map->end();
    for (;the_pair != end_pair;++ the_pair) {
      const int& new_rank = the_pair->first;
      const int& old_rank = the_pair->second.old_rank;
      the_pair->second.is_dummy = _forest->is_dummy(*geo);
      const bool& is_dummy = the_pair->second.is_dummy;

      os << new_rank << old_rank << is_dummy;
    }
  }

  template <int GDIM> void
  unpack_set_new_rank(HGeometry<GDIM,dow> * geo,
                      int remote_rank,
                      AFEPack::istream<>& is) {
    assert ((_forest->get_shared_info(*geo) != NULL));

    rank_map_t * p_map = get_rank_map(*geo);
    if (p_map == NULL) p_map = new_rank_map(*geo);

    int new_rank;
    std::size_t n;
    is >> n;
    for (u_int i = 0;i < n;++ i) {
      rank_map_val remote_rank_map;
      is >> new_rank 
         >> remote_rank_map.old_rank
         >> remote_rank_map.is_dummy;
      typename rank_map_t::iterator 
        the_pair = p_map->find(new_rank);
      if (the_pair == p_map->end()) {
        (*p_map)[new_rank] = remote_rank_map;
      } else {
        if (the_pair->second.is_dummy == remote_rank_map.is_dummy) {
          if (the_pair->second.old_rank > remote_rank_map.old_rank) {
            the_pair->second.old_rank = remote_rank_map.old_rank;
          }
        } else if (! remote_rank_map.is_dummy) {
          the_pair->second = remote_rank_map;
        }
      }
    }
  }

  private:
  template <class GEO> void
    geometry_set_new_rank(const GEO& geo, int new_rank) const {
    rank_map_t * p_map = get_rank_map(geo);
    if (p_map == NULL) {
      p_map = new_rank_map(geo);
    } else if (p_map->find(new_rank) != p_map->end()) {
      return; 
    }

    (*p_map)[new_rank].old_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);
      }
    }
  }

  template <class GEO> void
    geometry_set_new_rank_up(const GEO& geo, int new_rank) const {
    rank_map_t * p_map = get_rank_map(geo);
    if (p_map == NULL) p_map = new_rank_map(geo);

    (*p_map)[new_rank].old_rank = _forest->rank(); 
    if (GEO::dim == 0) return;
    
    for (u_int i = 0;i < geo.n_vertex;++ i) {
      geometry_set_new_rank_up(*geo.vertex[i], new_rank);
    }
    for (u_int i = 0;i < geo.n_boundary;++ i) {
      geometry_set_new_rank_up(*geo.boundary[i], new_rank);
    }
    GEO* const& p_parent = geo.parent;
    if (p_parent != NULL) {
      if (! is_on_this_new_rank(*p_parent, new_rank)) {
        geometry_set_new_rank_up(*p_parent, new_rank);
      }

      for (u_int i = 0;i < p_parent->n_child;++ i) {
        GEO* const& p_sib = p_parent->child[i];
        if (p_sib == &geo) continue; 
        if (this->is_on_this_new_rank(*p_sib, new_rank)) continue;
        geometry_set_new_rank_up(*p_sib, new_rank);
      }
    }
  }

  void set_new_rank_up() {
    typename forest_t::RootIterator
    the_ele = _forest->beginRootElement(),
    end_ele = _forest->endRootElement();
    for (;the_ele != end_ele;++ the_ele) {
      rank_map_t * p_map = get_rank_map(*the_ele);
      assert (p_map != NULL);
      typename rank_map_t::iterator
        the_pair = p_map->begin(),
        end_pair = p_map->end();
      for (;the_pair != end_pair;++ the_pair) {
        geometry_set_new_rank_up(*the_ele, the_pair->first);
      }
    }
  }

  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) {
      HGeometry<dim,dow> * p_geo = &*the_ele;
      do {
        if (p_geo->parent != NULL) {
          p_geo = p_geo->parent;
        } else break;
      } while (true);
      geometry_global_index(*p_geo, 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) {
      HGeometry<dim,dow> * p_geo = &*the_ele;
      do {
        if (p_geo->parent != NULL) {
          p_geo = p_geo->parent;
        } else break;
      } while (true);
      geometry_global_index(*p_geo, 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) {
    assert ((_forest->get_shared_info(*geo) != NULL));

    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; 

    rank_map_t * p_map = get_rank_map(geo);
    if (_forest->get_shared_info(geo) != NULL || 
        (p_map != 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);
    _forest->rootElement().swap(_nre); 
    _nre.clear();
    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) {
    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) {
    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) {
      const HGeometry<dim,dow> * p_geo = &(*the_ele);
      do {
        if (p_geo->parent != NULL) {
          p_geo = p_geo->parent;
        } else break;
      } while (true);
      set_shared_info_sent(*p_geo);
    }
  }

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

    for (u_int i = 0;i < geo.n_boundary;++ i) {
      this->set_dummy_flag(*geo.boundary[i]);
    }
    if (geo.isRefined()) {
      for (u_int i = 0;i < geo.n_child;++ i) {
        this->set_dummy_flag(*geo.child[i]);
      }
    }
  }

  template <class GEO>
    void clear_dummy_flag(GEO& geo) const {
    if (_forest->is_dummy(geo)) {
      _forest->undummy(geo);
    }

    for (u_int i = 0;i < geo.n_boundary;++ i) {
      this->clear_dummy_flag(*geo.boundary[i]);
    }
    if (geo.isRefined()) {
      for (u_int i = 0;i < geo.n_child;++ i) {
        this->clear_dummy_flag(*geo.child[i]);
      }
    }
  }

  void set_dummy_flag() {
    typename forest_t::RootIterator
      the_ele = _forest->beginRootElement(),
      end_ele = _forest->endRootElement();
    for (;the_ele != end_ele;++ the_ele) {
      const HGeometry<dim,dow> * p_geo = &(*the_ele);
      do {
        if (p_geo->parent != NULL) {
          p_geo = p_geo->parent;
        } else break;
      } while (true);
      set_dummy_flag(*p_geo);
    }

    the_ele = _forest->beginRootElement();
    for (;the_ele != end_ele;++ the_ele) {
      clear_dummy_flag(*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) {
          if (! _global_pointer_to_merge_later.empty()) {
            std::cerr << "Rank " << _forest->rank() << ": ";
            std::map<unsigned long, std::list<void *> >::iterator
              the_entry = _global_pointer_to_merge_later.begin(),
              end_entry = _global_pointer_to_merge_later.end();
            for (;the_entry != end_entry;++ the_entry) {
              std::cerr << "(" << the_entry->first 
                        << "->" << the_entry->second.size()
                << ")";
            }
            std::cerr << std::endl;
            assert (false);
          }

          this->share_global_pointer();
          this->coarse_root_element();
          this->set_shared_info_sent();
          this->set_dummy_flag();

          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_ULoadBalance_h__