* [126] NA49 Collaboration, C. Alt et al., “Energy dependence of $\Lambda$ and $\Xi$ production in central Pb+Pb collisions at A-20, A-30, A-40, A-80, and A-158 GeV measured at the CERN Super Proton Synchrotron,” Phys. Rev. C 78 (2008) 034918, arXiv:0804.3770 [nucl-ex].
  * [127] NA61/SHINE Collaboration, A. Acharya et al., “Measurements of $\pi^{\pm}$, $K^{\pm}$, $p$ and $\bar{p}$ spectra in 7Be+9Be collisions at beam momenta from 19$A$ to 150$A$ GeV/$c$ with the NA61/SHINE spectrometer at the CERN SPS,” Eur. Phys. J. C 81 no. 1, (2021) 73, arXiv:2010.01864 [hep-ex].
  * [128] T. T. Böhlen, F. Cerutti, M. P. W. Chin, A. Fassò, A. Ferrari, P. G. Ortega, A. Mairani, P. R. Sala, G. Smirnov, and V. Vlachoudis, “The FLUKA Code: Developments and Challenges for High Energy and Medical Applications,” Nucl. Data Sheets 120 (2014) 211.
  * [129] A. Ferrari, P. R. Sala, A. Fasso, and J. Ranft, “FLUKA: A multi-particle transport code (Program version 2005),” Tech. Rep. CERN-2005-010, SLAC-R-773, INFN-TC-05-11, 2005.  http://cds.cern.ch/record/898301.
  * [130] GEANT4 Collaboration, S. Agostinelli et al., “GEANT4–a simulation toolkit,” Nucl. Instrum. Meth. A 506 (2003) 250–303.
  * [131] J. Allison et al., “Geant4 developments and applications,” IEEE Trans. Nucl. Sci. 53 (2006) 270.
  * [132] G. Battistoni et al., “Overview of the FLUKA code,” Annals Nucl. Energy 82 (2015) 10–18.
  * [133] Particle Data Group, M. Tanabashi, et al., “Review of Particle Physics,” Phys. Rev. D98 (2018) 030001.
  * [134] S. Alioli, P. Nason, C. Oleari, and E. Re, “A general framework for implementing NLO calculations in shower Monte Carlo programs: the POWHEG BOX,” JHEP 06 (2010) 043, arXiv:1002.2581 [hep-ph].
  * [135] D. J. Lange, “The EvtGen particle decay simulation package,” Nucl. Instrum. Meth. A462 (2001) 152.
  * [136] NA49 Collaboration, S. V. Afanasiev et al., “Energy dependence of pion and kaon production in central Pb+Pb collisions,” Phys. Rev. C66 (2002) 054902, arXiv:nucl-ex/0205002 [nucl-ex].
  * [137] NA49 Collaboration, C. Alt et al., “Energy and centrality dependence of anti-p and p production and the anti-Lambda/anti-p ratio in Pb+Pb collisions between 20 AGeV and 158 AGeV,” Phys. Rev. C73 (2006) 044910.
  * [138] C. Lourenço and H. K. Wöhri, “Heavy flavour hadro-production from fixed-target to collider energies,” Phys. Rept. 433 (2006) 127, arXiv:hep-ph/0609101 [hep-ph].
  * [139] Hard Probe Collaboration, R. Vogt, “The A dependence of open charm and bottom production,” Int. J. Mod. Phys. E12 (2003) 211, arXiv:hep-ph/0111271 [hep-ph].
  * [140] ALICE Collaboration, “Letter of intent for ALICE 3: A next-generation heavy-ion experiment at the LHC,” arXiv:2211.02491 [physics.ins-det].
  * [141] NA49 Collaboration, C. Alt et al., “Energy dependence of $\phi$ meson production in central Pb+Pb collisions at $\sqrt{s_{\rm NN}}=6$ to 17 GeV,” Phys. Rev. C 78 (2008) 044907, arXiv:0806.1937 [nucl-ex].
  * [142] NA49 Collaboration, M. Mitrovski, “Omega and anti-Omega production in Pb+Pb and p + p collisions at 30-A-Gev, 40-A-Gev and 158-A-Gev,” J. Phys. G 30 (2004) S357–S362, arXiv:nucl-ex/0305029.
  * [143] S. B. Kotsiantis, “Decision trees: a recent overview,” Artificial Intelligence Review 39 no. 4, (2013) 261–283.
  * [144] “Xgboost documentation.” https://xgboost.readthedocs.io/en/latest/.
  * [145] J. Podolanski and R. Armenteros, “Analysis of V-events,” The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 45 no. 360, (1954) 13–30.
  * [146] D. Lissauer and E. V. Shuryak, “K meson modification in hot hadronic matter may be detected via $\phi$ meson decays,” Physics Letters B 253 no. 1-2, (1991) 15–18.
  * [147] F. Klingl, T. Waas, and W. Weise, “Modification of the $\phi$ meson spectrum in nuclear matter,” Phys. Lett. B 431 (1998) 254–262, arXiv:hep-ph/9709210.
  * [148] NA49 Collaboration, V. Friese, “Production of strange resonances in C + C and Pb + Pb collisions at 158-A-GeV,” Nucl. Phys. A 698 (2002) 487–490.
  * [149] NA49 Collaboration, C. Alt et al., “System-size dependence of strangeness production in nucleus-nucleus collisions at $\sqrt{s_{\rm NN}}=17.3$ GeV,” Phys. Rev. Lett. 94 (2005) 052301, arXiv:nucl-ex/0406031.
  * [150] NA50 Collaboration, B. Alessandro et al., “$\phi$ production in Pb-Pb collisions at 158-GeV/c per nucleon incident momentum,” Phys. Lett. B 555 (2003) 147–155. [Erratum: Phys.Lett.B 561, 294–294 (2003)].
  * [151] NA50 Collaboration, D. Jouan, “Phi and omega - rho production in d-C, d-U, S-U and Pb-Pb at SPS energies,” J. Phys. G 35 (2008) 104163.
  * [152] ALICE Collaboration, B. Abelev et al., “Technical Design Report for the Upgrade of the ALICE Inner Tracking System,” J. Phys. G41 (2014) 087002.
  * [153] FINUDA Collaboration, M. Agnello et al., “New results on Mesonic Weak Decay of p-shell Lambda-Hypernuclei,” Phys. Lett. B 681 (2009) 139–146, arXiv:0905.0623 [nucl-ex].
  * [154] V. Vovchenko and H. Stoecker, “Thermal-FIST: A package for heavy-ion collisions and hadronic equation of state,” Comput. Phys. Commun. 244 (2019) 295–310, arXiv:1901.05249 [nucl-th].
  * [155] V. Vovchenko, V. Koch, and C. Shen, “Proton number cumulants and correlation functions in Au-Au collisions at sNN=7.7–200 GeV from hydrodynamics,” Phys. Rev. C 105 no. 1, (2022) 014904, arXiv:2107.00163 [hep-ph].
  * [156] D. Gazda, A. Pérez-Obiol, A. Gal, and E. Friedman, “Lifetime of the hypertriton,” EPJ Web Conf. 271 (2022) 01002, arXiv:2211.17122 [nucl-th].
  * [157] F. Becattini, J. Manninen, and M. Gazdzicki, “Energy and system size dependence of chemical freeze-out in relativistic nuclear collisions,” Phys. Rev. C73 (2006) 044905, arXiv:hep-ph/0511092 [hep-ph].
  * [158] HADES Collaboration, S. Harabasz, “Multi-differential pattern of low-mass $\mathrm{e}^{+}\mathrm{e}^{-}$ excess from $\sqrt{s_{\rm NN}}=2.4$ GeV Au$+$Au collisions with HADES,” Nucl. Phys. A982 (2019) 771.
  * [159] NA57 Collaboration, F. Antinori et al., “Multiplicity of charged particles in Pb Pb collisions at SPS energies,” J. Phys. G 31 (2005) 321–335.
  * [160] R. Rapp, “Dilepton Spectroscopy of QCD Matter at Collider Energies,” Adv. High Energy Phys. 2013 (2013) 148253, arXiv:1304.2309 [hep-ph].
  * [161] R. Vogt, “$\mathrm{J}/\psi$ production and suppression,” Phys. Rept. 310 (1999) 197.
  * [162] NA50 Collaboration, B. Alessandro et al., “Charmonia and Drell-Yan production in proton–nucleus collisions at the CERN SPS,” Phys. Lett. B553 (2003) 167.
  * [163] NA50 Collaboration, B. Alessandro et al., “$psi^{\prime}$ production in Pb–Pb collisions at 158 GeV/nucleon,” Eur. Phys. J. C49 (2007) 559, arXiv:nucl-ex/0612013 [nucl-ex].
  * [164] K. Banicz, F. Bergsma, and H. K. Wöhri, “The PT7 magnet,” Tech. Rep. NA60-Note-2001-2, 2001.  http://na60web.ca.infn.it/www/experiment/notes/note_2001_2.pdf.
  * [165] G. Aglieri Rinella et al., “The ALPIDE pixel sensor chip for the upgrade of the ALICE Inner Tracking System,” Nucl. Instrum. Meth. A845 (2017) 583. Proceedings of the Vienna Conference on Instrumentation 2016.
  * [166] ALICE Collaboration, “Expression of Interest for an ALICE ITS Upgrade in LS3,”. https://cds.cern.ch/record/2644611/.
  * [167] “SoftwAnsys Engineering Software®.”  https://www.ansys.com/.
  * [168] “COMSOL Multiphysics®.”  https://www.comsol.it/.
  * [169] ALICE Collaboration, P. Antonioli, A. Kluge, and W. Riegler, “Upgrade of the ALICE Readout; Trigger System,” tech. rep., 2013.  https://cds.cern.ch/record/1603472.
  * [170] F. Sauli, “GEM: A new concept for electron amplification in gas detectors,” Nucl. Instrum. Meth. A 386 (1997) 531–534.
  * [171] B. Ketzer, Q. Weitzel, S. Paul, F. Sauli, and L. Ropelewski, “Performance of triple GEM tracking detectors in the COMPASS experiment,” Nucl. Instrum. Meth. A 535 (2004) 314–318. [Erratum: Nucl.Instrum.Meth.A 648, 293 (2011)].
  * [172] K. Gnanvo, X. Bai, C. Gu, N. Liyanage, V. Nelyubin, and Y. Zhao, “Performance in test beam of a large-area and light-weight GEM detector with 2D stereo-angle (U–V) strip readout,” Nucl. Instrum. Meth. A 808 (2016) 83–92, arXiv:1509.03875 [physics.ins-det].
  * [173] A. Zhang, V. Bhopatkar, E. Hansen, M. Hohlmann, S. Khanal, M. Phipps, E. Starling, J. Twigger, and K. Walton, “Performance of a Large-area GEM Detector Read Out with Wide Radial Zigzag Strips,” Nucl. Instrum. Meth. A 811 (2016) 30–41, arXiv:1508.07046 [physics.ins-det].
  * [174] A. Gasparian et al., “PRad-II: A new upgraded high precision measurement of the proton charge radius,” Jefferson Lab 2020 Proposals PAC48 (2020) . https://www.jlab.org/exp_prog/proposals/20/PR12-20-004_Proposal.pdf.
  * [175] H. Klest, “Overview and design of the sPHENIX TPC,” J. Phys. Conf. Ser. 1498 (2020) 012025, arXiv:2110.02082 [physics.ins-det].
  * [176] MOLLER Collaboration, J. Benesch et al., “The MOLLER Experiment: An Ultra-Precise Measurement of the Weak Mixing Angle Using Möller Scattering,” arXiv:1411.4088 [nucl-ex].
  * [177] H. Danielsson et al., “NA60+ Demonstrator, design, construction, simulation, experimental investigation, and lessons learned,” Tech. Rep. EDMS doc. 2694487 v.1, CERN, Geneva, 2022.  https://edms.cern.ch/document/2694487/1.
  * [178] LHCb Collaboration Collaboration, S. Amato et al., LHCb magnet: Technical Design Report.  Technical design report. LHCb. CERN, Geneva, 2000.  https://cds.cern.ch/record/424338.
  * [179] ALICE Collaboration, M. Puccio, “Tracking in high-multiplicity events,” PoS Vertex2016 (2017) 043.
  * [180] A. Gerbershagen, C. Ahdida, J. Bernhard, V. Clerc, S. Girod, E. Scomparin, G. Usai, and H. Vincke, “Design of beam optics and radiation protection concept for NA60+ heavy-ion experiment at CERN,” Nucl. Instrum. Meth. A 1047 (2023) 167887, arXiv:2210.17527 [physics.acc-ph].
  * [181] D. Forkel-Wirth, “Zonage radiologique au CERN,” EDMS 810149 .
  * [182] C. Ahdida et al., “New Capabilities of the FLUKA Multi-Purpose Code,” Front. in Phys. 9 (2022) 788253.
  * [183] V. Vlachoudis, “Flair: A powerful but user friendly graphical interface for FLUKA,” in International Conference on Mathematics, Computational Methods & Reactor Physics 2009, pp. 790–800.  2009\.
  * [184] H. Vincke and C. Theis, “ActiWiz – optimizing your nuclide inventory at proton accelerators with a computer code,” Prog. Nucl. Sci. Tech. 4 (2014) 228–232.
  * [185] Swiss Federal Council, “Radiological Protection Ordinance 814.501,”.
  * [186] F. Aberle, “Radiation monitoring layout - SPS North Area,” EDMS 2211938 .
  * [187] H. Boukabache, M. Pangallo, G. Ducos, N. Cardines, A. Bellotta, C. Toner, D. Perrin, and D. Forkel-Wirth, “Towards a Novel Modular Architecture for CERN Radiation Monitoring,” Radiat. Prot. Dosim. 173 no. 1-3, (2017) 240–244.

