Frederic Mazenc
email: frederic.mazenc@l2s.centralesupelec.fr

Vincent Andrieu
email: vincent.andrieu@gmail.com

Michael Malisoff
email: mal-isoff@lsu.edu

‡ Andrieu

Continuous-Discrete Observers for Time-Varying Nonlinear Systems: A Tutorial on Recent Results *

Continuous-discrete systems can occur when the plant state evolves in continuous time but the output values are only available at discrete instants. Continuous-discrete observers have the valuable property that the observation error between the true state of the system and the observer state converges to zero in a uniform way. The design of continuousdiscrete observers can often be done by building framers, which provide componentwise upper and lower bounds for the plant state. This paper is a tutorial on these approaches, highlighting recent results in the literature, and also providing previously unpublished, original results which are not being simultaneously submitted elsewhere.

Introduction

The search for more effective designs for observers for nonlinear systems has led to a substantial and complex literature [START_REF] Andrieu | Extended-Kalman-filter-like observers for continuous time systems with discrete time measurements[END_REF][START_REF] Andrieu | Observer design for Lipschitz systems with discrete-time measurements[END_REF][START_REF] Arcak | A framework for nonlinear sampled-data observer design via approximate discrete-time models and emulation[END_REF][START_REF] Dinh | Continuous-discrete time observer design for Lipschitz systems with sampled measurements[END_REF][START_REF] Karafyllis | From continuous-time design to sampled-data design of observers[END_REF][START_REF] Zemouche | On LMI conditions to design observers for Lipschitz nonlinear systems[END_REF][START_REF] Zemouche | Observers for a class of Lipschitz systems with extension to H ∞ performance analysis[END_REF]. The work is motivated by a plethora of real world applications where it may be difficult or impossible to measure the state variables. Then the goal is to use output measurements to design an observer for the state such that the observation error between the observer and the state converges to 0 as time goes to infinity. The well known literature on observers is largely for systems with continuous measurements; see, e.g., [START_REF] Zemouche | Observers for a class of Lipschitz systems with extension to H ∞ performance analysis[END_REF] for results based on writing the differential equation satisfied by the estimation error as a linear parameter varying system.

However, in real world applications, the output is often unavailable for continuous measurement. Since the dynamics are in continuous time and the output is only available for measurement at discrete instants, this produces a continuous-discrete system. There is now a large literature on observer designs for continuousdiscrete systems that spans over forty years. See, e.g., the work [START_REF] Jazwinski | Stochastic Processes and Filtering Theory[END_REF] of Jazwinski, which used a continuousdiscrete Kalman filter to solve a filtering problem for stochastic continuous-discrete time systems. The high gain observer approach in [START_REF] Gauthier | A simple observer for nonlinear systems applications to bioreactors[END_REF] was extended to continuous-discrete systems in [START_REF] Deza | High gain estimation for nonlinear systems[END_REF], where the impulsive correction gain is found using a continuousdiscrete Riccati equation. The robustness of observers under discretization was studied in [START_REF] Arcak | A framework for nonlinear sampled-data observer design via approximate discrete-time models and emulation[END_REF], and [START_REF] Ahmed-Ali | Continuous-discrete observer for state affine systems with sampled and delayed measurements[END_REF][START_REF] Farza | Continuous-discrete time observers for a class of MIMO nonlinear systems[END_REF][START_REF] Karafyllis | From continuous-time design to sampled-data design of observers[END_REF] used output predictors to design observers; see also the works [START_REF] Andrieu | Observer design for Lipschitz systems with discrete-time measurements[END_REF][START_REF] Deza | High gain estimation for nonlinear systems[END_REF][START_REF] Hammouri | Constant gain observer for continuous-discrete time uniformly observable systems[END_REF][START_REF] Mazenc | Construction of interval observers for continuous-time systems with discrete measurements[END_REF][START_REF] Tellez-Anguiano | Nonlinear continuous-discrete observer application to distillation columns[END_REF]. The paper [START_REF] Ahmed-Ali | Continuous-discrete adaptive observers for state affine systems[END_REF] designed continuous-discrete observers for nonlinear continuous time systems, where the input of the system satisfies a persistent excitation condition, and [START_REF] Nadri | Design of a continuous-discrete observer for state affine systems[END_REF] covers systems that are linear in the state and have known inputs. For an interesting application, see [START_REF] Astorga | Nonlinear continuousdiscrete observers: Applications to emulsion polymerization reactors[END_REF] for continuous-discrete observers for emulsion polymerization reactors.

In the recent work [START_REF] Mazenc | Design of continuous-discrete observers for time-varying nonlinear systems[END_REF], we built on [START_REF] Andrieu | Observer design for Lipschitz systems with discrete-time measurements[END_REF], by finding a new class of continuous-discrete observers for continuous time Lipschitz systems with discrete measurements. As in [START_REF] Andrieu | Observer design for Lipschitz systems with discrete-time measurements[END_REF] and [START_REF] Deza | High gain estimation for nonlinear systems[END_REF], the continuous-discrete observers in [START_REF] Mazenc | Design of continuous-discrete observers for time-varying nonlinear systems[END_REF] are obtained in two steps. First, when the output is not available for measurement, the state estimate is computed by integrating the model. Then, when a measurement occurs, the observer makes an impulsive correction to the estimated state.

The works [START_REF] Andrieu | Observer design for Lipschitz systems with discrete-time measurements[END_REF] and [START_REF] Dinh | Continuous-discrete time observer design for Lipschitz systems with sampled measurements[END_REF] used this two step approach to show that when no measurement occurs, the estimation error is a solution of an appropriate unknown linear parameter varying system. This led to a construction of a framer, meaning, an upper and lower bound for the solution vector, in a vector sense, which made it possible to design correction terms that ensure that the estimation error asymptotically converges to zero. However, [START_REF] Andrieu | Observer design for Lipschitz systems with discrete-time measurements[END_REF] and [START_REF] Dinh | Continuous-discrete time observer design for Lipschitz systems with sampled measurements[END_REF] find the framer by integrating a system with commutation, which does not lead to an explicit analytic expression for the framer, and [START_REF] Mazenc | Construction of interval observers for continuous-time systems with discrete measurements[END_REF] is limited to linear systems. By contrast, [START_REF] Mazenc | Design of continuous-discrete observers for time-varying nonlinear systems[END_REF] used an approach from [START_REF] Cacace | A new approach to the internally positive representation of linear MIMO systems[END_REF] on cooperative systems [START_REF] Haddad | Nonnegative and Compartmental Dynamical Systems[END_REF] to get analytic constructions of framers, which can be useful for applications where explicit expressions are needed.

This paper provides a tutorial that explains some of the preceding advances precisely, including the motivation for the assumptions and methods and the value added by our recent contributions [START_REF] Andrieu | Observer design for Lipschitz systems with discrete-time measurements[END_REF][START_REF] Dinh | Continuous-discrete time observer design for Lipschitz systems with sampled measurements[END_REF][START_REF] Mazenc | Design of continuous-discrete observers for time-varying nonlinear systems[END_REF], while also stating and proving some previously unpublished, original results on framers that are not being simultaneously submitted elsewhere. We believe that our tutorial will fill an important void in the literature, and increase the control community's appreciation for, and understanding of, continuous-discrete observers. In the next section, we provide the relevant definitions. In Section 3, we discuss the work [START_REF] Zemouche | Observers for a class of Lipschitz systems with extension to H ∞ performance analysis[END_REF] of Zemouche and others in the continuous time case, which illustrates one of the recurring themes in this article, namely, the possibility of using linear matrix inequalities (or LMIs), and therefore also LMI solvers, to design asymptotic observers for nonlinear systems. Then in Section 4, we discuss an extension of [START_REF] Zemouche | Observers for a class of Lipschitz systems with extension to H ∞ performance analysis[END_REF] to the discrete time measurement case, based on computing a reachable set for controlled systems and solving LMIs.

In Section 5, we discuss our alternative approach to continuous-discrete observers, which is based on designing framers, including results that have not appeared before. In Section 6, we use our framers to formulate our latest theorem on observers. In Section 7, we show how our closed form expressions for the framers allow us to check the assumptions of our theorem using linear matrix inequalities. In Section 8, we summarize our work and suggest future research topics. For novel applications of some of the theory in this paper to pendulum and robotic DC motor dynamics, see [START_REF] Mazenc | Design of continuous-discrete observers for time-varying nonlinear systems[END_REF].

Notation, Definitions, and Basic Result

Throughout the sequel, we omit arguments of functions, when they are clear from the context, and the dimensions of the matrices are arbitrary. We set N = {1, 2, . . .}. The k × n matrix all of whose entries are 0 will also be denoted by 0, and we use A = [a ij ] to indicate that an arbitrary matrix A ∈ R k×n has a ij in its ith row and jth column for each i ∈ {1, 2, . . . , k} and j ∈ {1, 2, . . . , n}. Also, I s is the identity matrix in any dimension s. The usual Euclidean norm

