Loading...
Searching...
No Matches
csr_mpi.cc
Go to the documentation of this file.
# include "rheolef/config.h"
#ifdef _RHEOLEF_HAVE_MPI
#include "rheolef/dis_macros.h"
#include "rheolef/csr.h"
#include "rheolef/asr_to_csr.h"
#include "rheolef/asr_to_csr_dist_logical.h"
#include "rheolef/csr_amux.h"
#include "rheolef/csr_cumul_trans_mult.h"
#include "rheolef/mpi_scatter_init.h"
#include "rheolef/mpi_scatter_begin.h"
#include "rheolef/mpi_scatter_end.h"
using namespace std;
namespace rheolef {
// ----------------------------------------------------------------------------
// useful class-functions
// ----------------------------------------------------------------------------
#include <algorithm> // for lower_bound
template <class Pair1, class Pair2, class RandomIterator>
struct op_ext2glob_t : unary_function<Pair1,Pair2> {
Pair2 operator()(const Pair1& x) const {
return Pair2(t[x.first], x.second); }
op_ext2glob_t(RandomIterator t1) : t(t1) {}
RandomIterator t;
};
template <class Pair1, class Pair2>
struct op_dia_t : unary_function<Pair1,Pair2> {
Pair2 operator()(const Pair1& x) const {
return Pair2(shift + x.first, x.second); }
typedef typename Pair1::first_type Size;
op_dia_t(Size s) : shift(s) {}
Size shift;
};
template <class Pair1, class Pair2, class RandomIterator>
struct op_dis_j2jext_t : unary_function<Pair1,Pair2> {
Pair2 operator()(const Pair1& x) const {
RandomIterator t = std::lower_bound(t1, t2, x.first);
assert_macro(*t == x.first, "problem in ditributed asr_to_csr");
return Pair2(distance(t1,t), x.second); }
op_dis_j2jext_t(RandomIterator u1, RandomIterator u2) : t1(u1), t2(u2) {}
RandomIterator t1, t2;
};
// ----------------------------------------------------------------------------
// allocators
// ----------------------------------------------------------------------------
template<class T>
csr_rep<T,distributed>::csr_rep()
: base(),
_ext(),
_jext2dis_j(0),
_dis_nnz(0),
_dis_ext_nnz(0),
_scatter_initialized(false),
_from(),
_to(),
_buffer()
{
}
template<class T>
csr_rep<T,distributed>::csr_rep(const csr_rep<T,distributed>& a)
: base(a),
_ext(a._ext),
_jext2dis_j(a._jext2dis_j),
_dis_nnz(a._dis_nnz),
_dis_ext_nnz(a._dis_ext_nnz),
_scatter_initialized(a._scatter_initialized),
_from(a._from),
_to(a._to),
_buffer(a._buffer)
{
// "physical copy of csr"
}
template<class T>
void
csr_rep<T,distributed>::resize (const distributor& row_ownership, const distributor& col_ownership, size_type nnz1)
{
base::resize (row_ownership, col_ownership, nnz1);
_ext.resize (row_ownership, col_ownership, 0); // note: the _ext part will be resized elsewhere
_scatter_initialized = false;
}
template<class T>
template<class A>
void
csr_rep<T,distributed>::build_from_asr (const asr<T,distributed,A>& a)
{
base::resize (a.row_ownership(), a.col_ownership(), 0);
_ext.resize (a.row_ownership(), a.col_ownership(), 0);
distributor::tag_type tag = distributor::get_new_tag();
typedef typename asr<T>::row_type row_type;
typedef typename row_type::value_type const_pair_type;
typedef pair<size_type,T> pair_type;
is_dia_t<size_type, const_pair_type>
is_dia(col_ownership().first_index(), col_ownership().last_index());
//
// step 1: compute pointers
//
set<size_type> colext;
size_type nnzext
= asr_to_csr_dist_logical (a.begin(), a.end(), is_dia, colext);
//
// step 2: resize and copy jext2dis_j
//
size_type nnzdia = a.nnz() - nnzext;
size_type ncoldia = col_ownership().size();
size_type ncolext = colext.size();
base::resize(a.row_ownership(), a.col_ownership(), nnzdia);
_ext.resize (a.row_ownership(), a.col_ownership(), nnzext);
_jext2dis_j.resize (ncolext);
copy (colext.begin(), colext.end(), _jext2dis_j.begin());
//
// step 3: copy values
// column indexes of a in 0..dis_ncol
// of dia in 0..ncoldia
// of dia in 0..ncolext
op_dia_t<const_pair_type,pair_type> op_dia(- col_ownership().first_index());
// step 3.a: copy dia part
asr_to_csr (
a.begin(),
a.end(),
is_dia,
op_dia,
base::begin(),
base::_data.begin());
// step 3.b: copy ext part
unary_negate<is_dia_t<size_type, const_pair_type> >
is_ext(is_dia);
op_dis_j2jext_t<const_pair_type, pair_type, typename vector<size_type>::const_iterator>
op_dis_j2jext(_jext2dis_j.begin(), _jext2dis_j.end());
asr_to_csr (
a.begin(), a.end(), is_ext, op_dis_j2jext,
_ext.begin(), _ext._data.begin());
// compute _dis_nnz:
_dis_nnz = a.dis_nnz();
_dis_ext_nnz = mpi::all_reduce (comm(), _ext.nnz(), std::plus<size_type>());
_scatter_initialized = false;
#ifdef TO_CLEAN
_dis_nnz = mpi::all_reduce (comm(), nnz() + _ext.nnz(), std::plus<size_type>());
check_macro (_dis_nnz == a.dis_nnz(), "build_from_asr: mistake: asr::dis_nnz="<<a.dis_nnz()
<< " while _dis_nnz="<<_dis_nnz);
#endif // TO_CLEAN
}
// messages building for A*x and trans(A)*x
template<class T>
void
csr_rep<T,distributed>::_scatter_init() const
{
vector<size_type> id (_jext2dis_j.size());
for (size_type i = 0, n = id.size(); i < n; i++) id[i] = i;
distributor::tag_type tag = distributor::get_new_tag();
mpi_scatter_init(
_jext2dis_j.size(),
_jext2dis_j.begin().operator->(),
id.size(),
id.begin().operator->(),
dis_ncol(),
col_ownership().begin().operator->(),
tag,
row_ownership().comm(),
_from,
_to);
_buffer.resize(_jext2dis_j.size());
}
// ----------------------------------------------------------------------------
// io
// ----------------------------------------------------------------------------
template<class T>
odiststream&
csr_rep<T,distributed>::put (odiststream& ops) const
{
// put all on io_proc
size_type io_proc = odiststream::io_proc();
size_type my_proc = comm().rank();
distributor a_row_ownership (dis_nrow(), comm(), (my_proc == io_proc ? dis_nrow() : 0));
distributor a_col_ownership (dis_ncol(), comm(), (my_proc == io_proc ? dis_ncol() : 0));
typedef std::allocator<T> A; // TODO: use heap_alloc for asr
asr<T,distributed,A> a (a_row_ownership, a_col_ownership);
size_type first_dis_i = row_ownership().first_index();
size_type first_dis_j = col_ownership().first_index();
if (nnz() != 0) {
const_iterator ia = begin();
for (size_type i = 0; i < nrow(); i++) {
for (const_data_iterator p = ia[i]; p < ia[i+1]; p++) {
const size_type& j = (*p).first;
const T& val = (*p).second;
a.dis_entry (i+first_dis_i, j+first_dis_j) += val;
}
}
}
if (_ext.nnz() != 0) {
const_iterator ext_ia = _ext.begin();
for (size_type i = 0; i < nrow(); i++) {
for (const_data_iterator p = ext_ia[i]; p < ext_ia[i+1]; p++) {
const size_type& j = (*p).first;
const T& val = (*p).second;
a.dis_entry (i+first_dis_i, _jext2dis_j[j]) += val;
}
}
}
a.dis_entry_assembly();
if (my_proc == io_proc) {
a.put_seq (ops);
}
return ops;
}
template<class T>
void
csr_rep<T,distributed>::dump (const string& name) const
{
odiststream ops;
std::string filename = name + std::to_string(comm().rank());
ops.open (filename, "mtx", comm());
check_macro(ops.good(), "\"" << filename << "[.mtx]\" cannot be created.");
ops << "%%MatrixMarket matrix coordinate real general" << std::endl
<< dis_nrow() << " " << dis_ncol() << " " << dis_nnz() << std::endl;
put(ops);
}
// ----------------------------------------------------------------------------
// blas2: basic linear algebra
// ----------------------------------------------------------------------------
template<class T>
void
csr_rep<T,distributed>::mult(
const vec<T,distributed>& x,
vec<T,distributed>& y)
const
{
check_macro (x.size() == ncol(), "csr*vec: incompatible csr("<<nrow()<<","<<ncol()<<") and vec("<<x.size()<<")");
y.resize (row_ownership());
distributor::tag_type tag = distributor::get_new_tag();
// initialize _from and _to scatter messages
_scatter_init_guard();
// send x to others
mpi_scatter_begin (
x.begin().operator->(),
_buffer.begin().operator->(),
_from,
_to,
details::generic_set_op(),
tag,
row_ownership().comm());
// y := dia*x
csr_amux (
base::begin(),
base::end(),
x.begin(),
details::generic_set_op(),
y.begin());
// receive tmp from others
mpi_scatter_end (
x.begin(),
_buffer.begin(),
_from,
_to,
details::generic_set_op(),
tag,
row_ownership().comm());
// y += ext*tmp
csr_amux (_ext.begin(), _ext.end(), _buffer.begin(), details::generic_set_plus_op(), y.begin());
}
template<class T>
void
csr_rep<T,distributed>::trans_mult(
const vec<T,distributed>& x,
vec<T,distributed>& y)
const
{
check_macro (x.size() == nrow(), "csr.trans_mult(vec): incompatible csr("<<nrow()<<","<<ncol()<<") and vec("<<x.size()<<")");
// initialize _from and _to scatter messages
_scatter_init_guard();
y.resize (col_ownership());
// y = dia*x
std::fill (y.begin(), y.end(), T(0));
csr_cumul_trans_mult (
base::begin(),
base::end(),
x.begin(),
details::generic_set_plus_op(),
y.begin());
// buffer = ext*x
std::fill (_buffer.begin(), _buffer.end(), T(0));
csr_cumul_trans_mult (
_ext.begin(),
_ext.end(),
x.begin(),
details::generic_set_plus_op(),
_buffer.begin());
// send buffer to others parts of y (+=)
distributor::tag_type tag = distributor::get_new_tag();
mpi_scatter_begin (
_buffer.begin().operator->(),
y.begin().operator->(),
_to, // reverse mode
_from,
details::generic_set_plus_op(),
tag,
col_ownership().comm());
// receive buffer from others
mpi_scatter_end (
_buffer.begin(),
y.begin(),
_to, // reverse mode
_from,
details::generic_set_plus_op(),
tag,
col_ownership().comm());
}
// ----------------------------------------------------------------------------
// blas3: basic linear algebra
// ----------------------------------------------------------------------------
template<class T>
csr_rep<T,distributed>&
csr_rep<T,distributed>::operator*= (const T& lambda)
{
base::operator*= (lambda);
_ext.operator*= (lambda);
return *this;
}
// ----------------------------------------------------------------------------
// expression c=a+b and c=a-b with a temporary c=*this
// ----------------------------------------------------------------------------
// NOTE: cet algo pourrait servir aussi au cas diag (dans csr_seq.cc)
// a condition de mettre des pseudo-renumerotations et d'enlever
// le set.insert dans la 1ere passe. Pas forcement plus lisible...
template<class T, class BinaryOp>
void
csr_ext_add (
const csr_rep<T,sequential>& a, const std::vector<typename csr<T>::size_type>& jext_a2dis_j,
const csr_rep<T,sequential>& b, const std::vector<typename csr<T>::size_type>& jext_b2dis_j,
csr_rep<T,sequential>& c, std::vector<typename csr<T>::size_type>& jext_c2dis_j,
BinaryOp binop)
{
typedef typename csr_rep<T,distributed>::size_type size_type;
typedef typename csr_rep<T,distributed>::iterator iterator;
typedef typename csr_rep<T,distributed>::const_iterator const_iterator;
typedef typename csr_rep<T,distributed>::data_iterator data_iterator;
typedef typename csr_rep<T,distributed>::const_data_iterator const_data_iterator;
typedef std::pair<size_type,T> pair_type;
//
// first pass: compute nnz_c and resize
//
size_type nnz_ext_c = 0;
const size_type infty = std::numeric_limits<size_type>::max();
const_iterator ia = a.begin();
const_iterator ib = b.begin();
std::set<size_type> jext_c_set;
for (size_type i = 0, n = a.nrow(); i < n; i++) {
for (const_data_iterator iter_jva = ia[i], last_jva = ia[i+1],
iter_jvb = ib[i], last_jvb = ib[i+1];
iter_jva != last_jva || iter_jvb != last_jvb; ) {
size_type dis_ja = iter_jva == last_jva ? infty : jext_a2dis_j [(*iter_jva).first];
size_type dis_jb = iter_jvb == last_jvb ? infty : jext_b2dis_j [(*iter_jvb).first];
if (dis_ja == dis_jb) {
jext_c_set.insert (dis_ja);
iter_jva++;
iter_jvb++;
} else if (dis_ja < dis_jb) {
jext_c_set.insert (dis_ja);
iter_jva++;
} else {
jext_c_set.insert (dis_jb);
iter_jvb++;
}
nnz_ext_c++;
}
}
c.resize (a.nrow(), b.ncol(), nnz_ext_c);
jext_c2dis_j.resize (jext_c_set.size());
std::copy (jext_c_set.begin(), jext_c_set.end(), jext_c2dis_j.begin());
//
// second pass: add and store in c
//
op_dis_j2jext_t<pair_type, pair_type, typename vector<size_type>::const_iterator>
op_dis_j2jext_c (jext_c2dis_j.begin(), jext_c2dis_j.end());
data_iterator iter_jvc = c._data.begin().operator->();
iterator ic = c.begin();
*ic++ = iter_jvc;
for (size_type i = 0, n = a.nrow(); i < n; i++) {
for (const_data_iterator iter_jva = ia[i], last_jva = ia[i+1],
iter_jvb = ib[i], last_jvb = ib[i+1];
iter_jva != last_jva || iter_jvb != last_jvb; ) {
size_type dis_ja = iter_jva == last_jva ? infty : jext_a2dis_j [(*iter_jva).first];
size_type dis_jb = iter_jvb == last_jvb ? infty : jext_b2dis_j [(*iter_jvb).first];
if (dis_ja == dis_jb) {
*iter_jvc++ = op_dis_j2jext_c (pair_type(dis_ja, binop((*iter_jva).second, (*iter_jvb).second)));
iter_jva++;
iter_jvb++;
} else if (dis_ja < dis_jb) {
*iter_jvc++ = op_dis_j2jext_c (pair_type(dis_ja, binop((*iter_jva).second, T(0))));
iter_jva++;
} else {
*iter_jvc++ = op_dis_j2jext_c (pair_type(dis_jb, binop(T(0),(*iter_jvb).second)));
iter_jvb++;
}
}
*ic++ = iter_jvc;
}
}
template<class T>
template<class BinaryOp>
void
csr_rep<T,distributed>::assign_add (
const csr_rep<T,distributed>& a,
const csr_rep<T,distributed>& b,
BinaryOp binop)
{
check_macro (a.dis_nrow() == b.dis_nrow() && a.dis_ncol() == b.dis_ncol(),
"a+b: invalid matrix a("<<a.dis_nrow()<<","<<a.dis_ncol()<<") and b("
<<b.dis_nrow()<<","<<b.dis_ncol()<<")");
check_macro (a.nrow() == b.nrow() && a.ncol() == b.ncol(),
"a+b: matrix local distribution mismatch: a("<<a.nrow()<<","<<a.ncol()<<") and b("
<<b.nrow()<<","<<b.ncol()<<")");
// 1) the diagonal part:
base::assign_add (a, b, binop);
// 2) the extra-diagonal part:
csr_ext_add (
a._ext, a._jext2dis_j,
b._ext, b._jext2dis_j,
_ext, _jext2dis_j,
binop);
_dis_nnz = mpi::all_reduce (comm(), nnz() + _ext.nnz(), std::plus<size_type>());
_dis_ext_nnz = mpi::all_reduce (comm(), _ext.nnz(), std::plus<size_type>());
_scatter_initialized = false;
#ifdef TO_CLEAN
// 3) scatter init for a*x :
vector<size_type> id(_jext2dis_j.size());
for (size_type i = 0; i < id.size(); i++) id[i] = i;
distributor::tag_type tag = distributor::get_new_tag();
_buffer.resize (_jext2dis_j.size());
mpi_scatter_init(
_jext2dis_j.size(),
_jext2dis_j.begin().operator->(),
id.size(),
id.begin().operator->(),
dis_ncol(),
col_ownership().begin().operator->(),
tag,
row_ownership().comm(),
_from,
_to);
#endif // TO_CLEAN
}
// ----------------------------------------------------------------------------
// trans(a)
// ----------------------------------------------------------------------------
template<class T>
void
csr_rep<T,distributed>::build_transpose (csr_rep<T,distributed>& b) const
{
//
// first: assembly all _ext parts of the a matrix in b=trans(a)
//
asr<T> b_ext (col_ownership(), row_ownership());
size_type first_i = row_ownership().first_index();
const_iterator ext_ia = ext_begin();
for (size_type i = 0, n = nrow(); i < n; i++) {
size_type dis_i = first_i + i;
for (const_data_iterator p = ext_ia[i], last_p = ext_ia[i+1]; p < last_p; p++) {
size_type dis_j = jext2dis_j ((*p).first);
const T& val = (*p).second;
b_ext.dis_entry (dis_j, dis_i) += val;
}
}
b_ext.dis_entry_assembly();
b.build_from_asr (b_ext);
//
// second: add all _diag parts
//
base::build_transpose (b);
//
// third: update dis_nnz by adding all diag nnz
//
b._dis_nnz = mpi::all_reduce (comm(), b.nnz() + b._ext.nnz(), std::plus<size_type>());
b._dis_ext_nnz = mpi::all_reduce (comm(), b._ext.nnz(), std::plus<size_type>());
b._scatter_initialized = false;
}
// ----------------------------------------------------------------------------
// set symmetry by check
// ----------------------------------------------------------------------------
template<class T>
void
csr_rep<T,distributed>::set_symmetry_by_check (const T& tol) const
{
if (dis_nrow() != dis_ncol()) {
set_symmetry(false);
return;
}
// check if |a-trans(a)| == 0 up to machine prec
csr_rep<T,distributed> at, d;
build_transpose (at);
d.assign_add (*this, at, std::minus<T>());
set_symmetry (d.max_abs() <= tol);
}
// ----------------------------------------------------------------------------
// expression c=a*b with a temporary c=*this
// ----------------------------------------------------------------------------
template<class T>
void
csr_rep<T,distributed>::assign_mult (
const csr_rep<T,distributed>& a,
const csr_rep<T,distributed>& b)
{
// TODO: distributed csr*csr could be simplified
// TODO: when a._jext2dis_j.size()==0 on all procs: no comms needed
// TODO: and the sequential algo could be used
check_macro (a.col_ownership() == b.row_ownership(),
"incompatible csr([0:"<<a.nrow()<<"|"<<a.dis_nrow()<<"[x"
<<"[0:"<<a.ncol()<<"|"<<a.dis_ncol()<<"[)"
"*csr([0:"<<b.nrow()<<"|"<<b.dis_nrow()<<"[x"
<<"[0:"<<b.ncol()<<"|"<<b.dis_ncol()<<"[)");
//
// 1) compress a non-empty column indexes numbering: creates a_zip_j numbering
//
size_type a_dia_ncol = a.col_ownership().size();
size_type a_ext_ncol = a._jext2dis_j.size(); // number of external columns; compressed
size_type first_a_dis_j = a.col_ownership().first_index();
std::map<size_type,size_type> a_dis_j2a_zip_j;
size_type a_zip_j = 0;
for (size_type jext = 0; jext < a_ext_ncol && a._jext2dis_j [jext] < first_a_dis_j; jext++, a_zip_j++) {
// row < local row index
size_type dis_j = a._jext2dis_j [jext];
a_dis_j2a_zip_j [dis_j] = a_zip_j;
}
size_type jext_up = a_zip_j;
for (size_type j = 0; j < a_dia_ncol; j++, a_zip_j++) {
// local rows
size_type dis_j = first_a_dis_j + j;
a_dis_j2a_zip_j [dis_j] = a_zip_j;
}
for (size_type jext = jext_up; jext < a_ext_ncol; jext++, a_zip_j++) {
// row > local row index
size_type dis_j = a._jext2dis_j [jext];
a_dis_j2a_zip_j [dis_j] = a_zip_j;
}
size_type a_zip_ncol = a_dia_ncol + a_ext_ncol;
//
// 2) create a seq matrix B_seq = submatrix of B by taking rows of B that equal to nonzero col of A
// requiers comms
//
typedef std::allocator<T> A; // TODO: use heap_alloc for asr
A alloc;
// 2.1) convert b in asr distributed format
asr<T,distributed,A> b_asr (b, alloc);
// 2.2) make available all external columns that are used in a
b_asr.set_dis_indexes (a._jext2dis_j);
// 2.3) convert local part of b_asr into b_asr_zip
distributor a_zip_col_ownership (distributor::decide, b.comm(), a_zip_ncol);
distributor b_zip_row_ownership = a_zip_col_ownership;
asr<T,sequential,A> b_asr_zip (b_zip_row_ownership, b.col_ownership());
size_type first_dis_j = b.row_ownership().first_index();
size_type j = 0;
for (typename asr<T,distributed,A>::const_iterator
iter = b_asr.begin(),
last = b_asr.end(); iter != last; ++iter, ++j) {
typedef typename asr<T,distributed,A>::row_type row_type;
size_type dis_j = first_dis_j + j;
size_type zip_j = a_dis_j2a_zip_j [dis_j];
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) {
size_type dis_k = (*row_iter).first;
const T& value = (*row_iter).second;
b_asr_zip.semi_dis_entry (zip_j, dis_k) = value;
}
}
// 2.4) convert external part of b_asr into b_asr_zip
typedef typename asr<T,distributed,A>::scatter_map_type ext_row_type;
const ext_row_type& b_ext_row_map = b_asr.get_dis_map_entries();
for (typename ext_row_type::const_iterator
iter = b_ext_row_map.begin(),
last = b_ext_row_map.end(); iter != last; ++iter) {
typedef typename asr<T,distributed,A>::row_type row_type;
size_type dis_j = (*iter).first;
size_type zip_j = a_dis_j2a_zip_j [dis_j];
const row_type& row = (*iter).second;
for (typename row_type::const_iterator
row_iter = row.begin(),
row_last = row.end(); row_iter != row_last; ++row_iter) {
size_type dis_k = (*row_iter).first;
const T& value = (*row_iter).second;
b_asr_zip.semi_dis_entry (zip_j, dis_k) = value;
}
}
b_asr_zip.dis_entry_assembly(); // recount nnz
// 2.5) convert b_asr_zip into csr: b_zip
csr<T,sequential> b_zip (b_asr_zip);
//
// 3) create a seq matrix a_zip = submatrix of A by compressing all local cols of A
// no comms requiered: a_zip := zip(a_dia, a_ext)
//
// 3.1) assembly a local asr matrix
asr<T,sequential> a_asr_zip (a.row_ownership(), a_zip_col_ownership);
size_type first_dis_i = a.row_ownership().first_index();
const_iterator dia_ia = a.begin();
const_iterator ext_ia = a.ext_begin();
for (size_type i = 0, n = a.nrow(); i < n; ++i) {
size_type dis_i = first_dis_i + i;
for (const_data_iterator p = dia_ia[i], last_p = dia_ia[i+1]; p != last_p; ++p) {
size_type j = (*p).first;
const T& value = (*p).second;
size_type dis_j = first_dis_j + j;
size_type zip_j = a_dis_j2a_zip_j [dis_j];
a_asr_zip.semi_dis_entry (i, zip_j) += value;
}
if (a.ext_nnz() == 0) continue;
for (const_data_iterator p = ext_ia[i], last_p = ext_ia[i+1]; p != last_p; ++p) {
size_type jext = (*p).first;
const T& value = (*p).second;
size_type dis_j = a._jext2dis_j [jext];
size_type zip_j = a_dis_j2a_zip_j [dis_j];
a_asr_zip.semi_dis_entry (i, zip_j) += value;
}
}
a_asr_zip.dis_entry_assembly();
// 3.2) convert to sequential csr
csr<T,sequential> a_zip (a_asr_zip);
//
// 4) compute the product sequentially (no comms requiered)
//
csr<T,sequential> c_zip = a_zip*b_zip;
//
// 5) create mpi matrix C by concatenating C_tmp on all procs
// no comms requiered: C = (C_dia,C_ext) = dispatch(C_tmp)
// 5.1) assembly in distributed asr
asr<T,distributed> c_asr_zip (a.row_ownership(), b.col_ownership());
const_iterator ic = c_zip.begin();
for (size_type i = 0, n = c_zip.nrow(); i < n; ++i) {
size_type dis_i = first_dis_i + i;
for (const_data_iterator p = ic[i], last_p = ic[i+1]; p != last_p; ++p) {
size_type dis_k = (*p).first;
const T& value = (*p).second;
c_asr_zip.semi_dis_entry (i, dis_k) += value;
}
}
c_asr_zip.dis_entry_assembly();
// 5.2) copy c_zip intro c_asr_zip
build_from_asr (c_asr_zip);
// 5.*) propagate flags ; cannot stat yet for symmetry and positive_definite
set_pattern_dimension (std::max(a.pattern_dimension(), b.pattern_dimension()));
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#ifndef _RHEOLEF_USE_HEAP_ALLOCATOR
#define _RHEOLEF_instanciate_class(T) \
template void csr_rep<T,distributed>::build_from_asr (const asr<T,distributed,std::allocator<T> >&);
#else // _RHEOLEF_USE_HEAP_ALLOCATOR
#define _RHEOLEF_instanciate_class(T) \
template void csr_rep<T,distributed>::build_from_asr (const asr<T,distributed,std::allocator<T> >&); \
template void csr_rep<T,distributed>::build_from_asr (const asr<T,distributed,heap_allocator<T> >&);
#endif // _RHEOLEF_USE_HEAP_ALLOCATOR
#define _RHEOLEF_istanciate(T) \
template class csr_rep<T,distributed>; \
template void csr_rep<T,distributed>::assign_add ( \
const csr_rep<T,distributed>&, \
const csr_rep<T,distributed>&, \
std::plus<T>); \
template void csr_rep<T,distributed>::assign_add ( \
const csr_rep<T,distributed>&, \
const csr_rep<T,distributed>&, \
std::minus<T>); \
_RHEOLEF_instanciate_class(T)
_RHEOLEF_istanciate(Float)
} // namespace rheolef
# endif // _RHEOLEF_HAVE_MPI
size_type
field::size_type size_type
Definition branch.cc:430
Float
see the Float page for the full documentation
rheolef::asr
Definition asr.h:49
rheolef::asr::semi_dis_entry
T & semi_dis_entry(size_type i, size_type dis_j)
Definition asr.h:183
rheolef::asr::dis_entry_assembly
void dis_entry_assembly()
Definition asr.h:112
rheolef::asr::row_type
pair_set< T, A > row_type
Definition asr.h:53
rheolef::asr::dis_entry
dis_reference dis_entry(size_type dis_i, size_type dis_j)
Definition asr.h:193
rheolef::csr_rep< T, distributed >::size_type
base::size_type size_type
Definition csr.h:191
rheolef::csr_rep< T, distributed >::const_iterator
base::const_iterator const_iterator
Definition csr.h:194
rheolef::csr_rep< T, distributed >::const_data_iterator
base::const_data_iterator const_data_iterator
Definition csr.h:196
rheolef::csr_rep< T, sequential >
Definition csr.h:84
rheolef::csr_rep< T, sequential >::pair_type
std::pair< size_type, T > pair_type
Definition csr.h:89
rheolef::csr_rep
Definition csr.h:81
rheolef::csr
see the csr page for the full documentation
Definition csr.h:317
rheolef::disarray::const_iterator
rep::base::const_iterator const_iterator
Definition disarray.h:503
rheolef::distributor
see the distributor page for the full documentation
Definition distributor.h:69
rheolef::distributor::get_new_tag
static tag_type get_new_tag()
returns a new tag
Definition distributor.cc:133
rheolef::distributor::decide
static const size_type decide
Definition distributor.h:83
rheolef::distributor::tag_type
int tag_type
Definition distributor.h:78
rheolef::odiststream
odiststream: see the diststream page for the full documentation
Definition diststream.h:137
rheolef::odiststream::open
void open(std::string filename, std::string suffix="", io::mode_type mode=io::out, const communicator &comm=communicator())
This routine opens a physical output file.
Definition diststream.cc:144
rheolef::odiststream::good
bool good() const
Definition diststream.cc:195
rheolef::odiststream::io_proc
static size_type io_proc()
Definition diststream.cc:79
rheolef::vec
see the vec page for the full documentation
Definition vec.h:79
rheolef::vec::resize
void resize(const distributor &ownership, const T &init_val=std::numeric_limits< T >::max())
Definition vec.h:199
_RHEOLEF_istanciate
#define _RHEOLEF_istanciate(T)
Definition csr_mpi.cc:701
assert_macro
#define assert_macro(ok_condition, message)
Definition dis_macros.h:113
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::asr_to_csr_dist_logical
Set::value_type asr_to_csr_dist_logical(InputPtrIterator iter_ptr_a, InputPtrIterator last_ptr_a, Predicate is_dia, Set &colext)
Definition asr_to_csr_dist_logical.h:73
rheolef::csr_amux
void csr_amux(InputIterator ia, InputIterator last_ia, InputRandomAcessIterator x, SetOperator set_op, OutputIterator y)
Definition csr_amux.h:77
rheolef::put
void put(std::ostream &out, std::string name, const tiny_matrix< T > &a)
Definition tiny_lu.h:155
rheolef::mpi_scatter_end
void mpi_scatter_end(InputIterator x, OutputIterator y, Message &from, Message &to, SetOp op, Tag tag, Comm comm)
Definition mpi_scatter_end.h:60
rheolef::mpi_scatter_begin
void mpi_scatter_begin(InputIterator x, OutputIterator y, Message &from, Message &to, SetOp op, Tag tag, Comm comm)
Definition mpi_scatter_begin.h:196
rheolef::mpi_scatter_init
void mpi_scatter_init(Size nidx, SizeRandomIterator1 idx, Size nidy, SizeRandomIterator2 idy, Size idy_maxval, SizeRandomIterator3 ownership, Tag tag, const distributor::communicator_type &comm, Message &from, Message &to)
Definition mpi_scatter_init.h:73
rheolef::csr_ext_add
void csr_ext_add(const csr_rep< T, sequential > &a, const std::vector< typename csr< T >::size_type > &jext_a2dis_j, const csr_rep< T, sequential > &b, const std::vector< typename csr< T >::size_type > &jext_b2dis_j, csr_rep< T, sequential > &c, std::vector< typename csr< T >::size_type > &jext_c2dis_j, BinaryOp binop)
Definition csr_mpi.cc:367
rheolef::csr_cumul_trans_mult
void csr_cumul_trans_mult(InputIterator1 ia, InputIterator1 last_ia, InputIterator3 x, SetOperator set_op, RandomAccessMutableIterator y)
Definition csr_cumul_trans_mult.h:67
rheolef::asr_to_csr
OutputPtrIterator asr_to_csr(InputPtrIterator iter_ptr_a, InputPtrIterator last_ptr_a, Predicate pred, Operation op, OutputPtrIterator iter_ptr_b, OutputDataIterator iter_data_b)
Definition asr_to_csr.h:70
rheolef::operator()
t operator()(const t &a, const t &b)
Definition space.cc:386
std
STL namespace.
d
size_t d
Definition neumann-laplace-lambda.cc:28
A
Definition diffusion_isotropic.h:25
lambda
Definition yield_slip_circle.h:34
p
Definition sphere.icc:25
rheolef::const_iterator
Definition field_expr_recursive.h:594
rheolef::is_dia_t
Definition asr_to_csr_dist_logical.h:99