## Appendix: NA60+ Collaboration

C. Ahdidacern&G. Aloccoucagliaricagliari&F. Antinoripadova&M. Arbacagliari&M.
Arestiucagliaricagliari&R. Arnalditorino&A. Baratto Roldancern&S.
Beolèutorinotorino&A. Beraudotorino&J. Bernhardcern&L. Bianchiutorinotorino&M.
Borysovaweizmannkyiv&S. Bresslerweizmann&S. Bufalinodisattorino&E.
Casulaucagliaricagliari&C. Cicalòcagliari&S. Colitorino&P.
Cortesepiemontetorino&A. Dainesepadova&H. Danielssoncern&A. De
Falcoucagliaricagliari&K. Dehmeltstonybrook&A. Dreesstonybrook&A.
Ferrettiutorinotorino&F. Fiondaucagliaricagliari&M. Gagliardiutorinotorino&A.
Gerbershagengroningen&F. Geurtsrice&V. Grecoucataniacatania&W. Lirice&M.P.
Lombardolfn&D. Marrascagliari&M. Maserautorinotorino&A. Masonicagliari&L.
Micheletticern&L. Mirasolaucagliaricagliari&F. Mazzaschicernutorino&M.
Mentinkcern&P. Mereutorino&A. Milovweizmann&A. Mulliriucagliaricagliari&L.
Musacern&C. Oppedisanotorino&B. Paulucagliaricagliari&M.
Pennisiutorinotorino&S. Plumariucatania&F. Prinotorino&M. Pucciocern&C.
Puggionicagliari&R. Rapptamu&I. Ravinovichweizmann&A. Rossipadova&V.
Sarritzuucagliaricagliari&B. Schmidtcern&E. Scomparintorino&S.
Siddhantacagliari&R. Shahoyancern&M. Tuvericagliari&A. Uraslyon&G.
Usaiucagliaricagliari&H. Vinckecern&I. Vorobyevcern

cernEuropean Organization for Nuclear Research (CERN), Geneva, Switzerland

ucagliariDipartimento di Fisica dell’Università di Cagliari, Cagliari, Italy

cagliariINFN, Sezione di Cagliari, Cagliari, Italy

padovaINFN, Sezione di Padova, Padova, Italy

torinoINFN, Sezione di Torino, Turin, Italy

utorinoDipartimento di Fisica dell Università di Torino, Turin, Italy

weizmannDepartment of Particle Physics and Astrophysics, Weizmann Insitute of
Science, Rehovot, Israel

kyivKyiv Institute for Nuclear Research (KINR), Natl. Acad. of Sci. of Ukraine
(NASU)

disatDipartimento DISAT del Politecnico di Torino, Turin, Italy

piemonteDipartimento di Scienze e Innovazione Tecnologica dell’Università del
Piemonte Orientale,
Alessandria, Italy

stonybrookDepartment of Physics and Astronomy, Stony Brook University, SUNY,
Stony Brook, New York, USA

groningenDepartment of Radiation Oncology, University of Groningen, Groningen,
The Netherlands

riceDepartment of Physics and Astronomy, Rice University, Houston, Texas, USA

ucataniaDipartimento di Fisica e Astronomia dell’Università di Catania,
Catania, Italy

cataniaINFN, Laboratori Nazionali del Sud, Catania, Italy

lfnINFN, Laboratori Nazionali di Frascati, Frascati, Italy

tamuCyclotron Institute and Department of Physics and Astronomy, Texas A&M
University, College Station, Texas, USA

lyonInstitut de Physique des 2 Infinis de Lyon, Université de Lyon,
CNRS/IN2P3, Lyon, France


Sequential and Shared-Memory Parallel Algorithms for Partitioned Local Depths
Aditya Devarakonda[Department of Computer Science, Wake Forest University.]

Grey Ballard[1]

In this work, we design, analyze, and optimize sequential and shared-memory parallel algorithms for partitioned local depths (PaLD).
Given a set of data points and pairwise distances, PaLD is a method for identifying strength of pairwise relationships based on relative distances, enabling the identification of strong ties within dense and sparse communities even if their sizes and within-community absolute distances vary greatly.
We design two algorithmic variants that perform community structure analysis through triplet comparisons of pairwise distances.
We present theoretical analyses of computation and communication costs and prove that the sequential algorithms are communication optimal, up to constant factors.
We introduce performance optimization strategies that yield sequential speedups of up to $29\times$ over a baseline sequential implementation and parallel speedups of up to $19.4\times$ over optimized sequential implementations using up to $32$ threads on an Intel multicore CPU.

§ INTRODUCTION.

Partitioned local depths (PaLD) is a method for revealing community structure in distance-based data [2].
Given pairwise distances (or dissimilarities) of a set of points, PaLD computes another pairwise measure called cohesion that measures closeness based on relative distances.
By relying on relative distance, PaLD is able to use a universal threshold to distinguish between strong and weak ties without defining neighborhoods by a single number of neighborhoods, neighborhood size, or absolute distance threshold.
In this way, PaLD can identify neighborhoods of varying size and density, making it useful for data where the relationships among points behave differently across the space.

