Processor-oriented dense matrix multiplication rank 2k symmetric updating method

Abstract

The invention discloses a processor-oriented dense matrix multiplication rank 2k symmetric updating method. The method comprises the following steps: 1) for a calculation task Cout = alpha A * BT + alpha B * AT + beta Cin, initializing the dimension A of a matrix to be N * K, the dimension B of the matrix to be K * N, and the dimensions of an output matrix Cout and an input matrix Cin to be N * N at a host end; copying the matrixes A, B and Cin from the host side memory to the equipment side memory; (2) the processor of the equipment side splits the calculation task into an alpha A * BT + alpha B * AT subtask and a beta * Cin subtask; the method comprises the following steps: firstly, calculating beta * Cin, and storing a calculation result into a memory of a device end; then, dividing an alpha A * BT + alpha B * AT subtask into batch triangular matrix multiplication of diagonal blocks and recursive batch rectangular matrix multiplication of non-diagonal blocks; 3) carrying out parallel calculation on staggered batch diagonal sub-blocks along the line; and 4) performing block recursion to call batch matrix multiplication to calculate a non-diagonal block.

Description

TECHNICAL FIELD

[0001] The present application belongs to the technical field of computer science, and particularly relates to a 3-level function of a BLAS mathematical library of a domestic GPU platform, and discloses a high-performance implementation method of rank-2k symmetric update and conjugate symmetric rank-2k update of dense matrix multiplication for a shuguang 7000 processor, which comprises an interleaved batch computation parallel method based on a fine-grained optimization algorithm and automatic code generation, and a calculation task division and data mapping parameter automatic optimization framework according to a structure characteristic of a domestic GPU. BACKGROUND

[0002] The terms related to the present application are explained as follows:

[0003] 1) BLAS (Basic Linear Algebra Subprograms): basic linear algebra operation mathematical library;

[0004] 2) CPU (Central Processing Unit): central processing unit;

[0005] 3) GPU (Graphics Processing Unit): graphics computing unit;

[0006] 4) GEMM (General Matrix Multiplications): general matrix multiplication, the calculation formula is C = alpha * op (A) * op (B) + beta * C, wherein op (X) is op (X) = X or op (X) = X T or op (X) = X H ;

[0007] 5) SYRK (Symmetric Rank-k Update): symmetric matrix rank-k update, the calculation formula is C = alpha * A * A T + beta * C or C = alpha * A T * A + beta * C;

[0008] 6) SYR2K (Symmetric Rank-2k Update): symmetric matrix rank-2k update, the calculation formula is C = alpha * A * B T + alpha * B * A T + beta * C or C = alpha * A T * B + alpha * B T * A + beta * C;

[0009] 7) HER2K (Hermitian Rank-2k Update): conjugate symmetric matrix rank-2k update, the calculation formula is C = alpha * A * B H + conjg (alpha) * B * A H+βC or C=αA H ×B+conjg(α)B H ×A+βC;

[0010] 8) MAGMA (Matrix Algebra on GPU and Multicore Architectures): An open-source, high-performance extended mathematics library that supports heterogeneous platforms, including GPU platforms;

[0011] 9) LDS (Local Data Share) Local shared memory.

[0012] 10) Tensile: A matrix multiplication automatic code generation tool suitable for domestic GPU platforms.

