You are on page 1of 13

An Experimental Study of a Parallel Shortest Path Algorithm

for Solving Large-Scale Graph Instances

Kamesh Madduri∗ David A. Bader∗ Jonathan W.Berry† Joseph R. Crobak‡

Abstract niques can be applied to obtain exact solutions for mem-

We present an experimental study of the single source short- ory and compute-intensive graph problems quickly. For
est path problem with non-negative edge weights (NSSP) on instance, recent experimental studies on Breadth-First
large-scale graphs using the ∆-stepping parallel algorithm. Search for large-scale graphs show that a parallel in-core
We report performance results on the Cray MTA-2, a mul- implementation is two orders of magnitude faster than
tithreaded parallel computer. The MTA-2 is a high-end an optimized external memory implementation [5, 2].
shared memory system offering two unique features that In this paper, we present an efficient parallel implemen-
aid the efficient parallel implementation of irregular algo- tation for the single source shortest paths problem that
rithms: the ability to exploit fine-grained parallelism, and can handle scale-free instances in the order of billions of
low-overhead synchronization primitives. Our implementa- edges. In addition, we conduct an experimental study of
tion exhibits remarkable parallel speedup when compared performance on several other graph families, also used
with competitive sequential algorithms, for low-diameter in the 9th DIMACS Implementation Challenge [17] on
sparse graphs. For instance, ∆-stepping on a directed scale- Shortest Paths. Please refer to our technical report [43]
free graph of 100 million vertices and 1 billion edges takes for additional performance details.
less than ten seconds on 40 processors of the MTA-2, with a Sequential algorithms for the single source shortest
relative speedup of close to 30. To our knowledge, these are path problem with non-negative edge weights (NSSP)
the first performance results of a shortest path problem on are studied extensively, both theoretically [20, 18, 23,
realistic graph instances in the order of billions of vertices 24, 54, 56, 33, 30, 46] and experimentally [19, 28, 27,
and edges. 15, 59, 29]. Let m and n denote the number of edges
and vertices in the graph respectively. Nearly all NSSP
1 Introduction algorithms are based on the classical Dijkstra’s [20] algo-
rithm. Using Fibonacci heaps [23], Dijkstra’s algorithm
We present an experimental study of the ∆-stepping
can be implemented in O(m + n log n) time. Thorup
parallel algorithm [49] for solving the single source
[56] presents an O(m + n) RAM algorithm for undi-
shortest path problem on large-scale graph instances.
rected graphs that differs significantly different from Di-
In addition to applications in combinatorial optimiza-
jkstra’s approach. Instead of visiting vertices in the or-
tion problems, shortest path algorithms are finding in-
der of increasing distance, it traverses a component tree.
creasing relevance in the domain of complex network
Meyer [47] and Goldberg [29] propose simple algorithms
analysis. Popular graph theoretic analysis metrics such
with linear average time for uniformly distributed edge
as betweenness centrality [25, 9, 39, 41, 32] are based
on shortest path algorithms. Our parallel implemen-
Parallel algorithms for solving NSSP are reviewed
tation targets graph families that are representative of
in detail by Meyer and Sanders [46, 49]. There are no
real-world, large-scale networks [7, 22, 12, 51, 50]. Real-
known PRAM algorithms that run in sub-linear time
world graphs are typically characterized by a low di-
and O(m + n log n) work. Parallel priority queues [21,
ameter, heavy-tailed degree distributions modeled by
11] for implementing Dijkstra’s algorithm have been de-
power laws, and self-similarity. They are often very
veloped, but these linear work algorithms have a worst-
large, with the number of vertices and edges ranging
case time bound of Ω(n), as they only perform edge
from several hundreds of thousands to billions. On cur-
relaxations in parallel. Several matrix-multiplication
rent workstations, it is not possible to do exact in-core
based algorithms [34, 26], proposed for the parallel All-
computations on these graphs due to the limited phys-
Pairs Shortest Paths (APSP), involve running time and
ical memory. In such cases, parallel computing tech-
efficiency trade-offs. Parallel approximate NSSP algo-
∗ Georgia rithms [40, 16, 55] based on the randomized Breadth-
Institute of Technology
† Sandia National Laboratories First search algorithm of Ullman and Yannakakis [58]
‡ Rutgers University
run in sub-linear time. However, it is not known how grained parallelism and low-overhead word-level syn-
to use the Ullman-Yannakakis randomized approach for chronization. The MTA-2 has no data cache; rather
exact NSSP computations in sub-linear time. than using a memory hierarchy to reduce latency, the
Meyer and Sanders give the ∆-stepping [49] NSSP MTA-2 processors use hardware multithreading to tol-
algorithm that divides Dijkstra’s algorithm into a num- erate the latency. The word-level synchronization sup-
ber of phases, each of which can be executed in par- port complements multithreading and makes perfor-
allel. For random graphs with uniformly distributed mance primarily a function of parallelism. Since graph
edge weights, this algorithm runs in sub-linear time with algorithms have an abundance of parallelism, yet often
linear average case work. Several theoretical improve- are not amenable to partitioning, the MTA-2 architec-
ments [48, 44, 45] are given for ∆-stepping (for instance, tural features lead to superior performance and scalabil-
finding shortcut edges, adaptive bucket-splitting), but ity. Our recent results highlight the exceptional perfor-
it is unlikely that they would be faster than the simple mance of the MTA-2 for implementations of key com-
∆-stepping algorithm in practice, as the improvements binatorial optimization and graph theoretic problems
involve sophisticated data structures that are hard to such as list ranking [3], connected components [3, 8],
implement efficiently. On a random d-regular graph in- subgraph isomorphism [8], Breadth-First Search and st -
stance (219 vertices and d = 3), Meyer and Sanders connectivity [5].
report a speedup of 9.2 on 16 processors of an Intel The main contributions of this paper are as follows:
Paragon machine, for a distributed memory implemen-
tation of the simple ∆-stepping algorithm. For the same • An experimental study of solving the single-source
graph family, we are able to solve problems three or- shortest paths problem in parallel using the ∆-
ders of magnitude larger with near-linear speedup on stepping algorithm. Prior studies have predomi-
the Cray MTA-2. For instance, we achieve a speedup of nantly focused on running sequential NSSP algo-
14.82 on 16 processors and 29.75 on 40 processors for a rithms on graph families that can be easily parti-
random d-regular graph of size 229 vertices and d set to tioned, whereas we also consider several arbitrary,
3. sparse graph instances. We also analyze perfor-
The literature contains few experimental studies on mance using machine-independent algorithmic op-
parallel NSSP algorithms [35, 52, 37, 57]. Prior im- eration counts.
plementation results on distributed memory machines
• Demonstration of the power of massive multithread-
resorted to graph partitioning [14, 1, 31], and running a
ing for graph algorithms on highly unstructured in-
sequential NSSP algorithm on the sub-graph. Heuristics
stances. We achieve impressive performance on
are used for load balancing and termination detection
low-diameter random and scale-free graphs.
[36, 38]. The implementations perform well for certain
graph families and problem sizes, but in the worst case, • Solving NSSP for large-scale realistic graph in-
there is no speedup. stances in the order of billions of edges. ∆-stepping
Implementations of PRAM graph algorithms for ar- on a synthetic directed scale-free graph of 100 mil-
bitrary sparse graphs are typically memory intensive, lion vertices and 1 billion edges takes 9.73 seconds
and the memory accesses are fine-grained and highly ir- on 40 processors of the MTA-2, with a relative
regular. This often leads to poor performance on cache- speedup of approximately 31. These are the first
based systems. On distributed memory clusters, few results that we are aware of, for solving instances
parallel graph algorithms outperform the best sequen- of this scale and also achieving near-linear speedup.
tial implementations due to long memory latencies and Also, the sequential performance of our implemen-
high synchronization costs [4, 3]. Parallel shared mem- tation is comparable to competitive NSSP imple-
ory systems are a more supportive platform. They offer mentations.
higher memory bandwidth and lower latency than clus-
ters, and the global shared memory greatly improves 2 Review of the ∆-stepping Algorithm
developer productivity. However, parallelism is depen-
Let G = (V, E) be a graph with n vertices and m
dent on the cache performance of the algorithm [53] and
edges. Let s ∈ V denote the source vertex. Each
scalability is limited in most cases.
edge e ∈ E is assigned a non-negative real weight
We present our shortest path implementation re-
by the length function l : E → R. Define the
sults on the Cray MTA-2, a massively multithreaded
weight of a path as the sum of the weights of its edges.
parallel machine. The MTA-2 is a high-end shared
The single source shortest paths problem with non-
memory system offering two unique features that aid
negative edge weights (NSSP) computes δ(v), the weight
considerably in the design of irregular algorithms: fine-
of the shortest (minimum-weighted) path from s to v.
δ(v) = ∞ if v is unreachable from s. We set δ(s) = 0.
Most shortest path algorithms maintain a tentative
distance value for each vertex, which are updated by
edge relaxations. Let d(v) denote the tentative distance Algorithm 1: ∆-stepping algorithm
of a vertex v. d(v) is initially set to ∞, and is an upper Input: G(V, E), source vertex s, length func-
bound on δ(v). Relaxing an edge hv, wi ∈ E sets d(w) tion l : E → R
to the minimum of d(w) and d(v) + l(v, w). Based on Output: δ(v), v ∈ V , the weight of the shortest
the manner in which the tentative distance values are path from s to v
updated, most shortest path algorithms can be classified
into two types: label-setting or label-correcting. Label- 1 foreach v ∈ V do

