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