Scalable High-Quality Hypergraph Partitioning

Lars Gottesbüren

Karlsruhe Institute of Technology
Am Fasanengarten 5

Tobias Heuer

Karlsruhe Institute of Technology

Nikolai Maas

Karlsruhe Institute of Technology

Peter Sanders

Karlsruhe Institute of Technology

Sebastian Schlag

Independent Researcher

Balanced hypergraph partitioning is an NP-hard problem with many applications, e.g., optimizing communication in distributed data placement problems.
The goal is to place all nodes across k different blocks of bounded size, such that hyperedges span as few parts as possible. This problem is well-studied
in sequential and distributed settings, but not in shared-memory. We close this gap by devising efficient and scalable shared-memory algorithms for all
components employed in the best sequential solvers without compromises with regards to solution quality.

This work presents the scalable and high-quality hypergraph partitioning framework Mt-KaHyPar. Its most important components are parallel
improvement algorithms based on the FM algorithm and maximum flows, as well as a parallel clustering algorithm for coarsening – which are used in a
multilevel scheme with $\log(n)$ levels. As additional components, we parallelize the $n$-level partitioning scheme, devise a deterministic version of
our algorithm, and present optimizations for plain graphs.

We evaluate our solver on more than 800 graphs and hypergraphs, and compare it with 25 different algorithms from the literature. Our fastest
configuration outperforms almost all existing hypergraph partitioners with regards to both solution quality and running time. Our highest-quality
configuration achieves the same solution quality as the best sequential partitioner KaHyPar, while being an order of magnitude faster
with ten threads. Thus, two of our configurations occupy all fronts of the Pareto curve for hypergraph partitioning.
Furthermore, our solvers exhibit good speedups, e.g., 29.6x in the geometric mean on 64 cores (deterministic), 22.3x ($\log(n)$-level), and 25.9x

<concept_desc>Theory of computation Shared memory algorithms</concept_desc>
<concept_desc>Mathematics of computing Hypergraphs</concept_desc>
<concept_desc>Theory of computation Graph algorithms analysis</concept_desc>

[500]Theory of computation Shared memory algorithms
[300]Mathematics of computing Hypergraphs
[100]Theory of computation Graph algorithms analysis

§ INTRODUCTION

The balanced hypergraph partitioning problem asks for a partition of the node set of a hypergraph into a fixed number of disjoint blocks with bounded size
such that an objective function defined on the hyperedges is minimized. The two most prominent objective functions
are the edge cut and connectivity metric. The former counts the number of hyperedges connecting more than
one block, while the latter additionally considers the number of blocks spanned by each hyperedge.
The problem has gained attraction in the field of very-large-scale-integration (VLSI) design already in the 1960s <cit.>.
Since then, it has been widely adopted in many other areas, such as minimizing the communication volume in parallel scientific
simulations <cit.>, storage sharding in distributed
databases <cit.>, simulations of distributed quantum circuits <cit.>,
and as a branching strategy in satisfiability solvers <cit.>.

Unfortunately, balanced partitioning is NP-hard <cit.> and hard to
approximate <cit.>.
Thus, heuristic solutions are used in practice  – with
the multilevel scheme emerging as the most successful method to achieve high solution quality in a reasonable
amount of time <cit.>.
Figure <ref> illustrates this
technique, which consists of three phases.
First, the hypergraph is coarsened to obtain a hierarchy of successively smaller and structurally similar
approximations of the input hypergraph by contracting pairs or clusters of highly-connected nodes.
Once the hypergraph is small enough, an initial partition into $k$ blocks is computed.
Subsequently, the contractions are reverted level-by-level, and, on each level, local search heuristics are
used to improve the partition from the previous level.

The multilevel paradigm.

There is a diverse landscape of algorithms that implement the multilevel framework with different time-quality trade-offs, as illustrated in
Figure <ref>.

The plot shows two major shortcomings of existing solvers: (i) higher solution
quality comes at the cost of higher running times often by several orders of magnitude, and (ii) parallel algorithms do
not achieve the same solution quality as the best sequential systems because they use comparatively
weaker components that are easier to parallelize (with the exception of our new solver Mt-KaHyPar).

Historically, the parallel partitioning community has focused on algorithms
for the distributed-memory model, which turned out to be not well-suited for the fine-grained parallelism required
to effectively parallelize high-quality techniques. However, as the number of cores and main-memory
capacity in modern machines increases, we believe that the shared-memory model has become a viable alternative for
processing large (hyper)graphs and can be used for closing the quality gap between sequential and parallel partitioning algorithms.

Solution quality and running times of existing algorithms for hypergraph partitioning (left, connectivity metric) and
graph partitioning (right, edge cut metric).
For the $y$-values in the plot (solution quality), we compute the ratios of the objective values of an algorithm relative to the best value produced by any algorithm for each instance and
aggregate them using the harmonic mean (similarly for running times on the $x$-axis).

Markers on the lower left side are considered better.
We run each parallel algorithm using 10 threads.
Instances are restricted to more than two million edges/pins to make parallelism worthwile (see  and  in Section <ref>).
Partially transparent markers indicate solvers producing more than 15$\%$ infeasible partitions (either imbalanced or timeout).

Main Contributions

Figure <ref> highlights the main contribution of this work: A shared-memory multilevel algorithm (Mt-KaHyPar) that achieves
the same solution quality as the best sequential codes, while being faster than most of the relevant parallel algorithms in its fastest configuration.

In particular, our Mt-KaHyPar solvers occupy all points on the Pareto frontier for hypergraphs (left) as well as the middle segment for graphs (right).

This is achieved by implementing parallel formulations for the core techniques used in the best sequential algorithms

without compromises in solution quality.
Our coarsening algorithm contracts a clustering of highly-connected nodes on each level and is guided by the community structure of the hypergraph.
The clustering algorithm uses a less restrictive locking protocol than a previous approach <cit.> and resolves conflicting clustering decisions
on-the-fly. Initial partitioning is done via parallel recursive bipartitioning
and a portfolio solver, leveraging work-stealing to account for load imbalances. The key feature distinguishing Mt-KaHyPar from previous parallel systems are the substantially stronger local search
algorithms. We present the first fully-parallel implementation of the FM algorithm and
a parallel version of flow-based refinement.
For these algorithms, we propose several novel and easy-to-implement solutions to overcome some fundamental parallelization challenges such as, for example, techniques to
(re)compute correct gain values for concurrent node moves.

We also present several extensions of the core multilevel algorithm.
We devise the first parallel formulation of the $n$-level partitioning scheme – the most extreme instantiation of the multilevel technique –
contracting only a single node on each level.  Correspondingly, in each refinement step, only a single node is uncontracted followed by a highly-localized
search for improvements around the uncontracted node, leading to more fine-grained refinement and ultimately better solution quality in a single run.
Furthermore, we present a deterministic version of our multilevel algorithm.
This offers reproducible results and thus also stable results, whereas previous algorithms may have large variance from repeated runs.
Furthermore, some applications even require deterministic results or value them highly (e.g., VLSI design due to manual post-processing).

Moreover, we present data structure optimizations that speed up our algorithm by a factor of two when running on plain graphs instead of hypergraphs.

In our extensive experimental evaluation, we compare Mt-KaHyPar to $25$ different sequential and parallel
graph and hypergraph partitioners on over $800$ graphs and hypergraphs with up to 2 billion edges/pins.
As of today and to the best of our knowledge, this is the most comprehensive comparison of partitioning algorithms in the literature.
As a main result, the highest-quality configuration of Mt-KaHyPar
produces partitions that are on par with KaHyPar <cit.> – the best sequential hypergraph partitioner – while being
almost an order of magnitude faster with only ten threads. The fastest configuration of Mt-KaHyPar achieves a self-relative speedup
of $22.3$ with 64 threads and computes partitions that are $23\%$ better
than those of Zoltan <cit.> (distributed-memory), while being a factor of $2.72$ faster on average.
Out of all evaluated algorithms, KaFFPa <cit.> (sequential) computes slightly better solutions, while
KaMinPar <cit.> (shared-memory) is faster than Mt-KaHyPar.

The work presents the main results of several conference publications <cit.> and summarizes
the dissertations of Gottesbüren <cit.> and Heuer <cit.>. The added value of the paper is the detailed
overview of the overall framework that contains the highest-quality and one of the fastest algorithms for partitioning (hyper)graphs.

This paper puts particular focus on our multilevel partitioning algorithm which provides the best time-quality trade-off.
We describe the algorithm with a greater level of detail compared to the corresponding conference version <cit.>.
Furthermore, the previously mentioned optimizations for graph partitioning are unpublished. Another key contribution is the large experimental
evaluation, going beyond the scope of the individual publications by including graph partitioning, breaking down the running times of individual components, and including even more competing baseline algorithms.
We included almost all publicly available multilevel graph
and hypergraph partitioning algorithms to provide a comprehensive overview on the landscape of partitioning tools.

Section <ref> introduces basic notation and definitions used throughout this work. We then start
the algorithm description with a high-level overview of the multilevel partitioning algorithm in Section <ref>.
The following sections are structured according to the different phases of the multilevel scheme: Section <ref>
and <ref> describe the coarsening and initial partitioning algorithm, while we discuss
different concurrent gain (re)computation techniques and the implementation of the parallel FM and flow-based
refinement algorithms in Section <ref>–<ref>. In Section <ref>, we
present the parallelization of the $n$-level partitioning scheme, and conclude the algorithmic part with our
data structure optimizations for graph partitioning and the deterministic version of the multilevel algorithm
in Section <ref> and <ref>. We then turn to the experimental evaluation
in Section <ref>.

Here, we evaluate the solution quality and scalability of the different configurations of
Mt-KaHyPar, and compare them to existing partitioning algorithms.
Section <ref> concludes the work and presents directions for future research.

As this paper covers a wide range of partitioning techniques, we review relevant literature in the
corresponding sections. For a comprehensive overview on (hyper)graph partitioning, we refer the reader
to existing surveys <cit.>
and the literature overviews in the theses of Lars Gottesbüren <cit.>, Tobias Heuer <cit.>,
and Sebastian Schlag <cit.>.

§ PRELIMINARIES

A weighted hypergraph $H=(V,E,c,\omega)$ is defined as a set of  $n$ nodes $V$ and a set of $m$
hyperedges $E$ (also called nets)
with node weights $c:V \rightarrow \mathbb{R}_{>0}$ and net weights $\omega:E \rightarrow \mathbb{R}_{>0}$, where each net $e$ is
a subset of the node set $V$. The nodes of a net are called its pins.
We extend $c$ and $\omega$ to sets in a natural way, i.e., $c(U) :=\sum_{u \in U} c(u)$ and $\omega(F) :=\sum_{e \in F} \omega(e)$.
A node $u$ is incident to a net $e$ if $u \in e$.
$\incnets(u) := \{e \mid u \in e\}$ is the set of all incident nets of $u$.
The set $\neighbors(u) := \{ v \mid \exists e \in E : \{u,v\} \subseteq e\}$ denotes the neighbors of $u$.
Two nodes $u$ and $v$ are adjacent if $v \in \neighbors(u)$.
The degree of a node $u$ is $d(u) := |\mathrm{I}(u)|$.
The size $|e|$ of a net $e$ is the number of its pins.
Nets of size one are called single-pin nets.
We denote the number of pins of a hypergraph with $p := \sum_{e \in E} |e| = \sum_{v \in V} d(v)$.

We call two nets $e_i$ and $e_j$ identical if $e_i = e_j$.
Given a subset $V' \subset V$, the subhypergraph
$\subhypergraph{V'}$ is defined as $\subhypergraph{V'}:=(V', \{e \cap V' \mid e \in E : e \cap V' \neq \emptyset \}, c, \omega')$
where $\omega'(e \cap V')$ is the weight of hyperedge $e$ in $H$.
The bipartite graph representation $G_x := (V \cup E, E_x)$ <cit.>
of an unweighted hypergraph $H = (V,E)$ contains the nodes and nets of $H$ as node set and for each pin $u \in e$, we add an undirected edge $\{u,e\}$ to $E_x$.
More formally, $E_x := \{\{u,e\} \mid \exists e \in E: u \in e\}$.

Clusterings and Partitions
A clustering $\cluster = \{C_1, \ldots, C_l\}$ of a hypergraph $H = (V,E,c,\omega)$ is a partition of the node set $V$ into disjoint subsets.
A cluster $C_i$ is called a singleton cluster if $|C_i| = 1$.
A node contained in a singleton cluster is called unclustered.
A $k$-way partition of a hypergraph $H$ is a clustering into a predefined number of disjoint blocks
$\Partition = \{V_1, \ldots, V_k\}$.
A $2$-way partition is also called a bipartition.
We denote the block to which a node $u$ is assigned by $\nodeblock{u}$.

For each net $e$, $\conset(e) := \{V_i \mid  V_i \cap e \neq \emptyset\}$ denotes the connectivity set of $e$.
The connectivity $\con(e)$ of a net $e$ is $\con(e) := |\conset(e)|$.
A net is called a cut net if $\con(e) > 1$.
A node $u$ that is incident to at least one cut net is called boundary node.
The number of pins of a net $e$ in block $V_i$ is denoted by $\pinsinpart(e,V_i) := |e \cap V_i|$.
We refer to $\pinsinpart(e,V_i)$ as the pin count value for a net $e$ and block $V_i$.
The set $E(V_i,V_j) := \{e \in E \mid \{V_i,V_j\} \subseteq \conset(e)\}$ represents the
cut nets connecting block $V_i$ and $V_j$. Two blocks $V_i$ and $V_j$ are adjacent if $E(V_i,V_j) \neq \emptyset$.
The quotient graph $\quotientgraph := (\Partition, E_{\Partition} := \{(V_i, V_j) \mid E(V_i,V_j) \neq \emptyset\})$
contains an edge between all adjacent blocks.

The Balanced Hypergraph Partitioning Problem

The balanced hypergraph partitioning problem is to
find a $k$-way partition $\Partition$ of a hypergraph $H$ that minimizes an objective function defined on the hyperedges where
each block $V' \in \Partition$ satisfies the balance constraint: $c(V') \leq \balancedconstraint{\max} := (1+\varepsilon) \lceil \frac{c(V)}{k} \rceil$[The
$\lceil\cdot\rceil$ in this definition ensures that there is always a feasible solution for inputs with unit node weights. However, this does not hold for general weighted inputs as finding
a balanced solution for its own is an NP-hard problem <cit.>. There exists several alternative definitions <cit.>, but no commonly
accepted way how to deal with feasibility. In this work, we use the original definition since our benchmark instances are unweighted.]
for some imbalance ratio $\varepsilon \in (0,1)$.
If $\Partition$ satisfies the balance constraint, we call $\Partition$ $\varepsilon$-balanced or just say balanced or feasible when $\varepsilon$
is clear from the context.
For $k = 2$, we refer to the problem as the bipartitioning problem.
The two most prominent objective functions are the cut-net metric $\ocut := \sum_{e \in \cutnets} \omega(e)$ (also called edge cut metric for graph partitioning)
and connectivity metric $\ocon(\Partition) := \sum_{e \in \cutnets} (\lambda(e) - 1) \cdot \omega(e)$ (also called $(\lambda - 1)$-metric)
where $\cutnets$ denotes the set of all cut nets. The cut-net metric directly generalizes the edge cut metric from graphs to
hypergraphs and minimizes the weight of all cut hyperedges.
The connectivity metric additionally considers the number of blocks connected by a net and thus more accurately models
the communication volume for parallel computations <cit.> (e.g., for the parallel sparse matrix-vector multiplication).
The hypergraph partitioning problem is NP-hard for both objective functions <cit.>.

