Romain Postoyan
email: romain.postoyan@univ-lorraine.fr

Marcos César Bragagnolo
email: marcos-cesar.bragagnolo@univ-lorraine.fr

Ernest Galbrun
email: ernest.galbrun@univ-lorraine.fr

Jamal Daafouz
email: jamal.daafouz@univ-lorraine.fr

Dragan Nesic
email: dnesic@unimelb.edu.au

Eugênio B Castelan
email: eugenio.castelan@das.ufsc.br

D Nešić

Event-triggered tracking control of unicycle mobile robots

Keywords: Event-triggered control, hybrid systems, nonholonomic systems, nonlinear systems

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The objective of the paper is to guarantee the tracking of given reference trajectories by unicycle mobile robots using a remote controller, while reducing the usage of the communication channel. In particular, we want to limit the amount of control input updates to save communication resources and to reduce the risk of packet losses and long transmission delays; noting that that also allows to curb the actuators wear and to reduce the energy consumption of the actuators. An event-triggered feedback law is designed for that purpose. The idea of event-triggered control is to transmit data between the controller and the plant whenever a state-dependent criterion is satisfied and not periodically as in traditional setups. In that way, the transmissions are adapted to the state of the system and these only occur when it is needed. Various techniques have been developed, see [START_REF] Arzén | A simple event-based PID controller[END_REF][START_REF] Aström | Comparison of Riemann and Lebesgue sampling for first order stochastic systems[END_REF][START_REF] Heemels | Analysis of event-driven controllers for linear systems[END_REF][START_REF] Postoyan | A framework for the event-triggered stabilization of nonlinear systems[END_REF][START_REF] Tabuada | Event-triggered real-time scheduling of stabilizing control tasks[END_REF] to mention a few. Most of them address the stabilization of equilibrium points, while very few controllers have been synthesized to stabilize time-varying trajectories, see e.g., [START_REF] Tallapragada | On event triggered tracking for nonlinear systems[END_REF]. It appears that tracking control induces additional difficulties as only approximate tracking can usually be ensured under communication constraints because of the time-varying component of the control law, see for more detail [START_REF] Postoyan | Tracking control for nonlinear networked control systems[END_REF]. As a consequence, available results on the stabilization of equilibrium points are not directly applicable in this context.

We follow an emulation-like approach (see e.g., [START_REF] Tabuada | Event-triggered real-time scheduling of stabilizing control tasks[END_REF][START_REF] Postoyan | A framework for the event-triggered stabilization of nonlinear systems[END_REF]). Thus, we first design the controller while ignoring the communication constraints and we derive the triggering strategy afterwards. We have selected the state-feedback controller of [START_REF] Jiang | Tracking control of mobile robots: a case study in backstepping[END_REF] among others, because the law is continuous and an explicit Lyapunov function is provided. The continuity of the feedback law is useful to guarantee the existence of a minimum amount of time between two transmissions, which is essential in practice as the hard-ware cannot transmit infinitely fast. The Lyapunov function is used to design the event-triggering condition. This Lyapunov function is weak in the sense that it does not strictly decrease along the solutions of the closed-loop system in the absence of communication constraints, see [START_REF] Jiang | Tracking control of mobile robots: a case study in backstepping[END_REF]. To overcome this issue, we impose conditions on the reference trajectories under which the Lyapunov function strictly decreases outside a given neighborhood of the origin, which can be made as small as desired by appropriately tuning the control law parameters. We then take into account the effects of the network and we model the problem as a hybrid system using the formalism of [START_REF] Goebel | Hybrid dynamical systems[END_REF] (as in e.g., [START_REF] Donkers | Outputbased event-triggered control with guaranteed L∞-gain and improved and decentralized event-triggering[END_REF][START_REF] Postoyan | A framework for the event-triggered stabilization of nonlinear systems[END_REF][START_REF] Seuret | Stability of nonlinear systems by means of event-triggered sampling algorithms[END_REF]). The event-triggering condition we construct is inspired by the technique in [START_REF] Tabuada | Event-triggered real-time scheduling of stabilizing control tasks[END_REF][START_REF] Forni | Eventtriggered transmission for linear control over communication channels[END_REF] where the Lyapunov function is forced to decrease at a certain rate, except in the aforementioned neighborhood of the origin, which we call a dead-zone. We prove that the solutions of the robot model practically converge towards the reference trajectories. The ultimate bound depends on a tuneable parameter which corresponds to the 'size' of the dead-zone of the triggering condition. It can be used to adjust the accuracy of the tracking error at the price of more transmissions. The approach we follow is similar to the one in [START_REF] Tallapragada | On event triggered tracking for nonlinear systems[END_REF], however we do not impose the same conditions on the reference trajectories and we propose a different triggering law. Finally, we have implemented the event-triggered controller on a benchmark. The controller sends its data to the robot over a wireless network according to the triggering condition and the measurements are periodically collected by cameras which are directly connected to the remote controller via a dedicated high bandwidth wired channel. The experimental results show that it is possible to significantly reduce the usage of the wireless network compared to periodic sampling, if we agree to slightly give up on the tracking accuracy.

A similar problem is investigated in [START_REF] Santos | Adaptive selftriggered control of a remotely operated robot[END_REF] where selftriggered controllers are developed. Beside the fact that we consider a nonlinear model of the robot and not a linear model as in [START_REF] Santos | Adaptive selftriggered control of a remotely operated robot[END_REF], we also envision a different setup. In our case, no local controller is implemented on the robot which helps saving batteries and the measurements are given by cameras which are connected via a wired channel to the controller unit (as opposed to odometry in [START_REF] Santos | Adaptive selftriggered control of a remotely operated robot[END_REF]), which justifies the choice of an eventtriggered implementation. Compared to our preliminary work in [START_REF] Postoyan | Nonlinear event-triggered tracking control of a mobile robot: design, analysis and experimental results[END_REF], the results rely on a simplified stability analysis and on different assumptions on the reference trajectories which cover a class of well-motivated cases.

In particular, we are able to guarantee the existence of a minimum amount of time between two transmissions for both of the reference trajectories we consider in the experimental part, which is not the case of [START_REF] Postoyan | Nonlinear event-triggered tracking control of a mobile robot: design, analysis and experimental results[END_REF]. Furthermore, we present new experimental results as we use cameras to measure the robot position as explained above and not odometry as in [START_REF] Postoyan | Nonlinear event-triggered tracking control of a mobile robot: design, analysis and experimental results[END_REF]; noting that the latter has the drawback to use the wireless network periodically to communicate with the controller.

The paper is organized as follows. Preliminaries are given in Section 2. In Section 3, the controller of [START_REF] Jiang | Tracking control of mobile robots: a case study in backstepping[END_REF] is recalled and the hybrid model is introduced. The event-triggering condition and the analytical results are presented in Section 4. Section 5 deals with the experimental results. The proofs are provided in the Appendix.

2 Preliminaries Let R = (-∞, ∞), R ≥0 = [0, ∞), R >0 = (0, ∞), Z ≥0 = {0, 1,
2, . . .}, and Z >0 = {1, 2, . . .}. For (x, y) ∈ R n+m , the notation (x, y) stands for [x T , y T ] T . A function γ : R ≥0 → R ≥0 is of class K if it is continuous, zero at zero and strictly increasing, and it is of class

K ∞ if in addition it is unbounded. A continuous function γ : R 2 ≥0 -→ R ≥0 is of class KL if for each t ∈ R ≥0 , γ(•, t) is of class K,
and, for each s ∈ R >0 , γ(s, •) is decreasing to zero. For a right-continuous function f : R → R n , f (t + ) stands for lim s→t, s>t f (s) for t ∈ R. We recall that the function sinc is defined from R to R by sinc(x) = sin x x when x = 0, and sinc(0) = 1, and it is twice continuously differentiable with sinc ′ (0) = 0 and sinc ′′ (0) = - 1 3 .

We will write the event-triggered controlled system as a hybrid system using the formalism of [START_REF] Goebel | Hybrid dynamical systems[END_REF]. In particular, we will consider a system of the form

ẋ = f (x) x ∈ C, x + = g(x) x ∈ D, (1)
where x ∈ R n is the state and C, D ⊂ R n are respectively the flow and the jump sets. We assume that f and g are continuous and that C and D are closed sets (this will be the case in the paper). We recall some basic definitions, see [START_REF] Goebel | Hybrid dynamical systems[END_REF].

A set E ⊂ R ≥0 × Z ≥0 is a compact hybrid time domain if E = ∪ j∈{0,...,J-1} ([t j , t j+1 ], j) for some finite sequence of times 0 = t 0 ≤ t 1 ≤ . . . ≤ t J . The set E is a hybrid time domain if for all (T, J) ∈ E, E ∩ ([0, T ] × {0, 1, . . . , J}) is a compact hybrid time domain.
A hybrid arc is a function φ defined on a hybrid time domain dom φ and such that, for each j ∈ Z ≥0 , t → φ(t, j) is locally absolutely continuous on 1) is said to be: nontrivial if dom φ contains at least two points; maximal if it cannot be extended; complete if dom φ is unbounded.

