Monsieur Raphaël Couturier

Julien Hascoet

Nicolas Brunie

Julien Lemaire

Romarik Jodin

Jérôme Reybert

Pierre Guironnet De Massas

Clément Leger

Michael Mercier

Van Toan Dao

Baptiste Jonglez

Elodie Morin

Henry-Joseph Audeoud

Pierre Brunisholz

G Van Zee

Robert A Van De Geijn

Devangi Parikh

Tyler Michael Smith

Keywords: ), BLAS Basic Linear Algebra Subprograms, BLIS BLAS-like Library Instantiation Software, CPU Central Processing Unit, DLA Dense Linear Algebra, GEMM General matrix multiplication, GPU Graphical Processing Unit, HPC High Performance Computing, HPL High Performance Linpack, KNL Knights Landing, LAPACK Linear Algebra PACKage, LBM Lattice Boltzmann method, MPI Message Passing Interface, MPPA Massively Parallel Processor Array, NUMA Non-uniform Memory Access, OpenCL Open Computing Language, OpenMP Open Multi-Processing, xvii Pthreads POSIX threads, ScaLAPACK Scalable LAPACK

Upcoming Exascale target in High Performance Computing (HPC) and disruptive achievements in artificial intelligence give emergence of alternative non-conventional many-core architectures, with energy efficiency typical of embedded systems, and providing the same software ecosystem as classic HPC platforms. A key enabler of energy-efficient computing on many-core architectures is the exploitation of data locality, specifically the use of scratchpad memories in combination with DMA engines in order to overlap computation and communication. Such software paradigm raises considerable programming challenges to both the vendor and the application developer. In this thesis, we tackle the memory transfer and performance issues, as well as the programming challenges of memory-and compute-intensive HPC applications on the Kalray MPPA many-core architecture.

With the first memory-bound use-case of the lattice Boltzmann method (LBM), we provide generic and fundamental techniques for decomposing three-dimensional iterative stencil problems onto clustered many-core processors fitted with scratchpad memories and DMA engines. The developed DMA-based streaming and overlapping algorithm delivers 33% performance gain over the default cache-based implementation. Highdimensional stencil computation suffers serious I/O bottleneck and limited on-chip memory space. We developed a new in-place LBM propagation algorithm, which reduces by half the memory footprint and yields 1.5 times higher performance-per-byte efficiency than the state-of-the-art out-of-place algorithm.

On the compute-intensive side with dense linear algebra computations, we build a matrix multiplication benchmark based on exploitation of scratchpad memory and efficient asynchronous DMA communication. This program delivers 350 GFLOPS, or 86% of theoretical performance of the MPPA. These techniques are then extended to a DMA module of the BLIS framework, which allows us to instantiate an optimized and portable level-3 BLAS numerical library on any DMA-based architecture, in less than 100 lines of code. We achieve 75% peak performance on the MPPA processor with the matrix multiplication operation (GEMM) from the standard BLAS library, without having to write thousands of lines of laboriously optimized code for the same result.

First, I would like to express my respectful thanks to my advisors: Bernard Tourancheau for his experience, patience and methodological research guidance; Christian Obrecht for his expertise in numerical simulation, relevant questions and thorough correction

Résumé

La prochaine cible de Exascale en calcul haute performance (High Performance Computing -HPC) et des récent accomplissements dans l'intelligence artificielle donnent l'émergence des architectures alternatives non conventionnelles, dont l'efficacité énergétique est typique des systèmes embarqués, tout en fournissant un écosystème de logiciel équivalent aux plateformes HPC classiques. Un facteur clé de performance de ces architectures à plusieurs coeurs est l'exploitation de la localité de données, en particulier l'utilisation de mémoire locale (scratchpad) en combinaison avec des circuits d'accès direct à la mémoire (Direct Memory Access -DMA) afin de chevaucher le calcul et la communication. Un tel paradigme soulève des défis de programmation considérables à la fois au fabricant et au développeur d'application. Dans cette thèse, nous abordons les problèmes de transfert et d'accès aux mémoires hiérarchiques, de performance de calcul, ainsi que les défis de programmation des applications HPC, sur l'architecture pluri-coeurs MPPA de Kalray.

Pour le premier cas d'application lié à la méthode de Boltzmann sur réseau (Lattice Boltzmann method -LBM), nous fournissons des techniques génériques et réponses fondamentales à la question de décomposition d'un domaine stencil itérative tridimensionnelle sur les processeurs clusterisés équipés de mémoires locales et de circuits DMA.

Nous proposons un algorithme de streaming et de recouvrement basé sur DMA, délivrant 33% de gain de performance par rapport à l'implémentation basée sur la mémoire cache par défaut. Le calcul de stencil multi-dimensionnel souffre d'un goulot d'étranglement important sur les entrées/sorties de données et d'espace mémoire sur puce limitée. Nous avons développé un nouvel algorithme de propagation LBM sur-place (in-place). Il consiste à travailler sur une seule instance de données, au lieu de deux, réduisant de moitié l'empreinte mémoire et cède une efficacité de performance-par-octet 1.5 fois meilleure par rapport à l'algorithme traditionnel dans l'état de l'art.

Du côté du calcul intensif avec l'algèbre linéaire dense, nous construisons un benchmark de multiplication matricielle, basé sur l'exploitation de la mémoire locale et la communication DMA asynchrone. Ce programme atteint 350 GFLOPS, soit 86% de la performance théorique de MPPA. Ces techniques sont ensuite étendues à un module DMA générique du framework BLIS, ce qui nous permet d'instancier une bibliothèque BLAS3 (Basic Linear Algebra Subprograms) portable et optimisée sur n'importe quelle architecture basée sur DMA, en moins de 100 lignes de code. Nous atteignons une performance maximale de 75% du théorique sur le processeur MPPA avec l'opération de multiplication de matrices (GEMM) de BLAS, sans avoir à écrire des milliers de lignes de code laborieusement optimisé pour le même résultat. The first four chapters in this manuscript introduce the current High Performance Computing (HPC) situation, the state-of-the-art and objectives of this thesis, and the target many-core platform. Contribution chapters in this manuscript are then presented as two main parts, organized in the thematic order and not in the chronological one. The first part presents approaches in optimizing data transfer and memory footprint of the threedimensional Lattice Boltzmann method (LBM), which belongs to the memory-bound category of applications. The second part focuses on Dense Linear Algebra (DLA) operations and associated numerical libraries, belonging to the compute-bound class.

Experimental results in this document are reported on the Kalray Massively Parallel

Processor Array (MPPA) many-core architecture, as well as other latest mainstream computing platforms such as NVIDIA Pascal Graphical Processing Unit (GPU), Intel

Xeon Haswell Non-uniform Memory Access (NUMA) Central Processing Unit (CPU)

and Intel Xeon Phi Knights Landing (KNL) processor.

Some contribution chapters in this manuscript contain principal materials from published and submitted papers, authored or co-authored by the writer of this manuscript as well. These papers are: Obrecht. BLIS-RDMA: A portable and high performance level-3 BLAS for DMAbased many-core architectures. In submission, 2017.

Chapter 1

High Performance Computing:

from Single-core to Many-core Knowledge has a beginning, but no end.

-Geeta S. Iyengar.

Introduction

Micro-processor architectures have made a considerable evolution since the first days of computer science. Several decades ago, when the processor clock was of the same order as the main memory speed, about hundreds of kilohertz to several megahertz, performance concerns were mostly on the computation cycles, rather than optimizing memory transfers. With advances in semiconductors, the transistor number and computing performance has exponentially increased over the years. Moore's law (Fig. 1.1) is a model of the semiconductor evolution, which depicts the increasing transistor count that doubles approximatively every 24 months for a constant circuit price. This evolution results in a parallel increase in computing power. Despite the limit in transistor size, this law has been nevertheless sustained by processor manufacturers since the last decade, by increasing the clock frequency and integrating more cores into a same silicon die. High Performance Computing (HPC) has become an essential field to guide and to be guided by the moving computing architectures. Nowadays, HPC is almost used in any scientific domain, from molecular dynamics simulation, bio-informatics, medical drug discovery, to computational fluid dynamics (CFD), oil and gas, image and signal processing, ocean simulation and weather forecast, and recently huge explosion in computing workload for astrophysics radio telescopes (Square Kilometer Array -SKA) or deep learning. Nonetheless, performance walls in exploiting parallelism are, were, and will be the main challenge to computer scientists on the road of developing future HPC systems.

Performance walls 1.2.1 Frequency-wall

The higher the clock frequency operates, the higher the performance will be. However, there is a physical limit between the clock frequency and the processor design. In synchronous circuits, there exists a maximal frequency, under which the processor still manages to synchronize its components properly with the clock signal and to maintain an operating state. Beyond that frequency, the various distance between the clock source to the working components introduces a micro-dephasing that is in the same order of a single clock rate, results in erroneous synchronization and unstable circuits. That maximal clock frequency is conventionally fixed at 10 GHz. In this case, a single clock rate is roughly the time for the light to travel three centimeters, therefore the maximum possible distance difference on an electronic circuit.

Energy consumption and power dissipation is also an issue. High-frequency working system produces heat and must be accordingly cooled down. Current leakage is proportional to the voltage and increases power consumption. World records in over-clocking end-user computers are often limited at 8 GHz. To reach this level, the system is often cooled down by liquid nitrogen. An alternative to boost performance was, instead, doubling the core count and introducing execution concurrency, which opened the era of multi-core processors.

Energy-wall

In the last years, the energy-wall has emerged as the main limiting factor in designing supercomputers. Today, we are in the age of Petascale (10 15 floating-point operations per second) within an energy budget of 20 MW. The first-ranked machine in the Top500 list (list of most powerful supercomputers in the world) delivers currently a power/energy ratio of about 10 GFLOPS/W. To reach Exascale (10 18 floating-point operations per second), predicted to appear in the 2020s (Fig. 1.2), we need a system capable to deliver more than 50 GFLOPS/W, with a global power load of less than 20 MW.

To lower power consumption and to increase flops count, the only solution is to embed a large number of low-frequency cores into a same processor, sometimes up to hundreds of cores on a same die. Such a processor would require disruptive memory and interconnection technologies to feed the core cluster. This arises the question of how to write efficient codes on those architectures, which programming model to design to expose the Source: https://www.top500.org.

massive parallelism, while ensuring ease of coding. It brings into light the challenges of many-core design and programming.

Memory-wall

There is an increasing gap between the computing performance and the memory speed.

The computing performance slope is much steeper than the memory one. Since the overall performance amounts to the weakest part of the whole system, memory bandwidth turns into a bottleneck for most applications. In order to better understand the former issue, we will use the concept of arithmetic intensity and the roofline model [6]. The roofline model provides an easy way to predict performance of an application or more precisely its computation kernel. The concept of arithmetic intensity (I) (ops/byte or flops/byte for floating-point) is defined by the number of arithmetic operations performed on a given quantity of data loaded. Each implementation of a numerical method possesses a proper I related to its computation workload and data access pattern. With a given compiler and predefined optimization flags, the application produces a binary code that reaches a performance π . Running a memory bandwidth benchmark on the system gives a sustained memory bandwidth β (π < π and β < β). The real performance P of the computation kernel is equal to min(π , β × I). When P = π , the kernel is considered compute-bound since its performance does not rely on the memory but on the computing capacity of the processor, the applied compilation options, and the parallelization or vectorization method. This performance is a constant for each configuration variant (horizontal line in the roofline figure, see Fig. 1.5). When P = β × I, P follows a sloping line with a slope β , where comes the name of roofline. In this case, the kernel is considered memory-bound, since its performance relies on the memory bandwidth.

Each processor architecture has its own roof lines, and every computation kernel running on that processor has a performance bounded by those two lines. This performance model, despite being simple to characterize, allows doing a rapid comparison of computation kernels on an architecture (which ones are memory-bound, which ones are compute-bound), or comparing a specific kernel across multiple architectures (if it is memory-bound, so it should be better to use a higher memory-bandwidth system). Using the roofline model, a developer or an integration engineer can make appropriate choice to use/port their code on a specific platform, to adopt a more aggressive optimization method, or even to make architectural decisions. To overcome the memory latency and keep cores busy, computer architectures implement data prefetching. As can be seen in Fig. 1.6, near-core memory levels (caches or scratchpad, of the order of kilobyte and megabyte) can prefetch data either by implicit hardware mechanism or explicit built-in instructions. They were demonstrated to work efficiently on contemporary CPUs. However, when off-chip memories (DDR, NVMe etc.) come up with supplementary physical media (NoC, PCIe), performing implicit hardware prefetching is much harder and not relevant anymore, as the memory scope is becoming too large (order of giga-byte and tera-byte). Instead, a software approach based on asynchronous RDMA communication libraries, despite requiring additional programming efforts, can improve data locality and deliver satisfactory performance.

Software-wall

Exposing parallelism while keeping the programming model simple is a hard question.

NUMA CPU processors can be programmed with Open Multi-Processing (OpenMP), POSIX threads (Pthreads), Intel Threading Building Blocks (TBB) or Intel Cilk. With the rise of heterogeneous systems, programming models must evolve as well. The OpenMP 3 for shared and NUMA architecture was revised to OpenMP 4 to support target devices. OpenACC [7], CUDA and Open Computing Language (OpenCL) [8] are also other programming models and APIs for accelerators. Programming language for distributed memory has longly been dominated by Message Passing Interface (MPI) [9] with incremental features leveraging the increasing node count: from two-sided communication in MPI 1 to one-sided communication since MPI 2, non-blocking collective operation since MPI 3 as well as a re-enforced one-sided specification.

Code and performance portability is also an issue due to the hardware diversity. To reduce the programming efforts for non-computer scientists, additional tools and metalanguages were introduced: domain specific languages (DSL) like Halide [START_REF] Ragan-Kelley | Halide: a language and compiler for optimizing parallelism, locality, and recomputation in image processing pipelines[END_REF], high-level language and automatic framework (SYCL [START_REF] Keryell | Khronos SYCL for OpenCL: a tutorial[END_REF]), data-flow and DAG-based (Directed Acyclic Graph) analysis tools (StarPU [START_REF] Augonnet | StarPU: a unified platform for task scheduling on heterogeneous multicore architectures[END_REF]). These high-level tools allow users to express their processing kernels and get their code automatically generated, compiled and deployed on multiple computation units.

We believe that future high performance systems will be a combination of various computing platforms. Applications (and their sub-modules), upon their arithmetic intensity range, will be deployed and run on the most suitable platform. Such a system will be highly heterogeneous and non-uniform in terms of performance and memory bandwidth of each sub-platform, where developers sometimes need to hand-tune the code to obtain the best performance, especially on embedded and non-conventional hardware.

Programming and running these all-in-one systems raises considerable complexity in design, scheduling, debugging, isolation and security.

Fault tolerance will also be a big concern on large-scale systems to resist against fail-stop failures, due to the fast decline of mean time between failure (MTBF) with the growing system size. Hardware-based fault tolerance mechanisms tend to be vendor-dependent at a certain degree. Without advances in research, these mechanisms may suffer significant overhead and would be difficult to optimize. Algorithm-Based Fault Tolerance (ABFT), since the last few years, has achieved important results by the research community [START_REF] Chen | Online-ABFT: An online algorithm based fault tolerance scheme for soft error detection in iterative methods[END_REF][START_REF] Du | Algorithmbased fault tolerance for dense matrix factorizations[END_REF][START_REF] Bosilca | Algorithm-based fault tolerance applied to high performance computing[END_REF][START_REF] Wunderlich | Efficacy and efficiency of algorithm-based fault-tolerance on GPUs[END_REF][START_REF] Rech | An efficient and experimentally tuned software-based hardening strategy for matrix multiplication on GPUs[END_REF]. However, relevance of following these approaches to the coming Exascale context is still unclear as of today.

Summary

In this chapter, we briefly present the high performance computing (HPC) and its main three performance walls: the frequency-wall, the energy-wall and the memory-wall. We also identify software challenges for the upcoming HPC applications, that we believe to be another obstacle to performance: the software-wall.

In the next chapter, we will introduce the first type of HPC application studied in this thesis: the lattice Boltzmann method, a stencil-based computation, known to be one of the most memory-bound applications.

Chapter 2

Lattice Boltzmann method (LBM)

If you can't explain it simply, you don't understand it well enough.

-Albert Einstein.

Background

Theory

Inspired from the lattice gas automata theory [START_REF] Frisch | Lattice-gas automata for the Navier-Stokes equation[END_REF] and first introduced by McNamara and Zanetti [START_REF] Guy | Use of the Boltzmann equation to simulate lattice-gas automata[END_REF], the lattice Boltzmann method has become widely used in computational fluid dynamics (CFD) as an alternative to the solving of Navier-Stokes equations.

Belonging to the structured grid-based discrete method, the LBM is known for its advantages such as straightforward meshing, ability to model complex geometries, and most of all its inherent parallelism, well-suited to massively parallel computing architectures.

An LBM model is characterized by a stencil type, denoted DdQq, where d is the number of space dimensions (one, two or three) and q is the number of particle distribution functions (PDFs) [START_REF] Cao | Physical symmetry and lattice symmetry in the lattice Boltzmann method[END_REF]. Particle distribution functions describe the interaction between a lattice node and its surrounding neighborhood. More precisely, q relates to the number of neighboring nodes that will be involved into interaction with the lattice node of interest.

In most cases, the node itself is taken into account and the number of neighboring nodes equals q -1. The most used stencil types for LBM are D2Q5 and D2Q9 for two-dimensional domains, or D3Q19 and D3Q27 for three-dimensional domains (see Fig. 2.1). Three-dimensional LBM often operates on D3Q19 or D3Q27 stencils. The LBM spatial domain is represented by a grid of nodes, discretized with a mesh size δx. The simulation duration is discretized in constant time steps δt. The LBM updating rule for each node at each time step is defined by the following equation:

f i (x + δtξ i , t + δt) -f i (x, t) = Ω f i (x, t) (2.1)
in which Ω is a pre-defined collision operator. The collision operator implements the time evolution of particle distribution functions f i (i ∈ {0, ..., q -1}) of a given node towards its nearest neighbors with respect to the ξ i velocities. For better presentation in LBM codes, Eq. 2.1 is often split into two sub-steps:

f * i (x, t + δt) = f i (x, t) + Ω f i (x, t) (2.
2)

|f i (x + δtξ i , t + δt) = |f * i (x, t + δt) (2.3)
in which, Eq. 2.2 applies the Ω operator to the current state f i (x, t) -a ket vector containing q PDFs of the lattice node. This reduces to local computations (also known as collision step), translated into floating-point arithmetic operations. Results of this sub-step are new PDFs of the next time step f * i (x, t + δt) , which temporarily remain within the local node. Eq. 2.3 then streams these PDFs into neighboring nodes (also known as streaming step), with the notation of spatial directions x + δtξ i . This sub-step is translated into memory load/store instructions.

Memory requirement

At each time step, the whole spatial domain must be updated before being able to start the next iteration. This spatio-temporal dependency of the LBM (shared by other stencil numerical schemes) compels developers, for the sake of code simplicity, to allocate two instances of the computational array, one as input of Eq. 2.2 (read-only) and one as output of Eq. 2.3 (write-only). This technique is usually known as two-lattice [START_REF] Massaioli | Achieving high performance in a LBM code using OpenMP[END_REF] (see the next section), whose the main drawback is the doubled memory consumption which significantly reduces the maximal reachable spatial resolution. It requires scientists to run their code on more machines with a larger aggregate memory space, thus resulting in larger cost and energy consumption.

From a programming point of view, LBM kernels are easy to implement and well-suited for parallelization on recent multi-/many-core platforms. However, lattice Boltzmann While recent architectures gain computing performance by increasing the clock speed and multiplying the number of cores, evolution of memory systems still cannot fetch enough data to keep cores busy. The dataset cannot always fit in caches and must be stored in the main (even remote) memory with much higher latency. The low arithmetic intensity of stencil kernels like LBM is thus the limit of performance, as well as their poor data-locality which reduces significantly the cache-reuse ratio. Previous studies in [START_REF] Mcintosh-Smith | On the performance portability of structured grid codes on many-core computer architectures[END_REF] and [START_REF] Obrecht | Performance Evaluation of an OpenCL Implementation of the Lattice Boltzmann Method on the Intel Xeon phi[END_REF] show that LBM implementations are memory-bound and hardly obtain good performance on CPU or Xeon Phi processors. GPU-based accelerators, thanks to their graphics-dedicated high-bandwidth memory, appear to be the most suitable platforms for LBM today.

Propagation algorithms

One-step two-lattice (OT)

One-step two-lattice (also known as two-lattice) is the most employed algorithm in LBM implementations on massively parallel architectures. The collision and streaming steps are fused into one compute kernel, either in pull or push scheme (see Fig. 2.3). This kernel loops on all lattice nodes and updates the whole domain at each time step. The two lattice arrays (A and B) which are swapped at the end of each time step, differ from each other by their access type within the compute kernel: one is read-only and one is write-only. Using non-temporal streaming store to perform the write operation in the streaming step is thus a good reason to use this algorithm. But this feature is not widely available on all architectures, due to its hardware cost.

One-step one-lattice (OO)

Different approaches, known as one-lattice algorithms, were introduced to reduce the memory footprint by working on only one lattice array and to improve data locality of the LBM. Most of them operate elaborate exchange of PDFs between neighboring lattice nodes in the parallel execution context. Pohl et al. [START_REF] Pohl | Optimization and profiling of the cache performance of parallel lattice Boltzmann codes[END_REF] proposed compressed-grid (also known as shift, see Fig. 2.4) approach to reduce the memory requirement of the two-lattice algorithm. With the same objective, Mattila et al. [START_REF] Mattila | An efficient swap algorithm for the lattice Boltzmann method[END_REF] developed a swap algorithm that requires almost half of memory space compared to the two-lattice algorithm. Comparisons of these algorithms with varying lattice-indexing and data layouts were carried out in [START_REF] Mattila | Comparison of implementations of the lattice-Boltzmann method[END_REF][START_REF] Wittmann | Comparison of different propagation steps for lattice Boltzmann methods[END_REF]. These studies show equivalent computational efficiency of compressed-grid and swap algorithm compared to the two-lattice approach, while consuming less memory. However, these two approaches both require definite iteration order and complex index calculations for shifting the two lattice grids (shift) or swapping distribution values between neighbors (swap). These dependencies make it very hard to implement shift and swap algorithms on highly parallel GPUs and accelerators in offloading mode (CUDA, OpenCL). Today, they are implemented only as sequential CPU code inside a subdomain and are scaled up by using MPI for inter-domain halo exchange [START_REF] Wittmann | Comparison of different propagation steps for lattice Boltzmann methods[END_REF]. This configuration yields satisfactory weak-scaling but cannot enable strong-scaling, since execution of each shared-memory subdomain cannot be parallelized by either OpenMP or Pthreads. Bailey et al. [START_REF] Bailey | Accelerating lattice Boltzmann fluid flow simulations using graphics processors[END_REF] presented the AA-pattern which overwrites read input PDFs by new collided data via two different kernels (even and odd time steps) (see Fig. 2.5). Geier and Schönherr [START_REF] Geier | Esoteric Twist: An Efficient in-Place Streaming Algorithmus for the Lattice Boltzmann Method on Massively Parallel Hardware[END_REF] introduced Esoteric twist (shortened to Esotwist) as an improvement of AA-pattern, by interacting only with neighboring nodes in positive xyz-direction. These two later algorithms work on one lattice array, are inherently asynchronous and thus are attractive for GPU and similar parallel architectures. However, they are more complex to implement than other algorithms mentioned above. Compressed-grid, AA-pattern and Esotwist need two kernels for even and odd time steps respectively. Esotwist can be implemented with one kernel, but requires imperatively the SoA (structure of arrays) storage layout to swap the Q pointers to their opposite direction after each collision, thus mostly only interesting for GPU architectures. locally in natural order. The even step performs reads of PDFs from neighboring nodes (pull), collides and writes back (push) to the same place on these nodes, with opposite PDFs.

