Mirosław Lachowicz

Martin Parisot
email: martin.parisot@inria.fr

Zuzanna Szymańska

Intracellular protein dynamics as a mathematical problem

. The model takes into account diffusive transport in the nucleus and cytoplasm, as well as active transport of protein molecules along microtubules in the cytoplasm. The model reproduces, at least in numerical simulations, the oscillatory changes in protein concentration observed in the experimental data. To our knowledge this is the first paper that, in the multidimensional case, deals with a rigorous mathematical analysis of a model of intracellular dynamics with active transport on microtubules. In particular, in the present paper we prove well-posedness of the model in any space dimension. The model is a complex system of nonlinear PDEs with specific boundary conditions. It may be adapted to other signaling pathways.

Introduction

Mathematical modeling of intracellular proteins dynamics is one of the most commonly undertaken challenges in mathematical biology. Of course, this follows from our belief that increased knowledge about these processes will allow us to better understand the functioning of a cell, especially if it is distorted, as it is in the case of many diseases like cancer. In the literature one can find many models describing the intracellular dynamics of different proteins, e.g. [START_REF] Gordon | Spatio-temporal modelling of the p53-mdm2 oscillatory system[END_REF][START_REF] Lipniacki | Mathematical model of nf-kappab regulatory module[END_REF][START_REF] Puszynski | Crosstalk between p53 and nuclear factor-b systems: pro-and anti-apoptotic functions of nf-b[END_REF][START_REF] Rybinski | Modelling the efficacy of hyperthermia treatment[END_REF][START_REF] Szymańska | Mathematical modelling of heat shock protein synthesis in response to temperature change[END_REF][START_REF] Sturrock | Spatiotemporal modelling of the Hes1 and p53-Mdm2 intracellular signalling pathways[END_REF]. However, most of these models, either ignores the spatial aspect, or treat it in a very simplified way, taking into account only the diffusion as a transport mechanism of molecules. Only recently new models that take into account directed transport that is held in the cytoplasm were published [START_REF] Szymańska | Mathematical modeling of the intracellular protein dynamics: the importance of active transport along microtubules[END_REF][START_REF] Cangiani | A spatial model of cellular molecular trafficking including active transport along microtubules[END_REF][START_REF] Sturrock | Influence of the nuclear membrane, active transport, and cell shape on the Hes1 and p53-Mdm2 pathways: Insights from spatio-temporal modelling[END_REF][START_REF] Chaplain | Spatio-Temporal Modelling of Intracellular Signalling Pathways: Transcription Factors, Negative Feedback Systems and Oscillations[END_REF].

One of few such models, was proposed in 2014 by Szymańska at al., the model of a simple intracellular protein concentration control system, that took into account active transport along microtubules [START_REF] Szymańska | Mathematical modeling of the intracellular protein dynamics: the importance of active transport along microtubules[END_REF]. The work focused not on particular protein dynamics but on the proper description of intracellular transport processes. The main goal was to provide a description that is fuller in terms of biology and correct from the mathematical point of view. Therefore the model is generic, however, corresponds to protein Hes1 simple regulatory system.

In the present paper, we show a well-posedness of a model of intracellular transport processes generalizing the system presented in [START_REF] Szymańska | Mathematical modeling of the intracellular protein dynamics: the importance of active transport along microtubules[END_REF] to the case of N couples mRNA/protein. The mathematical analysis is perform for an arbitrary space dimension and with non-linear parameters satisfying classical mathematical properties. To our knowledge, this is the first work that in multidimensional case deals with rigorous mathematical analysis of the model of intracellular protein dynamics including active transport along microtubules. The model can be describe as a system of hyperbolic and parabolic equation coupled through non-linear parameters and boundary condition. The mathematical analysis of such system is not trivial and we propose a analysis based on classical results for, on the one hand linear hyperbolic equation, on the second linear parabolic equation. The key point of the proof of well-posedness for the coupled nonlinear system lies in a priori estimations leading to contraction argument. Last but not least, some properties such as global evolution and positivity are given to ensure the biological relevance of the model.

For other paper with mathematical analysis of a gene regulatory network model in a 1-dimensional domain including linear stability analysis of the steady states we refer the reader to recent work of Chaplain at al. [START_REF] Chaplain | Hopf bifurcation in a gene regulatory network model: Molecular movement causes oscillations[END_REF]. The authors showed that a diffusion coefficient acts as a bifurcation parameter and give rise to a Hopf bifurcation. They however did not study directly the influence of microtubules.

The present paper is organized as follows. In Sections 2 we briefly introduce the mathematical model. More details of the biological background as well as its discussion may be found in [START_REF] Szymańska | Mathematical modeling of the intracellular protein dynamics: the importance of active transport along microtubules[END_REF]. Section 4 is devoted to the mathematical analysis of the model. Finally we discuss about further improvements in Section 5.

Intracellular dynamics

Here we present a very brief description of the underlying biological process. If interested in more details, please refer to our work Szymanska at al. [START_REF] Szymańska | Mathematical modeling of the intracellular protein dynamics: the importance of active transport along microtubules[END_REF] or, to a textbook of molecular biology by Alberts at al. [START_REF] Alberts | Essential Cell Biology[END_REF].

In simplest terms, an eukaryotic cell consists of a nucleus and cytoplasm, and the nucleus is located inside the cytoplasm, see Fig. 1. It means that considered domain, i.e. cell volume denoted by Ω, is composed of two separated parts -nucleus denoted by N and cytoplasm denoted by C . The entire cell is limited by a so-called cellular membrane here denoted by Γ c . The nucleus is separated from the cytoplasm by a membrane, which is called the nuclear envelope that is denoted by Γ n . Both on the cell membrane and the nuclear envelope we assume boundary conditions. In the cytoplasm there are so called microtubules, which are thin polymers arranged concentrically from the centre to the edges of the cell. They are used to give shape to the cell, and to facilitate intracellular transport of molecules and organelles. In the model the microtubules density is described by ξ function. The nucleus contains DNA, which is used for the production of mRNA molecules. The newly produced mRNA molecules escape from the nucleus to the cytoplasm. After passing the nuclear envelope, with help of the microtubules, they are moving into the cytoplasm, where they are released. Thus, in the present paper the mRNA in the cytoplasm is divided into associated with microtubules and free denoted by m i and m i , i = 1, . . . , N, respectively. In the cytoplasm there are many molecular machineries called ribosomes where, on the basis of mRNA, proteins are produced. Newly synthesised proteins diffuse in the cytoplasm where they perform its appropriate functions. Some of the proteins are needed in the nucleus in order to regulate the process of mRNA production. In such a case they bind to the microtubules and are transported along them to the nucleus. As in the case of mRNA, the protein can be divided into the associated with microtubules denoted by p i and free, denoted by p i , i = 1, . . . , N. In the model, the protein produced in ribosomes in the cytoplasm acts also as its own inhibitor. That means protein molecules are transported into the nucleus where they act as a suppressor for its own mRNA production. Concentration of mRNA in the nucleus is denoted by m i and concentration of protein there is denoted by p i , i = 1, . . . , N. This mechanism is called negative feedback and allows to maintain the limited protein concentration.

The system of equations

The cell volume, defined by Ω is a nonempty convex bounded and smooth domain in R d with d ≥ 1. In what follows we do not impose any restriction on the size of the space dimension d. The boundary, denoted by Γ c , corresponds to the cell membrane and c is the unit normal vector outward the cell. We defined the nucleus of the cell by N , a nonempty bounded and smooth sub-domain in Ω. Its boundary Γ n corresponds to the nuclear envelope, n is the unit normal vector outward the nucleus and we assume that Γ n ∩ Γ c = ∅. The cytoplasm C is defined as the complement of the nucleus in the cell, i.e. C = Ω \ N . Here and in what follows we adhere to the convention that given set U we denote by U T = [0, T ] × U, and T is arbitrary but fixed. In addition, the cytoplasm of eukaryote cells presents numerous microtubules which are thin filaments going from the nucleus envelope to the cell membrane. To model the microtubules, we introduce the unit vector η : C → R d oriented from the nucleus envelope to the cell membrane such that the microtubules are tangent to η. From now on, we denote by t > 0 the time and by x ∈ Ω the position in the cell.

In the present paper we consider the system of nonlinear PDEs with specific boundary conditions that can be referred to specific signaling pathway composed by N couples of mRNA and protein. The main novelty of the Γ n -nuclear envelope; ξ -concentration of microtubules; m i -mRNA molecules in the nucleus; p i -protein molecules in the nucleus; m i -mRNA molecules in the cytoplasm that are linked to the microtubules; m i -mRNA molecules in the cytoplasm that are free, i.e. not linked to the microtubules; p i -protein molecules in the cytoplasm that are free, i.e. not linked to the microtubules, i = 1, . . . , N;

model presented in [START_REF] Szymańska | Mathematical modeling of the intracellular protein dynamics: the importance of active transport along microtubules[END_REF] comes from the distinction of two ways of motion in the cytoplasm.We propose to model the distribution of microtubules through the function ξ on C and such that

0 ≤ ξ (t, x) ≤ 1 .
The density of microtubules is supposed to be given at any time thus it is not an unknown of our model. The movement of molecules in the cytoplasm depends on whether the molecules are bound to the microtubules or not. Therefore, in the cytoplasm, both species, protein and mRNA are divided into those associated with microtubules, namely linked molecules, and those not associated, namely free molecules. For a given t and x we denote by m i (t, x) and p i (t, x) the concentration of mRNA linked to the microtubules and the concentration of protein linked to the microtubules, respectively. Analogously, by m i (t, x) and p i (t, x) we denote the concentration of mRNA and the concentration of protein, respectively, that are free to the microtubules. Thus, we consider the following equations in the cytoplasm

on C T          ∂ t m i -∇ • k i ∇ m i = -i m i + b i m i , (1a) ∂ t m i + ∇ • (v i m i η) = -i m i -b i m i , (1b) ∂ t p i -∇ • κ i ∇ p i = -λ i p i -β i p i + α i , (1c)
∂ t p i -∇ • ν i p i η = -λ i p i + β i p i , (1d)
with the nonlinear parameters that are local functions of the unknown concentrations and the density of microtubules. More precisely, the diffusion parameters k i (ξ; m c ; p c ) and κ i (ξ; m c ; p c ) are square matrices, the speed parameters v i (ξ) and ν i (ξ), the parameters i (ξ; m c ; p c ) and λ i (ξ; m c ; p c ), the exchange parameters b i (ξ; m c ; p c ) and β i (ξ; m c ; p c ), the transcription rate a i (ξ; m c ; p c ) and the translation rate α i (ξ; m c ; p c ) are positive scalar function of the vectorial unknowns