The input to PaLD is a distance matrix, and the output is a cohesion matrix.
As detailed in <ref>, computing cohesion requires determining the size of the local neighborhood of each pair of points and then computing contributions to cohesion values based on neighborhood sizes.
In each case, the fundamental operation is a comparison of the pairwise distances among triplets of points.
Given $n$ points, this yields an arithmetic complexity of $O(n^3)$.
The goal of this paper is to develop efficient sequential and shared-memory parallel algorithms for scaling PaLD to datasets of size up to $O(10^5)$, making it computationally feasible to analyze ones that fit in memory on a single server.
<Ref> presents the structure of the PaLD computation and our two main algorithmic approaches, which we call pairwise and triplet, respectively.
As an $O(n^3)$ computation, PaLD shares many similarities with dense matrix multiplication (GEMM), and our algorithmic design borrows from ideas of cache-efficient algorithms for GEMM [3, 9, 18].
For example, the basic computation is a comparison between distances of points $x,y,z$, which involves distance matrix entries $d_{xy}$, $d_{yz}$, and $d_{xz}$ and has an access pattern similar to the fused multiply-adds (FMA) within GEMM.
There are a few key differences between PaLD and GEMM.
First, because of symmetric distances, the order of the points is irrelevant, so rather than requiring consideration of all $n^3$ possible values of $x,y,z$, we need consider only $\binom{n}{3}\approx n^3/6$ unique triplets.
Second, while the memory access of distances is regular, the updates of the cohesion requires branching based on distance comparisons.
Finally, the computation requires two passes because cohesion updates depend on the sizes of local neighborhoods.
Each pass requires a varying mix of integer and floating point operations in addition to the branching.
The pairwise and triplet approaches navigate a tradeoff between exploiting symmetry and achieving regular data access and parallelization.

In <ref> we prove a lower bound on the cache efficiency of any PaLD algorithm, and we show that both of our algorithms achieve optimal cache performance, up to constant factors.
By exploiting symmetry and applying cache blocking, we obtain data locality in cache and minimize the number of reads and writes of matrix values.
<Ref> details our low-level optimizations of the two PaLD algorithms. We show that branch avoidance has the highest impact on sequential performance given the high cost of branch misprediction [10, 13, 14].
Along with other optimizations including cache blocking and vectorization, we show performance improvements over naive sequential code of up to $29\times$.
In <ref> we design, optimize, and evaluate OpenMP parallel versions of the two PaLD algorithms.
We show that the pairwise algorithm enables regular data access patterns and loop-based parallelism that can largely avoid write conflicts.
The triplet algorithm exploits more symmetry to reduce arithmetic operations but requires task-based parallelism due to more complicated data access patterns and write conflicts.
We also apply Non-Uniform Memory Access (NUMA) optimizations when scaling across sockets.
We achieve strong scaling speedups up to $19.4\times$ for pairwise and $13.2\times$ for triplet over their optimized sequential versions on $32$ threads.
Finally, we describe a text analysis application in <ref>, demonstrating the utility of PaLD on larger datasets than previously considered, and we show a parallel speedup of $16.7\times$ on a task with $n=2712$ using $32$ threads.

§ BACKGROUND.

Given a set of points and a pairwise distance metric, partitioned local depth (PaLD) algorithms determine the pairwise cohesion between all pairs of points in a dataset [2].
Assuming that the dataset comprises sufficiently separated subsets, cohesion values are invariant to contraction and dilation of distances within subset distances.
Community structure revealed by cohesion values capture the concept of near neighbors based on relative positioning, adapting to varying density.
This approach is more flexible than standard cluster labeling or nearest neighbor approaches.
Density-based approaches (e.g. DBSCAN) [5, 6, 7] that attempt to combine points into high- and low-density groups based on pairwise distances include thresholding (tuning) parameters to reflect locality and cluster size.
Likewise, $k$-nearest neighbor (KNN) approaches [8] attempt to group points via comparisons against their $k$ nearest neighbors (using absolute distances).
The tuning parameter, $k$, controls the neighborhood size for a given point and is often fixed for all points.
Cohesion values depend on triplet distance comparisons (as opposed to absolute distances) which require only measures of relative similarity and can be more reliable than exact numerical distances for analyzing high dimensional, non-Euclidean data.
PaLD requires $O(n^3)$ operations to compute cohesion values without assumptions on underlying probability distribution or tuning parameters.

Given a set of points $\SSS$, the local focus of a pair of points $x,y\in \SSS$ is the set of all points within distance $d_{xy}$ of either $x$ or $y$, where $d_{xy}$ is the distance between $x$ and $y$: $ {\cal U}_{xy} = \{z \in \SSS \; | \; d_{xz} \leq d_{xy} \; \Or \;  d_{yz} \leq d_{xy}\}.$
We let $u_{xy} = |{\cal U}_{xy}|$ denote the size of the local focus.

The local depth of a point $x \in \SSS$ is the probability that, given a uniformly chosen random second point $Y\in \cal S$ and a random third point $Z$ chosen uniformly from the local focus ${\cal U}_{xY}$, $Z$ is closer to $x$ than $Y$:

\begin{equation}
\label{eq:LD}
\ell_x = \Pr\left[d_{Zx} < d_{ZY} \; | \; Y \sim \unif(\SSS {\setminus} \{x\}), Z \sim \unif({\cal U}_{xY})\right].
\end{equation}

The cohesion of a point $z$ to another point $x$ is a part of the local depth $\ell_x$ and is defined as

\begin{equation}
\label{eq:cohesion}
c_{xz} = \Pr\left[Z = z \; \tAnd \; d_{Zx} < d_{ZY} \right].
\end{equation}

The random variables $Y$ and $Z$ in <ref> are chosen from the same distributions as in <ref>; we drop the notation here and later.
This implies that $\ell_x = \sum_{z\in \SSS} c_{xz}$, or that cohesion is partitioned local depth.
The cohesion matrix, $C$, can be used to analyze community structure.
For example, two points have particularly strong cohesion if the impact of one of the points to the other is more than that expected from a random focus point of another random point.

§ PALD ALGORITHMS DESIGN.

In order to compute the cohesion of all pairs of points, we can again use the law of total probability to partition $c_{xz}$ across all points $y\in {\cal S}$:

\begin{equation*}
c_{xz} = \sum_{y\in \SSS} \Pr\left[Y = y \; \tAnd \; Z = z \; \tAnd \; d_{Zx} < d_{ZY} \right].
\end{equation*}

Using the law of conditional probability, this becomes

\begin{multline*}
c_{xz} = \sum_{y\in \SSS} \Pr\left[d_{zx} < d_{zy} \; | \; Y = y, Z = z \right] \\
\cdot \Pr\left[Z=z \; | \; Y = y\right] \cdot \Pr\left[Y=y\right]
\end{multline*}

which implies

\begin{equation}
\label{eq:fxz}
c_{xz} = \sum_{y\in \SSS} \indic{d_{xz} \leq d_{yz}} \cdot \frac{\indic{d_{xz}\leq d_{xy}}}{u_{xy}} \cdot \frac{1}{n-1}
= \frac{1}{n-1} \sum_{y\in \SSS} g_{xyz},
\end{equation}

where $\indic{}$ is the indicator function and we have defined

\begin{equation}
\label{eq:gxyz}
g_{xyz} = \indic{d_{xz} \leq d_{yz}} \cdot \indic{d_{xz}\leq d_{xy}} \, / \, u_{xy}.
\end{equation}

The task is then to compute $g_{xyz}$ for all $x,y,z\in\SSS$, a total of $n^3$ values.
However, only approximately 1/3rd of the $g_{xyz}$ values are nonzero because, given three points with unique pairwise distance values, only one pair has the minimum distance.
For example, given points $x,y,z\in\SSS$, if $x$ and $y$ are the closest pair, then $g_{xzy}$ and $g_{yzx}$ are nonzero, but $g_{xyz}=g_{yxz}=g_{zxy}=g_{zyz}=0$.
To compute the nonzero values $g_{xzy}$ and $g_{yzx}$, we need the values $u_{xz}$ and $u_{yz}$.
The size of any given local focus can be computed as $u_{xy} = \sum_{z\in\SSS} \indic{d_{xz} \leq d_{xy} \,\Or d_{yz} \leq d_{xy}}.$

We consider two algorithmic approaches to computing the local focus sizes and the final cohesion matrix that take advantage of the symmetry.
The first approach, which we call the pairwise algorithm, considers all $\binom{n}{2}$ pairs of points, and for each pair, first determines the size of its local focus and then computes contributions to the cohesion matrix from all points within the local focus.
The second approach, which we call the triplet algorithm, considers all $\binom{n}{3}$ triplets of points, and for each triplet, determines which of the two local foci the triplet contributes to and then (in a second pass) determines which of the two cohesion matrix entries the triplet contributes to.
We analyze and compare the two algorithms in <ref>.

§.§ Pairwise Algorithm.

The entry-wise pairwise algorithm is given as <ref>.
The idea is to perform the computations for each pair of points $x$ and $y$.
To compute $g_{xyz}$ for each third point $z$, we first must compute the size of the local focus, $u_{xy}$.
This requires a pass over all $n$ points with two comparisons and a possible integer increment.
A second pass over all $n$ points determines, for points in the local focus, which of the points $x$ or $y$ the third point supports, and the cohesion matrix is updated accordingly.
Note that only one local focus size need be stored at any one time, requiring minimal temporary memory.

Pairwise Sequential Algorithm

[1]
$D \in \mathbb{R}^{n \times n},$ Distance Matrix
$C \in \mathbb{R}^{n \times n},$ Cohesion Matrix
$x = 1$ $n - 1$
$y = x + 1$ $n$

$u_{xy} = 0$
$z = 1$ $n$
$d_{xz} < d_{xy}$ $d_{yz} < d_{xy}$
$u_{xy} = u_{xy}+1$

