Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 09:37:16

0001 // Copyright (C) 2004-2006 The Trustees of Indiana University.
0002 
0003 // Use, modification and distribution is subject to the Boost Software
0004 // License, Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0005 // http://www.boost.org/LICENSE_1_0.txt)
0006 
0007 //  Authors: Brian Barrett
0008 //           Douglas Gregor
0009 //           Andrew Lumsdaine
0010 #ifndef BOOST_GRAPH_PARALLEL_CC_PS_HPP
0011 #define BOOST_GRAPH_PARALLEL_CC_PS_HPP
0012 
0013 #ifndef BOOST_GRAPH_USE_MPI
0014 #error "Parallel BGL files should not be included unless <boost/graph/use_mpi.hpp> has been included"
0015 #endif
0016 
0017 #include <boost/assert.hpp>
0018 #include <boost/property_map/property_map.hpp>
0019 #include <boost/property_map/parallel/parallel_property_maps.hpp>
0020 #include <boost/graph/parallel/algorithm.hpp>
0021 #include <boost/pending/indirect_cmp.hpp>
0022 #include <boost/graph/graph_traits.hpp>
0023 #include <boost/graph/overloading.hpp>
0024 #include <boost/graph/distributed/concepts.hpp>
0025 #include <boost/graph/parallel/properties.hpp>
0026 #include <boost/graph/parallel/process_group.hpp>
0027 #include <boost/optional.hpp>
0028 #include <algorithm>
0029 #include <vector>
0030 #include <queue>
0031 #include <limits>
0032 #include <map>
0033 #include <boost/graph/parallel/container_traits.hpp>
0034 #include <boost/graph/iteration_macros.hpp>
0035 
0036 
0037 // Connected components algorithm based on a parallel search.
0038 //
0039 // Every N nodes starts a parallel search from the first vertex in
0040 // their local vertex list during the first superstep (the other nodes
0041 // remain idle during the first superstep to reduce the number of
0042 // conflicts in numbering the components).  At each superstep, all new
0043 // component mappings from remote nodes are handled.  If there is no
0044 // work from remote updates, a new vertex is removed from the local
0045 // list and added to the work queue.
0046 //
0047 // Components are allocated from the component_value_allocator object,
0048 // which ensures that a given component number is unique in the
0049 // system, currently by using the rank and number of processes to
0050 // stride allocations.
0051 //
0052 // When two components are discovered to actually be the same
0053 // component, a mapping is created in the collisions object.  The
0054 // lower component number is prefered in the resolution, so component
0055 // numbering resolution is consistent.  After the search has exhausted
0056 // all vertices in the graph, the mapping is shared with all
0057 // processes, and they independently resolve the comonent mapping (so
0058 // O((N * NP) + (V * NP)) work, in O(N + V) time, where N is the
0059 // number of mappings and V is the number of local vertices).  This
0060 // phase can likely be significantly sped up if a clever algorithm for
0061 // the reduction can be found.
0062 namespace boost { namespace graph { namespace distributed {
0063   namespace cc_ps_detail {
0064     // Local object for allocating component numbers.  There are two
0065     // places this happens in the code, and I was getting sick of them
0066     // getting out of sync.  Components are not tightly packed in
0067     // numbering, but are numbered to ensure each rank has its own
0068     // independent sets of numberings.
0069     template<typename component_value_type>
0070     class component_value_allocator {
0071     public:
0072       component_value_allocator(int num, int size) :
0073         last(0), num(num), size(size)
0074       {
0075       }
0076 
0077       component_value_type allocate(void)
0078       {
0079         component_value_type ret = num + (last * size);
0080         last++;
0081         return ret;
0082       }
0083 
0084     private:
0085       component_value_type last;
0086       int num;
0087       int size;
0088     };
0089 
0090 
0091     // Map of the "collisions" between component names in the global
0092     // component mapping.  TO make cleanup easier, component numbers
0093     // are added, pointing to themselves, when a new component is
0094     // found.  In order to make the results deterministic, the lower
0095     // component number is always taken.  The resolver will drill
0096     // through the map until it finds a component entry that points to
0097     // itself as the next value, allowing some cleanup to happen at
0098     // update() time.  Attempts are also made to update the mapping
0099     // when new entries are created.
0100     //
0101     // Note that there's an assumption that the entire mapping is
0102     // shared during the end of the algorithm, but before component
0103     // name resolution.
0104     template<typename component_value_type>
0105     class collision_map {
0106     public:
0107       collision_map() : num_unique(0)
0108       {
0109       }
0110 
0111       // add new component mapping first time component is used.  Own
0112       // function only so that we can sanity check there isn't already
0113       // a mapping for that component number (which would be bad)
0114       void add(const component_value_type &a) 
0115       {
0116         BOOST_ASSERT(collisions.count(a) == 0);
0117         collisions[a] = a;
0118       }
0119 
0120       // add a mapping between component values saying they're the
0121       // same component
0122       void add(const component_value_type &a, const component_value_type &b)
0123       {
0124         component_value_type high, low, tmp;
0125         if (a > b) {
0126           high = a;
0127           low = b;
0128         } else {
0129           high = b;
0130           low = a;
0131         }
0132 
0133         if (collisions.count(high) != 0 && collisions[high] != low) {
0134           tmp = collisions[high];
0135           if (tmp > low) {
0136             collisions[tmp] = low;
0137             collisions[high] = low;
0138           } else {
0139             collisions[low] = tmp;
0140             collisions[high] = tmp;
0141           }
0142         } else {
0143           collisions[high] = low;
0144         }
0145 
0146       }
0147 
0148       // get the "real" component number for the given component.
0149       // Used to resolve mapping at end of run.
0150       component_value_type update(component_value_type a)
0151       {
0152         BOOST_ASSERT(num_unique > 0);
0153         BOOST_ASSERT(collisions.count(a) != 0);
0154         return collisions[a];
0155       }
0156 
0157       // collapse the collisions tree, so that update is a one lookup
0158       // operation.  Count unique components at the same time.
0159       void uniqify(void)
0160       {
0161         typename std::map<component_value_type, component_value_type>::iterator i, end;
0162 
0163         end = collisions.end();
0164         for (i = collisions.begin() ; i != end ; ++i) {
0165           if (i->first == i->second) {
0166             num_unique++;
0167           } else {
0168             i->second = collisions[i->second];
0169           }
0170         }
0171       }
0172 
0173       // get the number of component entries that have an associated
0174       // component number of themselves, which are the real components
0175       // used in the final mapping.  This is the number of unique
0176       // components in the graph.
0177       int unique(void)
0178       {
0179         BOOST_ASSERT(num_unique > 0);
0180         return num_unique;
0181       }
0182 
0183       // "serialize" into a vector for communication.
0184       std::vector<component_value_type> serialize(void)
0185       {
0186         std::vector<component_value_type> ret;
0187         typename std::map<component_value_type, component_value_type>::iterator i, end;
0188 
0189         end = collisions.end();
0190         for (i = collisions.begin() ; i != end ; ++i) {
0191           ret.push_back(i->first);
0192           ret.push_back(i->second);
0193         }
0194 
0195         return ret;
0196       }
0197 
0198     private:
0199       std::map<component_value_type, component_value_type> collisions;
0200       int num_unique;
0201     };
0202 
0203 
0204     // resolver to handle remote updates.  The resolver will add
0205     // entries into the collisions map if required, and if it is the
0206     // first time the vertex has been touched, it will add the vertex
0207     // to the remote queue.  Note that local updates are handled
0208     // differently, in the main loop (below).
0209 
0210       // BWB - FIX ME - don't need graph anymore - can pull from key value of Component Map.
0211     template<typename ComponentMap, typename work_queue>
0212     struct update_reducer {
0213       BOOST_STATIC_CONSTANT(bool, non_default_resolver = false);
0214 
0215       typedef typename property_traits<ComponentMap>::value_type component_value_type;
0216       typedef typename property_traits<ComponentMap>::key_type vertex_descriptor;
0217 
0218       update_reducer(work_queue *q,
0219                      cc_ps_detail::collision_map<component_value_type> *collisions, 
0220                      processor_id_type pg_id) :
0221         q(q), collisions(collisions), pg_id(pg_id)
0222       {
0223       }
0224 
0225       // ghost cell initialization routine.  This should never be
0226       // called in this imlementation.
0227       template<typename K>
0228       component_value_type operator()(const K&) const
0229       { 
0230         return component_value_type(0); 
0231       }
0232 
0233       // resolver for remote updates.  I'm not entirely sure why, but
0234       // I decided to not change the value of the vertex if it's
0235       // already non-infinite.  It doesn't matter in the end, as we'll
0236       // touch every vertex in the cleanup phase anyway.  If the
0237       // component is currently infinite, set to the new component
0238       // number and add the vertex to the work queue.  If it's not
0239       // infinite, we've touched it already so don't add it to the
0240       // work queue.  Do add a collision entry so that we know the two
0241       // components are the same.
0242       component_value_type operator()(const vertex_descriptor &v,
0243                                       const component_value_type& current,
0244                                       const component_value_type& update) const
0245       {
0246         const component_value_type max = (std::numeric_limits<component_value_type>::max)();
0247         component_value_type ret = current;
0248 
0249         if (max == current) {
0250           q->push(v);
0251           ret = update;
0252         } else if (current != update) {
0253           collisions->add(current, update);
0254         }
0255 
0256         return ret;
0257       }                                    
0258 
0259       // So for whatever reason, the property map can in theory call
0260       // the resolver with a local descriptor in addition to the
0261       // standard global descriptor.  As far as I can tell, this code
0262       // path is never taken in this implementation, but I need to
0263       // have this code here to make it compile.  We just make a
0264       // global descriptor and call the "real" operator().
0265       template<typename K>
0266       component_value_type operator()(const K& v, 
0267                                       const component_value_type& current, 
0268                                       const component_value_type& update) const
0269       {
0270           return (*this)(vertex_descriptor(pg_id, v), current, update);
0271       }
0272 
0273     private:
0274       work_queue *q;
0275       collision_map<component_value_type> *collisions;
0276       boost::processor_id_type pg_id;
0277     };
0278 
0279   } // namespace cc_ps_detail
0280 
0281 
0282   template<typename Graph, typename ComponentMap>
0283   typename property_traits<ComponentMap>::value_type
0284   connected_components_ps(const Graph& g, ComponentMap c)
0285   {
0286     using boost::graph::parallel::process_group;
0287 
0288     typedef typename property_traits<ComponentMap>::value_type component_value_type;
0289     typedef typename graph_traits<Graph>::vertex_iterator vertex_iterator;
0290     typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;
0291     typedef typename boost::graph::parallel::process_group_type<Graph>
0292       ::type process_group_type;
0293     typedef typename process_group_type::process_id_type process_id_type;
0294     typedef std::queue<vertex_descriptor> work_queue;
0295 
0296     static const component_value_type max_component = 
0297       (std::numeric_limits<component_value_type>::max)();
0298     typename property_map<Graph, vertex_owner_t>::const_type
0299       owner = get(vertex_owner, g);
0300 
0301     // standard who am i? stuff
0302     process_group_type pg = process_group(g);
0303     process_id_type id = process_id(pg);
0304 
0305     // Initialize every vertex to have infinite component number
0306     BGL_FORALL_VERTICES_T(v, g, Graph) put(c, v, max_component);
0307 
0308     vertex_iterator current, end;
0309     boost::tie(current, end) = vertices(g);
0310 
0311     cc_ps_detail::component_value_allocator<component_value_type> cva(process_id(pg), num_processes(pg));
0312     cc_ps_detail::collision_map<component_value_type> collisions;
0313     work_queue q;  // this is intentionally a local data structure
0314     c.set_reduce(cc_ps_detail::update_reducer<ComponentMap, work_queue>(&q, &collisions, id));
0315 
0316     // add starting work
0317     while (true) {
0318         bool useful_found = false;
0319         component_value_type val = cva.allocate();
0320         put(c, *current, val);
0321         collisions.add(val);
0322         q.push(*current);
0323         if (0 != out_degree(*current, g)) useful_found = true;
0324         ++current;
0325         if (useful_found) break;
0326     }
0327 
0328     // Run the loop until everyone in the system is done
0329     bool global_done = false;
0330     while (!global_done) {
0331 
0332       // drain queue of work for this superstep
0333       while (!q.empty()) {
0334         vertex_descriptor v = q.front();
0335         q.pop();
0336         // iterate through outedges of the vertex currently being
0337         // examined, setting their component to our component.  There
0338         // is no way to end up in the queue without having a component
0339         // number already.
0340 
0341         BGL_FORALL_ADJ_T(v, peer, g, Graph) {
0342           component_value_type my_component = get(c, v);
0343 
0344           // update other vertex with our component information.
0345           // Resolver will handle remote collisions as well as whether
0346           // to put the vertex on the work queue or not.  We have to
0347           // handle local collisions and work queue management
0348           if (id == get(owner, peer)) {
0349             if (max_component == get(c, peer)) {
0350               put(c, peer, my_component);
0351               q.push(peer);
0352             } else if (my_component != get(c, peer)) {
0353               collisions.add(my_component, get(c, peer));
0354             }
0355           } else {
0356             put(c, peer, my_component);
0357           }
0358         }
0359       }
0360 
0361       // synchronize / start a new superstep.
0362       synchronize(pg);
0363       global_done = all_reduce(pg, (q.empty() && (current == end)), boost::parallel::minimum<bool>());
0364 
0365       // If the queue is currently empty, add something to do to start
0366       // the current superstep (supersteps start at the sync, not at
0367       // the top of the while loop as one might expect).  Down at the
0368       // bottom of the while loop so that not everyone starts the
0369       // algorithm with something to do, to try to reduce component
0370       // name conflicts
0371       if (q.empty()) {
0372         bool useful_found = false;
0373         for ( ; current != end && !useful_found ; ++current) {
0374           if (max_component == get(c, *current)) {
0375             component_value_type val = cva.allocate();
0376             put(c, *current, val);
0377             collisions.add(val);
0378             q.push(*current);
0379             if (0 != out_degree(*current, g)) useful_found = true;
0380           }
0381         }
0382       }
0383     }
0384 
0385     // share component mappings
0386     std::vector<component_value_type> global;
0387     std::vector<component_value_type> mine = collisions.serialize();
0388     all_gather(pg, mine.begin(), mine.end(), global);
0389     for (size_t i = 0 ; i < global.size() ; i += 2) {
0390       collisions.add(global[i], global[i + 1]);
0391     }
0392     collisions.uniqify();
0393 
0394     // update the component mappings
0395     BGL_FORALL_VERTICES_T(v, g, Graph) {
0396       put(c, v, collisions.update(get(c, v)));
0397     }
0398 
0399     return collisions.unique();
0400   }
0401 
0402 } // end namespace distributed
0403 
0404 } // end namespace graph
0405 
0406 } // end namespace boost
0407 
0408 #endif // BOOST_GRAPH_PARALLEL_CC_HPP