Louis Ledoux
email: louis.ledoux@bsc.es

Marc Casas
email: marc.casas@bsc.es

A Generator of Numerically-Tailored and High-Throughput Accelerators for Batched GEMMs

We propose a hardware generator of GEMM accelerators. Our generator produces vendor-agnostic HDL describing highly customizable systolic arrays guided by accuracy and energy efficiency goals. The generated arrays have three main novel aspects. First, the accelerators handle a large variety of computer number formats using intermediate representations based on our Sign Scale Significand (S3) format. Second, the processing elements perform all intermediate dot-product arithmetic operations required by the GEMM kernel without any intermediate rounding, which makes it possible to deliver better energy efficiency than state-of-the-art approaches while offering more accuracy and reproducible results. Third, our accelerators feature the Half-Speed Sink Down (HSSD) mechanism, which maximizes the overlap of host-accelerator data transfers with GEMM computations.

We evaluate our automatically generated designs in a cuttingedge setup composed of a POWER9 host, CAPI (Coherent Accelerator Processor Interface) link, and a Virtex Ultrascale Plus FPGA. Arrays can operate at the speed of the link and saturate it to reach a 13GB/s throughput. Our fine-grain customization approach allows to cover a wide range of accuracy versus efficiency scenarios and can reach 0.65GOps/s/W while producing 1024 accurate bits or 148.7GOps/s/W with 6 accurate bits. Our configurations achieve up to 1613GOps/s system performance and power efficiencies of up to 240GOps/s/W for the FPGA. This automatic generator is the first being able to produce such a variety of designs. We improve the single-precision energy efficiency of state-of-the-art FPGA GEMM accelerators by 1.86×.

I. INTRODUCTION

A large variety of scientific computing workloads rely on the GEneral Matrix-Matrix (GEMM) multiplication kernel. The requirements of these workloads in terms of numerical accuracy or performance are very heterogeneous. Indeed, emerging workloads such as Deep Neural Networks are naturally resilient to arithmetic tinkering [START_REF] Johnson | Rethinking floating point for deep learning[END_REF], even though they embed millions of parameters [START_REF] Szegedy | Going Deeper with Convolutions[END_REF] [START_REF] He | Deep Residual Learning for Image Recognition[END_REF]. Conversely, ill-conditioned linear systems require much higher precision arithmetic [START_REF] Beliakov | A Parallel Algorithm for Calculation of Large Determinants with High Accuracy for GPUs and MPI clusters[END_REF] [68] [START_REF] Ellis | One-loop amplitudes for W+3 jet production in hadron collisions[END_REF] [26] [START_REF] Lake | From Sir Isaac to the Sloan Survey Calculating the Structure and Chaos Owing to Gravity in the Universe[END_REF] [10] [START_REF] Iakymchuk | Reproducible and Accurate Matrix Multiplication for GPU Accelerators[END_REF]. Despite this large variety of accuracy requirements, commodity hardware just supports integers or 32-and 64-bits IEEE-754 arithmetic [START_REF]Ieee standard for binary floating-point arithmetic[END_REF], which often results in a waste of computing resources.

To address these issues, computer number formats alternative to the IEEE-754 standard, like the 16-bit Brainfloat (Bfloat16) [START_REF] Abadi | TensorFlow: Large-Scale Machine Learning on Heterogeneous Distributed Systems[END_REF] [START_REF] Kalamkar | A Study of BFLOAT16 for Deep Learning Training[END_REF], the Tapered Floating-Point (TFP) [START_REF] Morris | Tapered floating point: A new floating-point representation[END_REF] or the Posit [START_REF] Gustafson | Beating Floating Point at its Own Game: Posit Arithmetic[END_REF], are attracting the attention of computer architects. There is a variety of previous works implementing arithmetic-logic units (ALUs) [START_REF] Jain | CLARINET: A RISC-V Based Framework for Posit Arithmetic Empiricism[END_REF] [4] [START_REF] Tiwari | PERI: A Configurable Posit Enabled RISC-V Core[END_REF] and bare arithmetic data paths [START_REF] Chaurasiya | Parameterized Posit Arithmetic Hardware Generator[END_REF] [67] [START_REF] Uguen | Evaluating the Hardware Cost of the Posit Number System[END_REF] for alternative formats, and comparing their area requirements to IEEE-754 units. In addition, some studies have analyzed the numerical stability of the Posit and the Bfloat16 arithmetics with respect to IEEE-754 [START_REF] De Dinechin | Posits: the good, the bad and the ugly[END_REF] [START_REF] Buoncristiani | Evaluating the Numerical Stability of Posit Arithmetic[END_REF]. Despite this extensive amount of work, there are no approaches considering the area, performance, power, and accuracy tradeoffs of these alternative computer number formats in the context of a real system implementing them. Also, there is no unified proposal leveraging the large flexibility of alternative number formats with respect to the bitwidth or the accumulator size, which covers an ample scope in terms of performance, accuracy, area, or power.

To handle these issues, this paper makes the following contributions beyond the state-of-the-art: • The automated generation of arithmetic datapaths composed of Sign Scale Significand (S3), a hardware-friendly floatingpoint format, in conjunction with extended or reduced internal precision accumulators. We propose a systolic array generator that exploits these automatically generated arithmetic datapaths to support a large variety of computer formats with customizable accumulators while producing reproducible results.

• The elimination of intermediate rounding, which provides large benefits in terms of energy efficiency, reproducibility, and accuracy.

• The Half-Speed Sink Down (HSSD) mechanism, a stallingfree output scheme that allows generated designs to maximize the use of processing elements and bandwidth.

