Benoît Pradelle
email: benoit.pradelle@inria.fr

Philippe Clauss
email: philippe.clauss@inria.fr

Vincent Loechner Camus

Adaptive Runtime Selection of Parallel Schedules in the Polytope Model

Keywords: Adaptive code selection, runtime performance evaluation, parallel loop nest, polytope model

There is often no unique version of a program that provides the best performance in all circumstances. Compilers should rely on an adaptive runtime decision to choose which optimizing and parallelizing transformations will lead to the best performance in any execution context. We present a new adaptive framework solving two drawbacks of existing methods: it is effective since the very first execution, and it handles slight variations of input data shape and size.

In our proposal, different code versions of parallel loop nests are statically generated by the compiler. At install time, each version is profiled in different execution contexts. At runtime, the execution time of each code version is predicted using the profiling results, the current input data shape and the number of available processor cores. The predicted best version is then run.

Our framework handles several versions of possibly tiled parallel loops, using the polytope model for both the profiling and the dynamic selection phases. We show on several benchmark programs that our runtime system selects one of the most efficient version with a very low runtime overhead. This quick and efficient selection leads to speedups compared to the usage of a unique version in every execution context.

INTRODUCTION

Program performance is very difficult to ensure in general with the ever extending and changing variety of processor architectures and execution contexts that an application can meet. This is particularly true with the proliferation of multicore architectures, where parallel execution of applications introduces an even higher level of performance uncertainty. Hence static compilers are unable to decide which optimizing code transformations have to be applied to take advantage of the underlying platform resources and reach good performance [START_REF] Pouchet | Iterative optimization in the polyhedral model: Part II, multidimensional time[END_REF]. Even worse, there is no unique optimized version providing the best performance in all situations [START_REF] Ansel | Petabricks: a language and compiler for algorithmic choice[END_REF][START_REF] Mars | Scenario based optimization: A framework for statically enabling online optimizations[END_REF]. Hence the best performance relies on multiple versions of codes having the same functional feature.

For instance, we observed on gemm, a matrix multiplication code using simple control and data structures, that depending on the input data, distinct parallel versions provide the best performance, while running on a single machine. Further, we observed that the version providing the best performance is also varying depending on the computer, even with the same input data. As a consequence, performing adaptive version selection is particularly important in a library including such a function, and even more important if the library is distributed as a binary on several platforms.

From several optimizing and parallelizing code transformation alternatives, a usual approach is to do iterative compilation: compiling and profiling the execution of different versions in order to select the most efficient one [START_REF] Pouchet | Iterative optimization in the polyhedral model: Part II, multidimensional time[END_REF][START_REF] Pouchet | Iterative optimization in the polyhedral model: Part I, one-dimensional time[END_REF]. However, results obtained from such a strategy are closely related to the execution context when profiling on the compilation platform. The execution context includes the target architecture, the input data size and shape, and the processor load.

Another solution is dynamic work distribution, as the dynamic schedule implemented in OpenMP [START_REF]The OpenMP API specification for parallel programming[END_REF] or work stealing methods such as Cilk [START_REF] Blumofe | Scheduling multithreaded computations by work stealing[END_REF]. Those methods are adapted to task parallelism, and when applied to data parallelism they often lead to a high and unpredictable overhead: the task granularity is much smaller in this latter case, increasing the relative cost of work distribution control.

Other studies [START_REF] Ansel | Petabricks: a language and compiler for algorithmic choice[END_REF][START_REF] Mars | Scenario based optimization: A framework for statically enabling online optimizations[END_REF][START_REF] Thomas | A framework for adaptive algorithm selection in STAPL[END_REF][START_REF] Voss | High-level adaptive program optimization with ADAPT[END_REF]] have proposed adaptive runtime selection between several algorithms, code extracts, or versions of a function. Those tools can handle any kind of codes however their granularity of adaptiveness is limited to many complete executions of the versions. We will discuss them further in the related work section.

The contributions of this work are: (1) a static/dynamic collaborative framework providing a fast selection process between nested loops characterized by distinct parallel versions, where the selection process is launched each time a loop nest is executed, thus providing very high adaptiveness; (2) the construction of a ranking table for each considered version and target multicore processor, obtained through an install time profiling run; (3) an implementation of the framework and an evaluation of its effectiveness with several benchmark programs and execution contexts.

We use the polytope model [START_REF] Bastoul | Code generation in the polyhedral model is easier than you think[END_REF] to generate multiple versions of a loop nest. This model provides many loop transformation techniques, implemented, for instance, in the automatic loop parallelization tools PLuTo [START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF] and LeT-SeE/PoCC [START_REF] Pouchet | Iterative optimization in the polyhedral model: Part II, multidimensional time[END_REF][START_REF] Pouchet | Iterative optimization in the polyhedral model: Part I, one-dimensional time[END_REF]. The main differences between the versions are the parallel schedules, the tile sizes and their dimension, and the loop unroll factors.

Typically, loops that fit the polytope model can be found in intensive scientific applications: exit conditions and array access functions are affine functions of iterators and parameters, and there is no test on data values. The polytope model provides an accurate frame allowing us to envisage automatic generation of the parallel code versions and to extract high level information from the source code.

Some previous related proposals have chosen to characterize more general applications using empirical measurements such as the observed behavior of previous executions [START_REF] Mars | Scenario based optimization: A framework for statically enabling online optimizations[END_REF] or machine learning techniques [START_REF] Tian | An inputcentric paradigm for program dynamic optimizations[END_REF]. They have to assume that the empirical measurements used are representative of all the execution contexts in the general case. In the polytope model, the behavior of the loop nests is fully deterministic and statically analyzable. We use this property to select very efficiently a code version without making such strong assumptions on the similarity between the current execution context and the previous ones.

Using a collection of benchmarks (from SPECOMP [START_REF] Aslot | SPEComp: A new benchmark suite for measuring parallel computer performance[END_REF], BLAS [START_REF] Blackford | An updated set of basic linear algebra subprograms (BLAS)[END_REF], and PLuTo [START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF]), we show that our runtime system almost always selects the most efficient out of three to seven distinct versions, and selects another good performing version in the other cases. The time overhead is in order of a millisecond on our test platforms, thus being quite negligible in most cases and allowing speedups compared to the best static version.

The rest of the paper is organized as follows. In the next section, some background on the polytope model is recalled. Section 3 describes our framework in detail. Experiments are presented in Section 4, and related work in Section 5. Finally, Section 6 gives our conclusions and some perspectives.

POLYTOPE MODEL

We consider loop-intensive codes that can be handled by the polytope model: all loop bounds are affine expressions of the outer loop iterators and parameters. Thus, we exclude for example while loops, or for loops with additional control in the body of the loops such as break or goto statements.

State of the art

The polytope model allows us to represent a nest of d loops as a d-dimensional polytope called an iteration domain. Each loop bound is an affine inequality on the variables, that defines a geometric half space. All of these half spaces, corresponding to all loop bounds, intersect as a convex polytope.

Each statement of a loop-intensive code is associated to an iteration domain, computed from the loops that surround this statement. If the loop bounds depend linearly on integer parameters, then the loop nest is parametric and the corresponding polytope is also parametric (which is often the case in intensive computation loop codes, taking as input an M × N matrix for example).

Using the polytope model, one can transform the iteration domains in order to generate a new scanning loop nest, changing the schedule of the statements [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part 1 : one dimensional time[END_REF][START_REF] Feautrier | Some efficient solutions to the affine scheduling problem, part 2 : multidimensional time[END_REF] and leading in the end to a new, possibly parallel, loop nest. Those transformations are expressed by transformation matrices.

In order to be valid, that is to say to respect the program semantics, the transformations must respect the dependencies of the program. The dependence analysis computes which iteration depends on which other. A set of affine constraints automatically built from the access functions of the data allows us to compute the dependence vectors [START_REF] Feautrier | Dataflow analysis of scalar and array references[END_REF] and to ensure the correctness of the transformations.

Tiling can be done by increasing the dimension of the iteration domains by the number of outer iterators to scan tiles, and are easily coupled to the original iterators with additional constraints [START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF]. Unrolling can also be done, after polyhedral transformation.

We use CLooG [START_REF] Bastoul | Code generation in the polyhedral model is easier than you think[END_REF] to generate the resulting code, from the statements iteration domains and the transformation matrices.

Generating different versions

Many performance issues result from the choice of the transformation matrices: degree of parallelism, load balancing, communications volume, cache locality, and vectorization capability depend on this choice. Our framework is based on the capability of generating many different valid versions in order to choose the most performing one at runtime.

For our benchmarks, we generated different schedules by hand. In order to automatically generate many versions, current automatic parallelizers such as PLuTo [START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF] and LeT-SeE/PoCC [START_REF] Pouchet | Iterative optimization in the polyhedral model: Part II, multidimensional time[END_REF][START_REF] Pouchet | Iterative optimization in the polyhedral model: Part I, one-dimensional time[END_REF] could be slightly modified: instead of searching for the unique best optimization using heuristics, apply all the possibly efficient optimizations. It could also be interesting to generate code versions considered as less efficient due to some hardware hypothesis. For example a code version expected to be inefficient due to its cache misuse could be the best one on a hardware accelerator without any cache (FPGA, GPU, ...).

Ehrhart polynomials

The runtime selector presented in Section 3 uses polytope integer points counting, in order to compute the number of iterations that a statement will perform: it is the number of integer points contained in the iteration domain.

This counting is achieved using Ehrhart polynomials [START_REF] Clauss | Parametric analysis of polyhedral iteration spaces[END_REF]: an Ehrhart polynomial associated to a parametric polytope is a piecewise pseudo-polynomial in the parameters of the iteration domains, that expresses the number of integer points contained in the polytope. It is generated at compile time using the Barvinok library [START_REF] Verdoolaege | Counting integer points in parametric polytopes using Barvinok's rational functions[END_REF], and evaluated at runtime, when the parameters values are known.

SELECTION FRAMEWORK

While other proposals select the best available code version after many executions, our framework aims to select it before any loop nest execution. To achieve this goal, we predict each code version execution time and run the version predicted to be the best one. This prediction is not necessarily exact, but it is accurate enough to rank the different code versions.

To predict those execution times, our framework first profiles each code version in various execution contexts and with different number of active threads. The profiled execution times are gathered in a parametric ranking table, see for instance Table 1. A runtime prediction step, executed before each execution of a loop nest, adapts those results according to the current context in order to predict an execution time for each code version. It is embedded in the application binary and is run through a simple function call replacing all the target loop nest executions. As in every version selection framework, all the versions are embedded in the binary file, leading to increase the application size. This increase is linear in the number of versions but is often negligible relatively to the total application size: in our experiments we measured a few kilo-bytes for a loop nest version.

Figure 1 shows a global overview of the framework. The different phases are represented vertically from the generation of the different code versions at compile time to the execution of the loop nest. Automatic generation of the versions is described in Subsection 2.2. The profiling and selection phases are detailed below.

Initial profiling

Accuracy

The profiling phase, run once on the target computer, typically at install time, is in charge of measuring the execution time of each code version in some representative contexts. However many factors can impact the execution time. Those factors have to be taken into account in order to measure execution times in contexts that can safely be considered as representative. We distinguish two main kinds of factors: static and dynamic factors.

The impact of static factors is constant during the execution and are naturally taken into account when measuring the execution time of any run. They are thus considered by our profiling step.

The second category consists in dynamic factors whose impact is not constant. We define the dynamic extern factors as factors depending on events provoked outside of the considered execution, usually other applications. Such factors are considered as out of the scope of the profiling step.

Dynamic intern factors depend on the considered code execution. Among all the potential factors, we identified two main dynamic intern factors relevant to current architectures that might have a major impact on execution time.

First, the overall memory performance of the application is a dynamic intern factor. To handle it, the domain size is doubled between two consecutive measurements until having less than one percent of difference in the consecutively measured execution times. The profiling result is then the last measurement. It allows us to measure the performance of each version when the data cannot fit in the first cache levels without needing any architecture-dependent information.

We observed that this first factor has a small impact compared to the load balance. When parallelizing a loop nest, a parallel loop iteration is a task mapped to a thread. Thus there is a strong link between the parallel loop bounds and the maximal number of simultaneously active threads. The loop bounds depend on which code transformation has been applied and often on the input data size. Thus, this factor is very variable and can lead to large performance loss.

To characterize the load balance impact, each version execution time is measured with a number of parallel iterations increasing from 1 to the number of processor cores during consecutive measurements. The profiling result is therefore parametrized by the maximum number of simultaneously active threads.

Equity

The other key point is to ensure that the comparison based on the predicted execution times is as fair as possible. It means that, despite their differences, all the code versions have to be profiled in a controlled and comparable execution context.

As explained before, the domain size is increased during consecutive profiling runs. Then different versions can reach different domain sizes depending on their characteristics. Thus the ranking table actually contains execution times per loop nest iteration instead of absolute execution times.

To be efficient, codes parallelized with polyhedral transformations are usually tiled [START_REF] Ramanujam | Tiling multidimensional iteration spaces for multicomputers[END_REF]. Significant performance variations can be observed between codes where tiles are fully executed and where only incomplete tile executions occur. To handle this variation, we have chosen to build the profiling domains exclusively with full tiles.

Building a profiling domain

The profiling code of each version is built such that those accuracy and equity requirements are met. The profiling domain has to be made of full tiles and its shape has to be precisely controlled. Iteration domains in the polytope model are defined by a set of linear constraints. To achieve our objectives, the framework removes all the iteration domain boundary constraints from the set of constraints. Then it bounds the parallel loop with a new parameter par sz and the other dimensions at the first tiling level with a new parameter others sz. The par sz parameter will be instantiated with the successive number of threads to repeat the measurements for all the possible number of active threads. The others sz parameter is increased during the profiling to control the domain size until reaching stable measurements. Only the first tiling level is constrained by this new parameter in order to execute full tiles when there are more than one tiling levels. If the domain is not tiled, all the dimensions are constrained.

Figures 2 and3 illustrate the profiling code generation. In Figure 2, we present a code summing the columns of array A in array S. The code is fully tiled and parallelized. Figure 3 shows the corresponding profiling code. The domain boundaries are eliminated and the domain is now controlled by the two new parameters par sz and others sz. While others sz is doubled until reaching a stable measurement, par sz is incremented from 1 to the number of available processor cores. The result is made of the last execution times per iteration for each considered parallel loop trip count.

Figure 3. Sample profiling code

For each value of both parameters, the execution time is measured using the regular operating system timing function. The expression used to compute the number of iterations (others sz * par sz * 64 * 64) is the Ehrhart polynomial representing the iteration domain size.

The full ranking table, built during the profiling step, is two dimensional. One dimension is made of the different versions while the second dimension represents the number of available processors. For a given version, and a given number of processors P available to the application, the ranking table gives the average execution time per iteration of this version when using P processor cores. A sample ranking table is presented in Table 1.

Runtime selection

The runtime selector is in charge of instantiating the parametric ranking tables resulting from the profiling step, to predict each version execution time.

# of cores version 1 version 2 version

Iteration count measurement

The profiling phase considers most of the code versions performance factors, however the load balance has to be evaluated at runtime: the data size and the number of available processor cores impact the maximum number of active threads. To evaluate this load balance, the runtime selector executes a very simplified copy of each loop nest version called prediction nest. The prediction nest counts how many iterations each thread executes in this code version. The number of iterations executed by each possible quantity of threads is then deduced and used to select the entries in the ranking tables.

To build the prediction nest, we consider a code version loop nest. All the loops and statements enclosed in the parallel loop define the computation achieved in a parallel iteration. The workload of this computation is evaluated as the Ehrhart polynomial counting the number of iterations of the loops enclosed in this parallel loop. Thus, at compile time, we syntactically replace the content of the parallel loop by a computation of this Ehrhart polynomial. The value defined by the polynomial is added to a counter specific to the thread which executes the current parallel iteration. At the end of the prediction nest execution, each counter will then contain the number of iterations that one thread has executed.

Figure 4 presents a sample code and its associated prediction nest. Observe that the innermost loop and the statement are replaced by a counting code. This counting code increments the thread-specific counter with the trip count of the removed loop.

It is assumed here that the mapping of the parallel iterations to the threads is the same in the prediction nest and in the corresponding code version, thus excluding dynamic mapping from the scope of our framework. To evaluate the impact of this drawback, we measured the performance of the most efficient version of each benchmark presented in Section 4 using the OpenMP dynamic mapping. We observed that, in average, this dynamic mapping leads to a slowdown compared to static mapping. This slowdown reaches more than 20 % on some architectures. The same experiment has been performed using Cilk [START_REF] Blumofe | Scheduling multithreaded computations by work stealing[END_REF]: the parallel loop is transformed into the cilk for control structure. We observed slowdowns of 45 % in average, reaching more than 100 % in some cases. This illustrates that those dynamic scheduling or work-stealing systems can suffer from high overheads and are not suited to deal with regular loop nests.

To deduce the number of parallel iterations executed by each thread quantity, the counters array cnt is sorted in descending order. For each position 1 ≤ i < nb threads, we can state that cnt[i-1] -cnt[i] iterations have been executed in parallel by exactly i threads. For example, consider a 4-cores computer where the following numbers of iterations have been measured by the prediction nest: We can deduce that one thread has executed (1000 -800) iterations while the others were idle, two threads have executed (800 -200) iterations in parallel and four threads have simultaneously executed 200 iterations. There are never exactly three active threads in this case.

Predicting the execution time

Once the iteration counts are known, the execution time can easily be deduced using the parametric ranking table. For example, consider the previous example where 200 iterations are executed by one thread, 600 by two threads and 200 by four threads. For this example, we consider the profiled execution times of the first version in Table 1. The execution time is then computed as 200 × 4 × 5 ms for the part of the iteration domain executed by four threads, 600 × 2 × 10 ms for the iterations executed by two threads plus 200 × 1 × 30 ms for the sequential iterations. The sum predicts an execution time of 22 seconds.

The prediction time computation is executed with every version of the loop nest in the current execution context, before running the best one.

EXPERIMENTS

We run our experiments on three computers with different multicore processors: an AMD Opteron 2431 processor with six cores , an AMD Phenom II 965 processor with four cores and an Intel Core i7 920 processor with four cores and hyperthreading activated. The experiments are run on Linux 2.6.35 systems.

The benchmark codes are 12 common polyhedral loop nests. The code 2mm is made of two matrix multiply (D = A × B ×C), adi is the ADI kernel provided for example with PLuTo, covariance is a covariance matrix computation, gemm and gemver are taken from BLAS [START_REF] Blackford | An updated set of basic linear algebra subprograms (BLAS)[END_REF], jacobi-1D and jacobi-2D are the 1D and 2D versions of the Jacobi kernel, lu is a LU decomposition kernel, matmul is a simple matrix multiply, matmul-init is a matrix multiply combined with the initialization of the result matrix, mgrid is a kernel extracted from the mgrid code in SPECOMP [START_REF] Aslot | SPEComp: A new benchmark suite for measuring parallel computer performance[END_REF] and seidel is a Gauss-Seidel kernel as provided with PLuTo.

Those kernels are typically put in libraries and called many times by applications cumulating the benefits of our system.

For each code, many versions are generated. One of those version is automatically generated by the PLuTo parallelizing compiler. Other versions have been designed by an expert, changing the number of tiling levels and the tile sizes, and making polyhedral loop transformations. As explained in Subsection 2.2, those versions could have been generated by automatic tools. The details about those versions are available on request. From each version, the profiling and prediction codes are generated by a set of fully automatic scripts in our implementation. The result of the prediction codes determines the framework choice in each execution context.

The runtime system is evaluated on each of the three computers, for each code with five different problem sizes leading to execution times near to a couple of minutes. Each measurement is repeated with a number of threads varying from one to the number of logical cores simulating different resource availability. The presented statistics enclose more than one thousand runs, each one being the median value out of five executions.

Performance variation

During our experiments, we have observed many performance variations. First, there are performance variations for a given code version among different computers for a fixed problem size and number of cores. For example with adi, when considering a problem size of 2448 and 5 available processor cores, version 2 is the best one with the Opteron processor but the worst one with the Core i7 processor.

We have also observed performance variations when the number of available cores on a given computer varies for a fixed problem size. For example, with the Core i7 processor for a problem size of 9996, the first version of gemver is the best one when only one processor core is available but the worst one when all the cores can be used.

Finally, when the problem size changes there are performance variations for a fixed computer and a fixed number of available cores. For example, with all the processor cores of the Phenom processor, the fourth version of mgrid is inefficient for a problem size of 400 but is the best one for a problem size of 496.

The best version is therefore not necessarily the same in all the execution contexts and our runtime system succeeds in selecting the best one in those presented cases.

We have observed that the transformations making a version perform better than another are the tile size, the number of tiling levels and the parallel schedules. All those transformations interact and have a complex impact on performance, changing with the execution context. This emphases the need of a dynamic system, such as our framework, able to consider all of those transformations together.

Evaluation results

We present in Table 2 a quantitative study of our runtime system efficiency. All the presented speedups are averages of all the execution contexts where the number of active processors and the data size vary. They take the dynamic system overhead into account. We measured this overhead, induced when predicting the execution times, to be less than a tenth of milliseconds per code version. We define the best static version as the version that performs best in average in all the considered execution contexts on an architecture. Such version could be the one statically chosen by a perfect offline system. If a best static version would be the fastest one in every execution context, no speedup can be expected from any runtime system.

First we show the speedup of our framework compared to the best static version. We can see that our system sometimes leads to a slowdown compared to this version. Those slowdowns are mostly due to incorrect choices in some execution contexts. Such incorrect choices are made when the approximations of the profiling step lead to mispredict the execution time of some versions in some of the execution contexts.

Two typical examples are gemver and jacobi-2d on the Core i7 platform. We can see however that our framework usually selects another well performing version in those cases as the slowdown remains very limited. It is important to note that this slowdown is related to a perfect static system which is extremely complex to build as it would require to statically select the best parallel version for each architecture. Such a system does not exist currently.

Our framework reaches some speedups in some cases (e.g. adi or gemm), meaning that the execution time is, in average, shorter than the execution time of the version that would have been selected by a perfect static system. This speedup is impossible to reach with only one version, illustrating the importance of using a dynamic selection system.

In the second column, we show the speedup of the best static version compared to the best version in each execution context. No speedup higher than 100 % can be expected here: it is impossible to have any version faster than the fastest version in a specific context. A speedup of 100 % here means that a unique version is the best one in all the tested execution contexts. Those speedups illustrate that even a perfect static system cannot reach the maximum performance in all the cases.

Considering those statistics, we can say that our framework is able to select a very well performing version in all the presented cases. In some cases, it is even able to outperform any hypothetical perfect offline system.

RELATED WORK

Iterative compilation and dynamic work distribution (work-stealing or dynamic scheduling) have been discussed in the introduction and in Subsection 3.2.1. This section focuses on adaptive runtime selection methods.

In the ATLAS project [START_REF] Whaley | Automated empirical optimizations of software and the ATLAS project[END_REF], empirical timings are used in order to choose the best method for a given architecture "in a matter of hours". Further it does not handle parallel programs.

With the ADAPT system [START_REF] Voss | High-level adaptive program optimization with ADAPT[END_REF], a specific language allows the user to describe optimizations and heuristics for applying these optimizations dynamically. However, the resulting optimizer is run on a free processor or on a remote machine across the network.

PetaBricks [START_REF] Ansel | Petabricks: a language and compiler for algorithmic choice[END_REF] provides a language and a compiler where having multiple implementations of multiple algorithms to solve a problem is the natural way of programming. The as-sociated runtime system uses a choice dependency graph to select one or another algorithm and implementation at different steps of the whole computation, thus resulting to an optimized hybrid algorithm. Such an approach is suitable for programs where it is obviously possible to switch from one algorithm to another while still making progress in the whole computation. Further, such a changing behavior must not induce overheads due to a compulsory cache flush for instance.

Mars and Hundt's static/dynamic SBO framework [START_REF] Mars | Scenario based optimization: A framework for statically enabling online optimizations[END_REF] consists in generating at compile-time several versions of a function that are related to different dynamic scenarios. These scenarios are identified at runtime thanks to the microprocessor event registers. Execution is dynamically rerouted to the code relevant to the current identified scenario. However, execution is not rerouted during a function call, but for the next calls. Further, it seems difficult to use this approach with parallel program since it is actually challenging with multicore processors to deduce accurate global multithreaded program behaviors from registers disseminated on the cores.

The STAPL adaptive selection framework [START_REF] Thomas | A framework for adaptive algorithm selection in STAPL[END_REF] runs a profiling execution at install time to extract architectural dependent information. This information is used at runtime combined with previous runs and training runs performance measurements through machine learning to select the best version. Their system requires many training runs before being able to take good decisions.

More recently Tian et al. [START_REF] Tian | An inputcentric paradigm for program dynamic optimizations[END_REF] propose an input-centric program behavior analysis for general programs, a statistical approach where program inputs have to be characterized differently depending on the target application -input size, data distribution, etc. -, and program behavior is represented by relations between some programs parameters, as for instance loop trip-counts. This modeling is used to achieve some dynamic version selection. Such relevant statistical relations seem difficult to be determined for any kind of programs. The approach necessarily suffers from approximations, hardware characteristics are not handled and overall, this work does not consider parallel and multicore programs.

CONCLUSION AND PERSPECTIVES

In this paper, we present a new framework for adaptive code selection of parallel loop nests at runtime. It handles varying input data and varying execution contexts, such as the number of cores or the cache architecture. We showed on a set of benchmarks that our framework is efficient: depending on the input data shape, on the target architecture and cores availability, the version performing best is not the same, and our method mostly selects the best performing one with a low runtime overhead. Further, it provides speedups that could not have been reached with one statically selected version when considering several calls with different input data shapes.

We are able to consider very different versions where a wide range of optimizations can be applied, from common compiling optimizations to complex polyhedral transformations. Further, our initial profiling phase is also very fast, lasting less than a couple of minutes in our experiments.

As a perspective, we plan to improve our parametric profiling phase by enriching the structure of the simplified nest in order to better consider more complex iteration domains. Another interesting extension would be to consider loop nests that are not strictly polyhedral. Finally, a complete dynamic system, i.e. without any initial profiling phase, could be built either by alleviating some mechanisms or by using actual runs to build the parametric ranking table.

Figure 1 .

 1 Figure 1. Framework overview

Figure 4 .

 4 Figure 4. Sample loop nest (top) and its corresponding prediction nest (bottom)

Table 1 .

 1 Sample parametric ranking table built on a 4-cores processor for 3 versions. The table content is made of execution times per iteration.

	3

Table 2 .

 2 Speedup of our framework compared to the best static version and speedup of the best static version compared to the best version in each context.

	Processor	Program	Framework Max static speedup speedup
		2mm	100.0 %	100.0 %
		adi	101.5 %	97.5 %
		covariance	100.0 %	99.9 %
		gemm	102.8 %	96.9 %
		gemver	100.0 %	99.8 %
	Phenom	jacobi-1d jacobi-2d	99.7 % 99.5 %	100.0 % 100.0 %
		lu	100.0 %	100.0 %
		matmul	100.0 %	100.0 %
		matmul-init	100.0 %	100.0 %
		mgrid	98.0 %	98.1 %
		seidel	100.1 %	99.0 %
		2mm	100.0 %	100.0 %
		adi	102.2 %	97.3 %
		covariance	100.0 %	99.8 %
		gemm	101.4 %	96.7 %
		gemver	99.9 %	99.8 %
	Opteron	jacobi-1d jacobi-2d	99.6 % 100.0 %	100.0 % 100.0 %
		lu	100.0 %	100.0 %
		matmul	100.0 %	100.0 %
		matmul-init	100.0 %	100.0 %
		mgrid	97.8 %	98.5 %
		seidel	100.6 %	98.3 %
		2mm	100.0 %	100.0 %
		adi	102.4 %	97.5 %
		covariance	96.9 %	99.7 %
		gemm	100.0 %	93.5 %
		gemver	89.5 %	91.6 %
	Core i7	jacobi-1d jacobi-2d	99.6 % 90.6 %	99.9 % 99.6 %
		lu	93.1 %	98.3 %
		matmul	100.0 %	98.5 %
		matmul-init	100.0 %	100.0 %
		mgrid	100.5 %	97.0 %
		seidel	99.9 %	99.6 %


Anis Rahman
email: anis.rahman@gipsa-lab.grenoble-inp.fr

Dominique Houzet
email: dominique.houzet@gipsa-lab.grenoble-inp.fr

Denis Pellerin
email: denis.pellerin@gipsa-lab.grenoble-inp.fr

Lionel Agud
email: lionel.agud@gipsa-lab.grenoble-inp.fr

GPU implementation of motion estimation for visual saliency

Visual attention is a complex concept that includes many processes to find the region of concentration in a visual scene. In this paper, we discuss a spatio-temporal visual saliency model where the visual information contained in videos is divided into two types: static and dynamic that are processed by two separate pathways. These pathways produce intermediate saliency maps that are merged together to get salient regions distinct from what surround them. Evidently, to realize a more robust model will involve inclusion of more complex processes. Likewise, the dynamic pathway of the model involves compute-intensive motion estimation,that when implemented on GPU resulted in a speedup of up to 40x against its sequential counterpart. The implementation involves a number of code and memory optimizations to get the performance gains, resultantly materializing real-time video analysis capability for the visual saliency model.

Introduction

The motivation behind the mimicking human vision system is to develop a complete computer vision system. The research in this field includes methods to extract useful information for proper characterization of a visual scene, thus making possible its analysis for a computer system. One of these many techniques are by estimating the motion to find the regions of attention. Motion estimation is a process to examine movement of an object by generating motion vectors for a sequence of images. The motion estimator implemented in this paper is proposed by Bruno et al. [START_REF] Bruno | Robust motion estimation using spatial gabor-like filters[END_REF], which afterwards is incorporated into the spatio-temporal visual saliency model developed by Marat et al. [START_REF] Marat | Modelling spatio-temporal saliency to predict gaze direction for short videos[END_REF]. The model proposed focuses on the information extracted from the visual field, without any top-down modalities, and in the end resulting in a region of attention. The saliency model is inspired from Feature Integration Theory [START_REF] Treisman | A feature-integration theory of attention[END_REF] that involves the decomposition of the visual field into various simple attributes. These are afterwards treated using several processes, and at the end are combined to find the regions that dominate others. The complete model is compute-intensive, hence it is a perfect candidate for parallelization to be used for real-time video analysis. GPU presents us with an inherently parallel architecture for image-based algorithms that are often complex and time-consuming. The desired performance can be extracted off these highly parallel GPUs by employing various methods to optimize the memory usage and thread allocation/synchronization. In this paper, we implement the dynamic pathway of the visual saliency model on the GPU. We concentrate on the optimization of the most computeintensive kernel, the motion estimator, to achieve the desired speedup. The paper is organized as follows: In section 2, we presented a short introduction to the two pathway spatio-temporal visual saliency model. In section 3, we summarize the dynamic pathway of the model with a brief description of the motion estimator used. Afterwards, the section summarizes our implemented algorithm with the different optimizations made. In section 4, we report the achieved speedups, profiling, and validity of these results is evaluated. In the end, a conclusion of the article and its future prospects are discussed.

Spatio-temporal visual saliency model

The bottom-up model [START_REF] Marat | Modelling spatio-temporal saliency to predict gaze direction for short videos[END_REF] illustrated in figure 1 is linearly modeled from the retina to the visual cortical cells. Moreover, the separation of useful information into two distinct signals makes the processing of the information relatively easier. Furthermore, the partial maps from the two pathways are fused together into a single saliency map. All these features contribute in developing a model based on the human visual system.

Static pathway

The static pathway is based on retinal filtering, which then is followed by a bank of Gabor filters.

The main modalities used here are frequencies and orientations. The retinal model imitates the horizontal, bipolar, and ganglion cells to expose more detail by increasing the luminance of higher frequencies in the visual input. The primary visual cortex is a model of simple cell receptive fields sensitive to visual signal orientations and spatial frequencies. This is modeled as a 2D Gabor filter bank that processes the visual information in different frequencies and orientations, resulting in 24 partial maps. These filters demonstrate optimal localization properties and good compromise of resolution between frequency and spatial domains. In primate visual system, the response of cells is dependent on their neuronal environment; its lateral connections. Thus, this activity is modeled as a linear combination of simple cells interacting with their neighbors. Afterwards, intermediate energy maps are strengthened by normalization and summed up to extract a saliency map for the static pathway.

Dynamic pathway On the other hand, the dynamic pathway finds a saliency map from a moving scene. As a preprocessing, this pathway uses camera motion compensation [START_REF] Odobez | Robust multiresolution estimation of parametric motion models applied to complex scenes[END_REF] for estimation of dominant motion of the salient regions according to their background. This estimation is followed by a 2D motion estimation [START_REF] Bruno | Robust motion estimation using spatial gabor-like filters[END_REF] to find local motion with respect to the background. The motion vectors are calculated with modalities of speed, orientation and direction, which are then treated with a temporal median filtering to remove noise and to produce a saliency map. Ultimately, the resulting maps from both pathways are fused together into a master map with salient regions; one's with highest energy.

Implementation of the dynamic pathway

CUDA [START_REF]NVIDIA CUDA Compute Unified Device Architecture -Pro-gramming Guide[END_REF] and OpenCL [START_REF] Group | Opencl -the open standard for parallel programming of heterogeneous systems[END_REF] threading model virtually launches a sea of light-weight threads onto the streaming multiprocessors on the graphics hardware. All in all, this provides fine-grain data parallelism across threads and coarse-grain data and task parallelism across thread blocks. Therefore, almost any algorithm exhibiting data-parallelism is suited to be ported on to graphics processors, and exploiting its raw computational power. The first step is to develop a sequential version, which can be used to select and code the compute-intensive functional units as kernels that are to be executed on graphics processors. This task is fairly easier as this version acts as a wellstructured reference. The main functions of the dynamic pathway of the visual saliency model include: applying the retinal filtering, convolving with the Gabor filter bank to decompose the rescaled images into sub-bands with different orientations, solving the over-determined system of equations, applying The algorithm 1 presents the various steps of the dynamic pathway, where the first step is camera motion compensation done on the host. This compensation is then followed by other functions implemented on GPU like retina filter, Gabor filter, motion estimator, Gaussian prefiltering, and many others. At the end of this pathway, we get the dynamic saliency map.

Motion estimator

The motion estimator [START_REF] Bruno | Robust motion estimation using spatial gabor-like filters[END_REF] presented here employs a differential method using Gabor filters to estimate local motion. This estimation is done by solving a robust system of equations of optical flow. The method employed is to find the pixels that conserve their luminance for some interval of time. The motion vector υ(p) = (υ x , υ y ) can be found using the equation of optical flow:

dI (p, t) dt = I (p, t) • υ (p) + ∂I (p, t) ∂t = 0
where I(x, y, t) represents the spatial gradients of luminance I(x, y, t). We can use the direction of the velocity component from these gradients to find the direction of the motion of an object.

The motion vectors are calculated by averaging the movement to its spatial neighborhood. This spatial continuity within the optical flow is achieved by the convolving the spatio-temporal image sequence with a Gabor filter bank:

υ x • ∂ (I * G i ) ∂ x + υ y • ∂ (I * G i ) ∂ y + ∂ (I * G i ) ∂ t = 0
the 2D filter bank comprises of N filters Gi with the same radial frequencies. Hence, resulting in a system of N equations for each pixel:

    Ω x 2 Ω y 1 Ω x 2 Ω y 2 ... Ω x n Ω y n     • υ x υ y =     Ω t 2 Ω t 2 ... Ω t n    
this over-determined system is solved using least squares method(Biweight Tuckey test) [START_REF] Bruno | Robust motion estimation using spatial gabor-like filters[END_REF].

The approximation begins with a sub-sampled image from the highest level of the pyramid. This multi-pattern approach allows robust estimation of motion for every pixel.

Input: Pyramid of sub-sampled images Output:

Velocity vectors V 1 foreach Level k of pyramid do 2 if k == K then 3
Gabor filtering I(x,y,t)

Different optimizations

Reducing global accesses After convolving with Gabor filters and applying the spatial and temporal gradients, we get an over-determined system of N solutions for each pixel of the image. This system is resolved using iterative weighted least-square estimations. This method involves the calculation of N (2N -1) intermediate values that are stored as local array variables. Furthermore, these values are accessed at different multiple times, which leads to higher global memory accesses because these arrays are placed on the device memory. To reduce these accesses, we use shared memory. However, the amount of shared memory is limited that makes it impossible to eliminate all the local array variables. As a result, this improvement is not visible because shared memory allocated affects the number of active thread blocks residing on a single multiprocessor.

Reducing register count

The estimator kernel is involving several steps that are complex, and causes the register count on the boundary. The higher register count results in a limited number of active thread blocks per streaming multiprocessor, thus the occupancy is low. In our first implementation, the register count is 22 that is reduced to 15 with the use of shared memory to store some of the local variables. As a result, the number of active thread blocks increases to 8 per streaming multiprocessor, and the device occupancy increases from 0.33 to 0.67. The use of constant memory would be a preferable solution, in case if each variable is accessed by all the threads simultaneously. Here, they happen to be used and modified in several instructions, thus the candidate variables to be placed in the shared region are carefully selected to minimize the need for synchronization barriers.

Using texture memory The estimation kernel comprises of two operations: one to find N(N-1) equations for the pixel, and then averaging them to get components of its motion vector. To reduce the complexity of the kernel, we cut and place these two operations into separate device kernels. This result is reduced complexity of the kernel for further optimizations and tweaks. On the other hand, it results in calculating and storing of the intermediate equations from the application of spatial and temporal gradients on to the device memory. These equations are used by the second kernel to approximate the motion. Here, the use of texture memory provides faster accesses to these prefetched values, rather than using local arrays that resulted in slower device memory accesses. Therefore, it helps to avoid global memory latency by taking advantage of caching mechanisms of the texture memory, hence, it results up to 10% performance improvement.

Multi-GPU solution

In such a multi-threaded model, there is complexity involved in selecting an efficient strategy for creation and destruction of threads, resource utilization, and load balancing. Furthermore, inter-GPU communication between GPUs might be an overhead affecting the overall performance. This overhead can be tackled by overlapping communications with computations using streaming available in the CUDA programming model. In our implementation of the saliency model, we find that the static pathway is faster than the dynamic pathway of the model. To test for more improvement, we employ a three GPU shared system, where the first GPU executes the static pathway while the dynamic pathway is divided into two halves with almost the same execution times as shown in figure 2. We launch the threads executing their portion of the calculation of the saliency model on their assigned GPU. The threads are initially in sleep mode, and they wait for a signal from the main thread that their inputs are available to be processed. The threads copy their inputs onto their assigned GPU's device memory, and then invoke the device kernels. After the completion of the kernels, the results are copied back to the host, where the next frame is waiting to be fetched for processing. The estimator kernel takes about 120ms, hence we cut the kernel into two almost equal portions of 60ms as shown in figure 2. The first kernel applies the retina and Gabor filtering, while the other performs the approximation of local motion using the resulting system of equations from the previous kernel. Evidently, there is an expensive device to device transfer within thread 1 and thread 2, which involves the data to be copied on to host memory. Afterwards, thread 2 is responsible for the estimation. In our multi-GPU implementation, we have three threads for the three available graphics devices each responsible for their portion of the visual saliency model as following: thread 0 calculates the static saliency map, thread 1 applies the two filters, while thread 2 uses the outputted system of equations from thread 1 to estimate the motion vectors. After the temporal filtering, thread two results in dynamic saliency map that is copied back to the host memory. Finally, the main thread fuses the resulting saliency maps from thread 0 and thread 2 into the final visual saliency map.

Results

All implementations are tested on a 2.67 GHz quad-core system with 10GB of main memory, and Windows 7 running on it. The parallel version is implemented using lat-est CUDA v3.0, and is evaluated on a NVIDIA GeForce GTX 285 graphics card. The static and dynamic pathways are evaluated with various image sizes, including standard 640×480 pixels, whereas the result of the estimator is evaluated using several datasets of image sequences with sizes ranging from 150×150 to 316×252 pixels.

Speedup

The CUDA implementation is tested against on the newest GeForce GTX 285. The device has 30 streaming multiprocessor with total 240 cores of clock rates 1.48 GHz each, providing 1.062 TFLOPS of single-precision and 89 GFLOPS of double precision computational power with memory bandwidth 159 GB/sec.The figure 3 shows the gains achieved for the static and dynamic pathway with different image sizes against their sequential versions. Figure 4 shows the timings for the static and dynamic pathway on a GTX 285 graphics card using various image sizes. The timing for both pathways starts up at about the same level but as the image size gets bigger the computation time for dynamic pathway increases. The reason is due to the increased complexity of the estimator algorithm involving more computations and device memory accesses.

Evaluation of results

To evaluate the correctness of the motion estimator, we calculate error between estimated and real optical flows us-

α e = arccos uu r + vv r + 1 √ u 2 + v 2 + 1 u 2 r + v 2 r + 1
where a e is the angular error for a given pixel with (u, v) the estimated and (u r , v r ) the real motion vectors. We used "treetran" and "treediv" image sequences for the evaluation, showing translational and divergent motion respectively [START_REF] Bruno | Robust motion estimation using spatial gabor-like filters[END_REF].

The results illustrated in table 3 are obtained using "treetran" and "treediv" image sequences of sizes 150 × 150 pixels.

Real-time solution

OpenCV [START_REF] Bradski | Learning OpenCV: Computer Vision with the OpenCV Library[END_REF] and OpenVIDIA [START_REF] Fung | Openvidia: Parallel gpu computer vision[END_REF] are open source platforms providing functions for interfacing to video input and output devices. It also facilitates the programming on GPU through a collection of utility functions, and various implementations of image processing and computer vision algorithms. These algorithms range to include feature detection and tracking, skin tone tracking, projective panoramas, and a lot more. After the parallel implementation of the spatio-temporal visual saliency model, we used OpenCV for interfacing with webcam and videos to demonstrate the real-time processing of the model. This demonstration resulted in execution of the model at 25 fps on GeForce 285GTX. Therefore, the performance gains achieved shows that GPUs can be a strong candidate for parallelization of computer vision algorithms.

Conclusion

The GPU implementation for dynamic pathway of the visual saliency model comprises of several compute-intensive kernels like recursive Gaussian filters, the motion estimator, and many others. After various optimizations and tweaks in the code, the final solution resulted in a speedup of up to 40 times against its sequential implementation. Hence, it is now possible to achieve real-time video analysis using the model. Furthermore, the acceleration allows us to add more attributes to enhance the quality of information extracted from a visual scene like the integration of depth proposed in [START_REF] Ouerhani | Visual attention: from bio-inspired modeling to real-time implementation[END_REF]. Moreover, it will be interesting to investigate the inclusion of top-down processes as proposed in [START_REF] Peters | Beyond bottom-up : Incorporating task-dependent influences into a computational model of spatial attention[END_REF]. Finally, the performance gains on GPU will enable our model to be used for various applications such as scene recognition for mobile robots [START_REF] Siagian | Rapid biologically-inspired scene classification using features shared with visual attention[END_REF], video compression [START_REF] Larson | Can visual fixation patterns improve image fidelity assessment[END_REF][START_REF] Ninassi | Does where you gaze on an image affect your perception of quality? applying visual attention to image quality[END_REF], computer graphics rendering [START_REF] Yee | In neurobiology of attention, chapter Attention for computer graphics rendering[END_REF][START_REF] Yee | Spatiotemporal sensitivity and visual attention for efficient rendering of dynamic environments[END_REF]. The idea, here is to determine a region with better saliency requiring good spatial resolution.

Figure 1 .

 1 Figure 1. Block diagram of spatio-temporal visual saliency model

49 12 ( 2 :

 122 Gabor filtering I(x,y,t-1)Gabor filtering I(x,y,t) 10 Gabor filtering I(x,y,t -1)11 end G x , G y , G t )← Calculate gradients 13Vector v ←Resolution I(x,y,t), I(x,y,t-1), G 14 if k == K then 15 Motion estimator algorithm

Figure 2 .

 2 Figure 2. Block diagram for partitioned dynamic pathway

Figure 3 .

 3 Figure 3. Speedups for two pathways versus sequential C implementation

Figure 4 .

 4 Figure 4. Timings for two pathways on GTX 285 for various image sizes

Table 1 . Profile of motion estimator kernel

 1

		No. of	Occupancy GPU
		registers		time
	Naive solution	18	0.75	32
	Shared memory	29	0.5	31.5
	Shared + prefetching	32	0.5	31.7
	Shared + textures	17	0.375	32

Table 2 illustrates

 2 GPU time for different steps of the dynamic pathways of the visual saliency model on a GTX 285 graphics card with test image sequences of size 640×480 pixels.

Table 2 . Computational cost of each step in dynamic pathway

 2

	GeForce GTX	GPU

Table 3 . Evaluating the M-estimator

 3

		Angular error
		treetran	treediv
		x	σ	x	σ
	Matlab 1.63 5.27 6.06 8.22
	C	1.10 0.99 4.15 2.69
	CUDA 1.19 1.00 5.73 3.91


Christine Opagiste

Maxime Leroux

Pierre Rodiere

G Garbarino

Sébastien Pairis

Pierre Bordet

Pascal Lejay

P Rodière

Lu 5 Ir 4 Si 10 whiskers: morphology, crystal structure, superconducting and charge density wave transition studies

Keywords: Morphology, crystal structure, superconducting and charge density wave transition studies A1. Crystal structure, A1. X-ray diffraction, A2. Growth from melt, A2. Single crystal growth, B1.Rare earth compounds, B2. superconducting materials PACS classification codes: 61.50.Ah, 61.05.C-, 81.10.Fq, 68.70.+w, 72.15.Eb, 74.25.F-

Introduction

The compounds of the intermetallic series R 5 Ir 4 Si 10 with R = rare-earth exhibit a behaviour where either superconductivity and Charge Density Wave (CDW) ordering or magnetic ordering and CDW ordering coexist in the same material depending on the magnetic or non magnetic nature of the rare earth element R [1 -6]. Among them, the Lu 5 Ir 4 Si 10 compound becomes superconducting at 4 K and shows a strongly coupled CDW transition around 80K [1 -4]. Like the other compounds of those series, Lu 5 Ir 4 Si 10 crystallizes in the tetragonal Sc 5 Co 4 Si 10 type structure with space group P4/mbm [7]. Single crystals of Lu 5 Ir 4 Si 10 are generally grown in a three-arc furnace using a modified Czochralski technique [3,8] or by Ref : J. Cryst. Growth 312, 3204 (2010) p2/8 the floating zone technique using an image furnace in Ar atmosphere [9]. In our process to obtain bulk single crystal grown with the modified Czochralski technique, we have observed a spontaneous growth of whiskers out of the under cooled melt. In this paper we describe the morphology, the crystallographic structure refinement and the transport property measurements performed on Lu 5 Ir 4 Si 10 whiskers.

Experimental

A polycrystalline sample of Lu 5 Ir 4 Si 10 was prepared in an induction furnace, in a cold copper crucible and in a high purity argon atmosphere, with a stoichiometric amount of Lu (99.96%, Ames Lab.), Ir (99.9%, Alfa Aesar) and Si (99.9999%, Alfa Aesar). The sample was melted several times to improve its homogeneity. Mass losses during these first steps were less than 0.1%. The polycrystalline sample was then placed in a three-arc furnace in a water cooled copper crucible and in a purified argon atmosphere. The cooling speed seems to be an important parameter for the growth of whiskers. Whiskers appear only when the melt is rapidly cooled in a fraction of a second, but not when the liquid phase is slowly cooled. Under convenient conditions, lots of whiskers can be seen emerging from the polycrystalline solid mass (fig. 1). A typical brittle needle-like single crystal is shown in the inset of figure 1. whiskers can be observed emerging from the polycrystalline sample when this one is thermally quenched. In inset: a typical brittle needle-like whisker of 1.5 mm long.

Results and discussion

Whiskers surface morphology, orientation and mosaicity

From Field Emission Scanning Electron Microscopy observations (FESEM -Zeiss Ultra, all images presented in this paper are secondary electron images) we find that whiskers exhibit a very regular shape of a plate like crystal with typical dimensions: 15 x 60 x 1000 µm 3 (fig. 2). We have never observed whiskers with a scroll shape. Observations of the surface reveal typical growth figures resulting of the formation of dendrites (fig. 3). This is not surprising since we have cooled the sample in a fraction of a second. As seen on figure 4-a a growing dendrite is observed on the surface and above a cavity.

The dendrite has one long branch lying on the surface with side branches in the same plane. We can distinguish the second and third orders of the dendrite branches (fig. 4-b). Typical whisker growth begins with a square pyramid (fig. 5) and has only been observed on the periphery of the polycrystalline thermally quenched sample. On the top surface of the polycrystalline sample we have never observed the formation of whiskers but only the formation of a multitude of square pyramidal concretions (fig 6). The extent to which the whiskers are not pointing out at the top depends on the mass transport and gravity and also on the temperature gradients. It is difficult to explain why the macroscopic structure of the whiskers looks like a plate like structure with a perfect rectangular cross section with sharp edges (fig. 3 The surface formed by the dendrites growth is oriented along the [110] direction. It has been impossible to get such information on whiskers, due to the imperfect state of their surface and to the low thickness probed by EBSD technique. In order to get the whiskers orientation, we have used X-ray diffraction measurements. We have found that the whiskers grow along the c axis direction and that each surface is oriented along the [110] direction (see fig. 2).

In the inset of figure 7 a zoom of a diffraction image of the reflections (110) and (001) is shown, as well as the angle dependence of its intensities in semi-log scale. The measurements have been developed in the ID27 beam line of the ESRF in a transmition configuration using a MARCCD detector with a focused X ray beam of

intensity is given in log scale). The half width half maximum of a Gaussian fit is represented by

Δ in both cases. It can be observed that the mosaicity of these two reflections is better than € 0.15 o . In the inset, a diffraction CCD image of the whisker is shown.

X-ray diffraction refinement

A whisker of dimensions 65 x 30 x 600 µm 3 4.190(2) Å. The reflections were integrated using the EvalCCD suite [13] and anisotropic effects were corrected using the Sadabs program [14] leading to a Rint of 8.47% for all 1650 reflections averaged from 21612 measured ones. This rather high value of Rint is probably due to the high absorption (µ=425.9 cm -1 ) and strongly anisotropic shape of the sample. Analysis of the systematic extinction conditions leads to the P4/mbm space group. The structure was solved by direct methods using the Sir97 program [15]. The refinement was carried out using the Jana2000

program [16]. All atoms were refined with anisotropic atomic displacement parameters and an isotropic extinction correction, using a 1/sig(Fo) 2 weighting scheme, where the sig(Fo) are the estimated standard deviation on observed structure factor moduli. The final agreement factors were Rwobs=3.18% and Rwall=3.39%, gofobs=1.86 and gofall=1.86 (reflections are observed for I>3.sig(I)). The positional and harmonic displacement parameters obtained are displayed in table 1 and 2.

These results are very close to those reported by Becker et al. [3] and obtained for a bulk single crystal.

Electrical resistivity measurements

We carried out measurements of electrical resistivity from 300K down to 2K on a whisker. We used a traditional four-wire technique with an automatic resistance measurement apparatus TRMC2 [START_REF]TRMC2 resistance bridge[END_REF]. The typical contact resistances were several Ohms. Around 77K, the calibration of the thermometer was checked with the boiling temperature of nitrogen at ambient Ref : J. Cryst. Growth 312, 3204 (2010) p7/8

pressure. Moreover, to check the absence of any self-heating, different current densities ranging from 0.02 A/cm 2 to 20 A/cm 2 were used.

The resistivity of the whisker is shown figure 8. This is quantitatively the same temperature dependence as in c-axis bulk samples reported in the literature [3].

Below 76.6K, a clear jump on the resistivity of 30% is due to the formation of the Charge Density Wave. A highly reproducible temperature hysteresis of 1K wide has been measured.

The same width has been measured with different speeds ranging from 50mK/min to 500mK/min. The hysteresis appears over the whole transition jump contrary to previous report on bulk sample [START_REF] Ramakrishnan | [END_REF]. This hysteresis is a clear evidence of the 1 st order nature of the transition previously reported [3]. At lower temperature a full superconducting transition is measured. The critical temperature is 4.1K (50% of the residual resistivity) and 35mK wide (10%-90% criterion). No significant changes were observed with current densities from 20A/cm 2 down to 0.02 A/cm 2 .

Figure 1 :

 1 Figure 1:

  Ref : J. Cryst. Growth 312, 3204 (2010) p3/8

Figure 2 :

 2 Figure 2: SEM image of a whisker. Orientation of the crystal axis is given. The face is oriented along the [110] direction.

figure 3 :

 3 figure 3: (a) Detail of the extremity of a whisker showing a surface with signs of chevronshaped growth and sharp edges (b and c).

  figure 4: (a) Dendrites observed on the surface of the polycrystalline sample. (b) Zoom

figure 7 :

 7 figure 7 (a) and (b)) we obtained that the mosaicity of the (110) and (001) reflections is better than

