arXiv:1704.03767v1 [cs.DC] 12 Apr 2017

Parallelized Kendall’s Tau Coefficient Computation via
SIMD Vectorized Sorting On Many-Integrated-Core
Processors
Yongchao Liua,∗, Tony Pana , Oded Greena , Srinivas Alurua,∗
a School

of Computational Science & Engineering, Georgia Institute of Technology, Atlanta,
GA 30332, USA

Abstract
Pairwise association measure is an important operation in data analytics. Kendall’s
tau coefficient is one widely used correlation coefficient identifying non-linear relationships between ordinal variables. In this paper, we investigated a parallel
algorithm accelerating all-pairs Kendall’s tau coefficient computation via single
instruction multiple data (SIMD) vectorized sorting on Intel Xeon Phis by taking advantage of many processing cores and 512-bit SIMD vector instructions.
To facilitate workload balancing and overcome on-chip memory limitation, we
proposed a generic framework for symmetric all-pairs computation by building
provable bijective functions between job identifier and coordinate space. Performance evaluation demonstrated that our algorithm on one 5110P Phi achieves
two orders-of-magnitude speedups over 16-threaded MATLAB and three ordersof-magnitude speedups over sequential R, both running on high-end CPUs. Besides, our algorithm exhibited rather good distributed computing scalability
with respect to number of Phis. Source code and datasets are publicly available
at http://lightpcc.sourceforge.net.
Keywords: Pairwise correlation; Kendall’s tau coefficient; all-pairs
computation; many integrated core; Xeon Phi

1. Introduction
Identifying interesting pairwise association between variables is an important
operation in data analytics. In bioinformatics and computational biology, one
typical application is to mine gene co-expression relationship via gene expression
data, which can be realized by query-based gene expression database search [1]
∗ Corresponding

author
Email addresses: yliu@cc.gatech.edu (Yongchao Liu), tpan7@gatech.edu (Tony Pan),
ogreen@gatech.edu (Oded Green), aluru@cc.gatech.edu (Srinivas Aluru)
1 Preliminary work was presented in the 28th International Symposium on Computer Architecture and High Performance Computing, Los Angeles, USA, 2016

Preprint submitted to Journal of Parallel and Distributed Computing

April 13, 2017

or gene co-expression network analysis [2]. For gene expression database search,
it targets to select the subject genes in the database that are co-expressed
with the query gene. One approach is to first define some pairwise correlation/dependence measure over gene expression profiles across multiple samples
(gene expression profiles for short) and then rank query-subject gene pairs by
their scores. For gene co-expression networks, nodes usually correspond to genes
and edges represent significant gene interactions inferred from the association
of gene expression profiles. To construct a gene co-expression network, allpairs computation over gene expression profiles is frequently conducted based
on linear (e.g. [3] [4] [5]) or non-linear (e.g. [6] [7] [8]) co-expression measures.
A variety of correlation/dependence measures have been proposed in the literature and among them, Pearson’s product-moment correlation coefficient [9]
(or Pearson’s r correlation) is the most widely used correlation measure [10].
However, this correlation coefficient is only applicable to linear correlations. In
contrast, Spearman’s rank correlation coefficient [11] (or Spearman’s ρ coefficient) and Kendall’s rank correlation coefficient [12] (or Kendall’s τ coefficient)
are two commonly used measures for non-linear correlations [13]. Spearman’s
ρ coefficient is based on Pearson’s r coefficient but applies to ranked variables,
while Kendall’s τ coefficient tests the association between ordinal variables.
These two rank-based coefficients were shown to play complementary roles in
the cases when Pearson’s r is not effective [14]. Among other non-linear measures, mutual information [15] [16] [17], Euclidean distance correlation [18] [19],
Hilbert-Schmidt information criterion [20], and maximal information criterion
[21] are frequently used as well. In addition, some unified frameworks for pairwise dependence assessments were proposed in the literature (e.g. [22]).
Kendall’s τ coefficient (τ coefficient for short) measures the ordinal correlation between two vectors of ordinal variables. Given two ordinal vectors
u = {u1 , u2 , ..., un } and v = {v1 , v2 , ..., vn }, where variables ui and vi are both
ordinal (0 ≤ i < n), the τ coefficient computes the correlation by counting
the number of concordant pairs nc and the number of discordant pairs nd by
treating ui and vi as a joint ordinal variable (ui , vi ). For the τ coefficient, a
pair of observations (ui , vi ) and (uj , vj ), where i 6= j, is deemed as concordant
if ui > uj and vi > vj or ui < uj and vi < vj , and discordant if ui > uj and
vi < vj or ui < uj and vi > vj . Note that if ui = uj or vi = vj , this pair is
considered neither concordant nor discordant.
In our study, we will consider two categories of τ coefficient, namely Tau-a
(denoted as τA ) and Tau-b (denoted as τB ). τA does not take into account tied
elements in each vector and is defined as
τA =

nc − nd
n0

(1)

where n0 = n(n − 1)/2. If all elements in each vector are distinct, we have
nc + nd = n0 and can therefore re-write Equation (1) as
τA =

2nd
n0 − 2nd
= 1−
n0
n0
2

(2)