• It provides a large evaluation campaign of our generated accelerators considering 64 different designs combining four different bitwidths, three accumulator sizes, and four computer number formats. We deploy the generated hardware on a highend FPGA (Virtex Ultrascale Plus) attached to a POWER9 system via IBM CAPI link. This campaign exhibits the wide variety of scenarios covered by the generated designs. When combining small bitwidths with exact accumulators, our designs deliver 240GOps/s/W and display relative numerical errors below 2 -4 . To the contrary, high-precision designs can deliver 1024 bits of accuracy and 0.65GOps/s/W . Section V shows that our 32-bit designs deliver 19.5GOps/s/W , which outperforms the FPGA state-of-the-art.

978-1-6654-8332-2/22/$31.00 ©2022 IEEE

II. RELATED WORK A. Matrix Multiplication Systolic Arrays

Early systolic arrays originated during the 1980s [START_REF] Kung | Systolic Arrays for (VLSI)[END_REF] and were mostly used for convolutions [START_REF] Kung | A Systolic 2-1) Convolution Chip[END_REF], linear systems [START_REF] Kung | Why systolic architectures?[END_REF] and matrix multiplications [START_REF] Quinton | A new matrix multiplication systolic array[END_REF]. Because it is not trivial to map a time-sequential algorithm into an ad-hoc space-time hardwired algorithm, previous works propose to formalize this process [START_REF] Kung | Optimal Systolic Design for the Transitive Closure and the Shortest Path Problems[END_REF] [5] [START_REF] Navarro | Computing sizeindependent matrix problems on systolic array processors[END_REF]. Systolic arrays' prominence and popularity rekindle due to the massive demands of computations brought by recent AI algorithms. In fact, commercialized ASICs for systolic arrays comprise cloud-and edge-deployed tapeouts with different array sizes. Notable examples are the Google's TPUs [START_REF]System Architecture | Cloud TPU[END_REF] for the cloud and Tesla FSD [START_REF] Lamert | Accelerated Mathematical Engine Tesla[END_REF] for the edge.

Gemmini [START_REF] Genc | Gemmini: An Agile Systolic Array Generator Enabling Systematic Evaluations of Deep-Learning Architectures[END_REF] is a systolic array generator targeting ASIC with a RISC-V host for edge devices. Gemmini generator has been utilized for the tapeouts of two System-On-Chips (SoC) in TSMC 16nm and Intel22FFL process technologies. The fact that this generator is written in Chisel, makes it data-type generic at the cost of suboptimal datapaths.

B. Arithmetic-level Optimizations

The IEEE-754 [START_REF]Ieee standard for binary floating-point arithmetic[END_REF] 1985 standard describes various floatingpoint formats with different bitwidths: half-, single-, and double-precision numbers, which occupy 16, 32, 64 bits of storage, respectively. This small set of general-purpose datatypes can not optimally handle the plethora of specific numerical requirements proposed by modern workloads. This is the reason why scientists propose emerging arithmetics, ALUs, and even mathematical frameworks to describe real numbers [START_REF] Lindstrom | Universal coding of the reals: Alternatives to ieee floating point[END_REF].

Bfloat16 [START_REF] Wang | BFloat16: The secret to high performance on Cloud TPUs[END_REF] is 16-bit wide and built from an IEEE-754 32bit truncated from its LSBs and therefore shares identical dynamic range as its 32-bit counterpart. It is becoming a popular format with software [START_REF] Wang | Distributed MLPerf ResNet50 Training on Intel Xeon Architectures with TensorFlow[END_REF] and hardware (FPGA [START_REF] Chromczak | Architectural Enhancements in Intel Agilex FPGAs[END_REF] [19] [START_REF] Petrov | Using of Bfloat16 Format in Deep Learning Embedded Accelerators based on FPGA with Limited Quantity of Dedicated Multipliers[END_REF], TPU [START_REF] Wang | Benchmarking TPU, GPU, and CPU Platforms for Deep Learning[END_REF], ARMv8 ISA [START_REF] Stephens | BFloat16 extensions for Armv8-A[END_REF]) support.

The Tapered Floating-Point (TFP) [START_REF] Morris | Tapered floating point: A new floating-point representation[END_REF] format uses variablesized entries to represent the mantissa and scale instead of the fixed-length entries of IEEE-754. We define TFPs ωS, ωF as numbers represented in scientific notations with a scale on ωS bits, a fraction on ωF bit, and a sign bit. We distinguish TFP from the IEEE-754 floating-point by its absence of denormal numbers.

The Posit number system [START_REF] Gustafson | Beating Floating Point at its Own Game: Posit Arithmetic[END_REF] is a particular case of tapered formats that improves entropy per bit by trading exponent for mantissa storage bits via a Golomb-Rice [START_REF] Golomb | Run-length encodings (Corresp.)[END_REF] encoding. One posit configuration is usually refered as "posit N, es ", where N is the bitwidth, and es the variable length exponent size.

III. SYSTOLIC ARRAY GENERATOR

Our hardware generator of GEMM accelerators relies on systolic arrays. The GEMM accelerators can target any FPGA or ASIC because we describe RTL logic neither correlated with timings nor hardware target. In fact, we leverage flopoco [START_REF] De Dinechin | Designing custom arithmetic data paths with FloPoCo[END_REF] [START_REF] Istoan | Automating the pipeline of arithmetic datapaths[END_REF] for the back-end translation, which internally builds up timing graphs and inserts necessary and sufficient clocked elements.

A. Architecture of generated Accelerators

