Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 09:37:19

0001 // Copyright (C) 2004-2008 The Trustees of Indiana University.
0002 
0003 // Use, modification and distribution is subject to the Boost Software
0004 // License, Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0005 // http://www.boost.org/LICENSE_1_0.txt)
0006 
0007 //  Authors: Nick Edmonds
0008 //           Douglas Gregor
0009 //           Andrew Lumsdaine
0010 #ifndef BOOST_GRAPH_DISTRIBUTED_SCC_HPP
0011 #define BOOST_GRAPH_DISTRIBUTED_SCC_HPP
0012 
0013 #ifndef BOOST_GRAPH_USE_MPI
0014 #error "Parallel BGL files should not be included unless <boost/graph/use_mpi.hpp> has been included"
0015 #endif
0016 
0017 // #define PBGL_SCC_DEBUG
0018 
0019 #include <boost/assert.hpp>
0020 #include <boost/property_map/property_map.hpp>
0021 #include <boost/property_map/parallel/parallel_property_maps.hpp>
0022 #include <boost/property_map/parallel/distributed_property_map.hpp>
0023 #include <boost/property_map/parallel/caching_property_map.hpp>
0024 #include <boost/graph/parallel/algorithm.hpp>
0025 #include <boost/graph/parallel/process_group.hpp>
0026 #include <boost/graph/distributed/queue.hpp>
0027 #include <boost/graph/distributed/filtered_graph.hpp>
0028 #include <boost/pending/indirect_cmp.hpp>
0029 #include <boost/graph/breadth_first_search.hpp>
0030 #include <boost/graph/graph_traits.hpp>
0031 #include <boost/graph/overloading.hpp>
0032 #include <boost/graph/distributed/concepts.hpp>
0033 #include <boost/graph/distributed/local_subgraph.hpp>
0034 #include <boost/graph/parallel/properties.hpp>
0035 #include <boost/graph/named_function_params.hpp>
0036 #include <boost/graph/random.hpp>
0037 #include <boost/graph/distributed/reverse_graph.hpp>
0038 #include <boost/optional.hpp>
0039 #include <boost/graph/distributed/detail/filtered_queue.hpp>
0040 #include <boost/graph/distributed/adjacency_list.hpp>
0041 #ifdef PBGL_SCC_DEBUG
0042   #include <iostream>
0043   #include <cstdlib>
0044   #include <iomanip>
0045   #include <sys/time.h>
0046   #include <boost/graph/distributed/graphviz.hpp> // for ostringstream
0047 #endif
0048 #include <vector>
0049 #include <map>
0050 #include <boost/graph/parallel/container_traits.hpp>
0051 
0052 #ifdef PBGL_SCC_DEBUG
0053 #  include <boost/graph/accounting.hpp>
0054 #endif /* PBGL_SCC_DEBUG */
0055 
0056 // If your graph is likely to have large numbers of small strongly connected
0057 // components then running the sequential SCC algorithm on the local subgraph
0058 // and filtering components with no remote edges may increase performance
0059 // #define FILTER_LOCAL_COMPONENTS
0060 
0061 namespace boost { namespace graph { namespace distributed { namespace detail {
0062 
0063 template<typename vertex_descriptor>
0064 struct v_sets{
0065   std::vector<vertex_descriptor> pred, succ, intersect, ps_union;
0066 };
0067 
0068 /* Serialize vertex set */
0069 template<typename Graph>
0070 void
0071 marshal_set( std::vector<std::vector<typename graph_traits<Graph>::vertex_descriptor> > in,
0072              std::vector<typename graph_traits<Graph>::vertex_descriptor>& out )
0073 {
0074   for( std::size_t i = 0; i < in.size(); ++i ) {
0075     out.insert( out.end(), graph_traits<Graph>::null_vertex() );
0076     out.insert( out.end(), in[i].begin(), in[i].end() );
0077   }
0078 }
0079 
0080 /* Un-serialize vertex set */
0081 template<typename Graph>
0082 void
0083 unmarshal_set( std::vector<typename graph_traits<Graph>::vertex_descriptor> in,
0084                std::vector<std::vector<typename graph_traits<Graph>::vertex_descriptor> >& out )
0085 {
0086   typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;
0087 
0088   while( !in.empty() ) {
0089     typename std::vector<vertex_descriptor>::iterator end 
0090       = std::find( in.begin(), in.end(), graph_traits<Graph>::null_vertex() );
0091 
0092     if( end == in.begin() )
0093       in.erase( in.begin() );
0094     else {
0095       out.push_back(std::vector<vertex_descriptor>());
0096       out[out.size() - 1].insert( out[out.size() - 1].end(), in.begin(), end );
0097       in.erase( in.begin(), end );
0098     }
0099   }
0100 }
0101 
0102 /* Determine if vertex is in subset */
0103 template <typename Set>
0104 struct in_subset {
0105   in_subset() : m_s(0) { }
0106   in_subset(const Set& s) : m_s(&s) { }
0107 
0108   template <typename Elt>
0109   bool operator()(const Elt& x) const {
0110     return ((*m_s).find(x) != (*m_s).end());
0111   }
0112 
0113 private:
0114   const Set* m_s;
0115 };
0116 
0117 template<typename T>
0118 struct vertex_identity_property_map
0119   : public boost::put_get_helper<T, vertex_identity_property_map<T> >
0120 {
0121   typedef T key_type;
0122   typedef T value_type;
0123   typedef T reference;
0124   typedef boost::readable_property_map_tag category;
0125 
0126   inline value_type operator[](const key_type& v) const { return v; }
0127   inline void clear() { }
0128 };
0129 
0130 template <typename T>
0131 inline void synchronize( vertex_identity_property_map<T> & ) { }
0132 
0133 /* BFS visitor for SCC */
0134 template<typename Graph, typename SourceMap>
0135 struct scc_discovery_visitor : bfs_visitor<>
0136 {
0137   scc_discovery_visitor(SourceMap& sourceM)
0138     : sourceM(sourceM) {}
0139 
0140   template<typename Edge>
0141   void tree_edge(Edge e, const Graph& g)
0142   {
0143     put(sourceM, target(e,g), get(sourceM, source(e,g)));
0144   }
0145 
0146  private:
0147   SourceMap& sourceM;
0148 };
0149 
0150 } } } } /* End namespace boost::graph::distributed::detail */
0151 
0152 namespace boost { namespace graph { namespace distributed {
0153     enum fhp_message_tags { fhp_edges_size_msg, fhp_add_edges_msg, fhp_pred_size_msg, 
0154                             fhp_pred_msg, fhp_succ_size_msg, fhp_succ_msg };
0155 
0156     template<typename Graph, typename ReverseGraph,
0157              typename VertexComponentMap, typename IsoMapFR, typename IsoMapRF,
0158              typename VertexIndexMap>
0159     void
0160     fleischer_hendrickson_pinar_strong_components(const Graph& g,
0161                                                   VertexComponentMap c,
0162                                                   const ReverseGraph& gr,
0163                                                   IsoMapFR fr, IsoMapRF rf,
0164                                                   VertexIndexMap vertex_index_map)
0165     {
0166       typedef typename graph_traits<Graph>::vertex_iterator vertex_iterator;
0167       typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;
0168       typedef typename graph_traits<ReverseGraph>::vertex_iterator rev_vertex_iterator;
0169       typedef typename graph_traits<ReverseGraph>::vertex_descriptor rev_vertex_descriptor;
0170       typedef typename boost::graph::parallel::process_group_type<Graph>::type
0171         process_group_type;
0172       typedef typename process_group_type::process_id_type process_id_type;
0173       typedef iterator_property_map<typename std::vector<vertex_descriptor>::iterator,
0174                                     VertexIndexMap> ParentMap;
0175       typedef iterator_property_map<typename std::vector<default_color_type>::iterator,
0176                                     VertexIndexMap> ColorMap;
0177       typedef iterator_property_map<typename std::vector<vertex_descriptor>::iterator,
0178                                     VertexIndexMap> Rev_ParentMap;
0179       typedef std::vector<std::pair<vertex_descriptor, vertex_descriptor> > VertexPairVec;
0180 
0181       typedef typename property_map<Graph, vertex_owner_t>::const_type
0182         OwnerMap;
0183 
0184       OwnerMap owner = get(vertex_owner, g);
0185 
0186       using boost::graph::parallel::process_group;
0187       process_group_type pg = process_group(g);
0188       process_id_type id = process_id(pg);
0189       int num_procs = num_processes(pg);
0190       int n = 0;
0191 
0192       int my_n = num_vertices(g);
0193       all_reduce(pg, &my_n, &my_n+1, &n, std::plus<int>());
0194 
0195       //
0196       // Initialization
0197       //
0198 
0199 #ifdef PBGL_SCC_DEBUG
0200   accounting::time_type start = accounting::get_time();
0201 #endif
0202 
0203       vertex_iterator vstart, vend;
0204       rev_vertex_iterator rev_vstart, rev_vend;
0205       std::vector<std::vector<vertex_descriptor> > vertex_sets, new_vertex_sets;
0206 
0207       vertex_sets.push_back(std::vector<vertex_descriptor>());
0208 
0209       // Remove vertices that do not have at least one in edge and one out edge
0210       new_vertex_sets.push_back(std::vector<vertex_descriptor>());
0211       for( boost::tie(vstart, vend) = vertices(g); vstart != vend; vstart++ )
0212         if( out_degree( get(fr, *vstart), gr) > 0 && out_degree(*vstart, g) > 0 )
0213           new_vertex_sets[0].push_back( *vstart );
0214 
0215       // Perform sequential SCC on local subgraph, filter all components with external
0216       // edges, mark remaining components and remove them from vertex_sets
0217 #ifdef FILTER_LOCAL_COMPONENTS  
0218       // This doesn't actually speed up SCC in connected graphs it seems, but it does work
0219       // and may help in the case where there are lots of small strong components.
0220       {
0221         local_subgraph<const Graph> ls(g);
0222         typedef typename property_map<local_subgraph<const Graph>, vertex_index_t>::type
0223           local_index_map_type;
0224         local_index_map_type local_index = get(vertex_index, ls);
0225 
0226         std::vector<int> ls_components_vec(num_vertices(ls));
0227         typedef iterator_property_map<std::vector<int>::iterator, local_index_map_type>
0228           ls_components_map_type;
0229         ls_components_map_type ls_component(ls_components_vec.begin(), local_index);
0230         int num_comp = boost::strong_components(ls, ls_component);
0231 
0232         // Create map of components
0233         std::map<int, std::vector<vertex_descriptor> > local_comp_map;
0234         typedef typename graph_traits<local_subgraph<const Graph> >::vertex_iterator ls_vertex_iterator;
0235         ls_vertex_iterator vstart, vend;
0236         for( boost::tie(vstart,vend) = vertices(ls); vstart != vend; vstart++ )
0237           local_comp_map[get(ls_component, *vstart)].push_back( *vstart );
0238 
0239         // Filter components that have no non-local edges
0240         typedef typename graph_traits<Graph>::adjacency_iterator adjacency_iterator;
0241         typedef typename graph_traits<ReverseGraph>::adjacency_iterator rev_adjacency_iterator;
0242         adjacency_iterator abegin, aend;
0243         rev_adjacency_iterator rev_abegin, rev_aend;
0244         for( std::size_t i = 0; i < num_comp; ++i ) {
0245           bool local = true;
0246           for( std::size_t j = 0; j < local_comp_map[i].size(); j++ ) {
0247             for( boost::tie(abegin,aend) = adjacent_vertices(local_comp_map[i][j], g);
0248                  abegin != aend; abegin++ )
0249               if( get(owner, *abegin) != id ) {
0250                 local = false;
0251                 break;
0252               }
0253 
0254             if( local )
0255               for( boost::tie(rev_abegin,rev_aend) = adjacent_vertices(get(fr, local_comp_map[i][j]), gr);
0256                    rev_abegin != rev_aend; rev_abegin++ )
0257                 if( get(owner, *rev_abegin) != id ) {
0258                   local = false;
0259                   break;
0260                 }
0261 
0262             if( !local ) break;
0263           }
0264 
0265           if( local ) // Mark and remove from new_vertex_sets
0266             for( std::size_t j = 0; j < local_comp_map[i].size(); j++ ) {
0267               put( c, local_comp_map[i][j], local_comp_map[i][0] );
0268               typename std::vector<vertex_descriptor>::iterator pos =
0269                 std::find(new_vertex_sets[0].begin(), new_vertex_sets[0].end(), local_comp_map[i][j]);
0270               if( pos != new_vertex_sets[0].end() )
0271                 new_vertex_sets[0].erase(pos);
0272             }
0273         }
0274       }
0275 #endif // FILTER_LOCAL_COMPONENTS
0276 
0277       all_gather( pg, new_vertex_sets[0].begin(), new_vertex_sets[0].end(), vertex_sets[0] );
0278       new_vertex_sets.clear();
0279 
0280 #ifdef PBGL_SCC_DEBUG
0281   accounting::time_type end = accounting::get_time();
0282   if(id == 0)
0283     std::cerr << "Trim local SCCs time = " << accounting::print_time(end - start) << " seconds.\n";
0284 #endif
0285 
0286       if( vertex_sets[0].empty() ) return;
0287 
0288       //
0289       // Recursively determine SCCs
0290       //
0291 
0292 #ifdef PBGL_SCC_DEBUG
0293   int iterations = 0;
0294 #endif
0295 
0296       // Only need to be able to map starting vertices for BFS from now on
0297       fr.clear();
0298 
0299       do {
0300 
0301 #ifdef PBGL_SCC_DEBUG
0302   if(id == 0) {
0303     printf("\n\nIteration %d:\n\n", iterations++);
0304 
0305     if( iterations > 1 ) {
0306       end = accounting::get_time();
0307       std::cerr << "Running main loop destructors time = " << accounting::print_time(end - start) << " seconds.\n";
0308     }
0309 
0310     start = accounting::get_time();
0311   }
0312 #endif
0313 
0314         // Get forward->reverse mappings for BFS start vertices
0315        for (std::size_t i = 0; i < vertex_sets.size(); ++i)
0316           get(fr, vertex_sets[i][0]);
0317         synchronize(fr);
0318 
0319         // Determine local vertices to start BFS from
0320         std::vector<vertex_descriptor> local_start;
0321         for( std::size_t i = id; i < vertex_sets.size(); i += num_procs )
0322           local_start.push_back(vertex_sets[i][0]);
0323 
0324         if( local_start.empty() )
0325           local_start.push_back(vertex_sets[0][0]);
0326 
0327 
0328         // Make filtered graphs
0329         typedef std::set<vertex_descriptor> VertexSet;
0330         typedef std::set<rev_vertex_descriptor> Rev_VertexSet;
0331 
0332         VertexSet filter_set_g;
0333         Rev_VertexSet filter_set_gr;
0334         typename VertexSet::iterator fs;
0335 
0336         int active_vertices = 0;
0337         for (std::size_t i = 0; i < vertex_sets.size(); i++)
0338           active_vertices += vertex_sets[i].size();
0339 
0340         // This is a completely random bound
0341         if ( active_vertices < 0.05*n ) {
0342           // TODO: This set insertion is ridiculously inefficient, make it an in-place-merge?
0343           for (std::size_t i = 0; i < vertex_sets.size(); i++)
0344             filter_set_g.insert(vertex_sets[i].begin(), vertex_sets[i].end());
0345 
0346           for (fs = filter_set_g.begin(); fs != filter_set_g.end(); ++fs )
0347             filter_set_gr.insert(get(fr, *fs));
0348         }
0349 
0350         filtered_graph<const Graph, keep_all, detail::in_subset<VertexSet> >
0351           fg(g, keep_all(), detail::in_subset<VertexSet>(filter_set_g));
0352 
0353         filtered_graph<const ReverseGraph, keep_all, detail::in_subset<VertexSet> >
0354           fgr(gr, keep_all(), detail::in_subset<VertexSet>(filter_set_gr));
0355 
0356         // Add additional starting vertices to BFS queue
0357         typedef filtered_queue<queue<vertex_descriptor>, boost::detail::has_not_been_seen<VertexIndexMap> >
0358           local_queue_type;
0359         typedef boost::graph::distributed::distributed_queue<process_group_type, OwnerMap, local_queue_type>
0360           queue_t;
0361 
0362         typedef typename property_map<ReverseGraph, vertex_owner_t>::const_type
0363           RevOwnerMap;
0364 
0365         typedef filtered_queue<queue<rev_vertex_descriptor>, boost::detail::has_not_been_seen<VertexIndexMap> >
0366           rev_local_queue_type;
0367         typedef boost::graph::distributed::distributed_queue<process_group_type, RevOwnerMap, rev_local_queue_type>
0368           rev_queue_t;
0369 
0370         queue_t Q(process_group(g),
0371                   owner,
0372                   make_filtered_queue(queue<vertex_descriptor>(),
0373                                       boost::detail::has_not_been_seen<VertexIndexMap>
0374                                       (num_vertices(g), vertex_index_map)),
0375                   false);
0376 
0377         rev_queue_t Qr(process_group(gr),
0378                        get(vertex_owner, gr),
0379                        make_filtered_queue(queue<rev_vertex_descriptor>(),
0380                                            boost::detail::has_not_been_seen<VertexIndexMap>
0381                                            (num_vertices(gr), vertex_index_map)),
0382                        false);
0383 
0384         for( std::size_t i = 1; i < local_start.size(); ++i ) {
0385           Q.push(local_start[i]);
0386           Qr.push(get(fr, local_start[i]));
0387         }
0388 
0389 #ifdef PBGL_SCC_DEBUG
0390   end = accounting::get_time();
0391   if(id == 0)
0392     std::cerr << "  Initialize BFS time = " << accounting::print_time(end - start) << " seconds.\n";
0393   start = accounting::get_time();
0394 #endif
0395 
0396 #ifdef PBGL_SCC_DEBUG
0397   accounting::time_type start2 = accounting::get_time();
0398 #endif
0399 
0400         // Forward BFS
0401         std::vector<default_color_type> color_map_s(num_vertices(g));
0402         ColorMap color_map(color_map_s.begin(), vertex_index_map);
0403         std::vector<vertex_descriptor> succ_map_s(num_vertices(g), graph_traits<Graph>::null_vertex());
0404         ParentMap succ_map(succ_map_s.begin(), vertex_index_map);
0405 
0406         for( std::size_t i = 0; i < vertex_sets.size(); ++i )
0407           put(succ_map, vertex_sets[i][0], vertex_sets[i][0]);
0408 
0409 #ifdef PBGL_SCC_DEBUG
0410   accounting::time_type end2 = accounting::get_time();
0411   if(id == 0)
0412     std::cerr << "  Initialize forward BFS time = " << accounting::print_time(end2 - start2) << " seconds.\n";
0413 #endif
0414 
0415         if (active_vertices < 0.05*n)
0416           breadth_first_search(fg, local_start[0], Q,
0417                                detail::scc_discovery_visitor<filtered_graph<const Graph, keep_all,
0418                                                                             detail::in_subset<VertexSet> >, ParentMap>
0419                                (succ_map),
0420                                color_map);
0421         else
0422           breadth_first_search(g, local_start[0], Q,
0423                                detail::scc_discovery_visitor<const Graph, ParentMap>(succ_map),
0424                                color_map);
0425 
0426 #ifdef PBGL_SCC_DEBUG
0427   start2 = accounting::get_time();
0428 #endif
0429 
0430         // Reverse BFS
0431         color_map.clear(); // reuse color map since g and gr have same vertex index
0432         std::vector<vertex_descriptor> pred_map_s(num_vertices(gr), graph_traits<Graph>::null_vertex());
0433         Rev_ParentMap pred_map(pred_map_s.begin(), vertex_index_map);
0434 
0435         for( std::size_t i = 0; i < vertex_sets.size(); ++i )
0436           put(pred_map, get(fr, vertex_sets[i][0]), vertex_sets[i][0]);
0437 
0438 #ifdef PBGL_SCC_DEBUG
0439   end2 = accounting::get_time();
0440   if(id == 0)
0441     std::cerr << "  Initialize reverse BFS time = " << accounting::print_time(end2 - start2) << " seconds.\n";
0442 #endif
0443 
0444         if (active_vertices < 0.05*n)
0445           breadth_first_search(fgr, get(fr, local_start[0]),
0446                                Qr,
0447                                detail::scc_discovery_visitor<filtered_graph<const ReverseGraph, keep_all,
0448                                                                             detail::in_subset<Rev_VertexSet> >, Rev_ParentMap>
0449                                (pred_map),
0450                                color_map);
0451         else
0452           breadth_first_search(gr, get(fr, local_start[0]),
0453                                Qr,
0454                                detail::scc_discovery_visitor<const ReverseGraph, Rev_ParentMap>(pred_map),
0455                                color_map);
0456 
0457 #ifdef PBGL_SCC_DEBUG
0458   end = accounting::get_time();
0459   if(id == 0)
0460     std::cerr << "  Perform forward and reverse BFS time = " << accounting::print_time(end - start) << " seconds.\n";
0461   start = accounting::get_time();
0462 #endif
0463 
0464         // Send predecessors and successors discovered by this proc to the proc responsible for
0465         // this BFS tree
0466         typedef struct detail::v_sets<vertex_descriptor> Vsets;
0467         std::map<vertex_descriptor, Vsets> set_map;
0468 
0469         std::map<vertex_descriptor, int> dest_map;
0470 
0471         std::vector<VertexPairVec> successors(num_procs);
0472         std::vector<VertexPairVec> predecessors(num_procs);
0473 
0474         // Calculate destinations for messages
0475         for (std::size_t i = 0; i < vertex_sets.size(); ++i)
0476           dest_map[vertex_sets[i][0]] = i % num_procs;
0477 
0478         for( boost::tie(vstart, vend) = vertices(g); vstart != vend; vstart++ ) {
0479           vertex_descriptor v = get(succ_map, *vstart);
0480           if( v != graph_traits<Graph>::null_vertex() ) {
0481             if (dest_map[v] == id)
0482               set_map[v].succ.push_back(*vstart);
0483             else
0484               successors[dest_map[v]].push_back( std::make_pair(v, *vstart) );
0485           }
0486         }
0487 
0488         for( boost::tie(rev_vstart, rev_vend) = vertices(gr); rev_vstart != rev_vend; rev_vstart++ ) {
0489           vertex_descriptor v = get(pred_map, *rev_vstart);
0490           if( v != graph_traits<Graph>::null_vertex() ) {
0491             if (dest_map[v] == id)
0492               set_map[v].pred.push_back(get(rf, *rev_vstart));
0493             else
0494               predecessors[dest_map[v]].push_back( std::make_pair(v, get(rf, *rev_vstart)) );
0495           }
0496         }
0497 
0498         // Send predecessor and successor messages
0499         for (process_id_type i = 0; i < num_procs; ++i) {
0500           if (!successors[i].empty()) {
0501             send(pg, i, fhp_succ_size_msg, successors[i].size());
0502             send(pg, i, fhp_succ_msg, &successors[i][0], successors[i].size());
0503           }
0504           if (!predecessors[i].empty()) {
0505             send(pg, i, fhp_pred_size_msg, predecessors[i].size());
0506             send(pg, i, fhp_pred_msg, &predecessors[i][0], predecessors[i].size());
0507           }
0508         }
0509         synchronize(pg);
0510 
0511         // Receive predecessor and successor messages and handle them
0512         while (optional<std::pair<process_id_type, int> > m = probe(pg)) {
0513           BOOST_ASSERT(m->second == fhp_succ_size_msg || m->second == fhp_pred_size_msg);
0514           std::size_t num_requests;
0515           receive(pg, m->first, m->second, num_requests);
0516           VertexPairVec requests(num_requests);
0517           if (m->second == fhp_succ_size_msg) {
0518             receive(pg, m->first, fhp_succ_msg, &requests[0],
0519                     num_requests);
0520 
0521             std::map<vertex_descriptor, int> added;
0522             for (std::size_t i = 0; i < requests.size(); ++i) {
0523               set_map[requests[i].first].succ.push_back(requests[i].second);
0524               added[requests[i].first]++;
0525             }
0526 
0527             // If order of vertex traversal in vertices() is std::less<vertex_descriptor>,
0528             // then the successor sets will be in order
0529             for (std::size_t i = 0; i < local_start.size(); ++i)
0530               if (added[local_start[i]] > 0)
0531                   std::inplace_merge(set_map[local_start[i]].succ.begin(),
0532                                      set_map[local_start[i]].succ.end() - added[local_start[i]],
0533                                      set_map[local_start[i]].succ.end(),
0534                                      std::less<vertex_descriptor>());
0535 
0536           } else {
0537             receive(pg, m->first, fhp_pred_msg, &requests[0],
0538                     num_requests);
0539 
0540             std::map<vertex_descriptor, int> added;
0541             for (std::size_t i = 0; i < requests.size(); ++i) {
0542               set_map[requests[i].first].pred.push_back(requests[i].second);
0543               added[requests[i].first]++;
0544             }
0545 
0546             if (boost::is_same<detail::vertex_identity_property_map<vertex_descriptor>, IsoMapRF>::value)
0547               for (std::size_t i = 0; i < local_start.size(); ++i)
0548                 if (added[local_start[i]] > 0)
0549                   std::inplace_merge(set_map[local_start[i]].pred.begin(),
0550                                      set_map[local_start[i]].pred.end() - added[local_start[i]],
0551                                      set_map[local_start[i]].pred.end(),
0552                                      std::less<vertex_descriptor>());
0553           }
0554         }
0555 
0556 #ifdef PBGL_SCC_DEBUG
0557   end = accounting::get_time();
0558   if(id == 0)
0559     std::cerr << "  All gather successors and predecessors time = " << accounting::print_time(end - start) << " seconds.\n";
0560   start = accounting::get_time();
0561 #endif
0562 
0563         //
0564         // Filter predecessor and successor sets and perform set arithmetic
0565         //
0566         new_vertex_sets.clear();
0567 
0568         if( std::size_t(id) < vertex_sets.size() ) { //If this proc has one or more unique starting points
0569 
0570           for( std::size_t i = 0; i < local_start.size(); ++i ) {
0571 
0572             vertex_descriptor v = local_start[i];
0573 
0574             // Replace this sort with an in-place merges during receive step if possible
0575             if (!boost::is_same<detail::vertex_identity_property_map<vertex_descriptor>, IsoMapRF>::value) 
0576               std::sort(set_map[v].pred.begin(), set_map[v].pred.end(), std::less<vertex_descriptor>());
0577 
0578             // Limit predecessor and successor sets to members of the original set
0579             std::vector<vertex_descriptor> temp;
0580 
0581             std::set_intersection( vertex_sets[id + i*num_procs].begin(), vertex_sets[id + i*num_procs].end(),
0582                                    set_map[v].pred.begin(), set_map[v].pred.end(),
0583                                    back_inserter(temp),
0584                                    std::less<vertex_descriptor>());
0585             set_map[v].pred.clear();
0586             std::swap(set_map[v].pred, temp);
0587 
0588             std::set_intersection( vertex_sets[id + i*num_procs].begin(), vertex_sets[id + i*num_procs].end(),
0589                                    set_map[v].succ.begin(), set_map[v].succ.end(),
0590                                    back_inserter(temp),
0591                                    std::less<vertex_descriptor>());
0592             set_map[v].succ.clear();
0593             std::swap(set_map[v].succ, temp);
0594 
0595             // Intersection(pred, succ)
0596             std::set_intersection(set_map[v].pred.begin(), set_map[v].pred.end(),
0597                                     set_map[v].succ.begin(), set_map[v].succ.end(),
0598                                     back_inserter(set_map[v].intersect),
0599                                     std::less<vertex_descriptor>());
0600 
0601             // Union(pred, succ)
0602             std::set_union(set_map[v].pred.begin(), set_map[v].pred.end(),
0603                            set_map[v].succ.begin(), set_map[v].succ.end(),
0604                            back_inserter(set_map[v].ps_union),
0605                            std::less<vertex_descriptor>());
0606 
0607             new_vertex_sets.push_back(std::vector<vertex_descriptor>());
0608             // Original set - Union(pred, succ)
0609             std::set_difference(vertex_sets[id + i*num_procs].begin(), vertex_sets[id + i*num_procs].end(),
0610                                 set_map[v].ps_union.begin(), set_map[v].ps_union.end(),
0611                                 back_inserter(new_vertex_sets[new_vertex_sets.size() - 1]),
0612                                 std::less<vertex_descriptor>());
0613 
0614             set_map[v].ps_union.clear();
0615 
0616             new_vertex_sets.push_back(std::vector<vertex_descriptor>());
0617             // Pred - Intersect(pred, succ)
0618             std::set_difference(set_map[v].pred.begin(), set_map[v].pred.end(),
0619                                 set_map[v].intersect.begin(), set_map[v].intersect.end(),
0620                                 back_inserter(new_vertex_sets[new_vertex_sets.size() - 1]),
0621                                 std::less<vertex_descriptor>());
0622 
0623             set_map[v].pred.clear();
0624 
0625             new_vertex_sets.push_back(std::vector<vertex_descriptor>());
0626             // Succ - Intersect(pred, succ)
0627             std::set_difference(set_map[v].succ.begin(), set_map[v].succ.end(),
0628                                 set_map[v].intersect.begin(), set_map[v].intersect.end(),
0629                                 back_inserter(new_vertex_sets[new_vertex_sets.size() - 1]),
0630                                 std::less<vertex_descriptor>());
0631 
0632             set_map[v].succ.clear();
0633 
0634             // Label SCC just identified with the 'first' vertex in that SCC
0635             for( std::size_t j = 0; j < set_map[v].intersect.size(); j++ )
0636               put(c, set_map[v].intersect[j], set_map[v].intersect[0]);
0637 
0638             set_map[v].intersect.clear();
0639           }
0640         }
0641 
0642 #ifdef PBGL_SCC_DEBUG
0643   end = accounting::get_time();
0644   if(id == 0)
0645     std::cerr << "  Perform set arithemetic time = " << accounting::print_time(end - start) << " seconds.\n";
0646   start = accounting::get_time();
0647 #endif
0648 
0649         // Remove sets of size 1 from new_vertex_sets
0650         typename std::vector<std::vector<vertex_descriptor> >::iterator vviter;
0651         for( vviter = new_vertex_sets.begin(); vviter != new_vertex_sets.end(); /*in loop*/ )
0652           if( (*vviter).size() < 2 )
0653             vviter = new_vertex_sets.erase( vviter );
0654           else
0655             vviter++;
0656 
0657         // All gather new sets and recur (gotta marshal and unmarshal sets first)
0658         vertex_sets.clear();
0659         std::vector<vertex_descriptor> serial_sets, all_serial_sets;
0660         detail::marshal_set<Graph>( new_vertex_sets, serial_sets );
0661         all_gather( pg, serial_sets.begin(), serial_sets.end(), all_serial_sets );
0662         detail::unmarshal_set<Graph>( all_serial_sets, vertex_sets );
0663 
0664 #ifdef PBGL_SCC_DEBUG
0665   end = accounting::get_time();
0666   if(id == 0) {
0667     std::cerr << "  Serialize and gather new vertex sets time = " << accounting::print_time(end - start) << " seconds.\n\n\n";
0668 
0669     printf("Vertex sets: %d\n", (int)vertex_sets.size() );
0670     for( std::size_t i = 0; i < vertex_sets.size(); ++i )
0671       printf("  %d: %d\n", i, (int)vertex_sets[i].size() );
0672   }
0673   start = accounting::get_time();
0674 #endif
0675 
0676         // HACK!?!  --  This would be more properly implemented as a topological sort
0677         // Remove vertices without an edge to another vertex in the set and an edge from another
0678         // vertex in the set
0679        typedef typename graph_traits<Graph>::out_edge_iterator out_edge_iterator;
0680        out_edge_iterator estart, eend;
0681        typedef typename graph_traits<ReverseGraph>::out_edge_iterator r_out_edge_iterator;
0682        r_out_edge_iterator restart, reend;
0683        for (std::size_t i = 0; i < vertex_sets.size(); ++i) {
0684          std::vector<vertex_descriptor> new_set;
0685          for (std::size_t j = 0; j < vertex_sets[i].size(); ++j) {
0686            vertex_descriptor v = vertex_sets[i][j];
0687            if (get(owner, v) == id) {
0688              boost::tie(estart, eend) = out_edges(v, g);
0689              while (estart != eend && find(vertex_sets[i].begin(), vertex_sets[i].end(),
0690                                            target(*estart,g)) == vertex_sets[i].end()) estart++;
0691              if (estart != eend) {
0692                boost::tie(restart, reend) = out_edges(get(fr, v), gr);
0693                while (restart != reend && find(vertex_sets[i].begin(), vertex_sets[i].end(),
0694                                                get(rf, target(*restart,gr))) == vertex_sets[i].end()) restart++;
0695                if (restart != reend)
0696                  new_set.push_back(v);
0697              }
0698            }
0699          }
0700          vertex_sets[i].clear();
0701          all_gather(pg, new_set.begin(), new_set.end(), vertex_sets[i]);
0702          std::sort(vertex_sets[i].begin(), vertex_sets[i].end(), std::less<vertex_descriptor>());
0703        }
0704 #ifdef PBGL_SCC_DEBUG
0705   end = accounting::get_time();
0706   if(id == 0)
0707     std::cerr << "  Trim vertex sets time = " << accounting::print_time(end - start) << " seconds.\n\n\n";
0708   start = accounting::get_time();
0709 #endif
0710 
0711       } while ( !vertex_sets.empty() );
0712 
0713 
0714       // Label vertices not in a SCC as their own SCC
0715       for( boost::tie(vstart, vend) = vertices(g); vstart != vend; vstart++ )
0716         if( get(c, *vstart) == graph_traits<Graph>::null_vertex() )
0717           put(c, *vstart, *vstart);
0718 
0719       synchronize(c);
0720     }
0721 
0722     template<typename Graph, typename ReverseGraph, typename IsoMap>
0723     void
0724     build_reverse_graph( const Graph& g, ReverseGraph& gr, IsoMap& fr, IsoMap& rf )
0725     {
0726       typedef typename graph_traits<Graph>::vertex_iterator vertex_iterator;
0727       typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;
0728       typedef typename graph_traits<Graph>::out_edge_iterator out_edge_iterator;
0729       typedef typename boost::graph::parallel::process_group_type<Graph>::type process_group_type;
0730       typedef typename process_group_type::process_id_type process_id_type;
0731       typedef std::vector<std::pair<vertex_descriptor, vertex_descriptor> > VertexPairVec;
0732 
0733       typename property_map<Graph, vertex_owner_t>::const_type
0734         owner = get(vertex_owner, g);
0735 
0736       process_group_type pg = process_group(g);
0737       process_id_type id = process_id(pg);
0738 
0739       int n;
0740       vertex_iterator vstart, vend;
0741       int num_procs = num_processes(pg);
0742 
0743       vertex_descriptor v;
0744       out_edge_iterator oestart, oeend;
0745       for( boost::tie(vstart, vend) = vertices(g); vstart != vend; vstart++ )
0746         {
0747           v = add_vertex(gr);
0748           put(fr, *vstart, v);
0749           put(rf, v, *vstart);
0750         }
0751 
0752       gr.distribution() = g.distribution();
0753 
0754       int my_n = num_vertices(g);
0755       all_reduce(pg, &my_n, &my_n+1, &n, std::plus<int>());
0756 
0757       for (int i = 0; i < n; ++i)
0758         get(fr, vertex(i,g));
0759       synchronize(fr);
0760 
0761       // Add edges to gr
0762       std::vector<std::pair<vertex_descriptor, vertex_descriptor> > new_edges;
0763       for( boost::tie(vstart, vend) = vertices(g); vstart != vend; vstart++ )
0764         for( boost::tie(oestart, oeend) = out_edges(*vstart, g); oestart != oeend; oestart++ )
0765           new_edges.push_back( std::make_pair(get(fr, target(*oestart,g)), get(fr, source(*oestart, g))) );
0766 
0767       std::vector<VertexPairVec> edge_requests(num_procs);
0768 
0769       typename std::vector<std::pair<vertex_descriptor, vertex_descriptor> >::iterator iter;
0770       for( iter = new_edges.begin(); iter != new_edges.end(); iter++ ) {
0771         std::pair<vertex_descriptor, vertex_descriptor> p1 = *iter;
0772         if( get(owner,  p1.first ) == id )
0773           add_edge( p1.first, p1.second, gr );
0774         else
0775           edge_requests[get(owner, p1.first)].push_back(p1);
0776       }
0777       new_edges.clear();
0778 
0779       // Send edge addition requests
0780       for (process_id_type p = 0; p < num_procs; ++p) {
0781         if (!edge_requests[p].empty()) {
0782           VertexPairVec reqs(edge_requests[p].begin(), edge_requests[p].end());
0783           send(pg, p, fhp_edges_size_msg, reqs.size());
0784           send(pg, p, fhp_add_edges_msg, &reqs[0], reqs.size());
0785         }
0786       }
0787       synchronize(pg);
0788 
0789       // Receive edge addition requests and handle them
0790       while (optional<std::pair<process_id_type, int> > m = probe(pg)) {
0791         BOOST_ASSERT(m->second == fhp_edges_size_msg);
0792         std::size_t num_requests;
0793         receive(pg, m->first, m->second, num_requests);
0794         VertexPairVec requests(num_requests);
0795         receive(pg, m->first, fhp_add_edges_msg, &requests[0],
0796                 num_requests);
0797         for( std::size_t i = 0; i < requests.size(); ++i )
0798           add_edge( requests[i].first, requests[i].second, gr );
0799       }
0800           synchronize(gr);
0801     }
0802 
0803     template<typename Graph, typename VertexComponentMap, typename ComponentMap>
0804     typename property_traits<ComponentMap>::value_type
0805     number_components(const Graph& g, VertexComponentMap r, ComponentMap c)
0806     {
0807       typedef typename boost::graph::parallel::process_group_type<Graph>::type process_group_type;
0808       typedef typename graph_traits<Graph>::vertex_iterator vertex_iterator;
0809       typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;
0810       typedef typename property_traits<ComponentMap>::value_type ComponentMapType;
0811       std::vector<vertex_descriptor> my_roots, all_roots;
0812       vertex_iterator vstart, vend;
0813 
0814       for( boost::tie(vstart, vend) = vertices(g); vstart != vend; vstart++ )
0815         if( find( my_roots.begin(), my_roots.end(), get(r, *vstart) ) == my_roots.end() )
0816           my_roots.push_back( get(r, *vstart) );
0817 
0818       all_gather( process_group(g), my_roots.begin(), my_roots.end(), all_roots );
0819 
0820       /* Number components */
0821       std::map<vertex_descriptor, ComponentMapType> comp_numbers;
0822       ComponentMapType c_num = 0;
0823 
0824       // Compute component numbers
0825       for (std::size_t i = 0; i < all_roots.size(); ++i )
0826         if ( comp_numbers.count(all_roots[i]) == 0 )
0827           comp_numbers[all_roots[i]] = c_num++;
0828 
0829       // Broadcast component numbers
0830       for( boost::tie(vstart, vend) = vertices(g); vstart != vend; vstart++ )
0831         put( c, *vstart, comp_numbers[get(r,*vstart)] );
0832 
0833       // Broadcast number of components
0834       if (process_id(process_group(g)) == 0) {
0835         typedef typename process_group_type::process_size_type
0836           process_size_type;
0837         for (process_size_type dest = 1, n = num_processes(process_group(g));
0838               dest != n; ++dest)
0839           send(process_group(g), dest, 0, c_num);
0840       }
0841 
0842       synchronize(process_group(g));
0843 
0844       if (process_id(process_group(g)) != 0) receive(process_group(g), 0, 0, c_num);
0845 
0846       synchronize(c);
0847       return c_num;
0848     }
0849 
0850 
0851     template<typename Graph, typename ComponentMap, typename VertexComponentMap,
0852              typename VertexIndexMap>
0853     typename property_traits<ComponentMap>::value_type
0854     fleischer_hendrickson_pinar_strong_components_impl
0855       (const Graph& g,
0856        ComponentMap c,
0857        VertexComponentMap r,
0858        VertexIndexMap vertex_index_map,
0859        incidence_graph_tag)
0860     {
0861       typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;
0862       typedef iterator_property_map<typename std::vector<vertex_descriptor>::iterator,
0863                                     VertexIndexMap> IsoMap;
0864       typename boost::graph::parallel::process_group_type<Graph>::type pg = process_group(g);
0865 
0866 #ifdef PBGL_SCC_DEBUG
0867   accounting::time_type start = accounting::get_time();
0868 #endif
0869 
0870       typedef adjacency_list<listS,
0871                              distributedS<typename boost::graph::parallel::process_group_type<Graph>::type, vecS>,
0872                              directedS > ReverseGraph;
0873 
0874       ReverseGraph gr(0, pg);
0875       std::vector<vertex_descriptor> fr_s(num_vertices(g));
0876       std::vector<vertex_descriptor> rf_s(num_vertices(g));
0877       IsoMap fr(fr_s.begin(), vertex_index_map);  // fr = forward->reverse
0878       IsoMap rf(rf_s.begin(), vertex_index_map); // rf = reverse->forward
0879 
0880       build_reverse_graph(g, gr, fr, rf);
0881 
0882 #ifdef PBGL_SCC_DEBUG
0883   accounting::time_type end = accounting::get_time();
0884   if(process_id(process_group(g)) == 0)
0885     std::cerr << "Reverse graph initialization time = " << accounting::print_time(end - start) << " seconds.\n";
0886 #endif
0887 
0888   fleischer_hendrickson_pinar_strong_components(g, r, gr, fr, rf, 
0889                                                 vertex_index_map);
0890 
0891       typename property_traits<ComponentMap>::value_type c_num = number_components(g, r, c);
0892 
0893       return c_num;
0894     }
0895 
0896     template<typename Graph, typename ComponentMap, typename VertexComponentMap,
0897              typename VertexIndexMap>
0898     typename property_traits<ComponentMap>::value_type
0899     fleischer_hendrickson_pinar_strong_components_impl
0900       (const Graph& g,
0901        ComponentMap c,
0902        VertexComponentMap r,
0903        VertexIndexMap vertex_index_map,
0904        bidirectional_graph_tag)
0905     {
0906       typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;
0907 
0908       reverse_graph<Graph> gr(g);
0909       detail::vertex_identity_property_map<vertex_descriptor> fr, rf;
0910 
0911       fleischer_hendrickson_pinar_strong_components(g, r, gr, fr, rf, 
0912                                                     vertex_index_map);
0913 
0914       typename property_traits<ComponentMap>::value_type c_num
0915         = number_components(g, r, c);
0916 
0917       return c_num;
0918     }
0919 
0920     template<typename Graph, typename ComponentMap, typename VertexIndexMap>
0921     inline typename property_traits<ComponentMap>::value_type
0922     fleischer_hendrickson_pinar_strong_components
0923       (const Graph& g,
0924        ComponentMap c,
0925        VertexIndexMap vertex_index_map)
0926     {
0927       typedef typename graph_traits<Graph>::vertex_descriptor
0928         vertex_descriptor;
0929       typedef iterator_property_map<typename std::vector<vertex_descriptor>::iterator,
0930                                     VertexIndexMap> VertexComponentMap;
0931       typename boost::graph::parallel::process_group_type<Graph>::type pg 
0932         = process_group(g);
0933 
0934       if (num_processes(pg) == 1) {
0935         local_subgraph<const Graph> ls(g);
0936         return boost::strong_components(ls, c);
0937       }
0938 
0939       // Create a VertexComponentMap for intermediate labeling of SCCs
0940       std::vector<vertex_descriptor> r_s(num_vertices(g), graph_traits<Graph>::null_vertex());
0941       VertexComponentMap r(r_s.begin(), vertex_index_map);
0942 
0943       return fleischer_hendrickson_pinar_strong_components_impl
0944                (g, c, r, vertex_index_map,
0945                 typename graph_traits<Graph>::traversal_category());
0946     }
0947 
0948     template<typename Graph, typename ComponentMap>
0949     inline typename property_traits<ComponentMap>::value_type
0950     fleischer_hendrickson_pinar_strong_components(const Graph& g,
0951                                                   ComponentMap c)
0952     {
0953       return fleischer_hendrickson_pinar_strong_components(g, c, get(vertex_index, g));
0954     }
0955 
0956 }  // end namespace distributed
0957 
0958 using distributed::fleischer_hendrickson_pinar_strong_components;
0959 
0960 } // end namespace graph
0961 
0962 template<class Graph, class ComponentMap, class P, class T, class R>
0963 inline typename property_traits<ComponentMap>::value_type
0964 strong_components
0965  (const Graph& g, ComponentMap comp,
0966   const bgl_named_params<P, T, R>&
0967   BOOST_GRAPH_ENABLE_IF_MODELS_PARM(Graph, distributed_graph_tag))
0968 {
0969   return graph::fleischer_hendrickson_pinar_strong_components(g, comp);
0970 }
0971 
0972 template<class Graph, class ComponentMap>
0973 inline typename property_traits<ComponentMap>::value_type
0974 strong_components
0975  (const Graph& g, ComponentMap comp
0976   BOOST_GRAPH_ENABLE_IF_MODELS_PARM(Graph, distributed_graph_tag))
0977 {
0978   return graph::fleischer_hendrickson_pinar_strong_components(g, comp);
0979 }
0980 
0981 } /* end namespace boost */
0982 
0983 #endif // BOOST_GRAPH_DISTRIBUTED_SCC_HPP