Keywords: numerical analysis, second order schemes, hyperbolic problems, residual distribution, C-property, unstructured grids

Residual Distribution for the Shallow Water equations

Résumé : In this paper we consider the discretization the Shallow Water equations by means of Residual Distribution (RD) schemes, and review the conditions allowing the exact preservation of some exact steady solutions. These conditions are shown to be related to both the type of spatial approximation and to the quadrature used to evaluate the cell residual. Numerical examples are shown to validate the theory.

Introduction and generalities

In this paper we study the application of Residual Distribution (RD) discretizations [START_REF] Deconinck | Residual distribution schemes: foundation and analysis[END_REF][START_REF] Abgrall | Residual distribution schemes: Current status and future trends[END_REF] to the Shallow Water equations. Our objective is to determine the conditions under which a high order (RD) scheme is able to preserve exactly particular steady solutions of the problem. The discussion presented here generalizes the work of [START_REF] Ricchiuto | Application of conservative residual distribution schemes to the solution of the shallow water equations on unstructured meshes[END_REF][START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF]. Similar results have been discussed e.g. in [START_REF] Audusse | A well-balanced positivity preserving second order scheme for shallow water flows on unstructured meshes[END_REF][START_REF] Kurganov | A Second-Order Well-Balanced Positivity Preserving Central-Upwind Scheme for the Saint-Venant System[END_REF][START_REF] Xing | High-order finite difference WENO schemes with the exact conservation property for the shallow-water equations[END_REF][START_REF] Xing | High-order well-balanced finite volume WENO schemes and discontinuous Galerkin methods for a class of hyperbolic systems with source terms[END_REF][START_REF] Noelle | High-order well-balanced finite volume WENO schemes for shallow water equations with moving water[END_REF][START_REF] Yee | High-order well-balanced schemes and applications to non-equilibrium flow[END_REF] for finite volume, WENO finite volume, and discontinuous Galerkin schemes. Here, as in [START_REF] Audusse | A well-balanced positivity preserving second order scheme for shallow water flows on unstructured meshes[END_REF][START_REF] Kurganov | A Second-Order Well-Balanced Positivity Preserving Central-Upwind Scheme for the Saint-Venant System[END_REF] the main focus is on two dimensional unstructured mesh approximations and, as in [START_REF] Noelle | High-order well-balanced finite volume WENO schemes for shallow water equations with moving water[END_REF], on problems non necessarily involving trivial equilibria.

Our methodology is based on a weighted residual approach where the element integral of the equation (viz. the element residual) is split into nodal residuals that, once assembled, form algebraic nodal equations. The main result is that, due to the direct use the multidimensional equation to write the algebraic equations, to the continuity of the spatial approximation, and making use of an equivalence between the element integral of any quasi-linear form of the equation and the integral of the conservative form, we are able to preserve exactly steady solutions provided that : the approximation is written directly in terms of a set of steady invariants, exact integration and source term expressions are used. We will show numerical evidence that this is indeed the case, at least whenever the data of the problem are smooth.

Mathematical problem and notation

We consider the numerical approximation of solution to the Shallow Water Equations (SWE) :

∂ t u + ∇ • F (u) + S(u, x, y) = 0 (1)
where, denoting by H the water depth, by q and v the discharge and the local water speed such that q = (q x , q y ) = H v = H(v x , v y ), u is the array of conserved quantities u = H q and F is denotes the conservative fluxes

F = [F x F y ]   q v ⊗ q + g H 2 2 I 2  
with I 2 the 2 × 2 identity matrix, and g the gravity acceleration. The source term S in (1) models the effects of the bathymetry B(x, y) on the mean flow :

S = gH 0 ∇B

Particular steady solutions

By simple manipulations, we can rewrite the SWE as

∂ t u + ∂ x q x + ∂ y q y v(∂ x q x + ∂ y q y ) + H∇E + q ⊥ (∂ y v x -∂ x v y ) = 0 (2)
where E denotes the total energy

E = g(H + B) + v • v 2
and q ⊥ denotes the orthogonal discharge q ⊥ = (-q y , q x ). This particular form of the equations allows to highlight a whole family of exact steady solutions involving homoenergetic, irrotational flow, with a solenoidal discharge field, namely

E = g(H + B) + v • v 2 = constant ∂ y v x -∂ x v y = 0 ∂ x q x + ∂ y q y = 0 (3)
Within this family of steady states we find the following well-known steady solutions :

Lake at rest This solution is completely characterized by the initial uniform value of the state vector

v = E q (4)
This means that ∀t ≥ 0 and ∀ (x, y) we have v = v 0 given by :

H + B = (H + B) 0 = 1 g E 0 = const q x = 0 q y = 0 (5)
Clearl, the last two conditions also imply v = 0, so that all of (3) are verified.

INRIA

Pseudo 1d flow A generalization of the lake at rest solution can be obtained by rewriting the SWE in quasi-linear form in terms of the state vector (4).

After a few manipulations one obtains :

∂ t u + A v ∂ x v + B v ∂ y v + S v (v, x, y) = 0 (6)
with

A v = ∂ v F x and B v = ∂ v F y .
The entries of A v and B v are readily obtained using the chain rule. Note however, that since v depends on the bathymetry, extra terms containg the derivatives of B arise when computing these matrices, so that the expression of the source term is modified. The form of the source term is now

S v (v, x, y) = 1 H - u • u g    0 u -v    q ⊥ • ∇B (7)
These expressions show that a steady solution is given by v = v 0 ∀ t and ∀ (x, y), provided that q ⊥ • ∇B = 0, that is provided that the bathymetry has only one dimensional variations along discharge lines. In particular, without loss of generality in the following we will refer to the pseudo 1d flow solutions as the ones for which

v = v 0 =   E 0 q 0 0   ∀ t and ∀ (x, y)
and

B = B(x) (8)
Note that these solutions basically involve one dimensional flow. Additionally, since in presence of shocks the relevant form of the SWE is the conservative form, relations (2) are no longer valid, so that the pseudo 1d solutions in general do not admit shocks. More general situations in which two different pseudo 1d flows are connected across a shock via the jump conditions are known, but will not be considered here.

Potential flow For completeness we add to the list of steady solutions the case in which the flow speed is obtained by a potential. In particular, let ϕ denote the scalar potential verifying ∆ϕ = 0

If we set v = (∂ y ϕ, -∂ x ϕ), one immediately sees that the second in (3) is always satisfied. Moreover, if we set H = ϕ+H ∞ , with H ∞ a constant, the last in (3) is also automatically true. At last, the first in (3) imposes the form of the bathymetry, finally leading to the following potential solutions :

∆ϕ = 0 v = (∂ y ϕ, -∂ x ϕ) H = ϕ + H ∞ B = 1 g E 0 -ϕ -H ∞ - ∇ϕ 2 2g (9)
In numerical applications one often requires the above steady solutions to be preserved exactly in some sense, so that small perturbations of these solutions can be resolved without the need of excessive mesh refinement. Schemes enjoying this property are often referred to as well balanced. The following more precice characterization can be given.

C-property A numerical scheme that preserves exactly initial solutions of the Lake at rest type (cf. equation ( 5)) is said to verify the C-property. A numerical scheme that preserves an initial lake at rest solution within an accuracy higher than that of its truncation error is said to verify the approximate C-property.

Generalized C-property A numerical scheme that preserves exactly initial solutions of the pseudo 1d flow type (cf. equation ( 8)) is said to verify the generalized C-property. A numerical scheme that preserves an initial pseudo 1d flow solution within an accuracy higher than that of its truncation error is said to verify the approximate generalized C-property.

