Jouhayna Harmouch
email: jouhayna.harmouch@inria.fr

Bernard Mourrain
email: bernard.mourrain@inria.fr

Houssam Khalil
email: houssam.khalil@ul.edu.lb

Decomposition of Low Rank Multi-Symmetric Tensor

est destinée au dépôt

Introduction

The decomposition of symmetric and multi-symmetric tensors has many applications in engineering disciplines such that signal processing [START_REF] De Lathauwer | Tensor-based techniques for the blind separation of ds-cdma signals[END_REF], scientific data analysis [START_REF] Jiang | Kruskal's permutation lemma and the identification of candecomp/parafac and bilinear models with constant modulus constraints[END_REF], [START_REF] Smilde | Multi-way analysis: applications in the chemical sciences[END_REF], statistics [START_REF] Sidiropoulos | Blind parafac receivers for ds-cdma systems[END_REF], in bioinformatics and spectroscopy [START_REF] Comon | Handbook of Blind Source Separation: Independent component analysis and applications[END_REF], in neuroscience, in phylogenetic .... For instance, the study of symmetric tensor decomposition gives an idea about the geometric structure of intersecting fibers in human brain using the Fibers Orientation Fibers Function described in [START_REF] Megherbi | Détection des croisements de fibre en IRM de diffusion par décomposition de tenseur : Approche analytique[END_REF], ou [START_REF] Connelly | Robust determination of the fiber orientation distribution in diffusion mri: Non-negativity constrained superresolved spherical deconvolution[END_REF] and [START_REF] Atkins | Symmetric positive semi-definite cartesian tensor fiber orientation distributions (ct-fod)[END_REF]. The decomposition of multi symmetric tensors of small rank appear in several other contexts, for learning latent variable models which are algebraic statistics models. This is the case for the analysis of phylogenetic trees model described in [START_REF] Roch | A short proof that phylogenetic tree reconstruction by maximum likelihood is hard[END_REF] or for the analysis of contents of web pages model described in [START_REF] Anandkumar | Tensor Decompositions for Learning Latent Variable Models (A Survey for ALT)[END_REF]. Here, the mixture model is a collection of all non-negative probability tensors of low rank.

The tensor decomposition problem is also very interesting from an algebraic geometric point of view [START_REF] Landsberg | Tensors: Geometry and Applications[END_REF]. Important efforts have been developed over the last decades to better understand the theoretical aspects, as well as the algorithmic aspects of this difficult problem. Some of the well-known decomposition methods use local optimization techniques such Alternate Least Square, Gradient Descents, Quasi-Newton, . . . to minimize the error between the tensor and its decomposition. Some other approaches exploit the algebraic structure associated to the tensor decomposition [START_REF] Brachat | Symmetric tensor decomposition[END_REF], [START_REF] Bernardi | General tensor decomposition, moment matrices and applications[END_REF]. Homotopy techniques have also been used recently to compute such decomposition [START_REF] Bernardi | Tensor decomposition and homotopy continuation[END_REF].

In this paper, we describe a direct method for the decomposition of multi-symmetric tensors, based on simple linear algebra tools. The decomposition algorithm applies to tensors of low enough rank. We follow the approach in [START_REF] Bernardi | General tensor decomposition, moment matrices and applications[END_REF] but directly apply numerically stable linear algebra tools on submatrices of the Hankel matrices to recover the decomposition. In particular, we show how to recover directly the points and weights from eigenvectors of multiplication operators of the quotient algebra associated to the decomposition. The algorithm does not require the solution of polynomial equations. The proposed method extends the techniques of [START_REF] Sanchez | Tensorial resolution: A direct trilinear decomposition[END_REF] to more general tensors and to tensors of higher rank. It is closely connected to the multivariate Prony method investigated in [START_REF] Mourrain | Polynomial-exponential decomposition from moments[END_REF] and to the structured low rank decomposition of Hankel matrix [START_REF] Harmouch | Structured low rank decomposition of multivariate hankel matrices[END_REF].

A multi-linear tensor is in correspondence with a multilinear map from a product of vector spaces to the coefficient field. A tensor symmetric tensor is a tensor whose components stay invariant by any permutation of indices. In the following, we study the general class of multi symmetric tensor decomposition problem, which contains these two classes. We show the correlation between the dual of a tensor, formal power series and then the Hankel matrices associated to them. We use the singular value decomposition of Hankel matrices to compute the decomposition of a tensor of low rank. We exploit the properties of Artinian Gorenstein Algebra to find out some multiplication matrices which help to know the eigen-structure of points associated to linear forms and their weights. We slice variables into bunches of sub-variables and we adapt the description of Artinian Gorenstein Algebra to this case. We adapt the method of decomposition of Hankel matrices of low rank described in [START_REF] Harmouch | Structured low rank decomposition of multivariate hankel matrices[END_REF] to a decomposition of multi linear tensors method which is based on the decomposition of a formal power series as a weighted sum of exponential described in [START_REF] Mourrain | Polynomial-exponential decomposition from moments[END_REF]. The computation of multiplication matrices depend on the dimension of tensor, and the number of given moments or coefficients. We describe the algorithm in 3 dimensional space and we give its numerical implementation using MAPLE. This description gives an idea about the constraints and difficulties of the problem in n dimensional space. We show a numerical example of the decomposition of a tensor of rank 3 with order one in each bunch of 3 variables in 3 dimensional space.

Contributions. We study the decomposition of multi-linear tensor T as a sum of product of powers of linear forms in correlation with the decomposition of its dual T * as a weighted sum of evaluations. T * is defined via a formal power series τ . We exploit the structure of the quotient algebra A τ of the ring of multivariate polynomials in bunches of sub-variables by the kernel of the Hankel operator H τ associated to τ . We choose two bases A 1 and A 2 of monomials such that all given moments of the tensor T appear in the matrix H τ associated to T in the bases A 1 and A 2 and we substitute x 0 by one. We compute the Singular Value Decomposition of the Hankel matrix associated to a chosen truncated bases of A 1 and A 2 such that the multiplication of the matrix by one fixed variable is well defined. We exploit the eigen-structure properties of multiplication operators to compute the sub-coordinates of points and their corresponding weights. We show the constraints which arise from the computation of all multiplication matrices in higher dimension spaces. We propose a new algorithm to compute the sub-coordinates of points using the eigenvalues of multiplication matrices and their transpose. We deduce weights from eigenvectors of a linear combination of multiplication matrices. This method is an adaptation of Structured Low Rank Decomposition of Multivariate Hankel Matrices method proposed in [START_REF] Harmouch | Structured low rank decomposition of multivariate hankel matrices[END_REF] which is the generalization of Prony method. We give a numerical interpretation of the decomposition of a multi-linear tensor of low rank in 3 dimensional space.

Structure of the paper. In the following section, we recall the definition of multi symmetric tensors of rank r and the affine decomposition theory of them. In section 2, we recall some important properties of Artinian Gorenstein Algebra that we adapt to solve the dual decomposition problem which is resumed by the computation of points and their weights. In section 3, we give a theoretical approach of the multi linear symmetric decomposition problem. In section 4, we propose a new algorithm to solve the decomposition problem in 3 dimensional space. Finally, we give an implementation of our algorithm for one example using MAPLE and we interpret the results.

Partial Symmetric Tensor Decomposition Problem

In this section we give the definition of a multi-symmetric tensor as a multi-homogeneous polynomial of a different positive degree at each collection of variables. This polynomial can be defined as well as multi symmetric array of coefficients. In the opposite, for a multi symmetric array of coefficients we can define a multi-homogeneous polynomial and then deshomogenize it. We recall the definition of minimal affine decomposition of a multi symmetric tensor as weighted sum of product of power of linear forms. We show the relationship between the dual of deshomogenized tensor and the formal power series associated to it using the apolar product. Then, after scaling by the linear form of the decomposition and multiplying the weights by the scaling factor we deduce by linearity that the dual of the Tensor can be decomposed as a weighted sum of evaluations.

Definition 1. Let (E j ) 1≤j≤k be a family of n j + 1 dimensional vector spaces, each one of them is of basis x j such that E j = ⟨x j ⟩ = ⟨x j , . . . , x j,nj ⟩. Definition 2. S δj (E j ) is the vector space of homogeneous polynomials in the variables x j of degree δ j . Definition 3. S δ1 (E 1 ) ⊗ S δ2 (E 2 ) ⊗ . . . ⊗ S δ k (E k ) is the vector space of multi-homogeneous polynomials of degree δ j in each subset of variables x j for j = 1, . . . , k, an element [T ] of this vector space is called a multi symmetric tensor. It is denoted hereafter as S δ (E). Definition 4. A multi symmetric tensor of S δ1 (E 1 ) ⊗ S δ2 (E 2 ) ⊗ . . . ⊗ S δ k (E k ) can be interpreted as a multi symmetric array of coefficients

[T ] = [t α ′ 1 ,α ′ 2 ,...,α ′ k ] α ′ j =δj α ′ j ∈N n j +1 such that each α ′ j = (α ′ j,pj ) 0≤pj ≤nj is a multi-index for 1 ≤ j ≤ k.
For α ∈ N n with α ≤ δ, we denote ᾱ = (δ -α , α 1 , . . . , α n ). The multi symmetric tensor is defined as

[T ] = [t ᾱ1, ᾱ1,..., ᾱk ] αj ≤δj αj ∈N n j .
Such tensor is identified with the multi-homogeneous polynomial

T (x 1 , x 2 , . . . , x k ) = ᾱj =δj ᾱj ∈N n j +1 t ᾱ1, ᾱ2,..., ᾱk (x 1 ) ᾱ1 (x 2 ) ᾱ2 . . . (x k ) ᾱk
If we let x j = 1 for j = 1, . . . , k we get

T (x 1 , x 2 , . . . , x k ) = αj ≤δj αj ∈N n j t α1,α2,...,α k (x 1 ) α1 (x 2 ) α2 . . . (x k ) α k
where x j = (x j,1 , . . . , x j,nj ) for j = 1, . . . , k because of (x j ) αj = (x j ) ᾱj for j = 1, . . . , k.

A multilinear tensor is defined when ᾱj = δ j = 1 for j = 1, . . . , k, then by abuse of notation we obtain ᾱj [i j ] = 1 for some 0 ≤ i j ≤ n j and 0 elsewhere, so that the multi symmetric array associated to that tensor is defined as

