dune-istl  2.2.1
aggregates.hh
Go to the documentation of this file.
1 // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set ts=8 sw=2 et sts=2:
3 // $Id: aggregates.hh 1762 2013-01-19 11:21:36Z mblatt $
4 #ifndef DUNE_AMG_AGGREGATES_HH
5 #define DUNE_AMG_AGGREGATES_HH
6 
7 
8 #include"parameters.hh"
9 #include"graph.hh"
10 #include"properties.hh"
11 #include"combinedfunctor.hh"
12 
13 #include<dune/common/timer.hh>
14 #include<dune/common/tuples.hh>
15 #include<dune/common/stdstreams.hh>
16 #include<dune/common/poolallocator.hh>
17 #include<dune/common/sllist.hh>
18 
19 #include<utility>
20 #include<set>
21 #include<algorithm>
22 #include<limits>
23 #include<ostream>
24 
25 namespace Dune
26 {
27  namespace Amg
28  {
29 
43  template<class T>
44  class AggregationCriterion : public T
45  {
46 
47  public:
51  typedef T DependencyPolicy;
52 
63  : T()
64  {}
65 
67  : T(parms)
68  {}
78  void setDefaultValuesIsotropic(std::size_t dim, std::size_t diameter=2)
79  {
80  this->setMaxDistance(diameter-1);
81  std::size_t csize=1;
82 
83  for(;dim>0;dim--){
84  csize*=diameter;
85  this->setMaxDistance(this->maxDistance()+diameter-1);
86  }
87  this->setMinAggregateSize(csize);
88  this->setMaxAggregateSize(static_cast<std::size_t>(csize*1.5));
89  }
90 
101  void setDefaultValuesAnisotropic(std::size_t dim,std::size_t diameter=2)
102  {
103  setDefaultValuesIsotropic(dim, diameter);
104  this->setMaxDistance(this->maxDistance()+dim-1);
105  }
106  };
107 
108  template<class T>
109  std::ostream& operator<<(std::ostream& os, const AggregationCriterion<T>& criterion)
110  {
111  os<<"{ maxdistance="<<criterion.maxDistance()<<" minAggregateSize="
112  <<criterion.minAggregateSize()<< " maxAggregateSize="<<criterion.maxAggregateSize()
113  <<" connectivity="<<criterion.maxConnectivity()<<" debugLevel="<<criterion.debugLevel()<<"}";
114  return os;
115  }
116 
120  template<class M, class N>
121  class Dependency : public Parameters
122  {
123  public:
127  typedef M Matrix;
128 
132  typedef N Norm;
133 
137  typedef typename Matrix::row_type Row;
138 
143 
144  void init(const Matrix* matrix);
145 
146  void initRow(const Row& row, int index);
147 
148  void examine(const ColIter& col);
149 
150  template<class G>
151  void examine(G& graph, const typename G::EdgeIterator& edge, const ColIter& col);
152 
153  bool isIsolated();
154 
155  Dependency(const Parameters& parms)
156  : Parameters(parms)
157  {}
158 
160  :Parameters()
161  {}
162 
163  protected:
165  const Matrix* matrix_;
171  int row_;
174  };
175 
179  template<class M, class N>
181  {
182  public:
186  typedef M Matrix;
187 
191  typedef N Norm;
192 
196  typedef typename Matrix::row_type Row;
197 
202 
203  void init(const Matrix* matrix);
204 
205  void initRow(const Row& row, int index);
206 
207  void examine(const ColIter& col);
208 
209  template<class G>
210  void examine(G& graph, const typename G::EdgeIterator& edge, const ColIter& col);
211 
212  bool isIsolated();
213 
214 
216  : Parameters(parms)
217  {}
219  :Parameters()
220  {}
221 
222  protected:
224  const Matrix* matrix_;
230  int row_;
233  };
234 
239  template<int N>
240  class Diagonal
241  {
242  public:
243  enum{ /* @brief We preserve the sign.*/
245  };
246 
251  template<class M>
252  typename M::field_type operator()(const M& m) const
253  {
254  return m[N][N];
255  }
256  };
257 
262  class FirstDiagonal : public Diagonal<0>
263  {};
264 
270  struct RowSum
271  {
272 
273  enum{ /* @brief We preserve the sign.*/
275  };
280  template<class M>
281  typename M::field_type operator()(const M& m) const
282  {
283  return m.infinity_norm();
284  }
285  };
286 
288  {
289 
290  enum{ /* @brief We preserve the sign.*/
292  };
297  template<class M>
298  typename M::field_type operator()(const M& m) const
299  {
300  return m.frobenius_norm();
301  }
302  };
304  {
305 
306  enum{ /* @brief We preserve the sign.*/
308  };
313  template<class M>
314  typename M::field_type operator()(const M& m) const
315  {
316  return 1;
317  }
318  };
328  template<class M, class Norm>
329  class SymmetricCriterion : public AggregationCriterion<SymmetricDependency<M,Norm> >
330  {
331  public:
334  {}
336  {}
337  };
338 
339 
351  template<class M, class Norm>
352  class UnSymmetricCriterion : public AggregationCriterion<Dependency<M,Norm> >
353  {
354  public:
356  : AggregationCriterion<Dependency<M,Norm> >(parms)
357  {}
359  {}
360  };
361  // forward declaration
362  template<class G> class Aggregator;
363 
364 
372  template<class V>
374  {
375  public:
376 
380  static const V UNAGGREGATED;
381 
385  static const V ISOLATED;
389  typedef V VertexDescriptor;
390 
395 
400  typedef PoolAllocator<VertexDescriptor,100> Allocator;
401 
406  typedef SLList<VertexDescriptor,Allocator> VertexList;
407 
412  {
413  public:
414  template<class EdgeIterator>
415  void operator()(const EdgeIterator& egde) const
416  {}
417  };
418 
419 
423  AggregatesMap();
424 
430  AggregatesMap(std::size_t noVertices);
431 
435  ~AggregatesMap();
436 
448  template<class M, class G, class C>
449  tuple<int,int,int,int> buildAggregates(const M& matrix, G& graph, const C& criterion,
450  bool finestLevel);
451 
471  template<bool reset, class G, class F, class VM>
472  std::size_t breadthFirstSearch(const VertexDescriptor& start,
473  const AggregateDescriptor& aggregate,
474  const G& graph,
475  F& aggregateVisitor,
476  VM& visitedMap) const;
477 
501  template<bool remove, bool reset, class G, class L, class F1, class F2, class VM>
502  std::size_t breadthFirstSearch(const VertexDescriptor& start,
503  const AggregateDescriptor& aggregate,
504  const G& graph, L& visited, F1& aggregateVisitor,
505  F2& nonAggregateVisitor,
506  VM& visitedMap) const;
507 
513  void allocate(std::size_t noVertices);
514 
518  std::size_t noVertices() const;
519 
523  void free();
524 
531 
537  const AggregateDescriptor& operator[](const VertexDescriptor& v) const;
538 
540 
542  {
543  return aggregates_;
544  }
545 
547  {
548  return aggregates_+noVertices();
549  }
550 
552 
554  {
555  return aggregates_;
556  }
557 
559  {
560  return aggregates_+noVertices();
561  }
562  private:
564  AggregatesMap(const AggregatesMap<V>& map)
565  {
566  throw "Auch!";
567  }
568 
570  AggregatesMap<V>& operator=(const AggregatesMap<V>& map)
571  {
572  throw "Auch!";
573  return this;
574  }
575 
579  AggregateDescriptor* aggregates_;
580 
584  std::size_t noVertices_;
585  };
586 
590  template<class G, class C>
591  void buildDependency(G& graph,
592  const typename C::Matrix& matrix,
593  C criterion,
594  bool finestLevel);
595 
600  template<class G, class S>
601  class Aggregate
602  {
603 
604  public:
605 
606  /***
607  * @brief The type of the matrix graph we work with.
608  */
609  typedef G MatrixGraph;
614 
619  typedef PoolAllocator<Vertex,100> Allocator;
620 
625  typedef SLList<Vertex,Allocator> VertexList;
626 
627 
632  typedef S VertexSet;
633 
635  typedef typename VertexList::const_iterator const_iterator;
636 
640  typedef std::size_t* SphereMap;
641 
649  Aggregate(const MatrixGraph& graph, AggregatesMap<Vertex>& aggregates,
650  VertexSet& connectivity);
651 
658  void reconstruct(const Vertex& vertex);
659 
663  void seed(const Vertex& vertex);
664 
668  void add(const Vertex& vertex);
669 
673  void clear();
674 
678  typename VertexList::size_type size();
682  typename VertexList::size_type connectSize();
683 
687  int id();
688 
690  const_iterator begin() const;
691 
693  const_iterator end() const;
694 
695  private:
699  VertexList vertices_;
700 
705  int id_;
706 
710  const MatrixGraph& graph_;
711 
715  AggregatesMap<Vertex>& aggregates_;
716 
720  VertexSet& connected_;
721  };
722 
726  template<class G>
727  class Aggregator
728  {
729  public:
730 
734  typedef G MatrixGraph;
735 
740 
743 
747  Aggregator();
748 
752  ~Aggregator();
753 
770  template<class M, class C>
771  tuple<int,int,int,int> build(const M& m, G& graph,
772  AggregatesMap<Vertex>& aggregates, const C& c,
773  bool finestLevel);
774  private:
779  typedef PoolAllocator<Vertex,100> Allocator;
780 
784  typedef SLList<Vertex,Allocator> VertexList;
785 
789  typedef std::set<Vertex,std::less<Vertex>,Allocator> VertexSet;
790 
794  typedef std::size_t* SphereMap;
795 
799  MatrixGraph* graph_;
800 
805 
809  VertexList front_;
810 
814  VertexSet connected_;
815 
819  int size_;
820 
824  class Stack
825  {
826  public:
827  static const Vertex NullEntry;
828 
829  Stack(const MatrixGraph& graph,
830  const Aggregator<G>& aggregatesBuilder,
831  const AggregatesMap<Vertex>& aggregates);
832  ~Stack();
833  bool push(const Vertex& v);
834  void fill();
835  Vertex pop();
836  private:
837  enum{ N = 256000 };
838 
840  const MatrixGraph& graph_;
842  const Aggregator<G>& aggregatesBuilder_;
844  const AggregatesMap<Vertex>& aggregates_;
846  int size_;
847  int maxSize_;
849  int head_;
850  int filled_;
851 
853  Vertex* vals_;
854 
855  void localPush(const Vertex& v);
856  };
857 
858  friend class Stack;
859 
870  template<class V>
871  void visitAggregateNeighbours(const Vertex& vertex, const AggregateDescriptor& aggregate,
872  const AggregatesMap<Vertex>& aggregates,
873  V& visitor) const;
874 
879  template<class V>
880  class AggregateVisitor
881  {
882  public:
886  typedef V Visitor;
894  AggregateVisitor(const AggregatesMap<Vertex>& aggregates, const AggregateDescriptor& aggregate,
895  Visitor& visitor);
896 
903  void operator()(const typename MatrixGraph::ConstEdgeIterator& edge);
904 
905  private:
907  const AggregatesMap<Vertex>& aggregates_;
909  AggregateDescriptor aggregate_;
911  Visitor* visitor_;
912  };
913 
917  class Counter
918  {
919  public:
921  Counter();
923  int value();
924 
925  protected:
927  void increment();
929  void decrement();
930 
931  private:
932  int count_;
933  };
934 
935 
940  class FrontNeighbourCounter : public Counter
941  {
942  public:
947  FrontNeighbourCounter(const MatrixGraph& front);
948 
949  void operator()(const typename MatrixGraph::ConstEdgeIterator& edge);
950 
951  private:
952  const MatrixGraph& graph_;
953  };
954 
959  int noFrontNeighbours(const Vertex& vertex) const;
960 
964  class TwoWayCounter : public Counter
965  {
966  public:
967  void operator()(const typename MatrixGraph::ConstEdgeIterator& edge);
968  };
969 
981  int twoWayConnections(const Vertex&, const AggregateDescriptor& aggregate,
982  const AggregatesMap<Vertex>& aggregates) const;
983 
987  class OneWayCounter : public Counter
988  {
989  public:
990  void operator()(const typename MatrixGraph::ConstEdgeIterator& edge);
991  };
992 
1004  int oneWayConnections(const Vertex&, const AggregateDescriptor& aggregate,
1005  const AggregatesMap<Vertex>& aggregates) const;
1006 
1013  class ConnectivityCounter : public Counter
1014  {
1015  public:
1022  ConnectivityCounter(const VertexSet& connected, const AggregatesMap<Vertex>& aggregates);
1023 
1024  void operator()(const typename MatrixGraph::ConstEdgeIterator& edge);
1025 
1026  private:
1028  const VertexSet& connected_;
1030  const AggregatesMap<Vertex>& aggregates_;
1031 
1032  };
1033 
1045  double connectivity(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates) const;
1053  bool connected(const Vertex& vertex, const AggregateDescriptor& aggregate,
1054  const AggregatesMap<Vertex>& aggregates) const;
1055 
1063  bool connected(const Vertex& vertex, const SLList<AggregateDescriptor>& aggregateList,
1064  const AggregatesMap<Vertex>& aggregates) const;
1065 
1073  class DependencyCounter: public Counter
1074  {
1075  public:
1079  DependencyCounter();
1080 
1081  void operator()(const typename MatrixGraph::ConstEdgeIterator& edge);
1082  };
1083 
1090  class FrontMarker
1091  {
1092  public:
1099  FrontMarker(VertexList& front, MatrixGraph& graph);
1100 
1101  void operator()(const typename MatrixGraph::ConstEdgeIterator& edge);
1102 
1103  private:
1105  VertexList& front_;
1107  MatrixGraph& graph_;
1108  };
1109 
1116  void markFront(const AggregatesMap<Vertex>& aggregates);
1117 
1121  void unmarkFront();
1122 
1137  int unusedNeighbours(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates) const;
1138 
1152  std::pair<int,int> neighbours(const Vertex& vertex,
1153  const AggregateDescriptor& aggregate,
1154  const AggregatesMap<Vertex>& aggregates) const;
1171  int aggregateNeighbours(const Vertex& vertex, const AggregateDescriptor& aggregate, const AggregatesMap<Vertex>& aggregates) const;
1172 
1180  bool admissible(const Vertex& vertex, const AggregateDescriptor& aggregate, const AggregatesMap<Vertex>& aggregates) const;
1181 
1188  void seedFromFront(Stack& stack, bool isolated);
1189 
1197  std::size_t distance(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates);
1198 
1207  Vertex mergeNeighbour(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates) const;
1208 
1217  void nonisoNeighbourAggregate(const Vertex& vertex,
1218  const AggregatesMap<Vertex>& aggregates,
1219  SLList<Vertex>& neighbours) const;
1220 
1228  template<class C>
1229  void growAggregate(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates, const C& c);
1230  template<class C>
1231  void growIsolatedAggregate(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates, const C& c);
1232  };
1233 
1234 #ifndef DOXYGEN
1235 
1236  template<class M, class N>
1237  inline void SymmetricDependency<M,N>::init(const Matrix* matrix)
1238  {
1239  matrix_ = matrix;
1240  }
1241 
1242  template<class M, class N>
1243  inline void SymmetricDependency<M,N>::initRow(const Row& row, int index)
1244  {
1245  maxValue_ = std::min(- std::numeric_limits<typename Matrix::field_type>::max(), std::numeric_limits<typename Matrix::field_type>::min());
1246  row_ = index;
1247  diagonal_ = norm_(matrix_->operator[](row_)[row_]);
1248  }
1249 
1250  template<class M, class N>
1251  inline void SymmetricDependency<M,N>::examine(const ColIter& col)
1252  {
1253  typename Matrix::field_type eij = norm_(*col);
1254  typename Matrix::field_type eji = norm_(matrix_->operator[](col.index())[row_]);
1255 
1256  // skip positive offdiagonals if norm preserves sign of them.
1257  if(!N::is_sign_preserving || eij<0 || eji<0)
1258  maxValue_ = std::max(maxValue_,
1259  eij /diagonal_ * eji/
1260  norm_(matrix_->operator[](col.index())[col.index()]));
1261  }
1262 
1263  template<class M, class N>
1264  template<class G>
1265  inline void SymmetricDependency<M,N>::examine(G& graph, const typename G::EdgeIterator& edge, const ColIter& col)
1266  {
1267  typename Matrix::field_type eij = norm_(*col);
1268  typename Matrix::field_type eji = norm_(matrix_->operator[](col.index())[row_]);
1269 
1270  // skip positve offdiagonals if norm preserves sign of them.
1271  if(!N::is_sign_preserving || (eij<0 || eji<0))
1272  if(eji / norm_(matrix_->operator[](edge.target())[edge.target()]) *
1273  eij/ diagonal_ > alpha() * maxValue_){
1274  edge.properties().setDepends();
1275  edge.properties().setInfluences();
1276 
1277  typename G::EdgeProperties& other = graph.getEdgeProperties(edge.target(), edge.source());
1278  other.setInfluences();
1279  other.setDepends();
1280  }
1281  }
1282 
1283  template<class M, class N>
1285  {
1286  return maxValue_ < beta();
1287  }
1288 
1289 
1290  template<class M, class N>
1291  inline void Dependency<M,N>::init(const Matrix* matrix)
1292  {
1293  matrix_ = matrix;
1294  }
1295 
1296  template<class M, class N>
1297  inline void Dependency<M,N>::initRow(const Row& row, int index)
1298  {
1299  maxValue_ = std::min(- std::numeric_limits<typename Matrix::field_type>::max(), std::numeric_limits<typename Matrix::field_type>::min());
1300  row_ = index;
1301  diagonal_ = norm_(matrix_->operator[](row_)[row_]);
1302  }
1303 
1304  template<class M, class N>
1305  inline void Dependency<M,N>::examine(const ColIter& col)
1306  {
1307  maxValue_ = std::max(maxValue_,
1308  -norm_(*col));
1309  }
1310 
1311  template<class M, class N>
1312  template<class G>
1313  inline void Dependency<M,N>::examine(G& graph, const typename G::EdgeIterator& edge, const ColIter& col)
1314  {
1315  if(-norm_(*col) >= maxValue_ * alpha()){
1316  edge.properties().setDepends();
1317  typedef typename G::EdgeDescriptor ED;
1318  ED e= graph.findEdge(edge.target(), edge.source());
1319  if(e!=std::numeric_limits<ED>::max())
1320  {
1321  typename G::EdgeProperties& other = graph.getEdgeProperties(e);
1322  other.setInfluences();
1323  }
1324  }
1325  }
1326 
1327  template<class M, class N>
1328  inline bool Dependency<M,N>::isIsolated()
1329  {
1330  return maxValue_ < beta() * diagonal_;
1331  }
1332 
1333  template<class G,class S>
1334  Aggregate<G,S>::Aggregate(const MatrixGraph& graph, AggregatesMap<Vertex>& aggregates,
1335  VertexSet& connected)
1336  : vertices_(), id_(-1), graph_(graph), aggregates_(aggregates),
1337  connected_(connected)
1338  {}
1339 
1340  template<class G,class S>
1341  void Aggregate<G,S>::reconstruct(const Vertex& vertex)
1342  {
1343  vertices_.push_back(vertex);
1344  typedef typename VertexList::const_iterator iterator;
1345  iterator begin = vertices_.begin();
1346  iterator end = vertices_.end();
1347  throw "Not yet implemented";
1348 
1349  while(begin!=end){
1350  //for();
1351  }
1352 
1353  }
1354 
1355  template<class G,class S>
1356  inline void Aggregate<G,S>::seed(const Vertex& vertex)
1357  {
1358  dvverb<<"Connected cleared"<<std::endl;
1359  connected_.clear();
1360  vertices_.clear();
1361  connected_.insert(vertex);
1362  dvverb << " Inserting "<<vertex<<" size="<<connected_.size();
1363  id_ = vertex;
1364  add(vertex);
1365  }
1366 
1367 
1368  template<class G,class S>
1369  inline void Aggregate<G,S>::add(const Vertex& vertex)
1370  {
1371  vertices_.push_back(vertex);
1372  aggregates_[vertex]=id_;
1373 
1374  typedef typename MatrixGraph::ConstEdgeIterator iterator;
1375  const iterator end = graph_.endEdges(vertex);
1376  for(iterator edge = graph_.beginEdges(vertex); edge != end; ++edge){
1377  dvverb << " Inserting "<<aggregates_[edge.target()];
1378  connected_.insert(aggregates_[edge.target()]);
1379  dvverb <<" size="<<connected_.size();
1380  }
1381  dvverb <<std::endl;
1382  }
1383  template<class G,class S>
1384  inline void Aggregate<G,S>::clear()
1385  {
1386  vertices_.clear();
1387  connected_.clear();
1388  id_=-1;
1389  }
1390 
1391  template<class G,class S>
1392  inline typename Aggregate<G,S>::VertexList::size_type
1394  {
1395  return vertices_.size();
1396  }
1397 
1398  template<class G,class S>
1399  inline typename Aggregate<G,S>::VertexList::size_type
1401  {
1402  return connected_.size();
1403  }
1404 
1405  template<class G,class S>
1406  inline int Aggregate<G,S>::id()
1407  {
1408  return id_;
1409  }
1410 
1411  template<class G,class S>
1412  inline typename Aggregate<G,S>::const_iterator Aggregate<G,S>::begin() const
1413  {
1414  return vertices_.begin();
1415  }
1416 
1417  template<class G,class S>
1418  inline typename Aggregate<G,S>::const_iterator Aggregate<G,S>::end() const
1419  {
1420  return vertices_.end();
1421  }
1422 
1423  template<class V>
1424  const V AggregatesMap<V>::UNAGGREGATED = std::numeric_limits<V>::max();
1425 
1426  template<class V>
1427  const V AggregatesMap<V>::ISOLATED = std::numeric_limits<V>::max()-1;
1428 
1429  template<class V>
1431  : aggregates_(0)
1432  {}
1433 
1434  template<class V>
1436  {
1437  if(aggregates_!=0)
1438  delete[] aggregates_;
1439  }
1440 
1441 
1442  template<class V>
1443  inline AggregatesMap<V>::AggregatesMap(std::size_t noVertices)
1444  {
1445  allocate(noVertices);
1446  }
1447 
1448  template<class V>
1449  inline std::size_t AggregatesMap<V>::noVertices() const
1450  {
1451  return noVertices_;
1452  }
1453 
1454  template<class V>
1455  inline void AggregatesMap<V>::allocate(std::size_t noVertices)
1456  {
1457  aggregates_ = new AggregateDescriptor[noVertices];
1458  noVertices_ = noVertices;
1459 
1460  for(std::size_t i=0; i < noVertices; i++)
1461  aggregates_[i]=UNAGGREGATED;
1462  }
1463 
1464  template<class V>
1465  inline void AggregatesMap<V>::free()
1466  {
1467  assert(aggregates_ != 0);
1468  delete[] aggregates_;
1469  aggregates_=0;
1470  }
1471 
1472  template<class V>
1473  inline typename AggregatesMap<V>::AggregateDescriptor&
1474  AggregatesMap<V>::operator[](const VertexDescriptor& v)
1475  {
1476  return aggregates_[v];
1477  }
1478 
1479  template<class V>
1480  inline const typename AggregatesMap<V>::AggregateDescriptor&
1481  AggregatesMap<V>::operator[](const VertexDescriptor& v) const
1482  {
1483  return aggregates_[v];
1484  }
1485 
1486  template<class V>
1487  template<bool reset, class G, class F,class VM>
1488  inline std::size_t AggregatesMap<V>::breadthFirstSearch(const V& start,
1489  const AggregateDescriptor& aggregate,
1490  const G& graph, F& aggregateVisitor,
1491  VM& visitedMap) const
1492  {
1493  VertexList vlist;
1494 
1495  DummyEdgeVisitor dummy;
1496  return breadthFirstSearch<true,reset>(start, aggregate, graph, vlist, aggregateVisitor, dummy, visitedMap);
1497  }
1498 
1499  template<class V>
1500  template<bool remove, bool reset, class G, class L, class F1, class F2, class VM>
1501  std::size_t AggregatesMap<V>::breadthFirstSearch(const V& start,
1502  const AggregateDescriptor& aggregate,
1503  const G& graph,
1504  L& visited,
1505  F1& aggregateVisitor,
1506  F2& nonAggregateVisitor,
1507  VM& visitedMap) const
1508  {
1509  typedef typename L::const_iterator ListIterator;
1510  int visitedSpheres = 0;
1511 
1512  visited.push_back(start);
1513  put(visitedMap, start, true);
1514 
1515  ListIterator current = visited.begin();
1516  ListIterator end = visited.end();
1517  std::size_t i=0, size=visited.size();
1518 
1519  // visit the neighbours of all vertices of the
1520  // current sphere.
1521  while(current != end){
1522 
1523  for(;i<size; ++current, ++i){
1524  typedef typename G::ConstEdgeIterator EdgeIterator;
1525  const EdgeIterator endEdge = graph.endEdges(*current);
1526 
1527  for(EdgeIterator edge = graph.beginEdges(*current);
1528  edge != endEdge; ++edge){
1529 
1530  if(aggregates_[edge.target()]==aggregate){
1531  if(!get(visitedMap, edge.target())){
1532  put(visitedMap, edge.target(), true);
1533  visited.push_back(edge.target());
1534  aggregateVisitor(edge);
1535  }
1536  }else
1537  nonAggregateVisitor(edge);
1538  }
1539  }
1540  end = visited.end();
1541  size = visited.size();
1542  if(current != end)
1543  visitedSpheres++;
1544  }
1545 
1546  if(reset)
1547  for(current = visited.begin(); current != end; ++current)
1548  put(visitedMap, *current, false);
1549 
1550 
1551  if(remove)
1552  visited.clear();
1553 
1554  return visitedSpheres;
1555  }
1556 
1557  template<class G>
1559  : graph_(0), aggregate_(0), front_(), connected_(), size_(-1)
1560  {}
1561 
1562  template<class G>
1564  {
1565  size_=-1;
1566  }
1567 
1568  template<class G, class C>
1569  void buildDependency(G& graph,
1570  const typename C::Matrix& matrix,
1571  C criterion, bool firstlevel)
1572  {
1573  // The Criterion we use for building the dependency.
1574  typedef C Criterion;
1575 
1576  // assert(graph.isBuilt());
1577  typedef typename C::Matrix Matrix;
1578  typedef typename G::VertexIterator VertexIterator;
1579 
1580  criterion.init(&matrix);
1581 
1582  for(VertexIterator vertex = graph.begin(); vertex != graph.end(); ++vertex){
1583  typedef typename Matrix::row_type Row;
1584 
1585  const Row& row = matrix[*vertex];
1586 
1587  // Tell the criterion what row we will examine now
1588  // This might for example be used for calculating the
1589  // maximum offdiagonal value
1590  criterion.initRow(row, *vertex);
1591 
1592  // On a first path all columns are examined. After this
1593  // the calculator should know whether the vertex is isolated.
1594  typedef typename Matrix::ConstColIterator ColIterator;
1595  ColIterator end = row.end();
1596  typename Matrix::field_type absoffdiag=0;
1597 
1598  if(firstlevel){
1599  for(ColIterator col = row.begin(); col != end; ++col)
1600  if(col.index()!=*vertex){
1601  criterion.examine(col);
1602  absoffdiag=std::max(absoffdiag, col->frobenius_norm());
1603  }
1604 
1605  if(absoffdiag==0)
1606  vertex.properties().setExcludedBorder();
1607  }
1608  else
1609  for(ColIterator col = row.begin(); col != end; ++col)
1610  if(col.index()!=*vertex)
1611  criterion.examine(col);
1612 
1613  // reset the vertex properties
1614  //vertex.properties().reset();
1615 
1616  // Check whether the vertex is isolated.
1617  if(criterion.isIsolated()){
1618  //std::cout<<"ISOLATED: "<<*vertex<<std::endl;
1619  vertex.properties().setIsolated();
1620  }else{
1621  // Examine all the edges beginning at this vertex.
1622  typedef typename G::EdgeIterator EdgeIterator;
1623  typedef typename Matrix::ConstColIterator ColIterator;
1624  EdgeIterator end = vertex.end();
1625  ColIterator col = matrix[*vertex].begin();
1626 
1627  for(EdgeIterator edge = vertex.begin(); edge!= end; ++edge, ++col){
1628  // Move to the right column.
1629  while(col.index()!=edge.target())
1630  ++col;
1631  criterion.examine(graph, edge, col);
1632  }
1633  }
1634 
1635  }
1636  }
1637 
1638 
1639  template<class G>
1640  template<class V>
1641  inline Aggregator<G>::AggregateVisitor<V>::AggregateVisitor(const AggregatesMap<Vertex>& aggregates,
1642  const AggregateDescriptor& aggregate, V& visitor)
1643  : aggregates_(aggregates), aggregate_(aggregate), visitor_(&visitor)
1644  {}
1645 
1646  template<class G>
1647  template<class V>
1648  inline void Aggregator<G>::AggregateVisitor<V>::operator()(const typename MatrixGraph::ConstEdgeIterator& edge)
1649  {
1650  if(aggregates_[edge.target()]==aggregate_)
1651  visitor_->operator()(edge);
1652  }
1653 
1654  template<class G>
1655  template<class V>
1656  inline void Aggregator<G>::visitAggregateNeighbours(const Vertex& vertex,
1657  const AggregateDescriptor& aggregate,
1658  const AggregatesMap<Vertex>& aggregates,
1659  V& visitor) const
1660  {
1661  // Only evaluates for edge pointing to the aggregate
1662  AggregateVisitor<V> v(aggregates, aggregate, visitor);
1663  visitNeighbours(*graph_, vertex, v);
1664  }
1665 
1666 
1667  template<class G>
1668  inline Aggregator<G>::Counter::Counter()
1669  : count_(0)
1670  {}
1671 
1672  template<class G>
1673  inline void Aggregator<G>::Counter::increment()
1674  {
1675  ++count_;
1676  }
1677 
1678  template<class G>
1679  inline void Aggregator<G>::Counter::decrement()
1680  {
1681  --count_;
1682  }
1683  template<class G>
1684  inline int Aggregator<G>::Counter::value()
1685  {
1686  return count_;
1687  }
1688 
1689  template<class G>
1690  inline void Aggregator<G>::TwoWayCounter::operator()(const typename MatrixGraph::ConstEdgeIterator& edge)
1691  {
1692  if(edge.properties().isTwoWay())
1693  Counter::increment();
1694  }
1695 
1696  template<class G>
1697  int Aggregator<G>::twoWayConnections(const Vertex& vertex, const AggregateDescriptor& aggregate,
1698  const AggregatesMap<Vertex>& aggregates) const
1699  {
1700  TwoWayCounter counter;
1701  visitAggregateNeighbours(vertex, aggregate, aggregates, counter);
1702  return counter.value();
1703  }
1704 
1705  template<class G>
1706  int Aggregator<G>::oneWayConnections(const Vertex& vertex, const AggregateDescriptor& aggregate,
1707  const AggregatesMap<Vertex>& aggregates) const
1708  {
1709  OneWayCounter counter;
1710  visitAggregateNeighbours(vertex, aggregate, aggregates, counter);
1711  return counter.value();
1712  }
1713 
1714  template<class G>
1715  inline void Aggregator<G>::OneWayCounter::operator()(const typename MatrixGraph::ConstEdgeIterator& edge)
1716  {
1717  if(edge.properties().isOneWay())
1718  Counter::increment();
1719  }
1720 
1721  template<class G>
1722  inline Aggregator<G>::ConnectivityCounter::ConnectivityCounter(const VertexSet& connected,
1723  const AggregatesMap<Vertex>& aggregates)
1724  : Counter(), connected_(connected), aggregates_(aggregates)
1725  {}
1726 
1727 
1728  template<class G>
1729  inline void Aggregator<G>::ConnectivityCounter::operator()(const typename MatrixGraph::ConstEdgeIterator& edge)
1730  {
1731  if(connected_.find(aggregates_[edge.target()]) == connected_.end() || aggregates_[edge.target()]==AggregatesMap<Vertex>::UNAGGREGATED)
1732  // Would be a new connection
1733  Counter::increment();
1734  else{
1735  Counter::increment();
1736  Counter::increment();
1737  }
1738  }
1739 
1740  template<class G>
1741  inline double Aggregator<G>::connectivity(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates) const
1742  {
1743  ConnectivityCounter counter(connected_, aggregates);
1744  double noNeighbours=visitNeighbours(*graph_, vertex, counter);
1745  return (double)counter.value()/noNeighbours;
1746  }
1747 
1748  template<class G>
1749  inline Aggregator<G>::DependencyCounter::DependencyCounter()
1750  : Counter()
1751  {}
1752 
1753  template<class G>
1754  inline void Aggregator<G>::DependencyCounter::operator()(const typename MatrixGraph::ConstEdgeIterator& edge)
1755  {
1756  if(edge.properties().depends())
1757  Counter::increment();
1758  if(edge.properties().influences())
1759  Counter::increment();
1760  }
1761 
1762  template<class G>
1763  int Aggregator<G>::unusedNeighbours(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates) const
1764  {
1765  return aggregateNeighbours(vertex, AggregatesMap<Vertex>::UNAGGREGATED, aggregates);
1766  }
1767 
1768  template<class G>
1769  std::pair<int,int> Aggregator<G>::neighbours(const Vertex& vertex,
1770  const AggregateDescriptor& aggregate,
1771  const AggregatesMap<Vertex>& aggregates) const
1772  {
1773  DependencyCounter unused, aggregated;
1774  typedef AggregateVisitor<DependencyCounter> Counter;
1775  typedef tuple<Counter,Counter> CounterTuple;
1776  CombinedFunctor<CounterTuple> visitors(CounterTuple(Counter(aggregates, AggregatesMap<Vertex>::UNAGGREGATED, unused), Counter(aggregates, aggregate, aggregated)));
1777  visitNeighbours(*graph_, vertex, visitors);
1778  return std::make_pair(unused.value(), aggregated.value());
1779 }
1780 
1781 
1782  template<class G>
1783  int Aggregator<G>::aggregateNeighbours(const Vertex& vertex, const AggregateDescriptor& aggregate, const AggregatesMap<Vertex>& aggregates) const
1784  {
1785  DependencyCounter counter;
1786  visitAggregateNeighbours(vertex, aggregate, aggregates, counter);
1787  return counter.value();
1788  }
1789 
1790  template<class G>
1791  std::size_t Aggregator<G>::distance(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates)
1792  {
1793  typename PropertyMapTypeSelector<VertexVisitedTag,G>::Type visitedMap = get(VertexVisitedTag(), *graph_);
1794  VertexList vlist;
1795  typename AggregatesMap<Vertex>::DummyEdgeVisitor dummy;
1796  return aggregates.template breadthFirstSearch<true,true>(vertex,
1797  aggregate_->id(), *graph_,
1798  vlist, dummy, dummy, visitedMap);
1799  }
1800 
1801  template<class G>
1802  inline Aggregator<G>::FrontMarker::FrontMarker(VertexList& front, MatrixGraph& graph)
1803  : front_(front), graph_(graph)
1804  {}
1805 
1806  template<class G>
1807  inline void Aggregator<G>::FrontMarker::operator()(const typename MatrixGraph::ConstEdgeIterator& edge)
1808  {
1809  Vertex target = edge.target();
1810 
1811  if(!graph_.getVertexProperties(target).front()){
1812  front_.push_back(target);
1813  graph_.getVertexProperties(target).setFront();
1814  }
1815  }
1816 
1817 
1818  template<class G>
1819  void Aggregator<G>::markFront(const AggregatesMap<Vertex>& aggregates)
1820  {
1821  assert(front_.size()==0);
1822  FrontMarker frontBuilder(front_, *graph_);
1823  typedef typename Aggregate<G,VertexSet>::const_iterator Iterator;
1824 
1825  for(Iterator vertex=aggregate_->begin(); vertex != aggregate_->end(); ++vertex)
1826  visitAggregateNeighbours(*vertex, AggregatesMap<Vertex>::UNAGGREGATED, aggregates, frontBuilder);
1827 
1828  }
1829 
1830  template<class G>
1831  inline bool Aggregator<G>::admissible(const Vertex& vertex, const AggregateDescriptor& aggregate, const AggregatesMap<Vertex>& aggregates) const
1832  {
1833  // Todo
1834  Dune::dvverb<<" Admissible not yet implemented!"<<std::endl;
1835  return true;
1836  //Situation 1: front node depends on two nodes. Then these
1837  // have to be strongly connected to each other
1838 
1839  // Iterate over all all neighbours of front node
1840  typedef typename MatrixGraph::ConstEdgeIterator Iterator;
1841  Iterator vend = graph_->endEdges(vertex);
1842  for(Iterator edge = graph_->beginEdges(vertex); edge != vend; ++edge){
1843  // if(edge.properties().depends() && !edge.properties().influences()
1844  if(edge.properties().isStrong()
1845  && aggregates[edge.target()]==aggregate)
1846  {
1847  // Search for another link to the aggregate
1848  Iterator edge1 = edge;
1849  for(++edge1; edge1 != vend; ++edge1){
1850  //if(edge1.properties().depends() && !edge1.properties().influences()
1851  if(edge1.properties().isStrong()
1852  && aggregates[edge.target()]==aggregate)
1853  {
1854  //Search for an edge connecting the two vertices that is
1855  //strong
1856  bool found=false;
1857  Iterator v2end = graph_->endEdges(edge.target());
1858  for(Iterator edge2 = graph_->beginEdges(edge.target()); edge2 != v2end; ++edge2){
1859  if(edge2.target()==edge1.target() &&
1860  edge2.properties().isStrong()){
1861  found =true;
1862  break;
1863  }
1864  }
1865  if(found)
1866  return true;
1867  }
1868  }
1869  }
1870  }
1871 
1872  // Situation 2: cluster node depends on front node and other cluster node
1874  vend = graph_->endEdges(vertex);
1875  for(Iterator edge = graph_->beginEdges(vertex); edge != vend; ++edge){
1876  //if(!edge.properties().depends() && edge.properties().influences()
1877  if(edge.properties().isStrong()
1878  && aggregates[edge.target()]==aggregate)
1879  {
1880  // Search for a link from target that stays within the aggregate
1881  Iterator v1end = graph_->endEdges(edge.target());
1882 
1883  for(Iterator edge1=graph_->beginEdges(edge.target()); edge1 != v1end; ++edge1){
1884  //if(edge1.properties().depends() && !edge1.properties().influences()
1885  if(edge1.properties().isStrong()
1886  && aggregates[edge1.target()]==aggregate)
1887  {
1888  bool found=false;
1889  // Check if front node is also connected to this one
1890  Iterator v2end = graph_->endEdges(vertex);
1891  for(Iterator edge2 = graph_->beginEdges(vertex); edge2 != v2end; ++edge2){
1892  if(edge2.target()==edge1.target()){
1893  if(edge2.properties().isStrong())
1894  found=true;
1895  break;
1896  }
1897  }
1898  if(found)
1899  return true;
1900  }
1901  }
1902  }
1903  }
1904  return false;
1905  }
1906 
1907  template<class G>
1908  void Aggregator<G>::unmarkFront()
1909  {
1910  typedef typename VertexList::const_iterator Iterator;
1911 
1912  for(Iterator vertex=front_.begin(); vertex != front_.end(); ++vertex)
1913  graph_->getVertexProperties(*vertex).resetFront();
1914 
1915  front_.clear();
1916  }
1917 
1918  template<class G>
1919  inline void
1920  Aggregator<G>::nonisoNeighbourAggregate(const Vertex& vertex,
1921  const AggregatesMap<Vertex>& aggregates,
1922  SLList<Vertex>& neighbours) const
1923  {
1924  typedef typename MatrixGraph::ConstEdgeIterator Iterator;
1925  Iterator end=graph_->beginEdges(vertex);
1926  neighbours.clear();
1927 
1928  for(Iterator edge=graph_->beginEdges(vertex); edge!=end; ++edge)
1929  {
1930  if(aggregates[edge.target()]!=AggregatesMap<Vertex>::UNAGGREGATED && graph_->getVertexProperties(edge.target()).isolated())
1931  neighbours.push_back(aggregates[edge.target()]);
1932  }
1933  }
1934 
1935  template<class G>
1936  inline typename G::VertexDescriptor Aggregator<G>::mergeNeighbour(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates) const
1937  {
1938  typedef typename MatrixGraph::ConstEdgeIterator Iterator;
1939 
1940  Iterator end = graph_->endEdges(vertex);
1941  for(Iterator edge = graph_->beginEdges(vertex); edge != end; ++edge){
1942  if(aggregates[edge.target()] != AggregatesMap<Vertex>::UNAGGREGATED &&
1943  graph_->getVertexProperties(edge.target()).isolated() == graph_->getVertexProperties(edge.source()).isolated()){
1944  if( graph_->getVertexProperties(vertex).isolated() ||
1945  ((edge.properties().depends() || edge.properties().influences())
1946  && admissible(vertex, aggregates[edge.target()], aggregates)))
1947  return edge.target();
1948  }
1949  }
1951  }
1952 
1953  template<class G>
1954  Aggregator<G>::FrontNeighbourCounter::FrontNeighbourCounter(const MatrixGraph& graph)
1955  : Counter(), graph_(graph)
1956  {}
1957 
1958  template<class G>
1959  void Aggregator<G>::FrontNeighbourCounter::operator()(const typename MatrixGraph::ConstEdgeIterator& edge)
1960  {
1961  if(graph_.getVertexProperties(edge.target()).front())
1962  Counter::increment();
1963  }
1964 
1965  template<class G>
1966  int Aggregator<G>::noFrontNeighbours(const Vertex& vertex) const
1967  {
1968  FrontNeighbourCounter counter(*graph_);
1969  visitNeighbours(*graph_, vertex, counter);
1970  return counter.value();
1971  }
1972  template<class G>
1973  inline bool Aggregator<G>::connected(const Vertex& vertex,
1974  const AggregateDescriptor& aggregate,
1975  const AggregatesMap<Vertex>& aggregates) const
1976  {
1977  typedef typename G::ConstEdgeIterator iterator;
1978  const iterator end = graph_->endEdges(vertex);
1979  for(iterator edge = graph_->beginEdges(vertex); edge != end; ++edge)
1980  if(aggregates[edge.target()]==aggregate)
1981  return true;
1982  return false;
1983  }
1984  template<class G>
1985  inline bool Aggregator<G>::connected(const Vertex& vertex,
1986  const SLList<AggregateDescriptor>& aggregateList,
1987  const AggregatesMap<Vertex>& aggregates) const
1988  {
1989  typedef typename SLList<AggregateDescriptor>::const_iterator Iter;
1990  for(Iter i=aggregateList.begin(); i!=aggregateList.end(); ++i)
1991  if(connected(vertex, *i, aggregates))
1992  return true;
1993  return false;
1994  }
1995 
1996  template<class G>
1997  template<class C>
1998  void Aggregator<G>::growIsolatedAggregate(const Vertex& seed, const AggregatesMap<Vertex>& aggregates, const C& c)
1999  {
2000  SLList<Vertex> connectedAggregates;
2001  nonisoNeighbourAggregate(seed, aggregates,connectedAggregates);
2002 
2003  while(aggregate_->size()< c.minAggregateSize() && aggregate_->connectSize() < c.maxConnectivity()){
2004  double maxCon=-1;
2005  std::size_t maxFrontNeighbours=0;
2006 
2008  unmarkFront();
2009  markFront(aggregates);
2010  typedef typename VertexList::const_iterator Iterator;
2011 
2012  for(Iterator vertex = front_.begin(); vertex != front_.end(); ++vertex){
2013  if(distance(*vertex, aggregates)>c.maxDistance())
2014  continue; // distance of proposes aggregate too big
2015 
2016  if(connectedAggregates.size()>0){
2017  // there is already a neighbour cluster
2018  // front node must be connected to same neighbour cluster
2019 
2020  if(!connected(*vertex, connectedAggregates, aggregates))
2021  continue;
2022  }
2023 
2024  double con = connectivity(*vertex, aggregates);
2025 
2026  if(con == maxCon){
2027  std::size_t frontNeighbours = noFrontNeighbours(*vertex);
2028 
2029  if(frontNeighbours >= maxFrontNeighbours){
2030  maxFrontNeighbours = frontNeighbours;
2031  candidate = *vertex;
2032  }
2033  }else if(con > maxCon){
2034  maxCon = con;
2035  maxFrontNeighbours = noFrontNeighbours(*vertex);
2036  candidate = *vertex;
2037  }
2038  }
2039 
2041  break;
2042 
2043  aggregate_->add(candidate);
2044  }
2045  }
2046 
2047  template<class G>
2048  template<class C>
2049  void Aggregator<G>::growAggregate(const Vertex& seed, const AggregatesMap<Vertex>& aggregates, const C& c)
2050  {
2051  while(aggregate_->size() < c.minAggregateSize()){
2052  int maxTwoCons=0, maxOneCons=0, maxNeighbours=-1;
2053  double maxCon=-1;
2054 
2056 
2057  unmarkFront();
2058  markFront(aggregates);
2059 
2060  typedef typename VertexList::const_iterator Iterator;
2061 
2062  for(Iterator vertex = front_.begin(); vertex != front_.end(); ++vertex){
2063  // Only nonisolated nodes are considered
2064  if(graph_->getVertexProperties(*vertex).isolated())
2065  continue;
2066 
2067  int twoWayCons = twoWayConnections(*vertex, aggregate_->id(), aggregates);
2068 
2069  /* The case of two way connections. */
2070  if( maxTwoCons == twoWayCons && twoWayCons > 0){
2071  double con = connectivity(*vertex, aggregates);
2072 
2073  if(con == maxCon){
2074  int neighbours = noFrontNeighbours(*vertex);
2075 
2076  if(neighbours > maxNeighbours){
2077  maxNeighbours = neighbours;
2078 
2079  std::size_t distance_ = distance(*vertex, aggregates);
2080 
2081  if(c.maxDistance() >= distance_){
2082  candidate = *vertex;
2083  }
2084  }
2085  }else if( con > maxCon){
2086  maxCon = con;
2087  maxNeighbours = noFrontNeighbours(*vertex);
2088  std::size_t distance_ = distance(*vertex, aggregates);
2089 
2090  if(c.maxDistance() >= distance_){
2091  candidate = *vertex;
2092  }
2093  }
2094  }else if(twoWayCons > maxTwoCons){
2095  maxTwoCons = twoWayCons;
2096  maxCon = connectivity(*vertex, aggregates);
2097  maxNeighbours = noFrontNeighbours(*vertex);
2098  std::size_t distance_ = distance(*vertex, aggregates);
2099 
2100  if(c.maxDistance() >= distance_){
2101  candidate = *vertex;
2102  }
2103 
2104  // two way connections preceed
2105  maxOneCons = std::numeric_limits<int>::max();
2106  }
2107 
2108  if(twoWayCons > 0)
2109  continue; // THis is a two-way node, skip tests for one way nodes
2110 
2111  /* The one way case */
2112  int oneWayCons = oneWayConnections(*vertex, aggregate_->id(), aggregates);
2113 
2114  if(oneWayCons==0)
2115  continue; // No strong connections, skip the tests.
2116 
2117  if(!admissible(*vertex, aggregate_->id(), aggregates))
2118  continue;
2119 
2120  if( maxOneCons == oneWayCons && oneWayCons > 0){
2121  double con = connectivity(*vertex, aggregates);
2122 
2123  if(con == maxCon){
2124  int neighbours = noFrontNeighbours(*vertex);
2125 
2126  if(neighbours > maxNeighbours){
2127  maxNeighbours = neighbours;
2128  std::size_t distance_ = distance(*vertex, aggregates);
2129 
2130  if(c.maxDistance() >= distance_){
2131  candidate = *vertex;
2132  }
2133  }
2134  }else if( con > maxCon){
2135  maxCon = con;
2136  maxNeighbours = noFrontNeighbours(*vertex);
2137  std::size_t distance_ = distance(*vertex, aggregates);
2138  if(c.maxDistance() >= distance_){
2139  candidate = *vertex;
2140  }
2141  }
2142  }else if(oneWayCons > maxOneCons){
2143  maxOneCons = oneWayCons;
2144  maxCon = connectivity(*vertex, aggregates);
2145  maxNeighbours = noFrontNeighbours(*vertex);
2146  std::size_t distance_ = distance(*vertex, aggregates);
2147 
2148  if(c.maxDistance() >= distance_){
2149  candidate = *vertex;
2150  }
2151  }
2152  }
2153 
2154 
2155  if(candidate == AggregatesMap<Vertex>::UNAGGREGATED)
2156  break; // No more candidates found
2157 
2158  aggregate_->add(candidate);
2159  }
2160  }
2161 
2162  template<typename V>
2163  template<typename M, typename G, typename C>
2164  tuple<int,int,int,int> AggregatesMap<V>::buildAggregates(const M& matrix, G& graph, const C& criterion,
2165  bool finestLevel)
2166  {
2167  Aggregator<G> aggregator;
2168  return aggregator.build(matrix, graph, *this, criterion, finestLevel);
2169  }
2170 
2171  template<class G>
2172  template<class M, class C>
2173  tuple<int,int,int,int> Aggregator<G>::build(const M& m, G& graph, AggregatesMap<Vertex>& aggregates, const C& c,
2174  bool finestLevel)
2175  {
2176  // Stack for fast vertex access
2177  Stack stack_(graph, *this, aggregates);
2178 
2179  graph_ = &graph;
2180 
2181  aggregate_ = new Aggregate<G,VertexSet>(graph, aggregates, connected_);
2182 
2183  Timer watch;
2184  watch.reset();
2185 
2186  buildDependency(graph, m, c, finestLevel);
2187 
2188  dverb<<"Build dependency took "<< watch.elapsed()<<" seconds."<<std::endl;
2189  int noAggregates, conAggregates, isoAggregates, oneAggregates;
2190  std::size_t maxA=0, minA=1000000, avg=0;
2191  int skippedAggregates;
2192  noAggregates = conAggregates = isoAggregates = oneAggregates =
2193  skippedAggregates = 0;
2194 
2195  while(true){
2196  Vertex seed = stack_.pop();
2197 
2198  if(seed == Stack::NullEntry)
2199  // No more unaggregated vertices. We are finished!
2200  break;
2201 
2202  // Debugging output
2203  if((noAggregates+1)%10000 == 0)
2204  Dune::dverb<<"c";
2205 
2206  aggregate_->seed(seed);
2207 
2208  if(graph.getVertexProperties(seed).excludedBorder()){
2209  aggregates[seed]=AggregatesMap<Vertex>::ISOLATED;
2210  ++skippedAggregates;
2211  continue;
2212  }
2213 
2214  if(graph.getVertexProperties(seed).isolated()){
2215  if(c.skipIsolated()){
2216  // isolated vertices are not aggregated but skipped on the coarser levels.
2217  aggregates[seed]=AggregatesMap<Vertex>::ISOLATED;
2218  ++skippedAggregates;
2219  // skip rest as no agglomeration is done.
2220  continue;
2221  }else
2222  growIsolatedAggregate(seed, aggregates, c);
2223  }else
2224  growAggregate(seed, aggregates, c);
2225 
2226 
2227  /* The rounding step. */
2228  while(!(graph.getVertexProperties(seed).isolated()) && aggregate_->size() < c.maxAggregateSize()){
2229 
2230  unmarkFront();
2231  markFront(aggregates);
2232 
2234 
2235  typedef typename VertexList::const_iterator Iterator;
2236 
2237  for(Iterator vertex = front_.begin(); vertex != front_.end(); ++vertex){
2238 
2239  if(graph.getVertexProperties(*vertex).isolated())
2240  continue; // No isolated nodes here
2241 
2242  if(twoWayConnections( *vertex, aggregate_->id(), aggregates) == 0 &&
2243  (oneWayConnections( *vertex, aggregate_->id(), aggregates) == 0 ||
2244  !admissible( *vertex, aggregate_->id(), aggregates) ))
2245  continue;
2246 
2247  std::pair<int,int> neighbourPair=neighbours(*vertex, aggregate_->id(),
2248  aggregates);
2249 
2250  //if(aggregateNeighbours(*vertex, aggregate_->id(), aggregates) <= unusedNeighbours(*vertex, aggregates))
2251  // continue;
2252 
2253  if(neighbourPair.first >= neighbourPair.second)
2254  continue;
2255 
2256  if(distance(*vertex, aggregates) > c.maxDistance())
2257  continue; // Distance too far
2258  candidate = *vertex;
2259  break;
2260  }
2261 
2262  if(candidate == AggregatesMap<Vertex>::UNAGGREGATED) break; // no more candidates found.
2263 
2264  aggregate_->add(candidate);
2265 
2266  }
2267 
2268  // try to merge aggregates consisting of only one nonisolated vertex with other aggregates
2269  if(aggregate_->size()==1 && c.maxAggregateSize()>1){
2270  if(!graph.getVertexProperties(seed).isolated()){
2271  Vertex mergedNeighbour = mergeNeighbour(seed, aggregates);
2272 
2273  if(mergedNeighbour != AggregatesMap<Vertex>::UNAGGREGATED){
2274  // assign vertex to the neighbouring cluster
2275  aggregates[seed] = aggregates[mergedNeighbour];
2276  }else{
2277  minA=std::min(minA,static_cast<std::size_t>(1));
2278  maxA=std::max(maxA,static_cast<std::size_t>(1));
2279  ++oneAggregates;
2280  ++conAggregates;
2281  }
2282  }else{
2283  ++avg;
2284  minA=std::min(minA,static_cast<std::size_t>(1));
2285  maxA=std::max(maxA,static_cast<std::size_t>(1));
2286  ++oneAggregates;
2287  ++isoAggregates;
2288  }
2289  ++avg;
2290  }else{
2291  avg+=aggregate_->size();
2292  minA=std::min(minA,aggregate_->size());
2293  maxA=std::max(maxA,aggregate_->size());
2294  if(graph.getVertexProperties(seed).isolated())
2295  ++isoAggregates;
2296  else
2297  ++conAggregates;
2298  }
2299  unmarkFront();
2300  markFront(aggregates);
2301  seedFromFront(stack_, graph.getVertexProperties(seed).isolated());
2302  unmarkFront();
2303  }
2304 
2305  Dune::dinfo<<"connected aggregates: "<<conAggregates;
2306  Dune::dinfo<<" isolated aggregates: "<<isoAggregates;
2307  if(conAggregates+isoAggregates>0)
2308  Dune::dinfo<<" one node aggregates: "<<oneAggregates<<" min size="
2309  <<minA<<" max size="<<maxA
2310  <<" avg="<<avg/(conAggregates+isoAggregates)<<std::endl;
2311 
2312  delete aggregate_;
2313  return make_tuple(conAggregates+isoAggregates,isoAggregates,
2314  oneAggregates,skippedAggregates);
2315  }
2316 
2317  template<class G>
2318  inline void Aggregator<G>::seedFromFront(Stack& stack_, bool isolated)
2319  {
2320  typedef typename VertexList::const_iterator Iterator;
2321 
2322  Iterator end= front_.end();
2323  int count=0;
2324  for(Iterator vertex=front_.begin(); vertex != end; ++vertex,++count)
2325  if(graph_->getVertexProperties(*vertex).isolated()==isolated)
2326  stack_.push(*vertex);
2327  /*
2328  if(MINIMAL_DEBUG_LEVEL<=2 && count==0 && !isolated)
2329  Dune::dverb<< " no vertices pushed for nonisolated aggregate!"<<std::endl;
2330  */
2331  }
2332 
2333  template<class G>
2334  Aggregator<G>::Stack::Stack(const MatrixGraph& graph, const Aggregator<G>& aggregatesBuilder,
2335  const AggregatesMap<Vertex>& aggregates)
2336  : graph_(graph), aggregatesBuilder_(aggregatesBuilder), aggregates_(aggregates), size_(0), maxSize_(0), head_(0), filled_(0)
2337  {
2338  vals_ = new Vertex[N];
2339  }
2340 
2341  template<class G>
2342  Aggregator<G>::Stack::~Stack()
2343  {
2344  Dune::dverb << "Max stack size was "<<maxSize_<<" filled="<<filled_<<std::endl;
2345  delete[] vals_;
2346  }
2347 
2348  template<class G>
2349  const typename Aggregator<G>::Vertex Aggregator<G>::Stack::NullEntry
2350  = std::numeric_limits<typename G::VertexDescriptor>::max();
2351 
2352  template<class G>
2353  inline bool Aggregator<G>::Stack::push(const Vertex& v)
2354  {
2355  if(aggregates_[v] == AggregatesMap<Vertex>::UNAGGREGATED){
2356  localPush(v);
2357  return true;
2358  }else
2359  return false;
2360  }
2361 
2362  template<class G>
2363  inline void Aggregator<G>::Stack::localPush(const Vertex& v)
2364  {
2365  vals_[head_] = v;
2366  size_ = std::min<int>(size_+1, N);
2367  head_ = (head_+N+1)%N;
2368  }
2369 
2370  template<class G>
2371  void Aggregator<G>::Stack::fill()
2372  {
2373  int isolated = 0, connected=0;
2374  int isoumin, umin;
2375  filled_++;
2376 
2377  head_ = size_ = 0;
2378  isoumin = umin = std::numeric_limits<int>::max();
2379 
2380  typedef typename MatrixGraph::ConstVertexIterator Iterator;
2381 
2382  const Iterator end = graph_.end();
2383 
2384  for(Iterator vertex = graph_.begin(); vertex != end; ++vertex){
2385  // Skip already aggregated vertices
2386  if(aggregates_[*vertex] != AggregatesMap<Vertex>::UNAGGREGATED)
2387  continue;
2388 
2389  if(vertex.properties().isolated()){
2390  isoumin = std::min(isoumin, aggregatesBuilder_.unusedNeighbours(*vertex, aggregates_));
2391  isolated++;
2392  }else{
2393  umin = std::min(umin, aggregatesBuilder_.unusedNeighbours(*vertex, aggregates_));
2394  connected++;
2395  }
2396  }
2397 
2398  if(connected + isolated == 0)
2399  // No unaggregated vertices.
2400  return;
2401 
2402  if(connected > 0){
2403  // Connected vertices have higher priority.
2404  for(Iterator vertex = graph_.begin(); vertex != end; ++vertex)
2405  if(aggregates_[*vertex] == AggregatesMap<Vertex>::UNAGGREGATED && !vertex.properties().isolated()
2406  && aggregatesBuilder_.unusedNeighbours(*vertex, aggregates_) == umin)
2407  localPush(*vertex);
2408  }else{
2409  for(Iterator vertex = graph_.begin(); vertex != end; ++vertex)
2410  if(aggregates_[*vertex] == AggregatesMap<Vertex>::UNAGGREGATED && vertex.properties().isolated()
2411  && aggregatesBuilder_.unusedNeighbours(*vertex, aggregates_) == isoumin)
2412  localPush(*vertex);
2413  }
2414  maxSize_ = std::max(size_, maxSize_);
2415  }
2416 
2417  template<class G>
2418  inline typename G::VertexDescriptor Aggregator<G>::Stack::pop()
2419  {
2420  while(size_>0){
2421  head_ = (head_ + N -1) % N;
2422  size_--;
2423  Vertex v = vals_[head_];
2424  if(aggregates_[v]==AggregatesMap<Vertex>::UNAGGREGATED)
2425  return v;
2426  }
2427  // Stack is empty try to fill it
2428  fill();
2429 
2430  // try again
2431  while(size_>0){
2432  head_ = (head_ + N -1) % N;
2433  size_--;
2434  Vertex v = vals_[head_];
2435  if(aggregates_[v]==AggregatesMap<Vertex>::UNAGGREGATED)
2436  return v;
2437  }
2438  return NullEntry;
2439  }
2440 
2441 #endif // DOXYGEN
2442 
2443  template<class V>
2444  void printAggregates2d(const AggregatesMap<V>& aggregates, int n, int m, std::ostream& os)
2445  {
2446  std::ios_base::fmtflags oldOpts=os.flags();
2447 
2448  os.setf(std::ios_base::right, std::ios_base::adjustfield);
2449 
2450  V max=0;
2451  int width=1;
2452 
2453  for(int i=0; i< n*m; i++)
2454  max=std::max(max, aggregates[i]);
2455 
2456  for(int i=10; i < 1000000; i*=10)
2457  if(max/i>0)
2458  width++;
2459  else
2460  break;
2461 
2462  for(int j=0, entry=0; j < m; j++){
2463  for(int i=0; i<n; i++, entry++){
2464  os.width(width);
2465  os<<aggregates[entry]<<" ";
2466  }
2467 
2468  os<<std::endl;
2469  }
2470  os<<std::endl;
2471  os.flags(oldOpts);
2472  }
2473 
2474 
2475  }// namespace Amg
2476 
2477 }// namespace Dune
2478 
2479 
2480 #endif