T Cluzeau
email: cluzeau@ensil.unilim.fr

V Dolean
email: dolean@unice.fr

F Nataf
email: nataf@ann.jussieu.fr

A Quadrat
email: alban.quadrat@inria.fr

Symbolic methods for developing new domain decomposition algorithms

Keywords: Systems of partial differential equations, domain decomposition methods, symbolic computation, systems theory, algebraic analysis, decoupling methods, Gröbner basis techniques Systèmes d'équations aux dérivées partielles, méthodes de décomposition de domaine, calcul formel, théorie des systèmes, analyse algébrique, méthodes de découplage, techniques

The purpose of this work is to show how algebraic and symbolic techniques such as Smith normal forms and Gröbner basis techniques can be used to develop new Schwarz-like algorithms and preconditioners for linear systems of partial differential equations.

Introduction

Some algorithmic aspects of systems of PDEs based simulations can be better clarified by means of symbolic computation techniques. Numerical simulations heavily rely on solving systems of partial differential equations (PDEs), which can be linear or non-linear, time or space-dependent or stationary. Non-linear problems are solved by fixed-point methods or Newton-type algorithms, which themselves lead to solving linear systems of equations. For three (or higher) dimensional problems, the direct numerical methods can be used to solve these systems only when the number of unknowns is small. That is because they have important memory requirement and computational costs. Hence, for the large-scale problems we deal with in today's standard applications, it is necessary to rely on iterative Krylov methods that are scalable (i.e., weakly dependent of the number of degrees of freedom and of the number of subdomains) and have limited memory requirements. They are preconditioned by domain decomposition methods, incomplete factorizations, multigrid preconditioners, . . . These techniques are well understood and efficient for scalar symmetric equations (e.g., Laplacian, biLaplacian) and to some extent for non-symmetric equations (e.g., convection-diffusion). But they have poor performances and lack robustness when used for symmetric systems of PDEs, and even more so for non-symmetric complex systems (e.g., fluid mechanics, porous media). As a general rule, the study of iterative solvers for systems of PDEs as opposed to scalar PDEs is an underdeveloped subject.

We aim at building new robust and efficient solvers, such as domain decomposition methods and preconditioners for some linear and well-known systems of PDEs. In particular, we shall concentrate on Neumann-Neumann and Finite Element Tearing and Interconnecting (FETI) type algorithms which are very popular and well-known for scalar symmetric positive definite second order problems. For instance, see [START_REF] Roeck | Analysis and test of a local domain decomposition preconditioner[END_REF][START_REF] Mandel | Balancing domain decomposition[END_REF][START_REF] Mandel | Balancing domain decomposition: Theory and performance in two and three dimensions[END_REF][START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF][START_REF] Mandel | An algebraic theory for primal and dual substructuring methods by constraints[END_REF], and to some point to different other problems, like the advection-diffusion equations [START_REF] Gerardo-Giorda | A Robin-Robin preconditioner for advectiondiffusion equations with discontinuous coefficients[END_REF][START_REF] Achdou | A domain decomposition preconditioner for an advection-diffusion problem[END_REF], plate and shell problems [START_REF] Tallec | A Neumann-Neumann domain decomposition algorithm for solving plate and shell problems[END_REF], or the Stokes equations [START_REF] Pavarino | Balancing Neumann-Neumann methods for incompressible Stokes equations[END_REF]35]. This work is motivated by the fact that in some sense these methods applied to systems of PDEs (such as Stokes, Oseen, linear elasticity) are less optimal than the domain decomposition methods for scalar problems. Indeed, in the case of two subdomains consisting of the two half planes, it is well-known that, for scalar equations like the Laplace problem, the Neumann-Neumann preconditioner is an exact preconditioner for the Schur complement equation. We recall that a preconditioner is called exact if the preconditioned operator simplies to the identity. Unfortunately, this does not hold in the vector case.

In order to achieve this goal, we use algebraic methods developed in constructive algebra, D-modules (differential modules) theory, and symbolic computation such as the so-called Smith or Jacobson normal forms and Gröbner basis techniques for transforming a linear system of PDEs into a set of independent scalar PDEs. Decoupling linear partial differential (PD) systems leads to the design of new numerical methods based on the efficient techniques dedicated to scalar PDEs (e.g., Laplace equation, advection-diffusion equation). Moreover, these algebraic and symbolic methods provide important intrinsic information about the linear system of PDEs to solve which need to be taken into account in the design of new numerical methods which can supersede the usual ones based on a direct extension of the classical scalar methods to the linear systems.

This algebraic approach enables an intrinsic analysis of linear systems of PDEs. In particular, we are interested in transforming the linear system of PDEs into a set of decoupled PDEs under certain types of invertible transformations. The problem of decoupling the equations of linear time-varying systems of ordinary differential equations (ODEs) or difference equations has been an important issue in the symbolic computation community (e.g., methods based on the so-called eigenring). For instance, the techniques based on the eigenring can be considered as a general-Inria ization for time-varying linear systems of ODEs of the classical diagonalization method used for solving time-invariant first order linear system of ODEs with constant coefficients. Nowadays, different algorithms and implementations in computer algebra systems use these techniques for integrating time-varying linear systems of ODEs in closed forms. These results recently have been extended in [START_REF] Cluzeau | Factoring and decomposing a class of linear functional systems[END_REF][START_REF] Cluzeau | On algebraic simplifications of linear functional systems[END_REF] for linear systems of PDEs based on D-modules and Gröbner bases techniques and implemented in the OreMorphisms package [START_REF] Cluzeau | OreMorphisms: a homological algebraic package for factoring, reducing and decomposing linear functional systems[END_REF]. The invertible transformations used in these methods are defined by unimodular matrices, namely, invertible matrices over the (non-commutative) ring of PD operators.

An alternative way for decoupling a linear system of ODEs is to use the so-called Smith normal form of the matrix of OD operators associated with the linear system. This normal form was introduced by H. J. S. Smith (1826-1883) for matrices with integer entries (see, e.g., [START_REF] Kailath | Linear Systems[END_REF] or Theorem 1.4 of [START_REF] Wloka | Boundary Value Problems for Elliptic Systems[END_REF]). The Smith normal form has already been successfully applied to open problems in the design of Perfectly Matched Layers (PML). The theory of PML for scalar equations was well-developed and the usage of the Smith normal form allowed one to extend these works to linear systems of PDEs. In [START_REF] Nataf | A new approach to perfectly matched layers for the linearized Euler system[END_REF], a general approach is proposed and applied to the particular case of the compressible Euler equations that model aero-acoustic phenomena.

For domain decomposition methods, several preliminary results have been obtained for the compressible Euler equations [START_REF] Dolean | A new domain decomposition method for the compressible Euler equations[END_REF], Stokes and Oseen systems [START_REF] Dolean | Deriving a new domain decomposition method for the Stokes equations using the Smith factorization[END_REF]. More precisely, it has been shown that the Stokes system is equivalent to a biLaplacian in 2D and to a Laplacian and a biLaplacian in 3D. Hence, the classical Neumann-Neumann algorithm for the Stokes system [35] can be recast into an algorithm for biLaplacian and Laplacian. The resulting algorithms are different from the classical algorithms for these scalar equations and are therefore not optimal. To fix this problem, we first used optimal algorithms for the biLaplacian and the Laplacian, and then via the Smith factorization, we back-transform them formally using only symbolic computations (derivation, linear combination of different equations) into optimal algorithms for the Stokes system. Interestingly enough, the classical algorithm for the Stokes system is based on solving either stress imposed or displacement imposed problems in the subdomains. In the new algorithm developed in [START_REF] Dolean | Deriving a new domain decomposition method for the Stokes equations using the Smith factorization[END_REF], the tangential stress and the normal displacement problems are solved in each subdomain. Preliminary tests show that the new algorithm can be twice as fast as the classical one. Moreover, when solving time-dependent problems, the new algorithm has an iteration count independent of the time step. This property is not true for the classical algorithm. The approach is quite general and has already been applied to the Oseen (linearized Navier-Stokes) equations and to the compressible Euler equations. For both equations, the new algorithms are the first genuine extensions of the Neumann-Neumann or FETI algorithms. This report is organized as follows: Section 2 recalls the notion of the Smith normal form of a matrix over a principal ideal domain and states the problem we are considering in what follows. We also introduce the examples that we are using to illustrate our results in the whole paper, namely the linear elasticity and Oseen/Stokes systems. In Section 3, we give optimal algorithms for the scalar operators that appears in the study of the latter linear elasticity or Oseen/Stokes systems. Section 4 states one of the main problem studied in this paper: the choice of Smith variables. We recall how this was done by hand calculations in previous studies and ask whether symbolic computation techniques can give an automatic way for constructing many possible Smith variables. Then, in Section 5, we recall useful notion of module theory that are needed to understand the sequel. We show that the problem of choosing relevant Smith variables can be reduced to a completion problem and give an algorithm for computing many relevant Smith variables. Section 6 shows how symbolic computation techniques (e.g., Gröbner or Janet basis computations) can also be useful to automatically reduce the interface conditions that appear in the algorithm developed in Section 3, and thus to speed up the computations.

In the first appendix (Section 7), we recall the important concept of Gröbner basis for ideals and modules over commutative polynomial rings, which is the main computational tool used in the study of linear systems of PDEs via module theory. Finally, Appendix 8 illustrates the Maple implementation of our algorithms in the Schwarz package which computes relevant Smith variables and reduces interface conditions. Finally, we illustrate this on the examples of the linear elasticity (Cauchy-Navier equations) and Oseen/Stokes systems.

Notations. If R is a ring, then R p×q is the set of p × q matrices with entries in R. Moreover, the general linear group of index p of R, denoted by GL p (R), is the group of invertible matrices of R p×p , namely:

GL p (R) = {E ∈ R p×p | ∃ F ∈ R p×p : E F = F E = I p }.
An element of GL p (R) is called an unimodular matrix. A diagonal matrix with elements d i 's on the diagonal will simply be denoted by diag(d 1 , . . . , d p ). If k is a field (e.g., k = Q, R, C), then k[x] is the commutative polynomial ring in x with coefficients in k. In particular, an element r ∈ k[x] has the form r = d i=0 a i x i , where the a i 's belong to k and d ∈ N = {0, 1, . . . , }. More generally, if k is a field, then k[x 1 , . . . , x n ] is the commutative ring of polynomials in x 1 , . . . , x n with coefficients in k. In what follows, k(x 1 , . . . , x n ) will denote the field of rational functions in x 1 , . . . , x n with coefficients in k. Finally, if r, r ′ ∈ R, then the notation r ′ | r means that r ′ divides r, i.e., that there exists r ′′ ∈ R such that r = r ′′ r ′ .

Smith normal forms of linear systems of PDEs

We first introduce the concept of Smith normal form of a matrix with polynomial entries (see, e.g., [START_REF] Kailath | Linear Systems[END_REF], [START_REF] Wloka | Boundary Value Problems for Elliptic Systems[END_REF], Theorem 1.4). This concept will play an important role in what follows. In particular, we can take

d i = m i /m i-1
, where m i is the greatest common divisor of all the i × i-minors of A (i.e., the determinants of all i × i-submatrices of A), with the convention that m 0 = 1.

The matrix D = diag(d 1 , . . . , d p ) ∈ R p×p is called a Smith normal form of A and the diagonal elements d i 's are the invariants factors of A.

Remark 1. We note that for R = k[s], E ∈ GL p (R) is equivalent to det(E) is an invertible polynomial, i.e., a nonzero element of k. The invariants factors of A are uniquely defined by A up to multiplication by nonzero elements of k. In particular, in what follows, we shall assume that the d i 's are monic polynomials, i.e., their leading coefficients are 1, which will allow us to call the matrix D = diag(d 1 , . . . , d p ) the Smith normal form of A. Finally, the relation A = E D F shows that the matrix E (resp., F ) operates on the columns (resp., rows) of D.

The Smith normal form is a mathematical concept which is classically used in the literature of symbolic computation, ordinary differential systems and control theory (see, e.g., [START_REF] Kailath | Linear Systems[END_REF]). Since its use in scientific computing is rather new, we give here a few comments:

Inria

• Smith was an English mathematician of the end of the 19 th century. He worked in number theory and considered the problem of factorizing matrices with integer entries. We gave here the polynomial version of his theorem in the special case where the matrix A is square but the result is more general and applies as well when the matrix A is rectangular with entries in a principal ideal domain (e.g., Euclidean domain such as

Z, k[x],
where k is a field). For more details, see [START_REF] Cox | Ideals, Varieties, and Algorithms[END_REF][START_REF] Eisenbud | Commutative Algebra with a View Toward Algebraic Geometry[END_REF][START_REF] Kailath | Linear Systems[END_REF].

• As stated in Remark 1, the invariants factors d i 's of A are uniquely defined by A up to nonzero elements of k. Hence, if the d i 's are chosen to be monic polynomials, then the Smith normal form D = diag(d 1 , . . . , d p ) of A is unique. But, this result does not hold for the matrices E and F , i.e., they are not uniquely defined by A.

• One of the main interests of Theorem 1 is the following. If det(A) = 0, then by Cramer's formula, A -1 ∈ k(s) p×p , i.e., the inverse of A is generally a matrix with rational function entries. Since det(E) and det(F ) are nonzero elements of the field k, the inverse of E and F are matrices with polynomial entries in s, i.e., E -1 , F -1 ∈ R p×p . Using the Smith normal form of A, we get

A -1 = F -1 D -1 E -1
, which shows that the rational part of the inverse of A only appears in the intrinsic diagonal matrix D -1 .

• The proof of Theorem 1 is constructive and gives an algorithm for computing the matrix D and matrices E and F such that A = E D F . The computation of Smith normal forms is available in many computer algebra systems such as Maple, Mathematica, Magma, . . .

Consider now the following model problem in R d with d = 2, 3:

L d (w) = g in R d , (1) |w
(x)| → 0 for |x| → ∞.
(2)

For instance, L d (w) = 0 can represent the Stokes/Oseen/linear elasticity equations in dimension d. Moreover, if we suppose that the inhomogeneous linear system of PDEs (1) has constant coefficients, then (1) can be rewritten as

A d w = g, (3)
where In what follows, we shall study the domain decomposition problem in which R d is divided into subdomains and the case A d ∈ R p×p . The direction normal to the interface of the subdomains is particularized and denoted by ∂

A d ∈ R q×p , R = k[∂ x , ∂ y ] if d = 2, or R = k[∂ x , ∂ y , ∂ z ] if d = 3,
x . If S = k(∂ y )[∂ x ] for d = 2 or S = k(∂ y , ∂ z )[∂ x ] for d = 3, then computing the Smith normal form of the matrix A d ∈ S p×p , we get A d = E D F , where D ∈ S p×p , E ∈ GL p (S)
and F ∈ GL p (S). In particular, the entries of the matrices E, D, F are polynomials in ∂ x , and E and F are unimodular matrices, i.e., det(E), det(

F ) ∈ k(∂ y ) \ {0} if d = 2, or det(E), det(F ) ∈ k(∂ y , ∂ z ) \ {0} if d = 3.
We recall that the matrices E and F are not unique contrary to D. Then, using the Smith normal form of A d , we have:

A d w = g ⇔ w s := F w, D w s = E -1 g. (4)
In other words, (3) is equivalent to the following uncoupled linear PD system:

D w s = E -1 g. (5)
Since E ∈ GL p (S) and F ∈ GL p (S), the entries of their inverses are still polynomial in ∂ x . Thus, applying E -1 to the right-hand side g amounts to taking linear combinations of derivatives of g with respect to x.

If R d is split into two subdomains R -× R d-1 and R + × R d-1 , where R -:= {x ∈ R | x < 0} and R + := {x ∈ R | x > 0}
, then the application of E -1 and F -1 to vectors can be done for each subdomain independently. No communication between the subdomains is necessary.

In conclusion, it is enough to find a domain decomposition algorithm for the uncoupled system [START_REF] Chyzak | Effective algorithms for parametrizing linear control systems over Ore algebras[END_REF] and then transform it back to the original one (3) by means of the matrix F . This technique can be applied to any linear system of PDEs once it is rewritten in a polynomial form. Moreover, the uncoupled system acts on the variables w s , which we shall further call Smith variables since they are issued from the Smith normal form.

Remark 2. Since the matrix F is used to transform (5) to (3) (see the first equation of the right-hand side of ( 4)) and the matrix F is not unique, we need to find a matrix F as simple as possible (e.g., F has minimal degree in ∂ x ) so that to obtain a final algorithm that can be used for practical computations.

Application to Cauchy-Navier equations

Consider the two or three dimensional elasticity operator

E d (u) := -µ ∆ u -(λ + µ) ∇ div u,
where λ, µ are the two Lamé constants.

We first study the two-dimensional case. Let us consider the commutative polynomial rings

R = Q(λ, µ)[∂ x , ∂ y ] and S = Q(λ, µ)(∂ y )[∂ x ] = Q(λ, µ, ∂ y )[∂ x ], ∆ = ∂ 2 x + ∂ 2 y , and
A 2 = (λ + 2 µ) ∂ 2 x + µ ∂ 2 y (λ + µ) ∂ x ∂ y (λ + µ) ∂ x ∂ y µ ∂ 2 x + (λ + 2 µ) ∂ 2 y ∈ R 2×2
the matrix of PD operators associated with E 2 , i.e., E 2 (u) = A 2 u. Then, the Smith normal form of A 2 ∈ S 2×2 is defined by:

D A2 = 1 0 0 ∆ 2 . ( 6
)
In particular, the matrices E, F ∈ GL 2 (S) can be chosen as follows:

E = (λ + µ) ∂ x ∂ y µ ∂ -2 y µ ∂ x + (λ + 2 µ) ∂ 2 y µ 2 (λ + µ) -1 ∂ x ∂ -3 y , F = (λ + µ) -1 (λ -µ) ∂ -3 y ∂ x 1 1 0 .
The particular form of D A2 shows that the system of equations for the linear elasticity in R 2 is algebraically equivalent over S to a biharmonic equation ∆ 2 v = 0. We now study the three-dimensional case. Let us consider the commutative polynomial rings

R = Q(λ, µ)[∂ x , ∂ y , ∂ z ] and S = Q(λ, µ)(∂ y , ∂ z )[∂ x ] = Q(λ, µ, ∂ y , ∂ z )[∂ x ], ∆ = ∂ 2 x + ∂ 2 y + ∂ 2 z , and
A 3 =    -(λ + µ) ∂ 2 x -µ ∆ -(λ + µ) ∂ x ∂ y -(λ + µ) ∂ x ∂ z -(λ + µ) ∂ x ∂ y -(λ + µ) ∂ 2 y -µ ∆ -(λ + µ) ∂ y ∂ z -(λ + µ) ∂ x ∂ z -(λ + µ) ∂ y ∂ z -(λ + µ) ∂ 2 z -µ ∆    ∈ R 3×3 (7)
Inria the matrix of PD operators associated with E 3 , i.e., E 3 (u) = A 3 u. Then, the Smith normal form of A 3 ∈ S 3×3 is:

D A3 =    1 0 0 0 ∆ 0 0 0 ∆ 2    . (8)
The particular form of D A3 shows that the system of equations for the linear elasticity in R 3 is algebraically equivalent over S to the uncoupled system formed by the Laplace equation ∆ v = 0 and the biharmonic equation ∆ 2 w = 0.

Application to Oseen and Stokes equations

Consider the two or three dimensional Oseen operator

O d (w) = O d (v, q) := (c v -ν ∆v + b • ∇v + ∇q, ∇ • v),
where b = (b j ) 1≤j≤d is the convection velocity, c the reaction coefficient, and ν the viscosity. If b = 0, then we obtain the Stokes operator

S d (w) = S d (v, q) := (c v -ν ∆v + ∇q, ∇ • v).
We start with the two-dimensional case. The spatial variables are denoted by x and y.

If we consider R = Q(b 1 , b 2 , c, ν)[∂ x , ∂ y ] and S = Q(b 1 , b 2 , c, ν)(∂ y )[∂ x ] = Q(b 1 , b 2 , c, ν, ∂ y )[∂ x ], ∆ = ∂ 2 x + ∂ 2 y , and
O 2 =    -ν ∆ + b 1 ∂ x + b 2 ∂ y + c 0 ∂ x 0 -ν ∆ + b 1 ∂ x + b 2 ∂ y + c ∂ y ∂ x ∂ y 0   
the matrix of PD operators associated with O 2 , i.e., O 2 (w) = O 2 w, then the Smith normal form of O 2 ∈ S 3×3 is defined by

D O2 =    1 0 0 0 1 0 0 0 ∆ L 2    , (9)
where

L 2 = c -ν ∆ + b 1 ∂ x + b 2 ∂ y .
The matrices E, F ∈ GL 3 (S) can be chosen as follows:

E =    ∂ x -∂ -1 y ∂ x L 2 -ν ∂ -2 y ∂ y 0 0 0 ∂ y 0    , F =    0 -∂ -1 y L 2 1 ∂ -1 y ∂ x 1 0 1 0 0    .
From the form of D O2 we can deduce that the two-dimensional Oseen equations can be mainly characterized by the scalar fourth order PD operator ∆ L 2 = ∆ (cν ∆ + b • ∇). This is not surprising since the stream function formulation of the Oseen problem gives the same PDE for the stream function in the two-dimensional case (see, e.g., [19]).

In the three-dimensional case, the spatial variables are now denoted by x, y and z.

Let us consider R = Q(b 1 , b 2 , c, ν)[∂ x , ∂ y , ∂ z ], S = Q(b 1 , b 2 , c, ν)(∂ y , ∂ z )[∂ x ] = Q(b 1 , b 2 , c, ν, ∂ y , ∂ z )[∂ x ], ∆ = ∂ 2 x + ∂ 2 y + ∂ 2 z , and
O 3 =     -ν ∆ + b1 ∂x + b2 ∂y + c 0 0 ∂x 0 -ν ∆ + b1 ∂x + b2 ∂y + c 0 ∂y 0 0 -ν ∆ + b1 ∂x + b2 ∂y + c ∂z ∂x ∂y ∂z 0     (10)
the matrix of PD operators defining O 3 (v, q), then the Smith normal form of O 3 is

D O3 =      1 0 0 0 0 1 0 0 0 0 L 3 0 0 0 0 ∆ L 3      , (11)
where:

L 3 = c -ν ∆ + b 1 ∂ x + b 2 ∂ y + b 3 ∂ z = c -ν ∆ + b • ∇. ( 12
)
Analogously to the two-dimensional case, we see that the Oseen operator O 3 is determined by the diagonal matrix D O3 . Therefore, it can be represented by the fourth order PD operator ∆ L 3 and the second order PD operator L 3 .

Remark 3. The above applications of the Smith normal forms suggest the following conclusion: one should design an optimal domain decomposition method for

• the biharmonic operator (resp., L 2 ∆) in the case of linear elasticity (resp., for the Oseen/Stokes equations) for the two-dimensional case,

• the system formed by the Laplace and biharmonic operators (resp., L 3 and L 3 ∆) in the case of linear elasticity (resp., for the Oseen/Stokes equations) in the three-dimensional case, and then back-transform it to the original system.

Optimal domain decomposition algorithms for scalar equations