[T ] = [t i1,i2,...,i k ] 0≤ij ≤nj 1≤j≤k
Given e j basis of E j for j = 1, . . . , k, the tensor [T ] in the basis e 1 ⊗ e 2 ⊗ . . . ⊗ e k is equal to T = ∑0≤i 1≤n1 0≤i2≤n2 ⋮ 0≤i k ≤n k t i1,i2,...,i k e 1,i1 ⊗ e 2,i2 ⊗ . . . ⊗ e k,i k , such a tensor can be identified with the multi-homogeneous poly-

nomial T (x 1 , x 2 , . . . , x k ) = ∑0≤i 1≤n1 0≤i2≤n2 ⋮ 0≤i k ≤n k t i1,i2,...,i k x 1,i1 x 2,i2 . . . x k,i k because of (x j ) ᾱj = x j,ij for some 0 ≤ i j ≤ n j
and for all 1 ≤ j ≤ k.

The dual of the tensor is

T * (y 1 , y 2 , . . . , y k ) = ∑0≤i 1≤n1 0≤i2≤n2 ⋮ 0≤i k ≤n k t i1,i2,...,i k y 1,i1 y 2,i2 . . . y k,i k because of (y j ) ᾱj = y j,ij
for some 0 ≤ i j ≤ n j and for all 1 ≤ j ≤ k.

We denote R δ1,δ2,...,δ k the space obtained by the deshomogeneisation of elements in S δ (E) by setting

x j = 1 for j = 1, . . . , k where R = C[x 1 , x 2 , . . . , x k ]
is the space of polynomials in the variables x j = (x j,1 , . . . , x j,nj ) for j = 1, . . . , k Definition 5. The tensor decomposition problem of T (x 1 , x 2 , . . . , x k ) is the decomposition of T as a sum of product of power of linear forms such that

T (x 1 , x 2 , . . . , x k ) = ∑ r p=1 ω p u δ1 p,1 (x 1 )u δ2 p,2 (x 2 ) . . . u δ k i,k (x k )
where u p,j (x j ) = u p,j x j + u p,j,1 x j,1 + . . . + u p,j,nj x j,nj and

u p = (u p,j,pj ) 0≤pj ≤nj 1≤j≤k = (u p,1 , u p,1,1 , . . . , u p,1,n1 , u p,2 , u p,2,1 , . . . , u p,2,n2 , . . . . . . , u p,k , u p,k,1 , . . . , u p,k,n k ) ∈ C ∑ k j=1 (nj +1)
is the coefficient vector associated to the linear forms u p,j (x j ) in the basis x j for j = 1, . . . , k. Definition 6. The minimal number of terms in a decomposition of T (x) is called the rank of T .

We say that T (x 1 , x 2 , . . . , x k ) has an affine minimal decomposition of the previous form if u p,j ≠ 0 for p = 1, . . . , r and j = 1, . . . , k where r is the rank of T .

Definition 7. For T = (t α1,α2,...,α k ) αj ≤δj αj ∈N n j ∈ S δ (E) we denote τ α1,α2,...,α k (T ) = τ α1,α2,...,α k = t α1,α2,...,α k δ1 α1 -1 δ2 α2 -1 . . . δ k α k -1
. The dual of the tensor T (x 1 , x 2 , . . . , x k ) ∈ S δ (E) is defined via the formal power series as τ (y 1 , y 2 , . . . ,

y k ) = T * (y 1 , y 2 , . . . , y k ) = ∑ αj ≤δj αj ∈N n j τ α1,α2,...,α k (y1) ᾱ1 ᾱ1! (y2) ᾱ2 ᾱ2! . . . (y k ) ᾱk ᾱk !
where (y j ) ᾱj = (y j , y j,1 , . . . , y j,nj ) (αj ,αj,1,...,αj,n j ) = ∏ nj pj =0 (y j,pj ) αj,p j for j = 1, . . . , k Definition 8. For a polynomial p ∈ R and a formal power series τ ∈ R * , we define the multiplication operator * such that

p * τ ∶ R → C q ↦ τ (p.q) Definition 9. Let T 1 (x 1 , x 2 , . . . , x k ) and T 2 (x 1 , x 2 , . . . , x k ) be two tensors of S δ (E). The apolar product of T 1 (x 1 , x 2 , . . . , x k ) and T 2 (x 1 , x 2 , . . . , x k ) is defined as ⟨T 1 (x 1 , x 2 , . . . , x k ), T 2 (x 1 , x 2 , . . . , x k )⟩ = ∑ αj ≤δj αj ∈N n j τ (1) α1,α2,...,α k τ (2) α1,α2,...,α k δ α where δ α = δ1 α1 δ2 α2 . . . δ k α k .
Definition 10. The dual operator of a tensor is defined as

T * ∶ (R δ1,δ2,...,δ k ) → (R δ1,δ2,...,δ k ) * (1)
T 2 ↦ T * (T 2 ) = ⟨T (x), T 2 (x)⟩ (2)
Lemma 1. By a generic change of coordinates in each E j , we may assume that u p,j ≠ 0 and that T has an affine decomposition. Then by scaling u p (x) and multiplying ω p by the d th power of the scaling factor we may assume that u p,j = 1 for p = 1, . . . , r and j = 1, . . . , k. Thus the polynomial

T (x) = ∑ r p=1 ω ′ p u ′ p δ (x) = ∑ r p=1 ω ′ i u ′ p,1 δ1 (x 1 )u ′ p,2 δ2 (x 2 ) . . . u ′ p,k δ k (x k )
Proposition 1. The dual of the product of powers of linear forms u δ1 1 u δ2 2 . . . u δ k k is the evaluation e u at u = (u 1 , u 2 , . . . , u k ).

Proof. For T = u δ1 1 u δ2 2 . . . u δ k k and any T ′ ∈ R δ1,δ2,...,δ k , we check that ⟨T (x), T ′ (x)⟩ = T ′ (u). This shows that T * coincides with the evaluation e u .

Thus if

T = ∑ i ω i u δ1 i,1 u δ2 i,2 . . . u δ k i,k
, then T * coincides with the weighted sum of evalautions T * = ∑ i ω i e bui on R δ1,δ2,...,δ k . We reduce the decomposition problem of T to the decomposition of T * as a weighted sum of evaluations T * = ∑ i ω i e bui .

Solving polynomial equations by eigenvector computation

We recall the definition of quotient algebra A of the ring of polynomials in a collection of variables x = [x 1 , x 2 , . . . , x k ] where x j = [x j,1 , x j,2 , . . . , x j,nj ] for j = 1, . . . , k by an ideal I of multivariate polynomials. We use eigen-structure properties of multiplication operators and their transpose in a chosen basis of A and its dual to compute the x's coordinates of points in the decomposition of the associated multi-linear tensor.

A quotient algebra A = C[x] I is Artinian if it is of finite dimension over C. In the case of partial symmetric tensor, the variables x and y are divided into bunches of sub-variables such that x = [x 1 , x 2 , . . . , x k ] and y = [y 1 , y 2 , . . . , y k ]. Hereafter, in the case of multi-symmetric tensor the ideal I defines a finite number of roots V(I) = {ξ 1 , ξ 2 , . . . , ξ r ′ } = {ξ ∈ C n ∀q ∈ I, q(ξ) = 0} where n = ∑ k j=1 (n j + 1) such that n j + 1 is the dimension of each vector space spanned by x j and we have a decomposition of A as a sum of sub-algebras:

A = C[x] I = A 1 ⊕ ⋯ ⊕ A r ′
where A p = u ξp A ∼ C[x] Q p and Q p is the primary component of I associated to the root ξ p ∈ C n . The elements u 1 , . . . , u r ′ satisfy the relations

u 2 ξp (x) ≡ u ξp (x), r i=1 u ξp (x) ≡ 1.
The polynomials u ξ1 , . . . , u ξ r ′ are called idempotents of A. The dimension of A p is the multiplicity of the point ξ p . For more details, see [START_REF] Elkadi | Introduction à la résolution des systèmes polynomiaux[END_REF][Chap. 4].

For g ∈ C[x], the multiplication operator M g is defined by

M g ∶ A → A h ↦ M g (h) = g h.
The transpose M T g of the multiplication operator M g is

M T g ∶ A * → A * Λ ↦ M T g (Λ) = Λ ○ M g = g ⋆ Λ.
The main property that we will use to recover the roots is the following [START_REF] Elkadi | Introduction à la résolution des systèmes polynomiaux[END_REF][Thm. 4.23]:

Proposition 2. Let I be an ideal of C[x] and suppose that V(I) = {ξ 1 , ξ 2 , . . . , ξ r ′ }. Then for all g ∈ A, the eigenvalues of M g and M T g are the values g(ξ 1 ), . . . , g(ξ r ′ ) of the polynomial g at the roots with multiplicities µ p = dim A p .

-The eigenvectors common to all M T g with g ∈ A are -up to a scalar -the evaluations e ξ1 , . . . , e ξ r ′ .

If B = {b 1 , . . . , b r } is a basis of A, then the coefficient vector of the evaluation e ξp in the dual basis of B is ⟨e ξp b j ⟩ β∈B = [b j (ξ p )] p=1...r = B(ξ p ). The previous proposition says that if M g is the matrix of M g in the basis B of A, then M T g B(ξ p ) = g(ξ p ) B(ξ p ).

If moreover the basis B contains the monomials 1, x 1,1 , x 1,2 , . . . , x 1,n1 , then the common eigenvectors of M T g are of the form v p = c [1, ξ p,1,1 , . . . , ξ p,1,n1 , . . .] and the x's coordinates of the root ξ p can be computed from the coefficients of v p by taking the ratio of the coefficients of the monomials x 1,1 , . . . , x 1,n1 by the coefficient of 1:

ξ p,1,i1 = vp,1,i 1 +1
vp,1,1 . Thus computing the common eigenvectors of all the matrices M T g for g ∈ A yield the x's coordinates of the roots ξ p (p = 1, . . . , r).

In practice, it is enough to compute the common eigenvectors of

