Loading...
Searching...
No Matches
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
_RHEOLEF_instanciation
#define _RHEOLEF_instanciation(T, M, A)
Definition asr.cc:223
Float
see the Float page for the full documentation
rheolef::asr
Definition asr.h:49
rheolef::asr::put_seq_sparse_matlab
odiststream & put_seq_sparse_matlab(odiststream &ops, size_type first_dis_i=0) const
Definition asr.cc:85
rheolef::asr::get
idiststream & get(idiststream &ips)
Definition asr.cc:183
rheolef::asr::row_type
pair_set< T, A > row_type
Definition asr.h:53
rheolef::asr::put_mpi
odiststream & put_mpi(odiststream &ops) const
Definition asr.cc:144
rheolef::asr::put
odiststream & put(odiststream &ops) const
Definition asr.cc:173
rheolef::asr::size_type
base::size_type size_type
Definition asr.h:55
rheolef::asr::build_from_csr
void build_from_csr(const csr_rep< T, M > &)
Definition asr.cc:37
rheolef::asr::_recompute_nnz
void _recompute_nnz()
Definition asr.cc:65
rheolef::asr::put_seq
odiststream & put_seq(odiststream &ops, size_type first_dis_i=0) const
Definition asr.cc:134
rheolef::asr::put_seq_matrix_market
odiststream & put_seq_matrix_market(odiststream &ops, size_type first_dis_i=0) const
Definition asr.cc:110
rheolef::csr_rep
Definition csr.h:81
rheolef::csr
see the csr page for the full documentation
Definition csr.h:317
rheolef::distributor
see the distributor page for the full documentation
Definition distributor.h:69
rheolef::distributor::decide
static const size_type decide
Definition distributor.h:83
rheolef::idiststream
idiststream: see the diststream page for the full documentation
Definition diststream.h:336
rheolef::idiststream::io_proc
static size_type io_proc()
This routine returns the rank of a process that can perform i/o.
Definition diststream.cc:64
rheolef::idiststream::is
std::istream & is()
Definition diststream.h:400
rheolef::idiststream::comm
const communicator & comm() const
Definition diststream.h:356
rheolef::iorheo::flag_type
std::bitset< last > flag_type
Definition iorheo.h:440
rheolef::iorheo::format_field
static flag_type format_field
Definition iorheo.h:445
rheolef::iorheo::matlab
@ matlab
Definition iorheo.h:386
rheolef::iorheo::sparse_matlab
@ sparse_matlab
Definition iorheo.h:387
rheolef::odiststream
odiststream: see the diststream page for the full documentation
Definition diststream.h:137
rheolef::odiststream::os
std::ostream & os()
Definition diststream.h:247
rheolef::odiststream::io_proc
static size_type io_proc()
Definition diststream.cc:79
T
Expr1::float_type T
Definition field_expr.h:230
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
rheolef
This file is part of Rheolef.
Definition compiler_eigen.h:39
rheolef::read_matrix_market_header
struct matrix_market read_matrix_market_header(idiststream &ips)
Definition mm_io.cc:30
p
Definition sphere.icc:25
rheolef::distributed
Definition communicator.h:128
rheolef::is_distributed
Definition communicator.h:70
rheolef::is_sequential
Definition communicator.h:66
rheolef::matrix_market
Definition mm_io.h:29
rheolef::matrix_market::general
static const format_type general
Definition mm_io.h:32
rheolef::matrix_market::symmetric
static const format_type symmetric
Definition mm_io.h:33
rheolef::matrix_market::skew_symmetric
static const format_type skew_symmetric
Definition mm_io.h:34
rheolef::matrix_market::hermitian
static const format_type hermitian
Definition mm_io.h:35
rheolef::matrix_market::format
format_type format
Definition mm_io.h:37
rheolef::sequential
Definition communicator.h:62