As opposed to τA , τB makes adjustments for ties and is computed as
nc − nd
τB = p
(3)
(n0 − n1 )(n0 − n2 )
P ′ ′
P ′ ′
′
′
where n1 =
i ui (ui − 1)/2 and n2 =
i vi (vi − 1)/2. ui (vi ) denotes the
cardinality of the i-th group of ties for vector u (v). Within a vector, each
distinct value defines one tie group and this value acts as the identifier of the
corresponding group. The cardinality of a tie group is equal to the number of
elements in the vector whose values are identical to the identifier of the tie group.
If all elements in either vector are distinct, τB will equal τA as n1 = n2 = 0.
This indicates τA is a special case of τB and an implementation of τB will cover
τA inherently. In addition, from the definitions of τA and τB , we can see that
the computation of the τ coefficient between u and v is commutable.
Besides the values of correlation coefficients, some applications need to calculate P -value statistics to infer statistical significance between variables. For this
purpose, one approach is permutation test [23]. However, a permutation test
may need a substantial number of pairwise τ coefficient computation [24] even
for moderately large n, thus resulting in prohibitively long times for sequential
execution. In the literature, parallelizing pairwise τ coefficient computation has
not yet been intensively explored. One recent work is from Wang et al. [25],
which accelerated the sequential τ coefficient computation in R [26] based on
Hadoop MapReduce [27] parallel programming model. In [25], the sequential
all-pairs τ coefficient implementation in R was shown extremely slow on largescale datasets. In our study, we further confirmed this observation through our
performance assessment (refer to section 5.2).
In this paper, we parallelized all-pairs τ coefficient computation on Intel
Xeon Phis based on Many-Integrated-Core (MIC) architecture, the first work
accelerating all-pairs τ coefficient computation on MIC processors to the best
of our knowledge. This work is a continuation from our previous parallelization of all-pairs Pearson’s r coefficient on Phi clusters [28] and further enriches
our LightPCC library (http://lightpcc.sourceforge.net) targeting parallel
pairwise association measures between variables in big data analytics. In this
work, we have investigated three variants, namely the naı̈ve variant, the generic
sorting-enabled (GSE) variant and the vectorized sorting-enabled (VSE) variant, built upon three pairwise τ coefficient kernels, i.e. the naı̈ve kernel, the GSE
kernel and the VSE kernel, respectively. Given two ordinal vectors u and v of n
elements each, the naı̈ve kernel enumerates all possible pairs of joint variables
(ui , vi ) (0 ≤ i < n) to obtain nc and nd , resulting in O(n2 ) time complexity. In
contrast, both the GSE and VSE kernels take sorting as the core and manage
to reduce the time complexity to O(n log n).
Given m vectors of n elements each, the overall time complexity would be
O(m2 n2 ) for the naı̈ve variant and O(m2 n log n) for the GSE and VSE variants. The VSE variant enhances the GSE one by exploiting 512-bit wide single
instruction multiple data (SIMD) vector instructions in MIC processors to implement fast SIMD vectorized pairwise merge of sorted subarrays. Furthermore,
3

to facilitate workload balancing and overcome on-chip memory limitation, we
investigated a generic framework for symmetric all-pairs computation by pioneering to build a provable, reversible and bijective relationship between job
identifier and coordinate space in a job matrix.
The performance of our algorithm was assessed using a collection of real
whole human genome gene expression datasets. Our experimental results demonstrates that the VSE variant performs best on both the multi-threaded CPU
and Phi systems, compared to the other two variants. We further compared
our algorithm with the all-pairs τ coefficient implementations in the widely
used MATLAB [29] and R [26], revealing that our algorithm on a single 5110P
Phi achieves up to 812 speedups over 16-threaded MATLAB and up to 1,166
speedups over sequential R, both of which were benchmarked on high-end CPUs.
In addition, our algorithm exhibited rather good distributed computing scalability with respect to number of Phis.
2. Intel Many-Integrated-Core (MIC) Architecture
Intel MIC architecture targets to combine many Intel processor cores onto
a single chip and has already led to the release of two generations of MIC processors. The first generation is code named as Knights Corner (KNC) and
the second generation code named as Knights Landing (KNL) [30]. KNC is
a PCI Express (PCIe) connected coprocessor that must be paired with Intel
Xeon CPUs. KNC is actually a shared-memory computer [31] with full cache
coherency over the entire chip and running a specialized Linux operating system over many cores. Each core adopts an in-order micro-architecture and has
four hardware threads offering four-way simultaneous multithreading. Besides
scalar processing, each core is capable of vectorized processing from a newly
designed vector processing unit (VPU) featuring 512-bit wide SIMD instruction
set architecture (ISA). For KNC, each core has only one VPU and this VPU is
shared by all active hardware threads running on the same core. Each 512-bit
vector register can be split to either 8 lanes with 64 bits each or 16 lanes with 32
bits each. Note that the VPU does not support legacy SIMD ISAs such as the
Streaming SIMD extensions (SSE) series. As for caches, each core has separate
L1 instruction and data caches of size 32 KB each, and a 512 KB L2 cache
interconnected via a bidirectional ring bus with the L2 caches of all other cores
to form a unified shared L2 cache over the chip. The cache line size is 64 bytes.
In addition, two usage models can be used to invoke KNC: offload model and
native model, where the former relies on compiler pragmas/directives to offload
highly-parallel parts of an application to KNC, while the latter treats KNC as
symmetric multiprocessing computers. While primarily focusing on KNC Phis
in this work, we note that our KNC-based implementations can be easily ported
onto KNL processors, as KNL implements a superset of KNC instruction sets.
We expect that our implementations will be portable to future Phis as well.

4

Table 1: A list of notions used

Notion
m
n
u
ui
nc
nd
n0
n1
n2
n3
τA
τB
SX
|SX |
vX

Description
number of vectors
number of elements per vector
vector u of n elements, likewise for v
i-th element of vector u, likewise for vi
number of concordant variable pairs
number of discordant variable pairs
n(n − 1)/2
P
′
′
′
Pi u′i (u′i − 1)/2, u′i is size of the i-th tie group in u
Pi vi (vi − 1)/2, vi is size of the i-th tie group in v
i wi (wi − 1)/2, wi is size of the i-th joint tie group for u and v
(nc − nd )/n
p0
(nc − nd )/ (n0 − n1 )(n0 − n2 )
a sorted subarray
length of the sorted subarray SX
a 512-bit SIMD vector with 16 32-bit integer lanes

3. Pairwise Correlation Coefficient Kernels
For the τ coefficient, we have investigated three pairwise τ coefficient kernels:
the naı̈ve kernel, the GSE kernel and the VSE kernel. From its definition,
it can be seen that the Kendall’s τ coefficient only depends on the order of
variable pairs. Hence, given two ordinal vectors, we can first order all elements
in each vector, then replace the original value of every element with its rank
in each vector, and finally conduct the τ coefficient computation on the rank
transformed new vectors. This rank transformation does not affect the resulting
coefficient value, but could streamline the computation, especially for ordinal
variables in complex forms of representation. Moreover, this transformation
needs to be done only once beforehand for each vector. Hence, we will assume
that all ordinal vectors have already been rank transformed in the following
discussions. For the convenience of discussion, Table 1 shows a list of notions
used across our study.
3.1. Naı̈ve Kernel
The naı̈ve kernel enumerates all possible combinations of joint variables
(ui , vi ) (0 ≤ i < n) and counts the number of concordant pairs nc as well
as the number of discordant pairs nd . As mentioned above, given two joint
variable pairs (ui , vi ) and (uj , vj ) (i 6= j), they are considered concordant if
ui > uj and vi > vj or ui < uj and vi < vj , discordant if ui > uj and vi < vj
or ui < uj and vi > vj , and neither concordant nor discordant if ui = uj or
vi = vj . Herein, we can observe that the two joint variables are concordant
if and only if the value of (ui − uj ) × (vi − vj ) is positive; discordant if and
only if the value of (ui − uj ) × (vi − vj ) is negative; and neither concordant nor
5

discordant if and only if the value of (ui − uj ) × (vi − vj ) is equal to zero. In this
case, in order to avoid branching in execution paths (particularly important for
processors without hardware branch prediction units), we compute the value of
nc − nd by examining the sign bit of the product of ui − uj and vi − vj (refer to
lines 2 and 11 in Algorithm 1).
Algorithm 1 Pseudocode of our naı̈ve kernel
1: function calc sign(v)
2:
return (v > 0) − (v < 0);
3: end function

⊲ return 1 if v > 0, -1 if v < 0 and 0, otherwise

4: function kendall tau a naı̈ve(u, v, n)
5:
norminator = 0;
6:
for i = 1; i < n; ++i do
7:
a = ui ; b = vi ;
8:
#pragma vector aligned
9:
#pragma simd reduction(+:nominator)
10:
for j = 0; j < i; ++j do
11:
nominator += calc sign((a − uj )) × (b − vj ));
12:
end for
13:
end for
;
14:
return nominator
n(n−1)/2
15: end function

⊲ norminator represents nc − nd

⊲ compute nc − nd

Algorithm 1 shows the pseudocode of the naı̈ve kernel. From the code, the
naı̈ve kernel has a quadratic time complexity in a function of n, but its runtime
is independent of the actual content of u and v, due to the use of function
calc sign. Meanwhile, the space complexity is O(1). Note that this naı̈ve
kernel is only used to compute τA .
3.2. Generic Sorting-enabled Kernel
Considering the close relationship between calculating τ and ordering a list of
variables, Knight [32] proposed a merge-sort-like divide-and-conquer approach
with O(n log n) time complexity, based on the assumption that no element tie
exists within any vector. As this assumption is not always the case, Knight did
mention this drawback and suggested an approximation method by averaging
counts, rather than propose an exact solution. In this subsection, we investigate
an exact sorting-enabled solution to address both cases: with or without element
ties within any vector, together in a unified manner.
Given two ordinal vectors u and v, this GSE kernel generally works in the
following five steps.
• Step1 sorts the list of joint variables (ui , vi ) (0 ≤ i < n) in ascending order,
where the joint variables are sorted first by the first element ui and secondarily by the second element vi . In this step, we used quicksort via the
standard qsort library routine, resulting in O(n log n) time complexity.
• Step2 performs a linear-time scan over the sorted list to compute n1 (refer
to Equation (3)) by counting the number of groups consisting of tied
values as well as the number of tied values in each group. Meanwhile, we
compute a new value n3 for joint ties, with respect to the pair (ui , vi ), as
6

P

i wi (wi − 1)/2, where wi represents the number of jointly tied values in
the i-th group of joint ties for u and v.

• Step3 counts the number of discordant pairs nd by re-sorting the sorted
list obtained in Step1 in ascending order of all elements in v via a merge
sort procedure that can additionally accumulate the number of discordant
joint variable pairs each time two adjacent sorted subarrays are merged.
The rationale is as follows. Firstly, when merging two adjacent sorted
subarrays, we count the number of discordant pairs by only performing
pairwise comparison between joint variables from distinct subarrays. In
this way, we can ensure that every pair of joint variables will be enumerated once and only once during the whole Step3 execution. Secondly,
given two adjacent sorted subarrays to merge, it is guaranteed that the
first value (corresponding to u) of every joint variable in the left subarray
(with the smaller indices) is absolutely less than or equal to the first value
of every joint variable in the right subarray (with the larger indices), due
to the sort conducted in Step1. In particular, when the first value is identical for the two joint variables from distinct subarrays, the second value
(corresponding to v) of the joint variable from the left subarray is also absolutely less than or equal to the second value of the joint variable from the
left subarray. This means that discordance occurs only if the second value
of a joint variable from the right subarray is less than the second value of
a joint variable from the left subarray. Therefore, the value of nd can be
gained by accumulating the number of occurrences of the aforementioned
discordance in every pairwise merge of subarrays. Algorithm 2 shows the
pseudocode of counting discordant pairs with out-of-place pairwise merge
of adjacent sorted subarrays, where the time complexity is O(n log n) and
the space complexity is O(n).
• Step4 performs a linear-time scan over the sorted list obtained in Step3
to compute n2 (see Equation (3)) in a similar way to Step2. This works
because the sorted list actually corresponds to a sorted list of all elements
in v.
• Step 5 computes the numerator nc − nd in Equations (1) and (3) as n0 −
n1 − n2 + n3 − 2nd . Note that if there is no tie in each vector, n1 , n2 and
n3 will all be zero. In this case, nc − nd will be equal to n0 − 2nd as shown
in Equation (2).
From the above workflow, we can see that the GSE kernel takes into account
tied elements within each vector. Unlike the naı̈ve kernel that computes τA ,
the GSE kernel targets the computation of τB . As mentioned above, τA is
actually a special case of τB and an algorithm for τB will inherently cover τA .
Therefore, our GSE kernel is able to calculate both τA and τB in a unified
manner. In addition, the GSE kernel has O(n log n) time complexity, since the
time complexity is O(n log n) for both Step1 and Step3, O(n) for both Step2
and Step4, and O(1) for Step 5.
7

Algorithm 2 Pseudocode of Step3 of our GSE kernel
1: function gse merge(in, out, lef t, mid, right)
2:
l = p = lef t; r = mid; nd = 0;
3:
while l < mid && r < right do
4:
if in[r].v < in[l].v then
5:
nd + = mid − l;
6:
out[p++] = in[r++];
7:
else
8:
out[p++] = in[l++];
9:
end if
10:
end while
11:
return nd ;
12: end function
13: function kendall tau b step3 gse(pairs, buf f er, n)
14:
nd = 0; in = pairs; out = buf f er;
15:
for s = 1; s < n; s *= 2 do
16:
for l = 0; l < n; l += 2 * s do
17:
m = min(l + s, n);
18:
r = min(l + 2 ∗ s, n);
19:
nd += gse merge(in, out, l, m, r);
20:
end for
21:
swap(in, out);
22:
end for
23:
if pairs != in then
24:
memcpy(pairs, in, n * sizeof(*pairs));
25:
end if
26:
return nd ;
27: end function