We define the generator parameters in a similar way as BLAS level-3 routines. Parameters n, m, and p correspond to the dimensions of the matrices involved in the GEMM operation, which we denote as A, B, and C. n, m, p ∈ N + , α, β ∈ R : C n,m ← αA n,p B p,m + βC n,m Fig. 1 illustrates the overall architecture of a 3 • 3 array. Input data and control come from the bus on the left. Our architectures align data with clocked registers and translate them using A2S3 units before data enter the actual array. Section III-C describes in detail the A2S3 units. Once in the array, every translated numerical value is expressed in the S3 hardware format, which is described in Section III-B2, and travels through the corresponding Processing Elements (PEs). Section III-C describes the PEs design. Start/End Of Block (SOB and EOB) are 2 bits generated by the main memory/host that control the array. SOB and EOB indicate when to flush to zero accumulated values and when accumulated values are valid to be carried out, respectively. This local and distributed control of the array allows efficient scaling. The last step is to round and normalize the output values via the normalization modules L2A to the output format. Section III-C describes L2A modules. Output and input formats can be different. To extract the in-place accumulated value from the outputstationary arrays, several options are possible with their respective trade-offs and requirements. Previous work [START_REF] Khan | [END_REF] exposes two of them: having global lines (buses) from PEs to output buses to form semi-broadcast arrays; or reusing PEs internal registers while stalling the incoming data. None of the induced constraints are suitable for our needs since stalling the incoming pipeline would seriously harm computation and communication overlap and deliver poor performance. Furthermore, given that we aim to produce generic-sized and fast systolic arrays, we can not accept the unpredictable delays caused by high fan-outs yield by global lines, and therefore, we do not consider them.

To remedy this problem, our generator leverages a stallingfree output scheme for the output matrix C that allows extracting final and intermediate results through a propagating Since our architecture is output-stationary, it does not need to preload any input matrix. Therefore, combined with the HSSD chain that allows the output matrix to exit simultaneously, PEs are constantly computing and are fully used after the initial traversal of the array. The HSSD chain allows all links between PEs to be local, thus inducing no global routes except for the clock and the reset, which totally matches FPGAs' building blocks. We describe our HSSD design in Section III-C2 at the PE scale level to clearly show how output data is transferred out of the systolic array. Our HSSD design makes it possible for the systolic arrays to perpetually pump data from the input borders while simultaneously feeding the output border where the normalizing and rounding units reside, as we show in terms of dotted lines in Fig. 1. Moreover, due to the non-stalling pipeline introduced by HSSD, all transfers can be pipelined and overlapped to hide the latency of the systolic array, which directly translates into increased throughput. As a consequence of the HSSD extraction capability, matrix blocks can be sent as batches without interruptions to perform many independent GEMM operations.

The generated designs fuse dot-product arithmetic operations so that no intermediate rounding occurs. The reason to postpone the rounding step and shift it to the bottom of the array is twofold: On the one hand, it allows removing the expensive rounding logic from the n • m PEs since it is just required by the m L2A modules (see Section III-C4). On the other hand, deferred rounding retains the full accuracy of accumulated values by avoiding cancellations when scales are drastically different. Section III-C3 describes in detail this deferred rounding mechanism.

B. Generator Configurability

The generator automatically produces GEMM accelerators featuring a wide range of options in terms of accuracy, area, performance, and energy consumption. These are the main reconfiguration parameters:

1) Array Dimensions: The height (h) and the width (w) of generated arrays define the number of PEs (number P Es = h • w). The more the h • w product increases, the more operations per second (Ops/s) the array can deliver since every PE performs a fused-dot-product at every clock edge.

2) Arithmetic: There are three aspects of our generator that define the arithmetic of the accelerators:

Sign Scale Significant (S3) format: The S3 format is a hardware-friendly floating-point format that maps well to FPGA internal basic elements. S3 is composed of the following fixed-size fields that code the fixed and variable size of input numbers:

• Sign on 1 bit: active high when the number is negative.

• Scale on ωS bits: unsigned scale with a bias. The bias is calculated differently depending on the arithmetic. • Significand on 1 + ωF bits: an implicit bit prepended to the unsigned significand of the number. The 0 value is coded with implicit and significand set to zeros. • NaN on 1 bit: Active high when the number is not a real number or when a problem occurs, such as overflows.

The fields form the N aN, sign, scale, implicit, f raction quintuple whose elements are concatenated in one bus of size ωS +ωF +3 bits. S3 can handle any computer number by just adjusting the N aN, sign, scale, implicit, f raction values. Table I provides examples illustrating how computer number formats can be expressed in terms of the S3 encoding.

Accumulator size: PEs can be coupled with variable-size accumulators, from exact to light configurations. The following parameters, which are inspired by previous work [START_REF] De Dinechin | An FPGA-specific approach to floating-point accumulation and sum-of-products[END_REF], define the accuracy and the size of the numerical fixed-point accumulators that our designs contain: • LSB: Is the weight of the Least Significant Bit of the accumulator. Adding bits to the right of the accumulator adds precision and controls the final accuracy. However, such bits are expensive because they augment the size of internal arithmetic modules (especially the barrel shifter, see Fig. 3).

• M SB: Is the expected weight of the summand. Similarly, lowering M SB can save hardware resources by decreasing the barrel shifter size.

•

OV F : Is the number of carry bits (see Fig. 3 OVF) prepended to the left as MSBs in order to prevent overflows due to accumulation. The accumulator bitwidth is ωLA = OV F + M SB -LSB + 1.

Normalizing units: Located at the bottom of the arrays, they translate the internal accumulations expressed in fixed-point format to the output computer number format. We propose three approaches to generate the rounding units: • Same: A standard approach is to perform rounding and normalizing of the output values into the format they entered the systolic array.

• Specific Arithmetic: By taking advantage of the hardware reconfigurability, our generator can normalize the two's complement accumulations into a different arithmetic format than the input one. This has many applications, such as promoting to a bigger data type to avoid the final roundoff error, sending to another hardware that supports a specific type only, or even demoting to a lower data type to save bandwidth.

• Exact: This option does not perform the final normalization step and, does not require the L2A modules. The accelerators implementing this option send the final data as fixed-point.

C. Modules designs

This section describes the Systolic Array (SA) and Systolic Array Kernel (SAK) components (see Fig. 1).

