Keywords: programmation générique, DSELs, programmation générative, C++, algèbre linéaire, GPU Generic programming, DSELs, Generative programming, C++, linear algebra, GPU

Titre : Méthodes de génération automatique de code appliquées à l'algèbre linéaire numérique dans le calcul haute performance

Table of contents Introduction

In recent years, architectures and programming languages have gone through major changes and diversification. The constant need for more computational power and more generic software has been the driving force behind these changes. Solving linear systems of equations Ax = b is part of many computational applications leading to a vast amount of research. The development of dense linear algebra has gone in hand with the evolution of computer languages and architectures. The High Performance Computing landscape has gone through many changes with the development of multicore processors and accelerators. One can safely predict that this trend will continue as current architectures are at a crossroad with the limits being reached for the lithography process. At the same time, computer languages have started to diversify with a focus on high level language concepts and functional programming paradigms. Recent languages such as Scala, Dotty, Rust or modern C++ are prime examples of this new direction. Such languages are also starting to become standard in High Performance Computing (HPC) to develop scalable and robust softwares. In this PhD thesis, we present a method of high level programming that takes into account the features of heterogeneous architectures and algorithmic properties to build numerical libraries.

The issues we consider throughout our work are related to software design and optimization for heterogeneous architectures. Software has historically being designed around solving problems for specific architectures, which limits portability and reuse. We propose an approach to develop generic interfaces that support architectural extensions. As the programming model can vary tremendously depending on the hardware, being able to give specifications when necessary and properly categorize hardware for code reuse is essential. This requires high level programming techniques that cannot be found in every language. The concept called Domain Specific Embedded Languages (DSEL) relies on a generic programming language that is flexible and expressive enough to enable the development of a sub-language with a definite syntax. This technique is known to be used in the languages mentioned above. In our work, we employ such programming techniques to design an architecture aware solver for dense linear systems. The resulting codes have been incorporated into a C++ library called NT2 1 . We also consider the issues around batched computations and tensor contractions. This work is done outside of the NT2 library. A detailed explanation on how to develop a simple high level memory container without relying on current software is given. We then develop an optimized matrix-matrix product for small sizes using a tuning approach. In our benchmarks, computations on the GPU are run using the IEEE 754 compliant mode [START_REF] Whitehead | Precision & Performance: Floating Point and IEEE 754 Compliance for NVIDIA GPUs[END_REF] for arithmetic operations and with CUDA ECC activated.

Introduction that takes into account both hardware and algorithmic features to select and generate at compile time the proper LAPACK/MAGMA routine from the high-level C++ code, mapping over 160 kernels.

-We present a multi-stage approach to alleviate the host/device programming model with a cost model for offloading code. This provides adaptable strategy to generate CUDA kernel directly from a single C++ source file containing NT 2 statements.

-Our method takes advantage of modern C++ design for hybrid computations and interface development without relying on current software.

-We describe a strategy using performance analysis based on hardware features and counters to implement efficient algorithms on specific architecture. This work is based on the implementation of a small batched matrix-matrix product that can perform better than state of the art libraries.

Chapter 1

The issue of adapting numerical libraries to parallel architectures Parallelism in todays computer architectures is ubiquitous whether it be in supercomputers or small portable devices such as smartphones or watches. While computers hardware is evolving, computing power as expressed by Moore's law is still the incentive for all manufacturers. If frequency scaling has been a temporary solution, power dissipation and the increasing complexity in CPU design has lead manufacturers to develop more sophisticated architectures. These new architectures that allow to increase the computational power vary from CPU multicores to GPUs, or even integrated system-on-chip (SOC). Exploiting the full performance of such systems for numerical problems has become very challenging.

High level Programming for Scientific Computing

