Conceptio › Archive › arXiv CS

arXiv CSopen access

A Simple Communication Scheme for Distributed Fast Multipole Methods

2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

Open Source ↗Direct PDF ↓

cloud distributed-computing parallel-computing

distributed computing, parallel computing, cloud

arXiv:2604.19243v1 [cs.DC] 21 Apr 2026

A Simple Communication Scheme for Distributed Fast Multipole Methods

The

International

Journal

High

Performance

Computing

Applications XX(X):2–40 ©The Author(s) 2025 Reprints and permission: sagepub.co.uk/journalsPermissions.nav DOI: 10.1177/ToBeAssigned www.sagepub.com/

SAGE

Srinath Kailasa1, 2

Abstract We present a simple hierarchical communication scheme for distributed Fast Multipole Methods (FMMs) based on MPI neighborhood collectives and uniform trees. The method targets the common case of extending an existing high-performance shared-memory uniform-tree FMM implementation to distributed memory with minimal redesign while preserving any shared memory optimizations optimizations. Benchmarks on the ARCHER2 supercomputer demonstrate that our method can scale to very large problem sizes, we demonstrate weak-scaling up to 3.2e10 uniformly distributed points on 512 nodes of the machine in our largest runs. Our simplifications based on uniform trees result asymptotically optimal communication under practical constraints that we document. We record departure from this scaling for non-uniform points, however we still obtain practically useful runtimes due to the ability to retain our shared memory optimizations.

Keywords Fast Multipole Method, FMM, MPI, Distributed Memory, ARCHER2, High Performance Computing, HPC, Distributed FMM Prepared using sagej.cls [Version: 2017/01/17 v1.20]

The International Journal of High Performance Computing Applications

XX(X)

Introduction

Fast algorithms for computing long-range interactions between large numbers of particles arise in a wide range of scientific applications, including gravitational N -body simulations, electrostatics, molecular dynamics, and boundary integral methods for elliptic PDEs. These problems typically involve computing sums of the form f (xi ) =

N ∑

K(xi , yj )sj ,

i = 1, . . . , M,

(1)

j=1 d where {yj }N j=1 ⊂ R are ‘source’ points with associated densities sj ∈ F, and d {xi }M i=1 ⊂ R are the ‘target’ points at which we wish to evaluate the field,

F ∈ {C, R} and d is the spatial dimension. The kernel K(·, ·) : Rd × Rd → F describes the interaction between source and target points. The prototypical kernel of interest in examining N -body problems is the three-dimensional Laplace kernel, K(·, ·) : R3 × R3 → R, K(x, y) =

1 (d = 3), ∥x − y∥

(2)

which models gravitational and electrostatic potentials. Direct evaluation of (1) requires O (N M ) operations, making it computationally prohibitive at large scales. Moreover, (1) implies a global data dependency between all source and target points. However, when the kernel admits a lowrank structure for well-separated point sets, hierarchical algorithms such as the Fast Multipole Method (FMM) (Greengard and Rokhlin 1987) or the closely related H-matrices (Hackbusch and Nowak 1989; Hackbusch 1999) can reduce the computational cost to O (N + M ) in the best case. These methods work by recursively partitioning the computational domain using hierarchical data structures (commonly octrees for d = 3 or quadtrees for d = 2), and splitting the sum in (1) into near-field and far-field contributions

1 Graphcore, United Kingdom 2 Department of Engineering, University of Cambridge, United Kingdom

Corresponding author: Srinath Kailasa, Graphcore, 11-19 Wine Street, Bristol, England BS1 2PH, UK. Email: [email protected]

Prepared using sagej.cls

Kailasa

for each box in the hierarchy containing target particles: f (xi ) =

∑

K(xi , yj )sj +

yj ∈Near(xi )

∑

K(xi , yj )sj .

(3)

yj ∈Far(xi )

The near-field interactions are evaluated directly, while the far-field interactions are approximated using compressed representations. Despite reducing algorithmic complexity, hierarchical methods do not eliminate the global data dependencies of the original N -body problem. In distributed memory, a naive implementation typically results in expensive all( ) to-all communication, with worst-case cost O P 2 , where P is the number of processes. Communication complexity therefore limits parallel scalability and must be considered when adapting N -body software to current and emerging exascale hardware. A landmark result by Lashuk et al. (Lashuk et al. 2009) demonstrated an approach (√ ) that reduced the runtime communication complexity for FMMs to O P , where P is the number of processes. More recently, Ibeid et al. (Ibeid et al. 2016) showed that communication can in principle scale as low as O (log P ) when the communication pattern mirrors the hierarchical structure of the algorithm. However, for practical distributed implementations, asymptotic communication complexity is only one part of the story. Communication schemes must also be simple to implement, portable across MPI libraries and machine topologies, and compatible with the shared-memory optimizations critical for single-node performance. This is particularly relevant on modern HPC systems, where hardware is deeply hierarchical (Kim et al. 2008) and overall performance is strongly controlled by the computational intensity of the underlying shared-memory kernels (Williams et al. 2009). In this setting, a method that relies only on common collective operations and a precomputed communication layout while preserving eﬀicient local kernels can be more valuable than one that pursues the strongest asymptotic bound at the cost of substantial implementation complexity. In this paper, we contribute: 1. A communication strategy for distributed FMM using only standard MPI4 collective operations and a precomputed communication layout, with neighborhood collectives for local exchange, designed to preserve sharedmemory optimizations developed for single-node FMMs and to align with existing software using uniform trees and spatial encoding. Prepared using sagej.cls

The International Journal of High Performance Computing Applications

XX(X) 2. A communication analysis of the distributed FMM phase of our ( ) implementation showing a favorable bound of O log P + (N/P )2/3 under explicit practical assumptions, and identifying the conditions under which the logarithmic global term is expected to be relevant. 3. Large-scale benchmarks on ARCHER2 studying weak scaling for both uniform and non-uniform distributions, together with a practical strongscaling study. We find our method capable of practically computing large problem sizes. Despite not adhering to asymptotic optimality in communication in all cases, we are able to compute the FMM over large numbers (3.2e10) of uniform/non-uniform points in practically useful times (of the order of 10s of seconds). We provide an overview of the FMM in distributed memory in Section 2, with a high-level overview of the sequential algorithm in Section 2.1, the importance of shared memory eﬀiciency for distributed implementations in Section 2.2, and a review of efforts to extend the algorithm to distributed memory in Section 2.3 to build the required context of our own contributions. Section 3 describes our approach, Section 4 presents a theoretical complexity analysis, and Section 5 reports our performance results on the CPU-based ARCHER2 system hosted at the Edinburgh Parallel Computing Centre (EPCC). In Section 6, we contextualize these results relative to existing work, and conclude with future directions in Section 7.

2 2.1

Fast Multipole Methods High Level Overview of the FMM

The FMM has come to embody a wide variety of closely related methods, which though sharing a common algorithmic structure, differ in their approach to field approximation and compression of the far-field in (3). Additionally, implementation details such as the use of hardware accelerators, parallel runtimes, and programming languages lead to a further diversity in variant FMMs. To avoid this complexity, we describe the structure common to all techniques, and defer to the literature and sources therein for details of both the algorithmic structure (Ying et al. 2004; Greengard and Rokhlin 1987; Fong and Darve 2009; Darve 2000; Anderson 1992; Martinsson and Rokhlin 2007) in addition to approaches for the optimization of the implementation, which often depend closely on the chosen algorithmic approach (Kailasa et al. 2024; Prepared using sagej.cls

Kailasa

Malhotra and Biros 2015; Lashuk et al. 2009; Wang et al. 2021; Yokota and Barba 2012; Ltaief and Yokota 2014; AbdulJabbar et al. 2014; Gumerov and Duraiswami 2008; Chandramowlishwaran et al. 2010; Blanchard et al. 2015; Messner et al. 2012; Bramas 2020). Most modern implementations of FMMs are built using hierarchical data structures called Hashed Octrees (HOTs) for problems in three spatial dimensions (or Hashed Quadtrees (HQTs) in two dimensions), pioneered by Warren and Salmon (Warren and Salmon 1993). This includes the most recent exafmm variant (Wang et al. 2021), pvfmm (Malhotra and Biros 2015), and our own kifmm-rs (Kailasa 2025a). Alternative approaches based on graph bisection (AbdulJabbar et al. 2014) and methods specialized for task-based runtimes (Blanchard et al. 2015) also exist. We do not discuss their communication schemes in detail because they are closely tied to the underlying data structures and runtimes used in those cases. For ease of reference, we employ the notation of Lashuk et al. (Lashuk et al. 2009), summarized in Table 1. Consider a problem in R3 . We enclose all source and target points in a cube of suﬀicient size, then construct an octree T by refining the cube into eight octants, or boxes, which are recursively subdivided until a user-defined depth, d. This corresponds to an FMM with uniform refinement, in which each leaf box is the same size. We associate a level, l ∈ [0, d], with each box β in T , where l = d is the deepest level and l = 0 corresponds to the root level occupied by the bounding cube. Adaptive refinement is also possible, in which adjacent leaf boxes (that is, boxes that share a vertex, edge, or face) may differ in size. In that case boxes are refined until no leaf box contains more than q sources/targets. We denote the set of leaf boxes by L ⊂ T . In either case, we associate with each β ∈ T four so-called interaction lists. These are lists of boxes whose interactions with each β are coupled during the course of the FMM algorithm, called the U, V, X and W lists in the terminology of (Ying et al. 2004), which we illustrate in Figure 1. In the uniform case, only the U and V lists are defined. We begin by defining the colleagues, C(β), of β as boxes which are adjacent to β and at the same level. If β is a leaf box such that it is not recursively subdivided into more boxes, its U list consists of its colleagues, C(β), in addition to β itself. For all boxes β ∈ T , the V list consists of boxes at the same level whose parents are colleagues of β’s parent. Prepared using sagej.cls

