PageRenderTime 202ms CodeModel.GetById 7ms app.highlight 175ms RepoModel.GetById 1ms app.codeStats 1ms

/Src/Dependencies/Boost/boost/graph/distributed/hohberg_biconnected_components.hpp

http://hadesmem.googlecode.com/
C++ Header | 1129 lines | 743 code | 169 blank | 217 comment | 149 complexity | 106ced9a971d843c2228588366a39eae MD5 | raw file
   1// Copyright 2005 The Trustees of Indiana University.
   2
   3// Use, modification and distribution is subject to the Boost Software
   4// License, Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
   5// http://www.boost.org/LICENSE_1_0.txt)
   6
   7//  Authors: Douglas Gregor
   8//           Andrew Lumsdaine
   9
  10// An implementation of Walter Hohberg's distributed biconnected
  11// components algorithm, from:
  12//
  13//   Walter Hohberg. How to Find Biconnected Components in Distributed
  14//   Networks. J. Parallel Distrib. Comput., 9(4):374-386, 1990.
  15//
  16#ifndef BOOST_GRAPH_DISTRIBUTED_HOHBERG_BICONNECTED_COMPONENTS_HPP
  17#define BOOST_GRAPH_DISTRIBUTED_HOHBERG_BICONNECTED_COMPONENTS_HPP
  18
  19#ifndef BOOST_GRAPH_USE_MPI
  20#error "Parallel BGL files should not be included unless <boost/graph/use_mpi.hpp> has been included"
  21#endif
  22
  23/* You can define PBGL_HOHBERG_DEBUG to an integer value (1, 2, or 3)
  24 * to enable debugging information. 1 includes only the phases of the
  25 * algorithm and messages as their are received. 2 and 3 add
  26 * additional levels of detail about internal data structures related
  27 * to the algorithm itself.
  28 *
  29 * #define PBGL_HOHBERG_DEBUG 1
  30*/
  31
  32#include <boost/graph/graph_traits.hpp>
  33#include <boost/graph/parallel/container_traits.hpp>
  34#include <boost/graph/parallel/process_group.hpp>
  35#include <boost/static_assert.hpp>
  36#include <boost/mpi/operations.hpp>
  37#include <boost/type_traits/is_convertible.hpp>
  38#include <boost/graph/graph_concepts.hpp>
  39#include <boost/graph/iteration_macros.hpp>
  40#include <boost/optional.hpp>
  41#include <utility> // for std::pair
  42#include <boost/assert.hpp>
  43#include <algorithm> // for std::find, std::mismatch
  44#include <vector>
  45#include <boost/graph/parallel/algorithm.hpp>
  46#include <boost/graph/distributed/connected_components.hpp>
  47
  48namespace boost { namespace graph { namespace distributed {
  49
  50namespace hohberg_detail {
  51  enum message_kind {
  52    /* A header for the PATH message, stating which edge the message
  53       is coming on and how many vertices will be following. The data
  54       structure communicated will be a path_header. */
  55    msg_path_header,
  56    /* A message containing the vertices that make up a path. It will
  57       always follow a msg_path_header and will contain vertex
  58       descriptors, only. */
  59    msg_path_vertices,
  60    /* A header for the TREE message, stating the value of gamma and
  61       the number of vertices to come in the following
  62       msg_tree_vertices. */
  63    msg_tree_header,
  64    /* A message containing the vertices that make up the set of
  65       vertices in the same bicomponent as the sender. It will always
  66       follow a msg_tree_header and will contain vertex descriptors,
  67       only. */
  68    msg_tree_vertices,
  69    /* Provides a name for the biconnected component of the edge. */
  70    msg_name
  71  };
  72
  73  // Payload for a msg_path_header message.
  74  template<typename EdgeDescriptor>
  75  struct path_header
  76  {
  77    // The edge over which the path is being communicated
  78    EdgeDescriptor edge;
  79
  80    // The length of the path, i.e., the number of vertex descriptors
  81    // that will be coming in the next msg_path_vertices message.
  82    std::size_t    path_length;
  83
  84    template<typename Archiver>
  85    void serialize(Archiver& ar, const unsigned int /*version*/)
  86    {
  87      ar & edge & path_length;
  88    }
  89  };
  90
  91  // Payload for a msg_tree_header message.
  92  template<typename Vertex, typename Edge>
  93  struct tree_header
  94  {
  95    // The edge over which the tree is being communicated
  96    Edge edge;
  97
  98    // Gamma, which is the eta of the sender.
  99    Vertex gamma;
 100
 101    // The length of the list of vertices in the bicomponent, i.e.,
 102    // the number of vertex descriptors that will be coming in the
 103    // next msg_tree_vertices message.
 104    std::size_t    bicomp_length;
 105
 106    template<typename Archiver>
 107    void serialize(Archiver& ar, const unsigned int /*version*/)
 108    {
 109      ar & edge & gamma & bicomp_length;
 110    }
 111  };
 112
 113  // Payload for the msg_name message.
 114  template<typename EdgeDescriptor>
 115  struct name_header
 116  {
 117    // The edge being communicated and named.
 118    EdgeDescriptor edge;
 119
 120    // The 0-based name of the component
 121    std::size_t name;
 122
 123    template<typename Archiver>
 124    void serialize(Archiver& ar, const unsigned int /*version*/)
 125    {
 126      ar & edge & name;
 127    }
 128  };
 129
 130  /* Computes the branch point between two paths. The branch point is
 131     the last position at which both paths are equivalent, beyond
 132     which the paths diverge. Both paths must have length > 0 and the
 133     initial elements of the paths must be equal. This is guaranteed
 134     in Hohberg's algorithm because all paths start at the
 135     leader. Returns the value at the branch point. */
 136  template<typename T>
 137  T branch_point(const std::vector<T>& p1, const std::vector<T>& p2)
 138  {
 139    BOOST_ASSERT(!p1.empty());
 140    BOOST_ASSERT(!p2.empty());
 141    BOOST_ASSERT(p1.front() == p2.front());
 142
 143    typedef typename std::vector<T>::const_iterator iterator;
 144
 145    iterator mismatch_pos;
 146    if (p1.size() <= p2.size())
 147      mismatch_pos = std::mismatch(p1.begin(), p1.end(), p2.begin()).first;
 148    else
 149      mismatch_pos = std::mismatch(p2.begin(), p2.end(), p1.begin()).first;
 150    --mismatch_pos;
 151    return *mismatch_pos;
 152  }
 153
 154  /* Computes the infimum of vertices a and b in the given path. The
 155     infimum is the largest element that is on the paths from a to the
 156     root and from b to the root. */
 157  template<typename T>
 158  T infimum(const std::vector<T>& parent_path, T a, T b)
 159  {
 160    using std::swap;
 161
 162    typedef typename std::vector<T>::const_iterator iterator;
 163    iterator first = parent_path.begin(), last = parent_path.end();
 164
 165#if defined(PBGL_HOHBERG_DEBUG) and PBGL_HOHBERG_DEBUG > 2
 166    std::cerr << "infimum(";
 167    for (iterator i = first; i != last; ++i) {
 168      if (i != first) std::cerr << ' ';
 169      std::cerr << local(*i) << '@' << owner(*i);
 170    }
 171    std::cerr << ", " << local(a) << '@' << owner(a) << ", "
 172              << local(b) << '@' << owner(b) << ") = ";
 173#endif
 174
 175    if (a == b) {
 176#if defined(PBGL_HOHBERG_DEBUG) && PBGL_HOHBERG_DEBUG > 2
 177      std::cerr << local(a) << '@' << owner(a) << std::endl;
 178#endif
 179      return a;
 180    }
 181
 182    // Try to find a or b, whichever is closest to the end
 183    --last;
 184    while (*last != a) {
 185      // If we match b, swap the two so we'll be looking for b later.
 186      if (*last == b) { swap(a,b); break; }
 187
 188      if (last == first) {
 189#if defined(PBGL_HOHBERG_DEBUG) && PBGL_HOHBERG_DEBUG > 2
 190        std::cerr << local(*first) << '@' << owner(*first) << std::endl;
 191#endif
 192        return *first;
 193      }
 194      else --last;
 195    }
 196
 197    // Try to find b (which may originally have been a)
 198    while (*last != b) {
 199      if (last == first) {
 200#if defined(PBGL_HOHBERG_DEBUG) and PBGL_HOHBERG_DEBUG > 2
 201        std::cerr << local(*first) << '@' << owner(*first) << std::endl;
 202#endif
 203        return *first;
 204      }
 205      else --last;
 206    }
 207
 208#if defined(PBGL_HOHBERG_DEBUG) && PBGL_HOHBERG_DEBUG > 2
 209    std::cerr << local(*last) << '@' << owner(*last) << std::endl;
 210#endif
 211    // We've found b; it's the infimum.
 212    return *last;
 213  }
 214} // end namespace hohberg_detail
 215
 216/* A message that will be stored for each edge by Hohberg's algorithm. */
 217template<typename Graph>
 218struct hohberg_message
 219{
 220  typedef typename graph_traits<Graph>::vertex_descriptor Vertex;
 221  typedef typename graph_traits<Graph>::edge_descriptor   Edge;
 222
 223  // Assign from a path message
 224  void assign(const std::vector<Vertex>& path)
 225  {
 226    gamma = graph_traits<Graph>::null_vertex();
 227    path_or_bicomp = path;
 228  }
 229
 230  // Assign from a tree message
 231  void assign(Vertex gamma, const std::vector<Vertex>& in_same_bicomponent)
 232  {
 233    this->gamma = gamma;
 234    path_or_bicomp = in_same_bicomponent;
 235  }
 236
 237  bool is_path() const { return gamma == graph_traits<Graph>::null_vertex(); }
 238  bool is_tree() const { return gamma != graph_traits<Graph>::null_vertex(); }
 239
 240  /// The "gamma" of a tree message, or null_vertex() for a path message
 241  Vertex gamma;
 242
 243  // Either the path for a path message or the in_same_bicomponent
 244  std::vector<Vertex> path_or_bicomp;
 245};
 246
 247
 248/* An abstraction of a vertex processor in Hohberg's algorithm. The
 249   hohberg_vertex_processor class is responsible for processing
 250   messages transmitted to it via its edges. */
 251template<typename Graph>
 252class hohberg_vertex_processor
 253{
 254  typedef typename graph_traits<Graph>::vertex_descriptor Vertex;
 255  typedef typename graph_traits<Graph>::edge_descriptor   Edge;
 256  typedef typename graph_traits<Graph>::degree_size_type  degree_size_type;
 257  typedef typename graph_traits<Graph>::edges_size_type   edges_size_type;
 258  typedef typename boost::graph::parallel::process_group_type<Graph>::type
 259    ProcessGroup;
 260  typedef std::vector<Vertex> path_t;
 261  typedef typename path_t::iterator path_iterator;
 262
 263 public:
 264  hohberg_vertex_processor()
 265    : phase(1),
 266      parent(graph_traits<Graph>::null_vertex()),
 267      eta(graph_traits<Graph>::null_vertex())
 268  {
 269  }
 270
 271  // Called to initialize a leader in the algorithm, which involves
 272  // sending out the initial path messages and being ready to receive
 273  // them.
 274  void initialize_leader(Vertex alpha, const Graph& g);
 275
 276  /// Handle a path message on edge e. The path will be destroyed by
 277  /// this operation.
 278  void
 279  operator()(Edge e, path_t& path, const Graph& g);
 280
 281  /// Handle a tree message on edge e. in_same_bicomponent will be
 282  /// destroyed by this operation.
 283  void
 284  operator()(Edge e, Vertex gamma, path_t& in_same_bicomponent,
 285             const Graph& g);
 286
 287  // Handle a name message.
 288  void operator()(Edge e, edges_size_type name, const Graph& g);
 289
 290  // Retrieve the phase
 291  unsigned char get_phase() const { return phase; }
 292
 293  // Start the naming phase. The current phase must be 3 (echo), and
 294  // the offset contains the offset at which this processor should
 295  // begin when labelling its bicomponents. The offset is just a
 296  // parallel prefix sum of the number of bicomponents in each
 297  // processor that precedes it (globally).
 298  void
 299  start_naming_phase(Vertex alpha, const Graph& g, edges_size_type offset);
 300
 301  /* Determine the number of bicomponents that we will be naming
 302   * ourselves.
 303   */
 304  edges_size_type num_starting_bicomponents(Vertex alpha, const Graph& g);
 305
 306  // Fill in the edge component map with biconnected component
 307  // numbers.
 308  template<typename ComponentMap>
 309  void fill_edge_map(Vertex alpha, const Graph& g, ComponentMap& component);
 310
 311 protected:
 312  /* Start the echo phase (phase 3) where we propagate information up
 313     the tree. */
 314  void echo_phase(Vertex alpha, const Graph& g);
 315
 316  /* Retrieve the index of edge in the out-edges list of target(e, g). */
 317  std::size_t get_edge_index(Edge e, const Graph& g);
 318
 319  /* Retrieve the index of the edge incidence on v in the out-edges
 320     list of vertex u. */
 321  std::size_t get_incident_edge_index(Vertex u, Vertex v, const Graph& g);
 322
 323  /* Keeps track of which phase of the algorithm we are in. There are
 324   * four phases plus the "finished" phase:
 325   *
 326   *   1) Building the spanning tree
 327   *   2) Discovering cycles
 328   *   3) Echoing back up the spanning tree
 329   *   4) Labelling biconnected components
 330   *   5) Finished
 331   */
 332  unsigned char phase;
 333
 334  /* The parent of this vertex in the spanning tree. This value will
 335     be graph_traits<Graph>::null_vertex() for the leader. */
 336  Vertex parent;
 337
 338  /* The farthest ancestor up the tree that resides in the same
 339     biconnected component as we do. This information is approximate:
 340     we might not know about the actual farthest ancestor, but this is
 341     the farthest one we've seen so far. */
 342  Vertex eta;
 343
 344  /* The number of edges that have not yet transmitted any messages to
 345     us. This starts at the degree of the vertex and decreases as we
 346     receive messages. When this counter hits zero, we're done with
 347     the second phase of the algorithm. In Hohberg's paper, the actual
 348     remaining edge set E is stored with termination when all edges
 349     have been removed from E, but we only need to detect termination
 350     so the set E need not be explicitly represented. */
 351  degree_size_type num_edges_not_transmitted;
 352
 353  /* The path from the root of the spanning tree to this vertex. This
 354     vector will always store only the parts of the path leading up to
 355     this vertex, and not the vertex itself. Thus, it will be empty
 356     for the leader. */
 357  std::vector<Vertex> path_from_root;
 358
 359  /* Structure containing all of the extra data we need to keep around
 360     PER EDGE. This information can not be contained within a property
 361     map, because it can't be shared among vertices without breaking
 362     the algorithm. Decreasing the size of this structure will drastically */
 363  struct per_edge_data
 364  {
 365    hohberg_message<Graph> msg;
 366    std::vector<Vertex> M;
 367    bool is_tree_edge;
 368    degree_size_type partition;
 369  };
 370
 371  /* Data for each edge in the graph. This structure will be indexed
 372     by the position of the edge in the out_edges() list. */
 373  std::vector<per_edge_data> edge_data;
 374
 375  /* The mapping from local partition numbers (0..n-1) to global
 376     partition numbers. */
 377  std::vector<edges_size_type> local_to_global_partitions;
 378
 379  friend class boost::serialization::access;
 380
 381  // We cannot actually serialize a vertex processor, nor would we
 382  // want to. However, the fact that we're putting instances into a
 383  // distributed_property_map means that we need to have a serialize()
 384  // function available.
 385  template<typename Archiver>
 386  void serialize(Archiver&, const unsigned int /*version*/)
 387  {
 388    BOOST_ASSERT(false);
 389  }
 390};
 391
 392template<typename Graph>
 393void
 394hohberg_vertex_processor<Graph>::initialize_leader(Vertex alpha,
 395                                                   const Graph& g)
 396{
 397  using namespace hohberg_detail;
 398
 399  ProcessGroup pg = process_group(g);
 400
 401  typename property_map<Graph, vertex_owner_t>::const_type
 402    owner = get(vertex_owner, g);
 403
 404  path_header<Edge> header;
 405  header.path_length = 1;
 406  BGL_FORALL_OUTEDGES_T(alpha, e, g, Graph) {
 407    header.edge = e;
 408    send(pg, get(owner, target(e, g)), msg_path_header, header);
 409    send(pg, get(owner, target(e, g)), msg_path_vertices, alpha);
 410  }
 411
 412  num_edges_not_transmitted = degree(alpha, g);
 413  edge_data.resize(num_edges_not_transmitted);
 414  phase = 2;
 415}
 416
 417template<typename Graph>
 418void
 419hohberg_vertex_processor<Graph>::operator()(Edge e, path_t& path,
 420                                            const Graph& g)
 421{
 422  using namespace hohberg_detail;
 423
 424  typename property_map<Graph, vertex_owner_t>::const_type
 425    owner = get(vertex_owner, g);
 426
 427#ifdef PBGL_HOHBERG_DEBUG
 428//  std::cerr << local(source(e, g)) << '@' << owner(source(e, g)) << " -> "
 429//            << local(target(e, g)) << '@' << owner(target(e, g)) << ": path(";
 430//  for (std::size_t i = 0; i < path.size(); ++i) {
 431//    if (i > 0) std::cerr << ' ';
 432//    std::cerr << local(path[i]) << '@' << owner(path[i]);
 433//  }
 434  std::cerr << "), phase = " << (int)phase << std::endl;
 435#endif
 436
 437  // Get access to edge-specific data
 438  if (edge_data.empty())
 439    edge_data.resize(degree(target(e, g), g));
 440  per_edge_data& edata = edge_data[get_edge_index(e, g)];
 441
 442  // Record the message. We'll need it in phase 3.
 443  edata.msg.assign(path);
 444
 445  // Note: "alpha" refers to the vertex "processor" receiving the
 446  // message.
 447  Vertex alpha = target(e, g);
 448
 449  switch (phase) {
 450  case 1:
 451    {
 452      num_edges_not_transmitted = degree(alpha, g) - 1;
 453      edata.is_tree_edge = true;
 454      parent = path.back();
 455      eta = parent;
 456      edata.M.clear(); edata.M.push_back(parent);
 457
 458      // Broadcast the path from the root to our potential children in
 459      // the spanning tree.
 460      path.push_back(alpha);
 461      path_header<Edge> header;
 462      header.path_length = path.size();
 463      ProcessGroup pg = process_group(g);
 464      BGL_FORALL_OUTEDGES_T(alpha, oe, g, Graph) {
 465        // Skip the tree edge we just received
 466        if (target(oe, g) != source(e, g)) {
 467          header.edge = oe;
 468          send(pg, get(owner, target(oe, g)), msg_path_header, header);
 469          send(pg, get(owner, target(oe, g)), msg_path_vertices, &path[0],
 470               header.path_length);
 471        }
 472      }
 473      path.pop_back();
 474
 475      // Swap the old path in, to save some extra copying. Nobody
 476      path_from_root.swap(path);
 477
 478      // Once we have received our place in the spanning tree, move on
 479      // to phase 2.
 480      phase = 2;
 481
 482      // If we only had only edge, skip to phase 3.
 483      if (num_edges_not_transmitted == 0)
 484        echo_phase(alpha, g);
 485      return;
 486    }
 487
 488  case 2:
 489    {
 490      --num_edges_not_transmitted;
 491      edata.is_tree_edge = false;
 492
 493      // Determine if alpha (our vertex) is in the path
 494      path_iterator pos = std::find(path.begin(), path.end(), alpha);
 495      if (pos != path.end()) {
 496        // Case A: There is a cycle alpha beta ... gamma alpha
 497        // M(e) <- {beta, gammar}
 498        edata.M.clear();
 499        ++pos;
 500        // If pos == path.end(), we have a self-loop
 501        if (pos != path.end()) {
 502          // Add beta
 503          edata.M.push_back(*pos);
 504          ++pos;
 505        }
 506        // If pos == path.end(), we have a self-loop or beta == gamma
 507        // (parallel edge). Otherwise, add gamma.
 508        if (pos != path.end()) edata.M.push_back(path.back());
 509      } else {
 510        // Case B: There is a cycle but we haven't seen alpha yet.
 511        // M(e) = {parent, path.back()}
 512        edata.M.clear();
 513        edata.M.push_back(path.back());
 514        if (parent != path.back()) edata.M.push_back(parent);
 515
 516        // eta = inf(eta, bra(pi_t, pi))
 517        eta = infimum(path_from_root, eta, branch_point(path_from_root, path));
 518      }
 519      if (num_edges_not_transmitted == 0)
 520        echo_phase(alpha, g);
 521      break;
 522    }
 523
 524  default:
 525//    std::cerr << "Phase is " << int(phase) << "\n";
 526    BOOST_ASSERT(false);
 527  }
 528}
 529
 530template<typename Graph>
 531void
 532hohberg_vertex_processor<Graph>::operator()(Edge e, Vertex gamma,
 533                                            path_t& in_same_bicomponent,
 534                                            const Graph& g)
 535{
 536  using namespace hohberg_detail;
 537
 538#ifdef PBGL_HOHBERG_DEBUG
 539  std::cerr << local(source(e, g)) << '@' << owner(source(e, g)) << " -> "
 540            << local(target(e, g)) << '@' << owner(target(e, g)) << ": tree("
 541            << local(gamma) << '@' << owner(gamma) << ", ";
 542  for (std::size_t i = 0; i < in_same_bicomponent.size(); ++i) {
 543    if (i > 0) std::cerr << ' ';
 544    std::cerr << local(in_same_bicomponent[i]) << '@'
 545              << owner(in_same_bicomponent[i]);
 546  }
 547  std::cerr << ", " << local(source(e, g)) << '@' << owner(source(e, g))
 548            << "), phase = " << (int)phase << std::endl;
 549#endif
 550
 551  // Get access to edge-specific data
 552  per_edge_data& edata = edge_data[get_edge_index(e, g)];
 553
 554  // Record the message. We'll need it in phase 3.
 555  edata.msg.assign(gamma, in_same_bicomponent);
 556
 557  // Note: "alpha" refers to the vertex "processor" receiving the
 558  // message.
 559  Vertex alpha = target(e, g);
 560  Vertex beta = source(e, g);
 561
 562  switch (phase) {
 563  case 2:
 564    --num_edges_not_transmitted;
 565    edata.is_tree_edge = true;
 566
 567    if (gamma == alpha) {
 568      // Case C
 569      edata.M.swap(in_same_bicomponent);
 570    } else {
 571      // Case D
 572      edata.M.clear();
 573      edata.M.push_back(parent);
 574      if (beta != parent) edata.M.push_back(beta);
 575      eta = infimum(path_from_root, eta, gamma);
 576    }
 577    if (num_edges_not_transmitted == 0)
 578      echo_phase(alpha, g);
 579    break;
 580
 581  default:
 582    BOOST_ASSERT(false);
 583  }
 584}
 585
 586template<typename Graph>
 587void
 588hohberg_vertex_processor<Graph>::operator()(Edge e, edges_size_type name,
 589                                            const Graph& g)
 590{
 591  using namespace hohberg_detail;
 592
 593#ifdef PBGL_HOHBERG_DEBUG
 594  std::cerr << local(source(e, g)) << '@' << owner(source(e, g)) << " -> "
 595            << local(target(e, g)) << '@' << owner(target(e, g)) << ": name("
 596            << name << "), phase = " << (int)phase << std::endl;
 597#endif
 598
 599  BOOST_ASSERT(phase == 4);
 600
 601  typename property_map<Graph, vertex_owner_t>::const_type
 602    owner = get(vertex_owner, g);
 603
 604  // Send name messages along the spanning tree edges that are in the
 605  // same bicomponent as the edge to our parent.
 606  ProcessGroup pg = process_group(g);
 607
 608  Vertex alpha = target(e, g);
 609
 610  std::size_t idx = 0;
 611  BGL_FORALL_OUTEDGES_T(alpha, e, g, Graph) {
 612    per_edge_data& edata = edge_data[idx++];
 613    if (edata.is_tree_edge
 614        && find(edata.M.begin(), edata.M.end(), parent) != edata.M.end()
 615        && target(e, g) != parent) {
 616      // Notify our children in the spanning tree of this name
 617      name_header<Edge> header;
 618      header.edge = e;
 619      header.name = name;
 620      send(pg, get(owner, target(e, g)), msg_name, header);
 621    } else if (target(e, g) == parent) {
 622      // Map from local partition numbers to global bicomponent numbers
 623      local_to_global_partitions[edata.partition] = name;
 624    }
 625  }
 626
 627  // Final stage
 628  phase = 5;
 629}
 630
 631template<typename Graph>
 632typename hohberg_vertex_processor<Graph>::edges_size_type
 633hohberg_vertex_processor<Graph>::
 634num_starting_bicomponents(Vertex alpha, const Graph& g)
 635{
 636  edges_size_type not_mapped = (std::numeric_limits<edges_size_type>::max)();
 637
 638  edges_size_type result = 0;
 639  std::size_t idx = 0;
 640  BGL_FORALL_OUTEDGES_T(alpha, e, g, Graph) {
 641    per_edge_data& edata = edge_data[idx++];
 642    if (edata.is_tree_edge
 643        && find(edata.M.begin(), edata.M.end(), parent) == edata.M.end()) {
 644      // Map from local partition numbers to global bicomponent numbers
 645      if (local_to_global_partitions[edata.partition] == not_mapped)
 646        local_to_global_partitions[edata.partition] = result++;
 647    }
 648  }
 649
 650#ifdef PBGL_HOHBERG_DEBUG
 651  std::cerr << local(alpha) << '@' << owner(alpha) << " has " << result
 652            << " bicomponents originating at it." << std::endl;
 653#endif
 654
 655  return result;
 656}
 657
 658template<typename Graph>
 659template<typename ComponentMap>
 660void
 661hohberg_vertex_processor<Graph>::
 662fill_edge_map(Vertex alpha, const Graph& g, ComponentMap& component)
 663{
 664  std::size_t idx = 0;
 665  BGL_FORALL_OUTEDGES_T(alpha, e, g, Graph) {
 666    per_edge_data& edata = edge_data[idx++];
 667    local_put(component, e, local_to_global_partitions[edata.partition]);
 668
 669#if defined(PBGL_HOHBERG_DEBUG) && PBGL_HOHBERG_DEBUG > 2
 670    std::cerr << "component("
 671              << local(source(e, g)) << '@' << owner(source(e, g)) << " -> "
 672              << local(target(e, g)) << '@' << owner(target(e, g)) << ") = "
 673              << local_to_global_partitions[edata.partition]
 674              << " (partition = " << edata.partition << " of "
 675              << local_to_global_partitions.size() << ")" << std::endl;
 676#endif
 677  }
 678}
 679
 680template<typename Graph>
 681void
 682hohberg_vertex_processor<Graph>::
 683start_naming_phase(Vertex alpha, const Graph& g, edges_size_type offset)
 684{
 685  using namespace hohberg_detail;
 686
 687  BOOST_ASSERT(phase == 4);
 688
 689  typename property_map<Graph, vertex_owner_t>::const_type
 690    owner = get(vertex_owner, g);
 691
 692  // Send name messages along the spanning tree edges of the
 693  // components that we get to number.
 694  ProcessGroup pg = process_group(g);
 695
 696  bool has_more_children_to_name = false;
 697
 698  // Map from local partition numbers to global bicomponent numbers
 699  edges_size_type not_mapped = (std::numeric_limits<edges_size_type>::max)();
 700  for (std::size_t i = 0; i < local_to_global_partitions.size(); ++i) {
 701    if (local_to_global_partitions[i] != not_mapped)
 702      local_to_global_partitions[i] += offset;
 703  }
 704
 705  std::size_t idx = 0;
 706  BGL_FORALL_OUTEDGES_T(alpha, e, g, Graph) {
 707    per_edge_data& edata = edge_data[idx++];
 708    if (edata.is_tree_edge
 709        && find(edata.M.begin(), edata.M.end(), parent) == edata.M.end()) {
 710      // Notify our children in the spanning tree of this new name
 711      name_header<Edge> header;
 712      header.edge = e;
 713      header.name = local_to_global_partitions[edata.partition];
 714      send(pg, get(owner, target(e, g)), msg_name, header);
 715    } else if (edata.is_tree_edge) {
 716      has_more_children_to_name = true;
 717    }
 718#if defined(PBGL_HOHBERG_DEBUG) && PBGL_HOHBERG_DEBUG > 2
 719    std::cerr << "M[" << local(source(e, g)) << '@' << owner(source(e, g))
 720              << " -> " << local(target(e, g)) << '@' << owner(target(e, g))
 721              << "] = ";
 722    for (std::size_t i = 0; i < edata.M.size(); ++i) {
 723      std::cerr << local(edata.M[i]) << '@' << owner(edata.M[i]) << ' ';
 724    }
 725    std::cerr << std::endl;
 726#endif
 727  }
 728
 729  // See if we're done.
 730  if (!has_more_children_to_name)
 731    // Final stage
 732    phase = 5;
 733}
 734
 735template<typename Graph>
 736void
 737hohberg_vertex_processor<Graph>::echo_phase(Vertex alpha, const Graph& g)
 738{
 739  using namespace hohberg_detail;
 740
 741  typename property_map<Graph, vertex_owner_t>::const_type
 742    owner = get(vertex_owner, g);
 743
 744  /* We're entering the echo phase. */
 745  phase = 3;
 746
 747  if (parent != graph_traits<Graph>::null_vertex()) {
 748    Edge edge_to_parent;
 749
 750#if defined(PBGL_HOHBERG_DEBUG) && PBGL_HOHBERG_DEBUG > 1
 751     std::cerr << local(alpha) << '@' << owner(alpha) << " echo: parent = "
 752               << local(parent) << '@' << owner(parent) << ", eta = "
 753               << local(eta) << '@' << owner(eta) << ", Gamma = ";
 754#endif
 755
 756    std::vector<Vertex> bicomp;
 757    std::size_t e_index = 0;
 758    BGL_FORALL_OUTEDGES_T(alpha, e, g, Graph) {
 759      if (target(e, g) == parent && parent == eta) {
 760        edge_to_parent = e;
 761        if (find(bicomp.begin(), bicomp.end(), alpha) == bicomp.end()) {
 762#if defined(PBGL_HOHBERG_DEBUG) && PBGL_HOHBERG_DEBUG > 1
 763          std::cerr << local(alpha) << '@' << owner(alpha) << ' ';
 764#endif
 765          bicomp.push_back(alpha);
 766        }
 767      } else {
 768        if (target(e, g) == parent) edge_to_parent = e;
 769
 770        per_edge_data& edata = edge_data[e_index];
 771
 772        if (edata.msg.is_path()) {
 773          path_iterator pos = std::find(edata.msg.path_or_bicomp.begin(),
 774                                        edata.msg.path_or_bicomp.end(),
 775                                        eta);
 776          if (pos != edata.msg.path_or_bicomp.end()) {
 777            ++pos;
 778            if (pos != edata.msg.path_or_bicomp.end()
 779                && find(bicomp.begin(), bicomp.end(), *pos) == bicomp.end()) {
 780#if defined(PBGL_HOHBERG_DEBUG) && PBGL_HOHBERG_DEBUG > 1
 781              std::cerr << local(*pos) << '@' << owner(*pos) << ' ';
 782#endif
 783              bicomp.push_back(*pos);
 784            }
 785          }
 786        } else if (edata.msg.is_tree() && edata.msg.gamma == eta) {
 787          for (path_iterator i = edata.msg.path_or_bicomp.begin();
 788               i != edata.msg.path_or_bicomp.end(); ++i) {
 789            if (find(bicomp.begin(), bicomp.end(), *i) == bicomp.end()) {
 790#if defined(PBGL_HOHBERG_DEBUG) && PBGL_HOHBERG_DEBUG > 1
 791              std::cerr << local(*i) << '@' << owner(*i) << ' ';
 792#endif
 793              bicomp.push_back(*i);
 794            }
 795          }
 796        }
 797      }
 798      ++e_index;
 799    }
 800#ifdef PBGL_HOHBERG_DEBUG
 801    std::cerr << std::endl;
 802#endif
 803
 804    // Send tree(eta, bicomp) to parent
 805    tree_header<Vertex, Edge> header;
 806    header.edge = edge_to_parent;
 807    header.gamma = eta;
 808    header.bicomp_length = bicomp.size();
 809    ProcessGroup pg = process_group(g);
 810    send(pg, get(owner, parent), msg_tree_header, header);
 811    send(pg, get(owner, parent), msg_tree_vertices, &bicomp[0],
 812         header.bicomp_length);
 813  }
 814
 815  // Compute the partition of edges such that iff two edges e1 and e2
 816  // are in different subsets then M(e1) is disjoint from M(e2).
 817
 818  // Start by putting each edge in a different partition
 819  std::vector<degree_size_type> parent_vec(edge_data.size());
 820  degree_size_type idx = 0;
 821  for (idx = 0; idx < edge_data.size(); ++idx)
 822    parent_vec[idx] = idx;
 823
 824  // Go through each edge e, performing a union() on the edges
 825  // incident on all vertices in M[e].
 826  idx = 0;
 827  BGL_FORALL_OUTEDGES_T(alpha, e, g, Graph) {
 828    per_edge_data& edata = edge_data[idx++];
 829
 830    // Compute union of vertices in M
 831    if (!edata.M.empty()) {
 832      degree_size_type e1 = get_incident_edge_index(alpha, edata.M.front(), g);
 833      while (parent_vec[e1] != e1) e1 = parent_vec[e1];
 834
 835      for (std::size_t i = 1; i < edata.M.size(); ++i) {
 836        degree_size_type e2 = get_incident_edge_index(alpha, edata.M[i], g);
 837        while (parent_vec[e2] != e2) e2 = parent_vec[e2];
 838        parent_vec[e2] = e1;
 839      }
 840    }
 841  }
 842
 843  edges_size_type not_mapped = (std::numeric_limits<edges_size_type>::max)();
 844
 845  // Determine the number of partitions
 846  for (idx = 0; idx < parent_vec.size(); ++idx) {
 847    if (parent_vec[idx] == idx) {
 848      edge_data[idx].partition = local_to_global_partitions.size();
 849      local_to_global_partitions.push_back(not_mapped);
 850    }
 851  }
 852
 853  // Assign partition numbers to each edge
 854  for (idx = 0; idx < parent_vec.size(); ++idx) {
 855    degree_size_type rep = parent_vec[idx];
 856    while (rep != parent_vec[rep]) rep = parent_vec[rep];
 857    edge_data[idx].partition = edge_data[rep].partition;
 858  }
 859
 860  // Enter the naming phase (but don't send anything yet).
 861  phase = 4;
 862}
 863
 864template<typename Graph>
 865std::size_t
 866hohberg_vertex_processor<Graph>::get_edge_index(Edge e, const Graph& g)
 867{
 868  std::size_t result = 0;
 869  BGL_FORALL_OUTEDGES_T(target(e, g), oe, g, Graph) {
 870    if (source(e, g) == target(oe, g)) return result;
 871    ++result;
 872  }
 873  BOOST_ASSERT(false);
 874}
 875
 876template<typename Graph>
 877std::size_t
 878hohberg_vertex_processor<Graph>::get_incident_edge_index(Vertex u, Vertex v,
 879                                                         const Graph& g)
 880{
 881  std::size_t result = 0;
 882  BGL_FORALL_OUTEDGES_T(u, e, g, Graph) {
 883    if (target(e, g) == v) return result;
 884    ++result;
 885  }
 886  BOOST_ASSERT(false);
 887}
 888
 889template<typename Graph, typename InputIterator, typename ComponentMap,
 890         typename VertexProcessorMap>
 891typename graph_traits<Graph>::edges_size_type
 892hohberg_biconnected_components
 893  (const Graph& g,
 894   ComponentMap component,
 895   InputIterator first, InputIterator last,
 896   VertexProcessorMap vertex_processor)
 897{
 898  using namespace boost::graph::parallel;
 899  using namespace hohberg_detail;
 900  using boost::parallel::all_reduce;
 901
 902  typename property_map<Graph, vertex_owner_t>::const_type
 903    owner = get(vertex_owner, g);
 904
 905  // The graph must be undirected
 906  BOOST_STATIC_ASSERT(
 907    (is_convertible<typename graph_traits<Graph>::directed_category,
 908                    undirected_tag>::value));
 909
 910  // The graph must model Incidence Graph
 911  function_requires< IncidenceGraphConcept<Graph> >();
 912
 913  typedef typename graph_traits<Graph>::edges_size_type edges_size_type;
 914  typedef typename graph_traits<Graph>::degree_size_type degree_size_type;
 915  typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;
 916  typedef typename graph_traits<Graph>::edge_descriptor edge_descriptor;
 917
 918  // Retrieve the process group we will use for communication
 919  typedef typename process_group_type<Graph>::type process_group_type;
 920  process_group_type pg = process_group(g);
 921
 922  // Keeps track of the edges that we know to be tree edges.
 923  std::vector<edge_descriptor> tree_edges;
 924
 925  // The leaders send out a path message to initiate the algorithm
 926  while (first != last) {
 927    vertex_descriptor leader = *first;
 928    if (process_id(pg) == get(owner, leader))
 929      vertex_processor[leader].initialize_leader(leader, g);
 930    ++first;
 931  }
 932  synchronize(pg);
 933
 934  // Will hold the number of bicomponents in the graph.
 935  edges_size_type num_bicomponents = 0;
 936
 937  // Keep track of the path length that we should expect, based on the
 938  // level in the breadth-first search tree. At present, this is only
 939  // used as a sanity check. TBD: This could be used to decrease the
 940  // amount of communication required per-edge (by about 4 bytes).
 941  std::size_t path_length = 1;
 942
 943  typedef std::vector<vertex_descriptor> path_t;
 944  typedef typename path_t::iterator path_iterator;
 945
 946  unsigned char minimum_phase = 5;
 947  do {
 948    while (optional<std::pair<int, int> > msg = probe(pg)) {
 949      switch (msg->second) {
 950      case msg_path_header:
 951        {
 952          // Receive the path header
 953          path_header<edge_descriptor> header;
 954          receive(pg, msg->first, msg->second, header);
 955          BOOST_ASSERT(path_length == header.path_length);
 956
 957          // Receive the path itself
 958          path_t path(path_length);
 959          receive(pg, msg->first, msg_path_vertices, &path[0], path_length);
 960
 961          edge_descriptor e = header.edge;
 962          vertex_processor[target(e, g)](e, path, g);
 963        }
 964        break;
 965
 966      case msg_path_vertices:
 967        // Should be handled in msg_path_header case, unless we're going
 968        // stateless.
 969        BOOST_ASSERT(false);
 970        break;
 971
 972      case msg_tree_header:
 973        {
 974          // Receive the tree header
 975          tree_header<vertex_descriptor, edge_descriptor> header;
 976          receive(pg, msg->first, msg->second, header);
 977
 978          // Receive the tree itself
 979          path_t in_same_bicomponent(header.bicomp_length);
 980          receive(pg, msg->first, msg_tree_vertices, &in_same_bicomponent[0],
 981                  header.bicomp_length);
 982
 983          edge_descriptor e = header.edge;
 984          vertex_processor[target(e, g)](e, header.gamma, in_same_bicomponent,
 985                                         g);
 986        }
 987        break;
 988
 989      case msg_tree_vertices:
 990        // Should be handled in msg_tree_header case, unless we're
 991        // going stateless.
 992        BOOST_ASSERT(false);
 993        break;
 994
 995      case msg_name:
 996        {
 997          name_header<edge_descriptor> header;
 998          receive(pg, msg->first, msg->second, header);
 999          edge_descriptor e = header.edge;
1000          vertex_processor[target(e, g)](e, header.name, g);
1001        }
1002        break;
1003
1004      default:
1005        BOOST_ASSERT(false);
1006      }
1007    }
1008    ++path_length;
1009
1010    // Compute minimum phase locally
1011    minimum_phase = 5;
1012    unsigned char maximum_phase = 1;
1013    BGL_FORALL_VERTICES_T(v, g, Graph) {
1014      minimum_phase = (std::min)(minimum_phase, vertex_processor[v].get_phase());
1015      maximum_phase = (std::max)(maximum_phase, vertex_processor[v].get_phase());
1016    }
1017
1018#ifdef PBGL_HOHBERG_DEBUG
1019    if (process_id(pg) == 0)
1020      std::cerr << "<---------End of stage------------->" << std::endl;
1021#endif
1022    // Compute minimum phase globally
1023    minimum_phase = all_reduce(pg, minimum_phase, boost::mpi::minimum<char>());
1024
1025#ifdef PBGL_HOHBERG_DEBUG
1026    if (process_id(pg) == 0)
1027      std::cerr << "Minimum phase = " << (int)minimum_phase << std::endl;
1028#endif
1029
1030    if (minimum_phase == 4
1031        && all_reduce(pg, maximum_phase, boost::mpi::maximum<char>()) == 4) {
1032
1033#ifdef PBGL_HOHBERG_DEBUG
1034      if (process_id(pg) == 0)
1035        std::cerr << "<---------Naming phase------------->" << std::endl;
1036#endif
1037      // Compute the biconnected component number offsets for each
1038      // vertex.
1039      std::vector<edges_size_type> local_offsets;
1040      local_offsets.reserve(num_vertices(g));
1041      edges_size_type num_local_bicomponents = 0;
1042      BGL_FORALL_VERTICES_T(v, g, Graph) {
1043        local_offsets.push_back(num_local_bicomponents);
1044        num_local_bicomponents +=
1045          vertex_processor[v].num_starting_bicomponents(v, g);
1046      }
1047
1048      synchronize(pg);
1049
1050      // Find our the number of bicomponent names that will originate
1051      // from each process. This tells us how many bicomponents are in
1052      // the entire graph and what our global offset is for computing
1053      // our own biconnected component names.
1054      std::vector<edges_size_type> all_bicomponents(num_processes(pg));
1055      all_gather(pg, &num_local_bicomponents, &num_local_bicomponents + 1,
1056                 all_bicomponents);
1057      num_bicomponents = 0;
1058      edges_size_type my_global_offset = 0;
1059      for (std::size_t i = 0; i < all_bicomponents.size(); ++i) {
1060        if (i == (std::size_t)process_id(pg)) 
1061          my_global_offset = num_bicomponents;
1062        num_bicomponents += all_bicomponents[i];
1063      }
1064
1065      std::size_t index = 0;
1066      BGL_FORALL_VERTICES_T(v, g, Graph) {
1067        edges_size_type offset = my_global_offset + local_offsets[index++];
1068        vertex_processor[v].start_naming_phase(v, g, offset);
1069      }
1070    }
1071
1072    synchronize(pg);
1073  } while (minimum_phase < 5);
1074
1075  // Number the edges appropriately.
1076  BGL_FORALL_VERTICES_T(v, g, Graph)
1077    vertex_processor[v].fill_edge_map(v, g, component);
1078
1079  return num_bicomponents;
1080}
1081
1082template<typename Graph, typename ComponentMap, typename InputIterator>
1083typename graph_traits<Graph>::edges_size_type
1084hohberg_biconnected_components
1085  (const Graph& g, ComponentMap component,
1086   InputIterator first, InputIterator last)
1087
1088{
1089  std::vector<hohberg_vertex_processor<Graph> >
1090    vertex_processors(num_vertices(g));
1091  return hohberg_biconnected_components
1092           (g, component, first, last,
1093            make_iterator_property_map(vertex_processors.begin(),
1094                                       get(vertex_index, g)));
1095}
1096
1097template<typename Graph, typename ComponentMap, typename ParentMap>
1098typename graph_traits<Graph>::edges_size_type
1099hohberg_biconnected_components(const Graph& g, ComponentMap component,
1100                               ParentMap parent)
1101{
1102  // We need the connected components of the graph, but we don't care
1103  // about component numbers.
1104  connected_components(g, dummy_property_map(), parent);
1105
1106  // Each root in the parent map is a leader
1107  typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;
1108  std::vector<vertex_descriptor> leaders;
1109  BGL_FORALL_VERTICES_T(v, g, Graph)
1110    if (get(parent, v) == v) leaders.push_back(v);
1111
1112  return hohberg_biconnected_components(g, component,
1113                                        leaders.begin(), leaders.end());
1114}
1115
1116template<typename Graph, typename ComponentMap>
1117typename graph_traits<Graph>::edges_size_type
1118hohberg_biconnected_components(const Graph& g, ComponentMap component)
1119{
1120  typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;
1121  std::vector<vertex_descriptor> parents(num_vertices(g));
1122  return hohberg_biconnected_components
1123           (g, component, make_iterator_property_map(parents.begin(),
1124                                                     get(vertex_index, g)));
1125}
1126
1127} } } // end namespace boost::graph::distributed
1128
1129#endif // BOOST_GRAPH_DISTRIBUTED_HOHBERG_BICONNECTED_COMPONENTS_HPP