A 2D Parallel Triangle Counting Algorithm for Distributed-Memory Architectures

TL;DR

This paper introduces a scalable 2D parallel triangle counting algorithm for distributed-memory systems, optimizing communication and computation to handle large graphs efficiently.

Contribution

The paper presents a novel 2D cyclic decomposition-based distributed triangle counting algorithm with key optimizations for sparsity and reduced memory overhead.

Findings

Achieves 3.24 to 7.22 speedup over baseline with 169 MPI ranks

Attains an average speedup of 10.2 times compared to previous algorithms

Effectively handles large synthetic and real-world graphs

Abstract

Triangle counting is a fundamental graph analytic operation that is used extensively in network science and graph mining. As the size of the graphs that needs to be analyzed continues to grow, there is a requirement in developing scalable algorithms for distributed-memory parallel systems. To this end, we present a distributed-memory triangle counting algorithm, which uses a 2D cyclic decomposition to balance the computations and reduce the communication overheads. The algorithm structures its communication and computational steps such that it reduces its memory overhead and includes key optimizations that leverage the sparsity of the graph and the way the computations are structured. Experiments on synthetic and real-world graphs show that our algorithm obtains an average relative speedup that range between 3.24 and 7.22 out of 10.56 across the datasets using 169 MPI ranks over the performance achieved by 16 MPI ranks. Moreover, we obtain an average speedup of 10.2 times on comparison with previously developed distributed-memory parallel algorithms.

Tables

Table 1. Table 1. Datasets used in the experiments.

Graph	#vertices	#edges	#triangles
twitter (Kwak et al., 2010)	41,652,230	1,202,513,046	34,824,916,864
friendster (Samsi et al., 2017)	119,432,957	1,799,999,986	191,716
g500-s26 (Nelson et al., 2014)	67,108,864	1,073,741,824	49,158,464,716
g500-s27 (Nelson et al., 2014)	134,217,728	2,147,483,648	106,858,898,940
g500-s28 (Nelson et al., 2014)	268,435,456	4,294,967,296	231,425,307,324
g500-s29 (Nelson et al., 2014)	536,870,912	8,589,934,592	499,542,556,876

Table 2. Table 2. Parallel Performance achieved using 16 - 169 MPI ranks.

		#nodes	expected		ppt	ppt		tct	tct		overall	overall
datasets	ranks	used	speedup		time	speedup		time	speedup		runtime	speedup
g500-s28	16	8			151.71			576.30			728.23
	25	7	1.56		113.80	1.33		408.21	1.41		522.31	1.39
	36	9	2.25		76.19	1.99		291.11	1.98		367.47	1.98
	49	9	3.06		63.96	2.37		222.64	2.59		286.78	2.54
	64	7	4.00		56.40	2.69		221.65	2.60		278.26	2.62
	81	9	5.06		46.12	3.29		160.30	3.60		206.61	3.52
	100	10	6.25		43.12	3.52		136.90	4.21		180.18	4.04
	121	11	7.56		40.22	3.77		121.33	4.75		161.67	4.50
	144	18	9.00		34.35	4.42		105.60	5.46		140.05	5.20
	169	17	10.56		30.69	4.94		79.82	7.22		110.51	6.59
g500-s29	16	16			323.50			1371.75			1695.57
	25	25	1.56		158.90	2.04		731.32	1.88		890.56	1.90
	36	18	2.25		163.41	1.98		697.33	1.97		861.21	1.97
	49	17	3.06		126.98	2.55		510.08	2.69		637.41	2.66
	64	16	4.00		100.59	3.22		473.49	2.90		574.34	2.95
	81	17	5.06		88.60	3.65		386.65	3.55		475.55	3.57
	100	17	6.25		72.13	4.48		280.05	4.90		352.36	4.81
	121	21	7.56		65.69	4.92		250.37	5.48		316.24	5.36
	144	24	9.00		63.06	5.13		232.09	5.91		295.33	5.74
	169	29	10.56		53.54	6.04		191.16	7.18		244.81	6.93
twitter	16	2			60.76			109.46			170.45
	25	2	1.56		39.59	1.53		64.73	1.69		104.50	1.63
	36	2	2.25		39.63	1.53		61.33	1.78		101.17	1.68
	49	3	3.06		33.45	1.82		45.31	2.42		79.04	2.16
	64	3	4.00		30.16	2.01		42.13	2.60		72.48	2.35
	81	4	5.06		29.08	2.09		30.46	3.59		59.68	2.86
	100	5	6.25		32.74	1.86		30.81	3.55		63.70	2.68
	121	6	7.56		32.64	1.86		24.75	4.42		57.50	2.96
	144	7	9.00		33.36	1.82		25.39	4.31		58.85	2.90
	169	8	10.56		31.62	1.92		18.52	5.91		50.29	3.39
friendster	16	3			91.54			95.41			187.21
	25	2	1.56		57.84	1.58		71.82	1.33		129.78	1.44
	36	3	2.25		48.51	1.89		64.29	1.48		112.98	1.66
	49	5	3.06		36.47	2.51		46.75	2.04		83.37	2.25
	64	4	4.00		35.80	2.56		45.61	2.09		81.66	2.29
	81	5	5.06		33.24	2.75		35.36	2.70		68.78	2.72
	100	5	6.25		35.56	2.57		35.24	2.71		71.04	2.64
	121	7	7.56		29.51	3.10		27.51	3.47		57.09	3.28
	144	6	9.00		38.62	2.37		37.65	2.53		76.53	2.45
	169	8	10.56		31.55	2.90		29.43	3.24		61.23	3.06

Table 3. Table 3. g500-s29 maximum runtime and load imbalance incurred per shift.

