Loading...
Searching...
No Matches
csr_seq.cc
Go to the documentation of this file.
#include "rheolef/csr.h"
#include "rheolef/asr.h"
#include "rheolef/asr_to_csr.h"
#include "rheolef/msg_util.h"
#include "rheolef/csr_amux.h"
#include "rheolef/csr_cumul_trans_mult.h"
using namespace std;
namespace rheolef {
// ----------------------------------------------------------------------------
// class member functions
// ----------------------------------------------------------------------------
template<class T>
csr_rep<T,sequential>::csr_rep(const distributor& row_ownership, const distributor& col_ownership, size_type nnz1)
: vector_of_iterator<pair_type>(row_ownership.size()+1),
_row_ownership (row_ownership),
_col_ownership (col_ownership),
_data (nnz1),
_is_symmetric (false),
_is_definite_positive (false),
_pattern_dimension (0)
{
}
template<class T>
void
csr_rep<T,sequential>::resize (const distributor& row_ownership, const distributor& col_ownership, size_type nnz1)
{
vector_of_iterator<pair_type>::resize (row_ownership.size()+1);
_row_ownership = row_ownership;
_col_ownership = col_ownership;
_data.resize (nnz1);
// first pointer points to the beginning of the data:
vector_of_iterator<pair_type>::operator[](0) = _data.begin().operator->();
}
template<class T>
csr_rep<T,sequential>::csr_rep(size_type loc_nrow1, size_type loc_ncol1, size_type loc_nnz1)
: vector_of_iterator<pair_type> (loc_nrow1+1),
_row_ownership (distributor::decide, communicator(), loc_nrow1),
_col_ownership (distributor::decide, communicator(), loc_ncol1),
_data (loc_nnz1),
_is_symmetric (false),
_is_definite_positive (false),
_pattern_dimension (0)
{
}
template<class T>
void
csr_rep<T,sequential>::resize (size_type loc_nrow1, size_type loc_ncol1, size_type loc_nnz1)
{
vector_of_iterator<pair_type>::resize (loc_nrow1+1);
_row_ownership = distributor (distributor::decide, communicator(), loc_nrow1);
_col_ownership = distributor (distributor::decide, communicator(), loc_ncol1);
_data.resize (loc_nnz1);
// first pointer points to the beginning of the data:
vector_of_iterator<pair_type>::operator[](0) = _data.begin().operator->();
}
template<class T>
csr_rep<T,sequential>::csr_rep(const csr_rep<T,sequential>& b)
: vector_of_iterator<pair_type>(b.nrow()+1),
_row_ownership (b.row_ownership()),
_col_ownership (b.col_ownership()),
_data(b._data),
_is_symmetric (b._is_symmetric),
_is_definite_positive (b._is_definite_positive),
_pattern_dimension (b._pattern_dimension)
{
// physical copy of csr
const_iterator ib = b.begin();
iterator ia = begin();
ia[0] = _data.begin().operator->();
for (size_type i = 0, n = b.nrow(); i < n; i++) {
ia [i+1] = ia[0] + (ib[i+1] - ib[0]);
}
}
template<class T>
template<class A>
void
csr_rep<T,sequential>::build_from_asr (const asr<T,sequential,A>& a)
{
resize (a.row_ownership(), a.col_ownership(), a.nnz());
typedef std::pair<size_type,T> pair_type;
typedef std::pair<size_type,T> const_pair_type;
//typedef typename asr<T>::row_type::value_type const_pair_type;
asr_to_csr (
a.begin(),
a.end(),
always_true<const_pair_type>(),
pair_identity<const_pair_type,pair_type>(),
vector_of_iterator<pair_type>::begin(),
_data.begin());
}
template<class T>
odiststream&
csr_rep<T,sequential>::put_matrix_market (odiststream& ods, size_type first_dis_i) const
{
std::ostream& os = ods.os();
os << setprecision (std::numeric_limits<T>::digits10)
<< "%%MatrixMarket matrix coordinate real general" << std::endl
<< nrow() << " " << ncol() << " " << nnz() << endl;
const_iterator ia = begin();
const size_type base = 1;
for (size_type i = 0, n = nrow(); i < n; i++) {
for (const_data_iterator iter = ia[i], last = ia[i+1]; iter != last; ++iter) {
os << i+first_dis_i+base << " "
<< (*iter).first+base << " "
<< (*iter).second << endl;
}
}
return ods;
}
template<class T>
odiststream&
csr_rep<T,sequential>::put_sparse_matlab (odiststream& ods, size_type first_dis_i) const
{
std::ostream& os = ods.os();
std::string name = iorheo::getbasename(ods.os());
if (name == "") name = "a";
os << setprecision (std::numeric_limits<T>::digits10)
<< name << " = spalloc(" << nrow() << "," << ncol() << "," << nnz() << ");" << endl;
const_iterator ia = begin();
const size_type base = 1;
for (size_type i = 0, n = nrow(); i < n; i++) {
for (const_data_iterator iter = ia[i], last = ia[i+1]; iter != last; ++iter) {
os << name << "(" << i+first_dis_i+base << "," << (*iter).first+base << ") = "
<< (*iter).second << ";" << endl;
}
}
return ods;
}
template<class T>
odiststream&
csr_rep<T,sequential>::put (odiststream& ods, size_type first_dis_i) const
{
iorheo::flag_type format = iorheo::flags(ods.os()) & iorheo::format_field;
if (format [iorheo::matlab] || format [iorheo::sparse_matlab])
{ return put_sparse_matlab (ods,first_dis_i); }
// default is matrix market format
return put_matrix_market (ods,first_dis_i);
}
template<class T>
void
csr_rep<T,sequential>::dump (const string& name, size_type first_dis_i) const
{
std::ofstream os (name.c_str());
std::cerr << "! file \"" << name << "\" created." << std::endl;
odiststream ods(os);
put (ods);
}
// ----------------------------------------------------------------------------
// basic linear algebra
// ----------------------------------------------------------------------------
template<class T>
void
csr_rep<T,sequential>::mult(
const vec<T,sequential>& x,
vec<T,sequential>& y)
const
{
check_macro (x.size() == ncol(), "csr*vec: incompatible csr("<<nrow()<<","<<ncol()<<") and vec("<<x.size()<<")");
csr_amux (
vector_of_iterator<pair_type>::begin(),
vector_of_iterator<pair_type>::end(),
x.begin(),
details::generic_set_op(),
y.begin());
}
template<class T>
void
csr_rep<T,sequential>::trans_mult(
const vec<T,sequential>& x,
vec<T,sequential>& y)
const
{
// reset y and then cumul
check_macro (x.size() == nrow(), "csr.trans_mult(vec): incompatible csr("<<nrow()<<","<<ncol()<<") and vec("<<x.size()<<")");
std::fill (y.begin(), y.end(), T(0));
csr_cumul_trans_mult (
vector_of_iterator<pair_type>::begin(),
vector_of_iterator<pair_type>::end(),
x.begin(),
details::generic_set_plus_op(),
y.begin());
}
template<class T>
csr_rep<T,sequential>&
csr_rep<T,sequential>::operator*= (const T& lambda)
{
iterator ia = begin();
for (data_iterator p = ia[0], last_p = ia[nrow()]; p != last_p; p++) {
(*p).second *= lambda;
}
return *this;
}
// ----------------------------------------------------------------------------
// expression c=a+b and c=a-b with a temporary c=*this
// ----------------------------------------------------------------------------
template<class T>
template<class BinaryOp>
void
csr_rep<T,sequential>::assign_add (
const csr_rep<T,sequential>& a,
const csr_rep<T,sequential>& b,
BinaryOp binop)
{
check_macro (a.nrow() == b.nrow() && a.ncol() == b.ncol(),
"incompatible csr add(a,b): a("<<a.nrow()<<":"<<a.ncol()<<") and "
"b("<<b.nrow()<<":"<<b.ncol()<<")");
//
// first pass: compute nnz_c and resize
//
size_type nnz_c = 0;
const size_type infty = std::numeric_limits<size_type>::max();
const_iterator ia = a.begin();
const_iterator ib = b.begin();
for (size_type i = 0, n = a.nrow(); i < n; i++) {
for (const_data_iterator iter_jva = ia[i], last_jva = ia[i+1],
iter_jvb = ib[i], last_jvb = ib[i+1];
iter_jva != last_jva || iter_jvb != last_jvb; ) {
size_type ja = iter_jva == last_jva ? infty : (*iter_jva).first;
size_type jb = iter_jvb == last_jvb ? infty : (*iter_jvb).first;
if (ja == jb) {
iter_jva++;
iter_jvb++;
} else if (ja < jb) {
iter_jva++;
} else {
iter_jvb++;
}
nnz_c++;
}
}
resize (a.row_ownership(), b.col_ownership(), nnz_c);
data_iterator iter_jvc = _data.begin().operator->();
iterator ic = begin();
*ic++ = iter_jvc;
//
// second pass: add and store in c
//
for (size_type i = 0, n = a.nrow(); i < n; i++) {
for (const_data_iterator iter_jva = ia[i], last_jva = ia[i+1],
iter_jvb = ib[i], last_jvb = ib[i+1];
iter_jva != last_jva || iter_jvb != last_jvb; ) {
size_type ja = iter_jva == last_jva ? infty : (*iter_jva).first;
size_type jb = iter_jvb == last_jvb ? infty : (*iter_jvb).first;
if (ja == jb) {
*iter_jvc++ = std::pair<size_type,T> (ja, binop((*iter_jva).second, (*iter_jvb).second));
iter_jva++;
iter_jvb++;
} else if (ja < jb) {
*iter_jvc++ = std::pair<size_type,T> (ja, binop((*iter_jva).second, T(0)));
iter_jva++;
} else {
*iter_jvc++ = std::pair<size_type,T> (jb, binop(T(0), (*iter_jvb).second));
iter_jvb++;
}
}
*ic++ = iter_jvc;
}
set_symmetry (a.is_symmetric() && b.is_symmetric());
set_pattern_dimension (std::max(a.pattern_dimension(), b.pattern_dimension()));
set_definite_positive (a.is_definite_positive() || b.is_definite_positive());
// perhaps too optimist for definite_positive: should be ajusted before calling a solver
}
// ----------------------------------------------------------------------------
// trans(a)
// ----------------------------------------------------------------------------
/*@!
\vfill \pagebreak \mbox{} \vfill \begin{algorithm}[h]
\caption{{\tt sort}: sort rows by increasing column order}
\begin{algorithmic}
\INPUT {the matrix in CSR format}
ia(0:nrow), ja(0:nnz-1), a(0:nnz-1)
\ENDINPUT
\OUTPUT {the sorted CSR matrix}
iao(0:nrow), jao(0:nnzl-1), ao(0:nnzl-1)
\ENDOUTPUT
\BEGIN
{\em first: reset iao} \\
\FORTO {i := 0} {nrow}
iao(i) := 0;
\ENDFOR
{\em second: compute lengths from pointers} \\
\FORTO {i := 0} {nrow-1}
\FORTO {p := ia(i)} {ia(i+1)-1}
iao (ja(p)+1)++;
\ENDFOR
\ENDFOR
{\em third: compute pointers from lengths} \\
\FORTO {i := 0} {nrow-1}
iao(i+1) += iao(i)
\ENDFOR
{\em fourth: copy values} \\
\FORTO {i := 0} {nrow-1}
\FORTO {p := ia(i)} {ia(i+1)-1}
j := ja(p) \\
q := iao(j) \\
jao(q) := i \\
ao(q) := a(p) \\
iao (j)++
\ENDFOR
\ENDFOR
{\em fiveth: reshift pointers} \\
\FORTOSTEP {i := nrow-1} {0} {-1}
iao(i+1) := iao(i);
\ENDFOR
iao(0) := 0
\END
\end{algorithmic} \end{algorithm}
\vfill \pagebreak \mbox{} \vfill
*/
// NOTE: transposed matrix has always rows sorted by increasing column indexes
// even if original matrix has not
template<class T>
void
csr_rep<T,sequential>::build_transpose (csr_rep<T,sequential>& b) const
{
b.resize (col_ownership(), row_ownership(), nnz());
// first pass: set ib(*) to ib(0)
iterator ib = b.begin();
for (size_type j = 0, m = b.nrow(); j < m; j++) {
ib[j+1] = ib[0];
}
// second pass: compute lengths of row(i) of a^T in ib(i+1)-ib(0)
const_iterator ia = begin();
for (size_type i = 0, n = nrow(); i < n; i++) {
for (const_data_iterator p = ia[i], last_p = ia[i+1]; p != last_p; p++) {
size_type j = (*p).first;
ib [j+1]++;
}
}
// third pass: compute pointers from lengths
for (size_type j = 0, m = b.nrow(); j < m; j++) {
ib [j+1] += (ib[j]-ib[0]);
}
// fourth pass: store values
data_iterator q0 = ib[0];
for (size_type i = 0, n = nrow(); i < n; i++) {
for (const_data_iterator p = ia[i], last_p = ia[i+1]; p != last_p; p++) {
size_type j = (*p).first;
data_iterator q = ib[j];
(*q).first = i;
(*q).second = (*p).second;
ib[j]++;
}
}
// fiveth: shift pointers
for (long int j = b.nrow()-1; j >= 0; j--) {
ib[j+1] = ib[j];
}
ib[0] = q0;
}
// ----------------------------------------------------------------------------
// set symmetry by check
// ----------------------------------------------------------------------------
template<class T>
void
csr_rep<T,sequential>::set_symmetry_by_check (const T& tol) const
{
if (nrow() != ncol()) {
set_symmetry(false);
return;
}
// check if a-trans(a) == 0 up to machine prec
csr_rep<T,sequential> at, d;
build_transpose (at);
d.assign_add (*this, at, std::minus<T>());
set_symmetry(d.max_abs() <= tol);
}
// ----------------------------------------------------------------------------
// expression c=a*b with a temporary c=*this
// ----------------------------------------------------------------------------
// compute the sparse size of c=a*b
template<class T>
typename csr_rep<T,sequential>::size_type
csr_csr_mult_size (
const csr_rep<T,sequential>& a,
const csr_rep<T,sequential>& b)
{
typedef typename csr_rep<T,sequential>::size_type size_type;
typename csr_rep<T,sequential>::const_iterator ia = a.begin(), ib = b.begin();
size_type nnz_c = 0;
for (size_type i = 0, n = a.nrow(); i < n; i++) {
std::set<size_type> y;
for (typename csr<T>::const_data_iterator ap = ia[i], last_ap = ia[i+1]; ap != last_ap; ++ap) {
size_type j = (*ap).first;
typename std::set<size_type>::iterator iter_y = y.begin();
for (typename csr<T>::const_data_iterator bp = ib[j], last_bp = ib[j+1]; bp != last_bp; ++bp) {
size_type k = (*bp).first;
iter_y = y.insert(iter_y, k);
}
}
nnz_c += y.size();
}
return nnz_c;
}
// compute the values of the sparse matrix c=a*b
template<class T>
static
void
csr_csr_mult (
const csr_rep<T,sequential>& a,
const csr_rep<T,sequential>& b,
csr_rep<T,sequential>& c)
{
typedef typename csr_rep<T,sequential>::size_type size_type;
typename csr_rep<T,sequential>::const_iterator ia = a.begin(), ib = b.begin();
typename csr_rep<T,sequential>::iterator ic = c.begin();
for (size_type i = 0, n = a.nrow(); i < n; i++) {
std::map<size_type,T> y;
for (typename csr<T>::const_data_iterator ap = ia[i], last_ap = ia[i+1]; ap != last_ap; ++ap) {
size_type j = (*ap).first;
const T& a_ij = (*ap).second;
typename std::map<size_type,T>::iterator prev_y = y.begin();
for (typename csr<T>::const_data_iterator bp = ib[j], last_bp = ib[j+1]; bp != last_bp; ++bp) {
size_type k = (*bp).first;
const T& b_jk = (*bp).second;
T value = a_ij*b_jk;
typename std::map<size_type,T>::iterator curr_y = y.find(k);
if (curr_y == y.end()) {
// create a new (i,k,value) entry in matrix C
y.insert (prev_y, std::pair<size_type,T>(k, value));
} else {
// increment an existing (i,k,value) entry in matrix C
(*curr_y).second += value;
prev_y = curr_y;
}
}
}
ic[i+1] = ic[i] + y.size();
// copy the y temporary into the i-th row of C:
typename std::map<size_type,T>::const_iterator iter_y = y.begin();
for (typename csr<T>::data_iterator cp = ic[i], last_cp = ic[i+1]; cp != last_cp; ++cp, ++iter_y) {
*cp = *iter_y;
}
}
// propagate flags ; cannot stat yet for symmetry and positive_definite
c.set_pattern_dimension (std::max(a.pattern_dimension(), b.pattern_dimension()));
}
template<class T>
void
csr_rep<T,sequential>::assign_mult (
const csr_rep<T,sequential>& a,
const csr_rep<T,sequential>& b)
{
size_type nnz_c = csr_csr_mult_size (a, b);
resize (a.row_ownership(), b.col_ownership(), nnz_c);
csr_csr_mult (a, b, *this);
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#ifndef _RHEOLEF_USE_HEAP_ALLOCATOR
#define _RHEOLEF_instanciate_class(T) \
template void csr_rep<T,sequential>::build_from_asr (const asr<T,sequential,std::allocator<T> >&);
#else // _RHEOLEF_USE_HEAP_ALLOCATOR
#define _RHEOLEF_instanciate_class(T) \
template void csr_rep<T,sequential>::build_from_asr (const asr<T,sequential,std::allocator<T> >&); \
template void csr_rep<T,sequential>::build_from_asr (const asr<T,sequential,heap_allocator<T> >&);
#endif // _RHEOLEF_USE_HEAP_ALLOCATOR
#define _RHEOLEF_instanciate(T) \
template class csr_rep<T,sequential>; \
template void csr_rep<T,sequential>::assign_add ( \
const csr_rep<T,sequential>&, \
const csr_rep<T,sequential>&, \
std::plus<T>); \
template void csr_rep<T,sequential>::assign_add ( \
const csr_rep<T,sequential>&, \
const csr_rep<T,sequential>&, \
std::minus<T>); \
_RHEOLEF_instanciate_class(T)
_RHEOLEF_instanciate(Float)
} // namespace rheolef
size_type
field::size_type size_type
Definition branch.cc:430
Float
see the Float page for the full documentation
rheolef::asr
Definition asr.h:49
rheolef::csr_rep< T, sequential >::const_data_iterator
vector_of_iterator< pair_type >::const_value_type const_data_iterator
Definition csr.h:94
rheolef::csr_rep< T, sequential >::iterator
vector_of_iterator< pair_type >::iterator iterator
Definition csr.h:90
rheolef::csr_rep< T, sequential >::pair_type
std::pair< size_type, T > pair_type
Definition csr.h:89
rheolef::csr_rep< T, sequential >::size_type
std::vector< T >::size_type size_type
Definition csr.h:86
rheolef::csr_rep< T, sequential >::data_iterator
vector_of_iterator< pair_type >::value_type data_iterator
Definition csr.h:93
rheolef::csr_rep< T, sequential >::const_iterator
vector_of_iterator< pair_type >::const_iterator const_iterator
Definition csr.h:91
rheolef::csr_rep
Definition csr.h:81
rheolef::csr
see the csr page for the full documentation
Definition csr.h:317
rheolef::distributor
see the distributor page for the full documentation
Definition distributor.h:69
rheolef::distributor::resize
void resize(size_type dis_size=0, const communicator_type &c=communicator_type(), size_type loc_size=decide)
Definition distributor.cc:30
rheolef::distributor::size
size_type size(size_type iproc) const
Definition distributor.h:170
rheolef::distributor::decide
static const size_type decide
Definition distributor.h:83
rheolef::iorheo::flag_type
std::bitset< last > flag_type
Definition iorheo.h:440
rheolef::iorheo::format_field
static flag_type format_field
Definition iorheo.h:445
rheolef::iorheo::matlab
@ matlab
Definition iorheo.h:386
rheolef::iorheo::sparse_matlab
@ sparse_matlab
Definition iorheo.h:387
rheolef::odiststream
odiststream: see the diststream page for the full documentation
Definition diststream.h:137
rheolef::odiststream::os
std::ostream & os()
Definition diststream.h:247
rheolef::vec
see the vec page for the full documentation
Definition vec.h:79
rheolef::vector_of_iterator
Definition vector_of_iterator.h:35
rheolef::vector_of_iterator::operator[]
const_reference operator[](size_type i) const
Definition vector_of_iterator.h:169
rheolef::vector_of_iterator::begin
iterator begin()
Definition vector_of_iterator.h:110
_RHEOLEF_instanciate
#define _RHEOLEF_instanciate(T)
Definition csr_seq.cc:491
T
Expr1::float_type T
Definition field_expr.h:230
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
mkgeo_ball.b
int b
Definition mkgeo_ball.sh:152
mkgeo_ball.c
int c
Definition mkgeo_ball.sh:153
mkgeo_ball.a
int a
Definition mkgeo_ball.sh:151
rheolef
This file is part of Rheolef.
Definition compiler_eigen.h:39
rheolef::csr_amux
void csr_amux(InputIterator ia, InputIterator last_ia, InputRandomAcessIterator x, SetOperator set_op, OutputIterator y)
Definition csr_amux.h:77
rheolef::put
void put(std::ostream &out, std::string name, const tiny_matrix< T > &a)
Definition tiny_lu.h:155
rheolef::csr_csr_mult_size
csr_rep< T, sequential >::size_type csr_csr_mult_size(const csr_rep< T, sequential > &a, const csr_rep< T, sequential > &b)
Definition csr_seq.cc:405
rheolef::put_matrix_market
void put_matrix_market(std::ostream &out, const Eigen::SparseMatrix< T, Eigen::RowMajor > &a)
Definition eigen_util.h:78
rheolef::csr_cumul_trans_mult
void csr_cumul_trans_mult(InputIterator1 ia, InputIterator1 last_ia, InputIterator3 x, SetOperator set_op, RandomAccessMutableIterator y)
Definition csr_cumul_trans_mult.h:67
rheolef::asr_to_csr
OutputPtrIterator asr_to_csr(InputPtrIterator iter_ptr_a, InputPtrIterator last_ptr_a, Predicate pred, Operation op, OutputPtrIterator iter_ptr_b, OutputDataIterator iter_data_b)
Definition asr_to_csr.h:70
std
STL namespace.
d
size_t d
Definition neumann-laplace-lambda.cc:28
lambda
Definition yield_slip_circle.h:34
p
Definition sphere.icc:25
rheolef::always_true
Definition msg_util.h:46
rheolef::pair_identity
Definition msg_util.h:89