[0013] BLAS is one of the most important extended mathematical library modules in the field of high-performance computing. It is widely used in high-performance computing, artificial intelligence, scientific and engineering applications, providing reliable computational support and significant performance improvements for upper-level applications. In particular, benchmark programs used for ranking world supercomputers, training artificial neural networks, and symmetric eigenvalue solvers all heavily rely on the performance of BLAS level 3 functions. Therefore, research on BLAS level 3 functions is of significant importance and practical value. The BLAS level 3 functions SYR2K and HER2K studied in this patent take a symmetric matrix or conjugate symmetric matrix as output and only calculates the rank-2K update by multiplying the lower or upper triangular parts of the output symmetric matrix. MAGMA is a powerful open-source high-performance linear algebra library that plays a crucial role in dense linear algebra computation, supporting heterogeneous platform computation with different architectures. The SYR2K function algorithm in MAGMA is based on the SYRK and GEMM functions. Figure 1 As shown (B=A), the untransposed SYRK function iteratively calls a series of GEMM functions along the diagonal of the output matrix to complete the rank-k update of the symmetric matrix. Assuming matrix A has a size of n×k, to calculate αA×A... T Update the lower triangular part of C, MAGMA will update matrix A. T Subdivide into blocks of size k×nb. Connect the submatrices of A with A by calling n / nb GEMM functions. T The blocks are multiplied one by one. The row dimension of the submatrix of A starts from n. Each call decreases by nb in the row direction from top to bottom. Along the diagonal of the output matrix C, the extra computational overhead of the upper triangular off-diagonal blocks of C is minimized as much as possible, and only nb × nb / 2 right-angled small triangle blocks are calculated each time. Figure 1 As shown, based on SYRK, let the size of matrices A and B be n×k, and then call αA×B respectively. T and αB×AT The function HER2K is implemented in the same way as SYR2K, only with one more step of conjugate processing on the result, so it is not described here.

[0014] From Figure 1 It can be seen that when implementing the rank-2k update of the lower triangular symmetric matrix, the series calls GEMM multiplication to write the upper triangular part that is actually not needed. This writing causes unnecessary memory access and writing. These writes do not contribute to the correctness of the result, but still consume memory bandwidth for memory access, reducing the efficiency of memory operations. At the same time, this writing will cause waste of cache resources, invalid data will occupy cache resources, so that the part that really needs to be calculated cannot be stored, causing cache pollution, resulting in further performance degradation. In addition, writing the upper triangular part that is not needed can easily cause security problems in actual use and calculation. First of all, the user-unfriendly design makes it difficult for users to allocate space or release memory for the upper triangular part that is not needed in advance if they do not understand the design scheme, and then it is easy to cause out-of-bound access, program crash, and memory leak problems, increasing the difficulty of actual use. Even if all these problems are solved, in order to ensure the correctness of the result, the output matrix still needs to be modified, that is, the unnecessary upper triangular content needs to be zeroed, causing additional post-processing overhead.

[0015] As can be seen from Figure 1 It can be seen that the submatrix block calculated by the sequence calling matrix multiplication is limited by the triangular shape, and the calculation dimension continuously shrinks along the diagonal line, resulting in a decrease in calculation density and the inability to fully utilize GPU performance.

[0016] From the above analysis, it can be seen that the existing scheme has the following problems:

[0017] 1. Although the two series calls GEMM have good performance, there are security risks, use risks, and hidden post-processing overhead.

[0018] 2.In the process of multiple calls of GEMM, due to the fact that the general matrix multiplication is rectangular, the upper triangular or lower triangular matrix required by the output matrix cannot be spliced by cutting, thus causing unavoidable non-diagonal block additional calculation overhead. At the same time, the block size is hard-coded inside the function, and cannot be adaptively optimized and selected for different precision and dimension functions. Assuming that the size of the A and B matrices is n*k, and the size of the cutting block is nb, the series of GEMM function calls will cause the actual calculation amount to be (1+nb / n) times the theoretical required calculation amount. Since the BLAS 3-level function is a compute-intensive function, there is usually a larger scale of calculation requirement, and when the size of the matrix itself is very large or nb is set too large, the additional calculation overhead cannot be ignored, which seriously affects the performance. When nb is too small, even if the shape of the matrix itself is large, the number of sub-blocks calculated by each call of the GEMM function is still small. At the same time, the calculation of the sub-matrix contracted along the diagonal will further reduce the calculation density. SUMMARY

[0019] In view of the problems in the prior art, the purpose of the present application is to provide a high-performance computing method for the interlaced batch block recursion of SYR2K and HER2K, which uses the flexibility of handwritten code and the wide applicability of Tensile automatically generated code to increase the calculation density and improve the calculation speed. At the same time, a parameter automatic tuning framework for calculation task division and data mapping based on the structure characteristics of GPU is designed, which removes the dimension hard coding restriction, so that the parameters can be instantiated and customized for different situations.