I j := {t : (t, j) ∈ dom φ}. A hybrid arc φ : dom φ → R n is a solution to (1) if: (i) φ(0, 0) ∈ C ∪ D; (ii) for all j ∈ Z ≥0 and almost all t ∈ I j , φ(t, j) ∈ C and φ(t, j) = f (φ(t, j)); (iii) for (t, j) ∈ dom φ such that (t, j + 1) ∈ dom φ, φ(t, j) ∈ D and φ(t, j + 1) = g(φ(t, j)). A solution to (
The definition below characterizes hybrid systems that generate solutions for which two successive jumps are spaced by (at least) a strictly positive amount of time which is uniform over the ball of initial conditions (see [START_REF] Postoyan | A framework for the event-triggered stabilization of nonlinear systems[END_REF]). Definition 1 The solutions to (1) have a uniform semiglobal dwell-time if for any ∆ ≥ 0, there exists η(∆) > 0 such that for any solution φ to (1) with

|φ(0, 0)| ≤ ∆, sup I j -inf I j ≥ η(∆) ∀j ∈ Z >0 with I j = ∅. (2) 2
We recall the following invariance definition (see [START_REF] Goebel | Hybrid dynamical systems[END_REF]).

Definition 2 Consider system (1), the set A ⊂ R n is strongly pre-forward invariant if for any ξ ∈ A each solution φ with initial condition ξ satisfies φ(t, j) ∈ A for all (t, j) ∈ dom φ. 2

3 System model

Jiang and Nijmeijer's controller

We consider a mobile robot for which the dynamics are defined by the following nonholonomic system

ẋ = v cos(θ), ẏ = v sin(θ), θ = w, (3)
where (x, y) are the Cartesian coordinates and θ is the angle between the heading direction and the x-axis, (v, w) denotes the control input vector, see Figure 1.

The objective is to make the solution to system (3) track a given trajectory t → (x r (t), y r (t), θ r (t)). We restrict our attention to the case where the reference trajectory satisfies the following dynamical equations ẋr = v r cos(θ r ), ẏr = v r sin(θ r ), θr = w r ,

where v r and w r are the feedforward inputs to system (4). In practice, the term w r is obtained by differentiating θ r . The term v r is deduced by using that v 2 r = ( ẋr ) 2 +( ẏr ) 2 and its sign is obtained by investigating the first equations of (4) for the given trajectory (x r , y r , θ r ).

We refer to (4) as the reference system. We use the following error coordinates (see [START_REF] Kanayama | A stable tracking control method for an autonomous mobile robot[END_REF])

    x e y e θ e     =     cos(θ) sin(θ) 0 -sin(θ) cos(θ) 0 0 0 1         x r -x y r -y θ r -θ     (5)
which leads to the error system ẋe = wy ev + v r cos(θ e ), ẏe = -wx e + v r sin(θ e ) θe = w rw.

(6) The change of coordinates in ( 5) is invertible and the transformation matrix and its inverse are bounded. As a consequence, if we guarantee that x e , y e , θ e converge to a neighborhood of the origin, then so do x r -x, y r -y, θ rθ. In [START_REF] Jiang | Tracking control of mobile robots: a case study in backstepping[END_REF], the following controller is shown to ensure the global convergence of x e , y e , θ e towards the origin (under some conditions on the feedforward inputs v r and w r )

v = v 1 + c 1 xe , w = w r + γy e v r
with c 3 > 0. We do not necessarily require to differentiate online ẇ to compute the term v 1 in (8) as we can replace ẇ by its algebraic expression by using (3), ( 4) and [START_REF] Heemels | Analysis of event-driven controllers for linear systems[END_REF].

Assumptions on the feedforward inputs

We make the following assumption on the feedforward inputs v r and w r .

Assumption 3 There exists M ≥ 0 such that the absolute values of v r (t), w r (t), vr (t) and ẇr (t) are bounded by M ≥ 0 for any t ∈ [0, ∞). Furthermore there exists

c ∈ R >0 such that v r (t) 2 ≥ c 2 for any t ≥ 0. 2
The first part of Assumption 3 corresponds to the first condition of Proposition 2 in [START_REF] Jiang | Tracking control of mobile robots: a case study in backstepping[END_REF]. We also assume that the sign of v r is always the same and that v r does not become arbitrarily close to the origin. The latter differs from [START_REF] Jiang | Tracking control of mobile robots: a case study in backstepping[END_REF][START_REF] Postoyan | Nonlinear event-triggered tracking control of a mobile robot: design, analysis and experimental results[END_REF] and is satisfied by various reference trajectories, as those studied in Section 5 for example.

Network implementation & hybrid model

We focus on the scenario where the controller (v, w) in [START_REF] Heemels | Analysis of event-driven controllers for linear systems[END_REF] communicates with the robot (3) via a network, see Figure 2. We assume that the controller has access to the reference trajectory and the feedforward inputs v r and w r at any time.

The sequence of transmission instants is t i , i ∈ Z ≥0 , and it is defined according to a state-dependent criterion which will be specified in Section 4.1. Thus, we have, for all t ∈ [t i , t i+1 ), ẋe = ŵy ev + v r cos(θ e ), ẏe =ŵx e + v r sin(θ e ) θe = w rŵ

where v and ŵ are respectively the networked versions of v and w. The controller is implemented using zeroorder-hold devices, i.e. v(t) = v(t i ) and ŵ(t) = w(t i ) for t ∈ [t i , t i+1 ). Hence the controller is, for i ∈ Z ≥0 ,

ẇ = 0 v = 0 t ∈ [t i , t i+1 ) ŵ(t + i ) = w(t i ) v(t + i ) = v(t i ). ( 10
)
We write the model using the formalism of [START_REF] Goebel | Hybrid dynamical systems[END_REF], for which a jump represents a transmission q = f (q) q ∈ C, q + = g(q) q ∈ D,

with q := (x e , y e , θ e , ŵ, v, τ ),

f (q) :=             ŵy e -v + v r cos(θ e ) -ŵx e + v r sin(θ e ) w r - ŵ 0 0 1             , g(q) :=             x e y e θ e v w τ             , ( 12
)
where τ is the time-variable, which is introduced to render the model autonomous by considering the feedforward inputs as time-dependent terms, and v, w are de-fined in [START_REF] Heemels | Analysis of event-driven controllers for linear systems[END_REF]. The sets C and D are defined by the triggering condition. Indeed, the set C corresponds to the region of the state-space where the triggering condition is not satisfied and D is defined as the region of the statespace where the triggering condition is verified (which does not mean that C ∩ D = ∅ as we will see). In that way, defining the triggering condition is equivalent to defining the sets C and D (see e.g., [START_REF] Postoyan | A framework for the event-triggered stabilization of nonlinear systems[END_REF]). The model ( 11) is a simplification of the experimental setup considered in Section 5 as we assume that the triggering condition is continuously evaluated in [START_REF] Postoyan | Nonlinear event-triggered tracking control of a mobile robot: design, analysis and experimental results[END_REF]. Nevertheless, system [START_REF] Postoyan | Nonlinear event-triggered tracking control of a mobile robot: design, analysis and experimental results[END_REF] is useful to design the event-triggered controller, and we will see that the stability property we ensure for system [START_REF] Postoyan | Nonlinear event-triggered tracking control of a mobile robot: design, analysis and experimental results[END_REF] can be maintained when the sampling of the triggering mechanism is sufficiently small. Our objective is to design a triggering condition, i.e. the sets C and D in [START_REF] Postoyan | Nonlinear event-triggered tracking control of a mobile robot: design, analysis and experimental results[END_REF], such that the errors x e , y e , θ e practically converge to the origin. The tuneable parameter will be given by the triggering condition.

Design & analytical results

Event-triggering condition

The Lyapunov-like candidate function below is considered in [START_REF] Jiang | Tracking control of mobile robots: a case study in backstepping[END_REF] V

(q) = 1 2 x2 e + 1 2 y 2 e + 1 2γ θ 2 e , (13)
where xe is defined in (8) and γ comes from [START_REF] Heemels | Analysis of event-driven controllers for linear systems[END_REF]. We derive, by following similar lines as in [START_REF] Jiang | Tracking control of mobile robots: a case study in backstepping[END_REF], that for any q ∈ R 6 , ∇V (q), f (q) = -Σ(q) + Λ(q) [START_REF] Santos | Adaptive selftriggered control of a remotely operated robot[END_REF] where Σ(q) := c 1 x2 e + c 2 θ 2 e + c 3 w 2 y 2 e and Λ(q) := xe (e w y ee vc 3 y e e w ρ(x e , y e , θ e , v r ) + c 3 we w x e )y e e w x e -1 γ θ e e w , with ρ(x e , y e , θ e , v r ) := -γx e v r sincθ eγy e v r sinc ′ (θ e )+ c 2 γ, e w := ŵw and e v : vv. We note that we recover (41) in [START_REF] Jiang | Tracking control of mobile robots: a case study in backstepping[END_REF] when there is no network, i.e. when e v = e w = 0 (as Λ(q) = 0 in this case).

We define the flow and the jump sets as follows

C = q ∈ Q : ∇V (q), f (q) ≤ -σΣ(q) or V (q) ≤ ε 2 D = q ∈ Q : ∇V (q), f (q) ≥ -σΣ(q) and V (q) ≥ ε 2
(15) where Q := R 5 ×R ≥0 , σ ∈ (0, 1) and ε > 0 are design parameters. The idea is the following. In [START_REF] Jiang | Tracking control of mobile robots: a case study in backstepping[END_REF], the function V is ensured to decrease with the rate Σ. Here, we degrade this decreasing property by the factor σ so that V still decreases but with rate σΣ with σ ∈ (0, 1), like in [START_REF] Tabuada | Event-triggered real-time scheduling of stabilizing control tasks[END_REF]. In addition, we impose that condition only when V (q) ≥ ε2 as otherwise we can neither guarantee the convergence of x e , y e , θ e to the origin, nor ensure the existence of a minimum amount of time between two transmission instants under Assumption 3. That is in agreement with [START_REF] Postoyan | Tracking control for nonlinear networked control systems[END_REF] where it is shown that the time-triggered emulation of a nonlinear tracking controller leads to approximate tracking when the feedforward inputs are subject to communication constraints, which is the case here, even with very fast sampling. Similar ideas are followed in [START_REF] Forni | Eventtriggered transmission for linear control over communication channels[END_REF][START_REF] Postoyan | A framework for the event-triggered stabilization of nonlinear systems[END_REF][START_REF] Seuret | Stability of nonlinear systems by means of event-triggered sampling algorithms[END_REF] to design the triggering condition. Nevertheless, none of the strategies proposed in these papers directly applies to our problem as the triggering condition is based on a weak Lyapunov function V (as opposed to strong Lyapunov functions) and these do not address tracking control.

Stability analysis

We first state the following lemma.

Lemma 4 Suppose Assumption 3 holds, then for any c 1 , c 2 , c 3 , ε > 0 and γ > 0 sufficiently large, there exists ϑ > 0 such that for any q ∈ Q

V (q) ≥ ε 2 ⇒ Σ(q) ≥ ϑ. ( 16
)
2

When q ∈ C and V (q) > ε 2 , we know from (15) that ∇V (q), f (q) ≤ -σΣ(q). Lemma 4 implies that in this case ∇V (q), f (q) ≤ -σϑ < 0, i.e. that V strictly decreases on this part of the set C. That is very useful for the stability analysis as it allows us avoiding invoking hybrid invariance principle as in [START_REF] Postoyan | Nonlinear event-triggered tracking control of a mobile robot: design, analysis and experimental results[END_REF]. Moreover, Lemma 4 also plays a key role to prove the existence of a minimum amount of time between two jumps (see the proof of Theorem 5). We are ready to state the main result of this section.

Theorem 5 Consider system [START_REF] Postoyan | Nonlinear event-triggered tracking control of a mobile robot: design, analysis and experimental results[END_REF] and suppose Assumption 3 holds. For any c 1 , c 2 , c 3 , ε > 0 and γ > 0 sufficiently large, the following holds.

(i) There exists β ∈ KL such that for any solution φ, V (φ(t, j)) ≤ β(V (φ(0, 0)), t+j)+ε 2 for any (t, j) ∈ dom φ. (ii) The solutions have a uniform semiglobal dwell-time. (iii) The maximal solutions are complete.

2

Items (i) and (iii) of Theorem 5 together with the boundedness of v r and w r stated in Assumption 3 ensure that the states of the robot converge to the states of the reference system up to an error of the order of ε. Item (ii) of Theorem 5 satisfies the requirement that any two successive transmissions cannot occur arbitrarily close-intime, which is important in practice.

Digital implementation

As already mentioned, the triggering condition is not continuously evaluated in practice but only periodically1 every T units of time, see Section 5. The sampling effect can be modeled as a time-varying delay δ ∈ [0, T ] which affects the control inputs. Indeed, the control input applied to the plant is sent at least after T units of time have elapsed since the triggering condition [START_REF] Seuret | Stability of nonlinear systems by means of event-triggered sampling algorithms[END_REF] has been violated. The proof of Theorem 5 indicates that in this case, for T sufficiently small, there exists σ ′ ∈ (0, σ) such that the triggering condition defined by ( 15) is satisfied with σ ′ instead of σ as explained in [START_REF] Tabuada | Event-triggered real-time scheduling of stabilizing control tasks[END_REF]. Note that, while it is possible to use the proof of Theorem 5 to compute an upper-bound of T , we do not provide it because the developments carried out in the appendix are subject to some conservatism which would typically lead to a small constant.

Experimental validation

We have implemented the proposed strategy on a benchmark. The event-triggered controller is located on a remote computer which communicates with the robot via a IEEE 802.11g wireless network. The state of the robot is measured by cameras which communicate with the control unit via a dedicated wired network. Because of the time-interval between each position update and the computation time, the triggering condition is periodically verified every T eval = 50 ms, and not continuously as in Section 3-4. Experimental results show that the period T eval is sufficiently small to maintain the properties ensured by Theorem 5, which is consistent with Section 4.3.

The reduction of the usage of the communication channel only concerns the transmissions from the controller to the robot, which corresponds to the dashed channel in Figure 2 (and not from the robot to the controller where periodic sampling is applied). The controller parameters in [START_REF] Heemels | Analysis of event-driven controllers for linear systems[END_REF] are c 1 = 30, c 2 = 0.8, c 3 = 0.9, γ = 1.5.

We focus on two types of reference trajectories whose paths respectively describe an ellipse and a lemniscate and which ensure the satisfaction of Assumption 3. The triggering condition [START_REF] Seuret | Stability of nonlinear systems by means of event-triggered sampling algorithms[END_REF] has been implemented with σ = 0.5 and we have studied the influence of ε. Simulations have been run over 2 T simulation = 60s. The usage of the communication network is quantified using the percentage (# stands for 'number' below) The performance of the controller is evaluated by the constant e that is the maximum value of |(xx r , yy r , θ -θ r )|, which is (heuristically) normalised according to the dimensions of the system, over the second half of the experiment, i.e. In that way, e quantifies the tracking accuracy of the robot after the transient period. The scaling coefficient 0.1 in (18) approximatively corresponds to the robot diameter, i.e. 100 mm. Table 1 shows that the controller only uses the communication channel 23.6% of the time, in the worst case, and that this quantity can be significantly decreased by varying ε. Thus, we see that to increase ε reduces the usage of the communication channel at the price of degraded tracking performances (in general).

r := 100 × # control updates × T eval T simulation . (17)
We then compared the event-triggered controller with a corresponding time-triggered one. The same controller [START_REF] Heemels | Analysis of event-driven controllers for linear systems[END_REF] is considered but the control input is transmitted to the robot every T sampling seconds. We are not aware of any analytical result which allows to derive an upper bound on T sampling . Hence, we have selected T sampling to be the average inter-transmission interval of the corresponding event-triggered controller. In that way, both controllers do generate the same amount of transmissions over an experiment. Figure 3 shows the average values of e for different ε obtained over 4 tests. We see that both controllers generate a similar tracking performance for small values of ε and that the event-triggered controller outperforms the time-triggered one for large values of ε, i.e. when the wireless network is used less often. These results suggest that it is possible to reduce the usage of the wireless network using the event-triggered controller provided that we agree to give up a bit on the tracking accuracy.

Conclusion

An event-triggered tracking controller has been designed and implemented for a standard model of mobile robots.

A global practical stability property together with the existence of a strictly positive amount of time between two transmissions have been proved under appropriate conditions on the reference trajectories to be tracked. Finally, the experimental results illustrate the potential of event-triggering techniques over periodic sampling and support the on-going research on nonlinear eventtriggered control in general, as these demonstrate the benefits of the paradigm on a real system.

A Appendix

Proof of Lemma 4. Let q ∈ Q such that V (q) ≥ ε 2 . In view of the definition of V in (13), the following three cases can occur: (i) 16) holds with ϑ = 2 3 ε 2 min{c 1 , c 2 γ} for cases (i) and (ii). Case (iii) requires more attention. Suppose that θ 2 e ≤ ε2 with ε > 0 sufficiently small such that 2 e ≥ ε2 and ( 16) holds with ϑ = c 2 ε2 ). We want to ensure w 2 ≥ 1 by appropriately selecting γ. Since w = w r + γ(y e v r sincθ e + c 2 θ e ), w 2 is a polynomial in γ. The idea is to first prove that the leading term (y e v r sincθ e + c 2 θ e ) 2 is bigger than a constant λ > 0. We can then select γ sufficiently large to obtain that w 2 ≥ 1. It holds that, in view of Assumption (ii) of Theorem 5 follows by induction by repeating the same arguments for t j+1t j with j ∈ Z ≥0 .

