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