# Distributed Discrete Morse Sandwich: Efficient Computation of Persistence Diagrams for Massive Scalar Data

Eve Le Guillou, Pierre Fortin, Julien Tierny

► To cite this version:

Eve Le Guillou, Pierre Fortin, Julien Tierny. Distributed Discrete Morse Sandwich: Efficient Computation of Persistence Diagrams for Massive Scalar Data. IEEE Transactions on Parallel and Distributed Systems, 2025, pp.1-18. 10.1109/TPDS.2025.3626047. [hal-05337555](https://HAL.science/hal-05337555v1)

HAL Id: hal-05337555

https://hal.science/hal-05337555v1

Submitted on 29 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Distributed Discrete Morse Sandwich: Efficient Computation of Persistence Diagrams for Massive Scalar Data

E. Le Guillou, P. Fortin, J. Tierny

**Abstract**—The persistence diagram, which describes the topological features of a dataset, is a key descriptor in Topological Data Analysis. The “Discrete Morse Sandwich” (DMS) method has been reported to be the most efficient algorithm for computing persistence diagrams of 3D scalar fields on a single node, using shared-memory parallelism. In this work, we extend DMS to distributed-memory parallelism for the efficient and scalable computation of persistence diagrams for massive datasets across multiple compute nodes. On the one hand, we can leverage the embarrassingly parallel procedure of the first and most time-consuming step of DMS (namely the discrete gradient computation). On the other hand, the efficient distributed computations of the subsequent DMS steps are much more challenging. To address this, we have extensively revised the DMS routines by contributing a new self-correcting distributed pairing algorithm, redesigning key data structures and introducing computation tokens to coordinate distributed computations. We have also introduced a dedicated communication thread to overlap communication and computation. Detailed performance analyses show the scalability of our hybrid MPI+thread approach for strong and weak scaling using up to 16 nodes of 32 cores (512 cores total). Our algorithm outperforms *DIPHA*, a reference method for the distributed computation of persistence diagrams, with an average speedup of ×8 on 512 cores. We show the practical capabilities of our approach by computing the persistence diagram of a public 3D scalar field of 6 billion vertices in 174 seconds on 512 cores. Finally, we provide a usage example of our open-source implementation at https://github.com/eve-le-guillou/DDMS-example.

**Index Terms**—Topological data analysis, high-performance computing, distributed-memory algorithms, persistent homology.

## I. INTRODUCTION

Modern datasets continue to grow, through constant improvements of data acquisition and numerical computations. Hence grows the need for efficient tools for the representation and analysis of those datasets. In the case of scalar data, *Topological Data Analysis* (TDA) has proven to provide practical solutions for capturing and extracting structural characteristics in high-resolution data. Topological methods have been used for data analysis in several domains, such as combustion [1]–[3], fluid dynamics [4], [5], nuclear energy [6], data science [7], [8], chemistry [9]–[11] or astrophysics [12], [13].

TDA includes a wide range of representations [14], [15], several of which have been extensively studied and efficiently

implemented, such as the persistence diagram [16], the Reeb graph [17], [18], the contour tree [19]–[21] or the Morse-Smale complex [22]–[25]. In this paper, we focus on the persistence diagram [14], a concise descriptor of the topological features of a dataset. This abstraction helps separate noise from the noteworthy structures of the data, enabling multi-scale data simplification [26], [27], segmentation [28], [29], compression [30], tracking [31]–[33] or ensemble summarization [34]–[36]. When considering *d*-dimensional scalar data, a persistence diagram for each dimension has to be computed, namely $D_0$, $D_1$ and $D_2$ for *d* = 3. Several algorithms already exist to compute these data structures, using either shared-memory or distributed-memory parallelism. In a shared-memory setting, the Discrete Morse Sandwich (DMS) algorithm has been reported to be the most efficient implementation, as shown in its associated benchmark [37]. However this setting limits the maximum size of the datasets that can be processed to the memory of a single node. This capacity is often insufficient with regard to the potentially massive size of modern datasets.

In this work, we introduce the Distributed Discrete Morse Sandwich (DDMS), an efficient algorithm for persistence diagram computation exploiting distributed-memory parallel architectures. DDMS builds upon the DMS algorithm. The distributed-memory parallelism enables DDMS to leverage greater computational power and memory capacity, allowing a user to perform computation faster than DMS and on larger datasets than ever before. The workflow of the proposed DDMS algorithm closely follows that of DMS. This enables us to benefit from the DMS efficiency within each node, by leveraging the embarrassingly parallel computation of the discrete gradient which is the first, most time-consuming, DMS step. The subsequent DMS steps are however less suitable for efficient distributed computations, in particular the pairing steps to extract $D_0$, $D_1$ and $D_2$. Like most TDA algorithms, they indeed provide a global view point on the data, which requires multiple global data traversals with little computation. This combination of factors challenges scalability due to the additional communications, synchronizations and computations required to ensure the correctness of the algorithm. For example, computing one pair of the diagram may require to access the data of multiple processes. This creates significant communication overhead. Furthermore, the pairing should happen in a specific order, creating strong data dependencies across all processes.

These subsequent steps require thus an important redesign for DDMS. Besides, DDMS relies on a hybrid MPI+thread

E. Le Guillou is with the CNRS, Sorbonne Université and University of Lille. E-mail: eve.le_guillou@sorbonne-universite.fr

J. Tierny is with the CNRS and Sorbonne Université. E-mail: julien.tierny@sorbonne-universite.fr

P. Fortin is with Univ. Lille, CNRS, Centrale Lille, UMR 9189 CRISTAL, F-59000 Lille, France. E-mail: pierre.fortin@univ-lille.fr
design combining the advantages of both shared and dis-
tributed memories: (i) larger total memory and (ii) reduced
communication and synchronization overheads, along with
improved intra-node dynamic load balancing. Even on a single
node, DDMS achieves average performance comparable to
DMS while requiring less memory, thereby rendering DDMS
advantageous for single-node execution as well. In this work,
though the original DMS was designed for both structured and
unstructured grids, we focus on structured grids.

Contributions

This paper makes the following new contributions:

1) *A hybrid distributed, shared-memory algorithm for the computation of persistence diagrams for 1D, 2D and 3D data:* Our work extends the fastest shared-memory approach (DMS) to the distributed setting, enabling an efficient computation of persistence diagrams for scalar fields of significant sizes. This is achieved thanks to the following novel procedures:

a) A novel self-correcting distributed algorithm for computing the $D_0$ and $D_2$ pairings, extending the related procedure of DMS [37] to the distributed setting. The new self-correction mechanism removes data dependencies, allows for parallel computation and reduces the need for synchronization.

b) A novel procedure for computing $D_1$, based on an efficient extension of the "PairCriticalSimplices" procedure [37]. This new procedure introduces specific data structures and computation tokens to handle the distribution of data across multiple processes and ensures the correctness of the algorithm. A dedicated communication thread is also introduced to overlap communication and computation, to reduce the idle time due to synchronizations and to improve communication reactivity.

The algorithm is output sensitive and provides substantial gains over the original DMS approach as well as *DIPHA*, the reference method for computing persistence diagrams in a distributed setting.

2) *An open-source implementation:* For reproduction purposes, we provide a C++ implementation of our approach using MPI+OpenMP and based on the Topology ToolKit (TTK) [38]–[40] integrated in ParaView.