Summary

The lattice Boltzmann method is one of the most memory-bound applications among other iterative stencil-based methods, such as image processing and computer vision.

Improving LBM performance lies on optimizing data-locality for better utilization of cache memories, as well as reducing the global memory footprint by inventing new implementation methods.

The fundamental challenge of any one-lattice algorithm is that memory accesses (read and write) must be performed carefully on the same lattice buffer to enforce the spatiotemporal dependency between nodes and time steps. More over, implementation often requires two versions of kernel code, adding more programming effort and reducing the maintainability of the application. On another hand, LBM boundary conditions on new physical models tend to be more and more elaborate. Typical LBM boundary conditions, such as simple bounce-back or interpolated bounce-back [START_REF] Obrecht | Efficient GPU implementation of the linearly interpolated bounce-back boundary condition[END_REF], imposes specific exchange rules of PDFs between adjacent nodes. Combining these conditions with existing one-lattice algorithms raises considerable complexity in implementation and validation, especially for 3D domains.

Other clustered many-core processors, despite a much lower global memory bandwidth with respect to GPUs, embed significant amount of fast local memories [START_REF] Dupont De Dinechin | A clustered manycore processor architecture for embedded and accelerated applications[END_REF][START_REF] Fu | The Sunway TaihuLight supercomputer: system and applications[END_REF]. They also provide more predictability in both computing time and data transfer. This enables using explicit and efficient user buffers for elaborate optimizations, such as software prefetching and streaming, based on local memories and asynchronous DMA engines.

Chapter 3

Basic Linear Algebra

Subprograms (BLAS)

Today, most software exists, not to solve a problem, but to interface with other software.

-I. O. Angell.

Background

Introduction

Since its first release in the 1980s, the Basic Linear Algebra Subprograms (BLAS) [START_REF] Chuck L Lawson | Basic linear algebra subprograms for Fortran usage[END_REF][START_REF] Jack J Dongarra | A set of level 3 basic linear algebra subprograms[END_REF] has been widely used as the de facto foundation of high-level dense linear algebra libraries such as Linear Algebra PACKage (LAPACK) [START_REF] Anderson | LAPACK: A portable linear algebra library for high-performance computers[END_REF] and Scalable LAPACK (ScaLA-PACK) [START_REF] Choi | ScaLAPACK: A scalable linear algebra library for distributed memory concurrent computers[END_REF], as well as in many inter-disciplinary scientific software. BLAS was designed to provide an unified and portable interface of numerical linear algebra operations on various computer architectures. Along other must-have software tools and libraries, BLAS is the first numerical API to be implemented and optimized on any architecture that targets high-performance computing.

The BLAS API1 defines three levels of numerical operations: (1) level-1 within or between scalar vectors, (2) level-2 between vector and matrix, and (3) level-3 between matrix and matrix. A typical example of a level-1 operation is the vector-vector addition (AXPY) (equivalent Triad in the STREAM benchmark [START_REF] John D Mccalpin | A survey of memory bandwidth and machine balance in current high performance computers[END_REF]). The level-2 performs, for instance, vector-matrix multiplication (GEMV), or the well-known solver of linear system of equations A • x = b (TRSV). For level-3, one should mention the General matrix

multiplication (GEMM) operation C ← α • A • B + β • C
, which is used as the core block of many computation-intensive benchmarks and applications. The High Performance Linpack (HPL) benchmark [START_REF] Petitet | HPL -A Portable Implementation of the High-Performance Linpack Benchmark for Distributed-Memory Computers[END_REF][START_REF] Jack J Dongarra | The LINPACK benchmark: past, present and future[END_REF] from which the Top500 list is built, as well as many other scientific applications are designed to map on this operation as much as possible to reach the maximal computing capacity of the target platform. In the latest years, GEMM has also become the kingpin of machine learning and deep learning advances.

Nevertheless, implementing and optimizing BLAS (typically level-3) on a given architecture has never been a trivial task. Straightforward implementations seldom deliver satisfying performance without advanced tiling and blocking techniques. To fill up the core pipeline, eliminate stall cycles and reach near-peak performance, developers compulsorily need to understand the low-level functionality of the hardware and write optimal assembly-level kernels. Multiplicity of BLAS parameters and their combinatorial cases yield up to several hundreds of assembly kernels to hand-tune and to maintain; in which to add extensions of instruction set architecture (ISA) and cache size evolution throughout processor generations. Developing and maintaining such an optimal library requires substantial time and expertise, that sometimes can only be afforded by the processor manufacturers or specialized research groups.

Conventional CPUs are shipped with proprietary libraries that implement BLAS, sparse BLAS routines, and BLAS-like extensions, such as: Intel Math Kernel Library (MKL) [START_REF]Intel. Intel Math Kernel Library[END_REF],

AMD Core Math Library (ACML) [START_REF]AMD Core Math Library[END_REF] and IBM Engineering and Scientific Subroutine Library (ESSL) [42]. Open-source options for BLAS-like functionality include the handoptimized GotoBLAS [START_REF] Goto | High-performance implementation of the level-3 BLAS[END_REF] [44] and its derivative OpenBLAS [START_REF] Xianyi | OpenBLAS[END_REF], auto-tuning solution such as Automatically Tuned Linear Algebra Software (ATLAS) [START_REF] Whaley | Automated empirical optimizations of software and the ATLAS project[END_REF], and projects that target modern multi-core processors such as Parallel Linear Algebra Software for Multicore Architectures (PLASMA) [START_REF] Yarkhan | Porting the plasma numerical library to the openmp standard[END_REF]. The BLAS-like Library Instantiation Software (BLIS) framework [START_REF] Field | BLIS: A framework for rapidly instantiating BLAS functionality[END_REF] is a recent development in the area of open-source and portable BLAS solutions for CPU-based architectures.

Vendors of GPU and other CPU accelerators also develop proprietary BLAS implementations adapted to their heterogeneous computing context: NVIDIA's CuBLAS [START_REF]GPU-accelerated standard BLAS library[END_REF],

AMD's clBLAS [START_REF]Advanced Micro Devices (AMD). A software library containing BLAS functions written in OpenCL[END_REF] and rocBLAS [START_REF]Next generation BLAS implementation for ROCm platform[END_REF]. Open-source projects for such architectures include MAGMA [START_REF] Tomov | Towards dense linear algebra for hybrid GPU accelerated manycore systems[END_REF], clBLAST [START_REF] Nugteren | CLBlast: A Tuned OpenCL BLAS Library[END_REF], and KAUST-BLAS [START_REF] Abdelfattah | KBLAS: An Optimized Library for Dense Matrix-Vector Multiplication on GPU Accelerators[END_REF]. These libraries, written in C-like languages (CUDA, OpenCL), rely on the vendor compiler and runtime API to generate executable code and offload computation onto the device. In order to abstract the hardware complexity and reduce programming effort, scheduling and memory management is hidden as much as possible to the developer and is managed by the deployment runtime and the device driver. Despite facilitating usage by non-expert users, the application portability and performance crucially depends on the vendor's or the open-source community's ability to implement and maintain an optimized BLAS library across multiple architectural generations.

General Matrix Multiplication (GEMM)

3.1.2.1 Basic implementation

C ← α • A • B + β • C.
Let us assume that matrices are square (m = n = k). The basic implementation of GEMM is based on the following the loop-based approach:

C i,j = α × n-1 t=0 A i,t × B t,j + β × C i,j i, j, t ∈ [0, n) (3.1)
Eq. 3.1 performs, per each C i,j element, 2n + 2 floating-point operations (flops),2 for a data traffic of 2n + 2 words. 3 Commonly, those floating-point operations can be performed within n+1 fused multiply-add (FMA) instructions. 4 This naive implementation has a poor arithmetic intensity. Spatial and temporal locality of data accesses in cache levels are not optimal either.

Blocked (Tiled) implementation

Goto et al. [START_REF] Goto | High-performance implementation of the level-3 BLAS[END_REF] revisited GEMM algorithms with multi-layer blocking (or tiling), in which each layer is corresponding to a memory level. For simplicity, let consider a two-layer configuration between a slow memory (DDR) and a fast memory (scratchpad or L1 cache). The matrix C on the slow memory is divided into

N × N blocks, each block is of size b × b (N = n b )
. The blocksize b is chosen so that the fast memory, whose size is S, can hold at least one block A and one block B (S ≥ 2b 2 words), or preferably one block C as well (S ≥ 3b 2 words). Let consider that three blocks A, B and C can fit into S, the communication cost between the slow and fast memory, the I/O traffic and arithmetic intensity (AI) of the blocked algorithm are written as follows:

I/O traffic = N 2 b 2 × N (read every block of A (N ) times) + N 2 b 2 × N (read every block of B (N ) times) + 2N 2 b 2 × 1 (read and write every block of C once) = 2N 3 b 2 + 2N 2 b 2 = 2n 3 b + 2n 2 = 2n 3 1 b + 1 n Complexity = 2n 3 AI = Complexity I/O traffic = nb n + b = b 1 + b n ≈ b (n b) (3.2)
As can be seen from Eq. 3.2, we can improve performance of the blocked algorithm between two successive memory levels by increasing b. Applying recursive tiling on contemporary CPU and GPU architectures is apparently the optimal approach for GEMM.

Furthermore, Goto et al. [START_REF] Goto | High-performance implementation of the level-3 BLAS[END_REF] also proposed an additional packing step before inner computations, which consists in reordering A and B data blocks into a pre-defined contiguous layout,5 in order to maximize cache hit ratio and minimize penalty of Translation Lookaside Buffer (TLB) misses. Through the packing step, the implementation will also be able to handle multiple parameter combinations (transa, transb, uplo, sidea etc.) with only several well-defined inner micro-kernels. This allows implementing BLAS functions without tuning hundreds of assembly kernels, produces well-structured and highly maintainable code, whose a successful example is the BLIS framework.

BLAS-like Library Instantiation Software (BLIS)

BLIS, stands for BLAS-like Library Instantiation Software [START_REF] Field | BLIS: A framework for rapidly instantiating BLAS functionality[END_REF], is a sub-project of libflame [START_REF] Van Zee | Introducing: The libflame library for dense matrix computations[END_REF]. Inspired from GotoBLAS [START_REF] Goto | High-performance implementation of the level-3 BLAS[END_REF], BLIS is designed with fundamental principles in dense linear algebra, including incremental and recursive construction of BLAS operations for code-reusability, data-packing and cache/register blocking for optimal locality, as well as ability to integrate platform-specific assembly code. rank-k update and constitutes the sixth loop. Other level-3 operations are then built on top of GEMM and this partioning method.

m R m R 1 += += += += += += n C n C k C k C m C m C 1 n R k C n R Pack A i → A i ~Pack B p → B p ñR A B j C j A p A i B p C j A i ~Bp B p Ci C i k C L3 cache
Main advantages of BLIS can be mentioned as follow:

• BLIS introduces a reduced set of micro-kernels (gemm, trsm, gemmtrsm), written in the portable C99 standard by default, as a reference implementation.

• BLIS defines a reduced set of execution parameters (cache sizes, memory alignment, memory allocator), largely used by the internal algorithms, but can differ significantly from one architecture to another.

• BLIS provides a user-defined configuration header which allows arbitrary modification of these execution parameters and easy plug-and-play of ISA-specific micro-kernels to generate a nearly optimized BLAS library on any cache-based architecture, without touching the core functions of BLIS. This abstract design enables a custom fit on any cache-based system implementing any instruction set with straightforward portability and high-performance. Cache-blocking and packing implementation techniques in BLIS have been proved to be analytically optimal [START_REF] Tze | Analytical modeling is enough for high-performance BLIS[END_REF] on multi-core and NUMA memory hierarchies, delivering competitive performance to other vendor libraries [START_REF] Field G Van Zee | The BLIS framework: Experiments in portability[END_REF].

Summary

New multi-and many-core architectures keep appearing and are moving fast. On these systems, writing library and software are becoming more and more challenging. Moreover, code-portability sometimes counters performance, due to the hardware diversity.

The BLIS framework has emerged as a promising solution for instantiating a light-weight and high-performance BLAS library for cache-based architectures.

However, both HPC design strategy and modern embedded and intelligent computing are coming up with more and more power-efficient and non-conventional architectures, on which, writing a BLAS library in pace with the hardware represents a big software challenge. They often do not have a hardware-assisted cache coherency, whereas a software-based cache protocol would suffer significant overhead. Hardware prefetcher, out-of-order execution and advanced branch prediction are commonly discarded to reduce power consumption (and fortunately avoid security vulnerabilities 6 ).

On those architectures, computation-intensive parts of code are expected to perform software-managed data-prefetching, by leveraging Direct Memory Access (DMA) engines and operating on scratchpad memories, instead of the traditional cache-based load/store scheme. Support of the asynchronous programming model based on Remote-DMA (RDMA), considered the key enabler of performance on DMA-based architectures, is currently the missing point of any BLAS-like library on these later platforms. -Arthur C. Clarke.

Introduction

Company

Kalray is a fabless semiconductor company, founded in France in 2008 after a spin-off from the CEA (French Alternative Energies and Atomic Energy Commission). Kalray is specialized and pioneering in developing a new family of many-core processors, namely Massively Parallel Processor Array (MPPA). The MPPA architecture offers unique parallel computing capacity, low latency and low-power consumption. The Very Long Instruction Word (VLIW) core architecture, distributed non-coherent memory system and dual control-plane and data-plane Network-on-chip enable time-predictability necessary to embedded, mixed-critical and real-time systems. The MPPA massively parallel architecture is also suitable for modern and energy-efficient HPC workloads, typically in the area of image processing, computer vision and autonomous vehicles.

MPPA architecture overview

