librheolef-dev 6.5-1+b1 / usr / include / rheolef / csr.h

# ifndef _RHEOLEF_CSR_H
# define _RHEOLEF_CSR_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/vec.h"
#include "rheolef/asr.h"
#include "rheolef/vector_of_iterator.h"
#include "rheolef/scatter_message.h"
#include "rheolef/pair_util.h"

#include <boost/bind.hpp>

namespace rheolef {

// -------------------------------------------------------------
// the sequential representation
// -------------------------------------------------------------
template<class T, class M> class csr_rep {};

template<class T>
class csr_rep<T,sequential> : public vector_of_iterator<std::pair<typename std::vector<T>::size_type,T> > {
public:
    typedef typename std::vector<T>::size_type                            size_type;
    typedef T                                                             element_type;
    typedef sequential                                                    memory_type;
    typedef typename std::pair<size_type,T>                               pair_type;
    typedef typename vector_of_iterator<pair_type>::iterator              iterator;
    typedef typename vector_of_iterator<pair_type>::const_iterator        const_iterator;
    typedef typename vector_of_iterator<pair_type>::difference_type       difference_type;
    typedef typename vector_of_iterator<pair_type>::value_type            data_iterator;
    typedef typename vector_of_iterator<pair_type>::const_value_type      const_data_iterator;

    csr_rep (size_type loc_nrow1 = 0, size_type loc_ncol1 = 0, size_type loc_nnz1 = 0);
    void resize (size_type loc_nrow1 = 0, size_type loc_ncol1 = 0, size_type loc_nnz1 = 0);
    csr_rep (const distributor& row_ownership, const distributor& col_ownership, size_type nnz1 = 0);
    void resize (const distributor& row_ownership, const distributor& col_ownership, size_type nnz1 = 0);
    csr_rep (const csr_rep<T,sequential>& a);
    template<class A> void build_from_asr (const asr<T,sequential,A>& a);
    template<class A> explicit csr_rep    (const asr<T,sequential,A>& a);
    template<class A> void build_from_diag (const array_rep<T,sequential,A>& d);

    const distributor& row_ownership() const { return _row_ownership; }
    const distributor& col_ownership() const { return _col_ownership; }
    const_iterator begin() const { return vector_of_iterator<pair_type>::begin(); }
    const_iterator end()   const { return vector_of_iterator<pair_type>::end(); }
          iterator begin()       { return vector_of_iterator<pair_type>::begin(); }
          iterator end()         { return vector_of_iterator<pair_type>::end(); }
    size_type nrow() const       { return vector_of_iterator<pair_type>::size()-1; }
    size_type ncol() const { return _col_ownership.size(); }
    size_type nnz() const { return _data.size(); }
    size_type dis_nrow() const { return nrow(); }
    size_type dis_ncol() const { return ncol(); }
    size_type dis_nnz() const { return nnz(); }
    T max_abs () const;
    bool is_symmetric() const { return _is_symmetric; }
    void set_symmetry (bool is_symm) const { _is_symmetric = is_symm; }
    void set_symmetry_by_check (const T& tol = std::numeric_limits<T>::epsilon()) const;
    size_type pattern_dimension() const { return _pattern_dimension; }
    void set_pattern_dimension(size_type dim) const { _pattern_dimension = dim; }
    size_type row_first_index () const { return 0; }
    size_type row_last_index () const { return nrow(); }
    size_type col_first_index () const { return 0; }
    size_type col_last_index () const { return ncol(); }
    idiststream& get (idiststream&);
    odiststream& put (odiststream&, size_type istart = 0) const;
    odiststream& put_matrix_market (odiststream&, size_type istart = 0) const;
    odiststream& put_sparse_matlab (odiststream&, size_type istart = 0) const;
    void dump (const std::string& name, size_type istart = 0) const;
    void mult       (const vec<T,sequential>& x, vec<T,sequential>& y) const;
    void trans_mult (const vec<T,sequential>& x, vec<T,sequential>& y) const;
    csr_rep<T,sequential>& operator*= (const T& lambda);
    template <class BinaryOp>
    void assign_add (const csr_rep<T,sequential>& a, const csr_rep<T,sequential>& b, BinaryOp binop);
    void build_transpose (csr_rep<T,sequential>& b) const;
    void assign_mult (const csr_rep<T,sequential>& a, const csr_rep<T,sequential>& b);

