parRSB: Exascale Spectral Element Mesh Partitioning

parRSB: Exascale Spectral Element Mesh Partitioning Thilina Ratnayaka

Paul Fischer

June 15, 2026

arXiv:2606.14659v1 [cs.DC] 12 Jun 2026

Abstract We introduce parRSB – a parallel, highly scalable graph partitioner for spectral element meshes that produce high quality partitions. parRSB is based on Recursive Spectral Bisection (RSB) algorithm implemented on the dual graph of the input mesh. RSB uses the Fiedler vector, which is the eigenvector associated with the smallest non-zero eigenvalue of the Laplacian matrix of the dual graph for making partitioning decisions and tries to minimize the communication volume between the partitions. We implemented two numerical methods: Lanczos, and Inverse iteration using Conjugate Gradient method to compute the Fiedler vector. We present partitioning results using parRSB on Summit and Frontier supercomputers at Oak Ridge National Laboratory to illustrate the quality of the partitions produced by parRSB and the scalability of our implementation. We also present results for some of the optimizations we did to speed up the partitioning process. 1

Introduction

High Performance Computing (HPC) landscape is changing rapidly from CPU based systems to accelerator based systems powered by General Purpose Graphics Processing Units (GPGPUs or simply GPUs). This is evident by the fact that 7 of the top 10 computers of the latest TOP500 list of supercomputers [1] are accelerator based systems. The fastest, Frontier supercomputer at Oak Ridge Leadership Computing Facility (OLCF) using AMD Instinct™ MI250 GPUs is the first and the only exascale system on the list so far. A second exascale system, Aurora at Argonne Leadership Computing Facility (ALCF) is expected to be deployed by the end of 2023 and will also be an accelerator based system with Intel® GPU Max series GPUs. These exascale and other pre-exascale systems like Summit (also at OLCF) have delivered an unprecedented level of peak performance and storage capacity thus enabling scientific community to solve larger and more complex problems faster than ever. Spectral Element Method (SEM) community has always been among the first to use the state of the art supercomputers in leadership comput-

ing facilities for numerical solution of Partial Differential Equations (PDEs). In this paper we are going to revisit an important aspect of high performance SEM simulations – problem decomposition or simply known as partitioning in the context of the changing HPC landscape in the exascale era. SEM simulations are performed on a mesh, which is a triangulation of the computational domain using smaller regions called elements. These elements usually have the shape of quadrilaterals in 2D and hexahedra in 3D. Each element is further discretized into a set of points called nodes or degrees of freedom (DOF). Numerical (or approximate) solution is represented by a set of basis functions defined on these nodal points which provide a polynomial approximation to the solution of the underlying PDE. Number of elements, E, and the order of basis functions, N , are selected based on the geometric complexity of the computational domain and the accuracy requirements of the simulation respectively. Size of the problem, n, measured in terms of the number of degrees of freedom (DOF) is equal to EN d where d (usually 2 or 3) is the physical dimension of the problem. In this paper, our focus is on d = 3 since it has the highest practical importance in science and engineering applications. Figure 1 shows a timeline of important events in high performance SEM simulations and the evolution of the problem size, n. Introduction of Message Passing Interface (MPI) [2] in 1994 was a major milestone in high performance scientific computing since it standardized the use of distributed memory parallel computing with explicit message passing. MPI has since become the de facto standard for parallel computing and is used in almost all leadership class computing facilities enabling users to utilize thousands of processors in parallel. In 2008, a SEM simulation using 520 million DOF was performed on 65536 processors on BG/P machine at ALCF [3] using Nek5000 [4] – a highly scalable spectral element solver optimized for the CPUs. In 2015, another Nek5000 simulation had 2 billion DOF on 1 million processors on ALCF’s Blue Gene/Q machine [5]. Once the Summit computer at OLCF was commissioned with NVIDIA V100 GPUs, NekRS – a GPU port of

B

= n

60

2B

2015

2020

=

=

2010

, 38 P = 0B 2 , P 764 = 8 72 00 0

1M =

,P

0M =

2005

n

2000

n

1995

n

M

PI

52

1. 0 To do : G

or d

,P

on

=

B

el l

65 53 6

Pr iz e

Nek5000 using OCCA [6] library was used to perform a 60 billion DOF simulation using 27648 GPUs [7]. Then in 2023, NekRS was used in a massive simulation with 380 billion DOF on 72000 GPUs on the first exascale system – Frontier at OLCF. It is clear that the problem sizes will soon reach Trillion DOF on future exascale systems with 109 – 1010 spectral elements. Also, it is worth noting that the local problem size, n/P , is significantly larger on the GPU based simulations compared to the CPU based simulations.

2025

Figure 1: Timeline of important events in high performance scientific computing and the evolution of the problem size measured in degrees of freedom (DOF).

the total time. Since the arithmetic time is proportional to the local problem size, the minimum of the maximum arithmetic time is achieved when the problem is evenly distributed among the processors. Thus, load balance is important to reduce the total time. Minimizing the communication time is more complicated since it depends on the communication volume and the number of messages. The cost of the two terms in Equation 1.2 is equal when α/β = Wp /Mp . It is easy to see m2 = α/β is the size of message that would take as twice as long to transmit as it would take to transmit a single word. If the average message size, Wp /Mp is less than m2 , then the communication time is dominated by the number of messages, otherwise, it is dominated by the communication volume. We will demonstrate that for the GPU based supercomputers, the average message size will be much greater than m2 more often than not and the communication time will be dominated by the communication volume. This is due to the fact that GPU based systems usually have a much larger local problem size and thus larger total message sizes. On traditional CPU systems where the local problem was much lower, the communication time was dominated by the number of messages. So a partitioning method which reduces the communication volume is essential for exascale SEM simulations. The remainder of the paper is organized as follows. We define the graph partitioning problem on a mesh in Section 2 and then look at the related work in Section 3 followed by a description of our approach Sections 4–7. We present results in Section 8 followed by a discussion of future directions in Section 9.

