Keywords: Parallel Random Numbers, Parallel random streams, High-Performance Computing, Random Number Generation, Stochastic Simulation

There is an increasing interest in the distribution of parallel random number streams in the high-performance computing community particularly, with the manycore shift. Even if we have at our disposal statistically sound random number generators according to the latest and thorough testing libraries, their parallelization can still be a delicate problem. Indeed, a set of recent publications shows it still has to be mastered by the scientific community. With the arrival of multi-core and manycore processor architectures on the scientist desktop, modelers who are non-specialists in parallelizing stochastic simulations need help and advice in distributing rigorously their experimental plans and replications according to the state of the art in pseudorandom numbers parallelization techniques. In this paper, we discuss the different partitioning techniques currently in use to provide independent streams with their corresponding software. In addition to the classical approaches in use to parallelize stochastic simulations on regular processors, this paper also presents recent advances in pseudo-random number generation for general-purpose graphical processing units. The state of the art given in this paper is written for simulation practitioners.

Introduction

With the arrival of multi-core and manycore processor architectures on the scientist desktop, modelers who are non-specialists in parallelizing stochastic simulations need help and advice in distributing rigorously their experimental plans and replications according to the state of the art in pseudo-random numbers partitioning techniques. New programming models and platforms are now available to a larger audience, particularly with the arrival of very powerful general-purpose graphics processing units (GP-GPUs). Stochastic simulation is now an essential tool for many research domains. Many applications are concerned, particularly when dealing with complex system science. Random number generators (RNGs) are mainly used in simulation to model the stochastic phenomena that appear in the formulation of a problem or to apply an iterative simulation-based problem solving technique (such as Monte Carlo simulation [START_REF] Gentle | Random Number Generation and Monte Carlo Methods (Statistics and Computing)[END_REF] or coupling between simulation and operations research optimization tools. Pseudo-random number generators (PRNGs) use deterministic algorithms to produce sequences of random numbers and have been studied in pertinent reviews, many from Michael Mascagni and Pierre l'Ecuyer, and also from other authors, including the well-known first edition of Knuth's Art of Computer Programming, volume 2 and its re-editions [START_REF] James | A review of pseudorandom number generators[END_REF]; [START_REF] Knuth | The Art of Computer Programming[END_REF]; [START_REF] Lazaro | Fully 3D monte carlo reconstruction in SPECT: a feasibility study[END_REF]; [START_REF] Mascagni | Parallel linear congruential generators with sophie-germain moduli[END_REF][START_REF] L'ecuyer | An object-oriented randomnumber package with many long streams and substreams[END_REF]. If some applications can cope with 'bad' or poor randomness according to current standards, we cannot afford producing biased results for nuclear physics or medicine for instance El [START_REF] El Bitar | Fully 3D monte carlo image reconstruction in SPECT using functionalregions[END_REF]; [START_REF] Li | Improving performance via computational replication on a large-scale computational grid[END_REF]; [START_REF] Maigne | Parallelization of monte carlo simulations and submission to a grid environment[END_REF]; [START_REF] Reuillon | Testing 65536 parallel pseudo-random number streams[END_REF]. For parallel stochastic simulations, there is an even more critical need for a large number of parallel and independent random sequences, each with good statistical properties (approximating as close as possible a truly random sequence). We cannot give a mathematical proof of independence between two random sequences, so each time this word is used, it should be considered as (pseudo) independence according to the best accepted practices in mathematics and statistics.

A category of random numbers called quasi-random numbers can improve the convergence speed (and accuracy) of Monte Carlo integration and also of some Markovian analysis [START_REF] Niederreiter | Quasi-Monte Carlo Methods[END_REF]. However, these numbers are not independent and can only be used for specific applications. Random numbers produced by deterministic algorithms are called pseudo-random numbers by simulationists and by most researchers, whereas this notion of 'pseudo'-random refers to polynomial-time unpredictability properties in the area of cryptology. In this paper, we do not consider RNGs for cryptographic usage but only for scientific simulations. Truly random numbers could be a very interesting source, but they are not reproducible L' Ecuyer et al. (2002). They can be supplied for instance by a physical external source such as a radioactive decay, but there are also many drawbacks in their use. True random numbers can be clumsy to produce and use, some devices are subject to partial breakdowns (their production has to be regularly statistically checked), some have biases that have to be corrected, and most of the time, we need to store the sequences to be able to reproduce exactly the same sequence several times. This is not possible for many high-performance computing (HPC) applications. New devices named Quantis, proposed by 'Id Quantique', have generated a strong interest over the last few years. Their principle is based on an optical quantum process as source of randomness, and they are fairly easy to use through PCI cards or USB sticks. This source can deliver a sequence of bits at a high bit rate of 4-16 Mbits/s. This sequence is then post-processed by the Quantis processing unit to remove the potential biases. Although for a supercomputer we could consider the acquisition through many devices, this is not possible with world wide computing grids (such as the European Grid Initiative, EGI) with heterogeneous computing centers. Once again, the major problem in our opinion is that sequences are not reproducible. For small-scale simulations, we can store the sequences or use a central server approach (detailed in Section 3.1.1) similar with [START_REF] Stevanovic | Quantum random bit generator service for monte carlo and other stochastic simulations[END_REF]. Still, HPC applications, such as nuclear medicine simulations for instance, cannot download, or store the hundreds of thousands of billions of numbers needed to reproduce the experiments El [START_REF] El Bitar | Fully 3D monte carlo image reconstruction in SPECT using functionalregions[END_REF]; [START_REF] Li | Improving performance via computational replication on a large-scale computational grid[END_REF]; [START_REF] Maigne | Parallelization of monte carlo simulations and submission to a grid environment[END_REF]; [START_REF] Reuillon | Testing 65536 parallel pseudo-random number streams[END_REF]. If we consider smaller needs in terms of numbers of draws, optimizations using memory mapping and unrolling of pre-stored numbers can however be very efficient for various applications [START_REF] Hill | Object-oriented pattern for distributed simulation of large scale ecosystems[END_REF]. Reproducing experiments is the essence of science, and even if this is not always necessary, in the case of stochastic simulation, the need for reproducible generators is essential for various purposes, including the analysis of results. To investigate and understand the results, we have to reproduce the same scenarios and find the same confidence intervals every time we run the same stochastic experiment. When debugging parallel stochastic applications, we need to reproduce the same control flow and the same result to correct an anomalous behavior. Reproducibility is also necessary for variance reduction techniques, for sensitivity analysis and many other statistical techniques [START_REF] Kleijnen | Statistical tools for simulation practitioners[END_REF]. In addition, for rigorous scientific applications, we want to obtain the same results if we run the application in parallel or in serial. Consequently, software random number generation remains the prevailing method for HPC, and we will see that specialists are warning us to be particularly careful when dealing with parallel stochastic simulations De Matteis and Pagnutti (1990b); [START_REF] Hill | Good random number generators are (not so) easy to find[END_REF]; [START_REF] Pawlikowski | Speeding up stochastic discrete-event simulation[END_REF]; [START_REF] Pawlikowski | On automatic partitioning, run-time control and output analysis methodology for massively parallel simulations[END_REF].

In this paper, we focus our contribution on the following aspects:

• We point out that parallelization of random streams has still to be mastered by the scientific community. Although a set of fine techniques exists, recent publications show that more attention has to be given to this aspect (Section 2).

• We survey the parallelization methods and their design, and we propose a possible classification for simulation practitioners (Section 3). Reliable initializations of PRNGs are also discussed.

• In Section 4, we present the major PRNG parallelization software and libraries, and then we discuss the empirical testing approaches (Section 5).

• In the last 2 sections, we present considerations raised by the emergence of hybrid computing with GP-GPUs.

Parallel Stochastic Simulations and Random Numbers Generation

Stochastic simulations can require a huge computational capacity particularly when we design experiments to explore large parameter spaces [START_REF] Kleijnen | Statistical tools for simulation practitioners[END_REF]. New technologies in HPC and networking, including hybrid computing, have very significantly improved our computing capacities. The major consequence is the increased interest for parallel and distributed simulation [START_REF] Fujimoto | Parallel and Distributed Simulation Systems[END_REF]. Major problems arise when dealing with model partitioning and with the stochastic aspects, but in many cases, effective solutions can be developed. Communications protocols (synchronous simulation versus asynchronous) have to be examined seriously to avoid deadlocks and to preserve the causality and determinism principles. When dealing with stochastic parallel simulation, we have to mainly consider two aspects: the quality of the PRNG and the parallelization technique. If stochastic sequential simulations always require a statistically sound PRNG, in the case of parallel simulations, this requirement is even more crucial. As a consequence, simple linear congruential generators (LCGs) should be definitively banned from modern HPC applications L' Ecuyer et al. (2002). The parallelization technique should ideally generate pseudorandom numbers in parallel, that is, each processing element (PE) should autonomously obtain either its own random sequence or its own subsequence of a global sequence (partitioning of a main sequence). If such independence is not guaranteed, the parallelism is affected. The sequences assigned to each PE must not depend on the number of processors: if the latter varies, we cannot reproduce the simulation experiments. For example, in the case of the central server approach, the PE access to the central RNG in an order that relies on scheduling policies and network communications, which does not ensure the same PE will have the same numbers every time. Such a configuration disables reproducibility, and it also often creates a bottleneck. Designers of parallel stochastic simulations always have to reply to this fundamental question: how can we make a safe RNG repartition to keep, on the one hand, efficiency, and on the other hand, a sound statistical quality of the simulation to obtain credible results? Indeed, the validation of such parallel simulations is a critical issue. Paul Coddington precisely states: 'Random number generators, particularly for parallel computers, should not be trusted [START_REF] Coddington | Techniques for empirical testing of parallel random number generators[END_REF]. Much research has been undertaken to design good sequential RNGs. Whatever the parallelization technique is, we have to rely on a 'good' generator according to a set of main principles proposed by specialists [START_REF] Coddington | Techniques for empirical testing of parallel random number generators[END_REF]; Hellekalek (1998a); [START_REF] Lagarias | Pseudorandom numbers[END_REF] (1988,1995); [START_REF] Durst | Using linear congruential generators for parallel random number generation[END_REF]; [START_REF] Eddy | Random number generators for parallel processors[END_REF]; [START_REF] Entacher | Parallel random number generation: long-range correlations among multiple processors[END_REF]; [START_REF] Percus | Random number generators for MIMD parallel processors[END_REF]; [START_REF] Srinivasan | Testing parallel random number generators[END_REF].

Assessing the quality of random streams remains a hard problem, and many widely used partition techniques have been shown to be inadequate for some specific applications. For instance, recent studies in networking and telecommunication question the relevance of many stochastic simulations because of the use of poor quality RNGs [START_REF] Hechenleitner | Defects in Random Number Routines of Well-Known Network Simulators and Appropriate Improvements[END_REF]; [START_REF] Pawlikowski | On automatic partitioning, run-time control and output analysis methodology for massively parallel simulations[END_REF]; Pawlikowski (2003a), as well as the mediocre parallelization techniques used in some software [START_REF] Entacher | Pitfalls when using parallel streams in OMNET++ simulations[END_REF]. The next section presents the main advantages and drawbacks of partitioning techniques.

Design of Parallel and Distributed Random Streams

There are many approaches to obtain parallel random numbers streams either by partitioning the main sequence (stream) of a given generator into subsequences (substreams) or by parameterizing one generator to have several sequences (multiple streams). Different techniques have been surveyed in [START_REF] Traore | The use of random number generation for stochastic distributed simulation: application to ecological modeling[END_REF] and more recently in [START_REF] Bauke | Random numbers for large-scale distributed monte carlo simulations[END_REF]. This section tries to give a thorough overview of the distribution techniques found in the literature. Here, we propose taxonomy of distribution techniques, depending on the type of original random streams they target. We also try to identify the scope of each of these techniques, so that any practitioner can refer to this section to figure out which technique best fits his application.

Partitioning a unique original stream

Any PRNG whose state vector is n-bits wide generates, from an initial state, a periodic random number sequence whose period is inevitably less than 2 n . Random numbers distribution techniques introduced in this section share a common principle referred to as cycle division in [START_REF] Mascagni | High-performance monte carlo tools[END_REF]. It splits random numbers from the considered cycle between the different PEs.

Techniques detailed hereafter differ from each other only by the method they use to split the main stream. Basically, we consider four main partitioning techniques.

Central server technique

This first approach is not truly parallel. It consists in using a central server, running an RNG and providing on-demand pseudo-random numbers to different PEs. This approach displays two major drawbacks. First, a simulation using this technique will not be reproducible because of scheduling policies that might change the order in which numbers will be provided to PEs. This is very problematic when trying to debug a simulation because the same numbers will not necessarily be assigned to the same PE at each run. Moreover, the results of a simulation cannot be checked by other scientists for the same reason: there are few chances that the same random sequence will be issued twice. Second, the central server approach will create a bottleneck if too many PEs are considered. As a consequence, this is suitable for serious games with limited parallelism but not for scientific applications.

Sequence splitting

The sequence splitting 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 substream. 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 | On shortcomings of the NS-2 random number generator[END_REF], Hechenleitner showed that in the OMNeT++, the spacing between sequences set to 1 million draws led to biased results (due to inter-sequence correlations) for processes using more random numbers. However, the determination of a good value for m is not the only difficulty. If overlapping can be easily avoided, long-range correlations in the initial RNG can lead to small-range correlations between the potential substreams De Matteis andPagnutti (1988, 1990a). The impact of this kind of correlation is problematic as shown in [START_REF] Reuillon | Testing 65536 parallel pseudo-random number streams[END_REF]; [START_REF] Srinivasan | Introduction to parallel RNGs[END_REF]. Figure 1 represents an original stream chunked through sequence splitting for 2 and 3 PEs. The schema also insists on the fact that an original stream is split in equally sized parts whose length may vary from one application to another.

Figure 1: Sequence splitting in a unique original stream considering 2 processing elements (PEs), then 3 PEs

Random spacing

The random spacing or indexed sequences 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 2006, Wu and Huang [START_REF] Wu | Parallel use of multiplicative congruential random number generators[END_REF] showed that the minimum distance between n statuses generated in this way is on average 1/n 2 multiplied by the period length. The risk is of course to have a bad initialization linked to the fact that two random statuses could be too close to each other, implying an overlapping of corresponding sequences.

For a PRNG with a period P , the probability that n sequences of length L, generated by a random spacing technique will overlap is equal to 1 -(1 -nL/(P -1)) n-1 , which is equivalent to n(n -1)L/P when nL/P is in the neighborhood of 0. A situation implying 3 PEs is sketched in Figure 2 where we can notice the non-equal gaps between two random numbers ranges. 2004), (L = 10 12 ), a hundred of such replications (n = 100) makes n(n-1)L on the order of 10 16 (i.e., more than 2 53 ). The probability to see an overlapping between two subsequences issued from a PRNG of period P far bigger than 2 53 becomes negligible. Nonetheless, from 20 LCGs PRNGs tested in L'Ecuyer and Leydold (2005), Kleijnen1986 laid out an overlapping probability greater than 99.9% (period P such that log 2 (P ) < 50.4). Widely spread PRNGs, such as Comblec88 of period P = 2 61 (combined LCGs), are on the acceptance borderline for this technique (with n = 100 and L = 10 12 , the overlapping probability is equal to 0.43%). They are shipped with several renowned software packages though, including RANLIB, CERNLIB, Boost, Octave and Scilab L' Ecuyer and Leydold (2005). Our advice would be however to avoid such LCGs or combined LCGs for modern simulations because of their structural weaknesses L'Ecuyer et al. (2002).