1 2 x2 e ≥ ε 2 3 ; (ii) 1 2γ θ 2 e ≥ ε 2 3 ; (iii) 1 2 y 2 e ≥ ε 2 3 . Since Σ(q) = c 1 x2 e + c 2 θ 2 e + c 3 w 2 y 2 e , (
Item (i) of Theorem 5. Let r > 0 and ε ′ > ε and consider a solution φ to [START_REF] Postoyan | Nonlinear event-triggered tracking control of a mobile robot: design, analysis and experimental results[END_REF] such that V (φ(0, 0)) ≤ r. Let (t, j) ∈ dom φ with j ≥ 1 and t 1 = inf I 1 . Since V (φ(0, 0)) ≤ r, there exists ∆ ≥ 0 such that |(φ xe (t 1 , 1), φ ye (t 1 , 1), φ θe (t 1 , 1), φ ŵ (t 1 , 1), φ v (t 1 , 1))| ≤ ∆ for (t 1 , 1) ∈ dom φ in view of [START_REF] Postoyan | Tracking control for nonlinear networked control systems[END_REF], Assumption 3 and (12). According to item (ii) of Theorem 5, we derive that tt 1 ≥ (j -1)η(r) where we write η(r) (instead of η(∆) with some abuse of notation). We notice that V (φ(t, j)) ≤ max{V (φ(0, 0))σϑt, ε 2 } ≤ max{rσϑt, ε 2 } in view of (A.3) and (A.4). Consequently, t + j ≥ T ′ with T ′ = (1 + η(r) -1 ) -1 ( r-ε ′2 σϑ -1), implies V (φ(t, j)) ≤ ε ′2 . Note that when t 1 = ∞, we have that dom φ = R ≥0 × {0} and we obtain t + j = t ≥ r-ε ′2 σϑ implies V (φ(t, j)) ≤ ε ′2 for (t, j) ∈ dom φ. We have thus proved that for any r > 0, ε ′ > ε, there exists T > 0 such that for any solution φ to [START_REF] Postoyan | Nonlinear event-triggered tracking control of a mobile robot: design, analysis and experimental results[END_REF], t + j ≥ T with (t, j) ∈ dom φ implies V (φ(t, j)) ≤ ε ′2 . We now use this uniform global attractivity property together with the first inequality in (A.7) to prove item (i) of Theorem 5. We follow similar arguments as in the proof of Theorem 3.40 in [START_REF] Goebel | Hybrid dynamical systems[END_REF] for that purpose. We consider β 0 (r, s) = max{sup{V (φ(t, j)) : φ ∈ S(ξ), V (ξ) ≤ r, t + j ≥ s}ε 2 , 0} + r exp(-s) for r, s ≥ 0 (where S(ξ) is the set of the solutions to [START_REF] Postoyan | Nonlinear event-triggered tracking control of a mobile robot: design, analysis and experimental results[END_REF] starting at ξ) which ensures V (φ(t, j))ε 2 ≤ β 0 (V (φ(0, 0)), t + j) for any solution φ to [START_REF] Postoyan | Nonlinear event-triggered tracking control of a mobile robot: design, analysis and experimental results[END_REF]. The β 0 is strictly increasing in r and decreasing in s and β 0 (r, s) ≤ max{rε 2 , 0} + r in view of (A.7). Hence lim r→0 + β 0 (r, s) = 0 for each s ≥ 0. Furthermore, lim s→∞ β(r, s) = 0 in view of the uniform global attractivity shown above. Since β 0 is not necessarily continuous, we use Lemma 3.41 in [START_REF] Goebel | Hybrid dynamical systems[END_REF] to upper-bound it by β ∈ KL. As a consequence, item (i) of Theorem 5 holds.

Completeness of the maximal solutions Let φ be a maximal solution to [START_REF] Postoyan | Nonlinear event-triggered tracking control of a mobile robot: design, analysis and experimental results[END_REF]. In view of [START_REF] Postoyan | A framework for the event-triggered stabilization of nonlinear systems[END_REF], φ is nontrivial as f (q) always belong to the tangent cone to C at q (see Definition 5.12 in [START_REF] Goebel | Hybrid dynamical systems[END_REF]) for any q ∈ C\D, see Proposition 6.10 in [START_REF] Goebel | Hybrid dynamical systems[END_REF]. Noting that C ∪ D = Q, the only potential obstacle to the completeness of φ is if it explodes in finite-time according to Proposition 6.10 in [START_REF] Goebel | Hybrid dynamical systems[END_REF], but that is not possible in view of (A.8) and the dynamics of the τ -variable. Hence φ is complete.

2

Fig. 1 .

 1 Fig. 1. Mobile robot.

  sincθ e + c 2 γθ e , (7) where c 1 , c 2 , γ > 0 are design parameters and xe := x ec 3 wy e v 1 := v r cos θ ec 3 ẇy e + c 3 w(wx ev r sin θ e )

Fig. 3 .

 3 Fig. 3. Experimental results when the reference describes an ellipse (top) and a lemniscate (bottom).

1 2

 1 εc (otherwise θ 2

We actually have a periodic event-triggered controller[START_REF] Heemels | Periodic event-triggered control for linear systems[END_REF][START_REF] Postoyan | Periodic event-triggered control for nonlinear systems[END_REF].

Tests have shown that the results do not significantly vary for bigger values of T simulation .

Acknowledgements

The authors are grateful to Lucas Bchini for his help for the experimental results reported in Section 5, and to the anonymous reviewers for their valuable comments which lead to significant improvements.

This work was supported by the ANR under the grant COM-PACS (ANR-13-BS03-0004-02), the European 7th Framework Network of Excellence Highly-complex and networked control systems (HYCON2 No. 257462) and the Australian Research Council under the Discovery Projects and Future Fellowship schemes. Corresponding author R. Postoyan.

(A.1) In view of the conditions assumed on θ e , (y e v r sincθ e + c 2 θ e ) 2 ≥ λ for some λ > 0. Consequently, for w r such that |w r | ≤ M (which is the case in our study according to Assumption 3),

2) It suffices to select γ sufficiently large for the right handside of the last inequality above to be bigger than 1 (using the fact that |y e v r sincθ e + c 2 θ e | ≥ √ λ). We have proved that [START_REF] Tabuada | Event-triggered real-time scheduling of stabilizing control tasks[END_REF]

Proof of Theorem 5. Let c 1 , c 2 , c 3 , ε > 0 and γ > 0 be such that (16) holds. Let q ∈ C. According to the definition of the set C, V (q) > ε 2 implies ∇V (q), f (q) ≤ -σΣ(q), which leads to, in view of Lemma 4,

On the other hand, it holds that, for any q ∈ D, V (g(q)) ≤ V (q). (A.4)

Boundedness of the solutions. We now derive a boundedness property for the solutions to [START_REF] Postoyan | Nonlinear event-triggered tracking control of a mobile robot: design, analysis and experimental results[END_REF] based on (A.3) and (A.4), which will be used to prove that any maximal solution is complete and that the solutions have a uniform semiglobal dwell-time. Let ∆ > 0 and φ = (φ xe , φ ye , φ θe , φ ŵ , φ v , φ τ ) be a solution to [START_REF] Postoyan | Nonlinear event-triggered tracking control of a mobile robot: design, analysis and experimental results[END_REF] Note that γ|(x e , θ e , y e )| 2 ≤ V (q) ≤ γ|(x e , θ e , y e )| 2 with γ = 1 2 min{1, 1 γ } and γ = 1 2 max{1, 1 γ } for any q ∈ Q. In that way, we derive from (A.3), (A.4) and (A.6) that, for (t, j) ∈ dom φ,

) from which we obtain |(φ xe (t, j), φ ye (t, j), φ θe (t, j))| ≤ ∆ V for some ∆ V > 0 which only depends on the controller parameters, ∆, M and ε. As a consequence, using Assumption 3, we deduce that w(t, j) and v(t, j) are bounded by a constant ∆ v,w > 0 which depends on the same constants as ∆ V , for any (t, j) ∈ dom φ in view of [START_REF] Heemels | Analysis of event-driven controllers for linear systems[END_REF]. In view of ( 12), |φ ŵ (t, j)| ≤ ∆ w,v for (t, j) ∈ dom φ with j ≥ 1 and |φ ŵ (0, 0)| ≤ ∆ (since |φ(0, 0)| ≤ ∆), we derive that φ ŵ is bounded by ∆ v,w + ∆. We similarly deduce that |φ v (t, j)| ≤ ∆ v,w + ∆ for any (t, j) ∈ dom φ. We have shown that for all (t, j) ∈ dom φ |(φ xe (t, j), φ ye (t, j), φ θe (t, j), φ ŵ (t, j), φ v (t, j))| ≤ ∆ (A.8) for some ∆ > 0 which only depends on the controller parameters, ∆, M and ε.

