Understanding performance: graph connected components

Question

I need to analyze directed graphs with 10M edges, 1M vertices, and 300K strongly connected components, so that the largest one contains 400K vertices.

I read some explanations of Leonid Shifrin here, here and here. Although I don't quite understand why his code merge is so blazingly fast, I learned two things: no recursions only iterations, linked lists are cool. Trying to mimic his approach, I prepared two function for finding strongly connected components. The first one is Kosaraju's algorithm:

 ClearAll[readEdge, SowDFS, DFS, kosaraju]

 readEdge[{a_, b_}] := (adjOut[a] = q[b, adjOut[a]]; adjIn[b] = q[a, adjIn[b]]);

 (* Sow as a postvisit vertex function *)
 SowDFS[q[v_?newQ, rest_]] := SowDFS[newQ[v] = False; adjIn[v] /. q[] -> q[Sow[v], rest]];
 SowDFS[q[v_, rest_]] := SowDFS[rest];

 (* deep first scan *)
 DFS[q[v_?newQ, rest_]] := DFS[Sow[v, tag]; newQ[v] = False; adjOut[v] /. q[] -> rest];
 DFS[q[v_, rest_]] := DFS[rest];

 (* finding connected components *)
 kosaraju[graphList_] := Block[{adjOut, adjIn, order, newQ, q, tag = 1, $IterationLimit = Infinity},
 (* construction of adjacency lists *)
 adjOut[v_] = q[]; adjIn[v_] = q[]; Scan[readEdge, graphList];
 SetAttributes[q, HoldAllComplete];
 (* the first scan: topological sort of the reverse graph *)
 newQ[v_] = True;
 order = Reverse[Reap[Scan[SowDFS[q[#, q[]]] &, DeleteDuplicates@Flatten[graphList]]][[2, 1]]];
 (* the second scan: finding components *)
 Clear[adjIn, newQ]; newQ[v_] = True;
 Last@Reap[Scan[(tag++; DFS[q[#, q[]]]) &, order], _, #2 &]
 ]

The second one is based on Tarjan's algorithm, as it is described in Wikipedia. I only adapted it in non-recursive way:

 (* reading edge *)
 readOut[{a_, b_}] := (adj[a] = q[b, adj[a]]);

 (* test to start component *)
 strongQ[v_] := index[v] == link[v];

 (* pop from stack *)
 pop[v_, q[v_, rest_]] := (inQ[Sow[v, v]] = False; rest);
 pop[v_, q[a_, rest_]] := pop[v, inQ[Sow[a, v]] = False; rest];

 (* deep first scan with two stacks *)
 biDFS[q[v_?newQ, rest_], stack_] := biDFS[
    (* previsit function *)
    newQ[v] = False; index[v] = link[v] = ++idx; {p[v], root} = {root, v};
    adj[v] /. q[] -> q[h[v], rest], inQ[v] = True; q[v, stack]];

 biDFS[q[h[v_?strongQ], rest_], stack_] := biDFS[rest, pop[v, stack]];
 biDFS[q[v_?inQ, rest_], stack_] := biDFS[link[root] = Min[link[root], index[v]]; rest, stack];
 biDFS[q[v_, rest_], stack_] := biDFS[rest, stack]

 (* postvisit vertex function *)
 h[a_] := (root = p[a]; link[root] = Min[link[root], link[a]]; a);

 (* start scan *)
 start[v_?newQ] := Block[{p, root = v, inQ, link, index, idx = 0},
 (* p is for parent *)p[a_] := a; inQ[a_] = False;
 biDFS[q[v, q[]], q[]]]

(* finding connected components *)
tarjan[graphList_] := Block[{adj, newQ, q, $IterationLimit = Infinity},
(* construction of adjacency out-lists *)
adj[v_] = q[]; Scan[readOut, graphList]; SetAttributes[q, HoldAllComplete];
(* finding components *)
newQ[v_] = True;
Last@Reap[Scan[start, DeleteDuplicates@Flatten[graphList]], _, #2 &]
]

I try not to use WM build in functions for graph analysis in order to have a clear comparison.

To test the performance I chose Google+ graphs from Stanford Large Network Dataset Collection. The renamed version of one of them is available here.

SetDirectory[NotebookDirectory[]];
graphList = Import["edges.dat"];

This graph contains 5126172 edges but only 16405 vertices. The structure of SCC is not completely trivial: a single connected component with 11064 vertices, another one — 11 vertices, one more — 5 vertices, three ones — 3 vertices each, ten components with 2 vertices, and 5296 vertices which are not strongly connected.

On the one hand, WM somehow requires a lot of time to form the graph:

result2 = ConnectedComponents[
Graph[DirectedEdge @@@ graphList]]; // AbsoluteTiming
(* {139.523577, Null} *)

and only half of a second to find the strongly connected components.

On the other hand, reading the graph takes 29 seconds in my function tarjan, and 12 seconds to find the components, thus 42 second in total:

result1 = tarjan[graphList]; // AbsoluteTiming
(* {42.207118, Null} *)

Hence, my question is it possible to speed up the tarjan function so that it will be at least 10 times slower than the built-in function ConnectedComponents?

My nb-file is available here.