5#ifndef DUNE_AMG_AGGREGATES_HH
6#define DUNE_AMG_AGGREGATES_HH
12#include "combinedfunctor.hh"
84 this->setMaxDistance(diameter-1);
89 this->setMaxDistance(this->maxDistance()+diameter-1);
91 this->setMinAggregateSize(csize);
92 this->setMaxAggregateSize(
static_cast<std::size_t
>(csize*1.5));
108 this->setMaxDistance(this->maxDistance()+dim-1);
113 std::ostream& operator<<(std::ostream& os,
const AggregationCriterion<T>& criterion)
115 os<<
"{ maxdistance="<<criterion.maxDistance()<<
" minAggregateSize="
116 <<criterion.minAggregateSize()<<
" maxAggregateSize="<<criterion.maxAggregateSize()
117 <<
" connectivity="<<criterion.maxConnectivity()<<
" debugLevel="<<criterion.debugLevel()<<
"}";
132 template<
class M,
class N>
156 void init(
const Matrix* matrix);
158 void initRow(
const Row& row,
int index);
160 void examine(
const ColIter& col);
163 void examine(G& graph,
const typename G::EdgeIterator& edge,
const ColIter& col);
180 typedef typename FieldTraits<field_type>::real_type real_type;
188 std::vector<real_type> vals_;
189 typename std::vector<real_type>::iterator valIter_;
194 template<
class M,
class N>
200 template<
class M,
class N>
201 inline void SymmetricMatrixDependency<M,N>::initRow(
const Row& row,
int index)
204 vals_.assign(row.size(), 0.0);
205 assert(vals_.size()==row.size());
206 valIter_=vals_.begin();
209 diagonal_=norm_(row[index]);
213 template<
class M,
class N>
214 inline void SymmetricMatrixDependency<M,N>::examine(
const ColIter& col)
218 real_type eij = norm_(*col);
219 if(!N::is_sign_preserving || eij<0)
221 *valIter_ = eij/diagonal_*eij/norm_(matrix_->operator[](col.index())[col.index()]);
222 maxValue_ =
max(maxValue_, *valIter_);
228 template<
class M,
class N>
230 inline void SymmetricMatrixDependency<M,N>::examine(G&,
const typename G::EdgeIterator& edge,
const ColIter&)
232 if(*valIter_ > alpha() * maxValue_) {
233 edge.properties().setDepends();
234 edge.properties().setInfluences();
239 template<
class M,
class N>
240 inline bool SymmetricMatrixDependency<M,N>::isIsolated()
244 valIter_=vals_.begin();
245 return maxValue_ < beta();
251 template<
class M,
class N>
275 void init(
const Matrix* matrix);
277 void initRow(
const Row& row,
int index);
279 void examine(
const ColIter& col);
282 void examine(G& graph,
const typename G::EdgeIterator& edge,
const ColIter& col);
299 typedef typename FieldTraits<field_type>::real_type real_type;
312 template<
class M,
class N>
336 void init(
const Matrix* matrix);
338 void initRow(
const Row& row,
int index);
340 void examine(
const ColIter& col);
343 void examine(G& graph,
const typename G::EdgeIterator& edge,
const ColIter& col);
360 typedef typename FieldTraits<field_type>::real_type real_type;
369 void initRow(
const Row& row,
int index,
const std::true_type&);
370 void initRow(
const Row& row,
int index,
const std::false_type&);
382 is_sign_preserving =
true
390 typename FieldTraits<typename M::field_type>::real_type
operator()(
const M& m,
393 typedef typename M::field_type field_type;
394 typedef typename FieldTraits<field_type>::real_type real_type;
395 static_assert( std::is_convertible<field_type, real_type >::value,
396 "use of diagonal norm in AMG not implemented for complex field_type");
409 typedef typename FieldTraits<M>::real_type real_type;
410 static_assert( std::is_convertible<M, real_type >::value,
411 "use of diagonal norm in AMG not implemented for complex field_type");
420 static T signed_abs(
const T & v)
427 static T signed_abs(
const std::complex<T> & v)
431 return csgn(v) * std::abs(v);
436 static T csgn(
const T & v)
438 return (T(0) < v) - (v < T(0));
443 static T csgn(std::complex<T> a)
445 return csgn(a.real())+(a.real() == 0.0)*csgn(a.imag());
466 is_sign_preserving =
false
473 typename FieldTraits<typename M::field_type>::real_type
operator()(
const M& m)
const
475 return m.infinity_norm();
483 is_sign_preserving =
false
490 typename FieldTraits<typename M::field_type>::real_type
operator()(
const M& m)
const
492 return m.frobenius_norm();
499 is_sign_preserving =
false
506 typename FieldTraits<typename M::field_type>::real_type
operator()(
const M& m)
const
517 template<
class M,
class Norm>
537 template<
class M,
class Norm>
600 template<
class EdgeIterator>
601 void operator()([[maybe_unused]]
const EdgeIterator& edge)
const
634 template<
class M,
class G,
class C>
635 std::tuple<int,int,int,int>
buildAggregates(
const M& matrix, G& graph,
const C& criterion,
655 template<
bool reset,
class G,
class F,
class VM>
660 VM& visitedMap)
const;
685 template<
bool remove,
bool reset,
class G,
class L,
class F1,
class F2,
class VM>
688 const G& graph, L& visited, F1& aggregateVisitor,
689 F2& nonAggregateVisitor,
690 VM& visitedMap)
const;
725 const_iterator begin()
const
730 const_iterator end()
const
750 AggregatesMap<V>& operator=(
const AggregatesMap<V>&) =
delete;
760 std::size_t noVertices_;
766 template<
class G,
class C>
768 const typename C::Matrix& matrix,
776 template<
class G,
class S>
785 typedef G MatrixGraph;
820 VertexSet& connectivity, std::vector<Vertex>& front_);
854 typename VertexSet::size_type
size();
901 std::vector<Vertex>& front_;
951 template<
class M,
class C>
952 std::tuple<int,int,int,int>
build(
const M& m, G& graph,
970 typedef std::set<Vertex,std::less<Vertex>,
Allocator> VertexSet;
975 typedef std::size_t* SphereMap;
990 std::vector<Vertex> front_;
995 VertexSet connected_;
1008 static const Vertex NullEntry;
1016 enum { N = 1300000 };
1050 const AggregatesMap<Vertex>& aggregates,
1058 class AggregateVisitor
1118 class FrontNeighbourCounter :
public Counter
1142 class TwoWayCounter :
public Counter
1165 class OneWayCounter :
public Counter
1183 const AggregatesMap<Vertex>& aggregates)
const;
1191 class ConnectivityCounter :
public Counter
1206 const VertexSet& connected_;
1251 class DependencyCounter :
public Counter
1283 std::vector<Vertex>& front_;
1398 template<
class M,
class N>
1404 template<
class M,
class N>
1405 inline void SymmetricDependency<M,N>::initRow(
const Row& row,
int index)
1407 initRow(row, index, std::is_convertible<field_type, real_type>());
1410 template<
class M,
class N>
1411 inline void SymmetricDependency<M,N>::initRow(
const Row& row,
int index,
const std::false_type&)
1413 DUNE_THROW(InvalidStateException,
"field_type needs to convertible to real_type");
1416 template<
class M,
class N>
1417 inline void SymmetricDependency<M,N>::initRow([[maybe_unused]]
const Row& row,
int index,
const std::true_type&)
1422 diagonal_ = norm_(matrix_->operator[](row_)[row_]);
1425 template<
class M,
class N>
1426 inline void SymmetricDependency<M,N>::examine(
const ColIter& col)
1429 real_type eij = norm_(*col);
1431 matrix_->operator[](col.index()).find(row_);
1432 if ( opposite_entry == matrix_->operator[](col.index()).end() )
1437 real_type eji = norm_(*opposite_entry);
1440 if(!N::is_sign_preserving || eij<0 || eji<0)
1441 maxValue_ =
max(maxValue_,
1442 eij /diagonal_ * eji/
1443 norm_(matrix_->operator[](col.index())[col.index()]));
1446 template<
class M,
class N>
1448 inline void SymmetricDependency<M,N>::examine(G& graph,
const typename G::EdgeIterator& edge,
const ColIter& col)
1450 real_type eij = norm_(*col);
1452 matrix_->operator[](col.index()).find(row_);
1454 if ( opposite_entry == matrix_->operator[](col.index()).end() )
1459 real_type eji = norm_(*opposite_entry);
1461 if(!N::is_sign_preserving || (eij<0 || eji<0))
1462 if(eji / norm_(matrix_->operator[](edge.target())[edge.target()]) *
1463 eij/ diagonal_ > alpha() * maxValue_) {
1464 edge.properties().setDepends();
1465 edge.properties().setInfluences();
1466 typename G::EdgeProperties& other = graph.getEdgeProperties(edge.target(), edge.source());
1467 other.setInfluences();
1472 template<
class M,
class N>
1473 inline bool SymmetricDependency<M,N>::isIsolated()
1475 return maxValue_ < beta();
1479 template<
class M,
class N>
1480 inline void Dependency<M,N>::init(
const Matrix* matrix)
1485 template<
class M,
class N>
1486 inline void Dependency<M,N>::initRow([[maybe_unused]]
const Row& row,
int index)
1491 diagonal_ = norm_(matrix_->operator[](row_)[row_]);
1494 template<
class M,
class N>
1495 inline void Dependency<M,N>::examine(
const ColIter& col)
1498 maxValue_ =
max(maxValue_, -norm_(*col));
1501 template<
class M,
class N>
1503 inline void Dependency<M,N>::examine(G& graph,
const typename G::EdgeIterator& edge,
const ColIter& col)
1505 if(-norm_(*col) >= maxValue_ * alpha()) {
1506 edge.properties().setDepends();
1507 typedef typename G::EdgeDescriptor ED;
1508 ED e= graph.findEdge(edge.target(), edge.source());
1511 typename G::EdgeProperties& other = graph.getEdgeProperties(e);
1512 other.setInfluences();
1517 template<
class M,
class N>
1518 inline bool Dependency<M,N>::isIsolated()
1520 return maxValue_ < beta() * diagonal_;
1523 template<
class G,
class S>
1525 VertexSet& connected, std::vector<Vertex>&
front)
1526 : vertices_(), id_(-1), graph_(graph), aggregates_(aggregates),
1527 connected_(connected), front_(
front)
1530 template<
class G,
class S>
1538 throw "Not yet implemented";
1546 template<
class G,
class S>
1549 dvverb<<
"Connected cleared"<<std::endl;
1552 connected_.insert(
vertex);
1553 dvverb <<
" Inserting "<<
vertex<<
" size="<<connected_.size();
1559 template<
class G,
class S>
1562 vertices_.insert(
vertex);
1565 front_.erase(std::lower_bound(front_.begin(), front_.end(),
vertex));
1569 const iterator end = graph_.endEdges(
vertex);
1570 for(iterator edge = graph_.beginEdges(
vertex); edge != end; ++edge) {
1571 dvverb <<
" Inserting "<<aggregates_[edge.target()];
1572 connected_.insert(aggregates_[edge.target()]);
1573 dvverb <<
" size="<<connected_.size();
1575 !graph_.getVertexProperties(edge.target()).front())
1577 front_.push_back(edge.target());
1578 graph_.getVertexProperties(edge.target()).setFront();
1582 std::sort(front_.begin(), front_.end());
1585 template<
class G,
class S>
1589 std::size_t oldsize = vertices_.size();
1591 typedef typename std::vector<Vertex>::iterator Iterator;
1593 typedef typename VertexSet::iterator SIterator;
1595 SIterator pos=vertices_.begin();
1596 std::vector<Vertex> newFront;
1597 newFront.reserve(front_.capacity());
1599 std::set_difference(front_.begin(), front_.end(), vertices.begin(), vertices.end(),
1600 std::back_inserter(newFront));
1605 pos=vertices_.insert(pos,*
vertex);
1606 vertices_.insert(*
vertex);
1607 graph_.getVertexProperties(*vertex).resetFront();
1608 aggregates_[*
vertex]=id_;
1611 const iterator end = graph_.endEdges(*
vertex);
1612 for(iterator edge = graph_.beginEdges(*
vertex); edge != end; ++edge) {
1613 dvverb <<
" Inserting "<<aggregates_[edge.target()];
1614 connected_.insert(aggregates_[edge.target()]);
1616 !graph_.getVertexProperties(edge.target()).front())
1618 front_.push_back(edge.target());
1619 graph_.getVertexProperties(edge.target()).setFront();
1621 dvverb <<
" size="<<connected_.size();
1625 std::sort(front_.begin(), front_.end());
1626 assert(oldsize+vertices.size()==vertices_.size());
1628 template<
class G,
class S>
1636 template<
class G,
class S>
1637 inline typename Aggregate<G,S>::VertexSet::size_type
1640 return vertices_.size();
1643 template<
class G,
class S>
1644 inline typename Aggregate<G,S>::VertexSet::size_type
1647 return connected_.size();
1650 template<
class G,
class S>
1656 template<
class G,
class S>
1659 return vertices_.begin();
1662 template<
class G,
class S>
1665 return vertices_.end();
1683 delete[] aggregates_;
1690 allocate(noVertices);
1703 noVertices_ = noVertices;
1705 for(std::size_t i=0; i < noVertices; i++)
1706 aggregates_[i]=UNAGGREGATED;
1712 assert(aggregates_ != 0);
1713 delete[] aggregates_;
1721 return aggregates_[v];
1728 return aggregates_[v];
1732 template<
bool reset,
class G,
class F,
class VM>
1735 const G& graph, F& aggregateVisitor,
1736 VM& visitedMap)
const
1740 DummyEdgeVisitor dummy;
1741 return breadthFirstSearch<true,reset>(start, aggregate, graph, vlist, aggregateVisitor, dummy, visitedMap);
1745 template<
bool remove,
bool reset,
class G,
class L,
class F1,
class F2,
class VM>
1750 F1& aggregateVisitor,
1751 F2& nonAggregateVisitor,
1752 VM& visitedMap)
const
1754 typedef typename L::const_iterator ListIterator;
1755 int visitedSpheres = 0;
1757 visited.push_back(start);
1758 put(visitedMap, start,
true);
1760 ListIterator current = visited.begin();
1761 ListIterator end = visited.end();
1762 std::size_t i=0,
size=visited.size();
1766 while(current != end) {
1768 for(; i<
size; ++current, ++i) {
1769 typedef typename G::ConstEdgeIterator EdgeIterator;
1770 const EdgeIterator endEdge = graph.endEdges(*current);
1772 for(EdgeIterator edge = graph.beginEdges(*current);
1773 edge != endEdge; ++edge) {
1775 if(aggregates_[edge.target()]==aggregate) {
1776 if(!
get(visitedMap, edge.target())) {
1777 put(visitedMap, edge.target(),
true);
1778 visited.push_back(edge.target());
1779 aggregateVisitor(edge);
1782 nonAggregateVisitor(edge);
1785 end = visited.end();
1786 size = visited.size();
1792 for(current = visited.begin(); current != end; ++current)
1793 put(visitedMap, *current,
false);
1799 return visitedSpheres;
1804 : graph_(0), aggregate_(0), front_(), connected_(), size_(-1)
1813 template<
class G,
class C>
1815 const typename C::Matrix& matrix,
1816 C criterion,
bool firstlevel)
1819 typedef typename C::Matrix Matrix;
1820 typedef typename G::VertexIterator VertexIterator;
1822 criterion.init(&matrix);
1827 const Row& row = matrix[*
vertex];
1832 criterion.initRow(row, *
vertex);
1837 ColIterator end = row.end();
1838 typename FieldTraits<typename Matrix::field_type>::real_type absoffdiag=0.;
1842 for(ColIterator col = row.begin(); col != end; ++col)
1843 if(col.index()!=*
vertex) {
1844 criterion.examine(col);
1845 absoffdiag =
max(absoffdiag, Impl::asMatrix(*col).frobenius_norm());
1849 vertex.properties().setExcludedBorder();
1852 for(ColIterator col = row.begin(); col != end; ++col)
1854 criterion.examine(col);
1860 if(criterion.isIsolated()) {
1862 vertex.properties().setIsolated();
1865 auto eEnd =
vertex.end();
1866 auto col = matrix[*
vertex].begin();
1868 for(
auto edge =
vertex.begin(); edge!= eEnd; ++edge, ++col) {
1870 while(col.index()!=edge.target())
1872 criterion.examine(graph, edge, col);
1882 inline Aggregator<G>::AggregateVisitor<V>::AggregateVisitor(
const AggregatesMap<Vertex>& aggregates,
1884 : aggregates_(aggregates), aggregate_(aggregate), visitor_(&visitor)
1891 if(aggregates_[edge.target()]==aggregate_)
1892 visitor_->operator()(edge);
1897 inline void Aggregator<G>::visitAggregateNeighbours(
const Vertex&
vertex,
1899 const AggregatesMap<Vertex>& aggregates,
1903 AggregateVisitor<V> v(aggregates, aggregate, visitor);
1909 inline Aggregator<G>::Counter::Counter()
1914 inline void Aggregator<G>::Counter::increment()
1920 inline void Aggregator<G>::Counter::decrement()
1925 inline int Aggregator<G>::Counter::value()
1933 if(edge.properties().isTwoWay())
1934 Counter::increment();
1939 const AggregatesMap<Vertex>& aggregates)
const
1941 TwoWayCounter counter;
1942 visitAggregateNeighbours(
vertex, aggregate, aggregates, counter);
1943 return counter.value();
1948 const AggregatesMap<Vertex>& aggregates)
const
1950 OneWayCounter counter;
1951 visitAggregateNeighbours(
vertex, aggregate, aggregates, counter);
1952 return counter.value();
1958 if(edge.properties().isOneWay())
1959 Counter::increment();
1963 inline Aggregator<G>::ConnectivityCounter::ConnectivityCounter(
const VertexSet& connected,
1964 const AggregatesMap<Vertex>& aggregates)
1965 : Counter(), connected_(connected), aggregates_(aggregates)
1974 Counter::increment();
1976 Counter::increment();
1977 Counter::increment();
1982 inline double Aggregator<G>::connectivity(
const Vertex&
vertex,
const AggregatesMap<Vertex>& aggregates)
const
1984 ConnectivityCounter counter(connected_, aggregates);
1986 return (
double)counter.value()/noNeighbours;
1990 inline Aggregator<G>::DependencyCounter::DependencyCounter()
1997 if(edge.properties().depends())
1998 Counter::increment();
1999 if(edge.properties().influences())
2000 Counter::increment();
2004 int Aggregator<G>::unusedNeighbours(
const Vertex&
vertex,
const AggregatesMap<Vertex>& aggregates)
const
2010 std::pair<int,int> Aggregator<G>::neighbours(
const Vertex&
vertex,
2012 const AggregatesMap<Vertex>& aggregates)
const
2014 DependencyCounter unused, aggregated;
2015 typedef AggregateVisitor<DependencyCounter> CounterT;
2016 typedef std::tuple<CounterT,CounterT> CounterTuple;
2019 return std::make_pair(unused.value(), aggregated.value());
2026 DependencyCounter counter;
2027 visitAggregateNeighbours(
vertex, aggregate, aggregates, counter);
2028 return counter.value();
2032 std::size_t Aggregator<G>::distance(
const Vertex&
vertex,
const AggregatesMap<Vertex>& aggregates)
2035 typename PropertyMapTypeSelector<VertexVisitedTag,G>::Type visitedMap =
get(VertexVisitedTag(), *graph_);
2037 typename AggregatesMap<Vertex>::DummyEdgeVisitor dummy;
2038 return aggregates.template breadthFirstSearch<true,true>(
vertex,
2039 aggregate_->
id(), *graph_,
2040 vlist, dummy, dummy, visitedMap);
2044 inline Aggregator<G>::FrontMarker::FrontMarker(std::vector<Vertex>&
front,
MatrixGraph& graph)
2045 : front_(
front), graph_(graph)
2051 Vertex target = edge.target();
2053 if(!graph_.getVertexProperties(target).front()) {
2054 front_.push_back(target);
2055 graph_.getVertexProperties(target).setFront();
2063 Dune::dvverb<<
" Admissible not yet implemented!"<<std::endl;
2070 Iterator vend = graph_->endEdges(
vertex);
2071 for(Iterator edge = graph_->beginEdges(
vertex); edge != vend; ++edge) {
2073 if(edge.properties().isStrong()
2074 && aggregates[edge.target()]==aggregate)
2077 Iterator edge1 = edge;
2078 for(++edge1; edge1 != vend; ++edge1) {
2080 if(edge1.properties().isStrong()
2081 && aggregates[edge.target()]==aggregate)
2086 Iterator v2end = graph_->endEdges(edge.target());
2087 for(Iterator edge2 = graph_->beginEdges(edge.target()); edge2 != v2end; ++edge2) {
2088 if(edge2.target()==edge1.target() &&
2089 edge2.properties().isStrong()) {
2105 vend = graph_->endEdges(
vertex);
2106 for(Iterator edge = graph_->beginEdges(
vertex); edge != vend; ++edge) {
2108 if(edge.properties().isStrong()
2109 && aggregates[edge.target()]==aggregate)
2112 Iterator v1end = graph_->endEdges(edge.target());
2114 for(Iterator edge1=graph_->beginEdges(edge.target()); edge1 != v1end; ++edge1) {
2116 if(edge1.properties().isStrong()
2117 && aggregates[edge1.target()]==aggregate)
2121 Iterator v2end = graph_->endEdges(
vertex);
2122 for(Iterator edge2 = graph_->beginEdges(
vertex); edge2 != v2end; ++edge2) {
2123 if(edge2.target()==edge1.target()) {
2124 if(edge2.properties().isStrong())
2141 void Aggregator<G>::unmarkFront()
2143 typedef typename std::vector<Vertex>::const_iterator Iterator;
2146 graph_->getVertexProperties(*vertex).resetFront();
2153 Aggregator<G>::nonisoNeighbourAggregate(
const Vertex&
vertex,
2154 const AggregatesMap<Vertex>& aggregates,
2155 SLList<Vertex>& neighbours)
const
2158 Iterator end=graph_->beginEdges(
vertex);
2161 for(Iterator edge=graph_->beginEdges(
vertex); edge!=end; ++edge)
2164 neighbours.push_back(aggregates[edge.target()]);
2169 inline typename G::VertexDescriptor Aggregator<G>::mergeNeighbour(
const Vertex&
vertex,
const AggregatesMap<Vertex>& aggregates)
const
2173 Iterator end = graph_->endEdges(
vertex);
2174 for(Iterator edge = graph_->beginEdges(
vertex); edge != end; ++edge) {
2176 graph_->getVertexProperties(edge.target()).isolated() == graph_->getVertexProperties(edge.source()).isolated()) {
2177 if( graph_->getVertexProperties(
vertex).isolated() ||
2178 ((edge.properties().depends() || edge.properties().influences())
2179 && admissible(
vertex, aggregates[edge.target()], aggregates)))
2180 return edge.target();
2187 Aggregator<G>::FrontNeighbourCounter::FrontNeighbourCounter(
const MatrixGraph& graph)
2188 : Counter(), graph_(graph)
2194 if(graph_.getVertexProperties(edge.target()).front())
2195 Counter::increment();
2199 int Aggregator<G>::noFrontNeighbours(
const Vertex&
vertex)
const
2201 FrontNeighbourCounter counter(*graph_);
2203 return counter.value();
2206 inline bool Aggregator<G>::connected(
const Vertex&
vertex,
2208 const AggregatesMap<Vertex>& aggregates)
const
2210 typedef typename G::ConstEdgeIterator iterator;
2211 const iterator end = graph_->endEdges(
vertex);
2212 for(iterator edge = graph_->beginEdges(
vertex); edge != end; ++edge)
2213 if(aggregates[edge.target()]==aggregate)
2218 inline bool Aggregator<G>::connected(
const Vertex&
vertex,
2219 const SLList<AggregateDescriptor>& aggregateList,
2220 const AggregatesMap<Vertex>& aggregates)
const
2223 for(Iter i=aggregateList.begin(); i!=aggregateList.end(); ++i)
2224 if(connected(
vertex, *i, aggregates))
2231 void Aggregator<G>::growIsolatedAggregate(
const Vertex& seed,
const AggregatesMap<Vertex>& aggregates,
const C& c)
2233 SLList<Vertex> connectedAggregates;
2234 nonisoNeighbourAggregate(seed, aggregates,connectedAggregates);
2236 while(aggregate_->
size()< c.minAggregateSize() && aggregate_->
connectSize() < c.maxConnectivity()) {
2238 std::size_t maxFrontNeighbours=0;
2242 typedef typename std::vector<Vertex>::const_iterator Iterator;
2245 if(distance(*
vertex, aggregates)>c.maxDistance())
2248 if(connectedAggregates.size()>0) {
2252 if(!connected(*
vertex, connectedAggregates, aggregates))
2256 double con = connectivity(*
vertex, aggregates);
2259 std::size_t frontNeighbours = noFrontNeighbours(*
vertex);
2261 if(frontNeighbours >= maxFrontNeighbours) {
2262 maxFrontNeighbours = frontNeighbours;
2265 }
else if(con > maxCon) {
2267 maxFrontNeighbours = noFrontNeighbours(*
vertex);
2275 aggregate_->
add(candidate);
2281 void Aggregator<G>::growAggregate(
const Vertex& seed,
const AggregatesMap<Vertex>& aggregates,
const C& c)
2285 std::size_t distance_ =0;
2286 while(aggregate_->
size() < c.minAggregateSize()&& distance_<c.maxDistance()) {
2287 int maxTwoCons=0, maxOneCons=0, maxNeighbours=-1;
2290 std::vector<Vertex> candidates;
2291 candidates.reserve(30);
2293 typedef typename std::vector<Vertex>::const_iterator Iterator;
2297 if(graph_->getVertexProperties(*vertex).isolated())
2300 int twoWayCons = twoWayConnections(*
vertex, aggregate_->
id(), aggregates);
2303 if( maxTwoCons == twoWayCons && twoWayCons > 0) {
2304 double con = connectivity(*
vertex, aggregates);
2307 int neighbours = noFrontNeighbours(*
vertex);
2309 if(neighbours > maxNeighbours) {
2310 maxNeighbours = neighbours;
2312 candidates.push_back(*
vertex);
2314 candidates.push_back(*
vertex);
2316 }
else if( con > maxCon) {
2318 maxNeighbours = noFrontNeighbours(*
vertex);
2320 candidates.push_back(*
vertex);
2322 }
else if(twoWayCons > maxTwoCons) {
2323 maxTwoCons = twoWayCons;
2324 maxCon = connectivity(*
vertex, aggregates);
2325 maxNeighbours = noFrontNeighbours(*
vertex);
2327 candidates.push_back(*
vertex);
2339 int oneWayCons = oneWayConnections(*
vertex, aggregate_->
id(), aggregates);
2344 if(!admissible(*
vertex, aggregate_->
id(), aggregates))
2347 if( maxOneCons == oneWayCons && oneWayCons > 0) {
2348 double con = connectivity(*
vertex, aggregates);
2351 int neighbours = noFrontNeighbours(*
vertex);
2353 if(neighbours > maxNeighbours) {
2354 maxNeighbours = neighbours;
2356 candidates.push_back(*
vertex);
2358 if(neighbours==maxNeighbours)
2360 candidates.push_back(*
vertex);
2363 }
else if( con > maxCon) {
2365 maxNeighbours = noFrontNeighbours(*
vertex);
2367 candidates.push_back(*
vertex);
2369 }
else if(oneWayCons > maxOneCons) {
2370 maxOneCons = oneWayCons;
2371 maxCon = connectivity(*
vertex, aggregates);
2372 maxNeighbours = noFrontNeighbours(*
vertex);
2374 candidates.push_back(*
vertex);
2379 if(!candidates.size())
2381 distance_=distance(seed, aggregates);
2382 candidates.resize(
min(candidates.size(), c.maxAggregateSize()-
2383 aggregate_->
size()));
2384 aggregate_->
add(candidates);
2388 template<
typename V>
2389 template<
typename M,
typename G,
typename C>
2393 Aggregator<G> aggregator;
2394 return aggregator.build(matrix, graph, *
this, criterion, finestLevel);
2398 template<
class M,
class C>
2399 std::tuple<int,int,int,int>
Aggregator<G>::build(
const M& m, G& graph, AggregatesMap<Vertex>& aggregates,
const C& c,
2405 Stack stack_(graph, *
this, aggregates);
2409 aggregate_ =
new Aggregate<G,VertexSet>(graph, aggregates, connected_, front_);
2416 dverb<<
"Build dependency took "<< watch.elapsed()<<
" seconds."<<std::endl;
2417 int noAggregates, conAggregates, isoAggregates, oneAggregates;
2418 std::size_t maxA=0, minA=1000000, avg=0;
2419 int skippedAggregates;
2420 noAggregates = conAggregates = isoAggregates = oneAggregates =
2421 skippedAggregates = 0;
2424 Vertex seed = stack_.pop();
2426 if(seed == Stack::NullEntry)
2431 if((noAggregates+1)%10000 == 0)
2435 if(graph.getVertexProperties(seed).excludedBorder()) {
2437 ++skippedAggregates;
2441 if(graph.getVertexProperties(seed).isolated()) {
2442 if(c.skipIsolated()) {
2445 ++skippedAggregates;
2449 aggregate_->
seed(seed);
2450 growIsolatedAggregate(seed, aggregates, c);
2453 aggregate_->
seed(seed);
2454 growAggregate(seed, aggregates, c);
2458 while(!(graph.getVertexProperties(seed).isolated()) && aggregate_->
size() < c.maxAggregateSize()) {
2460 std::vector<Vertex> candidates;
2461 candidates.reserve(30);
2463 typedef typename std::vector<Vertex>::const_iterator Iterator;
2467 if(graph.getVertexProperties(*vertex).isolated())
2470 if(twoWayConnections( *
vertex, aggregate_->
id(), aggregates) == 0 &&
2471 (oneWayConnections( *
vertex, aggregate_->
id(), aggregates) == 0 ||
2472 !admissible( *
vertex, aggregate_->
id(), aggregates) ))
2475 std::pair<int,int> neighbourPair=neighbours(*
vertex, aggregate_->
id(),
2481 if(neighbourPair.first >= neighbourPair.second)
2484 if(distance(*
vertex, aggregates) > c.maxDistance())
2486 candidates.push_back(*
vertex);
2490 if(!candidates.size())
break;
2492 candidates.resize(
min(candidates.size(), c.maxAggregateSize()-
2493 aggregate_->
size()));
2494 aggregate_->
add(candidates);
2499 if(aggregate_->
size()==1 && c.maxAggregateSize()>1) {
2500 if(!graph.getVertexProperties(seed).isolated()) {
2501 Vertex mergedNeighbour = mergeNeighbour(seed, aggregates);
2505 aggregates[seed] = aggregates[mergedNeighbour];
2506 aggregate_->invalidate();
2509 minA=
min(minA,
static_cast<std::size_t
>(1));
2510 maxA=
max(maxA,
static_cast<std::size_t
>(1));
2516 minA=
min(minA,
static_cast<std::size_t
>(1));
2517 maxA=
max(maxA,
static_cast<std::size_t
>(1));
2523 avg+=aggregate_->
size();
2524 minA=
min(minA,aggregate_->
size());
2525 maxA=
max(maxA,aggregate_->
size());
2526 if(graph.getVertexProperties(seed).isolated())
2534 Dune::dinfo<<
"connected aggregates: "<<conAggregates;
2535 Dune::dinfo<<
" isolated aggregates: "<<isoAggregates;
2536 if(conAggregates+isoAggregates>0)
2537 Dune::dinfo<<
" one node aggregates: "<<oneAggregates<<
" min size="
2538 <<minA<<
" max size="<<maxA
2539 <<
" avg="<<avg/(conAggregates+isoAggregates)<<std::endl;
2542 return std::make_tuple(conAggregates+isoAggregates,isoAggregates,
2543 oneAggregates,skippedAggregates);
2548 Aggregator<G>::Stack::Stack(
const MatrixGraph& graph,
const Aggregator<G>& aggregatesBuilder,
2549 const AggregatesMap<Vertex>& aggregates)
2550 : graph_(graph), aggregatesBuilder_(aggregatesBuilder), aggregates_(aggregates), begin_(graph.begin()), end_(graph.end())
2556 Aggregator<G>::Stack::~Stack()
2567 inline typename G::VertexDescriptor Aggregator<G>::Stack::pop()
2573 typename G::VertexDescriptor current=*begin_;
2583 void printAggregates2d(
const AggregatesMap<V>& aggregates,
int n,
int m, std::ostream& os)
2587 std::ios_base::fmtflags oldOpts=os.flags();
2589 os.setf(std::ios_base::right, std::ios_base::adjustfield);
2594 for(
int i=0; i< n*m; i++)
2595 maxVal=
max(maxVal, aggregates[i]);
2597 for(
int i=10; i < 1000000; i*=10)
2603 for(
int j=0, entry=0; j < m; j++) {
2604 for(
int i=0; i<n; i++, entry++) {
2606 os<<aggregates[entry]<<
" ";
A class for temporarily storing the vertices of an aggregate in.
Definition: aggregates.hh:778
A Dummy visitor that does nothing for each visited edge.
Definition: aggregates.hh:598
Class providing information about the mapping of the vertices onto aggregates.
Definition: aggregates.hh:560
Base class of all aggregation criterions.
Definition: aggregates.hh:49
Class for building the aggregates.
Definition: aggregates.hh:909
Dependency policy for symmetric matrices.
Definition: aggregates.hh:253
Norm that uses only the [N][N] entry of the block to determine couplings.
Definition: aggregates.hh:379
Norm that uses only the [0][0] entry of the block to determine couplings.
Definition: aggregates.hh:455
Iterator over all edges starting from a vertex.
Definition: graph.hh:95
The vertex iterator type of the graph.
Definition: graph.hh:209
The (undirected) graph of a matrix.
Definition: graph.hh:51
M::size_type VertexDescriptor
The vertex descriptor.
Definition: graph.hh:73
EdgeIteratorT< const MatrixGraph< Matrix > > ConstEdgeIterator
The constant edge iterator type.
Definition: graph.hh:298
All parameters for AMG.
Definition: parameters.hh:416
Criterion taking advantage of symmetric matrices.
Definition: aggregates.hh:519
Dependency policy for symmetric matrices.
Definition: aggregates.hh:314
Dependency policy for symmetric matrices.
Definition: aggregates.hh:134
Criterion suitable for unsymmetric matrices.
Definition: aggregates.hh:539
derive error class from the base class in common
Definition: istlexception.hh:19
A generic dynamic dense matrix.
Definition: matrix.hh:561
typename Imp::BlockTraits< T >::field_type field_type
Export the type representing the underlying field.
Definition: matrix.hh:565
row_type::const_iterator ConstColIterator
Const iterator for the entries of each row.
Definition: matrix.hh:589
MatrixImp::DenseMatrixBase< T, A >::window_type row_type
The type implementing a matrix row.
Definition: matrix.hh:574
An allocator managing a pool of objects for reuse.
Definition: poolallocator.hh:223
A single linked list.
Definition: sllist.hh:44
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:74
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
constexpr GeometryType vertex
GeometryType representing a vertex.
Definition: type.hh:492
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:484
constexpr auto min
Function object that returns the smaller of the given values.
Definition: hybridutilities.hh:506
Matrix::ConstColIterator ColIter
Constant column iterator of the matrix.
Definition: aggregates.hh:273
Matrix::ConstColIterator ColIter
Constant column iterator of the matrix.
Definition: aggregates.hh:154
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:586
int id()
Get the id identifying the aggregate.
Norm norm_
The functor for calculating the norm.
Definition: aggregates.hh:302
MatrixGraph::VertexDescriptor Vertex
The vertex identifier.
Definition: aggregates.hh:920
AggregationCriterion()
Constructor.
Definition: aggregates.hh:66
const Matrix * matrix_
The matrix we work on.
Definition: aggregates.hh:357
auto operator()(const M &m, typename std::enable_if_t< Dune::IsNumber< M >::value > *sfinae=nullptr) const
Compute the norm of a scalar.
Definition: aggregates.hh:406
M Matrix
The matrix type we build the dependency of.
Definition: aggregates.hh:258
G MatrixGraph
The matrix graph type used.
Definition: aggregates.hh:915
Norm norm_
The functor for calculating the norm.
Definition: aggregates.hh:363
M Matrix
The matrix type we build the dependency of.
Definition: aggregates.hh:319
real_type diagonal_
The norm of the current diagonal.
Definition: aggregates.hh:187
N Norm
The norm to use for examining the matrix entries.
Definition: aggregates.hh:263
int row_
index of the currently evaluated row.
Definition: aggregates.hh:185
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:329
const Matrix * matrix_
The matrix we work on.
Definition: aggregates.hh:296
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:334
~AggregatesMap()
Destructor.
Matrix::field_type field_type
The current max value.
Definition: aggregates.hh:179
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:1064
std::size_t * SphereMap
Type of the mapping of aggregate members onto distance spheres.
Definition: aggregates.hh:809
Matrix::row_type Row
Constant Row iterator of the matrix.
Definition: aggregates.hh:268
VertexSet::size_type connectSize()
Get the number of connections to other aggregates.
N Norm
The norm to use for examining the matrix entries.
Definition: aggregates.hh:324
Norm norm_
The functor for calculating the norm.
Definition: aggregates.hh:183
VertexSet::const_iterator const_iterator
Const iterator over a vertex list.
Definition: aggregates.hh:804
MatrixGraph::VertexDescriptor AggregateDescriptor
The type of the aggregate descriptor.
Definition: aggregates.hh:923
real_type diagonal_
The norm of the current diagonal.
Definition: aggregates.hh:367
void add(const Vertex &vertex)
Add a vertex to the aggregate.
T DependencyPolicy
The policy for calculating the dependency graph.
Definition: aggregates.hh:55
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:506
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:82
V AggregateDescriptor
The aggregate descriptor type.
Definition: aggregates.hh:580
static const V ISOLATED
Identifier of isolated vertices.
Definition: aggregates.hh:571
int row_
index of the currently evaluated row.
Definition: aggregates.hh:365
DependencyCounter()
Constructor.
real_type diagonal_
The norm of the current diagonal.
Definition: aggregates.hh:306
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:105
AggregatesMap(std::size_t noVertices)
Constructs with allocating memory.
Matrix::field_type field_type
The current max value.
Definition: aggregates.hh:359
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:592
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:304
M Matrix
The matrix type we build the dependency of.
Definition: aggregates.hh:139
const Matrix * matrix_
The matrix we work on.
Definition: aggregates.hh:177
S VertexSet
The type of a single linked list of vertex descriptors.
Definition: aggregates.hh:801
FieldTraits< typenameM::field_type >::real_type operator()(const M &m) const
compute the norm of a matrix.
Definition: aggregates.hh:490
static const V UNAGGREGATED
Identifier of not yet aggregated vertices.
Definition: aggregates.hh:566
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:298
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:144
FieldTraits< typenameM::field_type >::real_type operator()(const M &m) const
compute the norm of a matrix.
Definition: aggregates.hh:473
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.
FieldTraits< typenameM::field_type >::real_type operator()(const M &m, typename std::enable_if_t<!Dune::IsNumber< M >::value > *sfinae=nullptr) const
compute the norm of a matrix.
Definition: aggregates.hh:390
MatrixGraph::VertexDescriptor Vertex
The vertex descriptor type.
Definition: aggregates.hh:789
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:575
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:149
PoolAllocator< Vertex, 100 > Allocator
The allocator we use for our lists and the set.
Definition: aggregates.hh:795
DVVerbType dvverb(std::cout)
stream for very verbose output.
Definition: stdstreams.hh:96
DInfoType dinfo(std::cout)
Stream for informative output.
Definition: stdstreams.hh:141
DVerbType dverb(std::cout)
Singleton of verbose debug stream.
Definition: stdstreams.hh:117
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
constexpr std::integral_constant< T, I0 > front(std::integer_sequence< T, I0, II... >)
Return the first entry of the sequence.
Definition: integersequence.hh:39
constexpr auto get(std::integer_sequence< T, II... >, std::integral_constant< std::size_t, pos >={})
Return the entry at position pos of the given sequence.
Definition: integersequence.hh:22
Parameter classes for customizing AMG.
An stl-compliant pool allocator.
Provides classes for handling internal properties in a graph.
Implements a scalar matrix view wrapper around an existing scalar.
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:463
Whether this type acts as a scalar in the context of (hierarchically blocked) containers.
Definition: typetraits.hh:194