x 2 1 + . . . + x 2 n of vectors and the induced norm of matrices are denoted by | • |. All inequalities and maxima are componentwise, i.e., if A = [a ij ] and B = [b ij ] are matrices of the same dimensions, then we use A ≤ B to mean that a ij ≤ b ij for all i and j, and max{A, B} is the matrix C = [c ij ] where c ij = max{a ij , b ij } for all i and j. A square matrix is cooperative or Metzler provided all of its off-diagonal entries are non-negative. We use to denote transpose. For each r ∈ N and each function F : [0, ∞) → R r , we use the left limits F(t -) = lim s→t,s<t F(s). A function ϕ * : R × R n → R n is uniformly Lipschitz in its second argument provided there is a constant L > 0 such that |ϕ * (t, x) -ϕ * (t, y)| ≤ L|x -y| holds for all t ∈ R, x ∈ R n , and y ∈ R n . For any square matrices A and B in R n×n , we use A B (resp., A B) to mean that X (A-B)X ≤ 0 for all x ∈ R n (resp., X (A-B)X > 0 for all X ∈ R n \ {0}). We use Conv to denote the closed convex hull.

Background on Continuous Time Observers

Consider the continuous-discrete system

(3.1) ẋ(t) = Ax(t) + φ(x(t), u(t))
with discrete output measurements at the known sample times {t k } ∞ k=1 given by (3.2)

y k = Cx(t k ) , t k+1 = t k + δ k
where A and C are constant matrices, and the δ k 's represent the sampling delays. The function u(t) can present an open or closed loop control, and is assumed to be continuous. We make this assumption on (3.1):

Assumption 1. The pair (A, C) ∈ R n×n × R p×n is observable.
Also, for each pair (i, j) of values in {1, 2, . . . , n}, there is a positive real number b ij such that

(3.3) ∂φi ∂xj (x, u) ≤ b ij for all (x, u) ∈ R n × R p .
By a continuous-discrete observer for (3.1), we mean a system of the form

(3.4) ẋ(t) = Ax(t) + φ(x(t), u(t)) , t ∈ [t k , t k+1 ) x(t k ) = x(t - k ) + K(y k -C x(t - k )) for all k ≥ 0, where (3.5) x(t - k ) = lim t→t k ,t<t k x(t).
Then the estimation problem is that of selecting K to ensure that

(3.6) lim t→+∞ |x(t) -x(t)| = 0
for all initial conditions. To motivate our search for K, we recall the results from [START_REF] Zemouche | Observers for a class of Lipschitz systems with extension to H ∞ performance analysis[END_REF] on the continuous time case where y(t) = Cx(t). Let R be the set of all matrices in R n×n such that for each matrix R = [r ij ] in R and each pair (i, j), the corresponding entry r ij is either b ij or -b ij . We then have the following result from [START_REF] Zemouche | Observers for a class of Lipschitz systems with extension to H ∞ performance analysis[END_REF]: Theorem 3.1. If there exist a positive definite symmetric matrix P ∈ R n×n and a matrix L ∈ R p×n such that the LMI

(3.7) (A + R) P + P (A + R) -C L -L C < 0 holds for all R ∈ R, then the system (3.8) ẋ(t) = Ax(t) + φ(x(t), u(t)) + P -1 L (y(t) -C x(t))
is an asymptotic observer for (3.1) in the continuous time case where y(t) = Cx(t), i.e., the limit (3.6) holds for all initial conditions for (3.1) and (3.8).

Although Theorem 3.1 does not cover discrete output observations, it has the important feature that the sufficient condition is stated in terms of LMIs, which makes it possible to check the assumptions using LMI solvers [START_REF] Zemouche | On LMI conditions to design observers for Lipschitz nonlinear systems[END_REF]. This contrasts with traditional LMI approaches, which are usually limited to linear time invariant systems. We next explain our approach from [START_REF] Andrieu | Observer design for Lipschitz systems with discrete-time measurements[END_REF] for extending the basic properties of the observer (3.8) to systems with discrete output measurements.

Discrete Output Measurements

To motivate our extension of [START_REF] Zemouche | Observers for a class of Lipschitz systems with extension to H ∞ performance analysis[END_REF] to the case of discrete output observations, let δ > 0 be any constant, and consider the error system      ė(t) = Ae(t) + ∆φ(x(t), u(t), e(t)) for all t ∈ [kδ, (k + 1)δ) and k ≥ 0 , e(kδ) = (I n -KC)e(kδ -) for all k ≥ 0.

where ∆φ(x, u, e) = φ(x, u) -φ(x -e, u) and e(kδ -) = lim t→kδ -e(t). Assuming that φ is C 1 in its first argument (i.e., the state) and that there are constants b ij > 0 such that (3.3) holds for all (x, u) ∈ R n × R p and all pairs (i, j), the Mean Value Theorem gives ∆φ i (x, u, e) = (∂φ i /∂x)(z i (x, e, u), u)e for all i for suitable points z i (x, e, u). Hence, between any two measurements, the error is a solution of ė(t) = Ae(t) + V (t)e(t), where

V (t) = [v ij (t)] ∈ R n×n is a matrix each of whose entries is bounded by b ij .
For each choice of e 0 ∈ R n and each real number δ > 0, let A δ (e) ⊂ R n denote the reachability set at time δ with the control constraint that (4.9)

|v ij (t)| ≤ b ij
for all i and j and t ≥ 0. This means that for each e 1 in A δ (e 0 ), there exists a function V (t) = [v ij (t)] such that (4.9) holds for all t in [0, δ] and all pairs (i, j), and such that the solution e(t) of ė(t) = Ae(t) + V (t)e(t) starting from e 0 satisfies e(δ) = e 1 . We then have the following discrete time extension from [START_REF] Andrieu | Observer design for Lipschitz systems with discrete-time measurements[END_REF]:

P M i (δ)(P -C W ) (P -W C)M i (δ) P
is positive definite for all i ∈ {1, 2, . . . , }, then the choice (4.12)

K = P -1 W
in our observer (3.8) achieves our estimation goal (3.6).

The preceding theorem can be summarized as follows. The first step is the computation of a reachable set for a controlled system. Exact computation of this set is not needed, but only the upper approximation (4.10) in terms of the M i 's for 1 ≤ i ≤ . For lower triangular systems, such an upper approximation was given in [START_REF] Andrieu | Observer design for Lipschitz systems with discrete-time measurements[END_REF]. The second step is to check the positive definiteness of the matrices (4.11). For applications of Theorem 4.1, including cases where the system is uniformly observable, see [START_REF] Andrieu | Observer design for Lipschitz systems with discrete-time measurements[END_REF]. A possible drawback of the preceding approach is that it does not lead to analytic formulas for framers. We next provide an alternative approach to continuousdiscrete observers that has the advantage of providing closed form expressions for framers.

Background on Framers

In this section, we present several results on framers for time-varying linear systems that we use in the next section to build our continuous-discrete observers for nonlinear systems. Take any linear time-varying system

(5.13) ẋ(t) = M(t)x(t)
with state space R n , where all entries of M : [0, ∞) → R n×n are continuous. Let (t, t 0 ) denote the fundamental solution of the system (5.13), meaning, (∂ /∂t)(t, t 0 ) = M(t) (t, t 0 ) and (t 0 , t 0 ) = I n hold for all t 0 ≥ 0 and t ≥ t 0 . In this section, we derive lower and upper bounds for the function Γ(t) = (t, 0). Note for later use that the solution φ of (5.14) ∂φ ∂t (t, x 0 ) = M(t)φ(t, x 0 ), φ(0, x 0 ) = x 0 satisfies φ(t, x 0 ) = Γ(t)x 0 for all t ≥ 0 and x 0 ∈ R n . Our next lemma on framers is a key ingredient needed to prove our main result on framers. It assumes: Assumption 2. There are two constant Metzler matrices M ∈ R n×n and M ∈ R n×n such that (5.15) M ≤ M(t) ≤ M for all t ≥ 0 .

Also, M : R → R n×n is continuous.

The following is shown in [START_REF] Mazenc | Design of continuous-discrete observers for time-varying nonlinear systems[END_REF]:

Lemma 5.1. If Assumption 2 holds, then exp(Mt) ≤ Γ(t) ≤ exp(Mt) hold for all t ≥ 0.

Next, we consider the system (5.13) under the following much weaker assumption than Assumption 2: Assumption 3. The matrix valued function M : [0, ∞) → R n×n is bounded and continuous.

Assumption 3 allows us to pick functions K : R → R n×n and L : R → R n×n , a constant matrix L ≥ 0, and constant Metzler matrices K and K such that (5. [START_REF] Mazenc | Design of continuous-discrete observers for time-varying nonlinear systems[END_REF])

M(t) = K(t) -L(t), 0 ≤ L(t) ≤ L,
and K ≤ K(t) ≤ K hold for all t ≥ 0. Since M is bounded, the decomposition (5.16) can be obtained by replacing each entry of

