Olivier Coulaud
email: olivier.coulaud@inria.fr

Pierre Fortin
email: fortin@labri.fr

Jean Roman
email: roman@labri.fr

High performance BLAS formulation of the multipole-to-local operator in the Fast Multipole Method (submitted)

The multipole-to-local (M2L) operator is the most time-consuming part of the far field computation in the Fast Multipole Method for Laplace equation. Its natural expression, though commonly used, does not respect a sharp error bound: we here first prove the correctness of a second expression. We then propose a matrix formulation implemented with BLAS (Basic Linear Algebra Subprograms) routines in order to speed up its computation for these two expressions. We also introduce special data storages in memory to gain greater computational efficiency. This BLAS scheme is finally compared, for uniform distributions, to other M2L improvements such as block FFT, rotations and plane wave expansions. When considering runtime, extra memory storage, numerical stability and common precisions for Laplace equation, the BLAS version appears as the best one.

Introduction

The N-body problem in numerical simulations describes the computation of all pairwise interactions among N bodies. The Fast Multipole Method (FMM) developed by Greengard & Rokhlin [START_REF] Greengard | A fast algorithm for particle simulations[END_REF] [START_REF] Greengard | The Rapid Evaluation of Potential Fields in Particle Systems[END_REF] solves this N-body problem for any given precision with O(N ) runtime complexity against O(N 2 ) for the direct computation. The potential field is indeed decomposed in a near field part, directly computed, and a far field part approximated thanks to multipole and local expansions. First introduced for gravitational potentials in astrodynamics or electrostatic (coulombic) potentials in molecular dynamics [START_REF] Greengard | A fast algorithm for particle simulations[END_REF] [START_REF] Elliott | Multipole algorithms for molecular dynamics simulation on high performance computers[END_REF], it has then been extended with different mathematical bases to electromagnetism [START_REF] Darve | Méthodes multipôles rapides : résolution des équations de Maxwell par formulations intégrales[END_REF], VLSI capacitance [START_REF] Nabors | Fastcap: A multipole accelerated 3-D capacitance extraction program[END_REF], radiosity [START_REF] Singh | Load balancing and data locality in adaptive hierarchical N-body methods: Barnes-Hut, fast multipole, and radiosity[END_REF], object modeling [START_REF] Carr | Reconstruction and representation of 3D objects with radial basis functions[END_REF] and many more. In this paper we will focus on the Laplace equation and we will therefore only consider gravitational and electrostatic potentials.

In the FMM computation for Laplace equation the multipole-to-local operator (M2L), which converts a multipole expansion into a local expansion, represents most of the runtime of the far field computation. The operation count of this operator is O(P 4 ) in the 3D FMM, where P is the maximum degree in the expansions. In order to fasten its computation, some schemes have been introduced that reduce the O(P 4 ) operation count to O(P 3 ) or O(P 2 ): (block) Fast Fourier Transform (FFT) [START_REF] Elliott | Fast Fourier Transform accelerated fast multipole algorithm[END_REF], rotations [START_REF] White | Rotating around the quartic angular momentum barrier in fast multipole method calculations[END_REF] and more recently plane wave expansions [START_REF] Greengard | A new version of the fast multipole method for the Laplace equation in three dimensions[END_REF] [START_REF] Cheng | A fast adaptive multipole algorithm in three dimensions[END_REF]. When considering operation count and numerical stability, the plane wave improvement appears to be the best one, but to our knowledge no other work has been published about this plane wave scheme 1 with Laplace equation1 since the article [START_REF] Cheng | A fast adaptive multipole algorithm in three dimensions[END_REF]. The FFT and rotation schemes are thus still used in this context: see [START_REF] Kurzaka | Communications overlapping in fast multipole particle dynamics methods[END_REF] and [START_REF] Gumerov | Recursions for the computation of multipole translation and rotation coefficients for the 3-D Helmholtz equation[END_REF].

However, on modern architectures the computation speed of the processor is much higher than the speed of memory access. This difference leads to important waiting times in order to access data in memory: the processor is then unused, which dramatically affects the CPU times of these simulations. This is why we propose here a different approach: we will use highly efficient implementation techniques such as BLAS (Basic Linear Algebra Subprograms) routines [START_REF] Dongarra | A set of level 3 Basic Linear Algebra Subprograms[END_REF] to improve the runtime of the FMM. Thanks to optimal use of the different layers of the hierarchical memory of the computer, so that the pipelines of the floating point units are filled at best, they indeed offer substantial runtime speedup on superscalar architectures. This speedup only affects the constant in the O(P 4 ) notation and we keep the O(P 4 ) operation count. But in molecular dynamic simulations for example, the required precisions for electrostatic potentials range usually between 10 -5 et 10 -7 , so P ranges from 3 to 15 (see [START_REF] Elliott | Multipole algorithms for molecular dynamics simulation on high performance computers[END_REF] or [START_REF] Rankin | Efficient parallel implementations of multipole based N-body algorithms[END_REF]); and for gravitational potentials that arise in astrophysics the required precisions and the P values are even lower [START_REF] Capuzzo-Dolcetta | A comparison between the fast multipole algorithm and the tree-code to evaluate gravitational forces in 3-D[END_REF]. All these values of P are quite low in terms of operation count, and the speedup obtained with the BLAS routines can thus exceed the one obtained with a scheme that offers a lower operation count.

We also distinguish two different expressions of the M2L operator. Indeed, the error bound of the 3D FMM has been historically ( [START_REF] Greengard | A fast algorithm for particle simulations[END_REF] [START_REF] Greengard | The Rapid Evaluation of Potential Fields in Particle Systems[END_REF]) presented for the evaluation of potential with either finite multipole expansions or finite local expansions. But, as mentioned by several authors ( [START_REF] Soelvason | Error estimates for the fast multipole method[END_REF], [START_REF] White | Derivation and efficient implementation of the fast multipole method[END_REF]), we have also to consider, when implementing the FMM, that the M2L operator acts on finite multipole expansions, which means that both multipole and local expansions are finite. When denoting P the maximum degree of the expansions, and n (respectively j) the degrees of the multipole (respectively local) expansion terms, two different kinds of M2L expressions can then be used. In the first one, both n and j fully range between 0 and P , whereas in the second M2L expression we use only terms with n + j ≤ P . While the first one is natural and commonly used ( [START_REF] Petersen | Error estimates for the fast multipole method. II. The three-dimensional case[END_REF], [START_REF] Soelvason | Error estimates for the fast multipole method[END_REF], [START_REF] White | Derivation and efficient implementation of the fast multipole method[END_REF]), no sharp error bound has yet been found, to our knowledge, for the corresponding summations. We will prove that the second one, though generally less efficient, respects such sharp error bound. For these two expressions, we will present the principles and the implementation features for two improvements of the M2L operator, namely FFT and rotations, as well as for our BLAS approach. We will then propose a detailed comparaison of the memory requirements, numerical precisions and runtimes, for uniform distributions sequentially computed, between these three schemes depending on the M2L expression used.

Thanks to the FMMPART3D code (version 1.0) distributed by the MadMax Optics company2 , we will also compare our BLAS version with the enhancement based on plane wave expansions. To our knowledge, this is the first comparison between this scheme with plane waves and another improvement of the M2L operator when considering the FMM for Laplace equation, since only comparisons with direct computation are given in [START_REF] Greengard | A new version of the fast multipole method for the Laplace equation in three dimensions[END_REF] and [START_REF] Cheng | A fast adaptive multipole algorithm in three dimensions[END_REF].

The rest of this article is organised as follows: in the next section we will introduce the formulae used in our implementation of the FMM. We will also discuss one error bound issue and present the two different versions of the M2L operator. We will expose and discuss the FFT enhancement in section 2.3, the use of rotations in section 2.4, the plane wave scheme in section 2.5 and the introduction of BLAS routines in section 3. The FFT, rotation and BLAS improvements have been implemented in a code named FMB for Fast Multipole with BLAS (or Fast Multipole in Bordeaux) thus allowing precise comparisons among them as presented in section 4. For almost all these sections, more details and explanations can be found in the research report [START_REF] Fortin | Multipole-to-local operator in the Fast Multipole Method: comparison of FFT, rotations and BLAS improvements[END_REF]. The current article has to be seen as a summary of this research report, except for the final comparison in section 4.6 between the scheme with plane waves and our BLAS approach.

Multipole to Local operator: presentation and current improvements

We present here the M2L operator, its mathematical formulae, the two possible expressions, their consequence on the error bound and two current improvements of its operation count that we have implemented: the FFT enhancement and the rotations. The plane wave scheme is also briefly presented.

FMM presentation

We focus in this article on uniform distributions sequentially computed and refer the reader to the articles of Greengard & Rokhlin for a presentation of the FMM and especially the notions of quadtree and octree, near and far fields, multipole and local expansions, upward and downward passes, nearest neighbors, interaction list, and well-separateness. We denote ws the wellseparateness criterion (see [START_REF] White | Derivation and efficient implementation of the fast multipole method[END_REF]) and use ws = 1 here, which means that two cells of the octree are well-separated provided they do not share a boundary point. The interaction list of a cell c is then defined as the set of all children of the nearest neighbors of the parent cell of c which are not themselves nearest neigbors of c; there is 189 members in each interaction list. We denote by M2M the operator that translates a multipole expansion. M2L denotes the conversion of a multipole expansion into a local expansion, and L2L the translation of a local expansion. We consider the root of the octree (or computational box) to be at level 0. The height of the octree is defined as its maximum level. Moreover we use Morton ordering for the indexing of each cell: this allows a fast access to all cells in the octree thanks to bit operations (see [START_REF] Elliott | Multipole algorithms for molecular dynamics simulation on high performance computers[END_REF]).

