dune-istl  2.8.0
aggregates.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #ifndef DUNE_AMG_AGGREGATES_HH
4 #define DUNE_AMG_AGGREGATES_HH
5 
6 
7 #include "parameters.hh"
8 #include "graph.hh"
9 #include "properties.hh"
10 #include "combinedfunctor.hh"
11 
12 #include <dune/common/timer.hh>
13 #include <dune/common/stdstreams.hh>
14 #include <dune/common/poolallocator.hh>
15 #include <dune/common/sllist.hh>
16 #include <dune/common/ftraits.hh>
17 #include <dune/common/scalarmatrixview.hh>
18 
19 #include <utility>
20 #include <set>
21 #include <algorithm>
22 #include <complex>
23 #include <limits>
24 #include <ostream>
25 #include <tuple>
26 
27 namespace Dune
28 {
29  namespace Amg
30  {
31 
45  template<class T>
46  class AggregationCriterion : public T
47  {
48 
49  public:
53  typedef T DependencyPolicy;
54 
65  : T()
66  {}
67 
69  : T(parms)
70  {}
80  void setDefaultValuesIsotropic(std::size_t dim, std::size_t diameter=2)
81  {
82  this->setMaxDistance(diameter-1);
83  std::size_t csize=1;
84 
85  for(; dim>0; dim--) {
86  csize*=diameter;
87  this->setMaxDistance(this->maxDistance()+diameter-1);
88  }
89  this->setMinAggregateSize(csize);
90  this->setMaxAggregateSize(static_cast<std::size_t>(csize*1.5));
91  }
92 
103  void setDefaultValuesAnisotropic(std::size_t dim,std::size_t diameter=2)
104  {
105  setDefaultValuesIsotropic(dim, diameter);
106  this->setMaxDistance(this->maxDistance()+dim-1);
107  }
108  };
109 
110  template<class T>
111  std::ostream& operator<<(std::ostream& os, const AggregationCriterion<T>& criterion)
112  {
113  os<<"{ maxdistance="<<criterion.maxDistance()<<" minAggregateSize="
114  <<criterion.minAggregateSize()<< " maxAggregateSize="<<criterion.maxAggregateSize()
115  <<" connectivity="<<criterion.maxConnectivity()<<" debugLevel="<<criterion.debugLevel()<<"}";
116  return os;
117  }
118 
130  template<class M, class N>
132  {
133  public:
137  typedef M Matrix;
138 
142  typedef N Norm;
143 
147  typedef typename Matrix::row_type Row;
148 
153 
154  void init(const Matrix* matrix);
155 
156  void initRow(const Row& row, int index);
157 
158  void examine(const ColIter& col);
159 
160  template<class G>
161  void examine(G& graph, const typename G::EdgeIterator& edge, const ColIter& col);
162 
163  bool isIsolated();
164 
165 
167  : Parameters(parms)
168  {}
170  : Parameters()
171  {}
172 
173  protected:
175  const Matrix* matrix_;
177  typedef typename Matrix::field_type field_type;
178  typedef typename FieldTraits<field_type>::real_type real_type;
183  int row_;
186  std::vector<real_type> vals_;
187  typename std::vector<real_type>::iterator valIter_;
188 
189  };
190 
191 
192  template<class M, class N>
193  inline void SymmetricMatrixDependency<M,N>::init(const Matrix* matrix)
194  {
195  matrix_ = matrix;
196  }
197 
198  template<class M, class N>
199  inline void SymmetricMatrixDependency<M,N>::initRow(const Row& row, int index)
200  {
201  using std::min;
202  vals_.assign(row.size(), 0.0);
203  assert(vals_.size()==row.size());
204  valIter_=vals_.begin();
205 
206  maxValue_ = min(- std::numeric_limits<real_type>::max(), std::numeric_limits<real_type>::min());
207  diagonal_=norm_(row[index]);
208  row_ = index;
209  }
210 
211  template<class M, class N>
213  {
214  using std::max;
215  // skip positive offdiagonals if norm preserves sign of them.
216  real_type eij = norm_(*col);
217  if(!N::is_sign_preserving || eij<0) // || eji<0)
218  {
219  *valIter_ = eij/diagonal_*eij/norm_(matrix_->operator[](col.index())[col.index()]);
220  maxValue_ = max(maxValue_, *valIter_);
221  }else
222  *valIter_ =0;
223  ++valIter_;
224  }
225 
226  template<class M, class N>
227  template<class G>
228  inline void SymmetricMatrixDependency<M,N>::examine(G&, const typename G::EdgeIterator& edge, const ColIter&)
229  {
230  if(*valIter_ > alpha() * maxValue_) {
231  edge.properties().setDepends();
232  edge.properties().setInfluences();
233  }
234  ++valIter_;
235  }
236 
237  template<class M, class N>
239  {
240  if(diagonal_==0)
241  DUNE_THROW(Dune::ISTLError, "No diagonal entry for row "<<row_<<".");
242  valIter_=vals_.begin();
243  return maxValue_ < beta();
244  }
245 
249  template<class M, class N>
250  class Dependency : public Parameters
251  {
252  public:
256  typedef M Matrix;
257 
261  typedef N Norm;
262 
266  typedef typename Matrix::row_type Row;
267 
272 
273  void init(const Matrix* matrix);
274 
275  void initRow(const Row& row, int index);
276 
277  void examine(const ColIter& col);
278 
279  template<class G>
280  void examine(G& graph, const typename G::EdgeIterator& edge, const ColIter& col);
281 
282  bool isIsolated();
283 
284  Dependency(const Parameters& parms)
285  : Parameters(parms)
286  {}
287 
289  : Parameters()
290  {}
291 
292  protected:
294  const Matrix* matrix_;
296  typedef typename Matrix::field_type field_type;
297  typedef typename FieldTraits<field_type>::real_type real_type;
302  int row_;
305  };
306 
310  template<class M, class N>
312  {
313  public:
317  typedef M Matrix;
318 
322  typedef N Norm;
323 
327  typedef typename Matrix::row_type Row;
328 
333 
334  void init(const Matrix* matrix);
335 
336  void initRow(const Row& row, int index);
337 
338  void examine(const ColIter& col);
339 
340  template<class G>
341  void examine(G& graph, const typename G::EdgeIterator& edge, const ColIter& col);
342 
343  bool isIsolated();
344 
345 
347  : Parameters(parms)
348  {}
350  : Parameters()
351  {}
352 
353  protected:
355  const Matrix* matrix_;
357  typedef typename Matrix::field_type field_type;
358  typedef typename FieldTraits<field_type>::real_type real_type;
363  int row_;
366  private:
367  void initRow(const Row& row, int index, const std::true_type&);
368  void initRow(const Row& row, int index, const std::false_type&);
369  };
370 
375  template<int N>
376  class Diagonal
377  {
378  public:
379  enum { /* @brief We preserve the sign.*/
380  is_sign_preserving = true
381  };
382 
387  template<class M>
388  typename FieldTraits<typename M::field_type>::real_type operator()(const M& m,
389  [[maybe_unused]] typename std::enable_if_t<!Dune::IsNumber<M>::value>* sfinae = nullptr) const
390  {
391  typedef typename M::field_type field_type;
392  typedef typename FieldTraits<field_type>::real_type real_type;
393  static_assert( std::is_convertible<field_type, real_type >::value,
394  "use of diagonal norm in AMG not implemented for complex field_type");
395  return m[N][N];
396  // possible implementation for complex types: return signed_abs(m[N][N]);
397  }
398 
403  template<class M>
404  auto operator()(const M& m,
405  typename std::enable_if_t<Dune::IsNumber<M>::value>* sfinae = nullptr) const
406  {
407  typedef typename FieldTraits<M>::real_type real_type;
408  static_assert( std::is_convertible<M, real_type >::value,
409  "use of diagonal norm in AMG not implemented for complex field_type");
410  return m;
411  // possible implementation for complex types: return signed_abs(m[N][N]);
412  }
413 
414  private:
415 
417  template<typename T>
418  static T signed_abs(const T & v)
419  {
420  return v;
421  }
422 
424  template<typename T>
425  static T signed_abs(const std::complex<T> & v)
426  {
427  // return sign * abs_value
428  // in case of complex numbers this extends to using the csgn function to determine the sign
429  return csgn(v) * std::abs(v);
430  }
431 
433  template<typename T>
434  static T csgn(const T & v)
435  {
436  return (T(0) < v) - (v < T(0));
437  }
438 
440  template<typename T>
441  static T csgn(std::complex<T> a)
442  {
443  return csgn(a.real())+(a.real() == 0.0)*csgn(a.imag());
444  }
445 
446  };
447 
452  class FirstDiagonal : public Diagonal<0>
453  {};
454 
460  struct RowSum
461  {
462 
463  enum { /* @brief We preserve the sign.*/
464  is_sign_preserving = false
465  };
470  template<class M>
471  typename FieldTraits<typename M::field_type>::real_type operator()(const M& m) const
472  {
473  return m.infinity_norm();
474  }
475  };
476 
478  {
479 
480  enum { /* @brief We preserve the sign.*/
481  is_sign_preserving = false
482  };
487  template<class M>
488  typename FieldTraits<typename M::field_type>::real_type operator()(const M& m) const
489  {
490  return m.frobenius_norm();
491  }
492  };
494  {
495 
496  enum { /* @brief We preserve the sign.*/
497  is_sign_preserving = false
498  };
503  template<class M>
504  typename FieldTraits<typename M::field_type>::real_type operator()(const M& m) const
505  {
506  return 1;
507  }
508  };
515  template<class M, class Norm>
516  class SymmetricCriterion : public AggregationCriterion<SymmetricDependency<M,Norm> >
517  {
518  public:
521  {}
523  {}
524  };
525 
526 
535  template<class M, class Norm>
536  class UnSymmetricCriterion : public AggregationCriterion<Dependency<M,Norm> >
537  {
538  public:
540  : AggregationCriterion<Dependency<M,Norm> >(parms)
541  {}
543  {}
544  };
545  // forward declaration
546  template<class G> class Aggregator;
547 
548 
556  template<class V>
558  {
559  public:
560 
564  static const V UNAGGREGATED;
565 
569  static const V ISOLATED;
573  typedef V VertexDescriptor;
574 
579 
584  typedef PoolAllocator<VertexDescriptor,100> Allocator;
585 
590  typedef SLList<VertexDescriptor,Allocator> VertexList;
591 
596  {
597  public:
598  template<class EdgeIterator>
599  void operator()([[maybe_unused]] const EdgeIterator& edge) const
600  {}
601  };
602 
603 
608 
614  AggregatesMap(std::size_t noVertices);
615 
620 
632  template<class M, class G, class C>
633  std::tuple<int,int,int,int> buildAggregates(const M& matrix, G& graph, const C& criterion,
634  bool finestLevel);
635 
653  template<bool reset, class G, class F, class VM>
654  std::size_t breadthFirstSearch(const VertexDescriptor& start,
655  const AggregateDescriptor& aggregate,
656  const G& graph,
657  F& aggregateVisitor,
658  VM& visitedMap) const;
659 
683  template<bool remove, bool reset, class G, class L, class F1, class F2, class VM>
684  std::size_t breadthFirstSearch(const VertexDescriptor& start,
685  const AggregateDescriptor& aggregate,
686  const G& graph, L& visited, F1& aggregateVisitor,
687  F2& nonAggregateVisitor,
688  VM& visitedMap) const;
689 
695  void allocate(std::size_t noVertices);
696 
700  std::size_t noVertices() const;
701 
705  void free();
706 
713 
720 
722 
724  {
725  return aggregates_;
726  }
727 
729  {
730  return aggregates_+noVertices();
731  }
732 
734 
736  {
737  return aggregates_;
738  }
739 
741  {
742  return aggregates_+noVertices();
743  }
744  private:
746  AggregatesMap(const AggregatesMap<V>&) = delete;
748  AggregatesMap<V>& operator=(const AggregatesMap<V>&) = delete;
749 
753  AggregateDescriptor* aggregates_;
754 
758  std::size_t noVertices_;
759  };
760 
764  template<class G, class C>
765  void buildDependency(G& graph,
766  const typename C::Matrix& matrix,
767  C criterion,
768  bool finestLevel);
769 
774  template<class G, class S>
775  class Aggregate
776  {
777 
778  public:
779 
780  /***
781  * @brief The type of the matrix graph we work with.
782  */
783  typedef G MatrixGraph;
788 
793  typedef PoolAllocator<Vertex,100> Allocator;
794 
799  typedef S VertexSet;
800 
802  typedef typename VertexSet::const_iterator const_iterator;
803 
807  typedef std::size_t* SphereMap;
808 
818  VertexSet& connectivity, std::vector<Vertex>& front_);
819 
820  void invalidate()
821  {
822  --id_;
823  }
824 
831  void reconstruct(const Vertex& vertex);
832 
836  void seed(const Vertex& vertex);
837 
841  void add(const Vertex& vertex);
842 
843  void add(std::vector<Vertex>& vertex);
847  void clear();
848 
852  typename VertexSet::size_type size();
856  typename VertexSet::size_type connectSize();
857 
861  int id();
862 
865 
868 
869  private:
873  VertexSet vertices_;
874 
879  int id_;
880 
884  MatrixGraph& graph_;
885 
889  AggregatesMap<Vertex>& aggregates_;
890 
894  VertexSet& connected_;
895 
899  std::vector<Vertex>& front_;
900  };
901 
905  template<class G>
907  {
908  public:
909 
913  typedef G MatrixGraph;
914 
919 
922 
927 
932 
949  template<class M, class C>
950  std::tuple<int,int,int,int> build(const M& m, G& graph,
951  AggregatesMap<Vertex>& aggregates, const C& c,
952  bool finestLevel);
953  private:
958  typedef PoolAllocator<Vertex,100> Allocator;
959 
963  typedef SLList<Vertex,Allocator> VertexList;
964 
968  typedef std::set<Vertex,std::less<Vertex>,Allocator> VertexSet;
969 
973  typedef std::size_t* SphereMap;
974 
978  MatrixGraph* graph_;
979 
984 
988  std::vector<Vertex> front_;
989 
993  VertexSet connected_;
994 
998  int size_;
999 
1003  class Stack
1004  {
1005  public:
1006  static const Vertex NullEntry;
1007 
1008  Stack(const MatrixGraph& graph,
1009  const Aggregator<G>& aggregatesBuilder,
1010  const AggregatesMap<Vertex>& aggregates);
1013  private:
1014  enum { N = 1300000 };
1015 
1017  const MatrixGraph& graph_;
1019  const Aggregator<G>& aggregatesBuilder_;
1021  const AggregatesMap<Vertex>& aggregates_;
1023  int size_;
1024  Vertex maxSize_;
1026  typename MatrixGraph::ConstVertexIterator begin_;
1027  typename MatrixGraph::ConstVertexIterator end_;
1028 
1030  Vertex* vals_;
1031 
1032  };
1033 
1034  friend class Stack;
1035 
1046  template<class V>
1047  void visitAggregateNeighbours(const Vertex& vertex, const AggregateDescriptor& aggregate,
1048  const AggregatesMap<Vertex>& aggregates,
1049  V& visitor) const;
1050 
1055  template<class V>
1056  class AggregateVisitor
1057  {
1058  public:
1062  typedef V Visitor;
1070  AggregateVisitor(const AggregatesMap<Vertex>& aggregates, const AggregateDescriptor& aggregate,
1071  Visitor& visitor);
1072 
1079  void operator()(const typename MatrixGraph::ConstEdgeIterator& edge);
1080 
1081  private:
1083  const AggregatesMap<Vertex>& aggregates_;
1085  AggregateDescriptor aggregate_;
1087  Visitor* visitor_;
1088  };
1089 
1093  class Counter
1094  {
1095  public:
1099  int value();
1100 
1101  protected:
1103  void increment();
1105  void decrement();
1106 
1107  private:
1108  int count_;
1109  };
1110 
1111 
1116  class FrontNeighbourCounter : public Counter
1117  {
1118  public:
1124 
1125  void operator()(const typename MatrixGraph::ConstEdgeIterator& edge);
1126 
1127  private:
1128  const MatrixGraph& graph_;
1129  };
1130 
1135  int noFrontNeighbours(const Vertex& vertex) const;
1136 
1140  class TwoWayCounter : public Counter
1141  {
1142  public:
1143  void operator()(const typename MatrixGraph::ConstEdgeIterator& edge);
1144  };
1145 
1157  int twoWayConnections(const Vertex&, const AggregateDescriptor& aggregate,
1158  const AggregatesMap<Vertex>& aggregates) const;
1159 
1163  class OneWayCounter : public Counter
1164  {
1165  public:
1166  void operator()(const typename MatrixGraph::ConstEdgeIterator& edge);
1167  };
1168 
1180  int oneWayConnections(const Vertex&, const AggregateDescriptor& aggregate,
1181  const AggregatesMap<Vertex>& aggregates) const;
1182 
1189  class ConnectivityCounter : public Counter
1190  {
1191  public:
1198  ConnectivityCounter(const VertexSet& connected, const AggregatesMap<Vertex>& aggregates);
1199 
1200  void operator()(const typename MatrixGraph::ConstEdgeIterator& edge);
1201 
1202  private:
1204  const VertexSet& connected_;
1206  const AggregatesMap<Vertex>& aggregates_;
1207 
1208  };
1209 
1221  double connectivity(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates) const;
1229  bool connected(const Vertex& vertex, const AggregateDescriptor& aggregate,
1230  const AggregatesMap<Vertex>& aggregates) const;
1231 
1239  bool connected(const Vertex& vertex, const SLList<AggregateDescriptor>& aggregateList,
1240  const AggregatesMap<Vertex>& aggregates) const;
1241 
1249  class DependencyCounter : public Counter
1250  {
1251  public:
1256 
1257  void operator()(const typename MatrixGraph::ConstEdgeIterator& edge);
1258  };
1259 
1266  class FrontMarker
1267  {
1268  public:
1275  FrontMarker(std::vector<Vertex>& front, MatrixGraph& graph);
1276 
1277  void operator()(const typename MatrixGraph::ConstEdgeIterator& edge);
1278 
1279  private:
1281  std::vector<Vertex>& front_;
1283  MatrixGraph& graph_;
1284  };
1285 
1289  void unmarkFront();
1290 
1305  int unusedNeighbours(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates) const;
1306 
1320  std::pair<int,int> neighbours(const Vertex& vertex,
1321  const AggregateDescriptor& aggregate,
1322  const AggregatesMap<Vertex>& aggregates) const;
1339  int aggregateNeighbours(const Vertex& vertex, const AggregateDescriptor& aggregate, const AggregatesMap<Vertex>& aggregates) const;
1340 
1348  bool admissible(const Vertex& vertex, const AggregateDescriptor& aggregate, const AggregatesMap<Vertex>& aggregates) const;
1349 
1357  std::size_t distance(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates);
1358 
1367  Vertex mergeNeighbour(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates) const;
1368 
1377  void nonisoNeighbourAggregate(const Vertex& vertex,
1378  const AggregatesMap<Vertex>& aggregates,
1379  SLList<Vertex>& neighbours) const;
1380 
1388  template<class C>
1389  void growAggregate(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates, const C& c);
1390  template<class C>
1391  void growIsolatedAggregate(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates, const C& c);
1392  };
1393 
1394 #ifndef DOXYGEN
1395 
1396  template<class M, class N>
1397  inline void SymmetricDependency<M,N>::init(const Matrix* matrix)
1398  {
1399  matrix_ = matrix;
1400  }
1401 
1402  template<class M, class N>
1403  inline void SymmetricDependency<M,N>::initRow(const Row& row, int index)
1404  {
1405  initRow(row, index, std::is_convertible<field_type, real_type>());
1406  }
1407 
1408  template<class M, class N>
1409  inline void SymmetricDependency<M,N>::initRow(const Row& row, int index, const std::false_type&)
1410  {
1411  DUNE_THROW(InvalidStateException, "field_type needs to convertible to real_type");
1412  }
1413 
1414  template<class M, class N>
1415  inline void SymmetricDependency<M,N>::initRow([[maybe_unused]] const Row& row, int index, const std::true_type&)
1416  {
1417  using std::min;
1418  maxValue_ = min(- std::numeric_limits<typename Matrix::field_type>::max(), std::numeric_limits<typename Matrix::field_type>::min());
1419  row_ = index;
1420  diagonal_ = norm_(matrix_->operator[](row_)[row_]);
1421  }
1422 
1423  template<class M, class N>
1424  inline void SymmetricDependency<M,N>::examine(const ColIter& col)
1425  {
1426  using std::max;
1427  real_type eij = norm_(*col);
1428  typename Matrix::ConstColIterator opposite_entry =
1429  matrix_->operator[](col.index()).find(row_);
1430  if ( opposite_entry == matrix_->operator[](col.index()).end() )
1431  {
1432  // Consider this a weak connection we disregard.
1433  return;
1434  }
1435  real_type eji = norm_(*opposite_entry);
1436 
1437  // skip positive offdiagonals if norm preserves sign of them.
1438  if(!N::is_sign_preserving || eij<0 || eji<0)
1439  maxValue_ = max(maxValue_,
1440  eij /diagonal_ * eji/
1441  norm_(matrix_->operator[](col.index())[col.index()]));
1442  }
1443 
1444  template<class M, class N>
1445  template<class G>
1446  inline void SymmetricDependency<M,N>::examine(G& graph, const typename G::EdgeIterator& edge, const ColIter& col)
1447  {
1448  real_type eij = norm_(*col);
1449  typename Matrix::ConstColIterator opposite_entry =
1450  matrix_->operator[](col.index()).find(row_);
1451 
1452  if ( opposite_entry == matrix_->operator[](col.index()).end() )
1453  {
1454  // Consider this as a weak connection we disregard.
1455  return;
1456  }
1457  real_type eji = norm_(*opposite_entry);
1458  // skip positve offdiagonals if norm preserves sign of them.
1459  if(!N::is_sign_preserving || (eij<0 || eji<0))
1460  if(eji / norm_(matrix_->operator[](edge.target())[edge.target()]) *
1461  eij/ diagonal_ > alpha() * maxValue_) {
1462  edge.properties().setDepends();
1463  edge.properties().setInfluences();
1464  typename G::EdgeProperties& other = graph.getEdgeProperties(edge.target(), edge.source());
1465  other.setInfluences();
1466  other.setDepends();
1467  }
1468  }
1469 
1470  template<class M, class N>
1472  {
1473  return maxValue_ < beta();
1474  }
1475 
1476 
1477  template<class M, class N>
1478  inline void Dependency<M,N>::init(const Matrix* matrix)
1479  {
1480  matrix_ = matrix;
1481  }
1482 
1483  template<class M, class N>
1484  inline void Dependency<M,N>::initRow([[maybe_unused]] const Row& row, int index)
1485  {
1486  using std::min;
1487  maxValue_ = min(- std::numeric_limits<real_type>::max(), std::numeric_limits<real_type>::min());
1488  row_ = index;
1489  diagonal_ = norm_(matrix_->operator[](row_)[row_]);
1490  }
1491 
1492  template<class M, class N>
1493  inline void Dependency<M,N>::examine(const ColIter& col)
1494  {
1495  using std::max;
1496  maxValue_ = max(maxValue_, -norm_(*col));
1497  }
1498 
1499  template<class M, class N>
1500  template<class G>
1501  inline void Dependency<M,N>::examine(G& graph, const typename G::EdgeIterator& edge, const ColIter& col)
1502  {
1503  if(-norm_(*col) >= maxValue_ * alpha()) {
1504  edge.properties().setDepends();
1505  typedef typename G::EdgeDescriptor ED;
1506  ED e= graph.findEdge(edge.target(), edge.source());
1507  if(e!=std::numeric_limits<ED>::max())
1508  {
1509  typename G::EdgeProperties& other = graph.getEdgeProperties(e);
1510  other.setInfluences();
1511  }
1512  }
1513  }
1514 
1515  template<class M, class N>
1516  inline bool Dependency<M,N>::isIsolated()
1517  {
1518  return maxValue_ < beta() * diagonal_;
1519  }
1520 
1521  template<class G,class S>
1522  Aggregate<G,S>::Aggregate(MatrixGraph& graph, AggregatesMap<Vertex>& aggregates,
1523  VertexSet& connected, std::vector<Vertex>& front)
1524  : vertices_(), id_(-1), graph_(graph), aggregates_(aggregates),
1525  connected_(connected), front_(front)
1526  {}
1527 
1528  template<class G,class S>
1529  void Aggregate<G,S>::reconstruct(const Vertex& vertex)
1530  {
1531  /*
1532  vertices_.push_back(vertex);
1533  typedef typename VertexList::const_iterator iterator;
1534  iterator begin = vertices_.begin();
1535  iterator end = vertices_.end();*/
1536  throw "Not yet implemented";
1537 
1538  // while(begin!=end){
1539  //for();
1540  // }
1541 
1542  }
1543 
1544  template<class G,class S>
1545  inline void Aggregate<G,S>::seed(const Vertex& vertex)
1546  {
1547  dvverb<<"Connected cleared"<<std::endl;
1548  connected_.clear();
1549  vertices_.clear();
1550  connected_.insert(vertex);
1551  dvverb << " Inserting "<<vertex<<" size="<<connected_.size();
1552  ++id_ ;
1553  add(vertex);
1554  }
1555 
1556 
1557  template<class G,class S>
1558  inline void Aggregate<G,S>::add(const Vertex& vertex)
1559  {
1560  vertices_.insert(vertex);
1561  aggregates_[vertex]=id_;
1562  if(front_.size())
1563  front_.erase(std::lower_bound(front_.begin(), front_.end(), vertex));
1564 
1565 
1566  typedef typename MatrixGraph::ConstEdgeIterator iterator;
1567  const iterator end = graph_.endEdges(vertex);
1568  for(iterator edge = graph_.beginEdges(vertex); edge != end; ++edge) {
1569  dvverb << " Inserting "<<aggregates_[edge.target()];
1570  connected_.insert(aggregates_[edge.target()]);
1571  dvverb <<" size="<<connected_.size();
1572  if(aggregates_[edge.target()]==AggregatesMap<Vertex>::UNAGGREGATED &&
1573  !graph_.getVertexProperties(edge.target()).front())
1574  {
1575  front_.push_back(edge.target());
1576  graph_.getVertexProperties(edge.target()).setFront();
1577  }
1578  }
1579  dvverb <<std::endl;
1580  std::sort(front_.begin(), front_.end());
1581  }
1582 
1583  template<class G,class S>
1584  inline void Aggregate<G,S>::add(std::vector<Vertex>& vertices)
1585  {
1586 #ifndef NDEBUG
1587  std::size_t oldsize = vertices_.size();
1588 #endif
1589  typedef typename std::vector<Vertex>::iterator Iterator;
1590 
1591  typedef typename VertexSet::iterator SIterator;
1592 
1593  SIterator pos=vertices_.begin();
1594  std::vector<Vertex> newFront;
1595  newFront.reserve(front_.capacity());
1596 
1597  std::set_difference(front_.begin(), front_.end(), vertices.begin(), vertices.end(),
1598  std::back_inserter(newFront));
1599  front_=newFront;
1600 
1601  for(Iterator vertex=vertices.begin(); vertex != vertices.end(); ++vertex)
1602  {
1603  pos=vertices_.insert(pos,*vertex);
1604  vertices_.insert(*vertex);
1605  graph_.getVertexProperties(*vertex).resetFront(); // Not a front node any more.
1606  aggregates_[*vertex]=id_;
1607 
1608  typedef typename MatrixGraph::ConstEdgeIterator iterator;
1609  const iterator end = graph_.endEdges(*vertex);
1610  for(iterator edge = graph_.beginEdges(*vertex); edge != end; ++edge) {
1611  dvverb << " Inserting "<<aggregates_[edge.target()];
1612  connected_.insert(aggregates_[edge.target()]);
1613  if(aggregates_[edge.target()]==AggregatesMap<Vertex>::UNAGGREGATED &&
1614  !graph_.getVertexProperties(edge.target()).front())
1615  {
1616  front_.push_back(edge.target());
1617  graph_.getVertexProperties(edge.target()).setFront();
1618  }
1619  dvverb <<" size="<<connected_.size();
1620  }
1621  dvverb <<std::endl;
1622  }
1623  std::sort(front_.begin(), front_.end());
1624  assert(oldsize+vertices.size()==vertices_.size());
1625  }
1626  template<class G,class S>
1627  inline void Aggregate<G,S>::clear()
1628  {
1629  vertices_.clear();
1630  connected_.clear();
1631  id_=-1;
1632  }
1633 
1634  template<class G,class S>
1635  inline typename Aggregate<G,S>::VertexSet::size_type
1637  {
1638  return vertices_.size();
1639  }
1640 
1641  template<class G,class S>
1642  inline typename Aggregate<G,S>::VertexSet::size_type
1644  {
1645  return connected_.size();
1646  }
1647 
1648  template<class G,class S>
1649  inline int Aggregate<G,S>::id()
1650  {
1651  return id_;
1652  }
1653 
1654  template<class G,class S>
1656  {
1657  return vertices_.begin();
1658  }
1659 
1660  template<class G,class S>
1661  inline typename Aggregate<G,S>::const_iterator Aggregate<G,S>::end() const
1662  {
1663  return vertices_.end();
1664  }
1665 
1666  template<class V>
1667  const V AggregatesMap<V>::UNAGGREGATED = std::numeric_limits<V>::max();
1668 
1669  template<class V>
1670  const V AggregatesMap<V>::ISOLATED = std::numeric_limits<V>::max()-1;
1671 
1672  template<class V>
1674  : aggregates_(0)
1675  {}
1676 
1677  template<class V>
1679  {
1680  if(aggregates_!=0)
1681  delete[] aggregates_;
1682  }
1683 
1684 
1685  template<class V>
1686  inline AggregatesMap<V>::AggregatesMap(std::size_t noVertices)
1687  {
1688  allocate(noVertices);
1689  }
1690 
1691  template<class V>
1692  inline std::size_t AggregatesMap<V>::noVertices() const
1693  {
1694  return noVertices_;
1695  }
1696 
1697  template<class V>
1698  inline void AggregatesMap<V>::allocate(std::size_t noVertices)
1699  {
1700  aggregates_ = new AggregateDescriptor[noVertices];
1701  noVertices_ = noVertices;
1702 
1703  for(std::size_t i=0; i < noVertices; i++)
1704  aggregates_[i]=UNAGGREGATED;
1705  }
1706 
1707  template<class V>
1708  inline void AggregatesMap<V>::free()
1709  {
1710  assert(aggregates_ != 0);
1711  delete[] aggregates_;
1712  aggregates_=0;
1713  }
1714 
1715  template<class V>
1716  inline typename AggregatesMap<V>::AggregateDescriptor&
1717  AggregatesMap<V>::operator[](const VertexDescriptor& v)
1718  {
1719  return aggregates_[v];
1720  }
1721 
1722  template<class V>
1723  inline const typename AggregatesMap<V>::AggregateDescriptor&
1724  AggregatesMap<V>::operator[](const VertexDescriptor& v) const
1725  {
1726  return aggregates_[v];
1727  }
1728 
1729  template<class V>
1730  template<bool reset, class G, class F,class VM>
1731  inline std::size_t AggregatesMap<V>::breadthFirstSearch(const V& start,
1732  const AggregateDescriptor& aggregate,
1733  const G& graph, F& aggregateVisitor,
1734  VM& visitedMap) const
1735  {
1736  VertexList vlist;
1737 
1738  DummyEdgeVisitor dummy;
1739  return breadthFirstSearch<true,reset>(start, aggregate, graph, vlist, aggregateVisitor, dummy, visitedMap);
1740  }
1741 
1742  template<class V>
1743  template<bool remove, bool reset, class G, class L, class F1, class F2, class VM>
1744  std::size_t AggregatesMap<V>::breadthFirstSearch(const V& start,
1745  const AggregateDescriptor& aggregate,
1746  const G& graph,
1747  L& visited,
1748  F1& aggregateVisitor,
1749  F2& nonAggregateVisitor,
1750  VM& visitedMap) const
1751  {
1752  typedef typename L::const_iterator ListIterator;
1753  int visitedSpheres = 0;
1754 
1755  visited.push_back(start);
1756  put(visitedMap, start, true);
1757 
1758  ListIterator current = visited.begin();
1759  ListIterator end = visited.end();
1760  std::size_t i=0, size=visited.size();
1761 
1762  // visit the neighbours of all vertices of the
1763  // current sphere.
1764  while(current != end) {
1765 
1766  for(; i<size; ++current, ++i) {
1767  typedef typename G::ConstEdgeIterator EdgeIterator;
1768  const EdgeIterator endEdge = graph.endEdges(*current);
1769 
1770  for(EdgeIterator edge = graph.beginEdges(*current);
1771  edge != endEdge; ++edge) {
1772 
1773  if(aggregates_[edge.target()]==aggregate) {
1774  if(!get(visitedMap, edge.target())) {
1775  put(visitedMap, edge.target(), true);
1776  visited.push_back(edge.target());
1777  aggregateVisitor(edge);
1778  }
1779  }else
1780  nonAggregateVisitor(edge);
1781  }
1782  }
1783  end = visited.end();
1784  size = visited.size();
1785  if(current != end)
1786  visitedSpheres++;
1787  }
1788 
1789  if(reset)
1790  for(current = visited.begin(); current != end; ++current)
1791  put(visitedMap, *current, false);
1792 
1793 
1794  if(remove)
1795  visited.clear();
1796 
1797  return visitedSpheres;
1798  }
1799 
1800  template<class G>
1802  : graph_(0), aggregate_(0), front_(), connected_(), size_(-1)
1803  {}
1804 
1805  template<class G>
1807  {
1808  size_=-1;
1809  }
1810 
1811  template<class G, class C>
1812  void buildDependency(G& graph,
1813  const typename C::Matrix& matrix,
1814  C criterion, bool firstlevel)
1815  {
1816  // assert(graph.isBuilt());
1817  typedef typename C::Matrix Matrix;
1818  typedef typename G::VertexIterator VertexIterator;
1819 
1820  criterion.init(&matrix);
1821 
1822  for(VertexIterator vertex = graph.begin(); vertex != graph.end(); ++vertex) {
1823  typedef typename Matrix::row_type Row;
1824 
1825  const Row& row = matrix[*vertex];
1826 
1827  // Tell the criterion what row we will examine now
1828  // This might for example be used for calculating the
1829  // maximum offdiagonal value
1830  criterion.initRow(row, *vertex);
1831 
1832  // On a first path all columns are examined. After this
1833  // the calculator should know whether the vertex is isolated.
1834  typedef typename Matrix::ConstColIterator ColIterator;
1835  ColIterator end = row.end();
1836  typename FieldTraits<typename Matrix::field_type>::real_type absoffdiag=0.;
1837 
1838  using std::max;
1839  if(firstlevel) {
1840  for(ColIterator col = row.begin(); col != end; ++col)
1841  if(col.index()!=*vertex) {
1842  criterion.examine(col);
1843  absoffdiag = max(absoffdiag, Impl::asMatrix(*col).frobenius_norm());
1844  }
1845 
1846  if(absoffdiag==0)
1847  vertex.properties().setExcludedBorder();
1848  }
1849  else
1850  for(ColIterator col = row.begin(); col != end; ++col)
1851  if(col.index()!=*vertex)
1852  criterion.examine(col);
1853 
1854  // reset the vertex properties
1855  //vertex.properties().reset();
1856 
1857  // Check whether the vertex is isolated.
1858  if(criterion.isIsolated()) {
1859  //std::cout<<"ISOLATED: "<<*vertex<<std::endl;
1860  vertex.properties().setIsolated();
1861  }else{
1862  // Examine all the edges beginning at this vertex.
1863  auto eEnd = vertex.end();
1864  auto col = matrix[*vertex].begin();
1865 
1866  for(auto edge = vertex.begin(); edge!= eEnd; ++edge, ++col) {
1867  // Move to the right column.
1868  while(col.index()!=edge.target())
1869  ++col;
1870  criterion.examine(graph, edge, col);
1871  }
1872  }
1873 
1874  }
1875  }
1876 
1877 
1878  template<class G>
1879  template<class V>
1880  inline Aggregator<G>::AggregateVisitor<V>::AggregateVisitor(const AggregatesMap<Vertex>& aggregates,
1881  const AggregateDescriptor& aggregate, V& visitor)
1882  : aggregates_(aggregates), aggregate_(aggregate), visitor_(&visitor)
1883  {}
1884 
1885  template<class G>
1886  template<class V>
1887  inline void Aggregator<G>::AggregateVisitor<V>::operator()(const typename MatrixGraph::ConstEdgeIterator& edge)
1888  {
1889  if(aggregates_[edge.target()]==aggregate_)
1890  visitor_->operator()(edge);
1891  }
1892 
1893  template<class G>
1894  template<class V>
1895  inline void Aggregator<G>::visitAggregateNeighbours(const Vertex& vertex,
1896  const AggregateDescriptor& aggregate,
1897  const AggregatesMap<Vertex>& aggregates,
1898  V& visitor) const
1899  {
1900  // Only evaluates for edge pointing to the aggregate
1901  AggregateVisitor<V> v(aggregates, aggregate, visitor);
1902  visitNeighbours(*graph_, vertex, v);
1903  }
1904 
1905 
1906  template<class G>
1907  inline Aggregator<G>::Counter::Counter()
1908  : count_(0)
1909  {}
1910 
1911  template<class G>
1912  inline void Aggregator<G>::Counter::increment()
1913  {
1914  ++count_;
1915  }
1916 
1917  template<class G>
1918  inline void Aggregator<G>::Counter::decrement()
1919  {
1920  --count_;
1921  }
1922  template<class G>
1923  inline int Aggregator<G>::Counter::value()
1924  {
1925  return count_;
1926  }
1927 
1928  template<class G>
1929  inline void Aggregator<G>::TwoWayCounter::operator()(const typename MatrixGraph::ConstEdgeIterator& edge)
1930  {
1931  if(edge.properties().isTwoWay())
1932  Counter::increment();
1933  }
1934 
1935  template<class G>
1936  int Aggregator<G>::twoWayConnections(const Vertex& vertex, const AggregateDescriptor& aggregate,
1937  const AggregatesMap<Vertex>& aggregates) const
1938  {
1939  TwoWayCounter counter;
1940  visitAggregateNeighbours(vertex, aggregate, aggregates, counter);
1941  return counter.value();
1942  }
1943 
1944  template<class G>
1945  int Aggregator<G>::oneWayConnections(const Vertex& vertex, const AggregateDescriptor& aggregate,
1946  const AggregatesMap<Vertex>& aggregates) const
1947  {
1948  OneWayCounter counter;
1949  visitAggregateNeighbours(vertex, aggregate, aggregates, counter);
1950  return counter.value();
1951  }
1952 
1953  template<class G>
1954  inline void Aggregator<G>::OneWayCounter::operator()(const typename MatrixGraph::ConstEdgeIterator& edge)
1955  {
1956  if(edge.properties().isOneWay())
1957  Counter::increment();
1958  }
1959 
1960  template<class G>
1961  inline Aggregator<G>::ConnectivityCounter::ConnectivityCounter(const VertexSet& connected,
1962  const AggregatesMap<Vertex>& aggregates)
1963  : Counter(), connected_(connected), aggregates_(aggregates)
1964  {}
1965 
1966 
1967  template<class G>
1968  inline void Aggregator<G>::ConnectivityCounter::operator()(const typename MatrixGraph::ConstEdgeIterator& edge)
1969  {
1970  if(connected_.find(aggregates_[edge.target()]) == connected_.end() || aggregates_[edge.target()]==AggregatesMap<Vertex>::UNAGGREGATED)
1971  // Would be a new connection
1972  Counter::increment();
1973  else{
1974  Counter::increment();
1975  Counter::increment();
1976  }
1977  }
1978 
1979  template<class G>
1980  inline double Aggregator<G>::connectivity(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates) const
1981  {
1982  ConnectivityCounter counter(connected_, aggregates);
1983  double noNeighbours=visitNeighbours(*graph_, vertex, counter);
1984  return (double)counter.value()/noNeighbours;
1985  }
1986 
1987  template<class G>
1988  inline Aggregator<G>::DependencyCounter::DependencyCounter()
1989  : Counter()
1990  {}
1991 
1992  template<class G>
1993  inline void Aggregator<G>::DependencyCounter::operator()(const typename MatrixGraph::ConstEdgeIterator& edge)
1994  {
1995  if(edge.properties().depends())
1996  Counter::increment();
1997  if(edge.properties().influences())
1998  Counter::increment();
1999  }
2000 
2001  template<class G>
2002  int Aggregator<G>::unusedNeighbours(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates) const
2003  {
2004  return aggregateNeighbours(vertex, AggregatesMap<Vertex>::UNAGGREGATED, aggregates);
2005  }
2006 
2007  template<class G>
2008  std::pair<int,int> Aggregator<G>::neighbours(const Vertex& vertex,
2009  const AggregateDescriptor& aggregate,
2010  const AggregatesMap<Vertex>& aggregates) const
2011  {
2012  DependencyCounter unused, aggregated;
2013  typedef AggregateVisitor<DependencyCounter> CounterT;
2014  typedef std::tuple<CounterT,CounterT> CounterTuple;
2015  CombinedFunctor<CounterTuple> visitors(CounterTuple(CounterT(aggregates, AggregatesMap<Vertex>::UNAGGREGATED, unused), CounterT(aggregates, aggregate, aggregated)));
2016  visitNeighbours(*graph_, vertex, visitors);
2017  return std::make_pair(unused.value(), aggregated.value());
2018  }
2019 
2020 
2021  template<class G>
2022  int Aggregator<G>::aggregateNeighbours(const Vertex& vertex, const AggregateDescriptor& aggregate, const AggregatesMap<Vertex>& aggregates) const
2023  {
2024  DependencyCounter counter;
2025  visitAggregateNeighbours(vertex, aggregate, aggregates, counter);
2026  return counter.value();
2027  }
2028 
2029  template<class G>
2030  std::size_t Aggregator<G>::distance(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates)
2031  {
2032  return 0;
2033  typename PropertyMapTypeSelector<VertexVisitedTag,G>::Type visitedMap = get(VertexVisitedTag(), *graph_);
2034  VertexList vlist;
2035  typename AggregatesMap<Vertex>::DummyEdgeVisitor dummy;
2036  return aggregates.template breadthFirstSearch<true,true>(vertex,
2037  aggregate_->id(), *graph_,
2038  vlist, dummy, dummy, visitedMap);
2039  }
2040 
2041  template<class G>
2042  inline Aggregator<G>::FrontMarker::FrontMarker(std::vector<Vertex>& front, MatrixGraph& graph)
2043  : front_(front), graph_(graph)
2044  {}
2045 
2046  template<class G>
2047  inline void Aggregator<G>::FrontMarker::operator()(const typename MatrixGraph::ConstEdgeIterator& edge)
2048  {
2049  Vertex target = edge.target();
2050 
2051  if(!graph_.getVertexProperties(target).front()) {
2052  front_.push_back(target);
2053  graph_.getVertexProperties(target).setFront();
2054  }
2055  }
2056 
2057  template<class G>
2058  inline bool Aggregator<G>::admissible(const Vertex& vertex, const AggregateDescriptor& aggregate, const AggregatesMap<Vertex>& aggregates) const
2059  {
2060  // Todo
2061  Dune::dvverb<<" Admissible not yet implemented!"<<std::endl;
2062  return true;
2063  //Situation 1: front node depends on two nodes. Then these
2064  // have to be strongly connected to each other
2065 
2066  // Iterate over all all neighbours of front node
2067  typedef typename MatrixGraph::ConstEdgeIterator Iterator;
2068  Iterator vend = graph_->endEdges(vertex);
2069  for(Iterator edge = graph_->beginEdges(vertex); edge != vend; ++edge) {
2070  // if(edge.properties().depends() && !edge.properties().influences()
2071  if(edge.properties().isStrong()
2072  && aggregates[edge.target()]==aggregate)
2073  {
2074  // Search for another link to the aggregate
2075  Iterator edge1 = edge;
2076  for(++edge1; edge1 != vend; ++edge1) {
2077  //if(edge1.properties().depends() && !edge1.properties().influences()
2078  if(edge1.properties().isStrong()
2079  && aggregates[edge.target()]==aggregate)
2080  {
2081  //Search for an edge connecting the two vertices that is
2082  //strong
2083  bool found=false;
2084  Iterator v2end = graph_->endEdges(edge.target());
2085  for(Iterator edge2 = graph_->beginEdges(edge.target()); edge2 != v2end; ++edge2) {
2086  if(edge2.target()==edge1.target() &&
2087  edge2.properties().isStrong()) {
2088  found =true;
2089  break;
2090  }
2091  }
2092  if(found)
2093  {
2094  return true;
2095  }
2096  }
2097  }
2098  }
2099  }
2100 
2101  // Situation 2: cluster node depends on front node and other cluster node
2103  vend = graph_->endEdges(vertex);
2104  for(Iterator edge = graph_->beginEdges(vertex); edge != vend; ++edge) {
2105  //if(!edge.properties().depends() && edge.properties().influences()
2106  if(edge.properties().isStrong()
2107  && aggregates[edge.target()]==aggregate)
2108  {
2109  // Search for a link from target that stays within the aggregate
2110  Iterator v1end = graph_->endEdges(edge.target());
2111 
2112  for(Iterator edge1=graph_->beginEdges(edge.target()); edge1 != v1end; ++edge1) {
2113  //if(edge1.properties().depends() && !edge1.properties().influences()
2114  if(edge1.properties().isStrong()
2115  && aggregates[edge1.target()]==aggregate)
2116  {
2117  bool found=false;
2118  // Check if front node is also connected to this one
2119  Iterator v2end = graph_->endEdges(vertex);
2120  for(Iterator edge2 = graph_->beginEdges(vertex); edge2 != v2end; ++edge2) {
2121  if(edge2.target()==edge1.target()) {
2122  if(edge2.properties().isStrong())
2123  found=true;
2124  break;
2125  }
2126  }
2127  if(found)
2128  {
2129  return true;
2130  }
2131  }
2132  }
2133  }
2134  }
2135  return false;
2136  }
2137 
2138  template<class G>
2139  void Aggregator<G>::unmarkFront()
2140  {
2141  typedef typename std::vector<Vertex>::const_iterator Iterator;
2142 
2143  for(Iterator vertex=front_.begin(); vertex != front_.end(); ++vertex)
2144  graph_->getVertexProperties(*vertex).resetFront();
2145 
2146  front_.clear();
2147  }
2148 
2149  template<class G>
2150  inline void
2151  Aggregator<G>::nonisoNeighbourAggregate(const Vertex& vertex,
2152  const AggregatesMap<Vertex>& aggregates,
2153  SLList<Vertex>& neighbours) const
2154  {
2155  typedef typename MatrixGraph::ConstEdgeIterator Iterator;
2156  Iterator end=graph_->beginEdges(vertex);
2157  neighbours.clear();
2158 
2159  for(Iterator edge=graph_->beginEdges(vertex); edge!=end; ++edge)
2160  {
2161  if(aggregates[edge.target()]!=AggregatesMap<Vertex>::UNAGGREGATED && graph_->getVertexProperties(edge.target()).isolated())
2162  neighbours.push_back(aggregates[edge.target()]);
2163  }
2164  }
2165 
2166  template<class G>
2167  inline typename G::VertexDescriptor Aggregator<G>::mergeNeighbour(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates) const
2168  {
2169  typedef typename MatrixGraph::ConstEdgeIterator Iterator;
2170 
2171  Iterator end = graph_->endEdges(vertex);
2172  for(Iterator edge = graph_->beginEdges(vertex); edge != end; ++edge) {
2173  if(aggregates[edge.target()] != AggregatesMap<Vertex>::UNAGGREGATED &&
2174  graph_->getVertexProperties(edge.target()).isolated() == graph_->getVertexProperties(edge.source()).isolated()) {
2175  if( graph_->getVertexProperties(vertex).isolated() ||
2176  ((edge.properties().depends() || edge.properties().influences())
2177  && admissible(vertex, aggregates[edge.target()], aggregates)))
2178  return edge.target();
2179  }
2180  }
2182  }
2183 
2184  template<class G>
2185  Aggregator<G>::FrontNeighbourCounter::FrontNeighbourCounter(const MatrixGraph& graph)
2186  : Counter(), graph_(graph)
2187  {}
2188 
2189  template<class G>
2190  void Aggregator<G>::FrontNeighbourCounter::operator()(const typename MatrixGraph::ConstEdgeIterator& edge)
2191  {
2192  if(graph_.getVertexProperties(edge.target()).front())
2193  Counter::increment();
2194  }
2195 
2196  template<class G>
2197  int Aggregator<G>::noFrontNeighbours(const Vertex& vertex) const
2198  {
2199  FrontNeighbourCounter counter(*graph_);
2200  visitNeighbours(*graph_, vertex, counter);
2201  return counter.value();
2202  }
2203  template<class G>
2204  inline bool Aggregator<G>::connected(const Vertex& vertex,
2205  const AggregateDescriptor& aggregate,
2206  const AggregatesMap<Vertex>& aggregates) const
2207  {
2208  typedef typename G::ConstEdgeIterator iterator;
2209  const iterator end = graph_->endEdges(vertex);
2210  for(iterator edge = graph_->beginEdges(vertex); edge != end; ++edge)
2211  if(aggregates[edge.target()]==aggregate)
2212  return true;
2213  return false;
2214  }
2215  template<class G>
2216  inline bool Aggregator<G>::connected(const Vertex& vertex,
2217  const SLList<AggregateDescriptor>& aggregateList,
2218  const AggregatesMap<Vertex>& aggregates) const
2219  {
2220  typedef typename SLList<AggregateDescriptor>::const_iterator Iter;
2221  for(Iter i=aggregateList.begin(); i!=aggregateList.end(); ++i)
2222  if(connected(vertex, *i, aggregates))
2223  return true;
2224  return false;
2225  }
2226 
2227  template<class G>
2228  template<class C>
2229  void Aggregator<G>::growIsolatedAggregate(const Vertex& seed, const AggregatesMap<Vertex>& aggregates, const C& c)
2230  {
2231  SLList<Vertex> connectedAggregates;
2232  nonisoNeighbourAggregate(seed, aggregates,connectedAggregates);
2233 
2234  while(aggregate_->size()< c.minAggregateSize() && aggregate_->connectSize() < c.maxConnectivity()) {
2235  double maxCon=-1;
2236  std::size_t maxFrontNeighbours=0;
2237 
2239 
2240  typedef typename std::vector<Vertex>::const_iterator Iterator;
2241 
2242  for(Iterator vertex = front_.begin(); vertex != front_.end(); ++vertex) {
2243  if(distance(*vertex, aggregates)>c.maxDistance())
2244  continue; // distance of proposes aggregate too big
2245 
2246  if(connectedAggregates.size()>0) {
2247  // there is already a neighbour cluster
2248  // front node must be connected to same neighbour cluster
2249 
2250  if(!connected(*vertex, connectedAggregates, aggregates))
2251  continue;
2252  }
2253 
2254  double con = connectivity(*vertex, aggregates);
2255 
2256  if(con == maxCon) {
2257  std::size_t frontNeighbours = noFrontNeighbours(*vertex);
2258 
2259  if(frontNeighbours >= maxFrontNeighbours) {
2260  maxFrontNeighbours = frontNeighbours;
2261  candidate = *vertex;
2262  }
2263  }else if(con > maxCon) {
2264  maxCon = con;
2265  maxFrontNeighbours = noFrontNeighbours(*vertex);
2266  candidate = *vertex;
2267  }
2268  }
2269 
2271  break;
2272 
2273  aggregate_->add(candidate);
2274  }
2275  }
2276 
2277  template<class G>
2278  template<class C>
2279  void Aggregator<G>::growAggregate(const Vertex& seed, const AggregatesMap<Vertex>& aggregates, const C& c)
2280  {
2281  using std::min;
2282 
2283  std::size_t distance_ =0;
2284  while(aggregate_->size() < c.minAggregateSize()&& distance_<c.maxDistance()) {
2285  int maxTwoCons=0, maxOneCons=0, maxNeighbours=-1;
2286  double maxCon=-1;
2287 
2288  std::vector<Vertex> candidates;
2289  candidates.reserve(30);
2290 
2291  typedef typename std::vector<Vertex>::const_iterator Iterator;
2292 
2293  for(Iterator vertex = front_.begin(); vertex != front_.end(); ++vertex) {
2294  // Only nonisolated nodes are considered
2295  if(graph_->getVertexProperties(*vertex).isolated())
2296  continue;
2297 
2298  int twoWayCons = twoWayConnections(*vertex, aggregate_->id(), aggregates);
2299 
2300  /* The case of two way connections. */
2301  if( maxTwoCons == twoWayCons && twoWayCons > 0) {
2302  double con = connectivity(*vertex, aggregates);
2303 
2304  if(con == maxCon) {
2305  int neighbours = noFrontNeighbours(*vertex);
2306 
2307  if(neighbours > maxNeighbours) {
2308  maxNeighbours = neighbours;
2309  candidates.clear();
2310  candidates.push_back(*vertex);
2311  }else{
2312  candidates.push_back(*vertex);
2313  }
2314  }else if( con > maxCon) {
2315  maxCon = con;
2316  maxNeighbours = noFrontNeighbours(*vertex);
2317  candidates.clear();
2318  candidates.push_back(*vertex);
2319  }
2320  }else if(twoWayCons > maxTwoCons) {
2321  maxTwoCons = twoWayCons;
2322  maxCon = connectivity(*vertex, aggregates);
2323  maxNeighbours = noFrontNeighbours(*vertex);
2324  candidates.clear();
2325  candidates.push_back(*vertex);
2326 
2327  // two way connections precede
2328  maxOneCons = std::numeric_limits<int>::max();
2329  }
2330 
2331  if(twoWayCons > 0)
2332  {
2333  continue; // THis is a two-way node, skip tests for one way nodes
2334  }
2335 
2336  /* The one way case */
2337  int oneWayCons = oneWayConnections(*vertex, aggregate_->id(), aggregates);
2338 
2339  if(oneWayCons==0)
2340  continue; // No strong connections, skip the tests.
2341 
2342  if(!admissible(*vertex, aggregate_->id(), aggregates))
2343  continue;
2344 
2345  if( maxOneCons == oneWayCons && oneWayCons > 0) {
2346  double con = connectivity(*vertex, aggregates);
2347 
2348  if(con == maxCon) {
2349  int neighbours = noFrontNeighbours(*vertex);
2350 
2351  if(neighbours > maxNeighbours) {
2352  maxNeighbours = neighbours;
2353  candidates.clear();
2354  candidates.push_back(*vertex);
2355  }else{
2356  if(neighbours==maxNeighbours)
2357  {
2358  candidates.push_back(*vertex);
2359  }
2360  }
2361  }else if( con > maxCon) {
2362  maxCon = con;
2363  maxNeighbours = noFrontNeighbours(*vertex);
2364  candidates.clear();
2365  candidates.push_back(*vertex);
2366  }
2367  }else if(oneWayCons > maxOneCons) {
2368  maxOneCons = oneWayCons;
2369  maxCon = connectivity(*vertex, aggregates);
2370  maxNeighbours = noFrontNeighbours(*vertex);
2371  candidates.clear();
2372  candidates.push_back(*vertex);
2373  }
2374  }
2375 
2376 
2377  if(!candidates.size())
2378  break; // No more candidates found
2379  distance_=distance(seed, aggregates);
2380  candidates.resize(min(candidates.size(), c.maxAggregateSize()-
2381  aggregate_->size()));
2382  aggregate_->add(candidates);
2383  }
2384  }
2385 
2386  template<typename V>
2387  template<typename M, typename G, typename C>
2388  std::tuple<int,int,int,int> AggregatesMap<V>::buildAggregates(const M& matrix, G& graph, const C& criterion,
2389  bool finestLevel)
2390  {
2391  Aggregator<G> aggregator;
2392  return aggregator.build(matrix, graph, *this, criterion, finestLevel);
2393  }
2394 
2395  template<class G>
2396  template<class M, class C>
2397  std::tuple<int,int,int,int> Aggregator<G>::build(const M& m, G& graph, AggregatesMap<Vertex>& aggregates, const C& c,
2398  bool finestLevel)
2399  {
2400  using std::max;
2401  using std::min;
2402  // Stack for fast vertex access
2403  Stack stack_(graph, *this, aggregates);
2404 
2405  graph_ = &graph;
2406 
2407  aggregate_ = new Aggregate<G,VertexSet>(graph, aggregates, connected_, front_);
2408 
2409  Timer watch;
2410  watch.reset();
2411 
2412  buildDependency(graph, m, c, finestLevel);
2413 
2414  dverb<<"Build dependency took "<< watch.elapsed()<<" seconds."<<std::endl;
2415  int noAggregates, conAggregates, isoAggregates, oneAggregates;
2416  std::size_t maxA=0, minA=1000000, avg=0;
2417  int skippedAggregates;
2418  noAggregates = conAggregates = isoAggregates = oneAggregates =
2419  skippedAggregates = 0;
2420 
2421  while(true) {
2422  Vertex seed = stack_.pop();
2423 
2424  if(seed == Stack::NullEntry)
2425  // No more unaggregated vertices. We are finished!
2426  break;
2427 
2428  // Debugging output
2429  if((noAggregates+1)%10000 == 0)
2430  Dune::dverb<<"c";
2431  unmarkFront();
2432 
2433  if(graph.getVertexProperties(seed).excludedBorder()) {
2434  aggregates[seed]=AggregatesMap<Vertex>::ISOLATED;
2435  ++skippedAggregates;
2436  continue;
2437  }
2438 
2439  if(graph.getVertexProperties(seed).isolated()) {
2440  if(c.skipIsolated()) {
2441  // isolated vertices are not aggregated but skipped on the coarser levels.
2442  aggregates[seed]=AggregatesMap<Vertex>::ISOLATED;
2443  ++skippedAggregates;
2444  // skip rest as no agglomeration is done.
2445  continue;
2446  }else{
2447  aggregate_->seed(seed);
2448  growIsolatedAggregate(seed, aggregates, c);
2449  }
2450  }else{
2451  aggregate_->seed(seed);
2452  growAggregate(seed, aggregates, c);
2453  }
2454 
2455  /* The rounding step. */
2456  while(!(graph.getVertexProperties(seed).isolated()) && aggregate_->size() < c.maxAggregateSize()) {
2457 
2458  std::vector<Vertex> candidates;
2459  candidates.reserve(30);
2460 
2461  typedef typename std::vector<Vertex>::const_iterator Iterator;
2462 
2463  for(Iterator vertex = front_.begin(); vertex != front_.end(); ++vertex) {
2464 
2465  if(graph.getVertexProperties(*vertex).isolated())
2466  continue; // No isolated nodes here
2467 
2468  if(twoWayConnections( *vertex, aggregate_->id(), aggregates) == 0 &&
2469  (oneWayConnections( *vertex, aggregate_->id(), aggregates) == 0 ||
2470  !admissible( *vertex, aggregate_->id(), aggregates) ))
2471  continue;
2472 
2473  std::pair<int,int> neighbourPair=neighbours(*vertex, aggregate_->id(),
2474  aggregates);
2475 
2476  //if(aggregateNeighbours(*vertex, aggregate_->id(), aggregates) <= unusedNeighbours(*vertex, aggregates))
2477  // continue;
2478 
2479  if(neighbourPair.first >= neighbourPair.second)
2480  continue;
2481 
2482  if(distance(*vertex, aggregates) > c.maxDistance())
2483  continue; // Distance too far
2484  candidates.push_back(*vertex);
2485  break;
2486  }
2487 
2488  if(!candidates.size()) break; // no more candidates found.
2489 
2490  candidates.resize(min(candidates.size(), c.maxAggregateSize()-
2491  aggregate_->size()));
2492  aggregate_->add(candidates);
2493 
2494  }
2495 
2496  // try to merge aggregates consisting of only one nonisolated vertex with other aggregates
2497  if(aggregate_->size()==1 && c.maxAggregateSize()>1) {
2498  if(!graph.getVertexProperties(seed).isolated()) {
2499  Vertex mergedNeighbour = mergeNeighbour(seed, aggregates);
2500 
2501  if(mergedNeighbour != AggregatesMap<Vertex>::UNAGGREGATED) {
2502  // assign vertex to the neighbouring cluster
2503  aggregates[seed] = aggregates[mergedNeighbour];
2504  aggregate_->invalidate();
2505  }else{
2506  ++avg;
2507  minA=min(minA,static_cast<std::size_t>(1));
2508  maxA=max(maxA,static_cast<std::size_t>(1));
2509  ++oneAggregates;
2510  ++conAggregates;
2511  }
2512  }else{
2513  ++avg;
2514  minA=min(minA,static_cast<std::size_t>(1));
2515  maxA=max(maxA,static_cast<std::size_t>(1));
2516  ++oneAggregates;
2517  ++isoAggregates;
2518  }
2519  ++avg;
2520  }else{
2521  avg+=aggregate_->size();
2522  minA=min(minA,aggregate_->size());
2523  maxA=max(maxA,aggregate_->size());
2524  if(graph.getVertexProperties(seed).isolated())
2525  ++isoAggregates;
2526  else
2527  ++conAggregates;
2528  }
2529 
2530  }
2531 
2532  Dune::dinfo<<"connected aggregates: "<<conAggregates;
2533  Dune::dinfo<<" isolated aggregates: "<<isoAggregates;
2534  if(conAggregates+isoAggregates>0)
2535  Dune::dinfo<<" one node aggregates: "<<oneAggregates<<" min size="
2536  <<minA<<" max size="<<maxA
2537  <<" avg="<<avg/(conAggregates+isoAggregates)<<std::endl;
2538 
2539  delete aggregate_;
2540  return std::make_tuple(conAggregates+isoAggregates,isoAggregates,
2541  oneAggregates,skippedAggregates);
2542  }
2543 
2544 
2545  template<class G>
2546  Aggregator<G>::Stack::Stack(const MatrixGraph& graph, const Aggregator<G>& aggregatesBuilder,
2547  const AggregatesMap<Vertex>& aggregates)
2548  : graph_(graph), aggregatesBuilder_(aggregatesBuilder), aggregates_(aggregates), begin_(graph.begin()), end_(graph.end())
2549  {
2550  //vals_ = new Vertex[N];
2551  }
2552 
2553  template<class G>
2554  Aggregator<G>::Stack::~Stack()
2555  {
2556  //Dune::dverb << "Max stack size was "<<maxSize_<<" filled="<<filled_<<std::endl;
2557  //delete[] vals_;
2558  }
2559 
2560  template<class G>
2561  const typename Aggregator<G>::Vertex Aggregator<G>::Stack::NullEntry
2562  = std::numeric_limits<typename G::VertexDescriptor>::max();
2563 
2564  template<class G>
2565  inline typename G::VertexDescriptor Aggregator<G>::Stack::pop()
2566  {
2567  for(; begin_!=end_ && aggregates_[*begin_] != AggregatesMap<Vertex>::UNAGGREGATED; ++begin_) ;
2568 
2569  if(begin_!=end_)
2570  {
2571  typename G::VertexDescriptor current=*begin_;
2572  ++begin_;
2573  return current;
2574  }else
2575  return NullEntry;
2576  }
2577 
2578 #endif // DOXYGEN
2579 
2580  template<class V>
2581  void printAggregates2d(const AggregatesMap<V>& aggregates, int n, int m, std::ostream& os)
2582  {
2583  using std::max;
2584 
2585  std::ios_base::fmtflags oldOpts=os.flags();
2586 
2587  os.setf(std::ios_base::right, std::ios_base::adjustfield);
2588 
2589  V maxVal=0;
2590  int width=1;
2591 
2592  for(int i=0; i< n*m; i++)
2593  maxVal=max(maxVal, aggregates[i]);
2594 
2595  for(int i=10; i < 1000000; i*=10)
2596  if(maxVal/i>0)
2597  width++;
2598  else
2599  break;
2600 
2601  for(int j=0, entry=0; j < m; j++) {
2602  for(int i=0; i<n; i++, entry++) {
2603  os.width(width);
2604  os<<aggregates[entry]<<" ";
2605  }
2606 
2607  os<<std::endl;
2608  }
2609  os<<std::endl;
2610  os.flags(oldOpts);
2611  }
2612 
2613 
2614  } // namespace Amg
2615 
2616 } // namespace Dune
2617 
2618 
2619 #endif
Provides classes for building the matrix graph.
Parameter classes for customizing AMG.
Provides classes for handling internal properties in a graph.
Col col
Definition: matrixmatrix.hh:349
Matrix::ConstColIterator ColIter
Constant column iterator of the matrix.
Definition: aggregates.hh:271
std::vector< real_type >::iterator valIter_
Definition: aggregates.hh:187
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:584
iterator begin()
Definition: aggregates.hh:735
int id()
Get the id identifying the aggregate.
Norm norm_
The functor for calculating the norm.
Definition: aggregates.hh:300
void operator()([[maybe_unused]] const EdgeIterator &edge) const
Definition: aggregates.hh:599
MatrixGraph::VertexDescriptor Vertex
The vertex identifier.
Definition: aggregates.hh:918
AggregationCriterion()
Constructor.
Definition: aggregates.hh:64
const Matrix * matrix_
The matrix we work on.
Definition: aggregates.hh:355
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:404
void initRow(const Row &row, int index)
SymmetricMatrixDependency(const Parameters &parms)
Definition: aggregates.hh:166
M Matrix
The matrix type we build the dependency of.
Definition: aggregates.hh:256
G MatrixGraph
The matrix graph type used.
Definition: aggregates.hh:913
FieldTraits< typename M::field_type >::real_type operator()(const M &m) const
compute the norm of a matrix.
Definition: aggregates.hh:488
const AggregateDescriptor & operator[](const VertexDescriptor &v) const
Get the aggregate a vertex belongs to.
Norm norm_
The functor for calculating the norm.
Definition: aggregates.hh:361
SymmetricCriterion()
Definition: aggregates.hh:522
Dependency()
Definition: aggregates.hh:288
void examine(const ColIter &col)
M Matrix
The matrix type we build the dependency of.
Definition: aggregates.hh:317
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:261
iterator end()
Definition: aggregates.hh:740
UnSymmetricCriterion(const Parameters &parms)
Definition: aggregates.hh:539
void operator()(const typename MatrixGraph::ConstEdgeIterator &edge)
void initRow(const Row &row, int index)
Definition: aggregates.hh:199
static const Vertex NullEntry
Definition: aggregates.hh:1006
std::tuple< int, int, int, int > build(const M &m, G &graph, AggregatesMap< Vertex > &aggregates, const C &c, bool finestLevel)
Build the aggregates.
std::ostream & operator<<(std::ostream &os, const AggregationCriterion< T > &criterion)
Definition: aggregates.hh:111
void examine(const ColIter &col)
Definition: aggregates.hh:212
Dependency(const Parameters &parms)
Definition: aggregates.hh:284
int row_
index of the currently evaluated row.
Definition: aggregates.hh:183
friend class Stack
Definition: aggregates.hh:1034
void operator()(const typename MatrixGraph::ConstEdgeIterator &edge)
FrontNeighbourCounter(const MatrixGraph &front)
Constructor.
Matrix::row_type Row
Constant Row iterator of the matrix.
Definition: aggregates.hh:327
FieldTraits< typename M::field_type >::real_type operator()(const M &m, [[maybe_unused]] typename std::enable_if_t<!Dune::IsNumber< M >::value > *sfinae=nullptr) const
compute the norm of a matrix.
Definition: aggregates.hh:388
const Matrix * matrix_
The matrix we work on.
Definition: aggregates.hh:294
void examine(G &graph, const typename G::EdgeIterator &edge, const ColIter &col)
AggregateVisitor(const AggregatesMap< Vertex > &aggregates, const AggregateDescriptor &aggregate, Visitor &visitor)
Constructor.
Matrix::ConstColIterator ColIter
Constant column iterator of the matrix.
Definition: aggregates.hh:332
~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:1062
std::size_t * SphereMap
Type of the mapping of aggregate members onto distance spheres.
Definition: aggregates.hh:807
AggregationCriterion(const Parameters &parms)
Definition: aggregates.hh:68
Matrix::row_type Row
Constant Row iterator of the matrix.
Definition: aggregates.hh:266
VertexSet::size_type connectSize()
Get tne number of connections to other aggregates.
std::vector< real_type > vals_
Definition: aggregates.hh:186
FieldTraits< typename M::field_type >::real_type operator()(const M &m) const
compute the norm of a matrix.
Definition: aggregates.hh:504
N Norm
The norm to use for examining the matrix entries.
Definition: aggregates.hh:322
void printAggregates2d(const AggregatesMap< V > &aggregates, int n, int m, std::ostream &os)
Definition: aggregates.hh:2581
void invalidate()
Definition: aggregates.hh:820
const_iterator begin() const
Definition: aggregates.hh:723
real_type maxValue_
Definition: aggregates.hh:179
Norm norm_
The functor for calculating the norm.
Definition: aggregates.hh:181
void operator()(const typename MatrixGraph::ConstEdgeIterator &edge)
VertexSet::const_iterator const_iterator
Const iterator over a vertex list.
Definition: aggregates.hh:802
SymmetricCriterion(const Parameters &parms)
Definition: aggregates.hh:519
MatrixGraph::VertexDescriptor AggregateDescriptor
The type of the aggregate descriptor.
Definition: aggregates.hh:921
real_type maxValue_
Definition: aggregates.hh:359
void init(const Matrix *matrix)
Definition: aggregates.hh:193
real_type diagonal_
The norm of the current diagonal.
Definition: aggregates.hh:365
AggregateDescriptor * iterator
Definition: aggregates.hh:733
SymmetricDependency(const Parameters &parms)
Definition: aggregates.hh:346
FieldTraits< field_type >::real_type real_type
Definition: aggregates.hh:358
~Aggregator()
Destructor.
void examine(G &graph, const typename G::EdgeIterator &edge, const ColIter &col)
AggregateDescriptor & operator[](const VertexDescriptor &v)
Get the aggregate a vertex belongs to.
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 add(std::vector< Vertex > &vertex)
void examine(const ColIter &col)
Aggregator()
Constructor.
void operator()(const typename MatrixGraph::ConstEdgeIterator &edge)
Examine an edge.
FrontMarker(std::vector< Vertex > &front, MatrixGraph &graph)
Constructor.
void init(const Matrix *matrix)
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
const_iterator end() const
Definition: aggregates.hh:728
V AggregateDescriptor
The aggregate descriptor type.
Definition: aggregates.hh:578
static const V ISOLATED
Identifier of isolated vertices.
Definition: aggregates.hh:569
int row_
index of the currently evaluated row.
Definition: aggregates.hh:363
SymmetricMatrixDependency()
Definition: aggregates.hh:169
real_type diagonal_
The norm of the current diagonal.
Definition: aggregates.hh:304
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:357
AggregatesMap()
Constructs without allocating memory.
int value()
Access the current count.
std::tuple< int, int, int, int > buildAggregates(const M &matrix, G &graph, const C &criterion, bool finestLevel)
Build the aggregates.
SLList< VertexDescriptor, Allocator > VertexList
The type of a single linked list of vertex descriptors.
Definition: aggregates.hh:590
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.
void init(const Matrix *matrix)
UnSymmetricCriterion()
Definition: aggregates.hh:542
void operator()(const typename MatrixGraph::ConstEdgeIterator &edge)
FieldTraits< typename M::field_type >::real_type operator()(const M &m) const
compute the norm of a matrix.
Definition: aggregates.hh:471
const AggregateDescriptor * const_iterator
Definition: aggregates.hh:721
int row_
index of the currently evaluated row.
Definition: aggregates.hh:302
Stack(const MatrixGraph &graph, const Aggregator< G > &aggregatesBuilder, const AggregatesMap< Vertex > &aggregates)
FieldTraits< field_type >::real_type real_type
Definition: aggregates.hh:178
void initRow(const Row &row, int index)
M Matrix
The matrix type we build the dependency of.
Definition: aggregates.hh:137
FieldTraits< field_type >::real_type real_type
Definition: aggregates.hh:297
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:799
static const V UNAGGREGATED
Identifier of not yet aggregated vertices.
Definition: aggregates.hh:564
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:296
void operator()(const typename MatrixGraph::ConstEdgeIterator &edge)
bool isIsolated()
Definition: aggregates.hh:238
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
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:787
void seed(const Vertex &vertex)
Initialize the aggregate with one vertex.
SymmetricDependency()
Definition: aggregates.hh:349
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:573
real_type maxValue_
Definition: aggregates.hh:298
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:793
G MatrixGraph
Definition: aggregates.hh:783
void operator()(const typename MatrixGraph::ConstEdgeIterator &edge)
@ is_sign_preserving
Definition: aggregates.hh:481
@ is_sign_preserving
Definition: aggregates.hh:380
@ is_sign_preserving
Definition: aggregates.hh:497
@ is_sign_preserving
Definition: aggregates.hh:464
Definition: allocator.hh:9
PropertyMapTypeSelector< Amg::VertexVisitedTag, Amg::PropertiesGraph< G, Amg::VertexProperties, EP, VM, EM > >::Type get([[maybe_unused]] const Amg::VertexVisitedTag &tag, Amg::PropertiesGraph< G, Amg::VertexProperties, EP, VM, EM > &graph)
Definition: dependency.hh:291
derive error class from the base class in common
Definition: istlexception.hh:17
A generic dynamic dense matrix.
Definition: matrix.hh:559
typename Imp::BlockTraits< T >::field_type field_type
Export the type representing the underlying field.
Definition: matrix.hh:563
row_type::const_iterator ConstColIterator
Const iterator for the entries of each row.
Definition: matrix.hh:587
MatrixImp::DenseMatrixBase< T, A >::window_type row_type
The type implementing a matrix row.
Definition: matrix.hh:572
Base class of all aggregation criterions.
Definition: aggregates.hh:47
Dependency policy for symmetric matrices.
Definition: aggregates.hh:132
Dependency policy for symmetric matrices.
Definition: aggregates.hh:251
Dependency policy for symmetric matrices.
Definition: aggregates.hh:312
Norm that uses only the [N][N] entry of the block to determine couplings.
Definition: aggregates.hh:377
Norm that uses only the [0][0] entry of the block to determine couplings.
Definition: aggregates.hh:453
Functor using the row sum (infinity) norm to determine strong couplings.
Definition: aggregates.hh:461
Definition: aggregates.hh:478
Definition: aggregates.hh:494
Criterion taking advantage of symmetric matrices.
Definition: aggregates.hh:517
Criterion suitable for unsymmetric matrices.
Definition: aggregates.hh:537
Class for building the aggregates.
Definition: aggregates.hh:907
Class providing information about the mapping of the vertices onto aggregates.
Definition: aggregates.hh:558
A Dummy visitor that does nothing for each visited edge.
Definition: aggregates.hh:596
A class for temporarily storing the vertices of an aggregate in.
Definition: aggregates.hh:776
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
Iterator over all edges starting from a vertex.
Definition: graph.hh:93
The vertex iterator type of the graph.
Definition: graph.hh:207
All parameters for AMG.
Definition: parameters.hh:391