To our knowledge, so far no scheme is able to preserve exactly potential solutions of type [START_REF] Dunavant | High degree efficient symmetrical gaussian quadrature rules for the triangle[END_REF]. Tha main difficulty lays in the need of preserving exactly the solenoidal and irrotational conditions.

Residual Distribution discretization

We consider discretizations of (1) based on a conservative Residual Distribution (RD) approach [START_REF] Csík | A conservative formulation of the multidimensional upwind residual distribution schemes for general nonlinear conservation laws[END_REF][START_REF] Ricchiuto | Residual distribution for general time dependent conservation laws[END_REF]. In particular, the schemes we analyze are those discussed in [START_REF] Ricchiuto | Application of conservative residual distribution schemes to the solution of the shallow water equations on unstructured meshes[END_REF][START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF]. The basic principles of the RD discretization procedure can be summarized as follows. Let T h be an unstructured triangulation of the two dimentional spatial domain, composed of non ovelapping triangles, h denoting a characteristic elemet size. Let T be the generic triangle of the mesh. On T h , let u h be a continuous piecewise polynomial approximation of u built starting from collocated (nodal) values u i , i ∈ T h . An example of such an approximation is given by standard P k Lagrange finite elements. In particular, we will denote by ψ i the basis function associated to the degree of freedom u i , so that

u h = i∈T h ψ i u i On a generic element T define the fluctuation φ(u h ) φ(u h ) = ∂T F h (u h ) • n dl + T S h (u h , x, y) dx dy (10)
having denoted by F h (u h ) and source h (u h , x, y) continuous numerical approximations of the flux and of the source term, to be defined. Next consider a discretization of the time derivative. Here, we only consider the second order Crank-Nicholson operator

∂ t u + ∇ • F (u) + S(u, x, y) ≈ u n+1 -u n ∆t + ∇ • F (u n+1/2 ) + S(u n+1/2 , x, y)
INRIA with ∆t the time step, and

u n+1/2 = (u n+1 + u n )/2. Wtih this notation owe define the cell residual Φ(u h ) Φ(u h ) = T u n+1 h -u n h ∆t + ∇ • F h (u n+1/2 h ) + S h (u n+1/2 h , x, y) dx dy = T u n+1 h -u n h ∆t dx dy+φ(u n+1/2 h ) (11)
The RD discrete counterpart of ( 1) is obtained by splitting the cell residual to all the degrees of freedom (nodes) of element T . In particular, if Φ j (u h ) is the amount of residual distributed to j ∈ T , then

j∈T Φ j (u h ) = Φ(u h ) (12)
Finally, a local discrete approximation of ( 1) is obtained by requiring

T ∈T h |i∈T Φ i (u h ) = 0 ∀ i ∈ T h (13)

Basic properties

The properties of scheme ( 12)-( 13) depend on the choices made in terms of approximation (viz. interpolation) and distribution (viz. definition of Φ j , given Φ). For a thourough discussion of these properties we refer the reader to [START_REF] Abgrall | Toward the ultimate conservative scheme : Following the quest[END_REF][START_REF] Abgrall | High order fluctuation schemes on triangular meshes[END_REF][START_REF] Csík | A conservative formulation of the multidimensional upwind residual distribution schemes for general nonlinear conservation laws[END_REF][START_REF] Ricchiuto | Application of conservative residual distribution schemes to the solution of the shallow water equations on unstructured meshes[END_REF][START_REF] Abgrall | Simplified stabilisation procedures for residual distribution schemes[END_REF][START_REF] Deconinck | Residual distribution schemes: foundation and analysis[END_REF]. The following two properties are relevant to our discussion.

Conservation For the schemes considered here, conservation is guaranteed by the use of boundary integration of the flux in the computation of the cell residual [START_REF] Csík | A conservative formulation of the multidimensional upwind residual distribution schemes for general nonlinear conservation laws[END_REF]. In particular, under some continuity assumptions on the Φ j s always verified in practice, scheme ( 12)-( 13) verifies a Lax-Wendroff theorem, as long as the discrete approximation of the flux F h used in [START_REF] Kurganov | A Second-Order Well-Balanced Positivity Preserving Central-Upwind Scheme for the Saint-Venant System[END_REF] is continuous across element edges [START_REF] Abgrall | High order fluctuation schemes on triangular meshes[END_REF].

Accuracy For continuous k-th order accurate approximations of u, and of the flux and of the source, scheme ( 12)-( 13) has a truncation error of order O(h k ) provided that whenever u h is the interpolant of a classical solution of (1), then (in 2d) [START_REF] Ricchiuto | Application of conservative residual distribution schemes to the solution of the shallow water equations on unstructured meshes[END_REF][START_REF] Abgrall | High order fluctuation schemes on triangular meshes[END_REF] Φ

j (u h ) = O(h k+1 )
In particular, schemes for which

Φ j (u h ) = β j Φ(u h )
with β j a uniformly bounded distribution matrix have a O(h k ) truncation error [START_REF] Ricchiuto | Application of conservative residual distribution schemes to the solution of the shallow water equations on unstructured meshes[END_REF][START_REF] Abgrall | High order fluctuation schemes on triangular meshes[END_REF].

In all the numerical applications shown in the paper, we will make use of the stabilized limited Lax-Friedrich's RD scheme presented in detail in [START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF]. Even though the analysis presented in this paper remains valid in the general case of P k interpolation, we will focus on a P 1 approximation in space, using the Crank-Nicholson scheme in time.

Preservation of steady solutions

The constraints given so far on the RD discretization allow for the construction of high order conservative schemes. The question we want to answer now is how to ensure the satisfaction of the C-property and/or of its generalized variant. In particular, two degrees of freedom are left : the choice of the spatial approximation of the unknown u, and that of flux and source needed for the evaluation of the fluctuation φ.

To begin with, let us have a look at the discrete equations [START_REF] Ricchiuto | Residual distribution for general time dependent conservation laws[END_REF]. These define a system of nonlinear algebraic equations for the nodal values u i . The solution of such system is obtained by means of some iterative procedure such as a standard Newton algorithm, or as an explicit fixed point algorithm as the one employed in [START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF]. As an example, we can consider the case of a dumped Newton procedure. If we set w = u n+1 , the dumped Newton loop writes

T ∈T h |i∈T j∈T δ ij ω i + ∂ uj Φ i k w k+1 j -w k j = - T ∈T h |i∈T Φ k i ,
where ω i is the dumping parameter. When initializing the loop with w 0 = u n we have on the right hand side, for a high order RD

Φ 0 i = β i Φ 0 = β i   T w 0 h -u n h ∆t dx dy + φ( w 0 h + u n h 2 )   = β i φ(u n h ) ,
so that we end with

T ∈T h |i∈T j∈T δ ij ω i + ∂ uj Φ i k w 1 j -u n j = - T ∈T h |i∈T β i φ i (u n h ) .
Clearly, if φ(u n ) = 0 ∀ T , the Newton algorithm converges in one iteration to the solution w = u n+1 = u n . Hence, our objective is to select the spatial approximation, and the continuous approximation of the flux and of the source, such that

φ(u h ) = ∂T F h (u h ) • n dl + T S h (u h , x, y) dx dy = 0 .

Lake at rest solutions

The preservation of the lake at rest solution is easily achieved in the RD context.