ranks	maximum runtime	average runtime	load imbalance
25	187.93	177.81	1.05
36	106.65	93.79	1.14

Table 4. Table 4. g500-s29 task count growth with respect to the number of ranks.

ranks	task	percent increase with respect
used	counts	to previous rank
$16$	$33907905131$
$25$	$42360246067$	$25\%$
$36$	$50801950709$	$20\%$

Table 5. Table 5. Comparisons with Havoq ’s  (Pearce, 2017 ) triangle counting runtime

	2core	directed wedge	our	speedup
dataset	time	counting time	runtime	obtained
g500-s26	1.59	239.64	20.35	11.9
g500-s27	3.37	576.45	41.93	13.7
g500-s28	7.32	1395.11	79.82	14.6
twitter	1.88	124.72	18.52	6.2
friendster	3.29	24.75	29.43	$-$

Table 6. Table 6. Twitter graph runtime contrasted against other distributed-memory triangle counting approaches.

algorithm	fastest runtime reported	cores used
Our work	51.7	169
AOP (Arifuzzaman et al., 2017)	564.0	200
Surrogate (Arifuzzaman et al., 2017)	739.8	200
OPT-PSP (Kanewala et al., 2018)	${23.14}^a$	2048

Equations

$$c_{i,j}=\sum _{k=j+1}^n{u}_{i,k}\times u_{j,k},$$

$$\sum _{a_{i,j}=1\wedge i<j}{c}_{i,j}.$$

$$c_{i,j}=\sum _{k=1}^n{u}_{i,k}\times l_{k,j},$$

$$C[U]=UL$$

$$C[U_{x,y}]=\sum _z{U}_{x,z}{L}_{z,y}.$$

$$C[U_{x,y}]=\sum _{z=0}^{\sqrt{p}-1}{U}_{x,(x+y+z)\%\sqrt{p}}{L}_{(x+y+z)\%\sqrt{p},y}.$$

$$T_{pre-processing}=p+\frac{m}{p}+\frac{n}{p}+\log (p)+d_{max}+d_{max}\log (p).$$

$$T_{triangle-counting}=d_{avg}\frac{n}{\sqrt{p}}\Bigg{(}\frac{d_{avg}}{\sqrt{p}}+1\Bigg{)}.$$

Full text

triangle counting, graph analytics, distributed-memory

††journalyear: 2019††copyright: acmcopyright††conference: 48th International Conference on Parallel Processing; August 5–8, 2019; Kyoto, Japan††booktitle: 48th International Conference on Parallel Processing (ICPP 2019), August 5–8, 2019, Kyoto, Japan††price: 15.00††doi: 10.1145/3337821.3337853††isbn: 978-1-4503-6295-5/19/08††ccs: High performance computing Parallel algorithms††ccs: High performance computing Graph analytics††ccs: High performance computing Triangle counting

1. Introduction

The use of graphs to model large scale real-world data is ubiquitous in our everyday lives. In order to analyze and study the relationships these graphs model, graph analytic operations such as finding patterns of interest, analyzing the community structure and connectivity, and determining influential entities in a given graph, are commonly used.

One such graph analytic operation is counting the number of triangles in a graph. The number of triangles in a graph is an important statistic that is used as an intermediary step in various applications. It is used in computing the clustering coefficient and the transitivity ratio of graphs, both of which are used in characterizing the tendency of the nodes in a graph to cluster together. Furthermore, the computations involved in triangle counting forms an important step in computing the $k$-truss decomposition of a graph, detecting community structures, studying motif occurrences, detecting spamming activities, and understanding the structure of biological networks (Girvan and Newman, 2002; Smith et al., 2017; Shrivastava et al., 2008; Watts and Strogatz, 1998).

In recent years, driven by the growing size of the graphs that needs to be analyzed, there has been significant research in improving the efficiency of parallel algorithms for computing the exact and approximate number of triangles. Parallel triangle counting algorithms have been specifically built for GPUs, external memory, shared-memory, and distributed-memory platforms (Green et al., 2014; Hu et al., 2018; Bisson and Fatica, 2018; Polak, 2016; Wang et al., 2016; Arifuzzaman et al., 2017; Parimalarangan et al., 2017; Shun and Tangwongsan, 2015; Yaşar et al., 2018). The shared-memory class of solutions are limited by the amount of memory that is available in a single processor, thus, limiting the size of the graphs that can be analyzed. Moreover, in many practical settings, such large graphs are stored in a distributed fashion in the aggregate memory that is available in a distributed-memory system. Being able to successfully analyze large graphs in such scenarios requires the need of developing distributed-memory algorithms for counting the number of triangles. However, despite the advantages of a distributed-memory system, these algorithms face higher costs as compared to shared-memory systems during communication and synchronization steps. Furthermore, distributed-memory graph algorithms also entail the problem of intelligent graph partitioning in order to reduce the costs involved in communicating with neighboring vertices.

We present an MPI-based distributed-memory algorithm for triangle counting using a set intersection based approach (Tom et al., 2017). The key difference between our algorithm and previously proposed approaches is that it utilizes a 2D decomposition of the data and associated computations, which increases the concurrency that can be exploited and reduces the overall communication cost. Furthermore, our algorithm moves the data among the processors by utilizing a sequence of communication steps that are similar to those used by Cannon’s parallel matrix-matrix multiplication algorithm. This ensures that our algorithm is memory scalable and faces low communication overhead. We also include key optimizations that leverage the sparsity of the graph and the way the computations are structured. Some of these optimizations include enumerating the triangles using an ordering that specifically leverages hash-maps in the set intersection computation, changing the hashing routine for vertices based on the density of its adjacency list, and eliminating unnecessary intersection operations.

