dune-istl  2.2.0
matrixredistribute.hh
Go to the documentation of this file.
1 // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set ts=8 sw=2 et sts=2:
3 #ifndef DUNE_MATRIXREDIST_HH
4 #define DUNE_MATRIXREDIST_HH
5 #include"repartition.hh"
6 #include<dune/common/exceptions.hh>
7 #include<dune/common/parallel/indexset.hh>
14 namespace Dune
15 {
16  template<typename T>
18  {
19  bool isSetup() const
20  {
21  return false;
22  }
23  template<class D>
24  void redistribute(const D& from, D& to) const
25  {}
26 
27  template<class D>
28  void redistributeBackward(D& from, const D& to) const
29  {}
30 
31  void resetSetup()
32  {}
33 
34  void setNoRows(std::size_t size)
35  {}
36 
37  void setNoCopyRows(std::size_t size)
38  {}
39 
40  void setNoBackwardsCopyRows(std::size_t size)
41  {}
42 
43  std::size_t getRowSize(std::size_t index) const
44  {
45  return -1;
46  }
47 
48  std::size_t getCopyRowSize(std::size_t index) const
49  {
50  return -1;
51  }
52 
53  std::size_t getBackwardsCopyRowSize(std::size_t index) const
54  {
55  return -1;
56  }
57 
58  };
59 
60 #if HAVE_MPI
61  template<typename T, typename T1>
63  {
64  public:
66 
68  : interface(), setup_(false)
69  {}
70 
71  RedistributeInterface& getInterface()
72  {
73  return interface;
74  }
75  template<typename IS>
76  void checkInterface(const IS& source,
77  const IS& target, MPI_Comm comm)
78  {
79  RemoteIndices<IS> *ri=new RemoteIndices<IS>(source, target, comm);
80  ri->template rebuild<true>();
81  Interface inf;
83  int rank;
84  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
85  inf.free();
86  inf.build(*ri, flags, flags);
87 
88 
89 #ifdef DEBUG_REPART
90  if(inf!=interface){
91 
92  int rank;
93  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
94  if(rank==0)
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;
98 
99  throw "autsch!";
100  delete ri;
101  }else
102 
103 #endif
104  delete ri;
105  }
106  void setSetup()
107  {
108  setup_=true;
109  interface.strip();
110  }
111 
112  void resetSetup()
113  {
114  setup_=false;
115  }
116 
117  template<class GatherScatter, class D>
118  void redistribute(const D& from, D& to) const
119  {
120  BufferedCommunicator communicator;
121  communicator.template build<D>(from,to, interface);
122  communicator.template forward<GatherScatter>(from, to);
123  communicator.free();
124  }
125  template<class GatherScatter, class D>
126  void redistributeBackward(D& from, const D& to) const
127  {
128 
129  BufferedCommunicator communicator;
130  communicator.template build<D>(from,to, interface);
131  communicator.template backward<GatherScatter>(from, to);
132  communicator.free();
133  }
134 
135  template<class D>
136  void redistribute(const D& from, D& to) const
137  {
138  redistribute<CopyGatherScatter<D> >(from,to);
139  }
140  template<class D>
141  void redistributeBackward(D& from, const D& to) const
142  {
143  redistributeBackward<CopyGatherScatter<D> >(from,to);
144  }
145  bool isSetup() const
146  {
147  return setup_;
148  }
149 
150  void reserve(std::size_t size)
151  {}
152 
153  std::size_t& getRowSize(std::size_t index)
154  {
155  return rowSize[index];
156  }
157 
158  std::size_t getRowSize(std::size_t index) const
159  {
160  return rowSize[index];
161  }
162 
163  std::size_t& getCopyRowSize(std::size_t index)
164  {
165  return copyrowSize[index];
166  }
167 
168  std::size_t getCopyRowSize(std::size_t index) const
169  {
170  return copyrowSize[index];
171  }
172 
173  std::size_t& getBackwardsCopyRowSize(std::size_t index)
174  {
175  return backwardscopyrowSize[index];
176  }
177 
178  std::size_t getBackwardsCopyRowSize(std::size_t index) const
179  {
180  return backwardscopyrowSize[index];
181  }
182 
183  void setNoRows(std::size_t rows)
184  {
185  rowSize.resize(rows, 0);
186  }
187 
188  void setNoCopyRows(std::size_t rows)
189  {
190  copyrowSize.resize(rows, 0);
191  }
192 
193  void setNoBackwardsCopyRows(std::size_t rows)
194  {
195  backwardscopyrowSize.resize(rows, 0);
196  }
197 
198  private:
199  std::vector<std::size_t> rowSize;
200  std::vector<std::size_t> copyrowSize;
201  std::vector<std::size_t> backwardscopyrowSize;
202  RedistributeInterface interface;
203  bool setup_;
204  };
205 
214  template<class M, class RI>
216  {
217  // Make the default communication policy work.
218  typedef typename M::size_type value_type;
219  typedef typename M::size_type size_type;
220 
226  CommMatrixRowSize(const M& m_, RI& rowsize_)
227  : matrix(m_), rowsize(rowsize_)
228  {}
229  const M& matrix;
230  RI& rowsize;
231 
232  };
233 
234 
243  template<class M, class I>
245  {
246  typedef typename M::size_type size_type;
247 
254  CommMatrixSparsityPattern(const M& m_, const Dune::GlobalLookupIndexSet<I>& idxset_, const I& aggidxset_)
255  : matrix(m_), idxset(idxset_), aggidxset(aggidxset_), rowsize()
256  {}
257 
265  CommMatrixSparsityPattern(const M& m_, const Dune::GlobalLookupIndexSet<I>& idxset_, const I& aggidxset_,
266  const std::vector<typename M::size_type>& rowsize_)
267  : matrix(m_), idxset(idxset_), aggidxset(aggidxset_), sparsity(aggidxset_.size()), rowsize(&rowsize_)
268  {}
269 
277  {
278  // insert diagonal to overlap rows
279  typedef typename Dune::GlobalLookupIndexSet<I>::const_iterator IIter;
280  typedef typename Dune::OwnerOverlapCopyCommunication<int>::OwnerSet OwnerSet;
281  std::size_t nnz=0;
282 #ifdef DEBUG_REPART
283  int rank;
284 
285  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
286 #endif
287  for(IIter i= aggidxset.begin(), end=aggidxset.end(); i!=end; ++i){
288  if(!OwnerSet::contains(i->local().attribute())){
289 #ifdef DEBUG_REPART
290  std::cout<<rank<<" Inserting diagonal for"<<i->local()<<std::endl;
291 #endif
292  sparsity[i->local()].insert(i->local());
293  }
294 
295  nnz+=sparsity[i->local()].size();
296  }
297  assert( aggidxset.size()==sparsity.size());
298 
299  if(nnz>0){
300  m.setSize(aggidxset.size(), aggidxset.size(), nnz);
301  m.setBuildMode(M::row_wise);
302  typename M::CreateIterator citer=m.createbegin();
303 #ifdef DEBUG_REPART
304  std::size_t idx=0;
305  bool correct=true;
306  Dune::GlobalLookupIndexSet<I> global(aggidxset);
307 #endif
308  typedef typename std::vector<std::set<size_type> >::const_iterator Iter;
309  for(Iter i=sparsity.begin(), end=sparsity.end(); i!=end; ++i, ++citer)
310  {
311  typedef typename std::set<size_type>::const_iterator SIter;
312  for(SIter si=i->begin(), send=i->end(); si!=send; ++si)
313  citer.insert(*si);
314 #ifdef DEBUG_REPART
315  if(i->find(idx)==i->end()){
316  const typename I::IndexPair* gi=global.pair(idx);
317  assert(gi);
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;
321  correct=false;
322  }
323  ++idx;
324 #endif
325  }
326 #ifdef DEBUG_REPART
327  if(!correct)
328  throw "bla";
329 #endif
330  }
331  }
332 
340  void completeSparsityPattern(std::vector<std::set<size_type> > add_sparsity)
341  {
342  typedef typename std::vector<std::set<size_type> >::const_iterator Iter;
343  for (unsigned int i = 0; i != sparsity.size(); ++i) {
344  if (add_sparsity[i].size() != 0){
345  typedef std::set<size_type> Set;
346  Set tmp_set;
347  std::insert_iterator<Set> tmp_insert (tmp_set, tmp_set.begin());
348  std::set_union(add_sparsity[i].begin(), add_sparsity[i].end(),
349  sparsity[i].begin(), sparsity[i].end(), tmp_insert);
350  sparsity[i].swap(tmp_set);
351  }
352  }
353  }
354 
355  const M& matrix;
356  typedef Dune::GlobalLookupIndexSet<I> LookupIndexSet;
357  const Dune::GlobalLookupIndexSet<I>& idxset;
358  const I& aggidxset;
359  std::vector<std::set<size_type> > sparsity;
360  const std::vector<size_type>* rowsize;
361  };
362 
363  template<class M, class I>
364  struct CommPolicy<CommMatrixSparsityPattern<M,I> >
365  {
367 
372  typedef typename I::GlobalIndex IndexedType;
373 
375  typedef VariableSize IndexedTypeFlag;
376 
377  static typename M::size_type getSize(const Type& t, std::size_t i)
378  {
379  if(!t.rowsize)
380  return t.matrix[i].size();
381  else
382  {
383  assert((*t.rowsize)[i]>0);
384  return (*t.rowsize)[i];
385  }
386  }
387  };
388 
395  template<class M, class I>
397  {
406  CommMatrixRow(M& m_, const Dune::GlobalLookupIndexSet<I>& idxset_, const I& aggidxset_)
407  : matrix(m_), idxset(idxset_), aggidxset(aggidxset_), rowsize()
408  {}
409 
413  CommMatrixRow(M& m_, const Dune::GlobalLookupIndexSet<I>& idxset_, const I& aggidxset_,
414  std::vector<size_t>& rowsize_)
415  : matrix(m_), idxset(idxset_), aggidxset(aggidxset_), rowsize(&rowsize_)
416  {}
423  {
424  typedef typename Dune::GlobalLookupIndexSet<I>::const_iterator Iter;
425  typedef typename Dune::OwnerOverlapCopyCommunication<int>::OwnerSet OwnerSet;
426 
427  for(Iter i= aggidxset.begin(), end=aggidxset.end(); i!=end; ++i)
428  if(!OwnerSet::contains(i->local().attribute())){
429  // Set to Dirchlet
430  typedef typename M::ColIterator CIter;
431  for(CIter c=matrix[i->local()].begin(), cend= matrix[i->local()].end();
432  c!= cend; ++c)
433  {
434  *c=0;
435  if(c.index()==i->local()){
436  typedef typename M::block_type::RowIterator RIter;
437  for(RIter r=c->begin(), rend=c->end();
438  r != rend; ++r)
439  (*r)[r.index()]=1;
440  }
441  }
442  }
443  }
445  M& matrix;
447  const Dune::GlobalLookupIndexSet<I>& idxset;
449  const I& aggidxset;
451  std::vector<size_t>* rowsize; // row sizes differ from sender side in overlap!
452  };
453 
454  template<class M, class I>
455  struct CommPolicy<CommMatrixRow<M,I> >
456  {
458 
463  typedef std::pair<typename I::GlobalIndex,typename M::block_type> IndexedType;
464 
466  typedef VariableSize IndexedTypeFlag;
467 
468  static std::size_t getSize(const Type& t, std::size_t i)
469  {
470  if(!t.rowsize)
471  return t.matrix[i].size();
472  else
473  {
474  assert((*t.rowsize)[i]>0);
475  return (*t.rowsize)[i];
476  }
477  }
478  };
479 
480  template<class M, class I, class RI>
482  {
484 
485  static const typename M::size_type gather(const Container& cont, std::size_t i)
486  {
487  return cont.matrix[i].size();
488  }
489  static void scatter(Container& cont, const typename M::size_type& rowsize, std::size_t i)
490  {
491  assert(rowsize);
492  cont.rowsize.getRowSize(i)=rowsize;
493  }
494 
495  };
496 
497  template<class M, class I, class RI>
499  {
501 
502  static const typename M::size_type gather(const Container& cont, std::size_t i)
503  {
504  return cont.matrix[i].size();
505  }
506  static void scatter(Container& cont, const typename M::size_type& rowsize, std::size_t i)
507  {
508  assert(rowsize);
509  if (rowsize > cont.rowsize.getCopyRowSize(i))
510  cont.rowsize.getCopyRowSize(i)=rowsize;
511  }
512 
513  };
514 
515  template<class M, class I>
517  {
518  typedef typename I::GlobalIndex GlobalIndex;
520  typedef typename M::ConstColIterator ColIter;
521 
522  static ColIter col;
524 
525  static const GlobalIndex& gather(const Container& cont, std::size_t i, std::size_t j)
526  {
527  if(j==0)
528  col=cont.matrix[i].begin();
529  else if (col!=cont.matrix[i].end())
530  ++col;
531 
532  //copy communication: different row sizes for copy rows with the same global index
533  //are possible. If all values of current matrix row are sent, send
534  //std::numeric_limits<GlobalIndex>::max()
535  //and receiver will ignore it.
536  if (col==cont.matrix[i].end()){
537  numlimits = std::numeric_limits<GlobalIndex>::max();
538  return numlimits;
539  }
540  else {
541  const typename I::IndexPair* index=cont.idxset.pair(col.index());
542  assert(index);
543  // Only send index if col is no ghost
544  if ( index->local().attribute() != 2)
545  return index->global();
546  else {
547  numlimits = std::numeric_limits<GlobalIndex>::max();
548  return numlimits;
549  }
550  }
551  }
552  static void scatter(Container& cont, const GlobalIndex& gi, std::size_t i, std::size_t j)
553  {
554  try{
555  if (gi != std::numeric_limits<GlobalIndex>::max()){
556  const typename I::IndexPair& ip=cont.aggidxset.at(gi);
557  assert(ip.global()==gi);
558  std::size_t col = ip.local();
559  cont.sparsity[i].insert(col);
560 
561  typedef typename Dune::OwnerOverlapCopyCommunication<int>::OwnerSet OwnerSet;
562  if(!OwnerSet::contains(ip.local().attribute()))
563  // preserve symmetry for overlap
564  cont.sparsity[col].insert(i);
565  }
566  }
567  catch(Dune::RangeError er){
568  // Entry not present in the new index set. Ignore!
569 #ifdef DEBUG_REPART
570  typedef typename Container::LookupIndexSet GlobalLookup;
571  typedef typename GlobalLookup::IndexPair IndexPair;
572  typedef typename Dune::OwnerOverlapCopyCommunication<int>::OwnerSet OwnerSet;
573 
574  GlobalLookup lookup(cont.aggidxset);
575  const IndexPair* pi=lookup.pair(i);
576  assert(pi);
577  if(OwnerSet::contains(pi->local().attribute())){
578  int rank;
579  MPI_Comm_rank(MPI_COMM_WORLD,&rank);
580  std::cout<<rank<<cont.aggidxset<<std::endl;
581  std::cout<<rank<<": row "<<i<<" (global="<<gi <<") not in index set for owner index "<<pi->global()<<std::endl;
582  throw er;
583  }
584 #endif
585  }
586  }
587 
588  };
589  template<class M, class I>
590  typename MatrixSparsityPatternGatherScatter<M,I>::ColIter MatrixSparsityPatternGatherScatter<M,I>::col;
591 
592  template<class M, class I>
593  typename MatrixSparsityPatternGatherScatter<M,I>::GlobalIndex MatrixSparsityPatternGatherScatter<M,I>::numlimits;
594 
595 
596  template<class M, class I>
598  {
599  typedef typename I::GlobalIndex GlobalIndex;
601  typedef typename M::ConstColIterator ColIter;
602  typedef typename std::pair<GlobalIndex,typename M::block_type> Data;
603  static ColIter col;
604  static Data datastore;
606 
607  static const Data& gather(const Container& cont, std::size_t i, std::size_t j)
608  {
609  if(j==0)
610  col=cont.matrix[i].begin();
611  else if (col!=cont.matrix[i].end())
612  ++col;
613  // copy communication: different row sizes for copy rows with the same global index
614  // are possible. If all values of current matrix row are sent, send
615  // std::numeric_limits<GlobalIndex>::max()
616  // and receiver will ignore it.
617  if (col==cont.matrix[i].end()){
618  numlimits = std::numeric_limits<GlobalIndex>::max();
619  datastore = std::make_pair(numlimits,*col);
620  return datastore;
621  }
622  else {
623  // convert local column index to global index
624  const typename I::IndexPair* index=cont.idxset.pair(col.index());
625  assert(index);
626  // Store the data to prevent reference to temporary
627  // Only send index if col is no ghost
628  if ( index->local().attribute() != 2)
629  datastore = std::make_pair(index->global(),*col);
630  else {
631  numlimits = std::numeric_limits<GlobalIndex>::max();
632  datastore = std::make_pair(numlimits,*col);
633  }
634  return datastore;
635  }
636  }
637  static void scatter(Container& cont, const Data& data, std::size_t i, std::size_t j)
638  {
639  try{
640  if (data.first != std::numeric_limits<GlobalIndex>::max()){
641  typename M::size_type column=cont.aggidxset.at(data.first).local();
642  cont.matrix[i][column]=data.second;
643  }
644  }
645  catch(Dune::RangeError er){
646  // This an overlap row and might therefore lack some entries!
647  }
648 
649  }
650  };
651 
652  template<class M, class I>
653  typename MatrixRowGatherScatter<M,I>::ColIter MatrixRowGatherScatter<M,I>::col;
654 
655  template<class M, class I>
656  typename MatrixRowGatherScatter<M,I>::Data MatrixRowGatherScatter<M,I>::datastore;
657 
658  template<class M, class I>
659  typename MatrixRowGatherScatter<M,I>::GlobalIndex MatrixRowGatherScatter<M,I>::numlimits;
660 
661  template<typename M, typename C>
662  void redistributeSparsityPattern(M& origMatrix, M& newMatrix, C& origComm, C& newComm,
664  {
665  typename C::CopySet copyflags;
666  typename C::OwnerSet ownerflags;
667  typedef typename C::ParallelIndexSet IndexSet;
668  typedef RedistributeInformation<C> RI;
669  std::vector<typename M::size_type> rowsize(newComm.indexSet().size(), 0);
670  std::vector<typename M::size_type> copyrowsize(newComm.indexSet().size(), 0);
671  std::vector<typename M::size_type> backwardscopyrowsize(origComm.indexSet().size(), 0);
672 
673  // get owner rowsizes
674  CommMatrixRowSize<M,RI> commRowSize(origMatrix, ri);
675  ri.template redistribute<MatrixRowSizeGatherScatter<M,IndexSet,RI> >(commRowSize,commRowSize);
676 
677  origComm.buildGlobalLookup();
678 
679  for (std::size_t i=0; i < newComm.indexSet().size(); i++){
680  rowsize[i] = ri.getRowSize(i);
681  }
682  // get sparsity pattern from owner rows
684  origsp(origMatrix, origComm.globalLookup(), newComm.indexSet());
686  newsp(origMatrix, origComm.globalLookup(), newComm.indexSet(), rowsize);
687 
688  ri.template redistribute<MatrixSparsityPatternGatherScatter<M,IndexSet> >(origsp,newsp);
689 
690  // build copy to owner interface to get missing matrix values for novlp case
691  if (origComm.getSolverCategory() == SolverCategory::nonoverlapping){
692  RemoteIndices<IndexSet> *ris = new RemoteIndices<IndexSet>(origComm.indexSet(),
693  newComm.indexSet(),
694  origComm.communicator());
695  ris->template rebuild<true>();
696 
697  ri.getInterface().free();
698  ri.getInterface().build(*ris,copyflags,ownerflags);
699 
700  // get copy rowsizes
701  CommMatrixRowSize<M,RI> commRowSize_copy(origMatrix, ri);
702  ri.template redistribute<MatrixCopyRowSizeGatherScatter<M,IndexSet,RI> >(commRowSize_copy,
703  commRowSize_copy);
704 
705  for (std::size_t i=0; i < newComm.indexSet().size(); i++){
706  copyrowsize[i] = ri.getCopyRowSize(i);
707  }
708  //get copy rowsizes for sender
709  ri.redistributeBackward(backwardscopyrowsize,copyrowsize);
710  for (std::size_t i=0; i < origComm.indexSet().size(); i++){
711  ri.getBackwardsCopyRowSize(i) = backwardscopyrowsize[i];
712  }
713 
714  // get sparsity pattern from copy rows
715  CommMatrixSparsityPattern<M,IndexSet> origsp_copy(origMatrix,
716  origComm.globalLookup(),
717  newComm.indexSet(),
718  backwardscopyrowsize);
719  CommMatrixSparsityPattern<M,IndexSet> newsp_copy(origMatrix, origComm.globalLookup(),
720  newComm.indexSet(), copyrowsize);
721  ri.template redistribute<MatrixSparsityPatternGatherScatter<M,IndexSet> >(origsp_copy,
722  newsp_copy);
723 
724  newsp.completeSparsityPattern(newsp_copy.sparsity);
725  newsp.storeSparsityPattern(newMatrix);
726  }
727  else
728  newsp.storeSparsityPattern(newMatrix);
729 
730 #ifdef DUNE_ISTL_WITH_CHECKING
731  // Check for symmetry
732  int ret=0;
733  typedef typename M::ConstRowIterator RIter;
734  for(RIter row=newMatrix.begin(), rend=newMatrix.end(); row != rend; ++row){
735  typedef typename M::ConstColIterator CIter;
736  for(CIter col=row->begin(), cend=row->end(); col!=cend; ++col)
737  {
738  try{
739  newMatrix[col.index()][row.index()];
740  }catch(Dune::ISTLError e){
741  std::cerr<<newComm.communicator().rank()<<": entry ("
742  <<col.index()<<","<<row.index()<<") missing! for symmetry!"<<std::endl;
743  ret=1;
744 
745  }
746 
747  }
748  }
749 
750  if(ret)
751  DUNE_THROW(ISTLError, "Matrix not symmetric!");
752 #endif
753  }
754 
755  template<typename M, typename C>
756  void redistributeMatrixEntries(M& origMatrix, M& newMatrix, C& origComm, C& newComm,
758  {
759  typedef typename C::ParallelIndexSet IndexSet;
760  typedef RedistributeInformation<C> RI;
761  typename C::OwnerSet ownerflags;
762  std::vector<typename M::size_type> rowsize(newComm.indexSet().size(), 0);
763  std::vector<typename M::size_type> copyrowsize(newComm.indexSet().size(), 0);
764  std::vector<typename M::size_type> backwardscopyrowsize(origComm.indexSet().size(), 0);
765 
766  for (std::size_t i=0; i < newComm.indexSet().size(); i++){
767  rowsize[i] = ri.getRowSize(i);
768  if (origComm.getSolverCategory() == SolverCategory::nonoverlapping){
769  copyrowsize[i] = ri.getCopyRowSize(i);
770  }
771  }
772 
773  for (std::size_t i=0; i < origComm.indexSet().size(); i++)
774  if (origComm.getSolverCategory() == SolverCategory::nonoverlapping)
775  backwardscopyrowsize[i] = ri.getBackwardsCopyRowSize(i);
776 
777 
778  if (origComm.getSolverCategory() == SolverCategory::nonoverlapping){
779  // fill sparsity pattern from copy rows
780  CommMatrixRow<M,IndexSet> origrow_copy(origMatrix, origComm.globalLookup(),
781  newComm.indexSet(), backwardscopyrowsize);
782  CommMatrixRow<M,IndexSet> newrow_copy(newMatrix, origComm.globalLookup(),
783  newComm.indexSet(),copyrowsize);
784  ri.template redistribute<MatrixRowGatherScatter<M,IndexSet> >(origrow_copy,
785  newrow_copy);
786  ri.getInterface().free();
787  RemoteIndices<IndexSet> *ris = new RemoteIndices<IndexSet>(origComm.indexSet(),
788  newComm.indexSet(),
789  origComm.communicator());
790  ris->template rebuild<true>();
791  ri.getInterface().build(*ris,ownerflags,ownerflags);
792  }
793 
795  origrow(origMatrix, origComm.globalLookup(), newComm.indexSet());
797  newrow(newMatrix, origComm.globalLookup(), newComm.indexSet(),rowsize);
798  ri.template redistribute<MatrixRowGatherScatter<M,IndexSet> >(origrow,newrow);
799  if (origComm.getSolverCategory() != static_cast<int>(SolverCategory::nonoverlapping))
800  newrow.setOverlapRowsToDirichlet();
801  if(newMatrix.N()>0&&newMatrix.N()<20)
802  printmatrix(std::cout, newMatrix, "redist", "row");
803  }
804 
821  template<typename M, typename C>
822  void redistributeMatrix(M& origMatrix, M& newMatrix, C& origComm, C& newComm,
824  {
825  ri.setNoRows(newComm.indexSet().size());
826  ri.setNoCopyRows(newComm.indexSet().size());
827  ri.setNoBackwardsCopyRows(origComm.indexSet().size());
828  redistributeSparsityPattern(origMatrix, newMatrix, origComm, newComm, ri);
829  redistributeMatrixEntries(origMatrix, newMatrix, origComm, newComm, ri);
830  }
831 #endif
832 }
833 #endif