The evaluation of pairwise interactions between N particles in R d , d ∈ Z + , commonly referred to as the N-body problem, is a classic problem in computational science with a wide variety of applications. Generally, the problem may be defined as the evaluation of the sum Φ(r i ) = N ∑ j=1 ψ(r i , r j )σ j , i = 1, . . . , N (1) where ψ(r, r ′ ) : R d × R d → C is called the kernel function, Φ(r) is the potential at the point r, and σ j is the coefficient of the jth source particle located at r j ∈ R d . Applications of this equation range from mathematical physics wherein the kernel is a Green's function used to compute fields from sources or modern machine learning/interpolation methods. Example kernel functions that are often used are given in Table 1. Several machine learning techniques, increasingly popular in recent years, and interpolation methods [1][2][3][4][5] also rely on the N-body calculation. As is well known, the principal difficulty with the N-body calculation is that it is an inherently dense computation, resulting in O(N 2 ) computational complexity for evaluation of (1). This operation may also be interpreted as a matrix-vector product Ax = b, with the kernel matrix given by A ij = ψ(r i , r j ) and vectors x j = σ j , b i = Φ(r i ). As N becomes very large the cost of the matrixvector product becomes prohibitive, and fast, error-controllable evaluation methods become necessary.
Given the vast applicability of the N-body calculation, a number of researchers have contributed to a robust and varied literature on fast algorithms for this problem. Perhaps the most impactful of these are the fast multipole methods (FMMs), first introduced by Greengard and Rokhlin for the Laplace potential [6,7]. In these methods, elements of the kernel matrix are not explicitly formed; instead, the FMM is a method for rapidly computing the application of the kernel matrix to a vector. One should also mention the Barnes-Hut (or treecode) algorithms [8], upon which the fundamental ideas of the FMM were based. For a discussion of modern FMMs, it is useful to define a coarse taxonomy of methods. Kernel functions of interest typically fall into one of two categories: oscillatory (kernels which encode phase information, e.g., Helmholtz potential) and non-oscillatory (kernels which do not encode phase, e.g., Laplace potential). Designing an FMM for oscillatory kernels requires particular care which will not be discussed here; the interested reader is instead referred to Refs. [9][10][11][12][13]. FMMs may additionally be categorized as either kernel-dependent FMMs, which are based on a kernel-dependent factorization of a particular kernel, or kernel-independent FMMs, which are not. When applicable, kernel-dependent FMMs may be considerably more efficient than their kernel-independent counterparts thanks to the optimal exploitation of particular properties of the factorization; see, for instance, the Ref. [9]. However, implementation of kernel-dependent FMMs is often considerably more involved than that of kernel-independent FMMs. Within the class of kernel-independent FMMs, there are several options, each with their own strengths and weaknesses. As each of these methods can in some sense be viewed as a local polynomial approximation, we shall denote the order of approximation as P. The KIFMM algorithm [14] offers acceleration with low O(NP 2 ) computational complexity thanks to a framework of surface integrals radiating equivalent source densities, requiring only a surface discretization of each cell. However, the use of the integral radiation operator restricts the applicability of this algorithm to kernels satisfying a radiation condition, e.g. O(r -ν ), ν ≥ 1. The black-box FMM (bbFMM) [15] employs tensor-product Chebyshev interpolation, requiring a volumetric discretization of each cell which implies an O(NP 3 ) cost due to dense matrix operations between cells with (P + 1) 3 interpolation nodes each. However, implementation is straightforward and the algorithm can be used with any smooth kernel. Additionally, the cost can be reduced via low-rank factorization or FFT acceleration [11,16].
Another method developed in 2007, is the accelerated Cartesian expansion (ACE) method [17] that uses the Taylor series expansion in Cartesian tensor form to systematically express addition theorems for any arbitrary smooth kernel. More importantly, ACE does not require kernel functions to satisfy the radiation condition. The structure of the method offers numerous computational advantages elucidated and demonstrated in Refs. [17][18][19][20]. Further, while its polynomial representation implies a cost of O(NP 3 ), further analysis reveals the constant is small and the number of degrees of freedom per cell is asymptotically about 1/6 of that of the bbFMM, rendering the algorithm highly efficient for low P. In some sense, ACE straddles the boundary between kernel-dependent and kernel-independent: While the KIFMM and bbFMM require only a function handle to evaluate the kernel, ACE requires knowledge of the kernel's derivatives. Fortunately, for many of the smooth kernels for which FMM-like methods are desired, these are often available analytically.
Parallelization of the aforementioned algorithms has also been the subject of intense focus over recent years. The advent of powerful supercomputers with massive numbers of CPU cores has driven interest in extremely large N-body calculations involving millions to billions of particles. Special-purpose algorithms for the Laplace kernel have been perhaps the most well-studied in this regard [21,22], but several parallel algorithms for kernelindependent FMMs have also been put forth and have exhibited similar scaling properties [14,23,24]. The sheer size of the problems of interest and variations in the density of the underlying particle distribution present challenges in parallelization.
This paper presents a scalable parallel algorithm for evaluating N-body sums involving arbitrary non-oscillatory potentials, including those that do not decay as r → ∞, via the ACE algorithm.
Building upon [18], our contributions are:
• Algorithms for efficiently constructing and load-balancing highly non-uniform trees;
• Controllably-accurate ACE operators for non-uniform trees;
• Optimizations to the matrix-vector product;
• Performance, accuracy, and scalability for non-decaying kernels over large and highly non-uniform distributions;
• Direct performance comparison of our method with a highly -optimized KIFMM implementation.
The paper is organized as follows. Section 2 gives a brief overview of the ACE algorithm and establishes notation. Section 3 describes the new parallel algorithm, and Section 4 gives numerical and performance results demonstrating the accuracy and efficiency of the presented approach, and a comparison of its performance against the KIFMM. Appendix A provides insight how efficiencies can be realized in applying ACE to evaluate Stokes potentials. Finally, Appendix B provides the theorems necessary to completely define ACE (for completeness purposes).
In this section, we give a succinct summary of the relevant aspects of the ACE algorithm and highlight some useful features. For more details, the reader is referred to Ref. [17]. Consider the N-body problem described in (1). Let Ω ∈ R 3 denote an axis-aligned cubical domain containing all of the particles. Like all FMM-style algorithms, ACE relies on a hierarchical octree decomposition of the space Ω, and the number of levels in the tree will be denoted by L, where level 1 represents the root of the tree and level L is the leaf level, or the finest level of refinement.
Consider two boxes Ω s and Ω o that contain source and observers, respectively, in the oct-tree decomposition of Ω, and let the parents of Ω s and Ω o be denoted using Ω p s and Ω p o . Centers of these boxes are denoted using r c s , r c o , r c s,p , and r c o,p . To set the stage, we start with some remarks and notations; (i) an nth rank totally symmetric tensor A (n) contains (n + 1)(n + 2)/2 independent components as opposed to 3 n components; (ii) in this compressed form this tensor can be represented using A (n) (n 1 , n 2 , n 3 ) where n = n 1 + n 2 + n 3 ; (iii) an example of such a tensor is the polyadic associated with r which is given by rr • • • rr n times = r n and can be represented in compressed form as r (n) = x n 1 y n 2 z n 3 ; (iv) an m-fold contraction between two tensors A (n) and B (m) is denoted using m) ; and (v) a direct product between two tensors can be written as
The heart of ACE is the Taylor series expansion, which provides a natural framework for developing addition theorems. Consider a source located at a point r ′ observed at a point r.
for |X| > |d|, can express an addition theorem for any sufficiently smooth, non-oscillatory kernel function. In practice, this sum is truncated at some finite P, with increasing accuracy for larger values of P. From this representation, the N-body sum (1) may be recast as
where Φ N (r i ) represents contributions to Φ(r i ) from particles in octree boxes adjacent to that containing the ith observation point r i . The remaining (far) contributions, given by Φ F (r i ), are calculated using the addition theorem.
The ACE algorithm follows the standard fast multipole procedure; for particulars, the reader is referred to [17]. First, multipole expansions of sources are computed for each leaf-level box and aggregated up the octree. Next, local expansions are formed for each octree box by accumulating translated multipole expansions in its interaction list. As implied by ( 2), a local expansion may be interpreted as the coefficients of the Taylor polynomial representing the influence of all particles exterior to the near-field of its octree box. Finally, local expansions are disaggregated in a topdown fashion to the leaf level, where they are used to compute the potential at each observer particle. For non-uniform particle distributions it is advantageous to use an adaptive octree which can conform to the local distribution density (see Section 3). Two additional steps involving interactions between octree boxes of different sizes are required [13,25]. Because the condition |d| < |X| is not satisfied in three dimensions for these interactions, the operators involved take the form of either (i) directly evaluating multipole expansions at observer particles or (ii) summing the influence of individual source particles onto local expansions; they are trivially derived by subsuming either the observer or source parts, respectively, of d into the vector X, yielding a valid addition theorem.
We now state several interesting and/or useful properties of the algorithm. First, it has been shown that the operators involved in aggregation and disaggregation are exact, i.e. they incur no numerical error beyond the specification of P and these are independent of the height of the tree. Second, the source-side multipole expansions do not involve the kernel function; hence, multiple potentials can be computed at once from the same source data, provided one can store a set of local expansions for each kernel considered. Third, operators for going up and down the tree (between any two consecutive levels) differ only by a constant, and are easily precomputed. The same is true for certain kernel functions. Fourth, the k-fold gradient of the potential is trivially evaluated using the formula [17] ∇
As we will demonstrate, this property is especially useful for evaluating multiple distinct potentials which derive from a common kernel function, or in some cases reduces the overall cost when a potential can be expressed in terms of differential operators acting upon a simple kernel function (like the Stokes potential). Pertinent details for application to the Stokes potential are presented in Appendix A.
In this section, we outline our parallel algorithms for constructing, adapting, and load balancing the tree, as well as parallel potential evaluation. While parallelization of FMMs is well-documented in the literature [3,23,26], as we will demonstrate in Section 4, redundant M2L computations on processes sharing the same local expansions in the commonly used locally essential tree (LET) based algorithms would significantly hamper the performance of our proposed framework, particularly for high-accuracy calculations. The M2L stage incurs the largest computational cost in ACE in terms of P, whereas in methods such as KIFMM this is not the case [23]. For this reason, our parallel evaluation strategy differs from the usual LET based implementations, and therefore we present our parallel implementation in detail.
Construction of the distributed octree starts by partitioning the input particles equally among all processes. On each process, the equi-distributed particles are first hashed into Morton keys, which is a binary encoding of the leaf boxes according to their space-filling Morton-Z curve ordering. A parallel bucket sort is then used to assign each process a distinct set of leaves contiguous in the Morton-Z ordering such that each process gets a roughly equal number of particles. This initial partitioning is performed at the finest refinement level, i.e., at level L, in preparation for the ensuing adaptive tree coarsening and load balancing stages, and is concluded with a local post-order traversal tree construction on each process from leaves all the way up to the root.
The above partitioning scheme necessarily incurs duplicate copies of internal tree nodes, corresponding to common ancestors of leaf nodes residing on different processes. Such nodes are referred to as plural nodes. Treatment of computations associated with plural nodes is one aspect of our parallelization scheme that differentiates it from the LET based implementations. Hence, we provide some terminology to aid the discussion. For any plural node, the process with the highest rank which owns a copy is designated as the resident process, and as such is deemed responsible for all its tree-based interactions. Other processes with a copy of the plural node are called users of the node, and are directed only to send and receive data to and from the resident process. A plural node is referred to as a shared node on its resident process, and the users' copies are called duplicate nodes. As a consequence of the post-order traversal tree construction, it can be shown that all duplicate nodes in a process' local sub-tree appear consecutively at the end of the post-order traversal sequence. Shared nodes may appear anywhere within the post-order sequence, though they tend to appear toward the beginning. As we will discuss, these properties allow effective overlapping of communication and computation in the potential evaluation stage.
Random (or homogeneous) distributions of source or observer particles are well-represented by a uniform tree structure. Within such distributions, the number of particles per box is approximately constant across the entire simulation domain, facilitating the linear scaling of the ACE algorithm or any other FMM implementation for that matter. For non-uniform distributions though, i.e., those with sub-regions of relatively high particle density, the minimum box size must be decreased to prevent near-field costs from dominating the computation. Hence, in a uniform tree framework, the densest region of discretization dictates the leaf box size throughout the entire tree. The result is large swaths of space with a poor work-to-particle ratio in sparsely populated regions, degrading the cost scaling of parallel ACE computations.
To prevent such inefficiencies, we adapt the tree structure to the particle distribution so as to restore the approximate uniformity in leaf box population throughout the tree. Adaptation of trees for N-body simulation has been well-studied and has become a standard feature of modern tree-based simulation methods for non-oscillatory kernels [24,[27][28][29]. We employ a bottom-up scheme for both constructing and merging the tree to minimize communication costs. Our parallel algorithm for merging the distributed octree is given in Algorithm 1. Simply put, this algorithm starts with a uniform tree constructed according to the description in Section 3.1 and merges each set of sibling leaf boxes into their parents, if the total number of particles in these leaf boxes do not exceed a pre-determined particle threshold M. This scheme ensures that leaf boxes in more sparsely-populated regions contain roughly as many points as those in the dense regions, thereby avoiding excessive tree interactions (in sparse regions) or expensive near field computations (in dense regions).
In adaptive octree algorithms, it is necessary to evaluate interactions between boxes of different sizes. Fig. 2(a) gives an illustrative example. As usual, these interactions are classified into U, V , W , X -lists based on the adjacency and relative size of source and observer box pairs [7,25]. For the sake of completeness, these lists are defined as follows:
• U-list: source and observer boxes are adjacent and are both leaves, i.e., they are within the near field of each other;
• V-list: source and observer boxes are not adjacent but their parents are, i.e., they are inside the far field of each other;
• X-list: source and observer boxes are not adjacent, but the observer box is (i) adjacent to the parent of the source box, (ii) a leaf, and (iii) larger than the source box, corresponding to a cross-level interaction with an observer node higher up in the octree;
• W-list: reciprocal interactions of those in the X -list, corresponding to a cross-level interaction with an observer node lower down in the octree.
For highly non-uniform distributions, unconstrained merging of the tree results in a tree with a matching degree of nonuniformity. While this is normally desirable from the perspective of reducing computational costs (less tree nodes means less number of interactions overall), it has a downside in terms of memory utilization. In a uniform tree where no X-list or W-list interactions exist, there can be at most 6 3 (total number of boxes in the neighborhood of a box) -3 3 (number of near field boxes) = 189 boxes in the far field of a box. Leveraging symmetries and scale invariance of certain kernels can significantly reduce the number of translation operators needed for V -list interactions. However, in regions containing sharp discontinuities in leaf box size, the number of X -and W -list interactions can become quite large, and the number of unique translation operators needed to perform the corresponding X -and W -list interactions can also be very large. The 2:1 balance constraint [24,29] remedies this storage problem by disallowing adjacent leaf boxes to differ in size by more than a factor of two, significantly reducing the number of different X and W interactions and the memory overhead for storing their particle-specific translation operators. Algorithm 1 for merging the tree includes the option to impose this constraint.
After the tree merging procedure under the 2:1 balance constraints is complete, each process is left with a subset of the original uniform tree. As such, the computational profile is sufficiently different from the original tree that re-balancing the load on processes is necessary. To accomplish this, we follow an empirical load balancing strategy where we estimate the work attendant to each node in the distributed tree and aim to assign them in a load-balanced manner.
The cost per leaf is determined as follows. As part of the initialization operations, we first time a set of dummy operations to obtain cost estimates for each of the U, V , W , and X -list interactions. The cost for each node in the tree is then determined by multiplying the number of each interaction by its corresponding cost estimate, and summing these costs across all interaction types. As with the initial partitioning, we partition the distributed tree by determining a set of separators between contiguous chunks of leaf boxes. Therefore we account for the costs of the interior nodes by percolating their estimated costs down to the leaf boxes and adding them to the extant cost estimate at each leaf. A variant of Algorithm 1 from [29] is then used to determine P -1 locations at which to split the leaf level Morton curve for a re-balancing of the computational load. While the overall strategy is similar to the costzones approach of [27,28], our strategy accounts for the work of the entire tree instead of just the leaves. This is important for surface geometries with non-uniform distributions, as the amount of work above the leaves is not negligible.
The parallel potential evaluation is performed in three stages: The upward pass (M2M), translation (M2L), and the downward pass (L2L). Algorithm 2 describes the M2M (upward pass) stage of our implementation, which essentially entails shifting the multipole data of each tree node to the center of its parent box and aggregating the shifted multipole data from all siblings. Each process starts processing the nodes in its local subtree from rightto-left in the post-order traversal sequence. By choosing to go from right-to-left, we ensure that duplicate nodes are encountered toward the beginning of the M2M stage, and shared nodes The translation stage commences on each process once the local upward pass is complete. This stage includes three substages for computations of X , V and W lists. First, X -list interactions are handled by sending the required multipole expansions to processes with X -list observers, where the observed potentials are computed. Next, V -list interactions are carried out. Typically, this is the most expensive stage of any fast multipole-like method. To hide communication overheads, the source multipole expansions to be exchanged are divided up into packets and communicated to processes that need this information. As each source expansion is received in full, all its V -list interactions are computed, making good use of temporal locality. Once all source expansions in the current packet are exhausted of work, the next packet is constructed and communicated using nonblocking primitives to facilitate overlapping of communication and computation. This process continues until all remote V -list interactions are completed. Local V -list interactions, i.e., those in which both the source and observer boxes belong to the same process, are computed next. Finally, W -list interactions are computed by exchanging the source weights, i.e. u i in (B.1), and applying the appropriate translation operators to form the local expansions of these sources.
To complete the evaluation, Algorithm 3 describes the L2L (downward pass) stage which involves re-centering and adding parents' local expansions to their childrens'. Contrary to the M2M stage, here we employ a pre-order traversal of the local subtree so that the shared nodes are now encountered at the beginning and the duplicates are encountered toward the end of the L2L stage. This way, by using MPI's non-blocking IBcast primitive, we ensure that broadcast for a shared node can be completed in the background, while L2L computations of other nodes are being performed. Also, the result of the V -list interaction evaluations for duplicate nodes is communicated to users during the downward pass, again using non-blocking primitives for overlapping communications with computations.
As described above, updating multipole and local expansions for nodes that are shared between processes are performed efficiently in our implementation using asynchronous communications. We note that at any level, the local subtree of any process can have at most two plural nodes -one shared, one duplicate. It follows, then, that we may create a total of 2L MPI communicators, each with groupings of shared nodes and their users, so that we can take advantage of non-blocking variants of MPI's reduce and broadcast operations. We note that our evaluation algorithm differs from that of the common LET-based implementations in that we avoid duplicating V -list interactions, which are typically the most computationally expensive part of any fast multipole-like algorithm (see Section 4). Each V -list interaction is computed exactly once, and the resident process is responsible for such computations.
In this section, we present an array of results demonstrating error convergence and performance of the presented algorithm and its parallel implementation. end if 20: end for
We first examine the error convergence of the ACE algorithm applied to the r -ν , Yukawa, and Stokes potentials (the Stokeslet, rotlet, and stresslet; see the Appendix for further details). The first two are computed using scalar sources, while the Stokes potentials are calculated using vector-valued sources. The error for Stokes potentials is compared to the analytical by taking an inner product of the tensor-valued potential with each observer's polarization vector. We first evaluate the r -ν and Yukawa potentials for 3.125 million points within a 0.64 m cube. The distribution is mapped onto a tree with a leaf box diameter of d 0 = 0.01 m, yielding a 7-level tree with 12 points per box on average. A single buffer box is used for the far-field. Fig. 3 shows convergence in the L 2 error for the far-field with increasing expansion order P = 1, 3, . . . , 19 for both kernels. Different parameters are used to control the growth or decay of the kernels as r increases.
We next consider the evaluation of error in Stokes potentials for a collection of points within a cube with both uniform and non-uniform spatial distributions. Both distributions are characterized by 3.125 million points within a 0.64 m cube. In the case of the uniform distribution, we employ a uniform tree with increasing box size as P is increased to minimize runtime. The nonuniform distribution is generated by generating random points within the unit cube and raising the generated x, y, z positions to the powers of 1.2, 0.7, and 1.7, respectively, before rescaling the results to fit into the 0.64 m cube. In this case, we use a non-uniform tree with 2:1 balance and increase the minimum box size with P to approximately minimize runtime according to the uniform-tree cost estimate. The error convergence with increasing P for both distributions is shown in Fig. 4. The error metric used here relies on randomly sampling the analytical fields at one observer on each process, ensuring a good spatial distribution of observers, and comparing it with the computed potential using ACE. For each selected particle, the error in the potential relative to the analytic solution is computed, and the average of these errors is reported.
Finally, we examine the performance of the algorithm in evaluating all three Stokes potentials simultaneously from the kernel ψ(r) = r using the formulas given in Appendix A. We empirically selected the box sizes for best performance. These results were obtained on the Haswell partition of the Cori supercomputer at the National Energy Research Scientific Computing Center (NERSC). This cluster comprises 2388 compute nodes with two sockets each, populated by 16-core Intel Xeon E5-2698 v3 ''Haswell'' CPUs running at 2.3 GHz and 64 GB DDR4 RAM at 2133 MHz per socket. The algorithm was implemented in Fortran 90 using double-precision arithmetic and parallelized strictly in distributed-memory fashion using MPI. The code was compiled using the Intel compiler version 18 with optimization and architecture-specific instructions using the -O3 -xHost flags. For these runs we use vector-valued sources in R 3 and evaluate the far-field only, selecting P = 7 to give O(10 -5 ) accuracy for the Stokeslet and O(10 -3 ) accuracy for the rotlet and stresslet. We next consider a uniform random distribution of 5 billion randomly-oriented vector sources with unit norm inside a cube of diameter 8 m. The minimum box size was set to 0.015 m, yielding an 11-level uniform tree with 33 particles per leaf box on average. Fig. 5(a) shows the parallel efficiency (strong scaling) of this distribution up to 16,384 processes with respect to 1024 processes. The increase in efficiency from 2048 to 4096 processes is due to the fact that the distribution is randomly re-generated for each run, given the extremely large number of particles. The efficiency is over 90% for all cases. Fig. 5(b) shows the timings for each process for the 16,384 process case. The C2M stage requires on average 0.54 s while the L2O stage, modified for the Stokes potentials, requires 6.05 s, about an 11.2X increase. Most of this increase is due to the calculation of the stresslet potential. The overall computational time is dominated by the V -list calculation, taking over 45 s of the 54 s, underscoring the importance of avoiding redundant V-list calculations related to duplicate nodes at upper levels of the tree.
We now consider uniformly and non-uniformly distributed sources on the surface of a sphere. Both distributions comprise 1.024 billion particles. For the uniform distribution, the sphere diameter is 2 m. A uniform tree is used with a leaf box size of 5 × 10 -4 m, yielding a 13-level tree with leaves containing 15 particles per box on average. The strong scaling and perprocess timings on 16,384 processes for the uniform distribution are shown in Fig. 6. A mock distribution is also shown inset in the strong scaling plot. In this case, the efficiency is over 77%. The choppiness of the L2L stage timings is due to the fact that the number of particles assigned to each process is uneven; instead, the load balancing algorithm aims to balance the overall computational load. As the bulk of the work lies in the V -list evaluation, Fig. 6 suggests the algorithm does a reasonably good job.
In the case of the non-uniform distribution, the sphere has a diameter of 8 m, and the particles are clustered around the north and south (±z) poles (as depicted inset in Fig. 7(a)). The minimum box size is chosen as 6.25 × 10 -7 m in diameter, resulting in a 25-level tree. We set s max = 50 particles per box, resulting in about 18 particles per box on average and distributing leaves over the bottom 16 levels of the tree. Despite the extremely non-uniform distribution of particles, the strong scaling shown in Fig. 7(a) is still as high as 78% efficient on 16,384 processes. In addition, the per-process timings presented in Fig. 7(b) for the far-field evaluation suggest that the load balancing algorithm does its job reasonably well. Again, the choppiness here is due to the fact that particles are not evenly distributed across processes, and processes with high particle counts spend considerable time in the L2O stage. Those with the most particles stick out in these plots.
In this section, we compare the performance of our parallel algorithm with that of the point-based kernel-independent FMM [14,23] implemented in the open-source PVFMM package [24]. We selected this code as a benchmark because it is perhaps the most well-known code of its kind. Its authors have extensively optimized the algorithm and implementation, and they have compared it with other contemporary codes [24,30,31]. The runs in this section were performed on a single compute node with an Intel (R) Xeon (R) Gold 6148 CPU at 2.40 GHz with 83 GB RAM. Because our implementation does not make use of multithreading, we considered strictly MPI parallelism. We did not do any detailed manual optimization of our implementation ψ(|r|) = e -γ r r with γ = 1. We ran both codes at two target accuracy settings (in the sense of the L 2 relative error), first at the lower accuracy of ε = 10 -4 and then at the higher accuracy of ε = 5 × 10 -7 . In each case, the number of points per leaf box was tuned to the best of our ability to optimize the execution time for both codes. Table 2 summarizes the simulation parameters for all cases. The KIFMM parameters m and s represent the multipole expansion order and the maximum number of particles per leaf box, respectively. The results of our performance comparison show that the algorithm of the present paper is competitive with the KIFMM algorithm for the geometry and target accuracies considered. Fig. 8(a) shows the time taken for a matrix-vector product in each code for both accuracies. For the lower-accuracy case, our algorithm is faster for all process counts, though the KIFMM code is faster for the higher accuracy case. This is expected, as the complexity estimates for ACE and KIFMM when the number of points per box is optimized for runtime are O(NP 3 ) and O(Nm), respectively [14,17]. As shown in Fig. 8(b), our algorithm scales reasonably well up to 20 MPI ranks, but not quite as well as the KIFMM code. We note that the KIFMM code ran out of memory for 20 MPI ranks in high-accuracy mode.
In this paper, we have presented a fast, adaptive, and scalable algorithm for evaluating pair potentials involving non-oscillatory kernels of arbitrary form. We introduced controllably-accurate operators for calculating interactions arising from the nonuniform tree structure and algorithms for building, load balancing, and efficiently traversing the resulting non-uniform tree structure. We demonstrated the merits of this algorithm in terms of both error control and parallel performance through a series of numerical example involving several different kernels, including non-decaying kernels, and both uniform and extremely nonuniform distributions of up to 5 billion particles on as many as 16,384 CPU cores. The presented methods have many potential applications in large-scale computational physics, machine learning, and beyond.
where M (n) (r c s ) is the rank-n multipole tensor about r c s .
Next, these multipoles can be re-expressed about r c s,p using These theorems collectively provide a framework for expressing addition theorems for any smooth and non-oscillatory kernel