The second generation of Kalray MPPA-256 processor, codenamed Bostan (see Fig.

Distributed-memory POSIX-C

In the distributed POSIX programming model, an MPPA2-256 processor is exposed as a distributed multi-process system. Each compute cluster is considered as an individual computing unit with a main function and a proper memory space (SMEM). Each main The MPPA asynchronous library (namely mppa-async) implements DMA-based asynchronous one-sided put/get, remote atomic operations, peek, poke and two-sided queues.

These features enable high-throughput data-plane communication, and low-latency controlplane signal and synchronization. This programming model requires high design and development efforts, but can offer parallel applications great performance improvement over conventional platforms.

Host-based OpenCL acceleration

The Kalray SDK also allows programming MPPA as a compute-accelerator based on the OpenCL 1.2 Data-Parallel model [8]. A software-based cache protocol, so called Distributed Shared Memory (DSM), is used to implement a globally coherent cache between the 256 cores, on which OpenCL kernels are deployed. As part of contributions of this thesis, an extended set of OpenCL asynchronous primitives, including general strided, 2D and 3D copy, was introduced and implemented in the Kalray OpenCL toolchain and is discussed in Chapter 8 on page 85.

Summary

The Kalray MPPA processor family is a highly parallel computing platform. Programming the MPPA processor requires careful segmentation of data and explicit transfers onto local memories. This section closes the introduction part of the dissertation. In this opening segment, we have presented the thesis context, the challenges of current and future HPC trends and the potential as well as problematics of clustered many-core architectures. Application focus was given to (1) the lattice Boltzmann method as a memory-bound context with high bandwidth requirement, and (2) the BLAS library as a compute-bound scenario (level-3) with its intrinsic complexity and portability issues.

In the following chapters, we present approaches and contributions to address these two typical fields of the HPC domain, revised and shaped for the clustered DMA-based many-core architectures. We also implement, demonstrate and report analysis of each of those solutions on the Kalray MPPA2-256 processor, as well as other mainstream computing architectures whenever applicable.

Contributions

Chapter 5

Optimizing 3D LBM on Many-core Processors However difficult life may seem, there is always something you can do and succeed at.

-Stephen Hawking.

Introduction

As discussed in Chapter 2, Section 2.1.2 (page 13), LBM implementations are memorybound and hardly obtain high performance on CPU or Xeon Phi processors. GPUbased accelerators, thanks to their graphics-dedicated high-bandwidth memory, appear to be the most suitable platforms for LBM today. However, their low capacity of local memory prevents from using optimization techniques for data prefetching (to reduce transfer time) and data sharing between cores (for stencil neighboring dependencies).

Although other clustered many-core processors have much less global memory bandwidth than GPUs do, they embed significant amount of fast local memory, see [START_REF] Dupont De Dinechin | A clustered manycore processor architecture for embedded and accelerated applications[END_REF] and [START_REF] Fu | The Sunway TaihuLight supercomputer: system and applications[END_REF], and provide more predictability in both computing time and data transfer. This enables using explicit and efficient user buffers for elaborate optimizations, such as software prefetching and streaming. This motivates our approach in developing a pipelined 3D LBM algorithm on the Kalray MPPA processor, based on local memory exploitation and asynchronous communications. Our algorithm is described in every detail and can be used on similar many-core architectures.

Our key contributions are as follows:

1. Introduction of a new parallel algorithm for decomposing and streaming 3D stencil domains on local-memory-centric clustered many-core processors, by user-buffers and asynchronous software-prefetching to build a pipelined 3D stencil kernel. The proposed approach is implemented from the LBM compute kernel of OPAL [START_REF] Obrecht | Performance Evaluation of an OpenCL Implementation of the Lattice Boltzmann Method on the Intel Xeon phi[END_REF] and delivers 33% performance gain compared to its original OpenCL code on the Kalray MPPA-256 Bostan many-core processor.

2. This work provides fundamental responses and methods to further domain decomposition algorithms on clustered many-core processors (2D/3D stencils, image processing). An API proposal is also given in designing simple 2D/3D asynchronous copy functions on DMA-based platforms.

Detailed description of the use of generic equations to calculate decomposition

indexes dynamically, subdomain dimension and halo size, usable with or without ghost layer as in [START_REF] Mattila | Comparison of implementations of the lattice-Boltzmann method[END_REF].

The remainder of this chapter is structured as follows. Section 5.2 presents some related works that are relevant for our contributions. Section 5.3 introduces some low-level asynchronous transfer primitives required for building 3D stencils streaming algorithm on the MPPA processor. Section 5.4 presents an overview and technical details of the new LBM streaming algorithm using these asynchronous transfers. Experimental results are presented in Section 5.5, and we conclude in Section 5.6.

Related work

The straightforward method for implementing LBM is to use two instances of the lattice grid. Collision is carried out on data read from the first grid and propagation consists in writing the new distribution values to the second one. At the next time step, the two grids are swapped and the same procedure is repeated. One-step two-lattice method with collision and propagation fused in a same kernel was first introduced by Massaioli and Amati [START_REF] Massaioli | Achieving high performance in a LBM code using OpenMP[END_REF]. In the fused kernel, propagation can be done either before (pull scheme) or after collision (push scheme). In spite of its implementation simplicity, the two-lattice method results in substantial memory allocations with large domains.

Most existing LBM implementations on GPU employ the fused two-lattice approach as the easiest and most computationally efficient method. In particular, OpenCL Processor Array LBM (OPAL) from [START_REF] Obrecht | Performance Evaluation of an OpenCL Implementation of the Lattice Boltzmann Method on the Intel Xeon phi[END_REF] implements a one-step two-lattice 3D LBM solver based on the D3Q19 stencil. OPAL is designed to be simple and portable on GPUs, accelerators and other OpenCL-enabled devices.

In the related work on porting a 3D seismic wave propagation on the MPPA processor, Castro et al. [START_REF] Castro | Energy efficient seismic wave propagation simulation on a low-power manycore processor[END_REF] developed a 2D-prefetching algorithm for anticipating data transfers between global memory and local memory. The 3D domain is decomposed in small 2D slices. These slices are copied to the local memory such that transfers overlap with computations. The authors observed important waiting time for data arrival without identifying clearly the DDR bandwidth limitation of the MPPA. The impingement of halo slices on data throughput was not studied either.

Raase and Nordstrom [START_REF] Raase | On the use of a many-core processor for computational fluid dynamics simulations[END_REF] presented a 2D and 3D LBM implementation on Epiphany, a clustered many-core architecture very similar to MPPA. The LBM domain is distributed on 16 cores with user local memory of 24KB per core. Subdomain distribution is done by static mapping of a 4x4 topology on the 16 cores. This 2D mapping is also used on the 3D problem where the third dimension of subdomains is assigned with one global domain dimension, giving rectangular parallelepiped subdomains. These choices allow simulating only very small problem sizes (e.g. 12 × 35 × 12) and cannot be scaled to large simulations. The authors declined using the DRAM memory to implement a streaming algorithm for large LBM domains. Neither memory bandwidth optimization nor possibility of using DMA to perform asynchronous transfer on Epiphany was discussed. In this work, we aim to provide a generic and scalable 3D decomposition with its cuboid distribution function and asynchronous subdomain streaming to reduce data transfer time. Such algorithm can be used as a reference point to implement further high performance LBM or stencil applications on clustered many-core architectures.

Nagar et al. [START_REF] Nagar | LBM-IB: A parallel library to solve 3D fluid-structure interaction problems on manycore systems[END_REF] implemented a similar cube-based decomposition and distribution function which maps on CPU threads in the shared-memory context of large-memory multi-socket systems. Halo exchange between threads is done by writing directly into the memory zone of the respective cube owners, protected by mutual locks thanks to the CPU cache system. This mechanism cannot be directly used onto clustered architectures like MPPA as it requires either: (1) explicit inter-cluster communications; or [START_REF] Hascoët | Asynchronous one-sided communications and synchronizations for a clustered manycore processor[END_REF] committing changes to the global memory then fetched by other clusters. Solution (1) is not relevant in our scope due to (small) local memories, numerous subdomains must be streamed continuously in the MPPA's compute clusters. Such streaming should be done preferably by a self-governing and synchronization-free algorithm. Thus, keeping data in the local memory and waiting for communication does not seem appropriate, not to mention the complexity of managing the inter-subdomain spatial data dependency. In this work, we choose to adapt solution (2), consisting in continuously committing changes of subdomains to the global memory and performing one global synchronization between clusters at each simulation time step.

To the best of our knowledge, there is no work yet on solving the challenges of simulating large LBM domains on clustered many-core architectures. However, large LBM domains cannot fit into on-chip memory and must be stored in the off-chip DDR memory, which has much higher latency. Hence, using DMA to perform asynchronous transfers between off-chip and on-chip memories becomes a key performance factor in order to mask the memory latency. This involves important code re-structuration, as well as new communication primitives and algorithms. In this work, all these problems are addressed and solved while keeping a clear abstraction level from the underlying target hardware for the sake of genericity.

Low-level 3D asynchronous API

In this section, we briefly present some essential primitives performing asynchronous 3D data transfers used to build our pipelined LBM algorithm onto the Kalray MPPA processor. As shown in Fig. 5.1, the mppa async point3d t type describes copy-position and dimensions of the global and local 3D buffers. The subdomain is represented by width×height×depth elements, the size of each element in bytes being denoted size. We take an example to illustrate this specification design. In a common image processing decomposition, one may need to copy a 2D sub-image of 16 × 16 pixels to a larger local buffer, allocated at 18 × 18 pixels for instance. In this case, one must deal with a local stride of two pixels between each data block. This is important when local buffers are declared as true multi-dimensional arrays in the C99 standard, a feature which particularly eases 2D and 3D stencil programming. With the convenient mppa async point[2|3]d t data type (see Fig. 5.2), arbitrary positions and copy-block dimensions are automatically taken into account inside the 2D/3D put and get functions, facilitating subdomain copy and computation.

A structure mppa async event t is also defined in the API to contain required information for performing an asynchronous transfer. In a put/get function, if the event structure is set, the function fills a pending transaction event and returns immediately (non-blocking paradigm). One can further come back and wait on this event by calling the mppa async event wait() function for job completion. Otherwise, when the event structure is NULL, the function blocks and returns whether the buffer is ready to be reused (put) or the data are received (get). In the following, we take the D3Q19 LBM kernel from OPAL [START_REF] Obrecht | Performance Evaluation of an OpenCL Implementation of the Lattice Boltzmann Method on the Intel Xeon phi[END_REF] as a reference point, from which we propose a generic 3D LBM streaming algorithm with domain decomposition, detailed index and halo size calculation in any configuration. The streaming method is used for updating the whole domain by one time step, then is repeated till the end of simulation duration. While we are focusing on optimizing LBM, our streaming method can also be generalized for other kinds of stencil codes, by adapting the compute kernel, and a suitable set of asynchronous transfer primitives (2D/3D).

The first step consists in re-writing the LBM kernel of OPAL from OpenCL-C to a standard C99 code to run on CCs. Given the similarity between OpenCL-C and standard C99, the porting process did not raise much difficulty. The one-step two-lattice method with pull scheme originally implemented in OPAL is re-applied. Two instances of the 3D lattice grid (LatticeEven, LatticeOdd ), each containing a macro-pipeline using asynchronous 3D put and get functions to overlap computation and communication. We also apply the two-lattice method on local memory, i.e. the number of buffer slots is doubled, one for fetching the pre-collision cuboid (S) from the first global lattice grid and one for storing the post-collision cuboid (S ) that will be put in the second lattice grid. The pre-collision cuboid is allocated for

L x × L y × L z nodes,
F x × F y × F z nodes,
while the post-collision cuboid only needs to store C x × C y × C z nodes. Fig. 5.3 only draws one global lattice grid for compactness, but it should be understood that the local post-collision cuboid will be put in the second grid. These two global grids are then swapped before starting the processing of the next time step.

Ideally, the algorithm should run on multiple compute clusters and exploit all processing cores (PEs) in each cluster (multi-cluster multi-PE). For instance, on MPPA, multithreading within a compute cluster is enabled by spawning up to 15 threads, one per PE, from the PE0 in the Pthreads fashion (create, join). As there are 16 compute clusters available on MPPA, each CC is then responsible for M 16 subdomains. Note that depending on the value of M , there might be K trailing subdomains (K ∈ [0..15]).

If K > 0, the algorithm must perform an extra step to copy, update and put back these K trailing subdomains by K compute clusters, while other clusters are waiting. A synchronization barrier at the end of each time step is needed between all CCs to avoid data races at the next time step. This procedure is then repeated as many times as the number of timesteps.

The double-buffering (2-depth) pipeline in Algorithm 1 is the most basic algorithm where communication is overlapped by only one compute-step. As computations are faster than data transfers, deeper pipelines such as triple-or quadruple-buffering (whose details are found in Fig. 5.4 ) provide better overlapping, but also require more local memory. Note that the time spent in GET and PUT is considered negligible (non-blocking) and transfers are executed in background. However, the time spent in COMPUTE depends Algorithm 1 Explicit macro-pipeline of 3D stencil updates using double-buffering within a time step. on core speed, while the WAIT time depends on how fast the memory system is serving transfer requests and how they are hidden entirely or partially by the COMPUTE function.

In the next sub-sections, we propose methods to solve the following questions that immediately arise from Algorithm 1:

• How can we distribute fairly and exclusively all subdomains across CCs with their proper subdomain-indexing?

• Which subdomain size and pipeline depth should we choose to fit with the local memory size and to obtain the best trade-off?

• How to manage copy indexes and halo size of any subdomain, with or without using ghost layer?

m=0 1 2 3 4 5 6 7 i=0 1 2 0 1 2 0 1 buffers[0] G WCP WG WCP WG WCP W buffers[1] G WCP WG WCP WG WCP W buffers[2] G WCP WG WCP W Prolog Epilog

Subdomain distribution

Given a CC identified by cc id ∈ [0..15], its working subdomains is indexed by a onedimensional range m as cc id × M 16 ≤ m < (cc id + 1) × M 16 (assume K = 0). Mapping bijectively this 1D domain (m) to a 3D one (m x , m y , m z ) for spatial cube indexing (see Fig. 5.3) was done by space filling curves, such as Morton or Hilbert curves, in [START_REF] Sagan | Space-filling curves[END_REF]. These curves have been efficiently implemented by bit-interleaving in [START_REF] David S Wise | Ahnentafel indexing into Morton-ordered arrays, or matrix locality for free[END_REF] or lookup-table in [START_REF] Thiyagalingam | Is Morton layout competitive for large two-dimensional arrays yet?[END_REF]. However, Morton, Hilbert and other curves are better suited for square or cubic grids where the number of elements in all dimensions is equal. In our 3D decomposition scheme, despite the fact that global lattice domain may be cubic (L x = L y = L z ), subdomains may be not (C x = C y = C z ) due to many reasons (see next section), thus these curves are not always suitable for subdomains, as M d can be different.

In order to solve this problem, we implement a simple alternative bijective function f : N → N 3 in Fig. 5.5. It follows the 3D-row-major layout, which is also a space-filling curve, to index subdomains. Each conversion of the 3D-row-major curve implemented by f takes less than 10 instructions and is as fast as Morton or Hilbert curves.

Local subdomain dimensions

In most of the cases, a cubic subdomain would be ideal for coding and optimizing.

However, the local memory of clustered many-core processors is usually limited but plays an important role. On each MPPA's compute cluster, 2 MB local memory is quite small and should also host an embedded operating system, services and the user application binary. A remaining space of about 1.5 MB is available for dynamic buffer allocations. Some auxiliary variables are also needed in LBM for macroscopic monitoring (velocity, density. . . ). The maximal allocatable space for local pre-collision and post-collision cuboids is around 1.4 MB. Halo copy also consumes memory bandwidth. Hereafter, we refer to halo bandwidth HBW as the bandwidth lost in fetching halo layers. The HBW ratio is defined as the quotient of the number of halo cells by the total number of copied cells (main block and halo). On small subdomains, this ratio can be significant. For example, given a cubic subdomain whose main block size is C x × C x × C x , its HBW ratio is calculated by the below formula and represented as in Fig. 5.6.

g(C x , C x , C x ) = (C x + 2) 3 -C 3 x (C x + 2) 3 (5.1)
Cube size (Cx = Cy = Cz) Halo bandwidth ratio

(C x ) = (C x + 2) 3 -C x 3 (C x + 2) 3
Find: (D, C x , C y , C z ) Maximize: C x × C y × C z (nodes updated per subdomain) Minimize: F x × F y × F z (per subdomain storage)
Minimize:

(F x ×F y ×F z )-(C x ×C y ×C z ) F x ×F y ×F z (HBW) Subject to: D × ((F x × F y × F z ) + (C x × C y × C z )) × 19 × 4 1024 2 ≤ 1.4 F d = C d + 2 ; C d ∈ {2 n } ; D, n ∈ N + (5.2)
For instance, using D ≥ 3 in order to have better overlapping than with a 2-depth HBW ratio. On the other hand, note that increasing pipeline depth is not relevant, because the higher D is, the smaller (C x , C y , C z ) will be, thus HBW will become unacceptable. Moreover, compute cores will switch between small subdomains more often.

The accumulated waiting time will also be more important due to the exponential number of DMA requests and processing overhead of the DDR asynchronous services.

Local and remote copy-index management

In this section, we present generic analytic formulae to process dynamically copy indexing, subdomain size computation and halo cutoff management depending on geometric position of the subdomain. Adding a ghost layer surrounding the computational domain is a common technique to simplify the implementation of the streaming step at boundary cells, see e.g. [START_REF] Mattila | Comparison of implementations of the lattice-Boltzmann method[END_REF]. However, we choose not to use this approach in our work, mainly to minimize global memory allocations and avoid wasting bandwidth/storage in moving ghost cells.

However, in our 3D decomposition algorithm, this decision requires careful calculation of copy parameters from subdomain indexes. It is important to note that as the pre-collision cuboid S embeds two additional halo layers for each dimension (F d ), its computational space begins at (1, 1, 1) and ends at (F x -2, F y -2, F z -2) included. When fetching a non-boundary subdomain (main block + halo) from global memory to S, the arrival point of data at the local buffer is set to (0, 0, 0), and the remote point is computed as the global beginning position of the subdomain minus one (back-off) in each dimension

((m d × C d ) -1).
As ghost layers are not used in our implementation, a boundary subdomain can have up to three missing sides, depending on its location (see Fig. 5.7). Consequently, the halo layer of these missing sides needs to be pruned from the copied cuboid. The remote read-point and local write-point must also be adjusted as well. In order to generalize the solution, we introduce here three parameters associated respectively to these three adjustments: halo cutoff, remote offset and local offset.

We present in the following, generic formulae which determines copied positions and halo cutoffs of a given 3D cuboid subdomain (m x , m y , m z ), generalized from the 2D representation of Fig. 5.7.

const A = (A x , A y , A z ) = (0, 0, 0) R = (R x , R y , R z ) = (m x × C x , m y × C y , m z × C z ) B ad = A d + local offset(m d , M d ) B rd = R d + remote offset(m d , M d ) S d = F d + halo cutoff (m d , M d ) (5.
3)

The point A = (A x , A y , A z ) = (0, 0, 0) is the start point of the local buffer. The local offset remote offset halo cutoff

point R = (R x , R y , R z ) = (m x × C x , m y × C y , m z × C z )
(from point A) (from point R) (from F d ) m d = 0 1 0 -1 0 <m d && m d < M d -1 0 -1 0
m d = M d -1 0 -1 -1
These position computations can also be applied on other implementations which use ghost layer, by setting all remote offset to -1, i.e. allowing to jump out of the computational domain, and all local offset, halo cutoff to zero, i.e. imposing to copy extended subdomain F d to A d , instead of copying S d to B ad . As one cannotice in Fig. 5.8, the OPAL async algorithm outperforms the OpenCL version by more than 30% on the single-DDR mode (from 12 MLUPS to 16 ± 1 MLUPS).

Results and discussions

We also see that the configuration with less HBW (3-depth, 36% HBW) delivers higher performance than the 4-depth configuration (43% HBW). While consuming memory bandwidth, halo cells are copied because of the read-dependency between neighbors.

This does not contribute to the final performance. Fig. 5.8 shows that the less memory bandwidth halo cells take up, the more performance we obtain. This leads to think that the HBW of 2D/3D stencil computations aimed to reach Exascale, like weather forecast, ocean simulation and CFD, should be lessened on future clustered many-core processors.

For this to happen, these many-core chips should embed bigger local memory on each compute unit to tear down the useless part of halo exchange due to domain decomposition. Finally, Fig. 5.8 also shows the expected 2× performance speedup by using two DDRs compared to the single-DDR mode.

Performance extrapolation

For a better understanding of the benefit of our streaming algorithm, we modified the OPAL async code to be able to work with arbitrary values of pipeline-depth. Different pipeline depths were then tried out [START_REF] Ho | MPI communication on MPPA many-core NoC: design, modeling and performance issues[END_REF][START_REF] Hascoët | Asynchronous one-sided communications and synchronizations for a clustered manycore processor[END_REF][START_REF] Ho | New parallel in-place update algorithm for better memory usage in 3D lattice Boltzmann algorithm[END_REF]6,8) to see if increasing the number of asynchronous buffers can improve the performance. The block size is thus reduced to 8 × 8 × 8 so that up to eight subdomains can be stored in the local memory. Moreover, instead of using all the 16 compute clusters, we now vary this number of clusters and set the domain size to 128 3 to study the strong scalability of the algorithm. We consider using only the double-DDR mode this time to obtain the best performance.

In Fig. 5.9, as expected, the 1-depth code (blue line) is slower than other version with communication-computation overlapping. However, we obtain exactly the same performance as the double-buffering case when using more than two buffers [START_REF] Ho | New parallel in-place update algorithm for better memory usage in 3D lattice Boltzmann algorithm[END_REF]6,8). The performance line scales from 1 cluster to 8 clusters, then reaches almost a stable value of between 20 to 22 MLUPS from 8 clusters to 16 clusters. To explain this, we added the sustained throughput of 3D transfer (red line) from the Kalray unit test dedicated to 3D asynchronous copy. This test only does some ping-pong copies to the DDR and does not perform any calculation (Arithmetic Intensity (AI) = 0 flops/byte). We observe that the native 3D copy reaches the maximum throughput with as few as four clusters (6 GB/s), then remains the same for higher numbers of clusters (which is the same trend as the performance of OPAL async.). Four clusters are thus enough to saturate the DDR bandwidth. Unlike the 3D unit test, our LBM code performs real computation on the copied data. Its AI is about 350/(2 × 19 × 4) = 2.3 flops/byte, which means that each CC spends more time working on a 3D data block. This explains in Fig. 5.9 the MLUPS performance which reaches its upper bound for 8 clusters, instead of 4 clusters of the 3D unit test.

Another precise way to interpret the performance of 22 MLUPS is to apply the performance estimation formula presented by McIntosh-Smith et al. [START_REF] Mcintosh-Smith | On the performance portability of structured grid codes on many-core computer architectures[END_REF]:

P = B × 10 9 19 × 2 × 4 × 10 6 (MLUPS) (5.4)
in which B is the effective memory bandwidth in GB/s. In order to take into account the additional cost of halo copy in our decomposition algorithm, we multiply P by

(1 -HBW ), the effective part of bandwidth (main block) which generates the real performance:

P h = 6.0 × 10 9 19 × 2 × 4 × 10 6 × 8 3 10 3 = 20.2 MLUPS (5.5)
This estimation P h , shows that there is seemingly a little performance gain to perform asynchronous transfers on clustered many-core processors (here MPPA as an example) as for today. This is not because the streaming algorithm is not good, but because the overlapping gain time is too small compared to the lengthy waiting time for data due to the DDR3 bottleneck. This also demonstrates the memory-bound property of general stencil computations and leads to think that newer memory technologies, such as DDR4 and others, will be a performance boost on these architectures.

Notice that the scale-down of the 3D throughput versus the peak 17GB/s of two DDRs is caused by the fact that strided copies (2D/3D) must read data from a lot of different DDR memory banks. Furthermore, these copies can unavoidably suffer bad alignments due to the access pattern of application (Q = 19 floats), thus bear an efficiency factor of 3D transfer compared to the linear copy. For instance, on the current MPPA Bostan platform, if the linear transfer factor is normalized at 1, the 3D factor lies often in between 0.35 and 0.42, depending on the copy layout (size of each contiguous block, alignment of strides and dataset).

A correlation, computed by the lm function in R, from 1 to 8 clusters gives the performance expectation of our streaming algorithm if we were not bounded by the memory bandwidth (gray line). These results confirm that our pipelined LBM algorithm is strongly scalable, but is quickly memory-bound on MPPA and that its performance heavily depends on the hardware memory bandwidth. Our results also show that the imbalance between computing power and data throughput is one of the largest drawbacks of actual clustered many-core processors, and demonstrate the interest of future high-bandwidth memory technologies.

Conclusions

We introduce a decomposition approach for generic 3D stencil problems with formulations for calculating dynamically copied position indexes, subdomain addresses, subdomain size and halo cells. These analytic expressions are valid with or without ghost layers and are also usable for 2D problems. Based on this decomposition, our new pipelined 3D LBM code outperforms the original OpenCL version by 33 %, by overlapping computation and communication.

We expect that anticipating data requests by asynchronous memory transfers would improve effective throughput and that we could overcome the memory bound of the studied LBM kernel, by introducing enough pipeline depth to hide the global memory access latency. In practice, performance results are still bound by memory bandwidth and increasing the number of buffers (pipeline depth) does not improve performance, as the DDR3 memory is already fully loaded. Moreover, reducing subdomain size to increase pipeline depth induces significant bandwidth consumption for halo copy. Furthermore, the impact of HBW on small local memories was also identified as a governing factor of performance in our algorithm. We find out that the best strategy is to have cubic subdomains as large as possible and that the double-buffering scheme is enough on the current generation of MPPA processor. We furthermore presented comprehensive linear-programming equations which give the best trade-off between these structuring parameters.

In the next chapter, we study a new LBM propagation method which performs inplace lattice update (one-step one-lattice). Such a method will reduce by half the local memory requirement, thus increase the subdomain size, trim down halo bandwidth and improve performance. Porting async work group copy {2D|3D} primitives to the next OpenCL specification is also under consideration, as this would considerably improve the exploitation of local memory on clustered many-core processors.

Chapter 6

In-place LBM Propagation

Algorithms

The biggest difference between time and space is that you can't reuse time.

-Merrick Furst.

Introduction

Besides the two-lattice propagation algorithm, several one-lattice algorithms were introduced to reduce the memory footprint (by working on only one lattice array) and also to improve data locality of the LBM. Most of them operate elaborate exchange of PDFs between neighboring lattice nodes in the parallel execution context. The fundamental challenge of any one-lattice algorithm is that memory accesses (read and write) must be performed carefully on the same lattice buffer to enforce the spatio-temporal dependency between nodes and time steps. More precisely, in a parallel algorithm, new PDFs of a node after collision should not be written directly into memory without special care about the old data that may have not been used yet by neighboring nodes.

On another hand, LBM boundary conditions on new physical models tend to be more and more elaborate in terms of interaction between solid and fluid cells. Typical LBM boundary conditions, such as simple bounce-back or interpolated bounce-back [START_REF] Obrecht | Efficient GPU implementation of the linearly interpolated bounce-back boundary condition[END_REF], imposes specific exchange rules of PDFs between adjacent nodes. Combining these conditions with existing one-lattice algorithms such as compressed-grid, swap, AA-pattern or Esoteric twist (see the next section), raises considerable complexity in implementation and code assessment, especially for 3D domains.

Our key contributions in this chapter are as follows:

1. Two novel algorithms, two-wall and three-wall, in the one-lattice LBM class and their implementation detail in OpenMP for shared memory context and OpenCL for heterogeneous memory systems.

2. Comparison and performance analysis of these two-/three-wall algorithms versus the AA-pattern (one-lattice) and the state-of-the-art two-lattice algorithm on current mainstream computing platforms: Intel Xeon NUMA CPU, Intel Xeon Phi Knights Landing MIC and NVIDIA Tesla Pascal GPU.

3. Promising results when using the proposed algorithms to implement complex LBM and CFD problems at high spatial resolution on current and future many-core processors.

The remainder of this chapter is structured as follows. Section 6.2 describes our approach, advantage and limitations of the two-wall algorithm. Section 6.3 presents the three-wall algorithm as an improvement of two-wall. Implementation details in OpenMP and OpenCL are given in Section 6.4. Experimental results and discussion on CPU, Xeon Phi KNL and GPU Pascal are given in Section 6.5 and we conclude in Section 6.6.

Before updating the wall z (a 2D YX-wall) from the time step t to t+1, the current wall is used for saving the pre-collision data of this wall before it is modified by the in-place post-collision propagation. Idem, the past wall contains the pre-collision data of the back wall z -1, as this wall in G has been updated to the time step t + 1 in the previous z-iteration, thus cannot be used here for the wall z at the time step t. Inside each wall, lattice nodes can be independently copied into the current wall. Then a memory-fence barrier (if parallel execution) is needed to ensure that all copies are done before performing the collision and the in-place propagation on the wall z of G. At the end of each wall-update, all working threads must perform another barrier to achieve global synchronization and to avoid data-races between concurrent threads. Pointers to current wall and past wall are then swapped before advancing in the next wall in zdirection. The current wall of the previous iteration z -1 now becomes the past wall of the current iteration z. Inversely, the past wall of the previous iteration z -1 is now the current wall of the current iteration z and can be recycled for storing the precollision state of the current wall z, before being overwritten by the upcoming collision.

Algorithm details are presented in Fig. 6.2 in which the D2Q9 lattice is used. From the Fig. 6.2, parallelism is therefore possible at two steps: (a) and merge of (b)+(c). These two steps are separated by a barrier to ensure the read-after-write requirement of the copy operation of current wall, since the collision will read data from the past wall and current wall for neighbors in z -1 and z walls. This easily yields an OpenMP implementation by two #pragma omp parallel for statements or an OpenCL one, consisting in two kernels (one for (a) and one for (b)+(c)). They will be presented later in Section 6.4.

Advantages and limitations

The main advantage of the two-wall algorithm is that adapting complex boundary conditions is trivial by replacing any access to z -1 and z cells by past wall and current wall, respectively. This means that there is only one kernel code to maintain, instead of two on AA-pattern or Esoteric twist (even and odd). Adding and testing new boundary conditions will require less effort. Data locality is also a good point. Memory accesses from different cores have more chance to share and hit in cache (2D walls) than other algorithms, where PDFs are read from or written to neighbors in three dimensions, thus increasing penalty of cache-miss and cache-thrashing.

