Axel Kr Öner
email: axel.kroener@oeaw.ac.at

Karl Kunisch
email: karl.kunisch@uni-graz.at

AND Hasnaa Zidani

OPTIMAL FEEDBACK CONTROL FOR UNDAMPED WAVE EQUATIONS BY SOLVING A HJB EQUATION

Keywords: optimal control, wave equation, Hamilton-Jacobi Bellman equation, spectral elements AMS subject classifications. 49J20, 35L05, 49J15

An optimal finite-time horizon feedback control problem for (semi-linear) wave equations is presented. The feedback law can be derived from the dynamic programming principle and requires to solve the evolutionary Hamilton-Jacobi-Bellman (HJB) equation. Classical discretization methods based on finite elements lead to approximated problems governed by ODEs in high dimensional spaces which makes the numerical resolution by the HJB approach infeasible. In the present paper, an approximation based on spectral elements is used to discretize the wave equation. The effect of noise is considered and numerical simulations are presented to show the relevance of the approach.

Introduction.

In this paper we consider optimal feedback control for onedimensional (semi-linear) wave equations. Let Ω ⊂ R 1 be an open interval, T > 0. We consider control problems of the following type for given distributed costs l : L 2 (Ω) × L 2 (Ω) → R, discount factor µ, set of admissible controls U , dynamics F given by a (semi-linear) wave equation, and initial point y 0 at time t = 0. The state y is defined as the solution of the wave equation under the control u. A precise definition is given in the next section.

Although there are many publications for numerical methods for open-loop optimal control problems with partial differential equations, there are only few results on numerical methods for closed-loop optimal control of partial differential equations which are conceived for nonlinear problems. Riccati-based methods and their numerical approximations for the linear regulator problem in the context of infinite dimensional systems have received a significant amount of attention in the past. The linear-quadratic regulator for infinite dimensional, second order (in time) linear oscillators is considered in [START_REF] Gibson | An analysis of optimal modal regulation: convergence and stability[END_REF]. In [START_REF] Banks | The linear regulator problem for parabolic systems[END_REF] an approximation framework for the computation (in finite dimensional spaces) of Riccati operators is presented that converges to the Riccati operator for the linear regulator problem for a large class of parabolic systems. In [START_REF]A numerical algorithm for optimal feedback gains in high dimensional linear quadratic regulator problems[END_REF] a hybrid method for computing feedback gains in linear quadratic regulator problems is analyzed which is in particular appropriate when used with large dimensional systems like control systems for partial differential equations and in [START_REF] Ito | On the regularity of solutions of an operator Riccati equation arising in linear quadratic optimal control problems for hereditary differential systems[END_REF] a regularity result for solutions of Riccati equation is shown which is essential to obtain a rate of convergence for numerical approximations. A review article on continuous and approximation theory for Riccati equations can be found in [START_REF] Lasiecka | Differential and Algebraic Riccati Equations with Applications to Boundary/Point Control Problems: Continuous Theory and Approximation Theory[END_REF].

For recent results on Riccati-based methods in the context of control of partial differential equations see, e.g., [START_REF] Mehrmann | Explicit solutions for a Riccati equation from transport theory[END_REF], where an explicit solution for Riccati equations arising in transport theory is considered. In [START_REF] Bänsch | Stabilization of incompressible flow problems by Riccati-based feedback, Constrained Optimization and Optimal Control for Partial Differential Equations[END_REF] the authors study a Riccati approach for the two dimensional Navier Stokes equation, and in [START_REF] Raymond | Feedback boundary stabilization of the two-dimensional Navier-Stokes equations[END_REF][START_REF] Buchot | Nonlinear feedback stabilization of a two dimensional Burgers equation[END_REF] the authors consider feedback stabilization for Burgers and Navier-Stokes equations.

For optimal feedback control of partial differential equations based on Hamilton-Jacobi Bellman (HJB) equations there exists only few publications. Here we mention some of them. In [START_REF] Ferretti | Internal approximation schemes for optimal control problems in Hilbert spaces[END_REF] internal approximation schemes for optimal control problems in Hilbert spaces are considered and conditions for convergence of the approximate value function are given. In [START_REF] Gombao | Hamilton-Jacobi equations for control problems of parabolic equations[END_REF] HJB equations for optimal control problems of semilinear parabolic equations are analyzed. For optimal feedback control of the Burgers equation using proper orthogonal decomposition (POD) we refer to [START_REF] Kunisch | HJB-POD based feedback design for the optimal control of evolution problems[END_REF][START_REF] Kunisch | POD-based feedback control of the Burgers equation by solving the evolutionary HJB equation[END_REF]. For an optimal feedback control approach of an advection-diffusion equation with an adaptive POD method see [START_REF] Alla | An adaptive POD approximation method for the control of advectiondiffusion equations[END_REF].

For existing results in the context of feedback control of the wave equation we refer to [START_REF] Gibson | An analysis of optimal modal regulation: convergence and stability[END_REF], for the strongly damped wave equation and the Timoshenko beam to [START_REF] Hernández | A locking-free scheme for the LQR control of a Timoshenko beam[END_REF][START_REF] Hernández | Numerical approximation of the LQR problem in a strongly damped wave equation[END_REF]. For open-loop control of the wave equation see [START_REF] Gerdts | Numerical optimal control of the wave equation: Optimal boundary control of a string to rest in finite time[END_REF][START_REF] Gugat | L ∞ -norm minimal control of the wave equation: on the weakness of the bang-bang principle[END_REF][START_REF] Kröner | A minimum effort optimal control problem for the wave equation[END_REF][START_REF] Kröner | Semi-smooth Newton methods for optimal control of the wave equation with control constraints[END_REF][START_REF] Kunisch | On time optimal control of the wave equation, its regularization and optimality system[END_REF] and the references cited therein.

In this paper we use a feedback law for the semi-discrete problem of (1.1) which is based on the solution of the corresponding HJB equation. We study the stabilizing effect of the feedback control with respect to noise effecting the dynamics or the observation for linear and nonlinear equations. For validation purposes we compare the results with a linearized Riccati approach. We observe that in the case of the linear wave equation the HJB and the Riccati approach lead to similar results. In case of a nonlinear equation the Riccati approach may fail while the HJB approach still stabilizes the system.

The discretization of partial differential equations in space usually leads to high dimensional dynamical systems and it is therefore very challenging to solve the corresponding HJB equation (so-called curse of dimensionality). In the last decades, several theoretical and numerical developments in HJB theory led to powerful and efficient numerical approaches that can be used for control problems up to 6-dimensional problems [START_REF] Altarovici | A general Hamilton-Jacobi framework for nonlinear state constrained control problems[END_REF][START_REF] Bokanowski | A discontinuous Galerkin scheme for front propagation with obstacles[END_REF][START_REF] Bokanowski | A general Hamilton-Jacobi framework for nonlinear state constrained control problems[END_REF][START_REF] Bokanowski | An efficient data structure and accurate scheme to solve front propagation problems[END_REF][START_REF] Falcone | Semi-Lagrangian Approximation Schemes for Linear and Hamilton-Jacobi Equations[END_REF][START_REF] Osher | Essentially nonoscillatory schemes for Hamilton-Jacobi equations[END_REF][START_REF] Sethian | Fast marching methods[END_REF]. For higher dimensional problems, various approaches have been studied in the literature, including model reduction or advanced numerical schemes as, e.g., sparse grids, see [START_REF] Bokanowski | An adaptive sparse grid semi-Lagrangian scheme for first order Hamilton-Jacobi Bellman Equations[END_REF]. Finally, we also point out that a general introduction for Hamilton-Jacobi Bellman equations and the numerical treatment can be found in the monograph [START_REF] Bardi | Optimal Control and Viscosity Solutions of Hamilton-Jacobi-Bellman Equations[END_REF].

In our approach we overcome the difficulty (curse of dimensionality) by discretizing the wave equation with spectral elements in space. In many situations already few basis functions are sufficient to obtain a useful approximation. This choice of basis functions leads to diagonal stiffness and mass matrices which allow in certain situations a decomposition of the originally problem in lower dimensional subproblems which can be solved in parallel.

The paper is organized as follows. In Section 2 we formulate the continuous control problem, in Section 3 we present a semi-discretization in space of (semi-linear) wave equations and formulate the DPP and the corresponding HJB equation, in Section 4 we describe the decomposition in subproblems in case of optimal control of the linear wave equation, in Section 5 we formulate the discrete scheme for solving the HJB equation, in Section 6 we formulate the numerical feedback law, and in Section 7 we present some numerical examples.

The control problem.

In this section we introduce the control problem in its functional analytic setting. Let

A : H 1 0 (Ω) × L 2 (Ω) → L 2 (Ω) × H -1 (Ω), A = 0 id c∆ 0
with the negative Laplacian (-∆) : H 1 0 (Ω) → H -1 (Ω), parameter c > 0, and identity map id : L 2 (Ω) → L 2 (Ω). We consider semi-linear wave equations (written as a first order system in time) of the following type (2.1)

   d ds y(s) -Ay(s) -g(y 1 (s)) = f (s), s ∈ (0, T ), y(0) = y 0 with solution y = (y 1 , y 2 ) T , t ∈ [0, T ), T > 0, f = (0, f (•)) T , f ∈ L 2 (L 2 (Ω)), initial point y 0 ∈ H 1 0 (Ω) × L 2 (Ω)
, and nonlinearity

g : L 2 (Ω) → { 0 } × L 2 (Ω), z → (0, g 2 (z)) T with g 2 (z) -g 2 (w) L 2 (Ω) ≤ L g z -w L 2 (Ω) ∀z, w ∈ L 2 (Ω) (2.2)
for Lipschitz constant L g > 0. We assume that there exists a unique bounded solution in C([0, T ]; H 1 0 (Ω) × L 2 (Ω)). Remark 2.1. For the linear wave equation (i.e. g ≡ 0) there exists a unique solution y ∈ C([0, T ]; H 1 0 (Ω) × L 2 (Ω)) (see [35, p. 275]). Existence and uniqueness for the semi-linear wave equation can be shown under further conditions on the nonlinear source term g, see [START_REF] Zuazua | Exact controllability for semilinear wave equations in one space dimension[END_REF]. The Lipschitz continuity assumption (2.2) of g will be used for the HJB approach in next section.

To set up a control problem we introduce the set of admissible controls