setting algorithms (for instance, Dijkstra’s algorithm) 2 heavy(v) ←− {hv, wi ∈ E : l(v, w) > ∆};
perform relaxations only from settled (d(v) = δ(v)) 3 light(v) ←− {hv, wi ∈ E : l(v, w) ≤ ∆};
vertices, and compute the shortest path from s to all 4 d(v) ←− ∞;
vertices in exactly m edge relaxations. Based on the 5 relax(s, 0);
values of d(v) and δ(v), at each iteration of a shortest 6 i ←− 0;
path algorithm, vertices can be classified into unreached 7 while B is not empty do
(d(v) = ∞), queued (d(v) is finite, but v is not settled) 8 S ←− φ;
or settled. Label-correcting algorithms (e.g., Bellman- 9 while B[i] 6= φ do
Ford) relax edges from unsettled vertices also, and may 10 Req ←− {(w, d(v) + l(v, w)) : v ∈ B[i] ∧
perform more than m relaxations. Also, all vertices hv, wi ∈ light(v)};
remain in a queued state until the final step of the 11 S ←− S ∪ B[i];
algorithm. ∆-stepping belongs to the label-correcting 12 B[i] ←− φ;
type of shortest path algorithms. 13 foreach (v, x) ∈ Req do
The ∆-stepping algorithm (see Alg. 1) is an “ap- 14 relax(v, x);
proximate bucket implementation of Dijkstra’s algo-
rithm” [49]. It maintains an array of buckets B such 15 Req ←− {(w, d(v) + l(v, w)) : v ∈ S ∧
that B[i] stores the set of vertices {v ∈ V : v is queued hv, wi ∈ heavy(v)};
and d(v) ∈ [i∆, (i + 1)∆)}. ∆ is a positive real number 16 foreach (v, x) ∈ Req do
that denotes the “bucket width”. 17 relax(v, x);
In each phase of the algorithm (the inner while loop 18 i ←− i + 1;
in Alg. 1, lines 9–14, when bucket B[i] is not empty), all 19 foreach v ∈ V do
vertices are removed from the current bucket, added to 20 δ(v) ←− d(v);
the set S, and light edges (l(e) ≤ ∆, e ∈ E) adjacent to
these vertices are relaxed (see Alg. 2). This may result
in new vertices being added to the current bucket, which
are deleted in the next phase. It is also possible that
vertices previously deleted from the current bucket may
be reinserted, if their tentative distance is improved.
Heavy edges (l(e) > ∆, e ∈ E) are not relaxed in
a phase, as they result in tentative values outside the
current bucket. Once the current bucket remains empty Algorithm 2: The relax routine in the ∆-stepping
after relaxations, all heavy edges out of the vertices algorithm
in S are relaxed at once (lines 15–17 in Alg. 1). The Input: v, weight request x
algorithm continues until all the buckets are empty. Output: Assignment of v to appropriate
Observe that edge relaxations in each phase can bucket
be done in parallel, as long as individual tentative
1 if x < d(v) then
distance values are updated atomically. The number
2 B [⌊d(v)/∆⌋] ← B [⌊d(v)/∆⌋] \{v};
of phases bounds the parallel running time, and the
3 B [⌊x/∆⌋] ← B [⌊x/∆⌋] ∪ {v};
number of reinsertions (insertions of vertices previously
4 d(v) ← x;
deleted) and rerelaxations (relaxation of their out-going
edges) costs an overhead over Dijkstra’s algorithm. The
performance of the algorithm also depends on the value
of the bucket-width ∆. For ∆ = ∞, the algorithm is
similar to the Bellman-Ford algorithm. It has a high We implement a timed pre-processing step to semi-
degree of parallelism, but is inefficient compared to sort the edges based on the value of ∆. All the light
Dijkstra’s algorithm. ∆-stepping tries to find a good edges adjacent to a vertex are identified in parallel
compromise between the number of parallel phases and and stored in contiguous virtual locations, and so we
the number of re-insertions. Theoretical bounds on the visit only light edges in a phase. The O(n) work pre-
number of phases and re-insertions, and the average case processing step scales well in parallel on the MTA-2.
analysis of the parallel algorithm are presented in [49]. We also support fast parallel insertions into the
We summarize the salient results. request set R. R stores hv, xi pairs, where v ∈ V and
Let dc denote the maximum shortest path weight, x is the requested tentative distance for v. We add a
and P∆ denote the set of paths with weight at most vertex v to R only if it satisfies the condition x < d(v).
∆. Define a parameter lmax , an upper bound on the We do not store duplicates in R. We use a sparse
maximum number of edges in any path in P∆ . The set representation similar to one used by Briggs and
following results hold true for any graph family. Torczon [10] for storing vertices in R. This sparse data
structure uses two arrays of size n: a dense array that
• The number of buckets in B is ⌈dc /∆⌉. contiguously stores the elements of the set, and a sparse
array that indicates whether the vertex is a member
• The total number of reinsertions is bounded by of the set. Thus, it is easy to iterate over the request
|P∆ |, and the total number of rerelaxations is set, and membership queries and insertions are constant
bounded by |P2∆ |. time. Unlike other Dijkstra-based algorithms, we do not
relax edges in one step. Instead, we inspect adjacencies
• The number of phases is bounded by d∆c lmax , i.e.,
(light edges) in each phase, construct a request set of
no bucket is expanded more than lmax times.
vertices, and then relax vertices in the relax step.
For graph families with random edge weights and All vertices in the request set R are relaxed in
a maximum degree of d, Meyer and Sanders [49] theo- parallel in the relax routine. In this step, we first delete
retically prove that ∆ = θ(1/d) is a good compromise a vertex from the old bucket, and then insert it into the
between work efficiency and parallelism. The sequen- new bucket. Instead of performing individual insertions,
tial algorithm performs O(dn) expected work divided we first determine the expansion factor of each bucket,
between O( d∆c · logloglogn n ) phases with high probability . expand the buckets, and add then all vertices into their
In practice, in case of graph families for which dc is new buckets in one step. Since there are no duplicates
O(log n) or O(1), the parallel implementation of ∆- in the request set, no synchronization is involved for
stepping yields sufficient parallelism for our parallel sys- updating the tentative distance values.
tem. To saturate the MTA-2 processors with work and to
obtain high system utilization, we need to minimize the
3 Parallel Implementation of ∆-stepping number of phases and non-empty buckets, and maximize
the request set sizes. Entering and exiting a parallel
The bucket array B is the primary data structure used
phase involves a negligible running time overhead in
by the parallel ∆-stepping algorithm. We implement
practice. However, if the number of phases is O(n),
individual buckets as dynamic arrays that can be re-
this overhead dominates the actual running time of
sized when needed and iterated over easily. To support
the implementation. Also, we enter the relax routine
constant time insertions and deletions, we maintain two
once every phase. The number of implicit barrier
auxiliary arrays of size n: a mapping of the vertex ID
synchronizations in the algorithm is proportional to the
to its current bucket, and a mapping from the vertex ID
number of phases. Our implementation reduces the
to the position of the vertex in the current bucket (see
number of barriers. Our source code for the ∆-stepping
Fig. 1 for an illustration). All new vertices are added to
implementation, along with the MTA-2 graph generator
the end of the array, and deletions of vertices are done
ports, is freely available online [42].
by setting the corresponding locations in the bucket and
the mapping arrays to −1. Note that once bucket i is
4 Experimental Setup
finally empty after a light edge relaxation phase, there
will be no more insertions into the bucket in subsequent 4.1 Platforms We report parallel performance re-
phases. Thus, the memory can be reused once we are sults on a 40-processor Cray MTA-2 system with 160
done relaxing the light edges in the current bucket. Also GB uniform shared memory. Each processor has a clock
observe that all the insertions are done in the relax rou- speed of 220 MHz and support for 128 hardware threads.
tine, which is called once in each phase, and once for The ∆-stepping code is written in C with MTA-2 spe-
relaxing the heavy edges. cific pragmas and directives for parallelization. We com-
Figure 1: Bucket array and auxiliary data structures