    // accessors, only for distributed (for interface compatibility)
    size_type ext_nnz() const                { return 0; }
    size_type dis_ext_nnz() const            { return 0; }
    const_iterator ext_begin() const         { return const_iterator(); }
    const_iterator ext_end()   const         { return const_iterator(); }
    size_type jext2dis_j (size_type jext) const { return 0; }

//protected:
    distributor          		                          _row_ownership;
    distributor          		                          _col_ownership;
    std::vector<std::pair<typename std::vector<T>::size_type,T> > _data;
    mutable bool						  _is_symmetric;
    mutable size_type						  _pattern_dimension; // e.g. FEM 3d-pattern
};
template<class T>
template<class A>
inline
csr_rep<T,sequential>::csr_rep(const asr<T,sequential,A>& a)
  : vector_of_iterator<pair_type>(a.nrow()+1),
   _row_ownership (a.row_ownership()),
   _col_ownership (a.col_ownership()),
   _data(a.nnz())
{
  build_from_asr (a);
}
template<class T>
inline
idiststream&
csr_rep<T,sequential>::get (idiststream& ids)
{
  typedef std::allocator<T> A; // TODO: use heap_allocator for asr
  asr<T,sequential,A> a;
  a.get(ids);
  build_from_asr (a);
  return ids;
}
// -------------------------------------------------------------
// the distributed representation
// -------------------------------------------------------------
#ifdef _RHEOLEF_HAVE_MPI
template<class T>
class csr_rep<T,distributed> : public csr_rep<T,sequential> {
public:
    typedef csr_rep<T,sequential>              base;
    typedef distributed                        memory_type;
    typedef typename base::size_type           size_type;
    typedef typename base::element_type        element_type;
    typedef typename base::iterator            iterator;
    typedef typename base::const_iterator      const_iterator;
    typedef typename base::data_iterator       data_iterator;
    typedef typename base::const_data_iterator const_data_iterator;

    csr_rep ();
    csr_rep (const csr_rep<T,distributed>& a);
    template<class A> explicit csr_rep (const asr<T,distributed,A>& a);
    template<class A> void build_from_asr  (const asr<T,distributed,A>& a);
    void resize (const distributor& row_ownership, const distributor& col_ownership, size_type nnz1 = 0);
    template<class A> void build_from_diag (const array_rep<T,distributed,A>& d);