3) *A reproducible example:* We provide a reference Python script for computing the persistence diagram with a dataset size that can be adjusted to fit the capacities of any system (publicly available at: https://github.com/eve-le-guillou/DDMS-example). The larger example on the Turbulent Channel Flow (Fig. 17) can be produced using the pipeline available in this script.

The rest of the paper is organized as follows. In Sec. II, we provide the necessary theoretical foundations, as well as a comparison to other existing methods. The original DMS algorithm is summarized in Sec. II-F. Sec. III offers a high-level overview of DDMS. In Sec. IV, we provide a detailed explanation of the distributed computation of $D_0$ and $D_2$, including our novel self-correcting procedure at the core

of this computation. In Sec. V, we present the distributed
computation of $D_1$ which has been designed by revisiting the
DMS procedure "PairCriticalSimplices". Our experiments are
presented in Sec. VI, both in strong and weak scaling. We also
demonstrate significant gain over *DIPHA* [41] (Sec. VI-E),
to our knowledge the only publicly available implementation
for this problem in a distributed context. Finally, we illustrate
the new distributed capabilities of DDMS by computing the
persistence diagram (Sec. VI-F) of a 3D scalar field of 6 billion
vertices distributed on 16 nodes of 32 cores each.

II. BACKGROUND

This section presents several theoretical elements of com-
putational topology. For a more comprehensive introduction to
this topic, we refer the reader to the reference textbook [14].

A. Input Data

The input domain $\mathcal{M}$ is a $d$-dimensional simplicial complex, with $d \le 3$ in our applications. A simplicial complex can be intuitively defined as a discretization of a topological space using small building blocks of different dimensions (a 0-simplex is a vertex, a 1-simplex an edge, a 2-simplex a triangle and a 3-simplex a tetrahedron). A face $\tau$ of a $d$-simplex $\sigma$ is a simplex defined by a non-empty subset of the $d+1$ vertices of $\sigma$. For example, the faces of a tetrahedron are its four triangles, its six edges and four vertices. Inversely, $\sigma$ is called the co-face of $\tau$. On the domain $\mathcal{M}$ is defined a piecewise linear (PL) scalar field $f: \mathcal{M} \to \mathbb{R}$. The input field $f$ is our subject of interest and is provided on the vertices of $\mathcal{M}$. $f$ is assumed to be injective on the vertices. This is enforced in practice by a procedure inspired from *Simulation of Simplicity* [42]. It consists in pre-sorting the vertices by increasing $f$ values (breaking ties according to the vertex offsets in memory) and by finally considering as input data the vertex orders in this global data sort.

When the input is provided on a regular grid, it is implic-
itly decomposed into a simplicial complex by following the
Freudenthal triangulation [43]. This method will subdivide
cubical cells into several simplices. For example, in 2D, a
square is divided in two triangles (as can be seen in Fig. 2).
This design choice is central to TTK and allows a genericity
of support which is also characterized by implementation
simplicity. Furthermore, this particular triangulation yields a
natural interpolation within the input grid which coincides with
the critical simplices identified by discrete Morse theory. The
comparison of two simplices is made by listing the vertex or-
ders of both simplices by decreasing values, and by comparing
the resulting lists with *lexicographic* comparison. This com-
parison provides a strict global order over all simplices of $\mathcal{M}$.
Let $\mathcal{M}^i$ be the union of the first $i$ simplices of $\mathcal{M}$ of this global
order. The *lexicographic filtration of* $\mathcal{M}$ is the nested sequence
of simplicial complexes $\emptyset \neq \mathcal{M}^0 \subset \mathcal{M}^1 \subset \cdots \subset \mathcal{M}^n = \mathcal{M}$
with $n$ the number of simplices of $\mathcal{M}$. Fig. 1 shows an example
of filtration. For a more detailed definition of the lexicographic
filtration, we refer the reader to [37]. The persistence diagram
(defined in Sec. II-E) captures changes in the topology of the
sub-complexes as the filtration evolves, making the filtration
central to its construction.
<figure><img src="image_2.png" /><figcaption>Fig. 1. Filtration example for a PL scalar field representing the elevation on a terrain. As simplices are added by increasing scalar value, the topology of the domain changes: in (a), the terrain is composed of two connected components. In (b), only one connected component remains, but a handle has appeared. In (c), the smaller hill at the front has been filled and there is still a handle. In (d), there are two handles. In (e), only one handle is left. Finally, in (f), the whole domain has been swept (one connected component).</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. Example of discrete gradient field for a given scalar field (left). The larger spheres represent critical simplices. The light yellow arrows are edge-triangle pairs. The light blue arrows are vertex-edge pairs. Two examples of v-paths are shown in transparent blue (right), going from a critical edge and ending in critical vertices. These form the unstable set of the critical edge.</figcaption></figure>

## B. Distributed Model

The domain $\mathcal{M}$ is divided between $n_p$ processes in $n_p$ disjoint blocks of data. $f$ is defined on those blocks. Each process $p \in \{0, \dots, n_p - 1\}$ is associated with a local block $f_p : \mathcal{M}_p \to \mathbb{R}$, such that:

* $\mathcal{M}_p \subset \mathcal{M}$: each block $\mathcal{M}_p$ is a $d$-dimensional simplicial complex and a subset of $\mathcal{M}$.

* A simplex $\sigma$ is exclusively owned by the process with the lowest identifier $p$ whose block $\mathcal{M}_p$ contains $\sigma$.

* A simplex $\sigma \in \mathcal{M}_p$ which is not exclusively owned by the process $p$ is called a *ghost simplex* (see also below).

* $\cup_{\mathcal{M}_p} = \mathcal{M}$: the union of the blocks is equal to the input.

A layer of additional *ghost simplices* is added to each block $\mathcal{M}_p$. We note $\mathcal{M}_p'$ the $d$-dimensional simplicial complex obtained by considering a layer of ghost simplices, i.e., by adding to $\mathcal{M}_p$ the $d$-simplices of $\mathcal{M}$ which share a face with a $d$-simplex of $\mathcal{M}_p$, along with all their $d'$-dimensional faces (with $d' \in \{0, \dots, d-1\}$). Overall, all the simplices added in this way to the block $\mathcal{M}_p$ to form the *ghosted block* $\mathcal{M}_p'$ are referred to as *ghost simplices*. In practice, the partitioning of our domain and the generation of ghosts is performed using ParaView [44], a well-known visualization software. This preprocessing step is considered as already performed on our input data.

Our work is implemented in TTK and makes use of its features, in particular, the triangulation. The triangulation is a data structure that allows for efficient query execution and data traversal. It has been recently ported to a distributed-memory setting [40], which induces specific preconditioning steps to instantiate the triangulation, such as the ghost layer generation or the computation of the local adjacency graph of $\mathcal{M}_p$. This preconditioning overhead is negligible compared to already existing distributed algorithms of TTK.

## C. Discrete Morse Theory

In this section, we introduce notions of *Discrete Morse Theory* (DMT) [45] which is central to DMS. A discrete vector

is a pair formed by a simplex $\sigma_i \in \mathcal{M}$ (of dimension $i$) and one of its co-facets $\sigma_{i+1}$ (i.e., one of its co-faces of dimension $i+1$), noted $\{\sigma_i < \sigma_{i+1}\}$. $\sigma_{i+1}$ is usually referred to as the head of the vector, while $\sigma_i$ is its tail. Examples of discrete vectors include a pair between a vertex and one of its incident edges, or a pair between an edge and a triangle containing it. A discrete vector field on $\mathcal{M}$ is then defined as a collection $\mathcal{V}$ of pairs $\{\sigma_i < \sigma_{i+1}\}$, such that each simplex of $\mathcal{M}$ is involved in at most one pair. A simplex $\sigma_i$ which is involved in no discrete vector $\mathcal{V}$ is called a *critical simplex*.

A *v-path* is a sequence of discrete vectors $\{{\sigma_i^0 < \sigma_{i+1}^0}, \dots, \{\sigma_i^k < \sigma_{i+1}^k\}}$, such that the tails of two consecutive vectors are distinct and the tail of a vector in the sequence is a face of the head of the previous vector (see Fig. 2, right). The collection of all the v-paths terminating in a given critical simplex $\sigma_i$ is called the *discrete stable set* of $\sigma_i$. Symmetrically, the collection of all the discrete v-paths starting at a given critical simplex $\sigma_i$ is called the *discrete unstable set* of $\sigma_i$. Finally, a *discrete gradient field* is a discrete vector field such that all its possible v-paths are loop-free (see Fig. 2, left). Intuitively, this means that all critical simplices can be connected by following discrete vectors without any loops in the paths. The dimension of the critical simplex corresponds to its index in the smooth setting [46], [47]. A critical 0-simplex (or vertex) is called a minimum, a 1-simplex (or edge) a 1-saddle, a 2-simplex (or triangle) a 2-saddle and a 3-simplex (or tetrahedron) a maximum.

Several algorithms have been proposed to compute such a discrete gradient field from an input PL scalar field. We consider in this work the algorithm by Robins et al. [48], given its proximity to the PL setting: each critical cell identified by this algorithm is guaranteed to be located in the star of a PL critical vertex. In practice the computation is performed through inspection of the local neighborhood of each vertex, which makes this step embarrassingly parallel.
<figure><img src="image_4.png" /><figcaption>Fig. 3. Examples of 1-cycles on a sub-complex of a filtration. The PL scalar field represents the elevation on a terrain. Both the orange and the blue cycles are 1-cycles: the computation of their boundaries results in summing all the vertices of each cycle twice (as all vertices are faces of exactly two edges of the cycle), yielding a boundary equal to 0 for each cycle. However, the two cycles are not homologous: it is not possible to go from one to the other by adding a p-boundary. The blue cycle can be chosen as the representative of its class, representing the hole at the top of the hill it surrounds.</figcaption></figure>

### D. Homology groups

In this section, we introduce key concepts of topology relative to homology groups used in the original DMS algorithm. Intuitively, homology relies on the concept of cycles, used to detect "holes" in the data. In 3D, "holes" correspond to connected components, handles and voids. Cycles are a sum of simplices, or a collection of simplices, that can "go around" the domain and end where they started. Fig. 3 provides a visual example of these notions. Here are a series of definitions to formalize these concepts. A p-chain *c* is a formal sum of modulo 2 coefficients of p-simplices $\sigma_i$ of $\mathcal{M}: \sum \alpha_i \sigma_i$ with $\alpha_i \in \{0,1\}$. Two p-chains can be summed together component-wise to form a new p-chain. A p-chain can be modeled with a bit mask where a simplex is present if it has been added an odd-number of times. For example, adding the two 0-chains $a = v_0 + v_1$ and $b = v_1 + v_2$ will yield $a+b = v_0 + v_2$. The boundary of a p-simplex $\sigma_i$ is noted $\partial\sigma_i$ and is defined as the sum of its faces of dimension ($p-1$). The boundary of a p-chain *c* is the sum of the boundaries of its simplices: $\partial c = \sum \alpha_i \partial\sigma_i$. A p-cycle *c* is a p-chain such that $\partial c = 0$. The group of all p-cycles is noted $\mathcal{Z}_p(\mathcal{M})$.

Cycles are not sufficient to detect the "holes" of a domain. One must look into the relationship between cycles and p-boundaries. A p-boundary is a p-chain that is itself the boundary of a ($p+1$)-chain. The group of p-boundaries of a simplicial complex $\mathcal{M}$ is the group of all p-boundaries of $\mathcal{M}$, noted $\mathcal{B}_p(\mathcal{M})$. The $p^{th}$ homology group of a simplicial complex $\mathcal{M}$ is the quotient group of its p-cycles modulo its p-boundaries: $\mathcal{H}_p(\mathcal{M}) = \mathcal{Z}_p(\mathcal{M})/\mathcal{B}_p(\mathcal{M})$. Specifically, two p-cycles *c* and $c'$ are called *homologous* if one can be transformed into the other by adding the boundary of a ($p+1$)-chain $c''$: $c = c' + \partial c''$. The set of all cycles that are homologous defines a *homology class*. Any cycle of the class can be chosen as a representative. Intuitively, a homology class gathers all the cycles surrounding the same set of "holes". The rank of the $p^{th}$ homology group $\mathcal{H}_p(\mathcal{M})$, also known as the $p^{th}$ Betti number, is the number of linearly independent classes of $\mathcal{H}_p(\mathcal{M})$ (called *generators*). Detecting the p-dimensional

<figure><img src="image_5.png" /><figcaption>Fig. 4. Persistence diagram for a clean (left) and noisy (right) dataset. Critical simplices are denoted by spheres (light blue: minima, light yellow: maxima, others: saddles). The persistence of each pair is the height of the bar. $D_0$ corresponds to minimum-saddle pairs and $D_2$ corresponds to saddle-maximum pairs. Salient features (long pairs) can easily be distinguished from the noise (short pairs).</figcaption></figure>

independent "holes" in a domain therefore means extracting generators for the $p^{th}$ homology group.

### E. Persistence diagrams

As we browse the filtration, the topology of the domain evolves as shown in Fig. 1. The persistence diagram of dimension *p*, noted $D_p$, tracks the evolution of $p^{th}$ homology generators during a filtration. In particular, $D_0$ tracks the evolution of connected components, $D_1$ of handles and $D_2$ of voids. Specifically, at the step *i* of the filtration, a new generator $\gamma$ can appear upon the insertion of a p-simplex $\sigma_i$. For instance, a new 1-cycle $\gamma$ can appear by adding an edge which connects two vertices of $M_i$ belonging to the same connected component. In that case, $\gamma$ is said to be *born* at step *i*. Further in the filtration, at step *j*, $\gamma$ can disappear upon the insertion of a ($p+1$)-simplex $\sigma_j$ which completes the filling of the *hole* created at step *i* (i.e., in $M^j$, $\partial\sigma_j$ is *homologous* to the cycle created by $\sigma_i$ in $M^i$). For instance, a 1-cycle $\gamma$ made of three edges can be filled upon the introduction of the triangle containing these three edges. Then $\gamma$ is said to *die* at step *j*. The simplices $\sigma_i$ and $\sigma_j$ form a p-dimensional *persistence pair*. As shown in Fig. 4, each p-dimensional persistence pair is embedded in the persistence diagram $D_p$, as a point in the so-called *birth/death plane*, at coordinates $(f(\sigma_i), f(\sigma_j))$, where $f(\sigma)$ is the highest data value among the vertices of $\sigma$. The *persistence* of a pair can be read as the height of the point to the diagonal (i.e., $f(\sigma_j) - f(\sigma_i)$). Salient features are characterized by a high persistence, with the points located far from the diagonal. Conversely, noise in the data produces points in the plane close to the diagonal.

### F. The original Discrete Morse Sandwich algorithm

Here is an overview of the original Discrete Morse Sandwich algorithm. For a more detailed description of the algorithm, we refer the reader to its introduction paper [37]. We consider here 3D datasets, where $D_0$, $D_1$ and $D_2$ have to be computed. First, the discrete gradient is computed in parallel using multi-threading. The critical simplices are then deduced from the gradient. This step is called the *zero-persistence skip*: each remaining, non-critical simplex forms a zero-persistence pair with the other simplex involved in its discrete vector (Sec. II-C). The rest of the algorithm will focus on pairing the obtained critical simplices, following a stratification strategy where $D_0$ and $D_2$ (special cases) are computed before $D_1$.
<figure><img src="image_6.png" /><figcaption>Fig. 5. Overview of the DMS algorithm for the computation of $\mathcal{D}_0$. First, the unstable set is computed from the vertices of the critical 1-simplex $\sigma$ by tracing v-paths (white curves, left). The set is then collapsed into an extremum graph $\mathcal{G}_0$ (middle). Each critical simplex is represented by a node in the graph, with the arcs representing the v-paths. Finally, the graph is processed with a Union-Find structure to produce the persistence pair in $\mathcal{D}_0$ (right).</figcaption></figure>

The diagram $\mathcal{D}_0$ is computed first. An overview of this algorithm is shown in Fig. 5. We start by building an extremum graph noted $\mathcal{G}_0$ by following, for the two vertices $v_0$ and $v_1$ of a critical edge $\sigma$, the gradient until a critical 0-simplex (or extremum) is reached ($t_0$ and $t_1$). Each critical edge is processed in parallel (multi-threading). If $t_0$ and $t_1$ are distinct, the triplet ($\sigma, t_0, t_1$) represents new elements of the graph $\mathcal{G}_0$, adding the three nodes $\sigma, t_0$ and $t_1$ (one per element of the triplet) and two arcs representing the v-paths between $\sigma$ and the two extrema. When all critical 1-simplices are processed, $\mathcal{G}_0$ is complete. $\mathcal{D}_0$ is then computed by visiting the edges of $\mathcal{G}_0$ in increasing order following the PairExtremaSaddles algorithm presented in Alg. 1. This step is intrinsically sequential, and relies on a Union-Find data structure for each node of $\mathcal{G}_0$. A Union-Find efficiently models connectivity in data through two primitives: *find()*, that returns the representative of the connected component containing the node, and *union()*, that merges together two components by unifying their representatives. Initially, each 0-simplex is its own representative. For each triplet ($\sigma, t_0, t_1$) of $\mathcal{G}_0$, the following procedure is applied: the representatives $r_0$ and $r_1$ of $t_0$ and $t_1$ are retrieved (using *find()*). $r_0$ is ensured to be the highest representative (l. 5-6 in Alg. 1), by swapping if necessary its value with $r_1$. The highest representative, $r_0$, is then paired with $\sigma$ and is assigned $r_1$ as representative. In Fig. 5, this step creates the pairing ($\sigma, t_0$), with $t_1$ the new representative of $t_0$. To speed up the computation, the triplets of $\mathcal{G}_0$ are collapsed as they are visited. This means that the representative of $t_0$ is also set to $r_1$ (l. 12), which is equivalent to a *path compression* in a Union-Find data structure. This concludes the computation of $\mathcal{D}_0$.

The diagram $\mathcal{D}_2$ is computed similarly on critical 2- and 3-simplices using a dual discrete gradient field, obtained by reversing every discrete vector of the discrete gradient on the domain's dual complex (see [37] for more details).

For 3D datasets, $\mathcal{D}_1$ is computed last. The first step of the construction of $\mathcal{D}_1$ consists in restricting its input to the unpaired critical 1- and 2-simplices: the critical 1- and 2-simplices already paired in $\mathcal{D}_0$ or $\mathcal{D}_2$ are hence not considered here. This stratification strategy greatly reduces the number of input simplices for $\mathcal{D}_1$. We then apply Alg. 2 using the Pair-CriticalSimplex procedure defined in Alg. 3. For each unpaired critical 2-simplex $\sigma$, a homologous propagation (See Fig. 6) is performed. It expands a boundary, initially equal to $\partial\sigma$,

**Algorithm 1 PairExtremaSaddles**

**Input:** An ordered set $C_1$ of triplets ($\sigma_j, t_0, t_1$) of $\mathcal{G}_0$

**Output:** The persistence diagram $\mathcal{D}_0$, a collection of pairs ($\sigma, t$)

1: **for** $j \in C_1$ **do** // Process the 1-simplex $\sigma_j$

2: $r_0 \leftarrow findRepresentative(t_0)$

3: $r_1 \leftarrow findRepresentative(t_1)$

4: **if** $r_0 \neq r_1$ **then**

5: **if** $r_0 < r_1$ **then**

6: **swap**($r_0, r_1$)

7: **end if**

8: **addPair**($\sigma_j, r_0$)

9: **Representative**[$r_0$] $\leftarrow r_1$

10: **Representative**[$t_0$] $\leftarrow r_1$ // The arc is collapsed

11: **end if**

12: **end for**

<figure><img src="image_7.png" /><figcaption>Fig. 6. A homologous propagation given a simplex $\sigma$. The boundary of $\sigma$, $\partial\sigma$, is iteratively expanded by selecting its highest 1-simplex $\tau$ and adding to the boundary of $\sigma$ the boundary of the 2-chain associated to $\tau$. The boundary is expanded until reaching a critical 1-simplex (here in dark blue). The pair $(\sigma, \tau)$ is then added to $\mathcal{D}_1$.</figcaption></figure>

by selecting the highest 1-simplex $\tau$ of the current boundary and adding to it the boundary of the 2-chain associated to $\tau$. For simplicity, we will refer in the remainder of this paper to the boundary initiated this way in $\sigma$ as the *boundary of* $\sigma$. The propagation stops at the first unpaired 1-simplex. The pair $(\tau, \sigma)$ can then be created. $\sigma$ is the death and $\tau$ the birth of the pair. Note that the propagation is expended in reverse relative to the filtration order: given a triangle $\sigma$, this process identifies the edge $\tau$ which created the latest 1-cycle $\gamma$ in the filtration which is *homologous* to $\partial\sigma$ (Sec. II-E). All this assumes that the input set of simplices is ordered. However, as introduced by Nigmetov et al. in [49], it is possible to process the simplices in a random manner, hence in parallel (multi-threading). An extra case has then to be considered: when performing the homologous propagation for a simplex $\sigma$, it is possible to reach a 1-simplex $\tau$ that has already been paired to a 2-simplex $\sigma_{\tau}$ through homologous propagation. In that case, there are two possibilities: either (i) $\sigma_{\tau}$ is lower than $\sigma$ and the propagation carries on by merging $\partial\sigma_{\tau}$ with $\partial\sigma$, or (ii) $\sigma_{\tau}$ is higher than $\sigma$. In that case, the pair $(\tau, \sigma_{\tau})$ is removed, $(\tau, \sigma)$ is added and the homologous propagation of $\sigma_{\tau}$ is resumed. Compare-And-Swap operations are used for thread-safe memory accesses as described by Nigmetov et al. [49]. $\mathcal{D}_1$ is finally created from the temporary pairs once all
2-simplices have been completely processed.

**Algorithm 2 PairCriticalSimplices**

**Input:** Set $C_2$ of unpaired critical 2-simplices $\sigma$

**Output:** The persistence diagram $\mathcal{D}_1$, a collection of pairs
(Pair($\sigma$), $\sigma$)

1: **for** $j \in C_2$ **in parallel (multi-threading) do**

2: PairCriticalSimplex($\sigma_j$)

3: **end for**

4: **for** $j \in C_2$ **do**

5: $\mathcal{D}_1 \leftarrow \mathcal{D}_1 \cup (Pair(\sigma_j), \sigma_j)$

6: **end for**

**Algorithm 3 PairCriticalSimplex (homologous propagation)**

**Input:** An unpaired critical 2-simplex $\sigma$

**Output:** A temporary pair of $\mathcal{D}_1$

1: **if** Boundary($\sigma$) == 0 **then**

2: Boundary($\sigma$) $\leftarrow \partial\sigma$

3: **end if**

4: **while** Boundary($\sigma$) $\neq$ 0 **do**

5: $\tau \leftarrow \max(Boundary(\sigma))$

6: **if** $\tau$ is not a critical simplex **then** // Expand boundary

7: Boundary($\sigma$) $\leftarrow$

8: Boundary($\sigma$) + Boundary(Pair($\tau$))

9: **else** // $\tau$ is critical

10: **if** Pair($\tau$) == $\emptyset$ **then** // $\tau$ is unpaired

11: addPair($\sigma$, $\tau$)

12: **break**

13: **else** // $\tau$ has already been paired to $\sigma_{\tau}$

14: $\sigma_{\tau} \leftarrow Pair(\tau)$

15: **if** $\sigma_{\tau} < \sigma$ **then** // Merge the boundaries

16: Boundary($\sigma$) $\leftarrow$

17: Boundary($\sigma$) + Boundary($\sigma_{\tau}$)

18: **else** // $\sigma$ is older and the true death of $\tau$

19: addPair($\sigma$, $\tau$)

20: Pair($\sigma_{\tau}$) $\leftarrow \emptyset$

21: PairCriticalSimplex($\sigma_{\tau}$) // Resume for $\sigma_{\tau}$

22: **end if**

23: **end if**

24: **end if**

25: **end while**

Regarding 2D datasets, only the diagrams $\mathcal{D}_0$ and $\mathcal{D}_1$ have to be computed. $\mathcal{D}_1$ is then computed in 2D on critical 1- and 2-simplices like $\mathcal{D}_2$ in 3D.

G. Related work

**Persistent homology:** Multiple research groups independently introduced persistent homology [16], [50]–[52]. In numerous data analysis applications, topological persistence rapidly emerged as a compelling measure of importance, helping in the identification of salient topological features within the data. The most common method for persistence homology computation relies on the reduction of the boundary matrix (that describes the facet/co-facet relations between the simplices of the input domain). DMS relies on a different

strategy, based on discrete Morse Theory [45], [48], [53], but there are several conceptual similarities to existing documented accelerations. Bauer et al. [54] introduced the idea of apparent pairs, which is similar to the zero-persistence skip pre-computation step of DMS. Furthermore, the stratification strategy used in DMS can be linked with the stratification used by Bauer et al. [55] through "Clearing" and "Compression". These operations enable one to discard simplices already involved in persistent pairs. Edelsbrunner et al. [16] observe that the persistence diagram for dimension 0 can be computed through a Union-Find data structure. They also observe that the 2-dimensional persistence diagram can be obtained, by symmetry, with a Union-Find structure. This was aggressively exploited by Guillou et al. [37] as they restricted the Union-Find to stable and unstable sets of 1 and 2 saddles.

Other methods have investigated Morse Theory to accelerate the computation of persistent homology [46], [47], specifically in a discrete setting [45]. In particular, Robins et al [48] introduced the discrete gradient employed in DMS and used it to accelerate the computation of persistence. Other approaches improved on this idea [56]–[58] or extended it to a more general setting not limited to regular grids [53]. For instance, Wagner introduced an out-of-core approach [58] capable of processing massive scalar datasets on commodity hardware, but at the expense of significantly long computations (typically hours of computation for several billions of data points). In contrast, our approach targets high-performance hardware in a distributed setting, with much faster computations (Sec. VI-F). Also, note that this out-of-core approach [58] is not included in the performance benchmark by Guillou et al. [37], which was published before. However, our preliminary experiments report a typical ×2 speedup in favor of DMS [37] on our hardware. Overall, in contrast to previous approaches based on discrete Morse theory, DMS significantly accelerates the process thanks to an aggressive stratification strategy, taking advantage of the specificities of the diagrams $\mathcal{D}_0$ and $\mathcal{D}_2$ and avoiding the computation of the full Morse complex.

In practice, there are numerous software packages available to produce persistence diagrams, such as PHAT [59], DIPHA [41], Gudhi [60], Ripser [54] or Eirene [61]. Each implementation however focuses on particular data structures, such as generic filtrations of cell complexes (for PHAT, DIPHA and Gudhi) or Rips filtration of high-dimensional point clouds (for Eirene and Ripser). Some of the listed software for persistence diagram computation are purely sequential (Eirene, Ripser), while others implement shared-memory parallelism (PHAT, DMS). Another implementation, Ripser++, builds upon Ripser to focus on GPU acceleration, also using apparent pairs in a pre-computation step to speed up the overall algorithm. C. Tóth et al. [62] also introduced GPU-accelerated computations of signature kernels, a machine learning tool relying on a dynamic sequence of persistence diagrams [63]. Like Ripser++, their implementation focused on Vietoris-Rips filtrations. Otter et al. [64] provided in 2017 a benchmark focusing on execution speed and memory footprint on one and multiple CPU cores. For Vietoris-Rips filtrations, Ripser has been identified as the fastest implementation, whereas for sub-level set filtrations on cubical complexes, DIPHA was recognized then for its
superior execution speed. More recently, DMS has been re-
ported to be the fastest implementation using shared-memory
parallelism, according to the benchmark provided with its
introductory paper [37].

*Distributed-memory algorithms:* There are few existing ap-
proaches relative to distributed-memory parallelism and per-
sistence diagrams. The most well-known method is *DIPHA*,
introduced by Bauer et al. in [41]. It computes the boundary
matrix of the domain and partitions the matrix into blocks
of contiguous columns. The boundary matrix represents the
relations between the simplices and their faces. The matrix
is then reduced using a variant of Gaussian elimination. In
particular, it is very similar to the spectral sequence algorithm
for persistent homology [14], with several adaptations to make
it correct and efficient in a distributed-memory setting. Each
block is first reduced locally on a process. When the blocks
have been reduced to the best of a process’s capabilities, the
unreduced columns are sent to the next process to the left.
These communication and computation steps are performed
until all columns are completely reduced. The persistence
pairs of the diagram can then be extracted from the columns
and rows of the reduced matrix. *DIPHA* offers good parallel
speedups, using only multi-process (MPI) parallelism (no
multi-threading), and allows for the analysis of larger datasets
than anterior work.

Another distributed method was introduced by Ceccaroni et al. in [65]. However, their parallelism is limited to the concurrent processing of multiple, smaller datasets (i.e., at least one full dataset per process). On the contrary, we aim in this paper at processing one voluminous dataset on multiple compute nodes. An approach recently introduced by Nigmetov et al. [66] also uses distributed-memory parallelization to produce a persistence diagram. The algorithm combines spatial and range partitioning by computing first a local reduction of the data and then switching to a global reduction. However, this algorithm relies on persistent co-homology while our algorithm uses homology. As stated by Nigmetov et al., experiments suggest that persistent homology may be more efficient for computing the persistence diagram. Experiments that we performed with *DIPHA* using both homology and co-homology confirm this. At the time of writing this paper, no public implementation has been found for the approach of Nigmetov et al. To our knowledge, *DIPHA* is thus the only public implementation for distributed-memory computation of persistent homology.

In contrast to *DIPHA* and the approach of Nigmetov et al. (which both consider a distribution of the data by scalar value), our algorithm employs a spatial distribution of the input data and entirely obviates the need to compute the global boundary matrix. This feature constitutes a practical advantage, especially in contexts where data distribution has already been determined upstream in the analysis pipeline. Also, the reliance on discrete Morse Theory yields challenges that intrinsically differ from the solutions offered by the two other approaches.

III. OVERVIEW

This section provides an overview of our approach. First, the
global order of vertices is computed. This step is called *Array*

*Preconditioning.* In a distributed-memory setting, a global order is necessary for comparing vertices owned by different processes. The computation is done in parallel on multiple processes in three steps: we start by creating locally a vector for all vertices of the elements to sort (comprised of the scalar value of the vertex, its global identifier and the rank of the process that owns it). A distributed sort is then performed on the vector using `psort` [67]. This particular implementation was chosen because it is lightweight, modifiable and relatively efficient. *DIPHA* also uses this implementation to construct its boundary matrix. Each process can then compute the global order of the elements that are present locally following the distributed sort. Finally, the global orders are sent back to the owner of the corresponding vertices. For simplices of higher dimension, comparisons are performed using the lexicographic comparison on their global vertex orders (Sec. II-A).

Second, the discrete gradient of the input data is computed
by using the algorithm described by Robins et al. [48], which
is embarrassingly parallel for both shared- and distributed-
memory contexts [40].

Third, the critical simplices are extracted from the gradient
and sorted in the step *Extract & sort*.

Fourth, the diagrams $D_0$ and $D_2$ are computed by processing
the unstable and stable sets of the 1-saddles and 2-saddles of
$f$ and applying a self-correcting pairing algorithm to extract
the pairs of the diagram Sec. IV.

Next, the diagram $D_1$ is computed from the unpaired 1- and
2-saddles using our novel algorithm "DistributedPairCritical-
Simplices" Sec. V.

Finally, the classes of infinite persistence are extracted by
collecting the remaining, unpaired critical simplices.

As the input to the steps *Extract & sort*, *D*₀, *D*₂ and *D*₁
are the critical simplices, our overall algorithm is output-
sensitive (the number of critical simplices directly determines
the output size, the number of persistence pairs). For grids
of identical dimension, the number of simplices involved in
these steps depends on the input scalar fields, potentially
ranging from very few to a large number. The more pairs
are present in the output persistence diagram, the greater the
workload. Consequently, two datasets of the same dimension
may yield persistence diagrams of markedly different sizes,
with correspondingly varying computation times.

IV. EXTREMUM-SADDLE PERSISTENCE PAIRS

In this section, we will describe the different modifications
to PairExtremaSaddles (see Alg. 1) that we contributed to
obtain a distributed-memory version. This algorithm is applied
to compute both $D_0$ and $D_2$. An overview of the algorithm is
provided in Fig. 7.

A. *Stable and unstable sets computation*

The unstable sets of all critical 1-simplices are computed as
follows: a v-path is extracted from each vertex of each critical
1-simplex. Then, two possibilities arise: a critical vertex is
encountered (i.e. a local minimum), ending the computation
there for this unstable set, or a ghost vertex is encountered, in
which case a message is stored to be later sent to the process
<figure><img src="image_8.png" /><figcaption>Fig. 7. Overview of our algorithm for extremum-saddle pairing for $\mathcal{D}_0$. First, the unstable set is computed from the vertices of the critical 1-simplices $s_0$ and $s_1$ (sub-figure a, the two thick edges with their vertices) by following the gradient. When the computations of sets reach ghost simplices, a message is sent to the relevant process to notify it to resume computation on its domain (sub-figure b, the newly computed part of the set is thicker). Rounds of computations and communications are performed until all sets are computed (sub-figure c). The set is then collapsed into the distributed graph $\mathcal{G}_0$ following the rules stated in Sec. IV-B (sub-figure d). The ownership of extremum $t_0$ is given to process 0 as it is the process with the lowest rank identifier that has a set that ends in $t_0$. The ownership of extremum $t_1$ is given to process 2 as it has an unstable set started in $M_2$, at $s_1$, that ends in $t_1$. The graph is then processed with a distributed Union-Find structure to produce the persistence pair in $\mathcal{D}_0$. The thicker arc ($s_0$, $t_1$) corresponds to the computed persistence pair.</figcaption></figure>

owning the ghost vertex so that the computation can resume there later. Once all computations either are completed or have generated a message, all processes exchange their stored messages and resume their computations on their block. These successive computation and communication steps run until no messages are sent on any process during a communication round. For stable sets, the computation is similar but is applied on 2-simplices as the start of the set and 3-simplices as the final simplex of the set. The gradient is followed in reverse to emulate the dual gradient without explicitly computing it.

## B. Distributed extremum graph construction

The previous step computes the stable and unstable sets. Now we collapse those sets into the distributed extremum graphs $\mathcal{G}_i$ with $i \in \{0, 2\}$. In the remainder, we focus on the case $i=0$, the case $i=2$ being symmetric. The nodes of the graph are the saddles and extrema of the previously computed sets (Sec. IV-A). The arcs of the graph represent the v-paths connecting these critical simplices. A triplet of the graph refers to the three nodes $(\sigma, t_0, t_1)$, where $\sigma$ is a critical saddle and $t_0$ and $t_1$ extrema linked to $\sigma$ by v-paths.

Nodes of this graph may be located on different processes. We therefore establish a few additional rules to fit the definitions of $\mathcal{G}_0$ to a distributed-memory setting. These rules were established to enable the computation of the next step of the algorithm (i.e. the self-correcting pairing).

* Saddles are present on only one process. This arises from the manner in which representatives are computed in distributed memory: the process is initiated at a saddle and subsequently proceeds from one extremum to another, thereby eliminating the need for ghost saddles. It is necessary to determine which process holds the next extremum of the representative computation, whereas this requirement does not apply to saddles.

* A saddle node in $\mathcal{G}_0$ is owned by a process $p$ if its associated critical simplex is exclusively owned by $p$ (Sec. II-B). This ownership limits computation and communication at this step, as graph saddles are created

where their originating saddles in the triangulation already exist.

* Extrema, however, can be present on multiple processes but they are owned by only one and are ghost on other processes. This design choice was driven by the need to determine to which process to resume the representative computation in the self-correcting algorithm. With the presence of ghosts, a necessary link is made between processes.

* The ownership of an extremum node is determined as follows. If an extremum simplex $e$ is exclusively owned by the process $p$ (Sec. II-B) and if there exists a saddle simplex $s$ also exclusively owned by $p$ such that one of its unstable sets terminates in $e$, then the extremum node of $\mathcal{G}_0$ associated to $e$ is owned by $p$. Otherwise, the node associated to $e$ is owned by the process with the lowest rank identifier that owns a saddle node whose unstable set ends in $e$. The ownership of the extrema was designed to limit the number of communications. However, as the construction of the graph and the self-correcting algorithm are saddle-centric, it is more convenient to ensure that an extremum is on the same process of at least one of the saddles it is linked to in the graph, hence the ownership design.

Fig. 7 shows an example of a distributed extremum graph and applied ownership rules. The local graph on a process $p$ is noted $\mathcal{G}_{0,p}$. The ghosted local graph of a process $p$ is noted $\mathcal{G}'_{0,p}$. An extremum node is called at the interface of two processes $p$ and $q$ if it is owned by either $p$ or $q$ and is a ghost on the other process.

The computation of the collapse of the sets to build a graph is fairly straightforward: once all the sets are computed, the processes possess lists of all sets ending on a local extremum they own with regard to their domain. Each process will then determine which process is the owner of the extremum in the graph and send back to the owner of the originating saddle the extremum and its new ownership. It will also send to the new owner of the node a list of all the processes on which the
<figure><img src="image_9.png" /><figcaption>Fig. 8. Illustration of our self-correcting pairing algorithm on an example of $G_0$. Sub-figure a represents a distributed extremum graph on 3 processes, 0 (in blue), 1 (in orange) and 2 (in yellow). Extrema at the interface of two processes intersect with the dashed grey lines separating the domain of the processes. Sub-figure b represents the first computation iteration. Each process computes on its local graph without taking other processes into consideration. On process 2, $s_2$ is paired with $t_4$ (shown here by a thicker arc) and $t_0$ is made the new representative of $t_4$ (shown here by a dashed arrow). Then the pairing of $s_3$ is performed: $s_3$ can be paired with either $t_3$ or the representative of $t_4$, $t_0$. Since $t_3$ is higher than $t_0$, $s_3$ is paired with $t_3$. On process 1, $s_1$ is also paired with $t_3$. After the computation, processes exchange data regarding relevant, shared pairings. Process 2 will tell process 1 it created the pairing $(s_3, t_3)$ with $t_0$ as the new representative, because $t_3$ is owned by process 1. Process 1 will tell process 2 about its owned pairing $((s_1, t_3),$ with $t_2$ as the new representative) because it knows that $t_3$ is present as a ghost on the graph of process 2. Process 2, upon receiving this message, will assess that the message is correct and its pairing $(s_3, t_3)$ is wrong. It will compute a new pairing: $(s_3, t_2)$ (as $t_2$ is the new representative of $t_3$) and tell process 1 about the pairing (sub-figure c). Process 1 will update $t_2$ in the pairing to its representative, $t_1$. As $t_1$ is owned by process 0, process 1 will tell it about the pairing. Process 0, knowing that $t_1$ is unpaired, will accept the pairing $(s_3, t_1)$. Then, it will update the representative of $t_1$ to $t_0$. Finally, it will send back the information to process 2 that the actual correct pairing is $(s_3, t_1)$ as shown in sub-figure d.</figcaption></figure>

extremum is a ghost with regard to the graph. Each process will then receive and build the parts of its local graph.

### C. Self-correcting distributed pairing

We now have to build the $D_0$ pairs from the local graph $G'_{0,p}$ (and similarly for $D_2$). This step was originally performed sequentially in DMS since its execution time was negligible compared to others. However, in a distributed-memory setting, a sequential execution would most likely greatly hurt the overall performance as it would prevent any speedup on multiple nodes for this particular step. A detailed analysis of the execution time associated with DMS reveals that, for both $D_0$ and $D_2$, this sequential stage can account for up to 8.6% of the total runtime for Magnetic Reconnection, one of our test datasets (see Sec. VI-A). According to Amdahl's law, the theoretical upper bound on the multi-node speedup under these conditions is $\frac{1}{0.086} \approx 11.6$. This limitation would therefore impose a significant constraint on the attainable efficiency of the algorithm as a whole. Also, and most importantly, sequential execution on one node would require to gather the intermediate data (i.e., the graph computed in Sec. IV-B) on the node. This would introduce additional communications and may create an additional memory bottleneck in an algorithm that has already high memory requirements.

We thus have to design a distributed pairing algorithm. The new algorithm heavily relies on algorithm Alg. 1, but introduces several components to face the challenges brought by distributed-memory parallelism. We incorporated a self-correction mechanism, enriched the representation of extrema with additional data, removed the path compression mechanism, employed maps and local identifiers in place of global identifiers, and introduced inter-process communications where necessary. The following paragraphs provide a detailed account of these modifications.

We use the idea of comparing saddles to see which is the oldest (or the youngest), but restrict the algorithm to the distributed extremum graph created at the previous step. Our algorithm design is also inspired by the self-correcting

mechanism of the parallel (multi-threaded) procedure Pair-CriticalSimplices, computing $D_1$ in the original DMS algorithm. This introduces incorrect pairings resulting in extra computations. Such overhead is however handled in parallel, ultimately resulting in overall performance gains (see Sec. VI). Comparisons to the original DMS algorithm in Fig. 14 and Tab. I also show that our method incurs limited overhead. Moreover, we expected other strategies relying on a more direct approach (i.e. without incorrect pairings) to induce too much synchronizations and idle time to perform well.

A key difference from the original DMS algorithm is that the representative of an extremum will now store two pieces of information: the representative and the identifier of the saddle that assigned the representative to the extremum. This enables, when computing the pairing of a saddle $\sigma$, to stop the computation of the representatives of its extrema nodes $t_0$ and $t_1$ when the loop reaches representatives assigned by saddles older than $\sigma$ (as this would not occur in sequential). The path compression mechanism of DMS mentioned in Sec. II-F (see Alg. 1) is no longer applied as the resulting representatives may be false, leading to potentially incorrect computations for other pairings during the computation of other representatives. When a wrong pairing is detected by saddle comparison, the computation of the representatives of its original triplet is restarted from the beginning using Alg. 4.

Here is a description of the overall self-correcting distributed algorithm. A practical example is shown in Fig. 8. First, each process executes *DistributedProcessTriplet* (See Alg. 4) for all its triplets (in a sorted manner, as it proved to be more efficient even though it is not required). The messages to be sent are stored until all computations are performed, then all processes exchange their messages. For each message $(\sigma, m_0, m_1)$ received, the process will first detect if it is a recomputation (encoded by $(\sigma, -1, -1)$) and will trigger it. If it is not a recomputation, the process will update the representatives of the message to $(\sigma, Representative(m_0, \sigma), Representative(m_1, \sigma))$. If the newly computed representatives belong to another process,
<caption>**Algorithm 4** DistributedProcessTriplet</caption>

**Input:** A triplet $(\sigma, t_0, t_1)$ of $\mathcal{G}_{0,p}$ on process $p$
**Output:** A temporary pair of $\mathcal{D}_0$
1: **for** $i$ in $\{0, 1\}$ **do** // Compute the representatives
2:     $r_i \leftarrow \text{findRepresentative}(t_i, \sigma)$
3:     $r_i \text{Paired} \leftarrow \text{Pair}(r_i) \neq \emptyset$
4:     $r_i \text{Invalid} \leftarrow r_i \text{Paired}$ and $\sigma < \text{Pair}(r_i)$
5:     **if** $r_i \text{Invalid}$ **then**
6:         $r_i \text{Paired} \leftarrow \text{false}$
7:     **end if**
8: **end for**
9: **if** $r_0 \notin \mathcal{G}'_{0,p}$ or $r_1 \notin \mathcal{G}'_{0,p}$ **then**
10:     Send $(\sigma, r_0, r_1)$ to owner
11: **end if**
12: **if** $r_0 == r_1$ **then**
13:     **if** $\text{Pair}(\sigma) \neq \emptyset$ **then** // $\sigma$ should not be paired
14:         $r \leftarrow \text{Pair}(\sigma)$
15:         $\text{Representative}[r] \leftarrow r$ // re-initialize the representative
16:         $\text{Pair}(\sigma), \text{Pair}(r) \leftarrow \emptyset, \emptyset$ // remove the invalid pair
17:     **end if**
18: **else**
19:     **if** $(r_0 < r_1$ or $r_0 \text{Paired}$) and $!(r_1 \text{Paired})$ **then**
20:         swap $r_0$ and $r_1$ data
21:     **end if**
22: **if** $!(r_0 \text{Paired})$ **then**
23:     **if** $r_0 \text{Invalid}$ **then** // remove the invalid pair
24:         $\sigma_k \leftarrow \text{Pair}(r_0)$
25:         $\text{Pair}(\sigma_k) \leftarrow \emptyset$
26:     **end if**
27:     $\text{Pair}(\sigma), \text{Pair}(r_0) \leftarrow r_0, \sigma$ // Add the new pair
28:     $\text{Representative}[r_0] \leftarrow (r_1, \sigma)$
29:     **if** $r_0$ is at the interface of $p$ and process $q$ **then**
30:         Send $(\sigma, r_0, r_1)$ to process $q$
31:     **end if**
32: **if** $r_0 \text{Invalid}$ **then** // Recompute the invalid $\sigma_k$
33:     **if** $\sigma_k \in \mathcal{G}_{0,p}$ **then**
34:         DistributedProcessTriplet($\sigma_k$)
35:     **else**
36:         Send recomputation signal to the owner of $\sigma_k$
37:     **end if**
38: **end if**
39: **end if**
40: **end if**

then the message is passed along to that other process. Otherwise, the process detects if the message should trigger a correction through saddle comparison. If so, the pairs and the representatives are updated and a recomputation is triggered for the invalid saddle. This cycle of communications and computations is repeated until no messages are sent on any process during a communication round.

The original DMS algorithm used vectors to store several variables relative to critical simplices in its *PairCriticalSimplices* algorithm. These vectors were defined for all simplices of the triangulation, even though they were used only for critical simplices. This allowed for fast memory accesses, as the index of the simplex in the triangulation was equal its index in the vectors. In a distributed-memory setting, this is no longer possible: the global number of simplices will most likely prevent the memory allocation of such large vectors. We therefore reduced the size of the vectors to the number of local critical simplices and used maps to convert a global simplex index to its index in these vectors.

## D. Shared-memory parallelism

The computations of stable and unstable sets are straightforwardly processed with multiple threads in each MPI process as all set computations are independent. The communications are however performed only by the OpenMP primary thread. The construction of the distributed graph is similarly parallelized.

The self-correcting pairing algorithm is not multi-threaded within each MPI process. Although feasible in principle, we made the choice not to explore shared-memory parallelism for this particular step. Our expectations were that multi-thread parallelization would yield only limited performance gains, an assessment informed by the performance results reported for a comparable algorithm by Smirnov et al. [68]. Furthermore, as shown by Guillou et al. [37] (Table 3, Appendix C), the computation of $\mathcal{D}_0$ and $\mathcal{D}_2$ with DMS is often limited in terms of computation time with regard to the other procedures. Since this step already benefits from distributed-memory parallelism, we considered that the additional complexity introduced by thread-level parallelization would most likely outweigh the potential performance benefits.

$\mathcal{D}_0$ and $\mathcal{D}_2$ being completely independent, we assigned each diagram computation to an OpenMP task that can itself generate threads in a nested manner.

## V. SADDLE-SADDLE PERSISTENCE PAIRS

In this section, we will describe the different modifications to the original *PairCriticalSimplex* (Alg. 3) and *PairCriticalSimplex* (Alg. 2) algorithms necessary for the efficient distributed computation of $\mathcal{D}_1$.

### A. Distributed-memory parallel algorithm

Our algorithm relies on a new data structure: the global-local boundary. For a 2-simplex $\sigma$, this structure is composed of two elements: the set of edges of the boundary initiated in $\sigma$ (called in the remainder, for simplicity, the *boundary of* $\sigma$), which are owned by the current process $p$ (called *local boundary*), and the highest boundary edges of all processes containing a part of the boundary of $\sigma$ (called the *global boundary*). The local boundary is identical to the boundary of the original DMS algorithm. The global boundary is updated by other processes all through the computation. At any given time, on a given process $p$, the highest edge $\sigma$ of its local boundary is always lower or equal to the highest edges reported by the global boundaries of the other processes.

A central idea in our distributed-memory algorithm is the notion of *computation token*. Each 2-simplex $\sigma$ for which a homologous propagation needs to be computed is associated with a token. At any time, the token of each propagation is present in only one process. Only the process owning the token is allowed to propagate the boundary of $\sigma$. This means that, taken individually, each propagation is carried out sequentially. However, all the propagations are computed in parallel similarly to the shared-memory context.

Here is a description of our algorithm *DistributedPairCriticalSimplex* as defined in Alg. 5 and illustrated in Fig. 9. This algorithm revisits *PairCriticalSimplex* Alg. 3 in order
<figure><img src="image_10.png" /><figcaption>Fig. 9. Distributed homologous propagation. Sub-figure a shows two critical simplices: $\sigma$, the death saddle, and $\tau$, the birth saddle. $\sigma$ is located on process 1 (light yellow) and $\tau$ is located on process 0 (light blue). Initially, the boundary of $\sigma$ is only located on process 1 (Sub-figure b). New edges are added to the boundary that are owned by process 1 (Sub-figure c). Then new edges are added to the boundary that are ghost for process 1 (Sub-figure d) and a message is sent to process 0 so that it adds the edges to its local boundary of $\sigma$. On process 1, a new piece of information is stored regarding the boundary: the highest edge for process 0 (equal to the highest edge of the two ghost edges). As the global highest edge is still located on process 1 (gray arrow), the computation continues on process 1 (Sub-figure d) and additional ghost edges are added (Sub-figure e). After adding two new ghost edges (Sub-figure f), the highest edge is now located on process 0. Process 1 will send the computation token to process 0. Upon receiving it, process 0 will resume the computation and propagate the boundary (Sub-figure f). Finally, the propagation ends in $\tau$ and the pair $(\sigma, \tau)$ is created on process 0.</figcaption></figure>

to compute a distributed homologous propagation for an unpaired critical 2-simplex $\sigma$. For an unpaired critical 2-simplex $\sigma \in M_p$, with $p$ a process, its local boundary is propagated by following the same rules as *PairCriticalSimplex*. However, when a ghost edge, owned by process $q$, needs to be added to the global-local boundary of $\sigma$, a message will be sent to $q$ so that $q$ adds this particular edge to its local boundary of $\sigma$. The propagation can trigger a merge between two boundaries. Local boundaries are merged similarly to *PairCriticalSimplex* (Alg. 3, 1.15-17). Global boundaries are merged by keeping the highest of the two edges for each process. A message will then be sent to notify all the relevant processes that a merge has occurred and should be performed by them as well. As soon as the highest edge of the local boundary is no longer the highest edge in the global boundary, the propagation on $p$ is stopped. Then, the computation token will be sent to the process owning the highest edge in the global boundary, to resume the propagation on its block.

At the end of the computation, the pair $(\tau, \sigma)$ is stored on the process that owns $\tau$, as the boundary of another propagation may reach $\tau$ and a comparison between the two originating 2-simplices may be required. The process that owns $\sigma$ does not have to be aware of which simplex completed the pair.

*DistributedPairCriticalSimplex* hence generates two types of messages: computation tokens and boundary updates. Boundary updates correspond to either a merge order between two global-local boundaries, an addition of an edge to a local boundary or an update of the highest edge in a global boundary. The received boundary updates have to be performed in a particular order, to ensure two properties: (i) updates from potentially multiple processes with regard to one particular propagation need to be received and processed in the same order they were created in, (ii) updates from one process with regard to potentially multiple propagations need also to be received and processed in the same order they were created in. Other orders may result in an incorrect outcome. For (i), the property can be ensured by following a round-by-round design with alternating communication and computation steps. Processing the boundary updates sequentially for each process will ensure property (ii). However, messages sent by different processes can be processed in any order as long as they involve different propagations.

**Algorithm 5 DistributedPairCriticalSimplex**

**Input:** An unpaired critical 2-simplex $\sigma$
**Output:** A temporary pair of $D_1$
1: **if** GlobalLocalBoundary($\sigma$) == 0 **then**
2:     addEdge(GlobalLocalBoundary($\sigma$), $\partial\sigma$)
3:     **end if**
4:     **while** GlobalLocalBoundary($\sigma$) $\neq$ 0 **do**
5:         $\tau$ $\leftarrow$ max(LocalBoundary($\sigma$))
6:         $\tau_p$ $\leftarrow$ max(GlobalBoundary($\sigma$))
7:         **if** max($\tau$, $\tau_p$) == $\tau_p$ **then** // $\tau_p$ is the highest edge
8:             UpdateMaxGlobal($\sigma$, $\tau$, GlobalBoundary($\sigma$))
9:             Mark computation token of $\sigma$ for sending to $p$
10:        **return**
11:    **end if**
12:    **if** $\tau$ is not a critical simplex **then**
13:        addEdge(GlobalLocalBoundary($\sigma$), $\partial(Pair(\tau))$)
14:    **else**
15:        **if** Pair($\tau$) == $\emptyset$ **then**
16:            addPair($\sigma$, $\tau$) // $\tau$ is unpaired
17:            UpdateMaxGlobal($\sigma$, $\tau$, GlobalBoundary($\sigma$))
18:            **break**
19:        **else**
20:            $\sigma_\tau$ $\leftarrow$ Pair($\tau$) // $\tau$ has already been paired to $\sigma_\tau$
21:            **if** $\sigma_\tau < \sigma$ **then**
22:                MergeGlobalLocalBoundaries($\sigma$, $\sigma_\tau$)
23:            **else**
24:                addPair($\sigma$, $\tau$) // $\sigma$ is older and the true death of $\tau$
25:                UpdateMaxGlobal($\sigma$, $\tau$, GlobalBoundary($\sigma$))
26:                Pair($\sigma_\tau$) $\leftarrow$ $\emptyset$
27:                PairCriticalSimplex($\sigma_\tau$) // Resume for $\sigma_\tau$
28:            **end if**
29:        **end if**
30:    **end if**
31: **end while**

The overall algorithm *DistributedPairCriticalSimplex* (Alg. 6) revisits *PairCriticalSimplex* (Alg. 2) with distributed processes, and following this round-by-round design. All the propagations are first computed locally using multi-thread parallelism. Once all propagations are either completed or their computation token needs to be sent to another process, the communications start. First, the boundary updates are exchanged and processed. Then the computation tokens are exchanged and the new propagations are computed in parallel using threads. These rounds of communications and computations are performed until all critical 2-simplices are paired.
<figure><img src="image_11.png" /><figcaption>Fig. 10. Example of anticipation of propagation computation for 4 processes (process 0 in light blue, process 1 in dark blue, process 2 in orange and process 3 in light yellow). A critical triangle (in yellow) is the starting point of the distributed homologous propagation. The boundary is propagated as described in Alg. 5, however the computation token is not sent when the global highest edge of the boundary becomes located on process 1 (Sub-figure b). Instead, the propagation is continued on process 3 until the number of propagation iterations (i.e. while loop iterations in Alg. 5) reaches a predefined counter. Then, the computation token is sent to the process owning the global highest edge (here, process 1). The boundary is propagated on process 1 (Sub-figure c) until an unpaired critical edge is reached (thick, dark, blue edge). At this point, the global highest edge is located on process 2 (orange). The blue critical edge is therefore not paired with the critical triangle and the computation token is sent to process 2 that resumes the propagation on its domain. When reaching its own critical edge, as the global highest edge is still located on the domain of process 2, the propagation ends here and the pair is created on process 2 (Sub-figure d). Instead of having to send the computation token back and forth between process 1, 2 and 3 to produce the final pairing, the anticipation of propagation enables to exchange the token only twice.</figcaption></figure>

<caption>**Algorithm 6 DistributedPairCriticalSimplices**</caption>

Input: Set $C_2$ of unpaired critical 2-simplices $\sigma$
Input: Set $C_1$ of unpaired critical 1-simplices $\tau$
Output: The persistence diagram $\mathcal{D}_1$, a collection of pairs $(\sigma, \tau)$
1: **for** $j \in C_2$ **in parallel (multi-threading) do**
2: &nbsp;&nbsp;&nbsp;&nbsp;DistributedPairCriticalSimplex($\sigma_j$)
3: &nbsp;&nbsp;&nbsp;&nbsp;**end for**
4: &nbsp;&nbsp;&nbsp;&nbsp;**while** Global number of terminated propagations < $|C_2|$ **do**
5: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; // Perform global boundary updates
6: &nbsp;&nbsp;&nbsp;&nbsp; Send boundary updates to other processes
7: &nbsp;&nbsp;&nbsp;&nbsp; Receive boundary updates from other processes
8: &nbsp;&nbsp;&nbsp;&nbsp; Update received boundaries
9: &nbsp;&nbsp;&nbsp;&nbsp; // Resume computations with tokens
10: Send computation tokens to other processes
11: Receive computation tokens from other processes
12: **for** all received tokens $\sigma$ **in parallel (multi-threading) do**
13: &nbsp;&nbsp;&nbsp;&nbsp; DistributedPairCriticalSimplex($\sigma$)
14: &nbsp;&nbsp;&nbsp;&nbsp; **end for**
15: **end while**
16: **for** $j \in C_1$ **do** // Extract pairs from boundary computation
17: &nbsp;&nbsp;&nbsp;&nbsp; $\mathcal{D}_1 \leftarrow \mathcal{D}_1 \cup (\sigma_j, \text{Pair}(\sigma_j))$
18: **end for**

Similarly to Sec. IV-C, our implementation reduced the size of data structures and vectors to the local number of critical simplices and used maps to convert a global simplex index to its index in vectors.

## B. Anticipation of propagation computation

The algorithm presented in the previous section presents a critical flaw: in the worst case scenario, for boundaries stretched out on multiple processes, the maximum may change process every time a simplex is added to the propagation, resulting in repeated exchanges of the computation token. This may generate an extremely high number of communications. Anticipating this back and forth is possible by changing slightly the original algorithm as shown in Fig. 10. Instead of stopping the computation and sending the computation token

to another process as soon as the highest edge of the boundary is located on another process, we further the computation regardless on the current process, until either the number of propagation iterations (i.e., while loop iterations in Alg. 5) reaches a predefined counter (arbitrarily equal to 0.01% of the number of triangles of $M_p'$) or until an unpaired critical simplex $c$ is reached. Only then is the computation token sent to the process that owns the highest edge. Not pairing the potential simplex $c$ ensures that the propagation never expands too far, leading to potentially incorrect pairs, which would have been difficult to detect and correct afterwards.

## C. Overlap of communication and computation

There are two limitations to the previous algorithm that we want to address in this section: thread idle time and cost of communication. Indeed, at the end of each computation round, when waiting for all work to complete, there is often just a few propagations being computed, resulting in significant idle time when using many threads. We aim at reducing this idle time by triggering a communication round before all computations are finished. On the other hand, the cost of communication can be reduced by effectively overlapping communications with computations at the MPI level.

A dedicated communication thread can solve both these problems. We preserve the round-by-round structure for the communications to ensure that the update of global-local boundaries is processed in the right order (see Sec. V-A), but the communication rounds are now triggered by the communication threads. In each MPI process, the communication thread sends and receives messages, updates boundary data and creates one OpenMP task for each propagation, while the other compute threads process the propagations.

Even if we aim at starting earlier each communication round, these are not triggered as soon as one message can be sent. Making the communication thread wait until reaching a
threshold and sending multiple messages in one MPI communication at once limits indeed the number of communications. It ensures that the MPI layer is not overloaded with numerous messages, and limits the number of OpenMP atomic operations performed by the communication thread: these atomic operations are required for a correct synchronization with the compute threads. Messages will only be sent if the number of messages waiting to be sent by the current process is above the following threshold. It is set dynamically to increase reactivity as the computation progresses. At first, it is equal to 0.01% of the local number of unpaired 2-simplices. Then, at every round of communications, it is updated using the remaining global number of unpaired critical 2-simplices to add reactivity to the communications. This value of 0.01% was determined empirically. In the case where there are no tasks left to be computed within the current process, messages are sent regardless of whether the threshold was reached. The communication thread also performs the update of global-local boundaries. This can be done in parallel of propagation computation as updating global-local boundaries from the current round will not interfere with the computation of propagations from previous rounds. It is because the updates are not directly related to the current computation tokens.

The compute threads can now continuously process the local and incoming propagations through a task pool filled by the communication thread, harnessing more efficiently the intra-node multi-core parallelism. The idle time of the threads is therefore significantly reduced as compute threads no longer have to wait at the end of each computation round and the cost of communication is effectively hidden with the overlap of communication with computation as usual with a communication thread [69], [70].

## VI. RESULTS

For the following results, we rely on Sorbonne Université's supercomputer, MCMeSU, which contains 48 nodes of 32 cores each. Each node is composed of 2 AMD EPYC 7313 Milan CPUs with 256GB of RAM. The nodes are interconnected with Mellanox Infiniband. In our tests, we use up to 16 nodes (512 cores total), with one MPI process and 32 threads per node to minimize MPI communications and synchronizations as well as the memory footprint. When using a communication thread (see Sec. V-C), we rely on 31 compute threads only. Our algorithm is implemented in C++ with MPI+OpenMP within TTK [38], [39]. The correctness of our implementation was checked for all test datasets, by comparing our outputs against those generated by DMS (which were already compared to DIPHA's for triangulated voxel data, see [37] for more details). Tests of strong and weak scaling are conducted to study the performance of our algorithm. The preconditioning time for TTK's distributed triangulation (Sec. II-B) is not accounted for in this paper (it is negligible for regular grids [40]). Specifically, when reporting execution times, we consider that the data is already distributed among the nodes, in the form of ghosted blocks (Sec. II-B), which is a standard input for analysis pipelines in distributed environments. Only the execution time of our algorithm and its direct preconditioning

are measured. Our algorithm is also compared to the original DMS algorithm as well as to DIPHA.

### A. Datasets

The performance of our software has been evaluated using multiple datasets, selected to demonstrate a broad spectrum of cases. These datasets are sourced from publicly available repositories [71], [72]. Datasets from the OpenSciVis collection can be prepared using the script available in the reproducible example (See Sec. I).

* Backpack: a CT scan of a backpack containing various items. This acquired dataset is spatially imbalanced with regard to its number and location of topological features, making it a good test case for workload balancing.

* Isabel: a simulation of hurricane Isabel that hit the east coast of the USA in 2003. This dataset is very smooth and will result in a persistence diagram with few but significant pairs.

* Elevation: a synthetic dataset of the altitude within a cube, with a unique maximum at one corner of the cube and a unique minimum at the opposite corner. This is a pathological case with only one class of infinite persistence in $D_0$.

* Wavelet: a synthetic dataset of 3D wavelets on a cube. This dataset is quite smooth and symmetric, which should result in a very good workload balance.

* Isotropic pressure: a direct numerical simulation of forced isotropic turbulence. This dataset is originally very big (4096³) and was down-sampled.

* Magnetic reconnection: A single time step from a computational simulation of magnetic reconnection, showing interactions between magnetic fields. This dataset is extremely noisy.

* Synthetic truss: A simulated CT scan of a truss with defects. This dataset possesses very rich and symmetric topological features.

* Random: a synthetic dataset of a random field on a cube. This dataset showcases a pathological case with a high number of spatially evenly distributed persistence pairs.

For the strong scaling benchmarks, all datasets were resampled to 512³ via trilinear interpolation, except for Random, that was resampled to 512² × 256 as the execution time for this dataset is particularly long for 512³ with all tested softwares, making it unpractical to manage. This smaller size still makes Random the dataset with the longest execution time, as it is our worst case scenario.

For the weak scaling benchmarks, the size of the input (number of vertices) doubles each times the number of nodes doubles (by doubling the number of vertices along one, alternating, dimension). The initial size on one node is the same as the strong scaling one (512² × 256 for Random, 512³ for the others). The datasets were re-sampled in different ways depending on the size of the original data: Isabel, Backpack, Magnetic Reconnection have been up-sampled, whereas Synthetic truss and Isotropic Pressure have been down-sampled. Random was generated for its biggest weak scaling case and then down-sampled to smaller datasets. Elevation was
<figure><img src="image_12.png" /><figcaption>Fig. 11. Performance impact of the different $D_1$ versions on the overall DDMS execution time for Backpack. Basic corresponds to the first version for $D_1$ described in Sec. V-A, Anticipation to the second one described in Sec. V-B, and Overlap to the final one described in Sec. V-C.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 12. Execution time of each step of DDMS, for strong (left) and weak (right) scalings for Backpack.</figcaption></figure>

generated for each size, so that it always has only one pair in its diagram. Due to its symmetry, Wavelet was generated for the largest weak scaling case and then resized by being cut in two along each dimension alternatively.

## B. Performance improvements

We start by assessing the performance improvements of our different versions for computing $D_1$, namely: Basic, the initial version (see Sec. V-A); Anticipation, that implements the anticipation of computation for $D_1$ (see Sec. V-B) and Overlap that iterates on Anticipation and adds the overlap of communication and computations thanks to the communication thread (see Sec. V-C). As shown in Fig. 11, the anticipation of computation dramatically improves the overall DDMS performance, making Anticipation 6 times faster than Basic on 16 nodes in strong scaling and over 12 times in weak scaling. Overlap also improves the performance, by adding reactivity to the execution. On 16 nodes, it reduces the overall execution time by 20% in weak scaling and 28% in strong scaling. These results validate and justify our modifications, which are hence necessary to efficiently deploy such TDA algorithms on multiple nodes.

From now on, we will only consider the Overlap version, whose detailed execution profile is presented in Fig. 12. For both strong and weak scalings, the Array Preconditioning step, which corresponds to the computation of global order of vertices (see Sec. III), is quite short and minority. In strong scaling, Array Preconditioning is marginally more expensive on two nodes than on one, after which execution time decreases and ultimately outperforms single-node performance.

<figure><img src="image_14.png" /><figcaption>Fig. 13. Parallel efficiency of DDMS for strong (left) and weak (right) scaling.</figcaption></figure>

It is due to the cost of the extra computation steps described in Sec. III, mainly the preparation of the vector and the computation of the global order from the sorted vector. Despite this increase, the computation of this step remains a minority.

The discrete gradient step scales very well as expected, which is an important source of overall performance gains. The computations of $D_0$ and $D_2$ scale also well, in both weak and strong scaling. The following step, Extract & Sort, corresponds to the extraction and local sort of critical simplices for all dimensions (see Sec. III). It is performed independently on each process and also scales well. Only the computation time of $D_1$ (which is the most intensive in terms of time complexity) remains constant in strong scaling and scales unfavorably in weak scaling. But, thanks to our successive improvements presented in Sec. V, there is no strong performance loss. Taking advantage of the embarrassingly parallel steps of our algorithm, we still manage to obtain overall significant performance gains when increasing the number of nodes.

## C. Strong scaling

The results for all datasets in strong scaling are shown in Fig. 13 (left) in terms of parallel efficiency with respect to the execution on one core (i.e., speedup divided by the number of cores, in percentage). The execution times are also available in Fig. 14. With the exception of Random, the efficiencies of all datasets fit within the range of 55% to 80% on one node and of 20% to 50% on 512 cores. This shows the scalability of our approach. Though the efficiencies decrease as the number of cores increases, Fig. 14 shows that the execution times continue to decrease even on 512 cores, with most datasets eventually requiring less than 20 seconds on 512 cores. These performance results demonstrate that the benefits of multi-node execution extend beyond increased memory availability: they also encompass the acceleration of computations that could already be performed on a single node.

Random behaves a bit worse than the other datasets, presenting the biggest drop in efficiency: from the best efficiency on one node (close to 100%) to one of the worst one on 512 cores (close to 26%). This is explained by the output-sensitivity of our algorithm (see Sec. III). Another factor is the spatial placement of the birth and death of a pair within the dataset. The further apart they are, the longer the computation will be. Random is one of the noisiest of our datasets (with Magnetic Reconnection and Synthetic Truss), however, unlike those two datasets its pairs are evenly distributed. Consequently, the
<figure><img src="image_15.png" /><figcaption>Fig. 14. Comparison between DMS (dotted lines) and DDMS (full lines) in terms of execution time (left) and per-node peak memory footprint (right).</figcaption></figure>

birth and death tend to be further apart spatially, requiring more work and more communications. This leads to very good efficiency on one node, but this quickly becomes a performance issue as the number of nodes increases, leading to such an efficiency drop.

## D. Weak scaling

We define the weak scaling efficiency for *p* processes as $e_p = \frac{t_1}{t_p} \times 100$, with $t_1$ and $t_p$ being the execution times on 1 and *p* nodes. The weak scaling results for all datasets are shown in Fig. 13 (right). The weak scaling efficiency is better than the strong scaling one for most datasets. This is partly due to the fact that doubling the dataset size through re-sampling often results in less than a twofold increase in the number of critical simplices, and hence in computational workload (given the output sensitivity of our algorithm). For most datasets, the efficiency is in the range of 35% to 80% on 16 nodes, which again shows the scalability of our approach. There are however two exceptions: Random and Magnetic Reconnection. For Random, the efficiency eventually drops lower to 17% for the same reasons as in strong scaling: its pairs are numerous and spatially stretched out. For Magnetic Reconnection, the efficiency largely exceeds 100%. This is due to the up-sampling of the original dataset that barely multiplies the number of pairs by a factor of 1.6 between 1 and 16 nodes. This is most likely because the topological features are already numerous and unevenly distributed across the dataset. While this also applies to other datasets, such as Isabel, the other specificity of Magnetic Reconnection is that it produces the most pairs out of all the datasets on 512³ (35 millions). The execution time for $D_0$, $D_1$ and $D_2$ is therefore substantial compared to the computation of other steps, such as the gradient. As the number of nodes increases, even though the size of the dataset increases, each process actually computes less and less pairs leading to lower execution times. The same applies for Isabel, but has no impact on the overall execution time, as the computation of $D_0$, $D_1$ and $D_2$ is originally negligible compared to the gradient computation.

## E. Performance comparison

**Comparison with DMS:** In Fig. 14 are shown the execution time and per-node peak memory footprint of our DDMS algorithm compared to the original DMS algorithm on one core and then one node. In Tab. I are shown the gains

<table><caption>TABLE I
Gain of DDMS over DMS (in percentage). A positive gain (bold) indicates that DDMS outperforms DMS in terms of execution speed.</caption><thead><tr><th>#cores</th><th>Synthetic</th><th>Wavelet</th><th>Elevation</th><th>Magnetic</th><th>Isabel</th><th>Random</th><th>Backpack</th><th>Isotropic</th><th>Average</th></tr></thead><tbody><tr><td>1</td><td>-20.38</td><td><b>2.14</b></td><td><b>2.95</b></td><td>-26.11</td><td><b>2.15</b></td><td>-16.54</td><td>-3.98</td><td>-3.95</td><td>-8.00</td></tr><tr><td>32</td><td>-0.83</td><td><b>23.39</b></td><td><b>24.87</b></td><td>-10.86</td><td><b>14.68</b></td><td><b>26.87</b></td><td><b>7.55</b></td><td><b>14.71</b></td><td><b>12.5</b></td></tr></tbody></table>

<figure><img src="image_16.png" /><figcaption>Fig. 15. Comparison based on execution time (left) between *DIPHA* (dotted lines) and DDMS (full lines) and performance gain (right) for a strong scaling setting. The performance gain on a given number of cores is defined as $t_{DIPHA}/t_{DDMS}$, with $t_{DIPHA}$ and $t_{DDMS}$ the execution times of *DIPHA* and DDMS respectively. A performance gain higher than 1 means that DDMS is faster than *DIPHA*.</figcaption></figure>

of DDMS over DMS for executions on 1 and 32 cores. A positive gain indicates that DDMS outperforms DMS in terms of execution speed. On 1 core, the gain is on average -8%. The datasets characterized by larger workloads (in terms of number of critical simplices), namely Random, Synthetic Truss, and Magnetic Reconnection, exhibit the lowest gains. This suggests that gain decreases as workload size increases. Then, on 32 cores, the gain increases for all datasets, yielding comparable execution times of both algorithms, the distributed DDMS algorithm executing slightly faster than its shared-memory counterpart for all datasets, except Synthetic Truss and Magnetic Reconnection. The DDMS extra cost for these two datasets is due to changes in the algorithm (for example, the removal of arc collapse in the $D_0$ and $D_2$ computations in Sec. IV-C) or to more costly data structures (Sec. IV-B), enabling MPI execution. The overhead is however very limited on average and executing DDMS on two nodes already outperforms DMS for all datasets.

In terms of peak memory footprint, DDMS uses significantly less memory for all datasets. This is due to our reduction of vector size to the number of critical simplices as mentioned in Sec. IV-C and Sec. V-A. This allowed DDMS to produce an overall smaller footprint on one node.

**Comparison with DIPHA:** We now compare our algorithm to *DIPHA*, to our knowledge the only publicly available MPI implementation (without multithreading) for persistence diagram computation. The *DIPHA* execution time is measured using the built-in *benchmark* mode and corresponds to its total execution without the I/O time. We start by comparing the execution times in strong scaling on Fig. 15. On one core, DDMS outperforms *DIPHA* only on the smoother datasets (Elevation, Wavelet, Isabel and Backpack), as these can really harness the preconditioning of the discrete gradient to speed up the rest of the computation. On multiple nodes, DDMS scales much better and hence outperforms *DIPHA* for all datasets starting from 4 nodes. Notice that on one and two nodes,
<figure><img src="image_17.png" /><figcaption>Fig. 16. Maximum over all nodes of the per-node peak memory footprint (left) for *DIPHA* (dotted lines) and DDMS (full lines) and memory gain (right) for a strong scaling setting. The memory gain on a given number of nodes is defined as $m_{DIPHA}/m_{DDMS}$, with $m_{DIPHA}$ and $m_{DDMS}$ the maximum over all node of the per-node peak memory footprint for *DIPHA* and DDMS respectively. A memory gain higher than 1 means that *DIPHA* uses more memory than DDMS.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 17. Persistence diagram on a subset of Turbulent Channel Flow (2048x1920x1536, 6 billion vertices). The dataset (left) is the pressure field of a direct simulation of a turbulent flow in a plane channel. The execution was performed on 16 nodes of 32 cores and 256GB of RAM each (512 cores and 4096GB of RAM total). The persistence diagram (right: $D_0$ in blue, $D_1$ in orange, $D_2$ in yellow) was computed in 174 seconds (19 million pairs).</figcaption></figure>

only one dataset out of eight (Magnetic Reconnection) is more efficiently processed by *DIPHA*. Moreover, considering both execution times and scaling, the worst case dataset for *DIPHA* (Synthetic Truss) scales relatively well for DDMS, whereas the worst case for DDMS (Random) is always processed faster (up to ×3) by DDMS on more than one core. Finally, the average speedup for all datasets is around ×8 on 512 cores, showing a substantial performance gain of DDMS over *DIPHA*.

Fig. 16 compares the memory consumption of both approaches. On the fewest number of nodes, *DIPHA* requires a lower footprint than DDMS for all datasets but the memory scalability is much better for DDMS. Hence, on the largest number of nodes the peak memory footprint of DDMS ends up being smaller than that of *DIPHA*. DDMS divides indeed by almost two its memory usage (with a small overhead due to the ghost cells) every time the number of cores is multiplied by two (Sec. II-B).

For *DIPHA*, its data distribution is more conducive to memory imbalance as the final reduced columns are stored on the process that completed their reduction, meaning that some process may store a significantly higher number of reduced columns than other processes. This imbalance is more likely as the number of processes increases, leading to *DIPHA* having a larger per-node peak memory footprint than DDMS.

### F. Example

We ran DDMS on a larger dataset (Turbulent Channel Flow [71]) to show our algorithm's capability to handle massive datasets. This dataset is a direct simulation of a fully developed flow at different Reynolds numbers in a plane channel. The scalar field is the three dimensional pressure field and has been converted to single-precision floating-point numbers for a lower memory consumption. The computation was run on a subset of the original dataset of size 2048x1920x1536 which is approximately 6 billion vertices. The memory bottleneck of our implementation mainly lies in the computation of the discrete gradient, which is the most memory-consuming step. For this example, we used TTK's compile option `TTK_ENABLE_DCG_OPTIMIZE_MEMORY` which optimizes the gradient memory footprint (by trading on its computation time). The execution was performed on 16 nodes of 32 cores

and 256GB of RAM each (512 cores and 4096GB of RAM total). The persistence diagram shown in Fig. 17 was computed in 174 seconds and contained a little under 19 million pairs. For comparison, out-of-core techniques [58] require several hours of computation for datasets of similar size, which further illustrates the practical interest of our work for high-performance computing contexts.

### G. Limitations

For completeness, we discuss here some limitations of our work. Our implementation of the DDMS algorithm only supports structured grids at the time of writing this paper, even though TTK support both structured and unstructured grids. Although the original DMS implementation has been validated on both, additional tests would be required to validate our DDMS implementation on unstructured grids. Certain aspects of practical use may be different when using unstructured meshes. For example, as unstructured grids require more memory, our tests and benchmarks may need some modifications. The performance may also be more memory-bound. Substantial additional tests are thus required to evaluate the performance of our algorithm on such grids. Additionally, several specialized domain representations popular in scientific computing such as adaptive mesh refinement (AMR) are not supported by TTK and therefore by our work.

The primary memory bottleneck of our implementation lies in the computation of the gradient. Indeed, this computation is not only costly in terms of execution time but also in terms of memory footprint, as discrete vectors are computed for all simplices of all dimensions. This generates a significant memory footprint, which may be addressed in the future by an improved compact encoding of the discrete vectors, exploiting the regular structures of structured grids.

Finally, similarly to DMS, since our work also exploits Robins' discrete gradient [48], it suffers from the same limitation regarding an extension to higher dimensions. In particular, strong guarantees on the restriction of critical cells upon homology changes of the sub-complexes are provided in up to three dimensions. Beyond, additional, spurious critical cells may appear, limiting the effectiveness of our approach.
## VII. CONCLUSION

This paper introduced a new algorithm for the efficient computation of the persistence diagram on scalar data in a distributed-memory setting using a hybrid MPI+thread parallelization. The performance of our algorithm was tested on a set of 8 datasets representing various use cases. Thanks to our algorithmic improvements for $D_0, D_1$ and $D_2$, and to the scalability of the discrete gradient computation, we have shown that our approach delivers significant speedups on up to 512 cores, with parallel efficiencies up to 50% in strong scaling and up to 80% in weak scaling. DDMS also outperforms *DIPHA*, the only publicly available implementation for computing the persistence diagram in distributed memory, by a factor of $\times 8$ on average. It produces a slightly larger memory footprint than *DIPHA* on few nodes, but a smaller one on 16 nodes.

Additionally, we showed that our algorithm is capable of processing massive datasets by running DDMS on a larger dataset (Turbulent Channel Flow [71]) of 6 billion vertices in under 3 minutes on 512 cores.

In the future, we will investigate strategies to reduce the memory footprint of the discrete gradient. Indeed, as previously stated, the primary memory bottleneck of our implementation lies in the computation of the gradient. Reducing its memory footprint while maintaining its execution speed would be beneficial for DDMS. An avenue for further research is deploying DDMS on heterogenous architectures. Given the gradient algorithm's embarrassingly parallel nature, CPU architectures with integrated GPUs, such as the AMD Instinct MI300A APU, may further accelerate DDMS by leveraging GPU cores for the gradient computation and CPU cores for the computation of $D_0, D_1$ and $D_2$.

Another research direction is out-of-core computation, which enables processing large datasets on a single compute server by loading data in and out of memory. Specifically, extensions of our self-correcting pairing (for $D_0$ and $D_2$, Sec. IV-C) as well as our anticipation strategy for homologous propagation (for $D_1$, Sec. V-B) could be considered for out-of-core contexts, but significant adaptations would be required to maintain the time performance of our approach.

## ACKNOWLEDGMENTS

This work is partially supported by the European Commission grant ERC-2019-COG "TORI" (ref. 863464, https://erc-tori.github.io/).

## REFERENCES

[1] D. E. Laney, P. Bremer, A. Mascarenhas, P. Miller, and V. Pascucci, "Understanding the Structure of the Turbulent Mixing Layer in Hydrodynamic Instabilities," *IEEE Trans. Vis. and Comp. Graph.* (*Proc. of IEEE VIS*), 2006.

[2] P. Bremer, G. Weber, J. Tierny, V. Pascucci, M. Day, and J. Bell, "Interactive Exploration and Analysis of Large Scale Simulations Using Topology-based Data Segmentation," *IEEE Trans. Vis. and Comp. Graph.*, 2011.

[3] A. Gyulassy, P. Bremer, R. Grout, H. Kolla, J. Chen, and V. Pascucci, "Stability of Dissipation Elements: A case study in combustion," *Euro-Vis: Proc. of Eurographics Conference on Visualization*, 2014.

[4] J. Kasten, J. Reininghaus, I. Hotz, and H. Hege, "Two-Dimensional Time-Dependent Vortex Regions Based on the Acceleration Magnitude," *IEEE Trans. Vis. and Comp. Graph.*, 2011.

[5] F. Nauleau, F. Vivodtzev, T. Bridel-Bertomeu, H. Beaugendre, and J. Tierny, "Topological Analysis of Ensembles of Hydrodynamic Turbulent Flows - An Experimental Study," in *IEEE Symposium on Large Data Analysis and Visualization*, 2022.

[6] D. Maljovec, B. Wang, P. Rosen, A. Alfonsi, G. Pastore, C. Rabiti, and V. Pascucci, "Topology-Inspired Partition-Based Sensitivity Analysis and Visualization of Nuclear Simulations," in *IEEE Pacific Vis*, 2016.

[7] F. Chazal, L. J. Guibas, S. Y. Oudot, and P. Skraba, "Persistence-Based Clustering in Riemannian Manifolds," *Journal of the ACM*, 2013.

[8] H. Doraiswamy, J. Tierny, P. J. S. Silva, L. G. Nonato, and C. T. Silva, "TopoMap: A 0-dimensional Homology Preserving Projection of High-Dimensional Data," *IEEE Trans. Vis. and Comp. Graph.* (*Proc. of IEEE VIS*), 2020.

[9] H. Bhatia, A. G. Gyulassy, V. Lordi, J. E. Pask, V. Pascucci, and P.-T. Bremer, "TopoMS: Comprehensive Topological Exploration for Molecular and Condensed-Matter Systems," *Journal of Computational Chemistry*, 2018.

[10] D. Guenther, R. Alvarez-Boto, J. Contreras-Garcia, J.-P. Piquemal, and J. Tierny, "Characterizing Molecular Interactions in Chemical Systems," *IEEE Trans. Vis. and Comp. Graph.* (*Proc. of IEEE VIS*), 2014.

[11] M. Olejniczak, A. S. P. Gomes, and J. Tierny, "A Topological Data Analysis Perspective on Non-Covalent Interactions in Relativistic Calculations," *International Journal of Quantum Chemistry*,

[12] T. Sousbie, "The Persistent Cosmic Web and its Filamentary Structure: Theory and Implementations," *Royal Astronomical Society*, 2011.

[13] N. Shivashankar, P. Pranav, V. Natarajan, R. van de Weygaert, E. P. Bos, and S. Rieder, "Felix: A topology based framework for visual exploration of cosmic filaments," *IEEE Trans. Vis. and Comp. Graph.*, 2016.

[14] H. Edelsbunner and J. Harer, *Computational Topology: An Introduction*. American Mathematical Society, 2009.

[15] C. Heine, H. Leitte, M. Hlawitschka, F. Iuricich, L. De Floriani, G. Scheuermann, H. Hagen, and C. Garth, "A Survey of Topology-based Methods in Visualization," *Computer Graphics Journal*, 2016.

[16] H. Edelsbunner, D. Letscher, and A. Zomorodian, "Topological persistence and simplification," *Discrete and Computational Geometry*, 2002.

[17] S. Parsa, "A Deterministic $o(m \log m)$ Time Algorithm for the Reeb Graph," in *Symposium on Computational Geometry*, 2012.

[18] C. Gueunet, P. Fortin, J. Jomier, and J. Tierny, "Task-based Augmented Reeb Graphs with Dynamic ST-Trees," in *Eurographics Symposium on Parallel Graphics and Visualization*, 2019.

[19] S. Tarasov and M. Vyali, "Construction of contour trees in 3D in $O(n \log n)$ steps," in *Symposium on Computational Geometry*, 1998.

[20] H. Carr, J. Snoeyink, and U. Axen, "Computing contour trees in all dimensions," in *Symp. on Dis. Alg.*, 2000.

[21] C. Gueunet, P. Fortin, J. Jomier, and J. Tierny, "Task-Based Augmented Contour Trees with Fibonacci Heaps," *IEEE Trans. Parallel Distrib. Syst.*, 2019.

[22] P. Bremer, H. Edelsbunner, B. Hamann, and V. Pascucci, "A Topological Hierarchy for Functions on Triangulated Surfaces," *IEEE Trans. Vis. and Comp. Graph.*, 2004.

[23] A. Gyulassy, V. Natarajan, V. Pascucci, P. Bremer, and B. Hamann, "A Topological Approach to Simplification of Three-Dimensional Scalar Functions," *IEEE Trans. Vis. and Comp. Graph.*, 2006.

[24] L. De Floriani, U. Fugacci, F. Iuricich, and P. Magillo, "Morse complexes for shape segmentation and homological analysis: discrete models and algorithms," *Computer Graphics Journal*, 2015.

[25] A. Gyulassy, P. Bremer, and V. Pascucci, "Shared-Memory Parallel Computation of Morse-Smale Complexes with Improved Accuracy," *IEEE Trans. Vis. and Comp. Graph.* (*Proc. of IEEE VIS*), 2018.

[26] J. Tierny and V. Pascucci, "Generalized topological simplification of scalar fields on surfaces," *IEEE Trans. Vis. and Comp. Graph.* (*Proc. of IEEE VIS*), 2012.

[27] J. Lukasczyk, C. Garth, R. Maciejewski, and J. Tierny, "Localized Topological Simplification of Scalar Data," *IEEE Trans. Vis. and Comp. Graph.* (*Proc. of IEEE VIS*), 2020.

[28] H. A. Carr, J. Snoeyink, and M. van de Panne, "Simplifying Flexible Isosurfaces Using Local Geometric Measures," in *IEEE VIS*, 2004.

[29] A. Bock, H. Doraiswamy, A. Summers, and C. T. Silva, "TopoAngler: Interactive Topology-Based Extraction of Fishes," *IEEE Trans. Vis. and Comp. Graph.* (*Proc. of IEEE VIS*), 2018.

[30] M. Soler, M. Plainchault, B. Conche, and J. Tierny, "Topologically Controlled Lossy Compression," in *IEEE Pacific Vis*, 2018.

[31] —, "Lifted Wasserstein Matcher for Fast and Robust Topology Tracking," in *IEEE Symposium on Large Data Analysis and Visualization*, 2018.
[32] M. Soler, M. Petitfrere, G. Darche, M. Plainchault, B. Conche, and J. Tierny, "Ranking Viscous Finger Simulations to an Acquired Ground Truth with Topology-Aware Matchings," in *IEEE Symposium on Large Data Analysis and Visualization*, 2019.

[33] J. Lukasczyk, C. Garth, G. H. Weber, T. Biedert, R. Maciejewski, and H. Leitte, "Dynamic Nested Tracking Graphs," *IEEE Trans. Vis. and Comp. Graph.* (*Proc. of IEEE VIS*), 2019.

[34] G. Favelier, N. Faraj, B. Summa, and J. Tierny, "Persistence Atlas for Critical Point Variability in Ensembles," *IEEE Trans. Vis. and Comp. Graph.* (*Proc. of IEEE VIS*), 2018.

[35] J. Vidal, J. Budin, and J. Tierny, "Progressive Wasserstein Barycenters of Persistence Diagrams," *IEEE Trans. Vis. and Comp. Graph.* (*Proc. of IEEE VIS*), 2019.

[36] M. Kontak, J. Vidal, and J. Tierny, "Statistical Parameter Selection for Clustering Persistence Diagrams," in *2019 IEEE/ACM UrgentHPC@SC*, 2019.

[37] P. Guillou, J. Vidal, and J. Tierny, "Discrete Morse Sandwich: Fast Computation of Persistence Diagrams for Scalar Data – An Algorithm and A Benchmark," *IEEE Trans. Vis. and Comp. Graph.*, 2023.

[38] T. Bin Masood, J. Budin, M. Falk, G. Favelier, C. Garth, C. Gueunet, P. Guillou, L. Hofmann, P. Hristov, A. Kamakshidasan, C. Kappe, P. Klacansky, P. Laurin, J. Levine, J. Lukasczyk, D. Sakurai, M. Soler, P. Steneteg, J. Tierny, W. Usher, J. Vidal, and M. Wozniak, "An Overview of the Topology ToolKit," in *TopolnVis*, 2019.

[39] J. Tierny, G. Favelier, J. A. Levine, C. Gueunet, and M. Michaux, "The Topology ToolKit," *IEEE Trans. Vis. and Comp. Graph.* (*Proc. of IEEE VIS*), 2017, https://topology-tool-kit.github.io/.

[40] E. Le Guillou, M. Will, P. Guillou, J. Lukasczyk, P. Fortin, C. Garth, and J. Tierny, "TTK is Getting MPI-Ready," *IEEE Trans. Vis. and Comp. Graph.*, 2024.

[41] U. Bauer, M. Kerber, and J. Reininghaus, "Distributed Computation of Persistent Homology," in *Proceedings of the Workshop on Algorithm Engineering and Experiments*, 2014.

[42] H. Edelsbrunner and E. P. Mucke, "Simulation of Simplicity: A Technique to Cope with Degenerate Cases in Geometric Algorithms," *ACM Trans. on Graphics*, 1990.

[43] H. Freudenthal, "Simplizialzerlegungen von beschränkter Flachheit," *Annals of Mathematics*, vol. 43, pp. 580–582, 1942.

[44] J. Ahrens, B. Geveci, and C. Law, "ParaView: An End-User Tool for Large-Data Visualization," *The Visualization Handbook*, 2005.

[45] R. Forman, "A User's Guide to Discrete Morse Theory," *AM*, 1998.

[46] J. Milnor, *Morse Theory*. Princeton University Press, 1963.

[47] M. Morse, "The Calculus of Variations in the Large," in *American Mathematical Society*, 1934.

[48] V. Robins, P. J. Wood, and A. P. Sheppard, "Theory and Algorithms for Constructing Discrete Morse Complexes from Grayscale Digital Images," *IEEE Trans. Pattern Anal. Mach. Intell.*, 2011.

[49] D. Morozov and A. Nigmetov, "Brief Announcement: Towards Lockfree Persistent Homology," in *Symposium on Parallelism in Algorithms and Architectures*, 2020.

[50] S. A. Barannikov, "The Framed Morse Complex and Its Invariants," in *ADVSOV*. American Mathematical Society, 1994.

[51] P. Frosini and C. Landi, "Size Theory as a Topological Tool for Computer Vision," *Pattern Recognition and Image Analysis*, 1999.

[52] V. Robins, "Toward Computing Homology from Finite Approximations," *Topology Proceedings*, 1999.

[53] K. Mischaikow and V. Nanda, "Morse Theory for Filtrations and Efficient Computation of Persistent Homology," *Discrete and Computational Geometry*, 2012.

[54] U. Bauer, "Rips: efficient computation of Vietoris-Rips persistence barcodes," https://github.com/Ripser/ripser/, 2019.

[55] U. Bauer, M. Kerber, and J. Reininghaus, "Clear and Compress: Computing Persistent Homology in Chunks," in *Topological Methods in Data Analysis and Visualization III*. Springer International Publishing, 2014, pp. 103–117.

[56] D. Günther, J. Reininghaus, H. Wagner, and I. Hotz, "Efficient computation of 3D Morse-Smale complexes and persistent homology using discrete Morse theory," *The Visual Computer*, 2012.

[57] F. Iuricich, "Persistence cycles for visual exploration of persistent homology," *IEEE Trans. Vis. and Comp. Graph.*, 2021.

[58] H. Wagner, "Slice, Simplify and Stitch: Topology-Preserving Simplification Scheme for Massive Voxel Data," in *Symposium on Computational Geometry*, 2023.

[59] U. Bauer, M. Kerber, J. Reininghaus, and H. Wagner, "Phat - Persistent Homology Algorithms Toolbox," *Journal of Symbolic Computation*, 2017.

[60] C. Maria, J. Boissonnat, M. Glisse, and M. Yvinec, "The Gudhi Library: Simplicial Complexes and Persistent Homology," in *Mathematical Software*, 2014, https://github.com/GUDHI/.

[61] G. Henselman-Petrasek, "Eirene.jl package for Homological Algebra," https://github.com/Eetion/Eirene.jl, 2018.

[62] C. Tóth, D. J. D. Cruz, and H. Oberhauser, "A user's guide to KSig: Gpu-accelerated computation of the signature kernel," *arXiv preprint arXiv:2501.07145*, 2025.

[63] D. Lee and H. Oberhauser, "The Signature Kernel," *arXiv preprint arXiv:2305.04625*, 2023.

[64] N. Otter, M. Porter, U. Tillmann, P. Grindrod, and H. Harrington, "A roadmap for the computation of persistent homology," *EPJ Data Science*, vol. 6, p. 17, 2017.

[65] R. Ceccaroni, L. Di Rocco, U. Ferraro Petrillo, and P. Brutti, "A Distributed Approach for Persistent Homology Computation on a Large Scale," *Journal of Symbolic Computation*, 2024.

[66] A. Nigmetov and D. Morozov, "Distributed Computation of Persistent Cohomology," 2024. [Online]. Available: https://arxiv.org/abs/2410.16553.

[67] D. Cheng, V. Shah, J. Gilbert, and A. Edelman, "A novel parallel sorting algorithm for contemporary architectures," 01 2007. [Online]. Available: https://github.com/DIPHA/dipha/tree/master/externals/psort-1.0

[68] D. Smirnov and D. Morozov, "Triplet Merge Trees," in *TopolnVis*, 2017.

[69] A. Denis and F. Trahay, "MPI Overlap: Benchmark and Analysis," in *International Conference on Parallel Processing*, 2016.

[70] G. Hager, G. Schubert, T. Schoenemeyer, and G. Wellein, "Prospects for Truly Asynchronous Communication with Pure MPI and Hybrid MPI/OpenMP on Current Supercomputing Platforms," in *Proceedings of The Cray User Group*, 2011.

[71] P. Klacansky, "Open Scientific Visualization Data Sets," https://klacansky.com/open-scvis-datasets/, 2020.

[72] TTK Contributors, "TTK Data," https://github.com/topology-tool-kit/ttk-data/tree/dev, 2020.

VIII. BIOGRAPHY SECTION

<figure><img src="image_19.png" /></figure>

Eve Le Guillou is a Ph.D. student at Sorbonne Université. She received a M.S. degree in 2020 in Computer Science from Cranfield University as well as an engineering degree in 2021 from École Centrale de Lille. She is an active contributor to the Topology ToolKit (TTK), an open source library for topological data analysis. Her notable contributions to TTK include the port to MPI of TTK's data structure and of several algorithms.

<figure><img src="image_20.png" /></figure>

Pierre Fortin received his Ph.D. degree in Computer Science from the University of Bordeaux in 2006. He was an Assistant Professor at Sorbonne University from 2007, and he is Professor at University of Lille since 2023. His research interests include parallel algorithmic and high performance scientific computing.

<figure><img src="image_21.png" /></figure>

Julien Tierny received the Ph.D. degree in Computer Science from the University of Lille in 2008. He is a CNRS research director at Sorbonne University. Prior to his CNRS tenure, he held a Fulbright fellowship and was a post-doctoral researcher at the University of Utah. His research expertise lies in topological methods for data analysis and visualization. He is the founder and lead developer of the Topology ToolKit (TTK), an open source library for topological data analysis.