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;
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> >&);
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> >&);
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)
void put(idiststream &in, odiststream &out, bool do_proj, bool do_lumped_mass, bool def_fill_opt, size_type extract_id, const Float &scale_value, const std::pair< Float, Float > &u_range, render_type render)
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
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)")
verbose clean transpose logscale grid shrink ball stereo iso volume skipvtk deformation fastfieldload lattice reader_on_stdin color format format format format format format format format format format format format format format format matlab
verbose clean transpose logscale grid shrink ball stereo iso volume skipvtk deformation fastfieldload lattice reader_on_stdin color format format format format format format format format format format format format format format format format format format dump
This file is part of Rheolef.
std::enable_if< details::is_rheolef_arithmetic< U >::value,ad3_basic< T > & >::type operator*=(ad3_basic< T > &a, const U &b)
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)
void csr_amux(InputIterator ia, InputIterator last_ia, InputRandomAcessIterator x, SetOperator set_op, OutputIterator y)