DUNE PDELab (2.8)

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