We evaluate the performance of our algorithm on various real-world and synthetically generated graphs and compare it against other existing state-of-the-art approaches. Our experiments show that we obtain a relative speedup that range between $3.24$ and $7.22$ out of $10.56$ across the datasets using 169 MPI ranks over the performance achieved by 16 MPI ranks. Moreover, the performance of our parallel algorithm compares favorably against those achieved by existing distributed memory algorithms that rely on 1D decomposition (Pearce et al., 2013; Arifuzzaman et al., 2017; Kanewala et al., 2018).

The rest of the paper is organized as follows. Section 2 introduces the notation and definitions used in the paper, Section 3 details the necessary background required, Section 4 describes competing approaches we use, Section 5 explains our 2D parallel algorithm and finally, Section 6 details several experiments which demonstrate our algorithm’s speedup and scaling capabilities.

2. Definitions and notations

We will assume that the graphs that we operate on are simple and undirected and are represented using the standard $G=(V,E)$ notation. We will use $\operatorname{Adj}(v_{i})$ to denote the adjacency list of $v_{i}$; i.e., the set of vertices that are adjacent to $v_{i}$. We will use $d(v_{i})$ to denote the degree of $v_{i}$, i.e., $d(v_{i})=|\operatorname{Adj}(v_{i})|$. We will use $\%$ to indicate a mod operation and $\div$ to indicate a divide operation.

We will use $A$ to denote the $n\times n$ adjacency matrix of a symmetric $n$ vertex graph, in which $a_{i,j}=1$, if there is an edge between $v_{i}$ and $v_{j}$, and 0, otherwise. Furthermore, we will use $U$ and $L$ to denote the upper and the lower triangular portion of the adjacency matrix of $G$. A triangle is a set of three vertices $\{v_{i},v_{j},v_{k}\}$ if the edges $(v_{i},v_{j})$, $(v_{i},v_{k})$, and $(v_{j},v_{k})$ exist in $E$. The problem of triangle counting is to compute the total number of unique triangles in $G$.

Lastly, let there be $p$ processors in the system, which can be arranged in a square grid of the form $\sqrt{p}\times\sqrt{p}$. The processor in the $i$th row and $j$th column in this square grid is denoted by $P_{i,j}$.

3. Background

3.1. Triangle counting

Triangle counting algorithms iterate over all the edges in a graph and for each edge, count the number of triangles that this edge is a part of. To ensure that they do not count the same triangle multiple times, they impose strict vertex enumeration rules and ignore any triangles that do not satisfy those. There are two such rules, which are referred to as $\langle i,j,k\rangle$ and $\langle j,i,k\rangle$, where $i<j<k$. The $\langle i,j,k\rangle$ rule dictates that the algorithm starts from the first vertex ($i=1$) and iterates over the non-zero entries of $U$ in a row-wise fashion (i.e., edges of the form $(v_{i},v_{j})$) and considers only the triangle-closing vertices $v_{k}$, where $j<k$. In contrast, the $\langle j,i,k\rangle$ rule dictates that the algorithm also starts from the first vertex ($i=1$) and considers the triangle-closing vertices $v_{k}$, where $j<k$, but, iterates over the non-zero entries of $U$ in a column-wise fashion (i.e., edges of the form $(v_{j},v_{i})$).

Using the $\langle i,j,k\rangle$ rule, it is easy to show that the number of triangles, $c_{i,j}$, that result from $U$’s $(v_{i},v_{j})$ edge is

[TABLE]

and the total number of triangles in the graph is

[TABLE]

The computation associated with $c_{i,j}$ can be viewed as a set-intersection operation between the adjacency lists of $v_{i}$ and $v_{j}$, where $i<j<k$, and can be modeled as a list-based intersection and a map-based intersection. For every edge $(v_{i},v_{j})$ in $U$, the list-based intersection involves jointly traversing the adjacency lists of vertex $v_{i}$ and $v_{j}$, and finding common vertices, $v_{k}$, between them, such that $j<k$. In the map-based approach, we use an auxiliary data structure and hash $U_{i,*}$ (or, $U_{j,*}$) within that. This auxiliary data structure is then used to lookup the vertices of $U_{j,*}$ (or, $U_{i,*}$). Each successful lookup accounts for a complete triangle. On comparing the performance of the two approaches, we observe that map-based approaches are faster than list-based approaches. This is because for each $v_{i}$, the hash-map used to store the adjacency list of $v_{i}$ can be reused for all $v_{j}\in U_{i,*}$. The list-based and map-based methods are further detailed in (Parimalarangan et al., 2017; Shun and Tangwongsan, 2015; Tom et al., 2017).

Note that an alternative way of writing Equation 1 is

[TABLE]

since all entries of $L$ for $k<j$ are by definition 0. Given this, we can conceptually view the computations associated with triangle counting as that of computing certain entries of the matrix that results from multiplying matrices $U$ and $L$. Specifically, the entries that we care about are those that correspond to the non-zero entries of $U$. We will use $C[U]$ to denote the non-zero entries of U, and we will use

[TABLE]

to denote this operation of computing the entries of $C[U]$.

Furthermore, studies (Arifuzzaman et al., 2017) have shown that ordering the vertices in non-decreasing degree before the triangle counting step leads to lower runtimes. This optimization has been incorporated in most triangle counting algorithms that have been developed so far. By using such a degree-based ordering, the length of $v_{i}$’s adjacency list will tend to be smaller than that of $v_{j}$’s for every edge $(v_{i},v_{j})$ with $i<j$. Therefore, if $v_{j}$’s adjacency list is considerably larger than that of $v_{i}$’s, then it is beneficial to employ the $\langle j,i,k\rangle$ enumeration scheme and hash $v_{j}$’s adjacency list. In our previous work on triangle counting algorithms for shared-memory platforms (Tom et al., 2017), the $\langle j,i,k\rangle$ enumeration scheme along with map-based intersection approaches gave us faster runtimes as compared to other approaches. In order to use the $\langle j,i,k\rangle$ enumeration scheme, the non-zero entries of $L$ are used and $C[L]$ is used to denote the non-zero entries of $L$. As in Equation 4, $C[L]$ is equal to $UL$.