pile it using the MTA-2 C compiler (Cray Programming • Scale-free graphs: We use the R-MAT graph model
Environment (PE) 2.0.3) with -O3 and -par flags. [13] for real-world networks to generate scale-free
The MTA-2 code also compiles and runs on sequen- graphs. We define the family ScaleFree4-n similar
tial processors without any modifications. Our test plat- to random graphs.
form for the sequential performance results is one pro-
cessor of a dual-core 3.2 GHz 64-bit Intel Xeon machine • Log-uniform weight distribution: The above graph
with 6GB memory, 1MB cache and running RedHat En- generators assume randomly distributed edge
terprise Linux 4 (linux kernel 2.6.9). We compare the weights. We report results for an additional log-
sequential performance of our implementation with the uniform distribution also. The generated integer
DIMACS reference solver [17]. Both the codes are com- edge weights are of the form 2i , where i is chosen
piled with the Intel C compiler (icc) Version 9.0, with from the uniform random distribution [1, log C] (C
the flags -O3. denotes the maximum integer edge weight). We
define Random4logUnif-n and ScaleFree4logUnif-n
4.2 Problem Instances We evaluate sequential and families for this weight distribution.
parallel performance on several graph families. Some of
the generators and graph instances are part of the DI- 4.3 Methodology For sequential runs, we report the
MACS Shortest Path Implementation Challenge bench- execution time of the reference DIMACS NSSP solver
mark package [17]: (an efficient implementation of Goldberg’s algorithm
[30], which has expected-case linear time for some in-
• Random graphs: Random graphs are generated by puts) and the baseline Breadth-First Search (BFS) on
first constructing a Hamiltonian cycle, and then every graph family. The BFS running time is a natural
adding m − n edges to the graph at random. The lower bound for NSSP codes and is a good indicator of
generator may produce parallel edges as well as how optimized the shortest path implementations are.
self-loops. We define the random graph family It is reasonable to directly compare the execution times
Random4-n such that n is varied, m n = 4, and the
of the reference code and our implementation: both use
edge weights are chosen from a uniform random a similar adjacency array representation for the graph,
distribution. are written in C, and compiled and run in identical ex-
perimental settings. Note that our implementation is
• Grid graphs: This synthetic generator produces optimized for the MTA-2 and we make no modifications
two-dimensional meshes with grid dimensions x to the code before running on a sequential machine. The
and y. Long-n√(x = 16 n
, y = 16) and Square-n time taken for semi-sorting and mechanisms to reduce
grid (x = y = n) families are defined, similar to memory contention on the MTA-2 both constitute over-
random graphs. head on a sequential processor. Also, our implementa-
tion assumes real-weighted edges, and we cannot use
• Road graphs: Road graph families with transit time fast bitwise operations. By default, we set the value of
(USA-road-t ) and distance (USA-road-d ) as the ∆ to m for all graph instances. We will show that this
length function. choice of ∆ may not be optimal for all graph classes and
weight distributions.
In addition, we also study the following families: On a sequential processor, we execute the BFS and
shortest path codes on all the core graph families, for
the recommended problem sizes. However, for parallel 3.0
runs, we only report results for sufficiently large graph DIMACS Reference
instances in case of the synthetic graph families. We 2.5 Delta-stepping