U = L 2 ((0, T ); U ) (2.3) with U = { u ∈ R m | u i a ≤ u i ≤ u i b for i = 1, . . . , m } (2.4) and u i a , u i b ∈ R, i = 1, . . . , m, m ∈ N.
Further, we define the operator

B : R m → { 0 } × L 2 (Ω), Bu = (0, ũ) T , ũ = m i=1 u i sin(πix) (2.5)
i.e. for every input u ∈ U , the function B is a control of the amplitudes of the basis modes. Now, consider the distributed cost given by:

l : L 2 (Ω) 2 × R N → R, l(z, u) = 1 2 Cz -y d 2 L 2 (Ω) + α 2 u 2 (2.6)
for parameter α > 0, C = (id, 0), and a given desired state y d ∈ L 2 (Ω), respectively.

Here, • denotes the Euclidean norm.

Then the control problem is defined as

(2.7)            min u∈U T 0 l(y(s), u(s))e -µs ds, d ds y(s) = F(y(s), u(s)), y(0) = y 0
for y 0 ∈ Y 1 , a discount factor µ ≥ 0, and operator F defined by

(2.8) F : Y 1 × L 2 (Ω) → L 2 (Ω) × L 2 (Ω), F(x, u) = Ax + Bu + g(x 1 ) with x = (x 1 , x 2 ) T ∈ Y 1 .
In the sequel, we assume that for each control u ∈ U , there exists a unique state y solution of (2.7).

3. Semi-discrete formulation and HJB equation. In this section a semidiscretization in space of the wave equation is presented. This allows to formulate the dynamic programming principle for the corresponding semi-discrete control problem and to derive a HJB equation whose solution is the value function of the semi-discrete problem.

The curse of dimensionality (i.e. the infeasibility to solve high dimensional HJB equations) requires to consider semi-discrete approximations with low dimensional discrete dynamical systems. This is the case when discretizing the wave equation by spectral elements in space. Then the obtained approximated control problem governed by ordinary differential equations has a reasonable number of state variables which allow the problem to be solved by a low dimensional HJB equation.

Discretization in space.

For the semi-discretization in space we make the following ansatz. Let the state be given by y h = (y h 1 , y h 2 ) with

y h 1 (x, t) = N k=1 y h 1,k (t)ϕ k (x), y h 2 (x, t) = N k=1 y h 2,k (t)ϕ k (x) (3.1)
for given basis functions ϕ k ∈ H 1 0 (Ω), k = 1, . . . , N , and set

y h := (y 1,1 , . . . , y 1,N , y 2,1 , . . . , y 2,N ) T ∈ L 2 ((0, T ); R 2N ).
As mentioned above it is important to choose a low order of approximation N , since the HJB equation that will be associated to the approximated control problem will be in a 2N -dimensional space, see Section 3.3 where we will specify the choice of the basis functions in more detail. Further, we define the stiffness and mass matrices (3.2) A h = ((∇ϕ i , ∇ϕ j ) L 2 (Ω) ) i=1,...,N,j=1,...,N , M h = ((ϕ i , ϕ j ) L 2 (Ω) ) i=1,...,N,j=1,...,N with respect to the basis ϕ i , i = 1, . . . , N . We introduce the discrete dynamics

F h : R 2N × R N → R 2N , F h (z, u) = M -1 h (A h z + B h u + g h (z)) (3.3) with matrices M h = M h 0 0 M h , A h = 0 M h -cA h 0 , B h = 0 M h ,
and nonlinearity

g h (z) = (0, g h,2 (z)) T , g h,2 (z) =   (g 2 ( N k=1 z k ϕ k ), ϕ i ) L 2 (Ω) i=1,...,N   .
and the discrete distributed cost

l h : R 2N × R N → R, l h (z, u) = 1 2 (C h z -y h d ) T M h (C h z -y h d ) + α 2 u T u (3.4) for α > 0, y h d ∈ R N
, where y h d is the vector containing the coordinates of y d with respect to the spectral basis functions, C h ∈ R 2N ×2N with C h = (id N , 0) T . Here, id N denotes the identity on R N . The Lipschitz continuity of g implies the Lipschitz continuity of g h . The semi-discrete control problem is given by

(3.5)            inf u∈U J h (u) = T 0 l h (y h
x,t (s), u(s))e -µs ds, s.t. d ds y h (s) = F h (y h (s), u(s)) for s ∈ (0, T ),

y h (0) = y h 0 ,
for initial value y h 0 ∈ R 2N . The initial state y h 0 is the projection of y 0 on the finite dimensional space V h × V h with V h = span(ϕ 1 , . . . , ϕ N ).

To simplify the notation we will henceforth omit the index h.

Dynamic programming principle and

Hamilton-Jacobi Bellman equation. For x ∈ R 2N and t ∈ [0, T ), we consider the parameterized control problem:

(P x,t ) inf u∈U T t l h (y x,t (s), u(s))e -µs ds, d ds y x,t (s) = F h (y x,t (s), u(s)) for s ∈ (t, T ), y x,t (t) = x.
In the above problem, the notation y x,t denotes the solution of the ODE system that starts in the position x at time t. To the family of problems (P x,t ), we consider the value function v : R 2N × [0, T ] → R that associates to each initial condition (x, t) the infimum value of the problem (P x,t ):

v(x, t) = inf(P x,t ).

The value function v satisfies the dynamic programming principle (DPP), see [6, p. 154-155]

(3.6) v(x, t) = inf u∈U τ t e µ(t-s) l h (y x,t (s), u(s))ds + e µ(t-τ ) v(y x,t (τ ), τ ) , for τ ∈ [t, T ], t ∈ [0, T ].
From the DPP we derive the value function v is given as the viscosity solution of the evolutionary Hamilton-Jacobi Bellman (HJB) equation, cf. [6, p. 155],

(3.7)

-∂ t v(x, t) + µv(x, t) + H(x, ∇ x v(x, t)) = 0 (x, t) ∈ R 2N × (0, T ), v(x, T ) = 0 x ∈ R 2N , with Hamiltonian (3.8)
H : R 2N × R 2N -→ R, H(x, p) = sup u∈U -F h (x, u) T p -l h (x, u) .
Moreover, from [START_REF] Ishii | Uniqueness of unbounded viscosity solution of Hamilton-Jacobi equations[END_REF] we know that v is the unique viscosity solution of (3.7) with quadratic growth. Finally, let us point out that by straightforward calculations, the Hamiltonian can be reformulated as:

(3.9) H(x, p) = -x T (M -1 h A h ) T p -(M -1 h g(x 1 )) T p -u(p) T (M -1 h B h ) T p -( 1 2 (C h x -y h d ) T M h (C h x -y h d )) - α 2 u(p) T u(p) with u(p) := argmax u∈U -F h (x, u) T p -l h (x, u) = P [ua,u b ] - 1 α M -1 h B T h p (3.10)
for any p ∈ R 2N , and projection P [ua,u b ] defined by

P [ua,u b ] (z) = ẑ, ẑi = min(u i b , max(u i a , z i ))
for z ∈ R 2N . In case of inactive control constraints this reduces to

H(x, p) = -x T (M -1 h A h ) T p -(M -1 h g(x 1 )) T p -( 1 2 (C h x -y h d ) T M h (C h x -y h d )) + 1 2α p T BB T p with B = M -1 h B h . 3.3. Spectral elements.
In general it is very challenging to solve the HJB (3.7) in high dimensions. To obtain a low dimensional problem we use spectral elements for the discretization of the wave equation in space. This allows to use only few basis functions to obtain an appropriate approximation. Spectral elements play an important role and are widely used for the numerical discretization of partial differential equations, see, e.g. [START_REF] Canuto | Spectral Methods: Evolution to Complex Geometries and Applications to Fluid Dynamics[END_REF][START_REF] Gottlieb | Numerical Analysis of Spectral Methods: Theory and Applications[END_REF][START_REF] Patera | A spectral element method for fluid dynamics: Laminar flow in a channel expansion[END_REF] and in the context of optimal control [START_REF] Gaudio | Spectral element discretization of optimal control problems, Spectral and High Order Methods for Partial Differential Equations[END_REF]. To discretize the wave equation by spectral elements we specify the domain Ω = (0, 1) and choose the basis functions by

ϕ k (x) = sin(kπx), k = 1, . . . , N, x ∈ Ω.
The corresponding mass and stiffness matrix are diagonal and are given by

M h = diag (1/2) k=1,...,N , A h = diag 1 2 (kπ) 2 k=1,...,N .
The following numerical example shows, that for the corresponding open-loop control problem only few spectral basis functions are necessary to obtain a good approximation.

Example 3.1. We consider problem (2.7) and allow that the distributed costs also depend on time. We choose T = 1, Ω = (0, 1), µ = 0, c = 1, and the distributed costs as

l : L 2 (Ω) × L 2 (Ω) × [0, T ] → R, l(z, u, t) = 1 2 z -y d (t) 2 L 2 (Ω) + α 2 u 2 L 2 (Ω)
with α = 1, and the desired state by

y d (z, t) = (t -T ) 4 z 3 (1 -z) 3 + 24z 3 (1 -z) 3 -12(t -T ) 2 (6z(1 -z) 3 -18z 2 (1 -z) 2 + 6z 3 (1 -z)) -12(t -T ) 2 (6z(1 -z) 3 -18z 2 (1 -z) 2 + 6z 3 (1 -z)) + (t -T ) 4 -18(1 -z) 2 + 216z(1 -z) -54(1 -z) 2 -72z 2 .
The initial values are given by

x(z) = (T 4 z 3 (1 -z) 3 , -4T 3 z 3 (1 -z) 3 ) T .
The constraints on the control are chosen in such a way that they remain inactive.

Then the exact state is given by

y * (z, t) = ((t -T ) 4 z 3 (1 -z) 3 , 4(t -T ) 3 z 3 (1 -z) 3 ) T .
For N=5 spectral basis functions the error between the exact solution y * and the nu-

merical approximation ỹ is ỹ -y * L 2 (L 2 (Ω)) ≈ 1.29 • 10 -4 .
For time discretization a Crank Nicolson scheme is applied with temporal mesh parameter k = 1.0 • 10 -4 . This example shows that already few spatial basis functions lead to a good approximation of the state.

Remark 3.2. As already known and mentioned earlier in the paper, the theory of the HJB equation approach is valid in any dimension. However, in practice, the implementation of the method in high dimension requires, in general, a huge numerical effort (data storage and computing time).