Leap frog

The leap frog (LF) is the way to partition a random stream in the manner of dealing cards to several players. Random numbers are allocated in turn to PEs, in the same way as cards are 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 [START_REF] Aluru | Lagged fibonacci random number generators for distributed memory parallel computers[END_REF].

Given the period P of the global sequence, the period of each stream is P/N . As with the splitting technique, the long-range correlations in the initial RNG can lead to small-range correlations between the potential substreams, particularly if we have a large number of PEs. In addition, Wu and Huang [START_REF] Wu | Parallel use of multiplicative congruential random number generators[END_REF] showed that depending on the interval used (i.e., the number of PEs and the length of the random sequences), poor spectral values could be observed. A case where the quality of the original RNG is seriously affected by the LF technique is shown in [START_REF] Hill | Good random number generators are (not so) easy to find[END_REF]. In addition, when this technique is used without jump ahead (Section 3.4), performances are divided by the number of PEs because of the bottleneck problem appearing in the central server technique.

This technique needs to be used with caution depending on the environment it is set up in. Thus, to preserve reproducibility between two executions of the same simulation, one should not use this technique when the parallelism grain is not fixed yet. As a matter of fact, different random streams will be assigned to PEs depending on the chosen grain. This situation is illustrated in Figure 3. For the classical LF approach, one must ensure that the original status can be shared between all the PEs when the underlying PRNG algorithm relies on a linear recurrence to produce the next numbers of the sequence. This limits the target architectures of the classic LF technique to shared memory architectures, where a PRNG state can be stored and accessed efficiently by the PEs. For distributed memory architectures, a smart implementation of the LF technique can be used when the number of PEs is a power of 2. In this case, using a unique Lagged Fibonacci generator, each PE can implement its own version of the generator that performs the right jump corresponding to the PE id [START_REF] Janowczyk | Fast, processor-cardinality agnostic PRNG with a tracking application[END_REF]. This approach is also discussed in Section 4.2, focusing on GP-GPU.