$z = 1$ $n$
$d_{xz} < d_{xy}$ $d_{yz} < d_{xy}$
$d_{xz} < d_{yz}$
$c_{xz} = c_{xz} + 1/u_{xy}$
$c_{yz} = c_{yz} + 1/u_{xy}$

To improve the cache locality, we block the algorithm as follows: instead of considering only one pair of points, we consider two sets of points $\mathcal{X}$ and $\mathcal{Y}$ and consider all the pairs $(x,y)\in \mathcal{X} \times \mathcal{Y}$.
In this way, we obtain locality on the distance matrix block $D_{\mathcal{X},\mathcal{Y}}$ and a temporary block of local focus sizes $U_{\mathcal{X},\mathcal{Y}}$.

As in the entry-wise algorithm, the blocked algorithm makes two passes over all $n$ third points.
The first pass computes $U_{\mathcal{X},\mathcal{Y}}$, a local focus size block, and the second pass makes updates to the cohesion matrix.

Dependency structure of the blocked pairwise algorithm. The highlighted regions represent quantities with temporal locality. Quantities in red correspond to reads and ones in green correspond to writes. Orange entries are computed and used in fast memory and then discarded.  Blue represents entry-wise dependencies within each matrix/vector.

<Ref> shows the dependencies among the distance, local focus, and cohesion matrices for the blocked ($b = 4$) pairwise algorithm.
The red blocks correspond to the entries of the distance matrix that are read and re-used while processing the pair of blocks $\mathcal{X}$ and $\mathcal{Y}$ (the pattern is the same in both passes, though $D_{\mathcal{X},\mathcal{Y}}$ remains in fast memory through both passes.
The orange blocks represent entries of the local focus matrix which are computed in fast memory during the first pass and used during the second pass.
The green blocks of the cohesion matrix are re-used during the second pass before being written back to slow memory.
The blue blocks represent dependencies between entries of $D, U,$ and $C$ for one entry-wise iteration.

§.§ Triplet Algorithm.

The entry-wise triplet algorithm is given as <ref>.
In <Ref>, if a third point $z$ is in the local focus of $x$ and $y$ and is closer to $x$, then only the support of $z$ for $x$ is recorded in $C$ ($c_{xz}$ is updated).
If $z$ is closer to $x$ in a focus with $y$, then $x$ is closer to $z$ in its focus with $y$.
The idea of the triplet algorithm is to minimize the number of distance comparisons.
By performing all the updates for each triplet of points, we can avoid redundant comparisons.
However, this method requires that the local focus sizes are pre-computed for all pairs of points within the triplet, so it requires more temporary memory.
Triplet Sequential Algorithm.

[1]
$D \in \mathbb{R}^{n \times n}$ Distance Matrix
$C \in \mathbb{R}^{n \times n}$ Cohesion Matrix

Initialize $U = \texttt{triu}(2*\texttt{ones}(n),1)$
$x = 1$ $n-1$
$y = x + 1$ $n$
$z = y + 1$ $n$
$d_{xy} < d_{xz}$ $d_{xy} < d_{yz}$
// $x,y$ is closest pair
$u_{xz} = u_{xz} + 1$
$u_{yz} = u_{yz} + 1$
$ d_{xz} < d_{yz}$
// $x,z$ is closest pair
$u_{xy} = u_{xy} + 1$
$u_{yz} = u_{yz} + 1$
// $y,z$ is closest pair
$u_{xy} = u_{xy} + 1$
$u_{xz} = u_{xz} + 1$

$x = 1$ $n-1$
$y = x + 1$ $n$
$z = y + 1$ $n$
$d_{xy} < d_{xz} $ $ d_{xy} < d_{yz}$

$c_{xy} = c_{xy} + 1/u_{xz}$
$c_{yx} = c_{yx} + 1/u_{yz}$
$ d_{xz} < d_{yz}$

$c_{xz} = c_{xz} + 1/u_{xy}$
$c_{zx} = c_{zx} + 1/u_{yz}$

$c_{yz} = c_{yz} + 1/u_{xy}$
$c_{zy} = c_{zy} + 1/u_{xz}$

We can also block the triplet algorithm to obtain better cache locality.
Instead of a single triplet of points, we consider three blocks $\mathcal{X},\mathcal{Y},\mathcal{Z}$ and all triplets $(x,y,z)\in \mathcal{X} \times \mathcal{Y} \times \mathcal{Z}$.
We obtain locality on cache blocks of all three matrices: distance, local focus, and cohesion.
Note that a first pass is required to compute the local focus matrix in its entirety, and then blocks of the local focus matrix are read from slow memory during the second pass as needed.

Dependency structure of the blocked triplet algorithm. The highlighted regions represent entries with temporal locality. Matrices in red correspond to reads and ones in green correspond to writes. Matrices in orange correspond to writes during the first pass and reads during the second pass. Blue represents the entry-wise dependencies within each matrix.

<Ref> illustrates the dependencies among the distance, local focus, and cohesion matrices for the (blocked) triplet algorithm.
In the first pass, the blocked triplet algorithm reads $3$ blocks from the distance matrix, corresponding to the triplet pairs: $(x,y), (x,z), (y,z)$, and writes to the corresponding $3$ blocks of the local focus matrix.
Note that the distance and local focus matrices are symmetric so only the upper triangular parts are required.
The cohesion matrix is not symmetric, thus in the second pass $6$ blocks must be updated by performing distance comparisons (by reading $D_{\mathcal{X},\mathcal{Y}}, D_{\mathcal{X},\mathcal{Z}}, D_{\mathcal{Y},\mathcal{Z}}$) and utilizing entries of the local focus matrix (by reading $U_{\mathcal{X},\mathcal{Y}}, U_{\mathcal{X},\mathcal{Z}}, U_{\mathcal{Y},\mathcal{Z}}$).

§ SEQUENTIAL ALGORITHM ANALYSIS.

We model performance using the model, $\gamma F + \beta W$, where $F$ and $W$ represent an algorithm's computation and bandwidth costs, respectively, and $\gamma$ (time per operation) and $\beta$ (time per word moved) represent hardware parameters.
We analyze communication cost assuming a two-level memory hierarchy, which contains fast memory (cache) of size $M$ words and slow memory (DRAM) with unbounded size.
We assume that computation can only be performed on operands residing in fast memory.
If operands are in slow memory, then they must first be read into fast memory.
We limit analysis in this section to a two-level memory hierarchy, but this memory model can be used to analyze communication for each adjacent pair of levels in a multi-level memory hierarchy.

§.§ Communication Lower Bounds.

We use the framework in [1] to derive communication lower bounds.
The lower bound of <cit.> applies to all three-nested-loops (3NL) computations as defined in that paper.
We reproduce the 3NL definition here using the same notation, with sets $S_a, S_b, S_c \subseteq [n] \times [n]$ where $[n] = \{1,2,\dots,n\}$ and mappings $\mathbf{a}:S_a \rightarrow \cal{M}$, $\mathbf{b}:S_b \rightarrow \cal{M}$, $\mathbf{c}:S_c \rightarrow \cal{M}$, where $\cal{M}$ is slow memory.
For each $(i,j)\in S_c$, we also have a set $S_{ij} \subseteq [n]$.

A computation is considered to be three-nested-loops (3NL) if it includes computing, for all $(i,j)\in S_c$ with $S_{ij}$,

Mem(𝐜(i,j)) = f_ij({g_ijk(Mem(𝐚(i,k)),Mem(𝐛(k,j))}_k∈S_ij),

(a) mappings $\mathbf{a}$, $\mathbf{b}$, $\mathbf{c}$ are all one-to-one into slow memory, and (b) functions $f_{ij}$ and $g_{ijk}$ depend nontrivially on their arguments.

We first verify that the cohesion matrix computation defined by <ref> is 3NL when the distance matrix is stored explicitly in memory.
To satisfy the first constraint, we define the mappings $\mathbf{a}$, $\mathbf{b}$, and $\mathbf{c}$ as all mapping onto the distance matrix (that is, each mapping is one-to-one but the three mappings are not disjoint).
Here $\mathbf{a}(x,y)$ maps to the distance matrix entry $d_{xy}$.
To satisfy the second constraint, we see that computing $g_{xyz}$ depends nontrivially on $\mathbf{a}(x,y)$ and $\mathbf{b}(y,z)$, as both values must be compared with $d_{xz}$ to evaluate the indicator functions, and computing $c_{xz}$ depends nontrivially on its arguments, as it computes the sum over all values.
As argued in <ref>, the number of 3NL operations is $\sum_{i,j} |S_{ij}| = O(n^3)$.
Then, by <cit.>, the bandwidth cost lower bound for PaLD is $W=\Omega(n^3/\sqrt M)$.

§.§ Cost Analysis.

The blocked algorithms are described in <ref> with memory reference patterns depicted in <ref>.
The loop structures of the blocked algorithms are shown (with OpenMP parallelization) in <ref>.
We focus on the sequential costs in this section and discuss parallelization in <ref>.
Since the algorithms require mixed comparison and arithmetic instructions, we explicitly define the hardware parameters $\gamma_{cmp}$ and $\gamma_{fma}$ to represent the time per instruction for floating-point comparisons and FMAs, respectively.
We ignore the cost of integer arithmetic.
<Ref> shows the loop structure of the blocked pairwise algorithm where inner loop computations match <ref>.
We use $b$ to represent the block size for the pairwise algorithm.

The blocked pairwise algorithm has the leading order computation and communication costs:

\begin{align*}
F &= (5\gamma_{cmp} + 1\gamma_{fma})\cdot n\binom{n}{2} \approx 3n^3~\text{flops.}\\
W &= 4\sqrt{2}~\frac{n^3}{\sqrt{M}} \approx 5.7~\frac{n^3}{\sqrt{M}}~\text{words moved.}
\end{align*}

The blocked pairwise algorithm selects $\binom{n/b+1}{2}$ unique sets of points $\mathcal{X}, \mathcal{Y}$ with $|\mathcal{X}| = |\mathcal{Y}| = b$.
A total of $nb^2$ iterations are required to determine if a third point, $z$, is in the local focus for each $(x,y) \in \mathcal{X} \times \mathcal{Y}$.
The local focus update requires $2$ floating-point comparisons followed by $1$ integer accumulate into $u_{xy}$.
The cohesion update requires $3$ floating-point comparisons and $1$ FMA, as the reciprocals of elements of $U_{\mathcal{X},\mathcal{Y}}$ can be pre-computed once.
When $\mathcal{X} = \mathcal{Y}$, only $n\binom{b}{2}$ iterations are required to perform local focus and cohesion updates.
There are $n/b$ such overlapping sets.
Multiplying over the iterations, summing the work over the local focus and cohesion update loops, and multiplying by $\gamma_{cmp}$ and $\gamma_{fma}$ yields the computation cost.

Each of the $\binom{n/b + 1}{2}$ possible combinations of $\mathcal{X} \times \mathcal{Y}$ points requires reading the $b\times b$ block $D_{\mathcal{X},\mathcal{Y}}$ from slow memory.
In the first pass to compute the local focus sizes, for each third point, $z$, we read the two $b \times 1$ vectors $D_{\mathcal{X},z}$ and $D_{\mathcal{Y},z}$ from slow memory.
The local focus block $U_{\mathcal{X}, \mathcal{Y}}$ is computed and remains resident in fast memory.
Similarly, each iteration of the second pass cohesion update requires reading the $b \times 1$ vectors $D_{\mathcal{X},z}, D_{\mathcal{Y},z}, C_{\mathcal{X},z}$ and $C_{\mathcal{Y},z}$ from slow memory.
After each iteration within the second pass, $C_{\mathcal{X},z}$ and $C_{\mathcal{Y},z}$ must be written to slow memory.
We must maintain $2 b^2$ words of data in fast memory for $D_{\mathcal{X},\mathcal{Y}}$ and $U_{\mathcal{X},\mathcal{Y}}$, along with a constant number of length-$b$ vectors, so $b \leq \sqrt{{M}/2}$ to leading order.
Multiplying and summing these reads and writes over all iterations yields the leading order communication cost $4 n^3/b$, and choosing $b \approx \sqrt{{M}/2}$ yields the result.

<Ref> shows the loop structure of the blocked triplet algorithm, and the inner loop computations match <ref>.
The local focus sizes and cohesion matrix updates are computed in two separate passes, and two block sizes $\hat b$ and $\tilde b$ can be tuned independently.

The blocked triplet algorithm has the leading order computation and communication costs:

\begin{align*}
F &= (6\gamma_{cmp} + 2\gamma_{fma})\cdot \binom{n}{3} \approx 1.33n^3~\text{flops.}\\
W &= \left(\sqrt{6} + 4\sqrt{3}\right)\frac{n^3}{\sqrt{M}}\approx 9.4 \frac{n^3}{\sqrt{M}}~\text{words moved.}
\end{align*}

The blocked local focus and cohesion matrix passes have the same loop structure, each selecting $\binom{n/b + 2}{3}$ triplets of sets $\mathcal{X}, \mathcal{Y},$ and $\mathcal{Z}$ each of size $b$ points, though the value of $b$ differs in the two passes.
The triplet algorithm contains $3$ types of symmetry: $\mathcal{X} = \mathcal{Y} = \mathcal{Z}$, $\mathcal{X} \neq \mathcal{Y} = \mathcal{Z}$, and $\mathcal{X} = \mathcal{Y} \neq \mathcal{Z}$.
While our implementation accounts for each type of symmetry, we ignore it in our leading order cost analysis.
The local focus and cohesion update inner iterations each require $3$ distance comparisons to determine the pair of points with minimum distance.
The cohesion update iteration additionally requires $2$ FMAs to update entries of the cohesion matrix.
Multiplying operations by their respective $\gamma$ terms and summing work over the two passes proves the computation cost.

There are $\binom{n/\hat{b} + 2}{3}$ possible combinations of triplet blocks in the local focus pass.
The local focus update must read 2 $\hat{b} \times \hat{b}$ blocks of $D$, read 2 $\hat{b}\times \hat{b}$ blocks of $U$, and write 2 $\hat{b}\times \hat{b}$ blocks of $U$ from/to slow memory.
Note that the block $D_{\mathcal{X},\mathcal{Y}}$ can be read and the block $U_{\mathcal{X},\mathcal{Y}}$ read and written only $\binom{n/\hat{b} + 1}{2}$ times since they remain fixed while blocks $\mathcal{Z}$ vary in the innermost loop.
The cohesion update requires reading 2 $\tilde{b} \times \tilde{b}$ blocks of $D$ and $U$, respectively, followed by reading and writing 4 $\tilde{b} \times \tilde{b}$ blocks of $C$.
The blocks $D_{\mathcal{X},\mathcal{Y}}$ and $U_{\mathcal{X},\mathcal{Y}}$ are read from slow memory and the blocks $C_{\mathcal{X},\mathcal{Y}}$ and $C_{\mathcal{Y},\mathcal{X}}$ can be read and written $\binom{n/\tilde{b} + 1}{2}$ times.
The total I/O cost is then $n^3/\hat b + 2 n^3/\tilde b$, assuming that all blocks can be stored in fast memory.
This requires that $\hat b \leq \sqrt{{M}/6}$ and $\tilde b \leq \sqrt{{M}/12}$ to leading order.
Choosing block sizes at their approximate maximum value yields the communication cost.

The constants for the communication cost in <Ref> can be improved by unblocking the innermost loop over $\mathcal{Z}$ for the local focus and cohesion update passes, which allows for a slightly larger block size.
We use this technique for the pairwise algorithm, and it is useful in practice for matrix multiplication as well [18].
However, incorporating this optimization did not allow for auto-vectorization during cohesion updates where some updates require a stride of $n$.
Blocking all three loops allowed for unit-stride for all cohesion updates. We provide more details in the following section.

We can conclude from <Ref> that the pairwise variant requires more computation than the triplet variant, but it moves less data.
Both sequential variants attain the 3NL lower bound of $\Omega({n^3}/{\sqrt{M}})$ and are communication-optimal within a constant factor.
We will show in the next section how additional performance optimizations can yield large speedups.
The optimized sequential algorithms serve as the baselines from which we derive efficient shared-memory parallel algorithms.

§ SEQUENTIAL PERFORMANCE OPTIMIZATION.

We study the performance improvements achieved by each optimization, the tuning parameters introduced, and performance tradeoffs between the pairwise and triplet variants.
All algorithms were written in C and compiled with the Intel C compiler () release 2021.06.
The code was compiled with the following compiler flags: .
Experiments are performed on a single-node, dual-socket platform with two Intel Xeon Gold 6226R CPUs (16 cores per socket).
We run 5 trials for each experiment and use the mean to compute speedups.
We observe low runtime variance across trials, so we omit error bars for simplicity.
We perform experiments on randomly generated distance matrices for powers of two $n \in \{128, \ldots,4096\}$.
Our code can handle arbitrary square matrix sizes, but we limit performance evaluation to powers of two.

Speedup achieved from various performance optimizations applied to the Pairwise and Triplet algorithms.
Speedups are arranged by algorithm and relative to the previous performance optimization attempted.
The naive implementations of pairwise and triplet have a speedup of $1$.

We begin performance tuning by applying one level of blocking to <ref> (naive pairwise) and <ref> (naive triplet).
We show speedups relative to the previous optimization tried in <ref> with a fixed $n = 2048$ matrix.
Overall speedup over naive pairwise (resp. naive triplet) may be obtained by multiplying speedups across all optimizations.
Naive triplet resulted in a speedup of $1.11\times$ over naive pairwise due to less computation.
Introducing one level of blocking to naive pairwise led to a speedup of $1.07\times$.
Applying blocking to the triplet variant led to speedups of $1.20\times$ over naive triplet ($1.33\times$ over naive pairwise).
<Ref> require branches to correctly update $U$ and $C$ based on distance comparisons.
Distance comparisons can be vectorized, but updates to $U$ and $C$ cannot due to branching.
We avoid branches in both algorithms by computing auxiliary mask variables and performing FMAs with these explicit masks.
For <ref>, we compute the masks: $r = d_{xz} < d_{xy}~\Or d_{yz} < d_{xy}$ and $s = d_{xz} < d_{yz}$.
The variable $r$ indicates that $z$ is in the $(x,y)$ local focus and $s$ determines the entry of $C$ to update.
$C$ can be updated via two FMAs: $c_{xz} = c_{xz} + r \cdot s \cdot (1/u_{xy})$ and $c_{yz} = c_{yz} + (r)(1 - s)(1/u_{xy})$.
Branch avoidance introduces a performance tradeoff by increasing computation (e.g. performing FMAs with explicit zeros) but eliminates branch misprediction overhead.
For <ref>, branch avoidance enables a fixed stride length for updates of $C$ and facilitates other compiler optimizations (e.g. auto-vectorization and loop unrolling).
Branch avoidance alone yielded a speedup of $1.7\times$ over naive pairwise.
While branch avoidance allows for vectorization, updates to $c_{xz}$ and $c_{yz}$ require a stride length of $n$.
After blocking, we reduce the stride length to $1$ by updating columns of $C$ instead (see <ref>).
The combination achieved speedups of $20.2\times$ over naive pairwise.

<Ref> must determine the closest pair of points from a triplet $(x,y,z)$.
We avoid branches in <ref> by computing three masks from three floating point comparisons: $r = d_{xy} < d_{xz}~\tAnd d_{xy} < d_{yz}$, $s = (1 - r)(d_{xz} < d_{yz})$, and $t = (1 - r)(1 - s)$.
$C$ can then be updated using six FMAs:

\begin{align*}
c_{xy} = c_{xy} + r\left(1/u_{xz}\right), \quad &
c_{yx} = c_{yx} + r\left(1/u_{yz}\right),\\
c_{xz} = c_{xz} + s\left(1/u_{xy}\right), \quad &
c_{zx} = c_{zx} + s\left(1/u_{yz}\right),\\
c_{yz} = c_{yz} + t\left(1/u_{xy}\right), \quad &
c_{zy} = c_{zy} + t\left(1/u_{xz}\right).
\end{align*}

Applying branch avoidance to the triplet algorithm yields a speedup of $0.98\times$ due to the stride-$n$ updates to $C$.
When combined with blocking, however, we attain speedups of $20\times$ over naive triplet.
Triplet with branch avoidance and blocking yields a speedup of $1.1\times$ over pairwise with the same optimizations.
We were able to extract additional speedup by replacing floating point operations with integer operations during local focus updates, and ignoring equality in pairwise/triplet distance comparisons.
Each entry of $U$ counts the number of points in the local focus based on distance comparisons, with results stored in a mask register.
If $U$ is stored as a floating point array, then each increment to update $U$ requires an expensive integer mask to 32-bit floating point cast operation.
We avoid this by storing $U$ as an integer array during the local focus computation.
This allowed us to combine casting with computing reciprocals prior to cohesion updates.

The theoretical formulation of PaLD [2] allows for ties in pairwise distances (e.g., $d_{xz} == d_{yz}$).
When ties occur, support is split between cohesion entries $c_{xz}$ and $c_{yz}$ (i.e. $c_{xz} = c_{xz} + r\cdot s \cdot \left(0.5/u_{xy}\right)$).
In finite arithmetic, floating point equality is unlikely due to round-off and truncation.
Avoiding ties is critical for <ref> which contains more distance tie permutations than pairwise.
Introducing these additional optimizations yields self-relative speedups (over naive) of $25.5\times$ and $26.2\times$ for pairwise and triplet, respectively.
Overall, optimized triplet achieves a speedup of $1.14\times$ over optimized pairwise for $n = 2048$.

Speedup achieved from block size tuning for pairwise (top) and triplet (bottom) for $n = 2048$.

We also perform block size tuning for each algorithm.
We experiment with (powers of two) block sizes in the range $[2^5, 2^{10}]$. Optimized pairwise attains a maximum speedup of $25.5\times$ for $n = 2048$ after tuning.

For optimized triplet, updates to $U$ require storing $3$ distinct blocks of $D$ and $3$ distinct blocks of $U$ in cache.
Updates to $C$ require $3$ distinct blocks of $D$, $3$ distinct blocks of $U$, and $6$ distinct blocks of $C$ in cache.
This suggests that different block sizes may be better than a fixed block size.
<Ref> (bottom) illustrates the speedups observed (over <ref>) for various block size combinations for the optimized triplet algorithm.
We observe a maximum speedup of $26.2\times$ over naive triplet with $\hat{b} = 256$ and $\tilde{b} = 128$.

$n$     Pairwise Optimized     Triplet Optimized
128     0.00117 (1.58$\times$)     0.00185
256      0.00497 (1.34$\times$)     0.00665
512      0.0188 (1.18$\times$)     0.0221
1024     0.1274     0.1208 (1.05$\times$)
2048     0.9942     0.8734  (1.14$\times$)
4096     8.3623     6.6111  (1.26$\times$)

Running time in seconds (and speedup) comparison of  pairwise and triplet algorithms.

In <ref> we compare running times (and speedups) of optimized pairwise and optimized triplet over a range of input matrix sizes.
For small matrix sizes, where $D,U$ and $C$ all fit in cache, optimized pairwise is fastest (e.g. speedup of $1.58\times$ over triplet at $n = 128$).
This is because $n/b$ is a small integer where lower order terms dominate (see <ref>).
For larger matrices, optimized triplet performs better (speedup of $1.26\times$ over pairwise at $n = 4096$) due to lower computation cost.
In practice, we expect triplet to be the better sequential variant for most applications of PaLD.
If distances ties must be handled correctly, then pairwise is the better variant due to fewer branches.

Finally, we note that optimized pairwise attains $27.7\%$ of hardware peak at $n = 2048$ and optimized triplet attains $28\%$ at $n = 8192$.
Our Intel CPU has a single-core, single precision peak of $249.6$ Gflops/sec.
Single precision comparisons on our CPU have a cycles-per-instruction (CPI) of $1$ while all other single precision ops have a CPI of $0.5$.
Thus, floating point comparisons are twice as expensive.
See <Ref> for details on percentage of peak calculations for each algorithm.

The combination of all optimizations achieves speedups of $25.5\times$ and $29\times$ for pairwise and triplet, respectively, over naive pairwise (for $n = 2048$).
We observe speedups of $23\times$ and $26.2\times$ over naive triplet.

§ SHARED-MEMORY PARALLEL ALGORITHMS

This section presents the OpenMP parallelization of the optimized sequential pairwise and triplet algorithms.
<Ref> shows the OpenMP version of the blocked pairwise algorithm.
The blocked pairwise algorithm first computes $U_{\mathcal{X},\mathcal{Y}}$ with a pass over all $n$ points $z$.
The local focus $z$-loop can be parallelized across $p$ threads using the OpenMP  construct.
All threads must write to $U_{\mathcal{X},\mathcal{Y}}$ so a sum-reduction is required to resolve write conflicts.
The cohesion update pass requires the quantities $1/u_{xy}~\forall~(x,y)~\in~\mathcal{X}{\times}\mathcal{Y}$, which can be parallelized without write conflicts.
Cohesion updates are within each column of $C$ to entries of $C_{\mathcal{X},z}$ and $C_{\mathcal{Y},z}$.
The cohesion pass can be parallelized without write conflicts by splitting the $z$-loop across $p$ threads.
<Ref> illustrates the write patterns for optimized OpenMP pairwise for $n = 16$, $b = 4$, and $p = 8$.
Updates to entries of $C$ requires corresponding entries from $D$, so $D$ can also be partitioned column-wise.
The pairwise algorithm is amenable to NUMA optimizations due to the regular data dependencies.

<Ref> shows the OpenMP version of the blocked triplet algorithm. The triplet approach requires reading all of $D$ for local focus and cohesion update passes.
Blocking is performed over triplets of points, $\mathcal{X},\mathcal{Y},\mathcal{Z}$, and updates to $U$ and $C$ become irregular.
We use the OpenMP tasking model [17] for parallelism.
Each triplet block, $\mathcal{X}\times\mathcal{Y}\times\mathcal{Z}$, is a new task that can be executed by any available thread.
Tasks in the local focus pass write to $3$ blocks of $U$. $C$ is not symmetric, so the cohesion update pass writes to $6$ blocks.
Write conflicts arise when multiple tasks need to update the same blocks of $U$ or $C$.
We resolve conflicts by annotating dependencies using the  clause with the  modifier.
<Ref> shows the write conflicts for the local focus pass.
Each vertex represents one of the $\binom{n/b + 2}{3}$ tasks and is labeled by $\mathcal{X},\mathcal{Y},\mathcal{Z}$ block values, and edges represent conflicts.
The degree for each vertex varies based on the symmetry in the block.
This leads to irregular dependencies which we will show in <ref> are not as amenable to NUMA optimizations.

for(xb = 0; xb < n/b; ++xb)
for(yb = 0; yb <= xb; ++yb)
!#pragma omp parallel for \!
for(z = 0; z < n; ++z)
for(x = 0; x < b; ++x)
y_start = (xb==yb) ? (x+1) : 0;
for(y = y_start; y < b; ++y)
// update !$u_{xy}$.!
!#pragma omp parallel for!
for(i = 0; i < b*b; ++i)
U[i] = 1/U[i];
!#pragma omp parallel for!
for(z = 0; z < n; ++z)
for(x = 0; x < b; ++x)
y_start = (xb==jb) ? (x+1) : 0;
for(y = y_start; y < b; ++y)
// update !$c_{xz}$ and $c_{yz}$.!

Blocked OpenMP pairwise Algorithm.

Distance matrix reads and Local Focus/Cohesion writes for parallel pairwise code with $n=16$, $b =4$, and $p = 8$.
All threads have write conflicts to the $U$ block for each pair $\mathcal{X},\mathcal{Y}$ (in red), so synchronization is required via reductions.
Only one $U$ block is needed in fast memory at any given point in time. Writes to $C$ are within one column, so column blocks can be partitioned across threads without write conflicts.

!#pragma omp single!
for(xb = 0; xb < n/b; ++xb)
for(yb = xb; yb < n/b; ++yb)
for(zb = yb; xb < n/b; ++zb)
x_end=(xb==yb        yb==zb)?(b-1):b
!#pragma omp task untied depend(inout, U[$\mathcal{X},\mathcal{Y}$],U[$\mathcal{X},\mathcal{Z}$],U[$\mathcal{X},\mathcal{Z}$])!
for(x = 0; x < x_end; ++x)
y_start=(xb==yb) ? (x+1) : 0;
for(y = y_start; y < b; ++y)
z_start=(yb==zb) ? (y+1) : 0;
for(z = z_start; z < zb; ++z)
// update !$u_{xy},~u_{xz},~u_{yz}$.!
!#pragma omp parallel for!
for(i = 0; i < n*n; ++i)
U[i] = 1/U[i];

!#pragma omp single!
for(xb = 0; xb < n/b; ++xb)
for(yb = xb; yb < n/b; ++yb)
for(zb = zb; xb < n/b; ++zb)
x_end=(xb==yb        yb==zb)?(b-1):b
!#pragma omp task untied depend(inout, C[$\mathcal{X},\mathcal{Y}$],C[$\mathcal{X},\mathcal{Z}$],C[$\mathcal{Y},\mathcal{Z}$]) \!
!depend(inout, C[$\mathcal{Y},\mathcal{X}$],C[$\mathcal{Z},\mathcal{X}$],C[$\mathcal{Z},\mathcal{Y}$])!
for(x = 0; x < xend; ++x)
y_start=(xb==yb) ? (i+1) : 0;
for(y = ystart; y < b; ++y)
z_start=(yb==zb) ? (y+1) : 0;
for(z = z_start; z < zb; ++z)
// update !$c_{xy},~c_{xz},~c_{yz}$,!
// update !$c_{yx},~c_{zx},~c_{zy}$.!

Blocked OpenMP triplet Algorithm.

Task diagram for parallel triplet with $n/b=4$, where nodes are labeled by their $\mathcal{X},\mathcal{Y},\mathcal{Z}$ block values. Edges represent write conflicts for $U$ between tasks.

§.§ OpenMP Performance.

We use OpenMP version 4.5 and test the OpenMP algorithms on randomly generated dense distance matrices with $n \in \{2048, 4096, 8192\}$. We incorporate NUMA optimizations into the pairwise algorithm by controlling thread affinity via the  and  environment variables. We map OpenMP threads to physical cores, by assigning OpenMP thread ids $0$ to $16$ to CPU $0$ and threads $17$ to $31$ to CPU $1$.

OpenMP pairwise speedup from NUMA optimizations with $n \in \{2048, 4096, 8192\}$ and $p = 32$.

A static loop schedule yields best performance due to the pairwise algorithm's regular dependencies.
Each thread reads columns of $D$ and $C$ from thread-local fast memory so updates to $C$ are spatially local.
Thread binding ensures that accesses are temporally local by assigning fixed column blocks of $D$/$C$ to threads.
OpenMP allocates memory pages using a first-touch policy by default.
If a single thread allocates $D$, then $D$ resides in the memory hierarchy of the thread's CPU.
$D$ is typically computed outside the scope of the OpenMP algorithms, so we also study the effects of partitioning $D$ across sockets (i.e. memory binding).

<Ref> shows the speedup achieved by introducing thread binding only and thread + memory binding into the OpenMP pairwise algorithm across three matrix sizes, $n \in \{2048, 4096, 8192\}$.
We use the OpenMP pairwise algorithm without NUMA-aware optimizations as our baseline and report speedups for $32$ OpenMP threads.
When we use thread binding only, we observe average speedups of $1.4\times, 1.5\times,$ and $1.13\times$ for $n = 2048, 4098,$ and $8192$, respectively.
Thread binding with memory binding yields average speedups of speedup of $1.7\times, 1.69\times,$ and $1.2\times$ over the baseline.
We did not perform TLB optimizations, therefore, we observe decreasing speedups for large matrix sizes.
We also found that NUMA optimizations are useful at smaller thread counts, $2 \leq p \leq 16$, by mapping half the threads to CPU $0$ and the other half to CPU $1$.
This mapping provides access to the fast memory hierarchies on both CPUs.
We observe speedups ranging from $1.05\times$ ($n = 4096, p = 2$) to $1.33\times$ ($n = 2048, p = 16$) when splitting threads (where $p \leq 16$) across sockets.
We experimented with thread binding for the OpenMP triplet algorithm but not memory binding due to the irregular data dependencies.
However, we did not observe significant performance improvements over the baseline, so we omit these results from <ref>.
We obtain best OpenMP scaling when using the  clause, which allows suspended tasks to be resumed on any available thread.
Suspended tasks may cause additional reads from slow memory after restart.
Hence, we do not expect NUMA optimizations to be helpful.

Self-relative strong scaling efficiency of OpenMP Pairwise (top) and Triplet (bottom).

We perform strong scaling experiments in <ref> of the OpenMP variants under the same settings as for <ref> and report self-relative efficiency achieved.
We report efficiencies with and without NUMA optimizations.
The pairwise algorithm without NUMA optimizations achieves efficiencies of $24.2\%, 33.5\%,$ and $50.6\%$ at $p = 32$ for $n = 2048, 4096$ and $8192$, respectively.
Including NUMA optimizations yields efficiencies of $42.9\%, 56.6\%,$ and $60.5\%$ for $p = 32$.
The triplet algorithm achieves efficiencies of $28.0\%, 29.2\%,$ and  $40.9\%$ without NUMA optimizations and $36.9\%, 34.9\%$, and $41.2\%$ with NUMA optimizations for $p = 32$.
The triplet algorithm is the faster sequential baseline, hence the OpenMP triplet efficiencies are lower than those reported for OpenMP pairwise.
We also study weak scaling of the two algorithms with and without NUMA optimizations.
We fix $n^3/p$ over the range of $p$ tested.
We use the matrix sizes $n_1 \in \{2048, 4096, 8192\}$, where $n_1$ is the matrix size at $p = 1$.
<Ref> shows the results of the weak scaling experiments.
The pairwise algorithm without NUMA optimizations attains weak scaling efficiencies of $30.6\%, 48.2\%$, and $61.4\%$ for $n_1 = 2048, 4096,$ and $8192$, respectively at $32$ threads.
With NUMA optimizations, the efficiencies increase to $59.1\%,~ 63.6\%$, and $65.6\%$ for each of the matrix size settings at $p = 32$.
Triplet without NUMA optimizations achieves weak scaling efficiencies of $44.2\%, 49.1\%,$ and $50.1\%$ and $47.6\%, 49.1\%,$ and $50.1\%$ with NUMA optimizations at $p = 32$.

Self-relative weak scaling efficiency of OpenMP Pairwise (top) and Triplet (bottom).

§ TEXT ANALYSIS APPLICATION

Word clouds from PaLD analysis (left column) and distance analysis (right column) of the words guilt and halt. Font size is proportional to cohesion values and inverse distances.

We demonstrate the utility of PaLD on larger datasets than previously considered [2] for semantic analysis of words extracted from Shakespeare sonnets [11].
Words are converted to vectors using the pre-trained fastText word embedding [4, 12], yielding a dataset of 2712 words.
We compute Euclidean distance between embedding vectors and generate the cohesion matrix $C$ using the OpenMP pairwise algorithm.
<Ref> shows words associated with guilt and halt obtained from PaLD and from analyzing only the distance matrix $D$.
PaLD is parameter-free, with strong ties determined by a universal threshold (see [2]), whereas analysis using $D$ requires a user-tuned distance or neighbor-count cutoff.
Note the differing sizes of strong-tie neighborhoods between the two words.
PaLD finds $20$ words with strong ties to guilt and 5 words for halt.
The $20$ closest words to guilt based on distance correspond to a cutoff of $2.26$.
We observe significant overlap between the two sets, though PaLD reports stronger ties to  expiate and conscience.
PaLD finds $5$ words with strong ties to halt.
To illustrate the pitfalls of tuning an absolute distance threshold, we apply the distance cutoff $2.26$ for halt, which yields $23$ words including several unrelated ones (e.g. just and say).
This suggests that absolute distance thresholds are not robust to varying density and distance scales within word neighborhoods.
A distance cutoff of $2.14$ is required for halt to match results obtained from PaLD.
Applying the cutoff to guilt identifies only $8$ related words, missing several words like expiate.
We attain a speedup of $16.7\times$ using the NUMA optimized OpenMP pairwise algorithm at $p = 32$ and an overall run time of $0.178$ seconds.

§ CONCLUSION

This paper presents several sequential and shared-memory parallel algorithms for PaLD [2].
We prove that sequential variants are communication-optimal, up to constant factors.
We illustrate that branch avoidance is critical to attaining high performance; achieving a speedup of up to $29\times$ over naive sequential variants.
Based on our theoretical and empirical studies, we conclude that the triplet variant is the faster sequential algorithm for large matrices due to less computation.
However, we show that the pairwise algorithm is more amenable to parallelization due to regular data dependencies and load balance.
We observe strong scaling speedups  up to $19.4\times$ ($60.5\%$ efficiency), and weak scaling efficiencies of up to $65.6\%$ at $p = 32$ after incorporating NUMA-aware optimizations.
With the performance achieved on the text analysis application, we show that PaLD can be scaled to nearly any dataset with a distance matrix that fits in the memory of a single server.

Acknowledgements. We would like to thank Kenneth S. Berenhaut for helpful feedback on the presentation of PaLD and discussions on applying PaLD to semantic analysis of word embedding in <Ref>. We would also like to thank Yixin Zhang for code contributions to preliminary versions of the pairwise algorithms. This work is supported by the National Science Foundation under Grant No. OAC-2106920 and the U.S. Department of Energy, Office of Science, Advanced Scientific Computing Research program under Award Number DE-SC-0023296.

[1]
G. Ballard, E. Carson, J. Demmel, M. Hoemmen, N. Knight, and O. Schwartz,
Communication lower bounds and optimal algorithms for numerical linear
algebra, Acta Numerica, 23 (2014), p. 1–155,

[2]
K. S. Berenhaut, K. E. Moore, and R. L. Melvin, A social perspective
on perceived distances reveals deep community structure, Proceedings of the
National Academy of Sciences, 119 (2022), p. e2003634119,

[3]
J. Bilmes, K. Asanovic, C.-W. Chin, and J. Demmel, Optimizing matrix
multiply using phipac: A portable, high-performance, ansi c coding
methodology, in ACM International Conference on Supercomputing, New York,
NY, USA, 1997, Association for Computing Machinery, p. 253–260,

[4]
P. Bojanowski, E. Grave, A. Joulin, and T. Mikolov, Enriching word
vectors with subword information, CoRR, abs/1607.04606 (2016),

[5]
R. J. G. B. Campello, D. Moulavi, and J. Sander, Density-based
clustering based on hierarchical density estimates, in Advances in Knowledge
Discovery and Data Mining, J. Pei, V. S. Tseng, L. Cao, H. Motoda, and G. Xu,
eds., Berlin, Heidelberg, 2013, Springer Berlin Heidelberg, pp. 160–172.

[6]
R. J. G. B. Campello, D. Moulavi, A. Zimek, and J. Sander, Hierarchical density estimates for data clustering, visualization, and
outlier detection, ACM Trans. Knowl. Discov. Data, 10 (2015),
<https://doi.org/10.1145/2733381>, <https://doi.org/10.1145/2733381>.

[7]
M. Ester, H.-P. Kriegel, J. Sander, and X. Xu, A density-based
algorithm for discovering clusters in large spatial databases with noise, in
Proceedings of the Second International Conference on Knowledge Discovery and
Data Mining, KDD'96, AAAI Press, 1996, p. 226–231.

[8]
J. Goldberger, G. E. Hinton, S. Roweis, and R. R. Salakhutdinov, Neighbourhood components analysis, in Advances in Neural Information
Processing Systems, L. Saul, Y. Weiss, and L. Bottou, eds., vol. 17, MIT
Press, 2004,

[9]
K. Goto and R. A. v. d. Geijn, Anatomy of high-performance matrix
multiplication, ACM Trans. Math. Softw., 34 (2008),

[10]
J. L. Hennessy and D. A. Patterson, Computer Architecture, Sixth
Edition: A Quantitative Approach, Morgan Kaufmann Publishers Inc., San
Francisco, CA, USA, 6th ed., 2017.

[11]
T. M. Inc., Text analytics toolbox version: 1.9 (r2022b), 2022,

[12]
A. Joulin, E. Grave, P. Bojanowski, and T. Mikolov, Bag of tricks
for efficient text classification, CoRR, abs/1607.01759 (2016),

[13]
N. P. Jouppi, Superscalar vs. superpipelined machines, SIGARCH
Comput. Archit. News, 16 (1988), pp. 71–80,

[14]
N. P. Jouppi and D. W. Wall, Available instruction-level parallelism
for superscalar and superpipelined machines, in ASPLOS-III Proceedings -
Third International Conference on Architectural Support for Programming
Languages and Operating Systems, Boston, Massachusetts, USA, April 3-6, 1989,
J. S. Emer and J. L. Hennessy, eds., ACM Press, 1989, pp. 272–282,

[15]
J. Leskovec, J. Kleinberg, and C. Faloutsos, Graph evolution:
Densification and shrinking diameters, ACM Trans. Knowl. Discov. Data, 1
(2007), p. 2–es, <https://doi.org/10.1145/1217299.1217301>,

[16]
J. Leskovec and A. Krevl, SNAP Datasets: Stanford large network
dataset collection.
<http://snap.stanford.edu/data>, June 2014.

[17]
OpenMP Architecture Review Board, OpenMP application program
interface version 5.2, Nov. 2021,

[18]
T. M. Smith, B. Lowery, J. Langou, and R. A. van de Geijn, A tight
I/O lower bound for matrix multiplication, Tech. Report 1702.02017, arXiv,
2019, <https://arxiv.org/abs/1702.02017>.

§ PERCENTAGE OF HARDWARE PEAK

§.§ Pairwise Algorithm.

This section details the operations counting for the optimized sequential pairwise and triplet algorithms used for percentage of peak calculations in <Ref>.
Unless otherwise specified, all operands are assumed to be in $32$-bit floating point format.
The optimized sequential pairwise algorithm requires $2$ comparisons during local focus update to determine if a point, $z$, is in the neighborhood.
The local focus matrix is incremented based on these comparisons.
However, since $U$ is stored in integer format, we ignore the cost of integer increments during the local focus pass.
The cohesion update requires $3$ comparisons: $2$ comparisons to compute mask $r$ which determines if a point $z$ is in the local focus and $1$ comparison to compute mask $s$ which determines the column entry of $C$ to update.
Since results of floating point comparisons are stored in unsigned (integer) format, $r$ and $s$ must be cast to $32$-bit floats before FMAs.
This requires $2$ unsigned int to floating point cast operations.
Finally, $2$ FMAs (each FMA requires two instructions) can be used to update $c_{xz}$ and $c_{yz}$.
We explicitly compute both entries and accumulate with an explicit zero, as this avoids branching.
Based on these operations, the total number of operations for sequential pairwise can be computed as follows:

\begin{align*}
F = \left(5 \gamma_{cmp} + 2\cdot 2\gamma_{fma} + 2\gamma_{cast}\right)\cdot n \binom{n}{2}
\end{align*}

On our Intel Xeon Gold 6226R CPU, floating point comparisons have a CPI of $1$ whereas FMAs and casting each have a CPI of $0.5$.
Since comparisons are twice as expensive, we normalize our operation count to be relative to FMA/cast.
After normalization, the total number of operations become:

\begin{align*}
F = 16\gamma\cdot n \binom{n}{2} \approx 8 n^3~\text{ops}.
\end{align*}

Finally, percentage of peak can be calculated by:

\begin{align}\label{eq:pct-peak}
\frac{1}{249.6}\cdot \frac{F}{10^9 \cdot t_n}
\end{align}

where $t_n$ is the runtime time (in seconds) obtained empirically from executing the optimized sequential pairwise algorithm on a matrix of size $n$ and $249.6$ Gflops/sec is the single precision, single core machine peak of our Intel CPU.
The setting $n = 2048$ and $t_n = 0.99422~\text{seconds}$ (averaged over $5$ trials) yields $27.7\%$ of peak, as reported in <Ref>.

§.§ Triplet Algorithm.

The optimized sequential triplet algorithm makes two passes: one to compute $U$ in its entirety and one to compute $C$.
The triplet algorithm, which ignores ties, requires $6$ comparisons across the two passes to uniquely determine the pair of points in a triplet with minimum pairwise distance.
The local focus pass and cohesion pass must compute these distances.
Once again, we ignore integer increments in the local focus pass.
The remaining instructions are $3$ casting and $6$ FMA operations to update entries of $C$.

\begin{align*}
F = \left(12 \gamma_{cmp} + 2\cdot 6\gamma_{fma} + 3\gamma_{cast}\right)\cdot& \binom{n}{3} \approx 6.5 n^3.
\end{align*}

Setting $F = 6.5n^3, n = 8192$ and $t_n = 51.15952$ seconds in (<ref>) yields $28\%$ of peak, as reported in <Ref>.

§ RUNTIME BREAKDOWN

Running time breakdown (as a fraction of total running time) of pairwise and triplet algorithms for a random dense, $n = 2048$ matrix.

<Ref> shows the running time breakdown of the pairwise and triplet algorithms, grouped by the algorithm.
We report the fraction of time taken to compute the local focus, cohesion update, and memory overhead (e.g. memcpy into explicit cache blocks).
The OpenMP pairwise algorithm requires a reduction during the local focus computation.
As $p$ increases, we see that the local focus computation becomes a barrier to scalability.
The OpenMP triplet algorithm however does not require explicit synchronization.
This results in better scalability for the local focus computation as $p$ increases.
The cohesion update, however, does not scale efficiently since it updates up to $6$ unique blocks of $C$.
Given the irregular task dependencies (see <ref>), the OpenMP triplet algorithm was not amenable to NUMA optimizations.

In addition, variance in task costs due to symmetries and a non-static task schedule were barriers to efficiently scaling the cohesion update pass.
Memory overhead is a negligible fraction of runtime for pairwise and triplet algorithms and does not increase with $p$.
At $32$ threads, the pairwise algorithm is faster than triplet algorithm for cohesion update.
The reverse is true for the local focus update.
This behavior may indicate that the two algorithms can be combined by utilizing the triplet approach for local focus update and the pairwise approach for cohesion update for additional speedup.

§ SCALING ON SNAP DATASETS.

Dataset     $n$     sequential     $p = 32$
ca-GrQc     5242     21.69    1.390 (15.6$\times$)
ca-HepPh     12008     259.9     13.16 (19.7$\times$)
ca-CondMat     23133     1913    91.89 (20.8$\times$)

Pairwise runtimes (in sec.) and maximum speedup over pairwise sequential on SNAP datasets.

We perform scaling experiments on large datasets obtained from the SNAP data repository [16] to illustrate PaLD scalability on collaboration networks.
We obtain distance matrices by computing all-pairs shortest path distances.
<Ref> reports the running times (in seconds) and speedup achieved at $p = 32$ for the pairwise algorithm.
We use the optimized sequential pairwise algorithm as our baseline.
We achieve speedups of $15.6\times$, $19.7\times$ and $20.8\times$ on the ca-GrQC, ca-HepPh, and ca-CondMat datasets [15], respectively.
For the largest dataset, ca-CondMat, we are able to reduce the running time of computing $C$ from $31$ minutes (optimized pairwise sequential) to $92$ seconds (OpenMP pairwise with $p = 32$).