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