We introduce now briefly the formulae used in our implementation of the FMM focusing on the M2L operator: more details and the complete formulae set can be found in [START_REF] Fortin | Multipole-to-local operator in the Fast Multipole Method: comparison of FFT, rotations and BLAS improvements[END_REF]. We have choosen the formulae of Epton & Dembart [START_REF] Epton | Multipole translation theory for the three-dimensional Laplace and Helmholtz equations[END_REF] for the clarity of the underlying theorems and their proofs, and also because they allow elementary formulation of the M2L operator, which represents the most time-consuming part of the far field computation, while obeying the symmetry property among the multipole and local expansion terms of opposite order (see lemma 1 at the end of the section).

In 3D space, θ denotes the co-latitudinal coordinate and φ the longitudinal coordinate. With m = (-1) m if m ≥ 0, and m = 1 otherwise, and P m l denoting the associated Legendre function, our unnormalized spherical harmonics Y m l of degree l and order m, with l ≥ 0 et -l ≤ m ≤ l, are defined by

Y m l (θ, φ) = m (l -m)! (l + m)! P m l (cos θ)e i.m.φ . (1)
We emphasize that our spherical harmonics are considered as identically null for l < 0 or |m| > l.

Like Epton & Dembart [START_REF] Epton | Multipole translation theory for the three-dimensional Laplace and Helmholtz equations[END_REF], we also respectively define the Outer and Inner functions by

O m l (r, θ, φ) = (-1) l i |m| A m l Y m l (θ, φ) 1 r l+1 ∀ (l, m) ∈ N × Z with 0 ≤ |m| ≤ l, (2)
and

I m l (r, θ, φ) = i -|m| A m l Y m l (θ, φ)r l ∀ (l, m) ∈ N × Z with 0 ≤ |m| ≤ l, (3)
where

A m l = (-1) l √ (l-m)!(l+m)! .
Letting X and X be two position vectors in 3D space, we now give the main theorems that the FMM requires: their proof can be found in [START_REF] Epton | Multipole translation theory for the three-dimensional Laplace and Helmholtz equations[END_REF].

Theorem 1 (Classical translation theorem) Under the assumption X > X , we have

1 X -X = +∞ n=0 n l=-n (-1) n I -l n (X )O l n (X).
The next theorem is used to establish M2M (Outer-to-Outer) and M2L (Outer-to-Inner) operators.

Theorem 2 (Outer-to-Outer, Outer-to-Inner Laplace translation theorem) Let X > X , then

O l n (X -X ) = +∞ j=0 j k=-j (-1) j I -k j (X )O l+k n+j (X) ∀ (n, l) ∈ N × Z with 0 ≤ |l| ≤ n.
The last theorem corresponds to the Third Addition Theorem in [START_REF] Greengard | The Rapid Evaluation of Potential Fields in Particle Systems[END_REF], and it is used for L2L (Inner-to-Inner) operator.

Theorem 3 (Inner-to-Inner Laplace translation theorem)

I l n (X -X ) = n j=0 j k=-j (-1) j I k j (X )I l-k n-j (X) ∀ (n, l) ∈ N × Z with 0 ≤ |l| ≤ n.
Thus, according to the classical translation theorem and the Outer-to-Inner Laplace translation theorem, we have the following expression for the M2L operator that converts multipole expansion terms into local expansion terms, as defined in [START_REF] Fortin | Multipole-to-local operator in the Fast Multipole Method: comparison of FFT, rotations and BLAS improvements[END_REF].

Proposition 1 (M2L operator) Let M l n , with n ≥ 0, |l| ≤ n, the multipole expansion terms centered in z 1 , the local expansion terms L k j centered in z 2 , write

L k j = +∞ n=0 n l=-n M l n O -k-l j+n (ρ, α, β) ∀ (j, k) ∈ N × Z with 0 ≤ |k| ≤ j
where (ρ, α, β) are the spherical coordinates of the M2L vector z 2 -z 1 .

The O -k-l j+n (ρ, α, β) are named the M2L transfer functions . The implementation of this formulae is named hereafter the classic M2L. Finally, we emphasize the following property which is also valid for the Outer and Inner functions, as well as for the multipole and local expansion terms: this enables the computation of only terms with positive orders in multipole and local expansions.

Error bound analysis

As exposed in the introduction, the M2L operator concretely uses finite multipole expansion to compute (finite) local expansions. When denoting P as the maximum degree of the expansions, whose terms have degrees therefore ranging from 0 to P , two different kinds of M2L expressions can then be used.

The first one, named M2L kernel3 with double height, corresponds to the outer sum on n, in the expression of the M2L operator given in proposition 1, stopping at n = P , since in this case the O k j terms range up to 2P . This one is natural and commonly used, but even if some interesting works have been done to estimate the behavior of the error induced by such M2L expression (see [START_REF] Petersen | Error estimates for the fast multipole method. II. The three-dimensional case[END_REF] [START_REF] Soelvason | Error estimates for the fast multipole method[END_REF]), no sharp error bound has been found yet.

In the second M2L expression the outer sum on n stops at n = Pj, so that the maximum degree of the O k j used is P . This is named M2L kernel with single height. This has been used for example in the DPMTA (Distributed Parallel Multipole Tree Algorithm) code developped at Duke University but the proof in the error bound given in [START_REF] Elliott | Multipole algorithms for molecular dynamics simulation on high performance computers[END_REF] is wrong 4 . We however recommend the reading of this appendix since it presents a worst case error analysis similar to the one used hereafter. In particular it shows that no additional error is introduced by the M2M operation. As explained in [START_REF] Greengard | The Rapid Evaluation of Potential Fields in Particle Systems[END_REF] the L2L operation does not introduce any error at all. We can thus focus on the M2L operation.

In the following we present a sharp error bound (similar to the one presented in [START_REF] Elliott | Multipole algorithms for molecular dynamics simulation on high performance computers[END_REF]) respected by an M2L operator with single height kernel. It has to be noted that this error bound for single height M2L kernel was briefly presented by White & Head-Gordon [START_REF] White | Derivation and efficient implementation of the fast multipole method[END_REF], but they have used double height kernel in their implementation.

In order to introduce this error bound, we study the potential in point Z due to one single unit charge located in Q as pictured in figure 1. We denote by B 1 (respectively B 2 ) the ball centered

X X X X r2 X2 X1 Q r1 Z x y z α Figure 1: Our problem configuration. in X 1 with radius r 1 (respectively X 2 and r 2 ). Q is enclosed in B 1 , Z in B 2 . Moreover, denoting R = X 1 -X 2 , we assume that R > r 1 + r 2 .
We first proceed as in [START_REF] Greengard | The Rapid Evaluation of Potential Fields in Particle Systems[END_REF]. Let α the angle between "vectors" Q and Z, r > (respectively r < ) the maximum (respectively minimum) between the norm of Q and the one of Z, and P n the Legendre polynomials, we have

Φ(Z) = 1 Z -Q = +∞ n=0 r n < r n+1 > P n (cos α). (4)
Defining Φ P (Z) the potential with finite summation up to P , we have, since for all x ∈ [-1, 1],

|P n (x)| ≤ 1, Φ(Z) -Φ P (Z) ≤ 1 r > -r < r < r > P +1
.

(

) 5
This leads to the following error bound

Proposition 2 ∀ Q ∈ B 1 , and ∀ Z ∈ B 2 , we have, under the condition R > r 1 + r 2 Φ(Z) -Φ P (Z) ≤ 1 R -(r 1 + r 2 ) r 1 + r 2 R P +1 . Proof. Since Z -Q = ((Z -X 2 ) + (X 1 -Q)) + (X 2 -X 1 ) and since (Z -X 2 ) + (X 1 -Q) ≤ r 1 + r 2 ,
this results directly from the inequality [START_REF] Nabors | Fastcap: A multipole accelerated 3-D capacitance extraction program[END_REF].

We can now decompose the potential in the same way the FMM does while maintaining this error bound, and we obtain the following theorem (whose proof is given in appendix A).

Theorem 4 The M2L operator with single height M2L kernel strictly respects the error bound in proposition 2.

Since the error in proposition 2 is higher than the one due to the multipole expansion terms (see [START_REF] Greengard | A new version of the fast multipole method for the Laplace equation in three dimensions[END_REF]), and since no additional error is introduced by the M2M and L2L operators, the error bound of the FMM with a single height kernel M2L operator is therefore the same as in proposition 2.

It has to be noted that the condition of well-separateness, with either ws = 1 or ws = 2, is in fact more strict than R > r 1 + r 2 in order to ensure a fast enough convergence. With ws = 1 for example, we have, for a cell of side l:

r 1 = r 2 = √ 3 2 and R ≥ 2. Thus Φ(Z) -Φ P (Z) ≤ 1 R -(r 1 + r 2 ) √ 3 2 P +1 = 1 R -(r 1 + r 2 ) (0.866) P +1 .
As a result, this error bound can safely be used in the case of single height kernel in order to determine the value of P according to the required precision. Moreover since the double height kernel, compared to the single height kernel, only adds terms in the M2L formula, it respects at least this error bound but even adds precision with additional cost at runtime. The problem is that we cannot predict this gain in the accuracy of the FMM, and therefore we cannot compare single and double height kernels according to the tradeoff between theoretical accuracy and runtime. Only comparisons with practical accuracies will be possible as described in section 4.

We will now present, for both single and double kernel heights, three current improvements of the M2L operator that reduce its theoretical operation count.

Fast Fourier Transform