M T x1,1 , . . . , M T x1,n 1 , since ∀g ∈ C[x 1 ], M T g = g(M T
x1,1 , . . . , M T x1,n 1 ). Therefore, the common eigenvectors M T x1,1 , . . . , M T x1,n 1 are also eigenvectors of any M T g . The multiplicity structure, that is the dual Q ⊥ p of each primary component Q p of I, also called the inverse system of the point ξ p can be deduced by linear algebra tools (see e.g. [START_REF] Mourrain | Isolated points, duality and residues[END_REF]).

In the case of simple roots, we have the following property [START_REF] Elkadi | Introduction à la résolution des systèmes polynomiaux[END_REF][Chap. 4]:

Proposition 3. If the roots {ξ 1 , ξ 2 , . . . , ξ r } of I are simple (i.e. µ p = dim A p = 1) then we have the following:

u = {u ξ1 , . . . , u ξr } is a basis of A.

-The polynomials u ξ1 , . . . , u ξr are interpolation polynomials at the roots ξ p : u ξp (ξ q ) = 1 if p = q and 0 otherwise.

-The matrix of M g in the basis u is the diagonal matrix diag(g(ξ 1 ), . . . , g(ξ r )).

This proposition tells us that if g is separating the roots, i.e. g(ξ p ) ≠ g(ξ q ) for p ≠ q, then the eigenvectors of M g are, up to a scalar, interpolation polynomials at the roots.

Artinian Gorenstein algebra of a multivariate Hankel operator

In this section, we detail the construction of the quotient algebra A τ by the kernel I τ of the Hankel operator H τ associated to the dual of the tensor T . We compute a basis of A τ such that the submatrix associated to it has a maximal non-zero minor of a truncated matrix of H τ . We recall how to compute the multiplication matrices in this associated basis and its dual using some shifted submatrices of H τ . We notice that not all of them are easy to compute. We benefit from properties of generalized eigenvalues of multiplication matrices by y ′ s to compute the x's coordinates of points. We show how to use the generalized eigenvectors of the multiplication matrices to compute the weights. We associate to a Hankel operator H τ , the quotient

A τ = C[x] I τ of the polynomial ring C[x] modulo the kernel I τ = {p ∈ C[x] ∀q ∈ R, ⟨τ pq⟩ = 0} of H τ .
We check that I τ is an ideal of C[x], so that A τ is an algebra.

As

A τ = C[x] I τ ∼ img H τ , the operator H τ is of finite rank r, if and only if, A τ is Artinian of dimension dim C A τ = r . A quotient algebra A is called Gorenstein if its dual A * = Hom C (A, C
) is a free A-module generated by one element.

In our context, we have the following equivalent properties [START_REF] Mourrain | Polynomial-exponential decomposition from moments[END_REF]:

-τ = ∑ r ′ p=1 ω p (y)e ξp (y) with ω p ∈ C[y], ξ p ∈ C n and ∑ r ′ p=1 µ(ω p ) = r where n = ∑ k j=1 (n j + 1) , -H τ is of rank r, -A τ is an Artinian Gorenstein algebra of dimension r.
Another property that will be helpful to determine a basis of A τ is the following:

Lemma 2. Let B = {b 1 , . . . , b r }, B ′ = {b ′ 1 , . . . , b ′ r } ⊂ C[x]. If the matrix H B,B ′ τ = (⟨τ b p b ′ q ⟩) 1≤p,q≤r is invertible, then B and B ′ are linearly independent in A τ .
By this Lemma, bases of A τ can be computed by identifying non-zero minors of maximal size of the matrix of H τ . Proposition 4. Let B, B ′ be basis of A τ and g ∈ C[x]. We have

H B,B ′ g⋆τ = (M B g ) T H B,B ′ τ = H B,B ′ τ M B ′ g . (3)
where M B g (resp. M B ′ g ) is the matrix of the multiplication by g in the basis B (resp. B ′ ) of A τ . We deduce the following property: Proposition 5. Let τ (y) = ∑ r p=1 ω p (y)e ξp (y) with ω p ∈ C[y] ∖ {0} and ξ p ∈ C n distinct and let B, B ′ be bases of A τ . We have the following properties:

-For g ∈ C[x], M B ′ g = (H B,B ′ τ ) -1 H B,B ′ g⋆τ , (M B g ) T = H B,B ′ g⋆τ (H B,B ′ τ ) -1 . -For g ∈ C[x], the generalized eigenvalues of (H B,B ′ g⋆τ , H B,B ′ τ
) are g(ξ p ) with multiplicity µ p = µ(ω p ), p = 1, . . . , r.

-The generalized eigenvectors common to all

(H B,B ′ g⋆τ , H B,B ′ τ ) for g ∈ C[x] are -up to a scalar -(H B,B ′ τ ) -1 B(ξ p ), p = 1, . . . , r.
Proof. The two first points are direct consequences of Propositions 4 and 2. The third point is also a consequence of Proposition 2, since the coordinate vector of the evaluation e ξp in the dual basis of B is B(ξ p ) for p = 1, . . . , r.

This proposition shows that the matrices of multiplication by an element g in A, and thus the roots {ξ 1 , . . . , ξ r } and their multiplicity structure, can be computed from truncated Hankel matrices, provided we can determine bases B, B ′ of A τ . In practice, it is enough to compute the generalized eigenvectors common to

(H B,B ′ x1,i 1 ⋆τ , H B,B ′ τ ) for i 1 = 1, . . . , n 1 to recover the roots. As H B,B ′ x1,i 1 ⋆τ = H x1,i 1 B,B ′ τ = H B,x1,i 1 B ′ τ
, the decomposition can be computed from sub-matrices of

H B,B ′+ τ or H B + ,B ′ τ where B + = B∪x 1,1 B∪⋯∪x 1,n1 B, B ′+ = B ′ ∪ x 1,1 B ′ ∪ ⋯ ∪ x 1,n1 B ′ .

Multilinear Tensor Decomposition Problem

In this section, we analyze the easiest case of multi symmetric tensor where it is of degree one at each bunch of sub-variables. Our goal is to decompose τ which is equal to T * as a weigthed sum of evaluations by computing the eigen-structure of A τ which is based on the computation of multiplication operators. We simplify notations by using subscripts of variables and coefficients instead of multi-index exponents. We compute the truncated Singular Value Decomposition of a generic linear combination of a shifted Hankel matrices by the first collection of variables. By linearity and properties of the multiplication operators by one variable described in section 2, we deduce the multiplication operators by more complex variables which could be used to compute weights and points.

We choose two monomial bases B 1 and B 2 indexing respectively rows and columns of the Hankel matrix H B1,B2 T * associated to the tensor T * , such that the set of monomials {B 1 * B 2 x j,ij , 0 ≤ i j ≤ n j , 1 ≤ j ≤ k} span the set of deshomogenized polynomials R δ1,δ2,...,δ k .

The matrix of the truncated Hankel operator in the basis B 1 and the dual basis

B 2 is H B1,B2 T * = [t i1,i2,...,i k ] 0≤i1≤n1 0≤i2≤n2 ⋮ 0≤i k ≤n k .
The Hankel matrix associated to the tensor x 1,i1 * T * is defined as

H 1,i1 = H B1,B2 x1,i 1 * T * = H x1,i 1 * B1,B2 T * = [t α+β ] α∈x1,i 1 * B1,β∈B2
, all the elements of the matrix are divisible in x 1,i1 and of degree δ.

For example, the Hankel matrix associated to x 1 * T * in the monomials basis B 1 and B 2 is denoted by H 0 . Let λ(x 1 ) = λ 0 + λ 1 x 1,1 + . . . + λ n1 x 1,n1 is a linear form with generic chosen coefficients λ i1 , i 1 = 0, . . . , n 1 , we build a linear combination of H 1,i1 , i 1 = 0, . . . , n 1 such that Ĥ0 = ∑ n1 i1=0 λ i1 H 1,i1 we compute the singular value decomposition of it.

Computing the singular value decomposition of Ĥ0 , we obtain

Ĥ0 = U SV T
where S is the diagonal matrix of all singular values of Ĥ0 arranged in a decreasing order, U is an unitary matrix whose columns are the left singular vectors of Ĥ0 , V is an unitary matrix whose columns are the right singular vectors of Ĥ0 . We denote by U H the hermitian transpose of U and V the conjugate of V . We denote by U r and V r the truncated matrices of the first r columns of U and V and S r the diagonal matrix of the first r rows and r columns of S.

We denote B 1 = ⟨1, x 1,1 , . . . , x 1,n1 ⟩ and B 2 = ⟨1, x k,1 , . . . , x k,n k ⟩. Let u i = [u α,i ] α∈B1 and v j = [v β,j ] β∈B2 be respectively the i th and j th columns of U H and V . We denote by u i (x 1 ) = u T i U H r and v j (x 1 ) = v T j V r the corresponding polynomials. The bases formed by these first r polynomials are denoted U H r ∶= (u i (x 1 )) i=1,...,r and V r ∶= (v j (x 1 )) j=1,...,r . We will also denote by U H r (resp. V r ) the corresponding coefficient matrix, formed by the first rows (resp. columns) of U H (resp. V ). We denote by S r the diagonal matrix of the first r rows and columns of S, formed by the first r singular values.

We denote by H r 0 ,H r 1,i1 and Ĥr 0 the matrices obtained by the truncated singular value decomposition of H 0 ,H i1 and Ĥ0 respectively.

We have the following property

H r i1 = (M U H r x1,i 1 ) T H r 0 = H r 0 M V r x1,i 1 * T where M U H r x1,i 1 (resp. M V r x1,i 1 )
is the multiplication matrix by x 1,i1 in the basis U H r (resp. V r ) and M V r x1,i 1 * T is the multiplication matrix by x 1,i1 * T in the basis V r . Then by linearity, we obtain Ĥr

0 = ∑ n1 i1=0 λ i1 H r 1,i1 = H r 0 ∑ n1 i1=0 λ i1 M V r x1,i 1 * T = H r 0 M V r λ(x1) * T . Then ( Ĥr 0 ) -1 = (M V r λ(x1) * T ) -1 (H r 0 ) -1
so multiplying by the first equation we get

( Ĥr 0 ) -1 H r 1,i1 = (M V r λ(x1) * T ) -1 M V r x1,i 1 * T = M V r (x1,i 1 λ(x1)) * T
We compute the eigenvalues and the eigenvectors of the multiplication matrices M V r (x1,i 1 λ(x1)) * T in order to obtain the weights and the points of the decomposition.

Algorithm

