dune-pdelab  2.5-dev
blockmatrixdiagonal.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #ifndef DUNE_PDELAB_BACKEND_ISTL_BLOCKMATRIXDIAGONAL_HH
4 #define DUNE_PDELAB_BACKEND_ISTL_BLOCKMATRIXDIAGONAL_HH
5 
6 // this is here for backwards compatibility and deprecation warnings, remove after 2.5.0
7 #include "ensureistlinclude.hh"
8 
13 
14 namespace Dune {
15  namespace PDELab {
16  namespace ISTL {
17 
18 #ifndef DOXYGEN
19 
20  // implementation namespace for matrix diagonal
21 
22  namespace diagonal {
23 
24  // TMP for determining the type of vector to use for the matrix diagonal
25  template<typename M>
26  struct matrix_element_vector;
27 
28  // At FieldMatrix level, we keep the whole matrix
29  template<typename E, int n, int m>
30  struct matrix_element_vector<
31  FieldMatrix<E,n,m>
32  >
33  {
34  typedef FieldMatrix<E,n,m> type;
35  };
36 
37  // At BCRSMatrix level, we use a BlockVector and recursively apply the
38  // TMP to the block type.
39  template<typename Block, typename Allocator>
40  struct matrix_element_vector<
41  Dune::BCRSMatrix<Block,Allocator>
42  >
43  {
44  typedef Dune::BlockVector<
45  typename matrix_element_vector<Block>::type,
46  Allocator
47  > type;
48  };
49 
50 
51  // Function for extracting the diagonal from a matrix.
52  // For the FieldMatrix, we just copy the complete matrix.
53  template<typename FieldMatrix>
54  void matrix_element_vector_from_matrix(tags::field_matrix, FieldMatrix& c, const FieldMatrix& matrix)
55  {
56  c = matrix;
57  }
58 
59  // For the BCRSMatrix, we recursively copy diagonal blocks.
60  template<typename BlockVector, typename BCRSMatrix>
61  void matrix_element_vector_from_matrix(tags::block_vector, BlockVector& c, const BCRSMatrix& m)
62  {
63  const std::size_t rows = m.N();
64  c.resize(rows,false);
65  for (std::size_t i = 0; i < rows; ++i)
66  matrix_element_vector_from_matrix(container_tag(c[i]),c[i],m[i][i]);
67  }
68 
69 
70  // Function for inverting the diagonal.
71  // The FieldMatrix supports direct inverson.
72  template<typename FieldMatrix>
73  void invert_blocks(tags::field_matrix, FieldMatrix& c)
74  {
75  c.invert();
76  }
77 
78  // For a BCRSMatrix, we recursively invert all diagonal entries.
79  template<typename BlockVector>
80  void invert_blocks(tags::block_vector, BlockVector& c)
81  {
82  const std::size_t rows = c.size();
83  for (std::size_t i = 0; i < rows; ++i)
84  invert_blocks(container_tag(c[i]),c[i]);
85  }
86 
87 
88  // Matrix-vector product between matrix consisting of only the
89  // diagonal and a matching vector.
90  // For the FieldMatrix, simply call its matrix-vector product.
91  template<typename FieldMatrix, typename X, typename Y>
92  void mv(tags::field_matrix, const FieldMatrix& c, const X& x, Y& y)
93  {
94  c.mv(x,y);
95  }
96 
97  // For the BCRSMatrix, recursively apply this function to the
98  // individual blocks.
99  template<typename BlockVector, typename X, typename Y>
100  void mv(tags::block_vector, const BlockVector& c, const X& x, Y& y)
101  {
102  const std::size_t rows = c.size();
103  for (std::size_t i = 0; i < rows; ++i)
104  mv(container_tag(c[i]),c[i],x[i],y[i]);
105  }
106 
107 
108  // We don't know the type of the container that stores the actual field values (double, complex,...)
109  // Moreover, there might be different containers in case of heterogeneous matrices (multidomain etc.)
110  // In order to communicate the matrix blocks, we need to stream the single row inside the lowest-level
111  // diagonal matrix block for each DOF.
112  // The following set of functions extracts this information from the container and provides a simple
113  // interface that uses pointers to the field type as iterators.
114  //
115  // WARNING: This assumes that matrix blocks at the lowest level are dense and stored in column-major format!
116  //
117 
118  template<typename FieldMatrix, typename CI>
119  std::size_t row_size(tags::field_matrix, const FieldMatrix& c, const CI& ci, int i)
120  {
121  return FieldMatrix::cols;
122  }
123 
124  template<typename FieldMatrix>
125  std::size_t row_size(tags::field_matrix, const FieldMatrix& c)
126  {
127  return FieldMatrix::cols;
128  }
129 
130  template<typename BlockVector, typename CI>
131  std::size_t row_size(tags::block_vector, const BlockVector& c, const CI& ci, int i)
132  {
133  return row_size(container_tag(c[0]),c[0]);
134  }
135 
136  template<typename BlockVector>
137  std::size_t row_size(tags::block_vector, const BlockVector& c)
138  {
139  return row_size(container_tag(c[0]),c[0]);
140  }
141 
142  // FieldMatrix with a single row is special because the last-level index isn't stored, so we have to
143  // manually extract row 0.
144  template<typename FieldMatrix, typename CI>
145  typename FieldMatrix::field_type* row_begin(tags::field_matrix_1_any, FieldMatrix& c, const CI& ci, int i)
146  {
147  assert(i == -1);
148  return &(*c[0].begin());
149  }
150 
151  template<typename FieldMatrix, typename CI>
152  typename FieldMatrix::field_type* row_begin(tags::field_matrix_n_any, FieldMatrix& c, const CI& ci, int i)
153  {
154  assert(i == 0);
155  return &(*c[ci[0]].begin());
156  }
157 
158 
159  // The end iterators are a little tricky: We want a pointer to the memory location directly after the last
160  // entry for the given row. In theory, we could get this location by dereferencing the end() iterator and
161  // then taking the address of that location, but we are not allowed to dereference an end() iterator. So
162  // we instead decrement the end() iterator by one, take the (valid) address of the element at that location
163  // and increment that pointer by 1. Yay for ugly hackery! :-P
164 
165  // With a 1x1 matrix, we can simply take the address directly following the begin() iterator's target.
166  template<typename FieldMatrix, typename CI>
167  typename FieldMatrix::field_type* row_end(tags::field_matrix_1_1, FieldMatrix& c, const CI& ci, int i)
168  {
169  assert(i == -1);
170  return &(*c[0].begin()) + 1;
171  }
172 
173  // For any other matrix, we perform the decrement iterator / increment address of target dance...
174  // Once for the optimized storage scheme of single row matrices...
175  template<typename FieldMatrix, typename CI>
176  typename FieldMatrix::field_type* row_end(tags::field_matrix_1_any, FieldMatrix& c, const CI& ci, int i)
177  {
178  assert(i == -1);
179  typename FieldMatrix::row_type::iterator it = c[0].end();
180  --it;
181  return &(*it) + 1;
182  }
183 
184  // ... and once for the general case.
185  template<typename FieldMatrix, typename CI>
186  typename FieldMatrix::field_type* row_end(tags::field_matrix_n_any, FieldMatrix& c, const CI& ci, int i)
187  {
188  assert(i == 0);
189  typename FieldMatrix::row_type::iterator it = c[ci[0]].end();
190  --it;
191  return &(*it) + 1;
192  }
193 
194 
195  // These are the standard begin() and end() methods for BlockVector. They recursvely call row_begin()
196  // to arrive at the lowest level block structure.
197 
198  template<typename BlockVector, typename CI>
199  typename BlockVector::field_type* row_begin(tags::block_vector, BlockVector& c, const CI& ci, std::size_t i)
200  {
201  return row_begin(container_tag(c[ci[i]]),c[ci[i]],ci,i-1);
202  }
203 
204  template<typename BlockVector, typename CI>
205  typename BlockVector::field_type* row_end(tags::block_vector, BlockVector& c, const CI& ci, std::size_t i)
206  {
207  return row_end(container_tag(c[ci[i]]),c[ci[i]],ci,i-1);
208  }
209 
210 
211  } // namespace diagonal
212 
213 #endif // DOXYGEN
214 
215 
216  template<typename M>
218  {
219 
221 
223  {
224 
225  typedef typename diagonal::matrix_element_vector<Matrix>::type Container;
226  typedef typename Container::field_type field_type;
227  typedef field_type* iterator;
228 
229  Container _container;
230 
232  {
233  diagonal::matrix_element_vector_from_matrix(container_tag(_container),_container,Backend::native(m));
234  }
235 
236  void invert()
237  {
238  diagonal::invert_blocks(container_tag(_container),_container);
239  }
240 
241  template<typename X, typename Y>
242  void mv(const X& x, Y& y) const
243  {
244  diagonal::mv(container_tag(_container),_container,Backend::native(x),Backend::native(y));
245  }
246 
247  template<typename ContainerIndex>
248  std::size_t row_size(const ContainerIndex& ci) const
249  {
250  return diagonal::row_size(container_tag(_container),_container,ci,ci.size()-1);
251  }
252 
253  template<typename ContainerIndex>
254  iterator row_begin(const ContainerIndex& ci)
255  {
256  return diagonal::row_begin(container_tag(_container),_container,ci,ci.size()-1);
257  }
258 
259  template<typename ContainerIndex>
260  iterator row_end(const ContainerIndex& ci)
261  {
262  return diagonal::row_end(container_tag(_container),_container,ci,ci.size()-1);
263  }
264 
265  };
266 
267 
268  template<typename GFS>
270  : public Dune::CommDataHandleIF<AddMatrixElementVectorDataHandle<GFS>,typename Matrix::field_type>
271  {
272 
273  public:
274 
275  typedef typename Matrix::field_type DataType;
276  typedef typename GFS::Traits::SizeType size_type;
277 
279  : _gfs(gfs)
280  , _index_cache(gfs)
281  , _v(v)
282  {}
283 
285  bool contains(int dim, int codim) const
286  {
287  return _gfs.dataHandleContains(codim);
288  }
289 
291  bool fixedsize(int dim, int codim) const
292  {
293  return _gfs.dataHandleFixedSize(codim);
294  }
295 
300  template<typename Entity>
301  size_type size(Entity& e) const
302  {
303  _index_cache.update(e);
304 
305  size_type s = 0;
306  for (size_type i = 0; i < _index_cache.size(); ++i)
307  s += _v.row_size(_index_cache.containerIndex(i));
308  return s;
309  }
310 
312  template<typename MessageBuffer, typename Entity>
313  void gather(MessageBuffer& buff, const Entity& e) const
314  {
315  _index_cache.update(e);
316  for (size_type i = 0; i < _index_cache.size(); ++i)
317  {
318  const CI& ci = _index_cache.containerIndex(i);
319  for (RowIterator it = _v.row_begin(ci),
320  end_it = _v.row_end(ci);
321  it != end_it;
322  ++it)
323  buff.write(*it);
324  }
325  }
326 
331  template<typename MessageBuffer, typename Entity>
332  void scatter(MessageBuffer& buff, const Entity& e, size_type n)
333  {
334  _index_cache.update(e);
335  for (size_type i = 0; i < _index_cache.size(); ++i)
336  {
337  const CI& ci = _index_cache.containerIndex(i);
338  for (RowIterator it = _v.row_begin(ci),
339  end_it = _v.row_end(ci);
340  it != end_it;
341  ++it)
342  {
343  DataType x;
344  buff.read(x);
345  *it += x;
346  }
347  }
348  }
349 
350  private:
351 
353  typedef typename IndexCache::ContainerIndex CI;
354  typedef typename MatrixElementVector::iterator RowIterator;
355 
356  const GFS& _gfs;
357  mutable IndexCache _index_cache;
359 
360  };
361 
362  };
363 
364  } // namespace ISTL
365  } // namespace PDELab
366 } // namespace Dune
367 
368 #endif // DUNE_PDELAB_BACKEND_ISTL_BLOCKMATRIXDIAGONAL_HH
static const int dim
Definition: adaptivity.hh:83
Definition: blockmatrixdiagonal.hh:217
const std::string s
Definition: function.hh:1101
bool fixedsize(int dim, int codim) const
returns true if size per entity of given dim and codim is a constant
Definition: blockmatrixdiagonal.hh:291
iterator row_begin(const ContainerIndex &ci)
Definition: blockmatrixdiagonal.hh:254
void invert()
Definition: blockmatrixdiagonal.hh:236
Container _container
Definition: blockmatrixdiagonal.hh:229
size_type size(Entity &e) const
how many objects of type DataType have to be sent for a given entity
Definition: blockmatrixdiagonal.hh:301
Matrix::field_type DataType
Definition: blockmatrixdiagonal.hh:275
void gather(MessageBuffer &buff, const Entity &e) const
pack data from user to message buffer
Definition: blockmatrixdiagonal.hh:313
MatrixElementVector(const M &m)
Definition: blockmatrixdiagonal.hh:231
void mv(const X &x, Y &y) const
Definition: blockmatrixdiagonal.hh:242
iterator row_end(const ContainerIndex &ci)
Definition: blockmatrixdiagonal.hh:260
Backend::Native< M > Matrix
Definition: blockmatrixdiagonal.hh:220
std::size_t row_size(const ContainerIndex &ci) const
Definition: blockmatrixdiagonal.hh:248
For backward compatibility – Do not use this!
Definition: adaptivity.hh:27
const Entity & e
Definition: localfunctionspace.hh:111
Container::field_type field_type
Definition: istl/vector.hh:39
Ordering::Traits::ContainerIndex ContainerIndex
Definition: entityindexcache.hh:23
std::enable_if< std::is_base_of< impl::WrapperBase, T >::value, Native< T > &>::type native(T &t)
Definition: backend/interface.hh:192
bool contains(int dim, int codim) const
returns true if data for this codim should be communicated
Definition: blockmatrixdiagonal.hh:285
Container::field_type field_type
Definition: blockmatrixdiagonal.hh:226
typename native_type< T >::type Native
Alias of the native container type associated with T or T itself if it is not a backend wrapper...
Definition: backend/interface.hh:176
GFS::Traits::SizeType size_type
Definition: blockmatrixdiagonal.hh:276
AddMatrixElementVectorDataHandle(const GFS &gfs, MatrixElementVector &v)
Definition: blockmatrixdiagonal.hh:278
diagonal::matrix_element_vector< Matrix >::type Container
Definition: blockmatrixdiagonal.hh:225
tags::container< T >::type container_tag(const T &)
Gets instance of container tag associated with T.
Definition: backend/istl/tags.hh:249
void scatter(MessageBuffer &buff, const Entity &e, size_type n)
unpack data from message buffer to user
Definition: blockmatrixdiagonal.hh:332
field_type * iterator
Definition: blockmatrixdiagonal.hh:227