Dina Irofti
email: dina.irofti@l2s.centralesupelec.fr

Islam Boussaada
email: islam.boussaada@l2s.centralesupelec.fr

Silviu-Iulian Niculescu
email: silviu.niculescu@l2s.centralesupelec.fr

Geometric vs. algebraic approach: A study of double imaginary characteristic roots in time-delay systems

Keywords: Characteristic roots, Spectral analysis, Time delay, Stability analysis, Distributed Parameter Systems

This paper studies double imaginary characteristic roots in the case of time-delay systems with two delays as parameters. We aim to identify the direction in which double roots cross the imaginary axis, when the delay parameters change. To determine what happens when parameters are under a small perturbation, we present two methods: the algebraic and geometric approaches. Taking a theoretical example, we show that, even if the two methods are conceptually different, the provided results are consistent.

INTRODUCTION AND PROBLEM STATEMENT

The stability of dynamical systems with time-delays has been a problem of great interest in the last few decades. It has been shown that the presence of delays is responsible for rather complex behaviours, and it might induce stability or instability, function of the domain in which the delays take value. For systems with one delay as parameter, methods of identifying all the stable delay intervals are given in [START_REF] Lee | On the taudecomposition method of stability analysis for retarded dynamical systems[END_REF] and [START_REF] Walton | Direct method for tds stability analysis[END_REF]. As for systems with two delays as parameters, a rich collection of stability charts (the parameter regions where the system is stable) are presented in [START_REF] Stépán | Retarded dynamical systems: stability and characteristic functions[END_REF].

In this paper, we study time-delay systems with two delays, and focus our attention on how the stability charts look like in a neighbourhood of a double (non semisimple) root of the characteristic equation. We present two methods, the algebraic approach in Section 2 and the geometric approach in Section 3, and briefly compare them in Section 4. The algebraic approach we present involves the computation of eigenvectors corresponding to the double root. Geometric approach was introduced by [START_REF] Gu | On stability crossing curves for general systems with two delays[END_REF][START_REF] Gu | Migration of double imaginary characteristic roots under small deviation of two delay parameters[END_REF] and is based on the continuity of characteristic roots as functions of parameters (the two delays). Both methods aim to judge characteristic roots direction of crossing (the imaginary axis) when parameters change.

Consider a linear time-delay system ẋ(t) = A 0 x(t) + A 1 x(t -τ 1 ) + A 2 x(t -τ 2 ),

(1) where x ∈ R n is the state vector, the constant delays τ 1 and τ 2 are real and positive, with A 0 , A 1 , A 2 ∈ R n×n . The characteristic matrix of the system (1) is given by

M (s, τ ) = sI -A 0 -A 1 e -sτ1 -A 2 e -sτ2 , (2)
where τ = (τ 1 τ 2 ) T , I is the n × n identity matrix, and s is the Laplace variable. Similar to the finite-dimensional case, the characteristic equation of system (1) is given by det M (s, τ ) = 0.

(3)

The roots of (3) are called characteristic roots of system (1). Thus, the finite-dimensional nonlinear eigenvalue problem of system (1) can be written as

M (s, τ )u = 0, (4)
where the vector u ∈ C n \ {0} is called a right eigenvector corresponding to the characteristic root s. In a similar manner, we can also construct a left eigenvector v T ∈ C n \ {0} corresponding to the eigenvalue s, satisfying v T M (s, τ ) = 0. Furthermore, the multiplicity of an eigenvalue as a root of the characteristic equation ( 3) is called algebraic multiplicity. An eigenvalue is called simple if its algebraic multiplicity is equal to one. Otherwise, the eigenvalue is called multiple. Note that, since matrices A 0 , A 1 and A 2 are real, eigenvalues are real or appear in complex conjugate pairs. If in the case of a simple eigenvalue there is a single corresponding eigenvector (up to a scaling factor), a multiple eigenvalue can have one or several corresponding eigenvectors. The maximal number of linearly independent eigenvectors corresponding to a multiple eigenvalue is called geometric multiplicity of the eigenvalue. In other words, the geometric multiplicity is equal to the dimension of the null space of the characteristic matrix (2). The geometric multiplicity is less than or equal to the algebraic multiplicity. A multiple eigenvalue is called semi-simple if the algebraic and geometric multiplicities are equal, otherwise the eigenvalue is called non semi-simple. If there is a single eigenvector corresponding to the eigenvalue, then the eigenvalue is called nonderogatory. For instance, the nonderogatory case of zero eigenvalue, s = 0, with algebraic multiplicity two and geometric multiplicity one is known as Bogdanov-Taken singularity, and it has been extensively studied in the literature, especially in models with time delays (see [START_REF] Bogdanov | Versal deformations of a singular point of a vector field on the plane in the case of zero eigenvalues[END_REF], [START_REF] Bogdanov | Bifurcations of a limit cycle for a family of vector fields on the plane[END_REF], [START_REF] Boussaada | Computing the codimension of the singularity at the origin for delay systems in the regular case: A vandermonde-based approach[END_REF], [START_REF] Campbell | Zero singularities of codimension two and three in delay differential equations[END_REF]Yuan (2008), Faria (2003), [START_REF] Takens | Forced ocillations and bifurcations[END_REF]).

In this paper, we suppose that real matrices A 0 , A 1 , A 2 are such that the characteristic function p(s, τ ) := det M (s, τ ) (also called quasi-polynomial) of system (1) has the form p(s, τ ) = p 0 (s) + p 1 (s)e -sτ1 + p 2 (s)e -sτ2 , (5) where p k (s), k ∈ {0, 1, 2} are polynomials of s with real coefficients. Remark 1. For n = 2, consider

A 0 = a 11 a 12 a 21 a 22 , A 1 = b 11 b 12 b 21 b 22 , A 2 = c 11 c 12 c 21 c 22 .
If at least one of the conditions are satisfied, the characteristic function of system (1) can be written of the form (5). Note that conditions (C1)-(C4) are satisfied when matrices A 1 and A 2 are of rank 1. 2 Remark 2. (An equivalent eigenvalue problem). The eigenvalue problem (4) corresponding to the functional differential equation ( 1) is nonlinear and finite-dimensional. By considering (as, for instance, in [START_REF] Hale | Introduction to Functional Differential Equations[END_REF]) a linear operator A = A 0 x(t) + A 1 x(t -τ 1 ) + A 2 x(t -τ 2 ), with the domain

D(A) = φ φ, dφ dθ ∈ C ([-max{τ 1 , τ 2 }, 0] , R n ) , φ(0) = A 0 φ(0) + A 1 φ(-τ 1 ) + A 2 φ(-τ 2 )} , such that Aφ = dφ dθ , where C ([-max{τ 1 , τ 2 }, 0] , R n
) is the space of continuous functions from [-max{τ 1 , τ 2 }, 0] to R n , we can rewrite system (1) as an abstract ordinary differential equation ẋt = Ax t , where

x t (φ) ∈ C ([-max{τ 1 , τ 2 }, 0] , R n ) is the function segment defined by x t (φ)(θ) = x(φ)(t + θ), θ ∈ [-max{τ 1 , τ 2 }, 0]
. Note that the characteristic roots are the eigenvalues of A. This allows us to write an infinite-dimensional and linear eigenvalue problem,equivalent to (4), (sI -A)u = 0, with u ∈ C ([-max{τ 1 , τ 2 }, 0] , C n ). In addition, the corresponding eigenfunctions of eigenvalue s take the form ue sθ , with θ ∈ [-max{τ 1 , τ 2 }, 0]. However, the algebraic approach described in Section 2 is based on the finite-dimensional and nonlinear form of eigenvalue problem (4). 2

This paper studies the behaviour of double imaginary characteristic roots (double roots s = ±iω = 0 of characteristic function ( 5)) under a small deviation of delay parameters (τ 1 and τ 2 ) and compares the geometric approach and the algebraic approach in the case of time-delay systems. By definition (see for instance [START_REF] Campbell | Zero singularities of codimension two and three in delay differential equations[END_REF]; [START_REF] Gu | Migration of double imaginary characteristic roots under small deviation of two delay parameters[END_REF]), we say that system (5) has a double root s 0 = ±iω 0 at ω = ω 0 for τ 0 = (τ 10 τ 20 )

T if and only if

p(s 0 , τ 0 ) = ∂p ∂s s=s0 τ =τ0 = 0 and ∂ 2 p ∂s 2 s=s0 τ =τ0 = 0.
In the sequel, we present the two approaches and illustrate their implications by considering the following example.

Example 3. (Case study) Let the characteristic matrix be cos(1) -1 1 -1 -e -sτ2 + 2e -sτ1 sin(1) s -1 .

M (s, τ ) = s + 2e -sτ1
We compute the characteristic function corresponding to the characteristic matrix (6), and obtain

p(s, τ ) = s 2 -2s + 2 + [(2 cos 1) s -2 (cos 1 + sin 1)] e -τ1s + e -τ2s . (7) For τ = τ 0 = (τ 10 , τ 20 ) T = (1, 2), the characteristic function (7) has a double imaginary root s = s 0 = ±iω 0 at ω 0 = 1. 2 2. ALGEBRAIC APPROACH
The eigenvalue s 0 = ±iω 0 = ±i in Example 3 is double and non semi-simple. This implies that there exist u 0 and v 0 , a right, respectively left eigenvectors, and u 1 and v 1 a right, respectively left generalized eigenvectors, such that conditions

M 0 u 0 = 0 (8) M 0 u 1 + M 0 1 u 0 = 0 (9) v T 0 M 0 = 0 (10) v T 1 M 0 + v T 0 M 0 1 = 0 (11)
are simultaneously satisfied (see [START_REF] Campbell | Zero singularities of codimension two and three in delay differential equations[END_REF], Chapter 7 in [START_REF] Hale | Introduction to Functional Differential Equations[END_REF]), where

M 0 = M (s 0 , τ 0 ), (12)
M 0 1 = ∂M (s, τ ) ∂s s=s0 τ =τ0 . ( 13
) Proposition 4. Eigenvectors u 0 , u 1 , v 0 and v 1 satisfy v T 0 M 0 1 u 0 = 0 (14) v T 1 M 0 1 u 0 = v T 0 M 0 1 u 1 = 0. ( 15
)
Proof. Multiply equation ( 9) by v T 0 on the left and obtain ( 14) in view of (10). Multiply equation ( 9) by v T 1 on the left and obtain

v T 1 M 0 u 1 = -v T 1 M 0 1 u 0 . Multiply equation (11) by u 1 on the right and obtain v T 1 M 0 u 1 = -v T 0 M 0 1 u 1 . Thus, we have (15) because v T 1 M 0 1 = 0 and M 0 1 u 1 = 0. 2
We consider a simple case, where the delay parameters are under a small perturbation ε > 0:

τ 1 (ε) = τ 10 + ετ 11 , (16)
τ 2 (ε) = τ 20 + ετ 21 . (17)
As τ (ε) is smooth, we can write a Taylor series expansion of the characteristic matrix (2). Moreover, provided that the eigenvalue is non semi-simple, we can write a Puiseux series expansion of it and of the corresponding right eigenvector (see, for instance, Kato (1995)):

s(ε) = s 0 + ε 1 2 s 1 + εs 2 + . . . (18) u 0 (ε) = u 0 + ε 1 2 w 1 + . . . (19)
We can now replace equations ( 16), ( 17), (18), and (19) in the eigenvalue problem (4) and write M (s(ε), τ (ε))u(ε) = 0. Note that we use series expansion of exponential functions, as ε → 0 (for instance, we write e -εs0τ10 = 1 -εs 0 τ 11 +

ε 2 s 2 0 τ 2 11 2
). We collect the terms of equal powers of ε, more precisely we collect only the first three orders (ε 0 , ε 1 2 , ε 1 ) and obtain

M 0 u 0 = 0, M 0 w 1 + s 1 M 0 1 u 0 = 0, ( 20
)
s 1 M 0 1 w 1 + s 2 1 2 M 0 2 + s 2 M 0 1 + M 1 1 u 0 = 0, ( 21
)
where the matrices M 0 , M 0 1 , M 0 2 , and M 1 1 are given by

M 0 = s 0 I -A 0 -A 1 e -s0τ10 -A 2 e -s0τ20 , (22)
M 0 1 = I + τ 10 A 1 e -s0τ10 + τ 20 A 2 e -s0τ20 , (23)
M 0 2 = -τ 2 10 A 1 e -s0τ10 -τ 2 20 A 2 e -s0τ20 , (24)
M 1 1 = τ 11 s 0 A 1 e -s0τ10 + τ 21 s 0 A 2 e -s0τ20 .
(25) Notice that equations ( 22) and ( 23) are explicit expressions of ( 12) and ( 13), respectively. For the sake of uniqueness of u 0 and without any loss of generality, we can use the normalization

v T 1 M 0 1 u 0 = 1 = v T 0 M 0 1 u 1 .
(26) Divide equation ( 20) by s 1 = 0 and obtain 1 s1 M 0 w 1 = M 0 u 1 in view of (9). Since M 0 = 0, we can write w 1 as w 1 = s 1 u 1 .

(27) Replace (27) in equality (21) and multiply by v T 0 on the left. We obtain

s 2 1 v T 0 M 0 1 u 1 + s 2 1 2 v T 0 M 0 2 u 0 + s 2 v T 0 M 0 1 u 0 + v T 0 M 1 1 u 0 = 0.
(28) Now we use the second equality from the normalization (26), and condition (14) from Proposition 4 in equation ( 28), and write the expression of s 1 :

s 1 = ± - v T 0 M 1 1 u 0 1 + 1 2 v T 0 M 0 2 u 0 . ( 29
)
Given the continuity properties of the spectrum (Chapter 1 in [START_REF] Michiels | Stability and Stabilization of Time-delay Systems. An Eigenvalue Based Approach[END_REF]) and the Puiseux series expansion (18), it is easy to see that under a small perturbation of delay parameters ( 16)-( 17) the double characteristic root s 0 = iω 0 splits up into two simple characteristic roots, which will move to the right or left halfplane function of the sign of Re(s( )). These results can be summarized in the following proposition and remarks. Proposition 5. Assume that a double and non semi-simple characteristic root of system ( 5) is located on the imaginary axis, but not at the origin. Then, under a small perturbation of delay parameters of the form ( 16)-( 17) this double characteristic root splits up into two simple roots, and each one of them will move towards stability (instability) if Re(s(ε)) < 0 (Re(s(ε)) > 0). 2 Remark 6. Provided that Re(s 0 ) = 0, we might suppose, in general, that the sign of Re(s(ε)) is given by the sign of Re(s 1 ) as defined in equation ( 29), when ε is very small. However, knowing the sign of Re(s 1 ) is not enough to conclude over a global tendency of the double root to move towards stability or instability, as we shall illustrate in the sequel. This is why we take a further step, multiply equation ( 21) with v T 1 on the left, and obtain the following expression of s 2 ,

s 2 = v T 0 M 1 1 u 0 1 + 1 2 v T 0 M 0 2 u 0 v T 1 M 0 1 u 1 + 1 2 v T 1 M 0 2 u 0 -v T 1 M 1 1 u 0 , (30)
in view of normalization (26). Other terms of s(ε) might be found in a similar way. 2 Remark 7. In order to approximate the value of s(ε), function of ε, τ 11 , and τ 21 , we use s(ε) defined as in equation ( 18), where s 1 is given by ( 29), s 2 is given by ( 30), with M 0 1 , M 0 2 , and M 1 1 given by ( 23)-( 25), u 0 and v T 0 are normalized right, respectively left eigenvectors, satisfying (26), and u 1 and v 1 are generalized right, respectively left eigenvectors, satisfying equations ( 8)-( 11). 2 Remark 8. Note that the double non semi-simple root will split up into two simple roots. The plus sign from equation ( 29) corresponds to one of these simple roots, and the minus sign corresponds to the other simple root. 2 Remark 9. Roughly speaking, the double root will split up into two simple roots that will follow one of the two tendencies: either one root moves towards a half-plane, and the other root towards the other half-plane, or both roots moves towards the same half-plane. This means that we have two types of qualitative behaviour when a perturbation ε arises. We illustrate both situations on Example 3. 2

In the sequel, we consider the characteristic matrix (6) in Example 3. A right eigenvector u 0 satisfying (8) is of the form

u 0 = α 1 + iα 2 -(cos(2) -i(-1 + sin(2))) (α 1 + iα 2 ) ,
where α 1 , α 2 ∈ R. We compute the generalized right eigenvector u 1 = (β 1 + iβ 2 u 12 ) T , verifying (9), with β 1 , β 2 ∈ R and

u 12 = (cos(2) -i sin(2))α 1 + (i cos(2) + sin(2))α 2 - -(i + cos(2) -i sin(2)) (β 1 + iβ 2 ) ,
In the same manner we write the left eigenvector v T 0 satisfying equation ( 10)

v 0 = (1 -i) (γ 1 + iγ 2 ) γ 1 + iγ 2 ,
and the left generalized eigenvector v T 1 such that condition (11) holds:

v 1 = -γ 1 -iγ 2 + (1 -i)δ 1 + (1 + i)δ 2 δ 1 + iδ 2 , with γ 1 , γ 2 , δ 1 , δ 2 ∈ R.
The normalization condition (26) leads to the following constraint on α 1 , α 2 , γ 1 , and γ 2 :

α 1 = 0, α 2 = - 1 γ 2 (cos 2 + sin 2 tan 2) , γ 1 = -γ 2 tan 2, γ 2 = 0.
Therefore, the normalized eigenvectors satisfy

u 0 = 1 γ 2 -i cos 2 cos 2(-1 + i cos 2 + sin 2) , u 12 = cos 2(-i cos 2 -sin 2) γ 2 - -(cos 2 -i(-1 + sin 2)) (β 1 + iβ 2 ) , v 0 = γ 2 (1 + i)(1 + i tan 2) -tan 2 + i , v 1 = (-i + tan 2)γ 2 + (1 -i) (δ 1 + iδ 2 ) δ 1 + iδ 2 .
Now we write M 0 2 and M 1 1 of the form

M 0 2 = ∂ 2 ∂s 2 M (s, τ ) s=s0 τ =τ0 , M 1 1 = τ 11 ∂ ∂τ 1 M (s, τ ) s=s0 τ =τ0 + τ 21 ∂ ∂τ 2 M (s, τ ) s=s0 τ =τ0
, and recover the explicit formulae ( 24) and ( 25), respectively. More precisely, we obtain

M 0 2 = 2e -i cos 1 0 -i -(4 -i)e -2i 0 , M 1 1 = -2ie -i cos 1τ 11 0 e -2i iτ 21 --1 + e 2i τ 11 0 .
We compute s 1 function of τ 11 and τ 21 , using equation ( 29):

s 1 = ± (2i + (4 + 2i)e 2i ) τ 11 -2iτ 21 -1 + (1 -2i)e 2i .
Using formula (30), we write s 2 function of τ 11 and τ 21 We notice that if we fix τ 1 and increase τ 2 , then the double root splits up into two simple roots moving towards instability (with positive real part), as depicted in Figure 1, left. On the other hand, if we fix τ 1 and decrease τ 2 , then the two simple roots will move one towards instability, and the other one towards stability, as illustrated in Figure 1, right. developed by [START_REF] Vyhlidal | Mapping based algorithm for large-scale computation of quasi-polynomial zeros[END_REF]. Yellow points s = ±i correspond to ε = 0, τ 1 = 1 and τ 2 = 2. Fix τ 1 = 1 (i.e. τ 11 = 0). For case (a), the characteristic roots of the quasi-polynomial are computed by using QPmR algorithm, for each value of τ 2 , from 2 to 3 using a step of 10 -3 . In Figure 1 left, characteristic roots s(ε) are represented by coloured points (yellow points corresponding to ε = 0). We remark that these points become red as τ 2 increases, i.e. as perturbation ε increases. We can see that double root s = ±i splits up into two simple roots that will move towards the right half-plane. This is coherent to our computation of s(ε) of the form (18), using equations ( 29) and (30). For case (b), we decrease τ 2 from 2 (yellow points) to 1.9 (red points) using a step of 10 -3 , and we notice in Figure 1 right that the double root splits up into two simple roots going in opposite directions, one to the left half-plane, and the other one to the right halfplane. This is also consistent with our prediction based on computing s(ε) as in ( 18).

s 2 = 3 80 + i 80 sec 2(cos 1 -i sin 1) (cos 1 -(1 -i) sin 1)((1 + i) cos 1 + sin 1) [(1 + i)((
Thus, we identify two type of tendencies when a double non semi-simple root lying on the imaginary axis is subject to a perturbation: either one root will go into a half-plane, and the other root into the other half-plane, or both roots move towards the same half-plane.

GEOMETRIC APPROACH

In order to use this approach as described in [START_REF] Gu | Migration of double imaginary characteristic roots under small deviation of two delay parameters[END_REF], the standing assumption

D = det Re ∂p ∂τ 1 Re ∂p ∂τ 2 Im ∂p ∂τ 1 Im ∂p ∂τ 2 s=s0 τ =τ0 = 0 (31)
has to be satisfied. It has been proven that there is a bijection between the complex plane (where the quasipolynomial (5) characteristic roots lie) and parameter space (τ 1 -τ 2 space). More precisely, the parameter space is divided by the stability crossing curve T (defined as in [START_REF] Gu | On stability crossing curves for general systems with two delays[END_REF]) into an S-sector and a G-sector, as depicted in Figure 2. Fig. 2. τ 1 -τ 2 parameter space is divided by the stability crossing curve into a small sector, called S-sector, and a large sector, named G-sector.

The point C in Figure 2 corresponds to the double root s 0 = ±iω 0 on the complex plane. This means that the coordinates of τ 1 and τ 2 corresponding to C are nothing else than τ 10 and τ 20 , for which the quasi-polynomial (3) has a double root at s 0 = ±iω 0 . Moreover, for τ 1 and τ 2 taking values on the curve A C B (on the stability crossing curve T ), the quasi-polynomial (3) has simple characteristic roots s = ±iω. Note that the stability crossing curve contains both positive and negative local stability crossing curves (A C and B C ). For instance, if A C is the positive (negative) local stability crossing curve, and (τ 1 , τ 2 ) takes values on this curve, then the characteristic equation ( 3) has at least one imaginary characteristic root s = iω, with ω > ω 0 (ω < ω 0 ). In addition, the stability crossing curve T is a separation curve that divides τ 1 -τ 2 parameter space into regions, such that the number of characteristic roots on the right half complex plane remain constant as the parameters vary within each such region.

Theorem 7 in [START_REF] Gu | Migration of double imaginary characteristic roots under small deviation of two delay parameters[END_REF] confirms what it has been discussed in the previous section, that when a perturbation occurs on the parameters τ 1 and τ 2 the characteristic roots move according to a specific pattern: either both roots move toward the same half-plane (corresponding to the case where the pair (τ 1 , τ 2 ) moves towards S-sector in the parameter space), or one root moves into the right half-plane and the other one into the left half-plane (corresponding to the case where the parameter set (τ 1 , τ 2 ) moves towards G-sector in the parameter space). However, judging towards which half-plane the two characteristic roots move when (τ 1 , τ 2 ) is in S-sector additionally requires knowing the sign of a parameter κ, defined as follows.

κ = Re ∂ 2 p ∂s 2 - ∂ 3 p ∂s 3 + 3 ∂ 2 p ∂τ 1 ∂s ∂ 2 τ 1 ∂u 2 + 3 ∂ 2 p ∂τ 2 ∂s ∂ 2 τ 2 ∂u 2 s=s0 τ1=τ 10 τ2=τ20 γ=i
, where u represents the perturbation on s. We can evaluate

∂ 2 τ 1 ∂u 2 and ∂ 2 τ 2
∂u 2 , by using

∂ 2 τ 1 ∂u 2 ∂ 2 τ 2 ∂u 2 = - Re ∂p ∂τ 1 Re ∂p ∂τ 2 Im ∂p ∂τ 1 Im ∂p ∂τ 2 -1 Re ∂ 2 p ∂s 2 γ 2
Im ∂ 2 p ∂s 2 γ 2 s=s 0 τ 1 =τ 10 τ 2 =τ 20 , (32) with γ = ±i. Notice that the inverse matrix in equation ( 32) exists due to hypotheses (31). Thus, the case where (τ 1 , τ 2 ) pair moves towards S-sector can be summarized in the following criterion. Criterion 10. If (τ 1 , τ 2 ) is in the S-sector in a sufficiently small neighbourhood of (τ 10 , τ 20 ), then the two characteristic roots in the neighborhood of s 0 are both in the left (right) half-plane if κ < 0 (κ > 0). 2

In the sequel, we use the geometric approach and apply Criterion 10 on Example 3.

To begin with, we verify whether the standing assumption (31) is satisfied, and compute D 1.74159 = 0. The parameter space is depicted in Figure 3. The stability crossing curve has been computed as in [START_REF] Gu | On stability crossing curves for general systems with two delays[END_REF]. We note that the neighbourhood of (τ 10 , τ 20 ) = (1, 2) is divided by the stability crossing curve T into an S-sector and a G-sector. The cusp at (1, 2) corresponds to s 0 = ±iω 0 .

Provided that D > 0, we compute κ 30.7082 > 0. Therefore, according to Criterion 10, both imaginary roots move to the right half-plane as (τ 1 , τ 2 ) moves to S-sector. According to Theorem 7 from [START_REF] Gu | Migration of double imaginary characteristic roots under small deviation of two delay parameters[END_REF], as (τ 1 , τ 2 ) moves to G-sector, one of the two imaginary roots moves to the right half-plane, the other moves to the left half-plane. In other words, as (τ 1 , τ 2 ) moves from the S-sector to the G-sector through the (1, 2), one root moves from the right half-plane to the left half-plane passing through the point i on the imaginary axis, another root on the right half-plane moves to touch the imaginary axis at i then return to the right half-plane.

CONCLUDING DISCUSSION ON THE TWO METHODS

In this paper we have presented two methods that can be used to analyse the behaviour of a double non semisimple characteristic root under a small perturbation of the parameters. More precisely, we have studied the case where the parameters are two delays of a system written of the form (1), for which the characteristic equation is given by (3).

The two methods are rather different and we denote them by algebraic and geometric approach, respectively.

However, both methods suggest that characteristic roots follow the same pattern when a perturbation occurs under the parameters.

The former method is based on eigenvectors computation. Even if in this paper we have only discussed the non semisimple double roots case, a similar study can be made for double semi-simple characteristic roots. The latter method involves the computation of partial derivatives of the quasi-polynomial (3) with respect to the delays, up to the third order. Using the geometric approach implies that condition (31) must be satisfied. Nonetheless, one difference with respect to the method presented in Section 2 is that the geometric approach (and Criterion 10) presented in this paper applies to both semi-simple and non semi-simple characteristic roots. In other words, the main limitation when using the algebraic approach is that we have to separately treat semi-simple and non semi-simple double characteristic roots, and the geometric approach main limitation is that assumption (31) has to be satisfied.

Both approaches show that there are typically two type of situations when a deviation on the delays occurs in the parameter space: either the two imaginary roots move to the same half-plane of the complex plane, or one of the imaginary roots moves to the left half-plane and the other one moves to the right half-plane. This is rather surprising in the case of double non semi-simple characteristic roots: the trajectory of such a root when a perturbation ε occurs is described by the equation ( 18), where s 0 lies on the imaginary axis, and s 1 is given by equation ( 29). Note that the double root splits up into two simple roots, the positive real part s 1 corresponding to a root, and the s 1 with negative real part corresponding to the other simple root. Some publications in the literature suggest that this is why non semi-simple double roots always splits up into two simple roots following only one type of tendency, which is to move towards different half-planes of the complex plane.

The main message of this paper is that we have to pay attention to the quantification of the small neighbourhood, which remains an open problem at least for the algebraic approach. We can briefly illustrate this on Example 3, as follows.

We consider a point in S-sector, (τ 1 , τ 2 ) = (1.1, 2.6), as depicted in Figure 3. We have seen in Section 3 that in this case both imaginary roots move to the right halfplane. We can verify that we obtain the same result by using algebraic approach, by choosing ε = 0.1, τ 11 = 1 and τ 21 = 6 in equations ( 16)-( 17). Thus, we compute s(ε) written of the form (18), yielding the value of the two roots s = 0.975521 -0.254494i and s = 0.826199 -1.58361i. As both roots have positive real part, they are both in the right half-plane. Of course, we can always find an extremely small value for ε such that we end up with two simple roots lying on different sides of the imaginary axis, but this choice rises a natural question regarding the relevance of considering such a value for the control community. We conclude by saying that when choosing values for ε, τ 11 and τ 21 we have to decide upon what a small neighbourhood means, which remains an open question.

(

  C1): {b 11 = 0, b 12 = 0, c 11 = 0, c 12 = 0} (C2): {b 11 = 0, b 21 = 0, c 11 = 0, c 21 = 0} (C3): {b 12 = 0, b 22 = 0, c 12 = 0, c 22 = 0} (C4): {b 21 = 0, b 22 = 0, c 21 = 0, c 22 = 0}

  115 + 82i) cos 1 + (87 + 14i) cos 3 + (61 -2i) cos 5 + (9 + 18i) cos 7-(73 -194i) sin 1 -(53 + 70i) sin 3 -(25 -70i) sin 5 -(45 + 18i) sin 7)τ 11 + 2(cos 1 -i sin 1)((42 -40i) -(21 -4i) cos 2 + (22 -40i) cos 4 -(3 -4i) cos 6 + (4 + 28i) sin 2 + (4 + 4i) sin 4 -(4 -24i) sin 6)τ 21 ]. Now that we have concrete expressions of s 1 and s 2 function of τ 11 and τ 21 for Example 3, we proceed with giving values to τ 11 , τ 21 , and ε in order to illustrate Proposition 5. Suppose γ 2 = β 1 = β 2 = γ 1 = γ 2 = 1. We consider two simple cases, as follows, and compute the value of s(ε) for each case: case τ

FigureFig. 1 .

 1 Fig. 1. Double charcteristic root behaviour for Example 3.Yellow points (lying on the imaginary axis) correspond to τ 10 = 1 and τ 20 = 2. Case (a): fix τ 1 = 1 and increase τ 2 from 2 to 3. Case (b): fix τ 1 = 1 and decrease τ 2 from 2 to 1.9.

Fig. 3 .

 3 Fig. 3. τ 1 -τ 2 parameter space for Example 3.


Thomas Lefeuvre

Emin Koray Kasnakli

Imen Fassi

Isabelle Puaut

Christoph Cullmann

Steven Derrien

Gernot Gebhard

Using Polyhedral Techniques to Tighten WCET Estimates of Optimized Code: A Case Study with Array Contraction

The ARGO H2020 European project aims at developing a Worst-Case Execution Time (WCET)-aware parallelizing compilation toolchain. This toolchain operates on Scilab and XCoS inputs, and targets ScratchPad memory (SPM)-based multi-cores. Data-layout and loop transformations play a key role in this flow as they improve SPM efficiency and reduce the number of accesses to shared main memory. In this paper 1 , we study how these transformations impact WCET estimates of sequential codes. We demonstrate that they can bring significant improvements of WCET estimates (up to 2.7×) provided that the WCET analysis process is guided with automatically generated flow annotations obtained using polyhedral counting techniques.

I. INTRODUCTION

Real-time systems are ubiquitous, and many of them play an important role in our daily life. In hard real-time systems, computing the correct results is not the only requirement. In addition, the results must be produced within predetermined timing constraints, typically deadlines. To obtain strong guarantees on the system temporal behavior, designers must compute upper bounds of the Worst-Case Execution Times (WCET) of the tasks composing the system, in order to finally guarantee that they meet deadlines. Standard static WCET estimation techniques [START_REF] Wilhelm | The worst-case execution-time problem -overview of methods and survey of tools[END_REF] compute such bounds from a static analysis of the machine code. Their goal is to obtain a safe and accurate estimation of a task execution time on a given hardware platform. The safety criterion ensures that the WCET holds for any possible execution of the task on the target platform, whereas accuracy avoids over-provisioning the system.

WCET analysis is confronted with two challenges: (i) extracting knowledge of the execution flow of an application from its machine code, and (ii) modeling the temporal behavior of the target platform. The former issue can benefit from the designer guidance through flow facts (loop bounds, infeasible paths) expressed using annotations. Multi-core platforms make the latter issue even more challenging, as interference caused by concurrent accesses to shared resources have also to be modeled. In this context, the goal of the ARGO H2020 project (http://www.argo-project.eu/) is to develop a WCET-aware parallelizing compilation toolchain, operating from Scilab and XCoS inputs specifications, and targeting scratchpad memory (SPM) based embedded multi-cores.

Obviously, accurate WCET analysis is facilitated by predictable hardware architectures, for which the hardware state during execution can be modeled without introducing too much pessimism in the timing analysis. For example, hardware platforms using ScratchPad Memories (SPMs) instead of caches are considered as more predictable. The reason is that accurately modeling the temporal behavior of data and instruction caches is very challenging, and leads to very pessimistic WCET estimates, for some cache designs (replacement policies, cache hierarchies, write policies, cache coherency).

SPMs are not handled by the hardware, such as caches, but are managed by the software. When SPM management is left to the programmer-managed, it makes them very difficult to use. Many researches have thus studied how to implement automatic SPM management at the compiler level, either as a standalone optimization [START_REF] Deverge | Wcet-directed dynamic scratchpad memory allocation of data[END_REF]- [START_REF] Avissar | An optimal memory allocation scheme for scratch-pad-based embedded systems[END_REF] or in combination with loop and array-layout transformations [START_REF] Baskaran | Automatic data movement and computation mapping for multi-level parallel architectures with explicitly managed memories[END_REF], [START_REF] Darte | Lattice-based memory allocation[END_REF]. For the latter, it has been shown that impressive average-case performance improvements could be obtained on compute intensive kernels, but their ability to reduce WCET estimates remains to be demonstrated. This is the topic we address in this work. More precisely, our contributions are the following:

• We demonstrate the ability of automatic polyhedral optimization techniques to reduce WCET estimates in the case of sequential codes, with a focus on locality improvement and array contraction. We show on representative real-time image processing use cases that they can bring significant improvements of WCET estimates (up to 40%) provided that the WCET analysis process is guided with automatically generated flow annotations. • We present an automatic source-to-source compiler toolchain (part of the ARGO flow) which circumvents some limitations of modern static WCET estimation tools in presence of aggressive code transformations. The approach helps WCET estimation tools with automaticallygenerated flow annotations.

The remainder of this paper is organized as follows. Section II provides background information on the ARGO project and recalls how polyhedral compilation techniques can be used to take advantage of SPMs to improve performance. Section III describes our approach based on automatic derivation of flow annotations, and Section IV provides preliminary results on two representative use cases (real-time image processing). Conclusion and future work are sketched in Section V.

II. BACKGROUND

In this Section, we first provide an overview of the ARGO project. Then, we describe the loop and array-layout transformations considered in this work.

A. The ARGO project

The goal of the ARGO project is to design a WCET-oriented parallelizing tool-chain for embedded multi-core architectures. The ARGO flow [START_REF] Derrien | Wcet-aware parallelization of modelbased applications for multi-cores: The ARGO approach[END_REF], operates on user applications expressed either as XCoS models (Simulink-like) or in the Scilab language. This input specification is then parallelized and optimized for WCET on an embedded multi-core architecture.

A fundamental requirement for obtaining WCET estimates using WCET analyzers, such as AbsInt's aiT [START_REF] Ferdinand | aiT: worst case execution time prediction by static program analysis[END_REF], is that the timing behavior of all processor cores is predictable. For example, predictability at the core level is compromised by speculative hardware mechanisms such as caches or branch prediction. On multi-core platforms, the additional problem of shared hardware resources for which the interference due to concurrent accesses have to be considered, arises. To keep the problem tractable, we only consider in ARGO hardware platforms with the following characteristics for each core:

• A multi-core platform composed of time-predictable processors (here Leon3 or Xentium processors), supported by tools such as AbsInt's aiT [START_REF] Ferdinand | aiT: worst case execution time prediction by static program analysis[END_REF] • Memory hierarchies based on private SPMs instead of data and instruction caches, in order to avoid challenging cache hit/miss analysis In this paper, we will restrict ourselves to one core of the InvasIC platform, a predictable configurable and tile-based multi-core platform based on Leon3 cores.

One of the challenges for ARGO is to be able to take advantage of the use SPMs as a replacement for data and instruction caches. Although SPM-based platforms are much easier to model from a WCET point of view, they raise many challenges from the compiler point of view, as scratchpad must be explicitly managed by the software. Because of their widespread use in embedded systems, there has been a lot of research on compiler support for dynamic scratchpad management [START_REF] Deverge | Wcet-directed dynamic scratchpad memory allocation of data[END_REF]- [START_REF] Avissar | An optimal memory allocation scheme for scratch-pad-based embedded systems[END_REF]. Most of this work aim at determining (i) which variable should be allocated to scratchpad (ii) when a given variable should be moved in scratchpad and copied back to shared memory. The problem is that in many cases (e.g real-time image processing), the arrays manipulated by the program are way too large to fit in a scratchpad memory.

It is however possible to circumvent this issue by resorting to array layout transformations. Such transformations, explored in this paper, can be used to reduce the array sizes (array contraction) and/or partition arrays into regions so that they can fit the scratchpad.

B. Polyhedral based SPM management

Polyhedral based loop and data-layout transformation techniques are very efficient when it comes to optimizing memory accesses for SPM-based platforms [START_REF] Baskaran | Automatic data movement and computation mapping for multi-level parallel architectures with explicitly managed memories[END_REF]. For example, they offer very efficient array contraction algorithms [START_REF] Darte | Lattice-based memory allocation[END_REF].

Such algorithms consists in reducing the footprint of all temporary/intermediate arrays. This is achieved through complex re-indexing transformations involving circular (i.e. modulo based) addressing functions. The goal is generally to find the re-indexing transformation leading to the smallest possible array. The technique is based on array cell live range analysis, which is made possible thanks to the polyhedral representation. The result of this analysis is then used to identify conflicting array cells locations (cells that are alive at the same time). From there, finding a legal contraction consists in obtaining a re-indexing function that does not remap two conflicting cells to the same memory location.

When used alone, array contraction brings few benefits, but, its efficiency can be significantly improved when combined with locality improving loop transformations (tiling and/or fusion). Thanks to steady improvements in the polyhedral compilation techniques, nowadays, we have algorithms that automatically find such program transformations [START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF].

Since a detailed discussion on array contraction is out of the scope of this work, we simply illustrate the effect of the transformation on an example, given in Figure 1. As visible in the figure, and further detailed in Section III, the optimized code now contains complex control flow, which may challenge WCET analysis tools. To address this issue, we enrich transformed code with annotations, that we describe in the following Section.

III. TIGHTENING WCET ESTIMATES USING AUTOMATICALLY GENERATED FLOW ANNOTATIONS

This Section is organized as follows. We first discuss the scalablity issue with static WCET estimation tools for complex loop nests, and then describe our approach based on source level flow annotations.

A. Scalability of static WCET analysis

As mentioned in the previous Section, polyhedral transformations often lead to complex loop structures (with guarded statements and non constant loop bounds for example). For average-case performance, this overhead is largely outweighed by others benefits (such as parallelism and locality). This is not true when dealing with worst-case performance, as this complexity can result in very pessimistic WCET estimates, or no WCET estimate at all when the estimation tool cannot derive loop bounds.

This complexity is illustrated in Figure 1, which shows a code excerpt taken from the optimized PIPS use case (see Section IV for more details). The interested reader may notice that, in this code snippet, loop L3 is guarded by a predicate expression based of the outer loop index c1. Similarly, statements S1 and S2 are guarded by complex combinations of guards. Because it cannot always determine precisely how often a guard will be valid, a static WCET tool may end-up overestimating the WCET by considering the worst-case path for every loop iteration. This will then result in a pessimistic WCET estimate.

As mentioned earlier, the ARGO project relies on the AbsInt aiT tool to perform binary level WCET analysis. As for all static WCET estimation methods [START_REF] Yau-Tsun | Performance analysis of embedded software using implicit path enumeration[END_REF], the tool must be able to determine an upper bound of the number of iterations for every loop to be able to estimate the WCET. Because inferring loop bounds from binary is challenging, aiT offers two ways to derive as precise as possible loop bounds.

• The first approach relies on annotations, which are used to specify loop bounds that could not otherwise be inferred by the tool. Such annotations have to be provided by the user/programmer, and because of this are errorprone and may lead to unsafe WCETs (due to incorrect annotations in the case of too complex source code for example). • The second one is based on virtual unrolling, in which several instances of a same loop/functions are considered during the analysis (each instance corresponding to a specific context). This approach is very efficient for nested loops operating on small/medium data-sets. In most cases, the use of aggressive virtual unrolling leads to very accurate WCET estimates, as the context associated to each loop instance can be used to eliminate dead guards in the execution path. As a matter of fact, in our case, the choice of a large enough unrolling factor would always lead to a tight WCET estimate, as we end-up performing a complete symbolic execution of the loop nest. Unfortunately, none of these two approaches can be used in the context of the ARGO project.

Since ARGO operates from very high-level models (Scilab, XCos), an approach based on user-provided loop bound annotations is not practical, because the objective is to hide details of code transformations to the end-user.

In addition, we also experimentally observed that virtual unrolling does not scale well for very large loop trip counts and large unrolling factors. In some cases, (see Section IV for details) the analysis had to be stopped after two hours, or failed due to memory exhaustion, which is not compatible with the user-in-the-loop approach as advocated by the ARGO toolset.

To address these shortcomings, we propose a meet-inthe-middle approach, which combines polyhedral-based static analysis (implemented at the C level in the Gecos sourceto-source compiler [START_REF] Floch | Gecos: A framework for prototyping custom hardware design flows[END_REF]), with flow annotations (used at the binary/assembly level in AbsInt aiT WCET analysis tool).

B. Polyhedral based flow analysis

Figure 2 summarizes the proposed loop and array transformation framework. Source-to-source transformations are first applied using a polyhedral representation of loop nests. The polyhedral transformations achieved in this Figure are loop fusion combined with loop shifting; the two for-loop nests are fused and the inner-most loop (j-loop) is shifted by a factor equal to 2 in order to respect the dependencies. Flow annotations are then generated and inserted in the source code using aiT flow annotation facilities.

Flow annotations are associated to basic blocks and can be used, among others, to specify the exact execution count of a given block along the worst-case execution path. As currently implemented in aiT, the annotation follows the syntax described below, where labelX is the address of a machine instruction (identified through a label) and cnt is the computed execution count. Similar annotations are provided to specify maximal execution counts of basic blocks.

ais2 {flow sum: point("labelX") == cnt; }.

In a polyhedral representation of a program, a domain is associated to every statement in the loop nest. This domain is defined using a set of affine constraints over loop indices and parameters. The dimension of the domain (i.e. the polyhedron) corresponds to the depth of the loop-nest it represents. In our use cases, we did not observe any scalability problem. Since the code regenerated after optimization also admits a polyhedral representation, we can reconstruct the domains of all the statements in the regenerated code, as we do from the original program. Because the program transformation is not bijective statement wise (a statement in the original program may be transformed into more than one statement because of loop peeling or loop unrolling for example) we cannot operate on the initial polyhedral representation, but need to reconstruct the domains of the transformed statements.

Deriving the flow information associated to a statement then amounts to counting the number of iterations in its corresponding domain. Counting the number of integral points in a polytope is a well studied problem [START_REF] Clauss | Counting solutions to linear and nonlinear constraints through ehrhart polynomials: Applications to analyze and transform scientific programs[END_REF], [START_REF] Verdoolaege | Counting integer points in parametric polytopes using barvinok's rational functions[END_REF], and most state-of-the-art polyhedral compilation toolboxes offer this ability. In the general case, the outcome of the counting process takes the form of a piecewise polynomial expression. In our case, since we currently only deal with constant sized domains, the result is always a constant.

For example, consider the basic-block containing the statement S1 in the transformed code example from Figure 2. After extracting the polyhedral representation of the transformed program, we obtain its domain D S1 defined by the following equalities:

D S1 = {i, j|i ≥ j ∧ 0 ≤ i < 640 ∧ 0 ≤ j < 480}
The flow annotation that we are interested in corresponds to the number of integral points in its associated domain (that we write card(D S1 )), in this case card(D S1 ) = 153280. The derived flow information remains attached to a source level statement (and its corresponding basic-block in the Gecos compiler intermediate representation), not to a machine instruction label as required by aiT. To bridge the gap between the two tools, we insert custom labels using inlined asm statements, as depicted in Figure 2, and make sure the resulting program is compiled without optimization which could make the annotation invalid, in case the C compiler applies loop transformations.

IV. EXPERIMENTS

As mentioned in the introduction, the goal of this work is to tighten WCET estimates through the combined use of SPMs, locality enhancing loop transformations and data layout optimizations. Our goal is to demonstrate that, a meetin-middle approach involving both source-level and binarylevel analysis is mandatory to observe WCET performance improvements for complex compute-intensive loop kernels.

A. Benchmarks and experimental setup

In our experiments, we considered two real-time image processing applications:

• The first one is a classical image processing pipeline which implements the classical Harris Corner detection algorithm. The implementation operates on high definition images, and consists of 4 stages. The implementation uses 10 distinct arrays of 4 MB each.

• The second application is used for post-processing raw data from a polarized image sensor, used in factory automation for product defect detection. This application, named PIPS hereafter, is an industrial use case from the ARGO project. Given its size and relative complexity, it was deemed relevant as a real world example code for our experiments. The application consists of several processing stages, each of them operating on its own intermediate version of the image. The implementation uses 15 distinct statically allocated arrays for storing intermediate results, with sizes ranging from 2.46 MB to 9.83 MB. In this work, we use the Leon3 core as target processor, as it is one of the building blocks for the InvasIC multi-core architectures designed in the ARGO project. The Leon3 core is an open-source implementation of the Sparc-V8 instruction set by Gaisler research (http://www.gaisler.com/). In our experiments, the processor does not include cache memory, but instead uses a 64 KB on-chip SPM. The compiler used for the experiments is gcc 4.9.4 for Sparc. In the following, we always consider the code to be executed from the SPM. We configured aiT for the InvasIC core, by specifying two address ranges: one for the SPM, and one for the external memory, with distinct access times. We considered 0 cycle delay penalty for SPM access (i.e a load/store instruction takes 1 cycle to execute) and used several distinct values (16, 32 and 64 cycles) for the access delay to shared main memory. The informed reader may observe that these values are higher than typical DRAM access latencies. However, although this paper concentrates on sequential codes, we are interested in the scope of ARGO in multi-cores, where tens of cores may have interfering accesses to this memory through its controller, explaining our choice. The same configuration for aiT is used in all experiments reported in this Section, except if the opposite is explicitly stated in the text.

B. Impact of compiler optimizations on WCET estimates

As a first experiment, we use aiT to estimate WCETs for the original version of our target algorithms under different levels of compiler optimizations (from -O0, no optimization, to -O3). The goal of this experiment is to observe the impact of optimizations on the tool ability to derive WCET estimates. We performed this experiment with various levels of virtual enrolling (from 2, default value, to 256 in order to test for scalability) to help the tool determining flow information. The results of this experiment are provided in Table I for both Harris and PIPS. Shared memory access time for this experiment is set to 32, and without array contraction, manipulated arrays are too large to fit in SPM, that is then used to store the code only.

Our results show that the ability of aiT to derive WCET estimates is highly dependent on the program control flow complexity. For example, aiT is able to derive bounds for all versions of the original Harris kernel. On the PIPS use case, the tool was unable to derive a WCET estimate for PIPS when optimizations are enabled, even with the help of virtual unrolling. As depicted in Table I

⊥ -O1 ∞ ∞ ∞ ∞ ⊥ -O2 ∞ ∞ ∞ ∞ ⊥ -O3 ∞ ∞ ∞ ∞ ⊥ TABLE I WCET ESTIMATES FOR ORIGINAL CODE (IN GCYCLES).
could not infer some loop bounds (marked as ∞) or causes memory exhaustion (marked as ⊥) after two hours of runtime.

C. Impact of polyhedral transformations

As a second experiment, we have used aiT to estimate the WCET on the transformed version of our examples, after loop fusion and array contraction is applied. To avoid costly modulos operation, we increased the size of all circular arrays so as to reach a power of two (modulos are then translated into simple bitmasking operations). SPM access access delay is set 0 cycle, whereas shared memory access time is set to 32 cycles.

The results of the transformation for our two use cases are summarized below:

• For the Harris use case, the size of all 10 intermediate arrays was reduced thanks to the combination of fusion and contraction (4 arrays reduced to scalars, 5 arrays sizes reduced by more than 98%). All intermediate arrays could then fit in the SPM. • For the PIPS use case, 6 arrays out of 15 benefited from our transformation (4 were reduced to scalars, one was reduced from 2.46 MB to 11.52 KB, and one from 7.3 MB to 24 Bytes. Here again, all intermediate arrays could fit in the SPM. We first exercised the aiT ability to estimate a WCET using its virtual unrolling feature enabled, and this without flow annotations. The results of this experiment are provided in table II for Harris and PIPS, for an external memory access latency of 32 cycles. In the following Tables, UF stands for aiT Unrolling Factor.

They show that for Harris, aiT is still able to derive a WCET estimate without annotations, for all optimization levels except -O1 (at the time of writing, we have no explanation for this observation). For PIPS, the tool can only find a WCET estimate at -O0 level. For higher optimization levels, the analysis fails because of some missing loop bounds or because of memory exhaustion. Virtual unrolling significantly tightens WCET estimates (by a factor 66 for PIPS). However, this tighter WCET estimate comes at the cost of increased analysis time (for PIPS compiled at -O0, the analysis time increased from 2 s for UF=2 to 23 mn for UF=64, still running after 2 hours for UF=256).

As level. A key issue is that flow annotations are only valid if the compiler does not alter the control flow exposed in the original source code. Although we observed that the use of assembly labels in the source prevents many gcc optimizations, and somewhat preserves the validity of the flow annotations, we currently have no systematic way of ensuring they remain valid other than compiling the code with the -O0 flag. Resorting to -O0 is probably too conservative, as there exist many optimizations that are flow-safe in the sense that they would not impact the validity of the source level flow annotations at assembly level [START_REF] Li | Traceability of flow information: Reconciling compiler optimizations and WCET estimation[END_REF]. However, identifying such flow-safe optimizations remains an open and challenging problem, especially given the complexity of current compiler frameworks (for example gcc has more than 150 user-exposed optimization flags). The results of the WCET analysis using flow annotations and compiled with optimization level -O0 is provided in Table III for both Harris and PIPS. They show that we are able to obtain WCET estimates in both cases, although the results for Harris are not as good as those obtained by aiT alone.

⊥ -O1 ∞ ∞ ∞ ∞ ⊥ -O2 ∞ ∞ ∞ ∞ ⊥ -O3 ∞ ∞ ∞ ∞ ⊥

D. Overall benefits and discussion

The question raised by this work is the following: are polyhedral compiler optimizations useful for reducing WCET estimates?

To provide a fair answer to this question, we need to compare the best WCET estimate for the original program with the best WCET estimate for the transformed program. In our context, the best WCET estimate corresponds to the lowest safe2 WCET estimate obtained using the various compiler optimization levels available (-O0, -O1, -O2 and -O3 the case of Harris, we compare the WCET estimate obtained for the original code at -O1 (surprisingly slightly lower than at -O2 and -O3) with the WCET estimate obtained for the transformed code at -O2 (slightly lower than at -O3). In the case of PIPS, we are only able to estimate the WCET at optimization level -O0 for both the original and transformed versions, and therefore compare these two results.

The WCET improvements observed for these two application are provided in Table IV and show improvement ratios between 1.2 and 2.7. For each example, we show the performance improvements for a given access delay to the shared external memory. As expected, for larger memory access times, the improvement offered by array contraction increases.

Of course, there are certainly examples where the locality enhancing polyhedral optimizations would end-up to be counter productive from a WCET point of view. A typical scenario would be a program for which aiT could automatically derive a WCET at -O3 before the polyhedral transformation, but would fail after transformation for all optimization levels but -O0 (thanks to flow annotations). However, we did not witness such a case so far.

Besides, in the context of the ARGO flow, the use of polyhedral loop transformations is not restricted to improve the efficiency of SPMs, it is also used to expose task level parallelism that can be later used by the scheduler/mapper to map and schedule tasks on cores. We therefore expect these transformations to also contribute in improving systemlevel WCET performance (i.e. the WCET of the parallel application).

V. CONCLUSION

The ARGO H2020 European project aims at designing a WCET-aware parallelizing compilation flow for predictable embedded multi-cores. Among the optimizations involved in such a flow, data-layout optimizations and parallelizing loop transformations obviously play a key role. In this paper, we study how such transformations effectively impact WCET analysis. We demonstrate on two use cases that polyhedral transformation can bring significant improvements of WCET estimates (up to 2.7 × provided they can also drive the WCET analysis process through flow annotations, that can be computed automatically using well known polyhedral counting techniques.

Another outcome of this work is the observation that understanding the impact of combinations of complex compiler optimizations on WCET is a very challenging problem. This complexity is also acknowledged in the context of averagecase performance, and has motivated research work on iterative compilation, which consists in exploring many combination of optimizations to pick-up the combination offering the best performance. We believe that a similar approach for WCET would be very relevant, especially in the context of complex toolchains such as ARGO, where both the compiler and WCET analyzers are used as grey boxes.

Fig. 1 .

 1 Fig. 1. Excerpt of the PIPS kernel after polyhedral transformation.

Fig. 2 .

 2 Fig. 2. Summary of the proposed WCET aware loop and array transformation framework.

  for the PIPS example, aiT either terminates, but fails to estimate the WCET because it

	kernel	opt	2	Virtual unrolling factor 4 16 64	256
	Harris	-O0 2.512 2.507 2.507 2.507 -O1 1.073 1.072 1.072 1.072 -O2 1.074 1.073 1.073 1.073	⊥ ⊥ ⊥
		-O3 1.173 1.173 1.173 1.173	⊥
		-O0 20.54 20.54 20.53 20.52
	PIPS

  explained earlier, we propose to guide the WCET analyzer by computing flow annotations at the source code

	kernel	opt	2	Virtual unrolling factor 4 16 64	256
	Harris	-O0 -O1 -O2	2.086 ∞ 0.563	2.078 ∞ 0.561	2.078 2.078 ∞ ∞ 0.561 0.561	⊥ ⊥ ⊥
		-O3	0.566	0.564	0.564 0.564	⊥
		-O0 1548.48 1545.53 26.62 23.33
	PIPS

The work presented in this paper is part of ARGO (http://www.argoproject.eu/), funded by the European Commission under

Horizon 2020 Research and Innovation Action, Grant Agreement Number 688131.

This excludes WCET estimates at -O3 needing additional flow annotations, because compiler optimizations may make the annotations invalid.


Y Kebukawa
email: kebukawa@ynu.ac.jp

M Ito

J Aoki

T Okada

Y Kawai

J Matsumoto

R Nakamura

H Yano

K Terada

M Toyoda

Mission Y Kebukawa

H Yabuta

N Grand

H Cottin

A Buch

C Briois

L Thirkell

O Mori

J Ka- Waguchi

IN-SITU ANALYSIS OF SURFACE AND SUBSURFACE SAMPLES FROM A JUPITER TROJAN ASTEROID USING A HIGH RESOLUTION MASS SPECTROMETER IN THE SOLAR POWER SAIL MISSION

d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction:

The Solar Power Sail (SPS) mission is one of candidates of the upcoming strategic middle-class space exploration to demonstrate the first outer Solar System journey of Japan. The mission concept includes unique instruments of in-situ high resolution mass spectrometry (HRMS) and possible sample return capability as well as remote sensing instruments. The current mission sequence proposes the launch in late 2020s, and rendezvous to a D or P type Trojan asteroid of ~20-30 km in diameter in 2030s after Jupiter flyby. The overview of this mission is presented by Okada et al. in this meeting [1].

Here we present specific scientific goals of in-situ analysis with a HRMS and sample return from the surface and subsurface (up to 1 m) materials of a Jupiter Trojan asteroid.

Scientific Goals:

The key questions for the Jupiter Trojan asteroid exploration are: (1) constraining planet formation/migration theories, (2) evolution and distribution of volatiles (water and organics) in the Solar System, (3) origin of Earth's water, and (4) surface processes of Jupiter Trojan asteroids.

(1) Constraining planet formation/migration theories. The classic model for the Solar System formation suggests that the Trojan asteroids are mainly survivors of building blocks of the Jupiter system [2], while recently proposed planetary migration models (e.g., Nice model) claim that they are intruders from outer regions after the planetary migration of the giant gas planets settled [3]. Thus, origin of Jupiter Trojans potentially contains a key to understand the planetary formation and migration. We can roughly consider that in the case of the former scenario, the Trojans have similar composition to primitive main belt asteroids (e.g., C, D, or P type), and the later scenario, the Trojans have similar composition to comets (Kuiper belt objects, KBOs). The apparent differences between primitive asteroids and comets are, in general, comets have not aqueously altered, and enriched in heavy isotopes (D and 15 N) than primitive asteroids due to isotopic fractionations in cold environments (e.g., molecular clouds and the outer protosolar disk) [4].

(2) Evolution and distribution of volatiles in the Solar System. It is first opportunity to in-situ investiga-tion and sample return of D/P type asteroid, and it could be a Missing link between comets and asteroids. Since the only known D/P type asteroid sample is the Tagish Lake meteorite (an ungrouped carbonaceous chondrite), our knowledge of these types of asteroids is very limited by observational studies. Consistent with aqueously-altered nature of the Tagish Lake, main belt D type asteroids show 3 μm bands which indicates phyllosilicate OH, but Jupiter Trojans do not show the 3 μm band [5]. Thus, we might expect the very early stage of aqueous alteration and/or water ice (at subsurface) from the Trojan samples.

What was origin and evolutional processes of extraterrestrial organic matter? Comets contain larger amount of volatile simple molecules than primitive asteroids, i.e., carbonaceous chondrites. Alexander et al. [6] hypothesized that insoluble organic matter (IOM) in various chondrites and possibly cometary refractory organics evolved from a common precursor, which has high H/C ratios and high δD. We might be able to confirm that theory with the analysis of a Trojan asteroid.

Finally, the stable isotopic compositions of H, C, N and O for volatiles in the Trojans give us an important insight into circulation, distributions and evolution of gas and solid materials within the Solar System.

(3) Origin of Earth's water. D/H ratio of water in the Solar System objects is interesting, with regards of origin of Earth's water. Oort cloud comets have higher δD than Earth's water, a Jupiter family comet (JFC) 103/P Hartley 2 shows similar δD to the Earth's water but recent result from Rosetta mission shows that water in 67P/Churyumov-Gerasimenko (a JFC) has much high δD [7]. While, water in chondrites are estimated to lower than the Earth's water [8]. How about the Trojan asteroid?

(4) Surface processes of Jupiter Trojan asteroids. Samplings from both surface and subsurface are good opportunity to study space weathering and surface evolution of the Jupiter Trojan asteroid. Noble gas isotopes are sensitive indicators to elucidate a history of irradiation from cosmic rays and/or solar wind on its surface, as shown by the Itokawa particles by the Hayabusa mission [9].

In-situ analytical sequences:

We plan to analyze volatile materials on the Jupitar Trojan, for their isotopic and elemental compositions using a HRMS with a combination of pyrolysis ovens and gas chromatography (GC) columns. This HRMS system allows to measure H, N, C, O isotopic compositions and elemental compositions of molecules prepared by various pre-MS procedures including stepwise heating up to 600ºC, pyrolysis-GC, and high-temperature pyrolysis with catalyst in order to decompose the samples into simple gaseous molecules (e.g., H 2 , CO, and N 2 ). The required mass resolution should be at least 30,000 for analyzing isotopic ratios (e.g., H 2 16 O, HD 16 O and H 2 18 O for H and O isotopic measurements) for simple gaseous molecules. For elemental compositions of molecules/ions, mass accuracy of ~10 ppm is required to determine elemental compositions for molecules with m/z up to 300 (as well as compound specific isotopic compositions for smaller molecules). Our planned analytical sequences consist of three runs for both surface and subsurface samples (Fig. 1). In addition, 'sniff mode' which simply introduces environmental gaseous molecules into a HRMS will be done by the system. The details of the analytical methods and apparatus are under developments.

Sampling Flash Heating 300 & 600ºC

Pyrolysis ~1000ºC

Low molecular weight gas CO, CO 2 , H 2 O, CH 4 , N 2 etc.

HRMS

Step heating Up to 600ºC

GC

Isotopes: D/H, 15 N/ 14 N, 13 The sample return from the Trojan asteroid: Analyses of returned samples from Moon [START_REF] Jolliff | New views of the MOON[END_REF], asteroid [11] and comet [12] were essential to understand their origin and nature as well as increasing our knowledge about the Solar System. The most recent returned sample was from the S-type asteroid Itokawa by Hayabusa mission in 2010. The results by series of analyses provided new insights for the connection to meteorite researches, space weathering processes, small asteroidal body formation in the Solar System [e.g. , 11]. JAXA Hayabusa 2 and NASA Osiris-REx are both current sample return missions from the organic-rich asteroids, Ryugu (C-type) and Bennu (B-type), respectively [13,14]. Both missions have complementary scientific goals that are to understand the Solar System evolution in the point of view of organics, water and associated minerals. We, therefore, are working on the possibility of the sample return from Trojan asteroid that is expected to contain primordial chemical information at the very beginning of Solar System formation. D/P-type Jupiter Trojan asteroids likely consist of dominant of organics (carbonaceous materials) and anhydrous silicates (hydrated silicates cannot be excluded), possibly with water (ice) in its interiors [START_REF] Guilbert-Lepoutre | Lunar and Planetary Science XLVIII[END_REF]. Beside in-situ HRMS analysis of isotopic ratios, elements and molecules in surface and subsurface samples on the Trojan asteroid, analysis of returned samples containing non-volatile materials (organics and minerals) as well as water (ice) will open a new insight of the detailed scientific objectives for the Solar System evolution. Since, in-situ analysis is limited in terms of sample preparations, lack of relationship among components, and mineralogical/petrological contexts, the state-of-the-art microanalysis techniques on the Earth will provide these additional information such as isotopic ratios of individual component (organics and associated minerals), trace amount of gaseous species (e.g., Noble gases, CO, CO 2 , NH 3 , CH 4 in the ice), and organic compounds that are hard to be detected under the current in-situ HRMS system (e.g., amino acids).

The details of the sample return capsule is not yet fixed but a cryo-system is highly encouraged. Thus, we will receive "extraterrestrial ice (water)" that has a pristine water at the Solar System which contains the information of nebular gas, formation of ice, reservoir of volatiles (water and organics), and the origin of the Earth's water.

Fig 1 :

 1 Fig 1: In-situ analytical sequences of surface and subsurface samples of a Jupiter Trojan asteroid using a high resolution mass spectrometer (HRMS).

2 , CH 4 , CO, N 2

	Organics, Water,
	Hydrous minerals
	+ Sniff mode
		Isotopes Molecules
		Volatile organics and
	Preheating	fragments of refractory organics Molecules
	300/600ºC
		Isotopes
	with nickel catalyst

C/

12

C, 18 O/ 16 O Molecules: Elemental compositions of gases molecules and fragments from refractory organics Refractory OM  Gases: H


César Marcon

Ney Calazans

Fernando Moraes

Altamiro Susin
email: susin@inf.ufrgs.br

Igor Reis

Fabiano Hessel
email: hessel@inf.pucrs.br

Exploring NoC Mapping Strategies: An Energy and Timing Aware Technique

Complex applications implemented as Systems on Chip (SoCs) demand extensive use of system level modeling and validation. Their implementation gathers a large number of complex IP cores and advanced interconnection schemes, such as hierarchical bus architectures or networks on chip (NoCs). Modeling applications involves capturing its computation and communication characteristics. Previously proposed communication weighted models (CWM) consider only the application communication aspects. This work proposes a communication dependence and computation model (CDCM) that can simultaneously consider both aspects of an application. It presents a solution to the problem of mapping applications on regular

NoCs while considering execution time and energy consumption. The use of CDCM is shown to provide estimated average reductions of 40% in execution time, and 20% in energy consumption, for current technologies.

Introduction

SoC design requires special communication resources to fulfill stringent design requirements. Deep sub-micron effects pose formidable physical design challenges for long wires and global on-chip communication. Many designers have proposed a change from the mainstream synchronous design paradigm to a globally asynchronous, locally synchronous (GALS) paradigm [1]. In a GALS design, the application is partitioned into synchronous domains. Each domain is locally synchronous and placed inside a limited region, usually called tile. An asynchronous communication resource provides the link between synchronous domains. A NoC is an intra-chip communication infrastructure, usually composed by a set of routers interconnected by point to point communication channels, implementing a chosen topology. The NoC channels can be designed to provide an asynchronous communication protocol between otherwise synchronous domains. NoCs can then be easily adapted to implement systems based on the GALS paradigm. Besides, they present high scalability, reusability, and reliability [2].

Consider a SoC implemented using the GALS para-digm and composed by n cores. Suppose this SoC employs a NoC as its internal sole communication resource.

The application mapping problem for this architecture consists in finding an association of each core to a tile (a mapping) such that some cost function is minimized.

In the most general case, the problem allows n! possible solutions. Given a future SoC with hundreds of tiles [3], exhaustive search of the problem solutions space will rapidly become unfeasible. Consequently, the optimal implementation of such SoCs requires efficient mapping strategies. Some mapping strategies have been proposed. For example, [4] and [5] propose a communication weighted model (CWM) to account for the overall communication volume of each channel. However, CWM does not consider communication timing. This paper proposes a communication dependence and computation model (CDCM) to capture both, the volume and timing of application communication. Comparing these models for a 0.07 technology, CDCM produced estimated average reductions of 40% and 20% in execution time and energy consumption w.r.t. CWM.

The remaining of the paper is organized as follows. Section 2 discusses related work. Section 3 gives a formulation of the mapping problem. Section 4 describes and compares CWM and CDCM algorithms and Section 5 presents experimental results. Section 6 presents some conclusions and directions for further work.

Related Work

Hu and Marculescu [4] showed that by using mapping algorithms it is possible to reduce by more than 60% the energy consumption when compared to random mapping solutions. The authors propose the use of an application characterization graph (APCG), a way of capturing CWMs. Murali and De Micheli [5], propose a solution similar to that in [4]. Their CWM is represented by a structure called core graph. The main goal of their work is to propose an algorithm that maps cores on mesh NoC architectures under bandwidth constraints, trying to minimize average communication delay.

Ye et al. [6], propose a model to evaluate the energy consumption in a communication infrastructure containing routers, internal buffers and interconnected wires. The same authors [7] describe the contention problem in NoCs, evaluating the associated performance reduction. They propose as solution a routing algorithm that minimizes energy consumption, by reducing the required buffers in the communication network.

The approach of the present paper proposes and uses CDCM, which models application packet dependence and computation time. It explores CDCM and CWM strategies to solve mapping problems fulfilling energy consumption and execution time requirements. It stresses that CWMs as the ones presented in [4] [5] abstract communication timing, an essential information to estimate execution time and energy consumption of the application. CDCMs lead to a better mapping solution if compared to CWMs, with low extra computational effort. Moreover, the dynamic energy model presented in [6] is extended here to include static energy consumption. Instead of evaluating different routing algorithms, as in [7], this work explores communication dependence graphs to reduce communication buffers, saving area, execution time and energy.

Problem Formulation

Each core behavior can be modeled by its computation and communication characteristics. Here a formulation of the mapping problem is proposed using graph structures, together with models, strategies, and algorithms for it.

Graph Definitions

Graph structures used in this work are defined here: CWG, CDCG, and CRG. Figure 1 W represents all communications between application cores, while CWG reveals information of application relative communication volume. This is similar to the definitions of APCG from [4] and core graph from [5].

Definition 2: Let C be a set of cores of a given application. The communication dependence and computation graph or CDCG of this application is a directed graph <P, D>. The set of vertices P contains all packets exchanged between every pair of cores communicating in an application. There are also two special vertices named Start and End. The set of edges D contains all communications dependences in an application. Elements of P are 4-tuples in the form p abq = (c a , c b , t aq , w abq ), where c a , c b C, and p abq is the q-th packet sent from c a to c b . This packet contains w abq bits and is transmitted after the computation time t aq of the originating core (c a ) has elapsed.

The set of all packets sent from c a to c b is P ab .

The CDCG represents the communication and computation for an application composed by an arbitrary number of cores. The direction of the edges in this graph denotes that the destination vertex computation depends on the computation of the origin vertex. In other words, the destination vertex presents a communication dependence w.r.t. the origin vertex.

CWM and CDCM are evaluated here using a mesh topology NoC using wormhole, deterministic XY routing algorithm. Other NoC topologies can be equally treated.

Definition 3:

A communication resource graph is a directed graph CRG = < , L>, where the vertex set is the set of tiles = { 1 , 2 ,…, n }, and the edge set

L = {( i , j ), i , j
} gives the set of paths from i to j . The value n is again the total number of tiles and is equal to , the product of the two NoC dimensions. CRG edges and vertices represent physical links and routers of the target architecture, respectively. The CRG definition is equivalent to the architecture characterization graph in [4] and to the NoC topology graph in [5].

Figure 1 illustrates the above definitions using a hypothetical application with four IP cores exchanging a total of six packets and a 2 2 NoC. Figure 1

(c) CRG 1 =<{ 1 , 2 , 3 , 4 }, {( 1 ,B),( 2 ,A),( 3 ,F),( 4 ,E)}> and (d) CRG 2 =<{ 1 , 2 , 3 , 4 },{( 1 ,B),( 2 ,E),( 3 ,F), ( 4 ,A)}>. Start End (b) (d) (a

Energy and Timing Model

Energy consumption originates from both IP cores and interconnection. This work focuses on NoC energy consumption, proposing models to estimate energy dissipation, used as cost functions to evaluate mappings. Dynamic energy consumption is proportional to switching activity, arising from packets moving across the NoC, dissipating energy on the interconnect wires and inside each router. Static energy consumption comes mainly from leakage current, and is proportional to the application execution time and to the number of gates. Normally, static energy contributes with the smallest part of the energy consumption. However, for deep-submicron technologies, leakage current cannot be neglected, and static energy consumption is a significant part of the total energy consumption, reaching up to 20% in new technologies [8].

While CWMs are only suitable to compute dynamic energy, CDCMs can be used to estimate the total energy consumption of the NoC, enabling the computation of static and dynamic power dissipation figures. The approach here is similar to those in [4] and [5], but some concepts are extended to consider static energy consumption. The concept of bit energy EBit [6] is used to estimate the dynamic energy consumption for each bit, when these flip its polarity from a previous value. EBit is split into three components: bit dynamic energy (ERbit) dissipated at the router (router wires, buffers and logic gates); bit dynamic energy dissipated on horizontal (ELHbit) and vertical (ELVbit) links between tiles; and bit dynamic energy (ECbit) dissipated on the links between the router and the core of the tile. The relationship between these quantities is expressed by Equation ( 1), which computes the dynamic energy consumption of a bit passing through a router and a vertical or horizontal link.

(1) EBit = ERbit + (ELHbit or ELVbit) + ECbit ERbit depends on the buffer structure and technology to estimate how many bit flips occur to write, read and preserve the information. ELbit is directly proportional to the tile dimension. For regular mesh NoCs with square tiles, it is reasonable to estimate that ELHbit and ELVbit have the same value. Therefore, ELHbit and ELVbit are represented simply by ELbit. Furthermore, for large tile dimensions, ECbit is negligible w.r.t. ELbit.

Equation ( 2) computes the dynamic energy consumed by a single bit traversing the NoC, from tile i to tile j , where corresponds to the number of routers through which the bit passes.

(2) EBit ij = ERbit + ( -1) ELbit Let i and j be the tiles to which c a and c b , are respectively mapped. Then, the dynamic energy consumed by a c a c b communication is given by EBit ab = w ab EBit ij . Equation (3) gives the total amount of NoC dynamic energy consumption (EDyNoC) for CWM, computing this for all bits of all communications (y) occurring in the NoC.

(3)

EDyNoC (CWM) = y 1 i (i) E ab Bit , c a ,c b C
The total NoC energy consumption (ENoC) using CWM corresponds to EDyNoC (CWM) , since CWM is a model that does not carry timing information. It should be clear from this reasoning that CWM is inappropriate to compute static energy consumption.

Let w abq be the total amount of bits of a packet p abq P ab , going from core c a to core c b . Then, the dynamic energy consumed by the q-th packet of a c a c b communication is given by EBit abq = w abq EBit ij . Hence, equation ( 4) gives EDyNoC for CDCM and k abi represents the number of packets of the i-th communication from core c a to core c b . ( 4)

EDyNoC (CDCM) = abi abq k 1 q y 1 i (i) Bit E c a ,c b C
The static power of each router (PSRouter) is proportional to the number of gates that compose the router and that can be estimated by electrical simulation. With n representing the number of tiles, equation ( 5) computes NoC static power consumption (PstNoC). ( 5) PStNoC = n PSRouter

The total packet delay of the wormhole routing algorithm is composed by the routing delay and by the packet delay. The routing delay is the time necessary to create the communication path, which is determined during the traffic of the first flit, where the header of the packet is placed. The packet delay depends on the number of remaining flits. Let n abq be the number of flits of the q-th packet from c a to c b , obtained by dividing w abq by the link width. Let be the period of a clock cycle, and let t r be the number of cycles needed for taking a routing decision inside a router. Also, let t l be the number of cycles needed to transmit a flit through a link (between tiles or between an IP core and a router). The routing delay (d Rijq ) and the packet delay (d Pijq ) of the q-th packet from i to j , are represented in Equations ( 6) and (7), considering that a packet goes through routers without contention. Contentions can only be determined at execution time. (6)

d Rijq = ( (t r + t l ) + t l ) (7) d Pijq = (t l (n abq -1))
The total packet delay (d ijq ), obtained from the sum of (d Rijq ) and (d Pijq ), is expressed by Equation ( 8). ( 8)

d ijq = ( (t r + t l ) + t l n abq )
The application execution time (texec) is obtained from both the computation of all t iq and d ijq , and the contention time. Static energy consumption is proportional to PstNoC and to texec. Thus, equation ( 9

Communication Algorithms Comparison

The FRW framework implements a simulated annealing search method to obtain mapping solutions for CWM and CDCM. Moreover, it can also execute an exhaustive search method to compare the quality of solutions against an absolute optimum solution, for small NoCs. Both algorithms, CWM and CDCM, start from an initial mapping, evaluate the mapping cost, and search for a new mapping that reduces the computed cost, until reaching a stop condition. For both algorithms, the mapping cost is stored inside cost variables of CRG edges and vertices. The sum of all cost variables determines the mapping objective function. Initially, all cores of C are randomly mapped onto the set of tiles and all cost variables are set to 0.

Given a mapping function, let tiles i and j be the respective images of cores c a and c b in this function. For the CWM algorithm, each total number of bits in all packets of c a c b communications (w ab ) is associated to the cost variable of the corresponding vertices and edges of CRG. These start at i , follow the path defined by the XY routing algorithm and end at j . Equation ( 3) is the objective function of CWM, used to evaluate the cost of each mapping. The cost variable of each CRG edge is used to compute the dynamic energy of a link by multiplying w ab by ELbit. The cost variable of each CRG vertex is used to compute the dynamic energy of a router by multiplying w ab by ERbit. The sum of all cost variables of CRG results in EDyNoC. This procedure is illustrated in Figure 2. CWM aims to find mappings reducing EDyNoC.

CDCG improves CWM, as it considers computation time and packet ordering. Dependent packets cannot be concurrent. However, independent packets can occur at the same time, and may consequently lead to package contention, if they share the same communication resource. Packet contention implies larger texec, leading to increased EStNoC. The CDCM algorithm searches for mappings that minimize the sharing of communication resources for concurrent packets. To allow this, each CRG edge and vertex are associated to a cost variable list, where each element of the list represents a packet that contains the number of bits, the absolute time interval that the packet is occupying the NoC resource, the source and target tiles.

The CDCM algorithm starts with all vertices pointed by the Start vertex (Figure 1(b)). Pointing to a CDCG vertex implies that the packet enclosed into the vertex may be executed onto the CRG. A vertex execution marks all of its output edges as free. A vertex can only be executed if all of its input edges are free. The algorithm searches for all free vertices following all paths, until all paths reach the End vertex. This procedure computes the time of each dependent path in the CRG. To execute a CDCG vertex that encloses the q-th packet onto the CRG, all bits of the q-th packet of c a c b communication (w abq ) are associated with the corresponding vertices and edges of CRG, starting from i , following the XY routing algorithm, and ending in j .

When two or more packets compete for the same resource, all concurrent packets have to be contained into router input buffers. Therefore, from the contention point until reaching the target tile, contention time is added to the elapsed time, enabling to estimate the total packet delay more accurately. This procedure is illustrated in Figure 3 and Figure 4, where the communication weight of A F is stored into the buffer of router 1 . The cost variable list of CRG edges and vertices computes the dynamic energy of links and routers by multiplying w abq by ELbit and by ERbit, respectively. For a given mapping, the sum of all cost variables generates the value EDyNoC. When all packets of CDCG are executed into CRG texec is obtained, enabling to compute ENoC through equation ( 10), the objective function to evaluate the mapping cost. The goal of CDCM is to find mappings that minimize ENoC.

For both, CWM and CDCM algorithms, if the mapping cost achieved with a new mapping is smaller than the previously stored, the current mapping and cost are saved for further comparison. If the stop condition has not been reached, a new mapping is chosen and cost re-evaluated.

CWM and CDCM are able to estimate EDyNoC accurately, since this value depends only on the bit traffic along the NoC. The essential difference between the models is that CWM is not appropriate to estimate texec, due to the absence of task computation time and the impossibility of evaluating packet contention. The main advantages of CWM are (i) easy extraction of application core graph (CWG), since this can be done by simulation techniques; and (ii) low computational complexity. The automatic extraction of the CDCM application core graph is a hard task, since simulations only allow extracting possible message orderings, and not the required message dependence. As a consequence, CDCGs are described by hand, increasing the error susceptibility. The greater complexity of CDCM is observed on its algorithmic implementation, which increases computation time and memory usage w.r.t. CWM. However, CDCM captures the message ordering, and thus allows estimating the absolute execution time, and consequently the EStNoC. The main drawback of CDCM evaluation is that in real embedded applications, the number of packets between cores is much larger than the number of cores. Since each vertex of CDCG represents a packet exchanged between two cores and each vertex of CWG represents a core, CDCGs are larger than CWGs. A comparison between the CPU time required by CDCM and CWM algorithms is presented in Section 5.

Energy and Timing Analysis

This Section shows the application of algorithms CWM and CDCM to the example application mappings of Figure 1. The CWM evaluation does not capture the differences between the two mappings. However, CDCM evaluation shows that these mappings imply different execution times and distinct energy consumption estimates.

Figure 2 illustrates how the internal model of the CWM algorithm represents the energy for the two mappings of Figure 1(c,d). For illustration purposes, the example assumes that ERbit = ELbit = 1 10 -12 J/bit. The mapping of CWG onto CRG results that each vertex and edge of CRG is annotated with ERbit and ELbit energies, respectively. For instance, the execution of E A communication of Figure 1(a) onto the mapping of Figure 2(a) implies 35 10 -12 J of energy consumption, which is computed in tiles 4 and 2 , and in the link between these tiles. Applying equation ( 3), over all communications of Figure 1(a), EDyNoC results in 390 10 -12 J for both mappings. Energy = 400 10 -12 J Execution time = 100 ns Energy = 399 10 -12 J Execution time = 90 ns 15(A B): [10,26] 40(B F): [11,52] *15(A F): [46,69] 15(F B): [83,99] 40(B F):[14,55] *15(A F): [56,72] 15(F B): [80,96] 15(A B): [START_REF] Dick | TGFF: task graphs for free[END_REF]24] 15(A F): [45,60] 40(B F):[13,53] *15(A F): [55,70] 0 0 0 15(A B): [START_REF] Ye | Packetization and routing analysis of on-chip multiprocessor networks[END_REF]23] 20(E A): [14,35] 15(E A): [60,76] 15(A F):[43,59] 0 2 1 15(A B): [START_REF] Ye | Analysis of power consumption on switch fabrics in network routers[END_REF]21] 20(E A): [16,36] 15(E A): [62,77] 15(A F): [42,57] B 15(A B): [12,27] 40(B F): [10,50] 15(F B): [85,100] A 20(E A): [11,32] 15(E A): [57,73] E 20(E A): [10,30] 15(E A): [56,71] 20(E A): [13,33] 15(E A): [59,74] 15(F B): [82,97] 40(B F): [16,56] *15(A F): [58,73] 15(F B):[79,94] F 4 3 15(A B): [13,29] 40(B F): [11,52] 15(F B): [73,89] 15(A B): [10,26] 40(B F): [14,55] 15(A F): [46,62] 15(F B): [70,86] 40(B F): [13,53] 0 0 0 20(E A): [11,32] 15(E A):[57,73] 2 1 15(A B): [START_REF] Ye | Analysis of power consumption on switch fabrics in network routers[END_REF]21] 20(E A): [16,36] 15(E A): [62,77] 15(A F): [42,57] B 15(A B): [15,30] 40(B F): [10,50] 15(F B):[75,90] E F 15(A B): [START_REF] Ye | Packetization and routing analysis of on-chip multiprocessor networks[END_REF]23] 20(E A): [14,35] 15(E A): [60,76] 15(A F): [43,59] A 20(E A): [10,30] 15(E A): [56,71] 20(E A): [13,33] 15(E A): [59,74] 15(A B): [12,27] 15(F B): [72,87] 40(B F): [16,56] 15(A F): [48,63] 15(F B): [69,84] 15(A B): [START_REF] Dick | TGFF: task graphs for free[END_REF]24] 15(A F): [45,60] To evaluate CDCM mappings, consider the same parameters stated before, t r = 2 clock cycles, t l = 1 clock cycle, = 1 ns, one-bit sized flits, and unbounded router buffers. Figure 3 shows the same mappings illustrated before, now evaluated with CDCM. Each edge and each vertex is annotated with the number of bits in a given time interval. For instance, in Figure 3(a) tile 4 is annotated with 20(E A):[11,32] and 15(E A):[57,73] which means that there are two packets from E to A. The first packet has 20 bits and occupies router 4 during the interval from 11 ns to 32 ns. This interval is obtained by adding the E computation time t E1 = 10 ns with 1 ns of the link time. The link from 4 to 2 is annotated with the same two packets, each one delayed by the router delay.

Figure 4 depicts the CDCM algorithm execution over the mapping of Figure 3(a), showing a timing diagram that illustrates all computations and all packet deliveries. During A F and B F packets transmission some contentions occur, since they compete for the same resources at the same time. When the B F packet uses the router of tile 1 , the A F packet is contained into the input buffer of 1 . This contention implies that the packet is delayed until the B F packet is transmitted to tile 3 . The effect of this delay is observed in the router variables of 1 and the variables of the link between the 3 router and IP core F. These variables are marked with '*' in Figure 3(a). Figure 5 shows the timing diagram corresponding to the mapping depicted in Figure 3(b). It is possible to observe that this mapping avoids contention, since there are no packets competing for the same link at the same time. Changing the mapping of Figure 3(a) to that of Figure 3(b) implies an execution time reduction of 11.1%, from 100 ns to 90 ns. As there are different execution times for both mappings, the consumed static energy is different, too. For instance, consider PstNoC=0.1 10 -12 J/ns. Applying equation (10), ENoC can be obtained for both mappings, showing that mapping (a) consumes 1 % more energy than (b). This difference cannot be computed by CWM, since it just captures the effects of dynamic energy.

Experimental Results

Table 1 summarizes the characteristics of 18 applica-tions mapped onto 8 different NoC sizes. There are 4 embedded applications (a distributed Romberg integration, an 8-point Fast Fourier Transform, and 2 image applications for object recognition and image encoding) with some variations, for a total of 8 embedded applications. The remaining applications are benchmarks randomly generated by a proprietary system, which is similar to TGFF [9]; however, the system describes benchmarks through CDCGs, representing message dependence and bit volume of each message. The NoC size is the number of CRG vertices, the number of cores corresponds to the number of CWG vertices, and the number of packets of all cores equals the number of CDCG vertices. For each application, the best mapping achieved with the CWM algorithm is compared to the best mapping achieved with the CDCM algorithm, and the results are summarized in Table 2. For both models exhaustive search (ES) and simulated annealing (SA) were applied, depending on the NoC size. ETR gives the average execu- tion time reduction, and ECS denotes the average energy consumption saving, for a given technology, when mappings achieved with the CDCM algorithm are compared to the ones achieved with the CWM algorithm. ECS 0.35 column refers to ECS values obtained from 0.35 technology, and ECS 0.07 column refers to ECS values obtained by estimating 0.07 technology [8]. First, for small NoC sizes (up to 3x4 or 2x5), both ES and SA methods reached the same results. For larger NoC sizes (8x8, 10x10 and 12x10), it is not possible to find optimum mappings with ES within a reasonable time. The ETR column shows that the CDCM algorithm results in 40% average reduction of execution time when compared to the CWM algorithm. The ECS 0.35 column illustrates a very small energy consumption saving, since the leakage current is not as important for this technology. However, for deep-submicron technologies, there is a significant reduction in energy consumption (20% in average), as observable in column ECS 0.07 . In addition, Table 2 shows a slight trend of energy consumption saving and execution time reduction when the NoC size increases.

The computational complexity of the CWM algorithm is proportional to the number of communications between cores (NCC) and the computational complexity of CDCM algorithm is proportional to the number of dependences and packets of all cores (NDP). In real embedded applications, NDP is larger than NCC. However, the increase in CPU time with the increase of the NDP/NCC ratio is approximately linear with a small slope. The worst case for CDCM took only 23% more CPU time than for CWM.

Conclusions and Future Work

This paper addressed the problem of mapping applications onto NoCs. A communication dependence and computation model (CDCM) has been introduced and compared to the communication weighted model (CWM). As a conclusion CDCM is able to reduce the application execution time and energy consumption when compared to CWM. Experimental results show an average of 40% in execution time reduction. CDCM also reduces energy consumption. For a 0.07 technology, an average of 20% in energy savings is obtained. Moreover, CDCM presents only a moderate increase in computational cost, when compared to CWM, with better mapping results.

Definition 1 :

 1 illustrates them. A communication weighted graph (CWG) is a directed graph <C, W>. The set of vertices C = {c 1 , c 2 …, c n } represents the set of cores in one application. Assuming w ab is the number of bits of all packets sent from a core c a to a core c b , then the set of edges W is {(c a , c b ) | c a , c b C and w ab 0}, and each edge is labeled with the value w ab .

  (a) shows a CWG where C = {A,B,E,F}, the edge labels are w AB = 15, w AF = 15, w BF = 40, w EA = 35, w FB = 15 and the set W can be extracted easily from the Figure.

Figure 1 (

 1 b) depicts one possible CDCG for the application, where P = {p EA1 = (E,A,10,20), p EA2 = (E,A,20,15), p AF1 = (A,F,6,15),…} and D = {(Start,p EA1 ), (p EA1 , p EA2 ), (p AB1 , p AF1 )…}. Figure 1(c) and Figure 1(d) depict two arbitrary mappings of C, each corresponding to a CRG as follows:

Figure 1 -

 1 Figure 1 -(a) CWG; (b) CDCG; (c,d) 2 mappings.

  ) computes NoC static energy consumption (EstNoC). (9) EStNoC = PStNoC texec Finally, equation (10) gives the overall static plus dynamic energy consumption at the NoC (ENoC) for CDCM.(10) ENoC (CDCM) = EStNoC + EDyNoC (CDCM)

Figure 2 -

 2 Figure 2 -Mappings for Figure 1(c, d) -estimation of energy consumption obtained with CWM.

Figure 3 -

 3 Figure 3 -Figure1(c, d) mappings, using CDCM.To evaluate CDCM mappings, consider the same parameters stated before, t r = 2 clock cycles, t l = 1 clock cycle, = 1 ns, one-bit sized flits, and unbounded router buffers. Figure3shows the same mappings illustrated before, now evaluated with CDCM. Each edge and each vertex is annotated with the number of bits in a given time interval. For instance, in Figure3(a) tile 4 is annotated with

Figure 4 -

 4 Figure 4 -Timing for Figure 3(a) mapping.Figure5shows the timing diagram corresponding to the mapping depicted in Figure3(b). It is possible to observe that this mapping avoids contention, since there are no packets competing for the same link at the same time. Changing the mapping of Figure3(a) to that of Figure3(b) implies an execution time reduction of 11.1%, from 100 ns to 90 ns. As there are different execution times for both mappings, the consumed static energy is different, too. For instance, consider PstNoC=0.1 10 -12 J/ns. Applying equation (10), ENoC can be obtained for both mappings, showing that mapping (a) consumes 1 % more energy than (b). This difference cannot be computed by CWM, since it just captures the effects of dynamic energy.

Figure 5 -

 5 Figure 5 -Timing for Figure 3(b) mapping.

Table 1 -Summary of NoC/application features.

 1

	NoC size Number of	Number of packets	Total volume of bits during
		cores	of all cores	application execution
	3 x 2	5; 6; 6	43; 17; 43	78,817; 174; 49,003
	2 x 4	5; 7;8	16; 33; 18	1,600; 23,235; 5,930
	3 x 3	7; 9; 9	16; 18; 32	1,600; 1,860; 43,120
	2 x 5	8; 9; 10	24; 51; 22	2,215; 23,244; 322,221
	3 x 4	10; 12; 14	15; 25; 88	3,100; 2,578,920; 115,778
	8 x 8	62	344	9,799,200
	10 x 10	93	415	562,565,990
	12 x 10	99	446	680,006,120

Table 2 -Average energy and execution time re- ductions for CWM and CDCM.

 2

	Algorithm	NoC size	ETR	ECS 0.35	ECS 0.07
		3 x 2	36 %	0,50 %	15 %
		2 x 4	27 %	0,43 %	13 %
		3 x 3	39 %	0,55 %	17 %
		2 x 5	42 %	0,72 %	23 %
		3 x 4	42 %	0,71 %	22 %
	Simulated an-nealing only	8 x 8 10 x 10 12 x 10	38 % 46 % 48 %	0,6 % 0,8 % 0,86 %	19 % 25 % 26 %
	Average		40 %	0,65 %	20 %

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05) 1530-1591/05 $ 20.00 IEEE