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/unused.hh>
17 #include <dune/common/ftraits.hh>
18 #include <dune/common/scalarmatrixview.hh>
83 this->setMaxDistance(diameter-1);
88 this->setMaxDistance(this->maxDistance()+diameter-1);
90 this->setMinAggregateSize(csize);
91 this->setMaxAggregateSize(
static_cast<std::size_t
>(csize*1.5));
107 this->setMaxDistance(this->maxDistance()+dim-1);
114 os<<
"{ maxdistance="<<criterion.maxDistance()<<
" minAggregateSize="
115 <<criterion.minAggregateSize()<<
" maxAggregateSize="<<criterion.maxAggregateSize()
116 <<
" connectivity="<<criterion.maxConnectivity()<<
" debugLevel="<<criterion.debugLevel()<<
"}";
131 template<
class M,
class N>
162 void examine(G& graph,
const typename G::EdgeIterator& edge,
const ColIter&
col);
179 typedef typename FieldTraits<field_type>::real_type
real_type;
188 typename std::vector<real_type>::iterator
valIter_;
193 template<
class M,
class N>
199 template<
class M,
class N>
203 vals_.assign(row.size(), 0.0);
204 assert(vals_.size()==row.size());
205 valIter_=vals_.begin();
207 maxValue_ = min(- std::numeric_limits<real_type>::max(), std::numeric_limits<real_type>::min());
208 diagonal_=norm_(row[index]);
212 template<
class M,
class N>
218 if(!N::is_sign_preserving || eij<0)
220 *valIter_ = eij/diagonal_*eij/norm_(matrix_->operator[](
col.index())[
col.index()]);
221 maxValue_ = max(maxValue_, *valIter_);
227 template<
class M,
class N>
231 if(*valIter_ > alpha() * maxValue_) {
232 edge.properties().setDepends();
233 edge.properties().setInfluences();
238 template<
class M,
class N>
243 valIter_=vals_.begin();
244 return maxValue_ < beta();
250 template<
class M,
class N>
281 void examine(G& graph,
const typename G::EdgeIterator& edge,
const ColIter&
col);
298 typedef typename FieldTraits<field_type>::real_type
real_type;
311 template<
class M,
class N>
342 void examine(G& graph,
const typename G::EdgeIterator& edge,
const ColIter&
col);
359 typedef typename FieldTraits<field_type>::real_type
real_type;
386 typename FieldTraits<typename M::field_type>::real_type
operator()(
const M& m,
387 typename std::enable_if_t<!Dune::IsNumber<M>::value>* sfinae =
nullptr)
const
389 typedef typename M::field_type field_type;
390 typedef typename FieldTraits<field_type>::real_type real_type;
391 static_assert( std::is_convertible<field_type, real_type >::value,
392 "use of diagonal norm in AMG not implemented for complex field_type");
403 typename std::enable_if_t<Dune::IsNumber<M>::value>* sfinae =
nullptr)
const
405 typedef typename FieldTraits<M>::real_type real_type;
406 static_assert( std::is_convertible<M, real_type >::value,
407 "use of diagonal norm in AMG not implemented for complex field_type");
416 static T signed_abs(
const T & v)
423 static T signed_abs(
const std::complex<T> & v)
427 return csgn(v) * std::abs(v);
432 static T csgn(
const T & v)
434 return (T(0) < v) - (v < T(0));
439 static T csgn(std::complex<T> a)
441 return csgn(a.real())+(a.real() == 0.0)*csgn(a.imag());
469 typename FieldTraits<typename M::field_type>::real_type
operator()(
const M& m)
const
471 return m.infinity_norm();
486 typename FieldTraits<typename M::field_type>::real_type
operator()(
const M& m)
const
488 return m.frobenius_norm();
502 typename FieldTraits<typename M::field_type>::real_type
operator()(
const M& m)
const
513 template<
class M,
class Norm>
533 template<
class M,
class Norm>
596 template<
class EdgeIterator>
599 DUNE_UNUSED_PARAMETER(edge);
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
1158 const AggregatesMap<Vertex>& aggregates)
const;
1163 class OneWayCounter :
public Counter
1181 const AggregatesMap<Vertex>& aggregates)
const;
1189 class ConnectivityCounter :
public Counter
1198 ConnectivityCounter(
const VertexSet& connected,
const AggregatesMap<Vertex>& aggregates);
1204 const VertexSet& connected_;
1206 const AggregatesMap<Vertex>& aggregates_;
1221 double connectivity(
const Vertex& vertex,
const AggregatesMap<Vertex>& aggregates)
const;
1230 const AggregatesMap<Vertex>& aggregates)
const;
1239 bool connected(
const Vertex& vertex,
const SLList<AggregateDescriptor>& aggregateList,
1240 const AggregatesMap<Vertex>& aggregates)
const;
1249 class DependencyCounter :
public Counter
1255 DependencyCounter();
1275 FrontMarker(std::vector<Vertex>& front,
MatrixGraph& graph);
1281 std::vector<Vertex>& front_;
1305 int unusedNeighbours(
const Vertex& vertex,
const AggregatesMap<Vertex>& aggregates)
const;
1320 std::pair<int,int> neighbours(
const Vertex& vertex,
1322 const AggregatesMap<Vertex>& aggregates)
const;
1339 int aggregateNeighbours(
const Vertex& vertex,
const AggregateDescriptor& aggregate,
const AggregatesMap<Vertex>& aggregates)
const;
1348 bool admissible(
const Vertex& vertex,
const AggregateDescriptor& aggregate,
const AggregatesMap<Vertex>& aggregates)
const;
1357 std::size_t distance(
const Vertex& vertex,
const AggregatesMap<Vertex>& aggregates);
1367 Vertex mergeNeighbour(
const Vertex& vertex,
const AggregatesMap<Vertex>& aggregates)
const;
1377 void nonisoNeighbourAggregate(
const Vertex& vertex,
1378 const AggregatesMap<Vertex>& aggregates,
1379 SLList<Vertex>& neighbours)
const;
1389 void growAggregate(
const Vertex& vertex,
const AggregatesMap<Vertex>& aggregates,
const C& c);
1391 void growIsolatedAggregate(
const Vertex& vertex,
const AggregatesMap<Vertex>& aggregates,
const C& c);
1396 template<
class M,
class N>
1402 template<
class M,
class N>
1406 DUNE_UNUSED_PARAMETER(row);
1407 maxValue_ = min(- std::numeric_limits<typename Matrix::field_type>::max(), std::numeric_limits<typename Matrix::field_type>::min());
1409 diagonal_ = norm_(matrix_->operator[](row_)[row_]);
1412 template<
class M,
class N>
1416 real_type eij = norm_(*
col);
1418 matrix_->operator[](
col.index()).find(row_);
1419 if ( opposite_entry == matrix_->operator[](
col.index()).end() )
1424 real_type eji = norm_(*opposite_entry);
1427 if(!N::is_sign_preserving || eij<0 || eji<0)
1428 maxValue_ = max(maxValue_,
1429 eij /diagonal_ * eji/
1430 norm_(matrix_->operator[](
col.index())[
col.index()]));
1433 template<
class M,
class N>
1437 real_type eij = norm_(*
col);
1439 matrix_->operator[](
col.index()).find(row_);
1441 if ( opposite_entry == matrix_->operator[](
col.index()).end() )
1446 real_type eji = norm_(*opposite_entry);
1448 if(!N::is_sign_preserving || (eij<0 || eji<0))
1449 if(eji / norm_(matrix_->operator[](edge.target())[edge.target()]) *
1450 eij/ diagonal_ > alpha() * maxValue_) {
1451 edge.properties().setDepends();
1452 edge.properties().setInfluences();
1453 typename G::EdgeProperties& other = graph.getEdgeProperties(edge.target(), edge.source());
1454 other.setInfluences();
1459 template<
class M,
class N>
1462 return maxValue_ < beta();
1466 template<
class M,
class N>
1472 template<
class M,
class N>
1476 DUNE_UNUSED_PARAMETER(row);
1477 maxValue_ = min(- std::numeric_limits<real_type>::max(), std::numeric_limits<real_type>::min());
1479 diagonal_ = norm_(matrix_->operator[](row_)[row_]);
1482 template<
class M,
class N>
1486 maxValue_ = max(maxValue_, -norm_(*
col));
1489 template<
class M,
class N>
1493 if(-norm_(*
col) >= maxValue_ * alpha()) {
1494 edge.properties().setDepends();
1495 typedef typename G::EdgeDescriptor ED;
1496 ED e= graph.findEdge(edge.target(), edge.source());
1497 if(e!=std::numeric_limits<ED>::max())
1499 typename G::EdgeProperties& other = graph.getEdgeProperties(e);
1500 other.setInfluences();
1505 template<
class M,
class N>
1508 return maxValue_ < beta() * diagonal_;
1511 template<
class G,
class S>
1513 VertexSet& connected, std::vector<Vertex>& front)
1514 : vertices_(), id_(-1), graph_(graph), aggregates_(aggregates),
1515 connected_(connected), front_(front)
1518 template<
class G,
class S>
1526 throw "Not yet implemented";
1534 template<
class G,
class S>
1537 dvverb<<
"Connected cleared"<<std::endl;
1540 connected_.insert(vertex);
1541 dvverb <<
" Inserting "<<vertex<<
" size="<<connected_.size();
1547 template<
class G,
class S>
1550 vertices_.insert(vertex);
1551 aggregates_[vertex]=id_;
1553 front_.erase(std::lower_bound(front_.begin(), front_.end(), vertex));
1557 const iterator end = graph_.endEdges(vertex);
1558 for(iterator edge = graph_.beginEdges(vertex); edge != end; ++edge) {
1559 dvverb <<
" Inserting "<<aggregates_[edge.target()];
1560 connected_.insert(aggregates_[edge.target()]);
1561 dvverb <<
" size="<<connected_.size();
1563 !graph_.getVertexProperties(edge.target()).front())
1565 front_.push_back(edge.target());
1566 graph_.getVertexProperties(edge.target()).setFront();
1570 std::sort(front_.begin(), front_.end());
1573 template<
class G,
class S>
1577 std::size_t oldsize = vertices_.size();
1579 typedef typename std::vector<Vertex>::iterator Iterator;
1581 typedef typename VertexSet::iterator SIterator;
1583 SIterator pos=vertices_.begin();
1584 std::vector<Vertex> newFront;
1585 newFront.reserve(front_.capacity());
1587 std::set_difference(front_.begin(), front_.end(), vertices.begin(), vertices.end(),
1588 std::back_inserter(newFront));
1591 for(Iterator vertex=vertices.begin(); vertex != vertices.end(); ++vertex)
1593 pos=vertices_.insert(pos,*vertex);
1594 vertices_.insert(*vertex);
1595 graph_.getVertexProperties(*vertex).resetFront();
1596 aggregates_[*vertex]=id_;
1599 const iterator end = graph_.endEdges(*vertex);
1600 for(iterator edge = graph_.beginEdges(*vertex); edge != end; ++edge) {
1601 dvverb <<
" Inserting "<<aggregates_[edge.target()];
1602 connected_.insert(aggregates_[edge.target()]);
1604 !graph_.getVertexProperties(edge.target()).front())
1606 front_.push_back(edge.target());
1607 graph_.getVertexProperties(edge.target()).setFront();
1609 dvverb <<
" size="<<connected_.size();
1613 std::sort(front_.begin(), front_.end());
1614 assert(oldsize+vertices.size()==vertices_.size());
1616 template<
class G,
class S>
1624 template<
class G,
class S>
1625 inline typename Aggregate<G,S>::VertexSet::size_type
1628 return vertices_.size();
1631 template<
class G,
class S>
1632 inline typename Aggregate<G,S>::VertexSet::size_type
1635 return connected_.size();
1638 template<
class G,
class S>
1644 template<
class G,
class S>
1647 return vertices_.begin();
1650 template<
class G,
class S>
1653 return vertices_.end();
1671 delete[] aggregates_;
1678 allocate(noVertices);
1691 noVertices_ = noVertices;
1693 for(std::size_t i=0; i < noVertices; i++)
1694 aggregates_[i]=UNAGGREGATED;
1700 assert(aggregates_ != 0);
1701 delete[] aggregates_;
1709 return aggregates_[v];
1716 return aggregates_[v];
1720 template<
bool reset,
class G,
class F,
class VM>
1723 const G& graph, F& aggregateVisitor,
1724 VM& visitedMap)
const
1728 DummyEdgeVisitor dummy;
1729 return breadthFirstSearch<true,reset>(start, aggregate, graph, vlist, aggregateVisitor, dummy, visitedMap);
1733 template<
bool remove,
bool reset,
class G,
class L,
class F1,
class F2,
class VM>
1738 F1& aggregateVisitor,
1739 F2& nonAggregateVisitor,
1740 VM& visitedMap)
const
1742 typedef typename L::const_iterator ListIterator;
1743 int visitedSpheres = 0;
1745 visited.push_back(start);
1746 put(visitedMap, start,
true);
1748 ListIterator current = visited.begin();
1749 ListIterator end = visited.end();
1750 std::size_t i=0, size=visited.size();
1754 while(current != end) {
1756 for(; i<size; ++current, ++i) {
1757 typedef typename G::ConstEdgeIterator EdgeIterator;
1758 const EdgeIterator endEdge = graph.endEdges(*current);
1760 for(EdgeIterator edge = graph.beginEdges(*current);
1761 edge != endEdge; ++edge) {
1763 if(aggregates_[edge.target()]==aggregate) {
1764 if(!
get(visitedMap, edge.target())) {
1765 put(visitedMap, edge.target(),
true);
1766 visited.push_back(edge.target());
1767 aggregateVisitor(edge);
1770 nonAggregateVisitor(edge);
1773 end = visited.end();
1774 size = visited.size();
1780 for(current = visited.begin(); current != end; ++current)
1781 put(visitedMap, *current,
false);
1787 return visitedSpheres;
1792 : graph_(0), aggregate_(0), front_(), connected_(), size_(-1)
1801 template<
class G,
class C>
1803 const typename C::Matrix& matrix,
1804 C criterion,
bool firstlevel)
1807 typedef typename C::Matrix Matrix;
1808 typedef typename G::VertexIterator VertexIterator;
1810 criterion.init(&matrix);
1812 for(VertexIterator vertex = graph.begin(); vertex != graph.end(); ++vertex) {
1815 const Row& row = matrix[*vertex];
1820 criterion.initRow(row, *vertex);
1825 ColIterator end = row.end();
1826 typename FieldTraits<typename Matrix::field_type>::real_type absoffdiag=0.;
1830 for(ColIterator
col = row.begin();
col != end; ++
col)
1831 if(
col.index()!=*vertex) {
1832 criterion.examine(
col);
1833 absoffdiag = max(absoffdiag, Impl::asMatrix(*col).frobenius_norm());
1837 vertex.properties().setExcludedBorder();
1840 for(ColIterator
col = row.begin();
col != end; ++
col)
1841 if(
col.index()!=*vertex)
1842 criterion.examine(
col);
1848 if(criterion.isIsolated()) {
1850 vertex.properties().setIsolated();
1853 typedef typename G::EdgeIterator EdgeIterator;
1855 EdgeIterator eEnd = vertex.end();
1856 ColIterator
col = matrix[*vertex].begin();
1858 for(EdgeIterator edge = vertex.begin(); edge!= eEnd; ++edge, ++
col) {
1860 while(
col.index()!=edge.target())
1862 criterion.examine(graph, edge,
col);
1872 inline Aggregator<G>::AggregateVisitor<V>::AggregateVisitor(
const AggregatesMap<Vertex>& aggregates,
1874 : aggregates_(aggregates), aggregate_(aggregate), visitor_(&visitor)
1881 if(aggregates_[edge.target()]==aggregate_)
1882 visitor_->operator()(edge);
1887 inline void Aggregator<G>::visitAggregateNeighbours(
const Vertex& vertex,
1889 const AggregatesMap<Vertex>& aggregates,
1893 AggregateVisitor<V> v(aggregates, aggregate, visitor);
1899 inline Aggregator<G>::Counter::Counter()
1904 inline void Aggregator<G>::Counter::increment()
1910 inline void Aggregator<G>::Counter::decrement()
1915 inline int Aggregator<G>::Counter::value()
1923 if(edge.properties().isTwoWay())
1924 Counter::increment();
1929 const AggregatesMap<Vertex>& aggregates)
const
1931 TwoWayCounter counter;
1932 visitAggregateNeighbours(vertex, aggregate, aggregates, counter);
1933 return counter.value();
1938 const AggregatesMap<Vertex>& aggregates)
const
1940 OneWayCounter counter;
1941 visitAggregateNeighbours(vertex, aggregate, aggregates, counter);
1942 return counter.value();
1948 if(edge.properties().isOneWay())
1949 Counter::increment();
1953 inline Aggregator<G>::ConnectivityCounter::ConnectivityCounter(
const VertexSet& connected,
1954 const AggregatesMap<Vertex>& aggregates)
1955 : Counter(), connected_(connected), aggregates_(aggregates)
1964 Counter::increment();
1966 Counter::increment();
1967 Counter::increment();
1972 inline double Aggregator<G>::connectivity(
const Vertex& vertex,
const AggregatesMap<Vertex>& aggregates)
const
1974 ConnectivityCounter counter(connected_, aggregates);
1976 return (
double)counter.value()/noNeighbours;
1980 inline Aggregator<G>::DependencyCounter::DependencyCounter()
1987 if(edge.properties().depends())
1988 Counter::increment();
1989 if(edge.properties().influences())
1990 Counter::increment();
1994 int Aggregator<G>::unusedNeighbours(
const Vertex& vertex,
const AggregatesMap<Vertex>& aggregates)
const
2000 std::pair<int,int> Aggregator<G>::neighbours(
const Vertex& vertex,
2002 const AggregatesMap<Vertex>& aggregates)
const
2004 DependencyCounter unused, aggregated;
2005 typedef AggregateVisitor<DependencyCounter> Counter;
2006 typedef std::tuple<Counter,Counter> CounterTuple;
2009 return std::make_pair(unused.value(), aggregated.value());
2014 int Aggregator<G>::aggregateNeighbours(
const Vertex& vertex,
const AggregateDescriptor& aggregate,
const AggregatesMap<Vertex>& aggregates)
const
2016 DependencyCounter counter;
2017 visitAggregateNeighbours(vertex, aggregate, aggregates, counter);
2018 return counter.value();
2022 std::size_t Aggregator<G>::distance(
const Vertex& vertex,
const AggregatesMap<Vertex>& aggregates)
2025 typename PropertyMapTypeSelector<VertexVisitedTag,G>::Type visitedMap =
get(VertexVisitedTag(), *graph_);
2027 typename AggregatesMap<Vertex>::DummyEdgeVisitor dummy;
2028 return aggregates.template breadthFirstSearch<true,true>(vertex,
2029 aggregate_->
id(), *graph_,
2030 vlist, dummy, dummy, visitedMap);
2034 inline Aggregator<G>::FrontMarker::FrontMarker(std::vector<Vertex>& front,
MatrixGraph& graph)
2035 : front_(front), graph_(graph)
2041 Vertex target = edge.target();
2043 if(!graph_.getVertexProperties(target).front()) {
2044 front_.push_back(target);
2045 graph_.getVertexProperties(target).setFront();
2050 inline bool Aggregator<G>::admissible(
const Vertex& vertex,
const AggregateDescriptor& aggregate,
const AggregatesMap<Vertex>& aggregates)
const
2053 Dune::dvverb<<
" Admissible not yet implemented!"<<std::endl;
2060 Iterator vend = graph_->endEdges(vertex);
2061 for(Iterator edge = graph_->beginEdges(vertex); edge != vend; ++edge) {
2063 if(edge.properties().isStrong()
2064 && aggregates[edge.target()]==aggregate)
2067 Iterator edge1 = edge;
2068 for(++edge1; edge1 != vend; ++edge1) {
2070 if(edge1.properties().isStrong()
2071 && aggregates[edge.target()]==aggregate)
2076 Iterator v2end = graph_->endEdges(edge.target());
2077 for(Iterator edge2 = graph_->beginEdges(edge.target()); edge2 != v2end; ++edge2) {
2078 if(edge2.target()==edge1.target() &&
2079 edge2.properties().isStrong()) {
2095 vend = graph_->endEdges(vertex);
2096 for(Iterator edge = graph_->beginEdges(vertex); edge != vend; ++edge) {
2098 if(edge.properties().isStrong()
2099 && aggregates[edge.target()]==aggregate)
2102 Iterator v1end = graph_->endEdges(edge.target());
2104 for(Iterator edge1=graph_->beginEdges(edge.target()); edge1 != v1end; ++edge1) {
2106 if(edge1.properties().isStrong()
2107 && aggregates[edge1.target()]==aggregate)
2111 Iterator v2end = graph_->endEdges(vertex);
2112 for(Iterator edge2 = graph_->beginEdges(vertex); edge2 != v2end; ++edge2) {
2113 if(edge2.target()==edge1.target()) {
2114 if(edge2.properties().isStrong())
2131 void Aggregator<G>::unmarkFront()
2133 typedef typename std::vector<Vertex>::const_iterator Iterator;
2135 for(Iterator vertex=front_.begin(); vertex != front_.end(); ++vertex)
2136 graph_->getVertexProperties(*vertex).resetFront();
2143 Aggregator<G>::nonisoNeighbourAggregate(
const Vertex& vertex,
2144 const AggregatesMap<Vertex>& aggregates,
2145 SLList<Vertex>& neighbours)
const
2148 Iterator end=graph_->beginEdges(vertex);
2151 for(Iterator edge=graph_->beginEdges(vertex); edge!=end; ++edge)
2154 neighbours.push_back(aggregates[edge.target()]);
2159 inline typename G::VertexDescriptor Aggregator<G>::mergeNeighbour(
const Vertex& vertex,
const AggregatesMap<Vertex>& aggregates)
const
2163 Iterator end = graph_->endEdges(vertex);
2164 for(Iterator edge = graph_->beginEdges(vertex); edge != end; ++edge) {
2166 graph_->getVertexProperties(edge.target()).isolated() == graph_->getVertexProperties(edge.source()).isolated()) {
2167 if( graph_->getVertexProperties(vertex).isolated() ||
2168 ((edge.properties().depends() || edge.properties().influences())
2169 && admissible(vertex, aggregates[edge.target()], aggregates)))
2170 return edge.target();
2177 Aggregator<G>::FrontNeighbourCounter::FrontNeighbourCounter(
const MatrixGraph& graph)
2178 : Counter(), graph_(graph)
2184 if(graph_.getVertexProperties(edge.target()).front())
2185 Counter::increment();
2189 int Aggregator<G>::noFrontNeighbours(
const Vertex& vertex)
const
2191 FrontNeighbourCounter counter(*graph_);
2193 return counter.value();
2196 inline bool Aggregator<G>::connected(
const Vertex& vertex,
2198 const AggregatesMap<Vertex>& aggregates)
const
2200 typedef typename G::ConstEdgeIterator iterator;
2201 const iterator end = graph_->endEdges(vertex);
2202 for(iterator edge = graph_->beginEdges(vertex); edge != end; ++edge)
2203 if(aggregates[edge.target()]==aggregate)
2208 inline bool Aggregator<G>::connected(
const Vertex& vertex,
2209 const SLList<AggregateDescriptor>& aggregateList,
2210 const AggregatesMap<Vertex>& aggregates)
const
2212 typedef typename SLList<AggregateDescriptor>::const_iterator Iter;
2213 for(Iter i=aggregateList.begin(); i!=aggregateList.end(); ++i)
2214 if(connected(vertex, *i, aggregates))
2221 void Aggregator<G>::growIsolatedAggregate(
const Vertex& seed,
const AggregatesMap<Vertex>& aggregates,
const C& c)
2223 SLList<Vertex> connectedAggregates;
2224 nonisoNeighbourAggregate(seed, aggregates,connectedAggregates);
2226 while(aggregate_->
size()< c.minAggregateSize() && aggregate_->
connectSize() < c.maxConnectivity()) {
2228 std::size_t maxFrontNeighbours=0;
2232 typedef typename std::vector<Vertex>::const_iterator Iterator;
2234 for(Iterator vertex = front_.begin(); vertex != front_.end(); ++vertex) {
2235 if(distance(*vertex, aggregates)>c.maxDistance())
2238 if(connectedAggregates.size()>0) {
2242 if(!connected(*vertex, connectedAggregates, aggregates))
2246 double con = connectivity(*vertex, aggregates);
2249 std::size_t frontNeighbours = noFrontNeighbours(*vertex);
2251 if(frontNeighbours >= maxFrontNeighbours) {
2252 maxFrontNeighbours = frontNeighbours;
2253 candidate = *vertex;
2255 }
else if(con > maxCon) {
2257 maxFrontNeighbours = noFrontNeighbours(*vertex);
2258 candidate = *vertex;
2265 aggregate_->
add(candidate);
2271 void Aggregator<G>::growAggregate(
const Vertex& seed,
const AggregatesMap<Vertex>& aggregates,
const C& c)
2275 std::size_t distance_ =0;
2276 while(aggregate_->
size() < c.minAggregateSize()&& distance_<c.maxDistance()) {
2277 int maxTwoCons=0, maxOneCons=0, maxNeighbours=-1;
2280 std::vector<Vertex> candidates;
2281 candidates.reserve(30);
2283 typedef typename std::vector<Vertex>::const_iterator Iterator;
2285 for(Iterator vertex = front_.begin(); vertex != front_.end(); ++vertex) {
2287 if(graph_->getVertexProperties(*vertex).isolated())
2290 int twoWayCons = twoWayConnections(*vertex, aggregate_->
id(), aggregates);
2293 if( maxTwoCons == twoWayCons && twoWayCons > 0) {
2294 double con = connectivity(*vertex, aggregates);
2297 int neighbours = noFrontNeighbours(*vertex);
2299 if(neighbours > maxNeighbours) {
2300 maxNeighbours = neighbours;
2302 candidates.push_back(*vertex);
2304 candidates.push_back(*vertex);
2306 }
else if( con > maxCon) {
2308 maxNeighbours = noFrontNeighbours(*vertex);
2310 candidates.push_back(*vertex);
2312 }
else if(twoWayCons > maxTwoCons) {
2313 maxTwoCons = twoWayCons;
2314 maxCon = connectivity(*vertex, aggregates);
2315 maxNeighbours = noFrontNeighbours(*vertex);
2317 candidates.push_back(*vertex);
2320 maxOneCons = std::numeric_limits<int>::max();
2329 int oneWayCons = oneWayConnections(*vertex, aggregate_->
id(), aggregates);
2334 if(!admissible(*vertex, aggregate_->
id(), aggregates))
2337 if( maxOneCons == oneWayCons && oneWayCons > 0) {
2338 double con = connectivity(*vertex, aggregates);
2341 int neighbours = noFrontNeighbours(*vertex);
2343 if(neighbours > maxNeighbours) {
2344 maxNeighbours = neighbours;
2346 candidates.push_back(*vertex);
2348 if(neighbours==maxNeighbours)
2350 candidates.push_back(*vertex);
2353 }
else if( con > maxCon) {
2355 maxNeighbours = noFrontNeighbours(*vertex);
2357 candidates.push_back(*vertex);
2359 }
else if(oneWayCons > maxOneCons) {
2360 maxOneCons = oneWayCons;
2361 maxCon = connectivity(*vertex, aggregates);
2362 maxNeighbours = noFrontNeighbours(*vertex);
2364 candidates.push_back(*vertex);
2369 if(!candidates.size())
2371 distance_=distance(seed, aggregates);
2372 candidates.resize(min(candidates.size(), c.maxAggregateSize()-
2373 aggregate_->
size()));
2374 aggregate_->
add(candidates);
2378 template<
typename V>
2379 template<
typename M,
typename G,
typename C>
2383 Aggregator<G> aggregator;
2384 return aggregator.build(matrix, graph, *
this, criterion, finestLevel);
2388 template<
class M,
class C>
2389 std::tuple<int,int,int,int>
Aggregator<G>::build(
const M& m, G& graph, AggregatesMap<Vertex>& aggregates,
const C& c,
2395 Stack stack_(graph, *
this, aggregates);
2399 aggregate_ =
new Aggregate<G,VertexSet>(graph, aggregates, connected_, front_);
2406 dverb<<
"Build dependency took "<< watch.elapsed()<<
" seconds."<<std::endl;
2407 int noAggregates, conAggregates, isoAggregates, oneAggregates;
2408 std::size_t maxA=0, minA=1000000, avg=0;
2409 int skippedAggregates;
2410 noAggregates = conAggregates = isoAggregates = oneAggregates =
2411 skippedAggregates = 0;
2414 Vertex seed = stack_.pop();
2416 if(seed == Stack::NullEntry)
2421 if((noAggregates+1)%10000 == 0)
2425 if(graph.getVertexProperties(seed).excludedBorder()) {
2427 ++skippedAggregates;
2431 if(graph.getVertexProperties(seed).isolated()) {
2432 if(c.skipIsolated()) {
2435 ++skippedAggregates;
2439 aggregate_->
seed(seed);
2440 growIsolatedAggregate(seed, aggregates, c);
2443 aggregate_->
seed(seed);
2444 growAggregate(seed, aggregates, c);
2448 while(!(graph.getVertexProperties(seed).isolated()) && aggregate_->
size() < c.maxAggregateSize()) {
2450 std::vector<Vertex> candidates;
2451 candidates.reserve(30);
2453 typedef typename std::vector<Vertex>::const_iterator Iterator;
2455 for(Iterator vertex = front_.begin(); vertex != front_.end(); ++vertex) {
2457 if(graph.getVertexProperties(*vertex).isolated())
2460 if(twoWayConnections( *vertex, aggregate_->
id(), aggregates) == 0 &&
2461 (oneWayConnections( *vertex, aggregate_->
id(), aggregates) == 0 ||
2462 !admissible( *vertex, aggregate_->
id(), aggregates) ))
2465 std::pair<int,int> neighbourPair=neighbours(*vertex, aggregate_->
id(),
2471 if(neighbourPair.first >= neighbourPair.second)
2474 if(distance(*vertex, aggregates) > c.maxDistance())
2476 candidates.push_back(*vertex);
2480 if(!candidates.size())
break;
2482 candidates.resize(min(candidates.size(), c.maxAggregateSize()-
2483 aggregate_->
size()));
2484 aggregate_->
add(candidates);
2489 if(aggregate_->
size()==1 && c.maxAggregateSize()>1) {
2490 if(!graph.getVertexProperties(seed).isolated()) {
2491 Vertex mergedNeighbour = mergeNeighbour(seed, aggregates);
2495 aggregates[seed] = aggregates[mergedNeighbour];
2499 minA=min(minA,
static_cast<std::size_t
>(1));
2500 maxA=max(maxA,
static_cast<std::size_t
>(1));
2506 minA=min(minA,
static_cast<std::size_t
>(1));
2507 maxA=max(maxA,
static_cast<std::size_t
>(1));
2513 avg+=aggregate_->
size();
2514 minA=min(minA,aggregate_->
size());
2515 maxA=max(maxA,aggregate_->
size());
2516 if(graph.getVertexProperties(seed).isolated())
2524 Dune::dinfo<<
"connected aggregates: "<<conAggregates;
2525 Dune::dinfo<<
" isolated aggregates: "<<isoAggregates;
2526 if(conAggregates+isoAggregates>0)
2527 Dune::dinfo<<
" one node aggregates: "<<oneAggregates<<
" min size="
2528 <<minA<<
" max size="<<maxA
2529 <<
" avg="<<avg/(conAggregates+isoAggregates)<<std::endl;
2532 return std::make_tuple(conAggregates+isoAggregates,isoAggregates,
2533 oneAggregates,skippedAggregates);
2538 Aggregator<G>::Stack::Stack(
const MatrixGraph& graph,
const Aggregator<G>& aggregatesBuilder,
2539 const AggregatesMap<Vertex>& aggregates)
2540 : graph_(graph), aggregatesBuilder_(aggregatesBuilder), aggregates_(aggregates), begin_(graph.begin()), end_(graph.end())
2546 Aggregator<G>::Stack::~Stack()
2554 = std::numeric_limits<typename G::VertexDescriptor>::max();
2557 inline typename G::VertexDescriptor Aggregator<G>::Stack::pop()
2563 typename G::VertexDescriptor current=*begin_;
2577 std::ios_base::fmtflags oldOpts=os.flags();
2579 os.setf(std::ios_base::right, std::ios_base::adjustfield);
2584 for(
int i=0; i< n*m; i++)
2585 maxVal=max(maxVal, aggregates[i]);
2587 for(
int i=10; i < 1000000; i*=10)
2593 for(
int j=0, entry=0; j < m; j++) {
2594 for(
int i=0; i<n; i++, entry++) {
2596 os<<aggregates[entry]<<
" ";