Dune Core Modules (2.3.1)

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_MATRIXREDIST_HH
4#define DUNE_MATRIXREDIST_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 int rank;
109 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
110 if(rank==0)
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;
114
115 throw "autsch!";
116 delete ri;
117 }else
118
119#endif
120 delete ri;
121 }
122 void setSetup()
123 {
124 setup_=true;
125 interface.strip();
126 }
127
128 void resetSetup()
129 {
130 setup_=false;
131 }
132
133 template<class GatherScatter, class D>
134 void redistribute(const D& from, D& to) const
135 {
136 BufferedCommunicator communicator;
137 communicator.template build<D>(from,to, interface);
138 communicator.template forward<GatherScatter>(from, to);
139 communicator.free();
140 }
141 template<class GatherScatter, class D>
142 void redistributeBackward(D& from, const D& to) const
143 {
144
145 BufferedCommunicator communicator;
146 communicator.template build<D>(from,to, interface);
147 communicator.template backward<GatherScatter>(from, to);
148 communicator.free();
149 }
150
151 template<class D>
152 void redistribute(const D& from, D& to) const
153 {
154 redistribute<CopyGatherScatter<D> >(from,to);
155 }
156 template<class D>
157 void redistributeBackward(D& from, const D& to) const
158 {
159 redistributeBackward<CopyGatherScatter<D> >(from,to);
160 }
161 bool isSetup() const
162 {
163 return setup_;
164 }
165
166 void reserve(std::size_t size)
167 {}
168
169 std::size_t& getRowSize(std::size_t index)
170 {
171 return rowSize[index];
172 }
173
174 std::size_t getRowSize(std::size_t index) const
175 {
176 return rowSize[index];
177 }
178
179 std::size_t& getCopyRowSize(std::size_t index)
180 {
181 return copyrowSize[index];
182 }
183
184 std::size_t getCopyRowSize(std::size_t index) const
185 {
186 return copyrowSize[index];
187 }
188
189 std::size_t& getBackwardsCopyRowSize(std::size_t index)
190 {
191 return backwardscopyrowSize[index];
192 }
193
194 std::size_t getBackwardsCopyRowSize(std::size_t index) const
195 {
196 return backwardscopyrowSize[index];
197 }
198
199 void setNoRows(std::size_t rows)
200 {
201 rowSize.resize(rows, 0);
202 }
203
204 void setNoCopyRows(std::size_t rows)
205 {
206 copyrowSize.resize(rows, 0);
207 }
208
209 void setNoBackwardsCopyRows(std::size_t rows)
210 {
211 backwardscopyrowSize.resize(rows, 0);
212 }
213
214 private:
215 std::vector<std::size_t> rowSize;
216 std::vector<std::size_t> copyrowSize;
217 std::vector<std::size_t> backwardscopyrowSize;
218 RedistributeInterface interface;
219 bool setup_;
220 };
221
230 template<class M, class RI>
231 struct CommMatrixRowSize
232 {
233 // Make the default communication policy work.
234 typedef typename M::size_type value_type;
235 typedef typename M::size_type size_type;
236
242 CommMatrixRowSize(const M& m_, RI& rowsize_)
243 : matrix(m_), rowsize(rowsize_)
244 {}
245 const M& matrix;
246 RI& rowsize;
247
248 };
249
250
259 template<class M, class I>
260 struct CommMatrixSparsityPattern
261 {
262 typedef typename M::size_type size_type;
263
270 CommMatrixSparsityPattern(const M& m_, const Dune::GlobalLookupIndexSet<I>& idxset_, const I& aggidxset_)
271 : matrix(m_), idxset(idxset_), aggidxset(aggidxset_), rowsize()
272 {}
273
281 CommMatrixSparsityPattern(const M& m_, const Dune::GlobalLookupIndexSet<I>& idxset_, const I& aggidxset_,
282 const std::vector<typename M::size_type>& rowsize_)
283 : matrix(m_), idxset(idxset_), aggidxset(aggidxset_), sparsity(aggidxset_.size()), rowsize(&rowsize_)
284 {}
285
292 void storeSparsityPattern(M& m)
293 {
294 // insert diagonal to overlap rows
295 typedef typename Dune::GlobalLookupIndexSet<I>::const_iterator IIter;
296 typedef typename Dune::OwnerOverlapCopyCommunication<int>::OwnerSet OwnerSet;
297 std::size_t nnz=0;
298#ifdef DEBUG_REPART
299 int rank;
300
301 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
302#endif
303 for(IIter i= aggidxset.begin(), end=aggidxset.end(); i!=end; ++i) {
304 if(!OwnerSet::contains(i->local().attribute())) {
305#ifdef DEBUG_REPART
306 std::cout<<rank<<" Inserting diagonal for"<<i->local()<<std::endl;
307#endif
308 sparsity[i->local()].insert(i->local());
309 }
310
311 nnz+=sparsity[i->local()].size();
312 }
313 assert( aggidxset.size()==sparsity.size());
314
315 if(nnz>0) {
316 m.setSize(aggidxset.size(), aggidxset.size(), nnz);
317 m.setBuildMode(M::row_wise);
318 typename M::CreateIterator citer=m.createbegin();
319#ifdef DEBUG_REPART
320 std::size_t idx=0;
321 bool correct=true;
322 Dune::GlobalLookupIndexSet<I> global(aggidxset);
323#endif
324 typedef typename std::vector<std::set<size_type> >::const_iterator Iter;
325 for(Iter i=sparsity.begin(), end=sparsity.end(); i!=end; ++i, ++citer)
326 {
327 typedef typename std::set<size_type>::const_iterator SIter;
328 for(SIter si=i->begin(), send=i->end(); si!=send; ++si)
329 citer.insert(*si);
330#ifdef DEBUG_REPART
331 if(i->find(idx)==i->end()) {
332 const typename I::IndexPair* gi=global.pair(idx);
333 assert(gi);
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;
337 correct=false;
338 }
339 ++idx;
340#endif
341 }
342#ifdef DEBUG_REPART
343 if(!correct)
344 throw "bla";
345#endif
346 }
347 }
348
356 void completeSparsityPattern(std::vector<std::set<size_type> > add_sparsity)
357 {
358 for (unsigned int i = 0; i != sparsity.size(); ++i) {
359 if (add_sparsity[i].size() != 0) {
360 typedef std::set<size_type> Set;
361 Set tmp_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(),
364 sparsity[i].begin(), sparsity[i].end(), tmp_insert);
365 sparsity[i].swap(tmp_set);
366 }
367 }
368 }
369
370 const M& matrix;
371 typedef Dune::GlobalLookupIndexSet<I> LookupIndexSet;
372 const Dune::GlobalLookupIndexSet<I>& idxset;
373 const I& aggidxset;
374 std::vector<std::set<size_type> > sparsity;
375 const std::vector<size_type>* rowsize;
376 };
377
378 template<class M, class I>
379 struct CommPolicy<CommMatrixSparsityPattern<M,I> >
380 {
381 typedef CommMatrixSparsityPattern<M,I> Type;
382
387 typedef typename I::GlobalIndex IndexedType;
388
390 typedef VariableSize IndexedTypeFlag;
391
392 static typename M::size_type getSize(const Type& t, std::size_t i)
393 {
394 if(!t.rowsize)
395 return t.matrix[i].size();
396 else
397 {
398 assert((*t.rowsize)[i]>0);
399 return (*t.rowsize)[i];
400 }
401 }
402 };
403
410 template<class M, class I>
411 struct CommMatrixRow
412 {
421 CommMatrixRow(M& m_, const Dune::GlobalLookupIndexSet<I>& idxset_, const I& aggidxset_)
422 : matrix(m_), idxset(idxset_), aggidxset(aggidxset_), rowsize()
423 {}
424
428 CommMatrixRow(M& m_, const Dune::GlobalLookupIndexSet<I>& idxset_, const I& aggidxset_,
429 std::vector<size_t>& rowsize_)
430 : matrix(m_), idxset(idxset_), aggidxset(aggidxset_), rowsize(&rowsize_)
431 {}
437 void setOverlapRowsToDirichlet()
438 {
439 typedef typename Dune::GlobalLookupIndexSet<I>::const_iterator Iter;
440 typedef typename Dune::OwnerOverlapCopyCommunication<int>::OwnerSet OwnerSet;
441
442 for(Iter i= aggidxset.begin(), end=aggidxset.end(); i!=end; ++i)
443 if(!OwnerSet::contains(i->local().attribute())) {
444 // Set to Dirchlet
445 typedef typename M::ColIterator CIter;
446 for(CIter c=matrix[i->local()].begin(), cend= matrix[i->local()].end();
447 c!= cend; ++c)
448 {
449 *c=0;
450 if(c.index()==i->local()) {
451 typedef typename M::block_type::RowIterator RIter;
452 for(RIter r=c->begin(), rend=c->end();
453 r != rend; ++r)
454 (*r)[r.index()]=1;
455 }
456 }
457 }
458 }
460 M& matrix;
462 const Dune::GlobalLookupIndexSet<I>& idxset;
464 const I& aggidxset;
466 std::vector<size_t>* rowsize; // row sizes differ from sender side in overlap!
467 };
468
469 template<class M, class I>
470 struct CommPolicy<CommMatrixRow<M,I> >
471 {
472 typedef CommMatrixRow<M,I> Type;
473
478 typedef std::pair<typename I::GlobalIndex,typename M::block_type> IndexedType;
479
481 typedef VariableSize IndexedTypeFlag;
482
483 static std::size_t getSize(const Type& t, std::size_t i)
484 {
485 if(!t.rowsize)
486 return t.matrix[i].size();
487 else
488 {
489 assert((*t.rowsize)[i]>0);
490 return (*t.rowsize)[i];
491 }
492 }
493 };
494
495 template<class M, class I, class RI>
496 struct MatrixRowSizeGatherScatter
497 {
498 typedef CommMatrixRowSize<M,RI> Container;
499
500 static const typename M::size_type gather(const Container& cont, std::size_t i)
501 {
502 return cont.matrix[i].size();
503 }
504 static void scatter(Container& cont, const typename M::size_type& rowsize, std::size_t i)
505 {
506 assert(rowsize);
507 cont.rowsize.getRowSize(i)=rowsize;
508 }
509
510 };
511
512 template<class M, class I, class RI>
513 struct MatrixCopyRowSizeGatherScatter
514 {
515 typedef CommMatrixRowSize<M,RI> Container;
516
517 static const typename M::size_type gather(const Container& cont, std::size_t i)
518 {
519 return cont.matrix[i].size();
520 }
521 static void scatter(Container& cont, const typename M::size_type& rowsize, std::size_t i)
522 {
523 assert(rowsize);
524 if (rowsize > cont.rowsize.getCopyRowSize(i))
525 cont.rowsize.getCopyRowSize(i)=rowsize;
526 }
527
528 };
529
530 template<class M, class I>
531 struct MatrixSparsityPatternGatherScatter
532 {
533 typedef typename I::GlobalIndex GlobalIndex;
534 typedef CommMatrixSparsityPattern<M,I> Container;
535 typedef typename M::ConstColIterator ColIter;
536
537 static ColIter col;
538 static GlobalIndex numlimits;
539
540 static const GlobalIndex& gather(const Container& cont, std::size_t i, std::size_t j)
541 {
542 if(j==0)
543 col=cont.matrix[i].begin();
544 else if (col!=cont.matrix[i].end())
545 ++col;
546
547 //copy communication: different row sizes for copy rows with the same global index
548 //are possible. If all values of current matrix row are sent, send
549 //std::numeric_limits<GlobalIndex>::max()
550 //and receiver will ignore it.
551 if (col==cont.matrix[i].end()) {
552 numlimits = std::numeric_limits<GlobalIndex>::max();
553 return numlimits;
554 }
555 else {
556 const typename I::IndexPair* index=cont.idxset.pair(col.index());
557 assert(index);
558 // Only send index if col is no ghost
559 if ( index->local().attribute() != 2)
560 return index->global();
561 else {
562 numlimits = std::numeric_limits<GlobalIndex>::max();
563 return numlimits;
564 }
565 }
566 }
567 static void scatter(Container& cont, const GlobalIndex& gi, std::size_t i, std::size_t j)
568 {
569 DUNE_UNUSED_PARAMETER(j);
570 try{
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();
575 cont.sparsity[i].insert(col);
576
577 typedef typename Dune::OwnerOverlapCopyCommunication<int>::OwnerSet OwnerSet;
578 if(!OwnerSet::contains(ip.local().attribute()))
579 // preserve symmetry for overlap
580 cont.sparsity[col].insert(i);
581 }
582 }
583 catch(Dune::RangeError er) {
584 // Entry not present in the new index set. Ignore!
585#ifdef DEBUG_REPART
586 typedef typename Container::LookupIndexSet GlobalLookup;
587 typedef typename GlobalLookup::IndexPair IndexPair;
588 typedef typename Dune::OwnerOverlapCopyCommunication<int>::OwnerSet OwnerSet;
589
590 GlobalLookup lookup(cont.aggidxset);
591 const IndexPair* pi=lookup.pair(i);
592 assert(pi);
593 if(OwnerSet::contains(pi->local().attribute())) {
594 int rank;
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;
598 throw er;
599 }
600#endif
601 }
602 }
603
604 };
605 template<class M, class I>
606 typename MatrixSparsityPatternGatherScatter<M,I>::ColIter MatrixSparsityPatternGatherScatter<M,I>::col;
607
608 template<class M, class I>
609 typename MatrixSparsityPatternGatherScatter<M,I>::GlobalIndex MatrixSparsityPatternGatherScatter<M,I>::numlimits;
610
611
612 template<class M, class I>
613 struct MatrixRowGatherScatter
614 {
615 typedef typename I::GlobalIndex GlobalIndex;
616 typedef CommMatrixRow<M,I> Container;
617 typedef typename M::ConstColIterator ColIter;
618 typedef typename std::pair<GlobalIndex,typename M::block_type> Data;
619 static ColIter col;
620 static Data datastore;
621 static GlobalIndex numlimits;
622
623 static const Data& gather(const Container& cont, std::size_t i, std::size_t j)
624 {
625 if(j==0)
626 col=cont.matrix[i].begin();
627 else if (col!=cont.matrix[i].end())
628 ++col;
629 // copy communication: different row sizes for copy rows with the same global index
630 // are possible. If all values of current matrix row are sent, send
631 // std::numeric_limits<GlobalIndex>::max()
632 // and receiver will ignore it.
633 if (col==cont.matrix[i].end()) {
634 numlimits = std::numeric_limits<GlobalIndex>::max();
635 datastore = std::make_pair(numlimits,*col);
636 return datastore;
637 }
638 else {
639 // convert local column index to global index
640 const typename I::IndexPair* index=cont.idxset.pair(col.index());
641 assert(index);
642 // Store the data to prevent reference to temporary
643 // Only send index if col is no ghost
644 if ( index->local().attribute() != 2)
645 datastore = std::make_pair(index->global(),*col);
646 else {
647 numlimits = std::numeric_limits<GlobalIndex>::max();
648 datastore = std::make_pair(numlimits,*col);
649 }
650 return datastore;
651 }
652 }
653 static void scatter(Container& cont, const Data& data, std::size_t i, std::size_t j)
654 {
655 DUNE_UNUSED_PARAMETER(j);
656 try{
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;
660 }
661 }
662 catch(Dune::RangeError er) {
663 // This an overlap row and might therefore lack some entries!
664 }
665
666 }
667 };
668
669 template<class M, class I>
670 typename MatrixRowGatherScatter<M,I>::ColIter MatrixRowGatherScatter<M,I>::col;
671
672 template<class M, class I>
673 typename MatrixRowGatherScatter<M,I>::Data MatrixRowGatherScatter<M,I>::datastore;
674
675 template<class M, class I>
676 typename MatrixRowGatherScatter<M,I>::GlobalIndex MatrixRowGatherScatter<M,I>::numlimits;
677
678 template<typename M, typename C>
679 void redistributeSparsityPattern(M& origMatrix, M& newMatrix, C& origComm, C& newComm,
680 RedistributeInformation<C>& ri)
681 {
682 typename C::CopySet copyflags;
683 typename C::OwnerSet ownerflags;
684 typedef typename C::ParallelIndexSet IndexSet;
685 typedef RedistributeInformation<C> RI;
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);
689
690 // get owner rowsizes
691 CommMatrixRowSize<M,RI> commRowSize(origMatrix, ri);
692 ri.template redistribute<MatrixRowSizeGatherScatter<M,IndexSet,RI> >(commRowSize,commRowSize);
693
694 origComm.buildGlobalLookup();
695
696 for (std::size_t i=0; i < newComm.indexSet().size(); i++) {
697 rowsize[i] = ri.getRowSize(i);
698 }
699 // get sparsity pattern from owner rows
700 CommMatrixSparsityPattern<M,IndexSet>
701 origsp(origMatrix, origComm.globalLookup(), newComm.indexSet());
702 CommMatrixSparsityPattern<M,IndexSet>
703 newsp(origMatrix, origComm.globalLookup(), newComm.indexSet(), rowsize);
704
705 ri.template redistribute<MatrixSparsityPatternGatherScatter<M,IndexSet> >(origsp,newsp);
706
707 // build copy to owner interface to get missing matrix values for novlp case
708 if (origComm.getSolverCategory() == SolverCategory::nonoverlapping) {
709 RemoteIndices<IndexSet> *ris = new RemoteIndices<IndexSet>(origComm.indexSet(),
710 newComm.indexSet(),
711 origComm.communicator());
712 ris->template rebuild<true>();
713
714 ri.getInterface().free();
715 ri.getInterface().build(*ris,copyflags,ownerflags);
716
717 // get copy rowsizes
718 CommMatrixRowSize<M,RI> commRowSize_copy(origMatrix, ri);
719 ri.template redistribute<MatrixCopyRowSizeGatherScatter<M,IndexSet,RI> >(commRowSize_copy,
720 commRowSize_copy);
721
722 for (std::size_t i=0; i < newComm.indexSet().size(); i++) {
723 copyrowsize[i] = ri.getCopyRowSize(i);
724 }
725 //get copy rowsizes for sender
726 ri.redistributeBackward(backwardscopyrowsize,copyrowsize);
727 for (std::size_t i=0; i < origComm.indexSet().size(); i++) {
728 ri.getBackwardsCopyRowSize(i) = backwardscopyrowsize[i];
729 }
730
731 // get sparsity pattern from copy rows
732 CommMatrixSparsityPattern<M,IndexSet> origsp_copy(origMatrix,
733 origComm.globalLookup(),
734 newComm.indexSet(),
735 backwardscopyrowsize);
736 CommMatrixSparsityPattern<M,IndexSet> newsp_copy(origMatrix, origComm.globalLookup(),
737 newComm.indexSet(), copyrowsize);
738 ri.template redistribute<MatrixSparsityPatternGatherScatter<M,IndexSet> >(origsp_copy,
739 newsp_copy);
740
741 newsp.completeSparsityPattern(newsp_copy.sparsity);
742 newsp.storeSparsityPattern(newMatrix);
743 }
744 else
745 newsp.storeSparsityPattern(newMatrix);
746
747#ifdef DUNE_ISTL_WITH_CHECKING
748 // Check for symmetry
749 int ret=0;
750 typedef typename M::ConstRowIterator RIter;
751 for(RIter row=newMatrix.begin(), rend=newMatrix.end(); row != rend; ++row) {
752 typedef typename M::ConstColIterator CIter;
753 for(CIter col=row->begin(), cend=row->end(); col!=cend; ++col)
754 {
755 try{
756 newMatrix[col.index()][row.index()];
757 }catch(Dune::ISTLError e) {
758 std::cerr<<newComm.communicator().rank()<<": entry ("
759 <<col.index()<<","<<row.index()<<") missing! for symmetry!"<<std::endl;
760 ret=1;
761
762 }
763
764 }
765 }
766
767 if(ret)
768 DUNE_THROW(ISTLError, "Matrix not symmetric!");
769#endif
770 }
771
772 template<typename M, typename C>
773 void redistributeMatrixEntries(M& origMatrix, M& newMatrix, C& origComm, C& newComm,
774 RedistributeInformation<C>& ri)
775 {
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);
781
782 for (std::size_t i=0; i < newComm.indexSet().size(); i++) {
783 rowsize[i] = ri.getRowSize(i);
784 if (origComm.getSolverCategory() == SolverCategory::nonoverlapping) {
785 copyrowsize[i] = ri.getCopyRowSize(i);
786 }
787 }
788
789 for (std::size_t i=0; i < origComm.indexSet().size(); i++)
790 if (origComm.getSolverCategory() == SolverCategory::nonoverlapping)
791 backwardscopyrowsize[i] = ri.getBackwardsCopyRowSize(i);
792
793
794 if (origComm.getSolverCategory() == SolverCategory::nonoverlapping) {
795 // fill sparsity pattern from copy rows
796 CommMatrixRow<M,IndexSet> origrow_copy(origMatrix, origComm.globalLookup(),
797 newComm.indexSet(), backwardscopyrowsize);
798 CommMatrixRow<M,IndexSet> newrow_copy(newMatrix, origComm.globalLookup(),
799 newComm.indexSet(),copyrowsize);
800 ri.template redistribute<MatrixRowGatherScatter<M,IndexSet> >(origrow_copy,
801 newrow_copy);
802 ri.getInterface().free();
803 RemoteIndices<IndexSet> *ris = new RemoteIndices<IndexSet>(origComm.indexSet(),
804 newComm.indexSet(),
805 origComm.communicator());
806 ris->template rebuild<true>();
807 ri.getInterface().build(*ris,ownerflags,ownerflags);
808 }
809
810 CommMatrixRow<M,IndexSet>
811 origrow(origMatrix, origComm.globalLookup(), newComm.indexSet());
812 CommMatrixRow<M,IndexSet>
813 newrow(newMatrix, origComm.globalLookup(), newComm.indexSet(),rowsize);
814 ri.template redistribute<MatrixRowGatherScatter<M,IndexSet> >(origrow,newrow);
815 if (origComm.getSolverCategory() != static_cast<int>(SolverCategory::nonoverlapping))
816 newrow.setOverlapRowsToDirichlet();
817 if(newMatrix.N()>0&&newMatrix.N()<20)
818 printmatrix(std::cout, newMatrix, "redist", "row");
819 }
820
837 template<typename M, typename C>
838 void redistributeMatrix(M& origMatrix, M& newMatrix, C& origComm, C& newComm,
839 RedistributeInformation<C>& ri)
840 {
841 ri.setNoRows(newComm.indexSet().size());
842 ri.setNoCopyRows(newComm.indexSet().size());
843 ri.setNoBackwardsCopyRows(origComm.indexSet().size());
844 redistributeSparsityPattern(origMatrix, newMatrix, origComm, newComm, ri);
845 redistributeMatrixEntries(origMatrix, newMatrix, origComm, newComm, ri);
846 }
847#endif
848}
849#endif
Dune::ConstArrayListIterator
A constant random access iterator for the Dune::ArrayList class.
Definition: arraylist.hh:380
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:507
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:280
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:244
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:258
indexset.hh
Provides a map between global and local indices.
Dune
Dune namespace.
Definition: alignment.hh:14
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:838
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:232
Dune::CommMatrixRowSize::CommMatrixRowSize
CommMatrixRowSize(const M &m_, RI &rowsize_)
Constructor.
Definition: matrixredistribute.hh:242
Dune::CommMatrixRow
Utility class for comunicating the matrix entries.
Definition: matrixredistribute.hh:412
Dune::CommMatrixRow::rowsize
std::vector< size_t > * rowsize
row size information for the receiving side.
Definition: matrixredistribute.hh:466
Dune::CommMatrixRow::matrix
M & matrix
The matrix to communicate the values of.
Definition: matrixredistribute.hh:460
Dune::CommMatrixRow::CommMatrixRow
CommMatrixRow(M &m_, const Dune::GlobalLookupIndexSet< I > &idxset_, const I &aggidxset_, std::vector< size_t > &rowsize_)
Constructor.
Definition: matrixredistribute.hh:428
Dune::CommMatrixRow::idxset
const Dune::GlobalLookupIndexSet< I > & idxset
Index set for the original matrix.
Definition: matrixredistribute.hh:462
Dune::CommMatrixRow::setOverlapRowsToDirichlet
void setOverlapRowsToDirichlet()
Sets the non-owner rows correctly as Dirichlet boundaries.
Definition: matrixredistribute.hh:437
Dune::CommMatrixRow::aggidxset
const I & aggidxset
Index set for the redistributed matrix.
Definition: matrixredistribute.hh:464
Dune::CommMatrixRow::CommMatrixRow
CommMatrixRow(M &m_, const Dune::GlobalLookupIndexSet< I > &idxset_, const I &aggidxset_)
Constructor.
Definition: matrixredistribute.hh:421
Dune::CommMatrixSparsityPattern
Utility class to communicate and build the sparsity pattern of a redistributed matrix.
Definition: matrixredistribute.hh:261
Dune::CommMatrixSparsityPattern::storeSparsityPattern
void storeSparsityPattern(M &m)
Creates and stores the sparsity pattern of the redistributed matrix.
Definition: matrixredistribute.hh:292
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:356
Dune::CommMatrixSparsityPattern::CommMatrixSparsityPattern
CommMatrixSparsityPattern(const M &m_, const Dune::GlobalLookupIndexSet< I > &idxset_, const I &aggidxset_)
Constructor for the original side.
Definition: matrixredistribute.hh:270
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:281
Dune::CommPolicy
Default policy used for communicating an indexed type.
Definition: communicator.hh:121
Dune::CommPolicy::IndexedType
V::value_type IndexedType
The type we get at each index with operator[].
Definition: communicator.hh:140
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:146
Dune::CommPolicy::Type
V Type
The type the policy is for.
Definition: communicator.hh:133
Dune::SolverCategory::nonoverlapping
@ nonoverlapping
Category for on overlapping solvers.
Definition: solvercategory.hh:24
Dune::VariableSize
Flag for marking indexed data structures where the data at each index may be a variable multiple of a...
Definition: communicator.hh:111
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 (Jul 15, 22:36, 2024)