The International Journal of High Performance Computing Applications

XX(X) (a) Adaptive Tree

U List

V List

(b) Uniform Tree

W List

X List

Figure 1. Interaction lists for a box β ∈ T , where T is a quadtree for a problem in R2 , shown for clarity. Panel (a) shows an adaptively refined tree and panel (b) a uniformly refined tree.

In the adaptive case, the U list of a leaf box, β, now consists of all adjacent boxes, not just colleagues, since adjacent boxes may not be at the same level. The definition of the V list remains the same. In addition we define W and X lists for β. The W list consists of boxes that are descendants of colleagues of β and are themselves not adjacent to β. The X list consists of boxes for which β appears in their own W list. Each box β ∈ T is associated with two vectors, u and d. With reference to Figure 2, the vector u is understood as a compressed representation of the field due to densities contained within the box, which is accurate outside the region covered by C(β). The vector d is understood as a compressed representation of the field due to densities located outside C(β), and therefore describes the field evaluated at target points contained within β due to source densities in its far-field. Additionally, leaf boxes, β ∈ L, are also associated with vectors x – the target points contained in β, y – the source points contained in β, s – the associated source densities, and f – the field evaluated at target points in β. When using HOTs or HQTs, boxes are encoded using a scheme that preserves spatial locality, most commonly Morton or Hilbert keys. These are typically stored as 64-bit integers that represent the relative positions of boxes within the domain, with three bits encoding the box’s location at each level of refinement. A direct mapping of these keys to memory addresses performs poorly for adaptive point distributions, since many boxes are empty. Therefore, keys corresponding to non-empty boxes are hashed so that adjacent hashes correspond to adjacent Prepared using sagej.cls

Kailasa

Figure 2. We illustrate a problem in R2 . The left plot shows the well-separated blue region in the exterior of an orange source box where the representation u for the source box is valid. The right plot shows the blue region of validity in the interior of a target box where d summarizes the field from well-separated source boxes in the orange exterior. The arrows represent the outgoing and incoming nature of the representations u and d respectively. Adapted from Figure 2.1 of (Kailasa 2025b).

memory addresses. We use Morton keys in our implementation (Sundar et al. 2008); when sorted, they correspond to a pre-order traversal of the tree and therefore provide an eﬀicient way to index data associated with each box during the recursive algorithm, while also ensuring that spatially local data corresponds to geometrically local boxes. The reduced computational complexity of the FMM is due to the interaction lists and recursive tree structure, which enable the field, f , to be evaluated at each leaf box, β ∈ L, by coupling interactions between β and only a bounded number of other boxes. Linear complexity is obtained as we refine the hierarchical data structure such that there are O (N ) leaf boxes. The sequential algorithm, designed for shared memory, is summarized for the uniform and adaptive cases in Algorithm 1. It relies on a number of linear operators, described in Table 1, which translate between the representations of the field due to s, described by u and d, and the final evaluated field f . The structure and construction of these operators is the main point of difference between various FMMs, with corresponding trade-offs in implementation, performance, and accuracy; some of these trade-offs are summarized in (Yokota 2013). Construction methods range from algebraic approaches, based on the construction and compression of blocks of a matrix corresponding to (1), to purely analytical techniques that employ potential theory to derive kernelspecific formulae. Our software kifmm-rs (Kailasa 2025a), which we use for the benchmarking reported in this paper is based on the algebraic Kernel Independent Fast Multipole Method (kiFMM) (Ying et al. 2004). Prepared using sagej.cls

The International Journal of High Performance Computing Applications

XX(X) Table 1: Summary of notation. Box Data β, α

Unique representation of a box

Target points (leaf boxes only)

Source points (leaf boxes only)

Source densities at source points (leaf boxes only)

Potential at target points (leaf boxes only)

Multipole representation

Local representation Operators

Sβ

Source-to-Up (S2U) or Particle-to-Multipole (P2M) operator

UP(β),β

Up-to-Up (U2U) or Multipole-to-Multipole (M2M) operator

Tβ,α

Up-to-Down (U2D) or Multipole-to-Local (M2L) operator

Dβ,P(β)

Down-to-Down (D2D) or Local-to-Local (L2L) operator

Qβ,α

Source-to-Down (S2D) or Particle-to-Local (P2L) operator

Rβ,α

Up-to-Target (U2T) or Multipole-to-Particle (M2P) operator

Eβ

Down-To-Target (D2T) or Local-to-Particle (L2P) operator

Kβ,α

Source-To-Target (S2T) or Particle-to-Particle (P2P) operator Tree Data

T L

Hierarchical Tree Leaf boxes of T where L ⊂ T ; boxes β ∈ L are defined by not being further subdivided

A(β)

All ancestor boxes of β ∈ T

L(β)

Tree level of a box β ∈ T

C(β)

Colleagues of a box β ∈ T

P(β)

Parent of a box β ∈ T

K(β)

Children of a box β ∈ T

N(β)

The near field of a box β ∈ T , defined by a halo of boxes α such that they are contained in C(β) where L(α) = L(β)

F(β)

The far field of a box β ∈ T , defined as a complement of T with respect to N(β)

U(β) V(β)

Boxes adjacent to box β, including β itself Boxes α ∈ T such that α ∈ K(C(P(β))), L(α) = L(β), and α and β are not adjacent

W(β)

Boxes α ∈ T such that α ∈ N(β) and the u of α converges within β

X(β)

Boxes α ∈ T such that β ∈ N(α), and the densities s are evaluated to construct d for β

Adapted from Table 1 (Lashuk et al. 2009).

Prepared using sagej.cls

Kailasa

Algorithm 1 (Shared Memory FMM): We adapt the compact specification of Algorithm 1 in (Lashuk et al. 2009), but highlight steps unique to the adaptive FMM. oct Input: {xi , yi , si }N i=1 - An octree with |T | = Noct octants. Nleaf Output: {fi }i=1 - the potentials evaluated at target points at each of |L| = Nleaf leaf octants. // (1) The upward pass, traverses the tree, T , in post-order - from finest to coarsest boxes, level by level ∀β ∈ L : uβ = Sβ sβ // (1a) S2U:source-to-up step ∀β ∈ T : uP (β) + = UP (β),β uβ // (1b) U2U: up-to-up step // (2) The downward pass, traverses the tree, T , in pre-order - from the coarsest to the finest boxes, level by level ∀β ∈ T : ∀α ∈ V(β): dβ + = Tβ,α uβ // (2a) VLI: V-list interactions if T is adaptive then ∀β ∈ T : ∀α ∈ X(β): dβ + = Qβ,α sα // (2b) XLI: X-list interactions end if ∀β ∈ T : dβ + = Dβ,P (β) dP (β) // (2c) D2D: down-to-down step // (3) The leaf level boxes are finally handled, where the far-field is accumulated, and the near field directly summed. if T is adaptive then ∀β ∈ L: ∀α ∈ W(β): fβ + = Rβ,α uα // (3a) WLI: W-list interactions end if ∀β ∈ L: fβ + = Eβ dβ // (3b) D2T: down-to-target step ∀β ∈ L: ∀α ∈ U (β): fβ + = Kβ,α sα // (3c) ULI: U-list interactions (direct kernel evaluations)

2.2

Shared Memory Optimizations and Their Role in Distributed FMMs

In this section, we review key memory and parallelism considerations, which we aim to preserve when scaling to distributed systems. Examining Algorithm 1, we highlight opportunities for shared memory optimizations. Firstly in terms of coarse-grained parallelism, we note that the (exact) ULI step is independent from the remainder of the upward/downward passes for approximate field calculations, and can be performed asynchronously. Furthermore the VLI and XLI steps are independent level-wise, after the upward pass. Within each VLI interaction for boxes β ∈ T at a given level l during the downward pass, many interactions across source and target boxes share translational symmetry and share a common Up-to-Down (U2D) operator during the VLI step, this is also true for U2D operations within the V-list of Prepared using sagej.cls

The International Journal of High Performance Computing Applications

