3#ifndef DUNE_GALERKIN_HH
4#define DUNE_GALERKIN_HH
70 void insert(
const typename M::size_type& index);
74 std::size_t minRowSize();
76 std::size_t maxRowSize();
78 std::size_t sumRowSize();
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;
96 class BaseGalerkinProduct
107 template<
class M,
class V,
class I,
class O>
109 const I& pinfo,
const O& copy);
114 class GalerkinProduct
115 :
public BaseGalerkinProduct
118 typedef T ParallelInformation;
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,
143 template<
class G,
class I,
class Set>
144 const OverlapVertex<typename G::VertexDescriptor>*
145 buildOverlapVertices(
const G& graph,
const I& pinfo,
148 std::size_t& overlapCount);
153 bool operator()(
const OverlapVertex<A>& o1,
const OverlapVertex<A>& o2)
155 return *o1.aggregate < *o2.aggregate;
161 class GalerkinProduct<SequentialInformation>
162 :
public BaseGalerkinProduct
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,
182 struct BaseConnectivityConstructor
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);
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);
202 template<
class G,
class S,
class V>
228 typedef typename Graph::VertexDescriptor
Vertex;
267 template<
class G,
class T>
268 struct ConnectivityConstructor :
public BaseConnectivityConstructor
270 typedef typename G::VertexDescriptor Vertex;
272 template<
class V,
class O,
class R>
273 static void examine(G& graph,
278 const OverlapVertex<Vertex>* overlapVertices,
279 const OverlapVertex<Vertex>* overlapEnd,
284 struct ConnectivityConstructor<G,SequentialInformation> :
public BaseConnectivityConstructor
286 typedef typename G::VertexDescriptor Vertex;
288 template<
class V,
class R>
289 static void examine(G& graph,
291 const SequentialInformation& pinfo,
292 const AggregatesMap<Vertex>& aggregates,
297 struct DirichletBoundarySetter
299 template<
class M,
class O>
300 static void set(M& coarse,
const T& pinfo,
const O& copy);
304 struct DirichletBoundarySetter<SequentialInformation>
306 template<
class M,
class O>
307 static void set(M& coarse,
const SequentialInformation& pinfo,
const O& copy);
310 template<
class R,
class G,
class V>
311 void BaseConnectivityConstructor::constructNonOverlapConnectivity(R& row, G& graph, V& visitedMap,
313 const typename G::VertexDescriptor& seed)
315 assert(row.index()==aggregates[seed]);
316 row.insert(aggregates[seed]);
318 typedef typename G::VertexDescriptor Vertex;
319 typedef std::allocator<Vertex> Allocator;
324 aggregates.template breadthFirstSearch<true,false>(seed,aggregates[seed], graph, vlist, dummy,
325 conBuilder, visitedMap);
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)
334 ConnectedBuilder<G,R,V> conBuilder(aggregates, graph, visitedMap, row);
335 const typename G::VertexDescriptor aggregate=*seed->aggregate;
337 if (row.index()==*seed->aggregate) {
338 while(seed != overlapEnd && aggregate == *seed->aggregate) {
339 row.insert(*seed->aggregate);
343 put(visitedMap, seed->vertex,
true);
349 template<
class G,
class S,
class V>
353 : aggregates_(aggregates), graph_(graph), visitedMap_(visitedMap), connected_(connected)
356 template<
class G,
class S,
class V>
362 connected_.insert(
vertex);
366 template<
class G,
class I,
class Set>
367 const OverlapVertex<typename G::VertexDescriptor>*
368 GalerkinProduct<T>::buildOverlapVertices(
const G& graph,
const I& pinfo,
371 std::size_t& overlapCount)
374 typedef typename G::ConstVertexIterator ConstIterator;
375 typedef typename I::GlobalLookupIndexSet GlobalLookup;
376 typedef typename GlobalLookup::IndexPair
IndexPair;
378 const ConstIterator end = graph.end();
381 const GlobalLookup& lookup=pinfo.globalLookup();
386 if(pair!=0 && overlap.contains(pair->
local().attribute()))
390 typedef typename G::VertexDescriptor Vertex;
392 OverlapVertex<Vertex>* overlapVertices =
new OverlapVertex<Vertex>[overlapCount=0 ? 1 : overlapCount];
394 return overlapVertices;
401 if(pair!=0 && overlap.contains(pair->
local().attribute())) {
402 overlapVertices[overlapCount].aggregate = &aggregates[pair->
local()];
403 overlapVertices[overlapCount].vertex = pair->
local();
408 dverb << overlapCount<<
" overlap vertices"<<std::endl;
410 std::sort(overlapVertices, overlapVertices+overlapCount, OVLess<Vertex>());
413 return overlapVertices;
416 template<
class G,
class T>
417 template<
class V,
class O,
class R>
418 void ConnectivityConstructor<G,T>::examine(G& graph,
421 const AggregatesMap<Vertex>& aggregates,
423 const OverlapVertex<Vertex>* overlapVertices,
424 const OverlapVertex<Vertex>* overlapEnd,
427 typedef typename T::GlobalLookupIndexSet GlobalLookup;
428 const GlobalLookup& lookup = pinfo.globalLookup();
430 typedef typename G::VertexIterator VertexIterator;
432 VertexIterator vend=graph.end();
434#ifdef DUNE_ISTL_WITH_CHECKING
435 std::set<Vertex> examined;
442 if(!get(visitedMap, *
vertex)) {
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]);
451 constructNonOverlapConnectivity(row, graph, visitedMap, aggregates, *
vertex);
456 if(overlapVertices != overlapEnd) {
458 constructOverlapConnectivity(row, graph, visitedMap, aggregates, overlapVertices, overlapEnd);
469 dvverb<<
"constructed "<<row.index()<<
" non-overlapping rows"<<std::endl;
473 while(overlapVertices != overlapEnd)
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]);
483 constructOverlapConnectivity(row, graph, visitedMap, aggregates, overlapVertices, overlapEnd);
491 template<
class V,
class R>
492 void ConnectivityConstructor<G,SequentialInformation>::examine(G& graph,
494 [[maybe_unused]]
const SequentialInformation& pinfo,
495 const AggregatesMap<Vertex>& aggregates,
498 typedef typename G::VertexIterator VertexIterator;
500 VertexIterator vend=graph.end();
502 if(!get(visitedMap, *
vertex)) {
503 constructNonOverlapConnectivity(row, graph, visitedMap, aggregates, *
vertex);
512 : row_(matrix.createbegin()),
513 minRowSize_(
std::numeric_limits<
std::size_t>::
max()),
514 maxRowSize_(0), sumRowSize_(0)
516#ifdef DUNE_ISTL_WITH_CHECKING
517 diagonalInserted =
false;
526 std::size_t SparsityBuilder<M>::minRowSize()
532 std::size_t SparsityBuilder<M>::sumRowSize()
537 void SparsityBuilder<M>::operator++()
539 sumRowSize_ += row_.size();
540 minRowSize_=
std::min(minRowSize_, row_.size());
541 maxRowSize_=
std::max(maxRowSize_, row_.size());
543#ifdef DUNE_ISTL_WITH_CHECKING
544 assert(diagonalInserted);
545 diagonalInserted =
false;
550 void SparsityBuilder<M>::insert(
const typename M::size_type& index)
553#ifdef DUNE_ISTL_WITH_CHECKING
554 diagonalInserted = diagonalInserted || row_.index()==index;
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,
567 typedef OverlapVertex<typename G::VertexDescriptor> OverlapVertex;
571 const OverlapVertex* overlapVertices = buildOverlapVertices(fineGraph,
576 typedef typename G::MutableMatrix M;
577 M* coarseMatrix =
new M(size, size, M::row_wise);
582 typedef typename G::VertexIterator Vertex;
583 Vertex vend = fineGraph.end();
589 typedef typename G::MutableMatrix M;
592 ConnectivityConstructor<G,T>::examine(fineGraph, visitedMap, pinfo,
595 overlapVertices+count,
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()
603 delete[] overlapVertices;
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)
616 typedef typename G::MutableMatrix M;
617 M* coarseMatrix =
new M(size, size, M::row_wise);
622 typedef typename G::VertexIterator Vertex;
623 Vertex vend = fineGraph.end();
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;
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)
643 coarse =
static_cast<typename M::field_type
>(0);
645 typedef typename M::ConstIterator RowIterator;
646 RowIterator endRow = fine.end();
648 for(RowIterator row = fine.begin(); row != endRow; ++row)
651 typedef typename M::ConstColIterator ColIterator;
652 ColIterator endCol = row->end();
654 for(ColIterator col = row->begin(); col != endCol; ++col)
657 coarse[aggregates[row.index()]][aggregates[col.index()]]+=*col;
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()];
679 template<
class M,
class O>
680 void DirichletBoundarySetter<T>::set(M& coarse,
const T& pinfo,
const O& copy)
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;
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)
705 for(++col; col != cend; ++col)
711 template<
class M,
class O>
712 void DirichletBoundarySetter<SequentialInformation>::set(M& coarse,
713 const SequentialInformation& pinfo,
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 ©)
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
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