← Back to branch summary

~ubuntu-branches/ubuntu/trusty/rheolef/trusty

« back to all changes in this revision

Viewing changes to skit/plib2/csr_concat.h

  • Committer: Package Import Robot
  • Author(s): Pierre Saramito
  • Date: 2012-04-06 09:12:21 UTC
  • mfrom: (1.1.5)
  • Revision ID: package-import@ubuntu.com-20120406091221-m58me99p1nxqui49
Tags: 6.0-1
http://bugs.debian.org/667356
http://bugs.debian.org/661689
* New upstream release 6.0 (major changes):
  - massively distributed and parallel support
  - full FEM characteristic method (Lagrange-Gakerkin method) support
  - enhanced users documentation 
  - source code supports g++-4.7 (closes: #667356)
* debian/control: dependencies for MPI distributed solvers added
* debian/rules: build commands simplified
* debian/librheolef-dev.install: man1/* to man9/* added
* debian/changelog: package description rewritted (closes: #661689)

Show diffs side-by-side

added added

removed removed

Lines of Context:
skit/plib2/csr_concat.h
 
1
#ifndef _RHEOLEF_CSR_CONCAT_H
 
2
#define _RHEOLEF_CSR_CONCAT_H
 
3
///
 
4
/// This file is part of Rheolef.
 
5
///
 
6
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
 
7
///
 
8
/// Rheolef is free software; you can redistribute it and/or modify
 
9
/// it under the terms of the GNU General Public License as published by
 
10
/// the Free Software Foundation; either version 2 of the License, or
 
11
/// (at your option) any later version.
 
12
///
 
13
/// Rheolef is distributed in the hope that it will be useful,
 
14
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
 
15
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
16
/// GNU General Public License for more details.
 
17
///
 
18
/// You should have received a copy of the GNU General Public License
 
19
/// along with Rheolef; if not, write to the Free Software
 
20
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
21
///
 
22
/// =========================================================================
 
23
// build csr from initializer list (c++ 2011)
 
24
//
 
25
#include "rheolef/csr.h"
 
26
 
 
27
namespace rheolef {
 
28
 
 
29
// =========================================================================
 
30
// 1rst case : one-line matrix initializer
 
31
//  A = {a, u};                 // matrix & vector
 
32
//  A = {trans(u), 3.2};        // trans(vector) & scalar
 
33
// =========================================================================
 
34
template <class T, class M>
 
35
struct vec_trans {
 
36
  vec<T,M> v;
 
37
  vec_trans (const vec<T,M>& w) : v(w) {}
 
38
};
 
39
template <class T, class M>
 
40
vec_trans<T,M>
 
41
trans (const vec<T,M>& w) {
 
42
  return vec_trans<T,M> (w);
 
43
}
 
44
 
 
45
template <class T, class M>
 
46
struct csr_concat_value {
 
47
// typedef:
 
48
 typedef enum { scalar, vector, vector_transpose, matrix} variant_type;
 
49
// allocators:
 
50
 csr_concat_value (const T& x)               : s(x), v(),     m(),  variant(scalar) {}
 
51
 csr_concat_value (const vec<T,M>& x)        : s(),  v(x),    m(),  variant(vector) {}
 
52
 csr_concat_value (const vec_trans<T,M>& vt) : s(),  v(vt.v), m(),  variant(vector_transpose) {}
 
53
 csr_concat_value (const csr<T,M>& x)        : s(),  v(),     m(x), variant(matrix) {}
 
54
// io/debug:
 
55
 friend std::ostream& operator<< (std::ostream& o, const csr_concat_value<T,M>& x) {
 
56
  if (x.variant == scalar) return o << "s";
 
57
  else if (x.variant == vector) return o << "v";
 
58
  else if (x.variant == vector_transpose) return o << "vt";
 
59
  else return o << "m";
 
60
 }
 
61
// data:
 
62
public:
 
63
 T            s;
 
64
 vec<T,M>     v;
 
65
 csr<T,M>     m;
 
66
 variant_type variant;
 
67
};
 
68
 
 
69
template <class T, class M>
 
70
struct csr_concat_line {
 
71
// typedef:
 
72
 typedef typename csr<T,M>::size_type   size_type;
 
73
 typedef csr_concat_value<T,M>          value_type;
 
74
// allocators:
 
75
 csr_concat_line () : _l() {}
 
76
 csr_concat_line (const std::initializer_list<value_type>& il) : _l() {
 
77
#ifdef _RHEOLEF_HAVE_STD_INITIALIZER_ITERATOR
 
78
    typedef typename std::initializer_list<value_type>::const_iterator const_iterator;
 
79
#else // _RHEOLEF_HAVE_STD_INITIALIZER_ITERATOR
 
80
    typedef const value_type* const_iterator;
 
81
#endif // _RHEOLEF_HAVE_STD_INITIALIZER_ITERATOR
 
82
    for(const_iterator iter = il.begin(); iter != il.end(); ++iter) {
 
83
        _l.push_back(*iter);
 
84
    }
 
85
 }
 
86
 void push_back (const value_type& x) { _l.push_back (x); }
 
87
 
 
88
 friend std::ostream& operator<< (std::ostream& o, const csr_concat_line<T,M>& x) {
 
89
    std::cout << "{";
 
90
    for(typename std::list<value_type>::const_iterator iter = x._l.begin(); iter != x._l.end(); ++iter) {
 
91
        std::cout << *iter << " ";
 
92
    }
 
93
    return std::cout << "}";
 
94
 }
 
95
// internals:
 
96
 static void set_sizes (
 
97
    std::pair<size_type,size_type>&         row_sizes,
 
98
    std::pair<size_type,size_type>&         col_sizes,
 
99
    const std::pair<size_type,size_type>&   new_row_sizes,
 
100
    const std::pair<size_type,size_type>&   new_col_sizes);
 
101
 
 
102
 static void finalize_sizes (
 
103
    std::pair<size_type,size_type>&         sizes,
 
104
    const communicator&                     comm);
 
105
 
 
106
 static void finalize_sizes (
 
107
    std::vector<std::pair<size_type,size_type> >&  sizes,
 
108
    const communicator&                            comm);
 
109
 
 
110
 void build_csr_pass0 (
 
111
    std::pair<size_type,size_type>&                row_sizes,
 
112
    std::vector<std::pair<size_type,size_type> >&  col_sizes,
 
113
    communicator&                                  comm) const;
 
114
 
 
115
 void build_csr_pass1 (
 
116
    const std::pair<size_type,size_type>&                row_sizes,
 
117
    const std::vector<std::pair<size_type,size_type> >&  col_sizes,
 
118
    const communicator&                                  comm,
 
119
    distributor&                                         row_ownership,
 
120
    distributor&                                         col_ownership,
 
121
    std::vector<distributor>&                            col_start_by_component) const;
 
122
 
 
123
 void build_csr_pass2 (
 
124
          asr<T,M>&                       a,
 
125
          const std::vector<distributor>& col_start_by_component,
 
126
          size_type                       row_start_comp_i = 0,
 
127
          size_type                       row_start_comp_dis_i = 0) const;
 
128
 
 
129
 csr<T,M> build_csr () const;
 
130
 
 
131
// data:
 
132
protected:
 
133
 std::list<value_type> _l;
 
134
};
 
135
// ------------------------
 
136
// sizes utilities
 
137
// ------------------------
 
138
template <class T, class M>
 
139
void
 
140
csr_concat_line<T,M>::set_sizes (
 
141
    std::pair<size_type,size_type>&         row_sizes,
 
142
    std::pair<size_type,size_type>&         col_sizes,
 
143
    const std::pair<size_type,size_type>&   new_row_sizes,
 
144
    const std::pair<size_type,size_type>&   new_col_sizes)
 
145
{
 
146
  size_type undef = std::numeric_limits<size_type>::max();
 
147
  if (row_sizes.first == undef) {
 
148
    row_sizes = new_row_sizes;
 
149
  } else if (new_row_sizes.first != undef) {
 
150
    check_macro (row_sizes.first == new_row_sizes.first,
 
151
            "matrix initializer list: matrix row argument [0:"
 
152
               << new_row_sizes.first << "|" << new_row_sizes.second << "["
 
153
               << " has incomptible size: expect [0:"
 
154
               << row_sizes.first  << "|" << row_sizes.second << "[");
 
155
  }
 
156
  if (col_sizes.first == undef) {
 
157
    col_sizes = new_col_sizes;
 
158
  } else if (new_col_sizes.first != undef) {
 
159
    check_macro (col_sizes.first == new_col_sizes.first,
 
160
            "matrix initializer list: matrix col argument [0:"
 
161
               << new_col_sizes.first << "|" << new_col_sizes.second << "["
 
162
               << " has incomptible size: expect [0:"
 
163
               << col_sizes.first  << "|" << col_sizes.second << "[");
 
164
  }
 
165
}
 
166
template <class T, class M>
 
167
void
 
168
csr_concat_line<T,M>::finalize_sizes (
 
169
  std::pair<size_type,size_type>& sizes,
 
170
  const communicator&             comm)
 
171
{
 
172
  size_type undef = std::numeric_limits<size_type>::max();
 
173
  if (sizes.first == undef) {
 
174
    sizes.first  = (comm.rank() == comm.size() - 1) ? 1 : 0;
 
175
    sizes.second = 1;
 
176
  }
 
177
}
 
178
template <class T, class M>
 
179
void
 
180
csr_concat_line<T,M>::finalize_sizes (
 
181
  std::vector<std::pair<size_type,size_type> >&  sizes,
 
182
  const communicator&                            comm)
 
183
{
 
184
  for (size_type i_comp = 0; i_comp < sizes.size(); i_comp++) {
 
185
    finalize_sizes (sizes [i_comp], comm);
 
186
  }
 
187
}
 
188
// ------------------------
 
189
// pass 0 : compute sizes
 
190
// ------------------------
 
191
template <class T, class M>
 
192
void
 
193
csr_concat_line<T,M>::build_csr_pass0 (
 
194
    std::pair<size_type,size_type>&                row_sizes,
 
195
    std::vector<std::pair<size_type,size_type> >&  col_sizes,
 
196
    communicator&                                  comm) const
 
197
{
 
198
  size_type undef = std::numeric_limits<size_type>::max();
 
199
  if (col_sizes.size() == 0) {
 
200
    col_sizes.resize (_l.size(), std::pair<size_type,size_type>(undef,undef));
 
201
  } else {
 
202
    check_macro (col_sizes.size() == _l.size(),
 
203
            "matrix initializer list: unexpected matrix line size "
 
204
               << _l.size() << ": size " << col_sizes.size()
 
205
               << " was expected");
 
206
  }
 
207
  size_type j_comp = 0;
 
208
  for (typename std::list<value_type>::const_iterator iter = _l.begin(); iter != _l.end(); ++iter, j_comp++) {
 
209
    const value_type& x = *iter;
 
210
    switch (x.variant) {
 
211
      case value_type::matrix: {
 
212
        size_type     nrow = x.m.row_ownership().size();
 
213
        size_type dis_nrow = x.m.row_ownership().dis_size();
 
214
        size_type     ncol = x.m.col_ownership().size();
 
215
        size_type dis_ncol = x.m.col_ownership().dis_size();
 
216
        set_sizes (row_sizes, col_sizes[j_comp], 
 
217
            std::make_pair (nrow, dis_nrow), std::make_pair (ncol, dis_ncol));
 
218
        comm = x.m.row_ownership().comm(); // TODO: check if compatible comm
 
219
        break;
 
220
      } 
 
221
      case value_type::vector: {
 
222
        size_type     nrow = x.v.ownership().size();
 
223
        size_type dis_nrow = x.v.ownership().dis_size();
 
224
        size_type     ncol = (comm.rank() == comm.size() - 1) ? 1 : 0;
 
225
        size_type dis_ncol = 1;
 
226
        set_sizes (row_sizes, col_sizes[j_comp], 
 
227
            std::make_pair (nrow, dis_nrow), std::make_pair (ncol, dis_ncol));
 
228
        comm = x.v.comm(); // TODO: check if compatible comm
 
229
        break;
 
230
      }
 
231
      case value_type::vector_transpose: {
 
232
        size_type     nrow = (comm.rank() == comm.size() - 1) ? 1 : 0;
 
233
        size_type dis_nrow = 1;
 
234
        size_type     ncol = x.v.ownership().size();
 
235
        size_type dis_ncol = x.v.ownership().dis_size();
 
236
        set_sizes (row_sizes, col_sizes[j_comp], 
 
237
            std::make_pair (nrow, dis_nrow), std::make_pair (ncol, dis_ncol));
 
238
        comm = x.v.comm(); // TODO: check if compatible comm
 
239
        break;
 
240
      }
 
241
      case value_type::scalar: {
 
242
        size_type     nrow = undef;
 
243
        size_type dis_nrow = undef;
 
244
        size_type     ncol = undef;
 
245
        size_type dis_ncol = undef;
 
246
#ifdef TODO
 
247
        if (x.s != 0) {
 
248
              nrow = (comm.rank() == comm.size() - 1) ? 1 : 0;
 
249
          dis_nrow = 1;
 
250
              ncol = (comm.rank() == comm.size() - 1) ? 1 : 0;
 
251
          dis_ncol = 1;
 
252
        }
 
253
#endif // TODO
 
254
        set_sizes (row_sizes, col_sizes[j_comp], 
 
255
            std::make_pair (nrow, dis_nrow), std::make_pair (ncol, dis_ncol));
 
256
        break;
 
257
      }
 
258
    }
 
259
  }
 
260
}
 
261
// ---------------------------
 
262
// pass 1 : compute ownership
 
263
// ---------------------------
 
264
template <class T, class M>
 
265
void
 
266
csr_concat_line<T,M>::build_csr_pass1 (
 
267
    const std::pair<size_type,size_type>&                row_sizes,
 
268
    const std::vector<std::pair<size_type,size_type> >&  col_sizes,
 
269
    const communicator&                                  comm,
 
270
    distributor&                                         row_ownership,
 
271
    distributor&                                         col_ownership,
 
272
    std::vector<distributor>&                            col_start_by_component) const
 
273
{
 
274
  row_ownership = distributor (row_sizes.second, comm, row_sizes.first);
 
275
  size_type col_comp_start_j     = 0;
 
276
  size_type col_comp_start_dis_j = 0;
 
277
  col_start_by_component.resize (col_sizes.size());
 
278
  for (size_type j_comp = 0; j_comp < col_sizes.size(); j_comp++) {
 
279
    col_start_by_component [j_comp] = distributor (col_comp_start_dis_j, comm, col_comp_start_j);
 
280
    col_comp_start_j     += col_sizes [j_comp].first;
 
281
    col_comp_start_dis_j += col_sizes [j_comp].second;
 
282
  }
 
283
  col_ownership = distributor (col_comp_start_dis_j, comm, col_comp_start_j);
 
284
}
 
285
// ------------------------
 
286
// pass 2 : copy
 
287
// ------------------------
 
288
template <class T, class M>
 
289
void
 
290
csr_concat_line<T,M>::build_csr_pass2 (
 
291
    asr<T,M>&                       a,
 
292
    const std::vector<distributor>& col_start_by_component,
 
293
    size_type                       row_start_comp_i,
 
294
    size_type                       row_start_comp_dis_i) const
 
295
{
 
296
  communicator comm = a.row_ownership().comm();
 
297
  size_type first_dis_i = a.row_ownership().first_index();
 
298
  size_type first_dis_j = a.col_ownership().first_index();
 
299
  size_type my_proc = comm.rank();
 
300
  size_type last_proc = comm.size() - 1;
 
301
  size_type col_start_comp_dis_j = 0;
 
302
  size_type j_comp = 0;
 
303
  for (typename std::list<value_type>::const_iterator iter = _l.begin(); iter != _l.end(); ++iter, j_comp++) {
 
304
    const value_type& x = *iter;
 
305
    switch (x.variant) {
 
306
      case value_type::matrix: {
 
307
        size_type col_start_comp_j = col_start_by_component[j_comp].size();
 
308
        typename csr<T,M>::const_iterator dia_ia = x.m.begin();
 
309
        typename csr<T,M>::const_iterator ext_ia = x.m.ext_begin();
 
310
        for (size_type i = 0, n = x.m.nrow(); i < n; i++) {
 
311
          size_type dis_i = first_dis_i + row_start_comp_i + i;
 
312
          for (typename csr<T,M>::const_data_iterator p = dia_ia[i]; p < dia_ia[i+1]; p++) {
 
313
            size_type j = (*p).first;
 
314
            size_type dis_j = first_dis_j + col_start_comp_j + j;
 
315
            a.dis_entry (dis_i, dis_j) = (*p).second;
 
316
          }
 
317
          if (is_distributed<M>::value && x.m.ext_nnz() != 0) {
 
318
            for (typename csr<T,M>::const_data_iterator p = ext_ia[i]; p < ext_ia[i+1]; p++) {
 
319
              size_type comp_dis_j = x.m.jext2dis_j ((*p).first);
 
320
              size_type iproc = x.m.col_ownership().find_owner (comp_dis_j);
 
321
              size_type first_comp_dis_j_iproc = x.m.col_ownership().first_index (iproc);
 
322
              assert_macro (comp_dis_j >= first_comp_dis_j_iproc, "invalid index");
 
323
              size_type j = comp_dis_j - first_comp_dis_j_iproc;
 
324
              size_type col_start_comp_j_iproc = col_start_by_component[j_comp].size(iproc);
 
325
              size_type first_dis_j_iproc      = a.col_ownership().first_index(iproc);
 
326
              size_type dis_j = first_dis_j_iproc + col_start_comp_j_iproc + j;
 
327
              a.dis_entry (dis_i, dis_j) = (*p).second;
 
328
            }
 
329
          }
 
330
        }
 
331
        col_start_comp_dis_j += x.m.col_ownership().dis_size();
 
332
        break;
 
333
      } 
 
334
      case value_type::vector: {
 
335
        size_type dis_i = row_start_comp_i       + first_dis_i;
 
336
        size_type dis_j = col_start_comp_dis_j;
 
337
        for (typename vec<T,M>::const_iterator iter = x.v.begin(), last = x.v.end(); iter != last; iter++, dis_i++) {
 
338
          if (*iter != 0) {
 
339
            a.dis_entry (dis_i, dis_j) = *iter;
 
340
          }
 
341
        }
 
342
        col_start_comp_dis_j++;
 
343
        break;
 
344
      }
 
345
      case value_type::vector_transpose: {
 
346
        size_type dis_i = row_start_comp_dis_i + 0;
 
347
        size_type dis_j = col_start_comp_dis_j + x.v.ownership().first_index();
 
348
        for (typename vec<T,M>::const_iterator iter = x.v.begin(), last = x.v.end(); iter != last; iter++, dis_j++) {
 
349
          if (*iter != 0) {
 
350
            a.dis_entry (dis_i, dis_j) = *iter;
 
351
          }
 
352
        }
 
353
        col_start_comp_dis_j += x.v.ownership().dis_size();
 
354
        break;
 
355
      }
 
356
      case value_type::scalar: {
 
357
        if (x.s != 0 && comm.rank() == comm.size() - 1) { 
 
358
          size_type dis_i = row_start_comp_dis_i + 0;
 
359
          size_type dis_j = col_start_comp_dis_j;
 
360
          a.dis_entry (dis_i, dis_j) = x.s;
 
361
        }
 
362
        col_start_comp_dis_j++;
 
363
        break;
 
364
      }
 
365
    }
 
366
  }
 
367
}
 
368
// ------------------------------
 
369
// merge passes 1 & 2
 
370
// ------------------------------
 
371
template <class T, class M>
 
372
csr<T,M>
 
373
csr_concat_line<T,M>::build_csr() const
 
374
{
 
375
  size_type undef = std::numeric_limits<size_type>::max();
 
376
  std::pair<size_type,size_type>               row_sizes = std::pair<size_type,size_type>(undef,undef);
 
377
  std::vector<std::pair<size_type,size_type> > col_sizes;
 
378
  communicator comm;
 
379
  build_csr_pass0 (row_sizes, col_sizes, comm);
 
380
  finalize_sizes (row_sizes, comm);
 
381
  finalize_sizes (col_sizes, comm);
 
382
  distributor              row_ownership;
 
383
  distributor              col_ownership;
 
384
  std::vector<distributor> col_start_by_component;
 
385
  build_csr_pass1 (row_sizes, col_sizes, comm, row_ownership, col_ownership, col_start_by_component);
 
386
  asr<T,M> a (row_ownership, col_ownership);
 
387
  build_csr_pass2 (a, col_start_by_component);
 
388
  a.dis_entry_assembly();
 
389
  return csr<T,M>(a);
 
390
}
 
391
// -------------------------------
 
392
// csr cstor from std::initializer
 
393
// -------------------------------
 
394
#ifdef _RHEOLEF_HAVE_STD_INITIALIZER_LIST
 
395
#define RHEOLEF_csr_cstor(M)                                            \
 
396
template <class T>                                                      \
 
397
inline                                                                  \
 
398
csr<T,M>::csr (const std::initializer_list<csr_concat_value<T,M> >& init_list)  \
 
399
{                                                                       \
 
400
  csr_concat_line<T,M> cc (init_list);                                  \
 
401
  csr<T,M>::operator= (cc.build_csr());                                 \
 
402
}
 
403
RHEOLEF_csr_cstor(sequential)
 
404
#ifdef _RHEOLEF_HAVE_MPI
 
405
RHEOLEF_csr_cstor(distributed)
 
406
#endif // _RHEOLEF_HAVE_MPI
 
407
#undef RHEOLEF_csr_cstor
 
408
#endif // _RHEOLEF_HAVE_STD_INITIALIZER_LIST
 
409
// =========================================================================
 
410
// 2nd case : multi-line matrix initializer
 
411
//  A = { {a,        u  },
 
412
//        {trans(u), 3.2} };
 
413
// =========================================================================
 
414
 
 
415
template <class T, class M>
 
416
struct csr_concat {
 
417
// typedef:
 
418
 typedef typename csr<T,M>::size_type   size_type;
 
419
 typedef csr_concat_line<T,M>           line_type;
 
420
 typedef csr_concat_value<T,M>          value_type;
 
421
// allocators:
 
422
 csr_concat () : _l() {}
 
423
 csr_concat (const std::initializer_list<line_type>& il) : _l() {
 
424
#ifdef _RHEOLEF_HAVE_STD_INITIALIZER_ITERATOR
 
425
    typedef typename std::initializer_list<line_type>::const_iterator const_iterator;
 
426
#else // _RHEOLEF_HAVE_STD_INITIALIZER_ITERATOR
 
427
    typedef const line_type* const_iterator;
 
428
#endif // _RHEOLEF_HAVE_STD_INITIALIZER_ITERATOR
 
429
    for(const_iterator iter = il.begin(); iter != il.end(); ++iter) {
 
430
        _l.push_back(*iter);
 
431
    }
 
432
 }
 
433
 void push_back (const line_type& line) { _l.push_back(line); }
 
434
 
 
435
 friend std::ostream& operator<< (std::ostream& o, const csr_concat<T,M>& x) {
 
436
    std::cout << "{";
 
437
    for(typename std::list<line_type>::const_iterator iter = x._l.begin(); iter != x._l.end(); ++iter) {
 
438
        std::cout << *iter << " ";
 
439
    }
 
440
    return std::cout << "}";
 
441
 }
 
442
// internals:
 
443
 csr<T,M> build_csr () const;
 
444
 
 
445
// data:
 
446
protected:
 
447
 std::list<line_type> _l;
 
448
};
 
449
template <class T, class M>
 
450
csr<T,M>
 
451
csr_concat<T,M>::build_csr() const
 
452
{
 
453
  // ---------------------------
 
454
  // pass 0 : compute sizes
 
455
  // ---------------------------
 
456
  size_type undef = std::numeric_limits<size_type>::max();
 
457
  std::vector<std::pair<size_type,size_type> > row_sizes (_l.size(), std::pair<size_type,size_type>(undef,undef));
 
458
  std::vector<std::pair<size_type,size_type> > col_sizes;
 
459
  communicator comm;
 
460
  size_type i_comp = 0;
 
461
  for (typename std::list<line_type>::const_iterator iter = _l.begin(); iter != _l.end(); ++iter, i_comp++) {
 
462
    const line_type& line = *iter;
 
463
    line.build_csr_pass0 (row_sizes[i_comp], col_sizes, comm);
 
464
  }
 
465
  csr_concat_line<T,M>::finalize_sizes (row_sizes, comm);
 
466
  csr_concat_line<T,M>::finalize_sizes (col_sizes, comm);
 
467
  // ---------------------------
 
468
  // pass 1 : compute ownership
 
469
  // ---------------------------
 
470
  size_type row_comp_start_i     = 0;
 
471
  size_type row_comp_start_dis_i = 0;
 
472
  for (size_type i_comp = 0; i_comp < row_sizes.size(); i_comp++) {
 
473
    row_comp_start_i     += row_sizes [i_comp].first;
 
474
    row_comp_start_dis_i += row_sizes [i_comp].second;
 
475
  }
 
476
  distributor row_ownership = distributor (row_comp_start_dis_i, comm, row_comp_start_i);
 
477
 
 
478
  size_type col_comp_start_j     = 0;
 
479
  size_type col_comp_start_dis_j = 0;
 
480
  std::vector<distributor> col_start_by_component (col_sizes.size());
 
481
  for (size_type j_comp = 0; j_comp < col_sizes.size(); j_comp++) {
 
482
    col_start_by_component [j_comp] = distributor (col_comp_start_dis_j, comm, col_comp_start_j);
 
483
    col_comp_start_j     += col_sizes [j_comp].first;
 
484
    col_comp_start_dis_j += col_sizes [j_comp].second;
 
485
  }
 
486
  distributor col_ownership = distributor (col_comp_start_dis_j, comm, col_comp_start_j);
 
487
  // ------------------------
 
488
  // pass 2 : copy
 
489
  // ------------------------
 
490
  asr<T,M> a (row_ownership, col_ownership);
 
491
  size_type row_start_comp_i = 0;
 
492
  size_type row_start_comp_dis_i = 0;
 
493
  i_comp = 0;
 
494
  for (typename std::list<line_type>::const_iterator iter = _l.begin(); iter != _l.end(); ++iter, ++i_comp) {
 
495
    const line_type& line = *iter;
 
496
    line.build_csr_pass2 (a, col_start_by_component, row_start_comp_i, row_start_comp_dis_i);
 
497
    row_start_comp_i     += row_sizes[i_comp].first;
 
498
    row_start_comp_dis_i += row_sizes[i_comp].second;
 
499
  }
 
500
  a.dis_entry_assembly();
 
501
  return csr<T,M>(a);
 
502
}
 
503
#ifdef _RHEOLEF_HAVE_STD_INITIALIZER_LIST
 
504
#define RHEOLEF_csr_cstor(M)                                            \
 
505
template <class T>                                                      \
 
506
inline                                                                  \
 
507
csr<T,M>::csr (const std::initializer_list<csr_concat_line<T,M> >& init_list) \
 
508
{                                                                       \
 
509
  csr_concat<T,M> cc (init_list);                                       \
 
510
  csr<T,M>::operator= (cc.build_csr());                                 \
 
511
}
 
512
RHEOLEF_csr_cstor(sequential)
 
513
#ifdef _RHEOLEF_HAVE_MPI
 
514
RHEOLEF_csr_cstor(distributed)
 
515
#endif // _RHEOLEF_HAVE_MPI
 
516
#undef RHEOLEF_csr_cstor
 
517
#endif // _RHEOLEF_HAVE_STD_INITIALIZER_LIST
 
518
 
 
519
} // namespace rheolef
 
520
#endif // _RHEOLEF_CSR_CONCAT_H