In order to take fully advantage of these new architectures, software libraries are developed to give user access to efficient linear algebra computations. The first major numerical libraries to appear include the Scientific Subroutine Package by IBM [START_REF]System/360 Scientific Subroutine Package[END_REF] in 1968, the Basic Linear Algebra Subprograms (BLAS [START_REF] Dongarra | Basic Linear Algebra Subprograms Technical Forum Standard[END_REF]) or the Linear Algebra Package (LINPACK [START_REF] Dongarra | LINPACK users' guide[END_REF]) in 1979. Each of these packages were already implemented using some form of modular programming for better flexibility. Since then, architecture specific libraries for dense linear algebra have been developed such as LAPACK [START_REF] Anderson | LAPACK Users' Guide[END_REF] for cache-based CPUs, ScaLAPACK [START_REF] Blackford | ScaLAPACK Users' Guide[END_REF] for distributed architectures, PLASMA [START_REF]PLASMA users' guide, parallel linear algebra software for multicore architectures[END_REF] for multicores introducing tiled algorithms, data layouts and more recently MAGMA [START_REF] Baboulin | Enhancing the performance of dense linear algebra solvers on gpus in the magma project[END_REF][START_REF] Tomov | Towards dense linear algebra for hybrid GPU accelerated manycore systems[END_REF] for accelerator-based architectures with GPUs or Xeon Phi.

These evolutions have gone in hand with programming techniques and languages. As we go back, LINPACK was written in Fortran 66 and LAPACK in Fortran 77. The more recent libraries were developed in C for PLASMA and C++ for MAGMA. In each new library the level of abstraction has been raised by using generic programming techniques. However, due to retro-compatibility concerns and the need to maintain a similar interface for the routines, the use of high level programming techniques is still very limited. Moreover, with the increasing parallelism and heterogeneity as well as the ever increasing data-communication cost, the algorithms were completely modified and redesigned to take advantage of each new architecture. The disparity between these libraries that target different architectures illustrates one of the major issues in designing optimized linear algebra software.

Parallel architectures are essential for scientific applications due to the high computing power they provide. Using efficiently these architectures is a difficult task due to the diversity of the hardware and the large amount of software solutions. Algorithms for CPU are still evolving to better adapt to architectures. For examples, the ReLAPACK [START_REF] Peise | Recursive algorithms for dense linear algebra: The relapack collection[END_REF] library has been developed using recursive programming techniques and can outperform LAPACK. It is in this context that we place our work. We focus on architecture aware software design for scientific computations.

In this chapter, we give an overview of the different architectures and programming techniques that are the premises for our work. After explaining why high level programming matters in Section 1.1, we present the algorithms and solvers for dense numerical libraries in Section 1.2. We then detail in Section 1.3 the different architectures and the impact they can have on the different algorithms. Following this, we present the programming tools and techniques in Section 1.1 that are used to develop high level libraries. Finally, we discuss in Section 1.5 about the language and library we used for our work.

High level Programming for Scientific Computing

The evolution of semiconductor technology is dramatically transforming the balance of future computer systems, producing changes at every level. From the point of view of numerical libraries, and the myriad of applications that depend on them, Chapter 1. The issue of adapting numerical libraries to parallel architectures having more computational power allows to work on bigger simulations. These simulations can take up to months of continual computation depending on size of the problem. For this reason, having robust frameworks has become essential. Developing fast applications in a reasonable amount of time requires a multidisciplinary effort incorporating Domain Specific Languages, code generation and optimization techniques, domain science and application-specific numerical algorithms. As computational power and architectures are growing each year, it has now become a major concern to properly design libraries. In 2008 the first Petaflop machine entered the Top 500 while a research article was already explaining the challenges for ExaScale computing [START_REF] Bergman | Exascale computing study: Technology challenges in achieving exascale systems[END_REF] and the emerging architectures. Modularity is now one of the main concern when writing software frameworks. Modern software such as GEANT4, a simulation toolkit for simulating the passage of particles through matter [4], or ROOT [START_REF] Brun | Root-an object oriented data analysis framework[END_REF], for High Energy and Nuclear Physics, are developed using high level languages like C++ and can integrate other languages such as Python or R.

Exploring the use of accelerators as a general purpose co-processor to increase the computational power has been an active research field since the release of the CUDA language in 2007. Dedicated libraries for such architectures have emerged such as MAGMA, CUBLAS by NVIDIA TM or ACML by AMD TM but these are architecture specific. The MAGMA library exists in different versions depending on the hardware vendor. For major frameworks such as GEANT4, we can find research studies explaining how to implement algorithms for new architectures such as Xeon Phi [START_REF] Schweitzer | Performance evaluation of multithreaded geant4 simulations using an intel xeon phi cluster[END_REF]. These changes require advanced methodologies to compile and execute any scientific program. The article also shows that applications need to be heavily tuned to increase performance. The scientific community has a wide range of domains that require dedicated scientific applications to solve the programming challenges. Having an adaptable library does not need to be reimplemented for each new architecture is one of the main challenge of Scientific Computing.

Solving dense linear systems

Many scientific applications whether dense or sparse require at some point to solve a linear system of equations Ax = b. Examples range from fields such as Fluid Dynamics or Electromagnetism using Boundary Integral equations to tensor-based simulations [START_REF] Kolda | Tensor decompositions and applications[END_REF]. Due to this, linear algebra algorithms must be robust, accurate and scalable to fulfill the needs of High Performance Computing.

This notion of numerical stability [START_REF] Higham | Accuracy and Stability of Numerical Algorithms[END_REF] plays a crucial role in designing algorithms as the result obtained by a numerical algorithm does not correspond to the solution of the problem in exact arithmetic. The different sources of errors may come from the input data of the algorithm, which may be caused by prior computations or measurement errors, which are themselves limited to the machine precision currently under the IEEE 754 norm [START_REF] Kahan | Ieee standard 754 for binary floating-point arithmetic[END_REF]. Also, errors may be caused by approximations made by the algorithm. Depending on the desired accuracy, several algorithms can be chosen to solve a linear system. We can find two major approaches to solve linear systems : direct or iterative methods. A direct method uses a finite sequence of operations to provide an exact solution x if there are no rounding errors. The common method used to solve such systems is called Gaussian elimination. It consists in adding the coefficients of one equation to the others in order to eliminate a variable and continue this process until only one variable is left. Iterative methods start with an approximation of the solution x 0 and successively compute a sequence of approximations x i to improve the solution. In our work we focus on direct methods which are generally used to solve dense linear systems. Direct methods provide higher numerical precision and a good granularity of computations making it easier to exploit the latent parallelism. We note the main decompositions that exist [START_REF] Golub | Matrix Computations[END_REF] : LU, Cholesky, QR, SVD, LDL T . LU: Decomposes a general matrix A as L × U where L is unit lower triangular and an U is upper triangular (about 2 × n 3 /3 flops).

Cholesky: A symmetric positive definite (SPD) matrix A is factored as A = LL T where L is lower triangular (about n 3 /3 flops).

QR: An m-by-n matrix A is factored as A = QR with Q an m-by-m orthogonal matrix and R an m-by-n upper triangular matrix (2n 2 (mn/3) flops).

SVD: Decomposition of an m-by-n matrix A as A = U ΣV T , U an m-by-m orthogonal matrix, Σ an m-by-n diagonal matrix (where the diagonal contains the singular values of A) and V an n-by-n orthogonal matrix V (the cost is about mn 2 flops to compute Σ).

LDL T : factorization (using symmetric pivoting) is P AP T = LDL T where P is a permutation matrix, A is a symmetric indefinite square matrix, L is unit lower triangular and D is block-diagonal, with blocks of size 1 × 1 or 2 × 2 (about n 3 /3 flops).

LU factorization

The LU factorization [57, p.111] is a modified Gaussian elimination without pivoting that decomposes a nonsingular matrix A into a lower triangular matrix, and an upper triangular matrix U. The number of floating point operations of the LU decomposition for a square matrix n is about 2n 3 /3. For dense matrices, the LU factorization is usually performed in place. This means that the output factors L and U overwrite the matrix A during the factorization.

Chapter 1. The issue of adapting numerical libraries to parallel architectures Algorithm 1 LU factorization in place, no pivoting Input: A is n × n matrix 1: for k = 1 : n -1 do 2:

for i = k + 1 : n do for i = k + 1 : n do 6:

for j = k + 1 : n do 7:

A(i, j) = A(i, j) -A(i, k) * A(k, j)

8:
end for 9:

end for 10: end for

The LU factorization is generally implemented using pivoting strategies because there are several important issues that can appear without pivoting. For example, if a zero is found on the diagonal of the matrix, a division by zero will occur leading to a numerical error. Also, if there are elements of small magnitude on the diagonal, the factors will end up growing too much also leading to erroneous results. We note P A = LU , with P a permutation matrix, when using partial pivoting in LU factorization. The stability of Gaussian elimination can be measured with the growth factor [START_REF] Higham | Accuracy and Stability of Numerical Algorithms[END_REF]. The growth factor measures the ratio between the entries of the matrix after the elimination steps and the original entries. The growth factor of a square matrix A under Gaussian elimination is defined as:

g n (A) = max i,j,k |a (k) ij | max i,j |a i,j | ,
where a

(k) ij is the element of index (i, j) after the step number k of the elimination [START_REF] Kincaid | Numerical analysis: mathematics of scientific computing[END_REF]. The most common pivoting strategy used is called partial pivoting. It consists in permuting the rows but not the columns to ensure that the pivot is the largest entry in its column. For a step j of the Gaussian elimination, partial pivoting permutes the rows such that |a j,j | |a i,j | for all i > j . This guarantees that ||L|| ∞ 1. Its growth factor upper bound is 2 n-1 which can be reached for certain problems [START_REF] Foster | Gaussian elimination with partial pivoting can fail in practice[END_REF].

for i = k + 1 : n do 4: if |A(i, k)| > |A(index, k)| then 5: index = i 6: end if 7:
end for

for i = k + 1 : n do 10: A(i, k) ← A(i, k)/A(k, k)
for i = k + 1 : n do 13:
for j = k + 1 : n do 14: The drawback of pivoting comes from the communication cost due to the comparisons performed to find the pivot and to the resulting swap of rows. For a square matrix, the partial pivoting requires O(n 2 ) comparisons. Other well know pivoting strategies are total pivoting (O(n 3 ) comparisons) and rook pivoting (O(1.5n 3/4 log n ) comparisons).

A(i, j) = A(i, j) -A(i, k) * A(k, j)
The LU factorization is a classical technique that is used with different implementations in high performance linear algebra libraries. We can find applications with different pivoting strategies, precision, with iterative refinement or even in mixed precision. For this reason, we have used the LU factorization and its various forms as an important part of our generic solver for dense linear systems that we are going to detail in Chapter 2.

QR factorization

The QR factorization [57, p.246] decomposes a matrix A into the product of an upper triangular matrix R and an orthogonal matrix Q. This is the main method used to solve the least square problem which consists in finding a minimal solution to a system Ax = b. This means solving ||Ax-b|| p , A being an m by n matrix, for a defined p. There are multiple methods to compute the QR factorization of a matrix. Here we will focus on factorization with the Householder transforms which can be found in numerical libraries. We will described in Chapter 2 how we implemented a versions of mixed-precision semi-normal equations based on this QR factorization.

During the QR factorization, applying the Householder transformations on the matrix A will remove every entry below the diagonal. By repeating this process n architectures times, we obtain the upper triangular matrix R such that :

H n H n-1 ...H 1 A = R.
We note the form of an orthogonal matrix computed durign the factorization as :

H = I -2 ||u|| 2 uu T .
If we set Q = H 1 ...H n , we can then obtain A = QR which means we have [v, β] = house(A(j : m, j))

Q T A = R.
3:

v j = [0...0, 0, 1, v j j+1 ....v j m ] T 4: β j = 2 1+||A(j+1:m,j|| 2 ) 5: A(j : m, j : n) = (I -βvv T )A(j : m, j : n) 6: if j < m then 7: A(j + 1 : m, j) = v(2 : m -j + 1) 8:
end if

9: end for

We note that this method is known to be numerically stable [START_REF] Wilkinson | Rounding erros in algebraic processes[END_REF].

Solvers in LAPACK and MAGMA

Solvers for linear systems can be found in LAPACK and MAGMA constitute the building blocks to write a generic library for dense linear algebra. The routines given in these libraries can be distinguished by the static properties of the matrix. These properties correspond to the structure of a dense matrix which are for instance general, band, symmetric, positive definite or tridiagonal. Other structures can be considered but they may require new implementations to be written. The more common class of linear systems correspond to general dense matrices. This study focuses on dense matrices which already cover a wide range of applications.

What we describe below is the guideline of the main algorithms that are implemented. This however does not represent the implementation since it will depend of the structure of the matrix and also the target architecture. LAPACK being for CPU, the implementation is more straightforward and matches the algorithms while MAGMA will tend to be more specific. Indeed, since MAGMA is designed for GPU, it is important that the transitions between CPU and GPU is smooth to ensure minimal communication cost.

Solving dense linear systems 13 1.2.3.1 General solvers

The solver used in LAPACK/MAGMA for general linear systems (xgesv, the letter x indicating the precision, either real (s), double-real (d), complex (c) or doublecomplex (z) ; ge standing for general) is based on the LU decomposition with partial row pivoting. This technique, called Gaussian elimination with partial pivoting is implemented in most numerical libraries including LAPACK and MAGMA and is used in the LINPACK benchmark for the TOP500 list.

To solve the linear system Ax = b, we compute a factorization P A = LU with P permutation matrix, L lower triangular and U upper triangular. Then the solution

x is computed by solving successively the triangular systems Ly = P b followed by U x = y. In the case of rectangular matrices, the routine xgelsy uses the QR decomposition of the matrix (possibly with column pivoting) to solve the linear least square problem [START_REF] Björck | Numerical methods for least squares problems[END_REF] : The interest of these factorizations is that we can reuse the obtained factors for other computations which can also be very expensive such as computing the condition number [START_REF] Higham | Accuracy and Stability of Numerical Algorithms[END_REF], or performing iterative refinement [START_REF] Demmel | Error bounds from extra-precise iterative refinement[END_REF]. In the case of a square matrix, we will identify three types of solvers that are available in LAPACK and will mention others that can added. In the case of a rectangular matrix we can identify a unique solver, which is the routine xgelsy.

min x∈R n Ax -b 2
The first one is the routine xgesv to compute a solution using the LU factorization. This can be seen as the standard case when there is no particular information about the system. Other cases include iterative refinement and mixed precision.

Chapter 1. The issue of adapting numerical libraries to parallel architectures

Iterative refinement

The second one is the routine xgesvx using iterative refinement which enables it to improve the accuracy of the computed solution. Iterative refinement is a method that produces a correction to the computed solution by iterating on it. Each m th iteration in this process consists of computing the residual r m = b -Ax m , solving the new system Ad m = r m , and adding the correction

x m+1 = x m + d m .
Algorithm 6 One iteration of iterative refinement

1: Compute the residual r = b -Ax. 2: Solve the new system Ad = r. 3: Add the correction x = x + d.
This method is computationally more expensive but can provide better solution for ill-conditioned systems.

Mixed precision algorithm

The third method is the solver dsgesv that uses mixed precision iterative refinement [START_REF] Baboulin | Accelerating scientific computations with mixed precision algorithms[END_REF]. This method takes advantage of the fact that computing in single precision is in general much faster than computing in double precision. The factorization (most expensive part of the computation) is computed in single precision (ε s ) (if the matrix is not ill-conditioned) and the solution is refined in double precision (ε d ) using the classical iterative refinement method.

Algorithm 7 Mixed precision, iterative refinement

1: Compute P A = LU (ε s ) 2: Solve Ly = P b (ε s ) 3: Solve U x 0 = y (ε s ) do k = 1,2,... 4: r k = b -Ax k-1 (ε d ) 5: Solve Ly = P r k (ε s ) 6: Solve U z k = y (ε s ) 7: x k = x k-1 + z k (ε d ) check convergence
These methods can be found in current open source libraries. There are also other interesting methods that can provide satisfactory performance. An example is a method that accelerates linear system solutions using randomization techniques [START_REF] Baboulin | Accelerating linear system solutions using randomization techniques[END_REF], in which the system is modified by a multiplicative preconditioning based on Random Butterfly Transformations [START_REF] Parker | Random butterfly transformations with applications in computational linear algebra[END_REF] that removes the need for pivoting (which is an expensive process). Another method is based on communicationavoiding algorithms that minimize the cost of pivoting [START_REF] Baboulin | A class of communication-avoiding algorithms for solving general dense linear systems on CPU/GPU parallel machines[END_REF][START_REF] Grigori | CALU: a communication optimal LU factorization algorithm[END_REF].

SIMD extensions

Single instruction multiple data (SIMD) architectures correspond to processing units that can perform an operation on multiple data simultaneously (see Figure 1.1). SIMD extensions are designed to help processors exploit the latent data parallelism available in many applications at a smaller cost than duplicating Arithmetic Logic Units (ALU).

Multiple constructors have implemented vector extensions. SIMD extensions started mostly with Sun Microsystems Visual Instruction Set (VIS) for the SPARC V9 [START_REF] Kohn | The visual instruction set (VIS) in UltraSPARC[END_REF] and Hewlett-Packard Multimedias Acceleration eXtensions (MAX) for the PARISC instruction set [START_REF] Lee | Accelerating multimedia with enhanced microprocessors[END_REF]. This lead other constructors to develop their own SIMD extension set. [START_REF] Peleg | MMX technology extension to the Intel architecture[END_REF] used floating point registers instead of having specific vector registers which disabled scalar computations. Later on, they developed their first processors with dedicated vector registers (Pentium III) with 64 bytes registers. Since then, manufacturers have continued to increase the size of dedicated SIMD registers which now go up to 512 bytes vectors for Intel Xeon Phi coprocessors [START_REF] Chrysos | Intel R xeon phi TM coprocessor-the architecture[END_REF] (which corresponds to a line of L1 cache).

Exploiting the parallelism offered by these extensions is critical for high perfor-Chapter 1. The issue of adapting numerical libraries to parallel architectures mance computing as it can theoretically increased performance up to 8 times with 512 bytes registers for double precision operations. Using SIMD extensions can also increase the bandwidth. Linear algebra can benefit greatly from SIMD extensions as BLAS operations which are fundamental for most algorithm are vector based.

Multicore systems

The processor manufacturers faced a technology limitation to their traditional approach of boosting clock speeds and increasing throughput due to power dissipation problems. This lead to the development of multicore architectures and hyperthreading as a solution which is a turning point in terms of software development.

While constructors could increase the frequency in processors freely, the performance of a program would naturally increase which each new generation of processor. This went on until the end of the Intel Pentium 4 processors in 2006 with frequencies around 3.6 GHz. By then multicores were already available on the market but were not yet generalized. The common concept behind a multicore processor is to let each core have an L1 independent cache and put a shared L2 cache on the die which is the interface to the main memory. Figure 1.2 show a multicore systems with two levels of cache. Modern CPUs generally have two levels of cache in each core and a shared L3 cache. Having multiple CPU cores on the same die allows for a highly efficient cache coherent system at a much higher clock rate than if the signal had to travel off chip. Front side bus Most devices ranging from desktop computers and supercomputers to smartphones are based on a multicore design (see Figure 1.2). The number of cores can vary greatly depending on the solution going up to 22 cores on recent server based Intel Xeon processors family. Other parallelization techniques based on simultaneous multithreading (SMT) to improve the efficiency of superscalar CPUs can be found. As an example, Hyper threading (HT) consists in sharing resources by having two logical core for one single physical core. Another solution is the Non Uniform Memory Access (NUMA) architecture. Such systems are generally built around having multiple processors with each processor having its own memory bank (see Figure 1.3). Each processor is able to access the memory bank of the others but latencies will vary depending on the locality of each processor. To support this system, the operating system provides a virtual memory address space for all the processors in the NUMA node.

Graphics Processing Units (GPU)

Applications are becoming more demanding in workload and complexity as architectures become more powerful. CPUs have shown their limits in terms of computational power as they tend to be more general. This technology gap has led to the development of new architectures that are more specialized and can therefore provide more computational power for specific applications. GPUs are an example of a new type of architecture that has been emerging during the last decade.

A GPU is a specialized hardware to accelerate image display and processing. The first programming language for such architectures was released in 2007 by NVIDIA. The Compute Unified Device Architecture (CUDA) is a platform designed to work with scientific programming languages such as C, C++ and Fortran. CUDA has a GPUs architectures are very different from multicore CPUs and provide more computational power and better energy efficiency. The number of processing units available in a GPU is much higher than in a CPU (see Figure 1.4). A Tesla K40 has 2880 processor cores available but each runs at a lower frequency (745 GHz). Each of the cores can be seen as simple ALU that can only be used for computations. Cores are grouped together in what is called streaming multiprocessor (SM) with 192 CUDA cores for Kepler cards. Each SM has registers, an shared L1 cache and has a shared L2 cache between each SM with a maximum of 15 SM for Kepler cards. Similarly to CPUs, the number of cache levels on the GPU has grown to mitigate the high latency from the GPUs main memory.

Contrary to the CPU where data goes through each cache level for Intel processors, the data on the GPU does not necessarily stay at each cache level. Some data management therefore needs to be done to optimize data movement. The GPU has a thread scheduler and control units in each SM(see Figure 1.4). The execution model of the GPU is based on deploying a large number of threads (in the order of millions) as they have no creation cost compared to a CPU. The scheduler will then distribute the threads minimizing the number of stall cycles. This creates a high parallel programming model where a high number of compute units will work in parallel. Communications with the CPU use the PCI express bus. The speed of the PCI bus can vary depending on the technology, with a maximum of 6 GB/s for PCI express gen 2 (16x) and 12 GB/s for PCI express gen 3. Depending on the application, the PCI express bus can be a bottleneck for memory-bound problems due to data transfers ( see Figure 1.6). With each new architecture generation, GPU architectures often go through many changes but manage to maintain the same programming model.

Accelerated Processing Units

GPUs and accelerator based hardware are a major component of todays systems. However, using accelerators to their full potential can be impossible due to the limitations of the PCI express bus. Furthermore, GPUs are not stand alone processors and requires a CPU to function. This has lead to the development of System On Chip (SOC) hardware where a CPU and a GPU are combined on a single die. NVIDIA and AMD were the first manufacturers to develop this type of hardware with the NVIDIA Tegra brand in 2008 and AMDs Accelerated Processing Unit (APU) in 2011. Later on, Intel released its Integrated Graphic Processors for the Chapter 1. The issue of adapting numerical libraries to parallel architectures Intel core processors (Figure 1.8). As there is no official term for such hardware, we will refer to them as Integrated Graphic Processors (IGP). IGPs do not require a PCI Express bus as the CPU and GPU are on the same die (see Figure 1.6).

L2 cache

Registers SM-0 As both hardware have a direct access to the main memory, the GPU has a higher memory bandwidth. The first implementations of IGPs by AMD did not have a unified memory system between the GPU and CPU. This meant that they could not access the same address space limiting the interactions between each processor. With it's second generation of IGP called Kaveri which was released in 2014, AMD implemented a zero-copy model where CPU and GPU can access the memory with the same address space. This means pointers can be freely passed between CPU and GPU, thus preventing actual deep copies. NVIDIA also released an IGP with zero-copy through the Tegra K1 and later on the X1. AMD IGPs are coupled with AMD CPUs while NVIDIA uses ARMs CPUs. It is important to note that current IGP have very poor support for double precision computations.

L1

Programming tools and techniques

Section 1.3 described the current architectures without giving too much details about their programming models or the tools available for such platforms. In the next section, we are going to describe some of these tools and techniques often used to develop code for parallel architectures. Each tool will be described by its support for different architecture and expressiveness while techniques will explain how to build a tool for more expressiveness and architectural support. In section 1.3.1, we introduced the concept of SIMD extensions. Exploiting SIMD units in todays code is essential to reach the peak performance of a modern CPU. The most common method to exploit these extensions is by using the function intrinsics available. These functions are available in C but have been ported to other languages even web oriented ones such as JavaScript. The low level intrinsics in C are still used in every day programming and tend to be quite verbose. Depending on the size of the SIMD registers or the architecture, the intrinsics are also not the same. This causes some API problems and render complex the generalization of SIMD extensions. In Listing 1.1 we can see an example of the intrinsic for the multiplication of two 256 bits vectors. The verbosity of SIMD instructions has been diminished in recent years and computational intrinsics such as addition, multiplication, etc.. can be replaced by its operator as seen in Listing 1.2. On ARM NEON, the SIMD maximum size is 128 bit. It is therefore not possible to do the same operation as above but we can do a similar one with 128 registers. It is also possible to simply call the operator instead of the full intrinsic (Listing 1.3). In these examples we do not show the load and store operations needed to fill SIMD registers. For Intel processors, data also needed to be aligned on older processors to properly use SIMD extensions. On modern processors, unaligned load and store are quite sufficient to reach good performance most of the time as long as data doesn't cross cache lines. With the introduction of ARM supercomputer based system, recent studies [START_REF] Rajovic | Tibidabo: Making the case for an arm-based hpc system[END_REF] have shown how to improve energy efficient and design of HPC systems on low-power cores. Studying how these ARM based system can perform is therefore of interest.

Compilers for low level languages such as C, C++ can generate SIMD code through their auto-vectorization techniques. The programmer can write a simple code that will be translated to SIMD code during the compilation phase. GCC [START_REF] Nuzman | Multi-platform auto-vectorization[END_REF] is an example of a compiler using a vectorizer. It is also possible to use Just In Time (JIT) compilation to vectorize code with tools such as VaporSIMD [START_REF] Nuzman | Vapor simd: Auto-vectorize once, run everywhere[END_REF] or use library based solutions such as Boost.SIMD [START_REF] Estérie | Boost.SIMD: generic programming for portable SIMDization[END_REF] or VC [START_REF] Kretz | Extending C++ for explicit data-parallel programming via SIMD vector types[END_REF]. The common problem with generating SIMD comes from how patterns are recognized. It is not always obvious how to vectorize a problem especially if data is not contiguous. Compiler can provide directives such as #pragma simd for ICC or #pragma GCC ivdep for GCC. These directives usually tells the compiler that their is no loop-carried dependencies which prevents consecutive execution of a loop with SIMD. Complex operations such as shuffle or intra-registers are hard for compiler to generate as they tend to be very specific. Some compilers also use SIMD extensions in their implementation of the standard library such as GCC.

Multicore support

As we have seen in section 1.3.2 modern CPUs contain multiple cores on the same die. Most programming languages do no support multicore programming by default meaning the programmer must do the effort. We can still cite a language such as Chapel by CRAY which is available on Github. Chapel [START_REF] Chamberlain | Parallel programmability and the chapel language[END_REF] is designed around multithreaded execution via high-level abstractions for data parallelism, task parallelism, concurrency and nested parallelism (see example in Listing 1.4).

u s e C y c l i c D i s t ; // u s e t h e C y c l i c d i s t r i b u t i o n l i b r a r y c o n f i g c o n s t n = 1 0 0 0 0 0 ; // o v e r r i d e d e f a u l t u s i n g . / a . out --n=<v a l > f o r a l l i i n { 1 . . n} dmapped C y c l i c ( s t a r t I d x =1) do w r i t e l n ( " H e l l o from i t e r a t i o n " , i , " o f " , n , " r u n n i n g on node " , h e r e . i d ) ;

Listing 1.4: Chapel Hello World

The main tools for parallel programming on multicore systems are Pthreads, OpenMP, TBB and MPI.

Shared memory multi-threading (pthreads) is an implementation of the POSIX 1003.1c Standard. It is based on a model where threads are spawned on the fly while each thread is given a task at its creation. The threads can then communicate through shared memory variables. The Pthread standard provides functions for thread synchronization -barriers, semaphores, and critical regions in its API. Pthread does not have a specific programming model and allows for task or data parallelism but requires a complex semantic and is error prone. It is a low-level explicit model where the programmer has to allocate threads, synchronize them, and then recompose them. Higher level implementation of Pthread can be found in programming language such as the C++11 std :: thread.

OpenMP is a high level parallel programming model. It uses pragmas to add parallelism on a sequential implementation of an application. This requires the code section to be run in parallel to be dependency free between iterations. OpenMP also allows for task based parallelism to express dependencies. It uses a runtime library for thread management and scheduling techniques. The standard is constantly evolving and since OpenMP 4.0 supports extensions for accelerators based architectures.

Intel Thread Building Blocks (TBB) [START_REF] Pheatt | Intel R threading building blocks[END_REF] is a C++ template library for data parallelism. It is based on the concept of tasks instead of threads which are higher level entities. TBB uses parallel skeletons constructs as building blocks. We note for example parallel f or, parallel r educe, parallel s can and parallel s ort. The task based abstraction is an efficient high level model and provide simple abstractions with the help of lambda expressions since C++11.

Message passing (MPI) is a popular programming model originally designed for distributed computing. It provides a low-level explicit model where each core of a architectures processor can be addressed independently having its own memory pool. MPI is more of a task-based parallelism model but also supports data parallelism by deploying similar MPI processes on multiple data. Recent MPI standard have added special communicators for shared memory programming.

Accelerator support

In section 1.3.3 we have explained how the GPU programming model is based on offloading data and using data parallelism. The standard programming tools used to program accelerators are CUDA and OpenCL.

CUDA is the programming language for NVIDIA graphic cards. It works for GPU as well as NVIDIA IGPs without having to change the device code. As explained before, the NVCC compiler supports a subset of the C++ language. To execute code on a GPU, the programming model includes a host code for the CPU and a device code which will be executed on the GPU. CUDA is an expressive language and simplifies the use of accelerators following the C++ standard evolutions. Unfortunately, it is not cross platform which can be a heavy restriction for portable code.

OpenCL is the portable solution for programming accelerators. It can work with CPUS and GPUs indistinctly but a single code cannot be highly efficient on both a CPU and a GPU. Similarly to CUDA, an OpenCL program has two parts. A kernel which will be executed on one or more OpenCL devices (similar syntax to CUDA ) , and a host program which is executed on the host. The Host code for OpenCL is very verbose as it is designed to work for every platform and several libraries try to alleviate this complexity. One of the main problem with OpenCL is that it is mostly C compliant and is therefore hard to mesh with high level programming without having to redesign around OpenCL.

Domain Specific Languages

By definition, a Domain-Specific Language (DSL) is a computer language designed to fit a specific application domain. This is contrary to a general-purpose language, which is widely applicable across multiple domains and thereof lacks features for a specialized domain. These languages focus on simplifying specific tasks for developers by removing the complexity of a programming language. Adding this layer of abstraction helps developers to focus on the specificities of the domain described by the language, reducing the gap between domain experts and simple developers. DSL also have the effect of increasing the maintainability and ease the validation of codes.

There are a multitude of DSLs available ranging from domain such as editing documents to expressing physical problems. Tex, a markup language, is an example of a widely used DSL for document editing created by Donald Knuth in 1977. Tex offers control over document formating with commands to help properly structure a text without having to focus too much on it. Other examples of well know DSLs are the Structured Query Language (SQL) to exploit databases or Hyper Text Markup Language (HTML) to create web pages.

In numerical computing, MATLAB

TM [START_REF]MATLAB User's Guide. The mathworks[END_REF], a multi-paradigm numerical environment, is the most popular DSL. The original purpose of MATLAB what to ease the use of the advanced LINPACK library without having to know any FORTRAN. MATLAB is mainly focused on matrix manipulations by exploiting optimized libraries such as LAPACK for CPU while providing utility tools for function plotting and data management. It is also interfacing with other languages such as C, C++, Java, FORTRAN or Python.

DSLs are usually part of a developing environment as they cannot be used as stand alone. To implement a DSL, one needs to write a compiler/interpreter for its language as it is not part of a general purpose computing language. This is a heavy cost as developing a DSL require more than just writing a language and can often suffer from performance issues due to the lack of optimizations. Another approach that we use in our work is based on Domain Specific Embedded Languages (DSEL) that we describe next.

Domain Specific Embedded Languages

Domain Specific Embedded Languages (DSELs) are a subclass of DSL that rely on an existing general-purpose language to host it. DSELs then reuse the host language syntax and tool ecosystem to be compiled or interpreted. This means inheriting from the languages constructs (structure, conditions, functions, etc) while adding primitives that allow programmer to have access to a higher level of abstraction. Differences and advantages of DSELs over DSL are discussed by Abrahams [3,p.226] There are several forms of DSEL available in different languages. Javascript is in a sense an embedded language in an HTML page, while HTML is a DSL. In scientific computing, DSEL rely on a generic language that is flexible and expressive enough to enable the development of a subset language with a definite syntax. This technique is know to be used in languages such as Lisp, Haskell, SmallTalk, C++, Scala and more. It is not always possible to reproduce a specific syntax in a general purpose language and sometimes a symbol needs to be replaced. For example, in C++ it is not possible to reuse the colon operator (:) as defined in MATLAB. This is one of the main shortcomings of a DSEL with the fact that it must be implemented inside a library and cannot be added to a language. Designing a DSEL is however often easier as they reuse existing compiler and rely on domain dependent analysis to generate efficient code.

We will describe how we can design a DSEL in C++ using meta-programming techniques in Chapter 2. Our work will focus on building a generic solver for dense linear algebra with multi-architectural support.

Multi-stage programming

Multi-stage programming (MSP) is another approach than DSL and DSEL to increase the expressiveness of a programming language. MSP is defined as a metaprogramming technique which consists in writing a program that will be able to generate new code. This means the state of the program will be modified during the compilation phase by the code meta-programmed.

MSP is usually applied by adding specific language extensions to trigger a new compilation phase. As an example, MetaOcaml [START_REF] Taha | Metaocaml-a compiled, type-safe, multi-stage programming language[END_REF] is a multi-stage language based on Ocaml [START_REF] Leroy | The objective caml system release 3.12. Documentation and user's manual[END_REF] with three basic constructs for runtime code optimization. It adds type code values and constructs to build them called quoting and splicing. The generated code with MSP can then be stored in a file or compiled and linked directly to a Ocaml program enabling run-time code optimization.

More recent adoptions of MSP for parallel computing include techniques like Lightweight modular staging (LMS) [START_REF] Rompf | Lightweight modular staging: a pragmatic approach to runtime code generation and compiled dsls[END_REF] in Scala. Instead of using quasi-quotation LMS is based on types and uses an environment for MSP. Similarly to DSEL, using MSP with a host language requires this language to be flexible enough to adapt to MSP. Other approaches include Language Virtualization [START_REF] Chafi | Language virtualization for heterogeneous parallel computing[END_REF] (base on LMS) or Terra [START_REF] Devito | Terra: a multistage language for high-performance computing[END_REF] using language extensions for HPC in Lua.

We will describe our work on MSP in Chapter 3 based on the C++ language and the CMake tool for the multi-level code generation.

Considered language and library

C++ language and programming techniques

Generic programming is a paradigm that aims at achieving reusable, adaptable libraries. Pioneered by Alexander Stepanov and David Musser, Generic Programming obtained its first major success when the Standard Template Library became part of the C++ Standard. The process followed in generic programming is to first find the similarities among different implementations of a same algorithm in the form of concepts [START_REF] Gregor | Concepts: linguistic support for generic programming in c++[END_REF] , and then provide an abstraction that can match all of these implementations. Concepts are a subtle mix of constraints and axioms used to represent the problem. The constraints define the statically evaluable predicates on the properties and syntax of the system, while axioms state the semantic informations about the types requirements, providing a better understanding of the problem. This process, called lifting, is then repeated until the generic algorithm has reached a suitable level of abstraction, where it provides maximal re-usability without sacrificing performance:

• Step 1 : Find similarities among different implementations • Step 2 : Provide an abstraction that can match all of those implementations • Step 3 : Iterate until a satisfactory level is reached 1.5. Considered language and library Generic Programming and Concepts provide support for Concept Checking, which allows for a limited support of parametrized types constraints [START_REF] Siek | Essential language support for generic programming[END_REF], simplifying error checking at compile-time. In C++, concepts are being studied to be integrated in the standard as extension to the template system. They have been published as an ISO Technical Specification ISO/IEC TS 19217:2015.

Generative Programming is about bringing the benefits of automation to software development. Following this paradigm, any complex software system can be broken down to a list of interchangeable components which tasks are clearly identified and a series of generator that combine components by following rules set by an a priori domain specific analysis. The application of generative programming to parallel programming and hardware synthesis can be seen as a way to leverage the classical limitations of classic, library based tools by encapsulating expertise at the source code level and not only at the binary level. One of the most classical techniques for generative programming is Template Meta-programming.

Template Meta-programming is a meta-programming technique in which templates are used by a compiler to generate temporary source code, which is merged with the rest of the source code and then finally compiled. The output of these templates includes compile-time constants, data structures, and complete functions. The use of templates can be thought of as compile-time execution. This technique is used by a number of languages, the most well-known being C++ [START_REF] Abrahams | C++ template metaprogramming: concepts, tools, and techniques from Boost and beyond[END_REF], D [START_REF] Bright | D language Templates revisited[END_REF], Haskell [START_REF] Sheard | Template Meta-programming for Haskell[END_REF] and OCaml [START_REF] Taha | A gentle introduction to multi-stage programming[END_REF].

The NT 2 library

NT 2 [START_REF] Estérie | The numerical template toolbox: A modern c++ design for scientific computing[END_REF] is a numerical computing C++ library implementing a subset of the Matlab language as a DSEL . The work that we described in Chapter 2 and 3 are extensions for this library. NT 2 simplifies the development of data-parallel applications on a large selection of architectures currently including multi-core systems [START_REF] Tran Tan | Automatic task-based code generation for high performance domain specific embedded language[END_REF] with SIMD extensions. Simply put, a Matlab program can be converted to NT 2 by copying the original code into a C++ file and performing minor cosmetic changes (defining variables, calling functions in place of certain operators). NT 2 also takes great care to provide numerical precision as close to Matlab as possible, ensuring that results between a Matlab and an NT 2 code are sensibly equal.

Internally, NT 2 is designed to leverage the well known Expression Templates [START_REF] Veldhuizen | Expression templates[END_REF] C++ idiom to build at compile time a flexible representation of the abstract syntax tree of any C++ expression containing at least one NT 2 component. This compile-time tree is then transformed in actual code to be executed on a parallel system. Contrary to other libraries based on the same technique, NT 2 relies on Boost.Proto , an external Expression Templates system to handle the creation and transformation of Abstract Syntax Trees (AST) [START_REF] Niebler | Proto : A compiler construction toolkit for DSELs[END_REF]. Boost.Proto allows us to replace the direct walk-through of the compile-time AST done in most Finally, the other main difference between NT 2 and similar tools is the fact that architectural specificities of the generated code are handled by an execution model based on Algorithmic Skeletons [START_REF] Cole | Bringing skeletons out of the closet: a pragmatic manifesto for skeletal parallel programming[END_REF]. Those skeletons simplify the extension of NT 2 for various hardware by separating the concerns of optimizing the AST code for different type of architecture. The main element of NT 2 is the table class. table is a template class that can be parametrized by its element type and an optional list of settings. Instances of table behave like Matlab array and supports the same operators and functions. NT 2 covers a very large subset of Matlab functionality, from standard arithmetic, exponential, hyperbolic and trigonometric functions, bitwise and boolean operations, IEEE related functions and of course linear algebra.

// Matlab : A = 1 : 1 0 0 0 ; t a b l e <double> A = _( 1 . , 1 0 0 0 . ) ; // Matlab : B = A + randn ( s i z e (A1) ) ; t a b l e <double> B = A + randn ( s i z e (A1) ) ; // Matlab : r = s q r t ( sum ( (A ( : )-B ( : ) ) . ^2 ) / numel (A) ) ; d o u b l e r = s q r t ( sum ( s q r (A(_)-B(_) ) ) / numel (A) ) ; Listing 1.5: NT 2 RMSD Computation

Conclusion 29

Listing 1.5 showcases some NT 2 basic features including the mapping of the colon function (:) to the _ object, various functions, a random number generator and some utility functions like numel or size.

In Chapter 2, we will explain how to build a hybrid programming model and how we have integrated ours inside NT 2 . In Chapter 3, we will tackle the problem of MSP and explain how we can integrate it inside a C++ library using expression templates and the CMake tool.

Conclusion

In this chapter, we have presented the linear algebra algorithms we will work with and their various versions in state of the art high performance dense linear algebra libraries. Then we have described the architectures and their different evolutions that have gone in hand with the development of new libraries to support these architectures.

We have also given a general description of the various programming tools available for these architectures. This is done as a reminder of how complex the architectural landscape and the various tools available is. We have then mentioned high level programming techniques to develop more generic and flexible software that can be adaptable to new architectures. Finally, we have described some of the C++ techniques we rely on and the library in which some of our code is implemented.

From what we have seen, it is clear that the diversity of architectures, which we have not all described, combined with the amount of libraries and programming languages available make for a very large landscape. In this manuscript, we will explain in detail techniques that can alleviate such problems. These techniques are implemented using the C++ language but are not exclusive to this language. As of now, C++ is a widespread high-level language in HPC. It is unclear which language will take this role in the future, but we can say that the programming techniques employed will have a place in it. It is in this context that we present our work.

In the next chapter, we will present one of our contribution in the form of a multi-architectural generic solver for dense linear algebra that is expendable to new architectures. We will show how to design simple GPU containers and a consistent memory model that will also be used in Chapter 3.

Chapter 2

High-level programming for dense linear algebra The increasing complexity of new parallel architectures has widened the gap between adaptability and efficiency of the codes. As high performance numerical libraries tend to focus more on performance, we wish to address this issue using a C++ library called NT 2 . By analyzing the properties of the linear algebra domain that can be extracted from numerical libraries and combining them with architectural features, we developed a generic approach to solve dense linear systems on various architectures including CPU and GPU.

The work presented in this chapter was accepted and presented at Ispa 2015 [START_REF] Masliah | Metaprogramming dense linear algebra solvers applications to multi and many-core architectures[END_REF]. The contributions of this chapter are :

• The definition of a generic hybrid programming model that can be used as a stand-alone

• An architecture aware binding between NT 2 and LAPACK/MAGMA based on type tags to dispatch between the different architectures and runtime backends in an extensible way.

• An implementation of linsolve (in reference to the MATLAB routine) that takes into account both hardware and algorithmic features to select and generate at compile time the proper LAPACK/MAGMA routine from the highlevel C++ code, mapping over 160 kernels. Note that the support of different factorizations (QR, Cholesky, LU, SVD) is also provided in NT 2 to facilitate the development of new solver.

• An application based on a linear least squares solver (using mixed precision) that uses of the already available routines for several architectures in NT 2 . This application can therefore be written with a single interface and automatically generated for CPU or GPU while showing a level of expressiveness similar to MATLAB.

This chapter is organized as follows: in Section 2.1 we give the software context and related work on high level libraries. Then, in Section 2.2, we describe various programming techniques that combine algorithmic and architectural features in libraries. The methods that we used in N T 2 are then introduced. They enable us to achieve re-use and adaptability of library codes while preserving performance. After that, we explain in Section 2.3 the concept behind our GPU container and how it should interact with the CPU. In Section 2.4.3, we detail our methods to develop efficient dense linear algebra software. We present in Section ??, as an example of an application, the code generation of a mixed-precision linear least squares solver. We give performance comparisons on CPU and GPU using respectively to the QR routines from LAPACK and MAGMA. To our knowledge such a solver does not exist in public domain libraries LAPACK, PLASMA [START_REF]PLASMA users' guide, parallel linear algebra software for multicore architectures[END_REF] and MAGMA. Concluding remarks are given in Section 2.5.

Related Work

A major concern when developing dense linear algebra software is to propose a user-friendly Application Programming Interface (API) that provides similar performance as BLAS-like optimized routines. Moreover, with the increasing parallelism and heterogeneity as well as the ever increasing data-communication costs, numerical libraries often require to be modified or redesigned in order to take advantage of new features in parallel architectures. In our study we consider the dense linear algebra libraries LAPACK (serial library for CPU processors) and MAGMA (for Graphics Processing Units). The disparity between these libraries that target different architectures illustrate one of the issues in designing optimized linear algebra software. While being able to maintain a similar interface for the routines, the code and structure of all algorithms ported from LAPACK to MAGMA has to be rewritten to match the architectural features and the programming language of the accelerators. Furthermore, these libraries are mostly implemented using low-level languages like C or FORTRAN and thus cannot provide a high-level interface that would be closer to the specification language of the numerical linear algebra practitioner without losing performance. This issue is represented by the abstraction/efficiency tradeoff problem where raising the abstraction level with object-oriented and generic programming techniques is obtained at the cost of performance (see Figure 2.1). However, performance inhibits the flexibility and adaptability of libraries. Some solutions have been proposed in recent years but they tend to solve partially the abstraction/efficiency trade-off problem. The method followed by the Formal Linear Algebra Methods Environment (FLAME) with the Libflame library [START_REF] Van Zee | The libflame library for dense matrix computations[END_REF] is a good example. It offers a framework to develop dense linear solvers using algorithmic skeletons and an API which is more user-friendly than LAPACK, giving satisfactory performance results.

Another method is the one used by code generation projects like Spiral [START_REF] Püschel | Spiral: A generator for platform-adapted libraries of signal processing alogorithms[END_REF] for signal processing or other linear algebra libraries such as Design by transformation [START_REF] Marker | Designing linear algebra algorithms by transformation: Mechanizing the expert developer[END_REF], Build to Order [START_REF] Siek | Build to order linear algebra kernels[END_REF], or Hydra [START_REF] Lewis | Hydra: a scalable proteomic search engine which utilizes the hadoop distributed computing framework[END_REF] that develop code generation through the use of a Domain Specific Language (or DSL) to express data dependencies. A more closely related work is a linear algebra compiler [START_REF] Fabregat-Traver | A domain-specific compiler for linear algebra operations[END_REF]. It is based on the Mathematica language and optimizes the algorithms by reusing variables and mapping to BLAS function calls using problem specific knowledge [START_REF] Fabregat-Traver | Application-tailored linear algebra algorithms: A search-based approach[END_REF]. By definition, a DSL is a computer language specialized for a particular application domain. This implies that the use of a DSL requires a preprocessor and a custom compiler or interpreter. Furthermore, only the DSL compiler is aware of the underlying compiler existence limiting the integration of the DSL with other components such as an Integrated development environment (IDE) [START_REF] Fowler | Language workbenches: The killer-app for domain specific languages[END_REF]. Since writing a compiler is very complex and time-consuming while often not being re-usable, we do not wish to do such a task.

A more generic approach is the one followed in recent years by C++ libraries built around expression templates or other generative programming [START_REF] Czarnecki | Generative programming[END_REF] principles to design a Domain Specific Embedded Language. DSEL are languages implemented inside another host language. Designing a DSEL is easier than a DSL as it reuses existing compilers and relies on domain dependent analysis to generate efficient code (see Section 1.4.3). Problems such as copy elision and return value optimization that are implemented by compilers are often not exploited by complex DSLs like MAT-LAB which incurs copy penalties and slows down algorithms. Examples of such libraries are Armadillo [START_REF] Conrad | Armadillo: An open source C++ linear algebra library for fast prototyping and computationally intensive experiments[END_REF] and MTL [START_REF] Gottschling | Representation-transparent matrix algorithms with scalable performance[END_REF]. Armadillo provides good performance with BLAS and LAPACK bindings and an API close to MATLAB for simplicity. However it does not provide a generic solver like the MATLAB routine linsolve that can analyze the matrix type and choose the correct routine to call from the LAPACK library. It also does not support GPU computations which are becoming mandatory for medium to large dense linear algebra problems. In a similar way, while MTL can topple the performance of vendor-tuned codes, it does neither linsolves-like implementation nor GPU support. Other examples of such libraries include Eigen [START_REF] Guennebaud | A C++ Linear Algebra Library[END_REF], Flens [START_REF] Lehn | FLENS[END_REF], Ublas [START_REF] Walter | The boost uBLAS library[END_REF] and Blaze [START_REF] Iglberger | Expression templates revisited: a performance analysis of current methodologies[END_REF].

Our objective in this chapter is to provide a solution to the problems of portability and adaptability on new computer architectures. To this end, we propose a hybrid solver for CPU and GPU architectures with a single interface to solve dense linear systems. Our solution is designed on top of N T 2 , an open-source scientific library written in C++ available at www.github.com/NumScale/nt2. NT 2 provides a MATLAB-inspired API and its implementation is based on a meta-programming technique known as "expression templates" [START_REF] Veldhuizen | Expression templates[END_REF].

Generative programming for designing numerical libraries

Optimization approaches based on a configuration space

As stated in Section 2.1, developing complex linear algebra software is a non trivial task due to the large amount of both algorithmic and architectural requirements. These combined factors create a configuration space containing the various configurations available for a given system. Choosing the correct combination of factors from a configuration space will then ensure optimal performance.

Compiler techniques based on iterative compilation [START_REF] Triantafyllis | Compiler optimization-space exploration[END_REF], where several optimizations from a configuration space are tested and the best one is selected, is a classical technique to improve performance.

An example of these methods can be found in the ATLAS [START_REF] Whaley | Automatically tuned linear algebra software[END_REF] library which is based on using optimized binaries. Each function's binary is generated during the installation phase with the iterative compilation technique. The generation process is accelerated by a hierarchical tuning system. In this system, the lower level functions are subject to a large selection process ensuring their optimal performance. High-level functions like BLAS 3 routines can then exploit feedback from the previous steps of the configuration process.

A second method is based on a performance analysis at runtime. For instance, a system like StarPU [START_REF] Augonnet | StarPU: a unified platform for task scheduling on heterogeneous multicore architectures[END_REF] uses a monitored runtime system in which the performance of each function on a given hardware configuration is monitored in real-time. This monitoring allows StarPU to select the most optimized version of the algorithm by changing its parameters (tiling size, number of iterations,. . . ) or the targeted architecture (CPU, GPU, hybrid).

Both methods described above are valid approaches in the field of high performance computing. However, in our case, we aim at providing a library level system for such exploration [START_REF] Veldhuizen | Active libraries: Rethinking the roles of compilers and libraries[END_REF] that will complement the compilers work. One way to do this is to use generative programming.

Generative programming in software development

Generative programming (see Section 1.5.1) consists of bringing the benefits of automation to software development. Following this paradigm, a model can be drawn to implement the different components [START_REF] Czarnecki | Components and generative programming[END_REF] of a system. It is then possible to build a generator that will combine these components based on a generative domain model. This generator (or configuration knowledge) will ensure the transition from a configuration space with domain-specific concepts and features to a solution space that encapsulates expertise at the source-code level. The code generation process will be hidden from the end-user by various meta-programming techniques which turn the user interface into a simple and clean API, where few to none details about the algorithms and structures are visible.

Template meta-programming is a classical generative programming technique in which templates are used by a compiler to generate temporary source code. It is then merged with the rest of the source code and finally compiled. The output of these templates includes compile-time constants, data structures, and functions. The use of templates can be thought of as a compile-time execution that enables us to implement domain-specific optimizations.

Domain engineering methods for active libraries

We call active libraries [START_REF] Czarnecki | Generative programming and active libraries[END_REF] a technique which combines a set of generative programming and meta-programming methods to solve the abstraction/efficiency trade-off problem mentioned in Section 2.1. The main idea is to perform high-level optimization based on a semantic analysis of the code before any real compilation process. Such informations and transformations are then carried on by a meta-language that allows the developer to embed meta-informations in the source code itself, helping compilers to generate a better code by using these semantic informations. Active libraries are often implemented as DSEL.

Czarnecki proposed a methodol called Domain Engineering Method for Reusable Algorithmic Libraries [START_REF] Czarnecki | Generative Programming: Methods, Tools, and Applications[END_REF] (or DEMRAL depicted in the blocks 1-3 of Figure 2.2) based on the techniques described previously. DEMRAL is a DSEL-based method where domain specific descriptors are used to represent the various states of the system (represented as Block 1 in Figure 2.2). The various combinations of descriptors represent all the possible configurations in the system. In N T 2 , these parameters are available at the API level for the user [START_REF] Falcou | Software Abstractions for Parallel Architectures[END_REF]. Once these configurations have been implemented, it is possible to program the parametric components that they represent (see Block 2 in Figure 2.2). In N T 2 , these would correspond to the various skeletons available and various kernels for the CPU. The final step is to build a generator that will take the various descriptors as parameters, choose the corresponding component and generate the concrete application at compile-time based on this (see Block 3 in Figure 2.2). In N T 2 , this corresponds to the solver we have implemented which will be described later on. DEMRAL can be seen as a specialization of a paradigm like object-oriented programming, aspect programming or model driven engineering [START_REF] Schmidt | Guest editor's introduction: Model-driven engineering[END_REF]. While we can find a large number of algorithms (N ) and implementations for distinct data structures (P ), the problem is that combining them can result in a large number of code to write (N * P ). Using DEMRAL, only N generic algorithms and P data structure descriptions are needed since the generator will choose the correct Chapter 2. High-level programming for dense linear algebra domain-specific implementation from the configuration space with the help of the configuration knowledge.

The DEMRAL method provides a high re-usability, allowing components to be customized while retaining the efficiency of statically configured code [START_REF] Musser | STL tutorial and reference guide: C++ programming with the standard template library[END_REF]. We extend it by adding an architectural layer in the design with a "Domain specific architecture description" (block 4 of Figure 2.2) and a specialized generator for GPUs ( block 5 of Figure 2.2) based on this description [START_REF] Estérie | Multi-Architectural Support: A Generic and Generative Approach[END_REF]. In NT 2 , this would represent an extension to the API with a GPU tag and the addition of GPU based skeletons and kernels based on the MAGMA library. This enables us to have a separate specific generator for accelerators that will create a generic component with the appropriate marker that can then be combined with the already existing ones.

Building a CPU/GPU programming model

GPU container concept

Integrating GPU support in a library that only supports CPUs requires to decide at which level and how this integration should take place. In this work we will focus on CUDA integration to exploit the MAGMA library. This can simplify our task if we focus on calling existing functions that have already been compiled by NVCC. To details this, a major problem with the CUDA environment is the need of a special compiler for GPU code. This means a library compiled with a standard C++ compiler will not be able to seamlessly support CUDA code as it will require a new compilation step. It is also not possible to compile modern C++ code with NVCC as the compiler does not support properly the standard. In the situation where we have access to a binary, it is possible to directly call the function from a code compiled in C++.

Our approach will then focus on creating a minimalist GPU vector that will encapsulate the pointer and basic functionalities. This will enable us to have a recognizable structure for GPU data while having no need to pass through NVCC. The Thrust library [START_REF] Bell | Thrust: A productivity-oriented library for CUDA[END_REF] uses an opposite model where it allows random data access on a GPU vector but requires the code the be in .cu file and compiled with NVCC. In our implementation, we declare an element access operator [ ] for conformity reasons but it is not usable (will make a static assert).

Integration in NT 2

If the single system computation model of NT 2 is rather classic, we needed a way to handle accelerators in a generic and extensible way. The first step is to provide a simple way to locate tables (see Section 1.5.2) on either the host system or on a device. This is simply done by using a couple of settings in the definition of NT 2 table instances. In the Section 2.3.1 we described a low level vector. Tables are high level entities that include a container with settings on the locality (which we added), allocator type and many more. Listing 2.1 shows the nt2::host_ and device_ settings that are used to specify if a table contains data stored on the host memory or device memory. The device flag will make it so the table uses our GPU vector as a back-end with its restrictions. Listing 2.1: NT 2 host and device specifications Note that a special function on_device can be used to specify on which device the memory must be allocated in the case where multiple devices are available.

Semantic of operations between host and device tables is quite straightforward as they will be carried on the proper memory segment of each table. When mixing tables of different location, memory transfers are implicitly performed. This means that assigning a host table to a device table is equivalent to performing a CUDA memory transfer. This can be used for example to simplify interaction with existing GPU kernels as shown in listing 2.2. As streams are not assigned to tables, this transfer will be synchronous. A copy function is also available to perform asynchronous memory transfers when a non-default stream is given.

// X i s a 1 e3 x 1 e3 m a t r i x f u l l o f 1 . Listing 2.2: NT 2 interaction with cuBLAS This semantic of transfer by assignment is a classical way of performing such operation transparently. It as been used by tools like Thrust or VexCL [START_REF] Demidov | Programming CUDA and OpenCL: a case study using modern C++ libraries[END_REF] and have been proved to be easy enough for the user while allowing for fine grain performance tuning.

The current GPU programming model often requires to use pinned memory on the CPU for optimal transfer time and computation overlap. For this reason, we have added a CUDA pinned allocator that will use the CUDA CPU allocation function instead of a standard CPU alloc (Listing 2.3).

// X i s a 1 e3 x 1 e3 m a t r i x f u l l o f 1 a l l o c a t e d i n CUDA p in ne d memory . To enable support for IGP, we have also extended this model with new functionalities. As there is no simple way as of now to detect if the system is an IGP or not, we have added a compile flag (NT2_CUDA_INTEGRATED=ON/OFF, OFF is by default). This will activate specific parts to the CUDA integrated part such as flags definitions ( Listing 2.4) using the const object setup idiom 1 . In case the pinned allocator is used for the CPU, it will automatically call the cudaHostAlloc function with the correct flag cudaHostAllocMapped. As on IGP we do not want to do transfers between CPU and GPU using the cudaMemCpy function, we have to provide functions that can satisfy these conditions. We have designed two functions called to_host and to_device that provide these functionalities. The function to_device will verify if the data of X is pinned and if the flag NT2_CUDA_INTEGRATED is activated (Listing 2.5). If that is the case, the function will return a view of the entry table using its semantic information. If one of these conditions are not met, the type returned will be a table on the GPU and data from X will automatically be transfered. The to_host function will also check if the architecture is an IGP and if the output table is allocated as pinned. If these conditions are met, it will apply no modifications to Y. Otherwise, data will be transfered from Y on the GPU to X on the CPU. It is important to note that for IGP the function to_host holds a double purpose. Due to the IGP model, it is possible to run computations from the same memory space on GPU and CPU simultaneously. This means that if the user is not careful enough he can compute on the same data concurrently removing consistency between results. The function to_host will therefore apply a synchronization to make sure the computation is finished. Another solution is to use the CUDA_LAUNCH_BLOCKING=ON flag.

Application to linear algebra solvers

In this section, we describe our approach to automatically generate linear algebra solvers on parallel architectures. Our solution stems from the programming techniques combined with a proper configuration space and smart containers for data management on GPU.

Linear system solvers

The first step to build linsolve is to identify the key properties of the configuration space and the proper way to represent them. Once this analysis is done, we can refine these properties into high-level abstractions that will parameterize linsolve. These abstractions will then be used to define the configuration knowledge necessary to ensure the transition to the solution space. These properties represent the informations necessary to dispatch on the various solvers that we can find in the numerical libraries LAPACK and MAGMA. Concerning dense linear systems, we can identify three main properties that need to be taken into account: matrix structure, condition number, and targeted architecture.

A matrix structure can be divided into subcategories that can be identified statically (data type, storage scheme, matrix type and storage format). The data type and storage scheme parameters are already identified through the problem domain, respectively being scalar entries (real, double or single/double complex) and a dense matrix. The storage format is defined by N T 2 and shares a common interface with FORTRAN77 (column-major arrays). The matrix types correspond to the different ones available in the numerical libraries LAPACK and MAGMA (e.g., general, symmetric, hermitian...).

The second domain corresponds to the conditioning of the system. In current numerical libraries, the linear solvers are usually based on LU or QR factorizations in fixed precision, or mixed-precision algorithms with iterative refinement. It is not possible to identify statically if a system is ill-conditioned since it requires expensive computations which are not manageable at compile-time. Furthermore, it would be too costly to estimate the condition number at runtime for mixed-precision routines since it requires the factored form of the matrix (the LAPACK function gecon estimates the reciprocal condition number but requires the LU form bringing the cost to θ(n 3 ) for an n * n matrix). However, since current dense linear algebra libraries propose mixed precision routines, it needs to be part of the configuration space.

The last key domain of our solver is the dispatch between different architectures. As explained in Section 2.1, the architectural features of a GPU result in a very different language compared to a CPU. The solution we used to solve this abstraction problem is to provide through the use of a DSEL (Section 2.2.2) a common syntax between CPU and GPU routines. Using architecture aware binding, we can then freely decide whether to call LAPACK or MAGMA routines by dispatching on the different back-ends in an extensible way. It is now possible to define a grammar that encapsulates these ideas into a configuration space (see Section 2.2.1). These parameters are not mutually exclusive and can be extended/combined. Most of the parameters we can access are defined by the user and therefore configurable at the API level. In MATLAB, the linsolve routine does not take into account the data type, and the matrix type needs to be defined in a parameter structure containing the different matrix properties recognized (lower/upper triangular, upper Hessenberg, symmetric, positive definite, rectangular). While creating a matrix in N T 2 , the user has the possibility to define the matrix and data type which are optimized as meta-data properties of the matrix, using the following instruction: When calling linsolve, he will then have the possibility to give additional information on the conditioning of the matrix either as a parameter of the system: It is also possible to ask for complementary information as output like the reciprocal condition number returned by LAPACK : Once this is done, linsolve will be able to parse the configuration space by reading out the nested domain-specific features while assigning default values to the unspecified ones.

Memory management for hybrid computation

When using GPU-based systems, we need to ensure data consistency between the different physical memories. In this section, we discuss how these techniques are used jointly with linsolve. We can discern the two most common approaches to CPU and GPU containers. The first one is to statically define the locality of the container which is done in the Thrust library, while the second uses a dynamic approach like in SkePU [START_REF] Enmyren | SkePU: a multi-backend skeleton programming library for multi-GPU systems[END_REF].

The memory management mechanism in a dynamic approach allows to change the locality of a container and reallocate the data. The container then needs to manage the memory and ensure consistency between data and locality. In this situation, it is not possible to statically define the locality of a container. Therefore, our approach consists of adding an architectural tag similarly to the matrix type tag on our container (default locality is CPU). The purpose of this method is to enable the user to write programs using GPU functions in a transparent way.

It is then possible to ensure the transitions from CPU to GPU memory by using explicitly the tag. This does not prevent the decision-making process of the solver when no locality tag is given by the user. The solver can generate a GPU code performing data transfers from CPU to GPU as well as the reverse. The generation process will choose the architecture based on a combination of factors, mainly the matrix size and the algorithm. The GPU tag can also hold complementary informations passed as template parameters of the tag.

The definition of container locality being static, it is easier to define a data efficient memory management unit. Let's use the following scenario as an example :

x = nt2::linsolve(a,b);

From here, we can apply different strategies depending on the locality of x, a and b. In a situation where all three containers are on GPU (respectively CPU) memory there will be no locality problem as various data are located on the same device. However, the scenario where x is on the GPU (respectively CPU) while a and b are on the CPU (resp. GPU) will generate a conflict. The rules to solve locality conflicts are static and do not depend on the runtime. Therefore, the priority will be given to the locality of the result x to ensure consistency between the data and container locality.

Experiments were carried out on a system using 2 sockets of Intel Xeon E5645 2.40GHz and a Tesla C2075. We consider single precision random square matrices of size 2000, 10000 and 20000 and we solve a system of linear equations Ax = b using the LU factorization. The light grey bars in Figure 2.4 correspond to the following call made in N T 2 through linsolve that can run either on CPU or GPU. The dark grey bars correspond to C++ calls to either the LAPACK function sgesv or the MAGMA function magma_sgesv. These results show that automatically generated routines do not exhibit any overhead compared to direct calls to LAPACK or MAGMA. Note that the performance of all others generated routines does not incur any overhead as well.

Application: example of a least squares solver

In this section we illustrate how the generative programming method described in Section 2.2 can be used to generate automatically new implementations of algorithms and achieve satisfactory performance. We consider the overdetermined full rank linear least squares (LLS) problem

min x∈R n Ax -b 2 , with A ∈ R m×n , m ≥ n and b ∈ R m .
The most classical methods for solving linear least squares problems are based on the QR factorization or the normal equations. The latter method is twice cheaper (mn 2 vs 2mn 2 operations) but the error is then proportional to cond(A) 2 [20, p. 49]). However if A can be saved, we can also use the semi-normal equations (SNE) method where we solve the system

