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"
37 _row_ownership (row_ownership),
38 _col_ownership (col_ownership),
40 _is_symmetric (false),
41 _is_definite_positive (false),
42 _pattern_dimension (0)
50 _row_ownership = row_ownership;
51 _col_ownership = col_ownership;
59 _row_ownership (
distributor::decide, communicator(), loc_nrow1),
60 _col_ownership (
distributor::decide, communicator(), loc_ncol1),
62 _is_symmetric (false),
63 _is_definite_positive (false),
64 _pattern_dimension (0)
74 _data.resize (loc_nnz1);
81 _row_ownership (
b.row_ownership()),
82 _col_ownership (
b.col_ownership()),
84 _is_symmetric (
b._is_symmetric),
85 _is_definite_positive (
b._is_definite_positive),
86 _pattern_dimension (
b._pattern_dimension)
91 ia[0] = _data.begin().operator->();
93 ia [i+1] = ia[0] + (ib[i+1] - ib[0]);
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;
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;
125 os << i+first_dis_i+base <<
" "
126 << (*iter).first+base <<
" "
127 << (*iter).second << endl;
136 std::ostream& os = ods.
os();
137 std::string
name = iorheo::getbasename(ods.
os());
139 os << setprecision (std::numeric_limits<T>::digits10)
140 <<
name <<
" = spalloc(" << nrow() <<
"," << ncol() <<
"," << nnz() <<
");" << endl;
145 os <<
name <<
"(" << i+first_dis_i+base <<
"," << (*iter).first+base <<
") = "
146 << (*iter).second <<
";" << endl;
155 iorheo::flag_type format = iorheo::flags(ods.
os()) & iorheo::format_field;
157 {
return put_sparse_matlab (ods,first_dis_i); }
165 std::ofstream os (
name.c_str());
166 std::cerr <<
"! file \"" <<
name <<
"\" created." << std::endl;
181 check_macro (x.size() == ncol(),
"csr*vec: incompatible csr("<<nrow()<<
","<<ncol()<<
") and vec("<<x.size()<<
")");
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));
220 template<
class BinaryOp>
228 "incompatible csr add(a,b): a("<<
a.nrow()<<
":"<<
a.ncol()<<
") and "
229 "b("<<
b.nrow()<<
":"<<
b.ncol()<<
")");
234 const size_type infty = std::numeric_limits<size_type>::max();
239 iter_jvb = ib[i], last_jvb = ib[i+1];
240 iter_jva != last_jva || iter_jvb != last_jvb; ) {
242 size_type ja = iter_jva == last_jva ? infty : (*iter_jva).first;
243 size_type jb = iter_jvb == last_jvb ? infty : (*iter_jvb).first;
247 }
else if (ja < jb) {
255 resize (
a.row_ownership(),
b.col_ownership(), nnz_c);
264 iter_jvb = ib[i], last_jvb = ib[i+1];
265 iter_jva != last_jva || iter_jvb != last_jvb; ) {
267 size_type ja = iter_jva == last_jva ? infty : (*iter_jva).first;
268 size_type jb = iter_jvb == last_jvb ? infty : (*iter_jvb).first;
270 *iter_jvc++ = std::pair<size_type,T> (ja, binop((*iter_jva).second, (*iter_jvb).second));
273 }
else if (ja < jb) {
274 *iter_jvc++ = std::pair<size_type,T> (ja, binop((*iter_jva).second,
T(0)));
277 *iter_jvc++ = std::pair<size_type,T> (jb, binop(
T(0), (*iter_jvb).second));
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());
346 b.resize (col_ownership(), row_ownership(), nnz());
363 ib [j+1] += (ib[j]-ib[0]);
372 (*q).second = (*p).second;
377 for (
long int j =
b.nrow()-1; j >= 0; j--) {
389 if (nrow() != ncol()) {
395 build_transpose (at);
396 d.assign_add (*
this, at, std::minus<T>());
397 set_symmetry(
d.max_abs() <= tol);
413 std::set<size_type> y;
416 typename std::set<size_type>::iterator iter_y = y.begin();
419 iter_y = y.insert(iter_y, k);
431 const csr_rep<T,sequential>&
a,
432 const csr_rep<T,sequential>&
b,
433 csr_rep<T,sequential>&
c)
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) {
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) {
446 const T& b_jk = (*bp).second;
448 typename std::map<size_type,T>::iterator curr_y = y.find(k);
449 if (curr_y == y.end()) {
451 y.insert (prev_y, std::pair<size_type,T>(k,
value));
454 (*curr_y).second +=
value;
459 ic[i+1] = ic[i] + y.size();
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) {
474 resize (
a.row_ownership(),
b.col_ownership(), nnz_c);
475 csr_csr_mult (
a,
b, *
this);
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
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>&, \
495 template void csr_rep<T,sequential>::assign_add ( \
496 const csr_rep<T,sequential>&, \
497 const csr_rep<T,sequential>&, \
499 _RHEOLEF_instanciate_class(T)