boost/graph/distributed/strong_components.hpp

   1 // Copyright (C) 2004-2008 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: Nick Edmonds
   8 //           Douglas Gregor
   9 //           Andrew Lumsdaine
  10 #ifndef BOOST_GRAPH_DISTRIBUTED_SCC_HPP
  11 #define BOOST_GRAPH_DISTRIBUTED_SCC_HPP
  12
  13 #ifndef BOOST_GRAPH_USE_MPI
  14 #error "Parallel BGL files should not be included unless <boost/graph/use_mpi.hpp> has been included"
  15 #endif
  16
  17 // #define PBGL_SCC_DEBUG
  18
  19 #include <boost/assert.hpp>
  20 #include <boost/property_map/property_map.hpp>
  21 #include <boost/property_map/parallel/distributed_property_map.hpp>
  22 #include <boost/property_map/parallel/caching_property_map.hpp>
  23 #include <boost/graph/parallel/algorithm.hpp>
  24 #include <boost/graph/parallel/process_group.hpp>
  25 #include <boost/graph/distributed/queue.hpp>
  26 #include <boost/graph/distributed/filtered_graph.hpp>
  27 #include <boost/pending/indirect_cmp.hpp>
  28 #include <boost/graph/breadth_first_search.hpp>
  29 #include <boost/graph/graph_traits.hpp>
  30 #include <boost/graph/overloading.hpp>
  31 #include <boost/graph/distributed/concepts.hpp>
  32 #include <boost/graph/distributed/local_subgraph.hpp>
  33 #include <boost/graph/parallel/properties.hpp>
  34 #include <boost/graph/named_function_params.hpp>
  35 #include <boost/graph/random.hpp>
  36 #include <boost/graph/distributed/reverse_graph.hpp>
  37 #include <boost/optional.hpp>
  38 #include <boost/graph/distributed/detail/filtered_queue.hpp>
  39 #include <boost/graph/distributed/adjacency_list.hpp>
  40 #ifdef PBGL_SCC_DEBUG
  41   #include <iostream>
  42   #include <cstdlib>
  43   #include <iomanip>
  44   #include <sys/time.h>
  45   #include <boost/graph/distributed/graphviz.hpp> // for ostringstream
  46 #endif
  47 #include <vector>
  48 #include <map>
  49 #include <boost/graph/parallel/container_traits.hpp>
  50
  51 #ifdef PBGL_SCC_DEBUG
  52 #  include <boost/graph/accounting.hpp>
  53 #endif /* PBGL_SCC_DEBUG */
  54
  55 // If your graph is likely to have large numbers of small strongly connected
  56 // components then running the sequential SCC algorithm on the local subgraph
  57 // and filtering components with no remote edges may increase performance
  58 // #define FILTER_LOCAL_COMPONENTS
  59
  60 namespace boost { namespace graph { namespace distributed { namespace detail {
  61
  62 template<typename vertex_descriptor>
  63 struct v_sets{
  64   std::vector<vertex_descriptor> pred, succ, intersect, ps_union;
  65 };
  66
  67 /* Serialize vertex set */
  68 template<typename Graph>
  69 void
  70 marshal_set( std::vector<std::vector<typename graph_traits<Graph>::vertex_descriptor> > in,
  71              std::vector<typename graph_traits<Graph>::vertex_descriptor>& out )
  72 {
  73   for( std::size_t i = 0; i < in.size(); ++i ) {
  74     out.insert( out.end(), graph_traits<Graph>::null_vertex() );
  75     out.insert( out.end(), in[i].begin(), in[i].end() );
  76   }
  77 }
  78
  79 /* Un-serialize vertex set */
  80 template<typename Graph>
  81 void
  82 unmarshal_set( std::vector<typename graph_traits<Graph>::vertex_descriptor> in,
  83                std::vector<std::vector<typename graph_traits<Graph>::vertex_descriptor> >& out )
  84 {
  85   typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;
  86
  87   while( !in.empty() ) {
  88     typename std::vector<vertex_descriptor>::iterator end
  89       = std::find( in.begin(), in.end(), graph_traits<Graph>::null_vertex() );
  90
  91     if( end == in.begin() )
  92       in.erase( in.begin() );
  93     else {
  94       out.push_back(std::vector<vertex_descriptor>());
  95       out[out.size() - 1].insert( out[out.size() - 1].end(), in.begin(), end );
  96       in.erase( in.begin(), end );
  97     }
  98   }
  99 }
 100
 101 /* Determine if vertex is in subset */
 102 template <typename Set>
 103 struct in_subset {
 104   in_subset() : m_s(0) { }
 105   in_subset(const Set& s) : m_s(&s) { }
 106
 107   template <typename Elt>
 108   bool operator()(const Elt& x) const {
 109     return ((*m_s).find(x) != (*m_s).end());
 110   }
 111
 112 private:
 113   const Set* m_s;
 114 };
 115
 116 template<typename T>
 117 struct vertex_identity_property_map
 118   : public boost::put_get_helper<T, vertex_identity_property_map<T> >
 119 {
 120   typedef T key_type;
 121   typedef T value_type;
 122   typedef T reference;
 123   typedef boost::readable_property_map_tag category;
 124
 125   inline value_type operator[](const key_type& v) const { return v; }
 126   inline void clear() { }
 127 };
 128
 129 template <typename T>
 130 inline void synchronize( vertex_identity_property_map<T> & ) { }
 131
 132 /* BFS visitor for SCC */
 133 template<typename Graph, typename SourceMap>
 134 struct scc_discovery_visitor : bfs_visitor<>
 135 {
 136   scc_discovery_visitor(SourceMap& sourceM)
 137     : sourceM(sourceM) {}
 138
 139   template<typename Edge>
 140   void tree_edge(Edge e, const Graph& g)
 141   {
 142     put(sourceM, target(e,g), get(sourceM, source(e,g)));
 143   }
 144
 145  private:
 146   SourceMap& sourceM;
 147 };
 148
 149 } } } } /* End namespace boost::graph::distributed::detail */
 150
 151 namespace boost { namespace graph { namespace distributed {
 152     enum fhp_message_tags { fhp_edges_size_msg, fhp_add_edges_msg, fhp_pred_size_msg,
 153                             fhp_pred_msg, fhp_succ_size_msg, fhp_succ_msg };
 154
 155     template<typename Graph, typename ReverseGraph,
 156              typename VertexComponentMap, typename IsoMapFR, typename IsoMapRF,
 157              typename VertexIndexMap>
 158     void
 159     fleischer_hendrickson_pinar_strong_components(const Graph& g,
 160                                                   VertexComponentMap c,
 161                                                   const ReverseGraph& gr,
 162                                                   IsoMapFR fr, IsoMapRF rf,
 163                                                   VertexIndexMap vertex_index_map)
 164     {
 165       typedef typename graph_traits<Graph>::vertex_iterator vertex_iterator;
 166       typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;
 167       typedef typename graph_traits<ReverseGraph>::vertex_iterator rev_vertex_iterator;
 168       typedef typename graph_traits<ReverseGraph>::vertex_descriptor rev_vertex_descriptor;
 169       typedef typename boost::graph::parallel::process_group_type<Graph>::type
 170         process_group_type;
 171       typedef typename process_group_type::process_id_type process_id_type;
 172       typedef iterator_property_map<typename std::vector<vertex_descriptor>::iterator,
 173                                     VertexIndexMap> ParentMap;
 174       typedef iterator_property_map<typename std::vector<default_color_type>::iterator,
 175                                     VertexIndexMap> ColorMap;
 176       typedef iterator_property_map<typename std::vector<vertex_descriptor>::iterator,
 177                                     VertexIndexMap> Rev_ParentMap;
 178       typedef std::vector<std::pair<vertex_descriptor, vertex_descriptor> > VertexPairVec;
 179
 180       typedef typename property_map<Graph, vertex_owner_t>::const_type
 181         OwnerMap;
 182
 183       OwnerMap owner = get(vertex_owner, g);
 184
 185       using boost::graph::parallel::process_group;
 186       process_group_type pg = process_group(g);
 187       process_id_type id = process_id(pg);
 188       int num_procs = num_processes(pg);
 189       int n = 0;
 190
 191       int my_n = num_vertices(g);
 192       all_reduce(pg, &my_n, &my_n+1, &n, std::plus<int>());
 193
 194       //
 195       // Initialization
 196       //
 197
 198 #ifdef PBGL_SCC_DEBUG
 199   accounting::time_type start = accounting::get_time();
 200 #endif
 201
 202       vertex_iterator vstart, vend;
 203       rev_vertex_iterator rev_vstart, rev_vend;
 204       std::vector<std::vector<vertex_descriptor> > vertex_sets, new_vertex_sets;
 205
 206       vertex_sets.push_back(std::vector<vertex_descriptor>());
 207
 208       // Remove vertices that do not have at least one in edge and one out edge
 209       new_vertex_sets.push_back(std::vector<vertex_descriptor>());
 210       for( boost::tie(vstart, vend) = vertices(g); vstart != vend; vstart++ )
 211         if( out_degree( get(fr, *vstart), gr) > 0 && out_degree(*vstart, g) > 0 )
 212           new_vertex_sets[0].push_back( *vstart );
 213
 214       // Perform sequential SCC on local subgraph, filter all components with external
 215       // edges, mark remaining components and remove them from vertex_sets
 216 #ifdef FILTER_LOCAL_COMPONENTS
 217       // This doesn't actually speed up SCC in connected graphs it seems, but it does work
 218       // and may help in the case where there are lots of small strong components.
 219       {
 220         local_subgraph<const Graph> ls(g);
 221         typedef typename property_map<local_subgraph<const Graph>, vertex_index_t>::type
 222           local_index_map_type;
 223         local_index_map_type local_index = get(vertex_index, ls);
 224
 225         std::vector<int> ls_components_vec(num_vertices(ls));
 226         typedef iterator_property_map<std::vector<int>::iterator, local_index_map_type>
 227           ls_components_map_type;
 228         ls_components_map_type ls_component(ls_components_vec.begin(), local_index);
 229         int num_comp = boost::strong_components(ls, ls_component);
 230
 231         // Create map of components
 232         std::map<int, std::vector<vertex_descriptor> > local_comp_map;
 233         typedef typename graph_traits<local_subgraph<const Graph> >::vertex_iterator ls_vertex_iterator;
 234         ls_vertex_iterator vstart, vend;
 235         for( boost::tie(vstart,vend) = vertices(ls); vstart != vend; vstart++ )
 236           local_comp_map[get(ls_component, *vstart)].push_back( *vstart );
 237
 238         // Filter components that have no non-local edges
 239         typedef typename graph_traits<Graph>::adjacency_iterator adjacency_iterator;
 240         typedef typename graph_traits<ReverseGraph>::adjacency_iterator rev_adjacency_iterator;
 241         adjacency_iterator abegin, aend;
 242         rev_adjacency_iterator rev_abegin, rev_aend;
 243         for( std::size_t i = 0; i < num_comp; ++i ) {
 244           bool local = true;
 245           for( std::size_t j = 0; j < local_comp_map[i].size(); j++ ) {
 246             for( boost::tie(abegin,aend) = adjacent_vertices(local_comp_map[i][j], g);
 247                  abegin != aend; abegin++ )
 248               if( get(owner, *abegin) != id ) {
 249                 local = false;
 250                 break;
 251               }
 252
 253             if( local )
 254               for( boost::tie(rev_abegin,rev_aend) = adjacent_vertices(get(fr, local_comp_map[i][j]), gr);
 255                    rev_abegin != rev_aend; rev_abegin++ )
 256                 if( get(owner, *rev_abegin) != id ) {
 257                   local = false;
 258                   break;
 259                 }
 260
 261             if( !local ) break;
 262           }
 263
 264           if( local ) // Mark and remove from new_vertex_sets
 265             for( std::size_t j = 0; j < local_comp_map[i].size(); j++ ) {
 266               put( c, local_comp_map[i][j], local_comp_map[i][0] );
 267               typename std::vector<vertex_descriptor>::iterator pos =
 268                 std::find(new_vertex_sets[0].begin(), new_vertex_sets[0].end(), local_comp_map[i][j]);
 269               if( pos != new_vertex_sets[0].end() )
 270                 new_vertex_sets[0].erase(pos);
 271             }
 272         }
 273       }
 274 #endif // FILTER_LOCAL_COMPONENTS
 275
 276       all_gather( pg, new_vertex_sets[0].begin(), new_vertex_sets[0].end(), vertex_sets[0] );
 277       new_vertex_sets.clear();
 278
 279 #ifdef PBGL_SCC_DEBUG
 280   accounting::time_type end = accounting::get_time();
 281   if(id == 0)
 282     std::cerr << "Trim local SCCs time = " << accounting::print_time(end - start) << " seconds.\n";
 283 #endif
 284
 285       if( vertex_sets[0].empty() ) return;
 286
 287       //
 288       // Recursively determine SCCs
 289       //
 290
 291 #ifdef PBGL_SCC_DEBUG
 292   int iterations = 0;
 293 #endif
 294
 295       // Only need to be able to map starting vertices for BFS from now on
 296       fr.clear();
 297
 298       do {
 299
 300 #ifdef PBGL_SCC_DEBUG
 301   if(id == 0) {
 302     printf("\n\nIteration %d:\n\n", iterations++);
 303
 304     if( iterations > 1 ) {
 305       end = accounting::get_time();
 306       std::cerr << "Running main loop destructors time = " << accounting::print_time(end - start) << " seconds.\n";
 307     }
 308
 309     start = accounting::get_time();
 310   }
 311 #endif
 312
 313         // Get forward->reverse mappings for BFS start vertices
 314        for (std::size_t i = 0; i < vertex_sets.size(); ++i)
 315           get(fr, vertex_sets[i][0]);
 316         synchronize(fr);
 317
 318         // Determine local vertices to start BFS from
 319         std::vector<vertex_descriptor> local_start;
 320         for( std::size_t i = id; i < vertex_sets.size(); i += num_procs )
 321           local_start.push_back(vertex_sets[i][0]);
 322
 323         if( local_start.empty() )
 324           local_start.push_back(vertex_sets[0][0]);
 325
 326
 327         // Make filtered graphs
 328         typedef std::set<vertex_descriptor> VertexSet;
 329         typedef std::set<rev_vertex_descriptor> Rev_VertexSet;
 330
 331         VertexSet filter_set_g;
 332         Rev_VertexSet filter_set_gr;
 333         typename VertexSet::iterator fs;
 334
 335         int active_vertices = 0;
 336         for (std::size_t i = 0; i < vertex_sets.size(); i++)
 337           active_vertices += vertex_sets[i].size();
 338
 339         // This is a completely random bound
 340         if ( active_vertices < 0.05*n ) {
 341           // TODO: This set insertion is ridiculously inefficient, make it an in-place-merge?
 342           for (std::size_t i = 0; i < vertex_sets.size(); i++)
 343             filter_set_g.insert(vertex_sets[i].begin(), vertex_sets[i].end());
 344
 345           for (fs = filter_set_g.begin(); fs != filter_set_g.end(); ++fs )
 346             filter_set_gr.insert(get(fr, *fs));
 347         }
 348
 349         filtered_graph<const Graph, keep_all, detail::in_subset<VertexSet> >
 350           fg(g, keep_all(), detail::in_subset<VertexSet>(filter_set_g));
 351
 352         filtered_graph<const ReverseGraph, keep_all, detail::in_subset<VertexSet> >
 353           fgr(gr, keep_all(), detail::in_subset<VertexSet>(filter_set_gr));
 354
 355         // Add additional starting vertices to BFS queue
 356         typedef filtered_queue<queue<vertex_descriptor>, boost::detail::has_not_been_seen<VertexIndexMap> >
 357           local_queue_type;
 358         typedef boost::graph::distributed::distributed_queue<process_group_type, OwnerMap, local_queue_type>
 359           queue_t;
 360
 361         typedef typename property_map<ReverseGraph, vertex_owner_t>::const_type
 362           RevOwnerMap;
 363
 364         typedef filtered_queue<queue<rev_vertex_descriptor>, boost::detail::has_not_been_seen<VertexIndexMap> >
 365           rev_local_queue_type;
 366         typedef boost::graph::distributed::distributed_queue<process_group_type, RevOwnerMap, rev_local_queue_type>
 367           rev_queue_t;
 368
 369         queue_t Q(process_group(g),
 370                   owner,
 371                   make_filtered_queue(queue<vertex_descriptor>(),
 372                                       boost::detail::has_not_been_seen<VertexIndexMap>
 373                                       (num_vertices(g), vertex_index_map)),
 374                   false);
 375
 376         rev_queue_t Qr(process_group(gr),
 377                        get(vertex_owner, gr),
 378                        make_filtered_queue(queue<rev_vertex_descriptor>(),
 379                                            boost::detail::has_not_been_seen<VertexIndexMap>
 380                                            (num_vertices(gr), vertex_index_map)),
 381                        false);
 382
 383         for( std::size_t i = 1; i < local_start.size(); ++i ) {
 384           Q.push(local_start[i]);
 385           Qr.push(get(fr, local_start[i]));
 386         }
 387
 388 #ifdef PBGL_SCC_DEBUG
 389   end = accounting::get_time();
 390   if(id == 0)
 391     std::cerr << "  Initialize BFS time = " << accounting::print_time(end - start) << " seconds.\n";
 392   start = accounting::get_time();
 393 #endif
 394
 395 #ifdef PBGL_SCC_DEBUG
 396   accounting::time_type start2 = accounting::get_time();
 397 #endif
 398
 399         // Forward BFS
 400         std::vector<default_color_type> color_map_s(num_vertices(g));
 401         ColorMap color_map(color_map_s.begin(), vertex_index_map);
 402         std::vector<vertex_descriptor> succ_map_s(num_vertices(g), graph_traits<Graph>::null_vertex());
 403         ParentMap succ_map(succ_map_s.begin(), vertex_index_map);
 404
 405         for( std::size_t i = 0; i < vertex_sets.size(); ++i )
 406           put(succ_map, vertex_sets[i][0], vertex_sets[i][0]);
 407
 408 #ifdef PBGL_SCC_DEBUG
 409   accounting::time_type end2 = accounting::get_time();
 410   if(id == 0)
 411     std::cerr << "  Initialize forward BFS time = " << accounting::print_time(end2 - start2) << " seconds.\n";
 412 #endif
 413
 414         if (active_vertices < 0.05*n)
 415           breadth_first_search(fg, local_start[0], Q,
 416                                detail::scc_discovery_visitor<filtered_graph<const Graph, keep_all,
 417                                                                             detail::in_subset<VertexSet> >, ParentMap>
 418                                (succ_map),
 419                                color_map);
 420         else
 421           breadth_first_search(g, local_start[0], Q,
 422                                detail::scc_discovery_visitor<const Graph, ParentMap>(succ_map),
 423                                color_map);
 424
 425 #ifdef PBGL_SCC_DEBUG
 426   start2 = accounting::get_time();
 427 #endif
 428
 429         // Reverse BFS
 430         color_map.clear(); // reuse color map since g and gr have same vertex index
 431         std::vector<vertex_descriptor> pred_map_s(num_vertices(gr), graph_traits<Graph>::null_vertex());
 432         Rev_ParentMap pred_map(pred_map_s.begin(), vertex_index_map);
 433
 434         for( std::size_t i = 0; i < vertex_sets.size(); ++i )
 435           put(pred_map, get(fr, vertex_sets[i][0]), vertex_sets[i][0]);
 436
 437 #ifdef PBGL_SCC_DEBUG
 438   end2 = accounting::get_time();
 439   if(id == 0)
 440     std::cerr << "  Initialize reverse BFS time = " << accounting::print_time(end2 - start2) << " seconds.\n";
 441 #endif
 442
 443         if (active_vertices < 0.05*n)
 444           breadth_first_search(fgr, get(fr, local_start[0]),
 445                                Qr,
 446                                detail::scc_discovery_visitor<filtered_graph<const ReverseGraph, keep_all,
 447                                                                             detail::in_subset<Rev_VertexSet> >, Rev_ParentMap>
 448                                (pred_map),
 449                                color_map);
 450         else
 451           breadth_first_search(gr, get(fr, local_start[0]),
 452                                Qr,
 453                                detail::scc_discovery_visitor<const ReverseGraph, Rev_ParentMap>(pred_map),
 454                                color_map);
 455
 456 #ifdef PBGL_SCC_DEBUG
 457   end = accounting::get_time();
 458   if(id == 0)
 459     std::cerr << "  Perform forward and reverse BFS time = " << accounting::print_time(end - start) << " seconds.\n";
 460   start = accounting::get_time();
 461 #endif
 462
 463         // Send predecessors and successors discovered by this proc to the proc responsible for
 464         // this BFS tree
 465         typedef struct detail::v_sets<vertex_descriptor> Vsets;
 466         std::map<vertex_descriptor, Vsets> set_map;
 467
 468         std::map<vertex_descriptor, int> dest_map;
 469
 470         std::vector<VertexPairVec> successors(num_procs);
 471         std::vector<VertexPairVec> predecessors(num_procs);
 472
 473         // Calculate destinations for messages
 474         for (std::size_t i = 0; i < vertex_sets.size(); ++i)
 475           dest_map[vertex_sets[i][0]] = i % num_procs;
 476
 477         for( boost::tie(vstart, vend) = vertices(g); vstart != vend; vstart++ ) {
 478           vertex_descriptor v = get(succ_map, *vstart);
 479           if( v != graph_traits<Graph>::null_vertex() ) {
 480             if (dest_map[v] == id)
 481               set_map[v].succ.push_back(*vstart);
 482             else
 483               successors[dest_map[v]].push_back( std::make_pair(v, *vstart) );
 484           }
 485         }
 486
 487         for( boost::tie(rev_vstart, rev_vend) = vertices(gr); rev_vstart != rev_vend; rev_vstart++ ) {
 488           vertex_descriptor v = get(pred_map, *rev_vstart);
 489           if( v != graph_traits<Graph>::null_vertex() ) {
 490             if (dest_map[v] == id)
 491               set_map[v].pred.push_back(get(rf, *rev_vstart));
 492             else
 493               predecessors[dest_map[v]].push_back( std::make_pair(v, get(rf, *rev_vstart)) );
 494           }
 495         }
 496
 497         // Send predecessor and successor messages
 498         for (process_id_type i = 0; i < num_procs; ++i) {
 499           if (!successors[i].empty()) {
 500             send(pg, i, fhp_succ_size_msg, successors[i].size());
 501             send(pg, i, fhp_succ_msg, &successors[i][0], successors[i].size());
 502           }
 503           if (!predecessors[i].empty()) {
 504             send(pg, i, fhp_pred_size_msg, predecessors[i].size());
 505             send(pg, i, fhp_pred_msg, &predecessors[i][0], predecessors[i].size());
 506           }
 507         }
 508         synchronize(pg);
 509
 510         // Receive predecessor and successor messages and handle them
 511         while (optional<std::pair<process_id_type, int> > m = probe(pg)) {
 512           BOOST_ASSERT(m->second == fhp_succ_size_msg || m->second == fhp_pred_size_msg);
 513           std::size_t num_requests;
 514           receive(pg, m->first, m->second, num_requests);
 515           VertexPairVec requests(num_requests);
 516           if (m->second == fhp_succ_size_msg) {
 517             receive(pg, m->first, fhp_succ_msg, &requests[0],
 518                     num_requests);
 519
 520             std::map<vertex_descriptor, int> added;
 521             for (std::size_t i = 0; i < requests.size(); ++i) {
 522               set_map[requests[i].first].succ.push_back(requests[i].second);
 523               added[requests[i].first]++;
 524             }
 525
 526             // If order of vertex traversal in vertices() is std::less<vertex_descriptor>,
 527             // then the successor sets will be in order
 528             for (std::size_t i = 0; i < local_start.size(); ++i)
 529               if (added[local_start[i]] > 0)
 530                   std::inplace_merge(set_map[local_start[i]].succ.begin(),
 531                                      set_map[local_start[i]].succ.end() - added[local_start[i]],
 532                                      set_map[local_start[i]].succ.end(),
 533                                      std::less<vertex_descriptor>());
 534
 535           } else {
 536             receive(pg, m->first, fhp_pred_msg, &requests[0],
 537                     num_requests);
 538
 539             std::map<vertex_descriptor, int> added;
 540             for (std::size_t i = 0; i < requests.size(); ++i) {
 541               set_map[requests[i].first].pred.push_back(requests[i].second);
 542               added[requests[i].first]++;
 543             }
 544
 545             if (boost::is_same<detail::vertex_identity_property_map<vertex_descriptor>, IsoMapRF>::value)
 546               for (std::size_t i = 0; i < local_start.size(); ++i)
 547                 if (added[local_start[i]] > 0)
 548                   std::inplace_merge(set_map[local_start[i]].pred.begin(),
 549                                      set_map[local_start[i]].pred.end() - added[local_start[i]],
 550                                      set_map[local_start[i]].pred.end(),
 551                                      std::less<vertex_descriptor>());
 552           }
 553         }
 554
 555 #ifdef PBGL_SCC_DEBUG
 556   end = accounting::get_time();
 557   if(id == 0)
 558     std::cerr << "  All gather successors and predecessors time = " << accounting::print_time(end - start) << " seconds.\n";
 559   start = accounting::get_time();
 560 #endif
 561
 562         //
 563         // Filter predecessor and successor sets and perform set arithmetic
 564         //
 565         new_vertex_sets.clear();
 566
 567         if( std::size_t(id) < vertex_sets.size() ) { //If this proc has one or more unique starting points
 568
 569           for( std::size_t i = 0; i < local_start.size(); ++i ) {
 570
 571             vertex_descriptor v = local_start[i];
 572
 573             // Replace this sort with an in-place merges during receive step if possible
 574             if (!boost::is_same<detail::vertex_identity_property_map<vertex_descriptor>, IsoMapRF>::value)
 575               std::sort(set_map[v].pred.begin(), set_map[v].pred.end(), std::less<vertex_descriptor>());
 576
 577             // Limit predecessor and successor sets to members of the original set
 578             std::vector<vertex_descriptor> temp;
 579
 580             std::set_intersection( vertex_sets[id + i*num_procs].begin(), vertex_sets[id + i*num_procs].end(),
 581                                    set_map[v].pred.begin(), set_map[v].pred.end(),
 582                                    back_inserter(temp),
 583                                    std::less<vertex_descriptor>());
 584             set_map[v].pred.clear();
 585             std::swap(set_map[v].pred, temp);
 586
 587             std::set_intersection( vertex_sets[id + i*num_procs].begin(), vertex_sets[id + i*num_procs].end(),
 588                                    set_map[v].succ.begin(), set_map[v].succ.end(),
 589                                    back_inserter(temp),
 590                                    std::less<vertex_descriptor>());
 591             set_map[v].succ.clear();
 592             std::swap(set_map[v].succ, temp);
 593
 594             // Intersection(pred, succ)
 595             std::set_intersection(set_map[v].pred.begin(), set_map[v].pred.end(),
 596                                     set_map[v].succ.begin(), set_map[v].succ.end(),
 597                                     back_inserter(set_map[v].intersect),
 598                                     std::less<vertex_descriptor>());
 599
 600             // Union(pred, succ)
 601             std::set_union(set_map[v].pred.begin(), set_map[v].pred.end(),
 602                            set_map[v].succ.begin(), set_map[v].succ.end(),
 603                            back_inserter(set_map[v].ps_union),
 604                            std::less<vertex_descriptor>());
 605
 606             new_vertex_sets.push_back(std::vector<vertex_descriptor>());
 607             // Original set - Union(pred, succ)
 608             std::set_difference(vertex_sets[id + i*num_procs].begin(), vertex_sets[id + i*num_procs].end(),
 609                                 set_map[v].ps_union.begin(), set_map[v].ps_union.end(),
 610                                 back_inserter(new_vertex_sets[new_vertex_sets.size() - 1]),
 611                                 std::less<vertex_descriptor>());
 612
 613             set_map[v].ps_union.clear();
 614
 615             new_vertex_sets.push_back(std::vector<vertex_descriptor>());
 616             // Pred - Intersect(pred, succ)
 617             std::set_difference(set_map[v].pred.begin(), set_map[v].pred.end(),
 618                                 set_map[v].intersect.begin(), set_map[v].intersect.end(),
 619                                 back_inserter(new_vertex_sets[new_vertex_sets.size() - 1]),
 620                                 std::less<vertex_descriptor>());
 621
 622             set_map[v].pred.clear();
 623
 624             new_vertex_sets.push_back(std::vector<vertex_descriptor>());
 625             // Succ - Intersect(pred, succ)
 626             std::set_difference(set_map[v].succ.begin(), set_map[v].succ.end(),
 627                                 set_map[v].intersect.begin(), set_map[v].intersect.end(),
 628                                 back_inserter(new_vertex_sets[new_vertex_sets.size() - 1]),
 629                                 std::less<vertex_descriptor>());
 630
 631             set_map[v].succ.clear();
 632
 633             // Label SCC just identified with the 'first' vertex in that SCC
 634             for( std::size_t j = 0; j < set_map[v].intersect.size(); j++ )
 635               put(c, set_map[v].intersect[j], set_map[v].intersect[0]);
 636
 637             set_map[v].intersect.clear();
 638           }
 639         }
 640
 641 #ifdef PBGL_SCC_DEBUG
 642   end = accounting::get_time();
 643   if(id == 0)
 644     std::cerr << "  Perform set arithemetic time = " << accounting::print_time(end - start) << " seconds.\n";
 645   start = accounting::get_time();
 646 #endif
 647
 648         // Remove sets of size 1 from new_vertex_sets
 649         typename std::vector<std::vector<vertex_descriptor> >::iterator vviter;
 650         for( vviter = new_vertex_sets.begin(); vviter != new_vertex_sets.end(); /*in loop*/ )
 651           if( (*vviter).size() < 2 )
 652             vviter = new_vertex_sets.erase( vviter );
 653           else
 654             vviter++;
 655
 656         // All gather new sets and recur (gotta marshal and unmarshal sets first)
 657         vertex_sets.clear();
 658         std::vector<vertex_descriptor> serial_sets, all_serial_sets;
 659         detail::marshal_set<Graph>( new_vertex_sets, serial_sets );
 660         all_gather( pg, serial_sets.begin(), serial_sets.end(), all_serial_sets );
 661         detail::unmarshal_set<Graph>( all_serial_sets, vertex_sets );
 662
 663 #ifdef PBGL_SCC_DEBUG
 664   end = accounting::get_time();
 665   if(id == 0) {
 666     std::cerr << "  Serialize and gather new vertex sets time = " << accounting::print_time(end - start) << " seconds.\n\n\n";
 667
 668     printf("Vertex sets: %d\n", (int)vertex_sets.size() );
 669     for( std::size_t i = 0; i < vertex_sets.size(); ++i )
 670       printf("  %d: %d\n", i, (int)vertex_sets[i].size() );
 671   }
 672   start = accounting::get_time();
 673 #endif
 674
 675         // HACK!?!  --  This would be more properly implemented as a topological sort
 676         // Remove vertices without an edge to another vertex in the set and an edge from another
 677         // vertex in the set
 678        typedef typename graph_traits<Graph>::out_edge_iterator out_edge_iterator;
 679        out_edge_iterator estart, eend;
 680        typedef typename graph_traits<ReverseGraph>::out_edge_iterator r_out_edge_iterator;
 681        r_out_edge_iterator restart, reend;
 682        for (std::size_t i = 0; i < vertex_sets.size(); ++i) {
 683          std::vector<vertex_descriptor> new_set;
 684          for (std::size_t j = 0; j < vertex_sets[i].size(); ++j) {
 685            vertex_descriptor v = vertex_sets[i][j];
 686            if (get(owner, v) == id) {
 687              boost::tie(estart, eend) = out_edges(v, g);
 688              while (estart != eend && find(vertex_sets[i].begin(), vertex_sets[i].end(),
 689                                            target(*estart,g)) == vertex_sets[i].end()) estart++;
 690              if (estart != eend) {
 691                boost::tie(restart, reend) = out_edges(get(fr, v), gr);
 692                while (restart != reend && find(vertex_sets[i].begin(), vertex_sets[i].end(),
 693                                                get(rf, target(*restart,gr))) == vertex_sets[i].end()) restart++;
 694                if (restart != reend)
 695                  new_set.push_back(v);
 696              }
 697            }
 698          }
 699          vertex_sets[i].clear();
 700          all_gather(pg, new_set.begin(), new_set.end(), vertex_sets[i]);
 701          std::sort(vertex_sets[i].begin(), vertex_sets[i].end(), std::less<vertex_descriptor>());
 702        }
 703 #ifdef PBGL_SCC_DEBUG
 704   end = accounting::get_time();
 705   if(id == 0)
 706     std::cerr << "  Trim vertex sets time = " << accounting::print_time(end - start) << " seconds.\n\n\n";
 707   start = accounting::get_time();
 708 #endif
 709
 710       } while ( !vertex_sets.empty() );
 711
 712
 713       // Label vertices not in a SCC as their own SCC
 714       for( boost::tie(vstart, vend) = vertices(g); vstart != vend; vstart++ )
 715         if( get(c, *vstart) == graph_traits<Graph>::null_vertex() )
 716           put(c, *vstart, *vstart);
 717
 718       synchronize(c);
 719     }
 720
 721     template<typename Graph, typename ReverseGraph, typename IsoMap>
 722     void
 723     build_reverse_graph( const Graph& g, ReverseGraph& gr, IsoMap& fr, IsoMap& rf )
 724     {
 725       typedef typename graph_traits<Graph>::vertex_iterator vertex_iterator;
 726       typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;
 727       typedef typename graph_traits<Graph>::out_edge_iterator out_edge_iterator;
 728       typedef typename boost::graph::parallel::process_group_type<Graph>::type process_group_type;
 729       typedef typename process_group_type::process_id_type process_id_type;
 730       typedef std::vector<std::pair<vertex_descriptor, vertex_descriptor> > VertexPairVec;
 731
 732       typename property_map<Graph, vertex_owner_t>::const_type
 733         owner = get(vertex_owner, g);
 734
 735       process_group_type pg = process_group(g);
 736       process_id_type id = process_id(pg);
 737
 738       int n;
 739       vertex_iterator vstart, vend;
 740       int num_procs = num_processes(pg);
 741
 742       vertex_descriptor v;
 743       out_edge_iterator oestart, oeend;
 744       for( boost::tie(vstart, vend) = vertices(g); vstart != vend; vstart++ )
 745         {
 746           v = add_vertex(gr);
 747           put(fr, *vstart, v);
 748           put(rf, v, *vstart);
 749         }
 750
 751       gr.distribution() = g.distribution();
 752
 753       int my_n = num_vertices(g);
 754       all_reduce(pg, &my_n, &my_n+1, &n, std::plus<int>());
 755
 756       for (int i = 0; i < n; ++i)
 757         get(fr, vertex(i,g));
 758       synchronize(fr);
 759
 760       // Add edges to gr
 761       std::vector<std::pair<vertex_descriptor, vertex_descriptor> > new_edges;
 762       for( boost::tie(vstart, vend) = vertices(g); vstart != vend; vstart++ )
 763         for( boost::tie(oestart, oeend) = out_edges(*vstart, g); oestart != oeend; oestart++ )
 764           new_edges.push_back( std::make_pair(get(fr, target(*oestart,g)), get(fr, source(*oestart, g))) );
 765
 766       std::vector<VertexPairVec> edge_requests(num_procs);
 767
 768       typename std::vector<std::pair<vertex_descriptor, vertex_descriptor> >::iterator iter;
 769       for( iter = new_edges.begin(); iter != new_edges.end(); iter++ ) {
 770         std::pair<vertex_descriptor, vertex_descriptor> p1 = *iter;
 771         if( get(owner,  p1.first ) == id )
 772           add_edge( p1.first, p1.second, gr );
 773         else
 774           edge_requests[get(owner, p1.first)].push_back(p1);
 775       }
 776       new_edges.clear();
 777
 778       // Send edge addition requests
 779       for (process_id_type p = 0; p < num_procs; ++p) {
 780         if (!edge_requests[p].empty()) {
 781           VertexPairVec reqs(edge_requests[p].begin(), edge_requests[p].end());
 782           send(pg, p, fhp_edges_size_msg, reqs.size());
 783           send(pg, p, fhp_add_edges_msg, &reqs[0], reqs.size());
 784         }
 785       }
 786       synchronize(pg);
 787
 788       // Receive edge addition requests and handle them
 789       while (optional<std::pair<process_id_type, int> > m = probe(pg)) {
 790         BOOST_ASSERT(m->second == fhp_edges_size_msg);
 791         std::size_t num_requests;
 792         receive(pg, m->first, m->second, num_requests);
 793         VertexPairVec requests(num_requests);
 794         receive(pg, m->first, fhp_add_edges_msg, &requests[0],
 795                 num_requests);
 796         for( std::size_t i = 0; i < requests.size(); ++i )
 797           add_edge( requests[i].first, requests[i].second, gr );
 798       }
 799           synchronize(gr);
 800     }
 801
 802     template<typename Graph, typename VertexComponentMap, typename ComponentMap>
 803     typename property_traits<ComponentMap>::value_type
 804     number_components(const Graph& g, VertexComponentMap r, ComponentMap c)
 805     {
 806       typedef typename boost::graph::parallel::process_group_type<Graph>::type process_group_type;
 807       typedef typename graph_traits<Graph>::vertex_iterator vertex_iterator;
 808       typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;
 809       typedef typename property_traits<ComponentMap>::value_type ComponentMapType;
 810       std::vector<vertex_descriptor> my_roots, all_roots;
 811       vertex_iterator vstart, vend;
 812
 813       for( boost::tie(vstart, vend) = vertices(g); vstart != vend; vstart++ )
 814         if( find( my_roots.begin(), my_roots.end(), get(r, *vstart) ) == my_roots.end() )
 815           my_roots.push_back( get(r, *vstart) );
 816
 817       all_gather( process_group(g), my_roots.begin(), my_roots.end(), all_roots );
 818
 819       /* Number components */
 820       std::map<vertex_descriptor, ComponentMapType> comp_numbers;
 821       ComponentMapType c_num = 0;
 822
 823       // Compute component numbers
 824       for (std::size_t i = 0; i < all_roots.size(); ++i )
 825         if ( comp_numbers.count(all_roots[i]) == 0 )
 826           comp_numbers[all_roots[i]] = c_num++;
 827
 828       // Broadcast component numbers
 829       for( boost::tie(vstart, vend) = vertices(g); vstart != vend; vstart++ )
 830         put( c, *vstart, comp_numbers[get(r,*vstart)] );
 831
 832       // Broadcast number of components
 833       if (process_id(process_group(g)) == 0) {
 834         typedef typename process_group_type::process_size_type
 835           process_size_type;
 836         for (process_size_type dest = 1, n = num_processes(process_group(g));
 837               dest != n; ++dest)
 838           send(process_group(g), dest, 0, c_num);
 839       }
 840
 841       synchronize(process_group(g));
 842
 843       if (process_id(process_group(g)) != 0) receive(process_group(g), 0, 0, c_num);
 844
 845       synchronize(c);
 846       return c_num;
 847     }
 848
 849
 850     template<typename Graph, typename ComponentMap, typename VertexComponentMap,
 851              typename VertexIndexMap>
 852     typename property_traits<ComponentMap>::value_type
 853     fleischer_hendrickson_pinar_strong_components_impl
 854       (const Graph& g,
 855        ComponentMap c,
 856        VertexComponentMap r,
 857        VertexIndexMap vertex_index_map,
 858        incidence_graph_tag)
 859     {
 860       typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;
 861       typedef iterator_property_map<typename std::vector<vertex_descriptor>::iterator,
 862                                     VertexIndexMap> IsoMap;
 863       typename boost::graph::parallel::process_group_type<Graph>::type pg = process_group(g);
 864
 865 #ifdef PBGL_SCC_DEBUG
 866   accounting::time_type start = accounting::get_time();
 867 #endif
 868
 869       typedef adjacency_list<listS,
 870                              distributedS<typename boost::graph::parallel::process_group_type<Graph>::type, vecS>,
 871                              directedS > ReverseGraph;
 872
 873       ReverseGraph gr(0, pg);
 874       std::vector<vertex_descriptor> fr_s(num_vertices(g));
 875       std::vector<vertex_descriptor> rf_s(num_vertices(g));
 876       IsoMap fr(fr_s.begin(), vertex_index_map);  // fr = forward->reverse
 877       IsoMap rf(rf_s.begin(), vertex_index_map); // rf = reverse->forward
 878
 879       build_reverse_graph(g, gr, fr, rf);
 880
 881 #ifdef PBGL_SCC_DEBUG
 882   accounting::time_type end = accounting::get_time();
 883   if(process_id(process_group(g)) == 0)
 884     std::cerr << "Reverse graph initialization time = " << accounting::print_time(end - start) << " seconds.\n";
 885 #endif
 886
 887   fleischer_hendrickson_pinar_strong_components(g, r, gr, fr, rf,
 888                                                 vertex_index_map);
 889
 890       typename property_traits<ComponentMap>::value_type c_num = number_components(g, r, c);
 891
 892       return c_num;
 893     }
 894
 895     template<typename Graph, typename ComponentMap, typename VertexComponentMap,
 896              typename VertexIndexMap>
 897     typename property_traits<ComponentMap>::value_type
 898     fleischer_hendrickson_pinar_strong_components_impl
 899       (const Graph& g,
 900        ComponentMap c,
 901        VertexComponentMap r,
 902        VertexIndexMap vertex_index_map,
 903        bidirectional_graph_tag)
 904     {
 905       typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;
 906
 907       reverse_graph<Graph> gr(g);
 908       detail::vertex_identity_property_map<vertex_descriptor> fr, rf;
 909
 910       fleischer_hendrickson_pinar_strong_components(g, r, gr, fr, rf,
 911                                                     vertex_index_map);
 912
 913       typename property_traits<ComponentMap>::value_type c_num
 914         = number_components(g, r, c);
 915
 916       return c_num;
 917     }
 918
 919     template<typename Graph, typename ComponentMap, typename VertexIndexMap>
 920     inline typename property_traits<ComponentMap>::value_type
 921     fleischer_hendrickson_pinar_strong_components
 922       (const Graph& g,
 923        ComponentMap c,
 924        VertexIndexMap vertex_index_map)
 925     {
 926       typedef typename graph_traits<Graph>::vertex_descriptor
 927         vertex_descriptor;
 928       typedef iterator_property_map<typename std::vector<vertex_descriptor>::iterator,
 929                                     VertexIndexMap> VertexComponentMap;
 930       typename boost::graph::parallel::process_group_type<Graph>::type pg
 931         = process_group(g);
 932
 933       if (num_processes(pg) == 1) {
 934         local_subgraph<const Graph> ls(g);
 935         return boost::strong_components(ls, c);
 936       }
 937
 938       // Create a VertexComponentMap for intermediate labeling of SCCs
 939       std::vector<vertex_descriptor> r_s(num_vertices(g), graph_traits<Graph>::null_vertex());
 940       VertexComponentMap r(r_s.begin(), vertex_index_map);
 941
 942       return fleischer_hendrickson_pinar_strong_components_impl
 943                (g, c, r, vertex_index_map,
 944                 typename graph_traits<Graph>::traversal_category());
 945     }
 946
 947     template<typename Graph, typename ComponentMap>
 948     inline typename property_traits<ComponentMap>::value_type
 949     fleischer_hendrickson_pinar_strong_components(const Graph& g,
 950                                                   ComponentMap c)
 951     {
 952       return fleischer_hendrickson_pinar_strong_components(g, c, get(vertex_index, g));
 953     }
 954
 955 }  // end namespace distributed
 956
 957 using distributed::fleischer_hendrickson_pinar_strong_components;
 958
 959 } // end namespace graph
 960
 961 template<class Graph, class ComponentMap, class P, class T, class R>
 962 inline typename property_traits<ComponentMap>::value_type
 963 strong_components
 964  (const Graph& g, ComponentMap comp,
 965   const bgl_named_params<P, T, R>&
 966   BOOST_GRAPH_ENABLE_IF_MODELS_PARM(Graph, distributed_graph_tag))
 967 {
 968   return graph::fleischer_hendrickson_pinar_strong_components(g, comp);
 969 }
 970
 971 template<class Graph, class ComponentMap>
 972 inline typename property_traits<ComponentMap>::value_type
 973 strong_components
 974  (const Graph& g, ComponentMap comp
 975   BOOST_GRAPH_ENABLE_IF_MODELS_PARM(Graph, distributed_graph_tag))
 976 {
 977   return graph::fleischer_hendrickson_pinar_strong_components(g, comp);
 978 }
 979
 980 } /* end namespace boost */
 981
 982 #endif // BOOST_GRAPH_DISTRIBUTED_SCC_HPP