3.2. Cannon’s parallel matrix multiplication algorithm

Cannon’s algorithm is a widely used method for multiplying two matrices in parallel (Cannon, 1969). It is a dense, memory-efficient matrix multiplication algorithm that performs the multiplication operation in $\sqrt{p}$ shifts. Let $A$ and $B$ be the input matrices that need to be multiplied, and $C$, be the resulting matrix, i.e., $C=AB$. The algorithm assumes a 2D block distribution of the input matrices and distributes these blocks across the square ($\sqrt{p}\times\sqrt{p}$) processor grid such that processor $P_{i,j}$ locally stores the blocks $A_{i,\small(j+i\small)\%\sqrt{p}}$ and $B_{\small(i+j\small)\%\sqrt{p},j}$. The algorithm multiplies these local blocks first. The local blocks of $A$ and $B$ are then shifted left along the row and up along the column, respectively, and multiplied again. For dense matrix multiplication, every processor is busy with computations after each shift, and the number of communications is bounded at $\sqrt{p}$.

4. Related work

Distributed memory algorithms are based on concepts similar to the ones described in Section 3.

Pearce et al. (Pearce, 2017) developed a triangle counting application over HavoqGT (Pearce et al., 2013), an asynchronous vertex-centric graph analytics framework for distributed-memory. Their approach includes distributing the undirected graph using distributed delegates (Pearce, 2017), followed by a two-core decomposition of the graph which removes the vertices that cannot be a part of any triangle. After this, they proceed with reordering their vertices based on degree to compute wedges. These wedges, partitioned using 1D decomposition, are henceforth queried for closure to check the existence of a triangle.

Arifuzzaman et al. (Arifuzzaman et al., 2017) developed two different approaches, one which avoids communication with overlapping partitions whereas, the other which optimizes on memory usage. In the communication-avoiding approach, the vertices of the graph are partitioned into $p$ disjoint subsets, where, $p$ is the number of processors. Each processor is responsible for storing one of these subsets and their corresponding adjacency lists. Additionally, in order for each processor to work independently, the adjacency lists of the adjacent vertices of these vertices are stored too. Since most real-world sparse graphs follow a power-law degree distribution, a naive partitioning of the vertices of such a graph will lead to high memory overheads as the lengths of adjacency lists will be very skewed. Furthermore, the triangle counting operation will also incur high load imbalance, which will negatively impact the performance. Arifuzzaman et al. have explored these aspects and have developed various partitioning schemes in order to load balance their computations. In order to reduce the memory overheads of the above approach, Arifuzzaman et al. have further developed a space-efficient method, which involves higher communication costs. In this approach, they partitioned the vertices across processors into disjoint subsets and only stored the adjacency list of these vertices. Subsequently, only one copy of the graph exists across all the processors. For every intersection operation, they follow a push-based mechanism, in which the adjacency list of a vertex is sent to processors which require this particular list for performing the intersection. However, this leads to high communication overheads.

Kanewala et al. (Kanewala et al., 2018) describes a distributed, shared-memory algorithm to triangle counting. Their algorithm explores different combinations of the upper triangular part of the adjacency matrix and the lower triangular part of the adjacency matrix to perform the set intersection operations (Refer to Section 3) between the adjacency lists. In order to parallelize the algorithm in a distributed setting, they perform a 1D decomposition of the adjacency matrix and send the adjacency list of a vertex to the rank which stores the adjacency lists of its adjacent vertex. However, in order to curb the number of messages generated, they block vertices and their adjacency lists and process them in blocks.

5. Methods

This section presents our parallel algorithm for triangle counting in distributed memory parallel systems. Our algorithm reads the input graph, preprocesses it to reorder the vertices in non-decreasing degree among other operations, and stores the graph using compressed spare row (CSR) format prior to triangle counting. Our implementation is based on the map-based triangle counting approach outlined in (Tom et al., 2017), which was briefly described in Section 3.

5.1. Parallelization

Task Decomposition and Data Distribution

In our algorithm, we treat the computations required for computing an entry of the $C[U]$ matrix (Equation 3) as an indivisible computational task. We decompose the computations among the $p$ processors, by mapping the tasks associated with $C[U]$ using a 2D decomposition of $C$. Specifically, the processors are organized in the form of a $\sqrt{p}\times\sqrt{p}$ processor grid and each processor is responsible for the elements of $C[U]$ that exist in the $n/\sqrt{p}\times n/\sqrt{p}$ entries of $C$ that were assigned to it.

However, as we only consider the upper triangular matrix, a naive 2D block partitioning will lead to load imbalance. Moreover, as the vertices are sorted in non-decreasing degree, the length of the adjacency list increases as the vertex id increases. Therefore, the tasks associated with the extreme right and the lower part of the $C[U]$ matrix will be computationally more expensive as they employ such vertices for the intersection. This further contributes to the load imbalance.

To address both issues and evenly distribute the work between the processors, we perform a cyclic distribution of $C$ over the $\sqrt{p}\times\sqrt{p}$ processor grid. Because of the degree-based ordering, successive rows/columns in the upper and lower triangular portions of the adjacency matrix will have similar number of non-zeros. Consequently, a cell-by-cell cyclic distribution will tend to assign a similar number of non-zeros (tasks) of $C[U]$ and at the same time, a similar number of light and heavy tasks to each processor.

