1 // Copyright (C) 2004-2006 The Trustees of Indiana University.
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)
7 // Authors: Brian Barrett
10 #ifndef BOOST_GRAPH_PARALLEL_CC_PS_HPP
11 #define BOOST_GRAPH_PARALLEL_CC_PS_HPP
13 #ifndef BOOST_GRAPH_USE_MPI
14 #error "Parallel BGL files should not be included unless <boost/graph/use_mpi.hpp> has been included"
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>
31 #include <boost/graph/parallel/container_traits.hpp>
32 #include <boost/graph/iteration_macros.hpp>
35 // Connected components algorithm based on a parallel search.
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.
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.
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 {
70 component_value_allocator(int num, int size) :
71 last(0), num(num), size(size)
75 component_value_type allocate(void)
77 component_value_type ret = num + (last * size);
83 component_value_type last;
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.
99 // Note that there's an assumption that the entire mapping is
100 // shared during the end of the algorithm, but before component
102 template<typename component_value_type>
103 class collision_map {
105 collision_map() : num_unique(0)
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)
114 assert(collisions.count(a) == 0);
118 // add a mapping between component values saying they're the
120 void add(const component_value_type &a, const component_value_type &b)
122 component_value_type high, low, tmp;
131 if (collisions.count(high) != 0 && collisions[high] != low) {
132 tmp = collisions[high];
134 collisions[tmp] = low;
135 collisions[high] = low;
137 collisions[low] = tmp;
138 collisions[high] = tmp;
141 collisions[high] = low;
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)
150 assert(num_unique > 0);
151 assert(collisions.count(a) != 0);
152 return collisions[a];
155 // collapse the collisions tree, so that update is a one lookup
156 // operation. Count unique components at the same time.
159 typename std::map<component_value_type, component_value_type>::iterator i, end;
161 end = collisions.end();
162 for (i = collisions.begin() ; i != end ; ++i) {
163 if (i->first == i->second) {
166 i->second = collisions[i->second];
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.
177 assert(num_unique > 0);
181 // "serialize" into a vector for communication.
182 std::vector<component_value_type> serialize(void)
184 std::vector<component_value_type> ret;
185 typename std::map<component_value_type, component_value_type>::iterator i, end;
187 end = collisions.end();
188 for (i = collisions.begin() ; i != end ; ++i) {
189 ret.push_back(i->first);
190 ret.push_back(i->second);
197 std::map<component_value_type, component_value_type> collisions;
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).
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);
213 typedef typename property_traits<ComponentMap>::value_type component_value_type;
214 typedef typename property_traits<ComponentMap>::key_type vertex_descriptor;
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)
223 // ghost cell initialization routine. This should never be
224 // called in this imlementation.
226 component_value_type operator()(const K&) const
228 return component_value_type(0);
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
244 const component_value_type max = (std::numeric_limits<component_value_type>::max)();
245 component_value_type ret = current;
247 if (max == current) {
250 } else if (current != update) {
251 collisions->add(current, update);
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().
264 component_value_type operator()(const K& v,
265 const component_value_type& current,
266 const component_value_type& update) const
268 return (*this)(vertex_descriptor(pg_id, v), current, update);
273 collision_map<component_value_type> *collisions;
274 boost::processor_id_type pg_id;
277 } // namespace cc_ps_detail
280 template<typename Graph, typename ComponentMap>
281 typename property_traits<ComponentMap>::value_type
282 connected_components_ps(const Graph& g, ComponentMap c)
284 using boost::graph::parallel::process_group;
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;
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);
301 // standard who am i? stuff
302 process_group_type pg = process_group(g);
303 process_id_type id = process_id(pg);
305 // Initialize every vertex to have infinite component number
306 BGL_FORALL_VERTICES_T(v, g, Graph) put(c, v, max_component);
308 vertex_iterator current, end;
309 boost::tie(current, end) = vertices(g);
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));
318 bool useful_found = false;
319 component_value_type val = cva.allocate();
320 put(c, *current, val);
323 if (0 != out_degree(*current, g)) useful_found = true;
325 if (useful_found) break;
328 // Run the loop until everyone in the system is done
329 bool global_done = false;
330 while (!global_done) {
332 // drain queue of work for this superstep
334 vertex_descriptor v = q.front();
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
341 BGL_FORALL_ADJ_T(v, peer, g, Graph) {
342 component_value_type my_component = get(c, v);
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);
352 } else if (my_component != get(c, peer)) {
353 collisions.add(my_component, get(c, peer));
356 put(c, peer, my_component);
361 // synchronize / start a new superstep.
363 global_done = all_reduce(pg, (q.empty() && (current == end)), boost::parallel::minimum<bool>());
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
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);
379 if (0 != out_degree(*current, g)) useful_found = true;
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]);
392 collisions.uniqify();
394 // update the component mappings
395 BGL_FORALL_VERTICES_T(v, g, Graph) {
396 put(c, v, collisions.update(get(c, v)));
399 return collisions.unique();
402 } // end namespace distributed
404 } // end namespace graph
406 } // end namespace boost
408 #endif // BOOST_GRAPH_PARALLEL_CC_HPP