Several works have dealt with model reduction of optimal control problems governed by partial differential equations leading to a reduced controlled system that can be studied numerically by using the HJB approach, see, e.g., [START_REF] Kunisch | POD-based feedback control of the Burgers equation by solving the evolutionary HJB equation[END_REF] and [START_REF] Kunisch | HJB-POD based feedback design for the optimal control of evolution problems[END_REF], where proper orthogonal decomposition (POD) is applied to reduce the dimension of the underlying dynamical system, and in [START_REF] Alla | An adaptive POD approximation method for the control of advectiondiffusion equations[END_REF], where for adaptively chosen small time horizons the dynamical system is reduced by POD.

Note also that, it would be interesting to reduce the dimension not in the definition of the control problem, but in the HJB equation associated to it. However, for such a nonlinear PDE the analysis of reduction techniques is still an open problem, see [START_REF] Bokanowski | An adaptive sparse grid semi-Lagrangian scheme for first order Hamilton-Jacobi Bellman Equations[END_REF] for a method based on sparse grids. 4. Decomposition in two dimensional subproblems. The discretization by spectral elements leads to diagonal mass and stiffness matrices, i.e. the set of basis functions is an orthogonal basis with respect to the L 2 (Ω)-and H 1 0 (Ω)-inner product. In certain situations this allows to decompose the 2N -dimensional HJB equation in N 2-dimensional decoupled HJB equations which can be solved in parallel. In this section we assume that the problem is of linear-quadratic type with possibly constraints on the control, i.e.

           min u∈U J h (x, t, u) = T t l h (y x,t (s), u(s))e -µs ds, s.t. d ds y x,t = M -1 h (A h y x,t + B h u) , y x,t (t) = x.
We introduce N subproblems defined, for i = 1, . . . , N by (4.2)

         v i (ẑ, t) = min u(•)∈Ui 1 2 T t 1 2 z 1 i (s) 2 + αu 2 (s) e -µs ds, żi (s) = A i z i (s) + Bu(s), z i (t) = ẑ, for z i = (z 1 i , z 2 i ) T and with U i = L 2 ((0, T ); U i ) for U i = { u ∈ R | u i a ≤ u ≤ u i b } and A i = 0 1 a i 0 , B = 0 1 ,
and where

a i = (M -1 h A h ) ii . For every i = 1, • • • , N , the value functions v i of the subproblem (4.
2) satisfies the following HJB equation:

(4.3) -∂ t v i (z, t) + µv i (z, t) + H i (z, ∇ z v i (z, t)) = 0 (z, t) ∈ R 2 × (0, T ), v i (z, T ) = 0 z ∈ R 2 with H i (z, p) := sup u∈Ui (-(A i z + Bu) T p -1 2 ( 1 2 z 2 + αu 2 )) for z, p ∈ R 2 .
There holds the following relation between the value function of problem (4.1) and the value functions of the subproblems (4.2).

Theorem 4.1. Let v i be the solution of the subproblem (4.2) for i = 1, . . . , N and let v be the solution of (4.1). Then there holds the relation

v(x, t) = N i=1 v i x i x i+N , t for every t ∈ (0, T ), x =    x 1 . . . x 2N    ∈ R 2N .
Proof. The arguments are straightforward since in the present case (when g ≡ 0) the controlled system can be split into N decoupled systems. Also the minimisation over the control set U can be split into N independent minimisation problems. The above result means that instead of solving (3.7), one can solve N independent 2-dimensional HJB equations (4.3) with different data. These computations can be done in parallel, which has as consequence a significant reduction of the computational effort (in term of data storage and of cpu time).

Remark 4.2. It is well-known that in general the stiffness matrix for spectral elements is ill-conditioned. By the decomposition the drawback of this property is reduced since in every subproblem only two dimensional matrices are considered.

5. Discretization of the HJB. Several discretization schemes for HJB equations have been studied in the literature. For an overview on finite differences schemes, we refer to [START_REF] Osher | Essentially nonoscillatory schemes for Hamilton-Jacobi equations[END_REF]. Another recent reference [START_REF] Falcone | Semi-Lagrangian Approximation Schemes for Linear and Hamilton-Jacobi Equations[END_REF] gives an overview on semi-Lagrangian schemes and the features of this class of methods. For a finite element approach we refer to [START_REF] Jensen | On the convergence of finite element methods for Hamilton-Jacobi-Bellman Equations[END_REF], and for discontinuous Galerkin methods to [START_REF] Bokanowski | A discontinuous Galerkin scheme for front propagation with obstacles[END_REF]. Another class of antidiffusive methods have also been studied in [START_REF] Bokanowski | An efficient data structure and accurate scheme to solve front propagation problems[END_REF][START_REF] Bokanowski | Convergence of a non-monotone scheme for hamilton-jacobi-bellman equations with discontinuous initial data[END_REF]. Here, we use a finite difference scheme method, more precisely, an ENO scheme for space discretization is coupled with a Runge-Kutta time discretization scheme of second order.

5.1. Numerical scheme.. Following [START_REF] Osher | Essentially nonoscillatory schemes for Hamilton-Jacobi equations[END_REF][START_REF] Bokanowski | A general Hamilton-Jacobi framework for nonlinear state constrained control problems[END_REF] we shortly recall the scheme for a given continuous Hamiltonian H : R 2N × R 2N → R. For given temporal mesh parameter ∆t and spatial mesh parameter ∆x = (∆x 1 , . . . , ∆x 2N ) we define a spatial mesh

G = { x I = (I 1 ∆x 1 , • • • , I 2N ∆x 2N ) for I ∈ Z 2N }
and temporal mesh

t = t 0 < t 1 < • • • < t M = T, t j+1 -t j = ∆t, j = 0, . . . , M T -1, M T ∈ N.
Further, we introduce difference quotients for x ∈ G by

D + v(x) = (D + x1 v(x), . . . , D + x 2N v(x)), D -v(x) = (D - x1 v(x), . . . , D - x 2N v(x)) with D ± xi v(x I ) = ± v(x I i,± ) -v(x I ) ∆x and I k,± = (i 1 , . . . , i k-1 , i k ± 1, i k+1 , . . . , i 2N ).
Then the Lax-Friedrichs scheme for the HJB equation reads as follows: v M I = 0 for every I ∈ Z 2N , and

v n-1 I = v n I -∆tµv n I -∆tH LF (x I , D + v n (x I ), D -v n (x I )) for n = M, • • • , 1
with Lax-Friedrichs Hamiltonian

H LF (x, p + , p -) = H(x, p + + p - 2 ) - 2N i=1 C i (x) 2 (p + i -p - i ) (5.1)
and the stabilizing functions

C i (x) satisfying max p∈R 2N ∂H ∂p i (x, p) ≤ C i (x)
for i = 1, . . . , 2N . It is well-known that the Lax-Friedrichs scheme is monotone and the convergence is of first order as far as the CFL condition

k max x∈G 2N i=1 C i (x) ∆x i ≤ 1 is satisfied.
An ENO scheme can be obtained by considering a variant of the LF scheme

v n-1 I = v n I -∆tµv n I -∆tH LF (x I , D + v n (x I ), D -v n (x I )),
where D ± v n (x I ) are higher approximations of the gradient ∂ xi v coupled with a Runge-Kutta time discretization scheme of second order, see [START_REF] Bokanowski | A general Hamilton-Jacobi framework for nonlinear state constrained control problems[END_REF][START_REF] Osher | Essentially nonoscillatory schemes for Hamilton-Jacobi equations[END_REF]. Convergence of this method is not proved but various numerical examples confirm its relevance.

5.2.

Boundary condition for HJB equation on the computational domain. Since the HJB is defined on the full space we have to define a bounded computational domain K × (0, T ), with K ⊂ R 2N , to solve the problem numerically. For simplicity K will be chosen as

K := Π 2N i=1 [c i , d i ],
and the considered grid is as follows:

G 0 = { x I = (c 1 + (I 1 -1)∆x 1 , • • • , c 2N + (I 2N -1)∆x 2N ) for I ∈ Π 2N i=1 [1, M i ] } ,
with ∆x i := di-ci Mi-1 . However, due the boundedness of K an appropriate choice of a boundary condition is required to get a well define numerical scheme.

When a strong invariance condition

F h (x, u) • n(x) < 0 ∀x ∈ ∂K, ∀u ∈ U,
(with exterior normal n to the boundary of K at x) is satisfied then all trajectories starting in D will remain within Ω forever, see [6, p. 485]. In this case, there is no need to impose a boundary condition. However, the strong invariance condition is very restrictive and cannot be insured even in the simple case when g ≡ 0.

Here, we introduce an artificial boundary condition by setting the second derivative on the boundary to zero:

∂ xx v(x, t) = 0 ∀(x, t) ∈ ∂Ω × (0, T ).
This means that at a time step t n , for n = M T , • • • , 1, if x I ∈ ∂Ω, then there exists i such that I i = 1 or I i = M i . In this case the upwind derivatives on the boundary should be defined as:

     D ± xi v(x I ) = v(x I i,+ ) -v(x I ) h i if I i = 0; D ± xi v(x I )) = - v(x I i,-) -v(x I ) h i if I i = M i .
Since the value function is continuous, this approximation is reasonable. For a more sophisticated treatment of the boundary condition using state constraints we refer to [START_REF] Altarovici | A general Hamilton-Jacobi framework for nonlinear state constrained control problems[END_REF].

6. Reconstruction of the trajectory. For the reconstruction of the trajectories we formulate a feedback law based on the HJB approach and present different numerical realizations. This approach is applicable on a large class of (nonlinear) problems without modifications of the general setting.

Notice that besides the HJB approach, to study optimal feedback control problems for nonlinear dynamical systems, it is possible in some situations to use a Riccati based approach. In this case, a Riccati feedback operator can be derived from the linearized state equation. This approach has been considered, for instance, in [START_REF] Raymond | Feedback boundary stabilization of the two-dimensional Navier-Stokes equations[END_REF] to study boundary control problems governed by Navier Stokes equation. In [START_REF] Buchot | Nonlinear feedback stabilization of a two dimensional Burgers equation[END_REF] the authors consider an infinite horizon problem and propose a method using formally a quadratic Taylor approximation of the solution of the HJB equation. For a non-autonomous system that appears in the internal stabilization of Navier Stokes equation to a nonstationary solution a feedback law satisfying a differential Riccati equation can be found in [START_REF] Barbu | Internal exponential stabilization to a nonstationary solution for 3d Navier-Stokes equation[END_REF].

According to the HJB approach and its associated dynamic programming principle, the optimal feedback control u * x for state x at time t ∈ [0, T ) satisfies

