dune-istl  2.2.1
indicescoarsener.hh
Go to the documentation of this file.
1 // $Id: indicescoarsener.hh 1871 2013-02-26 11:21:45Z mblatt $
2 #ifndef DUNE_AMG_INDICESCOARSENER_HH
3 #define DUNE_AMG_INDICESCOARSENER_HH
4 
5 #include<dune/common/parallel/indicessyncer.hh>
6 #include<vector>
7 #include"renumberer.hh"
8 
9 #if HAVE_MPI
11 #endif
12 
13 #include"pinfo.hh"
14 
15 namespace Dune
16 {
17  namespace Amg
18  {
19 
31  template<typename T, typename E>
33  {
34  };
35 
36 
37 #if HAVE_MPI
38 
39  template<typename T, typename E>
41  {
42  public:
46  typedef E ExcludedAttributes;
47 
52 
54 
58  typedef typename ParallelIndexSet::GlobalIndex GlobalIndex;
59 
63  typedef typename ParallelIndexSet::LocalIndex LocalIndex;
64 
68  typedef typename LocalIndex::Attribute Attribute;
69 
73  typedef Dune::RemoteIndices<ParallelIndexSet> RemoteIndices;
74 
86  template<typename Graph, typename VM>
87  static typename Graph::VertexDescriptor
88  coarsen(ParallelInformation& fineInfo,
89  Graph& fineGraph,
90  VM& visitedMap,
92  ParallelInformation& coarseInfo);
93 
94  private:
95  template<typename G, typename I>
96  class ParallelAggregateRenumberer: public AggregateRenumberer<G>
97  {
98  typedef typename G::VertexDescriptor Vertex;
99 
100  typedef I GlobalLookupIndexSet;
101 
102  typedef typename GlobalLookupIndexSet::IndexPair IndexPair;
103 
104  typedef typename IndexPair::GlobalIndex GlobalIndex;
105 
106  public:
107  ParallelAggregateRenumberer(AggregatesMap<Vertex>& aggregates, const I& lookup)
108  : AggregateRenumberer<G>(aggregates), isPublic_(false), lookup_(lookup),
109  globalIndex_(std::numeric_limits<GlobalIndex>::max())
110  {}
111 
112 
113  void operator()(const typename G::ConstEdgeIterator& edge)
114  {
116  const IndexPair* pair= lookup_.pair(edge.target());
117  if(pair!=0){
118  globalIndex(pair->global());
119  attribute(pair->local().attribute());
120  isPublic(pair->local().isPublic());
121  }
122  }
123 
124  Vertex operator()(const GlobalIndex& global)
125  {
126  Vertex current = this->number_;
127  this->operator++();
128  return current;
129  }
130 
131  bool isPublic()
132  {
133  return isPublic_;
134  }
135 
136  void isPublic(bool b)
137  {
138  isPublic_ = isPublic_ || b;
139  }
140 
141  void reset()
142  {
143  globalIndex_ = std::numeric_limits<GlobalIndex>::max();
144  isPublic_=false;
145  }
146 
147  void attribute(const Attribute& attribute)
148  {
149  attribute_=attribute;
150  }
151 
152  Attribute attribute()
153  {
154  return attribute_;
155  }
156 
157  const GlobalIndex& globalIndex() const
158  {
159  return globalIndex_;
160  }
161 
162  void globalIndex(const GlobalIndex& global)
163  {
164  globalIndex_ = global;
165  }
166 
167  private:
168  bool isPublic_;
169  Attribute attribute_;
170  const GlobalLookupIndexSet& lookup_;
171  GlobalIndex globalIndex_;
172  };
173 
174  template<typename Graph, typename VM, typename I>
175  static void buildCoarseIndexSet(const ParallelInformation& pinfo,
176  Graph& fineGraph,
177  VM& visitedMap,
179  ParallelIndexSet& coarseIndices,
180  ParallelAggregateRenumberer<Graph,I>& renumberer);
181 
182  template<typename Graph,typename I>
183  static void buildCoarseRemoteIndices(const RemoteIndices& fineRemote,
185  ParallelIndexSet& coarseIndices,
186  RemoteIndices& coarseRemote,
187  ParallelAggregateRenumberer<Graph,I>& renumberer);
188 
189  };
190 
194  template<typename G, typename L, typename E>
196  : public ParallelIndicesCoarsener<OwnerOverlapCopyCommunication<G,L>,E>
197  {};
198 
199 
200 #endif
201 
208  template<typename E>
210  {
211  public:
212  template<typename Graph, typename VM>
213  static typename Graph::VertexDescriptor
214  coarsen(const SequentialInformation& fineInfo,
215  Graph& fineGraph,
216  VM& visitedMap,
218  SequentialInformation& coarseInfo);
219  };
220 
221 #if HAVE_MPI
222  template<typename T, typename E>
223  template<typename Graph, typename VM>
224  inline typename Graph::VertexDescriptor
226  Graph& fineGraph,
227  VM& visitedMap,
229  ParallelInformation& coarseInfo)
230  {
231  ParallelAggregateRenumberer<Graph,typename ParallelInformation::GlobalLookupIndexSet> renumberer(aggregates, fineInfo.globalLookup());
232  buildCoarseIndexSet(fineInfo, fineGraph, visitedMap, aggregates,
233  coarseInfo.indexSet(), renumberer);
234  buildCoarseRemoteIndices(fineInfo.remoteIndices(), aggregates, coarseInfo.indexSet(),
235  coarseInfo.remoteIndices(), renumberer);
236 
237  return renumberer;
238  }
239 
240  template<typename T, typename E>
241  template<typename Graph, typename VM, typename I>
243  Graph& fineGraph,
244  VM& visitedMap,
246  ParallelIndexSet& coarseIndices,
247  ParallelAggregateRenumberer<Graph,I>& renumberer)
248  {
249  // fineGraph is the local subgraph corresponding to the vertices the process owns.
250  // i.e. no overlap/copy vertices can be visited traversing the graph
251  typedef typename Graph::ConstVertexIterator Iterator;
252  typedef typename ParallelInformation::GlobalLookupIndexSet GlobalLookupIndexSet;
253 
254  Iterator end = fineGraph.end();
255  const GlobalLookupIndexSet& lookup = pinfo.globalLookup();
256 
257  coarseIndices.beginResize();
258 
259  // Setup the coarse index set and renumber the aggregate consecutively
260  // ascending from zero according to the minimum global index belonging
261  // to the aggregate
262  for(Iterator index = fineGraph.begin(); index != end; ++index){
264  // Isolated vertices will not be represented on the next level.
265  // These should only be there if skipIsolated is activiated in
266  // the coarsening criterion as otherwise they will be aggregated
267  // and should have real aggregate number in the map right now.
268  if(!get(visitedMap, *index)){
269  // This vertex was not visited by breadthFirstSearch yet.
270  typedef typename GlobalLookupIndexSet::IndexPair IndexPair;
271  const IndexPair* pair= lookup.pair(*index);
272 
273  renumberer.reset(); // reset attribute and global index.
274  if(pair!=0){
275  // vertex is in the index set. Note that not all vertices have
276  // to be in the index set, just the ones where communication
277  // will happen.
278  assert(!ExcludedAttributes::contains(pair->local().attribute()));
279  renumberer.attribute(pair->local().attribute());
280  renumberer.isPublic(pair->local().isPublic());
281  renumberer.globalIndex(pair->global());
282  }
283 
284  // Reconstruct aggregate and mark vertices as visited
285  aggregates.template breadthFirstSearch<false>(*index, aggregates[*index],
286  fineGraph, renumberer, visitedMap);
287 
288  typedef typename GlobalLookupIndexSet::IndexPair::GlobalIndex GlobalIndex;
289 
290  if(renumberer.globalIndex()!=std::numeric_limits<GlobalIndex>::max()){
291  // vertex is in the index set.
292  //std::cout <<" Adding global="<< renumberer.globalIndex()<<" local="<<static_cast<std::size_t>(renumberer)<<std::endl;
293  coarseIndices.add(renumberer.globalIndex(),
294  LocalIndex(renumberer, renumberer.attribute(),
295  renumberer.isPublic()));
296  }
297 
298  aggregates[*index] = renumberer;
299  ++renumberer;
300  }
301  }
302 
303  coarseIndices.endResize();
304 
305  assert(static_cast<std::size_t>(renumberer) >= coarseIndices.size());
306 
307  // Reset the visited flags
308  for(Iterator vertex=fineGraph.begin(); vertex != end; ++vertex)
309  put(visitedMap, *vertex, false);
310  }
311 
312  template<typename T, typename E>
313  template<typename Graph, typename I>
316  ParallelIndexSet& coarseIndices,
317  RemoteIndices& coarseRemote,
318  ParallelAggregateRenumberer<Graph,I>& renumberer)
319  {
320  std::vector<char> attributes(static_cast<std::size_t>(renumberer));
321 
322  GlobalLookupIndexSet<ParallelIndexSet> coarseLookup(coarseIndices, static_cast<std::size_t>(renumberer));
323 
324  typedef typename RemoteIndices::const_iterator Iterator;
325  Iterator end = fineRemote.end();
326 
327  for(Iterator neighbour = fineRemote.begin();
328  neighbour != end; ++neighbour){
329  int process = neighbour->first;
330 
331  assert(neighbour->second.first==neighbour->second.second);
332 
333  // Mark all as not known
334  typedef typename std::vector<char>::iterator CIterator;
335 
336  for(CIterator iter=attributes.begin(); iter!= attributes.end(); ++iter)
337  *iter = std::numeric_limits<char>::max();
338 
339  typedef typename RemoteIndices::RemoteIndexList::const_iterator Iterator;
340  Iterator riEnd = neighbour->second.second->end();
341 
342  for(Iterator index = neighbour->second.second->begin();
343  index != riEnd; ++index){
344  if(!E::contains(index->localIndexPair().local().attribute()) &&
345  aggregates[index->localIndexPair().local()] !=
347  {
348  assert(aggregates[index->localIndexPair().local()]<attributes.size());
349  if (attributes[aggregates[index->localIndexPair().local()]] != 3)
350  attributes[aggregates[index->localIndexPair().local()]] = index->attribute();
351  }
352  }
353 
354  // Build remote index list
355  typedef RemoteIndexListModifier<ParallelIndexSet,typename RemoteIndices::Allocator,false> Modifier;
356  typedef typename RemoteIndices::RemoteIndex RemoteIndex;
357  typedef typename ParallelIndexSet::const_iterator IndexIterator;
358 
359  Modifier coarseList = coarseRemote.template getModifier<false,true>(process);
360 
361  IndexIterator iend = coarseIndices.end();
362  for(IndexIterator index = coarseIndices.begin(); index != iend; ++index)
363  if(attributes[index->local()] != std::numeric_limits<char>::max()){
364  // remote index is present
365  coarseList.insert(RemoteIndex(Attribute(attributes[index->local()]), &(*index)));
366  }
367  //std::cout<<coarseRemote<<std::endl;
368  }
369 
370  // The number of neighbours should not change!
371  assert(coarseRemote.neighbours()==fineRemote.neighbours());
372 
373  // snyc the index set and the remote indices to recompute missing
374  // indices
375  IndicesSyncer<ParallelIndexSet> syncer(coarseIndices, coarseRemote);
376  syncer.sync(renumberer);
377 
378  }
379 
380 #endif
381 
382  template<typename E>
383  template<typename Graph, typename VM>
384  typename Graph::VertexDescriptor
386  Graph& fineGraph,
387  VM& visitedMap,
389  SequentialInformation& coarseInfo)
390  {
391  typedef typename Graph::VertexDescriptor Vertex;
392  AggregateRenumberer<Graph> renumberer(aggregates);
393  typedef typename Graph::VertexIterator Iterator;
394 
395  for(Iterator vertex=fineGraph.begin(), endVertex=fineGraph.end();
396  vertex != endVertex; ++vertex)
397  if(aggregates[*vertex]!=AggregatesMap<Vertex>::ISOLATED &&
398  !get(visitedMap, *vertex)){
399 
400  aggregates.template breadthFirstSearch<false>(*vertex, aggregates[*vertex],
401  fineGraph, renumberer, visitedMap);
402  aggregates[*vertex] = renumberer;
403  ++renumberer;
404  }
405 
406  for(Iterator vertex=fineGraph.begin(), endVertex=fineGraph.end();
407  vertex != endVertex; ++vertex)
408  put(visitedMap, *vertex, false);
409 
410  return renumberer;
411  }
412 
413  } //namespace Amg
414 } // namespace Dune
415 #endif