A GPU-Accelerated Fast Summation Method Based on Barycentric Lagrange Interpolation and Dual Tree Traversal

Tree-based fast summation methods such as the DTT, TC, and FMM can be described as having two phases. In the precompute phase, the particles are partitioned into a hierarchical tree of clusters, and assuming the particle distribution is homogeneous, this phase generally scales like $O(N\log N)$. In the compute phase, well-separated particle-particle interactions are approximated using for example monopole approximations [3, 4] or higher order multipole expansions [5, 10, 11]. Some differences in the compute phase of these methods are noted as follows. The TC traverses the tree for each target particle to identify well-separated particle-cluster pairs, while the DTT traverses two copies of the tree simultaneously to identify well-separated cluster-cluster pairs; both methods use a multipole acceptance criterion (MAC) for this purpose. The FMM passes information from one level to the next, with a uniform definition of the interaction list at each level; an upward pass computes cluster moments, and a downward pass computes potentials using multipole-to-local and local-to-local translations. The compute phase of the TC scales like $O(N\log N)$, while the compute phase of the FMM [5, 6] and DTT [12, 11] scale like $O(N)$. Next we review developments relevant to the present work in more detail, including dual tree traversal methods, kernel-independent methods, and GPU implementations.

The DDT methods employ cluster-cluster approximations. Early work by Appel developed a DTT-style method using monopole approximations [3]. A later algorithm independently developed by Dehnen for stellar dynamics used Taylor approximations [11]. Subsequent large-scale astrophysical simulations implemented parallel DTTs on many-core architectures [13, 14, 15, 16], and the DTT with Taylor approximations has been applied in molecular dynamics simulations of periodic condensed phase systems [17] and for polarizable force fields [18].

Tree-based fast summation methods originally relied on analytic series expansions specific to a given kernel $G({\bf x},{\bf y})$. For example, in the case of the Coulomb and Yukawa kernels, multipole expansions were used in the FMM [6, 19], and Cartesian Taylor expansions were used in the TC [20, 21]. Alternative approximations for the Coulomb kernel were introduced involving numerical discretization of the Poisson integral formula [22] and multipole expansions at pseudoparticles [23]. Eventually kernel-independent methods were developed that require only kernel evaluations and are suitable for a large class of kernels. Among these, the kernel-independent FMM (KIFMM) uses equivalent densities defined on proxy surfaces [24], the black-box FMM (bbFMM) uses polynomial interpolation and SVD compression [25], and the barycentric Lagrange treecode (BLTC) uses barycentric Lagrange interpolation [26, 27]. A number of related proxy point methods have recently been developed using skeletonized interpolation [28] and interpolative decomposition [29, 30]. Several of the kernel-independent fast summation methods have been parallelized for multi-core CPU systems [31, 32, 33, 34, 35, 36, 30].

Many modern HPC systems rely heavily on graphics processing units (GPUs) for high-throughput arithmetic. The direct sum in Eq. (1) is well suited for GPU computing; this is because the kernel evaluations $G({\bf x}_{i},{\bf x}_{j})$ can be computed concurrently, where the targets ${\bf x}_{i}$ provide an outer level of parallelism, while the sources ${\bf x}_{j}$ provide an inner level of parallelism. Early GPU implementations of direct summation achieved a 25$\times$ speedup over an optimized CPU implementation [37] and a 250$\times$ speedup over a portable C implementation [38]. However, these codes still scale like $O(N^{2})$ and there is great interest in implementing the sub-quadratic scaling tree-based methods on GPUs, although this is challenging due to the complexity of these methods in comparison with direct summation. In this direction there have been several GPU implementations of the TC and FMM [39, 40, 41, 42, 36, 43, 44, 45, 46, 47], but we know of only a few GPU implementations of the DTT [48, 49].

The present work contributes a GPU-accelerated tree-based fast summation method called barycentric Lagrange dual tree traversal (BLDTT). The BLDTT employs several techniques from previous tree-based methods, including (1) the DTT algorithmic structure [3, 11], (2) barycentric Lagrange interpolation [26, 27, 47], and (3) upward and downward passes similar to those in the FMM, but adapted to the context of barycentric Lagrange interpolation. The BLDTT has several additional features that should be noted. The algorithm replaces well-separated particle-particle interactions by adaptively chosen particle-cluster, cluster-particle, and cluster-cluster approximations, where the clusters are represented by proxy particles at Chebyshev grid points. The approximations are done with barycentric Lagrange interpolation and require only kernel evaluations, hence the BLDTT is kernel-independent. Similar to other tree-based fast summation methods, the BLDTT has a precompute phase and a compute phase; the precompute phase scales like $O(N\log N)$ with a small prefactor, while the compute phase scales like $O(N)$, so the observed scaling of the BLDTT is essentially $O(N)$.

