Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
csr_concat.cc
Go to the documentation of this file.
1
21// build csr from initializer list (c++ 2011)
22//
23#include "rheolef/csr_concat.h"
24
25namespace 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// ------------------------
35template <class T, class M>
36void
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}
62template <class T, class M>
63void
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}
76template <class T, class M>
77void
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// ------------------------
89template <class T, class M>
90void
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// ---------------------------
210template <class T, class M>
211void
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// -----------------------------------
234template <class T, class M>
235void
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
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);
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;
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++) {
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++) {
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;
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// ------------------------------
351template <class T, class 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);
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// =========================================================================
376template <class T, class 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) \
435template class csr_concat_line<T,M>; \
436template class csr_concat<T,M>;
437
439#ifdef _RHEOLEF_HAVE_MPI
441#endif // _RHEOLEF_HAVE_MPI
442
443}} // namespace rheolef::details
#define _RHEOLEF_instanciation(T, M, A)
Definition asr.cc:223
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
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
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
csr_concat_value< T, M > value_type
Definition csr_concat.h:105
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)
Definition csr_concat.cc:64
csr< T, M > build_csr() const
csr_concat_value< T, M >::size_type size_type
Definition csr_concat.h:180
csr_concat_line< T, M > line_type
Definition csr_concat.h:187
see the distributor page for the full documentation
Definition distributor.h:69
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
base::const_iterator const_iterator
Definition vec.h:92
#define trace_macro(message)
Definition dis_macros.h:111
#define assert_macro(ok_condition, message)
Definition dis_macros.h:113
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.
Definition sphere.icc:25