rheolef  6.6
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  _is_definite_positive (false),
19  _pattern_dimension (0)
20 {
21 }
22 template<class T>
23 void
24 csr_rep<T,sequential>::resize (const distributor& row_ownership, const distributor& col_ownership, size_type nnz1)
25 {
26  vector_of_iterator<pair_type>::resize (row_ownership.size()+1);
27  _row_ownership = row_ownership;
28  _col_ownership = col_ownership;
29  _data.resize (nnz1);
30  vector_of_iterator<pair_type>::operator[](0) = _data.begin().operator->();
31 }
32 template<class T>
34  : vector_of_iterator<pair_type> (loc_nrow1+1),
35  _row_ownership (distributor::decide, communicator(), loc_nrow1),
36  _col_ownership (distributor::decide, communicator(), loc_ncol1),
37  _data (loc_nnz1),
38  _is_symmetric (false),
39  _is_definite_positive (false),
40  _pattern_dimension (0)
41 {
42 }
43 template<class T>
44 void
46 {
48  _row_ownership = distributor (distributor::decide, communicator(), loc_nrow1);
49  _col_ownership = distributor (distributor::decide, communicator(), loc_ncol1);
50  _data.resize (loc_nnz1);
51  vector_of_iterator<pair_type>::operator[](0) = _data.begin().operator->();
52 }
53 template<class T>
55  : vector_of_iterator<pair_type>(b.nrow()+1),
56  _row_ownership (b.row_ownership()),
57  _col_ownership (b.col_ownership()),
58  _data(b._data),
59  _is_symmetric (b._is_symmetric),
60  _pattern_dimension (b._pattern_dimension)
61 {
62  const_iterator ib = b.begin();
63  iterator ia = begin();
64  ia[0] = _data.begin().operator->();
65  for (size_type i = 0, n = b.nrow(); i < n; i++) {
66  ia [i+1] = ia[0] + (ib[i+1] - ib[0]);
67  }
68 }
69 template<class T>
70 template<class A>
71 void
73 {
74  resize (a.row_ownership(), a.col_ownership(), a.nnz());
75  typedef std::pair<size_type,T> pair_type;
76  typedef std::pair<size_type,T> const_pair_type;
77  //typedef typename asr<T>::row_type::value_type const_pair_type;
78  asr_to_csr (
79  a.begin(),
80  a.end(),
84  _data.begin());
85 }
86 template<class T>
89 {
90  std::ostream& os = ods.os();
91  os << setprecision (std::numeric_limits<T>::digits10)
92  << "%%MatrixMarket matrix coordinate real general" << std::endl
93  << nrow() << " " << ncol() << " " << nnz() << endl;
94  const_iterator ia = begin();
95  const size_type base = 1;
96  for (size_type i = 0, n = nrow(); i < n; i++) {
97  for (const_data_iterator iter = ia[i], last = ia[i+1]; iter != last; ++iter) {
98  os << i+first_dis_i+base << " "
99  << (*iter).first+base << " "
100  << (*iter).second << endl;
101  }
102  }
103  return ods;
104 }
105 template<class T>
108 {
109  std::ostream& os = ods.os();
110  std::string name = iorheo::getbasename(ods.os());
111  if (name == "") name = "a";
112  os << setprecision (std::numeric_limits<T>::digits10)
113  << name << " = spalloc(" << nrow() << "," << ncol() << "," << nnz() << ");" << endl;
114  const_iterator ia = begin();
115  const size_type base = 1;
116  for (size_type i = 0, n = nrow(); i < n; i++) {
117  for (const_data_iterator iter = ia[i], last = ia[i+1]; iter != last; ++iter) {
118  os << name << "(" << i+first_dis_i+base << "," << (*iter).first+base << ") = "
119  << (*iter).second << ";" << endl;
120  }
121  }
122  return ods;
123 }
124 template<class T>
127 {
129  if (format [iorheo::matlab] || format [iorheo::sparse_matlab])
130  { return put_sparse_matlab (ods,first_dis_i); }
131  return put_matrix_market (ods,first_dis_i);
132 }
133 template<class T>
134 void
135 csr_rep<T,sequential>::dump (const string& name, size_type first_dis_i) const
136 {
137  std::ofstream os (name.c_str());
138  std::cerr << "! file \"" << name << "\" created." << std::endl;
139  odiststream ods(os);
140  put (ods);
141 }
142 
143 template<class T>
144 void
146  const vec<T,sequential>& x,
148  const
149 {
150  csr_amux (
153  x.begin(),
154  set_op<T,T>(),
155  y.begin());
156 }
157 template<class T>
158 void
160  const vec<T,sequential>& x,
162  const
163 {
164  std::fill (y.begin(), y.end(), T(0));
168  x.begin(),
169  set_add_op<T,T>(),
170  y.begin());
171 }
172 template<class T>
175 {
176  iterator ia = begin();
177  for (data_iterator p = ia[0], last_p = ia[nrow()]; p != last_p; p++) {
178  (*p).second *= lambda;
179  }
180  return *this;
181 }
182 template<class T>
183 template<class BinaryOp>
184 void
186  const csr_rep<T,sequential>& a,
187  const csr_rep<T,sequential>& b,
188  BinaryOp binop)
189 {
190  check_macro (a.nrow() == b.nrow() && a.ncol() == b.ncol(),
191  "incompatible csr add(a,b): a("<<a.nrow()<<":"<<a.ncol()<<") and "
192  "b("<<b.nrow()<<":"<<b.ncol()<<")");
193  size_type nnz_c = 0;
195  const_iterator ia = a.begin();
196  const_iterator ib = b.begin();
197  for (size_type i = 0, n = a.nrow(); i < n; i++) {
198  for (const_data_iterator iter_jva = ia[i], last_jva = ia[i+1],
199  iter_jvb = ib[i], last_jvb = ib[i+1];
200  iter_jva != last_jva || iter_jvb != last_jvb; ) {
201 
202  size_type ja = iter_jva == last_jva ? infty : (*iter_jva).first;
203  size_type jb = iter_jvb == last_jvb ? infty : (*iter_jvb).first;
204  if (ja == jb) {
205  iter_jva++;
206  iter_jvb++;
207  } else if (ja < jb) {
208  iter_jva++;
209  } else {
210  iter_jvb++;
211  }
212  nnz_c++;
213  }
214  }
215  resize (a.row_ownership(), b.col_ownership(), nnz_c);
216  data_iterator iter_jvc = _data.begin().operator->();
217  iterator ic = begin();
218  *ic++ = iter_jvc;
219  for (size_type i = 0, n = a.nrow(); i < n; i++) {
220  for (const_data_iterator iter_jva = ia[i], last_jva = ia[i+1],
221  iter_jvb = ib[i], last_jvb = ib[i+1];
222  iter_jva != last_jva || iter_jvb != last_jvb; ) {
223 
224  size_type ja = iter_jva == last_jva ? infty : (*iter_jva).first;
225  size_type jb = iter_jvb == last_jvb ? infty : (*iter_jvb).first;
226  if (ja == jb) {
227  *iter_jvc++ = std::pair<size_type,T> (ja, binop((*iter_jva).second, (*iter_jvb).second));
228  iter_jva++;
229  iter_jvb++;
230  } else if (ja < jb) {
231  *iter_jvc++ = std::pair<size_type,T> (ja, binop((*iter_jva).second, T(0)));
232  iter_jva++;
233  } else {
234  *iter_jvc++ = std::pair<size_type,T> (jb, binop(T(0), (*iter_jvb).second));
235  iter_jvb++;
236  }
237  }
238  *ic++ = iter_jvc;
239  }
240  set_symmetry (a.is_symmetric() && b.is_symmetric());
241  set_pattern_dimension (std::max(a.pattern_dimension(), b.pattern_dimension()));
242 }
243 /*@!
244  \vfill \pagebreak \mbox{} \vfill \begin{algorithm}[h]
245  \caption{{\tt sort}: sort rows by increasing column order}
246  \begin{algorithmic}
247  \INPUT {the matrix in CSR format}
248  ia(0:nrow), ja(0:nnz-1), a(0:nnz-1)
249  \ENDINPUT
250  \OUTPUT {the sorted CSR matrix}
251  iao(0:nrow), jao(0:nnzl-1), ao(0:nnzl-1)
252  \ENDOUTPUT
253  \BEGIN
254  {\em first: reset iao} \\
255  \FORTO {i := 0} {nrow}
256  iao(i) := 0;
257  \ENDFOR
258 
259  {\em second: compute lengths from pointers} \\
260  \FORTO {i := 0} {nrow-1}
261  \FORTO {p := ia(i)} {ia(i+1)-1}
262  iao (ja(p)+1)++;
263  \ENDFOR
264  \ENDFOR
265 
266  {\em third: compute pointers from lengths} \\
267  \FORTO {i := 0} {nrow-1}
268  iao(i+1) += iao(i)
269  \ENDFOR
270 
271  {\em fourth: copy values} \\
272  \FORTO {i := 0} {nrow-1}
273  \FORTO {p := ia(i)} {ia(i+1)-1}
274  j := ja(p) \\
275  q := iao(j) \\
276  jao(q) := i \\
277  ao(q) := a(p) \\
278  iao (j)++
279  \ENDFOR
280  \ENDFOR
281 
282  {\em fiveth: reshift pointers} \\
283  \FORTOSTEP {i := nrow-1} {0} {-1}
284  iao(i+1) := iao(i);
285  \ENDFOR
286  iao(0) := 0
287  \END
288  \end{algorithmic} \end{algorithm}
289  \vfill \pagebreak \mbox{} \vfill
290 */
291 
292 template<class T>
293 void
295 {
296  b.resize (col_ownership(), row_ownership(), nnz());
297 
298  iterator ib = b.begin();
299  for (size_type j = 0, m = b.nrow(); j < m; j++) {
300  ib[j+1] = ib[0];
301  }
302  const_iterator ia = begin();
303  for (size_type i = 0, n = nrow(); i < n; i++) {
304  for (const_data_iterator p = ia[i], last_p = ia[i+1]; p != last_p; p++) {
305  size_type j = (*p).first;
306  ib [j+1]++;
307  }
308  }
309  for (size_type j = 0, m = b.nrow(); j < m; j++) {
310  ib [j+1] += (ib[j]-ib[0]);
311  }
312  data_iterator q0 = ib[0];
313  for (size_type i = 0, n = nrow(); i < n; i++) {
314  for (const_data_iterator p = ia[i], last_p = ia[i+1]; p != last_p; p++) {
315  size_type j = (*p).first;
316  data_iterator q = ib[j];
317  (*q).first = i;
318  (*q).second = (*p).second;
319  ib[j]++;
320  }
321  }
322  for (long int j = b.nrow()-1; j >= 0; j--) {
323  ib[j+1] = ib[j];
324  }
325  ib[0] = q0;
326 }
327 template<class T>
328 void
330 {
331  if (nrow() != ncol()) {
332  set_symmetry(false);
333  return;
334  }
336  build_transpose (at);
337  d.assign_add (*this, at, std::minus<T>());
338  set_symmetry(d.max_abs() <= tol);
339 }
340 template<class T>
343  const csr_rep<T,sequential>& a,
344  const csr_rep<T,sequential>& b)
345 {
347  typename csr_rep<T,sequential>::const_iterator ia = a.begin(), ib = b.begin();
348  size_type nnz_c = 0;
349  for (size_type i = 0, n = a.nrow(); i < n; i++) {
350  std::set<size_type> y;
351  for (typename csr<T>::const_data_iterator ap = ia[i], last_ap = ia[i+1]; ap != last_ap; ++ap) {
352  size_type j = (*ap).first;
353  typename std::set<size_type>::iterator iter_y = y.begin();
354  for (typename csr<T>::const_data_iterator bp = ib[j], last_bp = ib[j+1]; bp != last_bp; ++bp) {
355  size_type k = (*bp).first;
356  iter_y = y.insert(iter_y, k);
357  }
358  }
359  nnz_c += y.size();
360  }
361  return nnz_c;
362 }
363 template<class T>
364 static
365 void
367  const csr_rep<T,sequential>& a,
368  const csr_rep<T,sequential>& b,
370 {
372  typename csr_rep<T,sequential>::const_iterator ia = a.begin(), ib = b.begin();
373  typename csr_rep<T,sequential>::iterator ic = c.begin();
374  for (size_type i = 0, n = a.nrow(); i < n; i++) {
375  std::map<size_type,T> y;
376  for (typename csr<T>::const_data_iterator ap = ia[i], last_ap = ia[i+1]; ap != last_ap; ++ap) {
377  size_type j = (*ap).first;
378  const T& a_ij = (*ap).second;
379  typename std::map<size_type,T>::iterator prev_y = y.begin();
380  for (typename csr<T>::const_data_iterator bp = ib[j], last_bp = ib[j+1]; bp != last_bp; ++bp) {
381  size_type k = (*bp).first;
382  const T& b_jk = (*bp).second;
383  T value = a_ij*b_jk;
384  typename std::map<size_type,T>::iterator curr_y = y.find(k);
385  if (curr_y == y.end()) {
386  y.insert (prev_y, std::pair<size_type,T>(k, value));
387  } else {
388  (*curr_y).second += value;
389  prev_y = curr_y;
390  }
391  }
392  }
393  ic[i+1] = ic[i] + y.size();
394  typename std::map<size_type,T>::const_iterator iter_y = y.begin();
395  for (typename csr<T>::data_iterator cp = ic[i], last_cp = ic[i+1]; cp != last_cp; ++cp, ++iter_y) {
396  *cp = *iter_y;
397  }
398  }
399 }
400 template<class T>
401 void
403  const csr_rep<T,sequential>& a,
404  const csr_rep<T,sequential>& b)
405 {
406  size_type nnz_c = csr_csr_mult_size (a, b);
407  resize (a.row_ownership(), b.col_ownership(), nnz_c);
408  csr_csr_mult (a, b, *this);
409 }
410 #define _RHEOLEF_instanciate(T) \
411 template class csr_rep<T,sequential>; \
412 template void csr_rep<T,sequential>::assign_add ( \
413  const csr_rep<T,sequential>&, \
414  const csr_rep<T,sequential>&, \
415  std::plus<T>); \
416 template void csr_rep<T,sequential>::assign_add ( \
417  const csr_rep<T,sequential>&, \
418  const csr_rep<T,sequential>&, \
419  std::minus<T>); \
420 template void csr_rep<T,sequential>::build_from_asr (const asr<T,sequential,std::allocator<T> >&); \
421 template void csr_rep<T,sequential>::build_from_asr (const asr<T,sequential,heap_allocator<T> >&);
422 
424 
425 } // namespace rheolef
distributor - data distribution table
Definition: distributor.h:19
const distributor & col_ownership() const
Definition: asr.h:76
vec< T, M > & operator*=(vec< T, M > &l, const int &r)
size_type size(size_type ip) const
Definition: distributor.h:110
idiststream, odiststream - distributed interface for large data streams
Definition: diststream.h:68
vector_of_iterator< pair_type >::const_iterator const_iterator
Definition: csr.h:24
const distributor & row_ownership() const
Definition: csr.h:38
void put(std::ostream &out, std::string name, const tiny_matrix< T > &a)
Definition: tiny_lu.h:106
void csr_cumul_trans_mult(InputIterator1 ia, InputIterator1 last_ia, InputIterator3 x, SetOperator set_op, RandomAccessMutableIterator y)
STL namespace.
irheostream, orheostream - large data streams
Definition: compiler.h:7
const_reference operator[](size_type i) const
std::bitset< last > flag_type
Definition: iorheo.h:329
vector_of_iterator< pair_type >::const_value_type const_data_iterator
Definition: csr.h:27
vector_of_iterator< pair_type >::value_type data_iterator
Definition: csr.h:26
void resize(size_type dis_size=0, const communicator_type &c=communicator_type(), size_type loc_size=decide)
Definition: distributor.cc:7
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:342
#define check_macro(ok_condition, message)
Definition: compiler.h:104
double Float
Definition: compiler.h:177
size_t d
size_type nrow() const
Definition: csr.h:44
size_type nnz() const
Definition: asr.h:70
const_iterator begin() const
Definition: csr.h:40
static const size_type decide
Definition: distributor.h:31
static void csr_csr_mult(const csr_rep< T, sequential > &a, const csr_rep< T, sequential > &b, csr_rep< T, sequential > &c)
Definition: csr_seq.cc:366
size_type ncol() const
Definition: csr.h:45
const distributor & row_ownership() const
Definition: asr.h:75
flag_type flags() const
Definition: iorheo.cc:147
std::vector< T >::size_type size_type
Definition: csr.h:19
void assign_add(const csr_rep< T, sequential > &a, const csr_rep< T, sequential > &b, BinaryOp binop)
Definition: csr_seq.cc:185
static flag_type format_field
Definition: iorheo.h:334
boost::proto::result_of::make_expr< boost::proto::tag::function, vec_domain, const vec_detail::max_, const int &, const vec< T, M > & >::type max(const int &x, const vec< T, M > &y)
vec - vector in distributed environment
Definition: vec.h:44
OutputPtrIterator asr_to_csr(InputPtrIterator iter_ptr_a, InputPtrIterator last_ptr_a, Predicate pred, Operation op, OutputPtrIterator iter_ptr_b, OutputDataIterator iter_data_b)
asr_to_csr – sequential sparse matrix convertion
Definition: asr_to_csr.h:51
size_type pattern_dimension() const
Definition: csr.h:56
const distributor & col_ownership() const
Definition: csr.h:39
bool is_symmetric() const
Definition: csr.h:51
reference_element_H::size_type size_type
csr - compressed sparse row matrix
Definition: asr.h:8
vector_of_iterator< pair_type >::iterator iterator
Definition: csr.h:23
void dump(size_t order, int orient, size_t shift)
asr - associative sparse matrix
Definition: asr.h:26
void resize(size_type loc_nrow1=0, size_type loc_ncol1=0, size_type loc_nnz1=0)
Definition: csr_seq.cc:45
std::ostream & os()
Definition: diststream.h:166
std::pair< size_type, T > pair_type
Definition: csr.h:22
#define _RHEOLEF_instanciate(T)
Definition: csr_seq.cc:410