u * x = argmin u∈U [t,t+τ ]
e -µτ v (y x,t (t + τ ), t + τ )) + t+τ t e µ(t-s) l h (y x,t (s), u(s))ds (6.1) for τ ≥ 0 with t + τ ≤ T , see [6, p. 149/155], and

U [t,t+τ ] = { u : [t, t + τ ] → R N | ∃w ∈ U : w| [t,t+τ ] = u } .
The feedback operator K is defined as (6.2)

K : R 2N × (t, T ) -→ U, (x, τ ) → K(x, τ ) = u * x (τ ).
Once the feedback operator is defined, to reconstruct the optimal trajectory starting in x 0 at the initial time 0, one needs just to solve the dynamical system y ′ (s) = F h (y(s), K(y(s), s)), y(0) = x 0 , s ∈ (0, T ].

To realize the feedback law (6.2) numerically there are several methods. We give an overview about three of them. Let v ∆ be the discrete value function computed by on G 0 by the ENO scheme described in Section 5.

1. The dynamics and distributed cost are given in a separable form, i.e.

F(z, u) = F (z) + Bu, l(z, u) = j(z) + α 2 u T u with F : R 2N → R 2N , matrix B ∈ R 2N ×N ,

K(x, t)

= P [ua,u b ] - 1 α B T ∇ x v(x, t) . (6.3)
Since only an approximation v ∆ of v can be computed, an approximation of the feedback can be given by:

K ∆ (x, t) = P [ua,u b ] - 1 α B T ∇ x v ∆ (x, t) , (6.4)
where the gradient ∇ x v ∆ can be approximated component-by-component by a central difference quotient. 2. The control minimizing the Hamiltonian H can be saved at every time point and in every mesh point of G 0 when solving the HJB equation. For all x I ∈ G 0 , k ∈ { 0, . . . , M T }, we set

K ∆ (x I , t k ) = argmax u∈U (-F h (x I , u) T ∇ x v ∆ (x I , t k ) + l h (x I , u)),
where v ∆ denotes the discrete value function. This leads to a set of discrete controls { u I } I∈Z 2N ,x I ∈K whose elements correspond to the mesh points. The advantage of this approach is that it gives directly the set of discrete controls without any additional effort. A drawback of this approach is that this set has to be saved and could be very large (as big as the set of all stored values of the value function). Furthermore, to compute a feedback control which corresponds to a state which is not a mesh point in G 0 one has to compute the control by an interpolation method. This procedure requires the feedback control operator K to be continuous, which usually is not the case. 3. In this approach the dynamic programming principle is directly used to reconstruct the trajectories, i.e. using a discrete version of (6.1). We choose the feedback control by (6.5)

u * x (t k ) = argmin u∈U [v ∆ ](x + ∆tF h (y x,t k (t k ), u)) + ∆tl h (y x I ,t k (t k ), u) , =: argmin u∈U G(x, y x,t k (t k ), u).
Here, [•] is an interpolation of the value function v ∆ on G 0 . A broadly used method to determine the argmin consists in choosing a finite set U finite ⊂ U and computing the argmin with respect to this set by evaluating the expression G(x I , y x,t k (t k ), u) for all u ∈ U finite and taking the control value which minimizes the expression, see [6, p. 475]. If the state space is low dimensional, the approach works well. But in case of higher dimensions or a large control set U finite the approach is only feasible in combination with a strategy which reduces the number of evaluations of G(x I , y x,t k (t k ), u) in each time step. If the problem can be decomposed as described in Section 4 the number of evaluations reduces strongly. In all the simulations performed in this paper, we apply the first of these three methods assuming that the value function is sufficiently smooth. This method gives us a direct formula for the feedback control and saves much computational time in comparison to the third option.

Numerical examples.

Several numerical examples are presented studying the stabilizing effect of the feedback control. For solving the HJB equation we use the software package ROC-HJ-Solver, see [START_REF] Bokanowski | A C++ library for solving HJ equations[END_REF], which is parallized with OpenMP. For the reconstruction of the feedback controls and the corresponding trajectories MATLAB is used.

In the first three examples we consider an optimal feedback control problem for the linear wave equation and in the fourth one an optimal control problem for a semi-linear wave equation. In all examples we set Ω = (0, 1), µ = 0, y d = 0. and α = 0.1. To validate the numerical algorithm we compare the trajectory we obtain via the HJB approach with the trajectory via the Riccati approach. The differential Riccati equation for the problem is given by

-P ′ (t) = A T P (t) + P (t)A -P (t)BR -1 B T P (t) + Q, P (T ) = 0 (7.2) with Q = M h , R = α id N , A = M -1
h A h , and B = M -1 h B h . The feedback operator for the Riccati approach which maps the optimal state at time t namely y * (t) to the corresponding optimal control u * (t) is given by

K(t) = -R -1 B T P (t).
For the reconstruction of the trajectories via the HJB approach we choose the first method presented in Section 6. The computational domain for the HJB equation is chosen to be K = [-1, 1] 10 . For the temporal discretization we choose the mesh size ∆t = 2 • 10 -3 . For the spatial discretization we choose a mesh uniform in every spatial direction with the following numbers of discretization points 10 × 16 × 30 × 50 × 60 × 10 × 16 × 30 × 50 × 60.

That means for the spatial directions which correspond to higher modes we use more grid points. Further, we decompose the equation in subproblems as described in Section 4. For the discretization of the dynamical system as well as the Riccati equation ( 7.2) we apply a explicit Euler scheme and set the temporal mesh parameter to k ds = 10∆t = 0.02.

In Figure 7.1 we see the first and fifth component of the state y 1 and velocity y 2 obtained by the HJB (blue curves) and the Riccati (red curves) approach and the corresponding components of the feedback control of these two approaches. The other components show a corresponding behaviour, therefore we skip the plots here. For the chosen discretization both approaches lead to nearly the same results. If we choose finer spatial and temporal mesh sizes when solving the HJB equation the trajectory of the HJB approach approximates the trajectory of the Riccati approach more closely. Further numerical experiments show that for a fixed and sufficiently small temporal mesh size the discretization of the subproblems for the higher modes require a finer spatial mesh to obtain a good approximation. The corresponding state and velocity of the wave equation which correspond to the computed trajectories are shown in Figure 7.2. The example confirms that the chosen boundary condition described in Section 5.2 leads to appropriate numerical results.

Linear wave equation -with noise.

We study the stabilizing effect of the feedback control on the trajectories in presence of random noise which either effects the dynamics or the observation. In this example we choose the parameter α = 0.01 which will lead to a stronger decrease of the trajectories to zero than in the previous example.

7.2.1.

Noise in the dynamics. The noise enters the dynamics, i.e. we consider

y t = F h (y, u) + M -1 h B h δ, δ = (δ 1 , . . . , δ 5 ) T (7.3) with random functions δ i : [0, T ] → [δ a , δ b ], i = 1, . . . , 5, for δ a , δ b ∈ R.
The value function is, of course, computed for the dynamical system without noise and the feedback law is applied to the state resulting from the perturbed system. For data as in (7.1) we compute the numerical solution of (7.3) for noise scaled to the following interval (7.4)

δ i (s) ∈ [-4, 4], i = 1, . . . , 5, s ∈ (0, T ].
For solving the HJB equation we proceed as in the previous example. The components of the corresponding trajectory and feedback control are presented in Figure 7.3. Again, the blue curves are the components of the trajectory computed via the HJB approach (the corresponding Riccati approach leads to quite similar results, why we omit the corresponding curves here). The red curves we obtain by the Riccati approach applied to the equation without noise. If we plug the resulting control in the perturbed system we obtain the green curves. In contrast to the latter control we see that the feedback control stabilizes the dynamical system. The components y 2,i , i = 1, . . . , 5, are more susceptible to noise than y 1,i , since the noise enters the corresponding right hand sides. In Figure 7.4(a) and 7.4(b) the state and velocity for the HJB and Riccati approach are shown, respectively. Both approaches drive the state and velocity to zero. Components of the trajectory and the feedback control via the HJB (blue) approach for the problem with noise and via the Riccati approach (red) for the equation without noise as well as the components of the trajectory (green) when applying the optimal control of the noise free system to the perturbed system, (α = 0.01).

We compare these results with the state and velocity when applying the optimal control for the noise free dynamics to the perturbed system, see Figure 7.4(c). In this case the control fails to drive state and velocity to zero.

7.2.2.

Noise in the observation. The noise enters the observation, i.e. we consider the disturbed feedback law u(s) = K ∆ (x(s) + δ(s)e 10 , s), e 10 = (1, . . . , 1) T with noise δ(s) : (0, T ] → R. In Figure 7.5 we see the components of the trajectory and the control. Although we apply the feedback operator on the perturbed state the resulting control stabilizes the trajectories. Numerical examples show that the amplitude of the oscillations increase with the amplitude of the noise. We have to choose the noise level smaller than in the first example namely δ(s) ∈ [0.01, 0.01]. The noise has a much stronger effect on the trajectory than in the previous examples. As in the previous example the HJB and Riccati approach for the system with noise lead to very similar curves. 7.3. Linear wave equation -control of high order system. In Example 3.1 we saw that we can expect a good approximation when discretizing the wave equation with 5 basis functions as we did in the previous examples. However, when considering nondiagonal operators B h or nonlinear problems the decomposition in subproblems is no longer possible since the equations of the corresponding dynamical system are coupled. In these cases we restrict the discretization to 3 basis functions for the state and for the velocity, because of the curse of dimensionality, which results in a six dimensional dynamical system. Nevertheless we can use the feedback control of this low order system to control a corresponding high order system. In the following we present an example, where we see the effects of the feedback control obtained from a low order system on a high order one.

For the high order system we use a finite element discretization. For given τ > 0 we define the spatial mesh

x i = iτ,
with i = 0, . . . , q, q = 1/τ and consider a standard linear finite element space

V H τ = { φ ∈ C 0 (Ω) | φ| Iq ∈ P 1 (I q ), I q = [
iτ, (i + 1)τ ], i = 0, . . . , q -1 } , where P 1 (I q ) consist of linear functions on I q . Further we define the stiffness matrix A FEM and mass matrix M FEM as in (3.2) with respect to the finite element basis and the corresponding discrete dynamics F FEM as in (3.3) with g h ≡ 0. Then the semi-discrete system is given by Components of the trajectory and the feedback control via the HJB approach (blue) for the system with noise and via the Riccati approach (red) for the equation without noise as well as the components of the trajectory (green) when applying the optimal control of the noise free system to the perturbed system, (α = 0.01) with noise δ : [0, T ] → R, e ∈ R q-1 , e(i) = 1 for i = 1, . . . , q -1, and initial point x ∈ R q-1 . For the time discretization we apply a Crank-Nicolson scheme and use the same temporal mesh as for the low order system.