With problem sizes getting larger and larger, an important metric for a user of a HPC system is how fast one can get to the solution i.e., time-to-solution. In the context of parallel SEM simulations using P processors, time spent on the local computation by processor p, Ta,p is proportional to the local problem size np on that processor and the arithmetic time ta – which is usually measured as the time taken to perform a single floating point operation (FLOP) in small matrix-matrix multiplications. If F FLOPs are performed per DOF by the numerical method to find the approximate solution, 2 Problem Definition then the arithmetic time is given by: Since we want to partition mesh elements to processors, we consider the dual graph G(V, E) of the original (1.1) Ta,p = F · np · ta mesh. Here V is the set of vertices of G where each PP where p=1 np = n. When it comes to communication vertex represents a spectral element in the original time, Tc,p , there are two main factors – total message mesh. E is the set of edges in G which represent size (or volume) Wp measured in 64-bit words and the connections between spectral elements. Edges exist the number of messages Mp . Then using the postal between elements which share a vertex, edge, or face in model [5], the communication time can be approximated the original mesh. Given a graph G(V, E) with non-negative edge as: weights ω : E → R>0 and P processors, the graph (1.2) Tc,p = αMp + βWp partitioning problem as it relates to mesh partitioning is to find a set of subsets V1 , V2 , . . . VP , such that: where α is the latency in seconds and β is the inverse bandwidth in seconds per word. Total time is then given (2.4) V ∪ V ∪ . . . V = V 1 2 P by the sum of the arithmetic and communication times: (2.5) Vi ∩ Vj = ∅, i ̸= j, 1 ≤ i, j ≤ P (1.3) T (n, P ) = max [Ta,p ] + max [Tc,p ] (2.6) max | Vi | − | Vj | = 1, 1 ≤ i, j ≤ P 1≤p≤P

1≤p≤P

Based on 1.3, it is clear that we have to minimize both where | Vi | is the size of the subset Vi . The first arithmetic and communication time in order to reduce two constraints ensure that subsets Vi , 1 ≤ i ≤ P

are a partition of V and that the subsets are disjoint (i.e., each element is assigned to exactly one processor). Our weight function ω is the number of shared vertices between two elements connected by the edge e ∈ E. This could be 1 if the two elements share a vertex, 2 if they share an edge, or 4 if they share a face in the original mesh. 3

Related Work

Finding the optimal partition (either bisection or a kway partition) is known to be NP-complete. Optimal solutions are rarely obtained and many heuristics have therefore been developed to address this problem, which is important in many fields [8]. These heuristics try to group elements based on similarity measures like spatial position of the mesh elements or how strongly the elements are connected. The former set of methods are known as geometric methods and the latter are known as combinatorial methods. Spectral methods on the other hand solve an optimization problem that tries to minimize the edge-cut between partitions using the global connectivity information of the graph. Geometric methods make the partitioning decision based on the physical coordinates of the graph nodes in space ignoring the graph connectivity. Such methods are applicable only when nodes of the graph under consideration are associated with a physical space or when the connectivity information can be used to embed the graph nodes in a k−dimensional space [8]. Since the SEM/FEM meshes represent a computational domain in physical space, it is straightforward to apply these methods. Because connectivity information is ignored, geometric methods tend to produce partitions with lower quality than combinatorial or spectral partitioning methods. Recursive Coordinate Bisection (RCB) seeks to reduce the communication by trying to minimize the length of the boundary-cut between partitions. To minimize the boundary, RCB finds the longest coordinate axis and bisects the graph normal to that axis. This approach avoids cutting along the longest axis. RCB is extremely fast due to the fact that we only need a fast sorting algorithm to bisect the graph in parallel. Apart from giving low quality partitions, RCB can only partition along the normal coordinate axes. This can affect the quality of partitions as the longest boundary may not be a normal coordinate axis, for example when the mesh is oriented at an angle to coordinate axes [8]. Recursive Inertial Bisection fixes this issue by calculating the inertial axis of the mesh, projecting the element coordinates to this axis and then ordering elements by this projected coordinate. Both of these schemes can be improved by rescaling the geometry such that, on

