{'fraction_non_alphanumeric': 0.028767645936665396, 'fraction_numerical': 0.007605239730382805, 'mean_word_length': 4.787722655674959, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
FSpGEMM: An OpenCL-based HPC Framework for Accelerating General Sparse Matrix-Matrix Multiplication on FPGAs

General sparse matrix-matrix multiplication (SpGEMM) is an integral part of many scientific computing, high-performance computing (HPC), and graph analytic applications. This paper presents a new compressed sparse vector (CSV) format for representing sparse matrices and FSpGEMM, an OpenCL-based HPC framework for accelerating general sparse matrix-matrix multiplication on FPGAs. The proposed FSpGEMM framework includes an FPGA kernel implementing a throughput-optimized hardware architecture based on Gustavson's algorithm and a host program implementing pre-processing functions for converting input matrices to the CSV format tailored for the proposed architecture. FSpGEMM utilizes a new buffering scheme tailored to Gustavson's algorithm. We compare FSpGEMM implemented on an Intel Arria 10 GX FPGA development board with Intel Math Kernel Library (MKL) implemented on an Intel Xeon E5-2637 CPU and cuSPARSE on an NVIDIA GTX TITAN X GPU, respectively, for multiplying a set of sparse matrices selected from SuiteSparse Matrix Collection. The experiment results show that the proposed FSpGEMM solution achieves on average 4.9x and 1.7x higher performance with 31.9x and 13.1x lower energy consumption per SpGEMM computation than the CPU and GPU implementations, respectively.

However, the existing CPU-and GPU-based solutions to SpGEMM have limited computational performance due to two main reasons. First, an efficient implementation of SpGEMM algorithms (e.g., inner product, outer product, and Gustavson's [14] algorithms) requires a tailored buffering scheme for storing and accessing intermediate results. The deep memory hierarchy and the fixed architecture of CPUs and GPUs make the implementation of such algorithmtailored buffering schemes hardly efficient. Therefore, computing SpGEMM on CPUs and GPUs suffers from limited performance due to poor cache locality and frequent offchip memory access for buffering intermediate results. ond, the complex operations (e.g., sparse vector addition) of SpGEMM algorithms impose strong loop-carried data dependency. The legacy architecture of CPUs and GPUs is able to exploit spatial parallelism but lacks the temporal/pipeline parallelism necessary for resolving such loopcarried data dependency, resulting in large loop initiation intervals that further limits the performance. Additionally, CPU-and GPU-based solutions to SpGEMM for HPC applications suffer from very high energy consumption due to long runtime (i.e. low performance) and high power consumption (e.g., 120W-300W). The high energy consumption of CPUs and GPUs inevitably causes thermal and environmental conservation concerns in HPC and data center computing.
With the recent adoption of FPGAs in data centers [15] as an emerging accelerator, FPGA computing offers an alternative solution to accelerating SpGEMM. An FPGA is a farm of configurable hardware resources whose functionality and interconnection can be redefined at runtime by programming its configuration memory [16]. A state-of-the-art FPGA carries an enormous amount of fine-and coarse-grained logic, computation, memory, and I/O resources. Upon the reconfiguration of these resources, an FPGA can implement any custom hardware architecture tailed to a specific algorithm, often achieving performance and/or energy efficiency gains over CPU and GPU counterparts [17], [18], [19], [20]. The fine-grained logic resources and the abundant on-chip memory and register resources on FPGA devices can be configured to implement a custom buffering scheme tailored for a specific SpGEMM algorithm to allow efficient on-chip storage of intermediate results and data movement among processing elements (PEs) to avoid off-chip memory access [21], [22]. In addition, the reconfigurable resources on FPGA devices can compose not only spatial but also temporal/pipeline parallelism both at a fine granularity and on a massive scale to resolve the complex loop-carried data dependency that exists in SpGEMM algorithms to minimize loop initiation intervals for improved performance [23], [24]. Furthermore, while providing higher performance for executing high-dependency algorithms, FPGAs generally have lower power consumption than CPUs and GPUs, which directly translates into lower energy consumption per task (higher energy efficiency). So far, there has been limited work on accelerating SpGEMM on FPGAs [25], [26].
Existing work all adopts the inner product algorithm [27] in their implementations. The inner product algorithm attempts to compute all zero and nonzero output values. In the case of sparse matrices, there are a considerable amount of computations that result in a zero output, consuming clock cycles that can be spared with domain knowledge embedded into the architecture. Additionally, the dot product operation between a row and a column of the input matrices requires index matching, which further contributes to the overhead of the SpGEMM algorithm. Therefore, the inner product algorithm is not suitable for SpGEMM.
In this paper, we propose a new compressed sparse vector (CSV) format for representing sparse matrices and FSpGEMM, an OpenCL-based HPC framework for accelerating SpGEMM on FPGAs. The proposed FSpGEMM framework consists of an FPGA kernel implementing a throughput-optimized hardware architecture and a host program implementing pre-processing functions for converting input matrices to the CSV format tailored for the proposed architecture. Different from prior work, we adopt Gustavson's algorithm to avoid zero output computation and reduce the synchronization overhead of computing partial products. Benefiting from the hardware flexibility of FPGAs, we propose a custom buffering scheme tailored to Gustavson's algorithm to improve the reuse of input matrices, thus largely reduce the amount of off-chip memory access. Additionally, we co-design FSpGEMM with the CSV format to transform the memory access pattern of input matrices from irregular to regular, which improves the memory bandwidth utilization and eliminates unnecessary memory stall. Overall, such synergies between the buffering scheme and Gustavson's algorithm as well as between the CSV format and FSpGEMM as a result of the co-design methodology significantly improve the computational performance. The contributions of this paper are summarized as follows. • We propose a new CSV format for representing and storing sparse matrices in a format that is codesigned with FSpGEMM for reducing unnecessary off-chip memory access and improving performance.

•
We propose an end-to-end OpenCL-based HPC framework for accelerating SpGEMM on FP-GAs, consisting of an energy-efficient and deeply pipelined FPGA kernel for accelerating the Gustavson's algorithm and a host program implementing pre-processing functions for converting and/or storing raw input matrix files into the CSV format required for the proposed architecture. • We propose a parameterized and scalable FPGA hardware architecture and an analytical method for determining the optimized architectural parameters for maximizing the run-time performance of SpGEMM given an arbitrary FPGA board specifications. • We propose a new buffering scheme to improve the data reuse of Gustavson's algorithm and reduce the amount of off-chip memory access, which improves the overall computational performance.

•
The proposed CSV format and framework can be potentially applied to all host-accelerator computing models with the SpGEMM computation offloaded to accelerator devices with customizable hardware (e.g., an FPGA or ASIC).

•
We evaluate the performance and energy efficiency of FSpGEMM based on an Intel Stratix 10 GX FPGA development board and compare it with Intel Math Kernel Library (MKL) and cuSPARSE, the state-of-the-art libraries for SpGEMM, running on an Intel Xeon E5-2637 CPU and an NVIDIA GTX TITAN X GPU, respectively, for multiplying the poisson3Da, 2cubes sphere, filter3D, cage12, scircuit, mac econ fwd500, offshore, and webbase-1M matrices selected from SuiteSparse Matrix Collection [28]. The experiment results show that the proposed FSpGEMM solution has on average 4.9× and 1.7× higher performance with 31.9× and 13.1× lower energy consumption per SpGEMM computation than the CPU and GPU implementations, respectively.
The rest of this paper is organized as follows. Section 2 provides background on different SpGEMM algorithms and summarizes the related work on accelerating linear algebra kernels on FPGAs as well as the state-of-the-art implementations of SpGEMM on CPUs, GPUs, FPGAs, and ASICs. Section 3 elaborates the proposed CSV format. Section 4 describes the proposed framework's buffering scheme, hardware architecture, and the host program. The evaluation of FSpGEMM is presented in section 5. Finally, the paper is concluded in Section 6. Figure 2) shows a compressed representation of a sparse matrix using Compressed Sparse Row (CSR) format [29]. In this format, nonzero values of the sparse matrix are laid out in the row-major orientation in the off-chip memory. The CSR format stores a sparse matrix using three arrays V , COL IN DEX, and ROW P T R representing nonzero values and column index of the nonzero elements, and the pointer to the first nonzero element in the first two arrays, respectively.

Sparse Matrix Formats
In Compressed Sparse Column (CSC) [29] format, the nonzero elements are stored in the column-major orientation using three arrays V , ROW IN DEX, and COL P T R for nonzero values, row index, and the pointer to the start of each column.
These data formats are not tailored to a SpGEMM algorithm or hardware architecture. Hence, they are not efficient for specific methods or hardware designs, which leads to a huge performance loss. Differently, the new CSV format makes the input data access in the proposed buffering scheme regular when utilizing the CSR or CSC formats.

SpGEMM Algorithms
There are three main methods for computing SpGEMM: inner product, outer product, and Gustavson's method. The differences among the three methods are twofold. First, these methods require different data formats for the input matrices to acquire contiguous access to off-chip memory. Second, some methods avoid wasted computations by calculating the nonzero elements of the output matrix.
The inner product algorithm [27] computes all the elements of the output matrix, including zero elements. In SpGEMM, most of the output elements are nonzeros. Additionally, computing each element involves a dot product operation, including index matching and multiply-accumulate (MAC). Jamro et al. [26] identifies index matching as a hardware-expensive and the most time-consuming operation of the inner product method. Therefore, the inner product algorithm inevitably causes its implementations to summer from both performance and energy consumption overheads.
The outer product algorithm [27] performs an outer product operation between a column (row) of the first input matrix and row (column) of the second input matrix. The result of each outer product operation is a large partial sum matrix with the same dimensions as the input matrices. The number of partial sums produced is equal to the number of rows of the input matrices that are often large. Therefore, buffering and accessing partial sums requires offchip memory access that incurs long access latency and consumes high energy. Moreover, the addition of large partial sums suffers from synchronization overhead. Consequently, the outer product algorithm suffers from undesired performance and energy consumption overhead. Figure 1 illustrates the row-wise Gustavson's method for multiplying two sparse matrices (A × B). In this method, each non-zero element in a row of the first input matrix (e.g., A(i, j) where i and j are the row and column indices, respectively) is multiplied by all non-zero elements of the corresponding row of the second input matrix (e.g., B(j, :) where : is the slice operation), resulting in a intermediate row of sparse partial products (e.g., C j,temp (i, :)). The addition of the sparse partial products from the multiplication of all nonzero elements in a row (e.g., the ith row) of the first input matrix with the corresponding rows of the second input matrix results in a final row of the output matrix (e.g., C(i, :) = j C j,temp (i, :)).  The column-wise Gustavson's method is similar to the column-wise one but with rows and columns switched. Gustavson's method does not require the hardware-expensive index-matching operation among the elements of input matrices. Also, in Gustavson's method, the addition operation of sparse partial products has low on-chip memory requirement and synchronization overhead. In contrast, since an entire intermediate matrix is dealt with using off-chip memory in the outer product method, utilizing Gustavson's method results in both improved performance and lower energy consumption. Therefore, we choose the row-wise Gustavson's method to develop FSpGEMM upon.
However, Gustavson's method maintains sub-par reuse of the second input matrix's data due to the irregular access pattern of rows that are read from off-chip memory. Such an access pattern is dependent on the order of non-zero column indices in the rows of the first input matrix, making the caching policy or buffering scheme for the second input matrix too complex and costly to implement.

Accelerating Dense GEMM on FPGAs
Parametrized FPGA implementations of dot-product and matrix-vector multiplication kernels are presented in [30]. The work in [30] also compared the proposed kernels with CPU and GPU implementations in terms of performance and energy efficiency. FBLAS [31] proposes scalable, modular, and OpenCL-based implementations of the Basic Linear Algebra Subprograms (BLAS) library to improve the reusability of the FPGA kernels. This work on accelerating the BLAS library targets dense vectors and matrices thus are not suitable nor efficient for accelerating SpGEMM.

Accelerating SpGEMM on FPGAs
The work in [25] and [26] proposes the design and implementation of the inner product algorithm for SpGEMM. They study the performance and energy consumption tradeoff of the design by tuning the architectural parameters (i.e., the number of PEs and the block size in blocking decomposition). Changing architectural parameters results in different FPGA resource utilization. All the existing work adopts the inner product method for SpGEMM, thus suffers from undesired performance and energy consumption overheads due to the costly index matching operation. Differently, we adopt Gustavson's algorithm eliminating expensive index matching operation, zero output computations, and storage of large intermediate values resulting in significant performance and consequently energy efficiency improvements.

Accelerating SpGEMM on Application-specific Integrated Circuits (ASICs)
There has been some recent work on accelerating SpGEMM on ASICs using either the outer product or Gustavson's method. OuterSPACE [32], [33] and SpArch [34] utilizes the outer product method. SpArch reduces the number of partial output matrices by matrix condensing to mitigate the overheads of synchronization and off-chip memory access. However, both OuterSPACE and SpArch still require off-chip memory access to store and retrieve intermediate results (partial output matrices). MatRaptor [27] uses the Gustavson's method to accelerate SpGEMM. Nonetheless, this work suffers from poor data reuse of the input matrices, leading to unnecessary off-chip memory access for loading input matrices repetitively with a large performance penalty. Moreover, the scalar-vector multiplication and the addition of sparse vectors in MatRaptor [27] are performed with very low computational parallelism (1 DSP slice per PE), leading to limited performance.

CPU and GPU Implementations of GEMM and SpGEMM
Intel Math Kernel Library (MKL) [35] is highly optimized for Intel processors and includes threaded kernels, such as GEMM and SpGEMM, suitable for use with computationally intensive applications. MKL uses AVX, AVX2, AVX-512 SIMD instructions of Intel processors to improve performance. CSparse [36] is a sparse matrix package included in SuiteSparse. SuiteSparse is a suite of sparse matrix algorithms [37] that consists of sparse matrix operations. However, its SpGEMM kernel is single-threaded and does not utilize SIMD instructions.
The cuSPARSE library [38] contains a set of basic linear algebra subroutines used for handling sparse matrices implemented on top of the NVIDIA CUDA runtime as part of the CUDA Toolkit and is designed to be invokable from C and C++. MKL and cuSPARSE contain the CPU and GPU implementations of SpGEMM, respectively, and are used as the reference SpGEMM computing solutions (in the context of HPC) for comparison in this study.

COMPRESSED SPARSE VECTOR (CSV) FORMAT
As GEMM and SpGEMM are primarily computed with spatial parallelism in CPU, GPU, and most custom hardware accelerators, the input matrices are often read from off-chip memory and fed into the computation units one vector at a time, where the sparse vector length is equal to the number of the computation units. In the row-wise Gustavson's method, reading the nonzero elements from multiple rows of the input matrix based on the CSR format in a vector fashion (see Figure 2) leads to a non-continuous off-chip memory access pattern. The same problem exists for the column-wise Gustavson's method when reading multiple columns of the input matrix using based on the CSC format [29]. A non-continuous off-chip memory access pattern results in a large performance penalty. Additionally, using the CSC and CSR formats for the row-and columnwise Gustavson's algorithms, respectively, requires having a very large lookup table (in the size of the input matrix dimension) to keep track of the locations of the last nonzero entry read from the input matrix, which is too costly to implement for large matrices. To address these issues, we propose the CSV format tailored to Gustavson's algorithm for reading input matrices in a vector fashion. Such a vector in the length of the number of computation units is referred to as a CSV vector in the context of the CSV format.
The CSV format uses three attributes for representing each nonzero element: the value V AL, the row index ROW IN DEX, and the column index COL IN DEX. Therefore, the location of the last nonzero entry read from the input matrix is always clearly indicated as (ROW IN DEX, COL IN DEX) without needing a loopup table for storage. In addition, the nonzero elements in the CSV format are stored in a new vector-major order (see Figure 2), which is a key difference from the CSR or CSC format that adopts a row-or column-major order. Storing nonzero elements in a vector-major order (that matches the number of computational units) is the key to assuring a continuous off-chip memory access pattern, which improves both the off-chip memory bandwidth utilization and computational performance.

FRAMEWORK DESIGN
The FSpGEMM framework consists of an FPGA kernel written in OpenCL and a host program running on an FPGA accelerator and a host CPU, respectively. FSpGEMM utilizes a new buffering scheme tailored to Gustavson's algorithm to reduce the amount of off-chip memory access. The kernel code implemented on the FPGA accelerator performs the SpGEMM computation based on Gustavson's algorithm. On the host side, we provide the utility functions for preprocessing and storing raw matrix files into the CSV format and the OpenCL API for provisioning tasks running on the FPGA accelerator.

Data Buffering Scheme
We adopt the row-wise Gustavson's method as the baseline algorithm. To improve the data reuse of the second input matrix and reduce the amount of off-chip memory access, we propose to process multiple rows of the first input matrix in parallel using multiple computing units (i.e., one row per computing unit) while sharing a row of the second input matrix among all computing units. To enable that, a CSV vector of nonzero values of the first input matrix is read. Based on the column index of the vector, the corresponding row from the second input matrix is then read and buffered in the on-chip memory. As a result, the access to the buffered row of the second input matrix is reused for multiple rows of the first input matrix.
The example in Figure 3 shows how the proposed data buffering scheme avoids one round of off-chip memory access to 4th row (since the grayed out CSV vector is in the 4th column of the first input matrix) of the second input matrix by reusing this buffered row for all computing units.
Since for each nonzero element of the first input matrix, a row of the second input matrix is read and reused in the burst mode, we define the total number of off-chip memory access to the second input matrix as the total number of nonzero elements in the first input matrix. Consequently, We define off-chip memory access reduction (OMAR) percentage as where v is the CSV vector index, V is a set of indices of all nonzero CSV vectors (shown as colored in Figure 2), A(v) is a nonzero CSV vector indexed by v, and nnz(·) is the operator that measures the number of nonzero elements in a sparse vector or matrix.
Note that the proposed data buffering scheme can be implemented on all customizable hardware (e.g., FPGA and ASIC devices). Also, it can potentially be applied to any fixed-architecture device with a compatible memory hierarchy. Figure 4 shows a high-level block diagram of the hardware architecture of the FPGA kernel, including three types of kernels: PEs, load and store kernels. The PEs are responsible for computing the nonzero elements of the output matrix based on Gustavson's algorithm. The load and store kernels read and write input and output data from and to the off-chip memory as well as feed and receive data to and from the PEs via FIFO channels, respectively. Separating the load and store kernels from PEs and connecting them using FIFO channels facilitates the data flow control for balancing the off-chip memory bandwidth and the data processing throughput by PEs. The depth of the FIFO channels is optimized by the offline compiler such that the load and store kernels are able to continuously read and write data from and to the off-chip memory, respectively, while pushing data throughout the kernel pipeline.

Architectural Overview
Each PE is responsible for calculating one row of the output matrix at a time. Multiple PEs process multiple rows independently in parallel. Each PE receives a nonzero scalar value of the first input matrix along with additional scheduling information through a QA channel. Also, each PE gets nonzero values from the corresponding row of the second input matrix from a QB channel. The output results are streamed to the store kernel via a QC channel.
To parallelize the computation of multiplications in a vector fashion, we take advantage of the loop unrolling optimization directive provided by OpenCL Offline Compiler. The FPGA kernel in FSpGEMM includes two types of OpenCL kernels: 1) autorun kernel, which starts executing automatically and does not need to be managed by a host program, and 2) host-driven kernel. As the offline compiler provided by Intel FPGA SDK for OpenCL does not need to implement any glue logic for communicating with the host program, using autorun kernels reduces logic utilization and allows for the mapping of more computation resources for improved performance [39]. The PEs are autorun kernels, and the load and store kernels are host-driven kernels since they access the off-chip memory that must be allocated first by a host program. Overall, the FPGA kernel is parameterized with two architectural parameters: SW (SIMD width of PEs) and N U M P E (number of PEs), which allows users to balance the trade-off between computational parallelism and FPGA resource usage.