The Neumann-Neumann or FETI methods are well-known for some symmetric scalar equations such as Laplace equations (see [START_REF] Roeck | Analysis and test of a local domain decomposition preconditioner[END_REF][START_REF] Mandel | Balancing domain decomposition[END_REF][START_REF] Mandel | Balancing domain decomposition: Theory and performance in two and three dimensions[END_REF][START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF][START_REF] Mandel | An algebraic theory for primal and dual substructuring methods by constraints[END_REF]). We can give an example of such a method in its iterative version. Consider a decomposition of the domain Ω = R 2 into two half planes

Ω 1 = R -× R and Ω 2 = R + × R.
Let the interface {0} × R be denoted by Γ and (n i ) i=1,2 the outward normal of (Ω i ) i=1,2 . We consider the following problem: Find u : R 2 → R such that

-∆u = g in R 2 , |u(x)| → 0 for |x| → ∞, ( 13
)
where g is a given right-hand side.

The following algorithm is optimal in the sense that converges in two iterations.

Algorithm 1. Let u n i be the local solution in the domain Ω i at iteration n. We choose the initial values u 0 1 and u 0 2 such that u 0 1 = u 0 2 on Γ. We obtain (u n+1 i ) i=1,2 from (u n i ) i=1,2 by the following iterative procedure: Correction step. We compute the corrections ũn+1 i i=1,2 as solutions of the following homogeneous local problems

Inria            -∆ ũn+1 1 = 0 in Ω 1 , lim |x|→∞ |ũ n+1 1 | = 0, ∂ ũn+1 1 ∂ n 1 = γ n on Γ,            -∆ ũn+1 2 = 0 in Ω 2 , lim |x|→∞ |ũ n+1 2 | = 0, ∂ ũn+1 2 ∂ n 2 = γ n on Γ, ( 14
)
where

γ n = - 1 2 ∂ u n 1 ∂ n 1 + ∂ u n 2 ∂ n 2 .
Updating step. We update (u n+1 i ) i=1,2 by solving the following local problems

       -∆ u n+1 1 = g in Ω 1 , lim |x|→∞ |u n+1 1 | = 0, u n+1 1 = u n 1 + δ n+1 on Γ,        -∆ u n+1 2 = g in Ω 2 , lim |x|→∞ |u n+1 2 | = 0, u n+1 2 = u n 2 + δ n+1 on Γ, ( 15
)
where

δ n+1 = 1 2 ũn+1 1 + ũn+1 2 .
Since the biharmonic operator (or its very similar form ∆ (cν ∆)) seems to play a key role in the design of a new algorithm for both Stokes and elasticity problem, inspired by the Neumann-Neumann algorithm for the Laplace equation (see Algorithm 1), we need to build an optimal algorithm (converging in two iterations) for it.

For the sake of simplicity, we only illustrate the example of the ∆ 2 operator, the conclusions for ∆ (cν ∆) being identical. We consider the following problem: Find φ : R 2 → R such that

∆ 2 φ = g in R 2 , |φ(x)| → 0 for |x| → ∞, ( 16
)
where g is a given right-hand side.

The domain Ω is decomposed into two half planes Ω 1 = R -× R and Ω 2 = R + × R. Let the interface {0} × R be denoted by Γ and (n i ) i=1,2 the outward normal of (Ω i ) i=1,2 . The algorithm we propose is given as follows:

Algorithm 2. We choose the initial values φ 0 1 and φ 0 2 such that φ

0 1 = φ 0 2 , ∆ φ 0 1 = ∆ φ 0 2 , on Γ. We obtain (φ n+1 i ) i=1,2 from (φ n i ) i=1
,2 by the following iterative procedure:

Correction step. We compute the corrections

φn+1 i i=1,2
as solutions of the following homogeneous local problems

                     ∆ 2 φn+1 1 = 0 in Ω 1 , lim |x|→∞ | φn+1 1 | = 0, ∂ φn+1 1 ∂ n 1 = γ n 1 on Γ, ∂ ∆ φn+1 1 ∂ n 1 = γ n 2 on Γ,                      ∆ 2 φn+1 2 = 0 in Ω 2 , lim |x|→∞ | φn+1 2 | = 0, ∂ φn+1 2 ∂ n 2 = γ n 1 on Γ, ∂ ∆ φn+1 2 ∂ n 2 = γ n 2 on Γ, (17)

RR n°7953

where

γ n 1 = - 1 2 ∂ φ n 1 ∂ n 1 + ∂ φ n 2 ∂ n 2 and γ n 2 = - 1 2 ∂ ∆ φ n 1 ∂ n 1 + ∂ ∆ φ n 2 ∂ n 2 .
Updating step. We update (φ n+1 i ) i=1,2 by solving the following local problems

           ∆ 2 φ n+1 1 = g in Ω 1 , lim |x|→∞ |φ n+1 1 | = 0, φ n+1 1 = φ n 1 + δ n+1 1 on Γ, ∆ φ n+1 1 = ∆ φ n 1 + δ n+1 2 on Γ,            ∆ 2 φ n+1 2 = g in Ω 2 , lim |x|→∞ |φ n+1 2 | = 0, φ n+1 2 = φ n 2 + δ n+1 1 on Γ, ∆ φ n+1 2 = ∆ φ n 2 + δ n+1 2 on Γ, ( 18
)
where δ n+1

1 = 1 2 φn+1 1 + φn+1 2 and δ n+1 2 = 1 2 ∆ φn+1 1 + ∆ φn+1 2 .
This algorithm has a remarkable property: it is optimal in the sense that converges in two iterations. The proof of its optimality can be easily done by Fourier transform techniques (see [START_REF] Dolean | Deriving a new domain decomposition method for the Stokes equations using the Smith factorization[END_REF]). This is a generalization of the Neumann-Neumann algorithm for the ∆ operator acting on a Smith variable. For instance, in the case of the two dimensional linear elasticity, φ represents the second component of the vector of Smith variables, that is, φ = (w s ) 2 = (F u) 2 , where u = (u, v) is the displacement field. We now need to replace φ with (F u) 2 into Algorithm 2 and then simplify it using algebraically admissible operations. Thus, one can obtain an optimal algorithm for the Stokes (resp., Cauchy-Navier) equations. But this depends on the form of F .

In the three dimensional case, we should rather consider the Laplace/biharmonic system. We thus need to study the following problem: Find φ : R 3 → R and ψ : R 3 → R such that

∆(ψ) = h in R 3 , |ψ(x)| → 0 for |x| → ∞, ∆ 2 (φ) = g in R 3 , |φ(x)| → 0 for |x| → ∞, (19)
where g and h are given right-hand sides. The domain Ω is decomposed into two half spaces

Ω 1 = R -× R 2 and Ω 2 = R + × R 2 .
Let the interface {0} × R 2 be denoted by Γ and (n i ) i=1,2 the outward normal of (Ω i ) i=1,2 . We have the following optimal algorithm:

Algorithm 3. We choose the initial values φ 0 1 , φ 0 2 , ψ 0 1 , ψ 0 2 such that φ 0 1 = φ 0 2 , ψ 0 1 = ψ 0 2 and ∆ φ 0 1 = ∆ φ 0 2 on Γ. We obtain (φ n+1 i ) i=1,
                                         ∆ ψn+1 1 = 0 in Ω 1 , ∆ 2 φn+1 1 = 0 in Ω 1 , lim |x|→∞ | ψn+1 1 | = 0, lim |x|→∞ | φn+1 1 | = 0, ∂ ψn+1 1 ∂ n 1 = δ n on Γ, ∂ φn+1 1 ∂ n 1 = γ n 1 on Γ, ∂ ∆ φn+1 1 ∂ n 1 = γ n 2 on Γ,                                          ∆ ψn+1 2 = 0 in Ω 2 , ∆ 2 φn+1 2 = 0 in Ω 2 , lim |x|→∞ | ψn+1 2 | = 0, lim |x|→∞ | φn+1 2 | = 0, ∂ ψn+1 2 ∂ n 2 = γ n on Γ, ∂ φn+1 2 ∂ n 2 = γ n 1 on Γ, ∂ ∆ φn+1 1 ∂ n 1 = γ n 2 on Γ, (20)
where:

γ n = - 1 2 ∂ ψ n 1 ∂ n 1 + ∂ ψ n 2 ∂ n 2 , γ n 1 = - 1 2 ∂ φ n 1 ∂ n 1 + ∂ φ n 2 ∂ n 2 , γ n 2 = - 1 2 ∂ ∆ φ n 1 ∂ n 1 + ∂ ∆ φ n 2 ∂ n 2 .
Updating step. We update (φ n+1 i ) i=1,2 and (ψ n+1 i ) i=1,2 by solving the following local problems:

                             ∆ ψ n+1 1 = h in Ω 1 , ∆ 2 φ n+1 1 = g in Ω 1 , lim |x|→∞ |ψ n+1 1 | = 0, lim |x|→∞ |φ n+1 1 | = 0, ψ n+1 1 = ψ n 1 + δ n+1 on Γ, φ n+1 1 = φ n 1 + δ n+1 1 on Γ, ∆ φ n+1 1 = ∆φ n 1 + δ n+1 2 on Γ,                              ∆ ψ n+1 2 = h in Ω 2 , ∆ 2 φ n+1 2 = g in Ω 2 , lim |x|→∞ |ψ n+1 2 | = 0, lim |x|→∞ |φ n+1 2 | = 0, ψ n+1 2 = ψ n 2 + δ n+1 on Γ, φ n+1 2 = φ n 2 + δ n+1 1 on Γ, ∆ φ n+1 2 = ∆ φ n 2 + δ n+1 2 on Γ, (21) where:
δ n+1 = 1 2 ψn+1 1 + ψn+1 2 , δ n+1 1 = 1 2 φn+1 1 + φn+1 2 , δ n+1 2 = 1 2 ∆ φn+1 1 + ∆ φn+1 2 .
In this case, (φ, ψ) = (w s ) 2,3 = (F u) 2,3 , where u = (u, v, w) is the displacement field in R 3 and (V ) 2,3 denotes the second and third entries of a vector V .

The same question arises both for the two and three dimensional cases: Is there a simple optimal domain decomposition algorithm for the Stokes (resp., Cauchy-Navier) equations? If so, compute the corresponding F . Therefore, we need to properly choose a matrix F from the Smith normal form in order to obtain a "good" algorithm for those systems of PDEs based on the optimal ones for the biharmonic operator and for the Laplace/biharmonic operator. In [START_REF] Dolean | A new domain decomposition method for the compressible Euler equations[END_REF][START_REF] Dolean | Deriving a new domain decomposition method for the Stokes equations using the Smith factorization[END_REF], that is for Euler and Stokes/Oseen equations, the computation of the Smith normal form was done by hand or using the computer algebra system Maple. Surprisingly, the computed F has provided a good algorithm for the Stokes system even if the approach was entirely heuristic.

An approach by hand calculations

In this section, we shall show how to choose F in the case of the two dimensional elasticity problem even if the construction is not yet entirely automatic. Very similar reasoning can be carried out for the Stokes systems and for the three-dimensional problems.

To begin with, we first show some properties of the system which will guide us in the choice of the variables. Lemma 1. Consider the two-dimensional linear elasticity system:

E 2 u := -µ ∆ u -(λ + µ) ∂ x div(u) = f u -µ ∆ v -(λ + µ) ∂ y div(u) = f v , u = (u, v). ( 22
)
1. The quantities ∂ y u -∂ x v and div(u) are "annihilated" by the operator ∆, i.e., they verify an equation of the type -∆ w = g, where g is a linear PD combination of f u and f v . The same property remains true for any linear PD combination of ∂ y u -∂ x v and div(u).

2. The variables u and v are "annihilated" by the operator ∆ 2 , as well as any linear PD combination of them.

The biharmonic operator annihilates the components of u (i.e., the displacements for the linear elasticity). Furthermore, as part of the Smith normal form, it acts on a Smith variable. A natural question is then to find an F such that the second component of the corresponding Smith variables F u is exactly u or v. This leads to a simple and natural form of the algorithm.

We distinguish here two cases. In both cases, one needs to reduce the boundary conditions given in [START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF] and [START_REF] Gerardo-Giorda | A Robin-Robin preconditioner for advectiondiffusion equations with discontinuous coefficients[END_REF] by rewritting them in terms of physical variables in order to obtain well-posed local problems. The authorized operations are:

• linear combinations between the interface equations (using only PD operators in the ydirection),

• the equations inside the domain (in this case all the polynomial operations with respect to ∂ x are authorized).

In the following, we introduce the stress tensor in order to write the algorithm in a more intrinsic form:

σ(u) =   (2µ + λ) ∂u ∂x + λ ∂v ∂y µ ∂u ∂y + ∂v ∂x µ ∂v ∂x + ∂u ∂y (2 µ + λ) ∂v ∂y + λ ∂u ∂x   , u = (u, v). (23)
For any vector u, its normal (resp., tangential) component on the interface is

u ni = u • n i (resp., u τ i = u • τ i ). We denote by σ i ni := (σ i (u i ) • n i ) • n i (resp., σ i τ i := (σ i (u i ) • n i ) • τ i ) the normal (resp., tangential) part of the normal stress tensor σ i (u i ) • n i .
Case 1. u is one of the Smith variables. Using the fact that F must be polynomial in ∂ x and unimodular, one can easily obtain by direct computations that:

F = -(λ + µ) -1 ∂ -3 y ∂ x (µ ∂ 2 x -λ ∂ 2 y ) 1 1 0 .

Inria

Correction step. One needs to rewrite the interface conditions given in [START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF]. The first one is equivalent to ∂ x u which cannot be simplified. The second interface condition can be reduced by using [START_REF] Janet | Leçons sur les systèmes d'équations aux dérivées partielles[END_REF] as follows:

∂ x (∆ u) = ∂ 3 x u + ∂ x ∂ 2 y u, = -µ λ+2 µ ∂ x ∂ 2 y u -λ+µ λ+2 µ ∂ 2 x ∂ y v + ∂ x ∂ 2 y u, = λ+µ λ+2 µ ∂ y (∂ x ∂ y u -∂ 2 x v), = λ+µ λ+2 µ ∂ y λ+2 µ µ ∂ y (∂ y v + ∂ x u) , = λ+µ µ ∂ 2 y (∂ y v + ∂ x u). (24)
Integrating twice along the interface (which means removing ∂ 2 y in the last equality of ( 24)), we get ∂ y v +∂ x u. Taking a linear combination with the first interface condition ∂ x u and integrating again along the interface, we finally obtain that the second interface condition can be reduced to v. Therefore, the two quantities to be imposed along the interface can be chosen as v and σ 11 (u). Since, in this particular case, the normal vector at the interface is n = (1, 0), this is equivalent to imposing u τ = u • τ and

σ n = (σ(u) • n) • n.
Update step. The first interface condition given by ( 18) is equivalent to imposing u = u • n, and the second one can be further simplified using formally the first equation of [START_REF] Janet | Leçons sur les systèmes d'équations aux dérivées partielles[END_REF] as follows:

∆u = ∂ 2 x u + ∂ 2 y u = -µ λ+2 µ ∂ 2 y u -λ+µ λ+2 µ + ∂ 2 y u = λ+µ λ+2 µ ∂ y (∂ y u -∂ x v). (25)
Proceeding as in the Correction step above, the two interface conditions are equivalent to imposing u n = u • n and σ τ = (σ(u)

• n) • τ .
Case 2. v is one of the Smith variables. Using the fact that F must be polynomial in ∂ x and unimodular, one can easily obtain by direct computations that:

F = 1 -(λ + µ) ∂ -3 y ∂ x ((3 µ + 2 λ) ∂ 2 y + (2 µ + λ) ∂ 2 x ) 0 1 .
Performing similar calculations as in Case 1, we obtain the following conclusion: for the Correction step, the two interface conditions are equivalent to imposing u n = u•n and σ τ = (σ(u)•n)•τ .

For the Update step, the two interface conditions are equivalent to imposing u τ = u • τ and

σ n = (σ(u) • n) • n.
In conclusion, depending on the choice of F , one obtains two different algorithms for the two dimensional linear elasticity system. Algorithm 4. (Case 1) Starting with an initial guess (u 0 i ) N i=0 satisfying u 0 i,ni = u 0 j,nj and σ i τ i (u 0 i ) = σ j τ j (u 0 j ) on Γ ij , ∀ i, j, i = j, the correction step is expressed as follows:

1 ≤ i ≤ N,            E 2 (ũ n+1 i ) = 0 in Ω i , ũn+1 i,τ i = - 1 2 (u n i,τ i + u n j,τ j ) on Γ ij , σ i ni (ũ n+1 i ) = - 1 2 (σ i ni (u n i ) + σ j nj (u n j )) on Γ ij , (26)
followed by an updating step:

1 ≤ i ≤ N,                  E 2 (u n+1 i ) = g in Ω i , u n+1 i,ni = u n i,ni + 1 2 (ũ n+1 i,ni + ũn+1 j,nj ) on Γ ij , σ i τ i (u n+1 i ) = σ i τ i (u n i ) + 1 2 (σ i τ i (ũ n+1 i ) + σ j τ j (ũ n+1 j )) on Γ ij . (27)
The boundary conditions in the update step involve the normal velocity and the tangential stress, whereas in the correction step the tangential velocity and the normal stress are involved.

Algorithm 5. (Case 2) Starting with an initial guess (u 0 i ) N i=0 satisfying u 0 i,τ i = u 0 j,τ j and σ i ni (u 0 i ) = σ j nj (u 0 j ) on Γ ij , ∀ i, j, i = j, the correction step is expressed as follows:

1 ≤ i ≤ N,            E 2 (ũ n+1 i ) = 0 in Ω i , ũn+1 i,ni = - 1 2 (u n i,ni + u n j,nj ) on Γ ij , σ i τ i (ũ n+1 i ) = - 1 2 (σ i τ i (u n i ) + σ j τ j (u n j )) on Γ ij , (28)
followed by an updating step:

1 ≤ i ≤ N,                  E 2 (u n+1 i ) = g in Ω i u n+1 i,τ i = u n i,τ i + 1 2 (ũ n+1 i,τ i + ũn+1 j,τ j ) on Γ ij , σ i ni (u n+1 i ) = σ i ni (u n i ) + 1 2 (σ i ni (ũ n+1 i ) + σ j nj (ũ n+1 j )) on Γ ij . (29)
The boundary conditions in the correction step involve the normal velocity and the tangential stress, whereas in the updating step the tangential velocity and the normal stress are involved. These two algorithms are completely symmetric with the same numerical complexity and the coupling interface conditions have a physical meaning. The above results thus lead to the following question: Are those algorithms the only ones? If not, study the dependency of the algorithms with respect to F .

An algorithmic approach

As we have seen in the above sections, the efficiency of our algorithms heavily relies on the simplicity of the Smith variables, that is of the rows of the matrix F already defined. In this section, within a constructive algebraic analysis approach (see [START_REF] Chyzak | Effective algorithms for parametrizing linear control systems over Ore algebras[END_REF][START_REF] Cluzeau | Factoring and decomposing a class of linear functional systems[END_REF][START_REF] Quadrat | An introduction to constructive algebraic analysis and its applications[END_REF]), we develop a method for constructing many possible Smith variables. Taking into account physical aspects, the user can then choose the simplest ones among them.

Let A ∈ R p×p be a matrix with entries in the ring R = k[∂ 1 , . . . , ∂ d ] of PD operators with coefficients in a field k. An element P ∈ R has the form P = 0≤|µ|≤r a µ ∂ µ , where a µ ∈ k,

Inria µ = (µ 1 . . . µ n ) ∈ N 1×n , ∂ µ = ∂ µ1 1 . . . ∂ µn n , ∂ i = ∂ ∂xi , i = 1, . . . , n, and |µ| = µ 1 + • • • + µ n .
Moreover, let us introduce the following commutative polynomial rings:

     R 1 = k(∂ 2 , . . . , ∂ d )[∂ 1 ], R i = k(∂ 1 , . . . , ∂ i-1 , ∂ i+1 , . . . , ∂ d )[∂ i ], i = 2, . . . , d -1, R d = k(∂ 1 , . . . , ∂ d-1 )[∂ d ]. ( 30
)
We recall that the Smith normal form of A with respect to the direction x i is the diagonal matrix D = diag(d 1 , . . . , d p ) formed by monic polynomials

d i of R i satisfying d 1 | d 2 | • • • | d p and A = E D F , for certain matrices E, F ∈ GL p (R i ).
In particular, there exists 0 ≤ r ≤ p, such that the diagonal matrix D can be written as:

D = I r 0 0 D 2 , D 2 = diag(d r+1 , . . . , d p ) ∈ R (p-r)×(p-r) i , d r+1 = 1.
As it was previously pointed out, the matrices E and F are not unique. Computer algebra systems such as Maple provide a pair of matrices (E, F ). The goal of this section is to explain how to provide many distinct F 's so that the results of Section 4 can be applied on the simplest F 's or on the F 's leading to Smith variables having a physical meaning. This way, the procedure described in Section 4 will be made automatic. To do that, we first need to recall a few definitions on module theory and state useful material. For more details on a module-theoretic approach to mathematical systems theory, see [START_REF] Quadrat | An introduction to constructive algebraic analysis and its applications[END_REF] and references therein.

A few results of module theory

If R is a commutative ring, then a R-module M is an abelian group (M, +) (see, e.g., [START_REF] Buchsbaum | What annihilates a module?[END_REF]) equipped with a scalar multiplication

R × M -→ M (r, m) -→ r m,
which satisfies the following properties:

1. r 1 (m 1 + m 2 ) = r 1 m 1 + r 1 m 2 , 2. (r 1 + r 2 ) m 1 = r 1 m 1 + r 2 m 1 , 3. r 2 (r 1 m 1 ) = (r 2 r 1 ) m 1 , 4. 1 m 1 = m 1 ,
for all r 1 , r 2 ∈ R and for all m 1 , m 2 ∈ M . Hence, the definition of a R-module M is similar to the one of a vector space but where the scalars are taken in a ring R and not in a field (e.g., Q, R, C) as for vector spaces. For instance, an abelian group is a Z-module and an ideal of R is a R-submodule of R.

If M and N are two R-modules, then a R-homomorphism f : M -→ N is a R-linear map from M to N , namely: ∀ r 1 , r 2 ∈ R, ∀ m 1 , m 2 ∈ M, f (r 1 m 1 + r 2 m 2 ) = r 1 f (m 1 ) + r 2 f (m 2 ).
We denote by hom R (M, N ) the abelian group formed by the R-homomorphisms from M to N . Since R is a commutative ring, hom R (M, N ) inherits a R-module structure:

∀ f ∈ hom R (M, N ), ∀ r ∈ R, ∀ m ∈ M, (r f )(m) = f (r m). Now, f ∈ hom R (M, N ) is called an isomorphism if f is a bijective (i.e.
, an injective and surjective) homomorphism. If there exists an isomorphism from M to N , then M and N are said to be isomorphic, which is denoted by M ∼ = N .

If A ∈ R q×p , then we can consider the following R-homomorphism:

.A : R 1×q -→ R 1×p r = (r 1 . . . r q ) -→ r A.

The kernel of the R-homomorphism .A is the R-module defined by:

ker R (.A) = {r ∈ R 1×q | r A = 0}.
The matrix A is said to have full row rank if ker R (.A) = 0, i.e., if the rows of A are R-linearly independent. The image im R (.A) of .A, simply denoted by R 1×q A, is the R-module defined by all the R-linear combinations of the rows of A. Moreover, the cokernel coker R (.A) of .A is the factor R-module defined by coker R (.A) := R 1×p /(R 1×q A). To simplify the notation, we shall denote this module by M , i.e., M = R 1×p /(R 1×q A). Two vectors r, r ′ ∈ R 1×p are said to be equivalent, denoted by r ∼ r ′ , if there exists s ∈ R 1×q such that r = r ′ + s A. We can easily check that ∼ is an equivalence class. The residue class of r ∈ R 1×p for this equivalence relation

