Dune Core Modules (2.5.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>& map)
736 {
737 throw "Auch!";
738 }
739
741 AggregatesMap<V>& operator=(const AggregatesMap<V>& map)
742 {
743 throw "Auch!";
744 return this;
745 }
746
750 AggregateDescriptor* aggregates_;
751
755 std::size_t noVertices_;
756 };
757
761 template<class G, class C>
762 void buildDependency(G& graph,
763 const typename C::Matrix& matrix,
764 C criterion,
765 bool finestLevel);
766
771 template<class G, class S>
773 {
774
775 public:
776
777 /***
778 * @brief The type of the matrix graph we work with.
779 */
780 typedef G MatrixGraph;
785
791
796 typedef S VertexSet;
797
799 typedef typename VertexSet::const_iterator const_iterator;
800
804 typedef std::size_t* SphereMap;
805
814 Aggregate(MatrixGraph& graph, AggregatesMap<Vertex>& aggregates,
815 VertexSet& connectivity, std::vector<Vertex>& front_);
816
817 void invalidate()
818 {
819 --id_;
820 }
821
828 void reconstruct(const Vertex& vertex);
829
833 void seed(const Vertex& vertex);
834
838 void add(const Vertex& vertex);
839
840 void add(std::vector<Vertex>& vertex);
844 void clear();
845
849 typename VertexSet::size_type size();
853 typename VertexSet::size_type connectSize();
854
858 int id();
859
862
865
866 private:
870 VertexSet vertices_;
871
876 int id_;
877
881 MatrixGraph& graph_;
882
886 AggregatesMap<Vertex>& aggregates_;
887
891 VertexSet& connected_;
892
896 std::vector<Vertex>& front_;
897 };
898
902 template<class G>
904 {
905 public:
906
910 typedef G MatrixGraph;
911
916
919
924
929
946 template<class M, class C>
947 std::tuple<int,int,int,int> build(const M& m, G& graph,
948 AggregatesMap<Vertex>& aggregates, const C& c,
949 bool finestLevel);
950 private:
956
961
965 typedef std::set<Vertex,std::less<Vertex>,Allocator> VertexSet;
966
970 typedef std::size_t* SphereMap;
971
975 MatrixGraph* graph_;
976
981
985 std::vector<Vertex> front_;
986
990 VertexSet connected_;
991
995 int size_;
996
1000 class Stack
1001 {
1002 public:
1003 static const Vertex NullEntry;
1004
1005 Stack(const MatrixGraph& graph,
1006 const Aggregator<G>& aggregatesBuilder,
1007 const AggregatesMap<Vertex>& aggregates);
1008 ~Stack();
1009 Vertex pop();
1010 private:
1011 enum { N = 1300000 };
1012
1014 const MatrixGraph& graph_;
1016 const Aggregator<G>& aggregatesBuilder_;
1018 const AggregatesMap<Vertex>& aggregates_;
1020 int size_;
1021 Vertex maxSize_;
1023 typename MatrixGraph::ConstVertexIterator begin_;
1025
1027 Vertex* vals_;
1028
1029 };
1030
1031 friend class Stack;
1032
1043 template<class V>
1044 void visitAggregateNeighbours(const Vertex& vertex, const AggregateDescriptor& aggregate,
1045 const AggregatesMap<Vertex>& aggregates,
1046 V& visitor) const;
1047
1052 template<class V>
1053 class AggregateVisitor
1054 {
1055 public:
1059 typedef V Visitor;
1068 Visitor& visitor);
1069
1076 void operator()(const typename MatrixGraph::ConstEdgeIterator& edge);
1077
1078 private:
1080 const AggregatesMap<Vertex>& aggregates_;
1082 AggregateDescriptor aggregate_;
1084 Visitor* visitor_;
1085 };
1086
1090 class Counter
1091 {
1092 public:
1096 int value();
1097
1098 protected:
1103
1104 private:
1105 int count_;
1106 };
1107
1108
1113 class FrontNeighbourCounter : public Counter
1114 {
1115 public:
1121
1122 void operator()(const typename MatrixGraph::ConstEdgeIterator& edge);
1123
1124 private:
1125 const MatrixGraph& graph_;
1126 };
1127
1132 int noFrontNeighbours(const Vertex& vertex) const;
1133
1137 class TwoWayCounter : public Counter
1138 {
1139 public:
1140 void operator()(const typename MatrixGraph::ConstEdgeIterator& edge);
1141 };
1142
1154 int twoWayConnections(const Vertex&, const AggregateDescriptor& aggregate,
1155 const AggregatesMap<Vertex>& aggregates) const;
1156
1160 class OneWayCounter : public Counter
1161 {
1162 public:
1163 void operator()(const typename MatrixGraph::ConstEdgeIterator& edge);
1164 };
1165
1177 int oneWayConnections(const Vertex&, const AggregateDescriptor& aggregate,
1178 const AggregatesMap<Vertex>& aggregates) const;
1179
1186 class ConnectivityCounter : public Counter
1187 {
1188 public:
1195 ConnectivityCounter(const VertexSet& connected, const AggregatesMap<Vertex>& aggregates);
1196
1197 void operator()(const typename MatrixGraph::ConstEdgeIterator& edge);
1198
1199 private:
1201 const VertexSet& connected_;
1203 const AggregatesMap<Vertex>& aggregates_;
1204
1205 };
1206
1218 double connectivity(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates) const;
1226 bool connected(const Vertex& vertex, const AggregateDescriptor& aggregate,
1227 const AggregatesMap<Vertex>& aggregates) const;
1228
1236 bool connected(const Vertex& vertex, const SLList<AggregateDescriptor>& aggregateList,
1237 const AggregatesMap<Vertex>& aggregates) const;
1238
1246 class DependencyCounter : public Counter
1247 {
1248 public:
1253
1254 void operator()(const typename MatrixGraph::ConstEdgeIterator& edge);
1255 };
1256
1263 class FrontMarker
1264 {
1265 public:
1272 FrontMarker(std::vector<Vertex>& front, MatrixGraph& graph);
1273
1274 void operator()(const typename MatrixGraph::ConstEdgeIterator& edge);
1275
1276 private:
1278 std::vector<Vertex>& front_;
1280 MatrixGraph& graph_;
1281 };
1282
1286 void unmarkFront();
1287
1302 int unusedNeighbours(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates) const;
1303
1317 std::pair<int,int> neighbours(const Vertex& vertex,
1318 const AggregateDescriptor& aggregate,
1319 const AggregatesMap<Vertex>& aggregates) const;
1336 int aggregateNeighbours(const Vertex& vertex, const AggregateDescriptor& aggregate, const AggregatesMap<Vertex>& aggregates) const;
1337
1345 bool admissible(const Vertex& vertex, const AggregateDescriptor& aggregate, const AggregatesMap<Vertex>& aggregates) const;
1346
1354 std::size_t distance(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates);
1355
1364 Vertex mergeNeighbour(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates) const;
1365
1374 void nonisoNeighbourAggregate(const Vertex& vertex,
1375 const AggregatesMap<Vertex>& aggregates,
1376 SLList<Vertex>& neighbours) const;
1377
1385 template<class C>
1386 void growAggregate(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates, const C& c);
1387 template<class C>
1388 void growIsolatedAggregate(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates, const C& c);
1389 };
1390
1391#ifndef DOXYGEN
1392
1393 template<class M, class N>
1394 inline void SymmetricDependency<M,N>::init(const Matrix* matrix)
1395 {
1396 matrix_ = matrix;
1397 }
1398
1399 template<class M, class N>
1400 inline void SymmetricDependency<M,N>::initRow(const Row& row, int index)
1401 {
1403 maxValue_ = std::min(- std::numeric_limits<typename Matrix::field_type>::max(), std::numeric_limits<typename Matrix::field_type>::min());
1404 row_ = index;
1405 diagonal_ = norm_(matrix_->operator[](row_)[row_]);
1406 }
1407
1408 template<class M, class N>
1409 inline void SymmetricDependency<M,N>::examine(const ColIter& col)
1410 {
1411 real_type eij = norm_(*col);
1412 typename Matrix::ConstColIterator opposite_entry =
1413 matrix_->operator[](col.index()).find(row_);
1414 if ( opposite_entry == matrix_->operator[](col.index()).end() )
1415 {
1416 // Consider this a weak connection we disregard.
1417 return;
1418 }
1419 real_type eji = norm_(*opposite_entry);
1420
1421 // skip positive offdiagonals if norm preserves sign of them.
1422 if(!N::is_sign_preserving || eij<0 || eji<0)
1423 maxValue_ = std::max(maxValue_,
1424 eij /diagonal_ * eji/
1425 norm_(matrix_->operator[](col.index())[col.index()]));
1426 }
1427
1428 template<class M, class N>
1429 template<class G>
1430 inline void SymmetricDependency<M,N>::examine(G& graph, const typename G::EdgeIterator& edge, const ColIter& col)
1431 {
1432 real_type eij = norm_(*col);
1433 typename Matrix::ConstColIterator opposite_entry =
1434 matrix_->operator[](col.index()).find(row_);
1435
1436 if ( opposite_entry == matrix_->operator[](col.index()).end() )
1437 {
1438 // Consider this as a weak connection we disregard.
1439 return;
1440 }
1441 real_type eji = norm_(*opposite_entry);
1442 // skip positve offdiagonals if norm preserves sign of them.
1443 if(!N::is_sign_preserving || (eij<0 || eji<0))
1444 if(eji / norm_(matrix_->operator[](edge.target())[edge.target()]) *
1445 eij/ diagonal_ > alpha() * maxValue_) {
1446 edge.properties().setDepends();
1447 edge.properties().setInfluences();
1448 typename G::EdgeProperties& other = graph.getEdgeProperties(edge.target(), edge.source());
1449 other.setInfluences();
1450 other.setDepends();
1451 }
1452 }
1453
1454 template<class M, class N>
1455 inline bool SymmetricDependency<M,N>::isIsolated()
1456 {
1457 return maxValue_ < beta();
1458 }
1459
1460
1461 template<class M, class N>
1462 inline void Dependency<M,N>::init(const Matrix* matrix)
1463 {
1464 matrix_ = matrix;
1465 }
1466
1467 template<class M, class N>
1468 inline void Dependency<M,N>::initRow(const Row& row, int index)
1469 {
1471 maxValue_ = std::min(- std::numeric_limits<real_type>::max(), std::numeric_limits<real_type>::min());
1472 row_ = index;
1473 diagonal_ = norm_(matrix_->operator[](row_)[row_]);
1474 }
1475
1476 template<class M, class N>
1477 inline void Dependency<M,N>::examine(const ColIter& col)
1478 {
1479 maxValue_ = std::max(maxValue_,
1480 -norm_(*col));
1481 }
1482
1483 template<class M, class N>
1484 template<class G>
1485 inline void Dependency<M,N>::examine(G& graph, const typename G::EdgeIterator& edge, const ColIter& col)
1486 {
1487 if(-norm_(*col) >= maxValue_ * alpha()) {
1488 edge.properties().setDepends();
1489 typedef typename G::EdgeDescriptor ED;
1490 ED e= graph.findEdge(edge.target(), edge.source());
1491 if(e!=std::numeric_limits<ED>::max())
1492 {
1493 typename G::EdgeProperties& other = graph.getEdgeProperties(e);
1494 other.setInfluences();
1495 }
1496 }
1497 }
1498
1499 template<class M, class N>
1500 inline bool Dependency<M,N>::isIsolated()
1501 {
1502 return maxValue_ < beta() * diagonal_;
1503 }
1504
1505 template<class G,class S>
1506 Aggregate<G,S>::Aggregate(MatrixGraph& graph, AggregatesMap<Vertex>& aggregates,
1507 VertexSet& connected, std::vector<Vertex>& front)
1508 : vertices_(), id_(-1), graph_(graph), aggregates_(aggregates),
1509 connected_(connected), front_(front)
1510 {}
1511
1512 template<class G,class S>
1513 void Aggregate<G,S>::reconstruct(const Vertex& vertex)
1514 {
1515 /*
1516 vertices_.push_back(vertex);
1517 typedef typename VertexList::const_iterator iterator;
1518 iterator begin = vertices_.begin();
1519 iterator end = vertices_.end();*/
1520 throw "Not yet implemented";
1521
1522 // while(begin!=end){
1523 //for();
1524 // }
1525
1526 }
1527
1528 template<class G,class S>
1529 inline void Aggregate<G,S>::seed(const Vertex& vertex)
1530 {
1531 dvverb<<"Connected cleared"<<std::endl;
1532 connected_.clear();
1533 vertices_.clear();
1534 connected_.insert(vertex);
1535 dvverb << " Inserting "<<vertex<<" size="<<connected_.size();
1536 ++id_ ;
1537 add(vertex);
1538 }
1539
1540
1541 template<class G,class S>
1542 inline void Aggregate<G,S>::add(const Vertex& vertex)
1543 {
1544 vertices_.insert(vertex);
1545 aggregates_[vertex]=id_;
1546 if(front_.size())
1547 front_.erase(std::lower_bound(front_.begin(), front_.end(), vertex));
1548
1549
1550 typedef typename MatrixGraph::ConstEdgeIterator iterator;
1551 const iterator end = graph_.endEdges(vertex);
1552 for(iterator edge = graph_.beginEdges(vertex); edge != end; ++edge) {
1553 dvverb << " Inserting "<<aggregates_[edge.target()];
1554 connected_.insert(aggregates_[edge.target()]);
1555 dvverb <<" size="<<connected_.size();
1556 if(aggregates_[edge.target()]==AggregatesMap<Vertex>::UNAGGREGATED &&
1557 !graph_.getVertexProperties(edge.target()).front())
1558 {
1559 front_.push_back(edge.target());
1560 graph_.getVertexProperties(edge.target()).setFront();
1561 }
1562 }
1563 dvverb <<std::endl;
1564 std::sort(front_.begin(), front_.end());
1565 }
1566
1567 template<class G,class S>
1568 inline void Aggregate<G,S>::add(std::vector<Vertex>& vertices)
1569 {
1570#ifndef NDEBUG
1571 std::size_t oldsize = vertices_.size();
1572#endif
1573 typedef typename std::vector<Vertex>::iterator Iterator;
1574
1575 typedef typename VertexSet::iterator SIterator;
1576
1577 SIterator pos=vertices_.begin();
1578 std::vector<Vertex> newFront;
1579 newFront.reserve(front_.capacity());
1580
1581 std::set_difference(front_.begin(), front_.end(), vertices.begin(), vertices.end(),
1582 std::back_inserter(newFront));
1583 front_=newFront;
1584
1585 for(Iterator vertex=vertices.begin(); vertex != vertices.end(); ++vertex)
1586 {
1587 pos=vertices_.insert(pos,*vertex);
1588 vertices_.insert(*vertex);
1589 graph_.getVertexProperties(*vertex).resetFront(); // Not a front node any more.
1590 aggregates_[*vertex]=id_;
1591
1592 typedef typename MatrixGraph::ConstEdgeIterator iterator;
1593 const iterator end = graph_.endEdges(*vertex);
1594 for(iterator edge = graph_.beginEdges(*vertex); edge != end; ++edge) {
1595 dvverb << " Inserting "<<aggregates_[edge.target()];
1596 connected_.insert(aggregates_[edge.target()]);
1597 if(aggregates_[edge.target()]==AggregatesMap<Vertex>::UNAGGREGATED &&
1598 !graph_.getVertexProperties(edge.target()).front())
1599 {
1600 front_.push_back(edge.target());
1601 graph_.getVertexProperties(edge.target()).setFront();
1602 }
1603 dvverb <<" size="<<connected_.size();
1604 }
1605 dvverb <<std::endl;
1606 }
1607 std::sort(front_.begin(), front_.end());
1608 assert(oldsize+vertices.size()==vertices_.size());
1609 }
1610 template<class G,class S>
1611 inline void Aggregate<G,S>::clear()
1612 {
1613 vertices_.clear();
1614 connected_.clear();
1615 id_=-1;
1616 }
1617
1618 template<class G,class S>
1619 inline typename Aggregate<G,S>::VertexSet::size_type
1621 {
1622 return vertices_.size();
1623 }
1624
1625 template<class G,class S>
1626 inline typename Aggregate<G,S>::VertexSet::size_type
1628 {
1629 return connected_.size();
1630 }
1631
1632 template<class G,class S>
1633 inline int Aggregate<G,S>::id()
1634 {
1635 return id_;
1636 }
1637
1638 template<class G,class S>
1640 {
1641 return vertices_.begin();
1642 }
1643
1644 template<class G,class S>
1646 {
1647 return vertices_.end();
1648 }
1649
1650 template<class V>
1651 const V AggregatesMap<V>::UNAGGREGATED = std::numeric_limits<V>::max();
1652
1653 template<class V>
1654 const V AggregatesMap<V>::ISOLATED = std::numeric_limits<V>::max()-1;
1655
1656 template<class V>
1658 : aggregates_(0)
1659 {}
1660
1661 template<class V>
1663 {
1664 if(aggregates_!=0)
1665 delete[] aggregates_;
1666 }
1667
1668
1669 template<class V>
1670 inline AggregatesMap<V>::AggregatesMap(std::size_t noVertices)
1671 {
1672 allocate(noVertices);
1673 }
1674
1675 template<class V>
1676 inline std::size_t AggregatesMap<V>::noVertices() const
1677 {
1678 return noVertices_;
1679 }
1680
1681 template<class V>
1682 inline void AggregatesMap<V>::allocate(std::size_t noVertices)
1683 {
1684 aggregates_ = new AggregateDescriptor[noVertices];
1685 noVertices_ = noVertices;
1686
1687 for(std::size_t i=0; i < noVertices; i++)
1688 aggregates_[i]=UNAGGREGATED;
1689 }
1690
1691 template<class V>
1692 inline void AggregatesMap<V>::free()
1693 {
1694 assert(aggregates_ != 0);
1695 delete[] aggregates_;
1696 aggregates_=0;
1697 }
1698
1699 template<class V>
1701 AggregatesMap<V>::operator[](const VertexDescriptor& v)
1702 {
1703 return aggregates_[v];
1704 }
1705
1706 template<class V>
1707 inline const typename AggregatesMap<V>::AggregateDescriptor&
1708 AggregatesMap<V>::operator[](const VertexDescriptor& v) const
1709 {
1710 return aggregates_[v];
1711 }
1712
1713 template<class V>
1714 template<bool reset, class G, class F,class VM>
1715 inline std::size_t AggregatesMap<V>::breadthFirstSearch(const V& start,
1716 const AggregateDescriptor& aggregate,
1717 const G& graph, F& aggregateVisitor,
1718 VM& visitedMap) const
1719 {
1720 VertexList vlist;
1721
1722 DummyEdgeVisitor dummy;
1723 return breadthFirstSearch<true,reset>(start, aggregate, graph, vlist, aggregateVisitor, dummy, visitedMap);
1724 }
1725
1726 template<class V>
1727 template<bool remove, bool reset, class G, class L, class F1, class F2, class VM>
1728 std::size_t AggregatesMap<V>::breadthFirstSearch(const V& start,
1729 const AggregateDescriptor& aggregate,
1730 const G& graph,
1731 L& visited,
1732 F1& aggregateVisitor,
1733 F2& nonAggregateVisitor,
1734 VM& visitedMap) const
1735 {
1736 typedef typename L::const_iterator ListIterator;
1737 int visitedSpheres = 0;
1738
1739 visited.push_back(start);
1740 put(visitedMap, start, true);
1741
1742 ListIterator current = visited.begin();
1743 ListIterator end = visited.end();
1744 std::size_t i=0, size=visited.size();
1745
1746 // visit the neighbours of all vertices of the
1747 // current sphere.
1748 while(current != end) {
1749
1750 for(; i<size; ++current, ++i) {
1751 typedef typename G::ConstEdgeIterator EdgeIterator;
1752 const EdgeIterator endEdge = graph.endEdges(*current);
1753
1754 for(EdgeIterator edge = graph.beginEdges(*current);
1755 edge != endEdge; ++edge) {
1756
1757 if(aggregates_[edge.target()]==aggregate) {
1758 if(!get(visitedMap, edge.target())) {
1759 put(visitedMap, edge.target(), true);
1760 visited.push_back(edge.target());
1761 aggregateVisitor(edge);
1762 }
1763 }else
1764 nonAggregateVisitor(edge);
1765 }
1766 }
1767 end = visited.end();
1768 size = visited.size();
1769 if(current != end)
1770 visitedSpheres++;
1771 }
1772
1773 if(reset)
1774 for(current = visited.begin(); current != end; ++current)
1775 put(visitedMap, *current, false);
1776
1777
1778 if(remove)
1779 visited.clear();
1780
1781 return visitedSpheres;
1782 }
1783
1784 template<class G>
1786 : graph_(0), aggregate_(0), front_(), connected_(), size_(-1)
1787 {}
1788
1789 template<class G>
1791 {
1792 size_=-1;
1793 }
1794
1795 template<class G, class C>
1796 void buildDependency(G& graph,
1797 const typename C::Matrix& matrix,
1798 C criterion, bool firstlevel)
1799 {
1800 // assert(graph.isBuilt());
1801 typedef typename C::Matrix Matrix;
1802 typedef typename G::VertexIterator VertexIterator;
1803
1804 criterion.init(&matrix);
1805
1806 for(VertexIterator vertex = graph.begin(); vertex != graph.end(); ++vertex) {
1807 typedef typename Matrix::row_type Row;
1808
1809 const Row& row = matrix[*vertex];
1810
1811 // Tell the criterion what row we will examine now
1812 // This might for example be used for calculating the
1813 // maximum offdiagonal value
1814 criterion.initRow(row, *vertex);
1815
1816 // On a first path all columns are examined. After this
1817 // the calculator should know whether the vertex is isolated.
1818 typedef typename Matrix::ConstColIterator ColIterator;
1819 ColIterator end = row.end();
1820 typename FieldTraits<typename Matrix::field_type>::real_type absoffdiag=0.;
1821
1822 if(firstlevel) {
1823 for(ColIterator col = row.begin(); col != end; ++col)
1824 if(col.index()!=*vertex) {
1825 criterion.examine(col);
1826 absoffdiag=std::max(absoffdiag, col->frobenius_norm());
1827 }
1828
1829 if(absoffdiag==0)
1830 vertex.properties().setExcludedBorder();
1831 }
1832 else
1833 for(ColIterator col = row.begin(); col != end; ++col)
1834 if(col.index()!=*vertex)
1835 criterion.examine(col);
1836
1837 // reset the vertex properties
1838 //vertex.properties().reset();
1839
1840 // Check whether the vertex is isolated.
1841 if(criterion.isIsolated()) {
1842 //std::cout<<"ISOLATED: "<<*vertex<<std::endl;
1843 vertex.properties().setIsolated();
1844 }else{
1845 // Examine all the edges beginning at this vertex.
1846 typedef typename G::EdgeIterator EdgeIterator;
1847 typedef typename Matrix::ConstColIterator ColIterator;
1848 EdgeIterator eEnd = vertex.end();
1849 ColIterator col = matrix[*vertex].begin();
1850
1851 for(EdgeIterator edge = vertex.begin(); edge!= eEnd; ++edge, ++col) {
1852 // Move to the right column.
1853 while(col.index()!=edge.target())
1854 ++col;
1855 criterion.examine(graph, edge, col);
1856 }
1857 }
1858
1859 }
1860 }
1861
1862
1863 template<class G>
1864 template<class V>
1865 inline Aggregator<G>::AggregateVisitor<V>::AggregateVisitor(const AggregatesMap<Vertex>& aggregates,
1866 const AggregateDescriptor& aggregate, V& visitor)
1867 : aggregates_(aggregates), aggregate_(aggregate), visitor_(&visitor)
1868 {}
1869
1870 template<class G>
1871 template<class V>
1872 inline void Aggregator<G>::AggregateVisitor<V>::operator()(const typename MatrixGraph::ConstEdgeIterator& edge)
1873 {
1874 if(aggregates_[edge.target()]==aggregate_)
1875 visitor_->operator()(edge);
1876 }
1877
1878 template<class G>
1879 template<class V>
1880 inline void Aggregator<G>::visitAggregateNeighbours(const Vertex& vertex,
1881 const AggregateDescriptor& aggregate,
1882 const AggregatesMap<Vertex>& aggregates,
1883 V& visitor) const
1884 {
1885 // Only evaluates for edge pointing to the aggregate
1886 AggregateVisitor<V> v(aggregates, aggregate, visitor);
1887 visitNeighbours(*graph_, vertex, v);
1888 }
1889
1890
1891 template<class G>
1892 inline Aggregator<G>::Counter::Counter()
1893 : count_(0)
1894 {}
1895
1896 template<class G>
1897 inline void Aggregator<G>::Counter::increment()
1898 {
1899 ++count_;
1900 }
1901
1902 template<class G>
1903 inline void Aggregator<G>::Counter::decrement()
1904 {
1905 --count_;
1906 }
1907 template<class G>
1908 inline int Aggregator<G>::Counter::value()
1909 {
1910 return count_;
1911 }
1912
1913 template<class G>
1914 inline void Aggregator<G>::TwoWayCounter::operator()(const typename MatrixGraph::ConstEdgeIterator& edge)
1915 {
1916 if(edge.properties().isTwoWay())
1917 Counter::increment();
1918 }
1919
1920 template<class G>
1921 int Aggregator<G>::twoWayConnections(const Vertex& vertex, const AggregateDescriptor& aggregate,
1922 const AggregatesMap<Vertex>& aggregates) const
1923 {
1924 TwoWayCounter counter;
1925 visitAggregateNeighbours(vertex, aggregate, aggregates, counter);
1926 return counter.value();
1927 }
1928
1929 template<class G>
1930 int Aggregator<G>::oneWayConnections(const Vertex& vertex, const AggregateDescriptor& aggregate,
1931 const AggregatesMap<Vertex>& aggregates) const
1932 {
1933 OneWayCounter counter;
1934 visitAggregateNeighbours(vertex, aggregate, aggregates, counter);
1935 return counter.value();
1936 }
1937
1938 template<class G>
1939 inline void Aggregator<G>::OneWayCounter::operator()(const typename MatrixGraph::ConstEdgeIterator& edge)
1940 {
1941 if(edge.properties().isOneWay())
1942 Counter::increment();
1943 }
1944
1945 template<class G>
1946 inline Aggregator<G>::ConnectivityCounter::ConnectivityCounter(const VertexSet& connected,
1947 const AggregatesMap<Vertex>& aggregates)
1948 : Counter(), connected_(connected), aggregates_(aggregates)
1949 {}
1950
1951
1952 template<class G>
1953 inline void Aggregator<G>::ConnectivityCounter::operator()(const typename MatrixGraph::ConstEdgeIterator& edge)
1954 {
1955 if(connected_.find(aggregates_[edge.target()]) == connected_.end() || aggregates_[edge.target()]==AggregatesMap<Vertex>::UNAGGREGATED)
1956 // Would be a new connection
1957 Counter::increment();
1958 else{
1959 Counter::increment();
1960 Counter::increment();
1961 }
1962 }
1963
1964 template<class G>
1965 inline double Aggregator<G>::connectivity(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates) const
1966 {
1967 ConnectivityCounter counter(connected_, aggregates);
1968 double noNeighbours=visitNeighbours(*graph_, vertex, counter);
1969 return (double)counter.value()/noNeighbours;
1970 }
1971
1972 template<class G>
1973 inline Aggregator<G>::DependencyCounter::DependencyCounter()
1974 : Counter()
1975 {}
1976
1977 template<class G>
1978 inline void Aggregator<G>::DependencyCounter::operator()(const typename MatrixGraph::ConstEdgeIterator& edge)
1979 {
1980 if(edge.properties().depends())
1981 Counter::increment();
1982 if(edge.properties().influences())
1983 Counter::increment();
1984 }
1985
1986 template<class G>
1987 int Aggregator<G>::unusedNeighbours(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates) const
1988 {
1989 return aggregateNeighbours(vertex, AggregatesMap<Vertex>::UNAGGREGATED, aggregates);
1990 }
1991
1992 template<class G>
1993 std::pair<int,int> Aggregator<G>::neighbours(const Vertex& vertex,
1994 const AggregateDescriptor& aggregate,
1995 const AggregatesMap<Vertex>& aggregates) const
1996 {
1997 DependencyCounter unused, aggregated;
1998 typedef AggregateVisitor<DependencyCounter> Counter;
1999 typedef std::tuple<Counter,Counter> CounterTuple;
2000 CombinedFunctor<CounterTuple> visitors(CounterTuple(Counter(aggregates, AggregatesMap<Vertex>::UNAGGREGATED, unused), Counter(aggregates, aggregate, aggregated)));
2001 visitNeighbours(*graph_, vertex, visitors);
2002 return std::make_pair(unused.value(), aggregated.value());
2003 }
2004
2005
2006 template<class G>
2007 int Aggregator<G>::aggregateNeighbours(const Vertex& vertex, const AggregateDescriptor& aggregate, const AggregatesMap<Vertex>& aggregates) const
2008 {
2009 DependencyCounter counter;
2010 visitAggregateNeighbours(vertex, aggregate, aggregates, counter);
2011 return counter.value();
2012 }
2013
2014 template<class G>
2015 std::size_t Aggregator<G>::distance(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates)
2016 {
2017 return 0;
2018 typename PropertyMapTypeSelector<VertexVisitedTag,G>::Type visitedMap = get(VertexVisitedTag(), *graph_);
2019 VertexList vlist;
2020 typename AggregatesMap<Vertex>::DummyEdgeVisitor dummy;
2021 return aggregates.template breadthFirstSearch<true,true>(vertex,
2022 aggregate_->id(), *graph_,
2023 vlist, dummy, dummy, visitedMap);
2024 }
2025
2026 template<class G>
2027 inline Aggregator<G>::FrontMarker::FrontMarker(std::vector<Vertex>& front, MatrixGraph& graph)
2028 : front_(front), graph_(graph)
2029 {}
2030
2031 template<class G>
2032 inline void Aggregator<G>::FrontMarker::operator()(const typename MatrixGraph::ConstEdgeIterator& edge)
2033 {
2034 Vertex target = edge.target();
2035
2036 if(!graph_.getVertexProperties(target).front()) {
2037 front_.push_back(target);
2038 graph_.getVertexProperties(target).setFront();
2039 }
2040 }
2041
2042 template<class G>
2043 inline bool Aggregator<G>::admissible(const Vertex& vertex, const AggregateDescriptor& aggregate, const AggregatesMap<Vertex>& aggregates) const
2044 {
2045 // Todo
2046 Dune::dvverb<<" Admissible not yet implemented!"<<std::endl;
2047 return true;
2048 //Situation 1: front node depends on two nodes. Then these
2049 // have to be strongly connected to each other
2050
2051 // Iterate over all all neighbours of front node
2052 typedef typename MatrixGraph::ConstEdgeIterator Iterator;
2053 Iterator vend = graph_->endEdges(vertex);
2054 for(Iterator edge = graph_->beginEdges(vertex); edge != vend; ++edge) {
2055 // if(edge.properties().depends() && !edge.properties().influences()
2056 if(edge.properties().isStrong()
2057 && aggregates[edge.target()]==aggregate)
2058 {
2059 // Search for another link to the aggregate
2060 Iterator edge1 = edge;
2061 for(++edge1; edge1 != vend; ++edge1) {
2062 //if(edge1.properties().depends() && !edge1.properties().influences()
2063 if(edge1.properties().isStrong()
2064 && aggregates[edge.target()]==aggregate)
2065 {
2066 //Search for an edge connecting the two vertices that is
2067 //strong
2068 bool found=false;
2069 Iterator v2end = graph_->endEdges(edge.target());
2070 for(Iterator edge2 = graph_->beginEdges(edge.target()); edge2 != v2end; ++edge2) {
2071 if(edge2.target()==edge1.target() &&
2072 edge2.properties().isStrong()) {
2073 found =true;
2074 break;
2075 }
2076 }
2077 if(found)
2078 {
2079 return true;
2080 }
2081 }
2082 }
2083 }
2084 }
2085
2086 // Situation 2: cluster node depends on front node and other cluster node
2088 vend = graph_->endEdges(vertex);
2089 for(Iterator edge = graph_->beginEdges(vertex); edge != vend; ++edge) {
2090 //if(!edge.properties().depends() && edge.properties().influences()
2091 if(edge.properties().isStrong()
2092 && aggregates[edge.target()]==aggregate)
2093 {
2094 // Search for a link from target that stays within the aggregate
2095 Iterator v1end = graph_->endEdges(edge.target());
2096
2097 for(Iterator edge1=graph_->beginEdges(edge.target()); edge1 != v1end; ++edge1) {
2098 //if(edge1.properties().depends() && !edge1.properties().influences()
2099 if(edge1.properties().isStrong()
2100 && aggregates[edge1.target()]==aggregate)
2101 {
2102 bool found=false;
2103 // Check if front node is also connected to this one
2104 Iterator v2end = graph_->endEdges(vertex);
2105 for(Iterator edge2 = graph_->beginEdges(vertex); edge2 != v2end; ++edge2) {
2106 if(edge2.target()==edge1.target()) {
2107 if(edge2.properties().isStrong())
2108 found=true;
2109 break;
2110 }
2111 }
2112 if(found)
2113 {
2114 return true;
2115 }
2116 }
2117 }
2118 }
2119 }
2120 return false;
2121 }
2122
2123 template<class G>
2124 void Aggregator<G>::unmarkFront()
2125 {
2126 typedef typename std::vector<Vertex>::const_iterator Iterator;
2127
2128 for(Iterator vertex=front_.begin(); vertex != front_.end(); ++vertex)
2129 graph_->getVertexProperties(*vertex).resetFront();
2130
2131 front_.clear();
2132 }
2133
2134 template<class G>
2135 inline void
2136 Aggregator<G>::nonisoNeighbourAggregate(const Vertex& vertex,
2137 const AggregatesMap<Vertex>& aggregates,
2138 SLList<Vertex>& neighbours) const
2139 {
2140 typedef typename MatrixGraph::ConstEdgeIterator Iterator;
2141 Iterator end=graph_->beginEdges(vertex);
2142 neighbours.clear();
2143
2144 for(Iterator edge=graph_->beginEdges(vertex); edge!=end; ++edge)
2145 {
2146 if(aggregates[edge.target()]!=AggregatesMap<Vertex>::UNAGGREGATED && graph_->getVertexProperties(edge.target()).isolated())
2147 neighbours.push_back(aggregates[edge.target()]);
2148 }
2149 }
2150
2151 template<class G>
2152 inline typename G::VertexDescriptor Aggregator<G>::mergeNeighbour(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates) const
2153 {
2154 typedef typename MatrixGraph::ConstEdgeIterator Iterator;
2155
2156 Iterator end = graph_->endEdges(vertex);
2157 for(Iterator edge = graph_->beginEdges(vertex); edge != end; ++edge) {
2158 if(aggregates[edge.target()] != AggregatesMap<Vertex>::UNAGGREGATED &&
2159 graph_->getVertexProperties(edge.target()).isolated() == graph_->getVertexProperties(edge.source()).isolated()) {
2160 if( graph_->getVertexProperties(vertex).isolated() ||
2161 ((edge.properties().depends() || edge.properties().influences())
2162 && admissible(vertex, aggregates[edge.target()], aggregates)))
2163 return edge.target();
2164 }
2165 }
2167 }
2168
2169 template<class G>
2170 Aggregator<G>::FrontNeighbourCounter::FrontNeighbourCounter(const MatrixGraph& graph)
2171 : Counter(), graph_(graph)
2172 {}
2173
2174 template<class G>
2175 void Aggregator<G>::FrontNeighbourCounter::operator()(const typename MatrixGraph::ConstEdgeIterator& edge)
2176 {
2177 if(graph_.getVertexProperties(edge.target()).front())
2178 Counter::increment();
2179 }
2180
2181 template<class G>
2182 int Aggregator<G>::noFrontNeighbours(const Vertex& vertex) const
2183 {
2184 FrontNeighbourCounter counter(*graph_);
2185 visitNeighbours(*graph_, vertex, counter);
2186 return counter.value();
2187 }
2188 template<class G>
2189 inline bool Aggregator<G>::connected(const Vertex& vertex,
2190 const AggregateDescriptor& aggregate,
2191 const AggregatesMap<Vertex>& aggregates) const
2192 {
2193 typedef typename G::ConstEdgeIterator iterator;
2194 const iterator end = graph_->endEdges(vertex);
2195 for(iterator edge = graph_->beginEdges(vertex); edge != end; ++edge)
2196 if(aggregates[edge.target()]==aggregate)
2197 return true;
2198 return false;
2199 }
2200 template<class G>
2201 inline bool Aggregator<G>::connected(const Vertex& vertex,
2202 const SLList<AggregateDescriptor>& aggregateList,
2203 const AggregatesMap<Vertex>& aggregates) const
2204 {
2205 typedef typename SLList<AggregateDescriptor>::const_iterator Iter;
2206 for(Iter i=aggregateList.begin(); i!=aggregateList.end(); ++i)
2207 if(connected(vertex, *i, aggregates))
2208 return true;
2209 return false;
2210 }
2211
2212 template<class G>
2213 template<class C>
2214 void Aggregator<G>::growIsolatedAggregate(const Vertex& seed, const AggregatesMap<Vertex>& aggregates, const C& c)
2215 {
2216 SLList<Vertex> connectedAggregates;
2217 nonisoNeighbourAggregate(seed, aggregates,connectedAggregates);
2218
2219 while(aggregate_->size()< c.minAggregateSize() && aggregate_->connectSize() < c.maxConnectivity()) {
2220 double maxCon=-1;
2221 std::size_t maxFrontNeighbours=0;
2222
2224
2225 typedef typename std::vector<Vertex>::const_iterator Iterator;
2226
2227 for(Iterator vertex = front_.begin(); vertex != front_.end(); ++vertex) {
2228 if(distance(*vertex, aggregates)>c.maxDistance())
2229 continue; // distance of proposes aggregate too big
2230
2231 if(connectedAggregates.size()>0) {
2232 // there is already a neighbour cluster
2233 // front node must be connected to same neighbour cluster
2234
2235 if(!connected(*vertex, connectedAggregates, aggregates))
2236 continue;
2237 }
2238
2239 double con = connectivity(*vertex, aggregates);
2240
2241 if(con == maxCon) {
2242 std::size_t frontNeighbours = noFrontNeighbours(*vertex);
2243
2244 if(frontNeighbours >= maxFrontNeighbours) {
2245 maxFrontNeighbours = frontNeighbours;
2246 candidate = *vertex;
2247 }
2248 }else if(con > maxCon) {
2249 maxCon = con;
2250 maxFrontNeighbours = noFrontNeighbours(*vertex);
2251 candidate = *vertex;
2252 }
2253 }
2254
2256 break;
2257
2258 aggregate_->add(candidate);
2259 }
2260 }
2261
2262 template<class G>
2263 template<class C>
2264 void Aggregator<G>::growAggregate(const Vertex& seed, const AggregatesMap<Vertex>& aggregates, const C& c)
2265 {
2266
2267 std::size_t distance_ =0;
2268 while(aggregate_->size() < c.minAggregateSize()&& distance_<c.maxDistance()) {
2269 int maxTwoCons=0, maxOneCons=0, maxNeighbours=-1;
2270 double maxCon=-1;
2271
2272 std::vector<Vertex> candidates;
2273 candidates.reserve(30);
2274
2275 typedef typename std::vector<Vertex>::const_iterator Iterator;
2276
2277 for(Iterator vertex = front_.begin(); vertex != front_.end(); ++vertex) {
2278 // Only nonisolated nodes are considered
2279 if(graph_->getVertexProperties(*vertex).isolated())
2280 continue;
2281
2282 int twoWayCons = twoWayConnections(*vertex, aggregate_->id(), aggregates);
2283
2284 /* The case of two way connections. */
2285 if( maxTwoCons == twoWayCons && twoWayCons > 0) {
2286 double con = connectivity(*vertex, aggregates);
2287
2288 if(con == maxCon) {
2289 int neighbours = noFrontNeighbours(*vertex);
2290
2291 if(neighbours > maxNeighbours) {
2292 maxNeighbours = neighbours;
2293 candidates.clear();
2294 candidates.push_back(*vertex);
2295 }else{
2296 candidates.push_back(*vertex);
2297 }
2298 }else if( con > maxCon) {
2299 maxCon = con;
2300 maxNeighbours = noFrontNeighbours(*vertex);
2301 candidates.clear();
2302 candidates.push_back(*vertex);
2303 }
2304 }else if(twoWayCons > maxTwoCons) {
2305 maxTwoCons = twoWayCons;
2306 maxCon = connectivity(*vertex, aggregates);
2307 maxNeighbours = noFrontNeighbours(*vertex);
2308 candidates.clear();
2309 candidates.push_back(*vertex);
2310
2311 // two way connections precede
2312 maxOneCons = std::numeric_limits<int>::max();
2313 }
2314
2315 if(twoWayCons > 0)
2316 {
2317 continue; // THis is a two-way node, skip tests for one way nodes
2318 }
2319
2320 /* The one way case */
2321 int oneWayCons = oneWayConnections(*vertex, aggregate_->id(), aggregates);
2322
2323 if(oneWayCons==0)
2324 continue; // No strong connections, skip the tests.
2325
2326 if(!admissible(*vertex, aggregate_->id(), aggregates))
2327 continue;
2328
2329 if( maxOneCons == oneWayCons && oneWayCons > 0) {
2330 double con = connectivity(*vertex, aggregates);
2331
2332 if(con == maxCon) {
2333 int neighbours = noFrontNeighbours(*vertex);
2334
2335 if(neighbours > maxNeighbours) {
2336 maxNeighbours = neighbours;
2337 candidates.clear();
2338 candidates.push_back(*vertex);
2339 }else{
2340 if(neighbours==maxNeighbours)
2341 {
2342 candidates.push_back(*vertex);
2343 }
2344 }
2345 }else if( con > maxCon) {
2346 maxCon = con;
2347 maxNeighbours = noFrontNeighbours(*vertex);
2348 candidates.clear();
2349 candidates.push_back(*vertex);
2350 }
2351 }else if(oneWayCons > maxOneCons) {
2352 maxOneCons = oneWayCons;
2353 maxCon = connectivity(*vertex, aggregates);
2354 maxNeighbours = noFrontNeighbours(*vertex);
2355 candidates.clear();
2356 candidates.push_back(*vertex);
2357 }
2358 }
2359
2360
2361 if(!candidates.size())
2362 break; // No more candidates found
2363 distance_=distance(seed, aggregates);
2364 candidates.resize(std::min(candidates.size(), c.maxAggregateSize()-
2365 aggregate_->size()));
2366 aggregate_->add(candidates);
2367 }
2368 }
2369
2370 template<typename V>
2371 template<typename M, typename G, typename C>
2372 std::tuple<int,int,int,int> AggregatesMap<V>::buildAggregates(const M& matrix, G& graph, const C& criterion,
2373 bool finestLevel)
2374 {
2375 Aggregator<G> aggregator;
2376 return aggregator.build(matrix, graph, *this, criterion, finestLevel);
2377 }
2378
2379 template<class G>
2380 template<class M, class C>
2381 std::tuple<int,int,int,int> Aggregator<G>::build(const M& m, G& graph, AggregatesMap<Vertex>& aggregates, const C& c,
2382 bool finestLevel)
2383 {
2384 // Stack for fast vertex access
2385 Stack stack_(graph, *this, aggregates);
2386
2387 graph_ = &graph;
2388
2389 aggregate_ = new Aggregate<G,VertexSet>(graph, aggregates, connected_, front_);
2390
2391 Timer watch;
2392 watch.reset();
2393
2394 buildDependency(graph, m, c, finestLevel);
2395
2396 dverb<<"Build dependency took "<< watch.elapsed()<<" seconds."<<std::endl;
2397 int noAggregates, conAggregates, isoAggregates, oneAggregates;
2398 std::size_t maxA=0, minA=1000000, avg=0;
2399 int skippedAggregates;
2400 noAggregates = conAggregates = isoAggregates = oneAggregates =
2401 skippedAggregates = 0;
2402
2403 while(true) {
2404 Vertex seed = stack_.pop();
2405
2406 if(seed == Stack::NullEntry)
2407 // No more unaggregated vertices. We are finished!
2408 break;
2409
2410 // Debugging output
2411 if((noAggregates+1)%10000 == 0)
2412 Dune::dverb<<"c";
2413 unmarkFront();
2414
2415 if(graph.getVertexProperties(seed).excludedBorder()) {
2416 aggregates[seed]=AggregatesMap<Vertex>::ISOLATED;
2417 ++skippedAggregates;
2418 continue;
2419 }
2420
2421 if(graph.getVertexProperties(seed).isolated()) {
2422 if(c.skipIsolated()) {
2423 // isolated vertices are not aggregated but skipped on the coarser levels.
2424 aggregates[seed]=AggregatesMap<Vertex>::ISOLATED;
2425 ++skippedAggregates;
2426 // skip rest as no agglomeration is done.
2427 continue;
2428 }else{
2429 aggregate_->seed(seed);
2430 growIsolatedAggregate(seed, aggregates, c);
2431 }
2432 }else{
2433 aggregate_->seed(seed);
2434 growAggregate(seed, aggregates, c);
2435 }
2436
2437 /* The rounding step. */
2438 while(!(graph.getVertexProperties(seed).isolated()) && aggregate_->size() < c.maxAggregateSize()) {
2439
2440 std::vector<Vertex> candidates;
2441 candidates.reserve(30);
2442
2443 typedef typename std::vector<Vertex>::const_iterator Iterator;
2444
2445 for(Iterator vertex = front_.begin(); vertex != front_.end(); ++vertex) {
2446
2447 if(graph.getVertexProperties(*vertex).isolated())
2448 continue; // No isolated nodes here
2449
2450 if(twoWayConnections( *vertex, aggregate_->id(), aggregates) == 0 &&
2451 (oneWayConnections( *vertex, aggregate_->id(), aggregates) == 0 ||
2452 !admissible( *vertex, aggregate_->id(), aggregates) ))
2453 continue;
2454
2455 std::pair<int,int> neighbourPair=neighbours(*vertex, aggregate_->id(),
2456 aggregates);
2457
2458 //if(aggregateNeighbours(*vertex, aggregate_->id(), aggregates) <= unusedNeighbours(*vertex, aggregates))
2459 // continue;
2460
2461 if(neighbourPair.first >= neighbourPair.second)
2462 continue;
2463
2464 if(distance(*vertex, aggregates) > c.maxDistance())
2465 continue; // Distance too far
2466 candidates.push_back(*vertex);
2467 break;
2468 }
2469
2470 if(!candidates.size()) break; // no more candidates found.
2471
2472 candidates.resize(std::min(candidates.size(), c.maxAggregateSize()-
2473 aggregate_->size()));
2474 aggregate_->add(candidates);
2475
2476 }
2477
2478 // try to merge aggregates consisting of only one nonisolated vertex with other aggregates
2479 if(aggregate_->size()==1 && c.maxAggregateSize()>1) {
2480 if(!graph.getVertexProperties(seed).isolated()) {
2481 Vertex mergedNeighbour = mergeNeighbour(seed, aggregates);
2482
2483 if(mergedNeighbour != AggregatesMap<Vertex>::UNAGGREGATED) {
2484 // assign vertex to the neighbouring cluster
2485 aggregates[seed] = aggregates[mergedNeighbour];
2486 aggregate_->invalidate();
2487 }else{
2488 ++avg;
2489 minA=std::min(minA,static_cast<std::size_t>(1));
2490 maxA=std::max(maxA,static_cast<std::size_t>(1));
2491 ++oneAggregates;
2492 ++conAggregates;
2493 }
2494 }else{
2495 ++avg;
2496 minA=std::min(minA,static_cast<std::size_t>(1));
2497 maxA=std::max(maxA,static_cast<std::size_t>(1));
2498 ++oneAggregates;
2499 ++isoAggregates;
2500 }
2501 ++avg;
2502 }else{
2503 avg+=aggregate_->size();
2504 minA=std::min(minA,aggregate_->size());
2505 maxA=std::max(maxA,aggregate_->size());
2506 if(graph.getVertexProperties(seed).isolated())
2507 ++isoAggregates;
2508 else
2509 ++conAggregates;
2510 }
2511
2512 }
2513
2514 Dune::dinfo<<"connected aggregates: "<<conAggregates;
2515 Dune::dinfo<<" isolated aggregates: "<<isoAggregates;
2516 if(conAggregates+isoAggregates>0)
2517 Dune::dinfo<<" one node aggregates: "<<oneAggregates<<" min size="
2518 <<minA<<" max size="<<maxA
2519 <<" avg="<<avg/(conAggregates+isoAggregates)<<std::endl;
2520
2521 delete aggregate_;
2522 return std::make_tuple(conAggregates+isoAggregates,isoAggregates,
2523 oneAggregates,skippedAggregates);
2524 }
2525
2526
2527 template<class G>
2528 Aggregator<G>::Stack::Stack(const MatrixGraph& graph, const Aggregator<G>& aggregatesBuilder,
2529 const AggregatesMap<Vertex>& aggregates)
2530 : graph_(graph), aggregatesBuilder_(aggregatesBuilder), aggregates_(aggregates), begin_(graph.begin()), end_(graph.end())
2531 {
2532 //vals_ = new Vertex[N];
2533 }
2534
2535 template<class G>
2536 Aggregator<G>::Stack::~Stack()
2537 {
2538 //Dune::dverb << "Max stack size was "<<maxSize_<<" filled="<<filled_<<std::endl;
2539 //delete[] vals_;
2540 }
2541
2542 template<class G>
2543 const typename Aggregator<G>::Vertex Aggregator<G>::Stack::NullEntry
2544 = std::numeric_limits<typename G::VertexDescriptor>::max();
2545
2546 template<class G>
2547 inline typename G::VertexDescriptor Aggregator<G>::Stack::pop()
2548 {
2549 for(; begin_!=end_ && aggregates_[*begin_] != AggregatesMap<Vertex>::UNAGGREGATED; ++begin_) ;
2550
2551 if(begin_!=end_)
2552 {
2553 typename G::VertexDescriptor current=*begin_;
2554 ++begin_;
2555 return current;
2556 }else
2557 return NullEntry;
2558 }
2559
2560#endif // DOXYGEN
2561
2562 template<class V>
2563 void printAggregates2d(const AggregatesMap<V>& aggregates, int n, int m, std::ostream& os)
2564 {
2565 std::ios_base::fmtflags oldOpts=os.flags();
2566
2567 os.setf(std::ios_base::right, std::ios_base::adjustfield);
2568
2569 V max=0;
2570 int width=1;
2571
2572 for(int i=0; i< n*m; i++)
2573 max=std::max(max, aggregates[i]);
2574
2575 for(int i=10; i < 1000000; i*=10)
2576 if(max/i>0)
2577 width++;
2578 else
2579 break;
2580
2581 for(int j=0, entry=0; j < m; j++) {
2582 for(int i=0; i<n; i++, entry++) {
2583 os.width(width);
2584 os<<aggregates[entry]<<" ";
2585 }
2586
2587 os<<std::endl;
2588 }
2589 os<<std::endl;
2590 os.flags(oldOpts);
2591 }
2592
2593
2594 } // namespace Amg
2595
2596} // namespace Dune
2597
2598
2599#endif
A class for temporarily storing the vertices of an aggregate in.
Definition: aggregates.hh:773
A Dummy visitor that does nothing for each visited edge.
Definition: aggregates.hh:581
Class providing information about the mapping of the vertices onto aggregates.
Definition: aggregates.hh:543
Base class of all aggregation criterions.
Definition: aggregates.hh:47
Class for building the aggregates.
Definition: aggregates.hh:904
Dependency policy for symmetric matrices.
Definition: aggregates.hh:249
Norm that uses only the [N][N] entry of the block to determine couplings.
Definition: aggregates.hh:372
Norm that uses only the [0][0] entry of the block to determine couplings.
Definition: aggregates.hh:432
Iterator over all edges starting from a vertex.
Definition: graph.hh:93
The vertex iterator type of the graph.
Definition: graph.hh:207
The (undirected) graph of a matrix.
Definition: graph.hh:49
M::size_type VertexDescriptor
The vertex descriptor.
Definition: graph.hh:71
EdgeIteratorT< const MatrixGraph< Matrix > > ConstEdgeIterator
The constant edge iterator type.
Definition: graph.hh:296
All parameters for AMG.
Definition: parameters.hh:391
Criterion taking advantage of symmetric matrices.
Definition: aggregates.hh:499
Dependency policy for symmetric matrices.
Definition: aggregates.hh:310
Dependency policy for symmetric matrices.
Definition: aggregates.hh:132
Criterion suited for unsymmetric matrices.
Definition: aggregates.hh:522
Definition: bvector.hh:650
derive error class from the base class in common
Definition: istlexception.hh:16
A generic dynamic dense matrix.
Definition: matrix.hh:555
T::field_type field_type
Export the type representing the underlying field.
Definition: matrix.hh:559
row_type::const_iterator ConstColIterator
Const iterator for the entries of each row.
Definition: matrix.hh:583
MatrixImp::DenseMatrixBase< T, A >::window_type row_type
The type implementing a matrix row.
Definition: matrix.hh:568
An allocator managing a pool of objects for reuse.
Definition: poolallocator.hh:247
A single linked list.
Definition: sllist.hh:42
Type traits to determine the type of reals (when working with complex numbers)
Provides classes for building the matrix graph.
SLListConstIterator< T, A > const_iterator
The constant iterator of the list.
Definition: sllist.hh:72
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
Matrix::ConstColIterator ColIter
Constant column iterator of the matrix.
Definition: aggregates.hh:269
Matrix::ConstColIterator ColIter
Constant column iterator of the matrix.
Definition: aggregates.hh:152
std::size_t breadthFirstSearch(const VertexDescriptor &start, const AggregateDescriptor &aggregate, const G &graph, L &visited, F1 &aggregateVisitor, F2 &nonAggregateVisitor, VM &visitedMap) const
Breadth first search within an aggregate.
PoolAllocator< VertexDescriptor, 100 > Allocator
The allocator we use for our lists and the set.
Definition: aggregates.hh:569
int id()
Get the id identifying the aggregate.
Norm norm_
The functor for calculating the norm.
Definition: aggregates.hh:298
FieldTraits< typenameM::field_type >::real_type operator()(const M &m) const
compute the norm of a matrix.
Definition: aggregates.hh:383
MatrixGraph::VertexDescriptor Vertex
The vertex identifier.
Definition: aggregates.hh:915
AggregationCriterion()
Constructor.
Definition: aggregates.hh:64
const Matrix * matrix_
The matrix we work on.
Definition: aggregates.hh:353
M Matrix
The matrix type we build the dependency of.
Definition: aggregates.hh:254
G MatrixGraph
The matrix graph type used.
Definition: aggregates.hh:910
Norm norm_
The functor for calculating the norm.
Definition: aggregates.hh:359
M Matrix
The matrix type we build the dependency of.
Definition: aggregates.hh:315
real_type diagonal_
The norm of the current diagonal.
Definition: aggregates.hh:185
N Norm
The norm to use for examining the matrix entries.
Definition: aggregates.hh:259
int row_
index of the currently evaluated row.
Definition: aggregates.hh:183
std::tuple< int, int, int, int > build(const M &m, G &graph, AggregatesMap< Vertex > &aggregates, const C &c, bool finestLevel)
Build the aggregates.
FrontNeighbourCounter(const MatrixGraph &front)
Constructor.
Matrix::row_type Row
Constant Row iterator of the matrix.
Definition: aggregates.hh:325
const Matrix * matrix_
The matrix we work on.
Definition: aggregates.hh:292
const AggregateDescriptor & operator[](const VertexDescriptor &v) const
Get the aggregate a vertex belongs to.
AggregateVisitor(const AggregatesMap< Vertex > &aggregates, const AggregateDescriptor &aggregate, Visitor &visitor)
Constructor.
Matrix::ConstColIterator ColIter
Constant column iterator of the matrix.
Definition: aggregates.hh:330
~AggregatesMap()
Destructor.
Matrix::field_type field_type
The current max value.
Definition: aggregates.hh:177
void decrement()
Decrement counter.
Aggregate(MatrixGraph &graph, AggregatesMap< Vertex > &aggregates, VertexSet &connectivity, std::vector< Vertex > &front_)
Constructor.
V Visitor
The type of the adapted visitor.
Definition: aggregates.hh:1059
std::size_t * SphereMap
Type of the mapping of aggregate members onto distance spheres.
Definition: aggregates.hh:804
Matrix::row_type Row
Constant Row iterator of the matrix.
Definition: aggregates.hh:264
VertexSet::size_type connectSize()
Get tne number of connections to other aggregates.
N Norm
The norm to use for examining the matrix entries.
Definition: aggregates.hh:320
Norm norm_
The functor for calculating the norm.
Definition: aggregates.hh:181
VertexSet::const_iterator const_iterator
Const iterator over a vertex list.
Definition: aggregates.hh:799
MatrixGraph::VertexDescriptor AggregateDescriptor
The type of the aggregate descriptor.
Definition: aggregates.hh:918
real_type diagonal_
The norm of the current diagonal.
Definition: aggregates.hh:363
~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:796
FieldTraits< typenameM::field_type >::real_type operator()(const M &m) const
compute the norm of a matrix.
Definition: aggregates.hh:467
static const V UNAGGREGATED
Identifier of not yet aggregated vertices.
Definition: aggregates.hh:549
std::size_t breadthFirstSearch(const VertexDescriptor &start, const AggregateDescriptor &aggregate, const G &graph, F &aggregateVisitor, VM &visitedMap) const
Breadth first search within an aggregate.
Matrix::field_type field_type
The current max value.
Definition: aggregates.hh:294
void allocate(std::size_t noVertices)
Allocate memory for holding the information.
N Norm
The norm to use for examining the matrix entries.
Definition: aggregates.hh:142
FieldTraits< typenameM::field_type >::real_type operator()(const M &m) const
compute the norm of a matrix.
Definition: aggregates.hh:450
void reconstruct(const Vertex &vertex)
Reconstruct the aggregat from an seed node.
const_iterator begin() const
get an iterator over the vertices of the aggregate.
MatrixGraph::VertexDescriptor Vertex
The vertex descriptor type.
Definition: aggregates.hh:784
void seed(const Vertex &vertex)
Initialize the aggregate with one vertex.
void clear()
Clear the aggregate.
void free()
Free the allocated memory.
void increment()
Increment counter.
void buildDependency(G &graph, const typename C::Matrix &matrix, C criterion, bool finestLevel)
Build the dependency of the matrix graph.
V VertexDescriptor
The vertex descriptor type.
Definition: aggregates.hh:558
std::tuple< int, int, int, int > buildAggregates(const M &matrix, G &graph, const C &criterion, bool finestLevel)
Build the aggregates.
Matrix::row_type Row
Constant Row iterator of the matrix.
Definition: aggregates.hh:147
PoolAllocator< Vertex, 100 > Allocator
The allocator we use for our lists and the set.
Definition: aggregates.hh:790
DVVerbType dvverb(std::cout)
stream for very verbose output.
Definition: stdstreams.hh:93
DInfoType dinfo(std::cout)
Stream for informative output.
Definition: stdstreams.hh:138
DVerbType dverb(std::cout)
Singleton of verbose debug stream.
Definition: stdstreams.hh:114
Front front
PartitionSet for the front partition.
Definition: partitionset.hh:229
Dune namespace.
Definition: alignment.hh:11
Parameter classes for customizing AMG.
An stl-compliant pool allocator.
Provides classes for handling internal properties in a graph.
Implements a singly linked list together with the necessary iterators.
Standard Dune debug streams.
Functor using the row sum (infinity) norm to determine strong couplings.
Definition: aggregates.hh:440
A simple timing class.
Definition of the DUNE_UNUSED macro for the case that config.h is not available.
#define DUNE_UNUSED_PARAMETER(parm)
A macro to mark intentionally unused function parameters with.
Definition: unused.hh:18
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 23, 23:29, 2024)