graph_parallel/include/boost/graph/distributed/hohberg_biconnected_components.hpp
Marcel Raad f3a3284f42 Fix for MSVC (#9)
MSVC doesn't know "and". It warns:
hohberg_biconnected_components.hpp(166) : warning C4067: unexpected tokens following preprocessor directive - expected a newline
hohberg_biconnected_components.hpp(201) : warning C4067: unexpected tokens following preprocessor directive - expected a newline
2016-05-01 15:29:13 -06:00

1129 lines
37 KiB
C++

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