R T Rx = A T b,
where R is the triangular factor from the QR factorization of A (this is a straightforward reformulation of the normal equations). It is shown in [START_REF] Björck | Stability analysis of the method of semi-normal equations for least squares problems[END_REF] that, similarly to the normal equations method, the forward error bound involves a factor cond(A) 2 , even if we use a R-factor that is of better quality than the Cholesky factor because it has been computed via a backward stable algorithm. However, as explained in [20, p. 126 and p. 250], the accuracy of the SNE method can be improved by using the corrected semi-normal equations method (CSNE) that consists in adding one step of fixed precision iterative refinement to the SNE as follows:

1. Let x solve R T Rx = A T b 2. Compute r = b -Ax 3. Solve Aw = r 4. Corrected solution y = x + w
It is shown in [69, p. 392] that, if cond(A) 2 u ≤ 1 (u being the unit roundoff), then the forward error bound for the CSNE method is similar to that of a backward stable method (and even smaller when r = Axb is small). In that case, the CSNE method is more satisfactory than the SNE method but this is not true for all A. In the following we propose to use the CSNE method to solve LLS in mixed precision.

The efficiency of mixed precision algorithms has been proved on linear systems based on the LU factorization with results that can reach up to 90% [START_REF] Baboulin | Accelerating scientific computations with mixed precision algorithms[END_REF] of floating point computational rate in the lowest precision on current architectures. The method to solve mixed precision CSNE (or MCSNE) consists of first performing the factorization in single precision (ε s ) (if the matrix is not ill-conditioned) with the computational cost of θ(mn 2 ) and then refine the solution in double precision (ε d ) where operations cost θ(n 2 ). Iterative refinement [START_REF] Demmel | Error bounds from extra-precise iterative refinement[END_REF] is a method that produces a correction to the computed solution by iterating on it. Each k th iteration in this process consists of computing the residual r k = b -Ax k-1 , solving the new system Ad k = r k , and adding the correction x k+1 = x k + d k . Mixed precision iterative refinement will work as long as the condition number of the least squares problem [START_REF] Baboulin | Computing the conditioning of the components of a linear least-squares solution[END_REF] is smaller than the inverse of the lower precision used (i.e. here 10 8 ).