∼ is denoted π(r), i.e., π(r) = π(r ′ ) if there exists s ∈ R 1×q such that r = r ′ + s A. The R-module structure of M is defined by: ∀ r, r ′ ∈ R 1×p , ∀ r ∈ R, π(r) + π(r ′ ) := π(r + r ′ ), r π(r) := π(r r). (31)
We can easily check that the above operations are well-defined, i.e., they do not depend on the choice of the representative r and r ′ of the elements π(r) and π(r ′ ) of M . Indeed, if t, t ′ ∈ R 1×p are such that π(r) = π(t) and π(r ′ ) = π(t ′ ), then there exist s, s ′ ∈ R 1×q such that r = t + s A and r ′ = t ′ + s ′ A, which yields:

π(r + r ′ ) = π(t + t ′ + (s + s ′ ) A) = π(t + t ′ ), π(r r) = π(r t + (r s) A) = π(r t).
Since the R-module M = R 1×p /(R 1×q A) plays a fundamental role in what follows, we describe it in terms of generators and relations. Let f j be the j th vector of the standard basis {f j } j=1,...,p of R 1×p , namely, f j is the row vector of length p defined by 1 at the j th position and 0 elsewhere. Moreover, let π : R 1×p -→ M = R 1×p /(R 1×q A) be the R-homomorphism sending r ∈ R 1×p to its residue class π(r) in M (see (31)). We claim that {m j = π(f j )} j=1,...,p is a family of generators of the R-module M . Indeed, for any m ∈ M , there exists r = (r 1 . . . r p ) ∈ R 1×p such that:

m = π(r) = π   p j=1 r j f j   = p j=1 r j π(f j ) = p j=1 r j m j .
Now, since f i A is the i th row of the matrix A, π(f i A) = 0, which yields:

∀ i = 1, . . . , q, π(f i A) = π   p j=1 A ij f j   = p j=1 A ij π(f j ) = p j=1 A ij m j = 0. ( 32
)
Thus, the family of generators {m j } j=1,...,p of M satisfies the relations p j=1 A ij m j = 0 for all i = 1, . . . , q. The R-module M is then said to be finitely presented by A ∈ R q×p and the matrix A is a presentation matrix of M . Note that the above remarks do not assume the commutativity Inria of the ring R. For more details, see [START_REF] Chyzak | Effective algorithms for parametrizing linear control systems over Ore algebras[END_REF][START_REF] Quadrat | An introduction to constructive algebraic analysis and its applications[END_REF] and the references therein. Module theory and homological algebra study properties of the R-module M which do not depend on the particular presentation M = R 1×p /(R 1×q A) of M (up to isomorphism, the R-module M can be presented by different matrices of different sizes).

Example 1. The finitely presented R-module M = R 1×p /(R 1×q A) is classically used in number theory and in algebraic geometry. For instance, Cauchy's definition of the field of complex numbers is C = R[x]/(R[x] (x 2 + 1)). In this case, we can take R = R[x], A = (x 2 + 1), and p = q = 1. Similarly, the ring of Gaussian numbers

Z[i] := {a + b i | a, b ∈ Z, i 2 = -1} can be written as Z[i] = Z[x]/(x 2 + 1), i.e., R = Z[x], A = (x 2 + 1), and p = q = 1. Now, if R = k[x 1 , . . . , x d ]
is a commutative ring of polynomials in x 1 , . . . , x d with coefficients in a field k, and I is an ideal of R generated by the polynomials P 1 , . . . , P q , i.e.,

I = q i=1 R P i = q i=1 r i P i | r i ∈ R, i = 1, . . . , q ,
then an important concept in algebraic geometry is the affine coordinate ring R/I associated with the affine algebraic variety

V (I) = {(x 1 , . . . , x d ) ∈ k d | ∀ i = 1, . . . , q, P i (x 1 , . . . , x d ) = 0}
defined by the common zeros of the polynomials P 1 , . . . , P q in k d . In this case, R = k[x 1 , . . . , x d ], A = (P 1 . . . P q ) T ∈ R q×1 , i.e., p = 1, and

R 1×q A = q i=1 R P i = I. Let us now state one interest of the use of the R-module M = R 1×p /(R 1×q A) in mathematical systems theory. If F is a R-module, then a result due to Malgrange ([25]) asserts that: ker F (A.) := {η ∈ F p | A η = 0} ∼ = hom R (M, F). (33)
Let us explain [START_REF] Roeck | Analysis and test of a local domain decomposition preconditioner[END_REF]. If {m j = π(f j )} j=1,...,p is the family of generators of the R-module M = R 1×p /(R 1×q A) defined above, then φ ∈ hom R (M, F) is defined by φ(m j ) = η j ∈ F for j = 1, . . . , p. Moreover, since by definition of a homomorphism, φ(0) = 0, using [START_REF] Quadrat | An introduction to constructive algebraic analysis and its applications[END_REF], we get

∀ i = 1, . . . , q, p j=1 A ij η j = p j=1 A ij φ(m j ) = φ   p j=1 A ij m j   = φ(0) = 0,
which shows that A η = 0, where η = (η 1 . . . η p ) T ∈ F p , and thus η ∈ ker F (A.). Conversely, if η ∈ ker F (A.), i.e., A η = 0, then we can consider the map φ η : M -→ F defined by φ η (π(r)) = r η for all r ∈ R 1×p . Indeed, this map is well-defined since if π(r) = π(r ′ ), then there exists s ∈ R 1×q such that r = r ′ + s A, which yields r η = r ′ η + s A η = r ′ η and shows that φ η (π(r)) = r η depends only on π(r) and not on any representative r ∈ R 1×p of π(r). Moreover,

∀ i = 1, . . . , q, φ η p i=1 A ij m j = p i=1 A ij φ η (m j ) = p i=1 A ij η j = 0,
which shows that φ η ∈ hom R (M, F) and proves the isomorphism since for every η ∈ ker F (A.), φ η (m j ) = η j for all j = 1, . . . , p, and for all φ

∈ hom R (M, F), φ = φ (φ(m1) ... φ(mp)) T . Example 2. If R = Q[∂ 1 , ∂ 2 ] is the commutative polynomial ring of PD operators with rational constant coefficients, ∆ = ∂ 2 1 + ∂ 2 2 ∈ R, M = R/(R ∆), and F = C ∞ (Ω) (resp., D ′ (Ω), S ′ (Ω)
) the R-module of smooth functions (resp., distributions, temperate distributions) on Ω, where Ω is an open subset of R 2 , then the isomorphism (33) yields:

ker F (∆.) = {η ∈ F | ∆ η = 0} ∼ = hom R (M, F).

Let us give a module-theoretic interpretation of the differential operators appearing in the Smith normal form of the presentation matrix

A of the R-module M (see Section 2). If M is a R-module, then the annihilator ann R (m) of m ∈ M is the ideal of R defined by ann R (m) = {r ∈ R | r m = 0}. More generally, the annihilator ann R (M ) of the R-module M is the ideal of R defined by: ann R (M ) := {r ∈ R | ∀ m ∈ M : r m = 0} = m∈M ann R (m). ( 34
)
Let M = R 1×p /(R 1×q A) be the R-module finitely presented by A ∈ R q×p . Let us compute the annihilator of M . If r ∈ ann R (M ), then r m = 0 for all m ∈ M . In particular, we get r m j = 0 for all j = 1, . . . , p, where {m j = π(f j )} j=1,...,p is a set of generators of M . Conversely, if r ∈ R is such that r m j = 0 for all j = 1, . . . , p, then r m = 0 for all m ∈ M since R is a commutative ring and every element m ∈ M has the form of m = p j=1 r j m j for certain r j ∈ R for j = 1, . . . , p. Now, for j = 1, . . . , p, π(r f j ) = r π(f j ) = r m j = 0, and thus there exists

B j ∈ R 1×q such that r f j = B j A, i.e., r I p = B A, where B = (B T 1 . . . B T p ) T ∈ R p×q ,

and thus:

ann

R (M ) = {r ∈ R | ∃ B ∈ R p×q : r I p = B A}.
Let us now focus on full row rank square matrices, a class which will be considered for applications in what follows. Hence, let A ∈ R p×p have full row rank, i.e., det(A) = 0. Similarly to what happens in linear algebra, we can check that A Adj(A) = Adj(A) A = det(A) I p hold, where Adj(A) ∈ R p×p is the adjugate matrix of A, namely the transpose of the matrix of cofactors of A (see, e.g., [START_REF] Eisenbud | Commutative Algebra with a View Toward Algebraic Geometry[END_REF]). We recall that the matrix of cofactors of A is the matrix whose (i, j) entry is (-1) i+j m ij , where m ij is the minor for entry A ij , namely the determinant of the (p-1)×(p-1)submatrix of A obtained by removing the i th row and the j th column of A. If r ∈ ann R (M ), i.e., r I p = B A for a certain matrix B ∈ R p×p , then r Adj(A) = B A Adj(A) = B det(A), i.e., r Adj(A) ij = B ij det(A) for i, j = 1, . . . , p, which yields (the converse implication holds since det(A) is a non-zero divisor of R):

r ∈ ann R (M ) ⇔ ∀ i, j = 1, . . . , p, r m ij ∈ R det(A) := {s det(A) | s ∈ R}. ( 35
)
Let us now introduce the important concept of Fitting's ideals in module theory ( [START_REF] Buchsbaum | What annihilates a module?[END_REF][START_REF] Eisenbud | Commutative Algebra with a View Toward Algebraic Geometry[END_REF]).

Definition 1. Let R be a commutative ring, A ∈ R q×p , and M = R 1×p /(R 1×q A). Then, the ideal of R generated by the (pi) × (pi)-minors of A is denoted Fitt i (M ) and is called the i th Fitting ideal of the R-module M . By convention, Fitt i (M ) = 0 if pi > q, i.e., i < pq, and Fitt i (M ) = R for i ≥ p.

It can be proved that the Fitting ideals depend only on M and not on the matrix A which presents M (see, e.g., [START_REF] Eisenbud | Commutative Algebra with a View Toward Algebraic Geometry[END_REF]), a fact explaining the notation Fitt i (M ).

Since every i × i-minor of A is a R-linear relation of (i -1) × (i -1)-minors of A, we get:

0 ⊆ Fitt 0 (M ) ⊆ Fitt 1 (M ) ⊆ Fitt 2 (M ) ⊆ • • • ⊆ Fitt i (M ) ⊆ Fitt i+1 (M ) ⊆ • • • ⊆ R. ( 36
)
Let A ∈ R p×p have full row rank, i.e., det(A) = 0, and M = R 1×p /(R 1×p A). Then:

Fitt 0 (M ) = R det(A), Fitt 1 (M ) = p i,j=1
R m ij .

Inria

If I and J are two ideals of R such that I ⊆ J, then I : J = {r ∈ R | r J ⊆ I} is the conductor of J into I, where r J := {r j | j ∈ J} is the principal ideal of R generated by r. We note that I : J = ann R (J/I).

Then, (35) shows that:

ann R (M ) = {r ∈ R | ∀ i, j = 1, . . . , p, r m ij ∈ R det(A)} = {r ∈ R | r Fitt 1 (M ) ⊆ Fitt 0 (M )} = Fitt 0 (M ) : Fitt 1 (M ). ( 37
)
Let us now explain how to explicitly compute the annihilator ann R (M ) of M when R is a commutative polynomial ring with coefficients in a computational field (e.g.,

k = Q). If Fitt 1 (M ) = p 2 i=1 R r i (i.e.
, r i is one of the m jk 's) and L = (r 1 . . . r p 2 ) ∈ R 1×p 2 , then using (37) and the commutativity of R, we get that r ∈ ann R (M ) iff there exist s i ∈ R for i = 1, . . . , p 2 such that r r i = s i det(A), i.e., iff:

L r = (s 1 . . . s p 2 ) det(A) ⇔ (r -s 1 . . . -s p 2 ) L det(A) I p 2 = 0.
Using Gröbner basis techniques (see Algorithm 9 of Section 7), we can compute

W ∈ R s×(1+p 2 ) such that ker R (.(L T det(A) I T p 2 ) T ) = R s×(1+p 2 )
W , and thus we get:

ann

R (M ) = s k=1 R W k1 .
Remark 4. Using (37), we get Fitt 0 (M ) ⊆ ann R (M ). More generally, we can prove (see [START_REF] Eisenbud | Commutative Algebra with a View Toward Algebraic Geometry[END_REF])

∀ j ≥ 1, ann R (M ) Fitt j (M ) ⊆ Fitt j-1 (M ), (38)
where the product I J of two ideals I and J of R is defined by:

I J = t i=1 a i b i | t ∈ N, a i ∈ I, b i ∈ J .
If the R-module M can be generated by r elements, then there exists a matrix B ∈ R s×r such that M ∼ = R 1×r /(R 1×s B), which shows that Fitt r (M ) = R since the ideal generated by the 0 × 0-minors of the matrix B is R and the Fitting ideals do not depend on the presentation of M (see [START_REF] Eisenbud | Commutative Algebra with a View Toward Algebraic Geometry[END_REF]). Moreover, (38) yields ann R (M ) r Fitt r (M ) ⊆ Fitt 0 (M ), i.e., ann R (M ) r ⊆ Fitt 0 (M ).

In particular, if M is a cyclic R-module, namely, M can be generated by one generator, i.e., r = 1, then ann R (M ) = Fitt 0 (M ).

An important application of the characterization of the annihilator ann

R (M ) of the R-module M = R 1×p /(R 1×p A) is the computation of the so-called characteristic variety char R (M ) of M . Definition 2. Let R = R[∂ 1 , . . . , ∂ d ] be the commutative polynomial ring of PD operators with coefficients in R, A ∈ R p×p , and M = R 1×p /(R 1×p A).
Then, the characteristic ideal of the R-module M is defined by:

I(M ) = ann R (M ) := {r ∈ R | ∀ m ∈ M, ∃ l ∈ N : r l m = 0}.
The characteristic variety of M is the complex affine algebraic variety defined by I(M ), i.e.:

char

R (M ) = {χ = (χ 1 , . . . , χ d ) ∈ C d | ∀ r ∈ I(M ) : r(χ) = 0}. Now, if M is a R-module and d ∈ R, then the annihilator ann M (d) is the R-submodule of M defined by ann M (d) = {m ∈ M | d m = 0}. If M = R 1×p /(R 1×p A
) and m ∈ ann M (d), then there exists r ∈ R 1×p such that m = π(r), which yields:

π(d r) = d π(r) = d m = 0 ⇔ ∃ s ∈ R 1×p : d r = s A. Since R is a commutative ring, d r = s A is equivalent to: (r -s) d I p A = 0.
If R is a commutative polynomial ring with coefficients in a computable field k, then we can find a family of generators of the R-module ker R (.(d I T p A T ) T ), i.e., there exists a matrix S ∈ R q×2 p such that ker R (.

(d I T p A T ) T ) = R 1×q S (see Algorithm 9 of Section 7). If S = (U V ), where U ∈ R q×p and V ∈ R q×p , then ann M (d) = q k=1 R π(U k• )
, where U k• is the k th row of the matrix U . Computing the normal forms of the vectors U k• 's with respect to a Gröbner basis of the R-module R 1×p A for a certain monomial order (for more details, see Section 7), we obtain

u k ∈ R 1×p satisfying ann M (d) = q k=1 R π(u k ).
All these computations can be handled with the OreModules package [START_REF] Chyzak | OreModules: A symbolic package for the study of multidimensional linear systems[END_REF].

Let us now show that the Fitting ideals Fitt

i (M )'s of the R-module M = R 1×p /(R 1×p A) are finer invariants than the invariant factors d i 's of the Smith normal form of A ∈ R p×p defined over the principal ideal domain R i . If R j ⊗ R M denotes the R j -module finitely presented by A ∈ R p×p , namely, R j ⊗ R M := R 1×p j /(R 1×p j A), (39)
then we can easily check (see, e.g., [START_REF] Eisenbud | Commutative Algebra with a View Toward Algebraic Geometry[END_REF]) that:

Fitt i (R j ⊗ R M ) = R j Fitt i (M ) := t k=1 r k s k | t ∈ N, r k ∈ R j , s k ∈ Fitt i (M ) . ( 40
)
An important property of the ring R j is that it is a principal ideal domain, namely, a ring for which every ideal is principal, i.e., can be generated by one element (see, e.g., [START_REF] Buchsbaum | What annihilates a module?[END_REF]). For instance, if

I = R j r 1 + • • • + R j r m = { m k=1 s k r k | s k ∈ R j } is
an ideal of R j generated by the elements r 1 , . . . , r m of R j , then I = R j t, where t is the greatest common divisor of the r i 's. The element t can be computed by means of the Euclidean algorithm. In particular, the ideals Fitt p-i (R j ⊗ R M ) are principal, i.e., Fitt p-i (R j ⊗ R M ) = R j m i , where m i is the greatest common divisor of the i × i-minors of A over the principal ideal domain R j . Hence, each ideal in the chain [START_REF] Wloka | Boundary Value Problems for Elliptic Systems[END_REF]

is principal and Fitt

p-i (R j ⊗ R M ) = R j m i ⊆ Fitt p-i+1 (R j ⊗ R M ) = R j m i-1 yields that m i-1 divides m i .
According to Theorem 1, we obtain that the invariant factors of A ∈ R p×p j over R j are exactly the generators of the following R i -modules:

∀ i = 1, . . . , p, Fitt p-i (R j ⊗ M )/Fitt p-i+1 (R j ⊗ M ) = R j (m i /m i-1 ) = R j d i .
Let us illustrate the above results on the Cauchy-Navier and Oseen equations.

Example 3. Let R = Q(λ, µ)[∂ x , ∂ y , ∂ z ]
be the commutative polynomial ring in ∂ x , ∂ y and ∂ z with coefficients in the field Q(λ, µ), the matrix A = A 3 ∈ R 3×3 defined by [START_REF] Chyzak | Non-commutative elimination in Ore algebras proves multivariate identities[END_REF], and the finitely presented R-module M = R 1×3 /(R 1×3 A). If {f j } j=1,2,3 is the standard basis of R 1×3 and π : R 1×3 -→ M the canonical projection onto M , then {m j = π(f j )} j=1,2,3 is a family of Inria generators of M which satisfies the relations 3 i=1 A ij m j = 0 for all i = 1, 2, 3. We can easily check that:

         Fitt 0 (M ) = R ∆ 3 , Fitt 1 (M ) = (R ∆)(R ∂ x ∂ y + R ∂ x ∂ z + R ∂ y ∂ z + R (µ ∂ 2 x + (λ + 2 µ) (∂ 2 y + ∂ 2 z )) +R ((λ + 2 µ) (∂ 2 x + ∂ 2 z ) + µ ∂ 2 y )) + R ((λ + 2 µ) (∂ 2 x + ∂ 2 y ) + µ ∂ 2 z ))), Fitt 2 (M ) = R. (41) Then, ann R (M ) = Fitt 0 (M ) : Fitt 1 (M ) = R ∆ 2 , i.e., every element m of M is annihilated by ∆ 2 , i.e., ∆ 2 m = 0, which shows that: char R (M ) = {(χ x , χ y , χ z ) ∈ C 3 | χ 2 x + χ 2 y + χ 2 z = 0}.
Let us now characterize the elements of M which are annihilated by ∆, i.e., let us compute ann M (∆). Using Gröbner basis techniques, we get that ker R (.(∆ I T 3 A T ) T ) = R 1×4 S, where S = (U V ) and:

U =      -µ ∂ z 0 µ ∂ x -µ ∂ y µ ∂ x 0 0 -µ ∂ z µ ∂ y (λ + 2 µ) ∂ x (λ + 2 µ) ∂ y (λ + 2 µ) ∂ z      , V =      -∂ z 0 ∂ x -∂ y ∂ x 0 0 -∂ z ∂ y ∂ x ∂ y ∂ z     
.

Thus, we obtain:

ann M (∆) =R (-∂ z m 1 + ∂ x m 3 ) + R (-∂ y m 1 + ∂ x m 2 ) + R (-∂ z m 2 + ∂ y m 3 ) + R (∂ x m 1 + ∂ y m 2 + ∂ z m 3 ).
Similarly, we can easily check again that ann

M (∆ 2 ) = R m 1 + R m 2 + R m 3 = M .
If we now consider the ring

R 1 = Q(λ, µ)(∂ y , ∂ z )[∂ x ] = Q(λ, µ, ∂ y , ∂ z )[∂ x ] of OD operators in ∂ x with coefficients in the field Q(λ, µ, ∂ y , ∂ z ) and R 1 ⊗ R M = R 1×3 1 /(R 1×3 1 A), then (40) and (41) give      Fitt 0 (R 1 ⊗ R M ) = R 1 ∆ 3 , Fitt 1 (R 1 ⊗ R M ) = R 1 ∆, Fitt 2 (R 1 ⊗ R M ) = R 1 , since the element ∂ y ∂ z appearing in Fitt 1 (M ) (see (41)) is an invertible element of R 1 , and thus, ann R1 (R 1 ⊗ R M ) = (∆ 3 : ∆) = R 1 ∆ 2 .
Finally, we can check again that the generators of the following principal ideals

     Fitt 2 (R 1 ⊗ R M ) = R 1 , Fitt 1 (R 1 ⊗ R M )/Fitt 2 (R 1 ⊗ R M ) = R 1 ∆, Fitt 0 (R 1 ⊗ R M )/Fitt 1 (R 1 ⊗ R M ) = R 1 ∆ 2 ,
are exactly the invariants factors of the matrix A (see [START_REF] Cluzeau | Factoring and decomposing a class of linear functional systems[END_REF]). [START_REF] Cluzeau | OreMorphisms: a homological algebraic package for factoring, reducing and decomposing linear functional systems[END_REF], and the R-module M = R 1×4 /(R 1×4 A) finitely presented by O 3 . If {f j } j=1,...,4 is the standard basis of R 1×4 and π : R 1×4 -→ M the canonical projection onto M , then {m j = π(f j )} j=1,...,4 is a family of generators of M satisfying the relations where L 3 is the PD operator defined by [START_REF] Cox | Using Algebraic Geometry[END_REF], and:

Now, we consider R = Q(b 1 , b 2 , c, ν)[∂ x , ∂ y , ∂ z ], the matrix A = O 3 ∈ R 4×4 defined by
char R (M ) ={(χ x , χ y , χ z ) ∈ C 3 | χ 2 x + χ 2 y + χ 2 z = 0} ∪ {(χ x , χ y , χ z ) ∈ C 3 | c -ν (χ 2 x + χ 2 y + χ 2 z ) + b 1 χ x + b 2 χ y + b 3 χ z = 0}.
Finally, we have:

ann M (∆) = R m 4 , ann M (L 3 ) = R (-∂ z m 1 + ∂ x m 3 ) + R (-∂ y m 1 + ∂ x m 2 ) + R (-∂ z m 2 + ∂ y m 3 ) +R (∂ z (b 1 m 1 + b 2 m 2 + b 3 m 3 + m 4 ) + c m 3 ) +R (b 1 ∂ y m 1 + (-ν ∂ 2 z + b 2 ∂ y + b 3 ∂ z + c) m 2 + ν ∂ y ∂ z m 3 + ∂ y m 4 ).
In this particular example, the computation of ann R (M ) allows us to obtain in an intrinsic way the PD operators with their multiplicities appearing in the Smith normal form of the presentation matrix A of M . See Section 2.