Load and Store Kernels
One should note that the first input matrix and the output matrix are represented by the CSV format, and the  second input matrix is represented by the CSR format in FSpGEMM. The load kernel iterates over the total number of nonzero elements in the first input matrix. In each iteration, the load kernel first sends a data structure (A DS) containing the value of a nonzero element of the first input matrix along with the other attribute data described in Table 1 using a QA channel. Specifically, RESET indicates whether the received nonzero element is the last one in the current row of the first input matrix. If so, RESET = 1 signals the PEs to finish the computations and reset to their initial states. The load kernel takes advantage of the CSV format and identifies the nonzero elements of the first input matrix in the same CSV vector by comparing the column indices of two consecutive nonzero elements. If the column index of the next nonzero element is different from the column index of the current nonzero element, the load kernel buffers the nonzero element of the corresponding (based on the current nonzero element current index) row of the second matrix and sends it to the corresponding PEs based on the row indices of nonzero elements of the first input matrix using the CSV format. The load kernel also calculates the number of nonzero elements in a row of the second input matrix that the PEs need to expect from the load kernel. Since the values from a row of the second input matrix are loaded in a vector fashion and get fed into all the PEs, the load kernel sends N U M B V EC vectors of nonzero values each in the size of SW in parallel. Additionally, the row index of the nonzero element from the first input matrix (A ROW IN D) is sent to the PEs to determine the row index of the output matrix.
Since the load kernel reads a row of the second input matrix at a time, the data are read from the off-chip memory using the CSR format to enable contiguous and regular access. Each vector of the second input matrix is represented by the data structure B DS as described in Table 2. V AL and B COL IN D are both a vector in the size of SW that contains the values and the column indices of the nonzero elements in a row of the second input matrix.
The store kernel reads the computation results of the output matrix from the PEs, which is represented by the data  structure C DS, including the nonzero value and its row and column indices as described in Table 3. Since the valid output values of different PEs are not necessarily produced at the same time, the FIFO channels QC are used to facilitate the data flow control and assure data integrity. Figure 4 shows the high-level block diagram of a PE. Each PE is responsible for two major operations for producing each row of the output matrix: 1) scalar-vector multiplication of a nonzero value of the first input matrix and a row of the second input matrix resulting in an intermediate row of sparse partial products and, 2) addition of sparse partial products. These functionalities are implemented by three main submodules: a vectorized multiplication (VecMult) unit for computing sparse partial products with a SIMD width of SW , a sort and merge (SM) unit for sorting and accumulating partial products, and a memory unit implementing a double buffering scheme for caching intermediate results. The data structure C T EM P DS of type C DS is used to store the temporary computed output elements.