⊲ merge two sorted subarrays
⊲ count discordant pairs

⊲ Perform out-of-place merge
⊲ swap in and out

3.3. Vectorized Sorting-enabled Kernel
The VSE kernel enhances the GSE kernel by employing 512-bit SIMD vector
instructions on MIC processors to implement vectorized pairwise merge of sorted
subarrays. In contrast with the GSE kernel, the VSE kernel has made the
following algorithmic changes. The first is packing a rank variable pair (ui , vi )
into a signed 32-bit integer. In this packed format, each variable is represented
by 15 bits in the 32-bit integer, with ui taking the most significant 16 bits and
vi the least significant 16 bits. In this case, the VSE kernel limits the maximum
allowable vector size n to 215 − 1 = 32, 767. The second is that due to packing,
we replace the generic variable pair quicksort in Step1 with an integer sorting
method, and re-implement the discordant pair counting algorithm in Step3 (see
Algorithm 2) based on integer representation. In Step1, sorting packed integers
is equivalent to sorting generic variable pairs (ui , vi ), since within any packed
integer ui sits in higher bits and vi in lower bits. In Step3, an additional
preprocessing procedure is needed to reset to zero the most significant 16 bits
of each packed integer (corresponding to u). In our implementation, we split a
512-bit SIMD vector into a 16-lane 32-bit-integer vector and then investigated
a vectorized merge sort for Step1 and a vectorized discordant pair counting
algorithm for Step3, both of which use the same pairwise merge method and
also follow a very similar procedure. Algorithm 3 shows the pseudocode of Step3
of our VSE kernel.
In the literature, some research has been done to accelerate sorting algorithms by means of SIMD vectorized pairwise merge of sorted subarrays. Hiroshi et al. [33] employed a SSE vectorized odd-even merge network [34] to
merge sorted subarrays in an out-of-core way. Chhugan et al. [35] adopted

8

Algorithm 3 Pseudocode of Step3 our VSE kernel
1: function vse merge(in, out, lef t, mid, right)
2:
l = lef t; r = mid; p = lef t; nd = 0;
3:
if mid − lef t < 16||right − mid < 16 then
4:
while l < mid && r < right do
5:
if input[r] < input[l] then
6:
nd + = mid − l; out[p++] = in[r++];
7:
else
8:
out[p + +] = in[l + +];
9:
end if
10:
end while
11:
else
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:

⊲ correspond to variable v

⊲ count discordant pairs
while l < mid && r < right do
if in[r] < in[l] then
nd + = mid − 1; r + +;
else
l + +;
end if
end while

⊲ correspond to variable v

⊲ merge two sorted subarrays
l = lef t; r = mid;
vMin = mm512 load epi32(in + l); vMax = mm512 load epi32(in + r);
l+ = 16; r+ = 16;
while true do
if mm512 reduce min epi32(vMin) ≥ mm512 reduce max epi32(vMax) then
mm512 store epi32(out + p, vMax); p+ = 16; vMax = vMin;
else if mm512 reduce min epi32(vMax) ≥ mm512 reduce max epi32(vMin) then
mm512 store epi32(out + p, vMin); p+ = 16;
else
⊲ invoke Algorithm 4
bitonic merge 16way(vMin, vMax);
mm512 store epi32(out + p, vMin); p+ = 16;
end if
if l + 16 ≥ mid||r + 16 ≥ right then
break;
end if
A = in[l]; B = in[r]; C = mm512 reduce max epi32(vMax);
if C ≤ A && C ≤ B then
mm512 store epi32(out + p, vMax); p+ = 16;
vMin = mm512 load epi32(in + l); l+ = 16;
vMax = mm512 load epi32(in + r); r+ = 16;
else if B < A then
vMin = mm512 load epi32(in + r); r+ = 16;
else
vMin = mm512 load epi32(in + l); l+ = 16;
end if
end while
end if
⊲ invoke Algorithm 5
bitonic merge 16way leftover(in, out, l, r, p, mid, right, vMax);
return nd ;
end function

49: function kendall tau b step3 vse(pairs, buf f er, n)
50:
nd = 0; in = pairs; out = buf f er;
51:
for s = 1; s < n; s *= 2 do
52:
for l = 0; l < n; l += 2 * s do
53:
m = min(l + s, n);
54:
r = min(l + 2 ∗ s, n);
55:
nd += vse merge(in, out, l, m, r);
56:
end for
57:
swap(in, out);
58:
end for
59:
if pairs != in then
60:
memcpy(pairs, in, n * sizeof(*pairs));
61:
end if
62:
return nd ;
63: end function

9

⊲ Perform out-of-place merge
⊲ swap in and out

u15
u14
u13
u12
u11
u10
u9
u8
u7
u6
u5
u4
u3
u2
u1
u0

v0
v1
v2
v3
v4
v5
v6
v7
v8
v9
v10
v11
v12
v13
v14
v15

L
H
L
H
L
H
L
H
L
H
L
H
L
H
L
H

L
H
L
H
L
H

L
H
L
H
L
H
L
H
L
H

a0
a1
a2
a3
a4
a5
a6
a7
a8
a9
a10
a11
a12
a13
a14
a15

b0
b1
b2
b3
b4
b5
b6
b7
b8
b9
b10
b11
b12
b13
b14
b15

L
H
L
H
L
H

L
H
L
H
L
H
L
H
L
H
L
H
L
H

L
H
L
H
L
H
L
H
L
H
L
H

a0
a1
a2
a3
a4
a5
a6
a7
a8
a9
a10
a11
a12
a13
a14
a15
b0
b1
b2
b3
b4
b5
b6
b7
b8
b9
b10
b11
b12
b13
b14
b15

L
H
L
H
L
H
L
H
L
H
L
H
L
H
L
H

L
H
L
H
L
H

L
H
L
H
L
H
L
H
L
H

a0
a1
a2
a3
a4
a5
a6
a7
a8
a9
a10
a11
a12
a13
a14
a15
b0
b1
b2
b3
b4
b5
b6
b7
b8
b9
b10
b11
b12
b13
b14
b15

L
H
L
H
L
H
L
H
L
H
L
H
L
H
L
H
L
H
L
H
L
H
L
H
L
H
L
H
L
H
L
H

a0
a1
a2
a3
a4
a5
a6
a7
a8
a9
a10
a11
a12
a13
a14
a15
b0
b1
b2
b3
b4
b5
b6
b7
b8
b9
b10
b11
b12
b13
b14
b15

L
H
L
H
L
H
L
H
L
H
L
H
L
H
L
H
L
H
L
H

L
H
L
H
L
H
L
H
L
H
L
H

a0
a1
a2
a3
a4
a5
a6
a7
a8
a9
a10
a11
a12
a13
a14
a15
b0
b1
b2
b3
b4
b5
b6
b7
b8
b9
b10
b11
b12
b13
b14
b15

Figure 1: 16-way bitonic merge network: (1) input data are stored in two 16-lane vectors and
(2) pipelining from left-to-right

similar ideas to [33], but combined a SSE vectorized in-register odd-even merge
sort [34] with an in-memory bitonic merge network. Chen et al. [36] absorbed
the merge-path idea of [37, 38] and extended the SSE vectorized work of [35]
to take advantage of 512-bit SIMD VPUs on KNC Phis and used a vectorized
in-register bitonic merge sort, instead of the in-register odd-even merge sort.
In our VSE kernel, we engineered a 512-bit SIMD vectorized in-register bitonic
merge network as the core of our pairwise merge procedure for sorted subarrays,
which is similar to [36], and further proposed a predict-and-skip mechanism to
reduce the number of comparisons during pairwise merge.
3.3.1. In-register bitonic merge network
The in-register bitonic merge network is the core of our vectorized pairwise
merge of sorted subarrays. In our algorithm, this network has 16 ways and
merges two sorted vectors vM in and vM ax (Figure 1 shows the computation
layout of the network), where all elements in each vector are placed in ascending
order from lane 0 to lane 15. In this case, to generate an input bitonic sequence
from the two vectors, we need to reverse the order of all elements in one and only
one vector (reverse vM ax in our case) and this order reversal is realized by one
permutation instruction, i.e. mm512 permutevar epi32(·). Having completed
the order reversal, we can complete the sorting of vectors vM in and vM ax in
log2 (32) = 5 steps. Algorithm 4 shows the pseudocode for our 16-way in-register
bitonic merge network. From the code, the function bitonic merge 16way is

10

composed of a fixed number of vector instructions and thus has a constant time
complexity.
Algorithm 4 Pseudocode of our 16-way in-register bitonic merge network
1: procedure bitonic merge 16way(vMin, vMax)
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:

⊲ constant vector variables and reused
//vReverse = mm512 set epi32(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15);
//vP ermIndex16 = mm512 set epi32(7, 6, 5, 4, 3, 2, 1, 0, 15, 14, 13, 12, 11, 10, 9, 8)
//vP ermIndex8 = mm512 set epi32(11, 10, 9, 8, 15, 14, 13, 12, 3, 2, 1, 0, 7, 6, 5, 4);
//vP ermIndex4 = mm512 set epi32(13, 12, 15, 14, 9, 8, 11, 10, 5, 4, 7, 6, 1, 0, 3, 2);
//vP ermIndex2 = mm512 set epi32(14, 15, 12, 13, 10, 11, 8, 9, 6, 7, 4, 5, 2, 3, 0, 1);
⊲ Reserve vector vMax
vMax = mm512 permutevar epi32(vReverse, vMax);
⊲ Level 1
vL1 = mm512 min epi32(vMin, vMax);
vH1 = mm512 max epi32(vMin, vMax);
⊲ Level 2
vT mp = mm512 permutevar epi32(vP ermIndex16, vL1);
vT mp2 = mm512 permutevar epi32(vP ermIndex16, vH1);
vL2 = mm512 mask min epi32(vL2, 0x00ff, vT mp, vL1);
vH2 = mm512 mask min epi32(vH2, 0x00ff, vT mp2, vH1);
vL2 = mm512 mask max epi32(vL2, 0xff00, vT mp, vL1);
vH2 = mm512 mask max epi32(vH2, 0xff00, vT mp2, vH1);
⊲ Level 3
vT mp = mm512 permutevar epi32(vP ermIndex8, vL2);
vT mp2 = mm512 permutevar epi32(vP ermIndex8, vH2);
vL3 = mm512 mask min epi32(vL3, 0x0f0f, vT mp, vL2);
vH3 = mm512 mask min epi32(vH3, 0x0f0f, vT mp2, vH2);
vL3 = mm512 mask max epi32(vL3, 0xf0f0, vT mp, vL2);
vH3 = mm512 mask max epi32(vH3, 0xf0f0, vT mp2, vH2);
⊲ Level 4
vT mp = mm512 permutevar epi32(vP ermIndex4, vL3);
vT mp2 = mm512 permutevar epi32(vP ermIndex4, vH3);
vL4 = mm512 mask min epi32(vL4, 0x3333, vT mp, vL3);
vH4 = mm512 mask min epi32(vH4, 0x3333, vT mp2, vH3);
vL4 = mm512 mask max epi32(vL4, 0xcccc, vT mp, vL3);
vH4 = mm512 mask max epi32(vH4, 0xcccc, vT mp2, vH3);
⊲ Level 5: vMin and vMax store the sorted sequence
vT mp = mm512 permutevar epi32(vP ermIndex2, vL4);
vT mp2 = mm512 permutevar epi32(vP ermIndex2, vH4);
vMin = mm512 mask min epi32(vMin, 0x5555, vT mp, vL4);
vMax = mm512 mask min epi32(vMax, 0x5555, vT mp2, vH4);
vMin = mm512 mask max epi32(vMin, 0xaaaa, vT mp, vL4);
vMax = mm512 mask max epi32(vMax, 0xaaaa, vT mp2, vH4);
end procedure

3.3.2. Vectorized pairwise merge of sorted subarrays
Our vectorized pairwise merge of sorted subarrays relies on the aforementioned 16-way in-register bitonic merge network and adopted a very similar
procedure to [33]. Given two sorted subarrays SA and SB , we assume that they
are aligned to 64 bytes and their lengths |SA | and |SB | are multiples of 16, for
the convenience of discussion. In this way, the vectorized pairwise merge works
as follows.
1. loads the smallest 16 elements of SA and SB to vM in and vM ax, respectively, and advances the pointer of each subarray.
2. invokes bitonic merge 16way(·) to sort vectors vM in and vM ax, and
then stores the content of vM in, the smallest 16 elements, to the resulting
output array.
11

Algorithm 5 Pseudocode of processing the leftovers in both subarrays
1: procedure bitonic merge 16way leftover(in, out, l, r, p, mid, right, vMax)
2:
vIndexInc = mm512 set epi32(15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0);
⊲ maximum signed integer value
3:
vMaxInt = mm512 set1 epi32(0x7fffffff)
4:
vMid = mm512 set1 epi32(mid);
5:
vRight = mm512 set1 epi32(right);
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:

⊲ initialize the vector mask for vMax
maskMax = 0xffff;
while l < mid||r < right do
if l < mid && r < right then
if in[r] < in[l] then
vT mp = mm512 add epi32( mm512 set1 epi32(r), vIndexInc)
maskMin = mm512 cmplt epi32 mask(vT mp, vRight);
vMin = mm512 mask load epi32(vMaxInt, maskMin, in + r); r+ = 16;
else
vT mp = mm512 add epi32( mm512 set1 epi32(l), vIndexInc)
maskMin = mm512 cmplt epi32 mask(vT mp, vMid);
vMin = mm512 mask load epi32(vMaxInt, maskMin, in + l); l+ = 16;
end if
else if l < mid then
vT mp = mm512 add epi32( mm512 set1 epi32(l), vIndexInc);
maskMin = mm512 cmplt epi32 mask(vT mp, vMid);
vMin = mm512 mask load epi32(vMaxInt, maskMin, in + l);
tmp = mm512 mask reduce max epi32(maskMax, vMax);
if tmp ≤ mm512 reduce min epi32(vMin) then
break;
end if
l+ = 16;
else if r < right then
vT mp = mm512 add epi32( mm512 set1 epi32(r), vIndexInc);
maskMin = mm512 cmplt epi32 mask(vT mp, vRight);
vMin = mm512 mask load epi32(vMaxInt, maskMin, in + r);
tmp = mm512 mask reduce max epi32(maskMax, vMax);
if tmp ≤ mm512 reduce min epi32(vMin) then
break;
end if
r+ = 16;
end if
⊲ invoke Algorithm 4
bitonic merge 16way(vMin, vMax);
⊲ calculate new vector masks
nb = mm countbits 32(maskMin) + mm countbits 32(maskMax);
maskMin = nb > 15? 0x0ffff : (1 << nb) − 1;
maskMax = nb < 17? 0 : (1 << (nb − 16)) − 1;
mm512 mask store epi32(out + p, maskMin, vMin);
p += mm countbits 32(maskMin);
end while
⊲ write out vMax
mm512 mask packstorelo epi32(out + p, maskMax, vMax);
mm512 mask packstorehi epi32(out + p + 16, maskMax, vMax);
p += mm countbits 32(maskMax);
⊲ copy out the rest
if l < mid then
memcpy(out + p, in + l, (mid − l)∗ sizeof(int));
else if r < right then
memcpy(out + p, in + r, (right − r)∗ sizeof(int));
end if
end procedure

12

3. compares the next element of SA and SB and loads 16 elements into vM in
from the subarray whose next element is the smaller, followed by pointer
advancing for the subarray.
4. jumps back to the second step and repeats the procedure until all elements
in both subarrays have been processed.
5. completes the merge of SA and SB by writing the content of vM ax to the
output array.
Note that if |SA | or |SB | is not a multiple of 16, we will have to implement
some special treatment to deal with the leftovers in SA and SB . This special
treatment is implemented as Algorithm 5 and invoked by Algorithm 3 (see line
46). From the above procedure, we can see that the time complexity of merging
SA and SB is linear and thus the VSE kernel has a time complexity of O(n log n).
As assumed that both SA and SB are aligned to 64 bytes, we can therefore
use only one mm512 load epi32(·) instruction to load 16 integer elements from
memory to a vector and only one mm512 store epi32(·) instruction to store 16
elements in a vector to memory. This is because memory address is guaranteed
to keep aligned to 64 bytes during pairwise merge. Defining Trd to denote
the latency of memory load (in clock cycles) and Twt to denote the latency of
memory store, the number of clock cycles per element load from memory can be
estimated as (Trd + 16)/16 = Trd/16 + 1 clock cycles and the number of clocks
cycles per element store to memory as (Twt + 16)/16 = Twt /16 + 1 clock cycles.
Moreover, bitonic merge 16way(·) is composed of 27 vector instructions (see
lines 7∼33 in Algorithm 4) and each of its invocations leads to the output of 16
elements. In these regards, the average number of instructions per element can
be estimated as 27/16 ≈ 1.69. Considering that we only used simple integer
and single-source permutation instructions, each of which has a two clock cycle
latency according to the Intel Xeon Phi processor software optimization manual
[39], the average computational cost per element can be estimated as ⌈2 ×
27/16⌉ = 4 clock cycles in the worst cases where consecutive instructions always
have data dependency. In order to reduce inter-instruction latency caused by
data dependency, we have manually tuned the order of instructions by means
of interleaving. Since each instruction has a latency of two clock cycles, this
interleaving could enable to execute one instruction in only one clock cycle,
thus further reducing the cost per element to ⌈1 × 27/16⌉ = 2 clock cycles
optimistically. In sum, the overall cost per element can be favorably estimated
to be (Trd + Twr )/16 + 4 clock cycles.
We have proposed a predict-and-skip scheme to determine (i) whether we
need to invoke the 16-way bitonic merge network and (ii) whether we should
load two vectors of new elements from SA and SB , respectively. For cases (i),
the rationale is if the minimum value in vM in (vM ax) is ≥ the maximum value
in vM ax (vM in), every value in vM in (vM ax) will be ≥ to all values in vM in
(vM ax). In this case, there is no need of invoking the bitonic merge network
(lines 23∼26 in Algorithm 3), since the order of all elements in both vectors
is already deterministic. For cases (ii), we compare the maximum value C in
vM ax with both the smallest element A in the rest of SA and the one B in
13

