DUNE PDELab (git)

matrix.hh
1// -*- tab-width: 4; indent-tabs-mode: nil -*-
2#ifndef DUNE_PDELAB_BACKEND_SIMPLE_MATRIX_HH
3#define DUNE_PDELAB_BACKEND_SIMPLE_MATRIX_HH
4
5#include <vector>
6#include <algorithm>
7#include <functional>
8#include <memory>
9#include <numeric>
10
13#include <dune/pdelab/backend/interface.hh>
14#include <dune/pdelab/backend/common/uncachedmatrixview.hh>
15#include <dune/pdelab/backend/simple/descriptors.hh>
16
17namespace Dune {
18 namespace PDELab {
19 namespace Simple {
20
21 template<typename GFSV, typename GFSU, typename C>
22 class MatrixContainer
23 : public Backend::impl::Wrapper<C>
24 {
25
26 friend Backend::impl::Wrapper<C>;
27
28 public:
29
30 typedef C Container;
31 typedef typename Container::value_type ElementType;
32 typedef ElementType E;
33
34 typedef ElementType field_type;
35 typedef typename Container::size_type size_type;
36
37 typedef GFSU TrialGridFunctionSpace;
38 typedef GFSV TestGridFunctionSpace;
39
40 typedef typename GFSV::Ordering::Traits::ContainerIndex RowIndex;
41 typedef typename GFSU::Ordering::Traits::ContainerIndex ColIndex;
42
43 template<typename RowCache, typename ColCache>
44 using LocalView = UncachedMatrixView<MatrixContainer,RowCache,ColCache>;
45
46 template<typename RowCache, typename ColCache>
47 using ConstLocalView = ConstUncachedMatrixView<const MatrixContainer,RowCache,ColCache>;
48
49
50 template<typename GO>
51 MatrixContainer(const GO& go)
52 : _rows(go.testGridFunctionSpace().size())
53 , _cols(go.trialGridFunctionSpace().size())
54 , _container(std::make_shared<Container>(_rows*_cols,E(0)))
55 {}
56
57 template<typename GO>
58 MatrixContainer(const GO& go, const E& e)
59 : _rows(go.testGridFunctionSpace().size())
60 , _cols(go.trialGridFunctionSpace().size())
61 , _container(std::make_shared<Container>(_rows*_cols,e))
62 {}
63
65 explicit MatrixContainer(Backend::unattached_container = Backend::unattached_container())
66 : _rows(0)
67 , _cols(0)
68 {}
69
71 explicit MatrixContainer(Backend::attached_container)
72 : _rows(0)
73 , _cols(0)
74 , _container(std::make_shared<Container>())
75 {}
76
77 MatrixContainer(const MatrixContainer& rhs)
78 : _rows(rhs._rows)
79 , _cols(rhs._cols)
80 , _container(std::make_shared<Container>(*(rhs._container)))
81 {}
82
83 MatrixContainer& operator=(const MatrixContainer& rhs)
84 {
85 if (this == &rhs)
86 return *this;
87 if (_rows == 0 && _cols == 0)
88 {
89 _rows = rhs._rows;
90 _cols = rhs._cols;
91 }
92 if (attached())
93 {
94 (*_container) = (*(rhs._container));
95 }
96 else
97 {
98 _container = std::make_shared<Container>(*(rhs._container));
99 }
100 return *this;
101 }
102
103 void detach()
104 {
105 _container.reset();
106 }
107
108 void attach(std::shared_ptr<Container> container)
109 {
110 _container = container;
111 }
112
113 bool attached() const
114 {
115 return bool(_container);
116 }
117
118 const std::shared_ptr<Container>& storage() const
119 {
120 return _container;
121 }
122
123 size_type N() const
124 {
125 return _rows;
126 }
127
128 size_type M() const
129 {
130 return _cols;
131 }
132
133 MatrixContainer& operator=(const E& e)
134 {
135 std::fill(_container->begin(),_container->end(),e);
136 return *this;
137 }
138
139 MatrixContainer& operator*=(const E& e)
140 {
141 using namespace std::placeholders;
142 std::transform(_container->begin(),_container->end(),_container->begin(),std::bind(std::multiplies<E>(),e,_1));
143 return *this;
144 }
145
146 template<typename V>
147 void mv(const V& x, V& y) const
148 {
149 auto rowit = _container->begin();
150 for (auto& v : y)
151 {
152 v = std::inner_product(rowit,rowit + _cols,x.begin(),E(0));
153 rowit += _cols;
154 }
155 }
156
157 template<typename V>
158 void usmv(const E alpha, const V& x, V& y) const
159 {
160 auto rowit = _container->begin();
161 for (auto& v : y)
162 {
163 v += alpha * std::inner_product(rowit,rowit + _cols,x.begin(),E(0));
164 rowit += _cols;
165 }
166 }
167
168 E& operator()(const RowIndex& ri, const ColIndex& ci)
169 {
170 return (*_container)[ri[0]*_cols + ci[0]];
171 }
172
173 const E& operator()(const RowIndex& ri, const ColIndex& ci) const
174 {
175 return (*_container)[ri[0]*_cols + ci[0]];
176 }
177
178 const Container& base() const
179 {
180 return *_container;
181 }
182
183 Container& base()
184 {
185 return *_container;
186 }
187
188 private:
189
190 const Container& native() const
191 {
192 return *_container;
193 }
194
195 Container& native()
196 {
197 return *_container;
198 }
199
200 public:
201
202 void flush()
203 {}
204
205 void finalize()
206 {}
207
208 void clear_row(const RowIndex& ri, const E& diagonal_entry)
209 {
210 std::fill(_container->begin() + ri[0]*_cols,_container->begin() + (ri[0]+1)*_cols,E(0));
211 (*this)(ri,ri) = diagonal_entry;
212 }
213
214 void clear_row_block(const RowIndex& ri, const E& diagonal_entry)
215 {
216 std::fill(_container->begin() + ri[0]*_cols,_container->begin() + (ri[0]+1)*_cols,E(0));
217 (*this)(ri,ri) = diagonal_entry;
218 }
219
220 private:
221
222 std::size_t _rows;
223 std::size_t _cols;
224 std::shared_ptr<Container> _container;
225 };
226
227 } // namespace Simple
228
229 } // namespace PDELab
230} // namespace Dune
231
232#endif // DUNE_PDELAB_BACKEND_SIMPLE_MATRIX_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
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
STL namespace.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 24, 23:30, 2024)