Jonathan Passerat-Palmbach

David R C Hill

How to Correctly Deal With Pseudorandom Numbers in Manycore Environments -Application to GPU programming with Shoverand

Keywords: pseudorandom numbers, manycore, GP-GPU, High Performance Computing, Shoverand, stochastic simulation

Stochastic simulations are often sensitive to the source of randomness that characterizes the statistical quality of their results. Consequently, we need highly reliable Random Number Generators (RNGs) to feed such applications. Recent developments try to shrink the computation time by relying more and more General Purpose Graphics Processing Units (GP-GPUs) to speed-up stochastic simulations. Such devices bring new parallelization possibilities, but they also introduce new programming difficulties. Since RNGs are at the base of any stochastic simulation, they also need to be ported to GP-GPU. There is still a lack of well-designed implementations of quality-proven RNGs on GP-GPU platforms. In this paper, we introduce ShoveRand, a framework defining common rules to generate random numbers uniformly on GP-GPU. Our framework is designed to cope with any GPU-enabled development platform and to expose a straightforward interface to users. We also provide an existing RNG implementation with this framework to demonstrate its efficiency in both development and ease of use.

Introduction

More than 50% of High Performance Computing (HPC) applications now require the use of Random Number Generators (RNGs). We have had RNGs adapted to intensive parallel computing at our disposal for more than a decade . But even though the subject is well mastered by specialists like Pierre L'Ecuyer or Makoto Matsumoto, this is not the case for many colleagues who are acknowledged as specialists in other domains but are not aware of the recent parallelization techniques than can be used with modern generators. Even if we now also have access to true random numbers with quantum devices, they are still clumsy to produce in parallel with satisfying rates and costs. In addition, such devices are still subject to partial breakdowns and even the latest ones do have some very little biases that have to be corrected (such processing is sometimes included inside the device). In addition, stochastic experiments have to be reproducible, not only for debugging purposes but also for variance reduction techniques, for sensitivity analysis and for many other statistical techniques [START_REF] Kleijnen | Statistical tools for simulation practitioners[END_REF]. This implies storing the sequences (in order to be able to reproduce exactly the same sequence several times). This last point leads to the fact that we cannot provide enough "true" numbers for many High Performance Computing applications such as parallel Monte Carlo for nuclear medicine [START_REF] Lazaro | Fully 3d monte carlo reconstruction in spect: a feasibility study[END_REF][START_REF] Bitar | Fully 3d monte carlo image reconstruction in spect using functionalregions[END_REF]. Consequently, software random number generation remains the prevailing method for HPC, and specialists have been warning us for years to be particularly careful when dealing with parallel stochastic simulations [START_REF] De Matteis | Controlling correlations in parallel monte carlo[END_REF][START_REF] Hellekalek | Don't trust parallel monte carlo![END_REF][START_REF] Pawlikowski | Towards credible and fast quantitative stochastic simulation[END_REF].

Recent developments try to shrink computation time by relying more and more on General Purpose Graphics Processing Units (GP-GPUs) to speed-up stochastic simulations. Such devices bring new parallelization possibilities, but they also introduce new programming difficulties. Since the introduction of Tesla boards, Nvidia, ATI and other manufacturers of GP-GPUs have changed the way we use our high computing performance resources. Since 2010, we have seen that the top supercomputers are now often hybrid. Since RNGs are at the base of any stochastic simulation, they also need to be ported to GP-GPU. There is still a lack of welldesigned implementations of quality-proven RNGs on GP-GPU platforms. Consequently, we also need a survey of the current Pseudo Random Numbers Generators (PRNG) available on GPU. We will discuss the recent Mersenne Twister for Graphics Processors (MTGP) that was initially released in 2010, but also more recent generators with cryptographic inspiration that have been presented lately at the annual supercomputing conference.

In this paper we intend to present the good practices when dealing with pseudorandom streams in parallel. After explaining the kind of problems we may encounter, the state of the art in terms of RNGs, testing suite and parallelization techniques is presented. Then, GPU considerations and hybrid computing constraints are exposed. Next, the Shoverand framework that facilitates a safe usage of pseudorandom streams for modern GPU hardware is presented, and we also show how to integrate your preferred PRNG in Shoverand if your RNG of choice is not already included in the framework.

Using Parallel Random Streams

Parallel stochastic simulation

Parallel and Distributed Simulation (PDS) is an extensive research domain where effective solutions have been developed. Deterministic communications protocols for synchronous simulation and asynchronous simulation have been studied to avoid deadlocks and to preserve the causality and determinism principles. When dealing with stochastic simulations, random numbers should be generated in parallel in order for each Processing Element (PE) to be able to autonomously get its own independent random number stream. If such autonomy is not guaranteed, the parallelism is affected [START_REF] Coddington | Random number generator for parallel computers[END_REF].

State of the art in terms of sequential generators

The most famous sequential generator currently available is Mersenne Twister (MT hereafter) [START_REF] Matsumoto | Mersenne twister: A 623-dimensionnally equidistributed uniform pseudorandom number generator[END_REF]. Since the initial proposition at the end of the nineties, a whole family of MT generators has been proposed. SFMT, a member of this family, is an SIMD-oriented version of the original MT generator [START_REF] Saito | Simd-oriented fast mersenne twister: a 128-bit pseudorandom number generator[END_REF]. SFMT proposed the following improvements: speed (twice as fast as MT), a better equidistribution, a quicker recovery from bad initialization (zero-excess in the initial state) and even an increased period length (ranging from 2 607 -1 to 2 216091 -1). A GP-GPU version of MT (named MTGP) was also proposed by Saito and Matsumoto; it comes with companion software for parallelization (MTGPDC) [START_REF] Saito | Variants of mersenne twister suitable for graphic processors[END_REF]. Based on similar principles (Generalized Feedback Shift Register), Panneton, in collaboration with L'Ecuyer and Matsumoto, proposed the WELL generators with even better statistical properties [START_REF] Panneton | Improved long-period generators based on linear reccurences modulo 2[END_REF].

Pierre L'Ecuyer suggests that multiple recurrence generators (MRGs) with much smaller periods (above 2 100 but under 2 200 ) like MRG32k3a [START_REF] Pierre | An object-oriented randomnumber package with many long streams and substreams[END_REF] can suffice for modern applications. This generator comes with very interesting statistical properties and is easy to parallelize. Smaller periods also authorize fast mathematical parallelization techniques. The MT family now proposes TinyMT which can benefit of this point. In 2011, Salmon et al recently introduced statistically sound counter based pseudorandom generators with relatively small periods [START_REF] Salmon | Parallel random numbers: as easy as 1, 2, 3[END_REF]. These proposals led to very interesting implementations on GP-GPUs.

With this large set of 'good' sequential RNGs, the main question is: how can we make a safe RNG repartition in order to keep, on the one hand, efficiency, and on the other hand a sound statistical quality of the simulation in order to obtain reliable results.

Partitionning techniques

Assigning random sequences to parallel processing elements (PEs) can be done with one of the two following approaches. The first one proposes to partition a unique random stream, whereas the second approach deals with multiple independent streams. In this case, independent streams are obtained by parameterizing a family of generators.

When dealing with the partitioning of a unique random source we find the following variants in literature. The Central Server approach runs a single PRNG and provides on demand pseudorandom numbers to different PEs. We can also cite Boolean cellular automata since they have been considered to generate parallel pseudorandom numbers, but we are not aware of any recent high performance stochastic simulations using this technique. Another technique is known as the Sequence Splitting method, sometimes named Blocking or Regular Spacing. The principle behind those names is simple: the main sequence is split into 'N' non-overlapping contiguous blocks. Until recently, the computing of a jump in the sequence was unavailable for modern generators with linear recurrences modulo 2 and huge periods (such as the families of MT & WELL generators). This problem was solved in [START_REF] Haramoto | Efficient jump ahead for f2-linear random number generators[END_REF], but despite its relative efficiency, it is still considered as slow by specialists (compared to what can be obtained for MRG32k3a with a smaller period). The Indexed Sequence, or Random Spacing Technique, is another technique that consists in initializing the same generator with different random statuses. The last known technique to partition a unique stream is known as the Leap Frog Technique. With this method, random numbers are distributed to processing elements like a deck of cards dealt to card players.

When a single stream is partitioned into many different substreams, we talk about Parameterization. This technique issues many instances from the same family of generators. Independent generators, with different parameters, are used for each processing element (PE). For the Mersenne Twister family this technique is named Dynamic Creator (DC) and generates by parameterization highly independent Mersenne Twister generators. Our best mathematicians state that PNRGs based on linear recurrences defined by matrices with characteristic polynomials relatively prime to each other are supposed to be highly independent. In addition to the original MT generator, TinyMT and MTGP for GP-GPUs benefit from the same approach. Even though it is considered safe according to the scientific community, we recommend testing the resulting generators. In case of failure we can give feedback to authors and this was recently done for MTGP [START_REF] Passerat-Palmbach | Reliable initialization of gpu-enabled parallel stochastic simulations using mersenne twister for graphics processors[END_REF].

