3#ifndef DUNE_ISTL_MATRIXREDISTRIBUTE_HH
4#define DUNE_ISTL_MATRIXREDISTRIBUTE_HH
11#include <dune/istl/paamg/pinfo.hh>
20 struct RedistributeInformation
27 void redistribute(
const D& from, D& to)
const
34 void redistributeBackward(D& from,
const D& to)
const
43 void setNoRows(std::size_t size)
48 void setNoCopyRows(std::size_t size)
53 void setNoBackwardsCopyRows(std::size_t size)
58 std::size_t getRowSize(std::size_t index)
const
64 std::size_t getCopyRowSize(std::size_t index)
const
70 std::size_t getBackwardsCopyRowSize(std::size_t index)
const
79 template<
typename T,
typename T1>
80 class RedistributeInformation<OwnerOverlapCopyCommunication<T,T1> >
83 typedef OwnerOverlapCopyCommunication<T,T1> Comm;
85 RedistributeInformation()
86 : interface(), setup_(false)
89 RedistributeInterface& getInterface()
94 void checkInterface(
const IS& source,
95 const IS& target, MPI_Comm comm)
97 auto ri = std::make_unique<RemoteIndices<IS> >(source, target, comm);
98 ri->template rebuild<true>();
100 typename OwnerOverlapCopyCommunication<int>::OwnerSet flags;
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>
132 void redistribute(
const D& from, D& to)
const
134 BufferedCommunicator communicator;
135 communicator.template build<D>(from,to, interface);
136 communicator.template forward<GatherScatter>(from, to);
139 template<
class GatherScatter,
class D>
140 void redistributeBackward(D& from,
const D& to)
const
143 BufferedCommunicator communicator;
144 communicator.template build<D>(from,to, interface);
145 communicator.template backward<GatherScatter>(from, to);
150 void redistribute(
const D& from, D& to)
const
152 redistribute<CopyGatherScatter<D> >(from,to);
155 void redistributeBackward(D& from,
const D& to)
const
157 redistributeBackward<CopyGatherScatter<D> >(from,to);
164 void reserve(std::size_t size)
167 std::size_t& getRowSize(std::size_t index)
169 return rowSize[index];
172 std::size_t getRowSize(std::size_t index)
const
174 return rowSize[index];
177 std::size_t& getCopyRowSize(std::size_t index)
179 return copyrowSize[index];
182 std::size_t getCopyRowSize(std::size_t index)
const
184 return copyrowSize[index];
187 std::size_t& getBackwardsCopyRowSize(std::size_t index)
189 return backwardscopyrowSize[index];
192 std::size_t getBackwardsCopyRowSize(std::size_t index)
const
194 return backwardscopyrowSize[index];
197 void setNoRows(std::size_t rows)
199 rowSize.resize(rows, 0);
202 void setNoCopyRows(std::size_t rows)
204 copyrowSize.resize(rows, 0);
207 void setNoBackwardsCopyRows(std::size_t rows)
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;
216 RedistributeInterface interface;
228 template<
class M,
class RI>
232 typedef typename M::size_type value_type;
233 typedef typename M::size_type size_type;
241 : matrix(m_), rowsize(rowsize_)
257 template<
class M,
class I>
260 typedef typename M::size_type size_type;
269 : matrix(m_), idxset(idxset_), aggidxset(aggidxset_), rowsize()
280 const std::vector<typename M::size_type>& rowsize_)
281 : matrix(m_), idxset(idxset_), aggidxset(aggidxset_), sparsity(aggidxset_.size()), rowsize(&rowsize_)
299 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
301 for(IIter i= aggidxset.begin(), end=aggidxset.end(); i!=end; ++i) {
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());
309 nnz+=sparsity[i->local()].size();
311 assert( aggidxset.size()==sparsity.size());
314 m.setSize(aggidxset.size(), aggidxset.size(), nnz);
315 m.setBuildMode(M::row_wise);
316 typename M::CreateIterator citer=m.createbegin();
322 typedef typename std::vector<std::set<size_type> >::const_iterator Iter;
323 for(Iter i=sparsity.begin(), end=sparsity.end(); i!=end; ++i, ++citer)
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(),
362 sparsity[i].begin(), sparsity[i].end(), tmp_insert);
363 sparsity[i].swap(tmp_set);
372 std::vector<std::set<size_type> > sparsity;
373 const std::vector<size_type>* rowsize;
376 template<
class M,
class I>
377 struct CommPolicy<CommMatrixSparsityPattern<M,I> >
379 typedef CommMatrixSparsityPattern<M,I>
Type;
390 static typename M::size_type
getSize(
const Type& t, std::size_t i)
393 return t.matrix[i].size();
396 assert((*t.rowsize)[i]>0);
397 return (*t.rowsize)[i];
408 template<
class M,
class I>
427 std::vector<size_t>& rowsize_)
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 typedef typename M::block_type::RowIterator RIter;
450 for(RIter r=c->begin(), rend=c->end();
467 template<
class M,
class I>
476 typedef std::pair<typename I::GlobalIndex,typename M::block_type>
IndexedType;
481 static std::size_t
getSize(
const Type& t, std::size_t i)
484 return t.matrix[i].size();
487 assert((*t.rowsize)[i]>0);
488 return (*t.rowsize)[i];
493 template<
class M,
class I,
class RI>
494 struct MatrixRowSizeGatherScatter
496 typedef CommMatrixRowSize<M,RI> Container;
498 static const typename M::size_type gather(
const Container& cont, std::size_t i)
500 return cont.matrix[i].size();
502 static void scatter(Container& cont,
const typename M::size_type& rowsize, std::size_t i)
505 cont.rowsize.getRowSize(i)=rowsize;
510 template<
class M,
class I,
class RI>
511 struct MatrixCopyRowSizeGatherScatter
513 typedef CommMatrixRowSize<M,RI> Container;
515 static const typename M::size_type gather(
const Container& cont, std::size_t i)
517 return cont.matrix[i].size();
519 static void scatter(Container& cont,
const typename M::size_type& rowsize, std::size_t i)
522 if (rowsize > cont.rowsize.getCopyRowSize(i))
523 cont.rowsize.getCopyRowSize(i)=rowsize;
528 template<
class M,
class I>
529 struct MatrixSparsityPatternGatherScatter
531 typedef typename I::GlobalIndex GlobalIndex;
532 typedef CommMatrixSparsityPattern<M,I> Container;
533 typedef typename M::ConstColIterator ColIter;
536 static GlobalIndex numlimits;
538 static const GlobalIndex& gather(
const Container& cont, std::size_t i, std::size_t j)
541 col=cont.matrix[i].begin();
542 else if (col!=cont.matrix[i].end())
549 if (col==cont.matrix[i].end()) {
550 numlimits = std::numeric_limits<GlobalIndex>::max();
554 const typename I::IndexPair* index=cont.idxset.pair(col.index());
557 if ( index->local().attribute() != 2)
558 return index->global();
560 numlimits = std::numeric_limits<GlobalIndex>::max();
565 static void scatter(Container& cont,
const GlobalIndex& gi, std::size_t i, std::size_t j)
569 if (gi != std::numeric_limits<GlobalIndex>::max()) {
570 const typename I::IndexPair& ip=cont.aggidxset.at(gi);
571 assert(ip.global()==gi);
572 std::size_t column = ip.local();
573 cont.sparsity[i].insert(column);
576 if(!OwnerSet::contains(ip.local().attribute()))
578 cont.sparsity[column].insert(i);
584 typedef typename Container::LookupIndexSet GlobalLookup;
585 typedef typename GlobalLookup::IndexPair IndexPair;
588 GlobalLookup lookup(cont.aggidxset);
589 const IndexPair* pi=lookup.pair(i);
591 if(OwnerSet::contains(pi->local().attribute())) {
593 MPI_Comm_rank(MPI_COMM_WORLD,&rank);
594 std::cout<<rank<<cont.aggidxset<<std::endl;
595 std::cout<<rank<<
": row "<<i<<
" (global="<<gi <<
") not in index set for owner index "<<pi->global()<<std::endl;
603 template<
class M,
class I>
604 typename MatrixSparsityPatternGatherScatter<M,I>::ColIter MatrixSparsityPatternGatherScatter<M,I>::col;
606 template<
class M,
class I>
607 typename MatrixSparsityPatternGatherScatter<M,I>::GlobalIndex MatrixSparsityPatternGatherScatter<M,I>::numlimits;
610 template<
class M,
class I>
611 struct MatrixRowGatherScatter
613 typedef typename I::GlobalIndex GlobalIndex;
614 typedef CommMatrixRow<M,I> Container;
615 typedef typename M::ConstColIterator ColIter;
616 typedef typename std::pair<GlobalIndex,typename M::block_type> Data;
618 static Data datastore;
619 static GlobalIndex numlimits;
621 static const Data& gather(
const Container& cont, std::size_t i, std::size_t j)
624 col=cont.matrix[i].begin();
625 else if (col!=cont.matrix[i].end())
631 if (col==cont.matrix[i].end()) {
632 numlimits = std::numeric_limits<GlobalIndex>::max();
633 datastore = std::make_pair(numlimits,*col);
638 const typename I::IndexPair* index=cont.idxset.pair(col.index());
642 if ( index->local().attribute() != 2)
643 datastore = std::make_pair(index->global(),*col);
645 numlimits = std::numeric_limits<GlobalIndex>::max();
646 datastore = std::make_pair(numlimits,*col);
651 static void scatter(Container& cont,
const Data& data, std::size_t i, std::size_t j)
655 if (data.first != std::numeric_limits<GlobalIndex>::max()) {
656 typename M::size_type column=cont.aggidxset.at(data.first).local();
657 cont.matrix[i][column]=data.second;
667 template<
class M,
class I>
668 typename MatrixRowGatherScatter<M,I>::ColIter MatrixRowGatherScatter<M,I>::col;
670 template<
class M,
class I>
671 typename MatrixRowGatherScatter<M,I>::Data MatrixRowGatherScatter<M,I>::datastore;
673 template<
class M,
class I>
674 typename MatrixRowGatherScatter<M,I>::GlobalIndex MatrixRowGatherScatter<M,I>::numlimits;
676 template<
typename M,
typename C>
677 void redistributeSparsityPattern(M& origMatrix, M& newMatrix, C& origComm, C& newComm,
678 RedistributeInformation<C>& ri)
680 typename C::CopySet copyflags;
681 typename C::OwnerSet ownerflags;
682 typedef typename C::ParallelIndexSet IndexSet;
683 typedef RedistributeInformation<C> RI;
684 std::vector<typename M::size_type> rowsize(newComm.indexSet().size(), 0);
685 std::vector<typename M::size_type> copyrowsize(newComm.indexSet().size(), 0);
686 std::vector<typename M::size_type> backwardscopyrowsize(origComm.indexSet().size(), 0);
689 CommMatrixRowSize<M,RI> commRowSize(origMatrix, ri);
690 ri.template redistribute<MatrixRowSizeGatherScatter<M,IndexSet,RI> >(commRowSize,commRowSize);
692 origComm.buildGlobalLookup();
694 for (std::size_t i=0; i < newComm.indexSet().size(); i++) {
695 rowsize[i] = ri.getRowSize(i);
698 CommMatrixSparsityPattern<M,IndexSet>
699 origsp(origMatrix, origComm.globalLookup(), newComm.indexSet());
700 CommMatrixSparsityPattern<M,IndexSet>
701 newsp(origMatrix, origComm.globalLookup(), newComm.indexSet(), rowsize);
703 ri.template redistribute<MatrixSparsityPatternGatherScatter<M,IndexSet> >(origsp,newsp);
707 RemoteIndices<IndexSet> *ris =
new RemoteIndices<IndexSet>(origComm.indexSet(),
709 origComm.communicator());
710 ris->template rebuild<true>();
712 ri.getInterface().free();
713 ri.getInterface().build(*ris,copyflags,ownerflags);
716 CommMatrixRowSize<M,RI> commRowSize_copy(origMatrix, ri);
717 ri.template redistribute<MatrixCopyRowSizeGatherScatter<M,IndexSet,RI> >(commRowSize_copy,
720 for (std::size_t i=0; i < newComm.indexSet().size(); i++) {
721 copyrowsize[i] = ri.getCopyRowSize(i);
724 ri.redistributeBackward(backwardscopyrowsize,copyrowsize);
725 for (std::size_t i=0; i < origComm.indexSet().size(); i++) {
726 ri.getBackwardsCopyRowSize(i) = backwardscopyrowsize[i];
730 CommMatrixSparsityPattern<M,IndexSet> origsp_copy(origMatrix,
731 origComm.globalLookup(),
733 backwardscopyrowsize);
734 CommMatrixSparsityPattern<M,IndexSet> newsp_copy(origMatrix, origComm.globalLookup(),
735 newComm.indexSet(), copyrowsize);
736 ri.template redistribute<MatrixSparsityPatternGatherScatter<M,IndexSet> >(origsp_copy,
739 newsp.completeSparsityPattern(newsp_copy.sparsity);
740 newsp.storeSparsityPattern(newMatrix);
743 newsp.storeSparsityPattern(newMatrix);
745#ifdef DUNE_ISTL_WITH_CHECKING
748 typedef typename M::ConstRowIterator RIter;
749 for(RIter row=newMatrix.begin(), rend=newMatrix.end(); row != rend; ++row) {
750 typedef typename M::ConstColIterator CIter;
751 for(CIter col=row->begin(), cend=row->end(); col!=cend; ++col)
754 newMatrix[col.index()][row.index()];
756 std::cerr<<newComm.communicator().rank()<<
": entry ("
757 <<col.index()<<
","<<row.index()<<
") missing! for symmetry!"<<std::endl;
766 DUNE_THROW(ISTLError,
"Matrix not symmetric!");
770 template<
typename M,
typename C>
771 void redistributeMatrixEntries(M& origMatrix, M& newMatrix, C& origComm, C& newComm,
772 RedistributeInformation<C>& ri)
774 typedef typename C::ParallelIndexSet IndexSet;
775 typename C::OwnerSet ownerflags;
776 std::vector<typename M::size_type> rowsize(newComm.indexSet().size(), 0);
777 std::vector<typename M::size_type> copyrowsize(newComm.indexSet().size(), 0);
778 std::vector<typename M::size_type> backwardscopyrowsize(origComm.indexSet().size(), 0);
780 for (std::size_t i=0; i < newComm.indexSet().size(); i++) {
781 rowsize[i] = ri.getRowSize(i);
783 copyrowsize[i] = ri.getCopyRowSize(i);
787 for (std::size_t i=0; i < origComm.indexSet().size(); i++)
789 backwardscopyrowsize[i] = ri.getBackwardsCopyRowSize(i);
794 CommMatrixRow<M,IndexSet> origrow_copy(origMatrix, origComm.globalLookup(),
795 newComm.indexSet(), backwardscopyrowsize);
796 CommMatrixRow<M,IndexSet> newrow_copy(newMatrix, origComm.globalLookup(),
797 newComm.indexSet(),copyrowsize);
798 ri.template redistribute<MatrixRowGatherScatter<M,IndexSet> >(origrow_copy,
800 ri.getInterface().free();
801 RemoteIndices<IndexSet> *ris =
new RemoteIndices<IndexSet>(origComm.indexSet(),
803 origComm.communicator());
804 ris->template rebuild<true>();
805 ri.getInterface().build(*ris,ownerflags,ownerflags);
808 CommMatrixRow<M,IndexSet>
809 origrow(origMatrix, origComm.globalLookup(), newComm.indexSet());
810 CommMatrixRow<M,IndexSet>
811 newrow(newMatrix, origComm.globalLookup(), newComm.indexSet(),rowsize);
812 ri.template redistribute<MatrixRowGatherScatter<M,IndexSet> >(origrow,newrow);
814 newrow.setOverlapRowsToDirichlet();
833 template<
typename M,
typename C>
835 RedistributeInformation<C>& ri)
837 ri.setNoRows(newComm.indexSet().size());
838 ri.setNoCopyRows(newComm.indexSet().size());
839 ri.setNoBackwardsCopyRows(origComm.indexSet().size());
840 redistributeSparsityPattern(origMatrix, newMatrix, origComm, newComm, ri);
841 redistributeMatrixEntries(origMatrix, newMatrix, origComm, newComm, ri);
846 void redistributeMatrixEntries(M& origMatrix, M& newMatrix,
847 Dune::Amg::SequentialInformation& origComm,
848 Dune::Amg::SequentialInformation& newComm,
849 RedistributeInformation<Dune::Amg::SequentialInformation>& ri)
851 DUNE_THROW(InvalidStateException,
"Trying to redistribute in sequential program!");
855 Dune::Amg::SequentialInformation& origComm,
856 Dune::Amg::SequentialInformation& newComm,
857 RedistributeInformation<Dune::Amg::SequentialInformation>& ri)
859 DUNE_THROW(InvalidStateException,
"Trying to redistribute in sequential program!");
A constant random access iterator for the Dune::ArrayList class.
Definition: arraylist.hh:379
A set consisting only of one item.
Definition: enumset.hh:60
Decorates an index set with the possibility to find a global index that is mapped to a specific local...
Definition: indexset.hh:506
derive error class from the base class in common
Definition: istlexception.hh:16
Default exception class for range errors.
Definition: exceptions.hh:252
A few common exception classes.
const IndexPair * pair(const std::size_t &local) const
Get the index pair corresponding to a local index.
#define DUNE_UNUSED_PARAMETER(parm)
A macro to mark intentionally unused function parameters with.
Definition: unused.hh:25
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
Provides a map between global and local indices.
Dune namespace.
Definition: alignedallocator.hh:10
void redistributeMatrix(M &origMatrix, M &newMatrix, C &origComm, C &newComm, RedistributeInformation< C > &ri)
Redistribute a matrix according to given domain decompositions.
Definition: matrixredistribute.hh:834
Classes providing communication interfaces for overlapping Schwarz methods.
Functionality for redistributing a parallel index set using graph partitioning.
Utility class to communicate and set the row sizes of a redistributed matrix.
Definition: matrixredistribute.hh:230
CommMatrixRowSize(const M &m_, RI &rowsize_)
Constructor.
Definition: matrixredistribute.hh:240
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:464
M & matrix
The matrix to communicate the values of.
Definition: matrixredistribute.hh:458
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:460
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:462
CommMatrixRow(M &m_, const Dune::GlobalLookupIndexSet< I > &idxset_, const I &aggidxset_)
Constructor.
Definition: matrixredistribute.hh:419
Utility class to communicate and build the sparsity pattern of a redistributed matrix.
Definition: matrixredistribute.hh:259
void storeSparsityPattern(M &m)
Creates and stores the sparsity pattern of the redistributed matrix.
Definition: matrixredistribute.hh:290
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
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
Default policy used for communicating an indexed type.
Definition: communicator.hh:127
V::value_type IndexedType
The type we get at each index with operator[].
Definition: communicator.hh:146
static int getSize(const V &, int index)
Get the number of primitve elements at that index.
SizeOne IndexedTypeFlag
Whether the indexed type has variable size or there is always one value at each index.
Definition: communicator.hh:152
V Type
The type the policy is for.
Definition: communicator.hh:139
@ nonoverlapping
Category for non-overlapping solvers.
Definition: solvercategory.hh:25
Flag for marking indexed data structures where the data at each index may be a variable multiple of a...
Definition: communicator.hh:117
Definition of the DUNE_UNUSED macro for the case that config.h is not available.