Uniform semiglobal dwell-times. We now show that the solutions to [START_REF] Postoyan | Nonlinear event-triggered tracking control of a mobile robot: design, analysis and experimental results[END_REF] have a uniform semiglobal dwell-time. Let ∆ ≥ 0, φ be a solution to [START_REF] Postoyan | Nonlinear event-triggered tracking control of a mobile robot: design, analysis and experimental results[END_REF] with |φ(0, 0)| ≤ ∆, (t 1 , 1) ∈ dom φ and (t 2 , 1) where t 1 = inf I 1 and t 2 = sup I 1 (note that if I 1 = ∅ or I 1 is unbounded, φ jumps no more than once and the desired result holds). Necessarily V (φ(t 1 , 1)) ≥ ε 2 as φ(t 1 , 0) ∈ D and so V (φ(t 1 , 0)) = V (φ(t 1 , 1)) ≥ ε 2 . If there exists t ∈ I 1 such that V (φ(t, 1)) ≤ ε 2 then t 2 = ∞ as the set {q : V (q) ≤ ε 2 } is strongly pre-forward invariant for system [START_REF] Postoyan | Nonlinear event-triggered tracking control of a mobile robot: design, analysis and experimental results[END_REF] in view of (A.3) and (A.4), in which case (2) holds for any η > 0. When for any t ∈ I 1 , V (φ(t, 1)) ≥ ε 2 , we deduce that t 2t 1 is bigger than the time T it takes for Λ(φ) to grow from 0 (its value at (t 1 , 1)) to (1σ)ϑ (as this constant lower bounds (1σ)Σ(φ) when V (φ) ≥ ε 2 according to Lemma 4). In view of (A.8) and since Λ is continuously differentiable with respect to q, we can upper bound | ∂Λ(q) ∂q f (q)| by a constant ∆ Λ which depends on ∆, M , ε and the controller parameters. Hence, t 2t 1 is lower bounded by the time T which verifies ∆ Λ T = (1-σ)ϑ, thus T = (1-σ)ϑ∆ -1 Λ > 0. Two jumps are separated by T > 0: (2) holds with η(∆) = T . Item


Jonathan Passerat-Palmbach

Claude Mazel

David R C Hill

Pseudo-Random Streams for Distributed and Parallel Stochastic Simulations on GP-GPU

Keywords: Stochastic Simulations, GP-GPU, Pseudo-Random Number Generators, Random Stream Parallelization

Random number generation is a key element of stochastic simulations. It has been widely studied for sequential applications purposes, enabling us to reliably use pseudo-random numbers in this case. Unfortunately, we cannot be so enthusiastic when dealing with parallel stochastic simulations. Many applications still neglect random stream parallelization, leading to potentially biased results. In particular parallel execution platforms, such as Graphics Processing Units (GPUs), add their constraints to those of Pseudo-Random Number Generators (PRNGs) used in parallel. This results in a situation where potential biases can be combined with performance drops when parallelization of random streams has not been carried out rigorously. Here, we propose criteria guiding the design of good GPU-enabled PRNGs. We enhance our comments with a study of the techniques aiming to parallelize random streams correctly, in the context of GPU-enabled stochastic simulations.

Introduction

Stochastic simulations can be computationally greedy applications and the tendency is to always tackle bigger problems following the increase of computer performances. Domain experts, trying to decrease over and over again their simulation execution time, look for new solutions. For about five years now, we have seen a growing interest in General-Purpose computation on Graphics Processing Units (GP-GPU) through the simulation community. GP-GPUs offer a startling amount of computational power at an incredibly low cost compared to computing clusters and supercomputers. However, the introduction of this new kind of architecture implies not only the rewriting of the simulation algorithms, but also of the tools they use. In this article, we have chosen to focus on Pseudo-Random Number Generators (PRNGs), which are the basis of any stochastic simulation.

