1 # ifndef _RHEOLEF_CSR_H
2 # define _RHEOLEF_CSR_H
70 #include "rheolef/vec.h"
71 #include "rheolef/asr.h"
72 #include "rheolef/vector_of_iterator.h"
73 #include "rheolef/scatter_message.h"
74 #include "rheolef/pair_util.h"
81 template<
class T,
class M>
class csr_rep {};
129 idiststream& get (idiststream&);
137 template <
class BinaryOp>
162 _row_ownership (
a.row_ownership()),
163 _col_ownership (
a.col_ownership()),
165 _is_symmetric(false),
166 _is_definite_positive(false),
167 _pattern_dimension(0)
176 typedef std::allocator<T> A;
185 #ifdef _RHEOLEF_HAVE_MPI
193 typedef typename base::iterator iterator;
195 typedef typename base::data_iterator data_iterator;
196 typedef typename base::const_data_iterator const_data_iterator;
205 const distributor& row_ownership()
const {
return base::_row_ownership; }
206 const distributor& col_ownership()
const {
return base::_col_ownership; }
207 const communicator& comm()
const {
return row_ownership().
comm(); }
208 const_iterator begin()
const {
return base::begin(); }
209 const_iterator end()
const {
return base::end(); }
210 iterator begin() {
return base::begin(); }
211 iterator end() {
return base::end(); }
212 size_type ext_nnz()
const {
return _ext.nnz(); }
213 const_iterator ext_begin()
const {
return _ext.begin(); }
214 const_iterator ext_end()
const {
return _ext.end(); }
215 iterator ext_begin() {
return _ext.begin(); }
216 iterator ext_end() {
return _ext.end(); }
217 size_type nrow()
const {
return base::nrow(); }
218 size_type ncol()
const {
return base::ncol(); }
219 size_type nnz()
const {
return base::nnz(); }
220 size_type dis_nrow()
const {
return row_ownership().dis_size(); }
221 size_type dis_ncol()
const {
return col_ownership().dis_size(); }
222 size_type dis_nnz()
const {
return _dis_nnz; }
223 size_type dis_ext_nnz()
const {
return _dis_ext_nnz; }
225 bool is_symmetric()
const {
return base::is_symmetric(); }
226 void set_symmetry (
bool is_symm)
const { base::set_symmetry(is_symm); }
228 bool is_definite_positive()
const {
return base::is_definite_positive(); }
229 void set_definite_positive (
bool is_defpos)
const { base::set_definite_positive(is_defpos); }
230 size_type pattern_dimension()
const {
return base::pattern_dimension(); }
231 void set_pattern_dimension(
size_type dim)
const { base::set_pattern_dimension(
dim); }
232 size_type row_first_index ()
const {
return row_ownership().first_index(); }
233 size_type row_last_index ()
const {
return row_ownership().last_index(); }
234 size_type col_first_index ()
const {
return col_ownership().first_index(); }
235 size_type col_last_index ()
const {
return col_ownership().last_index(); }
237 idiststream& get (idiststream&);
238 odiststream&
put (odiststream&)
const;
239 void dump (
const std::string&
name)
const;
240 void mult (
const vec<T,distributed>& x, vec<T,distributed>& y)
const;
241 void trans_mult (
const vec<T,distributed>& x, vec<T,distributed>& y)
const;
243 template <
class BinaryOp>
244 void assign_add (
const csr_rep<T,distributed>&
a,
const csr_rep<T,distributed>&
b, BinaryOp binop);
245 void build_transpose (csr_rep<T,distributed>&
b)
const;
246 void assign_mult (
const csr_rep<T,distributed>&
a,
const csr_rep<T,distributed>&
b);
251 csr_rep<T,sequential> _ext;
252 std::vector<size_type> _jext2dis_j;
257 mutable bool _scatter_initialized;
258 mutable scatter_message<std::vector<T> > _from;
259 mutable scatter_message<std::vector<T> > _to;
260 mutable std::vector<T> _buffer;
262 void _scatter_init()
const;
263 void _scatter_init_guard()
const {
264 if (_scatter_initialized)
return;
265 _scatter_initialized =
true;
272 csr_rep<T,distributed>::jext2dis_j (
size_type jext)
const
274 check_macro (jext < _jext2dis_j.size(),
"jext2dis_j: jext="<<jext<<
" is out of range [0:"<<_jext2dis_j.size()<<
"[");
275 return _jext2dis_j [jext];
280 csr_rep<T,distributed>::csr_rep(
const asr<T,distributed,A>&
a)
281 : csr_rep<
T,sequential>(),
286 _scatter_initialized(false),
296 csr_rep<T,distributed>::get (idiststream& ips)
298 typedef std::allocator<T> A;
299 asr<T,distributed,A>
a;
308 template <
class T,
class M>
class csr_concat_value;
309 template <
class T,
class M>
class csr_concat_line;
316 template <
class T,
class M = rheo_default_memory_model>
344 { base::data().resize(loc_nrow1, loc_ncol1, loc_nnz1); }
346 { base::data().resize(row_ownership, col_ownership, nnz1); }
363 void set_symmetry (
bool is_symm)
const { base::data().set_symmetry(is_symm); }
365 { base::data().set_symmetry_by_check(); }
370 T max_abs ()
const {
return base::data().max_abs(); }
402 base::data().mult (x,y);
410 base::data().trans_mult (x,y);
424 base::data().operator*= (
lambda);
458 a.data().build_transpose (
b.data());
463 #ifdef _RHEOLEF_HAVE_MPI
466 class csr<
T,distributed> :
public smart_pointer<csr_rep<T,distributed> > {
471 typedef csr_rep<T,distributed> rep;
472 typedef smart_pointer<rep> base;
476 typedef typename rep::iterator iterator;
477 typedef typename rep::const_iterator const_iterator;
478 typedef typename rep::data_iterator data_iterator;
479 typedef typename rep::const_data_iterator const_data_iterator;
483 csr() : base(new_macro(rep())) {}
485 explicit csr(
const asr<T,memory_type,A>&
a) : base(new_macro(rep(
a))) {}
486 void resize (
const distributor& row_ownership,
const distributor& col_ownership,
size_type nnz1 = 0)
487 { base::data().resize(row_ownership, col_ownership, nnz1); }
491 csr (
const std::initializer_list<details::csr_concat_value<T,distributed> >& init_list);
492 csr (
const std::initializer_list<details::csr_concat_line<T,distributed> >& init_list);
497 const distributor& row_ownership()
const {
return base::data().row_ownership(); }
498 const distributor& col_ownership()
const {
return base::data().col_ownership(); }
499 size_type dis_nrow ()
const {
return row_ownership().dis_size(); }
500 size_type dis_ncol ()
const {
return col_ownership().dis_size(); }
501 size_type dis_nnz ()
const {
return base::data().dis_nnz(); }
502 size_type dis_ext_nnz ()
const {
return base::data().dis_ext_nnz(); }
503 bool is_symmetric()
const {
return base::data().is_symmetric(); }
504 void set_symmetry (
bool is_symm)
const { base::data().set_symmetry(is_symm); }
506 { base::data().set_symmetry_by_check(); }
507 bool is_definite_positive()
const {
return base::data().is_definite_positive(); }
508 void set_definite_positive (
bool is_defpos)
const { base::data().set_definite_positive(is_defpos); }
509 size_type pattern_dimension()
const {
return base::data().pattern_dimension(); }
510 void set_pattern_dimension(
size_type dim)
const { base::data().set_pattern_dimension(
dim); }
511 T max_abs ()
const {
return base::data().max_abs(); }
514 size_type nrow ()
const {
return base::data().nrow(); }
515 size_type ncol ()
const {
return base::data().ncol(); }
516 size_type nnz ()
const {
return base::data().nnz(); }
519 size_type row_first_index ()
const {
return base::data().row_first_index(); }
520 size_type row_last_index ()
const {
return base::data().row_last_index(); }
521 size_type col_first_index ()
const {
return base::data().col_first_index(); }
522 size_type col_last_index ()
const {
return base::data().col_last_index(); }
524 const_iterator begin()
const {
return base::data().begin(); }
525 const_iterator end()
const {
return base::data().end(); }
526 iterator begin_nonconst() {
return base::data().begin(); }
527 iterator end_nonconst() {
return base::data().end(); }
530 size_type ext_nnz()
const {
return base::data().ext_nnz(); }
531 const_iterator ext_begin()
const {
return base::data().ext_begin(); }
532 const_iterator ext_end()
const {
return base::data().ext_end(); }
533 iterator ext_begin_nonconst() {
return base::data().ext_begin(); }
534 iterator ext_end_nonconst() {
return base::data().ext_end(); }
537 int constraint_process_rank()
const;
542 void mult (
const vec<element_type,distributed>& x, vec<element_type,distributed>& y)
const {
543 base::data().mult (x,y);
545 vec<element_type,distributed>
operator* (
const vec<element_type,distributed>& x)
const {
546 vec<element_type,distributed> y (row_ownership(),
element_type());
550 void trans_mult (
const vec<element_type,distributed>& x, vec<element_type,distributed>& y)
const {
551 base::data().trans_mult (x,y);
553 vec<element_type,distributed> trans_mult (
const vec<element_type,distributed>& x)
const {
554 vec<element_type,distributed> y (col_ownership(),
element_type());
559 csr<T,distributed>
operator+ (
const csr<T,distributed>&
b)
const;
560 csr<T,distributed>
operator- (
const csr<T,distributed>&
b)
const;
561 csr<T,distributed>
operator* (
const csr<T,distributed>&
b)
const;
565 base::data().operator*= (
lambda);
570 void dump (
const std::string&
name)
const { base::data().dump(
name); }
578 csr<T,distributed>
b =
a;
594 trans (
const csr<T,distributed>&
a)
596 csr<T,distributed>
b;
597 a.data().build_transpose (
b.data());
603 template<
class T,
class M,
class Function>
612 if (
a.ext_nnz() != 0) {
619 template<
class T,
class M,
class Function>
627 template<
class T,
class M>
629 diag (
const vec<T,M>&);
634 template <
class T,
class M>
639 return x.data().get(s);
641 template <
class T,
class M>
646 return x.data().put(s);
656 c.data().assign_add (this->data(),
b.data(), std::plus<T>());
664 c.data().assign_add (this->data(),
b.data(), std::minus<T>());
667 #ifdef _RHEOLEF_HAVE_MPI
673 c.data().assign_add (this->data(),
b.data(), std::plus<T>());
680 csr<T,distributed>
c;
681 c.data().assign_add (this->data(),
b.data(), std::minus<T>());
694 c.data().assign_mult (this->data(),
b.data());
697 #ifdef _RHEOLEF_HAVE_MPI
703 c.data().assign_mult (this->data(),
b.data());
719 val = std::max (val, abs((*iter).second));
723 #ifdef _RHEOLEF_HAVE_MPI
730 val = mpi::all_reduce (comm(), val, mpi::maximum<T>());
field::size_type size_type
size_type dis_nnz() const
rep::const_data_iterator const_data_iterator
rep::element_type element_type
size_type jext2dis_j(size_type jext) const
void resize(const distributor &row_ownership, const distributor &col_ownership, size_type nnz1=0)
rep::const_iterator const_iterator
void resize(size_type loc_nrow1=0, size_type loc_ncol1=0, size_type loc_nnz1=0)
size_type row_last_index() const
const_iterator begin() const
size_type pattern_dimension() const
bool is_definite_positive() const
size_type col_last_index() const
void mult(const vec< element_type, sequential > &x, vec< element_type, sequential > &y) const
void set_pattern_dimension(size_type dim) const
iterator ext_begin_nonconst()
const distributor & col_ownership() const
csr(const std::initializer_list< details::csr_concat_line< T, sequential > > &init_list)
iterator ext_end_nonconst()
const_iterator ext_end() const
size_type dis_ncol() const
void set_symmetry(bool is_symm) const
void set_symmetry_by_check(const T &tol=std::numeric_limits< T >::epsilon()) const
iterator begin_nonconst()
vec< element_type, sequential > trans_mult(const vec< element_type, sequential > &x) const
size_type ext_nnz() const
csr(const asr< T, sequential, A > &a)
smart_pointer< rep > base
bool is_symmetric() const
csr(const std::initializer_list< details::csr_concat_value< T, sequential > > &init_list)
size_type row_first_index() const
size_type dis_ext_nnz() const
const distributor & row_ownership() const
rep::memory_type memory_type
void set_definite_positive(bool is_defpos) const
const_iterator end() const
void trans_mult(const vec< element_type, sequential > &x, vec< element_type, sequential > &y) const
int constraint_process_rank() const
rep::data_iterator data_iterator
size_type col_first_index() const
size_type dis_nrow() const
csr_rep< T, sequential > rep
const_iterator ext_begin() const
void dump(const std::string &name) const
size_type dis_nnz() const
size_type _pattern_dimension
size_type jext2dis_j(size_type jext) const
vector_of_iterator< pair_type >::const_value_type const_data_iterator
size_type row_last_index() const
const_iterator begin() const
size_type pattern_dimension() const
bool is_definite_positive() const
size_type col_last_index() const
vector_of_iterator< pair_type >::iterator iterator
void set_pattern_dimension(size_type dim) const
const distributor & col_ownership() const
const_iterator ext_end() const
std::vector< std::pair< typename std::vector< T >::size_type, T > > _data
size_type dis_ncol() const
void set_symmetry(bool is_symm) const
distributor _row_ownership
std::pair< size_type, T > pair_type
size_type ext_nnz() const
bool is_symmetric() const
vector_of_iterator< pair_type >::difference_type difference_type
std::vector< T >::size_type size_type
vector_of_iterator< pair_type >::value_type data_iterator
size_type row_first_index() const
size_type dis_ext_nnz() const
const distributor & row_ownership() const
bool _is_definite_positive
void set_definite_positive(bool is_defpos) const
const_iterator end() const
vector_of_iterator< pair_type >::const_iterator const_iterator
distributor _col_ownership
size_type col_first_index() const
size_type dis_nrow() const
const_iterator ext_begin() const
see the csr page for the full documentation
see the distributor page for the full documentation
const communicator_type & comm() const
odiststream: see the diststream page for the full documentation
see the smart_pointer page for the full documentation
see the vec page for the full documentation
V::difference_type difference_type
const_value_type *const const_iterator
const T * const_value_type
const_iterator end() const
geo_element_auto< heap_allocator< geo_element::size_type > > element_type
typename Expr1::memory_type memory_type
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 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)
std::istream & operator>>(std::istream &is, const catchmark &m)
csr< T, sequential > trans(const csr< T, sequential > &a)
trans(a): see the form page for the full documentation
std::enable_if< details::is_rheolef_arithmetic< U >::value,ad3_basic< T >>::type operator+(const U &a, const ad3_basic< T > &b)
csr< T, sequential > operator-(const csr< T, sequential > &a)
csr< T, M > diag(const vec< T, M > &d)
csr< T, sequential > operator*(const T &lambda, const csr< T, sequential > &a)
void put_matrix_market(std::ostream &out, const Eigen::SparseMatrix< T, Eigen::RowMajor > &a)
OutputPairIterator pair_transform_second(InputPairIterator first, InputPairIterator last, OutputPairIterator result, UnaryOperation unary_op)
std::ostream & operator<<(std::ostream &os, const catchmark &m)
csr< T, M > csr_apply(Function f, const csr< T, M > &a)