We describe now the algorithm to recover the sum T * (x, y, z) = ∑ r p=1 ω p e up (x, y, z),

ω p ∈ C ∖ {0}, u p ∈ C ∑ 3 l=1 (n l +1)
, from the moments of degree at most one at each bunch of coordinates (t i,j,k ) 0≤i≤n1 0≤j≤n2 0≤k≤n3

of the formal power series. To simplify, we change notations to better understand the nine dimensional multivariate space seen as three dimensional space. We only use 3 bunches of variables such that x 1,i1 by x i and x 2,i2 by y j and x 3,i3 by z k .

Algorithm 3.1: Decomposition of polynomial-exponential series with constant weights Input: the moments (t i,j,k ) 0≤i≤n1 0≤j≤n2 0≤k≤n3

of T * .

Compute the monomial sets

A 1 = (x i y j ) 0≤i≤n1 0≤j≤n2
and A 2 = (z 0 , z 1 , . . . , z n3 ) and substitute the x 0 , y 0 and z 0 by 1 to define B 1 and B 2 .

Compute the Hankel matrix

H B1,B2 T * = [t i,j,k ] 0≤i≤n1 0≤j≤n2 0≤k≤n3
for the monomial sets B 1 and B 2 .

Compute the singular value decomposition of H

B1,B2 T * = U SV T where B 1 = ⟨1, x 1 , . . . , x n1 ⟩ and B 2 = ⟨1, z 1 , . . . , z n3 ⟩ with singular values s 1 ≥ s 2 ≥ ⋯ ≥ s m ≥ 0. 4.
Determine its numerical rank, that is, the largest integer r such that sr s1 ≥ . 5. Form the multiplication matrices by y j in the basis V r , M V r yj = S -1 r U H r H B1,B2 yj * T * V r where H B1,B2 yj * T * is the Hankel matrix associated to y j ⋆ T * for j = 1, . . . , n 2 . 6. Compute the eigenvectors v p of ∑ n2 j=1 l j M V r yj such that l j ≤ 1, j = 1, . . . , n 2 and for each p = 1, . . . , r do the following:

-The y ′ s coordinates of the u p are the eigenvalues of the multiplication matrices by y j . Use the formula M V r yj v p = u p,2,j v p for p = 1, . . . , r and j = 1, . . . , n 2 and deduce the u p,2,j . -Write the matrix H B1,B2 T * in the basis of interpolation polynomials(ie. the eigenvectors v p ) and use the corresponding matrix T to compute the z ′ s coordinates. Divide the k th row on the first row of the matrix T to obtain the values of u p,3,k for p = 1, . . . , r and k = 1, . . . , n 3 .

-The x ′ s coordinates of u p are computed using the eigenvectors of the transpose of the matrix M V r yj . They -are up to scalar-the evaluations, they are represented by vectors of the form The cost of the SVD computation is in O(s 3 ) where s ≥ r is the maximal size of the Hankel matrix H 0 and r the rank of the decomposition. The computation of each multiplication matrice is in O(r 2 ) and the eigencomputation is in O(r 3 ). This yields a complexity bound in O(s 3 + n r 2 ) for the complete algorithm, where n = max(n 1 , n 2 , n 3 ) is a bound on the dimension of the spaces. This complexity bound extends to the decomposition of general multi-symmetric tensors, provided r = rank H 0 .

v * p = µ p [1, u p,

Example

In this section, we illustrate the decomposition algorithm on a multi-linear tensor of degree one at each bunch of 3 variables and of rank 3.

If δ l = 1 for all l = 1, . . . , k, k > 1 and n l = n, let k = 3, n l = n = 2, r = 3 and δ l = 1 then we have x = (x 0 , x 1 , x 2 ), y = (y 0 , y 1 , y 2 ) and z = (z 0 , z 1 , z 2 ). For ᾱ ∈ N 3 , ᾱ = 1 ⇒ ᾱ = (1), (0, 1), (0, 1) ⇒ x ᾱ = x i , i = 0, . . . , 2 β ∈ N 3 , β = 1 ⇒ β = (1), (0, 1), (0, 1) ⇒ y β = y j , j = 0, . . . , 2 γ ∈ N 3 , γ = 1 ⇒ γ = (1), (0, 1), (0, 1) ⇒ z γ = z k , k = 0, . . . , 2 The multi symmetric tensor is defined by a multi symmetric array of coefficients such that t α,β,γ ∶= t ᾱ, β,γ = t i,j,k then T (x, y, z) = ∑0≤i≤2

0≤j≤2 0≤k≤2 t i,j,k x i y j z k = 0.4461757334x 0 y 0 z 0 -0.2262004083x 0 y 0 z 1 +0.4427031740x 0 y 0 z 2 - 0.2756785277x 0 y 1 z 0 +0.1612318550x 0 y 1 z 1 -0.3100164212x 0 y 1 z 2 -0.1209490221x 0 y 2 z 0 +0.1465160338x 0 y 2 z 1 - 0.1169341103x 0 y 2 z 2 -0.01239930649x 1 y 0 z 0 -0.05189330981x 1 y 0 z 1 +0.01803564422x 1 y 0 z 2 -0.01336683543x 1 y 1 z 0 + 0.02632784503x 1 y 1 z 1 -0.02598209626x 1 y 1 z 2 -0.3195263612x 1 y 2 z 0 +0.09311605022x 1 y 2 z 1 -0.1116610246x 1 y 2 z 2 - 0.1460187051x 2 y 0 z 0 +0.06557223848x 2 y 0 z 1 -0.1734312692x 2 y 0 z 2 +0.1010145926x 2 y 1 z 0 -0.05743078561x 2 y 1 z 1 + 0.1238292801x 2 y 1 z 2 -0.1485221955x 2 y 2 z 0 + 0.03231415762x 2 y 2 z 1 -0.376925099e -2x 2 y 2 z 2 . Let x 0 = y 0 = z 0 = 1 then T (x, y, z) = ∑1≤i≤2 1≤j≤2 1≤k≤2 t i,j,k x i y j z k
Then the tensor decomposition problem is T (x, y, z) = ∑ r p=1 ω p u p,1 (x)u p,2 (y)u p,3 (z). Given all the moments of degree at most one at each bunch of coordinates

(t i,j,k ) 0≤i≤2 0≤j≤2 0≤k≤2 , We create two sets A 1 = (x ᾱy β ) ᾱ =1 β =1 = (x i y j ) 0≤i≤2 0≤j≤2 and A 2 = (z γ ) γ =1 = (z k ) 0≤k≤2 so that A 1 = (x 0 y 0 , x 0 y 1 , x 0 y 2 , x 1 y 0 , x 1 y 1 , x 1 y 2 , x 2 y 0 , x 2 y 1 , x 2 y 2 ) and A 2 = (z 0 , z 1 , z 2 ). For x 0 = y 0 = z 0 = 1 then B 1 = (1, y 1 , y 2 , x 1 , x 1 y 1 , x 2 , x 2 y 1 , x 2 y 2 ) and B 2 = (1, z 1 , z 2 )
, the Hankel matrix associated to the tensor in the monomial basis B 1 and B 2 is . We denote by U H the Hermitian transpose of U and V the conjugate of V .

H B1,B2 T * = [t ᾱ+ β+γ ] ᾱ =1 β =1 γ =1 = 1 y1 y2 x1 x1y1 x1y2 x2 x2y1 x2y2 ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ t 0 t 1 t 2 t 0,1 t 0,2 t 1,1 t 1,2 t 2,1 t 2,2 1 t 0,1 t 1,1 t 2,1 t 0,1,1 t 0,2,1 t 1,1,1 t 1,2,1 t 2,1,1 t 2,2,1 z1 t 0,2 t 1,2 t 2,2 t 1,2 t 0,2,2 t 1,1,2 t 1,2,2 t 2,1,2 t 2,2,2 z2 = 1 y1 y2 x1 x1y1 x1y2 x2 x2y1 x2y2 ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ 0
H B1,B2 T * = U SV T = 1 x1 x2 ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ t 0 t 1 t 2 1 t 0,1 t 1,1 t 2,1 z1 t 0,2 t 1,2 t 2,2 z2 = 1 x1 x2 ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ 0
Let v i = [v α,i ] α∈B1 and w j = [w β,j ] β∈B2 be respectively the i th and j th columns of U H and V . We denote by v i (x) = v T i U H r and w j (z) = w T j V r the corresponding polynomials. The bases formed by these first r polynomials are denoted U H r ∶= (v i (x)) i=1,...,r and V r ∶= (w j (z)) j=1,...,r . We will also denote by U H r (resp. V r ) the corresponding coefficient matrix, formed by the first rows (resp. columns) of U H (resp. V ). We denote by S r the diagonal matrix of the first r rows and columns of S, formed by the first r singular values. To compute the multiplication matrices M V r y1 and M V r y2 we need to compute the following matrices Then we compute M V r y1 = S -1 r U H r H B1,B2 y1 * T * V r and M V r y2 = S -1 r U H r H B1,B2 y2 * T * V r , and the eigenvectors v p of ∑ 2 j=1 l j M V r yj such that l j ≤ 1, j = 1, . . . , 2. To recover the points u p ∈ C n * k for p = 1, . . . , r of the form We do the following: In general we have M V r xj,i j v i = u i,j,ij v i , for i = 1, . . . , r, j = 1, . . . , k, i j = 1, , n j so in this case we get M V r y1 v p = u p,2,1 v p and M V r y2 v p = u p,2,2 v p for p = 1, . . . , 3 so we compute u p,2,1 and u p,2,2 for p = 1, . . . , 3. So that we get

H B1,B2 y1 * T * = y1 y1x1 x2y1 ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ t 0,
u p = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ u p,

1 , 1 p [ 1 ]

 111 , . . . , u p,1,n1 ]. Compute v * p as the p th column of the transpose of the inverse of the matrixV = [v 1 , . . . , v r ] for p = 1, . . . , r and deduce u p,1,i = v * p [i+1] v * for p = 1, . . . , r and i = 1, . . . , n 1 . -Compute ω p = ⟨T * vp⟩ vp(up) . Output: r ∈ N, ω p ∈ C (0), u p ∈ C ∑ 3 l=1 (n l +1) , p = 1, . . . , r such that T * (x, y, z) = ∑ r p=1ω p e up (x, y, z) up to degree one at each bunch of coordinates.

