dune-pdelab  2.7-git
eigen/matrix.hh
Go to the documentation of this file.
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 
15 #include <Eigen/Sparse>
16 
17 namespace 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
const Entity & e
Definition: localfunctionspace.hh:121
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
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.