parallelize the synthetic core graph generators and port

Execution TIme (seconds)

them to run on the MTA-2. 2.0

Our implementations accept both directed and

undirected graphs. For all the synthetic graph in- 1.5

stances, we report execution times on directed graphs in

this paper. The road networks are undirected graphs. 1.0

We also assume the edge weights to be distributed in

[0, 1] in the ∆-stepping implementation. So we have a 0.5

pre-processing step to scale the integer edge weights in

the core problem families to the interval [0, 1], dividing 0.0
Random4-n Long-n Square-n FLA (d) FLA (t)
the integer weights by the maximum edge weight.
Graph Family
On the MTA-2, we compare our implementation
running time with the execution time of a multithreaded
Figure 2: Sequential performance of our ∆-stepping
level-synchronized breadth-first search [6], optimized for
implementation on the core graph families. All the
low-diameter graphs. The multithreaded BFS scales as
synthetic graphs are directed, with 220 vertices and m
n ≈
well as ∆-stepping for all the graph instances consid-
4. FLA(d) and FLA(t) are road networks corresponding
ered, and the execution time serves as a lower bound
to Florida, with 1070376 vertices and 2712768 edges
for the shortest path running time.
The first run on the MTA-2 is usually slower than
subsequent ones (by about 10% for a typical ∆-stepping
machine-independent algorithm operation counts. The
run). So we report the average running time for 10
parallel performance is dependent on the value of ∆,
successive runs. We run the code from three randomly
the number of phases, the size of the request set in each
chosen source vertices and average the running time.
phase. Fig. 3 plots the size of the light request set in
We found that using three sources consistently gave
each phase, for different graph families. By default, ∆
us execution time results with little variation on both
is set to 0.25 for all runs. If the request set size is
the MTA-2 and the reference sequential platform. We
less than 10, it is not plotted. Consider the random
tabulate the sequential and parallel performance metrics
graph family (Fig. 3(a)). It executes in 84 phases, and
in [43], and report execution time in seconds.
the request set sizes vary from 0 to 27,000. Observe
the recurring pattern of three bars stacked together in
5 Results and Analysis
the plot. This indicates that all the light edges in a
5.1 Sequential Performance First we present the bucket are relaxed in roughly three phases, and the
performance results of our implementation on the ref- bucket then becomes empty. The size of the relax set is
erence sequential platform, experimenting with various relatively high for several phases, which provides scope
graph families. Fig. 2 compares the execution time for exploiting multithreaded parallelism. The relax set
across graph instances of the same size, but from dif- sizes of a similar problem instance from the Long grid
ferent families. The DIMACS reference code is about family (Fig. 3(b)) stands in stark contrast to that of
1.5 to 2 times faster than our implementation for large the random graph. It takes about 200,000 phases to
problem instances in each family. The running time execute, and the maximum request size is only 15. Both
on the Random4-n is slightly higher than the rest of of these values indicate that our implementation would
the families. For additional details such as performance fare poorly on long grid graphs (e.g. meshes with a
as we vary the problem size for BFS, ∆-stepping, and very high aspect ratio). On square grids (Fig. 3(c)), ∆-
the DIMACS implementation execution times, please stepping takes fewer phases, and the request set sizes
refer to Section B.1 of [43]. Our key observation is that go up to 500. For a road network instance (NE USA-
the ratio of the ∆-stepping execution time to the BFS road-d, Fig. 3(d)), the algorithm takes 23,000 phases to
time varies between 3 and 10 across different problem execute, and only a few phases (about 30) have request
instances. sets greater than 1000 in size.
Fig. 4 plots several key ∆-stepping operation counts
5.2 ∆-stepping analysis To better understand the for various graph classes. All synthetic graphs are
algorithm performance across graph families, we use roughly of the same size. Fig. 4(a) plots the average
(a) Random4-n family, n = 220 . (b) Long-n family, n = 220 .