Algorithm 8 Mixed-Precision

CSNE Compute A = QR (ε s ) Solve R T x = A T b (ε s ) Solve Rx 0 = x (ε s ) do k = 1,2,... r k = b -Ax k-1 (ε d ) Solve R T x = r k (ε s ) Solve Rd k = x (ε s ) x k = x k-1 + d k (ε d ) check convergence
The first step of Algorithm 8 in NT 2 is implemented in Listing 2.6 with triu(qr) , while the second and third steps are performed by the following calls to mtimes and linsolve. Note that the code is similar in terms of syntax and number of instructions to what would be written in MATLAB. Once the solver for MCSNE is implemented using NT 2 , it becomes possible to add it to linsolve as a dispatch case of mixed precision solver for overdetermined linear systems. This would result in the following call :

x = nt2::linsolve(a,b,nt2::mixed_precision_);

Here in Listing 2.7 we show the corresponding code without using linsolve and implementing directly with MAGMA. It produces the same result as our code described in Listing 2.6 but is much more verbose and dependent on MAGMA functions. We still use our simple cuda_buffer to allocate simply data on the GPU and for automatic deallocation. Listing 2.7: Implementation of MCSNE with MAGMA Benchmarks were carried out using 2 sockets of Intel Xeon E5645 2.40GHz (peak Gflop/s is 230) and a Tesla C2075 (peak Gflop/s is 1030.4). We used Intel MKL [START_REF]Intel. Math Kernel Library (MKL)[END_REF] version 10.2.3, MAGMA 1.9 with CUDA 7.0 [START_REF]NVIDIA CUDA C Programming Guide[END_REF] and gcc 4.8 [START_REF] Gough | An Introduction to GCC[END_REF]. We compare ourselves with Eigen 3.2.5 and Armadillo 5.200.2. As a reminder, computations on the GPU are run using the IEEE 754 compliant mode for arithmetic operations and with CUDA ECC activated. The random test problems were generated using the method described in [START_REF] Paige | LSQR: An algorithm for sparse linear equations and sparse least squares[END_REF]. Performance results include data transfers between CPU and GPU and data.

. . . 1 0 0 0 0 / 1 0 0 0 1 3 0 0 0 / 3 0 0 0 1 5 0 0 0 / 5 0 0 0 1 8 0 0 0 / 8 0 0 0 2 0 0 0 0 / 1 0 0 0 0 In Figure 2.5 we compare the performance of CSNE_mixed_prec on the CPU with the LAPACK routine xgels that solves the LLS problem with a QR factorization without column pivoting. The results show that performance of MCSNE is only 10% less than that of sgels. Note that this represents around 75% of the peak performance of a matrix-matrix multiply in single precision (routine SGEMM). We use random matrices and the iterative refinement converged in less than 4 iterations.

If we compare our results with Armadillo, we can see that the implementation of MCSNE with Armadillo (Listing 2.8) is slower than dgels. While the performance of dgels in Armadillo is the same as in N T 2 , it is much harder to obtain high performance when not relying entirely on the BLAS. For examples, the extraction of the Q,R matrices are very expensive and the cast is not vectorized.However, as Armadillo aims toward a good balance between speed and ease of use, this is justified. The interface is intuitive and the documentation is well structured which makes writing code quite appreciable. Eigen (Listing 2.9) on the other side manages to reach better performance but is still slower than NT 2 . This is due to better memory management in NT 2 and some optimization functionalities that are not present in Eigen. The Eigen is however harder to use to implement a version of MCSNE and also less readable. Some of the code for MCSNE had to be hand-written as Eigen could not properly extract the upper triangular matrix R. Having an implementation where the solver is a member function of the matrix/solver makes the code look less readable and lose some of its generality and flexibility.

On the GPU, the performance of MCSNE in NT 2 is depicted in Figure 2.6. It also approaches 90% of the performance of magma_sgels on GPU (QR_single_prec) while being near twice faster than the routine in double precision.

The behavior of MCSNE when compared with QR solvers in double and single precision is similar to what was observed in [8, p. 15] for the LU factorization. Neither Armadillo nor Eigen offer proper support for GPU to allow us to write an MCSNE version.

Listing 2.9: CPU Implementation of MCSNE with Eigen

In Section 1.3.4, we described the concept behind the IGP architecture and current available architectures. We also note that current IGPs are not adapted for double precision. In Listing 2.7, we can see that the DGEMM routine on the NVIDIA Tegra X1 has very poor performance compared to the SGEMM. This is due double precision computations being emulated on the GPU.

An algorithm using mixed-precision will therefore be performing much better than in double precision since most computations are done in single precision. This is confirmed in Figure 2.7. The implementation of MCSNE that achieves good performance does not use directly MAGMA routines. The reason is that MAGMA does not have proper support for IGP and most computational routines are hybrid. MAGMA routines end up doing memory transfers between the CPU and GPU instead of using a zero-copy model. On IGPs, this is very costly, thus limiting the performance of the mixed-precision MAGMA routine for QR to 1 Gflop/s. It is also interesting to note that the MAGMA and CUBLAS implementations of SGEMM do not perform the same, with one being more efficient than the other depending on the size of the problem. This clearly shows that the size of the problem has an impact on how the code should be written for current IGPs.

Conclusion

Combining the large number of algorithms available in numerical libraries and architectural requirements in a generic solver for dense linear systems is a complex task. We showed that generative programming is a valid software development approach for addressing these issues while maintaining a high level of performance.

Our contribution furthers the work in active libraries by providing a viable way to make our software architecture-aware. Performance results illustrate that for both existing routines like those in linsolve and new ones such as MCSNE, the delivered performance is close to what state of the art libraries achieve.

The other interesting result is that software like NT 2 can quickly prototype new algorithms while providing support for various architectures. With NT 2 , we reach a good combination of high-level codes for linear algebra problems that gives good speedups and offers the users enough expressiveness to describe the problem in the most efficient way.

Future work includes support for more architectures like Intel Xeon Phi, with work on new algorithms that provide good performances while not being available in numerical libraries like randomized algorithms [START_REF] Baboulin | A Parallel Tiled Solver for Dense Symmetric Indefinite Systems on Multicore Architectures[END_REF][START_REF] Baboulin | Accelerating linear system solutions using randomization techniques[END_REF] or communication-avoiding algorithms [START_REF] Baboulin | A class of communication-avoiding algorithms for solving general dense linear systems on CPU/GPU parallel machines[END_REF] for dense linear systems. Moving to sparse problems is also a possibility where libraries like Cups [START_REF] Bell | Cusp: Generic parallel algorithms for sparse matrix and graph computations[END_REF] or VexCL provide an interesting approach. Raising the level of expressiveness stays a major concern while trying to add content in NT 2 .

In the next chapter, we will explain how we have build a multi-stage programming tool which is included in NT 2 for GPU operations. Developing large applications in a simple, fast and efficient way has always been an issue for software developers. As computing hardware complexity rose with the advent of SIMD, multi-processor, multi-core systems and more recently accelerators like GPUs [START_REF] Owens | A Survey of General-Purpose Computation on Graphics Hardware[END_REF] or Intel Xeon Phi [START_REF] Cramer | OpenMP Programming on Intel Xeon Phi Coprocessors: An Early Performance Comparison[END_REF], software design methods did not undergo the same amount of changes. This complicates the exploitation of such hardware in mainstream applications.

Designing DSL has been presented as a solution to these issues. As DSLs allow solutions to be expressed in their programming idiom with a level of abstraction equivalent to the problem domain, the maintainability and quality of code is increased. One of the most popular examples is MATLAB which provides a large selection of toolboxes that allow a direct expression of high-level algebraic and numerical constructs in a easy-to-use imperative language. In this scope, Domain Specific Embedded Languages (or DSELs ) [START_REF] Hudak | Building domain-specific embedded languages[END_REF][START_REF] Tratt | Model transformations and tool integration[END_REF] are languages implemented inside a general-purpose, host language [START_REF] Czarnecki | DSL Implementation in MetaOCaml, Template Haskell, and C++[END_REF]. without the requirement of a dedicated compiler or interpreter as they are designed as a library-like component.

In Chapter 2, we presented a method to take advantage of heterogeneous architectures using DSEL and providing support current HPC libraries. The problem of not being able to write simple CUDA kernels was however not addressed. Here, we wish to present an approach based on multi-stage programming that can circumvent these problems.

NT 2 -The Numerical Template Toolbox -is such a DSEL using C++ template meta-programming to provide a Matlab -inspired API while supporting a large selection of parallel architectures and keeping a high level of expressiveness. NT 2 was designed to support architectural features like SIMD extensions and multicore programming. However, the support for accelerators like GPGPUs was limited as GPU kernel compilers where unable to process NT 2 C++ 11 based implementation of C++ based DSEL . In chapter 2, we addressed the first part of the problem by providing a method for GPU based computing for dense linear algebra kernels.

In this chapter, we present a new method and extension for NT 2 that takes care of such accelerators, especially CUDA based GPUs through multi-stage programming [START_REF] Eckhardt | Implicitly heterogeneous multi-stage programming[END_REF] (or MSP ). MSP consists in doing multiple compilation phases allowing for type-safe program generation. The work presented in this chapter was submitted and accepted at ATMG 2016 [START_REF] Masliah | Meta-programming and multi-stage programming for gpgpus[END_REF]. Our contributions include:

• A multi-stage method to alleviate the host/device programming model with a cost model for offloading code

• An adaptable strategy to generate CUDA kernel directly from a single C++ source file containing NT 2 statements

• The integration of this kernel generator with existing CUDA libraries like cuBLAS or MAGMA.

The purpose of this system is to provide the user some leeway on how to distribute the data between the host and device through a simple mechanism. As having a perfect cost model for load balancing is very complex to put in place and costly, letting the user provide some insight on data locality is beneficial.

After reviewing the concurrent state of the art software libraries (Section 3.1), we then describe the kernel generator process and how it can be integrated with an existing library (Section 3.2). Finally, we present benchmarks assessing the generated code quality (Section 3.4) and conclude on the future work regarding NT 2 and accelerators.

Related Work

Software libraries for GPU computing try to simplify the new programming paradigm brought by many-core based systems. The general trend followed in recent years by C++ libraries is to provide a high-level interface through template meta-programming techniques. This goal is reached by providing device containers with architectureaware generic parallel algorithms and/or a code generation based process. We will give a detailed explanation of both type of libraries that support OpenCL/CUDA or both.

Thrust [START_REF] Hoberock | Thrust: A parallel template library[END_REF] is a header only library providing a similar interface to the C++ Standard Template Library. Its high-level interface is based on meta-programming and traits to provide efficient parallel skeletons that can run on either CPU or GPU. Container locality is expressed through an explicit container declaration limiting the abstraction but allowing for easier transfers between host and device. Locality for functions is defined by default for device vectors and can be extended with tag dispatching. Overall, it is a well rounded utility library which can be combined with CUDA Toolkit libraries such as CUBLAS, CUFFT and NPP. However it lacks of code generation features and does not support OpenCL.

C++AMP [START_REF] Gregory | C++ amp: accelerated massive parallelism with microsoft visual c++[END_REF] is a framework developed by Microsoft that includes a C++ library and an associated compiler. The C++ AMP programming model includes multidimensional arrays, indexing, memory transfer, tiling, and a mathematical function library. You can use C++ AMP language extensions to control how data is moved from the CPU to the GPU and back. It relies on parallel skeletons such as parallel_for to iterate on data. It does not provide any AST optimization or MSP and is more of an API based approach like Thrust with OpenCL.

SyCLs [START_REF]Khronos OpenCL Working Group -SYCL subgroup. Opencl sycl provisional specification[END_REF] aim is to simplify programming manycore processors by offering a modern C++ interface on top of OpenCL. SyCL is similar to C++AMP in the sense that it provides single-source development and inlined kernels. That means kernel code can be embedded within C++ code. They have added support for templates in kernel code which is not possible with OpenCL. It does not provide a model using the CUDA language.

VexCL [START_REF] Demidov | Programming CUDA and OpenCL: a case study using modern C++ libraries[END_REF] is an expression template library for OpenCL/CUDA. It provides a high level generic interface that is suitable for both back-ends with static parameters defined within a DSEL for linear algebra. The expression template mechanism allows for code generation by lazy evaluation of vectors and elementary operations within the AST. Similarly to Thrust, it provides STL-like functions on containers that have a defined locality. It is also possible for the user to define custom functions on device that will be dynamically generated for CUDA. However, the transform used for the generation process requires a unique data locality limiting hybrid algorithms.

ViennaCL [START_REF] Rupp | ViennaCL-a high level linear algebra library for GPUs and multi-core CPUs[END_REF] is also an expression template library for OpenCL/CUDA. This library strictly follows the uBLAS programming interface and STL like algorithms for easier integration with other softwares. It has implementations for BLAS kernels and high level solvers for sparse and dense computation that provide good performance. Through the mechanism of expression templates it can evaluate basic linear algebra operations with operator overloading. ViennaCL focuses more on OpenCL due to the necessity for separate compilation with CUDA limiting its support through the OpenCL language. It is however possible to generate CUDA code it does not have any real support for code generation or a device API yet.

The purpose of the previously described libraries is usually to provide a wrapper over the C++ language for GPGPU computing. For this reason, the DSELs based on expression templates or Boost.Proto are usually lightweight and consist mostly of overloading elementary operations for containers. The code generation phase then ends up doing a dynamic compilation of the CUDA code/OpenCL kernel which adds a significant overhead on small sized or low arithmetic intensity problems.

Furthermore, as mentioned previously, it is not possible to compile a NT 2 source code with NVCC even with the latest version (7.5). To address these issues, we have added a two step compilation in NT 2 using a serialization mechanism based on Boost.Serialization. Before describing this process, we will first detail our concept behind MSP.

Multi-stage programming in C++

As described in Section 3.1, code generation for accelerators in C++ is based on Meta-programming techniques. This process is however limited by the fact that C++ does not natively support language extensions for runtime code generation and program execution [START_REF] Taha | A gentle introduction to multi-stage programming[END_REF]. Creating a multi-stage paradigm is therefore necessary to remove the runtime cost of compiling CUDA code.

Designing a multi-stage software

We have described in Section 1.4.4 some of the core concepts behind MSP. MSP has been used extensively in fields based on stochastic problems such as energy planning [START_REF] Pereira | Multi-stage stochastic optimization applied to energy planning[END_REF] as it relies on two-stage problems. LMS described in Section 1.4.4 is conceptually similar to Meta-programming techniques in C++ applied to the Domain Engineering Method for Reusable Algorithmic Libraries [START_REF] Czarnecki | Generative programming -methods, tools and applications[END_REF] which we described in Section 2.2.3. Language visualization described also in Section 1.4.4 goes further by combining DSEL techniques with multi-stage programming to cover multiple DSEL (Probabilities, machine learning, Scripting ...) and combine them with heterogeneous hardware. We wish to stay within the same DSEL of linear algebra but use multi-stage to remove compatibility problems with accelerator based systems due to compiler incompatibilities. In Chapter 2, we have shown that by relying on compiled libraries such as MAGMA, we can ease the programming model by not needing NVCC. This however limits code development inside a library as it is not possible to directly write CUDA code in NT 2 . What we want to provide is model that supports using existing optimized libraries while being still able to extend them for more flexibility.

MSP is an ideal solution to such as problem as it allows type-safe program generation and is expendable. We differentiate ourselves from other standard MSP approaches by wanting a second step compilation. Most MSP approaches add language extensions to a program which is compiled before the actual program itself. What we wish is a MSP layer that will safely generate accelerator code based on our algorithm skeletons for CPU. We want to perform code generation by providing CMake macro to generate all the required calls to our kernel generation system. CMake is a cross-platform tool for managing the build process of software using a compiler independent method. It has been widely adopted in recently year instead of Makefiles as it provide more abstraction and interoperability. The CMake layer is very lightweight and can be replaced by any generation tool or even simple scripts.

Multi-stage programming in C++

The way we use MSP is not based on language extensions but on a method we developed as explained in Section 3.2.1. It is possible to create a multi-stage compilation of a program using this representation. In C++ , we base this process on doing a compilation phase with an incomplete link to generate only an object file. This is depicted in Phases 1-2 of Figure 3.2. Generating an executable program using a compiler is a multi-stage process divided into two main components: compilation and linking.

By definition, compilation refers to the processing of source code files to create object files. This process does not create an executable file for the user. This object file contains the compiled code in binary form of the symbols defined in the source code.

Linking refers to the creation of an executable or a library from an object file set. During the compilation, if the compiler cannot find the definition of a specific function in a file, it will assume that this function is defined in another file. As long as the source code is well-formed, the compiler will not complain. It is then to the linker to replace the references to undefined symbols with the correct addresses. Remaining undefined symbol are left untouched if not found and the linking process becomes incomplete. This will generate a new object file where remaining undefined symbol are left untouched.

It is then possible to use a demangling tool like cppfilt or nm to decode the C++ ABI names (steps 2-3 of Figure 3.2). Each demangled symbol will correspond to the internal C++ representation with a complete prototype of the incomplete function. By parsing this representation we can generate the CUDA/OpenCL kernel code with a corresponding host code to complete the link phase. Figure 3.2 describes this multi-stage process. To benefit from this paradigm, it is necessary to include a two stage compilation that cannot be deployed in a header only library. To implemented this process we use the CMake tool. It is also important to develop a readable intermediate representation that can be easily parsed when demangled. A specific tool called symbol/code converter (Figure 3.2, Part 3) coupled with a serialization process was developed as an add-on to NT 2 to solve the issues described previously. The serialization process to solve the readability of the abstract representation is based on Boost.Serialization . This representation is similar to the expression-template AST based on Boost.Proto that is parsed by the C++ compiler. This enables us to get the semantic information of the container described in chapter 2 like its data locality, data type or matrix shape. It is however not possible to have run-time informations like the container size. The role of symbol/code converter is two-fold. First, it must parse the demangled symbols obtained. This is done with Boost.Spirit [START_REF] De Guzman | The boost spirit parser generator framework[END_REF], a C++ library to parse expressions and generate outputs based on them. It is implemented as a DSEL using Expression templates and Meta-programming techniques. With this library, we can generate outputs corresponding to the semantic information of the containers and operators. Secondly, it must generate the device and host code (Figure 3.2, Part 4-5) from the output generated by Boost.Spirit. In order to achieve this, we must deserialize the abstract representation in which the semantic informations are represented by a tag.