To set up an algorithm to compute the feedback control for the full order system we introduce the projection on the low dimensional space

P L : V H τ → R 3 , φ → P L φ with (P L φ) i = 2 Ω φ(x) sin(iπx)dx, i = 1, 2, 3
and the pointwise evaluation operator

E : C(Ω) → V H τ , E(φ) = φ, φi = φ(x i ), i = 1, . . . , q -1, φ ∈ C(Ω). (7.6)
The feedback rule for the full order system is presented in Algorithm 1.

For the numerical example we choose the data as follows c = 0.025, T = 8, y 0 (z) = (1.5 sin(πz) + sin(2πz) + sin(3πz), 0) T , z ∈ Ω (7.7) with α = 0.1. We solve the HJB equation on K = [-2.5, 2.5] 6 with ∆t = 2 • 10 -3 . For the discretization in space we choose a spatial mesh with the following numbers of discretization points

4 × 8 × 20 × 4 × 8 × 20.
For the discretization of the wave equation we set the temporal step size to k ds = 0.02 and τ = 1/32.

In Figure 7.6 we compare the components of the trajectory and control for the HJB and Riccati approach. For the chosen discretization of the HJB equation we observe very similar results for both approaches. In Figure 7.7 we see the optimal Algorithm 1 Feedback control for full order system 1: Set n = 0. 2: Initial data:

y L 0 = P L (y H 0 ) for initial data y H 0 ∈ V H τ . 3: Repeat 4:
Compute the low order feedback control u L n = K ∆ (y L n , nk ds ).

5:

Compute the high order feedback control

u H n = E 3 i=1 u L n,i ϕ i . 6:
Time step for high order system

y H n+1 = y H n + k ds 2 (F FEM τ (y H n+1 , u H n ) + F FEM τ (y H n , u H n )).

7:

Project y H n+1 in the low dimensional space y L n+1 = P L (y H n+1 ).

8:

n = n + 1. 9: Until nk ds = T . state and velocity of the low order system and the full order system without noise. As one can expect the results are the same. For the case with noise δ(s) ∈ [-10, 10], s ∈ (0, T ], we see in Figure 7.8 the state and velocity for the uncontrolled full order system and when applying the feedback control. In contrast to the uncontrolled case the state and velocity are stabilized by the feedback control.

If the initial data of the full order problem has components not lying in the low dimensional subspace they cannot be controlled as we see in Figure 7.9 where we choose the initial data as y 0 (z) = (1.5 sin(πz) + sin(2πz) + sin(3πx) + sin(4πz), 0) T , z ∈ Ω.

We observe that the component of the solution with respect to the fourth basis function remains in the system.

Next, let the control act only on an open subset ω of Ω, i.e.

B ω = χ ω B with characteristic function χ ω of ω. We choose ω = (0, 0.75) and obtain a discrete operator given by This leads to a dynamics given by The functions p i (x) = χ ω sin(iπx), i = 1, 2, . . . , are not eigenfunctions of the Laplacian. Therefore, since in every time step we project the state in the low dimensional space there will always remain components which cannot be controlled as wee see in Figure 7.10 where the state, velocity and control for the low and the full order system are shown for data as given (7.7).

F h (x, u) = M -1 h (A h x + B ω h u) .
7.4. Semilinear wave equation -with noise. In this example we consider an optimal control problem for a semilinear wave equation. Let the dynamics F be   with y h 1 = (y 1,1 , y 1,2 , y 1,3 ) T . We compare the results with a linearized Riccati approach. We solve the differential Riccati equation for the linearized state equation (cf. Section 6), i.e. we consider z t -M -1 h A h z -M -1 h 3ȳ 2 z = M -1 h B h u. and linearize around the target ȳ = 0.

For the discretization of the HJB equation we choose the computational domain as We consider two different initial values ỹ0 (z) = 1.8 sin(πz) + 0.8 sin(2πz) + 0.5 sin(3πz) 1.8 sin(πz) + 0.5 sin(2πz) + 0.8 sin(3πz) , ŷ0 (z) = 1.5 sin(πz) + 0.5 sin(2πz) + 0.5 sin(3πz) 1.5 sin(πz) for z ∈ Ω.

7.4.1. Stabilization via HJB approach for initial data ỹ0 (linearized Riccati approach fails). The linearized Riccati approach fails to stabilize the system for initial data ỹ0 . It leads to a blow up in finite time as we can see in Figure 7.11(a). This is due to the fact that the initial data is not in the set [-1, 1] 6 , which leads to a stronger dominance of the nonlinear term. Nevertheless, the HJB approach stabilizes the system, see For noise δ i : [0, T ] → [-3, 3], i = 1, 2, 3, effecting the dynamics as

y t = F h (y, u) + M -1
h B h δ, δ = (δ 1 , δ 2 , δ 3 ) T (7.9) Figure 7.12 shows the components of the trajectory. 7.4.2. Stabilization via HJB approach for initial data ŷ0 (linearized Riccati approach works). Finally, we consider the problem with noise in the dynamics for initial point ŷ0 . The components of the trajectory are shown in Figure 7.13. For t > 2.5 the state is already close to zero and in particular smaller than 1 which implies that the linear term becomes dominant and the nonlinear very small. The influence of the nonlinear term is hence very small and thus the two feedback laws lead to similar curves. Thus in this case the initial point is sufficiently close to zero so that the Riccati approach still works.

  (s), u(s))e -µs ds, d ds y(s) = F(y(s), u(s)), y(0) = y 0

  and j : U → R as it is the case in the problem under consideration. If the value function satisfy v ∈ C 1 (R 2N ) the feedback operator is determined by the projection of the gradient of the value function on the set of admissible controls (cf. (3.10)), i.e.

7. 1 .

 1 Linear wave equation -without noise. In this example we consider the control problem under consideration for the linear wave equation with data c = 0.025,

5 Fig. 7 . 1 :

 571 Fig. 7.1: Components of the trajectory and the feedback control via HJB (blue) and Riccati (red) approach, (α = 0.1)

Fig. 7 . 2 :

 72 Fig. 7.2: State and velocity of the wave equation by the HJB and Riccati approach, (α = 0.1)

5 Fig. 7 . 3 :

 573 Fig.7.3: Noise in the dynamics: Components of the trajectory and the feedback control via the HJB (blue) approach for the problem with noise and via the Riccati approach (red) for the equation without noise as well as the components of the trajectory (green) when applying the optimal control of the noise free system to the perturbed system, (α = 0.01).

Fig. 7 . 4 :

 74 Fig. 7.4: Noise in the dynamics: State and velocity of the wave equation via the HJB andRiccati approach as well as the state and velocity of a noisy dynamical system controlled by the optimal control of the noise free problem, (α = 0.01).

5 Fig. 7 . 5 :

 575 Fig.7.5: Noise in the observation: Components of the trajectory and the feedback control via the HJB approach (blue) for the system with noise and via the Riccati approach (red) for the equation without noise as well as the components of the trajectory (green) when applying the optimal control of the noise free system to the perturbed system, (α = 0.01)

3 Fig. 7 . 6 :Fig. 7 . 7 :

 37677 Fig. 7.6: Components of the trajectory and the control of the HJB (blue curves) and Riccati (red curves) approach

Fig. 7 . 8 :

 78 Fig. 7.8: Uncontrolled and controlled state and velocity for full order system with noise in the dynamics

Fig. 7 . 9 :y 3 )

 793 Fig.7.9: State and velocity of the full order system with initial data x(z) = (1.5 sin(πz) + sin(2πz) + sin(4πz), 0) T not lying in the low order subspace

Fig. 7 . 10 :

 710 Fig. 7.10: Local control of the low and high order system

[- 2 , 2 ]

 22 × [-1, 1] × [-1, 1] × [-4, 4] × [-1, 1] × [-1, 1]and the numbers of discretization points as 10 × 8 × 16 × 10 × 8 × 16.For time discretization we set ∆t = 3.2 • 10 -3 . For the discretization of the wave equation we choose k ds = 3.2 • 10 -2 .

Figure 7 .

 7 11(b).

Fig. 7 . 11 :

 711 Fig.7.11: Components of the trajectory for the linearized Riccati approach and HJB approach for nonlinear dynamics with initial point y0 = ỹ0.

