1 #ifndef DUNE_PDELAB_EIGENMATRIXBACKEND_HH
2 #define DUNE_PDELAB_EIGENMATRIXBACKEND_HH
11 #include <Eigen/Sparse>
21 struct MatrixPatternInserter
23 MatrixPatternInserter(M & mat)
27 template<
typename RI,
typename CI>
28 void add_link(
const RI& ri,
const CI& ci)
30 _matrix.coeffRef(ri.back(),ci.back()) = 0.0;
36 template<
typename GFSV,
typename GFSU,
typename ET,
int _Options>
42 typedef Eigen::SparseMatrix<ET,_Options> Container;
43 typedef ET ElementType;
45 typedef ElementType field_type;
46 typedef typename Container::Index size_type;
47 typedef typename Container::Index index_type;
49 typedef GFSU TrialGridFunctionSpace;
50 typedef GFSV TestGridFunctionSpace;
52 typedef typename GFSV::Ordering::Traits::ContainerIndex RowIndex;
53 typedef typename GFSU::Ordering::Traits::ContainerIndex ColIndex;
55 template<
typename RowCache,
typename ColCache>
56 using LocalView = UncachedMatrixView<MatrixContainer,RowCache,ColCache>;
58 template<
typename RowCache,
typename ColCache>
59 using ConstLocalView = ConstUncachedMatrixView<const MatrixContainer,RowCache,ColCache>;
62 typename GFSV::Ordering::Traits::DOFIndex,
63 typename GFSV::Ordering::Traits::ContainerIndex
67 typename GFSU::Ordering::Traits::DOFIndex,
68 typename GFSU::Ordering::Traits::ContainerIndex
71 typedef MatrixPatternInserter<Container> Pattern;
74 MatrixContainer(
const GO& go,
int avg_per_row)
75 : _container(make_shared<Container>())
77 allocate_matrix(_container, go, avg_per_row, ElementType(0));
81 MatrixContainer(
const GO& go,
int avg_per_row,
const ElementType&
e)
82 : _container(make_shared<Container>())
84 allocate_matrix(_container, go, avg_per_row, e);
88 explicit MatrixContainer(tags::unattached_container = tags::unattached_container())
92 explicit MatrixContainer(tags::attached_container)
93 : _container(make_shared<Container>())
96 MatrixContainer(
const MatrixContainer& rhs)
97 : _container(make_shared<Container>(*(rhs._container)))
100 MatrixContainer& operator=(
const MatrixContainer& rhs)
110 _container = make_shared<Container>(rhs.base());
120 void attach(shared_ptr<Container> container)
122 _container = container;
125 bool attached()
const
127 return bool(_container);
130 const shared_ptr<Container>& storage()
const
137 return _container->rows();
142 return _container->cols();
145 MatrixContainer& operator=(
const ElementType&
e)
147 if(!_container->isCompressed()) _container->makeCompressed();
148 for (
int i=0; i<_container->nonZeros(); i++)
149 _container->valuePtr()[i] =
e;
155 MatrixContainer& operator*=(
const ElementType& e)
162 void mv(
const V& x, V& y)
const
164 y.base() = base() * x.base();
168 void usmv(
const ElementType alpha,
const V& x, V& y)
const
170 y.base() += alpha * (base() * x.base());
173 ElementType& operator()(
const RowIndex& ri,
const ColIndex& ci)
175 return _container->coeffRef(ri[0],ci[0]);
178 const ElementType operator()(
const RowIndex& ri,
const ColIndex& ci)
const
180 return _container->coeffRef(ri[0],ci[0]);
183 const Container& base()
const
199 void clear_row(
const RowIndex& ri,
const ElementType& diagonal_entry)
201 _container->middleRows(ri[0],1) *= 0.0;
202 _container->coeffRef(ri[0],ri[0]) = diagonal_entry;
206 template<
typename GO>
207 static void allocate_matrix(shared_ptr<Container> & c,
const GO & go,
int avg_per_row,
const ElementType& e)
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));
216 go.fill_pattern(pattern);
221 shared_ptr< Container > _container;
const E & e
Definition: interpolate.hh:172