dune-pdelab 2.7-git
Loading...
Searching...
No Matches
/construction/math/dune-pdelab/dune-pdelab-20c7085389d3eb4f8ca99e1bc60f7fa6036536c8/dune/pdelab/backend/simple/sparse.hh

Simple backend for CSR matrices.

Simple backend for CSR matrices

Matrix stored as a: nnz Number of nonzero elements data CSR format data array of the matrix indices CSR format index array of the matrix rowptr CSR format index pointer array of the matrix

Consider the following 3x3 matrix [1, 0, 2] [0, 0, 3] [4, 5, 6] the data would be stored as nnz=6 data=[1 2 3 4 5 6] indices=[0 2 2 0 1 2] rowptr=[0 2 3 6]

// -*- tab-width: 4; indent-tabs-mode: nil -*-
#ifndef DUNE_PDELAB_BACKEND_SIMPLE_SPARSE_HH
#define DUNE_PDELAB_BACKEND_SIMPLE_SPARSE_HH
#include <vector>
#include <algorithm>
#include <functional>
#include <numeric>
#include <memory>
#include <unordered_set>
#include <dune/common/typetraits.hh>
namespace Dune {
namespace PDELab {
namespace Simple {
class SparseMatrixPattern
: public std::vector< std::unordered_set<std::size_t> >
{
public:
typedef std::unordered_set<std::size_t> col_type;
template<typename RI, typename CI>
void add_link(const RI& ri, const CI& ci)
{
(*this)[ri.back()].insert(ci.back());
}
SparseMatrixPattern(std::size_t rows)
: std::vector< std::unordered_set<std::size_t> >(rows)
{}
};
template<template<typename> class C, typename ET, typename I>
struct SparseMatrixData
{
typedef ET ElementType;
typedef I index_type;
typedef std::size_t size_type;
std::size_t _rows;
std::size_t _cols;
std::size_t _non_zeros;
C<ElementType> _data;
C<index_type> _colindex;
C<index_type> _rowoffset;
};
template<typename GFSV, typename GFSU, template<typename> class C, typename ET, typename I>
class SparseMatrixContainer
: public Backend::impl::Wrapper<SparseMatrixData<C,ET,I> >
{
public:
typedef SparseMatrixData<C,ET,I> Container;
private:
friend Backend::impl::Wrapper<Container>;
public:
typedef ET ElementType;
typedef typename Container::size_type size_type;
typedef I index_type;
typedef GFSU TrialGridFunctionSpace;
typedef GFSV TestGridFunctionSpace;
typedef typename GFSV::Ordering::Traits::ContainerIndex RowIndex;
typedef typename GFSU::Ordering::Traits::ContainerIndex ColIndex;
template<typename RowCache, typename ColCache>
using LocalView = UncachedMatrixView<SparseMatrixContainer,RowCache,ColCache>;
template<typename RowCache, typename ColCache>
using ConstLocalView = ConstUncachedMatrixView<const SparseMatrixContainer,RowCache,ColCache>;
typedef SparseMatrixPattern Pattern;
template<typename GO>
SparseMatrixContainer(const GO& go)
: _container(std::make_shared<Container>())
{
}
template<typename GO>
SparseMatrixContainer(const GO& go, const ElementType& e)
: _container(std::make_shared<Container>())
{
}
explicit SparseMatrixContainer(Backend::unattached_container = Backend::unattached_container())
{}
explicit SparseMatrixContainer(Backend::attached_container)
: _container(std::make_shared<Container>())
{}
SparseMatrixContainer(const SparseMatrixContainer& rhs)
: _container(std::make_shared<Container>(*(rhs._container)))
{}
SparseMatrixContainer& operator=(const SparseMatrixContainer& rhs)
{
if (this == &rhs)
return *this;
if (attached())
{
(*_container) = (*(rhs._container));
}
else
{
_container = std::make_shared<Container>(*(rhs._container));
}
return *this;
}
void detach()
{
_container.reset();
}
void attach(std::shared_ptr<Container> container)
{
_container = container;
}
bool attached() const
{
return bool(_container);
}
const std::shared_ptr<Container>& storage() const
{
return _container;
}
size_type N() const
{
return _container->_rows;
}
size_type M() const
{
return _container->_cols;
}
SparseMatrixContainer& operator=(const ElementType& e)
{
std::fill(_container->_data.begin(),_container->_data.end(),e);
return *this;
}
SparseMatrixContainer& operator*=(const ElementType& e)
{
using namespace std::placeholders;
std::transform(_container->_data.begin(),_container->_data.end(),_container->_data.begin(),std::bind(std::multiplies<ET>(),e,_1));
return *this;
}
template<typename V>
void mv(const V& x, V& y) const
{
assert(y.N() == N());
assert(x.N() == M());
for (std::size_t r = 0; r < N(); ++r)
{
y.base()[r] = sparse_inner_product(r,x);
}
}
template<typename V>
void usmv(const ElementType alpha, const V& x, V& y) const
{
assert(y.N() == N());
assert(x.N() == M());
for (std::size_t r = 0; r < N(); ++r)
{
y.base()[r] += alpha * sparse_inner_product(r,x);
}
}
ElementType& operator()(const RowIndex& ri, const ColIndex& ci)
{
// entries are in ascending order
auto begin = _container->_colindex.begin() + _container->_rowoffset[ri[0]];
auto end = _container->_colindex.begin() + _container->_rowoffset[ri[0]+1];
auto it = std::lower_bound(begin,end,ci[0]);
assert (it != end);
return _container->_data[it - _container->_colindex.begin()];
}
const ElementType& operator()(const RowIndex& ri, const ColIndex& ci) const
{
// entries are in ascending order
auto begin = _container->_colindex.begin() + _container->_rowoffset[ri[0]];
auto end = _container->_colindex.begin() + _container->_rowoffset[ri[0]+1];
auto it = std::lower_bound(begin,end,ci[0]);
assert(it != end);
return _container->_data[it - _container->_colindex.begin()];
}
const Container& base() const
{
return *_container;
}
{
return *_container;
}
private:
const Container& native() const
{
return *_container;
}
Container& native()
{
return *_container;
}
public:
void flush()
{}
void finalize()
{}
void clear_row(const RowIndex& ri, const ElementType& diagonal_entry)
{
std::fill(
_container->_data.begin() + _container->_rowoffset[ri[0]],
_container->_data.begin() + _container->_rowoffset[ri[0]+1], ElementType(0));
(*this)(ri,ri) = diagonal_entry;
}
void clear_row_block(const RowIndex& ri, const ElementType& diagonal_entry)
{
std::fill(
_container->_data.begin() + _container->_rowoffset[ri[0]],
_container->_data.begin() + _container->_rowoffset[ri[0]+1], ElementType(0));
(*this)(ri,ri) = diagonal_entry;
}
protected:
template<typename GO>
static void allocate_matrix(std::shared_ptr<Container> & c, const GO & go, const ElementType& e)
{
typedef typename Pattern::col_type col_type;
Pattern pattern(go.testGridFunctionSpace().ordering().blockCount());
go.fill_pattern(pattern);
c->_rows = go.testGridFunctionSpace().size();
c->_cols = go.trialGridFunctionSpace().size();
// compute row offsets
c->_rowoffset.resize(c->_rows+1);
auto calc_offset = [=](const col_type & entry) mutable -> size_t { offset += entry.size(); return offset; };
std::transform(pattern.begin(), pattern.end(),
c->_rowoffset.begin()+1,
calc_offset);
// compute non-zeros
c->_non_zeros = c->_rowoffset.back();
// allocate col/data vectors
c->_data.resize(c->_non_zeros, e);
c->_colindex.resize(c->_non_zeros);
// copy pattern
auto colit = c->_colindex.begin();
c->_rowoffset[0] = 0;
for (auto & row : pattern)
{
auto last = std::copy(row.begin(),row.end(),colit);
std::sort(colit, last);
colit = last;
}
}
template<typename V>
ElementType sparse_inner_product (std::size_t row, const V & x) const {
std::size_t begin = _container->_rowoffset[row];
std::size_t end = _container->_rowoffset[row+1];
auto accu = [](const ElementType & a, const ElementType & b) -> ElementType { return a+b; };
auto mult = [=](const ElementType & a, const index_type & i) -> ElementType { return a * x.base()[i]; };
return std::inner_product(
&_container->_data[begin],
&_container->_data[end],
&_container->_colindex[begin],
accu, mult
);
}
std::shared_ptr< Container > _container;
};
} // namespace Simple
} // namespace PDELab
} // namespace Dune
#endif // DUNE_PDELAB_BACKEND_SIMPLE_SPARSE_HH
const std::size_t offset
Definition: localfunctionspace.hh:75
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
std::unordered_set< std::size_t > col_type
Definition: sparse.hh:28
void add_link(const RI &ri, const CI &ci)
Definition: sparse.hh:31
ET ElementType
Definition: sparse.hh:45
C< index_type > _rowoffset
Definition: sparse.hh:53
C< ElementType > _data
Definition: sparse.hh:51
std::size_t _rows
Definition: sparse.hh:48
I index_type
Definition: sparse.hh:46
std::size_t size_type
Definition: sparse.hh:47
C< index_type > _colindex
Definition: sparse.hh:52
std::size_t _non_zeros
Definition: sparse.hh:50
std::size_t _cols
Definition: sparse.hh:49
size_type N() const
Definition: sparse.hh:173
Container::size_type size_type
Definition: sparse.hh:94
void clear_row_block(const RowIndex &ri, const ElementType &diagonal_entry)
Definition: sparse.hh:276
void flush()
Definition: sparse.hh:262
void mv(const V &x, V &y) const
Definition: sparse.hh:197
ElementType field_type
Definition: sparse.hh:93
ET ElementType
Definition: sparse.hh:91
UncachedMatrixView< SparseMatrixContainer, RowCache, ColCache > LocalView
Definition: sparse.hh:104
SparseMatrixContainer & operator=(const SparseMatrixContainer &rhs)
Definition: sparse.hh:138
void clear_row(const RowIndex &ri, const ElementType &diagonal_entry)
Definition: sparse.hh:268
void attach(std::shared_ptr< Container > container)
Definition: sparse.hh:158
SparseMatrixPattern Pattern
Definition: sparse.hh:109
std::shared_ptr< Container > _container
Definition: sparse.hh:332
void detach()
Definition: sparse.hh:153
const Container & base() const
Definition: sparse.hh:238
ElementType & operator()(const RowIndex &ri, const ColIndex &ci)
Definition: sparse.hh:218
const std::shared_ptr< Container > & storage() const
Definition: sparse.hh:168
GFSV::Ordering::Traits::ContainerIndex RowIndex
Definition: sparse.hh:100
ElementType sparse_inner_product(std::size_t row, const V &x) const
Definition: sparse.hh:318
bool attached() const
Definition: sparse.hh:163
GFSU::Ordering::Traits::ContainerIndex ColIndex
Definition: sparse.hh:101
GFSU TrialGridFunctionSpace
Definition: sparse.hh:97
SparseMatrixData< C, ET, I > Container
Definition: sparse.hh:83
void finalize()
Definition: sparse.hh:265
static void allocate_matrix(std::shared_ptr< Container > &c, const GO &go, const ElementType &e)
Definition: sparse.hh:286
void usmv(const ElementType alpha, const V &x, V &y) const
Definition: sparse.hh:208
ConstUncachedMatrixView< const SparseMatrixContainer, RowCache, ColCache > ConstLocalView
Definition: sparse.hh:107
GFSV TestGridFunctionSpace
Definition: sparse.hh:98
size_type M() const
Definition: sparse.hh:178
SparseMatrixContainer & operator*=(const ElementType &e)
Definition: sparse.hh:189
Various tags for influencing backend behavior.