[0020] The key points of the present application are as follows:

[0021] 1. An interlaced batch block recursion calculation parallel method based on fine-grained tuning algorithm and automatic code generation

[0022] The present application provides a diagonal interlaced batch calculation parallel method based on fine-grained tuning algorithm. Figure 1 It can be seen that the performance constraint of SYR2K and HER2K comes from the fact that the series of GEMM calls cannot meet the triangular calculation task when calling the GEMM, causing additional overhead. Therefore, according to Figure 2 , the present application proposes a fine-grained tuning interlaced batch triangular block calculation parallel algorithm, which updates the diagonal block part of the output symmetric matrix C by batch interlaced calculation of alpha AxB T and alpha BxA T ; the block recursion method is used to call the handwritten batch matrix multiplication operator or the kernel function generated by Tensile to update the output symmetric matrix C in the non-diagonal block part, so as to complete the rank-2K update. Through this design, the additional calculation overhead caused by calling the out-of-bound GEMM is effectively eliminated, and the calculation efficiency is significantly improved.Figure 2 As shown, assuming the dimension of output matrix C is n, and the dimension of split triangular sub-block is nb, the diagonal part of output matrix can be decomposed into nb diagonal sub-blocks. If n%nb≠0, there will be some remaining sub-blocks, which will be processed after the interleaved computation of αA×B T and αB×A T is completed. According to Figure 3 , after the diagonal part is processed, the remaining part can be split into different size block matrix multiplication to complete the computation. Therefore, in this part, the present application adopts block recursion to call batch matrix multiplication, and the block size increases by steps in the recursion process until the matrix block computation is completed, and the recursion terminates. Starting from the reference size nb×nb of the sub-block of output matrix C, two kernel functions are started to calculate the corresponding sub-blocks of αA×B T and αB×A T respectively each time. The sub-block is gradually expanded in a recursive manner, and the size of the sub-block is expanded to four times (i.e. 2nb×2nb) each time. This recursive expansion allows the present application to increase the computation density by increasing the batch size when calling batch matrix multiplication with smaller dimensions, and the impact of the reduced batch size on the overall computation density is smaller when processing larger matrix sub-blocks, thereby effectively utilizing the parallel computation capability and hardware resources of matrix multiplication. When the dimension of the matrix is not an integer multiple of nb, some parts cannot be completely blocked, and these residual parts will be marked in the recursion process and processed separately through additional matrix multiplication calls in the last stage of this round of loop, ensuring that all matrix elements are calculated. Since batch matrix multiplication of different sizes has different data characteristics, the present application combines the flexibility of handwritten kernel functions and the general applicability of automatic code generation to select kernel functions for different batch dimension matrix multiplication. In order to reduce the use overhead brought by Tensile, the search for the generated code library is pruned and optimized.

[0023] The core idea of the present application is to decompose the complex triangular update problem into multiple smaller sub-problems through blocking and interleaved computation, avoiding the additional computation overhead and safety problems caused by the boundary GEMM access. At the same time, it can be seen that Figure 2 the batch interleaved processing of diagonal block computation of αA×B T and αB×A T makes the input matrices A and B continuous when reading during the computation process, which is different from directly complete computation of αA×B T + αB×A TCompared with reducing the global memory jump access, the cache efficiency is improved. In addition, unlike the general matrix multiplication (GEMM), the diagonal sub-blocks of the output matrix C are usually discontinuously stored in memory. Therefore, in order to ensure the correctness of the calculation, the application additionally maintains the offset information of the input matrix A, B and the output matrix C in the kernel function. The maintenance of the offset enables the algorithm to flexibly access the non-continuous data block in the memory, and ensures the correct calculation of each sub-block. At the same time, in order to process the batch update of all diagonal sub-blocks, the application adds a parameter batch for recording the total number of diagonal sub-blocks (i.e. ) For triangular batch multiplication, the application designs a multi-stage pipeline based on double buffering and loop unrolling to realize fine-grained optimization of the kernel function. Through the combination of interleaved calculation and batch processing, the algorithm fully utilizes the parallel computing capability of the GPU, avoids the frequent memory access and out-of-bound storage in the traditional method, improves the cache utilization, and thus realizes higher performance.