m c = m i 1≤i≤N , (m i ) 1≤i≤N , p c = p i 1≤i≤N , p i 1≤i≤N .
In the nucleus, i.e. if x ∈ N , the concentrations of the i th couple are denoted by m i (t, x) for mNRA and p i (t, x) for protein. We assume that the motion in the nucleus is govern by the diffusion and the evolution of the concentrations are given by the following reaction-diffusion equations

on N T ∂ t m i -∇ • (k 0i ∇m i ) = -0i m i + a 0i , (2a) ∂ t p i -∇ • (κ 0i ∇p i ) = -λ 0i p i , (2b)
with the nonlinear parameters that are local functions of the unknown concentrations. More precisely, the diffusion operators k 0i (m n ; p n ) and κ 0i (m n ; p n ) are square matrices, the decay parameters 0i (m n ; p n ) and λ 0i (m n ; p n ) and the transcription rate a 0i (m n ; p n ) are positive scalar function of the vectorial unknowns m n = (m i ) 1≤i≤N and p n = (p i ) 1≤i≤N . At the nuclear envelope, i.e. if x ∈ Γ n , we suppose that the exchange is only one-way. More precisely the mRNA concentration m i leaving the nucleus is split into the subpopulations with concentrations m i and m i accordingly to the density of microtubules ξ at the nuclear membrane. The equations read as

on Γ n T    n • (k 0i ∇m i ) = -u i m i , (3a) -n • k i ∇ m i = (1 -ξ) u i m i , (3b) n • η v i m i = ξu i m i . (3c)
The protein concentration in the nucleus p i results from the contribution of two ways of motion, i.e.

on Γ n T n • (κ 0i ∇p i ) = µ i p i + n • η ν i p i , (3d) -n • κ i ∇ p i = -µ i p i , (3e)
with the speed passing through the nuclear nuclear envelope of mRNA u i (ξ) and of protein µ i (ξ).

Finally we consider no exchange of particles between the cell and its environment. It leads to the following boundary condition at the cell membrane

on Γ c T      c • k i ∇ m i = c • η v i m i (4a) -c • κ i ∇ p i = 0, (4b) c • η ν i p i = 0. ( 4c
)
Note that the boundary conditions (4) imply that the mRNA particles reaching the cell membrane via the advection are released into the cytoplasm and change its way of motion. We refer to [START_REF] Szymańska | Mathematical modeling of the intracellular protein dynamics: the importance of active transport along microtubules[END_REF] for details of the biological interpretation of Eqs. ( 1)-( 2) and the boundary conditions ( 3)-(4).

Analysis

In the following section, we give the main result of the paper. We prove the global well-posedness of the nonlinear and multi-dimensional model. The model is a system of non-linear coupled hyperbolic and parabolic equations and its well-posedness, especially globally in time, is not obvious. The proof is based on an adapted iterative approximation process which permits to use the classical results for the linear hyperbolic equations as well as those for the linear parabolic equations. The convergence of the approximations leading to the global well-posedness is obtained using classical contraction arguments.

Theorem 1. Let d ≥ 1 and T > 0 be arbitrary and fixed. Assume that for any

1 ≤ i ≤ N (T1h1) ξ ∈ C 1 (C T );
(T1h2) There exists ε > 0 such that for any φ ∈ R N , ψ ∈ R N , the eigenvalues of k 0i (φ; ψ) and κ 0i (φ; ψ) are larger than ε, and for any ξ ∈ [0, 1], φ ∈ R 2N and ψ ∈ R 2N the eigenvalues of k i (ξ; φ; ψ) and κ i (ξ; φ; ψ) are larger than ε;

(T1h3) For any ξ ∈ [0, 1], there exists U < ∞ such that 0 ≤ u i (ξ) ≤ U and 0 ≤ µ i (ξ) ≤ U; (T1h4) For any ξ ∈ [0, 1], we have v i (ξ) ≥ 0 with v i ∈ C 1 ([0, 1]) and ν i (ξ) ≥ 0 with ν i ∈ C 1 ([0, 1]); (T1h5) The following parameters are bounded 0i < ∞ , λ 0i < ∞ , a 0i < ∞ ,
where 0i = 0i (0; 0), λ 0i = λ 0i (0; 0), a 0i = a 0i (0; 0); and for any ξ

∈ [0, 1] i < ∞ , λ i < ∞ , α i < ∞ , b i < ∞ , β i < ∞ ∀ ξ ∈ [0, 1] ,
where i (ξ) = i (ξ; 0; 0), λ i (ξ) = λ i (ξ; 0; 0), α i (ξ) = α 0i (ξ; 0; 0), b i (ξ) = b i (ξ; 0; 0) and β i (ξ) = β i (ξ; 0; 0);

(T1h6) There exists A > 0 such that for any (m n ; p n ;

n n ; q n ) ∈ R N 4 | 0i (m n ; p n ) m i -0i (n n ; q n ) n i | ≤ A (m n -n n ; p n -q n ) N , |λ 0i (m n ; p n ) m i -λ 0i (n n ; q n ) n i | ≤ A (m n -n n ; p n -q n ) N , |a 0i (m n ; p n ) -a 0i (n n ; q n )| ≤ A (m n -n n ; p n -q n ) N ,
where (m n ; p n ) N = N j=1 m j + p j and for any ξ ∈ [0, 1] and (m c ; p c ;

n c ; q c ) ∈ R 2N 4 i (ξ; m c ; p c ) m i -i (ξ; n c ; q c ) n i ≤ A (m c -n c ; p c -q c ) C , | i (ξ; m c ; p c ) m i -i (ξ; n c ; q c ) n i | ≤ A (m c -n c ; p c -q c ) C , |b i (ξ; m c ; p c ) m i -b i (ξ; n c ; q c ) n i | ≤ A (m c -n c ; p c -q c ) C , λ i (ξ; m c ; p c ) p i -λ i (ξ; n c ; q c ) q i ≤ A (m c -n c ; p c -q c ) C , λ i (ξ; m c ; p c ) p i -λ i (ξ; n c ; q c ) q i ≤ A (m c -n c ; p c -q c ) C , β i (ξ; m c ; p c ) p i -β i (ξ; n c ; q c ) q i ≤ A (m c -n c ; p c -q c ) C , |α i (ξ; m c ; p c ) -α i (ξ; n c ; q c )| ≤ A (m c -n c ; p c -q c ) C ,
where (m c ; p c ) C = N j=1 m j + m j + p j + p j .

Then, for all initial data such that

m 0 i ∈ H 1 (N ) , m 0 i ∈ H 1 (C ) , m 0 i ∈ L 2 (C ) , p 0 i ∈ H 1 (N ) , p 0 i ∈ H 1 (C ) , p 0 i ∈ L 2 (C ) ,
there exists an unique global weak solution of (1)-( 2) with the boundary condition (3) and (4) such that

m i ∈ L ∞ 0, T ; L 2 (N ) ∩ L 2 0, T ; H 1 (N ) , m i ∈ L ∞ 0, T ; L 2 (C ) ∩ L 2 0, T ; H 1 (C ) , m i ∈ L ∞ 0, T ; L 2 (C ) ∩ L 2 Γ c T ; v i |η • c| dσdt , p i ∈ L ∞ 0, T ; L 2 (N ) ∩ L 2 0, T ; H 1 (N ) , p i ∈ L ∞ 0, T ; L 2 (C ) ∩ L 2 0, T ; H 1 (C ) , p i ∈ L ∞ 0, T ; L 2 (C ) ∩ L 2 Γ n T ; ν i |η • n| dσdt .
Assumptions of Theorem 1 seem complex, in particular for the transcription rate a 0i , the translation rate α i and the decay parameters 0i , λ 0i , i and λ i . However, several models of signalling pathways with negative feedback proposed in the literature satisfies those conditions. For instance, this is the case of Hes1 gene regulatory network (Hes1 is a mammalian protein that suppress transcription) and p53/mdm2 regulatory network (p53 regulates the cell cycle and is one of the mayor tumour suppressors) [START_REF] Sturrock | Influence of the nuclear membrane, active transport, and cell shape on the Hes1 and p53-Mdm2 pathways: Insights from spatio-temporal modelling[END_REF].

In order to prove Theorem 1, we start with the following results:

Lemma 1.
Let Ω ⊂ R d be a smooth bounded domain and be n the outward normal to its boundary ∂Ω. Assume that (L1h1) There exists ε > 0 such that for any (t, x) ∈ Ω T , the eigenvalues of the square matrix K (t, x) are larger than ε;

(L1h2) For any (t, x) ∈ Ω T , we have M (t, x) ≥ 0;

(L1h3) F ∈ L 2 (Ω T ) and G ∈ L 2 (∂Ω T ).
Then, for any initial data φ 0 ∈ L 2 (Ω), there exists an unique weak solution φ in L ∞ 0, T ; L 2 (Ω) ∩ L 2 0, T ; H 1 (Ω) of the following parabolic equation with the Robin boundary condition

∂ t φ -∇ • (K∇φ) + εφ = F on Ω T n • (K∇φ) + Mφ = G on ∂Ω T (5)
Proof. Let (ψ n ) n≥0 be a dense set in H 1 (Ω) with ψ 0 = φ 0 and H 1 n = Vect (ψ i ) 0≤i≤n . We consider the following approximation of the variational formulation of (5): For any initial condition in φ 0 ∈ H 1 (Ω), find φ n ∈ C 1 0, T ; H 1 n such that for any ψ ∈ H 1 n we have

Ω ψ∂ t φ n dx + Ω ∇ψ • K∇φ n dx + ε Ω ψφ n dx + ∂Ω Mψφ n dσ = Ω ψFdx + ∂Ω ψGdσ .
The approximated problem is a system of linear ordinary differential equations in the finite dimensional space H 1 n . We conclude that it is well-posed for any time T > 0. By setting ψ = φ n , using (L1h2), the Cauchy-Schwarz's inequality, the trace theorem and the binomial expansion, we obtain the following a priori estimate

∂ t φ n 2 L 2 (Ω) + ε φ n 2 H 1 (Ω) ≤ 2 F 2 L 2 (Ω) + C G 2 L 2 (∂Ω) ε ( 6
)
where C is the constant of the trace theorem, i.e.

φ L 2 (∂Ω) ≤ √ C φ H 1 (Ω) .
We also conclude using (L1h3) that the sequence (φ n ) n≥0 is uniformly bounded in L ∞ 0, T ; L 2 (Ω) ∩L 2 0, T ; H 1 (Ω) .

Then there exists a subsequence weakly converging in L 2 0, T ; H 1 (Ω) and weakly-in L ∞ 0, T ; L 2 (Ω) to the weak solution φ of the initial problem (5) satisfying the a priori estimate [START_REF] Szymańska | Mathematical modelling of heat shock protein synthesis in response to temperature change[END_REF]. Since the problem is linear, the difference between two solutions satisfies the same problem ( 5) with a vanishing source term F 0 and incoming flux G = 0. Using the a priori estimate (6), we conclude the uniqueness of the solution.