Sequential PRNGs have been studied for a long time L'Ecuyer (2010), and finding a good quality PRNG to use in a sequential application has not been a problem for more than a decade. Many works have been accomplished to characterize the statistical quality of PRNGs, leading to several testing libraries. Nowadays, reference testing suites are well known and are used to assess PRNGs. However, a PRNG should always be considered in relationship with the scope of the application it feeds. For instance, cryptographic applications developers ought to base their choice on the NIST testing battery [START_REF] Rukhin | A statistical test suite for random and pseudorandom number generators for cryptographic applications[END_REF], whereas simulationists should use TestU01 L' Ecuyer and Simard (2007). These two testing batteries can be considered as a standard for their respective domains at the time of writing.

According to [START_REF] Coddington | Random number generators for parallel computers[END_REF]; [START_REF] Hellekalek | Good random number generators are (not so) easy to find[END_REF], a PRNG should perform well on a single processor before being parallelized. Yet, statistical quality is a necessary but not sufficient condition when selecting a PRNG to use in a parallel context: indeed, parallel streams should be independent. Thus, providing high quality random numbers becomes even more difficult when dealing with parallel architectures. We have to take into account the parallelization technique: how will we partition random streams among parallel processing elements (threads or processors for instance)? How will we ensure the independence between parallel streams in order to prevent the simulations involved from producing biased results? The major problem concerning independence between random streams is that no mathematical proof exists to ensure it. However, some studies lay out well-known techniques to spread random streams through parallel applications [START_REF] Coddington | Random number generators for parallel computers[END_REF]; Hellekalek (1998a); [START_REF] Traore | The use of random number generation for stochastic distributed simulation: application to ecological modeling[END_REF]; [START_REF] Hill | Practical distribution of random streams for stochastic high performance computing[END_REF]; [START_REF] Reuillon | Parallel stochastic simulations with rigorous distribution of pseudo-random numbers with distme: Application to life science simulations[END_REF]. They try to ensure the maximum independence between random streams using different strategies. We will consider in this article whether or not, and how, these techniques can be implemented on GPU.

Apart from the parallelization technique, another point relies directly on the architecture where the involved stochastic simulations run. If we consider a GP-GPU environment, a new difficulty comes into play: harnessing the power of the device requires a rather good knowledge of GPUs. With recent programming frameworks like CUDA (Compute Unified Device Architecture) or OpenCL (Open Computing Language), almost anyone can develop applications for GPUs, but obtaining the announced performance gain implies a higher level of understanding. Most of the work and considerations exposed here rely principally on NVIDIA CUDA solutions. We are also working with the emerging OpenCL standard Khronos OpenCL Working Group (2010), but the latter is still not robust enough in our opinion. Its current performances are slower than what you could obtain with CUDA Karimi et al. (2010). However, we feel that this standard deserves our interest, and should take on an important part of our future work.

In this article, we will name as processing elements those effectively computing data in parallel. In a CUDA GP-GPU environment, threads will be regarded as these elements, since this framework relies on a thread level logic referred to as SIMT (Single Instruction, Multiple Threads). The latter abstracts a much more standard designation known as SIMD (Single Instruction, Multiple Data). GPUs are based on this kind of parallel architecture. Here, each thread must be given different data that will be computed by an identical operation. In the case of parallel stochastic simulations, we need to furnish an independent stochastic stream to each thread, in order to prevent potential biases that could be introduced otherwise. Now, the two subjects we are tackling in this study are:

• How are these random streams produced on GP-GPU?

• How can we ensure that they are independent?

These two main problems led us to share our experience of using stochasticity in GP-GPUenabled simulations. We focus on the parallelization techniques of pseudo-random streams used to directly feed parallel simulation programs running on GPU (called kernels in the CUDA language). Before stating what we will present in this paper, let us point out first, that we will not propose any new PRNG, and second, that our study is not tied to the parallelization of random number generation algorithms, albeit we do survey some parallel PRNG algorithms. The main point of this study is to give guidelines that will hopefully help developers to use reliable parallelization techniques of random streams to ensure their independence, and that are, in addition, well adapted to GPU architectures particularities. Some works have attempted to speed-up generation using the GPU before retrieving random numbers back onto the host. However, current CPU-running PRNGs display fairly good performances thanks to dedicated compiler optimizations. For instance, Mersenne Twister for Graphics Processors (MTGP) [START_REF] Saito | A variant of mersenne twister suitable for graphics processors[END_REF], the recent GPU implementation of the well-known Mersenne Twister Matsumoto and Nishimura (1998), is only announced as being 6 times faster than the CPU reference SFMT (SIMD-oriented Fast Mersenne Twister) [START_REF] Saito | Simd-oriented fast mersenne twister: a 128-bit pseudorandom number generator[END_REF], which is already very efficient in terms of performance. In previous studies, we showed that the time spent in generating pseudo-random numbers consumes at most 30% of CPU time for some "stochastic-intensive" nuclear simulations [START_REF] Maigne | Parallelization of monte carlo simulations and submission to a grid environment[END_REF], but they are very scarce. For less intensive simulations, when less than a billion numbers are needed, there is no real need for parallelization and unrolling is still the most efficient technique [START_REF] Hill | Urng: A portable optimization technique for software applications requiring pseudo-random numbers[END_REF]. Considering the small part of the execution time used by most stochastic simulations to generate random numbers, it is not worth limiting GPUs usage at the generation task. To use GP-GPUs at their full potential, we are more interested in providing random numbers to GPU-running applications that will consume them directly on the device.

In this study, we will:

• Propose GP-GPU specific criteria for PRNGs design;

• Survey random streams parallelization techniques of random streams in any distributed environment;

• Suggest requirements to implement parallelization techniques of random streams on GP-GPU;

• Study, according to the previously introduced requirements, the suitability of well-known PRNGs and random streams parallelization techniques for GP-GPU architectures.

Pseudorandom Number Generation on GPU

This section will survey the major propositions that can be found in the literature about PRNGs implementations on GPU platforms. Historically, a common way to deal with random numbers on GPU was to generate them on CPU before transferring them on the graphics processor. This solution had to face the well-known bottleneck of data transfer between the CPU host and the GPU device. Even with nowadays PCI Express 16X running at 8GB/s, this approach is considerably limiting the throughput of high performance applications.

Let us first focus on PRNGs implementations using GPUs as an hardware accelerator only, to provide random numbers quicker to host applications. Actually, such implementations can mostly be found during the genesis of GP-GPU, when the underlying architecture and programming languages features narrowed applications' scope. Until recently, designing a PRNG for GPUenabled platforms could be very tricky as it forced programmers to deal with graphics Application Programming Interfaces (APIs). Some implementations are presented in [START_REF] Sussman | Pseudorandom number generation on the gpu[END_REF]. The authors especially list the limitations of these GPU dedicated PRNGs due to the past weaknesses of the hardware. Limited output per thread or untruthful operations were part of the restrictions that made these PRNGs feeble for High Performance Computing (HPC) applications.

Since 2008 and the recent advances from NVIDIA, new GPU software and hardware architectures offer the precision and speed needed by many HPC applications. Langdon presents a minimal implementation of the standard Park Miller PRNG [START_REF] Park | Random number generators: Good ones are hard to find[END_REF] on an NVIDIA 8800 GTX GPU in its paper from 2008[START_REF] Kirk | A fast high quality pseudo random number generator for graphics processing units[END_REF]. He announces a speed up of more than 40 compared to his Intel 2.40 GHz CPU. One year later, he increased again the speed of his application by four [START_REF] Langdon | A fast high quality pseudo random number generator for nvidia cuda[END_REF] by using the NVIDIA CUDA technology NVIDIA (2010b) with a Tesla T10 GPU. Nevertheless, we do not advise the use of the old Park Miller generator that has many known flaws, though it was still employed until recently in some widespread networking simulation software [START_REF] Entacher | Pitfalls when using parallel streams in omnet++ simulations[END_REF].

CUDA has been designed to allow developers to easily harness the computation power of GPUs. In his first implementation, Langdon had to deal with a complex and unadapted graphics API. With CUDA, developers can program GPUs without wasting their time making algorithms and their data fit into graphics dedicated data structures, such as pixels shaders. Furthermore, CUDA does not propose a new programming language but only some C extensions, making it easier to learn for C familiars. CUDA-enabled graphics boards fulfil the requirements noted in the conclusion of the previously cited [START_REF] Sussman | Pseudorandom number generation on the gpu[END_REF], with for instance an implementation of the IEEE 754-2008 floating point numbers standard. The current generation of boards based upon the Fermi architecture is now proposing configurable L1 cache, ECC memory and a considerable increase of performance in double precision, while owning twice as much cores as the earlier mentioned T10 processor. At the time of writing, MTGP Saito (2010) is to our knowledge the sole parallel PRNG that has been specifically designed to run on GPU. The algorithm is intrinsically parallel, and targets the first goal of random number generation on GPU: speeding-up numbers output.

Recent GPU architectures, such as Fermi, opened new development perspectives. Having larger and faster memory areas available per thread, and being able to use object-oriented features, applications have become more and more ambitious, so have PRNGs implementations. GPUs are now considered as a fully capable platform, able to run entire applications by themselves. To do so, developers need PRNG implementations for GPU that allow their applications to directly consume the issued numbers. Such a tendency can be observed with the increasing number of available libraries for CUDA. Pseudo-Random numbers generation follow the same tendency, and we have noticed several contributions in the last two years.

We noted two main proposals in this domain: CURAND and Thrust::random. They both aim to provide a straightforward interface to generate random numbers on GPU. Introduced in the latest version, at the time of writing, of the CUDA framework, CURAND NVIDIA (2010a) has been designed to generate random numbers in a straightforward way on CUDA-enabled GPUs. The main advantage of CURAND is that it is able to produce both quasi-random and pseudorandom sequences, either on GPU or on CPU. The quasi-RNG and pseudo-RNG implemented are respectively Sobol and XorShift [START_REF] Marsaglia | Xorshift rngs[END_REF]. The API of the library stays the same no matter which kind of RNG is selected and the platform on which the application is run on.

