High-performance implementation method and device for dense matrix multiplication for SW39000 processor

Abstract

The present invention discloses a high-performance implementation method of dense matrix multiplication for SW39000 processor, comprising: performing task block division on non-transposed matrix A, non-transposed matrix B and problem matrix C at the main memory level, generating a plurality of sub-matrices δA il , submatrix δB lj With the submatrix δC ij ; When the submatrix δC ij When updating, at the LDM level of the slave core array, the submatrix δC ij , submatrix δA il and submatrix δB lj Divide the task into blocks and generate sub-blocks εC uv , sub-block εA uw and sub-block εB wv Map to the corresponding slave core; Submatrix δA il The tth column block and the submatrix δB lj Multiply the t-th row block by the matrix δC t ; Based on the submatrix δC ij With each matrix δC t , update the submatrix δC in main memory ij , to obtain the output result of the problem matrix C. The present invention fully exploits the parallel computing capability of the slave core array through reasonable data layout on the main memory, the LDM of the slave core and the register, uses triple buffering technology to hide the DMA memory access overhead and RMA communication overhead, and uses instruction reordering technology to improve the floating-point computing efficiency of GEMM.

Description

Technical Field

[0001] The present invention relates to the field of basic linear algebra library BLAS (Basic Linear Algebra Subprograms) implementation, and in particular to a high-performance implementation method and device for dense matrix multiplication (GEMM) for SW39000 processor. Background Art

[0002] BLAS (Basic Linear Algebra Subprograms) is a basic linear algebra subprogram library, which mainly includes basic operations of vectors and matrices. It is one of the most basic and important mathematical libraries and is widely used in scientific computing, weather forecasting, astrophysics and other fields. The BLAS library is the core of many professional software. Among them, the high-performance implementation of dense matrix multiplication (GEMM) in BLAS level 3 functions is particularly important. It is the most important part of the BLAS library and is widely used. Both the traditional numerical computing field and the emerging artificial intelligence field rely heavily on GEMM. In addition, many large-scale scientific computing applications in high-performance fields such as biopharmaceuticals, ocean circulation, nuclear energy development, etc., rely on highly optimized GEMM at the bottom to obtain high operating efficiency. In particular, the performance of Linpack, a benchmark program used for the world's TOP500 supercomputer ranking, relies heavily on the performance of GEMM. Practices in numerical matrix analysis, deep learning, etc. have shown that GEMM is of vital significance to improving the computing speed of applications and giving full play to the performance of high-performance computers.

[0003] Dense matrix multiplication (GEMM) mainly implements: C = α × op(A) × op(B) + β × C, supporting four data types: real single / double precision and complex single / double precision. A, B, C are matrices of size m×k, k×n, and m×n, respectively. α and β are scalar parameters. op(X) can be X, X T (transpose of X), X H (conjugate transpose of X). GEMM is computationally intensive, and its performance depends heavily on whether the hardware characteristics of the machine can be effectively utilized. That is, how to reasonably arrange data in the storage structure, improve performance by computing hidden memory access, and make full use of the parallel computing resources of high-performance computers to increase data reuse rate. This is a great challenge for the high-performance implementation of GEMM.

[0004] There is a lot of parallel optimization work for GEMM on high-performance computing platforms. For example, WangQian et al. studied the optimization technology of dense linear algebra subroutines for mainstream multi-core processors, which can directly generate highly optimized assembly code for GEMM. Zhang et al. optimized GEMM on Loongson 3A CPU. Ali M et al. highly optimized GEMM on multi-core DSP using OpenMP. The above are all works on GEMM on multi-core processors. With the rapid development of GPU accelerators, the optimization work of GEMM on GPU has attracted more and more attention. For example, LiangJuanjuan et al. designed and implemented the core function GEMM of the level 3 BLAS library on GPU using CUDA. Nath R et al. gave a parallel algorithm for implementing high-performance GEMM on Nvidia GPU.