We will first describe the generation of the CUDA kernel. Each representation obtained from the AST corresponds either to an elementary expression of containers such as a = b + c (or a = b + c + d . . . ), a sequence of operations if a fused operator (tie) is called or to an algorithmic skeleton (reduce, scan...). The parsing is separated between the left and right hand-side of the computation. The left hand-side will check if the expression is a terminal or an AST and generate the left part of the CUDA kernel expression. If it is a fused operator, it will generate a sequence of left operators in the kernel. The right hand-side will parse the operator and generate from the NT 2 CUDA back-end the corresponding operation with its parameters. Similarly, a fused operator will generate a sequence of right-hand side.

The generation of the host code consists in creating the source file corresponding to the functions with missing symbols. This amounts to adding the includes for the currents back-end, the CUDA kernel call and streaming data if necessary. As the device locality is available in the AST under the tag nt2 :: device_ (see Section 2.3), we can stream the data to the GPU only if needed. Chaining non-fused operations that are on the host will obviously trigger back and forth data transfers. The overhead generated by such a process is diminished by the streaming and multi-stream process enabled for each operation.

We just described the key concepts of the symbol/code converter tool implemented in NT 2 for multi-stage programming. In the next section we will give concrete examples of generated kernels from symbol/code converter in NT 2 and how it supports hybrid computation.

Software implementation for NT 2

To describe the generation process, we first use the Triad kernel for element-wise operations which consists in doing a fused multiply-add (or fma :

a = b + c * d ).
The resulting code in NT 2 is : During the compilation phase, the operation in Listing 3.1 is replaced with a call to the CUDA transform skeleton. This skeleton will then call the device kernel that is not implemented yet resulting in the incomplete link sequence(Part 1-2, Listing 3.1). The equivalent code for transform is detailed in Listing 3.2. This transform will analyze the informations on the input matrices and the architecture to decide if it should proceed with the generation process. We can then demangle the symbols resulting in the code described in Listing 3.4. This code corresponds to the Boost.Proto AST in NT 2 (see Figure 3.3)that we parse to generate the host and device code. The architectural tag of the machine is depicted in purple in Listing 3.4 and the fma computation node in red.

v o i d nt2 : : e x t e r n a l _ k e r n e l <nt2 : : t a g : : transform_ , nt2 : : t a g : : cuda_<b o o s t : : simd : : t a g : : avx_> > :: c a l l <nt2 : : c o n t a i n e r : : t a b l e <f l o a t , nt2 : s e t t i n g s ( ) >, nt2 : : c o n t a i n e r : : view<nt2 : : t a g : : table_ , f l o a t c o n s t , nt2 : : s e t t i n g s ( ) >, nt2 : : c o n t a i n e r : : view<nt2 : : t a g : : table_ , f l o a t c o n s t , nt2 : : s e t t i n g s ( )> >, 3 l >, nt2 : : memory : : c o n t a i n e r <nt2 : : t a g : : table_ , f l o a t , nt2 : : s e t t i n g s ( nt2 :

: of_size_ <-1l , -1l , -1l , -1l >)> > c o n s t &)
Listing 3.4: NT 2 CUDA Triad demangled

The function tag (nt2::tag::transform_) corresponds to each element-wise operations in NT 2 . Any function can be constructed using external kernel if necessary. As an example, algorithmic skeletons (scan, reduce, zip ...) will call external kernel when the CUDA back-end is activated. Then, we can generate a thrust code in a .cu file that will then be compiled with NVCC. Using this technique, we can overcome a weakness of the thrust library while benefiting from its optimized routines. The generated code for the .cu file is described in Listing 3.5 for Kepler cards. It corresponds to the AST representation with additional semantic information specific to the CUDA language. A function wrapper that calls the CUDA kernel (fma4_wrapper ) is used to separate the compilation of the .cu file with NVCC from the corresponding C++ host code. The triad kernel directly calls the fma function from CUDA as the NT 2 AST recognizes every occurrence of an fma and replaces it by its function (Listing 3.4, see boost::simd::tag::fma_). This optimization in itself is not essential since it can be done by NVCC but is still interesting as it demonstrate the potential of code generation. As NT 2 already optimizes the AST by replacing patterns with corresponding functions, we can benefit from this analysis. We can then call the corresponding CUDA function if available or our own implementation for each pattern. The host code is partially described in Listing 3.6 and 3.7. Listing 3.6 corresponds to the initialization of parameters and data before doing the actual computation on the device. The blockSize and stream number is determined during the generation process depending on the number of parameters and architecture. The blockSize is usually generated by measuring the bandwidth of transfers from host to device in a range and choosing the most optimal one.