Finally, since Fitt 0 (M ) ann R (M ), Remark 4 shows that Cauchy-Navier equations and Oseen equations do not define cyclic R-modules, and thus they are not equivalent to a linear system of PDEs with constant coefficients in only one unknown. Definition 3. A R-module M is said to be free if M admits a basis, namely, an independent family of generators of M . If M is a finitely generated R-module, then it means that M is isomorphic to R l for a certain l ∈ N. Then, l is called the rank of M and it is denoted by rank R (M ).

Theorem 2 ( [START_REF] Fabiańska | Applications of the Quillen-Suslin theorem to multidimensional systems theory[END_REF]). Let R be a commutative integral domain, A ∈ R q×p a full row rank matrix, i.e., ker R (.A) = 0, and M = R 1×p /(R 1×q A) the R-module finitely presented by A.

1. M is a free R-module iff there exist three matrices

B ∈ R p×q , C ∈ R (p-q)×p , D ∈ R p×(p-q)
such that the following two Bézout identities hold:

A C (B D) = I p , (B D) A C = I p . (42)
If π : R 1×p -→ M is the canonical projection onto M , then {π(C i• )} i=1,...,p-q is a basis of the free R-module M of rank pq.

2. If k is a field and R = k[∂ 1 , . . . , ∂ d ] a commutative polynomial ring with coefficients in k, then M is a free R-module of rank pq iff the matrix A admits a right inverse, namely, iff there exists a matrix B ∈ R p×q such that A B = I q (Quillen-Suslin theorem).

We first note that ( 42) is equivalent to the following identities:

A B = I q , A D = 0, C B = 0, C D = I p-q , B A + D C = I p . ( 43
)
The first part of Theorem 2 shows that the R-module M = R 1×p /(R 1×q A), finitely presented by a full row rank matrix A, is free of rank pq iff A can be completed to a square unimodular matrix, namely, iff there exists C ∈ R (p-q)×p such that:

A C ∈ GL p (R). ( 44
)

Inria

The problem of completing the matrix A to a unimodular matrix is called a completion problem.

The second part of Theorem 2 states that every full row rank matrix A with entries in a commutative polynomial ring D = k[∂ 1 , . . . , ∂ d ] over a field k, which admits a right inverse can be completed to a unimodular matrix, i.e., there exists C ∈ R (p-q)×p such that (44) holds. Despite of the simplicity of the formulation of the Quillen-Suslin theorem, this result, first raised by Serre in 1955, was only proved in 1976. The computation of such a matrix C (i.e., of a basis {π(C i• )} i=1,...,p-q of the free D-module M of rank pq) is generally a difficult issue which requires a constructive proof of the Quillen-Suslin theorem (see [START_REF] Fabiańska | Applications of the Quillen-Suslin theorem to multidimensional systems theory[END_REF] and the references therein).

The only simple case is

d = 1, i.e., R = k[∂ 1 ]
, where the computation of the Smith normal form of A gives a matrix C satisfying the completion problem (see [START_REF] Fabiańska | Applications of the Quillen-Suslin theorem to multidimensional systems theory[END_REF][START_REF] Quadrat | An introduction to constructive algebraic analysis and its applications[END_REF]). Let us recall this result. Let R be a principal left ideal domain and A ∈ R q×p a matrix admitting a right inverse B ∈ R p×q . Then, A has necessarily full row rank since:

λ A = 0 ⇒ λ = λ (A B) = 0. Computing a Smith normal form of A over R, we get F ∈ GL q (R) and G ∈ GL p (R) satisfying A = F D G, where D = (D 1 0) ∈ R q×p , D 1 = diag(d 1 , . . . , d q ) and d i ∈ R for i = 1, . . . , q.
Since A has full row rank, so is D 1 , i.e., d i = 0 for i = 1, . . . , q. Now, (F D G) B = A B = I q , and thus D (G B) = F -1 , i.e., D (G B F ) = I q , which proves that D admits the right inverse

E = G B F = (E T 1 E T 2 ) T
, where E 1 ∈ R q×q and E 2 ∈ R (p-q)×q . Thus, we get D 1 E 1 = I q , which shows that D 1 ∈ GL q (R), and thus all the d i 's are invertible elements of R, which can be assumed without loss of generality to be equal to 1. Hence, we can assume that D 1 = I q and:

A = F (I q 0) G. If r = p -q, G = (G T 1 G T 2 ) T ∈ GL p (R), where G 1 ∈ R q×p , G 2 ∈ R r×p , and G -1 = (H 1 H 2 )
, where H 1 ∈ R p×q , H 2 ∈ R p×r , then we obtain A = F G 1 , i.e., G 1 = F -1 A, and

F -1 A G 2 G -1 = I p ⇒ F -1 0 0 I r A G 2 G -1 = I p , ⇒ A G 2 G -1 F -1 0 0 I r = I p ⇒ A G 2 (H 1 F -1 H 2 ) = I p ,
which shows that we can take

C = G 2 and (B D) = (H 1 F -1 H 2 ) ∈ GL p (R). Hence, if R is principal ideal domain (e.g., R = k[∂ 1 ],
where k is a field, or if R = R i defined by ( 30)), then the completion problem can be solved by means of a Smith normal form computation.

. The R = k[∂ x , ∂ y , ∂ z ]-module M = R 1×2 /(R (∆ -1)
) is defined by the two generators m 1 = π(f 1 ) and m 2 = π(f 2 ), where {f 1 , f 2 } is the standard basis of R 1×2 and π : R 1×2 -→ M is the canonical projection onto M , and one R-linear relation ∆ m 1m 2 = 0. By 1 of Theorem 2, the R-module is free of rank 1 since:

∆ -1 1 0 ∈ GL 2 (R).
In particular, the residue class m 1 of (1 0) in M is a basis of M , i.e., M = R m 1 .

RR n°7953

A more complicated example is given by the

R = k[∂ x , ∂ y , ∂ z ]-module M = R 1×3 /(R A), where A = (∂ x ∂ y -1 ∂ 2 x ∂ 2 y ). We can check that    ∂ x ∂ y -1 ∂ 2 x ∂ 2 y ∂ 2 y ∂ x ∂ y + 1 0 0 0 1    ∈ GL 3 (R),
which shows that the R-module M defined by the generators {m j } j=1,2,3 satisfying the R-linear

relation (∂ x ∂ y -1) m 1 + ∂ 2 x m 2 + ∂ 2 y m 3 =
0 is free of rank 2 and a basis of M is defined by

p 1 = ∂ 2 y m 1 + (∂ x ∂ y + 1) m 2 and p 2 = m 3 , i.e., M = R p 1 + R p 2 .
Indeed, computing the matrix D defined in [START_REF] Blinkov | The MAPLE Package "Janet": I. Polynomial Systems[END_REF], we then obtain:

     m 1 = ∂ 2 x p 1 + (∂ x ∂ y + 1) ∂ 2 y p 2 , m 2 = -(∂ x ∂ y -1) p 1 -∂ 4 y p 2 , m 3 = p 2 .

Computation of relevant Smith variables

We shall now show that the problem of finding Smith variables can be reduced to a particular completion problem.

Let M = R 1×p /(R 1×p A) be the R = k[∂ 1 , . . . , ∂ d ]-
module finitely presented by the matrix of PD operators A ∈ R p×p . With the notation [START_REF] Pavarino | Balancing Neumann-Neumann methods for incompressible Stokes equations[END_REF], computing the Smith normal form D ∈ R p×p i of the matrix A over the principal ideal domain R i , we obtain two matrices E ∈ GL p (R i ) and

F ∈ GL p (R i ) such that A = E D F . Now, if d 1 = • • • = d r = 1, d r+1 = 1 and
D 2 = diag(d r+1 , . . . , d p ) ∈ R (p-r)×(p-r) i , F = F 1 F 2 , E -1 = G 1 G 2 ,
where

F 1 ∈ R r×p i , F 2 ∈ R (p-r)×p i , G 1 ∈ R r×p i and G 2 ∈ R (p-r)×p i
, then:

A = E D F ⇔ F 1 = G 1 A, D 2 F 2 = G 2 A. ( 45
)
Lemma 2. With the above notations, if R i ⊗ R ann M (d j ) denotes the R i -module generated by the elements of ann M (d j ), then we have

R i ⊗ R ann M (d j ) = ann Ri⊗ R M (d j ). Proof. The first inclusion R i ⊗ R ann M (d j ) ⊆ ann Ri⊗RM (d j /1) is clear since for m ∈ M , d j m = 0 implies (d j m
)/1 = 0, i.e., (d j /1) (m/1) = 0, and thus m/1 belongs to ann Ri⊗ R M (d j ). Conversely, if m/1 belongs to ann Ri⊗M (d j /1), this means that (d j /1) (m/1) = (d j m)/1 = 0, and then it exists [START_REF] Buchsbaum | What annihilates a module?[END_REF]). Then,