The use of Fast Fourier Transform (FFT) in order to speed up the M2L computation has been introduced by William D. Elliott ([3] or [START_REF] Elliott | Fast Fourier Transform accelerated fast multipole algorithm[END_REF]). This work was performed only for single height M2L kernel and the block version used to prevent numerical instability as written in DPMTA has restricted the possible values of P to multiples of the block size. In [START_REF] Fortin | Multipole-to-local operator in the Fast Multipole Method: comparison of FFT, rotations and BLAS improvements[END_REF] we have generalized this FFT improvement to both single and double height M2L kernels and our block version can handle any value for P .

The principle behind the use of FFT in FMM is to view the M2L operator as one 2D convolution (or correlation) between two sequences. This convolution is faster computed thanks to one forward 2D FFT, one point-wise product and one backward 2D FFT, provided that the two sequences are periodic. As presented in [START_REF] Epton | Multipole translation theory for the three-dimensional Laplace and Helmholtz equations[END_REF] for M2M, we have to use additional null terms in order to build periodic sequences out of the multipole and local expansions and the M2L transfer function: this process is named zero-padding. However, due to large varying magnitude in the norm of the expansion terms, this FFT computation results in numerical instabilities when P grows. This issue has partly been resolved in [START_REF] Elliott | Fast Fourier Transform accelerated fast multipole algorithm[END_REF] with the use of a block FFT that implies a large grain convolution among the blocks for the point-wise product. More details can be found in [START_REF] Fortin | Multipole-to-local operator in the Fast Multipole Method: comparison of FFT, rotations and BLAS improvements[END_REF] and in our implementation we have used FFTW [START_REF] Frigo | The design and implementation of FFTW3[END_REF] as an efficient FFT implementation.

As for operation counts the most time-consuming part of the FFT scheme is the point-wise product. While each point-wise product runs in O(P 2 ) in the non-block version, the large-grain convolution of the block version leads to an operation count in O(P 3 ). This part of the computation does not match any of the standard BLAS calls (see section 3) and offers no data reuse: no special speedup can be expected from its implementation.

The block version does not however ensure complete numerical stability as shown in figure 2. For growing values of P this may be explained by the influence of the orders in the magnitude of the norm of the terms used: the block version is indeed only applied in the degree dimension. But numerical instabilities are also subject to the octree height used: an increasing height indeed decreases the distance between two cell centers in the M2L transfer function terms, as well as the distance between a cell center and a particle location in the multipole expansion terms. Numerical instabilities hence arise when the order of magnitude of these distances is becomming too small. Of course this behaviour depends on the size of the computational box that encloses all particles, but with a bigger computational box the numerical instabilities would nevertheless appear for greater heights. These instabilities can however be reduced with lower block size: while, for single height M2L kernel, the FFT with block size 4 becomes unstable for P higher than 14 (see figure 2), a FFT with block size 3 is stable until 27, but is of course slower to compute because of the higher operation count due to the higher number of blocks. When running DPMTA (version 3.1.2p3) the instabilities appear for the same P values and octree heights. We have also compared in [START_REF] Fortin | Multipole-to-local operator in the Fast Multipole Method: comparison of FFT, rotations and BLAS improvements[END_REF] the runtime of the FFT improvement in DPMTA (with fixed interaction list with ws = 1) and in our FMB code: similar CPU times were measured which validates the efficiency of our implementation.

With double height M2L kernel, the maximum degree of the O k j is 2P : the sizes of the arrays used in the FFT and in the point-wise product are thus roughly multiplicated by 4 (possibly thanks to zero-padding). But since O k j have now degrees up to 2P , the range of magnitude is even greater and they become unstable for even lower values of P : this imposes even more the block version. As in the non-block version, the instabilities in the block FFT appear however for lower values of P than in the single height case. Finally as for operation counts, we focus in the block FFT on the large-grain convolution which is the most time-consuming part: the ratio between single and double height kernels is here 8, which clearly does not favor the double height kernel compared to the ratio of 6 in classic M2L computation5 .

Rotations

The use of rotations has already been introduced in several articles: see [START_REF] Cheng | A fast adaptive multipole algorithm in three dimensions[END_REF], [START_REF] Choi | Rapid and stable determination of rotation matrices between spherical harmonics by direct recursion[END_REF], [START_REF] Greengard | A new version of the fast multipole method for the Laplace equation in three dimensions[END_REF], [START_REF] Gumerov | Recursions for the computation of multipole translation and rotation coefficients for the 3-D Helmholtz equation[END_REF] and [START_REF] White | Rotating around the quartic angular momentum barrier in fast multipole method calculations[END_REF]. This improvement allows the computation of the M2L operator6 in O(P 3 ), against O(P 4 ) for their classic version. We have chosen to use the formulae detailed by Gumerov and Duraiswami in [START_REF] Gumerov | Recursions for the computation of multipole translation and rotation coefficients for the 3-D Helmholtz equation[END_REF] since: they use the same definition for spherical harmonics, they use symmetries to speed up the computation and they focus only on the needed rotations. Moreover the recurrence is performed on real numbers, and not complex ones, and is simple to initiate. However no proof is given on the numerical stability of the formulae used.

The improvement in the use of rotations for M2L operator is based on the fact that the cost of an M2L operation performed along the z axis is O(P 3 ) against O(P 4 ) for general M2L. This is due to a property of the associated Legendre functions and this is valid for both single and double height M2L kernels. For a general M2L operation whose M2L vector is not aligned with the z axis, we have first to rotate the Cartesian coordinate system so that so that the M2L vector is along the z axis, then we perform the M2L operation, and finally we rotate back the coordinate system. The whole procedure is worthwhile since the first rotation on multipole expansion and the second on local expansion are both performed with O(P 3 ) operations. These rotations are performed thanks to rotation coefficients applied to the spherical harmonics and precomputed at each level of the octree during the downward pass. More details about the implementation of this scheme can be found in [START_REF] Fortin | Multipole-to-local operator in the Fast Multipole Method: comparison of FFT, rotations and BLAS improvements[END_REF].

The numerical stability of this scheme has been checked through whole FMM computations: no difference was detected in [START_REF] Fortin | Multipole-to-local operator in the Fast Multipole Method: comparison of FFT, rotations and BLAS improvements[END_REF] between the classic M2L scheme and the M2L scheme with rotations, for values of P up to 30, for both single and double height kernels and for several octree heights. When considering detailed operation count the ratio between single and double heights for M2L kernel appears to be roughly 20 7 which favors the double height kernel compared to the ratio of 6 of the classic M2L implementation (see footnote 5). Moreover, the memory requirements are very low since the number of rotation coefficient arrays equals the number of M2L transfer functions (316 with ws = 1).

In [START_REF] Fortin | Multipole-to-local operator in the Fast Multipole Method: comparison of FFT, rotations and BLAS improvements[END_REF] we have also discussed the possibility of introducing the BLAS routines (presented in section 3) in order to speed up the M2L computation with rotations, but no efficient solution was found: while the rotations of multipole and local expansions cannot be written as matrixvector products, the matrix-vector product corresponding to the M2L operation along the z axis computes too many useless terms and cannot be extended to matrix-matrix products.

Plane waves

Another computation scheme for the M2L operator is based on plane wave expansions. It has been presented in [START_REF] Greengard | A new version of the fast multipole method for the Laplace equation in three dimensions[END_REF] and improved in [START_REF] Cheng | A fast adaptive multipole algorithm in three dimensions[END_REF]. We now briefly present this scheme, and we refer the reader to these articles for more details.

The plane wave scheme is based on the following rewriting of the potential between two points X = (x, y, z) and X = (x , y , z ), with z > z :

1 X -X = 1 2π ∞ 0 e -λ(z-z )
2π 0 e iλ((x-x ) cos α+(y-y ) sin α) dα dλ,

where the double integral is estimated by quadrature formulae that depends on the required precision. Four precisions are used in [START_REF] Cheng | A fast adaptive multipole algorithm in three dimensions[END_REF] (see section 4.6 for details), and new quadrature formulae have to be written for each new precision. These precisions correspond to theoretical error bounds ensured by the FMM. More details about the quadrature formulae can be found in [START_REF] Cheng | Nonlinear optimization, quadrature, and interpolation[END_REF] and [START_REF] Yarvin | Generalized gaussian quadratures and singular value decompositions of integral operators[END_REF].

In order to compute a given M2L operation, the source multipole expansion is first converted into a plane wave expansion thanks to the quadrature formulae. The plane wave expansion is then translated from the center of the multipole expansion to the center of the target local expansion, where it is finally converted into a local expansion. This is however valid only if the error introduced by the plane wave expansions corresponds to the error due to multipole and local expansions. The number of quadrature nodes and plane waves is thus tied up to P and P 2 (see [START_REF] Cheng | A fast adaptive multipole algorithm in three dimensions[END_REF]). The translation of the plane wave expansion is then performed in O(P 2 ), which justifies the interest of this computation scheme, and the total operation count for M2L operation grows like O(P 3 ) + O(P 2 ) (due to conversions and possible additional rotations).

In summary of section 2, we have shown that the single height M2L kernel, contrary to the double height one, respects the FMM error bound which justifies its implementation. Both block FFT and rotations schemes reduce the operation count, but the use of the block FFT, contrary to the use of rotations, requires important extra memory storage and may result in unpredictable numerical instabilities. In addition, the plane wave scheme is stable for all values of P and presents the lowest operation count. We will now propose an alternative to speed up the M2L computation.

Multipole to Local operator: implementation with BLAS routines

