dune-istl  2.6-git
matrixredistribute.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #ifndef DUNE_ISTL_MATRIXREDISTRIBUTE_HH
4 #define DUNE_ISTL_MATRIXREDISTRIBUTE_HH
5 #include <memory>
6 #include "repartition.hh"
7 #include <dune/common/exceptions.hh>
8 #include <dune/common/parallel/indexset.hh>
9 #include <dune/common/unused.hh>
11 #include <dune/istl/paamg/pinfo.hh>
17 namespace Dune
18 {
19  template<typename T>
21  {
22  bool isSetup() const
23  {
24  return false;
25  }
26  template<class D>
27  void redistribute(const D& from, D& to) const
28  {
29  DUNE_UNUSED_PARAMETER(from);
30  DUNE_UNUSED_PARAMETER(to);
31  }
32 
33  template<class D>
34  void redistributeBackward(D& from, const D& to) const
35  {
36  DUNE_UNUSED_PARAMETER(from);
37  DUNE_UNUSED_PARAMETER(to);
38  }
39 
40  void resetSetup()
41  {}
42 
43  void setNoRows(std::size_t size)
44  {
45  DUNE_UNUSED_PARAMETER(size);
46  }
47 
48  void setNoCopyRows(std::size_t size)
49  {
50  DUNE_UNUSED_PARAMETER(size);
51  }
52 
53  void setNoBackwardsCopyRows(std::size_t size)
54  {
55  DUNE_UNUSED_PARAMETER(size);
56  }
57 
58  std::size_t getRowSize(std::size_t index) const
59  {
60  DUNE_UNUSED_PARAMETER(index);
61  return -1;
62  }
63 
64  std::size_t getCopyRowSize(std::size_t index) const
65  {
66  DUNE_UNUSED_PARAMETER(index);
67  return -1;
68  }
69 
70  std::size_t getBackwardsCopyRowSize(std::size_t index) const
71  {
72  DUNE_UNUSED_PARAMETER(index);
73  return -1;
74  }
75 
76  };
77 
78 #if HAVE_MPI
79  template<typename T, typename T1>
81  {
82  public:
84 
86  : interface(), setup_(false)
87  {}
88 
90  {
91  return interface;
92  }
93  template<typename IS>
94  void checkInterface(const IS& source,
95  const IS& target, MPI_Comm comm)
96  {
97  auto ri = std::make_unique<RemoteIndices<IS> >(source, target, comm);
98  ri->template rebuild<true>();
99  Interface inf;
101  int rank;
102  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
103  inf.free();
104  inf.build(*ri, flags, flags);
105 
106 
107 #ifdef DEBUG_REPART
108  if(inf!=interface) {
109 
110  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
111  if(rank==0)
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;
115 
116  throw "autsch!";
117  }
118 #endif
119  }
120  void setSetup()
121  {
122  setup_=true;
123  interface.strip();
124  }
125 
126  void resetSetup()
127  {
128  setup_=false;
129  }
130 
131  template<class GatherScatter, class D>
132  void redistribute(const D& from, D& to) const
133  {
134  BufferedCommunicator communicator;
135  communicator.template build<D>(from,to, interface);
136  communicator.template forward<GatherScatter>(from, to);
137  communicator.free();
138  }
139  template<class GatherScatter, class D>
140  void redistributeBackward(D& from, const D& to) const
141  {
142 
143  BufferedCommunicator communicator;
144  communicator.template build<D>(from,to, interface);
145  communicator.template backward<GatherScatter>(from, to);
146  communicator.free();
147  }
148 
149  template<class D>
150  void redistribute(const D& from, D& to) const
151  {
152  redistribute<CopyGatherScatter<D> >(from,to);
153  }
154  template<class D>
155  void redistributeBackward(D& from, const D& to) const
156  {
157  redistributeBackward<CopyGatherScatter<D> >(from,to);
158  }
159  bool isSetup() const
160  {
161  return setup_;
162  }
163 
164  void reserve(std::size_t size)
165  {}
166 
167  std::size_t& getRowSize(std::size_t index)
168  {
169  return rowSize[index];
170  }
171 
172  std::size_t getRowSize(std::size_t index) const
173  {
174  return rowSize[index];
175  }
176 
177  std::size_t& getCopyRowSize(std::size_t index)
178  {
179  return copyrowSize[index];
180  }
181 
182  std::size_t getCopyRowSize(std::size_t index) const
183  {
184  return copyrowSize[index];
185  }
186 
187  std::size_t& getBackwardsCopyRowSize(std::size_t index)
188  {
189  return backwardscopyrowSize[index];
190  }
191 
192  std::size_t getBackwardsCopyRowSize(std::size_t index) const
193  {
194  return backwardscopyrowSize[index];
195  }
196 
197  void setNoRows(std::size_t rows)
198  {
199  rowSize.resize(rows, 0);
200  }
201 
202  void setNoCopyRows(std::size_t rows)
203  {
204  copyrowSize.resize(rows, 0);
205  }
206 
207  void setNoBackwardsCopyRows(std::size_t rows)
208  {
209  backwardscopyrowSize.resize(rows, 0);
210  }
211 
212  private:
213  std::vector<std::size_t> rowSize;
214  std::vector<std::size_t> copyrowSize;
215  std::vector<std::size_t> backwardscopyrowSize;
216  RedistributeInterface interface;
217  bool setup_;
218  };
219 
228  template<class M, class RI>
230  {
231  // Make the default communication policy work.
232  typedef typename M::size_type value_type;
233  typedef typename M::size_type size_type;
234 
240  CommMatrixRowSize(const M& m_, RI& rowsize_)
241  : matrix(m_), rowsize(rowsize_)
242  {}
243  const M& matrix;
244  RI& rowsize;
245 
246  };
247 
248 
257  template<class M, class I>
259  {
260  typedef typename M::size_type size_type;
261 
268  CommMatrixSparsityPattern(const M& m_, const Dune::GlobalLookupIndexSet<I>& idxset_, const I& aggidxset_)
269  : matrix(m_), idxset(idxset_), aggidxset(aggidxset_), rowsize()
270  {}
271 
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_)
282  {}
283 
291  {
292  // insert diagonal to overlap rows
293  typedef typename Dune::GlobalLookupIndexSet<I>::const_iterator IIter;
294  typedef typename Dune::OwnerOverlapCopyCommunication<int>::OwnerSet OwnerSet;
295  std::size_t nnz=0;
296 #ifdef DEBUG_REPART
297  int rank;
298 
299  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
300 #endif
301  for(IIter i= aggidxset.begin(), end=aggidxset.end(); i!=end; ++i) {
302  if(!OwnerSet::contains(i->local().attribute())) {
303 #ifdef DEBUG_REPART
304  std::cout<<rank<<" Inserting diagonal for"<<i->local()<<std::endl;
305 #endif
306  sparsity[i->local()].insert(i->local());
307  }
308 
309  nnz+=sparsity[i->local()].size();
310  }
311  assert( aggidxset.size()==sparsity.size());
312 
313  if(nnz>0) {
314  m.setSize(aggidxset.size(), aggidxset.size(), nnz);
315  m.setBuildMode(M::row_wise);
316  typename M::CreateIterator citer=m.createbegin();
317 #ifdef DEBUG_REPART
318  std::size_t idx=0;
319  bool correct=true;
320  Dune::GlobalLookupIndexSet<I> global(aggidxset);
321 #endif
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)
324  {
325  typedef typename std::set<size_type>::const_iterator SIter;
326  for(SIter si=i->begin(), send=i->end(); si!=send; ++si)
327  citer.insert(*si);
328 #ifdef DEBUG_REPART
329  if(i->find(idx)==i->end()) {
330  const typename I::IndexPair* gi=global.pair(idx);
331  assert(gi);
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;
335  correct=false;
336  }
337  ++idx;
338 #endif
339  }
340 #ifdef DEBUG_REPART
341  if(!correct)
342  throw "bla";
343 #endif
344  }
345  }
346 
354  void completeSparsityPattern(std::vector<std::set<size_type> > add_sparsity)
355  {
356  for (unsigned int i = 0; i != sparsity.size(); ++i) {
357  if (add_sparsity[i].size() != 0) {
358  typedef std::set<size_type> Set;
359  Set tmp_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);
364  }
365  }
366  }
367 
368  const M& matrix;
369  typedef Dune::GlobalLookupIndexSet<I> LookupIndexSet;
370  const Dune::GlobalLookupIndexSet<I>& idxset;
371  const I& aggidxset;
372  std::vector<std::set<size_type> > sparsity;
373  const std::vector<size_type>* rowsize;
374  };
375 
376  template<class M, class I>
377  struct CommPolicy<CommMatrixSparsityPattern<M,I> >
378  {
380 
385  typedef typename I::GlobalIndex IndexedType;
386 
388  typedef VariableSize IndexedTypeFlag;
389 
390  static typename M::size_type getSize(const Type& t, std::size_t i)
391  {
392  if(!t.rowsize)
393  return t.matrix[i].size();
394  else
395  {
396  assert((*t.rowsize)[i]>0);
397  return (*t.rowsize)[i];
398  }
399  }
400  };
401 
408  template<class M, class I>
410  {
419  CommMatrixRow(M& m_, const Dune::GlobalLookupIndexSet<I>& idxset_, const I& aggidxset_)
420  : matrix(m_), idxset(idxset_), aggidxset(aggidxset_), rowsize()
421  {}
422 
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_)
429  {}
436  {
437  typedef typename Dune::GlobalLookupIndexSet<I>::const_iterator Iter;
438  typedef typename Dune::OwnerOverlapCopyCommunication<int>::OwnerSet OwnerSet;
439 
440  for(Iter i= aggidxset.begin(), end=aggidxset.end(); i!=end; ++i)
441  if(!OwnerSet::contains(i->local().attribute())) {
442  // Set to Dirchlet
443  typedef typename M::ColIterator CIter;
444  for(CIter c=matrix[i->local()].begin(), cend= matrix[i->local()].end();
445  c!= cend; ++c)
446  {
447  *c=0;
448  if(c.index()==i->local()) {
449  typedef typename M::block_type::RowIterator RIter;
450  for(RIter r=c->begin(), rend=c->end();
451  r != rend; ++r)
452  (*r)[r.index()]=1;
453  }
454  }
455  }
456  }
458  M& matrix;
460  const Dune::GlobalLookupIndexSet<I>& idxset;
462  const I& aggidxset;
464  std::vector<size_t>* rowsize; // row sizes differ from sender side in overlap!
465  };
466 
467  template<class M, class I>
468  struct CommPolicy<CommMatrixRow<M,I> >
469  {
471 
476  typedef std::pair<typename I::GlobalIndex,typename M::block_type> IndexedType;
477 
479  typedef VariableSize IndexedTypeFlag;
480 
481  static std::size_t getSize(const Type& t, std::size_t i)
482  {
483  if(!t.rowsize)
484  return t.matrix[i].size();
485  else
486  {
487  assert((*t.rowsize)[i]>0);
488  return (*t.rowsize)[i];
489  }
490  }
491  };
492 
493  template<class M, class I, class RI>
495  {
497 
498  static const typename M::size_type gather(const Container& cont, std::size_t i)
499  {
500  return cont.matrix[i].size();
501  }
502  static void scatter(Container& cont, const typename M::size_type& rowsize, std::size_t i)
503  {
504  assert(rowsize);
505  cont.rowsize.getRowSize(i)=rowsize;
506  }
507 
508  };
509 
510  template<class M, class I, class RI>
512  {
514 
515  static const typename M::size_type gather(const Container& cont, std::size_t i)
516  {
517  return cont.matrix[i].size();
518  }
519  static void scatter(Container& cont, const typename M::size_type& rowsize, std::size_t i)
520  {
521  assert(rowsize);
522  if (rowsize > cont.rowsize.getCopyRowSize(i))
523  cont.rowsize.getCopyRowSize(i)=rowsize;
524  }
525 
526  };
527 
528  template<class M, class I>
530  {
531  typedef typename I::GlobalIndex GlobalIndex;
533  typedef typename M::ConstColIterator ColIter;
534 
535  static ColIter col;
536  static GlobalIndex numlimits;
537 
538  static const GlobalIndex& gather(const Container& cont, std::size_t i, std::size_t j)
539  {
540  if(j==0)
541  col=cont.matrix[i].begin();
542  else if (col!=cont.matrix[i].end())
543  ++col;
544 
545  //copy communication: different row sizes for copy rows with the same global index
546  //are possible. If all values of current matrix row are sent, send
547  //std::numeric_limits<GlobalIndex>::max()
548  //and receiver will ignore it.
549  if (col==cont.matrix[i].end()) {
550  numlimits = std::numeric_limits<GlobalIndex>::max();
551  return numlimits;
552  }
553  else {
554  const typename I::IndexPair* index=cont.idxset.pair(col.index());
555  assert(index);
556  // Only send index if col is no ghost
557  if ( index->local().attribute() != 2)
558  return index->global();
559  else {
560  numlimits = std::numeric_limits<GlobalIndex>::max();
561  return numlimits;
562  }
563  }
564  }
565  static void scatter(Container& cont, const GlobalIndex& gi, std::size_t i, std::size_t j)
566  {
567  DUNE_UNUSED_PARAMETER(j);
568  try{
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);
574 
575  typedef typename Dune::OwnerOverlapCopyCommunication<int>::OwnerSet OwnerSet;
576  if(!OwnerSet::contains(ip.local().attribute()))
577  // preserve symmetry for overlap
578  cont.sparsity[column].insert(i);
579  }
580  }
581  catch(Dune::RangeError er) {
582  // Entry not present in the new index set. Ignore!
583 #ifdef DEBUG_REPART
584  typedef typename Container::LookupIndexSet GlobalLookup;
585  typedef typename GlobalLookup::IndexPair IndexPair;
586  typedef typename Dune::OwnerOverlapCopyCommunication<int>::OwnerSet OwnerSet;
587 
588  GlobalLookup lookup(cont.aggidxset);
589  const IndexPair* pi=lookup.pair(i);
590  assert(pi);
591  if(OwnerSet::contains(pi->local().attribute())) {
592  int rank;
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;
596  throw er;
597  }
598 #endif
599  }
600  }
601 
602  };
603  template<class M, class I>
605 
606  template<class M, class I>
608 
609 
610  template<class M, class I>
612  {
613  typedef typename I::GlobalIndex GlobalIndex;
615  typedef typename M::ConstColIterator ColIter;
616  typedef typename std::pair<GlobalIndex,typename M::block_type> Data;
617  static ColIter col;
618  static Data datastore;
619  static GlobalIndex numlimits;
620 
621  static const Data& gather(const Container& cont, std::size_t i, std::size_t j)
622  {
623  if(j==0)
624  col=cont.matrix[i].begin();
625  else if (col!=cont.matrix[i].end())
626  ++col;
627  // copy communication: different row sizes for copy rows with the same global index
628  // are possible. If all values of current matrix row are sent, send
629  // std::numeric_limits<GlobalIndex>::max()
630  // and receiver will ignore it.
631  if (col==cont.matrix[i].end()) {
632  numlimits = std::numeric_limits<GlobalIndex>::max();
633  datastore = std::make_pair(numlimits,*col);
634  return datastore;
635  }
636  else {
637  // convert local column index to global index
638  const typename I::IndexPair* index=cont.idxset.pair(col.index());
639  assert(index);
640  // Store the data to prevent reference to temporary
641  // Only send index if col is no ghost
642  if ( index->local().attribute() != 2)
643  datastore = std::make_pair(index->global(),*col);
644  else {
645  numlimits = std::numeric_limits<GlobalIndex>::max();
646  datastore = std::make_pair(numlimits,*col);
647  }
648  return datastore;
649  }
650  }
651  static void scatter(Container& cont, const Data& data, std::size_t i, std::size_t j)
652  {
653  DUNE_UNUSED_PARAMETER(j);
654  try{
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;
658  }
659  }
660  catch(Dune::RangeError er) {
661  // This an overlap row and might therefore lack some entries!
662  }
663 
664  }
665  };
666 
667  template<class M, class I>
669 
670  template<class M, class I>
672 
673  template<class M, class I>
675 
676  template<typename M, typename C>
677  void redistributeSparsityPattern(M& origMatrix, M& newMatrix, C& origComm, C& newComm,
679  {
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);
687 
688  // get owner rowsizes
689  CommMatrixRowSize<M,RI> commRowSize(origMatrix, ri);
690  ri.template redistribute<MatrixRowSizeGatherScatter<M,IndexSet,RI> >(commRowSize,commRowSize);
691 
692  origComm.buildGlobalLookup();
693 
694  for (std::size_t i=0; i < newComm.indexSet().size(); i++) {
695  rowsize[i] = ri.getRowSize(i);
696  }
697  // get sparsity pattern from owner rows
699  origsp(origMatrix, origComm.globalLookup(), newComm.indexSet());
701  newsp(origMatrix, origComm.globalLookup(), newComm.indexSet(), rowsize);
702 
703  ri.template redistribute<MatrixSparsityPatternGatherScatter<M,IndexSet> >(origsp,newsp);
704 
705  // build copy to owner interface to get missing matrix values for novlp case
706  if (origComm.getSolverCategory() == SolverCategory::nonoverlapping) {
707  RemoteIndices<IndexSet> *ris = new RemoteIndices<IndexSet>(origComm.indexSet(),
708  newComm.indexSet(),
709  origComm.communicator());
710  ris->template rebuild<true>();
711 
712  ri.getInterface().free();
713  ri.getInterface().build(*ris,copyflags,ownerflags);
714 
715  // get copy rowsizes
716  CommMatrixRowSize<M,RI> commRowSize_copy(origMatrix, ri);
717  ri.template redistribute<MatrixCopyRowSizeGatherScatter<M,IndexSet,RI> >(commRowSize_copy,
718  commRowSize_copy);
719 
720  for (std::size_t i=0; i < newComm.indexSet().size(); i++) {
721  copyrowsize[i] = ri.getCopyRowSize(i);
722  }
723  //get copy rowsizes for sender
724  ri.redistributeBackward(backwardscopyrowsize,copyrowsize);
725  for (std::size_t i=0; i < origComm.indexSet().size(); i++) {
726  ri.getBackwardsCopyRowSize(i) = backwardscopyrowsize[i];
727  }
728 
729  // get sparsity pattern from copy rows
730  CommMatrixSparsityPattern<M,IndexSet> origsp_copy(origMatrix,
731  origComm.globalLookup(),
732  newComm.indexSet(),
733  backwardscopyrowsize);
734  CommMatrixSparsityPattern<M,IndexSet> newsp_copy(origMatrix, origComm.globalLookup(),
735  newComm.indexSet(), copyrowsize);
736  ri.template redistribute<MatrixSparsityPatternGatherScatter<M,IndexSet> >(origsp_copy,
737  newsp_copy);
738 
739  newsp.completeSparsityPattern(newsp_copy.sparsity);
740  newsp.storeSparsityPattern(newMatrix);
741  }
742  else
743  newsp.storeSparsityPattern(newMatrix);
744 
745 #ifdef DUNE_ISTL_WITH_CHECKING
746  // Check for symmetry
747  int ret=0;
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)
752  {
753  try{
754  newMatrix[col.index()][row.index()];
755  }catch(Dune::ISTLError e) {
756  std::cerr<<newComm.communicator().rank()<<": entry ("
757  <<col.index()<<","<<row.index()<<") missing! for symmetry!"<<std::endl;
758  ret=1;
759 
760  }
761 
762  }
763  }
764 
765  if(ret)
766  DUNE_THROW(ISTLError, "Matrix not symmetric!");
767 #endif
768  }
769 
770  template<typename M, typename C>
771  void redistributeMatrixEntries(M& origMatrix, M& newMatrix, C& origComm, C& newComm,
773  {
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);
779 
780  for (std::size_t i=0; i < newComm.indexSet().size(); i++) {
781  rowsize[i] = ri.getRowSize(i);
782  if (origComm.getSolverCategory() == SolverCategory::nonoverlapping) {
783  copyrowsize[i] = ri.getCopyRowSize(i);
784  }
785  }
786 
787  for (std::size_t i=0; i < origComm.indexSet().size(); i++)
788  if (origComm.getSolverCategory() == SolverCategory::nonoverlapping)
789  backwardscopyrowsize[i] = ri.getBackwardsCopyRowSize(i);
790 
791 
792  if (origComm.getSolverCategory() == SolverCategory::nonoverlapping) {
793  // fill sparsity pattern from copy rows
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,
799  newrow_copy);
800  ri.getInterface().free();
801  RemoteIndices<IndexSet> *ris = new RemoteIndices<IndexSet>(origComm.indexSet(),
802  newComm.indexSet(),
803  origComm.communicator());
804  ris->template rebuild<true>();
805  ri.getInterface().build(*ris,ownerflags,ownerflags);
806  }
807 
809  origrow(origMatrix, origComm.globalLookup(), newComm.indexSet());
811  newrow(newMatrix, origComm.globalLookup(), newComm.indexSet(),rowsize);
812  ri.template redistribute<MatrixRowGatherScatter<M,IndexSet> >(origrow,newrow);
813  if (origComm.getSolverCategory() != static_cast<int>(SolverCategory::nonoverlapping))
814  newrow.setOverlapRowsToDirichlet();
815  }
816 
833  template<typename M, typename C>
834  void redistributeMatrix(M& origMatrix, M& newMatrix, C& origComm, C& newComm,
836  {
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);
842  }
843 #endif
844 
845 template<typename M>
846  void redistributeMatrixEntries(M& origMatrix, M& newMatrix,
850  {
851  DUNE_THROW(InvalidStateException, "Trying to redistribute in sequential program!");
852  }
853  template<typename M>
854  void redistributeMatrix(M& origMatrix, M& newMatrix,
858  {
859  DUNE_THROW(InvalidStateException, "Trying to redistribute in sequential program!");
860  }
861 }
862 #endif
void reserve(std::size_t size)
Definition: matrixredistribute.hh:164
Dune::GlobalLookupIndexSet< I > LookupIndexSet
Definition: matrixredistribute.hh:369
std::size_t getBackwardsCopyRowSize(std::size_t index) const
Definition: matrixredistribute.hh:192
CommMatrixRowSize< M, RI > Container
Definition: matrixredistribute.hh:496
void setNoRows(std::size_t size)
Definition: matrixredistribute.hh:43
void setNoRows(std::size_t rows)
Definition: matrixredistribute.hh:197
Definition: matrixredistribute.hh:494
void redistribute(const D &from, D &to) const
Definition: matrixredistribute.hh:132
CommMatrixRow< M, I > Type
Definition: matrixredistribute.hh:470
EnumItem< AttributeSet, OwnerOverlapCopyAttributeSet::owner > OwnerSet
Definition: owneroverlapcopy.hh:192
Definition: allocator.hh:7
void storeSparsityPattern(M &m)
Creates and stores the sparsity pattern of the redistributed matrix.
Definition: matrixredistribute.hh:290
RedistributeInterface & getInterface()
Definition: matrixredistribute.hh:89
std::vector< size_t > * rowsize
row size information for the receiving side.
Definition: matrixredistribute.hh:464
std::pair< GlobalIndex, typename M::block_type > Data
Definition: matrixredistribute.hh:616
const std::vector< size_type > * rowsize
Definition: matrixredistribute.hh:373
Definition: matrixredistribute.hh:529
Functionality for redistributing a parallel index set using graph partitioning.
I::GlobalIndex GlobalIndex
Definition: matrixredistribute.hh:613
std::size_t & getCopyRowSize(std::size_t index)
Definition: matrixredistribute.hh:177
void redistribute(const D &from, D &to) const
Definition: matrixredistribute.hh:150
void redistributeSparsityPattern(M &origMatrix, M &newMatrix, C &origComm, C &newComm, RedistributeInformation< C > &ri)
Definition: matrixredistribute.hh:677
static const M::size_type gather(const Container &cont, std::size_t i)
Definition: matrixredistribute.hh:498
M::size_type size_type
Definition: matrixredistribute.hh:233
RI & rowsize
Definition: matrixredistribute.hh:244
CommMatrixRow< M, I > Container
Definition: matrixredistribute.hh:614
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:476
static M::size_type getSize(const Type &t, std::size_t i)
Definition: matrixredistribute.hh:390
static void scatter(Container &cont, const GlobalIndex &gi, std::size_t i, std::size_t j)
Definition: matrixredistribute.hh:565
static void scatter(Container &cont, const typename M::size_type &rowsize, std::size_t i)
Definition: matrixredistribute.hh:519
I::GlobalIndex GlobalIndex
Definition: matrixredistribute.hh:531
VariableSize IndexedTypeFlag
Each row varies in size.
Definition: matrixredistribute.hh:479
static ColIter col
Definition: matrixredistribute.hh:617
CommMatrixSparsityPattern(const M &m_, const Dune::GlobalLookupIndexSet< I > &idxset_, const I &aggidxset_)
Constructor for the original side.
Definition: matrixredistribute.hh:268
Category for non-overlapping solvers.
Definition: solvercategory.hh:25
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
CommMatrixRowSize(const M &m_, RI &rowsize_)
Constructor.
Definition: matrixredistribute.hh:240
M::size_type size_type
Definition: matrixredistribute.hh:260
std::size_t getBackwardsCopyRowSize(std::size_t index) const
Definition: matrixredistribute.hh:70
static GlobalIndex numlimits
Definition: matrixredistribute.hh:619
static const M::size_type gather(const Container &cont, std::size_t i)
Definition: matrixredistribute.hh:515
CommMatrixRow(M &m_, const Dune::GlobalLookupIndexSet< I > &idxset_, const I &aggidxset_, std::vector< size_t > &rowsize_)
Constructor.
Definition: matrixredistribute.hh:426
const I & aggidxset
Index set for the redistributed matrix.
Definition: matrixredistribute.hh:462
M::ConstColIterator ColIter
Definition: matrixredistribute.hh:533
M & matrix
The matrix to communicate the values of.
Definition: matrixredistribute.hh:458
bool isSetup() const
Definition: matrixredistribute.hh:159
I::GlobalIndex IndexedType
The indexed type we send. This is the global index indentitfying the column.
Definition: matrixredistribute.hh:385
Utility class to communicate and set the row sizes of a redistributed matrix.
Definition: matrixredistribute.hh:229
static void scatter(Container &cont, const Data &data, std::size_t i, std::size_t j)
Definition: matrixredistribute.hh:651
static void scatter(Container &cont, const typename M::size_type &rowsize, std::size_t i)
Definition: matrixredistribute.hh:502
const M & matrix
Definition: matrixredistribute.hh:368
Definition: matrixredistribute.hh:511
const I & aggidxset
Definition: matrixredistribute.hh:371
static std::size_t getSize(const Type &t, std::size_t i)
Definition: matrixredistribute.hh:481
std::size_t getCopyRowSize(std::size_t index) const
Definition: matrixredistribute.hh:64
static Data datastore
Definition: matrixredistribute.hh:618
void checkInterface(const IS &source, const IS &target, MPI_Comm comm)
Definition: matrixredistribute.hh:94
CommMatrixSparsityPattern< M, I > Container
Definition: matrixredistribute.hh:532
void redistribute(const D &from, D &to) const
Definition: matrixredistribute.hh:27
void redistributeBackward(D &from, const D &to) const
Definition: matrixredistribute.hh:140
CommMatrixRow(M &m_, const Dune::GlobalLookupIndexSet< I > &idxset_, const I &aggidxset_)
Constructor.
Definition: matrixredistribute.hh:419
void resetSetup()
Definition: matrixredistribute.hh:40
M::size_type value_type
Definition: matrixredistribute.hh:232
static const Data & gather(const Container &cont, std::size_t i, std::size_t j)
Definition: matrixredistribute.hh:621
const Dune::GlobalLookupIndexSet< I > & idxset
Index set for the original matrix.
Definition: matrixredistribute.hh:460
derive error class from the base class in common
Definition: istlexception.hh:16
const M & matrix
Definition: matrixredistribute.hh:243
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
void setNoCopyRows(std::size_t size)
Definition: matrixredistribute.hh:48
Definition: matrixredistribute.hh:611
Definition: pinfo.hh:25
std::vector< std::set< size_type > > sparsity
Definition: matrixredistribute.hh:372
std::size_t getCopyRowSize(std::size_t index) const
Definition: matrixredistribute.hh:182
void redistributeBackward(D &from, const D &to) const
Definition: matrixredistribute.hh:34
void setNoBackwardsCopyRows(std::size_t rows)
Definition: matrixredistribute.hh:207
bool isSetup() const
Definition: matrixredistribute.hh:22
CommMatrixSparsityPattern< M, I > Type
Definition: matrixredistribute.hh:379
std::size_t & getBackwardsCopyRowSize(std::size_t index)
Definition: matrixredistribute.hh:187
CommMatrixRowSize< M, RI > Container
Definition: matrixredistribute.hh:513
Definition: repartition.hh:265
void redistributeBackward(D &from, const D &to) const
Definition: matrixredistribute.hh:155
std::size_t getRowSize(std::size_t index) const
Definition: matrixredistribute.hh:172
VariableSize IndexedTypeFlag
Each row varies in size.
Definition: matrixredistribute.hh:388
Col col
Definition: matrixmatrix.hh:349
A class setting up standard communication for a two-valued attribute set with owner/overlap/copy sema...
Definition: owneroverlapcopy.hh:171
Classes providing communication interfaces for overlapping Schwarz methods.
const Dune::GlobalLookupIndexSet< I > & idxset
Definition: matrixredistribute.hh:370
std::size_t getRowSize(std::size_t index) const
Definition: matrixredistribute.hh:58
Utility class to communicate and build the sparsity pattern of a redistributed matrix.
Definition: matrixredistribute.hh:258
static GlobalIndex numlimits
Definition: matrixredistribute.hh:536
Definition: matrixredistribute.hh:20
void setOverlapRowsToDirichlet()
Sets the non-owner rows correctly as Dirichlet boundaries.
Definition: matrixredistribute.hh:435
void setNoCopyRows(std::size_t rows)
Definition: matrixredistribute.hh:202
void setNoBackwardsCopyRows(std::size_t size)
Definition: matrixredistribute.hh:53
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
static ColIter col
Definition: matrixredistribute.hh:535
M::ConstColIterator ColIter
Definition: matrixredistribute.hh:615
static const GlobalIndex & gather(const Container &cont, std::size_t i, std::size_t j)
Definition: matrixredistribute.hh:538
void redistributeMatrixEntries(M &origMatrix, M &newMatrix, C &origComm, C &newComm, RedistributeInformation< C > &ri)
Definition: matrixredistribute.hh:771
Utility class for comunicating the matrix entries.
Definition: matrixredistribute.hh:409
std::size_t & getRowSize(std::size_t index)
Definition: matrixredistribute.hh:167
OwnerOverlapCopyCommunication< T, T1 > Comm
Definition: matrixredistribute.hh:83