Erwan Lenormand
email: erwan.lenormand@cea.fr

Thierry Goubier
email: thierry.goubier@cea.fr

Loic Cudennec
email: loic.cudennec@intradef.gouv.fr

Henri-Pierre Charles
email: henri-pierre.charles@cea.fr

Loïc Cudennec

Data management

Keywords: Distributed Shared Memory, Field Programmable Gate Array, Irregular application

come

Introduction

As a response to the power wall problem, High Performance Computing (HPC) systems heterogeneity is gradually increasing. Part of this heterogeneity comes from the association of processors with co-processors. These latter allow the execution of computational-intensive portions of applications, called compute kernels, with high FLOP/W efficiency. Graphics Processing Units (GPUs) have become very popular for performing high arithmetic intensity processing on dense datasets. However, as illustrated by the efficiency on the High Performance Conjugate Gradient (HPCG) benchmark [START_REF]High-Performance Conjugate Gradient (HPCG) Benchmark results[END_REF], these systems achieve only a fraction of their theoretical peak performance and do not show efficiency gains from their heterogeneity for irregular applications. This poor performance is due on one hand, to the random data access patterns generated by these applications, and on the other hand, to the complexity of porting irregular compute kernels to GPUs. Data accesses cause saturation effects on memory subsystems and interconnection systems. Because GPUs are optimized to perform dense data vector computations, running irregular compute kernels reduces their efficiency. On the contrary, thanks to their reconfigurable architecture, Field-programmable gate arrays (FPGAs) are particularly suitable for processing irregular compute kernels [START_REF] Escobar | Suitability Analysis of FPGAs for Heterogeneous Platforms in HPC[END_REF]. The attractiveness of FPGAs for HPC systems is growing by means of their increasing computing power and the improvement of High Level Synthesis (HLS) tools. However, porting irregular compute kernels to FPGA remains a challenging task, because random data access patterns limit the abilities of HLS tools. Thus, designers must deal with low-level kernel design, optimization of data structures for FPGA memory systems and orchestration of distributed data transfers.

To address this issue, we propose a data management model for irregular compute kernels targeting heterogeneous distributed systems with reconfigurable accelerators.The latter is based on a shared-memory provided by a Software-Distributed Shared Memory (S-DSM). The application datasets are sliced in chunks managed by the S-DSM. The integration of reconfigurable accelerators in the S-DSM allows devices to initiate accesses to chunks. In this way, all the processing units can make fine-grained random data accesses. This unified data access model simplifies programming and meets the needs of irregular applications. By abstracting the data structure, chunk partitioning enables to prefetch the data as streams of chunks. This prefetching should make it possible to hide high data access latencies by implicitly overlaying the transferred data flow with the processed data flow. The efficiency of the proposed data management model relies on its ability to hide latencies. To assess this efficiency, we have used two case studies: General Sparse Matrix-Matrix Multiplication (SpGEMM) and a tsunami simulation code. These two applications generate a lot of irregular memory accesses, which are complex to optimize because they are data-dependent.

The paper is organized as follows: Section 2 presents the data management model, Section 3 describes the experiments conducted to validate the model, Section 4 gives some references on related work, finally, Section 5 concludes this paper.

Data management model