u p, 2 , 1 =T * (v 1 )λ 1 ω 1 λ 2 ω 2 λ 3 ω 3 1 λ 1

 21111 -0.746329870878 -0.293761776025 -0.304898408788 u p,2,2 = 1.40328849510 -0.336304368405 -3.59031087599 .The eigenvectors v p ∈ ⟨1, x 1 , x 2 ⟩ for p = 1, . . . , 3 are up to a scalar the interpolation polynomials at the roots so that if the dual of the tensor has an affine decompositionT * (x, y, z) = ∑ r p=1 ω p e up (x, y, z) then T * (v p ) = ∑ r p=1 ω p e up (v p ) = λ p ω p , T * (z 1 v p ) = ∑ r p=1 ω p e up (z 1 v p ) = λ p ω p u p,3,1 and T * (z 2 v p ) = ∑ r p=1ω p e up (z 2 v p ) = λ p ω p u p,3,2 , for p = 1, . . . , 3. Then the values of u p,3,1 and u p,3,2 for p = 1, . . . , 3 come from the computation of the matrix:T * (v 2 ) T * (v 3 ) 1 T * (z 1 v 1 ) T * (z 1 v 2 ) T * (z 1 v 3 ) z1 T * (z 2 v 1 ) T * (z 2 v 2 ) T * (z 2 v 3 ) z2= ω 1 u 1,3,1 λ 2 ω 2 u 2,3,1 λ 3 ω 3 u 3,3,1 z1 λ 1 ω 1 u 1,3,2 λ 2 ω 2 u 2,3,2 λ 3 ω 3 u 3,3,2 z2

  All the entries of this matrix are known, we choose B 1 = ⟨1, x 1 , x 2 ⟩ and B 2 = ⟨1, z 1 , z 2 ⟩ to be able to multiply by y 1 and to compute the multiplication matrix. Computing the singular value decomposition of

	.4461757 0.2756785 0.1209490 0.01239930 0.01336683 0.3195263 0.14601870 0.1010145 0.1485221
	-0.2262004 0.1612318 0.1465160 0.05189330 0.02632784 0.09311605 0.06557223 0.05743078 0.03231415
	0.4427031 0.3100164 0.1169341 .01803564 0.02598209 0.1116610 0.1734312 0.1238292 0.00376925
	H B1,B2 T *	, we obtain

  1 t 1,1 t 2,1 1 t 0,1,1 t 1,1,1 t 2,1,1 z1 t 0,1,2 t 1,1,2 t 2,1,2 z2 t 1,2 t 2,2 1 t 0,2,1 t 1,2,1 t 2,2,1 z1 t 0,2,2 t 1,2,2 t 2,2,2 z2

			y1	y1x1	x2y1
			= ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ -0.2756785 -0.01336683 0.1010145 0.1612318 0.02632784 -0.05743078 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ -0.3100164 -0.02598209 0.1238292
		y2 y2x1 x2y2	y2	y2x1	x2y2
	H B1,B2 y2 * T * = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣	⎤ ⎥ ⎥ ⎥ ⎥ ⎦ t 0,2 = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ -0.1209490 -0.3195263 -0.1485221 0.1465160 0.09311605 0.03231415 -0.1169341 -0.1116610 -0.00376925	⎤ ⎥ ⎥ ⎥ ⎥ ⎦

Therefore the value of u p,3,1 (resp. u p,3,2 ) comes from the ratio of the second row (resp. the third row) and the first row of the matrix for p = 1, . . . , 3. So that we get u p,3,1 = -0.655842579065 0.0321233423462 -0.520955291 u p,3,2 = 1.24749588143 0.403506877499 0.242728128570 .

The common eigenvectors of all (M V r yj ) T -are up to scalar-the evaluations, they are represented by vectors of the form v * p = µ p [1, u p,1,1 , u p,1,2 ] in the dual basis of B 1 = ⟨1, x 1 , x 2 ⟩ then the computation of the coordinates of u p,1,1 and u p,1,2 come from the eigenvectors of the transpose of multiplication operators which are obtained by transposing the inverse of the matrix V of vectors of M V r yj for j = 1, . . . , 2, therefore the value of u p,1,1 (resp. u p,1,2 ) comes from the ratio of the second element of v * p (resp. the third element) and the first element of it, so that u p,1,1 = 0.114279629148 -1.08600705528 1.23814628617 u p,1,2 = -0.405714894278 -0.567603220082 0.873482418287 .

Notice that the computation of ω p , p = 1, . . . , 3 can be done using the following formula


Minh Quan Ho

Christian Obrecht

Bernard Tourancheau

Benoît Dupont De Dinechin

Julien Hascoet

Improving 3D Lattice Boltzmann Method stencil with asynchronous transfers on many-core processors

CPU-based many-core processors present an alternative to multicore CPU and GPU processors. In particular, the 93-Petaflops Sunway supercomputer, built from clustered many-core processors, has opened a new era for high performance computing that does not rely on GPU acceleration. However, memory bandwidth remains the main challenge for these architectures. This motivates our endeavor for optimizing one of the most data-intensive kind of stencil computations, namely the three-dimensional applications of the lattice Boltzmann method (LBM). We propose optimizations on many-cores processors by using local memory and asynchronous software-prefetching on a representative 3D LBM solver as an example. We achieve 33 % performance gain on the Kalray MPPA-256 manycore processor by actively streaming data from/to local memory, compared to the "passive" OpenCL programming model.

Introduction

Since the last decade, the lattice Boltzmann method (LBM) has become widely used in computational fluid dynamics for incompressible and weakly compressible flows, see e.g. [START_REF] Succi | The lattice Boltzmann equation: for fluid dynamics and beyond[END_REF]. An LBM model is characterized by its stencil type, denoted DdQq, where d is the number of space dimensions (one, two or three) and q is the number of particle distribution functions (PDFs), see e.g. [START_REF] Cao | Physical symmetry and lattice symmetry in the lattice boltzmann method[END_REF]. Physically, an LBM time step on a lattice node consists of a collision and a propagation step (also known as streaming step). The collision applies a pre-defined physical model on the lattice distribution vector of q values. The propagation then updates these new distribution values to the node itself and q -1 of its neighboring nodes. The most used stencil types are D2Q5, D2Q9 and D3Q19 (see Fig. 1), or D3Q27.

From a programming point of view, LBM kernels are easy to implement and wellsuited for parallelization on recent multi-/many-core platforms. However, lattice Boltzmann methods are known for their low arithmetic intensity and particularly high memory bandwidth requirement. Taking the example of a basic LBM solver, depending on collision operator, between 200 and 400 floating-point operations are performed on a lattice node per time step. Most D3Q19 LBM implementations require storing all the 19 distribution values for each lattice node. A lattice domain L × L × L contains 19 × L 3 single-or double-precision floating-point numbers. Updating this lattice grid in a single time-step requires 19 × 2 × L 3 load/store memory operations for less than 400 × L 3 arithmetic operations. Thus, simulating the whole lattice domain through T time-steps will generate a huge amount of data movement of 19 × 2 × L 3 × T floating-point numbers for 400 × L 3 × T floating-point operations. While recent architectures gain computing performance by increasing the clock speed and multiplying the number of cores, evolution of memory systems still cannot fetch enough data to keep cores busy. The dataset cannot always fit in caches and must be stored in the main (even remote) memory with much higher latency. The low arithmetic intensity of stencil kernels like LBM is thus the limit of performance, as well as their poor data-locality which reduces significantly the cachereuse ratio. Previous studies by [START_REF] Mcintosh-Smith | On the performance portability of structured grid codes on many-core computer architectures[END_REF] and [START_REF] Obrecht | Performance evaluation of an opencl implementation of the lattice boltzmann method on the intel xeon phi[END_REF] show that LBM implementations are memorybound and hardly obtain good performance on CPU or Xeon Phi processors. GPU-based accelerators, thanks to their graphics-dedicated high-bandwidth memory, appear to be the most suitable platforms for LBM today. However, their low capacity of local memory inhibits optimization techniques for data prefetching (to reduce transfer time) and data sharing between cores (for stencil neighboring dependencies).

Although other clustered many-core processors have much less global memory bandwidth than GPUs do, they embed significant amount of fast local memory, see [START_REF] Dupont De Dinechin | A clustered manycore processor architecture for embedded and accelerated applications[END_REF] and [START_REF] Fu | The sunway taihulight supercomputer: system and applications[END_REF], and provide more predictability in both computing time and data transfer. This enables using explicit and efficient user buffers for elaborate optimizations, such as software prefetching and streaming. This motivates our approach in developing a pipelined 3D LBM algorithm on the Kalray MPPA processor, based on local memory exploitation and asynchronous communications. Our algorithm is described in every detail and can be used on similar many-core architectures. Our key contributions are as follows:

1) Introduction of a new parallel algorithm for decomposing and streaming 3D stencil domains on local-memory-centric clustered many-core processors, by user-buffers and asynchronous software-prefetching to build a pipelined 3D stencil kernel. The proposed approach is implemented from the LBM compute kernel of OPAL [START_REF] Obrecht | Performance evaluation of an opencl implementation of the lattice boltzmann method on the intel xeon phi[END_REF] and delivers 33% performance gain compared to its original OpenCL code on the Kalray MPPA-256 Bostan many-core processor. 2) This work provides fundamental responses and methods to further domain-decomposition algorithms on clustered many-core processors (2D/3D stencils, image processing). An API proposal is also given in designing user-friendly 2D/3D asynchronous copy on DMA-enabled platforms. 3) Detailed description of the use of generic equations to calculate decomposition indexes dynamically, subdomain dimension and halo size, usable with or without ghost layer as in [START_REF] Mattila | Comparison of implementations of the lattice-boltzmann method[END_REF]. The remainder of this paper is structured as follows. Section 2 presents some related works that are relevant for our contributions. Section 3 introduces the main characteristics of the MPPA architecture and some low-level asynchronous transfer primitives required for building 3D stencils streaming algorithm. Section 4 presents an overview and technical details of the new LBM streaming algorithm using these asynchronous transfers. Experimental results are presented in Section 5, and we conclude in Section 6.

Related work