In return, there are two performance challenges of the two-wall algorithm. The first is the data movement of 4q (see Table 6.1) per node per time step versus 2q or 3q of two-lattice and others. This difference can be visible on a bandwidth-sensitive platform with non-persistent cache between kernel executions (like OpenCL on GPU). Secondly, updating a wall requires two barriers of all threads, resulting a total of 2L z barrier count for a grid G of L z walls. Multiplied to the number of time step T , the algorithm faces a very high number of barriers. As a matter of fact, two-wall requires fast barrier support from the hardware and/or an efficient event-driven programming language to launch and synchronize successive kernels at least cost, ideally via persistent threads and asynchronous runtime in background. Subsequently, in the next section, we introduce three-wall, an extended version of two-wall which breaks down the read-after-write dependency within each wall and reduces the number of barriers from 2L z to L z .

Three-wall propagation algorithm

In two-wall, the memory-fence barrier after the wall-copy is mandatory because the current wall is needed in the following collision step (read-after-write). This leads to think that if the wall z has been copied beforehand (for example at the iteration z -1), the three steps (a) (b) and (c) in Fig. 6.2 can be merged into one unique kernel.

This implies that the new kernel must perform the copy operation of the next wall, which we will refer to as future wall. Appropriately, three-wall introduces a third wall buffer, called future wall. The algorithm sketch of three-wall is presented in Fig. 6.3.

During execution of the wall z, the compute kernel of three-wall copies the next wall (Grid[z + 1]) into future wall. On the next iteration, the future wall becomes the current wall, the current wall turns back to the past wall, and the past wall switches up to the future wall and the process starts again. Now, since current wall (formerly future wall) has already been copied in the previous iteration (or in prolog), there is no more need to ensure ordering of cell updates of the wall of interest, thus the memory barrier is no longer required for three-wall. The total barrier count is therefore reduced to L z , with the same quantity of data movement. Table 6.1 depicts the identical quantity of data movement between two-wall and three-wall algorithm, in either directand indirect-addressing, but with a different number of barrier count.

Table 6.1: Estimation of data requirement for updating one node with both two-wall and three-wall algorithm, and then the barrier count of each algorithm.

Direct addressing Indirect addressing # Type

Comment # Type Comment q PDF Copy: loads q PDF Copy: loads q PDF Copy: stores q PDF Copy: stores q PDF Collision: loads q PDF Collision: loads q PDF Collision: stores q -1 IDX Collision: lookup q PDF Collision: stores Σ = 4q PDFs Σ = 4q PDFs + (q -1) IDXs

Barrier count Two-wall

# barriers = 2L z Three-wall # barriers = L z Moreover, as the copy of future wall is completely independent from the computation (see Fig. 6.3), three-wall exposes the possibility of using streaming load/store instructions to perform cache-bypassing memory operation, avoid polluting useful data in different cache levels of computation. It is also possible to perform asynchronous copy of the future wall on DMA-enabled architectures, so that CPU cores only need to focus on the computation task. Note that the extra memory allocation for temporary wall buffers is considered negligible compared to the main lattice buffer. For instance, threewall buffers represents 2.3% of additional memory on a 128 3 domain ((3 × 128 2 )/128 3 ), or 1.1% on a 256 3 domain. Furthermore, on a non-cubic or complex-geometry domain, the developer can make a choice of loop direction (either in L z , L y or L x ) that satisfies a specific trade-off between performance and memory consumption. Our algorithm can also be extended to more than three walls to fill the performance gap due to the barrier cost, in exchange of the additional memory footprint.

Implementations

We present in this section two short versions of our D3Q19 LBM implementation, the first one within the shared-memory OpenMP paradigm and the second in OpenCL on heterogeneous platforms including accelerators. written (not read) in the current iteration. This removes the read-after-write dependency of two-wall, and exposes possibility of using streaming stores for further optimization. The black and the blue cells are then independently read (b), collided and stored in-place (c).

OpenMP

Both two-wall and three-wall possess a for-loop in one dimension. In this work, we choose to pin the loop on the z-direction and to allocate the grid buffer as float

Grid[L z ][L y ][L x ][Q].
For the sake of simplicity, computation within each L z wall is based on #pragma omp parallel for statements which embrace two inner-loops in L y and L x direction. Barrier in two-wall is implicitly done by two separate OpenMP pragma which involve fork/join from the master thread at each time. Further optimization can use #pragma omp barrier between persistent threads which can even embrace the z-loop.

Wall buffers are allocated as a circular buffer float Walls [START_REF] Ho | Improving 3D lattice Boltzmann method stencil with asynchronous transfers on manycore processors[END_REF]

[L y ][L x ][Q].
On each z wall, threads (or rather the master thread) compute(s) indexes of the past wall, current wall and future wall (if three-wall) as i p , i c , and i f respectively by the modulo operation based on z:

i p = z mod 3 i c = (z + 1) mod 3 i f = (z + 2) mod 3 (6.1)
Here, instead of assigning (z -1) mod 3 to i p , z mod 3 to i c and (z + 1) mod 3 to i f , which makes natural sense, we choose to begin at z mod 3 until (z +2) mod 3 for i p , i c , and i f respectively. This avoids the negative first-case where z = 0, while i p , i c and i f have distinct values ranging in {0, 1, 2}. Note that on two-wall, this range is {0, 1} since mod 2 is used instead of mod 3. Implementation details of two-wall and three-wall are shown in Fig. 6.4 in pseudo-code and Fortran-like OpenMP syntax. Fig. 6.2 and Fig. 6.3 graphically describe the code data accesses patterns. We note that on threewall, a preliminary action (prologue) needs to be taken. It consists in performing copy of the first future wall (z = (0 + 1) mod 3) before entering the for-loop, so that the first iteration can directly access its current wall. (a) Two-wall.

OpenCL

In the implementation of OpenCL, the two #pragma omp parallel statements of the OpenMP code of two-wall is replaced by two 2D OpenCL kernels: one for wall-copy and one for wall-compute. Regarding three-wall, similarly to OpenMP, only one OpenCL kernel is needed inside each z-iteration which includes (a) copy of future wall, (b) collision and (c) in-place update. Implementation of these two kernels in OpenCL is mostly identical to ones of the OpenMP code in the previous section, with the loops in L y and L x broken into two-dimensional L y × L x global work-items. Fig. 6.5 illustrates the OpenCL implementation of two-wall and three-wall, in which we depict the code from the host side (the for-loop in L z direction) which enqueues OpenCL kernels to the device.

Regarding the synchronization between kernels, we employ the cl event management provided by the OpenCL API [8]. Kernels enqueued to the device can be set to a specific execution order by their event-dependency. This allows us to explicitly enforce the FIFO order independently from any specific OpenCL driver implementation, which for a more efficient utilization of the hardware, enables the out-of-order execution by default.

1: for z in 0 .. Lz-1 do 2:

Set kernel argument CopyWall.z = z;

3:

Enqueue FIFO kernel CopyWall to device;

4: 5:

Set kernel argument Compute2Wall.z = z;

6:

Enqueue FIFO kernel Compute2Wall to device;

7: end for 8:

(a) Two-wall.

7:

Enqueue FIFO kernel Compute3Wall to device;

8: end for (b) Three-wall.

Results and discussions

