Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
csr_seq.cc
Go to the documentation of this file.
1
21
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"
28using namespace std;
29namespace rheolef {
30// ----------------------------------------------------------------------------
31// class member functions
32// ----------------------------------------------------------------------------
33
34template<class T>
35csr_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}
45template<class T>
46void
47csr_rep<T,sequential>::resize (const distributor& row_ownership, const distributor& col_ownership, size_type nnz1)
48{
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}
56template<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}
67template<class T>
68void
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}
78template<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}
96template<class T>
97template<class A>
98void
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}
113template<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}
132template<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}
151template<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}
161template<class T>
162void
163csr_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
174template<class T>
175void
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}
189template<class T>
190void
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}
206template<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// ----------------------------------------------------------------------------
219template<class T>
220template<class BinaryOp>
221void
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
342template<class T>
343void
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// ----------------------------------------------------------------------------
385template<class T>
386void
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
403template<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
427template<class T>
428static
429void
430csr_csr_mult (
431 const csr_rep<T,sequential>& a,
432 const csr_rep<T,sequential>& b,
433 csr_rep<T,sequential>& c)
434{
435 typedef typename csr_rep<T,sequential>::size_type size_type;
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}
469template<class T>
470void
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) \
484template 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) \
487template void csr_rep<T,sequential>::build_from_asr (const asr<T,sequential,std::allocator<T> >&); \
488template 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) \
492template class csr_rep<T,sequential>; \
493template void csr_rep<T,sequential>::assign_add ( \
494 const csr_rep<T,sequential>&, \
495 const csr_rep<T,sequential>&, \
496 std::plus<T>); \
497template 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)
size_type size(size_type iproc) const
static const size_type decide
Definition distributor.h:83
std::bitset< last > flag_type
Definition iorheo.h:440
static flag_type format_field
Definition iorheo.h:445
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
#define _RHEOLEF_instanciate(T)
Definition csr_seq.cc:491
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)")
This file is part of Rheolef.
void csr_amux(InputIterator ia, InputIterator last_ia, InputRandomAcessIterator x, SetOperator set_op, OutputIterator y)
Definition csr_amux.h:77
void put(std::ostream &out, std::string name, const tiny_matrix< T > &a)
Definition tiny_lu.h:155
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
STL namespace.
Definition sphere.icc:25