The straightforward method for implementing LBM is to use two instances of the lattice grid. Collision is carried out on data read from the first grid and propagation consists in writing the new distribution values to the second one. At the next time step, the two grids are swapped and the same procedure is repeated. One-step two-lattice method with collision and propagation fused in a same kernel was first introduced by Massaioli and Amati [START_REF] Massaioli | Achieving high performance in a lbm code using openmp[END_REF]. In the fused kernel, propagation can be done either before (pull scheme) or after collision (push scheme). In spite of its implementation simplicity, the two-lattice method results in substantial memory allocations with large domains.

Most existing LBM implementations on GPU employ the fused two-lattice approach as the easiest and most computationally efficient method. In particular, OpenCL Processor Array LBM (OPAL) from [START_REF] Obrecht | Performance evaluation of an opencl implementation of the lattice boltzmann method on the intel xeon phi[END_REF] implements a one-step two-lattice 3D LBM solver based on the D3Q19 stencil. OPAL is designed to be simple and portable on GPUs, accelerators and other OpenCL-enabled devices.

In the related work on porting a 3D seismic wave propagation on the MPPA processor, Castro et al. [START_REF] Castro | Energy efficient seismic wave propagation simulation on a low-power manycore processor[END_REF] developed a 2D-prefetching algorithm for anticipating data transfers between global memory and local memory. The 3D domain is decomposed in small 2D slices. These slices are copied to the local memory such that transfers overlap with computations. The authors observed important waiting time for data arrival without identifying clearly the DDR bandwidth limitation of the MPPA. The impingement of halo slices on data throughput was not studied either.

Raase and Nordstrom [START_REF] Raase | On the use of a many-core processor for computational fluid dynamics simulations[END_REF] presented a 2D and 3D LBM implementation on Epiphany, a clustered many-core architecture very similar to MPPA. The LBM domain is distributed on 16 cores with user local memory of 24KB per core. Subdomain distribution is done by static mapping of a 4x4 topology on the 16 cores. This 2D mapping is also used on the 3D problem where the third dimension of subdomains is assigned with one global domain dimension, giving rectangular parallelepiped subdomains. These choices allow simulating only very small problem sizes (e.g 12x35x12) and can not be scaled to large simulations. The authors declined using the DRAM memory to implement a streaming algorithm for large LBM domains. Neither memory bandwidth optimization nor possibility of using DMA to perform asynchronous transfer on Epiphany was discussed. In this work, we aim to provide a generic and scalable 3D decomposition with its cuboid distribution function and asynchronous subdomain streaming to reduce data transfer time. Such algorithm can be used as a reference point to implement further high performance LBM or stencil applications on clustered many-core architectures.

Nagar et al. [START_REF] Nagar | Lbm-ib: A parallel library to solve 3d fluid-structure interaction problems on manycore systems[END_REF] implemented a similar cube-based decomposition and distribution function which maps on CPU threads in the shared-memory context of large-memory multi-socket systems. Halo exchange between threads is done by writing directly into the memory zone of the respective cube owners, protected by mutual locks thanks to the CPU cache system. This mechanism cannot be directly used onto clustered architectures like MPPA as it requires either: (1) explicit inter-cluster communications; or (2) committing changes to the global memory then fetched by other clusters. Solution (1) is not relevant in our scope due to (small) local memories, numerous subdomains must be streamed continuously in the MPPA's compute clusters. Such streaming should be done preferably by a self-governing and synchronization-free algorithm. Thus, keeping data in the local memory and waiting for communication does not seem appropriate, not to mention the complexity of managing the inter-subdomain spatial data dependency. In this work, we choose to adapt solution (2), consisting in continuously committing changes of subdomains to the global memory and performing one global synchronization between clusters at each simulation time step.

To the best of our knowledge, there is no work yet on solving the challenges of simulating large LBM domains on clustered many-core architectures. However, large LBM domains cannot fit into on-chip memory and must be stored in the off-chip DDR memory, which has much higher latency. Hence, using DMA to perform asynchronous transfers between off-chip and on-chip memories becomes a key performance factor in order to mask the memory latency. This involves important code re-structuration, as well as new communication primitives and algorithms. In this work, all these problems are addressed and solved while keeping a clear abstraction level from the underlaying target hardware for the sake of genericity. The second generation of Kalray MPPA-256 processor named Bostan (see Fig. [START_REF] Saidi | The shift to multicores in real-time and safety-critical systems[END_REF]) embeds 256 VLIW compute cores grouped into 16 compute clusters (CC) and 16 system cores in two unified I/O subsystems (IOS). The processor delivers peak performance of 634 GFlops in single precision and 317 GFlops in double precision within a consumption of 20 W. Each compute cluster owns 2 MB of local memory (SMEM) shared between 16 user cores (Processing Elements-PEs) running at maximum frequency of 600 MHz. One system core, known as Resource Manager (RM), is reserved for running operating system and resources management. DMA engines of each compute cluster and I/O subsystem provide the system with high bandwidth and low latency transfers between SMEM-SMEM (symmetric inter-cluster) and SMEM-DDR memory (asymmetric cluster-IO).

Kalray OpenCL

The MPPA platform supports OpenCL 1.2 Data Parallel programming model. Such a model exposes each PE as a work-group with 64 KiB of __local memory. As defined in the OpenCL specification, these work-groups cannot share their __local memory to form somehow a large common memory. Moreover, we will see in further discussions that the global three-dimensional domain has to be cut down into smaller subdomains in order to fit in the local memory of clusters. The overhead of copying halo cells of small cuboids is more important than larger ones. In fact, subdomains should be sized to be cubic and as large as possible. Given the maximum 64 KiB of __local memory per work-group, implementing a pipelined 3D LBM algorithm could not be efficient on the current Kalray OpenCL programming model.

On the other hand, in the POSIX programming model, each compute cluster appears as an independent 16-core CPU with 2 MB of shared memory. Thus large buffers can be allocated in the on-chip memories and this multi-core can be programmed using either OpenMP or Pthread multi-threading. This avoids useless data replication for each core and makes it possible to work on a larger common buffer inside the multi-core; therefore, increasing data reuse and reduces halo copy overhead. Such a model is relevant to our scope for efficient on-chip memory usage.

POSIX low-level 3D asynchronous API

In this section, we briefly present some essential primitives performing asynchronous 3D data transfers used to build our pipelined LBM algorithm onto the Kalray MPPA processor. As shown in Fig. 3, the mppa_async_point3d_t type describes copyposition and dimensions of the global and local 3D buffers. The subdomain is represented by width × height × depth elements, whose each element has size in bytes. We take an example to illustrate this specification design. In a common image processing decomposition, one may need to copy a 2D sub-image of 16 × 16 pixels to a larger local buffer, allocated at 18×18 pixels for instance. In this case, users must to explicitly deal with a local stride of two pixels between each data block, since the local buffer is sparse and the data should not be written contiguously to it. This is important when local buffers are declared as true multi-dimensional arrays in the C99 standard, a feature which particularly eases 2D and 3D stencil programming. With the convenient mppa_async_point[2|3]d_t data type (see Fig. 4), arbitrary positions and copy-block dimensions are automatically taken into account inside the 2D/3D put and get functions, facilitating subdomain copy and computation.

A structure mppa_async_event_t is also defined in the API to contain required information for performing an asynchronous transfer. In a put/get function, if the event structure is set, the function fills a pending transaction event and returns immediately (non-blocking paradigm). One can further come back and wait on this event by calling the mppa_async_event_wait() function for job completion. Otherwise, when the event structure is NULL, the function blocks and returns whether the buffer is ready to be reused (put) or the data are received (get). In the following, we take the D3Q19 LBM kernel from OPAL [START_REF] Obrecht | Performance evaluation of an opencl implementation of the lattice boltzmann method on the intel xeon phi[END_REF] as a reference point, from which we propose a generic 3D LBM streaming algorithm with domain decomposition, detailed index and halo size calculation in any configuration. The streaming method is used for updating the whole domain by one time step, then is repeated till the end of simulation duration. While we are focusing on optimizing LBM, our streaming method can also be generalized for other kinds of stencil codes, by adapting the compute kernel, some auxiliary settings (DdQq, halo size, number of time steps T ), and a suitable set of asynchronous transfer primitives (2D/3D).

The first step consists in re-writing the LBM kernel of OPAL from OpenCL-C to a standard C99 code to run on CCs. Given the similarity between OpenCL-C and standard C99, the porting process did not raise much difficulty. The one-step two-lattice method with pull scheme originally implemented in OPAL is re-applied. Two instances of the 3D lattice grid (LatticeEven, LatticeOdd), each containing L x ×L y ×L z nodes, are allocated on the global DDR memory and are accessed in node-wise layout, i.e distribution values of a lattice are stored consecutively. The second step divides the lattice domain into subdomains (see Fig. 5), then copies and computes them one by one on the CC local memory. Each subdomain is defined as a C x × C y × C z cuboid. To avoid repetitions, we use the subscript d as a symbol for the three Cartesian coordinates (x, y, z). Any variable or equation whose variables are subscripted by d should be interpreted as three variables or equations with x-/y-/z-subscripted terms respectively.

For the sake of simplicity, we assume that L d and C d are powers of two and define M d the number of subdomains in each dimension (M d = L d /C d ). The total number of subdomains is the product of the number of subdomains in each dimension M = M x × M y × M z . Besides, we denote the constant F d = C d + h to be the extended subdomain size with halo layers (h) added 1 . Thus, updating a subdomain of C x × C y × C z nodes fetches an extended cuboid F = F x × F y × F z nodes to the local memory. This requirement is true for most cases (non-boundary subdomains -e.g subdomain 4 of Fig. 8). On boundary subdomains (e.g subdomains 0, 1, 2, 3 of Fig. 8), the extended cuboid should be adjusted by applying a halo cutoff to deal with solid nodes. A local subdomain slot must therefore be allocated for F x × F y × F z nodes to match any cases. Algorithm 1 sums up the mono-cluster context where the compute cluster 0 (CC0) is updating M subdomains within an LBM time step. These subdomains are organized in a macro-pipeline using asynchronous 3D put and get functions to overlap computation and communication. We also apply the two-lattice method on local memory, i.e the number of buffer slots is doubled, one for fetching the pre-collision cuboid (S) from the first global lattice grid and one for storing the post-collision cuboid (S ) that will be put in the second lattice grid. The pre-collision cuboid is allocated for F x × F y × F z nodes, while the postcollision cuboid only needs to store C x × C y × C z nodes. Fig. 5 only draws one global lattice grid for compactness, but it should be understood that the local post-collision cuboid will be put in the second grid. These two global grids are then swapped before starting the processing of the next time step.

