Igor Chollet
email: igor.chollet@inria.fr

Xavier Claeys
email: claeys@ann.jussieu.fr

Pierre Fortin
email: pierre.fortin@univ-lille.fr

Laura Grigori
email: laura.grigori@inria.fr

A Directional Equispaced interpolation-based Fast Multipole Method for oscillatory kernels

Fast Multipole Methods (FMMs) based on the oscillatory Helmholtz kernel can reduce the cost of solving N-body problems arising from Boundary Integral Equations (BIEs) in acoustic or electromagnetics. However, their cost strongly increases in the high-frequency regime. This paper introduces a new directional FMM for oscillatory kernels (defmm -directional equispaced interpolation-based fmm), whose precomputation and application are FFT-accelerated due to polynomial interpolations on equispaced grids. We demonstrate the consistency of our FFT approach, and show how symmetries can be exploited in the Fourier domain. We also describe the algorithmic design of defmm, well-suited for the BIE non-uniform particle distributions, and present performance optimizations on one CPU core. Finally, we exhibit important performance gains on all test cases for defmm over a state-of-the-art FMM library for oscillatory kernels.

where G : R d × R d → C. Such an N-body problem appears in the numerical solving of Boundary Integral Equations (BIEs). We are especially concerned by the oscillatory Helmholtz kernel involved in BIEs applied to acoustic or electromagnetics G(x, y) := e iκ|x-y| 4π|x -y| , where κ ∈ R + is named the wavenumber, i denotes the complex number and | • | refers to the Euclidian norm. Computing p using hierarchical methods can achieve an O(N log N ) complexity on large surface particle distributions (i.e. point clouds), that are the distributions generated when discretizing BIEs. However, the cost of the hierarchical methods can still be a bottleneck in the highfrequency regime, i.e. when κD ≫ 1, D denoting the side length of the smallest box encompassing X and Y . We aim here at mitigating this cost in the high-frequency regime for the non-uniform particle distributions that are typical of BIE problems.

Related work

The direct evaluation of the N-body problem (1) has a complexity O(N 2 ). This complexity can be reduced using a hierarchical method, such as the Fast Multipole Method (FMM) [START_REF] Cheng | A wideband fast multipole method for the Helmholtz equation in three dimensions[END_REF][START_REF] Greengard | A fast algorithm for particle simulations[END_REF]. Thanks to a hierarchical decomposition, such as a 2 d -tree1 representation of the (cubical) computational domain B, and to a multipole acceptance criterion, the computation of ( 1) is indeed divided into two parts: a near field one which is computed directly and a far field one which is approximated through multipole and local expansions. Nodes of such a 2 d -tree are named cells and correspond to cubical subdomains of B whose radii depend on the tree level (see Figure 1). There mainly exist two approaches for FMMs dealing with the high-frequency Helmholtz kernel. First, the kernel-explicit methods consist in an explicit diagonalization of large far-field matrix blocks using analytic kernel expansions [START_REF] Cecka | Fourier-based Fast Multipole Method for the Helmholtz equation[END_REF][START_REF] Cheng | A wideband fast multipole method for the Helmholtz equation in three dimensions[END_REF], considering that the N -body problem in equation 1 is interpreted as a matrix-vector product. These methods perform differently depending on the particle distribution, with complexities ranging from O(N ) on quasi-uniform distributions to O(N logN ) or O(N 3/2 ) on homogeneous surface distributions (see [START_REF] Cecka | Fourier-based Fast Multipole Method for the Helmholtz equation[END_REF]), and possibly even to O(N2 ) on highly non-uniform distributions (see [START_REF] Messner | Fast directional multilevel summation for oscillatory kernels based on Chebyshev interpolation[END_REF]). Second, the kernel-independent methods (see for instance [START_REF] Fong | The black-box fast multipole method[END_REF][START_REF] Ying | A kernel-independent adaptive fast multipole algorithm in two and three dimensions[END_REF]) allow to use the same formalism, hence the same implementation, for a given class of kernels. Based on a specific definition of the well-separateness (i.e. the subsets of B × B on which the far-field approximations can be applied), the directional kernel-independent approaches exploit the low-rank property of particular blocks of a modified kernel [START_REF] Brandt | Multilevel computations of integral transforms and particle interactions with oscillatory kernels[END_REF][START_REF] Engquist | Fast directional multilevel algorithms for oscillatory kernels[END_REF][START_REF] Messner | Fast directional multilevel summation for oscillatory kernels based on Chebyshev interpolation[END_REF]. As opposed to the kernelexplicit methods, the kernel-independent directional approaches allow to derive fast algorithms with complexity O(N log N ) regardless of the particle distribution. According to [START_REF] Messner | Fast directional multilevel summation for oscillatory kernels based on Chebyshev interpolation[END_REF], this results in superior performance for a directional FMM on highly non-uniform distributions, a kernel-explicit FMM being better-suited on quasi-uniform particle distributions and on homogeneous surface ones.

The directional FMM algorithms described in the literature exploit low-rank approximations to compress matrices representing the far-field contributions [START_REF] Engquist | Fast directional multilevel algorithms for oscillatory kernels[END_REF][START_REF] Messner | Optimized M2L kernels for the Chebyshev interpolation based fast multipole method[END_REF][START_REF] Messner | Fast directional multilevel summation for oscillatory kernels based on Chebyshev interpolation[END_REF]. However, some highly efficient kernel-independent methods for low-frequency and translation-invariant kernels use Fast Fourier Transforms (FFTs) to efficiently process these matrices in diagonal form [START_REF] Schobert | A multilevel interpolating fast integral solver with fast Fourier Transform acceleration[END_REF]. These FFT-based methods exploit equispaced grids (i.e. cartesian grids or tensorization of the same uniform sampling d times). In FMM formulations based on polynomial interpolation, the use of FFT techniques was shown to be more efficient than low-rank approximations in the low-frequency regime [START_REF] Blanchard | Fast hierarchical algorithms for generating Gaussian random fields[END_REF], especially for high approximation orders. To our knowledge, directional FFT-based FMMs for the high-frequency regime have not been investigated in the past. Moreover, in such polynomial interpolation-based FMMs, the matrices associated to some FMM operators have to be precomputed, but their number drastically increases in the high-frequency regime (see [START_REF] Messner | Optimized M2L kernels for the Chebyshev interpolation based fast multipole method[END_REF]), which strongly increases the precomputation cost. FFT techniques can also help to reduce such precomputation cost.

Contributions

In this paper, we aim at efficiently solving BIE N-body problems with oscillatory kernels in the highfrequency regime. We thus present the first directional interpolation-based FMM taking advantage of polynomial interpolation on equispaced grids to allow fast FFT-based evaluations and precomputations of the far-field matrices. Our FMM retain the overall O (N log N ) complexity of directional kernel-independent FMMs in the high-frequency regime, and the usual linear FMM complexity in the low-frequency regime: we thus rely on the FFTs to lower the computation cost in practice in the high-frequency regime. A consistency proof for the polynomial Lagrange interpolation on such grids and on well-separated sets when dealing with asymptotically smooth kernels is provided. This gives a rigorous explanation of the practical convergence of our approximation process despite the wellknown stability issue of such polynomial interpolation. We also show how to extend the use of some symmetries in the Fourier domain in order to minimize the number of precomputed matrices. FMMs can rely on different 2 d -tree traversals: we adapt here a specific one (namely the dual tree traversal) to directional FMMs and show its relevance, along with a suitable 2 d -tree data structure, for the non-uniform particle distributions typical of BIE problems. We provide a new, publicly available 2 , C++ library implementing our FMM, named as defmm (directional equispaced interpolation-based fmm) and highly optimized on one CPU core. We detail the vectorization process of our near-field direct computation using high-level programming, and we show how to efficiently process the numerous and small FFTs required in our method. We also improve the performance of the BLAS-based operators in the multipole and local expansion translations. Finally, we detail a comparison with a state-of-the-art directional polynomial interpolation-based FMM library (namely dfmm [START_REF] Messner | Optimized M2L kernels for the Chebyshev interpolation based fast multipole method[END_REF][START_REF] Messner | Fast directional multilevel summation for oscillatory kernels based on Chebyshev interpolation[END_REF]), exhibiting important performance gains for defmm in all the test cases.

The paper is organized as follows. In section 2, we first recall the mathematical bases of the FMMs involved in this paper. We then provide in section 3 a consistency proof of our interpolation process on equispaced grids. We detail in section 4 the defmm algorithmic design and we show how to extend and exploit the tree symmetries in the Fourier domain. In section 5, we present various HPC optimizations on one CPU core for defmm, and finally in section 6 we provide numerical results, including a detailed performance comparison.

Presentation of Fast Multipole Methods

We focus in this article on the FMM formulation using polynomial interpolation and its variants [START_REF] Fong | The black-box fast multipole method[END_REF][START_REF] Messner | Fast directional multilevel summation for oscillatory kernels based on Chebyshev interpolation[END_REF][START_REF] Schobert | A multilevel interpolating fast integral solver with fast Fourier Transform acceleration[END_REF], which are here briefly recalled.

Directional FMMs

Suppose that the kernel G can be factorized into

G(x, y) = e iκ|x-y| K(x, y), ( 2
)
where K is an asymptotically smooth non-oscillatory kernel. Directional approaches rely on the directional parabolic separation condition (DPSC) [START_REF] Engquist | Fast directional multilevel algorithms for oscillatory kernels[END_REF], imposing conditions on well-separated subsets of B × B (i.e. pairs of cells in practice) such that the term G u in the expression

G(x, y) = e iκ⟨x,u⟩ e iκ⟨x-y, x-y |x-y| -u⟩ K(x, y) =:Gu(x,y) e -iκ⟨y,u⟩ (3)
does not oscillate, where u ∈ S 2 (referred to as a direction) is a well-chosen approximation of (x -y)/|x -y| on the unit sphere. Since the only term in equation ( 3) depending on both x and y is G u , fast summation methods can be derived. The DPSC is verified when x ∈ t, y ∈ s, where t, s are two cells such that κw 2 ≤ η dist(t, s) and t ⊂ C(s),