[0005] SW39000 is a multi-core processor with a heterogeneous architecture. On the new generation of Shenwei supercomputers based on the SW39000 multi-core processor architecture, there is currently no customized high-performance GEMM function deployed, and the existing open source mathematical library has low performance on this platform and cannot provide effective performance support for applications. Therefore, it is urgent to design and implement a high-performance GEMM function for this multi-core platform to fully utilize the hardware characteristics of the Shenwei multi-core processor and meet the urgent needs of upper-level applications for high-performance GEMM on the Shenwei multi-core platform. Summary of the invention

[0006] The present invention provides a high-performance implementation method and device of dense matrix multiplication (GEMM) for SW39000 processor, so as to meet the demand for dense matrix multiplication (GEMM) on SW39000 multi-core processor and solve the problem of low performance of existing open source mathematical library. Dense matrix multiplication (GEMM) mainly implements: C = α × op (A) × op (B) + β × C, supports four data types of real single / double precision and complex single / double precision, where A, B, C are matrices of size m × k, k × n, m × n respectively, α, β are scalar parameters, op (X) can be X, X T (transpose of X), X H (Conjugate transpose of X). The storage system of the SW39000 processor includes a three-level storage structure of main memory, slave core LDM, and register from bottom to top. Limited by the space size of each level of storage structure, the present invention designs a three-level block parallel matrix multiplication algorithm to transfer data to the slave core and use the computing resources of the slave core for calculation. The slave core array of the SW39000 processor includes 64 slave cores organized into a grid structure of 8 rows and 8 columns.

[0007] The technical solution of the present invention includes:

[0008] A high-performance implementation method for dense matrix multiplication for SW39000 processor, the steps comprising:

[0009] 1) At the main memory level, the non-transposed matrix A, the non-transposed matrix B and the problem matrix C are divided into task blocks to generate several sub-matrices δA il , submatrix δB lj With the submatrix δC ij , where 0≤i≤M, 0≤j≤N, 0≤l≤K, The sizes of the non-transposed matrix A, non-transposed matrix B and problem matrix C are m×k, k×n and m×n respectively, and the submatrix δA il , submatrix δB lj With the submatrix δC ij The size is divided into b m ×b k 、b k ×b n With b m ×b n ;

[0010] 2) When the submatrix δC ij When updating, at the LDM level of the slave core array, the submatrix δC ij , submatrix δA il and submatrix δB lj Divide the task into blocks and generate sub-blocks εC uv , sub-block εA uw and sub-block εB wv Mapped to the corresponding slave core, where 0≤u≤T, 0≤w≤T, 0≤v≤T, and the parameter T is set according to the number of slave cores;

[0011] 3) Sub-matrix δA il The tth column block and submatrix δB lj Multiply the t-th row block by the matrix δC t , and based on the submatrix δC ij With each matrix δC t , and get the submatrix δC ij The update result δC′ ij ;

[0012] 4) Based on the update result δC′ ij , update the submatrix δC in main memory ij , to get the output result of the problem matrix C.

[0013] Furthermore, the data types of the non-transposed matrix A, the non-transposed matrix B and the problem matrix C include: real single precision, real double precision, complex single precision or complex double precision.

[0014] Further, obtain the submatrix δC to be updated ij The method includes: using the NMK loop order, traversing each submatrix δC in turn ij .

[0015] Furthermore, the submatrix δC ij , submatrix δA il and submatrix δB lj The method of transferring data to the LDM of the slave core array includes: using a DMA data transfer method.

[0016] Further, the matrix δC is generated by the following steps: t :

[0017] 1) According to the submatrix δA il Column block neutron block εA uw With the submatrix δB lj Subblock εB in row block wv , get the sub-block εC uv The core where it is located;

[0018] 2) Sub-block εC uv The slave core where it is located obtains the sub-block εA uw With sub-block εB wv ;

[0019] 3) At the register level, sub-block εC uv , sub-block εA uw With sub-block εB wv Divide the tasks separately to obtain corresponding computing sub-blocks;

[0020] 4) Based on sub-block εA uw With sub-block εB wv The computational sub-block uses floating-point vector multiplication and addition instructions and loop unrolling, and uses the sub-block εC uv The calculation sub-block is updated as a unit to obtain the sub-block εC uv The update result of