As will be shown, the barycentric Lagrange approximations resemble the direct sum in Eq. (1) and they can be efficiently mapped onto GPUs. Based on this observation we present an OpenACC GPU implementation of the BLDTT with MPI remote memory access for distributed memory parallelization. The performance of the BLDTT is documented for calculations with different problem sizes, particle distributions, geometric domains, and interaction kernels, unequal target and source particles. Comparison with our earlier particle-cluster barycentric Lagrange treecode (BLTC) shows the superior performance of the BLDTT. In particular, on a single GPU for problem sizes ranging from $N$=1E5 to 1E8, the BLTC has $O(N\log N)$ scaling, while the BLDTT has $O(N)$ scaling. In addition, MPI strong scaling results are presented for the BLTC and BLDTT using $N$=64E6 particles on up to 32 GPUs. The BLDTT code is part of the BaryTree library for fast summation of particle interactions available on GitHub [50].

The remainder of the paper is organized as follows. Section 2 describes the barycentric Lagrange dual tree traversal (BLDTT) fast summation method. Section 3 describes our implementation of the BLDTT using MPI remote memory access for distributed memory parallelization and OpenACC for GPU acceleration. Section 4 presents numerical results for several test cases. Section 5 gives the conclusions.

We briefly review the barycentric Lagrange form of polynomial interpolation in 1d [26]. Given a function $f(x)$ and $n+1$ points $s_{k},k=0:n$, the Lagrange form of the interpolating polynomial is

where the Lagrange polynomial $L_{k}(x)$ in barycentric form is

This work employs Chebyshev points of the second kind,

in which case the interpolation weights are given by

where $\delta_{k}=1/2$ if $k=0$ or $n$, and $\delta_{k}=1$ otherwise [51, 26]. The algorithms described below use barycentric Lagrange interpolation in 3D rectangular boxes with interpolation points ${\bf s}_{\bf k}=(s_{k_{1}},s_{k_{2}},s_{k_{3}})$ given by a Cartesian tensor product grid of Chebyshev points adapted to the box. Depending on the context, the interpolation points ${\bf s}_{\bf k}$ may also be referred to as proxy particles.

The BLDTT fast summation method computes the potential at a set of $M$ target particles due to interactions with a set of $N$ source particles; Eq. (1) is a special case in which the targets and sources refer to the same set of particles. First, two hierarchical trees of particle clusters are built, one for the target particles and one for the source particles, where each cluster is a rectangular box; clusters in the target tree are denoted $C_{t}$ and clusters in the source tree are denoted $C_{s}$. The computed potential at a target particle ${\bf x}_{i}$ has contributions from four types of interactions as determined by the dual tree traversal algorithm described below. The four types are direct particle-particle (PP) interactions of nearby particles, and particle-cluster (PC), cluster-particle (CP), and cluster-cluster (CC) approximations of well-separated clusters.

Algorithm 1 is a high-level overview of the BLDTT. Lines 1-4 describe the input consisting of target and source particle data, interpolation degree $n$, MAC parameter $\theta$, and the maximum number of particles in the leaves of each tree, $M_{0},N_{0}$. Line 5 describes the output potentials. Line 6 builds the target tree and source tree containing target clusters $C_{t}$ and source clusters $C_{s}$. Line 7 is the upward pass to compute proxy charges $q_{\bf k}$ at proxy particles ${\bf s}_{\bf k}$ in source clusters $C_{s}$. Line 8 is the dual tree traversal to compute $PP,PC,CP,CC$ interactions. Line 9 is the downward pass to interpolate potentials from proxy particles ${\bf t}_{\bm{\ell}}$ to target particles ${\bf x}_{i}$ in target clusters $C_{t}$. The steps will be described in detail below. We summarize the notation used in presenting the BLDTT in Table 1.