[0024] 2. An automatic optimization framework based on the structural characteristics of the domestic GPU.

[0025] The application provides an automatic optimization framework for task division and data mapping parameters based on the structural characteristics of a domestic GPU. Due to different shapes and a large number of calculation blocks, this block division mode not only limits the performance of block recursion calculation, but also has requirements for batch matrix multiplication functions of different sizes. If the calculation block is set too small, the calculation problem will degenerate into an irregular matrix multiplication with a large K dimension and small M and N dimensions, which cannot make good use of the performance of the DCU, has low GEMM performance of the diagonal block and the non-diagonal block, and causes resource waste. However, if the calculation block is set too large, although the batch matrix multiplication function called by the non-diagonal block is obviously improved, the calculation of the diagonal block will restrict the overall calculation performance. In addition, there is an additional overhead when calling the batch matrix multiplication function automatically generated by the platform, and this overhead cannot be avoided.

[0026] Therefore, for the SYR2K and HER2K functions of different precisions, the application realizes the optimization of the block size and the kernel function parameters based on the structural characteristics of the domestic GPU, so as to balance the diagonal calculation block and the matrix multiplication calculation block. The application tests different block parameters and different kernel function parameters for functions of different precisions, and records the best performance under different conditions. The parameters include a global reference block, an LDS dimension, a thread block dimension, and selection of handwritten code or automatically generated code of the batch matrix multiplication. A CUDAC++ template is used to provide a unified code library for the SYR2K and HER2K kernel functions. Each set of optimization parameters is abstracted as an integer array, and the template can enable the corresponding specific kernel function to be called after the matrix input of different precisions and different dimensions, so as to improve the overall function performance.

[0027] The technical scheme of the present application is:

[0028] A dense matrix multiplication rank 2k symmetric update method for a processor, comprising the steps of:

[0029] 1) For a computing task C out = αA×B T + αB×A T + βC in , the host initializes the dimensions of matrix A and matrix B as N*K and K*N, and initializes the output matrix C out and the input matrix C in , both of which have a dimension of N*N; matrix A, matrix B and input matrix C in are copied from the host memory to the device memory;

[0030] 2) The processor of the device splits the computing task into αA×B T + αB×A T subtasks and β*C in subtasks; first, calculate β*C in and store the calculation result in the device memory where the input matrix C in is located; then, divide αA×B T + αB×A T subtasks into batch triangular matrix multiplication of diagonal blocks and recursive batch rectangular matrix multiplication of non-diagonal blocks, and take the dimension nb of the global reference block as the starting subblock dimension size for batch calculation of diagonal blocks and non-diagonal blocks;

[0031] 3) Parallel computing of staggered batch diagonal subblocks: first, get the batch number of diagonal block subblocks for batch calculation through , calculate the access storage location of the i-th matrix subblock A i of matrix A and the i-th matrix subblock B i of matrix B, allocate a grid, a thread block and a thread for calculation according to the access storage location, the batch and the subblock size, and update the batch corresponding matrix blocks C out in the output matrix C i ; then, start the diagonal subblock batch calculation kernel functions of and

[0032] respectively, and batch calculate batch triangular subblocks that can be divided by nb, and update the matrix blocks C i ;

[0033] 4) Block recursive batch matrix multiplication to calculate non-diagonal blocks: (1) assuming that the corresponding calculation starting position of the output matrix C out is n start, each round of recursive calculation, by Get this batch of non-diagonal block sub-block number batch1, staggered start assembly matrix multiplication to update And Non-diagonal sub-blocks of matrix A j The jth matrix sub-block of matrix A, B i The jth matrix sub-block of matrix B; (2) start the next round of recursive calculation, offset the starting position n start +=nb, and increase the dimension of the starting sub-block nb=nb×2; (3) repeat steps (1) ~ (2) until the starting position is out of bounds, and calculate the output matrix C out .

