librheolef-dev 6.6-1build2 / usr / include / rheolef / array.h

#ifndef _RHEO_ARRAY_H
#define _RHEO_ARRAY_H
///
/// 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/distributed.h"
#include "rheolef/distributor.h"
#include "rheolef/diststream.h"
#include "rheolef/heap_allocator.h"
#include "rheolef/msg_util.h"
#include "rheolef/container_traits.h"
#ifdef _RHEOLEF_HAVE_MPI
#include "rheolef/mpi_pair_datatype.h"
#endif // _RHEOLEF_HAVE_MPI

#include <boost/mpl/bool.hpp>
#include <boost/type_traits/remove_const.hpp>

namespace rheolef {
/// @brief array element output helper
template <class T>
struct _array_put_element_type {
  std::ostream& operator() (std::ostream& os, const T& x) { return os << x; }
};
template <class T>
struct _array_put_matlab_type {
  std::ostream& operator() (std::ostream& os, const T& x) { return os << x << ";"; }
};
/// @brief array element input helper
template <class T>
struct _array_get_element_type {
  std::istream& operator() (std::istream& is, T& x) { return is >> x; }
};
} // namespace rheolef
// -------------------------------------------------------------
// the sequential representation
// -------------------------------------------------------------
namespace rheolef {

template <class T, class M, class A> class array_rep {};

template <class T, class A>
class array_rep<T,sequential,A> : public std::vector<T,A> {
public:
    typedef T                                         value_type;
    typedef A                                         allocator_type;
    typedef typename A::difference_type               difference_type;
    typedef std::vector<T,A>                          base;
    typedef typename base::size_type                  size_type;
    typedef typename base::iterator                   iterator;
    typedef typename base::const_iterator             const_iterator;
    typedef typename base::const_reference            const_reference;
    typedef typename base::reference                  reference;
    typedef reference	                              dis_reference;
    typedef distributor::communicator_type            communicator_type;
    typedef sequential                                memory_type;

    explicit array_rep (const A& alloc = A());
    array_rep (const distributor& ownership, const T& init_val = T(), const A& alloc = A());
    void resize   (const distributor& ownership, const T& init_val = T());
    array_rep (size_type loc_size = 0,       const T& init_val = T(), const A& alloc = A());
    void resize   (size_type loc_size = 0,       const T& init_val = T());
    array_rep (const array_rep<T,sequential,A>& x);

    A get_allocator() const { return base::get_allocator(); }
    size_type size() const { return base::size(); }
    iterator begin() { return base::begin(); }
    const_iterator begin() const { return base::begin(); }
    iterator end() { return base::end(); }
    const_iterator end() const { return base::end(); }
    const distributor& ownership() const { return _ownership; }

    reference       operator[] (size_type i)       { return base::operator[] (i); }
    const_reference operator[] (size_type i) const { return base::operator[] (i); }
    const_reference dis_at (size_type dis_i) const { return operator[] (dis_i); }
   
    size_type dis_size () const { return base::size(); }
    size_type first_index () const { return 0; }
    size_type last_index () const { return base::size(); }
    reference dis_entry     (size_type dis_i) { return operator[](dis_i); }
    void reset_dis_indexes() const {}
    template<class SetOp = typename default_set_op<T>::type>
    void dis_entry_assembly_begin (SetOp = SetOp()) {}
    template<class SetOp = typename default_set_op<T>::type>
    void dis_entry_assembly_end (SetOp = SetOp()) {}
    void repartition (					      // old_numbering for *this
        const array_rep<size_type,sequential,A>& partition,	      // old_ownership
        array_rep<T,sequential,A>&               new_array,	      // new_ownership (created)
        array_rep<size_type,sequential,A>&       old_numbering,      // new_ownership
        array_rep<size_type,sequential,A>&       new_numbering) const // old_ownership
    {
	error_macro ("not yet");
    }
    template<class A2>
    void reverse_permutation (                                // old_ownership for *this=iold2dis_inew
        array_rep<size_type,sequential,A2>& inew2dis_iold) const;    // new_ownership

