5 #ifndef DUNE_AMG_AGGREGATES_HH
6 #define DUNE_AMG_AGGREGATES_HH
14 #include <dune/common/timer.hh>
15 #include <dune/common/stdstreams.hh>
16 #include <dune/common/poolallocator.hh>
17 #include <dune/common/sllist.hh>
18 #include <dune/common/ftraits.hh>
19 #include <dune/common/scalarmatrixview.hh>
84 this->setMaxDistance(diameter-1);
89 this->setMaxDistance(this->maxDistance()+diameter-1);
91 this->setMinAggregateSize(csize);
92 this->setMaxAggregateSize(
static_cast<std::size_t
>(csize*1.5));
108 this->setMaxDistance(this->maxDistance()+dim-1);
115 os<<
"{ maxdistance="<<criterion.maxDistance()<<
" minAggregateSize="
116 <<criterion.minAggregateSize()<<
" maxAggregateSize="<<criterion.maxAggregateSize()
117 <<
" connectivity="<<criterion.maxConnectivity()<<
" debugLevel="<<criterion.debugLevel()<<
"}";
132 template<
class M,
class N>
163 void examine(G& graph,
const typename G::EdgeIterator& edge,
const ColIter&
col);
180 typedef typename FieldTraits<field_type>::real_type
real_type;
189 typename std::vector<real_type>::iterator
valIter_;
194 template<
class M,
class N>
200 template<
class M,
class N>
204 vals_.assign(row.size(), 0.0);
205 assert(vals_.size()==row.size());
206 valIter_=vals_.begin();
208 maxValue_ = min(- std::numeric_limits<real_type>::max(), std::numeric_limits<real_type>::min());
209 diagonal_=norm_(row[index]);
213 template<
class M,
class N>
219 if(!N::is_sign_preserving || eij<0)
221 *valIter_ = eij/diagonal_*eij/norm_(matrix_->operator[](
col.index())[
col.index()]);
222 maxValue_ = max(maxValue_, *valIter_);
228 template<
class M,
class N>
232 if(*valIter_ > alpha() * maxValue_) {
233 edge.properties().setDepends();
234 edge.properties().setInfluences();
239 template<
class M,
class N>
244 valIter_=vals_.begin();
245 return maxValue_ < beta();
251 template<
class M,
class N>
299 typedef typename FieldTraits<field_type>::real_type
real_type;
312 template<
class M,
class N>
360 typedef typename FieldTraits<field_type>::real_type
real_type;
369 void initRow(
const Row& row,
int index,
const std::true_type&);
370 void initRow(
const Row& row,
int index,
const std::false_type&);
390 typename FieldTraits<typename M::field_type>::real_type
operator()(
const M& m,
391 [[maybe_unused]]
typename std::enable_if_t<!Dune::IsNumber<M>::value>* sfinae =
nullptr)
const
393 typedef typename M::field_type field_type;
394 typedef typename FieldTraits<field_type>::real_type real_type;
395 static_assert( std::is_convertible<field_type, real_type >::value,
396 "use of diagonal norm in AMG not implemented for complex field_type");
407 typename std::enable_if_t<Dune::IsNumber<M>::value>* sfinae =
nullptr)
const
409 typedef typename FieldTraits<M>::real_type real_type;
410 static_assert( std::is_convertible<M, real_type >::value,
411 "use of diagonal norm in AMG not implemented for complex field_type");
420 static T signed_abs(
const T & v)
427 static T signed_abs(
const std::complex<T> & v)
431 return csgn(v) * std::abs(v);
436 static T csgn(
const T & v)
438 return (T(0) < v) - (v < T(0));
443 static T csgn(std::complex<T> a)
445 return csgn(a.real())+(a.real() == 0.0)*csgn(a.imag());
473 typename FieldTraits<typename M::field_type>::real_type
operator()(
const M& m)
const
475 return m.infinity_norm();
490 typename FieldTraits<typename M::field_type>::real_type
operator()(
const M& m)
const
492 return m.frobenius_norm();
506 typename FieldTraits<typename M::field_type>::real_type
operator()(
const M& m)
const
517 template<
class M,
class Norm>
537 template<
class M,
class Norm>
548 template<
class G>
class Aggregator;
600 template<
class EdgeIterator>
601 void operator()([[maybe_unused]]
const EdgeIterator& edge)
const
634 template<
class M,
class G,
class C>
635 std::tuple<int,int,int,int>
buildAggregates(
const M& matrix, G& graph,
const C& criterion,
655 template<
bool reset,
class G,
class F,
class VM>
660 VM& visitedMap)
const;
685 template<
bool remove,
bool reset,
class G,
class L,
class F1,
class F2,
class VM>
688 const G& graph, L& visited, F1& aggregateVisitor,
689 F2& nonAggregateVisitor,
690 VM& visitedMap)
const;
760 std::size_t noVertices_;
766 template<
class G,
class C>
768 const typename C::Matrix& matrix,
776 template<
class G,
class S>
820 VertexSet& connectivity, std::vector<Vertex>& front_);
845 void add(std::vector<Vertex>& vertex);
854 typename VertexSet::size_type
size();
901 std::vector<Vertex>& front_;
951 template<
class M,
class C>
952 std::tuple<int,int,int,int>
build(
const M& m, G& graph,
960 typedef PoolAllocator<Vertex,100> Allocator;
965 typedef SLList<Vertex,Allocator> VertexList;
970 typedef std::set<Vertex,std::less<Vertex>,Allocator> VertexSet;
975 typedef std::size_t* SphereMap;
990 std::vector<Vertex> front_;
995 VertexSet connected_;
1016 enum { N = 1300000 };
1058 class AggregateVisitor
1118 class FrontNeighbourCounter :
public Counter
1137 int noFrontNeighbours(
const Vertex& vertex)
const;
1142 class TwoWayCounter :
public Counter
1165 class OneWayCounter :
public Counter
1191 class ConnectivityCounter :
public Counter
1206 const VertexSet& connected_;
1241 bool connected(
const Vertex& vertex,
const SLList<AggregateDescriptor>& aggregateList,
1251 class DependencyCounter :
public Counter
1283 std::vector<Vertex>& front_;
1322 std::pair<int,int> neighbours(
const Vertex& vertex,
1379 void nonisoNeighbourAggregate(
const Vertex& vertex,
1381 SLList<Vertex>& neighbours)
const;
1398 template<
class M,
class N>
1404 template<
class M,
class N>
1407 initRow(row, index, std::is_convertible<field_type, real_type>());
1410 template<
class M,
class N>
1413 DUNE_THROW(InvalidStateException,
"field_type needs to convertible to real_type");
1416 template<
class M,
class N>
1420 maxValue_ = min(- std::numeric_limits<typename Matrix::field_type>::max(), std::numeric_limits<typename Matrix::field_type>::min());
1422 diagonal_ = norm_(matrix_->operator[](row_)[row_]);
1425 template<
class M,
class N>
1429 real_type eij = norm_(*
col);
1431 matrix_->operator[](
col.index()).find(row_);
1432 if ( opposite_entry == matrix_->operator[](
col.index()).end() )
1437 real_type eji = norm_(*opposite_entry);
1440 if(!N::is_sign_preserving || eij<0 || eji<0)
1441 maxValue_ = max(maxValue_,
1442 eij /diagonal_ * eji/
1443 norm_(matrix_->operator[](
col.index())[
col.index()]));
1446 template<
class M,
class N>
1450 real_type eij = norm_(*
col);
1452 matrix_->operator[](
col.index()).find(row_);
1454 if ( opposite_entry == matrix_->operator[](
col.index()).end() )
1459 real_type eji = norm_(*opposite_entry);
1461 if(!N::is_sign_preserving || (eij<0 || eji<0))
1462 if(eji / norm_(matrix_->operator[](edge.target())[edge.target()]) *
1463 eij/ diagonal_ > alpha() * maxValue_) {
1464 edge.properties().setDepends();
1465 edge.properties().setInfluences();
1466 typename G::EdgeProperties& other = graph.getEdgeProperties(edge.target(), edge.source());
1467 other.setInfluences();
1472 template<
class M,
class N>
1475 return maxValue_ < beta();
1479 template<
class M,
class N>
1485 template<
class M,
class N>
1489 maxValue_ = min(- std::numeric_limits<real_type>::max(), std::numeric_limits<real_type>::min());
1491 diagonal_ = norm_(matrix_->operator[](row_)[row_]);
1494 template<
class M,
class N>
1498 maxValue_ = max(maxValue_, -norm_(*
col));
1501 template<
class M,
class N>
1505 if(-norm_(*
col) >= maxValue_ * alpha()) {
1506 edge.properties().setDepends();
1507 typedef typename G::EdgeDescriptor ED;
1508 ED e= graph.findEdge(edge.target(), edge.source());
1509 if(e!=std::numeric_limits<ED>::max())
1511 typename G::EdgeProperties& other = graph.getEdgeProperties(e);
1512 other.setInfluences();
1517 template<
class M,
class N>
1520 return maxValue_ < beta() * diagonal_;
1523 template<
class G,
class S>
1525 VertexSet& connected, std::vector<Vertex>& front)
1526 : vertices_(), id_(-1), graph_(graph), aggregates_(aggregates),
1527 connected_(connected), front_(front)
1530 template<
class G,
class S>
1538 throw "Not yet implemented";
1546 template<
class G,
class S>
1549 dvverb<<
"Connected cleared"<<std::endl;
1552 connected_.insert(vertex);
1553 dvverb <<
" Inserting "<<vertex<<
" size="<<connected_.size();
1559 template<
class G,
class S>
1562 vertices_.insert(vertex);
1563 aggregates_[vertex]=id_;
1565 front_.erase(std::lower_bound(front_.begin(), front_.end(), vertex));
1569 const iterator end = graph_.endEdges(vertex);
1570 for(iterator edge = graph_.beginEdges(vertex); edge != end; ++edge) {
1571 dvverb <<
" Inserting "<<aggregates_[edge.target()];
1572 connected_.insert(aggregates_[edge.target()]);
1573 dvverb <<
" size="<<connected_.size();
1575 !graph_.getVertexProperties(edge.target()).front())
1577 front_.push_back(edge.target());
1578 graph_.getVertexProperties(edge.target()).setFront();
1582 std::sort(front_.begin(), front_.end());
1585 template<
class G,
class S>
1589 std::size_t oldsize = vertices_.size();
1591 typedef typename std::vector<Vertex>::iterator Iterator;
1593 typedef typename VertexSet::iterator SIterator;
1595 SIterator pos=vertices_.begin();
1596 std::vector<Vertex> newFront;
1597 newFront.reserve(front_.capacity());
1599 std::set_difference(front_.begin(), front_.end(), vertices.begin(), vertices.end(),
1600 std::back_inserter(newFront));
1603 for(Iterator vertex=vertices.begin(); vertex != vertices.end(); ++vertex)
1605 pos=vertices_.insert(pos,*vertex);
1606 vertices_.insert(*vertex);
1607 graph_.getVertexProperties(*vertex).resetFront();
1608 aggregates_[*vertex]=id_;
1611 const iterator end = graph_.endEdges(*vertex);
1612 for(iterator edge = graph_.beginEdges(*vertex); edge != end; ++edge) {
1613 dvverb <<
" Inserting "<<aggregates_[edge.target()];
1614 connected_.insert(aggregates_[edge.target()]);
1616 !graph_.getVertexProperties(edge.target()).front())
1618 front_.push_back(edge.target());
1619 graph_.getVertexProperties(edge.target()).setFront();
1621 dvverb <<
" size="<<connected_.size();
1625 std::sort(front_.begin(), front_.end());
1626 assert(oldsize+vertices.size()==vertices_.size());
1628 template<
class G,
class S>
1636 template<
class G,
class S>
1637 inline typename Aggregate<G,S>::VertexSet::size_type
1640 return vertices_.size();
1643 template<
class G,
class S>
1644 inline typename Aggregate<G,S>::VertexSet::size_type
1647 return connected_.size();
1650 template<
class G,
class S>
1656 template<
class G,
class S>
1659 return vertices_.begin();
1662 template<
class G,
class S>
1665 return vertices_.end();
1683 delete[] aggregates_;
1690 allocate(noVertices);
1703 noVertices_ = noVertices;
1705 for(std::size_t i=0; i < noVertices; i++)
1706 aggregates_[i]=UNAGGREGATED;
1712 assert(aggregates_ != 0);
1713 delete[] aggregates_;
1721 return aggregates_[v];
1728 return aggregates_[v];
1732 template<
bool reset,
class G,
class F,
class VM>
1735 const G& graph, F& aggregateVisitor,
1736 VM& visitedMap)
const
1740 DummyEdgeVisitor dummy;
1741 return breadthFirstSearch<true,reset>(start, aggregate, graph, vlist, aggregateVisitor, dummy, visitedMap);
1745 template<
bool remove,
bool reset,
class G,
class L,
class F1,
class F2,
class VM>
1750 F1& aggregateVisitor,
1751 F2& nonAggregateVisitor,
1752 VM& visitedMap)
const
1754 typedef typename L::const_iterator ListIterator;
1755 int visitedSpheres = 0;
1757 visited.push_back(start);
1758 put(visitedMap, start,
true);
1760 ListIterator current = visited.begin();
1761 ListIterator end = visited.end();
1762 std::size_t i=0, size=visited.size();
1766 while(current != end) {
1768 for(; i<size; ++current, ++i) {
1769 typedef typename G::ConstEdgeIterator EdgeIterator;
1770 const EdgeIterator endEdge = graph.endEdges(*current);
1772 for(EdgeIterator edge = graph.beginEdges(*current);
1773 edge != endEdge; ++edge) {
1775 if(aggregates_[edge.target()]==aggregate) {
1776 if(!
get(visitedMap, edge.target())) {
1777 put(visitedMap, edge.target(),
true);
1778 visited.push_back(edge.target());
1779 aggregateVisitor(edge);
1782 nonAggregateVisitor(edge);
1785 end = visited.end();
1786 size = visited.size();
1792 for(current = visited.begin(); current != end; ++current)
1793 put(visitedMap, *current,
false);
1799 return visitedSpheres;
1804 : graph_(0), aggregate_(0), front_(), connected_(), size_(-1)
1813 template<
class G,
class C>
1815 const typename C::Matrix& matrix,
1816 C criterion,
bool firstlevel)
1819 typedef typename C::Matrix Matrix;
1820 typedef typename G::VertexIterator VertexIterator;
1822 criterion.init(&matrix);
1824 for(VertexIterator vertex = graph.begin(); vertex != graph.end(); ++vertex) {
1827 const Row& row = matrix[*vertex];
1832 criterion.initRow(row, *vertex);
1837 ColIterator end = row.end();
1838 typename FieldTraits<typename Matrix::field_type>::real_type absoffdiag=0.;
1842 for(ColIterator
col = row.begin();
col != end; ++
col)
1843 if(
col.index()!=*vertex) {
1844 criterion.examine(
col);
1845 absoffdiag = max(absoffdiag, Impl::asMatrix(*col).frobenius_norm());
1849 vertex.properties().setExcludedBorder();
1852 for(ColIterator
col = row.begin();
col != end; ++
col)
1853 if(
col.index()!=*vertex)
1854 criterion.examine(
col);
1860 if(criterion.isIsolated()) {
1862 vertex.properties().setIsolated();
1865 auto eEnd = vertex.end();
1866 auto col = matrix[*vertex].begin();
1868 for(
auto edge = vertex.begin(); edge!= eEnd; ++edge, ++
col) {
1870 while(
col.index()!=edge.target())
1872 criterion.examine(graph, edge,
col);
1882 inline Aggregator<G>::AggregateVisitor<V>::AggregateVisitor(
const AggregatesMap<Vertex>& aggregates,
1884 : aggregates_(aggregates), aggregate_(aggregate), visitor_(&visitor)
1891 if(aggregates_[edge.target()]==aggregate_)
1892 visitor_->operator()(edge);
1897 inline void Aggregator<G>::visitAggregateNeighbours(
const Vertex& vertex,
1899 const AggregatesMap<Vertex>& aggregates,
1903 AggregateVisitor<V> v(aggregates, aggregate, visitor);
1909 inline Aggregator<G>::Counter::Counter()
1914 inline void Aggregator<G>::Counter::increment()
1920 inline void Aggregator<G>::Counter::decrement()
1925 inline int Aggregator<G>::Counter::value()
1933 if(edge.properties().isTwoWay())
1934 Counter::increment();
1939 const AggregatesMap<Vertex>& aggregates)
const
1941 TwoWayCounter counter;
1942 visitAggregateNeighbours(vertex, aggregate, aggregates, counter);
1943 return counter.value();
1948 const AggregatesMap<Vertex>& aggregates)
const
1950 OneWayCounter counter;
1951 visitAggregateNeighbours(vertex, aggregate, aggregates, counter);
1952 return counter.value();
1958 if(edge.properties().isOneWay())
1959 Counter::increment();
1963 inline Aggregator<G>::ConnectivityCounter::ConnectivityCounter(
const VertexSet& connected,
1964 const AggregatesMap<Vertex>& aggregates)
1965 : Counter(), connected_(connected), aggregates_(aggregates)
1974 Counter::increment();
1976 Counter::increment();
1977 Counter::increment();
1982 inline double Aggregator<G>::connectivity(
const Vertex& vertex,
const AggregatesMap<Vertex>& aggregates)
const
1984 ConnectivityCounter counter(connected_, aggregates);
1986 return (
double)counter.value()/noNeighbours;
1990 inline Aggregator<G>::DependencyCounter::DependencyCounter()
1997 if(edge.properties().depends())
1998 Counter::increment();
1999 if(edge.properties().influences())
2000 Counter::increment();
2004 int Aggregator<G>::unusedNeighbours(
const Vertex& vertex,
const AggregatesMap<Vertex>& aggregates)
const
2010 std::pair<int,int> Aggregator<G>::neighbours(
const Vertex& vertex,
2012 const AggregatesMap<Vertex>& aggregates)
const
2014 DependencyCounter unused, aggregated;
2015 typedef AggregateVisitor<DependencyCounter> CounterT;
2016 typedef std::tuple<CounterT,CounterT> CounterTuple;
2019 return std::make_pair(unused.value(), aggregated.value());
2024 int Aggregator<G>::aggregateNeighbours(
const Vertex& vertex,
const AggregateDescriptor& aggregate,
const AggregatesMap<Vertex>& aggregates)
const
2026 DependencyCounter counter;
2027 visitAggregateNeighbours(vertex, aggregate, aggregates, counter);
2028 return counter.value();
2032 std::size_t Aggregator<G>::distance(
const Vertex& vertex,
const AggregatesMap<Vertex>& aggregates)
2035 typename PropertyMapTypeSelector<VertexVisitedTag,G>::Type visitedMap =
get(VertexVisitedTag(), *graph_);
2037 typename AggregatesMap<Vertex>::DummyEdgeVisitor dummy;
2038 return aggregates.template breadthFirstSearch<true,true>(vertex,
2039 aggregate_->
id(), *graph_,
2040 vlist, dummy, dummy, visitedMap);
2044 inline Aggregator<G>::FrontMarker::FrontMarker(std::vector<Vertex>& front,
MatrixGraph& graph)
2045 : front_(front), graph_(graph)
2051 Vertex target = edge.target();
2053 if(!graph_.getVertexProperties(target).front()) {
2054 front_.push_back(target);
2055 graph_.getVertexProperties(target).setFront();
2060 inline bool Aggregator<G>::admissible(
const Vertex& vertex,
const AggregateDescriptor& aggregate,
const AggregatesMap<Vertex>& aggregates)
const
2063 Dune::dvverb<<
" Admissible not yet implemented!"<<std::endl;
2070 Iterator vend = graph_->endEdges(vertex);
2071 for(Iterator edge = graph_->beginEdges(vertex); edge != vend; ++edge) {
2073 if(edge.properties().isStrong()
2074 && aggregates[edge.target()]==aggregate)
2077 Iterator edge1 = edge;
2078 for(++edge1; edge1 != vend; ++edge1) {
2080 if(edge1.properties().isStrong()
2081 && aggregates[edge.target()]==aggregate)
2086 Iterator v2end = graph_->endEdges(edge.target());
2087 for(Iterator edge2 = graph_->beginEdges(edge.target()); edge2 != v2end; ++edge2) {
2088 if(edge2.target()==edge1.target() &&
2089 edge2.properties().isStrong()) {
2105 vend = graph_->endEdges(vertex);
2106 for(Iterator edge = graph_->beginEdges(vertex); edge != vend; ++edge) {
2108 if(edge.properties().isStrong()
2109 && aggregates[edge.target()]==aggregate)
2112 Iterator v1end = graph_->endEdges(edge.target());
2114 for(Iterator edge1=graph_->beginEdges(edge.target()); edge1 != v1end; ++edge1) {
2116 if(edge1.properties().isStrong()
2117 && aggregates[edge1.target()]==aggregate)
2121 Iterator v2end = graph_->endEdges(vertex);
2122 for(Iterator edge2 = graph_->beginEdges(vertex); edge2 != v2end; ++edge2) {
2123 if(edge2.target()==edge1.target()) {
2124 if(edge2.properties().isStrong())
2141 void Aggregator<G>::unmarkFront()
2143 typedef typename std::vector<Vertex>::const_iterator Iterator;
2145 for(Iterator vertex=front_.begin(); vertex != front_.end(); ++vertex)
2146 graph_->getVertexProperties(*vertex).resetFront();
2153 Aggregator<G>::nonisoNeighbourAggregate(
const Vertex& vertex,
2154 const AggregatesMap<Vertex>& aggregates,
2155 SLList<Vertex>& neighbours)
const
2158 Iterator end=graph_->beginEdges(vertex);
2161 for(Iterator edge=graph_->beginEdges(vertex); edge!=end; ++edge)
2164 neighbours.push_back(aggregates[edge.target()]);
2169 inline typename G::VertexDescriptor Aggregator<G>::mergeNeighbour(
const Vertex& vertex,
const AggregatesMap<Vertex>& aggregates)
const
2173 Iterator end = graph_->endEdges(vertex);
2174 for(Iterator edge = graph_->beginEdges(vertex); edge != end; ++edge) {
2176 graph_->getVertexProperties(edge.target()).isolated() == graph_->getVertexProperties(edge.source()).isolated()) {
2177 if( graph_->getVertexProperties(vertex).isolated() ||
2178 ((edge.properties().depends() || edge.properties().influences())
2179 && admissible(vertex, aggregates[edge.target()], aggregates)))
2180 return edge.target();
2187 Aggregator<G>::FrontNeighbourCounter::FrontNeighbourCounter(
const MatrixGraph& graph)
2188 : Counter(), graph_(graph)
2194 if(graph_.getVertexProperties(edge.target()).front())
2195 Counter::increment();
2199 int Aggregator<G>::noFrontNeighbours(
const Vertex& vertex)
const
2201 FrontNeighbourCounter counter(*graph_);
2203 return counter.value();
2206 inline bool Aggregator<G>::connected(
const Vertex& vertex,
2208 const AggregatesMap<Vertex>& aggregates)
const
2210 typedef typename G::ConstEdgeIterator iterator;
2211 const iterator end = graph_->endEdges(vertex);
2212 for(iterator edge = graph_->beginEdges(vertex); edge != end; ++edge)
2213 if(aggregates[edge.target()]==aggregate)
2218 inline bool Aggregator<G>::connected(
const Vertex& vertex,
2219 const SLList<AggregateDescriptor>& aggregateList,
2220 const AggregatesMap<Vertex>& aggregates)
const
2222 typedef typename SLList<AggregateDescriptor>::const_iterator Iter;
2223 for(Iter i=aggregateList.begin(); i!=aggregateList.end(); ++i)
2224 if(connected(vertex, *i, aggregates))
2231 void Aggregator<G>::growIsolatedAggregate(
const Vertex& seed,
const AggregatesMap<Vertex>& aggregates,
const C& c)
2233 SLList<Vertex> connectedAggregates;
2234 nonisoNeighbourAggregate(seed, aggregates,connectedAggregates);
2236 while(aggregate_->
size()< c.minAggregateSize() && aggregate_->
connectSize() < c.maxConnectivity()) {
2238 std::size_t maxFrontNeighbours=0;
2242 typedef typename std::vector<Vertex>::const_iterator Iterator;
2244 for(Iterator vertex = front_.begin(); vertex != front_.end(); ++vertex) {
2245 if(distance(*vertex, aggregates)>c.maxDistance())
2248 if(connectedAggregates.size()>0) {
2252 if(!connected(*vertex, connectedAggregates, aggregates))
2256 double con = connectivity(*vertex, aggregates);
2259 std::size_t frontNeighbours = noFrontNeighbours(*vertex);
2261 if(frontNeighbours >= maxFrontNeighbours) {
2262 maxFrontNeighbours = frontNeighbours;
2263 candidate = *vertex;
2265 }
else if(con > maxCon) {
2267 maxFrontNeighbours = noFrontNeighbours(*vertex);
2268 candidate = *vertex;
2275 aggregate_->
add(candidate);
2281 void Aggregator<G>::growAggregate(
const Vertex& seed,
const AggregatesMap<Vertex>& aggregates,
const C& c)
2285 std::size_t distance_ =0;
2286 while(aggregate_->
size() < c.minAggregateSize()&& distance_<c.maxDistance()) {
2287 int maxTwoCons=0, maxOneCons=0, maxNeighbours=-1;
2290 std::vector<Vertex> candidates;
2291 candidates.reserve(30);
2293 typedef typename std::vector<Vertex>::const_iterator Iterator;
2295 for(Iterator vertex = front_.begin(); vertex != front_.end(); ++vertex) {
2297 if(graph_->getVertexProperties(*vertex).isolated())
2300 int twoWayCons = twoWayConnections(*vertex, aggregate_->
id(), aggregates);
2303 if( maxTwoCons == twoWayCons && twoWayCons > 0) {
2304 double con = connectivity(*vertex, aggregates);
2307 int neighbours = noFrontNeighbours(*vertex);
2309 if(neighbours > maxNeighbours) {
2310 maxNeighbours = neighbours;
2312 candidates.push_back(*vertex);
2314 candidates.push_back(*vertex);
2316 }
else if( con > maxCon) {
2318 maxNeighbours = noFrontNeighbours(*vertex);
2320 candidates.push_back(*vertex);
2322 }
else if(twoWayCons > maxTwoCons) {
2323 maxTwoCons = twoWayCons;
2324 maxCon = connectivity(*vertex, aggregates);
2325 maxNeighbours = noFrontNeighbours(*vertex);
2327 candidates.push_back(*vertex);
2330 maxOneCons = std::numeric_limits<int>::max();
2339 int oneWayCons = oneWayConnections(*vertex, aggregate_->
id(), aggregates);
2344 if(!admissible(*vertex, aggregate_->
id(), aggregates))
2347 if( maxOneCons == oneWayCons && oneWayCons > 0) {
2348 double con = connectivity(*vertex, aggregates);
2351 int neighbours = noFrontNeighbours(*vertex);
2353 if(neighbours > maxNeighbours) {
2354 maxNeighbours = neighbours;
2356 candidates.push_back(*vertex);
2358 if(neighbours==maxNeighbours)
2360 candidates.push_back(*vertex);
2363 }
else if( con > maxCon) {
2365 maxNeighbours = noFrontNeighbours(*vertex);
2367 candidates.push_back(*vertex);
2369 }
else if(oneWayCons > maxOneCons) {
2370 maxOneCons = oneWayCons;
2371 maxCon = connectivity(*vertex, aggregates);
2372 maxNeighbours = noFrontNeighbours(*vertex);
2374 candidates.push_back(*vertex);
2379 if(!candidates.size())
2381 distance_=distance(seed, aggregates);
2382 candidates.resize(min(candidates.size(), c.maxAggregateSize()-
2383 aggregate_->
size()));
2384 aggregate_->
add(candidates);
2388 template<
typename V>
2389 template<
typename M,
typename G,
typename C>
2393 Aggregator<G> aggregator;
2394 return aggregator.build(matrix, graph, *
this, criterion, finestLevel);
2398 template<
class M,
class C>
2399 std::tuple<int,int,int,int>
Aggregator<G>::build(
const M& m, G& graph, AggregatesMap<Vertex>& aggregates,
const C& c,
2405 Stack stack_(graph, *
this, aggregates);
2409 aggregate_ =
new Aggregate<G,VertexSet>(graph, aggregates, connected_, front_);
2416 dverb<<
"Build dependency took "<< watch.elapsed()<<
" seconds."<<std::endl;
2417 int noAggregates, conAggregates, isoAggregates, oneAggregates;
2418 std::size_t maxA=0, minA=1000000, avg=0;
2419 int skippedAggregates;
2420 noAggregates = conAggregates = isoAggregates = oneAggregates =
2421 skippedAggregates = 0;
2424 Vertex seed = stack_.pop();
2426 if(seed == Stack::NullEntry)
2431 if((noAggregates+1)%10000 == 0)
2435 if(graph.getVertexProperties(seed).excludedBorder()) {
2437 ++skippedAggregates;
2441 if(graph.getVertexProperties(seed).isolated()) {
2442 if(c.skipIsolated()) {
2445 ++skippedAggregates;
2449 aggregate_->
seed(seed);
2450 growIsolatedAggregate(seed, aggregates, c);
2453 aggregate_->
seed(seed);
2454 growAggregate(seed, aggregates, c);
2458 while(!(graph.getVertexProperties(seed).isolated()) && aggregate_->
size() < c.maxAggregateSize()) {
2460 std::vector<Vertex> candidates;
2461 candidates.reserve(30);
2463 typedef typename std::vector<Vertex>::const_iterator Iterator;
2465 for(Iterator vertex = front_.begin(); vertex != front_.end(); ++vertex) {
2467 if(graph.getVertexProperties(*vertex).isolated())
2470 if(twoWayConnections( *vertex, aggregate_->
id(), aggregates) == 0 &&
2471 (oneWayConnections( *vertex, aggregate_->
id(), aggregates) == 0 ||
2472 !admissible( *vertex, aggregate_->
id(), aggregates) ))
2475 std::pair<int,int> neighbourPair=neighbours(*vertex, aggregate_->
id(),
2481 if(neighbourPair.first >= neighbourPair.second)
2484 if(distance(*vertex, aggregates) > c.maxDistance())
2486 candidates.push_back(*vertex);
2490 if(!candidates.size())
break;
2492 candidates.resize(min(candidates.size(), c.maxAggregateSize()-
2493 aggregate_->
size()));
2494 aggregate_->
add(candidates);
2499 if(aggregate_->
size()==1 && c.maxAggregateSize()>1) {
2500 if(!graph.getVertexProperties(seed).isolated()) {
2501 Vertex mergedNeighbour = mergeNeighbour(seed, aggregates);
2505 aggregates[seed] = aggregates[mergedNeighbour];
2509 minA=min(minA,
static_cast<std::size_t
>(1));
2510 maxA=max(maxA,
static_cast<std::size_t
>(1));
2516 minA=min(minA,
static_cast<std::size_t
>(1));
2517 maxA=max(maxA,
static_cast<std::size_t
>(1));
2523 avg+=aggregate_->
size();
2524 minA=min(minA,aggregate_->
size());
2525 maxA=max(maxA,aggregate_->
size());
2526 if(graph.getVertexProperties(seed).isolated())
2534 Dune::dinfo<<
"connected aggregates: "<<conAggregates;
2535 Dune::dinfo<<
" isolated aggregates: "<<isoAggregates;
2536 if(conAggregates+isoAggregates>0)
2537 Dune::dinfo<<
" one node aggregates: "<<oneAggregates<<
" min size="
2538 <<minA<<
" max size="<<maxA
2539 <<
" avg="<<avg/(conAggregates+isoAggregates)<<std::endl;
2542 return std::make_tuple(conAggregates+isoAggregates,isoAggregates,
2543 oneAggregates,skippedAggregates);
2548 Aggregator<G>::Stack::Stack(
const MatrixGraph& graph,
const Aggregator<G>& aggregatesBuilder,
2549 const AggregatesMap<Vertex>& aggregates)
2550 : graph_(graph), aggregatesBuilder_(aggregatesBuilder), aggregates_(aggregates), begin_(graph.begin()), end_(graph.end())
2556 Aggregator<G>::Stack::~Stack()
2564 = std::numeric_limits<typename G::VertexDescriptor>::max();
2567 inline typename G::VertexDescriptor Aggregator<G>::Stack::pop()
2573 typename G::VertexDescriptor current=*begin_;
2587 std::ios_base::fmtflags oldOpts=os.flags();
2589 os.setf(std::ios_base::right, std::ios_base::adjustfield);
2594 for(
int i=0; i< n*m; i++)
2595 maxVal=max(maxVal, aggregates[i]);
2597 for(
int i=10; i < 1000000; i*=10)
2603 for(
int j=0, entry=0; j < m; j++) {
2604 for(
int i=0; i<n; i++, entry++) {
2606 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:351
Matrix::ConstColIterator ColIter
Constant column iterator of the matrix.
Definition: aggregates.hh:273
std::vector< real_type >::iterator valIter_
Definition: aggregates.hh:189
Matrix::ConstColIterator ColIter
Constant column iterator of the matrix.
Definition: aggregates.hh:154
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:586
iterator begin()
Definition: aggregates.hh:737
int id()
Get the id identifying the aggregate.
Norm norm_
The functor for calculating the norm.
Definition: aggregates.hh:302
void operator()([[maybe_unused]] const EdgeIterator &edge) const
Definition: aggregates.hh:601
MatrixGraph::VertexDescriptor Vertex
The vertex identifier.
Definition: aggregates.hh:920
AggregationCriterion()
Constructor.
Definition: aggregates.hh:66
const Matrix * matrix_
The matrix we work on.
Definition: aggregates.hh:357
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:406
void initRow(const Row &row, int index)
SymmetricMatrixDependency(const Parameters &parms)
Definition: aggregates.hh:168
M Matrix
The matrix type we build the dependency of.
Definition: aggregates.hh:258
G MatrixGraph
The matrix graph type used.
Definition: aggregates.hh:915
FieldTraits< typename M::field_type >::real_type operator()(const M &m) const
compute the norm of a matrix.
Definition: aggregates.hh:490
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:363
SymmetricCriterion()
Definition: aggregates.hh:524
Dependency()
Definition: aggregates.hh:290
void examine(const ColIter &col)
M Matrix
The matrix type we build the dependency of.
Definition: aggregates.hh:319
real_type diagonal_
The norm of the current diagonal.
Definition: aggregates.hh:187
N Norm
The norm to use for examining the matrix entries.
Definition: aggregates.hh:263
iterator end()
Definition: aggregates.hh:742
UnSymmetricCriterion(const Parameters &parms)
Definition: aggregates.hh:541
void operator()(const typename MatrixGraph::ConstEdgeIterator &edge)
void initRow(const Row &row, int index)
Definition: aggregates.hh:201
static const Vertex NullEntry
Definition: aggregates.hh:1008
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:113
void examine(const ColIter &col)
Definition: aggregates.hh:214
Dependency(const Parameters &parms)
Definition: aggregates.hh:286
int row_
index of the currently evaluated row.
Definition: aggregates.hh:185
friend class Stack
Definition: aggregates.hh:1036
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:329
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:390
const Matrix * matrix_
The matrix we work on.
Definition: aggregates.hh:296
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:334
~AggregatesMap()
Destructor.
Matrix::field_type field_type
The current max value.
Definition: aggregates.hh:179
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:1064
std::size_t * SphereMap
Type of the mapping of aggregate members onto distance spheres.
Definition: aggregates.hh:809
AggregationCriterion(const Parameters &parms)
Definition: aggregates.hh:70
Matrix::row_type Row
Constant Row iterator of the matrix.
Definition: aggregates.hh:268
VertexSet::size_type connectSize()
Get tne number of connections to other aggregates.
std::vector< real_type > vals_
Definition: aggregates.hh:188
FieldTraits< typename M::field_type >::real_type operator()(const M &m) const
compute the norm of a matrix.
Definition: aggregates.hh:506
N Norm
The norm to use for examining the matrix entries.
Definition: aggregates.hh:324
void printAggregates2d(const AggregatesMap< V > &aggregates, int n, int m, std::ostream &os)
Definition: aggregates.hh:2583
void invalidate()
Definition: aggregates.hh:822
const_iterator begin() const
Definition: aggregates.hh:725
real_type maxValue_
Definition: aggregates.hh:181
Norm norm_
The functor for calculating the norm.
Definition: aggregates.hh:183
void operator()(const typename MatrixGraph::ConstEdgeIterator &edge)
VertexSet::const_iterator const_iterator
Const iterator over a vertex list.
Definition: aggregates.hh:804
SymmetricCriterion(const Parameters &parms)
Definition: aggregates.hh:521
MatrixGraph::VertexDescriptor AggregateDescriptor
The type of the aggregate descriptor.
Definition: aggregates.hh:923
real_type maxValue_
Definition: aggregates.hh:361
void init(const Matrix *matrix)
Definition: aggregates.hh:195
real_type diagonal_
The norm of the current diagonal.
Definition: aggregates.hh:367
AggregateDescriptor * iterator
Definition: aggregates.hh:735
SymmetricDependency(const Parameters &parms)
Definition: aggregates.hh:348
FieldTraits< field_type >::real_type real_type
Definition: aggregates.hh:360
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:55
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:82
const_iterator end() const
Definition: aggregates.hh:730
V AggregateDescriptor
The aggregate descriptor type.
Definition: aggregates.hh:580
static const V ISOLATED
Identifier of isolated vertices.
Definition: aggregates.hh:571
int row_
index of the currently evaluated row.
Definition: aggregates.hh:365
SymmetricMatrixDependency()
Definition: aggregates.hh:171
DependencyCounter()
Constructor.
real_type diagonal_
The norm of the current diagonal.
Definition: aggregates.hh:306
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:105
AggregatesMap(std::size_t noVertices)
Constructs with allocating memory.
Matrix::field_type field_type
The current max value.
Definition: aggregates.hh:359
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:592
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:544
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:473
const AggregateDescriptor * const_iterator
Definition: aggregates.hh:723
int row_
index of the currently evaluated row.
Definition: aggregates.hh:304
Stack(const MatrixGraph &graph, const Aggregator< G > &aggregatesBuilder, const AggregatesMap< Vertex > &aggregates)
FieldTraits< field_type >::real_type real_type
Definition: aggregates.hh:180
void initRow(const Row &row, int index)
M Matrix
The matrix type we build the dependency of.
Definition: aggregates.hh:139
FieldTraits< field_type >::real_type real_type
Definition: aggregates.hh:299
const Matrix * matrix_
The matrix we work on.
Definition: aggregates.hh:177
S VertexSet
The type of a single linked list of vertex descriptors.
Definition: aggregates.hh:801
static const V UNAGGREGATED
Identifier of not yet aggregated vertices.
Definition: aggregates.hh:566
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:298
void operator()(const typename MatrixGraph::ConstEdgeIterator &edge)
bool isIsolated()
Definition: aggregates.hh:240
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:144
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:789
void seed(const Vertex &vertex)
Initialize the aggregate with one vertex.
SymmetricDependency()
Definition: aggregates.hh:351
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:575
real_type maxValue_
Definition: aggregates.hh:300
Matrix::row_type Row
Constant Row iterator of the matrix.
Definition: aggregates.hh:149
PoolAllocator< Vertex, 100 > Allocator
The allocator we use for our lists and the set.
Definition: aggregates.hh:795
G MatrixGraph
Definition: aggregates.hh:785
void operator()(const typename MatrixGraph::ConstEdgeIterator &edge)
@ is_sign_preserving
Definition: aggregates.hh:483
@ is_sign_preserving
Definition: aggregates.hh:382
@ is_sign_preserving
Definition: aggregates.hh:499
@ is_sign_preserving
Definition: aggregates.hh:466
Definition: allocator.hh:11
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:293
derive error class from the base class in common
Definition: istlexception.hh:19
A generic dynamic dense matrix.
Definition: matrix.hh:561
typename Imp::BlockTraits< T >::field_type field_type
Export the type representing the underlying field.
Definition: matrix.hh:565
row_type::const_iterator ConstColIterator
Const iterator for the entries of each row.
Definition: matrix.hh:589
MatrixImp::DenseMatrixBase< T, A >::window_type row_type
The type implementing a matrix row.
Definition: matrix.hh:574
Base class of all aggregation criterions.
Definition: aggregates.hh:49
Dependency policy for symmetric matrices.
Definition: aggregates.hh:134
Dependency policy for symmetric matrices.
Definition: aggregates.hh:253
Dependency policy for symmetric matrices.
Definition: aggregates.hh:314
Norm that uses only the [N][N] entry of the block to determine couplings.
Definition: aggregates.hh:379
Norm that uses only the [0][0] entry of the block to determine couplings.
Definition: aggregates.hh:455
Functor using the row sum (infinity) norm to determine strong couplings.
Definition: aggregates.hh:463
Definition: aggregates.hh:480
Definition: aggregates.hh:496
Criterion taking advantage of symmetric matrices.
Definition: aggregates.hh:519
Criterion suitable for unsymmetric matrices.
Definition: aggregates.hh:539
Class for building the aggregates.
Definition: aggregates.hh:909
Class providing information about the mapping of the vertices onto aggregates.
Definition: aggregates.hh:560
A Dummy visitor that does nothing for each visited edge.
Definition: aggregates.hh:598
A class for temporarily storing the vertices of an aggregate in.
Definition: aggregates.hh:778
M::size_type VertexDescriptor
The vertex descriptor.
Definition: graph.hh:73
EdgeIteratorT< const MatrixGraph< Matrix > > ConstEdgeIterator
The constant edge iterator type.
Definition: graph.hh:298
Iterator over all edges starting from a vertex.
Definition: graph.hh:95
The vertex iterator type of the graph.
Definition: graph.hh:209
All parameters for AMG.
Definition: parameters.hh:393