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>
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>
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);
164 void reserve(std::size_t size)
169 return rowSize[index];
172 std::size_t
getRowSize(std::size_t index)
const
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;
216 RedistributeInterface interface;
228 template<
class M,
class RI>
229 struct CommMatrixRowSize
232 typedef typename M::size_type value_type;
233 typedef typename M::size_type size_type;
240 CommMatrixRowSize(
const M& m_, RI& rowsize_)
241 : matrix(m_), rowsize(rowsize_)
257 template<
class M,
class I>
258 struct CommMatrixSparsityPattern
260 typedef typename M::size_type size_type;
268 CommMatrixSparsityPattern(
const M& m_,
const Dune::GlobalLookupIndexSet<I>& idxset_,
const I& aggidxset_)
269 : matrix(m_), idxset(idxset_), aggidxset(aggidxset_), rowsize()
279 CommMatrixSparsityPattern(
const M& m_,
const Dune::GlobalLookupIndexSet<I>& idxset_,
const I& aggidxset_,
280 const std::vector<typename M::size_type>& rowsize_)
281 : matrix(m_), idxset(idxset_), aggidxset(aggidxset_), sparsity(aggidxset_.size()), rowsize(&rowsize_)
290 void storeSparsityPattern(M& m)
293 typedef typename Dune::GlobalLookupIndexSet<I>::const_iterator IIter;
294 typedef typename Dune::OwnerOverlapCopyCommunication<int>::OwnerSet OwnerSet;
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();
320 Dune::GlobalLookupIndexSet<I> global(aggidxset);
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;
354 void completeSparsityPattern(std::vector<std::set<size_type> > add_sparsity)
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);
369 typedef Dune::GlobalLookupIndexSet<I> LookupIndexSet;
370 const Dune::GlobalLookupIndexSet<I>& idxset;
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;
385 typedef typename I::GlobalIndex IndexedType;
388 typedef VariableSize IndexedTypeFlag;
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>
419 CommMatrixRow(M& m_,
const Dune::GlobalLookupIndexSet<I>& idxset_,
const I& aggidxset_)
420 : matrix(m_), idxset(idxset_), aggidxset(aggidxset_), rowsize()
426 CommMatrixRow(M& m_,
const Dune::GlobalLookupIndexSet<I>& idxset_,
const I& aggidxset_,
427 std::vector<size_t>& rowsize_)
428 : matrix(m_), idxset(idxset_), aggidxset(aggidxset_), rowsize(&rowsize_)
435 void setOverlapRowsToDirichlet()
437 typedef typename Dune::GlobalLookupIndexSet<I>::const_iterator Iter;
438 typedef typename Dune::OwnerOverlapCopyCommunication<int>::OwnerSet OwnerSet;
440 for(Iter i= aggidxset.begin(), end=aggidxset.end(); i!=end; ++i)
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;
466 std::vector<size_t>* rowsize;
469 template<
class M,
class I>
470 struct CommPolicy<CommMatrixRow<M,I> >
472 typedef CommMatrixRow<M,I> Type;
478 typedef std::pair<typename I::GlobalIndex,typename M::block_type> IndexedType;
481 typedef VariableSize IndexedTypeFlag;
483 static std::size_t getSize(
const Type& t, std::size_t i)
486 return t.matrix[i].size();
489 assert((*t.rowsize)[i]>0);
490 return (*t.rowsize)[i];
495 template<
class M,
class I,
class RI>
496 struct MatrixRowSizeGatherScatter
498 typedef CommMatrixRowSize<M,RI> Container;
500 static const typename M::size_type gather(
const Container& cont, std::size_t i)
502 return cont.matrix[i].size();
504 static void scatter(Container& cont,
const typename M::size_type& rowsize, std::size_t i)
507 cont.rowsize.getRowSize(i)=rowsize;
512 template<
class M,
class I,
class RI>
513 struct MatrixCopyRowSizeGatherScatter
515 typedef CommMatrixRowSize<M,RI> Container;
517 static const typename M::size_type gather(
const Container& cont, std::size_t i)
519 return cont.matrix[i].size();
521 static void scatter(Container& cont,
const typename M::size_type& rowsize, std::size_t i)
524 if (rowsize > cont.rowsize.getCopyRowSize(i))
525 cont.rowsize.getCopyRowSize(i)=rowsize;
530 template<
class M,
class I>
531 struct MatrixSparsityPatternGatherScatter
533 typedef typename I::GlobalIndex GlobalIndex;
534 typedef CommMatrixSparsityPattern<M,I> Container;
535 typedef typename M::ConstColIterator ColIter;
538 static GlobalIndex numlimits;
540 static const GlobalIndex& gather(
const Container& cont, std::size_t i, std::size_t j)
543 col=cont.matrix[i].begin();
544 else if (
col!=cont.matrix[i].end())
551 if (
col==cont.matrix[i].end()) {
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();
567 static void scatter(Container& cont,
const GlobalIndex& gi, std::size_t i, std::size_t j)
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();
575 cont.sparsity[i].insert(column);
577 typedef typename Dune::OwnerOverlapCopyCommunication<int>::OwnerSet OwnerSet;
578 if(!OwnerSet::contains(ip.local().attribute()))
580 cont.sparsity[column].insert(i);
583 catch(
const Dune::RangeError&) {
586 typedef typename Container::LookupIndexSet GlobalLookup;
587 typedef typename GlobalLookup::IndexPair IndexPair;
588 typedef typename Dune::OwnerOverlapCopyCommunication<int>::OwnerSet OwnerSet;
590 GlobalLookup lookup(cont.aggidxset);
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>
609 typename MatrixSparsityPatternGatherScatter<M,I>::GlobalIndex MatrixSparsityPatternGatherScatter<M,I>::numlimits;
612 template<
class M,
class I>
613 struct MatrixRowGatherScatter
615 typedef typename I::GlobalIndex GlobalIndex;
616 typedef CommMatrixRow<M,I> Container;
617 typedef typename M::ConstColIterator ColIter;
618 typedef typename std::pair<GlobalIndex,typename M::block_type> Data;
620 static Data datastore;
621 static GlobalIndex numlimits;
623 static const Data& gather(
const Container& cont, std::size_t i, std::size_t j)
626 col=cont.matrix[i].begin();
627 else if (
col!=cont.matrix[i].end())
633 if (
col==cont.matrix[i].end()) {
634 numlimits = std::numeric_limits<GlobalIndex>::max();
635 datastore = std::make_pair(numlimits,*
col);
640 const typename I::IndexPair* index=cont.idxset.pair(
col.index());
644 if ( index->local().attribute() != 2)
645 datastore = std::make_pair(index->global(),*
col);
647 numlimits = std::numeric_limits<GlobalIndex>::max();
648 datastore = std::make_pair(numlimits,*
col);
653 static void scatter(Container& cont,
const Data& data, std::size_t i, std::size_t j)
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>
673 typename MatrixRowGatherScatter<M,I>::Data MatrixRowGatherScatter<M,I>::datastore;
675 template<
class M,
class I>
676 typename MatrixRowGatherScatter<M,I>::GlobalIndex MatrixRowGatherScatter<M,I>::numlimits;
678 template<
typename M,
typename C>
679 void redistributeSparsityPattern(M& origMatrix, M& newMatrix, C& origComm, C& newComm,
680 RedistributeInformation<C>& ri)
682 typename C::CopySet copyflags;
683 typename C::OwnerSet ownerflags;
684 typedef typename C::ParallelIndexSet IndexSet;
685 typedef RedistributeInformation<C> RI;
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);
691 CommMatrixRowSize<M,RI> commRowSize(origMatrix, ri);
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++) {
697 rowsize[i] = ri.getRowSize(i);
700 CommMatrixSparsityPattern<M,IndexSet>
701 origsp(origMatrix, origComm.globalLookup(), newComm.indexSet());
702 CommMatrixSparsityPattern<M,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);
718 CommMatrixRowSize<M,RI> commRowSize_copy(origMatrix, ri);
719 ri.template redistribute<MatrixCopyRowSizeGatherScatter<M,IndexSet,RI> >(commRowSize_copy,
722 for (std::size_t i=0; i < newComm.indexSet().size(); i++) {
723 copyrowsize[i] = ri.getCopyRowSize(i);
726 ri.redistributeBackward(backwardscopyrowsize,copyrowsize);
727 for (std::size_t i=0; i < origComm.indexSet().size(); i++) {
728 ri.getBackwardsCopyRowSize(i) = backwardscopyrowsize[i];
732 CommMatrixSparsityPattern<M,IndexSet> origsp_copy(origMatrix,
733 origComm.globalLookup(),
735 backwardscopyrowsize);
736 CommMatrixSparsityPattern<M,IndexSet> newsp_copy(origMatrix, origComm.globalLookup(),
737 newComm.indexSet(), copyrowsize);
738 ri.template redistribute<MatrixSparsityPatternGatherScatter<M,IndexSet> >(origsp_copy,
741 newsp.completeSparsityPattern(newsp_copy.sparsity);
742 newsp.storeSparsityPattern(newMatrix);
745 newsp.storeSparsityPattern(newMatrix);
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>
774 RedistributeInformation<C>& ri)
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++) {
783 rowsize[i] = ri.getRowSize(i);
785 copyrowsize[i] = ri.getCopyRowSize(i);
789 for (std::size_t i=0; i < origComm.indexSet().size(); i++)
791 backwardscopyrowsize[i] = ri.getBackwardsCopyRowSize(i);
796 CommMatrixRow<M,IndexSet> origrow_copy(origMatrix, origComm.globalLookup(),
797 newComm.indexSet(), backwardscopyrowsize);
798 CommMatrixRow<M,IndexSet> newrow_copy(newMatrix, origComm.globalLookup(),
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);
810 CommMatrixRow<M,IndexSet>
811 origrow(origMatrix, origComm.globalLookup(), newComm.indexSet());
812 CommMatrixRow<M,IndexSet>
813 newrow(newMatrix, origComm.globalLookup(), newComm.indexSet(),rowsize);
814 ri.template redistribute<MatrixRowGatherScatter<M,IndexSet> >(origrow,newrow);
816 newrow.setOverlapRowsToDirichlet();
835 template<
typename M,
typename C>
837 RedistributeInformation<C>& ri)
839 ri.setNoRows(newComm.indexSet().size());
840 ri.setNoCopyRows(newComm.indexSet().size());
841 ri.setNoBackwardsCopyRows(origComm.indexSet().size());
842 redistributeSparsityPattern(origMatrix, newMatrix, origComm, newComm, ri);
853 DUNE_THROW(InvalidStateException,
"Trying to redistribute in sequential program!");
861 DUNE_THROW(InvalidStateException,
"Trying to redistribute in sequential program!");