diff --git a/doc/transitive_closure.html b/doc/transitive_closure.html
new file mode 100644
index 00000000..d28af4d2
--- /dev/null
+++ b/doc/transitive_closure.html
@@ -0,0 +1,66 @@
+
+
+
+Boost Graph Library: Transitive Closure
+
+
+
+
+
+
+
+
+
+template <class Graph>
+void transitive_closure(Graph& g)
+
+
+The transitive closure of a graph G = (V,E) is a graph G' =
+(V,E') such that E' contains an edge (u,v) if and
+only if G contains a path from u to
+v. The transitive_closure() function transforms the
+input graph g into the transitive closure of itself.
+
+
+
+Where Defined
+
+
+boost/graph/transitive_closure.hpp
+
+
Parameters
+
+IN/OUT: Graph& g
+
+ A directed graph, where the Graph type must model the
+ Vertex List Graph
+ and Edge Mutable Graph concepts.
+
+
+
+
+
+
+
+
+
+
diff --git a/include/boost/graph/transitive_closure.hpp b/include/boost/graph/transitive_closure.hpp
new file mode 100644
index 00000000..656307e6
--- /dev/null
+++ b/include/boost/graph/transitive_closure.hpp
@@ -0,0 +1,429 @@
+// Copyright (C) 2001 Vladimir Prus
+// Permission to copy, use, modify, sell and distribute this software is
+// granted, provided this copyright notice appears in all copies and
+// modified version are clearly marked as such. This software is provided
+// "as is" without express or implied warranty, and with no claim as to its
+// suitability for any purpose.
+
+#include
+#include
+#include
+#include
+
+namespace boost { namespace graph {
+
+ namespace detail {
+ // These classes really don't belong here!
+ // They should be moved somewhere
+ template
+ struct const_subscript_t : std::unary_function
+ {
+ const_subscript_t(const Container& c) : container(&c) {};
+ const RT& operator()(const ST& i) const
+ { return (*container)[i]; }
+
+ protected:
+ const Container* container;
+ };
+
+ template
+ struct subscript_t : std::unary_function
+ {
+ subscript_t(Container& c) : container(&c) {};
+ RT& operator()(const ST& i) const
+ { return (*container)[i]; }
+
+ protected:
+ Container *container;
+ };
+
+ struct subscript_vb_t : std::unary_function
+ {
+ subscript_vb_t(std::vector& c) : container(&c) {};
+ std::vector::reference operator()(std::size_t i) const
+ { return (*container)[i]; }
+
+ protected:
+ std::vector *container;
+ };
+
+ template
+ const_subscript_t
+ subscript(const Container& c)
+ { return const_subscript_t(c); }
+
+
+ template
+ subscript_t
+ subscript(Container &c)
+ { return subscript_t(c); }
+
+ inline
+ subscript_vb_t
+ subscript(std::vector& v)
+ {
+ return subscript_vb_t(v);
+ }
+
+ template
+ struct truth : std::unary_function
+ {
+ bool operator()(const T&) { return true; }
+ };
+
+ }
+
+ namespace detail {
+ // The following functions are used in implementation of the algorithm
+ // They are not generic in any way, since it would require much work
+ // for no benefit.
+
+ /** Makes adjacencent vertices of any given vertex to appear in
+ topological order.
+ */
+ void topologically_sort_edges(vector< vector >& g,
+ const vector top_num)
+ {
+ for (size_t i = 0; i < g.size(); ++i) {
+
+ sort(g[i].begin(), g[i].end(),
+ compose_f_gx_hy(less(),
+ subscript(top_num),
+ subscript(top_num)));
+ }
+ }
+
+ /** Computed decomposition of vertices into disjoint chains,
+ where chains are pathes in graph such that vertices occur
+ in them in topological order. The averange number of chains
+ for this algorithm in G(n, p) graphs is O(ln(np)/p).
+ */
+ void compute_chains(const vector< vector >& g,
+ const vector& top_order,
+ vector< vector >& chains)
+ {
+ vector in_chains(g.size());
+
+ for (size_t i = 0; i < top_order.size(); ++i) {
+
+ int v = top_order[i];
+ if (!in_chains[v]) {
+
+ chains.resize(chains.size()+1);
+ vector& chain = chains.back();
+
+ for(;;) {
+
+ chain.push_back(v);
+ in_chains[v] = true;
+
+ vector::const_iterator next;
+ next = find_if(g[v].begin(), g[v].end(),
+ not1(subscript(in_chains)));
+ if (next != g[v].end())
+ v = *next;
+ else
+ break;
+ }
+ }
+ }
+ }
+
+ void acyclic_closure(vector< vector >& g)
+ {
+ // The simplest algorithm would be
+ // for u \in reverse_topological_order
+ // for (v in adj(u))
+ // if (v \notin succ(u))
+ // succ(u) = succ(u) | { v } \cup succ(v)
+ // where succ(v) is the set of all transitive succesors of v
+ // if succ union takes linear time, this is clearly
+ // 0(n*e) algorithm.
+ // Different representation for succ impvored averange
+ // performance, but not worst case behaviour.
+ // In this case, chain decomposition is used -- division
+ // of V into disjoint pathes C, where vertices in each
+ // path are topologically ordered. So, succ(u) will
+ // be a vector, telling for each chain index of element
+ // of chain with the smallest topological number.
+
+ using namespace boost;
+ using namespace std;
+
+ vector top_order;
+ vector top_num(g.size());
+ {
+ identity_property_map id;
+ topological_sort(g, back_inserter(top_order),
+ vertex_index_map(id));
+ reverse(top_order.begin(), top_order.end());
+
+ for (size_t i = 0; i < top_order.size(); ++i)
+ top_num[top_order[i]] = i;
+ }
+
+ topologically_sort_edges(g, top_num);
+
+ vector< vector > chains;
+ vector chain_no(g.size()), index_in_chain(g.size());
+ {
+ compute_chains(g, top_order, chains);
+
+ for (size_t i = 0; i < chains.size(); ++i)
+ for (size_t j = 0; j < chains[i].size(); ++j) {
+ int v = chains[i][j];
+ chain_no[v] = i;
+ index_in_chain[v] = j;
+ }
+ }
+
+#ifdef PRINTS
+ cout << "Chains found are :\n" << multiline << chains << endl;
+#endif
+
+ int inf = numeric_limits::max();
+ vector< vector > successors;
+ {
+ successors.resize(g.size(), vector(chains.size(), inf));
+
+ for (int i = top_order.size()-1; i >= 0; --i) {
+ int v = top_order[i];
+
+
+ for (size_t j = 0; j < g[v].size(); ++j) {
+ int av = g[v][j];
+
+ if (top_num[av] < successors[v][chain_no[av]]) {
+
+ for (size_t k = 0; k < chains.size(); ++k)
+ successors[v][k] = min(successors[v][k],
+ successors[av][k]);
+ successors[v][chain_no[av]] = top_num[av];
+ }
+ }
+ }
+ }
+
+#ifdef PRINTS
+ cout << "Successors sets\n" << multiline << successors << endl;
+#endif
+
+
+ for (size_t i = 0; i < g.size(); ++i)
+ g[i].clear();
+
+ for (size_t i = 0; i < g.size(); ++i)
+ for (size_t j = 0; j < chains.size(); ++j) {
+ int s = successors[i][j];
+
+ if (s < inf) {
+ int v = top_order[s];
+ for (size_t k = index_in_chain[v]; k < chains[j].size(); ++k)
+ g[i].push_back(chains[j][k]);
+ }
+ }
+ }
+ }
+
+
+ /** Computes transitive closure of graph g. The algorithm used has
+ worst-case complexity 0(ne) and averange complexity on G(n, p)
+ graphs is O(n ln(np)/p). Note: algorithm temporary removed all
+ the edges, so any information assosiated with edges is lost.
+ */
+ template
+ void transitive_closure(G& g)
+ {
+ using namespace detail;
+
+ function_requires< AdjacencyGraphConcept >();
+ function_requires< MutableGraphConcept >();
+
+ typedef typename graph_traits::vertex_descriptor vertex;
+ typedef typename graph_traits::edge_descriptor edge;
+ typedef typename graph_traits::vertex_iterator vertex_iterator;
+ typedef typename graph_traits::adjacency_iterator adjacency_iterator;
+
+ // Find SCCs
+ vector component_no(g.size());
+ vector< vector > components;
+ {
+ // AAA!
+ identity_property_map id;
+ int n = strong_components(g, &component_no[0], vertex_index_map(id));
+ components.resize(n);
+ for (size_t i = 0; i < g.size(); ++i)
+ components[component_no[i]].push_back(i);
+ }
+
+#ifdef PRINTS
+ cout << "Components are\n" << multiline << components << endl;
+#endif
+
+ // Construct a condensation graph
+ // G=(V,E) -> G'=(V',E')
+ // V' -- set of strong components in G
+ // E' = { (s_1, s_2) : \exists (u \in s_1, v \in s_2) : (u,v) \in E }
+ // Note that G' is acyclic
+ vector< vector > cg(components.size());
+ {
+ for (size_t i = 0; i < components.size(); ++i) {
+
+ vector targets;
+ {
+ for (size_t j = 0; j < components[i].size(); ++j) {
+ int v = components[i][j];
+ for (size_t k = 0; k < g[v].size(); ++k) {
+ int t = component_no[g[v][k]];
+ if (t != i) // Avoid loops in the condensation graph
+ targets.push_back(t);
+ }
+ }
+ }
+ sort(targets.begin(), targets.end());
+ vector::iterator di = unique(targets.begin(), targets.end());
+ if (di != targets.end())
+ targets.erase(di, targets.end());
+
+ cg[i] = targets;
+ }
+ }
+
+#ifdef PRINTS
+ cout << "Condensation graph is\n" << multiline << cg << endl;
+#endif
+
+ detail::acyclic_closure(cg);
+
+#ifdef PRINTS
+ cout << "Its transitive closure is\n" << multiline << cg << endl;
+#endif
+
+ // Create needed edges in original graph.
+ // The complexity of the following code is due to fact that we don't
+ // want *reflexive* transitive closure to be computed, but rather
+ // ordinary one.
+
+ vector looped_vertices;
+ {
+ vertex_iterator vb, ve;
+ for (boost::tie(vb, ve) = vertices(g); vb != ve; ++vb) {
+ adjacency_iterator ab, ae;
+ for (boost::tie(ab, ae) = adjacent_vertices(*vb, g);
+ ab != ae;
+ ++ab)
+ if (*vb == *ab)
+ looped_vertices.push_back(*vb);
+ }
+ }
+
+ remove_edge_if(truth(), g);
+
+ for (size_t i = 0; i < cg.size(); ++i)
+ for (size_t j = 0; j < cg[i].size(); ++j) {
+ int s = i;
+ int t = cg[i][j];
+ for (size_t k = 0; k < components[s].size(); ++k)
+ for (size_t l = 0; l < components[t].size(); ++l)
+ add_edge(components[s][k], components[t][l], g);
+ }
+
+ // Since cg's closure is not reflexive closure, we need to process SCC's
+ for (size_t i = 0; i < components.size(); ++i)
+ if (components[i].size() > 1)
+ for (size_t k = 0; k < components[i].size(); ++k)
+ for (size_t l = 0; l < components[i].size(); ++l)
+ add_edge(components[i][k], components[i][l], g);
+
+ for (size_t i = 0; i < looped_vertices.size(); ++i) {
+ vertex v = looped_vertices[i];
+ if (components[component_no[v]].size() == 1)
+ add_edge(v, v, g);
+ }
+ }
+
+ template
+ void warshall_transitive_closure(G& g)
+ {
+ using namespace boost;
+ typedef typename graph_traits::vertex_descriptor vertex;
+ typedef typename graph_traits::vertex_iterator vertex_iterator;
+
+ function_requires< AdjacencyMatrixConcept >();
+ function_requires< MutableGraphConcept >();
+
+ // Matrix form:
+ // for k
+ // for i
+ // if A[i,k]
+ // for j
+ // A[i,j] = A[i,j] | A[k,j]
+ vertex_iterator ki, ke, ii, ie, ji, je;
+ for (boost::tie(ki, ke) = vertices(g); ki != ke; ++ki)
+ for (boost::tie(ii, ie) = vertices(g); ii != ie; ++ii)
+ if (edge(*ii, *ki, g).second)
+ for (boost::tie(ji, je) = vertices(g); ji != je; ++ji)
+ if (!edge(*ii, *ji, g).second &&
+ edge(*ki, *ji, g).second)
+ {
+ add_edge(*ii, *ji, g);
+ }
+ }
+
+ template
+ void warren_transitive_closure(G& g)
+ {
+ using namespace boost;
+ typedef typename graph_traits::vertex_descriptor vertex;
+ typedef typename graph_traits::vertex_iterator vertex_iterator;
+
+ function_requires< AdjacencyMatrixConcept >();
+ function_requires< MutableGraphConcept >();
+
+ // Make sure second loop will work
+ if (num_vertices(g) == 0)
+ return;
+
+ // for i = 2 to n
+ // for k = 1 to i - 1
+ // if A[i,k]
+ // for j = 1 to n
+ // A[i,j] = A[i,j] | A[k,j]
+
+ vertex_iterator ic, ie, jc, je, kc, ke;
+ for (boost::tie(ic, ie) = vertices(g), ++ic; ic != ie; ++ic)
+ for (boost::tie(kc, ke) = vertices(g); *kc != *ic; ++kc)
+ if (edge(*ic, *kc, g).second)
+ for (boost::tie(jc, je) = vertices(g); jc != je; ++jc)
+ if (!edge(*ic, *jc, g).second &&
+ edge(*kc, *jc, g).second)
+ {
+ add_edge(*ic, *jc, g);
+ }
+
+ // for i = 1 to n - 1
+ // for k = i + 1 to n
+ // if A[i,k]
+ // for j = 1 to n
+ // A[i,j] = A[i,j] | A[k,j]
+
+ for (boost::tie(ic, ie) = vertices(g), --ie; ic != ie; ++ic)
+ for (kc = ic, ke = ie, ++kc; kc != ke; ++kc)
+ if (edge(*ic, *kc, g).second)
+ for (boost::tie(jc, je) = vertices(g); jc != je; ++jc)
+ if (!edge(*ic, *jc, g).second &&
+ edge(*kc, *jc, g).second)
+ {
+ add_edge(*ic, *jc, g);
+ }
+ }
+
+}
+}
+
+
+
+
+
diff --git a/test/transitive_closure_test.cpp b/test/transitive_closure_test.cpp
new file mode 100644
index 00000000..0b1f4afa
--- /dev/null
+++ b/test/transitive_closure_test.cpp
@@ -0,0 +1,167 @@
+// Copyright (C) 2001 Vladimir Prus
+// Permission to copy, use, modify, sell and distribute this software is
+// granted, provided this copyright notice appears in all copies and
+// modified version are clearly marked as such. This software is provided
+// "as is" without express or implied warranty, and with no claim as to its
+// suitability for any purpose.
+
+
+#include "vector_as_graph_plug.hpp"
+#include "transitive_closure.hpp"
+
+#include
+
+#include
+#include
+
+#include
+#include
+#include
+using namespace std;
+
+using namespace boost;
+using namespace boost::graph;
+
+
+void generate_graph(int n, double p, vector< vector >& r1,
+ vector >& r2)
+{
+ static class {
+ public:
+ double operator()() {
+ return double(rand())/RAND_MAX;
+ }
+ } gen;
+
+ r1.clear();
+ r2.clear();
+
+ r1.resize(n);
+ r2.resize(n, vector(n));
+
+ for (int i = 0; i < n; ++i)
+ for (int j = 0; j < n; ++j)
+ if (gen() < p) {
+ r1[i].push_back(j);
+ r2[i][j] = true;
+ }
+}
+
+bool compare_graphs(const vector< vector >& g1,
+ const vector< vector >& g2)
+{
+ int ne1(0);
+ {
+ for (size_t i = 0; i < g1.size(); ++i)
+ ne1 += g1[i].size();
+ }
+ int ne2(0);
+ {
+ for (size_t i = 0; i < g2.size(); ++i)
+ for (size_t j = 0; j < g2.size(); ++j)
+ if (g2[i][j])
+ ++ne2;
+ }
+ if (ne1 == ne2) {
+
+ for (size_t i = 0; i < g1.size(); ++i)
+ for (size_t j = 0; j < g1[i].size(); ++j)
+ if (!g2[i][g1[i][j]])
+ return false;
+
+ } else {
+ return false;
+ }
+ return true;
+}
+
+void warshall_transitive_closure(vector< vector >& g)
+{
+ for (size_t k = 0; k < g.size(); ++k)
+ for (size_t i = 0; i < g.size(); ++i)
+ for (size_t j = 0; j < g.size(); ++j)
+ g[i][j] = g[i][j] || (g[i][k] && g[k][j]);
+}
+
+
+
+void matrix_strong_components(const vector< vector >& g, vector& scc_no)
+{
+ // SCCs are equivalence classes for equivivalence relation defined as
+ // u ~ v <==> (v is reachable from u) & (u is reachable from v)
+
+ vector< vector > reachable(g);
+ warshall_transitive_closure(reachable);
+
+ vector< vector > ec;
+ for (size_t i = 0; i < g.size(); ++i) {
+ size_t j;
+ for (j = 0; j < ec.size(); ++j) {
+ int representative = ec[j][0];
+ if (reachable[i][representative] && reachable[representative][i])
+ ec[j].push_back(i);
+ }
+ if (j == ec.size()) {
+ ec.resize(ec.size()+1);
+ ec.back().push_back(i);
+ }
+ }
+
+ for (size_t i = 0; i < ec.size(); ++i)
+ for (size_t j = 0; j < ec[i].size(); ++j)
+ scc_no[ec[i][j]] = i;
+}
+
+
+
+bool test(int n, double p)
+{
+ vector< vector > g1;
+ vector< vector > g2;
+
+
+ generate_graph(n, p, g1, g2);
+ cout << "Created graph with " << n << " vertices.\n";
+
+ vector< vector > g1_c(g1);
+
+ {
+ progress_timer t;
+ cout << "transitive_closure" << endl;
+ transitive_closure(g1);
+ cout << "warshall_transitive_closure" << endl;
+ warshall_transitive_closure(g2);
+ }
+
+ if(compare_graphs(g1, g2))
+ return true;
+ else {
+ //cout << "Original graph was " << multiline << g1_c << endl;
+ //cout << "Result is " << multiline << g1 << endl;
+ return false;
+ }
+}
+
+
+int main()
+{
+ srand(time(0));
+ static class {
+ public:
+ double operator()() {
+ return double(rand())/RAND_MAX;
+ }
+ } gen;
+
+
+ for (size_t i = 0; i < 100; ++i) {
+ int n = 0 + int(20*gen());
+ double p = gen();
+ if (!test(n, p)) {
+ cout << "Failed." << endl;
+ return 1;
+ }
+ }
+ cout << "Passed." << endl;
+}
+