The BLAS (Basic Linear Algebra Subprograms) (see [START_REF] Dongarra | A set of level 3 Basic Linear Algebra Subprograms[END_REF] [27] [START_REF] Lawson | Basic Linear Algebra Subprograms for FORTRAN usage[END_REF]) are a standard interface for some usual linear algebra operations such as a dot product of 2 vectors (level 1 BLAS), a matrixvector product (level 2 BLAS) or a matrix-matrix product (level 3 BLAS). Optimized for the pipelines of the floating point units and for the hierarchical memory of the computer, the BLAS routines offer an efficient implementation of these operations, and the higher the level of the BLAS used, the higher the speedup they reach. BLAS routines have already been used for hierarchical O(N ) N-body algorithms by Hu and Johnsson [START_REF] Hu | Implementing O(N) N-body algorithms efficiently in data-parallel languages[END_REF] with Anderson's method [START_REF] Anderson | An implementation of the fast multipole method without multipoles[END_REF] which uses different expansions than the FMM, and hence translation/conversion operators, but has the same algorithm for the upward and downward passes. Here we propose the first BLAS implementation for the M2L operator of the FMM for both single and double height kernel. In order to achieve the highest efficiency, we also detail a scheme with recopies that enables to use level 3 BLAS routines, and we show how to avoid these recopies thanks to new data storage in memory.

While the original FMM formulae of Greengard & Rokhlin (see [START_REF] Greengard | The Rapid Evaluation of Potential Fields in Particle Systems[END_REF]) do not allow a rewriting of their corresponding operators (M2M, M2L or L2L) as matrix-vector products, this can be easily done with the simpler formulae of [START_REF] Elliott | Multipole algorithms for molecular dynamics simulation on high performance computers[END_REF], [START_REF] Epton | Multipole translation theory for the three-dimensional Laplace and Helmholtz equations[END_REF] or [START_REF] White | Derivation and efficient implementation of the fast multipole method[END_REF]. The full FMM algorithm has also been rewritten in terms of matrix operations in [START_REF] Sun | A matrix version of the fast multipole method[END_REF]: here, we only focus on the rewriting of the M2L operator as a matrix-vector product.

As we compute only terms with positive orders in the local expansion, we can write each M2L operation as the following matrix-vector product, here written for P = 2 with single height M2L kernel:

O 0 0 O 1 1 O 0 1 O -1 1 O 2 2 O 1 2 O 0 2 O -1 2 O -2 2 O 0 1 O 1 2 O 0 2 O -1 2 0 0 0 0 0 O -1 1 O 0 2 O -1 2 O -2 2 0 0 0 0 0 O 0 2 0 0 0 0 0 0 0 0 O -1 2 0 0 0 0 0 0 0 0 O -2 2 0 0 0 0 0 0 0 0 3
M 0 0 M -1 1 M 0 1 M 1 1 M -2 2 M -1 2 M 0 2 M 1 2 M 2
The computed vector is named local vector. The matrix is named M2Ltransfer matrix and denoted T M 2L . With double height M2L kernel the matrix is dense: terms O k j with j > P do not vanish as in the single height case. The building of this matrix, detailed in [START_REF] Fortin | Multipole-to-local operator in the Fast Multipole Method: comparison of FFT, rotations and BLAS improvements[END_REF], is rather straightforward. The vector used is named multipole vector: terms with negative orders have to be stored.

Implementation with level 2 BLAS routines

Our matrices are stored in column storage mode, and we use the transfer matrix in its transposed form: it is indeed generally more efficient to compute C ← A T .B in row storage mode for A, than C ← A.B with column storage mode for A; indeed, when considering the BLAS implementation (see [START_REF]Fortran 77 reference implementation source code of the blas from netlib[END_REF] for a default implementation), the first solution leads to less writings, with however more readings, than the second one.

Double height M2L kernel

The matrix T M 2L is dense: the use of level 2 BLAS ZGEM V routine is therefore obvious, and this routine will automatically optimize the computation of the matrix-vector product to the underlying superscalar architecture. In order to differentiate it from the other BLAS method that will be used later, we name it full blas.

Single height M2L kernel

The matrix T M 2L is now sparse: we therefore have to split the sparse matrix-vector product in several dense block products.

We refer the reader to the BLAS literature (see [START_REF] Dongarra | Numerical Linear Algebra for High-Performance Computers (Software, Environments, Tools)[END_REF] for example) for the underlying techniques used to fill at best the pipelines of the floating point execution units in order to reach peak performance of the processor. These techniques are mainly: loop ordering for best spatial and temporal locality among data, loop blocking in order to provide maximum cache reuse, temporary copies in local arrays for problems with "leading dimensions" of the matrices (see BLAS routine interfaces), loop unrolling, register blocking, data prefetch, etc.

In [START_REF] Ling | GEMM-based level 3 BLAS: high-performance model implementations and performance evaluation benchmark[END_REF], it has been shown that level 3 BLAS routines can efficiently be implemented only with the GEM M level 3 BLAS routine and a few level 2 BLAS routines. For portability purpose, as well as for simplicity, we prefer to adopt such approach. We will thus decompose the sparse matrix-vector product corresponding to M2L operator in several ZGEM V block products in the most efficient way. ZGEM V routine is used here since we treat a matrix-vector product, but in section 3.2 we will see how to use matrix-matrix products, and ZGEM M routine will then be used as in [START_REF] Ling | GEMM-based level 3 BLAS: high-performance model implementations and performance evaluation benchmark[END_REF]. All the optimizations recalled above will be used but left as much as possible to the underlying (and machine dependent) ZGEM V /ZGEM M BLAS routine called. We point out now one important fact: in our implementation of the FMM, we are free in the memory storage of the blocks of T M 2L since its construction is precomputed at each level of the octree in the downward pass of the FMM. Therefore, most of problems that arise when considering the leading dimension of a matrix will be irrelevant with our blocks. With these principles in mind, we will present an efficient block decomposition, named block blas, for the matrix T M 2L with single height M2L kernel.

In [START_REF] Ling | GEMM-based level 3 BLAS: high-performance model implementations and performance evaluation benchmark[END_REF], triangular matrix-matrix product was performed with ZGEM M routines thanks to a decomposition of the triangular matrix in strips. First, we proceed similarly using strips in the horizontal direction since the transfer matrix is stored by rows, and since this leads to several traversals of the multipole vector (one for each strip) against one single traversal on the local vector. In other words, with strips in the row direction, the resulting blocks of the local vector are updated one after the other, whereas the corresponding data of the multipole vector may be reloaded several times during the block matrix-vector product. As the local vector is traversed for updating (i.e. both reading and writing) and the multipole expansion for reading, this is more efficient than vertical strips which lead to the contrary. In practice, our strips strictly 7 respect the underlying structure of T M 2L as pictured on figure 3(a): each strip is separately stored and separately computed with one dedicated call to ZGEM V routine. One obvious problem is that, as P grows, the first strips are too thin and too long while the last ones are too thick and too short: this can prevent the BLAS routine to internally decompose those strips according to the hierarchical memory of the computer.

But the underlying structure of T M 2L can be considered as recursive for a given P . Indeed when isolating the "biggest upper left square" (in the number of underlying sub-blocks) submatrix, with half of the sub-blocks in both the row and the column dimensions, as pictured in figures 3(b) and 3(c), we are left with one sub-matrix "on the right" and another one "below" whose sub-blocks are disposed in the same way as T M 2L for P 2 . Detailed expression of this recursion can be found in [START_REF] Fortin | Multipole-to-local operator in the Fast Multipole Method: comparison of FFT, rotations and BLAS improvements[END_REF]. 7 Another decomposition that uses additional zeros to obtain strips with optimal number of rows is discussed in [START_REF] Fortin | Multipole-to-local operator in the Fast Multipole Method: comparison of FFT, rotations and BLAS improvements[END_REF]. The "biggest upper left square" sub-matrix of T M 2L for P = 5 is pictured with the darker gray value. (c) In the remaining sub-matrices "on the right" and "below", whose sub-blocks are disposed in the same way as T M 2L for P = 2 (see equation ( 6)), we can recursively isolate a "biggest upper left square" sub-matrix.

Hence we use as much as possible dense matrices (i.e. the upper left sub-matrices) which are efficiently treated by ZGEM V routine. With small values for P , which corresponds also to the terminal cases of the recursion, it may be not worth isolating this dense matrix: in these cases we use our strips. The inefficiency due to the shape of the first and the last strips is therefore minimized since we reserve the strips for the terminal cases. The terminal value for the recursion, hereafter named s block , will have to be determined experimentally.

Level 3 BLAS routines

Since matrix-matrix product requires O(N 2 ) memory storage relatively to O(N 3 ) operation count, it is easier to overlap memory latency with computation of the floating point execution units with level 3 BLAS than with level 2 BLAS, and thus to reach peak performance of the processor. We will therefore try to group multiple M2L operations in one single matrix-matrix product.

During the downward pass, at a given level of the octree, all M2L operations that have the same M2L vector (see proposition 1) share the same M2L transfer matrix. When considering all pairs of multipole and local expansions that share the same M2L vector, it is thus possible to concatenate all their multipole vectors as columns of one single multipole matrix M M . The local vectors are also concatenated according to the same order to form one single local matrix M L . Then the matrix-matrix product M L = T M 2L . M M computes the corresponding M2L operations all at once as described in figure 4. The concatenation is easily achieved thanks to the column storage of our matrices.

First, we will roughly consider that the best efficiency is obtained with level 3 BLAS when the maximum number of multipole and local expansions are concatenated each time (see section 3.2.3 for revisions of this assertion). In the following, we will thus see how to concatenate the maximum number of multipole and local expansions for one given M2L transfer matrix, and we will then present the implementation of this matrix-matrix product thanks to level 3 BLAS calls.

