Dune Core Modules (2.4.2)

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 "repartition.hh"
6#include <dune/common/exceptions.hh>
7#include <dune/common/parallel/indexset.hh>
8#include <dune/common/unused.hh>
9#include <dune/istl/owneroverlapcopy.hh>
15namespace Dune
16{
17 template<typename T>
18 struct RedistributeInformation
19 {
20 bool isSetup() const
21 {
22 return false;
23 }
24 template<class D>
25 void redistribute(const D& from, D& to) const
26 {
27 DUNE_UNUSED_PARAMETER(from);
28 DUNE_UNUSED_PARAMETER(to);
29 }
30
31 template<class D>
32 void redistributeBackward(D& from, const D& to) const
33 {
34 DUNE_UNUSED_PARAMETER(from);
35 DUNE_UNUSED_PARAMETER(to);
36 }
37
38 void resetSetup()
39 {}
40
41 void setNoRows(std::size_t size)
42 {
43 DUNE_UNUSED_PARAMETER(size);
44 }
45
46 void setNoCopyRows(std::size_t size)
47 {
48 DUNE_UNUSED_PARAMETER(size);
49 }
50
51 void setNoBackwardsCopyRows(std::size_t size)
52 {
53 DUNE_UNUSED_PARAMETER(size);
54 }
55
56 std::size_t getRowSize(std::size_t index) const
57 {
58 DUNE_UNUSED_PARAMETER(index);
59 return -1;
60 }
61
62 std::size_t getCopyRowSize(std::size_t index) const
63 {
64 DUNE_UNUSED_PARAMETER(index);
65 return -1;
66 }
67
68 std::size_t getBackwardsCopyRowSize(std::size_t index) const
69 {
70 DUNE_UNUSED_PARAMETER(index);
71 return -1;
72 }
73
74 };
75
76#if HAVE_MPI
77 template<typename T, typename T1>
78 class RedistributeInformation<OwnerOverlapCopyCommunication<T,T1> >
79 {
80 public:
81 typedef OwnerOverlapCopyCommunication<T,T1> Comm;
82
83 RedistributeInformation()
84 : interface(), setup_(false)
85 {}
86
87 RedistributeInterface& getInterface()
88 {
89 return interface;
90 }
91 template<typename IS>
92 void checkInterface(const IS& source,
93 const IS& target, MPI_Comm comm)
94 {
95 RemoteIndices<IS> *ri=new RemoteIndices<IS>(source, target, comm);
96 ri->template rebuild<true>();
97 Interface inf;
98 typename OwnerOverlapCopyCommunication<int>::OwnerSet flags;
99 int rank;
100 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
101 inf.free();
102 inf.build(*ri, flags, flags);
103
104
105#ifdef DEBUG_REPART
106 if(inf!=interface) {
107
108 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
109 if(rank==0)
110 std::cout<<"Interfaces do not match!"<<std::endl;
111 std::cout<<rank<<": redist interface new :"<<inf<<std::endl;
112 std::cout<<rank<<": redist interface :"<<interface<<std::endl;
113
114 throw "autsch!";
115 delete ri;
116 }else
117
118#endif
119 delete ri;
120 }
121 void setSetup()
122 {
123 setup_=true;
124 interface.strip();
125 }
126
127 void resetSetup()
128 {
129 setup_=false;
130 }
131
132 template<class GatherScatter, class D>
133 void redistribute(const D& from, D& to) const
134 {
135 BufferedCommunicator communicator;
136 communicator.template build<D>(from,to, interface);
137 communicator.template forward<GatherScatter>(from, to);
138 communicator.free();
139 }
140 template<class GatherScatter, class D>
141 void redistributeBackward(D& from, const D& to) const
142 {
143
144 BufferedCommunicator communicator;
145 communicator.template build<D>(from,to, interface);
146 communicator.template backward<GatherScatter>(from, to);
147 communicator.free();
148 }
149
150 template<class D>
151 void redistribute(const D& from, D& to) const
152 {
153 redistribute<CopyGatherScatter<D> >(from,to);
154 }
155 template<class D>
156 void redistributeBackward(D& from, const D& to) const
157 {
158 redistributeBackward<CopyGatherScatter<D> >(from,to);
159 }
160 bool isSetup() const
161 {
162 return setup_;
163 }
164
165 void reserve(std::size_t size)
166 {}
167
168 std::size_t& getRowSize(std::size_t index)
169 {
170 return rowSize[index];
171 }
172
173 std::size_t getRowSize(std::size_t index) const
174 {
175 return rowSize[index];
176 }
177
178 std::size_t& getCopyRowSize(std::size_t index)
179 {
180 return copyrowSize[index];
181 }
182
183 std::size_t getCopyRowSize(std::size_t index) const
184 {
185 return copyrowSize[index];
186 }
187
188 std::size_t& getBackwardsCopyRowSize(std::size_t index)
189 {
190 return backwardscopyrowSize[index];
191 }
192
193 std::size_t getBackwardsCopyRowSize(std::size_t index) const
194 {
195 return backwardscopyrowSize[index];
196 }
197
198 void setNoRows(std::size_t rows)
199 {
200 rowSize.resize(rows, 0);
201 }
202
203 void setNoCopyRows(std::size_t rows)
204 {
205 copyrowSize.resize(rows, 0);
206 }
207
208 void setNoBackwardsCopyRows(std::size_t rows)
209 {
210 backwardscopyrowSize.resize(rows, 0);
211 }
212
213 private:
214 std::vector<std::size_t> rowSize;
215 std::vector<std::size_t> copyrowSize;
216 std::vector<std::size_t> backwardscopyrowSize;
217 RedistributeInterface interface;
218 bool setup_;
219 };
220
229 template<class M, class RI>
230 struct CommMatrixRowSize
231 {
232 // Make the default communication policy work.
233 typedef typename M::size_type value_type;
234 typedef typename M::size_type size_type;
235
241 CommMatrixRowSize(const M& m_, RI& rowsize_)
242 : matrix(m_), rowsize(rowsize_)
243 {}
244 const M& matrix;
245 RI& rowsize;
246
247 };
248
249
258 template<class M, class I>
259 struct CommMatrixSparsityPattern
260 {
261 typedef typename M::size_type size_type;
262
269 CommMatrixSparsityPattern(const M& m_, const Dune::GlobalLookupIndexSet<I>& idxset_, const I& aggidxset_)
270 : matrix(m_), idxset(idxset_), aggidxset(aggidxset_), rowsize()
271 {}
272
280 CommMatrixSparsityPattern(const M& m_, const Dune::GlobalLookupIndexSet<I>& idxset_, const I& aggidxset_,
281 const std::vector<typename M::size_type>& rowsize_)
282 : matrix(m_), idxset(idxset_), aggidxset(aggidxset_), sparsity(aggidxset_.size()), rowsize(&rowsize_)
283 {}
284
291 void storeSparsityPattern(M& m)
292 {
293 // insert diagonal to overlap rows
294 typedef typename Dune::GlobalLookupIndexSet<I>::const_iterator IIter;
295 typedef typename Dune::OwnerOverlapCopyCommunication<int>::OwnerSet OwnerSet;
296 std::size_t nnz=0;
297#ifdef DEBUG_REPART
298 int rank;
299
300 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
301#endif
302 for(IIter i= aggidxset.begin(), end=aggidxset.end(); i!=end; ++i) {
303 if(!OwnerSet::contains(i->local().attribute())) {
304#ifdef DEBUG_REPART
305 std::cout<<rank<<" Inserting diagonal for"<<i->local()<<std::endl;
306#endif
307 sparsity[i->local()].insert(i->local());
308 }
309
310 nnz+=sparsity[i->local()].size();
311 }
312 assert( aggidxset.size()==sparsity.size());
313
314 if(nnz>0) {
315 m.setSize(aggidxset.size(), aggidxset.size(), nnz);
316 m.setBuildMode(M::row_wise);
317 typename M::CreateIterator citer=m.createbegin();
318#ifdef DEBUG_REPART
319 std::size_t idx=0;
320 bool correct=true;
321 Dune::GlobalLookupIndexSet<I> global(aggidxset);
322#endif
323 typedef typename std::vector<std::set<size_type> >::const_iterator Iter;
324 for(Iter i=sparsity.begin(), end=sparsity.end(); i!=end; ++i, ++citer)
325 {
326 typedef typename std::set<size_type>::const_iterator SIter;
327 for(SIter si=i->begin(), send=i->end(); si!=send; ++si)
328 citer.insert(*si);
329#ifdef DEBUG_REPART
330 if(i->find(idx)==i->end()) {
331 const typename I::IndexPair* gi=global.pair(idx);
332 assert(gi);
333 std::cout<<rank<<": row "<<idx<<" is missing a diagonal entry! global="<<gi->global()<<" attr="<<gi->local().attribute()<<" "<<
334 OwnerSet::contains(gi->local().attribute())<<
335 " row size="<<i->size()<<std::endl;
336 correct=false;
337 }
338 ++idx;
339#endif
340 }
341#ifdef DEBUG_REPART
342 if(!correct)
343 throw "bla";
344#endif
345 }
346 }
347
355 void completeSparsityPattern(std::vector<std::set<size_type> > add_sparsity)
356 {
357 for (unsigned int i = 0; i != sparsity.size(); ++i) {
358 if (add_sparsity[i].size() != 0) {
359 typedef std::set<size_type> Set;
360 Set tmp_set;
361 std::insert_iterator<Set> tmp_insert (tmp_set, tmp_set.begin());
362 std::set_union(add_sparsity[i].begin(), add_sparsity[i].end(),
363 sparsity[i].begin(), sparsity[i].end(), tmp_insert);
364 sparsity[i].swap(tmp_set);
365 }
366 }
367 }
368
369 const M& matrix;
370 typedef Dune::GlobalLookupIndexSet<I> LookupIndexSet;
371 const Dune::GlobalLookupIndexSet<I>& idxset;
372 const I& aggidxset;
373 std::vector<std::set<size_type> > sparsity;
374 const std::vector<size_type>* rowsize;
375 };
376
377 template<class M, class I>
378 struct CommPolicy<CommMatrixSparsityPattern<M,I> >
379 {
380 typedef CommMatrixSparsityPattern<M,I> Type;
381
386 typedef typename I::GlobalIndex IndexedType;
387
389 typedef VariableSize IndexedTypeFlag;
390
391 static typename M::size_type getSize(const Type& t, std::size_t i)
392 {
393 if(!t.rowsize)
394 return t.matrix[i].size();
395 else
396 {
397 assert((*t.rowsize)[i]>0);
398 return (*t.rowsize)[i];
399 }
400 }
401 };
402
409 template<class M, class I>
410 struct CommMatrixRow
411 {
420 CommMatrixRow(M& m_, const Dune::GlobalLookupIndexSet<I>& idxset_, const I& aggidxset_)
421 : matrix(m_), idxset(idxset_), aggidxset(aggidxset_), rowsize()
422 {}
423
427 CommMatrixRow(M& m_, const Dune::GlobalLookupIndexSet<I>& idxset_, const I& aggidxset_,
428 std::vector<size_t>& rowsize_)
429 : matrix(m_), idxset(idxset_), aggidxset(aggidxset_), rowsize(&rowsize_)
430 {}
436 void setOverlapRowsToDirichlet()
437 {
438 typedef typename Dune::GlobalLookupIndexSet<I>::const_iterator Iter;
439 typedef typename Dune::OwnerOverlapCopyCommunication<int>::OwnerSet OwnerSet;
440
441 for(Iter i= aggidxset.begin(), end=aggidxset.end(); i!=end; ++i)
442 if(!OwnerSet::contains(i->local().attribute())) {
443 // Set to Dirchlet
444 typedef typename M::ColIterator CIter;
445 for(CIter c=matrix[i->local()].begin(), cend= matrix[i->local()].end();
446 c!= cend; ++c)
447 {
448 *c=0;
449 if(c.index()==i->local()) {
450 typedef typename M::block_type::RowIterator RIter;
451 for(RIter r=c->begin(), rend=c->end();
452 r != rend; ++r)
453 (*r)[r.index()]=1;
454 }
455 }
456 }
457 }
459 M& matrix;
461 const Dune::GlobalLookupIndexSet<I>& idxset;
463 const I& aggidxset;
465 std::vector<size_t>* rowsize; // row sizes differ from sender side in overlap!
466 };
467
468 template<class M, class I>
469 struct CommPolicy<CommMatrixRow<M,I> >
470 {
471 typedef CommMatrixRow<M,I> Type;
472
477 typedef std::pair<typename I::GlobalIndex,typename M::block_type> IndexedType;
478
480 typedef VariableSize IndexedTypeFlag;
481
482 static std::size_t getSize(const Type& t, std::size_t i)
483 {
484 if(!t.rowsize)
485 return t.matrix[i].size();
486 else
487 {
488 assert((*t.rowsize)[i]>0);
489 return (*t.rowsize)[i];
490 }
491 }
492 };
493
494 template<class M, class I, class RI>
495 struct MatrixRowSizeGatherScatter
496 {
497 typedef CommMatrixRowSize<M,RI> Container;
498
499 static const typename M::size_type gather(const Container& cont, std::size_t i)
500 {
501 return cont.matrix[i].size();
502 }
503 static void scatter(Container& cont, const typename M::size_type& rowsize, std::size_t i)
504 {
505 assert(rowsize);
506 cont.rowsize.getRowSize(i)=rowsize;
507 }
508
509 };
510
511 template<class M, class I, class RI>
512 struct MatrixCopyRowSizeGatherScatter
513 {
514 typedef CommMatrixRowSize<M,RI> Container;
515
516 static const typename M::size_type gather(const Container& cont, std::size_t i)
517 {
518 return cont.matrix[i].size();
519 }
520 static void scatter(Container& cont, const typename M::size_type& rowsize, std::size_t i)
521 {
522 assert(rowsize);
523 if (rowsize > cont.rowsize.getCopyRowSize(i))
524 cont.rowsize.getCopyRowSize(i)=rowsize;
525 }
526
527 };
528
529 template<class M, class I>
530 struct MatrixSparsityPatternGatherScatter
531 {
532 typedef typename I::GlobalIndex GlobalIndex;
533 typedef CommMatrixSparsityPattern<M,I> Container;
534 typedef typename M::ConstColIterator ColIter;
535
536 static ColIter col;
537 static GlobalIndex numlimits;
538
539 static const GlobalIndex& gather(const Container& cont, std::size_t i, std::size_t j)
540 {
541 if(j==0)
542 col=cont.matrix[i].begin();
543 else if (col!=cont.matrix[i].end())
544 ++col;
545
546 //copy communication: different row sizes for copy rows with the same global index
547 //are possible. If all values of current matrix row are sent, send
548 //std::numeric_limits<GlobalIndex>::max()
549 //and receiver will ignore it.
550 if (col==cont.matrix[i].end()) {
551 numlimits = std::numeric_limits<GlobalIndex>::max();
552 return numlimits;
553 }
554 else {
555 const typename I::IndexPair* index=cont.idxset.pair(col.index());
556 assert(index);
557 // Only send index if col is no ghost
558 if ( index->local().attribute() != 2)
559 return index->global();
560 else {
561 numlimits = std::numeric_limits<GlobalIndex>::max();
562 return numlimits;
563 }
564 }
565 }
566 static void scatter(Container& cont, const GlobalIndex& gi, std::size_t i, std::size_t j)
567 {
568 DUNE_UNUSED_PARAMETER(j);
569 try{
570 if (gi != std::numeric_limits<GlobalIndex>::max()) {
571 const typename I::IndexPair& ip=cont.aggidxset.at(gi);
572 assert(ip.global()==gi);
573 std::size_t column = ip.local();
574 cont.sparsity[i].insert(column);
575
576 typedef typename Dune::OwnerOverlapCopyCommunication<int>::OwnerSet OwnerSet;
577 if(!OwnerSet::contains(ip.local().attribute()))
578 // preserve symmetry for overlap
579 cont.sparsity[column].insert(i);
580 }
581 }
582 catch(Dune::RangeError er) {
583 // Entry not present in the new index set. Ignore!
584#ifdef DEBUG_REPART
585 typedef typename Container::LookupIndexSet GlobalLookup;
586 typedef typename GlobalLookup::IndexPair IndexPair;
587 typedef typename Dune::OwnerOverlapCopyCommunication<int>::OwnerSet OwnerSet;
588
589 GlobalLookup lookup(cont.aggidxset);
590 const IndexPair* pi=lookup.pair(i);
591 assert(pi);
592 if(OwnerSet::contains(pi->local().attribute())) {
593 int rank;
594 MPI_Comm_rank(MPI_COMM_WORLD,&rank);
595 std::cout<<rank<<cont.aggidxset<<std::endl;
596 std::cout<<rank<<": row "<<i<<" (global="<<gi <<") not in index set for owner index "<<pi->global()<<std::endl;
597 throw er;
598 }
599#endif
600 }
601 }
602
603 };
604 template<class M, class I>
605 typename MatrixSparsityPatternGatherScatter<M,I>::ColIter MatrixSparsityPatternGatherScatter<M,I>::col;
606
607 template<class M, class I>
608 typename MatrixSparsityPatternGatherScatter<M,I>::GlobalIndex MatrixSparsityPatternGatherScatter<M,I>::numlimits;
609
610
611 template<class M, class I>
612 struct MatrixRowGatherScatter
613 {
614 typedef typename I::GlobalIndex GlobalIndex;
615 typedef CommMatrixRow<M,I> Container;
616 typedef typename M::ConstColIterator ColIter;
617 typedef typename std::pair<GlobalIndex,typename M::block_type> Data;
618 static ColIter col;
619 static Data datastore;
620 static GlobalIndex numlimits;
621
622 static const Data& gather(const Container& cont, std::size_t i, std::size_t j)
623 {
624 if(j==0)
625 col=cont.matrix[i].begin();
626 else if (col!=cont.matrix[i].end())
627 ++col;
628 // copy communication: different row sizes for copy rows with the same global index
629 // are possible. If all values of current matrix row are sent, send
630 // std::numeric_limits<GlobalIndex>::max()
631 // and receiver will ignore it.
632 if (col==cont.matrix[i].end()) {
633 numlimits = std::numeric_limits<GlobalIndex>::max();
634 datastore = std::make_pair(numlimits,*col);
635 return datastore;
636 }
637 else {
638 // convert local column index to global index
639 const typename I::IndexPair* index=cont.idxset.pair(col.index());
640 assert(index);
641 // Store the data to prevent reference to temporary
642 // Only send index if col is no ghost
643 if ( index->local().attribute() != 2)
644 datastore = std::make_pair(index->global(),*col);
645 else {
646 numlimits = std::numeric_limits<GlobalIndex>::max();
647 datastore = std::make_pair(numlimits,*col);
648 }
649 return datastore;
650 }
651 }
652 static void scatter(Container& cont, const Data& data, std::size_t i, std::size_t j)
653 {
654 DUNE_UNUSED_PARAMETER(j);
655 try{
656 if (data.first != std::numeric_limits<GlobalIndex>::max()) {
657 typename M::size_type column=cont.aggidxset.at(data.first).local();
658 cont.matrix[i][column]=data.second;
659 }
660 }
661 catch(Dune::RangeError er) {
662 // This an overlap row and might therefore lack some entries!
663 }
664
665 }
666 };
667
668 template<class M, class I>
669 typename MatrixRowGatherScatter<M,I>::ColIter MatrixRowGatherScatter<M,I>::col;
670
671 template<class M, class I>
672 typename MatrixRowGatherScatter<M,I>::Data MatrixRowGatherScatter<M,I>::datastore;
673
674 template<class M, class I>
675 typename MatrixRowGatherScatter<M,I>::GlobalIndex MatrixRowGatherScatter<M,I>::numlimits;
676
677 template<typename M, typename C>
678 void redistributeSparsityPattern(M& origMatrix, M& newMatrix, C& origComm, C& newComm,
679 RedistributeInformation<C>& ri)
680 {
681 typename C::CopySet copyflags;
682 typename C::OwnerSet ownerflags;
683 typedef typename C::ParallelIndexSet IndexSet;
684 typedef RedistributeInformation<C> RI;
685 std::vector<typename M::size_type> rowsize(newComm.indexSet().size(), 0);
686 std::vector<typename M::size_type> copyrowsize(newComm.indexSet().size(), 0);
687 std::vector<typename M::size_type> backwardscopyrowsize(origComm.indexSet().size(), 0);
688
689 // get owner rowsizes
690 CommMatrixRowSize<M,RI> commRowSize(origMatrix, ri);
691 ri.template redistribute<MatrixRowSizeGatherScatter<M,IndexSet,RI> >(commRowSize,commRowSize);
692
693 origComm.buildGlobalLookup();
694
695 for (std::size_t i=0; i < newComm.indexSet().size(); i++) {
696 rowsize[i] = ri.getRowSize(i);
697 }
698 // get sparsity pattern from owner rows
699 CommMatrixSparsityPattern<M,IndexSet>
700 origsp(origMatrix, origComm.globalLookup(), newComm.indexSet());
701 CommMatrixSparsityPattern<M,IndexSet>
702 newsp(origMatrix, origComm.globalLookup(), newComm.indexSet(), rowsize);
703
704 ri.template redistribute<MatrixSparsityPatternGatherScatter<M,IndexSet> >(origsp,newsp);
705
706 // build copy to owner interface to get missing matrix values for novlp case
707 if (origComm.getSolverCategory() == SolverCategory::nonoverlapping) {
708 RemoteIndices<IndexSet> *ris = new RemoteIndices<IndexSet>(origComm.indexSet(),
709 newComm.indexSet(),
710 origComm.communicator());
711 ris->template rebuild<true>();
712
713 ri.getInterface().free();
714 ri.getInterface().build(*ris,copyflags,ownerflags);
715
716 // get copy rowsizes
717 CommMatrixRowSize<M,RI> commRowSize_copy(origMatrix, ri);
718 ri.template redistribute<MatrixCopyRowSizeGatherScatter<M,IndexSet,RI> >(commRowSize_copy,
719 commRowSize_copy);
720
721 for (std::size_t i=0; i < newComm.indexSet().size(); i++) {
722 copyrowsize[i] = ri.getCopyRowSize(i);
723 }
724 //get copy rowsizes for sender
725 ri.redistributeBackward(backwardscopyrowsize,copyrowsize);
726 for (std::size_t i=0; i < origComm.indexSet().size(); i++) {
727 ri.getBackwardsCopyRowSize(i) = backwardscopyrowsize[i];
728 }
729
730 // get sparsity pattern from copy rows
731 CommMatrixSparsityPattern<M,IndexSet> origsp_copy(origMatrix,
732 origComm.globalLookup(),
733 newComm.indexSet(),
734 backwardscopyrowsize);
735 CommMatrixSparsityPattern<M,IndexSet> newsp_copy(origMatrix, origComm.globalLookup(),
736 newComm.indexSet(), copyrowsize);
737 ri.template redistribute<MatrixSparsityPatternGatherScatter<M,IndexSet> >(origsp_copy,
738 newsp_copy);
739
740 newsp.completeSparsityPattern(newsp_copy.sparsity);
741 newsp.storeSparsityPattern(newMatrix);
742 }
743 else
744 newsp.storeSparsityPattern(newMatrix);
745
746#ifdef DUNE_ISTL_WITH_CHECKING
747 // Check for symmetry
748 int ret=0;
749 typedef typename M::ConstRowIterator RIter;
750 for(RIter row=newMatrix.begin(), rend=newMatrix.end(); row != rend; ++row) {
751 typedef typename M::ConstColIterator CIter;
752 for(CIter col=row->begin(), cend=row->end(); col!=cend; ++col)
753 {
754 try{
755 newMatrix[col.index()][row.index()];
756 }catch(Dune::ISTLError e) {
757 std::cerr<<newComm.communicator().rank()<<": entry ("
758 <<col.index()<<","<<row.index()<<") missing! for symmetry!"<<std::endl;
759 ret=1;
760
761 }
762
763 }
764 }
765
766 if(ret)
767 DUNE_THROW(ISTLError, "Matrix not symmetric!");
768#endif
769 }
770
771 template<typename M, typename C>
772 void redistributeMatrixEntries(M& origMatrix, M& newMatrix, C& origComm, C& newComm,
773 RedistributeInformation<C>& ri)
774 {
775 typedef typename C::ParallelIndexSet IndexSet;
776 typename C::OwnerSet ownerflags;
777 std::vector<typename M::size_type> rowsize(newComm.indexSet().size(), 0);
778 std::vector<typename M::size_type> copyrowsize(newComm.indexSet().size(), 0);
779 std::vector<typename M::size_type> backwardscopyrowsize(origComm.indexSet().size(), 0);
780
781 for (std::size_t i=0; i < newComm.indexSet().size(); i++) {
782 rowsize[i] = ri.getRowSize(i);
783 if (origComm.getSolverCategory() == SolverCategory::nonoverlapping) {
784 copyrowsize[i] = ri.getCopyRowSize(i);
785 }
786 }
787
788 for (std::size_t i=0; i < origComm.indexSet().size(); i++)
789 if (origComm.getSolverCategory() == SolverCategory::nonoverlapping)
790 backwardscopyrowsize[i] = ri.getBackwardsCopyRowSize(i);
791
792
793 if (origComm.getSolverCategory() == SolverCategory::nonoverlapping) {
794 // fill sparsity pattern from copy rows
795 CommMatrixRow<M,IndexSet> origrow_copy(origMatrix, origComm.globalLookup(),
796 newComm.indexSet(), backwardscopyrowsize);
797 CommMatrixRow<M,IndexSet> newrow_copy(newMatrix, origComm.globalLookup(),
798 newComm.indexSet(),copyrowsize);
799 ri.template redistribute<MatrixRowGatherScatter<M,IndexSet> >(origrow_copy,
800 newrow_copy);
801 ri.getInterface().free();
802 RemoteIndices<IndexSet> *ris = new RemoteIndices<IndexSet>(origComm.indexSet(),
803 newComm.indexSet(),
804 origComm.communicator());
805 ris->template rebuild<true>();
806 ri.getInterface().build(*ris,ownerflags,ownerflags);
807 }
808
809 CommMatrixRow<M,IndexSet>
810 origrow(origMatrix, origComm.globalLookup(), newComm.indexSet());
811 CommMatrixRow<M,IndexSet>
812 newrow(newMatrix, origComm.globalLookup(), newComm.indexSet(),rowsize);
813 ri.template redistribute<MatrixRowGatherScatter<M,IndexSet> >(origrow,newrow);
814 if (origComm.getSolverCategory() != static_cast<int>(SolverCategory::nonoverlapping))
815 newrow.setOverlapRowsToDirichlet();
816 if(newMatrix.N()>0&&newMatrix.N()<20)
817 printmatrix(std::cout, newMatrix, "redist", "row");
818 }
819
836 template<typename M, typename C>
837 void redistributeMatrix(M& origMatrix, M& newMatrix, C& origComm, C& newComm,
838 RedistributeInformation<C>& ri)
839 {
840 ri.setNoRows(newComm.indexSet().size());
841 ri.setNoCopyRows(newComm.indexSet().size());
842 ri.setNoBackwardsCopyRows(origComm.indexSet().size());
843 redistributeSparsityPattern(origMatrix, newMatrix, origComm, newComm, ri);
844 redistributeMatrixEntries(origMatrix, newMatrix, origComm, newComm, ri);
845 }
846#endif
847}
848#endif
Dune::ConstArrayListIterator
A constant random access iterator for the Dune::ArrayList class.
Definition: arraylist.hh:379
Dune::EnumItem
A set consisting only of one item.
Definition: enumset.hh:60
Dune::GlobalLookupIndexSet
Decorates an index set with the possibility to find a global index that is mapped to a specific local...
Definition: indexset.hh:506
Dune::ISTLError
derive error class from the base class in common
Definition: istlexception.hh:16
Dune::RangeError
Default exception class for range errors.
Definition: exceptions.hh:279
exceptions.hh
A few common exception classes.
Dune::GlobalLookupIndexSet::pair
const IndexPair * pair(const std::size_t &local) const
Get the index pair corresponding to a local index.
DUNE_THROW
#define DUNE_THROW(E, m)
Definition: exceptions.hh:243
Dune::printmatrix
void printmatrix(std::ostream &s, const M &A, std::string title, std::string rowtext, int width=10, int precision=2)
Print a generic block matrix.
Definition: io.hh:259
indexset.hh
Provides a map between global and local indices.
Dune
Dune namespace.
Definition: alignment.hh:10
Dune::redistributeMatrix
void redistributeMatrix(M &origMatrix, M &newMatrix, C &origComm, C &newComm, RedistributeInformation< C > &ri)
Redistribute a matrix according to given domain decompositions.
Definition: matrixredistribute.hh:837
owneroverlapcopy.hh
Classes providing communication interfaces for overlapping Schwarz methods.
repartition.hh
Functionality for redistributing a parallel index set using graph partitioning.
Dune::CommMatrixRowSize
Utility class to communicate and set the row sizes of a redistributed matrix.
Definition: matrixredistribute.hh:231
Dune::CommMatrixRowSize::CommMatrixRowSize
CommMatrixRowSize(const M &m_, RI &rowsize_)
Constructor.
Definition: matrixredistribute.hh:241
Dune::CommMatrixRow
Utility class for comunicating the matrix entries.
Definition: matrixredistribute.hh:411
Dune::CommMatrixRow::rowsize
std::vector< size_t > * rowsize
row size information for the receiving side.
Definition: matrixredistribute.hh:465
Dune::CommMatrixRow::matrix
M & matrix
The matrix to communicate the values of.
Definition: matrixredistribute.hh:459
Dune::CommMatrixRow::CommMatrixRow
CommMatrixRow(M &m_, const Dune::GlobalLookupIndexSet< I > &idxset_, const I &aggidxset_, std::vector< size_t > &rowsize_)
Constructor.
Definition: matrixredistribute.hh:427
Dune::CommMatrixRow::idxset
const Dune::GlobalLookupIndexSet< I > & idxset
Index set for the original matrix.
Definition: matrixredistribute.hh:461
Dune::CommMatrixRow::setOverlapRowsToDirichlet
void setOverlapRowsToDirichlet()
Sets the non-owner rows correctly as Dirichlet boundaries.
Definition: matrixredistribute.hh:436
Dune::CommMatrixRow::aggidxset
const I & aggidxset
Index set for the redistributed matrix.
Definition: matrixredistribute.hh:463
Dune::CommMatrixRow::CommMatrixRow
CommMatrixRow(M &m_, const Dune::GlobalLookupIndexSet< I > &idxset_, const I &aggidxset_)
Constructor.
Definition: matrixredistribute.hh:420
Dune::CommMatrixSparsityPattern
Utility class to communicate and build the sparsity pattern of a redistributed matrix.
Definition: matrixredistribute.hh:260
Dune::CommMatrixSparsityPattern::storeSparsityPattern
void storeSparsityPattern(M &m)
Creates and stores the sparsity pattern of the redistributed matrix.
Definition: matrixredistribute.hh:291
Dune::CommMatrixSparsityPattern::completeSparsityPattern
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:355
Dune::CommMatrixSparsityPattern::CommMatrixSparsityPattern
CommMatrixSparsityPattern(const M &m_, const Dune::GlobalLookupIndexSet< I > &idxset_, const I &aggidxset_)
Constructor for the original side.
Definition: matrixredistribute.hh:269
Dune::CommMatrixSparsityPattern::CommMatrixSparsityPattern
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:280
Dune::CommPolicy
Default policy used for communicating an indexed type.
Definition: communicator.hh:120
Dune::CommPolicy::IndexedType
V::value_type IndexedType
The type we get at each index with operator[].
Definition: communicator.hh:139
Dune::CommPolicy::getSize
static int getSize(const V &, int index)
Get the number of primitve elements at that index.
Dune::CommPolicy::IndexedTypeFlag
SizeOne IndexedTypeFlag
Whether the indexed type has variable size or there is always one value at each index.
Definition: communicator.hh:145
Dune::CommPolicy::Type
V Type
The type the policy is for.
Definition: communicator.hh:132
Dune::SolverCategory::nonoverlapping
@ nonoverlapping
Category for on overlapping solvers.
Definition: solvercategory.hh:23
Dune::VariableSize
Flag for marking indexed data structures where the data at each index may be a variable multiple of a...
Definition: communicator.hh:110
unused.hh
Definition of the DUNE_UNUSED macro for the case that config.h is not available.
DUNE_UNUSED_PARAMETER
#define DUNE_UNUSED_PARAMETER(parm)
A macro to mark intentional unused function parameters with.
Definition: unused.hh:18
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Dec 21, 23:30, 2024)