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

# ifndef _RHEO_NEW_ASR_H
# define _RHEO_NEW_ASR_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/array.h"
#include "rheolef/pair_set.h"
#include "rheolef/diststream.h"
// -------------------------------------------------------------
// the asr class
// -------------------------------------------------------------
namespace rheolef {

template<class T, class M> class csr;
template<class T, class M> class csr_rep;

/*Class:asr
NAME:  @code{asr} - associative sparse matrix (@PACKAGE@-@VERSION@)
SYNOPSYS:       
  Associative sparse matrix container, used during FEM assembling process.
IMPLEMENTATION NOTE:
  Elements are stored row by row using the pair_set class bqsed on the STL map class.
  Implementation of asr uses array<pair_set>
AUTHORS:
    LMC-IMAG, 38041 Grenoble cedex 9, France
   | Pierre.Saramito@imag.fr
DATE:   6 january 1999, last update 20 may 2012.
End:
*/
//<verbatim:
template<class T, class M = rheo_default_memory_model, class A = std::allocator<T> >
class asr : public array<pair_set<T,A>, M, A> {
public:
// typedefs:

    typedef pair_set<T,A>          	     row_type;
    typedef array<row_type,M,A>              base;
    typedef typename base::size_type         size_type;
    typedef M                                memory_type;

    struct dis_reference {
      dis_reference (typename base::dis_reference row_dis_i, size_type dis_j)
       : _row_dis_i(row_dis_i), _dis_j(dis_j) {}

      dis_reference& operator+= (const T& value) {
        _row_dis_i += std::pair<size_type,T>(_dis_j,value);
        return *this;
      }
      typename base::dis_reference  _row_dis_i;
      size_type                     _dis_j;
    };

// allocators/deallocators:

    asr (const A& alloc = A())
      : base(distributor(), row_type(alloc), alloc), _col_ownership(), _nnz(0), _dis_nnz(0) {}
 
    asr (const distributor& row_ownership, const distributor& col_ownership, const A& alloc = A())
      : base(row_ownership, row_type(alloc), alloc), _col_ownership(col_ownership), _nnz(0), _dis_nnz(0) {}
 
    asr (const csr_rep<T,M>&, const A& alloc = A());
    asr (const csr<T,M>&,     const A& alloc = A());
    void build_from_csr (const csr_rep<T,M>&);
 
    void resize (const distributor& row_ownership, const distributor& col_ownership) 
    {
       base::resize (row_ownership);
       _col_ownership = col_ownership;
       _nnz = _dis_nnz = 0;
    }

// accessors:

    const communicator& comm() const { return base::comm(); }

    size_type nrow () const { return base::size(); }
    size_type ncol () const { return _col_ownership.size(); }
    size_type nnz () const { return _nnz; }
    
    size_type dis_nrow () const { return base::dis_size(); }
    size_type dis_ncol () const { return _col_ownership.dis_size(); }
    size_type dis_nnz () const  { return _dis_nnz; }
    const distributor& row_ownership() const { return base::ownership(); }
    const distributor& col_ownership() const { return _col_ownership; }

// modifiers:

    T operator()            (size_type     i, size_type dis_j) const;
    T&       semi_dis_entry (size_type     i, size_type dis_j);
    dis_reference dis_entry (size_type dis_i, size_type dis_j);

    // dis_entry_assembly_end is redefined in order to recompute _nnz and _dis_nnz
    void dis_entry_assembly_begin() { base::dis_entry_assembly_begin (pair_set_add_op<row_type>()); }
    void dis_entry_assembly_end()   { base::dis_entry_assembly_end   (pair_set_add_op<row_type>()); _recompute_nnz(); }
    void dis_entry_assembly()       { dis_entry_assembly_begin(); dis_entry_assembly_end(); }

// io:
    odiststream& put (odiststream& ops) const;
    idiststream& get (idiststream& ips);

// internal:
    odiststream& put_mpi               (odiststream& ops) const;
    odiststream& put_seq               (odiststream& ops, size_type first_dis_i = 0) const;
    odiststream& put_seq_sparse_matlab (odiststream& ops, size_type first_dis_i = 0) const;
    odiststream& put_seq_matrix_market (odiststream& ops, size_type first_dis_i = 0) const;
protected:
    void _recompute_nnz();
// data:
    distributor _col_ownership;
    size_type   _nnz;
    size_type   _dis_nnz;
};
//>verbatim:
// ----------------------------------------------------------------------------
// inlined
// ----------------------------------------------------------------------------
template <class T, class M, class A>
inline
asr<T,M,A>::asr (const csr<T,M>& a, const A& alloc)
 : base(a.row_ownership(), row_type(alloc), alloc),
   _col_ownership(a.col_ownership()),
   _nnz(a.nnz()),
   _dis_nnz(a.dis_nnz())
{
  build_from_csr (a.data());
}
template <class T, class M, class A>
inline
asr<T,M,A>::asr (const csr_rep<T,M>& a, const A& alloc)
 : base(a.row_ownership(), row_type(alloc), alloc),
   _col_ownership(a.col_ownership()),
   _nnz(a.nnz()),
   _dis_nnz(a.dis_nnz())
{
  build_from_csr (a);
}
template <class T, class M, class A>
inline
idiststream&
operator>> (idiststream& s,  asr<T,M,A>& x)
{
    return x.get(s);
}
template <class T, class M, class A>
inline
odiststream&
operator<< (odiststream& s, const asr<T,M,A>& x)
{
    return x.put(s);
}
template <class T, class M, class A>
inline
T
asr<T,M,A>::operator() (size_type i, size_type dis_j) const
{
  typename row_type::const_iterator pos_aij  = base::operator[](i).find(dis_j);
  if (pos_aij != base::operator[](i).end()) {
    return (*pos_aij).second;
  } else {
    return T(0);
  }
}
template <class T, class M, class A>
inline
T&
asr<T,M,A>::semi_dis_entry (size_type i, size_type dis_j)
{
   row_type& row_i = base::operator[](i);
   std::pair<typename row_type::iterator,bool> status
    = row_i.insert (std::pair<size_type,T>(dis_j,T(0)));
   return (*(status.first)).second;
}
template <class T, class M, class A>
inline
typename asr<T,M,A>::dis_reference
asr<T,M,A>::dis_entry (size_type dis_i, size_type dis_j)
{
  assert_macro (dis_i < dis_nrow() && dis_j < dis_ncol(), 
	  "indexes ("<<dis_i<<" "<<dis_j<<") out of range [0:"
	  << dis_nrow() << "[x[0:" << dis_ncol() << "[");
  return dis_reference (base::dis_entry (dis_i), dis_j);
}

} // namespace rheolef
# endif // _RHEO_NEW_ASR_H