Thrust::random is part of a GPU-enabled general purpose library called Thrust [START_REF] Hoberock | Thrust: A parallel template library[END_REF]. This open-source project intends to provide a GPU-enabled library equivalent to standard general-purpose C++ libraries, such as STL or Boost. Classes are split through several namespaces, of which Thrust::random is an example. The latter contains all classes and methods related to random numbers generation on GP-GPU. Thrust::random implements three PRNGs, each through a different C++ template class. We find a Linear Congruential Generator (LCG), a Linear Feedback Shift (LFS) [START_REF] Tausworthe | Random numbers generated by linear recurrence modulo two[END_REF] and a Subtract With Borrow (SWB) [START_REF] Marsaglia | A random number generator for pc's[END_REF]; [START_REF] Marsaglia | A new class of random number generators[END_REF]. Although the latter PRNG is mentioned as Subtract With Carry in Thrust::random documentation, Marsaglia's original proposition is known as SWB. In spite of the known flaws laid out by all these generators, the library offers simple ways to combine them into better quality randomness sources, like L'Ecuyer's Tausworthe combined generators L'Ecuyer (1996).

We recently proposed our own pseudorandom number generation toolkit for GPU named ShoveRand [START_REF] Passerat-Palmbach | Shoverand: a model-driven framework to easily generate random numbers on gp-gpu[END_REF]. It distinguishes from its counterparts by introducing a meta-model that enables the description of every PRNG characteristics. Moreover, this meta-model is implemented exclusively through C++ compile-time template mechanisms, thus introducing no overhead at runtime. ShoveRand offers a common API to users, whichever PRNG they select, and it also guides developers who would like to integrate a new generator to the framework. The latter performs compile-time analysis on the provided source code to ensure that only PRNG implementations which public interface matches our guidelines will compile successfully.

Although GPUs bring much more peak performances than CPUs, they must be carefully programmed to deliver the expected power, and most PRNGs have to be rethought from scratch to leverage GPUs characteristics. In fact, GPU architectures combine a manycore approach with SIMD vector cores. As vector processors do, GPU-enabled algorithms need to repeat the same operation on different data to correctly exploit the device. This is the main reason of the recent dedicated PRNGs proposals.

Other quality implementations can be found in [START_REF] Bradley | GPU Computing Gems Emerald Edition, chapter 16 -Parallelization Techniques for Random Numbers Generators[END_REF]. This study does not propose any new PRNG but details how a small set of reference quasi-random and pseudorandom number generators (Sobol, MRG32k3a and Mersenne Twister) have been successfully ported to GPU through CUDA. We draw your attention to the fact that these fine RNGs led to different GPU implementations, depending on their characteristics. For example, the small memory footprint of MRG32k3a allows an instance per thread whereas Mersenne Twister's large data structure conduct to a block of threads implementation level. We will detail the different options offered when implementing a PRNG on GPU in a further section, but the work of [START_REF] Bradley | GPU Computing Gems Emerald Edition, chapter 16 -Parallelization Techniques for Random Numbers Generators[END_REF] implicitly distinguishes two RNG groups: those which CPU code can directly be ported to GPU, with very few changes, and those who need to be redesigned to fit with GPU constraints.

Implementing PRNGs in a way to draw numbers directly on GPU led us to think about the best design of such pieces of software, considering both PRNG characteristics and GPU constraints.

Implementing PRNGs on GPUs

GP-GPU specific criteria for PRNGs design

As a result of the SIMD parallelism between threads and of their graphics processors legacy, GP-GPUs are not equivalent to a set of standard processors used in parallel. These particularities introduce some constraints that need to be satisfied if we do not want to see the overall simulation performance drop significantly. Thus, we will introduce in this section the requirements we find compulsory for a PRNG to run efficiently on a GP-GPU architecture.

