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 14 #include <dune/pdelab/ordering/orderingbase.hh> 15 #include <dune/pdelab/backend/istl/tags.hh> 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 #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"); 58 59 typedef typename recursive_build_matrix_type<E,RV,CV,RV::blocklevel>::type type; 60 #else 61 static_assert(static_cast<int>(blockLevel<RV>()) == static_cast<int>(blockLevel<CV>()),"Both vectors must have identical blocking depth"); 62 63 typedef typename recursive_build_matrix_type<E,RV,CV,blockLevel<RV>()>::type type; 64 #endif 65 66 }; 67 68 template<typename RowOrdering, typename ColOrdering, typename SubPattern_ = void> 69 class Pattern 70 : public std::vector<std::unordered_map<std::size_t,SubPattern_> > 71 { 72 73 public: 74 75 typedef SubPattern_ SubPattern; 76 77 template<typename RI, typename CI> add_link(const RI & ri,const CI & ci)78 void add_link(const RI& ri, const CI& ci) 79 { 80 recursive_add_entry(ri.view(),ci.view()); 81 } 82 83 template<typename RI, typename CI> recursive_add_entry(const RI & ri,const CI & ci)84 void recursive_add_entry(const RI& ri, const CI& ci) 85 { 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()); 89 } 90 Pattern(const RowOrdering & row_ordering,const ColOrdering & col_ordering)91 Pattern(const RowOrdering& row_ordering, const ColOrdering& col_ordering) 92 : _row_ordering(row_ordering) 93 , _col_ordering(col_ordering) 94 {} 95 96 private: 97 98 const RowOrdering& _row_ordering; 99 const ColOrdering& _col_ordering; 100 101 }; 102 103 template<typename RowOrdering, typename ColOrdering> 104 class Pattern<RowOrdering,ColOrdering,void> 105 : public std::vector<std::unordered_set<std::size_t> > 106 { 107 108 public: 109 110 typedef void SubPattern; 111 112 template<typename RI, typename CI> add_link(const RI & ri,const CI & ci)113 void add_link(const RI& ri, const CI& ci) 114 { 115 recursive_add_entry(ri,ci); 116 } 117 118 template<typename RI, typename CI> recursive_add_entry(const RI & ri,const CI & ci)119 void recursive_add_entry(const RI& ri, const CI& ci) 120 { 121 this->resize(_row_ordering.blockCount()); 122 (*this)[ri.back()].insert(ci.back()); 123 } 124 Pattern(const RowOrdering & row_ordering,const ColOrdering & col_ordering)125 Pattern(const RowOrdering& row_ordering, const ColOrdering& col_ordering) 126 : _row_ordering(row_ordering) 127 , _col_ordering(col_ordering) 128 {} 129 130 private: 131 132 const RowOrdering& _row_ordering; 133 const ColOrdering& _col_ordering; 134 135 }; 136 137 #if DUNE_VERSION_LT(DUNE_ISTL,2,8) 138 template<typename M, int blocklevel = M::blocklevel> 139 #else 140 template<typename M, int blocklevel = blockLevel<M>()> 141 #endif 142 struct requires_pattern 143 { 144 static const bool value = requires_pattern<typename M::block_type,blocklevel-1>::value; 145 }; 146 147 template<typename M> 148 struct requires_pattern<M,0> 149 { 150 static const bool value = false; 151 }; 152 153 template<typename B, typename A, int blocklevel> 154 struct requires_pattern<Dune::BCRSMatrix<B,A>,blocklevel> 155 { 156 static const bool value = true; 157 }; 158 159 template<typename M, typename RowOrdering, typename ColOrdering, bool pattern> 160 struct _build_pattern_type 161 { 162 typedef void type; 163 }; 164 165 template<typename M, typename RowOrdering, typename ColOrdering> 166 struct _build_pattern_type<M,RowOrdering,ColOrdering,true> 167 { 168 typedef Pattern<RowOrdering,ColOrdering,typename _build_pattern_type<typename M::block_type,RowOrdering,ColOrdering,requires_pattern<typename M::block_type>::value>::type> type; 169 }; 170 171 template<typename M, typename GFSV, typename GFSU, typename Tag> 172 struct build_pattern_type 173 { 174 175 typedef OrderingBase< 176 typename GFSV::Ordering::Traits::DOFIndex, 177 typename GFSV::Ordering::Traits::ContainerIndex 178 > RowOrdering; 179 180 typedef OrderingBase< 181 typename GFSU::Ordering::Traits::DOFIndex, 182 typename GFSU::Ordering::Traits::ContainerIndex 183 > ColOrdering; 184 185 typedef typename _build_pattern_type<M,RowOrdering,ColOrdering,requires_pattern<M>::value>::type type; 186 }; 187 188 template<typename M, typename GFSV, typename GFSU> 189 struct build_pattern_type<M,GFSV,GFSU,FlatContainerAllocationTag> 190 { 191 typedef Pattern<typename GFSV::Ordering, typename GFSU::Ordering> type; 192 }; 193 194 195 template<typename RI, typename CI, typename Block> 196 typename Block::field_type& access_matrix_element(tags::field_matrix_1_1,Block & b,const RI & ri,const CI & ci,int i,int j)197 access_matrix_element(tags::field_matrix_1_1, Block& b, const RI& ri, const CI& ci, int i, int j) 198 { 199 assert(i == -1); 200 assert(j == -1); 201 return b[0][0]; 202 } 203 204 template<typename RI, typename CI, typename Block> 205 typename Block::field_type& access_matrix_element(tags::field_matrix_n_m,Block & b,const RI & ri,const CI & ci,int i,int j)206 access_matrix_element(tags::field_matrix_n_m, Block& b, const RI& ri, const CI& ci, int i, int j) 207 { 208 assert(i == 0); 209 assert(j == 0); 210 return b[ri[0]][ci[0]]; 211 } 212 213 template<typename RI, typename CI, typename Block> 214 typename Block::field_type& access_matrix_element(tags::field_matrix_1_m,Block & b,const RI & ri,const CI & ci,int i,int j)215 access_matrix_element(tags::field_matrix_1_m, Block& b, const RI& ri, const CI& ci, int i, int j) 216 { 217 assert(i == -1); 218 assert(j == 0); 219 return b[0][ci[0]]; 220 } 221 222 template<typename RI, typename CI, typename Block> 223 typename Block::field_type& access_matrix_element(tags::field_matrix_n_1,Block & b,const RI & ri,const CI & ci,int i,int j)224 access_matrix_element(tags::field_matrix_n_1, Block& b, const RI& ri, const CI& ci, int i, int j) 225 { 226 assert(i == 0); 227 assert(j == -1); 228 return b[ri[0]][0]; 229 } 230 231 template<typename RI, typename CI, typename Block> 232 typename Block::field_type& access_matrix_element(tags::bcrs_matrix,Block & b,const RI & ri,const CI & ci,int i,int j)233 access_matrix_element(tags::bcrs_matrix, Block& b, const RI& ri, const CI& ci, int i, int j) 234 { 235 return access_matrix_element(container_tag(b[ri[i]][ci[j]]),b[ri[i]][ci[j]],ri,ci,i-1,j-1); 236 } 237 238 239 template<typename RI, typename CI, typename Block> 240 const typename Block::field_type& access_matrix_element(tags::field_matrix_1_1,const Block & b,const RI & ri,const CI & ci,int i,int j)241 access_matrix_element(tags::field_matrix_1_1, const Block& b, const RI& ri, const CI& ci, int i, int j) 242 { 243 assert(i == -1); 244 assert(j == -1); 245 return b[0][0]; 246 } 247 248 template<typename RI, typename CI, typename Block> 249 const typename Block::field_type& access_matrix_element(tags::field_matrix_n_m,const Block & b,const RI & ri,const CI & ci,int i,int j)250 access_matrix_element(tags::field_matrix_n_m, const Block& b, const RI& ri, const CI& ci, int i, int j) 251 { 252 assert(i == 0); 253 assert(j == 0); 254 return b[ri[0]][ci[0]]; 255 } 256 257 template<typename RI, typename CI, typename Block> 258 const typename Block::field_type& access_matrix_element(tags::field_matrix_1_m,const Block & b,const RI & ri,const CI & ci,int i,int j)259 access_matrix_element(tags::field_matrix_1_m, const Block& b, const RI& ri, const CI& ci, int i, int j) 260 { 261 assert(i == -1); 262 assert(j == 0); 263 return b[0][ci[0]]; 264 } 265 266 template<typename RI, typename CI, typename Block> 267 const typename Block::field_type& access_matrix_element(tags::field_matrix_n_1,const Block & b,const RI & ri,const CI & ci,int i,int j)268 access_matrix_element(tags::field_matrix_n_1, const Block& b, const RI& ri, const CI& ci, int i, int j) 269 { 270 assert(i == 0); 271 assert(j == -1); 272 return b[ri[0]][0]; 273 } 274 275 template<typename RI, typename CI, typename Block> 276 const typename Block::field_type& access_matrix_element(tags::bcrs_matrix,const Block & b,const RI & ri,const CI & ci,int i,int j)277 access_matrix_element(tags::bcrs_matrix, const Block& b, const RI& ri, const CI& ci, int i, int j) 278 { 279 return access_matrix_element(container_tag(b[ri[i]][ci[j]]),b[ri[i]][ci[j]],ri,ci,i-1,j-1); 280 } 281 282 283 284 template<typename RI, typename Block> clear_matrix_row(tags::field_matrix_1_any,Block & b,const RI & ri,int i)285 void clear_matrix_row(tags::field_matrix_1_any, Block& b, const RI& ri, int i) 286 { 287 assert(i == -1); 288 b[0] = 0; 289 } 290 291 template<typename RI, typename Block> clear_matrix_row(tags::field_matrix_n_any,Block & b,const RI & ri,int i)292 void clear_matrix_row(tags::field_matrix_n_any, Block& b, const RI& ri, int i) 293 { 294 assert(i == 0); 295 b[ri[0]] = 0; 296 } 297 298 template<typename RI, typename Block> clear_matrix_row(tags::bcrs_matrix,Block & b,const RI & ri,int i)299 void clear_matrix_row(tags::bcrs_matrix, Block& b, const RI& ri, int i) 300 { 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); 305 } 306 307 308 template<typename RI, typename Block> clear_matrix_row_block(tags::field_matrix_1_1,Block & b,const RI & ri,int i)309 void clear_matrix_row_block(tags::field_matrix_1_1, Block& b, const RI& ri, int i) 310 { 311 assert(i == -1); 312 b = 0; 313 } 314 315 template<typename RI, typename Block> clear_matrix_row_block(tags::field_matrix_1_any,Block & b,const RI & ri,int i)316 void clear_matrix_row_block(tags::field_matrix_1_any, Block& b, const RI& ri, int i) 317 { 318 DUNE_THROW(Dune::Exception,"Should never get here!"); 319 } 320 321 template<typename RI, typename Block> clear_matrix_row_block(tags::field_matrix_n_any,Block & b,const RI & ri,int i)322 void clear_matrix_row_block(tags::field_matrix_n_any, Block& b, const RI& ri, int i) 323 { 324 assert(i == 0); 325 b = 0; 326 } 327 328 template<typename RI, typename Block> clear_matrix_row_block(tags::bcrs_matrix,Block & b,const RI & ri,int i)329 void clear_matrix_row_block(tags::bcrs_matrix, Block& b, const RI& ri, int i) 330 { 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); 335 } 336 337 338 339 template<typename T, typename RI, typename CI, typename Block> 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)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) 341 { 342 assert(i == -1); 343 assert(j == -1); 344 b[0][0] = v; 345 } 346 347 template<typename T, typename RI, typename CI, typename Block> 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)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) 349 { 350 assert(i == 0); 351 assert(j == 0); 352 b[ri[0]][ci[0]] = v; 353 } 354 355 template<typename T, typename RI, typename CI, typename Block> 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)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) 357 { 358 assert(i == -1); 359 assert(j == 0); 360 b[0][ci[0]] = v; 361 } 362 363 template<typename T, typename RI, typename CI, typename Block> 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)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) 365 { 366 assert(i == 0); 367 assert(j == -1); 368 b[ri[0]][0] = v; 369 } 370 371 template<typename T, typename RI, typename CI, typename Block> write_matrix_element_if_exists(const T & v,tags::bcrs_matrix,Block & b,const RI & ri,const CI & ci,int i,int j)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) 373 { 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); 376 } 377 378 379 380 381 template<typename T, typename RI, typename CI, typename Block> 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)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) 383 { 384 assert(i == -1); 385 assert(j == -1); 386 b[0][0] = v; 387 } 388 389 template<typename T, typename RI, typename CI, typename Block> 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)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) 391 { 392 assert(i == 0); 393 assert(j == 0); 394 for (std::size_t i = 0; i < b.rows; ++i) 395 b[i][i] = v; 396 } 397 398 template<typename T, typename RI, typename CI, typename Block> 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)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) 400 { 401 DUNE_THROW(Dune::Exception,"Should never get here!"); 402 } 403 404 template<typename T, typename RI, typename CI, typename Block> 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)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) 406 { 407 DUNE_THROW(Dune::Exception,"Should never get here!"); 408 } 409 410 template<typename T, typename RI, typename CI, typename Block> 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)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) 412 { 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); 415 } 416 417 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 422 >::type allocate_matrix(const OrderingV & ordering_v,const OrderingU & ordering_u,const Pattern & p,Container & c)423 allocate_matrix(const OrderingV& ordering_v, 424 const OrderingU& ordering_u, 425 const Pattern& p, 426 Container& c) 427 { 428 c.setSize(ordering_v.blockCount(),ordering_u.blockCount(),false); 429 c.setBuildMode(Container::random); 430 431 for (std::size_t i = 0; i < c.N(); ++i) 432 c.setrowsize(i,p[i].size()); 433 c.endrowsizes(); 434 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); 438 c.endindices(); 439 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) 442 { 443 allocate_matrix(ordering_v.childOrdering(i), 444 ordering_u.childOrdering(cit->first), 445 cit->second, 446 c[i][cit->first]); 447 } 448 } 449 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 454 >::type allocate_matrix(const OrderingV & ordering_v,const OrderingU & ordering_u,const Pattern & p,Container & c)455 allocate_matrix(const OrderingV& ordering_v, 456 const OrderingU& ordering_u, 457 const Pattern& p, 458 Container& c) 459 { 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) 462 { 463 allocate_matrix(ordering_v.childOrdering(i), 464 ordering_u.childOrdering(cit->first), 465 cit->second, 466 c[i][cit->first]); 467 } 468 } 469 470 template<typename OrderingV, typename OrderingU, typename Pattern, typename Container> 471 typename std::enable_if< 472 std::is_same<typename Pattern::SubPattern,void>::value 473 >::type allocate_matrix(const OrderingV & ordering_v,const OrderingU & ordering_u,const Pattern & p,Container & c)474 allocate_matrix(const OrderingV& ordering_v, 475 const OrderingU& ordering_u, 476 const Pattern& p, 477 Container& c) 478 { 479 c.setSize(ordering_v.blockCount(),ordering_u.blockCount(),false); 480 c.setBuildMode(Container::random); 481 482 for (std::size_t i = 0; i < c.N(); ++i) 483 c.setrowsize(i,p[i].size()); 484 c.endrowsizes(); 485 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) 488 c.addindex(i,*cit); 489 c.endindices(); 490 } 491 492 } // namespace ISTL 493 494 #endif // DOXYGEN 495 496 } // namespace PDELab 497 } // namespace Dune 498 499 #endif // DUNE_PDELAB_BACKEND_ISTL_MATRIXHELPERS_HH 500