3 #ifndef DUNE_MATRIXREDIST_HH
4 #define DUNE_MATRIXREDIST_HH
6 #include<dune/common/exceptions.hh>
7 #include<dune/common/parallel/indexset.hh>
61 template<
typename T,
typename T1>
68 : interface(), setup_(false)
76 void checkInterface(
const IS& source,
77 const IS& target, MPI_Comm comm)
79 RemoteIndices<IS> *ri=
new RemoteIndices<IS>(source, target, comm);
80 ri->template rebuild<true>();
84 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
86 inf.build(*ri, flags, flags);
93 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
95 std::cout<<
"Interfaces do not match!"<<std::endl;
96 std::cout<<rank<<
": redist interface new :"<<inf<<std::endl;
97 std::cout<<rank<<
": redist interface :"<<interface<<std::endl;
117 template<
class GatherScatter,
class D>
120 BufferedCommunicator communicator;
121 communicator.template build<D>(from,to, interface);
122 communicator.template forward<GatherScatter>(from, to);
125 template<
class GatherScatter,
class D>
129 BufferedCommunicator communicator;
130 communicator.template build<D>(from,to, interface);
131 communicator.template backward<GatherScatter>(from, to);
138 redistribute<CopyGatherScatter<D> >(from,to);
143 redistributeBackward<CopyGatherScatter<D> >(from,to);
150 void reserve(std::size_t size)
155 return rowSize[
index];
160 return rowSize[
index];
165 return copyrowSize[
index];
170 return copyrowSize[
index];
175 return backwardscopyrowSize[
index];
180 return backwardscopyrowSize[
index];
185 rowSize.resize(rows, 0);
190 copyrowSize.resize(rows, 0);
195 backwardscopyrowSize.resize(rows, 0);
199 std::vector<std::size_t> rowSize;
200 std::vector<std::size_t> copyrowSize;
201 std::vector<std::size_t> backwardscopyrowSize;
214 template<
class M,
class RI>
243 template<
class M,
class I>
266 const std::vector<typename M::size_type>& rowsize_)
279 typedef typename Dune::GlobalLookupIndexSet<I>::const_iterator IIter;
285 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
288 if(!OwnerSet::contains(i->local().attribute())){
290 std::cout<<rank<<
" Inserting diagonal for"<<i->local()<<std::endl;
292 sparsity[i->local()].insert(i->local());
301 m.setBuildMode(M::row_wise);
302 typename M::CreateIterator citer=m.createbegin();
306 Dune::GlobalLookupIndexSet<I> global(
aggidxset);
308 typedef typename std::vector<std::set<size_type> >::const_iterator Iter;
311 typedef typename std::set<size_type>::const_iterator SIter;
312 for(SIter si=i->begin(), send=i->end(); si!=send; ++si)
315 if(i->find(idx)==i->end()){
316 const typename I::IndexPair* gi=global.pair(idx);
318 std::cout<<rank<<
": row "<<idx<<
" is missing a diagonal entry! global="<<gi->global()<<
" attr="<<gi->local().attribute()<<
" "<<
319 OwnerSet::contains(gi->local().attribute())<<
320 " row size="<<i->size()<<std::endl;
342 for (
unsigned int i = 0; i !=
sparsity.size(); ++i) {
343 if (add_sparsity[i].size() != 0){
344 typedef std::set<size_type> Set;
346 std::insert_iterator<Set> tmp_insert (tmp_set, tmp_set.begin());
347 std::set_union(add_sparsity[i].begin(), add_sparsity[i].end(),
356 const Dune::GlobalLookupIndexSet<I>&
idxset;
362 template<
class M,
class I>
376 static typename M::size_type getSize(
const Type& t, std::size_t i)
379 return t.
matrix[i].size();
394 template<
class M,
class I>
405 CommMatrixRow(M& m_,
const Dune::GlobalLookupIndexSet<I>& idxset_,
const I& aggidxset_)
412 CommMatrixRow(M& m_,
const Dune::GlobalLookupIndexSet<I>& idxset_,
const I& aggidxset_,
413 std::vector<size_t>& rowsize_)
423 typedef typename Dune::GlobalLookupIndexSet<I>::const_iterator Iter;
427 if(!OwnerSet::contains(i->local().attribute())){
429 typedef typename M::ColIterator CIter;
430 for(CIter c=
matrix[i->local()].begin(), cend=
matrix[i->local()].end();
434 if(c.index()==i->local()){
435 typedef typename M::block_type::RowIterator RIter;
436 for(RIter r=c->begin(), rend=c->end();
446 const Dune::GlobalLookupIndexSet<I>&
idxset;
453 template<
class M,
class I>
462 typedef std::pair<typename I::GlobalIndex,typename M::block_type>
IndexedType;
467 static std::size_t getSize(
const Type& t, std::size_t i)
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();
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
col = ip.local();
561 if(!OwnerSet::contains(ip.local().attribute()))
566 catch(Dune::RangeError er){
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>
592 typename MatrixSparsityPatternGatherScatter<M,I>::GlobalIndex MatrixSparsityPatternGatherScatter<M,I>::numlimits;
595 template<
class M,
class I>
601 typedef typename std::pair<GlobalIndex,typename M::block_type>
Data;
617 numlimits = std::numeric_limits<GlobalIndex>::max();
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(Dune::RangeError er){
651 template<
class M,
class I>
654 template<
class M,
class I>
655 typename MatrixRowGatherScatter<M,I>::Data MatrixRowGatherScatter<M,I>::datastore;
657 template<
class M,
class I>
658 typename MatrixRowGatherScatter<M,I>::GlobalIndex MatrixRowGatherScatter<M,I>::numlimits;
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,
723 newsp.completeSparsityPattern(newsp_copy.sparsity);
724 newsp.storeSparsityPattern(newMatrix);
727 newsp.storeSparsityPattern(newMatrix);
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;
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);
798 newrow.setOverlapRowsToDirichlet();
799 if(newMatrix.N()>0&&newMatrix.N()<20)
800 printmatrix(std::cout, newMatrix,
"redist",
"row");
819 template<
typename M,
typename C>