Furthermore, in order to map the input blocks $U$ and $L$ with the tasks owned by the processors, we decompose the matrices using cyclic distribution over the processor grid. After the decomposition, as the number of vertices assigned to a processor is not contiguous anymore, the adjacency list of a vertex $v_{i}$ is accessed using the transformed index $v_{i}\div\sqrt{p}$ in the per-processor CSR representation. Let the blocks $U_{x,y}$ and $L_{x,y}$ be the respective decomposition of $U$ and $L$ over the $\sqrt{p}\times\sqrt{p}$ grid such that processor $P_{x,y}$ is responsible for those blocks.

Orchestration of Computation and Communication

Consider the set of tasks $C[U_{x,y}]$ that processor $P_{x,y}$ is responsible for. For every task $(v_{i},v_{j})$ in that block, in order to apply Equation 4, it requires the adjacency lists of the set of $v_{i}$ and $v_{j}$ vertices. Thus, $P_{x,y}$ needs the blocks $U_{x,*}$ and $L_{*,y}$ to determine the number of triangles in $C[U_{x,y}]$, and following the convention of Equation 4, the associated computations can be written as

[TABLE]

Though this can be done by having each processor first collect the necessary rows and column blocks of matrices $U$ and $L$, respectively, and then proceed to perform the required computations, such an approach will increase the memory overhead of the algorithm. We address this problem by realizing that the above summation can be appropriately performed by the communication pattern of Cannon’s 2D parallel matrix multiplication, and we utilize the same pattern in our algorithm (refer Section 3). In terms of this communication pattern, Equation 2 can be rewritten as (recall that $\%$ indicates a mod operation)

[TABLE]

In each shift, for each non-zero element $(v_{i},v_{j})$ that exists in block $C[U_{x,y}]$, we hash the row $v_{i}$ that exists in block $U_{x,(x+y+z)\%\sqrt{p}}$ and lookup the $v_{k}$ vertices that exist in the column $v_{j}$ in block $L_{{(x+y+z)\%{\sqrt{p}}},y}$ to find the number of triangles incident on the edge $(v_{i},v_{j})$.

The initial shifts of Cannon’s algorithm sends the block $U_{x,y}$ to processor $(x,(y+x)\%\sqrt{p})$ and the block $L_{x,y}$ to processor $((x+y)\%\sqrt{p},y)$. After performing the triangle counting operation on the blocks associated with that shift, the block $U_{x,y}$ is sent left to $P_{x,y-1}$ and the block of $L_{x,y}$ is sent up over to $P_{x-1,y}$ in the next $\sqrt{p}-1$ shifts, and the triangle counting operation is performed as before. Every processor accumulates this count of triangles corresponding to the tasks stored in $C[U_{x,y}]$ over the $\sqrt{p}$ shifts and, the same is globally reduced over all the processors in the grid in the end.

Finally, as discussed in Section 3, to leverage the benefits of enumerating the triangles using the $\langle j,i,k\rangle$ scheme in the map-based triangle counting approach, we define a task based on the non-zero elements in $L$ instead of $U$, as $L$ contains the incidence list for each vertex $v_{j}$. Therefore, $L$, instead of $U$, is cyclically distributed to construct a task block, denoted by $C[L_{x,y}]$.

5.2. Optimizations

We include several optimizations which leverage the characteristics of sparse graphs to further increase the performance of our distributed-memory algorithm. These are detailed below.

Modifying the hashing routine for sparser vertices

Due to the 2D decomposition and the fact that we perform the required computations by operating on blocks of $U$ and $L$, the lengths of the adjacency lists that are being intersected will tend to be smaller (on the average, they should be smaller by a factor of $\sqrt{p}$). A direct consequence of this is that even with a moderately sized hashmap, the number of collisions will tend to be smaller. In order to take advantage of this, before we hash the adjacency of a vertex within the triangle counting routine, we heuristically determine if the vertex is involved in collisions by utilizing the length of the adjacency list of the vertex. If the length of the adjacency list is less than the size of the hash-map, then those vertices will face no collision while being hashed. Such vertices are hashed by performing a direct bitwise AND operation without involving any probing.

Doubly sparse traversal of the CSR structure

As we are performing a block cyclic distribution of the upper and the lower triangular portion of the adjacency matrix, multiple vertices allocated to a processor may not contain any adjacent vertices. This is because each processor will have roughly $1/\sqrt{p}$ of the adjacency lists and if the vertices in the adjacency list has a degree in $U$ that is less than the degree of the vertex itself, then the adjacency list of that vertex could be rendered empty. However, these vertices can not be directly eliminated from the CSR structure due to the indexing scheme we use to avoid maintaining offsets. Therefore, in order to eliminate this unnecessary looping over these vertices while performing the set intersection operation, we use a data structure that is inspired by the doubly compressed sparse row structure (Buluç and Gilbert, 2010) to store the task matrix, as well as the upper and the lower triangular portion of the adjacency matrix. In our algorithm, while creating the CSR structure for each processor, we also associate with it a list of vertices that contain non-empty adjacency lists. We use this list of vertices to index into the CSR structure; thus, we avoid any vertices that have empty adjacency lists.

Eliminating unnecessary intersection operations

While performing the intersection operation between the two adjacency lists, we only need to consider those triangle-closing vertices, $k$, that satisfy $k>j$. However, as the adjacency lists are split over $\sqrt{p}$ processors, it is possible that some of the participating vertices are assigned to other processors, and performing an intersection operation with these entries will result in no triangles. Therefore, to overcome this problem, while performing the hashmap lookups, we traverse row $v_{i}$ backwards, that is, from the last adjacent vertex of $v_{i}$, and then break out of the loop as soon as we encounter an adjacent vertex id that is lesser than the last adjacent vertex id in column $v_{j}$. With this, we weed out all the unnecessary intersection operations that would otherwise result in no common vertices. Since the adjacency lists stored in the CSR structure are not required to be sorted for any given vertex, this optimization requires an initial sort before we start the triangle counting phase. However, the cost incurred by this sorting is amortized over the many set intersection operations that take place in this phase.

