DUNE PDELab (git)

sparse.hh
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
14#include <dune/pdelab/backend/interface.hh>
15#include <dune/pdelab/backend/common/uncachedmatrixview.hh>
16#include <dune/pdelab/backend/simple/descriptors.hh>
17
18namespace Dune {
19 namespace PDELab {
20 namespace Simple {
21
22 class SparseMatrixPattern
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>
43 struct SparseMatrixData
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>
77 class SparseMatrixContainer
78 : public Backend::impl::Wrapper<SparseMatrixData<C,ET,I> >
79 {
80
81 public:
82
83 typedef SparseMatrixData<C,ET,I> Container;
84
85 private:
86
87 friend Backend::impl::Wrapper<Container>;
88
89 public:
90
91 typedef ET ElementType;
92
93 typedef ElementType field_type;
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>
104 using LocalView = UncachedMatrixView<SparseMatrixContainer,RowCache,ColCache>;
105
106 template<typename RowCache, typename ColCache>
107 using ConstLocalView = ConstUncachedMatrixView<const SparseMatrixContainer,RowCache,ColCache>;
108
109 typedef SparseMatrixPattern Pattern;
110
111 template<typename GO>
112 SparseMatrixContainer(const GO& go)
113 : _container(std::make_shared<Container>())
114 {
115 allocate_matrix(_container, go, ElementType(0));
116 }
117
118 template<typename GO>
119 SparseMatrixContainer(const GO& go, const ElementType& e)
120 : _container(std::make_shared<Container>())
121 {
122 allocate_matrix(_container, go, e);
123 }
124
126 explicit SparseMatrixContainer(Backend::unattached_container = Backend::unattached_container())
127 {}
128
130 explicit SparseMatrixContainer(Backend::attached_container)
131 : _container(std::make_shared<Container>())
132 {}
133
134 SparseMatrixContainer(const SparseMatrixContainer& rhs)
135 : _container(std::make_shared<Container>(*(rhs._container)))
136 {}
137
138 SparseMatrixContainer& operator=(const SparseMatrixContainer& rhs)
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
183 SparseMatrixContainer& operator=(const ElementType& e)
184 {
185 std::fill(_container->_data.begin(),_container->_data.end(),e);
186 return *this;
187 }
188
189 SparseMatrixContainer& operator*=(const ElementType& e)
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
243 Container& base()
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
Various tags for influencing backend behavior.
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
STL namespace.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 24, 23:30, 2024)