Rheolef  7.2
an efficient C++ finite element environment
csr_seq.cc
Go to the documentation of this file.
1 
22 #include "rheolef/csr.h"
23 #include "rheolef/asr.h"
24 #include "rheolef/asr_to_csr.h"
25 #include "rheolef/msg_util.h"
26 #include "rheolef/csr_amux.h"
27 #include "rheolef/csr_cumul_trans_mult.h"
28 using namespace std;
29 namespace rheolef {
30 // ----------------------------------------------------------------------------
31 // class member functions
32 // ----------------------------------------------------------------------------
33 
34 template<class T>
35 csr_rep<T,sequential>::csr_rep(const distributor& row_ownership, const distributor& col_ownership, size_type nnz1)
36  : vector_of_iterator<pair_type>(row_ownership.size()+1),
37  _row_ownership (row_ownership),
38  _col_ownership (col_ownership),
39  _data (nnz1),
40  _is_symmetric (false),
41  _is_definite_positive (false),
42  _pattern_dimension (0)
43 {
44 }
45 template<class T>
46 void
47 csr_rep<T,sequential>::resize (const distributor& row_ownership, const distributor& col_ownership, size_type nnz1)
48 {
49  vector_of_iterator<pair_type>::resize (row_ownership.size()+1);
50  _row_ownership = row_ownership;
51  _col_ownership = col_ownership;
52  _data.resize (nnz1);
53  // first pointer points to the beginning of the data:
54  vector_of_iterator<pair_type>::operator[](0) = _data.begin().operator->();
55 }
56 template<class T>
58  : vector_of_iterator<pair_type> (loc_nrow1+1),
59  _row_ownership (distributor::decide, communicator(), loc_nrow1),
60  _col_ownership (distributor::decide, communicator(), loc_ncol1),
61  _data (loc_nnz1),
62  _is_symmetric (false),
63  _is_definite_positive (false),
64  _pattern_dimension (0)
65 {
66 }
67 template<class T>
68 void
70 {
72  _row_ownership = distributor (distributor::decide, communicator(), loc_nrow1);
73  _col_ownership = distributor (distributor::decide, communicator(), loc_ncol1);
74  _data.resize (loc_nnz1);
75  // first pointer points to the beginning of the data:
76  vector_of_iterator<pair_type>::operator[](0) = _data.begin().operator->();
77 }
78 template<class T>
80  : vector_of_iterator<pair_type>(b.nrow()+1),
81  _row_ownership (b.row_ownership()),
82  _col_ownership (b.col_ownership()),
83  _data(b._data),
84  _is_symmetric (b._is_symmetric),
85  _is_definite_positive (b._is_definite_positive),
86  _pattern_dimension (b._pattern_dimension)
87 {
88  // physical copy of csr
89  const_iterator ib = b.begin();
90  iterator ia = begin();
91  ia[0] = _data.begin().operator->();
92  for (size_type i = 0, n = b.nrow(); i < n; i++) {
93  ia [i+1] = ia[0] + (ib[i+1] - ib[0]);
94  }
95 }
96 template<class T>
97 template<class A>
98 void
100 {
101  resize (a.row_ownership(), a.col_ownership(), a.nnz());
102  typedef std::pair<size_type,T> pair_type;
103  typedef std::pair<size_type,T> const_pair_type;
104  //typedef typename asr<T>::row_type::value_type const_pair_type;
105  asr_to_csr (
106  a.begin(),
107  a.end(),
111  _data.begin());
112 }
113 template<class T>
116 {
117  std::ostream& os = ods.os();
118  os << setprecision (std::numeric_limits<T>::digits10)
119  << "%%MatrixMarket matrix coordinate real general" << std::endl
120  << nrow() << " " << ncol() << " " << nnz() << endl;
121  const_iterator ia = begin();
122  const size_type base = 1;
123  for (size_type i = 0, n = nrow(); i < n; i++) {
124  for (const_data_iterator iter = ia[i], last = ia[i+1]; iter != last; ++iter) {
125  os << i+first_dis_i+base << " "
126  << (*iter).first+base << " "
127  << (*iter).second << endl;
128  }
129  }
130  return ods;
131 }
132 template<class T>
135 {
136  std::ostream& os = ods.os();
137  std::string name = iorheo::getbasename(ods.os());
138  if (name == "") name = "a";
139  os << setprecision (std::numeric_limits<T>::digits10)
140  << name << " = spalloc(" << nrow() << "," << ncol() << "," << nnz() << ");" << endl;
141  const_iterator ia = begin();
142  const size_type base = 1;
143  for (size_type i = 0, n = nrow(); i < n; i++) {
144  for (const_data_iterator iter = ia[i], last = ia[i+1]; iter != last; ++iter) {
145  os << name << "(" << i+first_dis_i+base << "," << (*iter).first+base << ") = "
146  << (*iter).second << ";" << endl;
147  }
148  }
149  return ods;
150 }
151 template<class T>
154 {
155  iorheo::flag_type format = iorheo::flags(ods.os()) & iorheo::format_field;
156  if (format [iorheo::matlab] || format [iorheo::sparse_matlab])
157  { return put_sparse_matlab (ods,first_dis_i); }
158  // default is matrix market format
159  return put_matrix_market (ods,first_dis_i);
160 }
161 template<class T>
162 void
163 csr_rep<T,sequential>::dump (const string& name, size_type first_dis_i) const
164 {
165  std::ofstream os (name.c_str());
166  std::cerr << "! file \"" << name << "\" created." << std::endl;
167  odiststream ods(os);
168  put (ods);
169 }
170 // ----------------------------------------------------------------------------
171 // basic linear algebra
172 // ----------------------------------------------------------------------------
173 
174 template<class T>
175 void
177  const vec<T,sequential>& x,
179  const
180 {
181  check_macro (x.size() == ncol(), "csr*vec: incompatible csr("<<nrow()<<","<<ncol()<<") and vec("<<x.size()<<")");
182  csr_amux (
185  x.begin(),
186  details::generic_set_op(),
187  y.begin());
188 }
189 template<class T>
190 void
192  const vec<T,sequential>& x,
194  const
195 {
196  // reset y and then cumul
197  check_macro (x.size() == nrow(), "csr.trans_mult(vec): incompatible csr("<<nrow()<<","<<ncol()<<") and vec("<<x.size()<<")");
198  std::fill (y.begin(), y.end(), T(0));
202  x.begin(),
203  details::generic_set_plus_op(),
204  y.begin());
205 }
206 template<class T>
209 {
210  iterator ia = begin();
211  for (data_iterator p = ia[0], last_p = ia[nrow()]; p != last_p; p++) {
212  (*p).second *= lambda;
213  }
214  return *this;
215 }
216 // ----------------------------------------------------------------------------
217 // expression c=a+b and c=a-b with a temporary c=*this
218 // ----------------------------------------------------------------------------
219 template<class T>
220 template<class BinaryOp>
221 void
223  const csr_rep<T,sequential>& a,
224  const csr_rep<T,sequential>& b,
225  BinaryOp binop)
226 {
227  check_macro (a.nrow() == b.nrow() && a.ncol() == b.ncol(),
228  "incompatible csr add(a,b): a("<<a.nrow()<<":"<<a.ncol()<<") and "
229  "b("<<b.nrow()<<":"<<b.ncol()<<")");
230  //
231  // first pass: compute nnz_c and resize
232  //
233  size_type nnz_c = 0;
234  const size_type infty = std::numeric_limits<size_type>::max();
235  const_iterator ia = a.begin();
236  const_iterator ib = b.begin();
237  for (size_type i = 0, n = a.nrow(); i < n; i++) {
238  for (const_data_iterator iter_jva = ia[i], last_jva = ia[i+1],
239  iter_jvb = ib[i], last_jvb = ib[i+1];
240  iter_jva != last_jva || iter_jvb != last_jvb; ) {
241 
242  size_type ja = iter_jva == last_jva ? infty : (*iter_jva).first;
243  size_type jb = iter_jvb == last_jvb ? infty : (*iter_jvb).first;
244  if (ja == jb) {
245  iter_jva++;
246  iter_jvb++;
247  } else if (ja < jb) {
248  iter_jva++;
249  } else {
250  iter_jvb++;
251  }
252  nnz_c++;
253  }
254  }
255  resize (a.row_ownership(), b.col_ownership(), nnz_c);
256  data_iterator iter_jvc = _data.begin().operator->();
257  iterator ic = begin();
258  *ic++ = iter_jvc;
259  //
260  // second pass: add and store in c
261  //
262  for (size_type i = 0, n = a.nrow(); i < n; i++) {
263  for (const_data_iterator iter_jva = ia[i], last_jva = ia[i+1],
264  iter_jvb = ib[i], last_jvb = ib[i+1];
265  iter_jva != last_jva || iter_jvb != last_jvb; ) {
266 
267  size_type ja = iter_jva == last_jva ? infty : (*iter_jva).first;
268  size_type jb = iter_jvb == last_jvb ? infty : (*iter_jvb).first;
269  if (ja == jb) {
270  *iter_jvc++ = std::pair<size_type,T> (ja, binop((*iter_jva).second, (*iter_jvb).second));
271  iter_jva++;
272  iter_jvb++;
273  } else if (ja < jb) {
274  *iter_jvc++ = std::pair<size_type,T> (ja, binop((*iter_jva).second, T(0)));
275  iter_jva++;
276  } else {
277  *iter_jvc++ = std::pair<size_type,T> (jb, binop(T(0), (*iter_jvb).second));
278  iter_jvb++;
279  }
280  }
281  *ic++ = iter_jvc;
282  }
283  set_symmetry (a.is_symmetric() && b.is_symmetric());
284  set_pattern_dimension (std::max(a.pattern_dimension(), b.pattern_dimension()));
285  set_definite_positive (a.is_definite_positive() || b.is_definite_positive());
286  // perhaps too optimist for definite_positive: should be ajusted before calling a solver
287 }
288 // ----------------------------------------------------------------------------
289 // trans(a)
290 // ----------------------------------------------------------------------------
291 /*@!
292  \vfill \pagebreak \mbox{} \vfill \begin{algorithm}[h]
293  \caption{{\tt sort}: sort rows by increasing column order}
294  \begin{algorithmic}
295  \INPUT {the matrix in CSR format}
296  ia(0:nrow), ja(0:nnz-1), a(0:nnz-1)
297  \ENDINPUT
298  \OUTPUT {the sorted CSR matrix}
299  iao(0:nrow), jao(0:nnzl-1), ao(0:nnzl-1)
300  \ENDOUTPUT
301  \BEGIN
302  {\em first: reset iao} \\
303  \FORTO {i := 0} {nrow}
304  iao(i) := 0;
305  \ENDFOR
306 
307  {\em second: compute lengths from pointers} \\
308  \FORTO {i := 0} {nrow-1}
309  \FORTO {p := ia(i)} {ia(i+1)-1}
310  iao (ja(p)+1)++;
311  \ENDFOR
312  \ENDFOR
313 
314  {\em third: compute pointers from lengths} \\
315  \FORTO {i := 0} {nrow-1}
316  iao(i+1) += iao(i)
317  \ENDFOR
318 
319  {\em fourth: copy values} \\
320  \FORTO {i := 0} {nrow-1}
321  \FORTO {p := ia(i)} {ia(i+1)-1}
322  j := ja(p) \\
323  q := iao(j) \\
324  jao(q) := i \\
325  ao(q) := a(p) \\
326  iao (j)++
327  \ENDFOR
328  \ENDFOR
329 
330  {\em fiveth: reshift pointers} \\
331  \FORTOSTEP {i := nrow-1} {0} {-1}
332  iao(i+1) := iao(i);
333  \ENDFOR
334  iao(0) := 0
335  \END
336  \end{algorithmic} \end{algorithm}
337  \vfill \pagebreak \mbox{} \vfill
338 */
339 
340 // NOTE: transposed matrix has always rows sorted by increasing column indexes
341 // even if original matrix has not
342 template<class T>
343 void
345 {
346  b.resize (col_ownership(), row_ownership(), nnz());
347 
348  // first pass: set ib(*) to ib(0)
349  iterator ib = b.begin();
350  for (size_type j = 0, m = b.nrow(); j < m; j++) {
351  ib[j+1] = ib[0];
352  }
353  // second pass: compute lengths of row(i) of a^T in ib(i+1)-ib(0)
354  const_iterator ia = begin();
355  for (size_type i = 0, n = nrow(); i < n; i++) {
356  for (const_data_iterator p = ia[i], last_p = ia[i+1]; p != last_p; p++) {
357  size_type j = (*p).first;
358  ib [j+1]++;
359  }
360  }
361  // third pass: compute pointers from lengths
362  for (size_type j = 0, m = b.nrow(); j < m; j++) {
363  ib [j+1] += (ib[j]-ib[0]);
364  }
365  // fourth pass: store values
366  data_iterator q0 = ib[0];
367  for (size_type i = 0, n = nrow(); i < n; i++) {
368  for (const_data_iterator p = ia[i], last_p = ia[i+1]; p != last_p; p++) {
369  size_type j = (*p).first;
370  data_iterator q = ib[j];
371  (*q).first = i;
372  (*q).second = (*p).second;
373  ib[j]++;
374  }
375  }
376  // fiveth: shift pointers
377  for (long int j = b.nrow()-1; j >= 0; j--) {
378  ib[j+1] = ib[j];
379  }
380  ib[0] = q0;
381 }
382 // ----------------------------------------------------------------------------
383 // set symmetry by check
384 // ----------------------------------------------------------------------------
385 template<class T>
386 void
388 {
389  if (nrow() != ncol()) {
390  set_symmetry(false);
391  return;
392  }
393  // check if a-trans(a) == 0 up to machine prec
395  build_transpose (at);
396  d.assign_add (*this, at, std::minus<T>());
397  set_symmetry(d.max_abs() <= tol);
398 }
399 // ----------------------------------------------------------------------------
400 // expression c=a*b with a temporary c=*this
401 // ----------------------------------------------------------------------------
402 // compute the sparse size of c=a*b
403 template<class T>
406  const csr_rep<T,sequential>& a,
407  const csr_rep<T,sequential>& b)
408 {
410  typename csr_rep<T,sequential>::const_iterator ia = a.begin(), ib = b.begin();
411  size_type nnz_c = 0;
412  for (size_type i = 0, n = a.nrow(); i < n; i++) {
413  std::set<size_type> y;
414  for (typename csr<T>::const_data_iterator ap = ia[i], last_ap = ia[i+1]; ap != last_ap; ++ap) {
415  size_type j = (*ap).first;
416  typename std::set<size_type>::iterator iter_y = y.begin();
417  for (typename csr<T>::const_data_iterator bp = ib[j], last_bp = ib[j+1]; bp != last_bp; ++bp) {
418  size_type k = (*bp).first;
419  iter_y = y.insert(iter_y, k);
420  }
421  }
422  nnz_c += y.size();
423  }
424  return nnz_c;
425 }
426 // compute the values of the sparse matrix c=a*b
427 template<class T>
428 static
429 void
430 csr_csr_mult (
431  const csr_rep<T,sequential>& a,
432  const csr_rep<T,sequential>& b,
433  csr_rep<T,sequential>& c)
434 {
436  typename csr_rep<T,sequential>::const_iterator ia = a.begin(), ib = b.begin();
437  typename csr_rep<T,sequential>::iterator ic = c.begin();
438  for (size_type i = 0, n = a.nrow(); i < n; i++) {
439  std::map<size_type,T> y;
440  for (typename csr<T>::const_data_iterator ap = ia[i], last_ap = ia[i+1]; ap != last_ap; ++ap) {
441  size_type j = (*ap).first;
442  const T& a_ij = (*ap).second;
443  typename std::map<size_type,T>::iterator prev_y = y.begin();
444  for (typename csr<T>::const_data_iterator bp = ib[j], last_bp = ib[j+1]; bp != last_bp; ++bp) {
445  size_type k = (*bp).first;
446  const T& b_jk = (*bp).second;
447  T value = a_ij*b_jk;
448  typename std::map<size_type,T>::iterator curr_y = y.find(k);
449  if (curr_y == y.end()) {
450  // create a new (i,k,value) entry in matrix C
451  y.insert (prev_y, std::pair<size_type,T>(k, value));
452  } else {
453  // increment an existing (i,k,value) entry in matrix C
454  (*curr_y).second += value;
455  prev_y = curr_y;
456  }
457  }
458  }
459  ic[i+1] = ic[i] + y.size();
460  // copy the y temporary into the i-th row of C:
461  typename std::map<size_type,T>::const_iterator iter_y = y.begin();
462  for (typename csr<T>::data_iterator cp = ic[i], last_cp = ic[i+1]; cp != last_cp; ++cp, ++iter_y) {
463  *cp = *iter_y;
464  }
465  }
466  // propagate flags ; cannot stat yet for symmetry and positive_definite
467  c.set_pattern_dimension (std::max(a.pattern_dimension(), b.pattern_dimension()));
468 }
469 template<class T>
470 void
472  const csr_rep<T,sequential>& a,
473  const csr_rep<T,sequential>& b)
474 {
475  size_type nnz_c = csr_csr_mult_size (a, b);
476  resize (a.row_ownership(), b.col_ownership(), nnz_c);
477  csr_csr_mult (a, b, *this);
478 }
479 // ----------------------------------------------------------------------------
480 // instanciation in library
481 // ----------------------------------------------------------------------------
482 #ifndef _RHEOLEF_USE_HEAP_ALLOCATOR
483 #define _RHEOLEF_instanciate_class(T) \
484 template void csr_rep<T,sequential>::build_from_asr (const asr<T,sequential,std::allocator<T> >&);
485 #else // _RHEOLEF_USE_HEAP_ALLOCATOR
486 #define _RHEOLEF_instanciate_class(T) \
487 template void csr_rep<T,sequential>::build_from_asr (const asr<T,sequential,std::allocator<T> >&); \
488 template void csr_rep<T,sequential>::build_from_asr (const asr<T,sequential,heap_allocator<T> >&);
489 #endif // _RHEOLEF_USE_HEAP_ALLOCATOR
490 
491 #define _RHEOLEF_instanciate(T) \
492 template class csr_rep<T,sequential>; \
493 template void csr_rep<T,sequential>::assign_add ( \
494  const csr_rep<T,sequential>&, \
495  const csr_rep<T,sequential>&, \
496  std::plus<T>); \
497 template void csr_rep<T,sequential>::assign_add ( \
498  const csr_rep<T,sequential>&, \
499  const csr_rep<T,sequential>&, \
500  std::minus<T>); \
501 _RHEOLEF_instanciate_class(T)
502 
504 
505 } // namespace rheolef
field::size_type size_type
Definition: branch.cc:430
see the Float page for the full documentation
vector_of_iterator< pair_type >::const_value_type const_data_iterator
Definition: csr.h:94
vector_of_iterator< pair_type >::iterator iterator
Definition: csr.h:90
std::pair< size_type, T > pair_type
Definition: csr.h:89
std::vector< T >::size_type size_type
Definition: csr.h:86
vector_of_iterator< pair_type >::value_type data_iterator
Definition: csr.h:93
vector_of_iterator< pair_type >::const_iterator const_iterator
Definition: csr.h:91
see the csr page for the full documentation
Definition: csr.h:317
see the distributor page for the full documentation
Definition: distributor.h:69
void resize(size_type dis_size=0, const communicator_type &c=communicator_type(), size_type loc_size=decide)
Definition: distributor.cc:30
size_type size(size_type iproc) const
Definition: distributor.h:170
static const size_type decide
Definition: distributor.h:83
odiststream: see the diststream page for the full documentation
Definition: diststream.h:137
std::ostream & os()
Definition: diststream.h:247
see the vec page for the full documentation
Definition: vec.h:79
const_reference operator[](size_type i) const
size_t size_type
Definition: basis_get.cc:76
rheolef::std value
Expr1::float_type T
Definition: field_expr.h:230
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
verbose clean transpose logscale grid shrink ball stereo iso volume skipvtk deformation fastfieldload lattice reader_on_stdin color format format format format format format format format format format format format format format format matlab
verbose clean transpose logscale grid shrink ball stereo iso volume skipvtk deformation fastfieldload lattice reader_on_stdin color format format format format format format format format format format format format format format format format format format dump
This file is part of Rheolef.
std::enable_if< details::is_rheolef_arithmetic< U >::value,ad3_basic< T > & >::type operator*=(ad3_basic< T > &a, const U &b)
Definition: ad3.h:367
void put(std::ostream &out, std::string name, const tiny_matrix< T > &a)
Definition: tiny_lu.h:155
_RHEOLEF_instanciate(Float, sequential) _RHEOLEF_instanciate(Float
csr_rep< T, sequential >::size_type csr_csr_mult_size(const csr_rep< T, sequential > &a, const csr_rep< T, sequential > &b)
Definition: csr_seq.cc:405
void put_matrix_market(std::ostream &out, const Eigen::SparseMatrix< T, Eigen::RowMajor > &a)
Definition: eigen_util.h:78
void csr_cumul_trans_mult(InputIterator1 ia, InputIterator1 last_ia, InputIterator3 x, SetOperator set_op, RandomAccessMutableIterator y)
OutputPtrIterator asr_to_csr(InputPtrIterator iter_ptr_a, InputPtrIterator last_ptr_a, Predicate pred, Operation op, OutputPtrIterator iter_ptr_b, OutputDataIterator iter_data_b)
Definition: asr_to_csr.h:70
void csr_amux(InputIterator ia, InputIterator last_ia, InputRandomAcessIterator x, SetOperator set_op, OutputIterator y)
Definition: csr_amux.h:77
Definition: sphere.icc:25