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