average, the elements are of the same diameter, which avoids having skewing the long axis by high-aspect-ratio elements [8–11]. Both RCB and RIB are based on the projection of coordinate values in a single axis. Space-filling curves try to improve over RCB/RIB algorithms by positioning center of mass of mesh elements in a continuous curve in higher dimensions (2D, 3D). After elements are placed in a curve, that ordering can be used to partition into the required number of subdomains [12]. Combinatorial methods try to use the connectivity information to group highly connected vertices, thus producing partitions with lower edge-cuts that tend to have a lower number of disconnected components compared to geometric methods. Levelized Nested Dissection (LND) starts with a single vertex and tries to find a subdomain which contains half of the graph nodes where all the element of the subdomain are connected. LND achieves this objective by gradually growing subdomains through the addition of adjacent vertices. Although LND produces better partitions, it is slower than geometric methods since the algorithm doesn’t have the same level of parallelism as the previous approaches [13]. Spectral methods try to solve a relaxed version of a discrete optimization problem. Let G(V, E) be a graph with | V | vertices. Graph bisection can be thought of as splitting | V | vertices into roughly equal sets A and B in size while minimizing edge-cuts. Let’s consider a vector x of size n where each entry in x is either +1 if ith vertex of graph belongs to A or −1 if it belongs to B. If L is the unweighted Laplacian associated with the graph (which we will define in Section 4), then xT Lx is equal to 4δ(A, B) where δ(A, B) is the number of edges shared between A and B [14]. Now minimizing edge-cuts is equal to minimizing xT Lx subject to the constraint xi = ±1. This minimization problem itself is NP-Complete and can’t be expected to be solved exactly. The problem can be relaxed by allowing each component xi to vary continuously and the minimizer for the relaxed problem is the eigenvector corresponding to the smallest positive eigenvalue of L [14]. Parallel implementations of these algorithms are available in several software packages. parMetis is a Message Passing Interface (MPI) based parallel library that implements a variety of algorithms for partitioning unstructured graphs [15]. The main algorithm in parMetis is a multilevel k-way partition algorithm which works by first coarsening the input graph to a smaller graph with a few hundred vertices and then projecting the computed k-way partition for this smaller graph to the original graph. PT-Scotch also uses a similar multilevel approach but with nested dissection at the top

level to partition the projected separator set rather than k-way partitioning [16]. They also use additional band refinement heuristics to refine the identified projected separators. Other notable software packages include JOSTLE [17] which also performs multilevel partitioning and Zoltan [18] which performs hypergraph partitioning. While there is apparently a wealth of partitioners available, not all of them provide the level of control and scalability required for targeted exascale simulations using the SEM. Because of the SEM granularity, where each element typically has 500 to 2000 degreesof-freedom, we need strict load balance. Furthermore, the partitioner must scale to P = 106 –107 ranks for current and future platforms such as Frontier, Aurora, etc. Users of the SEM, generally engineers and physicists, need an end-to-end solution where they provide minimal data input (e.g., a mesh); derived quantities such as partitions should be dealt with automatically. The partitioner must support multi-material domains, for example, such as in conjugate heat transfer where there is an expensive (i.e., flow) region that dominates the work per element, coupled to a typically less expensive region where one only needs to solve for a scalar field (temperature). Partitioning both these regions to achieve reasonable load balance and minimal communication requires detailed control over the partitioning process. There are multiple reasons, therefore, to revisit the parallel partitioning question on modern architectures. In this paper we introduce parRSB, a parallel graph partitioner based on Recursive Spectral Bisection (RSB) described in Section 4 to achieve load balance, minimal communication volume and detailed control over the whole partitioning process. Our main contributions include space-efficient storage and evaluation of the Laplacian based on a scalable gather-scatter operator, augmented preconditioned projection for the inverse iteration, approximate Krylov-subspace projection of the inverse iterates, algebraic multigrid preconditioner (AMG) for solving the Laplacian system and use of RCB as a pre-partitioner to speed up the partitioning process. 4

Recursive Spectral Bisection

The basis of Recursive Spectral Bisection (RSB) is a graph bisection algorithm first developed by Pothen et al. [19]. The main idea of their work is to make the bisection decision based on the eigenvectors of the Laplacian matrix L of the graph G(V, E) where G is the dual graph defined in Section 1. The dimension of L is |V |×|V | where |V | is the number of spectral elements in the mesh. The unweighted Laplacian matrix L = (lij )

is defined as follows:  if (vi , vj ) ∈ E  −1 deg(vi ) if i = j (4.7) lij =  0 otherwise Here, deg(vi ) is the degree of vertex vi in the graph. In other terms, it is the number of neighbors of corresponding element ei in the mesh. Definition (4.7) is for the unweighted Laplacian since connections between elements are not weighted, i.e., each connection between two elements shows up in the Laplacian matrix as a −1 irrespective of whether they share a vertex, edge, or face. In other words, we have a constant (unit) weight function in this case. In SEM/FEM meshes, it makes sense to put more weight to the latter connection since if those two elements end up on different processors, the communication volume between the processors is higher than if the elements shared just a single vertex. Hence, we use the following weighted form of the Laplacian which tends to give better partitions in practice:  if (vi , vj ) ∈ E  −ω P ij if i = j (4.8) lij = i̸=j ωij  0 otherwise where ωij is the weight associated with the edge (vi , vj ). We define the weight function ω : E → R>0 associated with the graph G(V, E) to be the number of vertices shared between elements (ei , ej ) corresponding to (vi , vj ) in G(V, E). For example, the weights for a 3D hex mesh would be 1, 2, or 4, depending on whether the connection is by vertex, edge, or face. The Laplacian matrix can be written as L = D−A, where D is the degree matrix and A is the adjacency matrix. D is a diagonal matrix that contains the degree of each vertex for the unweighted Laplacian and the sum of weights for the weighted Laplacian. For the unweighted Laplacian, ij th entry of A is 1 if there is an edge between the elements ei and ej and zero otherwise. For the weighted Laplacian, the value of these non-zero elements of A is equal to the weight. The Laplacian matrix reveals interesting spectral properties of the corresponding graph due to its close relation to the Laplacian operator in partial differential equations (PDEs) [20]. Spectral properties of the Laplacian matrix have been studied by many authors, notably by Fiedler in [21, 22]. It is easy to see that the Laplacian matrix is positive semi-definite and hence all the eigenvalues are greater than or equal to zero. The smallest eigenvalue of the Laplacian matrix is zero with a corresponding eigenvector of all ones (1). This is because the row-sum of both weighted and unweighted