The target and source trees are constructed by the same routines, described here for the target tree. The maximum number of particles per leaf is a user-specified parameter, $M_{0}$ for the target tree and $N_{0}$ for the source tree. The root cluster is the minimal bounding box containing all target particles. The root is recursively divided into child clusters, terminating when a cluster contains fewer than $M_{0}$ particles. Division occurs at the midpoint of the cluster; in general the cluster is bisected in all three coordinate directions, resulting in eight child clusters, with two exceptions. First, a cluster is divided into only two or four children in order to maintain a good aspect ratio, that is, a ratio of longest to shortest side lengths no greater than $\sqrt{2}$. Second, a cluster is divided into only two or four children to avoid creating leaf clusters with fewer than $M_{0}/2$ particles on average; in particular, if a cluster contains between $M_{0}$ and $2M_{0}$ particles, it is divided into two children, and if it contains between $2M_{0}$ and $4M_{0}$ particles, it is divided into four children. Upon creation, each cluster is shrunk to the minimal bounding box containing its particles, and a tensor product grid of Chebyshev points adapted to the box is created; these are also referred to as proxy particles. After building the trees, the BLDTT performs the upward pass, dual tree traversal, and downward pass, but before discussing these steps, the next subsection describes the four types of interactions which are eventually combined to compute potentials.

Figure 1 depicts the four types of interactions between a target cluster $C_{t}$ (left, blue) and a source cluster $C_{s}$ (right, red), where dots are target/source particles ${\bf x}_{i},{\bf y}_{j}$, and crosses are target/source proxy particles ${\bf t}_{\bm{\ell}},{\bf s}_{\bf k}$. Also shown are the target/source cluster radii $r_{t},r_{s}$, and the target-source cluster distance $R$. These diagrams depict 2D versions of the interactions; in practice the particles are distributed in 3D and the clusters are rectangular boxes. In general, clusters can interact via their particles (dots) or their proxy particles (crosses). Figure 1 shows the four cases: (a) $C_{t}$ and $C_{s}$ use particles (PP), (b) $C_{t}$ uses particles and $C_{s}$ uses proxy particles (PC), (c) $C_{t}$ uses proxy particles and $C_{s}$ uses particles (CP), (d) $C_{t}$ and $C_{s}$ both use proxy particles (CC). The interactions are described in detail below; to simplify notation we present the interactions in 1D and note that the extension to 3D is straightforward using tensor products.

Particle-particle interaction. Figure 1a depicts a PP interaction. In this case the PP potential at a target particle $x_{i}\in C_{t}$ due to direct interaction with the source particles $y_{j}\in C_{s}$ is denoted by

Particle-cluster approximation. Figure 1b depicts a PC interaction. The kernel is approximated by holding $x_{i}$ fixed and interpolating with respect to $y_{j}$ at the proxy particles $s_{k}$ in $C_{s}$,

Substituting into the particle-particle interaction and rearranging terms yields the PC potential,

where the proxy charges $\widehat{q}_{k}$ of the proxy particles $s_{k}$ are

Equation (8) uses $C_{t},\widehat{C}_{s}$ to indicate that the target particles $x_{i}$ interact with the proxy source particles $s_{k}$.

Cluster-particle approximation. Figure 1c depicts a CP interaction. The kernel is approximated by interpolating with respect to $x_{i}$ at the proxy particles $t_{\ell}$ in $C_{t}$ and holding $y_{j}$ fixed,

Substituting into the particle-particle interaction and rearranging terms yields the CP potential,

where the CP proxy potential $\phi(t_{\ell},\widehat{C}_{t},C_{s})$ is

Equation (11) and Eq. (12) use $\widehat{C}_{t},C_{s}$ to indicate that the proxy target particles of $C_{t}$ interact with the source particles $y_{j}$. Equation (11) interpolates from the proxy target particles $t_{\ell}$ to the target particles $x_{i}$.

Cluster-cluster interaction. Figure 1d depicts a CC interaction. The kernel is interpolated with respect to $x_{i}$ at the proxy particles $t_{\ell}$ in $C_{t}$ and with respect to $y_{j}$ at the proxy particles $s_{k}$ in $C_{s}$,

Substituting into the particle-particle interaction and rearranging terms yields the CC potential,

where the CC proxy potential $\phi(t_{\ell},\widehat{C}_{t},\widehat{C}_{s})$ is

and the proxy charges $\widehat{q}_{k}$ were defined in Eq. (9). Equation (14) and Eq. (15) use $\widehat{C}_{t},\widehat{C}_{s}$ to indicate that the proxy target particles $t_{\ell}$ interact with the proxy source particles $s_{k}$. Equation (14) interpolates from the proxy target particles $t_{\ell}$ to the target particles $x_{i}$.

The following three subsections describe the rest of Algorithm 1, comprising the upward pass, dual tree traversal, and downward pass.

The upward pass computes the proxy charges $\widehat{q}_{k}$ defined in Eq. (9) for the proxy particles $s_{k}$ in each source cluster $\widehat{C}_{s}$ in the source tree, as required in the PC and CC approximations in Eqs. (8) and (15). Note that each source particle $y_{j}$ contributes to the proxy charges $\widehat{q}_{k}$ of exactly one cluster at a given level of the tree. Hence with $N$ source particles and tree depth $O(\log N)$, computing the proxy charges directly by Eq. (9) requires $O(N\log N)$ operations.

