librheolef-dev 6.7-6 / usr / include / rheolef / hack_array_mpi.icc

///
/// This file is part of Rheolef.
///
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
///
/// Rheolef is free software; you can redistribute it and/or modify
/// it under the terms of the GNU General Public License as published by
/// the Free Software Foundation; either version 2 of the License, or
/// (at your option) any later version.
///
/// Rheolef is distributed in the hope that it will be useful,
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
/// GNU General Public License for more details.
///
/// You should have received a copy of the GNU General Public License
/// along with Rheolef; if not, write to the Free Software
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
/// 
/// =========================================================================
#include "rheolef/config.h"
#ifdef _RHEOLEF_HAVE_MPI

#include "rheolef/hack_array.h"

namespace rheolef {

// ===============================================================
// allocator
// ===============================================================
template<class T, class A>
hack_array_mpi_rep<T,A>::hack_array_mpi_rep (const A& alloc)
  : base(alloc),
    _stash(),
    _send(),
    _receive(),
    _receive_max_size(),
    _ext_x()
{
}
template<class T, class A>
hack_array_mpi_rep<T,A>::hack_array_mpi_rep (const distributor& ownership, const parameter_type& param, const A& alloc)
  : base(ownership,param,alloc),
    _stash(),
    _send(),
    _receive(),
    _receive_max_size(),
    _ext_x()
{
}
template <class T, class A>
void
hack_array_mpi_rep<T,A>::resize (const distributor& ownership, const parameter_type& param)
{
    base::resize (ownership, param);
    _stash.clear();
    _send.waits.clear();
    _send.data.clear();
    _receive.waits.clear();
    _receive.data.clear();
    _receive_max_size = 0;
}
// ===============================================================
// assembly
// ===============================================================
template <class T, class A>
void
hack_array_mpi_rep<T,A>::set_dis_entry (size_type dis_i, const generic_value_type& val)
{
  if (ownership().is_owned (dis_i)) {
    size_type i = dis_i - ownership().first_index();
    operator[](i) = val;
    return;
  }
  assert_macro (dis_i < ownership().dis_size(), "index "<<dis_i
	<< " is out of range [0:" << ownership().dis_size() << "[");
  // loop on the raw data vector (fixedsize, run-time known)
  size_type dis_iraw = base::_value_size*dis_i + 1;
  typename generic_value_type::const_iterator iter = val._data_begin();
  for (size_type loc_iraw = 0, loc_nraw = val._data_size(); loc_iraw < loc_nraw; loc_iraw++, iter++) {
    _stash.insert (std::pair<const size_type,raw_type>(dis_iraw + loc_iraw, *iter));
  }
}
template <class T, class A>
void
hack_array_mpi_rep<T,A>::dis_entry_assembly_begin ()
{
    _receive_max_size = mpi_assembly_begin (
        _stash,
        make_apply_iterator(_stash.begin(), first_op<typename stash_map_type::value_type>()),
        make_apply_iterator(_stash.end(),   first_op<typename stash_map_type::value_type>()),
        raw_base::ownership(),
        _receive,
        _send);

    _stash.clear();
}
template <class T, class A>
void
hack_array_mpi_rep<T,A>::dis_entry_assembly_end()
{
    mpi_assembly_end (
        _receive,
        _send,
        _receive_max_size,
        disarray_make_store(
            raw_base::begin() - raw_base::ownership().first_index(),
            set_op<raw_type,raw_type>(),
            size_type(0),
            std::false_type()));

    _send.waits.clear();
    _send.data.clear();
    _receive.waits.clear();
    _receive.data.clear();
    _receive_max_size = 0;
}
// ===============================================================
// scatter
// ===============================================================
/** @brief get values from ext_idx_set, that are managed by another proc
 * new version: instead of sending automatic_data_type, send data_size*raw_value_type to boost::mpi
 * => should work with mpi::boost and use simple MPI_Datatype instead of mallocated complex one
 */
template <class T, class A>
template <class Set, class Map>
void
hack_array_mpi_rep<T,A>::append_dis_entry (const Set& ext_idx_set, Map& ext_idx_map) const
{
    // 0) declare the local context for raw data
    scatter_message<std::vector<raw_type> > raw_from;
    scatter_message<std::vector<raw_type> > raw_to;

    // 1) convert set to vector, for direct acess & multiply by data_size
    std::vector<size_type> raw_ext_idx (base::_data_size*ext_idx_set.size());
    typename Set::const_iterator iter = ext_idx_set.begin();
    for (size_type i = 0; i < ext_idx_set.size(); i++, iter++) {
      size_type idx = *iter;
      for  (size_type l = 0; l < base::_data_size; l++) {
        size_type raw_i   = base::_data_size*i + l;
        size_type raw_idx = base::_value_size*idx + l + (base::_value_size - base::_data_size);
        raw_ext_idx [raw_i] = raw_idx;
      }
    }
    // 2) declare id[i]=i for scatter
    std::vector<size_type> raw_id (raw_ext_idx.size());
    for (size_type i = 0; i < raw_id.size(); i++) raw_id[i] = i;

    // 3) init scatter
    size_type raw_dis_size = base::_value_size*ownership().dis_size();
    size_type raw_size     = base::_value_size*ownership().size();
    distributor raw_ownership (raw_dis_size, ownership().comm(), raw_size);
    distributor::tag_type tag_init = distributor::get_new_tag();
    mpi_scatter_init(
        raw_ext_idx.size(),
        raw_ext_idx.begin().operator->(),
        raw_id.size(),
        raw_id.begin().operator->(),
        raw_ownership.dis_size(),
        raw_ownership.begin().operator->(),
        tag_init,
        raw_ownership.comm(),
        raw_from,
        raw_to);

    // 4) begin scatter: send local data to others and get ask for missing data
    std::vector<raw_type> raw_buffer (raw_ext_idx.size());
    distributor::tag_type tag = distributor::get_new_tag();
    // access to an itrator to raw_data (behaves as a pointer on raw_type)
    typedef typename base::base raw_base;
    typename raw_base::const_iterator raw_begin = raw_base::begin();
    mpi_scatter_begin (
        raw_begin,
        raw_buffer.begin().operator->(),
        raw_from,
        raw_to,
        set_op<raw_type,raw_type>(),
        tag,
        raw_ownership.comm());

    // 5) end scatter: receive missing data
    mpi_scatter_end (
        raw_begin,
        raw_buffer.begin().operator->(),
        raw_from,
        raw_to,
        set_op<raw_type,raw_type>(),
        tag,
        raw_ownership.comm());

    // 6) unpack raw data: build the associative container: pair (ext_idx ; data)
    iter = ext_idx_set.begin();
    for (size_type i = 0; i < ext_idx_set.size(); i++, iter++) {
      size_type idx = *iter;
      automatic_value_type value (base::_parameter);
      typename automatic_value_type::iterator p = value._data_begin();
      for  (size_type l = 0; l < base::_data_size; l++, p++) {
        size_type raw_i   = base::_data_size*i + l;
        *p = raw_buffer[raw_i];
      }
      ext_idx_map.insert (std::make_pair (idx, value));
    }
}
template <class T, class A>
typename hack_array_mpi_rep<T,A>::const_reference
hack_array_mpi_rep<T,A>::dis_at (const size_type dis_i) const
{
    if (ownership().is_owned(dis_i)) {
      size_type i = dis_i - ownership().first_index();
      return operator[](i);
    }
    typename scatter_map_type::const_iterator iter = _ext_x.find (dis_i);
    check_macro (iter != _ext_x.end(), "unexpected external index="<<dis_i);
    return (*iter).second;
}
template <class T, class A>
void
hack_array_mpi_rep<T,A>::update_dis_entries() const
{  
  std::set<size_type> ext_i;
  for (typename scatter_map_type::const_iterator
	iter = _ext_x.begin(),
	last = _ext_x.end(); iter != last; ++iter) {
    ext_i.insert ((*iter).first);
  }
  get_dis_entry (ext_i, _ext_x);
}
// ===============================================================
// put & get
// ===============================================================
template <class T, class A>
template <class PutFunction>
odiststream&
hack_array_mpi_rep<T,A>::put_values (odiststream& ops, PutFunction put_element) const
{
  distributor::tag_type tag = distributor::get_new_tag();
  std::ostream& os = ops.os();
  
  // determine maximum message to arrive
  size_type max_size;
  mpi::reduce(comm(), size(), max_size, mpi::maximum<size_type>(), 0);

  size_type io_proc = odiststream::io_proc();  
  if (ownership().process() == io_proc) {
    base::put_values (ops, put_element);
    // then, receive and print messages
    std::vector<raw_type> raw_values (max_size*base::_data_size, std::numeric_limits<raw_type>::max());
    for (size_type iproc = 0; iproc < ownership().n_process(); iproc++) {
      if (iproc == io_proc) continue; 
      size_type loc_sz_i = ownership().size(iproc);
      if (loc_sz_i == 0) continue;
      mpi::status status = comm().recv(iproc, tag, raw_values.begin().operator->(), raw_values.size());
      boost::optional<int> n_data_opt = status.count<raw_type>();
      check_macro (n_data_opt, "receive failed");
      size_type n_data = n_data_opt.get();
      check_macro (n_data == loc_sz_i*base::_data_size, "unexpected receive message size");
      typename T::automatic_type tmp (base::_parameter);
      for (size_type i = 0; i < loc_sz_i; i++) {
	typename T::iterator p = tmp._data_begin();
        for (size_type iloc = 0; iloc < base::_data_size; iloc++, p++) {
	  *p = raw_values [i*base::_data_size + iloc];
        }
	put_element (os, tmp);
	os << std::endl;
      }
    }
    os << std::flush;
  } else {
    if (size() != 0) {
      std::vector<raw_type> raw_values (size()*base::_data_size, std::numeric_limits<raw_type>::max());
      for (size_type i = 0, n = size(); i < n; i++) {
        for (size_type j = 0; j < base::_data_size; j++) {
	  raw_values [i*base::_data_size + j] = raw_base::operator[] (i*base::_value_size + j+(base::_value_size - base::_data_size));
        }
      }
      comm().send(io_proc, tag, raw_values.begin().operator->(), raw_values.size());
    }
  }
  return ops;
}
template <class T, class A>
odiststream&
hack_array_mpi_rep<T,A>::put_values (odiststream& ops) const
{
  return put_values (ops, _disarray_put_element_type<generic_value_type>());
}
template <class T, class A>
template <class GetFunction>
idiststream&
hack_array_mpi_rep<T,A>::get_values (idiststream& ips, GetFunction get_element)
{
  distributor::tag_type tag = distributor::get_new_tag();
  std::istream& is = ips.is();
  size_type my_proc = ownership().process();
  size_type io_proc = odiststream::io_proc();
  size_type size_max = 1;
  for (size_type iproc = 0; iproc < ownership().n_process(); iproc++) {
    size_max = std::max (size_max, ownership().size(iproc));
  }
  distributor io_ownership (size_max, comm(), (my_proc == io_proc) ? size_max : 0);
  hack_array_seq_rep<T,A> data_proc_j (io_ownership, base::_parameter);
  if (my_proc == io_proc) {
    // load first chunk associated to proc 0
    check_macro (load_chunk (is, begin(), end(), get_element), "read failed on input stream.");
    if (ownership().n_process() > 1) {
      // read in other chuncks and send to other processors
      // determine maximum chunck owned by other
      std::vector<raw_type> raw_values (size_max*base::_data_size, std::numeric_limits<raw_type>::max());
      typename hack_array_seq_rep<T,A>::iterator start_j = data_proc_j.begin();
      // bizarre qu'on lise ts les blocs dans la meme zone de memoire
      // et qu'on attende pas que ce soit envoye pour ecraser par le suivant ?
      for (size_type jproc = 0; jproc < ownership().n_process(); jproc++) {
        if (jproc == io_proc) continue; 
        // load first chunk associated to proc j
        size_type loc_sz_j = ownership().size(jproc);
	if (loc_sz_j == 0) continue;
        typename hack_array_seq_rep<T,A>::iterator last_j = start_j + loc_sz_j;
        check_macro (load_chunk (is, start_j, last_j, get_element), "read failed on input stream.");
        for (size_type i = 0, n = loc_sz_j; i < n; i++) {
          for (size_type j = 0; j < base::_data_size; j++) {
	    raw_values [i*base::_data_size + j] = data_proc_j.raw_base::operator[] (i*base::_value_size + j+(base::_value_size - base::_data_size));
          }
        }
  	comm().send (jproc, tag, raw_values.begin().operator->(), loc_sz_j*base::_data_size);
      }
    }
  } else {
    if (size() != 0) {
      std::vector<raw_type> raw_values (size()*base::_data_size, std::numeric_limits<raw_type>::max());
      comm().recv (io_proc, tag, raw_values.begin().operator->(), size()*base::_data_size);
      for (size_type i = 0, n = size(); i < n; i++) {
        for (size_type j = 0; j < base::_data_size; j++) {
          raw_base::operator[] (i*base::_value_size + j+(base::_value_size - base::_data_size))
	    = raw_values [i*base::_data_size + j];
        }
      }
    }
  }
  return ips;
}
template <class T, class A>
idiststream&
hack_array_mpi_rep<T,A>::get_values (idiststream& ips)
{
  return get_values (ips, _disarray_get_element_type<generic_value_type>());
}
template <class T, class A>
template <class PutFunction, class Permutation>
odiststream&
hack_array_mpi_rep<T,A>::permuted_put_values (
  odiststream&                       ops,
  const Permutation&                 perm,
  PutFunction                        put_element) const
{
  // NOTE: could be merged with disarray::permuted_put_value : same code exactly
  assert_macro (perm.size() == size(), "permutation size does not match");
  size_type io_proc = odiststream::io_proc();
  size_type my_proc = comm().rank();
  distributor io_ownership (dis_size(), comm(), (my_proc == io_proc) ? dis_size() : 0);
  hack_array_mpi_rep<T,A>  perm_x (io_ownership, base::_parameter);
  for (size_type i = 0, n = size(); i < n; i++) {
    perm_x.dis_entry (perm[i]) = operator[](i);
  }
  perm_x.dis_entry_assembly();
  perm_x.hack_array_seq_rep<T,A>::put_values (ops, put_element);
  return ops;
}
// ===============================================================
// repartition
// ===============================================================
template <class T, class A>
template <class A2>
void
hack_array_mpi_rep<T,A>::repartition (			// old_numbering for *this
        const disarray_rep<size_type,distributed,A2>&         partition,	// old_ownership
        hack_array_mpi_rep<T,A>&                new_array,	// new_ownership
        disarray_rep<size_type,distributed,A2>&    old_numbering,	// new_ownership
        disarray_rep<size_type,distributed,A2>&    new_numbering) const // old_ownership
{
  using namespace std;
  communicator comm = ownership().comm();
  size_type nproc   = comm.size();
  size_type my_proc = comm.rank();
  vector<size_type> send_local_elt_size (nproc, 0);
  typename disarray_rep<size_type,distributed,A2>::const_iterator iter_part = partition.begin();
  for (size_type ie = 0; ie < partition.size(); ie++, iter_part++) {
    send_local_elt_size [*iter_part]++;
  }
  vector<size_type> recv_local_elt_size (nproc, 0);
  all_to_all (comm, send_local_elt_size, recv_local_elt_size);
  vector<size_type> recv_local_elt_start (nproc+1);
  recv_local_elt_start [0] = 0;
  for (size_type iproc = 0; iproc < nproc; iproc++) {
    recv_local_elt_start [iproc+1] = recv_local_elt_start [iproc] + recv_local_elt_size[iproc];
  }
  vector<size_type> send_local_elt_start (nproc);
  all_to_all (comm, recv_local_elt_start.begin().operator->(), send_local_elt_start.begin().operator->());
  size_type new_local_n_elt = recv_local_elt_start [nproc];
  size_type global_n_elt = dis_size();

  // re-distribute data:
  distributor new_elt_ownership (global_n_elt, comm, new_local_n_elt);
  new_array.resize     (new_elt_ownership, base::_parameter); // CHANGE FROM ARRAY here
  old_numbering.resize (new_elt_ownership, numeric_limits<size_type>::max());
  new_numbering.resize (ownership(), numeric_limits<size_type>::max());
  iter_part = partition.begin();
  const_iterator iter_elt = begin();
  typename disarray_rep<size_type,distributed,A2>::iterator iter_new_num_elt = new_numbering.begin();
  for (size_type ie = 0, ne = partition.size(); ie < ne; ie++, iter_part++, iter_elt++, iter_new_num_elt++) {
    size_type iproc      = *iter_part;
    const generic_value_type& x = *iter_elt; // CHANGE FROM ARRAY here
    size_type new_global_ie = new_elt_ownership[iproc] + send_local_elt_start[iproc];
    new_array.dis_entry (new_global_ie) = x;
    *iter_new_num_elt = new_global_ie;
    size_type old_global_ie = ownership()[my_proc] + ie;
    old_numbering.dis_entry (new_global_ie) = old_global_ie;
    send_local_elt_start[iproc]++;
  }
  new_array.dis_entry_assembly();
  old_numbering.template dis_entry_assembly<typename default_set_op<size_type>::type>();
}

} // namespace rheolef
#endif // _RHEOLEF_HAVE_MPI