2#ifndef DUNE_PDELAB_BACKEND_SIMPLE_SPARSE_HH
3#define DUNE_PDELAB_BACKEND_SIMPLE_SPARSE_HH
10#include <unordered_set>
14#include <dune/pdelab/backend/interface.hh>
15#include <dune/pdelab/backend/common/uncachedmatrixview.hh>
16#include <dune/pdelab/backend/simple/descriptors.hh>
22 class SparseMatrixPattern
23 :
public std::vector< std::unordered_set<std::size_t> >
28 typedef std::unordered_set<std::size_t> col_type;
30 template<
typename RI,
typename CI>
31 void add_link(
const RI& ri,
const CI& ci)
33 (*this)[ri.back()].insert(ci.back());
36 SparseMatrixPattern(std::size_t rows)
37 :
std::vector<
std::unordered_set<
std::size_t> >(rows)
42 template<
template<
typename>
class C,
typename ET,
typename I>
43 struct SparseMatrixData
45 typedef ET ElementType;
47 typedef std::size_t size_type;
50 std::size_t _non_zeros;
52 C<index_type> _colindex;
53 C<index_type> _rowoffset;
76 template<
typename GFSV,
typename GFSU,
template<
typename>
class C,
typename ET,
typename I>
77 class SparseMatrixContainer
78 :
public Backend::impl::Wrapper<SparseMatrixData<C,ET,I> >
83 typedef SparseMatrixData<C,ET,I> Container;
87 friend Backend::impl::Wrapper<Container>;
91 typedef ET ElementType;
93 typedef ElementType field_type;
94 typedef typename Container::size_type size_type;
97 typedef GFSU TrialGridFunctionSpace;
98 typedef GFSV TestGridFunctionSpace;
100 typedef typename GFSV::Ordering::Traits::ContainerIndex RowIndex;
101 typedef typename GFSU::Ordering::Traits::ContainerIndex ColIndex;
103 template<
typename RowCache,
typename ColCache>
104 using LocalView = UncachedMatrixView<SparseMatrixContainer,RowCache,ColCache>;
106 template<
typename RowCache,
typename ColCache>
107 using ConstLocalView = ConstUncachedMatrixView<const SparseMatrixContainer,RowCache,ColCache>;
109 typedef SparseMatrixPattern Pattern;
111 template<
typename GO>
112 SparseMatrixContainer(
const GO& go)
113 : _container(
std::make_shared<Container>())
115 allocate_matrix(_container, go, ElementType(0));
118 template<
typename GO>
119 SparseMatrixContainer(
const GO& go,
const ElementType& e)
120 : _container(
std::make_shared<Container>())
122 allocate_matrix(_container, go, e);
126 explicit SparseMatrixContainer(Backend::unattached_container = Backend::unattached_container())
130 explicit SparseMatrixContainer(Backend::attached_container)
131 : _container(
std::make_shared<Container>())
134 SparseMatrixContainer(
const SparseMatrixContainer& rhs)
135 : _container(
std::make_shared<Container>(*(rhs._container)))
138 SparseMatrixContainer& operator=(
const SparseMatrixContainer& rhs)
144 (*_container) = (*(rhs._container));
148 _container = std::make_shared<Container>(*(rhs._container));
158 void attach(std::shared_ptr<Container> container)
160 _container = container;
163 bool attached()
const
165 return bool(_container);
168 const std::shared_ptr<Container>& storage()
const
175 return _container->_rows;
180 return _container->_cols;
183 SparseMatrixContainer& operator=(
const ElementType& e)
185 std::fill(_container->_data.begin(),_container->_data.end(),e);
189 SparseMatrixContainer& operator*=(
const ElementType& e)
191 using namespace std::placeholders;
192 std::transform(_container->_data.begin(),_container->_data.end(),_container->_data.begin(),std::bind(std::multiplies<ET>(),e,
_1));
197 void mv(
const V& x, V& y)
const
199 assert(y.N() == N());
200 assert(x.N() == M());
201 for (std::size_t r = 0; r < N(); ++r)
203 y.base()[r] = sparse_inner_product(r,x);
208 void usmv(
const ElementType alpha,
const V& x, V& y)
const
210 assert(y.N() == N());
211 assert(x.N() == M());
212 for (std::size_t r = 0; r < N(); ++r)
214 y.base()[r] += alpha * sparse_inner_product(r,x);
218 ElementType& operator()(
const RowIndex& ri,
const ColIndex& ci)
221 auto begin = _container->_colindex.begin() + _container->_rowoffset[ri[0]];
222 auto end = _container->_colindex.begin() + _container->_rowoffset[ri[0]+1];
223 auto it = std::lower_bound(begin,end,ci[0]);
225 return _container->_data[it - _container->_colindex.begin()];
228 const ElementType& operator()(
const RowIndex& ri,
const ColIndex& ci)
const
231 auto begin = _container->_colindex.begin() + _container->_rowoffset[ri[0]];
232 auto end = _container->_colindex.begin() + _container->_rowoffset[ri[0]+1];
233 auto it = std::lower_bound(begin,end,ci[0]);
235 return _container->_data[it - _container->_colindex.begin()];
238 const Container& base()
const
250 const Container& native()
const
268 void clear_row(
const RowIndex& ri,
const ElementType& diagonal_entry)
271 _container->_data.begin() + _container->_rowoffset[ri[0]],
272 _container->_data.begin() + _container->_rowoffset[ri[0]+1], ElementType(0));
273 (*this)(ri,ri) = diagonal_entry;
276 void clear_row_block(
const RowIndex& ri,
const ElementType& diagonal_entry)
279 _container->_data.begin() + _container->_rowoffset[ri[0]],
280 _container->_data.begin() + _container->_rowoffset[ri[0]+1], ElementType(0));
281 (*this)(ri,ri) = diagonal_entry;
285 template<
typename GO>
286 static void allocate_matrix(std::shared_ptr<Container> & c,
const GO & go,
const ElementType& e)
288 typedef typename Pattern::col_type col_type;
289 Pattern pattern(go.testGridFunctionSpace().ordering().blockCount());
290 go.fill_pattern(pattern);
292 c->_rows = go.testGridFunctionSpace().size();
293 c->_cols = go.trialGridFunctionSpace().size();
295 c->_rowoffset.resize(c->_rows+1);
296 size_type offset = 0;
297 auto calc_offset = [=](
const col_type & entry)
mutable ->
size_t { offset += entry.size();
return offset; };
298 std::transform(pattern.begin(), pattern.end(),
299 c->_rowoffset.begin()+1,
302 c->_non_zeros = c->_rowoffset.back();
304 c->_data.resize(c->_non_zeros, e);
305 c->_colindex.resize(c->_non_zeros);
307 auto colit = c->_colindex.begin();
308 c->_rowoffset[0] = 0;
309 for (
auto & row : pattern)
311 auto last = std::copy(row.begin(),row.end(),colit);
312 std::sort(colit, last);
318 ElementType sparse_inner_product (std::size_t row,
const V & x)
const {
319 std::size_t begin = _container->_rowoffset[row];
320 std::size_t end = _container->_rowoffset[row+1];
321 auto accu = [](
const ElementType & a,
const ElementType & b) -> ElementType {
return a+b; };
322 auto mult = [=](
const ElementType & a,
const index_type & i) -> ElementType {
return a * x.base()[i]; };
323 return std::inner_product(
324 &_container->_data[begin],
325 &_container->_data[end],
326 &_container->_colindex[begin],
332 std::shared_ptr< Container > _container;
Traits for type conversions and type information.
constexpr index_constant< 1 > _1
Compile time index with value 1.
Definition: indices.hh:55
Dune namespace.
Definition: alignedallocator.hh:13