In this work, we implement the lid-driven cavity use-case based on the OPAL LBM kernel presented in [START_REF] Obrecht | Performance Evaluation of an OpenCL Implementation of the Lattice Boltzmann Method on the Intel Xeon phi[END_REF], originally written in OpenCL-C. Multiple-relaxation-time (MRT) collision from [START_REF] Humières | Multiple-relaxation-time lattice Boltzmann models in three dimensions[END_REF] is fused with the pull propagation and the half-way bounce-back boundary condition. The collision kernel performs about 350 floating-point arithmetic operations per lattice update. Four propagation algorithms are implemented for comparison: (

Two-lattice, (2) AA-pattern, (3) Two-wall and (4) Three-wall. Two-lattice is considered as the state-of-the-art algorithm. AA-pattern is chosen as an efficient algorithm which requires only one instance of the lattice buffer. Then two-wall and three-wall are implemented for comparison to the two former algorithms.

For comparison between the four propagation algorithms, we will use the customary performance metrics in LBM, i.e. million lattice-node updates per second (MLUPS) and the memory occupancy efficiency in terms of LUPS per byte, referred to as the perf-mem ratio. The main reason of comparing the perf-mem ratio is to understand the performance efficiency of algorithms on a given amount of memory, in a similar way as GFLOPS per watt account for algorithmic energy efficiency. Results of each implementation are checked for correctness against the reference OPAL code on a GPU.

OpenMP

The OpenMP code runs on two shared memory platforms:

• 24-core 12x2 NUMA CPU Intel Xeon Haswell E5-2680v3 at 2.5 GHz For both two-lattice and AA-pattern, the three nested for-loops in z, y, x are organized and collapsed as follow:

#pragma omp parallel for private(x,y,z) collapse( 2)

for(z = 0; z < Lz; z++)

for(y = 0; y < Ly; y++)

for(x = 0; x < Lx; x++)

This configuration leaves the possibility to the icc compiler of vectorizing the x-loop and gives the best performance in our tests. Regarding the computation code, the LBM kernel of OPAL is rewritten from OpenCL-C to the multi-dimensional array C99 standard. Given the similarity between OpenCL-C and C99, the porting process did not raise much difficulty, which consists mostly in copying the propagation and collision code blocks inside the three nested zyx-loops. Performance in MLUPS of the four propagation algorithms on CPU is shown in Fig. 6.6a.

On CPU, AA-pattern delivers similar performance to the one ot two-lattice algorithm, that shows its benefit over two-lattice in terms of memory requirement. Performance of two-wall and three-wall are between 20 to 40% lower than one of two-lattice and AA-pattern, which conform to 4q of data movement per lattice update. As expected, three-wall is faster than two-wall on CPU, but by only 2-9%, certainly thanks to the 1 The kona01 node from the PlaFRIM system which is configured with flat memory-mode and quadrant cluster-mode [START_REF] Codreanu | Best Practice Guide -Knights Landing[END_REF]. On both processors, the compiler version of Intel 2017 update2-knl is used. Compilation flags are given in the Tab. 6.2, by following the recommendations from [START_REF] Codreanu | Best Practice Guide -Knights Landing[END_REF]. Default OpenMP affinity environment variables like number of threads and thread locality are also given the Tab. 6.2. We observe that using the same number of OpenMP threads as the number of physical cores yields better performance, as it avoids additional cost related to thread-switching and cache-pollution. Thread-binding is also carefully chosen to optimize cache locality within NUMA cores on CPU or the shared common L2 cache of a tile (2 cores) on KNL.

efficient OpenMP runtime on the CPU cache system. However, we observe that twowall turns to be faster than three-wall from problem sizes larger than 384 3 . This is explained by the data size of three-wall buffers that exceeds the L3 cache size of CPU (30 MB) from 384 3 domain: 3 × (384 2 × 19 × 4) ≈ 32 MB. Performance of both two-wall and three-wall follows an oscillating line that fits the L2 and L3 cache, at 128 3 and 384 3 respectively. Performance efficiency of algorithms in terms of LUPS-per-byte is shown in Fig. 6.6c, whose values are normalized to the two-lattice algorithm. The AApattern, which consumes twice less memory than two-lattice while delivering the same performance, is two times more efficiency in terms of memory occupation. Besides, threewall and two-wall, although they do not reach the same efficiency as the AA-pattern, are always better than two-lattice.

Moreover, we observe interesting results on KNL. As seen from Fig. 6.6b, two-lattice considerably outperforms three other algorithms (300 MLUPS versus 225 MLUPS of AA-pattern). This performance gap is justified by the fact that icc has succeeded at generating non-temporal streaming stores [START_REF] Ronald W Green | Memory movement and initialization: Optimization and control[END_REF] (movntdq) for the two-lattice propagation. Instruction generation is confirmed by reading the executable binary of two-lattice with the objdump command. This instruction is introduced on Xeon Phi generations for bypassing cache levels when writing data to the memory (the second lattice buffer of two-lattice). This avoids expensive read-for-ownership (RFO) operations and saves memory bandwidth on the Xeon Phi architecture, which is known to be more sensitive to cache-pollution because of the write-allocate policy. Interestingly, however, performance of three-wall this time competes with AA-pattern and outperforms two-wall by 15-20% for problem sizes larger than 384 3 (≈ 32 MB of last-level L2 cache). This result is promising because it shows that one can use three-wall as an alternative to AA-pattern to implement complex LBM boundary conditions on the KNL architecture, while obtaining the same performance and perf-mem ratio (see Fig. 6.6d). Performance of two-wall on KNL is not stable on small and medium problem sizes between 64 3 and 384 3 (zigzag line). There are two possible reasons to this phenomenon: (1) important number of inter-thread barriers (2L z ) and ( 2) dependence of memory alignment for twowall buffers with respect to problem size. These two issues are likely to be amplified by the latency from the directory-based L2 coherence protocol on KNL and would require further investigation to explain the phenomenon.

As a result, we prove that three-wall and two-wall can be implemented on a CPU system with higher memory-efficiency than two-lattice. Moreover, using three-wall on CPU does not offer significant gain over two-wall, because the cache coherence protocol can still manage to deliver satisfactory latency and bandwidth in a low-concurrency context (24 cores). However, it is promising to employ three-wall, which outperforms clearly two-wall, on a high-core count architecture like KNL (versus traditional CPU) as the post-collision in-place update of PDFs for each of these three algorithms. OPAL is implemented in a way that one can easily switch the memory layout depending on the target architecture, typically between AoS (favorable for CPU) and SoA (favorable for GPU). OpenCL libraries used in our tests are summarized in Tab. 6.3, in which global and local work-items configuration of kernels are also given. Similar to OpenMP, in OpenCL we can have either 2D or 3D kernels corresponding to the four propagation algorithms (2D for three-wall/two-wall and 3D for two-lattice and AA-pattern). We study performances using 32 as the problem size increment. It should be mentioned that setting OpenCL work-group size to 32 × 1 × 1 already gives satisfactory performance on target platforms. Incrementing work-group size to for instance 64 × 1 × 1 does not yield important gains, while constraints the problem size to be multiple of 64. Surprisingly in Fig. 6.7a, despite using the same CPU, the AA-pattern OpenCL code delivers significantly higher performance than the OpenMP version (290 MLUPS versus 210 MLUPS). This represents up to 2.75 times higher in perf-mem efficiency with respect to the normalized two-lattice ratio (see Fig. 6.7c) (versus ×2 in OpenMP). This leads us to think that the Intel OpenCL SDK, by some advanced SIMD optimizations, performs much better than the icc compiler on OpenMP pragma's of the AA-pattern propagation.

The Intel OpenCL runtime might have detected that AA-pattern always writes data back to the same location of earlier reads and generates the corresponding optimized code. More fine-tuned compilation options for icc could be the missing point in our tests. However, determining the corresponding flags or using advanced OpenMP hints and extensions is not the main goal of this work. For the sake of portability and proofof-correctness of algorithms, we do not include any platform-specific code annotation and only use default compilation flags as recommended in [START_REF] Codreanu | Best Practice Guide -Knights Landing[END_REF]. Moreover, we see that implementing two-/three-wall in OpenCL on CPU is less favorable (125 MLUPS) than directly using OpenMP (180 MLUPS), with respect to the performance and perf-mem ratio. The oscillating line, seen in the OpenMP code, is replaced by a normal lower performance line, probably due to the OpenCL runtime overhead and cache flush after each kernel completion.

Furthermore, the latest Pascal GPU, with its HBM2 memory, delivers gratifying performance with more than 3000 MLUPS on two-lattice and 2700 MLUPS on AA-pattern, execution failure of the AA-pattern algorithm for problem sizes larger than 352 3 on the Tesla P100 GPU, which can be a drawback in applying AA-pattern on GPU-based platforms.

Conclusions

In this chapter, we introduce two-wall and three-wall, two novel algorithms of the LBM one-lattice algorithm class which use a single lattice instance. Two-wall and three-wall are easy to implement and can be used on complex boundary conditions.

Our algorithms do not deliver the highest performance on CPU and GPU compared to other algorithms, but offer better memory efficiency (LUPS per byte) than the stateof-the-art two-lattice algorithm. Especially on a high-core-count shared-memory architecture like KNL, three-wall gives competitive performance and memory efficiency compared to AA-pattern, one of the most efficient propagation algorithms in literature.

Employing three-wall on future many-core architectures which embed more local memory and private cache as well as the non-temporal streaming store instruction could lead to further performance improvements. Furthermore, hardware architecture in near future is likely to favor large on-chip coherent memory associated to many-core designs.

For such evolutions, our approach will become mandatory to increase performance with a straightforward OpenMP implementation. Using a single lattice instance also allows to increase spatial resolution and thus to obtain more reliable solutions of the physical problem under consideration.

Introduction

In this chapter, we propose the design of an MPI Message-Passing library [9] for the intra communication on many-core processors, using the vendor support library (MP-PAIPC [START_REF] Dupont De Dinechin | A Distributed Run-Time Environment for the Kalray MPPA R -256 Integrated Manycore Processor[END_REF]) as the transfer-fabric to build MPI protocols from scratch, while porting any of existing MPI implementations such as MPICH or OpenMPI would not be possible due to limited on-chip memory of most recent many-core processors.

Based on studied MPPA hardware specifications presented in [START_REF] Dupont De Dinechin | A clustered manycore processor architecture for embedded and accelerated applications[END_REF][START_REF] Dupont De Dinechin | A Distributed Run-Time Environment for the Kalray MPPA R -256 Integrated Manycore Processor[END_REF], this chapter does a brief hardware summary and focuses on an MPI design over (but not limited to) the MPPA architecture, with detailed implementation algorithms and formulated models following vendor-hardware characteristics (K, h) and different optimizing approaches (Lazy, Eager). These studies are generic enough to be compared/ported to other verysimilar architectures, such as Tilera [START_REF]Tile processor architecture overview for the Tilepro series[END_REF], STHORM [START_REF] Mottin | The STHORM Platform[END_REF] or Neo chip [START_REF] Hemsoth | The Tiny Chip That Could Disrupt Exascale Computing[END_REF], on which doing/optimizing MPI communication over Network-on-chip is still a challenging or never-posed question.

The remainder of this chapter is organized as follows: Related MPI-oriented works on other many-core platforms are compared in section 7. Using the MPPA-MPI library, the HPL benchmark [START_REF] Petitet | HPL -A Portable Implementation of the High-Performance Linpack Benchmark for Distributed-Memory Computers[END_REF][START_REF] Jack J Dongarra | The LINPACK benchmark: past, present and future[END_REF] was ported on MPPA with the support of the standard BLAS-Netlib [START_REF] Chuck L Lawson | Basic linear algebra subprograms for Fortran usage[END_REF][START_REF] Jack J Dongarra | A set of level 3 basic linear algebra subprograms[END_REF] (mono-threaded) and OpenBLAS [START_REF] Xianyi | OpenBLAS, version 0.2. 8[END_REF],

an OpenMP optimized implementation. These benchmark results are summarized in section 7.7 and conclusions are given in section 7.8.

Related works

Our design is similar to the co-processor-only MPI model on the Intel Xeon-Phi platform [START_REF] Jeffers | Intel Xeon Phi Coprocessor High Performance Programming[END_REF], with support of OpenMP for hybrid multi-threaded programming. Besides MPI ranks running on CC, we introduce an MPI I/O rank running on an I/O subsystem of the chip as bridge to communicate with the host through the PCI-e interface, while there are no direct communication link between the host and the compute clusters on MPPA. Along the way, some collective MPI functions were also implemented (MPI Comm split, MPI Bcast, MPI Reduce, MPI Allreduce and MPI Barrier).

Our message-trigger handling mechanism using the RM core was inspired by the similar work of Prylli and Tourancheau [START_REF] Prylli | BIP: a new protocol designed for high performance networking on myrinet[END_REF] [77] implementing the BIP protocol for an optimized MPI implementation over the Myrinet network, taking advantage of its dedicated hardware, an extra core like the MPPA RM core.

Today, there exist other multi-/many-core processors similar to the MPPA. Some of them has an MPI implementation, others do not. This section reviews MPI-oriented libraries on other many-core architectures and their performance related to our work.

Raw processor [START_REF] Bedford | The Raw microprocessor: A computational fabric for software circuits and general-purpose programs[END_REF], designed by the Computer Science Laboratory at MIT, combines 16 identical compute units, called tile. The 16 tiles are connected by one static NoC and two dynamic NoC. The static network is used for predefined memory access pattern at compile time, the dynamic ones are used for communication scheme at runtime. Psota

and James [START_REF] Ryan | rMPI: An MPI-compliant message passing library for tiled architectures[END_REF] propose rMPI, the first MPI library over the Raw achitecture by inheriting some design aspects from MPICH and LAM/MPI also other specific implementation belonging to the Raw hardware. The highest throughput obtained on the ping-pong test of the Raw processor is about 150 MB/s with buffer size of 3.2 MB (100K words) [START_REF] Psota | rMPI: message passing on multicore processors with on-chip interconnect[END_REF].

Tilera processors [START_REF]Tile processor architecture overview for the Tilepro series[END_REF] (KNC), proposes three MPI programming models [START_REF] Jeffers | Intel Xeon Phi Coprocessor High Performance Programming[END_REF] which are ( 1) offload (host-only),

(2) co-processor-only and ( 3) symmetric (both host and co-processor). The MPI communication in the intra-MPPA context corresponds to the co-processor-only intra-MIC case. Potluri et al. [START_REF] Potluri | MVAPICH2-MIC: A High Performance MPI Library for Xeon Phi Clusters with InfiniBand[END_REF] studied the communication throughput of the MVAPICH2 library on KNC and their results show that a MIC-optimized MVAPICH2 library can delivers more than 9 GB/s of uni-direction throughput for messages up to 1 MB.

MPPA-MPI design

In the MPPA context, each CC is referred to as an MPI rank. Thus, the MPPA-256 MPPA processor (i.e. 17 ranks). The MPPA-MPI architecture on each rank is then divided in two layers:

processor

MPI-inter-process Control (MPIC)

Each MPI transaction begins by exchanging control messages at the MPIC layer between MPI ranks. Control messages are used for:

• information exchange about MPI transaction type (send/receive, communicator split, etc.).

• synchronization point in case of rendez-vous protocol.

We implement an RQueue-based active message server [START_REF] Von Eicken | Active messages: a mechanism for integrated communication and computation[END_REF] In an MPI send/receive, The Tx rank posts a Request-To-Send (send post t) to the Rx rank; idem, the Rx rank sends back a Clear-To-Send (recv post t) containing its allocated dnoc tag, to which the Tx rank will send data. Beforehand, this dnoc tag needs to be configured and linked to the receive buffer to enable remote writing.

MPI-inter-process Data-Transfer (MPIDT)

MPIDT is a light-weight wrapper of MPPAIPC Portal primitives. Once the Tx rank has got a matching control message, it configures a data transfer using received information (e.g. dnoc tag). Data can then be sent in either blocking or non-blocking mode dependent on the calling MPI function, using appropriate Portal primitives (see Tab. 7.1 for detailed function mapping).

MPPA-MPI implementation

As mentioned above, on-flight control messages carry essential information depending on their purpose. We present now their usage as well as algorithms of the two communication scenarios in our work:

(1) synchronous blocking send (MPI Send) and

(2) asynchronous non-blocking send (MPI Isend)

MPI Send -MPI Recv

Most well-known and optimized MPI libraries contain many (combined) techniques to perform the MPI Send call. In the first time, we chose to implement this function with rendez-vous blocking behavior, in order to avoid extra buffer space and minimize memory usage. This choice certainly adds more synchronization cost but does not change the functionality of the send/receive transaction. Some optimization approaches will be presented in the coming sections. Algorithms 2 and 3 summarize the implementation of MPI Send and MPI Recv.

Algorithm 2 MPI Recv(buf, count, datatype, source, tag, comm, status) The implementation of MPI Irecv is the same as the one of MPI Recv, except that the function returns right after having posted the receive to the sender, and the completion of reading (mppa aio wait) is done in MPI Wait.

Algorithm 3 MPI Send(buf, count, datatype, dest, tag, comm) recv post ← find recv post(count, datatype, tag, ...); 7: until recv post = NULL 8: /* send data (MPIDT layer), using mppa pwrite(s) */ 9: send data(buf, count, datatype, recv post.dnoc tag); 10: return MPI SUCCESS;

MPI Isend -MPI Recv

The implementation of MPI Isend uses non-blocking Portal primitives on both PE and RM on the Tx side. When the Tx rank (PE0) reaches MPI Isend in its execution without having received any matching recv post, it creates a non-started pending isend request containing related information (buffer pointer, dest, count, tag etc.) and returns. On arrival of the matching recv post, the RM core (callback handler) reads the previous pending isend request and triggers a non-blocking data send (to the recv post.dnoc tag of the Rx rank). The request is then set to started state to be distinguished from other non-started requests.

On the other hand, when the recv post arrives before MPI Isend, the RM core saves it into the internal buffer. The PE core executing MPI Isend later reads this recv post, performs a non-blocking send and marks the pending isend request as started. This propriety ensures that the transfer is performed only once for each transaction, either by the PE core (in MPI Isend) or by the RM core (in callback handler). At the end, started requests will be finished and cleaned by MPI Wait. Algorithms 4 and 5 present in more details the implementation of MPI Isend and of the callback handler.

Optimization

Eager send optimization

The idea is to pack any MPI message which can fit into a 120-byte space, as a controlmessage and send it directly to the Rx active server. In reality, the maximum data payload is about 96 bytes (24 bytes is used for control header). An eager buffer needs to be allocated on each MPI rank and can be defined by the EAGER BUFFER LENGTH macro in main.c. This approach is synchronization-free when the MPI Send call can Algorithm 4 MPI Isend(buf, count, datatype, dest, tag, comm, request) save recv post(recv post); 10: end if 11: return ; return before a matching receive is posted (non-local). It also leads to an improvement of about 6 to 10% in performance for the HPL benchmark on MPPA (see Section 7.7).

For longer messages, using several eager sends introduces segmentation and reassembly costs. A test case is implemented with messages split into eager packets in order to determine the best communication trade-off in Fig. 7.4. Such segmentation however consumes buffering memory and therefore increases RM workload.

Lazy send optimization

Lazy send consists in copying medium-size message into a local buffer and returns. The RM is then responsible for sending it to the destination. Unlike eager buffer on the Rx side, lazy buffer is allocated on the Tx side and can be tuned via some macros (LAZY THRESHOLD, LAZY BUFFER LENGTH).

This approach must be used with care because bad communication scheduling may lead to buffer wasting and lazy messages remaining for too long. Inversely, a dense communication scheme should neither be set to lazy mode in order to be able to send data directly rather than spending time doing memcpy in local memory.

DMA thread usage

MPI Isend uses Portal non-blocking primitive to configure a Tx DMA thread for data sending. The DMA engine implements a fetch instruction that loads the next cache line while pushing the current line into the NoC. This fetch is nowadays not available on PE cores, meaning that outbound throughput using PE is 4 times lower than using DMA engine (1 B/cycle vs. 4 B/cycle). Thus, tuning to use non-blocking DMA on MPI Send for messages of size greater or equal to DMA THRESHOLD will maximize the transfer performance.

MPPA-MPI Throughput modeling

The MPPA-256 Network-on-chip [START_REF] Inc | MPPA-256 Cluster and I/O Subsystem Architecture[END_REF] is designed so that any path linking two CCs always contains less than eight hops (including two local hops -one at sender and one at receiver). The average switching time on a NoC router is 7 cycles, then it takes the packet at most 8 cycles to reach the next hop. In the worst case, the link distance (time a packet spends on NoC to reach its destination) is 112 cycles (7×8+8×(8-1)). However, the necessary time to send a buffer (transmission time -t) is about O(N ) cycles [START_REF]MPPAIPC Performance[END_REF] (where N the buffer size in bytes), which is much longer than the link distance [START_REF] Kumar | A network on chip architecture and design methodology[END_REF].

As a result, we describe the transmission time t as a function of the buffer size N , a constant transfer ratio K and a default overhead h (aka. the cost of sending an empty buffer). This default overhead presents the initial cost of MPI implementation management (ID mapping, metadata setup, synchronization, error checking ...) and/or configuring the peripherals (cache, DMA) to prepare for data sending. This cost is paid on each MPI call and is independant to the subsequent data-sending process (which is presented by a data-transfer factor K). The ping-pong round-trip time (RTT) is approximately the sum of the transmission time on both sides, as the propagation time is negligible.

T ransmissionT ime : t = K × N + h (cycles) (7.1) RT T ≈ 2 × t = 2 × (K × N + h) (cycles) (7.2)
T hroughput : T = 2 × N RT T ≈ N K × N + h = 1 K + h N (bytes/cycle) (7.3) lim N →∞ T ≈ lim N →∞ 1 K + h N = 1 K (bytes/cycle) = 400 × K -1 (M B/s) (at f requency 400 M Hz) (7.4)
The constant K is a value specific to each send function with its own underlying transport primitive. For example, the MPI Isend which uses the DMA engine with peak throughput of 4 B/cycle, would have its transfer ratio K of about 0.25. The MPI Send, with default peak throughput of 1 B/cycle (no DMA engine), should obtain a transfer ratio K around 1.

Results and Discussion

Using the MPPA Developer platform [START_REF] Inc | Kalray Platforms and Boards[END_REF] with the first-generation MPPA processor, we set up ping-pong tests between:

( In each case, the same MPI send function is used on both sides (MPI Send or MPI Isend).

At the first time, all tests are run without any optimization in order to calibrate the proper throughput of each context (Fig. 7.3). At the second time, we enable all optimization on the MPI Send test and compare our optimization approaches in terms of latency, throughput and messages sent per second (Fig. 7.4).

Each ping-pong is repeated 50 times. We assume that there is no waiting time inside the MPI send function, since all ranks start at the same time and run at the same clock speed. Hence, the duration of the MPI send function can be considered as the transmission time. Depending on the send context, the measured transmission time is fitted into a linear correlation K × N + h presented in Tab. 7.2. The standard deviation from all obtained results is always less than 0.2%.

From

CC-IOS communication

Contrary to the symmetric communication performance between CCs, the transmission rate on I/O subsystem relies on the DDR bandwidth, which is much lower than the on-chip memory on CCs. We observe higher K values and much more considerable overhead h on the IOS 128, showing that the communication link from IOS to CCs might be the bottleneck on the MPPA. It is then difficult ignoring h N in this case. Keeping on our throughput estimation by 400 × K + h N on Fig. 7.3b, where the performance gap between the CC 15 and the IOS 128 is also illustrated.

Optimization comparison

We focus now on finding, on a given message size, the best send method among the four (Normal, Eager, Lazy and DMA) to use on MPI Send, in order to obtain lowest latency (round-trip-time) and/or highest ping-pong throughput, by enabling all optimizations and re-running our experiments between CCs. We also evaluate the number of messages sent per second in each approach by dividing the clock frequency (400 MHz) by the duration of the MPI Send call (in cycles). As the message will now be eagerly sent or lazily buffered and MPI Send returns right afterward, this duration on Eager(-splitting) or Lazy could no longer be evaluated as the transmission time in the Tab. 7.2, but respectively by :

E × N 96 + 1 (cycles, E ≈ 3800) (7.5)
O memcpy (N ) = 1.28 × N + 5300 (cycles) (

where E is the constant necessary cost to send 1 eager-split and O memcpy is a linear function of memcpy cost. Note that in the Lazy approach, the message is sent in background by the RM.

Hence, Eager and Lazy methods provide lower latency and higher message rate on short buffers, since they were designed to get rid of two-sided synchronization and the buffer size is still small enough not to be outperformed by the DMA high-throughput capacity. Fig. 7.4a shows that the ping-pong latency from 1 to 256 bytes using eager-splitting is reduced by half compared to DMA or Normal. Otherwise, this latency increases radically as soon as its transmission time, despite being smaller at the beginning, getting repeated as many times as split segmentation N 96 + 1 . On the other hand, using DMA on large buffers optimizes bandwidth utilization compared to Normal (using PE) or Lazy (using RM) methods. (Fig. 7.4b). Fig. 7.5a shows the HPL result on MPPA-256 using BLAS and OpenBLAS. Note that 70 GFLOPS of annouced theoretical performance is for all the 256 cores, while the best benchmark score (R max = 1.2 GFLOPS) was achieved using only one core per CC (i.e. 16 cores in total) and MPI eager send. Also, we have seen no performance change by enabling MPI lazy optimization. This can be explained by well-scheduled HPL overlapping [START_REF] Jack J Dongarra | The LINPACK benchmark: past, present and future[END_REF] in which, either MPI processes arrive to the communication step at the same time, or all heavy sends are done asynchronously by MPI Isend, while lazy sending only shows its advantage in bad-scheduled MPI Send. Furthermore, multithreading on MPI compute ranks (OpenMP on CCs) did never give better HPL result, because of the small working set and the OpenMP overhead.

Limitations and conclusions

In this chapter, we have introduced the design and performance issues of an MPI implementation on the Kalray MPPA-256. The MPPA-MPI library provides 1.2 GB/s of throughput for any inter compute-cluster point-to-point communication and this performance depends on the underlying MPPAIPC library. Optimization ideas such as eager send and lazy message are proposed, implemented and compared to determine the best approach based on message size. A synthetic model is also presented for each approach to evaluate their communication latency and throughput. We also learn that supporting MPI programming model is not an easy task on recent many-core processors, including MPPA, since MPI has become a large API with high-level abstractions and many-core hardware is taking more diversity and complexity. Thus, optimizing an MPI implementation on each of these platforms is even more not trivial.

Despite the fact that the HPL benchmark was also successfully ported on MPPA as a validation test of our MPI library, we think that mapping an MPI rank on a single compute cluster of MPPA is not a long-term solution. This topology mapping can provide some ease of programming in the short-term, but will encounter a scalability issue in the future when we will want to connect hundreds or thousands of MPPAs together. The granularity will be too fine and such mapping will face a significant General Matrix Multiplication (GEMM) on Many-core

Processors

The best performance improvement is the transition from the non-working state to the working state.

-John Ousterhout.

Introduction

The largest difference of a DMA-based many-core architecture from other conventional platforms is the lack of a data cache system. On a DMA-based architecture, all data must be accessed and copied from the off-chip global memory to the on-chip local memory. This not only improves data locality and latency, but also allows computationcommunication overlap by asynchronous transfers.

Our objective in the present chapter is to implement a fast GEMM on DMA-based many-core processors, typically the MPPA2-256. Consequently, we are focusing on asynchronous algorithms to reduce memory latency and an MPPA assembly micro-kernel to obtain the highest performance. Therefore, we do not plan to support a full range of GEMM parameters from the BLAS API in this chapter.

8.2 GEMM in POSIX-C

Algorithm

In this work, we use blocksize naming convention from Matsumoto et al. [START_REF] Matsumoto | Performance tuning of matrix multiplication in OpenCL on different GPUs and CPUs[END_REF], which implement an auto-tuning blocked GEMM implementation on OpenCL devices, based on configurable runtime parameters. The matrix C is divided in blocks of size M wg × N wg (C i ). Blocks of A are partitioned as M wg × K wg (A i ), and blocks of B as K wg × N wg (B i ).

On MPPA2, the two-layer memory configuration matches. A large and slow DDR memory versus a small and fast scratchpad memory on each Compute Cluster (CC). Computing a block C i requires multiplying a block-row of A i to a block-column of B i (see To compute a block C i within each cluster, blocks A i and B i are streamed from the DDR to scratchpad in the double-buffering overlapping scheme. At any time, there is always

Assembly-level GEMM micro-kernel

In the above section, we have described the asynchronous GEMM algorithm between DDR and scratchpad memory. Once blocks A i , B i and C i are copied into scratchpad and zero-padded if necessary, the local computation can proceed. of B p to N wg , instead of respectively m r and n r as defined in BLIS. 3 The MPPA ISA provides a streaming-load instruction which performs prefetching from the scratchpad memory into registers, bypassing the cache L1. Using this prefetch instruction allows another finer-granularity overlapping level: in-core computation and register-scratchpad communication, versus the DDR-scratchpad overlapping described in the previous section. Combining these two overlapping schemes enables parallel and highly efficient GEMM computation on a DMA-based, cache-free architecture like MPPA. Implementation of the corresponding micro-kernel is summarized in Algorithm 6. Regarding the two-and three-dimensional asynchronous copy functions, we implement a set of extended asynchronous primitives of the ones defined in the OpenCL specification.

GEMM in

Load C i [m r ][n r ] → registers C r [m r ][n r ]
They are: (1) a general-strided copy function which supports remote and local strides, respectively in global and local memory (the original one does not support stride on local memory), (2) 2D asynchronous copy and (3) 3D asynchronous copy. These two later primitives are mapped one-to-one over the underlying 2D and 3D ones from the mppa-async library. Similarly, these new OpenCL primitives facilitate the writing of the GEMM-async algorithm in OpenCL.

GEMM in OpenCL POSIX-like (Task-Parallel)

In the above section, we have identified the two following performance limitations of GEMM-async in OpenCL data-parallel, compared to the POSIX-C model:

• Small scratchpad memory per compute-unit. This execution and memory mapping reduces DMA blocksizes, results to I/O bound and performance loss.

• Lack of support for assembly code or pre-compiled objects. A carefully handtuned BLIS micro-kernel could deliver higher efficiency than the built-in OpenCL compiler, especially on the MPPA VLIW architecture.

Results

The three GEMM-async implementations are run on the MPPA2-256 processor oper-

Conclusions

In this chapter, we introduce a set of GEMM implementations over existing programming models, as well as new proposed extensions on the MPPA2 processor. The best performance is obtained on the POSIX-C model with large development efforts (≈ 4000 lines of code). The default OpenCL data-parallel mode, despite reducing the programming cost, encounters the scratchpad memory limitation due to the execution mapping. The new OpenCL POSIX-like extension combines assets of each one of the two above models and provides an elegant and high performance programming mode for the MPPA2 processor.

Within a reduced development cost of about 2500 lines of code, the OpenCL POSIX-like yields a performance three times higher than the OpenCL data-parallel mode and 59% compared to the POSIX-C one, with rooms for improvement by using two DDRs and the largest local memory space.

Nevertheless, those implementations are only for benchmark and study purpose. They are not usable in system-wide software applications or production frameworks which, for the sake of portability, commonly invoke only standard APIs and libraries. For instance, convolution algorithms in deep learning or simulation solvers in physics, biology and chemistry frequently bind their compute kernels on dense and sparse linear algebra operations. Providing an optimized sparse and dense BLAS library has always been a decisive condition to any processor architecture. However, achieving that task is not always manageable.

Chapter 9

Portable and Optimized BLAS Library on Many-core Processors ...I would be panicked if I were in industry. Now I'm forced into an approach that I haven't laid the groundwork for, it requires a lot more software leverage ..., and the microprocessor manufacturers don't control the software business.

-John L. Hennessy, Stanford University.

Introduction

Besides traditional CPU and GPU platforms, recent initiatives and strategies for energyefficient HPC [START_REF] Rajovic | Supercomputing with commodity CPUs: Are mobile SoCs ready for HPC[END_REF] [93] also look at low-power system-on-chip (SoC) processors as the building blocks. These SoC processors integrate a reduced instruction set CPU accelerated with a high number of in-order DSP, VLIW or vector cores operating on scratchpad memories. To strip off power-consuming sources, accelerator cores do not implement outof-order execution, hardware prefetcher, or globally coherent caches. Applications are expected to focus on local memory to overcome the high latency of the main memory and to leverage Direct Memory Access (DMA) engines to overlap communication and computation. Integrating DMA capabilities into a BLAS library represents a design and development challenge. This combines all the complexity of the BLAS library development, known to be elaborate, expertise-requiring and time-consuming, with the asynchronous aspects of DMA transfers, the multiplicity of memory spaces and the high execution concurrency. Failure to design a portable DMA-based BLAS will result in repeating the same time-consuming process on each new architecture. Our key contributions are as follows:

1. Asynchronous implementation of the level-3 BLAS in the BLIS framework as a flexible and fully-compatible module to leverage DMA capabilities on embedded platforms with minimal and deterministic memory footprint. 1 2. Definition of a generic and user-friendly DMA back-end interface of six functions, inspired by the classic asynchronous one-sided put/get RDMA operations for interconnection networks, firstly proposed in the SHMEM paradigm [START_REF] Pakin | Fast Messages: Efficient, portable communication for workstation clusters and MPPs[END_REF]. The purpose of the back-end interface is to unify low-level communication libraries from different platforms and to provide a coherent mapping of BLIS primitives on any DMA-based architecture. Contributed code is released under the same BSD license as BLIS. 2 3. Validation of the principles on the Kalray MPPA2-256 many-core processor, with mapping of the back-end DMA interface over the Kalray DMA library as an evidence of portability. The implementation delivers 75% of peak performance on a single-core execution with a memory footprint of 480 KB. This reference implementation passed out of the box the BLIS test suite 3 on CPU as another evidence of portability and correctness. 1 We refer to memory footprint as the scratchpad memory space needed for DMA and packing buffers. This footprint depends on the block size parameters, the overlapping scheme of each matrix as well as the parallelization of different loops. These aspects will be discussed later in the chapter.

2 https://github.com/hominhquan/blis/tree/rdma 3 https://travis-ci.org/hominhquan/blis

The remainder of this chapter is structured as follows. Section 9.2 presents some related work about using DMA transfers to optimize BLAS on different architectures. Section 9.3 introduces a design of a portable DMA back-end for the BLIS framework, and motivates the key technical decisions which respect to embedded memory constraints.

Section 9.4 shows obtained performance on the Kalray MPPA2-256 processor with detailed observation and explanation of the results. We conclude and outline future work in Section 9.5.

Related work

Exploiting DMA transfers has been the primary option when developing applications on Digital Signal Processors (DSP), Field Programmable Gate Array (FPGA) and other accelerator architectures fitted with local memories. Early works on developing DLA library were presented for FPGA [START_REF] Dou | 64-bit floating-point FPGA matrix multiplication[END_REF], DSP [START_REF] Francisco D Igual | Unleashing the high-performance and low-power of multi-core DSPs for generalpurpose HPC[END_REF] and the CELL processor On the Kalray MPPA2-256 processor, there are 32 MB of Local Memory distributed across sixteen compute-clusters [START_REF] Dupont De Dinechin | A clustered manycore processor architecture for embedded and accelerated applications[END_REF]. In the remainder of this chapter, we will refer to these local memories in general as scratchpad memory (SMEM). A scratchpad memory is defined as a user-allocatable and DMA-accessible memory space, close to the cores in the memory hierarchy and thus accessible at lower latency and higher bandwidth than global memory.

On a specific architecture, the scratchpad memory may reside in L1/L2/L3 locked cache partitions (if accessible to DMA) or in a physically distinct memory region. It should be mentioned that in the case of the TI DSP, different cache levels can be configured as multiple levels of SMEM. However, at the present work, we consider only one level of SMEM. The main reason to that is that we do not know of any other hardware in the market with memory characteristics similar to the TI DSP. Another reason is that supporting multi-level SMEMs would add considerable complexity in the BLIS-RDMA control trees [START_REF] Field | BLIS: A framework for rapidly instantiating BLAS functionality[END_REF], which we strive to keep as simple as possible in the first instance, in order to be used as a reference point for any further platform-specific customization. following approaches proposed in [START_REF] Goto | High-performance implementation of the level-3 BLAS[END_REF]. Algorithm of each of these two routines is implemented as a control tree following the block-panel-based loop organization described in [START_REF] Field | Implementing high-performance complex matrix multiplication via the 3m and 4m methods[END_REF]. In order to develop a comprehensive DMA support with as few modification as possible, we re-apply the block-panel-based algorithm by instantiating two new DMA control trees for GEMM and TRSM respectively.

Pseudo-code to integrate DMA transfers into the control tree is summarized in Fig. 9.2, where code modifications for the DMA extension are highlighted using colored text.

Sub-partitions of A, B and C are copied by DMA from the main memory into SMEM DMA buffers, then packed in situ in an appropriate contiguous memory layout [START_REF] Goto | Anatomy of high-performance matrix multiplication[END_REF] ( A and B). Packing facilities are used unchanged. Computation is then performed on these SMEM packed buffers A, B and C dma . C dma is then copied back to the main memory as described in Fig. 9.2. In GEMM, as matrices A and B are read-only, double-buffering is used for copying blocks/panels of A and B into SMEM. Matrix C, which is used for both input and output, is traversed under a triple-buffering scheme: one buffer in computation, one buffer in prefetch (get) and one buffer in writing back (put).

Memory management

Control trees of BLIS-RDMA require SMEM dynamic allocation to perform on-demand asynchronous copies between the main memory and the SMEM. Conveniently, BLIS comes with a memory broker managing a pool of memory blocks, mainly used for the packing process. However, those blocks are often inflated to match some given algorithms or a pre-defined alignment that are needed by the underlying macro-and micro-kernel.

As the DMA copy occurs before the packing and does not interface directly with either macro-or micro-kernel (see Fig. 9.1b), the DMA buffer can be allocated as a dense and contiguous block in order to reduce the SMEM footprint. We thus added the ability of using a scratchpad allocator into BLIS-RDMA so that the control trees can allocate exactly the needed amount of SMEM for DMA transfers. This allocator is mapped on malloc by default and can be re-defined by user. Recent architectures tend to embed fast on-chip memory with a dedicated allocator, such as the MCDRAM on Intel KNL and the hbwmalloc library [START_REF] Sodani | Knights Landing (KNL): 2nd Generation Intel R Xeon Phi processor[END_REF], which match our requirement of having a dedicated SMEM allocator.

c × k c from A[0][p c ] → A dma DMA block m c × n c from C[0][j c ] → C dma for i c in 0, ..., m -1 steps m c DMA next block m c × k c from A[i c + m c ][p c ] → A dma Wait A dma and pack A dma to A DMA next block m c × n c from C[i c + m c ][j c ] → C dma Wait C dma for j r in 0, ...,n c -1 steps n r // Macro-kernel for i r in 0, ...,m c -1 steps m r // m r × n r rank-k c micro-kernel : gemm ukr( A [i r ][0], B [0][j r ], C dma [i r ][j r ],
By instantiating control trees from the default ones, BLIS-RDMA inherits the same multi-threading mechanism as was enabled in the cache-based version, through five thread-count variables (e.g. BLIS JC NT, BLIS IC NT and so on) corresponding to the five loops around the micro-kernel (see Fig. 9.2). Upon available SMEM capacity, one can enable or disable multi-threading in some or a specific layer to fit the hardware constraints. For instance, setting BLIS JC NT = 2 will trigger two simultaneous DMA control flows in the outermost j c loop and double the SMEM footprint when each DMA flow is traversing the four underlying loops. A good parallelization scheme of BLIS-RDMA on embedded platforms would be to enable multi-threading in the inner loops to expose data-sharing. Typically, threads within the i c loop will share the same B dma panels (hence same B). More deeply, threads in the j r loop will share both B dma and A dma (hence same B and A respectively) as well as the C dma blocks.

Detailed memory footprint calculation of a top-down single-threaded DMA flow is given in Tab. 9.1. To take into account per-loop parallelization, let f i and t i denote respectively the unitary memory footprint and the thread count (BLIS * NT) of the i-th loop.

The total memory footprint F of a DMA control tree is determined by the following expression: F = t 5 (f 5 + t 4 (f 4 + t 3 (f 3 + t 2 (f 2 + t 1 f 1 )))) + f c (9.1)

We note f c the footprint of BLIS internal data and control trees (≤ 80 KiB) and take them in account in the calculation as these structures may be allocated in the SMEM for the execution performance as well. They are used only by the control code and are not duplicated like other user buffers (DMA, packing). The current DMA implementation has f 5 , f 2 and f 1 equal to zero. There is further possibility to implement fine-grained multi-level SMEMs by enabling the f 2 and f 1 footprints to fit the L2 and L1 cache for instance. Regarding f 5 , a simple multi-threading via the BLIS JC NT variable would be enough to enable coarse-grained parallelization.

Asynchronous DMA back-end interface

We present in this section a new back-end interface providing DMA support to perform asynchronous copy operations. The main idea behind this interface is inspired from the abstract design of BLIS, enabling the use of assembly-tuned micro-kernels written for a given architecture in order to reach near-peak performance. Likewise, Fig. 9.4 introduces a vendor-defined data structure and six functions of the back-end DMA interface. The Once the dma event t datatype and the back-end interface is mapped and implemented over the platform-specific DMA library, the BLIS-RDMA control trees will then appro- Upper-squarization, while ensuring the intended performance by keeping the same k c block size, leads to excessive SMEM allocation that can quickly exceed the capacity of the embedded hardware. On the contrary, lower-squarization results in under-utilization of the reserved SMEM by lowering the default k c block size, hence is exposed to performance degradation, especially on low-bandwidth memory systems. Mid-squarization, despite being more complex to implement, appears to be the most relevant option as it can yield acceptable performance while being reasonable in memory consumption, by maintaining a trade-off between computation and memory footprint. In the current state of our work, lower-squarization is employed for the sake of simplicity and proof of concept, and affects only the performance of SYMM. Implementing the mid-squarization is part of our roadmap.

Experimental results

Hardware configuration

We use the Kalray MPPA2-256 many-core processor as the target platform to validate our development. The cluster local memory of 2 MB is equivalent to the SMEM in BLIS-RDMA terminology. Multi-threading within a cluster is enabled using Pthreads. Default clock frequency of cores is set to 500 MHz with the peak performance of 2 GFLOPS/core in SGEMM and 1 GFLOPS/core in DGEMM. The DDR3 main memory is configured at SGEMM or DGEMM, multi-threaded or not). Any figure in which both SGEMM and DGEMM is presented will have the top border equal to the SGEMM peak and another lower line precising the DGEMM peak. 9.4.2 Multi-core level-3 BLAS and double-precision (FP64). The left-side lower-part non-transpose case for SYMM was intentionally chosen in our benchmark to measure the performance loss of the lowersquarization handling as previously discussed in section 9.3.4. Typically, single-core SSYMM and DSYMM achieve respectively 61% and 58% of peak. This represents approximately 13 to 16% of performance degradation between SYMM compared to other level-3

routines. However, we note that this degradation also depends on the block size configuration. It is significant on the MPPA2 processor because the default k c is twice larger than m c . Lower-squarization will thus divide the k c block size by two, hence double the required data bandwidth as well as the DMA control-flow overhead, which explains the performance loss on the MPPA2 processor. This leads to think that larger (and square) block sizes will be less exposed (or even not at all) to this issue when they are large enough to cross the memory bar of the system and pass to the compute-bound side.

Nevertheless, it is not seemingly the case for embedded platforms where local (on-chip) memory is often the most constrained resource.

Multi-threaded performance using 4 cores and 8 cores is presented in Fig. 9.6b and Fig. 9.6c. We observe that the efficiency decreases with increasing core count. We identify the following factors to the multi-threaded performance:

• Small block size due to limited local memory capacity. This constraint not only reduces the strong scalability of the parallelized j r loop, but also increases the I/O cost.

• Local memory bank conflicts. Micro-kernels written in MPPA2 assembly language employ a cache-bypass streaming load instruction between SMEM and core registers. This instruction delivers higher single-core throughput, but is sensitive to the local memory bank conflicts (local memory is composed of 16 banks interleaved every 64 bytes on the MPPA2 processor).

• Barrier and DMA synchronization. Thread barriers incur a significant overhead in multi-core context, especially in BLIS-RDMA where threads must synchronize between each other and wait for termination of DMA transfers.

Introduction 1 1 10 2 16 3

 11016 High Performance Computing: from Single-core to Many-core 3 1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2 Performance walls . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.2.1 Frequency-wall . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.2.2 Energy-wall . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.2.3 Memory-wall . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.2.4 Software-wall . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Lattice Boltzmann method (LBM) 11 2.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.1.1 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.1.2 Memory requirement . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.2 Propagation algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.2.1 One-step two-lattice (OT) . . . . . . . . . . . . . . . . . . . . . . . 14 2.2.2 One-step one-lattice (OO) . . . . . . . . . . . . . . . . . . . . . . . 15 2.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Basic Linear Algebra Subprograms (BLAS) 19 3.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3.1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3.1.2 General Matrix Multiplication (GEMM) . . . . . . . . . . . . . . . 21 3.2 BLAS-like Library Instantiation Software (BLIS) . . . . . . . . . . . . . . 23List of Figures1.1 Moore's law illustration . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Projected performance of the Top500 list . . . . . . . . . . . . . . . . . . 1.3 Processor-memory performance gap . . . . . . . . . . . . . . . . . . . . . 1.4 Arithmetic intensity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Roofline model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6 Different levels and latencies in memory hierarchy. . . . . . . . . . . . . . 2.1 LBM D3Q19 stencil. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 D3Q19 LBM applied on MPPA2 roofline model. . . . . . . . . . . . . . . 2.3 Propagation schemes of LBM in the D2Q9 sketch. . . . . . . . . . . . . . 2.4 Compressed-grid (shift) propagation algorithm. . . . . . . . . . . . . . . . 2.5 D2Q9 version of AA-pattern. . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Matrix multiplication. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Cache-based layer design of BLIS. . . . . . . . . . . . . . . . . . . . . . . 4.1 MPPA-256 processor overview. . . . . . . . . . . . . . . . . . . . . . . . . 4.2 OpenCL Data-Parallel: execution mapping. . . . . . . . . . . . . . . . . . 4.3 OpenCL Data-Parallel: memory mapping. . . . . . . . . . . . . . . . . . . 5.1 A part of MPPA Async API for 3D transfer. . . . . . . . . . . . . . . . . 5.2 Illustration of point2d t datatype for 2D copy. . . . . . . . . . . . . . . . 5.3 3D LBM/stencil decomposition. . . . . . . . . . . . . . . . . . . . . . . . . 5.4 3-depth pipeline (triple-buffering). . . . . . . . . . . . . . . . . . . . . . . 5.5 3D Row-major subdomain-indexing. . . . . . . . . . . . . . . . . . . . . . 5.6 Halo bandwidth ratio. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.7 Local/Remote copied index in 2D. . . . . . . . . . . . . . . . . . . . . . . 5.8 OPAL async vs. OPAL OpenCL on MPPA. . . . . . . . . . . . . . . . . . 5.9 Performance extrapolation of OPAL async. . . . . . . . . . . . . . . . . . 6.1 In-place updates on a grid using two-wall algorithm. . . . . . . . . . . . . 6.2 D2Q9 version of Two-wall in pull scheme. . . . . . . . . . . . . . . . . . . 6.3 D2Q9 version of three-wall in pull scheme. . . . . . . . . . . . . . . . . . . 6.4 Shared-memory OpenMP pseudo-code of Two-wall and Three-wall in D3Q19. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.5 Heterogeneous memory OpenCL implementation of Two-wall and Threewall in D3Q19. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.6 Comparison of different propagation algorithms in OpenMP on CPU and KNL. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Introduction This manuscript is the fruit of a three-year PhD enduring work on optimizing scientific applications on many-core processors. This thesis is founded by the CIFRE collaboration (French term of Convention Industrielle de Formation par la Recherche -Industrial convention of research-driven training, by the French Ministry of Research and Innovation) between the Kalray corporate and the Grenoble Informatics Laboratory (LIG) from the University of Grenoble Alps (UGA) and Centre for Energy and Thermal Sciences of Lyon (CETHIL) from the National Institute of Applied Sciences of Lyon (INSA Lyon), the National Center for Scientific Research (CNRS) and the University Claude Bernard Lyon 1.

Figure 1 . 1 :

 11 Figure 1.1: Moore's law illustration at the Heinz Nixdorf Museum. Credit: Paul Townend.

Figure 1 . 2 :

 12 Figure 1.2: Projected performance of the Top500 list.

Figure 1 . 3 :

 13 Figure 1.3: Processor-memory performance gap. Source: Computer Architecture: A Quantitative Approach -David A. Patterson and John L. Hennessy.

Fig. 1 .

 1 4 depicts the range of arithmetic intensity of several well-known computation kernels, such as sparse and dense linear algebra (DLA), lattice Boltzmann method (LBM) and fast Fourier transform (FFT).

Figure 1 . 4 :

 14 Figure 1.4: Arithmetic intensity. Source: https://crd.lbl.gov/departments/computer-science/PAR/research/roofline/.

Figure 1 . 5 :

 15 Figure 1.5: Roofline model. Source: https://alchetron.com/Roofline-model.

Figure 1 . 6 :

 16 Figure 1.6: Different levels and latencies in memory hierarchy.

Figure 2 . 1 :

 21 Figure 2.1: LBM D3Q19 stencil.

  methods are known for their low arithmetic intensity and particularly high memory bandwidth requirement. Taking the example of a basic LBM solver, depending on collision operator, between 200 and 400 floating-point operations are performed on a lattice node per time step. Most D3Q19 LBM implementations require storing all the 19 distribution values for each lattice node. A lattice domain L × L × L contains 19 × L 3 single-or double-precision floating-point numbers. Updating this lattice grid in a single time-step requires 19 × 2 × L 3 load/store memory operations for less than 400 × L 3 arithmetic operations. Thus, simulating the whole lattice domain through T time-steps will generate a huge amount of data movement of 19 × 2 × L 3 × T floating-point numbers for 400 × L 3 × T floating-point operations. Fig. 2.2 illustrates the memory-bound aspect of a D3Q19 LBM model.

Figure 2 . 2 :

 22 Figure 2.2: D3Q19 LBM applied on MPPA2 roofline model.

( a )

 a Push propagation. (b) Pull propagation.

Figure 2 . 3 :

 23 Figure 2.3: Propagation schemes of LBM in the D2Q9 sketch.

( a )

 a Even step. (b) Odd step.

Figure 2 . 4 :

 24 Figure 2.4: Compressed-grid (shift) propagation algorithm. Source: Wittmann et al. [27].

Figure 2 . 5 :

 25 Figure 2.5: D2Q9 version of AA-pattern with two kernels at odd and even time steps. The odd step reads local PDFs in their opposite order, collides and stores backlocally in natural order. The even step performs reads of PDFs from neighboring nodes (pull), collides and writes back (push) to the same place on these nodes, with opposite PDFs.

Figure 3 . 1 :

 31 Figure 3.1: Matrix multiplication.

Fig. 3 .

 3 2 depicts the global algorithm of the GEMM operation in BLIS. The three matrices A, B and C are partitioned and traversed through five loops around a micro-kernel. The micro-kernel performs a 4 th loop around micro-kernel 5 th loop around micro-kernel 3 rd loop around micro-kernel

Figure 3 . 2 :

 32 Figure 3.2: Cache-based layer design of BLIS. Courtesy of: Field G. Van Zee and Robert A. Van de Geijn.

6

  https://meltdownattack.com/ Any sufficiently advanced technology is indistinguishable from magic.

4 . 1 )

 41 embeds 256 VLIW compute cores grouped into 16 compute clusters (CC) and 16 system cores in two unified I/O subsystems (IOS). The sixteen compute clusters are organized in a 4 x 4 grid connected by a 2D torus Network-on-Chip (NoC). The processor delivers peak performance of 634 GFLOPS in single precision and 317 GFLOPS in double precision at a frequency 600 MHz, within a power consumption of 20 W.

Figure 4 . 1 :

 41 Figure 4.1: MPPA2-256 processor overview. (Source: Kalray).

4. 2

 2 Programming modelsThe complete Kalray Software Development Kit (SDK) features standard C/C++ and Fortran compilers, OpenMP/Pthreads support, low-level programming and communication libraries, as well as an OpenCL offloading runtime. They are divided into two main programming model: 1. A low-level distributed and inter-process communication (IPC) POSIX environment, supporting Fortran and C/C++. 2. A high-level host-based acceleration runtime based on the standard OpenCL 1.2 specification.

  function can be written in standard C/C++ or Fortran. It has direct access to the CC's SMEM in load/store model, and to the SMEM of other CCs or the DDR of IOS by explicit NoC communication, mainly based on software-triggered DMA engines to perform asynchronous transfers. Programming multiple clusters on MPPA in the POSIX model is similar to the one-sided message-passing model, largely used in HPC applications via the Message-Passing interface (MPI2) [9].

Figure 4 . 2 :

 42 Figure 4.2: OpenCL Data-Parallel: execution mapping. (Source: Kalray).

Figure 4 . 3 :

 43 Figure 4.3: OpenCL Data-Parallel: memory mapping. (Source: Kalray).

Fig. 4 .

 4 Fig. 4.2 and Fig. 4.3 depict respectively the execution mapping and memory hierarchy of the OpenCL Data-Parallel model on MPPA cores. Work-items within a workgroup is linearized and mapped on one PE. Collective operations like barrier() and async work group copy*() are accordingly handled by the runtime compiler to comply with desired execution order. Each work-group has load/store access to the global DDR memory with coherence based on DSM, and a private local memory ( local),

  t y p e d e f s t r u c t { i n t xpos ; i n t ypos ; i n t z p o s ; / * copy i n d e x * / i n t xdim ; i n t ydim ; i n t zdim ; / * b u f f e r d i m e n s i o n s * / } m p p a a s y n c p o i n t 3 d t ; / * 3D a s y n c h r o n o u s t r a n s f e r from remote t o l o c a l * / i n t mppa async memsget block3d ( v o i d * l o c a l , c o n s t v o i d * g l o b a l , s i z e t s i z e , i n t width , i n t h e i g h t , i n t depth , c o n s t m p p a a s y n c p o i n t 3 d t * l o c a l p o i n t , c o n s t m p p a a s y n c p o i n t 3 d t * r e m o t e p o i n t , m p p a a s y n c e v e n t t * e v e n t ) ;

Figure 5 . 1 :Figure 5 . 2 :

 5152 Figure 5.1: A part of MPPA Async API for 3D transfer. Prototype of get and put are similar.

5. 4

 4 Pipelined 3D LBM stencil on clustered many-core processors 5.4.1 Global algorithm

  Fig. 5.7). On boundary subdomains (e.g. subdomains 0, 1, 2, 3 of Fig. 5.7), the extended cuboid should be adjusted by applying a halo cutoff to deal with solid nodes. A local subdomain slot must therefore be allocated for F x × F y × F z nodes to match any cases.

Figure 5 . 3 :

 53 Figure 5.3: 3D LBM/stencil decomposition where Main-block subdomain (green) is copied with its surrounding halo layers (if exists) and one extra subdomain (blue) is needed to store post-collision state.

Figure 5 . 4 :

 54 Figure 5.4: 3-depth pipeline (triple-buffering) which allows 2-step distance between GET and WAIT, but only 1-step distance between PUT and WAIT, thus the PUT transfer will not be well overlapped (m: index of subdomain to compute, i: index of local buffer slot; G = GET; P = PUT; W = WAIT; C = COMPUTE; WCP = {WAIT + COMPUTE + PUT}; WG = {WAIT + GET}).

Figure 5 . 5 :

 55 Figure 5.5: 3D Row-major subdomain-indexing f : N → N 3 .

Figure 5 . 6 :

 56 Figure 5.6: Halo bandwidth ratio.Since the halo cell are needed for spatial dependency due to domain decomposition, they do not change the total number of updated cells nor the overall performance. As can be seen in Fig.5.6, the halo bandwidth costs up to 29% of data throughput on block size 16 × 16 × 16, but is enthusiastically reduced to as few as 10% on block size 64 × 64 × 64.This leads to think that further many-core architectures with larger local memories can noticeably improve 3D LBM performance by enlarging the subdomain size. The best performance is so achieved when the volume of the main block (C x × C y × C z ) is maximized and the HBW is minimized. Likewise, local storage should be reduced as much as possible. Let's assume single-precision floating-point representation, applying a D-depth pipeline for the two-local-cuboid method described above must fit into 1.4 MB of local memory and satisfy the linear-programming formulation below:

Figure 5 . 7 :

 57 Figure 5.7: Local/Remote copied index in 2D (in lattice node) with A: begin of the local buffer = (0,0); R: begin of the remote main block cuboid (without halo); B: begin of the copied cuboid (S), represented by: B a : index of S on local memory (from A) and B r : index of S on global memory (from R).

5. 5 . 1

 51 Pipelined 3D LBM stencil on MPPA We implement the pipelined 3D LBM algorithm on the MPPA-Bostan platform using the POSIX programming model and asynchronous 3D primitives from the MPPA Asynchronous One-Sided library. By default, MPPA-256 cores are set to run at 400 MHz and LP-DDR3 frequency is configured at 1066 MHz, i.e. ∼8.5 GB/s peak per DDR. Note that MPPA embeds two DDR interfaces (North and South) and the current OpenCL runtime only uses one DDR and exposes 1 GB of available global device memory, while the MPPA Asynchronous One-Sided library exposes both single and double DDR modes. Different cubic cavity sizes, varying from 64 to 224 are used in our tests, with some exceptions. Problem sizes larger than 160 cannot be run in OpenCL on MPPA due to the 1 GB device memory limit. Local work-group size in OPAL OpenCL is always set to 32 × 1 × 1, as it delivers the best performance in most of the cases.In single-DDR mode (POSIX and OpenCL), both LatticeEven and LatticeOdd are allocated on the North DDR. In double-DDR mode (POSIX-only), the LatticeEven buffer is allocated on the North DDR and the LatticeOdd is on the South DDR. The effective throughput of the double-DDR mode can be considered as twice as one of the single-DDR mode, thus 2× performance is expected. We present here results of the OPAL kernel rewritten with our new POSIX pipelined algorithm on the MPPA-256, called OPAL async, in 3-depth and 4-depth pipelines and following the local two-lattice method (S and S ) on various cavity sizes. These tests are further run in both singleand double-DDR modes. All these runs are checked for correctness against the original OPAL code on GPU.

Figure 5 . 8 :

 58 Figure 5.8: OPAL async vs. OPAL OpenCL on MPPA for duration = 1000 steps.Performance in MLUPS TODO

Figure 5 . 9 :

 59 Figure 5.9: Performance extrapolation of OPAL async on 8 × 8 × 8 subdomains with the first eight clusters correlation represented by a gray line for 1000 timesteps and cavity size 128.

Figure 6 . 2 :

 62 Figure 6.2: D2Q9 version of Two-wall in pull scheme: Copy operation (a) is represented by green cells. A barrier is needed to respect the read-after-write dependency between (a) and (b). Then, any node in the z-wall (black and blue cells) can be read (b), collided and stored in-place (c) independently.

Figure 6 . 3 :

 63 Figure 6.3: D2Q9 version of three-wall in pull scheme: Copy operation (a) is represented by green cells. The barrier is no longer needed, since FutureWall is only written (not read) in the current iteration. This removes the read-after-write dependency of two-wall, and exposes possibility of using streaming stores for further optimization. The black and the blue cells are then independently read (b), collided and stored in-place (c).

19 :

 19 end for (b) Three-wall.

Figure 6 . 4 :

 64 Figure 6.4: Shared-memory OpenMP pseudo-code of one time step with Two-wall and Three-wall in D3Q19. We assume AoS (Array of Structures) storage on CPU architecture. Fortran-like $OMP PARALLEL pragma is optional.

Figure 6 . 5 :

 65 Figure 6.5: Heterogeneous memory OpenCL implementation of one time step in Two-wall and Three-wall in D3Q19. Each $OMP PARALLEL in Fig. 6.4 is replaced by an L y × L x OpenCL kernel.

2 .

 2 Section 7.3 describes our MPI architecture design. Section 7.4 resumes our MPI implementation in pseudo-codes of blocking and non-blocking communication (MPI Send and MPI Isend). Some optimization ideas are then proposed and developed in this section such as (1) synchronizationfree eager send and (2) implicit local-buffered lazy send for short and medium sized messages respectively. A throughput estimation model based on the data transmission time is also introduced in section 7.5 to evaluate the communication performance. Section 7.6 presents our results for the ping-pong test following two scenarios, either symmetric ranks (MPI compute node -MPI compute node) or asymmetric ranks (MPI compute node -MPI I/O), corresponding on MPPA to CC-CC and CC-I/O subsystem respectively. Different optimization approaches are also tested and compared.

  supports up to 16 MPI ranks. Each MPI rank owns a private memory space of 2MB. Moreover, a hybrid MPI I/O rank is introduced running on the North IOS and manages the off-chip DDR memory. This MPI I/O rank is started from the host via the k1-mpirun command and is responsible for spawning MPI compute ranks on CCs subsequently. To keep the portability of any MPI legacy code, this extra MPI I/O rank is not listed in the MPI COMM WORLD. Any communication with this rank can be achieved through a local communicator (MPI COMM LOCAL ) that groups all MPI ranks within an

Figure 7 . 1 :

 71 Figure 7.1: Control message structures

Fig. 7 .Figure 7 . 2 :

 772 Fig.7.2 illustrates the structure of our MPPA-MPI implementation. Each rank emits control-message to its involved partner at each MPI call. The active server runs on the RM core and processes incoming control-messages. Furthermore, depending on MPI transactions and their status at runtime, the server can decide whether to perform a data send if this has not been or could not be done by the main thread, especially in case of a pending MPI Isend request or a matching registered lazy message.

4 : 5 :

 45 /* configure/start a non-blocking write (MPIDT layer) */ aio request ← configure aio write (req→buf, recv post.dnoc tag, ...);

  MPI rank 0 (CC 0) -MPI rank 15 (CC 15) and (2) MPI I/O 128 (IOS 128) -MPI rank 15 (CC 15). All MPI cluster ranks run at the same clock frequency of 400 MHz. The North IOS running the MPI I/O rank is configured to use the DDR controller at the default frequency of 600 MHz.

7. 6 . 1 15 Figure 7 . 3 :

 611573 Figure 7.3: Ping-pong throughput MPI Send (PE core) vs. MPI Isend (DMA).

Fig. 7 .

 7 Fig. 7.4c illustrates the message-rate of the four send methods. Not only this kind of measure gives user a high-level point of view about the implementation capacity to support communication load, but it shows interesting advantages of Eager and Lazy methods in tuning MPI applications, thanks to their fast sending time for short messages and synchronization-free algorithm.

( a )Figure 7 . 4 :

 a74 Figure 7.4: Optimization approaches comparison.

( a )

 a HPL benchmark on MPPA-256 against number of cores: BLAS vs.OpenBLAS.(b) HPL score extrapolation with increasing problem size.

Figure 7 . 5 :

 75 Figure 7.5: HPL current performance (a) and extrapolation (b)

  communication overhead. A more coarse-grained mapping, for instance, one MPI rank per MPPA would be a versatile solution, with the larger DDR space, PCIe and Ethernet interfaces. Each MPI rank then can be accelerated by the 16 compute clusters via OpenMP 4 or OpenCL in the MPI+X hybrid model. Chapter 8

Fig. 8 .

 8 Fig. 8.1). Each block-row and block-column contains a same number of K Kwg A i or B i blocks. 1 Our DMA-based GEMM algorithm, called GEMM-async, is given in Fig. 8.1. Each block C i is computed by one CC, equipped with 1.5 MiB of local memory (scratchpad) and 16 PEs. There are totally L = M Mwg × N Nwg blocks C to be distributed over 16 CCs. Each CC is responsible for at least L 16 blocks C. Whenever L is not multiple of 16, certain clusters must process one more block to fulfill remaining blocks (L mod 16).