Shared memory is a convenient programming paradigm to develop multi-threaded applications, which randomly access data. Software-Distributed Shared Memory can be used to aggregate distributed physical memories into a shared logical space. In this work, we consider a S-DSM for heterogeneous micro-server that has been proposed in previous work [START_REF] Cudennec | Software-Distributed Shared Memory over heterogeneous microserver architecture[END_REF]. The latter is organized as a semi-structured super-peer network, where a set of clients are connected to a peer-to-peer network of servers. Clients execute the user code and servers manage the shared data and related metadata. The proposed data management model relies on the integration of reconfigurable accelerators in the S-DSM. The goal is to enable compute kernels to request access to distributed data. Obviously, this way to access the data can lead to high access latencies. To deal with this problem, the data management model aims to hide data access latencies by overlaying the transferred data flow with the processed data flow. This relies on the ability to access data as continuous streams. To do this we use chunks. The chunk is a common object in computer science, widely used in peer-to-peer programs. This concept is to use metadata to describe the data stored in it. Chunks are the atomic piece of data managed by the S-DSM. Each one has an unique identifier (chunk ID) and their maximum size can be set by the application. We use them to represent irregular data structures, as they are convenient objects for data management in distributed systems and their metadata allow to abstract the stored data. From the point of view of the compute kernels, the role of the S-DSM is to transparently provide the data and metadata corresponding to chunk ID. By partitioning the data structures according to the access granularity of the applications, data streams can be generated from sequences of chunk ID. Adapting the size of chunks to the granularity of access allows to avoid the transfer of unnecessary data. We have chosen two data structures widely used in irregular applications to illustrate the data management model: sparse matrix and unstructured meshes.

Sparse Linear Algebra (SLA) consists in performing linear operations on matrices (or vectors) for which the majority of the elements are equal to zero. Sparse matrices are compressed to reduce their memory footprint and to accelerate access to their nonzero elements (NNZ). The compressed sparse row format (CSR), shown in Figure 1b, is one of the most used sparse matrix representations. The column indices and the values of NNZ are stored in row-major order in the arrays Col and Val. RP [i] indicates the position of the first element of row i in the

0 0 A 0,2 0 0 A 1,1 0 A 1,3 0 0 0 0 A 3,0 0 0 0 (a) Dense format. Val A 0,2 A 1,1 A 1,3 A 3,0 Col 2 1 3 0 RP 0 1 3 3 4
(b) CSR format.

(2,A 0,2 ) (a) Hilbert space-filling curve over the mesh. 1c, we have adapted the CSR format to the use of chunks. We colocalize the value and the column index of an element to form a pair. The set of pairs representing a row is stored in a chunk. Then we use chunks metadata to indicate the number of elements contained in the row. This structure reduces the number of memory accesses required to read or write a matrix row. It can be easily adapted to another compressed format (e.g. compressed sparse column format). Reading or writing a matrix involves to request access to each row and to request the transfer of rows data between the memory and the compute kernel. Decoupling the access request and the transfer request allows the prefetching of the data into the FPGA memory and thus hides the access latency. Considering that the kernel is developed as a pipeline of stages, which are separated by FIFOs, then the prefetching speed is implicitly limited by the size of the FIFOs. This prevents the FPGA memory from being overloaded due to too early data prefetching. In an ideal case prefetching speed corresponds to the speed of data consumption of the following stages in the pipeline.

(1,A 1,1 ) (3,A 1,3 ) (0,A 3,0 ) ID count 0 1 3 1 2 1 (c) Chunk-based CSR format
Shallow Water Equations (SWE) is a Finite Element Method (FEM) applications consists of loops iterating over elements, nodes, or edges of meshes. Many HPC applications work on unstructured meshes with triangular elements in 2D and tetrahedral elements in 3D. Typical kernels on such unstructured meshes proceed by mesh updates -updating all elements of the mesh according to a function of each element neighbourhood values, applying a convolution or stencil and hence following indirections to both iterate over the mesh and to access neighbourhood information. Consequently, the topology of the mesh and indexing of data has a significant impact on data access locality and therefore application performance. Space-filling curves (SFC) allow to improve data access locality of the mesh [START_REF] Bader | Space-Filling Curves: An Introduction with Applications in Scientific Computing[END_REF]. We use this technique to do an efficient chunk partitioning of the mesh. As shown in Figure 2a, a SFC is drawn in the geometric space of the mesh. Vertices are indexed according the order in which they meet the curve. As illustrated in Figure 2b, we apply a basic partitioning along the curve, which consists in grouping the values of nodes of consecutive indices in chunks of constant size. The elements are numbered in order of the smallest index of their vertices. By following the path of the curve, most of the data of the mesh could be accessed through a sliding window, whose size would not be dependent of the mesh size. Thus, traversal of the mesh would be done through a continuous flow of data, where the majority of chunks would be accessed only once. Only the data of the elements located at the junction zones between the different spaces of the curve would not be accessible through the window. By following the curve, it is possible to identify the corresponding chunks. In this way, a chunk ID sequence corresponding to these data can be generated. We use these observations to design kernels iterating over unstructured meshes. The sliding window is implemented with an addressable FIFO. Buffers are used to access data not accessible through the sliding window.

Data management model validation

To validate the proposed data management model and assess its ability to hide data access latencies, we have conducted experiments with a simulation tool. This tool makes it possible to evaluate the performance of the system from high level modeling without requiring a full FPGA synthesis. The experiments focused on sparse matrix-matrix multiplication and a tsunami simulation code.

Simulation methodology

To conduct the experiments, we have chosen to use a simulation tool that we have developed [START_REF] Lenormand | A combined fast/cycle accurate simulation tool for reconfigurable accelerator evaluation: application to distributed data management[END_REF]. The objective was to evaluate the performance of the system from a high-level modeling. The irregularity of the applications we are studying and the distributed nature of the system we are targeting imply high and variable data access latencies. Thus, the main objective of this tool is to evaluate the effects of latency on the ability to speed up compute kernels using our data management model. Performance evaluation is based on the generation of data access latencies relating to the activity of the compute kernel. The tool uses a hybrid method: the activity of the compute kernel is generated by a simulation engine and latencies are produced by measuring the real latencies of S-DSM requests executed on the physical architecture, in order to produce faithful latencies. The simulation engine and the S-DSM server can be run on different nodes. This makes it possible to study different topologies associated with different latency profiles. In the rest of the section three topologies are used: No Latency which corresponds to the ideal case where all the data is stored in the FPGA memory, Local which corresponds to the case where the FPGA is located on the node running the S-DSM server and storing the data, and Remote which corresponds to the case where the FPGA is located on a node remote from the one running the S-DSM server. Local node latencies are medium (383 µs for a read request and 207 µs for a write request) and remote node latencies are high (1311 µs for read and 533 µs for write). We have used a Xilinx Virtex VC707 as a reference FPGA to set up the simulation engine. Thus, the clock frequency was set to 200 MHz and the theoretical peak memory bandwidth between the DDR and compute kernels was 12.8 Gb/s. The simulation being non-deterministic, the results presented are the median values of 10 runs.

Case study 1: General Sparse Matrix-Matrix Multiplication

Request A i & C i Read A i Request B k0,: A i,k0 Read B k0,: Compute A i,k0 × B k0,: PE0 Request B k1,: A i,k1 Read B k1,: Compute A i,k1 × B k1,: PE1 Request B k2,: A i,k2 Read B k2,: Compute A i,k2 × B k2,: PE2 Ci Cp i,: Cp i,:
Cp i,:

Write C i Fig. 3: Dataflow of the SpGEMM compute kernel using 3 PEs.

SpGEMM is widely used to study acceleration methods for sparse linear algebra. This application generates irregular memory access patterns that makes it complex to optimize, with usually a low efficiency in terms of floating point operations per unit of time. We have designed the compute kernel by using the row-wise sparse matrix-matrix multiplication algorithm formulated by Gustavson [START_REF] Gustavson | Two Fast Algorithms for Sparse Matrices: Multiplication and Permuted Transposition[END_REF]. Thanks to the row-wise traversal of the matrices, this algorithm is well suited to dataflow processing and is quite straightforward to parallelize. As illustrated in Figure 3, to parallelize the computations, the kernel is implemented with several processing elements (PEs). The first stages of the kernel access the nonzero elements of the first input matrix and distribute them to the PEs. Each PE multiplies the elements received by the corresponding rows of the second input matrix. Finally, the last stages sum the partial results computed by the PEs and write the result matrix. The indices and the values of the matrix are encoded with 4 bytes (Single precision computations).