Early results showing the potential of the approach can be found in [START_REF] Brufau | Unsteady free surface flow simulation over complex topography with a multidimensional upwind technique[END_REF]. The approach discussed here is a generalization of the construction proposed in the last reference. The idea is to exploit the fact that we have a set of steady invariants v. We would like to express the discrete equations in terms of differentials of v. In our case, this means that the approximation should be

u h = u(v h ) with v h = i∈T h ψ i g(H i + B i ) q i .

INRIA

On the lake at rest solution we have v h = v 0 . If we now look at the fluctuation φ we have :

φ = ∂T 1 u h q h •n dl+ ∂T   0 g H 2 h 2   n dl+ T 0 gH h ∇B h dx dy = ∂T   0 g H 2 h 2   n dl+ T 0 gH h ∇B h dx dy .
At this point we still have to choose how to represent the bathymetry. Note that, independently on this choice, for exact integration we can apply Gauss' theorem to get

φ = T 0 H h ∇E h dx dy = T 0 H h ∇E 0 dx dy = 0 ,
which shows that Proposition 4.1 (C-property -Exact integration). When using exact integration, a high order RD scheme verifies the C-property, provided that the fluctuation is evaluated using the approximation written in terms of the steady invariants g(H

+ B) q t .
Exact integration is of course never used in practice. A logical choice is to assume for the bathymetry the same variation used for the invariants, namely B h = i∈T h ψ i B i . In this case we recover the constructions of [START_REF] Brufau | Unsteady free surface flow simulation over complex topography with a multidimensional upwind technique[END_REF] and [START_REF] Ricchiuto | Application of conservative residual distribution schemes to the solution of the shallow water equations on unstructured meshes[END_REF][START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF], which are very similar to the hydrostatic reconstruction used e.g. in [START_REF] Audusse | A well-balanced positivity preserving second order scheme for shallow water flows on unstructured meshes[END_REF]. In particular, we have H h = E h /g -B h which belongs now to the same polynomial space used for the approximation. In this case, we can indeed use exact integration with respect to this polynomial variation obtaining, as before, the exact preservation of the lake at rest state. Proposition 4.2 (C-property -Same approximation for depth and bathymetry). When using the same polynomial approximation for H h and B h and exact integration with respect to this assumption, a high order RD scheme verifies the C-property.

The last possibility is that of using for B its exact expression, while evaluating φ with approximate integration. In this case, provided that B is smooth enough, the accuracy obtained will depend on the quadrature formula used in practice. An estimate of the asymptotic behavior of the error can be obtained by studying the explicit fixed point resolution technique used in [START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF]. In this approach, the Newton loop is replaced by the following update

w k+1 i = w k i -ω i T ∈T h |i∈T (β i Φ) k ,
with ω i a relaxation parameter, that can be easily shown to be bounded by O(h -1 ) for stability reasons [START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF]. In the following we assume ω i = O(h -1 ). Suppose now that the bathymetry has a smooth variation, and that we evaluate the line integrals in the fluctuation with formulas exact for polynomials of degree p l , and the surface integrals with formulas exact for polynomials of degree p s .

On the lake at rest solution we have

φ = ∂T   0 g H 2 h 2   n dl + T 0 gH h ∇B h dx dy = T 0 H h ∇E 0 dx dy + max(ǫ p l , ǫ ps ) = max(ǫ p l , ǫ ps ) = max(O(h p l +2 ), O(h ps+3 ))
,

where the last estimates are readily obtained by e.g. bounding truncated Taylor developments of the integrand. The first explicit fixed point iteration will give

w 1 i = u n i + O(h -1 )O (φ(v n h )) = u n i + max(O(h p l +1 ), O(h ps+2 )) .
Starting from the second iteration we have to take into account the time variation of the solution. Neglecting variations on the asymptotic order of φ, we can write

Φ i = β i Φ = β i O |T | ∆t O(w 1 -u n ) + O(φ) = β i O(h) max(O(h p l +1 ), O(h ps+2 )) + max(O(h p l +2 ), O(h ps+3 )) = max(O(h p l +2 ), O(h ps+3 ))
, which easily leads to the result that for the k-th iterate we have

w k+1 i = w k +max(O(h p l +1 ), O(h ps+2 )) = . . . = u n i +max(O(h p l +1 ), O(h ps+2 )) . (14)
This shows that Proposition 4.3 (Approximate C-property). When using approximate quadrature and the exact variation of the bathymetry, high order RD verify the approximate C-property, provided that the bathymetry variations are smooth enough. A coarse estimate of the error on the solution is ǫ = max(O(h p l +1 ), O(h ps+2 )) if the line quadrature formulas used are exact for polynomials of degree p l and the surface quadrature formulas are exact for polynomials of degree p s .

Pseudo 1d flow solutions

In the case of pseudo 1d flow solutions, we are unfortunately not able to provide a simple and general strategy using the same polynomial for the depth and the bathymetry. Things are much easier in a purely one dimensional context, in which our RD approach has very close similarities to the regular/singular residual decomposition approach discussed in [14, ?]. When looking at the problem from a multidimensional perspective, which is a must on unstructured grids, we can however repeat the asymptotic analysis discussed for the lake at rest state.

In particular, when evaluating the fluctuation we set u h = u(v h ), where v is now given by [START_REF] Deconinck | Residual distribution schemes: foundation and analysis[END_REF], with E as in [START_REF] Abgrall | Simplified stabilisation procedures for residual distribution schemes[END_REF]. Also, we set F h = F (v h ), and S h = S(v h , x, y), and use the exact expression of the bathymetry. Since all the quantities are continuous and differentiable within an element (provided INRIA that the bathymetry is), when using exact integration we can use Gauss' theorem and chain rule differentiation to get (cf. section 2.1, equations ( 6) and ( 7))

φ(v h ) = T (∂ t u + A v (u h )∂ x v h + B v (v h )∂ y v h + S v (v h , x, y)) dx dy = T 1 H h - u h • u h g    0 u h -v h    q ⊥ h •∇B dx dy .
Since, by hypothesis we are using the exact representation of the bathymetry, we will have (cf. equation ( 8)) ∇B = (∂ x B, 0). Since q ⊥ h = (0, q x h ), the source term vanishes, and so does the fluctuation. Hence Proposition 4.4 (Generalized C-property -Exact integration). When using exact integration, high order RD verify the C-property, provided that the fluctuation is evaluated using the exact bathymetry and the approximation written in terms of the steady invariants

E q t .
As in the case of the lake at rest solution, exact integration is in practice replaced by numerical quadrature. Also for the pseudo 1d solutions, if the regularity of the bathymetry is high enough, we can estimate the quadrature error, and the perturbation introduced in the numerical solution. The analysis is identical to the one performed for the lake at rest and again leads to the final estimate [START_REF] Noelle | High-order well-balanced finite volume WENO schemes for shallow water equations with moving water[END_REF]. Proposition 4.5 (Approximate generalized C-property). When using approximate quadrature and the exact variation of the bathymetry, high order RD verify the approximate generalized C-property, provided that the bathymetry variations are smooth enough. A coarse estimate of the error on the solution is ǫ = max(O(h p l +1 ), O(h ps+2 )) if the line quadrature formulas used integrate exactly polynomials of degree p l and the surface quadrature formulas are exact for polynomials of degree p s .

Numerical examples

This section provides a numerical verification of the theory discussed. Results showing the respect of the C-property when using the same approximation for the depth and for the bathymetry (cf. proposition 4.2) have been published already in [START_REF] Ricchiuto | Application of conservative residual distribution schemes to the solution of the shallow water equations on unstructured meshes[END_REF][START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF]. We refer to these references for a thorough numerical assessment of this approach.

