dune-pdelab  2.0.0
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 
10 
11 namespace Dune {
12  namespace PDELab {
13  namespace istl {
14 
15 #ifndef DOXYGEN
16 
17  // implementation namespace for matrix diagonal
18 
19  namespace diagonal {
20 
21  // TMP for determining the type of vector to use for the matrix diagonal
22  template<typename M>
23  struct matrix_element_vector;
24 
25  // At FieldMatrix level, we keep the whole matrix
26  template<typename E, int n, int m>
27  struct matrix_element_vector<
28  FieldMatrix<E,n,m>
29  >
30  {
31  typedef FieldMatrix<E,n,m> type;
32  };
33 
34  // At BCRSMatrix level, we use a BlockVector and recursively apply the
35  // TMP to the block type.
36  template<typename Block, typename Allocator>
37  struct matrix_element_vector<
38  BCRSMatrix<Block,Allocator>
39  >
40  {
41  typedef BlockVector<
42  typename matrix_element_vector<Block>::type,
43  Allocator
44  > type;
45  };
46 
47 
48  // Function for extracting the diagonal from a matrix.
49  // For the FieldMatrix, we just copy the complete matrix.
50  template<typename FieldMatrix>
51  void matrix_element_vector_from_matrix(tags::field_matrix, FieldMatrix& c, const FieldMatrix& matrix)
52  {
53  c = matrix;
54  }
55 
56  // For the BCRSMatrix, we recursively copy diagonal blocks.
57  template<typename BlockVector, typename BCRSMatrix>
58  void matrix_element_vector_from_matrix(tags::block_vector, BlockVector& c, const BCRSMatrix& m)
59  {
60  const std::size_t rows = m.N();
61  c.resize(rows,false);
62  for (std::size_t i = 0; i < rows; ++i)
63  matrix_element_vector_from_matrix(container_tag(c[i]),c[i],m[i][i]);
64  }
65 
66 
67  // Function for inverting the diagonal.
68  // The FieldMatrix supports direct inverson.
69  template<typename FieldMatrix>
70  void invert_blocks(tags::field_matrix, FieldMatrix& c)
71  {
72  c.invert();
73  }
74 
75  // For a BCRSMatrix, we recursively invert all diagonal entries.
76  template<typename BlockVector>
77  void invert_blocks(tags::block_vector, BlockVector& c)
78  {
79  const std::size_t rows = c.size();
80  for (std::size_t i = 0; i < rows; ++i)
81  invert_blocks(container_tag(c[i]),c[i]);
82  }
83 
84 
85  // Matrix-vector product between matrix consisting of only the
86  // diagonal and a matching vector.
87  // For the FieldMatrix, simply call its matrix-vector product.
88  template<typename FieldMatrix, typename X, typename Y>
89  void mv(tags::field_matrix, const FieldMatrix& c, const X& x, Y& y)
90  {
91  c.mv(x,y);
92  }
93 
94  // For the BCRSMatrix, recursively apply this function to the
95  // individual blocks.
96  template<typename BlockVector, typename X, typename Y>
97  void mv(tags::block_vector, const BlockVector& c, const X& x, Y& y)
98  {
99  const std::size_t rows = c.size();
100  for (std::size_t i = 0; i < rows; ++i)
101  mv(container_tag(c[i]),c[i],x[i],y[i]);
102  }
103 
104 
105  template<typename FieldMatrix, typename CI>
106  std::size_t row_size(tags::field_matrix, const FieldMatrix& c, const CI& ci, int i)
107  {
108  return FieldMatrix::cols;
109  }
110 
111  template<typename FieldMatrix>
112  std::size_t row_size(tags::field_matrix, const FieldMatrix& c)
113  {
114  return FieldMatrix::cols;
115  }
116 
117  template<typename BlockVector, typename CI>
118  std::size_t row_size(tags::block_vector, const BlockVector& c, const CI& ci, int i)
119  {
120  return row_size(container_tag(c[0]),c[0]);
121  }
122 
123  template<typename BlockVector>
124  std::size_t row_size(tags::block_vector, const BlockVector& c)
125  {
126  return row_size(container_tag(c[0]),c[0]);
127  }
128 
129  template<typename FieldMatrix, typename CI>
130  typename FieldMatrix::field_type* row_begin(tags::field_matrix_1_any, FieldMatrix& c, const CI& ci, int i)
131  {
132  assert(i == -1);
133  return &(*c[0].begin());
134  }
135 
136  template<typename FieldMatrix, typename CI>
137  typename FieldMatrix::field_type* row_begin(tags::field_matrix_n_any, FieldMatrix& c, const CI& ci, int i)
138  {
139  assert(i == 0);
140  return &(*c[ci[0]].begin());
141  }
142 
143  template<typename FieldMatrix, typename CI>
144  typename FieldMatrix::field_type* row_end(tags::field_matrix_1_any, FieldMatrix& c, const CI& ci, int i)
145  {
146  assert(i == -1);
147  return &(*c[0].end());
148  }
149 
150  template<typename FieldMatrix, typename CI>
151  typename FieldMatrix::field_type* row_end(tags::field_matrix_n_any, FieldMatrix& c, const CI& ci, int i)
152  {
153  assert(i == 0);
154  return &(*c[ci[0]].end());
155  }
156 
157 
158  template<typename BlockVector, typename CI>
159  typename BlockVector::field_type* row_begin(tags::block_vector, BlockVector& c, const CI& ci, std::size_t i)
160  {
161  return row_begin(container_tag(c[ci[i]]),c[ci[i]],ci,i-1);
162  }
163 
164  template<typename BlockVector, typename CI>
165  typename BlockVector::field_type* row_end(tags::block_vector, BlockVector& c, const CI& ci, std::size_t i)
166  {
167  return row_end(container_tag(c[ci[i]]),c[ci[i]],ci,i-1);
168  }
169 
170 
171  } // namespace diagonal
172 
173 #endif // DOXYGEN
174 
175 
176  template<typename M>
178  {
179 
180  typedef typename raw_type<M>::type Matrix;
181 
183  {
184 
185  typedef typename diagonal::matrix_element_vector<Matrix>::type Container;
186  typedef typename Container::field_type field_type;
188 
190 
192  {
193  diagonal::matrix_element_vector_from_matrix(container_tag(_container),_container,raw(m));
194  }
195 
196  void invert()
197  {
198  diagonal::invert_blocks(container_tag(_container),_container);
199  }
200 
201  template<typename X, typename Y>
202  void mv(const X& x, Y& y) const
203  {
204  diagonal::mv(container_tag(_container),_container,raw(x),raw(y));
205  }
206 
207  template<typename ContainerIndex>
208  std::size_t row_size(const ContainerIndex& ci) const
209  {
210  return diagonal::row_size(container_tag(_container),_container,ci,ci.size()-1);
211  }
212 
213  template<typename ContainerIndex>
214  iterator row_begin(const ContainerIndex& ci)
215  {
216  return diagonal::row_begin(container_tag(_container),_container,ci,ci.size()-1);
217  }
218 
219  template<typename ContainerIndex>
220  iterator row_end(const ContainerIndex& ci)
221  {
222  return diagonal::row_end(container_tag(_container),_container,ci,ci.size()-1);
223  }
224 
225  };
226 
227 
228  template<typename GFS>
230  : public Dune::CommDataHandleIF<AddMatrixElementVectorDataHandle<GFS>,typename Matrix::field_type>
231  {
232 
233  public:
234 
235  typedef typename Matrix::field_type DataType;
236  typedef typename GFS::Traits::SizeType size_type;
237 
239  : _gfs(gfs)
240  , _index_cache(gfs)
241  , _v(v)
242  {}
243 
245  bool contains(int dim, int codim) const
246  {
247  return _gfs.dataHandleContains(codim);
248  }
249 
251  bool fixedsize(int dim, int codim) const
252  {
253  return _gfs.dataHandleFixedSize(codim);
254  }
255 
260  template<typename Entity>
261  size_type size(Entity& e) const
262  {
263  _index_cache.update(e);
264 
265  size_type s = 0;
266  for (size_type i = 0; i < _index_cache.size(); ++i)
267  s += _v.row_size(_index_cache.containerIndex(i));
268  return s;
269  }
270 
272  template<typename MessageBuffer, typename Entity>
273  void gather(MessageBuffer& buff, const Entity& e) const
274  {
275  _index_cache.update(e);
276  for (size_type i = 0; i < _index_cache.size(); ++i)
277  {
278  const CI& ci = _index_cache.containerIndex(i);
279  for (RowIterator it = _v.row_begin(ci),
280  end_it = _v.row_end(ci);
281  it != end_it;
282  ++it)
283  buff.write(*it);
284  }
285  }
286 
291  template<typename MessageBuffer, typename Entity>
292  void scatter(MessageBuffer& buff, const Entity& e, size_type n)
293  {
294  _index_cache.update(e);
295  for (size_type i = 0; i < _index_cache.size(); ++i)
296  {
297  const CI& ci = _index_cache.containerIndex(i);
298  for (RowIterator it = _v.row_begin(ci),
299  end_it = _v.row_end(ci);
300  it != end_it;
301  ++it)
302  {
303  DataType x;
304  buff.read(x);
305  *it += x;
306  }
307  }
308  }
309 
310  private:
311 
312  typedef EntityIndexCache<GFS> IndexCache;
313  typedef typename IndexCache::ContainerIndex CI;
314  typedef typename MatrixElementVector::iterator RowIterator;
315 
316  const GFS& _gfs;
317  mutable IndexCache _index_cache;
319 
320  };
321 
322  };
323 
324  } // namespace istl
325  } // namespace PDELab
326 } // namespace Dune
327 
328 #endif // DUNE_PDELAB_BACKEND_ISTL_BLOCKMATRIXDIAGONAL_HH
void update(const Entity &e)
Definition: entityindexcache.hh:49
bool fixedsize(int dim, int codim) const
returns true if size per entity of given dim and codim is a constant
Definition: blockmatrixdiagonal.hh:251
Container _container
Definition: blockmatrixdiagonal.hh:189
void mv(const X &x, Y &y) const
Definition: blockmatrixdiagonal.hh:202
AddMatrixElementVectorDataHandle(const GFS &gfs, MatrixElementVector &v)
Definition: blockmatrixdiagonal.hh:238
const CI & containerIndex(size_type i) const
Definition: entityindexcache.hh:64
size_type size(Entity &e) const
how many objects of type DataType have to be sent for a given entity
Definition: blockmatrixdiagonal.hh:261
void scatter(MessageBuffer &buff, const Entity &e, size_type n)
unpack data from message buffer to user
Definition: blockmatrixdiagonal.hh:292
static const int dim
Definition: adaptivity.hh:82
diagonal::matrix_element_vector< Matrix >::type Container
Definition: blockmatrixdiagonal.hh:185
std::size_t row_size(const ContainerIndex &ci) const
Definition: blockmatrixdiagonal.hh:208
C type
Definition: backend/istl/utility.hh:42
size_type size() const
Definition: entityindexcache.hh:74
V & raw(V &v)
Returns the raw ISTL object associated with v, or v itself it is already an ISTL object.
Definition: backend/istl/utility.hh:26
void invert()
Definition: blockmatrixdiagonal.hh:196
field_type * iterator
Definition: blockmatrixdiagonal.hh:187
iterator row_begin(const ContainerIndex &ci)
Definition: blockmatrixdiagonal.hh:214
raw_type< M >::type Matrix
Definition: blockmatrixdiagonal.hh:180
MatrixElementVector(const M &m)
Definition: blockmatrixdiagonal.hh:191
GFS::Traits::SizeType size_type
Definition: blockmatrixdiagonal.hh:236
void gather(MessageBuffer &buff, const Entity &e) const
pack data from user to message buffer
Definition: blockmatrixdiagonal.hh:273
bool contains(int dim, int codim) const
returns true if data for this codim should be communicated
Definition: blockmatrixdiagonal.hh:245
tags::container< T >::type container_tag(const T &)
Gets instance of container tag associated with T.
Definition: backend/istl/tags.hh:247
Container::field_type field_type
Definition: blockmatrixdiagonal.hh:186
Matrix::field_type DataType
Definition: blockmatrixdiagonal.hh:235
Definition: blockmatrixdiagonal.hh:177
const E & e
Definition: interpolate.hh:172
iterator row_end(const ContainerIndex &ci)
Definition: blockmatrixdiagonal.hh:220
const std::string s
Definition: function.hh:1103
Ordering::Traits::ContainerIndex ContainerIndex
Definition: entityindexcache.hh:23