Reducing overheads associated with communication

In order to eliminate the cost of serializing and deserializing memory during MPI communication steps, we allocate the memory associated with all of the information for a sparse matrix as a single blob, and “allocate” from within that blob the various arrays of the sparse matrix that are required for the processing of the algorithm. Specifically, we convert and store the blocks $U_{x,y}$ and $L_{x,y}$ as a blob of bytes before the shifts begin. This gives us some savings with respect to the amount of time spent in communication.

5.3. Preprocessing

Our algorithm assumes that the graph is initially stored using a 1D distribution, in which each processor has $n/p$ vertices and its associated adjacency lists. Given that, it proceeds to perform a series of preprocessing steps whose goals are to (i) perform a cyclic distribution of the graph followed by a relabeling of the vertices, (ii) reorder the vertices of the graph in non-decreasing degree, (iii) redistribute the graph based on the 2D cyclic distribution required by our parallel formulation, and, (iv) create the upper and the lower triangular portion of the adjacency matrix. The rational behind these pre-processing steps and our corresponding parallel formations are described in the rest of this section.

Initial redistribution

In some cases the initial distribution of the graph may be such that even though each processor had $n/p$ vertices, it may have adjacency lists with significant variation in the number of non-zeros elements. In order to reduce the load imbalance incurred while preprocessing the graph datasets that contain localized sets of highly dense vertices, we perform an initial cyclic distribution of the graph and relabel the vertices in the adjacency list accordingly.

Reorder vertices based on non-decreasing degree

Recall from the discussion in Section 3, prior to triangle counting, the graph vertices are reordered in non-decreasing degree as this significantly reduces the time required to do triangle counting. In order to achieve the same, the vertices of the graph and the vertices in their respective adjacency lists are relabeled based on its position after sorting it in non-decreasing degree. To perform this relabeling efficiently, we use distributed counting sort. Note that a side effect of this reordering is on the utilization of the available locality. Hash-map based triangle counting routines could potentially make use of the locality obtained by processing the vertices in an order such that vertices with similar adjacency lists are processed together. However, although the degree-based ordering destroys this locality, the gains achieved by such a reordering results in faster runtimes and is proven in (Arifuzzaman et al., 2017).

Create the upper and the lower triangular portion of the adjacency matrix

Recall from Subsection 5.1 that the upper and the lower triangular portions of the adjacency matrix are used in order to perform the computations involved in the triangle counting phase. The algorithm first processes its local chunk of vertices and the associated adjacency lists to determine the vertices that will need to be distributed to remote nodes to form a 2D cyclic distribution. Once each node has received its portion of the adjacency matrix, a single scan through the adjacency list of each vertex in the chunk is performed to create the upper and the lower triangular portions. To convert the matrix into these triangular portions, the degree of a vertex is compared with the degree of an adjacent vertex, and the adjacent vertex is placed in the upper portion if its degree is greater than that of the former. If not, the vertex is placed in the lower portion of the adjacency matrix. Moreover, since the vertices are reordered based on non-decreasing degree, the global position of the vertices can be used to compare the degrees of the vertices. However, in many scenarios, the position of the adjacent vertex is not locally available. Thus, this requires us to perform a communication step with all nodes which adds to the overheads in the parallel algorithm.

5.4. Cost Analysis

In order to analyze our algorithm, we derive the parallel time complexity for the computation and the communication steps in the pre-processing and the triangle counting phases of our algorithm. Let $n$ be the total number of vertices of the graph, $m$ be the total number of edges of the graph, and $p$ the total number of ranks used. Moreover, let $d_{avg}$ and $d_{max}$ be the average degree and the maximum degree of the graph, respectively. The computations in the pre-processing phase, as discussed above, involve multiple scans over the chunk of the adjacency matrix owned by each rank which amounts to a computation time of $m/p$. Moreover, the communication step in the pre-processing phase requires an all-to-all personalized communication operation. Since we implemented this communication step using $p$ point-to-point send and receive operations, its complexity is lower bounded by $p+m/p$. The pre-processing phase also includes the distributed counting sort for relabeling the vertices in non-decreasing degree order. The computations associated with this sort includes two scans of the local vertices to determine the local maximum degree and the local positions, an overall reduction to compute the global maximum degree and, a scan with a cost of $d_{max}$ to determine the new labels of the vertices by computing the new positions in the distributed system. Furthermore, to find the new positions, we perform a prefix sum which incurs a communication time of $d_{max}\log(p)$. Thus, the total pre-processing phase takes time

[TABLE]

In the triangle counting phase, for each shift, the amount of computations is on the average $(n/\sqrt{p})\times(d^{2}_{avg}/p)$, whereas the communication cost is $(n/\sqrt{p})\times(d_{avg}/\sqrt{p})$. As a result, the overall amount of time, across the $\sqrt{p}$ shifts, spent in computation and communication for the triangle counting phase is

[TABLE]

6. Experimental methodology

6.1. Datasets

We used six large, real-world and synthetic graphs with varying degree distributions to evaluate the performance of our algorithms. Various statistics related to these graphs and the sources from where they were obtained are shown in Table 1. twitter and friendster datasets are sparse, social networks, g500-s26, g500-s27, g500-s28 and g500-s29 were generated using the graph500 generator provided in (Nelson et al., 2014). These follow the RMAT graph specifications (Graph500, 2018). Our algorithm creates these synthetic graphs as input to each run prior to calling our triangle counting routine. This way, we avoid reading the big graphs from the disk. We converted all the graph datasets to undirected, simple graphs.

