Dune Core Modules (2.8.0)

galerkin.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_GALERKIN_HH
4#define DUNE_GALERKIN_HH
5
6#include "aggregates.hh"
7#include "pinfo.hh"
10#include <set>
11#include <limits>
12#include <algorithm>
13
14namespace Dune
15{
16 namespace Amg
17 {
29 template<class T>
30 struct OverlapVertex
31 {
35 typedef T Aggregate;
36
40 typedef T Vertex;
41
46
51 };
52
53
54
59 template<class M>
61 {
62 public:
68 SparsityBuilder(M& matrix);
69
70 void insert(const typename M::size_type& index);
71
72 void operator++();
73
74 std::size_t minRowSize();
75
76 std::size_t maxRowSize();
77
78 std::size_t sumRowSize();
79 std::size_t index()
80 {
81 return row_.index();
82 }
83 private:
85 typename M::CreateIterator row_;
87 std::size_t minRowSize_;
89 std::size_t maxRowSize_;
90 std::size_t sumRowSize_;
91#ifdef DUNE_ISTL_WITH_CHECKING
92 bool diagonalInserted;
93#endif
94 };
95
96 class BaseGalerkinProduct
97 {
98 public:
107 template<class M, class V, class I, class O>
108 void calculate(const M& fine, const AggregatesMap<V>& aggregates, M& coarse,
109 const I& pinfo, const O& copy);
110
111 };
112
113 template<class T>
114 class GalerkinProduct
115 : public BaseGalerkinProduct
116 {
117 public:
118 typedef T ParallelInformation;
119
129 template<class G, class V, class Set>
130 typename G::MutableMatrix* build(G& fineGraph, V& visitedMap,
131 const ParallelInformation& pinfo,
133 const typename G::Matrix::size_type& size,
134 const Set& copy);
135 private:
136
143 template<class G, class I, class Set>
144 const OverlapVertex<typename G::VertexDescriptor>*
145 buildOverlapVertices(const G& graph, const I& pinfo,
147 const Set& overlap,
148 std::size_t& overlapCount);
149
150 template<class A>
151 struct OVLess
152 {
153 bool operator()(const OverlapVertex<A>& o1, const OverlapVertex<A>& o2)
154 {
155 return *o1.aggregate < *o2.aggregate;
156 }
157 };
158 };
159
160 template<>
161 class GalerkinProduct<SequentialInformation>
162 : public BaseGalerkinProduct
163 {
164 public:
174 template<class G, class V, class Set>
175 typename G::MutableMatrix* build(G& fineGraph, V& visitedMap,
176 const SequentialInformation& pinfo,
177 const AggregatesMap<typename G::VertexDescriptor>& aggregates,
178 const typename G::Matrix::size_type& size,
179 const Set& copy);
180 };
181
182 struct BaseConnectivityConstructor
183 {
184 template<class R, class G, class V>
185 static void constructOverlapConnectivity(R& row, G& graph, V& visitedMap,
186 const AggregatesMap<typename G::VertexDescriptor>& aggregates,
187 const OverlapVertex<typename G::VertexDescriptor>*& seed,
188 const OverlapVertex<typename G::VertexDescriptor>* overlapEnd);
189
193 template<class R, class G, class V>
194 static void constructNonOverlapConnectivity(R& row, G& graph, V& visitedMap,
195 const AggregatesMap<typename G::VertexDescriptor>& aggregates,
196 const typename G::VertexDescriptor& seed);
197
198
202 template<class G, class S, class V>
204 {
205 public:
209 typedef G Graph;
213 typedef typename Graph::ConstEdgeIterator ConstEdgeIterator;
214
218 typedef S Set;
219
223 typedef V VisitedMap;
224
228 typedef typename Graph::VertexDescriptor Vertex;
229
237 ConnectedBuilder(const AggregatesMap<Vertex>& aggregates, Graph& graph,
238 VisitedMap& visitedMap, Set& connected);
239
244 void operator()(const ConstEdgeIterator& edge);
245
246 private:
250 const AggregatesMap<Vertex>& aggregates_;
251
252 Graph& graph_;
253
257 VisitedMap& visitedMap_;
258
262 Set& connected_;
263 };
264
265 };
266
267 template<class G, class T>
268 struct ConnectivityConstructor : public BaseConnectivityConstructor
269 {
270 typedef typename G::VertexDescriptor Vertex;
271
272 template<class V, class O, class R>
273 static void examine(G& graph,
274 V& visitedMap,
275 const T& pinfo,
276 const AggregatesMap<Vertex>& aggregates,
277 const O& overlap,
278 const OverlapVertex<Vertex>* overlapVertices,
279 const OverlapVertex<Vertex>* overlapEnd,
280 R& row);
281 };
282
283 template<class G>
284 struct ConnectivityConstructor<G,SequentialInformation> : public BaseConnectivityConstructor
285 {
286 typedef typename G::VertexDescriptor Vertex;
287
288 template<class V, class R>
289 static void examine(G& graph,
290 V& visitedMap,
291 const SequentialInformation& pinfo,
292 const AggregatesMap<Vertex>& aggregates,
293 R& row);
294 };
295
296 template<class T>
297 struct DirichletBoundarySetter
298 {
299 template<class M, class O>
300 static void set(M& coarse, const T& pinfo, const O& copy);
301 };
302
303 template<>
304 struct DirichletBoundarySetter<SequentialInformation>
305 {
306 template<class M, class O>
307 static void set(M& coarse, const SequentialInformation& pinfo, const O& copy);
308 };
309
310 template<class R, class G, class V>
311 void BaseConnectivityConstructor::constructNonOverlapConnectivity(R& row, G& graph, V& visitedMap,
313 const typename G::VertexDescriptor& seed)
314 {
315 assert(row.index()==aggregates[seed]);
316 row.insert(aggregates[seed]);
317 ConnectedBuilder<G,R,V> conBuilder(aggregates, graph, visitedMap, row);
318 typedef typename G::VertexDescriptor Vertex;
319 typedef std::allocator<Vertex> Allocator;
320 typedef SLList<Vertex,Allocator> VertexList;
321 typedef typename AggregatesMap<Vertex>::DummyEdgeVisitor DummyVisitor;
322 VertexList vlist;
323 DummyVisitor dummy;
324 aggregates.template breadthFirstSearch<true,false>(seed,aggregates[seed], graph, vlist, dummy,
325 conBuilder, visitedMap);
326 }
327
328 template<class R, class G, class V>
329 void BaseConnectivityConstructor::constructOverlapConnectivity(R& row, G& graph, V& visitedMap,
331 const OverlapVertex<typename G::VertexDescriptor>*& seed,
332 const OverlapVertex<typename G::VertexDescriptor>* overlapEnd)
333 {
334 ConnectedBuilder<G,R,V> conBuilder(aggregates, graph, visitedMap, row);
335 const typename G::VertexDescriptor aggregate=*seed->aggregate;
336
337 if (row.index()==*seed->aggregate) {
338 while(seed != overlapEnd && aggregate == *seed->aggregate) {
339 row.insert(*seed->aggregate);
340 // Walk over all neighbours and add them to the connected array.
341 visitNeighbours(graph, seed->vertex, conBuilder);
342 // Mark vertex as visited
343 put(visitedMap, seed->vertex, true);
344 ++seed;
345 }
346 }
347 }
348
349 template<class G, class S, class V>
351 Graph& graph, VisitedMap& visitedMap,
352 Set& connected)
353 : aggregates_(aggregates), graph_(graph), visitedMap_(visitedMap), connected_(connected)
354 {}
355
356 template<class G, class S, class V>
358 {
359 const Vertex& vertex = aggregates_[edge.target()];
362 connected_.insert(vertex);
363 }
364
365 template<class T>
366 template<class G, class I, class Set>
367 const OverlapVertex<typename G::VertexDescriptor>*
368 GalerkinProduct<T>::buildOverlapVertices(const G& graph, const I& pinfo,
370 const Set& overlap,
371 std::size_t& overlapCount)
372 {
373 // count the overlap vertices.
374 typedef typename G::ConstVertexIterator ConstIterator;
375 typedef typename I::GlobalLookupIndexSet GlobalLookup;
376 typedef typename GlobalLookup::IndexPair IndexPair;
377
378 const ConstIterator end = graph.end();
379 overlapCount = 0;
380
381 const GlobalLookup& lookup=pinfo.globalLookup();
382
383 for(ConstIterator vertex=graph.begin(); vertex != end; ++vertex) {
384 const IndexPair* pair = lookup.pair(*vertex);
385
386 if(pair!=0 && overlap.contains(pair->local().attribute()))
387 ++overlapCount;
388 }
389 // Allocate space
390 typedef typename G::VertexDescriptor Vertex;
391
392 OverlapVertex<Vertex>* overlapVertices = new OverlapVertex<Vertex>[overlapCount=0 ? 1 : overlapCount];
393 if(overlapCount==0)
394 return overlapVertices;
395
396 // Initialize them
397 overlapCount=0;
398 for(ConstIterator vertex=graph.begin(); vertex != end; ++vertex) {
399 const IndexPair* pair = lookup.pair(*vertex);
400
401 if(pair!=0 && overlap.contains(pair->local().attribute())) {
402 overlapVertices[overlapCount].aggregate = &aggregates[pair->local()];
403 overlapVertices[overlapCount].vertex = pair->local();
404 ++overlapCount;
405 }
406 }
407
408 dverb << overlapCount<<" overlap vertices"<<std::endl;
409
410 std::sort(overlapVertices, overlapVertices+overlapCount, OVLess<Vertex>());
411 // due to the sorting the isolated aggregates (to be skipped) are at the end.
412
413 return overlapVertices;
414 }
415
416 template<class G, class T>
417 template<class V, class O, class R>
418 void ConnectivityConstructor<G,T>::examine(G& graph,
419 V& visitedMap,
420 const T& pinfo,
421 const AggregatesMap<Vertex>& aggregates,
422 const O& overlap,
423 const OverlapVertex<Vertex>* overlapVertices,
424 const OverlapVertex<Vertex>* overlapEnd,
425 R& row)
426 {
427 typedef typename T::GlobalLookupIndexSet GlobalLookup;
428 const GlobalLookup& lookup = pinfo.globalLookup();
429
430 typedef typename G::VertexIterator VertexIterator;
431
432 VertexIterator vend=graph.end();
433
434#ifdef DUNE_ISTL_WITH_CHECKING
435 std::set<Vertex> examined;
436#endif
437
438 // The aggregates owned by the process have lower local indices
439 // then those not owned. We process them in the first pass.
440 // They represent the rows 0, 1, ..., n of the coarse matrix
441 for(VertexIterator vertex = graph.begin(); vertex != vend; ++vertex)
442 if(!get(visitedMap, *vertex)) {
443 // In the first pass we only process owner nodes
444 typedef typename GlobalLookup::IndexPair IndexPair;
445 const IndexPair* pair = lookup.pair(*vertex);
446 if(pair==0 || !overlap.contains(pair->local().attribute())) {
447#ifdef DUNE_ISTL_WITH_CHECKING
448 assert(examined.find(aggregates[*vertex])==examined.end());
449 examined.insert(aggregates[*vertex]);
450#endif
451 constructNonOverlapConnectivity(row, graph, visitedMap, aggregates, *vertex);
452
453 // only needed for ALU
454 // (ghosts with same global id as owners on the same process)
455 if (SolverCategory::category(pinfo) == static_cast<int>(SolverCategory::nonoverlapping)) {
456 if(overlapVertices != overlapEnd) {
457 if(*overlapVertices->aggregate!=AggregatesMap<Vertex>::ISOLATED) {
458 constructOverlapConnectivity(row, graph, visitedMap, aggregates, overlapVertices, overlapEnd);
459 }
460 else{
461 ++overlapVertices;
462 }
463 }
464 }
465 ++row;
466 }
467 }
468
469 dvverb<<"constructed "<<row.index()<<" non-overlapping rows"<<std::endl;
470
471 // Now come the aggregates not owned by use.
472 // They represent the rows n+1, ..., N
473 while(overlapVertices != overlapEnd)
474 if(*overlapVertices->aggregate!=AggregatesMap<Vertex>::ISOLATED) {
475
476#ifdef DUNE_ISTL_WITH_CHECKING
477 typedef typename GlobalLookup::IndexPair IndexPair;
478 const IndexPair* pair = lookup.pair(overlapVertices->vertex);
479 assert(pair!=0 && overlap.contains(pair->local().attribute()));
480 assert(examined.find(aggregates[overlapVertices->vertex])==examined.end());
481 examined.insert(aggregates[overlapVertices->vertex]);
482#endif
483 constructOverlapConnectivity(row, graph, visitedMap, aggregates, overlapVertices, overlapEnd);
484 ++row;
485 }else{
486 ++overlapVertices;
487 }
488 }
489
490 template<class G>
491 template<class V, class R>
492 void ConnectivityConstructor<G,SequentialInformation>::examine(G& graph,
493 V& visitedMap,
494 [[maybe_unused]] const SequentialInformation& pinfo,
495 const AggregatesMap<Vertex>& aggregates,
496 R& row)
497 {
498 typedef typename G::VertexIterator VertexIterator;
499
500 VertexIterator vend=graph.end();
501 for(VertexIterator vertex = graph.begin(); vertex != vend; ++vertex) {
502 if(!get(visitedMap, *vertex)) {
503 constructNonOverlapConnectivity(row, graph, visitedMap, aggregates, *vertex);
504 ++row;
505 }
506 }
507
508 }
509
510 template<class M>
512 : row_(matrix.createbegin()),
513 minRowSize_(std::numeric_limits<std::size_t>::max()),
514 maxRowSize_(0), sumRowSize_(0)
515 {
516#ifdef DUNE_ISTL_WITH_CHECKING
517 diagonalInserted = false;
518#endif
519 }
520 template<class M>
522 {
523 return maxRowSize_;
524 }
525 template<class M>
526 std::size_t SparsityBuilder<M>::minRowSize()
527 {
528 return minRowSize_;
529 }
530
531 template<class M>
532 std::size_t SparsityBuilder<M>::sumRowSize()
533 {
534 return sumRowSize_;
535 }
536 template<class M>
537 void SparsityBuilder<M>::operator++()
538 {
539 sumRowSize_ += row_.size();
540 minRowSize_=std::min(minRowSize_, row_.size());
541 maxRowSize_=std::max(maxRowSize_, row_.size());
542 ++row_;
543#ifdef DUNE_ISTL_WITH_CHECKING
544 assert(diagonalInserted);
545 diagonalInserted = false;
546#endif
547 }
548
549 template<class M>
550 void SparsityBuilder<M>::insert(const typename M::size_type& index)
551 {
552 row_.insert(index);
553#ifdef DUNE_ISTL_WITH_CHECKING
554 diagonalInserted = diagonalInserted || row_.index()==index;
555#endif
556 }
557
558 template<class T>
559 template<class G, class V, class Set>
560 typename G::MutableMatrix*
561 GalerkinProduct<T>::build(G& fineGraph, V& visitedMap,
562 const ParallelInformation& pinfo,
564 const typename G::Matrix::size_type& size,
565 const Set& overlap)
566 {
567 typedef OverlapVertex<typename G::VertexDescriptor> OverlapVertex;
568
569 std::size_t count;
570
571 const OverlapVertex* overlapVertices = buildOverlapVertices(fineGraph,
572 pinfo,
573 aggregates,
574 overlap,
575 count);
576 typedef typename G::MutableMatrix M;
577 M* coarseMatrix = new M(size, size, M::row_wise);
578
579 // Reset the visited flags of all vertices.
580 // As the isolated nodes will be skipped we simply mark them as visited
581
582 typedef typename G::VertexIterator Vertex;
583 Vertex vend = fineGraph.end();
584 for(Vertex vertex = fineGraph.begin(); vertex != vend; ++vertex) {
587 }
588
589 typedef typename G::MutableMatrix M;
590 SparsityBuilder<M> sparsityBuilder(*coarseMatrix);
591
592 ConnectivityConstructor<G,T>::examine(fineGraph, visitedMap, pinfo,
593 aggregates, overlap,
594 overlapVertices,
595 overlapVertices+count,
596 sparsityBuilder);
597
598 dinfo<<pinfo.communicator().rank()<<": Matrix ("<<coarseMatrix->N()<<"x"<<coarseMatrix->M()<<" row: min="<<sparsityBuilder.minRowSize()<<" max="
599 <<sparsityBuilder.maxRowSize()<<" avg="
600 <<static_cast<double>(sparsityBuilder.sumRowSize())/coarseMatrix->N()
601 <<std::endl;
602
603 delete[] overlapVertices;
604
605 return coarseMatrix;
606 }
607
608 template<class G, class V, class Set>
609 typename G::MutableMatrix*
610 GalerkinProduct<SequentialInformation>::build(G& fineGraph, V& visitedMap,
611 const SequentialInformation& pinfo,
613 const typename G::Matrix::size_type& size,
614 [[maybe_unused]] const Set& overlap)
615 {
616 typedef typename G::MutableMatrix M;
617 M* coarseMatrix = new M(size, size, M::row_wise);
618
619 // Reset the visited flags of all vertices.
620 // As the isolated nodes will be skipped we simply mark them as visited
621
622 typedef typename G::VertexIterator Vertex;
623 Vertex vend = fineGraph.end();
624 for(Vertex vertex = fineGraph.begin(); vertex != vend; ++vertex) {
627 }
628
629 SparsityBuilder<M> sparsityBuilder(*coarseMatrix);
630
631 ConnectivityConstructor<G,SequentialInformation>::examine(fineGraph, visitedMap, pinfo,
632 aggregates, sparsityBuilder);
633 dinfo<<"Matrix row: min="<<sparsityBuilder.minRowSize()<<" max="
634 <<sparsityBuilder.maxRowSize()<<" average="
635 <<static_cast<double>(sparsityBuilder.sumRowSize())/coarseMatrix->N()<<std::endl;
636 return coarseMatrix;
637 }
638
639 template<class M, class V, class P, class O>
640 void BaseGalerkinProduct::calculate(const M& fine, const AggregatesMap<V>& aggregates, M& coarse,
641 const P& pinfo, [[maybe_unused]] const O& copy)
642 {
643 coarse = static_cast<typename M::field_type>(0);
644
645 typedef typename M::ConstIterator RowIterator;
646 RowIterator endRow = fine.end();
647
648 for(RowIterator row = fine.begin(); row != endRow; ++row)
649 if(aggregates[row.index()] != AggregatesMap<V>::ISOLATED) {
650 assert(aggregates[row.index()]!=AggregatesMap<V>::UNAGGREGATED);
651 typedef typename M::ConstColIterator ColIterator;
652 ColIterator endCol = row->end();
653
654 for(ColIterator col = row->begin(); col != endCol; ++col)
655 if(aggregates[col.index()] != AggregatesMap<V>::ISOLATED) {
656 assert(aggregates[row.index()]!=AggregatesMap<V>::UNAGGREGATED);
657 coarse[aggregates[row.index()]][aggregates[col.index()]]+=*col;
658 }
659 }
660
661 // get the right diagonal matrix values on copy lines from owner processes
662 typedef typename M::block_type BlockType;
663 std::vector<BlockType> rowsize(coarse.N(),BlockType(0));
664 for (RowIterator row = coarse.begin(); row != coarse.end(); ++row)
665 rowsize[row.index()]=coarse[row.index()][row.index()];
666 pinfo.copyOwnerToAll(rowsize,rowsize);
667 for (RowIterator row = coarse.begin(); row != coarse.end(); ++row)
668 coarse[row.index()][row.index()] = rowsize[row.index()];
669
670 // don't set dirichlet boundaries for copy lines to make novlp case work,
671 // the preconditioner yields slightly different results now.
672
673 // Set the dirichlet border
674 //DirichletBoundarySetter<P>::template set<M>(coarse, pinfo, copy);
675
676 }
677
678 template<class T>
679 template<class M, class O>
680 void DirichletBoundarySetter<T>::set(M& coarse, const T& pinfo, const O& copy)
681 {
682 typedef typename T::ParallelIndexSet::const_iterator ConstIterator;
683 ConstIterator end = pinfo.indexSet().end();
684 typedef typename M::block_type Block;
685 Block identity=Block(0.0);
686 for(typename Block::RowIterator b=identity.begin(); b != identity.end(); ++b)
687 b->operator[](b.index())=1.0;
688
689 for(ConstIterator index = pinfo.indexSet().begin();
690 index != end; ++index) {
691 if(copy.contains(index->local().attribute())) {
692 typedef typename M::ColIterator ColIterator;
693 typedef typename M::row_type Row;
694 Row row = coarse[index->local()];
695 ColIterator cend = row.find(index->local());
696 ColIterator col = row.begin();
697 for(; col != cend; ++col)
698 *col = 0;
699
700 cend = row.end();
701
702 assert(col != cend); // There should be a diagonal entry
703 *col = identity;
704
705 for(++col; col != cend; ++col)
706 *col = 0;
707 }
708 }
709 }
710
711 template<class M, class O>
712 void DirichletBoundarySetter<SequentialInformation>::set(M& coarse,
713 const SequentialInformation& pinfo,
714 const O& overlap)
715 {}
716
717 } // namespace Amg
718} // namespace Dune
719#endif
Provides classes for the Coloring process of AMG.
Class providing information about the mapping of the vertices onto aggregates.
Definition: aggregates.hh:558
Visitor for identifying connected aggregates during a breadthFirstSearch.
Definition: galerkin.hh:204
Functor for building the sparsity pattern of the matrix using examineConnectivity.
Definition: galerkin.hh:61
A pair consisting of a global and local index.
Definition: indexset.hh:83
A single linked list.
Definition: sllist.hh:42
Classes for building sets out of enumeration values.
LocalIndex & local()
Get the local index.
constexpr GeometryType vertex
GeometryType representing a vertex.
Definition: type.hh:504
G Graph
The type of the graph.
Definition: galerkin.hh:209
void operator()(const ConstEdgeIterator &edge)
Process an edge pointing to another aggregate.
Definition: galerkin.hh:357
T Aggregate
The aggregate descriptor.
Definition: galerkin.hh:35
SparsityBuilder(M &matrix)
Constructor.
Definition: galerkin.hh:511
ConnectedBuilder(const AggregatesMap< Vertex > &aggregates, Graph &graph, VisitedMap &visitedMap, Set &connected)
Constructor.
Definition: galerkin.hh:350
Graph::ConstEdgeIterator ConstEdgeIterator
The constant edge iterator.
Definition: galerkin.hh:213
T Vertex
The vertex descriptor.
Definition: galerkin.hh:40
V VisitedMap
The type of the map for marking vertices as visited.
Definition: galerkin.hh:223
int visitNeighbours(const G &graph, const typename G::VertexDescriptor &vertex, V &visitor)
Visit all neighbour vertices of a vertex in a graph.
static const Vertex ISOLATED
Identifier of isolated vertices.
Definition: aggregates.hh:569
S Set
The type of the connected set.
Definition: galerkin.hh:218
Aggregate * aggregate
The aggregate the vertex belongs to.
Definition: galerkin.hh:45
Vertex vertex
The vertex descriptor.
Definition: galerkin.hh:50
Graph::VertexDescriptor Vertex
The vertex descriptor of the graph.
Definition: galerkin.hh:228
void calculate(const M &fine, const AggregatesMap< V > &aggregates, M &coarse, const I &pinfo, const O &copy)
Calculate the galerkin product.
auto min(ADLTag< 0 >, const V &v1, const V &v2)
implements binary Simd::min()
Definition: defaults.hh:87
auto max(ADLTag< 0 >, const V &v1, const V &v2)
implements binary Simd::max()
Definition: defaults.hh:79
DVVerbType dvverb(std::cout)
stream for very verbose output.
Definition: stdstreams.hh:93
DInfoType dinfo(std::cout)
Stream for informative output.
Definition: stdstreams.hh:138
DVerbType dverb(std::cout)
Singleton of verbose debug stream.
Definition: stdstreams.hh:114
Dune namespace.
Definition: alignedallocator.hh:11
STL namespace.
An stl-compliant pool allocator.
@ nonoverlapping
Category for non-overlapping solvers.
Definition: solvercategory.hh:25
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
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 13, 23:29, 2024)