    idiststream& get_values (idiststream& s);
    odiststream& put_values (odiststream& s) const;
    odiststream& put_matlab (odiststream& s) const;
    template <class GetFunction> idiststream& get_values (idiststream& ips, GetFunction get_element);
    template <class PutFunction> odiststream& put_values (odiststream& ops, PutFunction put_element) const;
    void dump (std::string name) const;
protected:
// data:
    distributor      _ownership;
};
// -------------------------------------------------------------
// the distributed representation
// -------------------------------------------------------------
#ifdef _RHEOLEF_HAVE_MPI
template <class T, class A>
class array_rep<T,distributed,A> : public array_rep<T,sequential,A> {
public:

// typedefs:

    typedef array_rep<T,sequential,A>     	base;
    typedef typename base::value_type     	value_type;
    typedef typename base::size_type      	size_type;
    typedef typename base::difference_type      difference_type;
    typedef typename base::reference      	reference;
    typedef typename base::const_reference      const_reference;
    typedef typename base::iterator       	iterator;
    typedef typename base::const_iterator 	const_iterator;
    typedef distributor::communicator_type      communicator_type;
    typedef distributed                         memory_type;
    typedef std::map <size_type, T, std::less<size_type>, heap_allocator<std::pair<size_type,T> > > 
						scatter_map_type;

    struct dis_reference {
      dis_reference (array_rep<T,distributed,A>& x, size_type dis_i)
       : _x(x), _dis_i(dis_i) {}

      dis_reference& operator= (const T& value) {
        _x.set_dis_entry (_dis_i, value);
        return *this;
      }
      template<class U>
      dis_reference& operator+= (const U& value) {
        _x.set_add_dis_entry (_dis_i, value);
        return *this;
      }
    // data:
    protected:
      array_rep<T,distributed,A>& _x;
      size_type           _dis_i;
    };

// allocators:

    array_rep (const distributor& ownership, const T&  init_val = T(), const A& alloc = A());
    void resize   (const distributor& ownership, const T&  init_val = T());
    array_rep (const array_rep<T,distributed,A>& x);

    A get_allocator() const        { return base::get_allocator(); }
    size_type size() const         { return base::size(); }
    const_iterator begin() const   { return base::begin(); }
    const_iterator end() const     { return base::end(); }
    iterator begin()               { return base::begin(); }
    iterator end()                 { return base::end(); }

    const distributor& ownership() const  { return base::_ownership; }
    const mpi::communicator& comm() const { return ownership().comm(); }
    size_type first_index () const        { return ownership().first_index(); }
    size_type last_index () const         { return ownership().last_index(); }
    size_type dis_size () const           { return ownership().dis_size(); }

    dis_reference dis_entry (size_type dis_i) { return dis_reference (*this, dis_i); }

    template<class SetOp = typename default_set_op<T>::type>
    void dis_entry_assembly_begin (SetOp my_set_op = SetOp());
    template<class SetOp = typename default_set_op<T>::type>
    void dis_entry_assembly_end   (SetOp my_set_op = SetOp());
    template<class SetOp = typename default_set_op<T>::type>
    void dis_entry_assembly       (SetOp my_set_op = SetOp())
   		{ dis_entry_assembly_begin (my_set_op); dis_entry_assembly_end (my_set_op); }

    template<class Set, class Map>
    void append_dis_entry (const Set& ext_idx_set, Map& ext_idx_map) const;

    template<class Set, class Map>
    void get_dis_entry (const Set& ext_idx_set, Map& ext_idx_map) const {
    	    ext_idx_map.clear();
    	    append_dis_entry (ext_idx_set, ext_idx_map);
	}

    template<class Set>
    void append_dis_indexes (const Set& ext_idx_set) const { append_dis_entry (ext_idx_set, _ext_x); }

    template<class Set>
    void set_dis_indexes    (const Set& ext_idx_set) const { get_dis_entry    (ext_idx_set, _ext_x); }
    void reset_dis_indexes() const;