With free-space boundary conditions8 (FBC), where all the space outside of the computational box is considered as empty, the cells located at the boundaries of the computational box, in each level, have an incomplete interaction list with less than 189 members for a well-separateness criterion ws = 1 (see section 2.1). In order to ease their computation, and for efficiency purpose, these cells with incomplete interaction list are computed separately when using level 3 BLAS for M2L.

¡ ¡ ¡ ¡ ¡ ¡ ¡ ¢ ¢ ¢ ¢ ¢ ¢ ¢ £ £ £ £ £ £ £ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¦ §¦ ¦ §¦ ¦ §¦ ¦ §¦ ¦ §¦ ¦ §¦ ¦ §¦ § § § § § § § ¨© §© § § § § § § § § § § § § § § § § § § § § § § § § § § § § § § § § § = § § § § § § § § § § ! §! §! ! §! §! " §" " §" # §# # §# $ §$ $ §$ % §% % §% & §& & §& ' §' ' §' ( §( ( §( ) §) ) §) 0 §0 0 §0

Computing the local expansions of cells with incomplete interaction list

With free-space boundary conditions, the "cells with incomplete interaction list" are all the cells whose parent have at least one neighbor that is outside of the octree. In order to treat all the M2L operations for the local expansions of such cells, we use copies of the multipole and local vectors. The M2L computation is thus performed as M L ← T M 2L . M M + M L and the columns of the local matrix are then recopied in the original local expansions. For a given ws value, the M2L vectors of the interaction list of a given cell are determined according to the type of child of the cell which describes the location of the cell center relatively to the center of its father. In 3D, there is eight different possible type of child. Some M2L vectors, but not all of them, are shared by all the types of child. We first loop on each possible M2L vector and secondly, inside this first loop, we loop on each type of child that matches the current M2L vector. This has been prefered to the opposite ordering of these two loops since it leads to a greater number of vectors treated per matrix-matrix product.

Computing the local expansions of cells with complete interaction list

All cells with complete interaction list have the same interaction list size, and all cells of the same type of child share exactly the same M2L vectors. This regularity allows well suited algorithms.

The most simple uses each time recopies as in [START_REF] Hu | Implementing O(N) N-body algorithms efficiently in data-parallel languages[END_REF] in order to treat altogether the maximum number of pairs of multipole and local expansions. Contrary to the cells with incomplete interaction list, we first loop on each type of child and secondly loop on each M2L vector corresponding to the current type of child: see algorithm 1.

Algorithm 1 Scheme with recopies. V M 2L denotes the set of all possible M2L vectors, T C the set of all different types of child and T M 2L (v) the M2L transfer matrix corresponding to the M2L vector v.

1: for all t ∈ T C do

2:

Copy all local vectors in M L ;

3:

for all v ∈ V M 2L corresponding to t do 4:
Copy the corresponding multipole vectors in M M ;

5:

Perform M L = T M 2L .M M ; 6:
end for

7:

Copy back all the local vectors ; 8: end for Thanks to the regularity, the recopies of the local expansions can be here performed out of the second loop, and the opposite order of the loops would lead to more copies for the local vectors while having the same number of vectors treated each time.

It is however possible to avoid the additional cost of recopies thanks to special data storages of our multipole and local vectors which are obtained through a "rearrangement" of these vectors in memory performed before the downward pass. In the first data storage, named row storage and illustrated in figure 5, we store consecutively in memory all the local vectors that belong to cells of the same type of child along a row in one given dimension of the 3D space. This is done for each row of the level, and the same storage is also done for the multipole vectors. Note however that for local vectors, only cells with complete interaction list have their local expansions rearranged in rows, while for the multipole expansions the rearrangement has to be performed for all cells of the level. With such data storage, we can call a level 3 BLAS routine starting at the local vector of the first cell of each row, with the corresponding M2L transfer matrix and the corresponding multipole vector. It can be noted that Morton ordering (see section 2.1) is used here in order to access cells according to the coordinates of their center. With FBC, at level l and for a given type of child, the local expansions are hence rearranged in (2 l-1 -2.ws) 2 rows of size 2 l-1 -2.ws. The size of the rows corresponds to the number of local vectors treated per level 3 BLAS call; this might be not enough to obtain the best efficiency: for example there is only 6 columns at level 4. That is why we propose a second data storage, named slice storage, that stores consecutively the rows in memory in order to form a slice and thus enables several rows to be treated with one single level 3 BLAS call. In order to have a correct correspondence between slices of local expansions and slices of multipole expansions, we have to insert blank boxes between rows of local expansions. These blank boxes correspond to cells not belonging to the octree whose local vectors are computed uselessly. They are actually required to compensate the absence of the cells with incomplete interaction list which are treated separately. With ws = 1, there is 1 blank box at the beginning and 1 at the end of each row, as pictured in figure 6. We can however skip the first blank box of the slice (that is to say, the first blank box of the first row), and also the last one.

With FBC, at level l, and for each type of child, the local expansions are rearranged in 2 [START_REF] Hu | Implementing O(N) N-body algorithms efficiently in data-parallel languages[END_REF] 29 @ 2@ 2@ @ 2@ 2@ @ 2@ 2@ @ 2@ 2@ A 2A 2A The rows of the row data storage mode are concatenated with blank boxes (marked with stripes) in between: the matrix-matrix product now involves multipole and local matrices with 6 columns, against 2 columns with row storage (see figure 5).

A 2A 2A A 2A 2A A 2A 2A B 2B B 2B B 2B B 2B C 2C C 2C C 2C C 2C D 2D 2D D 2D 2D D 2D 2D D 2D 2D E 2E 2E E 2E 2E E 2E 2E E 2E 2E F 2F F 2F F 2F F 2F G 2G G 2G G 2G G 2G H 2H H 2H H 2H H 2H H 2H I 2I I 2I I 2I I 2I I 2I P 2P 2P P 2P 2P P 2P 2P P 2P 2P P 2P 2P Q 2Q 2Q Q 2Q 2Q Q 2Q 2Q Q 2Q 2Q Q 2Q 2Q R 2R 2R R 2R 2R R 2R 2R R 2R 2R R 2R 2R S 2S 2S S 2S 2S S 2S 2S S 2S 2S S 2S 2S T 2T T 2T T 2T T 2T T 2T U 2U U 2U U 2U U 2U U 2U V 2V V 2V V 2V V 2V V 2V V 2V V 2V V 2V V 2V V 2V V 2V V 2V V 2V V 2V W W W W W W W W W W W W W W X 2X X 2X X 2X X 2X X 2X X 2X X 2X X 2X X 2X X 2X X 2X X 2X X 2X X 2X Y Y Y Y Y Y Y Y Y Y Y Y Y Y = = =
slices of size (2 l-1 -2.ws) 2 + 2 * (2 l-1 -2.ws) -2. If slice storage allows a better efficiency with the level 3 BLAS, its clear drawback is that it also performs useless extra work due to the blank boxes. This can even be extended to a third data storage, named level storage, resulting in concatenation of slices so that the whole level can be computed in one single BLAS call: in this case rows of blank boxes have to be added between each slice.

Implementation with level 3 BLAS routines

As in section 3.1.1, the dense matrix-matrix product for double height M2L kernel can be directly implemented with one single level 3 BLAS call: the ZGEM M routine.

For single height M2L kernel the decomposition proposed in section 3.1.2 is also valid with matrix-matrix products when replacing ZGEM V routine by ZGEM M one. However, we have to face a new constraint since we have several level 3 BLAS calls per matrix-matrix product: if all the columns of the multipole or local matrix cannot be stored in a level of the hierarchical memory (cache L1, L2 or even L3), or need more memory pages than the TLB (Translation Lookaside Buffer) can address, the whole matrix would have to be reloaded at each BLAS call. A matrix with N bExp expansions is therefore split in sub-matrices with a constant number of columns, namely N bExp max , and the same number of rows as the original matrix as described in figure 7: each sub-matrix is then treated separately. This decomposition of the matrix-matrix product is hereafter named "N bExp max decomposition".

Optimal values for N bExp max will be experimentally determined in section 3.3.1.

Tests and comparisons

Performance tests have been performed in order to validate the BLAS implementation, along with its parameters, that leads to the fastest M2L computation. KB, L2 cache size: 1.41 MB, L3 cache size: 32 MB, and 8 GB of memory), both at LaBRI9 . The number of particles used for the simulation is not usually precised since the runtimes here measured depend only on the height of the unifom octree and on P , but of course the height used implies a range in the number of particles that balances the near field and the far field computations.

Decomposition for single height M2L kernel

In this section, we will see how we have established, for each architecture (Power3 or Power4), the best N bExp max value, as well as the best s block value for the block blas routine (see section 3.1.2). As use of level 3 BLAS calls always improves the performances over level 2 ones (see [START_REF] Fortin | Multipole-to-local operator in the Fast Multipole Method: comparison of FFT, rotations and BLAS improvements[END_REF]), these values have been determined when using matrix-matrix products. First we need to determine the optimal N bExp max values. These are clearly machine dependent, and they depend on P too since P determines the number of rows in the multipole matrix and in the local matrix. Depending on the level 3 BLAS scheme used and on the level in the octree, the number of multipole and local expansions concatenated, denoted by N bExp, differs significantly. For performance tests with different N bExp max values, we have only considered N bExp values for cells with complete interaction list since this corresponds to the majority of the matrix-matrix products performed. The table 1 shows N bExp values according to our different schemes for cells with complete interaction list.

In practice, tests for optimal N bExp max values will be performed, according to P , for one single "big enough" N bExp value. Indeed, as confirmed by more complete tests, higher N bExp values will have the same optimal N bExp max values and lower ones will not need such decomposition. Generally 2744 or 958 are sufficient: see for example figure 8. Moreover, results have shown that N bExp max is mainly independent from s block : when splitting both multipole and local matrices according to N bExp max , the sub-matrices size depends indeed only on P (for the number of rows) and N bExp max (for the number of columns): see figure 7.