(c) Square-n family, n = 220 . (d) USA-road-d family, Northeast USA (NE). n = 1524452, m =

Figure 3: ∆-stepping algorithm: Size of the light request set at the end of each phase, for the core graph families.
Request set sizes less than 10 are not plotted.
shortest path weight for various graph classes. Scale- We define the speedup on p processors of the MTA-
free and Long grid graphs are on the two extremes. 2 as the ratio of the execution time on 1 processor
A log-uniform edge weight distribution also results in to the execution time on p processors. In all graph
low average edge weight. The number of phases (see classes except long grids, there is sufficient parallelism to
Fig. 4(b)) is highest for Long grid graphs. The number saturate a single processor of the MTA-2 for reasonably
of buckets shows a similar trend as the average shortest large problem instances.
path weight. Fig. 4(d) plots the total number of As expected, ∆-stepping performs best for low-
insertions for each graph family. The number of vertices diameter random and scale-free graphs with randomly
is 220 for all graph families (slightly higher for the road distributed edge weights (see Fig. 6(a) and 6(b)). We
network), and so ∆-stepping results in roughly 20% achieve a speedup of approximately 31 on 40 processors
overhead in insertions for all the graph families with for a directed random graph of nearly a billion edges,
random edge weights. Note the number of insertions and the ratio of the BFS and ∆-stepping execution time
for graphs with log-uniform weight distributions. ∆- is a constant factor (about 3-5) throughout. The imple-
stepping performs a significant amount of excess work mentation performs equally well for scale-free graphs,
for these families, because the value of ∆ is quite high that are more difficult for partitioning-based parallel
for this particular distribution. computing models to handle due to the irregular degree
We next evaluate the performance of the algorithm distribution. The execution time on 40 processors of the
as ∆ is varied (tables in Section B.2). Fig. 5 plots MTA-2 for the scale-free graph instance is within 9% (a
the execution time of various graph instances on a difference of less than one second) of the running time
sequential machine, and one processor of the MTA- for a random graph and the speedup is approximately
2. ∆ is varied from 0.1 to 10 in each case. We find 30 on 40 processors. We have already shown that the
that the absolute running times on a 3.2 GHz Xeon execution time for smaller graph instances on a sequen-
processor and the MTA-2 are comparable for random, tial machine is comparable to the DIMACS reference
square grid and road network instances. However, on implementation, a competitive NSSP algorithm. Thus,
long grid graphs (Fig. 5(b)), the MTA-2 execution time achieving a speedup of 30 for a realistic scale-free graph
is two orders of magnitude greater than the sequential instance of one billion edges (Fig. 6(b)) is a substan-
time. The number of phases and the total number of tial result. To our knowledge, these are the first results
relaxations vary as ∆ is varied (See Section B.2 in [43]). to demonstrate near-linear speedup for such large-scale
On the MTA-2, the running time is not only dependent unstructured graph instances.
on the work done, but also on the number of phases and In case of all the graph families, the relative speedup
the average number of relax requests in a phase. For increases as the problem size is increased (for e.g., on 40
instance, in the case of long grids (see Fig. 5(b), with processors, the speedup for a Random4-n instance with
execution time plotted on a log scale), the running time n = 221 is just 3.96, whereas it is 31.04 for 228 vertices).
decreases significantly as the value of ∆ is decreased, as This is because there is insufficient parallelism in a
the number of phases reduce. On a sequential processor, problem instance of size 221 to saturate 40 processors
however, the running time is only dependent on the of the MTA-2. As the problem size increases, the ratio
work done (number of insertions). If the value of ∆ of ∆-stepping execution time to multithreaded BFS
is greater than the average shortest path weight, we running time decreases. On an average, ∆-stepping is 5
perform excess work and the running time noticeably times slower than BFS for this graph family.
increases (observe the execution time for ∆ = 5, 10 on For random graphs with a log-uniform weight dis-
the random graph and the road network). The optimal tribution, ∆ set to m results in a significant amount of
value of ∆ (and the execution time on the MTA-2) is additional work. The ∆-stepping to BFS ratio is typi-
also dependent on the number of processors. For a cally 40 in this case, about 8 times higher than the cor-
particular ∆, it may be possible to saturate a single responding ratio for random graphs with random edge
processor of the MTA-2 with the right balance of work weights. However, the execution time scales well with
and phases. The execution time on a 40-processor run the number of processors for large problem sizes.
may not be minimal with this value of ∆. In case of Long-n graphs and ∆ set to m , there is
insufficient parallelism to fully utilize even a single pro-
5.3 Parallel Performance We present the parallel cessor of the MTA-2. The execution time of the level-
scaling of the ∆-stepping algorithm in detail. We synchronous BFS also does not scale with the number
ran ∆-stepping and the level-synchronous parallel BFS of processors. In fact, the running time goes up in case
on all graph instances described in Section 4.2 (see of multiprocessor runs, as the parallelization overhead
[43] for complete tabular results from all experiments). becomes significant. Note that the execution time on a
(a) Average shortest path weight ( n (b) No. of phases
∗ v∈V δ(v))