Laplacians is zero. Let the eigenvalues of L be ordered Algorithm 1 Recursive Spectral Bisection (RSB) as follows: 1: Find y2 using Lanczos or Inverse Power Iteration 2: Sort mesh elements according to y2 (4.9) λ1 = 0 < λ2 ≤ λ3 . . . ≤ λk 3: Assign 1st half of elements to 1st half of processors 4: Assign 2nd half of elements to 2nd half of processors where k = |V |. The multiplicity of λ1 is equal to the 5: Recurse on two halves till they are entirely connumber of connected components of the graph. Since all tained in a single processor the SEM/FEM meshes are connected, the multiplicity of λ1 is 1. λ2 , known as the algebraic connectivity, is related to the edge and vertex connectivities [19]. λ2 is a good measure of how well the graph is connected [23]. Let the eigenvector associated with λ2 be y 2 . This is known as the Fiedler vector. Fiedler studied the partitions generated by the sorted components of y 2 and proved that vertices of the graph G can be partitioned into connected subgraphs based on the order of component values. Based on this result, a graph bisection is obtained by first sorting the components of y 2 and then assigning the elements corresponding to the first half of the components to the first half of processors and the rest of the elments to other half of processors. y 2 can be found by a suitable numerical method for solving the following eigenvalue problem and selecting the eigenvector associated with the second smallest eigenvalue: (4.10)

Lx = λx

Since we know that the eigenvector associated with the smallest eigenvalue is all ones (1), we can restate (4.10) as finding the eigenvector corresponding to the smallest eigenvalue which satisfies the following conditions: (4.11)

Lx = λx 1T x = 0

We could use inverse power iteration to solve (4.11) with an initial vector orthogonalized with respect to 1. We can use a preconditioned conjugate gradient iteration to solve the Laplacian system at each iteration of inverse power iteration with LAMG as the preconditioner [24]. Inverse power iteration is expensive since it involves solving a linear system at each iteration. An alternative method is to use Lanczos, which can be used to find the extreme eigenvalues and associated eigenvectors. We have implemented both of these methods, with multigrid-preconditioned conjugate gradient iteration to solve Laplacian system in inverse power iteration. Once y 2 is found, the graph is bisected as described above. Partitioning continues recursively in each half till each partition is entirely contained in a single processor. Algorithm 1 outlines the RSB algorithm. Going down the recursive tree in Algorithm 1, if we run into a situation where the number of processors is not divisible by two, we can create two partitions with number of processors in each of them differ by 1. In

general, If P is the number of processors at any given level, we can create two partitions with sizes ⌊ P2 ⌋ and ⌈ P2 ⌉. Elements in the mesh will also be distributed between two partitions proportional to the number of processors in each partition. 5

Evaluating the Laplacian

Both Lanczos and inverse power iteration require performing the action of the Laplacian on an input vector. The most straightforward way to do this is constructing the Laplacian matrix explicitly and store it in a sparse matrix format like Compressed Sparse Row (CSR), with rows distributed among different processors. Then a sparse matrix vector product can be performed simply by gathering required vector components in each processor. The other option is to partially assemble Laplacian simply by keeping track which values need to be summed up while performing the action of Laplacian on the vector. We call this the gather-scatter approach. In SEM/FEM applications, the gather-scatter operation is a central part of the matrix-assembly and load-vector assembly processes. In the PDE context, for a given mesh with m elements, one has mL = m × v local vertices that, for linear elements, correspond to local degrees-of-freedom. These local degrees-of-freedom allow for easy, element-based, evaluation of derivatives and integrals. However, for continuous finite element approximations, the true (or global) degrees-of-freedom are associated with the shared nodes—local nodal values are simply copies of these global values. Suppose there are mu of the unique global nodal values, u = {uj }, enumerated j = 1, . . . , mu . Then a map (or copy) from the global to representation can nominally be effected by the matrix-vector product, (5.12)

uL

=

Qu,

where Q is a sparse mL × mu Boolean matrix that has a single nonzero in each row. Each column of Q corresponds to a single entry in u and the number of nonzeros in that column is equal to the number of elements that share the global id associated with that column. We refer to (5.12) as the scatter operation—data is scattered to the elements from the global degrees-offreedom. Equally important in SEM/FEM applications

Figure 2: Application of (left) QT P and (right) P T Q. Blue arrows indicate the application of P and P T , which is local to each element (and hence to each processor), while red arrows indicate the application of Q and QT , which require inter-processor communication in cases where adjacent elements are not on the same processor. At the end, one has y = P T Qw = P T QQT P x. Q and P are copy operations, while QT and P T involve summation. Global vertex numbers are indicated near the v = 4 vertices of each element. is the gather operation, w = QT uL . In this case, local values are summed to their global counterparts. Figure 2 illustrates the application of QT (left) and Q (right) in red. In the gather operation (QT ) values on local nodes that are shared by more than one element are mapped, with summation, to their global counterparts (denoted by ⊕ in the left figure). If this operation is followed by a scatter (Q), these summed values (wj in the right figure) are copied back to the originating local vertices. For the graph Laplacian on the dual of the vertex graph (i.e., on the graph associated with element centroids), we need to map one input value, xe , from each element to one output value, ye , on each element. This operation can be expressed with a local copy matrix, P that first maps xe for element e to its v vertices. At that point, one applies the nonlocal gather-scatter, QQT , followed by P T which sums the v values back to ye . QT P is illustrated on the left and P T Q on the right of Fig. 2. For the dual graph, the weighted adjacency matrix can be represented as