Lemma 2. Let Ω ⊂ R d be a smooth bounded domain and n be the outward normal to its boundary ∂Ω and

∂Ω -= {x ∈ ∂Ω | V (t, x) • n (x) < 0} and ∂Ω + = {x ∈ ∂Ω | V (t, x) • n (x) > 0} .
Assume that (L2h1) For any time t ∈ [0, T ], we have V (t, .) ∈ C 1 (Ω);

(L2h2) F ∈ L 2 (Ω T ) and G ∈ L 2 ∂Ω - T ; |V • n| dσdt .
Then, for any initial data φ 0 ∈ L 2 (Ω), there exists an unique weak solution φ in L ∞ 0, T ; L 2 (Ω) ∩L 2 (∂Ω T ; |V • n| dσdt) of the following hyperbolic equation

   ∂ t φ + ∇ • (φV) = F on Ω T φ = G on ∂Ω - T φ (0, x) = φ 0 (x) on Ω (7)
Proof. Since C 1 (Ω T ) is dense in L 2 (Ω T ), we can define the sequences F n ∈ C 0 Ω T and G n ∈ C 1 ∂Ω - T ), bounded in L 2 (Ω T ) and converging to F in L 2 (Ω T ) and G in L 2 ∂Ω - T , respectively. Then for any initial condition

φ 0 n ∈ C 1 Ω , there exists solution φ in C 1 (Ω T ) ∩ C 0 Ω T of    ∂ t φ n + V • ∇φ n + φ n ∇ • V = F n on Ω T φ n = G n on ∂Ω - T φ n (0, x) = φ 0 n (x) on Ω (8)
To prove existence, it is enough to see that the following function is a solution

φ n (t, x) = t τ(t,x) F n (s, χ (s, t, x)) e - t s ∇•V(θ,χ(θ,t,x))dθ +1 τ(t,x)=0 φ 0 n (χ (0, t, x)) e - t 0 ∇•V(θ,χ(θ,t,x))dθ +1 τ(t,x)>0 G n (τ (t, x) , χ (τ (t, x) , t, x)) e - t τ(t,x) ∇•V(θ,χ(θ,t,x))dθ
with the flow χ (s, t, x) defined as the solution of

∂ s χ (s, t, x) = V (s, χ (s, t, x)) , χ (t, t, x) = x.
The flow χ (s, t, x) is well-defined accordingly to the Picard-Lindelöf theorem and (L2h1). The time 0 ≤ τ (t, x) ≤ t is defined such that χ (τ (t, x) , t, x) ∈ ∂Ω -or τ (t, x) = 0. Then we define a sequence of initial data φ 0 n ∈ C 1 Ω , bounded in L 2 (Ω) and converging to φ 0 in L 2 (Ω). Multiplying Eq. ( 8) by the solution φ n we get the a priori estimate

∂ t φ n 2 L 2 (Ω) + φ n 2 L 2 (∂Ω + ;|V•n|dσ) ≤ F n 2 L 2 (Ω) + G n 2 L 2 (∂Ω -;|V•n|dσ) + 1 + ∇ • V 2 L ∞ (Ω) φ n 2 L 2 (Ω) . ( 9
)
Neglecting the second term of the LHS and using the Grönwall lemma, (L2h1) and (L2h2), we conclude that the sequence φ n is uniformly bounded in L ∞ 0, T ; L 2 (Ω) . Then, there exists a subsequence weakly converging in L ∞ 0, T ; L 2 (Ω) to the solution of [START_REF] Sturrock | Spatiotemporal modelling of the Hes1 and p53-Mdm2 intracellular signalling pathways[END_REF]. Since the problem is linear, the difference between two solutions satisfies the same problem [START_REF] Sturrock | Spatiotemporal modelling of the Hes1 and p53-Mdm2 intracellular signalling pathways[END_REF] with the vanishing source term F = 0 and incoming flux G = 0. We conclude therefore uniqueness of the solution φ of [START_REF] Sturrock | Spatiotemporal modelling of the Hes1 and p53-Mdm2 intracellular signalling pathways[END_REF]. Finally using the a priori estimate ( 9), we also get the estimation of the solution at the outgoing boundary ∂Ω + T .

Proof of Theorem 1. In the first step of the proof we use Lemma 1 and Lemma 2 to show the existence and uniqueness of the solution of the linearized model of ( 1)-( 2)-( 3)-( 4), more precisely assuming that the parameters are not function of the unknown functions. In the second step, we use contraction argument to show the convergence of the sequence of successive approximations and conclude the global existence of solution of the nonlinear system ( 1)-( 2)-( 3)-( 4).

Well-posedness of the linearized system: In this step we introduce successive approximations (m j i ), ( m j i ), (m j i ), (p j i ), ( p j i ) and (p j i ) of the solution of the linearized system, obtained by taking the source and decay terms with the previous approximations, and the boundary condition with the new approximation. More precisely, we take

k j 0i (t, x) = k 0i m j n ; p j n , κ j 0i (t, x) = κ 0i m j n ; p j n , j 0i (t, x) = 0i m j n ; p j n , λ j 0i (t, x) = λ 0i m j n ; p j n , a j 0i (t, x) = a 0i m j n ; p j n , α j i (t, x) = α i ξ; m j c ; p j c , k j i (t, x) = k i ξ; m j c ; p j c , κ j i (t, x) = κ i ξ; m j c ; p j c , j i (t, x) = i ξ; m j c ; p j c , λ j i (t, x) = λ i ξ; m j c ; p j c , b j i (t, x) = b i ξ; m j c ; p j c , β j i (t, x) = β i ξ; m j c ; p j c .
Assume that the previous approximations are in the corresponding Banach space

m j i ∈ L ∞ 0, T ; L 2 (N ) ∩ L 2 0, T ; H 1 (N ) , m j i ∈ L ∞ 0, T ; L 2 (C ) ∩ L 2 0, T ; H 1 (C ) , m j i ∈ L ∞ 0, T ; L 2 (C ) ∩ L 2 Γ c T ; v i |η • c| dσdt , p j i ∈ L ∞ 0, T ; L 2 (N ) ∩ L 2 0, T ; H 1 (N ) , p j i ∈ L ∞ 0, T ; L 2 (C ) ∩ L 2 0, T ; H 1 (C ) , p j i ∈ L ∞ 0, T ; L 2 (C ) ∩ L 2 Γ n T ; ν i |η • n| dσdt .
Note that accordingly to the assumptions of Theorem, the initial data, considered as constant functions of time, are in the corresponding Banach space. We define the next approximation m j+1 i by

       ∂ t m j+1 i -∇ • k j 0i ∇m j+1 i + εm j+1 i = a j 0i + ε -j 0i m j i on N T , n • k j 0i ∇m j+1 i + u i m j+1 i = 0 on Γ n T , m j+1 i (0, x) = m 0 i (x) on N . (10)
We are going to show that the RHS is in L 2 (N T ). We note that in particular the estimate of term j 0i m j i = 0i m j n , p j n m j i is not trivial. In order to this end we use the Lipschitz continuity of the function 0i (m n , p n ) m i (assumed by (T1h6)).

However, the assumption that the function 0i (m n , p n ) m i vanishes when the unknowns vanish, i.e. 0 = 0i (0, 0) 0 = 0i 0, given by (T1h5), is needed. Then we obtain

j 0i m j i = j 0i m j i -0i 0 ≤ A m j n -0; p j n -0 N ≤ A m j n ; p j n N .
Similarly, the term a j 0i = a 0i m j n , p j n is estimated adding the term when the unknowns vanish, i.e. a 0i (0, 0) = a 0i . We have

a j 0i 2 ≤ 2 a j 0i -a 0i 2 + 2 a 0i 2 ≤ 2A 2 m j n -0; p j n -0 2 N + 2 a 0i 2 ≤ 2A 2 m j n ; p j n 2 N + 2 a 0i 2 .
This yields the following estimate of the RHS

a j 0i + ε -j 0i m j i 2 ≤ 8A 2 m j n ; p j n 2 N + 4 a 0i 2 + 4ε 2 m j i 2 .
This is clearly integrable since the previous approximations m j i and p j i are in L 2 (N T ). It leads that the auxiliary problem [START_REF] Chaplain | Spatio-Temporal Modelling of Intracellular Signalling Pathways: Transcription Factors, Negative Feedback Systems and Oscillations[END_REF] satisfies Assumption (L1h3). We can easily check that Eq. ( 10) satisfies the other assumptions of Lemma 1. We conclude that there exists the unique global solution

m j+1 i ∈ L ∞ 0, T ; L 2 (N ) ∩ L 2 0, T ; H 1 (N ) .
Then, we consider the approximations given by

     ∂ t m j+1 i + ∇ • v i m j+1 i η = -j i + b j i m j i on C T , n • η v i m j+1 i = ξu i m j+1 i on Γ n T , m j+1 i (0, x) = m 0 i (x) on C . (11)
Similarly as for the auxiliary problem [START_REF] Chaplain | Spatio-Temporal Modelling of Intracellular Signalling Pathways: Transcription Factors, Negative Feedback Systems and Oscillations[END_REF], we show that the RHS is in L 2 (C T ). In addition since m j+1 i is in L 2 0, T ; H 1 (N ) and using the Trace theorem we have

m j+1 i ∈ L 2 Γ c T . It leads ξu i m j+1 i ∈ L 2 Γ c
T since u i is bounded (T1h3) and we conclude applying Lemma 2 that there exists the unique global solution

m j+1 i ∈ L ∞ 0, T ; L 2 (N ) ∩ L 2 Γ c T ; v i |η • c| dσdt .
Moreover, we consider

             ∂ t m j+1 i -∇ • k j i ∇ m j+1 i + ε m j+1 i = b j i m j i + ε -j i m j i on N T , -n • k j i ∇ m j+1 i = (1 -ξ) u i m j+1 i on Γ n T , c • k j i ∇ m j+1 i = c • η v i m j+1 i on Γ c T , m j+1 i (0, x) = m 0 i (x) on N . (12)
Using the properties of m j+1 i and m j+1 i , by Lemma 1 we conclude that there exists the unique global solution

m j+1 i ∈ L ∞ 0, T ; L 2 (C ) ∩ L 2 0, T ; H 1 (C ) .
Similarly we define the approximations by

             ∂ t p j+1 i -∇ • κ j i ∇ p j+1 i + ε p j+1 i = α j i + ε -λ j i -β j i p j i on N T , -n • κ j i ∇ p j+1 i + µ i p j+1 i = 0 on Γ n T , c • κ j i ∇ p j+1 i = 0 on Γ c T , p j+1 i (0, x) = p 0 i (x) on N . (13)
