3#ifndef DUNE_AMG_AGGREGATES_HH
4#define DUNE_AMG_AGGREGATES_HH
10#include "combinedfunctor.hh"
82 this->setMaxDistance(diameter-1);
87 this->setMaxDistance(this->maxDistance()+diameter-1);
89 this->setMinAggregateSize(csize);
90 this->setMaxAggregateSize(
static_cast<std::size_t
>(csize*1.5));
106 this->setMaxDistance(this->maxDistance()+dim-1);
111 std::ostream& operator<<(std::ostream& os,
const AggregationCriterion<T>& criterion)
113 os<<
"{ maxdistance="<<criterion.maxDistance()<<
" minAggregateSize="
114 <<criterion.minAggregateSize()<<
" maxAggregateSize="<<criterion.maxAggregateSize()
115 <<
" connectivity="<<criterion.maxConnectivity()<<
" debugLevel="<<criterion.debugLevel()<<
"}";
130 template<
class M,
class N>
154 void init(
const Matrix* matrix);
156 void initRow(
const Row& row,
int index);
158 void examine(
const ColIter& col);
161 void examine(G& graph,
const typename G::EdgeIterator& edge,
const ColIter& col);
178 typedef typename FieldTraits<field_type>::real_type real_type;
186 std::vector<real_type> vals_;
187 typename std::vector<real_type>::iterator valIter_;
192 template<
class M,
class N>
198 template<
class M,
class N>
199 inline void SymmetricMatrixDependency<M,N>::initRow(
const Row& row,
int index)
202 vals_.assign(row.size(), 0.0);
203 assert(vals_.size()==row.size());
204 valIter_=vals_.begin();
207 diagonal_=norm_(row[index]);
211 template<
class M,
class N>
212 inline void SymmetricMatrixDependency<M,N>::examine(
const ColIter& col)
216 real_type eij = norm_(*col);
217 if(!N::is_sign_preserving || eij<0)
219 *valIter_ = eij/diagonal_*eij/norm_(matrix_->operator[](col.index())[col.index()]);
220 maxValue_ =
max(maxValue_, *valIter_);
226 template<
class M,
class N>
228 inline void SymmetricMatrixDependency<M,N>::examine(G&,
const typename G::EdgeIterator& edge,
const ColIter&)
230 if(*valIter_ > alpha() * maxValue_) {
231 edge.properties().setDepends();
232 edge.properties().setInfluences();
237 template<
class M,
class N>
238 inline bool SymmetricMatrixDependency<M,N>::isIsolated()
242 valIter_=vals_.begin();
243 return maxValue_ < beta();
249 template<
class M,
class N>
273 void init(
const Matrix* matrix);
275 void initRow(
const Row& row,
int index);
277 void examine(
const ColIter& col);
280 void examine(G& graph,
const typename G::EdgeIterator& edge,
const ColIter& col);
297 typedef typename FieldTraits<field_type>::real_type real_type;
310 template<
class M,
class N>
334 void init(
const Matrix* matrix);
336 void initRow(
const Row& row,
int index);
338 void examine(
const ColIter& col);
341 void examine(G& graph,
const typename G::EdgeIterator& edge,
const ColIter& col);
358 typedef typename FieldTraits<field_type>::real_type real_type;
367 void initRow(
const Row& row,
int index,
const std::true_type&);
368 void initRow(
const Row& row,
int index,
const std::false_type&);
380 is_sign_preserving =
true
388 typename FieldTraits<typename M::field_type>::real_type
operator()(
const M& m,
391 typedef typename M::field_type field_type;
392 typedef typename FieldTraits<field_type>::real_type real_type;
393 static_assert( std::is_convertible<field_type, real_type >::value,
394 "use of diagonal norm in AMG not implemented for complex field_type");
407 typedef typename FieldTraits<M>::real_type real_type;
408 static_assert( std::is_convertible<M, real_type >::value,
409 "use of diagonal norm in AMG not implemented for complex field_type");
418 static T signed_abs(
const T & v)
425 static T signed_abs(
const std::complex<T> & v)
429 return csgn(v) * std::abs(v);
434 static T csgn(
const T & v)
436 return (T(0) < v) - (v < T(0));
441 static T csgn(std::complex<T> a)
443 return csgn(a.real())+(a.real() == 0.0)*csgn(a.imag());
464 is_sign_preserving =
false
471 typename FieldTraits<typename M::field_type>::real_type
operator()(
const M& m)
const
473 return m.infinity_norm();
481 is_sign_preserving =
false
488 typename FieldTraits<typename M::field_type>::real_type
operator()(
const M& m)
const
490 return m.frobenius_norm();
497 is_sign_preserving =
false
504 typename FieldTraits<typename M::field_type>::real_type
operator()(
const M& m)
const
515 template<
class M,
class Norm>
535 template<
class M,
class Norm>
598 template<
class EdgeIterator>
599 void operator()([[maybe_unused]]
const EdgeIterator& edge)
const
632 template<
class M,
class G,
class C>
633 std::tuple<int,int,int,int>
buildAggregates(
const M& matrix, G& graph,
const C& criterion,
653 template<
bool reset,
class G,
class F,
class VM>
658 VM& visitedMap)
const;
683 template<
bool remove,
bool reset,
class G,
class L,
class F1,
class F2,
class VM>
686 const G& graph, L& visited, F1& aggregateVisitor,
687 F2& nonAggregateVisitor,
688 VM& visitedMap)
const;
723 const_iterator begin()
const
728 const_iterator end()
const
748 AggregatesMap<V>& operator=(
const AggregatesMap<V>&) =
delete;
758 std::size_t noVertices_;
764 template<
class G,
class C>
766 const typename C::Matrix& matrix,
774 template<
class G,
class S>
783 typedef G MatrixGraph;
818 VertexSet& connectivity, std::vector<Vertex>& front_);
852 typename VertexSet::size_type
size();
899 std::vector<Vertex>& front_;
949 template<
class M,
class C>
950 std::tuple<int,int,int,int>
build(
const M& m, G& graph,
968 typedef std::set<Vertex,std::less<Vertex>,
Allocator> VertexSet;
973 typedef std::size_t* SphereMap;
988 std::vector<Vertex> front_;
993 VertexSet connected_;
1006 static const Vertex NullEntry;
1014 enum { N = 1300000 };
1048 const AggregatesMap<Vertex>& aggregates,
1056 class AggregateVisitor
1116 class FrontNeighbourCounter :
public Counter
1140 class TwoWayCounter :
public Counter
1163 class OneWayCounter :
public Counter
1181 const AggregatesMap<Vertex>& aggregates)
const;
1189 class ConnectivityCounter :
public Counter
1204 const VertexSet& connected_;
1249 class DependencyCounter :
public Counter
1281 std::vector<Vertex>& front_;
1396 template<
class M,
class N>
1402 template<
class M,
class N>
1403 inline void SymmetricDependency<M,N>::initRow(
const Row& row,
int index)
1405 initRow(row, index, std::is_convertible<field_type, real_type>());
1408 template<
class M,
class N>
1409 inline void SymmetricDependency<M,N>::initRow(
const Row& row,
int index,
const std::false_type&)
1411 DUNE_THROW(InvalidStateException,
"field_type needs to convertible to real_type");
1414 template<
class M,
class N>
1415 inline void SymmetricDependency<M,N>::initRow([[maybe_unused]]
const Row& row,
int index,
const std::true_type&)
1420 diagonal_ = norm_(matrix_->operator[](row_)[row_]);
1423 template<
class M,
class N>
1424 inline void SymmetricDependency<M,N>::examine(
const ColIter& col)
1427 real_type eij = norm_(*col);
1429 matrix_->operator[](col.index()).find(row_);
1430 if ( opposite_entry == matrix_->operator[](col.index()).end() )
1435 real_type eji = norm_(*opposite_entry);
1438 if(!N::is_sign_preserving || eij<0 || eji<0)
1439 maxValue_ =
max(maxValue_,
1440 eij /diagonal_ * eji/
1441 norm_(matrix_->operator[](col.index())[col.index()]));
1444 template<
class M,
class N>
1446 inline void SymmetricDependency<M,N>::examine(G& graph,
const typename G::EdgeIterator& edge,
const ColIter& col)
1448 real_type eij = norm_(*col);
1450 matrix_->operator[](col.index()).find(row_);
1452 if ( opposite_entry == matrix_->operator[](col.index()).end() )
1457 real_type eji = norm_(*opposite_entry);
1459 if(!N::is_sign_preserving || (eij<0 || eji<0))
1460 if(eji / norm_(matrix_->operator[](edge.target())[edge.target()]) *
1461 eij/ diagonal_ > alpha() * maxValue_) {
1462 edge.properties().setDepends();
1463 edge.properties().setInfluences();
1464 typename G::EdgeProperties& other = graph.getEdgeProperties(edge.target(), edge.source());
1465 other.setInfluences();
1470 template<
class M,
class N>
1471 inline bool SymmetricDependency<M,N>::isIsolated()
1473 return maxValue_ < beta();
1477 template<
class M,
class N>
1478 inline void Dependency<M,N>::init(
const Matrix* matrix)
1483 template<
class M,
class N>
1484 inline void Dependency<M,N>::initRow([[maybe_unused]]
const Row& row,
int index)
1489 diagonal_ = norm_(matrix_->operator[](row_)[row_]);
1492 template<
class M,
class N>
1493 inline void Dependency<M,N>::examine(
const ColIter& col)
1496 maxValue_ =
max(maxValue_, -norm_(*col));
1499 template<
class M,
class N>
1501 inline void Dependency<M,N>::examine(G& graph,
const typename G::EdgeIterator& edge,
const ColIter& col)
1503 if(-norm_(*col) >= maxValue_ * alpha()) {
1504 edge.properties().setDepends();
1505 typedef typename G::EdgeDescriptor ED;
1506 ED e= graph.findEdge(edge.target(), edge.source());
1509 typename G::EdgeProperties& other = graph.getEdgeProperties(e);
1510 other.setInfluences();
1515 template<
class M,
class N>
1516 inline bool Dependency<M,N>::isIsolated()
1518 return maxValue_ < beta() * diagonal_;
1521 template<
class G,
class S>
1523 VertexSet& connected, std::vector<Vertex>&
front)
1524 : vertices_(), id_(-1), graph_(graph), aggregates_(aggregates),
1525 connected_(connected), front_(
front)
1528 template<
class G,
class S>
1536 throw "Not yet implemented";
1544 template<
class G,
class S>
1547 dvverb<<
"Connected cleared"<<std::endl;
1550 connected_.insert(
vertex);
1551 dvverb <<
" Inserting "<<
vertex<<
" size="<<connected_.size();
1557 template<
class G,
class S>
1560 vertices_.insert(
vertex);
1563 front_.erase(std::lower_bound(front_.begin(), front_.end(),
vertex));
1567 const iterator end = graph_.endEdges(
vertex);
1568 for(iterator edge = graph_.beginEdges(
vertex); edge != end; ++edge) {
1569 dvverb <<
" Inserting "<<aggregates_[edge.target()];
1570 connected_.insert(aggregates_[edge.target()]);
1571 dvverb <<
" size="<<connected_.size();
1573 !graph_.getVertexProperties(edge.target()).front())
1575 front_.push_back(edge.target());
1576 graph_.getVertexProperties(edge.target()).setFront();
1580 std::sort(front_.begin(), front_.end());
1583 template<
class G,
class S>
1587 std::size_t oldsize = vertices_.size();
1589 typedef typename std::vector<Vertex>::iterator Iterator;
1591 typedef typename VertexSet::iterator SIterator;
1593 SIterator pos=vertices_.begin();
1594 std::vector<Vertex> newFront;
1595 newFront.reserve(front_.capacity());
1597 std::set_difference(front_.begin(), front_.end(), vertices.begin(), vertices.end(),
1598 std::back_inserter(newFront));
1603 pos=vertices_.insert(pos,*
vertex);
1604 vertices_.insert(*
vertex);
1605 graph_.getVertexProperties(*vertex).resetFront();
1606 aggregates_[*
vertex]=id_;
1609 const iterator end = graph_.endEdges(*
vertex);
1610 for(iterator edge = graph_.beginEdges(*
vertex); edge != end; ++edge) {
1611 dvverb <<
" Inserting "<<aggregates_[edge.target()];
1612 connected_.insert(aggregates_[edge.target()]);
1614 !graph_.getVertexProperties(edge.target()).front())
1616 front_.push_back(edge.target());
1617 graph_.getVertexProperties(edge.target()).setFront();
1619 dvverb <<
" size="<<connected_.size();
1623 std::sort(front_.begin(), front_.end());
1624 assert(oldsize+vertices.size()==vertices_.size());
1626 template<
class G,
class S>
1634 template<
class G,
class S>
1635 inline typename Aggregate<G,S>::VertexSet::size_type
1638 return vertices_.size();
1641 template<
class G,
class S>
1642 inline typename Aggregate<G,S>::VertexSet::size_type
1645 return connected_.size();
1648 template<
class G,
class S>
1654 template<
class G,
class S>
1657 return vertices_.begin();
1660 template<
class G,
class S>
1663 return vertices_.end();
1681 delete[] aggregates_;
1688 allocate(noVertices);
1701 noVertices_ = noVertices;
1703 for(std::size_t i=0; i < noVertices; i++)
1704 aggregates_[i]=UNAGGREGATED;
1710 assert(aggregates_ != 0);
1711 delete[] aggregates_;
1719 return aggregates_[v];
1726 return aggregates_[v];
1730 template<
bool reset,
class G,
class F,
class VM>
1733 const G& graph, F& aggregateVisitor,
1734 VM& visitedMap)
const
1738 DummyEdgeVisitor dummy;
1739 return breadthFirstSearch<true,reset>(start, aggregate, graph, vlist, aggregateVisitor, dummy, visitedMap);
1743 template<
bool remove,
bool reset,
class G,
class L,
class F1,
class F2,
class VM>
1748 F1& aggregateVisitor,
1749 F2& nonAggregateVisitor,
1750 VM& visitedMap)
const
1752 typedef typename L::const_iterator ListIterator;
1753 int visitedSpheres = 0;
1755 visited.push_back(start);
1756 put(visitedMap, start,
true);
1758 ListIterator current = visited.begin();
1759 ListIterator end = visited.end();
1760 std::size_t i=0, size=visited.size();
1764 while(current != end) {
1766 for(; i<size; ++current, ++i) {
1767 typedef typename G::ConstEdgeIterator EdgeIterator;
1768 const EdgeIterator endEdge = graph.endEdges(*current);
1770 for(EdgeIterator edge = graph.beginEdges(*current);
1771 edge != endEdge; ++edge) {
1773 if(aggregates_[edge.target()]==aggregate) {
1774 if(!get(visitedMap, edge.target())) {
1775 put(visitedMap, edge.target(),
true);
1776 visited.push_back(edge.target());
1777 aggregateVisitor(edge);
1780 nonAggregateVisitor(edge);
1783 end = visited.end();
1784 size = visited.size();
1790 for(current = visited.begin(); current != end; ++current)
1791 put(visitedMap, *current,
false);
1797 return visitedSpheres;
1802 : graph_(0), aggregate_(0), front_(), connected_(), size_(-1)
1811 template<
class G,
class C>
1813 const typename C::Matrix& matrix,
1814 C criterion,
bool firstlevel)
1817 typedef typename C::Matrix Matrix;
1818 typedef typename G::VertexIterator VertexIterator;
1820 criterion.init(&matrix);
1825 const Row& row = matrix[*
vertex];
1830 criterion.initRow(row, *
vertex);
1835 ColIterator end = row.end();
1836 typename FieldTraits<typename Matrix::field_type>::real_type absoffdiag=0.;
1840 for(ColIterator col = row.begin(); col != end; ++col)
1841 if(col.index()!=*
vertex) {
1842 criterion.examine(col);
1843 absoffdiag =
max(absoffdiag, Impl::asMatrix(*col).frobenius_norm());
1847 vertex.properties().setExcludedBorder();
1850 for(ColIterator col = row.begin(); col != end; ++col)
1852 criterion.examine(col);
1858 if(criterion.isIsolated()) {
1860 vertex.properties().setIsolated();
1863 auto eEnd =
vertex.end();
1864 auto col = matrix[*
vertex].begin();
1866 for(
auto edge =
vertex.begin(); edge!= eEnd; ++edge, ++col) {
1868 while(col.index()!=edge.target())
1870 criterion.examine(graph, edge, col);
1880 inline Aggregator<G>::AggregateVisitor<V>::AggregateVisitor(
const AggregatesMap<Vertex>& aggregates,
1882 : aggregates_(aggregates), aggregate_(aggregate), visitor_(&visitor)
1889 if(aggregates_[edge.target()]==aggregate_)
1890 visitor_->operator()(edge);
1895 inline void Aggregator<G>::visitAggregateNeighbours(
const Vertex&
vertex,
1897 const AggregatesMap<Vertex>& aggregates,
1901 AggregateVisitor<V> v(aggregates, aggregate, visitor);
1907 inline Aggregator<G>::Counter::Counter()
1912 inline void Aggregator<G>::Counter::increment()
1918 inline void Aggregator<G>::Counter::decrement()
1923 inline int Aggregator<G>::Counter::value()
1931 if(edge.properties().isTwoWay())
1932 Counter::increment();
1937 const AggregatesMap<Vertex>& aggregates)
const
1939 TwoWayCounter counter;
1940 visitAggregateNeighbours(
vertex, aggregate, aggregates, counter);
1941 return counter.value();
1946 const AggregatesMap<Vertex>& aggregates)
const
1948 OneWayCounter counter;
1949 visitAggregateNeighbours(
vertex, aggregate, aggregates, counter);
1950 return counter.value();
1956 if(edge.properties().isOneWay())
1957 Counter::increment();
1961 inline Aggregator<G>::ConnectivityCounter::ConnectivityCounter(
const VertexSet& connected,
1962 const AggregatesMap<Vertex>& aggregates)
1963 : Counter(), connected_(connected), aggregates_(aggregates)
1972 Counter::increment();
1974 Counter::increment();
1975 Counter::increment();
1980 inline double Aggregator<G>::connectivity(
const Vertex&
vertex,
const AggregatesMap<Vertex>& aggregates)
const
1982 ConnectivityCounter counter(connected_, aggregates);
1984 return (
double)counter.value()/noNeighbours;
1988 inline Aggregator<G>::DependencyCounter::DependencyCounter()
1995 if(edge.properties().depends())
1996 Counter::increment();
1997 if(edge.properties().influences())
1998 Counter::increment();
2002 int Aggregator<G>::unusedNeighbours(
const Vertex&
vertex,
const AggregatesMap<Vertex>& aggregates)
const
2008 std::pair<int,int> Aggregator<G>::neighbours(
const Vertex&
vertex,
2010 const AggregatesMap<Vertex>& aggregates)
const
2012 DependencyCounter unused, aggregated;
2013 typedef AggregateVisitor<DependencyCounter> CounterT;
2014 typedef std::tuple<CounterT,CounterT> CounterTuple;
2017 return std::make_pair(unused.value(), aggregated.value());
2024 DependencyCounter counter;
2025 visitAggregateNeighbours(
vertex, aggregate, aggregates, counter);
2026 return counter.value();
2030 std::size_t Aggregator<G>::distance(
const Vertex&
vertex,
const AggregatesMap<Vertex>& aggregates)
2033 typename PropertyMapTypeSelector<VertexVisitedTag,G>::Type visitedMap = get(VertexVisitedTag(), *graph_);
2035 typename AggregatesMap<Vertex>::DummyEdgeVisitor dummy;
2036 return aggregates.template breadthFirstSearch<true,true>(
vertex,
2037 aggregate_->
id(), *graph_,
2038 vlist, dummy, dummy, visitedMap);
2042 inline Aggregator<G>::FrontMarker::FrontMarker(std::vector<Vertex>&
front,
MatrixGraph& graph)
2043 : front_(
front), graph_(graph)
2049 Vertex target = edge.target();
2051 if(!graph_.getVertexProperties(target).front()) {
2052 front_.push_back(target);
2053 graph_.getVertexProperties(target).setFront();
2061 Dune::dvverb<<
" Admissible not yet implemented!"<<std::endl;
2068 Iterator vend = graph_->endEdges(
vertex);
2069 for(Iterator edge = graph_->beginEdges(
vertex); edge != vend; ++edge) {
2071 if(edge.properties().isStrong()
2072 && aggregates[edge.target()]==aggregate)
2075 Iterator edge1 = edge;
2076 for(++edge1; edge1 != vend; ++edge1) {
2078 if(edge1.properties().isStrong()
2079 && aggregates[edge.target()]==aggregate)
2084 Iterator v2end = graph_->endEdges(edge.target());
2085 for(Iterator edge2 = graph_->beginEdges(edge.target()); edge2 != v2end; ++edge2) {
2086 if(edge2.target()==edge1.target() &&
2087 edge2.properties().isStrong()) {
2103 vend = graph_->endEdges(
vertex);
2104 for(Iterator edge = graph_->beginEdges(
vertex); edge != vend; ++edge) {
2106 if(edge.properties().isStrong()
2107 && aggregates[edge.target()]==aggregate)
2110 Iterator v1end = graph_->endEdges(edge.target());
2112 for(Iterator edge1=graph_->beginEdges(edge.target()); edge1 != v1end; ++edge1) {
2114 if(edge1.properties().isStrong()
2115 && aggregates[edge1.target()]==aggregate)
2119 Iterator v2end = graph_->endEdges(
vertex);
2120 for(Iterator edge2 = graph_->beginEdges(
vertex); edge2 != v2end; ++edge2) {
2121 if(edge2.target()==edge1.target()) {
2122 if(edge2.properties().isStrong())
2139 void Aggregator<G>::unmarkFront()
2141 typedef typename std::vector<Vertex>::const_iterator Iterator;
2144 graph_->getVertexProperties(*vertex).resetFront();
2151 Aggregator<G>::nonisoNeighbourAggregate(
const Vertex&
vertex,
2152 const AggregatesMap<Vertex>& aggregates,
2153 SLList<Vertex>& neighbours)
const
2156 Iterator end=graph_->beginEdges(
vertex);
2159 for(Iterator edge=graph_->beginEdges(
vertex); edge!=end; ++edge)
2162 neighbours.push_back(aggregates[edge.target()]);
2167 inline typename G::VertexDescriptor Aggregator<G>::mergeNeighbour(
const Vertex&
vertex,
const AggregatesMap<Vertex>& aggregates)
const
2171 Iterator end = graph_->endEdges(
vertex);
2172 for(Iterator edge = graph_->beginEdges(
vertex); edge != end; ++edge) {
2174 graph_->getVertexProperties(edge.target()).isolated() == graph_->getVertexProperties(edge.source()).isolated()) {
2175 if( graph_->getVertexProperties(
vertex).isolated() ||
2176 ((edge.properties().depends() || edge.properties().influences())
2177 && admissible(
vertex, aggregates[edge.target()], aggregates)))
2178 return edge.target();
2185 Aggregator<G>::FrontNeighbourCounter::FrontNeighbourCounter(
const MatrixGraph& graph)
2186 : Counter(), graph_(graph)
2192 if(graph_.getVertexProperties(edge.target()).front())
2193 Counter::increment();
2197 int Aggregator<G>::noFrontNeighbours(
const Vertex&
vertex)
const
2199 FrontNeighbourCounter counter(*graph_);
2201 return counter.value();
2204 inline bool Aggregator<G>::connected(
const Vertex&
vertex,
2206 const AggregatesMap<Vertex>& aggregates)
const
2208 typedef typename G::ConstEdgeIterator iterator;
2209 const iterator end = graph_->endEdges(
vertex);
2210 for(iterator edge = graph_->beginEdges(
vertex); edge != end; ++edge)
2211 if(aggregates[edge.target()]==aggregate)
2216 inline bool Aggregator<G>::connected(
const Vertex&
vertex,
2217 const SLList<AggregateDescriptor>& aggregateList,
2218 const AggregatesMap<Vertex>& aggregates)
const
2221 for(Iter i=aggregateList.begin(); i!=aggregateList.end(); ++i)
2222 if(connected(
vertex, *i, aggregates))
2229 void Aggregator<G>::growIsolatedAggregate(
const Vertex& seed,
const AggregatesMap<Vertex>& aggregates,
const C& c)
2231 SLList<Vertex> connectedAggregates;
2232 nonisoNeighbourAggregate(seed, aggregates,connectedAggregates);
2234 while(aggregate_->
size()< c.minAggregateSize() && aggregate_->
connectSize() < c.maxConnectivity()) {
2236 std::size_t maxFrontNeighbours=0;
2240 typedef typename std::vector<Vertex>::const_iterator Iterator;
2243 if(distance(*
vertex, aggregates)>c.maxDistance())
2246 if(connectedAggregates.size()>0) {
2250 if(!connected(*
vertex, connectedAggregates, aggregates))
2254 double con = connectivity(*
vertex, aggregates);
2257 std::size_t frontNeighbours = noFrontNeighbours(*
vertex);
2259 if(frontNeighbours >= maxFrontNeighbours) {
2260 maxFrontNeighbours = frontNeighbours;
2263 }
else if(con > maxCon) {
2265 maxFrontNeighbours = noFrontNeighbours(*
vertex);
2273 aggregate_->
add(candidate);
2279 void Aggregator<G>::growAggregate(
const Vertex& seed,
const AggregatesMap<Vertex>& aggregates,
const C& c)
2283 std::size_t distance_ =0;
2284 while(aggregate_->
size() < c.minAggregateSize()&& distance_<c.maxDistance()) {
2285 int maxTwoCons=0, maxOneCons=0, maxNeighbours=-1;
2288 std::vector<Vertex> candidates;
2289 candidates.reserve(30);
2291 typedef typename std::vector<Vertex>::const_iterator Iterator;
2295 if(graph_->getVertexProperties(*vertex).isolated())
2298 int twoWayCons = twoWayConnections(*
vertex, aggregate_->
id(), aggregates);
2301 if( maxTwoCons == twoWayCons && twoWayCons > 0) {
2302 double con = connectivity(*
vertex, aggregates);
2305 int neighbours = noFrontNeighbours(*
vertex);
2307 if(neighbours > maxNeighbours) {
2308 maxNeighbours = neighbours;
2310 candidates.push_back(*
vertex);
2312 candidates.push_back(*
vertex);
2314 }
else if( con > maxCon) {
2316 maxNeighbours = noFrontNeighbours(*
vertex);
2318 candidates.push_back(*
vertex);
2320 }
else if(twoWayCons > maxTwoCons) {
2321 maxTwoCons = twoWayCons;
2322 maxCon = connectivity(*
vertex, aggregates);
2323 maxNeighbours = noFrontNeighbours(*
vertex);
2325 candidates.push_back(*
vertex);
2337 int oneWayCons = oneWayConnections(*
vertex, aggregate_->
id(), aggregates);
2342 if(!admissible(*
vertex, aggregate_->
id(), aggregates))
2345 if( maxOneCons == oneWayCons && oneWayCons > 0) {
2346 double con = connectivity(*
vertex, aggregates);
2349 int neighbours = noFrontNeighbours(*
vertex);
2351 if(neighbours > maxNeighbours) {
2352 maxNeighbours = neighbours;
2354 candidates.push_back(*
vertex);
2356 if(neighbours==maxNeighbours)
2358 candidates.push_back(*
vertex);
2361 }
else if( con > maxCon) {
2363 maxNeighbours = noFrontNeighbours(*
vertex);
2365 candidates.push_back(*
vertex);
2367 }
else if(oneWayCons > maxOneCons) {
2368 maxOneCons = oneWayCons;
2369 maxCon = connectivity(*
vertex, aggregates);
2370 maxNeighbours = noFrontNeighbours(*
vertex);
2372 candidates.push_back(*
vertex);
2377 if(!candidates.size())
2379 distance_=distance(seed, aggregates);
2380 candidates.resize(
min(candidates.size(), c.maxAggregateSize()-
2381 aggregate_->
size()));
2382 aggregate_->
add(candidates);
2386 template<
typename V>
2387 template<
typename M,
typename G,
typename C>
2391 Aggregator<G> aggregator;
2392 return aggregator.build(matrix, graph, *
this, criterion, finestLevel);
2396 template<
class M,
class C>
2397 std::tuple<int,int,int,int>
Aggregator<G>::build(
const M& m, G& graph, AggregatesMap<Vertex>& aggregates,
const C& c,
2403 Stack stack_(graph, *
this, aggregates);
2407 aggregate_ =
new Aggregate<G,VertexSet>(graph, aggregates, connected_, front_);
2414 dverb<<
"Build dependency took "<< watch.elapsed()<<
" seconds."<<std::endl;
2415 int noAggregates, conAggregates, isoAggregates, oneAggregates;
2416 std::size_t maxA=0, minA=1000000, avg=0;
2417 int skippedAggregates;
2418 noAggregates = conAggregates = isoAggregates = oneAggregates =
2419 skippedAggregates = 0;
2422 Vertex seed = stack_.pop();
2424 if(seed == Stack::NullEntry)
2429 if((noAggregates+1)%10000 == 0)
2433 if(graph.getVertexProperties(seed).excludedBorder()) {
2435 ++skippedAggregates;
2439 if(graph.getVertexProperties(seed).isolated()) {
2440 if(c.skipIsolated()) {
2443 ++skippedAggregates;
2447 aggregate_->
seed(seed);
2448 growIsolatedAggregate(seed, aggregates, c);
2451 aggregate_->
seed(seed);
2452 growAggregate(seed, aggregates, c);
2456 while(!(graph.getVertexProperties(seed).isolated()) && aggregate_->
size() < c.maxAggregateSize()) {
2458 std::vector<Vertex> candidates;
2459 candidates.reserve(30);
2461 typedef typename std::vector<Vertex>::const_iterator Iterator;
2465 if(graph.getVertexProperties(*vertex).isolated())
2468 if(twoWayConnections( *
vertex, aggregate_->
id(), aggregates) == 0 &&
2469 (oneWayConnections( *
vertex, aggregate_->
id(), aggregates) == 0 ||
2470 !admissible( *
vertex, aggregate_->
id(), aggregates) ))
2473 std::pair<int,int> neighbourPair=neighbours(*
vertex, aggregate_->
id(),
2479 if(neighbourPair.first >= neighbourPair.second)
2482 if(distance(*
vertex, aggregates) > c.maxDistance())
2484 candidates.push_back(*
vertex);
2488 if(!candidates.size())
break;
2490 candidates.resize(
min(candidates.size(), c.maxAggregateSize()-
2491 aggregate_->
size()));
2492 aggregate_->
add(candidates);
2497 if(aggregate_->
size()==1 && c.maxAggregateSize()>1) {
2498 if(!graph.getVertexProperties(seed).isolated()) {
2499 Vertex mergedNeighbour = mergeNeighbour(seed, aggregates);
2503 aggregates[seed] = aggregates[mergedNeighbour];
2504 aggregate_->invalidate();
2507 minA=
min(minA,
static_cast<std::size_t
>(1));
2508 maxA=
max(maxA,
static_cast<std::size_t
>(1));
2514 minA=
min(minA,
static_cast<std::size_t
>(1));
2515 maxA=
max(maxA,
static_cast<std::size_t
>(1));
2521 avg+=aggregate_->
size();
2522 minA=
min(minA,aggregate_->
size());
2523 maxA=
max(maxA,aggregate_->
size());
2524 if(graph.getVertexProperties(seed).isolated())
2532 Dune::dinfo<<
"connected aggregates: "<<conAggregates;
2533 Dune::dinfo<<
" isolated aggregates: "<<isoAggregates;
2534 if(conAggregates+isoAggregates>0)
2535 Dune::dinfo<<
" one node aggregates: "<<oneAggregates<<
" min size="
2536 <<minA<<
" max size="<<maxA
2537 <<
" avg="<<avg/(conAggregates+isoAggregates)<<std::endl;
2540 return std::make_tuple(conAggregates+isoAggregates,isoAggregates,
2541 oneAggregates,skippedAggregates);
2546 Aggregator<G>::Stack::Stack(
const MatrixGraph& graph,
const Aggregator<G>& aggregatesBuilder,
2547 const AggregatesMap<Vertex>& aggregates)
2548 : graph_(graph), aggregatesBuilder_(aggregatesBuilder), aggregates_(aggregates), begin_(graph.begin()), end_(graph.end())
2554 Aggregator<G>::Stack::~Stack()
2565 inline typename G::VertexDescriptor Aggregator<G>::Stack::pop()
2571 typename G::VertexDescriptor current=*begin_;
2581 void printAggregates2d(
const AggregatesMap<V>& aggregates,
int n,
int m, std::ostream& os)
2585 std::ios_base::fmtflags oldOpts=os.flags();
2587 os.setf(std::ios_base::right, std::ios_base::adjustfield);
2592 for(
int i=0; i< n*m; i++)
2593 maxVal=
max(maxVal, aggregates[i]);
2595 for(
int i=10; i < 1000000; i*=10)
2601 for(
int j=0, entry=0; j < m; j++) {
2602 for(
int i=0; i<n; i++, entry++) {
2604 os<<aggregates[entry]<<
" ";
A class for temporarily storing the vertices of an aggregate in.
Definition: aggregates.hh:776
A Dummy visitor that does nothing for each visited edge.
Definition: aggregates.hh:596
Class providing information about the mapping of the vertices onto aggregates.
Definition: aggregates.hh:558
Base class of all aggregation criterions.
Definition: aggregates.hh:47
Class for building the aggregates.
Definition: aggregates.hh:907
Dependency policy for symmetric matrices.
Definition: aggregates.hh:251
Norm that uses only the [N][N] entry of the block to determine couplings.
Definition: aggregates.hh:377
Norm that uses only the [0][0] entry of the block to determine couplings.
Definition: aggregates.hh:453
Iterator over all edges starting from a vertex.
Definition: graph.hh:93
The vertex iterator type of the graph.
Definition: graph.hh:207
The (undirected) graph of a matrix.
Definition: graph.hh:49
M::size_type VertexDescriptor
The vertex descriptor.
Definition: graph.hh:71
EdgeIteratorT< const MatrixGraph< Matrix > > ConstEdgeIterator
The constant edge iterator type.
Definition: graph.hh:296
All parameters for AMG.
Definition: parameters.hh:391
Criterion taking advantage of symmetric matrices.
Definition: aggregates.hh:517
Dependency policy for symmetric matrices.
Definition: aggregates.hh:312
Dependency policy for symmetric matrices.
Definition: aggregates.hh:132
Criterion suitable for unsymmetric matrices.
Definition: aggregates.hh:537
derive error class from the base class in common
Definition: istlexception.hh:17
A generic dynamic dense matrix.
Definition: matrix.hh:559
typename Imp::BlockTraits< T >::field_type field_type
Export the type representing the underlying field.
Definition: matrix.hh:563
row_type::const_iterator ConstColIterator
Const iterator for the entries of each row.
Definition: matrix.hh:587
MatrixImp::DenseMatrixBase< T, A >::window_type row_type
The type implementing a matrix row.
Definition: matrix.hh:572
An allocator managing a pool of objects for reuse.
Definition: poolallocator.hh:223
A single linked list.
Definition: sllist.hh:42
Type traits to determine the type of reals (when working with complex numbers)
Provides classes for building the matrix graph.
SLListConstIterator< T, A > const_iterator
The constant iterator of the list.
Definition: sllist.hh:72
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
constexpr GeometryType vertex
GeometryType representing a vertex.
Definition: type.hh:504
Matrix::ConstColIterator ColIter
Constant column iterator of the matrix.
Definition: aggregates.hh:271
Matrix::ConstColIterator ColIter
Constant column iterator of the matrix.
Definition: aggregates.hh:152
std::size_t breadthFirstSearch(const VertexDescriptor &start, const AggregateDescriptor &aggregate, const G &graph, L &visited, F1 &aggregateVisitor, F2 &nonAggregateVisitor, VM &visitedMap) const
Breadth first search within an aggregate.
PoolAllocator< VertexDescriptor, 100 > Allocator
The allocator we use for our lists and the set.
Definition: aggregates.hh:584
int id()
Get the id identifying the aggregate.
Norm norm_
The functor for calculating the norm.
Definition: aggregates.hh:300
MatrixGraph::VertexDescriptor Vertex
The vertex identifier.
Definition: aggregates.hh:918
AggregationCriterion()
Constructor.
Definition: aggregates.hh:64
const Matrix * matrix_
The matrix we work on.
Definition: aggregates.hh:355
auto operator()(const M &m, typename std::enable_if_t< Dune::IsNumber< M >::value > *sfinae=nullptr) const
Compute the norm of a scalar.
Definition: aggregates.hh:404
M Matrix
The matrix type we build the dependency of.
Definition: aggregates.hh:256
G MatrixGraph
The matrix graph type used.
Definition: aggregates.hh:913
Norm norm_
The functor for calculating the norm.
Definition: aggregates.hh:361
M Matrix
The matrix type we build the dependency of.
Definition: aggregates.hh:317
real_type diagonal_
The norm of the current diagonal.
Definition: aggregates.hh:185
N Norm
The norm to use for examining the matrix entries.
Definition: aggregates.hh:261
int row_
index of the currently evaluated row.
Definition: aggregates.hh:183
std::tuple< int, int, int, int > build(const M &m, G &graph, AggregatesMap< Vertex > &aggregates, const C &c, bool finestLevel)
Build the aggregates.
FrontNeighbourCounter(const MatrixGraph &front)
Constructor.
Matrix::row_type Row
Constant Row iterator of the matrix.
Definition: aggregates.hh:327
const Matrix * matrix_
The matrix we work on.
Definition: aggregates.hh:294
const AggregateDescriptor & operator[](const VertexDescriptor &v) const
Get the aggregate a vertex belongs to.
AggregateVisitor(const AggregatesMap< Vertex > &aggregates, const AggregateDescriptor &aggregate, Visitor &visitor)
Constructor.
Matrix::ConstColIterator ColIter
Constant column iterator of the matrix.
Definition: aggregates.hh:332
~AggregatesMap()
Destructor.
Matrix::field_type field_type
The current max value.
Definition: aggregates.hh:177
void decrement()
Decrement counter.
Aggregate(MatrixGraph &graph, AggregatesMap< Vertex > &aggregates, VertexSet &connectivity, std::vector< Vertex > &front_)
Constructor.
V Visitor
The type of the adapted visitor.
Definition: aggregates.hh:1062
std::size_t * SphereMap
Type of the mapping of aggregate members onto distance spheres.
Definition: aggregates.hh:807
Matrix::row_type Row
Constant Row iterator of the matrix.
Definition: aggregates.hh:266
VertexSet::size_type connectSize()
Get tne number of connections to other aggregates.
N Norm
The norm to use for examining the matrix entries.
Definition: aggregates.hh:322
Norm norm_
The functor for calculating the norm.
Definition: aggregates.hh:181
VertexSet::const_iterator const_iterator
Const iterator over a vertex list.
Definition: aggregates.hh:802
MatrixGraph::VertexDescriptor AggregateDescriptor
The type of the aggregate descriptor.
Definition: aggregates.hh:921
real_type diagonal_
The norm of the current diagonal.
Definition: aggregates.hh:365
void add(const Vertex &vertex)
Add a vertex to the aggregate.
T DependencyPolicy
The policy for calculating the dependency graph.
Definition: aggregates.hh:53
void operator()(const typename MatrixGraph::ConstEdgeIterator &edge)
Examine an edge.
FrontMarker(std::vector< Vertex > &front, MatrixGraph &graph)
Constructor.
FieldTraits< typenameM::field_type >::real_type operator()(const M &m) const
compute the norm of a matrix.
Definition: aggregates.hh:504
int visitNeighbours(const G &graph, const typename G::VertexDescriptor &vertex, V &visitor)
Visit all neighbour vertices of a vertex in a graph.
void setDefaultValuesIsotropic(std::size_t dim, std::size_t diameter=2)
Sets reasonable default values for an isotropic problem.
Definition: aggregates.hh:80
V AggregateDescriptor
The aggregate descriptor type.
Definition: aggregates.hh:578
static const V ISOLATED
Identifier of isolated vertices.
Definition: aggregates.hh:569
int row_
index of the currently evaluated row.
Definition: aggregates.hh:363
DependencyCounter()
Constructor.
real_type diagonal_
The norm of the current diagonal.
Definition: aggregates.hh:304
std::size_t noVertices() const
Get the number of vertices.
void setDefaultValuesAnisotropic(std::size_t dim, std::size_t diameter=2)
Sets reasonable default values for an aisotropic problem.
Definition: aggregates.hh:103
AggregatesMap(std::size_t noVertices)
Constructs with allocating memory.
Matrix::field_type field_type
The current max value.
Definition: aggregates.hh:357
AggregateDescriptor & operator[](const VertexDescriptor &v)
Get the aggregate a vertex belongs to.
AggregatesMap()
Constructs without allocating memory.
int value()
Access the current count.
SLList< VertexDescriptor, Allocator > VertexList
The type of a single linked list of vertex descriptors.
Definition: aggregates.hh:590
ConnectivityCounter(const VertexSet &connected, const AggregatesMap< Vertex > &aggregates)
Constructor.
VertexSet::size_type size()
Get the size of the aggregate.
const_iterator end() const
get an iterator over the vertices of the aggregate.
int row_
index of the currently evaluated row.
Definition: aggregates.hh:302
M Matrix
The matrix type we build the dependency of.
Definition: aggregates.hh:137
const Matrix * matrix_
The matrix we work on.
Definition: aggregates.hh:175
S VertexSet
The type of a single linked list of vertex descriptors.
Definition: aggregates.hh:799
FieldTraits< typenameM::field_type >::real_type operator()(const M &m) const
compute the norm of a matrix.
Definition: aggregates.hh:488
static const V UNAGGREGATED
Identifier of not yet aggregated vertices.
Definition: aggregates.hh:564
std::size_t breadthFirstSearch(const VertexDescriptor &start, const AggregateDescriptor &aggregate, const G &graph, F &aggregateVisitor, VM &visitedMap) const
Breadth first search within an aggregate.
Matrix::field_type field_type
The current max value.
Definition: aggregates.hh:296
void allocate(std::size_t noVertices)
Allocate memory for holding the information.
N Norm
The norm to use for examining the matrix entries.
Definition: aggregates.hh:142
FieldTraits< typenameM::field_type >::real_type operator()(const M &m) const
compute the norm of a matrix.
Definition: aggregates.hh:471
void reconstruct(const Vertex &vertex)
Reconstruct the aggregat from an seed node.
const_iterator begin() const
get an iterator over the vertices of the aggregate.
FieldTraits< typenameM::field_type >::real_type operator()(const M &m, typename std::enable_if_t<!Dune::IsNumber< M >::value > *sfinae=nullptr) const
compute the norm of a matrix.
Definition: aggregates.hh:388
MatrixGraph::VertexDescriptor Vertex
The vertex descriptor type.
Definition: aggregates.hh:787
void seed(const Vertex &vertex)
Initialize the aggregate with one vertex.
void clear()
Clear the aggregate.
void free()
Free the allocated memory.
void increment()
Increment counter.
void buildDependency(G &graph, const typename C::Matrix &matrix, C criterion, bool finestLevel)
Build the dependency of the matrix graph.
V VertexDescriptor
The vertex descriptor type.
Definition: aggregates.hh:573
std::tuple< int, int, int, int > buildAggregates(const M &matrix, G &graph, const C &criterion, bool finestLevel)
Build the aggregates.
Matrix::row_type Row
Constant Row iterator of the matrix.
Definition: aggregates.hh:147
PoolAllocator< Vertex, 100 > Allocator
The allocator we use for our lists and the set.
Definition: aggregates.hh:793
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
constexpr auto front(const HybridTreePath< T... > &tp) -> decltype(treePathEntry< 0 >(tp))
Returns a copy of the first element of the HybridTreePath.
Definition: treepath.hh:267
Dune namespace.
Definition: alignedallocator.hh:11
Parameter classes for customizing AMG.
An stl-compliant pool allocator.
Provides classes for handling internal properties in a graph.
Implements a scalar matrix view wrapper around an existing scalar.
Implements a singly linked list together with the necessary iterators.
Standard Dune debug streams.
Functor using the row sum (infinity) norm to determine strong couplings.
Definition: aggregates.hh:461
Whether this type acts as a scalar in the context of (hierarchically blocked) containers.
Definition: typetraits.hh:194