Processing Elements (PEs)
Algorithm 1 describes the OpenCL pseudo-code implementation of the SM unit. The corresponding indices for C T EM P V EC are available in B V EC IN D. The SM unit iterates over all the elements of C T EM P V EC. In each iteration, it first checks if there is a valid value in the buffers to add C T EM P V EC with C T EM P ROW by comparing the head pointer (i.e., M EM P T R HEAD[S]) and the tail pointer (i.e., M EM P T R T AIL[S]). In this case, the SM unit compares the index of the current element of C T EM P V EC with the index of stored values indicate by M EM P T R HEAD [S]. Since the C T EM P V EC and C T EM P ROW are already sorted, the value with the smaller index has the smallest index. When the smaller index exists in the buffer, the corresponding head pointer will be increased by one. In the case that the indices from C T EM P V EC and C T EM P ROW are equal, the merge operation occurs where the corresponding values from C T EM P V EC and C T EM P ROW are added. In the consequent iterations, if the head pointer is equal to the tail pointer, the SM unit has checked all the elements stored in the buffer and the values of C DS are set based on the current C T EM P V EC.
All operations described in Algorithm 1 happens B N U M V EC times. At the end of each iteration, the double buffer selector (S) is switched between 0 and 1. Finally, there might be some remaining elements in the current buffer determined by S. This can be checked by comparing the head and tail pointers. All C DS elements are set equal to the remaining values stored in the buffer with no comparison. Whenever finished processing on the vectors of the second input matrix or the remaining elements in the buffer, the PE decides whether to store C DS in the other buffer or send the final value to a QC channel based on the value of A DS.RESET .

