Rheolef  7.1
an efficient C++ finite element environment
csr_concat.cc
Go to the documentation of this file.
1 // build csr from initializer list (c++ 2011)
22 //
23 #include "rheolef/csr_concat.h"
24 
25 namespace rheolef { namespace details {
26 
27 // =========================================================================
28 // 1rst case : one-line matrix initializer
29 // A = {a, u}; // matrix & vector
30 // A = {trans(u), 3.2}; // trans(vector) & scalar
31 // =========================================================================
32 // ------------------------
33 // sizes utilities
34 // ------------------------
35 template <class T, class M>
36 void
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)
42 {
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 << "[");
51  }
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 << "[");
60  }
61 }
62 template <class T, class M>
63 void
65  std::pair<size_type,size_type>& sizes,
66  const communicator& comm)
67 {
68  // when "0" has still unknown block size, set it to 1
69  if (sizes.first == undef) {
70  size_type my_proc = comm.rank();
71  size_type iproc0 = constraint_process_rank (comm);
72  sizes.first = (my_proc == iproc0) ? 1 : 0;
73  sizes.second = 1;
74  }
75 }
76 template <class T, class M>
77 void
79  std::vector<std::pair<size_type,size_type> >& sizes,
80  const communicator& comm)
81 {
82  for (size_type i_comp = 0; i_comp < sizes.size(); i_comp++) {
83  finalize_sizes (sizes [i_comp], comm);
84  }
85 }
86 // ------------------------
87 // pass 0 : compute sizes
88 // ------------------------
89 template <class T, class M>
90 void
92  std::pair<size_type,size_type>& row_sizes,
93  std::vector<std::pair<size_type,size_type> >& col_sizes,
94  communicator& comm) const
95 {
96  size_type my_proc = comm.rank();
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));
100  } else {
101  check_macro (col_sizes.size() == _l.size(),
102  "matrix initializer list: unexpected matrix line size "
103  << _l.size() << ": size " << col_sizes.size()
104  << " was expected");
105  }
106  size_type j_comp = 0;
107  for (typename std::list<value_type>::const_iterator iter = _l.begin(); iter != _l.end(); ++iter, j_comp++) {
108  const value_type& x = *iter;
109  switch (x.variant) {
110  case value_type::empty: {
111  const sizes_type& nrows = x.e.first;
112  const sizes_type& ncols = x.e.second;
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);
116  break;
117  }
118  case value_type::scalar: {
119  size_type nrow = undef;
120  size_type dis_nrow = undef;
121  size_type ncol = undef;
122  size_type dis_ncol = undef;
123 #ifdef TODO
124  // TODO: non-zero value => force size to one ?
125  if (x.s != 0) {
126  nrow = (my_proc == iproc0) ? 1 : 0;
127  dis_nrow = 1;
128  ncol = (my_proc == iproc0) ? 1 : 0;
129  dis_ncol = 1;
130  }
131 #endif // TODO
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);
135  break;
136  }
137  case value_type::vector: {
138  // vertical vec
139  size_type nrow = x.v.ownership().size();
140  size_type dis_nrow = x.v.ownership().dis_size();
141  size_type ncol = (my_proc == iproc0) ? 1 : 0;
142  size_type dis_ncol = 1;
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);
146  comm = x.v.comm(); // TODO: check if compatible comm
147  break;
148  }
149  case value_type::vector_transpose: {
150  // horizontal vec
151  size_type ncol = x.v.ownership().size();
152  size_type dis_ncol = x.v.ownership().dis_size();
153  size_type nrow = (my_proc == iproc0) ? 1 : 0;
154  size_type dis_nrow = 1;
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);
157  comm = x.v.comm(); // TODO: check if compatible comm
158  comm = x.v.comm(); // TODO: check if compatible comm
159  break;
160  }
161  case value_type::vector_vec: {
162  // horizontal vector<vec>
163  size_type ncol = undef;
164  size_type dis_ncol = undef;
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(); // TODO: check if compatible comm
169  }
170  size_type nrow = (my_proc == iproc0) ? x.vv.size() : 0;
171  size_type dis_nrow = x.vv.size();
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);
174  break;
175  }
176  case value_type::vector_vec_transpose: {
177  // vertical vector<vec>
178  size_type nrow = undef;
179  size_type dis_nrow = undef;
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(); // TODO: check if compatible comm
184  }
185  size_type ncol = (my_proc == iproc0) ? x.vv.size() : 0;
186  size_type dis_ncol = x.vv.size();
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);
189  break;
190  }
191  case value_type::matrix: {
192  size_type nrow = x.m.row_ownership().size();
193  size_type dis_nrow = x.m.row_ownership().dis_size();
194  size_type ncol = x.m.col_ownership().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(); // TODO: check if compatible comm
199  trace_macro ("block col j_comp="<<j_comp<<": M: "<<row_sizes.first<<"x"<<col_sizes[j_comp].first);
200  break;
201  break;
202  }
203  }
204  trace_macro ("block col j_comp="<<j_comp<<" done");
205  }
206 }
207 // ---------------------------
208 // pass 1 : compute ownership
209 // ---------------------------
210 template <class T, class M>
211 void
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,
216  distributor& row_ownership,
217  distributor& col_ownership,
218  std::vector<distributor>& col_start_by_component) const
219 {
220  row_ownership = distributor (row_sizes.second, comm, row_sizes.first);
221  size_type col_comp_start_j = 0;
222  size_type col_comp_start_dis_j = 0;
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;
228  }
229  col_ownership = distributor (col_comp_start_dis_j, comm, col_comp_start_j);
230 }
231 // -----------------------------------
232 // pass 2 : copy at the right location
233 // -----------------------------------
234 template <class T, class M>
235 void
237  asr<T,M>& A,
238  const std::pair<size_type,size_type>& row_sizes,
239  const distributor& row_start_by_component,
240  const std::vector<std::pair<size_type,size_type> >& col_sizes,
241  const std::vector<distributor>& col_start_by_component) const
242 {
243  communicator comm = A.row_ownership().comm();
244  size_type my_proc = comm.rank();
245  size_type iproc0 = constraint_process_rank (comm);
246  size_type j_comp = 0;
247  for (typename std::list<value_type>::const_iterator iter = _l.begin(); iter != _l.end(); ++iter, j_comp++) {
248  const value_type& x = *iter;
249  switch (x.variant) {
250  case value_type::empty: {
251  // no value to insert
252  break;
253  }
254  case value_type::scalar: {
255  // do not add eplicit zero values
256  size_type dis_i = A.row_ownership().first_index(iproc0)
257  + row_start_by_component.size(iproc0);
258  size_type dis_j = A.col_ownership().first_index(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;
261  break;
262  }
263  case value_type::vector: {
264  // vertical vec : moved column to iproc0, but still distributed by lines
265  size_type dis_i = A.row_ownership().first_index()
266  + row_start_by_component.size();
267  size_type dis_j = A.col_ownership().first_index(iproc0)
268  + col_start_by_component[j_comp].size(iproc0);
269  for (typename vec<T,M>::const_iterator iter = x.v.begin(), last = x.v.end(); iter != last; iter++, dis_i++) {
270  if (*iter != 0) A.dis_entry (dis_i, dis_j) += *iter;
271  }
272  break;
273  }
274  case value_type::vector_transpose: {
275  // horizontal vec: row integrally hosted by iproc=0 ; col still distributed
276  size_type dis_i = A.row_ownership().first_index(iproc0)
277  + row_start_by_component.size(iproc0);
278  size_type dis_j = A.col_ownership().first_index()
279  + col_start_by_component[j_comp].size();
280  for (typename vec<T,M>::const_iterator iter = x.v.begin(), last = x.v.end(); iter != last; iter++, dis_j++) {
281  if (*iter != 0) A.dis_entry (dis_i, dis_j) += *iter;
282  }
283  break;
284  }
285  case value_type::vector_vec: {
286  // horizontal vector<vec>: rows integrally hosted by iproc=0 ; col still distributed
287  size_type n = x.vv.size();
288  for (size_type k = 0; k < n; k++) {
289  size_type dis_i = A.row_ownership().first_index(iproc0)
290  + row_start_by_component.size(iproc0)
291  + k;
292  size_type dis_j = A.col_ownership().first_index()
293  + col_start_by_component[j_comp].size();
294  for (typename vec<T,M>::const_iterator iter = x.vv[k].begin(), last = x.vv[k].end(); iter != last; iter++, dis_j++) {
295  if (*iter != 0) A.dis_entry (dis_i, dis_j) += *iter;
296  }
297  }
298  break;
299  }
300  case value_type::vector_vec_transpose: {
301  // vertical vector<vec>: moved column to iproc0, but still distributed by lines
302  size_type n = x.vv.size();
303  for (size_type k = 0; k < n; k++) {
304  size_type dis_i = A.row_ownership().first_index()
305  + row_start_by_component.size();
306  size_type dis_j = A.col_ownership().first_index(iproc0)
307  + col_start_by_component[j_comp].size(iproc0)
308  + k;
309  for (typename vec<T,M>::const_iterator iter = x.vv[k].begin(), last = x.vv[k].end(); iter != last; iter++, dis_i++) {
310  if (*iter != 0) A.dis_entry (dis_i, dis_j) += *iter;
311  }
312  }
313  break;
314  }
315  case value_type::matrix: {
316  typename csr<T,M>::const_iterator dia_ia = x.m.begin();
317  typename csr<T,M>::const_iterator ext_ia = x.m.ext_begin();
318  for (size_type i = 0, n = x.m.nrow(); i < n; i++) {
319  size_type dis_i = A.row_ownership().first_index()
320  + row_start_by_component.size()
321  + i;
322  for (typename csr<T,M>::const_data_iterator p = dia_ia[i]; p < dia_ia[i+1]; p++) {
323  size_type j = (*p).first;
324  size_type dis_j = A.col_ownership().first_index()
325  + col_start_by_component[j_comp].size()
326  + j;
327  A.dis_entry (dis_i, dis_j) += (*p).second;
328  }
329  if (is_distributed<M>::value && x.m.ext_nnz() != 0) {
330  for (typename csr<T,M>::const_data_iterator p = ext_ia[i]; p < ext_ia[i+1]; p++) {
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;
336  size_type dis_j = A.col_ownership().first_index(iproc)
337  + col_start_by_component[j_comp].size(iproc)
338  + j;
339  A.dis_entry (dis_i, dis_j) += (*p).second;
340  }
341  }
342  }
343  break;
344  }
345  }
346  }
347 }
348 // ------------------------------
349 // merge passes 1 & 2
350 // ------------------------------
351 template <class T, class M>
352 csr<T,M>
354 {
355  communicator comm;
356  std::pair<size_type,size_type> row_sizes = std::pair<size_type,size_type>(undef,undef);
357  distributor row_start_by_component (0, comm, 0);
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);
363  distributor row_ownership;
364  distributor col_ownership;
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);
368  A.dis_entry_assembly();
369  return csr<T,M>(A);
370 }
371 // =========================================================================
372 // 2nd case : multi-line matrix initializer
373 // A = { {a, u },
374 // {trans(u), 3.2} };
375 // =========================================================================
376 template <class T, class M>
377 csr<T,M>
379 {
380  // ---------------------------
381  // pass 0 : compute sizes
382  // ---------------------------
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;
385  communicator comm;
386  size_type i_comp = 0;
387  for (typename std::list<line_type>::const_iterator iter = _l.begin(); iter != _l.end(); ++iter, i_comp++) {
388  const line_type& line = *iter;
389  trace_macro ("block line i_comp="<<i_comp<<"...");
390  line.build_csr_pass0 (row_sizes[i_comp], col_sizes, comm);
391  trace_macro ("block line i_comp="<<i_comp<<" done");
392  }
393  csr_concat_line<T,M>::finalize_sizes (row_sizes, comm);
394  csr_concat_line<T,M>::finalize_sizes (col_sizes, comm);
395  // ---------------------------
396  // pass 1 : compute ownership
397  // ---------------------------
398  size_type row_comp_start_i = 0;
399  size_type row_comp_start_dis_i = 0;
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;
406  }
407  distributor row_ownership = row_start_by_component[row_sizes.size()];
408 
409  size_type col_comp_start_j = 0;
410  size_type col_comp_start_dis_j = 0;
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;
417  }
418  distributor col_ownership = col_start_by_component[col_sizes.size()];
419  // ------------------------
420  // pass 2 : copy
421  // ------------------------
422  asr<T,M> a (row_ownership, col_ownership);
423  i_comp = 0;
424  for (typename std::list<line_type>::const_iterator iter = _l.begin(); iter != _l.end(); ++iter, ++i_comp) {
425  const line_type& line = *iter;
426  line.build_csr_pass2 (a, row_sizes[i_comp], row_start_by_component[i_comp], col_sizes, col_start_by_component);
427  }
428  a.dis_entry_assembly();
429  return csr<T,M>(a);
430 }
431 // ----------------------------------------------------------------------------
432 // instanciation in library
433 // ----------------------------------------------------------------------------
434 #define _RHEOLEF_instanciation(T,M) \
435 template class csr_concat_line<T,M>; \
436 template class csr_concat<T,M>;
437 
438 _RHEOLEF_instanciation(Float,sequential)
439 #ifdef _RHEOLEF_HAVE_MPI
440 _RHEOLEF_instanciation(Float,distributed)
441 #endif // _RHEOLEF_HAVE_MPI
442 
443 }} // namespace rheolef::details
see the Float page for the full documentation
void dis_entry_assembly()
Definition: asr.h:112
const distributor & col_ownership() const
Definition: asr.h:101
dis_reference dis_entry(size_type dis_i, size_type dis_j)
Definition: asr.h:193
const distributor & row_ownership() const
Definition: asr.h:100
see the csr page for the full documentation
Definition: csr.h:317
csr< T, M > build_csr() const
Definition: csr_concat.cc:353
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
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
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
csr_concat_value< T, M >::size_type size_type
Definition: csr_concat.h:98
csr_concat_value< T, M >::sizes_type sizes_type
Definition: csr_concat.h:99
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
static void finalize_sizes(std::pair< size_type, size_type > &sizes, const communicator &comm)
Definition: csr_concat.cc:64
std::vector< vec< T, M > > vv
Definition: csr_concat.h:87
csr< T, M > build_csr() const
Definition: csr_concat.cc:378
csr_concat_value< T, M >::size_type size_type
Definition: csr_concat.h:180
see the distributor page for the full documentation
Definition: distributor.h:62
size_type size(size_type iproc) const
Definition: distributor.h:163
size_type first_index(size_type iproc) const
global index range and local size owned by ip-th process
Definition: distributor.h:151
const communicator_type & comm() const
Definition: distributor.h:145
see the vec page for the full documentation
Definition: vec.h:79
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
This file is part of Rheolef.
_RHEOLEF_instanciation(Float) _RHEOLEF_instanciation_evaluate(Float
Definition: sphere.icc:25