3#ifndef DUNE_AMG_AGGREGATES_HH
4#define DUNE_AMG_AGGREGATES_HH
10#include "combinedfunctor.hh"
82 this->setMaxDistance(diameter-1);
87 this->setMaxDistance(this->maxDistance()+diameter-1);
89 this->setMinAggregateSize(csize);
90 this->setMaxAggregateSize(
static_cast<std::size_t
>(csize*1.5));
106 this->setMaxDistance(this->maxDistance()+dim-1);
111 std::ostream& operator<<(std::ostream& os,
const AggregationCriterion<T>& criterion)
113 os<<
"{ maxdistance="<<criterion.maxDistance()<<
" minAggregateSize="
114 <<criterion.minAggregateSize()<<
" maxAggregateSize="<<criterion.maxAggregateSize()
115 <<
" connectivity="<<criterion.maxConnectivity()<<
" debugLevel="<<criterion.debugLevel()<<
"}";
130 template<
class M,
class N>
154 void init(
const Matrix* matrix);
156 void initRow(
const Row& row,
int index);
158 void examine(
const ColIter& col);
161 void examine(G& graph,
const typename G::EdgeIterator& edge,
const ColIter& col);
178 typedef typename FieldTraits<field_type>::real_type real_type;
186 std::vector<real_type> vals_;
187 typename std::vector<real_type>::iterator valIter_;
192 template<
class M,
class N>
198 template<
class M,
class N>
199 inline void SymmetricMatrixDependency<M,N>::initRow(
const Row& row,
int index)
201 vals_.assign(row.size(), 0.0);
202 assert(vals_.size()==row.size());
203 valIter_=vals_.begin();
205 maxValue_ = std::min(- std::numeric_limits<real_type>::max(), std::numeric_limits<real_type>::min());
206 diagonal_=norm_(row[index]);
210 template<
class M,
class N>
211 inline void SymmetricMatrixDependency<M,N>::examine(
const ColIter& col)
214 real_type eij = norm_(*col);
215 if(!N::is_sign_preserving || eij<0)
217 *valIter_ = eij/diagonal_*eij/norm_(matrix_->operator[](col.index())[col.index()]);
218 maxValue_ = std::max(maxValue_, *valIter_);
224 template<
class M,
class N>
226 inline void SymmetricMatrixDependency<M,N>::examine(G&,
const typename G::EdgeIterator& edge,
const ColIter&)
228 if(*valIter_ > alpha() * maxValue_) {
229 edge.properties().setDepends();
230 edge.properties().setInfluences();
235 template<
class M,
class N>
236 inline bool SymmetricMatrixDependency<M,N>::isIsolated()
240 valIter_=vals_.begin();
241 return maxValue_ < beta();
247 template<
class M,
class N>
271 void init(
const Matrix* matrix);
273 void initRow(
const Row& row,
int index);
275 void examine(
const ColIter& col);
278 void examine(G& graph,
const typename G::EdgeIterator& edge,
const ColIter& col);
295 typedef typename FieldTraits<field_type>::real_type real_type;
308 template<
class M,
class N>
332 void init(
const Matrix* matrix);
334 void initRow(
const Row& row,
int index);
336 void examine(
const ColIter& col);
339 void examine(G& graph,
const typename G::EdgeIterator& edge,
const ColIter& col);
356 typedef typename FieldTraits<field_type>::real_type real_type;
375 is_sign_preserving =
true
383 typename FieldTraits<typename M::field_type>::real_type
operator()(
const M& m)
const
385 typedef typename M::field_type field_type;
386 typedef typename FieldTraits<field_type>::real_type real_type;
387 static_assert( std::is_convertible<field_type, real_type >::value,
388 "use of diagonal norm in AMG not implemented for complex field_type");
397 static T signed_abs(
const T & v)
404 static T signed_abs(
const std::complex<T> & v)
408 return csgn(v) * std::abs(v);
413 static T csgn(
const T & v)
415 return (T(0) < v) - (v < T(0));
420 static T csgn(std::complex<T> a)
422 return csgn(a.real())+(a.real() == 0.0)*csgn(a.imag());
443 is_sign_preserving =
false
450 typename FieldTraits<typename M::field_type>::real_type
operator()(
const M& m)
const
452 return m.infinity_norm();
460 is_sign_preserving =
false
467 typename FieldTraits<typename M::field_type>::real_type
operator()(
const M& m)
const
469 return m.frobenius_norm();
476 is_sign_preserving =
false
483 typename FieldTraits<typename M::field_type>::real_type
operator()(
const M& m)
const
497 template<
class M,
class Norm>
520 template<
class M,
class Norm>
583 template<
class EdgeIterator>
584 void operator()(
const EdgeIterator& edge)
const
619 template<
class M,
class G,
class C>
620 std::tuple<int,int,int,int>
buildAggregates(
const M& matrix, G& graph,
const C& criterion,
642 template<
bool reset,
class G,
class F,
class VM>
647 VM& visitedMap)
const;
672 template<
bool remove,
bool reset,
class G,
class L,
class F1,
class F2,
class VM>
675 const G& graph, L& visited, F1& aggregateVisitor,
676 F2& nonAggregateVisitor,
677 VM& visitedMap)
const;
712 const_iterator begin()
const
717 const_iterator end()
const
737 AggregatesMap<V>& operator=(
const AggregatesMap<V>&) =
delete;
747 std::size_t noVertices_;
753 template<
class G,
class C>
755 const typename C::Matrix& matrix,
763 template<
class G,
class S>
772 typedef G MatrixGraph;
807 VertexSet& connectivity, std::vector<Vertex>& front_);
841 typename VertexSet::size_type
size();
888 std::vector<Vertex>& front_;
938 template<
class M,
class C>
939 std::tuple<int,int,int,int>
build(
const M& m, G& graph,
957 typedef std::set<Vertex,std::less<Vertex>,
Allocator> VertexSet;
962 typedef std::size_t* SphereMap;
977 std::vector<Vertex> front_;
982 VertexSet connected_;
995 static const Vertex NullEntry;
1003 enum { N = 1300000 };
1037 const AggregatesMap<Vertex>& aggregates,
1045 class AggregateVisitor
1105 class FrontNeighbourCounter :
public Counter
1129 class TwoWayCounter :
public Counter
1152 class OneWayCounter :
public Counter
1170 const AggregatesMap<Vertex>& aggregates)
const;
1178 class ConnectivityCounter :
public Counter
1193 const VertexSet& connected_;
1238 class DependencyCounter :
public Counter
1270 std::vector<Vertex>& front_;
1385 template<
class M,
class N>
1391 template<
class M,
class N>
1392 inline void SymmetricDependency<M,N>::initRow(
const Row& row,
int index)
1395 maxValue_ = std::min(- std::numeric_limits<typename Matrix::field_type>::max(), std::numeric_limits<typename Matrix::field_type>::min());
1397 diagonal_ = norm_(matrix_->operator[](row_)[row_]);
1400 template<
class M,
class N>
1401 inline void SymmetricDependency<M,N>::examine(
const ColIter& col)
1403 real_type eij = norm_(*col);
1405 matrix_->operator[](col.index()).find(row_);
1406 if ( opposite_entry == matrix_->operator[](col.index()).end() )
1411 real_type eji = norm_(*opposite_entry);
1414 if(!N::is_sign_preserving || eij<0 || eji<0)
1415 maxValue_ = std::max(maxValue_,
1416 eij /diagonal_ * eji/
1417 norm_(matrix_->operator[](col.index())[col.index()]));
1420 template<
class M,
class N>
1422 inline void SymmetricDependency<M,N>::examine(G& graph,
const typename G::EdgeIterator& edge,
const ColIter& col)
1424 real_type eij = norm_(*col);
1426 matrix_->operator[](col.index()).find(row_);
1428 if ( opposite_entry == matrix_->operator[](col.index()).end() )
1433 real_type eji = norm_(*opposite_entry);
1435 if(!N::is_sign_preserving || (eij<0 || eji<0))
1436 if(eji / norm_(matrix_->operator[](edge.target())[edge.target()]) *
1437 eij/ diagonal_ > alpha() * maxValue_) {
1438 edge.properties().setDepends();
1439 edge.properties().setInfluences();
1440 typename G::EdgeProperties& other = graph.getEdgeProperties(edge.target(), edge.source());
1441 other.setInfluences();
1446 template<
class M,
class N>
1447 inline bool SymmetricDependency<M,N>::isIsolated()
1449 return maxValue_ < beta();
1453 template<
class M,
class N>
1454 inline void Dependency<M,N>::init(
const Matrix* matrix)
1459 template<
class M,
class N>
1460 inline void Dependency<M,N>::initRow(
const Row& row,
int index)
1463 maxValue_ = std::min(- std::numeric_limits<real_type>::max(), std::numeric_limits<real_type>::min());
1465 diagonal_ = norm_(matrix_->operator[](row_)[row_]);
1468 template<
class M,
class N>
1469 inline void Dependency<M,N>::examine(
const ColIter& col)
1471 maxValue_ = std::max(maxValue_,
1475 template<
class M,
class N>
1477 inline void Dependency<M,N>::examine(G& graph,
const typename G::EdgeIterator& edge,
const ColIter& col)
1479 if(-norm_(*col) >= maxValue_ * alpha()) {
1480 edge.properties().setDepends();
1481 typedef typename G::EdgeDescriptor ED;
1482 ED e= graph.findEdge(edge.target(), edge.source());
1483 if(e!=std::numeric_limits<ED>::max())
1485 typename G::EdgeProperties& other = graph.getEdgeProperties(e);
1486 other.setInfluences();
1491 template<
class M,
class N>
1492 inline bool Dependency<M,N>::isIsolated()
1494 return maxValue_ < beta() * diagonal_;
1497 template<
class G,
class S>
1499 VertexSet& connected, std::vector<Vertex>& front)
1500 : vertices_(), id_(-1), graph_(graph), aggregates_(aggregates),
1501 connected_(connected), front_(
front)
1504 template<
class G,
class S>
1512 throw "Not yet implemented";
1520 template<
class G,
class S>
1523 dvverb<<
"Connected cleared"<<std::endl;
1526 connected_.insert(
vertex);
1527 dvverb <<
" Inserting "<<
vertex<<
" size="<<connected_.size();
1533 template<
class G,
class S>
1536 vertices_.insert(
vertex);
1539 front_.erase(std::lower_bound(front_.begin(), front_.end(),
vertex));
1543 const iterator end = graph_.endEdges(
vertex);
1544 for(iterator edge = graph_.beginEdges(
vertex); edge != end; ++edge) {
1545 dvverb <<
" Inserting "<<aggregates_[edge.target()];
1546 connected_.insert(aggregates_[edge.target()]);
1547 dvverb <<
" size="<<connected_.size();
1549 !graph_.getVertexProperties(edge.target()).front())
1551 front_.push_back(edge.target());
1552 graph_.getVertexProperties(edge.target()).setFront();
1556 std::sort(front_.begin(), front_.end());
1559 template<
class G,
class S>
1563 std::size_t oldsize = vertices_.size();
1565 typedef typename std::vector<Vertex>::iterator Iterator;
1567 typedef typename VertexSet::iterator SIterator;
1569 SIterator pos=vertices_.begin();
1570 std::vector<Vertex> newFront;
1571 newFront.reserve(front_.capacity());
1573 std::set_difference(front_.begin(), front_.end(), vertices.begin(), vertices.end(),
1574 std::back_inserter(newFront));
1579 pos=vertices_.insert(pos,*
vertex);
1580 vertices_.insert(*
vertex);
1581 graph_.getVertexProperties(*vertex).resetFront();
1582 aggregates_[*
vertex]=id_;
1585 const iterator end = graph_.endEdges(*
vertex);
1586 for(iterator edge = graph_.beginEdges(*
vertex); edge != end; ++edge) {
1587 dvverb <<
" Inserting "<<aggregates_[edge.target()];
1588 connected_.insert(aggregates_[edge.target()]);
1590 !graph_.getVertexProperties(edge.target()).front())
1592 front_.push_back(edge.target());
1593 graph_.getVertexProperties(edge.target()).setFront();
1595 dvverb <<
" size="<<connected_.size();
1599 std::sort(front_.begin(), front_.end());
1600 assert(oldsize+vertices.size()==vertices_.size());
1602 template<
class G,
class S>
1610 template<
class G,
class S>
1611 inline typename Aggregate<G,S>::VertexSet::size_type
1614 return vertices_.size();
1617 template<
class G,
class S>
1618 inline typename Aggregate<G,S>::VertexSet::size_type
1621 return connected_.size();
1624 template<
class G,
class S>
1630 template<
class G,
class S>
1633 return vertices_.begin();
1636 template<
class G,
class S>
1639 return vertices_.end();
1657 delete[] aggregates_;
1664 allocate(noVertices);
1677 noVertices_ = noVertices;
1679 for(std::size_t i=0; i < noVertices; i++)
1680 aggregates_[i]=UNAGGREGATED;
1686 assert(aggregates_ != 0);
1687 delete[] aggregates_;
1695 return aggregates_[v];
1702 return aggregates_[v];
1706 template<
bool reset,
class G,
class F,
class VM>
1709 const G& graph, F& aggregateVisitor,
1710 VM& visitedMap)
const
1714 DummyEdgeVisitor dummy;
1715 return breadthFirstSearch<true,reset>(start, aggregate, graph, vlist, aggregateVisitor, dummy, visitedMap);
1719 template<
bool remove,
bool reset,
class G,
class L,
class F1,
class F2,
class VM>
1724 F1& aggregateVisitor,
1725 F2& nonAggregateVisitor,
1726 VM& visitedMap)
const
1728 typedef typename L::const_iterator ListIterator;
1729 int visitedSpheres = 0;
1731 visited.push_back(start);
1732 put(visitedMap, start,
true);
1734 ListIterator current = visited.begin();
1735 ListIterator end = visited.end();
1736 std::size_t i=0, size=visited.size();
1740 while(current != end) {
1742 for(; i<size; ++current, ++i) {
1743 typedef typename G::ConstEdgeIterator EdgeIterator;
1744 const EdgeIterator endEdge = graph.endEdges(*current);
1746 for(EdgeIterator edge = graph.beginEdges(*current);
1747 edge != endEdge; ++edge) {
1749 if(aggregates_[edge.target()]==aggregate) {
1750 if(!get(visitedMap, edge.target())) {
1751 put(visitedMap, edge.target(),
true);
1752 visited.push_back(edge.target());
1753 aggregateVisitor(edge);
1756 nonAggregateVisitor(edge);
1759 end = visited.end();
1760 size = visited.size();
1766 for(current = visited.begin(); current != end; ++current)
1767 put(visitedMap, *current,
false);
1773 return visitedSpheres;
1778 : graph_(0), aggregate_(0), front_(), connected_(), size_(-1)
1787 template<
class G,
class C>
1789 const typename C::Matrix& matrix,
1790 C criterion,
bool firstlevel)
1793 typedef typename C::Matrix Matrix;
1794 typedef typename G::VertexIterator VertexIterator;
1796 criterion.init(&matrix);
1801 const Row& row = matrix[*
vertex];
1806 criterion.initRow(row, *
vertex);
1811 ColIterator end = row.end();
1812 typename FieldTraits<typename Matrix::field_type>::real_type absoffdiag=0.;
1815 for(ColIterator col = row.begin(); col != end; ++col)
1816 if(col.index()!=*
vertex) {
1817 criterion.examine(col);
1818 absoffdiag=std::max(absoffdiag, col->frobenius_norm());
1822 vertex.properties().setExcludedBorder();
1825 for(ColIterator col = row.begin(); col != end; ++col)
1827 criterion.examine(col);
1833 if(criterion.isIsolated()) {
1835 vertex.properties().setIsolated();
1838 typedef typename G::EdgeIterator EdgeIterator;
1840 EdgeIterator eEnd =
vertex.end();
1841 ColIterator col = matrix[*
vertex].begin();
1843 for(EdgeIterator edge =
vertex.begin(); edge!= eEnd; ++edge, ++col) {
1845 while(col.index()!=edge.target())
1847 criterion.examine(graph, edge, col);
1857 inline Aggregator<G>::AggregateVisitor<V>::AggregateVisitor(
const AggregatesMap<Vertex>& aggregates,
1859 : aggregates_(aggregates), aggregate_(aggregate), visitor_(&visitor)
1866 if(aggregates_[edge.target()]==aggregate_)
1867 visitor_->operator()(edge);
1872 inline void Aggregator<G>::visitAggregateNeighbours(
const Vertex&
vertex,
1874 const AggregatesMap<Vertex>& aggregates,
1878 AggregateVisitor<V> v(aggregates, aggregate, visitor);
1884 inline Aggregator<G>::Counter::Counter()
1889 inline void Aggregator<G>::Counter::increment()
1895 inline void Aggregator<G>::Counter::decrement()
1900 inline int Aggregator<G>::Counter::value()
1908 if(edge.properties().isTwoWay())
1909 Counter::increment();
1914 const AggregatesMap<Vertex>& aggregates)
const
1916 TwoWayCounter counter;
1917 visitAggregateNeighbours(
vertex, aggregate, aggregates, counter);
1918 return counter.value();
1923 const AggregatesMap<Vertex>& aggregates)
const
1925 OneWayCounter counter;
1926 visitAggregateNeighbours(
vertex, aggregate, aggregates, counter);
1927 return counter.value();
1933 if(edge.properties().isOneWay())
1934 Counter::increment();
1938 inline Aggregator<G>::ConnectivityCounter::ConnectivityCounter(
const VertexSet& connected,
1939 const AggregatesMap<Vertex>& aggregates)
1940 : Counter(), connected_(connected), aggregates_(aggregates)
1949 Counter::increment();
1951 Counter::increment();
1952 Counter::increment();
1957 inline double Aggregator<G>::connectivity(
const Vertex&
vertex,
const AggregatesMap<Vertex>& aggregates)
const
1959 ConnectivityCounter counter(connected_, aggregates);
1961 return (
double)counter.value()/noNeighbours;
1965 inline Aggregator<G>::DependencyCounter::DependencyCounter()
1972 if(edge.properties().depends())
1973 Counter::increment();
1974 if(edge.properties().influences())
1975 Counter::increment();
1979 int Aggregator<G>::unusedNeighbours(
const Vertex&
vertex,
const AggregatesMap<Vertex>& aggregates)
const
1985 std::pair<int,int> Aggregator<G>::neighbours(
const Vertex&
vertex,
1987 const AggregatesMap<Vertex>& aggregates)
const
1989 DependencyCounter unused, aggregated;
1990 typedef AggregateVisitor<DependencyCounter> Counter;
1991 typedef std::tuple<Counter,Counter> CounterTuple;
1994 return std::make_pair(unused.value(), aggregated.value());
2001 DependencyCounter counter;
2002 visitAggregateNeighbours(
vertex, aggregate, aggregates, counter);
2003 return counter.value();
2007 std::size_t Aggregator<G>::distance(
const Vertex&
vertex,
const AggregatesMap<Vertex>& aggregates)
2010 typename PropertyMapTypeSelector<VertexVisitedTag,G>::Type visitedMap = get(VertexVisitedTag(), *graph_);
2012 typename AggregatesMap<Vertex>::DummyEdgeVisitor dummy;
2013 return aggregates.template breadthFirstSearch<true,true>(
vertex,
2014 aggregate_->
id(), *graph_,
2015 vlist, dummy, dummy, visitedMap);
2019 inline Aggregator<G>::FrontMarker::FrontMarker(std::vector<Vertex>& front,
MatrixGraph& graph)
2020 : front_(
front), graph_(graph)
2026 Vertex target = edge.target();
2028 if(!graph_.getVertexProperties(target).front()) {
2029 front_.push_back(target);
2030 graph_.getVertexProperties(target).setFront();
2038 Dune::dvverb<<
" Admissible not yet implemented!"<<std::endl;
2045 Iterator vend = graph_->endEdges(
vertex);
2046 for(Iterator edge = graph_->beginEdges(
vertex); edge != vend; ++edge) {
2048 if(edge.properties().isStrong()
2049 && aggregates[edge.target()]==aggregate)
2052 Iterator edge1 = edge;
2053 for(++edge1; edge1 != vend; ++edge1) {
2055 if(edge1.properties().isStrong()
2056 && aggregates[edge.target()]==aggregate)
2061 Iterator v2end = graph_->endEdges(edge.target());
2062 for(Iterator edge2 = graph_->beginEdges(edge.target()); edge2 != v2end; ++edge2) {
2063 if(edge2.target()==edge1.target() &&
2064 edge2.properties().isStrong()) {
2080 vend = graph_->endEdges(
vertex);
2081 for(Iterator edge = graph_->beginEdges(
vertex); edge != vend; ++edge) {
2083 if(edge.properties().isStrong()
2084 && aggregates[edge.target()]==aggregate)
2087 Iterator v1end = graph_->endEdges(edge.target());
2089 for(Iterator edge1=graph_->beginEdges(edge.target()); edge1 != v1end; ++edge1) {
2091 if(edge1.properties().isStrong()
2092 && aggregates[edge1.target()]==aggregate)
2096 Iterator v2end = graph_->endEdges(
vertex);
2097 for(Iterator edge2 = graph_->beginEdges(
vertex); edge2 != v2end; ++edge2) {
2098 if(edge2.target()==edge1.target()) {
2099 if(edge2.properties().isStrong())
2116 void Aggregator<G>::unmarkFront()
2118 typedef typename std::vector<Vertex>::const_iterator Iterator;
2121 graph_->getVertexProperties(*vertex).resetFront();
2128 Aggregator<G>::nonisoNeighbourAggregate(
const Vertex&
vertex,
2129 const AggregatesMap<Vertex>& aggregates,
2130 SLList<Vertex>& neighbours)
const
2133 Iterator end=graph_->beginEdges(
vertex);
2136 for(Iterator edge=graph_->beginEdges(
vertex); edge!=end; ++edge)
2139 neighbours.push_back(aggregates[edge.target()]);
2144 inline typename G::VertexDescriptor Aggregator<G>::mergeNeighbour(
const Vertex&
vertex,
const AggregatesMap<Vertex>& aggregates)
const
2148 Iterator end = graph_->endEdges(
vertex);
2149 for(Iterator edge = graph_->beginEdges(
vertex); edge != end; ++edge) {
2151 graph_->getVertexProperties(edge.target()).isolated() == graph_->getVertexProperties(edge.source()).isolated()) {
2152 if( graph_->getVertexProperties(
vertex).isolated() ||
2153 ((edge.properties().depends() || edge.properties().influences())
2154 && admissible(
vertex, aggregates[edge.target()], aggregates)))
2155 return edge.target();
2162 Aggregator<G>::FrontNeighbourCounter::FrontNeighbourCounter(
const MatrixGraph& graph)
2163 : Counter(), graph_(graph)
2169 if(graph_.getVertexProperties(edge.target()).front())
2170 Counter::increment();
2174 int Aggregator<G>::noFrontNeighbours(
const Vertex&
vertex)
const
2176 FrontNeighbourCounter counter(*graph_);
2178 return counter.value();
2181 inline bool Aggregator<G>::connected(
const Vertex&
vertex,
2183 const AggregatesMap<Vertex>& aggregates)
const
2185 typedef typename G::ConstEdgeIterator iterator;
2186 const iterator end = graph_->endEdges(
vertex);
2187 for(iterator edge = graph_->beginEdges(
vertex); edge != end; ++edge)
2188 if(aggregates[edge.target()]==aggregate)
2193 inline bool Aggregator<G>::connected(
const Vertex&
vertex,
2194 const SLList<AggregateDescriptor>& aggregateList,
2195 const AggregatesMap<Vertex>& aggregates)
const
2198 for(Iter i=aggregateList.begin(); i!=aggregateList.end(); ++i)
2199 if(connected(
vertex, *i, aggregates))
2206 void Aggregator<G>::growIsolatedAggregate(
const Vertex& seed,
const AggregatesMap<Vertex>& aggregates,
const C& c)
2208 SLList<Vertex> connectedAggregates;
2209 nonisoNeighbourAggregate(seed, aggregates,connectedAggregates);
2211 while(aggregate_->
size()< c.minAggregateSize() && aggregate_->
connectSize() < c.maxConnectivity()) {
2213 std::size_t maxFrontNeighbours=0;
2217 typedef typename std::vector<Vertex>::const_iterator Iterator;
2220 if(distance(*
vertex, aggregates)>c.maxDistance())
2223 if(connectedAggregates.size()>0) {
2227 if(!connected(*
vertex, connectedAggregates, aggregates))
2231 double con = connectivity(*
vertex, aggregates);
2234 std::size_t frontNeighbours = noFrontNeighbours(*
vertex);
2236 if(frontNeighbours >= maxFrontNeighbours) {
2237 maxFrontNeighbours = frontNeighbours;
2240 }
else if(con > maxCon) {
2242 maxFrontNeighbours = noFrontNeighbours(*
vertex);
2250 aggregate_->
add(candidate);
2256 void Aggregator<G>::growAggregate(
const Vertex& seed,
const AggregatesMap<Vertex>& aggregates,
const C& c)
2259 std::size_t distance_ =0;
2260 while(aggregate_->
size() < c.minAggregateSize()&& distance_<c.maxDistance()) {
2261 int maxTwoCons=0, maxOneCons=0, maxNeighbours=-1;
2264 std::vector<Vertex> candidates;
2265 candidates.reserve(30);
2267 typedef typename std::vector<Vertex>::const_iterator Iterator;
2271 if(graph_->getVertexProperties(*vertex).isolated())
2274 int twoWayCons = twoWayConnections(*
vertex, aggregate_->
id(), aggregates);
2277 if( maxTwoCons == twoWayCons && twoWayCons > 0) {
2278 double con = connectivity(*
vertex, aggregates);
2281 int neighbours = noFrontNeighbours(*
vertex);
2283 if(neighbours > maxNeighbours) {
2284 maxNeighbours = neighbours;
2286 candidates.push_back(*
vertex);
2288 candidates.push_back(*
vertex);
2290 }
else if( con > maxCon) {
2292 maxNeighbours = noFrontNeighbours(*
vertex);
2294 candidates.push_back(*
vertex);
2296 }
else if(twoWayCons > maxTwoCons) {
2297 maxTwoCons = twoWayCons;
2298 maxCon = connectivity(*
vertex, aggregates);
2299 maxNeighbours = noFrontNeighbours(*
vertex);
2301 candidates.push_back(*
vertex);
2304 maxOneCons = std::numeric_limits<int>::max();
2313 int oneWayCons = oneWayConnections(*
vertex, aggregate_->
id(), aggregates);
2318 if(!admissible(*
vertex, aggregate_->
id(), aggregates))
2321 if( maxOneCons == oneWayCons && oneWayCons > 0) {
2322 double con = connectivity(*
vertex, aggregates);
2325 int neighbours = noFrontNeighbours(*
vertex);
2327 if(neighbours > maxNeighbours) {
2328 maxNeighbours = neighbours;
2330 candidates.push_back(*
vertex);
2332 if(neighbours==maxNeighbours)
2334 candidates.push_back(*
vertex);
2337 }
else if( con > maxCon) {
2339 maxNeighbours = noFrontNeighbours(*
vertex);
2341 candidates.push_back(*
vertex);
2343 }
else if(oneWayCons > maxOneCons) {
2344 maxOneCons = oneWayCons;
2345 maxCon = connectivity(*
vertex, aggregates);
2346 maxNeighbours = noFrontNeighbours(*
vertex);
2348 candidates.push_back(*
vertex);
2353 if(!candidates.size())
2355 distance_=distance(seed, aggregates);
2356 candidates.resize(std::min(candidates.size(), c.maxAggregateSize()-
2357 aggregate_->
size()));
2358 aggregate_->
add(candidates);
2362 template<
typename V>
2363 template<
typename M,
typename G,
typename C>
2367 Aggregator<G> aggregator;
2368 return aggregator.build(matrix, graph, *
this, criterion, finestLevel);
2372 template<
class M,
class C>
2373 std::tuple<int,int,int,int>
Aggregator<G>::build(
const M& m, G& graph, AggregatesMap<Vertex>& aggregates,
const C& c,
2377 Stack stack_(graph, *
this, aggregates);
2381 aggregate_ =
new Aggregate<G,VertexSet>(graph, aggregates, connected_, front_);
2388 dverb<<
"Build dependency took "<< watch.elapsed()<<
" seconds."<<std::endl;
2389 int noAggregates, conAggregates, isoAggregates, oneAggregates;
2390 std::size_t maxA=0, minA=1000000, avg=0;
2391 int skippedAggregates;
2392 noAggregates = conAggregates = isoAggregates = oneAggregates =
2393 skippedAggregates = 0;
2396 Vertex seed = stack_.pop();
2398 if(seed == Stack::NullEntry)
2403 if((noAggregates+1)%10000 == 0)
2407 if(graph.getVertexProperties(seed).excludedBorder()) {
2409 ++skippedAggregates;
2413 if(graph.getVertexProperties(seed).isolated()) {
2414 if(c.skipIsolated()) {
2417 ++skippedAggregates;
2421 aggregate_->
seed(seed);
2422 growIsolatedAggregate(seed, aggregates, c);
2425 aggregate_->
seed(seed);
2426 growAggregate(seed, aggregates, c);
2430 while(!(graph.getVertexProperties(seed).isolated()) && aggregate_->
size() < c.maxAggregateSize()) {
2432 std::vector<Vertex> candidates;
2433 candidates.reserve(30);
2435 typedef typename std::vector<Vertex>::const_iterator Iterator;
2439 if(graph.getVertexProperties(*vertex).isolated())
2442 if(twoWayConnections( *
vertex, aggregate_->
id(), aggregates) == 0 &&
2443 (oneWayConnections( *
vertex, aggregate_->
id(), aggregates) == 0 ||
2444 !admissible( *
vertex, aggregate_->
id(), aggregates) ))
2447 std::pair<int,int> neighbourPair=neighbours(*
vertex, aggregate_->
id(),
2453 if(neighbourPair.first >= neighbourPair.second)
2456 if(distance(*
vertex, aggregates) > c.maxDistance())
2458 candidates.push_back(*
vertex);
2462 if(!candidates.size())
break;
2464 candidates.resize(std::min(candidates.size(), c.maxAggregateSize()-
2465 aggregate_->
size()));
2466 aggregate_->
add(candidates);
2471 if(aggregate_->
size()==1 && c.maxAggregateSize()>1) {
2472 if(!graph.getVertexProperties(seed).isolated()) {
2473 Vertex mergedNeighbour = mergeNeighbour(seed, aggregates);
2477 aggregates[seed] = aggregates[mergedNeighbour];
2478 aggregate_->invalidate();
2481 minA=std::min(minA,
static_cast<std::size_t
>(1));
2482 maxA=std::max(maxA,
static_cast<std::size_t
>(1));
2488 minA=std::min(minA,
static_cast<std::size_t
>(1));
2489 maxA=std::max(maxA,
static_cast<std::size_t
>(1));
2495 avg+=aggregate_->
size();
2496 minA=std::min(minA,aggregate_->
size());
2497 maxA=std::max(maxA,aggregate_->
size());
2498 if(graph.getVertexProperties(seed).isolated())
2506 Dune::dinfo<<
"connected aggregates: "<<conAggregates;
2507 Dune::dinfo<<
" isolated aggregates: "<<isoAggregates;
2508 if(conAggregates+isoAggregates>0)
2509 Dune::dinfo<<
" one node aggregates: "<<oneAggregates<<
" min size="
2510 <<minA<<
" max size="<<maxA
2511 <<
" avg="<<avg/(conAggregates+isoAggregates)<<std::endl;
2514 return std::make_tuple(conAggregates+isoAggregates,isoAggregates,
2515 oneAggregates,skippedAggregates);
2520 Aggregator<G>::Stack::Stack(
const MatrixGraph& graph,
const Aggregator<G>& aggregatesBuilder,
2521 const AggregatesMap<Vertex>& aggregates)
2522 : graph_(graph), aggregatesBuilder_(aggregatesBuilder), aggregates_(aggregates), begin_(graph.begin()), end_(graph.end())
2528 Aggregator<G>::Stack::~Stack()
2536 = std::numeric_limits<typename G::VertexDescriptor>::max();
2539 inline typename G::VertexDescriptor Aggregator<G>::Stack::pop()
2545 typename G::VertexDescriptor current=*begin_;
2555 void printAggregates2d(
const AggregatesMap<V>& aggregates,
int n,
int m, std::ostream& os)
2557 std::ios_base::fmtflags oldOpts=os.flags();
2559 os.setf(std::ios_base::right, std::ios_base::adjustfield);
2564 for(
int i=0; i< n*m; i++)
2565 max=std::max(max, aggregates[i]);
2567 for(
int i=10; i < 1000000; i*=10)
2573 for(
int j=0, entry=0; j < m; j++) {
2574 for(
int i=0; i<n; i++, entry++) {
2576 os<<aggregates[entry]<<
" ";
A class for temporarily storing the vertices of an aggregate in.
Definition: aggregates.hh:765
A Dummy visitor that does nothing for each visited edge.
Definition: aggregates.hh:581
Class providing information about the mapping of the vertices onto aggregates.
Definition: aggregates.hh:543
Base class of all aggregation criterions.
Definition: aggregates.hh:47
Class for building the aggregates.
Definition: aggregates.hh:896
Dependency policy for symmetric matrices.
Definition: aggregates.hh:249
Norm that uses only the [N][N] entry of the block to determine couplings.
Definition: aggregates.hh:372
Norm that uses only the [0][0] entry of the block to determine couplings.
Definition: aggregates.hh:432
Iterator over all edges starting from a vertex.
Definition: graph.hh:93
The vertex iterator type of the graph.
Definition: graph.hh:207
The (undirected) graph of a matrix.
Definition: graph.hh:49
M::size_type VertexDescriptor
The vertex descriptor.
Definition: graph.hh:71
EdgeIteratorT< const MatrixGraph< Matrix > > ConstEdgeIterator
The constant edge iterator type.
Definition: graph.hh:296
All parameters for AMG.
Definition: parameters.hh:391
Criterion taking advantage of symmetric matrices.
Definition: aggregates.hh:499
Dependency policy for symmetric matrices.
Definition: aggregates.hh:310
Dependency policy for symmetric matrices.
Definition: aggregates.hh:132
Criterion suited for unsymmetric matrices.
Definition: aggregates.hh:522
derive error class from the base class in common
Definition: istlexception.hh:16
A generic dynamic dense matrix.
Definition: matrix.hh:555
T::field_type field_type
Export the type representing the underlying field.
Definition: matrix.hh:559
row_type::const_iterator ConstColIterator
Const iterator for the entries of each row.
Definition: matrix.hh:583
MatrixImp::DenseMatrixBase< T, A >::window_type row_type
The type implementing a matrix row.
Definition: matrix.hh:568
An allocator managing a pool of objects for reuse.
Definition: poolallocator.hh:247
A single linked list.
Definition: sllist.hh:42
Type traits to determine the type of reals (when working with complex numbers)
Provides classes for building the matrix graph.
SLListConstIterator< T, A > const_iterator
The constant iterator of the list.
Definition: sllist.hh:72
#define DUNE_UNUSED_PARAMETER(parm)
A macro to mark intentionally unused function parameters with.
Definition: unused.hh:25
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
constexpr GeometryType vertex
GeometryType representing a vertex.
Definition: type.hh:727
Matrix::ConstColIterator ColIter
Constant column iterator of the matrix.
Definition: aggregates.hh:269
Matrix::ConstColIterator ColIter
Constant column iterator of the matrix.
Definition: aggregates.hh:152
std::size_t breadthFirstSearch(const VertexDescriptor &start, const AggregateDescriptor &aggregate, const G &graph, L &visited, F1 &aggregateVisitor, F2 &nonAggregateVisitor, VM &visitedMap) const
Breadth first search within an aggregate.
PoolAllocator< VertexDescriptor, 100 > Allocator
The allocator we use for our lists and the set.
Definition: aggregates.hh:569
int id()
Get the id identifying the aggregate.
Norm norm_
The functor for calculating the norm.
Definition: aggregates.hh:298
FieldTraits< typenameM::field_type >::real_type operator()(const M &m) const
compute the norm of a matrix.
Definition: aggregates.hh:383
MatrixGraph::VertexDescriptor Vertex
The vertex identifier.
Definition: aggregates.hh:907
AggregationCriterion()
Constructor.
Definition: aggregates.hh:64
const Matrix * matrix_
The matrix we work on.
Definition: aggregates.hh:353
M Matrix
The matrix type we build the dependency of.
Definition: aggregates.hh:254
G MatrixGraph
The matrix graph type used.
Definition: aggregates.hh:902
Norm norm_
The functor for calculating the norm.
Definition: aggregates.hh:359
M Matrix
The matrix type we build the dependency of.
Definition: aggregates.hh:315
real_type diagonal_
The norm of the current diagonal.
Definition: aggregates.hh:185
N Norm
The norm to use for examining the matrix entries.
Definition: aggregates.hh:259
int row_
index of the currently evaluated row.
Definition: aggregates.hh:183
std::tuple< int, int, int, int > build(const M &m, G &graph, AggregatesMap< Vertex > &aggregates, const C &c, bool finestLevel)
Build the aggregates.
FrontNeighbourCounter(const MatrixGraph &front)
Constructor.
Matrix::row_type Row
Constant Row iterator of the matrix.
Definition: aggregates.hh:325
const Matrix * matrix_
The matrix we work on.
Definition: aggregates.hh:292
const AggregateDescriptor & operator[](const VertexDescriptor &v) const
Get the aggregate a vertex belongs to.
AggregateVisitor(const AggregatesMap< Vertex > &aggregates, const AggregateDescriptor &aggregate, Visitor &visitor)
Constructor.
Matrix::ConstColIterator ColIter
Constant column iterator of the matrix.
Definition: aggregates.hh:330
~AggregatesMap()
Destructor.
Matrix::field_type field_type
The current max value.
Definition: aggregates.hh:177
void decrement()
Decrement counter.
Aggregate(MatrixGraph &graph, AggregatesMap< Vertex > &aggregates, VertexSet &connectivity, std::vector< Vertex > &front_)
Constructor.
V Visitor
The type of the adapted visitor.
Definition: aggregates.hh:1051
std::size_t * SphereMap
Type of the mapping of aggregate members onto distance spheres.
Definition: aggregates.hh:796
Matrix::row_type Row
Constant Row iterator of the matrix.
Definition: aggregates.hh:264
VertexSet::size_type connectSize()
Get tne number of connections to other aggregates.
N Norm
The norm to use for examining the matrix entries.
Definition: aggregates.hh:320
Norm norm_
The functor for calculating the norm.
Definition: aggregates.hh:181
VertexSet::const_iterator const_iterator
Const iterator over a vertex list.
Definition: aggregates.hh:791
MatrixGraph::VertexDescriptor AggregateDescriptor
The type of the aggregate descriptor.
Definition: aggregates.hh:910
real_type diagonal_
The norm of the current diagonal.
Definition: aggregates.hh:363
void add(const Vertex &vertex)
Add a vertex to the aggregate.
T DependencyPolicy
The policy for calculating the dependency graph.
Definition: aggregates.hh:53
void operator()(const typename MatrixGraph::ConstEdgeIterator &edge)
Examine an edge.
FrontMarker(std::vector< Vertex > &front, MatrixGraph &graph)
Constructor.
FieldTraits< typenameM::field_type >::real_type operator()(const M &m) const
compute the norm of a matrix.
Definition: aggregates.hh:483
int visitNeighbours(const G &graph, const typename G::VertexDescriptor &vertex, V &visitor)
Visit all neighbour vertices of a vertex in a graph.
void setDefaultValuesIsotropic(std::size_t dim, std::size_t diameter=2)
Sets reasonable default values for an isotropic problem.
Definition: aggregates.hh:80
V AggregateDescriptor
The aggregate descriptor type.
Definition: aggregates.hh:563
static const V ISOLATED
Identifier of isolated vertices.
Definition: aggregates.hh:554
int row_
index of the currently evaluated row.
Definition: aggregates.hh:361
DependencyCounter()
Constructor.
real_type diagonal_
The norm of the current diagonal.
Definition: aggregates.hh:302
std::size_t noVertices() const
Get the number of vertices.
void setDefaultValuesAnisotropic(std::size_t dim, std::size_t diameter=2)
Sets reasonable default values for an aisotropic problem.
Definition: aggregates.hh:103
AggregatesMap(std::size_t noVertices)
Constructs with allocating memory.
Matrix::field_type field_type
The current max value.
Definition: aggregates.hh:355
AggregateDescriptor & operator[](const VertexDescriptor &v)
Get the aggregate a vertex belongs to.
AggregatesMap()
Constructs without allocating memory.
int value()
Access the current count.
SLList< VertexDescriptor, Allocator > VertexList
The type of a single linked list of vertex descriptors.
Definition: aggregates.hh:575
ConnectivityCounter(const VertexSet &connected, const AggregatesMap< Vertex > &aggregates)
Constructor.
VertexSet::size_type size()
Get the size of the aggregate.
const_iterator end() const
get an iterator over the vertices of the aggregate.
int row_
index of the currently evaluated row.
Definition: aggregates.hh:300
M Matrix
The matrix type we build the dependency of.
Definition: aggregates.hh:137
const Matrix * matrix_
The matrix we work on.
Definition: aggregates.hh:175
S VertexSet
The type of a single linked list of vertex descriptors.
Definition: aggregates.hh:788
FieldTraits< typenameM::field_type >::real_type operator()(const M &m) const
compute the norm of a matrix.
Definition: aggregates.hh:467
static const V UNAGGREGATED
Identifier of not yet aggregated vertices.
Definition: aggregates.hh:549
std::size_t breadthFirstSearch(const VertexDescriptor &start, const AggregateDescriptor &aggregate, const G &graph, F &aggregateVisitor, VM &visitedMap) const
Breadth first search within an aggregate.
Matrix::field_type field_type
The current max value.
Definition: aggregates.hh:294
void allocate(std::size_t noVertices)
Allocate memory for holding the information.
N Norm
The norm to use for examining the matrix entries.
Definition: aggregates.hh:142
FieldTraits< typenameM::field_type >::real_type operator()(const M &m) const
compute the norm of a matrix.
Definition: aggregates.hh:450
void reconstruct(const Vertex &vertex)
Reconstruct the aggregat from an seed node.
const_iterator begin() const
get an iterator over the vertices of the aggregate.
MatrixGraph::VertexDescriptor Vertex
The vertex descriptor type.
Definition: aggregates.hh:776
void seed(const Vertex &vertex)
Initialize the aggregate with one vertex.
void clear()
Clear the aggregate.
void free()
Free the allocated memory.
void increment()
Increment counter.
void buildDependency(G &graph, const typename C::Matrix &matrix, C criterion, bool finestLevel)
Build the dependency of the matrix graph.
V VertexDescriptor
The vertex descriptor type.
Definition: aggregates.hh:558
std::tuple< int, int, int, int > buildAggregates(const M &matrix, G &graph, const C &criterion, bool finestLevel)
Build the aggregates.
Matrix::row_type Row
Constant Row iterator of the matrix.
Definition: aggregates.hh:147
PoolAllocator< Vertex, 100 > Allocator
The allocator we use for our lists and the set.
Definition: aggregates.hh:782
DVVerbType dvverb(std::cout)
stream for very verbose output.
Definition: stdstreams.hh:93
DInfoType dinfo(std::cout)
Stream for informative output.
Definition: stdstreams.hh:138
DVerbType dverb(std::cout)
Singleton of verbose debug stream.
Definition: stdstreams.hh:114
constexpr Front front
PartitionSet for the front partition.
Definition: partitionset.hh:279
Dune namespace.
Definition: alignedallocator.hh:10
Parameter classes for customizing AMG.
An stl-compliant pool allocator.
Provides classes for handling internal properties in a graph.
Implements a singly linked list together with the necessary iterators.
Standard Dune debug streams.
Functor using the row sum (infinity) norm to determine strong couplings.
Definition: aggregates.hh:440
Definition of the DUNE_UNUSED macro for the case that config.h is not available.