Ideally, the algorithm should run on multiple compute clusters and exploit all processing cores (PEs) in each cluster (multi-cluster multi-PE). For instance, on MPPA, multithreading within a compute cluster is enabled by spawning up to 15 threads, one per PE, from the PE0 in the POSIX-pthread fashion (create, join). As there are 16 compute clusters available on MPPA, each CC is then responsible for M 16 subdomains. Note that depending on the value of M , there might be K trailing subdomains (K ∈ [0..15]). If K > 0, the algorithm must perform an extra step to copy, update and put back these K trailing subdomains by K compute clusters, while other clusters are waiting. A synchronization barrier at the end of each time step is needed between all CCs to avoid data races at the next time step. This procedure is then repeated as many times as the number of timesteps.

Algorithm 1 Explicit macro-pipeline of 3D stencil updates using double-buffering within a time step. The double-buffering (2-depth) pipeline in Algorithm 1 is the most basic algorithm where communication is overlapped by only one compute-step. As computations are faster than data transfers, deeper pipelines such as triple-or quadruple-buffering (whose details are found in Fig. 6 ) provide better overlapping, but also require more local memory. Note that the time spent in GET and PUT is considered negligible (non-blocking) and transfers are executed in background. However, the time spent in COMPUTE depends on core speed, while the WAIT time depends on how fast the memory system is serving transfer requests and how they are hidden entirely or partially by the COMPUTE function.

In the next sub-sections, we propose methods to solve the following questions that immediately arise from Algorithm 1:

• How can we distribute fairly and exclusively all subdomains across CCs with their proper subdomain-indexing? • Which subdomain size and pipeline depth should we choose to fit with the local memory size and to obtain the best trade-off? • How to manage copy indexes and halo size of any subdomain, with or without using ghost layer?

Subdomain distribution

Given a CC identified by cc id ∈ [0..15], its working subdomains is indexed by a onedimensional range m as cc id × M 16 ≤ m < (cc id + 1) × M 16 (assume K = 0). Mapping bijectively this 1D domain (m) to a 3D one (m x , m y , m z ) for spatial cube indexing (see Fig. 5) was done by space filling curves, such as Morton or Hilbert curves, in [START_REF] Sagan | Space-filling curves[END_REF]. These curves have been efficiently implemented by bit-interleaving in [START_REF] Wise | Ahnentafel indexing into morton-ordered arrays, or matrix locality for free[END_REF] or lookup-table in [START_REF] Thiyagalingam | Is morton layout competitive for large two-dimensional arrays yet?[END_REF]. However, Morton, Hilbert and other curves are better suited for square or cubic grids where the number of elements in all dimensions is equal. In our 3D decomposition scheme, despite the fact that global lattice domain may be cubic (L x = L y = L z ), subdomains may be not (C x = C y = C z ) due to many reasons (see next section), thus these curves are not always suitable for subdomains, as M d can be different.

In order to solve this problem, we implement a simple alternative bijective function f : N → N 3 in Fig. 7. It follows the 3D-row-major layout, which is also a space-filling curve, to index subdomains. Each conversion of the 3D-row-major curve implemented by f takes less than 10 instructions and is as fast as Morton or Hilbert curves.

Local subdomain dimensions

In most of the cases, a cubic subdomain would be ideal for coding and optimizing. However, the local memory of clustered many-core processors is usually limited but plays an important role. On each MPPA's compute cluster, 2 MB local memory is quite small and should also host an embedded operating system, services and the user application binary.

A remaining space of about 1.5 MB is available for dynamic buffer allocations. Some auxiliary variables are also needed in LBM for macroscopic monitoring (velocity, density. . . ). The maximal allocatable space for local pre-collision and post-collision cuboids is around 1.4 MB. Halo copy also consumes memory bandwidth. Hereafter, we refer to "halo bandwidth" HBW as the bandwidth lost in fetching halo layers. The HBW ratio is defined as the quotient of the number of halo cells by the total number of copied cells (main block and halo). On small subdomains, this ratio can be significant. For example, given a cubic subdomain whose main block size is C x × C x × C x , its HBW ratio is:

g(C x , C x , C x ) = (C x + 2) 3 -C 3 x (C x + 2) 3 = O( 1 C x ) lim C x →∞ g(C x , C x , C x ) = 0 Example: g(16, 16, 16) = 18 3 -16 3 18 3 ≈ 0.29 (1)
The best performance is achieved when the volume of the main block (C x × C y × C z ) is maximized and the HBW is minimized. Likewise, local storage should be reduced as much as possible. Let's assume single-precision floating-point representation, applying a D-depth pipeline for the two-local-cuboid method described above must fit into 1.4 MB of local memory and satisfy the linear-programming formulation below:

Find: (D, C x , C y , C z ) Maximize: C x × C y × C z (nodes updated per subdomain) Minimize: F x × F y × F z (per subdomain storage) Minimize: (F x ×F y ×F z )-(C x ×C y ×C z ) F x ×F y ×F z (HBW)
Subject to:

D×((F x ×F y ×F z )+(C x ×C y ×C z ))×19×4 1024 2 ≤ 1.4 F d = C d + 2 ; C d ∈ {2 n } ; D, n ∈ N + (2)
For instance, using D ≥ 3 in order to have better overlapping than with a 2-depth pipeline, restricts to a very small search domain (C d ≤ 128) that can be resolved by running the branch-and-bound algorithm in a script. Solutions can either be (D, C x , C y , C z ) = (3, 16, 8, 16) with 36% HBW ratio or (D, C x , C y , C z ) = (4, 8, 8, 16) with 43% HBW ratio. A permutation of C x , C y , C z also gives other satisfactory solutions, with the same HBW ratio. On the other hand, note that increasing pipeline depth is not relevant, because the higher D is, the smaller (C x , C y , C z ) will be, thus HBW will become unacceptable. Moreover, compute cores will switch between small subdomains more often. The accumulated waiting time will also be more important due to the exponential number of DMA requests and processing overhead of the DDR asynchronous services.

Local and remote copy-index management

In this section, we present generic analytic formulae to process dynamically copy indexing, subdomain size computation and halo cutoff management depending on geometric position of the subdomain. Adding a ghost layer surrounding the computational domain is a common technique to simplify the implementation of the streaming step at boundary cells, see e.g. [START_REF] Mattila | Comparison of implementations of the lattice-boltzmann method[END_REF]. However, we choose not to use this approach in our work, mainly to minimize global memory allocations and avoid wasting bandwidth/storage in moving ghost cells.

However, in our 3D decomposition algorithm, this decision requires careful calculation of copy parameters from subdomain indexes. It is important to note that as the pre-collision cuboid S embeds two additional halo layers for each dimension (F d ), its computational space begins at (1, 1, 1) and ends at (F x -2, F y -2, F z -2) included. When fetching a non-boundary subdomain (main block + halo) from global memory to S, the arrival point of data at the local buffer is set to (0, 0, 0), and the remote point is computed as the global beginning position of the subdomain minus one (back-off) in each dimension ((m

d × C d ) -1).
As ghost layers are not used in our implementation, a boundary subdomain can have up to three missing sides, depending on its location (see Fig. 8). Consequently, the halo layer of these missing sides needs to be pruned from the copied cuboid. The remote read-point and local write-point must also be adjusted as well. In order to generalize the solution, we introduce here three parameters associated respectively to these three adjustments: halo_cutoff, remote_offset and local_offset.

We present in the following, generic formulae which determines copied positions and halo cutoffs of a given 3D cuboid subdomain (m x , m y , m z ), generalized from the 2D representation of Fig. 8.

const A = (A x , A y , A z ) = (0, 0, 0) R = (R x , R y , R z ) = (m x × C x , m y × C y , m z × C z ) B ad = A d + local offset(m d , M d ) B rd = R d + remote offset(m d , M d ) S d = F d + halo cutoff(m d , M d ) (3) The point A = (A x , A y , A z ) = (0, 0, 0) is the start point of the local buffer. The point R = (R x , R y , R z ) = (m x × C x , m y × C y , m z × C z )
is the start point of the remote subdomain, without its halo layers. The fetched cuboid S is sized at

S d = F d + halo cutoff(m d , M d ).
It is read from the remote position B rd = R d + remote offset(m d , M d ) and written to the local position

B ad = A d + local offset(m d , M d ).
The collision is performed on the main block of S and the result is then written to S for the propagation step. However, managing copied parameters of S is simpler than on S. Since S contains exactly the main block of the subdomain, updated data from a collision can be written to (0, 0, 0), which is also the local copied position for sending to remote memory R = (R x , R y , R z ). The parameters halo_cutoff, remote_offset and local_offset are implemented as macros with rules in Tab. 1. local offset remote offset halo cutoff

(from point A) (from point R) (from F d ) m d = 0 1 0 -1 0 <m d && m d < M d -1 0 -1 0 m d = M d -1 0 -1 -1
These position computations can also be applied on other implementations which use ghost layer, by setting all remote_offset to -1, i.e. allowing to jump out of the computational domain, and all local_offset, halo_cutoff to zero, i.e. imposing to copy extended subdomain F d to A d , instead of copying S d to B ad .

Results and discussions

Pipelined 3D LBM stencil on MPPA

We implement the pipelined 3D LBM algorithm on the MPPA-Bostan platform using the POSIX programming model and asynchronous 3D primitives from the MPPA Asynchronous One-Sided library. By default, MPPA-256 cores are set to run at 400 MHz and LP-DDR3 frequency is configured at 1066 MHz, i.e ∼8.5 GB/s peak per DDR. Note that MPPA embeds two DDR interfaces (North and South) and the current OpenCL runtime only uses one DDR and exposes 1 GB of available global device memory, while the MPPA Asynchronous One-Sided library exposes both single and double DDR modes. Different cubic cavity sizes, varying from 64 to 224 are used in our tests, with some exceptions.

Problem sizes larger than 160 can not be run in OpenCL on MPPA due to the 1 GB device memory limit. Local work-group size in OPAL OpenCL is always set to 32 × 1 × 1, as it delivers the best performance in most of the cases.