Here we are more interested in the numerical verification of propositions 4.3 and 4.5. In particular, our main interest is the verification of the estimate [START_REF] Noelle | High-order well-balanced finite volume WENO schemes for shallow water equations with moving water[END_REF]. To do this we will consider a simple test case involving a pseudo-one dimensional variation of the bathymetry. All the test cases involve the spatial domain [0, 25] 2 , over which we impose three different variations of the bathymetry : Case T1 In this case we take B = B(x), with given by

B(x) = 0.2 -0.05(x -10) 2 if x ∈ [8, 12] 0 otherwise .
Case T2 In this case we take B = B(x) given by

B(x) = 0.2 sin(0.25π(x -8)) 4 if x ∈ [8, 12] 0 otherwise .
Case T3 In this case we take B = B(x) given by

B(x) = 0.2 sin(0.25π(x -8)) 6 if x ∈ [8, 12] 0 otherwise .
The objective is to have bathymetries with different regularity in order to be able to test the influence of this parameter on the validity of ( 14). All the computations have been run on unstructured grids as the one reported in figure 1. The unstructured nature of the meshes introduces a multidimensional character in the tests, since no preferential direction is present. In particular, all the grid refinement studies have been performed by conformally refining the mesh of figure 1. Three levels of refinement have been considered. The coarsest mesh size is h ≈ 25/40.

We have considered three quadrature strategies :

Strategy Q1 We use 2 line quadrature points corresponding to the trapezium rule, and 1 surface quadrature point, corresponding to mid-point quadrature. Both formulas integrate exactly polynomials of degree 1, hence the expected accuracy, according to ( 14) is

ǫ = O max(h 2 , h 3 ) = O(h 2 ) .
Strategy Q2 We use a standard 2 points Gaussian line quadrature formula, and a 6 points surface quadrature formula. The Gaussian formula integrates exactly polynomials of degree 3, while the surface formula integrates exactly polynomials of degree 4. Hence the expected accuracy, according to ( 14) is

ǫ = O max(h 4 , h 6 ) = O(h 4 ) .
Strategy Q3 We use a standard 3 points Gaussian line quadrature formula, and a 6 points surface quadrature formula. The Gaussian formula integrates exactly polynomials of degree 5, while the surface formula integrates exactly polynomials of degree 4. Hence the expected accuracy, according to ( 14) is

ǫ = O max(h 6 , h 6 ) = O(h 6 ) .
The 6 point surface quadrature formula has been taken from [START_REF] Dunavant | High degree efficient symmetrical gaussian quadrature rules for the triangle[END_REF]. We consider both an initial lake at rest state, and a pseudo 1d flow state. In both cases the initial value for the total energy is taken as E 0 = 22.06605.

In the lake at rest solution we obviously set q 0 = 0, while for the pseudo 1d solution we take q 0 = 4.42. Time dependent computations have been run with the LLFs scheme described in detail in [START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF]. In particular, we set the initial solution to the exact steady one and let the scheme compute an unsteady solution until the final The results for the lake at rest case are reported on figures 2, 3, and 4. In the figures we plot the grid convergence of the error for the cases T1, T2, and T3, comparing the three quadrature strategies with the prediction of [START_REF] Noelle | High-order well-balanced finite volume WENO schemes for shallow water equations with moving water[END_REF].

From figure 2 we see that all the strategies fail to reach even the level of truncation error of the scheme, which is second order. This is definitely due to the lack of regularity of the bathymetry. Indeed, in the case T1, the function B(x) is only continuous, hence all formulas basically give first order of accuracy. The slope observed is indeed between 1 and 2, as shown in the picture.

More interesting is the result of figure 3. In this case, the slopes observed are in perfect agreement with [START_REF] Noelle | High-order well-balanced finite volume WENO schemes for shallow water equations with moving water[END_REF] for the strategies Q1 and Q2. For the quadrature Q3 we only observe fourth order of convergence, instead of the sixth order predicted by [START_REF] Noelle | High-order well-balanced finite volume WENO schemes for shallow water equations with moving water[END_REF]. However, in this case the bathymetry only has three continuous derivatives, which explains why we cannot obtain more than fourth order of accuracy, which is indeed what we get with the strategy Q3.

Finally, the results for the T3 test case, shown in figure 4, confirm the validity of [START_REF] Noelle | High-order well-balanced finite volume WENO schemes for shallow water equations with moving water[END_REF]. In this case, the bathymetry has enough continuous derivatives to achieve the theoretical sixth order of convergence when using the Q3 strategy, confirming the validity of proposition 4.3. The results obtained for the pseudo 1d flow problem are identical in nature to the ones discussed for the lake at rest state. The grid convergence plots comparing the different quadrature strategies are reported in figures 5, 6, and 7.

As before, the bathymetry in T1 is only continuous, and the accuracy observed is between one and two, as we can see on figure 5. As we can see on figure 6, the scheme does better in case T2, in which we can get both the second and fourth orders of accuracy predicted by [START_REF] Noelle | High-order well-balanced finite volume WENO schemes for shallow water equations with moving water[END_REF] for the quadrature strategies Q1 and Q2. In the case of quadrature Q3 we get again a fourth order slope, due to the lack of regularity of B(x). The fact that for the same bathymetry we get the same behavior in the lake at rest and in the pseudo 1d flow case is encouraging. It confirms our theoretical analysis that the error is only related to quadrature, hence to the regularity of the data.

At last, we can see on figure 7 that also for the pseudo 1d flow solution when the data is regular enough, equation ( 14) predicts correctly the asymptotic error, thus confirming the validity of proposition 4.5.

INRIA

Final remarks

In this paper we have analyzed the satisfaction of the C-property and of the generalized C-property when using high order residual distribution schemes to discretize the shallow Water equations. We have provided simple arguments to prove that provided that the approximation is written in the correct set of variables, one can obtain at least the approximate version of these properties on arbitrary triangulations, provided that the exact bathymetry is used when discretizing the formulas. An estimate of the asymptotic error has been given and confirmed by numerical computations.

From the practical point of view the need of using the exact data is a drawback, since these data are often not available in the form of an equation. Moreover the dependence of the accuracy on the regularity or the bathymetry is also a limitation, the variation of the bed slope being often quite irregular.

However, we judge our exercice interesting in the sense that it points out in 2d the need of a truly multidimensional approach, as well as the need of more information. In particular, in the case of the pseudo 1d flow, the form of the source term in [START_REF] Csík | A conservative formulation of the multidimensional upwind residual distribution schemes for general nonlinear conservation laws[END_REF] suggests that the variations of the bathymetry in the cross stream direction are important and, hence, a truly multidimensional treatment is needed if one wants to resolve better these non-trivial equilibria.

Other approaches in which the bathymetry is reconstructed with polynomials of higher degree than the solution, or in which one adopts a local stream-aligned reference frame to decompose different componedns of ∇B are of course also possible. The real challenge is also to be able to handle more complex 2d solutions such as the potential solution [START_REF] Dunavant | High degree efficient symmetrical gaussian quadrature rules for the triangle[END_REF]. All these topics will be explored in the future.

Aknowledgements

The author would like to express his gratitude to the organizers of the conference "Numerical approximations of hyperbolic systems with source terms and applications" for their warm welcome. The discussions with Remi Abgrall, Andreas Bollermann, Smadar Karni and Sebastian Noelle have been very helpful in the development of this work.

INRIA

Figure 1 :

 1 Figure 1: Unstructured triangulation used in the numerical tests