the rest of SB (lines 34∼38 in Algorithm 3). If C ≤ A and C ≤ B, none
of the elements in vM ax will be greater than any element in the rest of SA
and SB , indicating that the absolute order of all vM ax elements has already
been determined. In this case, we write vM ax to the resulting output array
and load two vectors of new elements from SA and SB to vM in and vM ax,
respectively. The updates in vM in and vM ax will be used for the next iteration
of comparison.
It is worth mentioning that we also developed a 32-way bitonic merge network to implement the VSE kernel. Unfortunately, this new kernel based on
32-way merge network yielded slightly inferior performance to the kernel with
the 16-way network through our tests. In this regard, we adopted the 16-way
merge network in our VSE kernel.
4. Parallel All-Pairs Computation
4.1. All-Pairs Computation Framework
We consider the m × m job matrix to be a 2-dimensional coordinate space
on the Cartesian plane, and define the left-top corner to be the origin, the
horizontal x-axis (corresponding to columns) in left-to-right direction and the
vertical y-axis (corresponding to rows) in top-to-bottom direction. For nonsymmetric all-pairs computation (non-commutative pairwise computation), the
workload distribution over processing elements (PEs), e.g. threads, processes,
cores and etc., would be relatively straightforward. This is because coordinates
in the 2-dimensional matrix corresponds to distinct jobs. Unlike non-symmetric
all-pairs computation, it suffices by only computing the upper-triangle (or lowertriangle) of the job matrix for symmetric all-pairs computation (commutative
pairwise computation). In this case, balanced workload distribution could be
more complex than non-symmetric all-pairs computation.
In this paper, we propose a generic framework for workload balancing in
symmetric all-pairs computation, since the computation of the τ coefficient between any u and v is commutable as mentioned above. This framework works
by assigning each job a unique identifier and then building a bijective relationship between a job identifier Jm (y, x) and its corresponding coordinate (y, x).
We refer to this as direct bijective mapping. While facilitating balanced
workload distribution, this mapping merely relies on bijective functions, which
is a prominent feature distinguished from existing methods. To the best of our
knowledge, in the literature bijective functions have not ever been proposed for
workload balancing in symmetric all-pairs computation. In [40], the authors
used a very similar job numbering approach to ours (shown in this study), but
did not derive a bijective function for symmetric all-pairs computation. Our
framework can be applied to cases with identical (e.g. using static workload
distribution) or varied workload per job (e.g. using a shared integer counter
to realize dynamic workload distribution via remote memory access operations
in MPI [41] and Unified Parallel C (UPC) programming models [42] [43]) and
is also particularly useful for parallel computing architectures with hardware

14

*"& (

"

0

1

2

3

8

9

10 11 12 13 14

4

5

6

7

! " # 0.5

15 16 17 18 19 20
21 22 23 24 25
26 27 28 29

30 31 32
!

)
!

! " # 0.5 #
"$ % " % 0.25 # 2&' % 1(
)

! " # 0.5 %
"$ % " % 0.25 # 2&' % 1(

33 34
35

(a)

(b)

(c)

Figure 2: (a) an example direct bijective mapping between job identifier and coordinate space,
(b) solution region for Equation (7), and (c) solution region for Equation (8)

schedulers such as GPUs and FPGAs. In the following, without loss of generality, we will interpret our framework relative to the upper triangle of the job
matrix by counting in the major diagonal. Nonetheless, this framework can be
easily adapted to the cases excluding the major diagonal.
Direct bijective mapping. Given a job (y, x) in the upper triangle, we compute
its integer job identifier Jm (y, x) as
Jm (y, x) = Fm (y) + x − y, 0 ≤ y ≤ x < m

(4)

for dimension size m. In this equation, Fm (y) is the total number of cells
preceding row y in the upper triangle and is computed as
Fm (y) =

y(2m − y + 1)
2

(5)

where y varies in [0, m]. In this way, we have defined Equation (4) based on our
job numbering policy, i.e. all job identifiers vary in [0, m(m + 1)/2) and jobs
are sequentially numbered left-to-right and top-to-bottom in the upper triangle
(see Fig. 2a for an example).
Reversely, given a job identifier J = Jm (y, x) (0 ≤ J < m(m+1)/2), we need
to compute the coordinate (y, x) in order to locate the corresponding variable
pair. As per our definition, we have

J ≥ Fm (y) ⇔ y 2 − (2m + 1)y + 2J ≥ 0
(6)
J ≤ Fm (y + 1) − 1 ⇔ y 2 − (2m − 1)y + 2(J + 1) − 2m ≤ 0
By solving J ≥ Fm (y) (see Figure 2b), we get
p
y ≤ m + 0.5 − m2 + m + 0.25 − 2J

(7)

while getting
y ≥ m − 0.5 −

p
m2 + m + 0.25 − 2(J + 1)
15

(8)

