/Src/Dependencies/Boost/boost/graph/distributed/hohberg_biconnected_components.hpp
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