Figure 2 :Figure 3 :

 23 Figure 2: Grid convergence for the case T1, lake at rest with E 0 = 22.06605

Figure 4 :

 4 Figure 4: Grid convergence for the case T3, lake at rest with E 0 = 22.06605. The dashed line next to each convergence curve represents the corresponding slope as predicted by equation (14).

Figure 5 :Figure 6 :

 56 Figure 5: Grid convergence for the case T1, pseudo-1d flow with q 0 = 4.42 and E 0 = 22.06605

Figure 7 :

 7 Figure 7: Grid convergence for the case T3, pseudo-1d flow with q 0 = 4.42 and E 0 = 22.06605. The dashed line next to each convergence curve represents the corresponding slope as predicted by equation (14).

RR n 7191


Brice Boyer

Jean-Guillaume Dumas
email: .dumas\@imag.fr

Pascal Giorgi
email: pascal.giorgi@lirmm.fr

\{brice Boyer

Claude Shannon

Exact Sparse Matrix-Vector Multiplication on GPU's and Multicore Architectures

We propose different implementations of the sparse matrix-dense vector multiplication (SpMV) for finite fields and rings Z /m Z. We take advantage of graphic card processors (GPU) and multi-core architectures. Our aim is to improve the speed of SpMV in the LinBox library, and henceforth the speed of its black box algorithms. Besides, we use this and a new parallelization of the sigma-basis algorithm in a parallel block Wiedemann rank implementation over finite fields.

Introduction

Nowadays, personal computers and laptops are often equipped with multicore architectures, as well as with more and more powerful graphic cards. The latter ones can be easily programmable for a general purpose computing usage (Nvidia Cuda, Ati Stream, OpenCL). Graphic processors offer nowadays quite frequently superior performance on a same budget as their CPU counterparts. However, programmers can also efficiently use many-core CPUs for parallelization e.g. with the OpenMP standard.

On the numerical side, several libraries automatically tune the sparse matrix kernels [START_REF] Vuduc | Oski: A library of automatically tuned sparse matrix kernels[END_REF][START_REF] Vuduc | Fast sparse matrix-vector multiplication by exploiting variable block structure[END_REF][START_REF] Tvrdik | A new approach for accelerating the sparse matrix-vector multiplication[END_REF] and recently some kernels have been proposed e.g. for GPU's [START_REF] Vazquez | The sparse matrix vector product on GPUs[END_REF][START_REF] Buluç | Parallel sparse matrix-vector and matrix-transpose-vector multiplication using compressed sparse blocks[END_REF][START_REF] Bell | Implementing sparse matrix-vector multiplication on throughput-oriented processors[END_REF]. In this paper we want to adapt those techniques for exact computations and we first mostly focused on Z /m Z rings, with m smaller that a machine word.

The first idea is to use the numerical methods in an exact way as has been done for dense matrix operations [START_REF] Dumas | Finite field linear algebra subroutines[END_REF]. For sparse matrices, however, the extraction of sparse matrices is slightly different. Also, over small fields some more dedicated optimizations (such as a separate format for ones and minus ones) can be useful. Finally, we want to be able to use both multi-cores and GPU's at the same time and the best format for a given matrix depends on the underlying architecture.

Therefore, we propose an architecture with hybrid data formats, user-specified or heuristically discovered dynamically. The idea is that a given matrix will have different parts in different formats adapted to its data or the resources. Also we present a "just-in-time" technique that allows to compile on the fly some parts of the matrix vector product directly with the values of the matrix.

We have efficiently implemented 1 "Sparse Matrix-Vector multiplication" (SpMV) on finite rings, together with the transpose product and iterative process to compute the power of a matrix times a vector, or a sequence of matrix products.

We also make use of this library to improve the efficiency of the block Wiedemann algorithm's of the LinBox2 library. Indeed, this kind of algorithm uses block "black box" [START_REF] Kaltofen | Factoring high-degree polynomials by the black box Berlekamp algorithm[END_REF] techniques: the core operation is a matrix-vector multiplication and the matrix is never modified. We use the new matrix-vector multiplication library, together with a new parallel version of the sigma-basis algorithm, used to compute minimal polynomials [START_REF] Giorgi | On the complexity of polynomial matrix computations[END_REF][START_REF] Dumas | Parallel computation of the rank of large sparse matrices from algebraic k-theory[END_REF].

In section 2 we present different approaches to the parallelization of the SpMV operation, with the adaptation of numerical libraries (section 2.3), new formats adapted to small finite rings (section 2.4) together with our new hybrid strategy and their iterative versions (section 2.5). Then in section 3 we propose a new parallelization of the block Wiedemann rank algorithm in LinBox, via the parallelization of the matrix-sequence generation (section 3.1) and the parallelization of the matrix minimal polynomial computation (section 3.2).

Sparse-Vector Matrix multiplication

We begin with introducing some notations. In this section, we will consider a matrix A; the element at row i, column j is A[i, j]. The number nbnz is the number of non zeros elements in matrix A, it has row lines and col columns. If x and y are vectors, then we perform here the operation y ← A x + y. The general operation y ← αA x +β y can then be done by pre-multiplying x and y by α and β respectively. The apply operation in black box algorithms, or y ← A x, is performed by first setting y elements to zero. For further use in block methods, we also provide the operation Y ← αA X +β Y where X and Y are sets of vectors.

Sparse Matrix Formats and Multiplication

Sparse matrices arise form various domains and their shapes can be very specific. Taking into consideration the structure of a sparse matrix can dramatically improve the performance of SpMV. However, there is no general storage format that is efficient for all kind of sparse matrices.

Among the most important storage formats is the COO (coordinate) format which stores triples. It consists of three vectors of size nbnz, named data, colid and rowid, such that data

[k]= A [rowid[k],colid[k]].
The CSR (compressed storage row) stores more efficiently the previous representation: the rowid field is replaced by a (row + 1)

long start vector such that if start[i] k < start[i + 1], then data[k]= A[i,colid[k]].
In other words, start indicates where a line starts and ends in the other two ordered fields.

The ELL (ELLpack) format stores data in a denser way: it has data and colid fields such that data[i, j 0 ]= A[i,colid[i, j 0 ]], where j 0 varies between 0 and the maximum number of non zero elements on a line of A. One notices that these fields can be stored in row-major or column-major order. A variant is the ELL_R format that adds a row long rownb vector that indicates how many non zero entries there are per line.

The DIA (DIAgonal) is used to store matrices with non zero elements grouped along diagonals. It stores these diagonals in an array along with offsets where they start. We refer to [START_REF] Bell | Implementing sparse matrix-vector multiplication on throughput-oriented processors[END_REF], [START_REF] Vazquez | The sparse matrix vector product on GPUs[END_REF] for more details on these formats.

This very schematic description of a few well-known formats shows that each of them has pros and cons. Our aim is to produce a more efficient implementation of the SpMV operation on finite fields than the one present in LinBox, taking first advantage of this variety of formats.

Finite field representation

We present now how the data is stored. We use data types such as float, int. Firstly, when doing modular linear algebra, we try to minimize the number of costly f mod (reduction) operation calls. For instance, we prefer if possible the left loop to the right one in the next figure: . The latter representation enables us to perform twice more operations before a reduction, but this reduction is slightly more expensive. Another trade-off consists in choosing a float representation instead of double (on the architectures that support double). Indeed, operations can be much faster on float that on double but the double representation lets us do more operations before reduction. This is particularly true on some GPU's. In figure 1, we present variations in efficiency due to the storage data type and the size of m on one core of a 3.2GHz Intel Xeon CPU and a Nvidia GTX280 GPU. The timings correspond to the average of 50 SpMV operations, where x and y are randomly generated on the CPU (which also takes into account every data transfer between the CPU and GPU). The measures corresponds to the number of million floating point operations per seconds (flops); a SpMV operation requires 2 * nbnz such operations. The performance correspond to the best ones achieved on the matrices 3 For further information about these techniques on these rings and fast arithmetic in Galois extensions, see e.g. [START_REF] Dumas | Finite field linear algebra subroutines[END_REF].

Adapting numerical libraries

Another speed-up consists in using existing numerical libraries. The ideas behind using them on the rings Z /m Z, is twofold. Firstly, we delay the modular reduction, secondly we can use highly optimized popular libraries and get instant speed-ups as compared to more naive self-written routines.

Just like BLAS libraries can be used to speed up modular linear algebra [START_REF] Dumas | Dense linear algebra over wordsize prime fields: the FFLAS and FFPACK packages[END_REF], we can use numerical libraries for our purposes, or get inspiration for our algorithms from their techniques. For instance, there is the Oski library [START_REF] Vuduc | Oski: A library of automatically tuned sparse matrix kernels[END_REF] for sequential numerical SpMV, or the GPU implementation of SpMV by Nathan Bell et al. in [START_REF] Bell | Implementing sparse matrix-vector multiplication on throughput-oriented processors[END_REF]. The BLAS specifications include sparse BLAS4 but they are seldom fully implemented on free BLAS implementations.

Unfortunately, they usually cannot be used as-is. We need to extract submatrices from the sparse matrices, which is more complicated than for its dense counterpart when the use of strides and dimensions suffices. For instance, if one can do b accumulations on y[i] before reducing and the line i of A has r i non zero elements. Then we want to split this line between ⌈r i /b⌉ matrices. Finally, we can use the numerical libraries on these submatrices we have created. The general algorithm reads like follows:

spmv(y,A,x){ foreach submatrix Ai in A do{ spmv_num(y,Ai,x); reduce(y,m); } }

New formats

Most of the formats implemented show a row-level parallelism, except COO that has element-wise parallelism. The COO case is not obvious to implement and generally much slower. The parallel efficiency of other formats will depend then on the length of the rows as well as the data regularity. Unbalanced rows on a GPU architecture will produce many idle threads. Two solutions exist: the vector approach of Bell (they split the rows into shorter chunks) or the rearranging of rows with permutations to sort the row according to their length. The latter will not work in e.g. a power distribution of the row lengths. The ELL format answers very well this problem because each row has the same length.

An other way to parallelize the SpMV operation is to split the matrix A along rows to get smaller submatrices and treat them in parallel. We took this approach on the CPU COO algorithm. Also, we are dealing with large matrices, used many times as black-boxes. Therefore there is a trade-off between the time spent on optimizing the matrix and how much faster these optimizations will make SpMV run.

Things to consider during preprocessing may include for instance: reordering row-columns to create denser parts, choosing best-fitting formats, cutting the matrix into efficient sub-matrices ( [START_REF] Vuduc | Fast sparse matrix-vector multiplication by exploiting variable block structure[END_REF], [START_REF] Tvrdik | A new approach for accelerating the sparse matrix-vector multiplication[END_REF]). . . The preprocessing approach is taken by Oski: if the expected number of SpMV is very high, optimizing the matrix deeper will prove efficient.

Base case: JIT

One idea to improve SpMV on a given matrix is to hard code this operation in a static library. We read the matrix file and create a library that will apply this matrix to input vectors. For instance the y ← y +A x operation on the matrix 2 1 0 3 would be translated to (if m = 27)

void spmv(float * y, const float * x) { y[0] += 2*x[0] ; y[0] += x[1] ; y[0] = fmod(y[0],27); y[1] += 3*x[1] ; y[1] = fmod(y[1],27); }
Then we compile this generated file a as static library and use dlopen to access its functions. As we can see in this example, one can implement various optimizations: rearranging the rows so that the work is more even (non implemented yet), replacing the occurrences of ±1 in the matrix by less costly additions or subtractions. We have better control on what the compiler will produce. However, large matrices take extremely long to compile. gcc cannot compile the library if the source holds in one huge file, so we divide the matrix into parts of 1000 non zero elements and compile them. Only then for instance, we could compile bibd_81_3 but it takes 63s on the same Xeon machine. Once it is compiled, the CPU version runs at 620 Mflops, which is reasonably fast.

Taking into account the ±1

The example of JIT and the observation that many matrices arising from different applications have a lot of ±1 F tends to draw our attention to this special case. Moreover, many matrices on a small fields also share this property. Thus we can extract two submatrices corresponding to the 1 and -1 from the rest of the matrix and replace multiplications by usually less expensive additions. Besides, the data field in most formats (except ELL, DIA) can be forgotten as we know they only consist of 1 or -1: this reduces the memory usage. Doing only additions as opposed to axpy also hugely delays reduction.

!"#

Figure 3: Speed improvement on one 3.2GHz Intel Xeon CPU and a Nvidia GTX280 GPU when segregating or not the ±1 Figure 3 shows a maximum 20% improvement on a matrix with only 1s and 15% on matrix with 50% of ±1.

Basic Formats

As evoked earlier, the matrix A can be split into smaller submatrices. These submatrices can have a format adapted to them and/or can be treated differently. For instance, we can split row-wise and distribute these matrices for parallelism, or split them column-wise as in the delaying case (figure 2. This makes (possibly) many matrices that we each want to optimize individually so we get better overall performance.

We start with some observations. The COO format is slow due to the many fmod calls, it is best used when the matrix is extremely sparse. The CSR format is denser and can let delayed reduction occur, but one has to ensure the row lengths are well balanced when parallelizing. The ELL formats are very efficient on matrices that have roughly the same number of non zeros per line. The ELL_R format ( [START_REF] Vazquez | The sparse matrix vector product on GPUs[END_REF]) is better for uneven rows lengths. One difference in the CPU and GPU architecture makes the ELL row-major on the CPU (for better cache use) and column-major on the GPU (for better coalescing). The following figure shows on one example (bibd_81_3) the variation of efficiency. The data is normalized so that CSR is 1 on the CPU or GPU.

Figure 4: Speed-ups for various formats on matrix bibd_81_3 both on one 3.2GHz Intel Xeon CPU and a Nvidia GTX280 GPU; reference is CSR on each architecture

Hybridization

The previous remarks lead us to combine these formats to take advantage of them. A hybrid format such as ELL(R)+COO or ELL(R)+CRS leads to good performance on the GPU. When the ELL part is taken out of a matrix, many rows can be left empty. Then, we use a format called COO_S that is a CSR format with pointers only to the non empty rows. It has data, colid same as in CSR and COO. The number rowid[k] corresponds to the k th non empty row that starts in data and colid at start[k]. This format could be avoided if we used row permutations and ordered the lines according to their weight.

Heuristic format chooser

The previous remarks show a great complexity in the formats and the cutting of the matrix. We have implemented a user-helped heuristic format chooser. For instance, the user can indicate if she wants to try and make use of ±1. If so, for each submatrix, the program tries to find an a priori efficient format for them or if it fails, does not separate the 1 or the -1 from the rest. She can also indicate what is the format she wants to fill in priority.

The hybridization of the matrix is usually done as follows. If the matrix is large enough and most of the lines are filled, it will try to fit a part of the matrix in an ELL or ELL_R format. This choice is supported by the observation that many matrices have a c + r row distribution where c is some constant and r ∈ Z varies and the fact that ELL is generally much faster that other formats for matrices with even row weight. The rest of the matrix will be put in a CSR, COO or COO_S format, according to the number of empty lines and the number of residual non zero elements. Parameters that decide when segregating the 1s, that choose the best length for ELL matrix, etc., vary according to the architecture of the computer and need some specific tuning. This tuning is not yet provided at compile time but some of it could be automatically performed at install time.

Experiments show that this heuristic often gives equal or better results that simple formats on the CPU and the GPU.

Block and iterative versions

Using multi-vectors

We have described the SpMV operation y ← A x where x and y are vectors. We also need x and y to be multi-vectors, for they may be used for block algorithms. There are at least two ways of representing them : row or column-major order. In the row-major order, we can use the standard SpMV many times (and align the vectors). In the column-major order, we can write dedicated versions that try and make use of the cache. Indeed, in this case, we traverse the matrix only once and x and y are read/written contiguously. On figure 5, we note that on the CPU, using column-major multivectors is a non negligible gain of speed. On the contrary, the GPU implementation fails to sustain good efficiency for blocks of more than 8 vectors and some large matrices start to reach the memory limit.

Performance issues

The GPU operation on a single SpMV call from the host point of view is very slow because we need to move the vectors between the host and the device. It is therefore only usable on operations that need no data moving between the host and the device. Examples include the computation y ← A n x or the 3 Parallel block Wiedemann algorithm Some of the most representative applications requiring efficient sparse matrixvector product are blackbox methods based on the Lanczos/Krylov approach. In particular, the method proposed by Wiedemann [START_REF] Wiedemann | Solving sparse linear equations over finite fields[END_REF] and its block version proposed by Coppersmith [START_REF] Coppersmith | Solving homogeneous linear equations over GF (2) via block Wiedemann algorithm[END_REF] are well suited to highlight efficiency of sparse matrix-vector product since the latter is quite often their bottleneck.

As an application, we propose to improve the implementation of the Block Wiedemann rank algorithm presented in [START_REF] Dumas | Parallel computation of the rank of large sparse matrices from algebraic k-theory[END_REF]. Let us first briefly recall the outline of this algorithm, we let the reader refer to e.g. [START_REF] Turner | A block Wiedemann rank algorithm[END_REF] for further details.

Let A ∈ F n×n be a matrix satisfying the preconditions of [START_REF] Kaltofen | On Wiedemann's method of solving sparse linear systems[END_REF]. Then the algorithm can be decomposed in three steps:

1. Compute the matrix sequence S i = Y T A i Y for i = 0..2n/s + O(1), with Y ∈ F n×s chosen at random 2. Compute the minimal matrix generator F A Y ∈ F s×s [x] of the matrix series S(x) = i S i x i 3. Return the rank r = deg(det(F A Y )) -codeg det(F A Y ).
Our approach is to separate the parallelization of each step. The first step is clearly related to sparse matrix-vector product and we will re-use our tools presented in previous sections. The second step needs the computation of a minimal matrix generator. This can be achieved by a σ-basis computation as explained in [8, section 2.2]. Finally, the last step reduces to computing the co-degree of the determinant of the σ-basis . The degree of the determinant being directly computed as the sum of the row degrees of F A Y since, due to the σ-basis properties, the matrix is already in Popov form.

Parallelization of the matrix sequence generation

The parallelization proposed in [START_REF] Dumas | Parallel computation of the rank of large sparse matrices from algebraic k-theory[END_REF] was to ship independent set of vector blocks of V to different cores and apply them in parallel. Then gather the results to compute the dense dot products by U T .

An alternative is to use the SpMV library and let it take care of the iteration with the algorithm of the preceding section.

In figure 7 we compare both approaches:

Figure 7: Speed up from the new SpMV library compared to the native LinBox implementation in the generation of the matrix sequence (2n iterations) on one core of a 2.33GHz Intel Xeon E5345 CPU

Parallelization of the σ-basis computation

One can efficiently compute σ-basis using the algorithm PM-Basis of [START_REF] Giorgi | On the complexity of polynomial matrix computations[END_REF]. This algorithm mainly reduces to polynomial matrix multiplication. Therefore a first parallelization approach is to parallelize the polynomial multiplication.

Parallel polynomial matrix multiplication

Let A, B ∈ F n×n [x] be two polynomial matrices of degree d. One can multiply A and B in O(n 3 d + n 2 d log d) operations in F assuming F has a d-th primitive root of unity [START_REF] Cantor | On fast multiplication of polynomials over arbitrary algebras[END_REF]. Assuming one has k processors such that k n 2 , one can perform this multiplication with a parallel complexity of O(

n 3 d k + n 2 d log d k
) operation in F. Let us now see the sequential fast polynomial matrix multiplication algorithm and how it achieves such a parallel complexity:

Fast Polynomial Matrix Multiplication:

Inputs: A, B ∈ F n×n [x] of degree d, ω a d-th primitive root of unity inf F. Outputs: A × B 1. Ā := DF T (A, [1, ω, ω 2 , ..., ω 2d ]) 2. B := DF T (B, [1, ω, ω 2 , ..., ω 2d ]) 3. C := Ā ⊗ B 4. C := 1 2d DF T ( C, [1, ω -1 , ω -2 , ..., ω -2d ]) return C.
Here, DF T (P, L) means the multi-points evaluation of the polynomial P on each points of L, while ⊗ means the point-wise product.

• step 1,2 and 4 can be accomplished by using Fast Fourier Transform on each matrix entries which gives n 2 × O(d log d) operations (see [START_REF] Gathen | Modern Computer Algebra[END_REF]Theorem 8.15]). This clearly can be distributed on k processors such that each processor achieves in parallel the FFT on n We report in figure 8 the performance of the implementation of this parallel algorithm in the LinBox 5 library. Our choice of using this parallel algorithm rather than another one, achieving a possible better parallel complexity, has been driven by the re-usability of efficient sequential components of the library (e.g. matrix multiplication) and the ease of use within the library itself (i.e. mostly the same code as sequential one, only some OpenMP pragmas have been added).

One can see on figure 8 that our coder does not completely match the theoretical parallel speedup. The best we can achieved with 16 processors is a speedup of 5.5, which is only one third of the theoretical optimality. Nevertheless, one can see that with less processors (e.g. less than 4) the speedup factor is closer to 75% of the optimality, which is quite fair. We think this phenomenon can be explained by the underlying many multi-core architecture (Quad-Core AMD Opteron), which may clearly suffers from cache effect if computation are done on same chip or not.

As expected, we can also point out from figure 8 that our implementation benefits at most from parallelism when matrices are larger. Since workload on each core is more important, this allows to hide the penalty from memory operations and threads management of OpenMP. This remarks also applies on the degree but the impact is less important.

Parallel σ-basis implementation

According to the reduction of PM-Basis to polynomial matrix multiplication, one can achieve a parallel complexity of O˜(

n 3 d k + n 2 d log d k
) operations in F with k processors for σ-basis calculation, assuming k n 2 . Therefore, it suffices to directly plug in our parallel polynomial matrix multiplication into the original code of the LinBox library to get a parallel σ-basis implementation.

We report in figure 9 the performance of the parallel version of PM-Basis algorithm within LinBox. Here again, the speedup factor of parallelism is quite low when compared to the theoretical optimality. At most we were able to obtain a speedup of 3 with 16 processors. However, this timings are consistent with the previous ones in figure 8 where the best speedup was 5.

One may notice that reduction to polynomial matrix multiplication of the PM-Basis algorithm relies on a divide a conquer approach on the degree of the approximation (see [11, theorem 2.4]). Therefore, the recursion calls are made with smaller and smaller approximation's degrees, which leads to use less efficient parallel multiplications. Moreover, when the degree is too small the use of the M-Basis algorithm of [START_REF] Giorgi | On the complexity of polynomial matrix computations[END_REF] should be prefered since it becomes more efficient in practice. We have not yet implemented a parallel implementation of this algorithm in LinBox and this clearly affects the performance of our implementation.

Parallel determinant co-degree

Here we just launch in parallel the evaluations of the matrix polynomial at different points, and the computation of the determinant of the obtained matrix at the given point, and gather the results sequentially with the Poly1CRT class of Givaro.

Parallel block Wiedemann performance

In table 2 we show the overall performance of our algorithm on an octo-processor Xeon E5345 CPU, 8 × 2.33GHz. *-LB shows the timings of the current LinBox implementation, where *-SpMV presents our new improvement, both in sequential and in parallel. The speed-up for SpMV between 1 and 8 processors is slightly larger than 5 for all the matrices where the speed-up for LinBox ranges from 4 to 4.9. Furthermore, the speed-up obtained with SpMV versus LinBox on the sequence generation seems scalable as it even improves when used in a parallel setting.

Conclusion

We have proposed a new SpMV library providing good results on Z /m Z rings.

To attain this efficiency it has been mandatory to augment the complexity of the SpMV algorithms, since OpenMP, Cuda et al. all manage differently the parallelization. Nonetheless, we provide new hybrid formats that improve the performance. Moreover we have also specialized it to the computation of a sequence of matrix-vector products together with a new parallelization of the As seen in 3.2.2, a first parallelization of the σ-basis computation has been achieved. Its efficiency is not matching the expected scalability and lot of work needs to be done to circumvent this problem. First, a deeper study on the parallelization of σ-basis computation has to be done. Beside the parallelization of PM-Basis and M-Basis algorithms themselves, we need to design new algorithms to avoid the numerous task dependencies, inherent to the existing methods. This will also enable an easier parallelization of early termination strategies (requiring to interleave the generation sequence and the σ-basis computation).

Another important task is to extend the sigma-basis algorithm to work on polynomial matrices over extension fields. Indeed the use of random projections U and V over extension fields might improve the probabilities to get the full minimal polynomial of the matrix [START_REF] Kaltofen | Analysis of Coppersmith's block Wiedemann algorithm for the parallel solution of sparse linear systems[END_REF][START_REF] Villard | A study of Coppersmith's block Wiedemann algorithm using matrix polynomials[END_REF][START_REF] Chen | Efficient matrix preconditioners for black box linear algebra[END_REF]. As shown in this paper and in [START_REF] Dumas | Parallel computation of the rank of large sparse matrices from algebraic k-theory[END_REF], σ-basis needs only a polynomial matrix multiplication implementation to work. In order to current LinBox's implementation to extension field, we will use the same technique as [START_REF] Dumas | Finite field linear algebra subroutines[END_REF]: first use Kronecker substitution to transform the extension field polynomial representation to an integer representation ; then use a Chinese remaindered version of the polynomial matrix multiplication to recover the resulting matrix polynomial over Z ; and finally convert back the integers using e.g. the REDQ inverse operation of [START_REF] Dumas | Q-adic transform revisited[END_REF].

The SpMV implementation also needs further work and other directions to be explored. For instance, we need to have dedicated implementations in Z /2 Z where x and y can be compressed. More formats, including dense submatrices, have yet to be explored, which is linked to spending some more time on pre-processing the matrix: for instance the use of Metis6 for partitioning and reordering A would also improve the performance. It will be interesting to deal with matrices such that A and ⊤ A cannot be simultaneously stored ( [START_REF] Buluç | Parallel sparse matrix-vector and matrix-transpose-vector multiplication using compressed sparse blocks[END_REF]). This problem indeed occurs on GPU's where on-chip memory is very limited. Finally, we will also provide multi-GPU and hybrid GPU/CPU implementations.

2 ,

 2 for (i=0 ; i<n ; ++i){ y += a[i] * b[i] ; } y = fmod(y,m); for (i=0 ; i<n ; ++i){ y += a[i] * b[i] ; y = fmod(y,m); } In this case, suppose y = 0 and a[i], b[i] are reduced modulo m at first, and M is the largest representable integer. Say that on Z /m Z, we represent the ring on 0, m -1 . Then we can do at most M/m 2 such accumulations before reducing. We can also represent the ring on -m-1

Figure 1 :

 1 Figure 1: float-double trade-off for different sizes of m, on the CPU and GPU

Figure 2 :

 2 Figure 2: Using numerical routines

!

Figure 5 :

 5 Figure 5: Matrix-multivector multiplication speed on one 3.2GHz Intel Xeon CPU (left) and a Nvidia GTX280 GPU (right) for column-major multi-vectors, with 1, 4, 8 and 16 vectors.

2 k + O( 1 )

 1 matrix entries. This gives a parallel complexity of O( n 2 d log d k ) operations in F. • step 3 requires the computation of 2d independent matrix multiplications of dimension n, which gives O(n 3 d) operations in F. One can easily see how to distribute this work on k processors such that each processor has a workload of O( n 3 d k ) operations.

Figure 8 :

 8 Figure 8: Scalability of parallel polynomial matrix multiplication with LinBox and OpenMP on a 16 core machine (based on Quad-Core AMD Opteron). n is the matrix dimension.

Figure 9 :

 9 Figure 9: Scalability of parallel σ-basis computation with LinBox and OpenMP on a 16 core machine (based on a Quad-Core AMD Opteron). n is the matrix dimension of the series.

Table 1 :

 1 presented in table 1. Matrices overview

	name mat1916 bibd_81_3 EX5	GL7d15 mpolyout2
	row 1916	3240	6545	460261	2410560
	col 1916	85320	6545	171375	2086560
	nbnz 195985	255960	295680 6080381 15707520
	rank 1916	3240	4740	132043	1352011

Table 2 :

 2 Rank modulo 65521 with OpenMP Parallel block Wiedemann on a Xeon E5345, 8 × 2.33GHz sigma-basis algorithm in order to enhance e.g. rank computations of very large sparse matrices.

	Matrix	mat1916		bibd_81_3	EX5
	Cores	1	8	1	8	1	8
	Seq-LB	15.09 3.08 47.73	12.41	84.21	20.22
	Seq-SpMV	5.02 0.91 41.28	7.56	49.66	7.36
	σ-basis	9.02 1.64 18.45	3.63	37.45	8.39
	Interpolation	0.37 0.29	1.07	0.82	2.29	1.75
	Total-LB	24.48 5.01 67.25	16.86 123.95	30.36
	Total-SpMV	14.41 2.84 60.80 12.01	89.40 17.50

https://ljkforge.imag.fr/projects/ffspmvgpu/

http://linalg.org

matrices available at http://www-ljk.imag.fr/membres/Jean-Guillaume.Dumas/simc. html

www.netlib.org/blas/blast-forum/chapter3.pdf

www.linalg.org

http://glaros.dtc.umn.edu/gkhome/metis/metis/overview

computation of the sequence A i x i∈ 0,m that are used in many of the black box methods.

On figure 6, we illustrate this differences, mostly reusing or not the data on the GPU, by comparing the performance of the following two pseudo-codes: We confirm on figure 6 that it is highly desirable to not move data on the GPU when avoidable.