Besides, the optimal s block values, which depend on P and on the machine used, were searched between 1 and P . For low values of P , the best s block values equal P , which means it is better to use only strips for the decomposition. But for higher values of P , the best s block values were lower than P which justifies the recursive decomposition: see figure 8.

Finally, with these optimal s block values, the relevance of the N bExp max decomposition is shown on figure 8 for the scheme with recopies. With growing octree height, the numbers of expansions N bExp treated with one single matrix-matrix product increases, and the gain offered by the N bExp max decomposition increases thus too. With row and slice storages, N bExp max decomposition will be likewise relevant when the N bExp value exceeds the optimal N bExp max value. Gain in percentage. Figure 8: Gain in percentage for the full downward pass CPU times offered by a computation with N bExp max decomposition relative to one without N bExp max decomposition. Tests performed on IBM Power4, using the block blas routine (with optimal s block ) and the scheme with recopies for an octree height equal to 5 (N bExp = 2744) and 6 (N bExp = 27000). The corresponding values for s block and N bExp max are also given.

Recopies and special data storages

In order to compare the scheme with recopies with row and slice storages, CPU times have been measured on the full downward pass of the FMM. As mentioned in [START_REF] Hu | Implementing O(N) N-body algorithms efficiently in data-parallel languages[END_REF], the cost of copying vectors relatively to the cost of the matrix-matrix product decreases for growing values of P : copies of multipoles and local expansions are indeed performed in O(N bExp × P 2 ) while the matrix-matrix product requires O(N bExp × P 4 ) operations. This is more obvious with double height kernel than with single height kernel since the amount of computation for the matrix-matrix product is much more costly with double height, while the cost of recopies is the same. Figure 9 shows the gain of row and slice storages relative to the scheme with recopies according to different values of P for an octree height of 6: row and slice storages offer thus better gains relative to recopies for low values of P , and these gains are always higher for single height kernel. These gains are also influenced by the height of the octree: with growing heights of the octree, the number of expansions treated with one matrix-matrix product increases for row and slice storages (see table 1) and the proportion of blank boxes decreases for slice storage. And for low heights of the octree, the use of recopies may sometimes be faster as in figure 9.

Slice storage is here more efficient than row storage because of the too small N bExp for row storage at this height: for greater heights, row storage is a little bit faster. Moreover slice storage needs very long consecutive memory areas that cannot be allocated for too high values of P or too high heights when memory allocation with row storage may succeed. For these reasons we prefer row storage, and since the N bExp values enabled by slice storage seem to be high enough, level storage has also been discarded.

At last, as mentioned in sections 2.3 and 2.4, no BLAS routine can be applied to the FFT and Rotations improvements, and attempts to apply the scheme with recopies or row data storage have failed in improving the performances of these computations. Figure 9: Gain in percentage offered by row and slice storages relative to the scheme with recopies for downard pass CPU times with an octree of height 6. Tests performed on IBM Power4. Now that we have determined the best BLAS implementations, namely block blas and full blas (depending on the M2L kernel height), with row data storage, we are able to practically compare them with the other M2L improvements.

Comparison of the M2L improvements

Thanks to the implementation of the schemes based on FFT, rotations and BLAS routines in our FMB code, we can now precisely compare all these improvements of the M2L operator. The comparison with the plane wave scheme will be done afterwards thanks to the FMMPART3D code.

Memory requirements

We first compare the memory requirements of the different schemes since it may be the first choice criterion. We here focus only on the memory used for the expansions: the memory used for the particles remains indeed unchanged when computing M2L operator differently. In a nutshell, M2L computation with rotations has very low extra memory needs (see section 2.4). BLAS computation needs multipole vectors, whose size is roughly twice the size of the multipole expansions, and extra memory for the M2L transfer matrices, especially with double height kernel since these matrices are dense in this case. slice data storage requires also extra memory for the blank boxes. At last M2L computation with block or non-block FFT uses bigger arrays for its multipole expansions and its M2L transfer functions: these arrays are roughly 4 times bigger with single height M2L kernel and 16 times bigger with double height kernel. Using more detailed theoretical estimations that can be found in [START_REF] Fortin | Multipole-to-local operator in the Fast Multipole Method: comparison of FFT, rotations and BLAS improvements[END_REF], we have plotted in figure 10, for each scheme, the ratio of its memory need to the classic M2L memory need, using an octree of height 6 and only multiples of the FFT block size.

As predicted, while the rotation requirements are almost unnoticeable, and the BLAS ones remain moderate, FFT extra memory appear as problematic, especially for double height M2L kernel. The same ratios apply for other FFT block sizes. Moreover the ratio of 2 for BLAS in the double height kernel with high values of P is due to the dense M2L transfer matrices: for greater octree heights, this ratio remains close to 1.5 since the number of M2L functions is constant while the number of cells in the octree grows exponentially. At last the additional cost of the blank boxes in slice storage over row storage is very small.

Single height

CPU times for single height M2L kernel

Figure 11 compares the different schemes for an octree of size 5. If the BLAS version outperforms the version with the classic M2L and the one with rotations, the block FFT is faster, especially for values of P greater than 10. However in this test, without rescaling of the particle coordinates, the FFT with block size 4 is unstable for P ≥ 14, and the one with block size 3 for P ≥ 23.

When focusing on low and medium precisions, we recall also that for growing heights of the octree, the BLAS computation with row or slice storages becomes more efficient since the length of rows and slices increases at the leaf level which represents most of the runtime: the number of expansions treated in one level 3 BLAS call thus increases, as well as the efficiency of the BLAS routines. Figure 12 shows that BLAS computation scheme is clearly competitive with the FFT scheme with block size 4, in an octree of height 7 and for low and medium precisions. Moreover the memory requirements of the FFT improvement are problematic here, since on figure 12, the brutal increase in runtime for FFT at P = 10 is due to swapping despite the 16 GB of memory available.

CPU times for double height M2L kernel

With double height kernel, the use of BLAS clearly outperforms all other methods. Figure 13 compares the different schemes for an octree of height 5. For the BLAS, we have plotted only the row data storage, but slice storage and the scheme with recopies have similar performances. For the FFT, the block version with blocks of size 4 have been plotted; in this test, it is only stable for P ≤ 8, and more stable FFT with lower block sizes are slower. Moreover, we have stop the tests at P = 19 since on the IBM Power3, the 2GB of memory were insufficient for all FFT block sizes with P > 19. However it must be noticed that at this octree height, the rotation scheme becomes faster than all the BLAS schemes for P ≥ 23: this is of course the aftermath of the lower operation count for the rotation scheme. Finally, as for single height kernel, the BLAS computation with row and slice storages become even more efficient for greater heights of the octree.

Computational efficiency

In order to illustrate the efficiency of the BLAS versions that makes them faster than the other schemes, we present the table 2 that shows the MFLOPS (Millions of Floating Point Operations per Seconds) when computing M2L. The results are shown as percentages of the peak performance of the IBM Power3 used: 1500 MFLOPS. They have been computed with the Hardware Performance Monitor (HPM) Toolkit (version 2.4.3) as the average MFLOPS rate over all M2L computations of a full FMM computation. For level 3 BLAS, we have used block blas and full blas routines for respectively single and double height kernels, with an octree of height 5 and with recopies (the additional cost of copying the expansions is not considered here since we focus on the BLAS routine efficiency). It takes into account all cells, i.e. with and without complete interaction list. The height of the octree does not matter for classic M2L, rotations and FFT. These results can however not be directly compared from one row to the other: we recall here that the operations count differs! However it clearly illustrates how efficiently the processor is used in the different implementations and why BLAS computations outperform the other schemes. This table also indicates that our decomposition of the matrix-matrix product for single height kernel in the block blas routine, though satisfactory, is not optimal when compared to the full blas routine efficiency of the double height kernel, especially for low values of P .

Single versus double height kernels

We have already seen (see theorem 4) that a sharp error bound has been theoretically proven only for the single M2L height kernel. The double height M2L kernel is certainly more precise, and for a given precision, it would therefore requires a lower value for P than the single height one. But since no precise error bound is available for double height kernel, we cannot a priori know what this lower P will be.

That is why we compare here CPU times with practical accuracies for both kernel heights. As already exposed (see for example [START_REF] White | Derivation and efficient implementation of the fast multipole method[END_REF] for double height M2L kernel, and [START_REF] Rankin | Efficient parallel implementations of multipole based N-body algorithms[END_REF] for single height kernel), these practical accuracies are better than the theoretical ones which correspond to worstcase errors. These worst-case errors are indeed obtained with spherical regions, while we use in fact smaller cubical cells because of the octree. Moreover, we use relative error instead of absolute one so that the results are problem independent, and these relative errors are computed for the potential and not for the force components because these latters can be almost null. Since some discontinuities appear in the potential error for particles crossing cell boundaries (see [START_REF] Rankin | Efficient parallel implementations of multipole based N-body algorithms[END_REF] ), we choose to study the more stable RMS average error instead of the maximum error. This RMS average error ε rms is computed as follows:

ε rms = 1 N N i=1 Φ Dir (X i ) -Φ F M M (X i ) Φ Dir (X i ) 2 .
We consider here a gravitational potential computed for an uniform distribution with 100000 bodies and an octree of height 4, which results in 25 bodies per leaf in the mean. We recall that when the number of bodies per leaf increases, the part of the near field (directly computed) in the potential becomes higher and the potential becomes thus more precise. Figure 14 presents, for each M2L scheme, the tradeoff between practical error and CPU times (downward pass only) with both single and double height kernels. The scales are logarithmic, and the values of P plotted for double height kernel range from 3 to 14 with a step of 1, while for single height kernel they range from 3 to 29 with step 2: P = 14 in double height and P = 29 in single height correspond indeed to the same accuracy (precisely the first accuracy below 1.0 × 10 -9 in our simulation).

