dune-pdelab 2.7-git
Loading...
Searching...
No Matches
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
18namespace 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
95 typedef I index_type;
96
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>
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
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 std::size_t offset
Definition: localfunctionspace.hh:75
STL namespace.
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
Container::size_type size_type
Definition: sparse.hh:94
Container & base()
Definition: sparse.hh:243
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
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
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 Container & base() const
Definition: sparse.hh:238
const ElementType & operator()(const RowIndex &ri, const ColIndex &ci) const
Definition: sparse.hh:228
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
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
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
SparseMatrixContainer & operator*=(const ElementType &e)
Definition: sparse.hh:189
Various tags for influencing backend behavior.