DUNE PDELab (git)

matrix.hh
1// -*- tab-width: 4; indent-tabs-mode: nil -*-
2// vi: set et ts=4 sw=2 sts=2:
3#ifndef DUNE_PDELAB_BACKEND_EIGEN_MATRIX_HH
4#define DUNE_PDELAB_BACKEND_EIGEN_MATRIX_HH
5
6#include <utility>
7#include <vector>
8#include <set>
9
10#if HAVE_EIGEN
11
12#include <dune/pdelab/backend/common/uncachedmatrixview.hh>
14#include <dune/pdelab/ordering/orderingbase.hh>
15#include <Eigen/Sparse>
16
17namespace Dune
18{
19 namespace PDELab
20 {
21 namespace Eigen
22 {
23
24 template<typename M>
25 struct MatrixPatternInserter
26 {
27 MatrixPatternInserter(M & mat)
28 : _matrix(mat)
29 {}
30
31 template<typename RI, typename CI>
32 void add_link(const RI& ri, const CI& ci)
33 {
34 _matrix.coeffRef(ri.back(),ci.back()) = 0.0;
35 }
36
37 M & _matrix;
38 };
39
46 template<typename GFSV, typename GFSU, typename ET, int _Options>
47 class MatrixContainer
48 : public Backend::impl::Wrapper<::Eigen::SparseMatrix<ET,_Options>>
49 {
50
51 public:
52
53 typedef ::Eigen::SparseMatrix<ET,_Options> Container;
54
55 private:
56
57 friend Backend::impl::Wrapper<Container>;
58
59 public:
60
61 typedef ET ElementType;
62
63 typedef ElementType field_type;
64 typedef typename Container::Index size_type;
65 typedef typename Container::Index index_type;
66
67 typedef GFSU TrialGridFunctionSpace;
68 typedef GFSV TestGridFunctionSpace;
69
70 typedef typename GFSV::Ordering::Traits::ContainerIndex RowIndex;
71 typedef typename GFSU::Ordering::Traits::ContainerIndex ColIndex;
72
73 template<typename RowCache, typename ColCache>
74 using LocalView = UncachedMatrixView<MatrixContainer,RowCache,ColCache>;
75
76 template<typename RowCache, typename ColCache>
77 using ConstLocalView = ConstUncachedMatrixView<const MatrixContainer,RowCache,ColCache>;
78
79 typedef OrderingBase<
80 typename GFSV::Ordering::Traits::DOFIndex,
81 typename GFSV::Ordering::Traits::ContainerIndex
82 > RowOrdering;
83
84 typedef OrderingBase<
85 typename GFSU::Ordering::Traits::DOFIndex,
86 typename GFSU::Ordering::Traits::ContainerIndex
87 > ColOrdering;
88
89 typedef MatrixPatternInserter<Container> Pattern;
90
91 template<typename GO>
92 MatrixContainer(const GO& go)
93 : _container(std::make_shared<Container>())
94 {
95 allocate_matrix(_container, go, ElementType(0));
96 }
97
98 template<typename GO>
99 MatrixContainer(const GO& go, const ElementType& e)
100 : _container(std::make_shared<Container>())
101 {
102 allocate_matrix(_container, go, e);
103 }
104
106 explicit MatrixContainer(Backend::unattached_container = Backend::unattached_container())
107 {}
108
110 explicit MatrixContainer(Backend::attached_container)
111 : _container(std::make_shared<Container>())
112 {}
113
114 MatrixContainer(const MatrixContainer& rhs)
115 : _container(std::make_shared<Container>(*(rhs._container)))
116 {}
117
118 MatrixContainer& operator=(const MatrixContainer& rhs)
119 {
120 if (this == &rhs)
121 return *this;
122 if (attached())
123 {
124 base() = rhs.base();
125 }
126 else
127 {
128 _container = std::make_shared<Container>(rhs.base());
129 }
130 return *this;
131 }
132
133 void detach()
134 {
135 _container.reset();
136 }
137
138 void attach(std::shared_ptr<Container> container)
139 {
140 _container = container;
141 }
142
143 bool attached() const
144 {
145 return bool(_container);
146 }
147
148 const std::shared_ptr<Container>& storage() const
149 {
150 return _container;
151 }
152
153 size_type N() const
154 {
155 return _container->rows();
156 }
157
158 size_type M() const
159 {
160 return _container->cols();
161 }
162
163 MatrixContainer& operator=(const ElementType& e)
164 {
165 if(!_container->isCompressed()) _container->makeCompressed();
166 for (int i=0; i<_container->nonZeros(); i++)
167 _container->valuePtr()[i] = e;
168 // we require a sufficiently new Eigen version to use setConstant (newer than in debian testing)
169 // _container->setConstant(e);
170 return *this;
171 }
172
173 MatrixContainer& operator*=(const ElementType& e)
174 {
175 base() *= e;
176 return *this;
177 }
178
179 template<typename V>
180 void mv(const V& x, V& y) const
181 {
182 y.base() = base() * x.base();
183 }
184
185 template<typename V>
186 void usmv(const ElementType alpha, const V& x, V& y) const
187 {
188 y.base() += alpha * (base() * x.base());
189 }
190
191 ElementType& operator()(const RowIndex& ri, const ColIndex& ci)
192 {
193 return _container->coeffRef(ri[0],ci[0]);
194 }
195
196 const ElementType operator()(const RowIndex& ri, const ColIndex& ci) const
197 {
198 return _container->coeffRef(ri[0],ci[0]);
199 }
200
201 const Container& base() const
202 {
203 return *_container;
204 }
205
206 Container& base()
207 {
208 return *_container;
209 }
210
211 private:
212
213 const Container& native() const
214 {
215 return *_container;
216 }
217
218 Container& native()
219 {
220 return *_container;
221 }
222
223 public:
224
225 void flush()
226 {}
227
228 void finalize()
229 {}
230
231 void clear_row(const RowIndex& ri, const ElementType& diagonal_entry)
232 {
233 _container->middleRows(ri[0],1) *= 0.0;
234 _container->coeffRef(ri[0],ri[0]) = diagonal_entry;
235 }
236
237 protected:
238 template<typename GO>
239 static void allocate_matrix(std::shared_ptr<Container> & c, const GO & go, const ElementType& e)
240 {
241 // guess size
242 int rows = go.testGridFunctionSpace().ordering().blockCount();
243 int cols = go.trialGridFunctionSpace().ordering().blockCount();
244 c->resize(rows,cols);
245 size_type nz = go.matrixBackend().avg_nz_per_row;
246 if (nz)
247 c->reserve(::Eigen::VectorXi::Constant(rows,nz));
248 // setup pattern
249 Pattern pattern(*c);
250 go.fill_pattern(pattern);
251 // compress matrix
252 c->makeCompressed();
253 }
254
255 std::shared_ptr< Container > _container;
256 };
257
258 } // end namespace EIGEN
259 } // namespace PDELab
260} // namespace Dune
261
262#endif /* HAVE_EIGEN */
263
264#endif // DUNE_PDELAB_BACKEND_EIGEN_MATRIX_HH
Various tags for influencing backend behavior.
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 12, 23:30, 2024)