3 #ifndef DUNE_GALERKIN_HH
4 #define DUNE_GALERKIN_HH
8 #include <dune/common/poolallocator.hh>
9 #include <dune/common/enumset.hh>
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;
107 template<
class M,
class V,
class I,
class O>
109 const I& pinfo,
const O& copy);
129 template<
class G,
class V,
class Set>
130 typename G::MutableMatrix*
build(G& fineGraph, V& visitedMap,
133 const typename G::Matrix::size_type& size,
143 template<
class G,
class I,
class Set>
145 buildOverlapVertices(
const G& graph,
const I& pinfo,
148 std::size_t& overlapCount);
174 template<
class G,
class V,
class Set>
175 typename G::MutableMatrix*
build(G& fineGraph, V& visitedMap,
178 const typename G::Matrix::size_type& size,
184 template<
class R,
class G,
class V>
193 template<
class R,
class G,
class V>
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>
270 typedef typename G::VertexDescriptor
Vertex;
272 template<
class V,
class O,
class R>
286 typedef typename G::VertexDescriptor
Vertex;
288 template<
class V,
class R>
299 template<
class M,
class O>
300 static void set(M& coarse,
const T& pinfo,
const O& copy);
306 template<
class M,
class O>
310 template<
class R,
class G,
class V>
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;
320 typedef SLList<Vertex,Allocator> VertexList;
324 aggregates.template breadthFirstSearch<true,false>(seed,aggregates[seed], graph, vlist, dummy,
325 conBuilder, visitedMap);
328 template<
class R,
class G,
class V>
335 const typename G::VertexDescriptor aggregate=*seed->
aggregate;
338 while(seed != overlapEnd && aggregate == *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>
359 const Vertex& vertex = aggregates_[edge.target()];
362 connected_.insert(vertex);
366 template<
class G,
class I,
class Set>
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();
383 for(ConstIterator vertex=graph.begin(); vertex != end; ++vertex) {
384 const IndexPair* pair = lookup.pair(*vertex);
386 if(pair!=0 && overlap.contains(pair->local().attribute()))
390 typedef typename G::VertexDescriptor Vertex;
394 return overlapVertices;
398 for(ConstIterator vertex=graph.begin(); vertex != end; ++vertex) {
399 const IndexPair* pair = lookup.pair(*vertex);
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>
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;
441 for(VertexIterator vertex = graph.begin(); vertex != vend; ++vertex)
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]);
456 if(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]);
491 template<
class V,
class R>
498 typedef typename G::VertexIterator VertexIterator;
500 VertexIterator vend=graph.end();
501 for(VertexIterator vertex = graph.begin(); vertex != vend; ++vertex) {
502 if(!
get(visitedMap, *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;
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;
553 #ifdef DUNE_ISTL_WITH_CHECKING
554 diagonalInserted = diagonalInserted || row_.index()==index;
559 template<
class G,
class V,
class Set>
560 typename G::MutableMatrix*
564 const typename G::Matrix::size_type& size,
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();
584 for(Vertex vertex = fineGraph.begin(); vertex != vend; ++vertex) {
589 typedef typename G::MutableMatrix M;
595 overlapVertices+count,
598 dinfo<<pinfo.communicator().rank()<<
": Matrix ("<<coarseMatrix->N()<<
"x"<<coarseMatrix->M()<<
" row: min="<<sparsityBuilder.
minRowSize()<<
" max="
600 <<
static_cast<double>(sparsityBuilder.
sumRowSize())/coarseMatrix->N()
603 delete[] overlapVertices;
608 template<
class G,
class V,
class Set>
609 typename G::MutableMatrix*
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();
624 for(Vertex vertex = fineGraph.begin(); vertex != vend; ++vertex) {
632 aggregates, sparsityBuilder);
633 dinfo<<
"Matrix row: min="<<sparsityBuilder.
minRowSize()<<
" max="
635 <<
static_cast<double>(sparsityBuilder.
sumRowSize())/coarseMatrix->N()<<std::endl;
639 template<
class M,
class V,
class P,
class O>
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>
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();
711 template<
class M,
class O>
Provides classes for the Coloring process of AMG.
Col col
Definition: matrixmatrix.hh:349
bool operator()(const OverlapVertex< A > &o1, const OverlapVertex< A > &o2)
Definition: galerkin.hh:153
static void constructNonOverlapConnectivity(R &row, G &graph, V &visitedMap, const AggregatesMap< typename G::VertexDescriptor > &aggregates, const typename G::VertexDescriptor &seed)
Construct the connectivity of an aggregate in the overlap.
Definition: galerkin.hh:311
G Graph
The type of the graph.
Definition: galerkin.hh:209
void operator++()
Definition: galerkin.hh:537
void operator()(const ConstEdgeIterator &edge)
Process an edge pointing to another aggregate.
Definition: galerkin.hh:357
void insert(const typename M::size_type &index)
Definition: galerkin.hh:550
static void constructOverlapConnectivity(R &row, G &graph, V &visitedMap, const AggregatesMap< typename G::VertexDescriptor > &aggregates, const OverlapVertex< typename G::VertexDescriptor > *&seed, const OverlapVertex< typename G::VertexDescriptor > *overlapEnd)
Definition: galerkin.hh:329
G::MutableMatrix * build(G &fineGraph, V &visitedMap, const ParallelInformation &pinfo, AggregatesMap< typename G::VertexDescriptor > &aggregates, const typename G::Matrix::size_type &size, const Set ©)
Calculates the coarse matrix via a Galerkin product.
Definition: galerkin.hh:561
std::size_t index()
Definition: galerkin.hh:79
T ParallelInformation
Definition: galerkin.hh:118
G::VertexDescriptor Vertex
Definition: galerkin.hh:270
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
G::VertexDescriptor Vertex
Definition: galerkin.hh:286
static void examine(G &graph, V &visitedMap, const T &pinfo, const AggregatesMap< Vertex > &aggregates, const O &overlap, const OverlapVertex< Vertex > *overlapVertices, const OverlapVertex< Vertex > *overlapEnd, R &row)
Definition: galerkin.hh:418
G::MutableMatrix * build(G &fineGraph, V &visitedMap, const SequentialInformation &pinfo, const AggregatesMap< typename G::VertexDescriptor > &aggregates, const typename G::Matrix::size_type &size, const Set ©)
Calculates the coarse matrix via a Galerkin product.
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.
S Set
The type of the connected set.
Definition: galerkin.hh:218
Aggregate * aggregate
The aggregate the vertex belongs to.
Definition: galerkin.hh:45
std::size_t sumRowSize()
Definition: galerkin.hh:532
Vertex vertex
The vertex descriptor.
Definition: galerkin.hh:50
std::size_t minRowSize()
Definition: galerkin.hh:526
static void set(M &coarse, const T &pinfo, const O ©)
Definition: galerkin.hh:680
Graph::VertexDescriptor Vertex
The vertex descriptor of the graph.
Definition: galerkin.hh:228
static void examine(G &graph, V &visitedMap, const SequentialInformation &pinfo, const AggregatesMap< Vertex > &aggregates, R &row)
void calculate(const M &fine, const AggregatesMap< V > &aggregates, M &coarse, const I &pinfo, const O ©)
Calculate the galerkin product.
std::size_t maxRowSize()
Definition: galerkin.hh:521
Definition: allocator.hh:9
PropertyMapTypeSelector< Amg::VertexVisitedTag, Amg::PropertiesGraph< G, Amg::VertexProperties, EP, VM, EM > >::Type get([[maybe_unused]] const Amg::VertexVisitedTag &tag, Amg::PropertiesGraph< G, Amg::VertexProperties, EP, VM, EM > &graph)
Definition: dependency.hh:291
Class providing information about the mapping of the vertices onto aggregates.
Definition: aggregates.hh:558
Definition: galerkin.hh:31
Functor for building the sparsity pattern of the matrix using examineConnectivity.
Definition: galerkin.hh:61
Definition: galerkin.hh:97
Definition: galerkin.hh:116
Definition: galerkin.hh:183
Visitor for identifying connected aggregates during a breadthFirstSearch.
Definition: galerkin.hh:204
Definition: galerkin.hh:269
Definition: galerkin.hh:298
@ 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