Using Lemma 1, we conclude that there exists the unique global solution

p j+1 i ∈ L ∞ 0, T ; L 2 (C ) ∩ L 2 0, T ; H 1 (C ) .
Then we consider

     ∂ t p j+1 i -∇ • ν i p j+1 i η = β j i p j i -λ j i p j i on C T , n • η ν i p j+1 i = 0 on Γ c T , p j+1 i (0, x) = p 0 i (x) on C . ( 14
)
and by Lemma 2 we conclude that there exists the unique global solution

p j+1 i ∈ L ∞ 0, T ; L 2 (N ) ∩ L 2 Γ n T ; ν i |η • n| dσdt .
Finally we consider

       ∂ t p j+1 i -∇ • κ j 0i ∇p j+1 i + εp j+1 i = ε -λ j 0i p j i on N T , -n • κ j 0i ∇p j+1 i = µ i p j+1 i + n • η ν i p j+1 i on Γ n T , p j+1 i (0, x) = p 0 i (x) on N . (15)
and by Lemma 1 we conclude that there exists the unique global solution

p j+1 i ∈ L ∞ 0, T ; L 2 (N ) ∩ L 2 0, T ; H 1 (N ) .
Global Well-posedness: This step is the proof of convergence of the sequences (m j i ), ( m j i ), (m j i ), (p j i ), ( p j i ) and (p j i ) when j tends to infinity. In the following, we use the sign function defined in a standard way

sgn (φ) =    1 , if φ > 0 , 0 , if φ = 0 , -1 , if φ < 0 . Let (m c ; p c ) L 1 (C ) = N i=1 m i L 1 (C ) + m i L 1 (C ) + p i L 1 (C ) + p i L 1 (C ) , (m n ; p n ) L 1 (N ) = N i=1 m i L 1 (N ) + p i L 1 (N ) and (m; p) L 1 (Ω) = (m n ; p n ) L 1 (N ) + (m c ; p c ) L 1 (C ) .
In addition we denote

δm j+1 i = m j+1 i -m j i , δ m j+1 i = m j+1 i -m j i , δm j+1 i = m j+1 i -m j i , δp j+1 i = p j+1 i -p j i , δ p j+1 i = p j+1 i -p j i , δp j+1 i = p j+1 i -p j i .
We estimate the above differences by multiplying each equation by the sign of the difference and integrating over the space domain. More precisely, we multiply Eq. ( 10) by sgn δm j+1 i and by Lipschitz-continuity of the RHS (T1h6), we obtain

∂ t δm j+1 i L 1 (N ) + Γ n u i δm j+1 i dσ ≤ (2A + ε) δm j n ; δp j n L 1 (N ) . (16)
Next we multiply [START_REF] Chaplain | Hopf bifurcation in a gene regulatory network model: Molecular movement causes oscillations[END_REF] by sgn δm j+1 i and we have

∂ t δm j+1 i L 1 (C ) + Γ c c • η v i δm j+1 i dσ ≤ 2A δm j c ; δp j c L 1 (C ) + Γ n ξu i δm j+1 i dσ . (17)
Multiplying ( 12) by sgn δ m j+1 i yields

∂ t δ m j+1 i L 1 (C ) ≤ (2A + ε) δm j c ; δp j c L 1 (C ) + Γ n (1 -ξ) u i δm j+1 i dσ + Γ c c • η v i δm j+1 i dσ . (18)
Summing up these estimates and using the boundaries balance we obtain

∂ t δm j+1 i L 1 (N ) + δm j+1 i L 1 (C ) + δ m j+1 i L 1 (C ) ≤ (4A + ε) δm j ; δp j L 1 (Ω) . (19)
Similarly we estimate the time evolution of the differences δp j+1 i , δp j+1 i , δ p j+1 i . Then summing with respect to 1

≤ i ≤ N yields ∂ t δm j+1 ; δp j+1 L 1 (Ω) ≤ 2N (4A + ε) δm j ; δp j L 1 (Ω) . (20)
Finally by induction we obtain

δm j+1 ; δp j+1 L 1 (Ω) (t) ≤ (2N (4A + ε) t) j j! δm 1 ; δp 1 L ∞ (0,T;L 1 (Ω)) . (21)
Since the RHS tends to zero when j goes to infinity, the sequences (m j i ), ( m j i ), (m j i ), (p j i ), ( p j i ) and (p j i ) converge. Then the global well-posedness of the nonlinear system (1)-( 2)-( 3)-( 4) follows.

We now the biological relevance of the model by the following positivity and conservation results: Proposition 1. [Global mass laws] The global concentration of each species satisfy the following equation

∂ t Ω M i dx = N a 0 i -0 i M i dx - C i M i dx and ∂ t Ω P i dx = - N λ 0 i P i dx + C (α i -λ i P i ) dx (22)
with the total concentration of mRNA defined by

M i (t, x) = m i (t, x) , if x ∈ N m i (t, x) + m i (t, x), if x ∈ C (23)
and the total concentration of proteins defined by

P i (t, x) = p i (t, x) , if x ∈ N p i (t, x) + p i (t, x), if x ∈ C . (24)
Proof. The result follows by direct estimates by the boundaries (3) and ( 4) balance.

The conservation law presented in Proposition 1 shows that the only processes that change of the global concentration of mRNA or protein are the transcription, the translation and the decay of molecules. Assuming the space homogeneous case, the model can be reduce to a ODEs system, see [START_REF] Agrawal | Computational models of the notch network elucidate mechanisms of context-dependent signaling[END_REF][START_REF] Momiji | Dissecting the dynamics of the hes1 genetic oscillator[END_REF][START_REF] Monk | Oscillatory expression of hes1, p53, and nf-k b driven by transcriptional time delays[END_REF]. However, the space distributions of the molecules are not homogeneous since the productions are located in a part of the cell. We performed the numerical simulations with parameters values relevant from biological point of view, that illustrate the heterogeneity of the space distribution, see [START_REF] Szymańska | Mathematical modeling of the intracellular protein dynamics: the importance of active transport along microtubules[END_REF].

Proposition 2.

[Positivity] Let the assumptions of Theorem 1 hold. Moreover, if for any (φ, ψ) ∈ R N 2 , 0i (φ, ψ) ≥ 0, λ 0i (φ, ψ) ≥ 0, and a 0i (φ, ψ) ≥ 0 ,

and for any ξ ∈

[0, 1], (φ, ψ) ∈ R 2N 2 i (ξ, φ, ψ) ≥ 0, λ i (ξ, φ, ψ) ≥ 0, α i (ξ, φ, ψ) ≥ 0, b i (ξ, φ, ψ) ≥ 0 and β i (ξ, φ, ψ) ≥ 0 , (26)
then for all initial data, almost everywhere nonnegative, the solution is almost everywhere nonnegative.

Proof. We use the positive and negative part function 2φ ± = |φ| ± φ. First we multiply (2a) by the negative part of the solution m - i and we integrate over the space domain. It follows

∂ t m - i 2 L 2 (N ) + 2 m - i 2 H 1 (N ) + 2 m - i 2 L 2 (Γ n ;u i dσ) ≤ -2 N a 0i m - i dx .
By the construction the negative part function is nonnegative. Since the transcription rate is nonnegative, it leads that the right hand side is nonpositive

∂ t m - i 2 L 2 (N ) ≤ 0 .
Since the initial data are nonnegative, we conclude that m i is nonnegative almost everywhere. In addition, the negative part of the outgoing flux vanishes m - i L 2 (Γ n ;u i dσ) = 0 .

Then we multiply (1b) by the negative part of the solution m i -and we integrate over the space domain. We have

∂ t m i -2 L 2 (C ) + m i -2 L 2 (Γ c ;v i |η•c|dσ) + m i -2 L 2 (C ;b i dx) ≤ ∇ • (v i η) L ∞ (C ) m i -2 L 2 (C ) + Γ n ξu i m - i m i -dσ - Γ n ξu i m + i m i -dσ . ( 28
)
Using the Cauchy-Schwarz inequality we obtain

Γ n ξu i m - i m i -dσ ≤ m - i L 2 (Γ n ;u i dσ) m i - L 2 (Γ n ;u i dσ) ,
which vanishes accordingly to (27). The last term of the RHS of (28) is clearly nonpositive, thus

∂ t m i -2 L 2 (C ) ≤ ∇ • (v i η) L ∞ (C ) m i -2 L 2 (C ) .
Using the Grönwall lemma, and since the initial data are nonnegative, we conclude that m i are nonnegative almost everywhere, i.e. m - i L 2 (C ) = 0 .

In addition, using the estimate (28), we conclude that the negative part of the outgoing flux and the negative part of the exchange term vanish, i.e. m i

- L 2 (Γ c ;v i |η•c|dσ) = m i - L 2 (C ;b i dx) = 0 . (30)
Finally, we multiply (1a) by m - i and integrate over the space domain. We obtain

∂ t m - i 2 L 2 (C ) + 2 m - i 2 H 1 (C ) ≤ -2 C b i m i m - i dx -2 Γ n (1 -ξ) u i m i m - i dσ -2 Γ c m i m - i v i |η • c| dσ .
We show that the positive part of the RHS vanishes using Cauchy-Schwarz inequality and the previous estimates (27),( 29) and (30). It follows ∂ t m - i 2 L 2 (C ) ≤ 0 , and since the initial data are nonnegative, we conclude that m i is nonnegative almost everywhere, i.e. m - i L 2 (C ) = 0. Similarly we show that p i and p i are nonnegative almost everywhere.

The non-negativity of the concentrations is clearly essential in the term of biological interpretation. The fact that the solution of the model can not be negative allowed to simulate and analyse extreme situation for example when some parameters tend to vanish or become very high.

Conclusion

In this paper we have presented the preliminary mathematical analysis of the non-linear, multidimensional and multi-species model proposed by Szymańska at al. in 2014 [START_REF] Szymańska | Mathematical modeling of the intracellular protein dynamics: the importance of active transport along microtubules[END_REF], including active transport along the microtubules. The main properties such as well-posedness, positivity and global evolution are highlighted. However there are still many open mathematical problems including the analysis with less regular parameters, in particular transcription rate or speed along the microtubules. This questions seem relevant in the biological point of view and required a trickier analysis. In addition, the existence of periodic solutions that has been observed numerically [START_REF] Szymańska | Mathematical modeling of the intracellular protein dynamics: the importance of active transport along microtubules[END_REF] are still an open problem and are an important issus on the understanding of cell homeostasis and many disease like cancer.