(c) Last non-empty bucket (d) Number of relax requests

Figure 4: ∆-stepping algorithm performance statistics for various graph classes. All synthetic graph instances
have n set to 220 and m ≈ 4n. Rand-rnd: Random graph with random edge weights, Rand-logU: Random graphs
with log-uniform edge weights, Scale-free: Scale-free graph with random edge weights, Lgrid: Long grid, SqGrid:
Square grid, USA NE: 1524452 vertices, 3897634 edges. Plots (a), (b), (c) are on a log scale, while (d) uses a
linear scale.
(a) Random4-n family. 220 vertices (b) Long-n family. 220 vertices

(c) Square-n family. 220 vertices (d) USA-road-d family, Florida (FLA). 1070376 vertices, 2712798

Figure 5: A comparison of the execution time on the reference sequential platform and a single MTA-2 processor,
as the bucket-width ∆ is varied.
(a) (b)

Figure 6: ∆-stepping execution time and relative speedup on the MTA-2 for Random4-n (left) and ScaleFree4-n
(right) graph instances (directed graph, n=228 vertices and m = 4n edges, random edge weights).

single processor of the MTA-2 is two orders of magni- algorithm. For high diameter graphs, there is a trade-off
tude slower than the reference sequential processor. In between the number of phases and the amount of work
case of square grid graphs, there is sufficient parallelism done (proportional to the number of bucket insertions).
to utilize up to 4 processors for a graph instance of 224 The execution time is dependent on the value of ∆ as
vertices. For all smaller instances, the running time well as the number of processors. We need to reduce
does not scale for multiprocessor runs. The ratio of the the number of phases for parallel runs and increase
running time to BFS is about 5 in this case, and the the system utilization by choosing an appropriate value
∆-stepping MTA-2 single processor time is comparable of ∆. Our parallel performance studies have been
to the sequential reference platform running time for restricted to the Cray MTA-2 in this paper. We
smaller instances. For the road networks, we note that have designed and have a preliminary implementation
the execution time does not scale well with the number of ∆-stepping for multi-core processors and symmetric
of processors, as the problem instances are quite small. multiprocessors (SMPs), and for future work we will
We observe better performance (lower execution time, analyze its performance.
better speedup) on USA-road-d graphs than on USA-
road-t graphs. Acknowledgments
This work was supported in part by NSF Grants CNS-
6 Conclusions and Future Work 0614915, CAREER CCF-0611589, ACI-00-93039, NSF
In this paper, we present an efficient implementation DBI-0420513, ITR ACI-00-81404, ITR EIA-01-21377,
of the parallel ∆-stepping NSSP algorithm along with Biocomplexity DEB-01-20709, ITR EF/BIO 03-31654,
an experimental evaluation. We study the perfor- and DARPA Contract NBCH30390004. We acknowl-
mance for several graph families on the Cray MTA-2, edge the algorithmic inputs from Bruce Hendrickson of
and observe that our implementation execution time Sandia National Laboratories. We would also like to
scales very well with number of processors for low- thank John Feo of Cray for helping us optimize the
diameter sparse graphs. Few prior implementations MTA-2 implementation. Sandia is a multipurpose lab-
achieve parallel speedup for NSSP, whereas we attain oratory operated by Sandia Corporation, a Lockheed-
near-linear speedup for several large-scale low-diameter Martin Company, for the United States Department of
graph families. We also analyze the performance using Energy under contract DE-AC04-94AL85000.
platform-independent ∆-stepping algorithm operation
counts such as the number of phases and the request set
sizes to explain performance across graph instances.
We intend to further study the dependence of the
bucket-width ∆ on the parallel performance of the
References [14] K.M. Chandy and J. Misra. Distributed computation
on graphs: Shortest path algorithms. Communications
of the ACM, 25(11):833–837, 1982.
[1] P. Adamson and E. Tick. Greedy partitioned algo- [15] B.V. Cherkassky, A.V. Goldberg, and T. Radzik.
rithms for the shortest path problem. Internat. J. Par- Shortest paths algorithms: theory and experimental
allel Program, 20(4):271–278, 1991. evaluation. Mathematical Programming, 73:129–174,
[2] D. Ajwani, R. Dementiev, and U. Meyer. A computa- 1996.
tional study of external-memory BFS algorithms. In [16] E. Cohen. Using selective path-doubling for paral-
Proc. 17th Ann. Symp. Discrete Algorithms (SODA- lel shortest-path computation. J. Algs., 22(1):30–56,
06), pages 601–610, Miami, FL, January 2006. ACM 1997.
Press. [17] C. Demetrescu, A. Goldberg, and D. Johnson. 9th
[3] D.A. Bader, G. Cong, and J. Feo. On the architec- DIMACS implementation challenge – Shortest Paths.
tural requirements for efficient execution of graph al-∼ challenge9/, 2006.
gorithms. In Proc. 34th Int’l Conf. on Parallel Process- [18] R.B. Dial. Algorithm 360: Shortest path forest with
ing (ICPP), Oslo, Norway, June 2005. IEEE Computer topological ordering. Communications of the ACM,
Society. 12:632–633, 1969.
[4] D.A. Bader, A.K. Illendula, B. M.E. Moret, and [19] R.B. Dial, F. Glover, D. Karney, and D. Klingman. A
N. Weisse-Bernstein. Using PRAM algorithms on computational analysis of alternative algorithms and
a uniform-memory-access shared-memory architecture. labeling techniques for finding shortest path trees.
In G.S. Brodal, D. Frigioni, and A. Marchetti- Networks, 9:215–248, 1979.
Spaccamela, editors, Proc. 5th Int’l Workshop on Al- [20] E.W. Dijkstra. A note on two problems in connexion
gorithm Engineering (WAE 2001), volume 2141 of Lec- with graphs. Numerische Mathematik, 1:269–271,
ture Notes in Computer Science, pages 129–144, Århus, 1959.
Denmark, 2001. Springer-Verlag. [21] J.R. Driscoll, H.N. Gabow, R. Shrairman, and R.E.
[5] D.A. Bader and K. Madduri. Designing multithreaded Tarjan. Relaxed heaps: An alternative to fibonacci
algorithms for breadth-first search and st-connectivity heaps with applications to parallel computation. Com-
on the Cray MTA-2. In Proc. 35th Int’l Conf. on munications of the ACM, 31(11):1343–1354, 1988.
Parallel Processing (ICPP), Columbus, OH, August [22] M. Faloutsos, P. Faloutsos, and C. Faloutsos. On
2006. IEEE Computer Society. power-law relationships of the Internet topology. In
[6] D.A. Bader and K. Madduri. Parallel algorithms for Proc. ACM SIGCOMM, pages 251–262, Cambridge,
evaluating centrality indices in real-world networks. In MA, August 1999. ACM.
Proc. 35th Int’l Conf. on Parallel Processing (ICPP), [23] M.L. Fredman and R.E. Tarjan. Fibonacci heaps
Columbus, OH, August 2006. IEEE Computer Society. and their uses in improved network optimization al-
[7] A.-L. Barabási and R. Albert. Emergence of scaling in gorithms. J. ACM, 34:596–615, 1987.
random networks. Science, 286(5439):509–512, 1999. [24] M.L. Fredman and D.E. Willard. Trans-dichotomous
[8] J.W. Berry, B. Hendrickson, S. Kahan, and P. Konecny. algorithms for minimum spanning trees and shortest
Graph software development and performance on the paths. J. Comput. System Sci., 48:533–551, 1994.
MTA-2 and Eldorado. In Proc. Cray User Group [25] L.C. Freeman. A set of measures of centrality based
meeting (CUG 2006), Lugano, Switzerland, May 2006. on betweenness. Sociometry, 40(1):35–41, 1977.
CUG Proceedings. [26] A.M. Frieze and L. Rudolph. A parallel algorithm for
[9] U. Brandes. A faster algorithm for betweenness cen- all-pairs shortest paths in a random graph. In Proc.
trality. J. Mathematical Sociology, 25(2):163–177, 22nd Allerton Conference on Communication, Control
2001. and Computing, pages 663–670, 1985.
[10] P. Briggs and L. Torczon. An efficient representation [27] G. Gallo and P. Pallottino. Shortest path algorithms.
for sparse sets. ACM Lett. Program. Lang. Syst., 2(1- Ann. Oper. Res., 13:3–79, 1988.
4):59–69, 1993. [28] F. Glover, R. Glover, and D. Klingman. Computa-
[11] G.S. Brodal, J.L. Träff, and C.D. Zaroliagis. A parallel tional study of an improved shortest path algorithm.
priority queue with constant time operations. Journal Networks, 14:23–37, 1984.
of Parallel and Distributed Computing, 49(1):4–21, [29] A.V. Goldberg. Shortest path algorithms: Engineering
1998. aspects. In ISAAC 2001: Proc. 12th Int’l Symp. on
[12] A. Broder, R. Kumar, F. Maghoul, P. Raghavan, Algorithms and Computation, pages 502–513, London,
S. Rajagopalan, R. Stata, A. Tomkins, and J. Wiener. UK, 2001. Springer-Verlag.
Graph structure in the web. Computer Networks, 33(1- [30] A.V. Goldberg. A simple shortest path algorithm with
6):309–320, 2000. linear average time. In 9th Ann. European Symp.
[13] D. Chakrabarti, Y. Zhan, and C. Faloutsos. R-MAT: A on Algorithms (ESA 2001), volume 2161 of Lecture
recursive model for graph mining. In Proc. 4th SIAM Notes in Computer Science, pages 230–241, Aachen,
Intl. Conf. on Data Mining (SDM), Orlando, FL, April Germany, 2001. Springer.
2004. SIAM. [31] D. Gregor and A. Lumsdaine. Lifting sequential graph
algorithms for distributed-memory parallel computa- tributed Processing Symp. (IPDPS), pages 1–8, Fort
tion. In Proc. 20th ACM SIGPLAN Conf. on Object Lauderdale, FL, April 2002. IEEE Computer Society.
Oriented Programming, Systems, Languages, and Ap- [46] U. Meyer. Design and Analysis of Sequential and
plications (OOPSLA), pages 423–437, New York, NY, Parallel Single-Source Shortest-Paths Algorithms. PhD
USA, 2005. ACM Press. thesis, Universität Saarlandes, Saarbrücken, Germany,
[32] R. Guimerà, S. Mossa, A. Turtschi, and L.A.N. October 2002.
Amaral. The worldwide air transportation network: [47] U. Meyer. Average-case complexity of single-source
Anomalous centrality, community structure, and cities’ shortest-paths algorithms: lower and upper bounds. J.
global roles. Proceedings of the National Academy of Algs., 48(1):91–134, 2003.
Sciences USA, 102(22):7794–7799, 2005. [48] U. Meyer and P. Sanders. Parallel shortest path for
[33] T. Hagerup. Improved shortest paths on the word arbitrary graphs. In Proc. 6th International Euro-Par
RAM. In 27th Colloquium on Automata, Languages Conference (Euro-Par 2000), volume 1900 of Lecture
and Programming (ICALP), volume 1853 of Lecture Notes in Computer Science, pages 461–470, Munich,
Notes in Computer Science, pages 61–72, Geneva, Germany, 2000. Springer-Verlag.
Switzerland, 2000. Springer-Verlag. [49] U. Meyer and P. Sanders. ∆-stepping: a parallelizable
[34] Y. Han, V. Pan, and J. Reif. Efficient parallel algo- shortest path algorithm. J. Algs., 49(1):114–152, 2003.
rithms for computing the all pair shortest paths in di- [50] M.E.J. Newman. Scientific collaboration networks: II.
rected graphs. Algorithmica, 17(4):399–415, 1997. shortest paths, weighted networks and centrality. Phys.
[35] M.R. Hribar and V.E. Taylor. Performance study of Rev. E, 64:016132, 2001.
parallel shortest path algorithms: Characteristics of [51] M.E.J. Newman. The structure and function of com-
good decomposition. In Proc. 13th Ann. Conf. Intel plex networks. SIAM Review, 45(2):167–256, 2003.
Supercomputers Users Group (ISUG), 1997. [52] M. Papaefthymiou and J. Rodrigue. Implementing
[36] M.R. Hribar, V.E. Taylor, and D.E. Boyce. Paral- parallel shortest-paths algorithms. DIMACS Series
lel shortest path algorithms: Identifying the factors in Discrete Mathematics and Theoretical Computer
that affect performance. Report CPDC-TR-9803-015, Science, 30:59–68, 1997.
Northwestern University, Evanston, IL, 1998. [53] J. Park, M. Penner, and V.K. Prasanna. Optimizing
[37] M.R. Hribar, V.E. Taylor, and D.E. Boyce. Reducing graph algorithms for improved cache performance. In
the idle time of parallel shortest path algorithms. Proc. Int’l Parallel and Distributed Processing Symp.
Report CPDC-TR-9803-016, Northwestern University, (IPDPS 2002), Fort Lauderdale, FL, April 2002. IEEE
Evanston, IL, 1998. Computer Society.
[38] M.R. Hribar, V.E. Taylor, and D.E. Boyce. Termi- [54] R. Raman. Recent results on single-source shortest
nation detection for parallel shortest path algorithms. paths problem. SIGACT News, 28:61–72, 1997.
Journal of Parallel and Distributed Computing, 55:153– [55] H. Shi and T.H. Spencer. Time-work tradeoffs of the
165, 1998. single-source shortest paths problem. J. Algorithms,
[39] H. Jeong, S.P. Mason, A.-L. Barabási, and Z.N. Oltvai. 30(1):19–32, 1999.
Lethality and centrality in protein networks. Nature, [56] M. Thorup. Undirected single-source shortest paths
411:41–42, 2001. with positive integer weights in linear time. Journal of
[40] P.N. Klein and S. Subramanian. A randomized parallel the ACM, 46(3):362–394, 1999.
algorithm for single-source shortest paths. J. Algs., [57] J. L. Träff. An experimental comparison of two dis-
25(2):205–220, 1997. tributed single-source shortest path algorithms. Paral-
[41] F. Liljeros, C.R. Edling, L.A.N. Amaral, H.E. Stanley, lel Computing, 21(9):1505–1532, 1995.
and Y. Åberg. The web of human sexual contacts. [58] J. Ullman and M. Yannakakis. High-probability par-
Nature, 411:907–908, 2001. allel transitive closure algorithms. In Proc. 2nd Ann.
[42] K. Madduri. 9th DIMACS implementation challenge: Symp. Parallel Algorithms and Architectures (SPAA-
Shortest Paths. ∆-stepping C/MTA-2 code. http:// 90), pages 200–209, Crete, Greece, July 1990. ACM.∼ kamesh/research/DIMACS-ch9, [59] F.B. Zhan and C.E. Noon. Shortest path algorithms:
2006. an evaluation using real road networks. Transp. Sci.,
[43] K. Madduri, D.A. Bader, J.W. Berry, and J.R. Crobak. 32:65–73, 1998.
Parallel shortest path algorithms for solving large-
scale instances. Technical report, Georgia Institute of
Technology, September 2006.
[44] U. Meyer. Heaps are better than buckets: parallel
shortest paths on unbalanced graphs. In Proc. 7th
International Euro-Par Conference (Euro-Par 2001),
pages 343–351, Manchester, United Kingdom, 2000.
[45] U. Meyer. Buckets strike back: Improved parallel
shortest-paths. In Proc. 16th Int’l Parallel and Dis-