[0034] Further, in step 3), if there is a triangular sub-block that cannot be divided by nb, it is combined and a kernel function is started to calculate To increase the calculation density, the calculation results are accumulated to the corresponding position of the input matrix C in .

[0035] Further, in step 4), if there is a triangular sub-block that cannot be divided by nb, start 2 batch matrix multiplication to complete And The calculation results are accumulated to the corresponding position of the input matrix C in .

[0036] Further, according to the performance of the current batch calculation sub-block, a fine-grained handwritten optimization kernel function is selected or Tensile is used to generate assembly code.

[0037] Further, Tensile is used to generate the code of the assembly matrix multiplication; and the overhead of calling the assembly matrix multiplication is reduced by optimizing the optimal code search.

[0038] Further, the data types in matrix A, matrix B and input matrix C in Root is single precision, double precision, single precision complex or double precision complex.

[0039] Further, the device end memory is H2D.

[0040] The advantages of the present application are as follows:

[0041] (1) The dense matrix rank 2k update of the domestic GPU platform shuguang 7000, the present application proposes an interlaced batch block recursive high-performance calculation method based on fine-grained optimization algorithm and automatic code generation. By batch interlaced calculation αA×B T And αB×A TThe diagonal block part is not only avoided from the extra overhead in the traditional GEMM method, but also makes the interface be more reasonably used, and avoids the safety problem.

[0042] (2) The present application is verified based on the dawn 7000 processor, compared with the MAGMA open source library and the optimized version of the present application, the single precision and double precision of the SYR2K function and the single precision complex and double precision complex of the HER2K function are tested, the average speed is tested by 10 times after removing the preheating, it is found that the average speed ratio is 1.09, the double precision SYR2K function can reach the average speed ratio 1.15, and the comparison effect is as shown in Figures 5 to 8 BRIEF DESCRIPTION OF DRAWINGS

[0043] Figure 1 is the rank-2K update algorithm of the open source library MAGMA symmetric matrix, and the following triangular output matrix is taken as an example.

[0044] Figure 2 is the staggered batch diagonal block parallel method proposed by the present application, and the lower triangular output matrix which cannot be divided by the global benchmark block is taken as an example.

[0045] Figure 3 is the staggered batch block recursive calculation SYR2K schematic diagram proposed by the present application, and the above triangular output matrix is taken as an example.

[0046] Figure 4 is the diagonal staggered batch block recursive algorithm flow chart on the domestic GPU.

[0047] Figure 5 is the effect comparison of SSYR2K MAGMABLAS and the optimized version (our-blas) of the present application on the domestic GPU,

[0048] Figure 6 ​is the effect comparison of DSYR2K MAGMA BLAS on a domestic GPU and our-blas after optimization.

[0049] Figure 7 is the effect comparison of CHER2K MAGMA BLAS on a domestic GPU and our-blas after optimization.

[0050] Figure 8 is the effect comparison of ZHER2K MAGMA BLAS on a domestic GPU and our-blas after optimization. DETAILED DESCRIPTION

[0051] The application will be further described in detail below with reference to the accompanying drawings, and the examples are only used to explain the application and not to limit the scope of the application.

[0052] As shown in Figure 4 , the application designs and implements an interleaved batch computation parallel method based on fine-grained optimization algorithm and automatic code generation, and a computing task partitioning and data mapping parameter automatic optimization framework based on the architecture of a domestic GPU according to the architecture characteristics of SYR2K and HER2K and the experimental platform. Figure 4 The tasks responsible by the host and the device are clearly shown in the application. The host is mainly responsible for task starting, initialization, memory allocation and the like, and the processor is located at the device and is responsible for all core computations. The specific steps are as follows:

[0053] 1. Initialization stage

[0054] (1) Taking the conventional SYR2K as an example, the computing task is C out = αA×B T + αB×A T + βC in , the input matrix A has a dimension of N*K, the input matrix B has a dimension of K*N, and the output matrix C out and the input matrix C in have a dimension of N*N. After setting the parameters of the automatic optimization framework, the input matrices A, B and C in are copied from the host memory to the device memory (H2D), and the data type can be single-precision (FP32), double-precision (FP64), single-precision complex number and double-precision complex number according to the user input;