Figure 1 :

 1 Figure 1: Scheme of the space decomposition of the cell. Designation: Ω -eukaryotic cell; C -cytoplasm; N -nucleus; Γ c -cell membrane;Γ n -nuclear envelope; ξ -concentration of microtubules; m i -mRNA molecules in the nucleus; p i -protein molecules in the nucleus; m i -mRNA molecules in the cytoplasm that are linked to the microtubules; m i -mRNA molecules in the cytoplasm that are free, i.e. not linked to the microtubules; p i -protein molecules in the cytoplasm that are free, i.e. not linked to the microtubules, i = 1, . . . , N;

Acknowledgements M.L. acknowledges a support from the University of KwaZulu-Natal Research Found (South Africa). Z.S. acknowledges a support from the National Science Centre Poland grant 2011/01/D/ST1/04133.


Sven Verdoolaege

Keywords: structures, function summary, dead code elimination, live-range reordering compilation polyhédrique, tableau de structures, sommaire de fonction, suppression du code mort, réorganisation de l'ordre

This report describes various changes to the polyhedral model extractor pet and the automatic parallelizer PPCG, including support for generic statements, arrays of structures, function summaries, dead code elimination and live-range reordering as well as support for annotations allowing the user to explicitly describe relations among program variables, kills and the absence of loop-carried dependences. Most of these changes are instrumental in supporting pencil in PPCG. The report also describes how synchronization is introduced by PPCG.

Key-words: polyhedral compilation, array of

Introduction

This document describes some of the changes that have been applied to pet (Verdoolaege and Grosser 2012) and PPCG (Verdoolaege et al. 2013b) since their initial descriptions, with a special focus on the changes that were needed to support pencil (Baghdadi et al. 2013b). Most of these changes have been available for some time already, but had never been described in detail before. The most recent changes have been made available in ppcg-0.03-202-gc63f95d. There are also some changes that are not covered by this document:

• OpenCL backend in PPCG This backend was initially contributed by Riyadh Baghdadi and has been further improved by Sven van Haastregt.

• support for while, break and continue

The direct support for these constructs in pet is mainly useful for constructing weakly dynamic polyhedral process networks [START_REF] Stefanov | Converting Weakly Dynamic Programs to Equivalent Process Network Specifications[END_REF]) through run-time dependent dataflow analysis (Verdoolaege et al. 2013a). These constructs are not yet supported in PPCG. Instead, they are encapsulated inside generic statements (see Section 3.2).

• switch to schedule trees pet and PPCG now use the schedule trees of [START_REF] Verdoolaege | Schedule Trees[END_REF] to represent schedules.

• multiple SCoPs pet now provides an interface for extracting multiple SCoPs from the same input file.

Terminology

An expression is called static-affine if it is a quasi-affine expression involving only symbolic constants and loop iterators.

The statement instance set is a set that contains an element for each dynamic execution of a statement in the analyzed program fragment.

Changes to pet

This section describes some of the changes that have been applied to pet since its initial description (Verdoolaege and Grosser 2012), with a special focus on the changes that were needed to support pencil (Baghdadi et al. 2013b).

Accesses

For each access in the analyzed fragment, the original version of pet (Verdoolaege and Grosser 2012) would only keep track of the following pieces of information.

• an isl_map representing the access relation RT n°457

• a bit indicating whether this access is a read

• a bit indicating whether this access is a write In order to support various extensions, the version of pet described by this document keeps track of the following pieces of information instead.

• an isl_multi_pw_aff representing the original index expression

• a reference identifier

• three isl_union_maps representing the may-reads, the may-writes and the must-writes

• a bit indicating whether this access is a read

• a bit indicating whether this access is a write

• a bit indicating whether this access is a kill

The reference identifier is needed to distinguish multiple accesses in the same statement. For example, support for live-range reordering (Baghdadi et al. 2013a) requires PPCG to keep track of which references cause any given dependence between statement instances. The reference identifiers can be used here to identify the references involved.

The index expression is mainly used for printing accesses in the transformed program. It is first rewritten and simplified in terms of the loop iterators of the transformed program and then optionally mapped to a shared or private memory tile (Verdoolaege et al. 2013b, Section 7). One advantage of an isl_multi_pw_aff representation is that the different index expressions can be processed individually. In particular, a disjunction appearing in one of them will not affect the others, whereas in an isl_map representation, the disjunction happens at the outer level such that it affects all indices. The main reason for introducing an explicit index expression, however, is that, as explained next, an isl_map no longer suffices to represent the access relations(s) and is therefore no longer available.

In the original version of pet, a slice of an array passed to a function was either completely read or completely written. With the support for function summaries (Section 3.6), pet can now construct a more accurate access relation by analyzing the function that is being called (or its summary function). This means that a distinction needs to be made between the elements that are read and those that are or may be written within a single array access. Furthermore, because the function may accept an array of structures (Section 3.3), the accessed elements may not all live in the same space. Instead of a single isl_map, three isl_union_maps are therefore required to express the access relations.

The kill bit is used to mark "accesses" that appear in kill statements (Section 3.8). These accesses do not actually access any data, but merely describe the set of elements that are being killed by the statement.

Generic statements

Whereas, originally, a statement in pet could only be an expression statement, it can now be any statement. In particular, the body of the statement is a RT n°457 for ( int x = 0; x < n ; ++ x ) { A : s = f (); B : while ( P (x , s )) s = g ( s ); C : h ( s ); } Listing 1: A program with a while loop for ( int i = 0; i < 100; ++ i ) { for ( int j = 0; j < 100

; ++ j ) { if ( A [ j ]) { x += 0; continue ; } x = 1; } }
Listing 2: A program with a continue pet_tree representing an AST. The support for such generic statements allows pet to (optionally) encapsulate any dynamic control inside its smallest containing statement. For example, the program in Listing 1 is considered to consist of three statements, with the second statement encapsulating the while loop and its body. The instance set is { A(x) : 0 ≤ x < n; B(x) : 0 ≤ x < n; C(x) : 0 ≤ x < n }. The program in Listing 2 is considered to only have a single statement, encapsulating the inner for loop. This inner loop has dynamic control because of the continue statement governed by a condition that is not static-affine. The statement instance set for this program is { S(i) : 0 ≤ i < 100 }.

The generic statement AST in the body of a polyhedral statement is not further analyzed by pet but is instead considered as one indivisible entity. pet does take into account that there may be some dynamic control inside such statements and therefore treats all writes inside such a statement as a maywrite. The dataflow analysis performed by PPCG is applied to these generic statements and their may-writes and it may therefore find dataflow from other statement instances to a given instance of a statement even if in reality the dataflow happens inside the statement. For example, the value of t read inside the if-statement in Listing 3 is clearly written inside the same instance of the if-statement. In practice, however, the result of the dataflow analysis will have the read depend on every previous iteration, as well as on the last write to t before the loop. This issue could be resolved by moving the dataflow analysis inside pet so that it could perform the dataflow analysis inside the generic statement before encapsulating it. As a work-around, the spurious dataflow can also be removed by inserting kills (Section 3.8).

RT n°457

for ( int i = 0; i < 100

; ++ i ) { if ( A [ i ] > 0) { t = A [ i ]; B [ i ] = t ; } } Listing 3: Dataflow inside an encapsulating statement struct s { struct { int a [10]; } f [10]; }; void f () { struct s s ;
for ( int i = 0; i < 10; ++ i ) for ( int j = 0; j < 10; ++ j ) s . f [ i ]. a [ j ] = i * j ; } Listing 4: Accesses to structure fields

Structures

Accesses to fields of structures are represented using wrapped spaces. In particular, the data space is a wrapped map with the domain space referring to the accessed array and the range space identifying the field of the structure. The wrapped space is itself given a name that is composed of the array name and the field name. In case of nested field accesses, the wrapped space corresponding to the outer field access is taken as the domain of the wrapped space of the inner field access. For example, the nested field access in Listing 4 is represented as

{ S(i, j) → s f a(s f(s() → f(i)) → a(j)) }.
When collecting all may-reads, may-writes or must-writes in a pet_scop, pet returns access relations to the innermost fields. For example, for the program in Listing 5, the index expressions are

{ S[] -> a_b[a[] -> b[]] }, { T[] -> a_a[a[] -> a[]] }, { U[i] -> b[(i)] } and { U[i] -> a[] }.
However, the collected read access relation is

{ U[i] -> a_b[a[] -> b[]]; U[i] -> a_a[a[] -> a[]] } (1)
{ S [] -> a_b [ a [] -> b []]; T [] -> a_a [ a [] -> a []]; U [ i ] -> b_a [ b [( i )] -> a []]; U [ i ] -> b_b [ b [( i )] -> b []] } (2)

Two-Phase Approach

The extraction of a polyhedral model now follows a two-phase approach. In the first phase, the clang AST is analyzed in a bottom-up fashion to find a region that is amenable for polyhedral analysis. At the same time, this part of the clang AST is also converted into an internal pet_tree representation. At this stage, pet has not yet detected any symbolic constants or outer loop iterators (because the tree is constructed bottom-up). All accesses are therefore represented as data dependent accesses (Verdoolaege and Grosser 2012, Section 4.1). For example, the read from D in Listing 6 is initially represented using the index expression

{ [[] -> [i0]] -> D[((i0) : i0 >= 0)] }, (3)
where i0 represents the (currently) unknown result of the expression i + m. Inside the representation of this expression, the two arguments are represented as accesses to the variables i and m with index expressions

{ [] -> i[] } and { [] -> m[] }.
In the second phase, a pet_scop is extracted from the pet_tree constructed in the first phase. During a top-down traversal of this pet_tree, pet keeps track of a context containing known variable assignments and the domain for affine expressions. In particular, the domain collects the outer loop iterators, some of which may be virtual, along with their known bounds. Virtual iterators occur in infinite loops or loops with an unsigned loop iterator that may wrap. The

RT n°457

for ( int i = 0; i < n ; i ++) { D

[ i ] = D [ i + m ]; }
Listing 6: Simple loop known variables assignments include symbolic constants, outer loop iterators and scalar variables that have been assigned a quasi-affine expression and that do not change within the current subtree. Each variable assignment maps an identifier to a piecewise quasi-affine expression. The domain of this expression is a prefix of the domain of the context. For example, when pet reaches the body of the loop in Listing 6, it has the following assignments:

i: [i] -> [(i)] [i] -> [(i)] m: [m] -> [] -> [(m)] n: [n] -> [] -> [(n)] (4)
Note that m and n are declared as parameters, while i is the first iterator in the context domain.

During the construction of a pet_scop from a pet_tree, each expression is evaluated in the context, meaning that the following steps are applied.

• Insert context domain.

• Plug in known values for read accesses.

• Try and plug in index expression arguments by converting them to affine expressions.

• Plug in function summaries. See Section 3.6.

Considering once more the loop body in Listing 6, inserting the context domain

