23 #include "rheolef/csr_concat.h"
25 namespace rheolef {
namespace details {
35 template <
class T,
class M>
38 std::pair<size_type,size_type>& row_sizes,
39 std::pair<size_type,size_type>& col_sizes,
40 const std::pair<size_type,size_type>& new_row_sizes,
41 const std::pair<size_type,size_type>& new_col_sizes)
43 if (row_sizes.first == undef) {
44 row_sizes = new_row_sizes;
45 }
else if (new_row_sizes.first != undef) {
46 check_macro (row_sizes.first == new_row_sizes.first,
47 "matrix initializer list: matrix row argument [0:"
48 << new_row_sizes.first <<
"|" << new_row_sizes.second <<
"["
49 <<
" has incompatible size: expect [0:"
50 << row_sizes.first <<
"|" << row_sizes.second <<
"[");
52 if (col_sizes.first == undef) {
53 col_sizes = new_col_sizes;
54 }
else if (new_col_sizes.first != undef) {
55 check_macro (col_sizes.first == new_col_sizes.first,
56 "matrix initializer list: matrix col argument [0:"
57 << new_col_sizes.first <<
"|" << new_col_sizes.second <<
"["
58 <<
" has incompatible size: expect [0:"
59 << col_sizes.first <<
"|" << col_sizes.second <<
"[");
62 template <
class T,
class M>
65 std::pair<size_type,size_type>& sizes,
66 const communicator& comm)
69 if (sizes.first == undef) {
71 size_type iproc0 = constraint_process_rank (comm);
72 sizes.first = (my_proc == iproc0) ? 1 : 0;
76 template <
class T,
class M>
79 std::vector<std::pair<size_type,size_type> >& sizes,
80 const communicator& comm)
82 for (
size_type i_comp = 0; i_comp < sizes.size(); i_comp++) {
83 finalize_sizes (sizes [i_comp], comm);
89 template <
class T,
class M>
92 std::pair<size_type,size_type>& row_sizes,
93 std::vector<std::pair<size_type,size_type> >& col_sizes,
94 communicator& comm)
const
97 size_type iproc0 = constraint_process_rank (comm);
98 if (col_sizes.size() == 0) {
99 col_sizes.resize (_l.size(), std::pair<size_type,size_type>(undef,undef));
102 "matrix initializer list: unexpected matrix line size "
103 << _l.size() <<
": size " << col_sizes.size()
107 for (
typename std::list<value_type>::const_iterator iter = _l.begin(); iter != _l.end(); ++iter, j_comp++) {
110 case value_type::empty: {
113 trace_macro (
"block col j_comp="<<j_comp<<
": e: "<<nrows.first<<
"x"<<ncols.first<<
"...");
114 set_sizes (row_sizes, col_sizes[j_comp], nrows, ncols);
115 trace_macro (
"block col j_comp="<<j_comp<<
": E: "<<row_sizes.first<<
"x"<<col_sizes[j_comp].first);
126 nrow = (my_proc == iproc0) ? 1 : 0;
128 ncol = (my_proc == iproc0) ? 1 : 0;
132 trace_macro (
"block col j_comp="<<j_comp<<
": s: "<<nrow<<
"x"<<ncol<<
"...");
133 set_sizes (row_sizes, col_sizes[j_comp], std::make_pair (nrow, dis_nrow), std::make_pair (ncol, dis_ncol));
134 trace_macro (
"block col j_comp="<<j_comp<<
": S: "<<row_sizes.first<<
"x"<<col_sizes[j_comp].first);
140 size_type dis_nrow = x.
v.ownership().dis_size();
141 size_type ncol = (my_proc == iproc0) ? 1 : 0;
143 trace_macro (
"block col j_comp="<<j_comp<<
": v: "<<nrow<<
"x"<<ncol<<
"...");
144 set_sizes (row_sizes, col_sizes[j_comp], std::make_pair (nrow, dis_nrow), std::make_pair (ncol, dis_ncol));
145 trace_macro (
"block col j_comp="<<j_comp<<
": V: "<<row_sizes.first<<
"x"<<col_sizes[j_comp].first);
149 case value_type::vector_transpose: {
152 size_type dis_ncol = x.
v.ownership().dis_size();
153 size_type nrow = (my_proc == iproc0) ? 1 : 0;
155 set_sizes (row_sizes, col_sizes[j_comp], std::make_pair (nrow, dis_nrow), std::make_pair (ncol, dis_ncol));
156 trace_macro (
"block col j_comp="<<j_comp<<
": VT: "<<row_sizes.first<<
"x"<<col_sizes[j_comp].first);
161 case value_type::vector_vec: {
165 if (x.
vv.size() > 0) {
166 ncol = x.
vv[0].ownership().size();
167 dis_ncol = x.
vv[0].ownership().dis_size();
168 comm = x.
vv[0].ownership().comm();
170 size_type nrow = (my_proc == iproc0) ? x.
vv.size() : 0;
172 set_sizes (row_sizes, col_sizes[j_comp], std::make_pair (nrow, dis_nrow), std::make_pair (ncol, dis_ncol));
173 trace_macro (
"block col j_comp="<<j_comp<<
": VV: "<<row_sizes.first<<
"x"<<col_sizes[j_comp].first);
176 case value_type::vector_vec_transpose: {
180 if (x.
vv.size() > 0) {
181 nrow = x.
vv[0].ownership().size();
182 dis_nrow = x.
vv[0].ownership().dis_size();
183 comm = x.
vv[0].ownership().comm();
185 size_type ncol = (my_proc == iproc0) ? x.
vv.size() : 0;
187 set_sizes (row_sizes, col_sizes[j_comp], std::make_pair (nrow, dis_nrow), std::make_pair (ncol, dis_ncol));
188 trace_macro (
"block col j_comp="<<j_comp<<
": VVT: "<<row_sizes.first<<
"x"<<col_sizes[j_comp].first);
191 case value_type::matrix: {
193 size_type dis_nrow = x.
m.row_ownership().dis_size();
195 size_type dis_ncol = x.
m.col_ownership().dis_size();
196 set_sizes (row_sizes, col_sizes[j_comp], std::make_pair (nrow, dis_nrow), std::make_pair (ncol, dis_ncol));
197 trace_macro (
"block col j_comp="<<j_comp<<
": m: "<<nrow<<
"x"<<ncol<<
"...");
198 comm = x.
m.row_ownership().comm();
199 trace_macro (
"block col j_comp="<<j_comp<<
": M: "<<row_sizes.first<<
"x"<<col_sizes[j_comp].first);
204 trace_macro (
"block col j_comp="<<j_comp<<
" done");
210 template <
class T,
class M>
213 const std::pair<size_type,size_type>& row_sizes,
214 const std::vector<std::pair<size_type,size_type> >& col_sizes,
215 const communicator& comm,
218 std::vector<distributor>& col_start_by_component)
const
220 row_ownership =
distributor (row_sizes.second, comm, row_sizes.first);
223 col_start_by_component.resize (col_sizes.size());
224 for (
size_type j_comp = 0; j_comp < col_sizes.size(); j_comp++) {
225 col_start_by_component [j_comp] =
distributor (col_comp_start_dis_j, comm, col_comp_start_j);
226 col_comp_start_j += col_sizes [j_comp].first;
227 col_comp_start_dis_j += col_sizes [j_comp].second;
229 col_ownership =
distributor (col_comp_start_dis_j, comm, col_comp_start_j);
234 template <
class T,
class M>
238 const std::pair<size_type,size_type>& row_sizes,
240 const std::vector<std::pair<size_type,size_type> >& col_sizes,
241 const std::vector<distributor>& col_start_by_component)
const
245 size_type iproc0 = constraint_process_rank (comm);
247 for (
typename std::list<value_type>::const_iterator iter = _l.begin(); iter != _l.end(); ++iter, j_comp++) {
250 case value_type::empty: {
257 + row_start_by_component.
size(iproc0);
259 + col_start_by_component[j_comp].size(iproc0);
260 if (x.
s != 0 && my_proc == iproc0) A.
dis_entry (dis_i, dis_j) += x.
s;
266 + row_start_by_component.
size();
268 + col_start_by_component[j_comp].size(iproc0);
270 if (*iter != 0) A.
dis_entry (dis_i, dis_j) += *iter;
274 case value_type::vector_transpose: {
277 + row_start_by_component.
size(iproc0);
279 + col_start_by_component[j_comp].size();
281 if (*iter != 0) A.
dis_entry (dis_i, dis_j) += *iter;
285 case value_type::vector_vec: {
290 + row_start_by_component.
size(iproc0)
293 + col_start_by_component[j_comp].size();
295 if (*iter != 0) A.
dis_entry (dis_i, dis_j) += *iter;
300 case value_type::vector_vec_transpose: {
305 + row_start_by_component.
size();
307 + col_start_by_component[j_comp].size(iproc0)
310 if (*iter != 0) A.
dis_entry (dis_i, dis_j) += *iter;
315 case value_type::matrix: {
320 + row_start_by_component.
size()
325 + col_start_by_component[j_comp].size()
327 A.
dis_entry (dis_i, dis_j) += (*p).second;
331 size_type comp_dis_j = x.
m.jext2dis_j ((*p).first);
332 size_type iproc = x.
m.col_ownership().find_owner (comp_dis_j);
333 size_type first_comp_dis_j_iproc = x.
m.col_ownership().first_index (iproc);
334 assert_macro (comp_dis_j >= first_comp_dis_j_iproc,
"invalid index");
335 size_type j = comp_dis_j - first_comp_dis_j_iproc;
337 + col_start_by_component[j_comp].size(iproc)
339 A.
dis_entry (dis_i, dis_j) += (*p).second;
351 template <
class T,
class M>
356 std::pair<size_type,size_type> row_sizes = std::pair<size_type,size_type>(undef,undef);
358 std::vector<std::pair<size_type,size_type> > col_sizes;
359 std::vector<distributor> col_start_by_component;
360 build_csr_pass0 (row_sizes, col_sizes, comm);
361 finalize_sizes (row_sizes, comm);
362 finalize_sizes (col_sizes, comm);
365 build_csr_pass1 (row_sizes, col_sizes, comm, row_ownership, col_ownership, col_start_by_component);
366 asr<T,M> A (row_ownership, col_ownership);
367 build_csr_pass2 (A, row_sizes, row_start_by_component, col_sizes, col_start_by_component);
376 template <
class T,
class M>
383 std::vector<std::pair<size_type,size_type> > row_sizes (_l.size(), std::pair<size_type,size_type>(undef,undef));
384 std::vector<std::pair<size_type,size_type> > col_sizes;
387 for (
typename std::list<line_type>::const_iterator iter = _l.begin(); iter != _l.end(); ++iter, i_comp++) {
391 trace_macro (
"block line i_comp="<<i_comp<<
" done");
400 std::vector<distributor> row_start_by_component (row_sizes.size()+1);
401 for (
size_type i_comp = 0; i_comp <= row_sizes.size(); i_comp++) {
402 row_start_by_component [i_comp] =
distributor (row_comp_start_dis_i, comm, row_comp_start_i);
403 if (i_comp == row_sizes.size())
break;
404 row_comp_start_i += row_sizes [i_comp].first;
405 row_comp_start_dis_i += row_sizes [i_comp].second;
407 distributor row_ownership = row_start_by_component[row_sizes.size()];
411 std::vector<distributor> col_start_by_component (col_sizes.size()+1);
412 for (
size_type j_comp = 0; j_comp <= col_sizes.size(); j_comp++) {
413 col_start_by_component [j_comp] =
distributor (col_comp_start_dis_j, comm, col_comp_start_j);
414 if (j_comp == col_sizes.size())
break;
415 col_comp_start_j += col_sizes [j_comp].first;
416 col_comp_start_dis_j += col_sizes [j_comp].second;
418 distributor col_ownership = col_start_by_component[col_sizes.size()];
422 asr<T,M> a (row_ownership, col_ownership);
424 for (
typename std::list<line_type>::const_iterator iter = _l.begin(); iter != _l.end(); ++iter, ++i_comp) {
426 line.
build_csr_pass2 (
a, row_sizes[i_comp], row_start_by_component[i_comp], col_sizes, col_start_by_component);
428 a.dis_entry_assembly();
434 #define _RHEOLEF_instanciation(T,M) \
435 template class csr_concat_line<T,M>; \
436 template class csr_concat<T,M>;
439 #ifdef _RHEOLEF_HAVE_MPI
see the Float page for the full documentation
void dis_entry_assembly()
const distributor & col_ownership() const
dis_reference dis_entry(size_type dis_i, size_type dis_j)
const distributor & row_ownership() const
see the csr page for the full documentation
csr< T, M > build_csr() const
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
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
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)
csr_concat_value< T, M >::size_type size_type
csr_concat_value< T, M >::sizes_type sizes_type
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
static void finalize_sizes(std::pair< size_type, size_type > &sizes, const communicator &comm)
std::vector< vec< T, M > > vv
csr< T, M > build_csr() const
csr_concat_value< T, M >::size_type size_type
see the distributor page for the full documentation
size_type size(size_type iproc) const
size_type first_index(size_type iproc) const
global index range and local size owned by ip-th process
const communicator_type & comm() const
see the vec page for the full documentation
#define trace_macro(message)
#define assert_macro(ok_condition, message)
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
_RHEOLEF_instanciation(Float, sequential) _RHEOLEF_instanciation(Float
This file is part of Rheolef.