3#ifndef DUNE_PDELAB_BACKEND_ISTL_MATRIXHELPERS_HH
4#define DUNE_PDELAB_BACKEND_ISTL_MATRIXHELPERS_HH
8#include <unordered_map>
9#include <unordered_set>
14#include <dune/pdelab/ordering/orderingbase.hh>
15#include <dune/pdelab/backend/istl/tags.hh>
24 template<
typename RV,
typename CV,
typename block_type>
25 struct matrix_for_vectors;
27 template<
typename B1,
typename A1,
typename B2,
typename A2,
typename block_type>
33 template<
typename B1,
int n1,
typename B2,
int n2,
typename block_type>
39 template<
typename E,
typename RV,
typename CV, std::
size_t blocklevel>
40 struct recursive_build_matrix_type
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;
45 template<
typename E,
typename RV,
typename CV>
46 struct recursive_build_matrix_type<E,RV,CV,1>
52 template<
typename E,
typename RV,
typename CV>
53 struct build_matrix_type
56#if DUNE_VERSION_LT(DUNE_ISTL,2,8)
57 static_assert(
static_cast<int>(RV::blocklevel) ==
static_cast<int>(CV::blocklevel),
"Both vectors must have identical blocking depth");
59 typedef typename recursive_build_matrix_type<E,RV,CV,RV::blocklevel>::type type;
61 static_assert(
static_cast<int>(blockLevel<RV>()) ==
static_cast<int>(blockLevel<CV>()),
"Both vectors must have identical blocking depth");
63 typedef typename recursive_build_matrix_type<E,RV,CV,blockLevel<RV>()>::type type;
68 template<
typename RowOrdering,
typename ColOrdering,
typename SubPattern_ =
void>
70 :
public std::vector<std::unordered_map<std::size_t,SubPattern_> >
75 typedef SubPattern_ SubPattern;
77 template<
typename RI,
typename CI>
78 void add_link(
const RI& ri,
const CI& ci)
80 recursive_add_entry(ri.view(),ci.view());
83 template<
typename RI,
typename CI>
84 void recursive_add_entry(
const RI& ri,
const CI& ci)
86 this->resize(_row_ordering.blockCount());
87 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()))));
88 r.first->second.recursive_add_entry(ri.back_popped(),ci.back_popped());
91 Pattern(
const RowOrdering& row_ordering,
const ColOrdering& col_ordering)
92 : _row_ordering(row_ordering)
93 , _col_ordering(col_ordering)
98 const RowOrdering& _row_ordering;
99 const ColOrdering& _col_ordering;
103 template<
typename RowOrdering,
typename ColOrdering>
104 class Pattern<RowOrdering,ColOrdering,void>
105 :
public std::vector<std::unordered_set<std::size_t> >
110 typedef void SubPattern;
112 template<
typename RI,
typename CI>
113 void add_link(
const RI& ri,
const CI& ci)
115 recursive_add_entry(ri,ci);
118 template<
typename RI,
typename CI>
119 void recursive_add_entry(
const RI& ri,
const CI& ci)
121 this->resize(_row_ordering.blockCount());
122 (*this)[ri.back()].insert(ci.back());
125 Pattern(
const RowOrdering& row_ordering,
const ColOrdering& col_ordering)
126 : _row_ordering(row_ordering)
127 , _col_ordering(col_ordering)
132 const RowOrdering& _row_ordering;
133 const ColOrdering& _col_ordering;
137#if DUNE_VERSION_LT(DUNE_ISTL,2,8)
138 template<
typename M,
int blocklevel = M::blocklevel>
140 template<typename M, int blocklevel = blockLevel<M>()>
142 struct requires_pattern
144 static const bool value = requires_pattern<
typename M::block_type,blocklevel-1>::value;
148 struct requires_pattern<M,0>
150 static const bool value =
false;
153 template<
typename B,
typename A,
int blocklevel>
154 struct requires_pattern<
Dune::BCRSMatrix<B,A>,blocklevel>
156 static const bool value =
true;
159 template<
typename M,
typename RowOrdering,
typename ColOrdering,
bool pattern>
160 struct _build_pattern_type
165 template<
typename M,
typename RowOrdering,
typename ColOrdering>
166 struct _build_pattern_type<M,RowOrdering,ColOrdering,true>
168 typedef Pattern<RowOrdering,ColOrdering,typename _build_pattern_type<typename M::block_type,RowOrdering,ColOrdering,requires_pattern<typename M::block_type>::value>::type> type;
171 template<
typename M,
typename GFSV,
typename GFSU,
typename Tag>
172 struct build_pattern_type
175 typedef OrderingBase<
176 typename GFSV::Ordering::Traits::DOFIndex,
177 typename GFSV::Ordering::Traits::ContainerIndex
180 typedef OrderingBase<
181 typename GFSU::Ordering::Traits::DOFIndex,
182 typename GFSU::Ordering::Traits::ContainerIndex
185 typedef typename _build_pattern_type<M,RowOrdering,ColOrdering,requires_pattern<M>::value>::type type;
188 template<
typename M,
typename GFSV,
typename GFSU>
189 struct build_pattern_type<M,GFSV,GFSU,FlatContainerAllocationTag>
191 typedef Pattern<typename GFSV::Ordering, typename GFSU::Ordering> type;
195 template<
typename RI,
typename CI,
typename Block>
196 typename Block::field_type&
197 access_matrix_element(tags::field_matrix_1_1, Block& b,
const RI& ri,
const CI& ci,
int i,
int j)
204 template<
typename RI,
typename CI,
typename Block>
205 typename Block::field_type&
206 access_matrix_element(tags::field_matrix_n_m, Block& b,
const RI& ri,
const CI& ci,
int i,
int j)
210 return b[ri[0]][ci[0]];
213 template<
typename RI,
typename CI,
typename Block>
214 typename Block::field_type&
215 access_matrix_element(tags::field_matrix_1_m, Block& b,
const RI& ri,
const CI& ci,
int i,
int j)
222 template<
typename RI,
typename CI,
typename Block>
223 typename Block::field_type&
224 access_matrix_element(tags::field_matrix_n_1, Block& b,
const RI& ri,
const CI& ci,
int i,
int j)
231 template<
typename RI,
typename CI,
typename Block>
232 typename Block::field_type&
233 access_matrix_element(tags::bcrs_matrix, Block& b,
const RI& ri,
const CI& ci,
int i,
int j)
235 return access_matrix_element(container_tag(b[ri[i]][ci[j]]),b[ri[i]][ci[j]],ri,ci,i-1,j-1);
239 template<
typename RI,
typename CI,
typename Block>
240 const typename Block::field_type&
241 access_matrix_element(tags::field_matrix_1_1,
const Block& b,
const RI& ri,
const CI& ci,
int i,
int j)
248 template<
typename RI,
typename CI,
typename Block>
249 const typename Block::field_type&
250 access_matrix_element(tags::field_matrix_n_m,
const Block& b,
const RI& ri,
const CI& ci,
int i,
int j)
254 return b[ri[0]][ci[0]];
257 template<
typename RI,
typename CI,
typename Block>
258 const typename Block::field_type&
259 access_matrix_element(tags::field_matrix_1_m,
const Block& b,
const RI& ri,
const CI& ci,
int i,
int j)
266 template<
typename RI,
typename CI,
typename Block>
267 const typename Block::field_type&
268 access_matrix_element(tags::field_matrix_n_1,
const Block& b,
const RI& ri,
const CI& ci,
int i,
int j)
275 template<
typename RI,
typename CI,
typename Block>
276 const typename Block::field_type&
277 access_matrix_element(tags::bcrs_matrix,
const Block& b,
const RI& ri,
const CI& ci,
int i,
int j)
279 return access_matrix_element(container_tag(b[ri[i]][ci[j]]),b[ri[i]][ci[j]],ri,ci,i-1,j-1);
284 template<
typename RI,
typename Block>
285 void clear_matrix_row(tags::field_matrix_1_any, Block& b,
const RI& ri,
int i)
291 template<
typename RI,
typename Block>
292 void clear_matrix_row(tags::field_matrix_n_any, Block& b,
const RI& ri,
int i)
298 template<
typename RI,
typename Block>
299 void clear_matrix_row(tags::bcrs_matrix, Block& b,
const RI& ri,
int i)
301 typedef typename Block::ColIterator col_iterator_type;
302 const col_iterator_type end = b[ri[i]].end();
303 for(col_iterator_type cit = b[ri[i]].begin(); cit != end; ++cit)
304 clear_matrix_row(container_tag(*cit),*cit,ri,i-1);
308 template<
typename RI,
typename Block>
309 void clear_matrix_row_block(tags::field_matrix_1_1, Block& b,
const RI& ri,
int i)
315 template<
typename RI,
typename Block>
316 void clear_matrix_row_block(tags::field_matrix_1_any, Block& b,
const RI& ri,
int i)
321 template<
typename RI,
typename Block>
322 void clear_matrix_row_block(tags::field_matrix_n_any, Block& b,
const RI& ri,
int i)
328 template<
typename RI,
typename Block>
329 void clear_matrix_row_block(tags::bcrs_matrix, Block& b,
const RI& ri,
int i)
331 typedef typename Block::ColIterator col_iterator_type;
332 const col_iterator_type end = b[ri[i]].end();
333 for(col_iterator_type cit = b[ri[i]].begin(); cit != end; ++cit)
334 clear_matrix_row_block(container_tag(*cit),*cit,ri,i-1);
339 template<
typename T,
typename RI,
typename CI,
typename Block>
340 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)
347 template<
typename T,
typename RI,
typename CI,
typename Block>
348 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)
355 template<
typename T,
typename RI,
typename CI,
typename Block>
356 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)
363 template<
typename T,
typename RI,
typename CI,
typename Block>
364 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)
371 template<
typename T,
typename RI,
typename CI,
typename Block>
372 void write_matrix_element_if_exists(
const T& v, tags::bcrs_matrix, Block& b,
const RI& ri,
const CI& ci,
int i,
int j)
374 if (b.exists(ri[i],ci[j]))
375 write_matrix_element_if_exists(v,container_tag(b[ri[i]][ci[j]]),b[ri[i]][ci[j]],ri,ci,i-1,j-1);
381 template<
typename T,
typename RI,
typename CI,
typename Block>
382 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)
389 template<
typename T,
typename RI,
typename CI,
typename Block>
390 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)
394 for (std::size_t i = 0; i < b.rows; ++i)
398 template<
typename T,
typename RI,
typename CI,
typename Block>
399 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)
404 template<
typename T,
typename RI,
typename CI,
typename Block>
405 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)
410 template<
typename T,
typename RI,
typename CI,
typename Block>
411 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)
413 if (b.exists(ri[i],ci[j]))
414 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);
418 template<
typename OrderingV,
typename OrderingU,
typename Pattern,
typename Container>
419 typename std::enable_if<
420 !std::is_same<typename Pattern::SubPattern,void>::value &&
421 requires_pattern<Container>::value
423 allocate_matrix(
const OrderingV& ordering_v,
424 const OrderingU& ordering_u,
428 c.setSize(ordering_v.blockCount(),ordering_u.blockCount(),
false);
429 c.setBuildMode(Container::random);
431 for (std::size_t i = 0; i < c.N(); ++i)
432 c.setrowsize(i,p[i].size());
435 for (std::size_t i = 0; i < c.N(); ++i)
436 for (
typename Pattern::value_type::const_iterator cit = p[i].begin(); cit != p[i].end(); ++cit)
437 c.addindex(i,cit->first);
440 for (std::size_t i = 0; i < c.N(); ++i)
441 for (
typename Pattern::value_type::const_iterator cit = p[i].begin(); cit != p[i].end(); ++cit)
443 allocate_matrix(ordering_v.childOrdering(i),
444 ordering_u.childOrdering(cit->first),
450 template<
typename OrderingV,
typename OrderingU,
typename Pattern,
typename Container>
451 typename std::enable_if<
452 !std::is_same<typename Pattern::SubPattern,void>::value &&
453 !requires_pattern<Container>::value
455 allocate_matrix(
const OrderingV& ordering_v,
456 const OrderingU& ordering_u,
460 for (std::size_t i = 0; i < c.N(); ++i)
461 for (
typename Pattern::value_type::iterator cit = p[i].begin(); cit != p[i].end(); ++cit)
463 allocate_matrix(ordering_v.childOrdering(i),
464 ordering_u.childOrdering(cit->first),
470 template<
typename OrderingV,
typename OrderingU,
typename Pattern,
typename Container>
471 typename std::enable_if<
472 std::is_same<typename Pattern::SubPattern,void>::value
474 allocate_matrix(
const OrderingV& ordering_v,
475 const OrderingU& ordering_u,
479 c.setSize(ordering_v.blockCount(),ordering_u.blockCount(),
false);
480 c.setBuildMode(Container::random);
482 for (std::size_t i = 0; i < c.N(); ++i)
483 c.setrowsize(i,p[i].size());
486 for (std::size_t i = 0; i < c.N(); ++i)
487 for (
typename Pattern::value_type::const_iterator cit = p[i].begin(); cit != p[i].end(); ++cit)
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:466
A vector of blocks with memory management.
Definition: bvector.hh:392
Base class for Dune-Exceptions.
Definition: exceptions.hh:96
A dense n x m matrix.
Definition: fmatrix.hh:117
vector space out of a tensor product of fields.
Definition: fvector.hh:91
Implementation of the BCRSMatrix class.
Implements a matrix constructed from a given type representing a field and compile-time given number ...
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
Dune namespace.
Definition: alignedallocator.hh:13