[0055] (2) After the above data copying is completed, the computing task is split into αA×B T + αB×A T subtasks and β*C in subtasks according to the above formula. The kernel function is first started to calculate β*C in , and the result is stored in Cin The device memory, other computing structure will be subsequently stored in the memory, so the output does not need to open additional device memory;

[0056] (3) Then, for αA×B T + αB×A T Subtask, will be in the processor side respectively αA×B T and αB×A T Calculation, and will be the results of parallel reduction to the output matrix C out . Where, since the output matrix is a symmetric matrix, only complete output matrix on the triangle or the triangle update can be. With the output matrix C out As the task division benchmark, in order to improve the parallelism and computing efficiency, block calculation is implemented. The diagonal part of the output matrix is calculated by using fine-grained optimization batch triangular matrix multiplication, and the off-diagonal block part is calculated by using handwritten and generated hybrid recursive batch rectangular matrix multiplication. The block size of the above task is in a recursive manner, and the global benchmark block dimension is set to nb, that is, the starting subblock dimension of the input matrix A is nb*K, the starting subblock dimension of the input matrix B is nb*K, and the starting subblock dimension of the output matrix C out is nb*nb. The global benchmark block dimension is used as the starting subblock dimension size of the diagonal block and the off-diagonal block batch calculation. The diagonal block is specifically represented by batch triangular matrix update, and the off-diagonal block is batch matrix multiplication, and the batch dimension is determined by the calculated C out subblock.

[0057] 2. Parallel computing staggered batch along the diagonal subblock

[0058] (1) By formula Get batch diagonal block subblock number batch, calculate the i-th subblock of input matrix A, B A i , B i Access storage location, device side will be according to the location, batch and subblock size to allocate grid, thread block and thread for calculation, update the batch corresponding matrix block C i of the output matrix.

[0059] (2) By selecting the optimal fine-grained optimization kernel function algorithm to complete the triangular matrix multiplication of batch input subblocks to complete the diagonal part calculation of the output matrix. That is, the diagonal subblock batch calculation kernel of and is started staggeredly, that is, batch triangular subblocks that can be divided by the global benchmark block dimension nb are calculated respectively, and after the calculation is completed, the batch diagonal block part C i of the output matrix can be updated, that is, the green part of the input matrix A, B and the output matrix C in Figure 3 .

[0060] (3) If there are remaining calculation tasks that cannot be divided by nb, since the amount of data is small, take the merging processing of the remaining sub-blocks that cannot be divided, and start a kernel function to complete Reduce the overhead of multiple kernel function launches, increase the calculation density, and accumulate the calculation results to the C in corresponding position.

[0061] 3. Block recursion batch matrix multiplication is used to calculate the non-diagonal blocks

[0062] (1) After the diagonal block part of the output matrix C out is updated, the calculation of the non-diagonal blocks is performed. The present application uses block recursion batch matrix multiplication to update the remaining part, and the batch calculation starting position is offset from the output matrix C out , that is, the calculation sub-block position for the diagonal block update is offset by nb along the row direction, that is Figure 3 , the calculation starting position of the output matrix is n start .The calculation process is as follows: first, calculate the number of this batch of non-diagonal block sub-blocks batch1 that can be obtained by dividing the formula , and use the handwritten code or the assembly matrix multiplication code generated by Tensile to update the non-diagonal blocks of and in batches. The code selection is based on which way has better performance for the current calculation task. When using the Tensile automatic code generation tool, Tensile needs to search for the optimal generation scheme according to the input matrix dimensions, shapes, etc. Since there are calculation tasks of the same dimensions for symmetric matrices, the same dimension sub-task search is reduced to reduce the call overhead. If there are remaining calculation tasks that cannot be divided, repeat the above steps to launch 2 batch matrix multiplications to complete and The calculation results are accumulated to the C in corresponding position mentioned in the initialization stage (2).

[0063] (2) After the above calculation is completed, the next round of recursive calculation is started, the starting position n start of the batch calculation sub-block is offset by nb, the sub-block dimension is increased by nb = nb x 2, and the above steps are repeated to complete the current round of recursive calculation.

[0064] (3) Repeat steps (1) and (2) until the calculation starting position is out of bounds, that is, n start >N, C out matrix calculation is completed.

[0065] While specific embodiments of the application have been disclosed in order to illustrate the application and to assist those skilled in the art in practicing the application, it is to be understood that various substitutions, modifications and changes can be made by those skilled in the art without departing from the spirit of the application and the scope of the appended claims. Accordingly, it is intended that the application not be limited, except by the scope of the claims.

Claims

1. A processor-oriented dense matrix-multiplication rank-2k symmetric update method, comprising the steps of: 1) For the calculation task C out = αAxB T + αBxA T + βC in , the host initializes the dimension of matrix A as N*K, the dimension of matrix B as K*N, and initializes the output matrix C out and the input matrix C in , both of which have the dimension of N*N; matrix A, matrix B and input matrix C in are copied from the host memory to the device memory; 2) The processor at the device end splits the computation task into αA x B T + αB x A T subtasks and βC x D in subtasks; first compute βC in and store the computation result in the memory at the device end where the input matrix C in is located. Then αA×B T + αB×A T The subtask is divided into batch triangular matrix multiplication of diagonal blocks and recursive batch rectangular matrix multiplication of non-diagonal blocks, and the dimension nb of the global reference block is taken as the starting subblock dimension size of batch calculation of diagonal blocks and non-diagonal blocks. 3) Parallel computing staggered batch of diagonal sub-blocks along the line: first by get the number of batch of diagonal sub-blocks batch, compute the i-th matrix sub-block A i , the i-th matrix sub-block B i of the matrix B, the access storage location, according to the access storage location, batch and sub-block size, the grid, thread block and thread will be allocated to calculate, update the corresponding matrix block C out in the output matrix C i ; Then staggered start and diagonal sub-block batch computation kernel function, respectively batch computation batch of triangular sub-blocks that can be divided by nb, update the matrix block C i ; 4) Block recursive call batch matrix multiplication to calculate off-diagonal blocks: (1) Set the output matrix C out The corresponding calculation starting position is n start , and each round of recursive calculation is performed by obtaining the number of batch1 subblocks for batch calculation of off-diagonal blocks, and starting the assembly matrix multiplication in batches to update the off-diagonal subblocks of and , wherein A j is the jth matrix subblock of matrix A, and B i is the jth matrix subblock of matrix B; (2) Start the next round of recursive calculation, offset the starting position n start of the batch calculation subblock by nb, and increase the dimension nb of the starting subblock by nb = nb x 2; (3) Repeat steps (1) and (2) until the calculation starting position is out of bounds, and the output matrix C out is calculated.

2. The method of claim 1, wherein, In step 3), if there is a triangular sub-block that cannot be divided by nb, it is merged and a kernel function is calculated To increase the calculation density, the calculation result is accumulated to the input matrix C in The corresponding position.

3. The method of claim 2, wherein, In Step 4), if there are triangular sub-blocks that cannot be divided by nb, 2 batch matrix multiplications are initiated to complete and The calculation results are all accumulated to the input matrix C in corresponding positions.

4. The method according to claim 1 or 2 or 3, characterized in that, selecting to use a fine-grained hand-written optimized kernel function or to use Tensile to generate assembly code according to the performance of a sub-block calculated in a current batch.

5. The method according to claim 1 or 2 or 3, characterized in that, generating code of the assembly matrix multiplication using Tensile; reducing the overhead of calling the assembly matrix multiplication by optimizing optimal code search.

6. The method according to claim 1 or 2 or 3, characterized in that, Matrix A, Matrix B and input matrix C in The data type root in the matrix A, the matrix B and the input matrix C is single precision, double precision, single complex or double complex.

7. The method of claim 1, wherein, the device-side memory is H2D.

Images