[0021] 5) Using each sub-block εC uv The update result generates the matrix δC t .

[0022] Furthermore, when the sub-block εA uw With sub-block εB wv Not in subblock εC uv When the slave core is located, the slave core obtains the sub-block εA based on the RMA communication method. uw With sub-block εB wv .

[0023] Furthermore, the sub-block εC uvThe size of the computational sub-block is q×r m × n , sub-block εA uw The size of the computational sub-block is q×r m ×p k , sub-block εB wv The size of the computational sub-block is p k × n , where q represents the number of elements that the vector operation instruction can run simultaneously, and sub-block εB wv The size of the computational sub-block is p k × n , r m Represents the storage sub-block εA uw The number of registers of the computational sub-block, r n Represents the storage sub-block εB wv The number of registers of the computation sub-block.

[0024] Further, update the result δC′ ij =δC ij +∑δC t .

[0025] A storage medium stores a computer program, wherein the computer program is configured to execute the above method when running.

[0026] An electronic device includes a memory and a processor, wherein the memory stores a program for executing the above method.

[0027] The present invention has achieved the following technical effects:

[0028] Based on the architecture features of the SW39000 multi-core processor, the present invention designs a multi-level block algorithm based on a data sharing strategy, uses a reasonable data layout on the main memory, the LDM of the slave core and the register, and fully exploits the parallel computing capabilities of the slave core array. During the calculation process, a data sharing strategy based on the RMA communication mechanism is designed to achieve data sharing between slave cores, provide data for the calculation process, and enhance data reusability. The triple buffering technology is used to hide the DMA memory access overhead and the RMA communication overhead, and the instruction reordering technology is used to improve the floating-point calculation efficiency of GEMM. BRIEF DESCRIPTION OF THE DRAWINGS

[0029] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative labor.

[0030] Figure 1 The figure is a schematic diagram of the overall flow of the high-performance implementation method of dense matrix multiplication (GEMM) for the SW39000 processor of the present invention.

[0031] Figure 2 Schematic diagram of the matrix block in the main memory.

[0032] Figure 3 Schematic diagram of submatrix partitioning from kernel LDM.

[0033] Figure 4 Schematic diagram of the mapping between matrix elements and slave kernel arrays.

[0034] Figure 5 8A il The first column block and δB lj Schematic diagram of multiplication of the first row block.

[0035] Figure 6 The broadcast sub-block εA from the core RMA row in column 1 u1 Schematic diagram.

[0036] Figure 7 The first row of the slave core RMA column broadcast sub-block εB 1v Schematic diagram.

[0037] Figure 8 It is the performance acceleration ratio of the GEMM of the present invention and the open source library GotoBLAS. DETAILED DESCRIPTION

[0038] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0039] The high-performance implementation method of dense matrix multiplication (GEMM) of the present invention comprises the following steps:

[0040] Step 1: Figure 2 As shown, at the main memory level, matrix C, matrix A and matrix B are divided into two matrices of size b m ×b n 、b m ×b k 、b k ×b n The sub-matrix δC ij , δA il and δB lj , where 0≤i≤M, 0≤j≤N, 0≤l≤K, δCij represents the submatrix of the i-th row and j-th column of matrix C, δA il represents the submatrix of the i-th row and l-th column of matrix A, δB lj represents the submatrix of the lth row and jth column of matrix B. Using the cyclic order of NMK, the submatrix δC ij Unit, serially update matrix C, where δC ij + = α × (δA il )×(δB lj ),in,

[0041] Step 2: Update submatrix δC ij When the DMA data transfer mechanism provided by SW39000 is used, the required δC ij , δA il and δB lj Transfer from main memory to the LDM of the slave core array. Figure 3 As shown, at the LDM level, the submatrix δC ij , δA il and δB lj Further divided into sub-blocks εC uv , εA uw and εB wv , and εC uv Mapped to the slave core at row u and column v, εA uw Mapped to the slave core at row u and column w, εB wv Mapped to the slave core of row w and column v. Among them, 0≤u≤7, 0≤w≤7, 0≤v≤7, εC uv Represents the submatrix δC ij The sub-block of row u and column v, εA uw Denotes the submatrix δA il The sub-block of the uth row and wth column, εB wv Denotes the submatrix δB ij The sub-block at the wth row and vth column of .