Recursive Bipartitioning vs Direct $k$-way Partitioning
A $k$-way partition of a hypergraph can be obtained either by recursive bipartitioning or
direct $k$-way partitioning. The former first computes a bipartition and then
calls the bipartitioning routine on both blocks recursively until the input hypergraph is divided into the desired number of blocks. The latter partitions
the hypergraph directly into $k$ blocks and applies $k$-way local search algorithms to improve the solution.

§ A BRIEF OVERVIEW OF THE PARTITIONING ALGORITHM

Algorithm <ref> shows the high-level structure of our multilevel partitioning algorithm.
While the pseudocode presented does not explicitly exhibit parallelism, it shows the algorithmic components for which we provide
parallel implementations.

The coarsening algorithm proceeds in rounds until the hypergraph is considered as small enough for initial partitioning.
In each round, we find a clustering of highly-connected nodes and subsequently contract the clustering in parallel.
The clustering algorithm iterates over the nodes in parallel and finds the best target cluster for a node according
to a rating function. Afterwards, the node joins its desired cluster for which we implement a novel locking protocol that detects
and resolves conflicting clustering decisions on-the-fly.

Initial partitioning is done via parallel
recursive bipartitioning using a novel work-stealing approach to account for load imbalances within the parallel bipartitioning
calls. To compute an initial bipartition, we use a portfolio of nine different bipartitioning techniques, which is run several
times in parallel. The best bipartition out of all runs is then used as initial solution.

In the uncoarsening phase, we project the partition onto the next hypergraph in the hierarchy
by assigning the nodes to the block of their corresponding constituent in the coarser representation.
Subsequently, we improve the partition using three different parallel refinement algorithms:
label propagation refinement (used in most of the existing parallel partitioning algorithms),
a highly-localized version of the FM algorithm (improves an existing implementation used in Mt-KaHIP <cit.>),
and a novel parallelization of flow-based refinement. The rationale behind the use of three different local search algorithms executed in this order
is that it allows for increasingly better solution quality at the cost of higher running times.

The following sections are structured according to the different phases of the multilevel scheme, and provide a more detailed
explanation of the different algorithmic components of Algorithm <ref>. In Section <ref>
and <ref>, we present our coarsening and initial partitioning algorithm. The
description of the uncoarsening phase is split into three separate sections: Section <ref> describes the
partition data structure and several concurrent gain (re)computation techniques, while Section <ref> and <ref> presents our parallel
FM and flow-based refinement algorithm.

Hypergraph $H = (V,E)$, number of blocks $k$
$k$-way partition $\Partition$ of $H$

The Multilevel Partitioning Algorithm

$H_1 \gets H\EndOfStatement \mathcal{H} \gets \langle H_1 \rangle \EndOfStatement n \gets 1$
$V_i$ has too many nodes
$\cluster \gets \FuncSty{ComputeClustering}(H_n)$
$H_{n + 1} \gets H_n.\FuncSty{Contract}(\cluster)\EndOfStatement \mathcal{H} \gets \mathcal{H} \cup \langle H_{n+1}\rangle \EndOfStatement \Increment n$

$\Partition \gets$ InitialPartition($H_n, k$)
$i = n - 1$ down to $1$
$\Partition \gets$ project $\Partition$ onto $H_i$
$\FuncSty{LabelPropagationRefinement}(H_i,\Partition)$ *[r]finds easy improvements by moving single nodes
$\FuncSty{FMRefinement}(H_i,\Partition)$   *[r]finds short and non-trivial move sets
$\FuncSty{FlowBasedRefinement}(H_i,\Partition)$ *[r]global optimization finding long and complex move sets

§ THE COARSENING PHASE

The goal of the coarsening phase is to find successively smaller and structurally similiar approximations of the input
hypergraph <cit.> such that initial partitioning can find a partition of high quality not significantly worse
than the partition that can be found on the input hypergraph <cit.>.
This can be achieved by grouping highly-connected nodes together and merging each group into a single node,
which can be done by computing either a matching <cit.> or clustering of the
nodes <cit.>.

The latter was shown to be more effective in reducing the
size of (hyper)graphs with highly-skewed node degree distributions <cit.> (e.g., social networks).
In the following, we present our parallel clustering-based coarsening algorithm that works similar to the shared-memory version
of PaToH's coarsening scheme <cit.>.
However, the algorithm of <cit.> excessively locks nodes when
evaluating the rating function and adding nodes to clusters. We therefore propose a less restrive locking protocol that completly
omits locking nodes when computing the best target cluster for a node.
Moreover, it detects and resolves conflicting clustering decisions on-the-fly, while previous approaches relied on a postprocessing step <cit.>.

§.§ The Clustering Algorithm

Our coarsening algorithm repeatedly finds a clustering $\cluster$ of the nodes and subsequently contracts it
until the hypergraph is small enough. We represent the clustering $\cluster$ using an array $\rep$ of size $n$. We then choose
one representative $v \in C$ for each cluster $C \in \cluster$ and store $\rep[u] = v$ for each node $u \in C$.
Initially, each node is unclustered (i.e., $\rep[u] = u$ for each node $u \in V$). The clustering algorithm
then iterates over the nodes in parallel and assigns each unclustered node to the best target cluster
according to a rating function, which we introduce in the subsequent paragraph.

Cluster Join Operation

A path (top-left) and cyclic conflict (top-middle), and a combination of both conflicts (top-right)
with their resolutions (bottom).

Once a node $u$ chooses its desired target cluster $C$ represented by a node $v$, we have to set $\rep[u] = v$.
Since several nodes can join clusters simultaneously, there may occur conflicts that must
be resolved. As illustrated in Figure <ref>, there are two types of conflicts:
path and cyclic conflicts. A path conflict involves several nodes $u_1, \ldots, u_l$ and occurs
when each node $u_i$ tries to join $u_{i+1}$. In a cyclic conflict, the last node $u_l$ additionally
tries to join $u_1$. It is also possible that a combination of both conflicts occurs, as illustrated in Figure <ref>
(right). We can resolve a path conflict when each node $u_i$ waits until $u_{i+1}$ has joined its desired cluster.
Afterwards, we can set $\rep[u_i] = \rep[u_{i+1}]$ to resolve the conflict. However, applying this resolution scheme
to cyclic conflicts would result in a deadlock. Therefore, the threads must agree on a cluster join operation that
breaks the cycle and reduces it to a path conflict.

Algorithm <ref> shows the pseudocode of our cluster join operation, which takes a node $u$ as input, and adds it to a cluster represented by a node $v$.
The algorithm associates each node with one of the following three states: unclustered, currently joining a cluster, or clustered.

Unclustered nodes ($\rep[u] = u$) can join clusters, while an already clustered node is not considered by the clustering algorithm
anymore and therefore its representative does not change. If a thread sets the state of a node $u$ from unclustered to joining via an
atomic compare-and-swap operation, it acquires exclusive ownership for modifying $\rep[u]$ and setting its state to clustered.
Thus, if we succeed in setting the state of $u$ and $v$ to joining or $v$ is already clustered, we can safely set $\rep[u] = \rep[v]$
(see Line <ref>–<ref>) since this guarantees that no other thread modifies $\rep[u]$ and $\rep[v]$.

Note that the representative of $v$ may have changed due to concurrent cluster join operations.
In that case, its representative is stored in $\rep[v]$. We therefore always set $\rep[u] = \rep[v]$ (instead of $\rep[u] = v$).

If another thread sets the state of $v$ to joining, we know that $v$ also tries to join a cluster.
To resolve the conflict, we spin in a busy-waiting loop until
the state of $v$ is updated to clustered (see Line <ref>), and then join its new cluster (path conflict).
In the busy-waiting loop, we additionally check if $u$ is part of a cycle of nodes trying to join each other.
To detect a cyclic conflict, each node writes its desired target cluster into a globally shared vector and checks if this induces a cycle.
If so, the node with the smallest ID in the cycle gets to join its desired cluster, thus breaking the cycle.

A node $u$ that wants to join $v$'s cluster

Cluster Join Operation

()  $\atomicfunc{compare-and-swap}(\statearray[u], \statevar{Unclustered}, \statevar{Joining})$

()  $\statearray[v] = \statevar{Clustered}$ or $\atomicfunc{compare-and-swap}(\statearray[v], \statevar{Unclustered}, \statevar{Joining})$
$\rep[u] \gets \rep[v]$

(*[f]Another thread tries to add $v$ to a cluster)
(*[f]busy-waiting loop)  $\statearray[v] = \statevar{Joining}$
cyclic conflict detected and $u$ is node with smallest ID in cycle
$\rep[u] \gets \rep[v]\EndOfStatement \statearray[u], \statearray[v] \gets \statevar{Clustered}\EndOfStatement \textbf{break}$

(*[f]resolves path conflicts)  $\statearray[u] = \statevar{Joining}$   $\rep[u] \gets \rep[v]$

$\statearray[u], \statearray[v] \gets \statevar{Clustered}$

Rating Function
A node $u$ joins the cluster $C$ maximizing the heavy-edge rating function

\[r(u,C) = \sum_{e \in \incnets(u) \cap \incnets(C)} \frac{\omega(e)}{|e| - 1}.\]

The rating function is commonly used in the partitioning literature <cit.> and prefers clusters connected to $u$ via
a large number of heavy nets with small size. We evaluate the rating function by iterating
over the incident nets $e \in \incnets(u)$ and aggregating the ratings to the representatives $\rep[v]$ of each pin $v \in e$
in a thread-local hash table.

Afterwards, we iterate over the aggregated ratings and determine the representative $\rep[v]$ that maximizes $r(u,\rep[v])$.
Ties are broken uniformly at random. Subsequently, we perform the cluster join operation that sets $\rep[u] = \rep[v]$.

To aggregate ratings, we use fixed-capacity linear probing hash tables with $2^{15}$ entries and resort to a larger hash table
if the fill ratio exceeds $\nicefrac{1}{3}$ of the capacity. This technique can considerably reduce the number of cache misses
since most neighborhoods are small in real-world hypergraphs.
We further note that the representative of a node
can change during the evaluation of the rating function since we do not lock the nodes.
However, it has already been shown that such conflicts rarely happen in practice <cit.>
and therefore have a negligible impact on the partitioning result.

Contraction Limit
We stop coarsening when the number of nodes in the smallest hypergraph reaches $160k$. This contraction limit
was chosen based on our prior research on sequential hypergraph partitioning <cit.>. In addition, we terminate
the clustering algorithm when the number of nodes would drop below $\frac{c(V)}{2.5}$ after the contraction step.
This prevents the coarsening process from reducing the size of the hypergraph too aggressively <cit.>.
Conversely, we also stop coarsening if the contraction step does not reduce the number of nodes by more than $1\%$, even if the
$160k$ node limit is not reached. This can happen since we enforce an upper weight limit $c_{\max}$

on the weight of the heaviest cluster (set to $\frac{c(V)}{160k}$ as in KaHyPar <cit.>),
which prevents highly-skewed node-weight distributions that would make it difficult for initial partitioning to find a balanced solution <cit.>.
When adding a node to a cluster $C \in \cluster$, we ensure that $c(C) \le c_{\max}$ by updating cluster weights via atomic fetch-and-add instructions.
If $c(C) > c_{\max}$ after the update, we reject the corresponding cluster join operation and revert the cluster weight update.
The cluster weight limit can lead to coarsening passes that do not sufficiently reduce the size of the hypergraph.

§.§ The Contraction Algorithm

The hypergraph data structure stores the incident nets $I(u)$ of each node $u \in V$ and the pin-lists of each net $e \in E$
using two adjacency arrays. Each node $u$ and net $e$ additionally
stores its weight $c(u)$ and $\omega(e)$.
Contracting a clustering $\cluster = \{C_1,\ldots,C_l\}$ replaces each
cluster $C_i$ with one supernode $u_i$ with weight $c(u_i) = \sum_{v \in C_i} c(v)$. For each net $e \in E$,
we replace each pin $v \in e$ with the node $u_i$ representing the cluster $C_i$ in which $v$ is contained ($\rep[v] = u_i$).
After the replacement, multiple occurrences of the same supernode in a net are discarded.

Our contraction algorithm consists of several simple, easily parallelizable operations including remapping node IDs to a consecutive range,
aggregating cluster weights and degrees using atomic fetch-and-add instructions,
eliminating duplicated entries in pin-lists, and using parallel prefix sum operations to construct the adjacency arrays of the contracted hypergraph.
As these steps are rather low level, we refer the reader to Ref. <cit.> for more details.

A challenging aspect is removing duplicates from the set of nets.
We identify groups of identical nets and remove all but one representative per group to which we assign their aggregate weight.
This can reduce the number of pins significantly and therefore accelerates the other algorithmic components.
A simple algorithm is to perform pair-wise comparisons between all nets, which is however too expensive in practice.
To this end, we parallelize the InrSrt algorithm of Aykanat et al. <cit.> for identical net detection.
It uses fingerprints $f(e) := \sum_{v \in e} v^2$ to eliminate unnecessary pairwise comparisons between nets, by grouping nets with equal fingerprints via sorting.
Nets with different fingerprints or different sizes cannot be identical.
We distribute the fingerprints and their associated nets to the threads using a hash function.
Each thread sorts the nets by their fingerprint and size, and then performs pairwise comparisons on the subranges of potentially identical nets.
We aggregate the weights of identical nets at a representative and mark the others as invalid in a bitset.
A parallel prefix sum over the bitset maps the hyperedge IDs to a consecutive range in the contracted hypergraph.
Note that we also remove nets that contain only a single pin since they do not contribute to the cut.

§.§ Community-Aware Coarsening

A popular approach to improve an existing $k$-way partition $\Partition$ is the iterated multilevel cycle technique <cit.> (also called V-cycle).
In the coarsening phase, the algorithm forbids contractions between nodes that are not in the same block in $\Partition$, thus preserving the already identified cut structure. While the technique can be effective, using it as a postprocessing step in a multilevel algorithm almost doubles the running time.
As a more lightweight alternative, Heuer and Schlag <cit.> proposed using a clustering of the nodes
computed via a community detection algorithm instead of an existing $k$-way partition.
Community detection still captures the sparse cut patterns that are often found in good $k$-way partitions.

The authors showed that this substantially improves the
quality of both the initial and the final partition, and only slightly increases the running time of the overall algorithm.

