boost/graph/distributed/hohberg_biconnected_components.hpp
// Copyright 2005 The Trustees of Indiana University. // Use, modification and distribution is subject to the Boost Software // License, Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at // http://www.boost.org/LICENSE_1_0.txt) // Authors: Douglas Gregor // Andrew Lumsdaine // An implementation of Walter Hohberg's distributed biconnected // components algorithm, from: // // Walter Hohberg. How to Find Biconnected Components in Distributed // Networks. J. Parallel Distrib. Comput., 9(4):374-386, 1990. // #ifndef BOOST_GRAPH_DISTRIBUTED_HOHBERG_BICONNECTED_COMPONENTS_HPP #define BOOST_GRAPH_DISTRIBUTED_HOHBERG_BICONNECTED_COMPONENTS_HPP #ifndef BOOST_GRAPH_USE_MPI #error "Parallel BGL files should not be included unless <boost/graph/use_mpi.hpp> has been included" #endif /* You can define PBGL_HOHBERG_DEBUG to an integer value (1, 2, or 3) * to enable debugging information. 1 includes only the phases of the * algorithm and messages as their are received. 2 and 3 add * additional levels of detail about internal data structures related * to the algorithm itself. * * #define PBGL_HOHBERG_DEBUG 1 */ #include <boost/graph/graph_traits.hpp> #include <boost/graph/parallel/container_traits.hpp> #include <boost/graph/parallel/process_group.hpp> #include <boost/static_assert.hpp> #include <boost/mpi/operations.hpp> #include <boost/type_traits/is_convertible.hpp> #include <boost/graph/graph_concepts.hpp> #include <boost/graph/iteration_macros.hpp> #include <boost/optional.hpp> #include <utility> // for std::pair #include <boost/assert.hpp> #include <algorithm> // for std::find, std::mismatch #include <vector> #include <boost/graph/parallel/algorithm.hpp> #include <boost/graph/distributed/connected_components.hpp> #include <boost/concept/assert.hpp> namespace boost { namespace graph { namespace distributed { namespace hohberg_detail { enum message_kind { /* A header for the PATH message, stating which edge the message is coming on and how many vertices will be following. The data structure communicated will be a path_header. */ msg_path_header, /* A message containing the vertices that make up a path. It will always follow a msg_path_header and will contain vertex descriptors, only. */ msg_path_vertices, /* A header for the TREE message, stating the value of gamma and the number of vertices to come in the following msg_tree_vertices. */ msg_tree_header, /* A message containing the vertices that make up the set of vertices in the same bicomponent as the sender. It will always follow a msg_tree_header and will contain vertex descriptors, only. */ msg_tree_vertices, /* Provides a name for the biconnected component of the edge. */ msg_name }; // Payload for a msg_path_header message. template<typename EdgeDescriptor> struct path_header { // The edge over which the path is being communicated EdgeDescriptor edge; // The length of the path, i.e., the number of vertex descriptors // that will be coming in the next msg_path_vertices message. std::size_t path_length; template<typename Archiver> void serialize(Archiver& ar, const unsigned int /*version*/) { ar & edge & path_length; } }; // Payload for a msg_tree_header message. template<typename Vertex, typename Edge> struct tree_header { // The edge over which the tree is being communicated Edge edge; // Gamma, which is the eta of the sender. Vertex gamma; // The length of the list of vertices in the bicomponent, i.e., // the number of vertex descriptors that will be coming in the // next msg_tree_vertices message. std::size_t bicomp_length; template<typename Archiver> void serialize(Archiver& ar, const unsigned int /*version*/) { ar & edge & gamma & bicomp_length; } }; // Payload for the msg_name message. template<typename EdgeDescriptor> struct name_header { // The edge being communicated and named. EdgeDescriptor edge; // The 0-based name of the component std::size_t name; template<typename Archiver> void serialize(Archiver& ar, const unsigned int /*version*/) { ar & edge & name; } }; /* Computes the branch point between two paths. The branch point is the last position at which both paths are equivalent, beyond which the paths diverge. Both paths must have length > 0 and the initial elements of the paths must be equal. This is guaranteed in Hohberg's algorithm because all paths start at the leader. Returns the value at the branch point. */ template<typename T> T branch_point(const std::vector<T>& p1, const std::vector<T>& p2) { BOOST_ASSERT(!p1.empty()); BOOST_ASSERT(!p2.empty()); BOOST_ASSERT(p1.front() == p2.front()); typedef typename std::vector<T>::const_iterator iterator; iterator mismatch_pos; if (p1.size() <= p2.size()) mismatch_pos = std::mismatch(p1.begin(), p1.end(), p2.begin()).first; else mismatch_pos = std::mismatch(p2.begin(), p2.end(), p1.begin()).first; --mismatch_pos; return *mismatch_pos; } /* Computes the infimum of vertices a and b in the given path. The infimum is the largest element that is on the paths from a to the root and from b to the root. */ template<typename T> T infimum(const std::vector<T>& parent_path, T a, T b) { using std::swap; typedef typename std::vector<T>::const_iterator iterator; iterator first = parent_path.begin(), last = parent_path.end(); #if defined(PBGL_HOHBERG_DEBUG) && PBGL_HOHBERG_DEBUG > 2 std::cerr << "infimum("; for (iterator i = first; i != last; ++i) { if (i != first) std::cerr << ' '; std::cerr << local(*i) << '@' << owner(*i); } std::cerr << ", " << local(a) << '@' << owner(a) << ", " << local(b) << '@' << owner(b) << ") = "; #endif if (a == b) { #if defined(PBGL_HOHBERG_DEBUG) && PBGL_HOHBERG_DEBUG > 2 std::cerr << local(a) << '@' << owner(a) << std::endl; #endif return a; } // Try to find a or b, whichever is closest to the end --last; while (*last != a) { // If we match b, swap the two so we'll be looking for b later. if (*last == b) { swap(a,b); break; } if (last == first) { #if defined(PBGL_HOHBERG_DEBUG) && PBGL_HOHBERG_DEBUG > 2 std::cerr << local(*first) << '@' << owner(*first) << std::endl; #endif return *first; } else --last; } // Try to find b (which may originally have been a) while (*last != b) { if (last == first) { #if defined(PBGL_HOHBERG_DEBUG) && PBGL_HOHBERG_DEBUG > 2 std::cerr << local(*first) << '@' << owner(*first) << std::endl; #endif return *first; } else --last; } #if defined(PBGL_HOHBERG_DEBUG) && PBGL_HOHBERG_DEBUG > 2 std::cerr << local(*last) << '@' << owner(*last) << std::endl; #endif // We've found b; it's the infimum. return *last; } } // end namespace hohberg_detail /* A message that will be stored for each edge by Hohberg's algorithm. */ template<typename Graph> struct hohberg_message { typedef typename graph_traits<Graph>::vertex_descriptor Vertex; typedef typename graph_traits<Graph>::edge_descriptor Edge; // Assign from a path message void assign(const std::vector<Vertex>& path) { gamma = graph_traits<Graph>::null_vertex(); path_or_bicomp = path; } // Assign from a tree message void assign(Vertex gamma, const std::vector<Vertex>& in_same_bicomponent) { this->gamma = gamma; path_or_bicomp = in_same_bicomponent; } bool is_path() const { return gamma == graph_traits<Graph>::null_vertex(); } bool is_tree() const { return gamma != graph_traits<Graph>::null_vertex(); } /// The "gamma" of a tree message, or null_vertex() for a path message Vertex gamma; // Either the path for a path message or the in_same_bicomponent std::vector<Vertex> path_or_bicomp; }; /* An abstraction of a vertex processor in Hohberg's algorithm. The hohberg_vertex_processor class is responsible for processing messages transmitted to it via its edges. */ template<typename Graph> class hohberg_vertex_processor { typedef typename graph_traits<Graph>::vertex_descriptor Vertex; typedef typename graph_traits<Graph>::edge_descriptor Edge; typedef typename graph_traits<Graph>::degree_size_type degree_size_type; typedef typename graph_traits<Graph>::edges_size_type edges_size_type; typedef typename boost::graph::parallel::process_group_type<Graph>::type ProcessGroup; typedef std::vector<Vertex> path_t; typedef typename path_t::iterator path_iterator; public: hohberg_vertex_processor() : phase(1), parent(graph_traits<Graph>::null_vertex()), eta(graph_traits<Graph>::null_vertex()) { } // Called to initialize a leader in the algorithm, which involves // sending out the initial path messages and being ready to receive // them. void initialize_leader(Vertex alpha, const Graph& g); /// Handle a path message on edge e. The path will be destroyed by /// this operation. void operator()(Edge e, path_t& path, const Graph& g); /// Handle a tree message on edge e. in_same_bicomponent will be /// destroyed by this operation. void operator()(Edge e, Vertex gamma, path_t& in_same_bicomponent, const Graph& g); // Handle a name message. void operator()(Edge e, edges_size_type name, const Graph& g); // Retrieve the phase unsigned char get_phase() const { return phase; } // Start the naming phase. The current phase must be 3 (echo), and // the offset contains the offset at which this processor should // begin when labelling its bicomponents. The offset is just a // parallel prefix sum of the number of bicomponents in each // processor that precedes it (globally). void start_naming_phase(Vertex alpha, const Graph& g, edges_size_type offset); /* Determine the number of bicomponents that we will be naming * ourselves. */ edges_size_type num_starting_bicomponents(Vertex alpha, const Graph& g); // Fill in the edge component map with biconnected component // numbers. template<typename ComponentMap> void fill_edge_map(Vertex alpha, const Graph& g, ComponentMap& component); protected: /* Start the echo phase (phase 3) where we propagate information up the tree. */ void echo_phase(Vertex alpha, const Graph& g); /* Retrieve the index of edge in the out-edges list of target(e, g). */ std::size_t get_edge_index(Edge e, const Graph& g); /* Retrieve the index of the edge incidence on v in the out-edges list of vertex u. */ std::size_t get_incident_edge_index(Vertex u, Vertex v, const Graph& g); /* Keeps track of which phase of the algorithm we are in. There are * four phases plus the "finished" phase: * * 1) Building the spanning tree * 2) Discovering cycles * 3) Echoing back up the spanning tree * 4) Labelling biconnected components * 5) Finished */ unsigned char phase; /* The parent of this vertex in the spanning tree. This value will be graph_traits<Graph>::null_vertex() for the leader. */ Vertex parent; /* The farthest ancestor up the tree that resides in the same biconnected component as we do. This information is approximate: we might not know about the actual farthest ancestor, but this is the farthest one we've seen so far. */ Vertex eta; /* The number of edges that have not yet transmitted any messages to us. This starts at the degree of the vertex and decreases as we receive messages. When this counter hits zero, we're done with the second phase of the algorithm. In Hohberg's paper, the actual remaining edge set E is stored with termination when all edges have been removed from E, but we only need to detect termination so the set E need not be explicitly represented. */ degree_size_type num_edges_not_transmitted; /* The path from the root of the spanning tree to this vertex. This vector will always store only the parts of the path leading up to this vertex, and not the vertex itself. Thus, it will be empty for the leader. */ std::vector<Vertex> path_from_root; /* Structure containing all of the extra data we need to keep around PER EDGE. This information can not be contained within a property map, because it can't be shared among vertices without breaking the algorithm. Decreasing the size of this structure will drastically */ struct per_edge_data { hohberg_message<Graph> msg; std::vector<Vertex> M; bool is_tree_edge; degree_size_type partition; }; /* Data for each edge in the graph. This structure will be indexed by the position of the edge in the out_edges() list. */ std::vector<per_edge_data> edge_data; /* The mapping from local partition numbers (0..n-1) to global partition numbers. */ std::vector<edges_size_type> local_to_global_partitions; friend class boost::serialization::access; // We cannot actually serialize a vertex processor, nor would we // want to. However, the fact that we're putting instances into a // distributed_property_map means that we need to have a serialize() // function available. template<typename Archiver> void serialize(Archiver&, const unsigned int /*version*/) { BOOST_ASSERT(false); } }; template<typename Graph> void hohberg_vertex_processor<Graph>::initialize_leader(Vertex alpha, const Graph& g) { using namespace hohberg_detail; ProcessGroup pg = process_group(g); typename property_map<Graph, vertex_owner_t>::const_type owner = get(vertex_owner, g); path_header<Edge> header; header.path_length = 1; BGL_FORALL_OUTEDGES_T(alpha, e, g, Graph) { header.edge = e; send(pg, get(owner, target(e, g)), msg_path_header, header); send(pg, get(owner, target(e, g)), msg_path_vertices, alpha); } num_edges_not_transmitted = degree(alpha, g); edge_data.resize(num_edges_not_transmitted); phase = 2; } template<typename Graph> void hohberg_vertex_processor<Graph>::operator()(Edge e, path_t& path, const Graph& g) { using namespace hohberg_detail; typename property_map<Graph, vertex_owner_t>::const_type owner = get(vertex_owner, g); #ifdef PBGL_HOHBERG_DEBUG // std::cerr << local(source(e, g)) << '@' << owner(source(e, g)) << " -> " // << local(target(e, g)) << '@' << owner(target(e, g)) << ": path("; // for (std::size_t i = 0; i < path.size(); ++i) { // if (i > 0) std::cerr << ' '; // std::cerr << local(path[i]) << '@' << owner(path[i]); // } std::cerr << "), phase = " << (int)phase << std::endl; #endif // Get access to edge-specific data if (edge_data.empty()) edge_data.resize(degree(target(e, g), g)); per_edge_data& edata = edge_data[get_edge_index(e, g)]; // Record the message. We'll need it in phase 3. edata.msg.assign(path); // Note: "alpha" refers to the vertex "processor" receiving the // message. Vertex alpha = target(e, g); switch (phase) { case 1: { num_edges_not_transmitted = degree(alpha, g) - 1; edata.is_tree_edge = true; parent = path.back(); eta = parent; edata.M.clear(); edata.M.push_back(parent); // Broadcast the path from the root to our potential children in // the spanning tree. path.push_back(alpha); path_header<Edge> header; header.path_length = path.size(); ProcessGroup pg = process_group(g); BGL_FORALL_OUTEDGES_T(alpha, oe, g, Graph) { // Skip the tree edge we just received if (target(oe, g) != source(e, g)) { header.edge = oe; send(pg, get(owner, target(oe, g)), msg_path_header, header); send(pg, get(owner, target(oe, g)), msg_path_vertices, &path[0], header.path_length); } } path.pop_back(); // Swap the old path in, to save some extra copying. Nobody path_from_root.swap(path); // Once we have received our place in the spanning tree, move on // to phase 2. phase = 2; // If we only had only edge, skip to phase 3. if (num_edges_not_transmitted == 0) echo_phase(alpha, g); return; } case 2: { --num_edges_not_transmitted; edata.is_tree_edge = false; // Determine if alpha (our vertex) is in the path path_iterator pos = std::find(path.begin(), path.end(), alpha); if (pos != path.end()) { // Case A: There is a cycle alpha beta ... gamma alpha // M(e) <- {beta, gammar} edata.M.clear(); ++pos; // If pos == path.end(), we have a self-loop if (pos != path.end()) { // Add beta edata.M.push_back(*pos); ++pos; } // If pos == path.end(), we have a self-loop or beta == gamma // (parallel edge). Otherwise, add gamma. if (pos != path.end()) edata.M.push_back(path.back()); } else { // Case B: There is a cycle but we haven't seen alpha yet. // M(e) = {parent, path.back()} edata.M.clear(); edata.M.push_back(path.back()); if (parent != path.back()) edata.M.push_back(parent); // eta = inf(eta, bra(pi_t, pi)) eta = infimum(path_from_root, eta, branch_point(path_from_root, path)); } if (num_edges_not_transmitted == 0) echo_phase(alpha, g); break; } default: // std::cerr << "Phase is " << int(phase) << "\n"; BOOST_ASSERT(false); } } template<typename Graph> void hohberg_vertex_processor<Graph>::operator()(Edge e, Vertex gamma, path_t& in_same_bicomponent, const Graph& g) { using namespace hohberg_detail; #ifdef PBGL_HOHBERG_DEBUG std::cerr << local(source(e, g)) << '@' << owner(source(e, g)) << " -> " << local(target(e, g)) << '@' << owner(target(e, g)) << ": tree(" << local(gamma) << '@' << owner(gamma) << ", "; for (std::size_t i = 0; i < in_same_bicomponent.size(); ++i) { if (i > 0) std::cerr << ' '; std::cerr << local(in_same_bicomponent[i]) << '@' << owner(in_same_bicomponent[i]); } std::cerr << ", " << local(source(e, g)) << '@' << owner(source(e, g)) << "), phase = " << (int)phase << std::endl; #endif // Get access to edge-specific data per_edge_data& edata = edge_data[get_edge_index(e, g)]; // Record the message. We'll need it in phase 3. edata.msg.assign(gamma, in_same_bicomponent); // Note: "alpha" refers to the vertex "processor" receiving the // message. Vertex alpha = target(e, g); Vertex beta = source(e, g); switch (phase) { case 2: --num_edges_not_transmitted; edata.is_tree_edge = true; if (gamma == alpha) { // Case C edata.M.swap(in_same_bicomponent); } else { // Case D edata.M.clear(); edata.M.push_back(parent); if (beta != parent) edata.M.push_back(beta); eta = infimum(path_from_root, eta, gamma); } if (num_edges_not_transmitted == 0) echo_phase(alpha, g); break; default: BOOST_ASSERT(false); } } template<typename Graph> void hohberg_vertex_processor<Graph>::operator()(Edge e, edges_size_type name, const Graph& g) { using namespace hohberg_detail; #ifdef PBGL_HOHBERG_DEBUG std::cerr << local(source(e, g)) << '@' << owner(source(e, g)) << " -> " << local(target(e, g)) << '@' << owner(target(e, g)) << ": name(" << name << "), phase = " << (int)phase << std::endl; #endif BOOST_ASSERT(phase == 4); typename property_map<Graph, vertex_owner_t>::const_type owner = get(vertex_owner, g); // Send name messages along the spanning tree edges that are in the // same bicomponent as the edge to our parent. ProcessGroup pg = process_group(g); Vertex alpha = target(e, g); std::size_t idx = 0; BGL_FORALL_OUTEDGES_T(alpha, e, g, Graph) { per_edge_data& edata = edge_data[idx++]; if (edata.is_tree_edge && find(edata.M.begin(), edata.M.end(), parent) != edata.M.end() && target(e, g) != parent) { // Notify our children in the spanning tree of this name name_header<Edge> header; header.edge = e; header.name = name; send(pg, get(owner, target(e, g)), msg_name, header); } else if (target(e, g) == parent) { // Map from local partition numbers to global bicomponent numbers local_to_global_partitions[edata.partition] = name; } } // Final stage phase = 5; } template<typename Graph> typename hohberg_vertex_processor<Graph>::edges_size_type hohberg_vertex_processor<Graph>:: num_starting_bicomponents(Vertex alpha, const Graph& g) { edges_size_type not_mapped = (std::numeric_limits<edges_size_type>::max)(); edges_size_type result = 0; std::size_t idx = 0; BGL_FORALL_OUTEDGES_T(alpha, e, g, Graph) { per_edge_data& edata = edge_data[idx++]; if (edata.is_tree_edge && find(edata.M.begin(), edata.M.end(), parent) == edata.M.end()) { // Map from local partition numbers to global bicomponent numbers if (local_to_global_partitions[edata.partition] == not_mapped) local_to_global_partitions[edata.partition] = result++; } } #ifdef PBGL_HOHBERG_DEBUG std::cerr << local(alpha) << '@' << owner(alpha) << " has " << result << " bicomponents originating at it." << std::endl; #endif return result; } template<typename Graph> template<typename ComponentMap> void hohberg_vertex_processor<Graph>:: fill_edge_map(Vertex alpha, const Graph& g, ComponentMap& component) { std::size_t idx = 0; BGL_FORALL_OUTEDGES_T(alpha, e, g, Graph) { per_edge_data& edata = edge_data[idx++]; local_put(component, e, local_to_global_partitions[edata.partition]); #if defined(PBGL_HOHBERG_DEBUG) && PBGL_HOHBERG_DEBUG > 2 std::cerr << "component(" << local(source(e, g)) << '@' << owner(source(e, g)) << " -> " << local(target(e, g)) << '@' << owner(target(e, g)) << ") = " << local_to_global_partitions[edata.partition] << " (partition = " << edata.partition << " of " << local_to_global_partitions.size() << ")" << std::endl; #endif } } template<typename Graph> void hohberg_vertex_processor<Graph>:: start_naming_phase(Vertex alpha, const Graph& g, edges_size_type offset) { using namespace hohberg_detail; BOOST_ASSERT(phase == 4); typename property_map<Graph, vertex_owner_t>::const_type owner = get(vertex_owner, g); // Send name messages along the spanning tree edges of the // components that we get to number. ProcessGroup pg = process_group(g); bool has_more_children_to_name = false; // Map from local partition numbers to global bicomponent numbers edges_size_type not_mapped = (std::numeric_limits<edges_size_type>::max)(); for (std::size_t i = 0; i < local_to_global_partitions.size(); ++i) { if (local_to_global_partitions[i] != not_mapped) local_to_global_partitions[i] += offset; } std::size_t idx = 0; BGL_FORALL_OUTEDGES_T(alpha, e, g, Graph) { per_edge_data& edata = edge_data[idx++]; if (edata.is_tree_edge && find(edata.M.begin(), edata.M.end(), parent) == edata.M.end()) { // Notify our children in the spanning tree of this new name name_header<Edge> header; header.edge = e; header.name = local_to_global_partitions[edata.partition]; send(pg, get(owner, target(e, g)), msg_name, header); } else if (edata.is_tree_edge) { has_more_children_to_name = true; } #if defined(PBGL_HOHBERG_DEBUG) && PBGL_HOHBERG_DEBUG > 2 std::cerr << "M[" << local(source(e, g)) << '@' << owner(source(e, g)) << " -> " << local(target(e, g)) << '@' << owner(target(e, g)) << "] = "; for (std::size_t i = 0; i < edata.M.size(); ++i) { std::cerr << local(edata.M[i]) << '@' << owner(edata.M[i]) << ' '; } std::cerr << std::endl; #endif } // See if we're done. if (!has_more_children_to_name) // Final stage phase = 5; } template<typename Graph> void hohberg_vertex_processor<Graph>::echo_phase(Vertex alpha, const Graph& g) { using namespace hohberg_detail; typename property_map<Graph, vertex_owner_t>::const_type owner = get(vertex_owner, g); /* We're entering the echo phase. */ phase = 3; if (parent != graph_traits<Graph>::null_vertex()) { Edge edge_to_parent; #if defined(PBGL_HOHBERG_DEBUG) && PBGL_HOHBERG_DEBUG > 1 std::cerr << local(alpha) << '@' << owner(alpha) << " echo: parent = " << local(parent) << '@' << owner(parent) << ", eta = " << local(eta) << '@' << owner(eta) << ", Gamma = "; #endif std::vector<Vertex> bicomp; std::size_t e_index = 0; BGL_FORALL_OUTEDGES_T(alpha, e, g, Graph) { if (target(e, g) == parent && parent == eta) { edge_to_parent = e; if (find(bicomp.begin(), bicomp.end(), alpha) == bicomp.end()) { #if defined(PBGL_HOHBERG_DEBUG) && PBGL_HOHBERG_DEBUG > 1 std::cerr << local(alpha) << '@' << owner(alpha) << ' '; #endif bicomp.push_back(alpha); } } else { if (target(e, g) == parent) edge_to_parent = e; per_edge_data& edata = edge_data[e_index]; if (edata.msg.is_path()) { path_iterator pos = std::find(edata.msg.path_or_bicomp.begin(), edata.msg.path_or_bicomp.end(), eta); if (pos != edata.msg.path_or_bicomp.end()) { ++pos; if (pos != edata.msg.path_or_bicomp.end() && find(bicomp.begin(), bicomp.end(), *pos) == bicomp.end()) { #if defined(PBGL_HOHBERG_DEBUG) && PBGL_HOHBERG_DEBUG > 1 std::cerr << local(*pos) << '@' << owner(*pos) << ' '; #endif bicomp.push_back(*pos); } } } else if (edata.msg.is_tree() && edata.msg.gamma == eta) { for (path_iterator i = edata.msg.path_or_bicomp.begin(); i != edata.msg.path_or_bicomp.end(); ++i) { if (find(bicomp.begin(), bicomp.end(), *i) == bicomp.end()) { #if defined(PBGL_HOHBERG_DEBUG) && PBGL_HOHBERG_DEBUG > 1 std::cerr << local(*i) << '@' << owner(*i) << ' '; #endif bicomp.push_back(*i); } } } } ++e_index; } #ifdef PBGL_HOHBERG_DEBUG std::cerr << std::endl; #endif // Send tree(eta, bicomp) to parent tree_header<Vertex, Edge> header; header.edge = edge_to_parent; header.gamma = eta; header.bicomp_length = bicomp.size(); ProcessGroup pg = process_group(g); send(pg, get(owner, parent), msg_tree_header, header); send(pg, get(owner, parent), msg_tree_vertices, &bicomp[0], header.bicomp_length); } // Compute the partition of edges such that iff two edges e1 and e2 // are in different subsets then M(e1) is disjoint from M(e2). // Start by putting each edge in a different partition std::vector<degree_size_type> parent_vec(edge_data.size()); degree_size_type idx = 0; for (idx = 0; idx < edge_data.size(); ++idx) parent_vec[idx] = idx; // Go through each edge e, performing a union() on the edges // incident on all vertices in M[e]. idx = 0; BGL_FORALL_OUTEDGES_T(alpha, e, g, Graph) { per_edge_data& edata = edge_data[idx++]; // Compute union of vertices in M if (!edata.M.empty()) { degree_size_type e1 = get_incident_edge_index(alpha, edata.M.front(), g); while (parent_vec[e1] != e1) e1 = parent_vec[e1]; for (std::size_t i = 1; i < edata.M.size(); ++i) { degree_size_type e2 = get_incident_edge_index(alpha, edata.M[i], g); while (parent_vec[e2] != e2) e2 = parent_vec[e2]; parent_vec[e2] = e1; } } } edges_size_type not_mapped = (std::numeric_limits<edges_size_type>::max)(); // Determine the number of partitions for (idx = 0; idx < parent_vec.size(); ++idx) { if (parent_vec[idx] == idx) { edge_data[idx].partition = local_to_global_partitions.size(); local_to_global_partitions.push_back(not_mapped); } } // Assign partition numbers to each edge for (idx = 0; idx < parent_vec.size(); ++idx) { degree_size_type rep = parent_vec[idx]; while (rep != parent_vec[rep]) rep = parent_vec[rep]; edge_data[idx].partition = edge_data[rep].partition; } // Enter the naming phase (but don't send anything yet). phase = 4; } template<typename Graph> std::size_t hohberg_vertex_processor<Graph>::get_edge_index(Edge e, const Graph& g) { std::size_t result = 0; BGL_FORALL_OUTEDGES_T(target(e, g), oe, g, Graph) { if (source(e, g) == target(oe, g)) return result; ++result; } BOOST_ASSERT(false); } template<typename Graph> std::size_t hohberg_vertex_processor<Graph>::get_incident_edge_index(Vertex u, Vertex v, const Graph& g) { std::size_t result = 0; BGL_FORALL_OUTEDGES_T(u, e, g, Graph) { if (target(e, g) == v) return result; ++result; } BOOST_ASSERT(false); } template<typename Graph, typename InputIterator, typename ComponentMap, typename VertexProcessorMap> typename graph_traits<Graph>::edges_size_type hohberg_biconnected_components (const Graph& g, ComponentMap component, InputIterator first, InputIterator last, VertexProcessorMap vertex_processor) { using namespace boost::graph::parallel; using namespace hohberg_detail; using boost::parallel::all_reduce; typename property_map<Graph, vertex_owner_t>::const_type owner = get(vertex_owner, g); // The graph must be undirected BOOST_STATIC_ASSERT( (is_convertible<typename graph_traits<Graph>::directed_category, undirected_tag>::value)); // The graph must model Incidence Graph BOOST_CONCEPT_ASSERT(( IncidenceGraphConcept<Graph> )); typedef typename graph_traits<Graph>::edges_size_type edges_size_type; typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor; typedef typename graph_traits<Graph>::edge_descriptor edge_descriptor; // Retrieve the process group we will use for communication typedef typename process_group_type<Graph>::type process_group_type; process_group_type pg = process_group(g); // Keeps track of the edges that we know to be tree edges. std::vector<edge_descriptor> tree_edges; // The leaders send out a path message to initiate the algorithm while (first != last) { vertex_descriptor leader = *first; if (process_id(pg) == get(owner, leader)) vertex_processor[leader].initialize_leader(leader, g); ++first; } synchronize(pg); // Will hold the number of bicomponents in the graph. edges_size_type num_bicomponents = 0; // Keep track of the path length that we should expect, based on the // level in the breadth-first search tree. At present, this is only // used as a sanity check. TBD: This could be used to decrease the // amount of communication required per-edge (by about 4 bytes). std::size_t path_length = 1; typedef std::vector<vertex_descriptor> path_t; unsigned char minimum_phase = 5; do { while (optional<std::pair<int, int> > msg = probe(pg)) { switch (msg->second) { case msg_path_header: { // Receive the path header path_header<edge_descriptor> header; receive(pg, msg->first, msg->second, header); BOOST_ASSERT(path_length == header.path_length); // Receive the path itself path_t path(path_length); receive(pg, msg->first, msg_path_vertices, &path[0], path_length); edge_descriptor e = header.edge; vertex_processor[target(e, g)](e, path, g); } break; case msg_path_vertices: // Should be handled in msg_path_header case, unless we're going // stateless. BOOST_ASSERT(false); break; case msg_tree_header: { // Receive the tree header tree_header<vertex_descriptor, edge_descriptor> header; receive(pg, msg->first, msg->second, header); // Receive the tree itself path_t in_same_bicomponent(header.bicomp_length); receive(pg, msg->first, msg_tree_vertices, &in_same_bicomponent[0], header.bicomp_length); edge_descriptor e = header.edge; vertex_processor[target(e, g)](e, header.gamma, in_same_bicomponent, g); } break; case msg_tree_vertices: // Should be handled in msg_tree_header case, unless we're // going stateless. BOOST_ASSERT(false); break; case msg_name: { name_header<edge_descriptor> header; receive(pg, msg->first, msg->second, header); edge_descriptor e = header.edge; vertex_processor[target(e, g)](e, header.name, g); } break; default: BOOST_ASSERT(false); } } ++path_length; // Compute minimum phase locally minimum_phase = 5; unsigned char maximum_phase = 1; BGL_FORALL_VERTICES_T(v, g, Graph) { minimum_phase = (std::min)(minimum_phase, vertex_processor[v].get_phase()); maximum_phase = (std::max)(maximum_phase, vertex_processor[v].get_phase()); } #ifdef PBGL_HOHBERG_DEBUG if (process_id(pg) == 0) std::cerr << "<---------End of stage------------->" << std::endl; #endif // Compute minimum phase globally minimum_phase = all_reduce(pg, minimum_phase, boost::mpi::minimum<char>()); #ifdef PBGL_HOHBERG_DEBUG if (process_id(pg) == 0) std::cerr << "Minimum phase = " << (int)minimum_phase << std::endl; #endif if (minimum_phase == 4 && all_reduce(pg, maximum_phase, boost::mpi::maximum<char>()) == 4) { #ifdef PBGL_HOHBERG_DEBUG if (process_id(pg) == 0) std::cerr << "<---------Naming phase------------->" << std::endl; #endif // Compute the biconnected component number offsets for each // vertex. std::vector<edges_size_type> local_offsets; local_offsets.reserve(num_vertices(g)); edges_size_type num_local_bicomponents = 0; BGL_FORALL_VERTICES_T(v, g, Graph) { local_offsets.push_back(num_local_bicomponents); num_local_bicomponents += vertex_processor[v].num_starting_bicomponents(v, g); } synchronize(pg); // Find our the number of bicomponent names that will originate // from each process. This tells us how many bicomponents are in // the entire graph and what our global offset is for computing // our own biconnected component names. std::vector<edges_size_type> all_bicomponents(num_processes(pg)); all_gather(pg, &num_local_bicomponents, &num_local_bicomponents + 1, all_bicomponents); num_bicomponents = 0; edges_size_type my_global_offset = 0; for (std::size_t i = 0; i < all_bicomponents.size(); ++i) { if (i == (std::size_t)process_id(pg)) my_global_offset = num_bicomponents; num_bicomponents += all_bicomponents[i]; } std::size_t index = 0; BGL_FORALL_VERTICES_T(v, g, Graph) { edges_size_type offset = my_global_offset + local_offsets[index++]; vertex_processor[v].start_naming_phase(v, g, offset); } } synchronize(pg); } while (minimum_phase < 5); // Number the edges appropriately. BGL_FORALL_VERTICES_T(v, g, Graph) vertex_processor[v].fill_edge_map(v, g, component); return num_bicomponents; } template<typename Graph, typename ComponentMap, typename InputIterator> typename graph_traits<Graph>::edges_size_type hohberg_biconnected_components (const Graph& g, ComponentMap component, InputIterator first, InputIterator last) { std::vector<hohberg_vertex_processor<Graph> > vertex_processors(num_vertices(g)); return hohberg_biconnected_components (g, component, first, last, make_iterator_property_map(vertex_processors.begin(), get(vertex_index, g))); } template<typename Graph, typename ComponentMap, typename ParentMap> typename graph_traits<Graph>::edges_size_type hohberg_biconnected_components(const Graph& g, ComponentMap component, ParentMap parent) { // We need the connected components of the graph, but we don't care // about component numbers. connected_components(g, dummy_property_map(), parent); // Each root in the parent map is a leader typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor; std::vector<vertex_descriptor> leaders; BGL_FORALL_VERTICES_T(v, g, Graph) if (get(parent, v) == v) leaders.push_back(v); return hohberg_biconnected_components(g, component, leaders.begin(), leaders.end()); } template<typename Graph, typename ComponentMap> typename graph_traits<Graph>::edges_size_type hohberg_biconnected_components(const Graph& g, ComponentMap component) { typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor; std::vector<vertex_descriptor> parents(num_vertices(g)); return hohberg_biconnected_components (g, component, make_iterator_property_map(parents.begin(), get(vertex_index, g))); } } } } // end namespace boost::graph::distributed #endif // BOOST_GRAPH_DISTRIBUTED_HOHBERG_BICONNECTED_COMPONENTS_HPP