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_ISTL_MATRIXREDISTRIBUTE_HH 4 #define DUNE_ISTL_MATRIXREDISTRIBUTE_HH 5 #include <memory> 6 #include "repartition.hh" 7 #include <dune/common/exceptions.hh> 8 #include <dune/common/parallel/indexset.hh> 9 #include <dune/istl/owneroverlapcopy.hh> 10 #include <dune/istl/paamg/pinfo.hh> 11 /** 12 * @file 13 * @brief Functionality for redistributing a sparse matrix. 14 * @author Markus Blatt 15 */ 16 namespace Dune 17 { 18 template<typename T> 19 struct RedistributeInformation 20 { isSetupDune::RedistributeInformation21 bool isSetup() const 22 { 23 return false; 24 } 25 template<class D> redistributeDune::RedistributeInformation26 void redistribute([[maybe_unused]] const D& from, [[maybe_unused]] D& to) const 27 {} 28 29 template<class D> redistributeBackwardDune::RedistributeInformation30 void redistributeBackward([[maybe_unused]] D& from, [[maybe_unused]]const D& to) const 31 {} 32 resetSetupDune::RedistributeInformation33 void resetSetup() 34 {} 35 setNoRowsDune::RedistributeInformation36 void setNoRows([[maybe_unused]] std::size_t size) 37 {} 38 setNoCopyRowsDune::RedistributeInformation39 void setNoCopyRows([[maybe_unused]] std::size_t size) 40 {} 41 setNoBackwardsCopyRowsDune::RedistributeInformation42 void setNoBackwardsCopyRows([[maybe_unused]] std::size_t size) 43 {} 44 getRowSizeDune::RedistributeInformation45 std::size_t getRowSize([[maybe_unused]] std::size_t index) const 46 { 47 return -1; 48 } 49 getCopyRowSizeDune::RedistributeInformation50 std::size_t getCopyRowSize([[maybe_unused]] std::size_t index) const 51 { 52 return -1; 53 } 54 getBackwardsCopyRowSizeDune::RedistributeInformation55 std::size_t getBackwardsCopyRowSize([[maybe_unused]] std::size_t index) const 56 { 57 return -1; 58 } 59 60 }; 61 62 #if HAVE_MPI 63 template<typename T, typename T1> 64 class RedistributeInformation<OwnerOverlapCopyCommunication<T,T1> > 65 { 66 public: 67 typedef OwnerOverlapCopyCommunication<T,T1> Comm; 68 RedistributeInformation()69 RedistributeInformation() 70 : interface(), setup_(false) 71 {} 72 getInterface()73 RedistributeInterface& getInterface() 74 { 75 return interface; 76 } 77 template<typename IS> checkInterface(const IS & source,const IS & target,MPI_Comm comm)78 void checkInterface(const IS& source, 79 const IS& target, MPI_Comm comm) 80 { 81 auto ri = std::make_unique<RemoteIndices<IS> >(source, target, comm); 82 ri->template rebuild<true>(); 83 Interface inf; 84 typename OwnerOverlapCopyCommunication<int>::OwnerSet flags; 85 int rank; 86 MPI_Comm_rank(MPI_COMM_WORLD, &rank); 87 inf.free(); 88 inf.build(*ri, flags, flags); 89 90 91 #ifdef DEBUG_REPART 92 if(inf!=interface) { 93 94 MPI_Comm_rank(MPI_COMM_WORLD, &rank); 95 if(rank==0) 96 std::cout<<"Interfaces do not match!"<<std::endl; 97 std::cout<<rank<<": redist interface new :"<<inf<<std::endl; 98 std::cout<<rank<<": redist interface :"<<interface<<std::endl; 99 100 throw "autsch!"; 101 } 102 #endif 103 } setSetup()104 void setSetup() 105 { 106 setup_=true; 107 interface.strip(); 108 } 109 resetSetup()110 void resetSetup() 111 { 112 setup_=false; 113 } 114 115 template<class GatherScatter, class D> redistribute(const D & from,D & to) const116 void redistribute(const D& from, D& to) const 117 { 118 BufferedCommunicator communicator; 119 communicator.template build<D>(from,to, interface); 120 communicator.template forward<GatherScatter>(from, to); 121 communicator.free(); 122 } 123 template<class GatherScatter, class D> redistributeBackward(D & from,const D & to) const124 void redistributeBackward(D& from, const D& to) const 125 { 126 127 BufferedCommunicator communicator; 128 communicator.template build<D>(from,to, interface); 129 communicator.template backward<GatherScatter>(from, to); 130 communicator.free(); 131 } 132 133 template<class D> redistribute(const D & from,D & to) const134 void redistribute(const D& from, D& to) const 135 { 136 redistribute<CopyGatherScatter<D> >(from,to); 137 } 138 template<class D> redistributeBackward(D & from,const D & to) const139 void redistributeBackward(D& from, const D& to) const 140 { 141 redistributeBackward<CopyGatherScatter<D> >(from,to); 142 } isSetup() const143 bool isSetup() const 144 { 145 return setup_; 146 } 147 reserve(std::size_t size)148 void reserve(std::size_t size) 149 {} 150 getRowSize(std::size_t index)151 std::size_t& getRowSize(std::size_t index) 152 { 153 return rowSize[index]; 154 } 155 getRowSize(std::size_t index) const156 std::size_t getRowSize(std::size_t index) const 157 { 158 return rowSize[index]; 159 } 160 getCopyRowSize(std::size_t index)161 std::size_t& getCopyRowSize(std::size_t index) 162 { 163 return copyrowSize[index]; 164 } 165 getCopyRowSize(std::size_t index) const166 std::size_t getCopyRowSize(std::size_t index) const 167 { 168 return copyrowSize[index]; 169 } 170 getBackwardsCopyRowSize(std::size_t index)171 std::size_t& getBackwardsCopyRowSize(std::size_t index) 172 { 173 return backwardscopyrowSize[index]; 174 } 175 getBackwardsCopyRowSize(std::size_t index) const176 std::size_t getBackwardsCopyRowSize(std::size_t index) const 177 { 178 return backwardscopyrowSize[index]; 179 } 180 setNoRows(std::size_t rows)181 void setNoRows(std::size_t rows) 182 { 183 rowSize.resize(rows, 0); 184 } 185 setNoCopyRows(std::size_t rows)186 void setNoCopyRows(std::size_t rows) 187 { 188 copyrowSize.resize(rows, 0); 189 } 190 setNoBackwardsCopyRows(std::size_t rows)191 void setNoBackwardsCopyRows(std::size_t rows) 192 { 193 backwardscopyrowSize.resize(rows, 0); 194 } 195 196 private: 197 std::vector<std::size_t> rowSize; 198 std::vector<std::size_t> copyrowSize; 199 std::vector<std::size_t> backwardscopyrowSize; 200 RedistributeInterface interface; 201 bool setup_; 202 }; 203 204 /** 205 * @brief Utility class to communicate and set the row sizes 206 * of a redistributed matrix. 207 * 208 * @tparam M The type of the matrix that the row size 209 * is communicated of. 210 * @tparam RI The type of class for redistribution information 211 */ 212 template<class M, class RI> 213 struct CommMatrixRowSize 214 { 215 // Make the default communication policy work. 216 typedef typename M::size_type value_type; 217 typedef typename M::size_type size_type; 218 219 /** 220 * @brief Constructor. 221 * @param m_ The matrix whose sparsity pattern is communicated. 222 * @param[out] rowsize_ RedistributeInformation object 223 */ CommMatrixRowSizeDune::CommMatrixRowSize224 CommMatrixRowSize(const M& m_, RI& rowsize_) 225 : matrix(m_), rowsize(rowsize_) 226 {} 227 const M& matrix; 228 RI& rowsize; 229 230 }; 231 232 233 /** 234 * @brief Utility class to communicate and build the sparsity pattern 235 * of a redistributed matrix. 236 * 237 * @tparam M The type of the matrix that the sparsity pattern 238 * is communicated of. 239 * @tparam I The type of the index set. 240 */ 241 template<class M, class I> 242 struct CommMatrixSparsityPattern 243 { 244 typedef typename M::size_type size_type; 245 246 /** 247 * @brief Constructor for the original side 248 * @param m_ The matrix whose sparsity pattern is communicated. 249 * @param idxset_ The index set corresponding to the local matrix. 250 * @param aggidxset_ The index set corresponding to the redistributed matrix. 251 */ CommMatrixSparsityPatternDune::CommMatrixSparsityPattern252 CommMatrixSparsityPattern(const M& m_, const Dune::GlobalLookupIndexSet<I>& idxset_, const I& aggidxset_) 253 : matrix(m_), idxset(idxset_), aggidxset(aggidxset_), rowsize() 254 {} 255 256 /** 257 * @brief Constructor for the redistruted side. 258 * @param m_ The matrix whose sparsity pattern is communicated. 259 * @param idxset_ The index set corresponding to the local matrix. 260 * @param aggidxset_ The index set corresponding to the redistributed matrix. 261 * @param rowsize_ The row size for the redistributed owner rows. 262 */ CommMatrixSparsityPatternDune::CommMatrixSparsityPattern263 CommMatrixSparsityPattern(const M& m_, const Dune::GlobalLookupIndexSet<I>& idxset_, const I& aggidxset_, 264 const std::vector<typename M::size_type>& rowsize_) 265 : matrix(m_), idxset(idxset_), aggidxset(aggidxset_), sparsity(aggidxset_.size()), rowsize(&rowsize_) 266 {} 267 268 /** 269 * @brief Creates and stores the sparsity pattern of the redistributed matrix. 270 * 271 * After the pattern is communicated this function can be used. 272 * @param m The matrix to build. 273 */ storeSparsityPatternDune::CommMatrixSparsityPattern274 void storeSparsityPattern(M& m) 275 { 276 // insert diagonal to overlap rows 277 typedef typename Dune::GlobalLookupIndexSet<I>::const_iterator IIter; 278 typedef typename Dune::OwnerOverlapCopyCommunication<int>::OwnerSet OwnerSet; 279 std::size_t nnz=0; 280 #ifdef DEBUG_REPART 281 int rank; 282 283 MPI_Comm_rank(MPI_COMM_WORLD, &rank); 284 #endif 285 for(IIter i= aggidxset.begin(), end=aggidxset.end(); i!=end; ++i) { 286 if(!OwnerSet::contains(i->local().attribute())) { 287 #ifdef DEBUG_REPART 288 std::cout<<rank<<" Inserting diagonal for"<<i->local()<<std::endl; 289 #endif 290 sparsity[i->local()].insert(i->local()); 291 } 292 293 nnz+=sparsity[i->local()].size(); 294 } 295 assert( aggidxset.size()==sparsity.size()); 296 297 if(nnz>0) { 298 m.setSize(aggidxset.size(), aggidxset.size(), nnz); 299 m.setBuildMode(M::row_wise); 300 typename M::CreateIterator citer=m.createbegin(); 301 #ifdef DEBUG_REPART 302 std::size_t idx=0; 303 bool correct=true; 304 Dune::GlobalLookupIndexSet<I> global(aggidxset); 305 #endif 306 typedef typename std::vector<std::set<size_type> >::const_iterator Iter; 307 for(Iter i=sparsity.begin(), end=sparsity.end(); i!=end; ++i, ++citer) 308 { 309 typedef typename std::set<size_type>::const_iterator SIter; 310 for(SIter si=i->begin(), send=i->end(); si!=send; ++si) 311 citer.insert(*si); 312 #ifdef DEBUG_REPART 313 if(i->find(idx)==i->end()) { 314 const typename I::IndexPair* gi=global.pair(idx); 315 assert(gi); 316 std::cout<<rank<<": row "<<idx<<" is missing a diagonal entry! global="<<gi->global()<<" attr="<<gi->local().attribute()<<" "<< 317 OwnerSet::contains(gi->local().attribute())<< 318 " row size="<<i->size()<<std::endl; 319 correct=false; 320 } 321 ++idx; 322 #endif 323 } 324 #ifdef DEBUG_REPART 325 if(!correct) 326 throw "bla"; 327 #endif 328 } 329 } 330 331 /** 332 * @brief Completes the sparsity pattern of the redistributed matrix with data 333 * from copy rows for the novlp case. 334 * 335 * After the pattern communication this function can be used. 336 * @param add_sparsity Sparsity pattern from the copy rows. 337 */ completeSparsityPatternDune::CommMatrixSparsityPattern338 void completeSparsityPattern(std::vector<std::set<size_type> > add_sparsity) 339 { 340 for (unsigned int i = 0; i != sparsity.size(); ++i) { 341 if (add_sparsity[i].size() != 0) { 342 typedef std::set<size_type> Set; 343 Set tmp_set; 344 std::insert_iterator<Set> tmp_insert (tmp_set, tmp_set.begin()); 345 std::set_union(add_sparsity[i].begin(), add_sparsity[i].end(), 346 sparsity[i].begin(), sparsity[i].end(), tmp_insert); 347 sparsity[i].swap(tmp_set); 348 } 349 } 350 } 351 352 const M& matrix; 353 typedef Dune::GlobalLookupIndexSet<I> LookupIndexSet; 354 const Dune::GlobalLookupIndexSet<I>& idxset; 355 const I& aggidxset; 356 std::vector<std::set<size_type> > sparsity; 357 const std::vector<size_type>* rowsize; 358 }; 359 360 template<class M, class I> 361 struct CommPolicy<CommMatrixSparsityPattern<M,I> > 362 { 363 typedef CommMatrixSparsityPattern<M,I> Type; 364 365 /** 366 * @brief The indexed type we send. 367 * This is the global index indentitfying the column. 368 */ 369 typedef typename I::GlobalIndex IndexedType; 370 371 /** @brief Each row varies in size. */ 372 typedef VariableSize IndexedTypeFlag; 373 getSizeDune::CommPolicy374 static typename M::size_type getSize(const Type& t, std::size_t i) 375 { 376 if(!t.rowsize) 377 return t.matrix[i].size(); 378 else 379 { 380 assert((*t.rowsize)[i]>0); 381 return (*t.rowsize)[i]; 382 } 383 } 384 }; 385 386 /** 387 * @brief Utility class for comunicating the matrix entries. 388 * 389 * @tparam M The type of the matrix. 390 * @tparam I The type of the ParallelIndexSet. 391 */ 392 template<class M, class I> 393 struct CommMatrixRow 394 { 395 /** 396 * @brief Constructor. 397 * @param m_ The matrix to communicate the values. That is the local original matrix 398 * as the source of the communication and the redistributed at the target of the 399 * communication. 400 * @param idxset_ The index set for the original matrix. 401 * @param aggidxset_ The index set for the redistributed matrix. 402 */ CommMatrixRowDune::CommMatrixRow403 CommMatrixRow(M& m_, const Dune::GlobalLookupIndexSet<I>& idxset_, const I& aggidxset_) 404 : matrix(m_), idxset(idxset_), aggidxset(aggidxset_), rowsize() 405 {} 406 407 /** 408 * @brief Constructor. 409 */ CommMatrixRowDune::CommMatrixRow410 CommMatrixRow(M& m_, const Dune::GlobalLookupIndexSet<I>& idxset_, const I& aggidxset_, 411 std::vector<size_t>& rowsize_) 412 : matrix(m_), idxset(idxset_), aggidxset(aggidxset_), rowsize(&rowsize_) 413 {} 414 /** 415 * @brief Sets the non-owner rows correctly as Dirichlet boundaries. 416 * 417 * This should be called after the communication. 418 */ setOverlapRowsToDirichletDune::CommMatrixRow419 void setOverlapRowsToDirichlet() 420 { 421 typedef typename Dune::GlobalLookupIndexSet<I>::const_iterator Iter; 422 typedef typename Dune::OwnerOverlapCopyCommunication<int>::OwnerSet OwnerSet; 423 424 for(Iter i= aggidxset.begin(), end=aggidxset.end(); i!=end; ++i) 425 if(!OwnerSet::contains(i->local().attribute())) { 426 // Set to Dirchlet 427 typedef typename M::ColIterator CIter; 428 for(CIter c=matrix[i->local()].begin(), cend= matrix[i->local()].end(); 429 c!= cend; ++c) 430 { 431 *c=0; 432 if(c.index()==i->local()) { 433 auto setDiagonal = [](auto&& scalarOrMatrix, const auto& value) { 434 auto&& matrixView = Dune::Impl::asMatrix(scalarOrMatrix); 435 for (auto rowIt = matrixView.begin(); rowIt != matrixView.end(); ++rowIt) 436 (*rowIt)[rowIt.index()] = value; 437 }; 438 setDiagonal(*c, 1); 439 } 440 } 441 } 442 } 443 /** @brief The matrix to communicate the values of. */ 444 M& matrix; 445 /** @brief Index set for the original matrix. */ 446 const Dune::GlobalLookupIndexSet<I>& idxset; 447 /** @brief Index set for the redistributed matrix. */ 448 const I& aggidxset; 449 /** @brief row size information for the receiving side. */ 450 std::vector<size_t>* rowsize; // row sizes differ from sender side in overlap! 451 }; 452 453 template<class M, class I> 454 struct CommPolicy<CommMatrixRow<M,I> > 455 { 456 typedef CommMatrixRow<M,I> Type; 457 458 /** 459 * @brief The indexed type we send. 460 * This is the pair of global index indentitfying the column and the value itself. 461 */ 462 typedef std::pair<typename I::GlobalIndex,typename M::block_type> IndexedType; 463 464 /** @brief Each row varies in size. */ 465 typedef VariableSize IndexedTypeFlag; 466 getSizeDune::CommPolicy467 static std::size_t getSize(const Type& t, std::size_t i) 468 { 469 if(!t.rowsize) 470 return t.matrix[i].size(); 471 else 472 { 473 assert((*t.rowsize)[i]>0); 474 return (*t.rowsize)[i]; 475 } 476 } 477 }; 478 479 template<class M, class I, class RI> 480 struct MatrixRowSizeGatherScatter 481 { 482 typedef CommMatrixRowSize<M,RI> Container; 483 gatherDune::MatrixRowSizeGatherScatter484 static const typename M::size_type gather(const Container& cont, std::size_t i) 485 { 486 return cont.matrix[i].size(); 487 } scatterDune::MatrixRowSizeGatherScatter488 static void scatter(Container& cont, const typename M::size_type& rowsize, std::size_t i) 489 { 490 assert(rowsize); 491 cont.rowsize.getRowSize(i)=rowsize; 492 } 493 494 }; 495 496 template<class M, class I, class RI> 497 struct MatrixCopyRowSizeGatherScatter 498 { 499 typedef CommMatrixRowSize<M,RI> Container; 500 gatherDune::MatrixCopyRowSizeGatherScatter501 static const typename M::size_type gather(const Container& cont, std::size_t i) 502 { 503 return cont.matrix[i].size(); 504 } scatterDune::MatrixCopyRowSizeGatherScatter505 static void scatter(Container& cont, const typename M::size_type& rowsize, std::size_t i) 506 { 507 assert(rowsize); 508 if (rowsize > cont.rowsize.getCopyRowSize(i)) 509 cont.rowsize.getCopyRowSize(i)=rowsize; 510 } 511 512 }; 513 514 template<class M, class I> 515 struct MatrixSparsityPatternGatherScatter 516 { 517 typedef typename I::GlobalIndex GlobalIndex; 518 typedef CommMatrixSparsityPattern<M,I> Container; 519 typedef typename M::ConstColIterator ColIter; 520 521 static ColIter col; 522 static GlobalIndex numlimits; 523 gatherDune::MatrixSparsityPatternGatherScatter524 static const GlobalIndex& gather(const Container& cont, std::size_t i, std::size_t j) 525 { 526 if(j==0) 527 col=cont.matrix[i].begin(); 528 else if (col!=cont.matrix[i].end()) 529 ++col; 530 531 //copy communication: different row sizes for copy rows with the same global index 532 //are possible. If all values of current matrix row are sent, send 533 //std::numeric_limits<GlobalIndex>::max() 534 //and receiver will ignore it. 535 if (col==cont.matrix[i].end()) { 536 numlimits = std::numeric_limits<GlobalIndex>::max(); 537 return numlimits; 538 } 539 else { 540 const typename I::IndexPair* index=cont.idxset.pair(col.index()); 541 assert(index); 542 // Only send index if col is no ghost 543 if ( index->local().attribute() != 2) 544 return index->global(); 545 else { 546 numlimits = std::numeric_limits<GlobalIndex>::max(); 547 return numlimits; 548 } 549 } 550 } scatterDune::MatrixSparsityPatternGatherScatter551 static void scatter(Container& cont, const GlobalIndex& gi, std::size_t i, [[maybe_unused]] std::size_t j) 552 { 553 try{ 554 if (gi != std::numeric_limits<GlobalIndex>::max()) { 555 const typename I::IndexPair& ip=cont.aggidxset.at(gi); 556 assert(ip.global()==gi); 557 std::size_t column = ip.local(); 558 cont.sparsity[i].insert(column); 559 560 typedef typename Dune::OwnerOverlapCopyCommunication<int>::OwnerSet OwnerSet; 561 if(!OwnerSet::contains(ip.local().attribute())) 562 // preserve symmetry for overlap 563 cont.sparsity[column].insert(i); 564 } 565 } 566 catch(const Dune::RangeError&) { 567 // Entry not present in the new index set. Ignore! 568 #ifdef DEBUG_REPART 569 typedef typename Container::LookupIndexSet GlobalLookup; 570 typedef typename GlobalLookup::IndexPair IndexPair; 571 typedef typename Dune::OwnerOverlapCopyCommunication<int>::OwnerSet OwnerSet; 572 573 GlobalLookup lookup(cont.aggidxset); 574 const IndexPair* pi=lookup.pair(i); 575 assert(pi); 576 if(OwnerSet::contains(pi->local().attribute())) { 577 int rank; 578 MPI_Comm_rank(MPI_COMM_WORLD,&rank); 579 std::cout<<rank<<cont.aggidxset<<std::endl; 580 std::cout<<rank<<": row "<<i<<" (global="<<gi <<") not in index set for owner index "<<pi->global()<<std::endl; 581 throw; 582 } 583 #endif 584 } 585 } 586 587 }; 588 template<class M, class I> 589 typename MatrixSparsityPatternGatherScatter<M,I>::ColIter MatrixSparsityPatternGatherScatter<M,I>::col; 590 591 template<class M, class I> 592 typename MatrixSparsityPatternGatherScatter<M,I>::GlobalIndex MatrixSparsityPatternGatherScatter<M,I>::numlimits; 593 594 595 template<class M, class I> 596 struct MatrixRowGatherScatter 597 { 598 typedef typename I::GlobalIndex GlobalIndex; 599 typedef CommMatrixRow<M,I> Container; 600 typedef typename M::ConstColIterator ColIter; 601 typedef typename std::pair<GlobalIndex,typename M::block_type> Data; 602 static ColIter col; 603 static Data datastore; 604 static GlobalIndex numlimits; 605 gatherDune::MatrixRowGatherScatter606 static const Data& gather(const Container& cont, std::size_t i, std::size_t j) 607 { 608 if(j==0) 609 col=cont.matrix[i].begin(); 610 else if (col!=cont.matrix[i].end()) 611 ++col; 612 // copy communication: different row sizes for copy rows with the same global index 613 // are possible. If all values of current matrix row are sent, send 614 // std::numeric_limits<GlobalIndex>::max() 615 // and receiver will ignore it. 616 if (col==cont.matrix[i].end()) { 617 numlimits = std::numeric_limits<GlobalIndex>::max(); 618 datastore = Data(numlimits,*col); 619 return datastore; 620 } 621 else { 622 // convert local column index to global index 623 const typename I::IndexPair* index=cont.idxset.pair(col.index()); 624 assert(index); 625 // Store the data to prevent reference to temporary 626 // Only send index if col is no ghost 627 if ( index->local().attribute() != 2) 628 datastore = Data(index->global(),*col); 629 else { 630 numlimits = std::numeric_limits<GlobalIndex>::max(); 631 datastore = Data(numlimits,*col); 632 } 633 return datastore; 634 } 635 } scatterDune::MatrixRowGatherScatter636 static void scatter(Container& cont, const Data& data, std::size_t i, [[maybe_unused]] std::size_t j) 637 { 638 try{ 639 if (data.first != std::numeric_limits<GlobalIndex>::max()) { 640 typename M::size_type column=cont.aggidxset.at(data.first).local(); 641 cont.matrix[i][column]=data.second; 642 } 643 } 644 catch(const Dune::RangeError&) { 645 // This an overlap row and might therefore lack some entries! 646 } 647 648 } 649 }; 650 651 template<class M, class I> 652 typename MatrixRowGatherScatter<M,I>::ColIter MatrixRowGatherScatter<M,I>::col; 653 654 template<class M, class I> 655 typename MatrixRowGatherScatter<M,I>::Data MatrixRowGatherScatter<M,I>::datastore; 656 657 template<class M, class I> 658 typename MatrixRowGatherScatter<M,I>::GlobalIndex MatrixRowGatherScatter<M,I>::numlimits; 659 660 template<typename M, typename C> redistributeSparsityPattern(M & origMatrix,M & newMatrix,C & origComm,C & newComm,RedistributeInformation<C> & ri)661 void redistributeSparsityPattern(M& origMatrix, M& newMatrix, C& origComm, C& newComm, 662 RedistributeInformation<C>& ri) 663 { 664 typename C::CopySet copyflags; 665 typename C::OwnerSet ownerflags; 666 typedef typename C::ParallelIndexSet IndexSet; 667 typedef RedistributeInformation<C> RI; 668 std::vector<typename M::size_type> rowsize(newComm.indexSet().size(), 0); 669 std::vector<typename M::size_type> copyrowsize(newComm.indexSet().size(), 0); 670 std::vector<typename M::size_type> backwardscopyrowsize(origComm.indexSet().size(), 0); 671 672 // get owner rowsizes 673 CommMatrixRowSize<M,RI> commRowSize(origMatrix, ri); 674 ri.template redistribute<MatrixRowSizeGatherScatter<M,IndexSet,RI> >(commRowSize,commRowSize); 675 676 origComm.buildGlobalLookup(); 677 678 for (std::size_t i=0; i < newComm.indexSet().size(); i++) { 679 rowsize[i] = ri.getRowSize(i); 680 } 681 // get sparsity pattern from owner rows 682 CommMatrixSparsityPattern<M,IndexSet> 683 origsp(origMatrix, origComm.globalLookup(), newComm.indexSet()); 684 CommMatrixSparsityPattern<M,IndexSet> 685 newsp(origMatrix, origComm.globalLookup(), newComm.indexSet(), rowsize); 686 687 ri.template redistribute<MatrixSparsityPatternGatherScatter<M,IndexSet> >(origsp,newsp); 688 689 // build copy to owner interface to get missing matrix values for novlp case 690 if (SolverCategory::category(origComm) == SolverCategory::nonoverlapping) { 691 RemoteIndices<IndexSet> *ris = new RemoteIndices<IndexSet>(origComm.indexSet(), 692 newComm.indexSet(), 693 origComm.communicator()); 694 ris->template rebuild<true>(); 695 696 ri.getInterface().free(); 697 ri.getInterface().build(*ris,copyflags,ownerflags); 698 699 // get copy rowsizes 700 CommMatrixRowSize<M,RI> commRowSize_copy(origMatrix, ri); 701 ri.template redistribute<MatrixCopyRowSizeGatherScatter<M,IndexSet,RI> >(commRowSize_copy, 702 commRowSize_copy); 703 704 for (std::size_t i=0; i < newComm.indexSet().size(); i++) { 705 copyrowsize[i] = ri.getCopyRowSize(i); 706 } 707 //get copy rowsizes for sender 708 ri.redistributeBackward(backwardscopyrowsize,copyrowsize); 709 for (std::size_t i=0; i < origComm.indexSet().size(); i++) { 710 ri.getBackwardsCopyRowSize(i) = backwardscopyrowsize[i]; 711 } 712 713 // get sparsity pattern from copy rows 714 CommMatrixSparsityPattern<M,IndexSet> origsp_copy(origMatrix, 715 origComm.globalLookup(), 716 newComm.indexSet(), 717 backwardscopyrowsize); 718 CommMatrixSparsityPattern<M,IndexSet> newsp_copy(origMatrix, origComm.globalLookup(), 719 newComm.indexSet(), copyrowsize); 720 ri.template redistribute<MatrixSparsityPatternGatherScatter<M,IndexSet> >(origsp_copy, 721 newsp_copy); 722 723 newsp.completeSparsityPattern(newsp_copy.sparsity); 724 newsp.storeSparsityPattern(newMatrix); 725 } 726 else 727 newsp.storeSparsityPattern(newMatrix); 728 729 #ifdef DUNE_ISTL_WITH_CHECKING 730 // Check for symmetry 731 int ret=0; 732 typedef typename M::ConstRowIterator RIter; 733 for(RIter row=newMatrix.begin(), rend=newMatrix.end(); row != rend; ++row) { 734 typedef typename M::ConstColIterator CIter; 735 for(CIter col=row->begin(), cend=row->end(); col!=cend; ++col) 736 { 737 try{ 738 newMatrix[col.index()][row.index()]; 739 }catch(const Dune::ISTLError&) { 740 std::cerr<<newComm.communicator().rank()<<": entry (" 741 <<col.index()<<","<<row.index()<<") missing! for symmetry!"<<std::endl; 742 ret=1; 743 744 } 745 746 } 747 } 748 749 if(ret) 750 DUNE_THROW(ISTLError, "Matrix not symmetric!"); 751 #endif 752 } 753 754 template<typename M, typename C> redistributeMatrixEntries(M & origMatrix,M & newMatrix,C & origComm,C & newComm,RedistributeInformation<C> & ri)755 void redistributeMatrixEntries(M& origMatrix, M& newMatrix, C& origComm, C& newComm, 756 RedistributeInformation<C>& ri) 757 { 758 typedef typename C::ParallelIndexSet IndexSet; 759 typename C::OwnerSet ownerflags; 760 std::vector<typename M::size_type> rowsize(newComm.indexSet().size(), 0); 761 std::vector<typename M::size_type> copyrowsize(newComm.indexSet().size(), 0); 762 std::vector<typename M::size_type> backwardscopyrowsize(origComm.indexSet().size(), 0); 763 764 for (std::size_t i=0; i < newComm.indexSet().size(); i++) { 765 rowsize[i] = ri.getRowSize(i); 766 if (SolverCategory::category(origComm) == SolverCategory::nonoverlapping) { 767 copyrowsize[i] = ri.getCopyRowSize(i); 768 } 769 } 770 771 for (std::size_t i=0; i < origComm.indexSet().size(); i++) 772 if (SolverCategory::category(origComm) == SolverCategory::nonoverlapping) 773 backwardscopyrowsize[i] = ri.getBackwardsCopyRowSize(i); 774 775 776 if (SolverCategory::category(origComm) == SolverCategory::nonoverlapping) { 777 // fill sparsity pattern from copy rows 778 CommMatrixRow<M,IndexSet> origrow_copy(origMatrix, origComm.globalLookup(), 779 newComm.indexSet(), backwardscopyrowsize); 780 CommMatrixRow<M,IndexSet> newrow_copy(newMatrix, origComm.globalLookup(), 781 newComm.indexSet(),copyrowsize); 782 ri.template redistribute<MatrixRowGatherScatter<M,IndexSet> >(origrow_copy, 783 newrow_copy); 784 ri.getInterface().free(); 785 RemoteIndices<IndexSet> *ris = new RemoteIndices<IndexSet>(origComm.indexSet(), 786 newComm.indexSet(), 787 origComm.communicator()); 788 ris->template rebuild<true>(); 789 ri.getInterface().build(*ris,ownerflags,ownerflags); 790 } 791 792 CommMatrixRow<M,IndexSet> 793 origrow(origMatrix, origComm.globalLookup(), newComm.indexSet()); 794 CommMatrixRow<M,IndexSet> 795 newrow(newMatrix, origComm.globalLookup(), newComm.indexSet(),rowsize); 796 ri.template redistribute<MatrixRowGatherScatter<M,IndexSet> >(origrow,newrow); 797 if (SolverCategory::category(origComm) != static_cast<int>(SolverCategory::nonoverlapping)) 798 newrow.setOverlapRowsToDirichlet(); 799 } 800 801 /** 802 * @brief Redistribute a matrix according to given domain decompositions. 803 * 804 * All the parameters for this function can be obtained by calling 805 * graphRepartition with the graph of the original matrix. 806 * 807 * @param origMatrix The matrix on the original partitioning. 808 * @param newMatrix An empty matrix to store the new redistributed matrix in. 809 * @param origComm The parallel information of the original partitioning. 810 * @param newComm The parallel information of the new partitioning. 811 * @param ri The remote index information between the original and the new partitioning. 812 * Upon exit of this method it will be prepared for copying from owner to owner vertices 813 * for data redistribution. 814 * @tparam M The matrix type. It is assumed to be sparse. E.g. BCRSMatrix. 815 * @tparam C The type of the parallel information, see OwnerOverlapCopyCommunication. 816 */ 817 template<typename M, typename C> redistributeMatrix(M & origMatrix,M & newMatrix,C & origComm,C & newComm,RedistributeInformation<C> & ri)818 void redistributeMatrix(M& origMatrix, M& newMatrix, C& origComm, C& newComm, 819 RedistributeInformation<C>& ri) 820 { 821 ri.setNoRows(newComm.indexSet().size()); 822 ri.setNoCopyRows(newComm.indexSet().size()); 823 ri.setNoBackwardsCopyRows(origComm.indexSet().size()); 824 redistributeSparsityPattern(origMatrix, newMatrix, origComm, newComm, ri); 825 redistributeMatrixEntries(origMatrix, newMatrix, origComm, newComm, ri); 826 } 827 #endif 828 829 template<typename M> redistributeMatrixEntries(M & origMatrix,M & newMatrix,Dune::Amg::SequentialInformation & origComm,Dune::Amg::SequentialInformation & newComm,RedistributeInformation<Dune::Amg::SequentialInformation> & ri)830 void redistributeMatrixEntries(M& origMatrix, M& newMatrix, 831 Dune::Amg::SequentialInformation& origComm, 832 Dune::Amg::SequentialInformation& newComm, 833 RedistributeInformation<Dune::Amg::SequentialInformation>& ri) 834 { 835 DUNE_THROW(InvalidStateException, "Trying to redistribute in sequential program!"); 836 } 837 template<typename M> redistributeMatrix(M & origMatrix,M & newMatrix,Dune::Amg::SequentialInformation & origComm,Dune::Amg::SequentialInformation & newComm,RedistributeInformation<Dune::Amg::SequentialInformation> & ri)838 void redistributeMatrix(M& origMatrix, M& newMatrix, 839 Dune::Amg::SequentialInformation& origComm, 840 Dune::Amg::SequentialInformation& newComm, 841 RedistributeInformation<Dune::Amg::SequentialInformation>& ri) 842 { 843 DUNE_THROW(InvalidStateException, "Trying to redistribute in sequential program!"); 844 } 845 } 846 #endif 847