Sparse matrix-vector multiplication acceleration method facing dawn GPU (Graphics Processing Unit)

Abstract

The invention discloses a sparse matrix-vector multiplication acceleration method for a dawn-bathing GPU, and belongs to the field of high-performance calculation. The method comprises the following steps: acquiring a sparse matrix A and a vector x in a CSR format; dividing A into 16 * 16 sub-blocks Tile according to the width of an operand of a Xichen GPU matrix multiplied instruction, and making a proportion eta < lt > for non-zero elements; a CSR sub-matrix rowbase list is generated by 10% of Tiles, and a Tile-Dense sub-matrix block base list is generated by the rest of the Tiles; rows of the CSR sub-matrix are divided into light rows, medium rows and weight rows according to the number of non-zero elements, three-level protocols are completed in the same GPU kernel through Thread Block Tile, on-chip shared memory and global memory atomic operation, and ycsr is obtained; an mfma instruction is called for the Tile-Dense sub-matrix, and accumulation is carried out to obtain ytile; and finally, accumulating the ycsr and the ytile according to rows to output a result y. According to the method, an mcSPARSI interface can be replaced in a binary level, and load unbalance is relieved and the matrix multiplication instruction utilization rate is improved on the premise that an upper frame is not changed.

Description

Technical Field

[0001] This invention relates to the field of high-performance computing, and in particular to a sparse matrix-vector multiplication acceleration method for Muxi GPUs, which can be widely applied to scientific computing, deep learning inference, graph neural networks and large-scale graph analysis. Background Technology

[0002] Sparse matrix-vector multiplication (SpMV) is a core computational operator in applications such as scientific computing, deep learning inference, graph neural networks, and large-scale graph analysis. With the exponential growth of problem size, SpMV accounts for a significantly larger share of time in end-to-end applications, becoming one of the key bottlenecks restricting overall system performance.

[0003] To fully leverage the parallel computing capabilities of graphics processing units (GPUs), academia and industry have proposed various SpMV optimization strategies, such as MERGE-CSR and ELL-C-SIGMA, based on general-purpose GPU architectures like NVIDIA and AMD. While these optimization methods have achieved good results on specific platforms, their high coupling with the underlying instruction set, memory hierarchy, and thread scheduling mechanism often leads to performance degradation or even functional abnormalities when directly migrated to Muxi GPUs designed for artificial intelligence and scientific computing scenarios.

[0004] In the sparse linear algebra library mcSPARSE provided by Muxi, the mcsparseCsrmvEx interface performs well in most typical scenarios. However, its performance degrades significantly when the distribution of non-zero elements in each row of the sparse matrix differs considerably. Analysis reveals that Muxi GPUs face two main challenges in implementing SpMV: Load imbalance problem: Sparse matrices, which are common in scientific computing, have extremely large variance in the distribution of non-zero elements in the row dimension. This leads to an extremely uneven distribution of tasks among the computing units inside the Muxi GPU many-core processor. Some threads are overloaded while a large number of threads are idle or waiting, which seriously affects parallel efficiency.

[0005] Low utilization of computing units: The sparse mode lacks spatial locality, making it difficult to trigger the matrix multiplication acceleration instructions that occupy the main computing power in the Muxi GPU, resulting in the core computing resources not being fully utilized and low computing power utilization.

[0006] Therefore, there is an urgent need for a SpMV acceleration method specifically designed for the characteristics of the Muxi GPU architecture, which can perceive the structural features of sparse matrices online, dynamically balance the computational load, and fully activate matrix multiplication acceleration instructions to achieve robust performance improvement across matrix distributions. Summary of the Invention

[0007] When the existing MXGPU executes the sparse matrix-vector multiplication (SpMV), due to the extremely uneven distribution of non-zero elements in the rows of the sparse matrix, the computational workload within the Thread Block Tile varies significantly, resulting in a low utilization rate of the matrix multiplication acceleration instruction mfma and a drastic fluctuation in the overall parallel efficiency with the change of matrix distribution.

