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
741 AggregatesMap<V>& operator=(
const AggregatesMap<V>& map)
755 std::size_t noVertices_;
761 template<
class G,
class C>
763 const typename C::Matrix& matrix,
771 template<
class G,
class S>
780 typedef G MatrixGraph;
815 VertexSet& connectivity, std::vector<Vertex>& front_);
840 void add(std::vector<Vertex>& vertex);
849 typename VertexSet::size_type
size();
896 std::vector<Vertex>& front_;
946 template<
class M,
class C>
947 std::tuple<int,int,int,int>
build(
const M& m, G& graph,
965 typedef std::set<Vertex,std::less<Vertex>,
Allocator> VertexSet;
970 typedef std::size_t* SphereMap;
985 std::vector<Vertex> front_;
990 VertexSet connected_;
1003 static const Vertex NullEntry;
1011 enum { N = 1300000 };
1045 const AggregatesMap<Vertex>& aggregates,
1053 class AggregateVisitor
1113 class FrontNeighbourCounter :
public Counter
1132 int noFrontNeighbours(
const Vertex& vertex)
const;
1137 class TwoWayCounter :
public Counter
1160 class OneWayCounter :
public Counter
1178 const AggregatesMap<Vertex>& aggregates)
const;
1186 class ConnectivityCounter :
public Counter
1201 const VertexSet& connected_;
1246 class DependencyCounter :
public Counter
1278 std::vector<Vertex>& front_;
1317 std::pair<int,int> neighbours(
const Vertex& vertex,
1374 void nonisoNeighbourAggregate(
const Vertex& vertex,
1393 template<
class M,
class N>
1399 template<
class M,
class N>
1400 inline void SymmetricDependency<M,N>::initRow(
const Row& row,
int index)
1403 maxValue_ = std::min(- std::numeric_limits<typename Matrix::field_type>::max(), std::numeric_limits<typename Matrix::field_type>::min());
1405 diagonal_ = norm_(matrix_->operator[](row_)[row_]);
1408 template<
class M,
class N>
1409 inline void SymmetricDependency<M,N>::examine(
const ColIter& col)
1411 real_type eij = norm_(*col);
1413 matrix_->operator[](col.index()).find(row_);
1414 if ( opposite_entry == matrix_->operator[](col.index()).end() )
1419 real_type eji = norm_(*opposite_entry);
1422 if(!N::is_sign_preserving || eij<0 || eji<0)
1423 maxValue_ = std::max(maxValue_,
1424 eij /diagonal_ * eji/
1425 norm_(matrix_->operator[](col.index())[col.index()]));
1428 template<
class M,
class N>
1430 inline void SymmetricDependency<M,N>::examine(G& graph,
const typename G::EdgeIterator& edge,
const ColIter& col)
1432 real_type eij = norm_(*col);
1434 matrix_->operator[](col.index()).find(row_);
1436 if ( opposite_entry == matrix_->operator[](col.index()).end() )
1441 real_type eji = norm_(*opposite_entry);
1443 if(!N::is_sign_preserving || (eij<0 || eji<0))
1444 if(eji / norm_(matrix_->operator[](edge.target())[edge.target()]) *
1445 eij/ diagonal_ > alpha() * maxValue_) {
1446 edge.properties().setDepends();
1447 edge.properties().setInfluences();
1448 typename G::EdgeProperties& other = graph.getEdgeProperties(edge.target(), edge.source());
1449 other.setInfluences();
1454 template<
class M,
class N>
1455 inline bool SymmetricDependency<M,N>::isIsolated()
1457 return maxValue_ < beta();
1461 template<
class M,
class N>
1462 inline void Dependency<M,N>::init(
const Matrix* matrix)
1467 template<
class M,
class N>
1468 inline void Dependency<M,N>::initRow(
const Row& row,
int index)
1471 maxValue_ = std::min(- std::numeric_limits<real_type>::max(), std::numeric_limits<real_type>::min());
1473 diagonal_ = norm_(matrix_->operator[](row_)[row_]);
1476 template<
class M,
class N>
1477 inline void Dependency<M,N>::examine(
const ColIter& col)
1479 maxValue_ = std::max(maxValue_,
1483 template<
class M,
class N>
1485 inline void Dependency<M,N>::examine(G& graph,
const typename G::EdgeIterator& edge,
const ColIter& col)
1487 if(-norm_(*col) >= maxValue_ * alpha()) {
1488 edge.properties().setDepends();
1489 typedef typename G::EdgeDescriptor ED;
1490 ED e= graph.findEdge(edge.target(), edge.source());
1491 if(e!=std::numeric_limits<ED>::max())
1493 typename G::EdgeProperties& other = graph.getEdgeProperties(e);
1494 other.setInfluences();
1499 template<
class M,
class N>
1500 inline bool Dependency<M,N>::isIsolated()
1502 return maxValue_ < beta() * diagonal_;
1505 template<
class G,
class S>
1507 VertexSet& connected, std::vector<Vertex>& front)
1508 : vertices_(), id_(-1), graph_(graph), aggregates_(aggregates),
1509 connected_(connected), front_(
front)
1512 template<
class G,
class S>
1520 throw "Not yet implemented";
1528 template<
class G,
class S>
1531 dvverb<<
"Connected cleared"<<std::endl;
1534 connected_.insert(vertex);
1535 dvverb <<
" Inserting "<<vertex<<
" size="<<connected_.size();
1541 template<
class G,
class S>
1544 vertices_.insert(vertex);
1545 aggregates_[vertex]=id_;
1547 front_.erase(std::lower_bound(front_.begin(), front_.end(), vertex));
1551 const iterator end = graph_.endEdges(vertex);
1552 for(iterator edge = graph_.beginEdges(vertex); edge != end; ++edge) {
1553 dvverb <<
" Inserting "<<aggregates_[edge.target()];
1554 connected_.insert(aggregates_[edge.target()]);
1555 dvverb <<
" size="<<connected_.size();
1557 !graph_.getVertexProperties(edge.target()).front())
1559 front_.push_back(edge.target());
1560 graph_.getVertexProperties(edge.target()).setFront();
1564 std::sort(front_.begin(), front_.end());
1567 template<
class G,
class S>
1571 std::size_t oldsize = vertices_.size();
1573 typedef typename std::vector<Vertex>::iterator Iterator;
1575 typedef typename VertexSet::iterator SIterator;
1577 SIterator pos=vertices_.begin();
1578 std::vector<Vertex> newFront;
1579 newFront.reserve(front_.capacity());
1581 std::set_difference(front_.begin(), front_.end(), vertices.begin(), vertices.end(),
1582 std::back_inserter(newFront));
1585 for(Iterator vertex=vertices.begin(); vertex != vertices.end(); ++vertex)
1587 pos=vertices_.insert(pos,*vertex);
1588 vertices_.insert(*vertex);
1589 graph_.getVertexProperties(*vertex).resetFront();
1590 aggregates_[*vertex]=id_;
1593 const iterator end = graph_.endEdges(*vertex);
1594 for(iterator edge = graph_.beginEdges(*vertex); edge != end; ++edge) {
1595 dvverb <<
" Inserting "<<aggregates_[edge.target()];
1596 connected_.insert(aggregates_[edge.target()]);
1598 !graph_.getVertexProperties(edge.target()).front())
1600 front_.push_back(edge.target());
1601 graph_.getVertexProperties(edge.target()).setFront();
1603 dvverb <<
" size="<<connected_.size();
1607 std::sort(front_.begin(), front_.end());
1608 assert(oldsize+vertices.size()==vertices_.size());
1610 template<
class G,
class S>
1618 template<
class G,
class S>
1619 inline typename Aggregate<G,S>::VertexSet::size_type
1622 return vertices_.size();
1625 template<
class G,
class S>
1626 inline typename Aggregate<G,S>::VertexSet::size_type
1629 return connected_.size();
1632 template<
class G,
class S>
1638 template<
class G,
class S>
1641 return vertices_.begin();
1644 template<
class G,
class S>
1647 return vertices_.end();
1665 delete[] aggregates_;
1672 allocate(noVertices);
1685 noVertices_ = noVertices;
1687 for(std::size_t i=0; i < noVertices; i++)
1688 aggregates_[i]=UNAGGREGATED;
1694 assert(aggregates_ != 0);
1695 delete[] aggregates_;
1703 return aggregates_[v];
1710 return aggregates_[v];
1714 template<
bool reset,
class G,
class F,
class VM>
1717 const G& graph, F& aggregateVisitor,
1718 VM& visitedMap)
const
1722 DummyEdgeVisitor dummy;
1723 return breadthFirstSearch<true,reset>(start, aggregate, graph, vlist, aggregateVisitor, dummy, visitedMap);
1727 template<
bool remove,
bool reset,
class G,
class L,
class F1,
class F2,
class VM>
1732 F1& aggregateVisitor,
1733 F2& nonAggregateVisitor,
1734 VM& visitedMap)
const
1736 typedef typename L::const_iterator ListIterator;
1737 int visitedSpheres = 0;
1739 visited.push_back(start);
1740 put(visitedMap, start,
true);
1742 ListIterator current = visited.begin();
1743 ListIterator end = visited.end();
1744 std::size_t i=0, size=visited.size();
1748 while(current != end) {
1750 for(; i<size; ++current, ++i) {
1751 typedef typename G::ConstEdgeIterator EdgeIterator;
1752 const EdgeIterator endEdge = graph.endEdges(*current);
1754 for(EdgeIterator edge = graph.beginEdges(*current);
1755 edge != endEdge; ++edge) {
1757 if(aggregates_[edge.target()]==aggregate) {
1758 if(!get(visitedMap, edge.target())) {
1759 put(visitedMap, edge.target(),
true);
1760 visited.push_back(edge.target());
1761 aggregateVisitor(edge);
1764 nonAggregateVisitor(edge);
1767 end = visited.end();
1768 size = visited.size();
1774 for(current = visited.begin(); current != end; ++current)
1775 put(visitedMap, *current,
false);
1781 return visitedSpheres;
1786 : graph_(0), aggregate_(0), front_(), connected_(), size_(-1)
1795 template<
class G,
class C>
1797 const typename C::Matrix& matrix,
1798 C criterion,
bool firstlevel)
1801 typedef typename C::Matrix Matrix;
1802 typedef typename G::VertexIterator VertexIterator;
1804 criterion.init(&matrix);
1806 for(VertexIterator vertex = graph.begin(); vertex != graph.end(); ++vertex) {
1809 const Row& row = matrix[*vertex];
1814 criterion.initRow(row, *vertex);
1819 ColIterator end = row.end();
1820 typename FieldTraits<typename Matrix::field_type>::real_type absoffdiag=0.;
1823 for(ColIterator col = row.begin(); col != end; ++col)
1824 if(col.index()!=*vertex) {
1825 criterion.examine(col);
1826 absoffdiag=std::max(absoffdiag, col->frobenius_norm());
1830 vertex.properties().setExcludedBorder();
1833 for(ColIterator col = row.begin(); col != end; ++col)
1834 if(col.index()!=*vertex)
1835 criterion.examine(col);
1841 if(criterion.isIsolated()) {
1843 vertex.properties().setIsolated();
1846 typedef typename G::EdgeIterator EdgeIterator;
1848 EdgeIterator eEnd = vertex.end();
1849 ColIterator col = matrix[*vertex].begin();
1851 for(EdgeIterator edge = vertex.begin(); edge!= eEnd; ++edge, ++col) {
1853 while(col.index()!=edge.target())
1855 criterion.examine(graph, edge, col);
1865 inline Aggregator<G>::AggregateVisitor<V>::AggregateVisitor(
const AggregatesMap<Vertex>& aggregates,
1867 : aggregates_(aggregates), aggregate_(aggregate), visitor_(&visitor)
1874 if(aggregates_[edge.target()]==aggregate_)
1875 visitor_->operator()(edge);
1880 inline void Aggregator<G>::visitAggregateNeighbours(
const Vertex& vertex,
1882 const AggregatesMap<Vertex>& aggregates,
1886 AggregateVisitor<V> v(aggregates, aggregate, visitor);
1892 inline Aggregator<G>::Counter::Counter()
1897 inline void Aggregator<G>::Counter::increment()
1903 inline void Aggregator<G>::Counter::decrement()
1908 inline int Aggregator<G>::Counter::value()
1916 if(edge.properties().isTwoWay())
1917 Counter::increment();
1922 const AggregatesMap<Vertex>& aggregates)
const
1924 TwoWayCounter counter;
1925 visitAggregateNeighbours(vertex, aggregate, aggregates, counter);
1926 return counter.value();
1931 const AggregatesMap<Vertex>& aggregates)
const
1933 OneWayCounter counter;
1934 visitAggregateNeighbours(vertex, aggregate, aggregates, counter);
1935 return counter.value();
1941 if(edge.properties().isOneWay())
1942 Counter::increment();
1946 inline Aggregator<G>::ConnectivityCounter::ConnectivityCounter(
const VertexSet& connected,
1947 const AggregatesMap<Vertex>& aggregates)
1948 : Counter(), connected_(connected), aggregates_(aggregates)
1957 Counter::increment();
1959 Counter::increment();
1960 Counter::increment();
1965 inline double Aggregator<G>::connectivity(
const Vertex& vertex,
const AggregatesMap<Vertex>& aggregates)
const
1967 ConnectivityCounter counter(connected_, aggregates);
1969 return (
double)counter.value()/noNeighbours;
1973 inline Aggregator<G>::DependencyCounter::DependencyCounter()
1980 if(edge.properties().depends())
1981 Counter::increment();
1982 if(edge.properties().influences())
1983 Counter::increment();
1987 int Aggregator<G>::unusedNeighbours(
const Vertex& vertex,
const AggregatesMap<Vertex>& aggregates)
const
1993 std::pair<int,int> Aggregator<G>::neighbours(
const Vertex& vertex,
1995 const AggregatesMap<Vertex>& aggregates)
const
1997 DependencyCounter unused, aggregated;
1998 typedef AggregateVisitor<DependencyCounter> Counter;
1999 typedef std::tuple<Counter,Counter> CounterTuple;
2002 return std::make_pair(unused.value(), aggregated.value());
2007 int Aggregator<G>::aggregateNeighbours(
const Vertex& vertex,
const AggregateDescriptor& aggregate,
const AggregatesMap<Vertex>& aggregates)
const
2009 DependencyCounter counter;
2010 visitAggregateNeighbours(vertex, aggregate, aggregates, counter);
2011 return counter.value();
2015 std::size_t Aggregator<G>::distance(
const Vertex& vertex,
const AggregatesMap<Vertex>& aggregates)
2018 typename PropertyMapTypeSelector<VertexVisitedTag,G>::Type visitedMap = get(VertexVisitedTag(), *graph_);
2020 typename AggregatesMap<Vertex>::DummyEdgeVisitor dummy;
2021 return aggregates.template breadthFirstSearch<true,true>(vertex,
2022 aggregate_->
id(), *graph_,
2023 vlist, dummy, dummy, visitedMap);
2027 inline Aggregator<G>::FrontMarker::FrontMarker(std::vector<Vertex>& front,
MatrixGraph& graph)
2028 : front_(
front), graph_(graph)
2034 Vertex target = edge.target();
2036 if(!graph_.getVertexProperties(target).front()) {
2037 front_.push_back(target);
2038 graph_.getVertexProperties(target).setFront();
2043 inline bool Aggregator<G>::admissible(
const Vertex& vertex,
const AggregateDescriptor& aggregate,
const AggregatesMap<Vertex>& aggregates)
const
2046 Dune::dvverb<<
" Admissible not yet implemented!"<<std::endl;
2053 Iterator vend = graph_->endEdges(vertex);
2054 for(Iterator edge = graph_->beginEdges(vertex); edge != vend; ++edge) {
2056 if(edge.properties().isStrong()
2057 && aggregates[edge.target()]==aggregate)
2060 Iterator edge1 = edge;
2061 for(++edge1; edge1 != vend; ++edge1) {
2063 if(edge1.properties().isStrong()
2064 && aggregates[edge.target()]==aggregate)
2069 Iterator v2end = graph_->endEdges(edge.target());
2070 for(Iterator edge2 = graph_->beginEdges(edge.target()); edge2 != v2end; ++edge2) {
2071 if(edge2.target()==edge1.target() &&
2072 edge2.properties().isStrong()) {
2088 vend = graph_->endEdges(vertex);
2089 for(Iterator edge = graph_->beginEdges(vertex); edge != vend; ++edge) {
2091 if(edge.properties().isStrong()
2092 && aggregates[edge.target()]==aggregate)
2095 Iterator v1end = graph_->endEdges(edge.target());
2097 for(Iterator edge1=graph_->beginEdges(edge.target()); edge1 != v1end; ++edge1) {
2099 if(edge1.properties().isStrong()
2100 && aggregates[edge1.target()]==aggregate)
2104 Iterator v2end = graph_->endEdges(vertex);
2105 for(Iterator edge2 = graph_->beginEdges(vertex); edge2 != v2end; ++edge2) {
2106 if(edge2.target()==edge1.target()) {
2107 if(edge2.properties().isStrong())
2124 void Aggregator<G>::unmarkFront()
2126 typedef typename std::vector<Vertex>::const_iterator Iterator;
2128 for(Iterator vertex=front_.begin(); vertex != front_.end(); ++vertex)
2129 graph_->getVertexProperties(*vertex).resetFront();
2136 Aggregator<G>::nonisoNeighbourAggregate(
const Vertex& vertex,
2137 const AggregatesMap<Vertex>& aggregates,
2138 SLList<Vertex>& neighbours)
const
2141 Iterator end=graph_->beginEdges(vertex);
2144 for(Iterator edge=graph_->beginEdges(vertex); edge!=end; ++edge)
2147 neighbours.push_back(aggregates[edge.target()]);
2152 inline typename G::VertexDescriptor Aggregator<G>::mergeNeighbour(
const Vertex& vertex,
const AggregatesMap<Vertex>& aggregates)
const
2156 Iterator end = graph_->endEdges(vertex);
2157 for(Iterator edge = graph_->beginEdges(vertex); edge != end; ++edge) {
2159 graph_->getVertexProperties(edge.target()).isolated() == graph_->getVertexProperties(edge.source()).isolated()) {
2160 if( graph_->getVertexProperties(vertex).isolated() ||
2161 ((edge.properties().depends() || edge.properties().influences())
2162 && admissible(vertex, aggregates[edge.target()], aggregates)))
2163 return edge.target();
2170 Aggregator<G>::FrontNeighbourCounter::FrontNeighbourCounter(
const MatrixGraph& graph)
2171 : Counter(), graph_(graph)
2177 if(graph_.getVertexProperties(edge.target()).front())
2178 Counter::increment();
2182 int Aggregator<G>::noFrontNeighbours(
const Vertex& vertex)
const
2184 FrontNeighbourCounter counter(*graph_);
2186 return counter.value();
2189 inline bool Aggregator<G>::connected(
const Vertex& vertex,
2191 const AggregatesMap<Vertex>& aggregates)
const
2193 typedef typename G::ConstEdgeIterator iterator;
2194 const iterator end = graph_->endEdges(vertex);
2195 for(iterator edge = graph_->beginEdges(vertex); edge != end; ++edge)
2196 if(aggregates[edge.target()]==aggregate)
2201 inline bool Aggregator<G>::connected(
const Vertex& vertex,
2202 const SLList<AggregateDescriptor>& aggregateList,
2203 const AggregatesMap<Vertex>& aggregates)
const
2206 for(Iter i=aggregateList.begin(); i!=aggregateList.end(); ++i)
2207 if(connected(vertex, *i, aggregates))
2214 void Aggregator<G>::growIsolatedAggregate(
const Vertex& seed,
const AggregatesMap<Vertex>& aggregates,
const C& c)
2216 SLList<Vertex> connectedAggregates;
2217 nonisoNeighbourAggregate(seed, aggregates,connectedAggregates);
2219 while(aggregate_->
size()< c.minAggregateSize() && aggregate_->
connectSize() < c.maxConnectivity()) {
2221 std::size_t maxFrontNeighbours=0;
2225 typedef typename std::vector<Vertex>::const_iterator Iterator;
2227 for(Iterator vertex = front_.begin(); vertex != front_.end(); ++vertex) {
2228 if(distance(*vertex, aggregates)>c.maxDistance())
2231 if(connectedAggregates.size()>0) {
2235 if(!connected(*vertex, connectedAggregates, aggregates))
2239 double con = connectivity(*vertex, aggregates);
2242 std::size_t frontNeighbours = noFrontNeighbours(*vertex);
2244 if(frontNeighbours >= maxFrontNeighbours) {
2245 maxFrontNeighbours = frontNeighbours;
2246 candidate = *vertex;
2248 }
else if(con > maxCon) {
2250 maxFrontNeighbours = noFrontNeighbours(*vertex);
2251 candidate = *vertex;
2258 aggregate_->
add(candidate);
2264 void Aggregator<G>::growAggregate(
const Vertex& seed,
const AggregatesMap<Vertex>& aggregates,
const C& c)
2267 std::size_t distance_ =0;
2268 while(aggregate_->
size() < c.minAggregateSize()&& distance_<c.maxDistance()) {
2269 int maxTwoCons=0, maxOneCons=0, maxNeighbours=-1;
2272 std::vector<Vertex> candidates;
2273 candidates.reserve(30);
2275 typedef typename std::vector<Vertex>::const_iterator Iterator;
2277 for(Iterator vertex = front_.begin(); vertex != front_.end(); ++vertex) {
2279 if(graph_->getVertexProperties(*vertex).isolated())
2282 int twoWayCons = twoWayConnections(*vertex, aggregate_->
id(), aggregates);
2285 if( maxTwoCons == twoWayCons && twoWayCons > 0) {
2286 double con = connectivity(*vertex, aggregates);
2289 int neighbours = noFrontNeighbours(*vertex);
2291 if(neighbours > maxNeighbours) {
2292 maxNeighbours = neighbours;
2294 candidates.push_back(*vertex);
2296 candidates.push_back(*vertex);
2298 }
else if( con > maxCon) {
2300 maxNeighbours = noFrontNeighbours(*vertex);
2302 candidates.push_back(*vertex);
2304 }
else if(twoWayCons > maxTwoCons) {
2305 maxTwoCons = twoWayCons;
2306 maxCon = connectivity(*vertex, aggregates);
2307 maxNeighbours = noFrontNeighbours(*vertex);
2309 candidates.push_back(*vertex);
2312 maxOneCons = std::numeric_limits<int>::max();
2321 int oneWayCons = oneWayConnections(*vertex, aggregate_->
id(), aggregates);
2326 if(!admissible(*vertex, aggregate_->
id(), aggregates))
2329 if( maxOneCons == oneWayCons && oneWayCons > 0) {
2330 double con = connectivity(*vertex, aggregates);
2333 int neighbours = noFrontNeighbours(*vertex);
2335 if(neighbours > maxNeighbours) {
2336 maxNeighbours = neighbours;
2338 candidates.push_back(*vertex);
2340 if(neighbours==maxNeighbours)
2342 candidates.push_back(*vertex);
2345 }
else if( con > maxCon) {
2347 maxNeighbours = noFrontNeighbours(*vertex);
2349 candidates.push_back(*vertex);
2351 }
else if(oneWayCons > maxOneCons) {
2352 maxOneCons = oneWayCons;
2353 maxCon = connectivity(*vertex, aggregates);
2354 maxNeighbours = noFrontNeighbours(*vertex);
2356 candidates.push_back(*vertex);
2361 if(!candidates.size())
2363 distance_=distance(seed, aggregates);
2364 candidates.resize(std::min(candidates.size(), c.maxAggregateSize()-
2365 aggregate_->
size()));
2366 aggregate_->
add(candidates);
2370 template<
typename V>
2371 template<
typename M,
typename G,
typename C>
2375 Aggregator<G> aggregator;
2376 return aggregator.build(matrix, graph, *
this, criterion, finestLevel);
2380 template<
class M,
class C>
2381 std::tuple<int,int,int,int>
Aggregator<G>::build(
const M& m, G& graph, AggregatesMap<Vertex>& aggregates,
const C& c,
2385 Stack stack_(graph, *
this, aggregates);
2389 aggregate_ =
new Aggregate<G,VertexSet>(graph, aggregates, connected_, front_);
2396 dverb<<
"Build dependency took "<< watch.elapsed()<<
" seconds."<<std::endl;
2397 int noAggregates, conAggregates, isoAggregates, oneAggregates;
2398 std::size_t maxA=0, minA=1000000, avg=0;
2399 int skippedAggregates;
2400 noAggregates = conAggregates = isoAggregates = oneAggregates =
2401 skippedAggregates = 0;
2404 Vertex seed = stack_.pop();
2406 if(seed == Stack::NullEntry)
2411 if((noAggregates+1)%10000 == 0)
2415 if(graph.getVertexProperties(seed).excludedBorder()) {
2417 ++skippedAggregates;
2421 if(graph.getVertexProperties(seed).isolated()) {
2422 if(c.skipIsolated()) {
2425 ++skippedAggregates;
2429 aggregate_->
seed(seed);
2430 growIsolatedAggregate(seed, aggregates, c);
2433 aggregate_->
seed(seed);
2434 growAggregate(seed, aggregates, c);
2438 while(!(graph.getVertexProperties(seed).isolated()) && aggregate_->
size() < c.maxAggregateSize()) {
2440 std::vector<Vertex> candidates;
2441 candidates.reserve(30);
2443 typedef typename std::vector<Vertex>::const_iterator Iterator;
2445 for(Iterator vertex = front_.begin(); vertex != front_.end(); ++vertex) {
2447 if(graph.getVertexProperties(*vertex).isolated())
2450 if(twoWayConnections( *vertex, aggregate_->
id(), aggregates) == 0 &&
2451 (oneWayConnections( *vertex, aggregate_->
id(), aggregates) == 0 ||
2452 !admissible( *vertex, aggregate_->
id(), aggregates) ))
2455 std::pair<int,int> neighbourPair=neighbours(*vertex, aggregate_->
id(),
2461 if(neighbourPair.first >= neighbourPair.second)
2464 if(distance(*vertex, aggregates) > c.maxDistance())
2466 candidates.push_back(*vertex);
2470 if(!candidates.size())
break;
2472 candidates.resize(std::min(candidates.size(), c.maxAggregateSize()-
2473 aggregate_->
size()));
2474 aggregate_->
add(candidates);
2479 if(aggregate_->
size()==1 && c.maxAggregateSize()>1) {
2480 if(!graph.getVertexProperties(seed).isolated()) {
2481 Vertex mergedNeighbour = mergeNeighbour(seed, aggregates);
2485 aggregates[seed] = aggregates[mergedNeighbour];
2486 aggregate_->invalidate();
2489 minA=std::min(minA,
static_cast<std::size_t
>(1));
2490 maxA=std::max(maxA,
static_cast<std::size_t
>(1));
2496 minA=std::min(minA,
static_cast<std::size_t
>(1));
2497 maxA=std::max(maxA,
static_cast<std::size_t
>(1));
2503 avg+=aggregate_->
size();
2504 minA=std::min(minA,aggregate_->
size());
2505 maxA=std::max(maxA,aggregate_->
size());
2506 if(graph.getVertexProperties(seed).isolated())
2514 Dune::dinfo<<
"connected aggregates: "<<conAggregates;
2515 Dune::dinfo<<
" isolated aggregates: "<<isoAggregates;
2516 if(conAggregates+isoAggregates>0)
2517 Dune::dinfo<<
" one node aggregates: "<<oneAggregates<<
" min size="
2518 <<minA<<
" max size="<<maxA
2519 <<
" avg="<<avg/(conAggregates+isoAggregates)<<std::endl;
2522 return std::make_tuple(conAggregates+isoAggregates,isoAggregates,
2523 oneAggregates,skippedAggregates);
2528 Aggregator<G>::Stack::Stack(
const MatrixGraph& graph,
const Aggregator<G>& aggregatesBuilder,
2529 const AggregatesMap<Vertex>& aggregates)
2530 : graph_(graph), aggregatesBuilder_(aggregatesBuilder), aggregates_(aggregates), begin_(graph.begin()), end_(graph.end())
2536 Aggregator<G>::Stack::~Stack()
2544 = std::numeric_limits<typename G::VertexDescriptor>::max();
2547 inline typename G::VertexDescriptor Aggregator<G>::Stack::pop()
2553 typename G::VertexDescriptor current=*begin_;
2563 void printAggregates2d(
const AggregatesMap<V>& aggregates,
int n,
int m, std::ostream& os)
2565 std::ios_base::fmtflags oldOpts=os.flags();
2567 os.setf(std::ios_base::right, std::ios_base::adjustfield);
2572 for(
int i=0; i< n*m; i++)
2573 max=std::max(max, aggregates[i]);
2575 for(
int i=10; i < 1000000; i*=10)
2581 for(
int j=0, entry=0; j < m; j++) {
2582 for(
int i=0; i<n; i++, entry++) {
2584 os<<aggregates[entry]<<
" ";
A class for temporarily storing the vertices of an aggregate in.
Definition: aggregates.hh:773
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:904
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
Definition: bvector.hh:650
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_THROW(E, m)
Definition: exceptions.hh:216
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:915
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:910
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:1059
std::size_t * SphereMap
Type of the mapping of aggregate members onto distance spheres.
Definition: aggregates.hh:804
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:799
MatrixGraph::VertexDescriptor AggregateDescriptor
The type of the aggregate descriptor.
Definition: aggregates.hh:918
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:796
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:784
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:790
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
Front front
PartitionSet for the front partition.
Definition: partitionset.hh:229
Dune namespace.
Definition: alignment.hh:11
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.
#define DUNE_UNUSED_PARAMETER(parm)
A macro to mark intentionally unused function parameters with.
Definition: unused.hh:18