by
Figure 2c). In this case, by defining ∆ =
p solving J ≤ Fm (y + 1) − 1 (see
√
2 + m + 0.25 − 2(J + 1), ∆′ =
m
m2 + m + 0.25 − 2J and z = m − 0.5 −
p
m2 + m + 0.25 − 2(J + 1), we can reformulate Equations (7) and (8) to be
z ≤ y ≤ z + 1 + ∆ − ∆′ . Because ∆ < ∆′ , we know that [z, z + 1 + ∆ − ∆′ ]
is a sub-range of [z, z + 1) and thereby have z ≤ y < z + 1. Based on our job
numbering policy stated before, as a function of integer y, Equation (6) definitely
has y solutions, meaning that at least one integer exists in [z, z + 1 + ∆ − ∆′ ],
which satisfies Equation (6). Meanwhile, it is known that there always exists
one and only one integer in [z, z + 1) (can be easily proved) and this integer
equals ⌈z⌉, regardless of the value of z. Since [z, z + 1 + ∆ − ∆′ ] is a sub-range
of [z, z + 1), we can conclude that Equation (6) has a unique solution y that is
computed as
l
m
p
y = ⌈z⌉ = m − 0.5 − m2 + m + 0.25 − 2(J + 1)
(9)
Having got y, we can compute the coordinate x as
x = J + y − Fm (y)

(10)

based on Equation (4).
4.2. Tiled Computing
Based on the direct bijective mapping, we adopted tiled computing with
the intention to benefit from L1/L2 caches. The rationale behind the tiled
computing is loading into cache a small subset of the bigger dataset and reusing
this block of data in cache for multiple passes. This technique partitions a
matrix into a non-overlapping set of equal-sized q × q tiles. In our case, we
partition the job matrix and produce a tile matrix of size w × w tiles, where
w = ⌈m/q⌉. In this way, all jobs in the upper triangle of the job matrix are still
fully covered by the upper triangle of the tile matrix. By treating a tile as a
unit, we can assign a unique identifier to each tile in the upper triangle of the
tile matrix and then build bijective functions between tile identifiers and tile
coordinates in the tile matrix, similarly as we do for the job matrix.
Because the tile matrix has an identical structure to the original job matrix,
we can directly apply our bijective mapping to the tile matrix. In this case,
given a coordinate (yq , xq ) (0 ≤ yq ≤ xq < w) in the upper triangle of the tile
matrix, we can compute a unique tile identifier Jw (yq , xq ) as
Jw (yq , xq ) = Fw (yq ) + xq − yq , 0 ≤ yq ≤ xq < w

(11)

where Fw (yq ) is defined similar to Equation (5) as
Fw (yq ) =

yq (2w − yq + 1)
2

(12)

Likewise, given a tile identifier Jw , such that 0 ≤ Jw < w(w + 1)/2, we can
16

reversely compute its unique vertical coordinate yq as
l
m
p
yq = w − 0.5 − w2 + w + 0.25 − 2(Jw + 1)

(13)

and subsequently its unique horizontal coordinate xq as
xq = Jw + yq − Fw (yq )

(14)

As title size is subject to cache sizes and input data, tuning tile size is
believed to be important for gaining high performance. This tuning process,
however, is tedious and has to be conducted case-by-case. For convenience, we
empirically set q to 8 for CPUs and to 4 for Phis. Nevertheless, users can feel
free to tune tile sizes to meet their needs.
4.3. Multithreading
4.3.1. Asynchronous kernel execution
When m is large, we may not have sufficient memory to reside the resulting
m × m correlation matrix Mτ entirely in memory. To overcome memory limitation, we adopted a multi-pass kernel execution model which partitions the tile
identifier range [0, w(w + 1)/2) into a set of non-overlapping sub-ranges (equalsized in our case) and finishes the computation one sub-range after another. In
this way, we do not need to allocate the whole memory for matrix Mτ . Instead,
we only need to allocate a small amount of memory to store the computing
results of one sub-range, thus considerably reducing memory footprint.
When using the KNC Phi, we need to transfer the newly computed results
from the Phi to the host after having completed each pass of kernel execution.
If kernel execution and data transfer is conducted in sequential, the Phi will
be kept idle during the interim of device-to-host data transfer. In this regard,
a more preferable solution would be to employ asynchronous kernel execution
by enabling concurrent execution of host-side tasks and accelerator-side kernel
execution. Fortunately, KNC Phis enable such a kind of asynchronous data
transfer and kernel execution associated with the offload model. This asynchronous execution can be realized by coupling the signal and wait clauses,
both of which are associated with each other via a unique identifier. That is,
a signal clause initiates an asynchronous operation such as data transfer and
kernel execution, while a wait clause blocks the current execution until the associated asynchronous operation is completed. Refer to our previous work [28]
for more details about the host-side asynchronous execution workflow proposed.
4.3.2. Workload balancing
For the variant using the naı̈ve kernel, all jobs have the same amount of
computation (refer to Algorithm 1). Thus, given a fixed number of threads, we
evenly distribute jobs over the set of threads by assigning a thread to process
one tile at a time. In contrast, for the two variants using the sorting-enabled
kernels, different jobs may have different amount of computation because of

17

the two rounds of sort used in Step1 and Step3. Typically, an OpenMP dynamic scheduling policy is supposed to be in favor to address workload irregularity, albeit having relatively heavier workload distribution overhead than static
scheduling. However, it is observed that dynamic scheduling produces slightly
inferior performance to static scheduling through our tests. In this regard, we
adopted static scheduling in our two sorting-enabled variants and assigned one
thread to process one tile at a time, same as the naı̈ve variant.
4.4. Distributed Computing
On KNC Phi clusters, we used MPI offload model which launches MPI
processes just as an ordinary CPU cluster does. In this model, one or more
Phis will be associated to a parental MPI process, which utilizes offload pragmas/directives to interact with the affiliated Phis, and inter-Phi communications
need to be explicitly managed by parental processes. In this sense, Phis may not
perceive the existence of remote inter-process communications. Our distributed
implementations require one-to-one correspondence between MPI processes and
Phis, and adopted a static workload distribution scheme based on tiled computation. This static distribution is inspired by our practice in multithreading
on single Phis. In our implementation, given p processes, we evenly distribute
tiles onto the p processes with the i-th (0 ≤ i < p) process assigned to compute
the tiles whose identifiers are in [i × ⌈ w(w+1)
⌉, (i + 1) × ⌈ w(w+1)
⌉). Within each
2p
2p
process, we adopted asynchronous execution workflow as well.
5. Performance Evaluation
We assessed the performance of our algorithm from four aspects: (i) comparison between our three variants, (ii) comparison with widely used counterparts:
MATLAB (version R2015b) and R (version 3.2.0), (iii) multithreading scalability on a single Phi, and (iv) distributed computing scalability on Phi clusters.
In these tests, we used four real whole human genome gene expression datasets
(refer to Table 2) produced by Affymetrix Human Genome U133 Plus 2.0 Array.
These datasets are publicly available in the GPL570 data collection of SEEK
[1], a query-based computational gene co-expression search engine over large
transcriptomic databases. In this study, unless otherwise specified, we compute
τ coefficients between genes, meaning that m is equal to the number of genes
and n equal to the number of samples for each dataset.
All tests are conducted on 8 compute nodes in CyEnce HPC Cluster (Iowa
State University), where each node has two high-end Intel E5-2650 8-core 2.0
GHz CPUs, two 5110P Phis (each has 60 cores and 8 GB memory) and 128 GB
memory. Our algorithm is compiled by Intel C++ compiler v15.0.1 with option
-fast enabled. In addition, for distributed computing scalability assessment,
when two processes are launched into the same node, we used the environment
variable I MPI PIN PROCESSOR LIST to guide Intel MPI runtime system to pin
the two processes within the node to distinct CPUs (recall that one node has
two CPUs).
18

Table 2: Information of whole human genome gene expression datasets used

Name
DS21050
DS19784
DS13070
DS3526

No. of Genes
17,941
17,941
17,941
17,941

No. of Samples
310
320
324
353

Platform
Affymetrix
Affymetrix
Affymetrix
Affymetrix

5.1. Assessment of Our Three Variants
All of the three variants work on CPUs, Phis and their clusters. For the
naı̈ve and GSE variants, they both used the same C++ core code for CPUand MIC-oriented instances. For the VSE variant, its 512-bit SIMD vectorization is only applicable to Phis. In this regard, to support CPUs, we further
developed a non-vectorized merge sort in Step1 and a non-vectorized discordant pair counting algorithm in Step3 instead, based on the aforementioned
packed integer representation. Considering that this non-vectorized version also
works on Phis, we have examined how much our 512-bit SIMD vectorization
contributes to speed improvement by comparing the vectorized version with the
non-vectorized one on the same 5110P Phi. Interestingly, by measuring the
correlation between genes using the datasets in Table 2, we observed that the
non-vectorized version is on a par with the vectorized one. This may be due to
the relatively small value of n, which is only 353 at maximum. Based on this
consideration, we further evaluated both versions by measuring the correlation
between samples, where m will be equal to the number of samples and n equal to
the number of genes. In this case, n becomes 17,941 for each dataset, more than
50× larger than before. In this context, our performance assessment exposed
that the vectorized version is consistently superior to the non-vectorized one,
yielding very stable speedups averaged to be 1.87 for all datasets. This result
was exciting, proving that our 512-bit SIMD vectorization did boost speed even
for moderately large n values. Hence, we have used the vectorized version of
the VSE variant for performance measurement all throughout our study.
5.1.1. On CPU
We first compared the performance of our three variants on multiple CPU
cores. Table 3 shows the performance of each variant on 16 CPU cores and
Figure 3 shows the speedup of the 16-threaded instance of each variant over its
single-threaded one.
For the naı̈ve variant, its 16-threaded instance achieves a roughly constant
speedup of 13.70 over its single-threaded one. This observation is consistent
with our expectation, as the runtime of the naı̈ve kernel is subject to the vector
size n but independent of actual vector content (refer to the implementation
shown in Algorithm 1). In contrast, the speed of the two sorting-enabled kernels
are sensitive to vector content to some degree. This can be explained by the
following three factors: variable pair sorting in Step1, discordant pair counting
based on pairwise merge of sorted subarrays in Step3, and linear-time scans
19

14.0

Naïve

GSE

VSE

13.9

Speedup

13.8
13.7
13.6

13.5
13.4
DS21050

DS19784

DS13070

DS3526

Figure 3: Speedups of 16-threaded instances over single-thread ones on multiple CPU cores

for determining tie groups in Step2 and Step4. Nevertheless, for both sortingenabled variants, their 16-threaded instances yield relatively consistent speedups
over their single-threaded ones, respectively. For the GSE (VSE) variant, its 16threaded instance runs 13.74× (13.79×) faster than its single-threaded one on
average, with a minimum speedup 13.63 (13.69) and a maximum speedup 13.90
(13.97). For each case, the two sorting-enabled variants both yield superior
performance to the naı̈ve variant, where the average speedup is 1.60 for the
GSE kernel and 2.72 for the VSE kernel.
Table 3: Performance comparison of three variants on multiple CPU cores

Instance

Single-threaded

16-threaded

Dataset
DS21050
DS19784
DS13070
DS3526
DS21050
DS19784
DS13070
DS3526

Naı̈ve
10,052
10,690
10,975
12,993
734
781
801
948

20

Time (s)
GSE
6,609
6,776
6,920
7,619
475
495
504
559

VSE
3,758
3,938
4,041
4,742
275
286
294
340

Speedup
GSE VSE
1.52 2.67
1.58 2.71
1.59 2.72
1.71 2.74
1.54 2.67
1.58 2.73
1.59 2.72
1.70 2.79

Table 4: Performance comparison of the three variants on Phi

Dataset
DS21050
DS19784
DS13070
DS3526

Naı̈ve
400
424
433
506

Time (s)
GSE
419
434
459
488

VSE
261
266
276
296

Speedup
GSE VSE
0.95 1.53
0.98 1.59
0.94 1.57
1.04 1.71

5.1.2. On Xeon Phi
Subsequently, we evaluated the three variants on a single 5110P Phi (see
Table 4). From Table 4, the VSE variant outperforms the naı̈ve variant for each
dataset by yielding an average speedup of 1.60 with the minimum speedup 1.53
and the maximum speedup 1.71. The GSE variant is superior to the naı̈ve one
on the DS3526 dataset, but inferior to the latter on the remaining three datasets.
Interestingly, by revisiting the CPU results shown in Table 3, we recalled that
the former actually performs consistently better than the latter for each dataset
on CPU. Since both the naı̈ve and the GSE variants use the same pairwise τ
coefficient kernel code for their corresponding CPU- and MIC-oriented implementations, this discordant performance ranking between the two variants could
owe to the architectural differences of the two types of processing unit (PU).
For instance, by enabling auto-vectorization, the naı̈ve kernel (see Algorithm
1) can concurrently process 16 integer elements by one 512-bit SIMD vector
instruction, in contrast with 4 integer elements by one 128-bit SSE instruction
on CPU. This fourfold increase in parallelism would enable the naı̈ve kernel to
further boost performance.
Figure 4 shows the speedups of each variant on the Phi over on multiple
CPU cores for each dataset. From Figure 4, it is observed that the Phi instance
of each variant outperforms its corresponding single-threaded and 16-threaded
CPU instances for each dataset. Specifically, the Phi instance of the naı̈ve
variant runs 25.34× faster on average than its single-threaded CPU instance,
with the maximum speedup 25.67, and 1.85× faster on average than the 16threaded CPU instance, with the maximum speedup 1.87. Compared to their
single-threaded CPU instances, the Phi instances of the GSE and VSE variants
produce the average speedups of 15.52 and 14.96 and the maximum speedups
of 15.77 and 16.00, respectively. Meanwhile, in comparison to their 16-threaded
CPU instances, their Phi instances performs 1.13× and 1.08× better on average
and 1.15× and 1.15× better at maximum, respectively.
5.2. Comparison With MATLAB and R
Secondly, we compared our algorithm with all-pairs τ coefficient implementations in the widely used MATLAB (version R2015b) and R (version 3.2.0).
MATLAB and R are executed in one 16-core compute node mentioned above,
where MATLAB runs 16 threads as it supports multithreading and R runs in
sequential. For fair comparison, both MATLAB and R merely compute all-pairs
21

Naïve

30

GSE

VSE

Naïve

2.0

GSE

VSE

1.8
25

1.6
1.4

Speedup

Speedup

20
15
10

1.2
1.0
0.8

0.6
0.4

5

0.2
0

0.0
DS21050

DS19784

DS13070

DS21050

DS3526

(a)

DS19784

DS13070

DS3526

(b)

Figure 4: Speedups of each variant on Phi over on CPU: (a) speedups over one CPU core and
(b) speedups over 16 CPU cores.
Table 5: Performance comparison with MATLAB and R

Dataset
DS21050
DS19784
DS13070
DS3526

Time (s)
MLAB R
182,343 266,001
205,181 284,215
205,063 287,635
240,778 345,668

Speedup over MLAB
Naı̈ve GSE VSE
456
435
699
484
473
772
473
446
743
476
494
812

Speedup over R
Naı̈ve GSE VSE
665
635
1,019
671
655
1,069
664
626
1,042
683
709
1,166

τ coefficient values, without conducting statistical tests, same as our algorithm
does. Table 5 shows the runtimes of the 16-threaded MATLAB and sequential
R instances as well as their comparison with our three variants that execute on
the 5110P Phi. From the table, each variant demonstrates excellent speedups,
i.e. two orders-of-magnitude over 16-threaded MATLAB and three orders-ofmagnitude over sequential R. Compared to 16-threaded MATLAB, the average
and maximum speedups are 472 and 484 for the naı̈ve variant, 462 and 494 for
the GSE variant and 756 and 812 for the VSE variant, respectively. In contrast,
the average and maximum speedups over sequential R are 671 and 683 for the
naı̈ve variant, 656 and 709 for the GSE variant and 1,074 and 1,166 for the VSE
variant, respectively.
5.3. Parallel Scalability Assessment
5.3.1. Multithreading
Thirdly, we evaluated the multithreading scalability of our variants on the
Phi with respect to number of threads. Figure 5 shows the parallel scalability
of the three variants. In this test, we used balanced thread affinity to ensure
that thread allocation is balanced over the cores and the threads allocated to
the same core have consecutive identifiers (i.e. neighbors of each other). This

22

59T

118T

177T

236T

1400

59T

118T

177T

1400

236T

1200

1200

1000

1000

1000

800
600

Time (s)

1200

Time (s)

Time (s)

1400

800
600

177T

236T

600

400

400

200

200

200

0

0

(a)

118T

800

400

DS21050 DS19784 DS13070 DS3526

59T

0
DS21050 DS19784 DS13070 DS3526

(b)

DS21050 DS19784 DS13070 DS3526

(c)

Figure 5: Multithreading scalability of different variants: (a) the naı̈ve variant, (b) the GSE
variant and (c) the VSE variant.

thread affinity configuration is attained by setting the environment variable
KMP AFFINITY to balanced.
From the figures, it is observed that each variant gets performance improved
as the number of active threads per core grows from 1, via 2 and 3, to 4 (corresponding to 59, 118, 177 and 236 threads, respectively). As each core is dual
issue and employs in-order execution, at least two threads per core are required
in order to saturate the compute capacity of each core. Meanwhile, when moving
from one thread per core to two threads per core, we expect that the speedup
could be close to 2. In this regard, we investigated the speedup of the instance
with two threads per core over the one with one thread per core for each variant, and found that the average speedups are 1.58, 1.61 and 1.56 for the naı̈ve,
GSE and VSE variants, respectively. This finding is close to our expectation
largely. Furthermore, since all threads per core share the dual in-order execution pipeline, deploying three or four threads per core may further improve
performance but normally with decreased parallel efficiency (with respect to
threads rather than cores). Note that for the 5110P Phi, only 59 out of 60 cores
are used for computing as one core is reserved for the operating system running
inside. Since each variant reaches peak performance at 236 threads, we used
this number of threads for performance evaluation in all tests, unless otherwise
stated.
5.3.2. Distributed computing
Finally, we measured the parallel scalability of our algorithm by varying the
number of Phis used in a distributed environment (see Figure 6). This test is
conducted in a cluster that consists of 16 Phis and is constituted by the aforementioned 8 compute nodes. Each variant used 236 threads, since this setting
leads to the best performance as mentioned above. From Figure 6, the naı̈ve
variant demonstrates nearly constant speedups for all datasets on a specific
number of Phis, while the two sorting-enabled variants exposed slight fluctuations under the same hardware configuration. This can be explained by the fact
that the runtime of our naı̈ve variant is independent of vector content, whereas
that of each sorting-enabled variant is sensitive to actual vector content to some
degree. Moreover, for each dataset, it is observed that every variant has demon23

18

18
DS21050
DS19784
DS13070
DS3526

14

Speedup

12

DS21050
DS19784
DS13070
DS3526

16
14
12

Speedup

16

10
8

10
8

6

6

4

4

2

2

0

0
2

4
8
Number of Phis

16

2

4
8
Number of Phis

(a)

(b)

DS21050
DS19784
DS13070
DS3526

12

GSE

VSE

500
400

Time (s)

14

Speedup

Naïve

600

18
16

16

10
8

300
200

6
4

100

2

0

0
2

4
8
Number of Phis

1

16

(c)

2

4
8
Number of Phis

16

(d)

Figure 6: Distributed computing scalability: (a) the naı̈ve variant, (b) the GSE variant and
(c) the VSE variant and (d) runtimes of each variant on the DS3526 dataset.

strated rather good scalability with respect to number of Phis. Concretely, the
naı̈ve variant achieves an average speedup of 2.00, 3.99, 7.96 and 15.81, the GSE
one yields an average speedup of 1.93, 3.82, 7.63 and 14.97, and the VSE one
produces an average speedup of 1.99, 3.94, 7.86 and 15.23, on 2, 4, 8 and 16
Phis, respectively.
6. Conclusion
Pairwise association measure is an important operation in searching for
meaningful insights within a dataset by examining potentially interesting relationships between data variables of the dataset. However, all-pairs association
computation are computationally demanding for a large volume of data and may
take a prohibitively long sequential runtime. This computational challenge for
big data motivated us to use parallel/high-performance computing architectures
to accelerate its execution.
In this paper, we have investigated the parallelization of all-pairs Kendall’s τ
coefficient computation using MIC processors as the accelerator. To the best of
our knowledge, this is for the first time that all-pairs τ coefficient computation
24

has been parallelized on MIC processors. For the τ coefficient, we have developed
three variants, namely the naı̈ve variant, the GSE variant and the VSE variant,
based on three pairwise τ coefficient kernels, i.e. the naı̈ve kernel, the GSE
kernel and the VSE kernel. The three variants all can execute on CPUs, Phis and
their clusters. The naı̈ve kernel has a time complexity of O(n2 ), while the other
two sorting-enabled kernels have an improved time complexity of O(n log n).
Furthermore, we have proposed a generic framework for workload balancing
in symmetric all-pairs computation and overcoming memory limitation on the
host or accelerator side. This framework assigns unique identifiers to jobs in the
upper triangle of the job matrix and builds provable bijective functions between
job identifier and coordinate space.
The performance of the three variants was evaluated using a collection of
real gene expression datasets produced from whole human genomes. Our experimental results demonstrated that on one 5110P Phi, the naı̈ve variant runs up
to 25.67× faster than its execution on a single CPU core, and up to 1.87× faster
than on 16 CPU cores. On the same Phi, the GSE and VSE variants run up to
15.77× and 16.00× faster than their executions on a single CPU core, and up
to 1.15× and 1.15× faster than on 16 CPU cores, respectively. Meanwhile, on
the same CPU/Phi hardware configuration, the VSE variant achieves superior
performance to the naı̈ve and GSE variants. Interestingly, the naı̈ve variant was
observed to be inferior to the GSE variant on CPUs for both single-threaded
and 16-threaded settings, but became superior to the latter on the 5110P Phi
for most benchmarking datasets used. As both variants use the same core
C++ code for their corresponding CPU- and MIC-oriented implementations as
mentioned above, this observation suggests that the architectural differences
between different types of PUs may make substantial impact on the resulting
performance. Subsequently, we compared our algorithm with the third-party
counterparts implemented in the popular MATLAB and R, observing that our
algorithm on a single 5110P Phi runs up to 812× faster than 16-threaded MATLAB and up to 1, 166× faster than sequential R, both of which are executed
on high-end Intel E5-2650 CPUs. We further assessed the parallel scalability
of our algorithm with respect to number of Phis in a distributed environment,
exposing that our algorithm demonstrated rather good distributed computing
scalability.
Finally, it is worth mentioning that the current version of our VSE kernel
constrains the largest n value to 215 − 1, due to 32-bit integer packing. One
solution to overcome this constraint is packing a variable pair (ui , vi ) into a
64-bit integer with each variable occupying 32 bits. In this case, we can split
a 512-bit vector into 8 lanes with each lane representing a 64-bit integer, and
then implement vectorized pairwise merge of sorted subarrays based on 64-bit
integers. Nonetheless, it should be noted that this 64-bit solution has twice less
parallelism than 32-bit. In the future, we plan to combine this non-linear measure with other linear or non-linear correlation/dependence measures to generate
fused co-expression networks for genome-wide gene expression data analysis at
population level. In addition, since R programming language is frequently used
in data analytics and the R implementation of all-pairs τ coefficient computation
25

is extremely slow for large-scale datasets, we also plan to release a R package
for public use based on our research in this study.
Acknowledgment
This research is supported in part by US National Science Foundation under
IIS-1416259 and an Intel Parallel Computing Center award. Conflict of interest :
none declared.
References
References
[1] Q. Zhu, A. K. Wong, A. Krishnan, M. R. Aure, A. Tadych, R. Zhang, D. C.
Corney, C. S. Greene, L. A. Bongo, V. N. Kristensen, et al., Targeted exploration and analysis of large cross-platform human transcriptomic compendia, Nature Methods 12 (3) (2015) 211–214.
[2] R. Steuer, J. Kurths, C. O. Daub, J. Weise, J. Selbig, The mutual information: detecting and evaluating dependencies between variables, Bioinformatics 18 (suppl 2) (2002) S231–S240.
[3] A. J. Butte, I. S. Kohane, Unsupervised knowledge discovery in medical
databases using relevance networks., in: Proceedings of the AMIA Symposium, American Medical Informatics Association, 1999, p. 711.
[4] M. Mutwil, S. Klie, T. Tohge, F. M. Giorgi, O. Wilkins, M. M. Campbell, A. R. Fernie, B. Usadel, Z. Nikoloski, S. Persson, Planet: combined
sequence and expression comparisons across plant networks derived from
seven species, The Plant Cell 23 (3) (2011) 895–910.
[5] A. Gobbi, G. Jurman, A null model for pearson coexpression networks,
PloS One 10 (6) (2015) e0128115.
[6] A. A. Margolin, I. Nemenman, K. Basso, C. Wiggins, G. Stolovitzky, R. D.
Favera, A. Califano, ARACNE: an algorithm for the reconstruction of gene
regulatory networks in a mammalian cellular context, BMC Bioinformatics
7 (Suppl 1) (2006) S7.
[7] M. Aluru, J. Zola, D. Nettleton, S. Aluru, Reverse engineering and analysis
of large genome-scale gene networks, Nucleic Acids Research 41 (1) (2013)
e24.
[8] A. Lachmann, F. M. Giorgi, G. Lopez, A. Califano, Aracne-ap: gene network reverse engineering through adaptive partitioning inference of mutual
information, Bioinformatics 32 (14) (2016) 2233–2235.
[9] J. Lee Rodgers, W. A. Nicewander, Thirteen ways to look at the correlation
coefficient, The American Statistician 42 (1) (1988) 59–66.
26

[10] L. Song, P. Langfelder, S. Horvath, Comparison of co-expression measures:
mutual information, correlation, and model based indices, BMC Bioinformatics 13 (1) (2012) 328.
[11] C. Spearman, Spearmans rank correlation coefficient, Amer J Psychol 15
(1904) 72–101.
[12] M. G. Kendall, Rank correlation methods, Griffin, 1948.
[13] Y. Wang, Y. Li, H. Cao, M. Xiong, Y. Y. Shugart, L. Jin, Efficient test
for nonlinear dependence of two continuous variables, BMC Bioinformatics
16 (1) (2015) 1.
[14] W. Xu, Y. Hou, Y. Hung, Y. Zou, A comparative analysis of spearman’s
rho and kendall’s tau in normal and contaminated normal models, Signal
Processing 93 (1) (2013) 261–276.
[15] G. A. Darbellay, I. Vajda, et al., Estimation of the information by an
adaptive partitioning of the observation space, IEEE Transactions on Information Theory 45 (4) (1999) 1315–1321.
[16] C. O. Daub, R. Steuer, J. Selbig, S. Kloska, Estimating mutual information
using b-spline functions–an improved similarity measure for analysing gene
expression data, BMC Bioinformatics 5 (1) (2004) 1.
[17] K.-C. Liang, X. Wang, Gene regulatory network reconstruction using conditional mutual information, EURASIP Journal on Bioinformatics and Systems Biology 2008 (1) (2008) 1–14.
[18] G. J. Székely, M. L. Rizzo, N. K. Bakirov, et al., Measuring and testing
dependence by correlation of distances, The Annals of Statistics 35 (6)
(2007) 2769–2794.
[19] M. R. Kosorok, On brownian distance covariance and high dimensional
data, The Annals of Applied Statistics 3 (4) (2009) 1266.
[20] A. Gretton, O. Bousquet, A. Smola, B. Schölkopf, Measuring statistical
dependence with hilbert-schmidt norms, in: International Conference on
Algorithmic Learning Theory, Springer, 2005, pp. 63–77.
[21] D. N. Reshef, Y. A. Reshef, H. K. Finucane, S. R. Grossman, G. McVean,
P. J. Turnbaugh, E. S. Lander, M. Mitzenmacher, P. C. Sabeti, Detecting
novel associations in large data sets, Science 334 (6062) (2011) 1518–1524.
[22] L. Song, A. Smola, A. Gretton, J. Bedo, K. Borgwardt, Feature selection via
dependence maximization, Journal of Machine Learning Research 13 (May)
(2012) 1393–1434.
[23] B. H. W. Chang, W. Tian, Gsa-lightning: ultra-fast permutation-based
gene set analysis, Bioinformatics 32 (19) (2016) 3029–3031.
27

[24] H. Abdi, The kendall rank correlation coefficient, Encyclopedia of Measurement and Statistics. Sage, Thousand Oaks, CA (2007) 508–510.
[25] S. Wang, I. Pandis, D. Johnson, I. Emam, F. Guitton, A. Oehmichen,
Y. Guo, Optimising parallel R correlation matrix calculations on gene expression data using mapreduce, BMC Bioinformatics 15 (1) (2014) 1.
[26] R. C. Team, et al., R: A language and environment for statistical computing, R Foundation for Statistical Computing, 2013.
[27] J. Dean, S. Ghemawat, Mapreduce: simplified data processing on large
clusters, Communications of the ACM 51 (1) (2008) 107–113.
[28] Y. Liu, T. Pan, S. Aluru, Parallel pairwise correlation computation on
intel xeon phi clusters, in: 28th International Symposium on Computer
Architecture and High Performance Computing, IEEE, 2016, pp. 141–149.
[29] Mathworks, https://www.mathworks.com/products/matlab.html (2015).
[30] A. Sodani, R. Gramunt, J. Corbal, H.-S. Kim, K. Vinod, S. Chinthamani,
S. Hutsell, R. Agarwal, Y.-C. Liu, Knights landing: Second-generation intel
xeon phi product, IEEE Micro 36 (2) (2016) 34–46.
[31] J. Jeffers, J. Reinders, Intel Xeon Phi coprocessor high-performance programming, Morgan Kaufmann, 2013.
[32] W. R. Knight, A computer method for calculating kendall’s tau with ungrouped data, Journal of the American Statistical Association 61 (314)
(1966) 436–439.
[33] H. Inoue, T. Moriyama, H. Komatsu, T. Nakatani, Aa-sort: A new parallel
sorting algorithm for multi-core simd processors, in: Proceedings of the
16th International Conference on Parallel Architecture and Compilation
Techniques, IEEE Computer Society, 2007, pp. 189–198.
[34] K. E. Batcher, Sorting networks and their applications, in: Proceedings of
the April 30–May 2, 1968, spring joint computer conference, ACM, 1968,
pp. 307–314.
[35] J. Chhugani, A. D. Nguyen, V. W. Lee, W. Macy, M. Hagog, Y.-K. Chen,
A. Baransi, S. Kumar, P. Dubey, Efficient implementation of sorting on
multi-core simd cpu architecture, Proceedings of the VLDB Endowment
1 (2) (2008) 1313–1324.
[36] T. Xiaochen, K. Rocki, R. Suda, Register level sort algorithm on multicore simd processors, in: Proceedings of the 3rd Workshop on Irregular
Applications: Architectures and Algorithms, ACM, 2013, p. 9.
[37] S. Odeh, O. Green, Z. Mwassi, O. Shmueli, Y. Birk, Merge path-cacheefficient parallel merge and sort, Tech. rep., Technical report, CCIT Report
(2012).
28

[38] O. Green, R. McColl, D. A. Bader, Gpu merge path: a gpu merging algorithm, in: Proceedings of the 26th ACM international conference on
Supercomputing, ACM, 2012, pp. 331–340.
[39] Intel, Intel xeon phi processor software optimization guide,
https://software.intel.com/en-us/articles/intel-xeon-phi-processorsoftware-optimization-guide (2016).
[40] T. Kiefer, P. B. Volk, W. Lehner, Pairwise element computation with
MapReduce, in: Proceedings of the 19th ACM International Symposium
on High Performance Distributed Computing, ACM, 2010, pp. 826–833.
[41] W. Gropp, T. Hoefler, R. Thakur, E. Lusk, Using advanced MPI: Modern
features of the message-passing interface, MIT Press, 2014.
[42] T. El-Ghazawi, L. Smith, Upc: unified parallel c, in: Proceedings of the
2006 ACM/IEEE conference on Supercomputing, ACM, 2006, p. 27.
[43] J. González-Domı́nguez, Y. Liu, B. Schmidt, Parallel and scalable shortread alignment on multi-core clusters using upc++, PLoS ONE 11 (1)
(2016) e0145490.

29