the graph Laplacian. To understand how this approach leads to the weighted Laplacian, consider a case with x1 =1 and x2 = x3 = 0. We see that z = Aw x will generate z2 = 2 because the value at x1 has two paths to reach element 2. Singletons, such as global vertex 1, will make no contribution to y = Lx because their contribution is cancelled out by the difference Dw − Aw . The advantage of the gather-scatter approach is that it can be implemented with minimal setup cost, which is particularly important given that, at the partitioning stage, the problem setup has yet to be well defined on the processors. The simulation starts with what is effectively a random distribution of elements assigned to each processor, where nothing is known about how the graph is mapped to the processors. (As discussed in Section 8, we have found it advantageous to use an initial pass of RCB to at least organize the data into contiguous subsets in order to minimize gather-scatter communication overhead during Lanczos or inverse iteration.) To set up the QQT communication, we use the gather-scatter library, gslib [25], which involves two calls, gs setup, and gs op. In the setup phase, each processor provides an mp × v list of global integer (long) Aw = P T QQT P. pointers corresponding to the global id of each vertex on each of the mp elements on rank p. In a discovThe weighted degree matrix Dw =diag(de ), with d = ery phase, gs setup discerns which ranks have commonly Aw 1 representing the vector containing the row-sums of shared vertices, which vertices are multiply represented of Aw . From these two matrices we have L = Dw −Aw as

within a rank, and which are singletons. Subsequently, one calls gs op with a list of mp × v doubles (or floats), to which QQT is applied. The setup phase executes in O(log P ) time and typically runs multiple trials of gs op using either a pairwise exchange, a generalized all-to-all (based on the crystal-router of [26]), or an all-reduce, choosing whichever is the fastest. Users can also prescribe a prefered communication algorithm if they wish to reduce the setup overhead. We illustrate the gs setup call for the the example of Fig. 2 under the assumption that elements 1 and 2 are on rank 0 while element 3 is on rank 1. On respective ranks, the setup differs only in the enumeration of the global ids,

there is a shared face. If we can number edges and faces of the mesh uniquely, we can setup the gather-scatter such that contributions from each neighbor element is counted exactly once. It turns out that it is very easy and fast to do this numbering as we have a global numbering for vertices already available. 6

Lanczos Iteration

The Lanczos algorithm is preferred for finding eigenvalues of real symmetric matrices since the number of floating point operations per iteration is smaller than other equivalent methods like Arnoldi [27]. We implemented Lanczos with restarts to find y 2 . Suppose Lanczos algorithm takes j number of iterations before it reaches one of terminal conditions (maxrank 0: global id=[ 1 2 3 4 2 5 4 6 ]; mL = 8; imum number of iterations or the tolerance). Lanczos mL = 4; algorithm produces as output a series of vectors q and rank 1: global id=[ 3 4 7 8 ]; i scalars αi and βi where 1 ≤ i ≤ j. The vectors q i are Thus defined, respective setup and execution of QQT known as Lanczos vectors and are orthonormal. Let S would be expressed on each processor by the commands be the n × j matrix which has the Lanczos vectors as the columns. gs handle = gs setup(global num,mL );   gs err = gs op(gs handle,w); | | | S =  q1 q2 . . . qj  The gs op call in this case generates the result w := | | | QQT w. (We remark that gslib has scaled to millions of ranks on the IBM BG/Qs, Mira and Sequoia, with The scalar series α form the diagonal of a symmetric j setup times typically measured in fractions of a second. tri-diagonal matrix T while β form the subdiagonal j For 1 billion global vertices, the setup time on 1048576 (and superdiagonal) of T. ranks of Mira was 1.24 seconds, including 20 trial runs   of gs op. Setup times for current examples on Summit α1 β2 are provided in Section 8.)  β2 α2 β3    Evaluating the unweighted Laplacian is more comT= . . . .. .. ..    plex than the weighted Laplacian since we must avoid βj αj adding contributions from the same neighbor element multiple times. In other words, connection between two The eigenvalues of T are known as Ritz values. If t , elements should be counted only once irrespective how 1 ≤ i ≤ j are the eigenvectors of T, the vectors y = Sti i i many vertices they share. The unique set of neighbors of are called the Ritz vectors. Ritz values and vectors are an element can be found by a simple counting principle. the approximations for the eigenvalues and eigenvectors In order to do so, we need to think of the connectivity of L. We can find an approximation to y by finding the 2 between two elements in terms of vertices, edges and eigenvector corresponding to minimum eigenvalue in T faces. and then multiplying it by S. When two elements are connected in a SEM mesh, they are either connected by a vertex, edge or a face 7 Inverse Iteration in 3D. In 2D, we only have vertices and edges. The weight of the connection is the total number of vertices Algorithm 2 Inverse Iteration they share. A neighbor which shares an edge is counted 1: Initialize b with random values twice and a neighbor which shares a face is counted four 2: for i ← 1 to j do times in the weighted Laplacian. We can adjust the 3: Orthogonalize b w.r.t 1 gather sum for edges by removing contributions once 4: b = b/∥b∥ for each shared edge. An undesirable side effect of this 5: Solve Ly2 = b subtraction is that the contributions from face neighbors 6: b = y 2 are eliminated since they have 4 vertices and 4 edges 7: end for in common. So we need to add back contributions if

