Families of Distributed Memory Parallel Graph Algorithms from Self-Stabilizing Kernels-An SSSP Case Study

TL;DR

This paper introduces the AGM model to transform self-stabilizing algorithms into scalable distributed graph algorithms, specifically for SSSP, demonstrating improved performance over standard methods.

Contribution

The paper presents the AGM model that converts self-stabilizing algorithms into distributed graph algorithms and extends it for architecture-specific optimization.

Findings

AGM model effectively converts self-stabilizing algorithms for distributed graph processing.

Derived SSSP variants outperform standard distributed algorithms in experiments.

Fine-grained orderings enable architecture-specific scalable algorithm generation.

Abstract

Self-stabilizing algorithms are an important because of their robustness and guaranteed convergence. Starting from any arbitrary state, a self-stabilizing algorithm is guaranteed to converge to a legitimate state.Those algorithms are not directly amenable to solving distributed graph processing problems when performance and scalability are important. In this paper, we show the "Abstract Graph Machine" (AGM) model that can be used to convert self-stabilizing algorithms into forms suitable for distributed graph processing. An AGM is a mathematical model of parallel computation on graphs that adds work dependency and ordering to self-stabilizing algorithms. Using the AGM model we show that some of the existing distributed Single Source Shortest Path (SSSP) algorithms are actually specializations of self-stabilizing SSSP. We extend the AGM model to apply more fine-grained orderings at different spatial levels to derive additional scalable variants of SSSP algorithms, essentially enabling the algorithm to be generated for a specific target architecture. Experimental results show that this approach can generate new algorithmic variants that out-perform standard distributed algorithms for SSSP.

Tables

Table 1. TABLE I : Mean timing results (in seconds) for buffer, threadq, nodeq and numaq when ran against listed real graphs.

Graph	$\left|V\right|$	$\left|E\right|$	$D$	Cores	AGM	buffer	${\sigma }_{buf}$	threadq	${\sigma }_{thread}$	nodeq	${\sigma }_{node}$	numaq	${\sigma }_{numa}$
SOC-Live [14]	4847571	68993773	16	64	$\Delta =3$	22.7	6.47	20.52	5.8	14.93	2.39	29.8	3.87
					$KLA,K=1$	90.39	5.78	38.66	4.33	25.08	1.94	35.96	2.55
					$Chaotic$	39.82	6.98	11.66	0.63	166.26	15.42	207.42	22.21
Wiki-Talk [15]	2394385	5021410	9	64	$\Delta =3$	2.27	0.28	2.26	0.57	3.44	0.28	9.71	0.5
					$KLA,K=1$	3.73	0.44	2.53	0.34	1.94	0.23	3.9	0.4
					$Chaotic$	41.46	6.54	1.41	0.05	8.34	1.29	5.97	0.78
roadNet-CA [14]	1965206	2766607	849	1024	$\Delta =1200$	24.53	6.49	22.28	6.09	12.23	1.2	19.18	2.05
					$KLA,K=10$	54.24	6.14	54.63	5.8	43.86	5.18	51.35	5.54
					$Chaotic$	44.62	7.4	2.68	0.21	44.17	2.08	23.76	2.1
Orkut [16]	3072441	117185083	9	1024	$\Delta =10$	3.72	0.43	3.12	0.4	4.29	0.36	6.77	0.3
					$KLA,K=5$	20.94	6.33	71.19	23.5	15.87	1.56	18.28	2.12
					$Chaotic$	64.84	12.19	2.97	1.14	56.81	5.43	50.41	3.99

Equations

$$\text{π}(w)=\bigcup _{s_i\in St}{s}_i(w)$$