The BLDTT algorithm uses an alternative approach as follows. Let $C_{s}$ denote a parent cluster with Lagrange polynomials $L_{k}(y)$ and interpolation points $s_{k}$, and let $C_{s}^{i},i=1:8$ denote the eight child clusters with Lagrange polynomials $L^{i}_{k_{i}}(y)$ and interpolation points $s_{k_{i}}$. As an alternative to Eq. (9), the proxy charges of the parent $\widehat{q}_{k}$ are computed from the proxy charges of the children $\widehat{q}_{k_{i}}$ by the expression

which is derived in B. The derivation uses the definitions of the parent and child proxy charges and the relation

which is a special case of Eq. (2).

The upward pass starts by computing the proxy charges of the leaves of the source tree using Eq. (9) and then ascending to the root by Eq. (16). This is analogous to the upward pass in the FMM [5] where the multipole moments of a parent cluster are obtained from the moments of its children. As with the FMM, computing the proxy charges this way requires $O(N)$ operations, which can be seen as follows. Computing the proxy charges for the leaves requires $O(n^{3}N)$ operations; each of the $N$ source particles contributes to one leaf, which contains $O(n^{3})$ proxy particles (in 3D). Then, evaluating the child-to-parent relation in Eq. (16) for each parent proxy charge requires $O(n^{6})$ operations and is independent of the number of source particles. Since there are $O(N)$ clusters in the tree, ascending the tree requires an additional $O(n^{6}N)$ operations. Hence the operation count for the BLDTT upward pass is $O(n^{3}N)+O(n^{6}N)=O(N)$, as it is for the FMM.

The dual tree traversal determines which pairs of clusters in the target and source trees interact by one of the four options described above (PP, PC, CP, CC). Before the traversal starts, two sets of potentials are initialized to zero, potentials $\phi(x_{i})$ at the target particles and potentials $\phi(t_{\ell})$ at the proxy target particles. In the course of the traversal, the potentials $\phi(x_{i})$ are incremented due to PP and PC interactions, and the potentials $\phi(t_{\ell})$ are incremented due to CP and CC interactions. Following the dual tree traversal, the $\phi(t_{\ell})$ are interpolated to the target particles $x_{i}$ and combined with the $\phi(x_{i})$ in the downward pass.

The dual tree traversal uses the recursive procedure DTT($C_{t},C_{s}$) described in Algorithm 2, which takes a target cluster $C_{t}$ and a source cluster $C_{s}$ as input. Initially the procedure is called for the root clusters of the target and source trees. In what follows, the clusters are considered to be well-separated if $(r_{t}+r_{s})/R<\theta$, where $r_{t},r_{s}$ are the target and source cluster radii and $R$ is the center-center distance between the clusters.

If $C_{t}$ and $C_{s}$ are well-separated (line 2), then they interact in one of four ways depending on the number of particles they contain relative to the number of proxy particles in a cluster, which is denoted by $n_{p}=(n+1)^{3}$ in 3D. If $C_{t}$ and $C_{s}$ are both large (lines 3-4), then the CC proxy potentials are incremented using Eq. (15); else if $C_{t}$ is large and $C_{s}$ is small (lines 5-6), then the CP proxy potentials are incremented using Eq. (12); else if $C_{t}$ is small and $C_{s}$ is large (lines 7-8), then the PC potentials are incremented using Eq. (8); else $C_{t}$ and $C_{s}$ are both small (line 9) and the PP potentials are incremented using Eq. (6).

If $C_{t}$ and $C_{s}$ are not well-separated, then the traversal continues as follows. If $C_{t}$ and $C_{s}$ are leaves (lines 11-12), then the PP potentials are incremented using Eq. (6). Otherwise if $C_{s}$ is a leaf, then it interacts recursively with the children of $C_{t}$ (line 13), while if $C_{t}$ is a leaf, then it interacts recursively with the children of $C_{s}$ (line 14). Finally if $C_{t}$ and $C_{s}$ are both not leaves, then the smaller cluster interacts recursively with the children of the larger cluster (lines 15-17).

The DTT yields potentials $\phi(x_{i})$ due to PP and PC interactions and proxy potentials $\phi(t_{\ell})$ due to CP and CC interactions. In the case of $N$ homogeneously distributed source and target particles, the operation count of the dual tree traversal has been shown to be $O(N)$ [12, 11].

The downward pass interpolates the proxy potentials $\phi(t_{\ell})$ to the target particles $x_{i}$ and increments the potentials $\phi(x_{i})$. This can be done in two ways as described below.