    const_reference dis_at (size_type dis_i) const;

    // get all external pairs (dis_i, values):
    const scatter_map_type& get_dis_map_entries() const { return _ext_x; }

    template<class A2>
    void repartition (					      // old_numbering for *this
        const array_rep<size_type,distributed,A2>& partition,	      // old_ownership
        array_rep<T,distributed,A>&                new_array,	      // new_ownership (created)
        array_rep<size_type,distributed,A2>&       old_numbering,     // new_ownership
        array_rep<size_type,distributed,A2>&       new_numbering) const; // old_ownership

    template<class A2>
    void permutation_apply (		 		       // old_numbering for *this
        const array_rep<size_type,distributed,A2>& new_numbering,      // old_ownership
        array_rep<T,distributed,A>&                new_array) const;   // new_ownership (already allocated)

    template<class A2>
    void reverse_permutation (                                // old_ownership for *this=iold2dis_inew
        array_rep<size_type,distributed,A2>& inew2dis_iold) const;    // new_ownership

    idiststream& get_values (idiststream& s);
    odiststream& put_values (odiststream& s) const;
    odiststream& put_matlab (odiststream& s) const;
    template <class GetFunction> idiststream& get_values (idiststream& ips, GetFunction get_element);
    template <class PutFunction> odiststream& put_values (odiststream& ops, PutFunction put_element) const;
    template <class PutFunction, class A2> odiststream& permuted_put_values (odiststream& ops, const array_rep<size_type,distributed,A2>& perm,
		PutFunction put_element) const;
    void dump (std::string name) const;
protected:
    void set_dis_entry (size_type dis_i, const T& val);
    template<class U>
    void set_add_dis_entry (size_type dis_i, const U& val);
// typedefs:
    /** 1) stash: store data before assembly() communications:
      *   select multimap<U> when T=set<U> and map<T> otherwise
      */
    template<class Pair>
    struct remove_const_in_pair {
        typedef Pair type;
    };
    template<class T1, class T2>
    struct remove_const_in_pair<std::pair<T1,T2> > {
        typedef std::pair<typename boost::remove_const<T1>::type,
                          typename boost::remove_const<T2>::type> type;
    };
    template<class U, class IsContainer> struct stash_traits {};
    template<class U>
    struct stash_traits<U,boost::mpl::false_> {
        typedef U mapped_type;
        typedef std::map <size_type, U, std::less<size_type>, heap_allocator<std::pair<size_type,U> > >   map_type;
    };
    template<class U>
    struct stash_traits<U,boost::mpl::true_> {
        typedef typename remove_const_in_pair<typename U::value_type>::type mapped_type;
        typedef std::multimap <size_type, mapped_type, std::less<size_type>, heap_allocator<std::pair<size_type,mapped_type> > >   map_type;
    };
    typedef typename is_container_of_mpi_datatype<T>::type is_container;
    typedef typename stash_traits<T,is_container>::mapped_type stash_value;
    typedef typename stash_traits<T,is_container>::map_type    stash_map_type;