As far as classic M2L is concerned, the single height kernel is more efficient for low precisions and the double height one for high precisions. And as theoretically justified in sections 2.3, 2.4 and 4.4, the rotation and BLAS computations generally favor the double height kernel, except for the low values of P where the single height kernel is more efficient. On the other hand, the FFT improvement is generally more efficient with single height kernel.

Therefore we still have to compare the best of each scheme: this is done in figure 15 where we select the best kernel height for different P values and present all the schemes on the same plot. The BLAS version is then clearly more efficient than all the other schemes. By examples, for a practical error below 1.0 × 10 -6 the gain of the BLAS over the FFT (respectively the rotations) is 35% (respectively 62%). This fully validates the relevance of our new BLAS version compared to these M2L improvements.

Comparison with plane waves

In order to compare our BLAS version with the scheme based on plane wave expansions, we will study the CPU times of both our FMB code and the FMMPART3D code (version 1.0) according to the practical accuracies. This comparison has been done on a Linux PC (Intel Pentium 4 processor at 2.6 Ghz, L1 cache size: 8 KB, L2 cache size: 512 KB, and 2 GB of memory). Our FMB code is compiled with gcc 4.0.4 and uses the ATLAS BLAS library (version 3.6.0) [START_REF]Atlas homepage[END_REF] [START_REF] Whaley | Minimizing development and maintenance costs in supporting persistently optimized BLAS[END_REF]. The FMMPART3D code is available as a precompiled library (Lahey/Fujitsu Fortran 95 compiler, version 6.2).

We start by outlining that the articles [START_REF] Greengard | A new version of the fast multipole method for the Laplace equation in three dimensions[END_REF] and [START_REF] Cheng | A fast adaptive multipole algorithm in three dimensions[END_REF] do not introduce the Outer and Inner functions (see section 2.1): their expressions of M2M, M2L and L2L operators are thus more complicated, and therefore more expensive than ours in terms of operation count. Rotations.

- As for previous comparisons, we only consider here uniform distributions of particles. Like in FMMPART3D, and in the articles [START_REF] Greengard | A new version of the fast multipole method for the Laplace equation in three dimensions[END_REF] and [START_REF] Cheng | A fast adaptive multipole algorithm in three dimensions[END_REF], we use here the L 2 norm as practical error measurement:

ε L2 = M i=1 (Φ Dir (X i ) -Φ F M M (X i )) 2 M i=1 Φ Dir (X i ) 2 .
The direct computation is then used to compute the exact potential Φ Dir of M = 1000 bodies (chosen randomly among the N bodies), and the error generated by FMMPART3D and by FMB for these M bodies is computed according to the L 2 norm.

The FMMPART3D code offers the four following theoretical precisions ensured by the plane wave scheme: 1.6 × 10 -3 (low), 1.3 × 10 -6 (medium), 1.1 × 10 -9 (high) and 1.0 × 10 -12 (very high).

The corresponding P values are presented in table 3. As usual with the FMM (see section 4.5), we notice in the following tests that the practical errors obtained with FMMPART3D are much lower than these theoretical error bounds.

For FMB, the values of P range from 1 to 28, with step 1, so that we cover the whole range of precisions available in FMMPART3D. The FMB code selects the best computation scheme with level 3 BLAS routines depending on the octree level: the row data storage mode when the level l in the octree is high enough (l ≥ 5), and otherwise the scheme with recopies (see section 3.3.2).

Moreover, some algorithmic differences between the two codes have to be precised.

• With FMB, we use single height and double height M2L kernels. The single height kernel is however presented here only up to P = 15 since, as noticed in section 4.5, the single height kernel is more efficient than the double height one only for the low values of P . In FMMPART3D and in the articles [START_REF] Greengard | A new version of the fast multipole method for the Laplace equation in three dimensions[END_REF] and [START_REF] Cheng | A fast adaptive multipole algorithm in three dimensions[END_REF], the choice of the kernel height is unprecised, which suggests a double height kernel (more casual and more common).

• As far as the octree height is concerned, this one is optimally set by the user in FMB in order to balance the near field and the far field computations. With growing values of P , the cost of the far field computation increases, and we may have to reduce the octree height. On the other hand, FMMPART3D is an adaptive FMM code, whose algorithm is presented in [START_REF] Cheng | A fast adaptive multipole algorithm in three dimensions[END_REF] and [START_REF] Carrier | A fast adaptive multipole algorithm for particle simulations[END_REF]. This algorithm imposes a maximal number of particles per leaf in the octree: this maximum value can not be tuned by the user in FMMPART3D.