M(t) = [m ij (t)
] by m ij (t) + B for a big enough constant B > 0 to produce the Metzler matrices K(t) for each t, and then letting L = L be the constant matrix with B as each entry, but other decompositions of the type (5.16) exist. In [START_REF] Mazenc | Design of continuous-discrete observers for time-varying nonlinear systems[END_REF], the following is shown:

Lemma 5.2. Let the system (5.13) satisfy Assumption 3, and let L, K, L ∈ R n×n , K ∈ R n×n , and K ∈ R n×n satisfy the preceding requirements. Define the C 1 functions Γ : [0, ∞) → R n×n and Γ : [0, ∞) → R n×n by (5.17) Γ(t) = e Kt + 1 2 e (K-L)t -e (K+L)t and Γ(t) = 1 2 e (K+L)t + e (K-L)t .

Then Γ(t) ≤ Γ(t) ≤ Γ(t) hold for all t ≥ 0.

The proof of Lemma 5.2 is based on applying Lemma 5.1 to the 2n dimensional system (5.18) Λ(t) = H(t)Λ(t), where (5. [START_REF] Tellez-Anguiano | Nonlinear continuous-discrete observer application to distillation columns[END_REF])

H(t) = K(t) L(t) L(t) K(t) .
Two key novel features of Lemma 5.1 are that it produces many framers, for different K's and L's, and that one only needs K, L, and K to build the framers.

However, it is useful to note that one can build framers by applying Lemma 5.1 to lower dimensional systems, where the embedding of the n dimensional system is into a larger system of dimension strictly less than 2n. We next present new original results in this direction that have not been submitted elsewhere. For each continuous matrix D(t) of size n × n having the fundamental matrix q(t, t 0 ), we set φ D (t) = q(t, 0). Then φ D (0) = I n and φ M (0) = (t, 0) = Γ(t). We also set A + = max{0, A} and A -= A + -A for any matrix A, so A + ≥ 0 and A -≥ 0. We prove: Proposition 5.1. Let Assumption 2 hold. Assume that there is a constant matrix P ∈ R n×n such that

(5.20) P M(t)P -1 = B(t) = B 1 (t) B 2 (t) 0 B 3 (t)
where B 1 is valued in R s×s for some s ∈ (0, n) and B 2 and Then

B 3 are of suitable dimensions. Set K = [I s , 0] ∈ R s×n , R = P -1 , G = R + K ,
(5.22) K(L a (t)G -L b (t)F ) ≤ φ B1 (t) ≤ K(L b (t)G -L a (t)F )
hold for all t ≥ 0.

Proof: Throughout the proof, all equalities and inequalities hold for all t ≥ 0. Our condition (5.20) implies that P (∂φ M /∂t)(t)P -1 = P M(t)P -1 P φ M (t)P -1 = B(t)P φ M (t)P -1 . Hence, by the uniqueness of solutions property, φ B (t) = P φ M (t)P

P + e Mt -P -e Mt R -≤ φ B (t)P R - ≤ P + e Mt -P -e Mt R -, (5.27)
since R + ≥ 0 and R -≥ 0. Therefore, subtracting the previous inequalities and recalling that R + -R -= R = P -1 , our choices (5.21) of L a (t) and L b (t) give

(5.28) L a (t)R + -L b (t)R -≤ φ B (t) ≤ L b (t)R + -L a (t)R -.
For a suitable matrix that we denote by * , we have

(5.29) φ B (t) = φ B1 (t) * 0 φ B3 (t) .
Hence, (5.22) follows by left multiplying (5.28) through by K, and right multiplying (5.28) through by K . Proposition 5.1 includes the framer result from Lemma 5.2 as a special case. To see why, define H as in (5.19), under the assumptions of Lemma 5.2. Then,

(5.30) 0 ≤ K 0 0 K M ≤ H(t) ≤ K L L K M ,
so we apply Proposition 5.1 with M = H. Note that (5.31)

P = I n 0 I n I n is such that P H(t)P -1 = I n 0 I n I n K(t) L(t) L(t) K(t) I n 0 -I n I n = K(t) -L(t) L(t) 0 K(t) + L(t) .
Then we choose B 1 (t) = K(t) -L(t), and (5.22) reads e Kte (K+L)t -e (K-L)t 2 ≤ φ B1 (t) ≤ e (K+L)t +e (K-L)t

2

(by [START_REF] Mazenc | Design of continuous-discrete observers for time-varying nonlinear systems[END_REF]Lemma A.1]), which agrees with the conclusion of Lemma 5.2. Therefore, Proposition 5.1 includes Lemma 5.2 as a special case. Proposition 5.1 also provides framers that are based on embedding an n dimensional system into a larger system of dimension strictly less than 2n, and which are therefore beyond the scope of [START_REF] Mazenc | Design of continuous-discrete observers for time-varying nonlinear systems[END_REF]. Here is an example where this occurs:

Example 1. Consider the case of the oscillator

(5.32) E(t) = ω(t) 0 1 -1 0 where ω(t) : [0, ∞) → [0, ∞) is a continuous function.
We choose the everywhere Metzler matrix

(5.33) M(t) = ω(t)   0 1 0 0 1 1 2 0 1   .
The eigenvalues of the matrix in (5.33) are 2, i, and -i. Therefore, there is constant matrix P ∈ R 3×3 such that

(5.34) P M(t)P -1 = ω(t)   0 1 0 -1 0 0 0 0 2   .
This gives estimates for φ E by embedding the dynamics in dimension 3, whereas [START_REF] Mazenc | Design of continuous-discrete observers for time-varying nonlinear systems[END_REF] only implies that one gets an estimate for φ E using a system in dimension 4.

The preceding approach for the oscillator can be generalized to many other systems that can be transformed into the triangular form from Proposition 5.1, using similarity transformations. To illustrate the basic ideas of how this can be done, we first consider the constant matrix case (5.35) ẋ = Ax for any matrix A ∈ R n×n . However, we can replace A by a bounded function A(t) which may be uncertain, if we allow the constants ζ and λ that follow to depend only on suitable bounds on A(t) (analogously to Lemma 5.2). We let V = (1.....1) ∈ R n be the column matrix whose entries are all 1's (so V V is a matrix of all 1's), and we denote the rows of A by A i for i = 1, 2, . . . , n, so the system can be written as ẋi = A i x for i = 1, 2, . . . , n.

Let ζ and λ be positive constants to be selected. Then

(5.36) ẋi = (A i + ζV )x -ζV x for 1 ≤ i ≤ n .

We can select ζ so that each entry of each row A i + ζV is nonnegative. Hence, A + ζV V is Metzler. Also, (5.37)

-V ẋ = -λV x + (λV -V A)x
Choosing λ large enough, it follows that all the entries of the vector λV -V A are nonnegative. Next, consider the system

(5.38) ẏi = (A i + ζV )y + ζz, 1 ≤ i ≤ n ż = λz + (λV -V A)y
with state space R n+1 . This system can be written in the form Ż = MZ for a constant Metzler matrix M, and we deduce from (5.36) and (5.37) that any solution of (5.35) is such that y i = x i and z = -V x provide a solution of (5.38). Then we can apply our approach from [START_REF] Mazenc | Design of continuous-discrete observers for time-varying nonlinear systems[END_REF] to designing framers for cooperative systems. Next, consider a bounded function A(t) valued in R n×n . A time-varying analog of the preceding approach provides positive constants ζ and λ such that (5.39)

M(t) = A(t) + ζV V ζV λV -V A(t) λ valued in R (n+1)×(n+1) is Metzler for all t ≥ 0. Let (5.40) P = I n 0 V 1 ∈ R (n+1)×(n+1)
We then define µ by (5.41)

µ(t) = P M(t)P -1 = A(t) ζV 0 λ + nζ .
Thus, we have an n + 1 dimensional upper triangular form that is covered by Proposition 5.1, so we get framers for the solutions of ẋ = A(t)x using solutions of Ẏ = µ(t)Y . In particular, we only need to add a onedimensional dynamic extension. In the next section, we use our framers to design continuous-discrete observers.

6 Framers and Continuous-Discrete Observers

We next discuss our new solution in [START_REF] Mazenc | Design of continuous-discrete observers for time-varying nonlinear systems[END_REF] to the problem of constructing exponentially stable continuous-discrete observers. Let ν > 0 and ν > ν be any two constants, and fix any sequences {t i } and {ν i } in [0, ∞) such that (6.42) t 0 = 0 , and t i+1 = t i + ν i and

ν i ∈ [ν, ν] for all i ∈ N.
The t i 's will serve as the measurement times for ( Using our bounds ν and ν from (6.42), we also assume:

Assumption 5. There exist a constant matrix K ∈ R n×p , a constant κ ∈ (0, 1), and a symmetric positive definite matrix Q ∈ R n×n such that for each constant matrix β ∈ R n×n satisfying (6.46)

β(ρ) ≤ β ≤ β(ρ) for all ρ ∈ [ν, ν],
the inequality β (I -KC) Q(I -KC)β κQ holds.

