rheolef  7.0
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  check_macro (x.size() == nrow(), "csr.trans_mult(vec): incompatible csr("<<nrow()<<","<<ncol()<<") and vec("<<x.size()<<")");
165  std::fill (y.begin(), y.end(), T(0));
169  x.begin(),
170  set_add_op<T,T>(),
171  y.begin());
172 }
173 template<class T>
176 {
177  iterator ia = begin();
178  for (data_iterator p = ia[0], last_p = ia[nrow()]; p != last_p; p++) {
179  (*p).second *= lambda;
180  }
181  return *this;
182 }
183 template<class T>
184 template<class BinaryOp>
185 void
187  const csr_rep<T,sequential>& a,
188  const csr_rep<T,sequential>& b,
189  BinaryOp binop)
190 {
191  check_macro (a.nrow() == b.nrow() && a.ncol() == b.ncol(),
192  "incompatible csr add(a,b): a("<<a.nrow()<<":"<<a.ncol()<<") and "
193  "b("<<b.nrow()<<":"<<b.ncol()<<")");
194  size_type nnz_c = 0;
195  const size_type infty = std::numeric_limits<size_type>::max();
196  const_iterator ia = a.begin();
197  const_iterator ib = b.begin();
198  for (size_type i = 0, n = a.nrow(); i < n; i++) {
199  for (const_data_iterator iter_jva = ia[i], last_jva = ia[i+1],
200  iter_jvb = ib[i], last_jvb = ib[i+1];
201  iter_jva != last_jva || iter_jvb != last_jvb; ) {
202 
203  size_type ja = iter_jva == last_jva ? infty : (*iter_jva).first;
204  size_type jb = iter_jvb == last_jvb ? infty : (*iter_jvb).first;
205  if (ja == jb) {
206  iter_jva++;
207  iter_jvb++;
208  } else if (ja < jb) {
209  iter_jva++;
210  } else {
211  iter_jvb++;
212  }
213  nnz_c++;
214  }
215  }
216  resize (a.row_ownership(), b.col_ownership(), nnz_c);
217  data_iterator iter_jvc = _data.begin().operator->();
218  iterator ic = begin();
219  *ic++ = iter_jvc;
220  for (size_type i = 0, n = a.nrow(); i < n; i++) {
221  for (const_data_iterator iter_jva = ia[i], last_jva = ia[i+1],
222  iter_jvb = ib[i], last_jvb = ib[i+1];
223  iter_jva != last_jva || iter_jvb != last_jvb; ) {
224 
225  size_type ja = iter_jva == last_jva ? infty : (*iter_jva).first;
226  size_type jb = iter_jvb == last_jvb ? infty : (*iter_jvb).first;
227  if (ja == jb) {
228  *iter_jvc++ = std::pair<size_type,T> (ja, binop((*iter_jva).second, (*iter_jvb).second));
229  iter_jva++;
230  iter_jvb++;
231  } else if (ja < jb) {
232  *iter_jvc++ = std::pair<size_type,T> (ja, binop((*iter_jva).second, T(0)));
233  iter_jva++;
234  } else {
235  *iter_jvc++ = std::pair<size_type,T> (jb, binop(T(0), (*iter_jvb).second));
236  iter_jvb++;
237  }
238  }
239  *ic++ = iter_jvc;
240  }
241  set_symmetry (a.is_symmetric() && b.is_symmetric());
242  set_pattern_dimension (std::max(a.pattern_dimension(), b.pattern_dimension()));
243 }
244 /*@!
245  \vfill \pagebreak \mbox{} \vfill \begin{algorithm}[h]
246  \caption{{\tt sort}: sort rows by increasing column order}
247  \begin{algorithmic}
248  \INPUT {the matrix in CSR format}
249  ia(0:nrow), ja(0:nnz-1), a(0:nnz-1)
250  \ENDINPUT
251  \OUTPUT {the sorted CSR matrix}
252  iao(0:nrow), jao(0:nnzl-1), ao(0:nnzl-1)
253  \ENDOUTPUT
254  \BEGIN
255  {\em first: reset iao} \\
256  \FORTO {i := 0} {nrow}
257  iao(i) := 0;
258  \ENDFOR
259 
260  {\em second: compute lengths from pointers} \\
261  \FORTO {i := 0} {nrow-1}
262  \FORTO {p := ia(i)} {ia(i+1)-1}
263  iao (ja(p)+1)++;
264  \ENDFOR
265  \ENDFOR
266 
267  {\em third: compute pointers from lengths} \\
268  \FORTO {i := 0} {nrow-1}
269  iao(i+1) += iao(i)
270  \ENDFOR
271 
272  {\em fourth: copy values} \\
273  \FORTO {i := 0} {nrow-1}
274  \FORTO {p := ia(i)} {ia(i+1)-1}
275  j := ja(p) \\
276  q := iao(j) \\
277  jao(q) := i \\
278  ao(q) := a(p) \\
279  iao (j)++
280  \ENDFOR
281  \ENDFOR
282 
283  {\em fiveth: reshift pointers} \\
284  \FORTOSTEP {i := nrow-1} {0} {-1}
285  iao(i+1) := iao(i);
286  \ENDFOR
287  iao(0) := 0
288  \END
289  \end{algorithmic} \end{algorithm}
290  \vfill \pagebreak \mbox{} \vfill
291 */
292 
293 template<class T>
294 void
296 {
297  b.resize (col_ownership(), row_ownership(), nnz());
298 
299  iterator ib = b.begin();
300  for (size_type j = 0, m = b.nrow(); j < m; j++) {
301  ib[j+1] = ib[0];
302  }
303  const_iterator ia = begin();
304  for (size_type i = 0, n = nrow(); i < n; i++) {
305  for (const_data_iterator p = ia[i], last_p = ia[i+1]; p != last_p; p++) {
306  size_type j = (*p).first;
307  ib [j+1]++;
308  }
309  }
310  for (size_type j = 0, m = b.nrow(); j < m; j++) {
311  ib [j+1] += (ib[j]-ib[0]);
312  }
313  data_iterator q0 = ib[0];
314  for (size_type i = 0, n = nrow(); i < n; i++) {
315  for (const_data_iterator p = ia[i], last_p = ia[i+1]; p != last_p; p++) {
316  size_type j = (*p).first;
317  data_iterator q = ib[j];
318  (*q).first = i;
319  (*q).second = (*p).second;
320  ib[j]++;
321  }
322  }
323  for (long int j = b.nrow()-1; j >= 0; j--) {
324  ib[j+1] = ib[j];
325  }
326  ib[0] = q0;
327 }
328 template<class T>
329 void
331 {
332  if (nrow() != ncol()) {
333  set_symmetry(false);
334  return;
335  }
337  build_transpose (at);
338  d.assign_add (*this, at, std::minus<T>());
339  set_symmetry(d.max_abs() <= tol);
340 }
341 template<class T>
344  const csr_rep<T,sequential>& a,
345  const csr_rep<T,sequential>& b)
346 {
348  typename csr_rep<T,sequential>::const_iterator ia = a.begin(), ib = b.begin();
349  size_type nnz_c = 0;
350  for (size_type i = 0, n = a.nrow(); i < n; i++) {
351  std::set<size_type> y;
352  for (typename csr<T>::const_data_iterator ap = ia[i], last_ap = ia[i+1]; ap != last_ap; ++ap) {
353  size_type j = (*ap).first;
354  typename std::set<size_type>::iterator iter_y = y.begin();
355  for (typename csr<T>::const_data_iterator bp = ib[j], last_bp = ib[j+1]; bp != last_bp; ++bp) {
356  size_type k = (*bp).first;
357  iter_y = y.insert(iter_y, k);
358  }
359  }
360  nnz_c += y.size();
361  }
362  return nnz_c;
363 }
364 template<class T>
365 static
366 void
368  const csr_rep<T,sequential>& a,
369  const csr_rep<T,sequential>& b,
371 {
373  typename csr_rep<T,sequential>::const_iterator ia = a.begin(), ib = b.begin();
374  typename csr_rep<T,sequential>::iterator ic = c.begin();
375  for (size_type i = 0, n = a.nrow(); i < n; i++) {
376  std::map<size_type,T> y;
377  for (typename csr<T>::const_data_iterator ap = ia[i], last_ap = ia[i+1]; ap != last_ap; ++ap) {
378  size_type j = (*ap).first;
379  const T& a_ij = (*ap).second;
380  typename std::map<size_type,T>::iterator prev_y = y.begin();
381  for (typename csr<T>::const_data_iterator bp = ib[j], last_bp = ib[j+1]; bp != last_bp; ++bp) {
382  size_type k = (*bp).first;
383  const T& b_jk = (*bp).second;
384  T value = a_ij*b_jk;
385  typename std::map<size_type,T>::iterator curr_y = y.find(k);
386  if (curr_y == y.end()) {
387  y.insert (prev_y, std::pair<size_type,T>(k, value));
388  } else {
389  (*curr_y).second += value;
390  prev_y = curr_y;
391  }
392  }
393  }
394  ic[i+1] = ic[i] + y.size();
395  typename std::map<size_type,T>::const_iterator iter_y = y.begin();
396  for (typename csr<T>::data_iterator cp = ic[i], last_cp = ic[i+1]; cp != last_cp; ++cp, ++iter_y) {
397  *cp = *iter_y;
398  }
399  }
400 }
401 template<class T>
402 void
404  const csr_rep<T,sequential>& a,
405  const csr_rep<T,sequential>& b)
406 {
407  size_type nnz_c = csr_csr_mult_size (a, b);
408  resize (a.row_ownership(), b.col_ownership(), nnz_c);
409  csr_csr_mult (a, b, *this);
410 }
411 #ifdef _RHEOLEF_DO_NOT_USE_HEAP_ALLOCATOR
412 #define _RHEOLEF_instanciate_class(T) \
413 template void csr_rep<T,sequential>::build_from_asr (const asr<T,sequential,std::allocator<T> >&);
414 #else // _RHEOLEF_DO_NOT_USE_HEAP_ALLOCATOR
415 #define _RHEOLEF_instanciate_class(T) \
416 template void csr_rep<T,sequential>::build_from_asr (const asr<T,sequential,std::allocator<T> >&); \
417 template void csr_rep<T,sequential>::build_from_asr (const asr<T,sequential,heap_allocator<T> >&);
418 #endif // _RHEOLEF_DO_NOT_USE_HEAP_ALLOCATOR
419 
420 #define _RHEOLEF_instanciate(T) \
421 template class csr_rep<T,sequential>; \
422 template void csr_rep<T,sequential>::assign_add ( \
423  const csr_rep<T,sequential>&, \
424  const csr_rep<T,sequential>&, \
425  std::plus<T>); \
426 template void csr_rep<T,sequential>::assign_add ( \
427  const csr_rep<T,sequential>&, \
428  const csr_rep<T,sequential>&, \
429  std::minus<T>); \
430 _RHEOLEF_instanciate_class(T)
431 
433 
434 } // namespace rheolef
distributor - data distribution table
Definition: distributor.h:19
space_mult_list< T, M > & operator*=(space_mult_list< T, M > &Xm, const space_basic< T, M > &Y)
Definition: space_mult.h:72
size_type ncol() const
Definition: csr.h:43
idiststream, odiststream - distributed interface for large data streams
Definition: diststream.h:68
vector_of_iterator< pair_type >::const_iterator const_iterator
Definition: csr.h:22
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)
static flag_type format_field
Definition: iorheo.h:338
STL namespace.
flag_type flags() const
Definition: iorheo.cc:148
const_iterator begin() const
Definition: csr.h:38
irheostream, orheostream - large data streams
Definition: compiler.h:7
const distributor & row_ownership() const
Definition: asr.h:75
const distributor & row_ownership() const
Definition: csr.h:36
std::bitset< last > flag_type
Definition: iorheo.h:333
vector_of_iterator< pair_type >::const_value_type const_data_iterator
Definition: csr.h:25
vector_of_iterator< pair_type >::value_type data_iterator
Definition: csr.h:24
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:343
const_reference operator[](size_type i) const
#define check_macro(ok_condition, message)
Definition: compiler.h:104
double Float
Definition: compiler.h:160
size_t d
size_type size(size_type ip) const
Definition: distributor.h:110
bool is_symmetric() const
Definition: csr.h:49
rheolef::std value
size_type nrow() const
Definition: csr.h:42
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:367
const distributor & col_ownership() const
Definition: asr.h:76
std::vector< T >::size_type size_type
Definition: csr.h:17
void assign_add(const csr_rep< T, sequential > &a, const csr_rep< T, sequential > &b, BinaryOp binop)
Definition: csr_seq.cc:186
vec - vector in distributed environment
Definition: vec.h:48
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
const distributor & col_ownership() const
Definition: csr.h:37
size_type pattern_dimension() const
Definition: csr.h:54
csr - compressed sparse row matrix
Definition: asr.h:8
vector_of_iterator< pair_type >::iterator iterator
Definition: csr.h:21
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:20
#define _RHEOLEF_instanciate(T)
Definition: csr_seq.cc:420
size_type nnz() const
Definition: asr.h:70