__global__ v o i d fma4 ( f l o a t * t0 , c o n s t f l o a t * t1 , c o n s t f l o a t * t2 , c o n s t f l o a t * t 3 ) { i n t i d x = b l o c k I d x .

// A l l o c a t i n g memory on t h e d e v i c e v a l u e ( a0 ) . s p e c i f i c s ( ) . a l l o c a t e ( b l o c k S i z e , nStreams , s i z e , t r u e ) ; v a l u e ( c h i l d _ c <0>(a1 ) ) . s p e c i f i c s ( ) . a l l o c a t e ( b l o c k S i z e , nStreams , s i z e ) ; v a l u e ( c h i l d _ c <1>(a1 ) ) . s p e c i f i c s ( ) . a l l o c a t e ( b l o c k S i z e , nStreams , s i z e ) ; v a l u e ( c h i l d _ c <2>(a1 ) ) . s p e c i f i c s ( ) . a l l o c a t e ( b l o c k S i z e , nStreams , s i z e ) ;

// c h e c k s redundancy between i n p u t s and o u t p u t s s t d : : u no rd ere d_ set <c o n s t f l o a t * > addr ; addr . i n s e r t ( c h i l d _ c <0>(a0 ) . data ( ) ) ; f o r ( s t d : : s i z e _ t i =0; i < nStreams ; ++i ) { cudaStreamCreate (& stream [ i ] ) ; } Listing 3.6: NT 2 CUDA Triad Host Code initialization Using our tuning framework, we observed that to benefit best from the Kepler GPU bandwidth a block size of 40000 is necessary for single precision values. Double precision computations would lead to a lower block size generated. Similarly, the blockDim used is the highest value available for the architecture as we generate element-wise operations. Depending on the number of parameters the size may be lowered to limit the allocations. The allocation process includes pinned memory and device memory allocation. If containers were defined on the GPU with nt2 :: device_, they would not appear in the allocation phase. As we have no information on the pointer for each container, we use an unordered set to limit redundancy in memory transfers. Listing 3.7 describes the computation phase. It relies on block streaming with transfers to the GPU only if NT 2 tables are on the host. This streaming process is based on the overlap data transfers concept described by NVIDIA. It consists in creating multiple streams (the number depends on the architecture and problem intensity/size) and launching for each stream a transfer host to device, the CUDA kernel and the transfers device to host for a block. As the host memory has already been allocated, we must first transfer the data to pinned memory with cudaHostAlloc to benefit from GPU optimizations. Since the difference in bandwidth between pinned and page-able memory only increases with new architectures, this optimization can give a speedup even with a mono-stream program.

f o r ( s t d : : s i z e _ t i = 0 ; i < n ; ++i ) { s t d : : s i z e _ t j = i % nStreams ; v a l u e ( a0 ) . s p e c i f i c s ( ) . t r a n s f e r _ h t d ( a0 , i , stream [ j ] , j ) ; v a l u e ( c h i l d _ c <0>(a1 ) ) . s p e c i f i c s ( ) . t r a n s f e r _ h t d ( c h i l d _ c <0>(a1 ) , i , stream [ j ] , j , addr ) ; v a l u e ( c h i l d _ c <1>(a1 ) ) . s p e c i f i c s ( ) . t r a n s f e r _ h t d ( c h i l d _ c <1>(a1 ) , i , stream [ j ] , j , addr ) ; v a l u e ( c h i l d _ c <2>(a1 ) ) . s p e c i f i c s ( ) . t r a n s f e r _ h t d ( c h i l d _ c <2>(a1 ) , i , stream [ j ] , j , addr ) ; fma4_wrapper ( v a l u e ( a0 ) . s p e c i f i c s ( ) . data ( j ) , v a l u e ( c h i l d _ c <0>(a1 ) ) . s p e c i f i c s ( ) . data ( j ) , v a l u e ( c h i l d _ c <1>(a1 ) ) . s p e c i f i c s ( ) . data ( j ) , v a l u e ( c h i l d _ c <2>(a1 ) ) . s p e c i f i c s ( ) . data ( j ) , dimGrid , blockDim , stream As stated in Section 2.3, computations on the device only occur if some conditions are met. As of now, these conditions are limited to the problem size and data locality but can be extended as the call to transform is automatic when NT 2 has defined that CUDA is available. Due to the hierarchical tag dispatching in NT 2 , a system with an Intel processor coupled with an NVIDIA card will have a tag similar to the following : cuda_ < openmp_ < simd_extension >> . Therefore, if conditions for dispatch on the GPU are not met we will call the next level of transform ( i.e.openmp). This enables us to use both the CPU and GPU in parallel depending on the problem which is a functionality rarely implemented in libraries.

To describe the generation process for algorithmic skeletons based on Thrust we will use a reduction. The resulting code in NT 2 is : Reductions are simpler to generate as they only require to call the correct underlying parallel skeleton and take care of memory transfers. The evolution of GPU architectures has an impact on the underlying implementation of the algorithm and its performance. Therefore, we prefer to rely on Thrust which is a standard NVIDIA library instead of implementing each parallel skeleton. It is simpler to use Thrust as a back-end in our MSP process rather than implementing complex algorithms such as fold,scan... and optimize them for each architecture. The demangled symbols corresponding to Listing 3.8 can be seen in Listing 3.9. The demangled code is similar to the one described previously in Listing 3.4. The function tag is changed to nt2 :: tag :: global_ which corresponds to the reduction skeleton. From there, it is easy to generate the corresponding host and kernel code. The kernel code described in Listing 3.10 is a wrapper for the thrust reduction function. There is no overhead compared to calling the function directly with Thrust. • GCC 4.9, CUDA 7.5

// D e f i n

Black & Scholes kernel

The Black & Scholes algorithm represents a mathematical model that gives a theoretical estimate of the price of European call and put options on a non-dividendpaying stock. It is a bandwidth bound algorithm for GPU if we take into account the memory transfers. Black & Scholes requires an important number of registers and high latency operations such as log or exp. The Black & scholes algorithm is a fusion of multiple transform skeletons requiring large entry sizes to obtain suitable performance.

The code is given in Listing 3.13 using the loop-fused technique described previously with the operator tie. The nt2 :: device_ tag is specific for accelerator enabled architectures. However, if we use the tag while no accelerator is available we will fall back to the default architecture which is nt2 :: host_. t i e ( d , d1 , d2 , r ) = t i e ( s q r t ( Ta ) , l o g ( S/X) +(fma ( s q r ( v ) , 0 . 5 f , r ) * Ta ) / ( v * d ) , fma(-v , d , d1 ) , S * normcdf ( d1 )-X * exp(-r * Ta ) * normcdf ( d2 ) ) ; r e t u r n r ; } Listing 3.13: NT 2 black and scholes This additional semantic information on memory locality can help to avoid useless memory transfers while staying simple enough for the user. This will result in the .cu file in Listing 3.14 generated for floating point values. Since the AST does not contain the name of the parameters, the kernel generator has to give a different name to each one. Even with duplicated parameters, the memory transfers are only done once and a single pointer can be passed multiple times. Memory allocation on device and data transfers between host and device do pointer checking in the host code (see triad example) to insure no redundant work is done incurring also negligible overhead. The fnms function is due to an NT 2 AST transformation and corresponds to the fused negated multiply-subtract of three values.

__global__ v o i d bs ( f l o a t * t 0 , f l o a t * t 1 , f l o a t * t 2 , f l o a t *
[ i ] = s q r t f ( t 4 [ i ] ) ; t 1 [ i ] = p l u s ( l o g f ( d i v i d e s ( t 5 [ i ] , t 6 [ i ] ) ) , d i v i d e s ( m u l t i p l i e s ( t7 , t 8 [ i ] ) , m u l t i p l i e s ( t9 , t 1 0 [ i ] ) ) ) ; t 2 [ i ] = fnms ( t11 , t 1 2 [ i ] , t 1 3 [ i ] ) ; t 3 [ i ] = fnms ( m u l t i p l i e s ( t 1 4 [ i ] , e x p f ( m u l t i p l i e s ( t15 , t 1 6 [ i ] ) ) ) , f a s t n o r m c d f ( t 1 7 [ i ] ) , m u l t i p l i e s ( t 1 8 [ i ] , f a s t n o r m c d f ( t 1 9 [ i ] )
) ) ; } Listing 3.14: NT 2 black and scholes fused CUDA kernel The Black & Scholes algorithm involves high latency and high register count operations. This will result in sub-optimal performance on SIMD for the CPU due to spilled registers while a Kepler GPU will not have any such problem. As seen in Figure 3.4, the execution time of the kernel on the GPU is negligible (2-3 ms) compared to the overall time of the optimized version with SIMD and OPENMP in NT This is due to the highly parallel nature of the GPU. As its computational power is much higher, simple elementwise operations will naturally fit the programming model while CPUs must be optimized. Most of the time is spent transferring the data between host and device memory which can be avoided with the right semantic information available ( nt2 :: device_). In the scenario presented in Listing 3.13, which represents the worse case, we still have better performance than the Boost.Compute version on which most C++ libraries are based as we use block streaming with pinned memory reaching a near-optimal throughput (average of 9.7 GB/s) for device transfers. As the bandwidth of a memcpy on the CPU (5 GB/s) is faster than page-able transfers (3GB/s) even without any overlap between transfers and computation we still increase the performance. This optimization can be disabled depending on the CUDA architecture. We have also compared with an optimized CUDA implementation and the performance was identical. This shows that our approach to generate element-wise operations does not incur any overhead on the kernel computation time.

As computations on the GPU are often performed in great number, if the user allocates the data on the GPU he will pay no extra transfer cost for the rest of the computations and in this case the CUDA kernel is up to twelve times faster than the hand-optimized version for CPU. Also, using the flag nt2 :: pinned_ if tables must stay on CPU allows us to gain up to a x3 factor in transfers time. The optimized Black & Scholes version we implemented in CUDA reached the same performance in kernel and transfers time when using the same type of memory allocation for data. On the Tegra X1, the GPU is much more efficient than the CPU. This is caused by the large difference in performance between the Maxwell GPU and the 4-core ARM CPU. The ARM cortex A57 has around 64 Gflops while the GPU achieves 500 Gflops. The code in Listing 3.15 corresponds to the one executed on the Tegra for the GPU. As the no-copy model from NVIDIA for IGPs require data to be pinned, we have to allocate CPU tables with the pinned_ flag. t i e ( d , d1 , d2 , r ) = t i e ( s q r t ( Ta ) , l o g ( S/X) +(fma ( s q r ( v ) , 0 . 5 f , r ) * Ta ) / ( v * d ) , fma(-v , d , d1 ) , S * normcdf ( d1 )-X * exp(-r * Ta ) * normcdf ( d2 ) ) ; r e t u r n r ; } Listing 3.15: NT 2 Black and Scholes for Tegra

The important performance difference between the K40m and the Tegra is due to the on device bandwidth difference between the two cards. On our system, the K40m has an on device bandwidth of 200 GB/s while the Tegra has around 30GB/s.

Linsolve kernel

Linsolve is the Matlab solver for linear systems. As NT 2 has its own implementation of linsolve with a LAPACK/MAGMA back-end, we can combine it with the code generation process. As an example, we consider the solution of a dense linear system using the LU factorization and apply one step of iterative refinement:

1. Compute r = b -Ax.
2. Solve Ad = r.

Update y

= x + d.
The equivalent NT 2 code is the following : t a b l e <T, device_> A, b , x , e ; t a b l e <T, s e t t i n g s ( device_ , u p p e r _ t r i a n g u l a r _ )> u ; t a b l e <T, s e t t i n g s ( device_ , l o w e r _ t r i a n g u l a r _ )> l ; If executed on the CPU, the code in Listing 3.16 will call the LAPACK routines. The semantic information upper_triangular_ allows linsolve to call the triangular solver instead of doing the classic linear solve. If executed on the GPU, the same optimizations will be applied and the iterative refinement process will trigger calls to transform for both element-wise operations.

The performance results in Figure 3.6 attest that the performance obtained with our model is relevant. The GPU version with MSP calls magma kernels using the CUBLAS DGEMM routine without doing any transfer and reaches near peak performance of a K40m GPU which corresponds to 1.40 Tflop/s. The version that does not use MSP is slower as transfers are done during the iterative refinement step. The CPU version quickly reaches the peak performance of both CPU which is 210 Gflop/s. As we can see, there is no performance loss while call the LAPACK/-MAGMA back-ends and if device pointers are passed to our code generator, there will be no memory transfers. Similar performance would also be reached using the other factorizations available in NT 2 .

Conclusion

The development of tools for simplifying accelerator programming has been an active topic since accelerators have become a mainstream element of high-performance computing systems. In this chapter, we proposed an extension of a high-level, dataparallel scientific computing library to specifically handle GPU accelerators. Our main objectives were to keep a similar level of expressiveness in the client code while supporting different use cases of accelerator programming.

To reach this goal, we have implemented a multi-stage system in which the initial C++ code is used to automatically generate the equivalent CUDA kernel by reusing our internal representation of this code. This representation is based on C++ Expression Templates and the Algorithmic Skeleton to identify and classify expressions based on the kind of loop nest that is required. Finally, we have illustrated on a selection of examples that the performance obtained is close to the hardware capability and exhibits benefits compared to other solutions.

Work is still on-going on this system, including support for more specific functions on the latest GPUs. Implementing a more thorough cost model to ensure better scheduling of computation between CPU and GPU is also being studied. The natural evolution of this work is to extend our approach to the Intel Xeon Phi coprocessor, the runtime-based CUDA compiler recently released, or OpenCL devices. An interesting track of research can be derived from the support of OpenCL by targeting OpenCL enabled FPGAs, as we could bridge the gap between highlevel C++ and hardware design.

In recent years, tensors have started to take an important place in HPC applications. This is due to tensors' frequent use in physics and engineering, where tensors provide a mathematical tool for brief and comprehensive formulations and solutions of problems in areas such as elasticity, fluid mechanics, multi-physics, quantum chemistry, general relativity, and many others [START_REF] Kolda | Tensor decompositions and applications[END_REF]. Tensors are multi-dimensional arrays that can be used to describe physical properties featuring multilinear relations. Well known mathematical objects like scalars, vectors, and matrices can be generalized to tensors that are of order zero, one, and two, respectively. Also, tensor transformations like flattening of a tensor to matrices or reshaping of matrices into tensors, can be used to link tensor computations to the developments in high-performance numerical linear algebra. Therefore, similar to many applications, tensor computations can also significantly benefit from representing their computations in terms of BLAS, as well as from other dense LA algorithms and techniques for multicore and GPU architectures.

This work is the result of a collaboration with Azzam Haidar, Stanimire Tomov, Ahmad Abdelfattah and Jack Dongarra from the Innovative Computing Laboratory in Tennessee. Parts of the work presented in this chapter have been accepted at ICCS 2016 [START_REF] Abdelfattah | High-performance tensor contractions for GPUs[END_REF] and Europar 2016 [START_REF] Masliah | High-performance matrix-matrix multiplications of very small matrices[END_REF]. As this project is meant to be part of the MAGMA library, we decided not to use NT2 to limit external dependencies. However, the techniques presented here stems from the same principles and have been integrated in NT 2 . The contributions of this chapter are the following :

• A method to take advantage of modern C++ design for hybrid computations and interface development

• A strategy using performance analysis based on hardware features and counters to implement efficient algorithms no specific architecture

• An implementation of a small batched matrix-matrix product that can perform better than state of the art libraries

We start by giving in Section 4.1 a small overview of some Tensor libraries for numerical algebra. Then, in Section 4.2 we explain the challenges we face while developing a tensor interface and how to exploit modern C++ (post C++11). The next Section 4.4 explains how we can optimize batched computations for small matrix-matrix products which are often found in tensor computations. We finish this chapter by giving concluding remarks in Section 4.5.

Tensors in numerical libraries

The use of tensors in physics and engineering applications has motivated an extensive research in both algorithms for tensor algebra and computational aspects for tensor-based simulations. See the survey [START_REF] Kolda | Tensor decompositions and applications[END_REF], and the references cited there, for an overview of the linear algebra aspects of tensors research, including higherorder tensor decompositions, their applications, and available software. A general overview with current and future direction in tensor research is presented in the Future directions in tensor-based computation and modeling workshop report [START_REF]Future directions in tensor-based computation and modeling[END_REF].

There has been an important number of software packages developed for tensors. Some are designed to be used through numerical computing mathematical environments like the Tensor Toolbox 1 for Matlab, GRTensor II 2 for Maple, or Ricci 3 for Mathematica. These packages tend to focus on providing proper tensor computations and do not provide highly optimized operations or target accelerators. Others contractions such as Tensor Transpositions Compiler [START_REF] Springer | Ttc: A high-performance compiler for tensor transpositions[END_REF] (TTC), an open-source parallel compiler for multidimensional tensor transpositions, focus on specific optimizations to generate efficient code. Tensor packages for specific applications such as quantum chemical computations [START_REF] Wenjing | Acceleration of Streamed Tensor Contraction Expressions on GPGPU-Based Clusters[END_REF] present ongoing work on using accelerators to accelerate tensor contractions on GPUs. The approach uses code generation techniques and is incorporated in the NW Chem computational chemistry suite. More recent work [START_REF] Nelson | Generating efficient tensor contractions for gpus[END_REF] also uses code generation techniques and autotuning (with a system called Baracuda, based on CUDA-CHiLL and a high-level module Optimizing Compiler with Tensor OPeration Intelligence (OCTOPI)) to report significant acceleration compared to NW Chem on particular tensor contractions.

While we use code generation and auto-tuning, our approach also relies on context knowledge, and in particular that tensor reshaping techniques can often transform tensor contractions from many applications into many GEMMs, right from machine learning. This transformation can not easily be detect automatically, and moreover, even if the transformation is somehow indicated, it is well known that general GEMM optimizations using only compiler techniques can not match in performance best-tailored solutions. For distributed memory system that are CPU only, [START_REF] Solomonik | A massively parallel tensor contraction framework for coupled-cluster computations[END_REF] executes contractions via GEMM on a properly ordered and structured tensor. As in the other approaches in quantum chemistry, large tensor contractions are targeted. In contrast, we target many but small contractions, that are often very small size (see next section), resulting in large-scale computations. Tensor reshuffle operations are done to cast the contractions to GEMMs, when possible, and a batched approach with custom-built BLAS kernels is used to efficiently execute them. Additional closely related efforts include [START_REF] Markidis | Openacc acceleration of the nek5000 spectral element code[END_REF][START_REF] Shin | Speeding up nek5000 with autotuning and specialization[END_REF] and [START_REF] Stock | Model-driven simd code generation for a multi-resolution tensor kernel[END_REF].

Applications and challenges

The evolution and diversification of computer hardwares and programming languages is transforming the balance of computer systems. This evolution is producing changes at every level from programming efficiently the hardware to developing high level programming tools. From the point of view of numerical libraries, and the important number of applications that depend on them, three challenges stand out:

• the need to exploit unprecedented amounts of parallelism;

• the need to maximize the use of data locality and vectorized operations; and

• the need to cope with component heterogeneity.

Below, we highlight our main contributions related to the algorithm's design and optimization strategies aimed at addressing these challenges on multicore CPU and GPUs.

Exploit Parallelism and Vector Instructions: Clock frequencies are expected to stay constant, or even decrease to conserve power; consequently, as we already see, the primary method of increasing computational capability of a chip will be to dramatically increase the number of processing units (cores), which in turn will require an increase of orders of magnitude in the amount of concurrency that routines must be able to utilize as well as increasing the computational capabilities of the floating point units by extending it to the classical Streaming SIMD Extensions set (SSE-1, to SSE-4) in the earlier 2000, and recently to Advanced Vector Extensions (AVX, AVX-2, AVX-512). We developed specific optimization techniques that demonstrate how to use the many cores (currently multi-socket 10 -22 cores for the Haswell CPU or a K40 GPU using MAGMA to get optimal performance. The techniques and kernels developed in this work are generic and can be used elsewhere.

Hierarchical Communication Techniques that Maximizes the use of Data Locality:

Recent reports (e.g., [START_REF] Fuller | Committee on Sustaining Growth in Computing Performance[END_REF]) have made it clear that time per flop, memory bandwidth, and communication latency are all improving, but at exponentially different rates. So computation on very small matrices, that can be considered as computation-bound on old processors, is, -today and in the future-communicationbound and depends from the communication between levels of the memory hierarchy. We demonstrate that, performance is indeed harder to get on new manycore architectures unless hierarchical communications and optimized memory management are considered in the design. We show that, only after we developed multilevel memory design, our implementations reach optimal performance.

Performance Analysis and Auto-tuning:

We demonstrate the theoretical maximal performance bounds that could be reached for computation on very small matrices. We studied various instructions and performance counters, as well as proposed a template design with different tunable parameters in order to evaluate the effectiveness of our implementation and optimize it to reach the theoretical limit.

Container concept

The design of our code is done using new features of C++ for better re-usability and adaptability of the code. By using advanced template techniques we can create high-level interfaces without adding any cost even for small matrix-matrix products. To do so, we have designed a batched structure which will contain a C++ vector for the data and static dimensions. By using the C++ constexpr keyword and integral constants we can make a generic batched code that will dispatch at compile time contractions the correct version depending on the size of matrices. We use this environment for each code sequence we generate.

To develop high-quality HPC software for tensor contractions, we impose the following three main requirements on our interface design:

Convenience of use: Interfaces of HPC libraries are extremely important for the libraries' adoption by the community. Interfaces must provide convenience and practicality of use, including ease of interoperability between libraries. Ideally, a tensor data type standard must be defined. The standard for a dense matrix is a pointer, sizes, leading dimension, and assumption for columnmajor data layout, e.g., as in BLAS and LAPACK. For tensors, motivated by reviewing interfaces in available libraries and the success of BLAS, we propose to represent a tensor by its dimensions and a data layout abstraction. The abstraction is to provide a choice of predefined layouts, and ways to switch it, or to easily add user-specified layouts, without changing the underlying algorithms. In general, a specific layout provides the formula of mapping the multi-dimensional tensor to the memory, e.g., a second order tensor can be stored as a column-major matrix A with leading dimension lda, in which case the abstraction maps A i,j to A[i + j * lda].

Readability: Numerical software must be understandable, which is needed for ease of maintenance, as well as code optimizations, and testing. While we can easily implement any interface, e.g., even expressing the interface and tensor APIs in a DSEL if needed (plus code generation techniques to translate the DSEL to a standard language; such as the Einstein notation used in tensor computations), we determined that it is better to provide implementations relying on a standard language and the code generation features provided within the language. The C++14 standard for example is expressive enough to allow us to implement readable and easy to use interfaces.

Performance: While we expect to obtain high performance mostly through algorithmic design and auto-tuning (see Section 4.4), we do not want to compromise on optimization opportunities like removing parameters checking and loop unrolling to eliminate jumps and loop count decrements. These optimizations are essential especially for the small computations that we target. Therefore, in our design we consider the use of compiler features related to code generation (e.g., templates, etc.), as further discussed below.

C++11/14 features

The development of programming languages and their use has dramatically changed in recent years leading to continuous evolution. C++ is an example of such a programming language. The standardization committee has decided to make a new standard every 3 years, with the next release being the C++17 standard. Compatibility with the C language has also starting to break with the introduction of the auto keyword. The cause of these changes is the need for higher level language that provide better idioms for generic and generative programming. Here we will discuss new features of the C++11/14 standard that we use to develop a high level tensor interface.

The first feature of the C+11 language that we will discuss is auto [START_REF] Järvi | Decltype and auto (revision 3)[END_REF]. Consider the following declaration in Listing 4.1 :

// x i s t h e t y p e o f 7 : i n t auto x = 7 ; Listing 4.1: C++ auto Here x will have the type int because it is the type of its initializer. In general, we can write // x i s t h e t y p e o f e x p r e s s i o n auto x = e x p r e s s i o n ; Listing 4.2: C++ generic auto and x will be of the type from the value expression in Listing 4.2. For any variable, auto specifies that the type of the variable that is being declared will be automatically deduced from its initializer. This allows to write high level complex code without having the burden of complex types that can appear. We can apply the auto keyword on several features of the C++ language. As an example, in Listing 4.3 we can iterate with a range for statement on a vector of a generic type without defining the type of the value x. Another important feature of the C++11 standard is the constexpr [START_REF] Reis | General constant expressions for system programming languages[END_REF] keyword. constexpr provides a mechanism that can guarantee that an initialization is done at compile time. It also allows constant expressions involving user-defined types.

In Listing 4.4, the fibonnaci function is guaranteed to be executed at compile time if the value passed x is available at compile time. An important addition to the C++ standard for extensible design is variadic templates [START_REF] Gregor | Variadic templates for c++[END_REF]. In computer programming, variadic templates are defined as templates that take a variable number of arguments. Prior to C++11, templates contractions could only take a fixed number of arguments which were specified at the declaration of the template. The function add will accept any number of arguments and will compile as long as the operator+ can be applied between the arguments. This is type-safe as the verification is done at compile time. The resolution will be done using C++'s template and overloading resolution rules. Putting the [...] before the parameter name, typename... Args is called a template parameter pack. Using the [...] after a definition Args... args is called function parameter pack. This is how the function add is defined : the first argument is peeled off from the template parameter pack into type T. With each subsequent call, the parameter pack gets shorter by one parameter until the last occurrence with only one parameter in the pack is attained.

// sum i s o f t y p e i n t auto sum = add ( 1 , 2 , 3 , 6 ) ; // sum1 i s o f t y p e s t r i n g s t d : : s t r i n g x = " s t r i n g " , y = " add " , z = " f u n c t i o n " ; auto sum1 = add ( x , y , z ) ; Using constexpr and the features described previously also allow for integral constants. Integral constants are part of the C++ standard and wrap a static constant of a specific type in a class. We show in Listing 4.9 a possible definition for them 4 . Listing 4.9: Possible integral constant definition These types of expression can be used as non-type template arguments, array sizes, and in other contexts that require constant expressions. They are often used in conjunction with meta-programming techniques. The most common library pre-C++11 that was based around integral constant is Boost.MPL5 . Integral constants can be declared in the following way : integral_constant < int, 25 > ic_25, with ic_25 being the type for the value 25. We can the use integral constants in conjunction with user literals. User literals allow integer, floating-point, character, and string literals to produce objects of a wanted type by defining a suffix. Listing 4.10 is the code to declare a user defined integral constant literal that we can then pass as function parameter. The keyword using is also a new feature of the C++11 standard that replaces typedef and allow for template parameters. t e m p l a t e <u n s i g n e d l o n g N> u s i n g int_ = s t d : : i n t e g r a l _ c o n s t a n t <u n s i g n e d l o n g , N>; c o n s t e x p r i n t t o _ i n t ( c h a r c ) { r e t u r n s t a t i c _ c a s t <i n t >(c ) - User defined literals give us the means to simply dispatch on these values or use them as static function parameters instead of using templates. This allows for a more natural interface for the user and a simpler code generation process. Using Listing 4.10, we can now just use 25_c to represent the value 25. It will be automatically converted into an integral constant. Since C++11, operator"" is used for user-defined literals allowing for specific implementations. The easiest way to parse such a literal taking a variadic template char argument is by unpacking it, hence the const char array as a parameter to the parse function.

The last important feature that we will talk about is lambda functions [START_REF] Järvi | C++ lambda expressions and closures[END_REF]. A lambda expression is a mechanism for specifying a function object. The primary use is to specify a simple action to be performed by some function. For example: s t d : : v e c t o r <i n t > in , in 1 , out ; s t d : : t r a n s f o r m ( i n . b e g i n ( ) , i n . end ( ) , i n 1 . b e g i n ( ) , out . b e g i n ( )

, [ ] ( i n t a , i n t b )->i n t { r e t u r n a + b ; } ) ; Listing 4.11: C++ lambda using std transform std :: transf orm will iterate on each value of in and in1 and apply the lambda function passed on each element. The resulting operation will be written in the out vector. The lambda starts with the [ ] identifier, called the capture specification. It serves to specify how local variable are passed (by value : = , by reference : &). Lambdas are similar to functors: the variables captured by reference are like private pointers, and the ones captured by value are like private variables.

Lambdas are lightweight, nameless functions that can be defined where they are used. Lambdas are now a part of the C++ standard template library and current and future algorithms such as <thread> or <future>, allow you to specify function arguments as lambdas beside C-like functions and functors. Parallel skeletons implemented in the C++ standard library and Thrust have added support lambda functions. Modern C++ libraries also support lambdas as seen in Section 3.1. An example for thrust is shown in Listing 4.12.

t h r u s t : : h o s t _ v e c t o r <i n t > in , in 1 , out ; t h r u s t : : t r a n s f o r m ( i n . b e g i n ( ) , i n . end ( ) , i n 1 . b e g i n ( ) , out . end

, [ ] ( i n t a , i n t b )->i n t { r e t u r n a + b ; } ) ; Listing 4.12: C++ lambda using thrust transform Since CUDA 7.5, NVIDIA has also added experimental support for device lambda using the _device_ flag. A proposal to use ranges6 for the C++ standard library instead of iterators is also in study.

Developing a tensor interface

To develop a high level tensor interface, we use the C++ features detailed previously while keeping in mind the challenges described in Section 4.2. Related to performance, a lost optimization opportunity is if static checking and compile time information is not provided as part of the scientific libraries used. As an example, LAPACK routines start by checking entry parameters dynamically, which results in extra time for small size computations. The type of algorithms that we intend to use, also require specific tuning whether it be on CPU or GPU with MAGMA batched [START_REF] Haidar | Batched matrix computations on hardware accelerators based on GPUs[END_REF] as performance will greatly vary depending on numerous parameters.

To avoid these performance drawbacks, and also benefit from optimization techniques like compiler loop unrolling for static dimension problems, we use features of the new C++14 standard. In particular, the constexpr specifier enables to evaluate the value of a function or variable at compile time. We use this feature with integral constants and template specialization to design our tensor dimensions layout: As seen in Listing 4.13, we propose several ways to specify dimensions using the function of_size which returns a layout containing the static or dynamic dimension. Each operator inside the layout uses the constexpr keyword which enables us to return sizes at compile time if possible. Using constexpr will give the compiler more optimization opportunities such as loop unrolling.

The layout is class based on variadic templates to have a parameterizable number of dimensions. At the initialization of the layout, parameters are automatically unpacked, giving the value for the nth dimension if it is available at compile time. Dynamically sized dimensions will get an initial value of -1 to say they have not been initialized. To enable this optimization, the number of dimensions in a layout is decided at compile time leaving only the size of each dimension to be changeable. In Listing 4.14, we give an example of the definition for a simple layout to better view the key concepts behind it. The DGEMM batched function can be dispatched depending on the locality, calling either our implementation for small matrix product with static sizes on CPU or the MAGMA/CUBLAS implementation if the tensor is defined on the GPU.

Optimizing small matrix-matrix products

Performance measures

To evaluate the efficiency of our algorithms we derive theoretical bounds for the maximum achievable performance P max = F/T min , where F is the number of operations needed by the computation and T min is the fastest time to solution. For simplicity, consider C = αAB + βC on square matrices of size n. We have F ≈ 2n 3 and T min = min T (T Read(A,B,C) + T Compute(C) + T W rite(C) ). Note that we have to read/write 4n 2 elements, or 32n 2 Bytes for double precision (DP) calculations. Thus, if the maximum achievable bandwidth is B (in Bytes/second), and we assume T Compute(C) → 0 for very small computation, then T min = T Read(A,B,C) + T W rite(C) = 4n 2 /B in DP. Note that this time is theoretically achievable if the computation totally overlaps the data transfer and does not disrupt the maximum rate B of read/write to the GPU memory. Thus,

P max = 2n 3 B 32n 2 = nB 16 in DP.
The achievable bandwidth can be obtained by benchmarks. For our measures, we used the STREAM benchmark [START_REF] Mccalpin | Memory bandwidth and machine balance in current high performance computers[END_REF] and the Intel memory latency checker 3.0 tool for CPU with gcc 5.2, and the NVIDIA's bandwidthTest for GPU. Our tests show that contractions the practical CPU bandwidth we are able to achieve using different benchmarks is about 44 GB/s per socket. On the K40 GPU with ECC on the peak is 180 GB/s, so in that case P max is 2.75 n GFlop/s per socket for the CPU and 11.25 n GFlop/s for the K40 GPU. Thus, when n = 16 for example, we expect a theoretical maximum performance of 180 GFlop/s in DP on the K40 GPU.

The implementation of a matrix-matrix products kernel for very small matrices for CPUs requires specific design and optimizations. As we can store three double precision matrices of size up to 32 × 32 in the L1 cache of an Intel Xeon E5-2650 v3 processor, one can expect that any implementation will not suffer from data cache misses. This can be seen on Figure 4.6 where the performance of an ijk implementation, which is not cache-aware and cannot be vectorized, is pretty close to the ikj one. For smaller sizes, the ijk implementation is even more efficient than the ikj one, as it optimizes the number of stores (Figure 4.2). To obtain a near optimal performance, we conduct an extensive study over the performance counters using the PAPI [START_REF] Weaver | Measuring energy and power with PAPI[END_REF] tools. Our analysis concludes that in order to achieve an efficient execution for such computation, we need to maximize the occupancy and minimize the data traffic while respecting the underlying hierarchical memory design while maximizing the use of SIMD registers. Unfortunately, today's compilers cannot introduce highly sophisticated cache/register/SIMD based loop transformations and, consequently, this kind of optimization effort should be studied and implemented by the developer [START_REF] Loshin | Efficient Memory Programming[END_REF]. This includes techniques like reordering the data so that it can be easily vectorized, reducing the number of instructions so that the processor spends less time in decoding them, prefetching the data that will be reused in registers, and using an optimal blocking strategy.

Data Access Optimizations and Loop Transformation Techniques

In our design, we propose to order the iterations of the nested loops in such a way that we increase locality and expose more parallelism for vectorization. We can then use AVX2 SIMD instructions to optimize to the fullest our code. The matrix-matrix product is an example of perfectly nested loops which means that all the assignment statements are in the innermost loop. Hence, loop unrolling, loop peeling, and loop interchange can be useful techniques for such algorithm [START_REF] Ahmed | Tiling imperfectly-nested loop nests[END_REF][START_REF] Bacon | Compiler transformations for high-performance computing[END_REF]. These transformations improve the locality and help to reduce the stride of an array based computation. In our approach, we propose to unroll the two innermost loops so that the accesses to matrix B are independent from the loop order, which also allows us to reorder the computations for continuous access and improved vectorization. This technique enables us to prefetch and hold some of the data of B into the SIMD registers. Here, we manage to take advantage from the knowledge of the algorithm, and based on the principle of locality of references [START_REF] Hennessy | Computer Architecture, Fifth Edition: A Quantitative Approach[END_REF], to optimize both the temporal and spatial data locality.

Register Data Reuse and Locality

Similarly to the blocking strategies for better cache reuse in numerically intensive operations (e.g., large matrix-matrix products), we focus on register blocking to increase the performance. Our study concludes that the register reuse ends up being the key factor for performance. The idea is that when data is loaded into SIMD register, it will be reused as much as possible before its replacement by new data. The amount of data that can be kept into registers becomes an important contractions tuning parameter. For example, an 8 × 8 matrix requires 16 256-bit AVX-2 registers to be completely loaded. As the targeted hardware consists of only 16 256-bit AVX-2 registers, one can expect that loading the whole B will not be optimal as we will have to reload the vectors for A and C. However, if we load only 8 registers for B, which is equal to 4 rows, we can compute a row of C at each iteration and reuse these 8 registers for each iteration. We propose an auto-tuning process to check all the possible scenarios and provide the best option. This reduces the number of load, store, and total instructions from O(n Each C ends up being loaded/stored once. We apply this strategy to matrix sizes ranging from 8 to 32 as for smaller sizes the whole matrix can fit in registers. Different blocking strategies (square versus rectangular) have been studied through Store Mi0, Mi1, Mi2, Mi3 14: end for our auto-tuning process in order to achieve the best performance. We generate each matrix-matrix product function at compile time with C++ templates. The matrix size is passed as a function parameter using C++ integral constants.

As explained in Section 1.4, using SIMD units is quite simple if the algorithm has been well designed. As we have taken into consideration from the start of the algorithmic design SIMD problems, we can see that Algorithm 9 is very similar to the code in Listing 4.18. The main differences come from the syntactic sugar of the Intel SIMD instruction set. The code written here can only be efficient using the constexpr keyword. In-be saturated by a few cores. We performed a set of experiments towards clarifying this behavior and illustrate our findings in Figure 4.4. As shown, the notion of perfect speed-up does not exist for a memory-bound algorithm, and adding more cores increases the performance slightly. We performed a bandwidth evaluation when varying the number of cores to find that a single core can achieve about 18 GB/s while 6 and 8 cores (over the available 10 cores) can reach about 88% and 93% of the practical peak bandwidth, which is about 44 GB/s. We also studied NUMA-socket (non-uniform memory access) [START_REF] Hager | Introduction to High Performance Computing for Scientists and Engineers[END_REF] when using two Xeon sockets as seen in Figure 4.5. A standard memory allocation puts all of the data in the memory slot associated to the first socket until it gets filled, then starts filling the second socket. Since the problem size we are targeting is very small, most of the data is allocated on one socket, and thus using extra 10 cores of the second socket will not increase the performance. This is due to the fact that the data required by the cores of the second socket goes through the memory bus of the first socket, and thus is limited by the bandwidth of one socket (44 GB/s). There are ways to overcome this issue. By using NUMA with the interleave=all option, which spreads the allocation over the two sockets by memory pages, we can improve the overall performance. However, for very small sizes, we observe that such solution remains far from the optimal bound since data is spread out over the memory of the two sockets without any rules that cores from socket 0 should only access data on socket 0, and vice versa. To further improve performance, we use a specific NUMA memory allocation, which allows us to allocate half of the matrices on each socket. As shown in

Application to ARM processor

The ARM processor that we use for this benchmark is the CPU of the Tegra X1, a 4-core Cortex A57. The problematics details earlier still apply to the Tegra but on a different scale. Indeed, the ARM intrinsics only support 128-bit vectors which severely limit the SIMD use for double precision computations. In Figure 4.7, we compare the performance between our generated code, an ijk code, an ikj code and OpenBLAS [START_REF] Xianyi | Openblas[END_REF] using the latest version available from the develop branch on Github (28/6/2016).

Results follow the same trend we saw on the Intel processors. On very small sizes, ijk and ikj versions are quite efficient as the arithmetic intensity is very low, limiting the usefulness of parallelism. With increased sizes, we start to see these version stale and reach a limit set around 3.5 Gflops. The OpenBLAS version provides good performance but is limited by its blocking model which is not adapted for very small sizes.

Efficient wrapper for GPU batched GEMM

As we have seen in Section 4.3.2, our tensor class supports GPU data. In Figure 4.8, we show the performance comparison between the MAGMA DGEMM batched routine, the read-only cache version (rocache) in black and the fully optimized version in red, and CUBLAS using CUDA 7.5. The most important factors to optimize GPU performance is to maximize the occupancy and the efficiency of global memory reads. The first one is the ratio between the number of active warps per active cycles and the maximum number of warps that can run on an SM. The second is defined as the ratio between the load throughput requested by the kernel, and the actual required throughput needed to fulfill the kernel load requests. As MAGMA manages to reach better performance than CUBLAS, we use the MAGMA routine for GPU batched. contractions

Conclusion

In this chapter, we provide a method to develop high level interface using modern programming techniques in C++. We also give a technique using architectural tools and performance analysis to implement efficient algorithms for a specific architecture. Furthermore, we show how to develop a flexible interface that can take advantage of optimized libraries such as MAGMA or CUBLAS with minimal effort. This work can easily be extended because the interface is generic. The use of variadic templates allow for simple extensions of the container as additional parameters can be added and parsed without modifying current code.

Our work is motivated by a large number of applications, ranging from machine learning to big data analytics, that require fast linear algebra on many independent problems that are of size 32 and smaller. For these applications, the use of batched GEMM for small matrices is essential to achieve performance.

Conclusion and Perspectives

Conclusion

Throughout this manuscript, we studied the challenges behind designing high-level software and optimized codes for modern architectures. It is currently one of the main challenge faced to develop high performance software. The new Top 500 7ranking of June 2016 has shown the appearance of the Chinese machine from the National Supercomputing Center in Wuxi. This machine is now ranked first in the Top 500 and has a new architecture and applications must be designed specifically for it [START_REF] Fu | The sunway taihulight supercomputer: system and applications[END_REF] 8 ). The software stack for this machine includes compiler support for C, C++ and Fortran. Currently, C++ is the most common high level language used for high performance computing. While algorithms must be reimplemented for this machine, a library can just be extended.

In the first part of this thesis, we presented a method of high level programming that takes into account the features of heterogeneous architectures and the properties of matrices to build a generic dense linear algebra solver. Our programming model supports both implicit or explicit data transfers for GPUs and IGPs. Combining the large number of algorithms available in numerical libraries and architectural requirements in a generic solver for dense linear systems is a complex task. We showed that using generative programming is a valid software development approach for addressing these issues while maintaining a high level of performance. We focused on supporting programming models for heterogeneous architectures and we addressed this issue. The resulting software has been integrated into a C++ library called NT 2 which is available on Github.

The second contribution of our thesis is to provide a multi-stage system that can alleviate interoperability problems between the CPU and GPU programming models. Our multistage approach is used to automatically generate GPU code for CPU-based element-wise expressions and parallel skeletons while allowing for typesafe program generation. This work has also been integrated into the NT 2 library and used conjointly with previous work. In this context, we provide an adaptable strategy to generate CUDA kernels directly from a single C++ source file containing NT 2 statements. We illustrated our approach on a selection of examples and the performance obtained is close to the hardware capability and exhibits benefits compared to other solutions. Finally, we showed that this approach can be extended to new architectures or implemented using other programming languages.

93

Our last contribution is to apply high level programming techniques to batched computations and tensor contractions using modern C++. Using the experience we gained from our previous work, we were able to show that by combining a high level programming approach and advanced parallel programming techniques, we can outperform state of the art numerical libraries without relying on current libraries. To do so, we studied the issues around batched computations, memory locality and code vectorization to implement a highly optimized matrix-matrix product for small sizes using SIMD instructions. To obtain a near optimal performance, we conducted an extensive study over the performance counters using the PAPI tool and the Intel memory checker V3.0 for memory bandwidth.

Perspectives

As future work, we are interested in developing support for the OpenCL 2.0 standard in both our solver and multi-stage programming model. As CUDA is not generic enough and works only for NVIDIA cards, it is a limiting factor to develop a multiplatform library. An implementation of MSP in NT 2 has been proposed for OpenCL but it is still not complete. The recent Tegra X1 card also has limited support for OpenCL which is a problem. As recent library have added a layer on top of OpenCL for C++, it can become interesting to study and use such an approach.

We have also started looking at applying similar high level programming techniques to sparse problems. In Chapter 4, we presented a first approach to tensors and an application to batched GEMMs. It would be interesting to keep digging in this direction and try to support more features and benchmark applications.

With the incoming new architectures such as the Xeon Phi Knight's Landing or the future NVIDIA GPUs with NVLINK and IBM processors, current software will need to be extended. Indeed, these new architectures will bring changes in the programming models and load-balancing strategies. Proving that a high level generic approach can be a safe and efficient way to program these architectures will be a challenge to face. Titre : Méthodes de génération automatique de code appliquées à l'algèbre linéaire numérique dans le calcul haute performance Title : Automatic code generation methods applied to numerical linear algebra in high performance computing Keywords : Generic programming, Generative programming, C++, linear algebra, GPU Abstract : Parallelism in today's computer architectures is ubiquitous whether it be in supercomputers, workstations or on portable devices such as smartphones. Exploiting efficiently these systems for a specific application requires a multidisciplinary effort that concerns Domain Specific Languages (DSL), code generation and optimization techniques as well as application-specific numerical algorithms. In this PhD thesis, we present a method of high level programming that takes into account the features of heterogeneous architectures and the properties of matrices to build a generic dense linear algebra solver. As GPUs have become an asset in high performance computing, incorporating their use in general solvers is an important issue. We extend our approach to a new multis-tage programming model that alleviates the interoperability problems between the CPU and GPU programming models. Our multistage approach is used to automatically generate GPU code for CPUbased element-wise expressions and parallel skeletons while allowing for type-safe program generation. Finally, we investigate how to apply high level programming techniques to batched computations and tensor contractions. We first explain how to design a simple data container using modern C++-14 programming techniques. Then, we study the issues around batched computations, memory locality and code vectorization to implement a highly optimized matrix-matrix product for small sizes using SIMD instructions.

1. 2 . 11 Algorithm 2

 2112 Solving dense linear systems LU factorization in place, partial pivoting Input: A is a n × n matrix 1: for k = 1 : n do

Algorithm 4 2 : 3 : 4 :

 4234 Solving Ax = b with LU factorization 1: Compute the factorization P A = LU with P permutation matrix, L lower triangular and U upper triangular. Solve Ly = P b. Solve U x = y. Solution is x.

Algorithm 5 1 : 2 : 3 : 4 :

 51234 Solving Ax = b with QR factorization Compute the factorization A = QR with Q orthogonal matrix and R upper triangular matrix. Form d = Q T b. Solve Rx = d by backward substitution. Solution is x.

1. 3 .

 3 Architectures : Design and evolution 15 1.3 Architectures : Design and evolution 1.3.

Figure 1 .

 1 Figure 1.1: SIMD concept

Figure 1 . 2 :

 12 Figure 1.2: Multicore concept

Figure 1 .

 1 Figure 1.3: NUMA concept

Figure 1 . 4 :

 14 Figure 1.4: GPU Kepler based architecture

Figure 1 . 5 :

 15 Figure 1.5: Work flow for a GPU

  Figure 1.6: IGP architecture

Figure 1 . 7 :

 17 Figure 1.7: Work flow for a zero copy IGP

  __m256 a , b , c ; c = _mm_mul_ps( b , c ) ; Listing 1.1: Intel AVX2 single precision multiplication

Listing 1 . 2 :

 12 Intel AVX2 single precision multiplication no intrinsic

Chapter 1 .Figure 1 . 8 :

 118 Figure 1.8: Intel's 3rd generation IGP architecture

Chapter 1 .Figure 1 . 9 :

 119 Figure 1.9: Expression Templates in NT 2

Figure 2 . 1 :

 21 Figure 2.1: Abstraction/Efficiency trade-off for languages

Figure 2 . 2 :

 22 Figure 2.2: Overview of the NT 2 skeleton based generation process

  // G e n e r a t e s a h o s t t a b l e by d e f a u l t t a b l e <double> A( o f _ s i z e ( 1 e3 , 1 e3 ) ) ; // G e n e r a t e s a h o s t t a b l e e x p l i c i t l y t a b l e <double , host_> A2( o f _ s i z e ( 1 e3 , 1 e3 ) ) ; // G e n e r a t e s a d e v i c e t a b l e t a b l e <double , device_> D( o f _ s i z e ( 1 e3 , 1 e3 ) ) ; // G e n e r a t e s a h o s t p i n e e d t a b l e with CUDA t a b l e <double , pinned_> D2( o f _ s i z e ( 1 e3 , 1 e3 ) ) ;

  t a b l e <double> X = o n e s ( 1 e3 , 1 e3 ) ; // T r a n s f e r t o d e v i c e t a b l e <double , device_> Y = X; // cuBLAS d i r e c t c a l l c u b l a s D s c a l ( Y. s i z e ( ) , 5 . , Y. data ( ) , 1 . ) ; // T r a n s f e r back t o h o s t X = Y;

  t a b l e <double , pinned_> X = o n e s ( 1 e3 , 1 e3 ) ; // F a s t e r t r a n s f e r t o t h e d e v i c e t a b l e <double , device_> Y = X; // cuBLAS d i r e c t c a l l c u b l a s D s c a l ( Y. s i z e ( ) , 5 . , Y. data ( ) , 1 . ) ; // F a s t e r T r a n s f e r back t o h o s t X = Y; Listing 2.3: NT 2 cuda CPU allocator

  #i f d e f i n e d (NT2_CUDA_INTEGRATED) s t r u c t i n t e g r a t e d _ c u d a _ s e t t i n g s { i n t e g r a t e d _ c u d a _ s e t t i n g s ( ) { c u d a S e t D e v i c e F l a g s ( cudaDeviceMapHost ) ; } } ; s t a t i c c o n s t i n t e g r a t e d _ c u d a _ s e t t i n g s p r o p e r _ s e t t i n g s= { } ; #e n d i f Listing 2.4: Initialize CUDA for IGP

  // X i s a 1 e3 x 1 e3 m a t r i x f u l l o f 1 a l l o c a t e d i n CUDA p in ne d memory . t a b l e <double , pinned_> X = o n e s ( 1 e3 , 1 e3 ) ; // W i l l c r e a t e a view o f X on IGPs auto Y = t o _ d e v i c e (X) ; // cuBLAS d i r e c t c a l l c u b l a s D s c a l ( Y. s i z e ( ) , 5 . , Y. data ( ) , 1 . ) ; // W i l l not do t h e t r a n s f e r on IGPs to_host (Y,X) ; Listing 2.5: NT 2 IGP functions

  nt2::table<double,nt2::symmetric_> a;

x

  = nt2::linsolve(a,b,nt2::ill_conditioned_); or of the matrix: x = nt2::linsolve(nt2::ill_conditioned_(a),b);

  nt2::tie(x,r) = nt2::linsolve(a,b);

Figure 3 . 7 .Figure 2 . 3 :

 3723 Figure 3.2 shows the various steps to perform a call to linsolve with a symmetric matrix (here of size 5000) in N T 2 linked with the MAGMA library. The user starts by defining the entry and output matrices with the correct description as seen in Part 7 of Figure 3.2. When the code is analyzed by the C++ compiler, the call corresponding to Part 1 of Figure 3.2 will end up triggering the generation phase.The routine in Part 6 is generated using a combination of the parameters mentioned in Parts 2 to 5. This routine can be a LAPACK/MAGMA kernel (if available, which is the case of Figure3.2) or a kernel directly implemented in N T 2 . After the C++ compilation phase, we obtain a code similar to the one given in Part 8.

xFigure 2 . 4 :

 24 Figure 2.4: Performance comparison between LAPACK/MAGMA routines and generated codes via NT 2 for general dense linear system solution

  t a b l e <double> mcsne ( t a b l e <double> c o n s t& A, t a b l e <double> c o n s t& B) { d o u b l e anrm = l a n g e (A, ' I ' ) ; d o u b l e c t e = anrm * Eps<double >() * nt2 : : s q r t ( width ( a ) ) ; t a b l e <f l o a t > SA = c a s t <f l o a t >(A) ; t a b l e <f l o a t , u p p e r _ t r i a n g u l a r _ > SR = t r i u ( q r (SA) ) ; t a b l e <f l o a t > SX = mtimes ( t r a n s (SA) , c a s t <f l o a t >(B) ) ) ; SX = l i n s o l v e ( t r a n s (SR) ,SX) ; SX = l i n s o l v e (SR ,SX) ; t a b l e <double> X = c a s t <double >(SX) ; t a b l e <double> E = Bmtimes (A,X) ; s t d : : s i z e _ t i = 0 ; do { SX = c a s t <f l o a t >(mtimes ( t r a n s (A) , c a s t <f l o a t >(E) ) ) ; SX = l i n s o l v e ( t r a n s (SR) ,SX) ; SX = l i n s o l v e (SR ,SX) ; E = c a s t <double >(SX) d o u b l e RNRM = maximum( abs (E(_) ) ) ; X += E ; d o u b l e XNRM = maximum( abs (X(_) ) ) ; E = Bmtimes (A,X) ; ++i ; } w h i l e ( not (RNRM < XNRM * c t e ) && ( i <max_iter ) ) ; r e t u r n X; } Listing 2.6: N T 2 implementation for MCSNE

  i n t lda ,m, na , nhrs , l d b // s t a n d a r d l a p a c k d e f i n i t i o n s i n t nb = magma_get_dgeqrf_nb (m) ; i n t lwkopt = ( 2 * s t d : : min (m, na )+ ( ( na +31) / 3 2 ) * 3 2 ) * nb ; nt2 : : c o n t a i n e r : : t a b l e <double> x ( nt2 : : o f _ s i z e ( ldb , n h r s ) ) ; nt2 : : c o n t a i n e r : : t a b l e <f l o a t > tau ( nt2 : : o f _ s i z e ( s t d : : min (m, na ) , 1 ) ) ; nt2 : : memory : : cuda_buffer <double> dA(m * na , a ) ; nt2 : : memory : : cuda_buffer <double> dB( l d b * nhrs , b ) ; nt2 : : memory : : cuda_buffer <double> dX( l d b * n h r s ) ; nt2 : : memory : : cuda_buffer <double> dE( l d b * n h r s ) ; nt2 : : memory : : cuda_buffer <double> temp ( l d b * n h r s ) ; nt2 : : memory : : cuda_buffer <f l o a t > dSR (m * na ) ; nt2 : : memory : : cuda_buffer <f l o a t > dSX( l d b * n h r s ) ; nt2 : : memory : : cuda_buffer <f l o a t > dT( lwkopt * n h r s ) ; d o u b l e e p s = b o o s t : : simd : : Eps<double >() ; d o u b l e anrm = magmablas_dlange ( MagmaInfNorm , m, na , dA . data ( ) , m, dE . data ( ) ) ; d o u b l e c t e = anrm * e p s * nt2 : : s q r t ( type_t ( na ) ) ; d o u b l e xnrm , rnrm ; magmablas_dlag2s (m, na , dA . data ( ) ,m, dSR . data ( ) ,m, &i n f o ) ; magmablas_dlag2s ( ldb , nhrs , dB . data ( ) , ldb , dSX . data ( ) , ldb ,& i n f o ) ; magmablas_sgemv ( MagmaTrans , m, na , one , dSR . data ( ) , m, dSX . data ( ) , one , 0 , dSX . data ( ) , one ) ; magma_sgeqrf2_gpu ( m, na , dSR . data ( ) , m, tau . data ( ) , &i n f o ) ; magmablas_strsm ( MagmaLeft , MagmaUpper , MagmaTrans , MagmaNonUnit , na , nhrs , ( f l o a t ) one , dSR . data ( ) ,m, dSX . data ( ) , na ) ; magmablas_strsm ( MagmaLeft , MagmaUpper , MagmaNoTrans , MagmaNonUnit , na , nhrs , ( f l o a t ) one , dSR . data ( ) ,m, dSX . data ( ) , na ) ; magmablas_slag2d ( ldb , nhrs , dSX . data ( ) , ldb , dX . data ( ) , ldb ,& i n f o ) ; magmablas_dlacpy ( MagmaFull , ldb , nhrs , dB . data ( ) , ldb , dE . data ( ) , l d b ) ; magmablas_dgemv ( MagmaNoTrans , ld a , na , negone , dA . data ( ) , lda , dX . data ( ) , one , one , dE . data ( ) , one ) ; f o r ( s i z e _ t i = 1 ; i <=10;++ i ) { magmablas_dgemv ( MagmaTrans , ld a , na , one , dA . data ( ) , l da , dE . data ( ) , one , 0 , temp . data ( ) , one ) ; magmablas_dlag2s ( ldb , nhrs , temp . data ( ) , ldb , dSX . data ( ) , ldb ,& i n f o ) ; magmablas_strsm ( MagmaLeft , MagmaUpper , MagmaTrans , MagmaNonUnit , na , nhrs , ( f l o a t ) one , dSR . data ( ) ,m, dSX . data ( ) , na ) ; magmablas_strsm ( MagmaLeft , MagmaUpper , MagmaNoTrans , MagmaNonUnit , na , nhrs , ( f l o a t ) one , dSR . data ( ) ,m, dSX . data ( ) , na ) ; magmablas_slag2d ( na , nhrs , dSX . data ( ) , ldb , dE . data ( ) , ldb ,& i n f o ) ; nt2 : : memory : : copy ( dE , x ) ; rnrm = nt2 : : maximum( nt2 : : abs ( x (_) ) ) ; magmablas_dgeadd ( na , nhrs , one , dE . data ( ) , ldb , dX . data ( ) , l d b ) ; nt2 : : memory : : copy (dX , x ) ; xnrm = nt2 : : maximum( nt2 : : abs ( x (_) ) ) ; magmablas_dlacpy ( MagmaFull , ldb , nhrs , dB . data ( ) , ldb , dE . data ( ) , l d b ) ; magmablas_dgemv ( MagmaNoTrans , ld a , na , negone , dA . data ( ) , lda , dX . data ( ) , one , one , dE . data ( ) , one ) ; i f ( rnrm < xnrm * c t e ) { b r e a k ; } } nt2 : : memory : : copy (dX , x ) ; r e t u r n x (_( 1 , na ) ) ;

Figure 2 . 5 :

 25 Figure 2.5: Performance results of generated code on CPU

  mat A = randu<mat>(h , w) ; v e c B = randu<vec >(h ) ; d o u b l e n = 0 , n1 = 0 ; fmat Q,R; Q. s e t _ s i z e ( h , w) ; R. s e t _ s i z e ( h , w) ; Q. z e r o s ( ) ; R. z e r o s ( ) ; fmat a = conv_to<fmat > : : from (A) ; qr_econ (Q, R, a ) ; f v e c b = conv_to<f v e c > : : from (B) ; f v e c sx = t r a n s ( a ) * b ; sx = s o l v e ( t r a n s ( t r i m a t u (R) ) , sx ) ; sx = s o l v e ( t r i m a t u (R) , sx ) ; v e c X = conv_to<vec > : : from ( sx ) ; v e c e = B -A * X ; f o r ( s t d : : s i z e _ t i =1 ; i <5 ; ++i ) { sx = conv_to<f v e c > : : from ( t r a n s (A) * e ) ; sx = s o l v e ( t r a n s ( t r i m a t u (R) ) , sx ) ; sx = s o l v e ( t r i m a t u (R) , sx ) ; e = conv_to<vec > : : from ( sx ) ; X = X + e ; e = B -A * X; } Listing 2.8: CPU Implementation of MCSNE with Armadillo

  MatrixXd A = MatrixXd : : Random ( h , w) ; VectorXd B = VectorXd : : Random ( h ) ; d o u b l e n = 0 , n1 = 0 ; MatrixXf a = A. c a s t <f l o a t >() ; HouseholderQR<MatrixXf> qr ( a ) ; MatrixXf r f = q r . matrixQR ( ) ; VectorXf b = B . c a s t <f l o a t >() ; VectorXf sx = a . t r a n s p o s e ( ) * b ; s i z e _ t c o l s = a . c o l s ( ) ; s i z e _ t rows = a . rows ( ) ; MatrixXf s r ( c o l s , c o l s ) ; s i z e _ t s i z e _ n = s i z e o f ( f l o a t ) * ( c o l s ) ; f o r ( s i z e _ t i =0; i <c o l s ; i ++) { s t d : : memcpy ( s r . data ( )+i * c o l s , r f . data ( )+i * rows , s i z e _ n ) ; } sx = s r . t r i a n g u l a r V i e w <Eigen : : Upper >() . t r a n s p o s e ( ) . s o l v e ( sx ) ; sx = s r . t r i a n g u l a r V i e w <Eigen : : Upper >() . s o l v e ( sx ) ; VectorXd X = sx . c a s t <double >() ; VectorXd e = B -A * X ; f o r ( s t d : : s i z e _ t i =1 ; i <5 ; ++i ) { sx = (A. t r a n s p o s e ( ) * e ) . c a s t <f l o a t >() ; sx = s r . t r i a n g u l a r V i e w <Eigen : : Upper >() . t r a n s p o s e ( ) . s o l v e ( sx ) ; sx = s r . t r i a n g u l a r V i e w <Eigen : : Upper >() . s o l v e ( sx ) ; e = sx . c a s t <double >() ; X = X + e ; e = B -A * X; }

50 Chapter 2 .

 502 High-level programming for dense linear algebra

Figure 2 . 6 :

 26 Figure 2.6: Performance results of generated code on GPU

Figure 2 . 7 :

 27 Figure 2.7: Performance results on the Tegra X1 GPU

/

  / D e f i n e h o s t t a b l e t a b l e <f l o a t > A, B, C,D; // Triad k e r n e l A = B + C * D; Listing 3.1: NT 2 Triad kernel The code in Listing 3.1 corresponds to Part 1 of Figure 3.2.

Listing 3 . 2 : NT 2 Listing 3 . 3 :

 32233 t a b l e <f l o a t > A, B, C,D; // Triad k e r n e l r e p l a c e with t r a n s f o r m c a l l t r a n s f o r m (A, B + C * D) ; Triad transformFrom there, once we compile our code and go through phase 2 of Figure3.2 we obtain the following mangled code for our transform as described in Listing 3.3. NT 2 CUDA Triad mangled

Figure 3 . 3 :

 33 Figure 3.3: Triad kernel transform AST

  u s i n g b o o s t : : p r o t o : : c h i l d _ c ; u s i n g b o o s t : : p r o t o : : v a l u e ; s t d : : s i z e _ t s i z e = numel ( b o o s t : : p r o t o : : c h i l d _ c <0>(a1 ) ) ; s t d : : s i z e _ t b l o c k S i z e = s t d : : min ( s t d : : s i z e _ t ( 4 0 0 0 0 ) , s i z e ) ; s t d : : s i z e _ t nStreams = s t d : : min ( s t d : : s i z e _ t ( 2 ) , s i z e / b l o c k S i z e ) ; s t d : : s i z e _ t n = s i z e / b l o c k S i z e ; s t d : : s i z e _ t l e f t o v e r = s i z e % b l o c k S i z e ; dim3 blockDim = s t d : : min ( s t d : : s i z e _ t ( 1 0 2 4 ) , s i z e ) ; dim3 dimGrid = b l o c k S i z e / s i z e _ t ( 1 0 2 4 ) ; cudaStream_t stream [ nStreams ] ;

  ) ; b o o s t : : p r o t o : : v a l u e ( a0 ) . s p e c i f i c s ( ) . t r a n s f e r _ d t h ( a0 , i , stream [ j ] , j ) ; } i f ( l e f t o v e r =0) { . . . } Listing 3.7: NT 2 CUDA Triad Host Code streaming

  e h o s t t a b l e t a b l e <f l o a t , device_> A; // D e f i n e t h e t y p e o f r e d u c t i o n f u n c t o r <sum_> sum ; // R e d u c t i o n u s i n g t h e sum o f e l e m e n t s auto r e s u l t = g l o b a l ( sum ,A) ; Listing 3.8: NT 2 Reduction kernel

U

  v o i d nt2 : : e x t e r n a l _ k e r n e l <nt2 : : t a g : : global_ , nt2 : : t a g : : cuda_<b o o s t : : d i s p a t c h : : t a g : : cpu_> > :: c a l l <nt2 : : c o n t a i n e r : : t a b l e <s h o r t , nt2 : : device_ >, b o o s t : : simd : : t a g : : plus_ , s h o r t >( nt2 : : c o n t a i n e r : : t a b l e <s h o r t , nt2 : : device_ >&, b o o s t : : simd : : t a g : : plus_ &, s h o r t &) Listing 3.9: NT 2 CUDA Reduction demangled

#i n c

  l u d e <t h r u s t / r e d u c e . h> #i n c l u d e <t h r u s t / f u n c t i o n a l . h> #i n c l u d e <t h r u s t / d e v i c e _ p t r . h> v o i d t h r u s t _ r e d u c e 1 d e 0 a 5 e b ( s h o r t * b e g i n , s h o r t * end , s h o r t & r e s u l t ) { t h r u s t : : d e v i c e _ p t r <s h o r t > b_ptr = t h r u s t : : d e v i c e _ p o i n t e r _ c a s t ( b e g i n ) ; t h r u s t : : d e v i c e _ p t r <s h o r t > e_ptr = t h r u s t : : d e v i c e _ p o i n t e r _ c a s t ( end ) ; r e s u l t = t h r u s t : : r e d u c e ( b_ptr , e_ptr , r e s u l t , t h r u s t : : p l u s <s h o r t >() ) ; } Listing 3.10: NT 2 Reduction kernel • GPU : Tesla K40m -Pageable host to device (HTD) : 3 GB/s -Pinned host to device : 9.8 GB/s • Memory : 65 GB with a memcpy bandwidth of 5GB/s

  t a b l e <T> b l a c k s c h o l e s ( t a b l e <T> c o n s t& S , t a b l e <T> c o n s t& X , t a b l e <T> c o n s t& Ta , T c o n s t r , T c o n s t v ) { auto s = e x t e n t ( Ta ) ; t a b l e <T , device_ > d ( s ) , d1 ( s ) , d2 ( s ) ; t a b l e <T> r ;

2

 2

  .

Figure 3 . 4 :

 34 Figure 3.4: Black and Scholes Performance Comparison (in ms)

Figure 3 . 5 :

 35 Figure 3.5: Performance comparison on Tegra X1 with NT 2 (time in ms)

  t a b l e <T> b l a c k s c h o l e s ( t a b l e <T, pinned_> c o n s t& S , t a b l e <T, pinned_> c o n s t& X , t a b l e <T, pinned_> c o n s t& Ta , T c o n s t r , T c o n s t v ) { auto s = e x t e n t ( Ta ) ; t a b l e <T , device_ > d ( s ) , d1 ( s ) , d2 ( s ) ; t a b l e <T, pinned_> r ;

  t i e ( l , u ) = l u (A) x = l i n s o l v e ( l , b ) ; // l o w e r t r i a n g u l a r s o l v e x = l i n s o l v e ( u , x ) ; // upper t r i a n g u l a r s o l v e // One-s t e p r e f i n e m e n t e = b -nt2 : : mtimes (A, x ) ; e = nt2 : : l i n s o l v e ( l , e ) ; e = nt2 : : l i n s o l v e ( u , e ) ; x = x + e ; Listing 3.16: NT 2 LU linear solve with iterative refinement

Figure 3 . 6 :

 36 Figure 3.6: Performance comparison of NT 2 linear solve (in Gflop/s)

  t e m p l a t e <typename T> v o i d f ( v e c t o r <T>& v ) { f o r ( auto& x : v ) ++x ; } Listing 4.3: C++ generic for auto

  c o n s t e x p r l o n g l o n g f i b o n a c c i ( c o n s t i n t x ) { r e t u r n x <= 1 ? 1 : f i b o n a c c i ( x -1 ) + f i b o n a c c i ( x -2 ) ; } Listing 4.4: C++ constexpr

/

  / The [ . . . ] d e f i n e s a v a r i a d i c t e m p l a t e t e m p l a t e <typename . . . Values> c l a s s t u p l e ; Listing 4.5: C++ Variadic template Listing 4.5 is an example of a variadic template class that takes a variable number of arguments. Listing 4.6 allows for a more in-depths description of variadic templates. t e m p l a t e <typename T> T add (T v ) { r e t u r n v ; } t e m p l a t e <typename T, typename . . . Args> T add (T f i r s t , Args . . . a r g s ) { r e t u r n f i r s t + add ( a r g s . . . ) ; } Listing 4.6: C++ Variadic template add

Listing 4 . 7 :

 47 C++ Variadic template add call Examples of how the function add can be called are shown in Listing 4.7. // Apply a f o n c t i o n on a l i s t o f arguments t e m p l a t e <c l a s s F , c l a s s . . . Ts> F f o r _ e a c h _ a r g s (F f , Ts & & . . . a ) { r e t u r n ( v o i d ) s t d : : i n i t i a l i z e r _ l i s t <i n t >{(( v o i d ) s t d : : r e f ( f ) ( s t d : : forward <Ts >(a ) ) , 0 ) . . . } , f ; } Listing 4.8: C++ Variadic template function applyVaradic templates are a very powerful tool that has some very useful tricks. For example, in Listing 4.8, we have defined a function that take another function and a list of arguments and applies the function on each argument.

  t e m p l a t e <c l a s s T, T v> s t r u c t i n t e g r a l _ c o n s t a n t { s t a t i c c o n s t e x p r T v a l u e = v ; t y p e d e f T value_type ; t y p e d e f i n t e g r a l _ c o n s t a n t t y p e ; c o n s t e x p r o p e r a t o r value_type ( ) c o n s t n o e x c e p t { r e t u r n v a l u e ; } c o n s t e x p r value_type o p e r a t o r ( ) ( ) c o n s t n o e x c e p t { r e t u r n v a l u e ; } } ;

78 Chapter 4 .

 784 4 8 ; } t e m p l a t e <s t d : : s i z e _ t N, typename I> c o n s t e x p r l o n g l o n g p a r s e ( c o n s t c h a r (& a r r ) [N] , I c o n s t &, i n t ) { l o n g l o n g number = 0 , b a s e = 1 ; f o r ( s t d : : s i z e _ t i = 0 ; i < N; ++i ) { number += t o _ i n t ( a r r [N -1 -i ] ) * b a s e ; b a s e * = 1 0 ; } r e t u r n number ; } t e m p l a t e <c h a r . . . c> c o n s t e x p r auto o p e r a t o r " " _c ( ) { r e t u r n int_<p a r s e <s i z e o f . . . ( c ) >({ c . . . } , int_<s i z e o f . . . ( c ) >{} ,1) >{}; } Listing 4.10: C++ User defined literals Application to batched computations and tensor contractions

/

  / Template s p e c i a l i z a t i o n c o n s t e x p r auto l a y o u t = o f _ s i z e <5 ,3 >() ; // Using I n t e g r a l c o n s t a n t c o n s t e x p r auto l a y o u t 1 = o f _ s i z e ( 5 _c , 3 _c ) ; // Using dynamic d i m e n s i o n s c o n s t e x p r auto l a y o u t 2 = o f _ s i z e ( 5 , 3 ) ; // Using s t a t i c and dynamic d i m e n s i o n s c o n s t e x p r auto l a y o u t 3 = o f _ s i z e ( 5 _c , 3 ) ; // A c c e s s Dimensions a t c o m p i l e time c o n s t e x p r auto dim1 = l a y o u t ( 1 ) ; // W i l l not c a u s e a c o m p i l e time e r r o r s t a t i c _ a s s e r t ( x4 ( 1 ) == 5 , " f i r s t d i m e n s i o n i s s t a t i c " ) ; Listing 4.13: Dimensions for Tensors

/

  / C r e a t e a batch t h a t w i l l c o n t a i n 15 t e n s o r s o f s i z e 5 , 3 , 6 c o n s t e x p r auto batch<f l o a t , gpu_> b = make_batch ( o f _ s i z e ( 5 _c , 3 _c , 6 _c ) , 1 5 ) ; // A c c e s s i n g a t e n s o r from t h e batch r e t u r n s a view on i t c o n s t e x p r auto view_b = b ( 0 ) ; // C r e a t e a g r o u p i n g o f t e n s o r s o f same s i z e t e n s o r s c o n s t e x p r auto group<f l o a t , gpu_> g ( o f _ s i z e ( 5 _c , 3 _c ) ) ; // Add a t e n s o r t o t h e group c o n s t e x p r auto t e n s o r <f l o a t , gpu_> d_ts ( o f _ s i z e ( 5 _c , 3 _c ) ) ; g . push_back ( d_ts ) ; Listing 4.16: Batched tensors Once we have defined these functions we can call the kernel to compute a batched DGEMM on tensors of dimension 2. c o n s t e x p r auto batch<f l o a t , gpu_> b = make_batch ( o f _ s i z e ( 5 _c , 3 _c ) , 1 5 ) ; c o n s t e x p r auto batch<f l o a t , gpu_> b1 = make_batch ( o f _ s i z e ( 5 _c , 3 _c ) , 1 5 ) ; // Product o f two t e n s o r b a t c h e d o f d i m e n s i o n 2 f o r m a t r i x p r o d u c t u s i n g C++ o p e r a t o r c o n s t e x p r auto c = b * b1 ; // Product u s i n g t h e batch dgemm f u n c t i o n t h a t can be s p e c i a l i z e d d e p e n d i n g on p a r a m e t e r s c o n s t e x p r auto c = batch_gemm ( b , b1 ) ; Listing 4.17: Tensor Operations

Figure 4 . 2 :

 42 Figure 4.1: # of load instructions

Figure 4 . 3 :

 43 Figure 4.3: Total CPU instruction count

4. 4 . 85 Algorithm 9 16 1: 2 : 5 : 6 :

 485916256 Optimizing small matrix-matrix products Generic matrix-matrix product applied to matrices of size 16 × Load B0, B1, B2, B3 Load α, β 3: S = 16 4: for i = 0, 1, ... , S-1 do Load A[i*S] Mi0 = A[i*S] * B0; ... Mi3 = A[i*S] *B3 Mi0 += A[i*S+u] * Bu0; ... Mi3 += A[i*S+u] *Bui3 Mi0 = α Mi0 + β (Load Ci0); ... Mi3 = α Mi3 + β (Load Ci3) 13:

Figure 4 . 4 :

 44 Figure 4.4: Effect of the number of CPU cores

Figure 4 . 5 :

 45 Figure 4.5: Effect of the NUMA memory management

Figure 4

 4

Figure 4 . 6 :

 46 Figure 4.6: Experimental results of the matrix-matrix multiplication on CPU

Figure 4 . 7 :

 47 Figure 4.7: Experimental results of the matrix-matrix multiplication on the Tegra X1

Figure 4 . 8 :

 48 Figure 4.8: GPU performance on K40

  Mots clefs : programmation générique, programmation générative, C++, algèbre linéaire, GPU Résumé : Les architectures parallèles sont aujourd'hui présentes dans tous les systèmes informatiques, allant des smartphones aux supercalculateurs en passant par les ordinateurs de bureau. Programmer efficacement ces architectures requiert un effort pluridisciplinaire portant sur les langages dédiés (Domain Specific Languages -DSL), les techniques de génération de code et d'optimisation, et les algorithmes numériques propres aux applications. Dans cette thèse, nous présentons une méthode de programmation haut niveau prenant en compte les caractéristiques des architectures hétérogènes ainsi que les propriétés des matrices pour produire un solveur générique d'algèbre linéaire dense. Dans la mesure où les GPUs sont devenus un outil important pour le calcul haute performance, il est essentiel de les utiliser dans les plateformes de calcul. Nous avons par la suite étendu nos travaux à un modèle de programmation multi-étape ("multistage") pour résoudre les problèmes d'interopérabilité entre les modèles de programmation CPU et GPU. Nous utilisons cette technique pour générer automatiquement du code pour accélérateur à partir d'un code effectuant des opérations point par point ou utilisant des squelettes algorithmiques. Enfin, nous montrons comment la programmation haut niveau peut être appliquée à des calculs groupés et des contractions de tenseurs. Tout d'abord, nous expliquons comment concevoir un modèle de conteneur en utilisant des techniques de programmation basées sur le C++ moderne (C++-14). Ensuite, nous avons implémenté un produit de matrices optimisé pour des matrices de petites tailles en utilisant des instructions SIMD.

  Contents1.1 High level Programming for Scientific Computing . . . . . .

7 1.2 Solving dense linear systems . . . . . . . . . . . . . . . . . . . 8 1.2.1 LU factorization . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.2.2 QR factorization . . . . . . . . . . . . . . . . . . . . . . . . . 11 1.2.3 Solvers in LAPACK and MAGMA . . . . . . . . . . . . . . . 12 1.3 Architectures : Design and evolution . . . . . . . . . . . . . 15 1.3.1 SIMD extensions . . . . . . . . . . . . . . . . . . . . . . . . . 15 1.3.2 Multicore systems . . . . . . . . . . . . . . . . . . . . . . . . 16 1.3.3 Graphics Processing Units (GPU) . . . . . . . . . . . . . . . 17 1.3.4 Accelerated Processing Units . . . . . . . . . . . . . . . . . . 19 1.4 Programming tools and techniques . . . . . . . . . . . . . . . 20 1.4.1 Low-level programming tools . . . . . . . . . . . . . . . . . . 21 1.4.2 Domain Specific Languages . . . . . . . . . . . . . . . . . . . 24 1.4.3 Domain Specific Embedded Languages . . . . . . . . . . . . 25 1.4.4 Multi-stage programming . . . . . . . . . . . . . . . . . . . . 26 1.5 Considered language and library . . . . . . . . . . . . . . . . 26 1.5.1 C++ language and programming techniques . . . . . . . . . 26 1.5.2 The NT 2 library . . . . . . . . . . . . . . . . . . . . . . . . . 27 1.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

  Algorithm 3 illustrates the QR factorization as described by Golub using the house function[57, p.237].

Algorithm 3 HouseHolder QR in place, no pivoting Input: A is a m × n matrix 1: for k = 1 : n do 2:

Table 2 .

 2 1: Configuration space parameter levels

	0-Matrix type	general | band | diagonal | symmetric | positive definite
	1-Data type	float | double | single/double complex
	2-Precision	fixed | mixed-precision
	3-Conditioning	no information | ill-conditioned
	4-Storage scheme	general | packed
	5-Architecture	CPU | GPU

  Related Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

	Chapter 3
	Multi-stage programming in C++
	for GPGPUs
	Contents
	3.1

3.2 Multi-stage programming in C++ . . . . . . . . . . . . . . . 56 3.2.1 Designing a multi-stage software . . . . . . . . . . . . . . . . 56 3.2.2 Multi-stage programming in C++ . . . . . . . . . . . . . . . 57 3.2.3 Multi-stage programming tool for NT 2 . . . . . . . . . . . . 58 3.3 Software implementation for NT 2 . . . . . . . . . . . . . . . 59 3.4 Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 3.4.1 Black & Scholes kernel . . . . . . . . . . . . . . . . . . . . . . 65 3.4.2 Linsolve kernel . . . . . . . . . . . . . . . . . . . . . . . . . . 68 3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

  2 ) to O(n), compared to a classical ijk or ikj implementation as depicted in Figures 4.1, 4.2, and 4.3, respectively.

			gen code
	Number of instructions	10 6 10 7 10 8	mkl code ijk code ikj code
		0	5	10	15	20	25	30	35
					Matrix Size