Partitioning multiple streams: 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 | High-performance monte carlo tools[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. An example of parameterization can be found in Figure 4, which states three independent parameterized PRNGs. In the case of LCGs and multiplicative congruential generators, this can rapidly lead to poor results De Matteis and Pagnutti (1990a); [START_REF] Wu | Parallel use of multiplicative congruential random number generators[END_REF] even when the parameters are very carefully checked. For instance Mascagni and Chi proposed that the modulus be Mersenne or Sophie Germain prime numbers [START_REF] Mascagni | Some methods of parallel pseudorandom number generation[END_REF] to improve the case of LCGs parallelization.

Other techniques

Now that the main partitioning methods have been presented, we can indulge in a bit of history of other approaches because the parallelization of pseudo-random numbers has been under study for at least 30 years.

Variants have been developed on the basis of classical methods such as the shuffling LF. Its principle is to parameterize both the seed and recursive functions of LCGs. One of the main contributions to this variant is that it results in a scalable period, that is, the number of different random numbers that can be used increases with the number of parallel streams. A parameterization method was used in [START_REF] Percus | Random number generators for MIMD parallel processors[END_REF] to obtain parallel streams from a LCG, by choosing for the j th stream a multiplier a(j) and an additive constant c(j). It is shown in this case that good results can be obtained if c(j) is the j th prime number less than (m/2) where m is the modulus. A specific way to generalize the partitioning methods was proposed in the 1980s [START_REF] Frederickson | Pseudo-random trees in monte carlo[END_REF]. It results in what has been called a 'Lehmer tree'. But the suggestion is apparently limited to LCGs, which we still discourage to use.

We can also think of hybrid approaches, but they are not distribution techniques in their own right. Yet they combine several standard distribution techniques to take advantage of their respective features. Please note that the resulting random stream of this combined approach needs to pass through the same test batteries than its parents to be validated. Indeed, combining several distribution techniques might not preserve the statistical quality of the original random stream.

In the end, let us add that Boolean Cellular Automata have been considered to generate parallel pseudo-random numbers by [START_REF] Sipper | Generating parallel random number generators by cellular programming[END_REF]; Tomassini Hill et. al (1999); [START_REF] Tomassini | Generating high-quality random numbers in parallel by cellular automata[END_REF], who tested the generated sequences for distributed environments. This technique is rather slow, and [START_REF] Ackermann | Parallel random number generator for inexpensive configurable hardware cells[END_REF] has considered its hardware implementation in programmable chips. We are not aware of any evaluation of this technique by any state-of-the-art testing battery; thus, we cannot advise the use of this technique.

A tool for partitioning: jump-ahead algorithms

A jump-ahead algorithm is a tool used by distribution techniques such as sequence splitting and LF to deal numbers efficiently. This technique enables the analytical computing of the generator state in advance after a huge number of cycles (generations) and corresponds to a jump ahead in the random stream. Knowing the recurrence formula s n+1 = f (s n ), where s n+1 is the generator status, the difficulty is to determine

f (n) = f •f •...•f (n-1 compositions of f by itself) because s n = f (n) (s 0 ).
Such a technique is interesting with generators that have very large periods. Among widespread generators, MT and WELL generators, on the basis of linear recurrences modulo 2 (F 2 -linear generators), embed an efficient software providing jump-ahead facilities. Matsumoto and L'Ecuyer teams joined [START_REF] Haramoto | Efficient jump ahead for f2-linear random number generators[END_REF] to develop a viable jump-ahead algorithm running in few milliseconds on current processors. In this case, f is a linear application in F w 2 , with a matrix A with w lines and w columns, where w represents the size of the state vector. The problem is to compute A n , and this can be achieved, thanks to the characteristic polynomial of A matrix. When an F 2 linear generator proposes a full period (2 w -1) and when the size of the w state vector is big, the computing of the characteristic polynomial will take time to compute L'Ecuyer et al. (2002).

More efficient algorithms exist, such as MRG32k3a from L'Ecuyer and Buist (2005). However, this generator is twice to three times as slow as MT on common 32-bit computers and also has a much shorter period length, although it can be sufficient for modern applications.

Joint use of the partitioning of a single stream and parameterization

The two main techniques of random streams distribution, presented in Sections 3.1 and 3.2, are not exclusive.

To simplify, we can say that any number xn of a sequence provided by a conventional PRNG is the result of a call x n = g(s n ), where s n is the internal state of the PRNG and g a function usually very simple. The statistical quality of these generators is therefore only in the complexity of the state transition function f that allows us to deduce the following state: s n+1 = f (s n ). However, the downside of this complexity is that it induces a sequential dependence between the successive states of the PRNG. The direct and rapid access to a state s n (required by the sequence splitting or LF techniques, for example) is possible only if one has a jump-ahead function (Section 3.4).

With counter-based PRNGs [START_REF] Salmon | Parallel random numbers: as easy as 1, 2, 3[END_REF], so named because f (s n ) = (s n + 1) mod 2 p (where p is the size of the state vector), this sequential dependence disappears because the numbers are generated by calls such as x n = g(n) where n is the state reduced to a simple counter. The parallelism inherent in the latter relationship allows immediate use of single stream partitioning techniques (Section 3.1).

The statistical quality of counter-based PRNGs are therefore due to the complexity of the g functions used (for PRNGs used in simulation, these functions are simplifications of cryptographic block ciphers, such as AES or Threefish, used in cryptographically secure PRNGs). These functions are indexed by keys, allowing a natural distribution of pseudo-random numbers by parameterization over the key space (Section 3.2).

Finally, the interest of the counter-based PRNGs is that the generation of numbers x n = g k (n) can be parallelized either by partitioning each stream according to the values of the counter n (Section 3.1) or by following a parameterization approach (Section 3.2) on the basis of the k keys. When the size of counter space and key space is sufficient (up to 2 128 and 2 64 , respectively, for counter-based PRNGs presented in [START_REF] Salmon | Parallel random numbers: as easy as 1, 2, 3[END_REF]), being able to use both techniques presented previously independently even allows, within each PE, to associate a PRNG to each software entity. For instance, in the case of an individualbased stochastic simulation with a very large number of individuals, the sequence splitting technique can be applied in the counter space to distribute the pseudo-random sequences on different replications of the simulation while setting on the keys used to assign each individual its own PRNG.

Discussions and taxonomy of random streams distribution techniques

For techniques such as sequence splitting and random spacing, we have seen that a common problem is overlapping, but we also have to consider the potential impact of the random initialization on the quality of the underlying PRNG. Recent history has shown that even some of the best RNG algorithms could fail when badly initialized. In the first version of the Mersenne Twister generator, if two initial states were too close with respect to the Hamming distance, then the corresponding output sequences were close to each other. Improvements have been proposed to overcome this problem1 . For the initialization problem, the remaining technique is to run empirical or statistical tests such as TestU01. In 2008, Reuillon proposed 1 million statuses for the first Mersenne Twister of period 2 19937 -1: he used a RNG with cryptographic qualities to propose independent and well-balanced bit statuses, knowing that when MT had a zero-excess initialization status, it could take a quite long number of draws to recover good statistical properties [START_REF] Reuillon | Testing 65536 parallel pseudo-random number streams[END_REF].

In Figure 5, we propose a Unified Modeling Language (UML) class diagram of the main parallelization techniques. The latter are basically ordered depending on the use of either a single stream or multiple streams. Then, the taxonomy is refined considering the strategies set up by distribution techniques: for instance, we distinguish techniques issuing blocks of contiguous numbers in opposition to those dealing numbers. Now, we have seen that techniques based upon unique original random streams might take advantage of a jump-ahead algorithm to improve their generation speed. Basically, this feature concerns the LF and sequence splitting approaches. However, we do not consider jump ahead as a relevant criterion to sort distribution techniques. That is why we decided to allow any UniqueOriginalStream instance to make use of jump ahead. Finally, please note that we chose not to distinguish a particular hybrid approach. Indeed, we consider as hybrid any technique combining at least two of its counterparts. This has been carried out, thanks to the composite design pattern on the right side of Figure 5, which allows us to describe the Hybrid Approach as a class that represents the combination of at least two other distribution techniques.

In a previous paper [START_REF] Reuillon | Declarative task delegation in OpenMOLE[END_REF], we have presented the DistMe software toolkit designed to help with the distribution of large parallel stochastic applications. It is a concrete implementation of the techniques and advice provided along this paper and is a concrete implementation of our experience in dealing with pseudo-random numbers for distributed simulations. This paper also presents examples of stochastic simulations for life science research where thousands of billions of pseudo-random numbers have been used. The evolution of this software Figure 5: Taxonomy of distribution techniques toolkit (OpenMOLE) has been presented in Reuillon (2008a). The preliminary design of the DistMe toolkit was achieved when tackling the distribution of a nuclear medicine application using the largest European computing grid (EGI) El Bitar et al. ( 2006); [START_REF] Li | Improving performance via computational replication on a large-scale computational grid[END_REF]; [START_REF] Maigne | Parallelization of monte carlo simulations and submission to a grid environment[END_REF]. At that time, we used sequence splitting and the dynamic creation of Mersenne Twister algorithms [START_REF] Reuillon | Parallel stochastic simulations with rigorous distribution of pseudo-random numbers with DistMe: application to life science simulations[END_REF]. Thanks to the EGI grid, the equivalent of a few years of computing was achieved in a few days. An interesting example of LF usage is given in [START_REF] Janowczyk | Fast, processor-cardinality agnostic PRNG with a tracking application[END_REF], and efficient usage of the jump-ahead technique is given by simulation software using Scalable Library for Pseudo-Random Number Generation (SPRNG) [START_REF] Mascagni | Parameterizing parallel multiplicative lagged-fibonacci generators[END_REF], SSJ L'Ecuyer ( 2001) and the variants of Rstream L' Ecuyer and Buist (2005).

Choosing the right distribution technique

To sum up the suggestions stated in this section, Table 1 proposes another taxonomy of distribution techniques. Here, we have focused on the environment, where the simulation is supposed to run, and the underlying PRNG, which original stream has to be split. This table should help practitioners to figure out which distribution technique they can use, provided their own environment constraints.

Random Numbers Parallelization Software

Techniques designed for CPUs

In the early 1990s, at the European Simulation Symposium, Pawlikowski and his colleagues proposed an interesting methodology for the parallelization and for the automatic partitioning of stochastic parallel simulations [START_REF] Hill | Towards credible and fast quantitative stochastic simulation[END_REF]; [START_REF] Pawlikowski | On credibility of simulation studies of telecommunication networks[END_REF]. As in [START_REF] Hill | URNG: a portable optimization technique for software applications requiring pseudo-random numbers[END_REF] and [START_REF] Pawlikowski | On automatic partitioning, run-time control and output analysis methodology for massively parallel simulations[END_REF], Li and Mascagni gave advice to achieve safe multiple replications in parallel L'Ecuyer ( 2010), but to our knowledge, the first reliable and sound library that takes into account the parallelization of pseudo-random numbers is SPRNG presented in [START_REF] Mascagni | Random number generation[END_REF]. We can find in the latter paper an overview of the mathematical grounds for random number generators, and there are also implementations of various parameterization techniques for different families of generators, which were presented in [START_REF] Mascagni | Parameterizing parallel multiplicative lagged-fibonacci generators[END_REF]. This library has been proved reliable for parallel Monte Carlo computations and also proposes a small test suite. An example of multiple replications in parallel is also given in [START_REF] Ewing | Akaroa-2: Exploiting network computing by distributing stochastic simulation[END_REF] dealing with performance evaluation studies of modern multimedia telecommunication networks. In [START_REF] Mascagni | SPRNG: a scalable library for pseudorandom number generation[END_REF], we find a short and interesting discussion around Monte Carlo tools for HPC at the end of the last millennium. As stated previously in the parameterization section, Matsumoto and Nishimura proposed the Dynamic Creator software to generate mutually independent Mersenne twister generators for parallel computing [START_REF] Matsumoto | Dynamic creation of pseudorandom number generators[END_REF]. This kind of approach is still a very good parallelization technique for MT generators even if we consider the recent development of efficient jump-ahead techniques for linear recurrences modulo 2 (MT and WELL generators) [START_REF] Haramoto | Efficient jump ahead for f2-linear random number generators[END_REF]. At the beginning of the millennium, L'Ecuyer and his team started to propose a package able to produce many long streams and substreams in C, C++, Java, and FORTRAN L' Ecuyer and Buist (2005); a version for R was proposed in 2005 L'Ecuyer (1996). In 2004, Coddington and Newell released a Java parallel random number library for HPC for the Java language [START_REF] Coddington | Random number generator for parallel computers[END_REF]. This library proposes three good generators (two by L'Ecuyer), parallelized by the sequence splitting or the indexed sequence technique. To our knowledge, the library is limited to shared memory machines. In our opinion, the current best library for Java is SSJ L'Ecuyer (2001). The SSJ package provides a RandomStream Java interface that defines the basic structures to handle multiple streams of uniform random numbers. Each stream of random numbers is a Java object for which the original sequence from a given period can be cut into adjacent streams (or segments) as in the sequence splitting approach. Efficient methods are proposed to move around streams.

The case of GP-GPU

Dealing with GP-GPU, we have only found a limited number of techniques and few high quality generators [START_REF] Bradley | Parallelization techniques for random numbers generators[END_REF]. The purpose here is not to list current PRNG implementations for GPUs. We have presented in Passerat-Palmbach et al. ( 2010) a survey of this kind, and in most cases, the purpose was to improve the generation speed, thanks to a GPU accelerator. If this point is interesting, it is not what simulationists are looking for. Our main concern is to be able to massively run independent stochastic functions on GP-GPUs (named 'kernels' in the NVIDIA CUDA terminology).

At the time of writing, we find interesting propositions through libraries, mainly CURAND and Thrust::random. CURAND has been introduced in the latest version of the CUDA framework and is 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 pseudo-random sequences either on GPU or on CPU. In this library, we find a Sobol quasirandom number generator, the pseudo-random number generator implemented being XorShift proposed by [START_REF] Marsaglia | A new class of random number generators[END_REF]. This kind of generator is fast, but it also presents statistical flaws as explained by Panneton in his PhD thesis [START_REF] Panneton | Improved long-period generators based on linear reccurences modulo 2[END_REF]. The interest of CURAND is that the API of the library stays the same no matter which kind of random sequence you select and the platform you run your application on. The second library found is named Thrust::random. It is part of a GPU-enabled general-purpose library called Thrust. This open-source project intends to provide a GPU-enabled library equivalent to classic general purpose C++ libraries such as STL or Boost. Classes are split through several namespaces, 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 an LCG, a linear feedback shift register and a subtract with borrow [START_REF] Marsaglia | DIEHARD: a battery of tests of randomness[END_REF][START_REF] Marsaglia | Xorshift RNGs[END_REF]. By the way, although the latter PRNG is mentioned as subtract with carry in Thrust::random documentation, the original Marsaglia's proposition is known as subtract with borrow. It is surprising that modern libraries propose such old generators, even if the library offers simple ways to combine them into better quality randomness sources, such as L'Ecuyer's Tausworthe combined generators L' Ecuyer and Simard (2007). In [START_REF] Janowczyk | Fast, processor-cardinality agnostic PRNG with a tracking application[END_REF], we can see a smart usage of the LF technique adapted to GP-GPU for particle filtering. This paper presents a Lagged Fibonacci generator adapted to the LF approach when the number of PEs is a power of 2. An efficient GP-GPU implementation is shown and used in many computer vision algorithms.

For modern GP-GPUs (starting from the introduction of NVIDIA Tesla boards), our opinion is to consider only recent generators that have already shown statistical strength according to modern and thorough testing. To propose multiple stochastic streams through a GP-GPU architecture, we suggest the use of the DC dedicated to the new MTGP generator proposed by Saito (as a GP-GPU implementation of Mersenne Twister [START_REF] Saito | A variant of mersenne twister suitable for graphics processors[END_REF]). This careful implementation benefits from GP-GPU shared and constant memory to improve generation speed. This is not the only aspect of the generator that retained our attention: its availability for both CPUs and GPUs architectures is very important for us. Indeed, MTGP has been designed to run on GPUs, but we can also find a CPU version on Matsumoto's Home Page (Saito MTGP part). We were able to test the PRNG on CPU-based hosts with reliable and well-known tools, and this will be discussed in the next section. This PRNG can be well suited for stochastic simulations following the hybrid computing paradigm: multi-scale and holistic simulation running parts of an application on a CPU host (possible 'manycores'), while the massive parallel side is executed on GP-GPU boards. With MTGP dynamic creation, we can rely on the same parallelization technique for both the CPU and the GPU. We have started to test this GPU version of Mersenne Twister, and it will be discussed in the next section. Still, this PRNG relies on a parameterization technique to provide independent random streams. When used on GPU, this technique needs to be coupled with another one because we cannot afford to store a Parameterized Status per thread. This drawback is the main reason why a 'tiny' version of MT has recently been released. It is promising for GP-GPU because it has a small status (16 bytes). Tested successfully with the BigCrush test battery (from TestU01) described in the next section, this generator comes with its DC software and also with an efficient jump-ahead implementation. This efficient implementation is achieved, thanks to its shorter period compared with the usual huge periods of MT family generators; this period is however large enough for random spacing (2 127 -1).

Besides Mersenne Twister-like and WELL algorithms, MRG32k3a L'Ecuyer and Buist (2005) from L'Ecuyer has also retained our attention for a GP-GPU use. In fact, it is designed to support natively jump ahead, making MRG32k3a a perfect match to partition a random stream following the sequence splitting strategy. In addition, these PRNG internal data structures display a small memory footprint (24 bytes) making it a definitely good candidate for GP-GPU-Hill et. al enabled stochastic applications. Please note that a CUDA version of this algorithm is described in [START_REF] Bradley | Parallelization techniques for random numbers generators[END_REF], although a valid implementation is available on GitHub2 . The counterbased random number generators are also recently available and very efficient for GP-GPUs [START_REF] Salmon | Parallel random numbers: as easy as 1, 2, 3[END_REF].

Statistical and Empirical Testing Software for Random Streams

Knuth proposed a set of statistical tests for random streams [START_REF] Knuth | The Art of Computer Programming[END_REF]. Marsaglia designed a testing suite, named DieHard, highly regarded for many years [START_REF] Marsaglia | A random number generator for PC's[END_REF]. The statistical test suite developed by the National Institute for Standards and Technology is also interesting, particularly when cryptographic qualities are required. As mentioned in the previous section, SPRNG is also providing a set of statistical tests. A detailed description of the main statistical tests for pseudo-random numbers was given in Rütti's thesis [START_REF] Rütti | A generic random number generator test suite[END_REF], and he also proposed a testing suite with some colleagues [START_REF] Rütti | A random number generator test suite for the C++ standard[END_REF]. The DieHarder testing suite is proposed and updated by [START_REF] Brown | DieHarder: a random number test suite[END_REF]. Although DieHarder is also an interesting testing suite, we highly recommend the TestU01 software library, which currently offers the most complete collection of utilities for the empirical statistical testing of uniform random number generators L' Ecuyer and Leydold (2005). In addition to the classical statistical tests for RNGs and the other tests previously cited and proposed in the literature, TestU01 proposes new original tests as well as predefined tests suites (Crush and BigCrush with more than a hundred tests). Test of bit sequences are included, and TestU01 also considers the size of the sample according to the period length and to the kind of test considered. The TestU01 software also proposes interleaving of random streams; this is particularly interesting in the context of parallel simulation, to test the influence of parallel substreams inter-correlations. Empirical techniques to test parallel random number generators have been proposed by Coddington and Ko Coddington and Newell (2004). In 2008, we considered the testing with TestU01 (Crush) of 65536 independent sequences created by DC Matsumoto and Nishimura (2000). This was achieved using the former EGEE (known as EGI) computing grid [START_REF] Reuillon | Parallel stochastic simulations with rigorous distribution of pseudo-random numbers with DistMe: application to life science simulations[END_REF]. In addition to testing software designed for sequential generators, Coddington proposed a set of criterion for 'parallel generators [START_REF] Coddington | Techniques for empirical testing of parallel random number generators[END_REF]: the generator should be able to work on any number of processors, the sequential sequence in use for each processor should satisfy the current statistical tests, the parallel sequences should be reproducible and lastly, as stated by many authors, parallel sequences or streams should be uncorrelated. To this set of criteria, we add the fact that each PE should possess its own sequence or subsequence to ensure the reproducibility of execution on each PE independent of the hardware architecture (PEs can be: threads, processes or processors). For instance, if we detect a bug during the execution on a specific PE, we have to be able to reproduce it with the same sequence and on any architecture. The need for mutual independence between the parallel sequences or streams has to be checked carefully to avoid long-range correlations De [START_REF] De Matteis | Parallelization of random number generators and longrange correlations[END_REF][START_REF] De Matteis | Controlling correlations in parallel monte carlo[END_REF], 1990a). The problem of long-range correlations has been identified in various simulation applications for many years, and techniques have been proposed to avoid them as far as possible because to our [START_REF] De Matteis | Long-range correlations in linear and non-linear random number generators[END_REF].

We achieved an analysis of MTGP by using the dynamic creation of this PRNG and facing the resulting generators up to the current most stringent TestU01 test battery: BigCrush. Few weaknesses were identified during these experiments when compared with the original MT. The purpose of our test was to obtain a large set of fine parameterized statuses allowing the initialization of MTGP without introducing any potential bias with a bad status. Only statuses that passed BigCrush were kept in this study [START_REF] Passerat-Palmbach | Reliable initialization of GPU-enabled parallel stochastic simulations using mersenne twister for graphics processors[END_REF]. However, problematic statuses are mostly embarrassed by the same two tests. We observed that in some cases with relatively small periods (2 3217 ), 30% of the parameters generated by dynamic creation led to MTGPs with failure to test 35 and 100 of TestU01. Thanks to precisions given by Makoto Matsumoto and Mutsuo Saito, we now know that test 35 discards the most significant 25 bits from the 32-bit words and then uses the next 5-bits. Matsumoto and Saito were extremely reactive and have already corrected the GPU version of MTGPDC to take this into account. Nonetheless, our study has shown that MTGP was particularly safe with longer periods, according to TestU01 criteria. Yet, the longer the period, the more space it needs to store its internal state vector used by its algorithm. And this point is currently noticeable because even the best current GP-GPU architectures propose a small size of fast shared memory. We have selected about 6000 fine MTGP parameterized statuses, with the lowest available period 2 3217 to reduce their GPU shared memory footprint. They have been provided to Matsumoto and Saito, and they will be publicly available.

Conclusion

We have tackled the use of random number generation for HPC and presented the main partitioning methods for stochastic parallel simulations. The current 'best' generators according to the latest test libraries have been discussed. We have also considered various existing tools because the use of random numbers in parallel stochastic simulations is still a challenging technical problem. A 'good' generator is one where statistical defaults are well hidden, although they are not inexistent because even the 'best' known generators can fail for a particular application or when they are badly initialized. Parallel simulations must first consider the choice of a nice sequential generator, but as usual, no statistical or practical test is universal because we cannot prove the statistical soundness of a generator. Test batteries are the only empirical way to ensure that the generator in use fits for a wide set of applications. Such batteries can be used to independently test parallel streams. In addition, for very sensitive simulations, it can be a good practice to test different generators and partitioning techniques to increase confidence in simulation results. In the particular case of using parallel random number streams, additional care should be taken to ensure that we avoid long-range correlations, even though it often implies an additional computing cost. In our opinion, it would be interesting to include random generator partitioning and advanced random streams testing facilities in HPC middleware for the scientific end users. An attempt to do so was proposed in the PhD thesis of a student of ours with the DistMe, DistRNG and DistTest toolkits [START_REF] Reuillon | Rigorous distribution of stochastic simulations using the DistMeToolkit[END_REF]. Another practical perspective could be to formalize an assessment process on the basis of the following tuple: the nature of the application, the random number generator selected, the partitioning method and the results of current test methods when available. A test method can be viewed as a rough model for a particular set of applications. Systematic testing of parallel sequences can be achieved once and for all for a particular generator, thanks to the current level of computing power available. As said previously, Reuillon has recently used the EGI Grid to achieve a test of 65536 parallel random streams for Mersenne Twister 19937 [START_REF] Reuillon | Parallel stochastic simulations with rigorous distribution of pseudo-random numbers with DistMe: application to life science simulations[END_REF]. The initial statuses have been generated with the DC proposed by Matsumoto and Nishimura. In the context of GP-GPUs, we have also recently tested the dynamic creation of GPU dedicated Mersenne Twister generators (MTGP), and the results presented in [START_REF] Passerat-Palmbach | Reliable initialization of GPU-enabled parallel stochastic simulations using mersenne twister for graphics processors[END_REF] led to a collaboration with Matsumoto's team. The use of GP-GPUs for intensive independent stochastic hybrid computing can be limited by the current hardware constraints, such as the small size of fast shared memory area, as well as the predominant SIMD paradigm. We have recently noticed the promising 'TinyMT' that is carefully crafted for the limited GP-GPU fast shared memory size. Furthermore, the new hybrid architectures proposed are precisely working on the two aforementioned hardware limitations. In conjunction with the increase of the number of 'cores' in the near future (Kepler and Maxwell boards from NVIDIA) and the improvement of GP-GPU programming friendliness, this could change the way many scientists will consider the use of 'desktop' HPC.

Figure 2 :

 2 Figure 2: Random spacing applied to 3 processing elements (PEs)

Figure 3 :

 3 Figure 3: Leap frog in a unique original stream considering 2 processing elements (PEs), then 4 PEs

Figure 4 :

 4 Figure 4: Example of Parameterization for 3 processing elements

http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/emt19937ar.html

https://github.com/jopasserat/ShoveRand

Distribution

Technique Where to use it?

Where to avoid it?


Keywords: Stochastic Simulation, Multiple Replications In Parallel (MRIP), GP-GPU, CUDA, Warp-Level Parallelism (WLP)

Stochastic simulations need multiple replications in order to build confidence intervals for their results. Even if we do not need a large amount of replications, it is a good practice to speed-up the whole simulation time using the Multiple Replications In Parallel (MRIP) approach. This approach usually supposes to have access to a parallel computer such as a symmetric multiprocessing machine (with many cores), a computing cluster or a computing grid. In this paper, we propose Warp-Level Parallelism (WLP), a GP-GPU-enabled solution to compute MRIP on GP-GPUs (General-Purpose Graphics Processing Units). These devices display a great amount of parallel computational power at low cost, but are tuned to process efficiently the same operation on several data, through different threads. Indeed, this paradigm is called Single Instruction, Multiple Threads (SIMT). Our approach proposes to rely on small threads groups, called warps, to perform independent computations such as replications. We have benchmarked WLP with three different models: it allows MRIP to be computed up to six times faster than with the SIMT computing paradigm.

INTRODUCTION

Replications are a widespread method to obtain confidence intervals for stochastic simulation results. It consists in running the same stochastic simulation with different random sources and averaging the results. According to the Central Limit Theorem, the average result is approximated in an accurate enough way by a Gaussian Law, for a number of replications greater than 30. Thus, for a number of replications greater than 30, we can obtain a confidence interval with a satisfactory precision.

There are many cases where a single simulation can last for a while, so 30 of them run sequentially may represent a very long computation time. Because of this overhead, 30 replications are hardly run in most simulations. Instead, a good practice is often to run 3 replications when debugging, and 10 replications are commonly used to compute a confidence interval. To maintain an acceptable computation time while running 30 or more replications, many scientists proposed to run in parallel these independent simulations. This approach has been named Multiple Replication in Parallel (MRIP) in the nineties [START_REF] Pawlikowski | Distributed stochastic discrete-event simulation in parallel time streams[END_REF]. As its name suggests, its main idea is to run each replication in parallel [START_REF] Hill | Object-oriented pattern for distributed simulation of large scale ecosystems[END_REF]; [START_REF] Pawlikowski | Towards credible and fast quantitative stochastic simulation[END_REF]. In addition, when we explore an experimental plan we have to run different sets of replications, with different factor levels according to the experimental framework [START_REF] Hill | Object-oriented analysis and simulation[END_REF]; [START_REF] Amblard | Mda compliant design of simexplorer a software tool to handle simulation experimental frameworks[END_REF]. In this paper, we will not consider any constraints that need to be satisfied when implementing MRIP. One of the main barriers that often prevents simulationists to achieve a decent amount of replications is, on the one hand, the lack of knowledge in the parallelization techniques and on the other hand the parallel computing facilities available. Our work tackles this problem by introducing a way to harness the computational power of GP-GPUs (General-Purpose Graphics Processing Units -GPUs hereafter), which are rather cheap compared to regular parallel computers, to process MRIP quicker than on a scalar CPU (Central Processing Unit).

GPUs deliver such an overwhelming power at a low cost that they now play an important role in the High Performance Computing world. However, this kind of devices display major constraints, tied to its intrinsic architecture. Basically, GPUs have been designed to deal with computation intensive applications such as image processing. One of their well-known limits is memory access. Indeed, since GPUs are designed to be efficient at computation, they badly cope with applications frequently accessing memory. Except by choosing the right applications, the only thing we can do to overcome this drawback is to wait for the hardware to evolve in such a way. Last NVIDIA GPUs generations, codenamed Fermi, show a move in this way by improving cache memories available on the GPU. This leads to better performances for most applications at no development cost, only by replacing the old hardware by the state-of-the-art one.

Now, what we can actually think about is the way we program GPUs. Whatever the programming language or architecture one chooses to develop his application with, CUDA (Compute Unified Device Architecture) or OpenCL, the underlying paradigm is the same: SIMT (Single Instruction, Multiple Threads). Thus, applications are tuned to exploit the hardware configuration, which is a particular kind of SIMD architecture (Single Instruction, Multiple Data). To obtain speed-ups, we must propose parallel applications that will be SIMD compliant. This point reduces the scope of GPU-enabled applications.

In the SIMT paradigm, threads are automatically grouped into 32-wide bundles called warps. Warps are the base unit used to schedule both computation on Arithmetic and Logic Units (ALUs) and memory accesses. Threads within the same warp follow the SIMD pattern, i.e. they are supposed to execute the same operation at a given clock cycle. If they do not, a different execution branch is created and executed sequentially every time a thread needs to compute differently from its neighbours. The latter phenomenon is called branch divergence, and leads to significant performance drops. However, threads contained in different warps do not suffer the same constraint. They are executed independently, since they belong to different warps.

In this paper, we introduce Warp-Level Parallelism (WLP), a paradigm to evaluate the approach of using GPUs to compute MRIP, using an independent warp for each replication. Our study will:

• Describe a mechanism to run MRIP on GPU;

• Propose an implementation of our approach: WLP;

• Benchmark WLP with three different simulation models.

GENERAL CONCEPTS OF GPU PROGRAMMING AND ARCHITECTURE

This section does a brief recall of the major concepts introduced by GPU programming and especially by CUDA. It also basically describes how a GPU architecture is organized, since these aspects are directly tied to our approach.

The Single Instruction Multiple Threads (SIMT) paradigm

SIMT is the underlying paradigm of any CUDA application. It is based on the well-known SIMD paradigm. While using SIMD, the same instruction is executed in parallel on multiple computational units, but take different data flows in input. Instead of viewing SIMT as a simple SIMD variant, one needs to understand that it has been created to simplify applications development on GPU. The main idea is first to allow developers to deal with a unique function, named a kernel, which is going to be run in parallel on the GPU. Second, developers manipulate threads in SIMT, which are a much more common tool nowadays than traditional vectors enabling SIMD parallelization.

In order to handle SIMT more easily, CUDA introduces different bundles of threads. As a matter of fact, threads are grouped into blocks, which size and 3D-geometry are defined by the user. The whole blocks of a kernel form a 2D grid. Each thread will be uniquely identified in the kernel thanks to an identifier computed from a combination of its own coordinates and of its belonging block's. More precisely, in addition to grid and blocks, CUDA devices automatically split threads into fixed-size bundles called warps. Currently, warps contain 32 threads. This group is extremely important in the low-level mechanisms running on a GPU.

As long as NVIDIA has defined both its GPU architecture and the SIMT paradigm, the latter is not only convenient, it also perfectly fits its host architecture. Its sole purpose is to be used on GPU architecture, which is quite different from other multi-core architectures, especially from CPU ones, as we will see in the next part.

Basic architecture of a GPU

While a CPU possesses few cores, each of them allowing the execution of one thread at a time, a GPU possesses a small number of Streaming Multiprocessors (SM) (for instance an NVIDIA Fermi C2050 has 14 SMs). Each SM embeds an important number of computational units (there are 32 floating point computational units -called Streaming Processor (SP) -on each SM of a Fermi C2050). In theory, the floating-point computation power of a GPU board is equal to the number of SMs multiplied by the number of SPs. Another figure that needs to be considered in the architecture is the number of warp schedulers. The latter are key elements of CUDA performance. In fact, memory accesses are done per warp. However, because of memory latency, the warps-schedulers select the warps that have their data ready to process. Consequently, the more warps can be scheduled, the more the memory latency can be hidden.

When the former generation of NVIDIA GPUs was issued with a single warp-scheduler per SM [START_REF] Kirk | Nvidia tesla : A unified graphics and computing architecture[END_REF], Fermi now owns two warp schedulers per SM [START_REF] Wittenbrink | Fermi gf100 gpu architecture[END_REF]. They are first employed when threads need to be scheduled on the SM they have been assigned to. In fact, threads within a warp also achieve memory accesses in parallel, before processing the same instruction on these data. To sum up, when threads are bound to each other, and must execute the same instructions according to SIMD machinery, warps are the smallest unit that run in parallel on the different SMs of a GPU, and are the smallest GPU element that is 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. Figure 1 shows a simplified representation of a SM of the Fermi architecture.

Blocks dispatching and warps scheduling in NVIDIA GPUs

Now that we have introduced the basic functioning of CUDA-enabled GPUs, let us detail the particular features that will help us to achieve MRIP on such architectures. We will see in this part how our GPU-enabled MRIP implementation relies on the scheduling features provided by NVIDIA CUDA devices. Previously, NVIDIA GPUs were only able to run a single kernel at a time. Thus, blocks of threads were dispatched through all the available SMs in a more or less logical way: SMs were activated in turn, striding indices four by four. When every SM had been activated, the process started again.

One of the key features of the cutting-edge Fermi architecture is the ability to run several kernels in parallel on the same device. To do so, the way blocks of threads are dispatched through the device has been redesigned in a new fashion. Now, every block of threads, no matter which kernel it belongs to, is first handled by a top-level scheduler referred to as the GigaThread Engine. It is supposed to dispatch blocks of threads to the Streaming Multiprocessors (SMs). The point is CUDA has always proposed asynchronous kernels calls to developers. Now that Fermi-enabled devices can run several kernels in parallel, GigaThread needs to take into account any potential upcoming kernel. Consequently, the dispatcher cannot reserve all the SMs to run a first kernel, given that a second one could be launched at any time. When the second kernel appears, some resources will still be available so that they can be assigned to the new kernel blocks.

Moreover, GigaThread enables immediate replacement of blocks on an SM when one completes executing. Since context switching has been fastened with Fermi, blocks of threads can fully take advantage of the hardware device thanks to GigaThread dispatching capabilities. From an external point of view, and since we do not have the real specifications, we have noted that the dispatching of blocks does not seem to be deterministic. NVIDIA uses a specific way to place blocks on SMs: indeed, SMs will not be enabled in order. SMs bearing non-consecutive identifiers will in fact run consecutively ordered blocks.

A WARP MECHANISM TO SPEED UP REPLICA-TIONS

Two problems arise when trying to port replications to GPU threads, considering a replication per thread. First, we generally compute few replications, whereas we have seen that GPUs needed to achieve large amounts of computations to hide their memory latency. Second, replications of stochastic simulations are not renowned for their SIMD-friendly behaviour. Usually, replications fed with different random sources will draw different random numbers at the same point of the execution. If a condition result is based on this draw, divergent execution paths are likely to appear, forcing threads within a same warp to be executed sequentially because of the intrinsic properties of the device.

The idea that we propose in this paper is to take advantage of the previously introduced warp mechanism to enable fast replications of a simulation. Instead of having to deal with Thread-Level Parallelism (TLP) and its constraints mentioned above, we place ourselves at a slightly higher scope to manipulate warps only. Let this paradigm be called Warp-Level Parallelism (WLP), as opposed to TLP. Now running only one replication per warp, it is possible to have each replication to execute different instructions without being faced to the branch divergence problem.

But to successfully enable easy development of simulation replications on GPU using one thread per warp, two mechanisms are needed.

First, it is necessary to restrict each warp to use only one valid thread. By doing so, we ensure not to have divergent paths within a warp. Moreover, we artificially increase the device's occupancy, and consequently, we take advantage of the quick context switching between warps to hide slow memory accesses. Theoretically, we should use the lowest block size maximizing occupancy. For instance, a C2050 board owns 14 SMs, and can schedule at most 8 blocks per SM. In this case, the optimal block size when running 50 replications would be 32 threads per block. This situation is represented in Figure 2, where we can see two warps running their respective first threads only. The 31 remaining threads are disabled, and will stall until the end of the kernel. Unfortunately, the GigaThread scheduler, introduced in the previous section, does not always enable a kernel to run on every available SM. In addition, SMs' memory constraints might compromise this ideal case by reducing the number of available blocks per SM. Second, there has to be an easy solution to get a unique index for each warp. TLP relies essentially on threads identifiers to retrieve or write data back. Thus, WLP needs to propose an equivalent mechanism so that warps can be distinguished to access and compute their own data.

Thanks to the two tools introduced in this section, it is possible to create a kernel where only one thread per warp will be valid, and where it will be easy to make each valid thread compute different instructions, or work on different data depending on the new index.

Although we could not figure out the real behavior of the GigaThread Engine dispatcher, the characteristics noticed in this part are sufficient to evaluate the performance of the dispatching policy. Furthermore, the new scheduling features introduced in Fermi significantly enhance the overall performance of our warp-based approach, given that it highly relies on warp scheduling and block dispatching.

IMPLEMENTATION

Now that we have defined our solution, we will propose an implementation in this section. To do so, we need to focus on two major constraints: first, we should keep a syntax close to C++ and CUDA, so that users are not confused when they use our approach. Second, we need to propose compile-time mechanisms as much as possible. Indeed, since WLP only exploits a restricted amount of the device's processing units, we have tried to avoid any overhead implied by our paradigm.

This paper intends to prove that our approach is up and running. Thus, this section will only introduce a restricted number of keywords used by WLP. As we have seen previously, we first have to be able to identify the different warps, in the same way SIMT does with threads. One way to obtain the warp identifier is to compute it at runtime. Indeed, we know that warps are formed by 32 threads in current architectures [NVIDIA2011a]. Thus, knowing the running kernel configuration thanks to CUDA defined data-structures, we are able to figure out the warp identifier with simple operations only, similarly to what have done [START_REF] Hong | Accelerating cuda graph algorithms at maximum warp[END_REF]. The definition of a warpIdx variable containing the warp's identifier can be written as in Figure 3 Conceptually, this definition is ideal because warpIdx is declared as a 'constant variable', and the warp identifier does not change during a kernel execution. This formula fits with the CUDA way to number threads, which first considers threads' x indices, then y and finally z, within a block. The same organization is applied to blocks numbering [START_REF] Kirk | Nvidia tesla : A unified graphics and computing architecture[END_REF]. Please note that the warpSize variable is provided by CUDA. This makes our implementation portable since warp sizes may evolve in future CUDA architectures.

Although this method introduces superfluous computations to figure out the kernel's configuration, we find it easier to understand for developers. Another way to compute the warp's identifier would have been to write CUDA PTX assemblyNVIDIA (2011). The latter is the Instruction Set Architecture (ISA) currently used by CUDA-enabled GPUs. CUDA enables developers to insert inlined PTX assembly into CUDA high-level code, as explained in [NVIDIA2011b]. However, this method is far less readable than ours, and would not be more efficient since we only compute warpIdx once: at initialization. This warp identifier will serve as a base in WLP. When classical CUDA parallelism makes a heavy use of the runtime-computed global thread identifier, WLP proposes warpIdx as an equivalent. Now that we are able to figure out threads' parent warps, let us restrain the execution of the kernel to a warp scope. Given that we need to determine whether or not the current thread is the first within its belonging warp, we will be faced to problems similar to those encountered when trying to determine the warp identifier. In fact, a straightforward solution reckoning on our knowledge of the architecture quickly appears. It consists in determining the global thread identifier within the block to ensure it is a multiple of the current warp size. Once again, the kernel configuration is issued by CUDA intrinsic data structures, but we still need a reliable way to get the warp size to take into account any potential evolution. Luckily, we can figure out this size at runtime thanks to the aforementioned warpSize variable. Consequently, here is how we begin a warp-scope kernel in WLP: We now own the bricks to perform WLP, but still lack a user-friendly API. Indeed, it would not be adapted to ask our users to directly use complex formulas without having wrapped them up before in higher-level calls. To do so, we chose to use macros, for the sake that they are compiletime mechanisms, thus not causing any runtime overhead, and that they are perfectly handled by nvcc, the CUDA compiler. Our previous investigations result in two distinct macros: WARP_BEGIN and WARP_INIT, which respectively mark the beginning of the warp-scope code portion, and correctly fill the warp identifier variable. When WARP_INIT presents no particularities, except the requirement to be called before any operations bringing into play warpIdx, WARP_BEGIN voluntary forgets the block-starting brace following the if statement. By doing so, we expect users to place both opening and closing braces of their WLP code if needed, just as they would do with any other block-initiating keyword.

To sum up, please note once again that this implementation mainly targets to validate our approach. Still, it lays the foundation of a more complete API dedicated to WLP. The efficient but not appealing intrinsic mechanisms are totally masked to users thanks to macros introduced in WLP.

RESULTS

In this part, we introduce three well-known stochastic simulation models in order to benchmark our solution. We have compared WLP's performances on a Tesla C2050 board to those of a state-of-the-art scalar CPU: an Intel Westemere running at 2.527 GHz. For all of the three following models, each replication runs in a different warp when considering the GPU, whereas the CPU runs the replications sequentially. The following implementations use L'Ecuyer's Tausworthe three-component PRNG, which is available on both CPU and GPU respectively through Boost.Random and Thrust.Random [START_REF] Hoberock | Thrust: A parallel template library[END_REF] libraries. Random streams issued from this PRNG are then split into several sub-sequences according to the Random Spacing distribution technique [START_REF] Hill | Practical distribution of random streams for stochastic high performance computing[END_REF].

Description of the models

First, we have a classical Monte Carlo simulation used to approximate the value of Pi. The application draws a succession of random points coordinates. The number of random points present in the quarter of a unit circle are counted and stored. At the end of the simulation, the Pi approximation corresponds to a ratio of the points in the quarter of a unit circle to the total number of drawn points. The output of the simulation is therefore an approximate of Pi value. This model takes two input parameters: the number of random points to draw and the number of replications to compute.

The second simulation is a M/M/1 queue. For each client, the time duration before its arrival and the service time is randomly drawn. All other statistics are computed from these values. The program outputs are the average idle time, the average time in queue of the clients and the average time spent by the clients in the system. Because it did not impact the performances, the parameters of the random distribution are static in our implementation. Only the number of clients in the system and the number of replications, which modify the execution time, can be specified when running the application.

The last simulation is an adaptation of the random walk tests for PNRGs exposed in Vattulainen and [START_REF] Vattulainen | Mission impossible: Find a random pseudorandom number generator[END_REF]. The idea is to simulate a walker moving randomly on a chessboard-like map. The original application tests the independence of multiple flows of the same PRNG. To achieve this, multiple random walkers are run with different initializations of a generator on identically configured maps. Basically, each walker computes a replication. In the end, we count the number of walkers in every area of the map. Depending on the PRNG quality, we should find an equivalent number of walkers in each area. When the original version splits the map in four quarters, our implementation uses 30 chunks to put the light on the opportunity of our approach when there are many divergent branches in an application.

Comparison CPU versus GPU warp

As we can see in Figure 5, the CPU computation time of the Monte Carlo application approximating the value of Pi grows linearly with the number of replications. The GPU computation time increases only by steps. This behaviour is due to the huge parallel capability of the device. Until the GPU card is fully used, adding another replication does not impact the computation time, because they are all done in parallel. So, when the board is full, any new iteration will increase the computation time. This only happens on the 65 th replication because the GPU saved some resources in case a new kernel would have to be computed simultaneously. The same mechanism explains that after this first overhead, a new threshold appears and so on.

Due to this behaviour, GPUs are less efficient than CPUs when the board is nearly empty. When less than 30 replications are used, more than two-thirds of the board computational power is idle. Because sequential computation on CPU is widely faster than sequential computation on GPU, if only a little of the parallel capability of the card is used, the GPU runs slower. But when the application uses more of the card parallel computation power, the GPU becomes more efficient than the CPU. The pattern is very similar for the second model: the M/M/1 queue (see Figure 6). When the board does not run enough warps in parallel, the CPU computation is faster than the GPU one. But with this model, the number of replications needed for the GPU approach to outperform the CPU is smaller than what we obtained with the previous simple model. The GPU computation is here faster as soon as 20 replications are performed, when it required 30 replications to show its efficiency with the first model. This can be explained by GPUs' architecture, where memory accesses are far more costly than floating point operations in terms of processing time. If the application has a better computational operations per memory accesses ratio, it will run more efficiently on GPU. Thus, the GPU approach will catch with the CPU one faster. This point is very important because it means that depending on the application characteristics, it can be adequate to use this approach from a certain number of replications, or not. A solution is to consider the warp approach only when the number of replications is big enough to guaranty that most of the applications will run faster.

Comparison GPU warp versus GPU thread

If the advantages of WLP-enabled replications compared to CPU ones in terms of computation time have been demonstrated with the previous examples, it is necessary to determine if WLP outperforms the classic TLP. This case study has been achieved using the last model introduced: our adaptation of the random walk. Figure 7 shows the computation time noticed for each approach: CPU, GPU with WLP and GPU with TLP (named thread in the caption). Obviously, CPU and WLP results confirm the previous pattern: the CPU computation time increases linearly when the WLP one increases by steps. TLP follows logically the same evolution shape as WLP. Although it is impossible to see it here because the number of replications is too small, it also evolves step by step, similarly to the warp approach. WLP consumes a whole warp for each replication. In the same time, TLP activates 32 threads per warp. Thus, the latter's steps will be 32 times as long as WLP's. Having said that, we easily conclude that the first step in TLP will occur after the 2048 th replication.

As we can see in Figure 7, the computation time needed by the thread approach is significantly more important than the computation time of the warp approach (about 6 times bigger for the first 64 replications). But WLP catches up with TLP when the number of replications increases. When more than 700 replications are performed, the benefit of using the warp approach is greatly reduced. The best use of the warp approach for this model is obtained when running between 20 and 700 replications. Please note that this perfectly matches our replications amount requirement. It even allows the user to run another set of replications according to an experimental plan, or to run another set of replications with a different high quality PRNG. The latter practice is a good way to ensure that the input pseudo-random streams do not bias the results.

These results are backed up by the output of the NVIDIA Compute Profiler for CUDA applications. The latter tool allows developers to visualize many data about their applications. In our case, we have studied the ratio between the time spent accessing global memory versus To explain this ratio, let us recall that computation time was lower for WLP. Since the same algorithm is computed by the two different approaches, we should have noticed the same amount of Global Memory accesses in the two cases. In the same way, the profiler indicates significant differences between Global Memory reads and writes for TLP and WLP. These figures are summed up in

CONCLUSION

This paper has shown that using GPUs to compute MRIP was both possible and relevant. Having depicted nowadays GPUs' architecture, we have detailed how warp scheduling was achieved on such devices, and especially how we could take advantage of this feature to process codes with a high rate of branch divergent parts. Our approach, WLP (Warp-Level Parallelism), intends to allow users to easily distribute their experimental plans with replications on GPU.

WLP has been implemented thanks to simple arithmetic operations. Consequently, WLP displays a minimalist impact on the overall runtime performance. For the sake of user-friendliness, the internal mechanisms enabling WLP have been wrapped in high-level macros. At the time of writing, our version is functional and allows users to create blocks of code that will be executed independently on the GPU. Each warp will run an independent replication of the same simulation, determined by the warp identifier figured out at runtime. By doing so, we prevent performances to drop as they would do in an SIMT environment confronted to branch-divergent execution paths. WLP also tackles the GPU underutilization problem by artificially increasing the occupancy.

To demonstrate our approach performances, we have compared the execution times of a sequence of independent replications for three different stochastic simulations. Results show that WLP is at least twice as fast as cutting-edge CPUs when asked to compute a reasonable amount of replications, that is to say more than 30 replications. This will always be the case when a stochastic simulation is studied with a design of experiments, where for each combination of deterministic factors we have to run at least 30 replications, according to the previously mentioned Central Limit Theorem. WLP also overcomes the traditional CUDA SIMT performances by up to 6 to compute the same set of replications. Here, SIMT suffers of an underutilized GPU, whereas WLP takes advantage of a quick warp scheduling.

Insofar performances of WLP increase with the recent Fermi architecture compared to Tesla, we can expect this approach to be even more efficient with future CUDA architectures. We will validate this approach with bigger simulation models. As a matter of fact, two parameters need to be considered to determine how WLP will scale. On the one hand, a bigger model will often be more complicated, and will consequently contain much more divergent branches. When our approach should benefit of this aspect, on the other hand, bigger models will also consume more memory, which is the bottleneck of GPU devices.

The current version of WLP forces users to distribute their replications with our keywords. The target audience of our approach should, for the moment, be familiar with CUDA or GPU development. To lower the level of technical difficulty, we are currently thinking about an automatic tool, taking a simulation model and the number of replications to process in input, and producing the WLP equivalent in output, thus fully automating MRIP on GPU.

Figure 1 :

 1 Figure 1: Simplified architecture of an NVIDIA Fermi C2050 streaming multiprocessor

Figure 2 :

 2 Figure 2: Representation of thread disabling to place the application at a warp-level

  Figure 3: Const-definition of warpIdx

Figure 4 :

 4 Figure 4: Directive enabling warp-scope execution

Figure 5 :

 5 Figure 5: Computation time versus number of replications for the Monte Carlo Pi approximation with 10000000 draws

Figure 6 :

 6 Figure 6: Computation time versus number of replications for a M/M/1 queue model with 10000 clients

Figure 7 :

 7 Figure 7: Computation time versus number of replications for a random walk model with 1000 steps (above: 100 replications, below: 1000 replications)

Figure 8 :

 8 Figure 8: Comparison of TLP and WLP ratio of the overall Global Memory access time versus computation time

Table 1 :

 1

		TLP WLP
	Reads	225	18
	Writes	302	104

Table 1 :

 1 Number of read and write accesses to Global Memory for TLP and WLP versions of the Random Walk

1 This author's PhD is partially funded by the Auvergne Regional Council and the FEDER