]> https://gitweb.dealii.org/ - dealii-svn.git/blob
8e4d286c3d1074127854eff020a2b8f5288d5b4f
[dealii-svn.git] /
1 // Copyright (C) 2004-2006 The Trustees of Indiana University.
2
3 // Use, modification and distribution is subject to the Boost Software
4 // License, Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
5 // http://www.boost.org/LICENSE_1_0.txt)
6
7 //  Authors: Brian Barrett
8 //           Douglas Gregor
9 //           Andrew Lumsdaine
10 #ifndef BOOST_GRAPH_PARALLEL_CC_PS_HPP
11 #define BOOST_GRAPH_PARALLEL_CC_PS_HPP
12
13 #ifndef BOOST_GRAPH_USE_MPI
14 #error "Parallel BGL files should not be included unless <boost/graph/use_mpi.hpp> has been included"
15 #endif
16
17 #include <boost/property_map/property_map.hpp>
18 #include <boost/graph/parallel/algorithm.hpp>
19 #include <boost/pending/indirect_cmp.hpp>
20 #include <boost/graph/graph_traits.hpp>
21 #include <boost/graph/overloading.hpp>
22 #include <boost/graph/distributed/concepts.hpp>
23 #include <boost/graph/parallel/properties.hpp>
24 #include <boost/graph/parallel/process_group.hpp>
25 #include <boost/optional.hpp>
26 #include <algorithm>
27 #include <vector>
28 #include <queue>
29 #include <limits>
30 #include <map>
31 #include <boost/graph/parallel/container_traits.hpp>
32 #include <boost/graph/iteration_macros.hpp>
33
34
35 // Connected components algorithm based on a parallel search.
36 //
37 // Every N nodes starts a parallel search from the first vertex in
38 // their local vertex list during the first superstep (the other nodes
39 // remain idle during the first superstep to reduce the number of
40 // conflicts in numbering the components).  At each superstep, all new
41 // component mappings from remote nodes are handled.  If there is no
42 // work from remote updates, a new vertex is removed from the local
43 // list and added to the work queue.
44 //
45 // Components are allocated from the component_value_allocator object,
46 // which ensures that a given component number is unique in the
47 // system, currently by using the rank and number of processes to
48 // stride allocations.
49 //
50 // When two components are discovered to actually be the same
51 // component, a mapping is created in the collisions object.  The
52 // lower component number is prefered in the resolution, so component
53 // numbering resolution is consistent.  After the search has exhausted
54 // all vertices in the graph, the mapping is shared with all
55 // processes, and they independently resolve the comonent mapping (so
56 // O((N * NP) + (V * NP)) work, in O(N + V) time, where N is the
57 // number of mappings and V is the number of local vertices).  This
58 // phase can likely be significantly sped up if a clever algorithm for
59 // the reduction can be found.
60 namespace boost { namespace graph { namespace distributed {
61   namespace cc_ps_detail {
62     // Local object for allocating component numbers.  There are two
63     // places this happens in the code, and I was getting sick of them
64     // getting out of sync.  Components are not tightly packed in
65     // numbering, but are numbered to ensure each rank has its own
66     // independent sets of numberings.
67     template<typename component_value_type>
68     class component_value_allocator {
69     public:
70       component_value_allocator(int num, int size) :
71         last(0), num(num), size(size)
72       {
73       }
74
75       component_value_type allocate(void)
76       {
77         component_value_type ret = num + (last * size);
78         last++;
79         return ret;
80       }
81
82     private:
83       component_value_type last;
84       int num;
85       int size;
86     };
87
88
89     // Map of the "collisions" between component names in the global
90     // component mapping.  TO make cleanup easier, component numbers
91     // are added, pointing to themselves, when a new component is
92     // found.  In order to make the results deterministic, the lower
93     // component number is always taken.  The resolver will drill
94     // through the map until it finds a component entry that points to
95     // itself as the next value, allowing some cleanup to happen at
96     // update() time.  Attempts are also made to update the mapping
97     // when new entries are created.
98     //
99     // Note that there's an assumption that the entire mapping is
100     // shared during the end of the algorithm, but before component
101     // name resolution.
102     template<typename component_value_type>
103     class collision_map {
104     public:
105       collision_map() : num_unique(0)
106       {
107       }
108
109       // add new component mapping first time component is used.  Own
110       // function only so that we can sanity check there isn't already
111       // a mapping for that component number (which would be bad)
112       void add(const component_value_type &a) 
113       {
114         assert(collisions.count(a) == 0);
115         collisions[a] = a;
116       }
117
118       // add a mapping between component values saying they're the
119       // same component
120       void add(const component_value_type &a, const component_value_type &b)
121       {
122         component_value_type high, low, tmp;
123         if (a > b) {
124           high = a;
125           low = b;
126         } else {
127           high = b;
128           low = a;
129         }
130
131         if (collisions.count(high) != 0 && collisions[high] != low) {
132           tmp = collisions[high];
133           if (tmp > low) {
134             collisions[tmp] = low;
135             collisions[high] = low;
136           } else {
137             collisions[low] = tmp;
138             collisions[high] = tmp;
139           }
140         } else {
141           collisions[high] = low;
142         }
143
144       }
145
146       // get the "real" component number for the given component.
147       // Used to resolve mapping at end of run.
148       component_value_type update(component_value_type a)
149       {
150         assert(num_unique > 0);
151         assert(collisions.count(a) != 0);
152         return collisions[a];
153       }
154
155       // collapse the collisions tree, so that update is a one lookup
156       // operation.  Count unique components at the same time.
157       void uniqify(void)
158       {
159         typename std::map<component_value_type, component_value_type>::iterator i, end;
160
161         end = collisions.end();
162         for (i = collisions.begin() ; i != end ; ++i) {
163           if (i->first == i->second) {
164             num_unique++;
165           } else {
166             i->second = collisions[i->second];
167           }
168         }
169       }
170
171       // get the number of component entries that have an associated
172       // component number of themselves, which are the real components
173       // used in the final mapping.  This is the number of unique
174       // components in the graph.
175       int unique(void)
176       {
177         assert(num_unique > 0);
178         return num_unique;
179       }
180
181       // "serialize" into a vector for communication.
182       std::vector<component_value_type> serialize(void)
183       {
184         std::vector<component_value_type> ret;
185         typename std::map<component_value_type, component_value_type>::iterator i, end;
186
187         end = collisions.end();
188         for (i = collisions.begin() ; i != end ; ++i) {
189           ret.push_back(i->first);
190           ret.push_back(i->second);
191         }
192
193         return ret;
194       }
195
196     private:
197       std::map<component_value_type, component_value_type> collisions;
198       int num_unique;
199     };
200
201
202     // resolver to handle remote updates.  The resolver will add
203     // entries into the collisions map if required, and if it is the
204     // first time the vertex has been touched, it will add the vertex
205     // to the remote queue.  Note that local updates are handled
206     // differently, in the main loop (below).
207
208       // BWB - FIX ME - don't need graph anymore - can pull from key value of Component Map.
209     template<typename ComponentMap, typename work_queue>
210     struct update_reducer {
211       BOOST_STATIC_CONSTANT(bool, non_default_resolver = false);
212
213       typedef typename property_traits<ComponentMap>::value_type component_value_type;
214       typedef typename property_traits<ComponentMap>::key_type vertex_descriptor;
215
216       update_reducer(work_queue *q,
217                      cc_ps_detail::collision_map<component_value_type> *collisions, 
218                      processor_id_type pg_id) :
219         q(q), collisions(collisions), pg_id(pg_id)
220       {
221       }
222
223       // ghost cell initialization routine.  This should never be
224       // called in this imlementation.
225       template<typename K>
226       component_value_type operator()(const K&) const
227       { 
228         return component_value_type(0); 
229       }
230
231       // resolver for remote updates.  I'm not entirely sure why, but
232       // I decided to not change the value of the vertex if it's
233       // already non-infinite.  It doesn't matter in the end, as we'll
234       // touch every vertex in the cleanup phase anyway.  If the
235       // component is currently infinite, set to the new component
236       // number and add the vertex to the work queue.  If it's not
237       // infinite, we've touched it already so don't add it to the
238       // work queue.  Do add a collision entry so that we know the two
239       // components are the same.
240       component_value_type operator()(const vertex_descriptor &v,
241                                       const component_value_type& current,
242                                       const component_value_type& update) const
243       {
244         const component_value_type max = (std::numeric_limits<component_value_type>::max)();
245         component_value_type ret = current;
246
247         if (max == current) {
248           q->push(v);
249           ret = update;
250         } else if (current != update) {
251           collisions->add(current, update);
252         }
253
254         return ret;
255       }                                    
256
257       // So for whatever reason, the property map can in theory call
258       // the resolver with a local descriptor in addition to the
259       // standard global descriptor.  As far as I can tell, this code
260       // path is never taken in this implementation, but I need to
261       // have this code here to make it compile.  We just make a
262       // global descriptor and call the "real" operator().
263       template<typename K>
264       component_value_type operator()(const K& v, 
265                                       const component_value_type& current, 
266                                       const component_value_type& update) const
267       {
268           return (*this)(vertex_descriptor(pg_id, v), current, update);
269       }
270
271     private:
272       work_queue *q;
273       collision_map<component_value_type> *collisions;
274       boost::processor_id_type pg_id;
275     };
276
277   } // namespace cc_ps_detail
278
279
280   template<typename Graph, typename ComponentMap>
281   typename property_traits<ComponentMap>::value_type
282   connected_components_ps(const Graph& g, ComponentMap c)
283   {
284     using boost::graph::parallel::process_group;
285
286     typedef typename property_traits<ComponentMap>::value_type component_value_type;
287     typedef typename graph_traits<Graph>::vertex_iterator vertex_iterator;
288     typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;
289     typedef typename boost::graph::parallel::process_group_type<Graph>
290       ::type process_group_type;
291     typedef typename process_group_type::process_id_type process_id_type;
292     typedef typename property_map<Graph, vertex_owner_t>
293       ::const_type vertex_owner_map;
294     typedef std::queue<vertex_descriptor> work_queue;
295
296     static const component_value_type max_component = 
297       (std::numeric_limits<component_value_type>::max)();
298     typename property_map<Graph, vertex_owner_t>::const_type
299       owner = get(vertex_owner, g);
300
301     // standard who am i? stuff
302     process_group_type pg = process_group(g);
303     process_id_type id = process_id(pg);
304
305     // Initialize every vertex to have infinite component number
306     BGL_FORALL_VERTICES_T(v, g, Graph) put(c, v, max_component);
307
308     vertex_iterator current, end;
309     boost::tie(current, end) = vertices(g);
310
311     cc_ps_detail::component_value_allocator<component_value_type> cva(process_id(pg), num_processes(pg));
312     cc_ps_detail::collision_map<component_value_type> collisions;
313     work_queue q;  // this is intentionally a local data structure
314     c.set_reduce(cc_ps_detail::update_reducer<ComponentMap, work_queue>(&q, &collisions, id));
315
316     // add starting work
317     while (true) {
318         bool useful_found = false;
319         component_value_type val = cva.allocate();
320         put(c, *current, val);
321         collisions.add(val);
322         q.push(*current);
323         if (0 != out_degree(*current, g)) useful_found = true;
324         ++current;
325         if (useful_found) break;
326     }
327
328     // Run the loop until everyone in the system is done
329     bool global_done = false;
330     while (!global_done) {
331
332       // drain queue of work for this superstep
333       while (!q.empty()) {
334         vertex_descriptor v = q.front();
335         q.pop();
336         // iterate through outedges of the vertex currently being
337         // examined, setting their component to our component.  There
338         // is no way to end up in the queue without having a component
339         // number already.
340
341         BGL_FORALL_ADJ_T(v, peer, g, Graph) {
342           component_value_type my_component = get(c, v);
343
344           // update other vertex with our component information.
345           // Resolver will handle remote collisions as well as whether
346           // to put the vertex on the work queue or not.  We have to
347           // handle local collisions and work queue management
348           if (id == get(owner, peer)) {
349             if (max_component == get(c, peer)) {
350               put(c, peer, my_component);
351               q.push(peer);
352             } else if (my_component != get(c, peer)) {
353               collisions.add(my_component, get(c, peer));
354             }
355           } else {
356             put(c, peer, my_component);
357           }
358         }
359       }
360
361       // synchronize / start a new superstep.
362       synchronize(pg);
363       global_done = all_reduce(pg, (q.empty() && (current == end)), boost::parallel::minimum<bool>());
364
365       // If the queue is currently empty, add something to do to start
366       // the current superstep (supersteps start at the sync, not at
367       // the top of the while loop as one might expect).  Down at the
368       // bottom of the while loop so that not everyone starts the
369       // algorithm with something to do, to try to reduce component
370       // name conflicts
371       if (q.empty()) {
372         bool useful_found = false;
373         for ( ; current != end && !useful_found ; ++current) {
374           if (max_component == get(c, *current)) {
375             component_value_type val = cva.allocate();
376             put(c, *current, val);
377             collisions.add(val);
378             q.push(*current);
379             if (0 != out_degree(*current, g)) useful_found = true;
380           }
381         }
382       }
383     }
384
385     // share component mappings
386     std::vector<component_value_type> global;
387     std::vector<component_value_type> mine = collisions.serialize();
388     all_gather(pg, mine.begin(), mine.end(), global);
389     for (size_t i = 0 ; i < global.size() ; i += 2) {
390       collisions.add(global[i], global[i + 1]);
391     }
392     collisions.uniqify();
393
394     // update the component mappings
395     BGL_FORALL_VERTICES_T(v, g, Graph) {
396       put(c, v, collisions.update(get(c, v)));
397     }
398
399     return collisions.unique();
400   }
401
402 } // end namespace distributed
403
404 } // end namespace graph
405
406 } // end namespace boost
407
408 #endif // BOOST_GRAPH_PARALLEL_CC_HPP

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.