Figure 8 . 2 :

 82 Figure 8.2: BLIS micro-kernel.

21 :else 23 : 1 (

 21231 C r [m r ][n r ] *= beta 22: Set registers C r [m r ][n r ] to zero 24: end if 25: acc[m r ][n r ] = alpha × acc[m r ][n r ] + C r [m r ][n r ] // m r × n r FMAs 26: Store acc[m r ][n r ] to C i [m r ][n r ] scratchpad In the Kalray OpenCL data-parallel mode, an independent compute-unit (work-group) possesses a maximal scratchpad of 64 KB ( local memory). We note that according to the OpenCL definition, scheduling of work-groups and their execution binding is controlled by the OpenCL runtime. It is not possible to synchronize work-groups or explicitly bind them on the same cluster to somehow build up and share a larger common local memory. Compared to the POSIX programming model, in which a compute-unit is mapped on a cluster with 1.5 MB of scratchpad memory (see Chapter 4, Section 4.2.page 29)), the OpenCL data-parallel mode has much less local memory available per compute-unit for DMA transfers. The same overlapping algorithm was ported from POSIX-C to OpenCL-C without any difficulty. However, blocksizes are reduced to M wg = N wg = 64, K wg = 48 for FP32, and M wg = N wg = K wg for FP64, due to the local memory constraint. Since the performance of the tiled algorithm relies strongly on the blocksize of the fast memory layer (see Chapter 3, Section 3.1.2.2 (page 22)), this memory constraint would suffer I/O bottleneck and deliver lower performance than the POSIX-C programming model. The BLIS micro-kernel was not usable either, since linking an assembly code or a user-compiled object to an OpenCL kernel was not supported at the time of our work.

Figure 8 . 3 :

 83 Figure 8.3: OpenCL POSIX-like: motivation.

8. 4 . 1

 41 Execution and memory mappingIn the POSIX-like mode, an OpenCL work-group is mapped on one compute cluster (CC) and operates 1024 KB of local memory. When a work-group is scheduled, its kernel is booted on the PE0 of the associated CC. The kernel, running on the PE0, then two OpenCL cl program's which are not self-standing. By performing a linking step with clLinkProgram(), the two binaries merge in a common and self-sufficient program, ready to be used to create an OpenCL cl kernel for the device.

Figure 8 . 5 :

 85 Figure 8.5: clLinkProgram(): linking a user-compiled object to an OpenCL kernel.

  ating at clock frequency of 400 MHz and DDR3 of 1066 MHz. At that frequency, the peak GEMM performance is 400 GFLOPS in FP32 and 200 GFLOPS in FP64. In this work, we present only performance of the FP32 precision to be in pace with an existing Kalray implementation of SGEMM, based on the name of another communication library MPPAIPC, considered as the baseline of our comparisons. The POSIX-C code is run in two configurations: (1) single-DDR with the three matrices A and B and C on one DDR, and (2) double-DRR with the matrices A and C on the DDR North and the matrix B on the DDR South. This configuration provides twice more DDR bandwidth to the A-B streaming algorithm. The OpenCL runtime presently enables the use of only one DDR for global memory, either in data-parallel or task-parallel mode. In the OpenCL task-parallel mode, we link the kernel to two BLIS micro-kernels, one using normal load through the L1 cache and one using streaming-load to prefetch data into registers. The both BLIS micro-kernels are written in assembly and compiled by gcc.

Tab. 8 .

 8 1 summarizes GEMM performance in various programming models and configurations on the MPPA2 many-core processor. We highlight the R max performance of 350 GFLOPS, (corresponding to 87% of R peak ) of the double-DDR POSIX-C implementation, combined with the BLIS micro-kernel that asynchronously reads data to registers by streaming-load and performs the 64-bit vectorized ffmawp instruction, the most intensive in the MPPA2 VLIW ISA. It yields 1.75x speedup over the baseline score of MPPAIPC, and 1.2x over the single-DDR POSIX-C configuration.The OpenCL data-parallel mode delivers 65 GFLOPS, due to the blocksize constraint as discussed in section 8.3. The OpenCL POSIX-like extension, with a large common scratchpad memory and capability of linking to a pre-compiled code, yields 207 GFLOPS which is three times better than the data-parallel mode and getting closer to the POSIX-C performance. Note that the current extension does not expose the same hardware resources as the POSIX-C model (one DDR versus two DDRs, 1024 KiB versus 1.5 MiB of local memory). An identical configuration will allows our OpenCL extension to reach the same GEMM performance of the POSIX model.

4 .

 4 Implementation of a reference DMA back-end with the memcpy function and Pthreads.

  [START_REF] Kurzak | Solving systems of linear equations on the CELL processor using Cholesky factorization[END_REF] [START_REF] Saxena | Optimization of BLAS on the Cell Processor[END_REF]. Lin et al.[START_REF] Lin | Optimizations of Two Compute-Bound Scientific Kernels on the SW26010 Many-Core Processor[END_REF] optimized the DGEMM operation on the Sunway many-core architecture. They achieved 88.7% peak efficiency by combining both asynchronous DMA between different memories and asynchronous RLC (Register-Level Communication) between cores within the same group. Tasende[START_REF] Tasende | Generation of the Single Precision BLAS library for the Parallella platform, with Epiphany co-processor acceleration, using the BLIS framework[END_REF] introduced a BLIS-instantiated Epiphany-accelerated SGEMM on the Parallela board, by offloading the micro-panel computation from the host processor to the Epiphany co-processor and overlapping the host-device communication.Each of these works comes with specific solutions suited to the architecture in terms of hardware topology, DMA characteristics and memory alignment, but the fundamental technique remains the same. While the main idea has always been using DMA engines to overlap communication and computation by streaming matrix sub-blocks to the local memory, the question of generalizing algorithms to a software API and a hardware abstraction layer is seldom brought up and tackled.Previous and comparable works to this chapter were on implementing the level-3 BLAS with DMA capabilities on Texas Instruments DSP by Igual et al.[START_REF] Francisco D Igual | Unleashing the high-performance and low-power of multi-core DSPs for generalpurpose HPC[END_REF] and Ali et al.[START_REF] Ali | Level-3 BLAS on the TI C6678 multi-core DSP[END_REF]. Later on, these techniques were applied on the BLIS framework[START_REF] Parikh | Implementation of Linear Algebra Libraries for Embedded Architectures Using BLIS[END_REF] [START_REF] Parikh | An Implementation of GEMM for DMA-enabled Architectures[END_REF] and released within the MCSDK HPC toolchain[START_REF]MCSDK HPC 3.x Linear Algebra Library[END_REF]. Likewise, Szydzik et al.[START_REF] Szydzik | Level-3 BLAS on Myriad multi-core media-processor SoC[END_REF] presented a level-3 BLIS port with DMA transfers on the Movidius Myriad, a cuttingedge low-power processor for computer vision, that appears quite similar to previous work on TI DSP.The aforementioned works illustrate the costly but fundamental need of providing an optimized BLAS library on novel and low-power many-core architectures. They leverage efficient DMA asynchronous communication between large (external, slow) and small (on-chip, fast) memories, which matches our first objective. However, the second objective of having a generic and portable DMA support for a BLAS library (e.g. BLIS framework) is not satisfied. Despite the fact that the TI's MCSDK HPC source code is public, development was done specifically for the C66x DSP without an abstraction layer. To the best of our knowledge, source code of the work on the Myriad architecture has not been released nor made open-source, as of today. None of other works proposes an open-source solution generic enough to be reused on another architecture. Finally, BLIS has been actively evolving since the last two years with significant changes and simplification in order to reach its current excellent maintainability. This provides an ideal opportunity for the design of a generic and portable DMA support in BLIS, so that developers of a new DMA-based architecture can quickly port a highly-optimized BLIS on their platform with minimal efforts. 9.3 Portable DMA support for level-3 BLIS 9.3.1 Algorithm overview To reduce the communication cost, DMA-based architectures feature one or several local memories close to computing cores. Each Synergistic Processing Elements (SPE) core of the IBM CELL platform possesses 256 KB of Local Store (LS) memory [98]. The TI KeyStone-II processor embeds 6 MB of Multicore Shared Memory (MSM) [101] between eight DSP cores and the host processor. The Movidius Myriad-2 processor implements 2 MB of shared Connection MatriX (CMX) memory between twelve SHAVE cores [105].

Fig. 9 .

 9 Fig. 9.1b depicts the global vision of the DMA support and how it is integrated into the layered design of the cache-based version (Fig. 9.1a) of BLIS. The three matrices A, B and C are partitioned and traversed through five loops around a micro-kernel, itself is a rank-k update constituting the sixth loop. Sub-partitions of A, B and C, before being used for computation by cores, are continuously streamed by DMA from the main memory to the SMEM. Level-3 in BLIS is currently built on top of GEMM (and TRSM)

  DMA-based BLIS layers.

Figure 9 . 1 :

 91 Figure 9.1: Cache-based and DMA-based layer design of BLIS. Image used and modified with permission.

Figure 9 . 3 :

 93 Figure 9.3: Illustration of a generic 2D copy. The receiving local buffer (right hand side) must be equal to (or can be larger than) the extent block (A dma or B dma ).

Figure 9 . 4 :

 94 Figure 9.4: BLIS-RDMA back-end interface, to be implemented by the hardware vendor in order to provide the generic 2D copy described in Fig.9.3. The dimension of the src (global) and dst (local) 2D buffer (xdim, ydim), the dimension of the extent block (width, height, element size) and its localization (xpos, ypos) within both src and dst buffer. The position offset, from which data is read/written, is then calculated from the start address of each src and dst 2D buffer.

Figure 9 . 5 :

 95 Figure 9.5: Illustration of DMA panel extension requirement for symmetrization in case of SYMM (left-side/lower/non-transpose A).

Fig. 9 .

 9 Fig.9.6 reports the performance of level-3 BLAS on a single MPPA2 compute cluster with different core counts. As can be seen from Fig.9.6a, most single-core operations (except SYMM) obtain between 68% and 75% of peak in both single-precision (FP32)

Figure 9 . 6 :

 96 Figure 9.6: Performance of multi-threaded level-3 BLAS on MPPA2-256, single compute-cluster, 500MHz.

9. 4 . 3

 43 Multi-core xGEMM in different shapesWe present in this section performance of GEMM in single and double precision in three matrix configurations: (1) square matrices (m = n = k), (2) various m = n and k = k c and (3) various k with constant m = n. Each configuration is run with 1, 4 and 8 cores and performance is normalized in terms of GFLOPS/core. Satisfactory strong scalability is achieved when the GFLOPS/core ratio does not depend from the number of cores used. In general, the multi-threaded performance on MPPA tends to yield lower GFLOPS/core ratio, for the reasons mentioned in the previous section. The only possible improvement in our opinion is to increase the SMEM footprint to enlarge the DMA panel block size, which is constrained by hardware limitations.Performance efficiency and scalability of configuration (1) and (2) are reported in Fig.9.7a-9.7b and Fig.9.7c-9.7d, in FP32 and FP64, respectively. We observe slightly higher and more stable performance of the k = k c configuration compared to the square (m = n = k) configuration. When k is not multiple of k c , the p c loop (4 th loop around the micro-kernel) must perform an extra iteration to process the trailing elements in the k-direction. This implies an additional cost in re-triggering the DMA transfers and overlapping mechanism in the inner loops to update a smaller number of elements. Fortunately, this overhead can be compensated by the computational load in the case of large k (Fig.9.7a-9.7b).

  SGEMM (m = n = k). DGEMM (m = n = k). DGEMM (k = k c = 130).

  SGEMM (m = n = 4000).

  DGEMM (m = n = 4000).

Figure 9 . 7 :

 97 Figure 9.7: Scalability of multi-threaded xGEMM in different shapes on MPPA2-256, single compute-cluster, 500MHz.

  In Gerhard R. Joubert, Hugh Leather, Mark Parsons, Frans J. Peters, and Mark Sawyer, editors, Parallel Computing: On the Road to Exascale, Proceedings of the International Conference on Parallel Computing, ParCo 2015, 1-4 September 2015, Edinburgh, Scotland, UK, volume 27 of Advances in Parallel Computing, pages 113-122. IOS Press, 2015. 2. Julien Hascoët, Benoît Dupont de Dinechin, Pierre Guironnet de Massas, and Minh Quan Ho. Asynchronous one-sided communications and synchronizations for a clustered manycore processor. In Proceedings of the 15th IEEE/ACM Symposium on Embedded Systems for Real-Time Multimedia, ESTImedia 2017, Seoul, Republic of Korea, October 15 -20, 2017, pages 51-60, 2017. 3. Minh Quan Ho, Christian Obrecht, Bernard Tourancheau, Benoit Dupont de Dinechin,

1. Minh Quan Ho, Bernard Tourancheau, Christian Obrecht, Benoît Dupont de Dinechin, and Jérôme Reybert. MPI communication on MPPA many-core NoC: design, modeling and performance issues. and Julien Hascoet. Improving 3D lattice Boltzmann method stencil with asynchronous transfers on many-core processors. In 2017 IEEE 36th International Performance Computing and Communications Conference (IPCCC) (IPCCC 2017), San Diego, USA, December 2017. 4. Minh Quan Ho, Christian Obrecht, and Bernard Tourancheau. New parallel inplace update algorithm for better memory usage in 3D lattice Boltzmann algorithm. In submission, 2017. 5. Minh Quan Ho, Benoit Dupont de Dinechin, Bernard Tourancheau, and Christian

  The libflame project implements LAPACK-like features and lies closely on BLIS for the BLAS support. Both BLIS and libflame are developed by the Science of High-Performance Computing (SHPC) group at the University of Texas at Austin.

	They are released under the open-source BSD 3-clause license, facilitating adoption by
	industry.
	During many years, the research community was missing a well-structured, open-source,
	light-weight, portable and high performance BLAS library. Proprietary implementa-
	tions are considered as black-boxes and platform-specific. Open-source libraries like
	ATLAS and OpenBLAS appear too bloated or difficult to port and to optimize on a
	new architecture. Researchers and vendors need an easy and extensible framework as
	an experimental tool, not only to implement and tune new linear algebra algorithms,
	but also to maintain and optimize easily BLAS functions on next-generation architec-
	tures. These crucial points have been tackled and successfully solved within the BLIS
	framework for conventional cache-based CPUs.

  pipeline, restricts to a very small search domain (C d ≤ 128) that can be resolved by running the branch-and-bound algorithm in a script. Solutions can either be (D, C x , C y , C z ) =[START_REF] Ho | Improving 3D lattice Boltzmann method stencil with asynchronous transfers on manycore processors[END_REF][START_REF] Wunderlich | Efficacy and efficiency of algorithm-based fault-tolerance on GPUs[END_REF] 8,[START_REF] Wunderlich | Efficacy and efficiency of algorithm-based fault-tolerance on GPUs[END_REF] with 36% HBW ratio or (D, C x , C y , C z ) = (4, 8, 8, 16) with 43% HBW ratio. A permutation of C x , C y , C z also gives other satisfactory solutions, with the same

Table 5 . 1 :

 51 Copied index offset and halo cutoff of a subdomain.

  It is worth mentioning that KNL embeds 16 GB of on-chip high-bandwidth MCDRAM accessible either as cache or flat memory. In this test, we explicitly allocate lattice buffers directly into the MCDRAM via the hbwmalloc library[START_REF] Sodani | Knights Landing (KNL): 2nd Generation Intel R Xeon Phi processor[END_REF].1

• 64-core MIC Intel Xeon Phi Knights Landing (KNL) 7230 at 1.30 GHz

Table 6 . 2 :

 62 Compilation flags and OpenMP context on CPU and KNL.

	Arch.	Compilation	OpenMP env.
	CPU	icc -O3 -qopenmp -align -fma -ftz -finline-functions	OMP_NUM_THREADS=24 OMP_PROC_BIND=close OMP_PLACES=cores
		icc -O3 -qopenmp -align	OMP_NUM_THREADS=64
	KNL	-fma -ftz -finline-functions	OMP_PROC_BIND=close
		-xMIC-AVX512	OMP_PLACES=cores

Table 6 . 3 :

 63 OpenCL drivers and configuration on CPU and GPU.

	Arch.	OpenCL driver	OpenCL config. in 2D (or 3D)
	CPU	Intel OpenCL 14.2	GlobalSize L x × L y (×L z ) LocalSize 32 × 1 (×1)
	GPU	CUDA 7.5.18	GlobalSize L

x × L y (×L z ) LocalSize 32 × 1 (×1)

  to evaluate the maximum communication latency[START_REF] Bruno D'ausbourg | Deterministic Execution on Many-Core Platforms: application to the SCC[END_REF] also incurs unordered packets, hence inappropriate for hard real-time applications.

	The SCC-specific MPI-like native communication library (RCCE) delivers peak through-
	put of 55 MB/s on the ping-pong test [83]. By the same test, Clauss et al. [84] presented
	iRCCE (an improved RCCE version) and SCC-MPICH (an MPICH-based implemen-
	tation over iRCCE) that reach respectively 150 MB/s and 120 MB/s of throughput.
	RCKMPI, an Intel MPI implementation for SCC is also bounded by the performance of
	the iRCCE layer. Our MPI library on MPPA was built from scratch over the MPPAIPC
	library, without any TCP/UDP layer, while an MPICH-based solution would not fit the
	cluster private memory space (2 MB). A such MPICH implementation on MPPA can
	be extrapolated to 1.0 GB/s by adding the same overhead of 20% of SCC-MPICH over
	to iRCCE.
	Intel Single-Chip Cloud (SCC) is a prototype aimed to promote many-core proces-
	sor. Its 48 cores are organized in 24 dual-core clusters with access to off-chip DRAM
	shared/private region for all/each core through a look-up table (LUT), also a dedicated
	shared on-chip Message Passing Buffer (MPB). This memory architecture gives extra ad-

are mainly used in high performance embedded systems such as networking and multimedia. The TilePro64 processor defines a flat 2D-mesh with 64 identical VLIW cores connected through the Tilera iMeshTM network-on-chip. Cache coherence on TilePro64 is guaranteed by a hardware mechanism called Dynamic Distributed Cache (DDC) [69]. Kang et al. [81] propose an MPI implementation on Tile64 processor which delivers up to 250 MB/s on MPI Send/MPI Recv communication, with the largest message size of 256 KB due to the limited memory per core. At this buffer size, our MPI implementation on MPPA delivers 400 MB/s on MPI Send or up to 1 GB/s using DMA. vantage for implementing quick message sending based on shared buffers. However, the use of the dynamic NoC routing on Intel SCC (instead of static NoC in MPPA) makes it difficult Intel Many Integrated Core (MIC), known as Intel Xeon Phi co-processor family, is a x86-based many-core architecture with native support of Linux operating system and standard software stack. The first Intel MIC generation, namely Knights Corner

  on each MPI rank (CC and IOS) to handle incoming control messages from all other ranks (including itself on loop-back). Upon control message arrival, a callback function is executed on the RM, consisting typically on saving it into an internal buffer which later will be read by MPI calls from the main function (PE0).

	Control messages exchanged in the MPIC layer contain either one of the structures
	defined in Fig. 7.1.
	i n t	m pi t a g ;	/ * MPI message t a g * /
	. . .
	} s e n d p o s t t ;
	/ * Message s e n t by Rx t o Tx ( Clear -To-Send ) * /
	t y p e d e f s t r u c t r e c v p o s t s {
	i n t	d n o c t a g ;	/ * DNoc a l l o c a t e d on Rx * /
	mppa pid t	r e a d e r i d ;	/ * ID o f Rx p r o c e s s * /
	i n t	m pi t a g ;	/ * MPI message t a g * /
	. . .
	} r e c v p o s t t ;