[0042] Step 3: Use the submatrix δA that has been transferred to LDM il and δB lj Update δC ij , 64 slave cores parallelly calculate δA il ×δB lj Each slave core is responsible for updating the sub-block εC mapped to the local LDM uv , where εC uv + = (εA uw )×(εB wv), each slave core obtains sub-blocks εA and εB on other slave cores through the RMA communication mechanism provided by SW39000. The slave core updates sub-block εC uv When using DMA, triple buffering technology is used to overlap DMA data transmission with RMA communication and calculation.

[0043] Step 4: Each slave core updates the sub-block εC uv When the calculation is mainly performed (εA uw )×(εB wv ), and the assembly instructions provided by SW39000 are used to perform the calculation. At the register level, the sub-block εC uv , εA uw and εB wv Further divided into computational sub-blocks, sub-block εC uv The size of the computational sub-block is 8×r m × n , sub-block εA uw The size of the computational sub-block is 8×r m ×p k , sub-block εB wv The size of the computational sub-block is p k × n . Use floating-point vector multiply-add instructions and loop unrolling to update the sub-block εC in units of computation sub-blocks uv Among them, r m Represents the storage sub-block εA uw The number of registers of the computational sub-block, r n Represents the storage sub-block εB wv The number of registers of the computational sub-block, r m × n Represents the storage sub-block εC uv The number of registers of the computation sub-block.

[0044] Step 5: Use the DMA data transmission mechanism provided by SW39000 to update the sub-matrix δC ij Write back to main memory from the LDM of the slave core array.

[0045] Furthermore, step one includes:

[0046] When the NMK cycle order is used to serially update the matrix C, δC ij is reused, i.e. uploaded to the δC in the slave core LDM ij Retained in LDM until no longer needed. m 、b n 、b k The selection is limited by the size of the slave core LDM space, that is, the corresponding storage space applied for cannot exceed the size of the SW39000 slave core LDM space.

[0047] Furthermore, step 2 includes:

[0048] Submatrix δC ij , δA il and δB lj Further divided into sub-blocks εC uv , εA uw and εB wv , the mapping method between the sub-block and the slave core array is as follows: Figure 4 As shown, the submatrix δC ij , δA il and δB lj The rows are evenly divided into 8 row blocks, with sizes p m ×b n 、p m ×b k 、p k ×b n , and map the 8 row blocks to the 8 rows of slave cores respectively, that is, the i-th row block is mapped to the i-th row of slave cores, and when each row block is mapped to the slave core of the corresponding row, each row block is evenly cut into 8 parts along the column direction and mapped to the 8 slave cores of the row, and each slave core obtains the corresponding sub-block εC uv , εA uw and εB wv , their sizes are p m ×p n 、p m ×p k 、p k ×p n , where p m =b m / 8, p n =b n / 8, p k =b k / 8.

[0049] Furthermore, step three includes:

[0050] δC ij The update requires 8 iterations, each of which is il The column block and δB lj The row block multiplication is implemented, for example, the t-th iteration is implemented by δA il The tth column block (0≤t≤7) and δB lj It is realized by multiplying the t-th row block of . Figure 5 As shown, it is the first iteration, that is, when t = 1, δA il The first column block and δB lj The first row block is multiplied to get δC ij A part of δA ilEach column block of contains 8 sub-blocks εA uw , δB lj Each row block of contains 8 sub-blocks εB wv , so, at each iteration, δA il The column block and δB lj The row block multiplication corresponds to the 64 sub-block multiplication εA uw ×εB wv , which is calculated in parallel by 64 slave cores, each of which is responsible for updating εC in the local LDM uv , execute εA uw ×εB wv After each iteration, the slave core accumulates the calculation results to the local LDM's εC uv In each iteration, each slave core calculates εA uw ×εB wv Calculate the required εA uw and εB wv It may come from the local LDM or need to be obtained from other slave cores through the RMA communication mechanism. For example, at the tth iteration, the slave cores in other columns need to obtain the sub-block εA of the slave core in the tth column. ut , the slave cores in other rows need to obtain the sub-block εB of the slave core in row t tv Therefore, the slave core in the tth column broadcasts the sub-block εA in the local LDM through RMA row broadcast. ut To other slave cores in the same row, the slave core in row t broadcasts the sub-block εB in the local LDM by RMA column broadcast tv To other slave cores in the same column. Figure 6 As shown, in the first iteration, the first column of the slave core RMA row broadcast sub-block εA u1 To other slave cores in the same industry, such as Figure 7 As shown, in the first iteration, the first row of the slave core RMA column broadcasts the sub-block εB 1v To other slave cores in the same column. The triple buffering technology is specifically to use the DMA data transmission mechanism to access the matrix A elements and matrix B elements required for the current calculation from the main memory to the slave core LDM, and use the RMA communication mechanism to broadcast the matrix A elements and matrix B elements required for the last calculation to complete the last matrix-matrix multiplication operation, that is, DMA data transmission, RMA communication, and calculation are performed in parallel, and the calculation time is slightly longer than the DMA data transmission time and RMA communication time.

[0051] Furthermore, step four includes:

[0052] εC uv The size of the computational sub-block is 8×r m × n , using r m registers to store the elements of matrix A, using rn registers to store the elements of matrix B, r m × n registers to store the accumulated matrix C elements. Limited by the number of registers of the slave core being 64, we have: m × n +r m +r n ≤64. Use floating-point vector multiply-add instructions to update the sub-block εC in units of calculation sub-blocks. uv , a three-layer loop is adopted, the innermost loop is expanded 8 times, the number of multiplication and addition operations in a single loop is increased, the dependency between floating-point vector multiplication and addition instructions and memory access instructions is avoided as much as possible, and the floating-point vector multiplication and addition instructions and memory access instructions without dependency are issued simultaneously.

[0053] The following example involves dense matrix multiplication (GEMM) involving matrices A and B in non-transposed form. It mainly implements: C = α × A × B + β × C, supporting four data types: real single / double precision and complex single / double precision, where A, B, and C are matrices of size m × k, k × n, and m × n, respectively, and α and β are scalar parameters. The specific implementation steps include:

[0054] Step 1: Task division:

[0055] like Figure 2 As shown, the present invention divides the matrix C, matrix A, and matrix B into task blocks in the main memory, and the matrix C generates a total of ((M+1)×(N+1)) sub-matrices δC ij , the matrix A generates a total of ((M+1)×(K+1)) sub-matrices δA il , the matrix B generates a total of ((K+1)×(N+1)) sub-matrices δB lj ,in, 0≤i≤M, 0≤j≤N, 0≤l≤K. Each submatrix δC ij Considered as a subproblem, it contains ((M+1)×(N+1)) subproblems in total, and each subproblem completes the following operation: δC ij + = α × (δA il )×(δB lj ).

[0056] The present invention adopts the NMK loop order to traverse each sub-problem in turn. Assuming that the current problem is sub-problem j, the present invention uses the DMA data transmission mechanism provided by SW39000 to transfer the δC required by the sub-problem ij , δA il and δB lj Transfer from main memory to the LDM of the slave core array. Figure 3 As shown, at the LDM level, the submatrix δC ij , δAil and δB lj Further divided into sub-blocks εC uv , εA uw and εB wv , where 0≤u≤7, 0≤w≤7, 0≤v≤7. uv It is considered as a subtask, which contains 64 subtasks in total and is completed in parallel by 64 slave cores. Each slave core completes the following operations: εC uv + = (εA uw )×(εB wv ).

[0057] Assuming that the slave core CPE (0≤CPE≤63) is currently processing the subtask, the present invention converts the sub-block εC uv , εA uw and εB wv Further divided into computational sub-blocks, sub-block εC uv The size of the computational sub-block is 8×r m × n , sub-block εA uw The size of the computational sub-block is 8×r m ×p k , sub-block εB wv The size of the computational sub-block is p k × n Among them, r m Represents the storage sub-block εA uw The number of registers of the computational sub-block, r n Represents the storage sub-block εB wv The number of registers of the computational sub-block, r m × n Represents the storage sub-block εC uv The number of registers of the computation sub-block.

[0058] Step 2: Resource allocation:

[0059] Assuming that the current problem is subproblem j, it completes the following operations: δC ij + = α × (δA il )×(δB lj ). Update δC ij When the DMA data transfer mechanism provided by SW39000 is used, the required δC ij , δA il and δB lj Transfer from main memory to the LDM of the slave core array. Figure 4 As shown, at the LDM level, the submatrix δC ij , δA il and δB ljFurther divided into sub-blocks εC uv , εA uw and εB wv , and εC uv Mapped to the slave core at row u and column v, εA uw Mapped to the slave core at row u and column w, εB wv Mapped to the slave core at row w and column v.

[0060] Step 3: Execute calculation:

[0061] Assuming the subproblem j is currently being processed, it completes the following operations: δC ij + = α × (δA il )×(δB lj ).like Figure 5 As shown, δC ij The update requires 8 iterations, each of which is il The column block and δB lj The row block multiplication is implemented, for example, in the first iteration, δA il The first column block and δB lj The first row block is multiplied to get δC ij A part of δA il Each column block of contains 8 sub-blocks εA uw , δB lj Each row block of contains 8 sub-blocks εB wv , so, at each iteration, δA il The column block and δB lj The row block multiplication corresponds to the multiplication of 64 sub-blocks, which are calculated in parallel by 64 slave cores. Each slave core is responsible for updating εC in the local LDM. uv , using floating-point vector multiply-add instructions and loop unrolling to complete the calculation: εA uw ×εB wv After each iteration, the slave core accumulates the calculation results to the local LDM's εC uv middle.

[0062] For example, δC ij The initial value of is 0, and its result needs 8 iterations to be obtained. In the first iteration, δA il The first column block and δB lj The result of multiplying the first row block is δC1. After the first iteration is completed: δC ij =δC ij +δC1, then we get δC ij The first partial sum; the second iteration, δA il The second column block and δB lj The result of multiplying the second row block is δC2. After the second iteration is completed: δC ij=δC ij +δC1+δC2, we get δC ij The second part of the sum; and so on, after each iteration, the result of multiplying the column block and the row block is accumulated to δC ij When all 8 iterations are completed, we get δC ij The result, δC ij =δC ij +∑δC t .

[0063] At each iteration, each slave core calculates εA uw ×εB wv When the calculation is completed based on the calculation sub-block, the required εA is calculated uw and εB wv It may come from the local LDM or need to be obtained from other slave cores through the RMA communication mechanism. Each slave core obtains the εA required for calculation from other slave cores. uw and εB wv After that, εA is calculated in units of calculation sub-blocks. uw ×εB wv εA uw ×εB wv A three-layer loop is adopted, and the innermost loop is expanded 8 times, which increases the number of multiplication and addition operations in a single loop. The dependency between floating-point vector multiplication and addition instructions and memory access instructions is avoided as much as possible, and the floating-point vector multiplication and addition instructions and memory access instructions without dependency are issued simultaneously.

[0064] This embodiment uses the GotoBLAS mathematical library to verify the performance acceleration effect of the present invention. The problem scale selected in the embodiment is 1024, 2048, 4096, and 8192 matrix scales, and the precision selected is real single precision, real double precision, complex single precision, and complex double precision. Compared with the performance acceleration ratio of GEMM of the open source library GotoBLAS, the average acceleration ratio of the present invention relative to GotoBLAS is 696.56, and the highest acceleration ratio is 978.59.

[0065] This embodiment is only for the SW39000 multi-core processor, and the contents of the present invention are transplanted to other platforms after simple transformation, or there is no creative improvement on the task segmentation and data sharing strategy of the present invention, or simple optimization is performed on the basis of the present invention during the execution calculation stage. In essence, it still does not deviate from the contents covered by the present invention and still falls within the scope of protection of the present invention.