    const distributor& row_ownership() const { return base::_row_ownership; }
    const distributor& col_ownership() const { return base::_col_ownership; }
    const communicator& comm() const { return row_ownership().comm(); }
    const_iterator begin() const { return base::begin(); }
    const_iterator end()   const { return base::end(); }
          iterator begin()       { return base::begin(); }
          iterator end()         { return base::end(); }
    size_type      ext_nnz()   const { return _ext.nnz(); }
    const_iterator ext_begin() const { return _ext.begin(); }
    const_iterator ext_end()   const { return _ext.end(); }
          iterator ext_begin()       { return _ext.begin(); }
          iterator ext_end()         { return _ext.end(); }
    size_type nrow() const { return base::nrow(); }
    size_type ncol() const { return base::ncol(); }
    size_type nnz() const { return base::nnz(); }
    size_type dis_nrow() const { return row_ownership().dis_size(); }
    size_type dis_ncol() const { return col_ownership().dis_size(); }
    size_type dis_nnz()     const { return _dis_nnz; }
    size_type dis_ext_nnz() const { return _dis_ext_nnz; }
    T max_abs () const;
    bool is_symmetric() const { return base::is_symmetric(); }
    void set_symmetry (bool is_symm) const { base::set_symmetry(is_symm); }
    void set_symmetry_by_check (const T& tol = std::numeric_limits<T>::epsilon()) const;
    size_type pattern_dimension() const { return base::pattern_dimension(); }
    void set_pattern_dimension(size_type dim) const { base::set_pattern_dimension(dim); }
    size_type row_first_index () const { return row_ownership().first_index(); }
    size_type row_last_index ()  const { return row_ownership().last_index(); }
    size_type col_first_index () const { return col_ownership().first_index(); }
    size_type col_last_index ()  const { return col_ownership().last_index(); }
    size_type jext2dis_j (size_type jext)  const;
    idiststream& get (idiststream&);
    odiststream& put (odiststream&) const;
    void dump (const std::string& name) const;
    void mult       (const vec<T,distributed>& x, vec<T,distributed>& y) const;
    void trans_mult (const vec<T,distributed>& x, vec<T,distributed>& y) const;
    csr_rep<T,distributed>& operator*= (const T& lambda);
    template <class BinaryOp>
    void assign_add (const csr_rep<T,distributed>& a, const csr_rep<T,distributed>& b, BinaryOp binop);
    void build_transpose (csr_rep<T,distributed>& b) const;
    void assign_mult (const csr_rep<T,distributed>& a, const csr_rep<T,distributed>& b);
protected:
// data:
    // diagonal part is the basic csr_rep<seq> type
    // extra-diagonal blocs are sequential csr also:
    csr_rep<T,sequential>      		     _ext;
    std::vector<size_type>    		     _jext2dis_j;
    size_type            		     _dis_nnz;
    size_type            		     _dis_ext_nnz;
    
    // A*x internal stuff: scatter and buffer (lazy initialization):
    mutable bool                             _scatter_initialized;
    mutable scatter_message<std::vector<T> > _from;
    mutable scatter_message<std::vector<T> > _to;
    mutable std::vector<T>    		     _buffer;
// internal:
    void _scatter_init() const;
    void _scatter_init_guard() const {
	if (_scatter_initialized) return;
	_scatter_initialized = true;
        _scatter_init();
    }
};
template<class T>
inline
typename csr_rep<T,distributed>::size_type
csr_rep<T,distributed>::jext2dis_j (size_type jext)  const
{
  check_macro (jext < _jext2dis_j.size(), "jext2dis_j: jext="<<jext<<" is out of range [0:"<<_jext2dis_j.size()<<"[");
  return _jext2dis_j [jext];
}
template<class T>
template<class A>
inline
csr_rep<T,distributed>::csr_rep(const asr<T,distributed,A>& a)
  : csr_rep<T,sequential>(),
    _ext (),
    _jext2dis_j(),
    _dis_nnz(0),
    _dis_ext_nnz(0),
    _scatter_initialized(false),
    _from(),
    _to(),
    _buffer()
{
    build_from_asr (a);
}
template<class T>
inline
idiststream&
csr_rep<T,distributed>::get (idiststream& ips)
{
    typedef std::allocator<T> A; // TODO: use heap_alloc for asr
    asr<T,distributed,A> a;
    a.get (ips);
    build_from_asr (a);
    return ips;
}
#endif // _RHEOLEF_HAVE_MPI

// these classes are used for allocator from the std::initializer_list
template <class T, class M> class csr_concat_value;
template <class T, class M> class csr_concat_line;

// -------------------------------------------------------------
// the basic class with a smart pointer to representation
// the user-level class with memory-model parameter
// -------------------------------------------------------------
/*Class:csr
NAME:  @code{csr} - compressed sparse row matrix (@PACKAGE@-@VERSION@)
SYNOPSYS:       
  Distributed compressed sparse matrix container stored row by row.
DESCRIPTION:
  Sparse matrix are compressed by rows. In distributed environment, the
  distribution follows the row distributor (see @ref{distributor class}).
ALGEBRA:
  Adding or substracting two matrices writes @code{a+b} and @code{a-b}, respectively,
  and multiplying a matrix by a scalar writes @code{lambda*x}.
  Thus, any linear combination of sparse matrices is available.
 
  Matrix-vector product writes @code{a*x} where @code{x} is a vector (see @ref{vec class}).
LIMITATIONS:
  Some basic linear algebra is still under development:
  @code{a.trans_mult(x)} matrix transpose vector product,
  @code{trans(a)} matrix transpose,
  @code{a*b} matrix product.
AUTHORS: Pierre.Saramito@imag.fr
DATE:   10 february 1999
METHODS: @csr
End:
*/
template <class T, class M = rheo_default_memory_model>
class csr {
public:
    typedef M memory_type;
};
//<verbatim:
template<class T>
class csr<T,sequential> : public smart_pointer<csr_rep<T,sequential> > {
public:

// typedefs:

