3 #ifndef DUNE_ISTL_MATRIXREDISTRIBUTE_HH
4 #define DUNE_ISTL_MATRIXREDISTRIBUTE_HH
7 #include <dune/common/exceptions.hh>
8 #include <dune/common/parallel/indexset.hh>
9 #include <dune/common/unused.hh>
29 DUNE_UNUSED_PARAMETER(from);
30 DUNE_UNUSED_PARAMETER(to);
36 DUNE_UNUSED_PARAMETER(from);
37 DUNE_UNUSED_PARAMETER(to);
45 DUNE_UNUSED_PARAMETER(size);
50 DUNE_UNUSED_PARAMETER(size);
55 DUNE_UNUSED_PARAMETER(size);
60 DUNE_UNUSED_PARAMETER(index);
66 DUNE_UNUSED_PARAMETER(index);
72 DUNE_UNUSED_PARAMETER(index);
79 template<
typename T,
typename T1>
86 : interface(), setup_(false)
95 const IS& target, MPI_Comm comm)
97 auto ri = std::make_unique<RemoteIndices<IS> >(source, target, comm);
98 ri->template rebuild<true>();
102 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
104 inf.build(*ri, flags, flags);
110 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
112 std::cout<<
"Interfaces do not match!"<<std::endl;
113 std::cout<<rank<<
": redist interface new :"<<inf<<std::endl;
114 std::cout<<rank<<
": redist interface :"<<interface<<std::endl;
131 template<
class GatherScatter,
class D>
134 BufferedCommunicator communicator;
135 communicator.template build<D>(from,to, interface);
136 communicator.template forward<GatherScatter>(from, to);
139 template<
class GatherScatter,
class D>
143 BufferedCommunicator communicator;
144 communicator.template build<D>(from,to, interface);
145 communicator.template backward<GatherScatter>(from, to);
152 redistribute<CopyGatherScatter<D> >(from,to);
157 redistributeBackward<CopyGatherScatter<D> >(from,to);
169 return rowSize[index];
174 return rowSize[index];
179 return copyrowSize[index];
184 return copyrowSize[index];
189 return backwardscopyrowSize[index];
194 return backwardscopyrowSize[index];
199 rowSize.resize(rows, 0);
204 copyrowSize.resize(rows, 0);
209 backwardscopyrowSize.resize(rows, 0);
213 std::vector<std::size_t> rowSize;
214 std::vector<std::size_t> copyrowSize;
215 std::vector<std::size_t> backwardscopyrowSize;
228 template<
class M,
class RI>
257 template<
class M,
class I>
280 const std::vector<typename M::size_type>& rowsize_)
293 typedef typename Dune::GlobalLookupIndexSet<I>::const_iterator IIter;
299 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
302 if(!OwnerSet::contains(i->local().attribute())) {
304 std::cout<<rank<<
" Inserting diagonal for"<<i->local()<<std::endl;
306 sparsity[i->local()].insert(i->local());
315 m.setBuildMode(M::row_wise);
316 typename M::CreateIterator citer=m.createbegin();
320 Dune::GlobalLookupIndexSet<I> global(
aggidxset);
322 typedef typename std::vector<std::set<size_type> >::const_iterator Iter;
325 typedef typename std::set<size_type>::const_iterator SIter;
326 for(SIter si=i->begin(), send=i->end(); si!=send; ++si)
329 if(i->find(idx)==i->end()) {
330 const typename I::IndexPair* gi=global.pair(idx);
332 std::cout<<rank<<
": row "<<idx<<
" is missing a diagonal entry! global="<<gi->global()<<
" attr="<<gi->local().attribute()<<
" "<<
333 OwnerSet::contains(gi->local().attribute())<<
334 " row size="<<i->size()<<std::endl;
356 for (
unsigned int i = 0; i !=
sparsity.size(); ++i) {
357 if (add_sparsity[i].size() != 0) {
358 typedef std::set<size_type> Set;
360 std::insert_iterator<Set> tmp_insert (tmp_set, tmp_set.begin());
361 std::set_union(add_sparsity[i].begin(), add_sparsity[i].end(),
370 const Dune::GlobalLookupIndexSet<I>&
idxset;
376 template<
class M,
class I>
390 static typename M::size_type
getSize(
const Type& t, std::size_t i)
393 return t.
matrix[i].size();
408 template<
class M,
class I>
419 CommMatrixRow(M& m_,
const Dune::GlobalLookupIndexSet<I>& idxset_,
const I& aggidxset_)
426 CommMatrixRow(M& m_,
const Dune::GlobalLookupIndexSet<I>& idxset_,
const I& aggidxset_,
427 std::vector<size_t>& rowsize_)
437 typedef typename Dune::GlobalLookupIndexSet<I>::const_iterator Iter;
441 if(!OwnerSet::contains(i->local().attribute())) {
443 typedef typename M::ColIterator CIter;
444 for(CIter c=
matrix[i->local()].begin(), cend=
matrix[i->local()].end();
448 if(c.index()==i->local()) {
449 auto setDiagonal = [](
auto&& scalarOrMatrix,
const auto& value) {
450 auto&&
matrix = Dune::Impl::asMatrix(scalarOrMatrix);
451 for (
auto rowIt =
matrix.begin(); rowIt !=
matrix.end(); ++rowIt)
452 (*rowIt)[rowIt.index()] = value;
462 const Dune::GlobalLookupIndexSet<I>&
idxset;
469 template<
class M,
class I>
478 typedef std::pair<typename I::GlobalIndex,typename M::block_type>
IndexedType;
486 return t.
matrix[i].size();
495 template<
class M,
class I,
class RI>
502 return cont.
matrix[i].size();
507 cont.
rowsize.getRowSize(i)=rowsize;
512 template<
class M,
class I,
class RI>
519 return cont.
matrix[i].size();
524 if (rowsize > cont.
rowsize.getCopyRowSize(i))
525 cont.
rowsize.getCopyRowSize(i)=rowsize;
530 template<
class M,
class I>
552 numlimits = std::numeric_limits<GlobalIndex>::max();
556 const typename I::IndexPair* index=cont.
idxset.pair(
col.index());
559 if ( index->local().attribute() != 2)
560 return index->global();
562 numlimits = std::numeric_limits<GlobalIndex>::max();
569 DUNE_UNUSED_PARAMETER(j);
571 if (gi != std::numeric_limits<GlobalIndex>::max()) {
572 const typename I::IndexPair& ip=cont.
aggidxset.at(gi);
573 assert(ip.global()==gi);
574 std::size_t column = ip.local();
578 if(!OwnerSet::contains(ip.local().attribute()))
583 catch(
const Dune::RangeError&) {
587 typedef typename GlobalLookup::IndexPair IndexPair;
591 const IndexPair* pi=lookup.pair(i);
593 if(OwnerSet::contains(pi->local().attribute())) {
595 MPI_Comm_rank(MPI_COMM_WORLD,&rank);
596 std::cout<<rank<<cont.
aggidxset<<std::endl;
597 std::cout<<rank<<
": row "<<i<<
" (global="<<gi <<
") not in index set for owner index "<<pi->global()<<std::endl;
605 template<
class M,
class I>
608 template<
class M,
class I>
612 template<
class M,
class I>
618 typedef typename std::pair<GlobalIndex,typename M::block_type>
Data;
634 numlimits = std::numeric_limits<GlobalIndex>::max();
640 const typename I::IndexPair* index=cont.
idxset.pair(
col.index());
644 if ( index->local().attribute() != 2)
647 numlimits = std::numeric_limits<GlobalIndex>::max();
655 DUNE_UNUSED_PARAMETER(j);
657 if (data.first != std::numeric_limits<GlobalIndex>::max()) {
658 typename M::size_type column=cont.
aggidxset.at(data.first).local();
659 cont.
matrix[i][column]=data.second;
662 catch(
const Dune::RangeError&) {
669 template<
class M,
class I>
672 template<
class M,
class I>
675 template<
class M,
class I>
678 template<
typename M,
typename C>
682 typename C::CopySet copyflags;
683 typename C::OwnerSet ownerflags;
684 typedef typename C::ParallelIndexSet IndexSet;
686 std::vector<typename M::size_type> rowsize(newComm.indexSet().size(), 0);
687 std::vector<typename M::size_type> copyrowsize(newComm.indexSet().size(), 0);
688 std::vector<typename M::size_type> backwardscopyrowsize(origComm.indexSet().size(), 0);
692 ri.template redistribute<MatrixRowSizeGatherScatter<M,IndexSet,RI> >(commRowSize,commRowSize);
694 origComm.buildGlobalLookup();
696 for (std::size_t i=0; i < newComm.indexSet().size(); i++) {
701 origsp(origMatrix, origComm.globalLookup(), newComm.indexSet());
703 newsp(origMatrix, origComm.globalLookup(), newComm.indexSet(), rowsize);
705 ri.template redistribute<MatrixSparsityPatternGatherScatter<M,IndexSet> >(origsp,newsp);
709 RemoteIndices<IndexSet> *ris =
new RemoteIndices<IndexSet>(origComm.indexSet(),
711 origComm.communicator());
712 ris->template rebuild<true>();
714 ri.getInterface().free();
715 ri.getInterface().build(*ris,copyflags,ownerflags);
719 ri.template redistribute<MatrixCopyRowSizeGatherScatter<M,IndexSet,RI> >(commRowSize_copy,
722 for (std::size_t i=0; i < newComm.indexSet().size(); i++) {
727 for (std::size_t i=0; i < origComm.indexSet().size(); i++) {
733 origComm.globalLookup(),
735 backwardscopyrowsize);
737 newComm.indexSet(), copyrowsize);
738 ri.template redistribute<MatrixSparsityPatternGatherScatter<M,IndexSet> >(origsp_copy,
747 #ifdef DUNE_ISTL_WITH_CHECKING
750 typedef typename M::ConstRowIterator RIter;
751 for(RIter row=newMatrix.begin(), rend=newMatrix.end(); row != rend; ++row) {
752 typedef typename M::ConstColIterator CIter;
753 for(CIter
col=row->begin(), cend=row->end();
col!=cend; ++
col)
756 newMatrix[
col.index()][row.index()];
758 std::cerr<<newComm.communicator().rank()<<
": entry ("
759 <<
col.index()<<
","<<row.index()<<
") missing! for symmetry!"<<std::endl;
768 DUNE_THROW(
ISTLError,
"Matrix not symmetric!");
772 template<
typename M,
typename C>
776 typedef typename C::ParallelIndexSet IndexSet;
777 typename C::OwnerSet ownerflags;
778 std::vector<typename M::size_type> rowsize(newComm.indexSet().size(), 0);
779 std::vector<typename M::size_type> copyrowsize(newComm.indexSet().size(), 0);
780 std::vector<typename M::size_type> backwardscopyrowsize(origComm.indexSet().size(), 0);
782 for (std::size_t i=0; i < newComm.indexSet().size(); i++) {
789 for (std::size_t i=0; i < origComm.indexSet().size(); i++)
797 newComm.indexSet(), backwardscopyrowsize);
799 newComm.indexSet(),copyrowsize);
800 ri.template redistribute<MatrixRowGatherScatter<M,IndexSet> >(origrow_copy,
802 ri.getInterface().free();
803 RemoteIndices<IndexSet> *ris =
new RemoteIndices<IndexSet>(origComm.indexSet(),
805 origComm.communicator());
806 ris->template rebuild<true>();
807 ri.getInterface().build(*ris,ownerflags,ownerflags);
811 origrow(origMatrix, origComm.globalLookup(), newComm.indexSet());
813 newrow(newMatrix, origComm.globalLookup(), newComm.indexSet(),rowsize);
814 ri.template redistribute<MatrixRowGatherScatter<M,IndexSet> >(origrow,newrow);
835 template<
typename M,
typename C>
853 DUNE_THROW(InvalidStateException,
"Trying to redistribute in sequential program!");
861 DUNE_THROW(InvalidStateException,
"Trying to redistribute in sequential program!");
Classes providing communication interfaces for overlapping Schwarz methods.
Functionality for redistributing a parallel index set using graph partitioning.
Col col
Definition: matrixmatrix.hh:349
Definition: allocator.hh:7
void redistributeMatrixEntries(M &origMatrix, M &newMatrix, C &origComm, C &newComm, RedistributeInformation< C > &ri)
Definition: matrixredistribute.hh:773
void redistributeSparsityPattern(M &origMatrix, M &newMatrix, C &origComm, C &newComm, RedistributeInformation< C > &ri)
Definition: matrixredistribute.hh:679
void redistributeMatrix(M &origMatrix, M &newMatrix, C &origComm, C &newComm, RedistributeInformation< C > &ri)
Redistribute a matrix according to given domain decompositions.
Definition: matrixredistribute.hh:836
derive error class from the base class in common
Definition: istlexception.hh:16
Definition: matrixredistribute.hh:21
void setNoBackwardsCopyRows(std::size_t size)
Definition: matrixredistribute.hh:53
void redistribute(const D &from, D &to) const
Definition: matrixredistribute.hh:27
void resetSetup()
Definition: matrixredistribute.hh:40
void setNoCopyRows(std::size_t size)
Definition: matrixredistribute.hh:48
bool isSetup() const
Definition: matrixredistribute.hh:22
void setNoRows(std::size_t size)
Definition: matrixredistribute.hh:43
void redistributeBackward(D &from, const D &to) const
Definition: matrixredistribute.hh:34
std::size_t getBackwardsCopyRowSize(std::size_t index) const
Definition: matrixredistribute.hh:70
std::size_t getRowSize(std::size_t index) const
Definition: matrixredistribute.hh:58
std::size_t getCopyRowSize(std::size_t index) const
Definition: matrixredistribute.hh:64
std::size_t getRowSize(std::size_t index) const
Definition: matrixredistribute.hh:172
void redistribute(const D &from, D &to) const
Definition: matrixredistribute.hh:150
void setNoBackwardsCopyRows(std::size_t rows)
Definition: matrixredistribute.hh:207
std::size_t & getBackwardsCopyRowSize(std::size_t index)
Definition: matrixredistribute.hh:187
RedistributeInformation()
Definition: matrixredistribute.hh:85
std::size_t getCopyRowSize(std::size_t index) const
Definition: matrixredistribute.hh:182
void setNoRows(std::size_t rows)
Definition: matrixredistribute.hh:197
RedistributeInterface & getInterface()
Definition: matrixredistribute.hh:89
void reserve(std::size_t size)
Definition: matrixredistribute.hh:164
OwnerOverlapCopyCommunication< T, T1 > Comm
Definition: matrixredistribute.hh:83
void setNoCopyRows(std::size_t rows)
Definition: matrixredistribute.hh:202
void redistributeBackward(D &from, const D &to) const
Definition: matrixredistribute.hh:155
void setSetup()
Definition: matrixredistribute.hh:120
std::size_t getBackwardsCopyRowSize(std::size_t index) const
Definition: matrixredistribute.hh:192
void redistribute(const D &from, D &to) const
Definition: matrixredistribute.hh:132
void resetSetup()
Definition: matrixredistribute.hh:126
void redistributeBackward(D &from, const D &to) const
Definition: matrixredistribute.hh:140
std::size_t & getCopyRowSize(std::size_t index)
Definition: matrixredistribute.hh:177
void checkInterface(const IS &source, const IS &target, MPI_Comm comm)
Definition: matrixredistribute.hh:94
bool isSetup() const
Definition: matrixredistribute.hh:159
std::size_t & getRowSize(std::size_t index)
Definition: matrixredistribute.hh:167
Utility class to communicate and set the row sizes of a redistributed matrix.
Definition: matrixredistribute.hh:230
M::size_type size_type
Definition: matrixredistribute.hh:233
M::size_type value_type
Definition: matrixredistribute.hh:232
RI & rowsize
Definition: matrixredistribute.hh:244
const M & matrix
Definition: matrixredistribute.hh:243
CommMatrixRowSize(const M &m_, RI &rowsize_)
Constructor.
Definition: matrixredistribute.hh:240
Utility class to communicate and build the sparsity pattern of a redistributed matrix.
Definition: matrixredistribute.hh:259
M::size_type size_type
Definition: matrixredistribute.hh:260
const Dune::GlobalLookupIndexSet< I > & idxset
Definition: matrixredistribute.hh:370
void storeSparsityPattern(M &m)
Creates and stores the sparsity pattern of the redistributed matrix.
Definition: matrixredistribute.hh:290
const I & aggidxset
Definition: matrixredistribute.hh:371
const std::vector< size_type > * rowsize
Definition: matrixredistribute.hh:373
void completeSparsityPattern(std::vector< std::set< size_type > > add_sparsity)
Completes the sparsity pattern of the redistributed matrix with data from copy rows for the novlp cas...
Definition: matrixredistribute.hh:354
CommMatrixSparsityPattern(const M &m_, const Dune::GlobalLookupIndexSet< I > &idxset_, const I &aggidxset_)
Constructor for the original side.
Definition: matrixredistribute.hh:268
const M & matrix
Definition: matrixredistribute.hh:368
CommMatrixSparsityPattern(const M &m_, const Dune::GlobalLookupIndexSet< I > &idxset_, const I &aggidxset_, const std::vector< typename M::size_type > &rowsize_)
Constructor for the redistruted side.
Definition: matrixredistribute.hh:279
std::vector< std::set< size_type > > sparsity
Definition: matrixredistribute.hh:372
Dune::GlobalLookupIndexSet< I > LookupIndexSet
Definition: matrixredistribute.hh:369
static M::size_type getSize(const Type &t, std::size_t i)
Definition: matrixredistribute.hh:390
CommMatrixSparsityPattern< M, I > Type
Definition: matrixredistribute.hh:379
I::GlobalIndex IndexedType
The indexed type we send. This is the global index indentitfying the column.
Definition: matrixredistribute.hh:385
VariableSize IndexedTypeFlag
Each row varies in size.
Definition: matrixredistribute.hh:388
Utility class for comunicating the matrix entries.
Definition: matrixredistribute.hh:410
std::vector< size_t > * rowsize
row size information for the receiving side.
Definition: matrixredistribute.hh:466
M & matrix
The matrix to communicate the values of.
Definition: matrixredistribute.hh:460
CommMatrixRow(M &m_, const Dune::GlobalLookupIndexSet< I > &idxset_, const I &aggidxset_, std::vector< size_t > &rowsize_)
Constructor.
Definition: matrixredistribute.hh:426
const Dune::GlobalLookupIndexSet< I > & idxset
Index set for the original matrix.
Definition: matrixredistribute.hh:462
void setOverlapRowsToDirichlet()
Sets the non-owner rows correctly as Dirichlet boundaries.
Definition: matrixredistribute.hh:435
const I & aggidxset
Index set for the redistributed matrix.
Definition: matrixredistribute.hh:464
CommMatrixRow(M &m_, const Dune::GlobalLookupIndexSet< I > &idxset_, const I &aggidxset_)
Constructor.
Definition: matrixredistribute.hh:419
std::pair< typename I::GlobalIndex, typename M::block_type > IndexedType
The indexed type we send. This is the pair of global index indentitfying the column and the value its...
Definition: matrixredistribute.hh:478
CommMatrixRow< M, I > Type
Definition: matrixredistribute.hh:472
static std::size_t getSize(const Type &t, std::size_t i)
Definition: matrixredistribute.hh:483
VariableSize IndexedTypeFlag
Each row varies in size.
Definition: matrixredistribute.hh:481
Definition: matrixredistribute.hh:497
static void scatter(Container &cont, const typename M::size_type &rowsize, std::size_t i)
Definition: matrixredistribute.hh:504
static const M::size_type gather(const Container &cont, std::size_t i)
Definition: matrixredistribute.hh:500
CommMatrixRowSize< M, RI > Container
Definition: matrixredistribute.hh:498
Definition: matrixredistribute.hh:514
static const M::size_type gather(const Container &cont, std::size_t i)
Definition: matrixredistribute.hh:517
static void scatter(Container &cont, const typename M::size_type &rowsize, std::size_t i)
Definition: matrixredistribute.hh:521
CommMatrixRowSize< M, RI > Container
Definition: matrixredistribute.hh:515
Definition: matrixredistribute.hh:532
M::ConstColIterator ColIter
Definition: matrixredistribute.hh:535
static void scatter(Container &cont, const GlobalIndex &gi, std::size_t i, std::size_t j)
Definition: matrixredistribute.hh:567
CommMatrixSparsityPattern< M, I > Container
Definition: matrixredistribute.hh:534
static GlobalIndex numlimits
Definition: matrixredistribute.hh:538
static ColIter col
Definition: matrixredistribute.hh:537
I::GlobalIndex GlobalIndex
Definition: matrixredistribute.hh:533
static const GlobalIndex & gather(const Container &cont, std::size_t i, std::size_t j)
Definition: matrixredistribute.hh:540
Definition: matrixredistribute.hh:614
I::GlobalIndex GlobalIndex
Definition: matrixredistribute.hh:615
static Data datastore
Definition: matrixredistribute.hh:620
static GlobalIndex numlimits
Definition: matrixredistribute.hh:621
static const Data & gather(const Container &cont, std::size_t i, std::size_t j)
Definition: matrixredistribute.hh:623
M::ConstColIterator ColIter
Definition: matrixredistribute.hh:617
std::pair< GlobalIndex, typename M::block_type > Data
Definition: matrixredistribute.hh:618
static void scatter(Container &cont, const Data &data, std::size_t i, std::size_t j)
Definition: matrixredistribute.hh:653
static ColIter col
Definition: matrixredistribute.hh:619
CommMatrixRow< M, I > Container
Definition: matrixredistribute.hh:616
A class setting up standard communication for a two-valued attribute set with owner/overlap/copy sema...
Definition: owneroverlapcopy.hh:172
EnumItem< AttributeSet, OwnerOverlapCopyAttributeSet::owner > OwnerSet
Definition: owneroverlapcopy.hh:192
Definition: repartition.hh:268
@ nonoverlapping
Category for non-overlapping solvers.
Definition: solvercategory.hh:25
static Category category(const OP &op, decltype(op.category()) *=nullptr)
Helperfunction to extract the solver category either from an enum, or from the newly introduced virtu...
Definition: solvercategory.hh:32