Loading...
Searching...
No Matches
csr_concat.cc
Go to the documentation of this file.
// build csr from initializer list (c++ 2011)
//
#include "rheolef/csr_concat.h"
namespace rheolef { namespace details {
// =========================================================================
// 1rst case : one-line matrix initializer
// A = {a, u}; // matrix & vector
// A = {trans(u), 3.2}; // trans(vector) & scalar
// =========================================================================
// ------------------------
// sizes utilities
// ------------------------
template <class T, class M>
void
csr_concat_line<T,M>::set_sizes (
std::pair<size_type,size_type>& row_sizes,
std::pair<size_type,size_type>& col_sizes,
const std::pair<size_type,size_type>& new_row_sizes,
const std::pair<size_type,size_type>& new_col_sizes)
{
if (row_sizes.first == undef) {
row_sizes = new_row_sizes;
} else if (new_row_sizes.first != undef) {
check_macro (row_sizes.first == new_row_sizes.first,
"matrix initializer list: matrix row argument [0:"
<< new_row_sizes.first << "|" << new_row_sizes.second << "["
<< " has incompatible size: expect [0:"
<< row_sizes.first << "|" << row_sizes.second << "[");
}
if (col_sizes.first == undef) {
col_sizes = new_col_sizes;
} else if (new_col_sizes.first != undef) {
check_macro (col_sizes.first == new_col_sizes.first,
"matrix initializer list: matrix col argument [0:"
<< new_col_sizes.first << "|" << new_col_sizes.second << "["
<< " has incompatible size: expect [0:"
<< col_sizes.first << "|" << col_sizes.second << "[");
}
}
template <class T, class M>
void
csr_concat_line<T,M>::finalize_sizes (
std::pair<size_type,size_type>& sizes,
const communicator& comm)
{
// when "0" has still unknown block size, set it to 1
if (sizes.first == undef) {
size_type my_proc = comm.rank();
size_type iproc0 = constraint_process_rank (comm);
sizes.first = (my_proc == iproc0) ? 1 : 0;
sizes.second = 1;
}
}
template <class T, class M>
void
csr_concat_line<T,M>::finalize_sizes (
std::vector<std::pair<size_type,size_type> >& sizes,
const communicator& comm)
{
for (size_type i_comp = 0; i_comp < sizes.size(); i_comp++) {
finalize_sizes (sizes [i_comp], comm);
}
}
// ------------------------
// pass 0 : compute sizes
// ------------------------
template <class T, class M>
void
csr_concat_line<T,M>::build_csr_pass0 (
std::pair<size_type,size_type>& row_sizes,
std::vector<std::pair<size_type,size_type> >& col_sizes,
communicator& comm) const
{
size_type my_proc = comm.rank();
size_type iproc0 = constraint_process_rank (comm);
if (col_sizes.size() == 0) {
col_sizes.resize (_l.size(), std::pair<size_type,size_type>(undef,undef));
} else {
check_macro (col_sizes.size() == _l.size(),
"matrix initializer list: unexpected matrix line size "
<< _l.size() << ": size " << col_sizes.size()
<< " was expected");
}
size_type j_comp = 0;
for (typename std::list<value_type>::const_iterator iter = _l.begin(); iter != _l.end(); ++iter, j_comp++) {
const value_type& x = *iter;
switch (x.variant) {
case value_type::empty: {
const sizes_type& nrows = x.e.first;
const sizes_type& ncols = x.e.second;
trace_macro ("block col j_comp="<<j_comp<<": e: "<<nrows.first<<"x"<<ncols.first<<"...");
set_sizes (row_sizes, col_sizes[j_comp], nrows, ncols);
trace_macro ("block col j_comp="<<j_comp<<": E: "<<row_sizes.first<<"x"<<col_sizes[j_comp].first);
break;
}
case value_type::scalar: {
size_type nrow = undef;
size_type dis_nrow = undef;
size_type ncol = undef;
size_type dis_ncol = undef;
#ifdef TODO
// TODO: non-zero value => force size to one ?
if (x.s != 0) {
nrow = (my_proc == iproc0) ? 1 : 0;
dis_nrow = 1;
ncol = (my_proc == iproc0) ? 1 : 0;
dis_ncol = 1;
}
#endif // TODO
trace_macro ("block col j_comp="<<j_comp<<": s: "<<nrow<<"x"<<ncol<<"...");
set_sizes (row_sizes, col_sizes[j_comp], std::make_pair (nrow, dis_nrow), std::make_pair (ncol, dis_ncol));
trace_macro ("block col j_comp="<<j_comp<<": S: "<<row_sizes.first<<"x"<<col_sizes[j_comp].first);
break;
}
case value_type::vector: {
// vertical vec
size_type nrow = x.v.ownership().size();
size_type dis_nrow = x.v.ownership().dis_size();
size_type ncol = (my_proc == iproc0) ? 1 : 0;
size_type dis_ncol = 1;
trace_macro ("block col j_comp="<<j_comp<<": v: "<<nrow<<"x"<<ncol<<"...");
set_sizes (row_sizes, col_sizes[j_comp], std::make_pair (nrow, dis_nrow), std::make_pair (ncol, dis_ncol));
trace_macro ("block col j_comp="<<j_comp<<": V: "<<row_sizes.first<<"x"<<col_sizes[j_comp].first);
comm = x.v.comm(); // TODO: check if compatible comm
break;
}
case value_type::vector_transpose: {
// horizontal vec
size_type ncol = x.v.ownership().size();
size_type dis_ncol = x.v.ownership().dis_size();
size_type nrow = (my_proc == iproc0) ? 1 : 0;
size_type dis_nrow = 1;
set_sizes (row_sizes, col_sizes[j_comp], std::make_pair (nrow, dis_nrow), std::make_pair (ncol, dis_ncol));
trace_macro ("block col j_comp="<<j_comp<<": VT: "<<row_sizes.first<<"x"<<col_sizes[j_comp].first);
comm = x.v.comm(); // TODO: check if compatible comm
comm = x.v.comm(); // TODO: check if compatible comm
break;
}
case value_type::vector_vec: {
// horizontal vector<vec>
size_type ncol = undef;
size_type dis_ncol = undef;
if (x.vv.size() > 0) {
ncol = x.vv[0].ownership().size();
dis_ncol = x.vv[0].ownership().dis_size();
comm = x.vv[0].ownership().comm(); // TODO: check if compatible comm
}
size_type nrow = (my_proc == iproc0) ? x.vv.size() : 0;
size_type dis_nrow = x.vv.size();
set_sizes (row_sizes, col_sizes[j_comp], std::make_pair (nrow, dis_nrow), std::make_pair (ncol, dis_ncol));
trace_macro ("block col j_comp="<<j_comp<<": VV: "<<row_sizes.first<<"x"<<col_sizes[j_comp].first);
break;
}
case value_type::vector_vec_transpose: {
// vertical vector<vec>
size_type nrow = undef;
size_type dis_nrow = undef;
if (x.vv.size() > 0) {
nrow = x.vv[0].ownership().size();
dis_nrow = x.vv[0].ownership().dis_size();
comm = x.vv[0].ownership().comm(); // TODO: check if compatible comm
}
size_type ncol = (my_proc == iproc0) ? x.vv.size() : 0;
size_type dis_ncol = x.vv.size();
set_sizes (row_sizes, col_sizes[j_comp], std::make_pair (nrow, dis_nrow), std::make_pair (ncol, dis_ncol));
trace_macro ("block col j_comp="<<j_comp<<": VVT: "<<row_sizes.first<<"x"<<col_sizes[j_comp].first);
break;
}
case value_type::matrix: {
size_type nrow = x.m.row_ownership().size();
size_type dis_nrow = x.m.row_ownership().dis_size();
size_type ncol = x.m.col_ownership().size();
size_type dis_ncol = x.m.col_ownership().dis_size();
set_sizes (row_sizes, col_sizes[j_comp], std::make_pair (nrow, dis_nrow), std::make_pair (ncol, dis_ncol));
trace_macro ("block col j_comp="<<j_comp<<": m: "<<nrow<<"x"<<ncol<<"...");
comm = x.m.row_ownership().comm(); // TODO: check if compatible comm
trace_macro ("block col j_comp="<<j_comp<<": M: "<<row_sizes.first<<"x"<<col_sizes[j_comp].first);
break;
break;
}
}
trace_macro ("block col j_comp="<<j_comp<<" done");
}
}
// ---------------------------
// pass 1 : compute ownership
// ---------------------------
template <class T, class M>
void
csr_concat_line<T,M>::build_csr_pass1 (
const std::pair<size_type,size_type>& row_sizes,
const std::vector<std::pair<size_type,size_type> >& col_sizes,
const communicator& comm,
distributor& row_ownership,
distributor& col_ownership,
std::vector<distributor>& col_start_by_component) const
{
row_ownership = distributor (row_sizes.second, comm, row_sizes.first);
size_type col_comp_start_j = 0;
size_type col_comp_start_dis_j = 0;
col_start_by_component.resize (col_sizes.size());
for (size_type j_comp = 0; j_comp < col_sizes.size(); j_comp++) {
col_start_by_component [j_comp] = distributor (col_comp_start_dis_j, comm, col_comp_start_j);
col_comp_start_j += col_sizes [j_comp].first;
col_comp_start_dis_j += col_sizes [j_comp].second;
}
col_ownership = distributor (col_comp_start_dis_j, comm, col_comp_start_j);
}
// -----------------------------------
// pass 2 : copy at the right location
// -----------------------------------
template <class T, class M>
void
csr_concat_line<T,M>::build_csr_pass2 (
asr<T,M>& A,
const std::pair<size_type,size_type>& row_sizes,
const distributor& row_start_by_component,
const std::vector<std::pair<size_type,size_type> >& col_sizes,
const std::vector<distributor>& col_start_by_component) const
{
communicator comm = A.row_ownership().comm();
size_type my_proc = comm.rank();
size_type iproc0 = constraint_process_rank (comm);
size_type j_comp = 0;
for (typename std::list<value_type>::const_iterator iter = _l.begin(); iter != _l.end(); ++iter, j_comp++) {
const value_type& x = *iter;
switch (x.variant) {
case value_type::empty: {
// no value to insert
break;
}
case value_type::scalar: {
// do not add eplicit zero values
size_type dis_i = A.row_ownership().first_index(iproc0)
+ row_start_by_component.size(iproc0);
size_type dis_j = A.col_ownership().first_index(iproc0)
+ col_start_by_component[j_comp].size(iproc0);
if (x.s != 0 && my_proc == iproc0) A.dis_entry (dis_i, dis_j) += x.s;
break;
}
case value_type::vector: {
// vertical vec : moved column to iproc0, but still distributed by lines
size_type dis_i = A.row_ownership().first_index()
+ row_start_by_component.size();
size_type dis_j = A.col_ownership().first_index(iproc0)
+ col_start_by_component[j_comp].size(iproc0);
for (typename vec<T,M>::const_iterator iter = x.v.begin(), last = x.v.end(); iter != last; iter++, dis_i++) {
if (*iter != 0) A.dis_entry (dis_i, dis_j) += *iter;
}
break;
}
case value_type::vector_transpose: {
// horizontal vec: row integrally hosted by iproc=0 ; col still distributed
size_type dis_i = A.row_ownership().first_index(iproc0)
+ row_start_by_component.size(iproc0);
size_type dis_j = A.col_ownership().first_index()
+ col_start_by_component[j_comp].size();
for (typename vec<T,M>::const_iterator iter = x.v.begin(), last = x.v.end(); iter != last; iter++, dis_j++) {
if (*iter != 0) A.dis_entry (dis_i, dis_j) += *iter;
}
break;
}
case value_type::vector_vec: {
// horizontal vector<vec>: rows integrally hosted by iproc=0 ; col still distributed
size_type n = x.vv.size();
for (size_type k = 0; k < n; k++) {
size_type dis_i = A.row_ownership().first_index(iproc0)
+ row_start_by_component.size(iproc0)
+ k;
size_type dis_j = A.col_ownership().first_index()
+ col_start_by_component[j_comp].size();
for (typename vec<T,M>::const_iterator iter = x.vv[k].begin(), last = x.vv[k].end(); iter != last; iter++, dis_j++) {
if (*iter != 0) A.dis_entry (dis_i, dis_j) += *iter;
}
}
break;
}
case value_type::vector_vec_transpose: {
// vertical vector<vec>: moved column to iproc0, but still distributed by lines
size_type n = x.vv.size();
for (size_type k = 0; k < n; k++) {
size_type dis_i = A.row_ownership().first_index()
+ row_start_by_component.size();
size_type dis_j = A.col_ownership().first_index(iproc0)
+ col_start_by_component[j_comp].size(iproc0)
+ k;
for (typename vec<T,M>::const_iterator iter = x.vv[k].begin(), last = x.vv[k].end(); iter != last; iter++, dis_i++) {
if (*iter != 0) A.dis_entry (dis_i, dis_j) += *iter;
}
}
break;
}
case value_type::matrix: {
typename csr<T,M>::const_iterator dia_ia = x.m.begin();
typename csr<T,M>::const_iterator ext_ia = x.m.ext_begin();
for (size_type i = 0, n = x.m.nrow(); i < n; i++) {
size_type dis_i = A.row_ownership().first_index()
+ row_start_by_component.size()
+ i;
for (typename csr<T,M>::const_data_iterator p = dia_ia[i]; p < dia_ia[i+1]; p++) {
size_type j = (*p).first;
size_type dis_j = A.col_ownership().first_index()
+ col_start_by_component[j_comp].size()
+ j;
A.dis_entry (dis_i, dis_j) += (*p).second;
}
if (is_distributed<M>::value && x.m.ext_nnz() != 0) {
for (typename csr<T,M>::const_data_iterator p = ext_ia[i]; p < ext_ia[i+1]; p++) {
size_type comp_dis_j = x.m.jext2dis_j ((*p).first);
size_type iproc = x.m.col_ownership().find_owner (comp_dis_j);
size_type first_comp_dis_j_iproc = x.m.col_ownership().first_index (iproc);
assert_macro (comp_dis_j >= first_comp_dis_j_iproc, "invalid index");
size_type j = comp_dis_j - first_comp_dis_j_iproc;
size_type dis_j = A.col_ownership().first_index(iproc)
+ col_start_by_component[j_comp].size(iproc)
+ j;
A.dis_entry (dis_i, dis_j) += (*p).second;
}
}
}
break;
}
}
}
}
// ------------------------------
// merge passes 1 & 2
// ------------------------------
template <class T, class M>
csr<T,M>
csr_concat_line<T,M>::build_csr() const
{
communicator comm;
std::pair<size_type,size_type> row_sizes = std::pair<size_type,size_type>(undef,undef);
distributor row_start_by_component (0, comm, 0);
std::vector<std::pair<size_type,size_type> > col_sizes;
std::vector<distributor> col_start_by_component;
build_csr_pass0 (row_sizes, col_sizes, comm);
finalize_sizes (row_sizes, comm);
finalize_sizes (col_sizes, comm);
distributor row_ownership;
distributor col_ownership;
build_csr_pass1 (row_sizes, col_sizes, comm, row_ownership, col_ownership, col_start_by_component);
asr<T,M> A (row_ownership, col_ownership);
build_csr_pass2 (A, row_sizes, row_start_by_component, col_sizes, col_start_by_component);
A.dis_entry_assembly();
return csr<T,M>(A);
}
// =========================================================================
// 2nd case : multi-line matrix initializer
// A = { {a, u },
// {trans(u), 3.2} };
// =========================================================================
template <class T, class M>
csr<T,M>
csr_concat<T,M>::build_csr() const
{
// ---------------------------
// pass 0 : compute sizes
// ---------------------------
std::vector<std::pair<size_type,size_type> > row_sizes (_l.size(), std::pair<size_type,size_type>(undef,undef));
std::vector<std::pair<size_type,size_type> > col_sizes;
communicator comm;
size_type i_comp = 0;
for (typename std::list<line_type>::const_iterator iter = _l.begin(); iter != _l.end(); ++iter, i_comp++) {
const line_type& line = *iter;
trace_macro ("block line i_comp="<<i_comp<<"...");
line.build_csr_pass0 (row_sizes[i_comp], col_sizes, comm);
trace_macro ("block line i_comp="<<i_comp<<" done");
}
csr_concat_line<T,M>::finalize_sizes (row_sizes, comm);
csr_concat_line<T,M>::finalize_sizes (col_sizes, comm);
// ---------------------------
// pass 1 : compute ownership
// ---------------------------
size_type row_comp_start_i = 0;
size_type row_comp_start_dis_i = 0;
std::vector<distributor> row_start_by_component (row_sizes.size()+1);
for (size_type i_comp = 0; i_comp <= row_sizes.size(); i_comp++) {
row_start_by_component [i_comp] = distributor (row_comp_start_dis_i, comm, row_comp_start_i);
if (i_comp == row_sizes.size()) break;
row_comp_start_i += row_sizes [i_comp].first;
row_comp_start_dis_i += row_sizes [i_comp].second;
}
distributor row_ownership = row_start_by_component[row_sizes.size()];
size_type col_comp_start_j = 0;
size_type col_comp_start_dis_j = 0;
std::vector<distributor> col_start_by_component (col_sizes.size()+1);
for (size_type j_comp = 0; j_comp <= col_sizes.size(); j_comp++) {
col_start_by_component [j_comp] = distributor (col_comp_start_dis_j, comm, col_comp_start_j);
if (j_comp == col_sizes.size()) break;
col_comp_start_j += col_sizes [j_comp].first;
col_comp_start_dis_j += col_sizes [j_comp].second;
}
distributor col_ownership = col_start_by_component[col_sizes.size()];
// ------------------------
// pass 2 : copy
// ------------------------
asr<T,M> a (row_ownership, col_ownership);
i_comp = 0;
for (typename std::list<line_type>::const_iterator iter = _l.begin(); iter != _l.end(); ++iter, ++i_comp) {
const line_type& line = *iter;
line.build_csr_pass2 (a, row_sizes[i_comp], row_start_by_component[i_comp], col_sizes, col_start_by_component);
}
a.dis_entry_assembly();
return csr<T,M>(a);
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciation(T,M) \
template class csr_concat_line<T,M>; \
template class csr_concat<T,M>;
_RHEOLEF_instanciation(Float,sequential)
#ifdef _RHEOLEF_HAVE_MPI
_RHEOLEF_instanciation(Float,distributed)
#endif // _RHEOLEF_HAVE_MPI
}} // namespace rheolef::details
_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::dis_entry_assembly
void dis_entry_assembly()
Definition asr.h:112
rheolef::asr::col_ownership
const distributor & col_ownership() const
Definition asr.h:101
rheolef::asr::dis_entry
dis_reference dis_entry(size_type dis_i, size_type dis_j)
Definition asr.h:193
rheolef::asr::row_ownership
const distributor & row_ownership() const
Definition asr.h:100
rheolef::csr
see the csr page for the full documentation
Definition csr.h:317
rheolef::details::csr_concat_line::build_csr
csr< T, M > build_csr() const
Definition csr_concat.cc:353
rheolef::details::csr_concat_line::build_csr_pass0
void build_csr_pass0(std::pair< size_type, size_type > &row_sizes, std::vector< std::pair< size_type, size_type > > &col_sizes, communicator &comm) const
Definition csr_concat.cc:91
rheolef::details::csr_concat_line::build_csr_pass2
void build_csr_pass2(asr< T, M > &a, const std::pair< size_type, size_type > &row_sizes, const distributor &row_start_by_component, const std::vector< std::pair< size_type, size_type > > &col_sizes, const std::vector< distributor > &col_start_by_component) const
Definition csr_concat.cc:236
rheolef::details::csr_concat_line::set_sizes
static void set_sizes(std::pair< size_type, size_type > &row_sizes, std::pair< size_type, size_type > &col_sizes, const std::pair< size_type, size_type > &new_row_sizes, const std::pair< size_type, size_type > &new_col_sizes)
Definition csr_concat.cc:37
rheolef::details::csr_concat_line::size_type
csr_concat_value< T, M >::size_type size_type
Definition csr_concat.h:98
rheolef::details::csr_concat_line::sizes_type
csr_concat_value< T, M >::sizes_type sizes_type
Definition csr_concat.h:99
rheolef::details::csr_concat_line::value_type
csr_concat_value< T, M > value_type
Definition csr_concat.h:105
rheolef::details::csr_concat_line::build_csr_pass1
void build_csr_pass1(const std::pair< size_type, size_type > &row_sizes, const std::vector< std::pair< size_type, size_type > > &col_sizes, const communicator &comm, distributor &row_ownership, distributor &col_ownership, std::vector< distributor > &col_start_by_component) const
Definition csr_concat.cc:212
rheolef::details::csr_concat_line::finalize_sizes
static void finalize_sizes(std::pair< size_type, size_type > &sizes, const communicator &comm)
Definition csr_concat.cc:64
rheolef::details::csr_concat::build_csr
csr< T, M > build_csr() const
Definition csr_concat.cc:378
rheolef::details::csr_concat::size_type
csr_concat_value< T, M >::size_type size_type
Definition csr_concat.h:180
rheolef::details::csr_concat::line_type
csr_concat_line< T, M > line_type
Definition csr_concat.h:187
rheolef::distributor
see the distributor page for the full documentation
Definition distributor.h:69
rheolef::distributor::size
size_type size(size_type iproc) const
Definition distributor.h:170
rheolef::distributor::first_index
size_type first_index(size_type iproc) const
global index range and local size owned by ip-th process
Definition distributor.h:158
rheolef::distributor::comm
const communicator_type & comm() const
Definition distributor.h:152
rheolef::vec::const_iterator
base::const_iterator const_iterator
Definition vec.h:92
trace_macro
#define trace_macro(message)
Definition dis_macros.h:111
assert_macro
#define assert_macro(ok_condition, message)
Definition dis_macros.h:113
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
p
Definition sphere.icc:25
rheolef::distributed
Definition communicator.h:128
rheolef::is_distributed
Definition communicator.h:70
rheolef::sequential
Definition communicator.h:62