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