As the arithmetic intensity of SpGEMM are strongly data-dependent, we have chosen matrices, presented in Table 1, with varying sizes, densities and patterns. Thus, the experiments allows to evaluate the capacity of the data management model model to adapt to irregularity. In order to limit the simulation time, the memory footprints of the matrices are smaller than a FPGA DRAM. To reproduce a situation where the capacity of the accelerator memory requires to transfer the data during the execution, we have adapted the simulated memory capacity accordingly to the dataset. Thus, the accelerator memory have been configured with 65536 locations of 1 kib (64 Mib). For each matrix, we have defined the theoretical peak computation speed by considering the processing time as the size of data transferred between the memory and the compute kernel divided by the FPGA memory bandwidth.

Table 1: Square matrices, from [START_REF] Davis | The University of Florida Sparse Matrix Collection[END_REF], used for simulations. NNZ and density refer to the source matrix. The memory footprint includes the three operand matrices. Fig. 5: Computation speed in GFLOP/s according the system topology (higher is better). The horizontal lines are the theoretical peak computation speeds.

For the first experiment, we varied the parallelism level of the compute kernel by implementing between 4 and 64 PEs on the local node. Figure 4a shows the computation speed obtained for this experiment. This shows that the increase in parallelism makes it possible to speed up computations, up to 16 PEs. The speed up obtained between 16 PEs and 32 PEs is low (between 1.04 and 1.22) or even negative. Between 32 and 64 PEs the speed ud is always negative. This efficiency limitation means that the PEs are under-exploited due to an insufficient supply of data (data starvation). The latter can be explained either by a data starvation in FPGA memory (due to excessive latencies), or by a FPGA memory bandwidth bottleneck. Figure 4b illustrates the occupancy rate of the FPGA memory controller. These results show that the controller is saturated for the configuration with 16 PEs. This information highlights that the FPGA memory bandwidth is the bottleneck for this kernel. This bandwidth limit is also one of the explanations for the nonlinear speed up between the configurations with 4 PEs and 16 PEs. The second experiment aimed to study the impact of topology on performance. For this we have used a configuration with 16 PEs, able to saturate the memory controller on the local node. The results obtained are illustrated in Figure 5. It shows that for the matrices consph, F2 and m t1 the performance gap with the ideal case for the local node (between 1% and 2%) and the remote node (between 2% and 6%) is very low. This small performance gap is mainly explained by the time to load the first data required to reach the nominal mode of the kernel. This shows the ability of the data management model to hide data access latencies. For the matrices cop20k A and s3dkt3m2, the performance gap is larger, but remains relatively small, respectively 38% and 6% for local node and 54% and 32% for remote node. For these two latter, the performance gap with the theoretical peak is also the largest, even for the ideal case. These results highlight a correlation between data density and the ability to speed up computations. Indeed, for the most sparse matrices, the memory accesses at row granularity do not use all the width of the data bus. Therefore, sparsity amplifies the effect of memory bandwidth bottleneck. Moreover, the more the rows are sparse, more processing time is short. This limits the ability to overlay the processing flow with the data transfer flow.

Case study 2: Shallow Water Equation

The Shallow Water Equations (SWE) are hyperbolic partial differential equations that describe a layer of fluid below a pressure surface. They can be solved with FEM. The code under study is the TsunAWI simulation code, a production code that implements the SWE with inundation, and whose results are used in the Indonesia Tsunami Early Warning System (InaTEWS), and under real-time constraints in the LEXIS European project [START_REF] Goubier | Real-Time Model of Computation over HPC/Cloud Orchestration -The LEXIS Approach[END_REF]. The base data structure is a 2D unstructured mesh. This code has been optimized for performance, especially concerning the mesh ordering [START_REF] Goubier | Fast Tsunami Simulations for a Real-Time Emergency Response Flow[END_REF].

In this code, we have designed a kernel to speed up the calculation of the gradient, one of the operations of the tsunami simulation code. This operation is an interpolation of the sea surface height via barycentric coordinates. The barycentric coordinates are precomputed for each vertex element. The processing of an element requires five floating point operations, resulting in a low FLOP/byte ratio. The kernel is implemented with several independent processing elements (PEs). Each PE process a different part of the mesh. The size of the sliding window is 1024 words. For the experiments, we have used five meshes presented in Table 2. The FPGA memory have been configured with 32 Mib. This was defined according to the maximum number of processing elements implemented (16 PEs), in order to allow each PE to prefetch up to 512 chunks of 1 kib per stream of data. Considering the memory footprint of each mesh and the maximum memory bandwidth of the FPGA, the theoretical peak computation speed is approximately 426 MElements/s for all meshes.