XX(X) a single target box, and techniques that exploit this symmetry are at the core of (Kailasa et al. 2024; Messner et al. 2012; Gumerov and Duraiswami 2008; Malhotra and Biros 2015) – as they enable optimizations designed to reduce cache misses through the sharing of data amongst multiple U2D calculations. As already mentioned, fine grained parallelism is exposed for the ULI step, which is well suited for modern SIMD or SIMT machines. Coalescing of the U2U and D2D operations is also possible across multiple sets of sibling octants (those which share a common parent), however is less significant for runtime as memory accesses are spatially local, corresponding to the geometric locality of the corresponding boxes involved in these operations. With reference to Figure 1 coupling interactions between a box β and boxes in its various interaction list naturally leads to cache-ineﬀicient memory accesses over non-contiguously stored data in a HOT approach. The eﬀicient handling of these memory accesses is often the greatest barrier to achieving high computational intensity in the shared memory portion of an FMM implementation. The memory hierarchies of modern processors impose significant costs on memory accesses that result from cache misses, and often determine the end-to-end performance of an implementation. As a result the VLI, XLI and WLI steps in Algorithm 1 require careful design to ensure that implementations are not limited by the slowness of ineﬀicient memory access (Messner et al. 2012; Kailasa et al. 2024; Malhotra and Biros 2015; Takahashi et al. 2012; Gumerov and Duraiswami 2008; Fong and Darve 2009). The FMM was originally introduced to avoid direct kernel evaluations using (1), and minimize the ULI step, however this step is highly suited to the structure of modern SIMD and SIMT processors – due to the high rate of data re-use available for each target particle in (1). Therefore, the performance of implementations is controlled by the depth of the associated tree. A deeper tree corresponding to a smaller number of q particles for β ∈ L, resulting in more VLI and XLI steps, which are bandwidth bound, against shallower trees which result in larger ULI steps which we can write with high operational intensity and are compute bound. Indeed, these steps, in particular the VLI and ULI steps, dominate the total runtime of an FMM, often exceeding 90 % of wall time (Chandramowlishwaran et al. 2010; Kailasa et al. 2024). In the uniform case, one can take advantage of stencil based methods to re-form the VLI step such that spatial and temporal locality is maintained in memory access, and have been explored in (Kailasa et al. 2024; Gumerov and Duraiswami 2008; Messner et al. 2012; Malhotra and Biros 2015; Takahashi Prepared using sagej.cls

Kailasa

et al. 2012). This results from the fact that the U and V are simply halos around each box and known a priori, and simply checked for existence at runtime (i.e. whether or not an expected box in U, or V, contains any particle data). This allows for significant memory optimizations, such as the pre-allocation of contiguous buffers to store box data required during the ULI and VLI steps and index maps to rapidly look up data associated with interactions lists for a given target box, this is crucial for achieving high-performance in shared memory (Kailasa et al. 2024). However, in the adaptive case interaction lists depend on the particular particle distribution — and must be calculated at runtime. Because of this, we favor shallow uniformly refined trees that minimize memory-bound interactions and delegate irregularity to compute-bound kernel evaluations, where SIMD acceleration is highly effective as shown in our previous work (Kailasa et al. 2024), or SIMT based techniques for GPUs as shown in (Lashuk et al. 2009). We justify this approach for the distributed memory case by experiment in Section 5.

2.3

Distributed Memory Parallelism

We review the communication requirements of the distributed FMM when HOTs/HQTs are used as the underlying hierarchical data structures. We describe the algorithm in three spatial dimensions using HOTs. The input to a distributed FMM is a set of N source points and N target points, which without loss of generality we take to be the same set, together with N densities associated with the source points. These points are input in arbitrary order, so each process i ∈ [0, P − 1] initially stores approximately N/P points, where P is the total number of processes used. By ‘process’ we mean a coarse-grained unit of computation with its own memory space and execution state, such as an MPI process. The output is a distributed vector f , where the portion fi on each process contains the evaluated field for the target points stored on process i.

2.3.1

Tree Construction and Domain Partitioning When using HOTs we follow

the approach for the construction of distributed linear octrees first described in (Sundar et al. 2008). Here, the computational domain Ω is partitioned across ∪P processes such that Ω = i=1 Ωi , where each subdomain Ωi is defined by the volume covered by the leaf boxes assigned to it. We illustrate such a tree in Figure 3. We define the Locally Essential Tree (LET) for each process as, Prepared using sagej.cls

The International Journal of High Performance Computing Applications

XX(X) Level 0 Level 1 Level 2

Level d

1 ..... 8

........

9 ..... 18

Process 1

...................

Process 2

8d-2

8d-1

Process P

Figure 3. We illustrate the partition of a computational domain Ω discretized with a uniform linear Hashed Octree (HOT) discretized to a depth d, such that each process i ∈ [0, P − 1] is responsible for a subdomain Ωi defined by the leaf boxes such that Ω ∪ is formed from their disjoint union Ω = P i=1 Ωi . We illustrate that a spatial encoding scheme that preserves geometric locality amongst boxes leads to a partition which reflects this spatial locality in the logical proximity of neighboring processes.

∪

LET(i) = β∈[Li

∪

I(β).

(4)

A(Li )]

This is the subtree of the global HOT containing the boxes it requires to compute its local contribution to the field, fi , at each process. It contains a union of data associated with the combined interaction lists, I(β), of all the leaves, Li , local to that process in addition to all of their collective ancestors A(Li ). The LET is seen to be the minimum set of box data required to evaluate fi . 2.3.2

Communication Phases in Distributed FMM When using HOTs, the leaf

boxes are found by performing a spatial encoding; here we describe a scheme based on Morton keys (Sundar et al. 2008). A similar approach would apply to alternative spatial encoding schemes. We assign each point a Morton key encoded at the maximum level, or depth d. This is pleasingly parallel over all processes, and as such requires O(N ) work. Once these keys have been found, there are five major communication-intensive steps, which can be distinguished by whether the communication is executed during setup or during runtime of the FMM algorithm:

Prepared using sagej.cls

Kailasa

Box 2.1: Distributed Communication Steps 1. (Setup): A distributed sort of the input points based on their spatial encoding, such that geometrically close points are assigned to logically close processes. 2. (Setup): The construction of LETs on each process to construct the interaction lists for boxes under their control. 3. (Setup): The communication of densities, s, required during the ULI step (in the uniform case), as well as the XLI step (in the adaptive case). This will necessarily be between logically close processes due to Step (1). 4. (Setup): Load balancing input to balance workload among partitions. 5. (Runtime): The communication of the u vectors calculated during the upward pass to processes that require this data during the VLI steps (in the uniform case), as well as the WLI steps (in the adaptive case). 2.3.3

Distributed Sorting Strategies A simple approach for Step (1) in Box 2.1

