3 #ifndef DUNE_MATRIXREDIST_HH
4 #define DUNE_MATRIXREDIST_HH
6 #include <dune/common/exceptions.hh>
7 #include <dune/common/parallel/indexset.hh>
8 #include <dune/common/unused.hh>
27 DUNE_UNUSED_PARAMETER(from);
28 DUNE_UNUSED_PARAMETER(to);
34 DUNE_UNUSED_PARAMETER(from);
35 DUNE_UNUSED_PARAMETER(to);
43 DUNE_UNUSED_PARAMETER(size);
48 DUNE_UNUSED_PARAMETER(size);
53 DUNE_UNUSED_PARAMETER(size);
58 DUNE_UNUSED_PARAMETER(index);
64 DUNE_UNUSED_PARAMETER(index);
70 DUNE_UNUSED_PARAMETER(index);
77 template<
typename T,
typename T1>
84 : interface(), setup_(false)
92 void checkInterface(
const IS& source,
93 const IS& target, MPI_Comm comm)
95 RemoteIndices<IS> *ri=
new RemoteIndices<IS>(source, target, comm);
96 ri->template rebuild<true>();
100 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
102 inf.build(*ri, flags, flags);
109 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
111 std::cout<<
"Interfaces do not match!"<<std::endl;
112 std::cout<<rank<<
": redist interface new :"<<inf<<std::endl;
113 std::cout<<rank<<
": redist interface :"<<interface<<std::endl;
133 template<
class GatherScatter,
class D>
136 BufferedCommunicator communicator;
137 communicator.template build<D>(from,to, interface);
138 communicator.template forward<GatherScatter>(from, to);
141 template<
class GatherScatter,
class D>
145 BufferedCommunicator communicator;
146 communicator.template build<D>(from,to, interface);
147 communicator.template backward<GatherScatter>(from, to);
154 redistribute<CopyGatherScatter<D> >(from,to);
159 redistributeBackward<CopyGatherScatter<D> >(from,to);
166 void reserve(std::size_t size)
171 return rowSize[
index];
176 return rowSize[
index];
181 return copyrowSize[
index];
186 return copyrowSize[
index];
191 return backwardscopyrowSize[
index];
196 return backwardscopyrowSize[
index];
201 rowSize.resize(rows, 0);
206 copyrowSize.resize(rows, 0);
211 backwardscopyrowSize.resize(rows, 0);
215 std::vector<std::size_t> rowSize;
216 std::vector<std::size_t> copyrowSize;
217 std::vector<std::size_t> backwardscopyrowSize;
230 template<
class M,
class RI>
259 template<
class M,
class I>
282 const std::vector<typename M::size_type>& rowsize_)
295 typedef typename Dune::GlobalLookupIndexSet<I>::const_iterator IIter;
301 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
304 if(!OwnerSet::contains(i->local().attribute())) {
306 std::cout<<rank<<
" Inserting diagonal for"<<i->local()<<std::endl;
308 sparsity[i->local()].insert(i->local());
317 m.setBuildMode(M::row_wise);
318 typename M::CreateIterator citer=m.createbegin();
322 Dune::GlobalLookupIndexSet<I> global(
aggidxset);
324 typedef typename std::vector<std::set<size_type> >::const_iterator Iter;
327 typedef typename std::set<size_type>::const_iterator SIter;
328 for(SIter si=i->begin(), send=i->end(); si!=send; ++si)
331 if(i->find(idx)==i->end()) {
332 const typename I::IndexPair* gi=global.pair(idx);
334 std::cout<<rank<<
": row "<<idx<<
" is missing a diagonal entry! global="<<gi->global()<<
" attr="<<gi->local().attribute()<<
" "<<
335 OwnerSet::contains(gi->local().attribute())<<
336 " row size="<<i->size()<<std::endl;
358 for (
unsigned int i = 0; i !=
sparsity.size(); ++i) {
359 if (add_sparsity[i].size() != 0) {
360 typedef std::set<size_type> Set;
362 std::insert_iterator<Set> tmp_insert (tmp_set, tmp_set.begin());
363 std::set_union(add_sparsity[i].begin(), add_sparsity[i].end(),
372 const Dune::GlobalLookupIndexSet<I>&
idxset;
378 template<
class M,
class I>
392 static typename M::size_type getSize(
const Type& t, std::size_t i)
395 return t.
matrix[i].size();
410 template<
class M,
class I>
421 CommMatrixRow(M& m_,
const Dune::GlobalLookupIndexSet<I>& idxset_,
const I& aggidxset_)
428 CommMatrixRow(M& m_,
const Dune::GlobalLookupIndexSet<I>& idxset_,
const I& aggidxset_,
429 std::vector<size_t>& rowsize_)
439 typedef typename Dune::GlobalLookupIndexSet<I>::const_iterator Iter;
443 if(!OwnerSet::contains(i->local().attribute())) {
445 typedef typename M::ColIterator CIter;
446 for(CIter c=
matrix[i->local()].begin(), cend=
matrix[i->local()].end();
450 if(c.index()==i->local()) {
451 typedef typename M::block_type::RowIterator RIter;
452 for(RIter r=c->begin(), rend=c->end();
462 const Dune::GlobalLookupIndexSet<I>&
idxset;
469 template<
class M,
class I>
478 typedef std::pair<typename I::GlobalIndex,typename M::block_type>
IndexedType;
483 static std::size_t getSize(
const Type& t, std::size_t i)
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();
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
col = ip.local();
578 if(!OwnerSet::contains(ip.local().attribute()))
583 catch(Dune::RangeError er) {
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>
609 typename MatrixSparsityPatternGatherScatter<M,I>::GlobalIndex MatrixSparsityPatternGatherScatter<M,I>::numlimits;
612 template<
class M,
class I>
618 typedef typename std::pair<GlobalIndex,typename M::block_type>
Data;
634 numlimits = std::numeric_limits<GlobalIndex>::max();
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(Dune::RangeError er) {
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>
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,
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;
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);
816 newrow.setOverlapRowsToDirichlet();
817 if(newMatrix.N()>0&&newMatrix.N()<20)
818 printmatrix(std::cout, newMatrix,
"redist",
"row");
837 template<
typename M,
typename C>