The interested reader can find the details of the different approaches in the following surveys [START_REF] David | Practical distribution of random streams for stochastic high performance computing[END_REF][START_REF] David | Distribution of random streams for simulation practitioners[END_REF], the latter being the most complete reference (a full journal paper with the latest updates), with discussions around advantages and drawbacks of each technique.

Testing sequential and paralell RNGs

The previous 'fine' generators were found reliable thanks to thorough statistical and empirical testing. Testing techniques were initially proposed by Knuth and the first testing software that became famous in the nineties was the DieHard testing suite proposed by Marsaglia. Another famous Statistical Testing Suite is proposed by the National Institute for Standards and Technology (NIST STS) and is particularly appreciated when cryptographic qualities are required. [START_REF] Brown | Dieharder: A random number test suite[END_REF] developed the DieHarder suite as an update of Marsaglia's work. Rütti and Troyer also presented a testing suite with Petersen [START_REF] Rütti | A generic random number generator test suite[END_REF]. But currently, the most complete collection of utilities for the empirical and statistical testing of uniform random number generators is in our opinion TestU01. This software library proposed by [START_REF] Pierre | Testu01: A c library for empirical testing of random number generators[END_REF] is now the reference for most scientists of the domain.

For parallel testing, Srinivasan and Mascagni gave interesting advice [START_REF] Srinivasan | Testing parallel random number generators[END_REF] based on the tests proposed by Mascagni in the Scalable library for Pseudorandom Number Generation (SPRNG) in 1997. Approximately at the same time, [START_REF] Coddington | Random number generator for parallel computers[END_REF] proposed a set of techniques for empirical testing of parallel random number generators. We have to be aware that we do not have at our disposal mathematical theorems or techniques to explicitly give a proof of independence between generators or between generated parallel random streams. To avoid long-range correlations we can have a look at interesting propositions made by [START_REF] De Matteis | Parallelization of random number gemerators and long-range correlations[END_REF][START_REF] De Matteis | Controlling correlations in parallel monte carlo[END_REF].

Software aid for the distribution of random streams

The reference software library dealing with parallelization of pseudorandom numbers is SPRNG (Scalable library for Pseudorandom Number Generation) [START_REF] Mascagni | Some methods of parallel pseudorandom number generation[END_REF]. As stated before, SPRNG also proposes a small test suite. The Dynamic Creator discussed above and proposed initially by Matsumoto and Nishimura is a software aid to provide mutually independent Mersenne twister generators for parallel computing, which was published two years later [START_REF] Matsumoto | Dynamic creation of pseudorandom number generators[END_REF]. L'Ecuyer and his team proposed a package able to produce many long streams and substreams in C, C++, Java, and FORTRAN [START_REF] Pierre | An object-oriented randomnumber package with many long streams and substreams[END_REF] and also in R in a later version. Coddington and Newell proposed the JaPara library (a Java Parallel Random Number Library) for High-Performance Computing [START_REF] Coddington | Japara-a java parallel random number generator library for high-performance computing[END_REF]. If we consider Java, our advice would be to use the SSJ package, still from L'Ecuyer's team. It provides the basic structures for handling multiple streams, with efficient methods to move streams around [START_REF] Ecuyer | Simulation in java with ssj[END_REF].

We also proposed higher level software tools to help in the distribution of stochastic simulations on local clusters and institutional computing Grids. The DistMe toolkit [START_REF] Reuillon | Rigorous distribution of stochastic simulations using the distmetoolkit[END_REF][START_REF] Reuillon | Parallel stochastic simulations with rigorous distribution of pseudo-random numbers with distme: Application to life science simulations[END_REF] and the OpenMole software based on declarative task delegation [START_REF] Reuillon | Declarative task delegation in openmole[END_REF] helps removing the burden of random streams distribution.

In the first part of this paper, we have given a shallow survey of distribution techniques, libraries and tools. The second part focuses on Shoverand: our library proposal integrating the guidelines introduced in the previous lines.

Purpose of Shoverand

As we have seen in the previous parts of this article, it is very important to deal carefully with pseudorandom numbers distribution when working in parallel environments such as GPUs. Still we cannot ask any user that wants to leverage GPUs' computational power in his stochastic simulations to be aware of every theoretical consideration that will prevent his simulation results from being biased. Thus, we introduced Shoverand [START_REF] Passerat-Palmbach | Shoverand: a model-driven framework to easily generate random numbers on gp-gpu[END_REF], a framework that provides Pseudorandom Number Generation (PRNG) facilities to CUDA-enabled GPU applications. Shoverand combines several aspects to ease developments of stochastic-enabled applications on GPU. First, its API is quite similar to what can be encountered when using high-level CPU languages like C++ or Java. Second, Shoverand's main goal is to handle the distribution of stochastic streams automatically without any intervention from the user. Finally, our framework also targets PRNG developers: indeed, Shoverand only integrates third-party PRNGs and focuses on unifying their interface. To do so, we integrate compile-time constraints that check whether the algorithm meets our guidelines.

PRNGs embedded in Shoverand

At the time of writing, Shoverand embeds several PRNGs. All these algorithms have been selected according to their intrinsic properties. We first consider their statistical properties in a sequential environment, because a PRNG could not cope with the requirements of parallel environments if its sequential version was poor. Consequently, every PRNG wrapped in Shoverand must satisfy the most stringent testing battery currently available, namely BigCrush from TestU01 [START_REF] Pierre | Testu01: A c library for empirical testing of random number generators[END_REF]. PRNGs that pass all those tests are referred to as "Crush-resistant" in [START_REF] Salmon | Parallel random numbers: as easy as 1, 2, 3[END_REF]. While being Crush-resistant cannot ensure a perfect randomness of the considered pseudorandom stream, it is a satisfying property that few PRNGs can be proud of.

Additionally, the retained algorithms must support a reliable technique to distribute numbers in a parallel environment. We have previously surveyed such techniques in [START_REF] David | Distribution of random streams for simulation practitioners[END_REF], but only some of them can be applied on a GPU platform [START_REF] Passerat-Palmbach | Pseudo-random streams for distributed and parallel stochastic simulations on gp-gpu[END_REF]. The chosen ones are then ideal candidates to be ported to GPU, if not available yet, and moreover to be integrated in Shoverand. We detail hereafter the PRNGs that are currently, or will soon be, embedded in Shoverand.

MRG32k3a

Introduced by Pierre L'Ecuyer in [START_REF] Pierre | Good parameters and implementations for combined multiple recursive generators[END_REF], MRG32k3a is particularly suited to parallelization among small computational elements such as threads thanks to its intrinsic properties. This PRNG's lightweight data structure only stores 6 integers to handle its state. The algorithm itself is quite short, and relies on simple operations to issue new random numbers. The parameters chosen for MRG32k3a are such that it has a full period of 2 191 numbers. This period is fairly enough since L'Ecuyer suggests that periods between 2 100 and 2 200 are highly sufficient even for large-scale simulations. MRG32k3a has been designed to produce independent streams and sub-streams from its original random sequence thanks to its parameters that enable safe Sequence Splitting [START_REF] David | Distribution of random streams for simulation practitioners[END_REF]. The internal parameters split the initial sequence into 2 64 adjacent streams of 2 127 random numbers, themselves divided into sub-streams containing 2 76 elements. Now considering the distribution aspect, we can assign a stream or a sub-stream to each computational element according to the Sequence Splitting technique. As long as we are focusing on parallel applications that are CUDA-enabled, we are dealing with fine-grained Single Instruction, Multiple Threads (SIMT) applications. It means that the computational elements are, in our case, the logical threads of a CUDA kernel and the principle of SIMT is to load the device with as much threads as possible. Still, we do not expect having to deal with more than 2 64 parallel threads, which is the total number of independent streams bearing 2 127 random numbers each that MRG32k3a can provide.

TinyMT

TinyMT is the latest offspring from the Mersenne Twister family. TinyMT is not described in any scientific article yet, but information about it can be found on its dedicated webpage [START_REF] Saito | [END_REF]. This PRNG matches the requirements we have formulated for a PRNG to be integrated into Shoverand: it is described as producing a good quality output, according to TestU01 statistical tests, and displays a long-enough period of 2 127 numbers. TinyMT leverages parameterization to provide highly independent streams, each stream being represented by a unique parameterized status.