We also integrate the approach into our partitioning algorithm. We run the algorithm
as a preprocessing step before the coarsening phase and then use the clustering to restrict contractions
to nodes that belong to the same cluster.
The algorithm consists of two steps: transforming the hypergraph into its bipartite graph representation and then running
the parallel Louvain method of Staudt and Meyerhenke <cit.> for modularity
maximization, a widely used objective function for community detection <cit.>.

§ THE INITIAL PARTITIONING PHASE

Partitioning algorithms based on the direct $k$-way partitioning scheme often use multilevel recursive
bipartitioning to obtain an initial $k$-way partition <cit.>, as this leads to partitions with significantly better solution quality than using flat (non-multilevel)
$k$-way partitioning methods.

Many parallel partitioners run sequential initial partitioning algorithms in parallel <cit.>. However, the sequential calls can become a bottleneck when the smallest
hypergraph is still large. A more scalable approach parallelizes the recursive calls after each bipartitioning operation <cit.>.
The common approach is to statically split the thread pool along with the subproblems. Since this can lead to load imbalance when processing hypergraphs with unequal densities
in the recursive partitioning calls ,we instead generate tasks that can be dynamically load balanced using work stealing.

Parallel Recursive Bipartitioning
We compute initial $k$-way partitions via parallel recursive bipartitioning using
Algorithm <ref> initialized with $k = 2$ (without flow-based refinement). For the bipartitioning case,
we replace the initial partitioning call with a portfolio of bipartitioning techniques.

Once we obtain a bipartition $\Partition = \{V_1,V_2\}$ of the input hypergraph $H$, we extract the subhypergraphs $\subhypergraph{V_1}$
and $\subhypergraph{V_2}$
and recurse on both in parallel by partitioning $\subhypergraph{V_1}$ into $\lceil \frac{k}{2} \rceil$
and $\subhypergraph{V_2}$ into $\lfloor \frac{k}{2} \rfloor$ blocks.
We ensure that the final $k$-way partition
obtained via recursive bipartitioning is $\varepsilon$-balanced by adapting the imbalance ratio for each bipartition individually <cit.>.
Let $\subhypergraph{V'}$ be a subhypergraph that should be recursively partitioned into $k' \le k$ blocks. Then,

\begin{equation}
\label{eq:adaptive_epsilon}
\varepsilon' := \biggl( (1 + \varepsilon) \frac{c(V)}{k} \cdot \frac{k'}{c(V')} \biggr)^{\frac{1}{\lceil \log_2{k'} \rceil}} - 1
\end{equation}

