dune-pdelab  2.7-git
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 
23  : public std::vector< std::unordered_set<std::size_t> >
24  {
25 
26  public:
27 
28  typedef std::unordered_set<std::size_t> col_type;
29 
30  template<typename RI, typename CI>
31  void add_link(const RI& ri, const CI& ci)
32  {
33  (*this)[ri.back()].insert(ci.back());
34  }
35 
36  SparseMatrixPattern(std::size_t rows)
37  : std::vector< std::unordered_set<std::size_t> >(rows)
38  {}
39 
40  };
41 
42  template<template<typename> class C, typename ET, typename I>
44  {
45  typedef ET ElementType;
46  typedef I index_type;
47  typedef std::size_t size_type;
48  std::size_t _rows;
49  std::size_t _cols;
50  std::size_t _non_zeros;
51  C<ElementType> _data;
52  C<index_type> _colindex;
53  C<index_type> _rowoffset;
54  };
55 
76  template<typename GFSV, typename GFSU, template<typename> class C, typename ET, typename I>
78  : public Backend::impl::Wrapper<SparseMatrixData<C,ET,I> >
79  {
80 
81  public:
82 
84 
85  private:
86 
87  friend Backend::impl::Wrapper<Container>;
88 
89  public:
90 
91  typedef ET ElementType;
92 
94  typedef typename Container::size_type size_type;
95  typedef I index_type;
96 
97  typedef GFSU TrialGridFunctionSpace;
98  typedef GFSV TestGridFunctionSpace;
99 
100  typedef typename GFSV::Ordering::Traits::ContainerIndex RowIndex;
101  typedef typename GFSU::Ordering::Traits::ContainerIndex ColIndex;
102 
103  template<typename RowCache, typename ColCache>
105 
106  template<typename RowCache, typename ColCache>
108 
110 
111  template<typename GO>
112  SparseMatrixContainer(const GO& go)
113  : _container(std::make_shared<Container>())
114  {
116  }
117 
118  template<typename GO>
119  SparseMatrixContainer(const GO& go, const ElementType& e)
120  : _container(std::make_shared<Container>())
121  {
123  }
124 
127  {}
128 
131  : _container(std::make_shared<Container>())
132  {}
133 
135  : _container(std::make_shared<Container>(*(rhs._container)))
136  {}
137 
139  {
140  if (this == &rhs)
141  return *this;
142  if (attached())
143  {
144  (*_container) = (*(rhs._container));
145  }
146  else
147  {
148  _container = std::make_shared<Container>(*(rhs._container));
149  }
150  return *this;
151  }
152 
153  void detach()
154  {
155  _container.reset();
156  }
157 
158  void attach(std::shared_ptr<Container> container)
159  {
160  _container = container;
161  }
162 
163  bool attached() const
164  {
165  return bool(_container);
166  }
167 
168  const std::shared_ptr<Container>& storage() const
169  {
170  return _container;
171  }
172 
173  size_type N() const
174  {
175  return _container->_rows;
176  }
177 
178  size_type M() const
179  {
180  return _container->_cols;
181  }
182 
184  {
185  std::fill(_container->_data.begin(),_container->_data.end(),e);
186  return *this;
187  }
188 
190  {
191  using namespace std::placeholders;
192  std::transform(_container->_data.begin(),_container->_data.end(),_container->_data.begin(),std::bind(std::multiplies<ET>(),e,_1));
193  return *this;
194  }
195 
196  template<typename V>
197  void mv(const V& x, V& y) const
198  {
199  assert(y.N() == N());
200  assert(x.N() == M());
201  for (std::size_t r = 0; r < N(); ++r)
202  {
203  y.base()[r] = sparse_inner_product(r,x);
204  }
205  }
206 
207  template<typename V>
208  void usmv(const ElementType alpha, const V& x, V& y) const
209  {
210  assert(y.N() == N());
211  assert(x.N() == M());
212  for (std::size_t r = 0; r < N(); ++r)
213  {
214  y.base()[r] += alpha * sparse_inner_product(r,x);
215  }
216  }
217 
218  ElementType& operator()(const RowIndex& ri, const ColIndex& ci)
219  {
220  // entries are in ascending order
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]);
224  assert (it != end);
225  return _container->_data[it - _container->_colindex.begin()];
226  }
227 
228  const ElementType& operator()(const RowIndex& ri, const ColIndex& ci) const
229  {
230  // entries are in ascending order
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]);
234  assert(it != end);
235  return _container->_data[it - _container->_colindex.begin()];
236  }
237 
238  const Container& base() const
239  {
240  return *_container;
241  }
242 
244  {
245  return *_container;
246  }
247 
248  private:
249 
250  const Container& native() const
251  {
252  return *_container;
253  }
254 
255  Container& native()
256  {
257  return *_container;
258  }
259 
260  public:
261 
262  void flush()
263  {}
264 
265  void finalize()
266  {}
267 
268  void clear_row(const RowIndex& ri, const ElementType& diagonal_entry)
269  {
270  std::fill(
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;
274  }
275 
276  void clear_row_block(const RowIndex& ri, const ElementType& diagonal_entry)
277  {
278  std::fill(
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;
282  }
283 
284  protected:
285  template<typename GO>
286  static void allocate_matrix(std::shared_ptr<Container> & c, const GO & go, const ElementType& e)
287  {
288  typedef typename Pattern::col_type col_type;
289  Pattern pattern(go.testGridFunctionSpace().ordering().blockCount());
290  go.fill_pattern(pattern);
291 
292  c->_rows = go.testGridFunctionSpace().size();
293  c->_cols = go.trialGridFunctionSpace().size();
294  // compute row offsets
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,
300  calc_offset);
301  // compute non-zeros
302  c->_non_zeros = c->_rowoffset.back();
303  // allocate col/data vectors
304  c->_data.resize(c->_non_zeros, e);
305  c->_colindex.resize(c->_non_zeros);
306  // copy pattern
307  auto colit = c->_colindex.begin();
308  c->_rowoffset[0] = 0;
309  for (auto & row : pattern)
310  {
311  auto last = std::copy(row.begin(),row.end(),colit);
312  std::sort(colit, last);
313  colit = last;
314  }
315  }
316 
317  template<typename V>
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],
327  ElementType(0),
328  accu, mult
329  );
330  }
331 
332  std::shared_ptr< Container > _container;
333  };
334 
335  } // namespace Simple
336  } // namespace PDELab
337 } // namespace Dune
338 
339 #endif // DUNE_PDELAB_BACKEND_SIMPLE_SPARSE_HH
const Entity & e
Definition: localfunctionspace.hh:123
const std::size_t offset
Definition: localfunctionspace.hh:75
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
Tag for requesting a vector or matrix container without a pre-attached underlying object.
Definition: backend/common/tags.hh:24
Tag for requesting a vector or matrix container with a pre-attached underlying object.
Definition: backend/common/tags.hh:28
Definition: uncachedmatrixview.hh:13
Definition: uncachedmatrixview.hh:167
std::unordered_set< std::size_t > col_type
Definition: sparse.hh:28
SparseMatrixPattern(std::size_t rows)
Definition: sparse.hh:36
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
ElementType & operator()(const RowIndex &ri, const ColIndex &ci)
Definition: sparse.hh:218
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
SparseMatrixContainer(Backend::attached_container)
Creates an SparseMatrixContainer with an empty underlying ISTL matrix.
Definition: sparse.hh:130
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
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
SparseMatrixContainer & operator*=(const ElementType &e)
Definition: sparse.hh:189
std::shared_ptr< Container > _container
Definition: sparse.hh:332
const Container & base() const
Definition: sparse.hh:238
SparseMatrixContainer(Backend::unattached_container=Backend::unattached_container())
Creates an SparseMatrixContainer without allocating an underlying ISTL matrix.
Definition: sparse.hh:126
void detach()
Definition: sparse.hh:153
const std::shared_ptr< Container > & storage() const
Definition: sparse.hh:168
GFSV::Ordering::Traits::ContainerIndex RowIndex
Definition: sparse.hh:100
SparseMatrixContainer & operator=(const SparseMatrixContainer &rhs)
Definition: sparse.hh:138
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
const ElementType & operator()(const RowIndex &ri, const ColIndex &ci) const
Definition: sparse.hh:228
SparseMatrixContainer(const SparseMatrixContainer &rhs)
Definition: sparse.hh:134
GFSU TrialGridFunctionSpace
Definition: sparse.hh:97
SparseMatrixData< C, ET, I > Container
Definition: sparse.hh:83
SparseMatrixContainer & operator=(const ElementType &e)
Definition: sparse.hh:183
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
Container & base()
Definition: sparse.hh:243
SparseMatrixContainer(const GO &go, const ElementType &e)
Definition: sparse.hh:119
void usmv(const ElementType alpha, const V &x, V &y) const
Definition: sparse.hh:208
GFSV TestGridFunctionSpace
Definition: sparse.hh:98
SparseMatrixContainer(const GO &go)
Definition: sparse.hh:112
size_type M() const
Definition: sparse.hh:178
Various tags for influencing backend behavior.