The theoretical aspect of this approach is very satisfying, but TinyMT parameterized statuses need to be precomputed by a piece of software called Dynamic Creator (DC), which is shipped with the PRNG as an open-source binary. The idea here is to initialize each computing elements with a different status, since DC can create over 2 32 × 2 16 independent statuses. However, memory footprint considerations forced us to propose a hybrid implementation where the same independent parameterized status is shared among all the threads of a CUDA block. Independence between random sources is achieved by feeding each thread with a sub-stream of the original stream, following the Sequence Splitting technique. To do so, the original stream is sliced in equal chunks whose starting point, the seed status, is indicated to threads depending on their identifier. As a consequence, each thread will always consume the same random sequence in different replications of the same execution, thus ensuring reproducibility of the experience.

Philox and Threefry

Philox and Threefry are counter-based PRNGs [START_REF] Salmon | Parallel random numbers: as easy as 1, 2, 3[END_REF] also relying on parameterization to solve random streams partition concerns. Like any other PRNG considered in this study, they are Crush-resistant and display good performance in regards to their low memory footprint and high numbers throughput. They appear to be better suited than TinyMT (or any other member of the Mersenne Twister family) to a straightforward GPU implementation since their parameters are formed by a single key that can be set at runtime according to each thread's unique identifier. Please note that the GPU implementation of these PRNGs is directly provided by their authors. Both CUDA and OpenCL implementations are proposed for Philox and Threefry.

Case study: generating pseudorandom numbers in a CUDA kernel with Shoverand

In this section, we describe a major aspect of Shoverand: its user-friendly interface. We will see that on both host and device sides, our API is very expressive while remaining quite concise. Shoverand competes with two major counterparts in the CUDA world, both coming from an NVIDIA initiative, named Thrust [START_REF] Hoberock | Thrust: A parallel template library[END_REF] and cuRand [START_REF] Spanier | CUDA CURAND Library[END_REF]. These two libraries are also providing random number generation features but vary from Shoverand on several points that we will compare in this section.

Host side: Initialization phase

From the end-user's point of view, Shoverand requires an initialization phase in order to allocate its internal data structures on the device and perform some initializations. As a matter of fact, depending on the chosen PRNG, initialization might involve external data to be read from a parameter file, as is the case with TinyMT for instance.

We previously saw that we needed to consider the distribution technique and the PRNG algorithm as a pair. As a consequence, distribution techniques vary from one PRNG to another in Shoverand, but their initialization phase require the same data, which is basically the number of CUDA blocks that the kernel using the PRNG will spawn. This data being stored prior to the kernel call, no superfluous parameter needs to be passed to the kernel. This feature allows users not to have their hands tied when designing their kernels, since Shoverand does not impact kernels' prototypes like other libraries do.

As a result, the host side initialization phase boils down to a single call to a static method named init, which must be provided by every PRNG implementation to satisfy Shoverand's rules.

Device side: Computation phase

Using the device side of Shoverand is even simpler than the host side. You only have to create an instance of the PRNG you want to use and let its class' constructor do the rest. Device side initializations are performed behind the scenes by the constructor so that users have nothing to do. Then, you can pick random numbers by calling the next() method on the previously created object. If you are used to random number generation facilities offered by high-level languages such as Java, making use of Shoverand in your kernel is really intuitive.

Comparison with Thrust and cuRand

Thrust and cuRand are two projects developed by NVIDIA fellows. While cuRand is part of the CUDA SDK, Thrust is an external open-source library that can be downloaded from an Internet repository. In the paper originally introducing Shoverand [START_REF] Passerat-Palmbach | Shoverand: a model-driven framework to easily generate random numbers on gp-gpu[END_REF], we had surveyed these two libraries and identified their major drawbacks that led us to design Shoverand. The following lines investigate the changes brought by the state-of-the-art versions of Thrust and cuRand.

cuRand is NVIDIA's solution to random number generation on GPU. In [START_REF] Passerat-Palmbach | Shoverand: a model-driven framework to easily generate random numbers on gp-gpu[END_REF], we mentioned that cuRand suffered from the poor statistical quality of the PRNGs it embedded. The last version of this library partially solves this problem by following the advice we made in our previous paper. Now, cuRand embeds renowned high-quality PRNGs such as MRG32k3a [START_REF] Pierre | Good parameters and implementations for combined multiple recursive generators[END_REF] and MTGP [START_REF] Saito | Variants of mersenne twister suitable for graphic processors[END_REF], whose pseudorandom streams are stated as "Crush-resistant".

On the other hand, cuRand's API remains poor and forces users to add extra parameters in the prototypes of every kernel taking advantage of the library. The C API is not generic and loses its consistency when non-default options are enabled: for instance, MTGP's initialization step is achieved through a dedicated call in cuRand that is totally irrelevant when used with another PRNG. This approach is not convenient when you want to quickly swap PRNGs to study the impact of various random sources on a given application.

Thrust is an open source GP-GPU-enabled general purpose library developed by NVIDIA fellows. This project aims at providing a GP-GPU-enabled library equivalent to some classic general-purpose C++ libraries, like STL or Boost. Classes are split through several namespaces, such as Thrust::random. 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++ class template. We find a Linear Congruential Generator (LCG), a Linear Feedback Shift (LFS) and a Subtract With Carry (SWC), which are widely reckoned as not adapted to High Performance Computing.

Still, Thrust offers a nice API that mirrors Boost's. The random namespace provides userfriendly features very close to Shoverand's abilities. For instance, neither explicit initializations nor parameters are required in order to benefit from random number generation facilities in a kernel.

As a conclusion, we have on the one hand cuRand, a library that is improving after having been criticized by the research community for the statistical characteristics of its embedded PRNGs, and that exposes a restrictive API; and on the other hand, we have Thrust and its nice "à la Boost" API, yet powered by poor quality PRNGs. Then, Shoverand is based upon the good achievements from these two libraries: indeed, it exposes a user-friendly API while integrating only Crush-resistant PRNGs.

Listing 1 shows off how to pick up pseudorandom numbers from Shoverand. In this code snippet, we consider both the initialization on the host side and the random number generation on the device side:

Case study: embedding a new PRNG into Shoverand

Shoverand is not only a library but also a framework that allows PRNG developers to insert their own proposals as long as they follow some rules. In order to help them in their task, Shoverand employs a mechanism called concept checking that lets us express constraints in Shoverand's code that will be checked at compile-time for all the classes that inherit from ours. Concept checking is a generic programming feature available through different implementations with C++. It was introduced by [START_REF] Siek | Concept checking: Binding parametric polymorphism in c++[END_REF] and implemented in the Boost Concept Check Library (BCCL). The application scope of concept checking is quite wide and goes from interface checking (i.e., verifying the presence of a given method within a class), to assessing the presence of a particular member from a given type in a class.

Such a mechanism forces developers to match our interface without having to introduce costly runtime techniques like polymorphism and consequently virtual methods.

Thanks to the BCCL, we are able to design constraints that will make any compilation attempt fail if they are not met. In Shoverand, we force every PRNG algorithm class to inherit from our base RNG class. Then, each user-defined subclass must define at least 3 methods: init() and release(), which deal with parameters allocation and initialization from the host side, and next(), which picks up the next random number within a kernel, on the device side. In the same way, BCCL also permits us to verify that developers have provided two members to their class: the seed and parameterized statuses. These two members respectively represent the current state of the PRNG and its initial parameters. Fig. 1 sketches a UML class diagram of the expected content of a PRNG class to be embedded in Shoverand:

Conclusion

In this paper, we demonstrate how parallel stochastic simulations can be seriously impacted by the quality of the underlying Random Number Generators (RNGs), and the partitioning techniques used to feed such applications [START_REF] David | Distribution of random streams for simulation practitioners[END_REF]. Still, we have at our disposal a set of very reliable sequential generators and a set of distribution techniques adapted to the different kinds of generators.

From a practitioner's point of view, we have described Shoverand: a CUDA library that embeds the PRNGs displaying the best statistical quality. Moreover, its user-friendly interface make it a very serious choice when faced to its counterparts.

Shoverand is freely available for download at the following URL: http://http://forge. clermont-universite.fr/projects/shoverand.

Listing 1 :

 1 Example of use of Shoverand's API u s i n g shoverand : :RNG; u s i n g shoverand : : MRG32k3a ; // k e r n e l u s i n g Shoverand __global__ v o i d f o o K e r n e l ( f l o a t * ddata ) { RNG < f l o a t , MRG32k3a > rng ; ddata [ blockDim . x * b l o c k I d x . x + t h r e a d I d x . x ] = rng . next ( ) ; } . . . RNG< f l o a t , MRG32k3a >:: i n i t ( block_num ) ; f o o K e r n e l <<< block_num , thread_num >>>(d_data ) ; RNG< f l o a t , MRG32k3a >:: r e l e a s e ( ) ;

Figure 1 :

 1 Figure 1: UML class diagram of the expected interface of a PRNG in Shoverand

Acknowledgements

The authors would like to thank Luc Touraille for his careful reading and useful suggestions on the draft. This work received financial support from the Auvergne Regional Council.


Jonathan Passerat-Palmbach

David R C Hill

OpenCL: A suitable solution to simplify and unify High Performance Computing developments

Keywords: OpenCL, Abstraction layers, Survey, Source-to-source transformation, Manycore

Manycore architectures are now available in a wide range of HPC systems. Going from CPUs to GPUs and FPGAs, modern hardware accelerators can be exploited using heterogeneous software technologies. In this chapter, we study the inputs that OpenCL offers to High Performance Computing applications, as a solution to unify developments. In order to overcome the lack of native OpenCL support for some architectures, we survey the third-party research works that propose a source-to-source approach to transform OpenCL into other parallel programming languages. We use FPGAs as a case study, because of their dramatic OpenCL support compared to GPUs for instance. These transformation approaches could also lead to potential works in the Model Driven Engineering (MDE) field that we conceptualize on this work. Moreover, OpenCL's standard API is quite rough, thus we also introduce several APIs from the simple high-level binder to the source code generator that intend to ease and boost the development process of any OpenCL application.

Introduction

At the end of the last decade we entered more deeply into the era of hybrid computing, where different types of processors are used in a common computing system. Modern accelerators are now available as COTS (Commercial Off The Shelf), we can buy and or design our own high performance hybrid architecture. Manycore CPUs (Central Processing Units), GPGPUs (General-Purpose Graphical Computing Units) and more recently the AMD APUs (Accelerated Processing Units) are gaining more and more interest. FPGAs have also been considered to be elements of a general purpose High Performance Computing system, they are used in some of the best architectures (like the one proposed by Cray Inc.) [START_REF] Craven | Examining the viability of FPGA supercomputing[END_REF]. Nowadays, the High Performance Computing (HPC) community tends to champion GPU architectures to speed-up applications, but according to [START_REF] Williams | Computational density of fixed and reconfigurable multi-core devices for application acceleration[END_REF], FPGAs offer better performances than GPUs when faced with some integer arithmetic operations. In addition to impressive performances, the FPGA approaches have many advantages including small footprint with very low energy costs. Most of FPGA logic is configurable, and this architecture is saving the processing power since we only design the specific features that will be used. If we can reuse the same silicon for different functions, we also have the advantage of an easy mapping of different operations to different silicon thus allowing massive parallelism. Long before the GPGPU boom, in 2004, companies joined to establish the FPGA High Performance Computing Alliance (FHPCA). This alliance is dedicated to the use of Xilink FPGAs for high end real-world industrial applications. Since the advent of the Virtex 4 (and now with much larger devices) the use of FPGAs for general purpose numerical computing has become a reality. Even if it was not really noticed in the Top500 at that time, the University of Edinburgh supercomputing centre proposed in 2007 an FPGA based supercomputer named "Maxwell".

The case of FPGAs in HPC is particularly significant. Given the advantages previously exposed, the question is why don't we have more FPGAs in HPC? The main cause in our opinion is that the programming model problem is neither mature nor accessible to the scientific community. Common HPC developers are not used to programming in Hardware Description Languages (HDLs, such as VHDL) for FPGAs, and the problem is not just learning a new language, since it is mainly a design issue. Indeed, we do not design circuits as we design programs. Even if predefined libraries are available, they do not compensate the fact that there is no real standard hardware. The point is that we need a technology that can be understood by the heterogeneous HPC community, allowing them to address the wide range of hardware accelerators at their disposal.

Since its introduction in June 2008, the OpenCL™ (Open Computing Language) standard has been specifically designed to address the challenge of executing HPC applications on heterogeneous systems consisting of multi-core processors, GPGPUs and FPGAs. OpenCL is based on the C programming language but it provides additional features to leverage highly-parallel hardware accelerators. OpenCL allows the programmer to explicitly specify and control parallelism.

Though numerous accelerator vendors started to support OpenCL for CPUs and GPUs, fewer efforts were given to support hybrid systems with FPGAs. Recent advances from the Altera Corporation, make us confident that the OpenCL programming model will propose advantageous solutions to common problems that refrained the adoption of FPGAs for HPC. While native tools are not reliable yet, we can still focus on third-party developments that act as converters from OpenCL to HDL without the need of any further development from FPGAs vendors.

Part of this work actually considers FPGAs as a case study of the ability to take advantage of OpenCL on platforms that do not officially support this technology. We have chosen to tackle this problem through FPGA-aimed proposals, since this community has been very productive in terms of solutions to detach from vendors. Still, the same conclusions could be applied to other architectures lacking OpenCL support, the most notable being Intel Sandy/Ivybridge GPU parts, when used under the Linux OS (Operating System).

Vendor-agnostic approaches should democratize the use of OpenCL for HPC developments, and promote GPU computing at the same time, as the main platform offering OpenCL support. However, OpenCL suffers from a highly verbose and constrained API. Thus, this study also surveys software solutions assisting OpenCL development.

In this chapter, we:

• State the need of an abstraction layer for OpenCL;

• Survey the propositions to provide OpenCL support for any architecture through a source-to-source transformations approach;

• Study the most significant Application Programming Interfaces (APIs) to ease OpenCL development and compare their features;

• Recall the inputs of the OpenCL technology in High Performance Computing applications, and especially for GPUs.

2 On the need for an abstraction layer for OpenCL

OpenCL in brief

OpenCL is a standard proposed by the Khronos group that aims to unify developments on various kinds of hardware accelerators architectures like CPUs, GPUs and FPGAs. It provides programming constructs based upon C99 to write the actual parallel code (called the kernel). Kernels execution is executed by several work-items, that will be mapped to different execution units depending on the target: for instance, GPUs will associate them to local threads. For scheduling purposes, work-items are then bundled into work-groups each containing an equivalent amount of work-items. They are enhanced by APIs (Application Programming Interface) used to control the device and the execution. At the time of writing, the latest version of the API is 1.2 [START_REF]The OpenCL Specification 1.2[END_REF] that has been released in November, 2011. OpenCL programs execution relies on specific drivers issued by the manufacturer of the hardware they run on. The point is OpenCL kernels are not compiled with the rest of the application, but on the fly at runtime. This allows specific tuning of the binary for the current platform.

OpenCL brings three major inputs to HPC developments: as a cross-platform standard, it allows developing applications once and for all for any supported architecture. It is also a great abstraction layer that lets developers concentrate on the parallelization of their algorithm, and leave the device specific mechanics to the driver.

OpenCL: a constrained API

In the introduction we have evoked the major problem of OpenCL: its complicated API. Concretely this lies in three major problems that will be detailed hereafter: bloated source code, double calls and risky constructs. First, when kernel functions that execute on the hardware accelerator remain concise and thus fully expressive, host API routines result in verbose source code where it is difficult for an external reader or for a non-regular developer, to determine the purpose of the application among all those lines. Second, some design choices of the API make it even more verbose when trying to write portable code that can run on several hosts without any change. For instance, OpenCL developers are used to calling most of the API query functions twice. Indeed, a first call is often required to determine the number of results to be expected, and a second actually gets the right amount of elements. We will designate these two subsequent invocations of the same routine as the double-call pattern hereafter. Finally, such verbose constructs discourage developers to check the results of each of their calls to the API. At the same time, the OpenCL API is very permissive with regards to the type of the parameters its functions accept. Now imagine that two parameters have been awkwardly swapped in an API call, it is very likely that the application keeps quiet about this error if the developer has not explicitly checked the error code returned by this call. In an heterogeneous environment such as OpenCL, where the parallel part of the computation will be relocated on hardware accelerators, runtime errors that only issues error codes are not the easiest bugs to get rid of.

All these drawbacks make it clear that although the OpenCL standard is a great tool offering portability in high performance computing, it does not meet the higher level APIs expected by parallel developers to help them avoid common mistakes, and to produce efficient applications in a reasonable lapse of time.

OpenCL support on various platforms: the FPGAs case study

Enabling OpenCL for a particular FPGA architecture can be achieved in two major ways. The first relies on vendors who need to provide an OpenCL driver implementation for their particular device. When both NVIDIA and ATi provide OpenCL drivers for their GPU products, Intel has recently enabled a large set of its processors to support OpenCL. However the driver implementation solution is still in its infancy across the FPGA community. Altera is known for their FPGA products, but is also an active member of the Khronos Group, responsible of the OpenCL standard. The company has announced the development of an OpenCL framework, and presents encouraging results of a Monte-Carlo Black Scholes simulation implementation in [START_REF]Implementing FPGA Design with the OpenCL Standard[END_REF].

Researchers' initiatives have led us to consider another way of thinking that does not involve any intervention from the FPGA's vendor, but relies on third-party developments to provide OpenCL support for FPGAs. The main advantage of these solutions is that they can be applied to any platform lacking OpenCL support. In this section, we study the most promising solutions from the literature, focusing on FPGAs, as long as this community has been in need of a unified programming model for a long time now, and proposed consequently more tools to meet their expectations.

Source-to-source transformation approach

Supporting the whole OpenCL standard can become a long and costly operation for a vendor, while its community might not perceive the benefits of switching to OpenCL developments, or simply do not require high performance. That being said, another strategy comes into play to bypass vendors decisions: offer OpenCL support for a range of devices from a third-party effort. Such projects are all based on the same principle that consists in transforming the OpenCL source code into a language already supported by the target platform. In this section, we will survey the most interesting proposals of source-to-source compiling aiming at similar goals.

At the beginning, there was CUDA

At the moment, OpenCL is mostly known to be a unified solution to program GPUs from various origins, and to compete with NVIDIA's CUDA technology in the GPU computing domain. Being almost two years older than OpenCL, CUDA has been targeted by research works earlier than its counterpart, particularly when considering the opportunities to have it running on other devices than NVIDIA GPUs. Most of these approaches make use of the same principle as the strategy we identified to provide OpenCL support for an FPGA architecture from a third-party initiative: i.e. source-to-source transformation.

In 2008, Stratton et al. proposed MCUDA, a tool intended to execute CUDA code on a multi-core CPU [START_REF] Stratton | MCUDA: An efficient implementation of CUDA kernels for multicore CPUs[END_REF]. This paper brings up the fundamental problems that will be encountered when attempting to set up such a code generation approach, and even puts forward efficient solutions. Let us recall that CUDA relies on the Single Instruction, Multiple Thread (SIMT) paradigm, so the application logic is placed at the thread level. Threads perform the same operations on different input data, but sometimes, an operation needs to access data that have been processed by another thread. In such a case, CUDA developers have to synchronize their threads at the particular step of the algorithm dealing with shared data to make sure the data they are processing is up to date.

Before parallelizing the application among CPU threads, MCUDA first serializes the CUDA input in order to analyse it more simply. MCUDA tries to identify independent code parts in terms of data within the newly serialized code. When such parts are detected, they are bounded by what is called a threadloop. Thread-loops induce implicit synchronization points at their boundaries. Thus, an explicit request to synchronize threads in the original CUDA source code will result in two thread-loops: one embracing the code before the synchronization point and the other one after.

Instead of basing the computation on logical threads, MCUDA rises to a block of logical threads level. This strategy excludes semantic losses since blocks are independent along a classical CUDA execution. They contain logical threads that will perform the same operations with no constraint on the order in which the blocks are processed. MCUDA behaves identically, since after being serialized, blocks are executed by OS threads created by OpenMP directives, and thus relying on the underlying scheduling policy for their order of execution to be determined. Now that the code is serialized and that the synchronization points are respected, MCUDA treats memory accesses. The most important point of this step is to understand that CUDA, such as OpenCL, defines several data areas mainly to take advantage of the different caching mechanisms available in a GPU architecture. As long as CPUs display a single cached memory space, MCUDA maps any kind of memory to the unique memory space of the CPU. Still, memory areas also bring information on the scope of the data they store, forcing MCUDA to apply another conversion step that again prevents semantic losses in the output code. For example, when variables shared by all the logical threads of a block can be replaced by a single variable, elements declared as private to a thread have to be duplicated. To do so, MCUDA introduces arrays containing the particular value of a given thread in the cell corresponding to its identifier. We will see further on that the same pondering about memory spaces must be observed when thinking about an FPGA implementation.

Basically, the philosophy of MCUDA lies in two main steps:

1. Group local-thread based computations and express them at a block level (SIMT → SIMB);

2. Parallelize the execution of independent blocks with a pool of worker threads.

More generally, the purpose is to first obtain independent elements before executing them concurrently following an arbitrary scheduling. This mechanism is classical and can also be found in other works involving CUDA and GPU architectures such as [START_REF] Passerat-Palmbach | Warp-Level Parallelism: Enabling Multiple Replications In Parallel on GPU[END_REF]. This implies two major choices when implementing an OpenCL converter: first, you need to figure out how the work-items will be mapped on the target architecture, and second, how the memory hierarchy will be represented on your FPGA.

One year later, some of the authors of MCUDA applied their conclusions to FPGA architectures and proposed FCUDA: a tool to compile CUDA kernels onto FPGAs. FCUDA's philosophy is quite similar to MCUDA's in the fact that its main purpose is to translate CUDA SIMT code into SIMB independent blocks that are ready to be transformed. Now, when MCUDA spreads those blocks on multiple CPU cores through OpenMP threads, FCUDA will synthesise FPGA cores according to the Register Transfer Level (RTL) design corresponding to the kernel function.

The way FCUDA deals with memory is quite interesting and thoroughly described in the paper [START_REF] Papakonstantinou | FCUDA: Enabling efficient compilation of CUDA kernels onto FPGAs[END_REF], wherein authors make concrete proposals concerning the memory locations of the FPGA assigned to their CUDA counterparts. For instance, they propose to take advantage of Direct Memory Access (DMA) burst mode in order to implement constant memory.

Additionally, the workflow allowing FCUDA to transform CUDA code into FPGA Hardware Description Language (HDL) requires an intermediate stage where C code is produced and passed to a third-party software. This step enables developers to tune the resulting code before it is synthesized into a RTL design. However, we need to note that FCUDA displays some drawbacks, indeed it only handles kernels but does not take into account the host API.

Then comes OpenCL

The MCUDA way to deal with parallel applications was pursued by several recent studies, this time generating HDL for a given FPGA from OpenCL source code. In 2011, Owaida et al. introduce SOpenCL (Silicon OpenCL) a source-to-source transformation process that converts OpenCL code into HDL to build FPGA-enabled hardware accelerators [START_REF] Owaida | Synthesis of platform architectures from opencl programs[END_REF]. Just like MCUDA, SOpenCL first applies a serialization step on the kernel in order to join logical threads from a given work-group into independent execution environments. SOpenCL uses the fact that work-groups are actually 3D containers of local-threads to insert a triple-nested loop that will pass on each local-thread thanks to its x, y, z coordinates in the kernel to compute the data it was assigned.

Barriers and synchronization points are also handled in a similar way. Whenever a synchronization point must be set up, the loop around the current statement is split into two sub-loops: the first deals with the part of the statement before the synchronization point, while the second handles the remaining part. This process is called loop fission.

Finally, memory areas are also treated in a way close to MCUDA. The idea is to figure out the lifetime of a variable to state whether it can be reused across loops or not. Concerning the variables' scope, each logical-thread is provided with a private copy of the variable only if the latter is meant to be used in several sub-loops, which would potentially overwrite its content. The final location of variables in memory is bound to the targeted FPGA, but we have seen previously that no semantic information was lost during this stage.

What makes SOpenCL original is its template-based approach when the time comes to generate HDL code. SOpenCL proposes a hardware accelerator template including both the Data Path and Stream Units that will be involved in the computation. Taking advantage of the LLVM compiler framework [START_REF] Lattner | LLVM: A compilation framework for lifelong program analysis & transformation[END_REF], the hardware generation flow results in synthesisable HDL for the accelerator following an architectural template, instead of issuing an unpredictable output code.

In 2010, Jääskeläinen et al. proposed another work inspired from MCUDA. However, this OpenCLbased design methodology [START_REF] Jääskeläinen | OpenCL-based design methodology for application-specific processors[END_REF] differs from other proposals by exploring the extension feature offered by the OpenCL standard. As we have seen previously, OpenCL allows vendors to add specific extensions from their own to exploit their hardware a bit further than the sole standard API would. The main advantage of this capacity is that it does not imply rewriting code when switching from an extension-enabled architecture to a standard one. Thus, vendors can allow to leverage the specificities of their hardware, while remaining compliant with others'. As long as extensions availability can be checked at source level, one can adapt its algorithm behaviour to take advantage of the extensions only when possible. In an FPGA approach, this feature appears crucial since it lets hardware designers provide custom architectures synthesized on an FPGA platform. Clients can then choose to benefit from these particular developments or not, and even try both solutions in the blink of an eye without having to revise their entire source.

In conclusion, we have seen that no matter the source-to-source tool we studied, the implementation choices are roughly the same. The process follows a common skeleton that is summed up hereafter:

1. Serialize work-groups by executing their local-threads sequentially; 2. Eliminate barriers and synchronization points thanks to loop fission; 3. Analyse variables livelihood to determine which can be shared and which should be copied in a private area for each local-thread;

4. Convert the resulting source code to an HDL, taking work-groups as the base unit for the computation.

4 High-Level APIs for OpenCL: Two Philosophies 4.1 Ease OpenCL development through high-level APIs

Standard API C++ Wrapper

The OpenCL standard not only describes the C API that any implementation should meet, but also the C++ wrapper that comes along with it. This wrapper provides both C++ bindings of the OpenCL calls, and also two restricted declinations of the Standard Template Library (STL) from genuine C++: the string and vector classes. These two classes have been rewritten for the purpose of the OpenCL wrapper as a subset of their counterparts but are mostly compliant. The idea is to get rid of STL's bloated classes, which std::string is nothing but the best example. As long as they display an interface close to the original classes, the OpenCL versions, stored in the cl namespace, can be swapped with the matching STL class thanks to a single macro.

Another C++ mechanism nicely leveraged by the OpenCL C++ wrapper is exceptions. OpenCL errors are handled in a low-level way with error codes to be compared to 0 to assess whether the previous call is completed successfully. Obviously this technique is not quite adapted to develop with higher level languages, because it inflates source codes with non-effective lines. Moreover, it forces developers to be very careful and to explicitly check the returning code of their invocations, and, in case of problems, to retrieve the corresponding error. In our case, this painful process is handled by the exception mechanism that forces developers to catch the exception and treat it consequently.

The remaining elements of this binding are dedicated to ease the OpenCL API for developers. The most significant example in this way is the double call pattern foreseen in our short introduction to OpenCL in section 2.2. Using the traditional OpenCL C API, one needs to perform two successive calls to the same function to first figure out the number of results to be stored in an array provided for this purpose, before actually filling it through a second call requesting the exact number of elements to store in the array of results subsequently. The C++ wrapper here highly facilitates the process since a single call is needed to obtain the results, while cl::vector is used instead of dynamically allocated arrays. Listing 1 shows an example of the wrapper syntax. Although efforts can be noticed to provide a simple API, the result remains quite verbose since Listing 1 only lists GPU devices and creates a dummy kernel from a source file in the context of the discovered device: Used to provide nice bindings for C++ development tools going from database drivers to concurrency, the Qt library developed by Nokia [START_REF] Nokia | QtOpenCL[END_REF] now also offers its own OpenCL wrapper as a set of extensions named QtOpenCL. At the time of writing, we would like to insist on the fact that QtOpenCL is neither included in the stable Qt 4.7 release, nor in the future 4.8 version according to Nokia's roadmap. However, QtOpenCL is freely available for download as an extension and is compatible with Qt 4.6.2 and Qt 4.7 frameworks.

As an OpenCL wrapper, QtOpenCL aims to break down three main barriers that we already identified as OpenCL drawbacks: initialization phase, memory management and kernel execution. To assist users in the initialization phase, and particularly to compile the external OpenCL source code, QtOpenCL handles automatically OpenCL source files reading through a single method. Memory management receives as much consideration so that buffers can be created more simply. Moreover, their object oriented interface allow users to call methods such as read() directly on buffers to retrieve results from a hardware accelerator.

The third point of the QtOpenCL intentions targets kernel executions. This aspect is handled "à la Qt" thanks to the QFuture class, well-known by developers using Qt for their CPU developments involving concurrency.

QFuture is an event returned when a kernel is run through the QtConcurrent::run method, and allows to wait for the kernel to complete. In addition, QFuture is compatible with the QFutureWatcher class that uses the signal/slot mechanism which Qt is based upon. The latter mechanism is an implementation of the Observer design pattern that causes a routine to be called when a particular event occurs. In our case, the event signals the completion of the associated kernel and can rise an update of the Graphical User Interface (GUI) of the application for example.

Listing 2 presents how to enqueue a dummy kernel that takes a vector of integers as a parameter and outputs another one, using QtOpenCL: Listing 2: GPU context creation kernel enqueueing using QtOpenCL / / c o n t e x t c r e a t i o n QCLContext c o n t e x t ; i f ( ! c o n t e x t . c r e a t e ( ) ) { s t d : : c e r r << " E r r o r i n c o n t e x t c r e a t i o n f o r t h e GPU" << s t d : : e n d l ; r e t u r n 1 ; } c o n s t i n t v e c t o r S i z e = 1 0 2 4 0 0 0 ; QCLVector < i n t > i n V e c t o r = c o n t e x t . c r e a t e V e c t o r < i n t > ( v e c t o r S i z e ) ; QCLVector < i n t > o u t V e c t o r = c o n t e x t . c r e a t e V e c t o r < i n t > ( v e c t o r S i z e ) ; f o r ( i n t i = 0 ; i < v e c t o r S i z e ; ++ i ) {

PyOpenCL

PyOpenCL is a research initiative developed at the same time as its CUDA counterpart PyCUDA [START_REF] Klöckner | PyCUDA and PyOpenCL: A Scripting-Based Approach to GPU Run-Time Code Generation[END_REF]. Both approaches rely on a concept called GPU Run Time Code Generation (RTCG) that intends to solve common problems encountered when harnessing hardware accelerators. In PyOpenCL, this translates in a flexible code generation mechanism that can adapt to new requirements transparently. In the end, programmers can summon kernels as if they were methods from the program instance they have just built.

Apart from those dynamic features, PyOpenCL displays the classical inputs from an OpenCL wrapper developed in a high level language such as Python. The latter being well-known for its concision and simplicity, PyOpenCL directly benefits from this characteristic. This is achieved widely because of Python dynamic typing system that delegates an important part of the work to the interpreter. Finally, OpenCL source code reading takes advantage of Python file handling capacities, so that a program can be read from source and built using no more than four lines, as it can be seen in Listing 3:

Generating OpenCL source code from high-level APIs

ScalaCL

ScalaCL is a project, part of the free and open-source Nativelibs4java initiative, led by Olivier Chafik [START_REF] Chafik | ScalaCL[END_REF]. The Nativelibs4java project is an ambitious bundle of libraries trying to allow users to take advantage of various native binaries in a Java environment.

From ScalaCL itself, two projects have recently emerged. The first one is named Scalaxy. It is a plugin for the Scala compiler that intends to optimize Scala code at compile time. Indeed, Scala functional constructs might run slower than their classical Java equivalents. Scalaxy deals with this problem by preprocessing Scala code to replace some constructions by more efficient ones. Basically, this plugin intends to transform Scala loop-like calls such as map or foreach by their while loops equivalents. The main advantage of this tool is that it is applicable to any Scala code, without relying on any hardware.

The second element resulting from this fork is the ScalaCL collections. It consists in a set of collections that support a restricted amount of Scala functions. However, these functions can be mapped at compile time to their OpenCL equivalents. A compiler plugin dedicated to OpenCL generation is called at compile time to normalize the code of the closure applied to the collection. Functional programming usually defines a closure as an anonymous function embedded in the body of another function. A closure can also access the variables from the calling host function. The resulting source is then converted to an OpenCL kernel. At runtime, another part of ScalaCL comes into play, since the rest of the OpenCL code, like the initializations, is coupled to the previously generated kernel to form the whole parallel application. The body of the closure will be computed by an OpenCL Processing Element (PE), which can be a thread or a core depending on the host where the program is being run. Listing 4 presents a simple closure that computes the apply the cosine function to every element of an array. Only the execution of the closure's body is deported to the hardware accelerator targeted by OpenCL: The obvious asset of ScalaCL is its ability to generate OpenCL at compile time. It means that we could enhance ScalaCL by adding an extra step to tune the issued OpenCL kernel at compile time and take advantage of a GPU vendor specific extensions for instance.

Aparapi

Aparapi [START_REF] Amd | Aparapi[END_REF] is a project initiated by AMD and recently freed under an MIT-like open source licence. It intends to provide a way to perform OpenCL actions directly from a pure Java source code. This process involves no upstream translation step, and is then wholly carried out at runtime. To do so, Aparapi relies on a Java Native Interface (JNI) wrapper of the OpenCL API that hides complexity to developers.

Basically, Aparapi proposes to implement the operations to be performed in parallel within a Kernel class. The kernel code takes place in an overridden run() abstract method of the Kernel class that sets the boundaries of the parallel computation. At the time of writing, only a subset of Java features are supported by Aparapi, it means that run() can only contain a restricted amount of features, data types, and mechanisms. Concretely, primitive data types, except char, are the sole data elements that can be manipulated within Aparapi's kernels.

For years, Java has been used to managing concurrency problems thanks to a Thread class whose implementation makes use of the native thread of the underlying Operating System where the Java Virtual Machine (JVM) runs. Thread implements the Runnable interface that only consists in providing a void run() method that will be called when the thread is launched. This mechanism is widely adopted in the Java world, and most of the frameworks offering an abstraction layer to Thread employ it in order to remain compliant between each other. At first, it seems that Aparapi follows the same path given that it designates a run() method to contain the parallel code. However, Aparapi's source code does not involve Runnable at any moment. A simple example of using Aparapi is set up in Listing 5 to square an array of 8 integers: Listing 5: Squaring an array of integers using Aparapi f i n a l i n t [ ] i n A r r a y = new i n t [ ] { 1 , 2 , 3 , 4 , 5 , 6 , 7 , 8 } ; f i n a l i n t [ ] o u t A r r a y = new f l o a t [ i n p u t . l e n g t h ] ; This design choice seems rather awkward when taking into consideration that the parallel tasks are assigned to a pool of CPU worker threads, the Java Thread Pool (JTP), when no accelerator can be found on the host platform. In fact, several implementations of a JTP are now shipping with Java Standard Development Kit (SDK) like Executors or Fork/Join, and have proved to be both efficient and userfriendly. Software design fosters reutilization as much as possible and being compliant with standard tools not only fastens development but it also strengthens it.

Additionally, it is interesting to note that Aparapi is not fully cross-platform, albeit being written in Java. The JNI bindings used to perform the calls to the OpenCL API make the Java package of Aparapi bound to native binaries. Thus, Aparapi cannot be shipped as a single package and depends on the ability of its underlying platform to run native code. This aspect can become a problem in terms of simplicity of use by clients, since it forces them to have a C++ tool-chain installed on each platform they want to run Aparapi on, so that the native part of the library can be recompiled for their particular platform. Still, native code dependency is not a major drawback for Aparapi given that there are few chances that a platform cannot execute the compiled-side of the JNI part of the library. This means that Aparapi can mostly be considered as a cross-platform tool that requires little effort from the client to be effective.

The first versions of Aparapi came with a slight limitation: AMD constrained the library to run on its devices only! Systems where AMD devices could not be found used to fall back to the JTP instead.

Thanks to the open-source licence that protects this tool, a third-party developer was able to remove that latter constraint. In fact, execution is locked on a particular set of devices, here AMD GPUs, only by comparing the OpenCL platform identifier to the string identifying an AMD OpenCL platform. The removal of this software restriction makes Aparapi a potential high-level API to design HPC application.

A complete solution: JavaCL

JavaCL is an open source initiative that targets a smooth integration of OpenCL in Java applications. Like ScalaCL that we have studied in the previous section, JavaCL is part of the Nativelibs4java project developed by Olivier Chafik. The JavaCL library displays two interesting aspects that we will describe hereafter.

The first input from this library is to ease OpenCL development from the host side. In order to be easily integrated in Java code, the OpenCL 1.1 API [START_REF]The OpenCL Specification 1.2[END_REF] has been fully wrapped in Java classes. Every element required to perform OpenCL computations can now be handled directly in Java. This allows us to write nothing but the kernel using the C OpenCL API.

JavaCL is issued as a single jar file containing all the dependencies needed to be executed on any platform, included native libraries in charge of being the bridge between Java and the OpenCL platform. This way, it is perfectly independent from its underlying platform, and can be easily shipped by several ways. The most convenient one in our opinion is a Maven repository. Maven [START_REF]Apache Maven Project[END_REF] is a build tool targeting Java applications that allows the projects it constructs to declare dependencies located in remote repositories. When a dependency is encountered for the first time, the maven client installed on the host system building the project downloads the dependency, and stores it in a local repository, so that it can be reused in a further build.

Functions intending to query an OpenCL installation have been designed so as to provide the most information in a single call. This completeness leads to an invocation pattern well-known by OpenCL developers: the double-call pattern foreseen in section 2.2. For instance, a function such as clGetPlat-formIDs will return a list of the available OpenCL platforms in an array passed as a parameter along with its size. This array must have been allocated upstream and its size is consequently the maximum amount of results it can store. In addition to filling the array with the available platforms, clGetPlatformIDs also returns the total number of platforms in the system. Consequently, an application needs to invoke clGet-PlatformIDs twice in order to figure out dynamically the amount of platforms on a given system: first, the function is called to determine the number of platforms and allocate an array according to this result, second, the function is summoned to actually fill the array. Such a design is widely used in the OpenCL API, and developers often have to call the same function twice in a row to make their code portable. The whole process results in bloated source files, whose real behaviour might become difficult to comprehend. JavaCL answers this awkwardness through its API that automatically issues a filled array as a result.

When platforms have been identified, it is time for an OpenCL application to list their belonging devices and to create a context combining some of them. As long as a high-performance computing application's main purpose is to speed up computation time, developers usually select the most efficient devices to form the OpenCL context in which their application will run. Once again, JavaCL takes care of this step by providing a createBestContext method that creates a context with the device containing the most work units.

The two previous features focused on increasing the ease of use of the OpenCL API. Moreover, JavaCL also enhances OpenCL development thanks to the Java language capabilities. Let us now examine the characteristics that make JavaCL really more than a simple OpenCL wrapper.

As we have seen in our introduction to OpenCL, the latter stores data that are to be treated on the device into buffers. The main concern with these buffers is that whatever the data they contain, they are represented by the sole cl_mem type. This is a potential source of harm since the compiler has no way to check that you are passing a wrong buffer type as parameter to a kernel for example. On the other side, object oriented languages such as Java are used to encapsulating data in dedicated containers. JavaCL does so by providing a different class for each data type that contains a buffer. Thanks to Java generics, a CLBuffer class can be parameterised with the primitive data type the buffer actually contains. Not only expressive containers is a nice feature for developers who have to correct a source code they did not write, but it is a particularly good point that helps compilers finding errors. Strong static typing prevents errors encountered when buffers from any type are accepted as parameters, which can appear dramatic when programming devices such as GPUs or FPGAs. To sum up, the JavaCL API allows the compiler to check little programming errors at compile time, before they turn into malicious bugs at runtime.

Being compiled at runtime, within the host program, OpenCL sources might lead to compiling errors as any program would. However, due to on-the-fly compiling, errors are a bit more tedious to take into account than in a usual program. Developers have to explicitly request for the compilation error log through the classical double call syntax inherent in the OpenCL API query system. This is not particularly suited to solve problems efficiently, and one has to make sure to correctly handle the output of the compiler anytime he builds a kernel through the C API. In the Java world, runtime errors are traditionally handled by exceptions. The exceptions mechanism is cleverly adapted by JavaCL so that any kernel build will rise a CLBuildException that can be caught, and whose content can be printed easily.

In the same way, enqueued tasks errors are also wrapped in JavaCL exceptions. An exception hierarchy has been designed to cover the scope of potential errors issued by an OpenCL-enabled application. Explicit exceptions matching the most spread OpenCL errors can be casted throughout the application lifetime; for example, memory-bound errors corrupting the command queue are covered by CLException.InvalidCommandQueue. Once again, this JavaCL feature avoids the developer having to deal with involved errors codes, while the way Java handles exceptions ensures that errors will be reported at one time or another.

Last but not least, the most promising part in our opinion is the Generator. All the previous features introduced in this section either wrapped OpenCL calls or enhanced the API with type and error checking capabilities for example. The interesting point with the JavaCL generator is its ability to parse a pure OpenCL source file containing one or several kernels, and to issue a corresponding Java source file. The latter contains a class that associates a method per kernel present in the source file. In order to smoothly integrate in the host application, the automatically generated Java class is named after the source file name, whereas each of its methods bears the name of a kernel parsed in the file. Two aspects need to be distinguished to fully understand the inputs the generator brings to OpenCL development.

First, it widely contributes to simplifying OpenCL development: while JavaCL wrappers already make OpenCL on-the-fly compilation comfortable, the Generator allows developers to skip this step by reading the source code at compile time and generating a corresponding class whose methods can be called to launch the associated kernel.

Second, applications safety is again enforced given that kernel parameters are now typed. Traditionally, kernel parameters need to be set once and for all before enqueueing the kernel. JavaCL brings a first improvement to this process by allowing parameters to be all passed at once as a list of Objects. However, this behaviour might lead to runtime errors since the type and order in which parameters are passed cannot be verified at compile time, each of them being identified by a reference to an Object. Here, the Generator enables an earlier detection of this kind of error since the prototype of the method acting as a proxy of the kernel to be enqueued will only match the right type and order of the parameters. Listing 6 displays the execution of an instance of the kernel MyKernel, which corresponding class was generated upstream by the JavaCL Generator: Listing 6: Squaring an array of integers using JavaCL CLContext c o n t e x t = JavaCL . c r e a t e B e s t C o n t e x t ( ) ; CLQueue q u e u e = c o n t e x t . c r e a t e D e f a u l t Q u e u e ( ) ; f i n a l i n t v e c t o r S i z e = 3 6 8 6 4 ; P o i n t e r < I n t e g e r > i n P t r = a l l o c a t e I n t s ( v e c t o r S i z e ) ; f o r ( i n t i = 0 ; i < v e c t o r S i z e ; ++ i ) { i n P t r . s e t ( i , i ) ; } C L B u f f e r < I n t e g e r > i n V e c t o r = c o n t e x t . c r e a t e I n t B u f f e r ( Usage . I n p u t , i n P t r ) ; / / C r e a t e an OpenCL o u t p u t b u f f e r CLBuffer < I n t e g e r > o u t V e c t o r = c o n t e x t . c r e a t e I n t B u f f e r ( Usage . O u t p u t , w i d t h * h e i g h t * 8 ) ; one meta-model to another. Thus, this modelling way would provide a base when enabling new targets to support source-to-source compilation of OpenCL applications, allowing developers to save lots of time and to simply follow the guidelines contained in the transformation rules.

Conclusion

In this chapter we have attempted to show that the OpenCL standard was a great tool offering portability in high performance computing and authorizing a wider usage of GPUs, FPGAs, and other manycore platforms. We have seen that a raw usage of OpenCL does not meet the comfort of higher level APIs that parallel developers expect to help them avoid common mistakes, and to produce efficient applications in a reasonable amount of time. We have presented different approaches enabling software reuse and even more portability with features that avoid the developed code to be bound to the intrinsic characteristics of hardware platforms. Thus developers can base their developments upon reliable software components that have been tested many times and across several architectures. This consideration opens a new scope on the benefits of OpenCL for HPC. Using a cross-platform standard gives new testing perspectives. Indeed, such a feature obviously profits to heterogeneous computing approaches where FPGA can now be integrated among CPUs and GPUs using OpenCL. Not only do we have at our disposal abstraction levels much more understandable for developers not used to HDL-enabled codes, but we have seen in the previous sections that lots of APIs intend to bind OpenCL to widespread high-level languages such as C++ or Java. At the same time, legacy codes written in these languages can now benefit from the horsepower of manycores, without having to mix several programming languages in order to enable an application to exploit new hardware accelerators. On the sole FPGA development side, OpenCL can also bring interesting features thanks to its ability to rely on the compiler, and the vendor driver to produce the hardware dependent part of the application. By doing so, migrations to a new FPGA versions do not rely on developers anymore but rather on the vendor's drivers. The involved application should consequently scale smoothly requiring no or few interventions from developers.

Last but not least, an abstraction layer that automatically generates low-level code changes the way FPGA applications are developed. Although several Intellectual Property (IP) cores are available for various purposes such as the support of Ethernet, they still need to be routed to be enabled. In the Passerat-Palmbach & Hill same way, developers might spend a lot of valuable time to design state machines that will rule the computation on the FPGA. Albeit being important for a successful development, these considerations distract developers from the principal purpose of any HPC application: speed-up! As a conclusion, delegating low-level tasks that can be automated to the backend of an abstraction layer should improve the quality of the parallel developments, since they will take up a larger place in the roadmap of an application.

Listing 1 :

 1 GPU devices listing and kernel creation using the C++ wrapper API (adapted from[START_REF] Scarpino | OpenCL in Action[END_REF]) t r y { / / P l a c e t h e GPU d e v i c e s o f t h e f i r s t p l a t f o r m i n t o a c o n t e x t c l : : v e c t o r < c l : : P l a t f o r m > p l a t f o r m s ; c l : : v e c t o r < c l : : Device > d e v i c e s ; c l : : P l a t f o r m : : g e t (& p l a t f o r m s ) ; p l a t f o r m s [ 0 ] . g e t D e v i c e s ( CL_DEVICE_TYPE_GPU , &d e v i c e s ) ; c l : : C o n t e x t c o n t e x t ( d e v i c e s ) ; / / C r e a t e and b u i l d p r o g r a m s t d : : i f s t r e a m p r o g r a m F i l e ( " k e r n e l . c l " ) ; s t d : : s t r i n g p r o g r a m S t r i n g ( s t d : : i s t r e a m b u f _ i t e r a t o r < c h a r > ( p r o g r a m F i l e ) , ( s t d : : i s t r e a m b u f _ i t e r a t o r < c h a r > ( ) ) ) ; c l : : Program : : S o u r c e s s o u r c e ( 1 , s t d : : m a k e _ p a i r ( p r o g r a m S t r i n g . c _ s t r ( ) , p r o g r a m S t r i n g . l e n g t h ( ) + 1 ) ) ; c l : : Program p r o g r a m ( c o n t e x t , s o u r c e ) ; p r o g r a m . b u i l d ( d e v i c e s ) ; / / Add k e r n e l t o progam c l : : K e r n e l f o o K e r n e l ( program , " f o o " ) ; / / C r e a t e k e r n e l c l : : v e c t o r < c l : : K e r n e l > a l l K e r n e l s ; p r o g r a m . c r e a t e K e r n e l s (& a l l K e r n e l s ) ; } c a t c h ( c l : : E r r o r e ) { s t d : : c e r r << e . what ( ) << s t d : : e n d l ; } 4.1.2 QtOpenCL

  i n V e c t o r [ i ] = i ; } / / k e r n e l b u i l d QCLProgram p r o g r a m = c o n t e x t . b u i l d P r o g r a m F r o m S o u r c e F i l e ( " . / k e r n e l . c l " ) ; QCLKernel k e r n e l = p r o g r a m . c r e a t e K e r n e l ( " f o o " ) ; / / e n q u e u e and r u n k e r n e l k e r n e l . s e t G l o b a l W o r k S i z e ( v e c t o r S i z e ) ; k e r n e l ( i n V e c t o r , o u t V e c t o r ) ;

Listing 3 :

 3 GPU program building using PyOpenCL p r o g r a m F i l e = open ( ' f o o . c l ' , ' r ' ) p r o g r a m T e x t = p r o g r a m F i l e . r e a d ( ) p r o g r a m = c l . Program ( c o n t e x t , p r o g r a m T e x t ) p r o g r a m . b u i l d ( )

Listing 4 :

 4 Computing cosine of the 1000000 first integers through ScalaCL import s c a l a c l . _ import s c a l a . math . _ i m p l i c i t v a l c o n t e x t = new S c a l a C L C o n t e x t v a l range = ( 0 t o 1 0 0 0 0 0 0 ) . c l v a l r a n g e A r r a y = r . t o C L A r r a y / / Runs a s y n c h r o n o u s l y on t h e GPU v i a OpenCL v a l m a p R e s u l t = r a n g e A r r a y . map ( v => c o s ( v ) . t o F l o a t )

  K e r n e l k e r n e l = new K e r n e l ( ) { p u b l i c v o i d r u n ( ) { o u t V e c t o r [ g e t G l o b a l I d ( ) ] = i n A r r a y [ g e t G l o b a l I d ( ) ] * i n A r r a y [ g e t G l o b a l I d ( ) ] ; } } k e r n e l . e x e c u t e ( i n A r r a y . l e g t h ) ;

Fig. 1 Figure 1 :

 11 Figure 1: Skeleton of an MDE approach to transform OpenCL into any programming language

Passerat-Palmbach & Hill

Acknowledgements

The authors would like to thank D. S. Hill, Khaled Benkrid and Wim Vanderbauwhede for their careful reading and useful suggestions on the draft. This work was partly funded by the Auvergne Regional Council.

funded by the Auvergne Regional Council and the FEDER

Finally, the combination of an enhanced wrapping API and a generation mechanism proves that JavaCL goes further than the other APIs studied here.

Summary table of the solutions

In Tab. 1, we compare the APIs studied in this work according to the availability of three features: highlevel wrapper, code generation facility, cross-platform portability of the resulting binary that will run on the host. Tab. 1 states that although C++ is one of the most spread language in the HPC community, it suffers from poor APIs to abstract OpenCL.

Perspectives

Source-to-source solutions might lead to interesting works in the field of Model Driven Engineering (MDE). MDE is a modelling approach that emphasises two main aspects: the first are models and metamodels that allow, for instance, to represent Domain Specific Modelling Languages (DSML), and are generally expressed in the UML formalism. The second element guiding MDE are transformation engines and code generators, they bring a dynamic part in MDE. Their role is to define how meta-models are bound together, and how specific code can be generated from the information they contain. We observe at least a three-level hierarchy often presented as a pyramid, where meta-models are at the top level (M2), classical models at the intermediate one (M1), and the last level represents the real world (M0). In our case, we can consider both OpenCL as a meta-model of kernel functions. We can say that kernels are conform to OpenCL. On the other side, programming constructs such as VHDL, NVIDIA PTX or Intel SSE instructions are three equivalent kinds of models describing concrete binaries or RTL designs, i.e.: the real world. These languages obviously display specificities but their final purpose is the same, so their common aspects can be factored in a meta-model to which they conform. Let us call this metamodel Meta-Language. Given that OpenCL and Meta-Language are two meta-models, MDE allows us to transform elements of the same level automatically, provided valid transformations are available from