First note that each target particle $x_{i}$ is contained in a chain of target clusters,

where the superscript denotes the level in the target tree. The target cluster $C_{t}^{m}$ at level $m$ in the chain has Lagrange polynomials $L^{m}_{k_{m}}(x_{i})$ and proxy potentials $\phi(t_{k_{m}}^{m})$ that contribute to $\phi(x_{i})$,

where $t^{m}_{k_{m}}$ are the proxy particles of $C_{t}^{m}$, and the $\mathrel{+}=$ indicates that the results are aggregated with the potentials $\phi(x_{i})$ due to PP and PC interactions previously computed in the DTT. In Eq. (19) the inner sum interpolates potential values from the proxy particles $t^{m}_{k_{m}}$ to the target particle $x_{i}$, and the outer sum accumulates the results from each level in the tree. Computing $\phi(x_{i})$ as indicated in Eq. (19) requires $O(M\log M)$ operations; the factor $M$ is the number of target particles $x_{i}$, the factor $\log M$ is the number of levels in the target tree, and the inner sum requires $O(n^{3})$ operations independent of $M$.

Rather than interpolating from the proxy particles $t_{m}^{k}$ directly to the target particle $x_{i}$ as in Eq. (19), we utilize a recursive alternative. In what follows, $C_{t}^{m}$ is a parent cluster at level $m$ and $C^{m-1}_{t}$ is a child cluster at level $m-1$. The procedure interpolates the parent proxy potentials $\phi(t_{k_{m}}^{m})$ to child proxy potentials $\phi(t_{k_{m-1}}^{m-1})$,

where the $\mathrel{{+}{=}}$ indicates that the results of the parent-to-child interpolation on the right are aggregated with the child proxy potentials $\phi(t_{k_{m-1}}^{m-1})$ due to CP and CC interactions previously computed in the DTT. This procedure starts with the root cluster of the target tree (level $m=L$) and descends to the leaves (level $m=1$). Upon reaching the leaves, the proxy potentials $\phi(t_{k_{1}}^{1})$ are interpolated to the target particles $x_{i}$ and aggregated with the PP and PC potentials previously computed in the DTT,

It should be noted that the expressions for $\phi(x_{i})$ in Eq. (19) and Eq. (21) are equivalent and this is shown in C.

The recursive form of the downward pass described here is similar to the local-to-local step in the FMM, where the local coefficients are shifted and accumulated, and as in that case the operation count is reduced to $O(M)$. In particular, the parent-to-child interpolation in Eq. (20) requires $O(n^{6})$ operations (in 3D), independent of the number of target particles $M$, and the tree contains $O(M)$ clusters, so interpolation from the root down to the leaves by Eq. (20) requires $O(n^{6}M)$ operations. Then the final interpolation from the leaf proxy particles to the target particles by Eq. (21) requires $O(n^{3})$ operations for each target, yielding complexity $O(n^{6}M)+O(n^{3}M)=O(M)$ for the downward pass.

We briefly describe our previous barycentric Lagrange treecode (BLTC) [27, 47] which has an algorithmic structure resembling the Barnes-Hut treecode [4]. Unlike the BLDTT which builds a tree on both the source and target particles, the BLTC builds a tree of clusters on the source particles and a set of batches on the target particles, where the batches correspond to the leaves of a target tree. Once the source tree is built, the BLTC computes the proxy charges for each source cluster directly from the source particles by Eq. (9). The source tree is then traversed for each target batch, starting at the root of the tree and checking whether a given source cluster and target batch are well-separated. If they are well-separated and the source cluster contains more particles than interpolation points, then the cluster and batch interact by the PC approximation in Eq. (8). If they are not well-separated, then the batch interacts with the children of the source cluster. Leaves in the source tree that are not well-separated from a given target batch, and source clusters that are well-separated but contain more interpolation points than particles interact directly with the target batch by the PP interaction in Eq. (6). For $M$ target particles and $N$ source particles, the BLTC operation count is $O(N\log N)$ + $O(M\log N)$, where the first term arises from the computation of the proxy charges and the second term arises from the source tree traversal.

We decouple dual tree traversal from the computation of the particle interactions. The dual tree traversal is performed on the CPU, creating four interaction lists for each cluster of the target tree, one for each type of interaction. Each list contains the indices of the source clusters that interact with the target cluster by the given interaction type. The interactions are then computed by directly iterating over the interaction lists; this improves the efficiency of GPU calculations because these lists can be iterated over rapidly and GPU compute kernels (described below) can be queued asynchronously.

