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