[0008] To solve the above problems, the present invention provides a SpMV acceleration method for MXGPU, which includes the following steps: Step S1: Obtain the sparse matrix in CSR format A ∈R^( M × K ) and the vector x ∈R^ K , and allocate video memory on the device to store the row offset rowPtr , column index colIdx , matrix value val , record the row number M , column number K and the number of non-zero elements nnz ; Step S2: According to the operand width of the MXGPU matrix multiplication instruction mfma.f32.16x16x16.f16, divide A evenly into 16×16 sub-block Tiles, define NNZ_PER_TILE as the number of non-zero elements in the Tile, calculate η = NNZ_PER_TILE / 256; if η < 10%, then move all the non-zero elements of this Tile into row_base_list and record it as the CSR sub-matrix, otherwise move it into block_base_list and record it as the Tile-Dense sub-matrix; Step S3: Define the row number of the sparse matrix r , and count the number of non-zero elements in each row NNZ( r ), set the thresholds th1 and th2 , divide the rows into: lightweight rows L1: NNZ(r) ≤ th1; medium-weight rows L2: th1 < NNZ(r) ≤ th2; heavyweight rows L3: NNZ(r) > th2; and generate L1_rows , L2_rows , L3_rows arrays that only save the row numbers; Step S4: Parallelly execute the CSR sub-matrix within the same GPU kernel: For L1 rows, taking the Thread Block Tile as the unit, read into the on-chip shared memory through merged memory access, use __shfl_xor_sync to complete register-level reduction, and write the result into y_csr r ​]; For line L2, each Thread Block Tile is responsible for one line, defining the thread number within the Thread Block Tile. tid ∈[0,1023], perform a binary tree reduction of depth 10 in shared memory, by tid =0 threads write y_csr [ r ]; Enable for L3 line T =NNZ(r) / 1024 Thread Block Tile cooperative local reduction, each Tile uses a global memory int32_t counter. C [ r ]Execute atomicAdd(&C[ r ],1), when C [ r ]= T The summation specification is completed and written by the last tile. y_csr [ r ]; Step S5: For the Tile-Dense submatrix, start the GPU kernel. Each Thread Block Tile is responsible for a 16×16 Tile. Load 256 float16 values ​​into the register, call the mfma_f32_16x16x16_f16 instruction to complete the local matrix-vector multiplication, and write the result to a temporary array. y_tile ; Step S6: Start the GPU kernel, and... y_tile and y_csr The final SpMV result is obtained by summing the row numbers. y . Attached Figure Description

[0009] The present invention will be further described below with reference to the accompanying drawings and embodiments.

[0010] Figure 1 This is a flowchart illustrating the implementation of a high-performance SpMV algorithm based on the Muxi GPU according to an embodiment of the present invention. Figure 2 A thumbnail of the sparse matrix impcol_e; Detailed Implementation

[0011] The present invention will be further described in detail below with reference to the accompanying drawings and optional embodiments. These embodiments are only used to fully disclose the technical solution of the present invention so that those skilled in the art can reproduce it, and do not constitute a limitation on the scope of protection of the present invention.

[0012] This embodiment uses the SuiteSparse Matrix Collection (https: / / sparse.tamu.edu / ) as the test set. This dataset contains 2,893 sparse matrices from real-world engineering, scientific, or commercial scenarios, spanning more than 40 application domains, and is a recognized reproducible benchmark library in the numerical linear algebra community. Figure 2 A thumbnail of the “impcol_e” matrix is ​​given. This matrix is ​​derived from power system simulation. It has a large number of non-zero row elements with large variance and lacks spatial locality. It can be used to verify the effectiveness of the present invention under extreme distributions.

[0013] The hardware platform was a Muxi MXC500 32 GB graphics card, the software stack was MXMACA SDK 2.3, driver version 5.2.0, and compiler MXCC 3.1. Unless otherwise specified, all parameters below were measured under this environment.

[0014] Figure 1 This is the overall flowchart of this embodiment. Steps S1-S6 complete the SpMV calculation of "load balancing + matrix multiplication instruction friendly" in one go, and can directly replace the mcsparseCsrmvEx interface of mcSPARSE at the binary level without modifying the source code of the upper-layer framework; Step S1: Obtain the sparse matrix.

[0015] S1.1 Reads the sparse matrix from host memory in CSR (Compressed Sparse Row) format. A ∈R^( M × K ) and vector x ∈R^ K .

[0016] S1.2 Allocate one-time video memory on the device side, record three arrays: row offset (rowPtr), column index (colIdx), and non-zero element value (val), and record the number of rows. M Column number K Non-zero element nnz Metadata such as COO or ELL. If the original format is COO or ELL, you can first call the mcSPARSE conversion API before entering S2.

[0017] Step S2: Divide the square matrix into sub-blocks.

[0018] S2.1 Based on the operand width requirement of the native matrix multiplication acceleration instruction mfma.f32.16x16x16.f16 of the Muxi GPU, divide A into several 16×16 sub-blocks Tile.

[0019] S2.2 Let each tile contain non-zero elements.NNZ_PER_TILE, Define the proportion of non-zero elements in a tile. η = ( NNZ_PER_ TILE ) / 256; If η If the percentage is less than 10%, then move all non-zero elements within that tile into the new tile. row_base_list And denote it as a CSR submatrix; otherwise, shift it into... block_base_list And denoted as the Tile-Dense submatrix.

[0020] After S2.3 is divided, block_base_list The data is stored in video memory as an array of structures, each containing 256 float16 values ​​and a 2-byte row offset, for subsequent direct mapping to the MXMACA instruction register layout. For other instruction specifications provided by the platform, such as mfma.f32.8x8x8.f16, the tile size can be changed to 8×8 accordingly, and the threshold can be adjusted accordingly. η To ensure that the computing power of the matrix multiplication unit is fully utilized.

[0021] Step S3: CSR sub-block row classification.

[0022] S3.1 Count the number of non-zero elements in each row NNZ( r ), r ∈[0, M -1], where r The line number.

[0023] S3.2 Set the threshold. th1 =64, which is the size of the thread bundle (64 threads). th2 =1024 means the maximum number of threads supported by a Thread BlockTile.

[0024] S3.3 Write rows into three types of arrays according to the following rules: Lightweight rows L1: NNZ( r ) ≤ th1 ;Medium row L2:th1 <NNZ( r ) ≤ th2 ; Weight line L3: NNZ ( r ) > th2 The three types of row arrays only store row numbers, occupying 4 bytes per row, and their video memory usage is negligible.

[0025] Step S4: Multiply CSR sub-blocks by vectors.

[0026] This step should be implemented in a three-stage manner within the same Muxi GPU kernel.

[0027] S4.1 Lightweight Row Computation. Each Thread Block Tile is responsible for computing multiple lightweight rows. The total number of non-zero elements in these lightweight rows should be less than the maximum number of threads supported by the Thread Block Tile, which is 1024. First, all threads within the Thread Block Tile cooperate to read the non-zero elements across rows to be computed into the programmable on-chip shared memory in a merged memory access manner. Second, within each thread bundle, the __shfl_xor_sync primitive performs register-level reduction. Finally, the result is written back. y_csr [ r ]. When NNZ( r When the number of threads is ≤32, smaller sub-warp (containing 2, 4, 8 or 16 threads) can be enabled to reduce the overhead of thread idle time.

[0028] S4.2 Medium-quantity line calculation. Each Thread Block Tile is responsible for one medium-quantity line, and the thread number within the Thread Block Tile is defined. tid ∈[0,1023]. First, the thread processes them sequentially. colIdx [ rowPtr [ r ]+ tid [Multiplication and addition; secondly, binary tree reduction is performed on the programmable shared memory, with a reduction depth of...] log 2(1024)=10; finally, by tid The thread with a value of 0 writes the reduction result to... y_csr [ r ].

[0029] S4.3 Weighted Row Calculation. A multi-threaded block tile collaborative calculation is used to calculate a weighted row. First, assume the row to be calculated requires... T Each Thread Block Tile participates in the specification and defines... T =NNZ( r ) / 1024; Secondly, after each Thread BlockTile completes its local reduction, it uses the int32_t counter C[ in global memory. r ]Execute atomicAdd(&C[ r ],1); Finally, when C[ r ]== T At that time, the last returning Thread Block Tile reads all local sums and completes the final reduction, then writes the result to... y_csr [ r ].

[0030] Step S5: Multiply vectors in Tile-Dense sub-blocks S5.1 starts the Muxi GPU kernel, with each Thread Block Tile responsible for a 16×16 Tile.

[0031] S5.2 Loads the 256 float16 values ​​in the tile into the register by calling the mfma_f32_16x16x16_f16 instruction. This completes the accumulation of A_tile × x, and the result is written back to the temporary array. y_tile .

[0032] S5.3 Repeat execution until block_base_list Traversal complete.

[0033] Step S6: Merge Results S6.1 starts the MuXi GPU kernel, and... y_tile and y_csr The final SpMV result is obtained by summing the row numbers. y .

[0034] S6.2 If reuse is still required in subsequent iterations A Then keep block_base_list , row_base_list The structure and preprocessing overhead can be amortized within 10-20 iterations. On average, net benefits are achieved after 12 iterations.

[0035] Performance verification: On the representative sparse matrix impcol_e, this embodiment achieves a speedup of 10.2 compared to mcsparseCsrmvEx. On the SuiteSparse sparse matrix dataset, compared to mcsparseCsrmvEx, this embodiment achieves an average speedup of 5.6, a peak speedup of 11.2, and a load balancing effect (the ratio of the maximum to minimum execution time of the Thread Block Tile) reduced from 12.3 to 1.3; the preprocessing time accounts for an average of 12% of the single SpMV time, and the preprocessing overhead can be amortized after 12 iterations.

[0036] The specific embodiments described herein are merely illustrative of the spirit of the invention. Those skilled in the art to which this invention pertains may make various modifications or additions to the described specific embodiments or use similar methods to substitute them, without departing from the spirit of the invention or exceeding the scope defined by the appended claims.

Claims

1. A sparse matrix-vector multiplication acceleration method for Muxi GPU, characterized in that, Includes the following steps: Step S1: Obtain the sparse matrix in Compressed Sparse Row (CSR) format. A ∈R^( M × K ) and vector x ∈R^ K And allocate video memory on the device side and store line offsets. rowPtr Column index colIdx Non-zero element values val Record the number of rows M Column number K Non-zero element nnz ; Step S2: Based on the operand width of the Muxi GPU matrix multiplication acceleration instruction mfma.f32.16×16×16.f16, multiply the sparse matrix... A Divide the area into 16×16 sub-tiles and define... NNZ_PER_TILE Non-zero elements within a Tile ， Calculate the proportion of non-zero elements for each tile. η = NNZ_PER_TILE / 256, of which; like η If the percentage is less than 10%, then move all non-zero elements within that tile into the new tile. row_base_list And denote it as a CSR submatrix, otherwise shift it in. block_base_list And denoted as the Tile-Dense submatrix; Step S3: Define r Given the row numbers of a sparse matrix, count the number of non-zero elements in each row: NNZ( r ), r ∈[0, M -1]; set up th1 =64, th2 =1024, classify rows according to the following rules: Lightweight rows L1: NNZ( r )≤ th1 ;Medium row L2: th1 <NNZ( r )≤ th2 ; Weight line L3: NNZ(r)> th2 And generate L1_rows, L2_rows, and L3_rows arrays that only store row numbers; Step S4: Execute the CSR submatrix in parallel within the same GPU kernel: For line L1, non-zero elements of the responsible line are read into on-chip shared memory via merged memory access, unit by unit, using the thread block tile primitive __shfl_xor_sync to perform register-level reduction, and the result is written to... y_csr [ r ]; For line L2, each Thread Block Tile is responsible for one line, defining the thread number within the Thread Block Tile. tid ∈[0,1023], perform binary tree reduction in shared memory, depth log 2(1024)=10, ultimately by tid =0 threads write y_csr [ r ]; Enable for L3 line T =NNZ( r 1024 Thread Block Tiles work together to complete the local reduction, with each Thread Block Tile using a global memory int32_t counter. C [ r ]Execute atomicAdd(& C [ r ],1), when C [ r ]= T The last Thread Block Tile completes the aggregation into a specification and writes it. y_csr [ r ]; Step S5: For the Tile-Dense submatrix, start the GPU kernel. Each Thread Block Tile is responsible for a 16×16 Tile. Load 256 float16 values ​​into the register, call the mfma_f32_16x16x16_f16 instruction to complete the local matrix-vector multiplication, and write the result to a temporary array. y_tile ; Step S6: Start the GPU kernel, and... y_tile and y_csr The final SpMV result is obtained by summing the row numbers. y。

Images