To implement distributed memory parallelization, we use locally essential trees (LET) [52]. Particles are partitioned by recursive coordinate bisection to create compact sub-domains on each MPI rank using Trilinos Zoltan [53, 54]. Each MPI rank constructs the local source tree and local target tree for its particles. The LET of a rank is the union of the rank’s local source tree and all source clusters from remote ranks interacting with its local target tree. Even though constructing the LETs requires an all-to-all communication, the amount of data acquired by each rank grows only logarithmically with the problem size [52]. The construction and communication required by the LETs are performed using MPI passive target synchronization remote memory access (RMA). RMA is a one-sided communication model within MPI in which an origin process can put data onto a target process or get data from a target process through specially declared memory windows, with no active involvement from the target process. This approach enables each rank to construct its LET asynchronously from other ranks.

The GPU implementation uses eight compute kernels, two for the upward pass, four for computing interactions determined by the dual tree traversal, and two for the downward pass. The compute kernels are generated with OpenACC directives, compiled with the PGI C compiler. We utilize asynchronous launch of kernels in multiple GPU streams to hide as much latency as possible. The approach described here generalizes to multiple GPUs in a straightforward manner, in which each GPU corresponds to one MPI rank.

The first upward pass kernel computes the proxy charges for a given leaf in the source tree. For each leaf, this kernel is launched asynchronously, and any further computation is blocked until all leaf proxy charges are computed. The second upward pass kernel computes the proxy charges of a parent cluster using the proxy charges of its children. For a given level of the tree, this kernel is launched asynchronously for each cluster at that level, and any further computation is blocked until all proxy charges at that level are computed. The proxy charges at a given level are computed after the computation of the proxy charges at the previous level is complete.

The four DTT kernels compute the interaction of a target cluster with a source cluster. Each PP, PC, CP, and CC interaction launches one compute kernel. All such kernels are launched asynchronously and further computation is blocked until these kernels complete. The four interaction kernels have a similar structure in which the outer loop is over the target particles or proxy particles in the target cluster, and the inner loop is over the source particles or proxy particles in the source cluster. Importantly, due to the Lagrange form of barycentric interpolation, the inner loop iterations are independent, unlike alternative approximation methods which are sequential. The outer loop is mapped to the gang construct in OpenACC and the inner loop is mapped to the vector construct. Conceptually, a member of a gang corresponds to a thread block, while a member of a vector corresponds to an individual thread.

The two downward pass compute kernels are similar in structure to the upward pass kernels. At each level of the target tree above the leaves, beginning with the root, the first downward pass compute kernel is launched asynchronously for each target cluster at that level to interpolate the proxy potentials of the cluster to its children. Further computation is blocked until all compute kernels at a given level complete. Finally, the second downward pass compute kernel is launched asynchronously for each target leaf to interpolate the proxy potentials to the target particles.

Figure 2 shows the compute time (s) for direct summation (green), BLTC (red), and BLDTT (blue) with $N$=1E5, 1E6, 1E7, 1E8 random uniformly distributed source and target particles in the $[-1,1]^{3}$ cube interacting by the Coulomb kernel. The BLTC and BLDTT use MAC parameter $\theta=0.7$ and interpolation degree $n=8$, yielding 7-8 digit accuracy. Figure 2(a) is a linear plot, showing that the BLDTT is about twice as fast as the BLTC, and both are much faster than direct summation. Figure 2(b) is a logarithmic plot with reference lines scaling as $O(N)$ (dashed), $O(N\log N)$ (dotted), and $O(N^{2})$ (dash-dotted), showing that as the problem size increases, the BLTC has asymptotic $O(N\log N)$ scaling, while the BLDTT has asymptotic $O(N)$ scaling, as expected. Table 2 records the compute time and error; the asymptotic scaling of the compute time can be quantitatively confirmed, and while there is a slight increase in the error with problem size, the BLDTT error is consistently smaller than the BLTC error.

In this subsection we compare the BLDTT running on a single NVIDIA P100 GPU to running on 8 CPU cores of an Intel Xeon E5-2680v3 processor at 2.50 GHz with MPI parallelization. We perform the same four calculations above, so the errors are the same as those in Table 2. Table 3 gives the compute times, showing that the BLDTT achieves 30-40$\times$ speedup on the GPU compared the 8 CPU cores. The efficiency of the BLDTT running on the GPU is due to the independent nature of the terms in the barycentric Lagrange approximation, which allows them to be computed concurrently.

We investigate the performance of the BLDTT for three different random particle distributions: (a) uniform particles in $[-1,1]^{3}$, (b) Gaussian particles with radial pdf $\frac{1}{\sqrt{6\pi}}\exp\left(-r^{2}/6\right)$, (c) Plummer particles [55, 56] with radial pdf $\frac{3}{4\pi}\left(1+r^{2}\right)^{-5/2}$ and cutoff at $\pm 100$ in all three Cartesian coordinates. The charges of the uniform and Gaussian particles are uniformly distributed in $[-1,1]$, while the Plummer particle charges are set to $1/N$, where $N$ is the total number of particles. To give a sense of the structure of the distributions, Figure 3 depicts the three distributions with $N$=4E5 particles. Compared to the uniform case (a), the Gaussian and Plummer distributions (b,c) are more highly concentrated near the origin, with the Gaussian decaying more rapidly away from the origin and the Plummer decaying less rapidly.

Figure 4 shows the compute time (s) versus relative $\ell_{2}$ error for the BLDTT (blue, solid) and BLTC (red, dashed) on these three distributions with $N$=2E7 particles. Each connected curve represents constant MAC with $\theta=0.5$ ($\times$), $\theta=0.7$ ($\circ$), $\theta=0.9$ ($*$), and the interpolation degree $n=1,2,4,6,8,10$ increases from right to left on each curve. For these parameter choices the errors span the range from 1 digit to 10 digit accuracy. Large $\theta$ is more efficient for low accuracy and small $\theta$ is more efficient for high accuracy. The results show that the BLDTT has consistently better performance than the BLTC and is less sensitive to the particle distribution.

To demonstrate the effect of including PC and CP interactions in the BLDTT, Figure 5 shows the compute time versus relative $\ell_{2}$ error for the BLDTT (blue, solid) as presented in this paper using PP, PC, CP and CC interactions, and a version of the BLDTT (red, dashed) using only CC and PP interactions. When only CC and PP interactions are used, the interaction between a target cluster and a source cluster is handled by PP interaction if either cluster contains more interpolation points than particles, whereas the flexibility to choose PC or CP interactions in those cases improves performance at higher interpolation degree, especially for the non-uniform particle distributions.

To further understand the effect of including PC and CP interactions, next we compare the number of pointwise interactions used by the two versions of the BLDTT, where by pointwise interaction we mean one evaluation of the kernel $G({\bf x},{\bf y})$. Results are shown for MAC $\theta=0.9$ and interpolation degree $n=1,2,4,6,8,10$, for the same three random distributions with $N$=2E7 particles as above. Figure 6 displays results for the four types of interactions in stacked bars, CC (blue), PP (orange), PC (yellow), CP (purple), from bottom to top, where the left bar in each pair is the BLDTT with CC and PP interactions only, and the right bar is the BLDTT with PP, PC, CP, and CC interactions. In this case a direct sum calculation would use 4E14 PP interactions, while the BLDTT calculations use less than 6E12 interactions. The results show that for high degree, introducing PC and CP interactions into the BLDTT significantly reduces the number of PP interactions, replacing them with a much smaller number of PC and CP interactions, and this effect is more prominent for the nonuniform particle distributions.

We demonstrate here the performance of the BLDTT on three examples with non-cubic particle domains depicted in Figure 7: (a) thin slab of dimensions $1\times 10\times 10$, (b) square rod of dimensions $1\times 1\times 10$, and (c) spherical surface of radius $1$. In all cases the particles are random uniformly distributed.

Figure 8 shows the compute time (s) versus the relative $\ell_{2}$ error for the BLDTT (blue, solid) and BLTC (red, dashed) on these three examples with $N$=2E7 particles, using MAC $\theta=0.5,0.7,0.9$ and interpolation degree $n=1,2,4,6,8,10$. The results show that the BLDTT has consistently better performance than the BLTC. Compared to the cubic domain results in Figure 4(a), the BLDTT achieves similar levels of error and runs somewhat faster for the non-cubic domains. Heuristically, the BLDTT run time depends on the complexity of the tree; in particular, the tree is an oct-tree for the cubic domain, close to a quad-tree for the thin slab and sphere surface, and close to a binary tree for the square rod. The results indicate that BLDTT automatically adapts to the complexity of the tree without requiring explicit reprogramming.

To demonstrate the performance of the BLDTT with unequal target and source particles, we consider a cluster-particle variant of the BLTC for comparison [57, 58]. The CP-BLTC builds a tree on the $M$ targets and a set of batches on the $N$ sources, and rather than using PP and PC interactions, it uses PP and CP interactions. Instead of an $O(N\log N)$ upward pass to compute proxy charges, the CP-BLTC has an $O(M\log M)$ downward pass to interpolate proxy potentials to targets. While the compute phase of the BLTC is $O(M\log N)$, the compute phase of the CP-BLTC is $O(N\log M)$. Since the compute phase is in general the most expensive part of the algorithm, we expect the BLTC to perform better than the CP-BLTC when $N>M$, and vice versa.

