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