1) A2S3: Denormalizing units are called A2S3 and are the modules responsible for the conversion from an incoming arithmetic word to a S3 quintuple. They take an N -bit bus and output a ωS + ωF + 3 bits bus. Any computer format 2) Processing Element (PE): Our generator hardwired algorithm consists in four operations that PEs execute: First, register the incoming S3s values (S3_A and S3_B in Fig. 2) arriving respectively from the left and the right, to present them to their corresponding opposite border (right for S3_A_Q and bottom for S3_B_Q). Second, register the control signals and expose them to the opposite border (SOB and EOB from top to bottom). This allows a distributed and local control of the array and avoids global lines, which are often a problem. Third, compute in augmented precision without rounding the dot-product of incoming S3s values in an output-stationary fashion. This is done by a dedicated module discussed in Section III-C3. Fourth, extract via the HSSD mechanism the output-stationary values of C vertically with two registers C1 and C2 and a multiplexer, as Fig. 2 shows.

3) S3 Fused Dot Product (S3FDP): S3FDP is the kernel of our generator as it is the module that performs the fused-dotproduct between values encoded with the S3 format. Fig. 3 depicts the three main stages of the S3 fused-dot-product. Firstly (Fig. 3 1 ), the exact product is performed. Secondly (Fig. 3 2 ), the mantissa of the product is aligned with regard to the accumulator weights depending on the scale of the product. Finally (Fig. 3 3 ), the fixed-point addition is performed. The last operation requires peculiar timinganalysis. In fact, as fixed-point accumulators require as many bits as the dynamic range of the product contains binades, the carry-rippling can quickly become the critical path. If such, our generator automatically generates Carry-Save-Adders (CSAs) with flopoco timing analysis to suggest the chunk-size.

4) Normalizing modules (L2A): Their role is to translate the fixed-point values arriving from the HSSD chains to the specified output format. As mentioned in Section III-B2, the output format can vary and even be the fixed-point output coming from the S3FDP itself. The logic to normalize and round comprises the following main steps: First, to count Leading Zero/One (sign) to compute the scale of the accumulator. Second, to shift the accumulator by its scale to get the significand and a pre-sticky bit. Third, to build exponent fields based on the computed scale and the accumulator sign. The second and third steps are merged thanks to the LZOCStickyShifter operator from flopoco. And, fourth, to build the final word taking into account possible exceptions and applying a rounding scheme.

IV. DESIGN SPACE EXPLORATION

We carry out a design space exploration of the designs on a Virtex Ultrascale Plus (VUP) speedgrade-2 FPGA without loss of generality. Although the specific results of this design exploration apply to the VUP FPGA, the whole methodology can be applied to a generic FPGA. In fact, modern FPGAs share common properties and the utilization ratio between different building blocks is analogous.

A. Evaluated Designs

We consider four computer number formats: IEEE-754, Bfloat16, TFP, Posit. All of them are described in Section II. We consider four different bitwidths (8, 16, 32, and 64) for the IEEE-754 and the TFP formats, that is, 4 different variations for each format. For the case of Bfloat16 we consider the 16-bit width, since it is the only one it admits by definition. For the case of Posit, we consider the following 13 different formats: the posit 4, 0 4-bit design, the posit 8, 0 , posit 8, 1 , posit 8, 2 8-bit designs; the Additionally, we propose three different accumulator configurations. We define them in terms of the parameters we describe in Section III-B2. We propose the following accumulator designs:

AI(α): The α accumulator targets AI workloads. Since previous work [START_REF] Tambe | AdaptivFloat: A Floating-point based Data Type for Resilient Deep Learning Inference[END_REF] shows that neural networks weights values distribution are relatively narrow (e.g. within the [-15, 21] interval for Resnet50 and transformer models) and resilient to arithmetic tinkering [START_REF] Johnson | Rethinking floating point for deep learning[END_REF], we set the msb parameter to be 5. We set the lsb parameter to build a small accumulator of factor-2 linear size with regard to the arithmetic bitwidth. We have the size ωLA = 2 • N , and

msb = 5, ovf = 2, lsb = 2 • N -1 -msb -ovf EXACT(β):
The β accumulator size is dictated by the size of the S3 fields to ensure that any possible result of the exact multiplication can be represented in the accumulator without accuracy loss. The Kulisch [START_REF] Kulisch | The arithmetic of the digital computer: A new approach[END_REF] and Quire [START_REF] Mercado | mightymercado/PySigmoid[END_REF] accumulators land in this category.

CONSTANT_100(γ): The γ accumulator configuration correlates neither with the arithmetic internal fields nor with any workload numerical specificity. The size ωLA = 1+9(ovf )+40(msb)+50(-lsb) allows all arithmetic to have the same internal accumulator, which enables comparisons to be made on other aspects besides the accumulator size.

We consider the three accumulators for each one of the 22 computer number formats we evaluate, except for the case of posit 4, 0 where we just consider the β accumulator since α and γ do not produce meaningful accuracy for this format. We have thus a total of 64 different combinations. 4 how the resource numbers grow vis-à-vis the arithmetic and the accumulator configuration. These numbers are obtained after place and route with "Explore" directives through the generation of all design checkpoints. We can see how the bitwidth of the computer number format largely impacts the resource requirements of S3FDP units. Table II reports these numbers for the IEEE-754 format and bitwidths 8, 16, 32, and 64. For instance, the smallest IEEE-754 β does not use DSPs as all arithmetic datapaths are emulated with LUTs. To the contrary, the 53 bits multiplication and the 4228 bits addition of IEEE-754 64β require 6 DSPS, 272 CARRY8s with more than 10k LUTs. We also compare to α-, γ-S3FDPs and accumulators built on conventional FMA blocks infered thanks to Xilinx IPs. Our configurations always save resources against traditional FMA hardware, and thus, save power consumption. The only exception is the IEEE-754 64β which requires to maintain more than 4000 binades in fixed-point format.

V. EVALUATION A. Experimental Setup

