librheolef-dev 6.7-6 / usr / include / rheolef / disarray_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/disarray.h"

# include "rheolef/mpi_assembly_begin.h"
# include "rheolef/mpi_assembly_end.h"
# include "rheolef/mpi_scatter_init.h"
# include "rheolef/mpi_scatter_begin.h"
# include "rheolef/mpi_scatter_end.h"
# include "rheolef/load_chunk.h"
# include "rheolef/disarray_store.h"
# include "rheolef/rheostream.h"

namespace rheolef {
// ----------------------------------------------------------------------------
// allocators
// ----------------------------------------------------------------------------
template <class T, class A>
disarray_rep<T,distributed,A>::disarray_rep (const disarray_rep<T,distributed,A>& x)
  : base(x),
    _stash(),
    _send(),
    _receive(),
    _receive_max_size(0),
    _ext_x(x._ext_x)
{
    assert_macro(x._stash.size() == 0 && 
		 x._send.data.size() == 0 && 
		 x._receive.data.size() == 0, 
		 "copy during assembly phase: should not be done");
}
template <class T, class A>
disarray_rep<T,distributed,A>::disarray_rep (
        const distributor& ownership,
        const T&  init_val,
        const A&  alloc)
  : base(ownership, init_val,alloc),
    _stash(),
    _send(),
    _receive(),
    _receive_max_size(0),
    _ext_x()
{
}
template <class T, class A>
void
disarray_rep<T,distributed,A>::resize (
        const distributor& ownership,
        const T&  init_val)
{
    base::resize(ownership, init_val);
    _stash.clear();
    _send.waits.clear();
    _send.data.clear();
    _receive.waits.clear();
    _receive.data.clear();
    _receive_max_size = 0;
}
// ===============================================================
// set dis_entry and assembly
// ===============================================================
#ifdef TO_CLEAN
template <class Map>
inline
typename Map::mapped_type&
_stash_create_new_entry (Map& stash, typename Map::size_type dis_i, std::false_type)
{
    typedef typename Map::size_type   size_type;
    typedef typename Map::mapped_type T;
    std::pair<typename Map::iterator,bool> status
     = stash.insert (std::pair<const size_type,T>(dis_i,T()));
    return (*(status.first)).second;
};
template <class MultiMap>
inline
typename MultiMap::mapped_type&
_stash_create_new_entry (MultiMap& stash, typename MultiMap::size_type dis_i, std::true_type)
{
    typedef typename MultiMap::size_type   size_type;
    typedef typename MultiMap::mapped_type U;
    std::pair<typename MultiMap::iterator,bool> status
     = stash.insert (std::pair<const size_type,U>(dis_i,U()));
    return (*(status.first)).second;
};
template <class T, class A>
typename disarray_rep<T,distributed,A>::reference
disarray_rep<T,distributed,A>::new_dis_entry (size_type dis_i)
{
  size_type      first_dis_i = ownership().first_index();
  size_type       last_dis_i = ownership().last_index();
  if (dis_i >= first_dis_i && dis_i < last_dis_i) {
    return disarray_rep<T,distributed,A>::operator[](dis_i - first_dis_i);
  } else {
    assert_macro (dis_i < ownership().dis_size(), "index "<<dis_i
	<< " is out of range [0:" << ownership().dis_size() << "[");
    return _stash_create_new_entry (_stash, dis_i, is_container());
  }
}
#endif // TO_CLEAN
template <class Map>
inline
void
_stash_set (Map& stash, typename Map::size_type dis_i, const typename Map::mapped_type& val, std::false_type)
{
    typedef typename Map::size_type   size_type;
    typedef typename Map::mapped_type T;
    std::pair<typename Map::iterator,bool> status
     = stash.insert (std::pair<const size_type,T>(dis_i,T()));
    (*(status.first)).second = val;
};
template <class Map>
inline
void
_stash_set_add (Map& stash, typename Map::size_type dis_i, const typename Map::mapped_type& val, std::false_type)
{
    typedef typename Map::size_type   size_type;
    typedef typename Map::mapped_type T;
    std::pair<typename Map::iterator,bool> status 
     = stash.insert (std::pair<const size_type,T>(dis_i,T()));
    (*(status.first)).second += val;
};
template <class MultiMap, class T>
inline
void
_stash_set (MultiMap& stash, typename MultiMap::size_type dis_i, const T& val, std::true_type)
{
    // check if a multi-value for dis_i exists and delete it; then insert the new multi-value
    typedef typename MultiMap::iterator    iterator;
    typedef typename MultiMap::size_type   size_type;
    typedef typename MultiMap::mapped_type U;
    std::pair<iterator, iterator> range_dis_i = stash.equal_range (dis_i);
    stash.erase (range_dis_i.first, range_dis_i.second);
    for (typename T::const_iterator iter = val.begin(), last = val.end(); iter != last; iter++) {
      stash.insert (std::pair<const size_type,U>(dis_i,*iter));
    }
};
// insert a "T value" in multi-map when T has fixed size
template <class MultiMap, class U>
inline
void
_stash_set_add_multi (MultiMap& stash, typename MultiMap::size_type dis_i, const U& val, std::false_type)
{
  typedef typename MultiMap::size_type   size_type;
  typedef typename MultiMap::mapped_type W;
  stash.insert (std::pair<const size_type,W>(dis_i,val));
}
// insert a "T value" in multi-map when T has variable size
template <class MultiMap, class U>
inline
void
_stash_set_add_multi (MultiMap& stash, typename MultiMap::size_type dis_i, const U& val, std::true_type)
{
    typedef typename MultiMap::size_type   size_type;
    typedef typename MultiMap::mapped_type W;
    for (typename U::const_iterator iter = val.begin(), last = val.end(); iter != last; iter++) {
      stash.insert (std::pair<const size_type,W>(dis_i,*iter));
    }
}
// insert a "U value" for simple, or variable-sized one, as index_set or pair_set types
template <class MultiMap, class U>
inline
void
_stash_set_add (MultiMap& stash, typename MultiMap::size_type dis_i, const U& val, std::true_type)
{
  _stash_set_add_multi (stash, dis_i, val, typename is_container<U>::type());
}
template <class T, class A>
void
disarray_rep<T,distributed,A>::set_dis_entry (size_type dis_i, const T& val)
{
    size_type      start = ownership().first_index();
    size_type      last  = ownership().last_index();
    if (dis_i >= start && dis_i < last) {
        disarray_rep<T,distributed,A>::operator[](dis_i - start) = val;
    } else {
        assert_macro (dis_i < ownership().dis_size(), "index "<<dis_i
		<< " is out of range [0:" << ownership().dis_size() << "[");
        _stash_set (_stash, dis_i, val, is_container());
    }
}
template <class T, class A>
template <class U>
void
disarray_rep<T,distributed,A>::set_add_dis_entry (size_type dis_i, const U& val)
{
    size_type      start = ownership().first_index();
    size_type      last  = ownership().last_index();
    if (dis_i >= start && dis_i < last) {
        base::operator[](dis_i - start) += val;
    } else {
        assert_macro (dis_i < ownership().dis_size(),
		"index "<<dis_i<<" is out of range [0:"<<ownership().dis_size());
        _stash_set_add (_stash, dis_i, val, is_container());
    }
}
template <class T, class A>
template <class SetOp>
void
disarray_rep<T,distributed,A>::dis_entry_assembly_begin (SetOp my_set_op)
{
  _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>()),
        ownership(),
        _receive,
        _send);

  _stash.clear();
}
template <class T, class A>
template <class SetOp>
void
disarray_rep<T,distributed,A>::dis_entry_assembly_end(SetOp my_set_op)
{
  mpi_assembly_end (
	_receive,
	_send,
	_receive_max_size,
        disarray_make_store(
	    begin() - ownership().first_index(),
	    my_set_op,
	    size_type(0),
            is_container()));
   
  _send.waits.clear();
  _send.data.clear();
  _receive.waits.clear();
  _receive.data.clear();
  _receive_max_size = 0;
}
// ===============================================================
// repartition
// ===============================================================
template <class T, class A>
template <class A2>
void
disarray_rep<T,distributed,A>::repartition (				// old_numbering for *this
        const disarray_rep<size_type,distributed,A2>& partition,	// old_ownership
        disarray_rep<T,distributed,A>&                new_disarray,	// 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_type 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_disarray.resize     (new_elt_ownership);
  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 T& x = *iter_elt;
    size_type new_global_ie = new_elt_ownership[iproc] + send_local_elt_start[iproc];
    new_disarray.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_disarray.template dis_entry_assembly<typename default_set_op<T>::type>();
  old_numbering.template dis_entry_assembly<typename default_set_op<size_type>::type>();
}
template <class T, class A>
template <class A2>
void
disarray_rep<T,distributed,A>::reverse_permutation (	         // old_ownership for *this=iold2dis_inew
        disarray_rep<size_type,distributed,A2>& inew2dis_iold) const // new_ownership
{
  check_macro (inew2dis_iold.dis_size() == dis_size(), "reverse permutation[0:"<<inew2dis_iold.dis_size()
	<<"[ has incompatible dis_range with oriinal permutation[0:"<<dis_size()<<"[");
  size_type first_dis_iold = ownership().first_index();
  for (size_type iold = 0; iold < size(); iold++) {
    size_type dis_iold = first_dis_iold + iold;
    size_type dis_inew = base::operator[] (iold);
    inew2dis_iold.dis_entry (dis_inew) = dis_iold;
  }
  inew2dis_iold.template dis_entry_assembly<typename default_set_op<T>::type>();
}
template <class T, class A>
template <class A2>
void
disarray_rep<T,distributed,A>::permutation_apply (			  // old_numbering for *this
        const disarray_rep<size_type,distributed,A2>& new_numbering,  // old_ownership
        disarray_rep<T,distributed,A>&                new_disarray) const // new_ownership
{
  check_macro (size() == new_numbering.size(),
	"permutation_apply: incompatible disarray("<<size()<<") and permutation("<<new_numbering.size()<<") sizes");
  check_macro (dis_size() == new_disarray.dis_size(),
	"permutation_apply: incompatible disarray("<<dis_size()<<") and permutation("<<new_disarray.dis_size()<<") dis_sizes");
  typename disarray_rep<size_type,distributed,A2>::const_iterator iter_dis_new_ie = new_numbering.begin();
  for (const_iterator iter = begin(), last = end(); iter != last; iter++, iter_dis_new_ie++) {
    size_type dis_new_ie = *iter_dis_new_ie;
    new_disarray.dis_entry (dis_new_ie) = *iter;
  }
  new_disarray.template dis_entry_assembly<typename default_set_op<T>::type>();
}
/// @brief get values from ext_idx_set, that are managed by another proc; here T=simple type
template <class T, class A>
template <class Set, class Map>
void
disarray_rep<T,distributed,A>::append_dis_entry (const Set& ext_idx_set, Map& ext_idx_map, std::false_type) const
{
    // 0) declare the local context
    scatter_message<std::vector<T,A> > from;
    scatter_message<std::vector<T,A> > to;

    // 1) convert set to vector, for direct acess:
    std::vector<size_type> ext_idx (ext_idx_set.size());
    std::copy (ext_idx_set.begin(), ext_idx_set.end(), ext_idx.begin());

    // 2) declare id[i]=i for scatter
    std::vector<size_type> id (ext_idx.size());
    for (size_type i = 0; i < id.size(); i++) id[i] = i;

    // 3) init scatter
    distributor::tag_type tag_init = distributor::get_new_tag();
    mpi_scatter_init(
        ext_idx.size(),
        ext_idx.begin().operator->(),
        id.size(),
        id.begin().operator->(),
        ownership().dis_size(),
        ownership().begin().operator->(),
        tag_init,
        ownership().comm(),
        from,
        to);

    // 4) begin scatter: send local data to others and get ask for missing data
    std::vector<T,A> buffer (ext_idx.size());
    distributor::tag_type tag = distributor::get_new_tag();
    mpi_scatter_begin (
        begin().operator->(),
        buffer.begin().operator->(),
        from,
        to,
        set_op<T,T>(),
        tag,
        ownership().comm());

    // 5) end scatter: receive missing data
    mpi_scatter_end (
        begin().operator->(),
        buffer.begin(),
        from,
        to,
        set_op<T,T>(),
        tag,
        ownership().comm());

    // 6) build the associative container: pair (ext_idx ; data)
    for (size_type i = 0; i < buffer.size(); i++) {
        ext_idx_map.insert (std::make_pair (ext_idx[i], buffer[i]));
    }
}
template <class T, class A>
void
disarray_rep<T,distributed,A>::reset_dis_indexes() const
{
    std::set<size_type> ext_idx_set;
    for (typename scatter_map_type::const_iterator iter = _ext_x.begin(), last = _ext_x.end(); iter != last; iter++) {
      ext_idx_set.insert ((*iter).first);
    }
    set_dis_indexes (ext_idx_set);
}
/// @brief get values from ext_idx_set, that are managed by another proc; here T=container type
template <class T, class A>
template <class Set, class Map>
void
disarray_rep<T,distributed,A>::append_dis_entry (const Set& ext_idx_set, Map& ext_idx_map, std::true_type) const
{
    typedef typename T::value_type S; // S is supposed to be MPI-simple, i.e. with fixed size

    // 0) declare the local context
    typedef scatter_message<std::vector<T,A>, true>    message_type;
    
    message_type from;
    message_type to;

    // 1) convert set to vector, for direct acess:
    std::vector<size_type> ext_idx (ext_idx_set.size());
    std::copy (ext_idx_set.begin(), ext_idx_set.end(), ext_idx.begin());

    // 2) declare id[i]=i for scatter
    std::vector<size_type> id (ext_idx.size());
    for (size_type i = 0; i < id.size(); i++) id[i] = i;

    // 3) init scatter
    distributor::tag_type tag_init = distributor::get_new_tag();
    mpi_scatter_init(
        ext_idx.size(),
        ext_idx.begin().operator->(),
        id.size(),
        id.begin().operator->(),
        ownership().dis_size(),
        ownership().begin().operator->(),
        tag_init,
        ownership().comm(),
        from.get_base(),
        to.get_base());

    // 4) copy size of multi-valued objects into a tmp
    std::vector<size_type> data_sizes (size());
    for (size_type i = 0, n = size(); i < n; i++) {
      data_sizes[i] = base::operator[](i).size();
    }
    // 5) begin scatter: send local data to others and get ask for missing data
    std::vector<size_type> buffer (ext_idx.size());
    distributor::tag_type tag = distributor::get_new_tag();
    mpi_scatter_begin (
        data_sizes.begin().operator->(),
        buffer.begin().operator->(),
        from.get_base(),
        to.get_base(),
        set_op<size_type,size_type>(),
        tag,
        ownership().comm());

    // 6) end scatter: receive missing data
    mpi_scatter_end (
        data_sizes.begin().operator->(),
        buffer.begin(),
        from.get_base(),
        to.get_base(),
        set_op<size_type,size_type>(),
        tag,
        ownership().comm());

    // 7) initialize multi-valued scatter
    from.multi_init();
    to.multi_init();

    // 8) begin multi-valued scatter
    std::vector<T,A> multi_buffer (ext_idx.size());
    distributor::tag_type multi_tag = distributor::get_new_tag();
    mpi_scatter_begin (
        begin().operator->(),
        multi_buffer.begin().operator->(),
        from,
        to,
        set_op<T,T>(),
        multi_tag,
        ownership().comm());

    // 9) end scatter: receive missing data
    mpi_scatter_end (
        begin().operator->(),
        multi_buffer.begin(),
        from,
        to,
        set_op<T,T>(),
        multi_tag,
        ownership().comm());

    // 10) build the associative container: pair (ext_idx ; data)
    for (size_type i = 0; i < multi_buffer.size(); i++) {
        ext_idx_map.insert (std::make_pair (ext_idx[i], multi_buffer[i]));
    }
}
/// @brief get values from ext_idx_set, that are managed by another proc
template <class T, class A>
template <class Set, class Map>
inline
void
disarray_rep<T,distributed,A>::append_dis_entry (const Set& ext_idx_set, Map& ext_idx_map) const
{
    append_dis_entry (ext_idx_set, ext_idx_map, is_container());
}
template <class T, class A>
typename disarray_rep<T,distributed,A>::const_reference
disarray_rep<T,distributed,A>::dis_at (const size_type dis_i) const
{
    if (dis_i >= ownership().first_index() && dis_i < ownership().last_index()) {
      size_type i = dis_i - ownership().first_index();
      return base::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;
}
// ===============================================================
// put & get
// ===============================================================
template <class T, class A>
template <class PutFunction>
odiststream&
disarray_rep<T,distributed,A>::put_values (odiststream& ops, PutFunction put_element) const
{
    distributor::tag_type tag = distributor::get_new_tag();
    std::ostream& s = ops.os();
  
    // determine maximum message to arrive
    size_type max_size = 0;
    mpi::reduce(comm(), size(), max_size, mpi::maximum<size_type>(), 0);
 
    size_type io_proc = odiststream::io_proc();
    if (ownership().process() == io_proc) {
        for (size_type i = 0; i < size(); i++) {
            put_element (s, base::operator[](i));
            s << std::endl;
        }
        // receive and print messages
        std::vector<T,A> values (max_size);
        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, values.begin().operator->(), max_size);
  	    boost::optional<int> n_data_opt = status.count<T>();
	    check_macro (n_data_opt, "receive failed");
  	    size_type n_data = n_data_opt.get();
  	    for (size_type i = 0; i < n_data; i++) {
                put_element (s, values[i]);
                s << std::endl;
            }          
        }
        s << std::flush;
    } else {
	if (size() != 0) {
  	    comm().send(io_proc, tag, begin().operator->(), size());
	}
    }
    return ops;
}
template <class T, class A>
odiststream&
disarray_rep<T,distributed,A>::put_values (odiststream& ops) const
{
  return put_values (ops, _disarray_put_element_type<T>());
}
template <class T, class A>
odiststream&
disarray_rep<T,distributed,A>::put_matlab (odiststream& ops) const
{
  ops << "[";
  put_values (ops, _disarray_put_matlab_type<T>());
  return ops << "];";
}
template <class T, class A>
template <class PutFunction, class A2>
odiststream&
disarray_rep<T,distributed,A>::permuted_put_values (
  odiststream&                       ops,
  const disarray_rep<size_type,distributed,A2>& perm,
  PutFunction                        put_element) const
{
  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);
  disarray_rep<T,distributed,A>  perm_x (io_ownership);
  for (size_type i = 0, n = size(); i < n; i++) {
    perm_x.dis_entry (perm[i]) = base::operator[](i);
  }
  perm_x.template dis_entry_assembly_begin<typename default_set_op<T>::type>();
  perm_x.template dis_entry_assembly_end  <typename default_set_op<T>::type>();
  return perm_x.base::put_values (ops, put_element);
}
template <class T, class A>
template <class GetFunction>
idiststream&
disarray_rep<T,distributed,A>::get_values (idiststream& ps, GetFunction get_element) {
    distributor::tag_type tag = distributor::get_new_tag();
    std::istream& s = ps.is();
    size_type io_proc = odiststream::io_proc();
    if (ownership().process() == io_proc) {
        // load first chunk associated to proc 0
        if (!load_chunk (s, begin(), end(), get_element))
            error_macro("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
            size_type size_max = 1;
            for (size_type iproc = 0; iproc < ownership().n_process(); iproc++) {
                size_max = std::max (size_max, ownership().size(iproc));
            }
            std::vector<T,A> data_proc_j (size_max);
            T *start_j = data_proc_j.begin().operator->();
	    // 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;
                T *last_j = start_j + loc_sz_j;
                if (!load_chunk (s, start_j, last_j, get_element))
                    error_macro("read failed on input stream.");
  	        comm().send (jproc, tag, start_j, loc_sz_j);
            }
	}
  } else {
	if (size() != 0) {
	    comm().recv(io_proc, tag, begin().operator->(), size());	
	}
  }
  return ps;
}
template <class T, class A>
idiststream&
disarray_rep<T,distributed,A>::get_values (idiststream& ips) 
{
  return get_values (ips, _disarray_get_element_type<T>());
}
template <class T, class A>
void
disarray_rep<T,distributed,A>::dump (std::string name) const
{
  base::dump (name + itos(comm().rank()));
}

} // namespace rheolef
# endif // _RHEOLEF_HAVE_MPI