/ * Message s e n t by Tx t o Rx ( Request-To-Send ) * / t y p e d e f s t r u c t s e n d p o s t s { mppa pid t s e n d e r i d ; / * ID o f Tx p r o c e s s * /

Table 7 . 1 :

 71 MPI send/receive implementation in MPIDT level

	MPPA-MPI	MPPAIPC Portal
	MPI Send, MPI Ssend	mppa pwrite,
		mppa pwrites
	MPI Isend, MPI Issend mppa aio write
	MPI Recv, MPI Irecv	mppa aio read
	MPI Wait	mppa aio wait

  /* send recv post to real source (MPIC layer) */ 12: send recv post(my rank, dnoc tag, real source, tag, comm, ...); 13: /* wait data (MPIDT layer) */ 14: mppa aio wait(aio request); 15: return MPI SUCCESS;

	7:	real source ← send post.source;
	8: else
	9:	real source ← source;
	10: end if
	11:

1: my rank ← get rank(comm); 2: dnoc tag ← allocate dnoc tag(); 3: /* configure to receive data on this dnoc tag */ 4: aio request ← configure aio read(buf, dnoc tag, ...); 5: if source == MPI ANY SOURCE then 6:

send post ← find send post(count, datatype, tag, ...);

1 :

 1 my rank ← get rank(comm); 2: /* send send post to dest */ 3: send send post(my rank, dest, tag, comm, ...);

	4: /* wait for matching recv post from Rx (MPIC layer) */
	5: repeat
	6:

1 :

 1 my rank ← get rank(comm); 2: /* send send post to dest */ 3: send send post(my rank, dest, tag, comm, ...); 4: req ← new request(buf, count, ..., PENDING ISEND); 5: /* look for a matching recv post (MPIC layer) */ 6: recv post ← find recv post(count, datatype, tag, ...); 7: if recv post = NULL then aio request ← configure aio write (buf, recv post.dnoc tag, ...);

	8:	/* configure/start a non-blocking write (MPIDT layer) */
	9:
	10:	req→status := STARTED;
	11:	req→aio request := aio request;
	12: else
	13:	/* Do nothing (request initialized NON STARTED) */
	14: end if
	15: request ← req;
	16: return MPI SUCCESS;
	Algorithm 5 callback recv post(recv post)

1: /* look for a matching pending isend */ 2: req ← find pending isend(recv post); 3: if req = NULL then

Table 7 . 2 :

 72 Transmission time (cycles).

		To	MPI Send	MPI Isend
	CC 0	CC 15	t = 0.98 × N + 31430	t = 0.27 × N + 33690
	CC 15	CC 0	t = 0.98 × N + 30240	t = 0.27 × N + 32850
	IOS 128 CC 15	t = 13.52 × N + 159544 t = 0.84 × N + 181300
	CC 15	IOS 128 t = 0.98 × N + 129200	t = 0.26 × N + 144500

  OpenCL Data-Parallel In OpenCL, blocks C i are distributed onto work-groups. Each work-group computes one block C i . Since GEMM is inherently parallel and a work-group is mapped on one PE on Prefetch first column A p [0*cs] → registers A r [m r ] 3: Prefetch first row B p [0*rs] → registers B r [n r ] 4: Set registers acc[m r ][n r ] to zero ][n r ] += A r [m r ] × B r [n r ] // m r × n r FMAs acc[m r ][n r ] += A r [m r ] × B r [n r ] // m r × n r FMAs

	Algorithm 6 BLIS gemm micro-kernel on MPPA2 with in-core computation and reg-
	ister prefetching.
	1: /* Prolog */
	2: 5:
	6: /* k c -loop */
	7: for k in 0 ... k c -2 do
	8:	Prefetch (k +1)-th column A p [(k +1)*cs] → registers A r [m r ]
	9:	Prefetch (k +1)-th row B p [(k +1)*rs] → registers B r [n r ]
	10: acc[m r 11: A r ← A r
	12:	B r ← B r
	13: end for
	14:
	15: /* Epilog */
	16: 17:
	18: /* alpha-beta post-processing */
	19: if (beta != zero) then
	20:
	MPPA2, we set the local work-group size to 1×1 and global work-items count to 256×1,
	so that there are exactly 256 work-groups scheduled on 256 PEs during execution. The
	asynchronous GEMM pipeline on each work-group will operate the longest path, necessary
	to amortize the sequential cost of prolog and epilog.

Table 8 . 1 :

 81 GEMM-async performance on MPPA2-256, FP32, frequency 400 MHz.

		MPPAIPC	POSIX	POSIX	OpenCL	OpenCL
		(baseline)			Data-Parallel Task-Parallel
	Matrix size	4096 × 4096	4096 × 4096	4096 × 4096	4096 × 4096	4096 × 4096
	Compiler	gcc (-O3)	gcc (-O3)	gcc (-O3)	clang (-O3)	clang + gcc
	#DDR	2	1	2	1	1
	Local memory	1.5 MiB	1.5 MiB	1.5 MiB	64 KiB	1024 KiB
	C standard	N/A	75 GFLOPS	76 GFLOPS	65 GFLOPS	N/A
	Assembly: ffmawp, 200 GFLOPS 227 GFLOPS	228 GFLOPS	N/A	192 GFLOPS
	cached-load
	Assembly: ffmawp,	N/A	290 GFLOPS 350 GFLOPS	N/A	207 GFLOPS
	streaming-load

  ...) Pseudo-code of BLIS-RDMA through the five layers. Green is used for matrix B. Blue is used for matrix A. Red is used for matrix C. Swap of A dma /B dma /C dma buffers within each iteration is not presented here.

	end for
	end for
	DMA block m c × n c C dma back to main memory
	end for
	end for
	end for
	Figure 9.2:

Table 9 . 1 :

 91 Single-threaded level-3 BLIS-RDMA SMEM footprint calculation.

http://www.netlib.org/blas/blasqr.pdf

n + 2 multiplications and n additions

2n loads for Ai,t and Bt,j, one load and one store for Ci,j

Included in the IEEE 754-2008 standard and largely available on modern CPU and GPU architectures. The FMA instruction has advantage of significantly reducing the number of CPU cycles and minimizing the accumulated error due to successive rounding steps.

This layout is similar to the row-major GEMM TN format: A transposed and B non-transposed.

The ceil() function rounds up to the next integer to deal with trailing matrix elements after decomposition.

An edge block is the block added by the ceil() function, often at the right and bottom of the decomposition grid, to cover trailing matrix elements.

In BLIS, a packing step is carried out on Ap and Bp to pack these panels into other temporary contiguous ones ( A and B). These packed panels reduce the memory spatial distance between each rank-1 update (as small as mr and nr), thus improve cache and TLB performance.

Acknowledgements

Two-wall propagation algorithm

In this section, we present a new parallel one-lattice LBM propagation algorithm, named two-wall algorithm, in which the lattice grid G is updated in-place. The code change for adapting this algorithm on complex boundary conditions remains simple, easy to understand and to implement.

Algorithm

We propose a new propagation algorithm enabling in-place update when using the pull scheme. Note that while we are assuming the pull scheme in this section, applying the push scheme is entirely possible as well. The main idea for the pull scheme is, before colliding and streaming new PDFs, outdated PDFs must be copied out and saved in a temporary buffer (small, 2D buffer). Any spatial dependency which requires reading data from a specific node in the grid G, should read in the temporary buffer instead. More in detail, let say that the grid G is allocated in the C99 multi-dimensional array convention:

We decompose the grid G into L z walls. Each wall is composed of L y × L x nodes. The lattice grid G is updated wall-by-wall. There will be a sequential for loop in the z-direction to sweep through the L z walls of G. At each ziteration, the wall z (referenced by Grid[z]) will be updated from time step t to t+1 and its new PDFs are written in-place into Grid [z]. We then introduce two wall buffers, a past wall and a current wall, each of L y ×L x nodes (float Walls [START_REF] Hascoët | Asynchronous one-sided communications and synchronizations for a clustered manycore processor[END_REF]

These two walls, swapped from one time step to another, are respectively used to store old PDFs of the wall z -1 and z before they were overridden by the new PDFs of the collision operation. Wall buffers can be allocated with or without ghost-layer, depending on the application. an alternative to other one-lattice algorithms. The barrier cost and OpenMP overhead of three-wall is acceptable in a high-concurrency context like on KNL, with 64 cores and a high-latency cache coherence protocol. Three-wall can also be optimized on the KNL architecture by generating non-temporal streaming stores via OpenMP hints and/or compiler flags, which we expect to yield higher performance than AA-pattern.

OpenCL

The OpenCL code is run on two platforms:

• CPU Intel Xeon Haswell E5-2680v3 (same as the OpenMP code)

The CPU is used as an OpenCL device by setting the device type to CL DEVICE TYPE CPU and similarly CL DEVICE TYPE GPU for the GPU. Two-lattice has been already implemented in the original version of OPAL [START_REF] Obrecht | Performance Evaluation of an OpenCL Implementation of the Lattice Boltzmann Method on the Intel Xeon phi[END_REF]. We took this code as a basis to write kernels for AA-pattern, two-wall and three-wall by adapting the propagation step and which is 10 times higher than on the CPU. Besides, we obtain about 1300 MLUPS for two-wall and 1380 MLUPS for three-wall, which is 6% better than two-wall (see Fig. 6.7b). For unknown reason, the AA-pattern implementation fails to execute on problem sizes larger than 352 3 and returns an XID 31 error. According to the NVIDIA documentation, this error is due to an illegal address access. It is possibly caused by the application code, but there could also be driver bugs or hardware bugs. As the issue does not appear on problem sizes smaller than 352 3 with the same code, we suspect a pathological corner-case of the AA-pattern algorithm on the target GPU architecture and/or the current driver version. Nonetheless, the original source of the error remains unknown for the moment. Regarding the LUPS per byte efficiency on GPU, in Fig. 6.7d, the perf-mem ratio of AA-pattern is two times higher than the one of two-lattice (same as in OpenMP). Furthermore, we see that two-wall and three-wall algorithms are less efficient than two-lattice, by their perf-mem ratio between 0.6 and 0.9, normalized to two-lattice. From these results, we conclude that two-/three-wall can also be easily implemented from an existing two-lattice kernel in OpenCL. While running correctly on CPU or GPU in OpenCL, two-/three-wall do not offer significant performance gain nor memoryefficiency compared to two-lattice and AA-pattern. In our testing environment, we notice Chapter 7

Message Passing Interface (MPI)

on Many-core Processors

The single biggest problem in communication is the illusion that it has taken place.

-George Bernard Shaw. two blocks in prefetch by DMA (one next A i and one next B i ), and other three blocks in computation (current A i , B i and C i ). Blocksize parameters (M wg , N wg and K wg ) are consequently set to M wg = N wg = K wg = 256 in FP32 and M wg = N wg = 256, K wg = 128 in FP64, so that the scratchpad memory can hold at least five blocks: two A i , two B i and one C i . We note that, in our algorithm, both single-buffering and double-buffering were implemented on blocks C i . The communication cost on C i (read at the beginning and write at the end), even being on the critical path with the single-buffering scheme (not overlapped), can be amortized by flops computation from blocks A i and B i when

On edge blocks 2 , since the matrix data does not entirely fulfill the block, we perform zero-padding on the remaining part (right and bottom) of the edge block, so that the block-computation kernel can be reused without change or polluting results.

In order to facilitate development process and to focus on algorithmic aspects, we de- From the ease-of-coding point of view, the OpenCL POSIX-like mode can be considered as an alternative deployment runtime to the POSIX programming model. The OpenCL kernel, once booted on the PE0, is identical to a main function in the POSIX-C model.

Developers no longer need to take care of writing compilation Makefile, nor launching scripts from the host processor onto the MPPA processor. Instead, they can directly provide a kernel wrapper to the OpenCL program and get it compiled and linked to an optimized object or library implementing the crucial part of performance (BLIS microkernel). The executable is then sent, scheduled and started on the MPPA (operations that are considered as non-productive work) by the OpenCL driver. The on-site runtime only need to provide essential features of the libc, such as malloc/free, printf and pthread *, as well as access to some vendor low-level tools like the mppa-async communication library or profiling APIs, for the kernel to be fully operational.

Integration of object code or library

Special cases handling

We discuss in this section the technical solutions that were implemented in BLIS-RDMA in order to comply with the limited local memory available on DMA-based platforms.

These restrictions may or may not apply on other (conventional) architectures where hardware resource is less critical. In any case, these operations can be disabled in the source code.

For TRMM (Triangular matrix multiply) and TRSM (Triangular equation system solving) operations, the long dimension of the micro-panel (k c ) is, if necessary, rounded up to be multiple of m r or n r (upon sidea) [START_REF] Field | BLIS: A framework for rapidly instantiating BLAS functionality[END_REF] so that the packing facility can manage edge cases by zero-padding. 4 Despite the fact that only a few elements are added to the DMA panel, this slightly increased k c size can exceed the available SMEM space reserved to the library, which will likely fail at execution. In BLIS-RDMA, k c will be rounded down to the lower nearest multiple of m r (or n r ), instead of the upper nearest value. As m r and n r are often small (4, 8 and less than 32) and k c is in the order of hundreds, this rounding-down mostly does not affect the global performance.

For SYMM (Symmetric matrix multiply) operations, let us denote A the symmetric matrix and consider the left-side lower-part non-transpose case with k c larger than m c as shown in Fig. 9.5. The packing process of a diagonal-intersected panel must perform a symmetrization -a mirror-copy of the lower part and symmetrically write to the upper part of the panel so that the gemm micro-kernel can be reused for computation. Since the SMEM address space is different from the main memory and the DMA panel of size m c × k c is physically detached and copied from the main memory to the SMEM, the packing routine will try to read data from an invalid region outside the DMA panel (the gray part in Fig. 9.5), which will probably cause a segmentation fault. There are three options to manage this case: (1) virtually increase the m c dimension of the DMA panel to be equal to k c (without changing the real m c step parameter) to cover the gray missing 4 Edge case occurs when a panel intersects the diagonal and thus is divided into two parts: a gemmlike part which is computation-ready and a trmm-or trsm-like part in which a small triangle needs to be zeroed out before the computation. The first part can be used as is as input of the gemm microkernel. Whereas on the second part, it is preferable to have kc multiple of mr (or nr) so that the trailing triangular sub-matrix is compatible with the trsm micro-kernel. The two parts can also be computed at once by the fused gemmtrsm micro-kernel as well.

1333 MHz with a theoretical bandwidth of 10.4 GB/s. For the time being, BLIS-RDMA is run on a single compute cluster and multi-threaded over eight cores. The eight remaining cores are not used, due to local memory capacity constraints, and will be enabled for BLIS when a planned more memory-efficient software code cache is released.

The back-end DMA interface is implemented on top of the MPPA2 asynchronous library in less than 70 lines of code (see Appendix A.2). In order to maximize data-sharing and minimize SMEM footprint, multi-threading is enabled only in the j r loop by setting the BLIS JR NT environment variable to the thread count. This is considered as the most modest mode in BLIS-RDMA, as all threads share the same A dma , B dma , A, B and C dma buffers. Blocksize parameters m c , n c , k c and register block size m r , n r used on MPPA2 are given in Tab. 9.2, on which we also calculate the SMEM footprint for compilation. This estimated SMEM footprint is calculated from formulas given in Tab. 9.1 and Eq. 9.1. This size bounds closely the real execution footprint, without the need of beforehand iterative run. Micro-kernels of gemm in single and double precision are written in the MPPA2 VLIW assembly language to maximize flops count, which we estimate to reach 97% efficiency. The fused gemmtrsm micro-kernels are combined from gemm micro-kernels and reference C99 trsm ones.

In the current state of work, the level-3 BLAS is run only in single and double precision on the MPPA2 processor, hence without c and z operations and Hermitian routines. With on-going software developments, we expect to be able to run the full BLIS testsuite on the MPPA2 processor in a near future. In the meantime, BLIS-RDMA has already passed the testsuite on CPU with the back-end DMA interface emulated by memcpy and Pthreads, which we believe to be a pertinent cross-validation. Comparing to the cache-based implementation on MPPA2, the DMA-based approach delivers between 5× to 10× speedup on small matrices (≤ 512). On larger matrices, the cache-based solution suffers severe performance degradation (up to 1000× slowdown), due to the increasing overhead of matrix-packing, caused by large strided-accesses and cache-thrashing on the software cache. In the case of small k however, as depicted in Fig. 9.7e-9.7f, when k is just slightly larger than k c , the performance drops by 8% on the first k c slump (130) and by 3% on the subsequent slumps (N ×k c ) This phenomenon, called divot, has been observed and cured by Smith et al. [START_REF] Tyler | Anatomy of high-performance many-threaded matrix multiplication[END_REF] in the cache-based version of BLIS on the Xeon Phi processor, by increasing the effective k c by a certain size at each k c -multiple point to create a "bridge" passing over the divot value. Interestingly enough, contrarily to the previous situation, this performance-throttling (or divot) on the MPPA2 processor cannot be quickly recovered to its previous sustained performance, but rather climb slowly up and re-reach the maximum performance exactly at the next k c -multiple point, and get throttled down again (like a TCP congestion line). We believe this is a representative case of the difference of overhead between a software DMA-based communication versus a hardware cache-based implementation, where a cache-miss takes often less than 50 cycles, while a "DMA-miss" is in the order of thousands of cycles.

Conclusions

Pursuing the previous chapter work on optimizing blocked GEMM on DMA-based manycore processor, in this chapter, we present BLIS-RDMA, a generic support for, not limited to, the level-3 BLAS in the BLIS framework. We achieve three objectives: [START_REF] Ho | MPI communication on MPPA many-core NoC: design, modeling and performance issues[END_REF] standard programmability by supporting the BLAS API, with (2) portability through a back-end interface, and (3) computation performance in leveraging the cumbersome DMA capabilities on an increasing range of embedded and non-conventional many-core architectures. The work are demonstrated on the MPPA2 many-core processor with satisfying performance (75% peak on single-core). The implementation, despite being presently ported only on MPPA and CPU, is designed having in mind to support any DMA-based architecture and to minimize the porting effort through a back-end interface of six functions, aimed to be easily implemented by hardware vendors. We observe that, even with a DMA engine that enables efficient overlap between communication and computation, it is still hard to obtain near-peak performance on embedded platforms that are constrained by the local memory capacity.

Conclusions

Home isn't where you're from, it's where you find light when all grows dark.

-Pierce Brown, Golden Son.

A few steps back

In this thesis, we walk through memory-bound and compute-bound applications on increasingly popular conventional and non-conventional many-core architectures. Main contributions of this thesis are summarized as follows.

We first propose approaches for improving the data bottleneck of the 3D lattice Boltzmann method on many-cores processors by using scratchpad memory and asynchronous DMA transfer. We achieve 33 % performance gain on the MPPA architecture by actively streaming data from and to local memory, unlike in the passive OpenCL programming model.

We then tackle the memory consumption of the LBM by proposing two novel algorithms which require only one lattice array. The proposed algorithms are implemented in OpenMP and OpenCL and offer more development facility than other algorithms among the one-lattice class. Results show the adequacy of using our algorithms on a large set of cores with a sophisticated cache and memory system, by obtaining the same performance as the AA-pattern algorithm on the Xeon Phi Knights Landing processor, and 1.5 times higher memory-efficiency than the state-of-the-art two-lattice algorithm.

We present a DMA-based matrix multiplication (GEMM) algorithm. We identify necessary key features of DMA-based programming models to obtain high performance and portability. We implement those approaches on MPPA in POSIX-C and OpenCL with a hand-tuned optimized VLIW assembly kernel. We achieve up to 350 GFLOPS of SGEMM (87 % peak performance).

We present BLIS-RDMA, a generic DMA support for the BLIS framework. BLIS-RDMA enables high-performance BLAS computations on DMA-based many-core processors.

BLIS-RDMA proposes a deterministic data footprint model for a controlled scratchpad allocation. BLIS-RDMA implements a portable asynchronous one-sided communication interface which allows porting of highly optimized BLAS library on any DMA-based architecture in less than 100 lines of code. Our implementation is validated on the MPPA processor, delivering over 75 % peak performance within a memory footprint of 480 KB.

These contributions are built around and for the DMA capabilities of non-conventional many-core processors in general, and the Kalray MPPA in particular. Providing a coherent application mapping on hardware and partitioning the dataset onto parallel compute-units with asynchronism and concurrency is the most challenging task.

Perspectives

Next step in our work plan is porting and improving proposed LBM in-place propagation algorithms with DMA engines to prefetch and copy wall data on the MPPA processor.

Extending them to a 2.5D approach and complex geometries is also a plausible direction.

Regarding BLIS-RDMA, we want to study the adaptation cost to run it on the whole MPPA processor, as well as the portability on other DMA-based platforms or high-end many-core processors. There will be implementation of the mid-squarization approach for optimal SYMM performance, as well as integration of the DMA capabilities to the lower (and easier) level-1 and level-2 in order to provide a full-blown BLIS-RDMA. We also plan to support more elaborate control trees for multi-level SMEMs upon hardware availability. We will also look at extending the DMA back-end to perform DMA-copyand-pack and remove the CPU-based packing step (at least for GEMM operation).

We believe future high performance systems will combine multiple computing architectures. Applications (and their sub-modules), upon their arithmetic intensity range, will be deployed and run on the most suitable or specific-designed architecture. Such a system will be heterogeneous and non-uniform in terms of memory technology, computing power and programming models. Abstraction layers and domain-specific frameworks are important to end-users who are not familiar to the hardware system. All the complexity is kept in the underlying level, where system developers will strive to deliver the highest performance.

The energy-efficient design of embedded many-core architectures is attractive to theoretical HPC power budget, but in reality may sacrifice ease of programming and sustainability of performance due to the reduced hardware feature set and memory constraints.

Opening minds to adopt new architectures, identifying their strength and weakness, designing portable libraries and programming models are steps to be taken to keep the HPC community moving forward. With this final message we conclude this thesis.