NAME
csr - compressed sparse row matrix (rheolef-7.0)
SYNOPSIS
Distributed compressed sparse matrix container stored row by row.
DESCRIPTION
Sparse matrix are compressed by rows. In distributed environment, the distribution follows the row distributor (see distributor(2)).
ALGEBRA
Adding or subtracting two matrices writes a+b and a-b, respectively, and multiplying a matrix by a scalar writes lambda*x. Thus, any linear combination of sparse matrices is available.
Matrix-vector product writes a*x where x is a vector (see vec(2)).
LIMITATIONS
Some basic linear algebra is still under development: a.trans_mult(x) matrix transpose vector product, trans(a) matrix transpose, a*b matrix product.
IMPLEMENTATION
template<class
T>
class csr<T,sequential> : public
smart_pointer<csr_rep<T,sequential> > {
public:
// typedefs:
typedef
csr_rep<T,sequential> rep;
typedef smart_pointer<rep> base;
typedef typename rep::memory_type memory_type;
typedef typename rep::size_type size_type;
typedef typename rep::element_type element_type;
typedef typename rep::iterator iterator;
typedef typename rep::const_iterator const_iterator;
typedef typename rep::data_iterator data_iterator;
typedef typename rep::const_data_iterator
const_data_iterator;
// allocators/deallocators:
csr() :
base(new_macro(rep())) {}
template<class A>
explicit csr(const asr<T,sequential,A>& a) :
base(new_macro(rep(a))) {}
void resize (size_type loc_nrow1 = 0, size_type loc_ncol1 =
0, size_type loc_nnz1 = 0)
{ base::data().resize(loc_nrow1, loc_ncol1, loc_nnz1); }
void resize (const distributor& row_ownership, const
distributor& col_ownership, size_type nnz1 = 0)
{ base::data().resize(row_ownership, col_ownership, nnz1);
}
// allocators from initializer list (c++ 2011):
#ifdef
_RHEOLEF_HAVE_STD_INITIALIZER_LIST
csr (const
std::initializer_list<csr_concat_value<T,sequential>
>& init_list);
csr (const
std::initializer_list<csr_concat_line<T,sequential>
>& init_list);
#endif // _RHEOLEF_HAVE_STD_INITIALIZER_LIST
// accessors:
// global sizes
const distributor& row_ownership() const { return
base::data().row_ownership(); }
const distributor& col_ownership() const { return
base::data().col_ownership(); }
size_type dis_nrow () const { return
row_ownership().dis_size(); }
size_type dis_ncol () const { return
col_ownership().dis_size(); }
size_type dis_nnz () const { return base::data().nnz(); }
size_type dis_ext_nnz () const { return 0; }
bool is_symmetric() const { return
base::data().is_symmetric(); }
void set_symmetry (bool is_symm) const {
base::data().set_symmetry(is_symm); }
void set_symmetry_by_check (const T& tol =
std::numeric_limits<T>::epsilon()) const
{ base::data().set_symmetry_by_check(); }
bool is_definite_positive() const { return
base::data().is_definite_positive(); }
void set_definite_positive (bool is_defpos) const {
base::data().set_definite_positive(is_defpos); }
size_type pattern_dimension() const { return
base::data().pattern_dimension(); }
void set_pattern_dimension(size_type dim) const {
base::data().set_pattern_dimension(dim); }
T max_abs () const { return base::data().max_abs(); }
// local sizes
size_type nrow () const { return base::data().nrow(); }
size_type ncol () const { return base::data().ncol(); }
size_type nnz () const { return base::data().nnz(); }
// range on
local memory
size_type row_first_index () const { return
base::data().row_first_index(); }
size_type row_last_index () const { return
base::data().row_last_index(); }
size_type col_first_index () const { return
base::data().col_first_index(); }
size_type col_last_index () const { return
base::data().col_last_index(); }
const_iterator
begin() const { return base::data().begin(); }
const_iterator end() const { return base::data().end(); }
iterator begin_nonconst() { return base::data().begin(); }
iterator end_nonconst() { return base::data().end(); }
// accessors,
only for distributed (for interface compatibility)
size_type ext_nnz() const { return 0; }
const_iterator ext_begin() const { return const_iterator();
}
const_iterator ext_end() const { return const_iterator(); }
iterator ext_begin_nonconst() { return iterator(); }
iterator ext_end_nonconst() { return iterator(); }
size_type jext2dis_j (size_type jext) const { return 0;
}
// algebra:
// y := a*x
void mult (const vec<element_type,sequential>& x,
vec<element_type,sequential>& y) const {
base::data().mult (x,y);
}
vec<element_type,sequential> operator* (const
vec<element_type,sequential>& x) const {
vec<element_type,sequential> y (row_ownership(),
element_type());
mult (x, y);
return y;
}
void trans_mult (const
vec<element_type,sequential>& x,
vec<element_type,sequential>& y) const {
base::data().trans_mult (x,y);
}
vec<element_type,sequential> trans_mult (const
vec<element_type,sequential>& x) const {
vec<element_type,sequential> y (col_ownership(),
element_type());
trans_mult (x, y);
return y;
}
// a+b, a-b, a*b
csr<T,sequential> operator+ (const
csr<T,sequential>& b) const;
csr<T,sequential> operator- (const
csr<T,sequential>& b) const;
csr<T,sequential> operator* (const
csr<T,sequential>& b) const;
// lambda*a
csr<T,sequential>& operator*= (const T&
lambda) {
base::data().operator*= (lambda);
return *this;
}
// output:
void dump
(const std::string& name) const {
base::data().dump(name); }
};
// lambda*a
template<class T>
inline
csr<T,sequential>
operator* (const T& lambda, const
csr<T,sequential>& a)
{
csr<T,sequential> b = a;
b.operator*= (lambda);
return b;
}
// -a
template<class T>
inline
csr<T,sequential>
operator- (const csr<T,sequential>& a)
{
return T(-1)*a;
}
// trans(a)
template<class T>
inline
csr<T,sequential>
trans (const csr<T,sequential>& a)
{
csr<T,sequential> b;
a.data().build_transpose (b.data());
return b;
}
IMPLEMENTATION
template<class
T>
class csr<T,distributed> : public
smart_pointer<csr_rep<T,distributed> > {
public:
// typedefs:
typedef
csr_rep<T,distributed> rep;
typedef smart_pointer<rep> base;
typedef typename rep::memory_type memory_type;
typedef typename rep::size_type size_type;
typedef typename rep::element_type element_type;
typedef typename rep::iterator iterator;
typedef typename rep::const_iterator const_iterator;
typedef typename rep::data_iterator data_iterator;
typedef typename rep::const_data_iterator
const_data_iterator;
// allocators/deallocators:
csr() :
base(new_macro(rep())) {}
template<class A>
explicit csr(const asr<T,memory_type,A>& a) :
base(new_macro(rep(a))) {}
void resize (const distributor& row_ownership, const
distributor& col_ownership, size_type nnz1 = 0)
{ base::data().resize(row_ownership, col_ownership, nnz1);
}
// allocators from initializer list (c++ 2011):
#ifdef
_RHEOLEF_HAVE_STD_INITIALIZER_LIST
csr (const
std::initializer_list<csr_concat_value<T,distributed>
>& init_list);
csr (const
std::initializer_list<csr_concat_line<T,distributed>
>& init_list);
#endif // _RHEOLEF_HAVE_STD_INITIALIZER_LIST
// accessors:
// global sizes
const distributor& row_ownership() const { return
base::data().row_ownership(); }
const distributor& col_ownership() const { return
base::data().col_ownership(); }
size_type dis_nrow () const { return
row_ownership().dis_size(); }
size_type dis_ncol () const { return
col_ownership().dis_size(); }
size_type dis_nnz () const { return base::data().dis_nnz();
}
size_type dis_ext_nnz () const { return
base::data().dis_ext_nnz(); }
bool is_symmetric() const { return
base::data().is_symmetric(); }
void set_symmetry (bool is_symm) const {
base::data().set_symmetry(is_symm); }
void set_symmetry_by_check (const T& tol =
std::numeric_limits<T>::epsilon()) const
{ base::data().set_symmetry_by_check(); }
bool is_definite_positive() const { return
base::data().is_definite_positive(); }
void set_definite_positive (bool is_defpos) const {
base::data().set_definite_positive(is_defpos); }
size_type pattern_dimension() const { return
base::data().pattern_dimension(); }
void set_pattern_dimension(size_type dim) const {
base::data().set_pattern_dimension(dim); }
T max_abs () const { return base::data().max_abs(); }
// local sizes
size_type nrow () const { return base::data().nrow(); }
size_type ncol () const { return base::data().ncol(); }
size_type nnz () const { return base::data().nnz(); }
// range on
local memory
size_type row_first_index () const { return
base::data().row_first_index(); }
size_type row_last_index () const { return
base::data().row_last_index(); }
size_type col_first_index () const { return
base::data().col_first_index(); }
size_type col_last_index () const { return
base::data().col_last_index(); }
const_iterator
begin() const { return base::data().begin(); }
const_iterator end() const { return base::data().end(); }
iterator begin_nonconst() { return base::data().begin(); }
iterator end_nonconst() { return base::data().end(); }
// accessors,
only for distributed
size_type ext_nnz() const { return base::data().ext_nnz(); }
const_iterator ext_begin() const { return
base::data().ext_begin(); }
const_iterator ext_end() const { return
base::data().ext_end(); }
iterator ext_begin_nonconst() { return
base::data().ext_begin(); }
iterator ext_end_nonconst() { return base::data().ext_end();
}
size_type jext2dis_j (size_type jext) const { return
base::data().jext2dis_j(jext); }
// algebra:
// y := a*x
void mult (const vec<element_type,distributed>& x,
vec<element_type,distributed>& y) const {
base::data().mult (x,y);
}
vec<element_type,distributed> operator* (const
vec<element_type,distributed>& x) const {
vec<element_type,distributed> y (row_ownership(),
element_type());
mult (x, y);
return y;
}
void trans_mult (const
vec<element_type,distributed>& x,
vec<element_type,distributed>& y) const {
base::data().trans_mult (x,y);
}
vec<element_type,distributed> trans_mult (const
vec<element_type,distributed>& x) const {
vec<element_type,distributed> y (col_ownership(),
element_type());
trans_mult (x, y);
return y;
}
// a+b, a-b, a*b
csr<T,distributed> operator+ (const
csr<T,distributed>& b) const;
csr<T,distributed> operator- (const
csr<T,distributed>& b) const;
csr<T,distributed> operator* (const
csr<T,distributed>& b) const;
// lambda*a
csr<T,distributed>& operator*= (const T&
lambda) {
base::data().operator*= (lambda);
return *this;
}
// output:
void dump
(const std::string& name) const {
base::data().dump(name); }
};
// lambda*a
template<class T>
inline
csr<T,distributed>
operator* (const T& lambda, const
csr<T,distributed>& a)
{
csr<T,distributed> b = a;
b.operator*= (lambda);
return b;
}
// -a
template<class T>
inline
csr<T,distributed>
operator- (const csr<T,distributed>& a)
{
return T(-1)*a;
}
// trans(a)
template<class T>
inline
csr<T,distributed>
trans (const csr<T,distributed>& a)
{
csr<T,distributed> b;
a.data().build_transpose (b.data());
return b;
}
#endif // _RHEOLEF_HAVE_MPI
// b = f(a); f
as a class-function or usual fct
template<class T, class M, class Function>
csr<T,M>
apply (Function f, const csr<T,M>& a)
{
csr<T,M> b = a;
typename csr<T,M>::size_type n = a.nrow();
typename csr<T,M>::const_iterator dia_ia = a.begin();
typename csr<T,M>::iterator dia_ib =
b.begin_nonconst();
pair_transform_second (dia_ia[0], dia_ia[n], dia_ib[0], f);
if (a.ext_nnz() != 0) {
typename csr<T,M>::const_iterator ext_ia =
a.ext_begin();
typename csr<T,M>::iterator ext_ib =
b.ext_begin_nonconst();
pair_transform_second (ext_ia[0], ext_ia[n], ext_ib[0], f);
}
return b;
}
template<class T, class M, class Function>
csr<T,M>
apply (T (*f)(const T&), const csr<T,M>& a)
{
return apply (std::ptr_fun(f), a);
}
SEE ALSO
COPYRIGHT
Copyright (C) 2000-2018 Pierre Saramito <Pierre.Saramito [AT] imag.fr> GPLv3+: GNU GPL version 3 or later <http://gnu.org/licenses/gpl.html>;. This is free software: you are free to change and redistribute it. There is NO WARRANTY, to the extent permitted by law.