• In FMMPART3D, the operations M2M and L2L are computed thanks to a scheme with rotations in O(P • We first notice with table 3 that for each precision available in FMMPART3D, the P value required to reach (or surpass) the corresponding practical precision (that is to say the error obtained in practice) with FMB is much lower than the P value used in FMMPART3D. Such a gap is probably due to a loss of the FMM accuracy because of the plane wave introduction. The error generated in practice with the plane wave expansions may be greater than the one due to the spherical harmonic expansions, even though the practical error of FMMPART3D remains of course below the theoretical error bounds. A decrease of this additional error due to the plane wave scheme could only be done at the expense of an increase in its operation count.

• We can also remark that FMB is more efficient, with respect to FMMPART3D, when both the potential and the force are computed (in comparison with the computation of the potential only). This can be explained by a better implementation of our direct computation for the force, or by a more efficient computation of the evaluation of the force from the local expansions. As this is not directly related to the M2L computation, we focus on the figures where only the potential is computed.

• As announced, the single height M2L kernel is more efficient than the double height kernel for the low values of P (concretely for 2 ≤ P ≤ 5 -6).

• Finally, if we compare the CPU times for a given precision and for the computation of the potential only, FMB is always better for the first precision (low), with for example a gain of 39% with P = 4 (single height M2L kernel) in the case of 1 million particles. For the two following precisions (medium and high), FMB is as fast as FMMPART3D, but slower for the highest precision (very high).

This perfectly illustrates the gain in performance of the BLAS computation scheme which keeps the operation count of the classic M2L computation but greatly reduces the underlying constant factor: for low and medium values of P the best gain is offered by the BLAS routines, and the gain in operation count offered by the plane wave scheme is crucial only for the highest values of P .

The breaks in the FMB plots correspond to decrements of the octree height. The octree heights being lower in the case of 100 000 particles, our computation scheme with BLAS is less efficient, but the same behaviour is observed. For the cases with 1 million and 3 million particles, the octree heights are the same (with however different decrements) and the differences in the plots are essentially due to the more important part of the direct computation with 3 million particles. Due to compiler and licence restrictions we could not run FMMPART3D on other architectures with more memory (like SMP nodes 10 ), but in this case we would have been able to run bigger test cases with greater octree height, where our BLAS scheme would be likely to be even more efficient (see section 3.3.2).

As exposed in the introduction of this paper, the targeted precisions are lower or equal to roughly 10 -7 . This corresponds to the range of P values where our FMB code is either faster than, or as fast as, the FMMPART3D code. By the way, this may explain why the plane wave scheme has not been used with Laplace equation since the article [START_REF] Cheng | A fast adaptive multipole algorithm in three dimensions[END_REF]: contrary to Maxwell and Helmholtz equations, where the P values are high, the values of P required in Laplace equation limit the gain offered by a the operation count in O(P 2 ) + O(P 3 ), and the underlying constant factor of this operation count is not negligeable in this case.

Conclusion

In this paper, we have presented an overall study of the most efficient implementation of the Fast Multipole Method for the serial computation of gravitational or electrostatic simulations of uniform distributions.

A detailed study of the error bound has lead us to two expressions of the M2L operator that converts a multipole expansion into a local expansion: while the double height M2L kernel is generally more efficient for medium and high precisions, the single height M2L kernel is faster for low precisions and is the only one to ensure a sharp error bound. For each M2L expression, we have efficiently implemented the block FFT and the rotation improvements. To our knowledge, this is the first implementation of the FFT enhancement for double height M2L kernel, and it has been shown that numerical instabilities remain even with the block FFT.

As an alternative, a BLAS (Basic Linear Algebra Subprograms) version has been proposed for the dense matrices of the double height kernel as well as for the sparse matrices of the single height one. A scheme with recopies has first made possible the use of level 3 BLAS resulting in impressive speedups. Special data storages for the expansions either by rows or by slices have then enabled us to avoid the additional cost of recopies, especially for low precisions.

Memory requirement estimations and CPU times from practical simulations have been used to precisely compare all these different schemes. It appears that the BLAS version and the FFT improvement with blocks are the most efficient ones. While the BLAS version is always faster in case of double height kernel, the block FFT is faster with single height kernel for high precisions. However the memory requirements of the block FFT as well as the remaining numerical instabilities, precisely for high precisions, limit severely the benefit it offers for runtime in this case. When comparing the tradeoff between practical accuracy and runtime for both M2L kernel heights, the BLAS version is then the most efficient, while introducing no numerical instabilities and low extra memory requirements.

Finally, we have also compared our BLAS computation scheme with the plane wave scheme thanks to the FMMPART3D code: in the range of precisions that interest us for Laplace equation in astrophysics and molecular dynamics, the gain offered by the BLAS routines is either greater than, or equal to, the one obtained with the plane wave scheme.

We have already extended this BLAS scheme to the adaptive version of the FMM [START_REF] Coulaud | High-performance BLAS formulation of the Adaptive Fast Multipole Method[END_REF] [39], and we are currently parallelizing it on shared and distributed memory architectures. This BLAS approach could also be extended to other potentials whose FMM operators can be written as matrix products.

where the L l n denote the local expansion terms, centered in X 2 , due to the unit charge located in Q, and obtained thanks to the M2L operator with single height M2L kernel applied to the multipole expansion terms M k j as

L l n = P -n j=0 j k=-j M k j O -l-k n+j (X 2 -X 1 ).

Lemma 1 (

 1 Symmetry among orders) Y -m l = Y m l where z denotes the complex conjugate of z ∈ C.

Figure 2 :

 2 Figure2: Logarithmic error for the potential (Φ) and the force coordinates (f x, f y, f z) according to P for single height M2L kernel with FFT with blocks of size 4. Tests performed on 10000 particles, uniformly distributed, with octrees of height 4 and 6: all particle coordinates are inside [0, 1.0]. The error plotted is the maximum absolute error over all the particles. Tests run on IBM Power3 and Power4 with double precision.

Figure 3 :

 3 Figure 3: block blas decomposition for P = 5. (a) Decomposition with only strips: the different gray values show the different strips. (b)The "biggest upper left square" sub-matrix of T M 2L for P = 5 is pictured with the darker gray value. (c) In the remaining sub-matrices "on the right" and "below", whose sub-blocks are disposed in the same way as T M 2L for P = 2 (see equation (6)), we can recursively isolate a "biggest upper left square" sub-matrix.

Figure 4 :

 4 Figure 4: Concatenation of multipole and local expansions in order to use level 3 BLAS. The 4 M2L operations that share the same M2L vector are computed together in one single matrix-matrix product: the corresponding multipole and local vectors are hence concatenated in one multipole matrix and in one local matrix, both with 4 columns. The transfer matrix represented here is for a single height M2L kernel.

Figure 5 :

 5 Figure 5: row storage (2D, ws = 1, level 3). The cells are indexed according to Morton ordering. Each of the 4 types of child has a different gray value. The expansions of the cells with same type of child along a given row are concatenated: the 4 M2L operations, whose M2L vectors are represented on the quadtree, can then be directly computed with 2 matrix-matrix products (level 3 BLAS) without recopies.

Figure 6 :

 6 Figure6: slice storage (2D, ws = 1, level 3): see also figure5. The rows of the row data storage mode are concatenated with blank boxes (marked with stripes) in between: the matrix-matrix product now involves multipole and local matrices with 6 columns, against 2 columns with row storage (see figure5).

  block blas routine (single height kernel).

  full blas routine (double height kernel).

Figure 10 :

 10 Figure 10: Memory requirements of the different M2L computation schemes for an height of 6 (ratios to classic M2L computation requirements).

Figure 11 :

 11 Figure 11: Logarithmic downward pass CPU times of the different M2L computation schemes, for an octree of height 5 and with single height M2L kernel. Tests are performed on an IBM Power3 with 2 GB of memory.

Figure 12 :

 12 Figure 12: Downward pass CPU times of the different M2L computation schemes, for an octree of height 7 and with single height M2L kernel. Tests are performed on an IBM Power3 with 16 GB of memory at CINES (Centre Informatique National de l'Enseignement Supérieur, Montpellier, France).

  RMS relative error for potential Downward pass (log) CPU times (seconds)

Figure 14 :

 14 Figure 14: Tradeoffs between practical accuracies and CPU times for single and double height kernels. The numerical instabilities of the FFT scheme with block size 4 prevent from reaching a precision below 10 -7 in these tests.

Figure 15 :

 15 Figure 15: Tradeoffs between practical accuracies and CPU times with the best kernel height for each accuracy. See figure 14.

Figure 16 :

 16 Figure 16: Comparison between FMMPART3D and FMB for an uniform distribution of 100 000 particles.

Figure 17 :

 17 Figure 17: Comparison between FMMPART3D and FMB for an uniform distribution of 1 million particles.

Figure 18 :

 18 Figure 18: Comparison between FMMPART3D and FMB for an uniform distribution of 3 million particles.

Table 1 :

 1 N bExp max decomposition (max denotes N bExp max ). This matrix-matrix product is computed as several matrix-matrix products with at most N bExp max columns in the multipole and local matrices. N bExp values according to our different schemes (with FBC, and ws = 1, see section 3.2.2).

	max

These tests have been performed either on one IBM Power3-II WH2+ (375 MHz, 1.5 GFLOPS, L1 cache size: 64 KB, L2 cache size: 4 MB, and 2 GB of memory) or on one IBM Power4+ (1454 MHz, ≈ 6 GFLOPS, L1 cache size: 64

Table 2 :

 2 Logarithmic downward pass CPU times of the different M2L computation schemes, for an octree of height 5 and with double height M2L kernel. Tests are performed on an IBM Power3 with 2 GB of memory. Computational efficiencies of the different M2L computation schemes.

	log time (seconds)	1000 10000 100000		classic FFT with block of size 4 Rotations Full blas with Row
		100
	PSfrag replacements
		10	4	6	8	10	12	14	16	18
						P
	Figure 13: Single height M2L kernel		Double height M2L kernel
				P = 7			P = 15		P = 7	P = 15
	classic			9.4 %			16.3 %		14.5 %	18.3 %
	rotations			5.4 %			7.8 %		8.4 %	10.9 %
	FFT with block size 4			4.5 %			13 %		12.4 %	18.4 %
	level 3 BLAS			46.4 %			74 %		85.9 %	89.2 %

Table 3 :

 3 3 ) (see section 2.4), whereas in FMB we use the classic computation in O(P 4 ). Nevertheless these operations are clearly minority in the total computation cost of the FMM for uniform distributions. For each of the four precisions available in FMMPART3D: P values used in FMMPART3D and P values in FMB that enable to reach (or surpass) the corresponding practical precision obtained with FMMPART3D.Figures 16, 17 and 18 present the tradeoff between CPU times and ε L2 error for both codes, and for uniform distributions of 100 000, 1 million, and 3 million particles. In each figure, the four dots for FMMPART3D correspond to the four available precisions. Several conclusions can be drawn from these figures.

	P	low	medium	high	very high
	FMMPART3D	10	19	29	40
	FMB (double height M2L kernel)	4	8	16	26-28
	FMB (single height M2L kernel)	4	14	-	-

The plane wave scheme has however been widely used for Maxwell and Helmholtz equations (see[START_REF] Darve | Efficient fast multipole method for low-frequency scattering[END_REF] for example) since the P values are much higher in these cases.

http://www.madmaxoptics.com/technology/products/FMMPART3D.html

The kernel name has been inspired by[START_REF] Elliott | Multipole algorithms for molecular dynamics simulation on high performance computers[END_REF].

The equation (C.23) in appendix C of[START_REF] Elliott | Multipole algorithms for molecular dynamics simulation on high performance computers[END_REF] is wrong because the composition of differential operators differs from the product of the derivatives.

For classic M2L computation with single height M2L kernel, we have thanks to lemma 1 an operation count like:[START_REF] Greengard | A fast algorithm for particle simulations[END_REF] 12 (P + 3)(P + 1)(P + 2) 2 . With double height kernel, we have: 1 2 (P + 2)(P + 1)[START_REF] Elliott | Multipole algorithms for molecular dynamics simulation on high performance computers[END_REF] . Therefore the ratio among the two heights is roughly

[START_REF] Singh | Load balancing and data locality in adaptive hierarchical N-body methods: Barnes-Hut, fast multipole, and radiosity[END_REF] This applies also to M2M and L2L operators.

Such specificity does not appear with "Periodic Boundary Conditions", where the computational box is periodically replicated in each dimension.

Laboratoire Bordelais de Recherche en Informatique, Talence, France.

On SMP nodes, vendor BLAS implementations are also likely to be faster than the ATLAS routines.

Acknowledgements

We are grateful to Guillaume Latu, ICPS-LSIIT (Laboratoire des Sciences de l'Image, de l'Informatique et de la Télédétection) and INRIA CALVI project, for several helpful discussions on the underlying techniques of BLAS and their integration in the FMM.

A Proof of the theorem 4.

Proof. With the notations and assumptions of section 2.2, we consider a multipole expansion, centered in X 1 , and a local expansion, centered in X 2 , to express the potential at Z due to the particle in Q. Since the condition R > r 1 + r 2 , which guarantees the convergence of both multipole and local expansions, implies (Q-X 1 )-(Z-X 2 ) < X 2 -X 1 , and since I l n (-X) = (-1) n I l n (X) the potential in Z writes, with the classical translation theorem (theorem 1):

Thanks to the Inner-to-Inner translation theorem (theorem 3), we obtain

We emphasize here that the equation ( 7) is just a rewriting of equation ( 4): that is why when troncating the series in equation [START_REF] Carr | Reconstruction and representation of 3D objects with radial basis functions[END_REF] for n > P , we obtain the same error bound as in proposition 2. Thus we have

After an inversion of the two finite summations on the degrees n and j and reindexing n-j → n and l + k → l we obtain

We remember here that I -l n (x) imposes -n ≤ l ≤ n . In the same way, I k j (x) imposes |k| ≤ j, that is to say j + k ≥ 0, and kj ≤ 0. Moreover,

Under these conditions and after reindexing -l → l and n → n and another inversion of summations on the degrees, we obtain

For all Q in B 1 , the multipole expansion terms M k j , ∀ (j, k) ∈ N × Z with 0 ≤ |k| ≤ j, being defined by M k j = (-1) j I k j (Q -X 1 ), we have thus proved that Φ P (Z) =