    typedef csr_rep<T,sequential>             rep;
    typedef smart_pointer<rep>                base;
    typedef typename rep::memory_type         memory_type;
    typedef typename rep::size_type           size_type;
    typedef typename rep::element_type        element_type;
    typedef typename rep::iterator            iterator;
    typedef typename rep::const_iterator      const_iterator;
    typedef typename rep::data_iterator       data_iterator;
    typedef typename rep::const_data_iterator const_data_iterator;

// allocators/deallocators:

    csr() : base(new_macro(rep())) {}
    template<class A>
    explicit csr(const asr<T,sequential,A>& a) : base(new_macro(rep(a))) {}
    void resize (size_type loc_nrow1 = 0, size_type loc_ncol1 = 0, size_type loc_nnz1 = 0)
        { base::data().resize(loc_nrow1, loc_ncol1, loc_nnz1); }
    void resize (const distributor& row_ownership, const distributor& col_ownership, size_type nnz1 = 0)
        { base::data().resize(row_ownership, col_ownership, nnz1); }

// allocators from initializer list (c++ 2011):

#ifdef _RHEOLEF_HAVE_STD_INITIALIZER_LIST
    csr (const std::initializer_list<csr_concat_value<T,sequential> >& init_list);
    csr (const std::initializer_list<csr_concat_line<T,sequential> >&  init_list);
#endif // _RHEOLEF_HAVE_STD_INITIALIZER_LIST

// accessors:

    // global sizes
    const distributor& row_ownership() const { return base::data().row_ownership(); }
    const distributor& col_ownership() const { return base::data().col_ownership(); }
    size_type dis_nrow () const              { return row_ownership().dis_size(); }
    size_type dis_ncol () const              { return col_ownership().dis_size(); }
    size_type dis_nnz () const               { return base::data().nnz(); }
    size_type dis_ext_nnz () const           { return 0; }
    bool is_symmetric() const                { return base::data().is_symmetric(); }
    void set_symmetry (bool is_symm) const   { base::data().set_symmetry(is_symm); }
    void set_symmetry_by_check (const T& tol = std::numeric_limits<T>::epsilon()) const
                                             { base::data().set_symmetry_by_check(); }
    size_type pattern_dimension() const      { return base::data().pattern_dimension(); }
    void set_pattern_dimension(size_type dim) const { base::data().set_pattern_dimension(dim); }
    T max_abs () const                       { return base::data().max_abs(); }

    // local sizes
    size_type nrow () const		     { return base::data().nrow(); }
    size_type ncol () const		     { return base::data().ncol(); }
    size_type nnz () const		     { return base::data().nnz(); }
    
    // range on local memory
    size_type row_first_index () const       { return base::data().row_first_index(); }
    size_type row_last_index () const        { return base::data().row_last_index(); }
    size_type col_first_index () const       { return base::data().col_first_index(); }
    size_type col_last_index () const        { return base::data().col_last_index(); }

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

