4 #ifndef DUNE_GALERKIN_HH
5 #define DUNE_GALERKIN_HH
9 #include <dune/common/poolallocator.hh>
10 #include <dune/common/enumset.hh>
11 #include <dune/common/unused.hh>
87 typename M::CreateIterator row_;
89 std::size_t minRowSize_;
91 std::size_t maxRowSize_;
92 std::size_t sumRowSize_;
93 #ifdef DUNE_ISTL_WITH_CHECKING
94 bool diagonalInserted;
109 template<
class M,
class V,
class I,
class O>
111 const I& pinfo,
const O& copy);
131 template<
class G,
class V,
class Set>
132 typename G::MutableMatrix*
build(G& fineGraph, V& visitedMap,
135 const typename G::Matrix::size_type& size,
145 template<
class G,
class I,
class Set>
147 buildOverlapVertices(
const G& graph,
const I& pinfo,
150 std::size_t& overlapCount);
176 template<
class G,
class V,
class Set>
177 typename G::MutableMatrix*
build(G& fineGraph, V& visitedMap,
180 const typename G::Matrix::size_type& size,
186 template<
class R,
class G,
class V>
195 template<
class R,
class G,
class V>
198 const typename G::VertexDescriptor& seed);
204 template<
class G,
class S,
class V>
230 typedef typename Graph::VertexDescriptor
Vertex;
269 template<
class G,
class T>
272 typedef typename G::VertexDescriptor
Vertex;
274 template<
class V,
class O,
class R>
288 typedef typename G::VertexDescriptor
Vertex;
290 template<
class V,
class R>
301 template<
class M,
class O>
302 static void set(M& coarse,
const T& pinfo,
const O& copy);
308 template<
class M,
class O>
312 template<
class R,
class G,
class V>
315 const typename G::VertexDescriptor& seed)
317 assert(row.index()==aggregates[seed]);
318 row.insert(aggregates[seed]);
320 typedef typename G::VertexDescriptor Vertex;
321 typedef std::allocator<Vertex> Allocator;
322 typedef SLList<Vertex,Allocator> VertexList;
326 aggregates.template breadthFirstSearch<true,false>(seed,aggregates[seed], graph, vlist, dummy,
327 conBuilder, visitedMap);
330 template<
class R,
class G,
class V>
337 const typename G::VertexDescriptor aggregate=*seed->
aggregate;
340 while(seed != overlapEnd && aggregate == *seed->
aggregate) {
345 put(visitedMap, seed->
vertex,
true);
351 template<
class G,
class S,
class V>
355 : aggregates_(aggregates), graph_(graph), visitedMap_(visitedMap), connected_(connected)
358 template<
class G,
class S,
class V>
361 typedef typename G::VertexDescriptor
Vertex;
362 const Vertex& vertex = aggregates_[edge.target()];
365 connected_.insert(vertex);
369 template<
class G,
class I,
class Set>
374 std::size_t& overlapCount)
377 typedef typename G::ConstVertexIterator ConstIterator;
378 typedef typename I::GlobalLookupIndexSet GlobalLookup;
379 typedef typename GlobalLookup::IndexPair IndexPair;
381 const ConstIterator end = graph.end();
384 const GlobalLookup& lookup=pinfo.globalLookup();
386 for(ConstIterator vertex=graph.begin(); vertex != end; ++vertex) {
387 const IndexPair* pair = lookup.pair(*vertex);
389 if(pair!=0 && overlap.contains(pair->local().attribute()))
393 typedef typename G::VertexDescriptor Vertex;
397 return overlapVertices;
401 for(ConstIterator vertex=graph.begin(); vertex != end; ++vertex) {
402 const IndexPair* pair = lookup.pair(*vertex);
404 if(pair!=0 && overlap.contains(pair->local().attribute())) {
405 overlapVertices[overlapCount].
aggregate = &aggregates[pair->local()];
406 overlapVertices[overlapCount].
vertex = pair->local();
411 dverb << overlapCount<<
" overlap vertices"<<std::endl;
413 std::sort(overlapVertices, overlapVertices+overlapCount, OVLess<Vertex>());
416 return overlapVertices;
419 template<
class G,
class T>
420 template<
class V,
class O,
class R>
430 typedef typename T::GlobalLookupIndexSet GlobalLookup;
431 const GlobalLookup& lookup = pinfo.globalLookup();
433 typedef typename G::VertexIterator VertexIterator;
435 VertexIterator vend=graph.end();
437 #ifdef DUNE_ISTL_WITH_CHECKING
438 std::set<Vertex> examined;
444 for(VertexIterator vertex = graph.begin(); vertex != vend; ++vertex)
445 if(!
get(visitedMap, *vertex)) {
447 typedef typename GlobalLookup::IndexPair IndexPair;
448 const IndexPair* pair = lookup.pair(*vertex);
449 if(pair==0 || !overlap.contains(pair->local().attribute())) {
450 #ifdef DUNE_ISTL_WITH_CHECKING
451 assert(examined.find(aggregates[*vertex])==examined.end());
452 examined.insert(aggregates[*vertex]);
459 if(overlapVertices != overlapEnd) {
472 dvverb<<
"constructed "<<row.index()<<
" non-overlapping rows"<<std::endl;
476 while(overlapVertices != overlapEnd)
479 #ifdef DUNE_ISTL_WITH_CHECKING
480 typedef typename GlobalLookup::IndexPair IndexPair;
481 const IndexPair* pair = lookup.pair(overlapVertices->
vertex);
482 assert(pair!=0 && overlap.contains(pair->local().attribute()));
483 assert(examined.find(aggregates[overlapVertices->
vertex])==examined.end());
484 examined.insert(aggregates[overlapVertices->
vertex]);
494 template<
class V,
class R>
501 DUNE_UNUSED_PARAMETER(pinfo);
502 typedef typename G::VertexIterator VertexIterator;
504 VertexIterator vend=graph.end();
505 for(VertexIterator vertex = graph.begin(); vertex != vend; ++vertex) {
506 if(!
get(visitedMap, *vertex)) {
516 : row_(matrix.createbegin()),
517 minRowSize_(std::numeric_limits<std::size_t>::max()),
518 maxRowSize_(0), sumRowSize_(0)
520 #ifdef DUNE_ISTL_WITH_CHECKING
521 diagonalInserted =
false;
543 sumRowSize_ += row_.size();
544 minRowSize_=std::min(minRowSize_, row_.size());
545 maxRowSize_=std::max(maxRowSize_, row_.size());
547 #ifdef DUNE_ISTL_WITH_CHECKING
548 assert(diagonalInserted);
549 diagonalInserted =
false;
557 #ifdef DUNE_ISTL_WITH_CHECKING
558 diagonalInserted = diagonalInserted || row_.index()==
index;
563 template<
class G,
class V,
class Set>
564 typename G::MutableMatrix*
568 const typename G::Matrix::size_type& size,
575 const OverlapVertex* overlapVertices = buildOverlapVertices(fineGraph,
580 typedef typename G::MutableMatrix M;
581 M* coarseMatrix =
new M(size, size, M::row_wise);
586 typedef typename G::VertexIterator Vertex;
587 Vertex vend = fineGraph.end();
588 for(Vertex vertex = fineGraph.begin(); vertex != vend; ++vertex) {
593 typedef typename G::MutableMatrix M;
599 overlapVertices+count,
602 dinfo<<pinfo.communicator().rank()<<
": Matrix ("<<coarseMatrix->N()<<
"x"<<coarseMatrix->M()<<
" row: min="<<sparsityBuilder.
minRowSize()<<
" max="
604 <<
static_cast<double>(sparsityBuilder.
sumRowSize())/coarseMatrix->N()
607 delete[] overlapVertices;
612 template<
class G,
class V,
class Set>
613 typename G::MutableMatrix*
617 const typename G::Matrix::size_type& size,
620 DUNE_UNUSED_PARAMETER(overlap);
621 typedef typename G::MutableMatrix M;
622 M* coarseMatrix =
new M(size, size, M::row_wise);
627 typedef typename G::VertexIterator Vertex;
628 Vertex vend = fineGraph.end();
629 for(Vertex vertex = fineGraph.begin(); vertex != vend; ++vertex) {
637 aggregates, sparsityBuilder);
638 dinfo<<
"Matrix row: min="<<sparsityBuilder.
minRowSize()<<
" max="
640 <<
static_cast<double>(sparsityBuilder.
sumRowSize())/coarseMatrix->N()<<std::endl;
644 template<
class M,
class V,
class P,
class O>
646 const P& pinfo,
const O& copy)
648 DUNE_UNUSED_PARAMETER(copy);
649 coarse =
static_cast<typename M::field_type
>(0);
651 typedef typename M::ConstIterator RowIterator;
652 RowIterator endRow = fine.end();
654 for(RowIterator
row = fine.begin();
row != endRow; ++
row)
657 typedef typename M::ConstColIterator ColIterator;
658 ColIterator endCol =
row->end();
660 for(ColIterator
col =
row->begin();
col != endCol; ++
col)
663 coarse[aggregates[
row.index()]][aggregates[
col.index()]]+=*
col;
668 typedef typename M::block_type BlockType;
669 std::vector<BlockType> rowsize(coarse.N(),BlockType(0));
670 for (RowIterator
row = coarse.begin();
row != coarse.end(); ++
row)
671 rowsize[
row.index()]=coarse[
row.index()][
row.index()];
672 pinfo.copyOwnerToAll(rowsize,rowsize);
673 for (RowIterator
row = coarse.begin();
row != coarse.end(); ++
row)
674 coarse[
row.index()][
row.index()] = rowsize[
row.index()];
685 template<
class M,
class O>
688 typedef typename T::ParallelIndexSet::const_iterator ConstIterator;
689 ConstIterator end = pinfo.indexSet().end();
690 typedef typename M::block_type Block;
691 Block identity=Block(0.0);
692 for(
typename Block::RowIterator b=identity.begin(); b != identity.end(); ++b)
693 b->operator[](b.index())=1.0;
695 for(ConstIterator
index = pinfo.indexSet().begin();
697 if(copy.contains(
index->local().attribute())) {
698 typedef typename M::ColIterator ColIterator;
699 typedef typename M::row_type Row;
701 ColIterator cend = row.find(
index->local());
702 ColIterator
col = row.begin();
703 for(; col != cend; ++
col)
711 for(++col; col != cend; ++
col)
717 template<
class M,
class O>
S Set
The type of the connected set.
Definition: galerkin.hh:220
SparsityBuilder(M &matrix)
Constructor.
Definition: galerkin.hh:515
int visitNeighbours(const G &graph, const typename G::VertexDescriptor &vertex, V &visitor)
Visit all neighbour vertices of a vertex in a graph.
void calculate(const M &fine, const AggregatesMap< V > &aggregates, M &coarse, const I &pinfo, const O ©)
Calculate the galerkin product.
std::size_t minRowSize()
Definition: galerkin.hh:530
static void set(M &coarse, const T &pinfo, const O ©)
Definition: galerkin.hh:686
ConnectedBuilder(const AggregatesMap< Vertex > &aggregates, Graph &graph, VisitedMap &visitedMap, Set &connected)
Constructor.
Definition: galerkin.hh:352
G::VertexDescriptor Vertex
Definition: galerkin.hh:272
Aggregate * aggregate
The aggregate the vertex belongs to.
Definition: galerkin.hh:47
Definition: galerkin.hh:270
V VisitedMap
The type of the map for marking vertices as visited.
Definition: galerkin.hh:225
void insert(const typename M::size_type &index)
Definition: galerkin.hh:554
void operator++()
Definition: galerkin.hh:541
Definition: galerkin.hh:116
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:565
Definition: galerkin.hh:98
Visitor for identifying connected aggregates during a breadthFirstSearch.
Definition: galerkin.hh:205
Col col
Definition: matrixmatrix.hh:347
Class providing information about the mapping of the vertices onto aggregates.
Definition: aggregates.hh:497
Definition: galerkin.hh:32
T Aggregate
The aggregate descriptor.
Definition: galerkin.hh:37
Functor for building the sparsity pattern of the matrix using examineConnectivity.
Definition: galerkin.hh:62
Provides classes for the Coloring process of AMG.
T ParallelInformation
Definition: galerkin.hh:120
Category for on overlapping solvers.
Definition: solvercategory.hh:24
G::VertexDescriptor Vertex
Definition: galerkin.hh:288
Graph::VertexDescriptor Vertex
The vertex descriptor of the graph.
Definition: galerkin.hh:230
Graph::ConstEdgeIterator ConstEdgeIterator
The constant edge iterator.
Definition: galerkin.hh:215
std::size_t index()
Definition: galerkin.hh:81
Vertex vertex
The vertex descriptor.
Definition: galerkin.hh:52
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:313
Row row
Definition: matrixmatrix.hh:345
Definition: galerkin.hh:184
std::size_t maxRowSize()
Definition: galerkin.hh:525
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:331
std::size_t index
Definition: matrixmarket.hh:561
std::size_t count
Definition: matrixmatrix.hh:215
Definition: galerkin.hh:299
std::size_t sumRowSize()
Definition: galerkin.hh:536
G Graph
The type of the graph.
Definition: galerkin.hh:211
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:421
bool operator()(const OverlapVertex< A > &o1, const OverlapVertex< A > &o2)
Definition: galerkin.hh:155
T Vertex
The vertex descriptor.
Definition: galerkin.hh:42
void operator()(const ConstEdgeIterator &edge)
Process an edge pointing to another aggregate.
Definition: galerkin.hh:359