with dist(t, s) := min x∈t,y∈s |x -y| the distance between t and s, w the radius3 of s (assuming that w > 1) and C(s) a cone directed by u with apex in the center of s and aperture µ κw [START_REF] Engquist | Fast directional multilevel algorithms for oscillatory kernels[END_REF], η, µ being two strictly positive constants. Such pairs (t, s) are said to be well-separated in the high-frequency regime.

Because the cone aperture in the inequality (4) decreases as the cell radius w increases, more wedges are needed to compute the far field for cells close to the 2 d -tree root than for deeper cells. Since the kernel approximations depend on the cone directions, the number of approximations required for each target cell increases from the 2 d -tree leaves to the root. To lower the far field computation cost, the set of directions at each tree level is chosen so that a nested property is verified [START_REF] Börm | Directional H2-matrix compression for high-frequency problems[END_REF][START_REF] Engquist | Fast directional multilevel algorithms for oscillatory kernels[END_REF]. This leads to a direction tree: each direction at a given 2 d -tree level E is the son of a unique direction at the 2 d -tree level E + 1.

Interpolation-based FMMs

The polynomial interpolation techniques for hierarchical methods [START_REF] Börm | Low-Rank Approximation of Integral Operators by Interpolation[END_REF][START_REF] Fong | The black-box fast multipole method[END_REF][START_REF] Gibermann | Multilevel Approximation of Boundary Integral Operators[END_REF] rely on approximations of G using Lagrange interpolation. We refer to [START_REF] Gasca | Polynomial interpolation in several variables[END_REF] for a description of the multivariate interpolation methods and to the Lagrange interpolation problem. In the following formula, we indicate in gray the terms only appearing when combining the polynomial interpolation techniques for FMM with the directional approach [START_REF] Messner | Optimized M2L kernels for the Chebyshev interpolation based fast multipole method[END_REF][START_REF] Messner | Fast directional multilevel summation for oscillatory kernels based on Chebyshev interpolation[END_REF]. Let Ξ s = {y 1 , ..., y #Ξs } and Ξ t = {x 1 , ..., x #Ξt } (where #Ξ m , m ∈ {t, s}, denotes the cardinal of Ξ m ) be two interpolation grids in the cells s and t respectively, where the pair (t, s) is well-separated, we have:

G(x, y) ≈ e iκ⟨x,u⟩ k S k [t](x) l G u (x k , y l )S l [s](y)e -iκ⟨y,u⟩ ≈ k e iκ⟨x,u⟩ S k [t](x)e -iκ⟨x k ,u⟩ =:S u k [t](x) l G(x k , y l ) e iκ⟨y l ,u⟩ S l [s](y)e -iκ⟨y,u⟩ =:S u l [s](y) , (5)
the polynomials S l [s] and S k [t] verifying for any y ∈ s and any

x ∈ t S l [s](y)      = 1 if y = y l = 0 if y ∈ Ξ s , y ̸ = y l ∈ R otherwise , S k [t](x)      = 1 if x = x k = 0 if x ∈ Ξ t , x ̸ = x k ∈ R otherwise .
Suppose that x ∈ t, t being a target cell. For any source cell s well-separated from t, we define

p t,s (x) := y∈s∩Y G(x, y)q(y) ≈ k S u k [t](x) l G(x k , y l ) y∈s∩Y S u l [s](y)q(y),
where the second line is obtained thanks to the approximation [START_REF] Blanchard | Fast hierarchical algorithms for generating Gaussian random fields[END_REF]. This interpolation process can be repeated on t ′ ∈ Sons(t) and s ′ ∈ Sons(s) recursively, leading in the end to a multilevel algorithm whose operators, named as in [START_REF] Abduljabbar | Extreme scale fmm-accelerated boundary integral equation solver for wave scattering[END_REF][START_REF] Fong | The black-box fast multipole method[END_REF], are the following:

P2M (particles-to-multipole): M u s ′ (y ′ r ) := y∈s∩Y S u r [s ′ ](y)q(y), ∀y ′ r ∈ Ξ s ′ ; M2M (multipole-to-multipole): M v s (y l ) := r S v l [s](y ′ r )M u s ′ (y ′ r ), ∀y l ∈ Ξ s ; M2L (multipole-to-local): L v t (x k ) := l G(x k , y l )M v s (y l ), ∀x k ∈ Ξ t ; L2L (local-to-local): L u t ′ (x ′ h ) := k S v k [t](x ′ h )L v t ′ (x k ), ∀x ′ h ∈ Ξ t ′ ; L2P (local-to-particles): p t,s (x) ≈ h S u h [t ′ ](x)L u t ′ (x ′ h ), ∀x ∈ X ∩ t; where: Ξ s ′ = {y ′ 1 , ..., y ′ #Ξ s ′ } and Ξ t ′ = {x ′ 1 , ..., x ′ #Ξ t ′
} are the interpolation grids in s ′ and t ′ respectively; v is a son of u in the direction tree such that u and v are the best approximations of (ctr(t) -ctr(s))/|ctr(t) -ctr(s)|; ctr(c) denoting the center of the cell c. The direct evaluation of p t,s for all x ∈ X ∩ t corresponds to the P2P (particles-to-particles) operator between t and s. The result of a P2M or M2M operation is a multipole expansion and the result of a M2L or L2L operation is a local expansion. The set of cells s such that p t,s is approximated (i.e. not directly evaluated through a P2P operator) is named the interaction list of t.

All the FMM operators we defined can be interpreted as matrix-vector products, where the vectors are local or multipole expansions. Hence, we shall for example use in the remainder M2L matrices to refer to the matrices involved in the M2L matrix-vector product (and similarly for the other operators).

Equispaced grids

When dealing with multivariate polynomial interpolations, the most natural approach relies on tensorizations of 1D interpolation grids in a cuboid [START_REF] Börm | Low-Rank Approximation of Integral Operators by Interpolation[END_REF][START_REF] Fong | The black-box fast multipole method[END_REF][START_REF] Gibermann | Multilevel Approximation of Boundary Integral Operators[END_REF]. Denoting B ⋆ = [0, 1] d the reference unit box, any cell c writes as

c = α c + β c • B ⋆ (6)
where β c = diam(B)/2 L with L = level(c), and α c ∈ R d . We focus here on equispaced grids, which allow to take advantage of the translation invariance of the kernel, i.e. G((x+z)-(y+z)) = G(x-y).

The interpolation grid Ξ c ⊂ R d in cell c is obtained from a reference grid Ξ ⋆ ⊂ R d of B ⋆ through ( 6)

Ξ c = α c + β c • Ξ ⋆ where Ξ ⋆ := L × • • • × L
where L := {ℓ/L, ℓ = 0, . . . , L}.

These transformations can be used to transfer all evaluations of the kernel G on the reference grid Ξ ⋆ . Take a source cell s and a target cell t, and assume that they belong to the same level of their 2 d -tree and thus admit the same size and β s = β t . Any x ∈ s can be written as x = α s + β s x and any y ∈ t can be expressed as y = α t + β t ŷ, so that we have

G(x -y) = G s,t (x -ŷ) with G s,t (z) := G(α s -α t + β s z). (8)
The function x, y → G(x -y) induces a linear map G(Ξ t , Ξ s ) :

C[Ξ s ] → C[Ξ t ]. Similarly p, q → G s,t (p -q) induces a linear map G s,t (Ξ ⋆ , Ξ ⋆ ) : C[Ξ ⋆ ] → C[Ξ ⋆
] and, by construction, these two linear maps can be identified to the same matrices under consistent orderings

G(Ξ t , Ξ s ) = G s,t (Ξ ⋆ , Ξ ⋆ ).
Rewriting the kernel as above, we explicitly take advantage of the fact that β s = β t . When examining interactions between s and t, the normalized kernel G s,t is evaluated over the difference grid

Ξ ♯ := {p -q, p, q ∈ Ξ ⋆ } ⊂ R d := {ℓ/L, ℓ ∈ N, -L ≤ ℓ ≤ +L} d . (9)
Obviously Ξ ⋆ ⊂ Ξ ♯ , and any vector u ∈ C[Ξ ⋆ ] can be considered as an element of C[Ξ ♯ ] by means of extension by zero: define χ(u) ∈ C[Ξ ♯ ] by χ(u)(p) := u(p) if p ∈ Ξ ⋆ and χ(u)(p) := 0 otherwise. The linear map χ simply reduces to a matrix with entries equal to 0 or 1, with at most one nonvanishing entry per row, so applying such a matrix to a vector is computationally harmless. The application of this matrix numerically corresponds to a zero-padding.

We have card(Ξ ♯ ) = T d with T := (2L + 1), i.e. Ξ ♯ contains T samples in each direction. The optimization we wish to exploit is based on Fast Fourier Transforms (FFTs), this is why we embed the action of (8) in a periodic setting. The kernel G s,t can be extended by T /L-periodicity, denoting

G ♯ s,t the unique element of C[Ξ ♯ + Z d ] satisfying G ♯ s,t (p + (T /L)τ ) = G s,t (p) for all p ∈ Ξ ♯ , τ ∈ Z d . (10) In particular G ♯ s,t (p) = G s,t (p) for p ∈ Ξ ♯ . Denote G ♯ s,t (Ξ ♯ , Ξ ♯ ) the matrix associated to the mapping from C[Ξ ♯ ] → C[Ξ ♯ ] associated to p, q → G ♯ s,t (p -q)
. This matrix admits a circulant form due to the periodicity of the corresponding kernel. The notations introduced above lead to a factorization of M2L interactions between cells s and t,

G(Ξ s , Ξ t ) = χ ⊤ • G ♯ s,t (Ξ ♯ , Ξ ♯ ) • χ. (11)
As we are considering tensorized interpolation grids, and since G ♯ s,t (Ξ ♯ , Ξ ♯ ) admits circulant form, we naturally make use of a multi-dimensional Discrete Fourier Transform (DFT) to process this matrix. Let us briefly explain how. Given the dual grid

Ξ ♯ = L • Ξ ♯ ⊂ R d , the d-dimensional DFT refers to the linear maps F : C[Ξ ♯ ] → C[ Ξ ♯ ] defined by the formula F(u)(ξ) := 1 T d/2 x∈Ξ ♯ u(x) exp(-2ıπ ξ • x) ξ ∈ Ξ ♯ . ( 12
)
The inverse DFT is simply the adjoint map defined by the formula

F * (û)(x) = T -d/2 ξ∈ Ξ ♯ û(ξ) exp(+2ıπ ξ• x) for x ∈ Ξ ♯ . Due to the circulant form of the matrix G ♯ s,t (Ξ ♯ , Ξ ♯ ), there exists a diagonal matrix D s,t such that G ♯ s,t (Ξ ♯ , Ξ ♯ ) = F * • D s,t • F hence G(Ξ s , Ξ t ) = χ ⊤ • F * • D s,t • F • χ. ( 13
)
The diagonal of D s,t = diag ξ∈ Ξ ♯ ( F(G s,t )(ξ) ) contains the Fourier coefficients of the kernel. From the factorized form above, we see that matrix-vector product for the M2L matrix G(Ξ s , Ξ t ) can be processed efficiently. Indeed the action of χ is inexpensive, and F can be applied efficiently with O(d L d log L) complexity by means of FFTs. Notice that the diagonal of D s,t (i.e. its non-zero entries) is obtained by DFT of the first column of G ♯ s,t (Ξ ♯ , Ξ ♯ ) and can also be efficiently processed using FFT.

Since the same Fourier basis is used for all M2L matrices, one can apply F and F * (as well as χ and χ ⊤ ) once to each multipole and local expansions respectively before and after all the M2L evaluations. This reduces M2L matrix evaluations to products with diagonal matrices D s,t (i.e. Hadamard products), each with a complexity of (2L -1) d flop. In comparison with the product of a vector by a rank-k matrix approximation of a M2L matrix, requiring 2kL d operations, the product with the diagonalization of the circulant embedding of this matrix is theoretically less costly if k ≥ 2 d-1 , which is valid except when requesting very low overall accuracies. However, low-rank approaches usually rely on highly optimized implementations (namely BLAS routines) to significantly reduce the computation times of the dense matrix-vector (or matrix-matrix) products (see e.g. [START_REF] Messner | Optimized M2L kernels for the Chebyshev interpolation based fast multipole method[END_REF]). Therefore, theoretical operation counts are insufficient to compare a FFT-based approach with a low-rank one. We will have to compare practical time measurements using implementations of both approaches to determine the best one. One can already notice that, for non-oscillatory kernels, FFT-based approaches offer better or similar performance (depending on the interpolation order) than low-rank compressions for the M2L matrix evaluations [START_REF] Blanchard | Fast hierarchical algorithms for generating Gaussian random fields[END_REF][START_REF] Chen | Fast algorithms for evaluating the stress field of dislocation lines in anisotropic elastic media[END_REF]: this encourages us to investigate the FFT-based approaches for the oscillatory kernels.

In this purpose, we recall first that the interpolation process on equispaced grids is known to be subject to the Runge phenomenon, which may cause the process to diverge, especially for high interpolation orders (see for instance [START_REF] Chen | Fast algorithms for evaluating the stress field of dislocation lines in anisotropic elastic media[END_REF]).

Consistency of the equispaced interpolation

Despite the asymptotic instability of the Lagrange polynomial interpolation process on equispaced grids, we show in this section that the approximation (5) actually converges on well-separated sets.

We define the Multipole Acceptance Criterion (or MAC) A as a boolean function that takes two cells in argument and returns 1 if and only if these two cells are well-separated. As in [START_REF] Blanchard | Fast hierarchical algorithms for generating Gaussian random fields[END_REF][START_REF] Fong | The black-box fast multipole method[END_REF], we only consider cells at the same 2 d -tree level as input parameters of the MAC. These cells t, s thus have the same radius a, i.e. they are translations of [-a, a] d in R d (see section 2.3). The interpolation grids in t, s can be expressed relatively to the cell centers, so the polynomial interpolation of G on t × s can be seen as an interpolation of G on [-a, a] 2d using G(x, y) := G(ctr(t) + x, ctr(s) + y).

Here, G depends on t × s.

We now introduce and prove the following theorem, regarding the interpolation process consistency on our equispaced grids.

Theorem 1. Suppose that the condition A(t, s) = 1 implies that G is analytic in each variable at any point in [-a, a] with a convergence radius R such that R > 2a e . Then the Lagrange interpolation of G on t × s using L equispaced points in each variable, denoted by

I t×s L [ G], verifies lim L→+∞ I t×s L [ G] - G L ∞ ([-a,a] 2d ) = 0. More specifically, if radius(s) + radius(t) ≤ dist(t, s) and if G is asymptotically smooth, the convergence is geometric, i.e. ∃ ν ∈ (0, 1) s.t. ||I t×s L [ G] -G|| L ∞([-a,a] 2d ) = O ν L .

Preliminary results

The idea of the proof of Theorem 1 consists in combining 1D convergence estimates with the interpolation error results provided in [START_REF] Mößner | Error bounds for polynomial tensor product interpolation[END_REF] (restricted to our asymptotically smooth kernel functions). Our convergence proof differs from the one of [START_REF] Fong | The black-box fast multipole method[END_REF][START_REF] Messner | Fast directional multilevel summation for oscillatory kernels based on Chebyshev interpolation[END_REF] because we cannot rely on the properties of the Chebyshev grids. We denote by C ∞ (Ω) the space of multivariate functions f with bounded derivatives

∂ β f := ∂ β1 f...∂ β 2 d f on Ω, ∀β ∈ N 2d such that max k=1,...,2d β k ≤ L for any domain Ω ⊂ R 2d . We denote by A := {α = (α 1 , ..., α 2d ) ∈ {0, 1} 2d , ∃j, α j = 1}. We also use the notation ||f || ∞ := sup z∈[-a,a] 2d |f ((z)|. Theorem 2. ( [26] Thm. 2.1) For f ∈ C ∞ ([-a, a] 2d ), the product interpolation I [-a,a] 2d L [f ] of f in [-a,
a] 2d with the same 1D rule with L interpolation nodes in each variable verifies

f -I [-a,a] 2d L [f ] ∞ ≤ α∈A ω ᾱ L ∂ αL f ∞ where αL = (α 1 L, ..., α 2d L), ᾱ := 2d k=1 α k and ω L := 1 L! L-1 k=0 (•-x k ) ∞ , x k being the k th interpolation point of the 1D rule.
For equispaced grids, the constant ω L can be bounded using the following lemma.

Lemma 1. Let {x k := -a+2ak/(L-1) | k ∈ [[0, L-1]]} ⊂ [-a, a] an interpolation grid of equispaced points. We have ω L ≤ 2a L-1 L /(4L). Proof. Let x ∈ [-a, a]. We have L-1 j=0 (x -x j ) = (2a) L L-1 j=0 x + a 2a - j (L -1) . Let y := x+a 2a ∈ [0, 1]. This leads to L-1 j=0 (x -x j ) ≤ 2a L -1 L L-1 j=0 y(L -1) -j . Because y(L -1) ∈ [0, (L -1)],
L-1 j=0 y(L -1) -j is maximal for y ∈ (0, 1/(L -1)). Using a simple recurrence, one may easily show that

L-1 j=0 y(L -1) -j ≤ (L -1)! 4 , which implies that ω L ≤ 1 L! 2a L-1 L (L -1)!/4 ≤ 2a L-1 L /(4L).
We now want to bound the partial derivatives of the interpolated function. This is the purpose of the following lemma.

Lemma 2. If f is analytic in all its variables at any point of [-a, a] with a convergence radius R > 0, we have

||∂ αL f || ∞ ≤ C r 2d L! r L ᾱ with 0 < r < R, ᾱ := 2d k=1 α k , C ∈ R * +
being a constant independent of L, and α being defined as in Thm. 2.

Proof. Since f is analytic in all its variables, we can apply the Cauchy integral formula (see [START_REF] Ahlfors | Complex analysis: An introduction to the theory of analytic functions of one complex variable[END_REF], Chapter 2, Theorem 6), allowing us to write:

f (p) = 1 2πi 2d Γ1 ... Γ 2d f (z) (z 1 -p 1 )...(z 2d -p 2d ) dz 1 ...dz 2d ,
where

Γ j := {ξ ∈ C | dist(ξ, [-a, a]) = r}.
We thus have:

∂ αL p f (p) ∞ ≤ (2π) -2d Γ1
...

Γ 2d f (z) ∂ αL p 1 (z 1 -p 1 )...(z 2d -p 2d ) ∞ dz 1 ...dz 2d ≤ |Γ 1 |...|Γ 2d | (2π) sup z∈Γ1×...×Γ 2d f (z) ∂ αL p 1 (z 1 -p 1 )...(z 2d -p 2d ) ∞ .
The term sup z∈Γ1×...×Γ 2d |f (z)| is bounded thanks to the analyticity of f on a neighborhood of [-a, a] 2d encompassing the Γ j 's (with a convergence radius equal to R). Hence, there exists a constant M (Γ 1 , ..., Γ 2d ) ∈ R + * such that, defining

C(Γ 1 , ..., Γ 2d ) := (2π) -2d |Γ 1 |...|Γ 2d |M (Γ 1 , ..., Γ 2d ),
where

|Γ j | denotes the length of the path Γ j , j ∈ [[1, 2d]], with ᾱ := 2d k=1
α k , we have

∂ αL f (p) ∞ ≤ C(Γ 1 , ..., Γ 2d ) ∂ αL 1 (z 1 -p 1 )...(z 2d -p 2d ) ∞ (Since α k = 0, 1) ≤ C(Γ 1 , ..., Γ 2d ) 2d k=1 L! α k (z k -p k ) -Lα k -1 ∞ ≤ C(Γ 1 , ..., Γ 2d ) 2d k=1 L! α k r -Lα k -1 =( L! r L ) ᾱ/(r 2d ) .

Proof of the main theorem

We can now prove Thm. 1.

Proof. Following Lem. 1 and Thm. 2, we have,

I t×s L [ G] -G ∞ ≤ α∈A 2d k=1 ω α k L ∂ αL G ∞ ≤ α∈A 2d k=1 2a L -1 L /(4L) α k ∂ αL G ∞ which becomes, thanks to Lem. 2 I t×s L [ G] -G ∞ ≤ C α∈A 2d k=1 2a L -1 L L! 4r L L α k r -1 .
Now, by applying Stirling's inequality L! ≤ e -(L-1) L L+1/2 , one obtains

I t×s L [ G] -G ∞ ≤ C α∈A 2d k=1 2a L -1 L e -(L-1) L L+1/2 4r L L α k r -1 ≤ C α∈A 2d k=1 2a L -1 L eL L+1/2 4(er) L L α k r -1 ≤ C α∈A 2d k=1 2a re L L L -1 L e 4 √ L α k r -1 .
For L ≥ 2, we have

L L-1 L
≤ 4, which allows to write

I t×s L [ G] -G L ∞ ([-a,a] 2d ) ≤ C α∈A 2d k=1 2a re L e √ L α k r -1 ≤ Cr -2d α∈A 2d k=1 2a re L e √ L α k =( 2a re ) ᾱL e √ L ᾱ using ᾱ = 2d k=1 α k .
The number of terms in the sum over A is finite and depends only on the dimension: there is indeed 2 2d -1 terms in this sum since α = (0, ..., 0) does not verify ||α|| ∞ = 1. This estimate thus tends to zero when

L tends to infinity if 2a re < 1. Since r < R, this is verified if 2a e < R. Indeed, e √ L ᾱ ≤ e √ 2 2d
since we assumed that L ≥ 2. In addition, each 2a re ᾱL can be bounded by 2a re L since 2a re < 1. We finally have

I t×s L [ G] -G L ∞ ([-a,a] 2d ) ≤ Cr -2d (2 2d -1) e √ 2 2d Does not depend on L 2a re L → L→+∞ 0 . (14)
Now, let t and s be two cells with radius a at the same 2 d -tree level such that dist(t, s) ≥ 2a = radius(t) + radius(s). For asymptotically smooth kernel G : R d × R d → C, G is analytic with convergence radius R = dist(t, s)/2 (only singular at x = y). Hence, R ≥ a. We can thus choose in inequality ( 14) r ∈ 2R e , R . Such a r exists since 2 e < 1. We thus have

2 e dist(t, s) 2 < r < dist(t, s) 2 ⇒ radius(t) + radius(s) dist(t, s) =:θ > 2a re > θ 2 e .
Finally, because θ ≤ 1, we conclude from inequality ( 14) that for any ν ∈ 2 e θ, θ

I t×s L [ G] -G L ∞ ([-a,a] 2d ) = O ν L . ( 15
)
This proof has a geometric interpretation. In the inequality ( 14), the term 2a re L somehow corresponds to a MAC such as in [START_REF] Dehnen | A hierarchical O(N) force calculation algorithm[END_REF]: a refers to the radius of an interacting cell and r is related to the distance between the interacting cells. The greater this distance, the greater r and the better this estimate. Another information we obtain from the inequality ( 14) is that the convergence should be geometric in the 1D interpolation order. Indeed, considering θ := radius(t)+radius(s) dist(t,s)

, Estimate [START_REF] Engquist | Fast directional multilevel algorithms for oscillatory kernels[END_REF] somehow indicates that the error of the interpolation process is O (cθ) L with c a constant smaller than 1. In practice, our new directional MAC (defined in section 4.3.2 below), as well as the MAC from [START_REF] Fong | The black-box fast multipole method[END_REF][START_REF] Messner | Optimized M2L kernels for the Chebyshev interpolation based fast multipole method[END_REF][START_REF] Messner | Fast directional multilevel summation for oscillatory kernels based on Chebyshev interpolation[END_REF] (referred to as the strict MAC4 in this article), both verify the assumptions of theorem 1: the interpolation process of G therefore converges on such well-separated sets. In Fig. 2, we provide numerical investigation of the relative error for our equispaced interpolation with respect to the interpolation order in each variable on well-separated sets for two kernels matching the scope of Thm. 1 (namely asymptotically smooth kernels) and for dimensions d = 1, 2, 3 (similar results were obtained for higher dimensions but are not reported here since they are less relevant for our application cases). To this end, we considered well separated sets X ⊂ [0, 1] d and Y ⊂ [2, 3] d with 30000 randomly sampled particles in each, and randomly chosen charges at each particle in Y . Then, we computed the maximum relative error between a direct N -body problem on these sets and the approximated one, based on an equispaced interpolation on [0, 1] d × [2, 3] d . These results numerically confirm, for any dimension, the expected geometric convergence of the equispaced interpolation applied on well-separated sets. Similar results (not shown here) were also obtained for other configurations of well-separated cells.

We conclude this section by two important remarks. First, the practical efficiency of the interpolation-based FMM relies on the implicit assumption that the constant C(Γ 1 , ..., Γ 2d ) is small. This is the case when using the strict MAC in the low-frequency regime (a justification can be found in [START_REF] Messner | Fast directional multilevel summation for oscillatory kernels based on Chebyshev interpolation[END_REF]) and our new directional MAC (see section 4.3.2) in the high-frequency one. Second, the ill-conditioning of such interpolation on equispaced grids may cause an exponential amplification of floating-point rounding errors at the boundaries of the interpolation domain (see [START_REF] Platte | Impossibility of fast stable approximation of analytic functions from equispaced samples[END_REF][START_REF] Smith | Lebesgue constants in polynomial interpolation[END_REF]). We thus cannot expect a numerical convergence for any order in practice. Nevertheless, according to our tests (not shown here), we can reach a relative error of 10 -12 on 3D particle distributions with double-precision arithmetic before facing numerical instabilities. Practical applications usually require much less accurate approximations. Figure 2: Relative error in maximum norm between matrix-vector product on well-separated cells using direct computation or using single equispaced interpolation with respect to the interpolation order in each variable. d denotes the dimension, Laplacian denotes the kernel |x -y| -1 , Gaussian denotes the kernel e -|x-y| 2 .

defmm: a directional equispaced interpolation-based FMM

In this section, we present our new defmm library, integrating polynomial interpolations on equispaced grids within a directional kernel-independent FMM.

Overall algorithm

Like directional FMM algorithms [START_REF] Engquist | Fast directional multilevel algorithms for oscillatory kernels[END_REF][START_REF] Messner | Fast directional multilevel summation for oscillatory kernels based on Chebyshev interpolation[END_REF], our FMM is divided into four mains steps:

1. a precomputation step that involves the 2 d -tree constructions and the precomputation of M2L matrices, 2. an upward pass consisting in applying the P2M and M2M operators along the source tree in order to compute the (possibly directional) multipole expansions, 3. a horizontal pass consisting in applying the M2L and P2P operators in order to transform (possibly directional) source multipole expansions into (possibly directional) target local expansions on well-separated pairs of cells, 4. a downward pass consisting in applying the L2L operators along the target tree in order to compute the leaf (possibly directional) local expansions that can be transformed into far-field contributions by means of L2P operators.

In the following sections, we detail how the efficient handling of equispaced interpolations as well as our algorithmic choices impact these different steps. Hence, in Sect. 4.2 we introduce new operators

P M 0 M0 M 1 M1 M 2 M2 P L0 L 0 L1 L 1 L2 L 2 P2P P2M M2M hfM2M L2P L2L hfL2L M2L (Fourier) M2L (Fourier) hfM2L (Fourier) M2F F2L M2F F2L hfM2F hfF2L hfP2M hfL2P
Low frequency regime High frequency regime that slightly modify the definition of the upward pass (step 2) and of the downward pass (step 4). In Sect. 4.3 we present how to efficiently perform the horizontal pass (step 3), as well as our strategy to minimize the number of precomputed M2L matrices (step 1). Finally, we show in Sect. 4.4 how the use of symmetries can further reduce the precomputation cost and how these techniques combine with the horizontal pass. The new final algorithm is provided in Sect. 4.5.

Directional FFT-based FMM

The integration of FFT techniques exploiting equispaced grids in a directional polynomial interpolationbased FMM requires operations between various directional expansions to be converted in the Fourier domain.

M2F and F2L operators

As described in section 2, the M2L application on equispaced grids involves three steps: the conversion of the involved multipole expansion (extended by zero-padding) into the Fourier domain, the application of the diagonal M2L operator in this domain, and then the conversion of the resulting (extended) local expansion in the Fourier domain into a local expansion. We thus introduce two extra operators, the M2F (multipole-to-Fourier) and F2L (Fourier-to-local) ones, applying the matrices Fχ to multipole expansions in the real domain and χ T F * to local expansions in the Fourier domain. Let T be a M2L matrix, we have (see section 2.3)

T = χ ⊤ F * F2L D s,t Diagonal M2L
(Fχ)

M2F

.

In Figure 3, we depict how the different operators are applied in defmm. The operators are non-directional (i.e. without the terms in grey in section 2.2) in the low-frequency regime, and directional in the high-frequency one.

Direction generation

In the directional FMM presented in [START_REF] Engquist | Fast directional multilevel algorithms for oscillatory kernels[END_REF] the symmetries within the sets of directions at each tree level are exploited to reduce the number of precomputed M2L matrices. This however constrains the way the directions are generated. On the contrary, using a directional interpolation-based approach the M2L matrices do not depend on the directions (see section 2.2). We thus exploit a similar construction algorithm as presented in [START_REF] Börm | Directional H2-matrix compression for high-frequency problems[END_REF][START_REF] Engquist | Fast directional multilevel algorithms for oscillatory kernels[END_REF]: starting from a given regular solid, the set of projected face centers on the unit sphere are chosen as directions at the first high-frequency level. Then, each face is subdivided into 2 d-1 other faces and the process is repeated recursively in order to obtain the direction tree. Engquist & Ying [START_REF] Engquist | Fast directional multilevel algorithms for oscillatory kernels[END_REF] relies on the d-cube for the direction generation: according to our tests, the d-cube provides indeed the best compromise between the number of directions and the obtained accuracy among the platonic solids. Since, a d-cube has d! faces, there are d!2 E(d-1) directions at the E th high-frequency level in defmm.

Algorithmic design

4.3.1 2 d -tree construction 2 d -trees can be built either by setting the maximum number of particles per leaf (denoted as N crit), or by setting the maximum possible level (denoted as M axDepth). For (highly) non-uniform particle distributions, the M axDepth strategy leads to (numerous) empty cells, which should not be stored and processed in practice: this requires advanced data structures (see e.g. [START_REF] Coulaud | High-performance BLAS formulation of the adaptive Fast Multipole Method[END_REF][START_REF] Hariharan | A Scalable Parallel Fast Multipole Method for Analysis of Scattering from Perfect Electrically Conducting Surfaces[END_REF][START_REF] Nabors | Preconditioned, Adaptive, Multipole-Accelerated Iterative Methods for Three-Dimensional First-Kind Integral Equations of Potential Theory[END_REF]). The M axDepth strategy may also generate cells with a greatly varying number of particles [START_REF] Pavlov | Porting and Verification of ExaFMM Library in MIC Architecture[END_REF]. Moreover, a comparison of hierarchical methods in [START_REF] Fortin | Comparisons of different codes for galactic N-body simulations[END_REF] has detailed how the N crit strategy better adapts to highly non-uniform distributions than the M axDepth one. Considering the highly non-uniform distributions used in BIE problems, we therefore choose to rely on the N crit strategy in defmm. This differs from dfmm [START_REF] Messner | Optimized M2L kernels for the Chebyshev interpolation based fast multipole method[END_REF][START_REF] Messner | Fast directional multilevel summation for oscillatory kernels based on Chebyshev interpolation[END_REF] which relies on the M axDepth strategy.

Concretely, our 2 d -tree construction algorithm is similar to the exafmm5 one: the particles are first sorted according to the Morton ordering, then each cell is built with a pointer on its first particle. The particles within a cell are thus stored consecutively in memory and followed by the particles of the next cell in the Morton ordering. This ensures that the data is localized for the P2P evaluations. In addition, since the particles of the sons of a given cell are stored consecutively, we can directly apply P2P operators on non-leaf cells. The charges and potentials are stored in two dedicated arrays whose entries correspond to the sorted particles for the same data-locality reason.

We also store all the directional expansions associated to a given cell c in the same array to enhance data locality during the translations of local and multipole expansions: the M2M and L2L evaluations regarding c are indeed performed alltogether (see section 5.2 for the optimizations of these steps).

Dual Tree Traversal

Using the N crit criterion however complexifies the interaction list structure (see lists U, V, W and X in e.g. [START_REF] Ying | A kernel-independent adaptive fast multipole algorithm in two and three dimensions[END_REF]). In directional FMMs, the DPSC for well-separateness (see equation ( 4)) further complicates the construction of the interaction lists (multiple lists for each target cell with varying sizes and shapes depending on the tree level). We therefore adapt here the Dual Tree Traversal (DTT) [START_REF] Dehnen | A hierarchical O(N) force calculation algorithm[END_REF][START_REF] Yokota | An FMM based on dual tree traversal for many-core architectures[END_REF] to directional FMMs. The DTT is a simple and recursive procedure that simultaneously traverses the target and source trees to evaluate on-the-fly the M2L and P2P operators by dynamically testing the MAC on pairs of target and source cells. Hence, the interaction lists are Algorithm 1 Dual Tree Traversal (DTT) between target cell t and source cell s end for 15: end if traversed and evaluated implicitly but never built. This is another difference with dfmm which relies on a DTT-like algorithm to explicitly build the interaction lists during a precomputation step, and then separately processes these lists for each target cell.

Since FFTs can only efficiently accelerate M2L operations between cells at the same level (because of the relatively small involved circulant embedding), our DTT differs from the original one [START_REF] Dehnen | A hierarchical O(N) force calculation algorithm[END_REF][START_REF] Yokota | An FMM based on dual tree traversal for many-core architectures[END_REF] as follows: (i) when the MAC fails, both cells are split (not the largest one); (ii) a P2P operation is required as soon as one of the two cells is a leaf (not both).

We also aim at one single DTT algorithm for both the low-and high-frequency regimes. This requires the MAC to depend only on the cell radii and distances, and not on directions. Fortunately, this is possible in a directional polynomial interpolation-based FMM, since the M2L matrices do not depend on the directions. Our DTT hence relies in the low-frequency regime on the strict MAC (see section 3.2), which is simple and allows the most far-field interactions for each target cell, and on the following new directional MAC in the high-frequency regime (for wavenumber κ) between cells t, s with radius w at the same tree level:

max{κw 2 , 2w} dist(t, s) ≤ η, (16)
where η > 0 (η = 1 in practice). This new directional MAC is similar to another directional MAC [START_REF] Börm | Directional H2-matrix compression for high-frequency problems[END_REF], but contrary to this latter our MAC does not depend on the directed wedge condition: we do not consider directions in this MAC but only ratios between cell radii and the cell distance. This means that the DTT performs in the high-frequency regime in a similar way than in the low-frequency one, without considering the directional aspects that are entirely hidden in the M2L applications. We emphasize that the interpolation process using equispaced grids is consistant according to section 3 on cells complying with these two MACs. When this directional MAC is satisfied, one then has to retrieve the best direction to select the directional multipole expansion (from source s) and the directional local expansion (from target t) which are relevant for the M2L operation between t and s. The search of the best direction may have a non-negligible cost, but there is one single direction associated to each M2L matrix. We thus precompute these best directions during the M2L matrix precomputation step (see section 4.3.3).

Algorithm 2 Blank Dual Tree Traversal (BDTT) between target cell t and source cell s end if 18: end if

In the end, we obtain the simple DTT algorithm presented in Algorithm 1 to differentiate in defmm the far-field M2L operations from the near-field P2P ones, in both low-and high-frequency regimes.

One may notice that list-based approaches (i.e. without DTT) allow to group or "stack" multiple M2L operations (written as matrix-vector products), into matrix-matrix products (see e.g. [START_REF] Coulaud | High-performance BLAS formulation of the adaptive Fast Multipole Method[END_REF][START_REF] Malhotra | Algorithm 967: A Distributed-Memory Fast Multipole Method for Volume Potentials[END_REF][START_REF] Messner | Optimized M2L kernels for the Chebyshev interpolation based fast multipole method[END_REF]). This enables indeed to benefit from the higher efficiency of level-3 BLAS routines, especially for uniform distributions. We refer to this technique as vector stacking. Here however, we have M2L operations corresponding to Hadamard products whose grouping cannot lead to more efficient level-3 BLAS operations. Such Hadamard products may be converted to matrix-vector products as shown in [START_REF] Malhotra | Algorithm 967: A Distributed-Memory Fast Multipole Method for Volume Potentials[END_REF], but this requires extra zeros for non-uniform particle distributions. Considering the highly non-uniform distributions typical of BIE problems, we believe that the potential gain would be too limited. We hence do not consider vector stacking for our Hadamard products, and we rather rely on the DTT to efficiently process the BIE non-uniform distributions.

Blank passes

On non-uniform distributions only a subset of all possible directional expansions will be used in practice due to the DPSC. Hence, defmm determines the actually needed directional expansions (named effective expansions) during a precomputation step and only computes and stores these effective expansions to save memory and computations. In this purpose, defmm performs first a blank DTT (see Algorithm 2) to compute each required M2L matrix and to mark all M2L-interacting cells with the corresponding direction (see line 5 in Algorithm 2). Then, thanks to a blank downward pass (see algorithm 3), the required directions are propagated down to the leaves of the 2 d -trees. end for 10: end if

Exploiting symmetries in the Fourier domain

To minimize the number of M2L matrices to be precomputed, in the case of centered expansions in 2 d -trees, it has been observed in the literature that symmetries can be used (see for instance [START_REF] Engquist | Fast directional multilevel algorithms for oscillatory kernels[END_REF][START_REF] Messner | Optimized M2L kernels for the Chebyshev interpolation based fast multipole method[END_REF]). The underlying symmetry group is actually the hyperoctahedral one [START_REF] Chollet | Symmetries and Fast Multipole Methods for Oscillatory Kernels[END_REF] (i.e. the octahedral group in 3D, corresponding to the symmetry group of the cube), denoted by D d in dimension d, that can be realized as a group of rotation matrices. The use of symmetries strongly reduces the precomputation cost, especially in the high-frequency regime. Hence, we want to also exploit these symmetries in defmm. To do so, one has to express these symmetries in the Fourier domain, since the (modified diagonal) M2L matrices are expressed in this domain.

Taking symmetries into account

Let us first detail how, in a given cell c of the (source or target) tree, a permutation matrix over the interpolation grid Ξ c can be associated to a symmetry. Take any rotation R : R d → R d that leaves the unit cube centered at 0 invariant. The set of rotations satisfying this property forms the so-called hyperoctahedral group. The translated grid Ξ c -ctr(c) is centered at the origin, so Ξ c -ctr(c) = R(Ξ c -ctr(c)) i.e. Ξ c = ctr(c) + R(Ξ c -ctr(c)). This transformation, represented in Figure 4, induces a linear map R c :

C[Ξ c ] → C[Ξ c ] defined by R c (u)(p) := u(ctr(c) + R(p -ctr(c))) ∀u ∈ C[Ξ c ], ∀p ∈ Ξ c . (17)
This matrix R c is a permutation of interpolation nodes, and the correspondence R ↔ R c implements a permutation representation of the hyperoctahedral group. This is the group representation that we propose to exploit to minimize the number of precomputed M2L matrices.

Denote

R ⋆ : C[Ξ ⋆ ] → C[Ξ ⋆ ]
the permutation induced on Ξ ⋆ by means of the correspondence described above. When the same rotation transformation is applied to both source and target cells s and t, the transformed M2L matrix writes

R ⊤ t G(Ξ t , Ξ s )R s = R ⊤ ⋆ G s,t (Ξ ⋆ , Ξ ⋆ )R ⋆ = R ⊤ ⋆ χ ⊤ F * D s,t FχR ⋆ . (18)
We are going to show that the product above actually writes χ ⊤ F * Ds,t Fχ where Ds,t is a diagonal

matrix mapping C[ Ξ ⋆ ] → C[ Ξ ⋆ ]
with coefficients obtained by permuting the diagonal of D s,t .

Figure 4: Left: Source cell s in the low-frequency regime (blue) with possible target cells t (red) such that s is in the interaction list of t. The corresponding M2L matrices for cells t in pale red can be all deduced (for instance) from permutations of the dark red ones. Symmetry axes of the square (2-cube) are represented in orange. Middle and right: M2L matrix between the red and blue well-separated cells (middle) as a rotation of the M2L matrix between the red and grey cells (right), that permutes the relative positions of the interpolation nodes (colored dots). Since the entries of the M2L matrices correspond to interpolation nodes, permutations matrices represent this process.

The underlying symmetry of D 2 is the reflection with regard to the line x = y.

In concrete computational terms, this means that, whenever two M2L interactions only differ by a rotation, one M2L matrix is deduced from the other by a simple permutation of the diagonal entries of the factor term D s,t , which reduces storage and precomputation cost.

Let c ⋆ = (1/2, . . . , 1/2) = ctr(Ξ ⋆ ) ∈ R d refer to the center of the normalized interpolation grid. Denote Ξ 0 = Ξ ⋆ -ctr(Ξ ⋆ ) the normalized interpolation grid translated so as to be centered at 0. Take two vectors u, v ∈ C[Ξ ⋆ ]. Then we have

v ⊤ R ⊤ ⋆ G s,t (Ξ ⋆ , Ξ ⋆ )R ⋆ u = p∈Ξ⋆ q∈Ξ⋆ G s,t (p -q)u(c ⋆ + R(p -c ⋆ ))v(c ⋆ + R(q -c ⋆ )) = x∈Ξ0 y∈Ξ0 G s,t (x -y)u(c ⋆ + R(x))v(c ⋆ + R(y)) = x ′ ∈Ξ⋆ y ′ ∈Ξ⋆ G s,t (R * (p -q))u(p)v(q). (19) Since u, v are arbitrarily chosen in C[Ξ ⋆ ], this can be rewritten in condensed form by R ⊤ ⋆ G s,t (Ξ ⋆ , Ξ ⋆ )R ⋆ = (G s,t • R * )(Ξ ⋆ , Ξ ⋆ )
. Now there only remains to return to the calculus presented in section 2.

3 that shows R ⊤ ⋆ G s,t (Ξ ⋆ , Ξ ⋆ )R ⋆ = χ ⊤ F * D[G s,t • R * ] Fχ with D[G s,t • R * ] = diag( F(G s,t • R * ) ). ( 20
)
To summarize, the only difference between R ⊤ ⋆ G s,t (Ξ ⋆ , Ξ ⋆ )R ⋆ and G s,t (Ξ ⋆ , Ξ ⋆ ) lies in the coefficients of the Fourier symbol in the central diagonal term of their factorized form. Let us examine how the rotation R * acts on the Fourier symbol of G s,t . According to [START_REF] Chollet | Symmetries and Fast Multipole Methods for Oscillatory Kernels[END_REF], and since R(Ξ ♯ ) = Ξ ♯ , we have

F(G s,t • R * )(ξ) = 1 T d/2 x∈Ξ ♯ G s,t (R * x) exp(-2ıπ ξ • x) = 1 T d/2 y∈Ξ ♯ G s,t (y) exp(-2ıπ R * (ξ) • y) = F(G s,t )(R * ξ) (21)
which again summarizes as F(G s,t • R * ) = F(G s,t ) • R * . In other words, the permutation on the grid Ξ ♯ that should be applied on the symbol of the periodized M2L operator is the permutation associated to the inverse rotation R * = R -1 .

Symmetries in practice

Usually, the multipole and local expansions are permuted respectively before and after the evaluation of a M2L operator (see e.g. [START_REF] Malhotra | Algorithm 967: A Distributed-Memory Fast Multipole Method for Volume Potentials[END_REF][START_REF] Messner | Optimized M2L kernels for the Chebyshev interpolation based fast multipole method[END_REF]). Here we can reduce this number of applied permutations to only one thanks to the diagonal form of the M2L matrices in the Fourier domain in defmm (see equation ( 20)). Indeed, the permutations can all be applied to these diagonal M2L matrices, resulting in a permutation of their diagonal entries. This reduces the number of permutations from two to only one.

In practice, the permutations induce memory indirections when performing the Hadamard products that may prevent the compiler auto-vectorization. We hence rely on an OpenMP6 directive to enforce the compiler vectorization (see [START_REF] Chollet | Symmetries and Fast Multipole Methods for Oscillatory Kernels[END_REF] for details).

Final algorithm

We here summarize the main steps of our final algorithm, along with our improvements. Following Sect. 2.2, we denote in grey the terms which are specific to the high-frequency regime. This high-level algorithm may be easily implemented using recursion for the different passes. We refer the interested reader to the public version of our code.

Optimizations

We here present algorithmic and programming optimizations regarding high performance computing on one CPU core.

Processing the FFTs

We rely on the state-of-the-art FFTW library [START_REF] Frigo | The design and implementation of FFTW3[END_REF] to efficiently perform all at once our numerous FFTs. Our 3D FFTs (considering d = 3) are however small: (2L -1) d elements, with usually L ≤ 7. This makes the FFTW "plan" creation, required to decide the most efficient FFT implementation, more time consuming than its execution (i.e. the actual FFT computation). Moreover, zero-padding is required for the circulant embedding (see section 4.4) of the expansions stored in the Fourier domain: this implies expansion copies before each FFT, as well as after each reverse FFT. One could first use the "batch" FFTW feature 7 to efficiently process our numerous small FFTs. This performs all FFTs with one single FFTW plan and can improve performance compared to multiple individual FFT calls. All expansion copies hence have to be performed all together before and after the batch FFTW call. We refer to this method as the batch one. Since all our expansions have the same size and the same memory alignment, we can also rely on one single FFTW plan for all our FFTs, and perform individually each FFT along with its expansion copy. This individual method benefits from the cache memory for our small-sized expansions, and hence avoids to load data twice from main memory (for FFTs and for copies, which are both memorybound operations) as done in the batch method. This is shown in figure 5, where the individual method always outperforms (up to 30%) the batch one. We thus use the individual method in defmm.

BLAS-based upward and downward pass operators

Because of the directional aspects of defmm, the cost of the upward and downward passes is linearithmic for surface meshes (see [START_REF] Engquist | Fast directional multilevel algorithms for oscillatory kernels[END_REF]). As opposed to FMMs for non-oscillatory kernels, these steps have a significant cost in practice, which justifies their careful optimization. We hence first adapt an optimization suggested in [START_REF] Agullo | Pipelining the Fast Multipole Method over a Runtime System[END_REF] to equispaced grids and to oscillatory kernels (see section 5.2.1), and then improve it for directional FMMs (see sections 5.2.2 and 5.2.3). We validate our optimizations in section 5.2.4.

Tensorized M2M and L2L operators

We here detail a fast evaluation scheme suggested in [START_REF] Agullo | Pipelining the Fast Multipole Method over a Runtime System[END_REF] for the M2M and L2L operators on interpolation grids using tensorized Chebyshev rules in the low-frequency regime. We extend this scheme to equispaced grids which are also tensorized grids. The bijection in definition 1 induces a node indexing in the interpolation grid G allowing to exploit this tensorized structure.

Definition 1. Let L ∈ N. I denotes the bijection from [[0, L d -1]] to [[0, L -1]] d such that I -1 (I) := d k=1 I k L k-1 , ∀ I := (I 1 , ..., I d ).
Thanks to this tensorized structure, the matrix representations of the M2M and L2L operators are tensorized matrices. The L2L case being obtained by transposition, we focus on the M2M case. Let M ∈ R L d ×L d be the matrix representation of a M2M operation, there exists therefore p) . Using definition 1, the following holds

M (p) ∈ R L×L , p ∈ [[1, d]], such that M = d ⊗ p=1 M (
d ⊗ k=1 M (k) i,j = d k=1 M (k) I(i) k ,I(j) k . For any v ∈ C L d , we have (Mv) i = d ⊗ k=1 M (k) v i = L d -1 j=0 d p=1 M (p) I(i)p,I(j)p v j = J∈[[0,L-1]] d d p=1 M (p) I(i)p,Jp v I -1 (J) = J∈[[0,L-1]] d d-1 p=1 M (p) I(i)p,Jp M (d) I(i) d ,J d v I -1 (J) = J∈[[0,L-1]] d J d =0 d-1 p=1 M (p) I(i)p,Jp L-1 q=0 M (d) I(i) d ,q v I -1 (J+qe d )
where the last sum over q matches a matrix-vector product of size L × L. For i varying, this matrixvector product is performed on L d-1 different restrictions of v, involving each time the same matrix. Hence, there exists a permutation

P ∈ R L d ×L d such that (Mv) i = J∈[[0,L-1]] d J d =0 d-1 p=1 M (p)
I(i)p,Jp diag(M (d) )P v I(i) d
where diag(M (d) ) is a block-diagonal matrix with all diagonal bocks equal to M (d) . This process can be repeated d times, leading to an overall complexity of O dL d+1 since the permutations are applied in O(L d ) operations. This compares favorably with the O L 2d complexity of a naive approach.

Since the same matrix is used for multiple vectors at each of the d iterations, matrix-vector products can be stacked into matrix-matrix products to benefit from the level-3 BLAS higher efficiency [START_REF] Agullo | Pipelining the Fast Multipole Method over a Runtime System[END_REF].

We will refer to this version as the tensorized (or t) method.

The extension of the tensorized method to oscillatory kernels is obtained by noting that the directional M2M matrix M(u) with direction u can be written

M (u) = D 0 (u) d ⊗ p=1 M (p) D 1 (u) (22)
with two diagonal matrices D 0 (u) and D 1 (u) composed of complex exponential evaluations (see section 2.2).

Directional stacking

The tensorized method can be further optimized in the high-frequency regime. Starting from equation 22, if we consider for a given cell two directions u and u ′ with corresponding directional multipole expansions v(u) and v(u ′ ), denoting by ⊙ the Hadamard product, the "stacking"

M(u)v(u) M(u ′ )v(u ′ )
can be expressed as

D 0 (u) D 0 (u ′ ) ⊙ d ⊗ p=1 M (p) • D 1 (u) D 1 (u ′ ) ⊙ v(u) v(u ′ ) .
⊗ d p=1 M (p) here applies to a matrix instead of a vector, which allows to further benefit from the level-3 BLAS efficiency. Contrary to the tensorized method where L d-1 vectors can be stacked, the number of vectors that can be stacked is now L d-1 multiplied by the number of effective directional multipole expansions in each cell. This method is referred to as the tensorized+stacking method (denoted t + s).

Benefiting from real matrix products

One may observe that the matrices M (p) are real since they are composed of evaluations of Lagrange polynomials, but are applied to complex vectors (at least in a directional method). For any v ∈ C L d , we thus have

d ⊗ p=1 M (p) v = d ⊗ p=1 M (p) Re{v} + i d ⊗ p=1 M (p) Im{v}. 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 ' t ( 3 ) t + r ( 3 ) ' t ( 4 ) t + r ( 4 ) ' t ( 5 ) t + r ( 5 ) ' t ( 6 ) t + r ( 6 ) Time (s) lfM2M
(a) Sphere κD = 0

0 0.5 1 1.5 2 2.5 3 3.5 4
(b) Sphere κD = 64

0 0.2 0.4 0.6 0.8 1 1.2 ' t ( 3 ) t + r ( 3 ) ' ' t ( 4 ) t + r ( 4 ) ' ' t ( 5 ) t + r ( 5 ) ' ' t ( 6 ) t + r ( 6 ) Time (s) lfM2M (c)
(d) Cube κD = 64

Figure 6: Timings of all M2M evaluations for the tensorized (t), tensorized+stacking (t + s) and tensorized+stacking+real (t + s + r) methods, for different distributions with 10 7 particles and for different frequency regimes. κD refers to the length of the particle distribution multiplied by the wavenumber κ. 1D interpolation orders are indicated inside parentheses. lfM2M (respectively hfM2M) refers to the low-(resp. high-) frequency M2M operations. See section 6.2 for a detailed test case description.

By deinterleaving the real and imaginary parts in the expansion vectors, one can obtain stacked matrices of real elements. Hence, by discarding the imaginary part of the (real) M2M matrices, we can halve the number of required arithmetic operations. This method is referred to as the tensorized+stacking+real method (denoted t + s + r). Since there is only one multipole and one local expansion per cell in the low-frequency regime, no stacking of the directional expansions is performed and the t + s + r method is reduced to t + r in this frequency regime.

Performance results

In figure 6 are detailed the M2M timings for these three optimizations. Regarding the sphere with κD = 64 (see figure 6b), the new t + s method fastens the M2M evaluations of the t method in the high-frequency regime, while our t + s + r method further reduces the M2M evaluation times of this t + s method in the two frequency regimes. In the end, our t + s + r method outperforms the original t method of [START_REF] Agullo | Pipelining the Fast Multipole Method over a Runtime System[END_REF] by a factor up to 2.16×. This is also valid for the cube with κD = 64 (see figure 6d), when considering low interpolation orders: for L = 3 the t + s + r method outperforms the t one by a factor of ≈ 4×. For larger interpolation orders, the t + s + r method is however less efficient than the t + s one, and even leads to a performance loss for L = 6. This is due to our small matrix sizes (one matrix being of size L × L) for which the t + s method already provides a high enough number of stacked vectors with the cube distribution: there is thus no benefit with a higher number of stacked vectors (as provided by the t + s + r method). Moreover, the deinterleaving cost is not here offset by the lower number of operations induced by the real matrix products.

In the low-frequency regime (see figures 6a,6b,6c), the deinterleaving still allows faster M2M evaluations than the original t method with performance gains up to 4x. Since we target surface particle distributions for BIE problems, we choose to rely on the t + s + r method in defmm for M2M evaluations, as well as for L2L ones (for which the same performance gains have been obtained).

Vectorization for direct computation

The direct computation involved in the P2P operator consists in two nested loops over the (possibly distinct) target and source particles. As usually performed in FMM codes, we aim at vectorizing this computation to benefit from the AVX2 or AVX-512 units. We target here the vectorization of the outer loop (on target particles), which leads to fewer reduction operations than with an inner loop vectorization. Based on the data locality among particles described in section 4.3.1, our P2P code is close to the one presented in [START_REF] Abduljabbar | Extreme scale fmm-accelerated boundary integral equation solver for wave scattering[END_REF]. As in [START_REF] Abduljabbar | Extreme scale fmm-accelerated boundary integral equation solver for wave scattering[END_REF], we load separately the real and imaginary parts of the complex numbers and we perform only real (not complex) mathematical operations: this is required for our compiler (Intel C++ compiler) to fully vectorize the code. There are nevertheless some differences between our code and the one presented in [START_REF] Abduljabbar | Extreme scale fmm-accelerated boundary integral equation solver for wave scattering[END_REF]:

• Since the inputs of our code are the particle point clouds, we have to be able to numerically handle the practical case of interactions between equally located particles on which the kernel function cannot be evaluated. These singularities (i.e. the interactions between a particle and itself) are resolved directly: if the distance between the particles is too small, we set one operand to 0.0 and we continue the interaction computation, which hence results in 0.0. This minimal test may be processed with masks in the vectorized code, hence leading to a minimal performance penalty.

• In order to vectorize our outer loop, an OpenMP compiler directive is required for the Intel C++ compiler (which targets the inner loop otherwise). It can be noticed that the GNU C++ compiler currently fails to vectorize our code (with or without the OpenMP directive).

• Since the particle positions are not updated from one FMM application to the next in BIE problems, we can rely on a structure-of-array data layout for the potentials and the charges, which favors efficient vector memory accesses.

In the end, using 64-bit floating-point numbers on 512-bit AVX-512 units (hence a maximum speedup of 8×), our vectorized P2P operator offers performance gains up to 7.6× over the original scalar implementation (see [START_REF] Chollet | Symmetries and Fast Multipole Methods for Oscillatory Kernels[END_REF] for details).

Numerical results

We now provide numerical results illustrating the convergence and the performance of defmm. All codes run sequentially on a Intel Xeon Gold 6152 CPU with AVX-512 units and 384 GB of RAM.

Relative error

We first check the overall accuracy of defmm on a (volumic) uniform cube test case composed of 125000 particles. We compute the error of a defmm approximation p of p (see equation ( 1)) for a norm || • || as ||p-p|| ||p|| . The charges are randomly chosen. As shown in figure 7, the convergence follows the estimate derived from the interpolation process on equispaced grids in section 3.2, i.e. a geometric convergence in the 1D interpolation order. However, the same interpolation order leads to lower accuracies in the high-frequency regime than in the low-frequency one: this is due to the partitioning of the direction set, which introduces additional approximations in the high-frequency regime compared to the non-directional interpolation in the low-frequency regime. We believe that this is a consequence of the small initial set of directions chosen for performance reasons in the defmm implementation (see section 4.2.2) and that increasing the number of possible directions at each tree level should allow to recover better accuracies. Similar results were obtained on a surface sphere distribution.

Performance comparison with dfmm

For this performance comparison, we consider the following test cases, all with 10 7 particles (±1%).

• The uniform cube test case (volumic distribution). This classical FMM test case is not typical for BIE problems, but is still interesting due to its high computation cost for directional FMMs (the number of effective expansions being maximum).

• The classical sphere test case where the particles are quasi-uniformly scattered on the surface of the unit sphere. This is also considered as a difficult case, with high memory requirements and high computation costs, for the directional methods on boundary particle distributions [START_REF] Engquist | Fast directional multilevel algorithms for oscillatory kernels[END_REF].

• The refined cube, where the particles are highly non-uniformly sampled along the surface of the unit cube with a higher concentration around the edges and corners. This mimics distributions obtained with mesh refinement.

• The ellipse, which is an elongated distribution favorable to directional FMMs due to the reduced number of directions in the high-frequency regime. Our ellipse distribution presents a higher concentration of particles on the two poles.

defmm is compiled using the Intel C++ compiler. All the BLAS calls use the Intel Math Kernel Library (MKL) and we rely on the FFTW3 library for the FFTs.

Regarding dfmm (only g++ supported), we consider its two best variants [START_REF] Messner | Optimized M2L kernels for the Chebyshev interpolation based fast multipole method[END_REF], both based on low-rank compressions and on interpolations using tensorized Chebyshev grids. The first variant (IAblk -Individual Approximation with symmetries and BLocKing) relies on vector stacking and symmetries to improve performance with level-3 BLAS routines. The second one (SArcmp -Single Approximation with ReCoMPression) uses global approximations of the M2L matrices and recompressions to reduce the numerical ranks of these matrices. According to [START_REF] Messner | Optimized M2L kernels for the Chebyshev interpolation based fast multipole method[END_REF], this can lead to faster M2L evaluations than with IAblk, but at the cost of longer precomputation times.

Following the methodology used in [START_REF] Messner | Fast directional multilevel summation for oscillatory kernels based on Chebyshev interpolation[END_REF], we fix the threshold in the low-rank approximations of dfmm to 10 -L . We checked that the accuracy of the Chebyshev interpolation in dfmm is then similar to the accuracy of the equispaced interpolation in defmm. For each test, the M axDepth (dfmm) and the N crit (defmm) values (see section 4.3.1) are tuned to minimize the FMM application time. Due to memory constraints, we were not able to run dfmm (with IAblk or SArcmp) on the refined cube for L = 6. The comparison results are given in figure 8.

Overall performance Considering the overall times (including precomputation and application times), defmm always outperforms dfmm, being 1.3× to 6.5× faster than its best variant (IAblk here). As detailed below, this is due to the different design of the two codes. Since the P2P operator is not vectorized in dfmm, we also present in figure 8 a defmm variant, denoted by defmm-S, with a scalar (i.e. non-vectorized) P2P implementation. This enables us to more precisely study the impact of the other differences between defmm and dfmm. defmm-S offers similar or better performance than the best dfmm variant on the uniform cube and on the sphere, and outperforms dfmm by factors 1.8× to 5.9× on more non-uniform distributions (i.e. the refined cube and the ellipse). This shows that regarding directional FMMs our FFT-based approach is competitive with or faster than the low-rank approximations used in dfmm.

Sensitivity to the particle distribution Since the distributions have the same number of particles (±1%), comparing two distributions (all other test parameters being identical) illustrates the sensitivity of the method to the particle distribution. Except for the uniform cube in the highfrequency regime (whose higher cost is due to its maximum number of effective directional expansions), the defmm performance is few sensitive to the distribution, contrary to the dfmm one. This is due to the differences in the tree construction and traversal (see section 4.3): the combination of the N crit strategy and of our specific DTT allows to better adapt in defmm to the particle distribution than the M axDepth strategy used in dfmm. Besides, the performance gap between IAblk and defmm-S is minimized for the (uniform) cube and the sphere (quasi-uniform on its surface) distributions. These two distributions benefit indeed most from the BLAS routines and from the vector stacking in IAblk.

Sensitivity to the wavenumber Since both codes are based on the same directional approach, their performance is similarly impacted by an increase in the wavenumber. Again, the uniform cube in the high-frequency regime is a special case (maximizing the number of effective expansions), where the dfmm features (level-3 BLAS routines, cache management) lower here the performance gap with defmm when the wavenumber increases: dfmm is however still outperformed by defmm in this case.

Sensitivity to the interpolation order When moving from L = 4 to L = 6, the average performance ratio of defmm-S over SArcmp moves from 1.58× to 2.40× (for the uniform cube, the sphere and the ellipse test cases, and for any frequency). This increasing performance gain with respect to the interpolation order is due to our FFT-based FMM. With respect to IAblk, the performance gain of defmm-S moves from 1.42× to 1.56×: here the FFT gain is counterbalanced by the BLAS performance gain (which is greater for larger matrices).

Precomputation costs Thanks to our extension of the symmetry usage to the Fourier domain (see section 4.4), we manage to have on the uniform cube and on the sphere precomputation costs in defmm as low as the dfmm ones (considering the IAblk variant which requires shorter precomputations than SArcmp). In addition, in highly non-uniform distributions such as the refined cube and the ellipse, the precomputation costs of defmm are drastically lower than the dfmm ones, partly thanks to our N crit-based 2 d -trees and to our blank passes (see section 4.3) which minimize the required precomputations. Figure 9: Timings of precomputations and applications of defmm on the sphere test case with respect to number of particles at fixed wavenumber κ (left), and at constant number of particles per wavelength (right).

Specific study of the sphere test case

In this section, we further detail performance results on the sphere test case, regarding various parameters.

Scalings

In Fig. 9, we provide timings results obtained with defmm on the sphere with respect to the number of particles in two situations:

1. using a fixed wavenumber, hence only increasing the number of particles without impacting the frequency;

2. using a constant (on average) number of particles per wavelength (i.e. increasing the wavenumber with the number of particles).

In the second situation (see Fig. 9 (right)), both the precomputation and the application timings scale as O(N log N ), which is in accordance with the theory. In the first situation (see Fig. 9 (left)), the results indicate that adding particles without increasing the wavenumber results in adding lowfrequency cells in our trees. Since our method scales as O(N ) in the low-frequency regime (as confirmed by the K = 0 case), we thus asymptotically recover this scaling even for non-zero but fixed wavenumbers (such as K = 64). The specific shape of the scaling curves for K = 64 is merely due to the lower cost of processing low-frequency cells compared to the high-frequency ones.

Comparison with another directional FMM

Finally, we provide an additional comparison with the kernel-independent directional FMM from [START_REF] Engquist | Fast directional multilevel algorithms for oscillatory kernels[END_REF]. This FMM, based on KIFMM [START_REF] Ying | A kernel-independent adaptive fast multipole algorithm in two and three dimensions[END_REF] and referred to as ki-dfmm here, exploits QR-based randomized low-rank approximations rather than polynomial interpolations to process the far-field contributions.

Using defmm, we have reproduced in Table 1 the same performance test as in [START_REF] Engquist | Fast directional multilevel algorithms for oscillatory kernels[END_REF], considering their most difficult test case (i.e. the sphere) on the exact same particle distributions. Direct timings were estimated using the same methodology as in [START_REF] Engquist | Fast directional multilevel algorithms for oscillatory kernels[END_REF]. All tests have been executed in sequential, but since there is at least a 10-year gap between the CPU from [START_REF] Engquist | Fast directional multilevel algorithms for oscillatory kernels[END_REF] (published in 2007) and our Intel Xeon Gold 6152 processor (released in 2017), we focus in the comparison on the performance ratios between direct and FMM computation times, rather than on the times in seconds. A denotes the ratio of the direct computation time over the accelerated computation time. Perfomance numbers given in italic are extracted from [START_REF] Engquist | Fast directional multilevel algorithms for oscillatory kernels[END_REF].

One can first notice that, due to memory constraints, defmm fails to run tests requiring a high precision. Indeed, the memory consumption of our FFT-based approach is far more important than the ki-dfmm one, especially on the sphere test case which implies a high number of directional expansions per cell. However, when the defmm execution succeeds, defmm offers each time a higher performance ratio (over the direct computation) than ki-dfmm, which confirms the interest of our approach.

Conclusion

In this paper we presented a new approach for the directional interpolation-based FMMs, using FFT techniques thanks to equispaced grids. We provided a consistency proof of the approximation process and showed how to extend the symmetries of the interpolation-based FMM to the Fourier domain. We detailed the algorithmic design of our defmm library, as well as its high-performance optimizations on one CPU core. Finally, comparisons two state-of-the-art libraries exhibited the superior performance of our library in all test cases (without memory constraints), the defmm performance being also few sensitive to the surface particle distribution. This shows the interest of our FFT-based approach for lowering the FMM cost in the high frequency regime, and for efficiently processing the non-uniform distributions of the BIE problems.

Future works will be dedicated to the defmm parallelization. In this purpose, we will be able to lean on the dual tree traversal, which is known to be highly efficient regarding shared-memory parallelism (see e.g. [START_REF] Abduljabbar | Extreme scale fmm-accelerated boundary integral equation solver for wave scattering[END_REF][START_REF] Yokota | An FMM based on dual tree traversal for many-core architectures[END_REF]), and on previous work (such as [START_REF] Benson | A parallel directional fast multipole method[END_REF]) regarding distributed-memory parallelism. We also plan to integrate defmm in an iterative solver in order to solve complex realistic boundary integral equations.

1 Introduction

 1 Considering two point clouds X, Y ⊂ R d with cardinal N ∈ N * , d ∈ N * and q ∈ C[Y ] (named charges), where C[Y ] denotes the set of application from Y to C, we are interested in the fast computation of p ∈ C[X] (referred to as potentials) such that p(x) := y∈Y G(x, y)q(y), ∀ x ∈ X,

Figure 1 :

 1 Figure 1: 2 d -tree (d = 2) representation of B with particles colored in blue.

Figure 3 :

 3 Figure 3: Operators and expansions in defmm. hf is added in front of the directional operators (used in the high-frequency regime). P: particles; M i : multipole expansions; L i : local expansions; Mi , Li : multipole and local expansions in the Fourier domain.

Algorithm 3

 3 Blank Downward Pass (BDP) for cell c 1: if c is in the high-frequency regime then 2: for each direction u with which c is marked do 3: for c ′ ∈ Sons(c) do 4: mark c ′ with F ather(u)

1 . 2 ) 4 .Figure 5 :

 1245 Figure5: Timings of all FFT applications with 1D interpolation order 4, using the batch and individual methods. Tests performed with 10 7 particles (see section 6 for details).

Figure 7 :

 7 Figure 7: Relative error of the FMM over a 3D uniform cube with 1.25 × 10 5 particles. Maximum norm: squares; l 1 norm: circles; l 2 norm: triangles.

  L = 6, κD = 64

Figure 8 :

 8 Figure 8: Timings for one FMM execution of the dfmm and defmm variants. The two tested variants of dfmm are IAblk and SArcmp. defmm-S corresponds to defmm without P2P vectorization.

1 :

 1 {A: strict (resp. directional) MAC in low-(resp. high-) frequency regime} 2: if A(t, s) then

	3:	apply the M2L operator between t and s
	4:	return
	5: else
	6:	if t is a leaf or s is a leaf then
	7:	apply the P2P operator between t and s
	8:	return
	9:	end if
	10:	for t ′ ∈ Sons(t) do
	11:	for s ′ ∈ Sons(s) do
	12:	DTT(t ′ , s ′ )
	13:	end for
	14:

Table 1 :

 1 Timings (in seconds) for estimated direct computation (T d ) and accelerated (FMM) application (T a ).

	K 2π , N, ϵ	T def mm a	T d	A def mm T ki-df mm a	T d	A ki-df mm
	(16,3.22e+5,1e-4) 1.58e+1 7.89e+3 5.00e+2 9.50e+1	1.21e+4	1.28e+2
	(32,1.29e+6,1e-4) 7.73e+1 1.28e+5 1.66e+3 4.28e+2	1.95e+5	4.55e+2
	(64,5.15e+6,1e-4) 4.01e+2 2.01e+6 5.03e+3 1.97e+3	3.04e+6	1.54e+3
	(16,3.22e+5,1e-6) 4.36e+1 7.85e+3 1.80e+2 2.42e+2	1.18e+4	4.86e+1
	(32,1.29e+6,1e-6) 2.38e+2 1.28e+5 5.37e+2 1.21e+3	1.87e+5	1.54e+2
	(64,5.15e+6,1e-6)	-	-	-	5.95e+3	3.13e+6	5.27e+2
	(32,3.22e+5,1e-8) 9.31e+1 7.89e+3 8.47e+2 5.11e+2	1.22e+4	2.39e+1
	(32,1.29e+6,1e-8)	-	-	-	2.62e+3	1.96e+5	7.51e+1
	(64,5.15e+6,1e-8)	-	-	-	1.25e+4	3.15e+6	2.52e+2

Binary tree (d=1), quadtree (d=2) or octree (d=3).

At: https://github.com/IChollet/defmm

i.e. the radius of the smallest ball containing the cell.

Two cells at the same 2 d -tree level comply with the strict MAC if the distance between them is greater or equal than their side length.

State-of-the-art parallel C++ (non-directional) FMM library: https://github.com/exafmm/exafmm

See: https://www.openmp.org/

See: http://www.fftw.org/fftw3_doc/Advanced-Complex-DFTs.html