DUNE-FEM (unstable)

matrixredistribute.hh
Go to the documentation of this file.
1// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
2// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
3// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
4// vi: set et ts=4 sw=2 sts=2:
5#ifndef DUNE_ISTL_MATRIXREDISTRIBUTE_HH
6#define DUNE_ISTL_MATRIXREDISTRIBUTE_HH
7#include <memory>
8#include "repartition.hh"
12#include <dune/istl/paamg/pinfo.hh>
18namespace Dune
19{
20 template<typename T>
21 struct RedistributeInformation
22 {
23 bool isSetup() const
24 {
25 return false;
26 }
27 template<class D>
28 void redistribute([[maybe_unused]] const D& from, [[maybe_unused]] D& to) const
29 {}
30
31 template<class D>
32 void redistributeBackward([[maybe_unused]] D& from, [[maybe_unused]]const D& to) const
33 {}
34
35 void resetSetup()
36 {}
37
38 void setNoRows([[maybe_unused]] std::size_t size)
39 {}
40
41 void setNoCopyRows([[maybe_unused]] std::size_t size)
42 {}
43
44 void setNoBackwardsCopyRows([[maybe_unused]] std::size_t size)
45 {}
46
47 std::size_t getRowSize([[maybe_unused]] std::size_t index) const
48 {
49 return -1;
50 }
51
52 std::size_t getCopyRowSize([[maybe_unused]] std::size_t index) const
53 {
54 return -1;
55 }
56
57 std::size_t getBackwardsCopyRowSize([[maybe_unused]] std::size_t index) const
58 {
59 return -1;
60 }
61
62 };
63
64#if HAVE_MPI
65 template<typename T, typename T1>
66 class RedistributeInformation<OwnerOverlapCopyCommunication<T,T1> >
67 {
68 public:
69 typedef OwnerOverlapCopyCommunication<T,T1> Comm;
70
71 RedistributeInformation()
72 : interface(), setup_(false)
73 {}
74
75 RedistributeInterface& getInterface()
76 {
77 return interface;
78 }
79 template<typename IS>
80 void checkInterface(const IS& source,
81 const IS& target, MPI_Comm comm)
82 {
83 auto ri = std::make_unique<RemoteIndices<IS> >(source, target, comm);
84 ri->template rebuild<true>();
85 Interface inf;
86 typename OwnerOverlapCopyCommunication<int>::OwnerSet flags;
87 int rank;
88 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
89 inf.free();
90 inf.build(*ri, flags, flags);
91
92
93#ifdef DEBUG_REPART
94 if(inf!=interface) {
95
96 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
97 if(rank==0)
98 std::cout<<"Interfaces do not match!"<<std::endl;
99 std::cout<<rank<<": redist interface new :"<<inf<<std::endl;
100 std::cout<<rank<<": redist interface :"<<interface<<std::endl;
101
102 throw "autsch!";
103 }
104#endif
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
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 for (unsigned int i = 0; i != sparsity.size(); ++i) {
343 if (add_sparsity[i].size() != 0) {
344 typedef std::set<size_type> Set;
345 Set tmp_set;
346 std::insert_iterator<Set> tmp_insert (tmp_set, tmp_set.begin());
347 std::set_union(add_sparsity[i].begin(), add_sparsity[i].end(),
348 sparsity[i].begin(), sparsity[i].end(), tmp_insert);
349 sparsity[i].swap(tmp_set);
350 }
351 }
352 }
353
354 const M& matrix;
355 typedef Dune::GlobalLookupIndexSet<I> LookupIndexSet;
356 const Dune::GlobalLookupIndexSet<I>& idxset;
357 const I& aggidxset;
358 std::vector<std::set<size_type> > sparsity;
359 const std::vector<size_type>* rowsize;
360 };
361
362 template<class M, class I>
363 struct CommPolicy<CommMatrixSparsityPattern<M,I> >
364 {
365 typedef CommMatrixSparsityPattern<M,I> Type;
366
371 typedef typename I::GlobalIndex IndexedType;
372
374 typedef VariableSize IndexedTypeFlag;
375
376 static typename M::size_type getSize(const Type& t, std::size_t i)
377 {
378 if(!t.rowsize)
379 return t.matrix[i].size();
380 else
381 {
382 assert((*t.rowsize)[i]>0);
383 return (*t.rowsize)[i];
384 }
385 }
386 };
387
394 template<class M, class I>
396 {
405 CommMatrixRow(M& m_, const Dune::GlobalLookupIndexSet<I>& idxset_, const I& aggidxset_)
406 : matrix(m_), idxset(idxset_), aggidxset(aggidxset_), rowsize()
407 {}
408
412 CommMatrixRow(M& m_, const Dune::GlobalLookupIndexSet<I>& idxset_, const I& aggidxset_,
413 std::vector<size_t>& rowsize_)
414 : matrix(m_), idxset(idxset_), aggidxset(aggidxset_), rowsize(&rowsize_)
415 {}
422 {
425
426 for(Iter i= aggidxset.begin(), end=aggidxset.end(); i!=end; ++i)
427 if(!OwnerSet::contains(i->local().attribute())) {
428 // Set to Dirchlet
429 typedef typename M::ColIterator CIter;
430 for(CIter c=matrix[i->local()].begin(), cend= matrix[i->local()].end();
431 c!= cend; ++c)
432 {
433 *c=0;
434 if(c.index()==i->local()) {
435 auto setDiagonal = [](auto&& scalarOrMatrix, const auto& value) {
436 auto&& matrixView = Dune::Impl::asMatrix(scalarOrMatrix);
437 for (auto rowIt = matrixView.begin(); rowIt != matrixView.end(); ++rowIt)
438 (*rowIt)[rowIt.index()] = value;
439 };
440 setDiagonal(*c, 1);
441 }
442 }
443 }
444 }
450 const I& aggidxset;
452 std::vector<size_t>* rowsize; // row sizes differ from sender side in overlap!
453 };
454
455 template<class M, class I>
456 struct CommPolicy<CommMatrixRow<M,I> >
457 {
458 typedef CommMatrixRow<M,I> Type;
459
464 typedef std::pair<typename I::GlobalIndex,typename M::block_type> IndexedType;
465
468
469 static std::size_t getSize(const Type& t, std::size_t i)
470 {
471 if(!t.rowsize)
472 return t.matrix[i].size();
473 else
474 {
475 assert((*t.rowsize)[i]>0);
476 return (*t.rowsize)[i];
477 }
478 }
479 };
480
481 template<class M, class I, class RI>
482 struct MatrixRowSizeGatherScatter
483 {
484 typedef CommMatrixRowSize<M,RI> Container;
485
486 static const typename M::size_type gather(const Container& cont, std::size_t i)
487 {
488 return cont.matrix[i].size();
489 }
490 static void scatter(Container& cont, const typename M::size_type& rowsize, std::size_t i)
491 {
492 assert(rowsize);
493 cont.rowsize.getRowSize(i)=rowsize;
494 }
495
496 };
497
498 template<class M, class I, class RI>
499 struct MatrixCopyRowSizeGatherScatter
500 {
501 typedef CommMatrixRowSize<M,RI> Container;
502
503 static const typename M::size_type gather(const Container& cont, std::size_t i)
504 {
505 return cont.matrix[i].size();
506 }
507 static void scatter(Container& cont, const typename M::size_type& rowsize, std::size_t i)
508 {
509 assert(rowsize);
510 if (rowsize > cont.rowsize.getCopyRowSize(i))
511 cont.rowsize.getCopyRowSize(i)=rowsize;
512 }
513
514 };
515
516 template<class M, class I>
517 struct MatrixSparsityPatternGatherScatter
518 {
519 typedef typename I::GlobalIndex GlobalIndex;
520 typedef CommMatrixSparsityPattern<M,I> Container;
521 typedef typename M::ConstColIterator ColIter;
522
523 static ColIter col;
524 static GlobalIndex numlimits;
525
526 static const GlobalIndex& gather(const Container& cont, std::size_t i, std::size_t j)
527 {
528 if(j==0)
529 col=cont.matrix[i].begin();
530 else if (col!=cont.matrix[i].end())
531 ++col;
532
533 //copy communication: different row sizes for copy rows with the same global index
534 //are possible. If all values of current matrix row are sent, send
535 //std::numeric_limits<GlobalIndex>::max()
536 //and receiver will ignore it.
537 if (col==cont.matrix[i].end()) {
539 return numlimits;
540 }
541 else {
542 const typename I::IndexPair* index=cont.idxset.pair(col.index());
543 assert(index);
544 // Only send index if col is no ghost
545 if ( index->local().attribute() != 2)
546 return index->global();
547 else {
549 return numlimits;
550 }
551 }
552 }
553 static void scatter(Container& cont, const GlobalIndex& gi, std::size_t i, [[maybe_unused]] std::size_t j)
554 {
555 try{
557 const typename I::IndexPair& ip=cont.aggidxset.at(gi);
558 assert(ip.global()==gi);
559 std::size_t column = ip.local();
560 cont.sparsity[i].insert(column);
561
563 if(!OwnerSet::contains(ip.local().attribute()))
564 // preserve symmetry for overlap
565 cont.sparsity[column].insert(i);
566 }
567 }
568 catch(const Dune::RangeError&) {
569 // Entry not present in the new index set. Ignore!
570#ifdef DEBUG_REPART
571 typedef typename Container::LookupIndexSet GlobalLookup;
572 typedef typename GlobalLookup::IndexPair IndexPair;
574
575 GlobalLookup lookup(cont.aggidxset);
576 const IndexPair* pi=lookup.pair(i);
577 assert(pi);
578 if(OwnerSet::contains(pi->local().attribute())) {
579 int rank;
580 MPI_Comm_rank(MPI_COMM_WORLD,&rank);
581 std::cout<<rank<<cont.aggidxset<<std::endl;
582 std::cout<<rank<<": row "<<i<<" (global="<<gi <<") not in index set for owner index "<<pi->global()<<std::endl;
583 throw;
584 }
585#endif
586 }
587 }
588
589 };
590 template<class M, class I>
591 typename MatrixSparsityPatternGatherScatter<M,I>::ColIter MatrixSparsityPatternGatherScatter<M,I>::col;
592
593 template<class M, class I>
594 typename MatrixSparsityPatternGatherScatter<M,I>::GlobalIndex MatrixSparsityPatternGatherScatter<M,I>::numlimits;
595
596
597 template<class M, class I>
598 struct MatrixRowGatherScatter
599 {
600 typedef typename I::GlobalIndex GlobalIndex;
601 typedef CommMatrixRow<M,I> Container;
602 typedef typename M::ConstColIterator ColIter;
603 typedef typename std::pair<GlobalIndex,typename M::block_type> Data;
604 static ColIter col;
605 static Data datastore;
606 static GlobalIndex numlimits;
607
608 static const Data& gather(const Container& cont, std::size_t i, std::size_t j)
609 {
610 if(j==0)
611 col=cont.matrix[i].begin();
612 else if (col!=cont.matrix[i].end())
613 ++col;
614 // copy communication: different row sizes for copy rows with the same global index
615 // are possible. If all values of current matrix row are sent, send
616 // std::numeric_limits<GlobalIndex>::max()
617 // and receiver will ignore it.
618 if (col==cont.matrix[i].end()) {
620 datastore = Data(numlimits,*col);
621 return datastore;
622 }
623 else {
624 // convert local column index to global index
625 const typename I::IndexPair* index=cont.idxset.pair(col.index());
626 assert(index);
627 // Store the data to prevent reference to temporary
628 // Only send index if col is no ghost
629 if ( index->local().attribute() != 2)
630 datastore = Data(index->global(),*col);
631 else {
633 datastore = Data(numlimits,*col);
634 }
635 return datastore;
636 }
637 }
638 static void scatter(Container& cont, const Data& data, std::size_t i, [[maybe_unused]] std::size_t j)
639 {
640 try{
641 if (data.first != std::numeric_limits<GlobalIndex>::max()) {
642 typename M::size_type column=cont.aggidxset.at(data.first).local();
643 cont.matrix[i][column]=data.second;
644 }
645 }
646 catch(const Dune::RangeError&) {
647 // This an overlap row and might therefore lack some entries!
648 }
649
650 }
651 };
652
653 template<class M, class I>
654 typename MatrixRowGatherScatter<M,I>::ColIter MatrixRowGatherScatter<M,I>::col;
655
656 template<class M, class I>
657 typename MatrixRowGatherScatter<M,I>::Data MatrixRowGatherScatter<M,I>::datastore;
658
659 template<class M, class I>
660 typename MatrixRowGatherScatter<M,I>::GlobalIndex MatrixRowGatherScatter<M,I>::numlimits;
661
662 template<typename M, typename C>
663 void redistributeSparsityPattern(M& origMatrix, M& newMatrix, C& origComm, C& newComm,
664 RedistributeInformation<C>& ri)
665 {
666 typename C::CopySet copyflags;
667 typename C::OwnerSet ownerflags;
668 typedef typename C::ParallelIndexSet IndexSet;
669 typedef RedistributeInformation<C> RI;
670 std::vector<typename M::size_type> rowsize(newComm.indexSet().size(), 0);
671 std::vector<typename M::size_type> copyrowsize(newComm.indexSet().size(), 0);
672 std::vector<typename M::size_type> backwardscopyrowsize(origComm.indexSet().size(), 0);
673
674 // get owner rowsizes
675 CommMatrixRowSize<M,RI> commRowSize(origMatrix, ri);
676 ri.template redistribute<MatrixRowSizeGatherScatter<M,IndexSet,RI> >(commRowSize,commRowSize);
677
678 origComm.buildGlobalLookup();
679
680 for (std::size_t i=0; i < newComm.indexSet().size(); i++) {
681 rowsize[i] = ri.getRowSize(i);
682 }
683 // get sparsity pattern from owner rows
684 CommMatrixSparsityPattern<M,IndexSet>
685 origsp(origMatrix, origComm.globalLookup(), newComm.indexSet());
686 CommMatrixSparsityPattern<M,IndexSet>
687 newsp(origMatrix, origComm.globalLookup(), newComm.indexSet(), rowsize);
688
689 ri.template redistribute<MatrixSparsityPatternGatherScatter<M,IndexSet> >(origsp,newsp);
690
691 // build copy to owner interface to get missing matrix values for novlp case
693 RemoteIndices<IndexSet> *ris = new RemoteIndices<IndexSet>(origComm.indexSet(),
694 newComm.indexSet(),
695 origComm.communicator());
696 ris->template rebuild<true>();
697
698 ri.getInterface().free();
699 ri.getInterface().build(*ris,copyflags,ownerflags);
700
701 // get copy rowsizes
702 CommMatrixRowSize<M,RI> commRowSize_copy(origMatrix, ri);
703 ri.template redistribute<MatrixCopyRowSizeGatherScatter<M,IndexSet,RI> >(commRowSize_copy,
704 commRowSize_copy);
705
706 for (std::size_t i=0; i < newComm.indexSet().size(); i++) {
707 copyrowsize[i] = ri.getCopyRowSize(i);
708 }
709 //get copy rowsizes for sender
710 ri.redistributeBackward(backwardscopyrowsize,copyrowsize);
711 for (std::size_t i=0; i < origComm.indexSet().size(); i++) {
712 ri.getBackwardsCopyRowSize(i) = backwardscopyrowsize[i];
713 }
714
715 // get sparsity pattern from copy rows
716 CommMatrixSparsityPattern<M,IndexSet> origsp_copy(origMatrix,
717 origComm.globalLookup(),
718 newComm.indexSet(),
719 backwardscopyrowsize);
720 CommMatrixSparsityPattern<M,IndexSet> newsp_copy(origMatrix, origComm.globalLookup(),
721 newComm.indexSet(), copyrowsize);
722 ri.template redistribute<MatrixSparsityPatternGatherScatter<M,IndexSet> >(origsp_copy,
723 newsp_copy);
724
725 newsp.completeSparsityPattern(newsp_copy.sparsity);
726 newsp.storeSparsityPattern(newMatrix);
727 }
728 else
729 newsp.storeSparsityPattern(newMatrix);
730
731#ifdef DUNE_ISTL_WITH_CHECKING
732 // Check for symmetry
733 int ret=0;
734 typedef typename M::ConstRowIterator RIter;
735 for(RIter row=newMatrix.begin(), rend=newMatrix.end(); row != rend; ++row) {
736 typedef typename M::ConstColIterator CIter;
737 for(CIter col=row->begin(), cend=row->end(); col!=cend; ++col)
738 {
739 try{
740 newMatrix[col.index()][row.index()];
741 }catch(const Dune::ISTLError&) {
742 std::cerr<<newComm.communicator().rank()<<": entry ("
743 <<col.index()<<","<<row.index()<<") missing! for symmetry!"<<std::endl;
744 ret=1;
745
746 }
747
748 }
749 }
750
751 if(ret)
752 DUNE_THROW(ISTLError, "Matrix not symmetric!");
753#endif
754 }
755
756 template<typename M, typename C>
757 void redistributeMatrixEntries(M& origMatrix, M& newMatrix, C& origComm, C& newComm,
758 RedistributeInformation<C>& ri)
759 {
760 typedef typename C::ParallelIndexSet IndexSet;
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);
769 copyrowsize[i] = ri.getCopyRowSize(i);
770 }
771 }
772
773 for (std::size_t i=0; i < origComm.indexSet().size(); i++)
775 backwardscopyrowsize[i] = ri.getBackwardsCopyRowSize(i);
776
777
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
794 CommMatrixRow<M,IndexSet>
795 origrow(origMatrix, origComm.globalLookup(), newComm.indexSet());
796 CommMatrixRow<M,IndexSet>
797 newrow(newMatrix, origComm.globalLookup(), newComm.indexSet(),rowsize);
798 ri.template redistribute<MatrixRowGatherScatter<M,IndexSet> >(origrow,newrow);
799 if (SolverCategory::category(origComm) != static_cast<int>(SolverCategory::nonoverlapping))
800 newrow.setOverlapRowsToDirichlet();
801 }
802
819 template<typename M, typename C>
820 void redistributeMatrix(M& origMatrix, M& newMatrix, C& origComm, C& newComm,
821 RedistributeInformation<C>& ri)
822 {
823 ri.setNoRows(newComm.indexSet().size());
824 ri.setNoCopyRows(newComm.indexSet().size());
825 ri.setNoBackwardsCopyRows(origComm.indexSet().size());
826 redistributeSparsityPattern(origMatrix, newMatrix, origComm, newComm, ri);
827 redistributeMatrixEntries(origMatrix, newMatrix, origComm, newComm, ri);
828 }
829#endif
830
831template<typename M>
832 void redistributeMatrixEntries(M& origMatrix, M& newMatrix,
833 Dune::Amg::SequentialInformation& origComm,
834 Dune::Amg::SequentialInformation& newComm,
835 RedistributeInformation<Dune::Amg::SequentialInformation>& ri)
836 {
837 DUNE_THROW(InvalidStateException, "Trying to redistribute in sequential program!");
838 }
839 template<typename M>
840 void redistributeMatrix(M& origMatrix, M& newMatrix,
841 Dune::Amg::SequentialInformation& origComm,
842 Dune::Amg::SequentialInformation& newComm,
843 RedistributeInformation<Dune::Amg::SequentialInformation>& ri)
844 {
845 DUNE_THROW(InvalidStateException, "Trying to redistribute in sequential program!");
846 }
847}
848#endif
A constant random access iterator for the Dune::ArrayList class.
Definition: arraylist.hh:368
A set consisting only of one item.
Definition: enumset.hh:61
Decorates an index set with the possibility to find a global index that is mapped to a specific local...
Definition: indexset.hh:507
derive error class from the base class in common
Definition: istlexception.hh:19
Default exception class for range errors.
Definition: exceptions.hh:254
A few common exception classes.
Provides a map between global and local indices.
Classes providing communication interfaces for overlapping Schwarz methods.
const IndexPair * pair(const std::size_t &local) const
Get the index pair corresponding to a local index.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:484
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
void redistributeMatrix(M &origMatrix, M &newMatrix, C &origComm, C &newComm, RedistributeInformation< C > &ri)
Redistribute a matrix according to given domain decompositions.
Definition: matrixredistribute.hh:820
Functionality for redistributing a parallel index set using graph partitioning.
Utility class to communicate and set the row sizes of a redistributed matrix.
Definition: matrixredistribute.hh:216
CommMatrixRowSize(const M &m_, RI &rowsize_)
Constructor.
Definition: matrixredistribute.hh:226
Utility class for comunicating the matrix entries.
Definition: matrixredistribute.hh:396
std::vector< size_t > * rowsize
row size information for the receiving side.
Definition: matrixredistribute.hh:452
M & matrix
The matrix to communicate the values of.
Definition: matrixredistribute.hh:446
CommMatrixRow(M &m_, const Dune::GlobalLookupIndexSet< I > &idxset_, const I &aggidxset_, std::vector< size_t > &rowsize_)
Constructor.
Definition: matrixredistribute.hh:412
const Dune::GlobalLookupIndexSet< I > & idxset
Index set for the original matrix.
Definition: matrixredistribute.hh:448
void setOverlapRowsToDirichlet()
Sets the non-owner rows correctly as Dirichlet boundaries.
Definition: matrixredistribute.hh:421
const I & aggidxset
Index set for the redistributed matrix.
Definition: matrixredistribute.hh:450
CommMatrixRow(M &m_, const Dune::GlobalLookupIndexSet< I > &idxset_, const I &aggidxset_)
Constructor.
Definition: matrixredistribute.hh:405
Utility class to communicate and build the sparsity pattern of a redistributed matrix.
Definition: matrixredistribute.hh:245
void storeSparsityPattern(M &m)
Creates and stores the sparsity pattern of the redistributed matrix.
Definition: matrixredistribute.hh:276
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:340
CommMatrixSparsityPattern(const M &m_, const Dune::GlobalLookupIndexSet< I > &idxset_, const I &aggidxset_)
Constructor for the original side.
Definition: matrixredistribute.hh:254
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:265
Default policy used for communicating an indexed type.
Definition: communicator.hh:128
V::value_type IndexedType
The type we get at each index with operator[].
Definition: communicator.hh:147
static int getSize(const V &, int index)
Get the number of primitive elements at that index.
SizeOne IndexedTypeFlag
Whether the indexed type has variable size or there is always one value at each index.
Definition: communicator.hh:153
V Type
The type the policy is for.
Definition: communicator.hh:140
@ nonoverlapping
Category for non-overlapping solvers.
Definition: solvercategory.hh:27
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:34
Flag for marking indexed data structures where the data at each index may be a variable multiple of a...
Definition: communicator.hh:118
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 12, 23:30, 2024)