Performance Optimization for Determining Architectural Parameters
We provide a guideline to derive the optimal architectural parameters (SW and N U M P E) for minimizing the runtime subject to the available resources on a given FPGA device. We formulate the runtime (R) as R = N Ops F ×P ×U , where N Ops , F , P , and U are the number of floating-point operations (FLOPs) to produce an output matrix, the clock frequency, the computational parallelism, and the spatialtemporal utilization factor (STUF), respectively. STUF is a statistical metric that measures on average how busy the available computation resources in a computing device are (in both space and time) for performing useful operations of a given algorithm. Thus, STUF is defined as the ratio of effective computational parallelism per clock cycle on average and ranges from 0 to 1. N Ops is an algorithmic parameter and depends on the dimensions and the number of nonzero elements of input matrices. Computational parallelism is determined by the total number of utilized DSP units calculated as SW × N U M P E. Therefore, R can also be expressed as where N Ops , F , and U are lumped as α. We define the constrained optimization problem as where f 1 (SW ) = sizeof (f loat) × SW × F and C 1 are the required and available off-chip memory bandwidth, respectively, with sizeof (f loat) = 32b. Function f 2 (SW, N U M P E) = β × P = β × SW × N U M P E and C 2 are the total required and available logic resources, respectively, where β is a linear fitting parameter that captures the proportionality of the actual logic resource usage over computational parallelism that can be derived based on the logic resource (f 2 ) usage reported by the offline compiler given a target FPGA device.
Based on R(SW ) = α SW ×N U M P E , the runtime R is minimized when SW and N U M P E are maximized. Thus, the constrained optimization problem defined in Equation 3 has an analytical solution that can be derived as follows.

Host Program
The host program running on the CPU includes two major parts: pre-processing utility functions and OpenCL API function calls. As the first input matrix needs to be transformed and stored in the CSV format prior to the computation by the FPGA kernel, the utility functions read in the raw matrix files in an existing sparse matrix format then convert and store the matrices in the CSV format. The preprocessing step only needs to be performed once. We have utilized OpenCL API functions to create memory buffers, enqueue kernels, and read the results back from the FPGA device.

Experiment Setup
We evaluate the performance and energy efficiency of FSpGEMM in terms of runtime (s) and energy consumption (J ) per SpGEMM computation. To evaluate the FSpGEMM framework, we select a set of sparse matrices from the publicly available SuiteSparse Matrix Collection [28] (formerly known as the University of Florida Sparse Matrix Collection), a collection of sparse matrices in real applications. The specifications of the matrices are summarized in Table 4 including matrix dimensions (the number of rows times the number of columns) and the density calculated as . Our design adopts the single-precision floating-point (32-bit) data format. We develop and implement FSpGEMM using OpenCL with Intel FPGA SDK for OpenCL with Quartus Prime Pro 20.1. We compare FSpGEMM with MKL and cuSPARSE libraries for the CPU and GPU implementations, respectively. We measure the performance of the CPU implementation of MKL on a dual-socket Intel Xeon E5-2637 v3 CPU [40], and the GPU implementation of cuSPARSE on   [41]. Table 5 summarizes the specifications of the CPU system, the GPU device, and the FPGA board used in the evaluation.