s ∈ k[∂ 1 , . . . , ∂ j-1 , ∂ j+1 , . . . , ∂ n ] \ {0} such that s d j m = 0 in M (see
d j (s m) = 0 and s m ∈ ann M (d j ). Thus, (1/s) (s m) = m/1 holds in R i ⊗ R ann M (d j ), i.e., m/1 ∈ R i ⊗ R ann M (d j )
which leads to the second inclusion and ends the proof.

The last equality of (45) implies that the residue class n j of the j th row of the matrix F 2 in R i ⊗ R M must be annihilated by the PD operator d j ∈ R i \ {1}. Thus, using Lemma 2, every row of the submatrix Then, a natural problem is to seek for matrices

F 2 ∈ R (p-r)×p i of F ∈ GL p (R i ) satisfying A = E D F , where D is the Smith normal form of A and E ∈ GL p (R i ), is a R i -
F = (F T 1 F T 2 ) T ∈ GL p (R i ) -
where the residue classes of the rows of F 2 in R i ⊗ M (see (39)) correspond to the chosen Smith variables n j 's (e.g., ones having a physical meaning) which satisfy d j n j = 0 for j = r + 1, . . . , psuch that A = E D F for a certain E ∈ GL p (R i ), where D = diag(1, . . . , 1, d r+1 , . . . , d p ) is the Smith normal form of A. In other words, given a matrix F 2 ∈ R (p-r)×p such that the residue class n j of the j th row of 45), the matrices F 1 's are necessarily of the form F 1 = Λ A for certain matrices Λ ∈ R r×p i . Moreover, using (43), we can easily check that a necessary condition for

F 2 in R i ⊗ R M is annihilated by the PD operator d j ∈ R i \ {1}, we search for matrices F 1 ∈ R r×p i such that F = (F T 1 F T 2 ) T ∈ GL p (R i ) and A = E D F . From (
F 2 ∈ R (p-r)×p i to be completed to a unimodular matrix F is that F 2 admits a right inverse S 2 ∈ R p×(p-r) i , i.e., F 2 S 2 = I p-r .
In particular, it yields that ker Ri (.F 2 ) = 0 since λ F 2 = 0 then implies λ = λ F 2 S 2 = 0, i.e., the rows of F 2 are R i -linearly independent, in other words, F 2 has full row rank. Since R i is a commutative polynomial ring over a field, 2 of Theorem 2 proves that this necessary condition is also sufficient. Therefore, our problem can be restated as follows: Given F 2 ∈ R (p-r)×p i which admits a right inverse over R i and such that the residue class n j of the j th row of

F 2 in R i ⊗ R M is annihilated by the PD operator d j ∈ R i \ {1}, does it exist a matrix Λ ∈ R r×p i such that: F = ((Λ A) T F T 2 ) T ∈ GL p (R i ), A = E D F. (46)
If so, compute matrices Λ's satisfying (46).

Since the matrix F 2 ∈ R (p-r)×p i admits a right inverse over the commutative polynomial ring R i with coefficients in a field, by 2 of Theorem 2, the R i -module

N = R 1×p i /(R 1×(p-r) i
F 2 ) is free of rank r. Thus, by 1 of Theorem 2, there exist

T ∈ R r×p i , Q ∈ R p×r i , S ∈ R p×(p-r) i such that: T F 2 (Q S) = I r 0 0 I p-r = I p , Q T + S F 2 = I p . ( 47
)
As explained in 1 of Theorem 2, the residue classes of the rows of the matrix T in N define a basis of the free R i -module N . It is clear that the matrix T is not uniquely defined since the matrix T ′ = T + Z F 2 will also defines the same basis of N for all Z ∈ R r×(p-r) i

. Hence, we have to find (when they exist) two matrices Z ∈ R r×(p-r) i and Λ ∈ R r×p i such that:

T = (Λ -Z) A F 2 .
This means that the rows of T must belong to the R i -module

R 1×p i A + R 1×(p-r) i F 2 ,
A =    -(λ + µ) ∂ 2 x -µ ∆ -(λ + µ) ∂ x ∂ y -(λ + µ) ∂ x ∂ z -(λ + µ) ∂ x ∂ y -(λ + µ) ∂ 2 y -µ ∆ -(λ + µ) ∂ y ∂ z -(λ + µ) ∂ x ∂ z -(λ + µ) ∂ y ∂ z -(λ + µ) ∂ 2 z -µ ∆    , where ∆ = ∂ 2 x + ∂ 2 y + ∂ 2 z , and the R = Q(λ, µ)[∂ x , ∂ y , ∂ z ]-module M = R 1×3 /(R 1×3 A). The Smith normal form of the matrix A over R 1 = Q(λ, µ)(∂ y , ∂ z )[∂ x ] is defined by: D =    1 0 0 0 ∆ 0 0 0 ∆ 2    .
With the above notation, it means that r = 1 and

D 2 = diag(∆, ∆ 2 ) ∈ R 2×2 . If we write F = (F T 1 F T 2 ) T , where F 1 ∈ R 1×3 1 and F 2 ∈ R 2×3 1
, then the explanations above imply that the residue class of the first (resp., second) row of F 2 must be an element of the R 1 -module R 1 ⊗ R M annihilated by ∆ ∈ R (resp., ∆ 2 ∈ R). Using Example 3, we find that families of generators of ann M (∆) and ann M (∆ 2 ) are respectively given by the residue classes in R 1 ⊗ R M of the rows of the matrices A ∆ = U and A ∆ 2 = I 3 , where the matrix U is defined in Example 3. Now, the first (resp., second) row of F 2 must be a R 1 -linear combination of the rows of A ∆ (resp., A ∆ 2 ). We thus have several choices and for each of them, we are reduced to a completion problem. For instance, choosing the 4 th row of A ∆ (resp., the 3 rd row of A ∆ 2 ) as first (resp., second) row of F 2 , we obtain the following matrix

F 2 = 0 -∂ z ∂ y 0 0 1 ,
which trivially admits a right inverse with entries in R 1 . To solve (46), we have to find if it exists a row vector

Λ ∈ R 1×3 1 such that F = ((Λ A) T F T 2 ) T ∈ GL 3 (R 1
). If we can find such a row vector F 1 = Λ A, then the corresponding matrix F provides a choice of relevant Smith variables.

Let us now state the main results of this section.

Theorem 3. Let L = R 1×p /(R 1×q B) be a free R-module of rank r, Q ∈ R p×r and Q ′ ∈ R p×r two injective parametrizations of L, i.e., ker R (.Q) = R 1×p B, T Q = I r , ker R (.Q ′ ) = R 1×p B and T ′ Q ′ = I r ,
for some matrices T, T ′ ∈ R r×p . Then, there exist two matrices V ∈ GL r (R) and Z ∈ R r×q such that:

Q ′ = Q V, T ′ = V -1 T + Z B. (48)
Proof. Let κ : R 1×p -→ L be the canonical projection onto L. Then, the map φ : L -→ R 1×r , defined by φ(κ(r)) = r Q for all r ∈ R 1×p , is well-defined: if κ(r) = κ(s), then κ(rs) = 0, i.e., there exists t ∈ R 1×q such that r = s + t B, which yields φ(κ(r

)) = (s + t B) Q = s Q = φ(κ(s)). Hence, φ is a R-homomorphism from L to R 1×r . Since ker R (.Q) = R 1×q B, we get ker φ = {κ(r) | r Q = 0} = {κ(r) = κ(u B) | u ∈ R 1×q } = 0, i.e., φ is an injective R- homomorphism. Now, let ϕ : R 1×r -→ L be the R-homomorphism defined by ϕ(u) = κ(u T ) for all u ∈ R 1×r . Then, (φ • ϕ)(u) = φ(κ(u T )) = (u T ) Q = u for all u ∈ R 1×r , which shows that φ • ϕ = id R 1×r
, and proves that φ is surjective. Thus, φ is an R-isomorphism from L to R 1×r and ϕ = φ -1 . We note that if {h k } k=1,...,r is the standard basis of the free R-module R 1×r , then {κ(h k T ) = κ(T k• )} k=1,...,r is a basis of the free R-module L of rank r.

Similarly, the R-homomorphism ψ : L -→ R 1×r , defined by ψ(κ(r)) = r Q ′ for all r ∈ R 1×p , is an R-isomorphism, and ψ -1 : R 1×r -→ L is defined by

ψ -1 (u) = κ(u T ′ ) for all u ∈ R 1×r . The R-isomorphism ψ • φ -1 : R 1×r -→ R 1×r is defined by (ψ • φ -1 )(u) = u (T Q ′ ) for all u ∈ R 1×r . Since ψ • φ -1 is surjective, then for every vector h k of the standard basis of R 1×r , Inria there exists U k ∈ R 1×r such that h k = U k (T Q ′ ), which shows that U = (U 1 T . . . U r T ) ∈ R r×r satisfies U (T Q ′ ) = I r . Thus, det(U ) det(T Q ′ ) = 1, which shows that det(T Q ′ ) is a unit of R, and thus T Q ′ ∈ GL r (R). Hence, if we note V = T Q ′ ∈ GL r (R) and if α : R 1×r -→ R 1×r is the R-automorphism of R 1×r defined by α(v) = v V for all v ∈ R 1×r , then ψ = α • φ, i.e.: ∀ r ∈ R 1×p , ψ(κ(r)) = α(r Q) = r Q V.
Since by definition, ψ(κ(r)) = r Q ′ for all r ∈ R 1×p , we then get that

Q ′ = Q V . Finally, we have ψ -1 = φ -1 • α -1 , which yields κ(u T ′ ) = φ -1 (u V -1 ) = κ(u V -1 T ) for all u ∈ R 1×r , i.e., κ(u (T ′ -V -1 T )) = 0 for all u ∈ R 1×r , which proves the existence of a matrix Z ∈ R r×q such that T ′ -V -1 T = Z B, i.e., T ′ = V -1 T + Z B,

and finally proves (48).

We now obtain the following useful corollary of Theorem 3.

Corollary 1. Let A ∈ R p×p , L = R 1×p /(R 1×q B) be a free R-module of rank r, Q ∈ R p×r and Q ′ ∈ R p×r two injective parametrizations of L, i.e., ker R (.Q) = R 1×p B, T Q = I r , ker R (.Q ′ ) = R 1×p B and T ′ Q ′ = I r ,
for some T, T ′ ∈ R r×p . Then, the following assertions are equivalent:

1. There exist two matrices X ∈ R r×q and Y ∈ R r×p satisfying T + X B = Y A.

There exist two matrices

X ′ ∈ R r×q and Y ′ ∈ R r×p satisfying T ′ + X ′ B = Y ′ A.
Hence, T can be left factorized by (B T A T ) T iff T ′ can be left factorized by (B T A T ) T .

Proof. Using (48), there exist V ∈ GL r (R) and Z ∈ R r×q such that T ′ = V -1 T + Z B. Let us suppose that 1 holds. Then, using the identity T = V T ′ + (-V Z) B, we get:

T + X B = Y A ⇔ T ′ + (V -1 X -Z) B = (V -1 Y ) A.
Similarly, if 2 holds, then we obtain:

T ′ + X ′ B = Y ′ A ⇔ T + (V (X ′ + Z)) B = (V Y ′ ) A.
Using Corollary 1, we deduce an algorithm that given a matrix A ∈ R p×p and a matrix F 2 ∈ R (p-r)×p such that the residue class of the i th row of

F 2 in R i ⊗ R M = R 1×p i /(R 1×p i A)
defines an element of ann Ri⊗ R M (d r+i ), where d r+i ∈ R is the numerator of the i th diagonal element the Smith normal form D of A over R i , find a completion if it exists. Algorithm 6. Input: A computable field k, R = k[∂ 1 , . . . , ∂ d ], the ring R i defined by [START_REF] Pavarino | Balancing Neumann-Neumann methods for incompressible Stokes equations[END_REF], a full row rank matrix A ∈ R p×p , and F 2 ∈ R (p-r)×p defined as explained above.

Output: A completion F of F 2 if it exists and [ ] if not. 1. Compute a right inverse of F 2 over R i .
2. If such a right inverse does not exist, then return [ ], else do:

(a) Compute an injective minimal parametrization Q ∈ R p×r i of the free R i -module L = R 1×p i /(R 1×(p-r) i
F =      1 - ∂y ∂x ∂y 2 +∂z 2 - ∂x (∂y 2 λ+2 ∂z 2 λ+∂x 2 λ+3 ∂z 2 µ+2 ∂y 2 µ+2 ∂x 2 µ) (λ+µ)(∂y 2 +∂z 2 )∂z 0 -∂ z ∂ y 0 0 1      ∈ GL 3 (R 1 ).
However, if we choose the matrix F 2 defined by

F 2 = 2 ∂ x µ + ∂ x λ ∂ y λ + 2 ∂ y µ ∂ z λ + 2 ∂ z µ 1 0 0 ,
then no completion exists for such F 2 .

Reduction of interface conditions

In the algorithms presented in Section 3, we have both equations in the domains Ω i and at the interface Γ. In Sections 4 and 5, we saw how to construct many possible Smith variables and how to use them to get simpler equations for the corresponding decomposition domain algorithm. But, to simplify and speed up the algorithm, we also need to reduce the interface conditions with respect to the equations in the domains. In the present section, we show how symbolic computation techniques can be used to perform such reductions. A natural idea consists in gathering all equations and computing a Gröbner basis (see Section 7). However, one has to keep in mind that the independent variables do not play the same role. More precisely, the interface conditions cannot be differentiated with respect to x 1 since the border of the interface is defined by the equation x 1 = 0.

Consequently, we need to develop an alternative approach for reducing interface conditions. This method has been implemented in Maple using the OreModules package ( [START_REF] Chyzak | OreModules: A symbolic package for the study of multidimensional linear systems[END_REF]). See Section 8.2 for many examples of computations. The approach can be sketched as follows:

1. Compute a Gröbner basis of the equations inside the domain for a relevant monomial order.

2. Compute the normal forms of the interface conditions w.r.t. the latter Gröbner basis.

Write this normal forms in jet variables

in x 1 .
4. Perform linear algebra manipulations to simplify the normal forms. See Section 7 for precise definitions of these concepts.

Without loss of generality and in order to simplify the notations, we shall assume that the equations only involve three dependent variables u, v, w in two independent variables x 1 and x 2 . We then set

R = k[∂ 1 , ∂ 2 ],
where k is a field, and

R 1 = k(∂ 2 )[∂ 1 ]. Equations inside the domains Inria have the form p 1 (∂ 1 , ∂ 2 ) u + p 2 (∂ 1 , ∂ 2 ) v + p 3 (∂ 1 , ∂ 2 ) w = 0,
where the p i 's are polynomials with coefficients in k. We consider the set of monomials in the four variables ∂ x , u, v, w, and endow it with the graded reverse lexicographic order defined in 3 of Example 7 below such that 1 ≺ degrevlex w ≺ degrevlex v ≺ degrevlex u ≺ degrevlex ∂ x . Then, we compute a Gröbner basis (see Section 7) of the equations inside the domain with respect to the latter admissible term order (note that ∂ 2 is then viewed as a parameter). Now, the interface conditions (computed from the Smith variables) are given by equations of the form q 1 (∂ 1 ) u + q 2 (∂ 1 ) v + q 3 (∂ 1 ) w = 0, where the q i 's are polynomials with coefficients in k(∂ 2 ). We can then compute their normal forms with respect to the Gröbner basis. Once this has been done, to keep on simplifying the interface conditions, we consider them all together and try to perform linear algebra simplifications between them. Here again, we must be careful since multiplications by ∂ 1 are not allowed. To avoid such multiplications, we first rewrite the normal forms of the interface conditions in jet variables in x 1 , i.e., we replace the successive derivatives with respect to x 1 by new variables. For instance, ∂ 1 u, ∂ 1 v, ∂ 1 w are replaced by the new variables u 1 , v 1 , w 1 . We can then perform linear algebra simplifications between the normal forms of the interface conditions to get simplified form of the interface conditions. This method is implemented in the Maple package Schwarz and is illustrated step by step on the example of the 2D linear elasticity system given in Section 8.2.1 below.

Appendix: An introduction to Gröbner basis techniques

For polynomial rings, Janet ( [START_REF] Janet | Leçons sur les systèmes d'équations aux dérivées partielles[END_REF]) and Buchberger ([3]) developed two constructive algorithms which compute new sets of generators of a (left) ideal or a (left) module, called a Janet or a Gröbner basis. Algorithms rewriting any element of that (left) ideal or (left) module in terms of the new generators were also obtained. More generally, normal forms of general elements can be computed with respect to the Janet or Gröbner bases (see, e.g., [START_REF] Buchberger | An algorithm for finding the basis elements in the residue class ring modulo a zero dimensional polynomial ideal[END_REF][START_REF] Cox | Ideals, Varieties, and Algorithms[END_REF][START_REF] Cox | Using Algebraic Geometry[END_REF][START_REF] Greuel | A Singular Introduction to Commutative Algebra[END_REF] and the references therein). More precisely, if R denotes a commutative polynomial ring over a computable field k (e.g., Q, Z/Z p, where p is a prime), and A ∈ R q×p , then the knowledge of a Janet or a Gröbner basis of the R-submodule R 1×q A of R 1×p allows one to compute the normal form of any element λ ∈ R 1×p with respect to the Janet or Gröbner basis, i.e., to compute a distinguished representative of the residue class of λ in the quotient R-module M = R 1×p /(R 1×q A). In particular, Janet or Gröbner basis techniques allow one to constructively work in the (left) ideals and (left) modules (e.g., computation of kernels, images, factorizations, left/right/generalized inverses of multivariate polynomial matrices). They are nowadays implemented in different computer algebra systems such as, e.g., Maple, Mathematica, Singular, Macaulay2, CoCoA, . . . for different classes of commutative and noncommutative polynomial rings.

We shortly introduce the basic definitions and results on Gröbner bases for ideals and modules over commutative polynomial rings. For more details, see [START_REF] Buchberger | An algorithm for finding the basis elements in the residue class ring modulo a zero dimensional polynomial ideal[END_REF][START_REF] Cox | Ideals, Varieties, and Algorithms[END_REF][START_REF] Cox | Using Algebraic Geometry[END_REF][START_REF] Greuel | A Singular Introduction to Commutative Algebra[END_REF][START_REF] Kandri-Rody | Non-commutative Gröbner bases in algebras of solvable type[END_REF].

In what follows, we explain how to compute a Gröbner basis in the case of a commutative polynomial ring

R = k[x 1 , . . . x n ],
where k is a computable field (e.g., Q, Z/Z p, where p is a prime). Every element P of R can uniquely be written in the form

P = |µ|=0,...,r a µ x µ , a µ ∈ k, x µ = x µ1 1 • • • x µn n ,
where µ = (µ 1 . . . µ n ) ∈ N 1×n . The set defined by Mon(R) = {x µ | µ ∈ N n } is a basis of the k-vector space R. The elements of Mon(R) are called monomials of R.

In order to effectively study systems over commutative polynomial rings, we first need to introduce monomial orders to compare polynomials (see, e.g., [START_REF] Buchberger | An algorithm for finding the basis elements in the residue class ring modulo a zero dimensional polynomial ideal[END_REF][START_REF] Cox | Ideals, Varieties, and Algorithms[END_REF][START_REF] Greuel | A Singular Introduction to Commutative Algebra[END_REF][START_REF] Kandri-Rody | Non-commutative Gröbner bases in algebras of solvable type[END_REF]). Definition 4.

1. A total order ≺ on the set Mon(R) is called a well-ordering if every nonempty subset of Mon(R) has a least element with respect to ≺.

2. A total well-ordering order ≺ on the set Mon(R) is called an admissible term order of Mon(R) if it satisfies the following two conditions:

(a) 1 is the least element of Mon(R), namely, for all u ∈ Mon(R), u = 1 ⇒ 1 ≺ u.

(b) ≺ is compatible with the product, namely, if u, v ∈ Mon(R) satisfy u ≺ v, then w u ≺ w v for all w ∈ Mon(R).

3. Given a nonzero polynomial P ∈ R = k[x 1 , . . . , x n ] and an admissible term order ≺ on Mon(R), we can compare the nonzero terms of P with respect to ≺. The greatest of these monomials is called the leading monomial of P and is denoted by lm(P ). The coefficient of lm(P ) is the leading coefficient of P , denoted by lc(P ), and the leading term lc(P ) lm(P ) of P denoted by lt(P ).

Let us give important examples of admissible term orders of Mon(R).

Example 7. We can identify an element x µ of Mon(R) with the multi-index µ ∈ N 1×n .

1. The pure lexicographical order on Mon(R) is defined by µ ≺ plex ν whenever the first nonzero entry of νµ is positive. In particular, it means that

x 3 ≺ plex x 2 ≺ plex x 1 . More generally, if we consider R = Q[x 1 , x 2 , x 3
], then we have:

1 ≺ plex x 3 ≺ plex x 2 3 ≺ plex x 2 ≺ plex x 2 x 3 ≺ plex x 2 2 ≺ plex x 1 ≺ plex x 1 x 3 ≺ plex x 1 x 2 ≺ plex x 2 1 .
2. The graded lexicographic order on Mon(R) is defined by µ ≺ grlex ν whenever |µ| < |ν| or if we have |µ| = |ν|, then the first nonzero entry of νµ is positive. For instance, if we consider R = Q[x 1 , x 2 , x 3 ], then we have:

1 ≺ grlex x 3 ≺ grlex x 2 ≺ grlex x 1 ≺ grlex x 2 3 ≺ grlex x 2 x 3 ≺ grlex x 2 2 ≺ grlex x 1 x 3 ≺ grlex x 1 x 2 ≺ grlex x 2 1 .
3. The graded reverse lexicographic order on Mon(R), also called degree reverse lexicographical order, is defined by µ ≺ degrevlex ν whenever |µ| < |ν| or if we have |µ| = |ν|, then the last nonzero entry of νµ is negative. It is also denoted by ≺ tdeg . For instance, if we consider R = Q[x 1 , x 2 , x 3 ], then we have:

1 ≺ degrevlex x 3 ≺ degrevlex x 2 ≺ degrevlex x 1 ≺ degrevlex x 2 3 ≺ degrevlex x 2 x 3 ≺ degrevlex x 1 x 3 ≺ degrevlex x 2 2 ≺ degrevlex x 1 x 2 ≺ degrevlex x 2 1 .
4. Let R be the polynomial ring over k with indeterminates x 1 , . . . , x n , y 1 , . . . , y m . Assume that admissible monomial orders ≺ x and ≺ y on the monomials that only contain respectively the x i 's and the y i 's are given. An elimination order is then defined by

u v ≺ w t ⇔ u ≺ x w or u = w and v ≺ y t,
where u, w (resp., v, t) are monomials containing only the x i 's (resp., y i 's). An elimination order serves to eliminate the x i 's. The elimination order which we shaill use in what follows is the one induced by the degree reverse lexicographical orders on x 1 , . . . , x n and y 1 , . . . , y m . This is a very common order called lexdeg. For instance, if we consider R = Q[x 1 , x 2 , x 3 ],

x = (x 1 , x 2 ), y = x 3 , ≺ x =≺ tdeg and ≺ y =≺ tdeg , then we have:

1 ≺ lexdeg x 3 ≺ lexdeg x 2 3 ≺ lexdeg x 2 ≺ lexdeg x 2 x 3 ≺ lexdeg x 1 ≺ lexdeg x 1 x 3 ≺ lexdeg x 2 2 ≺ lexdeg x 1 x 2 ≺ lexdeg x 2 1 .
Inria Definition 5. We call monomial ideal of R = k[x 1 , . . . , x n ] an ideal of R generated by monomials.

If J is a subset of N 1×n , then x J denotes the set of monomials {x µ | µ ∈ J} and (x J ) the ideal of R = k[x 1 , . . . , x n ] generated by the elements of x J .

Example 8. Let us consider the commutative polynomial ring R = Q[x 1 , x 2 ] and the subset J = {(2, 0), (0, 2)} of N 1×2 . Then, the monomial ideal (x 2 1 , x 2 2 ) of R generated by x 2 1 and x 2 2 is defined by elements of the form P 1 x 2 1 + P 2 x 2 2 , where P 1 and P 2 are two arbitrary polynomials of R. In particular, the monomials of polynomial multiples of x 2 1 correspond to the integer points of the translated first quadrant at the point (2, 0). Similarly, the monomials of polynomial multiples of x 2 2 correspond to the integer points of the translated first quadrant at the point (0, 2). Hence, the monomials of any element of (x 2 1 , x 2 2 ) belong to those two translated first quadrants.

The following lemma is called Dickson's lemma.

Lemma 3. Every monomial ideal (x J ) of R = k[x 1 , . . . ,
x n ] is generated by a finite set of monomials of x J .

A proof of Lemma 3 can be found in any textbooks on Gröbner bases (see, e.g., [START_REF] Cox | Ideals, Varieties, and Algorithms[END_REF][START_REF] Greuel | A Singular Introduction to Commutative Algebra[END_REF]).

If we now consider a non-empty subset J of N 1×n and the corresponding non-empty set x J of monomials, then according to Lemma 3, the monomial ideal (x J ) is generated by a finite number of elements of x J , say L = {x α , x β , . . . , x θ }. Let us now consider a total order ≺ on the set Mon(R) satisfying a and b of 2 of Definition 4, and µ ∈ J. Then, x µ ∈ (x α , x β , . . . , x θ ) and thus there exists x γ ∈ L and ν ∈ N n such that x µ = x ν x γ . Hence, we either have x ν = 1 which yields x µ = x γ ∈ L or x ν = 1 which, according to 1 ≺ x ν and b of 2 of Definition 4, implies x γ ≺ x γ x ν = x µ . In other words, we have x µ ∈ L or x µ is greater than an element of L, which shows that any element x µ of x J is greater or equal to the least element of L, which finally proves that a total order on the set Mon(R) satisfying a and b of 2 of Definition 4 is a well-ordering order, and thus an admissible term order.

We can now introduce the concept of a Gröbner basis which plays an important role for the computational issues in mathematical systems theory. Definition 6. Let R = k[x 1 , . . . , x n ] be a commutative polynomial ring over a computable field k, ≺ an admissible term order on Mon(R), and I an ideal of R. A set of nonzero polynomials G = {Q i } i=1,...,t ⊂ I is called a Gröbner basis for I if for all nonzero element P in I, there exists i ∈ {1, . . . , t} such that lm(Q i ) divides lm(P ), i.e., lm(P ) = M lm(Q i ), for a certain M ∈ Mon(R).

If we denote by lm(I) the monomial ideal of R = k[x 1 , . . . , x n ] generated by the leading terms of the elements of I, i.e., lm(I) = (lm(P )) P ∈I\{0} , then a family G = {Q i } i=1,...,t of elements of I is a Gröbner basis of I iff lm(I) = (lm(Q 1 ), . . . , lm(Q t )). We note that the key point in Definition 6 is that G must be finite because, otherwise, we can always take G = I \ {0}.

The existence of a Gröbner basis of I is a straightforward consequence of Lemma 3. Indeed, according to Lemma 3, there exists a finite set {Q i } i=1,...,t of polynomials of I such that lm(I) = (lm(Q 1 ), . . . , lm(Q t )). Hence, an interesting issue is to explicitly compute those elements Q i 's from a given set of generators of I and an admissible term order ≺ on Mon(R).

Before explaining how to compute Gröbner bases, we first want to explain their most important property. To do that, we first need to introduce the concept of reduction. In what follows, we fix an admissible term order ≺ on Mon(R). Let P, Q ∈ R be given. If lm(Q) divides lm(P ), then P can be reduced modulo Q, namely, the leading term of P can be eliminated by considering the new polynomial S defined by

S = P - lc(P ) lm(P ) lc(Q) lm(Q) Q = P - lt(P ) lt(Q) Q,
which satisfies lm(S) ≺ lm(P ). Let us illustrate the concept of reduction.

Example 9. Let R = Q[x 1 ,
x 2 ] be the polynomial ring with the degree reverse lexicographical order ≺ degrevlex and the two polynomials

P = 2 x 2 1 + 3 x 1 x 7 2 and Q = 6 x 3 2 -x 1 .
In particular, we have 1

≺ x 2 ≺ x 1 ≺ x 2 2 ≺ x 1 x 2 ≺ x 2 1 ≺ x 3 2 ≺ x 1 x 2 2 ≺ x 2 1 x 2 ≺ x 3 1 ≺ . . .,
where, for simplicity reasons, we have used the notation ≺ for ≺ degrevlex . Then, we have lm(

P ) = x 1 x 7 2 , lm(Q) = x 3 2
and lm(Q) divides lm(P ), which shows that P can be reduced modulo Q as follows:

S = P - 3 x 1 x 7 2 6 x 3 2 Q = P - 1 2 x 1 x 4 2 Q = 1 2 x 2 1 x 4 2 + 2 x 2 1 .
Now, we can reduce P modulo a finite set {Q i } i=1,...,t of polynomials as follows.

Algorithm 7.

• Input:

P ∈ R = k[x 1 , . . . , x n ], a finite set {Q i } i=1,.
..,t of elements of R, and admissible term order ≺ on Mon(R).

• Output: Elements d 1 , . . . , d t and S of R such that

P = t i=1 d i Q i + S,
and none of the monomials of S can be reduced by an element of {Q i } i=1,...,t .

1. Set S = 0 and d i = 0 for i = 1, . . . , t.

While P = 0, do:

• If there exists i ∈ {1, . . . , t} such that lm(Q i ) divides lm(P ), then do:

d i ←-d i + lt(P ) lt(Q i ) , P ←-P - lt(P ) lt(Q i ) Q i ,
• Else, do S ←-S + lt(P ) and P ←-Plt(P ).

3. Return d 1 , . . . , d t and S.

Remark 5. We point out that the output of Algorithm 7 depends on the particular choices of i ∈ {1, . . . , t} such that lm(Q i ) divides lm(P ).

Algorithm 7 is usually called the division algorithm in R and S is the remainder of P on division by {Q i } i=1,...,t for the admissible term order ≺ or simply, the remainder of P modulo {Q i } i=1,...,t . In particular, it generalizes the Euclidean division for a multivariate commutative polynomial ring

R = k[x 1 , . . . , x n ].
An important property of a Gröbner basis G = {Q i } i=1,...,t of an ideal I of R is that each polynomial P ∈ I is reduced to 0 modulo G, i.e., by subtracting suitable multiples of the Q i from P , we get the zero polynomial. Indeed, applying Algorithm 7 to P and to the Gröbner basis G = {Q i } i=1,...,t of I, we obtain polynomials d 1 , . . . , d t and S such that

P = t i=1 d i Q i + S,

Inria

and none of the monomials of S can be reduced by an element of {Q i } i=1,...,t . But, S = P -t i=1 d i Q i ∈ I, and if S = 0, then by definition of G, lm(S) ∈ lm(I) = (lm(Q 1 ), . . . , lm(Q t )), i.e., there exist b 1 , . . . , b t ∈ R such that lm(S) = t i=1 b i lm(Q i ), which yields that there exists j ∈ 1, . . . , t such that lm(Q j ) divides lm(S), which contradicts the hypothesis on S, and thus S = 0. Hence, every element P ∈ I can be written as P = t i=1 d i Q i , i.e., P ∈ (Q 1 , . . . , Q t ), and thus I = (Q 1 , . . . , Q t ) because the Q i 's belong to I, which shows that a Gröbner basis G = {Q i } i=1,...,t of an ideal I defines a set of generators of I, i.e., I = (Q 1 , . . . , Q t ).

Proposition 1. Given a Gröbner basis G of an ideal I of R = k[x 1 , . . . , x n ] for an admissible term order ≺, P ∈ I iff the remainder of P modulo G is zero.

The knowledge of a Gröbner basis G of I allows us to solve the membership problem, namely, to check whether or not a given polynomial P belongs to I: P belongs to I iff the remainder S returned in the output of Algorithm 7 for a Gröbner basis G = {Q i } i=1,...,t of I is zero.

Let us now shortly explain how to compute a Gröbner basis of an ideal. For more details, we refer to [START_REF] Buchberger | An algorithm for finding the basis elements in the residue class ring modulo a zero dimensional polynomial ideal[END_REF][START_REF] Cox | Ideals, Varieties, and Algorithms[END_REF][START_REF] Greuel | A Singular Introduction to Commutative Algebra[END_REF][START_REF] Kandri-Rody | Non-commutative Gröbner bases in algebras of solvable type[END_REF]. Given P, Q ∈ I, we can let L = lcm(lm(P ), lm(Q)) to be the least common multiple of the leading monomials of P and Q. One can obviously find another element of the ideal I by computing the following so-called S-polynomial of P and Q:

S(P, Q) = L lt(P ) P - L lt(Q) Q.
Example 10. We consider the polynomial ring R = Q[x 1 , x 2 , x 3 ] with the degree reverse lexicographical order ≺ degrevlex and denote by I the ideal of R generated by

P = 2 x 1 x 3 -x 3 and Q = x 2 1 x 2 -1. Then, we have lm(P ) = x 1 x 3 , lc(P ) = 2, lm(Q) = x 2 1
x 2 , and lc(Q) = 1. Hence, the S-polynomial S(P, Q) is defined by:

S(P, Q) = x 2 1 x 2 x 3 2 x 1 x 3 P - x 2 1 x 2 x 3 x 2 1 x 2 Q = 1 2 x 1 x 2 P -x 3 Q = - 1 2 x 1 x 2 x 3 + x 3 .
Now, since lm(S(P, Q)) = x 1 x 2 x 3 is divisible by lm(P ), we obtain that S(P, Q) can be reduced modulo P as follows:

T = S(P, Q) - -1 2 x 1 x 2 x 3 2 x 1 x 3 P = S(P, Q) -- 1 4 x 2 P = x 3 - 1 4 x 2 x 3 .
Of course, by construction, both S(P, Q) and T belong to the ideal I.

As previously noticed, S-polynomials of elements of I belong to I.

Theorem 4. Let I be an ideal of R = k[x 1 , . . . , x n ] and an admissible term order ≺ on Mon(R). Then, G = {Q i } i=1,...,t is a Gröbner basis of I for ≺ iff for all pairs

(Q i , Q j ) of distinct elements of G, the remainder of S(Q i , Q j ) modulo G is zero.
According to Theorem 4, Buchberger's algorithm ( [START_REF] Buchberger | An algorithm for finding the basis elements in the residue class ring modulo a zero dimensional polynomial ideal[END_REF]) takes a generating set {P j } j=1,...,r of the ideal I as input and constructs a Gröbner basis of I by starting with G = {P j } j=1,...,r and computing the S-polynomials of pairs of distinct polynomials in G. The S-polynomials are reduced as long as possible modulo polynomials in G. The result of this sequence of reductions are called the remainders of the S-polynomials modulo G. Every remainder different from zero is then added to the set G and this procedure is iterated as long as there exist nonzero remainders of S-polynomials. This algorithm terminates with a Gröbner basis G of the ideal I for the admissible term order ≺.

Example 12. Given an ideal I of R = k[x 1 , . . . , x n , y 1 , . . . , y m ], we obtain a Gröbner basis of the left ideal I ∩ k[y 1 , . . . , y m ] by computing the Gröbner basis G of I with respect to an elimination order (see 4 of Example 7) and intersecting G with k[y 1 , . . . , y m ], which merely amounts to omitting all polynomials in G that involve any x i .

Finally, we shortly explain how we can extend the previous results from ideals to modules over a commutative polynomial algebra R (see, e.g., [START_REF] Cox | Using Algebraic Geometry[END_REF][START_REF] Greuel | A Singular Introduction to Commutative Algebra[END_REF]).

Let us denote by {f j } j=1,...,p the standard basis of the following free left R-module

R 1×p = {(λ 1 . . . λ p ) | λ i ∈ R, i = 1, . . . , p},
namely, the k th component of f j is 1 if k = j and 0 otherwise. First, we need to extend the term order ≺ on Mon(R) to the set of monomials of the form u f j , where u ∈ Mon(R) and j = 1, . . . , p. When it does not lead to any confusion, we still denote by ≺ the extension of ≺ to Mon(R 1×p ) = p j=1 Mon(R) f j , where ≺ is asked to satisfy the following two conditions:

1. ∀ w ∈ Mon(R) : u f i ≺ v f j ⇒ w u f i ≺ w v f j . 2. ∀ j = 1, . . . , p : u ≺ v ⇒ u f j ≺ v f j .
Without loss of generality, we let

f p ≺ f p-1 ≺ • • • ≺ f 1 .
Then, there are two natural extensions of a term order to Mon(R 1×p ) defined in the next definition. Definition 8. Let ≺ be an admissible term order on Mon(R), u, v ∈ Mon(R), and {f j } j=1,...,p the standard basis of the R-module R 1×p .

1. The term over position order on Mon(R 1×p ) induced by ≺ is defined by:

u f i ≺ v f j ⇔ u ≺ v or u = v and f i ≺ f j .
2. The position over term order on Mon(R 1×p ) induced by ≺ is defined by:

u f i ≺ v f j ⇔ f i ≺ f j or f i = f j and u ≺ v.
The term over position order is of more computational value with regard to efficiency, whereas the position over term order can be used to eliminate components: Buchberger's algorithm applied to a generating set of a R-submodule L of R 1×p with respect to a position over term order finds (if it exists) an element of the form P f j ∈ L, where P ∈ R.

If an admissible monomial order on Mon(R 1×p ) is fixed, then leading monomials, leading coefficients, leading terms of elements in R 1×p are defined as in the case of ideals. Moreover, Buchberger's algorithm carries over to R 1×p . For more details, we refer, e.g., to [START_REF] Cox | Using Algebraic Geometry[END_REF][START_REF] Greuel | A Singular Introduction to Commutative Algebra[END_REF]. We only give here an example.

Example 13. We endow the R = Q[x 1 , x 2 ]-module R 1×2 with the term over position order induced by ≺ degrevlex on Mon(R) and denote by L the R-submodule of R 1×2 generated by

P = (x 2 1 + x 2 ) f 1 + x 1 x 2 f 2 and Q = (x 2 + 1) f 1 + f 2 . The S-polynomial S(P, Q), defined by S(P, Q) = x 2 P -x 2 1 Q = (x 2 2 -x 2 1 ) f 1 + (x 1 x 2 2 -x 2 1
) f 2 , which can be reduced by P and Q to:

T = S(P, Q) + P -x 2 Q = (x 1 x 2 2 -x 2 1 + x 1 x 2 -x 2 ) f 2 .
Then, P can still be replaced as a generator of L by P -Q = (x 2 1 -1) f 1 + (x 1 x 2 -1) f 2 , but no further nonzero S-polynomials can be computed from P -Q, Q and T , since the leading term of T lies in the second component, whereas the leading terms of P -Q and Q lie in the first one. Thus, Buchberger's algorithm terminates with a Gröbner basis

G = {P -Q, Q, T } of L.
Thereafter, non-experts on Gröbner basis techniques can only consider a Gröbner basis as a kind of "black box" which allows one to compute much important information on polynomial systems. For instance, using them, we can constructively solve the following problems:

1. Computation of kernels: Given a matrix A ∈ R q×p , find B ∈ R r×q satisfying:

ker R (.A) := {λ = (λ 1 . . . λ q ) ∈ R 1×q | λ A = 0} = R 1×r B := {µ B | µ ∈ R 1×r }.
This is done by the following algorithm: Algorithm 9.

• Input: A commutative polynomial ring R and a finitely generated R-submodule L of R 1×p defined by a matrix A ∈ R q×p , i.e., L = R 1×q A.

• Output: A matrix B ∈ R r×q such that ker R (.A) = R 1×r B.
(a) Introduce the indeterminates η 1 , . . . , η p , ζ 1 , . . . , ζ q over R and define the following set:

P =    p j=1 A ij η j -ζ i | i = 1, . . . , q    .
(b) Compute a Gröbner basis G of P in the free R-module generated by the η j 's and the ζ i 's for j = 1, . . . , p and i = 1, . . . , q, namely,

p j=1 R η j ⊕ q i=1 R ζ i ,
with respect to a term order which eliminates the η j 's.

(c) Compute the intersection G ∩ ( q i=1 R ζ i ) = { q i=1 B ki ζ i | k = 1, . . . , r} by selecting the elements of G containing only the ζ i 's, and return B = (B ij ) ∈ R r×q .
2. Computation of factorizations: Given two matrices A ∈ R q×p and A ′ ∈ R q ′ ×p , find a matrix A ′′ ∈ R q×q ′ , when it exists, satisfying

A = A ′′ A ′ .
This can be done by the following algorithm:

Algorithm 10.

• Input: A commutative polynomial ring R, two matrices A ∈ R q×p , and A ′ ∈ R q ′ ×p .

• Output: A matrix A ′′ ∈ R q×q ′ such that A = A ′′ A ′ if A ′′ exists and [ ] otherwise.
(a) Introduce the indeterminates λ j , j = 1, . . . , p and µ i , i = 1, . . . , q ′ , over R, and define:

P =    p j=1 A ′ ij λ j -µ i | i = 1, . . . , q ′    . (b) Compute the Gröbner basis G of P in p j=1 R λ j ⊕ q ′
i=1 R µ i with respect to a term order which eliminates the λ j 's.

(c) Define the following set:

Q =    p j=1 A kj λ j | k = 1, . . . , q    . (d) Reduce each element of Q w.r.t. G (see Algorithm 7) and call H = {H i | i = 1, . . . , q}
the set obtained.

Inria

From the two above simple examples, it clearly appears that the corresponding computations cannot generally be obtained without the use of Gröbner basis techniques and without their implementations in a computer algebra system. If an implementation of Buchberger's algorithm is only at hand for R and not for R 1×p , where p > 1, then one can overcome this problem by introducing new variables f 1 , . . . , f p and ignoring any polynomial which contains any nonlinear term in the f j 's in the result (we note that the input for Buchberger's algorithm consists of polynomials which are homogeneous in the f j 's and reduction keeps homogeneity). In particular, the Maple package Groebner actually does not form any product with the f j 's when the monomial order is properly defined. This fact is used in the development of the OreModules package ( [START_REF] Chyzak | OreModules: A symbolic package for the study of multidimensional linear systems[END_REF]) built upon the Maple packages Groebner and Ore -algebra ( [START_REF] Chyzak | Non-commutative elimination in Ore algebras proves multivariate identities[END_REF]). In particular, OreModules allows one to effectively study the previously listed problems (and more) for the class of Ore algebras of functional operators available in Ore -algebra. The OreModules package is a main computational tool for studying the properties of finitely presented modules over Ore algebras appearing in mathematical systems theory. The package Janet ([2]) follows the same philosophy but is based on the concept of Janet bases ([2, [START_REF] Janet | Leçons sur les systèmes d'équations aux dérivées partielles[END_REF]31]).

Appendix: Maple computations 8.1 Completion problem

In this section, we study the completion problem considered in Section 5.2.

Elasticity 3D

> restart: > with(linalg): > with(OreModules): > with(Schwarz):

We consider the elastostatic equations, i.e., the Navier-Cauchy equations in R 3 . This system is defined by R y = 0, where R ∈ A 3×3 is a matrix with entries in the commutative polynomial ring A = Q(λ, µ)[dx, dy, dz] of PD operators in dx = ∂/∂x, dy = ∂/∂y and dz = ∂/∂z with coefficients in Q(λ, µ), where λ and µ are the two Lamé constants.

We first define A and R.

R :=     -2 dx 2 µ -dx 2 λ -dy 2 µ -dz 2 µ -dx dy (λ + µ) -dx dz (λ + µ) -dx dy (λ + µ) -dx 2 µ -2 dy 2 µ -dy 2 λ -dz 2 µ -dy dz (λ + µ) -dx dz (λ + µ) -dy dz (λ + µ) -dx 2 µ -dy 2 µ -2 dz 2 µ -dz 2 λ     Inria Let ∆ = dx 2 + dy 2 + dz 2 .
We can check that the Smith normal form of R is the diagonal matrix having 1, ∆ and ∆ 2 as diagonal entries > S:=map(factor,smith(R,dx));

S :=     1 0 0 0 dx 2 + dy 2 + dz 2 0 0 0 dx 2 + dy 2 + dz 2 2     i.e., R = E S F , where E, F ∈ GL 3 (B) and B = Q(λ, µ, dy, dz)[dx].
As a consequence, the residue class of the second (resp., third) row of the matrix F in the B-module B 1×3 /(B 1×3 R) must be annihilated by the PD operator ∆ (resp., ∆ 2 ). To find a set of possible F 's (i.e., possible "Smith variables"), we first compute families of generators of the elements of the Amodule M = A 1×3 /(A 1×3 R) respectively annihilated by the PD operators ∆ and ∆ 2 . This can be done by means of the AnnOp procedure.

> Delta:=dx^2+dy^2+dz^2: > F2:=AnnOp(Delta,R,A);

F 2 :=        -dz µ 0 dx µ -dy µ dx µ 0 2 dx µ + dx λ dy λ + 2 dy µ dz λ + 2 dz µ 0 -dz µ dy µ        > F3:=AnnOp(Delta^2,R,A); F 3 :=     2 µ 2 + λ µ 0 0 0 2 µ 2 + λ µ 0 0 0 2 µ 2 + λ µ    
Then, the second (resp., third) row of F can be chosen among the B-linear combinations of the rows of F 2 (resp., F 3). Once these two rows have been chosen, we can try to complement them to a unimodular matrix F which further satisfies that its first row (the one that we need to find) is annihilated by 1, i.e., by the first entry of the diagonal matrix S. As it is illustrated below, this complement does not always exist since the choices made for the second and third row of F may not lead to a unimodular matrix whose first row is annihilated by 1.

The procedure SmithVariablesCompletion takes as inputs the matrix R, the lower part of the Smith form (here diag(∆, ∆ 2 )), the chosen rows of F (here a row of F 2 and a row of F 3) and B. It decides whether or not the given choice can be complemented as previously explained. If so, then it returns the corresponding matrices E, S and F , and [ ] otherwise. For instance, if we choose the fourth row of F 2 and the third row of F 3 (up to constant) > F23:=stackmatrix(row(F2,4),row(F3,3)/(2*mu^2+lambda*mu));

F 23 := 0 -dz µ dy µ 0 0 1
and run the procedure SmithVariablesCompletion:

RR n°7953 > B:=DefineOreAlgebra(diff=[dx,x],polynom=[x]
,comm=[dy,dz,lambda,mu]): > SVC:=SmithVariablesCompletion(R,diag(Delta,Delta^2),F23,B):

> E:=SVC[1]; > S:=SVC[2]; > F:=SVC[3]; E :=       (-2 µ -λ) dx 2 -dy 2 + dz 2 µ (2 µ+λ)dy dx (dy 2 +dz 2 )dz µ (-4 µ 2 -4 λ µ-λ 2 )dx dz (dy 2 λ+dz 2 µ+dz 2 λ+dy 2 µ)
-dx dy (λ + µ)

dz 2 µ+2 dy 2 µ+dy 2 λ (dy 2 +dz 2 )dz µ -(2 µ+λ)dy (dy 2 +dz 2 )dz -dx dz (λ + µ) (λ+µ)dy (dy 2 +dz 2 )µ -2 µ-λ dy 2 +dz 2       S :=     1 0 0 0 dx 2 + dy 2 + dz 2 0 0 0 dx 2 + dy 2 + dz 2 2     F :=      1 -dy dx dy 2 +dz 2 - dx (dy 2 λ+2 dz 2 λ+dx 2 λ+3 dz 2 µ+2 dy 2 µ+2 dx 2 µ) (λ+µ)(dy 2 +dz 2 )dz 0 -dz µ dy µ 0 0 1     
We find that this particular choice can be complemented to a unimodular matrix F whose first row is annihilated by 1. However, choosing the third row of F 2 and the first row of F 3, we obtain > F23p:=stackmatrix(row(F2,3),row(F3,1)/(2*mu^2+lambda*mu));

F 23p := 2 dx µ + dx λ dy λ + 2 dy µ dz λ + 2 dz µ 1 0 0
and running again the procedure, we find that this choice cannot be complemented, which yields an empty output:

> SVC:=SmithVariablesCompletion(R,diag(Delta,Delta^2),F23p,B); []
Once a "good choice" has been made, we can run the procedure ReducedInterfaceConditions (see Section 8.2) to reduce the interface conditions obtained with the corresponding F . For instance, with the "good choice" previously given, we can reduce the interface conditions of the update step of the algorithm defined by means of the PD operators 1, 1 and ∆, i.e., apply ReducedInterfaceConditions to U , V and ∆ W , where U, V, W denote the system variables.

> ReducedInterfaceConditions(R,SVC[3],A,[1,1,Delta],[U,V,W]); [V = 0, W = 0, U x = 0]
To obtain a longer list of reduced interface conditions (e.g., for the interface conditions defined by means of 1, 1 and ∆ for the update step, and by dx, dx, dx ∆ for the correction step) computed with distinct matrices F constructed as previously explained, we loop over all the possible combinations made with one row of F 2 and one row of F 3, check whether or not this candidate is "good" (i.e., whether or not the corresponding matrix can be complemented as previously explained), and if so, compute the corresponding F , and then run the procedure ReducedInterfaceConditions.

Inria

> for i from 1 to rowdim(F2) do > for j from 1 to rowdim(F3) do > F23:=stackmatrix(row(F2,i),row(F3,j)/(2*mu^2+lambda*mu)): > SVC:=SmithVariablesCompletion(R,diag(Delta,Delta^2),F23,B):

> if SVC=[] then > print([[i,j],[]]); > else > print([[i,j], > ReducedInterfaceConditions(R,SVC[3],A,[1,1,Delta],[U,V,W]), > ReducedInterfaceConditions(R,SVC[3],A,[dx,dx,dx*Delta],[U,V,W])]); > fi: > od: > od: [[1, 1], []] [[1, 2], [V = 0, U x = -W dz , W x = dz U ], [U x = -dy dz V λ+dy dz V µ+W dy 2 µ+2 W dz 2 µ+W dz 2 λ dz (2 µ+λ) , V x = 0, W x = (dy 2 +dz 2 )U dz ]] [[1, 3], [W = 0, U x = -dy V, W x = dz U ], [U x = -dy dz V λ+dy dz V µ+W dy 2 µ+2 W dz 2 µ+W dz 2 λ dz (2 µ+λ) , V x = (dy 2 +dz 2 )U dy , W x = 0]] [[2, 1], []] [[2, 2], [V = 0, U x = -W dz , V x = U dy], [U x = -2 V dy 2 µ+V dy 2 λ+V dz 2 µ+dy dz W λ+dy dz W µ (2 µ+λ)dy , V x = 0, W x = (dy 2 +dz 2 )U dz ]] [[2, 3], [W = 0, U x = -dy V, V x = U dy], [U x = -2 V dy 2 µ+V dy 2 λ+V dz 2 µ+dy dz W λ+dy dz W µ (2 µ+λ)dy , V x = (dy 2 +dz 2 )U dy , W x = 0]] [[3, 1], []] [[3, 2], []] [[3, 3], []] [[4, 1], [U = 0, W = dz V dy , V x = -Wxdz dy ], [V = -W dz dy , U x = 0, W x = dz Vx dy ]] [[4, 2], [V = 0, W = 0, U x = 0], [U = 0, V x = 0, W x = 0]] [[4, 3], [V = 0, W = 0, U x = 0], [U = 0, V x = 0, W x = 0]]
With this method, we can find the Smith variables (i.e., the F 's) that lead to simple reduced interface conditions (for both the update and the correction steps). From the above computations, the simplest reduced interface conditions were obtained for the matrix F formed by the fourth row of F 2 and the second (or third) row of F 3.

We note that we do not test all the possible unimodular matrices F such that R = E S F , where E unimodular, since the second (resp., third) row of F can be chosen as a B-linear combination of the rows of F 2 (resp., F 3). Consequently, it might happen that another choice of B-linear combinations provides a simpler form for the reduced interface conditions.

Stokes 2D

We consider the Stokes equations in R 2 . This system is defined by R y = 0, where R ∈ A 3×3 is a matrix with entries in the commutative polynomial ring A = Q(ν, c)[dx, dy] of PD operators in dx = ∂/∂x and dy = ∂/∂y with coefficients in Q(ν, c), where ν is the viscosity and c the reaction coefficient.

We first define A and R. > Delta:=dx^2+dy^2: > L:=-nu*Delta+c: > F3:=AnnOp(Delta*L,R,A);

F 3 :=        1 0 0 dy -dx 0 -ν dx -ν dy 1 0 1 0       
We then obtain different choices for the last row of F . Using the SmithVariablesCompletion procedure, we can try to complement each of the rows of F 3 to a unimodular matrix F whose first two rows are annihilated by 1. This yields distinct choices for F from which we can run the ReducedInterfaceConditions procedure to reduce the interface conditions defined by means of the PD operators 1 and L for the update step, and by dx and dx L for the correction step. > for i from 1 to rowdim(F3) do > F31:=row(F3,i): > SVC:=SmithVariablesCompletion(R,diag(Delta*L),F31,B):

> if SVC=[] then > print([[i],[]]); > else > print([[i], > ReducedInterfaceConditions(R,SVC[3],A,[1,L],[U,V,P]), > ReducedInterfaceConditions(R,SVC[3],A,[dx,dx*L],[U,V,P])]); > fi; > od: [[1], [U = 0, P x = 0], [P = 0, V = 0]] [[2], []] [[3], []]
[ [START_REF] Buchsbaum | What annihilates a module?[END_REF], [P = 0, V = 0], [U = 0, P x = 0]]

Stokes 3D

We consider the Stokes equations in R 3 . This system is defined by R y = 0, where R ∈ A 4×4 is a matrix with entries in the commutative polynomial ring A = Q(ν, c)[dx, dy, dz] of PD operators in dx = ∂/∂x, dy = ∂/∂y and dz = ∂/∂z with coefficients in Q(ν, c), where ν is the viscosity and c the reaction coefficient.

We first define A and R.

R :=        -ν dx 2 + dy 2 + dz 2 + c 0 0 dx 0 -ν dx 2 + dy 2 + dz 2 + c 0 dy 0 0 -ν dx 2 + dy 2 + dz 2 + c dz dx dy dz 0       
Let ∆ = dx 2 + dy 2 + dz 2 and L = -ν ∆ + c. We can check that the Smith normal form S of R is the diagonal matrix which entries are 1, 1, L and L ∆ > S:=map(factor,smith(R,dx));

S :=         1 0 0 0 0 1 0 0 0 0 --ν dx 2 -ν dz 2 -ν dy 2 +c ν 0 0 0 0 - (-ν dx 2 -ν dz 2 -ν dy 2 +c)(dx 2 +dy 2 +dz 2 ) ν        
i.e., R = E S F , where E, F ∈ GL 4 (B) and B = Q(ν, c, dy, dz) [dx]. As a consequence, the residue class of the third (resp., fourth) row of F in the B-module B 1×4 /(B 1×4 R) must be annihilated by the PD operator L (resp., L ∆). Using the AnnOp procedure, we compute the families of generators of the elements of the A-module M = A 1×4 /(A 1×4 R) respectively annihilated by L and by L ∆.

> Delta:=dx^2+dy^2+dz^2: > L:=-nu*Delta+c: > F3:=AnnOp(L,R,A);

F 3 :=                 -dz 0 dx 0 -dy dx 0 0 ν dx dz ν dy dz ν dz 2 -c -dz ν dx dy ν dy 2 + ν dz 2 -c 0 -dy dx dy dz 0 ν dx 2 + ν dy 2 + ν dz 2 -c 0 0 -dx 0 -dz dy 0                
> F4:=AnnOp(L*Delta,R,A);

F 4 :=                 1 0 0 0 dz 0 -dx 0 dy -dx 0 0 -ν dx -ν dy -ν dz 1 0 1 0 0 0 dz -dy 0 0 0 1 0                
We obtain different choices for the last two rows of F . Using the SmithVariablesCompletion procedure, we can try to complement the matrix formed by one row of F 3 and one row of F 4 to a unimodular matrix F whose first two rows are annihilated by 1. This yields distinct choices for F from which we can run the ReducedInterfaceConditions procedure to reduce the interface conditions defined by means of the PD operators 1, 1 and L for the update step, and by dx, dx and dx L for the correction step. We only print out the result when the choice of the last rows of F is a "good" one, i.e., when it can be complemented the matrix formed by the proposed two rows to a unimodular matrix whose first two rows are annihilated by 1.

Inria > B:=DefineOreAlgebra(diff=[dx,x],polynom=[x]
,comm=[dy,dz,nu,c]): > for i from 1 to rowdim(F3) do > for j from 1 to rowdim(F4) do > F34:=stackmatrix(row(F3,i),row(F4,j)): > SVC:=SmithVariablesCompletion(R,diag(L,L*Delta),F34,B):

> if SVC=[] then > else > print([[i,j], > ReducedInterfaceConditions(R,SVC[3],A,[1,1,L],[U,V,W,P]), > ReducedInterfaceConditions(R,SVC[3],A,[dx,dx,dx*L],[U,V,W,P])]); > fi: > od: > od: [[1, 5], [P = 0, V = 0, W x = dz U ], [P = --W ν dy 2 +W c+ν dz dy V dz , P x = 0, W x = - U (-ν dy 2 -ν dz 2 +c) ν dz ]]
[ [START_REF] Achdou | A domain decomposition preconditioner for an advection-diffusion problem[END_REF][START_REF] Chyzak | Non-commutative elimination in Ore algebras proves multivariate identities[END_REF],

[P = 0, W = 0, W x = dz U ], [P = --W ν dy 2 +W c+ν dz dy V dz , P x = 0, W x = 0]] [[2, 5], [P = 0, V = 0, P x = U ν dz 2 -U c -ν W x dz ], [P = -ν dy W dz -V ν dz 2 +V c dy , P x = 0, W x = - U (-ν dy 2 -ν dz 2 +c) ν dz ]] [[2, 7], [P = 0, W = 0, P x = U ν dz 2 -U c -ν W x dz ], [P = -ν dy W dz -V ν dz 2 +V c dy , P x = 0, W x = 0]] [[3, 1], [P = -W c dz , U = 0, P x = 0], [P = 0, V = -W dz dy , W x = -dz Px c ]]
[ [START_REF] Buchberger | An algorithm for finding the basis elements in the residue class ring modulo a zero dimensional polynomial ideal[END_REF][START_REF] Chyzak | Effective algorithms for parametrizing linear control systems over Ore algebras[END_REF], [P = 0, V = 0, W = 0], [U = 0, P x = 0, W x = 0]]

[ [START_REF] Buchsbaum | What annihilates a module?[END_REF][START_REF] Achdou | A domain decomposition preconditioner for an advection-diffusion problem[END_REF], [P = -ν dy W dz -V ν dz 2 +V c dy , U = 0, P x = 0], [P = 0, V = -W dz dy , P

x = -U c + U ν dz 2 -ν W x dz ]]
[ [START_REF] Buchsbaum | What annihilates a module?[END_REF][START_REF] Chyzak | Effective algorithms for parametrizing linear control systems over Ore algebras[END_REF], [P = 0, V = 0, W = 0], [U = 0, P x = 0, W x = 0]]

[ [START_REF] Buchsbaum | What annihilates a module?[END_REF][START_REF] Chyzak | Non-commutative elimination in Ore algebras proves multivariate identities[END_REF], [P = 0, V = 0, W = 0], [U = 0, P x = 0, W x = 0]] [ [START_REF] Chyzak | Non-commutative elimination in Ore algebras proves multivariate identities[END_REF][START_REF] Achdou | A domain decomposition preconditioner for an advection-diffusion problem[END_REF], [U = 0, V = dy W dz , P x = 0], [P = 0, V = -W dz dy , P x = --dz U ν dy 2 -U ν dz 3 +dz U c+Wxν dy 2 +Wxν dz 2 dz ]]

[ [START_REF] Chyzak | Non-commutative elimination in Ore algebras proves multivariate identities[END_REF][START_REF] Chyzak | Effective algorithms for parametrizing linear control systems over Ore algebras[END_REF], [P = 0, V = 0, W = 0], [U = 0, P x = 0, W x = 0]] [ [START_REF] Chyzak | Non-commutative elimination in Ore algebras proves multivariate identities[END_REF][START_REF] Chyzak | Non-commutative elimination in Ore algebras proves multivariate identities[END_REF], [P = 0, V = 0, W = 0], [U = 0, P x = 0,

W x = 0]]
With this method, we can find the Smith variables (i.e., the F 's) that lead to simple reduced interface conditions (for both the update and the correction steps). From the above computations, the simplest reduced interface conditions were obtained, for instance, for a matrix F formed by the seventh row of F 3 and the seventh row of F 4.

Oseen 2D

We consider the Oseen equations in R 2 . This system is defined by R y = 0, where R ∈ A 3×3 is a matrix with entries in the commutative polynomial ring A = Q(ν, c, b1, b2)[dx, dy] of PD operators in dx = ∂/∂x and dy = ∂/∂y with coefficients in Q(ν, c, b1, b2), where ν is the viscosity, c the reaction coefficient, and (b1, b2) the convection velocity.

We first define A and R.

F 3 :=        1 0 0 dy -dx 0 b1 -ν dx -ν dy + b2 1 0 1 0       
We then obtain different choices for the last row of F . Using the SmithVariablesCompletion procedure, we can try to complement each row of F 3 to a unimodular matrix F whose first two rows are annihilated by 1. This yields distinct choices for F which can be used as an input of the ReducedInterfaceConditions procedure to reduce the interface conditions in the different following cases:

Case 1 Let Robin = ν dx -b1/2. The interface conditions are given by the PD operators 1 and L for the update step, and by Robin and dx L for the correction step.

Inria

> Robin:=nu*dx-b1/2: > for i from 1 to rowdim(F3) do > F31:=stackmatrix(row(F3,i)): > SVC:=SmithVariablesCompletion(R,diag(L*Delta),F31,B):

> if SVC=[] then > print([[i],[]]); > else > print([[i], > ReducedInterfaceConditions(R,SVC[3],A,[1,L],[U,V,P]), > ReducedInterfaceConditions(R,SVC[3],A,[Robin,dx*L],[U,V,P])]): > fi: > od: [[1], [U = 0, P x = 0], [P = 0, U = -2 ν dy V b1 ]] [[2], []] [[3], [P = 0, U = -V b2 b1 ], [P = -dy U b1 2 +2 V b1 ν dy 2 -dy V b1 b2 -2 b2 U ν dy 2 +2 b2 U c+2 U b2 2 dy b1 dy , P x = 0]] [[4], [P = 0, V = 0], [V = 2 (-ν dy 2 +c+b2 dy)U b1 dy , P x = 0]]
With this method, we can find the Smith variables (i.e., the F 's) that lead to simple reduced interface conditions (for both the update and the correction steps). From the above computations, the simplest reduced interface conditions were obtained for the matrix formed by the first row of F 3.

Case 2

The interface conditions are given by the PD operators 1 and L for the update step, and by dx and Robin ∆ for the correction step.

> for i from 1 to rowdim(F3) do > F31:=stackmatrix(row(F3,i)): > SVC:=SmithVariablesCompletion(R,diag(L*Delta),F31,B):

> if SVC=[] then > print([[i],[]]); > else > print([[i],
> ReducedInterfaceConditions(R,SVC [START_REF] Buchberger | An algorithm for finding the basis elements in the residue class ring modulo a zero dimensional polynomial ideal[END_REF],A, [1,L],[U,V,P]), > ReducedInterfaceConditions(R,SVC [START_REF] Buchberger | An algorithm for finding the basis elements in the residue class ring modulo a zero dimensional polynomial ideal[END_REF],A,[dx,Robin*Delta],[U,V,P])]): > fi: > od:

[[1], [U = 0, P x = 0], [V = 0, P x = 2 ν dy 2 P -b2 dy U b1 +2 ν dy 2 U b1 -cU b1 b1 ]] [[2], []] [[3], [P = 0, U = -V b2 b1 ], [P = X Y , P x = -V b1 ν dy 2 -b2 U ν dy 2 +b2 U c+U b2 2 dy-dy V b1 b2 -ν dy+b2 ]]
where

X = -3 ν cdy V b1 b2 + 2 ν 2 U b1 2 dy 3 -ν b1 2 dy U c -2 b2 ν 2 cU dy 2 + 2 b2 ν c 2 U + 2 cV b1 ν 2 dy 2 +2 ν 2 V b2 dy 3 b1 -3 ν b1 b2 2 dy 2 V -3 ν b1 2 dy 2 U b2 + b1 b2 3 dy V + 2 ν dy b2 2 U c + b1 2 b2 U c +cV b1 b2 2 + b1 2 b2 2 U dy,
and:

Y = -dy b1 2 ν 2 dy 2 -3 ν dy b2 + b2 2 .

[ [START_REF] Buchsbaum | What annihilates a module?[END_REF], [P = 0, V = 0], [P = -dy V b1 b2 -b1 V c+2 ν cU dy+2 b2 U ν dy 2 -2 ν 2 dy 3 U b1 dy

,

P x = U ν dy 2 -U c -U b2 dy + b1 dy V ]]
With this method, we can find the Smith variables (i.e., the F 's) that lead to simple reduced interface conditions (for both the update and the correction steps). From the above computations, the simplest reduced interface conditions were obtained for the matrix formed by the first row of F 3.

Case 3

The interface conditions are given by the PD operators L and Robin for the update step, and by 1 and dx L for the correction step.

> for i from 1 to rowdim(F3) do > F31:=stackmatrix(row(F3,i)): > SVC:=SmithVariablesCompletion(R,diag(L*Delta),F31,B):

> if SVC=[] then > print([[i],[]]); > else > print([[i], > ReducedInterfaceConditions(R,SVC[3],A,[L,Robin],[U,V,P]), > ReducedInterfaceConditions(R,SVC[3],A,[1,dx*L],[U,V,P])]): > fi: > od: [[1], [U = -2 ν dy V b1 , P x = 0], [P = 0, U = 0]] [[2], []] [[3], [P = 0, P x = -1/2 dy b1 2 U -2 b2 U ν dy 2 +2 b2 U c+2 U b2 2 dy-dy V b1 b2 +2 V b1 ν dy 2 -ν dy+b2 ], [P = -U b1 -V b2 , P x = 0]] [[4], [P = 0, P x = U ν dy 2 -U c -U b2 dy + 1/2 b1 dy V ], [V = 0, P x = 0]]
With this method, we can find the Smith variables (i.e., the F 's) that lead to simple reduced interface conditions (for both the update and the correction steps). From the above computations, the simplest reduced interface conditions were obtained for the matrix formed by the first row of F 3.

Case 4

The interface conditions are given by the PD operators dx and ∆ for the update step, and by 1 and Robin ∆ for the correction step.

Inria

> for i from 1 to rowdim(F3) do > F31:=stackmatrix(row(F3,i)): > SVC:=SmithVariablesCompletion(R,diag(L*Delta),F31,B):

> if SVC=[] then > print([[i],[]]); > else > print([[i], > ReducedInterfaceConditions(R,SVC[3],A,[dx,Delta],[U,V,P]), > ReducedInterfaceConditions(R,SVC[3],A,[1,Robin*Delta],[U,V,P])]): > fi: > od: [[1], [U = -Px b2 dy+c , V = 0], [U = 0, P x = dy (2 ν dy P +2 cν V +2 ν dy V b2 +b1 2 V ) b1 ]] [[2], []] [[3], [P = -b2 dy U b1 +cb2 V +b2 2 dy V +cU b1 b2 dy , P x = -V b1 ν dy 2 -b2 U ν dy 2 +b2 U c+U b2 2 dy-dy V b1 b2 -ν dy+b2 ], [P = -U b1 -V b2 , P x = X Y ]] [[4], [P = -V (b2 dy+c) dy , P x = U ν dy 2 -U c -U b2 dy + b1 dy V ], [P = Z T , V = 0]
] where X, Y, Z and T are four rather involved polynomials which are not printing here.

With this method, we can find the Smith variables (i.e., the F 's) that lead to simple reduced interface conditions (for both the update and the correction steps). From the above computations, the simplest reduced interface conditions were obtained from the first row of F 3.

Oseen 3D

We consider the Oseen equations in R 3 . This system is defined by R y = 0, where R ∈ A 4×4 is a matrix with entries in the commutative polynomial ring A = Q(ν, c, b1, b2, b3)[dx, dy, dz] of PD operators in dx = ∂/∂x, dy = ∂/∂y and dz = ∂/∂z with coefficients in Q(ν, c, b1, b2), where ν is the viscosity, c the reaction coefficient, and (b1, b2, b3) the convection velocity.

We first define A and R.

L := -ν dx 2 + b1 dx -ν dy 2 -ν dz 2 + c + b2 dy + b3 dz R :=        L 0 0 dx 0 L 0 dy 0 0 L dz dx dy dz 0       
Let ∆ = dx 2 + dy 2 + dz 2 and L = -ν ∆ + b1 dx + b2 dy + b3 dz + c. We can check that the Smith normal form S of R is the diagonal matrix which entries are 1, 1, L and L ∆.

> S:=map(factor,smith(R,dx)); > Delta:=dx^2+dy^2+dz^2:

RR n°7953 S :=         1 0 0 0 0 1 0 0 0 0 -L ν 0 0 0 0 - L(dx 2 +dy 2 +dz 2 ) ν         i.e.,
> L:=-nu*Delta+b1*dx+b2*dy+b3*dz+c: > F3:=AnnOp(L,R,A); F 3 :=                 -dz 0 dx 0 -dy dx 0 0 dz ν dx -dz b1 -b2 dz + ν dy dz -c -b3 dz + ν dz 2 -dz dx ν dy -b1 dy ν dy 2 + ν dz 2 -b3 dz -c -b2 dy 0 -dy dx dy dz 0 ν dx 2 -b1 dx + ν dy 2 + ν dz 2 -b3 dz -c -b2 dy 0 0 -dx 0 -dz dy 0                 > F4:=AnnOp(L*Delta,R,A); F 4 :=                 1 0 0 0 dz 0 -dx 0 dy -dx 0 0 b1 -ν dx -ν dy + b2 b3 -dz ν 1 0 1 0 0 0 dz -dy 0 0 0 1 0                
We then obtain different choices for the last rows of F . Using the SmithVariablesCompletion procedure, we can try to complement a matrix formed by one row of F 3 and one row of F 4 to a unimodular matrix F whose first two rows are annihilated by 1. This yields distinct choices for F which can be used as an input of the ReducedInterfaceConditions procedure to reduce the interface conditions in the following distinct cases: Case 1 Let Robin = ν dx -b1/2. The interface conditions are defined by means of the PD operators 1, L and 1 for the update step, and by Robin, Robin and dx L for the correction step. We run the algorithms for all the distinct choices for F . > Robin:=nu*dx-b1/2: > for i from 1 to rowdim(F3) do > for j from 1 to rowdim(F4) do > F34:=stackmatrix(row(F3,i),row(F4,j)); > SVC:=SmithVariablesCompletion(R,diag(L,L*Delta),F34,B);

> if SVC=[] then > s:=[[i,j],[]]: > else > s:=[[i,j],ReducedInterfaceConditions(R,SVC[3],A,[1,L,1],[U,V,W,P]), > ReducedInterfaceConditions(R,SVC[3],A,[Robin,Robin,dx*L],[U,V,W,P])]: > fi: > od: > od:
Since the results are quite large, we only give the simplest reduced interface conditions (for both the update and the correction steps of the algorithm). It is obtained by choosing the first row of F 3 and the seventh row of F 4:

[[P = 0, W = 0, W x = dz U ], [P = -1/2 Wxb1 2 +b1 2 dz U +2 b1 dz ν dy V -4 ν 2 dy 2 Wx+4 b3 dz ν Wx+4 cν Wx+4 b2 dy ν Wx dz b1 , W = 2 ν Wx b1 , P x = 0]]
Case 2 The interface conditions are defined by means of the PD operators 1, L and 1 for the update step, and by Robin, Robin ∆ and dx for the correction step. We run the algorithms for all the distinct choices for F .

> for i from 1 to rowdim(F3) do > for j from 1 to rowdim(F4) do > F34:=stackmatrix(row(F3,i),row(F4,j)); > SVC:=SmithVariablesCompletion(R,diag(L,L*Delta),F34,B);

> if SVC=[] then > s:=[[i,j],[]]: > else > s:=[[i,j],ReducedInterfaceConditions(R,SVC[3],A,[1,L,1],[U,V,W,P]), > ReducedInterfaceConditions(R,SVC[3],A,[Robin,Robin*Delta,dx],[U,V,W,P])]: > fi: > od: > od:
Since the results are quite large, we only give the simplest reduced interface conditions (for both the update and the correction steps of the algorithm). It is obtained by choosing the first row of F 3 and the seventh row of F 4:

[[P = 0, W = 0, W x = dz U ], [P = -2 ν dz Px-2 W ν dy 2 b1 -2 ν W b1 dz 2 +b1 W b3 dz +b1 W c+b1 W b2 dy dz b1 , U = -2 ν (-2 dz Px+W dy 2 b1 +2 W b1 dz 2 +dy b1 dz V ) dz b1 2 , W x = 0]] Case 3
The interface conditions are defined by means of the PD operators 1, L and Robin for the update step, and by Robin, dx L and 1 for the correction step. We run the algorithms for all the distinct choices for F .

> for i from 1 to rowdim(F3) do > for j from 1 to rowdim(F4) do > F34:=stackmatrix(row(F3,i),row(F4,j)); > SVC:=SmithVariablesCompletion(R,diag(L,L*Delta),F34,B);

> if SVC=[] then > s:=[[i,j],[]]: > else > s:=[[i,j],ReducedInterfaceConditions(R,SVC[3],A,[1,L,Robin],[U,V,W,P]), > ReducedInterfaceConditions(R,SVC[3],A,[Robin,dx*L,1],[U,V,W,P])]: > fi: > od: > od:
Since the results are quite large, we only give the simplest reduced interface conditions (for both the update and the correction steps of the algorithm). It is obtained by choosing the first row of F 3 and the seventh row of F 4:

[[P = 0, W = 2 ν dz U b1 , W x = dz U ], [W = 0, P x = 0, W x = -dz (U b1 +2 P +2 ν dy V ) b1 ]]
Case 4 The interface conditions are defined by means of the PD operators 1, ∆ and dx for the update step, and by Robin, Robin ∆ and 1 for the correction step. We run the algorithms for all the distinct choices for F .

> for i from 1 to rowdim(F3) do > for j from 1 to rowdim(F4) do > F34:=stackmatrix(row(F3,i),row(F4,j)); > SVC:=SmithVariablesCompletion(R,diag(L,L*Delta),F34,B);

> if SVC=[] then > s:=[[i,j],[]]: > else > s:=[[i,j], > ReducedInterfaceConditions(R,SVC[3],A,[1,Delta,dx],[U,V,W,P]), > ReducedInterfaceConditions(R,SVC[3],A,[Robin,Robin*Delta,1],[U,V,W,P])]: > fi: > od: > od:
Since the results are quite large, we only give the simplest reduced interface conditions (for both the update and the correction steps of the algorithm). It is obtained by choosing the first row of F 3 and the seventh row of F 4:

[P = -(b3 dz +c+b2 dy)W dz , U = 0, W x = 0] [W = 0, P x = X 2 b1 ν , W x = -dz (U b1 +2 P +2 ν dy V ) b1 ]]
where:

X =b1 2 P + b1 3 U + 2 ν b3 dz U b1 + 2 b1 ν cU + 2 ν b1 2 dy V + 2 ν dy b1 b2 U + 4 cV ν 2 dy + 4 ν 2 dy 2 V b2 + 4 b3 dz V ν 2 dy + 4 ν b2 P dy + 4 ν cP + 4 ν b3 dz P.
The computations were done by the Maple Schwarz package built upon OreModules ( [START_REF] Chyzak | OreModules: A symbolic package for the study of multidimensional linear systems[END_REF]).

Inria

The unimodular matrices U and V over the ring C = Q(λ, µ, dy)[dx] are such that U R V = S. Equivalently defining E = U -1 and F = V -1 , we have E S F = R. In our case, the matrices E and F are the following ones:

> E:=inverse(U); > F:=inverse(V);

E :=   (λ + µ) dx dy µ dy 2 2 µ dy 2 + dy 2 λ + µ dx 2 µ 2 dx dy 3 (λ+µ)   F :=   - (-dy 2 λ+µ dx 2 )dx dy 3 (λ+µ) 1 1 0  
Given two operators Op1 and Op2, for instance, Op1 = 1 and Op2 = dx 2 + dy 2 , the interface conditions are computed from the matrix F (here only from the second row of F ) as follows:

> Op1:=1: > IC_1:=Mult(Op1,linalg[submatrix](F,2..2,1..2),A); IC 1 := 1 0 > Op2:=dx^2+dy^2: > IC_2:=Mult(Op2,linalg[submatrix](F,2..2,1..2),A); IC 2 := dy 2 + dx 2 0
We then reduce the interface conditions with respect to the Gröbner basis of the equations of the system. We take the first interface condition. We then do the same with the second interface condition:

> IC2:=evalm(IC_2&*evalm([[u],[v]])); IC2 := dy 2 + dx 2 u > NIC2:='OreModules/normal_form'(IC2[1,1],GB,mTord); N IC2 := -(λ+µ)dx dy v 2 µ+λ + dy 2 u(λ+µ) 2 µ+λ
Finally, we perform linear algebra simplifications to the system formed by the two normal forms N IC1 and N IC2. To avoid any multiplication by dx in the computation, we first replace dx u and dx v by the jet variables u x and v x in N IC1 and N IC2 and subtract them by right hand sides f 1 and f 2 .

> M1:=subs(u=u[x],v=v[x],coeff(NIC1,dx,1))+coeff(NIC1,dx,0)-f[1]; M 1 := u -f 1 > M2:=subs(u=u[x],v=v[x],coeff(numer(NIC2),dx,1))+coeff(numer(NIC2),dx,0)-f[2]; Inria M 2 := -(λ + µ) dy v x + dy 2 u (λ + µ) -f 2
We finally solve {M 1 = 0, M 2 = 0} in the unknowns u x , v x , u, v, and put f 1 = f 2 = 0 in the result to obtain the reduced interface conditions.

> Sols:=subs(f1=0,f2=0,solve(M1,M2,u

[x],v[x],u,v)); Sols := {u = 0, u x = u x , v = v, v x = 0}
We obtain the following reduced interface conditions:

u(x, y) = 0, ∂ v(x, y) ∂ x = 0.
The result can be directly obtained using the procedure ReducedInterfaceConditions. The procedure takes as inputs the matrix R of the system, a unimodular matrix F corresponding to the Smith normal form of R, the commutative polynomial ring A, the PD operators Op1 = 1 and Op2 = dx 2 + dy 2 defining the interface conditions and the names of the variables u and v. The ouput contains the reduced interface conditions.

> ReducedInterfaceConditions(R,F,A,[1,dx^2+dy^2],[u,v]); [u = 0, v x = 0]
Given the matrix R, the unimodular matrices E and F satisfying R = E S F , where S is the Smith normal form of R, are not unique. Moreover, as we have noticed before, the interface conditions depend on the choice of F . Thus, we can obtained distinct reduced interface conditions for distinct choices of F . For instance, we can consider another F :

F 1 :=   1 - dx (3 µ dy 2 +2 dy 2 λ+2 µ dx 2 +dx 2 λ) dy 3 (λ+µ) 0 1  
We then obtain the following reduced interface conditions:

> ReducedInterfaceConditions(R,F1,A, [1,dx^2+dy^2], [u,v]); [v = 0, u x = 0]

Elasticity 3D

We consider the elastostatic equations (i.e., the Navier-Cauchy equations) in R 3 . This system is defined by R y = 0, where R ∈ A 3×3 is a matrix with entries in the commutative polynomial ring A = Q(λ, µ)[dx, dy, dz] of PD operators in dx, dy, dz with coefficients in Q(λ, µ), where λ and µ are the two Lamé constants, dx = d/dx, dy = d/dy and dz = d/dz are the derivations with respect to x respectively y and z.

We first define A and R. The Smith normal form of R with respect to the variable dx is given by: The reduced interface conditions depend on the unimodular matrices E and F satisfying R = E S F . Let us first consider the F returned by the Maple procedure smith.

> F:=inverse(V);

F :=       - dx (
2 dx 2 µ 2 +µ dx 2 λ+dz 2 µ 2 -dy 2 λ 2 -2 µ dy 2 λ) dy (λ+µ)(2 dy 2 µ+dy 2 λ+dz 2 µ)

1 -dz (dx 2 µ-dy 2 λ-dy 2 µ) dy (2 dy 2 µ+dy 2 λ+dz 2 µ) dx dz dy 2 +dz 2 0 1

--dy 2 -dz 2 +dx dz dy 2 +dz 2 0 -1

     
Choosing this particular F , the interface conditions are the following: We reduce them using the ReducedInterfaceConditions procedure. This yields:

> ReducedInterfaceConditions(R,F,A,[dx,dx,dx*Delta],[u,v,w]);

[v = -dz w dy , u x = 0, v x = --wxdy 2 λ+dz udy 2 µ-2 wxdy 2 µ+udz 3 µ-wxdz 2 µ dy dz (λ+µ) ] However, choosing others F provides distinct (simpler) reduced interface conditions. With [v = -dz w dy , u x = 0, v x = dy wx dz ]

Stokes 2D

We consider the Stokes equations in R 2 . This system is defined by R y = 0, where R ∈ A 3×3 is a matrix with entries in the commutative polynomial ring A = Q(ν, c)[dx, dy] of PD operators in dx, dy with coefficients in Q(ν, c), where ν is the viscosity and c the reaction coefficient, dx = d/dx and dy = d/dy are the derivations with respect to x respectively y.

We first define A and R. The Smith normal form of R with respect to the variable dx is given by: We first define A and R. The Smith normal form of R with respect to the variable dx is given by: > S:=smith(R,dx,'U','V');

S :=         1 0 0 0 0 1 0 0 0 0 dx 2 --ν dy 2 -ν dz 2 +c ν 0 0 0 0 - (-ν dx 2 -ν dy 2 -ν dz 2 +c)(dx 2 +dy 2 +dz 2 ) ν        
The reduced interface conditions depend on the unimodular matrices E and F satisfying R = E S F . Let us consider, for instance, the F returned by the Maple procedure smith. [p = 0, v = -dz w dy , p x = --dz uν dy 2 -uν dz 3 +dz uc+wxν dy 2 +wxν dz 2 dz ] 8.2.5 Oseen 2D

We consider the Oseen equations in R 2 . This system is defined by R y = 0, where R ∈ A 3×3 is a matrix with entries in the commutative polynomial ring A = Q(ν, c, b1, b2)[dx, dy] of PD operators in dx, dy with coefficients in Q(ν, c, b1, b2), where ν is the viscosity, c the reaction coefficient, (b1, b2) the convection velocity, dx = d/dx and dy = d/dy are the derivations with respect to x respectively y.

We first define A and R. The Smith normal form of R with respect to the variable dx is given by: The reduced interface conditions depend on the unimodular matrices E and F satisfying R = E S F . Let us consider, for instance, the F returned by the Maple procedure smith. Choosing this particular F , we can compute the reduced interface conditions using the Re-ducedInterfaceConditions procedure.

Four cases of interface conditions have to be distinguished. To write these four cases, we introduce the following PD operators: The reduced interface conditions depend on the unimodular matrices E and F satisfying R = E S F . Let us consider, for instance, the F returned by the Maple procedure smith. Choosing this particular F , we can compute the reduced interface conditions using the Re-ducedInterfaceConditions procedure.

Four cases of interface conditions have to be distinguished. Here, we shall only give the computed reduced interface conditions for the first case: the other ones can be obtained in a similar way but we then get huge expressions that are not readable. To write these four cases, we define the following PD operators: [u = 0, w = dz v dy , p x = 0] The computations were done by the Maple Schwarz package built upon OreModules ( [START_REF] Chyzak | OreModules: A symbolic package for the study of multidimensional linear systems[END_REF]).

Theorem 1 .

 1 Let k be a field, R = k[s], p a positive integer, and A ∈ R p×p . Then, there exist two matrices E ∈ GL p (R) and F ∈ GL p (R) such that A = E D F, where D = diag(d 1 , . . . , d p ) and the d i 's are elements of R satisfying d 1 | d 2 | • • • | d p .

  and k is a field.

Example 4 .

 4 Let us consider the underdetermined linear PD equation ∆ fg = 0 in the two unknowns f and g

  linear combination of the representatives of the generators of ann M (d j ) for j = r + 1, . . . , p. The computation of the ann M (d j )'s provides Inria Smith variables with a physical meaning contrary to the direct computation of ann Ri⊗ R M (d j ) by means of Smith normal forms.

Example 5 .

 5 a fact which can be constructively checked by means of Gröbner basis techniques (see Algorithm 10 of Section 7). If such a factorization exists, then we get a matrix Z ∈ R r×(p-r) i, and thus the matrix F 1 = T + Z F 2 = Λ A satisfies (46). Let us consider the elastostatic equations or Navier-Cauchy equations in R 3 , i.e.,

F 2 )Example 6 .

 26 of rank r (e.g., by computing a Smith normal form of F 2 as explained at the end of Section 5.1). (b) Compute a left inverse T ∈ R r×p i of Q (see the end of Section 5.1).(c) Factorize T with respect to(F T 2 A T ) T over R i . RR n°7953 (d) If such a factorization does not exist, then return [ ], else if T = -Z F 2 + Λ A, where Z ∈ R r×(p-r) i and Λ ∈ R r×p i , then return the matrix F = T + Z F 2 F 2 .Algorithm 6 is implemented in the Maple package Schwarz built upon OreModules ([START_REF] Chyzak | OreModules: A symbolic package for the study of multidimensional linear systems[END_REF]) and is demonstrated in Appendix 8. Consider again the elostatic equations introduced in Example 5. For the choice of F 2 given at the end of Example 5, we can find a completion of F 2 defined by:

>

  A:=DefineOreAlgebra(diff=[dx,x],diff=[dy,y],diff=[dz,z],polynom=[x,y,z], > comm=[lambda,mu]): > R := matrix(3, 3,[-2*dx^2*mu-dx^2*lambda-dy^2*mu-dz^2*mu,-dx*dy*(lambda+mu), > -dx*dz*(lambda+mu),-dx*dy*(lambda+mu),-dx^2*mu-2*dy^2*mu-dy^2*lambda-dz^2*mu, > -dy*dz*(lambda+mu),-dx*dz*(lambda+mu),-dy*dz*(lambda+mu), > -dx^2*mu-dy^2*mu-2*dz^2*mu-dz^2*lambda]);

>

  A:=DefineOreAlgebra(diff=[dx,x],diff=[dy,y],polynom=[x,y],comm=[nu,c]): > R:=evalm([[-nu*(dx^2+dy^2)+c,0,dx],[0,-nu*(dx^2+dy^2)+c,dy],[dx,dy,0]]); 2 + dy 2 + c 0 dx 0 -ν dx 2 + dy 2 + c dx 2 + dy 2 and L = -ν ∆ + c.We can check that the Smith normal form S of R is the diagonal matrix which entries are 1, 1 and L ∆ > S:=map(factor,smith(R,dx)); (dx 2 +dy 2 )(-ν dx 2 -ν dy 2 +c) ., R = E S F , where E, F ∈ GL 3 (B) and B = Q(ν, c, dy, dz)[dx]. As a consequence, the residue class of the third row of F in the B-module B 1×3 /(B 1×3 R) must be annihilated by the PD operator L ∆. We compute a family of generators of the elements of the A-module M = A 1×3 /(A 1×3 R) annihilated by L ∆ using the AnnOp procedure.

Inria

>

  B:=DefineOreAlgebra(diff=[dx,x],polynom=[x],comm=[dy,nu,c]):

>

  A:=DefineOreAlgebra(diff=[dx,x],diff=[dy,y],diff=[dz,z],polynom=[x,y,z], > comm=[nu,c]): > R:=evalm([[-nu*(dx^2+dy^2+dz^2)+c,0,0,dx],[0,-nu*(dx^2+dy^2+dz^2)+c,0,dy], > [0,0,-nu*(dx^2+dy^2+dz^2)+c,dz],[dx,dy,dz,0]]);

>

  A:=DefineOreAlgebra(diff=[dx,x],diff=[dy,y],polynom=[x,y],comm=[nu,c,b1,b2]): > R:=evalm([[-nu*(dx^2+dy^2)+c+b1*dx+b2*dy,0,dx], > [0,-nu*(dx^2+dy^2)+c+b1*dx+b2*dy,dy],[dx,dy,0]]); 2 + dy 2 + c + b1 dx + b2 dy 0 dx 0 -ν dx 2 + dy 2 + c + b1 dx + b2 dy dx 2 + dy 2 and L = -ν ∆ + b1 dx + b2 dy + c. We can check that the Smith normal form S of R is the diagonal matrix which entries are 1, 1 and L ∆ > S:=map(factor,smith(R,dx)); (dx 2 +dy 2 )(-ν dx 2 -ν dy 2 +c+b1 dx +b2 dy) ., R = E S F , where E, F ∈ GL 3 (B) and B = Q(ν, c, b1, b2, dy)[dx]. As a consequence, the residue class of the third row of F in the B-module B 1×3 /(B 1×3 R) must be annihilated by the PD operator L ∆. Using the AnnOp procedure, we compute the family of generators of the elements of the A-module M = A 1×3 /(A 1×3 R) annihilated by L ∆. > Delta:=dx^2+dy^2: > L:=-nu*Delta+b1*dx+b2*dy+c: > F3:=AnnOp(L*Delta,R,A);

>

  A:=DefineOreAlgebra(diff=[dx,x],diff=[dy,y],diff=[dz,z],polynom=[x,y,z], > comm=[nu,c,b1,b2,b3]): > L:=-nu*dx^2-nu*dy^2-nu*dz^2+b1*dx+b2*dy+b3*dz+c; > R:=evalm([[L,0,0,dx],[0,L,0,dy],[0,0,L,dz],[dx,dy,dz,0]]);

  R = E S F , where E, F ∈ GL 4 (B) and B = Q(ν, c, b1, b2, b3, dy, dz)[dx]. As a consequence, the residue class of the third (resp., fourth) row of F in the B-module B 1×4 /(B 1×4 R) must be annihilated by the PD operator L (resp., L ∆). Using the AnnOp procedure, we compute the families of generators of the elements of the A-module M = A 1×4 /(A 1×4 R) respectively annihilated by L and L ∆.

Inria

>

  B:=DefineOreAlgebra(diff=[dx,x],polynom=[x],comm=[dy,dz,nu,c,b1,b2,b3]):

>

  IC1:=evalm(IC_1&*evalm([[u],[v]])); IC1 := u Then we compute its normal form with respect to GB. > NIC1:='OreModules/normal_form'(IC1[1,1],GB,mTord); N IC1 := u

>- 2

 2 A:=DefineOreAlgebra(diff=[dx,x],diff=[dy,y],diff=[dz,z],polynom=[x,y,z], > comm=[lambda,mu]): > R := matrix(3, 3,[-2*dx^2*mu-dx^2*lambda-dy^2*mu-dz^2*mu,-dx*dy*(lambda+mu), > -dx*dz*(lambda+mu),-dx*dy*(lambda+mu),-dx^2*mu-2*dy^2*mu-dy^2*lambda-dz^2*mu, > -dy*dz*(lambda+mu),-dx*dz*(lambda+mu),-dy*dz*(lambda+mu), > -dx^2*mu-dy^2*mu-2*dz^2*mu-dz^2*lambda]); dx 2 µdx 2 λdy 2 µdz 2 µ -dx dy (λ + µ) -dx dz (λ + µ)-dx dy (λ + µ) -dx 2 µ -2 dy 2 µdy 2 λdz 2 µ -dy dz (λ + µ)-dx dz (λ + µ) -dy dz (λ + µ) -dx 2 µdy 2 µ -2 dz 2 µdz 2 λ

> 4

 4 S:=smith(R,dx,'U','V'); + 2 dx 2 dy 2 + 2 dy 2 dz 2 + 2 dz 2 dx 2 + dx 4 + dz4

>

  Delta:=dx^2+dy^2+dz^2: > IC1:=Mult(dx,linalg[submatrix](F,2..2,1..3),A); > IC2:=Mult(dx,linalg[submatrix](F,3..3,1..3),A); > IC3:=Mult(dx*Delta,linalg[submatrix](F,3..3,1..3),A);IC1 := dz dx 2 dy 2 +dz 2 0 dx IC2 := -dx (-dy 2 -dz 2 +dx dz ) dy 2 +dz 2 0 -dx IC3 := -dx (dy 2 +dz 2 +dx 2 )(-dy 2 -dz 2 +dx dz )dy 2 +dz 2 0 -dx dz 2dx dy 2dx 3

F 1 ( 3 > 1 ->

 131 dy 2 µ+2 dy 2 λ+2 dx 2 µ+2 dz 2 µ+dx 2 λ+dz 2 λ) dy (dy 2 +dz 2 )(λ+µ) ReducedInterfaceConditions(R,F1,A,[dx,dx,dx*Delta],[u,v,w]); [u = 0, v x = 0, w x = 0dx (3 dy 2 µ+2 dy 2 λ+2 dx 2 µ+2 dz 2 µ+dx 2 λ+dz 2 λ) dy (dy 2 +dz 2 )(λ+µ) ReducedInterfaceConditions(R,F2,A,[dx,dx,dx*Delta],[u,v,w]);

>

  A:=DefineOreAlgebra(diff=[dx,x],diff=[dy,y],polynom=[x,y],comm=[nu,c]): > R:=evalm([[-nu*(dx^2+dy^2)+c,0,dx],[0,-nu*(dx^2+dy^2)+c,dy],[dx,dy,0]]);

>

  S:=smith(R,dx,'U','V'); -2 dy 2 ν dx 2 -dy 4 ν+dy 2 c-ν dx 4 +dx 2 c ν The reduced interface conditions depend on the unimodular matrices E and F satisfying R = E S F . Let us consider, for instance, the F returned by the Maple procedure smith.

  Choosing this particular F , the interface conditions are the following:RR n°7953> L:=-nu*(dx^2+dy^2)+c: > IC1:=Mult(dx,linalg[submatrix](F,3..3,1..3),A); > IC2:=Mult(dx*L,linalg[submatrix](F,3..3,1..3),A);IC1 := dx 0 0 IC2 := -ν dx 3dx ν dy 2 + dx c 0 0We reduce them using the ReducedInterfaceConditions procedure. This yields:> ReducedInterfaceConditions(R,F,A,[dx,dx*L],[u,v,p]); [p = 0, v = 0]8.2.4 Stokes 3DWe consider the Stokes equations in R 3 . This system is defined by R y = 0, where R ∈ A 4×4 is a matrix with entries in the commutative polynomial ring A = Q(ν, c)[dx, dy, dz] of PD operators in dx, dy, dz with coefficients in Q(ν, c), where ν is the viscosity and c the reaction coefficient, dx = d/dx, dy = d/dy and dz = d/dz are the derivations with respect to x respectively y and z.

>

  A:=DefineOreAlgebra(diff=[dx,x],diff=[dy,y],diff=[dz,z],polynom=[x,y,z], > comm=[nu,c]): > R:=evalm([[-nu*(dx^2+dy^2+dz^2)+c,0,0,dx],[0,-nu*(dx^2+dy^2+dz^2)+c,0,dy], > [0,0,-nu*(dx^2+dy^2+dz^2)+c,dz],[dx,dy,dz,0]]); 2 + dy 2 + dz 2 + c 0 0 dx 0 -ν dx 2 + dy 2 + dz 2 + c 0 dy 0 0 -ν dx 2 + dy 2 + dz 2 + c

-Inria

  2 -ν dy 2 -ν dz 2 +c dy -dy 2 -dz 2 +dz dx dy 2 +dz 2 Choosing this particular F , the interface conditions are the following: > L:=-nu*(dx^2+dy^2+dz^2)+c: > IC1:=Mult(dx,linalg[submatrix](F,3..3,1..4),A); > IC2:=Mult(dx,linalg[submatrix](F,4..4,1..4),A); > IC3:=Mult(dx*L,linalg[submatrix](F,4..4,1..4),A);IC1 := dz dx 2 dy 2 +dz 2 0 dx 0 IC2 := -dx (-dy 2 -dz 2 +dx dz ) dy 2 +dz 2 0 -dx 0 IC3 := -(-dy 2 -dz 2 +dx dz )dx (-ν dy 2 -ν dz 2 -ν dx 2 +c) dy 2 +dz 2 0 dx ν dy 2 + dx ν dz 2 + ν dx 3dx c 0We reduce them using the ReducedInterfaceConditions procedure. This yields:> ReducedInterfaceConditions(R,F,A,[dx,dx,L*dx],[u,v,w,p]);

>

  A:=DefineOreAlgebra(diff=[dx,x],diff=[dy,y],polynom=[x,y],comm=[nu,c,b1,b2]): > R:=evalm([[-nu*(dx^2+dy^2)+c+b1*dx+b2*dy,0,dx], > [0,-nu*(dx^2+dy^2)+c+b1*dx+b2*dy,dy],[dx,dy,0]]); 2 + dy 2 + c + b1 dx + b2 dy 0 dx 0 -ν dx 2 + dy 2 + c + b1 dx + b2 dy dy

>

  S:=map(factor,smith(R,dx,'U','V')); (dx 2 +dy 2 )(-ν dx 2 +b1 dx -ν dy 2 +c+b2 dy)

>

  F:=inverse(V); dx 2 +b1 dx -ν dy 2 +c+b2 dy dy

>Inria

  Delta:=dx^2+dy^2: > L2:=-nu*(dx^2+dy^2)+b1*dx+b2*dy+c: > Robin:=nu*dx-b1/2: Case 1 Correction step The interface conditions are the following:> IC1:=Mult(Robin,linalg[submatrix](F,3..3,1..3),A); > IC2:=Mult(L2*dx,linalg[submatrix](F,3..3,1..3),A); IC1 := ν dx -1/2 b1 0 0 IC2 := -ν dx 3 + b1 dx 2dx ν dy 2 + dx c + dx b2 dy 0 0We reduce them:> ReducedInterfaceConditions(R,F,A,[Robin,L2*dx],[u,v,p]); [p = 0, u = -2 ν dy v b1 ]Update step The interface conditions are the following: > IC1:=Mult(L2,linalg[submatrix](F,3..3,1..3),A); > IC2:=Mult(1,linalg[submatrix](F,3..3,1..3),A); IC1 := -ν dx 2 + b1 dxν dy 2 + c + b2 dy 0 0 IC2 := 1 0 0We reduce them:> ReducedInterfaceConditions(R,F,A,[L2,1],[u,v,p]); [u = 0, p x = 0]Case 2Correction step The interface conditions are the following: > IC1:=Mult(Robin*Delta,linalg[submatrix](F,3..3,1..3),A); > IC2:=Mult(dx,linalg[submatrix](F,3..3,1..3),A); IC1 := ν dx 3 + dx ν dy 2 -1/2 b1 dx 2 -1/2 dy 2 b1 0 0 IC2 := dx 0 0

  2 +b1 dx -ν dy 2 -ν dz 2 +c+b2 dy+b3 dz dy

>-2 ν 2 wdz dy 2 b1+-2 ν 2 dy vdz 2 b1+ 1 / 2 ( 4 -2 ν 2 wdz 3 b1+ 2 b3 ν wdz 2 b1+ 1 ]0 - 1 0

 221243211 Delta:=dx^2+dy^2+dz^2: > Robin:=nu*dx-b1/2: Case 1 Correction step The interface conditions are the following: > IC1:=Mult(Robin,linalg[submatrix](F,3..3,1..4),A); > IC2:=Mult(Robin,linalg[submatrix](F,4..4,1..4),A); > IC3:=Mult(dx*L,linalg[submatrix](F,4..4,1..4),A);IC1 := -1/2 (-2 ν dx +b1 )dz dx dy 2 +dz 2 0 ν dx -1/2 b1 0 IC2 := 1/2 (-2 ν dx +b1 )(-dy 2 -dz 2 +dx dz ) dy 2 +dz 2 0 -ν dx + 1/2 b1 0 IC3 := -(-dy 2 -dz 2 +dx dz )dx (-ν dx 2 -ν dy 2 -ν dz 2 +b1 dx +b2 dy+b3 dz +c) dy 2 +dz 2 0 -dx L 0We reduce them:> ReducedInterfaceConditions(R,F,A,[Robin,Robin,dx*L],[u,v,w,p]); RR n°7953 [ p = 0, u = -2 ν dy v b1 -2 ν dz w b1 , p x = -2 ν 2 dy 3 v b1 b2 ν w+4 ν b3 v)dz dy b1 + 1/2 (vb1 2 +4 ν cv)dy b1 /2 (2 b1 2 w-2 ν wxb1 +4 cν w)dz b1Update step The interface conditions are the following:> IC1:=Mult(1,linalg[submatrix](F,3..3,1..4),A); > IC2:=Mult(L,linalg[submatrix](F,4..4,1..4),A); > IC3:=Mult(1,linalg[submatrix](F,4..4,1..4),A); -(-ν dx 2 -ν dy 2 -ν dz 2 +b1 dx +b2 dy+b3 dz +c)(-dy 2 -dz 2 +dx dz ) dy 2 +dz 2 0 -L 0 IC3 := --dy 2 -dz 2 +dx dz dy 2 +dz 2We reduce them:> ReducedInterfaceConditions(R,F,A,[1,L,1],[u,v,w,p]);

RR n°7953

i=1 A ij m j = 0 for i = 1, . . . , 4. We can check thatFitt0 (M ) = R det(A) = R (∆ L 2 3 ), ann R (M ) = R (∆ L 3 ), RR n°7953

Work supported by the PEPS Maths-ST2I SADDLEShttp://www-math.unice.fr/ ~dolean/saddles/

Algorithm 8.

• Input: A finite set {P j } j=1,...,r of polynomials of R = k[x 1 , . . . , x n ] and an admissible term order ≺ on Mon(R).

• Output: A Gröbner basis G = {Q i } i=1,...,t of the ideal I = (P 1 , . . . , P r ) of R.

1. Set G = {P j } j=1,...,r and denote by P the set of distinct pairs of elements of G.

2. While P = ∅, do:

• Choose a pair (P i , P j ) of P and remove it from P.

• Compute S(P i , P j ) and reduce it modulo G.

• If the remainder S is nonzero, then:

-Add to P all the pairs of the form (P, S), where P ∈ G.

-Add S to G.

Return G.

Let us illustrate Algorithm 8 with an explicit example.

Example 11. We consider again Example 10. Let us compute a Gröbner basis of the ideal I. We first set G = {P, Q}. In Example 10, we computed S(P, Q) = -1 2 x 1 x 2 x 3 + x 3 and reduced it by P to obtain the following polynomial:

Then, we set G = {P, Q, T }. Now, we need to compute the following S-polynomials:

Reducing these new S-polynomials modulo the elements of G, we obtain:

Hence, Buchberger's algorithm terminates with the Gröbner basis G = {P, Q, T } of I.

Definition 7.

A Gröbner basis G = {Q i } i=1,...,t of an ideal I of R is called minimal if the following two conditions are satisfied:

1. For all i = 1, . . . , t, lc(Q i ) = 1.

For all

Proposition 2. Every nonzero ideal I of R = k[x 1 , . . . , x n ] admits a minimal Gröbner basis for any given admissible term order ≺.

Remark 6. For computational purposes, the degree reverse lexicographical order is much more feasible than the purely lexicographical one. The latter serves elimination purposes: Buchberger's algorithm, applied to a generating set of an ideal with respect to the lexicographical order, finds a polynomial of the ideal I (if it exists) which contains only the ≺-least variable.

Inria (e) If any H i contains any λ i , then return [ ], else let H i = q ′ j=1 A ′′ ij µ j , for i = 1, . . . , q, and return

3. Computation of Bézout identities: Given a matrix A ∈ R q×p , find (if it exists) a left inverse

The left inverse can be computed as follows:

Algorithm 11.

• Input: A commutative polynomial ring R and a matrix A ∈ R q×p .

• Output:

(a) Introduce the indeterminates λ j , j = 1, . . . , p and µ i , i = 1, . . . , q, over R and define:

R µ i with respect to a term order which eliminates the λ j 's.

(c) Remove from G the elements which do not contain any λ i and call H this new set.

The computation of the right inverse of a matrix over a commutative polynomial ring can be done by applying the previous algorithm to the transpose of the matrix and transposing the obtained matrix.

Let us illustrate the use of Gröbner bases to solve some of the previous problems.

Example 14. Let us consider again Examples 10 and 11. In Example 10, we found:

Substituting T into the formulas for S(P, T ) and S(Q, T ) obtained in Example 11, we get:

Then, the identities (49) become the following two relations

)) because the second equation of the previous system is a multiple of the first one.

Finally, using the Gröbner basis G = {P, Q, T } of the ideal I = (P, Q) of R generated by P and Q, we can prove that the ring R/I is 1-dimensional, which means that the dimension of the complex solutions V (I) of the polynomial system {P = 0, Q = 0} is 1, which can be checked by direct computation:

Reduction of interface conditions

In this section, we illustrate the problem studied in Section 6.

Elasticity 2D

> restart:

> with(linalg):

> with(OreModules):

> with(Schwarz):

We consider the system of the linear elasticity equations in R 2 . This system is defined by R y = 0, where R ∈ A 2×2 is a matrix with entries in the commutative polynomial ring A = Q(λ, µ)[dx, dy] of PD operators in dx, dy with coefficients in Q(λ, µ), where λ and µ are the two Lamé constants, dx = d/dx and dy = d/dy are the derivations with respect to x respectively y.

We define A and R.

The equations can be written: We then define an appropriate term order to sort the indeterminates of the commutative polynomial ring B and compute a Gröbner basis of the set of equations G with respect to this term order.

> mTord:='OreModules/term_order'(B [START_REF] Achdou | A domain decomposition preconditioner for an advection-diffusion problem[END_REF],tdeg(dx,u,v,x),[u,v]);

> GB:='OreModules/gb'(G,mTord); GB := [dx dy uλ + dx dy uµ + 2 vµ dy 2 + vdy 2 λ + vµ dx 2 , 2 uµ dx 2 + udx 2 λ + uµ dy 2 + dx dy vλ + dx dy vµ]

We now show how to reduce the interface conditions with respect to the Gröbner basis GB. We first need to define the interface conditions. To achieve this, we compute the Smith normal form of R with respect to the variable dx.

> S:=map(factor,smith(R,dx,'U','V'));

RR n°7953

We reduce them:

Update step The interface conditions are the following:

We reduce them:

Correction step The interface conditions are the following: We reduce them:

Oseen 3D

We consider the Oseen equations in R 3 . This system is defined by R y = 0, where R ∈ A 4×4 is a matrix with entries in the commutative polynomial ring A = Q(ν, c, b1, b2, b3)[dx, dy, dz] of PD operators in dx, dy, dz with coefficients in Q(ν, c, b1, b2), where ν is the viscosity, c the reaction coefficient, (b1, b2, b3) is convection velocity, dx = d/dx, dy = d/dy and dz = d/dz are the derivations with respect to x respectively y and z.

We first define A and R. The Smith normal form of R with respect to the variable dx is given by: > S:=map(factor,smith(R,dx,'U','V')); Inria