We implemented inverse iteration (Algorithm 2) using a preconditioned flexible conjugate gradient (flexcg) method, augmented with full projection to solve for Fiedler vector once we reached larger problem sizes for which the Lanczos didn’t converge. A key point of our flexcg implementation is that the initial search direction is not preconditioned. This choice is motivated by the fact that we are solving for y 2 , which is an eigenvector of L and that the preceding iterate is used as the right-hand side of Lx = b. As b −→ y 2 , the Krylov subspace (in L, but not in M −1 L) will be invariant and the modified flexcg will return in a single iteration, at which point we stop the outer inverse-iteration. We used an aggregation-based algebraic multigrid (AMG) preconditioner inspired by the Lean Algebraic Multigrid [24] method. Pseudocode for the multigrid preconditioner is shown in Algorithm 3. Algorithm 3 V-Cycle Input: σ, nsmooth , i, r, L Output: u 1: D = 1.0/diag(L) 2: u = σ Dr 3: r = r − Lu 4: for s ← 1 to nsmooth do 5: u = u + σDr 6: r = r − LσDr 7: end for i+1 i 8: Lc = Ji LJi+1 i+1 9: rc = Ji r 10: ec = V-Cycle(σ, nsmooth , i, rc , Lc ) i 11: u = u + Ji+1 ec 12: for s ← 1 to nsmooth do 13: r = r − Lu 14: u = u + σDr 15: end for The finest level of multigrid, level l0 , uses the Laplacian corresponding to the actual element mesh, which is evaluated using either the CSR or gather-scatter representation. For coarser levels (l1 , l2 , . . .), the gatherscatter based Laplacian evaluation is not straightforward to setup as we don’t have a global vertex numbering readily available once the elements are coarsened. For those levels, we forget about the geometric aspect of coarsening and view it as simply applying a restriction operator Jii+1 going from level i to i + 1 and its transi i pose Ji+1 to generate Li+1 = Jii+1 Li Ji+1 , where Jii+1 is a Boolean matrix corresponding to piecewise-constant prolongation. It is easy to see that this preserves the qualities of the Laplacian. We bootstrap the prolongation operator from an RCB ordering of the mesh elements. Before setting up

