dune-istl
2.2.0
|
00001 // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*- 00002 // vi: set ts=8 sw=2 et sts=2: 00003 // $Id: aggregates.hh 1532 2012-02-13 16:38:54Z christi $ 00004 #ifndef DUNE_AMG_AGGREGATES_HH 00005 #define DUNE_AMG_AGGREGATES_HH 00006 00007 00008 #include"parameters.hh" 00009 #include"graph.hh" 00010 #include"properties.hh" 00011 #include"combinedfunctor.hh" 00012 00013 #include<dune/common/timer.hh> 00014 #include<dune/common/tuples.hh> 00015 #include<dune/common/stdstreams.hh> 00016 #include<dune/common/poolallocator.hh> 00017 #include<dune/common/sllist.hh> 00018 00019 #include<utility> 00020 #include<set> 00021 #include<algorithm> 00022 #include<limits> 00023 #include<ostream> 00024 00025 namespace Dune 00026 { 00027 namespace Amg 00028 { 00029 00043 template<class T> 00044 class AggregationCriterion : public T 00045 { 00046 00047 public: 00051 typedef T DependencyPolicy; 00052 00062 AggregationCriterion() 00063 : T() 00064 {} 00065 00066 AggregationCriterion(const Parameters& parms) 00067 : T(parms) 00068 {} 00078 void setDefaultValuesIsotropic(std::size_t dim, std::size_t diameter=2) 00079 { 00080 this->setMaxDistance(diameter-1); 00081 std::size_t csize=1; 00082 00083 for(;dim>0;dim--){ 00084 csize*=diameter; 00085 this->setMaxDistance(this->maxDistance()+diameter-1); 00086 } 00087 this->setMinAggregateSize(csize); 00088 this->setMaxAggregateSize(static_cast<std::size_t>(csize*1.5)); 00089 } 00090 00101 void setDefaultValuesAnisotropic(std::size_t dim,std::size_t diameter=2) 00102 { 00103 setDefaultValuesIsotropic(dim, diameter); 00104 this->setMaxDistance(this->maxDistance()+dim-1); 00105 } 00106 }; 00107 00108 template<class T> 00109 std::ostream& operator<<(std::ostream& os, const AggregationCriterion<T>& criterion) 00110 { 00111 os<<"{ maxdistance="<<criterion.maxDistance()<<" minAggregateSize=" 00112 <<criterion.minAggregateSize()<< " maxAggregateSize="<<criterion.maxAggregateSize() 00113 <<" connectivity="<<criterion.maxConnectivity()<<" debugLevel="<<criterion.debugLevel()<<"}"; 00114 return os; 00115 } 00116 00120 template<class M, class N> 00121 class Dependency : public Parameters 00122 { 00123 public: 00127 typedef M Matrix; 00128 00132 typedef N Norm; 00133 00137 typedef typename Matrix::row_type Row; 00138 00142 typedef typename Matrix::ConstColIterator ColIter; 00143 00144 void init(const Matrix* matrix); 00145 00146 void initRow(const Row& row, int index); 00147 00148 void examine(const ColIter& col); 00149 00150 template<class G> 00151 void examine(G& graph, const typename G::EdgeIterator& edge, const ColIter& col); 00152 00153 bool isIsolated(); 00154 00155 Dependency(const Parameters& parms) 00156 : Parameters(parms) 00157 {} 00158 00159 Dependency() 00160 :Parameters() 00161 {} 00162 00163 protected: 00165 const Matrix* matrix_; 00167 typename Matrix::field_type maxValue_; 00169 Norm norm_; 00171 int row_; 00173 typename Matrix::field_type diagonal_; 00174 }; 00175 00179 template<class M, class N> 00180 class SymmetricDependency : public Parameters 00181 { 00182 public: 00186 typedef M Matrix; 00187 00191 typedef N Norm; 00192 00196 typedef typename Matrix::row_type Row; 00197 00201 typedef typename Matrix::ConstColIterator ColIter; 00202 00203 void init(const Matrix* matrix); 00204 00205 void initRow(const Row& row, int index); 00206 00207 void examine(const ColIter& col); 00208 00209 template<class G> 00210 void examine(G& graph, const typename G::EdgeIterator& edge, const ColIter& col); 00211 00212 bool isIsolated(); 00213 00214 00215 SymmetricDependency(const Parameters& parms) 00216 : Parameters(parms) 00217 {} 00218 SymmetricDependency() 00219 :Parameters() 00220 {} 00221 00222 protected: 00224 const Matrix* matrix_; 00226 typename Matrix::field_type maxValue_; 00228 Norm norm_; 00230 int row_; 00232 typename Matrix::field_type diagonal_; 00233 }; 00234 00239 template<int N> 00240 class Diagonal 00241 { 00242 public: 00243 enum{ /* @brief We preserve the sign.*/ 00244 is_sign_preserving = true 00245 }; 00246 00251 template<class M> 00252 typename M::field_type operator()(const M& m) const 00253 { 00254 return m[N][N]; 00255 } 00256 }; 00257 00262 class FirstDiagonal : public Diagonal<0> 00263 {}; 00264 00270 struct RowSum 00271 { 00272 00273 enum{ /* @brief We preserve the sign.*/ 00274 is_sign_preserving = false 00275 }; 00280 template<class M> 00281 typename M::field_type operator()(const M& m) const 00282 { 00283 return m.infinity_norm(); 00284 } 00285 }; 00286 00287 struct FrobeniusNorm 00288 { 00289 00290 enum{ /* @brief We preserve the sign.*/ 00291 is_sign_preserving = false 00292 }; 00297 template<class M> 00298 typename M::field_type operator()(const M& m) const 00299 { 00300 return m.frobenius_norm(); 00301 } 00302 }; 00303 struct AlwaysOneNorm 00304 { 00305 00306 enum{ /* @brief We preserve the sign.*/ 00307 is_sign_preserving = false 00308 }; 00313 template<class M> 00314 typename M::field_type operator()(const M& m) const 00315 { 00316 return 1; 00317 } 00318 }; 00328 template<class M, class Norm> 00329 class SymmetricCriterion : public AggregationCriterion<SymmetricDependency<M,Norm> > 00330 { 00331 public: 00332 SymmetricCriterion(const Parameters& parms) 00333 : AggregationCriterion<SymmetricDependency<M,Norm> >(parms) 00334 {} 00335 SymmetricCriterion() 00336 {} 00337 }; 00338 00339 00351 template<class M, class Norm> 00352 class UnSymmetricCriterion : public AggregationCriterion<Dependency<M,Norm> > 00353 { 00354 public: 00355 UnSymmetricCriterion(const Parameters& parms) 00356 : AggregationCriterion<Dependency<M,Norm> >(parms) 00357 {} 00358 UnSymmetricCriterion() 00359 {} 00360 }; 00361 // forward declaration 00362 template<class G> class Aggregator; 00363 00364 00372 template<class V> 00373 class AggregatesMap 00374 { 00375 public: 00376 00380 static const V UNAGGREGATED; 00381 00385 static const V ISOLATED; 00389 typedef V VertexDescriptor; 00390 00394 typedef V AggregateDescriptor; 00395 00400 typedef PoolAllocator<VertexDescriptor,100> Allocator; 00401 00406 typedef SLList<VertexDescriptor,Allocator> VertexList; 00407 00411 class DummyEdgeVisitor 00412 { 00413 public: 00414 template<class EdgeIterator> 00415 void operator()(const EdgeIterator& egde) const 00416 {} 00417 }; 00418 00419 00423 AggregatesMap(); 00424 00430 AggregatesMap(std::size_t noVertices); 00431 00435 ~AggregatesMap(); 00436 00448 template<class M, class G, class C> 00449 tuple<int,int,int,int> buildAggregates(const M& matrix, G& graph, const C& criterion, 00450 bool finestLevel); 00451 00471 template<bool reset, class G, class F, class VM> 00472 std::size_t breadthFirstSearch(const VertexDescriptor& start, 00473 const AggregateDescriptor& aggregate, 00474 const G& graph, 00475 F& aggregateVisitor, 00476 VM& visitedMap) const; 00477 00501 template<bool remove, bool reset, class G, class L, class F1, class F2, class VM> 00502 std::size_t breadthFirstSearch(const VertexDescriptor& start, 00503 const AggregateDescriptor& aggregate, 00504 const G& graph, L& visited, F1& aggregateVisitor, 00505 F2& nonAggregateVisitor, 00506 VM& visitedMap) const; 00507 00513 void allocate(std::size_t noVertices); 00514 00518 std::size_t noVertices() const; 00519 00523 void free(); 00524 00530 AggregateDescriptor& operator[](const VertexDescriptor& v); 00531 00537 const AggregateDescriptor& operator[](const VertexDescriptor& v) const; 00538 00539 typedef const AggregateDescriptor* const_iterator; 00540 00541 const_iterator begin() const 00542 { 00543 return aggregates_; 00544 } 00545 00546 const_iterator end() const 00547 { 00548 return aggregates_+noVertices(); 00549 } 00550 00551 typedef AggregateDescriptor* iterator; 00552 00553 iterator begin() 00554 { 00555 return aggregates_; 00556 } 00557 00558 iterator end() 00559 { 00560 return aggregates_+noVertices(); 00561 } 00562 private: 00564 AggregatesMap(const AggregatesMap<V>& map) 00565 { 00566 throw "Auch!"; 00567 } 00568 00570 AggregatesMap<V>& operator=(const AggregatesMap<V>& map) 00571 { 00572 throw "Auch!"; 00573 return this; 00574 } 00575 00579 AggregateDescriptor* aggregates_; 00580 00584 std::size_t noVertices_; 00585 }; 00586 00590 template<class G, class C> 00591 void buildDependency(G& graph, 00592 const typename C::Matrix& matrix, 00593 C criterion, 00594 bool finestLevel); 00595 00600 template<class G, class S> 00601 class Aggregate 00602 { 00603 00604 public: 00605 00606 /*** 00607 * @brief The type of the matrix graph we work with. 00608 */ 00609 typedef G MatrixGraph; 00613 typedef typename MatrixGraph::VertexDescriptor Vertex; 00614 00619 typedef PoolAllocator<Vertex,100> Allocator; 00620 00625 typedef SLList<Vertex,Allocator> VertexList; 00626 00627 00632 typedef S VertexSet; 00633 00635 typedef typename VertexList::const_iterator const_iterator; 00636 00640 typedef std::size_t* SphereMap; 00641 00649 Aggregate(const MatrixGraph& graph, AggregatesMap<Vertex>& aggregates, 00650 VertexSet& connectivity); 00651 00658 void reconstruct(const Vertex& vertex); 00659 00663 void seed(const Vertex& vertex); 00664 00668 void add(const Vertex& vertex); 00669 00673 void clear(); 00674 00678 typename VertexList::size_type size(); 00682 typename VertexList::size_type connectSize(); 00683 00687 int id(); 00688 00690 const_iterator begin() const; 00691 00693 const_iterator end() const; 00694 00695 private: 00699 VertexList vertices_; 00700 00705 int id_; 00706 00710 const MatrixGraph& graph_; 00711 00715 AggregatesMap<Vertex>& aggregates_; 00716 00720 VertexSet& connected_; 00721 }; 00722 00726 template<class G> 00727 class Aggregator 00728 { 00729 public: 00730 00734 typedef G MatrixGraph; 00735 00739 typedef typename MatrixGraph::VertexDescriptor Vertex; 00740 00742 typedef typename MatrixGraph::VertexDescriptor AggregateDescriptor; 00743 00747 Aggregator(); 00748 00752 ~Aggregator(); 00753 00770 template<class M, class C> 00771 tuple<int,int,int,int> build(const M& m, G& graph, 00772 AggregatesMap<Vertex>& aggregates, const C& c, 00773 bool finestLevel); 00774 private: 00779 typedef PoolAllocator<Vertex,100> Allocator; 00780 00784 typedef SLList<Vertex,Allocator> VertexList; 00785 00789 typedef std::set<Vertex,std::less<Vertex>,Allocator> VertexSet; 00790 00794 typedef std::size_t* SphereMap; 00795 00799 MatrixGraph* graph_; 00800 00804 Aggregate<MatrixGraph,VertexSet>* aggregate_; 00805 00809 VertexList front_; 00810 00814 VertexSet connected_; 00815 00819 int size_; 00820 00824 class Stack 00825 { 00826 public: 00827 static const Vertex NullEntry; 00828 00829 Stack(const MatrixGraph& graph, 00830 const Aggregator<G>& aggregatesBuilder, 00831 const AggregatesMap<Vertex>& aggregates); 00832 ~Stack(); 00833 bool push(const Vertex& v); 00834 void fill(); 00835 Vertex pop(); 00836 private: 00837 enum{ N = 256000 }; 00838 00840 const MatrixGraph& graph_; 00842 const Aggregator<G>& aggregatesBuilder_; 00844 const AggregatesMap<Vertex>& aggregates_; 00846 int size_; 00847 int maxSize_; 00849 int head_; 00850 int filled_; 00851 00853 Vertex* vals_; 00854 00855 void localPush(const Vertex& v); 00856 }; 00857 00858 friend class Stack; 00859 00870 template<class V> 00871 void visitAggregateNeighbours(const Vertex& vertex, const AggregateDescriptor& aggregate, 00872 const AggregatesMap<Vertex>& aggregates, 00873 V& visitor) const; 00874 00879 template<class V> 00880 class AggregateVisitor 00881 { 00882 public: 00886 typedef V Visitor; 00894 AggregateVisitor(const AggregatesMap<Vertex>& aggregates, const AggregateDescriptor& aggregate, 00895 Visitor& visitor); 00896 00903 void operator()(const typename MatrixGraph::ConstEdgeIterator& edge); 00904 00905 private: 00907 const AggregatesMap<Vertex>& aggregates_; 00909 AggregateDescriptor aggregate_; 00911 Visitor* visitor_; 00912 }; 00913 00917 class Counter 00918 { 00919 public: 00921 Counter(); 00923 int value(); 00924 00925 protected: 00927 void increment(); 00929 void decrement(); 00930 00931 private: 00932 int count_; 00933 }; 00934 00935 00940 class FrontNeighbourCounter : public Counter 00941 { 00942 public: 00947 FrontNeighbourCounter(const MatrixGraph& front); 00948 00949 void operator()(const typename MatrixGraph::ConstEdgeIterator& edge); 00950 00951 private: 00952 const MatrixGraph& graph_; 00953 }; 00954 00959 int noFrontNeighbours(const Vertex& vertex) const; 00960 00964 class TwoWayCounter : public Counter 00965 { 00966 public: 00967 void operator()(const typename MatrixGraph::ConstEdgeIterator& edge); 00968 }; 00969 00981 int twoWayConnections(const Vertex&, const AggregateDescriptor& aggregate, 00982 const AggregatesMap<Vertex>& aggregates) const; 00983 00987 class OneWayCounter : public Counter 00988 { 00989 public: 00990 void operator()(const typename MatrixGraph::ConstEdgeIterator& edge); 00991 }; 00992 01004 int oneWayConnections(const Vertex&, const AggregateDescriptor& aggregate, 01005 const AggregatesMap<Vertex>& aggregates) const; 01006 01013 class ConnectivityCounter : public Counter 01014 { 01015 public: 01022 ConnectivityCounter(const VertexSet& connected, const AggregatesMap<Vertex>& aggregates); 01023 01024 void operator()(const typename MatrixGraph::ConstEdgeIterator& edge); 01025 01026 private: 01028 const VertexSet& connected_; 01030 const AggregatesMap<Vertex>& aggregates_; 01031 01032 }; 01033 01045 double connectivity(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates) const; 01053 bool connected(const Vertex& vertex, const AggregateDescriptor& aggregate, 01054 const AggregatesMap<Vertex>& aggregates) const; 01055 01063 bool connected(const Vertex& vertex, const SLList<AggregateDescriptor>& aggregateList, 01064 const AggregatesMap<Vertex>& aggregates) const; 01065 01073 class DependencyCounter: public Counter 01074 { 01075 public: 01079 DependencyCounter(); 01080 01081 void operator()(const typename MatrixGraph::ConstEdgeIterator& edge); 01082 }; 01083 01090 class FrontMarker 01091 { 01092 public: 01099 FrontMarker(VertexList& front, MatrixGraph& graph); 01100 01101 void operator()(const typename MatrixGraph::ConstEdgeIterator& edge); 01102 01103 private: 01105 VertexList& front_; 01107 MatrixGraph& graph_; 01108 }; 01109 01116 void markFront(const AggregatesMap<Vertex>& aggregates); 01117 01121 void unmarkFront(); 01122 01137 int unusedNeighbours(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates) const; 01138 01152 std::pair<int,int> neighbours(const Vertex& vertex, 01153 const AggregateDescriptor& aggregate, 01154 const AggregatesMap<Vertex>& aggregates) const; 01171 int aggregateNeighbours(const Vertex& vertex, const AggregateDescriptor& aggregate, const AggregatesMap<Vertex>& aggregates) const; 01172 01180 bool admissible(const Vertex& vertex, const AggregateDescriptor& aggregate, const AggregatesMap<Vertex>& aggregates) const; 01181 01188 void seedFromFront(Stack& stack, bool isolated); 01189 01197 std::size_t distance(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates); 01198 01207 Vertex mergeNeighbour(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates) const; 01208 01217 void nonisoNeighbourAggregate(const Vertex& vertex, 01218 const AggregatesMap<Vertex>& aggregates, 01219 SLList<Vertex>& neighbours) const; 01220 01228 template<class C> 01229 void growAggregate(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates, const C& c); 01230 template<class C> 01231 void growIsolatedAggregate(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates, const C& c); 01232 }; 01233 01234 #ifndef DOXYGEN 01235 01236 template<class M, class N> 01237 inline void SymmetricDependency<M,N>::init(const Matrix* matrix) 01238 { 01239 matrix_ = matrix; 01240 } 01241 01242 template<class M, class N> 01243 inline void SymmetricDependency<M,N>::initRow(const Row& row, int index) 01244 { 01245 maxValue_ = std::min(- std::numeric_limits<typename Matrix::field_type>::max(), std::numeric_limits<typename Matrix::field_type>::min()); 01246 row_ = index; 01247 diagonal_ = norm_(matrix_->operator[](row_)[row_]); 01248 } 01249 01250 template<class M, class N> 01251 inline void SymmetricDependency<M,N>::examine(const ColIter& col) 01252 { 01253 typename Matrix::field_type eij = norm_(*col); 01254 typename Matrix::field_type eji = norm_(matrix_->operator[](col.index())[row_]); 01255 01256 // skip positive offdiagonals if norm preserves sign of them. 01257 if(!N::is_sign_preserving || eij<0 || eji<0) 01258 maxValue_ = std::max(maxValue_, 01259 eij /diagonal_ * eji/ 01260 norm_(matrix_->operator[](col.index())[col.index()])); 01261 } 01262 01263 template<class M, class N> 01264 template<class G> 01265 inline void SymmetricDependency<M,N>::examine(G& graph, const typename G::EdgeIterator& edge, const ColIter& col) 01266 { 01267 typename Matrix::field_type eij = norm_(*col); 01268 typename Matrix::field_type eji = norm_(matrix_->operator[](col.index())[row_]); 01269 01270 // skip positve offdiagonals if norm preserves sign of them. 01271 if(!N::is_sign_preserving || (eij<0 || eji<0)) 01272 if(eji / norm_(matrix_->operator[](edge.target())[edge.target()]) * 01273 eij/ diagonal_ > alpha() * maxValue_){ 01274 edge.properties().setDepends(); 01275 edge.properties().setInfluences(); 01276 01277 typename G::EdgeProperties& other = graph.getEdgeProperties(edge.target(), edge.source()); 01278 other.setInfluences(); 01279 other.setDepends(); 01280 } 01281 } 01282 01283 template<class M, class N> 01284 inline bool SymmetricDependency<M,N>::isIsolated() 01285 { 01286 return maxValue_ < beta(); 01287 } 01288 01289 01290 template<class M, class N> 01291 inline void Dependency<M,N>::init(const Matrix* matrix) 01292 { 01293 matrix_ = matrix; 01294 } 01295 01296 template<class M, class N> 01297 inline void Dependency<M,N>::initRow(const Row& row, int index) 01298 { 01299 maxValue_ = std::min(- std::numeric_limits<typename Matrix::field_type>::max(), std::numeric_limits<typename Matrix::field_type>::min()); 01300 row_ = index; 01301 diagonal_ = norm_(matrix_->operator[](row_)[row_]); 01302 } 01303 01304 template<class M, class N> 01305 inline void Dependency<M,N>::examine(const ColIter& col) 01306 { 01307 maxValue_ = std::max(maxValue_, 01308 -norm_(*col)); 01309 } 01310 01311 template<class M, class N> 01312 template<class G> 01313 inline void Dependency<M,N>::examine(G& graph, const typename G::EdgeIterator& edge, const ColIter& col) 01314 { 01315 if(-norm_(*col) >= maxValue_ * alpha()){ 01316 edge.properties().setDepends(); 01317 typedef typename G::EdgeDescriptor ED; 01318 ED e= graph.findEdge(edge.target(), edge.source()); 01319 if(e!=std::numeric_limits<ED>::max()) 01320 { 01321 typename G::EdgeProperties& other = graph.getEdgeProperties(e); 01322 other.setInfluences(); 01323 } 01324 } 01325 } 01326 01327 template<class M, class N> 01328 inline bool Dependency<M,N>::isIsolated() 01329 { 01330 return maxValue_ < beta() * diagonal_; 01331 } 01332 01333 template<class G,class S> 01334 Aggregate<G,S>::Aggregate(const MatrixGraph& graph, AggregatesMap<Vertex>& aggregates, 01335 VertexSet& connected) 01336 : vertices_(), id_(-1), graph_(graph), aggregates_(aggregates), 01337 connected_(connected) 01338 {} 01339 01340 template<class G,class S> 01341 void Aggregate<G,S>::reconstruct(const Vertex& vertex) 01342 { 01343 vertices_.push_back(vertex); 01344 typedef typename VertexList::const_iterator iterator; 01345 iterator begin = vertices_.begin(); 01346 iterator end = vertices_.end(); 01347 throw "Not yet implemented"; 01348 01349 while(begin!=end){ 01350 //for(); 01351 } 01352 01353 } 01354 01355 template<class G,class S> 01356 inline void Aggregate<G,S>::seed(const Vertex& vertex) 01357 { 01358 dvverb<<"Connected cleared"<<std::endl; 01359 connected_.clear(); 01360 vertices_.clear(); 01361 connected_.insert(vertex); 01362 dvverb << " Inserting "<<vertex<<" size="<<connected_.size(); 01363 id_ = vertex; 01364 add(vertex); 01365 } 01366 01367 01368 template<class G,class S> 01369 inline void Aggregate<G,S>::add(const Vertex& vertex) 01370 { 01371 vertices_.push_back(vertex); 01372 aggregates_[vertex]=id_; 01373 01374 typedef typename MatrixGraph::ConstEdgeIterator iterator; 01375 const iterator end = graph_.endEdges(vertex); 01376 for(iterator edge = graph_.beginEdges(vertex); edge != end; ++edge){ 01377 dvverb << " Inserting "<<aggregates_[edge.target()]; 01378 connected_.insert(aggregates_[edge.target()]); 01379 dvverb <<" size="<<connected_.size(); 01380 } 01381 dvverb <<std::endl; 01382 } 01383 template<class G,class S> 01384 inline void Aggregate<G,S>::clear() 01385 { 01386 vertices_.clear(); 01387 connected_.clear(); 01388 id_=-1; 01389 } 01390 01391 template<class G,class S> 01392 inline typename Aggregate<G,S>::VertexList::size_type 01393 Aggregate<G,S>::size() 01394 { 01395 return vertices_.size(); 01396 } 01397 01398 template<class G,class S> 01399 inline typename Aggregate<G,S>::VertexList::size_type 01400 Aggregate<G,S>::connectSize() 01401 { 01402 return connected_.size(); 01403 } 01404 01405 template<class G,class S> 01406 inline int Aggregate<G,S>::id() 01407 { 01408 return id_; 01409 } 01410 01411 template<class G,class S> 01412 inline typename Aggregate<G,S>::const_iterator Aggregate<G,S>::begin() const 01413 { 01414 return vertices_.begin(); 01415 } 01416 01417 template<class G,class S> 01418 inline typename Aggregate<G,S>::const_iterator Aggregate<G,S>::end() const 01419 { 01420 return vertices_.end(); 01421 } 01422 01423 template<class V> 01424 const V AggregatesMap<V>::UNAGGREGATED = std::numeric_limits<V>::max(); 01425 01426 template<class V> 01427 const V AggregatesMap<V>::ISOLATED = std::numeric_limits<V>::max()-1; 01428 01429 template<class V> 01430 AggregatesMap<V>::AggregatesMap() 01431 : aggregates_(0) 01432 {} 01433 01434 template<class V> 01435 AggregatesMap<V>::~AggregatesMap() 01436 { 01437 if(aggregates_!=0) 01438 delete[] aggregates_; 01439 } 01440 01441 01442 template<class V> 01443 inline AggregatesMap<V>::AggregatesMap(std::size_t noVertices) 01444 { 01445 allocate(noVertices); 01446 } 01447 01448 template<class V> 01449 inline std::size_t AggregatesMap<V>::noVertices() const 01450 { 01451 return noVertices_; 01452 } 01453 01454 template<class V> 01455 inline void AggregatesMap<V>::allocate(std::size_t noVertices) 01456 { 01457 aggregates_ = new AggregateDescriptor[noVertices]; 01458 noVertices_ = noVertices; 01459 01460 for(std::size_t i=0; i < noVertices; i++) 01461 aggregates_[i]=UNAGGREGATED; 01462 } 01463 01464 template<class V> 01465 inline void AggregatesMap<V>::free() 01466 { 01467 assert(aggregates_ != 0); 01468 delete[] aggregates_; 01469 aggregates_=0; 01470 } 01471 01472 template<class V> 01473 inline typename AggregatesMap<V>::AggregateDescriptor& 01474 AggregatesMap<V>::operator[](const VertexDescriptor& v) 01475 { 01476 return aggregates_[v]; 01477 } 01478 01479 template<class V> 01480 inline const typename AggregatesMap<V>::AggregateDescriptor& 01481 AggregatesMap<V>::operator[](const VertexDescriptor& v) const 01482 { 01483 return aggregates_[v]; 01484 } 01485 01486 template<class V> 01487 template<bool reset, class G, class F,class VM> 01488 inline std::size_t AggregatesMap<V>::breadthFirstSearch(const V& start, 01489 const AggregateDescriptor& aggregate, 01490 const G& graph, F& aggregateVisitor, 01491 VM& visitedMap) const 01492 { 01493 VertexList vlist; 01494 01495 DummyEdgeVisitor dummy; 01496 return breadthFirstSearch<true,reset>(start, aggregate, graph, vlist, aggregateVisitor, dummy, visitedMap); 01497 } 01498 01499 template<class V> 01500 template<bool remove, bool reset, class G, class L, class F1, class F2, class VM> 01501 std::size_t AggregatesMap<V>::breadthFirstSearch(const V& start, 01502 const AggregateDescriptor& aggregate, 01503 const G& graph, 01504 L& visited, 01505 F1& aggregateVisitor, 01506 F2& nonAggregateVisitor, 01507 VM& visitedMap) const 01508 { 01509 typedef typename L::const_iterator ListIterator; 01510 int visitedSpheres = 0; 01511 01512 visited.push_back(start); 01513 put(visitedMap, start, true); 01514 01515 ListIterator current = visited.begin(); 01516 ListIterator end = visited.end(); 01517 std::size_t i=0, size=visited.size(); 01518 01519 // visit the neighbours of all vertices of the 01520 // current sphere. 01521 while(current != end){ 01522 01523 for(;i<size; ++current, ++i){ 01524 typedef typename G::ConstEdgeIterator EdgeIterator; 01525 const EdgeIterator endEdge = graph.endEdges(*current); 01526 01527 for(EdgeIterator edge = graph.beginEdges(*current); 01528 edge != endEdge; ++edge){ 01529 01530 if(aggregates_[edge.target()]==aggregate){ 01531 if(!get(visitedMap, edge.target())){ 01532 put(visitedMap, edge.target(), true); 01533 visited.push_back(edge.target()); 01534 aggregateVisitor(edge); 01535 } 01536 }else 01537 nonAggregateVisitor(edge); 01538 } 01539 } 01540 end = visited.end(); 01541 size = visited.size(); 01542 if(current != end) 01543 visitedSpheres++; 01544 } 01545 01546 if(reset) 01547 for(current = visited.begin(); current != end; ++current) 01548 put(visitedMap, *current, false); 01549 01550 01551 if(remove) 01552 visited.clear(); 01553 01554 return visitedSpheres; 01555 } 01556 01557 template<class G> 01558 Aggregator<G>::Aggregator() 01559 : graph_(0), aggregate_(0), front_(), connected_(), size_(-1) 01560 {} 01561 01562 template<class G> 01563 Aggregator<G>::~Aggregator() 01564 { 01565 size_=-1; 01566 } 01567 01568 template<class G, class C> 01569 void buildDependency(G& graph, 01570 const typename C::Matrix& matrix, 01571 C criterion, bool firstlevel) 01572 { 01573 // The Criterion we use for building the dependency. 01574 typedef C Criterion; 01575 01576 // assert(graph.isBuilt()); 01577 typedef typename C::Matrix Matrix; 01578 typedef typename G::VertexIterator VertexIterator; 01579 01580 criterion.init(&matrix); 01581 01582 for(VertexIterator vertex = graph.begin(); vertex != graph.end(); ++vertex){ 01583 typedef typename Matrix::row_type Row; 01584 01585 const Row& row = matrix[*vertex]; 01586 01587 // Tell the criterion what row we will examine now 01588 // This might for example be used for calculating the 01589 // maximum offdiagonal value 01590 criterion.initRow(row, *vertex); 01591 01592 // On a first path all columns are examined. After this 01593 // the calculator should know whether the vertex is isolated. 01594 typedef typename Matrix::ConstColIterator ColIterator; 01595 ColIterator end = row.end(); 01596 typename Matrix::field_type absoffdiag=0; 01597 01598 if(firstlevel){ 01599 for(ColIterator col = row.begin(); col != end; ++col) 01600 if(col.index()!=*vertex){ 01601 criterion.examine(col); 01602 absoffdiag=std::max(absoffdiag, col->frobenius_norm()); 01603 } 01604 01605 if(absoffdiag==0) 01606 vertex.properties().setExcludedBorder(); 01607 } 01608 else 01609 for(ColIterator col = row.begin(); col != end; ++col) 01610 if(col.index()!=*vertex) 01611 criterion.examine(col); 01612 01613 // reset the vertex properties 01614 //vertex.properties().reset(); 01615 01616 // Check whether the vertex is isolated. 01617 if(criterion.isIsolated()){ 01618 //std::cout<<"ISOLATED: "<<*vertex<<std::endl; 01619 vertex.properties().setIsolated(); 01620 }else{ 01621 // Examine all the edges beginning at this vertex. 01622 typedef typename G::EdgeIterator EdgeIterator; 01623 typedef typename Matrix::ConstColIterator ColIterator; 01624 EdgeIterator end = vertex.end(); 01625 ColIterator col = matrix[*vertex].begin(); 01626 01627 for(EdgeIterator edge = vertex.begin(); edge!= end; ++edge, ++col){ 01628 // Move to the right column. 01629 while(col.index()!=edge.target()) 01630 ++col; 01631 criterion.examine(graph, edge, col); 01632 } 01633 } 01634 01635 } 01636 } 01637 01638 01639 template<class G> 01640 template<class V> 01641 inline Aggregator<G>::AggregateVisitor<V>::AggregateVisitor(const AggregatesMap<Vertex>& aggregates, 01642 const AggregateDescriptor& aggregate, V& visitor) 01643 : aggregates_(aggregates), aggregate_(aggregate), visitor_(&visitor) 01644 {} 01645 01646 template<class G> 01647 template<class V> 01648 inline void Aggregator<G>::AggregateVisitor<V>::operator()(const typename MatrixGraph::ConstEdgeIterator& edge) 01649 { 01650 if(aggregates_[edge.target()]==aggregate_) 01651 visitor_->operator()(edge); 01652 } 01653 01654 template<class G> 01655 template<class V> 01656 inline void Aggregator<G>::visitAggregateNeighbours(const Vertex& vertex, 01657 const AggregateDescriptor& aggregate, 01658 const AggregatesMap<Vertex>& aggregates, 01659 V& visitor) const 01660 { 01661 // Only evaluates for edge pointing to the aggregate 01662 AggregateVisitor<V> v(aggregates, aggregate, visitor); 01663 visitNeighbours(*graph_, vertex, v); 01664 } 01665 01666 01667 template<class G> 01668 inline Aggregator<G>::Counter::Counter() 01669 : count_(0) 01670 {} 01671 01672 template<class G> 01673 inline void Aggregator<G>::Counter::increment() 01674 { 01675 ++count_; 01676 } 01677 01678 template<class G> 01679 inline void Aggregator<G>::Counter::decrement() 01680 { 01681 --count_; 01682 } 01683 template<class G> 01684 inline int Aggregator<G>::Counter::value() 01685 { 01686 return count_; 01687 } 01688 01689 template<class G> 01690 inline void Aggregator<G>::TwoWayCounter::operator()(const typename MatrixGraph::ConstEdgeIterator& edge) 01691 { 01692 if(edge.properties().isTwoWay()) 01693 Counter::increment(); 01694 } 01695 01696 template<class G> 01697 int Aggregator<G>::twoWayConnections(const Vertex& vertex, const AggregateDescriptor& aggregate, 01698 const AggregatesMap<Vertex>& aggregates) const 01699 { 01700 TwoWayCounter counter; 01701 visitAggregateNeighbours(vertex, aggregate, aggregates, counter); 01702 return counter.value(); 01703 } 01704 01705 template<class G> 01706 int Aggregator<G>::oneWayConnections(const Vertex& vertex, const AggregateDescriptor& aggregate, 01707 const AggregatesMap<Vertex>& aggregates) const 01708 { 01709 OneWayCounter counter; 01710 visitAggregateNeighbours(vertex, aggregate, aggregates, counter); 01711 return counter.value(); 01712 } 01713 01714 template<class G> 01715 inline void Aggregator<G>::OneWayCounter::operator()(const typename MatrixGraph::ConstEdgeIterator& edge) 01716 { 01717 if(edge.properties().isOneWay()) 01718 Counter::increment(); 01719 } 01720 01721 template<class G> 01722 inline Aggregator<G>::ConnectivityCounter::ConnectivityCounter(const VertexSet& connected, 01723 const AggregatesMap<Vertex>& aggregates) 01724 : Counter(), connected_(connected), aggregates_(aggregates) 01725 {} 01726 01727 01728 template<class G> 01729 inline void Aggregator<G>::ConnectivityCounter::operator()(const typename MatrixGraph::ConstEdgeIterator& edge) 01730 { 01731 if(connected_.find(aggregates_[edge.target()]) == connected_.end() || aggregates_[edge.target()]==AggregatesMap<Vertex>::UNAGGREGATED) 01732 // Would be a new connection 01733 Counter::increment(); 01734 else{ 01735 Counter::increment(); 01736 Counter::increment(); 01737 } 01738 } 01739 01740 template<class G> 01741 inline double Aggregator<G>::connectivity(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates) const 01742 { 01743 ConnectivityCounter counter(connected_, aggregates); 01744 double noNeighbours=visitNeighbours(*graph_, vertex, counter); 01745 return (double)counter.value()/noNeighbours; 01746 } 01747 01748 template<class G> 01749 inline Aggregator<G>::DependencyCounter::DependencyCounter() 01750 : Counter() 01751 {} 01752 01753 template<class G> 01754 inline void Aggregator<G>::DependencyCounter::operator()(const typename MatrixGraph::ConstEdgeIterator& edge) 01755 { 01756 if(edge.properties().depends()) 01757 Counter::increment(); 01758 if(edge.properties().influences()) 01759 Counter::increment(); 01760 } 01761 01762 template<class G> 01763 int Aggregator<G>::unusedNeighbours(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates) const 01764 { 01765 return aggregateNeighbours(vertex, AggregatesMap<Vertex>::UNAGGREGATED, aggregates); 01766 } 01767 01768 template<class G> 01769 std::pair<int,int> Aggregator<G>::neighbours(const Vertex& vertex, 01770 const AggregateDescriptor& aggregate, 01771 const AggregatesMap<Vertex>& aggregates) const 01772 { 01773 DependencyCounter unused, aggregated; 01774 typedef AggregateVisitor<DependencyCounter> Counter; 01775 typedef tuple<Counter,Counter> CounterTuple; 01776 CombinedFunctor<CounterTuple> visitors(CounterTuple(Counter(aggregates, AggregatesMap<Vertex>::UNAGGREGATED, unused), Counter(aggregates, aggregate, aggregated))); 01777 visitNeighbours(*graph_, vertex, visitors); 01778 return std::make_pair(unused.value(), aggregated.value()); 01779 } 01780 01781 01782 template<class G> 01783 int Aggregator<G>::aggregateNeighbours(const Vertex& vertex, const AggregateDescriptor& aggregate, const AggregatesMap<Vertex>& aggregates) const 01784 { 01785 DependencyCounter counter; 01786 visitAggregateNeighbours(vertex, aggregate, aggregates, counter); 01787 return counter.value(); 01788 } 01789 01790 template<class G> 01791 std::size_t Aggregator<G>::distance(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates) 01792 { 01793 typename PropertyMapTypeSelector<VertexVisitedTag,G>::Type visitedMap = get(VertexVisitedTag(), *graph_); 01794 VertexList vlist; 01795 typename AggregatesMap<Vertex>::DummyEdgeVisitor dummy; 01796 return aggregates.template breadthFirstSearch<true,true>(vertex, 01797 aggregate_->id(), *graph_, 01798 vlist, dummy, dummy, visitedMap); 01799 } 01800 01801 template<class G> 01802 inline Aggregator<G>::FrontMarker::FrontMarker(VertexList& front, MatrixGraph& graph) 01803 : front_(front), graph_(graph) 01804 {} 01805 01806 template<class G> 01807 inline void Aggregator<G>::FrontMarker::operator()(const typename MatrixGraph::ConstEdgeIterator& edge) 01808 { 01809 Vertex target = edge.target(); 01810 01811 if(!graph_.getVertexProperties(target).front()){ 01812 front_.push_back(target); 01813 graph_.getVertexProperties(target).setFront(); 01814 } 01815 } 01816 01817 01818 template<class G> 01819 void Aggregator<G>::markFront(const AggregatesMap<Vertex>& aggregates) 01820 { 01821 assert(front_.size()==0); 01822 FrontMarker frontBuilder(front_, *graph_); 01823 typedef typename Aggregate<G,VertexSet>::const_iterator Iterator; 01824 01825 for(Iterator vertex=aggregate_->begin(); vertex != aggregate_->end(); ++vertex) 01826 visitAggregateNeighbours(*vertex, AggregatesMap<Vertex>::UNAGGREGATED, aggregates, frontBuilder); 01827 01828 } 01829 01830 template<class G> 01831 inline bool Aggregator<G>::admissible(const Vertex& vertex, const AggregateDescriptor& aggregate, const AggregatesMap<Vertex>& aggregates) const 01832 { 01833 // Todo 01834 Dune::dvverb<<" Admissible not yet implemented!"<<std::endl; 01835 return true; 01836 //Situation 1: front node depends on two nodes. Then these 01837 // have to be strongly connected to each other 01838 01839 // Iterate over all all neighbours of front node 01840 typedef typename MatrixGraph::ConstEdgeIterator Iterator; 01841 Iterator vend = graph_->endEdges(vertex); 01842 for(Iterator edge = graph_->beginEdges(vertex); edge != vend; ++edge){ 01843 // if(edge.properties().depends() && !edge.properties().influences() 01844 if(edge.properties().isStrong() 01845 && aggregates[edge.target()]==aggregate) 01846 { 01847 // Search for another link to the aggregate 01848 Iterator edge1 = edge; 01849 for(++edge1; edge1 != vend; ++edge1){ 01850 //if(edge1.properties().depends() && !edge1.properties().influences() 01851 if(edge1.properties().isStrong() 01852 && aggregates[edge.target()]==aggregate) 01853 { 01854 //Search for an edge connecting the two vertices that is 01855 //strong 01856 bool found=false; 01857 Iterator v2end = graph_->endEdges(edge.target()); 01858 for(Iterator edge2 = graph_->beginEdges(edge.target()); edge2 != v2end; ++edge2){ 01859 if(edge2.target()==edge1.target() && 01860 edge2.properties().isStrong()){ 01861 found =true; 01862 break; 01863 } 01864 } 01865 if(found) 01866 return true; 01867 } 01868 } 01869 } 01870 } 01871 01872 // Situation 2: cluster node depends on front node and other cluster node 01874 vend = graph_->endEdges(vertex); 01875 for(Iterator edge = graph_->beginEdges(vertex); edge != vend; ++edge){ 01876 //if(!edge.properties().depends() && edge.properties().influences() 01877 if(edge.properties().isStrong() 01878 && aggregates[edge.target()]==aggregate) 01879 { 01880 // Search for a link from target that stays within the aggregate 01881 Iterator v1end = graph_->endEdges(edge.target()); 01882 01883 for(Iterator edge1=graph_->beginEdges(edge.target()); edge1 != v1end; ++edge1){ 01884 //if(edge1.properties().depends() && !edge1.properties().influences() 01885 if(edge1.properties().isStrong() 01886 && aggregates[edge1.target()]==aggregate) 01887 { 01888 bool found=false; 01889 // Check if front node is also connected to this one 01890 Iterator v2end = graph_->endEdges(vertex); 01891 for(Iterator edge2 = graph_->beginEdges(vertex); edge2 != v2end; ++edge2){ 01892 if(edge2.target()==edge1.target()){ 01893 if(edge2.properties().isStrong()) 01894 found=true; 01895 break; 01896 } 01897 } 01898 if(found) 01899 return true; 01900 } 01901 } 01902 } 01903 } 01904 return false; 01905 } 01906 01907 template<class G> 01908 void Aggregator<G>::unmarkFront() 01909 { 01910 typedef typename VertexList::const_iterator Iterator; 01911 01912 for(Iterator vertex=front_.begin(); vertex != front_.end(); ++vertex) 01913 graph_->getVertexProperties(*vertex).resetFront(); 01914 01915 front_.clear(); 01916 } 01917 01918 template<class G> 01919 inline void 01920 Aggregator<G>::nonisoNeighbourAggregate(const Vertex& vertex, 01921 const AggregatesMap<Vertex>& aggregates, 01922 SLList<Vertex>& neighbours) const 01923 { 01924 typedef typename MatrixGraph::ConstEdgeIterator Iterator; 01925 Iterator end=graph_->beginEdges(vertex); 01926 neighbours.clear(); 01927 01928 for(Iterator edge=graph_->beginEdges(vertex); edge!=end; ++edge) 01929 { 01930 if(aggregates[edge.target()]!=AggregatesMap<Vertex>::UNAGGREGATED && graph_->getVertexProperties(edge.target()).isolated()) 01931 neighbours.push_back(aggregates[edge.target()]); 01932 } 01933 } 01934 01935 template<class G> 01936 inline typename G::VertexDescriptor Aggregator<G>::mergeNeighbour(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates) const 01937 { 01938 typedef typename MatrixGraph::ConstEdgeIterator Iterator; 01939 01940 Iterator end = graph_->endEdges(vertex); 01941 for(Iterator edge = graph_->beginEdges(vertex); edge != end; ++edge){ 01942 if(aggregates[edge.target()] != AggregatesMap<Vertex>::UNAGGREGATED && 01943 graph_->getVertexProperties(edge.target()).isolated() == graph_->getVertexProperties(edge.source()).isolated()){ 01944 if( graph_->getVertexProperties(vertex).isolated() || 01945 ((edge.properties().depends() || edge.properties().influences()) 01946 && admissible(vertex, aggregates[edge.target()], aggregates))) 01947 return edge.target(); 01948 } 01949 } 01950 return AggregatesMap<Vertex>::UNAGGREGATED; 01951 } 01952 01953 template<class G> 01954 Aggregator<G>::FrontNeighbourCounter::FrontNeighbourCounter(const MatrixGraph& graph) 01955 : Counter(), graph_(graph) 01956 {} 01957 01958 template<class G> 01959 void Aggregator<G>::FrontNeighbourCounter::operator()(const typename MatrixGraph::ConstEdgeIterator& edge) 01960 { 01961 if(graph_.getVertexProperties(edge.target()).front()) 01962 Counter::increment(); 01963 } 01964 01965 template<class G> 01966 int Aggregator<G>::noFrontNeighbours(const Vertex& vertex) const 01967 { 01968 FrontNeighbourCounter counter(*graph_); 01969 visitNeighbours(*graph_, vertex, counter); 01970 return counter.value(); 01971 } 01972 template<class G> 01973 inline bool Aggregator<G>::connected(const Vertex& vertex, 01974 const AggregateDescriptor& aggregate, 01975 const AggregatesMap<Vertex>& aggregates) const 01976 { 01977 typedef typename G::ConstEdgeIterator iterator; 01978 const iterator end = graph_->endEdges(vertex); 01979 for(iterator edge = graph_->beginEdges(vertex); edge != end; ++edge) 01980 if(aggregates[edge.target()]==aggregate) 01981 return true; 01982 return false; 01983 } 01984 template<class G> 01985 inline bool Aggregator<G>::connected(const Vertex& vertex, 01986 const SLList<AggregateDescriptor>& aggregateList, 01987 const AggregatesMap<Vertex>& aggregates) const 01988 { 01989 typedef typename SLList<AggregateDescriptor>::const_iterator Iter; 01990 for(Iter i=aggregateList.begin(); i!=aggregateList.end(); ++i) 01991 if(connected(vertex, *i, aggregates)) 01992 return true; 01993 return false; 01994 } 01995 01996 template<class G> 01997 template<class C> 01998 void Aggregator<G>::growIsolatedAggregate(const Vertex& seed, const AggregatesMap<Vertex>& aggregates, const C& c) 01999 { 02000 SLList<Vertex> connectedAggregates; 02001 nonisoNeighbourAggregate(seed, aggregates,connectedAggregates); 02002 02003 while(aggregate_->size()< c.minAggregateSize() && aggregate_->connectSize() < c.maxConnectivity()){ 02004 double maxCon=-1; 02005 std::size_t maxFrontNeighbours=0; 02006 02007 Vertex candidate=AggregatesMap<Vertex>::UNAGGREGATED; 02008 unmarkFront(); 02009 markFront(aggregates); 02010 typedef typename VertexList::const_iterator Iterator; 02011 02012 for(Iterator vertex = front_.begin(); vertex != front_.end(); ++vertex){ 02013 if(distance(*vertex, aggregates)>c.maxDistance()) 02014 continue; // distance of proposes aggregate too big 02015 02016 if(connectedAggregates.size()>0){ 02017 // there is already a neighbour cluster 02018 // front node must be connected to same neighbour cluster 02019 02020 if(!connected(*vertex, connectedAggregates, aggregates)) 02021 continue; 02022 } 02023 02024 double con = connectivity(*vertex, aggregates); 02025 02026 if(con == maxCon){ 02027 std::size_t frontNeighbours = noFrontNeighbours(*vertex); 02028 02029 if(frontNeighbours >= maxFrontNeighbours){ 02030 maxFrontNeighbours = frontNeighbours; 02031 candidate = *vertex; 02032 } 02033 }else if(con > maxCon){ 02034 maxCon = con; 02035 maxFrontNeighbours = noFrontNeighbours(*vertex); 02036 candidate = *vertex; 02037 } 02038 } 02039 02040 if(candidate==AggregatesMap<Vertex>::UNAGGREGATED) 02041 break; 02042 02043 aggregate_->add(candidate); 02044 } 02045 } 02046 02047 template<class G> 02048 template<class C> 02049 void Aggregator<G>::growAggregate(const Vertex& seed, const AggregatesMap<Vertex>& aggregates, const C& c) 02050 { 02051 while(aggregate_->size() < c.minAggregateSize()){ 02052 int maxTwoCons=0, maxOneCons=0, maxNeighbours=-1; 02053 double maxCon=-1; 02054 02055 Vertex candidate = AggregatesMap<Vertex>::UNAGGREGATED; 02056 02057 unmarkFront(); 02058 markFront(aggregates); 02059 02060 typedef typename VertexList::const_iterator Iterator; 02061 02062 for(Iterator vertex = front_.begin(); vertex != front_.end(); ++vertex){ 02063 // Only nonisolated nodes are considered 02064 if(graph_->getVertexProperties(*vertex).isolated()) 02065 continue; 02066 02067 int twoWayCons = twoWayConnections(*vertex, aggregate_->id(), aggregates); 02068 02069 /* The case of two way connections. */ 02070 if( maxTwoCons == twoWayCons && twoWayCons > 0){ 02071 double con = connectivity(*vertex, aggregates); 02072 02073 if(con == maxCon){ 02074 int neighbours = noFrontNeighbours(*vertex); 02075 02076 if(neighbours > maxNeighbours){ 02077 maxNeighbours = neighbours; 02078 02079 std::size_t distance_ = distance(*vertex, aggregates); 02080 02081 if(c.maxDistance() >= distance_){ 02082 candidate = *vertex; 02083 } 02084 } 02085 }else if( con > maxCon){ 02086 maxCon = con; 02087 maxNeighbours = noFrontNeighbours(*vertex); 02088 std::size_t distance_ = distance(*vertex, aggregates); 02089 02090 if(c.maxDistance() >= distance_){ 02091 candidate = *vertex; 02092 } 02093 } 02094 }else if(twoWayCons > maxTwoCons){ 02095 maxTwoCons = twoWayCons; 02096 maxCon = connectivity(*vertex, aggregates); 02097 maxNeighbours = noFrontNeighbours(*vertex); 02098 std::size_t distance_ = distance(*vertex, aggregates); 02099 02100 if(c.maxDistance() >= distance_){ 02101 candidate = *vertex; 02102 } 02103 02104 // two way connections preceed 02105 maxOneCons = std::numeric_limits<int>::max(); 02106 } 02107 02108 if(twoWayCons > 0) 02109 continue; // THis is a two-way node, skip tests for one way nodes 02110 02111 /* The one way case */ 02112 int oneWayCons = oneWayConnections(*vertex, aggregate_->id(), aggregates); 02113 02114 if(oneWayCons==0) 02115 continue; // No strong connections, skip the tests. 02116 02117 if(!admissible(*vertex, aggregate_->id(), aggregates)) 02118 continue; 02119 02120 if( maxOneCons == oneWayCons && oneWayCons > 0){ 02121 double con = connectivity(*vertex, aggregates); 02122 02123 if(con == maxCon){ 02124 int neighbours = noFrontNeighbours(*vertex); 02125 02126 if(neighbours > maxNeighbours){ 02127 maxNeighbours = neighbours; 02128 std::size_t distance_ = distance(*vertex, aggregates); 02129 02130 if(c.maxDistance() >= distance_){ 02131 candidate = *vertex; 02132 } 02133 } 02134 }else if( con > maxCon){ 02135 maxCon = con; 02136 maxNeighbours = noFrontNeighbours(*vertex); 02137 std::size_t distance_ = distance(*vertex, aggregates); 02138 if(c.maxDistance() >= distance_){ 02139 candidate = *vertex; 02140 } 02141 } 02142 }else if(oneWayCons > maxOneCons){ 02143 maxOneCons = oneWayCons; 02144 maxCon = connectivity(*vertex, aggregates); 02145 maxNeighbours = noFrontNeighbours(*vertex); 02146 std::size_t distance_ = distance(*vertex, aggregates); 02147 02148 if(c.maxDistance() >= distance_){ 02149 candidate = *vertex; 02150 } 02151 } 02152 } 02153 02154 02155 if(candidate == AggregatesMap<Vertex>::UNAGGREGATED) 02156 break; // No more candidates found 02157 02158 aggregate_->add(candidate); 02159 } 02160 } 02161 02162 template<typename V> 02163 template<typename M, typename G, typename C> 02164 tuple<int,int,int,int> AggregatesMap<V>::buildAggregates(const M& matrix, G& graph, const C& criterion, 02165 bool finestLevel) 02166 { 02167 Aggregator<G> aggregator; 02168 return aggregator.build(matrix, graph, *this, criterion, finestLevel); 02169 } 02170 02171 template<class G> 02172 template<class M, class C> 02173 tuple<int,int,int,int> Aggregator<G>::build(const M& m, G& graph, AggregatesMap<Vertex>& aggregates, const C& c, 02174 bool finestLevel) 02175 { 02176 // Stack for fast vertex access 02177 Stack stack_(graph, *this, aggregates); 02178 02179 graph_ = &graph; 02180 02181 aggregate_ = new Aggregate<G,VertexSet>(graph, aggregates, connected_); 02182 02183 Timer watch; 02184 watch.reset(); 02185 02186 buildDependency(graph, m, c, finestLevel); 02187 02188 dverb<<"Build dependency took "<< watch.elapsed()<<" seconds."<<std::endl; 02189 int noAggregates, conAggregates, isoAggregates, oneAggregates; 02190 std::size_t maxA=0, minA=1000000, avg=0; 02191 int skippedAggregates; 02192 noAggregates = conAggregates = isoAggregates = oneAggregates = 02193 skippedAggregates = 0; 02194 02195 while(true){ 02196 Vertex seed = stack_.pop(); 02197 02198 if(seed == Stack::NullEntry) 02199 // No more unaggregated vertices. We are finished! 02200 break; 02201 02202 // Debugging output 02203 if((noAggregates+1)%10000 == 0) 02204 Dune::dverb<<"c"; 02205 02206 aggregate_->seed(seed); 02207 02208 if(graph.getVertexProperties(seed).excludedBorder()){ 02209 aggregates[seed]=AggregatesMap<Vertex>::ISOLATED; 02210 ++skippedAggregates; 02211 continue; 02212 } 02213 02214 if(graph.getVertexProperties(seed).isolated()){ 02215 if(c.skipIsolated()){ 02216 // isolated vertices are not aggregated but skipped on the coarser levels. 02217 aggregates[seed]=AggregatesMap<Vertex>::ISOLATED; 02218 ++skippedAggregates; 02219 throw "huch"; 02220 // skip rest as no agglomeration is done. 02221 continue; 02222 }else 02223 growIsolatedAggregate(seed, aggregates, c); 02224 }else 02225 growAggregate(seed, aggregates, c); 02226 02227 02228 /* The rounding step. */ 02229 while(!(graph.getVertexProperties(seed).isolated()) && aggregate_->size() < c.maxAggregateSize()){ 02230 02231 unmarkFront(); 02232 markFront(aggregates); 02233 02234 Vertex candidate = AggregatesMap<Vertex>::UNAGGREGATED; 02235 02236 typedef typename VertexList::const_iterator Iterator; 02237 02238 for(Iterator vertex = front_.begin(); vertex != front_.end(); ++vertex){ 02239 02240 if(graph.getVertexProperties(*vertex).isolated()) 02241 continue; // No isolated nodes here 02242 02243 if(twoWayConnections( *vertex, aggregate_->id(), aggregates) == 0 && 02244 (oneWayConnections( *vertex, aggregate_->id(), aggregates) == 0 || 02245 !admissible( *vertex, aggregate_->id(), aggregates) )) 02246 continue; 02247 02248 std::pair<int,int> neighbourPair=neighbours(*vertex, aggregate_->id(), 02249 aggregates); 02250 02251 //if(aggregateNeighbours(*vertex, aggregate_->id(), aggregates) <= unusedNeighbours(*vertex, aggregates)) 02252 // continue; 02253 02254 if(neighbourPair.first >= neighbourPair.second) 02255 continue; 02256 02257 if(distance(*vertex, aggregates) > c.maxDistance()) 02258 continue; // Distance too far 02259 candidate = *vertex; 02260 break; 02261 } 02262 02263 if(candidate == AggregatesMap<Vertex>::UNAGGREGATED) break; // no more candidates found. 02264 02265 aggregate_->add(candidate); 02266 02267 } 02268 02269 // try to merge aggregates consisting of only one nonisolated vertex with other aggregates 02270 if(aggregate_->size()==1 && c.maxAggregateSize()>1){ 02271 if(!graph.getVertexProperties(seed).isolated()){ 02272 Vertex mergedNeighbour = mergeNeighbour(seed, aggregates); 02273 02274 if(mergedNeighbour != AggregatesMap<Vertex>::UNAGGREGATED){ 02275 // assign vertex to the neighbouring cluster 02276 aggregates[seed] = aggregates[mergedNeighbour]; 02277 }else{ 02278 minA=std::min(minA,static_cast<std::size_t>(1)); 02279 maxA=std::max(maxA,static_cast<std::size_t>(1)); 02280 ++oneAggregates; 02281 ++conAggregates; 02282 } 02283 }else{ 02284 ++avg; 02285 minA=std::min(minA,static_cast<std::size_t>(1)); 02286 maxA=std::max(maxA,static_cast<std::size_t>(1)); 02287 ++oneAggregates; 02288 ++isoAggregates; 02289 } 02290 ++avg; 02291 }else{ 02292 avg+=aggregate_->size(); 02293 minA=std::min(minA,aggregate_->size()); 02294 maxA=std::max(maxA,aggregate_->size()); 02295 if(graph.getVertexProperties(seed).isolated()) 02296 ++isoAggregates; 02297 else 02298 ++conAggregates; 02299 } 02300 unmarkFront(); 02301 markFront(aggregates); 02302 seedFromFront(stack_, graph.getVertexProperties(seed).isolated()); 02303 unmarkFront(); 02304 } 02305 02306 Dune::dinfo<<"connected aggregates: "<<conAggregates; 02307 Dune::dinfo<<" isolated aggregates: "<<isoAggregates; 02308 if(conAggregates+isoAggregates>0) 02309 Dune::dinfo<<" one node aggregates: "<<oneAggregates<<" min size=" 02310 <<minA<<" max size="<<maxA 02311 <<" avg="<<avg/(conAggregates+isoAggregates)<<std::endl; 02312 02313 delete aggregate_; 02314 return make_tuple(conAggregates+isoAggregates,isoAggregates, 02315 oneAggregates,skippedAggregates); 02316 } 02317 02318 template<class G> 02319 inline void Aggregator<G>::seedFromFront(Stack& stack_, bool isolated) 02320 { 02321 typedef typename VertexList::const_iterator Iterator; 02322 02323 Iterator end= front_.end(); 02324 int count=0; 02325 for(Iterator vertex=front_.begin(); vertex != end; ++vertex,++count) 02326 if(graph_->getVertexProperties(*vertex).isolated()==isolated) 02327 stack_.push(*vertex); 02328 /* 02329 if(MINIMAL_DEBUG_LEVEL<=2 && count==0 && !isolated) 02330 Dune::dverb<< " no vertices pushed for nonisolated aggregate!"<<std::endl; 02331 */ 02332 } 02333 02334 template<class G> 02335 Aggregator<G>::Stack::Stack(const MatrixGraph& graph, const Aggregator<G>& aggregatesBuilder, 02336 const AggregatesMap<Vertex>& aggregates) 02337 : graph_(graph), aggregatesBuilder_(aggregatesBuilder), aggregates_(aggregates), size_(0), maxSize_(0), head_(0), filled_(0) 02338 { 02339 vals_ = new Vertex[N]; 02340 } 02341 02342 template<class G> 02343 Aggregator<G>::Stack::~Stack() 02344 { 02345 Dune::dverb << "Max stack size was "<<maxSize_<<" filled="<<filled_<<std::endl; 02346 delete[] vals_; 02347 } 02348 02349 template<class G> 02350 const typename Aggregator<G>::Vertex Aggregator<G>::Stack::NullEntry 02351 = std::numeric_limits<typename G::VertexDescriptor>::max(); 02352 02353 template<class G> 02354 inline bool Aggregator<G>::Stack::push(const Vertex& v) 02355 { 02356 if(aggregates_[v] == AggregatesMap<Vertex>::UNAGGREGATED){ 02357 localPush(v); 02358 return true; 02359 }else 02360 return false; 02361 } 02362 02363 template<class G> 02364 inline void Aggregator<G>::Stack::localPush(const Vertex& v) 02365 { 02366 vals_[head_] = v; 02367 size_ = std::min<int>(size_+1, N); 02368 head_ = (head_+N+1)%N; 02369 } 02370 02371 template<class G> 02372 void Aggregator<G>::Stack::fill() 02373 { 02374 int isolated = 0, connected=0; 02375 int isoumin, umin; 02376 filled_++; 02377 02378 head_ = size_ = 0; 02379 isoumin = umin = std::numeric_limits<int>::max(); 02380 02381 typedef typename MatrixGraph::ConstVertexIterator Iterator; 02382 02383 const Iterator end = graph_.end(); 02384 02385 for(Iterator vertex = graph_.begin(); vertex != end; ++vertex){ 02386 // Skip already aggregated vertices 02387 if(aggregates_[*vertex] != AggregatesMap<Vertex>::UNAGGREGATED) 02388 continue; 02389 02390 if(vertex.properties().isolated()){ 02391 isoumin = std::min(isoumin, aggregatesBuilder_.unusedNeighbours(*vertex, aggregates_)); 02392 isolated++; 02393 }else{ 02394 umin = std::min(umin, aggregatesBuilder_.unusedNeighbours(*vertex, aggregates_)); 02395 connected++; 02396 } 02397 } 02398 02399 if(connected + isolated == 0) 02400 // No unaggregated vertices. 02401 return; 02402 02403 if(connected > 0){ 02404 // Connected vertices have higher priority. 02405 for(Iterator vertex = graph_.begin(); vertex != end; ++vertex) 02406 if(aggregates_[*vertex] == AggregatesMap<Vertex>::UNAGGREGATED && !vertex.properties().isolated() 02407 && aggregatesBuilder_.unusedNeighbours(*vertex, aggregates_) == umin) 02408 localPush(*vertex); 02409 }else{ 02410 for(Iterator vertex = graph_.begin(); vertex != end; ++vertex) 02411 if(aggregates_[*vertex] == AggregatesMap<Vertex>::UNAGGREGATED && vertex.properties().isolated() 02412 && aggregatesBuilder_.unusedNeighbours(*vertex, aggregates_) == isoumin) 02413 localPush(*vertex); 02414 } 02415 maxSize_ = std::max(size_, maxSize_); 02416 } 02417 02418 template<class G> 02419 inline typename G::VertexDescriptor Aggregator<G>::Stack::pop() 02420 { 02421 while(size_>0){ 02422 head_ = (head_ + N -1) % N; 02423 size_--; 02424 Vertex v = vals_[head_]; 02425 if(aggregates_[v]==AggregatesMap<Vertex>::UNAGGREGATED) 02426 return v; 02427 } 02428 // Stack is empty try to fill it 02429 fill(); 02430 02431 // try again 02432 while(size_>0){ 02433 head_ = (head_ + N -1) % N; 02434 size_--; 02435 Vertex v = vals_[head_]; 02436 if(aggregates_[v]==AggregatesMap<Vertex>::UNAGGREGATED) 02437 return v; 02438 } 02439 return NullEntry; 02440 } 02441 02442 #endif // DOXYGEN 02443 02444 template<class V> 02445 void printAggregates2d(const AggregatesMap<V>& aggregates, int n, int m, std::ostream& os) 02446 { 02447 std::ios_base::fmtflags oldOpts=os.flags(); 02448 02449 os.setf(std::ios_base::right, std::ios_base::adjustfield); 02450 02451 V max=0; 02452 int width=1; 02453 02454 for(int i=0; i< n*m; i++) 02455 max=std::max(max, aggregates[i]); 02456 02457 for(int i=10; i < 1000000; i*=10) 02458 if(max/i>0) 02459 width++; 02460 else 02461 break; 02462 02463 for(int j=0, entry=0; j < m; j++){ 02464 for(int i=0; i<n; i++, entry++){ 02465 os.width(width); 02466 os<<aggregates[entry]<<" "; 02467 } 02468 02469 os<<std::endl; 02470 } 02471 os<<std::endl; 02472 os.flags(oldOpts); 02473 } 02474 02475 02476 }// namespace Amg 02477 02478 }// namespace Dune 02479 02480 02481 #endif