All experiments run on an AC922 POWER9 system with attached FPGAs as accelerators. The POWER9 system has 2 sockets (chips), each with 20 cores and 4 threads/core for a total of 160 threads running at 2.3GHz. The FPGA is a xcvu3p- The physical link between the host and the FPGA is PCIe gen4-8lanes, whose 16GT/s provide up to 15.754 GB/s [START_REF]PCI Express[END_REF]. All communications are done on top of the IBM Coherent Accelerator Processor Interface 2(CAPI2 [START_REF] Morgan | Opening Up The Server Bus For Coherent Acceleration[END_REF] [47] [START_REF] Wikipedia | Coherent Accelerator Processor Interface[END_REF]) protocol with the IBM SNAP framework [START_REF] Ibm | CAPI SNAP Framework Hardware and Software. Contribute to open-power/snap development by creating an account on GitHub[END_REF]. To avoid Clock Domain Crossing (CDC), we generate all the arrays for an operating frequency of 250M Hz, which is the SNAP buses clock.

B. Performance, Energy Efficiency, Accuracy, and Throughput

We evaluate the throughput, performance, energy efficiency, and accuracy of generated designs by considering the 64 configurations we describe in Section IV-A running at 250M Hz. The number of PEs varies depending on the bitwidths: 64 • 63 for 4-bits, 32 • 31 for 8-bits, 16 • 15 for 16-bits, 8 • 7 for 32bits, and 4 • 3 for 64-bit widths. For instance, Fig. 5 shows placed and routed arrays of 32•31 and 64•63 PEs with the additional communication logic, which represents 30% of the total floorplaning.

All metrics we consider depend on the matrix dimensions (see Eq. III-A) and the number of matrix batches (b). We propose a Payload Size metric, which combines both matrix dimensions and the number of matrix batches. This metric is defined in terms of Bytes as follows:

P ayload size = N • n • m • p • b • 2/8.
The Payload Size is equivalent to the DMA size. For throughput, performance, and energy efficiency, the number of accumulations per batch does not impact the results. However, accuracy is strongly dependant on the number of numerical accumulations per batch, as we discuss in Section V-B4. We display results setting b = 1, except if we precise the contrary. All displayed values result from an average of 10 runs. We generate the coefficients of our matrices by using a uniformly distributed random variable within the interval [-1, 1].

1) Physical Link Throughput: Fig. 6 shows the average throughput achieved by the 64 arrays and the peak bandwidth offered by the 8-lanes PCIe gen4 link. The x-axis shows the payload size of the GEMM problem in terms of Bytes, and the y-axis shows the measured throughput in GB/s. All the measured throughput are averaged and plotted as one curve as the standard deviation is never more than 5% of the expected saturation throughput. The maximum measured throughput reaches ∼ 83% of the peak link bandwidth, which corresponds to payload sizes of 2 24 Bytes or larger. As previous work Fig. 7: Measured and averaged vs theoretical Performance (GOps/s) for different array sizes and Payload Sizes (Bytes). describes [START_REF] Lawley | Understanding Performance of PCI Express Systems[END_REF], the hardware PCIe protocol, which contains several layers and encoding schemes, makes it impossible to achieve the peak link bandwidth. In fact, our designs saturate data link bandwidth at their disposal. Fig. 8 compares throughputs at application level for different array sizes with and without our HSSD extraction chain. The experiment uses a constant payload size of 2 30 Bytes with a decreasing number of batches b, and therefore, an increasing size of the common dimension p (shown by x-axis of Fig. 8). Without HSSD, and thus, without overlaping incoming and outgoing DMAs, an array stalling of n clock edges happens every p clock edges to be able to read out the data. This experiment exhibits that stalling the pipeline for small matrices hinders performances up to two order of magnitude.

2) Performance: Fig. 7 shows the performance achieved by the 64 designs in terms of GOps/s. The x-axis shows data payload in Bytes, and the y-axis displays performance in terms of GOps/s. Our experiments clearly show how the designs experiment 5 possible behaviors in terms of performance depending on the bitwidth. The reason why bitwidths define the performance behavior of the designs is the fact that the number of PEs in the array is defined by the data format bitwidth, as we explain in Section V-B. The designs achieve a performance peak of 4.8, 22.5, 96.7, 399.7, and 1613 GOps/s for 64-, 32-, 16-, 8-and 4-bit numerical formats, respectively. We observe these peaks for payload sizes equal to or larger than 2 24 Bytes. They are very close the maximum theoretical values, which are 6, 28, 120, 496, and 2016 Gops/s for 64-, 32-, 16-, 8-, and 4-bit numerical formats, respectively, which correspond to arrays with 4 • 3, 8 • 7, 16 • 15, 32 • 31, and 64 • 63 PEs. The theoretical peaks are not achieved due to the physical link bandwidth saturation that we observe for payload sizes equal or larger than 2 24 Bytes in Fig. 6.

3) Energy Efficiency: We report the rate of operations delivered by the system per watt of power consumed in Fig. 9. We also compare to energy efficiencies delivered by equivalent arrays but with intermediate roundings, thus, mimicking the behavior of traditional FMA hardware. The than 2 28 Bytes. The corresponding β-designs, posit 8, 0 β, IEEE-754 8β, and TFP8β, deliver maximum energy efficiencies of 98.3, 84.9, 98.2 GOps/s/W, respectively. The difference between the power efficiencies delivered by these two groups arises from the size of the accumulator, which significantly impacts the power cost of arithmetic operations. This trend is confirmed for other arithmetic bitwidths. For instance, all β accumulators of 32-and 64-bits arithmetics are less efficient than α and γ accumulators due to the quadratic scaling of internal operators with regard to the arithmetic bitwidth. The most efficient 16-, 32-, and 64-bit designs are TFP16α, TFP32α, and TFP64γ, which deliver maximum efficiencies of 47.9, 19.5, and 5.3 GOps/s/W, respectively. In terms of comparing the energy efficiency of the four formats we consider, we can see how TFP and IEEE-754 formats deliver very similar efficiency when using the same bitwidth, while posit can be either better or worse than TFP and IEEE-754 for the same bitwidth depending on the es parameter, which is defined in Section II. The energy efficiency of designs with PEs including rounding units is much worse than the ones without intermediate rounding. For example, IEEE-32β with rounding units reaches energy efficiencies up to 14.5GOps/s/W , while its non intermediate rounding counterpart reaches 19.5GOps/s/W . In fact, considering all the de-signs, we observe an average loss of 70% of energy efficiency when the rounding unit is present in the accumulation step. 4) Accuracy: We report in Fig. 9 the final accuracy for different payload sizes. Workloads are composed of a single batch, so the number of numerical accumulations per PE is the same as the payload size. We compute the log2 of the relative error, which indicates the number of exact significant bits with respect to a high-precision oracle obtained with the mpfr library [START_REF] Fousse | MPFR: A multiple-precision binary floating-point library with correct rounding[END_REF]. The oracle consists in the same computation made in software with a fixed-point of 1024 significant bits as fractional part, and 1000 bits to code the integer part.