OMAR Evaluation
As mentioned in Section 3, a naive implementation of Gustavson's algorithm suffers from poor reuse of the second input matrix. To address this issue, we propose to process multiple rows of the first input matrix in parallel and share the rows of the second among multiple PEs. We calculate the off-chip memory access reduction percentage (OMAR %) using Equation 1 for the matrices summarized in Table 4. Figure 6 shows the OMAR percentage that can be achieved by the proposed data buffering scheme with respect to different numbers of PEs (i.e., parallel rows) and matrices. The results show that the amount of OMAR that can be archived by the proposed data buffering scheme monotonically improves as the number of PEs increases. 1.7%-24.8%, 6.0%-38.6%, 15.9%-46.5%, 28.1%-51.3%, and 39.2%-54.0% OMAR can be achieved at the PE number of 2, 4, 8, 16, and 32, respectively, across the selected sparse matrices. The improvement with increasing the number of PEs is because a row of the second input matrix is buffered and shared for processing more rows of the first input matrix. Thus, the amount of access to the row of the second input matrix is decreased.

Experiment Results
The optimal architectural parameters of the FPGA kernel implementation derived from the performance optimization method presented in Section 4.2.4 are SW = 16, N U M P E = 32. The FPGA kernel implementation runs at 236 MHz. Table 8 shows the resource utilization for the implemented FPGA kernel. As projected by our performance and architecture models, the FPGA kernel implementation fully utilizes the FPGA resources (limited by logic  utilization), and further increasing SW beyond 16 does not improve performance any further due to the saturation of off-chip memory bandwidth. Table 7 shows the performance comparison of FSpGEMM in terms of runtime in ms with the CPU implementation of MKL and the GPU implementation of cuSPARSE for performing SpGEMM computation. A lower runtime indicates higher performance. Figure 7 shows the performance speedup achieved by FSpGEMM over the CPU and GPU implementations with respect to different sparse matrices. Overall, the speedup over the CPU and GPU implementations ranges from 1.1×-9.7× and 0.6×-3.0×, respectively. On average, FSpGEMM achieves 4.9× and 1.7× higher performance than the CPU and GPU counterparts, respectively. It should be noted that such speedup ratios are achieved at the condition that the FPGA implementation runs at a 5-15× lower clock frequency than the CPU and GPU implementations. The speedup mainly stems from the offchip memory access reduction due to the customized data buffer scheme tailored to Gustavson's algorithm, the improved off-chip memory bandwidth utilization as a result of the architectural co-design with the proposed CSV format, the temporal/pipeline parallelism offered by FPGAs for resolving the strong data dependency in the Gustavson's algorithm, and the throughput-optimized hardware architecture tailored to the Gustavson's algorithm, all of which improves the overall STUF of the computation resources on the FPGA at run time.

