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>
18 #include <dune/common/unused.hh>
81 this->setMaxDistance(diameter-1);
86 this->setMaxDistance(this->maxDistance()+diameter-1);
88 this->setMinAggregateSize(csize);
89 this->setMaxAggregateSize(static_cast<std::size_t>(csize*1.5));
105 this->setMaxDistance(this->maxDistance()+dim-1);
110 std::ostream& operator<<(std::ostream& os, const AggregationCriterion<T>& criterion)
112 os<<
"{ maxdistance="<<criterion.maxDistance()<<
" minAggregateSize="
113 <<criterion.minAggregateSize()<<
" maxAggregateSize="<<criterion.maxAggregateSize()
114 <<
" connectivity="<<criterion.maxConnectivity()<<
" debugLevel="<<criterion.debugLevel()<<
"}";
129 template<
class M,
class N>
160 void examine(G& graph,
const typename G::EdgeIterator& edge,
const ColIter&
col);
183 std::vector<typename Matrix::field_type>
vals_;
184 typename std::vector<typename Matrix::field_type>::iterator
valIter_;
189 template<
class M,
class N>
195 template<
class M,
class N>
198 vals_.assign(row.size(), 0.0);
199 assert(vals_.size()==row.size());
200 valIter_=vals_.begin();
202 maxValue_ = std::min(- std::numeric_limits<typename Matrix::field_type>::max(), std::numeric_limits<typename Matrix::field_type>::min());
203 diagonal_=norm_(row[index]);
207 template<
class M,
class N>
212 if(!N::is_sign_preserving || eij<0)
214 *valIter_ = eij/diagonal_*eij/norm_(matrix_->operator[](col.index())[col.index()]);
215 maxValue_ = std::max(maxValue_, *valIter_);
221 template<
class M,
class N>
225 if(*valIter_ > alpha() * maxValue_) {
226 edge.properties().setDepends();
227 edge.properties().setInfluences();
232 template<
class M,
class N>
237 valIter_=vals_.begin();
238 return maxValue_ < beta();
244 template<
class M,
class N>
275 void examine(G& graph,
const typename G::EdgeIterator& edge,
const ColIter&
col);
303 template<
class M,
class N>
334 void examine(G& graph,
const typename G::EdgeIterator& edge,
const ColIter&
col);
407 return m.infinity_norm();
424 return m.frobenius_norm();
452 template<
class M,
class Norm>
475 template<
class M,
class Norm>
538 template<
class EdgeIterator>
541 DUNE_UNUSED_PARAMETER(edge);
574 template<
class M,
class G,
class C>
575 tuple<int,int,int,int>
buildAggregates(
const M& matrix, G& graph,
const C& criterion,
597 template<
bool reset,
class G,
class F,
class VM>
602 VM& visitedMap)
const;
627 template<
bool remove,
bool reset,
class G,
class L,
class F1,
class F2,
class VM>
630 const G& graph, L& visited, F1& aggregateVisitor,
631 F2& nonAggregateVisitor,
632 VM& visitedMap)
const;
696 AggregatesMap<V>& operator=(
const AggregatesMap<V>& map)
710 std::size_t noVertices_;
716 template<
class G,
class C>
718 const typename C::Matrix& matrix,
726 template<
class G,
class S>
770 VertexSet& connectivity, std::vector<Vertex>& front_);
795 void add(std::vector<Vertex>& vertex);
804 typename VertexSet::size_type
size();
851 std::vector<Vertex>& front_;
901 template<
class M,
class C>
902 tuple<int,int,int,int>
build(
const M& m, G& graph,
910 typedef PoolAllocator<Vertex,100> Allocator;
915 typedef SLList<Vertex,Allocator> VertexList;
920 typedef std::set<Vertex,std::less<Vertex>,Allocator> VertexSet;
925 typedef std::size_t* SphereMap;
940 std::vector<Vertex> front_;
945 VertexSet connected_;
966 enum { N = 1300000 };
1008 class AggregateVisitor
1068 class FrontNeighbourCounter :
public Counter
1087 int noFrontNeighbours(
const Vertex& vertex)
const;
1092 class TwoWayCounter :
public Counter
1110 const AggregatesMap<Vertex>& aggregates)
const;
1115 class OneWayCounter :
public Counter
1133 const AggregatesMap<Vertex>& aggregates)
const;
1141 class ConnectivityCounter :
public Counter
1150 ConnectivityCounter(
const VertexSet& connected,
const AggregatesMap<Vertex>& aggregates);
1156 const VertexSet& connected_;
1158 const AggregatesMap<Vertex>& aggregates_;
1173 double connectivity(
const Vertex& vertex,
const AggregatesMap<Vertex>& aggregates)
const;
1182 const AggregatesMap<Vertex>& aggregates)
const;
1191 bool connected(
const Vertex& vertex,
const SLList<AggregateDescriptor>& aggregateList,
1192 const AggregatesMap<Vertex>& aggregates)
const;
1201 class DependencyCounter :
public Counter
1207 DependencyCounter();
1227 FrontMarker(std::vector<Vertex>& front,
MatrixGraph& graph);
1233 std::vector<Vertex>& front_;
1257 int unusedNeighbours(
const Vertex& vertex,
const AggregatesMap<Vertex>& aggregates)
const;
1272 std::pair<int,int> neighbours(
const Vertex& vertex,
1274 const AggregatesMap<Vertex>& aggregates)
const;
1291 int aggregateNeighbours(
const Vertex& vertex,
const AggregateDescriptor& aggregate,
const AggregatesMap<Vertex>& aggregates)
const;
1300 bool admissible(
const Vertex& vertex,
const AggregateDescriptor& aggregate,
const AggregatesMap<Vertex>& aggregates)
const;
1309 std::size_t distance(
const Vertex& vertex,
const AggregatesMap<Vertex>& aggregates);
1319 Vertex mergeNeighbour(
const Vertex& vertex,
const AggregatesMap<Vertex>& aggregates)
const;
1329 void nonisoNeighbourAggregate(
const Vertex& vertex,
1330 const AggregatesMap<Vertex>& aggregates,
1331 SLList<Vertex>& neighbours)
const;
1341 void growAggregate(
const Vertex& vertex,
const AggregatesMap<Vertex>& aggregates,
const C& c);
1343 void growIsolatedAggregate(
const Vertex& vertex,
const AggregatesMap<Vertex>& aggregates,
const C& c);
1348 template<
class M,
class N>
1354 template<
class M,
class N>
1357 DUNE_UNUSED_PARAMETER(row);
1358 maxValue_ = std::min(- std::numeric_limits<typename Matrix::field_type>::max(), std::numeric_limits<typename Matrix::field_type>::min());
1360 diagonal_ = norm_(matrix_->operator[](row_)[row_]);
1363 template<
class M,
class N>
1370 if(!N::is_sign_preserving || eij<0 || eji<0)
1371 maxValue_ = std::max(maxValue_,
1372 eij /diagonal_ * eji/
1373 norm_(matrix_->operator[](col.index())[col.index()]));
1376 template<
class M,
class N>
1384 if(!N::is_sign_preserving || (eij<0 || eji<0))
1385 if(eji / norm_(matrix_->operator[](edge.target())[edge.target()]) *
1386 eij/ diagonal_ > alpha() * maxValue_) {
1387 edge.properties().setDepends();
1388 edge.properties().setInfluences();
1390 typename G::EdgeProperties& other = graph.getEdgeProperties(edge.target(), edge.source());
1391 other.setInfluences();
1396 template<
class M,
class N>
1399 return maxValue_ < beta();
1403 template<
class M,
class N>
1409 template<
class M,
class N>
1412 DUNE_UNUSED_PARAMETER(row);
1413 maxValue_ = std::min(- std::numeric_limits<typename Matrix::field_type>::max(), std::numeric_limits<typename Matrix::field_type>::min());
1415 diagonal_ = norm_(matrix_->operator[](row_)[row_]);
1418 template<
class M,
class N>
1421 maxValue_ = std::max(maxValue_,
1425 template<
class M,
class N>
1429 if(-norm_(*col) >= maxValue_ * alpha()) {
1430 edge.properties().setDepends();
1431 typedef typename G::EdgeDescriptor ED;
1432 ED e= graph.findEdge(edge.target(), edge.source());
1433 if(e!=std::numeric_limits<ED>::max())
1435 typename G::EdgeProperties& other = graph.getEdgeProperties(e);
1436 other.setInfluences();
1441 template<
class M,
class N>
1444 return maxValue_ < beta() * diagonal_;
1447 template<
class G,
class S>
1449 VertexSet& connected, std::vector<Vertex>& front)
1450 : vertices_(), id_(-1), graph_(graph), aggregates_(aggregates),
1451 connected_(connected), front_(front)
1454 template<
class G,
class S>
1462 throw "Not yet implemented";
1470 template<
class G,
class S>
1473 dvverb<<
"Connected cleared"<<std::endl;
1476 connected_.insert(vertex);
1477 dvverb <<
" Inserting "<<vertex<<
" size="<<connected_.size();
1483 template<
class G,
class S>
1486 vertices_.insert(vertex);
1487 aggregates_[vertex]=id_;
1489 front_.erase(std::lower_bound(front_.begin(), front_.end(), vertex));
1493 const iterator end = graph_.endEdges(vertex);
1494 for(iterator edge = graph_.beginEdges(vertex); edge != end; ++edge) {
1495 dvverb <<
" Inserting "<<aggregates_[edge.target()];
1496 connected_.insert(aggregates_[edge.target()]);
1497 dvverb <<
" size="<<connected_.size();
1499 !graph_.getVertexProperties(edge.target()).front())
1501 front_.push_back(edge.target());
1502 graph_.getVertexProperties(edge.target()).setFront();
1506 std::sort(front_.begin(), front_.end());
1509 template<
class G,
class S>
1513 std::size_t oldsize = vertices_.size();
1515 typedef typename std::vector<Vertex>::iterator Iterator;
1517 typedef typename VertexSet::iterator SIterator;
1519 SIterator pos=vertices_.begin();
1520 std::vector<Vertex> newFront;
1521 newFront.reserve(front_.capacity());
1523 std::set_difference(front_.begin(), front_.end(), vertices.begin(), vertices.end(),
1524 std::back_inserter(newFront));
1527 for(Iterator vertex=vertices.begin(); vertex != vertices.end(); ++vertex)
1529 pos=vertices_.insert(pos,*vertex);
1530 vertices_.insert(*vertex);
1531 graph_.getVertexProperties(*vertex).resetFront();
1532 aggregates_[*vertex]=id_;
1535 const iterator end = graph_.endEdges(*vertex);
1536 for(iterator edge = graph_.beginEdges(*vertex); edge != end; ++edge) {
1537 dvverb <<
" Inserting "<<aggregates_[edge.target()];
1538 connected_.insert(aggregates_[edge.target()]);
1540 !graph_.getVertexProperties(edge.target()).front())
1542 front_.push_back(edge.target());
1543 graph_.getVertexProperties(edge.target()).setFront();
1545 dvverb <<
" size="<<connected_.size();
1549 std::sort(front_.begin(), front_.end());
1550 assert(oldsize+vertices.size()==vertices_.size());
1552 template<
class G,
class S>
1560 template<
class G,
class S>
1561 inline typename Aggregate<G,S>::VertexSet::size_type
1564 return vertices_.size();
1567 template<
class G,
class S>
1568 inline typename Aggregate<G,S>::VertexSet::size_type
1571 return connected_.size();
1574 template<
class G,
class S>
1580 template<
class G,
class S>
1583 return vertices_.begin();
1586 template<
class G,
class S>
1589 return vertices_.end();
1593 const V AggregatesMap<V>::UNAGGREGATED = std::numeric_limits<V>::max();
1596 const V AggregatesMap<V>::ISOLATED = std::numeric_limits<V>::max()-1;
1607 delete[] aggregates_;
1614 allocate(noVertices);
1627 noVertices_ = noVertices;
1629 for(std::size_t i=0; i < noVertices; i++)
1630 aggregates_[i]=UNAGGREGATED;
1636 assert(aggregates_ != 0);
1637 delete[] aggregates_;
1642 inline typename AggregatesMap<V>::AggregateDescriptor&
1645 return aggregates_[v];
1649 inline const typename AggregatesMap<V>::AggregateDescriptor&
1652 return aggregates_[v];
1656 template<
bool reset,
class G,
class F,
class VM>
1659 const G& graph, F& aggregateVisitor,
1660 VM& visitedMap)
const
1664 DummyEdgeVisitor dummy;
1665 return breadthFirstSearch<true,reset>(start, aggregate, graph, vlist, aggregateVisitor, dummy, visitedMap);
1669 template<
bool remove,
bool reset,
class G,
class L,
class F1,
class F2,
class VM>
1674 F1& aggregateVisitor,
1675 F2& nonAggregateVisitor,
1676 VM& visitedMap)
const
1678 typedef typename L::const_iterator ListIterator;
1679 int visitedSpheres = 0;
1681 visited.push_back(start);
1682 put(visitedMap, start,
true);
1684 ListIterator current = visited.begin();
1685 ListIterator end = visited.end();
1686 std::size_t i=0, size=visited.size();
1690 while(current != end) {
1692 for(; i<size; ++current, ++i) {
1693 typedef typename G::ConstEdgeIterator EdgeIterator;
1694 const EdgeIterator endEdge = graph.endEdges(*current);
1696 for(EdgeIterator edge = graph.beginEdges(*current);
1697 edge != endEdge; ++edge) {
1699 if(aggregates_[edge.target()]==aggregate) {
1700 if(!
get(visitedMap, edge.target())) {
1701 put(visitedMap, edge.target(),
true);
1702 visited.push_back(edge.target());
1703 aggregateVisitor(edge);
1706 nonAggregateVisitor(edge);
1709 end = visited.end();
1710 size = visited.size();
1716 for(current = visited.begin(); current != end; ++current)
1717 put(visitedMap, *current,
false);
1723 return visitedSpheres;
1728 : graph_(0), aggregate_(0), front_(), connected_(), size_(-1)
1737 template<
class G,
class C>
1739 const typename C::Matrix& matrix,
1740 C criterion,
bool firstlevel)
1743 typedef typename C::Matrix Matrix;
1744 typedef typename G::VertexIterator VertexIterator;
1746 criterion.init(&matrix);
1748 for(VertexIterator vertex = graph.begin(); vertex != graph.end(); ++vertex) {
1751 const Row& row = matrix[*vertex];
1756 criterion.initRow(row, *vertex);
1761 ColIterator end = row.end();
1765 for(ColIterator col = row.begin(); col != end; ++
col)
1766 if(col.index()!=*vertex) {
1767 criterion.examine(col);
1768 absoffdiag=std::max(absoffdiag, col->frobenius_norm());
1772 vertex.properties().setExcludedBorder();
1775 for(ColIterator col = row.begin(); col != end; ++
col)
1776 if(col.index()!=*vertex)
1777 criterion.examine(col);
1783 if(criterion.isIsolated()) {
1785 vertex.properties().setIsolated();
1788 typedef typename G::EdgeIterator EdgeIterator;
1790 EdgeIterator eEnd = vertex.end();
1791 ColIterator col = matrix[*vertex].begin();
1793 for(EdgeIterator edge = vertex.begin(); edge!= eEnd; ++edge, ++
col) {
1795 while(col.index()!=edge.target())
1797 criterion.examine(graph, edge, col);
1807 inline Aggregator<G>::AggregateVisitor<V>::AggregateVisitor(
const AggregatesMap<Vertex>& aggregates,
1809 : aggregates_(aggregates), aggregate_(aggregate), visitor_(&visitor)
1816 if(aggregates_[edge.target()]==aggregate_)
1817 visitor_->operator()(edge);
1822 inline void Aggregator<G>::visitAggregateNeighbours(
const Vertex& vertex,
1824 const AggregatesMap<Vertex>& aggregates,
1828 AggregateVisitor<V> v(aggregates, aggregate, visitor);
1834 inline Aggregator<G>::Counter::Counter()
1839 inline void Aggregator<G>::Counter::increment()
1845 inline void Aggregator<G>::Counter::decrement()
1850 inline int Aggregator<G>::Counter::value()
1858 if(edge.properties().isTwoWay())
1859 Counter::increment();
1864 const AggregatesMap<Vertex>& aggregates)
const
1866 TwoWayCounter counter;
1867 visitAggregateNeighbours(vertex, aggregate, aggregates, counter);
1868 return counter.value();
1873 const AggregatesMap<Vertex>& aggregates)
const
1875 OneWayCounter counter;
1876 visitAggregateNeighbours(vertex, aggregate, aggregates, counter);
1877 return counter.value();
1883 if(edge.properties().isOneWay())
1884 Counter::increment();
1888 inline Aggregator<G>::ConnectivityCounter::ConnectivityCounter(
const VertexSet& connected,
1889 const AggregatesMap<Vertex>& aggregates)
1890 : Counter(), connected_(connected), aggregates_(aggregates)
1899 Counter::increment();
1901 Counter::increment();
1902 Counter::increment();
1907 inline double Aggregator<G>::connectivity(
const Vertex& vertex,
const AggregatesMap<Vertex>& aggregates)
const
1909 ConnectivityCounter counter(connected_, aggregates);
1911 return (
double)counter.value()/noNeighbours;
1915 inline Aggregator<G>::DependencyCounter::DependencyCounter()
1922 if(edge.properties().depends())
1923 Counter::increment();
1924 if(edge.properties().influences())
1925 Counter::increment();
1929 int Aggregator<G>::unusedNeighbours(
const Vertex& vertex,
const AggregatesMap<Vertex>& aggregates)
const
1935 std::pair<int,int> Aggregator<G>::neighbours(
const Vertex& vertex,
1937 const AggregatesMap<Vertex>& aggregates)
const
1939 DependencyCounter unused, aggregated;
1940 typedef AggregateVisitor<DependencyCounter> Counter;
1941 typedef tuple<Counter,Counter> CounterTuple;
1944 return std::make_pair(unused.value(), aggregated.value());
1949 int Aggregator<G>::aggregateNeighbours(
const Vertex& vertex,
const AggregateDescriptor& aggregate,
const AggregatesMap<Vertex>& aggregates)
const
1951 DependencyCounter counter;
1952 visitAggregateNeighbours(vertex, aggregate, aggregates, counter);
1953 return counter.value();
1957 std::size_t Aggregator<G>::distance(
const Vertex& vertex,
const AggregatesMap<Vertex>& aggregates)
1960 typename PropertyMapTypeSelector<VertexVisitedTag,G>::Type visitedMap =
get(VertexVisitedTag(), *graph_);
1962 typename AggregatesMap<Vertex>::DummyEdgeVisitor dummy;
1963 return aggregates.template breadthFirstSearch<true,true>(vertex,
1964 aggregate_->
id(), *graph_,
1965 vlist, dummy, dummy, visitedMap);
1969 inline Aggregator<G>::FrontMarker::FrontMarker(std::vector<Vertex>& front,
MatrixGraph& graph)
1970 : front_(front), graph_(graph)
1976 Vertex target = edge.target();
1978 if(!graph_.getVertexProperties(target).front()) {
1979 front_.push_back(target);
1980 graph_.getVertexProperties(target).setFront();
1985 inline bool Aggregator<G>::admissible(
const Vertex& vertex,
const AggregateDescriptor& aggregate,
const AggregatesMap<Vertex>& aggregates)
const
1988 Dune::dvverb<<
" Admissible not yet implemented!"<<std::endl;
1995 Iterator vend = graph_->endEdges(vertex);
1996 for(Iterator edge = graph_->beginEdges(vertex); edge != vend; ++edge) {
1998 if(edge.properties().isStrong()
1999 && aggregates[edge.target()]==aggregate)
2002 Iterator edge1 = edge;
2003 for(++edge1; edge1 != vend; ++edge1) {
2005 if(edge1.properties().isStrong()
2006 && aggregates[edge.target()]==aggregate)
2011 Iterator v2end = graph_->endEdges(edge.target());
2012 for(Iterator edge2 = graph_->beginEdges(edge.target()); edge2 != v2end; ++edge2) {
2013 if(edge2.target()==edge1.target() &&
2014 edge2.properties().isStrong()) {
2030 vend = graph_->endEdges(vertex);
2031 for(Iterator edge = graph_->beginEdges(vertex); edge != vend; ++edge) {
2033 if(edge.properties().isStrong()
2034 && aggregates[edge.target()]==aggregate)
2037 Iterator v1end = graph_->endEdges(edge.target());
2039 for(Iterator edge1=graph_->beginEdges(edge.target()); edge1 != v1end; ++edge1) {
2041 if(edge1.properties().isStrong()
2042 && aggregates[edge1.target()]==aggregate)
2046 Iterator v2end = graph_->endEdges(vertex);
2047 for(Iterator edge2 = graph_->beginEdges(vertex); edge2 != v2end; ++edge2) {
2048 if(edge2.target()==edge1.target()) {
2049 if(edge2.properties().isStrong())
2066 void Aggregator<G>::unmarkFront()
2068 typedef typename std::vector<Vertex>::const_iterator Iterator;
2070 for(Iterator vertex=front_.begin(); vertex != front_.end(); ++vertex)
2071 graph_->getVertexProperties(*vertex).resetFront();
2078 Aggregator<G>::nonisoNeighbourAggregate(
const Vertex& vertex,
2079 const AggregatesMap<Vertex>& aggregates,
2080 SLList<Vertex>& neighbours)
const
2083 Iterator end=graph_->beginEdges(vertex);
2086 for(Iterator edge=graph_->beginEdges(vertex); edge!=end; ++edge)
2089 neighbours.push_back(aggregates[edge.target()]);
2094 inline typename G::VertexDescriptor Aggregator<G>::mergeNeighbour(
const Vertex& vertex,
const AggregatesMap<Vertex>& aggregates)
const
2098 Iterator end = graph_->endEdges(vertex);
2099 for(Iterator edge = graph_->beginEdges(vertex); edge != end; ++edge) {
2101 graph_->getVertexProperties(edge.target()).isolated() == graph_->getVertexProperties(edge.source()).isolated()) {
2102 if( graph_->getVertexProperties(vertex).isolated() ||
2103 ((edge.properties().depends() || edge.properties().influences())
2104 && admissible(vertex, aggregates[edge.target()], aggregates)))
2105 return edge.target();
2112 Aggregator<G>::FrontNeighbourCounter::FrontNeighbourCounter(
const MatrixGraph& graph)
2113 : Counter(), graph_(graph)
2119 if(graph_.getVertexProperties(edge.target()).front())
2120 Counter::increment();
2124 int Aggregator<G>::noFrontNeighbours(
const Vertex& vertex)
const
2126 FrontNeighbourCounter counter(*graph_);
2128 return counter.value();
2131 inline bool Aggregator<G>::connected(
const Vertex& vertex,
2133 const AggregatesMap<Vertex>& aggregates)
const
2135 typedef typename G::ConstEdgeIterator iterator;
2136 const iterator end = graph_->endEdges(vertex);
2137 for(iterator edge = graph_->beginEdges(vertex); edge != end; ++edge)
2138 if(aggregates[edge.target()]==aggregate)
2143 inline bool Aggregator<G>::connected(
const Vertex& vertex,
2144 const SLList<AggregateDescriptor>& aggregateList,
2145 const AggregatesMap<Vertex>& aggregates)
const
2147 typedef typename SLList<AggregateDescriptor>::const_iterator Iter;
2148 for(Iter i=aggregateList.begin(); i!=aggregateList.end(); ++i)
2149 if(connected(vertex, *i, aggregates))
2156 void Aggregator<G>::growIsolatedAggregate(
const Vertex& seed,
const AggregatesMap<Vertex>& aggregates,
const C& c)
2158 SLList<Vertex> connectedAggregates;
2159 nonisoNeighbourAggregate(seed, aggregates,connectedAggregates);
2161 while(aggregate_->
size()< c.minAggregateSize() && aggregate_->
connectSize() < c.maxConnectivity()) {
2163 std::size_t maxFrontNeighbours=0;
2167 typedef typename std::vector<Vertex>::const_iterator Iterator;
2169 for(Iterator vertex = front_.begin(); vertex != front_.end(); ++vertex) {
2170 if(distance(*vertex, aggregates)>c.maxDistance())
2173 if(connectedAggregates.size()>0) {
2177 if(!connected(*vertex, connectedAggregates, aggregates))
2181 double con = connectivity(*vertex, aggregates);
2184 std::size_t frontNeighbours = noFrontNeighbours(*vertex);
2186 if(frontNeighbours >= maxFrontNeighbours) {
2187 maxFrontNeighbours = frontNeighbours;
2188 candidate = *vertex;
2190 }
else if(con > maxCon) {
2192 maxFrontNeighbours = noFrontNeighbours(*vertex);
2193 candidate = *vertex;
2200 aggregate_->
add(candidate);
2206 void Aggregator<G>::growAggregate(
const Vertex& seed,
const AggregatesMap<Vertex>& aggregates,
const C& c)
2209 std::size_t distance_ =0;
2210 while(aggregate_->
size() < c.minAggregateSize()&& distance_<c.maxDistance()) {
2211 int maxTwoCons=0, maxOneCons=0, maxNeighbours=-1;
2214 std::vector<Vertex> candidates;
2215 candidates.reserve(30);
2217 typedef typename std::vector<Vertex>::const_iterator Iterator;
2219 for(Iterator vertex = front_.begin(); vertex != front_.end(); ++vertex) {
2221 if(graph_->getVertexProperties(*vertex).isolated())
2224 int twoWayCons = twoWayConnections(*vertex, aggregate_->
id(), aggregates);
2227 if( maxTwoCons == twoWayCons && twoWayCons > 0) {
2228 double con = connectivity(*vertex, aggregates);
2231 int neighbours = noFrontNeighbours(*vertex);
2233 if(neighbours > maxNeighbours) {
2234 maxNeighbours = neighbours;
2236 candidates.push_back(*vertex);
2238 candidates.push_back(*vertex);
2240 }
else if( con > maxCon) {
2242 maxNeighbours = noFrontNeighbours(*vertex);
2244 candidates.push_back(*vertex);
2246 }
else if(twoWayCons > maxTwoCons) {
2247 maxTwoCons = twoWayCons;
2248 maxCon = connectivity(*vertex, aggregates);
2249 maxNeighbours = noFrontNeighbours(*vertex);
2251 candidates.push_back(*vertex);
2254 maxOneCons = std::numeric_limits<int>::max();
2263 int oneWayCons = oneWayConnections(*vertex, aggregate_->
id(), aggregates);
2268 if(!admissible(*vertex, aggregate_->
id(), aggregates))
2271 if( maxOneCons == oneWayCons && oneWayCons > 0) {
2272 double con = connectivity(*vertex, aggregates);
2275 int neighbours = noFrontNeighbours(*vertex);
2277 if(neighbours > maxNeighbours) {
2278 maxNeighbours = neighbours;
2280 candidates.push_back(*vertex);
2282 if(neighbours==maxNeighbours)
2284 candidates.push_back(*vertex);
2287 }
else if( con > maxCon) {
2289 maxNeighbours = noFrontNeighbours(*vertex);
2291 candidates.push_back(*vertex);
2293 }
else if(oneWayCons > maxOneCons) {
2294 maxOneCons = oneWayCons;
2295 maxCon = connectivity(*vertex, aggregates);
2296 maxNeighbours = noFrontNeighbours(*vertex);
2298 candidates.push_back(*vertex);
2303 if(!candidates.size())
2305 distance_=distance(seed, aggregates);
2306 candidates.resize(std::min(candidates.size(), c.maxAggregateSize()-
2307 aggregate_->
size()));
2308 aggregate_->
add(candidates);
2312 template<
typename V>
2313 template<
typename M,
typename G,
typename C>
2317 Aggregator<G> aggregator;
2318 return aggregator.build(matrix, graph, *
this, criterion, finestLevel);
2322 template<
class M,
class C>
2323 tuple<int,int,int,int>
Aggregator<G>::build(
const M& m, G& graph, AggregatesMap<Vertex>& aggregates,
const C& c,
2327 Stack stack_(graph, *
this, aggregates);
2331 aggregate_ =
new Aggregate<G,VertexSet>(graph, aggregates, connected_, front_);
2338 dverb<<
"Build dependency took "<< watch.elapsed()<<
" seconds."<<std::endl;
2339 int noAggregates, conAggregates, isoAggregates, oneAggregates;
2340 std::size_t maxA=0, minA=1000000, avg=0;
2341 int skippedAggregates;
2342 noAggregates = conAggregates = isoAggregates = oneAggregates =
2343 skippedAggregates = 0;
2346 Vertex seed = stack_.pop();
2348 if(seed == Stack::NullEntry)
2353 if((noAggregates+1)%10000 == 0)
2357 if(graph.getVertexProperties(seed).excludedBorder()) {
2359 ++skippedAggregates;
2363 if(graph.getVertexProperties(seed).isolated()) {
2364 if(c.skipIsolated()) {
2367 ++skippedAggregates;
2371 aggregate_->
seed(seed);
2372 growIsolatedAggregate(seed, aggregates, c);
2375 aggregate_->
seed(seed);
2376 growAggregate(seed, aggregates, c);
2380 while(!(graph.getVertexProperties(seed).isolated()) && aggregate_->
size() < c.maxAggregateSize()) {
2382 std::vector<Vertex> candidates;
2383 candidates.reserve(30);
2385 typedef typename std::vector<Vertex>::const_iterator Iterator;
2387 for(Iterator vertex = front_.begin(); vertex != front_.end(); ++vertex) {
2389 if(graph.getVertexProperties(*vertex).isolated())
2392 if(twoWayConnections( *vertex, aggregate_->
id(), aggregates) == 0 &&
2393 (oneWayConnections( *vertex, aggregate_->
id(), aggregates) == 0 ||
2394 !admissible( *vertex, aggregate_->
id(), aggregates) ))
2397 std::pair<int,int> neighbourPair=neighbours(*vertex, aggregate_->
id(),
2403 if(neighbourPair.first >= neighbourPair.second)
2406 if(distance(*vertex, aggregates) > c.maxDistance())
2408 candidates.push_back(*vertex);
2412 if(!candidates.size())
break;
2414 candidates.resize(std::min(candidates.size(), c.maxAggregateSize()-
2415 aggregate_->
size()));
2416 aggregate_->
add(candidates);
2421 if(aggregate_->
size()==1 && c.maxAggregateSize()>1) {
2422 if(!graph.getVertexProperties(seed).isolated()) {
2423 Vertex mergedNeighbour = mergeNeighbour(seed, aggregates);
2427 aggregates[seed] = aggregates[mergedNeighbour];
2431 minA=std::min(minA,static_cast<std::size_t>(1));
2432 maxA=std::max(maxA,static_cast<std::size_t>(1));
2438 minA=std::min(minA,static_cast<std::size_t>(1));
2439 maxA=std::max(maxA,static_cast<std::size_t>(1));
2445 avg+=aggregate_->
size();
2446 minA=std::min(minA,aggregate_->
size());
2447 maxA=std::max(maxA,aggregate_->
size());
2448 if(graph.getVertexProperties(seed).isolated())
2456 Dune::dinfo<<
"connected aggregates: "<<conAggregates;
2457 Dune::dinfo<<
" isolated aggregates: "<<isoAggregates;
2458 if(conAggregates+isoAggregates>0)
2459 Dune::dinfo<<
" one node aggregates: "<<oneAggregates<<
" min size="
2460 <<minA<<
" max size="<<maxA
2461 <<
" avg="<<avg/(conAggregates+isoAggregates)<<std::endl;
2464 return make_tuple(conAggregates+isoAggregates,isoAggregates,
2465 oneAggregates,skippedAggregates);
2470 Aggregator<G>::Stack::Stack(
const MatrixGraph& graph,
const Aggregator<G>& aggregatesBuilder,
2471 const AggregatesMap<Vertex>& aggregates)
2472 : graph_(graph), aggregatesBuilder_(aggregatesBuilder), aggregates_(aggregates), begin_(graph.begin()), end_(graph.end())
2478 Aggregator<G>::Stack::~Stack()
2485 const typename Aggregator<G>::Vertex Aggregator<G>::Stack::NullEntry
2486 = std::numeric_limits<typename G::VertexDescriptor>::max();
2489 inline typename G::VertexDescriptor Aggregator<G>::Stack::pop()
2495 typename G::VertexDescriptor current=*begin_;
2507 std::ios_base::fmtflags oldOpts=os.flags();
2509 os.setf(std::ios_base::right, std::ios_base::adjustfield);
2514 for(
int i=0; i< n*m; i++)
2515 max=std::max(max, aggregates[i]);
2517 for(
int i=10; i < 1000000; i*=10)
2523 for(
int j=0, entry=0; j < m; j++) {
2524 for(
int i=0; i<n; i++, entry++) {
2526 os<<aggregates[entry]<<
" ";
std::size_t * SphereMap
Type of the mapping of aggregate members onto distance spheres.
Definition: aggregates.hh:759
const_iterator end() const
Definition: aggregates.hh:672
Definition: aggregates.hh:411
row_type::const_iterator ConstColIterator
Const iterator for the entries of each row.
Definition: matrix.hh:53
Matrix::row_type Row
Constant Row iterator of the matrix.
Definition: aggregates.hh:146
const_iterator end() const
get an iterator over the vertices of the aggregate.
A Dummy visitor that does nothing for each visited edge.
Definition: aggregates.hh:535
const AggregateDescriptor * const_iterator
Definition: aggregates.hh:665
void buildDependency(G &graph, const typename C::Matrix &matrix, C criterion, bool finestLevel)
Build the dependency of the matrix graph.
Matrix::field_type maxValue_
The current max value.
Definition: aggregates.hh:350
M::field_type operator()(const M &m) const
compute the norm of a matrix.
Definition: aggregates.hh:376
Definition: aggregates.hh:431
int visitNeighbours(const G &graph, const typename G::VertexDescriptor &vertex, V &visitor)
Visit all neighbour vertices of a vertex in a graph.
void init(const Matrix *matrix)
Dependency()
Definition: aggregates.hh:283
MatrixGraph::VertexDescriptor Vertex
The vertex descriptor type.
Definition: aggregates.hh:739
void initRow(const Row &row, int index)
G MatrixGraph
Definition: aggregates.hh:735
M::size_type VertexDescriptor
The vertex descriptor.
Definition: graph.hh:72
const Matrix * matrix_
The matrix we work on.
Definition: aggregates.hh:289
void initRow(const Row &row, int index)
Iterator over all edges starting from a vertex.
Definition: graph.hh:93
VertexSet::size_type connectSize()
Get tne number of connections to other aggregates.
Matrix::field_type diagonal_
The norm of the current diagonal.
Definition: aggregates.hh:297
void invalidate()
Definition: aggregates.hh:772
static const V ISOLATED
Identifier of isolated vertices.
Definition: aggregates.hh:509
Criterion taking advantage of symmetric matrices.
Definition: aggregates.hh:453
void examine(const ColIter &col)
SLList< VertexDescriptor, Allocator > VertexList
The type of a single linked list of vertex descriptors.
Definition: aggregates.hh:530
void printAggregates2d(const AggregatesMap< V > &aggregates, int n, int m, std::ostream &os)
Definition: aggregates.hh:2505
const Matrix * matrix_
The matrix we work on.
Definition: aggregates.hh:348
AggregateDescriptor * iterator
Definition: aggregates.hh:677
M Matrix
The matrix type we build the dependency of.
Definition: aggregates.hh:136
T::field_type field_type
Export the type representing the underlying field.
Definition: matrix.hh:29
Class for building the aggregates.
Definition: aggregates.hh:486
SymmetricCriterion(const Parameters &parms)
Definition: aggregates.hh:456
Norm that uses only the [0][0] entry of the block to determine couplings.
Definition: aggregates.hh:386
tuple< int, int, int, int > build(const M &m, G &graph, AggregatesMap< Vertex > &aggregates, const C &c, bool finestLevel)
Build the aggregates.
SymmetricCriterion()
Definition: aggregates.hh:459
V VertexDescriptor
The vertex descriptor type.
Definition: aggregates.hh:513
Aggregate(MatrixGraph &graph, AggregatesMap< Vertex > &aggregates, VertexSet &connectivity, std::vector< Vertex > &front_)
Constructor.
Base class of all aggregation criterions.
Definition: aggregates.hh:45
The (undirected) graph of a matrix.
Definition: graph.hh:49
Functor using the row sum (infinity) norm to determine strong couplings.
Definition: aggregates.hh:394
VertexSet::size_type size()
Get the size of the aggregate.
SymmetricDependency()
Definition: aggregates.hh:342
UnSymmetricCriterion(const Parameters &parms)
Definition: aggregates.hh:479
void init(const Matrix *matrix)
void setDefaultValuesAnisotropic(std::size_t dim, std::size_t diameter=2)
Sets reasonable default values for an aisotropic problem.
Definition: aggregates.hh:102
void clear()
Clear the aggregate.
void init(const Matrix *matrix)
Definition: aggregates.hh:190
void free()
Free the allocated memory.
Col col
Definition: matrixmatrix.hh:347
V Visitor
The type of the adapted visitor.
Definition: aggregates.hh:1014
int row_
index of the currently evaluated row.
Definition: aggregates.hh:180
Matrix::row_type Row
Constant Row iterator of the matrix.
Definition: aggregates.hh:261
Class providing information about the mapping of the vertices onto aggregates.
Definition: aggregates.hh:497
Matrix::ConstColIterator ColIter
Constant column iterator of the matrix.
Definition: aggregates.hh:151
M::field_type operator()(const M &m) const
compute the norm of a matrix.
Definition: aggregates.hh:405
SymmetricDependency(const Parameters &parms)
Definition: aggregates.hh:339
std::vector< typename Matrix::field_type > vals_
Definition: aggregates.hh:183
N Norm
The norm to use for examining the matrix entries.
Definition: aggregates.hh:256
AggregationCriterion(const Parameters &parms)
Definition: aggregates.hh:67
static const Vertex NullEntry
Definition: aggregates.hh:958
Matrix::row_type Row
Constant Row iterator of the matrix.
Definition: aggregates.hh:320
Provides classes for building the matrix graph.
SymmetricMatrixDependency()
Definition: aggregates.hh:168
SymmetricMatrixDependency(const Parameters &parms)
Definition: aggregates.hh:165
G MatrixGraph
The matrix graph type used.
Definition: aggregates.hh:865
MatrixGraph::VertexDescriptor Vertex
The vertex identifier.
Definition: aggregates.hh:870
Dependency policy for symmetric matrices.
Definition: aggregates.hh:130
Dependency policy for symmetric matrices.
Definition: aggregates.hh:245
Definition: aggregates.hh:427
Definition: aggregates.hh:368
V AggregateDescriptor
The aggregate descriptor type.
Definition: aggregates.hh:518
const Matrix * matrix_
The matrix we work on.
Definition: aggregates.hh:174
PoolAllocator< VertexDescriptor, 100 > Allocator
The allocator we use for our lists and the set.
Definition: aggregates.hh:524
void reconstruct(const Vertex &vertex)
Reconstruct the aggregat from an seed node.
MatrixGraph::VertexDescriptor AggregateDescriptor
The type of the aggregate descriptor.
Definition: aggregates.hh:873
N Norm
The norm to use for examining the matrix entries.
Definition: aggregates.hh:141
Matrix::field_type diagonal_
The norm of the current diagonal.
Definition: aggregates.hh:356
Criterion suited for unsymmetric matrices.
Definition: aggregates.hh:476
void seed(const Vertex &vertex)
Initialize the aggregate with one vertex.
PoolAllocator< Vertex, 100 > Allocator
The allocator we use for our lists and the set.
Definition: aggregates.hh:745
T DependencyPolicy
The policy for calculating the dependency graph.
Definition: aggregates.hh:52
M::field_type operator()(const M &m) const
compute the norm of a matrix.
Definition: aggregates.hh:422
void add(const Vertex &vertex)
Add a vertex to the aggregate.
int id()
Get the id identifying the aggregate.
Norm norm_
The functor for calculating the norm.
Definition: aggregates.hh:293
N Norm
The norm to use for examining the matrix entries.
Definition: aggregates.hh:315
The vertex iterator type of the graph.
Definition: graph.hh:207
AggregateDescriptor & operator[](const VertexDescriptor &v)
Get the aggregate a vertex belongs to.
Row row
Definition: matrixmatrix.hh:345
AggregatesMap()
Constructs without allocating memory.
A class for temporarily storing the vertices of an aggregate in.
Definition: aggregates.hh:727
void initRow(const Row &row, int index)
Definition: aggregates.hh:196
tuple< int, int, int, int > buildAggregates(const M &matrix, G &graph, const C &criterion, bool finestLevel)
Build the aggregates.
VariableBlockVector< T, A >::window_type row_type
The type implementing a matrix row.
Definition: matrix.hh:38
Parameter classes for customizing AMG.
~AggregatesMap()
Destructor.
Dependency(const Parameters &parms)
Definition: aggregates.hh:279
S VertexSet
The type of a single linked list of vertex descriptors.
Definition: aggregates.hh:751
std::size_t index
Definition: matrixmarket.hh:561
derive error class from the base class in common
Definition: istlexception.hh:16
Provides classes for handling internal properties in a graph.
M Matrix
The matrix type we build the dependency of.
Definition: aggregates.hh:310
int row_
index of the currently evaluated row.
Definition: aggregates.hh:354
Matrix::field_type maxValue_
The current max value.
Definition: aggregates.hh:176
M Matrix
The matrix type we build the dependency of.
Definition: aggregates.hh:251
void allocate(std::size_t noVertices)
Allocate memory for holding the information.
Matrix::field_type maxValue_
The current max value.
Definition: aggregates.hh:291
bool isIsolated()
Definition: aggregates.hh:233
M::field_type operator()(const M &m) const
compute the norm of a matrix.
Definition: aggregates.hh:438
Definition: aggregates.hh:415
Definition: aggregates.hh:398
int row_
index of the currently evaluated row.
Definition: aggregates.hh:295
std::vector< typename Matrix::field_type >::iterator valIter_
Definition: aggregates.hh:184
iterator end()
Definition: aggregates.hh:684
iterator begin()
Definition: aggregates.hh:679
const_iterator begin() const
Definition: aggregates.hh:667
Norm norm_
The functor for calculating the norm.
Definition: aggregates.hh:352
static const V UNAGGREGATED
Identifier of not yet aggregated vertices.
Definition: aggregates.hh:504
std::size_t noVertices() const
Get the number of vertices.
Definition: bvector.hh:584
AggregationCriterion()
Constructor.
Definition: aggregates.hh:63
VertexSet::const_iterator const_iterator
Const iterator over a vertex list.
Definition: aggregates.hh:754
void examine(const ColIter &col)
Definition: aggregates.hh:208
const_iterator begin() const
get an iterator over the vertices of the aggregate.
void examine(const ColIter &col)
Dependency policy for symmetric matrices.
Definition: aggregates.hh:304
std::size_t breadthFirstSearch(const VertexDescriptor &start, const AggregateDescriptor &aggregate, const G &graph, F &aggregateVisitor, VM &visitedMap) const
Breadth first search within an aggregate.
EdgeIteratorT< const MatrixGraph< Matrix > > ConstEdgeIterator
The constant edge iterator type.
Definition: graph.hh:297
friend class Stack
Definition: aggregates.hh:986
void operator()(const EdgeIterator &edge) const
Definition: aggregates.hh:539
void setDefaultValuesIsotropic(std::size_t dim, std::size_t diameter=2)
Sets reasonable default values for an isotropic problem.
Definition: aggregates.hh:79
Matrix::ConstColIterator ColIter
Constant column iterator of the matrix.
Definition: aggregates.hh:325
UnSymmetricCriterion()
Definition: aggregates.hh:482
Norm norm_
The functor for calculating the norm.
Definition: aggregates.hh:178
All parameters for AMG.
Definition: parameters.hh:390
Matrix::ConstColIterator ColIter
Constant column iterator of the matrix.
Definition: aggregates.hh:266
Definition: example.cc:34
Matrix::field_type diagonal_
The norm of the current diagonal.
Definition: aggregates.hh:182
Norm that uses only the [N][N] entry of the block to determine couplings.
Definition: aggregates.hh:364