6.2. Experimental setup

Our experiments were conducted on up to 29 dual-socket nodes with 24 cores of Intel Haswell E5-2680v3, each with 64 GB memory. More details about the system is detailed in (Institute, 2018). Our programs were developed using C and OpenMPI (v3.1.2), and compiled using GCC (v8.1.0) with -O3 optimization. We ran our MPI programs with the option --bind-to core. We ran our experiments with ranks ranging from 16 to 169, such that the number of ranks is a perfect square which forms a $\sqrt{p}\times\sqrt{p}$ processor grid. The number of nodes used for the different number of ranks is detailed in Table 2. In order to best utilize the available resources, we resort to minimizing the maximum number of nodes used such that the aggregate memory that exists in the selected number of nodes satisfies the memory requirement of our algorithm. Moreover, we have bound cores to socket as well. For example, with 36 ranks, we bind one core per socket to get better performance while ensuring our program does not run out of memory. We obtain runtimes starting with 16 ranks, as the largest graph in our experiments, i.e., g500-s29, required all the memory provided by 16 nodes for the processing.

Our competing approach, Havoq, was executed on 48 nodes, with all 24 cores being used across the node. The program was compiled using GCC (v5.4.0) and OpenMPI (v3.1.2). We had to resort to this particular version of GCC since Havoq ran successfully on this version. Moreover, Havoq produced faster runtimes on this version than on GCC (v8.1.0).

6.3. Performance metrics

We make use of the following two performance metrics for assessing the performance of our parallel formulation

1. (1)

Speedup - Speedup is computed by dividing the runtime obtained by the baseline algorithm against the runtime obtained by the parallel algorithm. We consider the runtime obtained with 16 ranks as the baseline and report the speedup of the algorithm computed against that. 2. (2)

Efficiency - Similar to above, we compute the efficiency obtained by using the runtime of the 16-rank case as the baseline. Specifically, if $p$ is the number of processors and $T_{p}$ is the parallel runtime of our algorithm (or one of its phases), then we use $16T_{16}/pT_{p}$.

7. Results

7.1. Parallel Performance

Table 2 shows the performance achieved by our algorithm on the two larger synthetic graphs and the two real-world graphs. From these results, we notice that as the number of ranks increase, the time required by the pre-processing and triangle counting times decrease. The overall speedup on 169 ranks relative to 16 ranks is in the range of $3.06$ to $6.93$ (compared to an expected speedup of $10.56$). Note that the synthetic graphs achieving better speedups over the real-world graphs is due to the fact that the synthetic graphs we experiment on are larger. Since the performance advantage of our algorithm tapers off after a certain number of ranks (which is also contingent on the size of the dataset), we restrict ourselves to measuring runtimes up until $13\times 13$ processor grid.

Comparing the scalability of the pre-processing and the triangle counting phases, we can see that the latter scales better. The relative speedup of the triangle counting phase on $169$ ranks is on the average $1.71$ times higher than that achieved by the pre-processing phase. This can also be seen in the efficiency plots of Figure 1, where the efficiency of the pre-processing phase decreases faster than the triangle counting phase when the number of ranks increase. In light of the analysis presented in Section 5.4, this scaling behavior was expected. The communication and the computation of the pre-processing phase is of the same order. However, the computation in the triangle counting phase is of the order $d_{avg}/\sqrt{p}$ more than the pre-processing phase. Thus, the triangle counting phase continues to scale better than the pre-processing phase with increasing number of ranks.

Moreover, for almost all graphs, the performance at $25$ ranks shows a super-linear speedup when compared to the runtimes obtained at $16$ ranks. We believe this happens because both the triangle counting and the pre-processing phase utilize caches better as the aggregate amount of cache memory increases with increasing number of ranks. This is further quantified in Figure 2, which plots the operation rate in kOps/second for both the phases for g500-s29. We can see that although the pre-processing phase continues to show higher operation rates with increasing ranks, the triangle counting phase shows its peak performance at $25$ ranks.

Finally, on analyzing the real-world graphs, we notice that twitter attains better speedups than friendster. We believe this happens because the triangle counting phase in twitter involves more work as opposed to friendster. We measure the average number of probes performed per shift in every rank for a $169$ rank run in both twitter and friendster, and we observe that the number of probes in twitter is $68$% more than that of friendster.

7.2. Sources of overhead

We analyze three different sources of parallel overheads in our triangle counting algorithm. The first has to do with load imbalance as a result of assigning a different number of computations to each rank during each one of the $\sqrt{p}$ steps of the algorithm. The second is due to the redundant work performed with increasing number of ranks. Finally, the third overhead we analyze is due to the time spent in communication as we increase the number of ranks.

Load imbalance

Recall from Section 5, that the computations during triangle counting is organized in $\sqrt{p}$ phases, and in each phase each processor performs a shift operation and processes the block of $U$ and $L$ that it just received. If the amount of computation associated with processing each pair of blocks is different, then this can lead to load imbalance. In order to quantify this load imbalance in the just the compute phase, we performed a series of experiments in which we measured the time obtained per shift for the computations involved in the triangle counting phase for 25 and 36 MPI ranks. The load imbalance was measured as the ratio of the maximum amount of time over all pairs of blocks over the average time. The results are shown in Table 3. For 25 MPI ranks, the load imbalance is $1.05$ and for 36 MPI ranks is $1.14$. We also quantify how distributing the data and the tasks contributes to the load imbalance. We count the number of non-zero tasks associated with each rank with increasing grid sizes, and compute the load imbalance. In general, the load imbalance that we observed was less than $6\%$, which can further explain the load imbalance observed over $25$ and $36$ ranks runtimes.