    /** 2) message: for communication during assembly_begin(), assembly_end()
      */
    struct message_type {
        std::list<std::pair<size_type,mpi::request>,A>    waits;
        std::vector<std::pair<size_type,stash_value>,A>   data;
    };
    /** 3) scatter (get_entry): specialized versions for T=container and T=simple type
      */
    template<class Set, class Map>
    void append_dis_entry (const Set& ext_idx_set, Map& ext_idx_map, boost::mpl::true_) const;
    template<class Set, class Map>
    void append_dis_entry (const Set& ext_idx_set, Map& ext_idx_map, boost::mpl::false_) const;

// data:
    stash_map_type   _stash;		// for assembly msgs:
    message_type     _send;
    message_type     _receive;
    size_type        _receive_max_size;
    mutable scatter_map_type _ext_x;		// for ext values (scatter)
};
#endif // _RHEOLEF_HAVE_MPI
// -------------------------------------------------------------
// the basic class with a smart pointer to representation
// the user-level class with memory-model parameter
// -------------------------------------------------------------
/*Class:array
NAME:  @code{array} - container in distributed environment (@PACKAGE@-@VERSION@)
SYNOPSYS:       
  STL-like vector container for a distributed memory machine model.
EXAMPLE:
   A sample usage of the class is:
   @example
     int main(int argc, char**argv) @{
        environment distributed(argc, argv);
        array<double> x(distributor(100), 3.14);
        dout << x << endl;
     @}
   @end example
   The array<T> interface is similar to those of the std::vector<T> with the 
   addition of some communication features in the distributed case:
   write accesses with entry/assembly and read access with dis_at.

DISTRIBUTED WRITE ACCESS:
   Loop on any @code{dis_i} that is not managed by the current processor:
   @example
	x.dis_entry (dis_i) = value;
   @end example
   and then, after loop, perform all communication:
   @example
	x.dis_entry_assembly();
   @end example
   After this command, each value is stored in the array, available the processor
   associated to @code{dis_i}.
DISTRIBUTED READ ACCESS:
   First, define the set of indexes:
   @example
 	std::set<size_t> ext_idx_set; 
   @end example
   Then, loop on @code{dis_i} indexes that are not managed by the current processor:
   @example
	ext_idx_set.insert (dis_i);
   @end example
   After the loop, performs the communications:
   @example
        x.set_dis_indexes (ext_idx_set);
   @end example
   After this command, each values associated to the @code{dis_i} index,
   and that belongs to the index set, is now available also on the
   current processor as:
   @example
        value = x.dis_at (dis_i);
   @end example
   For convenience, if @code{dis_i} is managed by the current processor, this
   function returns also the value.
NOTE:
  The class takes two template parameters: one for the type T and the second
  for the memory model M, that could be either M=distributed or M=sequential.
  The two cases are associated to two diferent implementations, but proposes
  exactly the same interface. The sequential interface propose also a supplementary
  constructor:
   @example
        array<double,sequential> x(local_size, init_val);
   @end example
   This constructor is a STL-like one but could be consufused in the distributed case,
   since there are two sizes: a local one and a global one. In that case, the use
   of the distributor, as a generalization of the size concept, clarify the situation
   (@pxref{distributor class}).
 
IMPLEMENTATION NOTE:
  "scatter" via "get_dis_entry".

  "gather" via "dis_entry(dis_i) = value"
  or "dis_entry(dis_i) += value". Note that += applies when T=idx_set where
  idx_set is a wrapper class of std::set<size_t> ; the += operator represents the
  union of a set. The operator= is used when T=double or others simple T types
  without algebra. If there is a conflict, i.e. several processes set the dis_i
  index, then the result of operator+= depends upon the order of the process at
  each run and is not deterministic. Such ambiguous behavior is not detected
  yet at run time.

AUTHOR: Pierre.Saramito@imag.fr
End:
*/
template <class T, class M = rheo_default_memory_model, class A = std::allocator<T> >
class array {
public:
    typedef M memory_type;
    typedef typename std::vector<T,A>::size_type      size_type;
    typedef typename std::vector<T,A>::iterator       iterator;
    typedef typename std::vector<T,A>::const_iterator const_iterator;
};
//<verbatim:
template <class T, class A>
class array<T,sequential,A> : public smart_pointer<array_rep<T,sequential,A> > {
public:

// typedefs:

    typedef array_rep<T,sequential,A>  	  rep;
    typedef smart_pointer<rep> 		  base;

    typedef sequential 			  memory_type;
    typedef typename rep::size_type 	  size_type;
    typedef typename rep::difference_type difference_type;
    typedef typename rep::value_type 	  value_type;
    typedef typename rep::reference 	  reference;
    typedef typename rep::dis_reference   dis_reference;
    typedef typename rep::iterator 	  iterator;
    typedef typename rep::const_reference const_reference;
    typedef typename rep::const_iterator  const_iterator;

// allocators:


    array       (size_type loc_size = 0,       const T& init_val = T(), const A& alloc = A());
    void resize (size_type loc_size = 0,       const T& init_val = T());
    array       (const distributor& ownership, const T& init_val = T(), const A& alloc = A());
    void resize (const distributor& ownership, const T& init_val = T());

// local accessors & modifiers:

    A get_allocator() const              { return base::data().get_allocator(); }
    size_type     size () const          { return base::data().size(); }
    size_type dis_size () const          { return base::data().dis_size(); }
    const distributor& ownership() const { return base::data().ownership(); }
    const communicator& comm() const     { return ownership().comm(); }

    reference       operator[] (size_type i)       { return base::data().operator[] (i); }
    const_reference operator[] (size_type i) const { return base::data().operator[] (i); }
    reference       operator() (size_type i)       { return base::data().operator[] (i); }
    const_reference operator() (size_type i) const { return base::data().operator[] (i); }
    const_reference dis_at (size_type dis_i) const { return operator[] (dis_i); }

          iterator begin()       { return base::data().begin(); }
    const_iterator begin() const { return base::data().begin(); }
          iterator end()         { return base::data().end(); }
    const_iterator end() const   { return base::data().end(); }

// global modifiers (for compatibility with distributed interface):

    dis_reference dis_entry (size_type dis_i) { return base::data().dis_entry(dis_i); }
    template<class SetOp = typename default_set_op<T>::type>
    void dis_entry_assembly (SetOp my_set_op = SetOp()) {}
    template<class SetOp = typename default_set_op<T>::type>
    void dis_entry_assembly_begin (SetOp my_set_op = SetOp()) {}
    template<class SetOp = typename default_set_op<T>::type>
    void dis_entry_assembly_end (SetOp my_set_op = SetOp()) {}

    void dis_entry_assembly_begin() {}
    void dis_entry_assembly_end()   {}
    void dis_entry_assembly()       {}

    void reset_dis_indexes() const {}
    template<class Set> void set_dis_indexes    (const Set& ext_idx_set) const {}
    template<class Set> void append_dis_indexes (const Set& ext_idx_set) const {}
    template<class Set, class Map> void append_dis_entry (const Set& ext_idx_set, Map& ext_idx_map) const {}
    template<class Set, class Map> void get_dis_entry    (const Set& ext_idx_set, Map& ext_idx_map) const {}

// apply a partition:
 
    template<class RepSize>
    void repartition (			             // old_numbering for *this
        const RepSize&         partition,	     // old_ownership
        array<T,sequential,A>& new_array,	     // new_ownership (created)
        RepSize&               old_numbering,	     // new_ownership
        RepSize&               new_numbering) const  // old_ownership
        { return base::data().repartition (partition, new_array, old_numbering, new_numbering); }

    template<class RepSize>
    void permutation_apply (                       // old_numbering for *this
        const RepSize&          new_numbering,     // old_ownership
        array<T,sequential,A>&  new_array) const   // new_ownership (already allocated)
        { return base::data().permutation_apply (new_numbering, new_array); }

    void reverse_permutation (                                 // old_ownership for *this=iold2dis_inew
        array<size_type,sequential,A>& inew2dis_iold) const   // new_ownership
        { base::data().reverse_permutation (inew2dis_iold.data()); }

// i/o:

    odiststream& put_values (odiststream& ops) const { return base::data().put_values(ops); }
    idiststream& get_values (idiststream& ips)       { return base::data().get_values(ips); }
    template <class GetFunction>
    idiststream& get_values (idiststream& ips, GetFunction get_element)       { return base::data().get_values(ips, get_element); }
    template <class PutFunction>
    odiststream& put_values (odiststream& ops, PutFunction put_element) const { return base::data().put_values(ops, put_element); }
    void dump (std::string name) const { return base::data().dump(name); }
};
//>verbatim:
template <class T, class A>
inline
array<T,sequential,A>::array (
    	size_type loc_size,
	const T&  init_val,
        const A&  alloc)
 : base(new_macro(rep(loc_size,init_val,alloc)))
{
}
template <class T, class A>
inline
array<T,sequential,A>::array (
    	const distributor& ownership,
	const T&           init_val,
        const A&           alloc)
 : base(new_macro(rep(ownership,init_val,alloc)))
{
}
template <class T, class A>
inline
void
array<T,sequential,A>::resize (
    	size_type loc_size,
	const T&  init_val)
{
  base::data().resize (loc_size,init_val);
}
template <class T, class A>
inline
void
array<T,sequential,A>::resize (
    	const distributor& ownership,
	const T&           init_val)
{
  base::data().resize (ownership,init_val);
}
#ifdef _RHEOLEF_HAVE_MPI
//<verbatim:
template <class T, class A>
class array<T,distributed,A> : public smart_pointer<array_rep<T,distributed,A> > {
public:

// typedefs:

    typedef array_rep<T,distributed,A>    rep;
    typedef smart_pointer<rep> 		  base;

    typedef distributed 		  memory_type;
    typedef typename rep::size_type 	  size_type;
    typedef typename rep::difference_type difference_type;
    typedef typename rep::value_type 	  value_type;
    typedef typename rep::reference 	  reference;
    typedef typename rep::dis_reference   dis_reference;
    typedef typename rep::iterator 	  iterator;
    typedef typename rep::const_reference const_reference;
    typedef typename rep::const_iterator  const_iterator;
    typedef typename rep::scatter_map_type scatter_map_type;

// allocators:

    array       (const distributor& ownership = distributor(), const T& init_val = T(), const A& alloc = A());
    void resize (const distributor& ownership = distributor(), const T& init_val = T());

// local accessors & modifiers:

    A get_allocator() const              { return base::data().get_allocator(); }
    size_type     size () const          { return base::data().size(); }
    size_type dis_size () const          { return base::data().dis_size(); }
    const distributor& ownership() const { return base::data().ownership(); }
    const communicator& comm() const     { return base::data().comm(); }

    reference       operator[] (size_type i)       { return base::data().operator[] (i); }
    const_reference operator[] (size_type i) const { return base::data().operator[] (i); }
    reference       operator() (size_type i)       { return base::data().operator[] (i); }
    const_reference operator() (size_type i) const { return base::data().operator[] (i); }

          iterator begin()       { return base::data().begin(); }
    const_iterator begin() const { return base::data().begin(); }
          iterator end()         { return base::data().end(); }
    const_iterator end() const   { return base::data().end(); }

// global accessor:

    template<class Set, class Map>
    void append_dis_entry (const Set& ext_idx_set, Map& ext_idx_map) const { base::data().append_dis_entry (ext_idx_set, ext_idx_map); }

    template<class Set, class Map>
    void get_dis_entry    (const Set& ext_idx_set, Map& ext_idx_map) const { base::data().get_dis_entry (ext_idx_set, ext_idx_map); }

    template<class Set>
    void append_dis_indexes (const Set& ext_idx_set) const { base::data().append_dis_indexes (ext_idx_set); }
    void reset_dis_indexes() const { base::data().reset_dis_indexes(); }

    template<class Set>
    void set_dis_indexes    (const Set& ext_idx_set) const { base::data().set_dis_indexes (ext_idx_set); }

    const T& dis_at (size_type dis_i) const { return base::data().dis_at (dis_i); }

    // get all external pairs (dis_i, values):
    const scatter_map_type& get_dis_map_entries() const { return base::data().get_dis_map_entries(); }

// global modifiers (for compatibility with distributed interface):

    dis_reference dis_entry (size_type dis_i) { return base::data().dis_entry(dis_i); }

