4 #ifndef DUNE_AMG_AGGREGATES_HH
5 #define DUNE_AMG_AGGREGATES_HH
13 #include<dune/common/timer.hh>
14 #include<dune/common/tuples.hh>
15 #include<dune/common/stdstreams.hh>
16 #include<dune/common/poolallocator.hh>
17 #include<dune/common/sllist.hh>
80 this->setMaxDistance(diameter-1);
85 this->setMaxDistance(this->maxDistance()+diameter-1);
87 this->setMinAggregateSize(csize);
88 this->setMaxAggregateSize(static_cast<std::size_t>(csize*1.5));
104 this->setMaxDistance(this->maxDistance()+dim-1);
109 std::ostream& operator<<(std::ostream& os, const AggregationCriterion<T>& criterion)
111 os<<
"{ maxdistance="<<criterion.maxDistance()<<
" minAggregateSize="
112 <<criterion.minAggregateSize()<<
" maxAggregateSize="<<criterion.maxAggregateSize()
113 <<
" connectivity="<<criterion.maxConnectivity()<<
" debugLevel="<<criterion.debugLevel()<<
"}";
120 template<
class M,
class N>
151 void examine(G& graph,
const typename G::EdgeIterator& edge,
const ColIter&
col);
179 template<
class M,
class N>
210 void examine(G& graph,
const typename G::EdgeIterator& edge,
const ColIter&
col);
283 return m.infinity_norm();
300 return m.frobenius_norm();
328 template<
class M,
class Norm>
351 template<
class M,
class Norm>
362 template<
class G>
class Aggregator;
414 template<
class EdgeIterator>
448 template<
class M,
class G,
class C>
449 tuple<int,int,int,int>
buildAggregates(
const M& matrix, G& graph,
const C& criterion,
471 template<
bool reset,
class G,
class F,
class VM>
476 VM& visitedMap)
const;
501 template<
bool remove,
bool reset,
class G,
class L,
class F1,
class F2,
class VM>
504 const G& graph, L& visited, F1& aggregateVisitor,
505 F2& nonAggregateVisitor,
506 VM& visitedMap)
const;
570 AggregatesMap<V>& operator=(
const AggregatesMap<V>& map)
584 std::size_t noVertices_;
590 template<
class G,
class C>
592 const typename C::Matrix& matrix,
600 template<
class G,
class S>
678 typename VertexList::size_type
size();
770 template<
class M,
class C>
771 tuple<int,int,int,int>
build(
const M& m, G& graph,
779 typedef PoolAllocator<Vertex,100> Allocator;
784 typedef SLList<Vertex,Allocator> VertexList;
789 typedef std::set<Vertex,std::less<Vertex>,Allocator> VertexSet;
794 typedef std::size_t* SphereMap;
814 VertexSet connected_;
833 bool push(
const Vertex& v);
855 void localPush(
const Vertex& v);
880 class AggregateVisitor
940 class FrontNeighbourCounter :
public Counter
959 int noFrontNeighbours(
const Vertex& vertex)
const;
964 class TwoWayCounter :
public Counter
982 const AggregatesMap<Vertex>& aggregates)
const;
987 class OneWayCounter :
public Counter
1005 const AggregatesMap<Vertex>& aggregates)
const;
1013 class ConnectivityCounter :
public Counter
1022 ConnectivityCounter(
const VertexSet& connected,
const AggregatesMap<Vertex>& aggregates);
1028 const VertexSet& connected_;
1030 const AggregatesMap<Vertex>& aggregates_;
1045 double connectivity(
const Vertex& vertex,
const AggregatesMap<Vertex>& aggregates)
const;
1054 const AggregatesMap<Vertex>& aggregates)
const;
1063 bool connected(
const Vertex& vertex,
const SLList<AggregateDescriptor>& aggregateList,
1064 const AggregatesMap<Vertex>& aggregates)
const;
1073 class DependencyCounter:
public Counter
1079 DependencyCounter();
1099 FrontMarker(VertexList& front,
MatrixGraph& graph);
1116 void markFront(
const AggregatesMap<Vertex>& aggregates);
1137 int unusedNeighbours(
const Vertex& vertex,
const AggregatesMap<Vertex>& aggregates)
const;
1152 std::pair<int,int> neighbours(
const Vertex& vertex,
1154 const AggregatesMap<Vertex>& aggregates)
const;
1171 int aggregateNeighbours(
const Vertex& vertex,
const AggregateDescriptor& aggregate,
const AggregatesMap<Vertex>& aggregates)
const;
1180 bool admissible(
const Vertex& vertex,
const AggregateDescriptor& aggregate,
const AggregatesMap<Vertex>& aggregates)
const;
1188 void seedFromFront(
Stack& stack,
bool isolated);
1197 std::size_t distance(
const Vertex& vertex,
const AggregatesMap<Vertex>& aggregates);
1207 Vertex mergeNeighbour(
const Vertex& vertex,
const AggregatesMap<Vertex>& aggregates)
const;
1217 void nonisoNeighbourAggregate(
const Vertex& vertex,
1218 const AggregatesMap<Vertex>& aggregates,
1219 SLList<Vertex>& neighbours)
const;
1229 void growAggregate(
const Vertex& vertex,
const AggregatesMap<Vertex>& aggregates,
const C& c);
1231 void growIsolatedAggregate(
const Vertex& vertex,
const AggregatesMap<Vertex>& aggregates,
const C& c);
1236 template<
class M,
class N>
1242 template<
class M,
class N>
1245 maxValue_ = std::min(- std::numeric_limits<typename Matrix::field_type>::max(), std::numeric_limits<typename Matrix::field_type>::min());
1247 diagonal_ = norm_(matrix_->operator[](row_)[row_]);
1250 template<
class M,
class N>
1257 if(!N::is_sign_preserving || eij<0 || eji<0)
1258 maxValue_ = std::max(maxValue_,
1259 eij /diagonal_ * eji/
1260 norm_(matrix_->operator[](col.index())[col.index()]));
1263 template<
class M,
class N>
1271 if(!N::is_sign_preserving || (eij<0 || eji<0))
1272 if(eji / norm_(matrix_->operator[](edge.target())[edge.target()]) *
1273 eij/ diagonal_ > alpha() * maxValue_){
1274 edge.properties().setDepends();
1275 edge.properties().setInfluences();
1277 typename G::EdgeProperties& other = graph.getEdgeProperties(edge.target(), edge.source());
1278 other.setInfluences();
1283 template<
class M,
class N>
1286 return maxValue_ < beta();
1290 template<
class M,
class N>
1296 template<
class M,
class N>
1299 maxValue_ = std::min(- std::numeric_limits<typename Matrix::field_type>::max(), std::numeric_limits<typename Matrix::field_type>::min());
1301 diagonal_ = norm_(matrix_->operator[](row_)[row_]);
1304 template<
class M,
class N>
1307 maxValue_ = std::max(maxValue_,
1311 template<
class M,
class N>
1315 if(-norm_(*col) >= maxValue_ * alpha()){
1316 edge.properties().setDepends();
1317 typedef typename G::EdgeDescriptor ED;
1318 ED e= graph.findEdge(edge.target(), edge.source());
1319 if(e!=std::numeric_limits<ED>::max())
1321 typename G::EdgeProperties& other = graph.getEdgeProperties(e);
1322 other.setInfluences();
1327 template<
class M,
class N>
1330 return maxValue_ < beta() * diagonal_;
1333 template<
class G,
class S>
1335 VertexSet& connected)
1336 : vertices_(), id_(-1), graph_(graph), aggregates_(aggregates),
1337 connected_(connected)
1340 template<
class G,
class S>
1341 void Aggregate<G,S>::reconstruct(
const Vertex& vertex)
1343 vertices_.push_back(vertex);
1344 typedef typename VertexList::const_iterator iterator;
1345 iterator begin = vertices_.begin();
1346 iterator end = vertices_.end();
1347 throw "Not yet implemented";
1355 template<
class G,
class S>
1356 inline void Aggregate<G,S>::seed(
const Vertex& vertex)
1358 dvverb<<
"Connected cleared"<<std::endl;
1361 connected_.insert(vertex);
1362 dvverb <<
" Inserting "<<vertex<<
" size="<<connected_.size();
1368 template<
class G,
class S>
1369 inline void Aggregate<G,S>::add(
const Vertex& vertex)
1371 vertices_.push_back(vertex);
1372 aggregates_[vertex]=id_;
1375 const iterator end = graph_.endEdges(vertex);
1376 for(iterator edge = graph_.beginEdges(vertex); edge != end; ++edge){
1377 dvverb <<
" Inserting "<<aggregates_[edge.target()];
1378 connected_.insert(aggregates_[edge.target()]);
1379 dvverb <<
" size="<<connected_.size();
1383 template<
class G,
class S>
1384 inline void Aggregate<G,S>::clear()
1391 template<
class G,
class S>
1392 inline typename Aggregate<G,S>::VertexList::size_type
1393 Aggregate<G,S>::size()
1395 return vertices_.size();
1398 template<
class G,
class S>
1399 inline typename Aggregate<G,S>::VertexList::size_type
1400 Aggregate<G,S>::connectSize()
1402 return connected_.size();
1405 template<
class G,
class S>
1406 inline int Aggregate<G,S>::id()
1411 template<
class G,
class S>
1412 inline typename Aggregate<G,S>::const_iterator Aggregate<G,S>::begin()
const
1414 return vertices_.begin();
1417 template<
class G,
class S>
1418 inline typename Aggregate<G,S>::const_iterator Aggregate<G,S>::end()
const
1420 return vertices_.end();
1424 const V AggregatesMap<V>::UNAGGREGATED = std::numeric_limits<V>::max();
1427 const V AggregatesMap<V>::ISOLATED = std::numeric_limits<V>::max()-1;
1430 AggregatesMap<V>::AggregatesMap()
1435 AggregatesMap<V>::~AggregatesMap()
1438 delete[] aggregates_;
1443 inline AggregatesMap<V>::AggregatesMap(std::size_t noVertices)
1445 allocate(noVertices);
1449 inline std::size_t AggregatesMap<V>::noVertices()
const
1455 inline void AggregatesMap<V>::allocate(std::size_t noVertices)
1458 noVertices_ = noVertices;
1460 for(std::size_t i=0; i < noVertices; i++)
1461 aggregates_[i]=UNAGGREGATED;
1465 inline void AggregatesMap<V>::free()
1467 assert(aggregates_ != 0);
1468 delete[] aggregates_;
1473 inline typename AggregatesMap<V>::AggregateDescriptor&
1474 AggregatesMap<V>::operator[](
const VertexDescriptor& v)
1476 return aggregates_[v];
1480 inline const typename AggregatesMap<V>::AggregateDescriptor&
1481 AggregatesMap<V>::operator[](
const VertexDescriptor& v)
const
1483 return aggregates_[v];
1487 template<
bool reset,
class G,
class F,
class VM>
1488 inline std::size_t AggregatesMap<V>::breadthFirstSearch(
const V& start,
1490 const G& graph, F& aggregateVisitor,
1491 VM& visitedMap)
const
1495 DummyEdgeVisitor dummy;
1496 return breadthFirstSearch<true,reset>(start, aggregate, graph, vlist, aggregateVisitor, dummy, visitedMap);
1500 template<
bool remove,
bool reset,
class G,
class L,
class F1,
class F2,
class VM>
1501 std::size_t AggregatesMap<V>::breadthFirstSearch(
const V& start,
1505 F1& aggregateVisitor,
1506 F2& nonAggregateVisitor,
1507 VM& visitedMap)
const
1509 typedef typename L::const_iterator ListIterator;
1510 int visitedSpheres = 0;
1512 visited.push_back(start);
1513 put(visitedMap, start,
true);
1515 ListIterator current = visited.begin();
1516 ListIterator end = visited.end();
1517 std::size_t i=0, size=visited.size();
1521 while(current != end){
1523 for(;i<size; ++current, ++i){
1524 typedef typename G::ConstEdgeIterator EdgeIterator;
1525 const EdgeIterator endEdge = graph.endEdges(*current);
1527 for(EdgeIterator edge = graph.beginEdges(*current);
1528 edge != endEdge; ++edge){
1530 if(aggregates_[edge.target()]==aggregate){
1531 if(!
get(visitedMap, edge.target())){
1532 put(visitedMap, edge.target(),
true);
1533 visited.push_back(edge.target());
1534 aggregateVisitor(edge);
1537 nonAggregateVisitor(edge);
1540 end = visited.end();
1541 size = visited.size();
1547 for(current = visited.begin(); current != end; ++current)
1548 put(visitedMap, *current,
false);
1554 return visitedSpheres;
1559 : graph_(0), aggregate_(0), front_(), connected_(), size_(-1)
1568 template<
class G,
class C>
1570 const typename C::Matrix& matrix,
1571 C criterion,
bool firstlevel)
1574 typedef C Criterion;
1577 typedef typename C::Matrix Matrix;
1578 typedef typename G::VertexIterator VertexIterator;
1580 criterion.init(&matrix);
1582 for(VertexIterator vertex = graph.begin(); vertex != graph.end(); ++vertex){
1585 const Row& row = matrix[*vertex];
1590 criterion.initRow(row, *vertex);
1595 ColIterator end = row.end();
1599 for(ColIterator col = row.begin(); col != end; ++
col)
1600 if(col.index()!=*vertex){
1601 criterion.examine(col);
1602 absoffdiag=std::max(absoffdiag, col->frobenius_norm());
1606 vertex.properties().setExcludedBorder();
1609 for(ColIterator col = row.begin(); col != end; ++
col)
1610 if(col.index()!=*vertex)
1611 criterion.examine(col);
1617 if(criterion.isIsolated()){
1619 vertex.properties().setIsolated();
1622 typedef typename G::EdgeIterator EdgeIterator;
1624 EdgeIterator end = vertex.end();
1625 ColIterator col = matrix[*vertex].begin();
1627 for(EdgeIterator edge = vertex.begin(); edge!= end; ++edge, ++
col){
1629 while(col.index()!=edge.target())
1631 criterion.examine(graph, edge, col);
1641 inline Aggregator<G>::AggregateVisitor<V>::AggregateVisitor(
const AggregatesMap<Vertex>& aggregates,
1643 : aggregates_(aggregates), aggregate_(aggregate), visitor_(&visitor)
1650 if(aggregates_[edge.target()]==aggregate_)
1651 visitor_->operator()(edge);
1656 inline void Aggregator<G>::visitAggregateNeighbours(
const Vertex& vertex,
1658 const AggregatesMap<Vertex>& aggregates,
1662 AggregateVisitor<V> v(aggregates, aggregate, visitor);
1668 inline Aggregator<G>::Counter::Counter()
1673 inline void Aggregator<G>::Counter::increment()
1679 inline void Aggregator<G>::Counter::decrement()
1684 inline int Aggregator<G>::Counter::value()
1692 if(edge.properties().isTwoWay())
1693 Counter::increment();
1698 const AggregatesMap<Vertex>& aggregates)
const
1700 TwoWayCounter counter;
1701 visitAggregateNeighbours(vertex, aggregate, aggregates, counter);
1702 return counter.value();
1707 const AggregatesMap<Vertex>& aggregates)
const
1709 OneWayCounter counter;
1710 visitAggregateNeighbours(vertex, aggregate, aggregates, counter);
1711 return counter.value();
1717 if(edge.properties().isOneWay())
1718 Counter::increment();
1722 inline Aggregator<G>::ConnectivityCounter::ConnectivityCounter(
const VertexSet& connected,
1723 const AggregatesMap<Vertex>& aggregates)
1724 : Counter(), connected_(connected), aggregates_(aggregates)
1733 Counter::increment();
1735 Counter::increment();
1736 Counter::increment();
1741 inline double Aggregator<G>::connectivity(
const Vertex& vertex,
const AggregatesMap<Vertex>& aggregates)
const
1743 ConnectivityCounter counter(connected_, aggregates);
1745 return (
double)counter.value()/noNeighbours;
1749 inline Aggregator<G>::DependencyCounter::DependencyCounter()
1756 if(edge.properties().depends())
1757 Counter::increment();
1758 if(edge.properties().influences())
1759 Counter::increment();
1763 int Aggregator<G>::unusedNeighbours(
const Vertex& vertex,
const AggregatesMap<Vertex>& aggregates)
const
1769 std::pair<int,int> Aggregator<G>::neighbours(
const Vertex& vertex,
1771 const AggregatesMap<Vertex>& aggregates)
const
1773 DependencyCounter unused, aggregated;
1774 typedef AggregateVisitor<DependencyCounter> Counter;
1775 typedef tuple<Counter,Counter> CounterTuple;
1778 return std::make_pair(unused.value(), aggregated.value());
1783 int Aggregator<G>::aggregateNeighbours(
const Vertex& vertex,
const AggregateDescriptor& aggregate,
const AggregatesMap<Vertex>& aggregates)
const
1785 DependencyCounter counter;
1786 visitAggregateNeighbours(vertex, aggregate, aggregates, counter);
1787 return counter.value();
1791 std::size_t Aggregator<G>::distance(
const Vertex& vertex,
const AggregatesMap<Vertex>& aggregates)
1793 typename PropertyMapTypeSelector<VertexVisitedTag,G>::Type visitedMap =
get(VertexVisitedTag(), *graph_);
1795 typename AggregatesMap<Vertex>::DummyEdgeVisitor dummy;
1796 return aggregates.template breadthFirstSearch<true,true>(vertex,
1797 aggregate_->
id(), *graph_,
1798 vlist, dummy, dummy, visitedMap);
1802 inline Aggregator<G>::FrontMarker::FrontMarker(VertexList& front,
MatrixGraph& graph)
1803 : front_(front), graph_(graph)
1809 Vertex target = edge.target();
1811 if(!graph_.getVertexProperties(target).front()){
1812 front_.push_back(target);
1813 graph_.getVertexProperties(target).setFront();
1819 void Aggregator<G>::markFront(
const AggregatesMap<Vertex>& aggregates)
1821 assert(front_.size()==0);
1822 FrontMarker frontBuilder(front_, *graph_);
1823 typedef typename Aggregate<G,VertexSet>::const_iterator Iterator;
1825 for(Iterator vertex=aggregate_->
begin(); vertex != aggregate_->
end(); ++vertex)
1831 inline bool Aggregator<G>::admissible(
const Vertex& vertex,
const AggregateDescriptor& aggregate,
const AggregatesMap<Vertex>& aggregates)
const
1834 Dune::dvverb<<
" Admissible not yet implemented!"<<std::endl;
1841 Iterator vend = graph_->endEdges(vertex);
1842 for(Iterator edge = graph_->beginEdges(vertex); edge != vend; ++edge){
1844 if(edge.properties().isStrong()
1845 && aggregates[edge.target()]==aggregate)
1848 Iterator edge1 = edge;
1849 for(++edge1; edge1 != vend; ++edge1){
1851 if(edge1.properties().isStrong()
1852 && aggregates[edge.target()]==aggregate)
1857 Iterator v2end = graph_->endEdges(edge.target());
1858 for(Iterator edge2 = graph_->beginEdges(edge.target()); edge2 != v2end; ++edge2){
1859 if(edge2.target()==edge1.target() &&
1860 edge2.properties().isStrong()){
1874 vend = graph_->endEdges(vertex);
1875 for(Iterator edge = graph_->beginEdges(vertex); edge != vend; ++edge){
1877 if(edge.properties().isStrong()
1878 && aggregates[edge.target()]==aggregate)
1881 Iterator v1end = graph_->endEdges(edge.target());
1883 for(Iterator edge1=graph_->beginEdges(edge.target()); edge1 != v1end; ++edge1){
1885 if(edge1.properties().isStrong()
1886 && aggregates[edge1.target()]==aggregate)
1890 Iterator v2end = graph_->endEdges(vertex);
1891 for(Iterator edge2 = graph_->beginEdges(vertex); edge2 != v2end; ++edge2){
1892 if(edge2.target()==edge1.target()){
1893 if(edge2.properties().isStrong())
1908 void Aggregator<G>::unmarkFront()
1910 typedef typename VertexList::const_iterator Iterator;
1912 for(Iterator vertex=front_.begin(); vertex != front_.end(); ++vertex)
1913 graph_->getVertexProperties(*vertex).resetFront();
1920 Aggregator<G>::nonisoNeighbourAggregate(
const Vertex& vertex,
1921 const AggregatesMap<Vertex>& aggregates,
1922 SLList<Vertex>& neighbours)
const
1925 Iterator end=graph_->beginEdges(vertex);
1928 for(Iterator edge=graph_->beginEdges(vertex); edge!=end; ++edge)
1931 neighbours.push_back(aggregates[edge.target()]);
1936 inline typename G::VertexDescriptor Aggregator<G>::mergeNeighbour(
const Vertex& vertex,
const AggregatesMap<Vertex>& aggregates)
const
1940 Iterator end = graph_->endEdges(vertex);
1941 for(Iterator edge = graph_->beginEdges(vertex); edge != end; ++edge){
1943 graph_->getVertexProperties(edge.target()).isolated() == graph_->getVertexProperties(edge.source()).isolated()){
1944 if( graph_->getVertexProperties(vertex).isolated() ||
1945 ((edge.properties().depends() || edge.properties().influences())
1946 && admissible(vertex, aggregates[edge.target()], aggregates)))
1947 return edge.target();
1954 Aggregator<G>::FrontNeighbourCounter::FrontNeighbourCounter(
const MatrixGraph& graph)
1955 : Counter(), graph_(graph)
1961 if(graph_.getVertexProperties(edge.target()).front())
1962 Counter::increment();
1966 int Aggregator<G>::noFrontNeighbours(
const Vertex& vertex)
const
1968 FrontNeighbourCounter counter(*graph_);
1970 return counter.value();
1973 inline bool Aggregator<G>::connected(
const Vertex& vertex,
1975 const AggregatesMap<Vertex>& aggregates)
const
1977 typedef typename G::ConstEdgeIterator iterator;
1978 const iterator end = graph_->endEdges(vertex);
1979 for(iterator edge = graph_->beginEdges(vertex); edge != end; ++edge)
1980 if(aggregates[edge.target()]==aggregate)
1985 inline bool Aggregator<G>::connected(
const Vertex& vertex,
1986 const SLList<AggregateDescriptor>& aggregateList,
1987 const AggregatesMap<Vertex>& aggregates)
const
1989 typedef typename SLList<AggregateDescriptor>::const_iterator Iter;
1990 for(Iter i=aggregateList.begin(); i!=aggregateList.end(); ++i)
1991 if(connected(vertex, *i, aggregates))
1998 void Aggregator<G>::growIsolatedAggregate(
const Vertex& seed,
const AggregatesMap<Vertex>& aggregates,
const C& c)
2000 SLList<Vertex> connectedAggregates;
2001 nonisoNeighbourAggregate(seed, aggregates,connectedAggregates);
2003 while(aggregate_->
size()< c.minAggregateSize() && aggregate_->
connectSize() < c.maxConnectivity()){
2005 std::size_t maxFrontNeighbours=0;
2009 markFront(aggregates);
2010 typedef typename VertexList::const_iterator Iterator;
2012 for(Iterator vertex = front_.begin(); vertex != front_.end(); ++vertex){
2013 if(distance(*vertex, aggregates)>c.maxDistance())
2016 if(connectedAggregates.size()>0){
2020 if(!connected(*vertex, connectedAggregates, aggregates))
2024 double con = connectivity(*vertex, aggregates);
2027 std::size_t frontNeighbours = noFrontNeighbours(*vertex);
2029 if(frontNeighbours >= maxFrontNeighbours){
2030 maxFrontNeighbours = frontNeighbours;
2031 candidate = *vertex;
2033 }
else if(con > maxCon){
2035 maxFrontNeighbours = noFrontNeighbours(*vertex);
2036 candidate = *vertex;
2043 aggregate_->
add(candidate);
2049 void Aggregator<G>::growAggregate(
const Vertex& seed,
const AggregatesMap<Vertex>& aggregates,
const C& c)
2051 while(aggregate_->
size() < c.minAggregateSize()){
2052 int maxTwoCons=0, maxOneCons=0, maxNeighbours=-1;
2058 markFront(aggregates);
2060 typedef typename VertexList::const_iterator Iterator;
2062 for(Iterator vertex = front_.begin(); vertex != front_.end(); ++vertex){
2064 if(graph_->getVertexProperties(*vertex).isolated())
2067 int twoWayCons = twoWayConnections(*vertex, aggregate_->
id(), aggregates);
2070 if( maxTwoCons == twoWayCons && twoWayCons > 0){
2071 double con = connectivity(*vertex, aggregates);
2074 int neighbours = noFrontNeighbours(*vertex);
2076 if(neighbours > maxNeighbours){
2077 maxNeighbours = neighbours;
2079 std::size_t distance_ = distance(*vertex, aggregates);
2081 if(c.maxDistance() >= distance_){
2082 candidate = *vertex;
2085 }
else if( con > maxCon){
2087 maxNeighbours = noFrontNeighbours(*vertex);
2088 std::size_t distance_ = distance(*vertex, aggregates);
2090 if(c.maxDistance() >= distance_){
2091 candidate = *vertex;
2094 }
else if(twoWayCons > maxTwoCons){
2095 maxTwoCons = twoWayCons;
2096 maxCon = connectivity(*vertex, aggregates);
2097 maxNeighbours = noFrontNeighbours(*vertex);
2098 std::size_t distance_ = distance(*vertex, aggregates);
2100 if(c.maxDistance() >= distance_){
2101 candidate = *vertex;
2105 maxOneCons = std::numeric_limits<int>::max();
2112 int oneWayCons = oneWayConnections(*vertex, aggregate_->
id(), aggregates);
2117 if(!admissible(*vertex, aggregate_->
id(), aggregates))
2120 if( maxOneCons == oneWayCons && oneWayCons > 0){
2121 double con = connectivity(*vertex, aggregates);
2124 int neighbours = noFrontNeighbours(*vertex);
2126 if(neighbours > maxNeighbours){
2127 maxNeighbours = neighbours;
2128 std::size_t distance_ = distance(*vertex, aggregates);
2130 if(c.maxDistance() >= distance_){
2131 candidate = *vertex;
2134 }
else if( con > maxCon){
2136 maxNeighbours = noFrontNeighbours(*vertex);
2137 std::size_t distance_ = distance(*vertex, aggregates);
2138 if(c.maxDistance() >= distance_){
2139 candidate = *vertex;
2142 }
else if(oneWayCons > maxOneCons){
2143 maxOneCons = oneWayCons;
2144 maxCon = connectivity(*vertex, aggregates);
2145 maxNeighbours = noFrontNeighbours(*vertex);
2146 std::size_t distance_ = distance(*vertex, aggregates);
2148 if(c.maxDistance() >= distance_){
2149 candidate = *vertex;
2158 aggregate_->
add(candidate);
2162 template<
typename V>
2163 template<
typename M,
typename G,
typename C>
2164 tuple<int,int,int,int> AggregatesMap<V>::buildAggregates(
const M& matrix, G& graph,
const C& criterion,
2167 Aggregator<G> aggregator;
2168 return aggregator.build(matrix, graph, *
this, criterion, finestLevel);
2172 template<
class M,
class C>
2173 tuple<int,int,int,int>
Aggregator<G>::build(
const M& m, G& graph, AggregatesMap<Vertex>& aggregates,
const C& c,
2177 Stack stack_(graph, *
this, aggregates);
2181 aggregate_ =
new Aggregate<G,VertexSet>(graph, aggregates, connected_);
2188 dverb<<
"Build dependency took "<< watch.elapsed()<<
" seconds."<<std::endl;
2189 int noAggregates, conAggregates, isoAggregates, oneAggregates;
2190 std::size_t maxA=0, minA=1000000, avg=0;
2191 int skippedAggregates;
2192 noAggregates = conAggregates = isoAggregates = oneAggregates =
2193 skippedAggregates = 0;
2196 Vertex seed = stack_.pop();
2198 if(seed == Stack::NullEntry)
2203 if((noAggregates+1)%10000 == 0)
2206 aggregate_->
seed(seed);
2208 if(graph.getVertexProperties(seed).excludedBorder()){
2210 ++skippedAggregates;
2214 if(graph.getVertexProperties(seed).isolated()){
2215 if(c.skipIsolated()){
2218 ++skippedAggregates;
2222 growIsolatedAggregate(seed, aggregates, c);
2224 growAggregate(seed, aggregates, c);
2228 while(!(graph.getVertexProperties(seed).isolated()) && aggregate_->
size() < c.maxAggregateSize()){
2231 markFront(aggregates);
2235 typedef typename VertexList::const_iterator Iterator;
2237 for(Iterator vertex = front_.begin(); vertex != front_.end(); ++vertex){
2239 if(graph.getVertexProperties(*vertex).isolated())
2242 if(twoWayConnections( *vertex, aggregate_->
id(), aggregates) == 0 &&
2243 (oneWayConnections( *vertex, aggregate_->
id(), aggregates) == 0 ||
2244 !admissible( *vertex, aggregate_->
id(), aggregates) ))
2247 std::pair<int,int> neighbourPair=neighbours(*vertex, aggregate_->
id(),
2253 if(neighbourPair.first >= neighbourPair.second)
2256 if(distance(*vertex, aggregates) > c.maxDistance())
2258 candidate = *vertex;
2264 aggregate_->
add(candidate);
2269 if(aggregate_->
size()==1 && c.maxAggregateSize()>1){
2270 if(!graph.getVertexProperties(seed).isolated()){
2271 Vertex mergedNeighbour = mergeNeighbour(seed, aggregates);
2275 aggregates[seed] = aggregates[mergedNeighbour];
2277 minA=std::min(minA,static_cast<std::size_t>(1));
2278 maxA=std::max(maxA,static_cast<std::size_t>(1));
2284 minA=std::min(minA,static_cast<std::size_t>(1));
2285 maxA=std::max(maxA,static_cast<std::size_t>(1));
2291 avg+=aggregate_->
size();
2292 minA=std::min(minA,aggregate_->
size());
2293 maxA=std::max(maxA,aggregate_->
size());
2294 if(graph.getVertexProperties(seed).isolated())
2300 markFront(aggregates);
2301 seedFromFront(stack_, graph.getVertexProperties(seed).isolated());
2305 Dune::dinfo<<
"connected aggregates: "<<conAggregates;
2306 Dune::dinfo<<
" isolated aggregates: "<<isoAggregates;
2307 if(conAggregates+isoAggregates>0)
2308 Dune::dinfo<<
" one node aggregates: "<<oneAggregates<<
" min size="
2309 <<minA<<
" max size="<<maxA
2310 <<
" avg="<<avg/(conAggregates+isoAggregates)<<std::endl;
2313 return make_tuple(conAggregates+isoAggregates,isoAggregates,
2314 oneAggregates,skippedAggregates);
2318 inline void Aggregator<G>::seedFromFront(
Stack& stack_,
bool isolated)
2320 typedef typename VertexList::const_iterator Iterator;
2322 Iterator end= front_.end();
2324 for(Iterator vertex=front_.begin(); vertex != end; ++vertex,++
count)
2325 if(graph_->getVertexProperties(*vertex).isolated()==isolated)
2326 stack_.push(*vertex);
2335 const AggregatesMap<Vertex>& aggregates)
2336 : graph_(graph), aggregatesBuilder_(aggregatesBuilder), aggregates_(aggregates), size_(0), maxSize_(0), head_(0), filled_(0)
2342 Aggregator<G>::Stack::~Stack()
2344 Dune::dverb <<
"Max stack size was "<<maxSize_<<
" filled="<<filled_<<std::endl;
2349 const typename Aggregator<G>::Vertex Aggregator<G>::Stack::NullEntry
2350 = std::numeric_limits<typename G::VertexDescriptor>::max();
2353 inline bool Aggregator<G>::Stack::push(
const Vertex& v)
2363 inline void Aggregator<G>::Stack::localPush(
const Vertex& v)
2366 size_ = std::min<int>(size_+1, N);
2367 head_ = (head_+N+1)%N;
2371 void Aggregator<G>::Stack::fill()
2373 int isolated = 0, connected=0;
2378 isoumin = umin = std::numeric_limits<int>::max();
2382 const Iterator end = graph_.end();
2384 for(Iterator vertex = graph_.begin(); vertex != end; ++vertex){
2389 if(vertex.properties().isolated()){
2390 isoumin = std::min(isoumin, aggregatesBuilder_.unusedNeighbours(*vertex, aggregates_));
2393 umin = std::min(umin, aggregatesBuilder_.unusedNeighbours(*vertex, aggregates_));
2398 if(connected + isolated == 0)
2404 for(Iterator vertex = graph_.begin(); vertex != end; ++vertex)
2406 && aggregatesBuilder_.unusedNeighbours(*vertex, aggregates_) == umin)
2409 for(Iterator vertex = graph_.begin(); vertex != end; ++vertex)
2411 && aggregatesBuilder_.unusedNeighbours(*vertex, aggregates_) == isoumin)
2414 maxSize_ = std::max(size_, maxSize_);
2418 inline typename G::VertexDescriptor Aggregator<G>::Stack::pop()
2421 head_ = (head_ + N -1) % N;
2432 head_ = (head_ + N -1) % N;
2446 std::ios_base::fmtflags oldOpts=os.flags();
2448 os.setf(std::ios_base::right, std::ios_base::adjustfield);
2453 for(
int i=0; i< n*m; i++)
2454 max=std::max(max, aggregates[i]);
2456 for(
int i=10; i < 1000000; i*=10)
2462 for(
int j=0, entry=0; j < m; j++){
2463 for(
int i=0; i<n; i++, entry++){
2465 os<<aggregates[entry]<<
" ";