asr.cc

Go to the documentation of this file.

#ifdef USE_NEW_ASR
#define asr asr
#endif // USE_NEW_ASR
 
#include "rheolef/asr.h"
#include "rheolef/csr.h"
#include "rheolef/iorheo.h"
#include "rheolef/mm_io.h"
 
namespace rheolef {
 
// ----------------------------------------------------------------------------
// csr2asr convertion
// ----------------------------------------------------------------------------
template<class T, class M, class A>
void
asr<T,M,A>::build_from_csr (const csr_rep<T,M>& a)
{
  typename csr_rep<T,M>::const_iterator
        dia_ia = a.begin(),
        ext_ia = a.ext_begin();
  size_type first_dis_i = a.row_ownership().first_index();
  size_type first_dis_j = a.col_ownership().first_index();
  for (size_type i = 0, n = a.nrow(); i < n; i++) {
    for (typename csr<T,M>::const_data_iterator p = dia_ia[i]; p < dia_ia[i+1]; ++p) {
      size_type dis_j = first_dis_j + (*p).first;
      const T&  value = (*p).second; 
      semi_dis_entry (i, dis_j) += value;
    }
    if (is_sequential<M>::value || a.ext_nnz() == 0) continue;
    size_type dis_i = first_dis_i + i;
    for (typename csr<T,M>::const_data_iterator p = ext_ia[i]; p < ext_ia[i+1]; p++) {
      size_type dis_j = a.jext2dis_j ((*p).first);
      const T&  value = (*p).second; 
      semi_dis_entry (i, dis_j) += value;
    }
  }
  dis_entry_assembly();
}
// ----------------------------------------------------------------------------
// after a asr::dis_entry_assembly, may recompute nnz and dis_nnz
// ----------------------------------------------------------------------------
template <class T, class M, class A>
void
asr<T,M,A>::_recompute_nnz ()
{
  _nnz = 0;
  for (typename base::const_iterator
          iter = base::begin(),
          last = base::end(); iter != last; ++iter) {
    _nnz += (*iter).size();
  }
  _dis_nnz = _nnz;
#ifdef _RHEOLEF_HAVE_MPI
  if (is_distributed<M>::value) {
    _dis_nnz = mpi::all_reduce (comm(), _nnz, std::plus<size_type>());
  }
#endif // _RHEOLEF_HAVE_MPI
}
// ----------------------------------------------------------------------------
// io
// ----------------------------------------------------------------------------
template <class T, class M, class A>
odiststream&
asr<T,M,A>::put_seq_sparse_matlab (odiststream& ods, size_type first_dis_i) const
{
  std::ostream& os = ods.os();
  std::string name = iorheo::getbasename(ods.os());
  if (name == "") name = "a";
  os << std::setprecision (std::numeric_limits<T>::digits10)
     << name << " = spalloc(" << nrow() << "," << ncol() << "," << nnz() << ");" << std::endl;
  if (_nnz == 0) return ods;
  size_type i = 0;
  for (typename base::const_iterator
          iter = base::begin(),
          last = base::end(); iter != last; ++iter, ++i) {
    const row_type& row = *iter;
    for (typename row_type::const_iterator
            row_iter = row.begin(),
            row_last = row.end(); row_iter != row_last; ++row_iter) {
      os << name << "(" << first_dis_i+i+1 << ","
         << (*row_iter).first+1 << ") = "
         << (*row_iter).second << ";" << std::endl;
    }
  }
  return ods;
}
template <class T, class M, class A>
odiststream&
asr<T,M,A>::put_seq_matrix_market (odiststream& ods, size_type first_dis_i) const
{
  std::ostream& os = ods.os();
  os << std::setprecision(std::numeric_limits<T>::digits10)
     << "%%MatrixMarket matrix coordinate real general" << std::endl
     << nrow() << " " << ncol() << " " << nnz() << std::endl;
  if (_nnz == 0) return ods;
  size_type i = 0;
  for (typename base::const_iterator
          iter = base::begin(),
          last = base::end(); iter != last; ++iter, ++i) {
    const row_type& row = *iter;
    for (typename row_type::const_iterator
            row_iter = row.begin(),
            row_last = row.end(); row_iter != row_last; ++row_iter) {
      os << first_dis_i+i+1 << " "
         << (*row_iter).first+1 << " "
         << (*row_iter).second << std::endl;
    }
  }
  return ods;
}
template <class T, class M, class A>
odiststream&
asr<T,M,A>::put_seq (odiststream& ods, size_type first_dis_i) const
{
  iorheo::flag_type format = iorheo::flags(ods.os()) & iorheo::format_field;
  if (format [iorheo::matlab] || format [iorheo::sparse_matlab])
                                      { return put_seq_sparse_matlab   (ods,first_dis_i); }
  // default is matrix market format
  return put_seq_matrix_market (ods,first_dis_i);
}
template <class T, class M, class A>
odiststream&
asr<T,M,A>::put_mpi (odiststream& ods) const
{
  // send all in a pseudo-sequential matrix on process 0
  size_type io_proc = odiststream::io_proc();
  size_type my_proc = comm().rank();
  distributor io_row_ownership (dis_nrow(), comm(), (my_proc == io_proc ? dis_nrow() : 0));
  distributor io_col_ownership (dis_ncol(), comm(), (my_proc == io_proc ? dis_ncol() : 0));
  asr<T> a (io_row_ownership, io_col_ownership);
  size_type first_dis_i = row_ownership().first_index();
  size_type i = 0;
  for (typename base::const_iterator
          iter = base::begin(),
          last = base::end(); iter != last; ++iter, ++i) {
    const row_type& row = *iter;
    for (typename row_type::const_iterator
            row_iter = row.begin(),
            row_last = row.end(); row_iter != row_last; ++row_iter) {
      a.dis_entry (first_dis_i + i, (*row_iter).first) += (*row_iter).second;
    }
  }
  a.dis_entry_assembly();
  if (my_proc == io_proc) {
    // print sequential matrix on io_proc
    a.put_seq (ods);
  }
  return ods;
}
template <class T, class M, class A>
odiststream&
asr<T,M,A>::put (odiststream& ods) const
{
  if (is_sequential<M>::value) {
    return put_seq (ods);
  } else {
    return put_mpi (ods);
  }
}
template <class T, class M, class A>
idiststream&
asr<T,M,A>::get (idiststream& ips) 
{
  // matrix market format:
  //   %%mm_header
  //   nrow ncol nnz
  //   {i, j, aij}*
  // we suppose nrow ncol already readed
  // and the matrix has good dimensions
  struct matrix_market mm = read_matrix_market_header (ips);
  check_macro (mm.format != matrix_market::hermitian, "Hermitian matrix not yet supported");
  size_type dis_nrow, dis_ncol, dis_nnz;
  ips >> dis_nrow >> dis_ncol >> dis_nnz;
  communicator comm = ips.comm();
  distributor row_ownership (dis_nrow, comm, distributor::decide);
  distributor col_ownership (dis_ncol, comm, distributor::decide);
  resize (row_ownership, col_ownership);
  size_type my_proc = comm.rank();
  if (my_proc == idiststream::io_proc()) {
    std::istream& is = ips.is();
    for (size_type p = 0; p < dis_nnz; p++) {
      size_type dis_i, dis_j;
      T value;
      is >> dis_i >> dis_j >> value;
      dis_i--; dis_j--; // 1..N convention -> 0..N-1
      dis_entry(dis_i,dis_j) += value;
      if (mm.format == matrix_market::general || dis_i == dis_j) continue;
      // when symmetries, add also the corespondings coefs, for simplicity
      if (mm.format == matrix_market::symmetric) {
         dis_entry(dis_j,dis_i) += value;
      } else if (mm.format == matrix_market::skew_symmetric) {
         dis_entry(dis_j,dis_i) += -value;
      }
    }
  }
  dis_entry_assembly();
  return ips;
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciation(T,M,A)                                               \
template class asr<T,M,A>;
 
_RHEOLEF_instanciation(Float,sequential,std::allocator<Float>)
#ifdef _RHEOLEF_HAVE_MPI
_RHEOLEF_instanciation(Float,distributed,std::allocator<Float>)
#endif // _RHEOLEF_HAVE_MPI
 
} // namespace rheolef