Figure 9 shows the compute time (s) versus relative $\ell_{2}$ error for the BLTC (red, dashed), CP-BLTC (green, dash-dotted), and BLDTT (blue, solid) with (a) $M$=2E7 targets, $N$=2E6 sources, (b) $M$=2E6 targets, $N$=2E7 sources, for MAC $\theta$ and interpolation degree $n$ as above. For (a) $M>N$, the CP-BLTC outperforms the BLTC, for (b) $N>M$, the BLTC outperforms the CP-BLTC for errors below 1E-5, while the BLDTT outperforms the two treecodes in both cases. Note however in (b) that for MAC $\theta=0.9$ and error larger than 1E-3, the CP-BLTC runs slightly faster than the BLDTT; this is due to the cost of the upward pass in the BLDTT, which is more expensive than the downward pass, and with low degree $n$, makes up a substantial portion of the compute time. The results demonstrate the ability of the BLDTT to efficiently adapt to the case of unequal targets and sources.

In previous sections we considered particles interacting through the Coulomb potential. Here we demomstrate the performance of the BLDTT on three other interaction kernels: (a) an oscillatory kernel, $\sin(\pi r)/r$, (b) a Yukawa kernel, $\exp(-0.5r)/r$, and (c) a regularized Coulomb kernel, $1/(r^{2}+\epsilon^{2})^{1/2}$, with $\epsilon=0.005$. Figure 10 shows the compute time (s) versus relative $\ell_{2}$ error for the BLDTT (blue, solid) and BLTC (red, dashed) for these three kernels on $N$=2E7 random uniformly distributed particles in the cube $[-1,1]^{3}$ for various values of the MAC $\theta$ and interpolation degree $n$. The results show that the BLDTT has consistently better performance than the BLTC. Comparing with the Coulomb potential resuts in Figure 4(a), we see that the BLDTT has similar performance for the various interaction kernels, reflecting the kernel-independent nature of the method.

Finally, we demonstrate the MPI strong scaling of the BLDTT up to 32 NVIDIA P100 GPUs with one MPI rank per GPU. The particles are partitioned into geometrically localized domains by Trilinos Zoltan [53, 54]. Figure 11 depicts a sample domain decomposition for $N$=1.6E6 random uniformly distributed particles in the cube $[-1,1]^{3}$, (a) across 8 ranks with 2E5 particles per rank, and (b) across 16 ranks with 1E5 particles per rank. Colors represent particles residing on different ranks.

We consider a problem with $N$=64E6 particles using MAC $\theta=0.7$ and interpolation degree $n=8$ yielding error $\approx$ 1E-8. Figure 12 shows the compute time versus the number of GPUs for (a) random uniform, (b) Gaussian, and (c) Plummer distributions for the BLDTT (blue) and BLTC (red), where dashed lines indicate ideal scaling and the boxed numbers show the parallel efficiency. As was shown earlier for one GPU in Figure 4, the BLDTT is consistently faster than the BLTC up to 32 GPUs, and the speedup improves for the nonuniform particle distributions. The BLDTT and BLTC have generally similar parallel efficiency for all three distributions; for example on 32 GPUs, the BLDTT has parallel efficiency 77%, 65%, and 66% for the uniform, Gaussian, and Plummer distributions, compared to 83%, 81%, and 64% for the BLTC.

The slightly lower parallel efficiency for the BLDTT compared to the BLTC can be explained by the greater efficiency of the BLDTT algorithm itself. Figure 13 shows the component breakdown as a percentage of run time (total wall clock time) of the (a) BLTC and (b) BLDTT for the uniform distribution results in Figure 12(a). The components shown are the upward pass (blue), compute due to local sources and source clusters (orange), compute due to remote sources and source clusters (yellow), downward pass (purple), LET construction and communication (green), and other (light blue), which includes tree building and interaction list building. The breakdown is based on timing results for the most expensive MPI rank in each case. Note that the LET construction accounts for an increasing percentage of the run time as the number of ranks increases. This is to be expected as more particles reside on remote ranks and must be communicated, and this is the primary factor that impedes ideal parallel scaling. The LET construction time is nearly identical for the BLTC and BLDTT, however since the BLDTT computations are more efficient, the LET construction accounts for a larger percentage of the run time and this results in the lower parallel efficiency seen in Figure 12 as the number of ranks increases.