dune-pdelab  2.4-dev
sparse.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil -*-
2 #ifndef DUNE_PDELAB_BACKEND_SIMPLE_SPARSE_HH
3 #define DUNE_PDELAB_BACKEND_SIMPLE_SPARSE_HH
4 
5 #include <vector>
6 #include <algorithm>
7 #include <functional>
8 #include <numeric>
9 #include <memory>
10 #include <unordered_set>
11 
12 #include <dune/common/typetraits.hh>
17 
18 namespace Dune {
19  namespace PDELab {
20  namespace simple {
21 
22  template<typename _RowOrdering, typename _ColOrdering>
24  : public std::vector< std::unordered_set<std::size_t> >
25  {
26 
27  public:
28 
29  typedef _RowOrdering RowOrdering;
30  typedef _ColOrdering ColOrdering;
31 
32  typedef std::unordered_set<std::size_t> col_type;
33 
34  template<typename RI, typename CI>
35  void add_link(const RI& ri, const CI& ci)
36  {
37  this->resize(_row_ordering.blockCount());
38  (*this)[ri.back()].insert(ci.back());
39  }
40 
41  SparseMatrixPattern(const RowOrdering& row_ordering, const ColOrdering& col_ordering)
42  : _row_ordering(row_ordering)
43  , _col_ordering(col_ordering)
44  {}
45 
46  private:
47 
48  const RowOrdering& _row_ordering;
49  const ColOrdering& _col_ordering;
50 
51  };
52 
53  template<template<typename> class C, typename ET, typename I>
55  {
56  typedef ET ElementType;
57  typedef I index_type;
58  typedef std::size_t size_type;
59  std::size_t _rows;
60  std::size_t _cols;
61  std::size_t _non_zeros;
62  C<ElementType> _data;
63  C<index_type> _colindex;
64  C<index_type> _rowoffset;
65  };
66 
87  template<typename GFSV, typename GFSU, template<typename> class C, typename ET, typename I>
89  : public Backend::impl::Wrapper<SparseMatrixData<C,ET,I> >
90  {
91 
92  public:
93 
95 
96  private:
97 
98  friend Backend::impl::Wrapper<Container>;
99 
100  public:
101 
102  typedef ET ElementType;
103 
104  typedef ElementType field_type;
105  typedef typename Container::size_type size_type;
106  typedef I index_type;
107 
109  typedef GFSV TestGridFunctionSpace;
110 
111  typedef typename GFSV::Ordering::Traits::ContainerIndex RowIndex;
112  typedef typename GFSU::Ordering::Traits::ContainerIndex ColIndex;
113 
114  template<typename RowCache, typename ColCache>
116 
117  template<typename RowCache, typename ColCache>
119 
120  typedef OrderingBase<
121  typename GFSV::Ordering::Traits::DOFIndex,
122  typename GFSV::Ordering::Traits::ContainerIndex
124 
125  typedef OrderingBase<
126  typename GFSU::Ordering::Traits::DOFIndex,
127  typename GFSU::Ordering::Traits::ContainerIndex
129 
131 
132  template<typename GO>
133  SparseMatrixContainer(const GO& go)
134  : _container(std::make_shared<Container>())
135  {
137  }
138 
139  template<typename GO>
140  SparseMatrixContainer(const GO& go, const ElementType& e)
141  : _container(std::make_shared<Container>())
142  {
143  allocate_matrix(_container, go, e);
144  }
145 
148  {}
149 
152  : _container(std::make_shared<Container>())
153  {}
154 
156  : _container(std::make_shared<Container>(*(rhs._container)))
157  {}
158 
160  {
161  if (this == &rhs)
162  return *this;
163  if (attached())
164  {
165  (*_container) = (*(rhs._container));
166  }
167  else
168  {
169  _container = std::make_shared<Container>(*(rhs._container));
170  }
171  return *this;
172  }
173 
174  void detach()
175  {
176  _container.reset();
177  }
178 
179  void attach(std::shared_ptr<Container> container)
180  {
181  _container = container;
182  }
183 
184  bool attached() const
185  {
186  return bool(_container);
187  }
188 
189  const std::shared_ptr<Container>& storage() const
190  {
191  return _container;
192  }
193 
194  size_type N() const
195  {
196  return _container->_rows;
197  }
198 
199  size_type M() const
200  {
201  return _container->_cols;
202  }
203 
204  SparseMatrixContainer& operator=(const ElementType& e)
205  {
206  std::fill(_container->_data.begin(),_container->_data.end(),e);
207  return *this;
208  }
209 
210  SparseMatrixContainer& operator*=(const ElementType& e)
211  {
212  using namespace std::placeholders;
213  std::transform(_container->_data.begin(),_container->_data.end(),_container->_data.begin(),std::bind(std::multiplies<ET>(),e,_1));
214  return *this;
215  }
216 
217  template<typename V>
218  void mv(const V& x, V& y) const
219  {
220  assert(y.N() == N());
221  assert(x.N() == M());
222  for (std::size_t r = 0; r < N(); ++r)
223  {
224  y.base()[r] = sparse_inner_product(r,x);
225  }
226  }
227 
228  template<typename V>
229  void usmv(const ElementType alpha, const V& x, V& y) const
230  {
231  assert(y.N() == N());
232  assert(x.N() == M());
233  for (std::size_t r = 0; r < N(); ++r)
234  {
235  y.base()[r] += alpha * sparse_inner_product(r,x);
236  }
237  }
238 
239  ElementType& operator()(const RowIndex& ri, const ColIndex& ci)
240  {
241  // entries are in ascending order
242  auto begin = _container->_colindex.begin() + _container->_rowoffset[ri[0]];
243  auto end = _container->_colindex.begin() + _container->_rowoffset[ri[0]+1];
244  auto it = std::lower_bound(begin,end,ci[0]);
245  assert (it != end);
246  return _container->_data[it - _container->_colindex.begin()];
247  }
248 
249  const ElementType& operator()(const RowIndex& ri, const ColIndex& ci) const
250  {
251  // entries are in ascending order
252  auto begin = _container->_colindex.begin() + _container->_rowoffset[ri[0]];
253  auto end = _container->_colindex.begin() + _container->_rowoffset[ri[0]+1];
254  auto it = std::lower_bound(begin,end,ci[0]);
255  assert(it != end);
256  return _container->_data[it - _container->_colindex.begin()];
257  }
258 
259  const Container& base() const
260  {
261  return *_container;
262  }
263 
264  Container& base()
265  {
266  return *_container;
267  }
268 
269  private:
270 
271  const Container& native() const
272  {
273  return *_container;
274  }
275 
276  Container& native()
277  {
278  return *_container;
279  }
280 
281  public:
282 
283  void flush()
284  {}
285 
286  void finalize()
287  {}
288 
289  void clear_row(const RowIndex& ri, const ElementType& diagonal_entry)
290  {
291  std::fill(
292  _container->_data.begin() + _container->_rowoffset[ri[0]],
293  _container->_data.begin() + _container->_rowoffset[ri[0]+1], ElementType(0));
294  (*this)(ri,ri) = diagonal_entry;
295  }
296 
297  protected:
298  template<typename GO>
299  static void allocate_matrix(std::shared_ptr<Container> & c, const GO & go, const ElementType& e)
300  {
301  typedef typename Pattern::col_type col_type;
302  Pattern pattern(go.testGridFunctionSpace().ordering(),go.trialGridFunctionSpace().ordering());
303  go.fill_pattern(pattern);
304 
305  c->_rows = go.testGridFunctionSpace().size();
306  c->_cols = go.trialGridFunctionSpace().size();
307  // compute row offsets
308  c->_rowoffset.resize(c->_rows+1);
309  size_type offset = 0;
310  auto calc_offset = [=](const col_type & entry) mutable -> size_t { offset += entry.size(); return offset; };
311  std::transform(pattern.begin(), pattern.end(),
312  c->_rowoffset.begin()+1,
313  calc_offset);
314  // compute non-zeros
315  c->_non_zeros = c->_rowoffset.back();
316  // allocate col/data vectors
317  c->_data.resize(c->_non_zeros, e);
318  c->_colindex.resize(c->_non_zeros);
319  // copy pattern
320  auto colit = c->_colindex.begin();
321  c->_rowoffset[0] = 0;
322  for (auto & row : pattern)
323  {
324  auto last = std::copy(row.begin(),row.end(),colit);
325  std::sort(colit, last);
326  colit = last;
327  }
328  }
329 
330  template<typename V>
331  ElementType sparse_inner_product (std::size_t row, const V & x) const {
332  std::size_t begin = _container->_rowoffset[row];
333  std::size_t end = _container->_rowoffset[row+1];
334  auto accu = [](const ElementType & a, const ElementType & b) -> ElementType { return a+b; };
335  auto mult = [=](const ElementType & a, const index_type & i) -> ElementType { return a * x.base()[i]; };
336  return std::inner_product(
337  &_container->_data[begin],
338  &_container->_data[end],
339  &_container->_colindex[begin],
340  ElementType(0),
341  accu, mult
342  );
343  }
344 
345  std::shared_ptr< Container > _container;
346  };
347 
348  } // namespace simple
349  } // namespace PDELab
350 } // namespace Dune
351 
352 #endif // DUNE_PDELAB_BACKEND_SIMPLE_SPARSE_HH
const std::size_t offset
Definition: localfunctionspace.hh:74
void detach()
Definition: sparse.hh:174
OrderingBase< typename GFSV::Ordering::Traits::DOFIndex, typename GFSV::Ordering::Traits::ContainerIndex > RowOrdering
Definition: sparse.hh:123
SparseMatrixContainer(const GO &go, const ElementType &e)
Definition: sparse.hh:140
Container & base()
Definition: sparse.hh:264
SparseMatrixContainer & operator=(const ElementType &e)
Definition: sparse.hh:204
SparseMatrixContainer(const SparseMatrixContainer &rhs)
Definition: sparse.hh:155
Container::size_type size_type
Definition: sparse.hh:105
OrderingBase< typename GFSU::Ordering::Traits::DOFIndex, typename GFSU::Ordering::Traits::ContainerIndex > ColOrdering
Definition: sparse.hh:128
Various tags for influencing backend behavior.
const Container & base() const
Definition: sparse.hh:259
SparseMatrixPattern< RowOrdering, ColOrdering > Pattern
Definition: sparse.hh:130
void clear_row(const RowIndex &ri, const ElementType &diagonal_entry)
Definition: sparse.hh:289
const E & e
Definition: interpolate.hh:172
STL namespace.
ElementType sparse_inner_product(std::size_t row, const V &x) const
Definition: sparse.hh:331
Definition: uncachedmatrixview.hh:158
static void allocate_matrix(std::shared_ptr< Container > &c, const GO &go, const ElementType &e)
Definition: sparse.hh:299
C< ElementType > _data
Definition: sparse.hh:62
void mv(const V &x, V &y) const
Definition: sparse.hh:218
I index_type
Definition: sparse.hh:106
SparseMatrixContainer(Backend::attached_container)
Creates an SparseMatrixContainer with an empty underlying ISTL matrix.
Definition: sparse.hh:151
GFSU::Ordering::Traits::ContainerIndex ColIndex
Definition: sparse.hh:112
void attach(std::shared_ptr< Container > container)
Definition: sparse.hh:179
std::size_t size_type
Definition: sparse.hh:58
C< index_type > _colindex
Definition: sparse.hh:63
std::unordered_set< std::size_t > col_type
Definition: sparse.hh:32
_RowOrdering RowOrdering
Definition: sparse.hh:29
ElementType field_type
Definition: sparse.hh:104
ElementType & operator()(const RowIndex &ri, const ColIndex &ci)
Definition: sparse.hh:239
void usmv(const ElementType alpha, const V &x, V &y) const
Definition: sparse.hh:229
std::size_t _non_zeros
Definition: sparse.hh:61
void flush()
Definition: sparse.hh:283
Tag for requesting a vector or matrix container without a pre-attached underlying object...
Definition: backend/common/tags.hh:23
Tag for requesting a vector or matrix container with a pre-attached underlying object.
Definition: backend/common/tags.hh:27
void finalize()
Definition: sparse.hh:286
std::size_t _rows
Definition: sparse.hh:59
ET ElementType
Definition: sparse.hh:102
Definition: uncachedmatrixview.hh:13
void add_link(const RI &ri, const CI &ci)
Definition: sparse.hh:35
SparseMatrixContainer & operator=(const SparseMatrixContainer &rhs)
Definition: sparse.hh:159
GFSV TestGridFunctionSpace
Definition: sparse.hh:109
Definition: adaptivity.hh:27
std::size_t _cols
Definition: sparse.hh:60
SparseMatrixContainer(Backend::unattached_container=Backend::unattached_container())
Creates an SparseMatrixContainer without allocating an underlying ISTL matrix.
Definition: sparse.hh:147
C< index_type > _rowoffset
Definition: sparse.hh:64
I index_type
Definition: sparse.hh:57
bool attached() const
Definition: sparse.hh:184
GFSV::Ordering::Traits::ContainerIndex RowIndex
Definition: sparse.hh:111
SparseMatrixData< C, ET, I > Container
Definition: sparse.hh:94
Definition: orderingbase.hh:21
std::shared_ptr< Container > _container
Definition: sparse.hh:345
const ElementType & operator()(const RowIndex &ri, const ColIndex &ci) const
Definition: sparse.hh:249
_ColOrdering ColOrdering
Definition: sparse.hh:30
SparseMatrixContainer(const GO &go)
Definition: sparse.hh:133
size_type M() const
Definition: sparse.hh:199
SparseMatrixPattern(const RowOrdering &row_ordering, const ColOrdering &col_ordering)
Definition: sparse.hh:41
SparseMatrixContainer & operator*=(const ElementType &e)
Definition: sparse.hh:210
const std::shared_ptr< Container > & storage() const
Definition: sparse.hh:189
ET ElementType
Definition: sparse.hh:56
size_type N() const
Definition: sparse.hh:194
GFSU TrialGridFunctionSpace
Definition: sparse.hh:108