    // accessors, only for distributed (for interface compatibility)
    size_type ext_nnz() const                { return 0; }
    const_iterator ext_begin() const         { return const_iterator(); }
    const_iterator ext_end()   const         { return const_iterator(); }
          iterator ext_begin_nonconst()      { return iterator(); }
          iterator ext_end_nonconst()        { return iterator(); }
    size_type jext2dis_j (size_type jext) const { return 0; }

// algebra:

    // y := a*x
    void mult (const vec<element_type,sequential>& x, vec<element_type,sequential>& y) const {
      base::data().mult (x,y);
    }
    vec<element_type,sequential> operator* (const vec<element_type,sequential>& x) const {
      vec<element_type,sequential> y (row_ownership(), element_type());
      mult (x, y);
      return y;
    }
    void trans_mult (const vec<element_type,sequential>& x, vec<element_type,sequential>& y) const {
      base::data().trans_mult (x,y);
    }
    vec<element_type,sequential> trans_mult (const vec<element_type,sequential>& x) const {
      vec<element_type,sequential> y (col_ownership(), element_type());
      trans_mult (x, y);
      return y;
    }
    // a+b, a-b, a*b
    csr<T,sequential> operator+ (const csr<T,sequential>& b) const;
    csr<T,sequential> operator- (const csr<T,sequential>& b) const;
    csr<T,sequential> operator* (const csr<T,sequential>& b) const;

    // lambda*a
    csr<T,sequential>& operator*= (const T& lambda) {
      base::data().operator*= (lambda);
      return *this;
    }
// output:

    void dump (const std::string& name) const { base::data().dump(name); }
};
// lambda*a
template<class T>
inline
csr<T,sequential>
operator* (const T& lambda, const csr<T,sequential>& a)
{
  csr<T,sequential> b = a;
  b.operator*= (lambda);
  return b;
}
// -a
template<class T>
inline
csr<T,sequential>
operator- (const csr<T,sequential>& a)
{
  return T(-1)*a;
}
// trans(a)
template<class T>
inline
csr<T,sequential>
trans (const csr<T,sequential>& a)
{
  csr<T,sequential> b;
  a.data().build_transpose (b.data());
  return b;
}
//>verbatim:

#ifdef _RHEOLEF_HAVE_MPI
//<verbatim:
template<class T>
class csr<T,distributed> : public smart_pointer<csr_rep<T,distributed> > {
public:

// typedefs:

    typedef csr_rep<T,distributed>                    rep;
    typedef smart_pointer<rep>                base;
    typedef typename rep::memory_type         memory_type;
    typedef typename rep::size_type           size_type;
    typedef typename rep::element_type        element_type;
    typedef typename rep::iterator            iterator;
    typedef typename rep::const_iterator      const_iterator;
    typedef typename rep::data_iterator       data_iterator;
    typedef typename rep::const_data_iterator const_data_iterator;

// allocators/deallocators:

    csr() : base(new_macro(rep())) {}
    template<class A>
    explicit csr(const asr<T,memory_type,A>& a) : base(new_macro(rep(a))) {}
    void resize (const distributor& row_ownership, const distributor& col_ownership, size_type nnz1 = 0)
        { base::data().resize(row_ownership, col_ownership, nnz1); }

// allocators from initializer list (c++ 2011):

#ifdef _RHEOLEF_HAVE_STD_INITIALIZER_LIST
    csr (const std::initializer_list<csr_concat_value<T,distributed> >& init_list);
    csr (const std::initializer_list<csr_concat_line<T,distributed> >&  init_list);
#endif // _RHEOLEF_HAVE_STD_INITIALIZER_LIST

// accessors:

    // global sizes
    const distributor& row_ownership() const { return base::data().row_ownership(); }
    const distributor& col_ownership() const { return base::data().col_ownership(); }
    size_type dis_nrow () const              { return row_ownership().dis_size(); }
    size_type dis_ncol () const              { return col_ownership().dis_size(); }
    size_type dis_nnz () const               { return base::data().dis_nnz(); }
    size_type dis_ext_nnz () const           { return base::data().dis_ext_nnz(); }
    bool is_symmetric() const                { return base::data().is_symmetric(); }
    void set_symmetry (bool is_symm) const   { base::data().set_symmetry(is_symm); }
    void set_symmetry_by_check (const T& tol = std::numeric_limits<T>::epsilon()) const
                                             { base::data().set_symmetry_by_check(); }
    size_type pattern_dimension() const      { return base::data().pattern_dimension(); }
    void set_pattern_dimension(size_type dim) const { base::data().set_pattern_dimension(dim); }
    T max_abs () const                       { return base::data().max_abs(); }

    // local sizes
    size_type nrow () const		     { return base::data().nrow(); }
    size_type ncol () const		     { return base::data().ncol(); }
    size_type nnz () const		     { return base::data().nnz(); }
    
    // range on local memory
    size_type row_first_index () const       { return base::data().row_first_index(); }
    size_type row_last_index () const        { return base::data().row_last_index(); }
    size_type col_first_index () const       { return base::data().col_first_index(); }
    size_type col_last_index () const        { return base::data().col_last_index(); }

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

    // accessors, only for distributed
    size_type ext_nnz() const                { return base::data().ext_nnz(); }
    const_iterator ext_begin() const         { return base::data().ext_begin(); }
    const_iterator ext_end()   const         { return base::data().ext_end(); }
          iterator ext_begin_nonconst()      { return base::data().ext_begin(); }
          iterator ext_end_nonconst()        { return base::data().ext_end(); }
    size_type jext2dis_j (size_type jext) const { return base::data().jext2dis_j(jext); }

// algebra:


    // y := a*x
    void mult (const vec<element_type,distributed>& x, vec<element_type,distributed>& y) const {
      base::data().mult (x,y);
    }
    vec<element_type,distributed> operator* (const vec<element_type,distributed>& x) const {
      vec<element_type,distributed> y (row_ownership(), element_type());
      mult (x, y);
      return y;
    }
    void trans_mult (const vec<element_type,distributed>& x, vec<element_type,distributed>& y) const {
      base::data().trans_mult (x,y);
    }
    vec<element_type,distributed> trans_mult (const vec<element_type,distributed>& x) const {
      vec<element_type,distributed> y (col_ownership(), element_type());
      trans_mult (x, y);
      return y;
    }
    // a+b, a-b, a*b
    csr<T,distributed> operator+ (const csr<T,distributed>& b) const;
    csr<T,distributed> operator- (const csr<T,distributed>& b) const;
    csr<T,distributed> operator* (const csr<T,distributed>& b) const;

    // lambda*a
    csr<T,distributed>& operator*= (const T& lambda) {
      base::data().operator*= (lambda);
      return *this;
    }
// output:

    void dump (const std::string& name) const { base::data().dump(name); }
};
// lambda*a
template<class T>
inline
csr<T,distributed>
operator* (const T& lambda, const csr<T,distributed>& a)
{
  csr<T,distributed> b = a;
  b.operator*= (lambda);
  return b;
}
// -a
template<class T>
inline
csr<T,distributed>
operator- (const csr<T,distributed>& a)
{
  return T(-1)*a;
}
// trans(a)
template<class T>
inline
csr<T,distributed>
trans (const csr<T,distributed>& a)
{
  csr<T,distributed> b;
  a.data().build_transpose (b.data());
  return b;
}
#endif // _RHEOLEF_HAVE_MPI

// b = f(a); f as a class-function or usual fct
template<class T, class M, class Function>
csr<T,M>
apply (Function f, const csr<T,M>& a)
{
  csr<T,M> b = a;
  typename csr<T,M>::size_type n = a.nrow(); 
  typename csr<T,M>::const_iterator dia_ia = a.begin(); 
  typename csr<T,M>::iterator       dia_ib = b.begin_nonconst(); 
  pair_transform_second (dia_ia[0], dia_ia[n], dia_ib[0], f);
  if (a.ext_nnz() != 0) {
    typename csr<T,M>::const_iterator ext_ia = a.ext_begin(); 
    typename csr<T,M>::iterator       ext_ib = b.ext_begin_nonconst(); 
    pair_transform_second (ext_ia[0], ext_ia[n], ext_ib[0], f);
  }
  return b;
}
template<class T, class M, class Function>
csr<T,M>
apply (T (*f)(const T&), const csr<T,M>& a)
{
  return apply (std::ptr_fun(f), a);
}
//>verbatim:

template<class T, class M>
csr<T,M>
diag (const vec<T,M>&);

// ------------------------------
// i/o
// ------------------------------
template <class T, class M>
inline
idiststream&
operator >> (idiststream& s,  csr<T,M>& x)
{ 
    return x.data().get(s); 
}
template <class T, class M> 
inline
odiststream&
operator << (odiststream& s, const csr<T,M>& x)
{
    return x.data().put(s); 
}
// ------------------------------
// a+b, a-b
// ------------------------------
template <class T>
inline
csr<T,sequential>
csr<T,sequential>::operator+ (const csr<T,sequential>& b) const {
  csr<T,sequential> c;
  c.data().assign_add (this->data(), b.data(), std::plus<T>());
  return c;
}
template <class T> 
inline
csr<T,sequential>
csr<T,sequential>::operator- (const csr<T,sequential>& b) const {
  csr<T,sequential> c;
  c.data().assign_add (this->data(), b.data(), std::minus<T>());
  return c;
}
#ifdef _RHEOLEF_HAVE_MPI
template <class T> 
inline
csr<T,distributed>
csr<T,distributed>::operator+ (const csr<T,distributed>& b) const {
  csr<T,distributed> c;
  c.data().assign_add (this->data(), b.data(), std::plus<T>());
  return c;
}
template <class T> 
inline
csr<T,distributed>
csr<T,distributed>::operator- (const csr<T,distributed>& b) const {
  csr<T,distributed> c;
  c.data().assign_add (this->data(), b.data(), std::minus<T>());
  return c;
}
#endif // _RHEOLEF_HAVE_MPI
  
// ------------------------------
// a*b
// ------------------------------
template <class T>
inline
csr<T,sequential>
csr<T,sequential>::operator* (const csr<T,sequential>& b) const {
  csr<T,sequential> c;
  c.data().assign_mult (this->data(), b.data());
  return c;
}
#ifdef _RHEOLEF_HAVE_MPI
template <class T> 
inline
csr<T,distributed>
csr<T,distributed>::operator* (const csr<T,distributed>& b) const {
  csr<T,distributed> c;
  c.data().assign_mult (this->data(), b.data());
  return c;
}
#endif // _RHEOLEF_HAVE_MPI

// ------------------------------
// a.max_abs
// ------------------------------
template <class T>
inline
T
csr_rep<T,sequential>::max_abs () const
{
  T val = 0;
  typename csr<T,sequential>::const_iterator ia = begin();
  for (typename csr_rep<T,sequential>::const_data_iterator iter = ia[0], last = ia[nrow()]; iter != last; ++iter) {
    val = std::max (val, abs((*iter).second));
  }
  return val;
}
#ifdef _RHEOLEF_HAVE_MPI
template <class T>
inline
T
csr_rep<T,distributed>::max_abs () const
{
  T val = csr_rep<T,sequential>::max_abs();
  val = mpi::all_reduce (comm(), val, mpi::maximum<T>());
  return val;
}
#endif // _RHEOLEF_HAVE_MPI
  
} // namespace rheolef
# endif // _RHEOLEF_CSR_H