The first experiment aimed to study the speed up efficiency of the kernel according the parallelism level. To do this, the kernel has been implemented with 1 to 16 PEs. The computation speeds obtained are represented in Figure 6a. It shows that the speed up efficiency decreases sharply beyond 2 PEs. The occupancy rate of the memory controller illustrated in the Figure 6b time of the first data as a percentage of simulation time. This way the space between the colored bars and the gray bars represents the part of time during which the controller could have been used more. For the four largest meshes the results show a saturation of the memory controller for configurations with 4 PEs or more. We have concluded that the bandwidth of the FPGA memory limits the performance scaling of the kernel. From 2 PEs, the controller occupancy rate is too high to efficiently speed up the processing by increasing the parallelism.

For the smallest mesh, the loading time of the first data represents a significant part of the total processing time, which reduces the computation speed.

For the second experiment, we have evaluated the kernel performance according to the topology with 4 PEs. Figure 7a and 7b respectively illustrate the computation speed obtained and the associated occupancy rate of the memory controller. These results show that increasing data access latency reduces computation speed and highlight a correlation between the size of the mesh and the slowdown. As shown in Figure 7b, this effect can be explained by the proportion of the processing time spent to load the first data. For the large datasets processing on the local node, where the load time is the least impacting, processing speeds almost reach those of the ideal case. For the remote node, where the read access latency is three times greater than on the local node, the slowing down of the computation speed is relatively low for the three largest meshes (from 10% to 17%). Finally, the performance of the No Latency configuration is close to the theoretical peak computation speed. The gap is due to the inability to access all the node from the sliding window, which requiring to read several times some data. Thus, we have concluded that this mesh traversal method is almost ideal.

Discussion

The experiments have evaluated a data management model where accelerator tasks initiate access to distributed data. The experimental scenario used was the most disadvantageous, as the FPGA memory was empty at startup and all data had to be transferred during runtime. The results showed that thanks to prefetching, the programming model can efficiently hide the latencies of distributed data access. Nevertheless, this efficiency depends on the workload of the compute kernel. In practice, the observed workloads are huge. The size of the sparse matrices used in scientific applications can exceed ten gigabytes. The size of the complete datasets used for the tsunami simulation are at least ten times larger than the data subset used for the calculation of the gradient only. For the complete simulation each element of the mesh involves a hundred floating operation per iteration. Thus, the processing time of the largest meshes on a high end processor can exceed several hours. This processing time is much too long for urgent computing applications as tsunami simulation. What motivates to distribute the processing. This work shows that an S-DSM could simplify the distributed data management thanks to chunk partitioning and that the presented data management model solves the data access latency issue. Experiments have shown how in this model an FPGA can be supplied with data. As each accelerator is master of its access to data, this model can be extended to a distributed system integrating several FPGAs.

Related Work

Prior work have been done to provide shared memory for distributed systems with accelerators. Willendberg et al. [START_REF] Willenberg | A Remote Memory Access Infrastructure for Global Address Space Programming Models in FPGAs[END_REF] have proposed an FPGA communication infrastructure compatible to GASNet. This enables processing elements implemented on an FPGA to initiate remote direct memory access to remote FPGAs. Unicorn [START_REF] Beri | The Unicorn Runtime: Efficient Distributed Shared Memory Programming for Hybrid CPU-GPU Clusters[END_REF] provides a distributed shared memory (DSM) for CPU-GPU clusters. This is achieved with transactional semantics and deferred bulk data synchronization. StarPU [START_REF] Augonnet | StarPU: a unified platform for task scheduling on heterogeneous multicore architectures[END_REF] uses a DSM to manage data replication for heterogeneous distributed systems, but this DSM is not directly exposed to users. Recent work has studied chunk partitioning applied to sparse matrix for acceleration of sparse linear algebra. Winter et al. [START_REF] Winter | Adaptive Sparse Matrix-Matrix Multiplication on the GPU[END_REF] have proposed an adaptive chunk-based SpGEMM for GPU. This approach uses chunks to store the partial results of multiplication, then uses the chunk metadata for the merge stage. Rubensson and Rudberg [START_REF] Rubensson | Chunks and Tasks: A programming model for parallelization of dynamic algorithms[END_REF] have proposed the Chunks and Tasks programming model for parallelization of irregular applications. In this model, matrices are represented by sparse quatrees of chunks. MatRaptor [START_REF] Srivastava | MatRaptor: A Sparse-Sparse Matrix Multiplication Accelerator Based on Row-Wise Product[END_REF] and REAP [START_REF] Soltaniyeh | Synergistic CPU-FPGA Acceleration of Sparse Linear Algebra[END_REF] uses a chunk-based CSR format adaptation and the row-wise product to implement SpGEMM kernel on FPGA. Barrio et al. [4] have proposed an unstructured mesh sorting algorithm to enabling stream processing for finite element method applications. This algorithm was applied to study the acceleration of scientific codes on CPU-FPGA platform.

Conclusion

Increasing the energy efficiency of HPC systems has become a major issue. Thanks to their reconfigurable architecture, FPGAs could increase power efficiency for HPC applications with irregular compute kernels. However, due to their complexity of use, FPGAs are underemployed in HPC systems. In this paper we have proposed a data management model for irregular compute kernel acceleration on FPGA integrated in distributed system. This model relies on a S-DSM to allow accelerators to initiate access to distributed data and on chunk partitioning to abstract the irregular structure of the datasets. We have shown how this data management model could be applied to compute kernels of sparse linear algebra and finite element method. We have conducted experiments with a hybrid simulation tool, which exploits the physical system to provide accurate data. These experiments have shown that the data management model enables to efficiently hide high data access latencies. Finally, experiments have shown that memory bandwidth is a bottleneck. This phenomenon is normal since the studied applications are memory bound. High Memory Bandwidth (HBM) technologies as available on current and future FPGAs should help to remove this bottleneck and improving performance of compute kernels.

Fig. 1 :

 1 Fig. 1: Matrix representation.

  Chunk partitioning of the mesh with chunks of two elements.

Fig. 2 :

 2 Fig. 2: Reordering and partitioning of a 2D unstructured mesh.

Fig. 4 :

 4 Fig. 4: Performance according to the number of processing elements (PEs).

  Computation speed in MElements/s (higher is better). P a d a n g _ C P a d a n g _ F C o q u im b o _ C C o q u im b o _ F M e d it e r r a n e Memory controller activity (occupancy percentage). Gray bars correspond to the loading time of the first data.

Fig. 6 :

 6 Fig. 6: Performance according to the number of processing elements (PEs).

  provides a better understanding of these results. In this figure, the colored bars represent the percentage of time the controller is active. The gray bars correspond to the loading P a d a n g _ C P a d a n g _ F C o q u im b o _ C C o q u im b o _ F M e d it e r r a n e Memory controller activity (occupancy percentage). Gray bars correspond to the loading time of the first data.

Fig. 7 :

 7 Fig. 7: Performance according the system topology.

Table 2 :

 2 Characteristics of the set of meshes used for the experiments.

	Region	Name	# Elements # Vertices Memory footprint
	Indian Ocean	Padang C Padang F	460 119 2 470 345 1 242 653 231 586	14 Mb 74 Mb
	Pacific Ocean	Coquibo C Coquimbo F	3 396 755 1 709 506 9 762 027 4 887 927	102 Mb 293 Mb
	Mediterranean Sea Mediterranean	9 917 645 4 999 404

This work was supported by the LEXIS project, funded by the EU's Horizon 2020 research and innovation programme (2014-2020) under grant agreement no. 825532.