NT2 is available at http://github.com/jfalcou/nt2

http://www.sandia.gov/∼tgkolda/TensorToolbox/

http://grtensor.phy.queensu.ca/

http://www.math.washington.edu/∼lee/Ricci/

http://en.cppreference.com/w/cpp/types/integral_constant

http://www.boost.org/doc/libs/1_61_0/libs/mpl/doc/index.html

https://ericniebler.github.io/std/wg21/D4128.html

http://www.top500.org/

http://www.netlib.org/utk/people/JackDongarra/PAPERS/sunway-report-2016.pdf

Université Paris-Saclay Espace Technologique / Immeuble Discovery Route de l'Orme aux Merisiers RD 128 / 91190 Saint-Aubin, France

LD/ST SFU LD/ST LD/ST LD/ST

with ViennaCL at runtime. [START_REF] Lutz | Boost.Compute[END_REF] is a header only C++ library based on the OpenCL standard. Similar to other libraries, it manages device memory through a designated container. It provides an interesting concept of future for asynchronous copy on the device allowing for more versatility. Boost.Compute also supports closures and adaptable structures for device. Similarly to thrust, it is a well rounded library based on OpenCL to simplify the coding process on accelerators. It however lacks support for numerical analysis and cannot generate CUDA code.

Boost.Compute

Eigen [START_REF] Guennebaud | A C++ Linear Algebra Library[END_REF] is a popular library to solve linear systems using expression templates. It should be able to support accelerator code with CUDA by writing Eigen code in a .cu file. It however does not provide any high-level interface or special containers for GPU computing. Since Eigen 3.3 it has a partial experimental CUDA support for Eigen's objects and algorithms within CUDA kernels. Eigen does not support OpenCL.

Kokkos [START_REF] Edwards | Kokkos: Enabling manycore performance portability through polymorphic memory access patterns[END_REF] is a C++ developed by Sandia that provides an API for hybrid computations. It does not follow the guidelines of standard C++ libraries and uses view based on a shared pointer. It uses parallel skeletons and lambda as their parallel model. The advantage of Kokkos is that their API is rather simple and can be directly compiled with nvcc simplifying integration. It however does not provide any AST and high level operations.

SciPAL [START_REF] Kramer | SciPAL: Expression templates and composition closure objects for high performance computational physics with cuda and openmp[END_REF] is an expression template library with a DSEL for dense and sparse linear algebra. It focuses mainly on recognizing GEMM calls by parsing the AST and regular operations on containers. Its code generation process for CUDA kernels is done mostly by explicitly writing the code in callable objects. The code generation process is then done at runtime for non-BLAS kernels. SciPAL does not support OpenCL. In Figure 3.5, we compare features between NT 2 and the previously described libraries. We did not include SciPAL since the code is not available, and Eigen as This method can be extended simply as we have shown. Any functions implementing external_kernel will end up triggering our MSP process. Each function being differentiated by the function tag used. It is then possible to extend the symbol/code converter to support this new function. Each function can be mixed with one another and also used the MAGMA back-end we described in Section 2.4. Listing 3.12 is an example of a code doing a reduction and an elementwise operations that will go through our MSP process. The demangled object will contain both missing symbol for global and transform. The symbol/code converter will parse each symbol and generate the corresponding code.

Experiments

In this section, we will show that our generation process produces satisfactory performances in most situations. The benchmarks are realized with the following components :

• CPU : 2 x 6 cores Intel Xeon E5-2620 15MB L3, AVX Chapter 4 We can then freely extract the dimensions (Listing 4.13) and use them to specify our CPU and GPU kernels at compile time. Our tensor model is based on our layout, the data type of the tensor and its locality. The memory buffer is based on vector from the Standard Template Library for CPU and Thrust for GPU. To generate a tensor, we need to pass a data type and locality as template parameter and the size to the constructor (Listing 4.15). // C r e a t e a rank 2 t e n s o r o f s i z e 5 , 3 on GPU c o n s t e x p r t e n s o r <f l o a t , gpu_> d_ts ( o f _ s i z e <5 ,3 >() ) ; // C r e a t e a rank 2 t e n s o r o f s i z e 5 , 3 on CPU c o n s t e x p r t e n s o r <f l o a t > t s ( o f _ s i z e <5 ,3 >() ) ; // Use t h r u s t t o f i l l d_ts with 9 t h r u s t : : f i l l ( d_ts . b e g i n ( ) , d_ts . end ( ) , 9 ) ; // Copy d_ts from GPU t o t s on CPU copy ( d_ts , t s ) ; Listing 4.15: Create Tensor and copy Transfers between CPU and GPU tensors can be expressed through the function copy (Listing 4.15). This function will use the stream 0 by default but a stream can be passed for asynchronous copy. Locality definitions are simple classes that contain the semantic information necessary. For example, gpu_ will be a class definition that will inherit from the locality class itself. It is then possible to iterate on each parameter of the layout looking for a locality class. Meta-functions are available in the C++ standard library (utility package) to do type checking. For example, std :: is_base_of verifies if a class inherits from another one. If we cannot find a corresponding locality in the class definition, we can then fall-back to the default locality we have decided.

Application to batched computations and tensor contractions

We have designed two models for batched computing (Listing 4.16). The first one is based on allocating a single memory block for all tensors to improve data locality while the other is a group of same size tensors.

Chapter 4. Application to batched computations and tensor contractions

deed, the compiler will only be able to unroll the loops properly at compile time if the number of iterations in the loops are fixed.The inner loop can only be factorized in such a way because we can trust the compiler. Most C++ compilers use the polyhedral model to optimize loops [START_REF] Pouchet | Iterative optimization in the polyhedral model: Part ii, multidimensional time[END_REF] and compile time information is needed for proper optimizations.

Effect of the Multi-threading

As described above, operating on matrices of very small sizes is memory-bound computation and thus, increasing the number of CPU cores may not always increase the performance since the performance will be limited by the bandwidth which can