the multigrid V -cycle, we partition the mesh using RCB. Denote the coarsest level as lN , which comprises all elements. The next level comprises the two sets arising from the initial coordinate bisection, such that the N −1 prolongation matrix is JN = [ 1 1 ]T . If I2 is the ranki−1 i 2 identity matrix, then Ji = I2 ⊗Ji+1 , up to the finest level, where the number of entries in L0 is typically not a power of 2. Note that L1 = J01 L0 J10 will be approximately half as large as L0 because aggregating element pairs at level l0 is equal to condensing (i.e., summing) the rows and columns of L0 corresponding to those pairs. Our aggregation strategy is local to a processor till we reach a single row per processor and we keep number of processors constant till we reach that threshold. Once we reach that threshold, we turn off the number of processors by half at each level since the size of the Laplacian matrix is reduced by half. This process is described in detail below. 1. Assume that at the level l, we have a nl × nl matrix and Pl active processors. Also, Laplacian Ll is represented using CSR format. 2. If nl = Pl , we set Pl+1 = Pl /2, otherwise Pl+1 = Pl . We will keep processors p1 , . . . pPl+1 active in level l + 1. 3. Gather rows 1, . . . nl of level l, distributed over Pl processors, to Pl+1 ≤ Pl processors. (If Pl = Pl+1 , no data movement takes place.) 4. Collapse these rows locally on Pl+1 processors to get nl+1 rows where nl+1 ≈ nl /2. (The orginal rows are preserved; the collapsed outputs are written to a new row id in level l + 1. 5. Repeat the above for columns as well. 6. These steps produce a mapping of row and column ids from level l to l + 1 which can be used to setup i Jii+1 and Ji+1 . With the above steps, we generate L0 , L1 , L2 , . . . as CSR matrices. 8

Results

We tested parRSB as a mesh partitioner for NekRS [28] – a highly scalable GPU-accelerated spectral element solver. Tests were done on Summit and Frontier supercomputers at OLCF where parRSB was run on the CPU. Running parRSB on CPU is not a bottleneck for NekRS since parRSB finished in less than a minute during all our tests. In the context of the SEM simulations where the setup time could easily take

minutes due to various preprocessing steps and the preconditioner setup, the time taken by parRSB is negligible. In SEM production runs, total setup time is only a fraction of the time taken by the solver since the actual simulation time could easily be hours or days. We made sure that the partitions generated by parRSB are load balanced as much as possible with the maximum load imbalanced being only a single element as mentioned in Equation 2.4. To evaluate the partition quality, we choose number of neighbors and the average communication volume per neighbor (i.e., total outgoing message size from a processor divided by the number of neighbors). As pointed out in Section 1, these two metrics are responsible for the communication. We used two pebble bed meshes with 13 million and 99 million spectral elements for comparing Lanczos and inverse iteration on Summit supercomputer. For both Lanczos and inverse iteration, we evaluated Laplacian operator using the gather-scatter approach listed in Section 5. Gather-scatter setup times for the 99 million element mesh was about 0.15s for 115 billion vertices using 37800 processors. For the 13 million mesh, it took about 0.037s for 15 billion vertices with 11340 processors. These setup times are a fraction of the time taken to run either Lanczos or inverse iteration. Timing data and number of neighbors for the 13 million element mesh is shown in Table 1 as the number of processors is increased from 4872 to 11340. First column of the table is the number of processors and second column is the time taken to finish the partitioning with Lanczos. We found that running a geometric partitioning scheme like RCB or RIB before RSB can significantly reduce time taken by the RSB algorithm. Lanczos runtime reduced approximately by 2× when RCB is used as a pre-partitioner. This gain is more pronounced at the higher levels of the bisection tree where the global partition sizes are larger. As we go down the bisection tree running RCB or RIB as a pre-partitioner has diminishing returns. Third column is the time taken to finish the partitioning with RCB as a pre-partitioner for Lanczos. Fourth and fifth columns show maximum number of neighbors and the average neighbors respectively. Table 2 shows the same information for preconditioned inverse iteration. These results are comparable to Lanczos (in Table 1) except when it comes to execution time. Main culprit for slowing down inverse iteration is the AMG preconditioner setup, which has to be done at each level in RSB tree, not the inverse iteration itself. For the mesh under consideration, for P = 4872, both Lanczos and inverse iteration spent about 2s finding the Fiedler vector for the first cut. Inverse iteration took about 6 iterations to find the cut compared to

Lanczos which reached the maximum number of restarts (50) permitted in our implementation. Table 3 list results for 99 million element pebble-bed mesh. We can see that the average number of neighbors stay more or less the same as the number of processors are increased. Run time for Lanczos for the 99 million element mesh is higher than the 13 Million element mesh due to the larger local problem sizes. Table 4 shows results for a weak scaling study on the Frontier supercomputer using cube meshes as the number of processors is increased from 8 to 65536 while the number of elements per processor is kept constant at ≈ 8000. We can see that parRSB generates partitions with maximum and average number of neighbors in the expected range (26) for a SEM mesh. What is more interesting is the average message size reported in the last column. Given that the m2 for Frontier is about 5000 (TODO: Cite ping pong tests), based on the numbers reported for the average message size, we can see that the most of these runs fall on the regime where communication volume is the dominant factor in communication time. 9

Future Work

It is not uncommon in graph Laplacians to have eigenvalues with multiplicty greater than 1. For example a checkerboard would have a Fiedler-vector pair that shares the same minimum nontrivial eigenvector. In fact, the space is much larger, because things only need to be topologically a checkerboard to have this property. So, a rectangular mesh with the same number of elements in x and y would have a matched pair of Fiedler vectors, and any linear combination of these eigenvectors is a valid Fiedler vector in the RSB algorithm. Unfortunately, the quality of the cut is not the same for all vectors. Cutting along the vertical or horizontal midline of an N × N graph will expose N faces; cutting on a 45o angle would expose ≈ 2N faces, which is suboptimal. Conceivably, one could explore the one dimensional parameter space spanned by θv 1 + (1 − θ)v 2 , which constitutes all linear combinations of the eigenvector pair to try to find the min-cut solution. Unfortunately, neither standard Lanczos nor inverse iteration will reveal a second eigenvector in this pair— what emerges is merely the same linear combination of those eigenvectors that was in the initial seed vector for the Lanczos/inverse iteration. To find the second eigenvector, one would need to run a block, or subspace, iteration comprising two vectors that are orthogonalized on each iteration. Fortunately, the overhead for performing such a block iteration is unlikely to be large given that these algorithms are dominated by relatively sparse (i.e., latency-limited) communication. Moreover, gslib

P 4872 6468 8106 9744 11340

Lanczos Time (s) 45.5 36.8 30.5 30.6 29.1

RCB + Lanczos Time(s) 24.8 20.4 17.2 16.6 16.4

Max Neighbors 25 26 27 30 27

Average Neighbors 16.0 16.5 17.0 16.8 17.0

Table 1: Partition time and number of neighbors for 13 Million elements mesh using Lanczos. P 4872 6468 8106 9744 11340

Inverse Time (s) 32.8 30.0 28.8 30.0 28.5

Max Neighbors 28 27 35 32 29

Average Neighbors 16.2 16.6 17.2 17.0 17.1

Table 2: Partition time and number of neighbors for 13 Million elements using preconditioned inverse iteration. already supports vector-based communication, that is, multiple values at each node of the graph. In cases limited by interprocessor latency, the cost of communicating extra values at each node is close to nil. In addition, a block algorithm would amortize the setup overhead of the graph and of the aggregation-based AMG for inverse iteration. We believe that such an algorithm might be of value in generating higher-quality RSB-based cuts and are pursuing these ideas with that goal in mind. There are multiple ways to reduce setup cost of the AMG preconditioner. 108 element meshes are just a start and future exascale computers will be running billion element meshes. For SEM/FEM simulations, CPU based supercomputers like Mira at Argonne Leadership Computing Facility used to run efficiently with 2000 grid pointers per rank. With GPGPUS, Summit at OLCF usually requires about 2 million grid points per rank in order to run at an acceptable efficiency. We expect this trend to continue and the local probem size in exascale systems most likely be larger than that of Summit. For a larger local problem size, it will take longer for parRSB to just run on CPUs. So porting parRSB to use accelerators is in our roadmap. References [1] Top500. Top500 list - june 2023, 2023. [2] W.D. Gropp, E. Lusk, and A. Skjellum. Using MPI: Portable Parallel Programming with the MessagePassing Interface, 2nd edition. MIT Press, Cambridge, MA, 1999. [3] P.F. Fischer, J. Lottes, W.D. Pointer, and A. Siegel. Petascale algorithms for reactor hydrodynamics. J. Phys. Conf. Series, 125:012076, 2008.

[4] Paul Fischer, James Lottes, and Henry Tufo. Nek5000. [Computer Software] https://doi.org/10.11578/dc. 20210416.29, jun 2007. [5] Paul F Fischer. Scaling limits for pde-based simulation. In 22nd AIAA Computational Fluid Dynamics Conference, page 3049, 2015. [6] David S Medina, Amik St-Cyr, and Tim Warburton. Occa: A unified approach to multi-threading languages. arXiv preprint arXiv:1403.0968, 2014. [7] Paul Fischer, Stefan Kerkemeier, Misun Min, YuHsiang Lan, Malachi Phillips, Thilina Rathnayake, Elia Merzari, Ananias Tomboulides, Ali Karakus, Noel Chalmers, et al. Nekrs, a gpu-accelerated spectral element navier–stokes solver. Parallel Computing, 114:102982, 2022. [8] Kirk Schloegel, George Karypis, and Vipin Kumar. Graph partitioning for high performance scientific simulations. 2000. [9] Michael T Heath and Padma Raghavan. A cartesian parallel nested dissection algorithm. SIAM Journal on Matrix Analysis and Applications, 16(1):235–253, 1995. [10] John R Gilbert, Gary L Miller, and Shang-Hua Teng. Geometric mesh partitioning: Implementation and experiments. SIAM Journal on Scientific Computing, 19(6):2091–2110, 1998. [11] Marsha J Berger and Shahid H Bokhari. A partitioning strategy for nonuniform problems on multiprocessors. IEEE Transactions on Computers, 36(05):570– 580, 1987. [12] John R Pilkington and Scott B Baden. Dynamic partitioning of non-uniform structured workloads with spacefilling curves. IEEE Transactions on parallel and Distributed Systems, 7(3):288–300, 1996. [13] Alan George and Joseph W Liu. Computer solution of large sparse positive definite. Prentice Hall Professional Technical Reference, 1981. [14] Alex Pothen. Graph partitioning algorithms with applications to scientific computing. In Parallel Numeri-

P 16212 21588 27006 37800

RCB + Lanczos Time (s) 67.2 57.7 50.4 43.6

Max Neighbors 31 27 26 29

Average Neighbors 17.2 16.3 16.6 16.8

Table 3: Partition time and number of neighbors for 99 Million elements mesh using Lanczos. P 8 512 1024 4096 8192 16384 32768 65536

RCB + Lanczos Time (s) 0.1 1.5 6.0 6.0 20.0 20.0 20.0 41.0

Max Neighbors 7 26 23 26 24 25 26 27

Average Neighbors 7.0 19.8 14.0 23.0 15.6 15.8 24.0 16.0

Average Message Size 8580 5600 9890 5130 8680 9092 4890 8427

Table 4: Partition time and quality for a weak scaling study with E/P ≈ 8000. cal Algorithms, pages 323–368. Springer, 1997. [15] George Karypis, Kirk Schloegel, and Vipin Kumar. Parmetis: Parallel graph partitioning and sparse matrix ordering library. 1997. [16] Cédric Chevalier and François Pellegrini. Pt-scotch: A tool for efficient parallel graph ordering. Parallel computing, 34(6-8):318–331, 2008. [17] Chris Walshaw and Mark Cross. Jostle: parallel multilevel graph-partitioning software–an overview. Mesh partitioning techniques and domain decomposition techniques, 10:27–58, 2007. [18] Karen D Devine, Erik G Boman, Robert T Heaphy, Rob H Bisseling, and Umit V Catalyurek. Parallel hypergraph partitioning for scientific computing. In Proceedings 20th IEEE International Parallel & Distributed Processing Symposium, pages 10–pp. IEEE, 2006. [19] Alex Pothen, Horst D Simon, and Kang-Pu Liou. Partitioning sparse matrices with eigenvectors of graphs. SIAM journal on matrix analysis and applications, 11(3):430–452, 1990. [20] Stephen T Barnard and Horst D Simon. Fast multilevel implementation of recursive spectral bisection for partitioning unstructured problems. Concurrency: Practice and experience, 6(2):101–117, 1994. [21] Miroslav Fiedler. Algebraic connectivity of graphs. Czechoslovak mathematical journal, 23(2):298–305, 1973. [22] Miroslav Fiedler. A property of eigenvectors of nonnegative symmetric matrices and its application to graph theory. Czechoslovak mathematical journal, 25(4):619– 633, 1975. [23] Nair Maria Maia De Abreu. Old and new results on algebraic connectivity of graphs. Linear algebra and its applications, 423(1):53–73, 2007. [24] Oren E Livne and Achi Brandt. Lean algebraic

multigrid (lamg): Fast graph laplacian linear solver. SIAM Journal on Scientific Computing, 34(4):B499– B522, 2012. [25] gslib sparse communication library. https://github. com/Nek5000/gslib. Accessed: 2021-10-03. [26] G. C. Fox, M. A. Johnson, G. A. Lyzenga, S. W. Otto, J. K. Salmon, and D. W. Walker. Solving Problems on Concurrent Processors. Prentice-Hall, Englewood Cliffs, NJ, 1988. [27] Kesheng Wu and Horst Simon. Thick-restart lanczos method for large symmetric eigenvalue problems. SIAM Journal on Matrix Analysis and Applications, 22(2):602–616, 2000. [28] Paul Fischer, Stefan Kerkemeier, Misun Min, YuHsiang Lan, Malachi Phillips, Thilina Rathnayake, Elia Merzari, Ananias Tomboulides, Ali Karakus, Noel Chalmers, et al. Nekrs, a gpu-accelerated spectral element navier-stokes solver. arXiv preprint arXiv:2104.05829, 2021.