is the imbalance ratio used for the bipartition of $H_{V'}$. If each bipartition is $\varepsilon'$-balanced, then it is guaranteed
that the final $k$-way partition is $\varepsilon$-balanced <cit.>.

Portfolio-Based Bipartitioning
We implemented the same portfolio of initial bipartitioning techniques as in KaHyPar <cit.>,
including seven different variants of (greedy) hypergraph growing <cit.>,
random assignment <cit.>, and label propagation initial
partitioning <cit.>.

We refer the reader to Ref. <cit.> for more details on their implementation.
We run each algorithm independently in parallel for at least $5$ and at most $20$ times.
After $5$ runs, we only run an algorithm again if it is likely to improve the best solution $\Partition^*$ found so far.
We estimate this based on the arithmetic mean $\mu$ and standard deviation $\sigma$ of the connectivity values achieved by that algorithm so far, using the $95\%$ rule.
Assuming the connectivity values follow a normal distribution, roughly $95\%$ of the runs will fall between $\mu - 2 \sigma$ and $\mu + 2 \sigma$.
If $\mu - 2 \sigma > \ocon(\Partition^*)$, we do not run the algorithm again.
Additionally, we refine each bipartition using sequential $2$-way FM refinement <cit.>.
We continue uncoarsening using the bipartition with the best connectivity value. In case of ties, we prefer the bipartition
with the best balance.

§ GAIN COMPUTATION TECHNIQUES

Local search algorithms greedily move nodes to different blocks according to a gain value.
The gain value reflects the change in the objective function for a particular node move. For the connectivity metric,
the gain $\gain{u}{t}$ of moving a node $u$ to a target block $V_t$

can be expressed as follows:

\[\gain{u}{t} := \omega(\{e \in \incnets(u) \mid \pinsinpart(e, \nodeblock{u}) = 1 \}) - \omega(\{e \in \incnets(u) \mid \pinsinpart(e, V_t) = 0 \}). \]

Moving node $u$ to block $V_t$ decreases the connectivity of all nets by one for which $u$ is the last remaining pin in its current block $\nodeblock{u}$.
Conversely, the move increases the connectivity of all nets $e \in \incnets(u)$ by one for which no pin $v \in e$ is assigned to the target block $V_t$.

Example of a move conflict when two nodes are moved simultaneously. Both threads assume that the individual node moves
removes the edge with weight $5$ from the cut. However, the edge is still cut after moving both nodes and
the edges with weight $2$ become cut edges.

To achieve meaningful speedups, parallel refinement algorithms need to move nodes concurrently.
The actual gain of a node move can change between the time it is initially calculated and the time it is applied to the partition, due to concurrent node moves in its neighborhood <cit.>.
As a consequence, two concurrent node moves can worsen the connectivity metric, even if their individual gains suggested an improvement, as illustrated in Figure <ref>.
Thus, correctly calculating gains is a fundamental challenge for parallel refinement algorithms.

These conflicts occur when two adjacent nodes change their blocks simultaneously. Common remedies include computing a node coloring and only moving nodes of the same color at a time <cit.>, scheduling $2$-way refinement
algorithms on block pairs that form a matching in the quotient graph in parallel <cit.>,
allowing only node moves from a block $V_s$ to $V_t$ if $s < t$ <cit.> (and vice versa in a second phase),
or following an optimistic strategy assuming that conflicts happen rarely in practice <cit.>.

The presented approaches still allow all individual node moves,
but combining arbitrary moves into a single move sequence might be not always possible.
This is problematic for parallelizing local search techniques as their sequential counterparts often identify a
set of moves that only yield an improvement if moved together.

While ignoring search conflicts appears to be the preferred approach, their impact on solution quality is unpredictable and deserves further consideration.

We therefore contribute several parallel gain computation techniques to compute accurate gain values and detect conflicts
between moves without restricting possible moves. We present a technique named attributed gains to double-check the gain of a node move
in Section <ref>, a concurrent gain table to accelerate gain calculations and communicate updates between threads in Section <ref>,
and a novel parallel algorithm for recomputing exact gains of a sequence of node moves in Section <ref>.
These techniques build on our concurrent partition data structure which we describe in the next section in more detail.

§.§ The Partition Data Structure

Our partition data structure stores and maintains the block assigments $\Partition$, the block weights $c(V_i)$,
the pin count values $\pinsinpart(e,V_i)$, and
connectivity sets $\conset(e)$ for each net $e \in E$ and block $V_i \in \Partition$.

The Move Node Operation

Algorithm <ref> shows the updates to the partition data structure when moving a node $u$ from its source block $V_s$ to a target block $V_t$.
We only perform a node move if it does not violate the balance constraint, which we ensure by adding the weight
of node $u$ to the weight of block $V_t$ via an atomic fetch-and-add instruction.
If the node move is feasible, we update the block assignment of node $u$ to block $V_t$ and subtract the node weight $u$ from
its previous block $V_s$.
If the move is infeasible, we subtract the weight again and reject the move.

A node $u$ that should be moved from its source block $V_s$ to a target block $V_t$
Attributed gain value $\attrgain$

The Move Node Operation

$c_t \gets \atomicfunc{fetch-and-add}(c(V_t), c(u))$
(*[f]Revert block weight update if balance constraint violated)  $c_t + c(u) > \balancedconstraint{\max}$
$c(V_t) \minuseqatomic c(u)$0

$\nodeblock{u} \gets V_t$$c(V_s) \minuseqatomic c(u)$$\attrgain \gets 0$
$e \in I(u)$
$\FuncSty{lock}(e)$$\Phi_s \gets \Decrement\pinsinpart(e,V_s)$$\Phi_t \gets \Increment\pinsinpart(e,V_t)$$\FuncSty{unlock}(e)$
*[l]Update connectivity set $\conset(e)$ and attributed gain value $\attrgain$
()  $\Phi_s = 0$   $\conset(e) \gets \conset(e) \setminus \{ V_s \}$$\attrgain \pluseq \omega(e)$
()  $\Phi_t = 1$   $\conset(e) \gets \conset(e) \cup \{ V_t \}$$\attrgain \minuseq \omega(e)$
$\FuncSty{UpdateGainTable}(e, \Phi_s, \Phi_t)$ *[r]see Section <ref>

Data Layout

The size of a pin count value is bounded by the size of the largest hyperedge.
To save memory, we use a packed representation with $\lceil \log({\max_{e \in E} |e|}) \rceil$ bits per entry
for the $\pinsinpart(e,V_i)$ values. Furthermore, we use a bitset of size $k$ to store the connectivity set
$\conset(e)$ of each hyperedge $e \in E$. We iterate over the connectivity set $\conset(e)$ by taking a snapshot
of its bitset and then use count-leading-zeroes instructions.
We compute the connectivity $\con(e) = |\conset(e)|$ of a hyperedge $e$ using pop-count instructions
(counts the number of $1$-bits in a machine word).
To add or remove a block from the connectivity set, we flip the corresponding bit using an atomic xor operation.
The move node operation can be made lock-free by updating $\pinsinpart(e, V_i)$ with atomic fetch-and-add instructions, but this requires one machine word per value.
We therefore use a spin-lock for each net $e$ due to the packed representation.

Attributed Gains

As the gain value of a node move can change between its initial calculation and actual execution due to concurrent node moves in its neighborhood,
we additionally compute an attributed gain value for each move based on the atomic updates
of the pin count values $\pinsinpart(e,V_s)$ and $\pinsinpart(e,V_t)$ in Line <ref>
of Algorithm <ref>.
We attribute a connectivity decrease by $\omega(e)$ to the move that reduces $\pinsinpart(e, V_s)$ to zero (see Line <ref>)
and an increase by $\omega(e)$ for increasing $\pinsinpart(e,V_t)$ to one (see Line <ref>).

Since we do not lock all incident nets $e \in \incnets(u)$ before moving a node $u$,
there is no guarantee on the order in which concurrent moves perform the pin count updates.
Hence, this scheme may distribute the connectivity reductions to different threads,
but the sum of the attributed gains of all node moves equals the overall connectivity reduction <cit.>.

Attributed Gains for Label Propagation Refinement

The most widely used refinement technique in parallel partitioning
algorithms is label propagation <cit.>.
The algorithm works in rounds. In each round, it iterates over all nodes in parallel, and whenever it
visits a node $u$, it moves it to the block $V_t$ maximizing its move gain $\gain{u}{t}$ (respecting the balance constraint).
The algorithm only performs moves with positive gain and therefore cannot escape from local optima.
However, we use it in our partitioning algorithm to find all simple node moves such that our more advanced
refinement techniques can focus on finding non-trivial improvements (for more technical details on its implementation,
see Ref. <cit.>).

Since the label propagation algorithm performs only positive gain moves, we immediately revert a node move
if it has negative attributed gain.
Note that reverting such a node move
does not guarantee to improve the connectivity metric again as other concurrent node moves may have
changed the pin count values of the corresponding nets in the meantime. However, reverting them directly after detection
decreases the likelihood of such conflicts.
Furthermore, we use attributed gains to track the value of the connectivity metric instead of recomputing it after each round.

§.§ The Gain Table

For our FM algorithm, we use a gain table which stores and maintains the gain values for all possible moves.
This enables repeatedly looking up gains in $O(1)$ time and is a globalized way of updating the gains of nodes owned by other threads.

Gain tables are not a new idea <cit.> but have gone out of fashion
due to their memory requirements <cit.>.
To the best of our knowledge, our introduction of parallel gain tables is novel.

We use atomic fetch-and-add instructions to update the gains as soon as nodes are moved.
Updates on some nodes become visible while the overall update procedure is still in flight.
Therefore, updates trickle in over time, and some outdated or inconsistent values may be read by other threads.
Still, with concurrent node moves this is the most accurate we can be.

Recall that the gain $\gain{u}{t}$ of moving a node $u$ to a target block $V_t$ can be expressed
as follows:

\[\gain{u}{t} := \omega(\{e \in \incnets(u) \mid \pinsinpart(e, \nodeblock{u}) = 1 \}) - \omega(\{e \in \incnets(u) \mid \pinsinpart(e, V_t) = 0 \}). \]

The first term $b(u) \coloneq \omega(\{e \in \incnets(u) \mid \pinsinpart(e, \nodeblock{u}) = 1 \})$ is the benefit of moving $u$ out of its block.
Conversely, the term $p(u, V_t) \coloneq \omega(\{e \in \incnets(u) \mid \pinsinpart(e, V_t) = 0 \})$ is the penalty for moving $u$ into $V_t$.

Update Rules

Instead of storing $\gain{u}{i}$, we store $b(u)$ and $p(u, V_i)$ separately for each node $u$, so that changes to $b(u)$ only require one update,
instead of updates to $k$ gain values.
This approach uses $(k+1)n$ memory words in total.
For each net $e \in I(u)$, we update $b(u)$ and $p(u,V_i)$ using atomic fetch-and-add instructions as follows.

If  (e, V_s) = 0  then  ∀v ∈e                do      p(v, V_s)          ω(e)

If  (e, V_s) = 1  then  ∀v ∈e ∩V_s       do      b(v)               ω(e)

If  (e, V_t) = 1  then  ∀v ∈e                do      p(v, V_t)          ω(e)

If  (e, V_t) = 2  then  ∀v ∈e ∩V_t       do      b(v)               ω(e)

The update conditions implement the UpdateGainTable procedure from Line <ref> in Algorithm <ref>.

Benefit Pecularities

There is a race condition on $\nodeblock{v}$ in the check $\nodeblock{v} = V_s$ (case 2) or $\nodeblock{v} = V_t$ (case 4).
When $\nodeblock{v}$ changes, we may perform a benefit update on $v$ that was also intended for a different pin of $e$ in the new $\nodeblock{v}$.
The penalty values are not affected since they are independent of the pin's current block.
Our FM algorithm is organized in rounds in which each node can be moved at most once.
Therefore, once $u$ gets moved, we do not read $b(u)$ for the rest of the round.
Due to the race condition it may still be updated, which is why we recalculate $b(u)$ after the round is finished instead of recalculating $b(u)$ for the new block immediately after the move.
We note that it is possible to correctly update benefits by using $k$ benefit values per node <cit.>.

Correctness and Complexity of Gain Updates

In the following, we prove that once all updates for a given set of moves are completed and no further moves are performed, the gain values are correct.

After performing all gain updates associated with a set of moves $M$ in parallel, each unmoved node $v \in V \setminus M$ has correct $b(v)$, and each $v \in V$ has correct $p(v, V_i)$ terms.

First, we note that the updates are correct in the sequential setting <cit.>.
Due to the atomic consistency of pin-count and gain updates, it suffices to prove correctness for arbitrary linearized (sequential) orders of updates.
The remaining difficulty is that different orders may yield different intermediate values.
However, due to commutativity we arrive at the same final $\pinsinpart(e,V_i)$ values.
Thus, it suffices to argue that gain updates triggered by $\pinsinpart(e, V_i) \pluseq 1$ cancel out those triggered by $\pinsinpart(e, V_i) \minuseq 1$.
This statement holds, as case 1 and 3 are complimentary, as well as case 2 and 4.
Therefore, the final $p(v, V_i)$ and $b(v)$ values only depend on the final $\pinsinpart(e, V_i)$ values.

The work of gain updates for moving all nodes once is $\Oh{\sum_{e \in E}|e| \cdot \min(k, |e|)} = \Oh{kp}$.

The core to the argument is that each of the update cases is only triggered a constant number of times per hyperedge and block <cit.>, and costs $\mathcal{O}(|e|)$ work per update.
This hinges on moving each node at most once.

Note that due to the $\min(k, |e|)$ term, this bound matches the $\mathcal{O}(m)$ bound on plain graphs.
On real-word hypergraphs, we observed work much closer to $\Oh{p}$ since most nets have small size or few pins per block.

§.§ The Parallel Gain Recalculation Algorithm

We now propose a parallel algorithm to recompute exact gain values of a sequence of node moves $M = \langle m_1, \dots, m_l \rangle$
if they are supposed to be performed in this order. Each move $m_i \in M$ is of the form $m_i = (u,V_s,V_t)$, which means
that node $u$ is moved from block $V_s$ to $V_t$.

Again, we assume that each node is moved at most once.
Recall that a move of a node $u$ from block $V_s$ to $V_t$ decreases the connectivity of a hyperedge $e$,
if $\pinsinpart(e,V_s)$ decreases to zero. Conversely, it increases the connectivity if $\pinsinpart(e,V_t)$
increases to one. The idea of the following algorithm is to iterate over the hyperedges in parallel, and
identify the node moves in $M$ that increase or decrease the connectivity of a hyperedge using Algorithm <ref>.

Consider a hyperedge $e$ and a block $V_i \in \Partition$. The first observation is that if we move a pin $v \in e$ to $V_i$,
then $\pinsinpart(e,V_i)$ cannot decrease to zero anymore since each node is moved at most once. In order to decrease $\pinsinpart(e,V_i)$
to zero, we have to move all pins $u \in e \cap V_i$ out of block $V_i$ before we move the first pin $v \in e \setminus V_i$ to block $V_i$.
In this case, the last pin $u \in e$ moved out of block $V_i$ decreases the connectivity of $e$ and the first pin $v \in e$
moved to block $V_i$ increases its connectivity again. Thus, we can decide whether or not a move increases or decreases the
connectivity of a hyperedge by simply comparing the indices of the node moves in $M$, which were last moved out and first moved
to a particular block. Additionally, we need to know if the move sequence $M$ moves all pins out of block $V_i$.
To do so, we count the number of non-moved pins $v \in e$ in each block. If the number of non-moved pins is zero for a block $V_i$, then
either $\pinsinpart(e,V_i)$ was zero before, or the move sequence $M$ moved all nodes out of block $V_i$.

Algorithm <ref> shows the pseudocode that identifies the node moves in $M$ that increase or decrease the connectivity
of a hyperedge $e$. The algorithm uses two loops, both iterating over the pins of hyperedge $e$.
The first loop computes the indices of the node moves
that first moved to and last moved out of each block $V_i \in \Partition$ (see Line <ref>), in addition to the number of pins in $e$ that were not moved (see Line <ref>).

The second loop then decides for each moved pin $u \in e$ whether or not it increases or decreases the connectivity of hyperedge $e$ by evaluating
the conditions shown in Lines <ref> and <ref>.
Let $m_i := (u,V_s,V_t)$ be the corresponding node move of pin $u \in e$ in $M$.
If $M$ moves all nodes out of block $V_s$ ($\matharray{non\_moved}[V_s] = 0$) and $u$ is the last pin moved out
of block $V_s$ ($\matharray{last\_out}[V_s] = i$), while the first move that moves a pin into block $V_s$ happens strictly after
$m_i$ ($i < \matharray{first\_in}[V_s]$), then $m_i$ reduces the connectivity metric by $\omega(e)$.
Conversely, if $M$ moves all nodes out of block $V_t$ ($\matharray{non\_moved}[V_t] = 0$) and $u$ is the first pin
moved into block $V_t$ ($\matharray{first\_in}[V_t] = i$), while the last move that moves a pin out of block $V_t$ happens strictly
before $m_i$ ($i > \matharray{last\_out}[V_t]$), then $m_i$ increases the connectivity metric by $\omega(e)$.
Since we run the algorithm for each hyperedge in parallel, several threads can modify the gain value $g_i$ of a node move $m_i$ simultaneously.
We therefore use atomic fetch-and-add instructions (see Line <ref> and <ref>).

To further reduce the complexity of the algorithm, we only process hyperedges containing moved nodes. To do so, we iterate
over the node moves in $M$ in parallel and run Algorithm <ref> only for incident edges
of moved nodes. We mark already processed hyperedges in a shared bitset using atomic test-and-set instructions.

Hyperedge $e$, a sequence of node moves $M = \langle m_1, \dots, m_l \rangle$ and
a shared gain vector $\mathcal{G} = \langle g_1, \dots, g_l \rangle$ representing the recalculated gain values

Parallel Gain Recalculation

$\matharray{first\_in} \gets [\infty,\ldots,\infty]\EndOfStatement \matharray{last\_out} \gets [-\infty,\ldots,-\infty]$ *[r]Arrays of size $k$
$\matharray{non\_moved} \gets [0,\ldots,0]$ *[r]Array of size $k$
$u \in e$
(*[f]moved nodes are marked in a bitset)  $u$ was moved
$m_i := (u,V_s,V_t) \gets$ find corresponding move in $M$
$\matharray{last\_out}[V_s] \gets \max(i, \matharray{last\_out}[V_s])$$\matharray{first\_in}[V_t] \gets \min(i, \matharray{first\_in}[V_t])$

$u \in e$
$u$ was moved
$m_i := (u,V_s,V_t) \gets$ find corresponding move in $M$
()  $\matharray{last\_out}[V_s] = i~\land~i < \matharray{first\_in}[V_s]~\land~\matharray{non\_moved}[V_s] = 0$
$g_i \pluseqatomic \omega(e)$

()  $\matharray{first\_in}[V_t] = i~\land~i > \matharray{last\_out}[V_t]~\land~\matharray{non\_moved}[V_t] = 0$
$g_i \minuseqatomic \omega(e)$

§ THE FIDUCCIA-MATTHEYSES ALGORITHM

The Fiduccia-Mattheyses (FM) algorithm <cit.> is the most widely used local search algorithm in sequential partitioning algorithms.
Most of the existing variants insert all possible moves or only the highest gain move for each boundary node into a PQ (boundary FM)
and then perform the following two steps:
(i) repeatedly perform the highest gain move subject to the balance constraint, followed by (ii) reverting moves back to the prefix
with the highest cumulative gain in the sequence of performed moves.

The revert is necessary, since moves with negative gains are allowed, so the algorithm is able to escape from local minima.
Unfortunately, calculating the same move sequence as FM is P-hard <cit.>, i.e., it is unlikely that a parallel algorithm with poly-log depth exists.

Sanders and Schulz <cit.> proposed a relaxed version that inserts only the highest gain move for a single seed node into a PQ and then
gradually expands around the node by claiming neighbors of moved node (localized FM). The algorithm not only produces better solutions
than boundary FM <cit.>, it is also highly amenable to parallelization as multiple FM searches can
run in parallel, each starting from a different seed node. In the following, we present our parallel implementation of the localized FM algorithm, and
discuss its main differences to an existing parallelization <cit.>.

The Parallel $k$-Way FM Algorithm
Algorithm <ref> shows the pseudocode of our parallel FM algorithm.
The algorithm proceeds in rounds, and each round starts with inserting all boundary nodes into a globally
shared task queue $Q$.
The threads then poll a fixed number of nodes ($= 25$) from $Q$ that they use as seed nodes for the
localized FM searches, which expand to neighbors of moved nodes.

The searches are non-overlapping, i.e., threads acquire exclusive ownership of nodes,

while hyperedges can touch multiple searches.
Node moves performed by the different searches are not visible to other threads, as they are performed locally
using thread-local hash tables.
However, once a thread finds an improvement, it immediately applies it to the global partition.
The local moves are atomically appended to a global move sequence (using one atomic fetch-and-add for all local moves).
We repeatedly start localized FM searches until the task queue is empty.
Note that we initialize the searches with multiple seed nodes instead of a single node as this substantially
accelerates the algorithm in practice without sacrifices in solution quality.

Once the task queue is empty, we proceed to the second phase, where we recalculate the gains of the global move sequence
(see Section <ref>) and then use a parallel prefix
sum and reduce operation on the recomputed gain values to identify and revert to the best seen solution.

We perform multiple rounds until a maximum number is reached or the connectivity
metric is not improved.

Localized $k$-Way FM Search

The localized FM search uses a single PQ storing the move with the highest gain for each inserted node.
We initialize the PQ with several seed nodes
and use the gain table to compute the initial best move for each node (see Line <ref>).
Then, we repeatedly select the move with the highest gain
and apply it to a thread-local partition $\DeltaPartition$.
Changes on $\DeltaPartition$ are not visible to other threads for now. However, we apply the move sequence
to the global partition $\Partition$ as soon as we find an improvement (see Line <ref>), then triggering gain updates in the global gain table.

When we move a node $u$ locally, we collect the nets $e \in \incnets(u)$ affected
by gain updates. We use them to update the gain values of nodes in the PQ –
combining global gain table and $\DeltaPartition$ data, thus gradually infusing updates from other threads into the search
– and expand the search to neighbors of moved nodes.
A localized search terminates when the PQ becomes empty or the adaptive stopping rule of Osipov and Sanders <cit.>
is triggered. The stopping rule assumes that the observed gain values follow a normal distribution and terminates a search
when it becomes unlikely to find further improvements.
We release the ownership of non-moved nodes at the end such that other searches
can acquire them again. We do not release the ownership of moved nodes to
ensure that each node is moved at most once during an FM pass.

Parallel $k$-Way FM Algorithm

() $\FuncSty{FMRefinement}($Hypergraph $H,$ $k$-way partition $\Partition)$
improvement found and maximum number of rounds not reached
$Q \gets$ initialize task queue with all boundary nodes
$Q$ not empty
$V_{\text{\footnotesize seed}} \gets$ poll $25$ seed nodes from $Q$
$\FuncSty{LocalizedFMRefinement}(H,\Partition,V_{\text{\footnotesize seed}})$

recompute gains of global move sequence and revert to best prefix *[r]see Section <ref>

$\FuncSty{LocalizedFMRefinement}(H,\Partition,V_{\text{\footnotesize seed}})$
(*[f]Initialize PQ with seed nodes)  $u \in V_{\text{\footnotesize seed}}$
$(g_t,V_t) \gets \FuncSty{ComputeMaxGainMove}(u)$PQ.$\FuncSty{Insert}(u, g_t, V_t)$

$\Delta \gets 0\EndOfStatement M \gets \emptyset$
PQ not empty and search should continue
$(u,g_t,V_t) \gets$ PQ.$\FuncSty{PopMaxGainMove}()$
move $u$ to $V_t$ in thread-local partition $\DeltaPartition$ with gain table update
$\Delta \gets \Delta + g_t\EndOfStatement M \gets M \cup \{(u,V_t)\}$
$\Delta > 0$ or ($\Delta = 0$ and move improved balance)

apply move sequence $M$ to global partition $\Partition$
$\Delta \gets 0\EndOfStatement M \gets \emptyset\EndOfStatement \DeltaPartition \gets \Partition$

()  all nets $e \in I(u)$ involved in a gain update
$v \in e$
$v$ is not marked
()  PQ.$\FuncSty{Contains}(v)$   update gain of $v$ in PQ
()  try to acquire node $v$
$(g_t,V_t) \gets \FuncSty{ComputeMaxGainMove}(v)$PQ.$\FuncSty{Insert}(v, g_t, V_t)$

mark $v$

unmark all nodes

We explicitly allow moves with negative gains, which will worsen the solution quality intermediately.
At the end of each localized search, we thus revert back to the best seen solution.
If we directly applied moves to the global partition, other searches could base their decisions on states that will later be reverted.
Therefore, we apply node moves to a thread-local partition $\DeltaPartition$ first and only perform them on the global partition if they lead to an improvement, i.e., will not be reverted for now.

The thread-local partition $\DeltaPartition$ stores changes relative to the global partition in a set of hash tables.
For example, we compute the weight of a block $V_i$ by calculating $c(V_i) + \Delta c(V_i)$ where $c(V_i)$
is the weight of block $V_i$ stored in the global partition data structure and $\Delta c(V_i)$ is the weight
of all nodes that locally moved to block $V_i$ minus the weight of nodes that moved out of block $V_i$. We maintain
the block ID, pin count values, as well as benefit and penalty terms of the gain table analogously.

Applying a move sequence to the global partition makes it immediately visible to the searches on other threads.
Since $\DeltaPartition$ stores local changes relative to the global partition, the block weights and pin count values
are still correct. However, some gain values may be incorrect since the gain table updates on the global partition
do not consider moves performed locally.
This is only a small issue since thread-local deltas are cleared after applying the moves to the global partition.
In practice, the scheme drastically reduces conflicts[We have found that the recomputed gain values of
the global move sequence match the observed gain values during the localized FM searches in most cases.].
Another reason for applying moves as soon as possible is to keep the memory footprint of the hash tables small.
The overall peak memory incurred by thread-local partition data is small, because the memory is proportional to the number of moves,
and long-running searches must find improvements to keep going.

Differences to Mt-KaHIP
The FM implementation in Mt-KaHIP <cit.> performs node moves only locally, which are therefore not
visible to other threads. At the end of an FM pass, the move sequences found by the different searches are concatenated to a global
move sequence, for which gains are recomputed sequentially. We improved the algorithm by making improvements immediately
visible to other threads using the thread-local partition and gain table data structure, leading to more accurate gain values.
Moreover, we removed the last sequential part of the algorithm with our parallel gain recomputation technique.

§ FLOW-BASED REFINEMENT

A major shortcoming of move-based local search algorithms is that they greedily move nodes to other
blocks based on a gain value considering only the block assignment of adjacent nodes.
Thus, the decision to apply a move depends only on local information, which may not be sufficient to find some non-trivial improvements <cit.>.
Maximum flows overcome this limitation by deriving a minimum cut
seperating two nodes <cit.> and therefore have a more global view on the partitioning problem.

Although it seems natural to use them as local search strategy in partitioning algorithms, maximum flows were long perceived as computationally expensive and it was unclear how to derive balanced partitions <cit.>.
This changed over the last two decades as flow-based refinement techniques were successfully implemented in the highest-quality
sequential graph and hypergraph partitioning algorithms <cit.>.
Today it is considered to be the most powerful improvement heuristics for (hyper)graph partitioning.
However, since they come at the cost of substantially higher running times, they can be impractical for partitioning very large hypergraphs.

Algorithm Overview

In this section, we present the first parallel formulation of the sequential flow-based refinement approach
used in KaHyPar <cit.>.
The high-level pseudocode of the algorithm is outlined in Algorithm <ref>.
Flow-based refinement works on bipartitions and can be scheduled on different block pairs to improve $k$-way partitions <cit.>.
We therefore start with a parallel scheduling scheme of adjacent block pairs based on the quotient graph in
Section <ref> (see Line <ref> and <ref>).
In Section <ref>, we describe the flow network construction algorithm that extracts a subhypergraph induced by a region $B \subseteq V$ around the boundary nodes
of two adjacent blocks, which then yields a flow network (see Line <ref> and <ref>).
On each network, we run the FlowCutter algorithm <cit.> to derive a balanced minimum cut using incremental maximum flow computations.
FlowCutter and its parallelization are discussed in Sections <ref> and <ref>.
We then convert the minimum cut into a set of moves $M$ and an expected connectivity reduction $\expgain$.
If FlowCutter claims an improvement, i.e., if $\expgain \ge 0$, we apply the moves to the global partition and compute the exact reduction $\attrgain$,
based on which we either mark the blocks for further refinement, or revert the moves (see Line <ref> and <ref>).

Hypergraph $H = (V,E,c,\omega)$ and $k$-way partition $\Partition$ of $H$

Parallel Flow-Based Refinement

$\quotientgraph \gets \FuncSty {BuildQuotientGraph}(H, \Partition)$ *[r]see Section <ref>
(*[f]see Section <ref>) $\exists$ active $(V_i,V_j) \in \quotientgraph$
$B := \gets \FuncSty {ConstructRegion}(H,V_i,V_j)$ *[r]see Section <ref>
$(\flowhypergraph,\source,\sink) \gets \FuncSty {ConstructFlowNetwork}(H,B)$ *[r]see Section <ref>
$(M, \expgain) \gets \FuncSty {FlowCutterRefinement}(\flowhypergraph,\source,\sink)$ *[r]see Section <ref> – <ref>
(*[f]potential improvement)  $\expgain \ge 0$
$\attrgain \gets \FuncSty {ApplyMoves}(H,\Partition,M)$ *[r]see Section <ref>
(*[f]found improvement)  $\attrgain > 0$
mark $V_i$ and $V_j$ as active

(*[f]no improvement) $\attrgain < 0$  $\FuncSty{RevertMoves}(H, \Partition, M)$

Maximum Flows
A flow network $\flownetwork = (\flownodeset,\flowedgeset,\capacity)$
is a directed graph with a dedicated source node
$\source \in \flownodeset$ and sink node $\sink \in \flownodeset$ in which
each edge $e \in \flowedgeset$ has capacity $\capacity(e) \ge 0$. An $\sourcesink$-flow is a function
$f: \flownodeset \times \flownodeset \rightarrow \mathbb{R}$
that satisfies the capacity constraint $\forall u,v \in \flownodeset: f(u,v) \le \capacity(u,v)$,
the skew symmetry constraint $\forall u,v \in \flownodeset: f(u,v) = -f(v,u)$ and the
flow conservation constraint $\forall u \in \flownodeset \setminus \{\source,\sink\}: \sum_{v \in \flownodeset} f(u,v) = 0$.
The value of a flow $|f| := \sum_{v \in \flownodeset} f(\source,v) = \sum_{v \in \flownodeset} f(v,\sink)$ is defined as the
total amout of flow transferred from $\source$ to $\sink$.
An $\sourcesink$-flow $f$ is a maximum $\sourcesink$-flow if there exists no other $\sourcesink$-flow $f'$ with $|f| < |f'|$.
The residual capacity is defined as $r_f(e) = \capacity(e) - f(e)$.
An edge $e$ is saturated if $r_f(e) = 0$.
The residual network $\flownetwork_f = (\flownodeset, \flowedgeset_f, r_f)$ with
$\flowedgeset_f := \{(u,v) \in \flownodeset \times \flownodeset \mid r_f(u,v) > 0\}$ contains all non-saturated edges.
The max-flow min-cut theorem states that the value $|f|$ of a maximum $\sourcesink$-flow equals the weight
of a minimum cut that separates $\source$ and $\sink$ <cit.>.
This is also called a minimum $\sourcesink$-cut.
The minimum $\sourcesink$-cut can be derived by exploring the nodes reachable from the source or sink via
residual edges ($r_f(e) > 0$), which is also called the source-side or sink-side cut.

§.§ Parallel Active Block Scheduling

Sanders and Schulz <cit.> propose the active block scheduling strategy to apply
their flow-based refinement algorithm for bipartitions on $k$-way partitions.
Their algorithm proceeds in rounds.
In each round, it schedules all pairs of adjacent blocks where at least one is marked as active.
Initially, all blocks are marked as active.
If a search on two blocks finds an improvement, both are marked as active for the next round.

Our parallel implementation schedules multiple flow computations on adjacent block pairs in parallel.
We do not enforce any constraints on the block pairs processed concurrently, e.g., there can
be multiple threads running on the same block and they can also share some of their nodes.
We use $\min(t,\tau \cdot k)$ threads to process the active block pairs in parallel,
where $t$ is the number of available threads in the system and
the parameter $\tau$ controls the available parallelism in the scheduler.
With higher values of $\tau$, more block pairs are scheduled in parallel, which
can lead to more interferences between searches that operate on overlapping regions.
Threads that are not involved in scheduling can join parallel flow computations.
In a parameter study <cit.>, we found that $\tau = 1$ offers a good
trade-off between conflicting searches and scalability.

Initially, we push all pairs of adjacent blocks into a concurrent FIFO queue $\activequeue$.
The threads then poll from $\activequeue$ and if a search finds an improvement on a block pair $(V_i,V_j)$,
we mark both as active using a seperate bitset for each round. If either $V_i$ or $V_j$ becomes active, we
push all adjacent blocks into $\activequeue$ if they are not contained yet.
Thus, active block pairs of different rounds are stored interleaved in $\activequeue$ and the end of a
round does not induce a synchronization point as in the original algorithm <cit.>.
A round ends when all of its block pairs have been processed and all prior rounds have ended.
If the relative improvement at the end of a round is less than $0.1\%$, we immediately terminate the algorithm.

Apply Moves
Since concurrently scheduled flow computations can operate on overlapping regions, there are three
conflict types that can occur when applying a sequence of node moves $M$ to the global partition $\Partition$:

balance constraint violations, $\expgain \neq \attrgain$ (i.e., the expected does not
match the actual connectivity reduction), and nodes in $M$ may already be moved by other searches.

In practice, the running time to apply a sequence of node moves is negligible compared to solving flow problems <cit.>.
Thus, we can afford to use a lock so that only one thread applies moves at a time to address these conflicts.
First, we remove all nodes from $M$ that are not in their expected block.
Afterwards, we compute the block weights as if all remaining moves were applied.
If the resulting partition is balanced, we perform the moves, during which we aggregate the attributed gains $\attrgain$ of
each move. If $\attrgain < 0$, we revert all moves.

§.§ Flow Network Construction

A hypergraph $\flowhypergraph$ (left) induced by a region $B := B_1 \cup B_2$ and the flow network $\flownetwork$ (right)
given by the Lawler expansion of $\flowhypergraph$.

To improve the cut of a bipartition $\Partition = \{V_1, V_2\}$, we grow a size-constrained
region $B := B_1 \cup B_2$ with $B_1 \subseteq V_1$ and $B_2 \subseteq V_2$
around the cut hyperedges of $\Partition$ via two breadth-first-searches (BFS) <cit.>.
The first BFS is initialized with all boundary nodes of block $V_1$ and continues to add nodes to $B_1$ as long as
$c(B_1) \le (1 + \alpha \varepsilon) \lceil \frac{c(V)}{2} \rceil - c(V_2)$, where $\alpha$ is an input parameter.
The second BFS that constructs $B_2$ proceeds analogously.
We then contract all nodes in $V_1 \setminus B$ to the source $\source$ and
$V_2 \setminus B$ to the sink $\sink$ <cit.>
and obtain a coarser hypergraph $\flowhypergraph = (\flownodeset, \flowedgeset)$.
The flow network $\flownetwork$ is then given by the Lawler expansion of $\flowhypergraph$ <cit.>, which is illustrated
in Figure <ref>.
For each hyperedge $e \in \flowedgeset$, we add two nodes $\innode$ and $\outnode$ and a bridging edge
$(\innode,\outnode)$ with capacity $\capacity(\innode,\outnode) = \omega(e)$ to $\flownetwork$. For each pin $u \in e$,
we add two edges $(u,\innode)$ and $(\outnode,u)$ with infinite capacity to $\flownetwork$.
Note that we do not construct $\flownetwork$ explicitly in our actual implementation, since our maximum flow algorithm runs
on $\flowhypergraph$ by implicitly exploiting the structure of the Lawler expansion.

The parameter $\alpha$ controls the size of the flow network.
For $\alpha = 1$, each flow computation yields a balanced bipartition
with a possibly smaller cut in the original hypergraph, since only nodes of $B$ can move
to the opposite block ($c(B_1) + c(V_2) \le (1 + \varepsilon)\lceil \frac{c(V)}{2} \rceil$ and vice versa for block $B_2$).
Larger values for $\alpha$ lead to larger flow problems with potentially smaller minimum cuts,
but also increase the likelihood of violating the balance constraint.
However, this is not a problem since the flow-based refinement routine guarantees balance through incremental minimum cut computations (see Section <ref>).
In practice, we use $\alpha = 16$ (also used in KaHyPar <cit.>).
We additionally restrict the distance of each
node $v \in B$ to the cut hyperedges to be smaller than or equal to a parameter $\delta$ ($= 2$). We observed that it
is unlikely that a node far way from the cut is moved to the opposite block by the flow-based refinement.

§.§ The FlowCutter Algorithm

In this section, we discuss the flow-based refinement on a bipartition.
We introduce the aforementioned FlowCutter algorithm <cit.>,
which parallelization is described in the next section.

To speed up convergence and make parallelism worthwhile, we propose an optimization named bulk piercing.

Algorithm Overview

FlowCutter solves a sequence of incremental maximum flow problems until a balanced bipartition is found.
Algorithm <ref> shows the pseudocode for the approach.
In each iteration, first the previous flow (initially zero) is augmented to a maximum flow regarding the current source set $S$ and sink set $T$.
Subsequently, the node sets $S_r, T_r \subset \flownodeset$ of the source- and sink-side cuts are derived.
This is done via residual (parallel) BFS (forward from $S$ for $S_r$, backward from $T$ for $T_r$).
The node sets induce two bipartitions $(S_r, \flownodeset \setminus S_r)$ and $(\flownodeset \setminus T_r, T_r)$.
If neither is balanced, all nodes on the side with smaller weight are transformed to a source (if $c(S_r) \leq c(T_r)$) or a sink otherwise.
As this would yield the same cut in the next iteration, we add one additional node, called piercing node, to the terminal set of the smaller side.
Thus, the bipartitions contributed by the currently smaller side will be more balanced with a possibly larger cut in future iterations.
Since the smaller side is grown, this process will converge to a balanced bipartition.

The FlowCutter Algorithm
Original hypergraph $H = (V,E,c,\omega)$, flow network $\flowhypergraph = (\flownodeset, \flowedgeset)$
and a source $\source \in \flownodeset$ and sink $\sink \in \flownodeset$
Balanced Bipartition of $\flowhypergraph$
$S \gets \{ \source \} , T \gets \{ \sink \}$ *[r]initialize source and sink set
initialize flow $f: V \times V \rightarrow \mathbb{R}_{\ge 0}$ with $\forall (u,v) \in V \times V: f(u,v) = 0$
() no balanced bipartition found
$f \gets \FuncSty{ParallelMaxPreflow}(\flowhypergraph, S, T, f)$ *[r]augment $f$ to a maximum preflow
$(S_r,T_r) \gets$ derive source- and sink-side cut $S_r, T_r \subset \flownodeset$
() $(S_r, \flownodeset \setminus S_r)$ is balanced  $(S_r, \flownodeset \setminus S_r)$
($\flownodeset \setminus T_r, T_r)$ is balanced   $(\flownodeset \setminus T_r, T_r)$
() $c(S_r) \leq c(T_r)$
$S \gets S_r \cup \FuncSty{selectPiercingNode}(S \cup S_r)$
$T \gets T_r \cup \FuncSty{selectPiercingNode}(T \cup T_r)$

For our purpose, there are two important piercing node selection heuristics: avoid augmenting paths <cit.> and distance from cut <cit.>.
Whenever possible, a node that is not reachable from the source or sink should be picked, i.e., $v \in \flownodeset \setminus (S_r \cup T_r)$.
Such nodes do not increase the weight of the cut, while improving balance <cit.>.
As a secondary criterion, larger distances from the original cut are preferred, to reconstruct parts of it.

Bulk Piercing Optimization

On larger instances, piercing only one node per iteration converges slowly. We therefore increase the amount of work in each iteration
by piercing multiple nodes, as long as we are far from balance.

To achieve a small number of iterations (e.g., poly-log) we set a goal on the weights of the sides of the bipartition, and pierce more aggressively the further we are from it.
Assume, we want to pierce the source side next and have already performed $r-1$ piercing iterations on it.
In the $r$-th iteration we want to add $\frac{1}{2^r} \cdot (\frac{c(V)}{2} - c(S))$ new weight to the source side, where $c(S)$ is the weight of the initial source-side terminals (before any piercing) and $\frac{c(V)}{2}$ is the weight of a perfectly balanced bipartition.
Thus the overall weight goal for the $r$-th iteration on the source side is set to $(\frac{c(V)}{2} - c(S)) \sum_{i=1}^{r} \frac{1}{2^i}$.
This is chosen such that we allow a lot of progress early on and become more careful as we get closer to a balanced bipartition.
We track the average weight added per node in previous iterations and from this estimate the number of required piercing nodes to reach the goal for the $r$-th iteration.
To boost measurement accuracy, we pierce only one node for the first few rounds, and then switch to bulk piercing.

§.§ Parallel Maximum Flow Algorithm

Maximum flow algorithms are notoriously difficult to parallelize
efficiently <cit.>.
The synchronous push-relabel approach of Baumstark  <cit.> is a recent algorithm that
sticks closely to sequential FIFO and thus shows good results.
We first describe the sequential push-relabel algorithm proposed by Goldberg and Tarjan <cit.> and
then briefly outline its parallelization.

We conclude with implementation details and intricacies of using FlowCutter with preflows.

Push-Relabel Algorithm
The push-relabel algorithm <cit.> stores a distance label $\distance(u)$ and an excess value $\excess(u) := \sum_{v \in \flownodeset} f(v,u)$ for each node.
It maintains a preflow <cit.> which is a flow where the conservation constraint is replaced by $\excess(u) \geq 0$.
The distance labels represent a lower bound for the distance of each node to the sink.
A node $u \in \flownodeset$ is active if $\excess(u) > 0$. An edge $(u,v) \in \flowedgeset$ is admissible
if $r_f(u,v) > 0$ and $\distance(u) = \distance(v) + 1$.
A $\text{push}(u,v)$ operation sends $\delta = \min(\excess(u),r_f(u,v))$ flow units over $(u,v)$.
It is applicable if $u$ is active and $(u,v)$ is admissible.
A $\text{relabel}(u)$ operation updates the distance label of $u$ to $\min(\{\distance(v) + 1 \mid r_f(u,v) > 0\})$, which is applicable if $u$ is active and has no admissible edges.
The distance labels are initialized to $\forall u \in \flownodeset \setminus \{ \source \}: \distance(u) = 0$ and $\distance(\source) = |V|$ and all source edges are saturated.
Efficient variants use the discharge routine, which repeatedly scans the edges of an active node until its excess is zero.
All admissible edges are pushed and at the end of a scan, the node is relabeled.

The global relabeling heuristic <cit.> frequently assigns exact
distance labels by performing a reverse BFS from the sink to reduce relabel work in practice.
Note that a maximum preflow already induces a minimum sink-side cut, so if only a minimum cut is required, the algorithm can already stop
once no active nodes with distance label $< n$ exist.

The parallel push-relabel algorithm of Baumstark  <cit.> proceeds in rounds in which all active nodes are discharged in parallel.
The flow is updated globally, the nodes are relabeled locally and the excess differences are aggregated in a second array using atomic instructions.
After all nodes have been discharged, the distance labels $\distance$ are updated to the local labels $\distance'$ and the excess deltas are applied.
The discharging operations thus use the labels and excesses from the previous round.
This is repeated until there are no nodes with $\excess(v) > 0$ and $\distance(v) < n$ left.
To avoid concurrently pushing flow on residual arcs in both directions (race condition on flow values), a deterministic winning criterion on the old distance labels is used to determine which direction to push, if both nodes are active.
If an arc cannot be pushed due to this, the discharge terminates after the current scan, as the node may not be relabeled in this round.
The rounds are interleaved with global relabeling <cit.>, after linear push and relabel work, using parallel reverse BFS in the residual network.
We additionally fixed an undocumented bug in the original algorithm (not source code) for which we refer the reader to Ref. <cit.>.

Intricacies with Preflows and FlowCutter
A maximum preflow only yields a sink-side cut via the reverse residual BFS, but we also need the source-side cut.
We can run flow decomposition <cit.> to push excesses back to the source.
However, flow decomposition is difficult to parallelize <cit.>.
Instead, we initialize the forward residual BFS with all active non-sink excess nodes.
This finds the reverse paths that carry flow from the source to the excess nodes, which is what we need.

Furthermore, when transforming a node with positive excess to a sink, its excess must be added to the flow value.
This only happens when piercing, as sink-side nodes have no excess.

Finally, we want to reuse the distance labels from the previous round to avoid re-initialization overheads.
However, as the labels are a lower bound on the distance from the sink, piercing on the sink side invalidates the labels.

In this case, we run global relabeling to fix the labels and collect the existing excess nodes, before starting the main discharge loop.
When piercing on the source side the labels remain valid and new excesses are created.
These are added to the active nodes and we do not run an additional global relabeling.
The existing excess nodes are collected during regular global relabel runs.

Implementation Details

Since $(\innode, \outnode)$ is the only outgoing edge of $\innode$ with non-zero capacity in the Lawler expansion (see Figure <ref>),
the flow on edges $(u, \innode)$ is also bounded by $\omega(e)$ (instead of $\infty$). Adding these capacities is a trivial optimization, but
significantly accelerates the algorithm and increases the available parallelism. This can be explained by the fact that a hypernode $u$ does not
immediately relieve all of its excess to one of its incident nets $e \in \incnets(u)$ during the discharge routine, which is later pushed back
due to the $\omega(e)$ bound. We set the capacities $c(u,\innode)$ to $\infty$ again when deriving the source- and sink-side cut, since only
bridging edges can be cut in the Lawler expansion.

Moreover, we observed that the number of active nodes follows a power-law distribution. Due to little work in later rounds, it takes
many rounds to trigger the global relabeling step that also terminates the algorithm when a maximum preflow is found.
Therefore, we perform additional relabeling if the flow value has not changed for some rounds (500),
and only few active nodes ($< 1500$) were available in each.

§ A PARALLELIZATION OF THE $N$-LEVEL PARTITIONING SCHEME

Our multilevel algorithm contracts a clustering of highly-connected nodes on each level, which
induces a hierarchy with a logarithmic number of levels. In contrast, KaHyPar <cit.> – the currently
best sequential partitioning algorithm with regards to solution quality – contracts only a single node on each level.
Correspondingly, in each refinement step, only a single node is uncontracted followed
by a highly-localized search for improvements around the uncontracted node.
This technique produces almost $n$ levels and is therefore known as the $n$-level partitioning scheme.
More levels provide more opportunities to refine the current solution <cit.>

at different
granularities but also increases the running time of multilevel algorithms. Therefore, KaHyPar is
the method of choice for computing high-quality partitions but comes at the cost of substantially higher running times than
other systems – prohibitively so for very large hypergraphs. Although $n$-level partitioning seems inherently sequential,
we present the first shared-memory parallelization of the technique which achieves good speedups and comparable
solution quality to KaHyPar in a fraction of its running time.

We start this section with a formal definition of the (un)contraction operation, and provide
a high-level overview of our $n$-level partitioning algorithm.
We then discuss and present solutions for the main challenges in this algorithm:
finding a parallel schedule of (un)contraction operations and performing them on a dynamic hypergraph data structure in parallel.
We conclude the algorithm description with a discussion on how the refinement algorithms from the previous section are
integrated into the $n$-level algorithm.

The Contraction and Uncontraction Operation
Contracting a node $v$ onto another node $u$ replaces $v$ with $u$ in all nets $e \in \incnets(v) \setminus \incnets(u)$
and removes $v$ from all nets $e \in \incnets(u) \cap \incnets(v)$. The weight of node $u$ is then $c(u) + c(v)$. We call
$u$ the representative of the contraction and $v$ its contraction partner. Uncontracting a node $v$ reverses
the corresponding contraction operation.

Hypergraph $H = (V,E)$, number of blocks $k$
$k$-way partition $\Partition$ of $H$

The $n$-Level Partitioning Algorithm

() $V$ has too many nodes
() $u \in V$ in random order
$v \gets \operatorname{arg\,max}_{v} \sum_{e \in \incnets(u) \cap \incnets(v)} \frac{\omega{(e)}}{|e| - 1}$ *[r]find best contraction partner for $u$
() $(v,u)$ is eligible for contraction
contract $v$ onto $u$

$\Partition \gets$ InitialPartition($H, k$)
$\batches = \langle B_1, \dots, B_l \rangle \gets \FuncSty {ConstructBatches}(\contractionforest)$
(*[f]$|B| \approx \maxbatchsize$) $B \in \batches$
() $v \in B$
uncontract $v$ from $\rep[v]$$\Partition[v] \gets \Partition[\rep[v]]$

Algorithm Overview
Algorithm <ref> shows the high-level pseudocode of our $n$-level partitioning
algorithm. In the coarsening phase, we iterate in parallel over all nodes and find
the best contraction partner $v$ for each node $u$ using the heavy-edge rating function <cit.> (similar as in our
multilevel algorithm). We then check whether $v$ can be contracted right away
onto $u$ or if there are any other pending contractions that must be performed before.
In the latter case, we transfer the responsibility
of contracting $v$ onto $u$ to the thread resolving the last dependency that defers the contraction.
Once the hypergraph is small enough, we compute an initial partition into $k$ blocks.

Uncontracting only a single node followed by a localized refinement step is
inherently sequential, which is why we have to relax the $n$-level idea in the uncoarsening phase.
We construct a sequence of batches $\batches = \langle B_1, \dots, B_l \rangle$ of contracted nodes,
such that $|B_i| \approx \maxbatchsize$ where $\maxbatchsize$ is an input parameter.
Batches are processed one after another, enabling the uncontraction of nodes in subsequent batches.
Nodes in the same batch are uncontracted in parallel.
The main challenge is to identify which nodes can or even must appear in the same batch.
After uncontracting each batch, we apply highly-localized refinement algorithms around the batched nodes.

A Forest-Based Scheduling of Contraction Operations
Let us consider a sequence of contractions $\contractions := \splitatcommas{\langle (v_1,u_1),\ldots, (v_n,u_n) \rangle}$
executed exactly in this order ($v_i$ is contracted onto $u_i$).
In this sequence, each node is contracted onto at most one representative, and there are no cyclic contraction dependencies.

Thus, the sequence of contractions form a forest $\contractionforest := (V,\contractions)$ if interpreted
as a graph with directed edges $(v_i,u_i) \in \contractions$.  We will refer to $\contractionforest$ as the
contraction forest. Our parallelization uses the observation that there exists several permutations of $\contractions$
leading to the same contraction forest $\contractionforest$. We can contract a node as soon as all of its
children in $\contractionforest$ have been contracted. To obtain parallelism, different subtrees and siblings
can be contracted independently, i.e., we traverse $\contractionforest$ in a bottom-up fashion in parallel.

In our actual algorithm, we do not know $\contractionforest$ in advance.
However, we show how to construct $\contractionforest$ dynamically and from that we derive a
parallel schedule of contraction operations. We call a contraction $(v,u)$ compatible with existing
contractions if it satisfies the following three conditions: (i) $v$ must be a root of $\contractionforest$, (ii) adding $(v,u)$ to
$\contractionforest$ must not induce a cycle, and (iii)
the contraction of $u$ onto its parent in $\contractionforest$ must not have started yet.
We represent $\contractionforest$ using an array $\rep$ of size $n$ storing for each node its representative ($u$ is a root
if $\rep[u] = u$). Additionally, we use a zero-initialized array $\pending$, where $\pending[u]$ stores the number of
children of $u$ whose contraction is not finished. If $\pending[u] = 0$ and $\rep[u] \neq u$, we assume that the contraction of
$u$ onto $\rep[u]$ has started and prevent further contractions onto $u$. The entries $\rep[u]$ and $\pending[u]$ are only modified while
holding a node-specific lock for $u$.

If we add a contraction $(v,u)$ to $\contractionforest$, we first lock $v$ and check if $v$ is a root.
If $\rep[v] \neq v$, we discard the contraction as another thread has already selected a representative for $v$.
Otherwise, we walk the path towards the root of $u$'s tree in $\contractionforest$ to find the lowest ancestor $w$
of $u$ whose contraction has not started yet ($rep[w] = w$ or $\pending[w] > 0$, in most cases $w = u$). If $v$ is found
on this path, the contraction is discarded, as it would induce a cycle in $\contractionforest$. If no cycle is found,
we lock $w$, and check $\rep[w]$ and $\pending[w]$ again. If they changed, we find a new suitable ancestor and perform the
cycle check again. Otherwise, we set $\rep[u] = w$ and increment $\pending[w]$ by one, and unlock $v$ and $w$.
We immediately contract $v$ onto $w$ if $\pending[v] = 0$ and subsequently reduce $\pending[w]$ by one. If this
reduces $\pending[w]$ to zero, we recursively apply this process to the contraction $(w, \rep[w])$ if $\rep[w] \neq w$.

Batch Uncontractions
For the uncoarsening phase, we construct a sequence of batches $\batches = \langle B_1, \ldots, B_l \rangle$ where each
batch $B_i$ contains roughly $\maxbatchsize$ contracted nodes that can be uncontracted independently in parallel.
The batch size $\maxbatchsize$ is an input parameter (set to $\maxbatchsize = 1000$ in our implementation) that interpolates between
scalability (high values) and the inherently sequential $n$-level scheme ($\maxbatchsize = 1$). Uncontracting a batch $B_i$
resolves the last dependencies required to uncontract the next batch $B_{i+1}$. After each batch uncontraction, we apply a highly-localized version of the
label propagation and FM algorithm searching for improvements in a small region around the uncontracted nodes.

Uncontracting $w$ increases the cut by one since we do not uncontract $v$ and $w$ in reverse order.

We construct the batches via a top-down traversal of the contraction forest $\contractionforest$.
There are two constraints that we need to consider when constructing the batches:
(i) a node must appear in a batch strictly after the batch containing its representative, and (ii) siblings in
$\contractionforest$ must be uncontracted in reverse order of contraction. The second condition prevents uncontractions
increasing the cut size as illustrated in Figure <ref>, which would violate a fundamental property of
the multilevel scheme. Since contractions can be performed at the same time, it is often not possible to define a strict order in which we have to
revert the contractions. We therefore associate each contraction $(v,u)$ with a time interval $[s_v, e_v]$ by atomically incrementing
a counter before starting ($s_v$) and after finishing ($e_v$) a contraction operation. If the time interval of two nodes overlap,
we assume they were contracted at the same time, otherwise one is strictly earlier than the other. Among siblings, we compute
the transitive closure of nodes with overlapping time intervals and order them decreasingly if one is stricly earlier than the other.
We then use this as the reverse order of contractions, while we add siblings with overlapping time intervals to the same batch.
As this is only a high-level description of the batch construction algorithm, we refer the reader to Ref. <cit.> for more details.

The Dynamic Hypergraph Data Structure
Figure <ref> illustrates the dynamic hypergraph data structure that stores the pin-list of each net $e$ and the incident
nets $\incnets(u)$ of each node $u$ using two seperate adjacency arrays. When we contract a node $v$ onto another node $u$, we iterate over the incident nets
of $v$ and search for $u$ and $v$ in the pin-list of each net $e \in \incnets(v)$. If we do not find $u$ in $e$, we replace $v$ with
$u$ ($e \in \incnets(v) \setminus \incnets(u)$). Otherwise, we swap $v$ to the end of $e$'s pin-list and decrement
the size of $e$ by one ($e \in \incnets(u) \cap \incnets(v)$), dividing its pin-list into an active and inactive part.
We use a seperate lock for each net to synchronize edits to the pin-lists. We further mark the nets $e \in \incnets(u) \cap \incnets(v)$
in a bitset $X$, which we then use to update the incident nets of $u$.

The key idea for updating the incident nets is to remove $\incnets(u) \cap \incnets(v)$ from $\incnets(v)$ and concatenate $u$
and $v$ in a doubly-linked list $L_u$. All nodes contracted onto $u$ are then stored in $L_u$. We can then iterate over the incident
nets of $u$ by iterating over all entries $w \in L_u$ and the modified $\incnets(w)$ arrays.
We associate each $\incnets(w)$ array with a counter $t_w$ and each entry $e \in \incnets(w)$ with a marker $t_{w,e}$ (initially set to zero).
Entries with markers $\ge t_w$ are active, i.e., were not remove yet.
For removing $\incnets(u) \cap \incnets(v)$ from $\incnets(v)$
(marked in bitset $X$), we iterate over all nodes $w \in L_v$ and increment $t_w$ by one. We then iterate over the previously
active entries of $\incnets(w)$ (now marked with $t_w - 1$) and if an entry is not in $X$, we set its marker to $t_w$.
Otherwise, we swap the entry to the end of the active part but keeping its marker at $t_w - 1$.
A simpler approach would be to represent the incident nets of each node as an adjacency list and add $\incnets(v) \setminus \incnets(u)$ to $\incnets(u)$,
as it is done in KaHyPar <cit.>. However, this could lead to quadratic memory usage, and is therefore
not practical for large hypergraphs.

When uncontracting a node $v$ from its representative $u$, we first restore $L_v$ from $L_u$. To do this, we additionally store the last node in
$L_v$ at the time $v$ is contracted onto $u$. We then iterate over all nodes $w \in L_v$ and decrement their counters $t_w$ by one. This
reactivates all nets $e \in \incnets(w)$ that became inactive due to contracting $v$, i.e., were part of $\incnets(u) \cap \incnets(v)$
before the contraction. In these nets, we swap $v$ to the active part of their pin-lists again. All previously active nets $e \in \incnets(w)$ (now marked
with $t_{w,e} > t_w$) were part of $\incnets(v) \setminus \incnets(u)$ before the contraction in which we then replace $u$ with $v$ again. Note that we sort all
pins in the inactive part of a pin-list by the batches in which they are uncontracted. Then, all pins of a net $e$ part of the current batch can be restored
simultaneously by appropriately incrementing the size of $e$. Only one thread that triggers the restore case on a net performs the restore operation,
which we ensure with an atomic test-and-set instruction.

Contraction operation applied on the dynamic hypergraph data structure.

Removing Single-Pin and Identical Nets
We remove single-pin nets and aggregate the weight of all identical nets at one representative after a pass over all nodes
in the coarsening phase using the same algorithm as already described in Section <ref>. This adds several synchronization points ($\approx \log{n}$)
at which we have to restore them in the uncoarsening phase. KaHyPar <cit.> removes these nets directly
after each contraction operation. However, doing this in the parallel setting would introduce additional dependencies for batches, which is why we
decided against it.

After uncontracting a batch, we run a highly-localized version of label propagation and FM refinement initialized
with the boundary nodes of the current batch. The searches then expand to a small region around the uncontracted nodes.
We complement the localized refinement with a refinement pass on the entire hypergraph after restoring single-pin and identical nets.
Here, we run FM (initialized with all boundary nodes) and flow-based refinement. Additionally, we implemented a concurrent
gain table update procedure for batch uncontractions for which we refer the reader to Ref. <cit.> for more details.

§ UNIFYING HYPERGRAPH AND GRAPH PARTITIONING

Hypergraph partitioning (HGP) is considered inherently more complicated <cit.> and therefore more complex in terms
of implementation and running time <cit.> than graph partitioning (GP). However, the high-level description of
partitioning algorithms often does not reveal any difference between the two. For example, label propagation refinement iterates over all nodes, and
moves each node to the block with the highest gain value. While the algorithm is widely used in GP and HGP, the main difference
in its implementation lies in the representation of the graph data structure and the computation of gain values. GP tools build on
data structures using one adjacency array to represent the neighbors of nodes, while HGP requires two adjacency arrays storing
the pin-lists of hyperedges and the incident nets of nodes. This results in a better cache utilization and faster access times for graph algorithms.
Moreover, the gain value of a node move for the edge cut metric depends on the block assignments of neighbors for GP,
while HGP tools have to maintain or compute the pin count values of hyperedges to decide whether or not it can be removed from the cut.
In the following, we present an optimized graph and partition data structure for GP implementing the interface of our hypergraph data
structure such that we can use them as a drop-in replacement in our partitioning algorithm. We focus on the multilevel algorithm
and refer the reader to Ref. <cit.> for a description of the $n$-level graph data structure.

An undirected and weighted graph $G = (V,E,c,\omega)$ can be considered as a hypergraph where each net contains only two pins (also called an edge).
Therefore, the definitions and notations for hypergraphs also apply to undirected graphs.
We define the weight of an edge $e = \{u,v\} \in E$ as $\omega(u,v) := \omega(e)$.
If $\{u,v\} \notin E$, then $\omega(u,v) = 0$.
An edge $\{u,u\} \in E$ is called a selfloop.
For a subset $V' \subseteq V$, $\omega(u, V') := \sum_{v \in V'} \omega(u,v)$
is the weight of all edges connecting node $u$ to $V'$.

§.§ Graph Data Structure

We use one adjacency array to represent an undirected graph.
The adjacency lists stores the directed edges $(u,v) \in \incnets(u)$ for each node $u$.
Note that we could reduce the memory overhead of the data structure by storing only the neighbors
$v \in \neighbors(u)$.
However, our partitioning algorithms are designed to operate on hypergraphs, which often request
the pin-list of a net using its ID (position in the adjacency array). This is only possible when storing
additional information about the source and target node for each edge.
Thus, our data structure requires twice as much memory as traditional graph data structures.

Our multilevel partitioning algorithm contracts a clustering of the nodes on each level.
The coarsening algorithm stores the clustering in an array $\rep$ where $\rep[u] = v$ stores the representative of $u$'s cluster.
For each representative $v$, we maintain the invariant that $\rep[v] = v$.

The contraction algorithm first remaps cluster IDs to a consecutive range by computing a parallel prefix sum on an array of size $n$ that
has a one at position $v$ if $v$ is a representative of a cluster and zero otherwise.
Then, we accumulate the weights and degrees of nodes in each cluster using atomic fetch-and-add instructions.
Afterwards, we copy the incident edges of each cluster to a consecutive range in a temporary adjacency array by computing
a parallel prefix sum over the cluster degrees.
We then iterate over the adjacency lists of each cluster in parallel, sort them, and remove selfloops and identical edges except for one representative
at which we aggregate their weights. Finally, we construct the adjacency array of the coarse graph by computing a parallel prefix
sum over the remaining cluster degrees.

§.§ The Partition Data Structure

For hypergraphs, our partition data structure stores the block assigments $\Partition$, the block weights $c(V_i)$,
the pin count values $\pinsinpart(e,V_i)$, and connectivity sets $\conset(e)$ for each net $e \in E$ and block
$V_i \in \Partition$. Since a graph edge connects only two nodes, we can remove the pin count values and
connectivity sets as we can calculate them on-the-fly. However, we used the synchronized writes to the pin count
values to update the gain table and compute the attributed gain values. We therefore present alternative approaches
for both techniques that exploit the properties of graphs.

The Gain Table

The connectivity metric reverts to the edge cut metric for plain graphs (since $\con(e) \le |e| = 2$). The gain value of moving
a node $u$ to another block $V_t$ is then defined as $\gain{u}{t} := \omega(u, V_t) - \omega(u, \nodeblock{u})$ (external minus internal edges).

Thus, the gain table for graph partitioning stores and maintains the $\omega(u,V_i)$ values for each node $u \in V$ and block $V_i \in \Partition$ ($n \cdot k$ entries).

If we move a node $u$ from block $V_s$ to $V_t$, we update the gain table by adding $\omega(u,v)$ to $\omega(v,V_t)$ and $-\omega(u,v)$
to $\omega(v,V_s)$ for each neighbor $v \in \neighbors(u)$ using atomic fetch-and-add instructions.
Hence, the complexity of the gain table updates when each node is moved at most once is $\sum_{u \in V} 2d(u) = \Oh{m}$.

Attributed Gains
For a node $u$ moved from block $V_s$ to $V_t$, we attribute a connectivity reduction or increase by
$\omega(e)$ to each net $e \in \incnets(u)$ based on the synchronized writes to $\pinsinpart(e,V_s)$
and $\pinsinpart(e,V_t)$. Since we do not longer maintain the pin count values,
we need another synchronization mechanism to decide if a node move removes an edge from the cut or makes
it a cut edge, and based on that attribute a reduction or an increase by the weight of the edge to the move.

We therefore use an array $B$ of size $m$ (initialized with $\perp$) to synchronize the node moves for each edge.
To compute the attributed gain of a node move $u$ from block $V_s$ to $V_t$, we iterate over all incident edges $e = \{u,v\} \in \incnets(u)$ and write
the target block $V_t$ of $u$ to $B[e]$ using an atomic compare-and-swap operation. If the operation
succeeds, no other thread has moved its neighbor $v$ yet. In this case, $e$ becomes an internal edge if $V_t = \nodeblock{v}$
(reduces the edge cut by $\omega(e)$) and a cut edge otherwise (increases the edge cut by $\omega(e)$).
If we do not succeed in setting $B[e]$ from $\perp$ to $V_t$, another thread has already moved or is currently
moving $v$ to another block. In both cases, its target block is $B[e]$ and we can compute the attributed gain value
for edge $e$ as before by comparing $V_t$ and $B[e]$. After calculating the attributed gain value for
each net $e \in \incnets(u)$, we set the block ID of $u$ to $V_t$. Note that the algorithm only works when
each node is moved at most once, as it is done in our refinement algorithms ($B[e]$ values are reset to $\perp$ after
each refinement round).

§ DETERMINISTIC PARTITIONING

A program is externally deterministic <cit.> if, given the same input, it produces the same output, each time it is run.

Sequential programs are usually deterministic by default, whereas parallel programs are non-deterministic by default due to randomness in scheduling.
Yet, researchers have advocated the benefits of deterministic parallel programs for several decades <cit.>.
It is easier to debug the program, to reason about performance and it yields reproducible results: in experiments and applications.
Unfortunately, with the exception of BiPart <cit.>, all published parallel partitioning algorithms so far are non-deterministic.
This stems from concurrently performed moves affecting other ongoing move decisions.

In this section, we present deterministic versions for a subset of the components in our multilevel framework: label propagation refinement, heavy-edge clustering for coarsening, and the Louvain community detection method <cit.> (optimizing the popular modularity metric), which we used to guide coarsening decisions.
These clustering algorithms all follow the local moving scheme.
Nodes are visited asynchronously in parallel and are moved to the best cluster in their neighborhood.

To achieve determinism, we use the synchronous local moving approach which is popular in distributed Louvain implementations for community detection <cit.>.
Moves are calculated but not applied until the end of a local moving round and thus do not influence one another.
The difficult part and difference to prior work is that not all calculated moves can be applied, for example due to the balance constraint. We must select a subset that is as profitable as possible.
We also break down each round into further sub-rounds, to trade off more frequent synchronization for more accurate gains.

Except for the use of non-internally deterministic sub-routines such as sorting, group-by, and emitting elements to a collection in parallel, our algorithms are internally deterministic, i.e., additionally pass through the same internal states on each run <cit.>.

Deterministic Label Propagation Refinement

In synchronous label propagation, we first calculate the highest gain move for each node in the current sub-round.
In a second step, we perform balance-preserving swaps between block pairs, prioritized by the gains of the calculated moves.
This generalizes a previous approach in SocialHash <cit.> to weighted hypergraphs, and thus allows the use in a multilevel framework.

For each block pair $(V_s, V_t)$, we sort the two move sequences $M_{st}$ from $V_s$ to $V_t$ and $M_{ts}$ from $V_t$ to $V_s$ by gain and then select a prefix
$M_{st}[0:i]$ and $M_{ts}[0:j]$, from each sequence to apply.
We use the node ID as tie-breaker for determinism.
Let $x(i,j) \coloneq \sum_{a = 0}^{i-1} c(M_{st}[a]) - \sum_{a = 0}^{j-1} c(M_{ts}[a])$ be the weight added to block $V_t$ and removed from block $V_s$ after swapping the nodes in the corresponding prefixes.
We call $i,j$ feasible if $- (\balancedconstraint{\max} - c(V_s)) \leq x(i,j) \leq \balancedconstraint{\max} - c(V_t)$, i.e., after the swaps the partition is still balanced.
To maximize gain, we look for the longest feasible prefixes.
This can be computed similar to merging two sorted arrays.
Keep two pointers $i,j$ to the current prefixes of $M_{st}[0:i], M_{ts}[0:j]$.
In each step advance the pointer of the sequence whose source block receives more weight, i.e., advance $j$ if $x(i,j) < 0$ or $i$ if $x(i,j) > 0$.
If $x(i,j) = 0$ advance either, if the end of the corresponding sequence is not yet reached.

The parallelization follows a common idea for parallel merging.
We first compute the cumulative gains of the sequences via parallel prefix sum operations.
Then the following algorithm is applied recursively to perform the selection.
We do binary search to find  the smallest index $q$ in the shorter sequence whose cumulative weight is not less than the cumulative weight of the middle of the longer sequence.
The two sub-sequence pairs to the left and right of the middle and $q$ can be searched independently in parallel.
Let $l$ denote the length of the longer sequence, then the algorithm does $\mathcal{O}(l)$ work and has $\mathcal{O}(\log^2(l))$ depth.
There are parallel merge algorithms with $\mathcal{O}(\log(l))$ depth, but these are more complicated and unlikely to yield faster running time in practice.

Additionally, we propose two optimizations that are helpful in practice but do not affect the theoretical running time.
If the right parts contain feasible prefixes we return them as we prefer longer prefixes, otherwise we return the result from the left parts.
If the prefixes at the splitting points are feasible, we can omit the left call.
Further, we can omit the right call if the cumulative weight at the middle of the longer sequence exceeds that at the end of the shorter sequence.

Deterministic Louvain Method

There is no weight constraint on clusters in the Louvain algorithm.
Therefore, we can apply all of the calculated moves.
However, there is an intricacy with floating point weights, which we need due to the edge weight model employed <cit.>.

In modularity optimization with the Louvain method, the cluster volume (weighted degree sum of the cluster) is part of the move decisions.
Usually, the volumes are updated in parallel after a node move <cit.>.
Unfortunaly, floating point arithmetic is not associative: different schedules will lead to slightly different rounded values, which actually resulted in non-deterministic outcomes.
One option is to recompute the volumes after each local moving subround.
However, this is substantially slower than only considering updates from moved nodes, particularly at later stages when fewer nodes are moved.
Therefore, we have to establish an order in which the volume updates of each cluster are aggregated.
First we group the updates by cluster, then sort by node ID. Subsequently, we perform a reduction on each group with static load balancing, which is needed for determinism.

Deterministic Clustering for Coarsening

During coarsening, we bound the weight of the heaviest cluster by an upper weight limit $c_{\max}$,
to ensure that initial partitioning can find a feasible solution.
The difference to refinement is that we have significantly more clusters, and only unclustered nodes (singletons) can move.
Therefore, the approach from refinement is not applicable here, which is why we use a simpler scheme.

Each unclustered node in a sub-round first determines its desired target cluster according to the heavy-edge rating function.
Then, we group the moves by the target cluster, and sort them in order of ascending node weight and use the node ID as tie-breaker for determinism.
For each group, we compute a prefix sum on the node weights, and apply the longest prefix that does not exceed the weight constraint, rejecting the remaining moves.

As an optimization to reduce the amount of work in the group-by stage (second largest bottleneck), we already sum up cluster weights during the target-cluster calculation step (main bottleneck).
If all moves into a cluster combined do not exceed the weight constraint, we simply approve them all and exclude the target cluster from the group-by stage.

For further details such as implementation details and group-by mechanisms used, as well as initial partitioning and contraction algorithms, we refer to <cit.>.

§ EXPERIMENTS

All presented algorithms have been made available in the Multi-Threaded Karlsruhe Hypergraph Partitioning framework
Mt-KaHyPar[Mt-KaHyPar is publicly available from <https://github.com/kahypar/mt-kahypar>].
It implements a parallel multilevel and $n$-level partitioning algorithm, as well as a deterministic version of the multilevel
algorithm and optimized data structures for graph partitioning. Mt-KaHyPar optimizes the connectivity metric for
hypergraph partitioning and the edge cut metric for graph partitioning.

The following experimental evaluation is structured as follows: We first decribe the four different benchmark sets composed of
over 800 graphs and hypergraphs, and discuss the experimental setup and methodology used in our experiments. We
then evaluate the time-quality trade-offs and speedups of Mt-KaHyPar's different partitioning configurations, and
analyze the running time of its algorithmic components in Section <ref>.
We conclude the evaluation by comparing Mt-KaHyPar to $25$ different sequential and parallel graph and hypergraph
partitioning algorithms in Section <ref>.

Instances We assembled four different benchmark sets. Two of the sets consist of graphs (G), while the other two consists of hypergraphs (HG).
The sets are further subdivided into medium-sized (M) and large instances (L). We abbreviate the name of a benchmark set, e.g., with .
The baseline denotes the size of the instances, while the subscript indicates whether it contains graphs or hypergraphs. We summarize the properties
of the instances contained in the benchmark sets in Figure <ref>[
We made all benchmark sets and detailed statistics of their properties publicly available
from <https://algo2.iti.kit.edu/heuer/talg/>.].
Note that all graphs and hypergraphs have unit node and (hyper)edge weights.

The hypergraph instances are derived from four sources encompassing three application domains:
the ISPD98 VLSI Circuit Benchmark Suite <cit.> (),
the DAC 2012 Routability-Driven Placement Contest <cit.> (),
the SuiteSparse Matrix Collection <cit.> (), and the International SAT Competition 2014 <cit.> ().
We interprete the rows and columns of a sparse matrix as nets and nodes, and a non-zero entry in a cell $(i,j)$ indicates whether
or not node $j$ is a pin of net $i$ <cit.>. We translate satisfiability formulas into three different hypergraph representations <cit.>.
The  resp.  representation interpretes the variables resp. literals as nodes, while the clauses form the hyperedges spanning the corresponding nodes.
The  representation models the clauses as nodes and variables as hyperedges.

Set  contains all 488 hypergraphs from the well-established benchmark set of Heuer and Schlag <cit.>
(18 , 10 , 184 , 276 , , and  instances). The benchmark  is composed of 94 large hypergraphs that
we selected in order to have more inputs where parallelization is important and useful. It contains the 8 largest  instances from  that we enhanced with 6 even larger
satisfiability formulas from the International SAT Competition 2014 <cit.> ($3 \cdot 14 = 42$ different hypergraph representations). We also included 42 sparse matrices with at least 15 million non-zeros, randomly sampled from
the SuiteSparse Matrix Collection <cit.>. Additionally,
we added all  instances from . The largest hypergraph of  has roughly two billion pins.

Our graph benchmark sets are composed of instances from the 10th DIMACS Implementation Challenge <cit.> (),
the Stanford Large Network Dataset Collection <cit.> and the Laboratory for Web Algorithms <cit.> (),
the DAC 2012 Routability-Driven Placement Contest <cit.> (),
the SuiteSparse Matrix Collection <cit.> (), and several randomly generated graphs <cit.> ().

Set  was initially assembled by Gottesbüren  <cit.> (195 graphs) from which we excluded the 39 largest graphs
and additionally added 16 social networks from the Stanford Large Network Dataset Collection <cit.> (114 , 30 , 15 , 3 , and
10  instances). The benchmark  contains 38 out of 42 instances from a graph collection assembled by Ahkremtsev <cit.> (four instances
were considered as too large as they were used to evaluate external memory algorithms). Additionally, we enhanced  with 15 graphs that we excluded from  and were
not contained in  yet (16 , 16 , 15 , and 6  instances). The largest graph of  has roughly two billion edges.

Summary of different properties for our benchmark sets. It shows for each
(hyper)graph (points), the number of nodes $|V|$, nets $|E|$ and pins $|P|$,
as well as the median and maximum net size ($\medsize$ and $\maxsize{e}$) and node degree ($\meddeg$ and $\maxsize{v}$).

Experimental Setup
Experiments on medium-sized instances ( and ) run on a cluster of Intel Xeon Gold 6230 processors (2 sockets with 20 cores each)
running at $2.1$ GHz with 96GB RAM.
In these experiments, we partition each (hyper)graph ten times using different random seeds into
$k \in \{2,4,8,16,32,64,128\}$ blocks with an allowed imbalance of $\varepsilon = 3\%$ and a time limit of eight hours.
Experiments on large instances ( and ) are done on an AMD EPYC 7702 processor (1 socket with 64 cores) running
at $2.0$–$3.5$ GHz with 1024GB RAM. Here, we partition each (hyper)graph three times using different random seeds into $k \in \{2,8,16,64\}$
blocks with an allowed imbalance of $\varepsilon = 3\%$ and a time limit of two hours. Note that we restrict the parameter space for experiments
on large instances due to limited computational resources. For graph partitioning, we configure the algorithms to optimize the
edge cut metric, while we focus on the connectivity metric for hypergraphs. We will also refer to both metrics as the solution
quality of a partition.

Aggregating Performance Numbers
We call a (hyper)graph partitioned into $k$ blocks an instance. For each instance, we aggregate running times
and the solution quality using the arithmetic mean over all seeds. To further aggregate over multiple instances, we use the geometric mean
for absolute running times and self-relative speedups. If all runs of an algorithm produced an imbalanced partition or ran into the time limit on
an instance, we consider the solution as infeasible. In plots, we mark imbalanced solutions with 55 and similarly instances
that timed out with . Runs with imbalanced partitions are not excluded from aggregated running times. For runs that exceeded the time limit,
we use the time limit itself in the aggregates. When comparing running times,
we say that an algorithm $\mathcal{A}$ is faster than $\mathcal{B}$ by a factor of $x$ on average if
$x := \nicefrac{\gmean{t_\mathcal{B}}}{\gmean{t_\mathcal{A}}} > 1$ where $\gmean{t_\mathcal{A}}$ and
$\gmean{t_\mathcal{B}}$ are geometric mean running times of $\mathcal{A}$ and $\mathcal{B}$.

Performance Profiles
Performance profiles can be used to compare the solution quality of different algorithms <cit.>.
Let $\mathcal{X}$ be the set of all algorithms, $\mathcal{I}$ the set of instances, and $q_{\mathcal{A}}(I)$ the quality of algorithm
$\mathcal{A} \in \mathcal{X}$ on instance $I \in \mathcal{I}$ ($q_{\mathcal{A}}(I)$ is the arithmetic mean over all seeds).
For each algorithm $\mathcal{A}$, performance profiles show the fraction of instances ($y$-axis) for which $q_\mathcal{A}(I) \leq \tau \cdot \text{Best}(I)$, where $\tau$ is on the $x$-axis
and $\text{Best}(I) := \min_{\mathcal{A}' \in \mathcal{X}}q_{\mathcal{A}'}(I)$ is the best solution produced by an algorithm $\mathcal{A}' \in \mathcal{X}$ for an instance $I \in \mathcal{I}$.
For $\tau = 1$, the $y$-value indicates the percentage of instances for which an algorithm $\mathcal{A} \in \mathcal{X}$ performs best.
Achieving higher fractions at smaller $\tau$ values is considered better.
The 55- and -tick indicates the fraction of instances for which all runs of that algorithm produced an imbalanced solution or timed out.
Note that these plots relate the quality of an algorithm to the best solution and thus do not permit a full ranking of
three or more algorithms.

Effectiveness Tests
Ahkremtsev  <cit.> introduce effectiveness tests to compare solution quality when two algorithms are given a
similar running time by performing additional repetitions with the faster algorithm.
Following this approach, we generate virtual instances that we compare using performance profiles.
Consider two algorithms $\mathcal{A}$ and $\mathcal{B}$, and an instance $I$.
We first sample one run of both algorithms for instance $I$.
Let $t_\mathcal{A}^1, t_\mathcal{B}^1$ be their running times and assume that $t_\mathcal{A}^1 \le t_\mathcal{B}^1$.
We then sample additional runs without replacement for $\mathcal{A}$ until their accumulated time exceeds $t_\mathcal{B}^1$ or all runs have been sampled.
Let $t_\mathcal{A}^2, \dots, t_\mathcal{A}^l$ denote their running times.
We accept the last run with probability $(t_\mathcal{B}^1 - \sum_{i = 1}^{l-1} t_\mathcal{A}^i) / t_\mathcal{A}^l$ so that the expected time for the sampled runs of $\mathcal{A}$ equals $t_\mathcal{B}^1$.
The solution quality is the minimum out of the sampled runs.
For each instance, we generate $10$ virtual instances.

Statiscal Significance Tests
We use the Wilcoxon signed-rank test <cit.> to determine whether the difference of the solutions produced by two algorithms
is statiscally significant. At a $1\%$ significance level ($p \le 0.01$), a Z-score with $|Z| > 2.576$ is deemed significant <cit.>.

§.§ Evaluation of Framework Configurations

In this section, we present a detailed evaluation of our shared-memory partitioning algorithm Mt-KaHyPar. We first describe its different
configurations and compare them regarding solution quality and running time. We then discuss the running times of the different
algorithmic components, present speedups, and evaluate the impact of our optimizations for plain graphs.

Framework Configurations

The Mt-KaHyPar framework provides a multilevel (Mt-KaHyPar-D, Default) and $n$-level partitioning algorithm
(Mt-KaHyPar-Q, Quality), as well as configurations extending them with flow-based refinement
(Mt-KaHyPar-D-F and Mt-KaHyPar-Q-F, Flows). It also implements a determistic version of the multilevel
algorithm (Mt-KaHyPar-SDet, Speed-Deterministic), which does not use the FM algorithm.
The code is written in $17$, parallelized using the  parallelization library <cit.>, and compiled using 9.2 with the
flags .
All of these algorithms have a large number of configuration options and were carefully tuned to provide the best trade-off between solution
quality and running time. However, a detailed parameter tuning study is beyond the scope of this paper.
We already mentioned specific choices for relevant parameters in the text and refer the reader to our conference publications <cit.> and
the dissertation of Heuer <cit.> for a detailed overview[Parameter tuning was done on a subset  of  that consists
of $100$ instances not contained in . We compared the quality produced by different partitioning algorithms on  and $\setminus$ using performance profiles and found
that they do not differ <cit.>. Thus, we decided to include the parameter tuning instances in the final evaluation to increase the evidence of the following
experimental results.].

Time-Quality Trade-Off

Performance profiles and running times comparing the different configurations
of Mt-KaHyPar executed with $10$ threads on .

Figure <ref> compares the solution quality of the partitions produced by the different configurations
of Mt-KaHyPar and their running times relative to Mt-KaHyPar-D on . The
configurations can be ranked from lowest to highest quality as follows: Mt-KaHyPar-SDet ( $1.25$s),
Mt-KaHyPar-D ($0.88$s), Mt-KaHyPar-Q ($2.99$s), Mt-KaHyPar-D-F ($2.73$s), and Mt-KaHyPar-Q-F ($5.08$s).
The ranking looks similar for running times except for Mt-KaHyPar-D which is faster than Mt-KaHyPar-SDet.
However, this changes when we compare their running times on the larger instances of .
Here, our deterministic configuration is faster than Mt-KaHyPar-D (Mt-KaHyPar-SDet: $3.14$s vs
Mt-KaHyPar-D: $4.65$s with 64 threads). For smaller instances, initial partitioning is the most time-consuming component
since we stop coarsening when we reach $160k$ nodes which can be close to the original number of nodes for some instances (e.g., $10240$ nodes for $k = 64$).
To reduce the running time of initial partitioning, Mt-KaHyPar-D adaptively adjusts the number
of repetitions of the different algorithms in the bipartitioning portfolio based on their success so far.
For larger instances, the smallest hypergraph is often significantly smaller than the input, and therefore
the running time of initial partitioning becomes negligible compared to the other phases.

The median improvement in solution quality of Mt-KaHyPar-D over Mt-KaHyPar-SDet is $6\%$, while flow-based refinement (Mt-KaHyPar-D-F)
improves Mt-KaHyPar-D by $4.2\%$ in the median at the cost of a $3$ times slower running time on average.
When we compare the multilevel (Mt-KaHyPar-D) and $n$-level partitioning algorithm (Mt-KaHyPar-Q), we see that $n$-level partitioning produces
partitions that are $1.9\%$ better than those produced by our multilevel algorithm in the median, but its running time is $3.4$ times slower on average.
The differences in solution quality and running time are less pronounced when both configurations use flow-based refinement (median improvement of Mt-KaHyPar-Q-F
over Mt-KaHyPar-D-F is $0.6\%$). Note that multilevel partitioning with flow-based refinement produces better partitions than our $n$-level
configuration ($2\%$), while it is also slightly faster.

We have seen that using stronger refinement algorithms leads to substantially better solution quality at the cost of higher running times.
Moreover, traditional multilevel algorithms can produce better partitions than $n$-level algorithms when flow-based refinement is used.

Effectiveness Tests

Our $n$-level algorithm computes better partitions than our multilevel algorithm without flow-based refinement, but is $3$ times slower on average.
When both configurations use flow-based refinement, the difference in solution quality becomes less pronounced.
We therefore use effectiveness tests to compare Mt-KaHyPar-D(-F) and Mt-KaHyPar-Q(-F) when both are given the same amount of time by
performing additional repetitions with the faster algorithm until the accumulated running time equals the running time of the slower algorithm.

Figure <ref> shows the results of these experiments. As we can see, the performance lines of
Mt-KaHyPar-D(-F) and Mt-KaHyPar-Q(-F) are almost identical in the performance profiles. This means
that Mt-KaHyPar-D(-F) computes partitions of comparable quality to its $n$-level counterpart when we give more time for
additional repetitions.

In contrast to the prevalent perception in the literature that more levels lead to better partitioning results <cit.>,