3 #ifndef DUNE_AMG_AGGREGATES_HH
4 #define DUNE_AMG_AGGREGATES_HH
12 #include <dune/common/timer.hh>
13 #include <dune/common/stdstreams.hh>
14 #include <dune/common/poolallocator.hh>
15 #include <dune/common/sllist.hh>
16 #include <dune/common/ftraits.hh>
17 #include <dune/common/scalarmatrixview.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);
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>
161 void examine(G& graph,
const typename G::EdgeIterator& edge,
const ColIter&
col);
178 typedef typename FieldTraits<field_type>::real_type
real_type;
187 typename std::vector<real_type>::iterator
valIter_;
192 template<
class M,
class N>
198 template<
class M,
class N>
202 vals_.assign(row.size(), 0.0);
203 assert(vals_.size()==row.size());
204 valIter_=vals_.begin();
206 maxValue_ = min(- std::numeric_limits<real_type>::max(), std::numeric_limits<real_type>::min());
207 diagonal_=norm_(row[index]);
211 template<
class M,
class N>
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>
230 if(*valIter_ > alpha() * maxValue_) {
231 edge.properties().setDepends();
232 edge.properties().setInfluences();
237 template<
class M,
class N>
242 valIter_=vals_.begin();
243 return maxValue_ < beta();
249 template<
class M,
class N>
297 typedef typename FieldTraits<field_type>::real_type
real_type;
310 template<
class M,
class N>
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&);
388 typename FieldTraits<typename M::field_type>::real_type
operator()(
const M& m,
389 [[maybe_unused]]
typename std::enable_if_t<!Dune::IsNumber<M>::value>* sfinae =
nullptr)
const
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");
405 typename std::enable_if_t<Dune::IsNumber<M>::value>* sfinae =
nullptr)
const
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());
471 typename FieldTraits<typename M::field_type>::real_type
operator()(
const M& m)
const
473 return m.infinity_norm();
488 typename FieldTraits<typename M::field_type>::real_type
operator()(
const M& m)
const
490 return m.frobenius_norm();
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>
546 template<
class G>
class Aggregator;
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;
758 std::size_t noVertices_;
764 template<
class G,
class C>
766 const typename C::Matrix& matrix,
774 template<
class G,
class S>
818 VertexSet& connectivity, std::vector<Vertex>& front_);
843 void add(std::vector<Vertex>& vertex);
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,
958 typedef PoolAllocator<Vertex,100> Allocator;
963 typedef SLList<Vertex,Allocator> VertexList;
968 typedef std::set<Vertex,std::less<Vertex>,Allocator> VertexSet;
973 typedef std::size_t* SphereMap;
988 std::vector<Vertex> front_;
993 VertexSet connected_;
1014 enum { N = 1300000 };
1056 class AggregateVisitor
1116 class FrontNeighbourCounter :
public Counter
1135 int noFrontNeighbours(
const Vertex& vertex)
const;
1140 class TwoWayCounter :
public Counter
1163 class OneWayCounter :
public Counter
1189 class ConnectivityCounter :
public Counter
1204 const VertexSet& connected_;
1239 bool connected(
const Vertex& vertex,
const SLList<AggregateDescriptor>& aggregateList,
1249 class DependencyCounter :
public Counter
1281 std::vector<Vertex>& front_;
1320 std::pair<int,int> neighbours(
const Vertex& vertex,
1377 void nonisoNeighbourAggregate(
const Vertex& vertex,
1379 SLList<Vertex>& neighbours)
const;
1396 template<
class M,
class N>
1402 template<
class M,
class N>
1405 initRow(row, index, std::is_convertible<field_type, real_type>());
1408 template<
class M,
class N>
1411 DUNE_THROW(InvalidStateException,
"field_type needs to convertible to real_type");
1414 template<
class M,
class N>
1418 maxValue_ = min(- std::numeric_limits<typename Matrix::field_type>::max(), std::numeric_limits<typename Matrix::field_type>::min());
1420 diagonal_ = norm_(matrix_->operator[](row_)[row_]);
1423 template<
class M,
class N>
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>
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>
1473 return maxValue_ < beta();
1477 template<
class M,
class N>
1483 template<
class M,
class N>
1487 maxValue_ = min(- std::numeric_limits<real_type>::max(), std::numeric_limits<real_type>::min());
1489 diagonal_ = norm_(matrix_->operator[](row_)[row_]);
1492 template<
class M,
class N>
1496 maxValue_ = max(maxValue_, -norm_(*
col));
1499 template<
class M,
class N>
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());
1507 if(e!=std::numeric_limits<ED>::max())
1509 typename G::EdgeProperties& other = graph.getEdgeProperties(e);
1510 other.setInfluences();
1515 template<
class M,
class N>
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);
1561 aggregates_[vertex]=id_;
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));
1601 for(Iterator vertex=vertices.begin(); vertex != vertices.end(); ++vertex)
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);
1822 for(VertexIterator vertex = graph.begin(); vertex != graph.end(); ++vertex) {
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)
1851 if(
col.index()!=*vertex)
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());
2022 int Aggregator<G>::aggregateNeighbours(
const Vertex& vertex,
const AggregateDescriptor& aggregate,
const AggregatesMap<Vertex>& aggregates)
const
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();
2058 inline bool Aggregator<G>::admissible(
const Vertex& vertex,
const AggregateDescriptor& aggregate,
const AggregatesMap<Vertex>& aggregates)
const
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;
2143 for(Iterator vertex=front_.begin(); vertex != front_.end(); ++vertex)
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
2220 typedef typename SLList<AggregateDescriptor>::const_iterator Iter;
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;
2242 for(Iterator vertex = front_.begin(); vertex != front_.end(); ++vertex) {
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;
2261 candidate = *vertex;
2263 }
else if(con > maxCon) {
2265 maxFrontNeighbours = noFrontNeighbours(*vertex);
2266 candidate = *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;
2293 for(Iterator vertex = front_.begin(); vertex != front_.end(); ++vertex) {
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);
2328 maxOneCons = std::numeric_limits<int>::max();
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;
2463 for(Iterator vertex = front_.begin(); vertex != front_.end(); ++vertex) {
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];
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()
2562 = std::numeric_limits<typename G::VertexDescriptor>::max();
2565 inline typename G::VertexDescriptor Aggregator<G>::Stack::pop()
2571 typename G::VertexDescriptor current=*begin_;
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]<<
" ";
Provides classes for building the matrix graph.
Parameter classes for customizing AMG.
Provides classes for handling internal properties in a graph.
Col col
Definition: matrixmatrix.hh:349
Matrix::ConstColIterator ColIter
Constant column iterator of the matrix.
Definition: aggregates.hh:271
std::vector< real_type >::iterator valIter_
Definition: aggregates.hh:187
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
iterator begin()
Definition: aggregates.hh:735
int id()
Get the id identifying the aggregate.
Norm norm_
The functor for calculating the norm.
Definition: aggregates.hh:300
void operator()([[maybe_unused]] const EdgeIterator &edge) const
Definition: aggregates.hh:599
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
void initRow(const Row &row, int index)
SymmetricMatrixDependency(const Parameters &parms)
Definition: aggregates.hh:166
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
FieldTraits< typename M::field_type >::real_type operator()(const M &m) const
compute the norm of a matrix.
Definition: aggregates.hh:488
const AggregateDescriptor & operator[](const VertexDescriptor &v) const
Get the aggregate a vertex belongs to.
Norm norm_
The functor for calculating the norm.
Definition: aggregates.hh:361
SymmetricCriterion()
Definition: aggregates.hh:522
Dependency()
Definition: aggregates.hh:288
void examine(const ColIter &col)
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
iterator end()
Definition: aggregates.hh:740
UnSymmetricCriterion(const Parameters &parms)
Definition: aggregates.hh:539
void operator()(const typename MatrixGraph::ConstEdgeIterator &edge)
void initRow(const Row &row, int index)
Definition: aggregates.hh:199
static const Vertex NullEntry
Definition: aggregates.hh:1006
std::tuple< int, int, int, int > build(const M &m, G &graph, AggregatesMap< Vertex > &aggregates, const C &c, bool finestLevel)
Build the aggregates.
std::ostream & operator<<(std::ostream &os, const AggregationCriterion< T > &criterion)
Definition: aggregates.hh:111
void examine(const ColIter &col)
Definition: aggregates.hh:212
Dependency(const Parameters &parms)
Definition: aggregates.hh:284
int row_
index of the currently evaluated row.
Definition: aggregates.hh:183
friend class Stack
Definition: aggregates.hh:1034
void operator()(const typename MatrixGraph::ConstEdgeIterator &edge)
FrontNeighbourCounter(const MatrixGraph &front)
Constructor.
Matrix::row_type Row
Constant Row iterator of the matrix.
Definition: aggregates.hh:327
FieldTraits< typename M::field_type >::real_type operator()(const M &m, [[maybe_unused]] typename std::enable_if_t<!Dune::IsNumber< M >::value > *sfinae=nullptr) const
compute the norm of a matrix.
Definition: aggregates.hh:388
const Matrix * matrix_
The matrix we work on.
Definition: aggregates.hh:294
void examine(G &graph, const typename G::EdgeIterator &edge, const ColIter &col)
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
AggregationCriterion(const Parameters &parms)
Definition: aggregates.hh:68
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.
std::vector< real_type > vals_
Definition: aggregates.hh:186
FieldTraits< typename M::field_type >::real_type operator()(const M &m) const
compute the norm of a matrix.
Definition: aggregates.hh:504
N Norm
The norm to use for examining the matrix entries.
Definition: aggregates.hh:322
void printAggregates2d(const AggregatesMap< V > &aggregates, int n, int m, std::ostream &os)
Definition: aggregates.hh:2581
void invalidate()
Definition: aggregates.hh:820
const_iterator begin() const
Definition: aggregates.hh:723
real_type maxValue_
Definition: aggregates.hh:179
Norm norm_
The functor for calculating the norm.
Definition: aggregates.hh:181
void operator()(const typename MatrixGraph::ConstEdgeIterator &edge)
VertexSet::const_iterator const_iterator
Const iterator over a vertex list.
Definition: aggregates.hh:802
SymmetricCriterion(const Parameters &parms)
Definition: aggregates.hh:519
MatrixGraph::VertexDescriptor AggregateDescriptor
The type of the aggregate descriptor.
Definition: aggregates.hh:921
real_type maxValue_
Definition: aggregates.hh:359
void init(const Matrix *matrix)
Definition: aggregates.hh:193
real_type diagonal_
The norm of the current diagonal.
Definition: aggregates.hh:365
AggregateDescriptor * iterator
Definition: aggregates.hh:733
SymmetricDependency(const Parameters &parms)
Definition: aggregates.hh:346
FieldTraits< field_type >::real_type real_type
Definition: aggregates.hh:358
void examine(G &graph, const typename G::EdgeIterator &edge, const ColIter &col)
AggregateDescriptor & operator[](const VertexDescriptor &v)
Get the aggregate a vertex belongs to.
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 add(std::vector< Vertex > &vertex)
void examine(const ColIter &col)
void operator()(const typename MatrixGraph::ConstEdgeIterator &edge)
Examine an edge.
FrontMarker(std::vector< Vertex > &front, MatrixGraph &graph)
Constructor.
void init(const Matrix *matrix)
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
const_iterator end() const
Definition: aggregates.hh:728
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
SymmetricMatrixDependency()
Definition: aggregates.hh:169
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
AggregatesMap()
Constructs without allocating memory.
int value()
Access the current count.
std::tuple< int, int, int, int > buildAggregates(const M &matrix, G &graph, const C &criterion, bool finestLevel)
Build the aggregates.
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.
void init(const Matrix *matrix)
UnSymmetricCriterion()
Definition: aggregates.hh:542
void operator()(const typename MatrixGraph::ConstEdgeIterator &edge)
FieldTraits< typename M::field_type >::real_type operator()(const M &m) const
compute the norm of a matrix.
Definition: aggregates.hh:471
const AggregateDescriptor * const_iterator
Definition: aggregates.hh:721
int row_
index of the currently evaluated row.
Definition: aggregates.hh:302
Stack(const MatrixGraph &graph, const Aggregator< G > &aggregatesBuilder, const AggregatesMap< Vertex > &aggregates)
FieldTraits< field_type >::real_type real_type
Definition: aggregates.hh:178
void initRow(const Row &row, int index)
M Matrix
The matrix type we build the dependency of.
Definition: aggregates.hh:137
FieldTraits< field_type >::real_type real_type
Definition: aggregates.hh:297
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
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 operator()(const typename MatrixGraph::ConstEdgeIterator &edge)
bool isIsolated()
Definition: aggregates.hh:238
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
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.
MatrixGraph::VertexDescriptor Vertex
The vertex descriptor type.
Definition: aggregates.hh:787
void seed(const Vertex &vertex)
Initialize the aggregate with one vertex.
SymmetricDependency()
Definition: aggregates.hh:349
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
real_type maxValue_
Definition: aggregates.hh:298
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
G MatrixGraph
Definition: aggregates.hh:783
void operator()(const typename MatrixGraph::ConstEdgeIterator &edge)
@ is_sign_preserving
Definition: aggregates.hh:481
@ is_sign_preserving
Definition: aggregates.hh:380
@ is_sign_preserving
Definition: aggregates.hh:497
@ is_sign_preserving
Definition: aggregates.hh:464
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
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
Base class of all aggregation criterions.
Definition: aggregates.hh:47
Dependency policy for symmetric matrices.
Definition: aggregates.hh:132
Dependency policy for symmetric matrices.
Definition: aggregates.hh:251
Dependency policy for symmetric matrices.
Definition: aggregates.hh:312
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
Functor using the row sum (infinity) norm to determine strong couplings.
Definition: aggregates.hh:461
Definition: aggregates.hh:478
Definition: aggregates.hh:494
Criterion taking advantage of symmetric matrices.
Definition: aggregates.hh:517
Criterion suitable for unsymmetric matrices.
Definition: aggregates.hh:537
Class for building the aggregates.
Definition: aggregates.hh:907
Class providing information about the mapping of the vertices onto aggregates.
Definition: aggregates.hh:558
A Dummy visitor that does nothing for each visited edge.
Definition: aggregates.hh:596
A class for temporarily storing the vertices of an aggregate in.
Definition: aggregates.hh:776
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
Iterator over all edges starting from a vertex.
Definition: graph.hh:93
The vertex iterator type of the graph.
Definition: graph.hh:207
All parameters for AMG.
Definition: parameters.hh:391