In single-DDR mode (POSIX and OpenCL), both LatticeEven and LatticeOdd are allocated on the North DDR. In double-DDR mode (POSIX-only), the LatticeEven buffer is allocated on the North DDR and the LatticeOdd is on the South DDR. The effective throughput of the double-DDR mode can be considered as twice as one of the single-DDR mode, thus 2× performance is expected. We present here results of the OPAL kernel rewritten with our new POSIX pipelined algorithm on the MPPA-256, called OPAL async, in 3-depth and 4-depth pipelines and following the local two-lattice method (S and S ) on various cavity sizes. These tests are further run in both single-and double-DDR modes. All these runs are checked for correctness against the original OPAL code on GPU.

As one can notice in Fig. 9, the OPAL async algorithm outperforms the OpenCL version by more than 30% on the single-DDR mode (from 12 MLUPS to 16 ± 1 MLUPS). We also see that the configuration with less HBW (3-depth, 36% HBW) delivers higher performance than the 4-depth configuration (43% HBW). While consuming memory bandwidth, halo cells are copied because of the read-dependency between neighbors. This does not contribute to the final performance. Fig. 9 shows that the less memory bandwidth halo cells take up, the more performance we obtain. This leads to think that the HBW of 2D/3D stencil computations aimed to reach Exascale, like weather forecast, ocean simulation and CFD, should be lessened on future clustered many-core processors. For this to happen, these many-core chips should embed bigger local memory on each compute unit to tear down the useless part of halo exchange due to domain decomposition. Finally, Fig. 9 also shows the expected 2× performance speedup by using two DDRs compared to the single-DDR mode.

Performance extrapolation

For a better understanding of the benefit of our streaming algorithm, we modified the OPAL async code to be able to work with arbitrary values of pipeline-depth. Different pipeline depths were then tried out [START_REF] Succi | The lattice Boltzmann equation: for fluid dynamics and beyond[END_REF][START_REF] Cao | Physical symmetry and lattice symmetry in the lattice boltzmann method[END_REF][START_REF] Obrecht | Performance evaluation of an opencl implementation of the lattice boltzmann method on the intel xeon phi[END_REF][START_REF] Fu | The sunway taihulight supercomputer: system and applications[END_REF][START_REF] Massaioli | Achieving high performance in a lbm code using openmp[END_REF] to see if increasing the number of asynchronous buffers can improve the performance. The block size is thus reduced to 8 × 8 × 8 so that up to eight subdomains can be stored in the local memory. Moreover, instead of using all the 16 compute clusters, we now vary this number of clusters and set the domain size to 128 3 to study the strong scalability of the algorithm. We consider using only the double-DDR mode this time to obtain the best performance.

In Fig. 10, as expected, the 1-depth code (blue line) is slower than other version with communication-computation overlapping. However, we obtain exactly the same performance as the double-buffering case when using more than two buffers [START_REF] Obrecht | Performance evaluation of an opencl implementation of the lattice boltzmann method on the intel xeon phi[END_REF][START_REF] Fu | The sunway taihulight supercomputer: system and applications[END_REF][START_REF] Massaioli | Achieving high performance in a lbm code using openmp[END_REF]. The performance line scales from 1 cluster to 8 clusters, then reaches almost a stable value of between 20-22 MLUPS from 8 clusters to 16 clusters. To explain this, we added the sustained throughput of 3D transfer (red line) from the Kalray unit test dedicated to 3D asynchronous copy. This test only does some ping-pong copies to the DDR and does not perform any calculation (Arithmetic Intensity (AI) = 0 flops/byte). We observe that the native 3D copy reaches the maximum throughput with as few as four clusters (6GB/s), then remains the same for higher numbers of clusters (which is the same trend as the performance of OPAL async.). Four clusters are thus enough to saturate the DDR bandwidth. Unlike the 3D unit test, our LBM code performs real computation on the copied data. Its AI is about 350/(2 * 19 * 4) = 2.3 flops/byte, which means that each CC spends more time working on a 3D data block. This explains in Fig. 10 the MLUPS performance which reaches its upper bound for 8 clusters, instead of 4 clusters of the 3D unit test.

Another precise way to interpret the performance of 20-22 MLUPS is to apply the performance estimation formula presented by McIntosh-Smith et al. [START_REF] Mcintosh-Smith | On the performance portability of structured grid codes on many-core computer architectures[END_REF]:

P = B × 10 9 19 × 2 × 4 × 10 6 (M LU P S) (4)
in which B is the effective memory bandwidth in GB/s. In order to take into account the additional cost of halo copy in our decomposition algorithm, we multiply P by (1 -HBW ), the effective part of bandwidth (main block) which generates the real performance:

P h = 6.0 × 10 9 19 × 2 × 4 × 10 6 × 8 3 10 3 = 20.2 M LU P S (5)
This estimation P h , shows that there is seemingly a little performance gain to perform asynchronous transfers on clustered many-core processors (here MPPA as an example) as for today. This is not because the streaming algorithm is not good, but because the overlapping gain time is too small compared to the lengthy waiting time for data due to the DDR3 bottleneck. This also demonstrates the memory-bound property of general stencil computations and leads to think that newer memory technologies, such as DDR4 and others, will be a performance boost on these architectures. Notice that the scale-down of the 3D throughput versus the peak 17GB/s of two DDRs is caused by the fact that strided copies (2D/3D) must read data from a lot of different DDR memory banks. Furthermore, these copies can unavoidably suffer bad alignments due to the access pattern of application (Q = 19 floats), thus bear an efficiency factor of 3D transfer compared to the linear copy. For instance, on the current MPPA Bostan platform, if the linear transfer factor is normalized at 1, the 3D factor lies often in between 0.35 and 0.42, depending on the copy layout (size of each contiguous block, alignment of strides and dataset). A correlation, computed by the lm function in R, from 1 to 8 clusters gives the performance expectation of our streaming algorithm if we were not bounded by the memory bandwidth (gray line). These results confirm that our pipelined LBM algorithm is strongly scalable, but is quickly memory-bound on MPPA and that its performance heavily depends on the hardware memory bandwidth. Our results also show that the imbalance between computing power and data throughput is one of the largest drawbacks of actual clustered many-core processors, and demonstrate the interest of future high-bandwidth memory technologies.

Conclusions

We introduced a decomposition approach for generic 3D stencil problems with formulations for calculating dynamically copied position indexes, subdomain addresses, subdomain size and halo cells. These analytic results work with or without using ghost layers and are also usable for 2D problems. Based on this decomposition, our new pipelined 3D LBM code outperforms the original OpenCL version by 33 %, by overlapping computation and communication.

We expected that anticipating data requests by asynchronous memory transfers would improve effective throughput and that we could overcome the memory bound of the studied LBM kernel, by introducing enough pipeline depth to hide the global memory access latency. In practice, performance results are still bound by memory bandwidth and increasing the number of buffers (pipeline depth) does not improve performance, as the DDR3 memory is already fully loaded. Moreover, reducing subdomain size to increase pipeline depth induces significant bandwidth consumption for halo copy. Furthermore, the impact of HBW on small local memories was also identified as a governing factor of performance in our algorithm. We found out that the best strategy is to have cubic subdomains as large as possible and that the double-buffering scheme is enough on the current generation of MPPA processor. We furthermore presented comprehensive linear-programming equations which give the best trade-off between these structuring parameters.

In the future, we plan to study a new LBM propagation method which performs inplace lattice update (one-step one-lattice). Such a method will reduce by half the local memory requirement, thus increase the subdomain size, trim down halo bandwidth and improve performance. Porting async_work_group_copy_{2D|3D} primitives to the next OpenCL specification is also under consideration, as this would considerably improve the exploitation of local memory on clustered many-core processors.

Figure 1 .

 1 Figure 1. LBM D3Q19 stencil.

Figure 2 .

 2 Figure 2. MPPA-256 processor overview (Source: Kalray).

Figure 4 .

 4 Figure 4. Illustration of point2d_t datatype for 2D copy. 3D copy is conceived by adding depth and Z fields.

4 .

 4 Pipelined 3D LBM stencil on clustered many-core processors 4.1. Global algorithm

Figure 5 .

 5 Figure 5. 3D LBM/stencil decomposition where Main-block subdomain (green) is copied with its surrounding halo layers (if exists) and one extra subdomain (blue) is needed to store post-collision state.

Figure 6 . 3 -

 63 Figure 6. 3-depth pipeline (triple-buffering) which allows 2-step distance between GET and WAIT, but only 1-step distance between PUT and WAIT, thus the PUT transfer will not be well overlapped (m: index of subdomain to compute, i: index of local buffer slot; G = GET; P = PUT; W = WAIT; C = COMPUTE; WCP = {WAIT + COMPUTE + PUT}; WG = {WAIT + GET}).

Figure 8 .

 8 Figure 8. Local/Remote copied index in 2D (in lattice node) with A: begin of the local buffer = (0,0); R: begin of the remote main block cuboid (without halo); B: begin of the copied cuboid (S), represented by: Ba: index of S on local memory (from A) and Br: index of S on global memory (from R).

Figure 9 .

 9 Figure 9. OPAL async vs. OPAL OpenCL on MPPA for duration = 1000 steps.

Figure 10 .

 10 Figure 10. Performance extrapolation of OPAL async on 8 × 8 × 8 subdomains with the first eight clusters correlation represented by a gray line for 1000 timesteps and cavity size 128.

  * mx = x ; / * o u t p u t s * / * my = y ; / * o u t p u t s * / * mz = z ; / * o u t p u t s * / } . Figure 7. 3D Row-major subdomain-indexing f : N → N 3 .

	v o i d c u b o i d i n d e x 1 t o 3 ( i n t m,	/

* i n p u t * / i n t * mx , i n t * my , i n t * mz ) / * o u t p u t s * / { i n t z = (m / (Mx * My) ) ; i n t y = (m -( z * (Mx * My) ) ) / Mx; i n t x = (m -( z * (Mx * My) ) ) -( y * Mx) ;

TABLE 1 .

 1 COPIED INDEX OFFSET AND HALO CUTOFF OF A SUBDOMAIN.

h = 2 with the D3Q19 stencil.

Acknowledgments

The authors would like to thank some members from the Kalray software team: Nicolas Brunie, Romarik Jodin, Pierre Guironnet de Massas and Clément Leger for interesting discussions and technical support during this work.