dune-istl
2.2.0
|
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