We observe that the smallest accumulators (all α and posit 4, 0 β) fail to maintain their initial accuracies while the number of accumulations increases. They are more statistically sensitive to overflows with respect to their larger counterparts, which can result in a N aN propagation (lines falling down in Fig. 9). All γ accumulators yield 50 accurate bits as it is equal to -lsb. These experiments show that the final rounding significantly impacts the accuracy as the internal precision can be lost to fit in a mantissa of general purpose size. We also observe how large accumulators are able to achieve good accuracy levels as accumulation counts increase. Some interesting 8 bits combinations are TFP8β and IEEE-754 8β, which sustain 10 and 15 bits of accuracy for all the considered accumulation counts and, according to Fig. 9, they deliver 98.2 and 84.9 GOps/s/Watt, respectively. More drastic combinations are posit 4, 0 β and TFP 64β. The TFP 64β always maintains the correct 1024 bits after millions of accumulated values while delivering 0.65GOps/s/W . Conversely, posit 4, 0 β delivers 240GOps/s/W but tends to overflow after 2 14 fused accumulations, which makes it lose its initial 4 accurate bits. The accurcy of designs with PEs including rounding units is much worse than the ones without intermediate rounding units. The most extreme cases are IEEE 64β, and any 64-bitγ, where versions with intermediate rounding lose 971 and 0 accurate bits, respectively, with respect to their counterparts without intermediate rounding.

C. Comparison with the State-of-the-art

Table III displays the most relevant characteristics of previous approaches implementing hardware accelerators for GEMM. The numbers are extracted from the respective manuscripts. Table columns show the year of publication, FPGA device, frequency, FPGA resource utilization, host system, physical link, performance/energy-efficiency they achieve for 64-and 32-bits numerical formats, the maximum number of accurate bits, the number of computer arithmetics they can use, and other aspects like numerical reproducibility, automatic generation, or vendor agnosticism. The column "reproducibility" specifies whether accumulation order hinders vector-dot-product results.

Moss et al. [START_REF] Moss | A Customizable Matrix Multiplication Framework for the Intel HARPv2 Xeon+FPGA Platform: A Deep Learning Case Study[END_REF] exploit the flexibility offered by the HPC hybrid Xeon/Arria 10 HARPv2 to offer a customizable Matrix Multiplication framework that supports 8 distinct data types. They leverage only data types that map directly into the DSPs of the Arria 10 target, and therefore support single-precision IEEE-754 and other fixed-point formats. The design is heavily specialized for Deep Neural Networks (DNNs) and delivers a performance of ~99GOps/s and 2.8GOps/s/W energy efficiency for the IEEE-754 32-bit format.

Similarly, Guan et al. [START_REF] Guan | FP-DNN: An Automated Framework for Mapping Deep Neural Networks onto FPGAs with RTL-HLS Hybrid Templates[END_REF] specialize their matrixmultiplication for DNNs by proposing FP-DNN. It is the first work to ease the development of DNNs for FPGAs. Indeed, it generates hybrid RTL-HLS and host source code from the outputs of a trained Tensorflow [START_REF] Abadi | TensorFlow: Large-Scale Machine Learning on Heterogeneous Distributed Systems[END_REF] model. When using the IEEE-754 32-bit format, this design delivers ~81GOps/s and an energy efficiency of ~3.24GOps/s/W .

De Fine Licht et al. [START_REF] De Fine Licht | Flexible Communication Avoiding Matrix Multiplication on FPGA with High-Level Synthesis[END_REF] propose a high level synthesis (HLS) open-source approach for minimizing bandwidth requirements by augmenting the operational intensity. They report performances of 122GOps/s for the floating-point 64bit, 409GOps/s for the 32-bit version, with energy efficiencies of 3.31GOps/s/W , 5.3GOps/s/W , respectively. The proposed designs use DDR4 modules to offload the matrices with DMAs between the HOST and the logic. This creates copies of the data in both directions (user space, kernel space, FPGA DDR4) and partially explains their lower energy efficiencies compared to ours. The authors only leverage the arithmetic data types proposed by Vivado HLS which comprise: FP64, FP32, FP16, uint32, uint16 and uint8. Additionally, since the authors use Vivado HLS, the floating-point types are not IEEE-754 compliant as they do not handle the subnormnals [START_REF] Hrica | Floating-Point Design with Vivado HLS[END_REF].

The fourth row of Table III corresponds to the generated designs running on the FPGA platform described in Section V-A. Sections V-B2 and V-B3 describe the performance and energy efficiency metrics present in this row. Our arrays deliver 19.5GOps/s/W when operating with 32-bit floatingpoint data types, which outperforms the single-precision energy efficiency of the best state-of-the-art FPGA-based GEMM accelerators by 1.86×. Thus, our generator is more energyefficient than state-of-the-art approaches while offering more arithmetic reconfigurability in terms of number of number formats, data bitwidths, and accumulator configurations.

VI. CONCLUSIONS