$$s_i(w)=\{\begin{array}{l}\{w_{new}|w_{new}\in N(w)\\ \text{ if }\text{C (w) is True \& U(w) is True}\}\\ \{\}else\end{array}$$

$$\text{π}(w)=\{\begin{array}{l}\{w_k|w_k\in <N(w)>,\\ <U(w)>,\\ <C(w)>\}\end{array}$$

$${\text{π}}^{SSSP}(w)=\{\begin{array}{l}\{w_k|w_k\in <w_k[0]\in neighbors(w[0])and\\ w_k[1]=w[1]+weight(w[0],w_k[0])>,\\ <distance(w[0])⟵w[1]>,\\ <ifw[1]<distance(w[0])>\\ \}\end{array}$$

$${\text{π}}^{kla}(w)=\{\begin{array}{l}\{w_k|w_k\in <w_k[0]\in neighbors(w[0])and\\ w_k[1]=w[1]+weight(w[0],w_k[0])\\ and{w}_k[2]=w[2]+1>,\\ <distance(w[0])⟵w[1]>,\\ <ifw[0]<distance(w[0])>\}\end{array}$$

Full text

I Introduction

Most of the existing parallel algorithms are developed based on Parallel Random Access Machine (PRAM) [1] model. While PRAM is a simple machine model it does not consider factors that are significant in distributed memory system (e.g., overhead of synchronization, remote message communication etc.). Further, algorithms designed for PRAM may assume global data structures and subgraph computations are efficient, hence their cost is not counted into algorithm performance. While these algorithms perform well in shared memory systems they tend to perform poorly on distributed memory systems.

Self-stabilizing graph algorithms rely on local information to solve graph problems. In self-stabilizing algorithms every vertex in the graph is associated with a state. Whenever, there is a state change, neighboring vertices are notified via edges. Self-stabilizing algorithms does not use global data structures and does not rely on operations such as subgraph computations (these operations are expensive in a distributed environment). The fact that, self-stabilizing algorithms rely only on local information and does not assume global data structures motivate us to investigate the applicability of self-stabilizing graph algorithms for large scale static graph processing.

A self-stabilizing algorithm consists of set of rules. Every rule has a condition. A rule is evaluated only if its condition evaluates to true. A self-stabilizing algorithm reaches a legitimate state irrespective of its initial state.

Self-stabilizing graph algorithms have been introduced for number of graph applications including Single Source Shortest Path, Breadth First Search (BFS), Spanning Tree Construction, Maximal Independent Set and Graph Coloring. A survey of self-stabilizing graph algorithms is presented in [2].

The strategy for updating vertex states is defined by an execution model. In self-stabilization, those execution models are called demons. We find three types of demons in self-stabilizing algorithms: central demon, synchronous demon and distributed demon. In a central demon algorithm, only one vertex can update the state at a time. While synchronous demon updates all the vertex states at the same time, a distributed demon select a subset of vertices to update states at the same time. Since our main focus is to reduce global synchronization and to rely on “local” data for processing, we only consider distributed demon algorithms.

A distributed demon self-stabilizing Single Source Shortest Path algorithm is presented in Algorithm 1 (discussed in  [3]). At stabilization this algorithm will have the minimum distance to every vertex, from a given source vertex. Self-stabilizing algorithms describe the algorithm using rules. A rule consists of a condition and an action, when condition evaluated true action is invoked. Rules in the Algorithm 1 have the format; current state $\longrightarrow$ new state.

Note that in the algorithm, $d(i)$ represents the state of vertex $i$ and $N(i)=\{i\in V|(i,j)\in E\}$ stands for set of all neighbours of vertex $i$. The pre assigned weight for an edge is denoted by $w(i,j)$.

The algorithm consists of two rules. The first rule is only applicable to the source. It says, if the source “distance state” is not 0, then the “distance state” should be transferred to 0. The second rule is activated only if the current distance of the vertex is not equal to the minimum of its neighbours distance plus the weight of the edge connecting the neighbour. So, the legitimate state of the system is $d(r)=0$ ($r$ is the given source) and $\forall i\neq r,d(i)=min_{j\in N(i)}(d(j)+w(i,j))$.  [3] proves that this algorithm ultimately stabilizes and $d(i)$ represents the shortest distance from the source, for each vertex $i$ at stabilization.

Another important aspect of self-stabilizing algorithms is the atomicity requirement. Algorithm 1 requires to querying vertex $i$’s state and its neighbors state and updating vertex $i$’s state (if the rule evaluates to true for vertex $i$) in a single atomic step. In self-stabilization, the requirement to query the current vertex state and its neighbours states and update the current vertex state in a single atomic step is called composite atomicity. The Read-Write atomicity is another form of atomicity discussed in detail in [4]. In this paper we consider self-stabilizing algorithms with composite atomicity.

Using self-stabilizing algorithms for static graph analysis is challenging due to several reasons: 1. Self-stabilizing algorithms are designed for streaming contexts in which those algorithms do not terminate after stabilizing; 2. To implement the composite atomicity requirement we need to lock the current vertex as well as its neighbors. This will generate a lot of synchronization regions and involves a significant number of locking/unlocking operations, hence it reduces performance especially in distributed execution; 3. A naive implementation of the self-stabilizing SSSP would require program to iterate through all the vertices and apply rules R0 and R1 until they reach the legitimate state. When processing a large graph, such an implementation would generate a massive amount of unnecessary work and would perform very poorly.

To overcome the challenges discussed above, we model self-stabilizing SSSP Algorithm 1 using an Abstract Graph Machine (AGM) [5]. An AGM is a model for designing distributed memory parallel graph algorithms. An AGM essentially converts the self-stabilizing algorithm to a stabilizing algorithm by adding work dependency and uses ordering to reduce the amount of work. A stabilizing algorithm starts from a specific initial state where as self-stabilizing algorithm can execute the algorithm from any initial state. The modeled stabilizing algorithm is in a suitable form for distributed graph processing. We also show that some of the existing well-known distributed SSSP algorithms are specific versions of the modeled self-stabilizing algorithm in 1.

We further enhance the AGM model to incorporate architecture dependent spatial characteristics to generate less synchronized orderings. The extended model is called Extended AGM (EAGM). Using EAGM we generate nine variations of SSSP algorithms and compare their performance to standard distributed SSSP algorithms under three different types of graph inputs. Our results show that some of the generated algorithm variations perform better compared to standard distributed SSSP algorithms.

II Self-Stabilizing SSSP & AGM

In the following we discuss several important aspects when modeling the self-stabilizing SSSP algorithm using AGM. More specifically, our discussion is focused on termination, composite atomicity and ordering.

Self-stabilizing algorithms are not iterative algorithms neither the AGM algorithms. While self-stabilizing algorithms does not terminate they go to an idle state when states are stabilized but AGM algorithms terminate at stabilization. Algorithm termination depends on the amount of active work available in the system. Our implementations use standard termination detection algorithms to count active work available in the system alive. As long as there are state changes there will be active work available. When active work is zero there are no more state changes and guarantees that the algorithm state is stabilized.

Under composite atomicity Algorithm 1 needs to query states from neighbors, query the current state and update current the vertex state in a single atomic step. This requires synchronization between neighboring vertices. However, for the SSSP algorithm in 1, the synchronization requirement can be alleviated due to the monotonic function “min”. In rule 1, neighbors states are processed in the “min” function. The “min” function selects the minimum distance + weight value for all neighbors, irrespective of the order states pushed in. Therefore we only need to maintain the state of the current vertex.

The AGM model makes state transitions based on work events. We define a unit of work to be a pair, $<v,s_{v}>$ where $v\in V$ and $s_{v}$ is the state associated with $v$. We call a unit of work a $\mathit{workitem}$ and we denote all the possible work items generated by an algorithm as the set $\mathit{WorkItem}$. In the AGM model, every-time a state is updated new work items are generated.

Before processing, new work items are ordered. In an AGM ordering is defined as a strict weak ordering. The strict weak ordering divide work items into different equivalence classes.

III Abstract Graph Machine (AGM)

An Abstract Graph Machine(AGM) consists of a definition of a $\mathit{WorkItem}$, an initial $\mathit{workitem}$ set, a set of states, a processing function and a strict weak ordering relation. The processing function takes a $\mathit{workitem}$ ($\in\mathit{WorkItem}$) as an input and may produce 0 or more work items. Further, the processing function may change values associated with states. The strict weak ordering relation order work items into a set of ordered (induced) equivalence classes.

The AGM model denotes a $\mathit{workitem}$ ($\in\mathit{WorkItem}$) as a tuple. A tuple’s first element stores a vertex and the rest of the positions store the state/s and ordering attribute values. For example, the self-stabilizing SSSP algorithm stores a vertex and a distance in a $\mathit{workitem}$. The $\mathit{workitem}$ values are populated by the processing function. The size of the tuple (i.e. the number of additional elements) is determined by the states in the algorithm and the ordering attributes used in the AGM formulation.

To access tuple elements we use the bracket operator; e.g., if $w\in\mathit{WorkItem}$ and if $w$ = <$v,p_{0},p_{1}\dots,p_{n}$>then $w[0]$ = v and $w[1]$ = $p_{0}$ and $w[2]$ = $p_{1}$, etc. The $\mathit{workitem}$ data (i.e. tuple elements) are read by the processing function. After reading values, the processing function can change the states associated with the vertex in the $\mathit{workitem}$.

An AGM maintains state values as mappings(functions). The domain of the state mappings is the set $V$. The co-domain depends on the possible values that can be held in states. For example, in the self-stabilizing SSSP algorithm the state mapping is $distance:V\longrightarrow\mathbb{R}$. In AGM terminology, accessing a state value associated with a vertex (or edge) “v” is denoted as “mapping_name[v]” (e.g., distance[v]).

In addition to state mappings, algorithms also use read-only properties. These properties usually hold graph data and are interpreted as functions, e.g., edge weights ($weights:E\longrightarrow\mathbb{R}$). In the abstraction we treat those read-only properties as part of the graph definition. In terms of the syntax, to read values, the AGM model does not distinguish between read-only properties and state mappings.

Both state mappings and properties are accessed within a processing function and the processing function only reads from the properties while processing function may modify state mappings when processing a $\mathit{workitem}$. The processing function updates the state value associated with a vertex in a single atomic step.

The logic inside the processing function is analogous to the code that runs infinitely on every process in a self-stabilizing algorithm. However, in the AGM model a single process may run multiple instances of processing functions (e.g., distributed shared memory model). The logic inside the processing function is based on the rules defined for self-stabilizing algorithms but is modified to work with work items. A processing function ($\pi$) takes a $\mathit{workitem}$ and may produce more work items based on the logic defined inside the $\pi$. In our formalism we treat states as side effects, in the sense they are not passed as explicit arguments but subject to change when executing $\pi$. We will denote set of states using $\mathit{Q}$ and we will use $\mathcal{P}(W)$ to denote the powerset of set $W$. Then, mathematically the $\pi$ is declared as \textit{\pi}:\mathit{WorkItem}\longrightarrow\mathcal{P}(\mathit{WorkItem}).

The processing function defines the basic logic of an algorithm. It consists of a set of statements ($St$). A statement specifies a condition based on input $\mathit{workitem}$ and/or states ($\mathit{\mathit{C}}:\mathit{WorkItem}\longrightarrow\{True,False\}$) and an update to states ($\mathit{\mathit{U}}:\mathit{WorkItem}\longrightarrow\{True,False\}$) and how a output work items ($w_{new}$) should be constructed ($\mathit{\mathit{N}}:\mathit{WorkItem}\longrightarrow\mathcal{P}(\mathit{WorkItem})$). The $\mathit{N}$ of a statement is evaluated only if its $\mathit{C}$ (condition) and $\mathit{U}$ are evaluates to $True$. We distinguish between $U$ and $C$, since $U$ may have side effects where it updates states but $C$ does not create any side effects.

An abstract version of $\pi$ is formally defined in Definition 1.

Definition 1

\mathit{\textit{\pi}}:\mathit{WorkItem}\longrightarrow\mathcal{P}(\mathit{WorkItem})**

[TABLE]

where; $\mathit{s_{i}}:\mathit{WorkItem}\longrightarrow\mathcal{P}(\mathit{WorkItem})$

[TABLE]

The output work items of a processing function are ordered according to the strict weak ordering defined on $\mathit{WorkItem}$. The ordered work items are again fed into the processing function in the order they appear. The interaction between the processing function and the ordering is graphically depicted in the Figure 1.

The strict weak ordering relation (denoted by $<_{\mathit{wis}}$) must satisfy the following properties;

For all $w\in\mathit{WorkItem}$, $w\nless_{\mathit{wis}}w$. 2. 2.

For all $w_{1},w_{2}\in\mathit{WorkItem}$, if $w_{1}<_{\mathit{wis}}w_{2}$, then $w_{2}\nless_{\mathit{wis}}w_{1}$. 3. 3.

For all $w_{1},w_{2},w_{3}\in\mathit{WorkItem}$, if $w_{1}<_{\mathit{wis}}w_{2}$ and $w_{2}<_{\mathit{wis}}w_{3}$, then $w_{1}<_{\mathit{wis}}w_{3}$. 4. 4.

For all $w_{1},w_{2},w_{3}\in\mathit{WorkItem}$, if $w_{1}$ not comparable with $w_{2}$ and $w_{2}$ not comparable with $w_{3}$, then $w_{1}$ is not comparable with $w_{3}$.

Properties 1 and 2 states that the strict weak ordering relation is not reflexive and is anti-symmetric. Property 3 denotes the transitivity of the “comparable work items” and Property 4 states that transitivity is preserved among non-comparable elements in the $\mathit{workitem}$ set. These properties give rise to an equivalence (i.e. non-comparable work items belong to the same equivalence class) relation defined on $\mathit{WorkItem}$, hence the strict weak ordering relation partitions the complete $\mathit{WorkItem}$. Since work items in different equivalence classes are comparable, the strict weak ordering relation defined on the set $\mathit{WorkItem}$ induces an ordering on generated equivalence classes. In general, there are several ways to define the induced ordering relation (denoted $<_{\mathit{WIS}}$). For our work we use the definition given in the Definition 2.

Definition 2

$<_{\mathit{WIS}}$* is a binary relation defined on $\mathcal{P}(\mathit{WorkItem})$, such that if $W_{1},W_{2}\in\mathcal{P}(\mathit{WorkItem})$ then; $W_{1}\leq_{\mathit{WIS}}W_{2}\;iff\\ \;forall\;w_{1}\in W_{1}\;and\;forall\;w_{2}\in W_{2};w_{1}<_{\mathit{wis}}w_{2}$.*

Having defined all supporting concepts we now give the definition of an AGM in Definition 3.

Definition 3

An Abstract Graph Machine(AGM) is a 6-tuple (G, $\mathit{WorkItem}$, Q, $\pi$, $<_{\mathit{wis}}$, S), where

G = (V, E) is the input graph, 2. 2.

\textbf{\mathit{WorkItem}}\subseteq(V\times P_{0}\times P_{1}\dots\times P_{n})* where each $P_{i}$ represents a state value or an ordering attribute,* 3. 3.

Q - Set of states represented as mappings, 4. 4.

\textit{\pi}:\mathit{WorkItem}\longrightarrow\mathcal{P}(\mathit{WorkItem})* is the processing function,* 5. 5.

$<_{\mathit{wis}}$* - Strict weak ordering relation defined on work items* 6. 6.

S ($\subseteq\mathit{WorkItem}$) - Initial $\mathit{workitem}$ set.

In the following we give the semantics of an AGM. An AGM starts execution with the initial $\mathit{workitem}$ set. The initial $\mathit{workitem}$ set is ordered according to the strict weak ordering relation. Next, the work items within the smallest equivalence class is fed to the processing function. If the processing function generates new work items, then they are separated into equivalence classes based on the strict weak ordering relation. The work items within a single equivalence class can execute the processing function in parallel. However, work items in two different equivalence classes must be ordered according to the induced relation (i.e. $<_{\mathit{WIS}}$). When executing work items in an equivalence class, it may generate new work items for the same equivalence class or to an equivalence class greater (as per $<_{\mathit{WIS}}$) than currently processing equivalence class. The AGM executes work items in the next equivalence class, once it finished executing all the work items in the current smallest equivalence class. An AGM terminates when it executes all the work items in all the equivalence classes.

The AGM is used to model graph algorithms related to graph applications such as SSSP, BFS, PageRank and Connected Components. Processing functions and orderings used in those algorithms is discussed detail in [5].

III-A SSSP Algorithms in AGM

In this subsection, we build AGM model for Algorithm 1. We also show that adding different orderings to the modeled algorithm reveals behaviours of existing distributed SSSP algorithms.

The SSSP algorithms discussed in this paper can be formulated using a single statement. For brevity we use following notation to represent the processing function.

Notation 1

\textit{\pi}:\mathit{WorkItem}\longrightarrow\mathcal{P}(\mathit{WorkItem})**

[TABLE]

In notation 1, $w_{k}$ is the new $\mathit{workitem}$ generated from $\mathit{N}$. As discussed previously $\mathit{U}$ and $\mathit{C}$ represents state update and condition.

To build the AGM for self-stabilizing  Algorithm 1, we need to define each tuple element in Definition 3 (i.e. (G, $\mathit{WorkItem}$, Q, $\pi$, $<_{\mathit{wis}}$, S)). As the input graph we use $G=(V,E,weight)$, where weight is a read-only property map that has weights associated to edges. As explained in Section II the set $\mathit{WorkItem}$ is defined based on the vertex state. For Algorithm 1, the state we are interested in is the distance from the source vertex. Therefore, we define, $\mathit{WorkItem}^{sssp}\subseteq(V\times Distance)$, where Distance $\subseteq\mathbb{R}_{+}^{*}$. The only state AGM needs is the distance from the source vertex and it is represented as distance. Values for distance state is assigned when the processing function is executed with work items.

The AGM processing function is defined based on the rules in Algorithm 1 and given in Definition 4. Since Rule 1 is only applied to the source vertex, we can move it to the initial $\mathit{workitem}$ set in the AGM representation. Rule 2 is encoded into the processing function in the format defined in Definition 1. The definition of the processing function uses a helper function called neighbors (Declared as $neighbors:V\longrightarrow\mathcal{P}(V)$), which operates on graph vertices.

Definition 4

\textit{\pi}^{sssp}:\mathit{WorkItem}^{sssp}\rightarrow\mathcal{P}(\mathit{WorkItem}^{sssp})**

[TABLE]

The next required parameter definition for the AGM model is the strict weak ordering relation ($<_{\mathit{wis}}$). If work items are not ordered in any form, then we will have a single large equivalence class of generated work items. We can have a single large equivalence class by defining strict weak ordering relation as in Definition 5.

Definition 5

$<_{\mathit{chaotic}}$* is a binary relation defined on $\mathit{WorkItem}^{sssp}$ where, $w_{1}<_{\mathit{chaotic}}w_{2}\;False\;\forall w_{1},w_{2}\in\mathit{WorkItem}^{sssp}$*

Basically, the binary relation $<_{\mathit{chaotic}}$ does not divide work items into any comparable equivalence classes. Using the strict weak ordering defined in Definition 5, we present the AGM model for Algorithm 1, in Proposition 1.

Proposition 1

$G=(V,E,weight)$* is the input graph,* 2. 2.

$\mathit{WorkItem}$* = $\mathit{WorkItem}^{sssp}$,* 3. 3.

Q = {distance} is the state mappings, 4. 4.

$\pi$* = \textit{\pi}^{SSSP},* 5. 5.

Strict weak ordering relation $<_{\mathit{wis}}$ = $<_{\mathit{chaotic}}$, 6. 6.

S = {<$v_{s}$, 0>}* where $v_{s}\in V$ and $v_{s}$ is the source vertex.*

Since the AGM presented in Proposition 1, does not perform any ordering on work items, we call it Chaotic AGM. However, the ordering on work items can be improved by defining strict weak ordering relations that generate smaller comparable equivalence classes. There are numerous ways to define strict weak orderings so that they generate smaller equivalence classes. Some of those strict weak orderings yield us, existing distributed SSSP algorithms such as Dijkstra’s SSSP [6] algorithm, $\Delta$-stepping SSSP [7] algorithm and KLA SSSP algorithm [8]. All those algorithms share almost the same processing function but uses different orderings on work items.

III-A1 Dijkstra’s Algorithm

Dijkstra’s SSSP algorithm is the work efficient SSSP algorithm. Dijkstra’s algorithm globally orders vertices by their associated distances and the shortest distance vertices are processed first. In the following, we define the ordering relation for Dijkstra’s algorithm and we instantiate Dijkstra’s algorithm using Abstract Graph Machine.

Definition 6

$<_{\mathit{dj}}$* is a binary relation defined on $\mathit{WorkItem}^{sssp}$ as follows; Let $w_{1},w_{2}\in\mathit{WorkItem}^{sssp}$, then; $w_{1}<_{\mathit{dj}}w_{2}$ iff $w_{1}[1]<w_{2}[1]$*

It can be proved that $<_{dj}$ is a strict weak ordering that satisfy constraints listed under Definition 6 (proof is omitted to save space). AGM instantiation for Dijkstra’s algorithm is given in  Proposition 2.

Proposition 2

Dijkstra’s Algorithm is an instance of AGM where;

$G=(V,E,weight)$* is the input graph,* 2. 2.

$\mathit{WorkItem}$* = $\mathit{WorkItem}^{sssp}$,* 3. 3.

Q = {distance} is the state mappings, 4. 4.

$\pi$* = \textit{\pi}^{SSSP},* 5. 5.

Strict weak ordering relation $<_{\mathit{wis}}$ = $<_{\mathit{dj}}$, 6. 6.

S = {<$v_{s}$, 0>}* where $v_{s}\in V$ and $v_{s}$ is the source vertex.*

III-A2 $\Delta$-Stepping Algorithm

$\Delta$-Stepping [7] SSSP algorithm arranges vertex-distance pairs into distance ranges (buckets) of size $\Delta(\in\mathbb{N})$ and execute buckets in order. Within a bucket, vertex-distance pairs are not ordered, and can be executed in any order. Processing a bucket may produce extra work for the same bucket or for a successive buckets. The strict weak ordering relation for $\Delta$-stepping algorithm is given in Definition 7.

Definition 7

$<_{\Delta}$* is a binary relation defined on $\mathit{WorkItem}^{sssp}$ as follows; Let $w_{1},w_{2}\in\mathit{WorkItem}^{sssp}$, then;

$w_{1}<_{\Delta}w_{2}$ iff $\lfloor w_{1}[1]/\Delta\rfloor<\lfloor w_{2}[1]/\Delta\rfloor$*

Instantiation of $\Delta$-Stepping algorithm in the AGM is same as in Proposition 2, except the strict weak ordering relation is $<_{\Delta}$ (= $<_{\mathit{wis}}$).

III-A3 KLA SSSP Algorithm

The $k$-level asynchronous (KLA) paradigm [8] bridges level-synchronous and asynchronous paradigms for processing graphs. In level-synchronous approach (E.g:- In level-synchronous BFS), vertices are processed level by level. KLA processes vertices up-to $k$ levels asynchronously and then moves to next $k$ levels.

In the following we model KLA SSSP with the AGM. The KLA approach order work items by the level of the resulting tree. Therefore, we need an additional ordering attribute in the $\mathit{WorkItem}$ definition. The KLA $\mathit{WorkItem}$ is defined as $\mathit{WorkItem}^{kla}\subseteq V\times\mathit{Distance}\times\mathit{Level}$ where $Level\subseteq\mathbb{N}$.

Further, the processing function also needs to be altered to populate value for the new ordering attribute. The processing function for KLA SSSP is defined in Definition 8. The processing function generates work items with an updated level, while changing distance state appropriately.

Definition 8

\textit{\pi}^{kla}:\mathit{WorkItem}^{kla}\longrightarrow\mathcal{P}(\mathit{WorkItem}^{kla})**

[TABLE]

KLA SSSP order work items based on the $k$ level. The work items within two consecutive $k$ levels can be executed in parallel and work items that are not in two consecutive $k$ levels must be ordered. The strict weak ordering relation for KLA SSSP is given in Definition 9.

Definition 9

$<_{kla}$* is a binary relation defined on

$\mathit{WorkItem}^{kla}$ as follows:

Let $w_{1},w_{2}\in\mathit{WorkItem}^{kla}$, then

$w_{1}<_{kla}w_{2}$ iff $\lfloor w_{1}[2]/k\rfloor<\lfloor w_{2}[2]/k\rfloor$*

The AGM formulation for KLA SSSP algorithm is given in Proposition 3.

Proposition 3

KLA SSSP Algorithm is an instance of AGM where;

$G=(V,E,weight)$* is the input graph* 2. 2.

$\mathit{WorkItem}$* = $\mathit{WorkItem}^{t}$* 3. 3.

Q = {distance} are the state mappings, 4. 4.

$\pi$* = \textit{\pi}^{kla}* 5. 5.

Strict weak ordering relation $<_{\mathit{wis}}$ = $<_{kla}$ 6. 6.

S = {<$v_{s}$, 0, 0>}* where $v_{s}\in V$ and $v_{s}$ is the source vertex and level starts with 0.*

Each AGM discussed above divides $\mathit{WorkItem}$ into equivalence classes differently (Figure 2) The Dijkstra’s AGM (Figure 2(a)) generates an equivalence class per each different distance value. The work items that have the same distance belong to the same equivalence class while work items of different distances go into different equivalence classes. The $\Delta$-stepping (Figure 2(b)) AGM also performs ordering based on the distance but equivalence classes are generated based on a $\Delta$ value and in general have more elements compared to Dijkstra’s AGM. All the work items within a single equivalence class are guaranteed to have distances between $\Delta*i$ and $\Delta*(i+1)$ for some $i\in\mathbb{N}$. Similar to $\Delta$-stepping algorithms, KLA, too, generates larger equivalence classes but uses level as the ordering attribute. Figure 2(c), shows how KLA algorithm arranges equivalence classes. As shown in Figure 2(d), the Chaotic version has a single large equivalence class containing all the work items.

Each of the above algorithms is different because of the way they generate equivalence classes and how those equivalence classes are ordered. Otherwise, they tend to share the same processing function that implements Rule 1 in Algorithm 1. Further, if two SSSP algorithms share the same ordering attribute, then they share the same processing function in AGM model. For example, both Dijkstra’s AGM and $\Delta$-stepping share the same processing function. If two algorithms use different ordering attributes, then processing functions differ only to update values of different ordering attributes.

IV Extended-AGM

Ordering in terms of distance reduces the amount of redundant work in SSSP algorithms. Out of the algorithms discussed in the previous section, Dijkstra’s algorithm performs the best ordering, so that it does the minimum amount of redundant work. However, the overhead of ordering in Dijkstra’s algorithm is significant in a parallel distributed run-time due to the frequent synchronization. In other words, when the AGM instance generates more equivalence classes, the global synchronization overhead increases. The other algorithms discussed above reduce overhead of ordering by chunking work items into larger equivalence classes. This reduces the number of equivalence classes generated. The Extended-AGM adds ordering to chunks but to reduce the overhead of ordering at AGM level, it applies ordering at different lower spatial levels.

A spatial level defines the amount of memory accessible to order work items. The highest spatial level is the accumulation of all the memory of the participating nodes (also called global memory). The next spatial level is the memory available at a node. The memory in a single node can be further subdivided into logical regions such as numa domains. Each numa domain may be shared by several threads. The lowest spatial level is the thread local memory. Such a spatial division is highly architecture-dependent and hierarchical.

The EAGM depicts a spatial hierarchy as a tree (Figure 3). Every node in the spatial hierarchy has an ordering attached. The ordering attached to the root represents the ordering defined by the relevant AGM. The example given inFigure 3, is an EAGM hierarchy derived from the $\Delta$-stepping AGM. As shown inFigure 3 the $\Delta$ ordering ( Definition 7) is attached to the root node. The orderings attached to the lower spatial levels are performed on work items available to the memory in the relevant spatial level. Therefore, the orderings attached to lower spatial levels are more relaxed than the ordering attached to the root. By default the EAGM spatial hierarchy assumes Chaotic (i.e. no ordering) ordering but specific orderings can be enforced. The example in Figure 3 enforces strict weak ordering $<_{dj}$ at numa level and uses Chaotic orderings at other levels.

An Extended AGM is an AGM but instead of consisting a strict weak ordering relation an EAGM has a spatial hierarchy with annotated orderings. An EAGM extends an AGM if and only if the EAGM generates same equivalence classes in the AGM at the root level of the spatial hierarchy. Therefore, an AGM can have multiple EAGMs, where each EAGM has the same ordering as the AGM at the root but different orderings at lower spatial levels. Each EAGM represents a variation of the algorithm modeled in the relevant AGM. If an AGM generates equivalence classes with many work items, then the EAGM has provision to perform fine-grain ordering at different spatial levels. However, if the AGM ordering is such that it generates equivalence classes with few work items, then the derived EAGM has less opportunity to perform ordering at spatial levels. For example, the $\Delta$-stepping AGM generates equivalence classes with many work items (provided $\Delta$ is sufficiently large). Variations of the $\Delta$-stepping AGM can be generated by applying ordering to work items at the process level, the numa level or at the thread level. However, for Dijkstra’s AGM the spatial orderings may not improve the overall performance of the algorithm because the equivalence classes generated by Dijkstra’s AGM has fewer work items on average.

Out of the algorithms discussed in the previous section, the fine-grained spatial ordering is effective to AGMs defined for $\Delta$-stepping, KLA and Chaotic. By considering the spatial hierarchy used in Figure 3, we apply Dijkstra’s strict weak ordering relation ( Definition 6) to spatial hierarchy levels PROCESS, NUMA and THREAD to derive EAGMs (Figure 4). Each EAGM generates a variation of the main algorithm defined by its corresponding AGM.

Figure 4(a), shows EAGM variations derived for $\Delta$-stepping algorithm. Figure 4(a)-(i), applies $<_{dj}$ ordering to THREAD level and Figure 4(a)-(ii) applies $<_{dj}$ ordering to NUMA level and Figure 4(a)-(iii) applies ordering to PROCESS level. EAGMs for KLA and Chaotic are derived in the same way.

For convenience, we will refer to the original AGM implementation as buffer, the variation that does THREAD level ordering as threadq, the variation that does PROCESS level ordering as nodeq and the variation that does NUMA level ordering as numaq.

V Implementation

We implemented EAGMs shown in Figure 4. Each implementation generated a variation of main algorithm. For implementation we used a light weight active messaging framework based on MPI. To represent the local graph data, we used compressed sparse row format. 1D distribution is used to distribute the graph. States and read only mappings are maintained as property maps indexed by vertices or edges and property maps are also distributed. We used concurrent priority queues (flat combining synchronous priority queue [9]) for node level and numa level ordering.

VI Results

We experimented with the performance of EAGM variations on synthetic graphs and on real-world graphs. For synthetic graphs we used:

• •

RMAT1: Graphs based on the current Graph500 [10] BFS benchmark specification with R-MAT [11] parameters $A=0.57$, $B=C=0.19$ and $D=0.05$ and random edge weights from $1$ to $100$.

• •

RMAT2: Graphs generated based on the proposed Graph500 [12] SSSP benchmark specification with R-MAT parameters $A=0.50$, $B=C=0.1$ and $D=0.3$ and random edge weights from $1$ to $255.$

For real world graphs, we use four graphs with varying parameter from the SNAP [13] repository with random edge weights from $1$ to $100.$

All our experiments were carried out on a Cray XE6/XK7 supercomputer, with 32 AMD Opteron Abu Dhabi CPUs 64 GB memory per node (4 NUMA domains), running Cray Linux Environment and the cray-mpich2 ver. 7.2.5 MPI implementation. The run-time architecture we used, aligned with the spatial hierarchy used in our EAGM implementations (Figure 4).

VI-A Scaling Results

We ran weak scaling experiments on Graph500 graphs from scale 19 ($2^{19}$ vertices) to scale 31. We compare each AGM algorithm (i.e. main SSSP algorithms) to their EAGM variations (see section IV). The implementation of the main AGM algorithm is represented in buffer for each algorithm ($\Delta$-Stepping, KLA and Chaotic). The thread-level ordering variation, node level ordering variation and numa level ordering variation implementations are denoted using threadq, nodeq and numaq, for each algorithm.

VI-A1 $\Delta$-Stepping Variations

In fig. 5, we present weak scaling results for $\Delta$-stepping variations on RMAT1 and RMAT2 synthetic graphs with $\Delta$ values 3, 5, and 7. The original $\Delta$-stepping (buffer) algorithm performs the best in-node. Since no communication is involved, the additional ordering provided by the other implementations does not provide a sufficient benefit for its overhead.

In general, the threadq variation is the fastest in the distributed setting. The nodeq and the numaq variations perform better with increasing deltas, and they are not competitive with the buffer implementation.

In summary, while in-node performance is dominated by the traditional $\Delta$-stepping algorithm (aka the buffer implementation) the distributed execution shows significant improvement with threadq variation. Though, the numaq and nodeq variations should provide better performance than the threadq variation by providing more ordering, the overhead of the concurrent ordering reduces the performance of numaq and nodeq.

VI-A2 KLA Variations

KLA variations show different performance characteristics than $\Delta$-stepping. Weak scaling results for KLA implementations on RMAT1 and RMAT2 graphs $K=1,2$ and $3$ are shown in fig. 6. For KLA, the nodeq and the numaq variations perform the best at scale, with $K=1.$ At greater $K$ values, the performance of threadq is comparable to nodeq and numaq, but in absolute terms, the performance at higher $K$ values is worse than at $K=1.$ As said in the previous section, the numaq and nodeq provide the best potential ordering by ordering the most items. The overheads are kept at bay because at $K=1$ all the writes to the next level’s queue occur before all the reads. The flat combining queue [9] that we use performs the best in this scenario. For higher $K$ values, writes and reads get more mixed, and the advantage of numaq and nodeq becomes less pronounced (when $K$ is higher more work items go into queues and concurrent ordering overhead become significant). In KLA, for both RMAT1 and RMAT2 inputs, all EAGM variations (threadq, nodeq and numaq) perform better compared to original KLA algorithm (buffer).

VI-A3 Chaotic Variations

Figure 7 shows results for chaotic variations with RMAT1 and RMAT2 input graphs. The chaotic AGM has a single large equivalence class and does not perform any form of ordering. Due to work explosion we were unable to run the chaotic algorithm except for small scales. For the same reason both nodeq and numaq end up having larger queue sizes, hence the overhead of ordering became significant as we increase the scale. However, the thread level ordering shows good performance, specially in distributed execution. For RMAT2, threadq achieves almost perfect weak scaling. Furthermore, the threadq chaotic variation is faster than all other variations in terms of absolute performance, demonstrating how the structured (E)AGM approach may result in new, highly performant algorithms.

VI-B Real World Graphs

The real world graphs we used in our experiments are listed in Table I along with their characteristics and relevant results. Most of the real world graphs are fairly small. For our experiments, we pick either 64 or 1024 cores, depending on the size of the graph. SOC-Live Journal [14] and Wiki-Talk [15] experiments are run on 64 cores, California Road Network [14], and Orkut [16] experiments are run on 1024 cores.

For Live Journal, nodeq showed better results out of the $\Delta$-stepping variations and the $KLA$ variations. Similarly to synthetic graph results, the chaotic variation show best results for threadq. KLA nodeq and chaotic threadq showed the best results on Wiki-Talk graph. California Road Network has the highest diameter, with the edge weights ranging from 0 to 100. Both $\Delta$-stepping and KLA show good results with the node level ordering and Chaotic showed good results for threadq. Orkut graph input shows minimum timing values for $\Delta$-stepping in their threadq implementations. In Chaotic variations, threadq implementation continued to perform better for social networking graphs.

VII Related Work

SSSP is a classic example of an irregular application. Parallel graph algorithms for SSSP being well-studied. $\Delta$-Stepping [7], KLA [8], Bellman Ford [17], Crauser’s SSSP [18] are popular algorithms that address distributed memory parallel SSSP problem.

A distributed SSSP algorithm connecting self-stabilizing was studied in [19]. The same algorithm is discussed for different run-time characteristics in [20]. Algorithm discussed in [19] is an EAGM instantiation of Chaotic AGM. Much of the work related to self-stabilization is already discussed in Section I.

AGM is an abstract model for designing distributed memory parallel graph algorithms. Most of the graph algorithms existing today are designed based on PRAM architecture. In PRAM we have a single shared memory and individual processors reading/writing from/to shared memory.

As discussed in Section I, PRAM algorithms suffer from performance issues in distributed memory settings. Bulk Synchronous Parallel (BSP) [21] is a model used for designing distributed memory parallel algorithms. In BSP we have super steps where in each super step we do computation, communication and barrier synchronization. Also, there are certain variations of BSP where computation and communication are overlapped to improve performance yet uses barriers. BSP is more like an extended version of PRAM and hence the graph algorithms designed for PRAM can be implemented in distributed settings using BSP. In addition, there are models that consider network parameters and considers communications (e.g., LogP [22] and its descendants).

Regarding spatial orderings, Galois scheduler; OBIM [23] considers spatial features when processing an irregular application but it is an implementation for shared memory systems rather than an abstract model.

VIII Conclusion

Using the AGM abstraction, we converted Algorithm 1 to a form suitable for distributed memory, parallel graph processing. We showed that existing distributed graph algorithms; Dijkstra’s SSSP, $\Delta$-stepping SSSP and KLA SSSP are variations of the converted algorithm. Those algorithms basically implement the Rule 2 of Algorithm 1 with different ordering based on distance state or based on a different ordering attribute such as level. We also showed, proposed EAGM model can generate more fine-grained orderings at less synchronized spatial levels. Results of our experiments showed that some of the generated algorithms perform better compared to standard distributed memory, parallel SSSP algorithms under different graph inputs.

Acknowledgment

This work is supported by the National Science Foundation under Grant No. 1319520. Further, this research was supported in part by Lilly Endowment, Inc., through its support for the Indiana University Pervasive Technology Institute, and in part by the Indiana METACyt Initiative. The Indiana METACyt Initiative at IU was also supported in part by Lilly Endowment, Inc.