The work of the second author was supported in part by the Fonds zur Förderung der wissenschaftlichen Forschung under SFB 32, Mathematical Optimization and Applications in the Biomedical Sciences. (Hasnaa.ZidaniThe work of the third author was partially supported by the EU under the 7th Framework Program Marie Curie Initial Training Network FP7-PEOPLE-2010-ITN, SADCO project, GA number 264735-SADCO.

1


Paul Feautrier

Éric Violard

Alain Ketterlin

Improving X10 Program Performances by Clock Removal

X10 is a promising recent parallel language designed specifically to address the challenges of productively programming a wide variety of target platforms. The sequential core of X10 is an object-oriented language in the Java family. This core is augmented by a few parallel constructs that create activities as a generalization of the well known fork/join model. Clocks are a generalization of the familiar barriers. Synchronization on a clock is specified by the advance() method call. Activities that execute advances stall until all existent activities have done the same, and then are released at the same (logical) time. This naturally raises the following question: are clocks strictly necessary for X10 programs? Surprisingly enough, the answer is no, at least for sufficiently regular programs. One assigns a date to each operation, denoting the number of advances that the activity has executed before the operation. Operations with the same date constitute a front, fronts are executed sequentially in order of increasing dates, while operations in a front are executed in parallel if possible. Depending on the nature of the program, this may entail some overhead, which can be reduced to zero for polyhedral programs. We show by experiments that, at least for the current X10 runtime, this transformation usually improves the performance of our benchmarks. Besides its theoretical interest, this transformation may be of interest for simplifying a compiler or runtime library.

Introduction

Due to physical limitations, today computers all have explicit parallelism. This is true over the whole power spectrum, from embedded systems to high performance number crunchers, in which millions of cores must contribute to a common task. Efficient programming of such architectures is one of the most important challenge of the next decade. Among the many solutions which have been proposed -parallel programming libraries, domain specific languages, automatic parallelization -one of the most interesting is the use of parallel programming languages: languages in which parallel constructs are first class citizens, on a par with standard control constructs like the sequence or the loop. This approach has two advantages. Firstly, it hides the intricate details of parallel programming at the hardware or operating system level. Second, and most importantly, the programmer can express the problem inherent parallelism. Such parallelism might be difficult to infer from a sequential implementation.

The recent years have seen the creation of many such languages, among which Titanium [START_REF] Yelick | Titanium: A high-performance Java dialect[END_REF], Chapel [START_REF] Chamberlain | Parallel programmability and the Chapel language[END_REF], Co-Array Fortran [START_REF] Numrich | Co-array Fortran for parallel programming[END_REF], UPC [START_REF] Consortium | UPC language specifications[END_REF], Habanero Java [START_REF] Cavé | Habanero-java: the new adventures of old X10[END_REF]. This paper deals with X103 which is being developed at IBM Research. However, we believe that our techniques -if not our results -can be adapted without difficulties to other languages. Basically, parallelism is expressed by syntactic constructions, async/finish in X10 or cobegin/coend in Chapel. For some algorithms, it is necessary to restrict temporarily the degree of parallelism, using synchronization objects, called clocks in X10 or phasers in Habanero Java. These primitives are somewhat redundant, and may be used interchangeably in some circumstances. The aim of this paper is to explore these redundancies for X10, and to evaluate their impact on program performance.

Our key contributions are:

we give a general scheme for clock elimination, which applies only to static control programs, -we show that this scheme is correct and does not lose parallelism, -for polyhedral programs, the control overhead of the target program can be reduced or even eliminated by loop transformations, -experiments show that for the latest version of the X10 compiler and runtime, the proposed transformation improves the running time for fine grain parallel programs.

The rest of the paper is structured as follows. We will first give as much information on X10 as necessary to understand our approach. We will then define the polyhedral subset of X10. While our approach is not limited to this subset, it gives the best results in the case of polyhedral programs. results.

The X10 Language

The Base Language X10 is an object oriented language of the Java family. It has classes and methods, assignments and method invocation, and the usual control constructs: conditionals and loops. Dealing with method invocation necessitates interprocedural analysis, and is beyond the scope of this paper. The exact shape of assignments is irrelevant in this work.

X10 has two kind of loops: the ordinary Java loop:

for(<initialization>; <tests>; <increment>) S and an enumerator loop:

for(x in <range>) S
where the type of the counter x is inferred from the type of the range.

Concurrency Concurrency is expressed in X10 by two constructs, async S and finish S, where S is an arbitrary statement or statement block. Such constructs can be arbitrarily nested, except that the whole program is always embedded in an implicit or explicit finish, which creates the main activity. The effect of async S is to create a new activity or lightweight thread, which executes S in parallel with the rest of the program. The effect of finish S is to launch the execution of S, then to wait until all activities which were created inside S have terminated. X10 also allows the distribution of work on a set of logical places (typically, various nodes of a compute cluster), in a way that is transparent to the organization of activities. This aspect is orthogonal to the work explained in this paper, and will not be further evoked.

Synchronization

In some cases, it may be necessary to synchronize several parallel activities. This can be achieved using clocks. Clocks are created by clocked finish constructs. Activities are registered to the clock associated to the innermost enclosing clocked finish if created by a clocked async construct. An activity deregisters itself from its associated clock when it terminates. Synchronization occurs when an activity executes the advance primitive. This activity is held until all registered activities have executed an advance, at which time all registered activities are released.

Clocks can be seen as generalization of the classical barriers. The main differences are that activities may be distributed among several clocks which work independently, and that this distribution is dynamic as it can change when an activity is created or terminated. Refer to Figure 1 for a sample X10 program.

Intuitively, it should be clear that an unclocked finish or a set of advances can be used interchangeably. In both cases, several activities are stalled until all of them have reached a synchronization point. If a clock is used, then all clocked activities are released for further processing, while in the case of a finish, they are destroyed. The aim of this paper is to explore this analogy, both from the point of view of expressiveness and from the point of view of performance.

The Polyhedral Subset of X10

In general, analysis of arbitrary programs in a high level language like X10 is difficult, due to the presence of dynamic constructs like while loops, tests, and method invocation. Hence, many authors [START_REF] Feautrier | The polyhedral model[END_REF] have defined the so-called polyhedral model, in which many analyzes can be done at compile time. The polyhedral subset of X10 has been defined in [START_REF] Yuki | Array dataflow analysis for polyhedral X10 programs[END_REF]. The present section is a summary of this work. An X10 program is in the polyhedral model if its data structures are arrays and its control structures are loops. An enumerator loop is polyhedral if the range is an integer interval. The bounds of the range and the array subscripts must be affine functions of surrounding loops counters and integer parameters.

If these conditions are met, one can define statement instances, iteration domains, and an order of execution or happens-before relation. Statement in-stances are named by position vectors, which are deduced from abstract syntax trees (AST).

Consider the following example and its AST: The position vector for an elementary statement S is obtained by following the unique path in the AST from the root to S. In the example, the position vector of S 1 is [f, 0, i, 0] and that of S 2 is [f, 0, i, 1, a, 0], where f stands for finish and a stands for async. Let x and y be two position vectors, and let us write x ≺ y for "instance x happens before instance y". To decide whether x ≺ y, first expand x ≪ y, where ≪ is the ordinary lexicographic order. Then, remove a term if, after elimination of a common prefix, the first letter one encounter on the left is an a. This rule reflects the fact that the only temporal effect of async S is to postpone the execution of S. The reader may care to check that in the above example, instances of S 2 are unordered, while S 1 (i) happens before

S 2 (i ′ ) if i < i ′ .
Another construction is necessary for programs that use clocks. The simplest case is that of one-clock programs (or of innermost clocked finishes). One must distinguish the unclocked happens-before relation, for which advances are treated as ordinary statements, and the clocked happens-before, noted ≺ ≺. Let A be the set of advances inside one clocked finish. The advance counter at operation u is defined as:

φ(u) = Card{u ′ ∈ A | u ′ ≺ u}.
When the effect of clocks is taken into account, one can prove that if φ(u) < φ(v), then u happens before v. As a consequence, the clocked happens-before relation is:

u ≺ ≺ v ≡ φ(u) < φ(v) ∨ u ≺ v.
Since for polyhedral programs A is a union of disjoint polyhedra, and u ′ ≺ u is a disjunction of affine inequalities, the set {u ′ ∈ A | u ′ ≺ u} is the set of integer points which belong to a union of polyhedra. The cardinal of this set can be computed in closed form using the theory of Ehrhart polynomial, for which there exists efficient libraries [START_REF] Verdoolaege | Counting integer points in parametric polytopes using Barvinok's rational functions[END_REF].

A Generic Transformation Strategy

Our goal is to remove clocks from X10 programs. To understand the idea of this transformation, consider Figure 2: the center graph depicts the execution of an imaginary X10 program, where activities are represented by vertical boxes that contain regular instruction executions and clock synchronization operations. These activities "align" on their calls to advance(). The code on the left side of the figure is one possible source of this program. The idea of the transformation is to extract "slices" (or phases) across activities, represented by horizontal dashed boxes on the graph. A possible corresponding program appears on the right of the figure: the usage of clocks has been replaced by the barrier ending finish blocks. We will prove in the next section that both programs execute the same operations in the same order, except for clocks and the number (and duration) of activities. This transformation, can be implemented by a straightforward technique. Starting with a given clocked finish block, the result of the transformation can be sketched as follows:

for (d ...) // where d is a monotonically increasing phase number finish S d // the original finish block restricted to phase d

Writing the transformed program this way assumes that it is possible 1) to determine the number of phases of the program, either statically or dynamically, 2) to execute the given block for a given phase only (the restriction of the block to that phase), and 3) to repeatedly execute the original block. The rest of this section explores these three issues.

Motivating Example

Our goal in this section is three-fold. First, it is important to understand what class of programs the transformation can be applied to. Second, the example will help pinpointing potential optimizations. And third, we want to empirically validate our intuition that managing clocks is more expensive than creating activities.

Our working example appears on the left of Figure 3. The finish block creates only two activities in addition to the main activity. Each of these execute a loop that does some work (in abstract instructions S0 and S1), and then conditionally synchronizes with the other. A set of input parameters, contained in arrays a0 and a1, drives the control of the program and the synchronization scheme. These parameters make it impossible to statically derive how many phases the program has, and how many executions of S0 and S1 are performed in each phase. The resulting program appears on the right of Figure 3. The transformation can be broken into four successive steps:

clocked finish { clocked async { for (i in 0..(N-1)) { S0(i); if (a0(i) > 0) advance(); } } clocked async { for (i in 0..(N-1)) { S1(i); if (a1(i) > 0) advance(); } } } 1 cont = true; 2 for (d=0 ; cont ; d++) { 3 cont = false; 4 finish {φ := 0; 5 async { φ0 := φ; 6 for (i in 0..(N-1)) { 7 if (d == φ0) S0(i); 8 if (a0(i) > 0) 9 ++ φ0;
1. The finish block is wrapped inside a loop over d, whose iterations represent the various phases of the execution (line 2 on Fig. 3). The exit condition is represented with a boolean, named cont, whose role is detailed in the fourth phase.

2. Every activity gets its own local "counter" (named φ, φ 0 and φ 1 in the example), 4 To evaluate the performance impact of the transformation, we still need to give some definition to S0(i) and S1(i). In the experiment below, we use some "dummy" code of the form:

for (t in 1..T) a(i) += garbage(k%4)
that is to say, two accesses to arrays plus two arithmetic operations (subject to optimization). The T parameter is used to control the amount of work performed by one call to either S0 and S1: on a recent laptop, we have observed that such a loop takes roughly T nanoseconds. To run either the original or the modified program, the arrays a0 and a1 are filled with randomly generated values with equiprobable signs. Figure 4 shows the execution times in milliseconds of both versions with N = 100 as a function of the parameter T. The original version uses clocks to synchronize both activities, whereas the modified version simply repeats the whole finish block as many times as necessary (therefore creating many more activities). These curves are surprisingly close to each other. For moderately heavy instruction grain (here between 10 and 100 µs per call to S0 or S1), it seems that the cost of handling clocks is approximately as high as executing around 50 instances of the block (including the creation of activities). This accomplishes our third goal, and validates our intuition that clocks are expensive.

Applicability and Correctness

There are two main aspects in the generic transformation:

1. guarding the instructions, so as to have them execute during the right phase; 2. maintaining phase numbers ("dates") during each iteration of the loop on d. Let us for a moment assume that the second aspect is enforced. In that case, it is easy to see that both versions of the program are equivalent. For, if two instructions of the original block are executed during different phases, then they will be executed inside different finish blocks in the transformed program. And if they are executed during the same phase, they will still be executed in the right order, since the transformed block is a copy of the original block, and thus faithfully reproduces program order. Therefore, correctness of the transformation is guaranteed if phase numbers can be correctly maintained.

Maintaining correct dates at all times, i.e., during every iteration of the dloop, however, is not possible for all programs. Here is a simple modification of our previous example, where the body of the i-loop inside the first async becomes:

if (a0(i) > 0) { a0(i) = -1; advance(); }
Here, code executed at date d updates a value that will be used at a later date (an iteration d ′ with d ′ > d). This means that later iterations of the d-loop will not be able to maintain phase numbers correctly, leading to an incorrect result.

The general criterion to distinguish programs that can be transformed correctly is the following: every iteration of the d-loop must perform exactly the same sequence of advance counter incrementations. To formally capture this notion, let us define control variables: a variable that is used in a conditional branch (including loop back-edges), or to update another control variable (with arrays considered as single variables). 6 Then, a program will be correctly transformed, if the history of each control variable is the same in each iteration of the d loop, i.e. if no control variable is live in at the beginning of the d loop. For sequential structured programs, this can be checked by many classical algorithms, including reaching definition analysis and transformation to SSA. These algorithms can be extended to parallel programs: see for instance [START_REF] Lee | Basic compiler algorithms for parallel programs[END_REF], where a Concurrent Static Single Assignment form is defined for programs with parallel construct similar to those of X10, including post / wait synchronization. The results of this analysis are approximate. For polyhedral X10 programs, it is possible to do better, as shown in [START_REF] Yuki | Array dataflow analysis for polyhedral X10 programs[END_REF].

Optimization Opportunities

The generic transformation described above uses a very costly strategy: reexecuting the original code again and again, inhibiting the execution of almost all instructions at each iteration. Even though this cost seems to be amortized even for moderately heavy computations, the whole structure of the transformed program is unsatisfactory. This section tries to highlight characteristics of the transformed program that may lead to further simplification. The various steps of the intuitive transformation provide important clues on classes of programs where this applies.

The first aspect is about the local counters that have to be maintained to model the "date" inside an activity. In some cases, the date may be available at compile time as a closed form function of loop counters. Or it can be precomputed to avoid repeated incrementations of the local counter. In our example, precomputation would fill two arrays d0 and d1, indexed on i and containing the date at iteration i. The code of the first activity becomes:

for (i in 0..(N-1)) if (d == d0(i)) S0(i);
In other cases, like the ones described in the next section, dates are functions of the enclosing loop counters, and do not need dedicated storage. The second aspect is very much related to the first, and relates to the upper bound of the enclosing loop. When dates are available, it is immediate to compute or memorize their maximal value (which is the upper bound on d) . This also removes the need of a boolean variable and tests at the end of activities.

The third aspect relates to the interplay between statement guards and loop bounds. In our example, we have reached a situation where a loop iterates from 0 to N -1, but where only a sub-range of this iteration space leads to actual execution. It is therefore possible to adjust the range of the loop to cover only the relevant sub-range. In our example, the first activity becomes: while (d0(i0) == d) { S0(i0); ++ i0; } where i0 is a global counter, suitably initialized and preserved across activities.

We have given an informal account on how a program with clocks removed can be further simplified and optimized. The next section describes a class of programs where these optimizations can be systematically applied, and details their implementation.

General Polyhedral Programs

Polyhedral programs with clocks have the property that a date can be computed directly for any instruction, by counting the number of calls to advance() performed before an instance of the given instruction. This removes the need to maintain explicit counters, and provides symbolic expressions involving loop counters and symbolic parameters. One can always evaluate the number of integer poins inside a (parametrized) polyhedron, and therefore assign a monotonically increasing rank to any instance of an instruction. Under reasonable assumptions, i.e., that loops have unit steps, such ranks are integer-valued polynomials with rational coefficients [START_REF] Clauss | Counting solutions to linear and nonlinear constraints through Ehrhart polynomials: applications to analyze and transform scientific programs[END_REF].

An example program appears on Figure 5, with date expressions placed in comments. Note that the counting happens in two phases: first, the starting date of an activity is computed, and second the date of instructions are computed relative to the activity's starting date. The transformation process starts by computing the maximal date at which an instruction of the original finish block executes. This can be done by maximizing for each instruction individually, and then taking the maximum of the results. This maximum, N -2 + M (M -1)/2 in our example, is the upper bound of the loop wrapped around the original block. Then, calls to advance() are removed, and guards are placed around statements. The result appears on the right part of Figure 5.

After having inserted guards around statements, the last step is to examine whether the guards have an impact on the bounds of the loops that enclose the statement. Our example illustrates this situation: after transformation, the innermost loop becomes:

for (k in 0..(j-1)) if (d == i+j*(j-1)/2+k) S1(i,j,k);
A trivial rewriting of the guard shows that even though the loop iterates over a range of values for k, the whole loop will actually execute S1(i,j,k) at most once. This construct can therefore be replaced by:

k = d -i -j*(j-1)/2; if (0<=k && k<=j-1) S1(i,j,k);
which tests whether the single value selected by the guard is inside the range of the loop. Note that we started with a depth three loop nest. Then a new loop level was added around this nest. And finally the deepest level is removed. This is likely to reduce the overhead introduced by the transformation. This optimization may be extended to loops containing several statements (at the same or different dates). However, it applies only when date expressions are linear in the nearest enclosing loop counter, which we think is a very common case. Actually, for this not to be the case, an innermost statement-bearing loop should also contain another loop 1) containing only calls to advance(), and 2) with a bound being a function of its parent loop counter. Here is the simplest example of such a construction: for (y in ...) { S(...); for (z in 0..y) advance(); } We think this pathological case and its variations are sufficiently infrequent not to cause real trouble in practice.

Note that after a loop is removed, the statement is still guarded, but with a condition involving inequalities. Therefore, there is no possibility of re-applying the same "iteration space collapsing", but nothing says that the new guard may not imply bound adjustments on the enclosing loops. The next section shows an example of such chained loop adjustments.

Polyhedral Programs with Affine Dates

We have seen that when polynomial dates are available, the resulting program can be optimized by combining guards and loop bounds. However, dealing with polynomials of high degree is difficult and may restricts how far optimizations can go. It is therefore interesting to consider the particular case of affine dates. In that case, all obstacles to optimization are lifted, and one can hope to be able to optimize the transformed program up to the point where it has the same complexity as the original program. 7Whenever the original program induces dates that are all affine forms in the enclosing loop counters (and parameters), we are guaranteed that the deepest loop level can be removed. In fact, this last level of loop contains only guarded statements, and the guards are of the form d = α 0 i 0 +. . .+α n i n , where i 0 , . . . , i n are the counters of the enclosing loops. Such a guard always determines at most one value of the counter of the nearest loop. This property appears in the program in Figure 6. The left part shows the original clocked finish block (dates appear in comments), whereas the right part shows the mechanically transformed program. Since dates are affine, one can immediately apply the "iteration space collapsing" optimization mentioned in the previous section. The first loop on j then becomes:

if (i<=d && d<=N-1) S0(i,d);
and the second loop on j can be transformed as well (note that we do not keep a variable to store the value of j, but rather substitute it immediately). The major advantage of having affine dates is the fact that the resulting program can be further optimized. We are going to illustrate these additional optimizations on the example program in Figure 6, and then we will show how the program transformations involved are strongly related to the problem of code generation from a polyhedral model of a program. We will then show, in the next section, how existing tools can be adapted to directly produce the optimized version.

Regarding the example of Figure 6, the first step is to replace constructs of the form async if (...) S(...) with if (...) async S(...), because there is no need to create an activity that does nothing. All these initial modifications lead to the following program:

for (d in 0..(2*N-2)) finish for (i in 0..(N-1)) { if (i<=d && d<=N-1) async S0(i,d);
S0 := [N]
->{[f,0,i,0,a,0,j,0,d]: 0<=i<N and i<=j<N and d=j}; 2. Representing the whole program simply consists in computing the union of the individual instruction polyhedra: P := S0+S1; Here P represents the set of instances of S0 and S1, in the order of the original program. 83. Iterating on dates first is performed by changing the order of the dimensions. This is written as: U := {[f,p0,i,p1,a,p2,j,p3,d]->[d,f,p0,i,p1,a,p2,j,p3]}(P); This actually doesn't do anything, but is an important indication to the next step. 4. Producing the final code is performed by generating a program scanning the resulting polyhedron U. We use the CLooG algorithm [START_REF] Bastoul | Code generation in the polyhedral model is easier than you think[END_REF], which produces a new loop nest with a loop scanning dates (d) first, and whose body contains various constructions (finish, async, loops, and instructions) in the order prescribed by the various other dimensions (the original CLooG algorithm had to be hacked to handle finish and async).

The final code (after trivial cosmetic post-processing) appears on Figure 9: CLooG has adjusted all loop bounds (even though it could have gone further). This code generation phase actually under-uses CLooG, which is able to apply a "scattering function" (taking its default value, the identity, in our case). The same result could be obtained by applying, e.g., some variation of Fourier-Motzkin elimination for bound adjustment [START_REF] Ancourt | Scanning polyhedra with DO loops[END_REF]. However, reconstructing the structure would still need additional work. CLooG does both iteration domain computation and code generation.

Experimental Results

To evaluate the effect of eliminating clocks on execution time, we have used eight different polyhedral programs with affine dates. All these programs are parametrized by a number N that determines the number of activities and the number of iterations of loops in various ways. Their execution is depicted on Figure 10 for N = 6: vertical lines represent activities, dots represent individual instruction executions, and horizontal dashed lines represent phases of execution. To compare the clocked and clockless versions of each program, we have measured their execution times (averaged over 20 executions). We have used a not-quite-recent X86-64 compatible AMD machine with 24 cores (two sockets of 12 cores). X10 programs were compiled with the official release of X10, version 2.3.1, available from http://x10-lang.org/. All programs have been run with N = 100. Because the elimination of clocks affects only the control of the program, and not its actual work, we have varied the time taken by a single instruction execution (a call of the form S k (i, j) in all cases) the same way we did in Section 2.1: a single parameter T controls how much time a single execution of any S k (i, j) takes. The goal of the experiment is therefore to measure the difference in execution times as a function of T .

The results are shown on Figure 11. Every graph shows the execution time of both versions. In all cases, the clockless version runs faster than the version with clocks. Rows of three graphs show the times of a given program for various values of T (the workload): every graph displays ten evenly spaced values of T , with one order of magnitude variation from one graph to the next. The vertical scales are different from one graph to the other, but all scales are zero-based.

Since X10 is not the only language allowing finish/async programming, we have also conducted preliminary experiments with Habanero-Java [START_REF] Cavé | Habanero-java: the new adventures of old X10[END_REF] (version 1.3.1), with results similar to those presented here.

There are several lessons to learn from Figure 11. First, eliminating clocks always has a positive impact on execution time. This validates our intuition that clocks are expensive to manage. At least their use is more expensive than launching more activities (by a factor O(N ) in our case). We acknowledge that this is fairly dependent on implementation issues, but we also think that it will be easier to optimize activity creation rather than clock synchronization. Future implementations of X10 (and related languages, like Habanero and Chapel [START_REF] Chamberlain | Parallel programmability and the Chapel language[END_REF]) may change this situation.

Examination of the leftmost column of Figure 11 shows the relatively irregular behavior of programs using clocks with fine-grain instructions: it looks as if the frequent calls to advance() make the actual time difficult to predict, whereas clockless programs display a smoother, quasi-linear curve. Again, this heavily depends on the implementation of the activity scheduler, but it seems clear that clockless programs are "easier" to schedule over an arbitrary number of threads.

The third lesson learned is that, as expected, the difference between versions vanishes when the workload is reasonably large, because the time spent in control becomes negligible compared to the time spent on computation. What is less obvious is that this happens for values of T around one million (which, on our machine, is about 1 millisecond). Considering the kind of programs we have used (basically loops over arrays, where every instruction accesses one element of one or more arrays), there is little chance that this workload is reached. This means that for fine-grain programs, the transformation is probably advantageous, providing significant speedup in most cases.

Related Work

There exists a large body of literature on barriers and clocks, their analysis, optimization and verification. Nearest to the subject of this paper is work on T from 100K to 1M Fig. 11. Execution times for clocked ( ) and clockless ( ) versions, for various scales of workload. Vertical scales vary from graph to graph, but are all zero-based.

optimal barrier placement [START_REF] Aiken | Barrier inference[END_REF][START_REF] Kamil | Concurrency analysis for parallel programs with textually aligned barriers[END_REF][START_REF] Darte | A linear-time algorithm for optimal barrier placement[END_REF] and verification [START_REF] Vasudevan | Compile-time analysis and specialization of clocks in concurrent programs[END_REF]. While apparently related to the present work, Chau-Weng Tseng paper [START_REF] Tseng | Compiler optimizations for eliminating barrier synchronization[END_REF] deals in fact with a completely different problem, namely how to distribute work among threads in order to minimize synchronization. Several authors have argued that barriers or clocks can be implemented more efficiently than task or activity creation, and have advocated algorithms for minimizing the number of tasks. To the best of our knowledge, the word SPMDization was coined by Padua and Paek in [START_REF] Padua | Compiling for scalable multiprocessors with Polaris[END_REF]. A recent discussion of the same idea is by Zaho et. al. [START_REF] Zhao | Reducing task creation and termination overhead in explicitly parallel programs[END_REF] in which an algorithm, which amount to moving parallel loops outside sequential loops with barrier insertion is proposed. Our contention here is that moving in the opposite direction may be beneficial in some cases. Choosing between the two solutions depends on many factors: the target system, the compiler and runtime and the source program. For instance, if the target is hardware, where it is almost impossible to dynamically create activities and synchronization is cheap, using clocks might be the best solution. Our work shows that the situation is exactly the reverse for software.

Conclusion

When one has to generate a parallel program, either manually or automatically, one has to choose between two extreme program shapes: one a sequence of parallel constructs, the other a parallel composition of sequential threads. Obviously, these two extreme cases can be combined to produce many intermediate solutions.

In the first approach, it is usually possible to restrict synchronization to one barrier after each parallel block. This corresponds to the exclusive use of async / finish in X10, and is especially suitable for vector or data-parallelism. In the second approach, it is usually not possible, except in the case of embarassingly parallel programs, to construct independent threads. Residual dependences must be satisfied using clock or phasers [START_REF] Cavé | Habanero-java: the new adventures of old X10[END_REF]. This work shows that deciding which approach gives the best performance is not obvious, and must be approached experimentally. Our main contribution is a systematic method for converting a large class of clocked programs into unclocked ones. Our algorithms can easily be automated, thus simplifying the comparison process.

This paper has introduced a program transformation that acts on an explicitly parallel program, an unusual characteristic in the polyhedral framework. Such an ability opens up a large space of new potential optimizations, extending the scope of automatic parallelization. Taking the cost of synchronization primitives into account must also be extended and further generalized, to cover cases where implementations may have different semantics and/or relative overheads. Also, the cost of synchronization is only one part of the picture, and more work is needed to combine synchronization costs with more "traditional" transformation objectives in the polyhedral framework, like, e.g., temporal and spatial locality. Finally, we plan to investigate the use of parallel-to-parallel program transformations in dynamic optimization frameworks, where switching between various versions of the same program can alleviate the variation of synchronization costs linked to changing runtime conditions.

Fig. 1 .

 1 Fig. 1. A Sample Program

Fig. 2 .

 2 Fig. 2. Parallelism and synchronization in X10, with and without clocks

10 } 11 if|Fig. 3 .

 10113 Fig. 3. An example program, before and after transformation

Fig. 4 .

 4 Fig. 4. Clocked and clockless execution times (in milliseconds)

Fig. 5 .

 5 Fig. 5. A polyhedral program with polynomial dates (in comments) on the left, and the result of the transformation, on the right.

Fig. 8 .

 8 Fig. 8. An AST for the program on Figure 6, and the corresponding polyhedra

Fig. 9 .

 9 Fig. 9. The final result, produced by CLooG

Fig. 10 .

 10 Fig. 10. Example iteration spaces, here for N = 6. All examples spawn O(N ) activities, last for O(N ) clock steps, and execute O(N 2 ) instructions.

  initialized at the start of the activity by capturing the value of the parent activity's counter if any (lines 4, 5, and 12). Local counters are maintained by replacing calls to advance() by an incrementation (lines 9 and 16). 3. All instructions that have an effect visible outside the finish block are guarded (lines 7 and 14), and the guard condition checks whether the value of the local counter matches the currently executed phase (given by d). 4. Finally, when any activity reaches its end, the value of the local counter has reached its maximum value for that activity. This maximum value is the index of the last phase for which this activity has work to do. A simple test decides whether the loop on d should continue iterating (lines 11, 18, and 19).5

x10.sourceforge.net/documentation/languagespec/x10-latest.pdf

The local counter of the activity executing the body of finish, named φ, is useless here and was left for completeness only.

Activities could be aborted once their local counter is above the value of d: this aspect is more or less orthogonal to our goal, and is ignored here.

In practice, the collection of control-variables can be restricted to programming constructs containing at least one call to advance().

Note that the complexity in terms of the number of executions of individual instructions is always the same on both versions. Here we refer to the complexity of the associated control, i.e., the number of times guards and loop exit conditions are evaluated.

For this union operation to have any meaning, the dimension lists of the various must coincide; this is trivially achieved by padding with zeros. No modification is necessary in our example.

if (i+1<=d && d<=2*i) async S1(i,d-i-1); } At this point, all remaining optimizations are made possible by the comparison of the various inequalities that apply to the individual instructions. Since our goal is to reduce the time taken by evaluating the guards, we are going to rearrange this code to eliminate useless guards and uselessly large bounds.

The first batch of useless evaluations of guards is caused by the d<=N-1 condition, because at this point d is supposed to iterate from zero to 2*N-2. This means that half of the values of d will simply fail to satisfy the condition. Eliminating these useless tests requires that the range of d is split into two sub-ranges, the first of which makes the condition trivially true, and the second which makes it false. The result appears in Figure 2.5. Range-splitting globally enlarges the code, but removes any occurrence of S0 in the loop iterating over the second sub-range. Note also that condition i+1<=d around S1 has become trivially true in the second loop, and is therefore also omitted.

The second set of unnecessary tests is caused by the remaining conditions, which in all cases are stricter than the surrounding loop bounds. The range of the first loop on i can be split into three sub-ranges, namely 0..(d-1), d, and (d+1)..(N-1): the first leads to the bulk of the work, the second selects only S0, and the third leads to nothing. The result of bound adjustment appears in Figure 2.5. What was just done on upper bounds can now be done on lower bounds as well: the condition d<=2*i appears twice inside loops whose lower bound on i is zero, for any value of d. Therefore, the lower bound can be adjusted as well. The details are left to the reader.

Polyhedral Implementation and Optimized Control

The approach we have taken in the previous section consists in a succession of elementary transformations: wrapping a loop around the original code, placing guards around elementary statements, and adjusting iteration domains according to the guards. In contrast, in the polyhedral model, all these transformations can be represented in a uniform framework, and polyhedral operations can be used to manipulate the program. A polyhedral model of an instruction (a polyhedron, for short) is made up from two distinct parts: first, an ordered list of dimensions, and second a set of constraints (inequalities) on the values of the various dimensions. There are three types of dimensions: 1) syntactic dimensions, which are usually constants, 2) loop iterators, and 3) parallel constructs indicators, which are the abstract symbols f(inish) and a(sync). Figure 8 displays the polyhedra corresponding to the instructions appearing in the original program of Figure 6. The left part of the figure shows an abstract syntax tree, which is convenient to read the various dimensions. The right part shows the polyhedra, using the notation of the iscc polyhedral calculator, part of the barvinok library [START_REF] Verdoolaege | Counting integer points in parametric polytopes using Barvinok's rational functions[END_REF]. Note that polyhedra can be parametrized (by N in our case), and that constant dimensions can be written literally, i.e., {[f,0,...]: ...} is equivalent to {[f,p_0,...]: p_0=0 and ...}.

All manipulations necessary for the elimination of clocks can now be formulated as operations on polyhedra:

1. Introducing dates is performed by adding a dimension, at the very end of the list of dimensions since the date may depend on any of the enclosing loop counters. For instance, the definition of S0 becomes: