BLAS3 structured operator acceleration computing system based on hopper architecture gpu

Abstract

The application provides a BLAS3 structured operator acceleration calculation system based on a Hopper architecture GPU, and relates to the technical field of computers.The system comprises: a calculation unit discrimination module, which determines the calculation unit used by a current operator during operation and estimates the maximum row dimension upper bound of the tensor core execution path of the current operator; an instruction-aware block parameter determination module, which dynamically determines the optimal block size and quantity of an input matrix in real time; a block matrix loading and alignment module, which divides the input matrix and the matrix to be updated into submatrices based on the block size as a basic block and completes the loading of the corresponding submatrices; an operator kernel function execution module, which completes the shared memory structured parallel loading and storage of the double-precision floating point array of the submatrices corresponding to the input matrix, and calls the tensor core to perform multiplication and accumulation calculation; and a pipeline and concurrent scheduling module, which adds the block calculation task to the corresponding task set and performs multi-flow concurrent scheduling on the task set.

Description

Technical Field

[0001] This invention relates to the field of computer technology, and in particular to a BLAS3 structured operator-accelerated computing system based on a Hopper architecture GPU. Background Technology

[0002] Matrix computation is a core mathematical tool in scientific and engineering computing, widely used in complex computational tasks such as large-scale numerical simulations, optimization problems, and system modeling. Efficient matrix operations are crucial for improving computational performance and accuracy. To achieve this, the Basic Linear Algebra Subroutine (BLAS) serves as a fundamental computational kernel, widely applied in deep learning and high-performance computing. BLAS comprises three levels of operations, with the third level (BLAS3) primarily used for matrix multiplication and encompassing both structured operators (such as SYRK and SYR2K) and unstructured operators (such as GEMM).

[0003] With the development of modern GPU architecture, NVIDIA's Hopper architecture GPU, launched in 2022, introduced the fourth-generation Tensor Core, specifically designed to accelerate matrix multiplication and accumulation operations. This innovation significantly improved the execution efficiency of BLAS3-level operations. However, current mainstream high-performance computing libraries, such as cuBLAS and MAGMA, still have some performance bottlenecks.

[0004] Although the cuBLAS library optimizes for the Hopper architecture, its overall matrix processing mode lacks explicit block optimization for ultra-large-scale matrices, leading to tensor core degradation. Furthermore, the cuBLAS library's technical solution binds operator execution to a single CUDA stream, lacking support for multi-stream concurrent computation, thus limiting the full utilization of GPU resources (such as streaming multiprocessors (SM), memory bandwidth, etc.) in structured large-scale matrix computation. Actual tests show that the performance of structured operators in cuBLAS is only about 20% of the theoretical peak performance of GPUs, failing to fully utilize the parallel computing capabilities of GPUs, resulting in low resource utilization and significantly longer runtime.

[0005] Meanwhile, many numerical linear algebra computation routines in the MAGMA library (such as matrix decomposition and eigenvalue solving) rely on cuBLAS structured operators. When dealing with large-scale matrices, this is limited by the overall matrix processing mode of cuBLAS and fails to effectively utilize the tensor core, further restricting the potential of high-performance computing, becoming a bottleneck for performance improvement, and affecting the application and promotion of large-scale scientific and engineering problems. Summary of the Invention

[0006] The technical problem to be solved by the present invention is to provide a BLAS3 structured operator acceleration computing system based on Hopper architecture GPU to overcome the shortcomings of the prior art and realize the accelerated computing of BLAS3 structured operators.

[0007] To solve the above-mentioned technical problems, the technical solution adopted by the present invention is: a BLAS3 structured operator accelerated computing system based on Hopper architecture GPU, including a computing unit discrimination module, an instruction-aware block parameter determination module, a block matrix loading and alignment module, an operator kernel function execution module, and a pipeline and concurrency scheduling module.

[0008] The computational unit discrimination module dynamically collects GPU tensor core related metrics in real time in the target sampling section of the operator through a performance acquisition tool, determines the computational unit used by the current operator during runtime, and uses a progressive adaptive probe window growth strategy to estimate the upper bound of the maximum row dimension of the current operator in the tensor core execution path online.

[0009] The instruction-aware block parameter determination module dynamically senses the hardware resources for the tensor block size used by the tensor core computation instruction and the dimension information of the input matrix. Through an integer mapping collaborative block generation strategy based on tensor upper bound constraints, it dynamically determines the optimal block size parameter of the input matrix in real time. and number of blocks parameter The tensor upper bound constraint is the upper bound of the maximum row dimension of the operator in the tensor core execution path;

[0010] The block matrix loading and alignment module divides the input matrix and the matrix to be updated into sub-matrices with the block size as the basic block, and creates and allocates double-precision floating-point arrays to complete the loading of the corresponding sub-matrices;

[0011] The operator kernel function execution module is used on the GPU with tensor core aligned operation unit blocks as the basic granularity. Based on the tensor mapping description structure, it completes the shared memory structured parallel loading and storage of the double-precision floating-point array of the submatrix corresponding to the input matrix through direct memory access. It coordinates asynchronous transport and computation by using "arrival-wait" synchronization barriers. At the register level, it completes row or column priority alignment through thread program numbers and calls the tensor core through inline assembly instructions to perform matrix block-level multiplication and addition accumulation calculations in a thread group cooperative manner, and accumulates the results into register segments. Then, it performs cooperative write-back of the results in slices according to the structure region awareness rules.

[0012] The pipeline and concurrent scheduling module, based on the structured characteristics of different operators and the block-level task attributes and computational data dependencies, adds block computing tasks to the corresponding task sets, performs multi-stream concurrent scheduling on the task sets, and designs a CUDA stream dynamic optimization strategy based on wave capacity to dynamically and adaptively configure the number of CUDA streams at runtime to maximize the utilization of streaming multiprocessors.

[0013] Furthermore, the target sampling segment of the operator is obtained by marking "segment start" at the entry point of the computation library call and marking "segment end" at the exit point of the computation library call through a source-level instrumentation mechanism;

[0014] The GPU tensor core-related metrics include the average number of tensor computation instructions executed and the average number of active cycles in the tensor pipeline. The average number of tensor computation instructions executed represents the average number of thread-beam-level tensor computation instructions executed by the tensor core within the operator target sampling segment. The average number of active cycles in the tensor pipeline represents the average number of clock cycles in which the tensor computation pipeline is active within the operator target sampling segment.

[0015] Furthermore, the specific method for estimating the upper bound of the maximum row dimension of the current operator in the tensor core execution path online using the progressively adaptive probe window growth strategy is as follows:

[0016] Step S11: Using thread block segments as the basic unit, construct a probe window W of size P according to an integer multiple sequence of their size, which is used to generate the block row dimension sequence;

[0017] Step S12: In the target sampling segment of the operator, calculate the weighted joint score S of the average number of tensor computation instructions executed and the average number of active cycles of the tensor pipeline for each row dimension within the probe window W;

[0018] Step S13: Detect window size P The growth step size is adaptively increased until the first degradation signal is triggered. When the degradation signal χ appears for the first time, the row dimension of the previous valid block is taken as the upper bound of the maximum row dimension of the current operator in the execution path of the tensor core. The degradation signal χ is the signal that the computing unit is degraded from the tensor core to the CUDA core. The weighted joint score S of the average number of tensor computing instructions executed and the average number of active cycles in the tensor pipeline is less than or equal to the threshold as the degradation signal χ.

[0019] Furthermore, the optimal block size parameter is dynamically determined in real time through an integral mapping collaborative block generation strategy based on tensor upper bound constraints. and number of blocks parameter The specific method is as follows:

[0020] Step S21: During operator runtime, a performance analysis tool is used to perform real-time performance profiling of the tensor core execution process on the Hopper GPU architecture, dynamically extracting the tensor core aligned operation unit block size parameter used for tensor core matrix multiplication and addition instructions; the tensor core aligned operation unit block size is the size of the tensor block used to perform calculations in the tensor core calculation instructions; this size consists of three dimensions, represented as tile_m, tile_n, and tile_k, where: tile_m represents the size of the tensor block in the first dimension; tile_n represents the size of the tensor block in the second dimension; and tile_k represents the size of the tensor block in the third dimension;

[0021] Step S22: Initialize the value range of the block size nb to the least common multiple of tile_m and tile_n. The initial value set of nb consisting of integer multiples thereof. The maximum value in the set does not exceed the upper bound of the maximum row dimension that triggers the downgrading of the tensor core to a CUDA core;

[0022] Step S23: Iterate through the initial value set of block size nb obtained in step S22, and calculate the number of blocks based on the row dimension n of input matrix A and input matrix B. The value, and retain it in the set such that The final set of nb values ​​is composed of integer values. .

[0023] Furthermore, for the DSYRK, DSYR2K, and DSYRKX operators, when the matrix submatrix loading and alignment module divides the submatrix using the block size as the basic block size, it will load the double-precision random matrix... and Divide the data along the row dimension using a block size of nb as the basic block, and obtain... There are nb×k row submatrices; the matrix to be updated It is divided into two dimensions along both the row and column dimensions using a block size of nb as the basic block, forming There are nb×nb block matrices; where the matrix C to be updated is a double-precision symmetric matrix;

[0024] For DTRMM and DSYMM operators, when the block matrix loading and alignment module divides the submatrix using the block size as the basic block size, it will load the double-precision symmetric matrix... The data is divided into two dimensions along both the row and column dimensions, using blocks of size nb as the basic blocks, resulting in... A block matrix of size nb×nb; a double-precision random matrix and the matrix to be updated Divide the data along the row dimension using a block size of nb as the basic block, and obtain... There are nb×k row submatrices, where the column dimension k remains the same as the original matrix; among them, the matrix to be updated is... It is a double-precision random matrix.

[0025] Furthermore, the operator kernel function execution module includes a normal block multiplication and accumulation kernel function, a symmetric block update kernel function, and a symmetric block multiplication and accumulation kernel function;

[0026] The ordinary block multiplication accumulator kernel function will combine the double-precision random matrices A and B and the matrix to be updated C or a double-precision symmetric matrix. Double-precision random matrix and the matrix to be updated The sub-matrix blocks, divided according to a predetermined block partitioning method, are used as input. Matrix multiplication and accumulation operations are performed on each pair of sub-matrix blocks to finally generate a double-precision output matrix D or The corresponding submatrix;

[0027] The symmetric block update kernel function takes the double-precision random matrices A and B and the matrix to be updated C as input, divided into sub-matrix blocks according to a predetermined block partitioning method, and performs a symmetric matrix multiplication and accumulation operation with a rank not exceeding 2k on each pair of sub-matrix blocks, finally generating the corresponding sub-matrix of the double-precision output matrix D.

[0028] The symmetric block multiplication accumulator kernel function double-precision symmetric matrix Double-precision random matrix and the matrix to be updated The sub-matrix blocks, divided according to a predetermined block partitioning method, are used as input. Symmetric matrix multiplication and accumulation operations are performed on each pair of sub-matrix blocks to finally generate a double-precision output matrix. The corresponding submatrix;

[0029] Among them, the symmetric block update kernel function and the symmetric block multiply-accumulate kernel function perform thread block-level structure scheduling and pruning on the upper / lower triangular region of the matrix to be updated or the redundant symmetric region of the input matrix during the loading and writing-back phases of matrix multiplication operations.

[0030] Furthermore, the pipeline and concurrent scheduling module includes a task set construction submodule, an update rule dynamic construction submodule, a wave capacity-aware multi-stream concurrent adaptive configuration submodule, and a multi-stream concurrent scheduling submodule;

[0031] The task set construction submodule is used to execute a configurable block-level traversal and task construction strategy according to the operator type and structural characteristics, and add block computing tasks to the corresponding task set according to the strategy.

[0032] The update rule dynamic construction submodule optimizes the data loading, write-back dependencies, and task invocation relationships in the task set through the update rule dynamic construction mechanism. The update rule dynamic construction mechanism dynamically constructs different data update rules based on the data dependencies between different task sets and automatically adjusts the scheduling order of computing tasks according to these rules.

[0033] The wave capacity-aware multi-stream concurrent adaptive configuration submodule is used to obtain in real time the usage of each thread block of each kernel function of the operator kernel function execution module on the hardware resources of the streaming multiprocessor during the execution process, and determine the total wave capacity and workload according to the online decision mechanism of concurrent stream number driven by wave capacity calibration. It makes an online decision on the optimal number of concurrent streams of various CUDA streams under the current device and sends it to the multi-stream concurrent scheduling submodule to improve the utilization of streaming multiprocessors.

[0034] The multi-stream concurrent scheduling submodule adaptively selects the pipeline mode from the policy library containing dual-pipeline and triple-pipeline modes, and instantiates it into a CUDA stream based on the task set and operator structure features.

[0035] Furthermore, the configurable block-level traversal and task construction strategy includes a diagonal-order block-level traversal strategy for the upper triangular region of the matrix to be updated C, and an input matrix. Diagonal sequential block-level traversal strategy for the upper triangular region;

[0036] The diagonal sequential block-level traversal strategy for the upper triangular region of the matrix C to be updated is as follows: Set the diagonal-level iteration variable t∈[0, ... -1], perform block-level traversal and update the upper triangular region of the symmetric matrix C in diagonal order; for each iteration t, for all row block indices i∈ With column block index j∈[0, -1] Block matrix C satisfying ji=t ij When i=j, add the task of updating this submatrix to the task set. When i ≠ j, add the task of updating this submatrix to the task set. ;

[0037] The input matrix The diagonal sequential block-level traversal strategy for the upper triangular region is as follows: Set the diagonal-level iteration variable t∈[0, ... -1], for double-precision random symmetric matrices The upper triangular region is updated by block-level traversal in diagonal order; for each iteration t, for the matrix... All row block indices i∈[0, -1] and column block index j∈[0, -1] Block matrix satisfying ji=t When i=j, add the task of calculating this submatrix to the task set. When i ≠ j, add the task of calculating this submatrix to the task set. .

[0038] Furthermore, the dynamic construction mechanism for the update rules specifically involves: for the task set... and By constructing bilateral update rules for the block matrix to be updated Perform dynamic updates for the task set. and By constructing a one-sided update rule for the submatrix of rows to be updated Perform dynamic updates;

[0039] The bilateral update rule is calculated based on two types of matrices. and Joint block matrix The update is performed; during the update process, the calculation results of the two types of matrices are weighted by coefficients. and Merge into middle;

[0040] The unilateral update rule is based on single-class matrix calculation. Treatment of updating submatrix To perform an update, during the update process, only... The calculation results are merged into middle.

[0041] Furthermore, the multi-stream concurrent scheduling submodule, in dual-pipeline mode, divides the task set... The task set is fixed and assigned to the first type of pipeline to execute the symmetric block update kernel function. Assign to the second type of pipeline and batch schedule non-main diagonal regions according to diagonal numbering. In each scheduling, the positions of the submatrix contained in the current batch scheduling are paired. calculate , , , and The base address and layout mapping is used to call the ordinary block multiply-accumulate kernel function to complete the two-stage multiply-accumulate and write back;

[0042] In the three-pipeline mode, the first type of pipeline calls the symmetric block multiplication-accumulation kernel function to process the task set. The matrix product terms in the pipeline are calculated, and completion events are registered in their respective computation flows; the second type of pipeline calls the ordinary block multiplication accumulator kernel function to process the task set. The matrix product terms in the process are calculated and a completion event is registered in the computation flow to which they belong. After all the above completion events are satisfied, the third type of pipeline schedules the accumulator kernel, aggregates the matrix product results to the target block and completes the write-back.

[0043] In the pipeline instantiation into CUDA stream design, concurrent CUDA streams are dynamically acquired and configured through the wave capacity-aware multi-stream concurrent adaptive configuration module.

[0044] The beneficial effects of adopting the above technical solution are as follows: The BLAS3 structured operator accelerated computing system based on Hopper architecture GPU provided by this invention solves the key problem of tensor core degradation by dynamically and reasonably dividing the row load of large-size matrices through a block matrix algorithm based on the input matrix size and tensor core hardware awareness dynamic adjustment strategy. In the kernel function design, it fully combines the new features of the fourth-generation tensor core of Hopper architecture (Tensor Core natively supports double-precision floating-point (FP64) matrix multiply-accumulate instruction MMA and tensor memory accelerator) and the structural characteristics of the matrix (upper triangular matrix and symmetric matrix); by designing a CUDA stream dynamic configuration strategy based on wave capacity, it estimates the maximum number of concurrent CUDA streams on Hopper architecture GPU, and determines the number of concurrent startups of each type of CUDA stream at runtime, thereby maximizing the utilization of streaming multiprocessors. Attached Figure Description

[0045] Figure 1 A structural block diagram of the BLAS3 structured operator accelerated computing system based on Hopper architecture GPU provided in an embodiment of the present invention;

[0046] Figure 2 A schematic diagram of the submatrix structure of the DSYRK, DSYR2K, and DSYRKX operators provided in the embodiments of the present invention;

[0047] Figure 3 A schematic diagram of the submatrix structure of the DSYMM and DTRMM operators provided in the embodiments of the present invention;

[0048] Figure 4 A flowchart of a CUDA stream dynamic configuration strategy based on wave capacity provided in an embodiment of the present invention;

[0049] Figure 5 This is a schematic diagram of a dual pipeline mode provided in an embodiment of the present invention;

[0050] Figure 6 This is a schematic diagram of a three-pipeline mode provided in an embodiment of the present invention;

[0051] Figure 7The figure shows the experimental results of the speedup ratio between the accelerated computing system of the present invention and cuBLAS, provided for embodiments of the present invention.

[0052] Figure 8 The figure shows the experimental results of the speedup ratio between the accelerated computing system of the present invention and the MAGMA baseline implementation in the application of symmetric matrix eigenvalue decomposition, as provided in the embodiments of the present invention. Detailed Implementation

[0053] The specific embodiments of the present invention will be described in further detail below with reference to the accompanying drawings and examples. The following examples are for illustrative purposes only and are not intended to limit the scope of the invention.

[0054] In this embodiment, a BLAS3 structured operator-accelerated computing system based on a Hopper architecture GPU is used, such as... Figure 1 As shown, it includes a computational unit discrimination module, an instruction-aware block parameter determination module, a block matrix loading and alignment module, an operator kernel function execution module, and a pipeline and concurrency scheduling module;

[0055] The computational unit discrimination module dynamically collects GPU tensor core related metrics in real time in the target sampling section of the operator through a performance acquisition tool, determines the computational unit used by the current operator during runtime, and uses a progressive adaptive probe window growth strategy to estimate the upper bound of the maximum row dimension of the current operator in the tensor core execution path online.

[0056] The target sampling segment of the operator is obtained by marking "segment start" at the entry point of the computation library call and "segment end" at the exit point of the computation library call through a source-level instrumentation mechanism;

[0057] The GPU tensor core related metrics include the average number of tensor computation instructions executed and the average number of active cycles in the tensor pipeline. The average number of tensor computation instructions executed is used to characterize the average number of thread-beam-level tensor computation instructions executed by the tensor core within the operator target sampling segment. The average number of active cycles in the tensor pipeline is used to characterize the average number of clock cycles in which the tensor computation pipeline is active within the operator target sampling segment.

[0058] Tensor core execution path refers to the computational process and data flow path in a GPU where operators use tensor cores to perform matrix multiplication or other tensor operations. Specifically, tensor core execution path includes the entire process from data input, tensor block loading, computation execution to result output.

[0059] In a preferred embodiment, the two GPU tensor core-related metrics may correspond to the average execution count of tensor core computation instructions (DMMA-type instructions) and the average active cycle count of the pipeline in the GPU hardware performance counters, or equivalent metrics. These two metrics correspond to the NVIDIA NsightCompute metrics: sm__inst_executed_pipe_tensor_op_dmma (average number of executions of thread-beam tensor pipeline DMMA instructions) and sm__pipe_tensor_op_dmma_cycles_active (average number of cycles in which tensor pipeline DMMA instructions are active).

[0060] In this embodiment, the specific method for estimating the upper bound of the maximum row dimension of the current operator in the tensor core execution path online using a progressively adaptive probe window growth strategy is as follows:

[0061] Step S11: Using thread block segments as the basic unit, construct a probe window W of size P according to an integer multiple sequence of their size, which is used to generate the block row dimension sequence;

[0062] The thread block is divided into the smallest matrix sub-regions that the thread block is responsible for calculating and updating. The resulting block row dimension sequence is shown in the following formula:

[0063] ;

[0064] ;

[0065] in, Let be the candidate value for the first block row dimension in the detection window W, where the value represents the number of rows in the block along the row direction. For the detection window W, the first Candidate values ​​for each block row dimension Divide the thread block into blocks;

[0066] Step S12: In the target sampling region of the operator, for each row dimension within the detection window W The weighted joint score S of the average number of tensor computation instructions executed and the average number of active cycles in the tensor pipeline is calculated separately, as shown in the following formula:

[0067] ;

[0068] in, To calculate the average number of instructions executed for tensors, This represents the average number of active cycles in the tensor channel. and The preset weight is 0.5 by default.

[0069] Step S13: Detect window size P The growth step size is adaptively increased until the first degradation signal is triggered; when the degradation signal χ first appears, the row dimension of the previous valid block is taken as the upper bound of the maximum row dimension of the current operator in the tensor core execution path: = ; where the growth step size e represents the number of times the probe window slides forward; the degradation signal χ is the signal that the computing unit is downgraded from a tensor core to a CUDA core. The weighted joint score S of the average number of tensor computing instructions executed and the average number of active cycles in the tensor pipeline is less than or equal to a threshold as the degradation signal χ; the threshold is determined by the size of the input matrix and the type of tensor core computing instructions; the input matrix is ​​the data matrix entered by the user when calling the computing system, and this matrix is ​​processed as input data during the computing process;

[0070] The instruction-aware block parameter determination module dynamically senses the hardware resources and input matrix dimension information of the tensor core aligned operation unit block size used by the tensor core calculation instruction. Through the integer mapping collaborative block generation strategy of tensor upper bound constraints, it dynamically determines the optimal block size parameters of the input matrix in real time. and number of blocks parameter This allows each submatrix to be efficiently computed using Tensor Core computation instructions. The Tensor Core aligned operation unit block size is the size of the tensor block used to perform computations in Tensor Core computation instructions, particularly Matrix Multiply-Accumulate (MMA) instructions. This size consists of three dimensions, denoted as tile_m, tile_n, and tile_k, where: tile_m represents the size of the tensor block in the first dimension (row dimension); tile_n represents the size of the tensor block in the second dimension (column dimension); and tile_k represents the size of the tensor block in the third dimension (depth dimension, typically corresponding to the inner product dimension in matrix multiplication).

[0071] In this embodiment, the optimal block size parameter is dynamically determined in real time through an integer mapping collaborative block generation strategy based on tensor upper bound constraints. and number of blocks parameter The specific method is as follows:

[0072] Step S21: During operator runtime, use performance analysis tools to perform real-time performance analysis of the tensor core execution process on the HopperGPU architecture, and dynamically extract the tensor core aligned operation unit block size parameters used for tensor core matrix multiplication and addition instructions;

[0073] Step S22: Initialize the value range of the block size nb to the least common multiple of tile_m and tile_n. The initial value set of nb consisting of integer multiples thereof. The maximum value in the set does not exceed the upper bound of the maximum row dimension TC_Max_row that triggers the downgrading of the tensor core to a CUDA core;

[0074] , ;

[0075] in, This is the upper bound of the maximum row dimension of the current operator in the execution path of the tensor core. Let u be the least common multiple of tile_m and tile_n, and let u be a subset of the set of positive integers. A positive integer, This represents the set of positive integers, that is, all integers greater than zero.

[0076] This step ensures that each input matrix can be fully mapped to the matrix multiplication and addition instruction tensor block structure during computation, avoiding computational core degradation caused by tensor block boundaries and row dimensions exceeding the upper bound supported by the tensor kernel, thereby improving the utilization efficiency of the tensor kernel.

[0077] Step S23: Iterate through the initial value set of block size nb obtained in step S22, and calculate the number of blocks based on the row dimension n of input matrix A and input matrix B. The value, and retain it in the set such that The final set of nb values ​​is composed of integer values. ;

[0078] ;

[0079] This strategy can completely map an input matrix of any dimension into a submatrix, ensuring that access to and updating of the submatrix in subsequent loading and computation stages has structural integrity and memory safety.

[0080] The block matrix loading and alignment module divides the input matrix and the matrix to be updated into sub-matrices with the block size as the basic block, and creates and allocates double-precision floating-point arrays to complete the loading of the corresponding sub-matrices;

[0081] like Figure 2 As shown, for the DSYRK, DSYR2K, and DSYRKX operators: the block matrix loading and alignment module divides the matrix into sub-matrices based on the block size, which is to load the input matrix (double-precision random matrix). and Divide the data along the row dimension using a block size of nb as the basic block, and obtain... There are nb×k row submatrices; where the column dimension k remains the same as the original matrix; the matrix to be updated is a double-precision symmetric matrix. It is divided into two dimensions along both the row and column dimensions using a block size of nb as the basic block, forming There are nb×nb block matrices, and the operator kernel function execution module only needs to update the matrix. The upper triangular block submatrix is ​​processed.

[0082] The design of keeping the row submatrices of input matrices A and B unchanged from the column dimension k of the two input matrices is to accommodate the computational characteristics of the tensor cores in the Hopper architecture GPU. Since the Hopper architecture's tensor cores possess higher computational power and execution frequency, their computational capabilities are only fully utilized in algorithms with high arithmetic density. Therefore, to fully leverage their hardware computational potential, it is essential to ensure high arithmetic density in the operators, achieving a high ratio between computation and memory access. To improve the utilization of the tensor cores under the Hopper architecture, this structural design ensures that the row submatrices have a large column dimension. This allows each submatrix, with the block size as its basic unit, to maintain high arithmetic density in the column dimension during subsequent computations, thereby improving the parallel computation efficiency and throughput performance of the tensor cores, making it suitable for large-scale matrix computations on the Hopper architecture.

[0083] The operator kernel function execution module only processes the upper triangular block submatrix of the matrix C to be updated. This takes advantage of the structural characteristics of symmetric matrices and only processes the necessary upper triangular region to eliminate repeated calculations of symmetric positions in subsequent algorithms, thereby reducing computational redundancy and global memory access.

[0084] For the DSYRK, DSYR2K, and DSYRKX operators: the loading method for each submatrix is ​​as follows:

[0085] By setting block index variables i and j, the loading of different row submatrices of input matrix A and input matrix B is controlled, respectively. Simultaneously, the loading of block submatrices at different row and column address indices of the matrix C to be updated is also controlled, respectively. The ranges of i and j are [0, ..., ... -1], load the row submatrix A[i×nb:(i+1)×nb,:] of input matrix A into the corresponding double-precision floating-point array SA, and store it in column-major order; load the row submatrix B[j×nb:(j+1)×nb,:] of input matrix B into the corresponding double-precision floating-point array SB, and store it in column-major order; load the block submatrix C[i×nb:(i+1)×nb,j×nb:(j+1)×nb] of the matrix C to be updated that satisfies the preset condition i≤j into the corresponding double-precision floating-point array SC, and store it in column-major order;

[0086] like Figure 3 As shown, for the DTRMM and DSYMM operators, when the matrix submatrix loading and alignment module divides the submatrix using the block size as the basic block size, it will load the double-precision symmetric matrix... The data is divided into two dimensions along both the row and column dimensions, using blocks of size nb as the basic blocks, resulting in... A block matrix of size nb×nb; a double-precision random matrix and the matrix to be updated Divide the data along the row dimension using a block size of nb as the basic block, and obtain... There are nb×k row submatrices, where the column dimension k remains the same as the original matrix; among them, the matrix to be updated is... It is a double-precision random matrix.

[0087] Submatrices for DTRMM and DSYMM operators: The loading method for each submatrix is ​​as follows:

[0088] The double-precision symmetric matrix is ​​controlled by setting the block index variables i and j respectively. Different row address indices and different column address indices, i also controls the double-precision random matrix. and Different row address indices, i and j range from [0, ... -1]; By iterating through the row index i and column index j of the submatrix in a nested loop, in each iteration, if i ≤ j, then the outer loop variable i is used to iterate along the matrix. The content from row i×nb to row (i+1)×nb is extracted along the row index direction; and under the control of the inner loop variable j, the data from column j×nb to column (j+1)×nb is extracted along the column index direction, thus obtaining the matrix. Block matrix in The corresponding data area is And load the submatrix into a double-precision floating-point array. In the matrix, data is stored in column-major order for subsequent calculations. The loop variable `i` also controls the matrix. Loading of row submatrices, sequentially extracting double-precision random matrices. The submatrix block from row i×nb to row (i+1)×nb That is, the corresponding data area is And load it into a double-precision floating-point array. The data is stored in column-major order; double-precision random matrices are extracted sequentially, controlled by the loop variable i. The submatrix block from row i×nb to row (i+1)×nb That is, the corresponding data area is And load it into a double-precision floating-point array. In this context, data is stored in column-major order for use with submatrices. Participate in block matrix multiplication and addition operations.

[0089] The operator kernel function execution module is used on the Hopper GPU with tensor core aligned operation unit blocks as the basic granularity. Based on the tensor mapping description structure, it completes the structured parallel loading and shared memory storage of the double-precision floating-point arrays corresponding to the submatrices of the input matrix through direct memory access. It coordinates asynchronous transport and computation by using an "arrival-wait" synchronization barrier. At the register level, it completes row or column priority alignment through thread program numbers and calls the tensor core through inline assembly instructions to perform matrix block-level multiplication and addition accumulation calculations in a thread group cooperative manner, and accumulates the results into register segments. Then, it performs cooperative write-back of the results in slices according to the structure region awareness rules.

[0090] The operator kernel function execution module includes a normal block multiplication and accumulation kernel function, a symmetric block update kernel function, and a symmetric block multiplication and accumulation kernel function.

[0091] The ordinary block multiplication accumulator kernel function will combine the double-precision random matrices A and B and the matrix to be updated C or a double-precision symmetric matrix. Double-precision random matrix and the matrix to be updated The sub-matrix blocks, divided according to a predetermined block partitioning method, are used as input. Matrix multiplication and accumulation operations are performed on each pair of sub-matrix blocks to finally generate a double-precision output matrix D or The corresponding submatrix;

[0092] The symmetric block update kernel function takes the double-precision random matrices A and B and the matrix to be updated C as input, divided into sub-matrix blocks according to a predetermined block partitioning method. It performs a matrix multiplication accumulation operation with a symmetric rank not exceeding 2k on each pair of sub-matrix blocks, and finally generates the corresponding sub-matrix of the double-precision output matrix D. The matrix multiplication accumulation operation with a symmetric rank not exceeding 2k is to multiply the sub-matrix block of A with the transpose of the sub-matrix block of B, and accumulate the result of multiplying the sub-matrix block of B with the transpose of the sub-matrix block of A to generate a symmetric matrix block.

[0093] The symmetric block multiplication accumulator kernel function double-precision symmetric matrix Double-precision random matrix and the matrix to be updated The sub-matrix blocks, divided according to a predetermined block partitioning method, are used as input. Symmetric matrix multiplication and accumulation operations are performed on each pair of sub-matrix blocks to finally generate a double-precision output matrix. The corresponding submatrix;

[0094] Among them, the symmetric block update kernel function and the symmetric block multiply-accumulate kernel function perform thread block-level structure scheduling and pruning on the lower triangle (or redundant symmetric region) during the loading and write-back phases to reduce redundant computation and memory access, and improve throughput and computation unit occupancy stability.

[0095] The ordinary block multiplication accumulator kernel function includes the following calculation steps:

[0096] Step C1: Construct and configure the tensor mapping description structure of the submatrices corresponding to input matrix A and matrix B. This structure describes the logical block division of the submatrices corresponding to input matrix A and matrix B in global memory and the data transfer path between them and shared memory, using tensor core aligned operation unit blocks as the unit. This supports the structured and parallel data loading operation implemented by the tensor memory access engine in the Hopper architecture in step C3.

[0097] The tensor mapping description structure describes the block partitioning method of each submatrix in the global memory's two-dimensional tensor core alignment operation unit blocks, and determines the mapping relationship with shared memory through preset step size, logical size and alignment strategy; wherein, the tensor mapping description structure is configured to be directly read through tensor transport instructions, and realize the structured loading of block matrix-level data without relying on thread index calculation.

[0098] Step C2: Initialize the shared memory array and synchronization barrier mechanism of the tensor core alignment operation unit block corresponding to the cached submatrix within the thread block, which is used to cache the tensor core alignment operation unit block data and ensure data consistency between threads;

[0099] The synchronization barrier mechanism includes one or more barrier objects used to coordinate the synchronization dependency between the loading of tensor core alignment operation unit blocks and subsequent computation processes. The mechanism employs arrive-wait control semantics, allowing each thread to non-block the recording of data copy states without obstructing the execution path, and to release the data uniformly after the transfer is complete, ensuring the integrity and access consistency of the tensor core alignment operation unit block data in the shared memory. Two synchronization barrier objects are initialized within each thread block, respectively controlling the access order of the tensor core alignment operation unit blocks of the submatrices of input matrices A and B. The two synchronization barrier objects are configured by thread 0 calling the `init()` function, and other threads synchronize the transfer completion state through `arrive()` and `wait()` operations.

[0100] Step C3: All tensor core alignment operation unit blocks are cyclically loaded along the row dimension of the submatrix of input matrix A and the column dimension of the submatrix of input matrix B. The data transfer of the tensor core alignment operation unit blocks is initiated by a single thread in the thread block. This thread obtains the mapping information of the current tensor core alignment operation unit block according to the preset tensor mapping description structure, and calls the asynchronous copy instruction to move the operation unit block from global memory to the shared memory cache through the tensor memory access engine using direct memory access. During data transfer, all threads register their arrival status through the synchronization barrier mechanism, and after the transfer is completed, they release the signal through the barrier object and uniformly proceed to step C4.

[0101] Step C4: Each thread within the thread block loads the tensor core aligned operation unit blocks corresponding to the submatrices of input matrices A and B from shared memory to registers using an index-based method, based on its thread procedure number identifier value within the thread bundle. To ensure consistency with the alignment format of tensor core calculation instructions in subsequent steps, the tensor core aligned operation unit blocks of the submatrices of matrix A are loaded into registers in row-major order, while those of the submatrices of matrix B are loaded into registers in column-major order. Through thread bundle mapping control and register array scheduling, tensor core aligned operation unit blocks are explicitly bound to thread registers, improving the efficiency of concurrent data loading at the thread bundle level.

[0102] Step C5: Use inline PTX assembly instructions to perform matrix multiplication and addition operations on the supported tensor dimensions using the double-precision tensor core matrix multiplication and addition instructions under the Hopper architecture. During execution, the tensor core aligned operation unit blocks corresponding to the submatrices of input matrices A and B come from the register fragments in step C4, and the accumulated results are stored in the register group c_frag corresponding to the tensor core aligned operation unit block of the block submatrix of the matrix C to be updated.

[0103] Step C6: Calculate the position of the data in the register set c_frag obtained in step C5 within the submatrix of the matrix C to be updated, based on the thread index, and write it back to global memory. The write-back process uses a thread-based cooperative approach, updating the corresponding operand matrix sub-segments in the submatrix of the matrix C to be updated in row-major order. To address the possibility of non-full matrix fragments, thread-level conditional checks ensure that only valid data is written back, avoiding out-of-bounds access and guaranteeing memory access safety.

[0104] The first three processing steps of the symmetric block update kernel function are the same as the processing steps C1-C3 of the ordinary block multiplication and accumulation kernel function. However, there are design differences in the subsequent steps: since the submatrices of the matrix C to be updated in the symmetric block update kernel function are symmetric, the kernel function design makes full use of symmetry to optimize GPU computing.

[0105] In step C4, thread block-level structure scheduling is first pruned: the grid block dimension configuration in the kernel function startup primitive is used to schedule only thread blocks that meet the condition that the row index of the tensor core alignment operation unit block is less than or equal to the column index of the tensor core alignment operation unit block, thereby avoiding the calculation and resource occupation of the lower triangular region of the submatrix of matrix C;

[0106] In step C6, a structure region awareness judgment is performed: before each write-back, a judgment condition is added: whether the row index coordinate of the current write-back position is less than or equal to the column index coordinate of the current write-back position. Only if this condition is met will a write-back be performed, ensuring that only the upper triangular region of the submatrix of the matrix to be updated C is written back. Simultaneously, a boundary judgment condition is combined: the row index coordinate of the current write-back position is less than the row dimension of the submatrix of the matrix to be updated C, and the column index coordinate of the current write-back position is less than the column dimension of the submatrix of the matrix to be updated C, ensuring that address access is legal. By sensing whether the write-back position is within the upper triangular region within the thread, the submatrix of the matrix to be updated C is written back according to the upper triangular structure, precisely avoiding invalid access to the lower triangular region.

[0107] The first three processing steps of the symmetric block multiplication-accumulation kernel function are the same as steps C1-C3 of the ordinary block multiplication-accumulation kernel function. However, there are design differences in the subsequent steps: because the input matrix of the symmetric block multiplication-accumulation kernel function... The submatrices of the kernel have symmetry, so the symmetry can be fully utilized in the kernel function design to optimize GPU computing.

[0108] In step C4, each thread calculates an index to access the matrix based on its position identifier (lane_id) within the thread bundle. and The corresponding submatrix is ​​loaded from shared memory into a register, for the matrix To ensure the symmetry of the submatrix, a preset condition is added: for thread blocks where the row index of the tensor core aligned operation unit block is less than or equal to the column index of the tensor core aligned operation unit block, only matrices satisfying the preset condition are loaded. The part of the submatrix.

[0109] In this embodiment, the calculation process of the ordinary block multiplication and accumulation kernel function is implemented on the GPU as follows:

[0110] Step 1: Use the CUtensorMap data structure provided by the NVIDIA CUDA Driver API to construct a tensor map description structure. This structure describes the logical block distribution of matrices in global memory and the data transfer path between global and shared memory, using tensor core aligned operation unit blocks as units. The CUtensorMap data structure includes the following fields:

[0111] baseAddress: Represents the starting address of the input tensor in global video memory;

[0112] globalShape: Defines the overall logical size of the tensor, i.e. the number of elements in each dimension;

[0113] globalStride: Represents the storage stride of a tensor in each dimension (in bytes), used to support non-contiguous memory access modes;

[0114] boxShape: Represents the logical shape of a single tensor core aligned operation unit block in each dimension;

[0115] boxStride: Represents the stride of the elements in each dimension within the tensor core alignment operation unit block. It is usually set to 1, indicating that the elements are stored contiguously in the tensor core alignment operation unit block.

[0116] Taking tensor core alignment operation unit block sizes tile_m=16, tile_k=16, tile_n=8 as an example, when the submatrix dimension of the input matrix A is nb×k, it can be divided into × There are n×k tensor core alignment operation unit blocks, where each tensor core alignment operation unit block corresponds to a 16×16 submatrix region; the submatrix of matrix B (with dimensions nb×k) can be divided into: × Tensor core alignment operation unit blocks with size tile_k=16 and tile_n=8 are used. For each tensor core alignment operation unit block, its mapping description structure is encoded into a GPU-recognizable format through the CUtensorMapEncode interface and copied to the device for subsequent asynchronous data transfer by the tensor memory access engine.

[0117] The construction of the tensor mapping structure described above enables high-throughput shared memory loading operations to be completed in units of tensor core aligned operation units under the Hopper architecture without explicit thread index management. This effectively improves the parallelism and memory access efficiency of tensor computation, reduces global memory access latency, and improves the overall execution performance of the operation.

[0118] Multiple shared memory arrays are allocated within the thread block to cache input matrix fragments, storing data from the corresponding tensor core aligned operation unit blocks in the submatrices of matrix A and matrix B. The shared memory arrays are constructed according to the size of the tensor core aligned operation unit blocks and are configured to meet the alignment boundaries required by the tensor memory access engine. Alignment constraints are implemented through memory alignment instructions to ensure that subsequent tensor core aligned operation unit block data can be directly loaded using a single direct memory access transfer. Preferably, the arrays are configured to meet 128-byte boundary alignment to match the hardware handling requirements of the GPU tensor memory access engine.

[0119] Step 2: Within each thread block, declare two shared memory arrays sA[n×k] and sB[n×k] for temporarily storing submatrix data; these correspond to the submatrix data of input matrix A and input matrix B, respectively. The arrays are aligned to 128 bytes using the alignas(128) instruction to meet the shared memory address alignment requirements of the tensor memory access engine during direct memory access operations.

[0120] Two thread-block-level barrier objects, `barA` and `barB`, are declared to control the asynchronous loading process of tensor core alignment operation unit blocks. These barriers control the copying progress of tensor core alignment operation unit blocks for submatrices of input matrices A and B, respectively. The declarations of the two barrier objects are as follows: `__shared__barrierbarA`; `__shared__barrierbarB`. Initialization operations are performed only by the thread with thread ID 0 within the thread block: `init(&barA, blockDim.x)` binds the barrier objects `barA` and `barB` to the current thread block, allowing all threads to record and synchronize. The `fence_proxy_async_shared_cta()` instruction forces synchronization of data visibility between temporary memory operations and shared memory, ensuring that data is promptly visible to computation threads after asynchronous data transfer. Finally, the thread-block-level synchronization instruction `__syncthreads()` is called once to unify the execution progress of all threads, ensuring that the shared memory space, barrier objects, and register states are initialized before entering the tensor core alignment operation unit block data transfer and computation phase, maintaining consistent visibility across threads.

[0121] Step 3: During the loading of tensor core alignment operation unit block data, the only thread in the thread block calls the tensor memory access engine based on the preset tensor mapping description structure and initiates an asynchronous transport operation on the current tensor core alignment operation unit block.

[0122] The data transfer process is executed using the GPU hardware's tensor memory access engine. It employs direct memory access to copy tensor core-aligned computation unit block data from global memory to the shared memory cache. During this process, other threads register their synchronization states non-blockingly using a synchronization barrier mechanism. Upon completion of the data copy, they receive a barrier release signal to trigger subsequent thread-based collaborative computation operations. Two `arrival_token` objects, `tokenA` and `tokenB`, are declared to track the current state of the tensor core-aligned computation unit block loading barrier. The token serves as a handle to the barrier controller's state and will be used later for precise blocking and releasing of synchronization operations. Then, all tensor core-aligned computation unit block submatrices are iteratively loaded along the middle dimension direction of matrix multiplication. The middle dimension index of the tensor core-aligned computation unit block is `block_k_iter`, and the total number of iterations is `num_blocks_k = k / tile_k`. Each iteration corresponds to loading the tensor core aligned operation unit block submatrices of input matrices A and B into shared memory. In each iteration, an asynchronous tensor transport instruction is initiated by the only thread within the thread block (threadIdx.x==0), using the tensor memory access engine to call the cp_async_bulk_tensor_2d_global_to_shared() interface to transport the data of the specified tensor core aligned operation unit block from global memory to shared memory addresses shared_A and shared_B. Here, tensorMapA and tensorMapB are tensor mapping description structures generated on the host side and copied to the device. The transport operation precisely locates the required data fragment according to the middle dimension index (block_k_iter) and the row and column indices (block_m_iter, block_n_iter) of the tensor core aligned operation unit block. After each transport by the Tensor Memory Access Engine, the arrival state of the corresponding barrier object is automatically updated, and a token for subsequent wait() operations is returned. Other threads within the thread block do not participate in the transport initiation during this iteration; they only call barA.arrive() and barB.arrive() to register the barrier state, ensuring consistent semantic behavior for subsequent synchronization operations. After the transport initiation or arrival registration, all threads uniformly call the barX.wait(tokenX) function to wait for the barrier to be released, ensuring that the Tensor Memory Access Engine operation has been completed and the data in the tensor core aligned operation unit block in the shared memory can be safely read. After the barrier is released, the thread synchronization primitive __syncthreads() is called once to implement a synchronization barrier for all threads in the thread block, ensuring that the shared memory is viewed as a barrier. Figure 1Therefore, the tensor core aligned operation unit block data can be used by each thread in subsequent multiply-accumulate calculations;

[0123] Step 4: The process of loading the corresponding tensor core aligned operation unit blocks from the submatrices of matrices A and B into the registers based on the thread index adopts the following loading strategy:

[0124] Each thread determines its coordinates (a_row, a_col) and (b_row, b_col) in the sub-fragment of the operand matrix by two-dimensional expansion based on its thread bundle position identifier lane_id. The sub-matrix of input matrix A is stored in row-major order, and the intra-thread mapping is: a_row = lane_id / 4, a_col = lane_id%4; the sub-matrix of input matrix B is stored in column-major order, and the intra-thread mapping is: b_row = lane_id%4, b_col = lane_id / 4. Based on this, each thread is responsible for loading multiple sets of tensor core aligned operation unit blocks from global memory and filling their register arrays a_frag[8] and b_frag[4] in sequence.

[0125] Submatrix loading strategy for matrix A: Map register array a_frag[0-7] to a 16×1 column fragment (spanning 16 rows, width 1), covering 16 row dimensions (tile_m=16) in the tensor core alignment operation unit block data input; By traversing variable l∈[0,7], the upper and lower row positions are distinguished by even and odd indices: When l is even, the value of register array comes from the upper half row region of the current tensor core alignment operation unit block, and the address offset is expressed as: A[global_row+a_row×K+(l×2)+tensor core alignment operation unit block_k+a_col]; where global_row represents the global row index, a_row represents the row index in the submatrix of matrix A, K is the dimension constant of the matrix, tensor core alignment operation unit block_k represents the alignment block position in the tensor core, and a_col is the column index of the submatrix of matrix A. When l is odd, the value of the register array comes from the lower half of the tensor core aligned operation unit block (i.e., shifted down 8 rows relative to a_row), and the address offset is expressed as: A[global_row+(a_row+8)×K+(l×2-2)+tensor core aligned operation unit block_k+a_col]; This design logically constructs a two-layer row expansion mechanism, loading the 16 row dimensions of the submatrix corresponding to matrix A into the register array in two complete segments, satisfying the requirement in mma.sync.m16n8k16 that the operands of matrix A are in row-major order and can be accessed across rows.

[0126] The loading strategy for the submatrix corresponding to matrix B is as follows: each thread loads 4 fragments into b_frag[0-3], which are arranged in a column-first layout of the B operands in the Tensor Core Alignment Operation Unit Block data. The loading address is expressed as: B[Tensor Core Alignment Operation Unit Block_k+b_row+(l×4)+global_col+b_col×K]; this address mapping realizes the sampling of the 16×8 submatrix of B, where the increment of l is 1 in each iteration, that is, each iteration steps 4 rows. Finally, each thread extracts its corresponding k-dimensional fragments from the 16 rows of B in column order. This layout is aligned with the column-major order format of B required by mma.sync, ensuring that the Tensor Core matrix multiplication and addition instructions accept continuous and structured matrix column fragments. This register loading mapping strategy maps the thread position identifier lane_id to two-dimensional coordinates in the global matrix. By using an interleaved row and column offset strategy to load the operand matrix sub-segments of A and B in segments, it achieves precise structural docking with the tensor core computation mode and provides data layout guarantee for the efficient execution of mma.sync.aligned.m16n8k16.row.col.f64.f64.f64.f64.f64.

[0127] In this embodiment, firstly, based on the lane_id value of the current thread within the thread bundle, the thread index mapping mechanism within the thread bundle is used: (lane_id / 4) controls the row offset, (lane_id%4) controls the column offset, to calculate the sub-elements of matrix A and matrix B that each thread needs to load, and store them in a_frag[8] and b_frag[4] respectively. To maximize register utilization efficiency, the following optimization strategy is adopted: the 16×16 tensor core aligned operation unit block data block of matrix A is mapped to a continuous register group, and the 8 double-precision elements required by the thread bundle are precisely loaded through 8 loops; the 16×8 tensor core aligned operation unit block of matrix B only uses 4 threads to complete 4 data extractions, and each time it is loaded into the b_frag register array;

[0128] Step 5: After loading, the double-precision tensor core calculation instruction mma.sync.aligned.m16n8k16.row.col.f64.f64.f64.f64.rn supported by the Hopper architecture is called. The register fragments a_frag and b_frag are used as the input of the instruction, respectively. The high-throughput thread bundle-level matrix multiplication and accumulation calculation is performed, and the initial result fragment is stored in the register group c_frag[4]. This implementation has the following characteristics: each execution can complete D=A×B+C, where A is 16×16, B is 16×8, and C / D is 16×8. The inline PTX explicitly controls the binding of the source register and the target register, thereby providing maximum flexibility for thread calculation within the thread bundle. During the calculation process, the tensor core pipeline is automatically used to realize hardware-level fused multiply-accumulate (FMA) operation, which significantly improves the computation efficiency within a unit cycle.

[0129] Step Six: Calculate the corresponding position of the data in register group c_frag obtained in Step Five in the global matrix C based on the thread index, and write it back to global memory. The write-back process adopts a thread-based cooperative approach, updating the corresponding fragments in the block sub-matrix C in row-major order. To address the possibility of non-full fragments in the matrix, thread-level conditional checks ensure that only valid data is written back, avoiding out-of-bounds access and guaranteeing memory access safety.

[0130] Calculate the corresponding position of the thread index in the global matrix C using the following mapping strategy: Each thread writes the calculation results in c_frag[0 - 3] into the global matrix C in a two-dimensional expansion manner according to its position identifier lane_id. The thread implements a two-dimensional mapping of the sub-slice results of the operand matrix through (row = lane_id / 4, col = lane_id % 4), where each thread corresponds to the position corresponding to the row and column (row, col) of a tensor core alignment operation unit block with a row size of 16 and a column size of 8 (tile_m = 16, tile_n = 8). Further, c_frag[0 - 3] is logically divided into upper and lower sub-matrices: for the index variable q satisfying the preset condition q < 2, it is written into the upper half row region of the corresponding global matrix C; for the index variable q satisfying the preset condition q ≥ 2, it is written into the lower half row region (i.e., offset 8 rows downward relative to global_row); this structure forms a logically two-layer nested write-back path, effectively mapping the output fragment results of the tensor core alignment operation unit block with the shape of tile_m × tile_n = 16 × 8 to the corresponding positions in the global matrix; ensuring conflict-free and continuous write-back, and adapting to the parallel execution model of the tensor core of the Hopper architecture. After the calculation is completed, the thread determines the updated position in the target matrix C according to its (lane_id / 4) and (lane_id % 4) index information. Each thread writes the calculation results in c_frag[0 - 3] back to C[row × N + col] in a two-dimensional expansion manner, where: through the two-layer nested write-back logic, the sub-slice mapping of the result matrix with tile_m = 16 and tile_n = 8 to the exact memory address of the global matrix is realized: C[(global_row + row) × N + global_col + col × 2 + (q & 0x1)] += c_frag[q] (q < 2) and C[(global_row + (row + 8)) × N + global_col + col × 2 + (q & 0x1)] += c_frag[q] (q ≥ 2); global_row and global_col represent the row and column indices in the global matrix C, which are calculated according to the alignment strategy of the tensor core, specifically: int global_row = tensor core alignment operation unit block_row × tensor core alignment operation unit block_M; int global_col = tensor core alignment operation unit block_col × tensor core alignment operation unit block_N; row and col are calculated through lane_id respectively: int row = lane_id / 4; int col = lane_id % 4; and by adding boundary judgment (row < M && col < N), the generality of the large matrix block calculation is enhanced and the security of memory access is ensured;

[0131] The pipeline and concurrent scheduling module is designed to address the structured characteristics of different operators. Based on block-level task attributes and computational data dependencies, it performs multi-stream concurrent scheduling of task sets. Furthermore, it employs a CUDA stream dynamic optimization strategy based on wave capacity to dynamically and adaptively configure the number of CUDA streams at runtime, maximizing the utilization of the streaming multiprocessor. This module includes the following sub-modules:

[0132] 1. Task Set Construction Submodule: This module executes configurable block-level traversal and task construction strategies based on operator type and structural characteristics, and adds block computation tasks to the corresponding task set according to these strategies.

[0133] The configurable block-level traversal and task construction strategy includes a block-level traversal strategy for the diagonal order of the upper triangular region of the update matrix C, and an input matrix. Diagonal sequential block-level traversal strategy for the upper triangular region;

[0134] The configurability is reflected in the following: For the DSYRK, DSRKX, and DSYR2K operators, the strategy is configured as: a block-level traversal strategy for the diagonal order of the upper triangular region of the update matrix C; for the DSYMM and DTRMM operators, the strategy is configured as: input matrix Diagonal sequential block-level traversal strategy for the upper triangular region;

[0135] The diagonal sequential block-level traversal strategy for the upper triangular region of the update matrix C is as follows: Set the diagonal-level iteration variable t∈[0, ... -1], perform block-level traversal and update the upper triangular region of the symmetric matrix C in diagonal order; for each iteration t, for all row block indices i∈[0, -1] and column block index j∈[0, -1] Block matrix C satisfying ji=t ij When i=j, add the task of updating this submatrix to the task set. When i ≠ j, add the task of updating this submatrix to the task set. ;

[0136] The input matrix The diagonal sequential block-level traversal strategy for the upper triangular region is as follows: Set the diagonal-level iteration variable t∈[0, ... -1], for double-precision random symmetric matrices The upper triangular region is updated by block-level traversal in diagonal order; for each iteration t, for the matrix... All row block indices i∈[0, -1] and column block index j∈[0, -1] Block matrix satisfying ji=t When i=j, add the task of calculating this submatrix to the task set. When i ≠ j, add the task of calculating this submatrix to the task set. .

[0137] The task set It has the following characteristics: block matrix Only need to be with row submatrix B i Perform calculations; the task set It has the following characteristics: block matrix Requires the row submatrix B i The calculation also requires the row submatrix B. j Perform calculations;

[0138] 2. Dynamic Rule Construction Submodule: This submodule optimizes data loading, write-back dependencies, and task invocation relationships within task sets through a dynamic rule construction mechanism. This mechanism dynamically constructs different data update rules based on data dependencies between different task sets and automatically adjusts the scheduling order of computation tasks according to these rules. In this way, the system ensures coordination and data consistency between tasks and effectively enables parallel task execution.

[0139] The dynamic construction mechanism for the update rules specifically involves: targeting the task set... and By constructing bilateral update rules, the block matrix C to be updated is treated as follows: ij Perform dynamic updates for the task set. and By constructing unilateral update rules, the block matrix to be updated is... Perform dynamic updates;

[0140] The bilateral update rule is calculated based on two types of matrices. and block matrix The update is performed. During the update process, the calculation results of the two types of matrices are weighted by coefficients. and merge into Specifically, for each block matrix update task that satisfies ji=t, the update rule is as follows:

[0141] ;

[0142] in, This is the updated block matrix, i.e., the output matrix;

[0143] The bilateral update rule clarifies two types of matrix multiplication ( and The operations do not have a specific order and can be computed in parallel. Furthermore, for the same submatrix to be updated... The results of multiple matrix multiplications are summed. and This approach has a clear sequence, concatenating addition operations to establish edge relationships, ensuring the correct writing order, and effectively handling tasks that update the same output block multiple times. The explicit management of computation order through the concatenated edge structure of the dependency graph avoids data conflicts in parallel computing, thereby improving computational efficiency and ensuring correct task execution. In matrix update tasks, the bilateral update rule allows for parallel computation of two matrix operations, and the concatenated edges of the dependency graph control the order of addition operations, reducing data conflicts and computational bottlenecks caused by parallel computing in traditional methods.

[0144] The unilateral update rule is based on single-class matrix calculation. Treatment of updating submatrix Perform an update. During the update process, only [the process involves...]. The calculation results; specifically, for those that meet the conditions For each block matrix update task, the update rule is as follows:

[0145] ;

[0146] 3. Wavelength Capacity Aware Multi-Stream Concurrency Adaptive Configuration Submodule: This module is used to obtain the real-time usage of hardware resources (registers, shared memory, etc.) of the streaming multiprocessor by each thread block during the execution of each kernel function of the operator kernel function execution module. Based on the online decision-making mechanism for the number of concurrent streams driven by the wavelength capacity calibration, it determines the total wavelength capacity and workload, and makes an online decision on the optimal number of concurrent streams for various CUDA streams under the current device. The decision is then sent to the multi-stream concurrency scheduling submodule to improve the utilization of the streaming multiprocessor.

[0147] like Figure 4 As shown, the specific method for determining the total wave capacity and workload based on the concurrent flow online decision-making mechanism driven by wave capacity calibration is as follows:

[0148] Step A1: Calculate the single streaming multiprocessor wave capacity of the symmetric block update kernel function, the ordinary block multiply-accumulate kernel function, and the symmetric block multiply-accumulate kernel function based on registers, shared memory, and thread constraints, and obtain the total wave capacity.

[0149] The single-stream multiprocessor wave capacity of the symmetric block update kernel function is defined as: in the Hopper architecture, the maximum number of thread blocks that the symmetric block update kernel function can execute concurrently on a single stream multiprocessor, that is, the upper limit of the number of thread bundles that can reside at any time, which depends on the usage of stream multiprocessor hardware resources (registers, shared memory, etc.) by each thread block of the symmetric block update kernel function.

[0150] The formula for calculating the single streaming multiprocessor wave capacity of the symmetric block update kernel function is as follows:

[0151] ;

[0152] in, For updating the kernel function of a symmetric block, the single streaming multiprocessor wave capacity, The maximum number of registers available for each streaming multiprocessor (SM). The maximum shared memory available for each streaming multiprocessor. The maximum number of threads available for each streaming multiprocessor (SM). The number of registers required to update a single thread block of the kernel function for a symmetric block. The shared memory required by a single thread block to update the kernel function for a symmetric block; these parameters , , , , The performance was obtained in real time during the execution of the symmetric block update kernel function in the operator kernel function execution module by using the NVIDIA NsightCompute performance analysis tool (ncu).

[0153] The single-stream multiprocessor wave capacity of the ordinary block multiply-accumulate kernel function is defined as follows: In the Hopper architecture, the maximum number of thread blocks that the ordinary block multiply-accumulate kernel function can execute concurrently on a single stream multiprocessor depends on the usage of stream multiprocessor hardware resources (registers, shared memory, etc.) by each thread block of the ordinary block multiply-accumulate kernel function.

[0154] The formula for calculating the single-stream multiprocessor wave capacity of the ordinary block multiply-accumulate kernel function is as follows:

[0155] ;

[0156] in, The single-stream multiprocessor wave capacity of the ordinary block multiply-accumulate kernel function. The number of registers required by a single thread block of the accumulator kernel function for a normal block. The shared memory required by a single thread block for multiplying a normal block by the accumulator kernel function; these parameters , The performance was obtained in real time during the execution of the ordinary block multiplication and accumulation kernel function in the operator kernel function execution module using the NVIDIA NsightCompute performance analysis tool (ncu).

[0157] The single-stream multiprocessor wave capacity of the symmetric block multiply-accumulate kernel function is defined as follows: In the Hopper architecture, the maximum number of thread blocks that the symmetric block multiply-accumulate kernel function can execute concurrently on a single stream multiprocessor depends on the usage of stream multiprocessor hardware resources (registers, shared memory, etc.) by each thread block of the symmetric block multiply-accumulate kernel function.

[0158] The formula for calculating the single streaming multiprocessor wave capacity of the symmetric block multiply-accumulate kernel function is as follows:

[0159] ;

[0160] in, Single streaming multiprocessor wave capacity of symmetric block multiply-accumulate kernel function The number of registers required for a single thread block of the symmetric block multiplication accumulator kernel function. The shared memory required for a single thread block of the symmetric block multiplication accumulator kernel function; these parameters , The NsightCompute performance analysis tool (ncu) was used to obtain real-time data during the execution of the symmetric block multiplication-accumulation kernel function in the operator kernel function execution module;

[0161] Step A2: Calculate the workload of the symmetric block update kernel function, the ordinary block multiply-accumulate kernel function, and the symmetric block multiply-accumulate kernel function calls, and determine the proportion of streaming multiprocessors and the total wave capacity allocated to each type of CUDA stream based on the pipeline mode selected by the multi-stream concurrent scheduling submodule according to the workload.

[0162] The workload calculation formula for the symmetric update kernel function call is as follows: ;

[0163] The workload calculation formula for the ordinary block multiplication accumulator kernel function call is as follows: ;

[0164] The workload calculation formula for the symmetric block multiply-accumulate kernel function call is as follows: ;

[0165] For dual-flow mode:

[0166] Assigned to the set of tasks to be executed The number of streaming multiprocessors for the first type of CUDA stream of the block matrix update task for:

[0167] ;

[0168] in, This refers to the total number of streaming multiprocessors (SMs) in a Hopper architecture GPU. The quota ratio corresponding to the first type of CUDA stream in the dual-pipeline mode is determined by the following formula:

[0169] ;

[0170] Assigned to the set of tasks to be executed The number of streaming multiprocessors for the second type of CUDA stream of the block matrix update task for:

[0171] ;

[0172] in, The quota ratio corresponding to the second type of CUDA stream in the dual-pipeline mode is determined by the following formula:

[0173] ;

[0174] Assigned to the set of tasks to be executed The first type of CUDA stream total wave capacity of the block matrix update task for:

[0175] ;

[0176] Assigned to the set of tasks to be executed Total wave capacity of the second type of CUDA stream in the block matrix update task for:

[0177] ;

[0178] For the three-flow model:

[0179] Assigned to the task set responsible for execution The number of streaming multiprocessors for the first type of CUDA stream in the block matrix update task for:

[0180] ;

[0181] in, This refers to the total number of streaming multiprocessors (SMs) in a Hopper architecture GPU. The quota ratio corresponding to the first type of CUDA flow under the three-stream model is determined by the following formula:

[0182] ;

[0183] Assigned to the set of tasks to be executed Multi-stream processor count of the second type of CUDA stream for block matrix update tasks for:

[0184] ;

[0185] in, The quota ratio corresponding to the second type of CUDA flow under the three-flow model is determined by the following formula:

[0186] ;

[0187] Assigned to the set of tasks to be executed The first type of CUDA stream total wave capacity of the block matrix update task for:

[0188] ;

[0189] Assigned to the set of tasks to be executed Total wave capacity of the second type of CUDA stream in the block matrix update task for:

[0190] ;

[0191] Step A3: Based on the pipeline mode selected by the multi-stream concurrent scheduling submodule, calculate the number of thread blocks required to start and execute a single CUDA stream for each type of task set under the current pipeline mode, and make an online decision on the number of CUDA streams to be started concurrently for each type of task set according to the allocated total wave capacity.

[0192] For dual pipeline mode:

[0193] Single execution task set Number of thread blocks required to start and execute CUDA streams for:

[0194] ;

[0195] in, This indicates the number of threads in the row direction of the thread block that updates the kernel function of the symmetric block. This indicates the number of threads in the column direction of the thread block that updates the kernel function of the symmetric block;

[0196] Single execution task set Number of thread blocks required to start and execute CUDA streams for:

[0197] ;

[0198] in, This indicates the number of threads in the line direction of the thread block in the ordinary block multiplication and accumulation kernel function. This indicates the number of threads in the column direction of the thread block of the ordinary block multiply-accumulate kernel function;

[0199] Execution task set Number of concurrent CUDA stream launches for: = ;

[0200] Execution task set Number of concurrent CUDA stream launches for: = ;

[0201] For the three-pipeline model:

[0202] Single execution task set Number of thread blocks required to start and execute CUDA streams for:

[0203] ;

[0204] in, This indicates the number of threads in the row direction of the symmetric block multiplication-accumulation kernel function thread block. This indicates the number of threads in the column direction of the symmetric block multiplication-accumulation kernel function thread block;

[0205] Single set of computing tasks Number of thread blocks required to start and execute CUDA streams for:

[0206] ;

[0207] in, This indicates the number of threads in the line direction of the thread block in the ordinary block multiplication and accumulation kernel function. This indicates the number of threads in the column direction of the thread block of the ordinary block multiply-accumulate kernel function;

[0208] Execution task set Number of concurrent CUDA stream launches for: = ;

[0209] Execution task set Number of concurrent CUDA stream launches for: = ;

[0210] 4. Multi-stream Concurrent Scheduling Submodule: From a strategy library containing dual-pipeline and triple-pipeline modes, this module adaptively selects the pipeline mode based on task set and operator structure characteristics and instantiates it into a CUDA stream. In dual-pipeline mode, task set ρ1 is fixedly assigned to the first type of pipeline to execute the symmetric block update kernel function, and task set ρ2 is assigned to the second type of pipeline. Non-main diagonal regions are batch-scheduled according to diagonal numbering. In each scheduling, A is calculated based on the position sequence pair (i,j) of the block submatrix contained in the current batch. i B j and C i,j The base address and layout mapping is used to call the ordinary block multiply-accumulate kernel function to complete the two-stage multiply-accumulate and write back;

[0211] In the three-pipeline mode, the first and second types of pipelines respectively call the symmetric block multiply-accumulate kernel function and the ordinary block multiply-accumulate kernel function to compute the task set. , The two types of matrix product terms are registered and completed within their respective streams. The third type of convergence stream schedules the accumulator kernel to converge to the target block and write it back after all events are satisfied. In the pipeline instantiation into CUDA stream design, concurrent CUDA streams are dynamically acquired and configured through the wave capacity-aware multi-stream concurrent adaptive configuration module.

[0212] In a preferred embodiment, the dual-pipeline mode is as described in the appendix to the specification. Figure 5 As shown, two types of CUDA streams are used to implement the execution of each pipeline. The first type of CUDA stream, StreamClass1, is responsible for executing the task set. The task of updating a block matrix corresponds to updating the block matrix C at the main diagonal position where i=j. ij The specific design is as follows: First, read the row submatrix A. i The corresponding double-precision floating-point array SA, read the row submatrix B. j The corresponding double-precision floating-point array SB is read and its transpose is performed; the block submatrix C of the matrix C to be updated is read. ij Once the corresponding double-precision floating-point array SC is completed, the symmetric block update kernel function is called to execute the block submatrix C. ij The update operation is performed, and the calculation result is written back to the block submatrix D of the output matrix D. ij The corresponding double-precision floating-point array SD. The second type of CUDA stream, StreamClass2, is responsible for executing the task set. Block matrix update tasks, which correspond to updating matrices C such that ji=t, where t∈[1, The block matrix C at the non-main diagonal positions of [-1] ijTo improve kernel function throughput and data concurrency loading efficiency to adapt to the high-throughput computation of Hopper architecture TensorCore, a batch processing approach is adopted, using the diagonal as the unit: In each CUDA stream scheduling, all submatrix update tasks that satisfy the current diagonal number t are collected into a task batch, and all block submatrices of the current diagonal are updated at once using a batch processing method. Specifically, according to the preset condition ji=t, where t∈[1, -1], first execute the block submatrix C ij of Update operation: For each value of t, calculate the memory address mapping relationship from the matrix block to the batch: row submatrix A i The starting memory address in global memory is A_base + i × nb, and the memory layout is nb × k. A_base represents the starting address of matrix A in global memory; row submatrix B j The starting address of the global memory is B_base + j × nb, and the memory layout is nb × k. B_base is the starting address of matrix B in global memory; the block submatrix C ij The starting address in global memory is C_base + i × nb × n + j × nb, and the memory layout is nb × nb. C_base is the starting address of matrix C in global memory; read the row submatrix A. i The corresponding double-precision floating-point arrays SA_1 and B j The corresponding double-precision floating-point array SB_1 is transposed, and the block submatrix C is read. ij The corresponding double-precision floating-point array SC_1, after being loaded and transposed, calls the ordinary block multiplication and accumulation kernel function to execute the block submatrix C. ij The update operation has a batch size of [size missing]. The calculation results are then written back to the submatrix D of the output matrix D. ij The corresponding double-precision floating-point array SD_1; then execute the block submatrix C. ij of Update operation: This operation uses a batch processing method similar to... The same. Through the above pipeline design, two different CUDA streams are responsible for the computation of different task types of the matrix C to be updated. This ensures that all tasks without data dependencies can be executed concurrently, achieving maximum overlap and pipeline between data loading, computation, and write-back, and maximizing the concurrent execution capability and resource utilization efficiency of the GPU.

[0213] Three-line model such as Figure 6 As shown, the first type of CUDA stream, StreamClass1, is responsible for executing the task set. The block matrix update task is specifically designed as follows: (based on the number of blocks) (For example) First, read the block matrix. The double-precision floating-point array SA_21 corresponding to (i=j=0) is used to read the row submatrix. The corresponding double-precision floating-point array SB_21; after loading, the symmetric block multiplication-accumulation kernel function is called to execute the row submatrix of the matrix C to be updated. The update operation is performed, and the calculation result is written back to the row submatrix of the output matrix D. The corresponding double-precision floating-point array SD_21; then read the block matrix. The double-precision floating-point array SA_22 corresponding to (i=j=1) is used to read the row submatrix. The corresponding double-precision floating-point array SB_22; after loading, the symmetric block multiplication-accumulation kernel function is called to execute the row submatrix of the matrix C to be updated. The update operation is performed, and the calculation result is written back to the row submatrix of the output matrix D. The corresponding double-precision floating-point array SD_22.

[0214] The second type of CUDA stream, StreamClass2, is responsible for executing task sets. To fully utilize the symmetric matrix structure of matrix A, the task of updating the block matrix involves loading the block matrix once. (i≠j) to the double-precision floating-point array SA_211 respectively call the ordinary block multiplication and accumulation kernel function to calculate the row submatrix. and The product values ​​are written back to the double-precision floating-point arrays SD_211 and SD_222 respectively. This implementation reduces the number of accesses to the double-precision floating-point array SA_211.

[0215] For each of the two types of CUDA streams corresponding to a product term, after the corresponding kernel function computation is completed in the stream, a computation completion event `event_k` is registered within the stream by calling the CUDA-provided event recording primitive API `cudaEventRecord(event_k, StreamClass_k)`. This event is used to indicate that the corresponding intermediate product submatrix has been computed and is in a state that can be read and accumulated.

[0216] The third type of CUDA stream, StreamClass3, is responsible for loading the row submatrices of the matrix C to be updated. The double-precision floating-point array SC_21 corresponding to (i=0) is used, and the matrix calculation results SD_21 and SD_211 of the first and second type CUDA streams are accumulated and written back to the double-precision array SD_21 corresponding to the row submatrix of the output matrix D; the row submatrix of the matrix C to be updated is loaded. The double-precision floating-point array SC_22 corresponding to (i=1); and the matrix calculation results of the first type of CUDA stream and the second type of CUDA stream are accumulated by SD_22 and SD_222, and written back to the double-precision array SD_22 corresponding to the row submatrix of the corresponding output matrix D;

[0217] The specific design is as follows: A third type of CUDA stream, StreamClass3, is created specifically for performing the final accumulation operation. Within StreamClass3, cudaStreamWaitEvent(StreamClass_k, event_k) is called sequentially to explicitly establish synchronization dependencies on the event_k in each StreamClass_k. This mechanism ensures that tasks in StreamClass3 can only start after all intermediate product calculations are complete, thus avoiding write conflicts and inconsistencies caused by concurrent data access. The synchronization dependency is satisfied when all intermediate product blocks tempC are completed. k The data is ready; the accumulation kernel function is scheduled and started in CUDA StreamClass3. This accumulation kernel function is responsible for loading the submatrix of the matrix to be updated. The original value, and all intermediate product blocks temp temp ,...,temp And perform the following accumulation operation:

[0218] ;

[0219] The output submatrix obtained after accumulation The value will be written back to the target storage area in global memory through this accumulation kernel function. Since all data dependencies have been synchronized via an event mechanism, the write-back operation is performed sequentially after all product calculations, ensuring data consistency and accumulation correctness.

[0220] In this embodiment, the BLAS3 structured operator acceleration computing system based on the Hopper architecture GPU of the present invention was tested on a Hopper architecture GPU platform (NVIDIA H100-SMX) using a large-scale double-precision random matrix (the number of rows in the matrix ranges from 49152 to 65536, while the number of columns is fixed at 1024 and 2048). The experimental results of the acceleration effect are shown in Table 1, with the runtime comparisons showing an acceleration of [missing information]. Figure 7As shown, compared with the existing CUDA Basic Linear Algebra Subroutines (cuBLAS), the accelerated computing system of this invention achieves a speedup of up to 14.99 times. From the performance experimental results, as shown in Table 2, the highest performance of the structured operator using the existing technology (cuBLAS) is 13.292 TFLOPs, only 20% of the theoretical peak performance of H100 (67 TFLOPs), while the highest performance of the accelerated computing system of this invention is 48.885 TFLOPs, 73% of the theoretical peak performance of H100. It is evident that the accelerated computing system of this invention achieves significant improvements in both acceleration effect and performance, meeting the application requirements for high-performance computing in the fields of artificial intelligence and scientific engineering computing.

[0221] To further verify the beneficial effects of the accelerated computing system of this invention in numerical linear algebra routines in the field of scientific and engineering computing, this embodiment compares the performance and runtime of the baseline implementation routine in the MAGMA library and the implementation routine of this invention in the baseline implementation with the underlying BLAS3 operator replaced by the baseline implementation routine in the same hardware and parameter settings, specifically in the applications of symmetric matrix tridiagonalization and symmetric matrix eigenvalue decomposition. Performance data is shown in Table 3. In the application of symmetric matrix tridiagonalization, the accelerated computing system of this invention achieved a maximum performance of 831.36 GFLOPs, significantly better than the 764.56 GFLOPs of the baseline implementation of the MAGMA technology. Runtime data is shown in Table 4. For example, the accelerated computing system... Figure 8 As shown, in the application of symmetric matrix eigenvalue decomposition, the accelerated computing system of the present invention achieves a 1.2-fold speedup compared to the MAGMA technology solution.

[0222] Table 1 Comparison of the running time of the accelerated computing system of the present invention and cuBLAS

[0223]

[0224] Table 2 Performance Comparison of the Accelerated Computing System of the Present Invention with cuBLAS

[0225]

[0226] Table 3 Performance comparison of the accelerated computing system of the present invention and the MAGMA baseline implementation in tridiagonalization applications.

[0227]

[0228] Table 4 Comparison of runtime between the accelerated computing system of this invention and the MAGMA baseline implementation in symmetric matrix eigenvalue decomposition applications.

[0229] 。

Claims

1. A BLAS3 structured operator-accelerated computing system based on Hopper architecture GPU, characterized in that: It includes a computational unit discrimination module, an instruction-aware block parameter determination module, a block matrix loading and alignment module, an operator kernel function execution module, and a pipeline and concurrency scheduling module; The computational unit discrimination module dynamically collects GPU tensor core related metrics in real time in the target sampling section of the operator through a performance acquisition tool, determines the computational unit used by the current operator during runtime, and uses a progressive adaptive probe window growth strategy to estimate the upper bound of the maximum row dimension of the current operator in the tensor core execution path online. The instruction-aware block parameter determination module dynamically senses the hardware resources for the tensor block size used by the tensor core computation instruction and the dimension information of the input matrix. Through an integer mapping collaborative block generation strategy based on tensor upper bound constraints, it dynamically determines the optimal block size parameter of the input matrix in real time. and number of blocks parameter The tensor upper bound constraint is the upper bound of the maximum row dimension of the operator in the tensor core execution path; The block matrix loading and alignment module divides the input matrix and the matrix to be updated into sub-matrices with the block size as the basic block, and creates and allocates double-precision floating-point arrays to complete the loading of the corresponding sub-matrices; The operator kernel function execution module is used on the GPU with tensor core aligned operation unit blocks as the basic granularity. Based on the tensor mapping description structure, it completes the shared memory structured parallel loading and storage of the double-precision floating-point arrays of the submatrices corresponding to the input matrix through direct memory access. It coordinates asynchronous transport and computation by using "arrival-wait" synchronization barriers. At the register level, it completes row or column priority alignment through thread program numbers and calls the tensor core through inline assembly instructions to perform matrix block-level multiplication and addition accumulation calculations in a thread group cooperative manner, and accumulates the results into register segments. Then, it performs cooperative write-back of the results in slices according to the structure region awareness rules. The pipeline and concurrent scheduling module, based on the structured characteristics of different operators and the block-level task attributes and computational data dependencies, adds block computing tasks to the corresponding task sets, performs multi-stream concurrent scheduling on the task sets, and designs a CUDA stream dynamic optimization strategy based on wave capacity to dynamically and adaptively configure the number of CUDA streams at runtime to maximize the utilization of streaming multiprocessors.

2. The BLAS3 structured operator accelerated computing system based on Hopper architecture GPU according to claim 1, characterized in that: The target sampling segment of the operator is obtained by marking "segment start" at the entry point of the computation library call and marking "segment end" at the exit point of the computation library call through a source-level instrumentation mechanism; The GPU tensor core-related metrics include the average number of tensor computation instructions executed and the average number of active cycles in the tensor pipeline. The average number of tensor computation instructions executed represents the average number of thread-beam-level tensor computation instructions executed by the tensor core within the operator target sampling segment. The average number of active cycles in the tensor pipeline represents the average number of clock cycles in which the tensor computation pipeline is active within the operator target sampling segment.

3. The BLAS3 structured operator accelerated computing system based on Hopper architecture GPU according to claim 2, characterized in that: The specific method for estimating the upper bound of the maximum row dimension of the current operator in the tensor core execution path online using the progressively adaptive probe window growth strategy is as follows: Step S11: Using thread block segments as the basic unit, construct a probe window W of size P according to an integer multiple sequence of their size, which is used to generate the block row dimension sequence; Step S12: In the target sampling segment of the operator, calculate the weighted joint score S of the average number of tensor computation instructions executed and the average number of active cycles of the tensor pipeline for each row dimension within the probe window W; Step S13: Detect window size P The growth step size is adaptively increased until the first degradation signal is triggered. When the degradation signal χ first appears, the row dimension of the previous valid block is taken as the upper bound of the maximum row dimension of the current operator in the execution path of the tensor core. The degradation signal χ is the signal that the computing unit is degraded from the tensor core to the CUDA core. The weighted joint score S of the average number of tensor computing instructions executed and the average number of active cycles of the tensor pipeline is less than or equal to the threshold as the degradation signal χ.

4. The BLAS3 structured operator accelerated computing system based on Hopper architecture GPU according to claim 3, characterized in that: The optimal block size parameter is dynamically determined in real time through an integral mapping collaborative block generation strategy based on tensor upper bound constraints. and number of blocks parameter The specific method is as follows: Step S21: During operator runtime, a performance analysis tool is used to perform real-time performance profiling of the tensor core execution process on the Hopper GPU architecture, dynamically extracting the tensor core aligned operation unit block size parameter used for tensor core matrix multiplication and addition instructions; the tensor core aligned operation unit block size is the size of the tensor block used to perform calculations in the tensor core calculation instructions; this size consists of three dimensions, represented as tile_m, tile_n, and tile_k, where: tile_m represents the size of the tensor block in the first dimension; tile_n represents the size of the tensor block in the second dimension; and tile_k represents the size of the tensor block in the third dimension; Step S22: Initialize the value range of the block size nb to the least common multiple of tile_m and tile_n. The initial value set of nb consisting of integer multiples thereof. The maximum value in the set does not exceed the upper bound of the maximum row dimension that triggers the downgrading of the tensor core to a CUDA core; Step S23: Iterate through the initial value set of block size nb obtained in step S22, and calculate the number of blocks based on the row dimension n of input matrix A and input matrix B. The value, and retain it in the set such that The final set of nb values ​​is composed of integer values. .

5. The BLAS3 structured operator accelerated computing system based on Hopper architecture GPU according to claim 4, characterized in that: For the DSYRK, DSYR2K, and DSYRKX operators, when the matrix submatrix loading and alignment module divides the submatrix using the block size as the basic block size, it will load the double-precision random matrix... and Divide the data along the row dimension using a block size of nb as the basic block, and obtain... There are nb×k row submatrices; the matrix to be updated It is divided into two dimensions along both the row and column dimensions using a block size of nb as the basic block, forming There are nb×nb block matrices; where the matrix C to be updated is a double-precision symmetric matrix; For DTRMM and DSYMM operators, when the block matrix loading and alignment module divides the submatrix using the block size as the basic block size, it will load the double-precision symmetric matrix... The data is divided into two dimensions along both the row and column dimensions, using blocks of size nb as the basic blocks, resulting in... A block matrix of size nb×nb; a double-precision random matrix and the matrix to be updated Divide the data along the row dimension using a block size of nb as the basic block, and obtain... There are nb×k row submatrices, where the column dimension k remains the same as the original matrix; among them, the matrix to be updated is... It is a double-precision random matrix.

6. The BLAS3 structured operator accelerated computing system based on Hopper architecture GPU according to claim 5, characterized in that: The operator kernel function execution module includes a normal block multiplication and accumulation kernel function, a symmetric block update kernel function, and a symmetric block multiplication and accumulation kernel function. The ordinary block multiplication accumulator kernel function will combine the double-precision random matrices A and B and the matrix to be updated C or a double-precision symmetric matrix. Double-precision random matrix and the matrix to be updated The sub-matrix blocks, divided according to a predetermined block partitioning method, are used as input. Matrix multiplication and accumulation operations are performed on each pair of sub-matrix blocks to finally generate a double-precision output matrix D or The corresponding submatrix; The symmetric block update kernel function takes the double-precision random matrices A and B and the matrix to be updated C as input, divided into sub-matrix blocks according to a predetermined block partitioning method, and performs a symmetric matrix multiplication and accumulation operation with a rank not exceeding 2k on each pair of sub-matrix blocks, finally generating the corresponding sub-matrix of the double-precision output matrix D. The symmetric block multiplication accumulator kernel function double-precision symmetric matrix Double-precision random matrix and the matrix to be updated The sub-matrix blocks, divided according to a predetermined block partitioning method, are used as input. Symmetric matrix multiplication and accumulation operations are performed on each pair of sub-matrix blocks to finally generate a double-precision output matrix. The corresponding submatrix; Among them, the symmetric block update kernel function and the symmetric block multiply-accumulate kernel function perform thread block-level structure scheduling and pruning on the upper / lower triangular region of the matrix to be updated or the redundant symmetric region of the input matrix during the loading and writing-back phases of matrix multiplication operations.

7. The BLAS3 structured operator accelerated computing system based on Hopper architecture GPU according to claim 6, characterized in that: The pipeline and concurrent scheduling module includes a task set construction submodule, an update rule dynamic construction submodule, a wave capacity-aware multi-stream concurrent adaptive configuration submodule, and a multi-stream concurrent scheduling submodule; The task set construction submodule is used to execute a configurable block-level traversal and task construction strategy according to the operator type and structural characteristics, and add block computing tasks to the corresponding task set according to the strategy. The update rule dynamic construction submodule optimizes the data loading, write-back dependencies, and task invocation relationships in the task set through the update rule dynamic construction mechanism. The update rule dynamic construction mechanism dynamically constructs different data update rules based on the data dependencies between different task sets and automatically adjusts the scheduling order of computing tasks according to these rules. The wave capacity-aware multi-stream concurrent adaptive configuration submodule is used to obtain in real time the usage of each thread block of each kernel function of the operator kernel function execution module on the hardware resources of the streaming multiprocessor during the execution process, and determine the total wave capacity and workload according to the online decision mechanism of concurrent stream number driven by wave capacity calibration. It makes an online decision on the optimal number of concurrent streams of various CUDA streams under the current device and sends it to the multi-stream concurrent scheduling submodule to improve the utilization of streaming multiprocessors. The multi-stream concurrent scheduling submodule adaptively selects the pipeline mode from the policy library containing dual-pipeline and triple-pipeline modes, and instantiates it into a CUDA stream based on the task set and operator structure features.

8. The BLAS3 structured operator accelerated computing system based on Hopper architecture GPU according to claim 7, characterized in that: The configurable block-level traversal and task construction strategy includes a diagonal-order block-level traversal strategy for the upper triangular region of the matrix to be updated C, and an input matrix. Diagonal sequential block-level traversal strategy for the upper triangular region; The diagonal sequential block-level traversal strategy for the upper triangular region of the matrix C to be updated is as follows: Set the diagonal-level iteration variable t∈[0, ... -1], perform block-level traversal and update the upper triangular region of the symmetric matrix C in diagonal order; for each iteration t, for all row block indices i∈ With column block index j∈[0, Block matrices satisfying ji=t When i=j, add the task of updating this submatrix to the task set. When i ≠ j, add the task of updating this submatrix to the task set. ; The input matrix The diagonal sequential block-level traversal strategy for the upper triangular region is as follows: Set the diagonal-level iteration variable t∈[0, ... -1], for double-precision random symmetric matrices The upper triangular region is updated by block-level traversal in diagonal order; for each iteration t, for the matrix... All row block indices i∈[0, ] and column block index j∈[0, Block matrices satisfying ji=t When i=j, add the task of calculating this submatrix to the task set. ; When i≠j, add the task of calculating this submatrix to the task set. .

9. The BLAS3 structured operator accelerated computing system based on Hopper architecture GPU according to claim 8, characterized in that: The dynamic construction mechanism for the update rules specifically involves: targeting the task set... and By constructing bilateral update rules for the block matrix to be updated Perform dynamic updates for the task set. and By constructing a one-sided update rule for the submatrix of rows to be updated Perform dynamic updates; The bilateral update rule is calculated based on two types of matrices. and Joint block matrix The update is performed; during the update process, the calculation results of the two types of matrices are weighted by coefficients. and Merge into middle; The unilateral update rule is based on single-class matrix calculation. Treatment of updating submatrix To perform an update, during the update process, only... The calculation results are merged into middle.

10. The BLAS3 structured operator accelerated computing system based on Hopper architecture GPU according to claim 9, characterized in that: The multi-stream concurrent scheduling submodule, in dual-pipeline mode, will schedule the task set... The task set is fixed and assigned to the first type of pipeline to execute the symmetric block update kernel function. Assign to the second type of pipeline and batch schedule non-main diagonal regions according to diagonal numbering. In each scheduling, the positions of the submatrix contained in the current batch scheduling are paired. calculate , , , and The base address and layout mapping is used to call the ordinary block multiply-accumulate kernel function to complete the two-stage multiply-accumulate and write back; In the three-pipeline mode, the first type of pipeline calls the symmetric block multiplication-accumulation kernel function to process the task set. The matrix product terms in the pipeline are calculated, and completion events are registered in their respective computation flows; the second type of pipeline calls the ordinary block multiplication accumulator kernel function to process the task set. The matrix product terms in the process are calculated and a completion event is registered in the computation flow to which they belong. After all the above completion events are satisfied, the third type of pipeline schedules the accumulator kernel, aggregates the matrix product results to the target block and completes the write-back. In the pipeline instantiation into CUDA stream design, concurrent CUDA streams are dynamically acquired and configured through the wave capacity-aware multi-stream concurrent adaptive configuration module.

Images