STUF Comparison with CPU and GPU Implementations
In section 4, we defined the runtime as R = N Ops F ×P ×U , where N Ops , F , P , and U are the number of FLOPs in the Gustavson's algorithm, clock frequency, computational parallelism, and STUF, respectively. We consider computational parallelism as the number of possible floating-point operations per clock cycle. A DSP in the FPGA and a CUDA core in the GPU [42] can perform one multiplication and one addition per clock cycle. Therefore, the computational parallelism for the FPGA and GPU is equal to twice the number of DSPs (1,518 DSPs) and twice the number of CUDA cores (3,584 cores), respectively. The CPU has two sockets with four cores per socket, and each core provides 32 FLOPs per cycle, enabled by Intel AVX2 instruction set extension [43]. Therefore, we can estimate the actual STUF of different solutions based on the runtime measurements as U = N Ops F ×P ×R .

Energy Efficiency Comparison with CPU and GPU Implementations
We measure the average power consumption of the Intel Arria 10 GX development board using the Power Monitor tool in the Board Test System (BTS) application provided by Intel [44] during FPGA kernel execution. BTS measures the supply voltage and the drawn current of the entire FPGA board by reads values from onboard sensors. For the power measurement of the CPU, we utilize the likwid-powermeter tool from Likwid [45] to access the  Running Average Power Limit (RAPL) counters on the Intel CPU. The RAPL interface is controlled through MSR registers [46]. For the power measurement of the GPU, we utilize the POWER query option [47] of the NVIDIA System Management Interface (nvidia-smi) [48] tool. Table 9 summarizes the energy consumption (J ) per SpGEMM computation of different implementations, which is calculated as product of runtime (s) and average power consumption (W ). Figure 8 shows the corresponding energy consumption reduction achieved by the FPGA implementation. Overall, the energy consumption reduction over the CPU and GPU implementations ranges from 7.6×-58.6× and 4.3×-22.3×, respectively. On average, FSpGEMM achieves 31.9× and 13.1× higher energy efficiency than the CPU and GPU counterparts, respectively.

CONCLUSION
In this paper, we present FSpGEMM, an OpenCL-based HPC framework for FPGA acceleration of SpGEMM. FSpGEMM includes a deeply pipelined and scalable FPGA kernel that accelerates the SpGEMM algorithm and a set of utility functions for pre-processing step of input matrices into the proposed CSV format.
We introduce a performance-optimized model to derive architectural parameters for the FPGA kernel subject to the available on-chip resources (DSPs and memory blocks) to map the design into the arbitrary FPGAs. The experiment results based on the Intel Arria 10 GX FPGA development board for accelerating SpGEMM of a set of sparse matrices from SuiteSparse Matrix Collection show on average one order of magnitude higher performance and lower energy consumption compared to the state-of-the-art implementations of SpGEMM on CPU and GPU.