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