See Section 7 for ways to verify Assumption 5. In [START_REF] Mazenc | Design of continuous-discrete observers for time-varying nonlinear systems[END_REF], we prove the following, where x(t - 0 ) = x(t 0 ): Theorem 6.1. Let the system (6.43) satisfy Assumptions 4-5 and choose the continuous-discrete system

   ẋ * (t) = A * x * (t) + ϕ * (t, x * (t)) for all t ∈ [t k , t k+1 ) x * (t k ) = x * (t - k ) + P -1 K[y * (t k ) -P x * (t - k )
] for all k ≥ 0. Then the dynamics for the observation error x * -x * is uniformly globally exponentially stable to 0. Remark 1. We can always rewrite the x * dynamics in (6.43) as ẋ * (t) = Ψ * (t, x * (t)), where Ψ * (t, x * ) = A * x * + ϕ * (t, x * ). Hence, we can replace Assumption 4 by the requirements that ϕ * is C 1 and (∂ϕ * /∂x)(t, x) is bounded and select A * = 0. However, different A * 's and P 's produce different conditions in Assumption 5, so it is helpful to consider different possible A * 's. The proof of Theorem 6.1 provides explicit formulas for the functions in the final exponential stability estimate.

Remark 2. The functions (6.45) are obtained by setting K = A -D, K = A + V , and L = V -D in (5.17). They correspond to choosing M(t) = K(t) -L(t) = A + V (t) in Lemma 5.2 for each i, where

K(t) = A V (t), L(t) = V q (t), A V (t) = A+D V (t)+V p (t), D V (t) = diag{v 11 , v 22 , . . . , v nn }, V (t) = [v ij (t)] =
w(t, P x * (t), P x * (t)) for any fixed solutions of the system and observer, V p (t) = max{V N (t), 0}, V q (t) = max{V N (t), 0}-V N (t), and V N (t) = V (t)-D V (t). They provide framers for the dynamics for the error variable x = P (x * -x * ) on [t i , t i+1 ) for all i ≥ 0. However, they can be replaced by our new framers from Proposition 5.1, using arguments from the preceding section. This produces different versions of Theorem 6.1 for different framers, which were not considered in [START_REF] Mazenc | Design of continuous-discrete observers for time-varying nonlinear systems[END_REF].

LMI Formalism

In applications, it can be convenient to check Assumption 5 using LMIs. To see why, we define the functions β and β by (6.45), and we introduce the sets of matrices

F(ρ) = β ∈ R n×n : β ij ∈ β ij (ρ), β ij[ (ρ)
for all i and j .

Following [START_REF] Dinh | Continuous-discrete time observer design for Lipschitz systems with sampled measurements[END_REF], this allows us to rewrite Assumption 5 as an LMI, as in the following result from [START_REF] Mazenc | Design of continuous-discrete observers for time-varying nonlinear systems[END_REF]: Proposition 7.1. Assume that there exist a symmetric positive definite matrix Q ∈ R n×n and a matrix W ∈ R n×p such that

(7.47) Q (Q -W C)β β (Q -W C) Q 0
holds for all β ∈ F(ρ) and all ρ in [ν, ν]. Then there is a constant κ ∈ (0, 1) such that Assumption 5 holds with

K = Q -1 W .
Moreover, if a condition ensuring the existence of an observer in the case of continuous time measurement holds, then Assumption 5 holds provided the ν k 's are small enough. This is made precise in the following proposition in [START_REF] Mazenc | Design of continuous-discrete observers for time-varying nonlinear systems[END_REF]: Proposition 7.2. Let Assumption 4 hold, and define V as in Section 6. Assume that there exist a matrix K 0 ∈ R n×p , a constant κ 0 ∈ (0, 1), and a symmetric positive definite matrix Q ∈ R n×n such that

(7.48) (β 0 -K 0 C) Q + Q(β 0 -K 0 C) -κ 0 Q
holds for all matrices β 0 ∈ R n×n such that A -V ≤ β 0 ≤ A + V . Then there exists a constant ν * > 0 such that for all ν ∈ (0, ν * ], Assumption 5 is satisfied with the choices ν = ν = ν and K = νK 0 .

Conclusions and Open Problems

This tutorial discussed several recent results on continuous-discrete observers. Such observers are useful for finding asymptotic estimators of the states of continuous time systems in situations where only discrete observations of an output of the system are available for measurement. While many results on observers have appeared in the literature, the authors' approaches are distinguished through their use of differential inclusions coupled with the fact that they provide framers for the original systems, i.e., explicit upper and lower bounds for the unknown states (in the vector sense) that hold for all times. The explicit framers in [START_REF] Mazenc | Design of continuous-discrete observers for time-varying nonlinear systems[END_REF] are found by embedding n dimensional systems as subsystems of larger 2n dimensional systems. This tutorial also included new original results on framers, which have not appeared before and which are not being submitted simultaneously for publication elsewhere. Our new framers include the framers in [START_REF] Mazenc | Design of continuous-discrete observers for time-varying nonlinear systems[END_REF] as a special case, but are more general because they make it possible to obtain frames by embedding n dimensional systems as subsystems of larger systems having dimension n + 1. One research direction worth pursuing is to compare the performances of observers that are obtained using the framers from our new Proposition 5.1 with the observers that are obtained from the framers in the previously reported Lemma 5.2. Having sample output values can be viewed as having delayed output measurements, with a timevarying delay. It would be interesting to generalize our approaches to allow delays in the original plant, or cases where the original dynamics is a PDE.

  and F = R -K , and set (5.21) L a (t) = P + e Mt -P -e Mt and L b (t) = P + e Mt -P -e Mt .

  r(b -a) + a)dr .Then the Fundamental Theorem of Calculus givesϕ(t, b) -ϕ(t, a) = w(t, a, b)(b -a) for all t ≥ 0, a ∈ R n , and b ∈ R n . Also, Assumption 4 provides positive constants v ij such that each entry of w = [w ij ] satisfies w ij (t, a, b) ∈ [-v ij , v ij ] for all t ≥ 0, a ∈ R n , b ∈ R n , i ∈ {1, ..., n}, and j ∈ {1, ..., n}. Set D = diag{v 11 , ..., v nn } ∈ R n×n , V = [v ij ] ∈ R n×n, and (6.45) β(ρ) = e (A-D)ρ + 1 2 e (A+D)ρ -e (A+2V -D)ρ and β(ρ) = 1 2 e (A+2V -D)ρ + e (A+D)ρ .

  Theorem 4.1. Assume that there are constants b ij > 0 such that (3.3) holds for all (x, u) ∈ R n × R p and all pairs (i, j) and that φ is C 1 . If there exist a finite set of matrix valued functions S = {M 1 , M 2 , . . . , M } mapping [0, ∞) into R n×n , a symmetric positive definite matrix P in R n×n , and a W ∈ R n×p such that

	(4.10)	A

δ (e) ⊆ Conv{M i (δ)e : 1 ≤ i ≤ } holds for all e ∈ R n , and such that the matrix inequality

(4.11)

  Mt ≤ P + Rφ B (t)P ≤ P + e Mt and P -e Mt ≤ P -Rφ B (t)P ≤ P -e Mt . Mt -P -e Mt ≤ φ B (t)P ≤ P + e Mt -P -e Mt and therefore also P + e Mt -P -e Mt R + ≤ φ B (t)P R

	Since R = P -1 , it follows that
	(5.25)	P + e

-1 , so Assumption 2 gives (5.23) e Mt ≤ Rφ B (t)P = φ M (t) ≤ e Mt . Consequently, since P + ≥ 0 and P -≥ 0, we have (5.24) P + e + ≤ P + e Mt -P -e Mt R + and (5.26)

  (t k ), t ∈ [t k , t k+1 ), k ∈ Nwith discrete measurements, where x * and y * are valued in R n and R p respectively. Assume: Assumption 4. There is an invertible matrix P ∈ R n×n such that the matrix A = P A * P -1 is Metzler. Also, ϕ * is C 1 , and (∂ϕ * /∂x)(t, x) is bounded.

	6.43) ẋ * (t) = A Set ϕ(t, x) = P ϕ * (t, P -1 x), C = C * P -1 , and
	(6.44)	w(t, a, b) =

* x * (t) + ϕ * (t, x * (t)) y * (t) = C * x *

* Supported by ANR LIMICOS Contract 12 BS03 005 01 (Andrieu) and NSF Grant 1102348 (Malisoff). † Mazenc is with EPI DISCO INRIA-Saclay, the


Denis Barthou

Albert Cohen
email: acohen@prism.uvsq.fr

Jean-Franc ¸ois Collard

Maximal Static Expansion

Keywords: Expansion of data structures, privatization, single assignment

Memory expansions are classical means to extract parallelism from imperative programs. However, for dynamic control programs with general memory accesses, such transformations either fail or require some run-time mechanism to restore the data flow. This paper presents an expansion framework for any type of data structure in any imperative program, without the need for dynamic data flow restoration. The key idea is to group together the write operations that participate in the flow of the same datum. We show that such an expansion boils down to mapping each group to a single memory cell. We give a practical algorithm for code transformation. This algorithm, however, is valid for (possibly non-affine) loops over arrays only.