    template<class SetOp = typename default_set_op<T>::type>
    void dis_entry_assembly_begin (SetOp my_set_op = SetOp()) { base::data().dis_entry_assembly_begin (my_set_op); }
    template<class SetOp = typename default_set_op<T>::type>
    void dis_entry_assembly_end   (SetOp my_set_op = SetOp()) { base::data().dis_entry_assembly_end   (my_set_op); }
    template<class SetOp = typename default_set_op<T>::type>
    void dis_entry_assembly       (SetOp my_set_op = SetOp()) { base::data().dis_entry_assembly       (my_set_op); }

    void dis_entry_assembly_begin() { base::data().template dis_entry_assembly_begin<typename default_set_op<T>::type>(); }
    void dis_entry_assembly_end()   { base::data().template dis_entry_assembly_end<typename default_set_op<T>::type>(); }
    void dis_entry_assembly()       { dis_entry_assembly_begin(); dis_entry_assembly_end(); }

// apply a partition:
 
    template<class RepSize>
    void repartition (			            // old_numbering for *this
        const RepSize&        partition,            // old_ownership
        array<T,distributed>& new_array,            // new_ownership (created)
        RepSize&              old_numbering,        // new_ownership
        RepSize&              new_numbering) const  // old_ownership
        { return base::data().repartition (partition.data(), new_array.data(), old_numbering.data(), new_numbering.data()); }

    template<class RepSize>
    void permutation_apply (                       // old_numbering for *this
        const RepSize&          new_numbering,     // old_ownership
        array<T,distributed,A>& new_array) const   // new_ownership (already allocated)
        { base::data().permutation_apply (new_numbering.data(), new_array.data()); }

    void reverse_permutation (                                 // old_ownership for *this=iold2dis_inew
        array<size_type,distributed,A>& inew2dis_iold) const   // new_ownership
        { base::data().reverse_permutation (inew2dis_iold.data()); }

// i/o:

    odiststream& put_values (odiststream& ops) const { return base::data().put_values(ops); }
    idiststream& get_values (idiststream& ips)       { return base::data().get_values(ips); }
    void dump (std::string name) const 	    { return base::data().dump(name); }

    template <class GetFunction>
    idiststream& get_values (idiststream& ips, GetFunction get_element)       { return base::data().get_values(ips, get_element); }
    template <class PutFunction>
    odiststream& put_values (odiststream& ops, PutFunction put_element) const { return base::data().put_values(ops, put_element); }
    template <class PutFunction, class A2> odiststream& permuted_put_values (
		odiststream& ops, const array<size_type,distributed,A2>& perm, PutFunction put_element) const 
								     { return base::data().permuted_put_values (ops, perm.data(), put_element); }
};
//>verbatim:
template <class T, class A>
inline
array<T,distributed,A>::array (
    	const distributor& ownership,
	const T&           init_val,
        const A&           alloc)
 : base(new_macro(rep(ownership,init_val,alloc)))
{
}
template <class T, class A>
inline
void
array<T,distributed,A>::resize (
    	const distributor& ownership,
	const T         &  init_val)
{
  base::data().resize (ownership,init_val);
}
#endif // _RHEOLEF_HAVE_MPI

// -------------------------------------------------------------
// i/o with operator<< & >>
// -------------------------------------------------------------
template <class T, class A>
inline
idiststream&
operator >> (idiststream& ips,  array<T,sequential,A>& x)
{ 
    return x.get_values(ips); 
}
template <class T, class A>
inline
odiststream&
operator << (odiststream& ops, const array<T,sequential,A>& x)
{
    return x.put_values(ops);
}
#ifdef _RHEOLEF_HAVE_MPI
template <class T, class A>
inline
idiststream&
operator >> (idiststream& ips,  array<T,distributed,A>& x)
{ 
    return x.get_values(ips); 
}
template <class T, class A>
inline
odiststream&
operator << (odiststream& ops, const array<T,distributed,A>& x)
{
    return x.put_values(ops);
}
#endif // _RHEOLEF_HAVE_MPI
} // namespace rheolef
// -------------------------------------------------------------
// not inlined : longer code
// -------------------------------------------------------------
#include "rheolef/array_seq.icc"
#include "rheolef/array_mpi.icc"
#endif // _RHEO_ARRAY_H