Redundant work

In Section 5, we discussed various optimizations that were designed to efficiently operate on the very sparse blocks of $U$ and $L$ (e.g., doubly sparse traversal), in order to eliminate redundant computations. However, those optimizations do not entirely eliminate the redundant computations. To measure how much extra work we do as we increase the number of ranks, we instrumented our code to count the number of tasks that result in the map-based set intersection operation throughout the execution of the triangle counting phase. The number of such tasks for g500-s29 on 16, 25, and 36 ranks are shown in Table 4. We see that as the number of ranks grow from 16 to 25 and 25 to 36, the number of tasks increases by $25\%$ and $20\%$, respectively. This extra work is responsible for some of the loss in the potential speedup observed in the results shown in Table 3.

Communication overheads

The fraction of time spent on communication over the entire runtime of the pre-processing and the triangle counting phase can be identified by looking at Figure 3. This plot shows that for both pre-processing and triangle counting, the bulk of the time is spent in computations for the largest graph in our testbed. However, the portion of the overall runtime that is attributed to the communication keeps increasing as we increase the number of ranks.

7.3. Quantifying the gains achieved by the optimizations

Recall from Section 5.2 that we introduced various optimizations in our triangle counting phase that leverage the sparsity that occur in the graphs and the structure of the computations. In order to quantify the reduction in the runtime of the triangle counting phase, we choose the first two optimizations, which we believe gave us maximum benefits: (i) using a doubly sparse traversal of the CSR structure, and, (ii) modifying the hashing routine for sparser vertices, and, recorded the runtime of the triangle counting phase without using these optimizations.

Based on the results that we obtained from these experiments on g500-s29, we observe that the doubly-sparse traversal of the vertices has reduced the runtime of the triangle counting phase by $10\%$ and $15\%$ for $16$ ranks and $100$ ranks, respectively. In a similar light, the modified hashing routine has reduced the runtime of the triangle counting phase by $1.2\%$ and $8.7\%$ for $16$ ranks and $100$ ranks, respectively. Moreover, we also recorded the improvement obtained by using the enumeration scheme as $\langle j,i,k\rangle$ as opposed to $\langle i,j,k\rangle$ in our algorithm. We observe that the triangle counting runtime decreased by $72.8\%$ when we use the $\langle j,i,k\rangle$ enumeration scheme as compared to the $\langle i,j,k\rangle$ enumeration scheme.

7.4. Comparison against other algorithms

Comparison against Havoq (Pearce, 2017)

As discussed in Section 4, various distributed memory parallel triangle counting algorithms have been developed. We perform two different evaluations. First, we perform a direct comparison with Havoq on the graphs detailed in Table 1. Table 5 compares the triangle counting runtime obtained by Havoq and the triangle counting time obtained by our approach. Havoq runtimes were obtained on 1152 cores (using 48 nodes) and our runtimes were obtained on 169 cores. On average, we get a speedup of $10.2$ times over their approach. In friendster, our approach is slower than Havoq. We believe this is because Havoq does an edge-based partitioning scheme (referred to as delegate partitioning), which leads to better scaling capability as compared to our method, which incurs more overheads as the number of ranks increase. Furthermore, Havoq required more number of nodes than what was available in our system for g500-s29 and we could not obtain the runtime for the same.

Comparison against other distributed-memory algorithm (Arifuzzaman

et al., 2017; Kanewala et al., 2018)

We also contrast the performance achieved by our algorithm against what was achieved by other previous approaches on only the twitter graph, since this was the common benchmark. For this second evaluation, we use the runtimes that were reported in the respective papers, which were obtained on different architectures and number of ranks. Thus, the comparisons presented with these approaches in Table 6 should be interpreted in view of this caveat. The performance achieved by the various algorithms on twitter is shown in Table 6. Algorithm with Overlapping Partitioning (AOP) (Arifuzzaman et al., 2017) was run on 200 cores. Their experimental setup included 64 computing nodes (QDR InfiniBand interconnect) with 16 processors (Sandy Bridge E5-2670, 2.6GHz) per node, has 4GB memory per processor, and uses the operating system CentOS Linux 6. Surrogate (Arifuzzaman et al., 2017) was also run on 200 cores and their experimental setup was the same as that of AOP. OPT-PSP algorithm (Kanewala et al., 2018) was run on 2048 cores. The experimental setup for this algorithm included a Cray XC system which has 2 Broadwell 22-core Intel Xeon processors, and the scaling experiments used only up to 16 cores to uniformly double the problem size. From these results we can see that our implementation is comparable to all previous approaches. Moreover, the relative performance advantage of our method still holds, when we account for the fact that some of the runtimes reported in Table 6 use more number of cores than those used in our experiments.

8. Conclusion

In this paper we presented a distributed memory formulation for triangle counting and evaluated its performance on real-word and synthetic graphs. Compared to prior parallel formulations, our formulation utilizes a 2D decomposition which increases the concurrency that it can exploit while reducing the overall communication overhead. The experimental results showed that these features lead to good scaling performance. Our analysis also identified areas that can benefit from further algorithmic and data structure improvements in order to better balance the work and reduce the amount of redundant computations. Moreover, we also note that this work can be easily extended to deal with rectangular processor grids using the SUMMA (Van De Geijn and Watts, 1997) algorithm.

Acknowledgements.

This work was supported in part by NSF (1447788, 1704074, 1757916, 1834251), Army Research Office (W911NF1810344), Intel Corp, and the Digital Technology Center at the University of Minnesota. Access to research and computing facilities was provided by the Digital Technology Center and the Minnesota Supercomputing Institute.