dune-istl  2.2.0
aggregates.hh
Go to the documentation of this file.
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