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>
26 void redistribute([[maybe_unused]]
const D& from, [[maybe_unused]] D& to)
const
45 std::size_t
getRowSize([[maybe_unused]] std::size_t index)
const
63 template<
typename T,
typename T1>
70 : interface(), setup_(false)
79 const IS& target, MPI_Comm comm)
81 auto ri = std::make_unique<RemoteIndices<IS> >(source, target, comm);
82 ri->template rebuild<true>();
86 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
88 inf.build(*ri, flags, flags);
94 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
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;
115 template<
class GatherScatter,
class D>
118 BufferedCommunicator communicator;
119 communicator.template build<D>(from,to, interface);
120 communicator.template forward<GatherScatter>(from, to);
123 template<
class GatherScatter,
class D>
127 BufferedCommunicator communicator;
128 communicator.template build<D>(from,to, interface);
129 communicator.template backward<GatherScatter>(from, to);
136 redistribute<CopyGatherScatter<D> >(from,to);
141 redistributeBackward<CopyGatherScatter<D> >(from,to);
153 return rowSize[index];
158 return rowSize[index];
163 return copyrowSize[index];
168 return copyrowSize[index];
173 return backwardscopyrowSize[index];
178 return backwardscopyrowSize[index];
183 rowSize.resize(rows, 0);
188 copyrowSize.resize(rows, 0);
193 backwardscopyrowSize.resize(rows, 0);
197 std::vector<std::size_t> rowSize;
198 std::vector<std::size_t> copyrowSize;
199 std::vector<std::size_t> backwardscopyrowSize;
212 template<
class M,
class RI>
241 template<
class M,
class I>
264 const std::vector<typename M::size_type>& rowsize_)
277 typedef typename Dune::GlobalLookupIndexSet<I>::const_iterator IIter;
283 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
286 if(!OwnerSet::contains(i->local().attribute())) {
288 std::cout<<rank<<
" Inserting diagonal for"<<i->local()<<std::endl;
290 sparsity[i->local()].insert(i->local());
299 m.setBuildMode(M::row_wise);
300 typename M::CreateIterator citer=m.createbegin();
304 Dune::GlobalLookupIndexSet<I> global(
aggidxset);
306 typedef typename std::vector<std::set<size_type> >::const_iterator Iter;
309 typedef typename std::set<size_type>::const_iterator SIter;
310 for(SIter si=i->begin(), send=i->end(); si!=send; ++si)
313 if(i->find(idx)==i->end()) {
314 const typename I::IndexPair* gi=global.pair(idx);
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;
340 for (
unsigned int i = 0; i !=
sparsity.size(); ++i) {
341 if (add_sparsity[i].size() != 0) {
342 typedef std::set<size_type> 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(),
354 const Dune::GlobalLookupIndexSet<I>&
idxset;
360 template<
class M,
class I>
374 static typename M::size_type
getSize(
const Type& t, std::size_t i)
377 return t.
matrix[i].size();
392 template<
class M,
class I>
403 CommMatrixRow(M& m_,
const Dune::GlobalLookupIndexSet<I>& idxset_,
const I& aggidxset_)
410 CommMatrixRow(M& m_,
const Dune::GlobalLookupIndexSet<I>& idxset_,
const I& aggidxset_,
411 std::vector<size_t>& rowsize_)
421 typedef typename Dune::GlobalLookupIndexSet<I>::const_iterator Iter;
425 if(!OwnerSet::contains(i->local().attribute())) {
427 typedef typename M::ColIterator CIter;
428 for(CIter c=
matrix[i->local()].begin(), cend=
matrix[i->local()].end();
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;
446 const Dune::GlobalLookupIndexSet<I>&
idxset;
453 template<
class M,
class I>
462 typedef std::pair<typename I::GlobalIndex,typename M::block_type>
IndexedType;
470 return t.
matrix[i].size();
479 template<
class M,
class I,
class RI>
486 return cont.
matrix[i].size();
491 cont.
rowsize.getRowSize(i)=rowsize;
496 template<
class M,
class I,
class RI>
503 return cont.
matrix[i].size();
508 if (rowsize > cont.
rowsize.getCopyRowSize(i))
509 cont.
rowsize.getCopyRowSize(i)=rowsize;
514 template<
class M,
class I>
536 numlimits = std::numeric_limits<GlobalIndex>::max();
540 const typename I::IndexPair* index=cont.
idxset.pair(
col.index());
543 if ( index->local().attribute() != 2)
544 return index->global();
546 numlimits = std::numeric_limits<GlobalIndex>::max();
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();
561 if(!OwnerSet::contains(ip.local().attribute()))
566 catch(
const Dune::RangeError&) {
570 typedef typename GlobalLookup::IndexPair IndexPair;
574 const IndexPair* pi=lookup.pair(i);
576 if(OwnerSet::contains(pi->local().attribute())) {
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;
588 template<
class M,
class I>
591 template<
class M,
class I>
595 template<
class M,
class I>
601 typedef typename std::pair<GlobalIndex,typename M::block_type>
Data;
617 numlimits = std::numeric_limits<GlobalIndex>::max();
623 const typename I::IndexPair* index=cont.
idxset.pair(
col.index());
627 if ( index->local().attribute() != 2)
630 numlimits = std::numeric_limits<GlobalIndex>::max();
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;
644 catch(
const Dune::RangeError&) {
651 template<
class M,
class I>
654 template<
class M,
class I>
657 template<
class M,
class I>
660 template<
typename M,
typename C>
664 typename C::CopySet copyflags;
665 typename C::OwnerSet ownerflags;
666 typedef typename C::ParallelIndexSet IndexSet;
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);
674 ri.template redistribute<MatrixRowSizeGatherScatter<M,IndexSet,RI> >(commRowSize,commRowSize);
676 origComm.buildGlobalLookup();
678 for (std::size_t i=0; i < newComm.indexSet().size(); i++) {
683 origsp(origMatrix, origComm.globalLookup(), newComm.indexSet());
685 newsp(origMatrix, origComm.globalLookup(), newComm.indexSet(), rowsize);
687 ri.template redistribute<MatrixSparsityPatternGatherScatter<M,IndexSet> >(origsp,newsp);
691 RemoteIndices<IndexSet> *ris =
new RemoteIndices<IndexSet>(origComm.indexSet(),
693 origComm.communicator());
694 ris->template rebuild<true>();
696 ri.getInterface().free();
697 ri.getInterface().build(*ris,copyflags,ownerflags);
701 ri.template redistribute<MatrixCopyRowSizeGatherScatter<M,IndexSet,RI> >(commRowSize_copy,
704 for (std::size_t i=0; i < newComm.indexSet().size(); i++) {
709 for (std::size_t i=0; i < origComm.indexSet().size(); i++) {
715 origComm.globalLookup(),
717 backwardscopyrowsize);
719 newComm.indexSet(), copyrowsize);
720 ri.template redistribute<MatrixSparsityPatternGatherScatter<M,IndexSet> >(origsp_copy,
729#ifdef DUNE_ISTL_WITH_CHECKING
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)
738 newMatrix[
col.index()][row.index()];
740 std::cerr<<newComm.communicator().rank()<<
": entry ("
741 <<
col.index()<<
","<<row.index()<<
") missing! for symmetry!"<<std::endl;
750 DUNE_THROW(
ISTLError,
"Matrix not symmetric!");
754 template<
typename M,
typename C>
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);
764 for (std::size_t i=0; i < newComm.indexSet().size(); i++) {
771 for (std::size_t i=0; i < origComm.indexSet().size(); i++)
779 newComm.indexSet(), backwardscopyrowsize);
781 newComm.indexSet(),copyrowsize);
782 ri.template redistribute<MatrixRowGatherScatter<M,IndexSet> >(origrow_copy,
784 ri.getInterface().free();
785 RemoteIndices<IndexSet> *ris =
new RemoteIndices<IndexSet>(origComm.indexSet(),
787 origComm.communicator());
788 ris->template rebuild<true>();
789 ri.getInterface().build(*ris,ownerflags,ownerflags);
793 origrow(origMatrix, origComm.globalLookup(), newComm.indexSet());
795 newrow(newMatrix, origComm.globalLookup(), newComm.indexSet(),rowsize);
796 ri.template redistribute<MatrixRowGatherScatter<M,IndexSet> >(origrow,newrow);
817 template<
typename M,
typename C>
835 DUNE_THROW(InvalidStateException,
"Trying to redistribute in sequential program!");
843 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:9
void redistributeMatrixEntries(M &origMatrix, M &newMatrix, C &origComm, C &newComm, RedistributeInformation< C > &ri)
Definition: matrixredistribute.hh:755
void redistributeSparsityPattern(M &origMatrix, M &newMatrix, C &origComm, C &newComm, RedistributeInformation< C > &ri)
Definition: matrixredistribute.hh:661
void redistributeMatrix(M &origMatrix, M &newMatrix, C &origComm, C &newComm, RedistributeInformation< C > &ri)
Redistribute a matrix according to given domain decompositions.
Definition: matrixredistribute.hh:818
derive error class from the base class in common
Definition: istlexception.hh:17
Definition: matrixredistribute.hh:20
void setNoBackwardsCopyRows(std::size_t size)
Definition: matrixredistribute.hh:42
void redistribute(const D &from, D &to) const
Definition: matrixredistribute.hh:26
void resetSetup()
Definition: matrixredistribute.hh:33
void setNoCopyRows(std::size_t size)
Definition: matrixredistribute.hh:39
bool isSetup() const
Definition: matrixredistribute.hh:21
void setNoRows(std::size_t size)
Definition: matrixredistribute.hh:36
void redistributeBackward(D &from, const D &to) const
Definition: matrixredistribute.hh:30
std::size_t getBackwardsCopyRowSize(std::size_t index) const
Definition: matrixredistribute.hh:55
std::size_t getRowSize(std::size_t index) const
Definition: matrixredistribute.hh:45
std::size_t getCopyRowSize(std::size_t index) const
Definition: matrixredistribute.hh:50
std::size_t getRowSize(std::size_t index) const
Definition: matrixredistribute.hh:156
std::size_t & getBackwardsCopyRowSize(std::size_t index)
Definition: matrixredistribute.hh:171
RedistributeInterface & getInterface()
Definition: matrixredistribute.hh:73
void redistribute(const D &from, D &to) const
Definition: matrixredistribute.hh:134
void setNoBackwardsCopyRows(std::size_t rows)
Definition: matrixredistribute.hh:191
std::size_t & getCopyRowSize(std::size_t index)
Definition: matrixredistribute.hh:161
RedistributeInformation()
Definition: matrixredistribute.hh:69
std::size_t getCopyRowSize(std::size_t index) const
Definition: matrixredistribute.hh:166
void setNoRows(std::size_t rows)
Definition: matrixredistribute.hh:181
void reserve(std::size_t size)
Definition: matrixredistribute.hh:148
OwnerOverlapCopyCommunication< T, T1 > Comm
Definition: matrixredistribute.hh:67
void setNoCopyRows(std::size_t rows)
Definition: matrixredistribute.hh:186
void redistributeBackward(D &from, const D &to) const
Definition: matrixredistribute.hh:139
void setSetup()
Definition: matrixredistribute.hh:104
std::size_t getBackwardsCopyRowSize(std::size_t index) const
Definition: matrixredistribute.hh:176
void redistribute(const D &from, D &to) const
Definition: matrixredistribute.hh:116
void resetSetup()
Definition: matrixredistribute.hh:110
void redistributeBackward(D &from, const D &to) const
Definition: matrixredistribute.hh:124
void checkInterface(const IS &source, const IS &target, MPI_Comm comm)
Definition: matrixredistribute.hh:78
bool isSetup() const
Definition: matrixredistribute.hh:143
std::size_t & getRowSize(std::size_t index)
Definition: matrixredistribute.hh:151
Utility class to communicate and set the row sizes of a redistributed matrix.
Definition: matrixredistribute.hh:214
M::size_type size_type
Definition: matrixredistribute.hh:217
M::size_type value_type
Definition: matrixredistribute.hh:216
RI & rowsize
Definition: matrixredistribute.hh:228
const M & matrix
Definition: matrixredistribute.hh:227
CommMatrixRowSize(const M &m_, RI &rowsize_)
Constructor.
Definition: matrixredistribute.hh:224
Utility class to communicate and build the sparsity pattern of a redistributed matrix.
Definition: matrixredistribute.hh:243
M::size_type size_type
Definition: matrixredistribute.hh:244
const Dune::GlobalLookupIndexSet< I > & idxset
Definition: matrixredistribute.hh:354
void storeSparsityPattern(M &m)
Creates and stores the sparsity pattern of the redistributed matrix.
Definition: matrixredistribute.hh:274
const I & aggidxset
Definition: matrixredistribute.hh:355
const std::vector< size_type > * rowsize
Definition: matrixredistribute.hh:357
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:338
CommMatrixSparsityPattern(const M &m_, const Dune::GlobalLookupIndexSet< I > &idxset_, const I &aggidxset_)
Constructor for the original side.
Definition: matrixredistribute.hh:252
const M & matrix
Definition: matrixredistribute.hh:352
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:263
std::vector< std::set< size_type > > sparsity
Definition: matrixredistribute.hh:356
Dune::GlobalLookupIndexSet< I > LookupIndexSet
Definition: matrixredistribute.hh:353
static M::size_type getSize(const Type &t, std::size_t i)
Definition: matrixredistribute.hh:374
CommMatrixSparsityPattern< M, I > Type
Definition: matrixredistribute.hh:363
I::GlobalIndex IndexedType
The indexed type we send. This is the global index indentitfying the column.
Definition: matrixredistribute.hh:369
VariableSize IndexedTypeFlag
Each row varies in size.
Definition: matrixredistribute.hh:372
Utility class for comunicating the matrix entries.
Definition: matrixredistribute.hh:394
std::vector< size_t > * rowsize
row size information for the receiving side.
Definition: matrixredistribute.hh:450
M & matrix
The matrix to communicate the values of.
Definition: matrixredistribute.hh:444
CommMatrixRow(M &m_, const Dune::GlobalLookupIndexSet< I > &idxset_, const I &aggidxset_, std::vector< size_t > &rowsize_)
Constructor.
Definition: matrixredistribute.hh:410
const Dune::GlobalLookupIndexSet< I > & idxset
Index set for the original matrix.
Definition: matrixredistribute.hh:446
void setOverlapRowsToDirichlet()
Sets the non-owner rows correctly as Dirichlet boundaries.
Definition: matrixredistribute.hh:419
const I & aggidxset
Index set for the redistributed matrix.
Definition: matrixredistribute.hh:448
CommMatrixRow(M &m_, const Dune::GlobalLookupIndexSet< I > &idxset_, const I &aggidxset_)
Constructor.
Definition: matrixredistribute.hh:403
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:462
CommMatrixRow< M, I > Type
Definition: matrixredistribute.hh:456
static std::size_t getSize(const Type &t, std::size_t i)
Definition: matrixredistribute.hh:467
VariableSize IndexedTypeFlag
Each row varies in size.
Definition: matrixredistribute.hh:465
Definition: matrixredistribute.hh:481
static void scatter(Container &cont, const typename M::size_type &rowsize, std::size_t i)
Definition: matrixredistribute.hh:488
static const M::size_type gather(const Container &cont, std::size_t i)
Definition: matrixredistribute.hh:484
CommMatrixRowSize< M, RI > Container
Definition: matrixredistribute.hh:482
Definition: matrixredistribute.hh:498
static const M::size_type gather(const Container &cont, std::size_t i)
Definition: matrixredistribute.hh:501
static void scatter(Container &cont, const typename M::size_type &rowsize, std::size_t i)
Definition: matrixredistribute.hh:505
CommMatrixRowSize< M, RI > Container
Definition: matrixredistribute.hh:499
Definition: matrixredistribute.hh:516
M::ConstColIterator ColIter
Definition: matrixredistribute.hh:519
static void scatter(Container &cont, const GlobalIndex &gi, std::size_t i, std::size_t j)
Definition: matrixredistribute.hh:551
CommMatrixSparsityPattern< M, I > Container
Definition: matrixredistribute.hh:518
static GlobalIndex numlimits
Definition: matrixredistribute.hh:522
static ColIter col
Definition: matrixredistribute.hh:521
I::GlobalIndex GlobalIndex
Definition: matrixredistribute.hh:517
static const GlobalIndex & gather(const Container &cont, std::size_t i, std::size_t j)
Definition: matrixredistribute.hh:524
Definition: matrixredistribute.hh:597
I::GlobalIndex GlobalIndex
Definition: matrixredistribute.hh:598
static Data datastore
Definition: matrixredistribute.hh:603
static GlobalIndex numlimits
Definition: matrixredistribute.hh:604
M::ConstColIterator ColIter
Definition: matrixredistribute.hh:600
static const Data & gather(const Container &cont, std::size_t i, std::size_t j)
Definition: matrixredistribute.hh:606
std::pair< GlobalIndex, typename M::block_type > Data
Definition: matrixredistribute.hh:601
static void scatter(Container &cont, const Data &data, std::size_t i, std::size_t j)
Definition: matrixredistribute.hh:636
static ColIter col
Definition: matrixredistribute.hh:602
CommMatrixRow< M, I > Container
Definition: matrixredistribute.hh:599
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:258
@ 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