Dune Core Modules (2.6.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>
15#include <dune/common/sllist.hh>
16#include <dune/common/unused.hh>
18
19#include <utility>
20#include <set>
21#include <algorithm>
22#include <complex>
23#include <limits>
24#include <ostream>
25#include <tuple>
26
27namespace Dune
28{
29 namespace Amg
30 {
31
45 template<class T>
46 class AggregationCriterion : public T
47 {
48
49 public:
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:
178 typedef typename FieldTraits<field_type>::real_type real_type;
179 real_type maxValue_;
183 int row_;
185 real_type diagonal_;
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 vals_.assign(row.size(), 0.0);
202 assert(vals_.size()==row.size());
203 valIter_=vals_.begin();
204
205 maxValue_ = std::min(- std::numeric_limits<real_type>::max(), std::numeric_limits<real_type>::min());
206 diagonal_=norm_(row[index]);
207 row_ = index;
208 }
209
210 template<class M, class N>
211 inline void SymmetricMatrixDependency<M,N>::examine(const ColIter& col)
212 {
213 // skip positive offdiagonals if norm preserves sign of them.
214 real_type eij = norm_(*col);
215 if(!N::is_sign_preserving || eij<0) // || eji<0)
216 {
217 *valIter_ = eij/diagonal_*eij/norm_(matrix_->operator[](col.index())[col.index()]);
218 maxValue_ = std::max(maxValue_, *valIter_);
219 }else
220 *valIter_ =0;
221 ++valIter_;
222 }
223
224 template<class M, class N>
225 template<class G>
226 inline void SymmetricMatrixDependency<M,N>::examine(G&, const typename G::EdgeIterator& edge, const ColIter&)
227 {
228 if(*valIter_ > alpha() * maxValue_) {
229 edge.properties().setDepends();
230 edge.properties().setInfluences();
231 }
232 ++valIter_;
233 }
234
235 template<class M, class N>
236 inline bool SymmetricMatrixDependency<M,N>::isIsolated()
237 {
238 if(diagonal_==0)
239 DUNE_THROW(Dune::ISTLError, "No diagonal entry for row "<<row_<<".");
240 valIter_=vals_.begin();
241 return maxValue_ < beta();
242 }
243
247 template<class M, class N>
248 class Dependency : public Parameters
249 {
250 public:
254 typedef M Matrix;
255
259 typedef N Norm;
260
264 typedef typename Matrix::row_type Row;
265
270
271 void init(const Matrix* matrix);
272
273 void initRow(const Row& row, int index);
274
275 void examine(const ColIter& col);
276
277 template<class G>
278 void examine(G& graph, const typename G::EdgeIterator& edge, const ColIter& col);
279
280 bool isIsolated();
281
282 Dependency(const Parameters& parms)
283 : Parameters(parms)
284 {}
285
286 Dependency()
287 : Parameters()
288 {}
289
290 protected:
295 typedef typename FieldTraits<field_type>::real_type real_type;
296 real_type maxValue_;
300 int row_;
302 real_type diagonal_;
303 };
304
308 template<class M, class N>
310 {
311 public:
315 typedef M Matrix;
316
320 typedef N Norm;
321
325 typedef typename Matrix::row_type Row;
326
331
332 void init(const Matrix* matrix);
333
334 void initRow(const Row& row, int index);
335
336 void examine(const ColIter& col);
337
338 template<class G>
339 void examine(G& graph, const typename G::EdgeIterator& edge, const ColIter& col);
340
341 bool isIsolated();
342
343
344 SymmetricDependency(const Parameters& parms)
345 : Parameters(parms)
346 {}
348 : Parameters()
349 {}
350
351 protected:
356 typedef typename FieldTraits<field_type>::real_type real_type;
357 real_type maxValue_;
361 int row_;
363 real_type diagonal_;
364 };
365
370 template<int N>
372 {
373 public:
374 enum { /* @brief We preserve the sign.*/
375 is_sign_preserving = true
376 };
377
382 template<class M>
383 typename FieldTraits<typename M::field_type>::real_type operator()(const M& m) const
384 {
385 typedef typename M::field_type field_type;
386 typedef typename FieldTraits<field_type>::real_type real_type;
387 static_assert( std::is_convertible<field_type, real_type >::value,
388 "use of diagonal norm in AMG not implemented for complex field_type");
389 return m[N][N];
390 // possible implementation for complex types: return signed_abs(m[N][N]);
391 }
392
393 private:
394
396 template<typename T>
397 static T signed_abs(const T & v)
398 {
399 return v;
400 }
401
403 template<typename T>
404 static T signed_abs(const std::complex<T> & v)
405 {
406 // return sign * abs_value
407 // in case of complex numbers this extends to using the csgn function to determine the sign
408 return csgn(v) * std::abs(v);
409 }
410
412 template<typename T>
413 static T csgn(const T & v)
414 {
415 return (T(0) < v) - (v < T(0));
416 }
417
419 template<typename T>
420 static T csgn(std::complex<T> a)
421 {
422 return csgn(a.real())+(a.real() == 0.0)*csgn(a.imag());
423 }
424
425 };
426
431 class FirstDiagonal : public Diagonal<0>
432 {};
433
439 struct RowSum
440 {
441
442 enum { /* @brief We preserve the sign.*/
443 is_sign_preserving = false
444 };
449 template<class M>
450 typename FieldTraits<typename M::field_type>::real_type operator()(const M& m) const
451 {
452 return m.infinity_norm();
453 }
454 };
455
456 struct FrobeniusNorm
457 {
458
459 enum { /* @brief We preserve the sign.*/
460 is_sign_preserving = false
461 };
466 template<class M>
467 typename FieldTraits<typename M::field_type>::real_type operator()(const M& m) const
468 {
469 return m.frobenius_norm();
470 }
471 };
472 struct AlwaysOneNorm
473 {
474
475 enum { /* @brief We preserve the sign.*/
476 is_sign_preserving = false
477 };
482 template<class M>
483 typename FieldTraits<typename M::field_type>::real_type operator()(const M& m) const
484 {
485 return 1;
486 }
487 };
497 template<class M, class Norm>
498 class SymmetricCriterion : public AggregationCriterion<SymmetricDependency<M,Norm> >
499 {
500 public:
501 SymmetricCriterion(const Parameters& parms)
503 {}
505 {}
506 };
507
508
520 template<class M, class Norm>
521 class UnSymmetricCriterion : public AggregationCriterion<Dependency<M,Norm> >
522 {
523 public:
524 UnSymmetricCriterion(const Parameters& parms)
526 {}
528 {}
529 };
530 // forward declaration
531 template<class G> class Aggregator;
532
533
541 template<class V>
543 {
544 public:
545
549 static const V UNAGGREGATED;
550
554 static const V ISOLATED;
559
564
570
576
581 {
582 public:
583 template<class EdgeIterator>
584 void operator()(const EdgeIterator& edge) const
585 {
587 }
588 };
589
590
595
602
607
619 template<class M, class G, class C>
620 std::tuple<int,int,int,int> buildAggregates(const M& matrix, G& graph, const C& criterion,
621 bool finestLevel);
622
642 template<bool reset, class G, class F, class VM>
643 std::size_t breadthFirstSearch(const VertexDescriptor& start,
644 const AggregateDescriptor& aggregate,
645 const G& graph,
646 F& aggregateVisitor,
647 VM& visitedMap) const;
648
672 template<bool remove, bool reset, class G, class L, class F1, class F2, class VM>
673 std::size_t breadthFirstSearch(const VertexDescriptor& start,
674 const AggregateDescriptor& aggregate,
675 const G& graph, L& visited, F1& aggregateVisitor,
676 F2& nonAggregateVisitor,
677 VM& visitedMap) const;
678
684 void allocate(std::size_t noVertices);
685
689 std::size_t noVertices() const;
690
694 void free();
695
702
709
710 typedef const AggregateDescriptor* const_iterator;
711
712 const_iterator begin() const
713 {
714 return aggregates_;
715 }
716
717 const_iterator end() const
718 {
719 return aggregates_+noVertices();
720 }
721
722 typedef AggregateDescriptor* iterator;
723
724 iterator begin()
725 {
726 return aggregates_;
727 }
728
729 iterator end()
730 {
731 return aggregates_+noVertices();
732 }
733 private:
735 AggregatesMap(const AggregatesMap<V>&) = delete;
737 AggregatesMap<V>& operator=(const AggregatesMap<V>&) = delete;
738
742 AggregateDescriptor* aggregates_;
743
747 std::size_t noVertices_;
748 };
749
753 template<class G, class C>
754 void buildDependency(G& graph,
755 const typename C::Matrix& matrix,
756 C criterion,
757 bool finestLevel);
758
763 template<class G, class S>
765 {
766
767 public:
768
769 /***
770 * @brief The type of the matrix graph we work with.
771 */
772 typedef G MatrixGraph;
777
783
788 typedef S VertexSet;
789
791 typedef typename VertexSet::const_iterator const_iterator;
792
796 typedef std::size_t* SphereMap;
797
806 Aggregate(MatrixGraph& graph, AggregatesMap<Vertex>& aggregates,
807 VertexSet& connectivity, std::vector<Vertex>& front_);
808
809 void invalidate()
810 {
811 --id_;
812 }
813
821
825 void seed(const Vertex& vertex);
826
830 void add(const Vertex& vertex);
831
832 void add(std::vector<Vertex>& vertex);
836 void clear();
837
841 typename VertexSet::size_type size();
845 typename VertexSet::size_type connectSize();
846
850 int id();
851
854
857
858 private:
862 VertexSet vertices_;
863
868 int id_;
869
873 MatrixGraph& graph_;
874
878 AggregatesMap<Vertex>& aggregates_;
879
883 VertexSet& connected_;
884
888 std::vector<Vertex>& front_;
889 };
890
894 template<class G>
896 {
897 public:
898
902 typedef G MatrixGraph;
903
908
911
916
921
938 template<class M, class C>
939 std::tuple<int,int,int,int> build(const M& m, G& graph,
940 AggregatesMap<Vertex>& aggregates, const C& c,
941 bool finestLevel);
942 private:
948
953
957 typedef std::set<Vertex,std::less<Vertex>,Allocator> VertexSet;
958
962 typedef std::size_t* SphereMap;
963
967 MatrixGraph* graph_;
968
973
977 std::vector<Vertex> front_;
978
982 VertexSet connected_;
983
987 int size_;
988
992 class Stack
993 {
994 public:
995 static const Vertex NullEntry;
996
997 Stack(const MatrixGraph& graph,
998 const Aggregator<G>& aggregatesBuilder,
999 const AggregatesMap<Vertex>& aggregates);
1000 ~Stack();
1001 Vertex pop();
1002 private:
1003 enum { N = 1300000 };
1004
1006 const MatrixGraph& graph_;
1008 const Aggregator<G>& aggregatesBuilder_;
1010 const AggregatesMap<Vertex>& aggregates_;
1012 int size_;
1013 Vertex maxSize_;
1015 typename MatrixGraph::ConstVertexIterator begin_;
1017
1019 Vertex* vals_;
1020
1021 };
1022
1023 friend class Stack;
1024
1035 template<class V>
1036 void visitAggregateNeighbours(const Vertex& vertex, const AggregateDescriptor& aggregate,
1037 const AggregatesMap<Vertex>& aggregates,
1038 V& visitor) const;
1039
1044 template<class V>
1045 class AggregateVisitor
1046 {
1047 public:
1051 typedef V Visitor;
1060 Visitor& visitor);
1061
1068 void operator()(const typename MatrixGraph::ConstEdgeIterator& edge);
1069
1070 private:
1072 const AggregatesMap<Vertex>& aggregates_;
1074 AggregateDescriptor aggregate_;
1076 Visitor* visitor_;
1077 };
1078
1082 class Counter
1083 {
1084 public:
1088 int value();
1089
1090 protected:
1095
1096 private:
1097 int count_;
1098 };
1099
1100
1105 class FrontNeighbourCounter : public Counter
1106 {
1107 public:
1113
1114 void operator()(const typename MatrixGraph::ConstEdgeIterator& edge);
1115
1116 private:
1117 const MatrixGraph& graph_;
1118 };
1119
1124 int noFrontNeighbours(const Vertex& vertex) const;
1125
1129 class TwoWayCounter : public Counter
1130 {
1131 public:
1132 void operator()(const typename MatrixGraph::ConstEdgeIterator& edge);
1133 };
1134
1146 int twoWayConnections(const Vertex&, const AggregateDescriptor& aggregate,
1147 const AggregatesMap<Vertex>& aggregates) const;
1148
1152 class OneWayCounter : public Counter
1153 {
1154 public:
1155 void operator()(const typename MatrixGraph::ConstEdgeIterator& edge);
1156 };
1157
1169 int oneWayConnections(const Vertex&, const AggregateDescriptor& aggregate,
1170 const AggregatesMap<Vertex>& aggregates) const;
1171
1178 class ConnectivityCounter : public Counter
1179 {
1180 public:
1187 ConnectivityCounter(const VertexSet& connected, const AggregatesMap<Vertex>& aggregates);
1188
1189 void operator()(const typename MatrixGraph::ConstEdgeIterator& edge);
1190
1191 private:
1193 const VertexSet& connected_;
1195 const AggregatesMap<Vertex>& aggregates_;
1196
1197 };
1198
1210 double connectivity(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates) const;
1218 bool connected(const Vertex& vertex, const AggregateDescriptor& aggregate,
1219 const AggregatesMap<Vertex>& aggregates) const;
1220
1228 bool connected(const Vertex& vertex, const SLList<AggregateDescriptor>& aggregateList,
1229 const AggregatesMap<Vertex>& aggregates) const;
1230
1238 class DependencyCounter : public Counter
1239 {
1240 public:
1245
1246 void operator()(const typename MatrixGraph::ConstEdgeIterator& edge);
1247 };
1248
1255 class FrontMarker
1256 {
1257 public:
1264 FrontMarker(std::vector<Vertex>& front, MatrixGraph& graph);
1265
1266 void operator()(const typename MatrixGraph::ConstEdgeIterator& edge);
1267
1268 private:
1270 std::vector<Vertex>& front_;
1272 MatrixGraph& graph_;
1273 };
1274
1278 void unmarkFront();
1279
1294 int unusedNeighbours(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates) const;
1295
1309 std::pair<int,int> neighbours(const Vertex& vertex,
1310 const AggregateDescriptor& aggregate,
1311 const AggregatesMap<Vertex>& aggregates) const;
1328 int aggregateNeighbours(const Vertex& vertex, const AggregateDescriptor& aggregate, const AggregatesMap<Vertex>& aggregates) const;
1329
1337 bool admissible(const Vertex& vertex, const AggregateDescriptor& aggregate, const AggregatesMap<Vertex>& aggregates) const;
1338
1346 std::size_t distance(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates);
1347
1356 Vertex mergeNeighbour(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates) const;
1357
1366 void nonisoNeighbourAggregate(const Vertex& vertex,
1367 const AggregatesMap<Vertex>& aggregates,
1368 SLList<Vertex>& neighbours) const;
1369
1377 template<class C>
1378 void growAggregate(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates, const C& c);
1379 template<class C>
1380 void growIsolatedAggregate(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates, const C& c);
1381 };
1382
1383#ifndef DOXYGEN
1384
1385 template<class M, class N>
1386 inline void SymmetricDependency<M,N>::init(const Matrix* matrix)
1387 {
1388 matrix_ = matrix;
1389 }
1390
1391 template<class M, class N>
1392 inline void SymmetricDependency<M,N>::initRow(const Row& row, int index)
1393 {
1395 maxValue_ = std::min(- std::numeric_limits<typename Matrix::field_type>::max(), std::numeric_limits<typename Matrix::field_type>::min());
1396 row_ = index;
1397 diagonal_ = norm_(matrix_->operator[](row_)[row_]);
1398 }
1399
1400 template<class M, class N>
1401 inline void SymmetricDependency<M,N>::examine(const ColIter& col)
1402 {
1403 real_type eij = norm_(*col);
1404 typename Matrix::ConstColIterator opposite_entry =
1405 matrix_->operator[](col.index()).find(row_);
1406 if ( opposite_entry == matrix_->operator[](col.index()).end() )
1407 {
1408 // Consider this a weak connection we disregard.
1409 return;
1410 }
1411 real_type eji = norm_(*opposite_entry);
1412
1413 // skip positive offdiagonals if norm preserves sign of them.
1414 if(!N::is_sign_preserving || eij<0 || eji<0)
1415 maxValue_ = std::max(maxValue_,
1416 eij /diagonal_ * eji/
1417 norm_(matrix_->operator[](col.index())[col.index()]));
1418 }
1419
1420 template<class M, class N>
1421 template<class G>
1422 inline void SymmetricDependency<M,N>::examine(G& graph, const typename G::EdgeIterator& edge, const ColIter& col)
1423 {
1424 real_type eij = norm_(*col);
1425 typename Matrix::ConstColIterator opposite_entry =
1426 matrix_->operator[](col.index()).find(row_);
1427
1428 if ( opposite_entry == matrix_->operator[](col.index()).end() )
1429 {
1430 // Consider this as a weak connection we disregard.
1431 return;
1432 }
1433 real_type eji = norm_(*opposite_entry);
1434 // skip positve offdiagonals if norm preserves sign of them.
1435 if(!N::is_sign_preserving || (eij<0 || eji<0))
1436 if(eji / norm_(matrix_->operator[](edge.target())[edge.target()]) *
1437 eij/ diagonal_ > alpha() * maxValue_) {
1438 edge.properties().setDepends();
1439 edge.properties().setInfluences();
1440 typename G::EdgeProperties& other = graph.getEdgeProperties(edge.target(), edge.source());
1441 other.setInfluences();
1442 other.setDepends();
1443 }
1444 }
1445
1446 template<class M, class N>
1447 inline bool SymmetricDependency<M,N>::isIsolated()
1448 {
1449 return maxValue_ < beta();
1450 }
1451
1452
1453 template<class M, class N>
1454 inline void Dependency<M,N>::init(const Matrix* matrix)
1455 {
1456 matrix_ = matrix;
1457 }
1458
1459 template<class M, class N>
1460 inline void Dependency<M,N>::initRow(const Row& row, int index)
1461 {
1463 maxValue_ = std::min(- std::numeric_limits<real_type>::max(), std::numeric_limits<real_type>::min());
1464 row_ = index;
1465 diagonal_ = norm_(matrix_->operator[](row_)[row_]);
1466 }
1467
1468 template<class M, class N>
1469 inline void Dependency<M,N>::examine(const ColIter& col)
1470 {
1471 maxValue_ = std::max(maxValue_,
1472 -norm_(*col));
1473 }
1474
1475 template<class M, class N>
1476 template<class G>
1477 inline void Dependency<M,N>::examine(G& graph, const typename G::EdgeIterator& edge, const ColIter& col)
1478 {
1479 if(-norm_(*col) >= maxValue_ * alpha()) {
1480 edge.properties().setDepends();
1481 typedef typename G::EdgeDescriptor ED;
1482 ED e= graph.findEdge(edge.target(), edge.source());
1483 if(e!=std::numeric_limits<ED>::max())
1484 {
1485 typename G::EdgeProperties& other = graph.getEdgeProperties(e);
1486 other.setInfluences();
1487 }
1488 }
1489 }
1490
1491 template<class M, class N>
1492 inline bool Dependency<M,N>::isIsolated()
1493 {
1494 return maxValue_ < beta() * diagonal_;
1495 }
1496
1497 template<class G,class S>
1498 Aggregate<G,S>::Aggregate(MatrixGraph& graph, AggregatesMap<Vertex>& aggregates,
1499 VertexSet& connected, std::vector<Vertex>& front)
1500 : vertices_(), id_(-1), graph_(graph), aggregates_(aggregates),
1501 connected_(connected), front_(front)
1502 {}
1503
1504 template<class G,class S>
1506 {
1507 /*
1508 vertices_.push_back(vertex);
1509 typedef typename VertexList::const_iterator iterator;
1510 iterator begin = vertices_.begin();
1511 iterator end = vertices_.end();*/
1512 throw "Not yet implemented";
1513
1514 // while(begin!=end){
1515 //for();
1516 // }
1517
1518 }
1519
1520 template<class G,class S>
1521 inline void Aggregate<G,S>::seed(const Vertex& vertex)
1522 {
1523 dvverb<<"Connected cleared"<<std::endl;
1524 connected_.clear();
1525 vertices_.clear();
1526 connected_.insert(vertex);
1527 dvverb << " Inserting "<<vertex<<" size="<<connected_.size();
1528 ++id_ ;
1529 add(vertex);
1530 }
1531
1532
1533 template<class G,class S>
1534 inline void Aggregate<G,S>::add(const Vertex& vertex)
1535 {
1536 vertices_.insert(vertex);
1537 aggregates_[vertex]=id_;
1538 if(front_.size())
1539 front_.erase(std::lower_bound(front_.begin(), front_.end(), vertex));
1540
1541
1542 typedef typename MatrixGraph::ConstEdgeIterator iterator;
1543 const iterator end = graph_.endEdges(vertex);
1544 for(iterator edge = graph_.beginEdges(vertex); edge != end; ++edge) {
1545 dvverb << " Inserting "<<aggregates_[edge.target()];
1546 connected_.insert(aggregates_[edge.target()]);
1547 dvverb <<" size="<<connected_.size();
1548 if(aggregates_[edge.target()]==AggregatesMap<Vertex>::UNAGGREGATED &&
1549 !graph_.getVertexProperties(edge.target()).front())
1550 {
1551 front_.push_back(edge.target());
1552 graph_.getVertexProperties(edge.target()).setFront();
1553 }
1554 }
1555 dvverb <<std::endl;
1556 std::sort(front_.begin(), front_.end());
1557 }
1558
1559 template<class G,class S>
1560 inline void Aggregate<G,S>::add(std::vector<Vertex>& vertices)
1561 {
1562#ifndef NDEBUG
1563 std::size_t oldsize = vertices_.size();
1564#endif
1565 typedef typename std::vector<Vertex>::iterator Iterator;
1566
1567 typedef typename VertexSet::iterator SIterator;
1568
1569 SIterator pos=vertices_.begin();
1570 std::vector<Vertex> newFront;
1571 newFront.reserve(front_.capacity());
1572
1573 std::set_difference(front_.begin(), front_.end(), vertices.begin(), vertices.end(),
1574 std::back_inserter(newFront));
1575 front_=newFront;
1576
1577 for(Iterator vertex=vertices.begin(); vertex != vertices.end(); ++vertex)
1578 {
1579 pos=vertices_.insert(pos,*vertex);
1580 vertices_.insert(*vertex);
1581 graph_.getVertexProperties(*vertex).resetFront(); // Not a front node any more.
1582 aggregates_[*vertex]=id_;
1583
1584 typedef typename MatrixGraph::ConstEdgeIterator iterator;
1585 const iterator end = graph_.endEdges(*vertex);
1586 for(iterator edge = graph_.beginEdges(*vertex); edge != end; ++edge) {
1587 dvverb << " Inserting "<<aggregates_[edge.target()];
1588 connected_.insert(aggregates_[edge.target()]);
1589 if(aggregates_[edge.target()]==AggregatesMap<Vertex>::UNAGGREGATED &&
1590 !graph_.getVertexProperties(edge.target()).front())
1591 {
1592 front_.push_back(edge.target());
1593 graph_.getVertexProperties(edge.target()).setFront();
1594 }
1595 dvverb <<" size="<<connected_.size();
1596 }
1597 dvverb <<std::endl;
1598 }
1599 std::sort(front_.begin(), front_.end());
1600 assert(oldsize+vertices.size()==vertices_.size());
1601 }
1602 template<class G,class S>
1603 inline void Aggregate<G,S>::clear()
1604 {
1605 vertices_.clear();
1606 connected_.clear();
1607 id_=-1;
1608 }
1609
1610 template<class G,class S>
1611 inline typename Aggregate<G,S>::VertexSet::size_type
1613 {
1614 return vertices_.size();
1615 }
1616
1617 template<class G,class S>
1618 inline typename Aggregate<G,S>::VertexSet::size_type
1620 {
1621 return connected_.size();
1622 }
1623
1624 template<class G,class S>
1625 inline int Aggregate<G,S>::id()
1626 {
1627 return id_;
1628 }
1629
1630 template<class G,class S>
1632 {
1633 return vertices_.begin();
1634 }
1635
1636 template<class G,class S>
1638 {
1639 return vertices_.end();
1640 }
1641
1642 template<class V>
1643 const V AggregatesMap<V>::UNAGGREGATED = std::numeric_limits<V>::max();
1644
1645 template<class V>
1646 const V AggregatesMap<V>::ISOLATED = std::numeric_limits<V>::max()-1;
1647
1648 template<class V>
1650 : aggregates_(0)
1651 {}
1652
1653 template<class V>
1655 {
1656 if(aggregates_!=0)
1657 delete[] aggregates_;
1658 }
1659
1660
1661 template<class V>
1662 inline AggregatesMap<V>::AggregatesMap(std::size_t noVertices)
1663 {
1664 allocate(noVertices);
1665 }
1666
1667 template<class V>
1668 inline std::size_t AggregatesMap<V>::noVertices() const
1669 {
1670 return noVertices_;
1671 }
1672
1673 template<class V>
1674 inline void AggregatesMap<V>::allocate(std::size_t noVertices)
1675 {
1676 aggregates_ = new AggregateDescriptor[noVertices];
1677 noVertices_ = noVertices;
1678
1679 for(std::size_t i=0; i < noVertices; i++)
1680 aggregates_[i]=UNAGGREGATED;
1681 }
1682
1683 template<class V>
1684 inline void AggregatesMap<V>::free()
1685 {
1686 assert(aggregates_ != 0);
1687 delete[] aggregates_;
1688 aggregates_=0;
1689 }
1690
1691 template<class V>
1693 AggregatesMap<V>::operator[](const VertexDescriptor& v)
1694 {
1695 return aggregates_[v];
1696 }
1697
1698 template<class V>
1699 inline const typename AggregatesMap<V>::AggregateDescriptor&
1700 AggregatesMap<V>::operator[](const VertexDescriptor& v) const
1701 {
1702 return aggregates_[v];
1703 }
1704
1705 template<class V>
1706 template<bool reset, class G, class F,class VM>
1707 inline std::size_t AggregatesMap<V>::breadthFirstSearch(const V& start,
1708 const AggregateDescriptor& aggregate,
1709 const G& graph, F& aggregateVisitor,
1710 VM& visitedMap) const
1711 {
1712 VertexList vlist;
1713
1714 DummyEdgeVisitor dummy;
1715 return breadthFirstSearch<true,reset>(start, aggregate, graph, vlist, aggregateVisitor, dummy, visitedMap);
1716 }
1717
1718 template<class V>
1719 template<bool remove, bool reset, class G, class L, class F1, class F2, class VM>
1720 std::size_t AggregatesMap<V>::breadthFirstSearch(const V& start,
1721 const AggregateDescriptor& aggregate,
1722 const G& graph,
1723 L& visited,
1724 F1& aggregateVisitor,
1725 F2& nonAggregateVisitor,
1726 VM& visitedMap) const
1727 {
1728 typedef typename L::const_iterator ListIterator;
1729 int visitedSpheres = 0;
1730
1731 visited.push_back(start);
1732 put(visitedMap, start, true);
1733
1734 ListIterator current = visited.begin();
1735 ListIterator end = visited.end();
1736 std::size_t i=0, size=visited.size();
1737
1738 // visit the neighbours of all vertices of the
1739 // current sphere.
1740 while(current != end) {
1741
1742 for(; i<size; ++current, ++i) {
1743 typedef typename G::ConstEdgeIterator EdgeIterator;
1744 const EdgeIterator endEdge = graph.endEdges(*current);
1745
1746 for(EdgeIterator edge = graph.beginEdges(*current);
1747 edge != endEdge; ++edge) {
1748
1749 if(aggregates_[edge.target()]==aggregate) {
1750 if(!get(visitedMap, edge.target())) {
1751 put(visitedMap, edge.target(), true);
1752 visited.push_back(edge.target());
1753 aggregateVisitor(edge);
1754 }
1755 }else
1756 nonAggregateVisitor(edge);
1757 }
1758 }
1759 end = visited.end();
1760 size = visited.size();
1761 if(current != end)
1762 visitedSpheres++;
1763 }
1764
1765 if(reset)
1766 for(current = visited.begin(); current != end; ++current)
1767 put(visitedMap, *current, false);
1768
1769
1770 if(remove)
1771 visited.clear();
1772
1773 return visitedSpheres;
1774 }
1775
1776 template<class G>
1778 : graph_(0), aggregate_(0), front_(), connected_(), size_(-1)
1779 {}
1780
1781 template<class G>
1783 {
1784 size_=-1;
1785 }
1786
1787 template<class G, class C>
1788 void buildDependency(G& graph,
1789 const typename C::Matrix& matrix,
1790 C criterion, bool firstlevel)
1791 {
1792 // assert(graph.isBuilt());
1793 typedef typename C::Matrix Matrix;
1794 typedef typename G::VertexIterator VertexIterator;
1795
1796 criterion.init(&matrix);
1797
1798 for(VertexIterator vertex = graph.begin(); vertex != graph.end(); ++vertex) {
1799 typedef typename Matrix::row_type Row;
1800
1801 const Row& row = matrix[*vertex];
1802
1803 // Tell the criterion what row we will examine now
1804 // This might for example be used for calculating the
1805 // maximum offdiagonal value
1806 criterion.initRow(row, *vertex);
1807
1808 // On a first path all columns are examined. After this
1809 // the calculator should know whether the vertex is isolated.
1810 typedef typename Matrix::ConstColIterator ColIterator;
1811 ColIterator end = row.end();
1812 typename FieldTraits<typename Matrix::field_type>::real_type absoffdiag=0.;
1813
1814 if(firstlevel) {
1815 for(ColIterator col = row.begin(); col != end; ++col)
1816 if(col.index()!=*vertex) {
1817 criterion.examine(col);
1818 absoffdiag=std::max(absoffdiag, col->frobenius_norm());
1819 }
1820
1821 if(absoffdiag==0)
1822 vertex.properties().setExcludedBorder();
1823 }
1824 else
1825 for(ColIterator col = row.begin(); col != end; ++col)
1826 if(col.index()!=*vertex)
1827 criterion.examine(col);
1828
1829 // reset the vertex properties
1830 //vertex.properties().reset();
1831
1832 // Check whether the vertex is isolated.
1833 if(criterion.isIsolated()) {
1834 //std::cout<<"ISOLATED: "<<*vertex<<std::endl;
1835 vertex.properties().setIsolated();
1836 }else{
1837 // Examine all the edges beginning at this vertex.
1838 typedef typename G::EdgeIterator EdgeIterator;
1839 typedef typename Matrix::ConstColIterator ColIterator;
1840 EdgeIterator eEnd = vertex.end();
1841 ColIterator col = matrix[*vertex].begin();
1842
1843 for(EdgeIterator edge = vertex.begin(); edge!= eEnd; ++edge, ++col) {
1844 // Move to the right column.
1845 while(col.index()!=edge.target())
1846 ++col;
1847 criterion.examine(graph, edge, col);
1848 }
1849 }
1850
1851 }
1852 }
1853
1854
1855 template<class G>
1856 template<class V>
1857 inline Aggregator<G>::AggregateVisitor<V>::AggregateVisitor(const AggregatesMap<Vertex>& aggregates,
1858 const AggregateDescriptor& aggregate, V& visitor)
1859 : aggregates_(aggregates), aggregate_(aggregate), visitor_(&visitor)
1860 {}
1861
1862 template<class G>
1863 template<class V>
1864 inline void Aggregator<G>::AggregateVisitor<V>::operator()(const typename MatrixGraph::ConstEdgeIterator& edge)
1865 {
1866 if(aggregates_[edge.target()]==aggregate_)
1867 visitor_->operator()(edge);
1868 }
1869
1870 template<class G>
1871 template<class V>
1872 inline void Aggregator<G>::visitAggregateNeighbours(const Vertex& vertex,
1873 const AggregateDescriptor& aggregate,
1874 const AggregatesMap<Vertex>& aggregates,
1875 V& visitor) const
1876 {
1877 // Only evaluates for edge pointing to the aggregate
1878 AggregateVisitor<V> v(aggregates, aggregate, visitor);
1879 visitNeighbours(*graph_, vertex, v);
1880 }
1881
1882
1883 template<class G>
1884 inline Aggregator<G>::Counter::Counter()
1885 : count_(0)
1886 {}
1887
1888 template<class G>
1889 inline void Aggregator<G>::Counter::increment()
1890 {
1891 ++count_;
1892 }
1893
1894 template<class G>
1895 inline void Aggregator<G>::Counter::decrement()
1896 {
1897 --count_;
1898 }
1899 template<class G>
1900 inline int Aggregator<G>::Counter::value()
1901 {
1902 return count_;
1903 }
1904
1905 template<class G>
1906 inline void Aggregator<G>::TwoWayCounter::operator()(const typename MatrixGraph::ConstEdgeIterator& edge)
1907 {
1908 if(edge.properties().isTwoWay())
1909 Counter::increment();
1910 }
1911
1912 template<class G>
1913 int Aggregator<G>::twoWayConnections(const Vertex& vertex, const AggregateDescriptor& aggregate,
1914 const AggregatesMap<Vertex>& aggregates) const
1915 {
1916 TwoWayCounter counter;
1917 visitAggregateNeighbours(vertex, aggregate, aggregates, counter);
1918 return counter.value();
1919 }
1920
1921 template<class G>
1922 int Aggregator<G>::oneWayConnections(const Vertex& vertex, const AggregateDescriptor& aggregate,
1923 const AggregatesMap<Vertex>& aggregates) const
1924 {
1925 OneWayCounter counter;
1926 visitAggregateNeighbours(vertex, aggregate, aggregates, counter);
1927 return counter.value();
1928 }
1929
1930 template<class G>
1931 inline void Aggregator<G>::OneWayCounter::operator()(const typename MatrixGraph::ConstEdgeIterator& edge)
1932 {
1933 if(edge.properties().isOneWay())
1934 Counter::increment();
1935 }
1936
1937 template<class G>
1938 inline Aggregator<G>::ConnectivityCounter::ConnectivityCounter(const VertexSet& connected,
1939 const AggregatesMap<Vertex>& aggregates)
1940 : Counter(), connected_(connected), aggregates_(aggregates)
1941 {}
1942
1943
1944 template<class G>
1945 inline void Aggregator<G>::ConnectivityCounter::operator()(const typename MatrixGraph::ConstEdgeIterator& edge)
1946 {
1947 if(connected_.find(aggregates_[edge.target()]) == connected_.end() || aggregates_[edge.target()]==AggregatesMap<Vertex>::UNAGGREGATED)
1948 // Would be a new connection
1949 Counter::increment();
1950 else{
1951 Counter::increment();
1952 Counter::increment();
1953 }
1954 }
1955
1956 template<class G>
1957 inline double Aggregator<G>::connectivity(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates) const
1958 {
1959 ConnectivityCounter counter(connected_, aggregates);
1960 double noNeighbours=visitNeighbours(*graph_, vertex, counter);
1961 return (double)counter.value()/noNeighbours;
1962 }
1963
1964 template<class G>
1965 inline Aggregator<G>::DependencyCounter::DependencyCounter()
1966 : Counter()
1967 {}
1968
1969 template<class G>
1970 inline void Aggregator<G>::DependencyCounter::operator()(const typename MatrixGraph::ConstEdgeIterator& edge)
1971 {
1972 if(edge.properties().depends())
1973 Counter::increment();
1974 if(edge.properties().influences())
1975 Counter::increment();
1976 }
1977
1978 template<class G>
1979 int Aggregator<G>::unusedNeighbours(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates) const
1980 {
1981 return aggregateNeighbours(vertex, AggregatesMap<Vertex>::UNAGGREGATED, aggregates);
1982 }
1983
1984 template<class G>
1985 std::pair<int,int> Aggregator<G>::neighbours(const Vertex& vertex,
1986 const AggregateDescriptor& aggregate,
1987 const AggregatesMap<Vertex>& aggregates) const
1988 {
1989 DependencyCounter unused, aggregated;
1990 typedef AggregateVisitor<DependencyCounter> Counter;
1991 typedef std::tuple<Counter,Counter> CounterTuple;
1992 CombinedFunctor<CounterTuple> visitors(CounterTuple(Counter(aggregates, AggregatesMap<Vertex>::UNAGGREGATED, unused), Counter(aggregates, aggregate, aggregated)));
1993 visitNeighbours(*graph_, vertex, visitors);
1994 return std::make_pair(unused.value(), aggregated.value());
1995 }
1996
1997
1998 template<class G>
1999 int Aggregator<G>::aggregateNeighbours(const Vertex& vertex, const AggregateDescriptor& aggregate, const AggregatesMap<Vertex>& aggregates) const
2000 {
2001 DependencyCounter counter;
2002 visitAggregateNeighbours(vertex, aggregate, aggregates, counter);
2003 return counter.value();
2004 }
2005
2006 template<class G>
2007 std::size_t Aggregator<G>::distance(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates)
2008 {
2009 return 0;
2010 typename PropertyMapTypeSelector<VertexVisitedTag,G>::Type visitedMap = get(VertexVisitedTag(), *graph_);
2011 VertexList vlist;
2012 typename AggregatesMap<Vertex>::DummyEdgeVisitor dummy;
2013 return aggregates.template breadthFirstSearch<true,true>(vertex,
2014 aggregate_->id(), *graph_,
2015 vlist, dummy, dummy, visitedMap);
2016 }
2017
2018 template<class G>
2019 inline Aggregator<G>::FrontMarker::FrontMarker(std::vector<Vertex>& front, MatrixGraph& graph)
2020 : front_(front), graph_(graph)
2021 {}
2022
2023 template<class G>
2024 inline void Aggregator<G>::FrontMarker::operator()(const typename MatrixGraph::ConstEdgeIterator& edge)
2025 {
2026 Vertex target = edge.target();
2027
2028 if(!graph_.getVertexProperties(target).front()) {
2029 front_.push_back(target);
2030 graph_.getVertexProperties(target).setFront();
2031 }
2032 }
2033
2034 template<class G>
2035 inline bool Aggregator<G>::admissible(const Vertex& vertex, const AggregateDescriptor& aggregate, const AggregatesMap<Vertex>& aggregates) const
2036 {
2037 // Todo
2038 Dune::dvverb<<" Admissible not yet implemented!"<<std::endl;
2039 return true;
2040 //Situation 1: front node depends on two nodes. Then these
2041 // have to be strongly connected to each other
2042
2043 // Iterate over all all neighbours of front node
2044 typedef typename MatrixGraph::ConstEdgeIterator Iterator;
2045 Iterator vend = graph_->endEdges(vertex);
2046 for(Iterator edge = graph_->beginEdges(vertex); edge != vend; ++edge) {
2047 // if(edge.properties().depends() && !edge.properties().influences()
2048 if(edge.properties().isStrong()
2049 && aggregates[edge.target()]==aggregate)
2050 {
2051 // Search for another link to the aggregate
2052 Iterator edge1 = edge;
2053 for(++edge1; edge1 != vend; ++edge1) {
2054 //if(edge1.properties().depends() && !edge1.properties().influences()
2055 if(edge1.properties().isStrong()
2056 && aggregates[edge.target()]==aggregate)
2057 {
2058 //Search for an edge connecting the two vertices that is
2059 //strong
2060 bool found=false;
2061 Iterator v2end = graph_->endEdges(edge.target());
2062 for(Iterator edge2 = graph_->beginEdges(edge.target()); edge2 != v2end; ++edge2) {
2063 if(edge2.target()==edge1.target() &&
2064 edge2.properties().isStrong()) {
2065 found =true;
2066 break;
2067 }
2068 }
2069 if(found)
2070 {
2071 return true;
2072 }
2073 }
2074 }
2075 }
2076 }
2077
2078 // Situation 2: cluster node depends on front node and other cluster node
2080 vend = graph_->endEdges(vertex);
2081 for(Iterator edge = graph_->beginEdges(vertex); edge != vend; ++edge) {
2082 //if(!edge.properties().depends() && edge.properties().influences()
2083 if(edge.properties().isStrong()
2084 && aggregates[edge.target()]==aggregate)
2085 {
2086 // Search for a link from target that stays within the aggregate
2087 Iterator v1end = graph_->endEdges(edge.target());
2088
2089 for(Iterator edge1=graph_->beginEdges(edge.target()); edge1 != v1end; ++edge1) {
2090 //if(edge1.properties().depends() && !edge1.properties().influences()
2091 if(edge1.properties().isStrong()
2092 && aggregates[edge1.target()]==aggregate)
2093 {
2094 bool found=false;
2095 // Check if front node is also connected to this one
2096 Iterator v2end = graph_->endEdges(vertex);
2097 for(Iterator edge2 = graph_->beginEdges(vertex); edge2 != v2end; ++edge2) {
2098 if(edge2.target()==edge1.target()) {
2099 if(edge2.properties().isStrong())
2100 found=true;
2101 break;
2102 }
2103 }
2104 if(found)
2105 {
2106 return true;
2107 }
2108 }
2109 }
2110 }
2111 }
2112 return false;
2113 }
2114
2115 template<class G>
2116 void Aggregator<G>::unmarkFront()
2117 {
2118 typedef typename std::vector<Vertex>::const_iterator Iterator;
2119
2120 for(Iterator vertex=front_.begin(); vertex != front_.end(); ++vertex)
2121 graph_->getVertexProperties(*vertex).resetFront();
2122
2123 front_.clear();
2124 }
2125
2126 template<class G>
2127 inline void
2128 Aggregator<G>::nonisoNeighbourAggregate(const Vertex& vertex,
2129 const AggregatesMap<Vertex>& aggregates,
2130 SLList<Vertex>& neighbours) const
2131 {
2132 typedef typename MatrixGraph::ConstEdgeIterator Iterator;
2133 Iterator end=graph_->beginEdges(vertex);
2134 neighbours.clear();
2135
2136 for(Iterator edge=graph_->beginEdges(vertex); edge!=end; ++edge)
2137 {
2138 if(aggregates[edge.target()]!=AggregatesMap<Vertex>::UNAGGREGATED && graph_->getVertexProperties(edge.target()).isolated())
2139 neighbours.push_back(aggregates[edge.target()]);
2140 }
2141 }
2142
2143 template<class G>
2144 inline typename G::VertexDescriptor Aggregator<G>::mergeNeighbour(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates) const
2145 {
2146 typedef typename MatrixGraph::ConstEdgeIterator Iterator;
2147
2148 Iterator end = graph_->endEdges(vertex);
2149 for(Iterator edge = graph_->beginEdges(vertex); edge != end; ++edge) {
2150 if(aggregates[edge.target()] != AggregatesMap<Vertex>::UNAGGREGATED &&
2151 graph_->getVertexProperties(edge.target()).isolated() == graph_->getVertexProperties(edge.source()).isolated()) {
2152 if( graph_->getVertexProperties(vertex).isolated() ||
2153 ((edge.properties().depends() || edge.properties().influences())
2154 && admissible(vertex, aggregates[edge.target()], aggregates)))
2155 return edge.target();
2156 }
2157 }
2159 }
2160
2161 template<class G>
2162 Aggregator<G>::FrontNeighbourCounter::FrontNeighbourCounter(const MatrixGraph& graph)
2163 : Counter(), graph_(graph)
2164 {}
2165
2166 template<class G>
2167 void Aggregator<G>::FrontNeighbourCounter::operator()(const typename MatrixGraph::ConstEdgeIterator& edge)
2168 {
2169 if(graph_.getVertexProperties(edge.target()).front())
2170 Counter::increment();
2171 }
2172
2173 template<class G>
2174 int Aggregator<G>::noFrontNeighbours(const Vertex& vertex) const
2175 {
2176 FrontNeighbourCounter counter(*graph_);
2177 visitNeighbours(*graph_, vertex, counter);
2178 return counter.value();
2179 }
2180 template<class G>
2181 inline bool Aggregator<G>::connected(const Vertex& vertex,
2182 const AggregateDescriptor& aggregate,
2183 const AggregatesMap<Vertex>& aggregates) const
2184 {
2185 typedef typename G::ConstEdgeIterator iterator;
2186 const iterator end = graph_->endEdges(vertex);
2187 for(iterator edge = graph_->beginEdges(vertex); edge != end; ++edge)
2188 if(aggregates[edge.target()]==aggregate)
2189 return true;
2190 return false;
2191 }
2192 template<class G>
2193 inline bool Aggregator<G>::connected(const Vertex& vertex,
2194 const SLList<AggregateDescriptor>& aggregateList,
2195 const AggregatesMap<Vertex>& aggregates) const
2196 {
2197 typedef typename SLList<AggregateDescriptor>::const_iterator Iter;
2198 for(Iter i=aggregateList.begin(); i!=aggregateList.end(); ++i)
2199 if(connected(vertex, *i, aggregates))
2200 return true;
2201 return false;
2202 }
2203
2204 template<class G>
2205 template<class C>
2206 void Aggregator<G>::growIsolatedAggregate(const Vertex& seed, const AggregatesMap<Vertex>& aggregates, const C& c)
2207 {
2208 SLList<Vertex> connectedAggregates;
2209 nonisoNeighbourAggregate(seed, aggregates,connectedAggregates);
2210
2211 while(aggregate_->size()< c.minAggregateSize() && aggregate_->connectSize() < c.maxConnectivity()) {
2212 double maxCon=-1;
2213 std::size_t maxFrontNeighbours=0;
2214
2216
2217 typedef typename std::vector<Vertex>::const_iterator Iterator;
2218
2219 for(Iterator vertex = front_.begin(); vertex != front_.end(); ++vertex) {
2220 if(distance(*vertex, aggregates)>c.maxDistance())
2221 continue; // distance of proposes aggregate too big
2222
2223 if(connectedAggregates.size()>0) {
2224 // there is already a neighbour cluster
2225 // front node must be connected to same neighbour cluster
2226
2227 if(!connected(*vertex, connectedAggregates, aggregates))
2228 continue;
2229 }
2230
2231 double con = connectivity(*vertex, aggregates);
2232
2233 if(con == maxCon) {
2234 std::size_t frontNeighbours = noFrontNeighbours(*vertex);
2235
2236 if(frontNeighbours >= maxFrontNeighbours) {
2237 maxFrontNeighbours = frontNeighbours;
2238 candidate = *vertex;
2239 }
2240 }else if(con > maxCon) {
2241 maxCon = con;
2242 maxFrontNeighbours = noFrontNeighbours(*vertex);
2243 candidate = *vertex;
2244 }
2245 }
2246
2248 break;
2249
2250 aggregate_->add(candidate);
2251 }
2252 }
2253
2254 template<class G>
2255 template<class C>
2256 void Aggregator<G>::growAggregate(const Vertex& seed, const AggregatesMap<Vertex>& aggregates, const C& c)
2257 {
2258
2259 std::size_t distance_ =0;
2260 while(aggregate_->size() < c.minAggregateSize()&& distance_<c.maxDistance()) {
2261 int maxTwoCons=0, maxOneCons=0, maxNeighbours=-1;
2262 double maxCon=-1;
2263
2264 std::vector<Vertex> candidates;
2265 candidates.reserve(30);
2266
2267 typedef typename std::vector<Vertex>::const_iterator Iterator;
2268
2269 for(Iterator vertex = front_.begin(); vertex != front_.end(); ++vertex) {
2270 // Only nonisolated nodes are considered
2271 if(graph_->getVertexProperties(*vertex).isolated())
2272 continue;
2273
2274 int twoWayCons = twoWayConnections(*vertex, aggregate_->id(), aggregates);
2275
2276 /* The case of two way connections. */
2277 if( maxTwoCons == twoWayCons && twoWayCons > 0) {
2278 double con = connectivity(*vertex, aggregates);
2279
2280 if(con == maxCon) {
2281 int neighbours = noFrontNeighbours(*vertex);
2282
2283 if(neighbours > maxNeighbours) {
2284 maxNeighbours = neighbours;
2285 candidates.clear();
2286 candidates.push_back(*vertex);
2287 }else{
2288 candidates.push_back(*vertex);
2289 }
2290 }else if( con > maxCon) {
2291 maxCon = con;
2292 maxNeighbours = noFrontNeighbours(*vertex);
2293 candidates.clear();
2294 candidates.push_back(*vertex);
2295 }
2296 }else if(twoWayCons > maxTwoCons) {
2297 maxTwoCons = twoWayCons;
2298 maxCon = connectivity(*vertex, aggregates);
2299 maxNeighbours = noFrontNeighbours(*vertex);
2300 candidates.clear();
2301 candidates.push_back(*vertex);
2302
2303 // two way connections precede
2304 maxOneCons = std::numeric_limits<int>::max();
2305 }
2306
2307 if(twoWayCons > 0)
2308 {
2309 continue; // THis is a two-way node, skip tests for one way nodes
2310 }
2311
2312 /* The one way case */
2313 int oneWayCons = oneWayConnections(*vertex, aggregate_->id(), aggregates);
2314
2315 if(oneWayCons==0)
2316 continue; // No strong connections, skip the tests.
2317
2318 if(!admissible(*vertex, aggregate_->id(), aggregates))
2319 continue;
2320
2321 if( maxOneCons == oneWayCons && oneWayCons > 0) {
2322 double con = connectivity(*vertex, aggregates);
2323
2324 if(con == maxCon) {
2325 int neighbours = noFrontNeighbours(*vertex);
2326
2327 if(neighbours > maxNeighbours) {
2328 maxNeighbours = neighbours;
2329 candidates.clear();
2330 candidates.push_back(*vertex);
2331 }else{
2332 if(neighbours==maxNeighbours)
2333 {
2334 candidates.push_back(*vertex);
2335 }
2336 }
2337 }else if( con > maxCon) {
2338 maxCon = con;
2339 maxNeighbours = noFrontNeighbours(*vertex);
2340 candidates.clear();
2341 candidates.push_back(*vertex);
2342 }
2343 }else if(oneWayCons > maxOneCons) {
2344 maxOneCons = oneWayCons;
2345 maxCon = connectivity(*vertex, aggregates);
2346 maxNeighbours = noFrontNeighbours(*vertex);
2347 candidates.clear();
2348 candidates.push_back(*vertex);
2349 }
2350 }
2351
2352
2353 if(!candidates.size())
2354 break; // No more candidates found
2355 distance_=distance(seed, aggregates);
2356 candidates.resize(std::min(candidates.size(), c.maxAggregateSize()-
2357 aggregate_->size()));
2358 aggregate_->add(candidates);
2359 }
2360 }
2361
2362 template<typename V>
2363 template<typename M, typename G, typename C>
2364 std::tuple<int,int,int,int> AggregatesMap<V>::buildAggregates(const M& matrix, G& graph, const C& criterion,
2365 bool finestLevel)
2366 {
2367 Aggregator<G> aggregator;
2368 return aggregator.build(matrix, graph, *this, criterion, finestLevel);
2369 }
2370
2371 template<class G>
2372 template<class M, class C>
2373 std::tuple<int,int,int,int> Aggregator<G>::build(const M& m, G& graph, AggregatesMap<Vertex>& aggregates, const C& c,
2374 bool finestLevel)
2375 {
2376 // Stack for fast vertex access
2377 Stack stack_(graph, *this, aggregates);
2378
2379 graph_ = &graph;
2380
2381 aggregate_ = new Aggregate<G,VertexSet>(graph, aggregates, connected_, front_);
2382
2383 Timer watch;
2384 watch.reset();
2385
2386 buildDependency(graph, m, c, finestLevel);
2387
2388 dverb<<"Build dependency took "<< watch.elapsed()<<" seconds."<<std::endl;
2389 int noAggregates, conAggregates, isoAggregates, oneAggregates;
2390 std::size_t maxA=0, minA=1000000, avg=0;
2391 int skippedAggregates;
2392 noAggregates = conAggregates = isoAggregates = oneAggregates =
2393 skippedAggregates = 0;
2394
2395 while(true) {
2396 Vertex seed = stack_.pop();
2397
2398 if(seed == Stack::NullEntry)
2399 // No more unaggregated vertices. We are finished!
2400 break;
2401
2402 // Debugging output
2403 if((noAggregates+1)%10000 == 0)
2404 Dune::dverb<<"c";
2405 unmarkFront();
2406
2407 if(graph.getVertexProperties(seed).excludedBorder()) {
2408 aggregates[seed]=AggregatesMap<Vertex>::ISOLATED;
2409 ++skippedAggregates;
2410 continue;
2411 }
2412
2413 if(graph.getVertexProperties(seed).isolated()) {
2414 if(c.skipIsolated()) {
2415 // isolated vertices are not aggregated but skipped on the coarser levels.
2416 aggregates[seed]=AggregatesMap<Vertex>::ISOLATED;
2417 ++skippedAggregates;
2418 // skip rest as no agglomeration is done.
2419 continue;
2420 }else{
2421 aggregate_->seed(seed);
2422 growIsolatedAggregate(seed, aggregates, c);
2423 }
2424 }else{
2425 aggregate_->seed(seed);
2426 growAggregate(seed, aggregates, c);
2427 }
2428
2429 /* The rounding step. */
2430 while(!(graph.getVertexProperties(seed).isolated()) && aggregate_->size() < c.maxAggregateSize()) {
2431
2432 std::vector<Vertex> candidates;
2433 candidates.reserve(30);
2434
2435 typedef typename std::vector<Vertex>::const_iterator Iterator;
2436
2437 for(Iterator vertex = front_.begin(); vertex != front_.end(); ++vertex) {
2438
2439 if(graph.getVertexProperties(*vertex).isolated())
2440 continue; // No isolated nodes here
2441
2442 if(twoWayConnections( *vertex, aggregate_->id(), aggregates) == 0 &&
2443 (oneWayConnections( *vertex, aggregate_->id(), aggregates) == 0 ||
2444 !admissible( *vertex, aggregate_->id(), aggregates) ))
2445 continue;
2446
2447 std::pair<int,int> neighbourPair=neighbours(*vertex, aggregate_->id(),
2448 aggregates);
2449
2450 //if(aggregateNeighbours(*vertex, aggregate_->id(), aggregates) <= unusedNeighbours(*vertex, aggregates))
2451 // continue;
2452
2453 if(neighbourPair.first >= neighbourPair.second)
2454 continue;
2455
2456 if(distance(*vertex, aggregates) > c.maxDistance())
2457 continue; // Distance too far
2458 candidates.push_back(*vertex);
2459 break;
2460 }
2461
2462 if(!candidates.size()) break; // no more candidates found.
2463
2464 candidates.resize(std::min(candidates.size(), c.maxAggregateSize()-
2465 aggregate_->size()));
2466 aggregate_->add(candidates);
2467
2468 }
2469
2470 // try to merge aggregates consisting of only one nonisolated vertex with other aggregates
2471 if(aggregate_->size()==1 && c.maxAggregateSize()>1) {
2472 if(!graph.getVertexProperties(seed).isolated()) {
2473 Vertex mergedNeighbour = mergeNeighbour(seed, aggregates);
2474
2475 if(mergedNeighbour != AggregatesMap<Vertex>::UNAGGREGATED) {
2476 // assign vertex to the neighbouring cluster
2477 aggregates[seed] = aggregates[mergedNeighbour];
2478 aggregate_->invalidate();
2479 }else{
2480 ++avg;
2481 minA=std::min(minA,static_cast<std::size_t>(1));
2482 maxA=std::max(maxA,static_cast<std::size_t>(1));
2483 ++oneAggregates;
2484 ++conAggregates;
2485 }
2486 }else{
2487 ++avg;
2488 minA=std::min(minA,static_cast<std::size_t>(1));
2489 maxA=std::max(maxA,static_cast<std::size_t>(1));
2490 ++oneAggregates;
2491 ++isoAggregates;
2492 }
2493 ++avg;
2494 }else{
2495 avg+=aggregate_->size();
2496 minA=std::min(minA,aggregate_->size());
2497 maxA=std::max(maxA,aggregate_->size());
2498 if(graph.getVertexProperties(seed).isolated())
2499 ++isoAggregates;
2500 else
2501 ++conAggregates;
2502 }
2503
2504 }
2505
2506 Dune::dinfo<<"connected aggregates: "<<conAggregates;
2507 Dune::dinfo<<" isolated aggregates: "<<isoAggregates;
2508 if(conAggregates+isoAggregates>0)
2509 Dune::dinfo<<" one node aggregates: "<<oneAggregates<<" min size="
2510 <<minA<<" max size="<<maxA
2511 <<" avg="<<avg/(conAggregates+isoAggregates)<<std::endl;
2512
2513 delete aggregate_;
2514 return std::make_tuple(conAggregates+isoAggregates,isoAggregates,
2515 oneAggregates,skippedAggregates);
2516 }
2517
2518
2519 template<class G>
2520 Aggregator<G>::Stack::Stack(const MatrixGraph& graph, const Aggregator<G>& aggregatesBuilder,
2521 const AggregatesMap<Vertex>& aggregates)
2522 : graph_(graph), aggregatesBuilder_(aggregatesBuilder), aggregates_(aggregates), begin_(graph.begin()), end_(graph.end())
2523 {
2524 //vals_ = new Vertex[N];
2525 }
2526
2527 template<class G>
2528 Aggregator<G>::Stack::~Stack()
2529 {
2530 //Dune::dverb << "Max stack size was "<<maxSize_<<" filled="<<filled_<<std::endl;
2531 //delete[] vals_;
2532 }
2533
2534 template<class G>
2535 const typename Aggregator<G>::Vertex Aggregator<G>::Stack::NullEntry
2536 = std::numeric_limits<typename G::VertexDescriptor>::max();
2537
2538 template<class G>
2539 inline typename G::VertexDescriptor Aggregator<G>::Stack::pop()
2540 {
2541 for(; begin_!=end_ && aggregates_[*begin_] != AggregatesMap<Vertex>::UNAGGREGATED; ++begin_) ;
2542
2543 if(begin_!=end_)
2544 {
2545 typename G::VertexDescriptor current=*begin_;
2546 ++begin_;
2547 return current;
2548 }else
2549 return NullEntry;
2550 }
2551
2552#endif // DOXYGEN
2553
2554 template<class V>
2555 void printAggregates2d(const AggregatesMap<V>& aggregates, int n, int m, std::ostream& os)
2556 {
2557 std::ios_base::fmtflags oldOpts=os.flags();
2558
2559 os.setf(std::ios_base::right, std::ios_base::adjustfield);
2560
2561 V max=0;
2562 int width=1;
2563
2564 for(int i=0; i< n*m; i++)
2565 max=std::max(max, aggregates[i]);
2566
2567 for(int i=10; i < 1000000; i*=10)
2568 if(max/i>0)
2569 width++;
2570 else
2571 break;
2572
2573 for(int j=0, entry=0; j < m; j++) {
2574 for(int i=0; i<n; i++, entry++) {
2575 os.width(width);
2576 os<<aggregates[entry]<<" ";
2577 }
2578
2579 os<<std::endl;
2580 }
2581 os<<std::endl;
2582 os.flags(oldOpts);
2583 }
2584
2585
2586 } // namespace Amg
2587
2588} // namespace Dune
2589
2590
2591#endif
A class for temporarily storing the vertices of an aggregate in.
Definition: aggregates.hh:765
A Dummy visitor that does nothing for each visited edge.
Definition: aggregates.hh:581
Class providing information about the mapping of the vertices onto aggregates.
Definition: aggregates.hh:543
Base class of all aggregation criterions.
Definition: aggregates.hh:47
Class for building the aggregates.
Definition: aggregates.hh:896
Dependency policy for symmetric matrices.
Definition: aggregates.hh:249
Norm that uses only the [N][N] entry of the block to determine couplings.
Definition: aggregates.hh:372
Norm that uses only the [0][0] entry of the block to determine couplings.
Definition: aggregates.hh:432
Iterator over all edges starting from a vertex.
Definition: graph.hh:93
The vertex iterator type of the graph.
Definition: graph.hh:207
The (undirected) graph of a matrix.
Definition: graph.hh:49
M::size_type VertexDescriptor
The vertex descriptor.
Definition: graph.hh:71
EdgeIteratorT< const MatrixGraph< Matrix > > ConstEdgeIterator
The constant edge iterator type.
Definition: graph.hh:296
All parameters for AMG.
Definition: parameters.hh:391
Criterion taking advantage of symmetric matrices.
Definition: aggregates.hh:499
Dependency policy for symmetric matrices.
Definition: aggregates.hh:310
Dependency policy for symmetric matrices.
Definition: aggregates.hh:132
Criterion suited for unsymmetric matrices.
Definition: aggregates.hh:522
derive error class from the base class in common
Definition: istlexception.hh:16
A generic dynamic dense matrix.
Definition: matrix.hh:555
T::field_type field_type
Export the type representing the underlying field.
Definition: matrix.hh:559
row_type::const_iterator ConstColIterator
Const iterator for the entries of each row.
Definition: matrix.hh:583
MatrixImp::DenseMatrixBase< T, A >::window_type row_type
The type implementing a matrix row.
Definition: matrix.hh:568
An allocator managing a pool of objects for reuse.
Definition: poolallocator.hh:247
A single linked list.
Definition: sllist.hh:42
Type traits to determine the type of reals (when working with complex numbers)
Provides classes for building the matrix graph.
SLListConstIterator< T, A > const_iterator
The constant iterator of the list.
Definition: sllist.hh:72
#define DUNE_UNUSED_PARAMETER(parm)
A macro to mark intentionally unused function parameters with.
Definition: unused.hh:25
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
constexpr GeometryType vertex
GeometryType representing a vertex.
Definition: type.hh:727
Matrix::ConstColIterator ColIter
Constant column iterator of the matrix.
Definition: aggregates.hh:269
Matrix::ConstColIterator ColIter
Constant column iterator of the matrix.
Definition: aggregates.hh:152
std::size_t breadthFirstSearch(const VertexDescriptor &start, const AggregateDescriptor &aggregate, const G &graph, L &visited, F1 &aggregateVisitor, F2 &nonAggregateVisitor, VM &visitedMap) const
Breadth first search within an aggregate.
PoolAllocator< VertexDescriptor, 100 > Allocator
The allocator we use for our lists and the set.
Definition: aggregates.hh:569
int id()
Get the id identifying the aggregate.
Norm norm_
The functor for calculating the norm.
Definition: aggregates.hh:298
FieldTraits< typenameM::field_type >::real_type operator()(const M &m) const
compute the norm of a matrix.
Definition: aggregates.hh:383
MatrixGraph::VertexDescriptor Vertex
The vertex identifier.
Definition: aggregates.hh:907
AggregationCriterion()
Constructor.
Definition: aggregates.hh:64
const Matrix * matrix_
The matrix we work on.
Definition: aggregates.hh:353
M Matrix
The matrix type we build the dependency of.
Definition: aggregates.hh:254
G MatrixGraph
The matrix graph type used.
Definition: aggregates.hh:902
Norm norm_
The functor for calculating the norm.
Definition: aggregates.hh:359
M Matrix
The matrix type we build the dependency of.
Definition: aggregates.hh:315
real_type diagonal_
The norm of the current diagonal.
Definition: aggregates.hh:185
N Norm
The norm to use for examining the matrix entries.
Definition: aggregates.hh:259
int row_
index of the currently evaluated row.
Definition: aggregates.hh:183
std::tuple< int, int, int, int > build(const M &m, G &graph, AggregatesMap< Vertex > &aggregates, const C &c, bool finestLevel)
Build the aggregates.
FrontNeighbourCounter(const MatrixGraph &front)
Constructor.
Matrix::row_type Row
Constant Row iterator of the matrix.
Definition: aggregates.hh:325
const Matrix * matrix_
The matrix we work on.
Definition: aggregates.hh:292
const AggregateDescriptor & operator[](const VertexDescriptor &v) const
Get the aggregate a vertex belongs to.
AggregateVisitor(const AggregatesMap< Vertex > &aggregates, const AggregateDescriptor &aggregate, Visitor &visitor)
Constructor.
Matrix::ConstColIterator ColIter
Constant column iterator of the matrix.
Definition: aggregates.hh:330
~AggregatesMap()
Destructor.
Matrix::field_type field_type
The current max value.
Definition: aggregates.hh:177
void decrement()
Decrement counter.
Aggregate(MatrixGraph &graph, AggregatesMap< Vertex > &aggregates, VertexSet &connectivity, std::vector< Vertex > &front_)
Constructor.
V Visitor
The type of the adapted visitor.
Definition: aggregates.hh:1051
std::size_t * SphereMap
Type of the mapping of aggregate members onto distance spheres.
Definition: aggregates.hh:796
Matrix::row_type Row
Constant Row iterator of the matrix.
Definition: aggregates.hh:264
VertexSet::size_type connectSize()
Get tne number of connections to other aggregates.
N Norm
The norm to use for examining the matrix entries.
Definition: aggregates.hh:320
Norm norm_
The functor for calculating the norm.
Definition: aggregates.hh:181
VertexSet::const_iterator const_iterator
Const iterator over a vertex list.
Definition: aggregates.hh:791
MatrixGraph::VertexDescriptor AggregateDescriptor
The type of the aggregate descriptor.
Definition: aggregates.hh:910
real_type diagonal_
The norm of the current diagonal.
Definition: aggregates.hh:363
~Aggregator()
Destructor.
void add(const Vertex &vertex)
Add a vertex to the aggregate.
T DependencyPolicy
The policy for calculating the dependency graph.
Definition: aggregates.hh:53
Aggregator()
Constructor.
void operator()(const typename MatrixGraph::ConstEdgeIterator &edge)
Examine an edge.
FrontMarker(std::vector< Vertex > &front, MatrixGraph &graph)
Constructor.
FieldTraits< typenameM::field_type >::real_type operator()(const M &m) const
compute the norm of a matrix.
Definition: aggregates.hh:483
int visitNeighbours(const G &graph, const typename G::VertexDescriptor &vertex, V &visitor)
Visit all neighbour vertices of a vertex in a graph.
void setDefaultValuesIsotropic(std::size_t dim, std::size_t diameter=2)
Sets reasonable default values for an isotropic problem.
Definition: aggregates.hh:80
V AggregateDescriptor
The aggregate descriptor type.
Definition: aggregates.hh:563
static const V ISOLATED
Identifier of isolated vertices.
Definition: aggregates.hh:554
int row_
index of the currently evaluated row.
Definition: aggregates.hh:361
real_type diagonal_
The norm of the current diagonal.
Definition: aggregates.hh:302
std::size_t noVertices() const
Get the number of vertices.
void setDefaultValuesAnisotropic(std::size_t dim, std::size_t diameter=2)
Sets reasonable default values for an aisotropic problem.
Definition: aggregates.hh:103
AggregatesMap(std::size_t noVertices)
Constructs with allocating memory.
Matrix::field_type field_type
The current max value.
Definition: aggregates.hh:355
AggregateDescriptor & operator[](const VertexDescriptor &v)
Get the aggregate a vertex belongs to.
AggregatesMap()
Constructs without allocating memory.
int value()
Access the current count.
SLList< VertexDescriptor, Allocator > VertexList
The type of a single linked list of vertex descriptors.
Definition: aggregates.hh:575
ConnectivityCounter(const VertexSet &connected, const AggregatesMap< Vertex > &aggregates)
Constructor.
VertexSet::size_type size()
Get the size of the aggregate.
const_iterator end() const
get an iterator over the vertices of the aggregate.
int row_
index of the currently evaluated row.
Definition: aggregates.hh:300
M Matrix
The matrix type we build the dependency of.
Definition: aggregates.hh:137
const Matrix * matrix_
The matrix we work on.
Definition: aggregates.hh:175
S VertexSet
The type of a single linked list of vertex descriptors.
Definition: aggregates.hh:788
FieldTraits< typenameM::field_type >::real_type operator()(const M &m) const
compute the norm of a matrix.
Definition: aggregates.hh:467
static const V UNAGGREGATED
Identifier of not yet aggregated vertices.
Definition: aggregates.hh:549
std::size_t breadthFirstSearch(const VertexDescriptor &start, const AggregateDescriptor &aggregate, const G &graph, F &aggregateVisitor, VM &visitedMap) const
Breadth first search within an aggregate.
Matrix::field_type field_type
The current max value.
Definition: aggregates.hh:294
void allocate(std::size_t noVertices)
Allocate memory for holding the information.
N Norm
The norm to use for examining the matrix entries.
Definition: aggregates.hh:142
FieldTraits< typenameM::field_type >::real_type operator()(const M &m) const
compute the norm of a matrix.
Definition: aggregates.hh:450
void reconstruct(const Vertex &vertex)
Reconstruct the aggregat from an seed node.
const_iterator begin() const
get an iterator over the vertices of the aggregate.
MatrixGraph::VertexDescriptor Vertex
The vertex descriptor type.
Definition: aggregates.hh:776
void seed(const Vertex &vertex)
Initialize the aggregate with one vertex.
void clear()
Clear the aggregate.
void free()
Free the allocated memory.
void increment()
Increment counter.
void buildDependency(G &graph, const typename C::Matrix &matrix, C criterion, bool finestLevel)
Build the dependency of the matrix graph.
V VertexDescriptor
The vertex descriptor type.
Definition: aggregates.hh:558
std::tuple< int, int, int, int > buildAggregates(const M &matrix, G &graph, const C &criterion, bool finestLevel)
Build the aggregates.
Matrix::row_type Row
Constant Row iterator of the matrix.
Definition: aggregates.hh:147
PoolAllocator< Vertex, 100 > Allocator
The allocator we use for our lists and the set.
Definition: aggregates.hh:782
DVVerbType dvverb(std::cout)
stream for very verbose output.
Definition: stdstreams.hh:93
DInfoType dinfo(std::cout)
Stream for informative output.
Definition: stdstreams.hh:138
DVerbType dverb(std::cout)
Singleton of verbose debug stream.
Definition: stdstreams.hh:114
constexpr Front front
PartitionSet for the front partition.
Definition: partitionset.hh:279
Dune namespace.
Definition: alignedallocator.hh:10
Parameter classes for customizing AMG.
An stl-compliant pool allocator.
Provides classes for handling internal properties in a graph.
Implements a singly linked list together with the necessary iterators.
Standard Dune debug streams.
Functor using the row sum (infinity) norm to determine strong couplings.
Definition: aggregates.hh:440
A simple timing class.
Definition of the DUNE_UNUSED macro for the case that config.h is not available.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 13, 23:29, 2024)