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->();
92 for (
size_type i = 0, n = b.nrow(); i < n; i++) {
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;
123 for (
size_type i = 0, n = nrow(); i < n; i++) {
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());
138 if (name ==
"") name =
"a";
139 os << setprecision (std::numeric_limits<T>::digits10)
140 << name <<
" = spalloc(" << nrow() <<
"," << ncol() <<
"," << nnz() <<
");" << endl;
143 for (
size_type i = 0, n = nrow(); i < n; i++) {
145 os << name <<
"(" << i+first_dis_i+base <<
"," << (*iter).first+base <<
") = "
146 << (*iter).second <<
";" << endl;
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()<<
")");
186 details::generic_set_op(),
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));
203 details::generic_set_plus_op(),
220template<
class BinaryOp>
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()<<
")");
234 const size_type infty = std::numeric_limits<size_type>::max();
237 for (
size_type i = 0, n = a.nrow(); i < n; i++) {
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);
262 for (
size_type i = 0, n = a.nrow(); i < n; i++) {
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());
350 for (
size_type j = 0, m = b.nrow(); j < m; j++) {
355 for (
size_type i = 0, n = nrow(); i < n; i++) {
362 for (
size_type j = 0, m = b.nrow(); j < m; j++) {
363 ib [j+1] += (ib[j]-ib[0]);
367 for (
size_type i = 0, n = nrow(); i < n; i++) {
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);
412 for (
size_type i = 0, n = a.nrow(); i < n; i++) {
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)
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) {
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) {
467 c.set_pattern_dimension (std::max(
a.pattern_dimension(),
b.pattern_dimension()));
476 resize (a.row_ownership(), b.col_ownership(), nnz_c);
477 csr_csr_mult (a, b, *
this);
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> >&);
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> >&);
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>&, \
497template void csr_rep<T,sequential>::assign_add ( \
498 const csr_rep<T,sequential>&, \
499 const csr_rep<T,sequential>&, \
501_RHEOLEF_instanciate_class(T)
field::size_type size_type
see the Float page for the full documentation
vector_of_iterator< pair_type >::const_value_type const_data_iterator
vector_of_iterator< pair_type >::iterator iterator
std::pair< size_type, T > pair_type
std::vector< T >::size_type size_type
vector_of_iterator< pair_type >::value_type data_iterator
vector_of_iterator< pair_type >::const_iterator const_iterator
see the csr page for the full documentation
see the distributor page for the full documentation
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
std::bitset< last > flag_type
static flag_type format_field
odiststream: see the diststream page for the full documentation
see the vec page for the full documentation
const_reference operator[](size_type i) const
#define _RHEOLEF_instanciate(T)
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)
void put(std::ostream &out, std::string name, const tiny_matrix< T > &a)
csr_rep< T, sequential >::size_type csr_csr_mult_size(const csr_rep< T, sequential > &a, const csr_rep< T, sequential > &b)
void put_matrix_market(std::ostream &out, const Eigen::SparseMatrix< T, Eigen::RowMajor > &a)
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)