{ [i] } into the index expression of the read (3) results in the updated index expression { [[i] -> [i0]] -> D[((i0) : i0 >= 0)] }. (5)
In the same way, the index expressions of the nested accesses are updated to

{ [i] -> i[] } and { [i] -> m[] }.
The latter two index expressions access scalars for which a known value is available in the context, so they are further updated to

{ [i] -> [(i)] } and [m] -> { [i] -> [(m)] }.
At this point, pet tries to convert the + operation in the argument of the read of D into a quasiaffine expression. This operation succeeds here and results in the expression

[m] -> { [i] -> [(i + m)] }.
This result can then be plugged into the index expression (5) resulting in

[m] -> { [i] -> D[((i + m) : i + m >= 0)] }. (6)
Note that as explained next, pet also takes into account signed integer overflow, so the actual result is slightly more complicated. Next to the context in which expressions are evaluated and which is used during the construction of a pet_scop, a pet_scop itself also keeps track of a set called "context". This context collects constraints on the symbolic constants

RT n°457 void foo ( int n , int m , int S , int D [ const restrict static S ]) { __pencil_assume ( m > n ); for ( int i = 0; i < n ; i ++) { D [ i ] = D [ i + m ]; } }
Listing 7: Explicit assumptions that can be assumed to hold because the code represented by the pet_scop could not be executed otherwise. These constraints are derived from the bounds on the integer types used to declare the variables corresponding to the symbolic constants, the absence of signed integer overflow and the absence of negative array sizes and array indices. Note that while negative array indices are allowed in C99, they are not allowed by pet because they may lead to aliasing. During the construction of a pet_scop, the context of the pet_scop does not only refer to the symbolic constants, but also to the outer loop iterators. Each of these outer loop iterators is projected out on the way back up the pet_tree, from innermost to outermost. In particular, when leaving a loop with loop iterator j and with outer loop iterators i, if dom(i, j) represents the known constraints on the outer and current loop iterators and valid(i, j) represents the current context inside the loop construct, then the context outside the loop is ∀j : dom(i, j) =⇒ valid(i, j).

Explicit Assumptions

As explained in the previous section, pet derives necessary conditions on the symbolic constants from the bounds of the corresponding integer types and the absence of signed integer overflow and negative array sizes and array indices. The user is allowed to add additional constraints through a "call" to __pencil_assume. The argument of this call is a boolean condition that is guaranteed to hold at that position by the user. While this boolean condition may have any form, only those that are quasi-affine expressions in the symbolic constants and outer loop iterators are effectively taken into account. Consider, for example, the program in Listing 7. Due to the explicit assumption, the final context satisfies the constraint m > n. In principle, this constraint could be used to prove that there are no loop-carried dependences in the program, but this context information is currently not exploited by PPCG during its dataflow analysis.

Function Summaries

Whenever pet comes across a function call during the first phase of the analysis, it checks if the body of the function is available. If so, a pet_scop is constructed for the entire function body. A failure to extract such a pet_scop is currently considered to be an error. For pencil input, this should not be a problem since the entire source file is supposed to conform to the pencil restrictions. Before

RT n°457 int f ( int i ); int maybe (); struct s { int a ; }; void set_odd ( int n , struct s A [ static n ]) { int t ; S : t = f (0); for ( int i = 1; i < n ; i += 2) if ( maybe ()) T : A [ i ]. a = t ; } void foo ( int n , struct s B [ static 2 * n ]) { # pragma scop set_odd (2 * n , B ); # pragma endscop }
Listing 8: Program with function call extracting a pet_scop, pet first analyzes the arguments of the function and groups them into three classes

• integer arguments, i.e., those with an integer type,

• array arguments, i.e., those that may access data from the caller, and

• other arguments.
The pet_scop is then extracted with the integer arguments as the only symbolic constants. During this extraction, the autodetect option is turned off to ensure that the entire body is extracted. After the extraction, a function summary is constructed that collects summary information about the function. In particular, for each argument, the function summary keeps track of its type and some type-dependent information. For integer arguments, this extra piece of information is the identifier. For array arguments, three isl_union_maps are stored corresponding to the may-reads, the may-writes and the must-writes. The domains of these relations are formed by a space with a coordinate for each integer argument. The access relations are obtained by collecting the sets of accessed elements from the pet_scop extracted from the function body, distributing them over the array arguments and then replacing the symbolic constants by the corresponding local domain coordinates. Consider, for example, the program in Listing 8. While analyzing the body of foo, pet runs across a call to set_odd. This function has two arguments,

RT n°457

n and A. The first is recognized as an integer argument, with identifier n, and the second argument is recognized as an array argument. The entire set_odd function body is therefore analyzed with n as the only symbolic constant. Assuming the option to encapsulate dynamic control is turned on, the body of the for loop is considered as a monolithic statement performing a may-write. The may-write access relation for the entire body is therefore

[ n ] -> { S [] -> t []; T [ i ] -> A_a [ A [ i ] -> a []] : exists ( e0 : 2 e0 = -1 + i and 1 <= i <= -1 + n ) } (7)
The must-write access relation is

[n] -> { S[] -> t[] } (8)
while the may-read access relation is

[ n ] -> { T [ i ] -> t []: exists ( e0 : 2 e0 = -1 + i and 1 <= i <= -1 + n ) } (9)
Projecting out the statement instances of the body (which are not relevant to the caller) from the may-write access relation, results in

[ n ] -> { t []; A_a [ A [ i ] -> a []] : exists ( e0 : 2 e0 = -1 + i and 1 <= i <= -1 + n ) } (10)
The array arguments are then considered in turn and for each of them the extent is determined in the context of the same symbolic constants. For the second argument, this extent is

[n] -> { A[i0] : i0 >= 0 and i0 <= -1 + n }. (11)
Since the access relations are all given in terms of the innermost fields, this extent first needs to be mapped to those innermost fields as well, resulting in

[n] -> { A_a[A[i0] -> a[]] : i0 >= 0 and i0 <= -1 + n }. ( 12
)
Intersecting this extent with the may-write access (10) results in the may-writes for this array argument

[ n ] -> { A_a [ A [ i ] -> a []] : exists ( e0 : 2 e0 = -1 + i and 1 <= i <= -1 + n ) } (13)
The symbolic constants in this set are then replaced by domain coordinates resulting in

{ [ n ] -> A_a [ A [ i ] -> a []] : exists ( e0 : 2 e0 = -1 + i and 1 <= i <= -1 + n ) } (14) RT n°457
pencil support in pet and PPCG

The may-reads and the must-writes for this array argument are empty.

During the construction of a pet_scop from the pet_tree representing the body of foo, the call expression is evaluated in its context, which performs the steps detailed in Section 3.4 and at the end plugs in the corresponding function summary, if any. Plugging in a function summary involves the following steps. The integer arguments in the function call are converted to affine expressions in the current context, if possible. Otherwise, additional arguments are added to the context domain, resulting in data dependent expressions (Verdoolaege and Grosser 2012, Section 4.1). For each array argument in the function summary, if the corresponding argument in the function call is effectively an access expression, then integer arguments in the function call are plugged into the access relations in the function summary and combined with the index expression of the array access in the function call to produce access relations for this array access.

Returning to the example of Listing 8, there is only a single integer argument. The actual argument in the function call can in this case be represented as an affine expression in terms of the context domain

[n] -> { [] -> [(2n)] }. ( 15
)
Since there is only one integer argument, collecting all integer arguments results in the same affine expression here. Plugging in this parameter binding in the may-write access relation ( 14), results in

[ n ] -> { [] -> A_a [ A [ i ] -> a []] : exists ( e0 : 2 e0 = -1 + i and 1 <= i <= -1 + 2 * n ) } (16)
Finally, this access relation is combined with the index expression { [] -> B[] }, resulting in the final may-write access relation

[ n ] -> { [] -> B_a [ B [ i ] -> a []] : exists ( e0 : 2 e0 = -1 + i and 1 <= i <= -1 + 2 * n ) } (17)
Combining the index expression of the caller with the access relations in the function summary involves dropping the name of the outer array in the access relation and pasting the index expressions of this outer array to those of the inner array of the index expression. In the example above, the dropped outer array name is A and its single index expression is added to the initially zero index expressions of B.

Summary Functions

In some cases, the access relations derived by pet from a called function may be too inaccurate. In the extreme case, no code may be available for the function such that pet can only assume that every element of the passed arrays is accessed. In order to obtain more accurate access relations, the user may tell pet to derive the access relations not from the actual function body, but from some other function with the same signature. Such a function is called a summary function and is treated in exactly the same way as in Section 3.6 for

RT n°457 int f ( int i ); int maybe (); struct s { int a ; }; void set_odd_summary ( int n , struct s A [ static n ]) { for ( int i = 1; i < n ; i += 2) if ( maybe ()) A [ i ]. a = 0; } __attribute__ (( pencil_access ( set_odd_summary ))) void set_odd ( int n , struct s A [ static n ]); void set_odd ( int n , struct s A [ static n ]) { for ( int i = 0; i < n ; ++ i ) A [2 * f ( i ) + 1]. a = i ; } void foo ( int n , struct s B [ static 2 * n ]) { # pragma scop set_odd (2 * n , B ); # pragma endscop }
Listing 9: Summary function the purpose of deriving access relations of the calling statement. In particular, pet only collects the accesses in the function body and does not care about the order in which they are executed. Listing 9 illustrates the use of summary functions. Since pet is unable to analyze the index expression in the body of the set_odd function, it would derive the may-write access relation

{ [ n ] -> A_a [ A [ i ] -> a []] : 1 <= i <= -1 + n } (18)
for its function summary, missing the fact that only the odd elements could possibly be accessed. By writing a summary function set_odd_summary where this fact is more explicitly available and by marking it as a summary function for set_odd through the pencil_access attribute, the user can make pet derive the more accurate may-write access relation (14) of the previous section.

RT n°457 void foo ( int n , int A [ n ][ n ] , int B [ n ][ n ]) {
for ( int i = 0; i < n ; ++ i ) # pragma pencil independent for ( int j = 0; j < n ; ++ j )

B [ i ][ A [ i ][ j ]] = i + j ; }
Listing 10: Explicitly marked absence of dependence

Kills

A kill statement in a pet_scop represents the fact that no dataflow on the killed data elements can pass through any instance of the statement. This information can be used during dataflow analysis to stop the search for potential sources on data elements killed by the statement. Whenever pet comes across a variable declaration, two kill statements that kill the entire array are introduced, one at the location of the variable declaration and one at the end of the block that contains the variable declaration. Note that, currently, only single variable declarations are supported by pet. Any variable declared inside a statement encapsulating dynamic control is killed before and after the statement.

A user can also introduce explicit kills by adding a "call" to __pencil_kill. A kill statement is introduced for each argument passed to __pencil_kill. Section 4.4 describes an example of the use of __pencil_kill in Listing 12.

Absence of Dependences

Whenever there is control or an index expression that is not static-affine, the may-read and may-write access relations may be strict overapproximations of the actual accesses. These may therefore result in dependences being detected that do not occur in practice. The user can exclude such dependences by introducing a #pragma pencil independent. When placed in front of a loop, this pragma means that, assuming that a variable that is declared inside the loop is considered private to any given iteration, the iterations of the loop may be freely reordered with respect to each other, including reorderings that result in (partial) overlaps of distinct iterations. In particular, the user asserts through this pragma that no dependences need to be introduced to prevent such reorderings. For example, in the program of Listing 10, the user may have additional information about the A array ensuring that each row has all distinct entries. This means that each iteration of the inner loop accesses a different element of the B array such that the iterations of this loop may be freely reordered.

Since pet itself does not perform any form of dataflow analysis, it needs to propagate the independence information to its users. In particular, for each #pragma pencil independent, it keeps track of all pairs of statement instances for which no dependence needs to be introduced as well as all variables that are local to the body of the loop. This information is currently only collected for for-loops with static-affine loop bounds and a constant increment. The representation is fairly inefficient because it is quadratic in the number of statements

RT n°457

inside the loop body. The support for these pragmas predates the introduction of schedule trees. Now that pet keeps track of schedule trees, it should be possible to come up with a more efficient representation. For example, a mark node could be introduced on top of the band node corresponding to the loop decorated with the pragma. Ultimately, however, it may be best for pet to perform dataflow analysis itself such that this information no longer needs to be propagated to its users.

Changes to PPCG

This section describes some of the changes that have been applied to PPCG since its initial description (Verdoolaege et al. 2013b), with a special focus on the changes that were needed to support pencil (Baghdadi et al. 2013b).

Dead Code Elimination

Dead code elimination is performed by constructing the set of statement instances that may need to be executed. Initially, this is the set of statement instances that may be the last to write to an array element. In particular, a dataflow analysis is performed where the "sinks" are the must-writes and the kills, while the "potential sources" are the may-writes. All the may-writes that appear as the source in the result of this dataflow analysis are definitely killed by the corresponding must-write or kill and therefore cannot possibly be the last to write to the corresponding array element. The statement instances performing the remaining set may-writes then form the initial set of statement instances that need to be executed. This set is extended with all the statement instances that perform a function call since pet currently does keep track of "pure" functions and so PPCG can only assume that each function call needs to be performed.

The statement instances in the initial set may depend on other statement instances through potential flow dependences, which therefore need to be executed as well. The dataflow dependence relation is therefore applied to the current set and the result is added to this set. This process then needs to be repeated until it stops producing additional elements. However, a naive application of this scheme could lead to an infinite loop in case of parametric bounds on the original instance set since each iteration of the scheme would consider ever increasing values of the symbolic constants. In order to avoid this trap, a form of widening [START_REF] Cousot | Comparing the Galois Connection and Widening/Narrowing Approaches to Abstract Interpretation[END_REF]) is applied. In particular, after each addition of new statement instances, the current set is replaced by its integer affine hull (Verdoolaege 2010) and then intersected with the complete set of statement instances. Each iteration therefore increases the dimension of the instance set of at least one statement. Since the number of statements is finite and since the dimension of each of the corresponding instance sets is finite, this process is guaranteed to terminate.

Structures

Support in PPCG is fairly basic at the moment. The fields of a structure are assumed not to alias with each other (or indeed anything else), so that dataflow RT n°457 analysis does not need to be adjusted to support structures. If some element of a structure is accessed on the device, then the entire outer array containing the structure is copied as a whole to the device and back from the device (if needed). Arrays of structures are currently not considered for mapping to shared or private memory.

Kills

The kills introduced in Section 3.8 are used by PPCG in three steps,

• the dataflow computation,

• the computation of the set of live-out accesses, and

• the computation of the may-persist set.

The first two are explained in Section 4.4. The third is explained here.

The may-persist set is used during the computation of which elements to copy in and/or out of the device. PPCG first computes the set of data elements that need to be copied out. This set is computed from the may-writes that occur within the kernel(s) for which the copying needs to be performed. Some of these may-writes may however only be needed within the given subtree of the schedule tree (which is entirely mapped to the device) and therefore may not need to be copied out. In particular, PPCG removes those may-writes that are only involved in local dataflow and that are not live-out. In the end, the copy-out set is approximated by the entire arrays for which at least one element needs to be copied out. The copy-in set is similarly derived from the may-reads. However, there may be elements that are copied out that are not necessarily written inside the relevant part of the schedule tree. That is, there may be may-writes that are not also must-writes or there may be elements in the approximation of the copy-out set that may not be possibly written at all. In order to preserve their values with respect to the host, these arrays therefore need to be added to the copy-in set. This is where the may-persist set comes in. The host code may not care about some of these elements and the copy-in would be a waste of time and energy. Instead of adding all elements that are copied out but that may not have been written to the copy-in set, PPCG only adds those that are also elements of the may-persist set or that may be involved in dataflow that crosses the schedule subtree at hand. The may-persist set then is computed from the extent of all arrays that are not locally declared inside the program fragment, mapped to the innermost fields. In particular, the may-persist set consists of all elements of those arrays except those elements that are definitely written or killed somewhere in the program fragment.

Note that the way the copy-in set is updated to include those elements of the copy-out set that may not have been written is similar to the way approximations are handled by Alias et al. 2011.

Live-range Reordering

Baghdadi et al. (2013a) describe a relaxed permutability criterion that can be used to decide whether a schedule band can be tiled. However, the criterion RT n°457

for ( i = 0; i < n ; i ++) for ( j = 0; j < n ; j ++) { S1 : The same criterion can however also be used during the construction of schedule bands and this is how it is implemented in PPCG.

t = A [ i ][ j ] * y_1 [ j ] ; S2 : x1 [ i ] = x1 [ i ] + t ; } R0 R1 R2 R3 R4 R5 Listing
The main idea of Baghdadi et al. (2013a) is that dependences that are only introduced to prevent live-ranges from overlapping may be ignored during tiling if the corresponding live-ranges are local to the schedule band that is being tiled. That is, if the source and the sink of a live-range are mapped to the same iteration of the schedule band, then this live-range may be freely reordered with respect to other live-ranges on the same array elements by the schedule band and therefore does not prevent tiling. PPCG exploits this idea during the very construction of permutable bands.

Since a dependence between two statement instances may be caused by several pairs of accesses inside those two statements, especially if these are generic statements encapsulating dynamic control (Section 3.2), it is important to keep track of the accesses at hand. PPCG therefore computes so-called "tagged" dependence relations, where domain and range do not simply refer to a statement instance, but rather to a pair of statement instance and reference identifier. Consider, for example, the program fragment of Baghdadi et al. (2013a, Fig. 7), reproduced in Listing 11 with reference identifier annotations. The tagged dataflow relation of this fragment is

[ n ] -> { [ S1 [i , j ] -> R0 []] -> [ S2 [i , j ] -> R5 []] : 0 <= i , j < n ; [ S2 ][ i , j ] -> R3 []] -> [ S2 [i , j +1] -> R4 []] : 0 <= i < n and 0 <= j < n -1 } (19
) While constructing a schedule band, the scheduler needs to make sure that either the source and the sink of a live-range are assigned the same value by the band or that some validity constraints are respected that ensure that the live-range does not overlap with any other live-range on the same data element. The mechanism used by isl for communicating such information is the conditional validity schedule constraint. Such a constraint consists of two parts, a condition and a (conditional) validity constraint. Both of these are tagged relations, i.e., they relate two pairs of statement instance and reference identifier. The scheduler ensures that for each element of the condition relation, if domain and range are scheduled apart, that then all adjacent conditional validity constraints are satisfied. Given an element { (S(i) → U ()) → (T (j) → V ()) } of the condition relation, an adjacent conditional validity constraint is one where either the domain is the element { (T (j) → V ()) } or the range is the element

RT n°457 { (S(i) → U ()) }.
In order to exploit this mechanism, PPCG computes the following dependences. First recall that the dataflow analysis of isl takes as input a schedule, a set of "sink" accesses, a set of "must-source" accesses and a set of "maysource" accesses. The resulting dataflow relation matches each sink access to the last must-source access (according to the schedule) that accesses the same data element (if any) as well as to any intermediate may-source that accesses the same data element. If there is no matching must-source, then the sink is matched to all previous may-sources.

• the flow dependences are computed using the isl dataflow analysis with the may-reads as sinks, the may-writes as may may-sources and as mustsources the union of the must-writes and the kills (Section 3.8). The kills are subsequently removed from the domain of the dataflow relation to only retain dataflow between actual writes and reads. Removing these kills is important for the use of the dataflow relation for dead code elimination (Section 4.1), as it should not add kill statement instances to the set of statement instances that need to be executed.

• the order dependences are the ones that need to ensure that live-ranges (i.e., dataflow dependences) do not overlap in case they are not local to a schedule band. In particular, the order dependences maintain the initial order of the live-ranges and relate any read or unmatched write to any later write. An unmatched write is one that does not occur as the source in any live-range. In principle, such writes should have been removed by dead code elimination (Section 4.1), but the write may occur in a statement that performs a call or it may have been introduced due to the widening step.

It is important to include such unmatched writes in the order dependences to prevent these writes from getting scheduled in the middle of a later liverange. The writes that do appear in a live-range do not need to be included since they are necessarily executed before the corresponding read(s) and these reads are prevented from getting scheduled in the middle of a later live-range.

• the forced dependences are validity dependences that need to be enforced even when live-range reordering is used and consist of two parts, the external dependences and the order dependences between sources with a shared sink. The external dependences ensure that a live-in read remains live-in and that a live-out write remains live-out. They are computed by relating any live-in read to any later may-write to the same data element and similarly any may-write to any later live-out may-write to the same data element. The order dependences between sources with a shared sink ensure that such sources retain their relative execution order. If there are any may-writes, then there may be several flow dependences with the same sink. Only one of the corresponding sources will be the actual source at run-time, but some of the earlier potential sources may also be executed and they need to be prevented from getting executed after the actual source in the transformed program. Dependences are therefore introduced between any source of a live-range and any later source with a shared sink.

RT n°457

• the false dependences are kept for compatibility with the mode of PPCG where live-ranges are not reordered. They are computed using the isl dataflow analysis with the may-writes as sinks, the union of the may-writes and the may-reads as may-sources and the must-writes as must-sources.

Using these dependence relations, the schedule constraints for the isl scheduler are set as follows.

• the validity constraints are set to the union of the flow and the forced dependences.

• the proximity constraints are set to the union of the flow and the false dependences. Again, the reason for adding the false dependences is compatibility with the mode of PPCG where live-ranges are not reordered.

• the coincidence constraints, i.e., the constraints that determine which members of a band are considered "coincident", are set to the union of the flow dependences, the forced dependences and those order dependences that are derived from non-scalars. The reason for adding these order dependences is that PPCG currently only supports privatization of scalars. If there is an order dependence between two iterations of a schedule dimension then this means that the two iterations may access the same data elements. If this schedule dimension is considered to be coincident, then this data element needs to be privatized. Since PPCG can only handle such privatization for scalars, the presence of an order constraint between non-scalars should therefore prevent the schedule dimension from being considered coincident.

• the conditional validity constraints have the flow dependences as condition and (all) order dependences as (conditional) validity constraints.

The order dependences on the scalars are used to check which scalars should forcibly be mapped to private memory. Despite the fact that tiling was the main motivation for Baghdadi et al. (2013a), it should be noted that live-range reordering as implemented in PPCG is also useful for other transformations. Consider, for example, the example program of Mehta (2014, Figure 5.4 (a)), reproduced in Listing 12, with a couple of minor changes. First, the lower bound on the k1-loop has been changed from 0 to 1 to avoid an out-of-bounds access in statement S2. Second, an explicit kill has been added to tell PPCG that the contents of a0 and am1 do not need to be preserved. Otherwise, PPCG would have to make sure that the final writes to these scalars remain the last writes, implying that a loop containing such a final write could not be considered to be coincident. The approach of Mehta (2014, Section 5.4) appears to be more restrictive than live-range reordering. In particular, it only allows the outer two loops in Listing 12 to be fused, while live-range reordering allows all three loops to be fused.

Without live-range reordering, the dependences on a0 prevent essentially any reordering and the computed schedule tree is shown in Figure 1 (with the root domain node omitted). The generated C code would be essentially the same as the input code. With live-range reordering, the loops can be fused as long as no live-ranges overlap. In this case, the schedule tree in Figure 2 is computed (again with the root domain node omitted). Without the kills on a0 and am1, none of for ( int c0 = 0; c0 < Nx ; c0 += 1) for ( int c1 = 0; c1 < Ny ; c1 += 1) for ( int c2 = 1; c2 < Nz ; c2 += 1) { a0

-> [ { S5 [ i2 , j2 , k2 ] -> [( i2 )]; S1 [ i1 , j1 , k1 ] -> [( i1 )]; S2 [ i1 , j1 , k1 ] -> [( i1 )]; S3 [ i1 , j1 , k1 ] -> [( i1 )]; S4 [ i2 , j2 , k2 ] -> [( i2 )] } , { S5 [ i2 , j2 , k2 ] -> [( j2 )]; S1 [ i1 , j1 , k1 ] -> [( j1 )]; S2 [ i1 , j1 , k1 ] -> [( j1 )]; S3 [ i1 , j1 , k1 ] -> [( j1 )]; S4 [ i2 , j2 , k2 ] -> [( j2 )] } , { S5 [ i2 , j2 , k2 ] -> [(1 + k2 )]; S1 [ i1 ,
= a [ c0 ][ c1 ][ c2 ]; am1 = a [ c0 ][ c1 ][ c2 -1]; x [ c0 ][ c1 ][ c2 ] = ( a0 + am1 ); a0 = a [ c0 ][ c1 ][ c2 -1]; rho [ c0 ][ c1 ][ c2 -1] = ( a0 + ( x [ c0 ][ c1 ][ c2 ] -x [ c0 ][ c1 ][ c2 -1]
)); } Listing 13: Plain C code generated from the schedule tree in Figure 2 the schedule dimensions in the band would be considered coincident. Since PPCG insists that the outer schedule dimension is coincident (when generating CUDA or OpenCL code), this schedule would be rejected without the kills. Generating plain C code for the schedule in Figure 2 produces the code in Listing 13.

Absence of Dependences

The information collected by pet from a #pragma pencil independent as explained in Section 3.9 is currently only exploited by PPCG when live-range reordering (Section 4.4) is enabled (which is the default). The reason is that the live-range reordering support already takes care of ensuring that the live-ranges of variables local to the marked loops either do not overlap or that they are mapped to private variables on the device.

Currently, the information about which pairs of statement instances do not need to be protected from reordering or overlapping is only used after the dependences have been computed. In particular, the effect on the schedule constraints is as follows:

• the independent pairs of statement instances of which the first is not a must-write are removed from the flow dependences. Those pairs where the first is a must-write cannot simply be removed because this mustwrite may have killed other dependences that would have to be reinstated if the dependence were removed.

Note that it would arguably be better to adjust the dataflow analysis such that the unwanted dependences would not be added in the first place rather than trying to remove them afterwards. Such an adjusted dataflow analysis would be similar to the one proposed by [START_REF] Collard | Array Dataflow Analysis for Explicitly Parallel Programs[END_REF].

• the independent pairs of statement instances are removed from the external dependences. Note that the order dependences between sources with a shared sink are computed from the flow dependences after the independent pairs of statement instances have been removed from these flow dependences. These dependences therefore do not require any further adjustment.

RT n°457

• the independent pairs of statement instances are removed from the order dependences on arrays that are not local to the corresponding loop. Those that are local to the corresponding loop still need to be protected from overlapping, but for the non-local arrays, the user has guaranteed that no such protection is required.

• the independent pairs of statement instances are removed from the flow and forced dependences for the purpose of computing the coincidence constraints. Recall that the coincidence constraints are formed by taking the union of the flow dependence, the forced dependence and the order dependences on non-scalars. The independent pairs of statement instances have already been removed from the order dependences on non-local arrays. In the end, this means that the independent pairs of statement instances are removed from all contributors to the coincidence constraints except the order dependences on the local arrays.

Synchronization

There are two reasons why synchronization may need to be introduced.

• If any data is copied to/from shared memory, then the thread that copies the data may be different from the thread that reads/writes the data. Synchronization therefore needs to be added when any such copying is being performed. In particular, if there is any copying from global memory to shared memory (before the core computation), then synchronization needs to be added between this copying (which writes to shared memory) and the core computation (which reads from shared memory) and also between the core computation and the next iteration of the copying. Similarly, if there is any copying from shared memory back to global memory (after the core computation), then synchronization needs to be added between the core computation (which writes to shared memory) and the copying (which reads from shared memory) and also between this copying and the next iteration of the core computation.

• If any data is written to global memory, then it may be read again by the same block in a later iteration and this subsequent read may be performed by a different thread. Synchronization therefore needs to be added after any write to global memory that may be read again by the same block (and within the same iteration of all outer bands).

Note that PPCG currently does not treat these different types of synchronization differently in any way. In particular, they are not mapped to different kinds of barriers when targeting OpenCL.

Let us now consider in more detail how PPCG adds synchronization. PPCG first collects all the write references that require any synchronization on global memory. These are write references for which one of the corresponding reads or a conflicting write may occur in the same block but in a different thread. They are computed from the potential flow dependences and the pairs of potential writes accessing the same memory location that have the same schedule prefix for the kernel launch, but a different schedule prefix for the node that introduces the thread identifiers. If any of these references that require synchronization is mapped to global memory, then synchronization is added after RT n°457 the core computation unless the schedule depth of the node that is mapped to threads is the same as the schedule depth of the kernel. In the case where these depths are the same, there can be no loops within the kernel but outside of the synchronization, meaning that there is nothing to synchronize against.

Whenever PPCG adds extra statements for copying to/from shared or private memory, it also checks if any further synchronization is needed. Such statements and potential synchronization are added for each reference group separately. The reference groups are currently considered in a fairly arbitrary order. It may be better to first consider those groups that need to have copy statements introduced deeper in the schedule tree such that that the synchronization introduced for these groups could be taken into account when considering synchronization at higher levels.

If the reference group is mapped to private memory, then synchronization is only added when data is copied back from private memory to global memory. There is no need to add synchronization between the core computation and the copy because the copy is necessarily performed by the same thread. There is also no synchronization required before any copying from global memory to private memory because synchronization is introduced after every write to global memory that may have a corresponding read in the same kernel. Furthermore, synchronization is only added if the reference group contains any of the writes that require synchronization and if the schedule depth where the copy statement is introduced is strictly greater than the schedule depth of the kernel.

If the reference group is mapped to shared memory, then synchronization may be needed both after the copying and after the core computation. In case of a read from global memory to shared memory, this means that synchronization may be added before and after the core computation. Before adding this synchronization, PPCG first checks if such synchronization has already been introduced at the schedule depth where the copying is being performed or any deeper level. If so, then no additional synchronization needs to be introduced. Furthermore, the synchronization after the core computation is only added if the copy depth is strictly greater than the kernel depth. In case of a write from shared memory back to global memory, the synchronization is usually inserted before and after the copying. The synchronization before the copying comes after the core computation and can then also be avoided if such synchronization has already been introduced at the same or a deeper level. The synchronization after the copying is only introduced if the copy depth is strictly greater than the kernel depth. If the reference group contains any of the writes that require synchronization, then the synchronization is added after the copying. Otherwise, no synchronization is needed to protect the write to global memory, but the read from shared memory still needs to be protected from being overwritten by a later write. Note that the only later writes of concern are those performed by the core computation since any potential later read from global memory to shared memory is performed by the same thread that is performing the write from shared memory to global memory. It is therefore sufficient to add synchronization before the core computation at the schedule depth where the copying is introduced. As usual, this synchronization can be skipped if synchronization has already been introduced at the same or a deeper level.

RT n°457

  for ( int i = 0; i < 10; ++ i ) U : b [ i ] = a ; } Listing 5: Accesses to structures while the collected write access relation is

11 :

 11 One loop from the mvt kernel, with reference identifiers is only used after the schedule band has already been constructed. In the implementation on top of Pluto (Bondhugula et al. 2008) described by Baghdadi et al. (2013a), this is achieved by reinterpreting the bands computed by Pluto.

Figure 1 :

 1 Figure 1: Schedule tree for Listing 12 without live-range reordering

Figure 2 :

 2 Figure 2: Schedule tree for Listing 12 with live-range reordering

RT n°457RESEARCH CENTRE PARIS -Domainede Voluceau, -Rocquencourt B.P. 105 -78153 Le Chesnay Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-0803

Acknowledgements

This work was partly funded by the European FP7 project CARP id. 287767.

for ( int i1 = 0; i1 < Nx ; i1 ++) { for ( int j1 = 0; j1 < Ny ; j1 ++) { for ( int k1 = 1; k1 < Nz ; k1 ++) { S1 :

for ( int j2 = 0; j2 < Ny ; j2 ++) { for ( int k2 = 0; k2 < Nz -1; k2 ++) { S4 : RT n°457