Our array generator exploits hardware reconfigurability to produce GEMM kernels tailored to and dictated by specific numerical requirements. We conduct an extensive design space exploration and a system-wise evaluation, which show that our generator improves the state-of-the-art in terms of energy efficiency. The generated designs can retain accuracy with large problem sizes (millions of accumulations). Ranging from 240GOps/s/W with 4 accurate significant bits to 0.65GOps/s/W with 1024 accurate bits, we adequately generate drastically different configurations targeting many application domains.

This paper illustrates how the increasing size of reconfigurable hardware, along with emerging communication links, brings many opportunities beyond the IEEE-754 standard, which may produce a paradigm shift in terms of computer arithmetic and hardware heterogeneity.

Fig. 1 :

 1 Fig. 1: Generated array with 3 • 3 Processing Elements.To extract the in-place accumulated value from the outputstationary arrays, several options are possible with their respective trade-offs and requirements. Previous work[START_REF] Khan | [END_REF] exposes two of them: having global lines (buses) from PEs to output buses to form semi-broadcast arrays; or reusing PEs internal registers while stalling the incoming data. None of the induced constraints are suitable for our needs since stalling the incoming pipeline would seriously harm computation and communication overlap and deliver poor performance. Furthermore, given that we aim to produce generic-sized and fast systolic arrays, we can not accept the unpredictable delays caused by high fan-outs yield by global lines, and therefore, we do not consider them.To remedy this problem, our generator leverages a stallingfree output scheme for the output matrix C that allows extracting final and intermediate results through a propagating

Fig. 2 :

 2 Fig.2: Schematic view of a Processing Element (PE) with the vertical carry chain (HSSD) highlighted in red-dashed lines can be translated to S3 and work with the array as long as the corresponding A2S3 unit exists.2) Processing Element (PE): Our generator hardwired algorithm consists in four operations that PEs execute: First, register the incoming S3s values (S3_A and S3_B in Fig.2) arriving respectively from the left and the right, to present them to their corresponding opposite border (right for S3_A_Q and bottom for S3_B_Q). Second, register the control signals and expose them to the opposite border (SOB and EOB from top to bottom). This allows a distributed and local control of the array and avoids global lines, which are often a problem. Third, compute in augmented precision without rounding the dot-product of incoming S3s values in an output-stationary fashion. This is done by a dedicated module discussed in Section III-C3. Fourth, extract via the HSSD mechanism the output-stationary values of C vertically with two registers C1 and C2 and a multiplexer, as Fig.2shows.3) S3 Fused Dot Product (S3FDP): S3FDP is the kernel of our generator as it is the module that performs the fused-dotproduct between values encoded with the S3 format. Fig.3depicts the three main stages of the S3 fused-dot-product. Firstly (Fig.3 1), the exact product is performed. Secondly (Fig.32 ), the mantissa of the product is aligned with regard to the accumulator weights depending on the scale of the product. Finally (Fig.33 ), the fixed-point addition is performed. The last operation requires peculiar timinganalysis. In fact, as fixed-point accumulators require as many bits as the dynamic range of the product contains binades, the carry-rippling can quickly become the critical path. If such, our generator automatically generates Carry-Save-Adders (CSAs) with flopoco timing analysis to suggest the chunk-size.4) Normalizing modules (L2A): Their role is to translate the fixed-point values arriving from the HSSD chains to the

Fig. 3 :

 3 Fig. 3: Circuit design of a generic S3FDP

  posit 16, 0 , posit 16, 1 , and posit 16, 2 16-bit designs; the posit 32, 0 , posit 32, 1 , and posit 32, 2 32-bit designs; and the posit 64, 1 , posit 64, 2 , and posit 64, 3 64bit designs. Our evaluation includes a total of 22 computers number formats.

Fig. 4 :

 4 Fig. 4: Resource utilization of routed S3FDPs at 250MHz for VU3P-2 FPGA in logarithmic scale.

2 )

 2 Cost of PE, SAK and SA: The costs of PE, SAK, and SA arise from the cost of S3FDP components. The extra hardware added to an S3FDP to form a PE does not depend on the frequency and thus can be calculated without the use of tools. This extra cost comes from the delaying of input coefficients of A and B matrices (2 S3-sized registers), control delaying (2 Flip-Flops for SOB and EOB), and HSSD (accumulator-sized mux 2:1 and two accumulator-sized registers). Similarly, with the PE cost we can compute the cost of the SAK that contains n * m PEs as follows: cost SAK (n, m, arith, ovf, msb, lsb) = n • m • cost P E (arith, ovf, msb, lsb).

Fig. 5 :

 5 Fig. 5: Floorplaning of 32•31 and 64•63 arrays with additional PCIe/CAPI/SNAP logic (blue).

Fig. 6 :

 6 Fig. 6: Measured and averaged vs theoretical Throughput (GBytes/s) for different Payload Sizes (Bytes). ffvc1517-2 [3] [66] which embeds 384k LUTs, 788k FFs, 2280 DSPs, 25.3Mb of BRAM, and 90Mb of UltraRAM (URAM).The physical link between the host and the FPGA is PCIe gen4-8lanes, whose 16GT/s provide up to 15.754 GB/s[START_REF]PCI Express[END_REF]. All communications are done on top of the IBM Coherent Accelerator Processor Interface 2(CAPI2[START_REF] Morgan | Opening Up The Server Bus For Coherent Acceleration[END_REF] [47][START_REF] Wikipedia | Coherent Accelerator Processor Interface[END_REF]) protocol with the IBM SNAP framework[START_REF] Ibm | CAPI SNAP Framework Hardware and Software. Contribute to open-power/snap development by creating an account on GitHub[END_REF]. To avoid Clock Domain Crossing (CDC), we generate all the arrays for an operating frequency of 250M Hz, which is the SNAP buses clock.