[0066] Parts of the present invention that are not described in detail belong to the common knowledge of those skilled in the art.

[0067] The embodiments described above are only descriptions of specific examples of the present invention and are not intended to limit the scope of the present invention. Without departing from the design spirit of the present invention, various modifications and improvements made to the technical solutions of the present invention by ordinary technicians in this field should all fall within the protection scope determined by the claims of the present invention.

Claims

1. A high-performance implementation method for dense matrix multiplication for SW39000 processor, the steps comprising: 1) At the main memory level, the non-transposed matrix A, the non-transposed matrix B and the problem matrix C are divided into task blocks to generate several sub-matrices δA il , submatrix δB lj With the submatrix δC ij , where 0≤i≤M, 0≤j≤N, 0≤l≤K, The sizes of the non-transposed matrix A, non-transposed matrix B and problem matrix C are m×k, k×n and m×n respectively, and the submatrix δA il , submatrix δB lj With the submatrix δC ij The size is divided into b m ×b k , b k ×b n With b m ×b n ; 2) When the submatrix δC ij When updating, at the LDM level of the slave core array, the submatrix δC ij , submatrix δA il and submatrix δB lj Divide the task into blocks and generate sub-blocks εC uv , sub-block εA uw and sub-block εB wv Mapped to the corresponding slave core, where 0≤u≤T, 0≤w≤T, 0≤v≤T, parameter T is set according to the number of slave cores, sub-block εC uv The size of the computational sub-block is q×r m × n , sub-block εA uw The size of the computational sub-block is q×r m ×p k , sub-block εB wv The size of the computational sub-block is p k × n , q represents the number of elements that the vector operation instruction can run simultaneously, r m Represents the storage sub-block εA uw The number of registers of the computational sub-block, r n Represents the storage sub-block εB wv The number of registers of the computation sub-block; 3) Sub-matrix δA il The tth column block and submatrix δB lj Multiply the t-th row block by the matrix δC t , and based on the submatrix δC ij With each matrix δC t , and get the submatrix δC ij The update result δC′ ij ; Wherein, the sub-matrix δA il The tth column block and submatrix δB lj Multiply the t-th row block by the matrix δC t ,include: According to the submatrix δA il Column block neutron block εA uw With the submatrix δB lj Subblock εB in row block wv , get the sub-block εC uv The core where it is located; From the sub-block εC uv The slave core where it is located obtains the sub-block εA uw With sub-block εB wv ; Among them, when the sub-block εA uw With sub-block εB wv Not in subblock εC uv When the slave core is located, the slave core obtains the sub-block εA based on the RMA communication method. uw With sub-block εB wv ; At the register level, the sub-block εC uv , sub-block εA uw With sub-block εB wv Divide the tasks separately to obtain corresponding computing sub-blocks; Based on sub-block εA uw With sub-block εB wv The computational sub-block uses floating-point vector multiplication and addition instructions and loop unrolling, and uses the sub-block εC uv The calculation sub-block is updated as a unit to obtain the sub-block εC uv The update result of Using each sub-block εC uv The update result generates the matrix δC t ; 4) Based on the update result δC′ ij , update the submatrix δC in main memory ij , to get the output result of the problem matrix C.

2. The method according to claim 1, characterized in that The data types of the non-transposed matrix A, the non-transposed matrix B and the problem matrix C include: real single precision, real double precision, complex single precision or complex double precision.

3. The method according to claim 1, characterized in that Get the submatrix δC to be updated ij The method includes: using the NMK loop order, traversing each submatrix δC in turn ij .

4. The method according to claim 1, characterized in that Sub-matrix δC ij , submatrix δA il and submatrix δB lj The method of transferring data to the LDM of the slave core array includes: using a DMA data transfer method.

5. The method according to claim 1, characterized in that Update result δC′ ij =δC ij +∑δC t .

6. A storage medium storing a computer program, wherein: The computer program is configured to execute the method according to any one of claims 1 to 5 when executed.

7. An electronic device, comprising a memory and a processor, wherein the memory stores a computer program, and the processor is configured to run the computer program to perform the method according to any one of claims 1 to 5.

Images