Figure 8 :

 8 Figure 8: Temperature dependence of the electrical resistivity of a whisker of Lu 5 Ir 4 Si 10 along caxis with a current density of 20A/cm 2 . a) from room temperature to 2K b) zoom on the superconducting transition at 4.1K c) zoom on the temperature hysteresis at the Charge Density Wave transition around 76.6K.

Table 1 :

 1 Atomic coordinates and isotropic displacement parameters for Lu 5 Ir 4 Si 10 .

	Atom	Wyckoff pos.	X	Y	Z	Uiso (Å 2 )
	Ir	8i	0.018812(18)	0.245589(17)	0	0.00530(5)
	Lu1	2a	0	0	0	0.00691(8)
	Lu2	4h	0.17411(2)	0.32589(2)	0.5	0.00760(6)
	Lu3	4h	-0.11531(2)	0.38469(2)	0.5	0.00815(6)
	Si1	4g	0.06628(15)	0.43372(15)	0	0.0064(4)
	Si2	8i	0.19999(16)	0.16501(16)	0	0.0095(4)
	Si3	8j	0.0035(2)	0.15636(17)	0.5	0.0157(5)

Table 2 :

 2 Harmonic displacement parameters for Lu 5 Ir 4 Si 10 .

	Atom	U11 (Å 2 )	U22 (Å 2 )	U33 (Å 2 )	U12 (Å 2 )
	Ir	0.00492(8)	0.00477(8)	0.00620(10)	-0.00036(5)
	Lu1	0.00597(11)	0.00597(11)	0.00879(19)	0
	Lu2	0.00812(9)	0.00812(9)	0.00656(14)	0.00196(10)
	Lu3	0.00802(9)	0.00802(9)	0.00841(14)	0.00024(10)
	Si1	0.0057(5)	0.0057(5)	0.0079(9)	0.0014(7)
	Si2	0.0076(7)	0.0052(6)	0.0158(8)	0.0016(5)
	Si3	0.0316(12)	0.0088(8)	0.0068(7)	-0.0048(8)

Acknowledgments:

We greatly acknowledge the technical assistance of Joël Balay and Abdellali Hadj-Azzem from the Crystal Growth Service at the Néel Institut.