is provided by Samplesort (Blelloch et al. 1991). Here, we identify P buckets associated with each process such that all keys in the i’th bucket are less than or equal to all keys in the (i + 1)’th bucket for i ∈ [1, P ]. We distribute the leaves L into appropriate buckets and sort the local buckets, Li , to end up with a globally sorted L. In our software (Kailasa 2025a), we use a randomized scheme for selecting P − 1 splitters to define these buckets. Here, the Morton keys of b points are selected from each of i ∈ [1, P ] processes and gathered at a nominated process. Once a local sort has been applied to this subset of entries, one can select P − 1 entries by choosing every b’th entry as a splitter. The splitters can be used to form an all-to-all communication between all P processes to redistribute local portions of point data into their appropriate bucket. Using modern implementations of MPI the dominant costs of the sort have a complexity of, ) ( N N O tr ( log ) + ti b log P + ti m̄P , P P

(5)

where the first term describes the cost of the local sorts at each node, each initialized with O (N/P ) keys, with tr describing the slowness of local memory Prepared using sagej.cls

The International Journal of High Performance Computing Applications

XX(X) accesses as tr = 1/bandwidth. The second term describes the cost of gathering samples of size b from each process using schemes based on binomial trees or recursive doubling, where ti describes the cost of communication over the interconnect as ti = 1/bandwidth. The final term describes the cost of the allto-all operation required to distribute entries into their respective buckets, where m̄ is the mean message size. This corresponds to an MPI_Alltoallv operation, and for large message sizes can be accomplished in O (P ) communication steps in modern MPI implementations (Kang et al. 2020). As such, this forms the most expensive part of our implementation. However, as our focus is on the runtime phases of the FMM in this paper we treat this cost as a fixed startup cost that is paid only once when setting up the tree. We acknowledge that asymptotically better algorithms exist Sundar et al. (2013). 2.3.4

Load Balancing Our implementation targets uniform and near uniform

point distributions, such that we can uniform trees. In which case The Central Limit Theorem (CLT) implies that the sample quantiles used as splitters the parallel sort step increasingly approximate the true quantiles as the sample size b grows. Specifically, the standard deviation of the splitter positions decreases √ proportionally to 1/ b. This ensures that the buckets formed by these splitters contain approximately equal numbers of entries, with deviations from the ideal size N/P diminishing with larger sample sizes. In the context of the uniform FMM the estimated buckets correspond to partitions of the points such that each bucket contains an approximately equal number of points. For uniform or approximately uniform point distributions, such as points distributed evenly throughout the problem domain or uniformly on a surface without sharp features, the proposed random sampling approach therefore leads to an approximately load balanced partition of points across MPI processes that improves with the number of samples. The uniformity of the distribution means that this also corresponds to an approximately load balanced distribution of leaf boxes Li ⊂ L across processes i ∈ [0, P − 1] This analysis does not to non-uniform point distributions, and in this case we expect a significant load imbalance. In our experiments we demonstrate that practically useful large scale computations are still achievable in practice despite this – without the need for further specialized load-balancing algorithms. 2.3.5

Runtime Communication Scheme of Lashuk et al Lashuk et al. (Lashuk

et al. 2009) split runtime communication in Step (5) in Box 2.1 into two parts, in the first part they perform the upward pass independently over each LET, Prepared using sagej.cls

Kailasa

acknowledging that for boxes close to the root of the tree their u vectors will only contain partial results, which must be synchronized across all processes. In order to avoid a global all-to-all, they use a customized ‘reduce-scatter’ scheme based on point-to-point communication, discussed in Algorithm 3 of (Lashuk et al. 2009), using a classical hypercube communication scheme. This algorithm has a communication complexity of, ( ) √ O tl log P + ti m̄(3 P − 2) ,

(6)

where tl is a constant associated with the latency of establishing communication between two processes over log P communication rounds, and ti is constant related to the bandwidth constraints of a HPC system of sending a mean of m̄ messages between processes during each round. The bandwidth limited nature of modern HPC systems leads to the dominant cost being √ ( )2/3 proportional to P . In the uniform case m̄ can be estimated as N , with an P argument relating the surface area to volume ratio of the subdomain controlled by each process, Ωi , which is proportional to the octant data required by other processes. This leads to an estimate of the dominant cost in communication complexity being, (( O

N P

)2/3

√

) (7)

Once completed the downward passes can be performed independently on each process to calculate fi for their local portion of target points. This scheme maintains any shared memory optimizations developed for both the upward and downward passes, however introduces redundant computations for d at higher levels of the tree during the VLI step – which are repeated across all processes which share many common ancestor boxes close to the root. Data dependencies are explicitly resolved at runtime, including the discovery of coupled interactions and the associated allocations required to store interaction list box data. 2.3.6

Runtime Communication Scheme of Ibeid et al Ibeid et al. (Ibeid et al.

2016) describe how the runtime communication complexity of the distributed FMM can be reduced to ( O log P +

Prepared using sagej.cls

(

N P

)2/3 ) ,

(8)

The International Journal of High Performance Computing Applications

XX(X) which has been shown to apply to both uniform and adaptive cases (Yokota et al. 2014). This reduction is achieved by splitting the HOT into two hierarchical components: (1) a global or process tree and (2) a local tree within each process. The idea is illustrated in Figure 4 for the uniform case. In this scheme, each leaf of the global tree corresponds to a process, and serves as the root of a finegrained local tree owned by that process. For adaptive trees, the leaves of the global tree are defined uniformly at the same level, while adaptivity is captured entirely within the local trees. To understand how this splitting reduces redundant computation near the root, consider the upward pass. After each process computes the upward representation u for the root of its local tree, the U2U (upward-to-upward) operation needs to be performed across local roots. However only one process needs to carry out this computation per group of eight processes – whose local roots correspond to sibling boxes (i.e. those sharing a common parent in the global tree). These processes can transmit their u data to a single nominated process responsible for the U2U at that level. This logic can be applied recursively up the global tree, involving only 8l active processes at level l, and yielding a total communication complexity of O (log P ), since each process participates in a constant number of communications and the height of the global tree can be seen to be log P from Figure 4. A similar pattern applies during the downward pass over the global tree. With reference to Figure 4, for each level l ∈ [2, dglobal − 1], the u data needed for the VLI step at a set of sibling boxes is determined by the interaction list K(C(P(α))), where α is their shared parent. The nominated process, which already holds the u data from the upward pass, can then broadcast or exchange with at most 26 processes whose data overlaps with this interaction list. Since each process communicates with a bounded number of neighbors, this pattern again yields a total communication cost of O (log P ) derived from the height of the global tree. Communication for the downward pass within the local trees involves only O (1) neighboring processes. Regardless of whether the tree is uniform or adaptive, the u and s ghost data required for the local U, V, W, and X list interactions is restricted to boxes contained within, or overlapping with, those defined by K(C(P(β))), where β is a box at levels [dglobal + 1, dglobal + dlocal ] within a local tree (see Figure 4). These interactions remain confined to local trees rooted at boxes adjacent to the root of the current local tree, and therefore involve communication with at most 26 neighboring processes. Prepared using sagej.cls

Kailasa

Level 0 Level 1 Level 2

1 ..... 8

........

9 ..... 18

Global Tree Level dglobal - 2 Level dglobal - 1 Level dglobal Level dglobal + 1

1 ..... 8

........

dglobal = log8(P)

........ ........

9 ..... 18

Local Tree Level dglobal + dlocal

dlocal = log8(N/P)

Process i

Figure 4. An illustration of octrees distributed over P processes using Ibeid et al.’s specification (Ibeid et al. 2016). Each node in the graph corresponds to a box, and each edge corresponds to a parent-child relationship, omitting some nodes for clarity. The root box of the global tree is labelled G, and corresponds to a box covering the problem domain. The root boxes of each local tree are labelled Pi for i ∈ [1, P ] where P is the total number of processes. Each node is labelled with its index when sorted in Morton order, at each level l ∈ [0, d], where d is the tree depth. Each process contains a non-overlapping subset of the leaf boxes. Adapted from Figure 3.7 in (Kailasa 2025b)

Algorithm

Our algorithm iterates on the communication pattern described by (Ibeid et al. 2016), with a specific focus on practical engineering concerns to ensure highperformance, portability by relying only on standard MPI collectives and ease of implementation in an existing single-node FMM software. We restrict ourselves to uniform trees, which simplifies the communication graph, and allows us to identify and resolve all data dependencies before any tree traversal. We delegate to a highly optimized ULI implementation to handle irregular point distributions. 1. We create a static communication graph at setup time, and describe it using MPI’s neighborhood communicators – introduced in MPI-3∗ . As a ∗ Introduced in MPI-3, neighborhood communicators allow processes to define logical neighbors

within a topology (e.g., a Cartesian grid or a graph-based partition tailored to an application), enabling eﬀicient, localized communication with those neighbors. While the communicators themselves define logical relationships between MPI processes, MPI implementations can optimize communication by leveraging information about the physical hardware topology,

Prepared using sagej.cls

The International Journal of High Performance Computing Applications

XX(X) result our implementation shares conceptual similarity with task-based runtime approaches. 2. We evaluate the upward and downward passes on the global portion of the tree on a single process, thus ensuring we maintain our crucial VLI optimizations developed for single-node performance, with the communication at the global level of the algorithm reduced to a single all-gather and all-scatter. 3. Except for the non-uniform case at the very largest scales evaluating 32e9 particles, communication scaling is broadly supported by our complexity analysis. The description is split into two parts. Section 3.1 describes the steps required to establish data dependencies between all boxes in the tree, and encode the static communication graph in Neighborhood communicators. Section 3.2 describes the runtime execution of our algorithm. We note that neighbourhood communication is particularly well suited to the static, sparse communication graph implicit in Ibeid et al.’s scheme. Each process has only a limited number of neighbors with which it must exchange data: at most 26 during the local and global portions of the downward pass, and 8 during the global portion of the upward pass. Because these exchanges arise from adjacent local roots in the tree, the communication pattern is spatially localized in the application domain.

3.1

Algorithm Setup – Precomputing Data Dependencies Box 3.1: Algorithm Setup (1) Tree Construction A user specifies the depth of the global tree dglobal , which is used to encode all locally contained points within a process. These are then sorted with this coarse grained Morton key – the resulting leaves correspond to the local roots. Each process, i ∈ [1, P ], can in principle be responsible for multiple local trees – a level of granularity that can be tuned by adjusting the splitters in the parallel sort. i,j , where Each process builds its own set of local trees, denoted Tlocal j ∈ [1, . . . , niP ] indexes the local roots at process i ∈ [1, P ]. The Locally

when available, to improve performance. Furthermore, neighborhood communicators support localized variants of collective operations.

Prepared using sagej.cls

Kailasa

Essential Tree (LET) of each local tree is constructed using (4), where boxes in the interaction list are defined from the halos that constitute the U and V interaction lists, those contained on ghost processes are checked for existence as we describe below. (2) Global Tree Layout The set of local roots across all processes defines the structure of the global tree, which we refer to as the layout. The layout maps MPI ranks to the local roots they own. In the uniform case, the number of local roots assigned to a process acts as a proxy for computational load, enabling coarse-grained load balancing by adjusting the number of local trees assigned to each process. To construct the layout, we perform an MPI_Allgatherv over a global communicator to collect local root metadata (e.g., Morton keys and owning ranks) from all processes. Each process then has global knowledge of the layout. i,j (3) Query Identification For each box β ∈ Tlocal , we have computed its U

and V interaction lists in Step (1). To determine whether an interaction list box is locally available or must be fetched from a remote process, each process checks whether the box lies within the subdomain defined by the layout. Since all trees descend from known local roots, this check takes O (P ) per box. In the worst case, the total cost of forming these queries is O (N ) per process, but this can be reduced. In practice, a box’s U and V interactions only involve trees rooted in the 26 neighboring subdomains. Thus, queries need only be formed for boxes potentially owned by these neighbors, reducing the cost to O (26 · (N/P )). (4) Forming Query Packets Each process now constructs two query packets: • The V-list query packet: Morton keys and associated process ranks of all non-local boxes required for VLI interactions. • The U-list query packet: Morton keys and associated process ranks of non-local boxes needed for ULI operations. (5) Building Neighborhood Communicators We now construct neighborhood communicators, one for each of these queries which we call the U-List Communicator and V-List Communicator. We use MPI_Dist_graph_create to encode inter-process dependencies identified in Step (4). Because interaction lists are spatially local in the problem

Prepared using sagej.cls

The International Journal of High Performance Computing Applications

XX(X) domain, the resulting communication graphs are sparse. Each process communicates with at most 26 neighbors – the maximum number of adjacent subdomains in three dimensions. This worst case occurs when each process owns a single local root. If multiple local roots are assigned to the same process, the process-local communication graph is formed from the union of the corresponding local-root halos, so grouping roots onto one process does not increase the worst-case neighbor count beyond the single-root case. Figure 5 illustrates this graph construction for the communicators. (6) Establish Data Dependencies 1. Query Existence Exchange: Each process calls MPI_Neighbor_alltoallv to send its U and V query packets to neighbors via their respective communicators. In the uniform case, each recipient checks whether the requested boxes exist locally (which can be accomplished in constant time if a set is used to store all boxes in the LET) and returns the keys of those that do, with another MPI_Neighbor_alltoallv. 2. U-list Data Exchange: Using the U-list communicator, each process sends the source densities for ghost boxes needed for ULI interactions with a further call to MPI_Neighbor_alltoallv. 3. V-list Buffer Allocation: Using the V-list existence responses, each process allocates buffers to receive the relevant u vectors at runtime. Internally, each process stores its local box metadata using hash-based Morton key indexing, allowing O (1) access to required data. We note that the after setup new source data can be added to existing trees by re-running the U-list data exchange, and clearing existing buffers.

3.2

Algorithm – Runtime Execution

Above we describe the process by which we establish data dependencies for the local trees, but not the global trees. A direct implementation of Ibeid et al.’s scheme (Ibeid et al. 2016) would require neighborhood communication among 8l MPI processes at level l during the upward and downward passes over the Prepared using sagej.cls

Kailasa

global portion of the tree. At higher levels, the communicating processes become increasingly (logically) distant from one another, while each active process contributes only a single u or d vector. Our approach both simplifies the communication steps, and allows us to compute u and d at all global tree boxes using our highly-optimized shared memory kernel implementations. We do this by centralizing the evaluation of the global tree on a single process. We use MPI_Gatherv to gather the u vectors associated with all local roots on a nominated process, run the upward and downward passes on the global portion of the tree using the existing highly optimized shared-memory implementation, and return the d vectors to their associated local roots using MPI_Scatterv. This introduces an artificial bottleneck, where the majority of the system is briefly idle while the global tree is traversed. However this cost is in practice small due to the eﬀiciency of our shared memory implementation (especially without an expensive ULI step required). Additionally, the implementation of the communication for the global tree is reduced to just two global collectives. Having performed the steps in Box 3.1, our runtime algorithm executes in local and global stages with reference to Figure 4: Box 3.2: Runtime Execution (1) Local Stage 1. Run the upward pass for levels l ∈ [dglobal − 1, dglobal + dlocal ] i,j independently on all local trees, Tlocal for i ∈ [1, P ] and j ∈ [1, niP ]

where P is the total number of MPI processes and niP is the number of local trees at process i ∈ [1, P ]. 2. Once the upward pass is complete, the u data required during the downward pass is available for all local trees via the V-list communicator. This data is exchanged using MPI_Neighbor_alltoallv with the V-list communicator and inserted into the buffers, both created during setup. (2) Global Stage 1. Select a nominated MPI process to manage the global tree computation. Gather u data associated with each local root to the nominated process using MPI_Gatherv over the whole network.

Prepared using sagej.cls

The International Journal of High Performance Computing Applications

XX(X) 2. On the nominated process, construct an explicit global tree, Tglobal , and perform the upward and downward passes entirely in shared memory (excluding the S2T operation) for levels l ∈ [0, dglobal − 1]. 3. Scatter the d vectors at the leaf level of the global tree, Tglobal , back to their associated processes using the layout, calculated during algorithm setup, with MPI_Scatterv over the whole network. These expansions correspond to the local roots of each local tree, i,j . Tlocal

4. Each MPI process can now perform the downward pass for its i,j , independently in parallel for levels l ∈ [dglobal − local trees, Tlocal

1, dglobal + dlocal ], with no further communication required.

Complexity Analysis

Our goal in this section is not to argue that asymptotic complexity alone determines the merit of the method, but rather to show that the communication structure proposed above remains compatible with favorable runtime communication bounds in the regimes targeted by our implementation. We therefore also account for the behavior of modern MPI implementations, which commonly use heuristics based on the number of MPI processes and the size of messages to select an appropriate communication algorithm. For example, Figure 6 illustrates the decision tree used by MPICH for MPI_Alltoallv; other implementations use similar strategies. 4.0.1

Algorithm Setup The tree construction during Step (1) in Box 3.1

is dictated by the communication complexity of the parallel sort, which we described in Section 2.3. We treat this sorting phase as a preprocessing step, separate from the distributed FMM phase analyzed below, and in this work we use Samplesort for simplicity of implementation rather than an asymptotically more specialized alternative. Step (2) in Box 3.1 constructs the layout - associating all local roots with MPI ranks and storing at each process. This involves a call to MPI_Allgatherv over the entire network. Each process contributes a small message: the metadata of its local roots, typically Morton keys and ranks, which can be encoded using 64-bit integers. If the total message size is less than ∼ 10 KB per-process, which corresponds to over 1000 local roots per process, modern MPI implementations Prepared using sagej.cls

Kailasa

3 4

Each box corresponds to a local root, labelled with their associated MPI process rank in the range [1, 8] which we color for clarity.

5 4 7

U/V List Communicators deﬁned by colleagues of each local root.

All interactions of 's children are described by children of its colleagues. Therefore U/V list communicators constructed with colleagues of local roots are suﬃcient to capture all child interactions. Consider interactions of 3 :

ULI Interactions (If 3 at leaf level)

VLI Interactions

Figure 5. The first row illustrates an example of a layout of a distributed tree over 8 MPI processes, each box corresponds to a local root, where we use a two dimensional tree for clarity. Colors correspond to MPI processes, each of which can contain multiple local roots. In the second row we highlight a target box β in process 4 and the edges and nodes that must be added to the neighborhood communication graphs at process 4 due to the interactions of β and its descendants. These are calculated from the layout, interaction lists, and nearest neighbors of τ̂ . In the third and fourth rows we consider the required communications for a child node of τ̂ , showing how the neighborhood communication graph created for its parent can be re-used. Prepared using sagej.cls

The International Journal of High Performance Computing Applications

XX(X) can use logarithmic algorithms such as Bruck’s algorithm by default, resulting in a communication complexity of O (log P ) for this step (Thakur et al. 2005). During step (6) in Box 3.1, sub-steps (1) and (2) involve neighborhood exchanges of query packets and ghost data between each process and its spatial neighbors. In the worst case this is bounded by 26 neighboring processes in three dimensions, and does not grow with P . The resulting communication complexity therefore matches that of the local portion of the downward pass described in Sections 5.1.3 and 5.1.4 of (Ibeid et al. 2016), given by ( ) O (N/P )2/3 . 4.0.2

Algorithm

4.0.3

Local Portion Step (2) during the local portion of Box 3.2 again requires

neighborhood communications between O (1) MPI processes for each process as during sub-steps (1) and (2) during step (6) of Box 3.1, resulting in the same complexity estimate. 4.0.4

Global Portion Communication during the global portion in Box 3.2

occurs in steps (1) and (3): a gather of u vectors and a scatter of d vectors associated with local roots. These steps are performed using MPI_Gatherv and MPI_Scatterv, respectively. As mentioned above, when message sizes are small MPI implementations default to logarithmic complexity communication schemes O (log P ). However, if message sizes exceed implementation thresholds, fallback algorithms with O (P ) complexity may be triggered. The message size per process depends on the number of local roots niP , the floating-point precision b ∈ {32, 64}, and the FMM ‘expansion P̃ – which tunes the accuracy ( order’ )

of the evaluated field. This is of O P̃ 2

in 3D FMMs, and in the kiFMM

implementation used by our software (Ying et al. 2004), each u or d vector has length 6(P̃ − 1)2 + 2. The total message size per process is thus approximately Message size ≈ niP · (6(P̃ − 1)2 + 2) · 4.0.5

b bytes. 8

Runtime Communication Bound Assuming message sizes remain below

threshold and spatial communication graphs remain sparse, the runtime communication complexity is ) ( O log P + (N/P ) 2/3 , Prepared using sagej.cls

(9)

Kailasa

Process Count < 8? Yes

Message Size <= 32768B? Yes Spread Out O(P)

Message Size <= 32768B?

No Pairwise O(P)

Yes

Message Size <= 256B? Yes Brucks O(log(P))

Pairwise O(P)

No Spread Out O(P)

Figure 6. The default decision tree in MPICH for choosing an algorithm for MPI_Alltoallv. Adapted from Figure 2 in (Netterville et al. 2022).

which matches the best-case bound of Ibeid et al. (Ibeid et al. 2016). We emphasize, however, that this result should be interpreted as support for the practicality of the design rather than its sole motivation: the main benefit of the method remains that it realizes hierarchical communication using only standard collectives and a precomputed communication layout while preserving the highperformance shared-memory kernels used locally. Achieving this bound in practice depends on the underlying hardware and software stack. In particular, interconnect technologies must support hierarchical collectives, and message sizes must be tuned to remain below system-specific thresholds. If not, global communication may exhibit linear complexity, resulting in an upper bound of ( ) O P + (N/P ) 2/3 .

Benchmarks

The ARCHER2 supercomputer is a HPE Cray EX system, consisting of 5680 compute nodes each containing two AMD EPYC 7742 64-core 2.25GHz CPUs, networked with a HPE Cray Slingshot interconnect in a Dragonfly topology. There are 8 NUMA regions per node, with 16 cores per NUMA region and up to 512 GB of memory per node, or equivalently 4 GB per core. In total there are 750080 CPU cores, supporting 1500160 software threads with hyperthreading. Each node consists of two CPU sockets, each of which contains 8 Core Complex Dies (CCDs); each CCD consists of 2 Core Complexes (CCXs), and each CCX owns 4 compute cores that share 16 MB of L3 cache. This hierarchy makes process placement an important part of performance tuning for Prepared using sagej.cls

The International Journal of High Performance Computing Applications

XX(X) our implementation, since the local portion of the FMM should remain within a shared cache or NUMA region whenever possible. We note that users are able to request at most 512 nodes of the system. We benchmark our software on two point distributions: a uniform random distribution in the volume, and points distributed randomly on the surface of a unit sphere. The spherical distribution is intentionally more challenging for a uniform tree, since many fine boxes are empty and the work associated with the local trees is therefore less balanced than in the uniform case. All results reported here use expansion order P̃ = 3 in single precision for the 3D Laplace problem (2). Each benchmark configuration is repeated five times, and the runtime and eﬀiciency plots report mean wall times together with bands showing one standard deviation across runs. Some run-to-run variability on ARCHER2 is expected even for a fixed configuration, since batch-scheduled runs exhibit system-level variability that is not directly controlled by the application. In the non-uniform case this system-level variability is compounded by the intrinsic imbalance of the workload and messages sent across the network, so the derived parallel-eﬀiciency curves are more variable and we use the runtime curves as the primary indicator of scaling behavior. These experiments are designed to test three claims of the method: first, that preserving shared-memory locality through an appropriate rank granularity matters strongly for performance; second, that a communication scheme based on standard collectives and a precomputed layout still delivers good large-scale runtime performance for the distributed FMM phase; and third, that the centralized global phase remains a secondary cost in practice. To study the effect of node-level granularity we compare three MPI process placement strategies on ARCHER2: 1. CCX pinning: each MPI process is pinned to a single CCX, corresponding to 4 CPU cores that share one L3 cache. 2. CCD pinning: each MPI process is pinned to a CCD, corresponding to 8 CPU cores, or two CCXs, within a single NUMA region. 3. Socket pinning: each MPI process is pinned to an entire socket, corresponding to one AMD EPYC 7742 64-core CPU. At 4e9 points we compare all three strategies directly. Although the algorithm can in principle assign multiple local trees to a rank, the reported ARCHER2 weak-scaling runs use one local tree per MPI rank. The benchmark configurations are chosen so that each node carries the same total number Prepared using sagej.cls

Kailasa

of points despite the differing rank granularity: CCX pinning uses 32 ranks per node with 250000 points per rank, CCD pinning uses 16 ranks per node with 500000 points per rank, and socket pinning uses 2 ranks per node with 4 × 106 points per rank. In all three cases this gives 8 × 106 points per node, so the comparison is made at matched node-level problem size even though the amount of work owned by each rank differs substantially across strategies. These results should therefore be interpreted as a comparison of complete nodelevel configurations rather than of hardware pinning in isolation. Within each pinning strategy, the total number of ranks then increases in factors of eight, so the global tree depth increments by one at each scaling step. The local depth is fixed within a given benchmark configuration, but chosen per strategy to keep the local trees at a reasonable resolution in terms of estimated points per leaf box: the 4e9 campaigns use local depth 4 for CCX and CCD pinning, and local depth 5 for socket pinning. Since all benchmark results use P̃ = 3 in single precision, each rank contributes a single local-root u or d vector of 26 floating-point values, corresponding to 104 bytes, so the global gather/scatter messages remain firmly in the small-message regime assumed in Section 4. For the largest 3.2e10-point experiments we report results for socket pinning, giving one MPI rank per socket and 1024 ranks on 512 nodes. These socket-pinned campaigns use local depth 5 for the uniform case and 6 for the spherical case, again reflecting a choice of local tree resolution rather than a change in the weak-scaling construction itself. This choice is deliberate: for the ARCHER2 weak-scaling results reported below, socket pinning provides the most robust large-scale behavior among the node-level configurations we study because it leaves a larger shared-memory workload within each rank and fewer ranks to coordinate during global collective phases, even though finer-grained placements remain useful when strong scaling fixed problem sizes.

5.1

Weak Scaling at 4e9 Points

Figure 7 compares weak scaling up to 4e9 points on 512 ARCHER2 nodes for all three pinning strategies using the hierarchy-aware construction described above. In these experiments, runtime is reported from wall time. These experiments expose the trade-off between fine-grained MPI parallelism and reuse of sharedmemory locality inside a rank. Across both distributions we observe useful weak scaling over the full machine range, with distributed FMM runtime remaining in the low single-digit seconds even at the largest runs shown here. Prepared using sagej.cls

The International Journal of High Performance Computing Applications

XX(X) The spherical distribution exhibits visibly larger error bars and more irregular eﬀiciency curves, reflecting a combination of workload imbalance from the nonuniform point set and the run-to-run system noise already noted above, so we interpret the eﬀiciency curves in conjunction with the runtime data rather than in isolation. The dominant costs between each process pinning strategy appear to be determined by the number of points-per-MPI process, determining the depth of the local trees to traverse within each rank. The socket pinning regime appears to have the best parallel eﬀiciency. However, this is likely a reflection of the performance of our global collectives which now coordinate fewer ranks for the same problem size. However, this also indicates that this is a regime we can use to scale to significantly larger problem sizes – if we accept the larger fixed costs per-rank with each rank being responsible for a deeper FMM tree (containing more points). This cost is acceptable if the shared-memory implementation of the FMM is suﬀiciently optimized – reflected in our low overall runtimes.

5.2

Weak Scaling at 3.2e10 Points

Figure 8 shows the socket-pinned weak-scaling construction extended to the largest experiments. Each rank owns a single local tree, responsible for 32e6 points each, each local tree is refined to depth 5 in the uniform case, and depth 6 in the spherical case. We require only 21 0 ranks due to the large numbers of points handled per-rank. This figure is the clearest demonstration of the regime we ultimately target: a coarse-grained distributed decomposition in which each rank still owns enough local work for the shared-memory kernels to remain the main determinant of runtime. The large error bars for given rank configurations are hard to attribute to a single cause from our data. We have little control over the exact rank placement, or noise on the system due to co-running programs. However, just looking at the mean runtimes we observe good parallel eﬀiciency, except for the most extreme case for spherical distributions where the load imbalance due to our imposition of uniform trees becomes visible despite the eﬀiciency of our shared memory kernels. There clearly exists a trade-off between regimes where our communication strategy begins to degrade in eﬀiciency and the maximum number of particles per-rank our shared memory implementation can handle. However, even in the Prepared using sagej.cls

Kailasa

29 Weak Scaling up to 4e9 Points on 512 nodes of ARCHER2 (Uniform) Weak Scaling (Runtime)

6000

CCX CCD Socket

100

4000 CCX CCD Socket

3000 2000

Parallel Efficiency (%)

5000 Runtime (ms)

Parallel Efficiency

120

80 60 40 20

1000 10

Total Number of Points (DoFs)

Number of Ranks

Weak Scaling up to 4e9 Points on 512 nodes of ARCHER2 (Sphere)

Weak Scaling (Runtime)

3000

CCX CCD Socket

100 Parallel Efficiency (%)

Runtime (ms)

2500

Parallel Efficiency

120

CCX CCD Socket

2000 1500 1000

80 60 40 20

500 10

Total Number of Points (DoFs)

Number of Ranks

Figure 7. Weak scaling to 4e9 points on 512 ARCHER2 nodes. Top: uniform random distribution. Bottom: spherical surface distribution. Blue circles, orange squares, and green triangles correspond to CCX, CCD, and socket pinning, respectively. Times are reported from wall time, error bars are reported from standard deviation over 5 runs.

least favorable cases with up to 3.2e10 non-uniform points exposing a significant load-imbalance in work per-rank, as well as the messages sent across the network, our observed mean wall times remain in the region of 10s of seconds.

5.3

Strong Scaling

Figure 9 reports strong-scaling experiments for a fixed problem of 128 × 106 points, using up to 64 ARCHER2 nodes. In contrast to the weak-scaling study, the total number of points is held fixed here while the resources are increased in factors of two. Within each pinning strategy the number of ranks per node is fixed by the placement choice, so doubling the node count also doubles the total number of ranks and halves the points owned by each rank. The local depth is then adjusted so that the estimated leaf occupancy remains comparable across Prepared using sagej.cls

The International Journal of High Performance Computing Applications

XX(X) Weak Scaling up to 3.2e10 Points on 512 nodes of ARCHER2 (Uniform) Weak Scaling (Runtime)

Socket

100 Parallel Efficiency (%)

11000 Runtime (ms)

Parallel Efficiency

120

Socket

12000 10000 9000 8000 7000

80 60 40 20

6000 10

Total Number of Points (DoFs)

Number of Ranks

Weak Scaling up to 3.2e10 Points on 512 nodes of ARCHER2 (Sphere)

25000

Weak Scaling (Runtime)

Socket

100 Parallel Efficiency (%)

20000 Runtime (ms)

Parallel Efficiency

120

Socket

15000 10000

80 60 40 20

Total Number of Points (DoFs)

Number of Ranks

Figure 8. Weak scaling to 3.2e10 points on 512 ARCHER2 nodes using socket pinning. Top: uniform random distribution. Bottom: spherical surface distribution.

the sequence. The resulting figure should therefore be interpreted as a practical strong-scaling study for this implementation, rather than as a strictly fixeddiscretization scaling test with an identical tree at every point. In particular, the reference lines should be read only as visual guides, since the tree structure is mildly adjusted across the sequence. These runs nevertheless complement the weak-scaling results by showing how far each process placement strategy can be pushed before the amount of work per rank becomes too small to maintain good eﬀiciency. In this setting the picture is different from the largest weakscaling runs: CCX and CCD pinning remain eﬀicient to finer granularities than socket pinning, particularly for the uniform distribution, because they expose Prepared using sagej.cls

Kailasa

31 Strong Scaling For 128e6 Points on 64 Nodes of ARCHER2 (Uniform) Strong Scaling (Runtime)

Parallel Efficiency

CCX CCD Socket

100 Parallel Efficiency (%)

Runtime (ms)

CCX CCD Socket

80 60 40 20

Points Per Rank

Strong Scaling For 128e6 Points on 64 Nodes of ARCHER2 (Sphere)

Strong Scaling (Runtime)

Parallel Efficiency

CCX CCD Socket

100 Parallel Efficiency (%)

Runtime (ms)

CCX CCD Socket

80 60 40 20

Points Per Rank

Figure 9. Practical strong-scaling study for a fixed problem size of 128 × 106 points, using up to 64 ARCHER2 nodes. Top: uniform random distribution. Bottom: spherical surface distribution. Blue circles, orange squares, and green triangles correspond to CCX, CCD, and socket pinning, respectively. Dashed lines indicate ideal strong-scaling trends for the corresponding pinning strategy and should be read only as visual references. The local depth is adjusted across the sequence to maintain comparable estimated leaf occupancy, so the tree is not identical at every point.

eﬀicient process-level parallelism once the per-rank work has become small – ;likely due to the eﬀiciency of our shared-memory kernels optimized for these configurations. For the spherical distribution, eﬀiciency drops earlier for all three strategies, consistent with the greater imbalance already seen in the weak-scaling runs. Together, the weak- and strong-scaling results therefore support a consistent interpretation: for the largest production-style weak-scaling runs on ARCHER2, socket pinning is the most effective granularity because it preserves a larger shared-memory workload per rank, while finer-grained pinning is most useful when aggressively strong scaling a fixed problem. Prepared using sagej.cls

The International Journal of High Performance Computing Applications

XX(X)

5.4

Setup Breakdown

Figure 10 decomposes the one-off setup cost for the socket-pinned weak-scaling configurations at 4e9 points, while Figure 11 shows the same breakdown for the largest runs up to 3.2e10 points. In both settings, the setup is dominated by tree construction together with the initial distributed sort, so each combined plot also reports a second panel in which this dominant phase is excluded to make the remaining setup stages visible. These right-hand panels isolate the communication-specific setup work discussed in Section 3.1: the domainbounds exchanges used to establish box bounds, the global layout exchange, communicator construction, and the U-list and V-list query/setup phases.

5.5

Runtime Communication Breakdown

Figure 12 partitions the total wall time into communication and computation for the same socket-pinned weak-scaling configuration considered above for scaling up to 3.2e10 particles. Communication time is measured as wall-clock time spent in the main runtime collectives, namely MPI_Neighbor_Alltoallv, MPI_Gatherv, and MPI_Scatterv in the runtime phase of our algorithm described in Section 3.2; computation time is the corresponding wall-clock time spent in the FMM kernel phases. The stacked bars in Figure 12 report mean communication and computation times over the five repeated runs, and Figures 12 and 13 are included to illustrate how the runtime and communication decompose under this same setup. In both distributions, total wall time is dominated by computation, while communication grows with scale. This is consistent with the design choice of centralizing the global tree traversal: although global collectives are introduced, they remain a minority of the total runtime throughout the uniform study and across most of the non-uniform study. The effect is more pronounced for the spherical distribution, where irregular occupancy increases communication costs and makes communication a much larger fraction of the runtime at the largest scales. These measurements also provide the clearest experimental connection to the complexity model of Section 4. There, the favorable bound relies on two practical conditions: that neighborhood communication remains sparse, and that the global gather/scatter steps remain in the small-message regime where modern MPI implementations use eﬀicient collective algorithms. Figure 12 is broadly consistent with this picture over the measured range for the uniform distribution. This begins to break down for the non-uniform study at the largest problem Prepared using sagej.cls

Kailasa

CCX

Configuration CCD

Socket

Total Setup

Distributed Setup Excluding Trees

70000

2500

50000

Setup Time (ms)

60000 40000 30000 20000

2000 1500 1000 500

10000 0

8e6

e9 4.096

64e6

512e

Total Number of Points (DoFs)

Bounds exchange U-list query + particle

CCX

Configuration CCD

Socket

Distributed Setup Excluding Trees

3500

120000

3000

100000

2500

Setup Time (ms)

e9 4.096

Layout exchange Comm. construction

Total Setup

140000

80000 60000

2000 1500

40000

1000

20000

500

512e

Total Number of Points (DoFs) Setup phase

Trees + sorting V-list query + buffers

64e6

8e6

64e6

6 512e

Total Number of Points (DoFs)

4.096

Trees + sorting V-list query + buffers

8e6

64e6

512e

Total Number of Points (DoFs) Setup phase

Bounds exchange U-list query + particle

4.096

Layout exchange Comm. construction

Figure 10. Setup-time breakdown for the socket-pinned weak-scaling configurations on ARCHER2 at 4e9 points. Top: uniform random distribution. Bottom: spherical surface distribution. In each plot, the left panel shows the total setup cost including tree construction and the initial distributed sort, while the right panel excludes that dominant phase to reveal the remaining distributed setup stages: domain-bounds exchange, global layout exchange, communicator construction, U-list query and particle exchange, and V-list query and buffer setup.

size measured, however still doesn’t quite reach the linear trend-line. It appears likely that increasing the number of MPI ranks further would see communication complexity move to a worse regime for non-uniform point distributions under our algorithm due to the imbalance in the global collective calls. However, as we show here, we are still able to compute extremely large problems with practically useful runtimes despite this loss of optimality. Prepared using sagej.cls

The International Journal of High Performance Computing Applications

XX(X) Total Setup

Distributed Setup Excluding Trees

17500

500000

15000 Setup Time (ms)

Setup Time (ms)

400000 300000 200000

6 512e

64e6

4.096

1e6

8e9

32.76

Trees + sorting U-list query + particle

64e6

512e

4.096

Total Number of Points (DoFs) Setup phase

V-list query + buffers Bounds exchange

8e9

32.76

Layout exchange Comm. construction

Total Setup

Distributed Setup Excluding Trees

1.2

30000

1.0

25000

0.8

Setup Time (ms)

7500 2500

Total Number of Points (DoFs)

0.6 0.4 0.2 0.0

10000 5000

100000 0

12500

20000 15000 10000 5000

64e6

512e

4.096

Total Number of Points (DoFs)

8e9

32.76

Trees + sorting V-list query + buffers

64e6

512e

4.096

Total Number of Points (DoFs) Setup phase

U-list query + particle Bounds exchange

8e9

32.76

Layout exchange Comm. construction

Figure 11. Setup-time breakdown for the largest socket-pinned weak-scaling configurations on ARCHER2, up to 3.2 × 1010 points. Top: uniform random distribution. Bottom: spherical surface distribution. As in Figure 10, the left panel shows the total setup cost including tree construction and the initial distributed sort, while the right panel excludes that dominant phase to expose the remaining distributed setup stages.

Discussion

Our method deliberately targets uniform octrees. This design choice simplifies the construction of static communication graphs and allows for fully precomputed interaction lists and data exchanges. As a result, we avoid runtime filtering, memory allocation, and communication planning. In doing so, we sacrifice adaptability to nonuniform trees in favor of maximal precomputation and memory reuse. While this limits generality in comparison to (Lashuk et al. 2009; Abduljabbar et al. 2017), it enables the full exploitation of shared-memory hardware on uniform problems – and moderately irregular point distributions Prepared using sagej.cls

Wall Time (ms)

Kailasa

~ log(p) scaling ~ O(p) scaling Communication Time Computation Time

128

Number of MPI Ranks

1024 ~ log(p) scaling ~ O(p) scaling Communication Time Computation Time

128

Number of MPI Ranks

1024

Figure 12. Aggregate communication and computation time for the same socket-pinned weak-scaling configurations used in the preceding scaling study. Top: uniform random distribution. Bottom: spherical surface distribution. The stacked bars report mean times over five repeated runs and separate communication from computation, while the dashed red and green curves indicate reference log P and P trends computed from the ranks in each figure and normalized with respect to the first measured value.

are also handled well by shared-memory optimizations for the ULI step – which is typically the most straightforward phase of the FMM to parallelize. Prepared using sagej.cls

The International Journal of High Performance Computing Applications

XX(X)

ScatterV GatherV Neighbor AllToAllV

350

Wall Time (ms)

300 250 200 150 100 50 0

128

Number of MPI Ranks

1024

1400

ScatterV GatherV Neighbor AllToAllV

Wall Time (ms)

1200 1000 800 600 400 200 0

128

Number of MPI Ranks

1024

Figure 13. Breakdown of communication time for the same socket-pinned weak-scaling configurations used above. Top: uniform random distribution. Bottom: spherical surface distribution. The stacked bars report mean communication times over five repeated runs and are plotted for illustration; from bottom to top they show MPI_Scatterv, MPI_Gatherv, and MPI_Neighbor_Alltoallv.

We also introduce a controlled bottleneck by executing the global portion of the FMM in a single shared-memory execution context. This is a conscious tradeoff: by centralizing global tree traversal, we preserve the performance Prepared using sagej.cls

Kailasa

benefits of shared-memory kernels, particularly for the VLI step. As shown in the ARCHER2 experiments of Section 5, this choice does not appear to hinder scalability over the ranges studied: we observe useful weak scaling together with useful practical strong-scaling behavior on smaller configurations. Our approach may lose optimal worst-case complexity with P depending on the point distribution, with a clear trade-off between the number of points handled per-rank (determined by the eﬀiciency of our shared memory implementation) and the total number of ranks. However, the clear strength of our approach is the ease with which it can be implemented in FMM software already developed and optimized for a the single-node case, using only standard collective operations.

Conclusion

We have presented a communication scheme for distributed FMMs based on standard MPI neighborhood collectives, uniform trees, and precomputed communication layouts. The central goal of the method is to simplify the distributed-memory extension of existing high-performance FMM software while preserving the shared-memory kernels that dominate runtime on modern machines. The ARCHER2 results show that this design delivers weak scaling for uniform and moderately non-uniform point distributions, together with useful practical strong-scaling behavior on smaller configurations. The distributed FMM runtime remains on the order of tens of seconds even at the largest scales studied. Although favorable runtime communication complexity can still be obtained in the best case, the main contribution of this work is to show that a distributed FMM can be realized using standard MPI collectives and precomputed communication layouts while still delivering good large-scale runtime performance. Future work would involve deploying our software on other leading HPC systems, to identify whether the issues we faced with system noise were unique to ARCHER2’s communication subsystem, and identify the with more granularity the reasons for the departure from optimal communication complexity in the case of non-uniform point distributions at large P . Prepared using sagej.cls

The International Journal of High Performance Computing Applications

XX(X) References Abduljabbar M, Markomanolis GS, Ibeid H, Yokota R and Keyes D (2017) Communication reducing algorithms for distributed hierarchical n-body problems with boundary distributions. In: International Conference on High Performance Computing. Springer, pp. 79–96. AbdulJabbar M, Yokota R and Keyes D (2014) Asynchronous execution of the fast multipole method using charm++. arXiv preprint arXiv:1405.7487 . Anderson CR (1992) An implementation of the fast multipole method without multipoles. SIAM Journal on Scientific and Statistical Computing 13(4): 923–947. Blanchard P, Bramas B, Coulaud O, Darve E, Dupuy L, Etcheverry A and Sylvand G (2015) Scalfmm: A generic parallel fast multipole library. In: SIAM Conference on Computational Science and Engineering (SIAM CSE 2015). Blelloch GE, Leiserson CE, Maggs BM, Plaxton CG, Smith SJ and Zagha M (1991) A comparison of sorting algorithms for the connection machine cm-2. In: Proceedings of the third annual ACM symposium on Parallel algorithms and architectures. pp. 3–16. Bramas B (2020) Tbfmm: A c++ generic and parallel fast multipole method library. Journal of Open Source Software 5(56): 2444. Chandramowlishwaran A, Williams S, Oliker L, Lashuk I, Biros G and Vuduc R (2010) Optimizing and tuning the fast multipole method for state-of-the-art multicore architectures. In: 2010 IEEE International Symposium on Parallel & Distributed Processing (IPDPS). IEEE, pp. 1–12. Darve E (2000) The fast multipole method: numerical implementation. Journal of Computational Physics 160(1): 195–240. Fong W and Darve E (2009) The black-box fast multipole method.

Journal of

Computational Physics 228(23): 8712–8725. Greengard L and Rokhlin V (1987) A fast algorithm for particle simulations. Journal of computational physics 73(2): 325–348. Gumerov NA and Duraiswami R (2008) Fast multipole methods on graphics processors. Journal of Computational Physics 227(18): 8290–8313. Hackbusch W (1999) A sparse matrix arithmetic based on-matrices. part i: Introduction to-matrices. Computing 62(2): 89–108. Hackbusch W and Nowak ZP (1989) On the fast matrix multiplication in the boundary element method by panel clustering. Numerische Mathematik 54(4): 463–491.

Prepared using sagej.cls

Kailasa

Ibeid H, Yokota R and Keyes D (2016) A performance model for the communication in fast multipole methods on high-performance computing platforms.

The

International Journal of High Performance Computing Applications 30(4): 423– 437. Kailasa S (2025a) kifmm-rs: A kernel-independent fast multipole framework in rust. Journal of Open Source Software 10(110): 7124. Kailasa S (2025b) Modern Resesarch Software For Fast Multipole Methods. PhD Thesis, UCL (University College London). Kailasa S, Betcke T and Kazdadi SE (2024) M2l translation operators for kernel independent fast multipole methods on modern architectures.

arXiv preprint

arXiv:2408.07436 . Kang Q, Ross R, Latham R, Lee S, Agrawal A, Choudhary A and Liao Wk (2020) Improving all-to-many personalized communication in two-phase i/o. In: SC20: International Conference for High Performance Computing, Networking, Storage and Analysis. IEEE, pp. 1–13. Kim J, Dally WJ, Scott S and Abts D (2008) Technology-driven, highly-scalable dragonfly topology. ACM SIGARCH Computer Architecture News 36(3): 77–88. Lashuk I, Chandramowlishwaran A, Langston H, Nguyen TA, Sampath R, Shringarpure A, Vuduc R, Ying L, Zorin D and Biros G (2009) A massively parallel adaptive fast-multipole method on heterogeneous architectures. In: Proceedings of the Conference on High Performance Computing Networking, Storage and Analysis. pp. 1–12. Ltaief H and Yokota R (2014) Data-driven execution of fast multipole methods. Concurrency and Computation: Practice and Experience 26(11): 1935–1946. Malhotra D and Biros G (2015) Pvfmm: A parallel kernel independent fmm for particle and volume potentials. Communications in Computational Physics 18(3): 808–830. Martinsson PG and Rokhlin V (2007) An accelerated kernel-independent fast multipole method in one dimension. SIAM Journal on Scientific Computing 29(3): 1160– 1178. Messner M, Bramas B, Coulaud O and Darve E (2012) Optimized m2l kernels for the chebyshev interpolation based fast multipole method.

arXiv preprint

arXiv:1210.7292 . Netterville N, Fan K, Kumar S and Gilray T (2022) A Visual Guide to MPI All-to-all.

Proceedings - 2022 IEEE 29th International Conference on High

Performance Computing, Data, and Analytics Workshops, HiPCW 2022 : 20– 2710.1177/ToBeAssigned10.1109/HiPCW57629.2022.00008.

Prepared using sagej.cls

The International Journal of High Performance Computing Applications

XX(X) Sundar H, Malhotra D and Biros G (2013) Hyksort: a new variant of hypercube quicksort on distributed memory architectures.

In: Proceedings of the 27th

international ACM conference on international conference on supercomputing. pp. 293–302. Sundar H, Sampath RS and Biros G (2008) Bottom-up construction and 2: 1 balance refinement of linear octrees in parallel. SIAM Journal on Scientific Computing 30(5): 2675–2708. Takahashi T, Cecka C, Fong W and Darve E (2012) Optimizing the multipole-to-local operator in the fast multipole method for graphical processing units. International Journal for Numerical Methods in Engineering 89(1): 105–133. Thakur R, Rabenseifner R and Gropp W (2005) Optimization of collective communication operations in MPICH. International Journal of High Performance Computing Applications 19(1): 49–66. 10.1177/ToBeAssigned10.1177/1094342005051521. Wang T, Yokota R and Barba LA (2021) Exafmm: a high-performance fast multipole method library with c++ and python interfaces. Journal of Open Source Software 6(61): 3145. Warren MS and Salmon JK (1993) A parallel hashed oct-tree n-body algorithm. In: Proceedings of the 1993 ACM/IEEE conference on Supercomputing. pp. 12–21. Williams S, Waterman A and Patterson D (2009) Roofline: an insightful visual performance model for multicore architectures. Communications of the ACM 52(4): 65–76. Ying L, Biros G and Zorin D (2004) A kernel-independent adaptive fast multipole algorithm in two and three dimensions. Journal of Computational Physics 196(2): 591–626. Yokota R (2013) An fmm based on dual tree traversal for many-core architectures. Journal of Algorithms & Computational Technology 7(3): 301–324. Yokota R and Barba LA (2012) A tuned and scalable fast multipole method as a preeminent algorithm for exascale systems. The International Journal of High Performance Computing Applications 26(4): 337–346. Yokota R, Turkiyyah G and Keyes D (2014) Communication complexity of the fast multipole method and its algebraic variants. arXiv preprint arXiv:1406.1974 .

Prepared using sagej.cls

Record · ID 124026 · SHA-256 93fc16b16f9e731f
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.