Rheolef  7.1
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  set_op<T,T>(),
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  set_add_op<T,T>(),
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 }
467 template<class T>
468 void
470  const csr_rep<T,sequential>& a,
471  const csr_rep<T,sequential>& b)
472 {
473  size_type nnz_c = csr_csr_mult_size (a, b);
474  resize (a.row_ownership(), b.col_ownership(), nnz_c);
475  csr_csr_mult (a, b, *this);
476 }
477 // ----------------------------------------------------------------------------
478 // instanciation in library
479 // ----------------------------------------------------------------------------
480 #ifdef _RHEOLEF_DO_NOT_USE_HEAP_ALLOCATOR
481 #define _RHEOLEF_instanciate_class(T) \
482 template void csr_rep<T,sequential>::build_from_asr (const asr<T,sequential,std::allocator<T> >&);
483 #else // _RHEOLEF_DO_NOT_USE_HEAP_ALLOCATOR
484 #define _RHEOLEF_instanciate_class(T) \
485 template void csr_rep<T,sequential>::build_from_asr (const asr<T,sequential,std::allocator<T> >&); \
486 template void csr_rep<T,sequential>::build_from_asr (const asr<T,sequential,heap_allocator<T> >&);
487 #endif // _RHEOLEF_DO_NOT_USE_HEAP_ALLOCATOR
488 
489 #define _RHEOLEF_instanciate(T) \
490 template class csr_rep<T,sequential>; \
491 template void csr_rep<T,sequential>::assign_add ( \
492  const csr_rep<T,sequential>&, \
493  const csr_rep<T,sequential>&, \
494  std::plus<T>); \
495 template void csr_rep<T,sequential>::assign_add ( \
496  const csr_rep<T,sequential>&, \
497  const csr_rep<T,sequential>&, \
498  std::minus<T>); \
499 _RHEOLEF_instanciate_class(T)
500 
502 
503 } // namespace rheolef
put
void put(idiststream &in, odiststream &out, bool do_proj, bool do_lumped_mass, bool def_fill_opt, size_type extract_id, const Float &scale_value, const std::pair< Float, Float > &u_range, render_type render)
Definition: branch.cc:500
rheolef::pair_identity
Definition: msg_util.h:117
rheolef::csr_rep
Definition: asr.h:32
rheolef::csr_rep< T, sequential >::const_iterator
vector_of_iterator< pair_type >::const_iterator const_iterator
Definition: csr.h:91
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
rheolef::put
void put(std::ostream &out, std::string name, const tiny_matrix< T > &a)
Definition: tiny_lu.h:155
rheolef::put_matrix_market
void put_matrix_market(std::ostream &out, const Eigen::SparseMatrix< T, Eigen::RowMajor > &a)
Definition: eigen_util.h:79
rheolef::csr_rep< T, sequential >::const_data_iterator
vector_of_iterator< pair_type >::const_value_type const_data_iterator
Definition: csr.h:94
rheolef::csr_rep< T, sequential >::size_type
std::vector< T >::size_type size_type
Definition: csr.h:86
dump
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
Definition: iorheo-members.h:119
rheolef::csr_amux
void csr_amux(InputIterator ia, InputIterator last_ia, InputRandomAcessIterator x, SetOperator set_op, OutputIterator y)
Definition: csr_amux.h:77
rheolef::value
rheolef::std value
rheolef::csr_cumul_trans_mult
void csr_cumul_trans_mult(InputIterator1 ia, InputIterator1 last_ia, InputIterator3 x, SetOperator set_op, RandomAccessMutableIterator y)
Definition: csr_cumul_trans_mult.h:67
rheolef::csr_rep< T, sequential >::data_iterator
vector_of_iterator< pair_type >::value_type data_iterator
Definition: csr.h:93
rheolef::distributor
see the distributor page for the full documentation
Definition: distributor.h:62
mkgeo_ball.c
c
Definition: mkgeo_ball.sh:153
rheolef::odiststream::os
std::ostream & os()
Definition: diststream.h:236
rheolef::distributor::decide
static const size_type decide
Definition: distributor.h:76
rheolef::size_type
size_t size_type
Definition: basis_get.cc:76
rheolef::vec
see the vec page for the full documentation
Definition: vec.h:79
rheolef::csr_rep< T, sequential >::pair_type
std::pair< size_type, T > pair_type
Definition: csr.h:89
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 matlab
Definition: iorheo-members.h:109
p
Definition: sphere.icc:25
rheolef::csr_csr_mult_size
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
rheolef::always_true
Definition: msg_util.h:51
rheolef::vector_of_iterator
Definition: vector_of_iterator.h:33
rheolef::csr_rep< T, sequential >
Definition: csr.h:84
a
Definition: diffusion_isotropic.h:25
rheolef::asr
Definition: asr.h:49
rheolef::csr
see the csr page for the full documentation
Definition: asr.h:31
mkgeo_sector.m
m
Definition: mkgeo_sector.sh:118
rheolef::vector_of_iterator::begin
iterator begin()
Definition: vector_of_iterator.h:110
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::odiststream
odiststream: see the diststream page for the full documentation
Definition: diststream.h:126
Float
see the Float page for the full documentation
mkgeo_ball.d
d
Definition: mkgeo_ball.sh:154
rheolef::csr_rep< T, sequential >::iterator
vector_of_iterator< pair_type >::iterator iterator
Definition: csr.h:90
mkgeo_ball.b
b
Definition: mkgeo_ball.sh:152
mkgeo_ball.n
n
Definition: mkgeo_ball.sh:150
rheolef::vector_of_iterator::operator[]
const_reference operator[](size_type i) const
Definition: vector_of_iterator.h:169
size_type
field::size_type size_type
Definition: branch.cc:425
rheolef::set_op
Definition: msg_util.h:56
rheolef::_RHEOLEF_instanciate
_RHEOLEF_instanciate(Float, sequential) _RHEOLEF_instanciate(Float
rheolef::set_add_op
Definition: msg_util.h:61
rheolef::asr_to_csr
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
rheolef::distributor::resize
void resize(size_type dis_size=0, const communicator_type &c=communicator_type(), size_type loc_size=decide)
Definition: distributor.cc:30
rheolef::operator*=
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
mkgeo_contraction.name
name
Definition: mkgeo_contraction.sh:133
rheolef::distributor::size
size_type size(size_type iproc) const
Definition: distributor.h:163
T
Expr1::float_type T
Definition: field_expr.h:218
lambda
Definition: yield_slip_circle.h:34