Introduction

Data dependences are known to hamper automatic parallelization of imperative programs and their efficient compilation on modern superscalar or VLIW processors. A general method to tackle this problem is to disambiguate memory accesses and to assign distinct memory cells to non-conflicting writes, i.e. to expand data structures. In parallel processing, expanding a datum also allows to place one copy of the datum on each processor, enhancing parallelism. This technique is known as array privatization [START_REF] Tu | Automatic array privatization[END_REF][START_REF] Maydan | Array dataflow analysis and its use in array privatization[END_REF][START_REF] Creusillet | Array Region Analyses and Applications[END_REF] and is extremely important to parallelizing and vectorizing compilers [START_REF] Levine | A comparative study of automatic vectorizing compilers[END_REF][START_REF] Lee | Parallelizing Compilers: Implementation and Effectiveness[END_REF]. A similar technique is register or variable renaming.

In the extreme case, each memory cell is written at most once, and the program is said to be in single assignment (SA) form. Unfortunately, when the control flow cannot be predicted at compiletime, some run-time computation is needed to preserve the original data flow: In the static single-assignment framework, φ-functions may be needed to "merge" multiple reaching definitions, i.e. possible data definitions due to several incoming control paths [START_REF] Cytron | An efficient method of computing static single assignment form[END_REF][START_REF] Cytron | Efficiently computing static single assignment form and the control dependence graph[END_REF]. Such φ-functions may be an overhead at run-time, especially for non-scalar data structures or when replicated data are distributed across processors. We are thus looking for a static expansion, i.e. an expansion of data structures that does not need a φ-function.

To be published in the proceedings of PoPL'98: The twenty-fifth annual ACM SIGACT-SIGPLAN symposium on Principles of Programming Languages, January 19-21 1998, San Diego, CA.

(Notice that according to our definition, an expansion in the static single assignment framework may not be static.) The goal of this paper is to automatically find the static expansion which expands all data structures as much as possible, i.e. the maximal static expansion. Maximal static expansion may be considered as a trade-off between parallelism and memory usage.

We present an algebraic framework to derive the maximal static expansion. The input of this framework is the (perhaps inaccurate) output of a data-flow analysis, so our method is "optimal" with respect to the precision of this analysis. Our framework is valid for any imperative program, without restriction-the only restrictions being those of your favorite data-flow analysis. We then present an algorithm to construct the maximal static expansion for programs with arrays only, but where subscripts and control structures are unrestricted.

The paper is organized as follows: Section 2 studies motivating examples showing what we want to achieve. Section 3 formally states what (maximal) static expansion is, and Section 4 presents a general framework to solve this problem. This framework is applied in Section 5 to derive an algorithm for maximal static expansion. Section 6 applies this algorithm to the motivating examples, before our conclusion.

Motivating Examples

The general framework presented in this paper is valid for any imperative programs. However, the three examples we study in this section are basically loop nests over arrays (mainly because our own analysis [START_REF] Barthou | Fuzzy array dataflow analysis[END_REF] is restricted to such programs).

Definitions

For any statement, the iteration vector is the vector built from surrounding loop counters. The iteration domain is the set of values the iteration vector takes during program execution. For instance, the iteration domain of T in Figure 1

is Dom(T ) = {i : 1 ≤ i ≤ N}.
Each iteration of a loop spawns instances, called operations, of statements included in the loop body. In the example program, the for loop on i yields N instances of T , denoted by 1, T , . . . , N, T . Moreover, we introduce artificial integer counters for while loops. For instance, operations of S in Figure 1 are labeled i, w, S , with 1 ≤ i ≤ N and w ∈ N.

The execution order on operations is denoted by ≪. Each operation i, T assigns a new value to variable x. In turn, statement S assigns x an undefined number of times (possibly zero). The value read in x by statement R is thus defined either by T , or by some instance of S, in the same iteration of the for loop (the same i). Therefore, if the expansion assigns distinct memory cells to i, T and to instances of i, w, S , how could operation i, R "know" which memory cell to read from?

To formalize this problem, we use a data-flow analysis to describe where values are defined and where they are used. The intuitive picture is that a datum flows from its source to the sink . We assume that the data-flow analysis works at operation level. Moreover the analysis may be more or less accurate: When the exact source of a read operation cannot be predicted at compile time, we suppose that it returns a set of possible sources for this read. This set is a conservative approximation of the source.

We may thus call σ the function mapping a read operation to its set of sources. Applied to the example in Figure 1, it tells us that the set of sources σ( i, w, S ) of an operation i, w, S is:

σ( i, w, S ) = if w > 0 then { i, w -1, S } else { i, T } (1)
And the set of sources σ( i, R ) of an operation i, R is:

σ( i, R ) = { i, T } ∪ { i, w, S : w ≥ 0}, ( 2
)
where w is the (arbitrary) counter of the while-loop.

Let us try to expand scalar x. One way is to convert the program into SA, making T write into x'[ i ] and S into x''[ i, w ]: Then, each memory cell is assigned to at most once, complying with the definition of SA. However, what should right-hand sides look like now? A brute-force application of (2) yields the program in Figure 2. While the right-hand side of S only depends on w, the right-hand side of R depends on the control flow, thus needing a function similar to a φ-function in the SSA framework (even if, on this introductory example, the φ-function would be very simple) [START_REF] Griebl | Generation of synchronous code for automatic parallelization of while loops[END_REF].

The aim of this paper is to expand x as much as possible in this program but without having to insert φ-functions.

A possible static expansion is to uniformly expand x into x[i] and avoid output dependencies between distinct iterations of the for loop. The resulting maximal static expansion of this example is given by Figure 3. It has the same degree of parallelism and is simpler than the program in single-assignment. Notice that it should be easy to adapt the array privatization techniques by Maydan et al. [START_REF] Maydan | Array dataflow analysis and its use in array privatization[END_REF] to handle the program in Figure 1; This would tell us that x can be privatized along i. However, we want to do more than privatization along loops, as illustrated in the following examples.

for i = 1 to N do T x'[i] = .

Second Example: Array Expansion

Let us give a more complex example; We would like to expand array A in the program in Figure 4.

Since T always executes when j equals N, a value read by i, j, S , j > N is never defined by an instance i ′ , j ′ , S of S with j ′ ≤ N. Figure 4 describes the data-flow relations between S instances: An arrow from (i ′ , j ′ ) to (i, j) means that instance (i ′ , j ′ ) defines a value that may reach (i, j).

Formally, the source of one instance of statement S is:

σ( i, j, S ) = if j ≤ N then { i ′ , j ′ , S : 1 ≤ i ′ ≤ 2N ∧ 1 ≤ j ′ < j ∧ i ′ -j ′ = i -j} else { i ′ , j ′ , S : 1 ≤ i ′ ≤ 2N ∧ N < j ′ < j ∧ i ′ -j ′ = i -j} ∪ { i ′ , N, T : 1 ≤ i ′ < i ∧ i ′ = i -j + N} (3)
Because sources are non-singleton sets, converting this program to SA form would require run-time computation of the memory location read by S.

However, we notice that the iteration domain of S may be split into disjoint subsets by grouping together operations involved in the same data flow. These subsets build a partition of the iteration domain. Each subset may have its own memory cell, a cell that will not be written nor read by operations outside the subset. The partition is given in Figure 5.a.

Using this property, we can duplicate only those elements of A that are used twice. These are all the array elements A The maximal static expansion is shown in Figure 6. Notice that this program has the same degree of parallelism as the corresponding single-assignment program, without the run-time overhead.

[ c ], 1 + N ≤ c ≤ 3N -1.

Third Example: Non-Affine Array Subscripts

Consider the program in Figure 7.a, where foo and bar are arbitrary subscripting functions 2 . Since all array elements are assigned by T , the value read by R at the i th iteration must have been produced by S or T at the same iteration. The data-flow graph is similar to the first example:

σ( i, R ) = { i, S } ∪ { i, j, T : 1 ≤ j ≤ N}. ( 4
)
The maximal static expansion adds a new dimension to A subscripted by i. It is sufficient to make the first loop parallel.

These examples show the need for an automatic static expansion technique. We present in the following section a formal definition of expansion and a general framework for maximal static expansion. We then describe an expansion algorithm for arrays that yields the expanded programs shown above. Notice that it is easy to recognize the original programs in their expanded counterparts, which is a practical property of our algorithm.

Related work

If the input program is built of nested for loops with affine bounds and accesses arrays with affine subscripts, one can find a static expansion which is also in single-assignment form. Feautrier [START_REF] Feautrier | Dataflow analysis of scalar and array references[END_REF] coined the term static control programs for this class of programs.

In the case of programs with general control and unrestricted arrays subscripts, array data-flow analyses are approximate [START_REF] Collard | Fuzzy array dataflow analysis[END_REF][START_REF] Barthou | Fuzzy array dataflow analysis[END_REF][START_REF] Wonnacott | Nonlinear array dependence analysis[END_REF][START_REF] Wonnacott | Constraint-Based Array Dependence Analysis[END_REF]: Several writes may be the unique definition of a given value, but the analysis cannot tell. [START_REF] Griebl | Generation of synchronous code for automatic parallelization of while loops[END_REF] describes how to obtain a single-assignment program to the price of dynamic restoration of data flow.

Many studies are related to array privatization. As hinted above, Maydan et al. [START_REF] Maydan | Array dataflow analysis and its use in array privatization[END_REF] proposed an algorithm to privatize arrays. However, their method only applies to static control programs. Tu and Padua [START_REF] Tu | Automatic array privatization[END_REF] proposed a privatization technique for a very large class of programs. But it resorts to dynamic restoration of data flow. Another accurate approach using array regions has been described by Creusillet [START_REF] Creusillet | Array Region Analyses and Applications[END_REF]. Her method avoids the cost of a dynamic restoration and copies back the privatized elements into the original arrays.

However, privatization only detects parallelism along the enclosing loops; It is thus less powerful than general array expansion techniques. Indeed, the example in Section 2.3 shows that our method not only may expand along diagonals in the iteration space but may also do some "blocking" along these diagonals.

Static Expansion

Let Ω be the set of all operations in the program, f the function mapping operations to memory cells they write into, and be the set of all writes. We still use σ to denote the function mapping a read operation to its set of possible sources. Notice that σ may also be seen as a relation between read and write operations. Let f ′ be the expansion, that is the new function, after program transformation, mapping operations to the memory cells they write into.

W ⊆ Ω real A[1..4*N-1,1..2] for i = 1 to 2*N do for j = 1 to 2*N do {expansion of statement S} if -2*N+1 <= i-j <= -N then if ... then S 1 A[i-j+2*N,1] = ... A[i-j+2*N,1] ... end if elsif -N+1 <= i-j <= N-1 then if j <= N then if ... then S 2 A[i-j+2*N,1] = ... A[i-j+2*N,1] ... end if else if ... then S 3 A[i-j+2*N,2] = ... A[i-j+2*N,2] ... end if end if else if ... then S 4 A[i-j+2*N,1] = ... A[i-j+2*N,1] ... end if end if {expansion of statement T } T if j = N then A[i+N,2] = ...

end if end for end for

Let us consider two operations u and v belonging to the same set of possible sources of some read r. If they both write in the same memory cell ( f (u) = f (v)) and if we assign two distinct memory cells to u and v ( f ′ (u) = f ′ (v)), then a φ-function is needed to restore the data flow since we do not know which of the two cells has the value needed by r. Static expansion enforces f ′ (u) = f ′ (v).

Definition 1 (Static expansion) A static expansion is a mapping f ′ from operations to memory cells such that

∀u, v : (∃r, u ∈ σ(r) ∧ v ∈ σ(r) ∧ f (u) = f (v)) =⇒ f ′ (u) = f ′ (v).
Because the sources of a read are mapped to the same memory cell by f ′ , static expansion preserves the original data-flow graph. Notice also that, according to this definition, even a constant function on W is a static expansion. Because we are interested in maximizing the memory expansion, the range of a "good" static expansion should be as large as possible. In other words, such an expansion should be constant on sets as small as possible:

Definition 2 (Maximal static expansion) A static expansion f ′ is maximal on the set of operations W if, for any static expansion f ′′ , ∀u, v ∈ W : f ′ (u) = f ′ (v) =⇒ f ′′ (u) = f ′′ (v).
Intuitively, if f ′ is maximal, then f ′′ cannot do better: it maps two writes to the same memory cell when f ′ does. We need to characterize the sets of operations on which a maximal static expansion f ′ is constant, i.e. the equivalence classes of

the relation {u, v ∈ W : f ′ (u) = f ′ (v)}. The set of theses classes is denoted by W f ′ .
The number of memory cells after maximal static expansion is thus equal to the cardinal of W f ′ . However, this hardly gives us an expansion scheme, because this result does not tell us how much each individual memory cell should be expanded. The purpose of Section 4 is to give a similar result for each memory cell c used in the original program. This result appears in Theorem 1. This theorem is then used to give a practical expansion scheme.

Expansion Scheme

Let us define the relation:

uR ⊑ ⇐⇒ ∃∇, ⊓ ∈ σ(∇) ∧ ⊑ ∈ σ(∇). (5)
σ is itself a relation on Ω × Ω and the reciprocal relation is denoted by σ

-1 . Therefore, uR ⊑ ⇐⇒ ⊓ ∈ σ(σ -∞ (⊑)), i.e., R = σ•σ -∞ . Relation R is obviously symmetric. Definition 1 requires that a static expansion ′ verifies f ′ (u) = f ′ (v) when f (u) = f (v) and uR ⊑. Given u, v and w in W , if f (u) = f (v) = f (w), uR ⊑ and vR ⊒ then f ′ (u) = f ′ (v) = f ′ (w). Therefore, given u ∈ W , f ′ is constant on the set of all v ∈ W such that f (u) = f (v) and uR * v,
R * being the transitive closure of R . We may give an equivalent definition of a static expansion:

Definition 3 A static expansion is a mapping f ′ from operations to memory cells such that ∀u, v : uR * v ∧ f (u) = f (v) =⇒ f ′ (u) = f ′ (v).
We now characterize any maximal static expansion in terms of R * and f :

Lemma 1 f ′ is a maximal static expansion if and only if ∀u, v ∈ W : uR * v ∧ f (u) = f (v) ⇐⇒ f ′ (u) = f ′ (v). (6)
Sufficient condition-the "if" part

Let f ′ be a mapping s.t. ∀u, v ∈ W : f ′ (u) = f ′ (v) ⇐⇒ uR * v∧ f (u) = f (v). By definition, f ′ is a static expan- sion. Let us show that f ′ is maximal. Suppose that for u, v ∈ W : f ′ (u) = f ′ (v). (6) implies uR * v and f (u) = f (v).
Thus, from Definition 3, any static expansion

f ′′ sat- isfies f ′′ (u) = f ′′ (v). Therefore, f ′ (u) = f ′ (v) =⇒ f ′′ (u) = f ′′ (v), so f ′ complies with Definition 2.
Necessary condition-the "only if" part Let f ′ be a maximal static expansion. Because f ′ is a static expansion, we only have to prove that ∀u, v ∈

W : f ′ (u) = f ′ (v) =⇒ uR * v ∧ f (u) = f (v).
On the one hand,

f ′ (u) = f ′ (v) =⇒ f (u) = f (v) be- cause f is a static expansion.
On the other hand, assume f ′ (u) = f ′ (v) and ¬uR * v. We show that it contradicts the maximality of f ′ : Let f ′′ (w) = f ′ (w) when ¬uR * w, and

f ′′ (w) = c when uR * w, with c = f ′ (u). f ′′ is a static expansion: By construction, f ′′ (u ′ ) = f ′′ (v ′ ) for any u ′ and v ′ such that u ′ R * v ′ . The contra- diction comes from the fact that f ′′ (u) = f ′′ (v).

✷

Let us define M = f (W ) the set of all memory cells accessed by write and for c M, W (c) = {u ∈ W : (u) = c} the set of operations writing into c. Given c ∈ M, the previous lemma entails that a static expansion f ′ is maximal iff

∀u, v ∈ W (c) : f ′ (u) = f ′ (v) ⇐⇒ uR * v.
classes of R * in W (c) are exactly the sets we are looking for:

Theorem 1 The sets on which a maximal static expansion f ′ is constant are described by:

W f ′ = c∈M W (c) R * (7)
The equivalence classes defined in this theorem gives the partition intuitively found in Section 2, and the expansion factor of each individual memory cell c is Card( W (c) R * ). Consider for instance A[0] in Figure 5.a. The instances of S that belong to W (A[0]) are on the main diagonal {(i, j) : 1 ≤ i, j ≤ 2N ∧ i = j = 0}. R * partitions these operations in exactly the two subsets depicted in the figure.

To generate the transformed code, one has to remember which equivalent class an operation belongs to: Let ϕ be the function mapping an operation u to a representative of its equivalence class. One may label each element of W (c) R * , or equivalently, label each element of ϕ(W (c)). Such a labeling scheme is obviously arbitrary, but all programs transformed using our method are equivalent up to a permutation of these labels. We denote by ν(u) the label we choose for the elements of ϕ(W ( f (u))). Then, f ′ = ( f , ν).

Our expansion scheme depends on the transitive closure calculator and on the part calculating W (c). We would like to stress the fact that the expansion produced is static and maximal with respect to the results yielded by these parts, whatever their accuracy:

• The exact transitive closure may be too complicated and may therefore be over-approximated. The expansion factor of a memory cell c is then lower than Card( W (c) R * ). However, the expansion remains static and is maximal with respect to the transitive closure given to the algorithm.

• The sets W (c) may not be known precisely at compile-time.

(For instance, when data structures are arrays with non-affine subscripts.) One may use some approximation W (c) instead, such that W (c) ⊆ W (c), and expand c into as many cells as elements in W (c) R * . However, an operation u may then beto two distinct classes of W (c) R * and W (c ′ ) R * , c = c ′ , that is, have several representatives and be associated to different class labels. To avoid this pitfall, we enforce the same labels for all classes including u: We first label all classes of W R * , which in turn gives labels to the classes of all W R * . The drawback of this method is that some memory cells not used during program execution may be allocated. The reasons are that we cannot know statically which cells will be referred to, and that the set of numbers labeling the classes of a given W (c) R * may not be dense.

The maximal static expansion scheme given above works for any imperative programs. More precisely, you may expand any imperative program using maximal static expansion, provided that a data-flow analysis technique can handle it (at operation level) and that transitive closure computation, relation composition, intersection, and so on, are feasible in your framework.

Expanding scalars and arrays is done by renaming the variables and adding new dimensions to arrays; However, no straightforward expansion exists for trees, graphs, dynamic data structures with pointers . . . In the general case, appropriate expansion "rules" must be defined-depending on both the data and control structures.

We give below an algorithm to construct expanded codes for loops nests and arrays only 3 .

Before giving the algorithm, we would like to focus on two important points:

• The algebraic view given in this section considered each memory cell c in turn. Obviously, since the number of memory cells brought into play in a program is often unknown or parameterized, a naive application of this view would not practical. Our method gives a solution parameterized by the identity of the cell c, so its complexity does not grow with Card(M).

• The definitions given in Section 3 and the expansion scheme are valid for any class of imperative programs. The only restrictions and limitations are those of the data-flow analysis and of the algorithm to compute transitive closures.

In the sequel, since we apply our own array data-flow analysis framework to maximal static expansion, we inherit its syntactical restrictions: Data structures are scalars and arrays; Pointers are not allowed. Loops, conditionals and array subscripts are unrestricted.

An Algorithm for Loop Nests

Using a data-flow analysis such as FADA [START_REF] Barthou | Fuzzy array dataflow analysis[END_REF], the data-flow graph is described by systems of affine inequalities over iteration variables and structure parameters. Our algorithm then reduces to well known transformations on affine integer polyhedra, most of them being implemented in Omega [START_REF] Pugh | The omega test: a fast and practical integer programming algorithm for dependence analysis[END_REF]. We present below the expansion algorithm for all accesses to a given array A.

Input: The data-flow graph as an affine relation σ between reads and their reaching definitions (the sources).

Output: The target expanded code.

1. Compute R = σ • σ -∞ . (This boils down to eliminating r in (5).)

2. If R is not transitive, compute R * with Omega's transitive closure operator. Because the transitive closure of an affine relation is not necessary affine, the result may be an upperapproximation. See [START_REF] Kelly | Transitive closure of infinite graphs and its applications[END_REF] for details. This approximation is a conservative one, but may hide an interesting possible static expansion. Using Omega, R * is described as a mapping from u to u ( u being the class of u for relation R * : u = {v ∈ W : uR * v}).

3.

In each class u, pick a single, arbitrary element. This chosen element is now considered as the representative ϕ(u). How do we pick this element? As long as the element we pick is unique, any method is fine. Let us choose the minimum according to lexicographical order (which a case of overkill).

4. Are all subscript functions affine?

Yes Let us consider c = A [ x ]. W (A [ x ]) is the union of { i, S : i ∈ Dom(S) ∧ f ( i, S ) = x} over all statements S writing into A. Compute ϕ(W (A [ x ])), which is a set of representa- tives of W (A[x]) R * .
Give a number to each element in the set of representatives. No Compute ϕ(W ). Give a number to each element in the set of representatives.

If an element in the set of representatives is itself a parameterized affine set of operations, labeling boils down to scanning exactly once all the integer points in the set, which can be done using classical techniques [START_REF] Ancourt | Scanning polyhedra with DO loop[END_REF][START_REF] Kelly | Code generation for multiple mappings[END_REF]. In both cases, u has a single representative and is therefore mapped to a unique label ν(u).

Code generation is then straightforward: any reference

A[ f (u) ] in the left hand side is transformed into A[ f (u), ν(u) ].
For any reference in the right hand side, one has to find the label of the source of the read. That is, any read

A[ g(u) ] is transformed into A[ g(u), ν(σ(u)) ].
(Recall that σ(u) is a set, mapped by construction of ν to a single label ν(σ(u)).) When ν is a conditional whose predicate is affine w.r.t. loop counters, then the conditional can be taken out of A's subscript.

The size declaration

A[ ... ] of A is transformed into A[ ..., max S max u∈Dom(S) ν(u) ] 4 .
Computing the Lexicographical Minimum Let us call u the equivalence class of u for relation R * . The lexicographical minimum of u is: min

≪ ( u) = v s.t. uR * v ∧ ( ∃w : uR * w ∧ w ≪ v)
This definition may be simplified in writing ≪ as a relation between operations: 4 Arrays usually have to be rectangular; Therefore A ν(u) [ f (u)] may be a better renaming. Consider for instance the expanded version of example 2: Expanding A into A1 and A2 would require 6N -2 array elements instead of 8N -2 in Figure 6.

≪ = {(u, v) : u ≪ v}.

Thus, min

≪ ( u) = (R * \ (≪ •R * ))(u) (8)
Complexity For each array in the source program, the algorithm proceeds as follows:

• Compute the reciprocal relation σ -1 of σ. This is different from computing the inverse of a function and barely consists in a swap of the two arguments of σ.

• Composing two relations σ and σ ′ boils down to eliminating y in xσy ∧ yσ ′ z.

• Computing the exact transitive closure of R is quite expensive. Kelly et al. [START_REF] Kelly | Transitive closure of infinite graphs and its applications[END_REF] do not give a formal bound on the complexity of their algorithm, but their implementation in the Omega toolkit proved to be efficient if not concise. Notice again that the exact transitive closure is not necessary for our expansion scheme to be correct.

Moreover, R happens to be often transitive in practice. In our implementation, this is first checked before triggering the computation of the closure by testing whether the difference

(R • R ) \ R is empty.
In all three examples, the relation is already transitive.

• In the algorithm above, ϕ is a lexicographical minimum. This clearly is a bad idea, because the expansion scheme just needs a way to pick one element per equivalence class. Computing the lexicographical minimum is expensive a priori, but was easy to implement in our first prototype.

• Finally, numbering classes is costly only when we have to scan a polyhedral set of representatives in dimension greater than 1. In practice, we only had intervals on the examples we tried.

Implementation The maximal static expansion is implemented in C++ on top of the Omega library. Figure 8 summaries the computation times for the three examples (on a Sun SPARCstation 5). These results do not include the computation times for data-flow analysis and code generation. Moreover, computing the class representatives is relatively fast; It validates our choice to compute function ϕ (mapping operations to their representatives) using a lexicographical minimum. The intuition behind these results is that the computation time mainly depends on the number of affine constraints in the data-flow analysis relation.

Our only concern so far would be to find a way to approximate the expressions of transitive closures when they become large.

Back to the examples

This section applies our algorithm to the motivating examples, using the Omega Calculator [START_REF] Pugh | The omega test: a fast and practical integer programming algorithm for dependence analysis[END_REF] as a tool to manipulate affine relations.

First Example

Let us consider the source program at Figure 1. Using the Omega Calculator text-based interface, we describe a step-by-step execution of the expansion algorithm. We have to code operations as integer-valued vectors. An operation i, S s is denoted by vector [i,..,s], where [..] possibly pads the vector with zeroes. We number T, S, R with 1, 2, 3 in this order, so i, T , i, j, S and i, R are written [i,0,1], [i,j,2] and[i,0,3], respectively.

From ( 1) and ( 2), we construct the source relation S:

# S := {[i,0,2]->[i,0,1] : 1<=i<=N} # union {[i,w,2]->[i,w-1,2] : 1<=i<=N && 1<=w} # union {[i,0,3]->[i,0,1] : 1<=i<=N} # union {[i,0,3]->[i,w,2] : 1<=i<=N && 0<=w};
Step 1. Computing R is straightforward:

# S' := inverse S; # R := S(S'); # R; {[i,0,1]->[i,0,1] : 1<=i<=N} union {[i,w,2]->[i,0,1] : 1<=i<=N && 0<=w} union {[i,0,1]->[i,w',2] : 1<=i<=N && 0<=w'} union {[i,w,2]->[i,w',2] : 1<=i<=N && 0<=w' && 0<=w}
Step 2. R is already transitive, no closure computation is thus necessary.

Step 3. Let us choose ϕ(u) as the first executed operation in class u (the least operation according to the sequential order): ϕ(u) = min ≪ ({u ′ : u ′ R * u}).

To compute the lexicographical minimum, let us rewrite its definition using the Omega Calculator syntax. We thus describe ϕ by a relation of the form:

ϕ([ i, w, s ]) = [ i ′ , w ′ , s ′ ] s.t. [ w, s ], [ i ′ , w ′ , s ′ ] ∈ W (x) ∧ [ i, w, s ] R * [ i ′ , w ′ , s ′ ] ∧ ( ∃ [ i ′′ , w ′′ , s ′′ ] ∈ W (x) : [ i, w, s ] R * [ i ′′ , w ′′ , s ′′ ] ∧ [ i ′ , w ′ , s ′ ] ≪ [ i ′′ , w ′′ , s ′′ ])
Since [ i, w, s ] ∈ W (x) is always verified in this example, we may simplify this expression in using (8):

ϕ([ i, w, s ]) = (R * \ (≪ •R * ))([ i, w, s ]). (9)
The result of this computation is:

∀i, w, 1 ≤ i ≤ N, w ≥ 1 : ϕ( i, T ) = i, T , ϕ( i, w, S ) = i, T .
Step 4. Since we have only one memory cell, W (x) = W .

Computing ϕ(W (x)) yields N operations of the form i, T . Maximal static expansion of accesses to variable x requires N memory cells. i is an obvious label:

∀i, w, 1 ≤ i ≤ N, w ≥ 1 : ν( i, w, S ) = ν( i, T ) = i. ( 10
)
Step 5. All left-hand side references to x are transformed into x[i]; All references to x in the right hand side are transformed into x[i] too since their sources are instances of S or T for the same i. The expanded code is thus exactly the one found intuitively in Figure 3.

Step 6. The size declaration of the new array is x[1..N].

Second Example

We now consider the source program in Figure 4. Operations i, j, S and i, N, T are denoted by [i,j,1] and [i,N,2], respectively. From (3), the source relation S is defined as:

# S := {[i,j,1]->[i', j',1] : 1<=i,i'<=2N # && 1<=j'<j<=N && i'-j'=i-j} # union {[i,j,1]->[i',N,2] : 1<=i,i'<=2N # && N<j<=2N && i'=i-j+N} # union {[i,j,1]->[i',j',1] : 1<=i,i'<=2N # && N<j'<j<=2N && i'-j'=i-j};
Step 1. As in the first example, we compute relation R using Omega:

# S' := inverse S; # R := S(S'); # R; {[i,j,1]->[i',j-i+i',1] : 1<=i<=2N-1 && 1<=j<N && 1<=i'<=2N-1 && i<j+i' && j+i'<N+i} union {[i,j,1]->[i',j-i+i',1] : N<j<=2N-1 && 1<=i<=2N-1 && 1<=i'<=2N-1 && N+i<j+i' && j+i'<2N+i} union {[i,N,2]->[i',N-i+i',1] : 1<=i<i'<=2N-1 && i'<N+i} union {[i,j,1]->[N+i-j,N,2] : N<j<=2N-1 && i<=2N-1 && j<N+i} union {[i,N,2]->[i,N,2] : 1<=i<=2N-1}
Step 2. We compute R * . Figure 5.a shows the equivalence classes of R * .

Step 3. We compute ϕ(u) as a relation similar to (9), using Omega. The result follows:

∀i, j, 1 ≤ i ≤ 2N, 1 ≤ j ≤ N, j -i ≥ 0 ϕ( i, j, S ) = 1, j -i + 1, S ∀i, j, 1 ≤ i ≤ 2N, 1 ≤ j ≤ N, j -i < 0 ϕ( i, j, S ) = i -j + 1, 1, S ∀i, j, 1 ≤ i ≤ 2N, N + 1 ≤ j ≤ 2N, j -i ≥ N ϕ( i, j, S ) = 1, j -i + 1, S ∀i, j, 1 ≤ i ≤ 2N, N + 1 ≤ j ≤ 2N, j -i < N ϕ( i, j, S ) = i -j + N, N, T ∀i, j, 1 ≤ i ≤ 2N ϕ( i, N, T ) = i, N, T Step 4. W (c) = { i, j, S : i -j + 2N = c} ∪ { c -N, N, T }.
Let us compute the representatives for W (c):

1 ≤ c ≤ N : ϕ(W (c)) = 1, 1 -c, S N + 1 ≤ c ≤ 3N -1 : ϕ(W (c)) = { c + 1, 1, S , c + N, N, T } 3N ≤ c ≤ 4N -1 : ϕ(W (c)) = c + 1, 1, S
This result shows three intervals of constant cardinality of W (c) R * ; They are described in Figure 5.b. A labeling can be found mechanically. If ij + 2N ≤ N or ij + 2N ≥ 3N, there is only one representative in ϕ(W (ij + 2N)), thus ν( i, j, S ) = 1. If N + 1 ≤ ij + 2N ≤ 3N -1, there are two representatives; Thus we define ν( i, j, S ) = 1 if j ≤ n, ν( i, j, S ) = 2 if j > N, and ν( i, N, T ) = 2.

Step 5. The static expansion code appears in Figure 6. As hinted in Section 5, conditionals in ν have been taken out of array subscripts.

Step 6. Array A is allocated as A[1..4*N-1,1..2].

Third Example: Non-Affine Array Subscripts

We come back to the program in Figure 7.a. Operations i, j, T , i, S and i, R are written [i,j,1], [i,0,2] and [i,0,3]. From (4), we build the source relation as follows:

# S := {[i,0,3]->[i,j,1] : 1<=i,j<=N} # union {[i,0,3]->[i,0,2] : 1<=i<=N};

Step 1. Step 2. R is already transitive: R = R * .

Step 3. We compute ϕ(u) as a relation similar to [START_REF] Kelly | Code generation for multiple mappings[END_REF].

∀i, 1 ≤ i ≤ N : ϕ( i, S ) = i, 1, T ∀i, j, 1 ≤ i ≤ N, 1 ≤ j ≤ N : ϕ( j, T ) = 1, T that every i, j, T operation is in relation with i, 1, T .

Step 4. Since some subscripts are not affine, we cannot compute at compile-time the exact sets W (A[ x ]) of operations writing in some cell A[ x ]. Therefore, we compute ϕ(W ):

ϕ(W ) = { i, 1, T }.
We can use i to label these representatives; Thus the resulting ν function is: ν( i, S ) = ν( i, j, T ) = i.

Step 5. Using this labeling, all left hand side references to A[...] become A[..., i] in the expanded code. Since the source of i, R is an instance of S or T at the same iteration i, the right hand side of R is expanded the same way. Expanding the code thus leads to the intuitive result given at Figure 7.b.

Step 6. The size declaration of A is now A[1..N,1..N].

Conclusion

Expanding data structures is a classical optimization to cut memorybased dependences. However, the generated code has to ensure that all reads refer to the correct memory cell. When control flow is dynamic, the main drawback of such methods is therefore that some run-time computation has to be done to decide the identity of the correct memory cell. This paper presented a new and general expansion framework: A cell can be expanded at most as many times as there are classes of independent (as far as data-flow is concerned) writes. A practical algorithm was given and applied to real-life loop nests accessing arrays.

Interestingly enough, the framework does not require any precision level of the data-flow analysis, nor does it require the closure computation to be exact. Conservative approximate results are fine as well, the only drawback being a probable loss in static expansion. However, we cannot do any better, and in accordance to our definition, the static expansion we derive is still maximal. When the data-flow analysis and/or the transitive closure tool give poor results, our expansion scheme does not fail but degrades gracefully.

Future work will study the application of the framework to a wider class of problems. We also intend to enhance the algorithm so as to handle pointer-based data structures and recursive programs.

1 2. 2 Figure 1 :

 121 Figure 1: First example.

Figure 4 :

 4 Figure 4: Second example.

Figure 5 Figure 5 :

 55 Figure 5.a.

Figure 6 :

 6 Figure 6: Maximal static expansion for the second example.

  Figure 7.a: Source program.

Figure 7 :

 7 Figure 7: Third Example.

Figure 8 :

 8 Figure 8: Computation times, in milliseconds.

  # S' := inverse S; # R := S(S'); # R;{[i,j,1]->[i,j',1] : 1<=i<=N && 1<=j<=N && 1<=j'<=N} union {[i,0,2]->[i,j',1] : 1<=i<=N && 1<=j'<=N} union {[i,j,1]->[i,0,2] : 1<=i<=N && 1<=j<=N} union {[i,0,2]->[i,0,2] : 1<=i<=N}

For instance, Horn and Schunck's algorithm to perform 3D Gaussian smoothing by separable convolution.

A[foo(i)] stands for an array subscript between 1 and N, "too complex" to be analyzed at compile-time.

This is mainly due to the fact that our implementation of the expansion scheme is based on our own data-flow analysis, which is restricted to such programs.

Acknowledgments The first two authors are supported by the French Ministre de l'Enseignement Suprieur et de la Recherche (MESR) and the third by the Centre National de la Recherche Scientifique (CNRS). All authors are, in addition, supported by INRIA project AAA and the German-French ProCoPe program.

We would like to thank Paul Feautrier, Mx Geigl, Martin Griebl and Vincent Lefebvre for fruitful discussions on this topic.