The main goal targeted when using GP-GPUs is to obtain greater speed-ups. But this kind of device has not been primarily designed to support general computations and is more inclined to perform some arithmetic operations. Nowadays GPUs still display different performances with single and double precision floating point numbers, in accordance with the IEEE 754-2008 standard. For instance, the previous generation of NVIDIA supercomputing-dedicated GPU, the Tesla T10, was known to deliver ten times less computational power when dealing with double precision floating point numbers rather than with simple precision. Even if the current cutting-edge GPU, the NVIDIA Fermi T20, has considerably reduced the gap between these two precisions (a factor 2 still exists), it is wise to remain cautious before using double precision operations on GPU. Most of the time, single precision floats are sufficient enough to handle random values contained in [0 ; 1[ and should consequently be favoured. This proposition leads us to our first criterion: single precision floating point numbers should be preferred throughout the GP-GPU random number generation algorithm.

Another legacy of graphics processors is the heterogeneous memory organization. To complete what has previously been said on this subject, let us recall that several memory areas are reachable by threads running on a GPU. The capabilities of these memories, i.e. their capacity and response time, depend on two characteristics. First, the more threads can reach a memory area, the slower it is. Indeed, registers allocated to a single thread are the fastest memory this thread will be able to communicate with. Close to the same speed, we find shared memory, reachable by a relatively small amount of threads, all belonging to the same block, and so running on the same core of the GPU. On the other hand, every thread running on the GPU, regardless of which core they are located on, can access a wide memory area, commonly called global memory. This latter area is far slower than its previous counterparts (a few hundreds of GPU cycles are necessary for a basic global memory access). Second, read-only memories are faster since they can fully benefit from cache mechanisms, contrary to read-write memories. In the light of these memory constraints, it is obvious that GPU PRNGs should be designed with particular attention to sparing costly memory accesses. Commonly, static parameters will take place in read-only areas, whereas dynamic elements such as state vectors will be handled at the thread or thread group level. In a more formal way, the following is another criterion of good design: the algorithm should be designed in a way to avoid global memory accesses.

Taking into account the particularities of GPUs and their architecture, we have proposed two new requirements for PRNGs to run efficiently on such devices. They are summed up hereafter:

1. Single precision floating point numbers should be preferred throughout the GP-GPU random number generation algorithm;

2. The algorithm should be designed in a way to avoid global memory accesses.

Location of PRNGs' internal data structures on GP-GPU

We focus here on the implementation level of PRNGs on the GPU, which describes the memory area where the algorithm's data is located. In Passerat-Palmbach et al. ( 2010), we identified three implementation levels, mapped on the CUDA thread hierarchy: threads, blocks of threads and grid of blocks. These strategies directly impact the PRNG implementation, so as the parallelization technique coupled with it. Let us introduce them to understand the technical prerequisites about the subjects we are tackling in this study.

Obviously, new PRNG algorithms have to take advantage of GPU intrinsic properties. For simplicity purposes, we will briefly introduce the major concepts that rule GPU architectures, that is to say: heterogeneous memory hierarchy and thread organization. Notions described in this paragraph are represented in Figure 1. On the one hand, the thread organization is meant to maximize the performances of the application. Threads are bundled into blocks of threads to be assigned to one of the Streaming Multiprocessor (SM) of the GPU. SMs can mostly be considered as GPU equivalents to CPU cores. The important point is that they have their own thread scheduler, which champions threads which data are available. Selected threads then run on Streaming Processors (SP): the processing units of SMs. The matching notion of heterogeneous memories comes into play at this point. Threads can access several memories displaying various capacities, which we will discuss more thoroughly in a further section. For now, we just need to keep in mind the hierarchy of memory areas: i.e. the classification of memories based on their response time and visibility from threads. From the fastest to the slowest, threads can access: registers, shared memory, local memory and global memory. This enumeration does not take into account any constant memory since they cannot be written from kernel programs. Thus, they would not be able to store the produced random numbers. When registers and local memory are dedicated to a single thread, shared memory is visible to all the threads within a common block (inside a SM), while every thread can reach global memory.

This memory organization highly affects the PRNG's performances. With the first considered implementation level, using a generator per thread, the internal state of the PRNG has to be saved in each thread's local memory. Most of the time, internal states are formed by several elements contained in arrays. Now, CUDA related works, like [START_REF] Kirk | A fast high quality pseudo random number generator for graphics processing units[END_REF], specify that arrays declared within a thread are stored in the local memory. Although its name seems to indicate a thread scope, please note that local memory is actually a subset of global memory, and suffers consequently from the same slowness. This area is also allocated to threads when their register set is depleted, since registers are available in limited quantities on GPUs. At the time of writing, Fermi, the latest CUDA architecture has partially solved this problem, thanks to a larger amount of registers per thread, thus allowing small-memory-footprint PRNGs to store all their data in the register and shared memory space. Equivalently, with the third cited implementation level, a PRNG for all threads, that is to say a grid-level scope, the global memory is solicited to store the state of the unique PRNG of the application. Each thread draws a number and updates its component of the state in global memory. In [START_REF] Zhmurov | Efficient pseudo-random number generators for biomolecular simulations on graphics processors[END_REF], authors present three basic generation algorithms working either with a single instance of the PRNG for the kernel or with an instance per thread. The three algorithms exposed are quite basic: Ran2, Hybrid Taus and a Lagged Fibonacci generator. In the same way, [START_REF] Langdon | A fast high quality pseudo random number generator for nvidia cuda[END_REF] chooses to generate a number per thread in its GPU version of the Park-Miller algorithm. These two approaches make a heavy use of global memory. This has the advantage of being persistent across kernel launches within the same application. It is however important to realize that this area used to be quite slow: it implied a 400 to 800 clock cycle latency because it was not cached NVIDIA (2010b) in previous CUDA architectures. Once again, Fermi plays its part and enhances memory access time. However, the global memory approach will still display a worse latency than its counterparts.

The second implementation scope, the block of threads level, is the only one left to discuss. Every thread in a block can access a shared memory area. PRNG algorithms can consequently store their internal state in this space, enabling every thread of the block to update it. Shared memory is implemented on-chip and is consequently announced as fast as registers. Thus, PRNGs implemented at a block level will not suffer from the memory latency induced by slow global memory accesses. For independence purpose, we could imagine a variant of this block of threads scope, with a PRNG per warp.

The concept of Warps, as introduced by NVIDIA, corresponds to a subgroup of threads dynamically formed by the device at runtime: threads within a warp achieve memory accesses in parallel. Warps are thus the smallest GPU units that are able to process independent code sections. Indeed, given that different warps either run on different SMs, or on the same but at different clock ticks, they are fully independent to each other. Because of memory latency, warps-schedulers select the warps that have their data ready to process. Consequently, the more warps can be scheduled, the better the memory latency can be hidden.

In order to implement a PRNG at a given level, the following requirements must be met: first we need a common memory area accessible by every member of the group. In the case of warps, the shared memory area assigned to their belonging block will perfectly suit. Second, in order to build their own random sequence, processing elements need to be able to distinguish their corresponding PRNG. There is no problem to do so when dealing with a PRNG implemented at either thread, block or kernel level, since CUDA provides us a way to uniquely identify each of the aforesaid element. Although warps identifiers are not directly available through CUDA keywords, we have shown in another study how a thread could fill a variable with the value of its parent warp's identifier [Passerat-Palmbach.etal.2011b].

The three main implementation scopes detailed in this section are sketched in Figure 2:

!"#$%&"'()" *+)(,-" !"#$%&"'()" *+)(,-" !"#$%&"./)",00"*+)(,-1" !"#$%&"'()"20/34" /."*+)(,-1"

!"# As a conclusion, no matter which strategy we choose to implement random number generation facilities on GPU, we will always have to deal with distribution techniques. Such techniques could be hidden in a well-designed library, or directly applied by the simulation developer.

$ %# $ &# $ '# $ (# $ )# $ *# $ +# $ ,# !

Requirements for distribution techniques of random streams on GP-GPU

GP-GPU specific requirements for random streams parallelization

The literature is full of references describing the profile of what a good usage of parallel PRNGs should be. For example, [START_REF] Coddington | Random number generators for parallel computers[END_REF]; [START_REF] Hellekalek | Good random number generators are (not so) easy to find[END_REF] list requirements that any sequential or parallel PRNG should meet. GPUs are particular parallel architectures, so any PRNG running on this kind of device should, at least, match the requirements enumerated in the previous references. In this section, we will successively check how these criteria can be adapted to GPU architectures.

Emphasizing parallel PRNG performances, [START_REF] Coddington | Random number generators for parallel computers[END_REF] noticed that each processor should generate its sequence independently of the other processors. We consider, indeed, that every processing element should have its own stochastic stream at its disposal. This condition must be satisfied first, not only for efficiency, but especially because GPU-enabled stochastic simulation parallelization principles rely on it. First, it is a necessary, but not sufficient, condition to fulfil to ensure a higher independence of stochastic streams feeding different replications of a simulation. Second, considering a single replication, the SIMT parallelism level leads threads to compute their own data sets, including their own stochastic stream. Thus, our first requirement concerning random streams parallelization can be expressed as follows: each thread should dispose of its own random sequence.

As we explained previously, GP-GPU programming frameworks offer a thread scope rather than a processor one. The threads in use for GP-GPU propose an abstraction of the underlying architecture. They are concurrently running on the same device and handle their own local memory area. Thread scheduling is at the basis of GPU performances. Memory accesses are the well-known bottleneck of this kind of device. Indeed, running a large amount of threads in turn allows GPUs to bypass memory latency. There should always be runnable threads while others are waiting for their input data. Disregarding the effective number of processors, we theoretically say that the more threads you have, the better your application will leverage the device. Applications need to be written to use the maximum number of threads, but also to scale up transparently when the next GPU generation will be able to run twice as many threads as today. So, in accordance with [START_REF] Coddington | Random number generators for parallel computers[END_REF] who advocates that the generator should work for any number of processors, our second GP-GPU specific requirement for parallelization techniques of random streams is that it must be usable for any number of GP-GPU threads.

Returning to the original requirements, we then find [START_REF] Coddington | Random number generators for parallel computers[END_REF] states that parallel random streams produced should be uncorrelated. This criterion is related to both PRNG intrinsic properties and to the parallelization technique set up. We previously stated that a PRNG candidate to parallelization should first perform well on a single processor. Thus, we will not take its intrinsic qualities into account here. However, no matter the worth of the used PRNG, the parallelization techniques must be used carefully, as it will be shown in the next section. Please note that this requirement is neither affected by GP-GPU architectures nor by programming frameworks. As a result, we will just recall it without modifying its expression. [START_REF] Coddington | Random number generators for parallel computers[END_REF] also noted that the same sequence of random numbers should be produced for different numbers of processors, and for the special case of a single processor. Here, we understand that the PRNG output on each processor should not depend on the number of processors used. This point is very important and must be treated carefully when choosing a parallelization technique. For example, in a distributed environment containing several processors, a scheduler can govern execution. Depending on the scheduler algorithm and on the global system charge, parallel executions of different parts of a simulation might not execute in the same order. In such a case, it is compulsory for the PRNG output to be independent from the order in which simulations parts may run. If this requirement is not met, reproducibility of simulations executions is no longer ensured. We can do this for games, but not for scientific applications. Reproducibility is needed when dealing with stochastic simulations, in order to debug a problematic case raised by a particular random stream for instance. We also think about design of experiments for simulations, where reproducibility is mandatory to isolate the impact of parameters variations on results.

In the case of GP-GPU, we find exactly the same problem at the thread level. These entities are also scheduled, not atomically but by bundles. Fortunately, both threads and their bundles own a unique identifier allowing us to distinguish them among executions. Thus, if a parallel random stream is only bound to the unique identifier of a thread, according to our first requirement, output will be reproducible through multiple executions. Therefore, we can obtain the necessary bijective relation between a T i thread and an SS i stochastic stream, as stated in Figure 3: Finally, we can write our last requirement in such a way: when the status of the PRNG is not modified, the sequence of random numbers generated for a given thread must be the same no matter the number of threads and no matter of threads scheduling.

! "# ! $# %# ! &# %# ! '# (( "# (( $# %# (( &# %# (( '#
Let us now sum up the requirements targeting GPUs we highlighted in this part, continuing the enumeration started in the previous section: 3. Each thread should dispose of its own random sequence; 4. The parallelization technique must be usable for any number of GP-GPU threads; 5. The parallel random streams produced should be uncorrelated; 6. When the status of the PRNG is not modified, the sequence of random; numbers generated for a given thread must be the same no matter the number of threads and no matter of threads scheduling.

Random streams parallelization techniques fitting GP-GPUs

This section presents how the main techniques used to distribute random streams between processing elements can be adapted to GPU architectures, depending on their ability to fulfil the previously introduced requirements.

Sequence Splitting

The Sequence Splitting (SS) method is also known as "Blocking" or "Regular Spacing". It consists in allocating non-overlapping, contiguous and equally sized blocks from the original random stream to form substreams. When partitioning a sequence (x i , i = 0, 1, . . .) into N streams, the j th stream is x k+(j-1)m , k = 0, . . . , m -1, where m is the length of each stream; m must be chosen so that each stream is long enough to achieve the stochastic simulation performed by the corresponding process. For instance in [START_REF] Hechenleitner | Defects in Random Number Routines of Well-Known Network Simulators and Appropriate Improvements[END_REF], Hechenleitner showed that in the OMNeT++ simulation software, the spacing between sequences set to 1 million draws led to biased results (due to inter-sequence correlations) for processes using more random numbers.

This technique implies the knowledge of how many numbers each thread will consume at most. Indeed, knowing that each thread consumes at most L random numbers, then the first L numbers will be attributed to the first thread, the L following to the second thread and so forth. Following the previous formalism, we have Y i = {X iL , X iL+1 , . . . , X (i+1)L-1 }. This numbers repartition is described in Figure 4.

!"#$%&'(')' !"#$%&'('*' !'#+,-"./' 0 1' 0 2' 0 3' 0 4' 0 )' 0 5' 0 *' 0 6' 0 7' 0 8' 0 21' 0 22' 0 1' 0 2' 0 3' 0 4' 0 )' 0 5' 0 *' 0 6' 0 7' 0 8' 0 21' 0 22' Efficient Sequence Splitting relies on a particular feature of the PRNG called Jump Ahead, or Skip Ahead. Here, we discern two categories of algorithms. Some PRNGs contain an algorithm able to perform Jump Ahead, thanks to intrinsic properties of the PRNG Haramoto et al. (2008). This allows us to reach any part of the sequence in equal time, regardless of the destination point. Let this technique be considered as Intrinsic Jump Ahead hereafter. The other solution is to emulate Skipping Ahead: to do so, we have to compute an advanced state by processing step by step the previous ones (for example, starting from X iL , X (i+1)L can be computed by running the PRNG L times in a sequential way in order to get X iL+1 , ..., X (i+1)L-1 and finally X (i+1)L ). In fact, whichever PRNG you use, you can unfold the sequence to the desired point, and store the state vector at this point in order to be able to load it later. Such a vector is named Seed Status in Passerat-Palmbach et al. ( 2010), since it is able to set a generator in a predefined state. Emulated Jump Ahead can become very costly though: indeed, the further you need to go in the sequence, the more time it takes to compute the Seed Status.

When an intrinsic Jump Ahead algorithm is available for the involved PRNG, Sequence Splitting is a very good approach for GP-GPUs. However, as far as we know there are few GP-GPU ports of algorithms with Jump Ahead features. At the time of writing, we are only aware of an MRG32k3a implementation detailed in [START_REF] Bradley | GPU Computing Gems Emerald Edition, chapter 16 -Parallelization Techniques for Random Numbers Generators[END_REF] and of the recent Tiny Mersenne Twister (TinyMT) [START_REF] Saito | Tiny mersenne twister (tinymt)[END_REF], available to download but not described in any scientific paper yet. On the other hand, Emulated Jump Ahead is not GPU-compliant because statuses computation is a purely sequential operation (we need X n to compute X n+1 ). As a consequence, different threads may lead to different computation time to process the next state to jump to, because of their sequential execution. Thus, the SIMD parallelism would be shrunk, making the overall gain decrease. To solve this problem, we propose to pre-compute substreams on the host side, store the Seed Status at each substream starting point and then transfer all these statuses to the device.

Random Spacing

The Random Spacing (RS) or Indexed Sequences (IS) method builds a partition of N streams by initializing the same generator with N random statuses. In the case of old LCGs it was named random seeding. For modern generators with a more complex status, the random statuses are generated with another RNG, and this technique is interesting when generators have a huge period. This technique is easy to set up. In [START_REF] Wu | Parallel use of multiplicative congruential random number generators[END_REF], the authors have given the minimum distance between N sub-sequences which status were randomly chosen: it is in average equal to 1/N 2 multiplied by the period length. More precisely, the probability of overlapping between N sequences of length L, issued from a PRNG of period P is equal to: 1 -(1 -N L/(P -1)) (N -1) . That is equivalent to N (N -1)L/P when N L/(P -1) is in the neighbourhood of 0.

Overlapping risks become sizeable with short period PRNGs. All the PRNGs tested by Pierre L'Ecuyer and Richard Simard in L'Ecuyer and Simard (2007) display periods P from 2 24 to 2 131072 . Most of them have log 2 P of the order of a few dozens. Now, given that nowadays longest simulations can consume up to several hundreds of billions of random numbers, (L = 10 11 ), a hundred of such replications (N = 100) leads to a N (N -1)L of approximately 10 15 , i.e. 2 50 . The probability to see an overlapping between two sub-sequences issued from a PRNG of period P far bigger than 2 50 is negligible.

Random Spacing initialization process consists to draw a random number (from another generator), and to set it as the seed of the considered PRNG. Consequently, it fits GP-GPUs well, since this operation can be done in parallel without any constraint. In Figure 5, we have sketched the use of Random Spacing to issue a random stream assigned to each of the 3 threads represented. Yet, the risk of overlapping between sub-sequences must be evaluated according to the amount and the length of the sub-sequences and to the period of the PRNG used. If we select generators with large periods, such as WELLs and Mersenne Twister, this risk is negligible.

!"#$%&'()" * +" * ," * -" * ." * /" * 0" * 1" * 2" * 3" * 4" * ,+" * ,," * ,-" * ,." * ,/" * ,0"

Leap Frog

Leap Frog (LF) is the way to partition a random stream like a deck of cards. Random numbers are allocated in turn to processors like cards would be dealt to players. Pragmatically, let each processor hold an i identifier. Every such PE will build a Y i substream from an X original random stream such as Y i = {X i , X i+N , . . . , X i+kN }, with N equal to the number of processors. This technique is not quite adapted to split random streams on GP-GPU since it does not satisfy the last constraint expressed in the previous section. In fact, if the number of threads changes, the subsequence assigned to each thread will be different. This situation is shown in Figure 6. Now, the number of threads for an application is bound to the underlying device: GPUs can run a different number of threads concurrently, depending on their architecture generation. As a result, we would not be able to ensure the reproducibility of a simulation from a GPU to another, even if we initialized the PRNG with the same parameters. Furthermore, some bugs induced by random draws might not appear on the developer's device, while they would on the user's.

The solution would be to implement the PRNG at a level with a constant number of threads. The CUDA framework handles constant-sized bundles called warps, which always contain 32 threads at the time of writing. They are in fact a subdivision of blocks of threads, and access consequently the same shared memory area. Nowadays, a great number of GPUs do not have enough shared memory to store a PRNG status per warp. However, the newest GPU generations offer larger shared memory spaces that could potentially enable us to use Leap Frog following this idea.

!"#$%&'()" * +" * ," * -" * ." * /" * 0" * 1" * 2" * 3" * 4" * ,+" * ,," -"56('78)" /"56('78)" * +" * ," * -" * ." * /" * 0" * 1" * 2" * 3" * 4" * ,+" * ,,"

Parameterization

Techniques presented so far tried to split a single stream into several substreams. Another approach consists in using several declinations of the same PRNG: each generator has the same structure and generation mechanism with a unique parameter set, called Parameterized Status hereafter.

Although no mathematical proof can establish this independence, some implementations of Parameterization are safe according to the current state of the art [START_REF] Mascagni | Parameterizing parallel multiplicative lagged-fibonacci generators[END_REF]. We especially think to the Dynamic Creator (DC) algorithm [START_REF] Matsumoto | Dynamic creation of pseudorandom number generators[END_REF] coming along with most of the generators from the Mersenne Twister family. DC integrates a unique identifier, which belongs to the Parameterized Status of the PRNG. This identifier becomes a part of the characteristic polynomial of the matrix that defines the recurrence of the PRNG. Two identifiers will consequently lead to two different Parameterized Statuses. Furthermore, DC ensures that the characteristic polynomials we obtain are mutually prime, and the authors assert that the random sequences generated with such distinct Parameterized Statuses will be highly independent, even if, as mentioned before, this fact cannot be mathematically proven.

This technique displays some constraints making it difficult to port to GPU. Unfortunately, few PRNGs propose it intrinsically. Some, such as LCGs, are even reckoned as bad candidates for Parameterization De [START_REF] De Matteis | Parallelization of random number gemerators and long-range correlations[END_REF]. In addition, storing Seed Statuses (basically the common seed given by the user to initialize a generator, or the internal state vector of the PRNG) being already problematic, we can scarcely imagine spending vast amounts of memory to store a Parameterized Status per thread. A Parameterization example is proposed in Figure 7:

!"#$%&'()" *+(+%','(-.'/"" *012"34" 5 6" 5 4" 5 7" 5 8" 5 9" 5 :" 5 ;" 5 <" 5 =" 5 >" 5 46" 5 44" *+(+%','(-.'/"" *012"37" ? 6" ? 4" ? 7" ? 8" ? 9" ? :" ? ;" ? <" ? =" ? >" ? 46" ? 44" *+(+%','(-.'/"" *012"38" @ 6" @ 4" @ 7" @ 8" @ 9" @ :" @ ;" @ <" @ =" @ >" @ 46" @ 44"

Figure 7: Three parameterized PRNGs producing three highly independent random sequences Every technique introduced so far, presents advantages and drawbacks. Most of them are related to the chosen PRNG. Depending on the application and environment you own, you might be forced to select a PRNG knowing it has some flaws in particular cases. In this way, Table 1 states PRNG kinds and parallelization techniques that work well together:

Conclusion

We have seen that more and more applications, and especially stochastic simulations, tend to take advantage of recent GP-GPU architectures in order to improve their performance. However GPU computing needs to offer the same tools as other platforms. High quality PRNGs belong to this category and have existed for more than a decade, although some recent publications dealing with GP-GPU implementations of PRNGs still propose old and weak generators. In this paper we have shown how difficult it could be to obtain good quality pseudo random sequences on GP-GPU. Indeed, it implies taking into consideration two different domains: GP-GPU programming and PRNG parallelization techniques. Issuing a PRNG that can produce independent stochastic streams when used in parallel is a first hurdle that not all PRNGs can get over. When you have at your disposal one that fulfils this requirement, it has to be ported to GP-GPU. It means that you need to think about a GPU implementation of your PRNG, if it is not already available.

This paper introduced the main problems that you will be faced with when trying to port your stochastic simulations to GP-GPU. To avoid these difficulties, and above all, errors and performance drops that could result from the use of a hazardous GPU-enabled PRNG, we first proposed PRNGs criteria dedicated to GPUs. In order to take advantage of the existing PRNG parallelization techniques on GPU, we first defined a set of requirements that should be met by any PRNG implementing a distribution technique on GPU. Finally, we studied the applicability of widespread distribution techniques on GPU, and determined which techniques best matched GPU constraints. These requirements and chosen techniques sum up the experience accumulated in our research team concerning GPU-enabled stochastic simulations.

We encountered lots of parameters brought up by PRNGs and GPU programming. As long as it can dramatically impact both the overall performance of the simulation and the quality of its results, it might be a good point to propose a straightforward API to use well-defined PRNGs on GPUs. In this way, libraries laid out in this paper represent in our mind an elegant manner to do so. They allow every user to easily call PRNG functions without wasting time on how parameters should be set. In the same way, they enable advanced users to adjust parameters in their own fashion. As a matter of fact, we introduced a recent framework named ShoveRand, which embed PRNGs following the presently established requirements into a GPUenabled library with a unified API. This library permanently evolves to integrate new PRNGs, such as the recent TinyMT from Mutsuo Saito, but also features proposed by the C++ Technical Report 1 to enhance random number generation.

Figure 1 :

 1 Figure 1: Simple representation of the major elements of a GPU

Figure 2 :

 2 Figure 2: PRNGs implementation scopes and their location in the different memory areas of a GPU

Figure 3 :

 3 Figure 3: Bijective relation between threads and stochastic streams

Figure 4 :

 4 Figure 4: Two random streams parallelizations based upon Sequence Splitting with two different sub-sequences lengths

Figure 5 :

 5 Figure 5: Random Spacing creation of three sub-sequences of equal length but differently spaced from each other

Figure 6 :

 6 Figure 6: Different threads numbers leading to different random substreams through the Leap Frog method

Table 1 :

 1 Summary of the potential PRNG/Parallelization technique associations

	Technique	Preferred	PRNGs
		PRNGs	to avoid
	Leap Frog	None
		(disable reproducibility) Linear generators
	Sequence Splitting	Intrinsic Jump-Ahead	Emulated
		compliant	Jump-Ahead
	Random Spacing	Large period	Short period
	Parameterization	MT family	LCG

partially funded by the Auvergne Regional Council and the FEDER