dune-istl  2.2.0
indicescoarsener.hh
Go to the documentation of this file.
00001 // $Id: indicescoarsener.hh 1444 2011-01-23 17:24:33Z rebecca $
00002 #ifndef DUNE_AMG_INDICESCOARSENER_HH
00003 #define DUNE_AMG_INDICESCOARSENER_HH
00004 
00005 #include<dune/common/parallel/indicessyncer.hh>
00006 #include<vector>
00007 #include"renumberer.hh"
00008 
00009 #if HAVE_MPI
00010 #include<dune/istl/owneroverlapcopy.hh>
00011 #endif
00012 
00013 #include"pinfo.hh"
00014 
00015 namespace Dune
00016 {
00017   namespace Amg
00018   {
00019     
00031     template<typename T, typename E>
00032     class IndicesCoarsener
00033     {
00034     };
00035 
00036     
00037 #if HAVE_MPI
00038 
00039     template<typename T, typename E>
00040     class ParallelIndicesCoarsener
00041     {
00042     public:
00046       typedef E ExcludedAttributes;
00047       
00051       typedef T ParallelInformation;
00052       
00053       typedef typename ParallelInformation::ParallelIndexSet ParallelIndexSet;
00054       
00058       typedef typename ParallelIndexSet::GlobalIndex GlobalIndex;
00059       
00063       typedef typename ParallelIndexSet::LocalIndex LocalIndex;
00064 
00068       typedef typename LocalIndex::Attribute Attribute;
00069       
00073       typedef Dune::RemoteIndices<ParallelIndexSet> RemoteIndices;
00074       
00086       template<typename Graph, typename VM>
00087       static typename Graph::VertexDescriptor
00088       coarsen(ParallelInformation& fineInfo,
00089               Graph& fineGraph,
00090               VM& visitedMap,
00091               AggregatesMap<typename Graph::VertexDescriptor>& aggregates,
00092               ParallelInformation& coarseInfo);
00093       
00094     private:
00095       template<typename G, typename I>
00096       class ParallelAggregateRenumberer: public AggregateRenumberer<G>
00097       {
00098         typedef typename G::VertexDescriptor Vertex;
00099 
00100         typedef I GlobalLookupIndexSet;
00101         
00102         typedef typename GlobalLookupIndexSet::IndexPair IndexPair;
00103         
00104         typedef typename IndexPair::GlobalIndex GlobalIndex;
00105         
00106       public:
00107         ParallelAggregateRenumberer(AggregatesMap<Vertex>& aggregates, const I& lookup)
00108           :  AggregateRenumberer<G>(aggregates),  isPublic_(false), lookup_(lookup), 
00109              globalIndex_(std::numeric_limits<GlobalIndex>::max())
00110         {}
00111                 
00112         
00113         void operator()(const typename G::ConstEdgeIterator& edge)
00114         {
00115           AggregateRenumberer<G>::operator()(edge);
00116           const IndexPair* pair= lookup_.pair(edge.target());
00117           if(pair!=0){
00118             globalIndex(pair->global());
00119             attribute(pair->local().attribute());
00120             isPublic(pair->local().isPublic());
00121           }
00122         }
00123         
00124         Vertex operator()(const GlobalIndex& global)
00125         {
00126           Vertex current = this->number_;
00127           this->operator++();
00128           return current;
00129         }
00130         
00131         bool isPublic()
00132         {
00133           return isPublic_;
00134         }
00135         
00136         void isPublic(bool b)
00137         {
00138           isPublic_ = isPublic_ || b;
00139         }
00140         
00141         void reset()
00142         {
00143           globalIndex_ = std::numeric_limits<GlobalIndex>::max();
00144           isPublic_=false;
00145         }
00146 
00147         void attribute(const Attribute& attribute)
00148         {
00149           attribute_=attribute;
00150         }
00151         
00152         Attribute attribute()
00153         {
00154           return attribute_;
00155         }
00156         
00157         const GlobalIndex& globalIndex() const
00158         {
00159           return globalIndex_;
00160         }
00161         
00162         void globalIndex(const GlobalIndex& global)
00163         {
00164           globalIndex_ = global;
00165         }
00166         
00167       private:
00168         bool isPublic_;
00169         Attribute attribute_;
00170         const GlobalLookupIndexSet& lookup_;
00171         GlobalIndex globalIndex_;
00172       };
00173       
00174     template<typename Graph, typename VM, typename I>
00175     static void buildCoarseIndexSet(const ParallelInformation& pinfo,
00176                                     Graph& fineGraph,
00177                                     VM& visitedMap,
00178                                     AggregatesMap<typename Graph::VertexDescriptor>& aggregates,
00179                                     ParallelIndexSet& coarseIndices,
00180                                     ParallelAggregateRenumberer<Graph,I>& renumberer);
00181     
00182       template<typename Graph,typename I>
00183       static void buildCoarseRemoteIndices(const RemoteIndices& fineRemote,
00184                                            const AggregatesMap<typename Graph::VertexDescriptor>& aggregates,
00185                                            ParallelIndexSet& coarseIndices,
00186                                            RemoteIndices& coarseRemote,
00187                                            ParallelAggregateRenumberer<Graph,I>& renumberer);
00188             
00189     };
00190 
00194     template<typename G, typename L, typename E>
00195     class IndicesCoarsener<OwnerOverlapCopyCommunication<G,L>,E>
00196       : public ParallelIndicesCoarsener<OwnerOverlapCopyCommunication<G,L>,E>
00197     {};
00198     
00199 
00200 #endif    
00201 
00208     template<typename E>
00209     class IndicesCoarsener<SequentialInformation,E>
00210     {
00211     public:
00212       template<typename Graph, typename VM>
00213       static typename Graph::VertexDescriptor 
00214       coarsen(const SequentialInformation& fineInfo,
00215               Graph& fineGraph,
00216               VM& visitedMap,
00217               AggregatesMap<typename Graph::VertexDescriptor>& aggregates,
00218               SequentialInformation& coarseInfo);
00219     };
00220 
00221 #if HAVE_MPI
00222     template<typename T, typename E>
00223     template<typename Graph, typename VM>
00224     inline typename Graph::VertexDescriptor 
00225     ParallelIndicesCoarsener<T,E>::coarsen(ParallelInformation& fineInfo,
00226                                    Graph& fineGraph,
00227                                    VM& visitedMap,
00228                                    AggregatesMap<typename Graph::VertexDescriptor>& aggregates,
00229                                    ParallelInformation& coarseInfo)
00230     {
00231       ParallelAggregateRenumberer<Graph,typename ParallelInformation::GlobalLookupIndexSet> renumberer(aggregates, fineInfo.globalLookup());
00232       buildCoarseIndexSet(fineInfo, fineGraph, visitedMap, aggregates, 
00233                           coarseInfo.indexSet(), renumberer);
00234       buildCoarseRemoteIndices(fineInfo.remoteIndices(), aggregates, coarseInfo.indexSet(), 
00235                                coarseInfo.remoteIndices(), renumberer);
00236 
00237       return renumberer;
00238     }
00239     
00240     template<typename T, typename E>
00241     template<typename Graph, typename VM, typename I>
00242     void ParallelIndicesCoarsener<T,E>::buildCoarseIndexSet(const ParallelInformation& pinfo,
00243                                                     Graph& fineGraph,
00244                                                     VM& visitedMap,
00245                                                     AggregatesMap<typename Graph::VertexDescriptor>& aggregates,
00246                                                     ParallelIndexSet& coarseIndices,
00247                                                     ParallelAggregateRenumberer<Graph,I>& renumberer)
00248     {
00249       // fineGraph is the local subgraph corresponding to the vertices the process owns.
00250       // i.e. no overlap/copy vertices can be visited traversing the graph
00251       typedef typename Graph::VertexDescriptor Vertex;       
00252       typedef typename Graph::ConstVertexIterator Iterator;
00253       typedef typename ParallelInformation::GlobalLookupIndexSet GlobalLookupIndexSet;
00254       
00255       Iterator end = fineGraph.end();
00256       const GlobalLookupIndexSet& lookup = pinfo.globalLookup();
00257       
00258       coarseIndices.beginResize();
00259 
00260       // Setup the coarse index set and renumber the aggregate consecutively
00261       // ascending from zero according to the minimum global index belonging
00262       // to the aggregate
00263       for(Iterator index = fineGraph.begin(); index != end; ++index){
00264         if(aggregates[*index]!=AggregatesMap<typename Graph::VertexDescriptor>::ISOLATED)
00265           // Isolated vertices will not be represented on the next level.
00266           // These should only be there if skipIsolated is activiated in
00267           // the coarsening criterion as otherwise they will be aggregated
00268           // and should have real aggregate number in the map right now.
00269           if(!get(visitedMap, *index)){
00270             // This vertex was not visited by breadthFirstSearch yet.
00271             typedef typename GlobalLookupIndexSet::IndexPair IndexPair;   
00272             const IndexPair* pair= lookup.pair(*index);
00273               
00274             renumberer.reset(); // reset attribute and global index.
00275             if(pair!=0){
00276               // vertex is in the index set. Note that not all vertices have
00277               // to be in the index set, just the ones where communication
00278               // will happen.
00279               assert(!ExcludedAttributes::contains(pair->local().attribute()));
00280               renumberer.attribute(pair->local().attribute());
00281               renumberer.isPublic(pair->local().isPublic());
00282               renumberer.globalIndex(pair->global());
00283             }
00284              
00285             // Reconstruct aggregate and mark vertices as visited
00286             aggregates.template breadthFirstSearch<false>(*index, aggregates[*index], 
00287                                                             fineGraph, renumberer, visitedMap);
00288             
00289             typedef typename GlobalLookupIndexSet::IndexPair::GlobalIndex GlobalIndex;
00290             
00291             if(renumberer.globalIndex()!=std::numeric_limits<GlobalIndex>::max()){
00292               // vertex is in the index set.
00293               //std::cout <<" Adding global="<< renumberer.globalIndex()<<" local="<<static_cast<std::size_t>(renumberer)<<std::endl;
00294               coarseIndices.add(renumberer.globalIndex(), 
00295                                 LocalIndex(renumberer, renumberer.attribute(), 
00296                                            renumberer.isPublic()));
00297             }
00298             
00299             aggregates[*index] = renumberer;
00300             ++renumberer;
00301           }
00302       }
00303 
00304       coarseIndices.endResize();
00305 
00306       assert(static_cast<std::size_t>(renumberer) >= coarseIndices.size());
00307       
00308       // Reset the visited flags      
00309       for(Iterator vertex=fineGraph.begin(); vertex != end; ++vertex)
00310         put(visitedMap, *vertex, false);      
00311     }
00312     
00313     template<typename T, typename E>
00314     template<typename Graph, typename I>
00315     void ParallelIndicesCoarsener<T,E>::buildCoarseRemoteIndices(const RemoteIndices& fineRemote,
00316                                                          const AggregatesMap<typename Graph::VertexDescriptor>& aggregates,
00317                                                          ParallelIndexSet& coarseIndices,
00318                                                          RemoteIndices& coarseRemote,
00319                                                          ParallelAggregateRenumberer<Graph,I>& renumberer)
00320     {
00321       std::vector<char> attributes(static_cast<std::size_t>(renumberer));
00322 
00323       GlobalLookupIndexSet<ParallelIndexSet> coarseLookup(coarseIndices, static_cast<std::size_t>(renumberer));
00324       
00325       typedef typename RemoteIndices::const_iterator Iterator;
00326       Iterator end = fineRemote.end();
00327       
00328       for(Iterator neighbour = fineRemote.begin();
00329           neighbour != end; ++neighbour){
00330         int process = neighbour->first;
00331         
00332         assert(neighbour->second.first==neighbour->second.second);
00333         
00334         // Mark all as not known
00335         typedef typename std::vector<char>::iterator CIterator;
00336 
00337         for(CIterator iter=attributes.begin(); iter!= attributes.end(); ++iter)
00338           *iter = std::numeric_limits<char>::max();
00339         
00340         typedef typename RemoteIndices::RemoteIndexList::const_iterator Iterator;
00341         Iterator riEnd = neighbour->second.second->end();
00342 
00343         for(Iterator index = neighbour->second.second->begin();
00344             index != riEnd; ++index){
00345           if(!E::contains(index->localIndexPair().local().attribute()) && 
00346              aggregates[index->localIndexPair().local()] != 
00347              AggregatesMap<typename Graph::VertexDescriptor>::ISOLATED)
00348           {
00349             assert(aggregates[index->localIndexPair().local()]<attributes.size());
00350             if (attributes[aggregates[index->localIndexPair().local()]] != 3)
00351               attributes[aggregates[index->localIndexPair().local()]] = index->attribute();
00352           }  
00353         }
00354         
00355         // Build remote index list
00356         typedef RemoteIndexListModifier<ParallelIndexSet,typename RemoteIndices::Allocator,false> Modifier;
00357         typedef typename RemoteIndices::RemoteIndex RemoteIndex;
00358         typedef typename ParallelIndexSet::const_iterator IndexIterator;
00359 
00360         Modifier coarseList = coarseRemote.template getModifier<false,true>(process);
00361         
00362         IndexIterator iend = coarseIndices.end();
00363         for(IndexIterator index = coarseIndices.begin(); index != iend; ++index)
00364           if(attributes[index->local()] != std::numeric_limits<char>::max()){
00365             // remote index is present
00366             coarseList.insert(RemoteIndex(Attribute(attributes[index->local()]), &(*index)));
00367           }
00368         //std::cout<<coarseRemote<<std::endl;
00369       }
00370       
00371       // The number of neighbours should not change!
00372       assert(coarseRemote.neighbours()==fineRemote.neighbours());
00373       
00374       // snyc the index set and the remote indices to recompute missing
00375       // indices
00376       IndicesSyncer<ParallelIndexSet> syncer(coarseIndices, coarseRemote);
00377       syncer.sync(renumberer);
00378       
00379     }
00380 
00381 #endif
00382 
00383     template<typename E>
00384     template<typename Graph, typename VM>
00385     typename Graph::VertexDescriptor 
00386     IndicesCoarsener<SequentialInformation,E>::coarsen(const SequentialInformation& fineInfo,
00387                                                        Graph& fineGraph,
00388                                                        VM& visitedMap,
00389                                                        AggregatesMap<typename Graph::VertexDescriptor>& aggregates,
00390                                                        SequentialInformation& coarseInfo)
00391     {
00392       typedef typename Graph::VertexDescriptor Vertex;
00393       AggregateRenumberer<Graph> renumberer(aggregates);
00394       typedef typename Graph::VertexIterator Iterator;
00395       
00396       for(Iterator vertex=fineGraph.begin(), endVertex=fineGraph.end();
00397           vertex != endVertex; ++vertex)
00398         if(aggregates[*vertex]!=AggregatesMap<Vertex>::ISOLATED && 
00399            !get(visitedMap, *vertex)){
00400 
00401           aggregates.template breadthFirstSearch<false>(*vertex, aggregates[*vertex], 
00402                                                         fineGraph, renumberer, visitedMap);
00403           aggregates[*vertex] = renumberer;
00404           ++renumberer;
00405         }
00406         
00407       for(Iterator vertex=fineGraph.begin(), endVertex=fineGraph.end();
00408           vertex != endVertex; ++vertex)
00409         put(visitedMap, *vertex, false);
00410 
00411       return renumberer;
00412     }
00413     
00414   } //namespace Amg
00415 } // namespace Dune
00416 #endif