Fig. 8 :

 8 Fig.7: Measured and averaged vs theoretical Performance (GOps/s) for different array sizes and Payload Sizes (Bytes). describes[START_REF] Lawley | Understanding Performance of PCI Express Systems[END_REF], the hardware PCIe protocol, which contains several layers and encoding schemes, makes it impossible to achieve the peak link bandwidth. In fact, our designs saturate data link bandwidth at their disposal. Fig.8compares throughputs at application level for different array sizes with and without our HSSD extraction chain. The experiment uses a constant payload size of 230 Bytes with a decreasing number of batches b, and therefore, an increasing size of the common dimension p (shown by x-axis of Fig.8). Without HSSD, and thus, without overlaping incoming and outgoing DMAs, an array stalling of n clock edges happens every p clock edges to be able to read out the data. This experiment exhibits that stalling the pipeline for small matrices hinders performances up to two order of magnitude.2) Performance: Fig.7shows the performance achieved by the 64 designs in terms of GOps/s. The x-axis shows data payload in Bytes, and the y-axis displays performance in terms of GOps/s. Our experiments clearly show how the designs experiment 5 possible behaviors in terms of performance depending on the bitwidth. The reason why bitwidths define the performance behavior of the designs is the fact that the number of PEs in the array is defined by the data format bitwidth, as we explain in Section V-B. The designs achieve a performance peak of 4.8, 22.5, 96.7, 399.7, and 1613 GOps/s for 64-, 32-, 16-, 8-and 4-bit numerical formats, respectively. We observe these peaks for payload sizes equal to or larger than 224 Bytes. They are very close the maximum theoretical values, which are 6, 28, 120, 496, and 2016 Gops/s for 64-, 32-, 16-, 8-, and 4-bit numerical formats, respectively, which correspond to arrays with 4 • 3, 8 • 7, 16 • 15, 32 • 31, and 64 • 63 PEs. The theoretical peaks are not achieved due to the physical link bandwidth saturation that we observe for payload sizes equal or larger than 224 Bytes in Fig.6.3) Energy Efficiency: We report the rate of operations delivered by the system per watt of power consumed in Fig.9. We also compare to energy efficiencies delivered by equivalent arrays but with intermediate roundings, thus, mimicking the behavior of traditional FMA hardware. The Figure is split into five categories corresponding to the five arithmetic bitwidths and array sizes that we evaluate. The most efficient design is the 64•63 array of posit 4, 0 β PEs, which delivers 240GOps/s/W . The 4 most performant α-designs are posit 8, 0 α, IEEE-754 8α, TFP8α, and posit 8, 1 α that deliver energy efficiencies of 148.7, 146.6, 146.7, and 136.6 GOps/s/W, respectively, for payload sizes equal or larger

Fig. 9 :

 9 Fig. 9: Energy efficiency (top) against Accuracy (bottom) for different arithmetic bitwidths for different Payload Sizes (Bytes).

TABLE I :

 I Examples of common computer number formats with their corresponding S3 translation.

	Computer arithmetic Value S3 ωS	bias	S3 ωF	S3 quintuple
	IEEE-754 16	0	5	2 ωS-1 -1 = 15	10	0, x, xxxxx, 0, 0000000000
	posit 8, 0	1	4	(N -2) • 2 es = 6	5	0, 0, 0110, 1, 00000
	posit 16, 2	N aN	7	(N -2) • 2 es = 56	11	1, x, xxxxxxx, x, xxxxxxxxxxx
	bfloat16	3.5	8	2 ωS-1 -1 = 127	7	0, 0, 10000000, 1, 1100000
	chain called the Half-Speed Sink Down (HSSD). HSSD comes
	with a resource cost of O(n • m • ωLA) registered elements,
	where ωLA is the bitwidth of the fixed-point accumulator.

TABLE II :

 II Area and power consumption of IEEE-754-α/β/γ S3FDPs and Xilinx FMA IP for VUP-2 FPGA @250MHz

	Bitwidth			8-bit			16-bit			32-bit			64-bit
	Config.	α	β	γ	F M A	α	β	γ	F M A	α	β	γ	F M A	α	β	γ	F M A
	LUTs	85	109	166	233	134	187	200	425	345	1039	416	769	1806 11414 1668	1785
	FFs	18	42	102	342	34	124	102	680	225	865	317	1278	1383	8996	1089	2884
	DSPs	0	0	0	1	1	1	1	1	2	2	2	2	6	6	6	10
	CARRY8s	2	6	13	17	4	16	13	22	12	73	19	37	38	272	34	70
	Power(PE)(W) 0.003 0.004 0.006 0.016 0.008 0.018 0.014 0.032 0.019 0.058 0.023 0.066 0.085 0.460 0.081 0.191

TABLE III :

 III Comparison with previous related works. Box with"-" means that the information is not available.

	work	year	device	freq	% util.	HOST	link	best perf. energy eff. accurate bits	#	reproducibility	generator vendor
			(process node)	(MHz)				(GOps/s) (GOps/s/W) (significant) data types (of computations)		agnostic
								64/32	64/32
	Guan [23]	2017	Stratix V(28nm) [11]	150	95	Xeon E5-2650v2 [31] PCI-e Gen2 x8 [65]	-/~81	-/3.24	max. 23	2	no	yes	-
	Moss [52]	2018 HARPv2(14/20nm) [17]	312.5	99	Xeon E5-26xxv2 [31]	QPI [30]	-/~99	-/~2.8	max. 23	8	no	no	no
	de Fine Licht [15] 2020	VU9P(16nm) [66]	[146-190] [69-90]	-	-	122/409	3.13/10.9	max. 52	6	no	no	no
	this work	2021	VU3P(16nm)	250	[68-92]	POWER9 [50]	CAPI2 [49]	4.8/22.5	5.3/19.5	1024	64	yes	yes	yes

This work has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 955606 Marc Casas is supported by Grant RYC-2017-23269 funded by MCIN/AEI/ 10.13039/501100011033 and by "ESF Investing in your future"