rheolef  6.5
csr_seq.cc
Go to the documentation of this file.
1 
2 #include "rheolef/csr.h"
3 #include "rheolef/asr.h"
4 #include "rheolef/asr_to_csr.h"
5 #include "rheolef/msg_util.h"
6 #include "rheolef/csr_amux.h"
7 #include "rheolef/csr_cumul_trans_mult.h"
8 using namespace std;
9 namespace rheolef {
10 
11 template<class T>
12 csr_rep<T,sequential>::csr_rep(const distributor& row_ownership, const distributor& col_ownership, size_type nnz1)
13  : vector_of_iterator<pair_type>(row_ownership.size()+1),
14  _row_ownership (row_ownership),
15  _col_ownership (col_ownership),
16  _data (nnz1),
17  _is_symmetric (false),
18  _pattern_dimension (0)
19 {
20 }
21 template<class T>
22 void
23 csr_rep<T,sequential>::resize (const distributor& row_ownership, const distributor& col_ownership, size_type nnz1)
24 {
25  vector_of_iterator<pair_type>::resize (row_ownership.size()+1);
26  _row_ownership = row_ownership;
27  _col_ownership = col_ownership;
28  _data.resize (nnz1);
29  vector_of_iterator<pair_type>::operator[](0) = _data.begin().operator->();
30 }
31 template<class T>
33  : vector_of_iterator<pair_type> (loc_nrow1+1),
34  _row_ownership (distributor::decide, communicator(), loc_nrow1),
35  _col_ownership (distributor::decide, communicator(), loc_ncol1),
36  _data (loc_nnz1),
37  _is_symmetric (false),
38  _pattern_dimension (0)
39 {
40 }
41 template<class T>
42 void
44 {
46  _row_ownership = distributor (distributor::decide, communicator(), loc_nrow1);
47  _col_ownership = distributor (distributor::decide, communicator(), loc_ncol1);
48  _data.resize (loc_nnz1);
49  vector_of_iterator<pair_type>::operator[](0) = _data.begin().operator->();
50 }
51 template<class T>
53  : vector_of_iterator<pair_type>(b.nrow()+1),
54  _row_ownership (b.row_ownership()),
55  _col_ownership (b.col_ownership()),
56  _data(b._data),
57  _is_symmetric (b._is_symmetric),
58  _pattern_dimension (b._pattern_dimension)
59 {
60  const_iterator ib = b.begin();
61  iterator ia = begin();
62  ia[0] = _data.begin().operator->();
63  for (size_type i = 0, n = b.nrow(); i < n; i++) {
64  ia [i+1] = ia[0] + (ib[i+1] - ib[0]);
65  }
66 }
67 template<class T>
68 template<class A>
69 void
71 {
72  resize (a.row_ownership(), a.col_ownership(), a.nnz());
73  typedef std::pair<size_type,T> pair_type;
74  typedef std::pair<size_type,T> const_pair_type;
75  //typedef typename asr<T>::row_type::value_type const_pair_type;
76  asr_to_csr (
77  a.begin(),
78  a.end(),
82  _data.begin());
83 }
84 template<class T>
87 {
88  std::ostream& os = ods.os();
89  os << setprecision (std::numeric_limits<T>::digits10)
90  << "%%MatrixMarket matrix coordinate real general" << std::endl
91  << nrow() << " " << ncol() << " " << nnz() << endl;
92  const_iterator ia = begin();
93  const size_type base = 1;
94  for (size_type i = 0, n = nrow(); i < n; i++) {
95  for (const_data_iterator iter = ia[i], last = ia[i+1]; iter != last; ++iter) {
96  os << i+first_dis_i+base << " "
97  << (*iter).first+base << " "
98  << (*iter).second << endl;
99  }
100  }
101  return ods;
102 }
103 template<class T>
106 {
107  std::ostream& os = ods.os();
108  std::string name = iorheo::getbasename(ods.os());
109  if (name == "") name = "a";
110  os << setprecision (std::numeric_limits<T>::digits10)
111  << name << " = spalloc(" << nrow() << "," << ncol() << "," << nnz() << ");" << endl;
112  const_iterator ia = begin();
113  const size_type base = 1;
114  for (size_type i = 0, n = nrow(); i < n; i++) {
115  for (const_data_iterator iter = ia[i], last = ia[i+1]; iter != last; ++iter) {
116  os << name << "(" << i+first_dis_i+base << "," << (*iter).first+base << ") = "
117  << (*iter).second << ";" << endl;
118  }
119  }
120  return ods;
121 }
122 template<class T>
125 {
127  if (format [iorheo::matlab] || format [iorheo::sparse_matlab])
128  { return put_sparse_matlab (ods,first_dis_i); }
129  return put_matrix_market (ods,first_dis_i);
130 }
131 template<class T>
132 void
133 csr_rep<T,sequential>::dump (const string& name, size_type first_dis_i) const
134 {
135  std::ofstream os (name.c_str());
136  std::cerr << "! file \"" << name << "\" created." << std::endl;
137  odiststream ods(os);
138  put (ods);
139 }
140 
141 template<class T>
142 void
144  const vec<T,sequential>& x,
146  const
147 {
148  csr_amux (
151  x.begin(),
152  set_op<T,T>(),
153  y.begin());
154 }
155 template<class T>
156 void
158  const vec<T,sequential>& x,
160  const
161 {
162  std::fill (y.begin(), y.end(), T(0));
166  x.begin(),
167  set_add_op<T,T>(),
168  y.begin());
169 }
170 template<class T>
173 {
174  iterator ia = begin();
175  for (data_iterator p = ia[0], last_p = ia[nrow()]; p != last_p; p++) {
176  (*p).second *= lambda;
177  }
178  return *this;
179 }
180 template<class T>
181 template<class BinaryOp>
182 void
184  const csr_rep<T,sequential>& a,
185  const csr_rep<T,sequential>& b,
186  BinaryOp binop)
187 {
188  check_macro (a.nrow() == b.nrow() && a.ncol() == b.ncol(),
189  "incompatible csr add(a,b): a("<<a.nrow()<<":"<<a.ncol()<<") and "
190  "b("<<b.nrow()<<":"<<b.ncol()<<")");
191  size_type nnz_c = 0;
193  const_iterator ia = a.begin();
194  const_iterator ib = b.begin();
195  for (size_type i = 0, n = a.nrow(); i < n; i++) {
196  for (const_data_iterator iter_jva = ia[i], last_jva = ia[i+1],
197  iter_jvb = ib[i], last_jvb = ib[i+1];
198  iter_jva != last_jva || iter_jvb != last_jvb; ) {
199 
200  size_type ja = iter_jva == last_jva ? infty : (*iter_jva).first;
201  size_type jb = iter_jvb == last_jvb ? infty : (*iter_jvb).first;
202  if (ja == jb) {
203  iter_jva++;
204  iter_jvb++;
205  } else if (ja < jb) {
206  iter_jva++;
207  } else {
208  iter_jvb++;
209  }
210  nnz_c++;
211  }
212  }
213  resize (a.row_ownership(), b.col_ownership(), nnz_c);
214  data_iterator iter_jvc = _data.begin().operator->();
215  iterator ic = begin();
216  *ic++ = iter_jvc;
217  for (size_type i = 0, n = a.nrow(); i < n; i++) {
218  for (const_data_iterator iter_jva = ia[i], last_jva = ia[i+1],
219  iter_jvb = ib[i], last_jvb = ib[i+1];
220  iter_jva != last_jva || iter_jvb != last_jvb; ) {
221 
222  size_type ja = iter_jva == last_jva ? infty : (*iter_jva).first;
223  size_type jb = iter_jvb == last_jvb ? infty : (*iter_jvb).first;
224  if (ja == jb) {
225  *iter_jvc++ = std::pair<size_type,T> (ja, binop((*iter_jva).second, (*iter_jvb).second));
226  iter_jva++;
227  iter_jvb++;
228  } else if (ja < jb) {
229  *iter_jvc++ = std::pair<size_type,T> (ja, binop((*iter_jva).second, T(0)));
230  iter_jva++;
231  } else {
232  *iter_jvc++ = std::pair<size_type,T> (jb, binop(T(0), (*iter_jvb).second));
233  iter_jvb++;
234  }
235  }
236  *ic++ = iter_jvc;
237  }
238  set_symmetry (a.is_symmetric() && b.is_symmetric());
239  set_pattern_dimension (std::max(a.pattern_dimension(), b.pattern_dimension()));
240 }
241 /*@!
242  \vfill \pagebreak \mbox{} \vfill \begin{algorithm}[h]
243  \caption{{\tt sort}: sort rows by increasing column order}
244  \begin{algorithmic}
245  \INPUT {the matrix in CSR format}
246  ia(0:nrow), ja(0:nnz-1), a(0:nnz-1)
247  \ENDINPUT
248  \OUTPUT {the sorted CSR matrix}
249  iao(0:nrow), jao(0:nnzl-1), ao(0:nnzl-1)
250  \ENDOUTPUT
251  \BEGIN
252  {\em first: reset iao} \\
253  \FORTO {i := 0} {nrow}
254  iao(i) := 0;
255  \ENDFOR
256 
257  {\em second: compute lengths from pointers} \\
258  \FORTO {i := 0} {nrow-1}
259  \FORTO {p := ia(i)} {ia(i+1)-1}
260  iao (ja(p)+1)++;
261  \ENDFOR
262  \ENDFOR
263 
264  {\em third: compute pointers from lengths} \\
265  \FORTO {i := 0} {nrow-1}
266  iao(i+1) += iao(i)
267  \ENDFOR
268 
269  {\em fourth: copy values} \\
270  \FORTO {i := 0} {nrow-1}
271  \FORTO {p := ia(i)} {ia(i+1)-1}
272  j := ja(p) \\
273  q := iao(j) \\
274  jao(q) := i \\
275  ao(q) := a(p) \\
276  iao (j)++
277  \ENDFOR
278  \ENDFOR
279 
280  {\em fiveth: reshift pointers} \\
281  \FORTOSTEP {i := nrow-1} {0} {-1}
282  iao(i+1) := iao(i);
283  \ENDFOR
284  iao(0) := 0
285  \END
286  \end{algorithmic} \end{algorithm}
287  \vfill \pagebreak \mbox{} \vfill
288 */
289 
290 template<class T>
291 void
293 {
294  b.resize (col_ownership(), row_ownership(), nnz());
295 
296  iterator ib = b.begin();
297  for (size_type j = 0, m = b.nrow(); j < m; j++) {
298  ib[j+1] = ib[0];
299  }
300  const_iterator ia = begin();
301  for (size_type i = 0, n = nrow(); i < n; i++) {
302  for (const_data_iterator p = ia[i], last_p = ia[i+1]; p != last_p; p++) {
303  size_type j = (*p).first;
304  ib [j+1]++;
305  }
306  }
307  for (size_type j = 0, m = b.nrow(); j < m; j++) {
308  ib [j+1] += (ib[j]-ib[0]);
309  }
310  data_iterator q0 = ib[0];
311  for (size_type i = 0, n = nrow(); i < n; i++) {
312  for (const_data_iterator p = ia[i], last_p = ia[i+1]; p != last_p; p++) {
313  size_type j = (*p).first;
314  data_iterator q = ib[j];
315  (*q).first = i;
316  (*q).second = (*p).second;
317  ib[j]++;
318  }
319  }
320  for (long int j = b.nrow()-1; j >= 0; j--) {
321  ib[j+1] = ib[j];
322  }
323  ib[0] = q0;
324 }
325 template<class T>
326 void
328 {
329  if (nrow() != ncol()) {
330  set_symmetry(false);
331  return;
332  }
334  build_transpose (at);
335  d.assign_add (*this, at, std::minus<T>());
336  set_symmetry(d.max_abs() <= tol);
337 }
338 template<class T>
341  const csr_rep<T,sequential>& a,
342  const csr_rep<T,sequential>& b)
343 {
345  typename csr_rep<T,sequential>::const_iterator ia = a.begin(), ib = b.begin();
346  size_type nnz_c = 0;
347  for (size_type i = 0, n = a.nrow(); i < n; i++) {
348  std::set<size_type> y;
349  for (typename csr<T>::const_data_iterator ap = ia[i], last_ap = ia[i+1]; ap != last_ap; ++ap) {
350  size_type j = (*ap).first;
351  typename std::set<size_type>::iterator iter_y = y.begin();
352  for (typename csr<T>::const_data_iterator bp = ib[j], last_bp = ib[j+1]; bp != last_bp; ++bp) {
353  size_type k = (*bp).first;
354  iter_y = y.insert(iter_y, k);
355  }
356  }
357  nnz_c += y.size();
358  }
359  return nnz_c;
360 }
361 template<class T>
362 static
363 void
365  const csr_rep<T,sequential>& a,
366  const csr_rep<T,sequential>& b,
368 {
370  typename csr_rep<T,sequential>::const_iterator ia = a.begin(), ib = b.begin();
371  typename csr_rep<T,sequential>::iterator ic = c.begin();
372  for (size_type i = 0, n = a.nrow(); i < n; i++) {
373  std::map<size_type,T> y;
374  for (typename csr<T>::const_data_iterator ap = ia[i], last_ap = ia[i+1]; ap != last_ap; ++ap) {
375  size_type j = (*ap).first;
376  const T& a_ij = (*ap).second;
377  typename std::map<size_type,T>::iterator prev_y = y.begin();
378  for (typename csr<T>::const_data_iterator bp = ib[j], last_bp = ib[j+1]; bp != last_bp; ++bp) {
379  size_type k = (*bp).first;
380  const T& b_jk = (*bp).second;
381  T value = a_ij*b_jk;
382  typename std::map<size_type,T>::iterator curr_y = y.find(k);
383  if (curr_y == y.end()) {
384  y.insert (prev_y, std::pair<size_type,T>(k, value));
385  } else {
386  (*curr_y).second += value;
387  prev_y = curr_y;
388  }
389  }
390  }
391  ic[i+1] = ic[i] + y.size();
392  typename std::map<size_type,T>::const_iterator iter_y = y.begin();
393  for (typename csr<T>::data_iterator cp = ic[i], last_cp = ic[i+1]; cp != last_cp; ++cp, ++iter_y) {
394  *cp = *iter_y;
395  }
396  }
397 }
398 template<class T>
399 void
401  const csr_rep<T,sequential>& a,
402  const csr_rep<T,sequential>& b)
403 {
404  size_type nnz_c = csr_csr_mult_size (a, b);
405  resize (a.row_ownership(), b.col_ownership(), nnz_c);
406  csr_csr_mult (a, b, *this);
407 }
408 #define _RHEOLEF_instanciate(T) \
409 template class csr_rep<T,sequential>; \
410 template void csr_rep<T,sequential>::assign_add ( \
411  const csr_rep<T,sequential>&, \
412  const csr_rep<T,sequential>&, \
413  std::plus<T>); \
414 template void csr_rep<T,sequential>::assign_add ( \
415  const csr_rep<T,sequential>&, \
416  const csr_rep<T,sequential>&, \
417  std::minus<T>); \
418 template void csr_rep<T,sequential>::build_from_asr (const asr<T,sequential,std::allocator<T> >&); \
419 template void csr_rep<T,sequential>::build_from_asr (const asr<T,sequential,heap_allocator<T> >&);
420 
422 
423 } // namespace rheolef
424