dune-pdelab  2.5-dev
matrixhelpers.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_MATRIXHELPERS_HH
4 #define DUNE_PDELAB_BACKEND_ISTL_MATRIXHELPERS_HH
5 
6 #include<utility>
7 #include<vector>
8 #include <unordered_map>
9 #include <unordered_set>
10 
11 #include <dune/common/fmatrix.hh>
12 #include <dune/istl/bcrsmatrix.hh>
13 
16 
17 namespace Dune {
18  namespace PDELab {
19 
20 #ifndef DOXYGEN
21 
22  namespace ISTL {
23 
24  template<typename RV, typename CV, typename block_type>
25  struct matrix_for_vectors;
26 
27  template<typename B1, typename A1, typename B2, typename A2, typename block_type>
28  struct matrix_for_vectors<Dune::BlockVector<B1,A1>,Dune::BlockVector<B2,A2>,block_type>
29  {
30  typedef Dune::BCRSMatrix<block_type> type;
31  };
32 
33  template<typename B1, int n1, typename B2, int n2, typename block_type>
34  struct matrix_for_vectors<Dune::FieldVector<B1,n1>,Dune::FieldVector<B2,n2>,block_type>
35  {
36  typedef Dune::FieldMatrix<block_type,n1,n2> type;
37  };
38 
39  template<typename E, typename RV, typename CV, std::size_t blocklevel>
40  struct recursive_build_matrix_type
41  {
42  typedef typename matrix_for_vectors<RV,CV,typename recursive_build_matrix_type<E,typename RV::block_type,typename CV::block_type,blocklevel-1>::type>::type type;
43  };
44 
45  template<typename E, typename RV, typename CV>
46  struct recursive_build_matrix_type<E,RV,CV,1>
47  {
48  typedef Dune::FieldMatrix<E,RV::dimension,CV::dimension> type;
49  };
50 
51 
52  template<typename E, typename RV, typename CV>
53  struct build_matrix_type
54  {
55 
56  static_assert(static_cast<int>(RV::blocklevel) == static_cast<int>(CV::blocklevel),"Both vectors must have identical blocking depth");
57 
58  typedef typename recursive_build_matrix_type<E,RV,CV,RV::blocklevel>::type type;
59 
60  };
61 
62  template<typename RowOrdering, typename ColOrdering, typename SubPattern_ = void>
63  class Pattern
64  : public std::vector<std::unordered_map<std::size_t,SubPattern_> >
65  {
66 
67  public:
68 
69  typedef SubPattern_ SubPattern;
70 
71  template<typename RI, typename CI>
72  void add_link(const RI& ri, const CI& ci)
73  {
74  recursive_add_entry(ri.view(),ci.view());
75  }
76 
77  template<typename RI, typename CI>
78  void recursive_add_entry(const RI& ri, const CI& ci)
79  {
80  this->resize(_row_ordering.blockCount());
81  std::pair<typename std::unordered_map<std::size_t,SubPattern>::iterator,bool> r = (*this)[ri.back()].insert(make_pair(ci.back(),SubPattern(_row_ordering.childOrdering(ri.back()),_col_ordering.childOrdering(ci.back()))));
82  r.first->second.recursive_add_entry(ri.back_popped(),ci.back_popped());
83  }
84 
85  Pattern(const RowOrdering& row_ordering, const ColOrdering& col_ordering)
86  : _row_ordering(row_ordering)
87  , _col_ordering(col_ordering)
88  {}
89 
90  private:
91 
92  const RowOrdering& _row_ordering;
93  const ColOrdering& _col_ordering;
94 
95  };
96 
97  template<typename RowOrdering, typename ColOrdering>
98  class Pattern<RowOrdering,ColOrdering,void>
99  : public std::vector<std::unordered_set<std::size_t> >
100  {
101 
102  public:
103 
104  typedef void SubPattern;
105 
106  template<typename RI, typename CI>
107  void add_link(const RI& ri, const CI& ci)
108  {
109  recursive_add_entry(ri,ci);
110  }
111 
112  template<typename RI, typename CI>
113  void recursive_add_entry(const RI& ri, const CI& ci)
114  {
115  this->resize(_row_ordering.blockCount());
116  (*this)[ri.back()].insert(ci.back());
117  }
118 
119  Pattern(const RowOrdering& row_ordering, const ColOrdering& col_ordering)
120  : _row_ordering(row_ordering)
121  , _col_ordering(col_ordering)
122  {}
123 
124  private:
125 
126  const RowOrdering& _row_ordering;
127  const ColOrdering& _col_ordering;
128 
129  };
130 
131  template<typename M, int blocklevel = M::blocklevel>
132  struct requires_pattern
133  {
134  static const bool value = requires_pattern<typename M::block_type,blocklevel-1>::value;
135  };
136 
137  template<typename M>
138  struct requires_pattern<M,0>
139  {
140  static const bool value = false;
141  };
142 
143  template<typename B, typename A, int blocklevel>
144  struct requires_pattern<Dune::BCRSMatrix<B,A>,blocklevel>
145  {
146  static const bool value = true;
147  };
148 
149  template<typename M, typename RowOrdering, typename ColOrdering, bool pattern>
150  struct _build_pattern_type
151  {
152  typedef void type;
153  };
154 
155  template<typename M, typename RowOrdering, typename ColOrdering>
156  struct _build_pattern_type<M,RowOrdering,ColOrdering,true>
157  {
159  };
160 
161  template<typename M, typename GFSV, typename GFSU, typename Tag>
162  struct build_pattern_type
163  {
164 
165  typedef OrderingBase<
166  typename GFSV::Ordering::Traits::DOFIndex,
167  typename GFSV::Ordering::Traits::ContainerIndex
168  > RowOrdering;
169 
170  typedef OrderingBase<
171  typename GFSU::Ordering::Traits::DOFIndex,
172  typename GFSU::Ordering::Traits::ContainerIndex
173  > ColOrdering;
174 
176  };
177 
178  template<typename M, typename GFSV, typename GFSU>
179  struct build_pattern_type<M,GFSV,GFSU,FlatContainerAllocationTag>
180  {
181  typedef Pattern<typename GFSV::Ordering, typename GFSU::Ordering> type;
182  };
183 
184 
185  template<typename RI, typename CI, typename Block>
186  typename Block::field_type&
187  access_matrix_element(tags::field_matrix_1_1, Block& b, const RI& ri, const CI& ci, int i, int j)
188  {
189  assert(i == -1);
190  assert(j == -1);
191  return b[0][0];
192  }
193 
194  template<typename RI, typename CI, typename Block>
195  typename Block::field_type&
196  access_matrix_element(tags::field_matrix_n_m, Block& b, const RI& ri, const CI& ci, int i, int j)
197  {
198  assert(i == 0);
199  assert(j == 0);
200  return b[ri[0]][ci[0]];
201  }
202 
203  template<typename RI, typename CI, typename Block>
204  typename Block::field_type&
205  access_matrix_element(tags::field_matrix_1_m, Block& b, const RI& ri, const CI& ci, int i, int j)
206  {
207  assert(i == -1);
208  assert(j == 0);
209  return b[0][ci[0]];
210  }
211 
212  template<typename RI, typename CI, typename Block>
213  typename Block::field_type&
214  access_matrix_element(tags::field_matrix_n_1, Block& b, const RI& ri, const CI& ci, int i, int j)
215  {
216  assert(i == 0);
217  assert(j == -1);
218  return b[ri[0]][0];
219  }
220 
221  template<typename RI, typename CI, typename Block>
222  typename Block::field_type&
223  access_matrix_element(tags::bcrs_matrix, Block& b, const RI& ri, const CI& ci, int i, int j)
224  {
225  return access_matrix_element(container_tag(b[ri[i]][ci[j]]),b[ri[i]][ci[j]],ri,ci,i-1,j-1);
226  }
227 
228 
229  template<typename RI, typename CI, typename Block>
230  const typename Block::field_type&
231  access_matrix_element(tags::field_matrix_1_1, const Block& b, const RI& ri, const CI& ci, int i, int j)
232  {
233  assert(i == -1);
234  assert(j == -1);
235  return b[0][0];
236  }
237 
238  template<typename RI, typename CI, typename Block>
239  const typename Block::field_type&
240  access_matrix_element(tags::field_matrix_n_m, const Block& b, const RI& ri, const CI& ci, int i, int j)
241  {
242  assert(i == 0);
243  assert(j == 0);
244  return b[ri[0]][ci[0]];
245  }
246 
247  template<typename RI, typename CI, typename Block>
248  const typename Block::field_type&
249  access_matrix_element(tags::field_matrix_1_m, const Block& b, const RI& ri, const CI& ci, int i, int j)
250  {
251  assert(i == -1);
252  assert(j == 0);
253  return b[0][ci[0]];
254  }
255 
256  template<typename RI, typename CI, typename Block>
257  const typename Block::field_type&
258  access_matrix_element(tags::field_matrix_n_1, const Block& b, const RI& ri, const CI& ci, int i, int j)
259  {
260  assert(i == 0);
261  assert(j == -1);
262  return b[ri[0]][0];
263  }
264 
265  template<typename RI, typename CI, typename Block>
266  const typename Block::field_type&
267  access_matrix_element(tags::bcrs_matrix, const Block& b, const RI& ri, const CI& ci, int i, int j)
268  {
269  return access_matrix_element(container_tag(b[ri[i]][ci[j]]),b[ri[i]][ci[j]],ri,ci,i-1,j-1);
270  }
271 
272 
273 
274  template<typename RI, typename Block>
275  void clear_matrix_row(tags::field_matrix_1_any, Block& b, const RI& ri, int i)
276  {
277  assert(i == -1);
278  b[0] = 0;
279  }
280 
281  template<typename RI, typename Block>
282  void clear_matrix_row(tags::field_matrix_n_any, Block& b, const RI& ri, int i)
283  {
284  assert(i == 0);
285  b[ri[0]] = 0;
286  }
287 
288  template<typename RI, typename Block>
289  void clear_matrix_row(tags::bcrs_matrix, Block& b, const RI& ri, int i)
290  {
291  typedef typename Block::ColIterator col_iterator_type;
292  const col_iterator_type end = b[ri[i]].end();
293  for(col_iterator_type cit = b[ri[i]].begin(); cit != end; ++cit)
294  clear_matrix_row(container_tag(*cit),*cit,ri,i-1);
295  }
296 
297 
298  template<typename RI, typename Block>
299  void clear_matrix_row_block(tags::field_matrix_1_1, Block& b, const RI& ri, int i)
300  {
301  assert(i == -1);
302  b = 0;
303  }
304 
305  template<typename RI, typename Block>
306  void clear_matrix_row_block(tags::field_matrix_1_any, Block& b, const RI& ri, int i)
307  {
308  DUNE_THROW(Dune::Exception,"Should never get here!");
309  }
310 
311  template<typename RI, typename Block>
312  void clear_matrix_row_block(tags::field_matrix_n_any, Block& b, const RI& ri, int i)
313  {
314  assert(i == 0);
315  b = 0;
316  }
317 
318  template<typename RI, typename Block>
319  void clear_matrix_row_block(tags::bcrs_matrix, Block& b, const RI& ri, int i)
320  {
321  typedef typename Block::ColIterator col_iterator_type;
322  const col_iterator_type end = b[ri[i]].end();
323  for(col_iterator_type cit = b[ri[i]].begin(); cit != end; ++cit)
324  clear_matrix_row_block(container_tag(*cit),*cit,ri,i-1);
325  }
326 
327 
328 
329  template<typename T, typename RI, typename CI, typename Block>
330  void write_matrix_element_if_exists(const T& v, tags::field_matrix_1_1, Block& b, const RI& ri, const CI& ci, int i, int j)
331  {
332  assert(i == -1);
333  assert(j == -1);
334  b[0][0] = v;
335  }
336 
337  template<typename T, typename RI, typename CI, typename Block>
338  void write_matrix_element_if_exists(const T& v, tags::field_matrix_n_m, Block& b, const RI& ri, const CI& ci, int i, int j)
339  {
340  assert(i == 0);
341  assert(j == 0);
342  b[ri[0]][ci[0]] = v;
343  }
344 
345  template<typename T, typename RI, typename CI, typename Block>
346  void write_matrix_element_if_exists(const T& v, tags::field_matrix_1_m, Block& b, const RI& ri, const CI& ci, int i, int j)
347  {
348  assert(i == -1);
349  assert(j == 0);
350  b[0][ci[0]] = v;
351  }
352 
353  template<typename T, typename RI, typename CI, typename Block>
354  void write_matrix_element_if_exists(const T& v, tags::field_matrix_n_1, Block& b, const RI& ri, const CI& ci, int i, int j)
355  {
356  assert(i == 0);
357  assert(j == -1);
358  b[ri[0]][0] = v;
359  }
360 
361  template<typename T, typename RI, typename CI, typename Block>
362  void write_matrix_element_if_exists(const T& v, tags::bcrs_matrix, Block& b, const RI& ri, const CI& ci, int i, int j)
363  {
364  if (b.exists(ri[i],ci[j]))
365  write_matrix_element_if_exists(v,container_tag(b[ri[i]][ci[j]]),b[ri[i]][ci[j]],ri,ci,i-1,j-1);
366  }
367 
368 
369 
370 
371  template<typename T, typename RI, typename CI, typename Block>
372  void write_matrix_element_if_exists_to_block(const T& v, tags::field_matrix_1_1, Block& b, const RI& ri, const CI& ci, int i, int j)
373  {
374  assert(i == -1);
375  assert(j == -1);
376  b[0][0] = v;
377  }
378 
379  template<typename T, typename RI, typename CI, typename Block>
380  void write_matrix_element_if_exists_to_block(const T& v, tags::field_matrix_n_m, Block& b, const RI& ri, const CI& ci, int i, int j)
381  {
382  assert(i == 0);
383  assert(j == 0);
384  for (std::size_t i = 0; i < b.rows; ++i)
385  b[i][i] = v;
386  }
387 
388  template<typename T, typename RI, typename CI, typename Block>
389  void write_matrix_element_if_exists_to_block(const T& v, tags::field_matrix_1_m, Block& b, const RI& ri, const CI& ci, int i, int j)
390  {
391  DUNE_THROW(Dune::Exception,"Should never get here!");
392  }
393 
394  template<typename T, typename RI, typename CI, typename Block>
395  void write_matrix_element_if_exists_to_block(const T& v, tags::field_matrix_n_1, Block& b, const RI& ri, const CI& ci, int i, int j)
396  {
397  DUNE_THROW(Dune::Exception,"Should never get here!");
398  }
399 
400  template<typename T, typename RI, typename CI, typename Block>
401  void write_matrix_element_if_exists_to_block(const T& v, tags::bcrs_matrix, Block& b, const RI& ri, const CI& ci, int i, int j)
402  {
403  if (b.exists(ri[i],ci[j]))
404  write_matrix_element_if_exists_to_block(v,container_tag(b[ri[i]][ci[j]]),b[ri[i]][ci[j]],ri,ci,i-1,j-1);
405  }
406 
407 
408  template<typename OrderingV, typename OrderingU, typename Pattern, typename Container>
409  typename std::enable_if<
412  >::type
413  allocate_matrix(const OrderingV& ordering_v,
414  const OrderingU& ordering_u,
415  const Pattern& p,
416  Container& c)
417  {
418  c.setSize(ordering_v.blockCount(),ordering_u.blockCount(),false);
419  c.setBuildMode(Container::random);
420 
421  for (std::size_t i = 0; i < c.N(); ++i)
422  c.setrowsize(i,p[i].size());
423  c.endrowsizes();
424 
425  for (std::size_t i = 0; i < c.N(); ++i)
426  for (typename Pattern::value_type::const_iterator cit = p[i].begin(); cit != p[i].end(); ++cit)
427  c.addindex(i,cit->first);
428  c.endindices();
429 
430  for (std::size_t i = 0; i < c.N(); ++i)
431  for (typename Pattern::value_type::const_iterator cit = p[i].begin(); cit != p[i].end(); ++cit)
432  {
433  allocate_matrix(ordering_v.childOrdering(i),
434  ordering_u.childOrdering(cit->first),
435  cit->second,
436  c[i][cit->first]);
437  }
438  }
439 
440  template<typename OrderingV, typename OrderingU, typename Pattern, typename Container>
441  typename std::enable_if<
444  >::type
445  allocate_matrix(const OrderingV& ordering_v,
446  const OrderingU& ordering_u,
447  const Pattern& p,
448  Container& c)
449  {
450  for (std::size_t i = 0; i < c.N(); ++i)
451  for (typename Pattern::value_type::iterator cit = p[i].begin(); cit != p[i].end(); ++cit)
452  {
453  allocate_matrix(ordering_v.childOrdering(i),
454  ordering_u.childOrdering(cit->first),
455  cit->second,
456  c[i][cit->first]);
457  }
458  }
459 
460  template<typename OrderingV, typename OrderingU, typename Pattern, typename Container>
461  typename std::enable_if<
463  >::type
464  allocate_matrix(const OrderingV& ordering_v,
465  const OrderingU& ordering_u,
466  const Pattern& p,
467  Container& c)
468  {
469  c.setSize(ordering_v.blockCount(),ordering_u.blockCount(),false);
470  c.setBuildMode(Container::random);
471 
472  for (std::size_t i = 0; i < c.N(); ++i)
473  c.setrowsize(i,p[i].size());
474  c.endrowsizes();
475 
476  for (std::size_t i = 0; i < c.N(); ++i)
477  for (typename Pattern::value_type::const_iterator cit = p[i].begin(); cit != p[i].end(); ++cit)
478  c.addindex(i,*cit);
479  c.endindices();
480  }
481 
482  } // namespace ISTL
483 
484 #endif // DOXYGEN
485 
486  } // namespace PDELab
487 } // namespace Dune
488 
489 #endif // DUNE_PDELAB_BACKEND_ISTL_MATRIXHELPERS_HH
tags::container< T >::type container_tag(const T &)
Gets instance of container tag associated with T.
Definition: backend/istl/tags.hh:246
static const unsigned int value
Definition: gridfunctionspace/tags.hh:139
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
const P & p
Definition: constraints.hh:147