Memory architecture oriented double precision general matrix multiplication optimization method and system

Abstract

The application belongs to the technical field of high-performance computing, and aims at solving the problems of limited computing power and access efficiency, and provides a double-precision general matrix multiplication optimization method and system for a memory architecture, wherein a main core receives an operation request, determines a unified blocking rule by using a multi-constraint optimization model, and decomposes a matrix into a plurality of sub-matrix blocks according to a core array topology; a slave core receives the sub-matrix blocks issued by the main core, divides the sub-matrix blocks into sub-matrix small blocks based on the unified blocking rule, loads the sub-matrix small blocks into local data memory independent buffers based on a DMA double-buffer protocol, divides the sub-matrix small blocks in the buffer into SIMD vectors according to the characteristics of a SIMD unit, and alternately switches vectorization calculation and cache operation through different independent buffers according to an iteration period; after all the slave cores complete the calculation, the main core collects the results written back to the main memory by the slave cores to obtain the final operation result, and the application realizes a double breakthrough of computing power and memory access efficiency.

Description

TECHNICAL FIELD

[0001] The present application belongs to the technical field related to high-performance computing, and particularly relates to a double-precision general matrix multiplication optimization method and system for a memory architecture. BACKGROUND

[0002] The statements in this section merely provide background information related to the present application and do not necessarily constitute prior art.

[0003] As the most core basic operator in the field of scientific computing, double-precision general matrix multiplication (DGEMM) usually accounts for 60% to 80% of the entire computing time in applications such as quantum chemistry, molecular dynamics, and density functional theory. With the development of large-scale scientific applications such as first-principle calculations to the level of millions of cores, the performance of double-precision general matrix multiplication has become a key bottleneck restricting the overall application efficiency. In particular, in the real-space density functional theory software, the eigenvalue solving process based on the Chebyshev filter subspace iteration algorithm includes projection Hamiltonian calculation and subspace rotation operations, which need to call double-precision general matrix multiplication thousands or even tens of thousands of times in a single self-consistent field iteration. In the computing scenario of thousands of atoms, the cumulative execution time of DGEMM accounts for as high as 54.6%, which seriously affects the overall performance of the software and the feasibility of large-scale applications.

[0004] The mathematical expression of double-precision general matrix multiplication operation is: wherein, , , is a double-precision floating-point matrix, and is a scalar coefficient. The floating-point operation amount of this operation is: . Among them is the multiplication and addition operation of matrix multiplication, is the scaling and addition operation of . When the matrix size is large, the dominant term is . For matrix operations with a size of , a single double-precision general matrix multiplication call needs to perform about floating-point operations. In the single SCF iteration of the CheFSI algorithm, double-precision general matrix multiplication is called 24525 times, and the total calculation amount is as high as floating-point operations, making double-precision general matrix multiplication optimization a key to improving overall performance.

[0005] The existing high-performance BLAS library optimization work is mainly aimed at the traditional hardware cache architecture design. Intel MKL, OpenBLAS and other libraries have achieved close to the hardware theoretical peak performance on mainstream architectures such as x86 and ARM through hierarchical blocking strategy and automatic tuning technology. These methods rely on the multi-level cache automatic replacement mechanism, hardware prefetch logic and cache coherence protocol to optimize data access, and the core optimization idea is based on the premise of automatic management of hardware cache.

[0006] The Godson SW26010pro many-core processor adopts a completely different design concept, and its computing processing core is equipped with a 256 KB local data memory to replace the traditional hardware cache. This software-managed temporary storage architecture requires programmers to explicitly control each data transfer between main memory and LDM, including the selection of transfer timing, the determination of data transfer granularity, the specification of target address, and the management of buffer life cycle. LDM does not have the automatic replacement mechanism, prefetch logic and cache coherence protocol of the hardware cache. Any data access that does not have the target data in the LDM will cause program errors or a sharp decline in performance. This architectural feature makes the traditional cache-oriented architecture optimized BLAS library unable to run directly on the Godson platform or severely degraded in performance, and must be specially redesigned and optimized. More seriously, the main memory access bandwidth of the Godson core is relatively limited, and the 64 slave cores of each core group need to share a memory channel with a theoretical bandwidth of only 51.2 GB / s. In dense matrix operations, if the memory access delay cannot be effectively hidden, a serious memory bottleneck will easily occur.

[0007] The real space DFT method solves the eigenvalue problem by using the Chebyshev filter subspace iteration algorithm, and the matrix operation presents a special mode different from the traditional scientific computing application. In each iteration of the CheFSI algorithm, the projection Hamiltonian calculation and the subspace rotation operation need to be performed, involving NoTrans-NoTrans, Trans-NoTrans, NoTrans-Trans, Trans-Trans and other matrix transpose combinations. Each transpose combination has completely different requirements for SIMD vectorization: NoTrans-NoTrans can use memory continuity for direct vector loading and FMA operation, while Trans-NoTrans needs to handle stride access and horizontal reduction problems. The existing BLAS library is usually deeply optimized for NoTrans-NoTrans, and for other transpose combinations, it often uses the detour strategy of calculating after transposition, which will introduce additional storage overhead and data movement overhead in the environment with limited LDM capacity. The CheFSI algorithm needs to call DGEMM thousands to tens of thousands of times in a single SCF iteration, and the matrix size of each call is usually thousands by hundreds. This frequent small and medium-sized call mode puts higher requirements on the initialization overhead and buffer reuse efficiency of the optimization method.

[0008] The prior art has the following systematic problems: first, the matrix block parameter selection faces the problem of multiple constraint coupling, which needs to meet the double buffering space, SIMD alignment, data reuse efficiency and DMA transmission granularity constraints, and the existing method lacks a systematic parameter determination method; second, the double buffering implementation has deficiencies in synchronous timing selection, buffer switching overhead and boundary processing, and the pipeline overlap effect is limited; third, the SIMD vectorization lacks a complete framework for uniformly processing the four transpose combinations of NN, TN, NT and TT, and multiple sets of code paths need to be maintained; fourth, the master-slave core heterogeneous model has technical problems in task allocation among 64 slave cores, data race avoidance and load balancing processing.

[0009] The above problems result in that the performance of double-precision general matrix multiplication is much lower than the theoretical peak value, which becomes a key bottleneck of real space DFT calculation on domestic supercomputers. There is an urgent need for a systematic optimization method for double-precision general matrix multiplication facing software managed memory architecture, which realizes collaborative optimization under multiple constraints. SUMMARY

[0010] To overcome the above problems of the prior art, the application provides a double-precision general matrix multiplication optimization method and system facing a memory architecture, which improves the computing power and memory efficiency and can quickly support the needs of fields such as scientific computing and engineering simulation which have high requirements for matrix operation efficiency.

[0011] To achieve the above purpose, the application adopts the following technical solutions:

[0012] In a first aspect, the application provides a double-precision general matrix multiplication optimization method for a memory architecture, comprising:

[0013] The master core receives a double-precision general matrix multiplication operation request, determines a unified blocking rule by using a multi-constraint optimization model, and decomposes the double-precision general matrix into a plurality of sub-matrix blocks according to the slave core array topology.

[0014] The slave core receives the sub-matrix blocks issued by the master core, divides the sub-matrix blocks into sub-matrix small blocks based on the unified blocking rule, loads the sub-matrix small blocks into independent buffer areas of the local data storage based on the DMA double-buffering protocol, and after the loading is completed, divides the sub-matrix small blocks into SIMD vectors for calculation according to the characteristics of the SIMD unit of the slave core; by alternately switching between different independent buffer areas, the pipeline parallelism of the SIMD vectorization calculation and the DMA data transmission is realized, and the slave core completes the calculation of all the allocated sub-matrix small blocks.

[0015] After all the slave cores complete the calculation, the master core synchronously waits and collects the results written back to the main memory by the slave cores, and obtains the final double-precision general matrix operation result.

[0016] In a second aspect, the application provides a double-precision general matrix multiplication optimization system for a memory architecture, comprising:

[0017] The master core is configured to receive a double-precision general matrix multiplication operation request, determine a unified blocking rule by using a multi-constraint optimization model, and decompose the double-precision general matrix into a plurality of sub-matrix blocks according to the slave core array topology.

[0018] The slave core is configured to receive the sub-matrix blocks issued by the master core, divide the sub-matrix blocks into sub-matrix small blocks based on the unified blocking rule, load the sub-matrix small blocks into independent buffer areas of the local data storage based on the DMA double-buffering protocol, and after the loading is completed, divide the sub-matrix small blocks into SIMD vectors for calculation according to the characteristics of the SIMD unit of the slave core; by alternately switching between different independent buffer areas, the pipeline parallelism of the SIMD vectorization calculation and the DMA data transmission is realized, and the slave core completes the calculation of all the allocated sub-matrix small blocks.

[0019] The master core is configured to, after all the slave cores complete the calculation, synchronously wait and collect the results written back to the main memory by the slave cores, and output the final double-precision general matrix operation result.

[0020] In a third aspect, the application provides an electronic device, comprising a memory, a processor, and computer instructions stored in the memory and running on the processor, wherein when the computer instructions are run by the processor, the method of the first aspect is completed.

[0021] In a fourth aspect, the present application provides a computer readable storage medium for storing computer instructions, which, when executed by a processor, complete the method of the first aspect.

[0022] The above one or more technical solutions have the following beneficial effects:

[0023] In the present application, through the combination optimization of master-slave core cooperation, three-level blocking, DMA buffering and SIMD vectorization, the dual breakthroughs of computing power and memory efficiency are realized, which can quickly support the field requirements of scientific computing, engineering simulation and other high requirements for matrix operation efficiency.

[0024] In the present application, a three-level blocking strategy is designed, the master core determines a unified blocking rule by using a multi-constraint optimization model, and the master core decomposes a double-precision general matrix into a plurality of sub-matrix blocks according to the slave core array topology; the slave core divides the sub-matrix blocks into sub-matrix small blocks based on the unified blocking rule; and the slave core divides the sub-matrix small blocks in the local data memory independent buffer into SIMD vectors according to the characteristics of the SIMD unit. Each layer of the three-level blocking precisely matches the hardware characteristics of the master-slave core, LDM and SIMD unit, avoiding resource waste caused by single blocking.

[0025] In the present application, a unified vectorization framework is designed for four matrix transpose combinations, and the strategy is automatically selected according to the data access mode; a three-level strategy is adopted for boundary processing, which avoids slowing down the overall efficiency while ensuring that no data is missed.

[0026] The advantages of the additional aspects of the present application will be partially given in the following description, partially become obvious from the following description, or be understood through the practice of the present application. BRIEF DESCRIPTION OF DRAWINGS

[0027] The drawings accompanying the specification of the present application serve to provide a further understanding of the present application, and the illustrative embodiments of the present application and their descriptions serve to explain the present application, and do not constitute an improper limitation on the present application.

[0028] Figure 1 It is a double-precision general matrix multiplication optimization system architecture diagram in the embodiment of the present application;

[0029] Figure 2 It is a matrix blocking decomposition schematic diagram in the embodiment of the present application;

[0030] Figure 3 It is a DMA double-buffering pipeline timing diagram in the embodiment of the present application;

[0031] Figure 4 It is a double-precision general matrix multiplication optimization method flowchart in the embodiment of the present application. DETAILED DESCRIPTION

[0032] It should be noted that the following detailed description is exemplary in nature and is intended to provide further description of the application. Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this application belongs.

[0033] It should be noted that the terminology used herein is for the purpose of describing particular embodiments only and is not intended to be limiting of example embodiments according to the present application.

[0034] In the case of no conflict, the embodiments in the present application and the features in the embodiments can be combined with each other.

[0035] Term explanation:

[0036] DGEMM refers to double-precision general matrix multiplication, which is one of the most basic and important linear algebra operations in scientific computing, and its operation form is C = aAB + bC, where A, B, and C are double-precision floating point matrices, and a and b are scalar coefficients.

[0037] LDM refers to local data memory, which is a software-managed on-chip high-speed memory, and is essentially different from the traditional hardware-managed cache. Programmers need to explicitly control the data transfer time, transmission granularity and space allocation between main memory and LDM through software.

[0038] DMA refers to direct memory access, which is a technology that allows hardware subsystems to access system memory independently of the central processor, and is used in the present application to realize asynchronous data transmission between main memory and LDM.

[0039] SIMD refers to single instruction multiple data, which is a parallel computing architecture that allows a single instruction to process multiple data elements simultaneously. The SIMD unit of the Godson core is 256 bits wide and can process 8 double-precision floating point numbers simultaneously.

[0040] MPE refers to management processing element, which adopts a traditional cache architecture and is responsible for task scheduling, parameter transmission and system I / O operation.

[0041] CPE refers to computing processing element, which is equipped with software-managed LDM and is specifically designed to execute compute-intensive tasks.

[0042] DFT refers to Density Functional Theory, which is an important theoretical method for studying the electronic structure of many-electron systems in quantum mechanics, and is widely used in the fields of material science and quantum chemistry.

[0043] CheFSI refers to Chebyshev Filtered Subspace Iteration, which is an efficient numerical method for solving large-scale symmetric eigenvalue problems, and is used to solve the eigenvalues and eigenvectors of the Kohn-Sham equation in real-space DFT calculations.

[0044] FMA refers to Fused Multiply-Add, which can complete multiplication and addition operations in one clock cycle, improving the throughput of floating-point operations. Arithmetic intensity refers to the number of floating-point operations that can be completed per byte of data transfer, which is used to measure the computational intensity of an algorithm and is an important indicator for determining whether the performance bottleneck of an algorithm is memory-bound or computation-bound.

[0045] Scratchpad Memory refers to on-chip memory that needs to be explicitly managed by software, in contrast to cache memory that is automatically managed by hardware.

[0046] FLOPs / Byte represents the number of floating-point operations per byte;

[0047] NN → NoTrans-NoTrans indicates that both matrices are not transposed;

[0048] NT → NoTrans-Trans indicates that the first matrix is not transposed and the second matrix is transposed;

[0049] TN → Trans-NoTrans indicates that the first matrix is transposed and the second matrix is not transposed;

[0050] TT → Trans-Trans indicates that both matrices are transposed;

[0051] DMA_iget() is a DMA asynchronous read function, used to asynchronously load data from main memory to LDM;

[0052] DMA_wait() is a DMA synchronous waiting function, used to wait for the completion of asynchronous transmission;

[0053] DMA_iput() is a DMA asynchronous write function, used to asynchronously write data from LDM back to main memory;

[0054] simd_load() is a SIMD vector load function, which loads vector data from memory to SIMD registers;

[0055] simd_vmad() is a SIMD vector fused multiply-add function, which performs vector-level a x b + c operation;

[0056] simd_reduce_add() is a SIMD vector reduction add function, which sums up each element of the vector to a scalar;

[0057] doublev8 is a vector data type of 8 double-precision floating-point numbers;

[0058] current_buf represents the current buffer index variable used;

[0059] next_buf represents the next buffer index variable to be used;

[0060] attribute_((aligned(64))) represents a compiler alignment attribute declaration, which specifies that the data is aligned according to the 64-byte boundary;

[0061] Buffer[0] calculates Block k while Buffer[1] loads Block k+1, which means that the 0th group of buffers performs Block k calculation while the 1st group of buffers loads Block k+1 data;

[0062] scatter-gather represents a scatter-gather instruction;

[0063] stride represents stride access;

[0064] gather represents a gather instruction;

[0065] padding represents a padding instruction.

[0066] Embodiment one

[0067] This embodiment aims at the performance bottleneck of DGEMM, i.e., double-precision general matrix multiplication, on the software-managed memory architecture of ShenWei many-core processors, and proposes a systematic collaborative optimization method that integrates LDM-aware blocking strategy, DMA double-buffering protocol and SIMD vectorization mechanism. The core technologies include:

[0068] (1) Under the condition of 256KB LDM capacity, a blocking parameter optimization method is established, which meets the constraints of double-buffering space, SIMD alignment, data reuse efficiency and DMA transmission granularity;

[0069] (2) Design and calculate the transmission pipeline overlap guarantee mechanism, ensure the effectiveness of double buffering through time analysis and adaptive synchronization strategy;

[0070] (3) Design a unified SIMD vectorization framework for four matrix transpose combinations, and select continuous vectorization or horizontal reduction strategy according to the transpose attribute;

[0071] (4) Realize efficient parallel task scheduling and load balancing of master-slave core heterogeneous architecture, avoid data competition and process irregular dimension.

[0072] The embodiment discloses a double-precision general matrix multiplication optimization method for a memory architecture, comprising:

[0073] The master core receives a double-precision general matrix multiplication operation request, determines a unified blocking rule by using a multi-constraint optimization model, and decomposes the double-precision general matrix into a plurality of sub-matrix blocks according to the slave core array topology;

[0074] The slave core receives the sub-matrix block issued by the master core, divides the sub-matrix block into a sub-matrix small block based on the unified blocking rule, loads the sub-matrix small block into an independent buffer area of the local data storage based on the DMA double-buffering protocol; after loading is completed, the slave core divides the sub-matrix small block into a SIMD vector for calculation according to the characteristics of the SIMD unit of the slave core; by alternately switching between different independent buffer areas, the pipeline parallelism of SIMD vectorization calculation and DMA data transmission is realized, until the slave core completes the calculation of all the allocated sub-matrix small blocks;

[0075] After all the slave cores complete the calculation, the master core synchronously waits and collects the results written back to the main memory by the slave cores, and obtains the final double-precision general matrix operation result.

[0076] In the embodiment, as shown in Figure 4 The master core MPE executes an initialization phase, the initialization phase includes parameter legality check, calculation blocking parameter, parameter structure construction and slave core start; the main loop phase is divided into master core MPE flow and slave core MPE flow; the master core MPE flow is waiting for the slave core, and the slave core MPE flow includes reading the slave core coordinates, circulating task allocation, LDM buffer area allocation, double-buffering pipeline in K-dimensional loop, SIMD calculation and asynchronous result writing. The master core MPE executes an end phase, the end phase includes waiting for the slave core to return, reading the complete result matrix C from the main memory, beta scaling processing and returning the result matrix C.

[0077] At the overall architecture level, a heterogeneous parallel model of master-slave core cooperation is adopted. The master core is responsible for global decomposition, parameter packaging and result collection, and 64 slave cores are started by the athread_spawn() function for SPMD (Single Program Multiple Data) parallel execution. The slave cores allocate independent buffers in their respective local data memories LDM through the_thread_local keyword to avoid data competition. The slave cores directly write back to the main memory through DMA after completing the calculation, and the master core waits through the athread_join() function.

[0078] The system architecture of the embodiment adopts a three-layer progressive design, as shown in Figure 1 The application interface layer provides a standard DGEMM interface and performs parameter legality checking. The task scheduling and allocation layer runs on the master core MPE and includes matrix block division, packaged parameter structure, slave core start and result collection. The parallel computing execution layer runs on the 64 slave cores CPE and includes task allocation, LDM buffer management, DMA data transmission, double-buffer scheduling, SIMD calculation and boundary processing.

[0079] In terms of the LDM-aware block strategy, the embodiment proposes a block parameter determination method based on a multi-constraint optimization model, which formalizes the parameter selection as a constraint satisfaction problem. Five key constraints are systematically considered: (1) LDM capacity constraint. The DGEMM double-buffer mechanism needs to store two sets of A matrix buffers, two sets of B matrix buffers and a single set of C matrix buffers in the LDM, and the five buffers are A matrix double buffer (A_buffer[0], A_buffer[1]), B matrix double buffer (B_buffer[0], B_buffer[1]), C matrix single buffer (C_buffer). A matrix and B matrix need double buffering to realize pipeline parallelism of calculation and transmission, and C matrix only needs single buffering as the accumulation result. Space is also reserved for DMA descriptors, loop counters and other local data. The capacity constraint is expressed as 2MbKb+2KbNb+MbNb+Overhead<=LDM_SIZE, where Mb, Nb and Kb are the three dimensions of the matrix block, Overhead is the auxiliary data space, and LDM_SIZE represents the total capacity of the local data memory.

[0080] (2) SIMD alignment constraint. To fully utilize the 512-bit SIMD vector unit of the Godson slave core, the dimensions of the matrix block must be an integer multiple of the vector width VLEN, and the double v8 type VLEN is 8. VLEN represents the vector length.

[0081] (3) Data multiplexing efficiency constraint, the algorithm is changed from memory limited to computation limited by improving arithmetic intensity;

[0082] (4) DMA transfer efficiency constraint, the efficiency of the Godson from core DMA engine significantly decreases when the transfer granularity is less than 512 bytes, and the block byte number 8MbKb>=512 bytes is required to avoid transfer fragmentation;

[0083] (5) Load balancing constraint, the block size is preferentially selected to be an integer multiple of the common matrix dimension to reduce the boundary block imbalance.

[0084] To simplify the solution and achieve inter-core load balancing, the embodiment adopts a cubic block strategy to set Mb=Nb=Kb=B, and the capacity constraint is simplified to 5B²+Overhead<=LDM_SIZE. For the Godson from core 256KB capacity, 32768 double-precision floating-point numbers can be stored, and the Overhead is about 2000 elements. It is derived that B²<=(32768-2000) / 5=6153.6, that is, B<=78.4. Combined with the SIMD alignment constraint that B must be a multiple of 8, the candidate values are 64 and 72. Considering the DMA transfer efficiency, when B=64, the single block byte number is 64x64x8=32768 bytes, which is much larger than the 512 byte threshold. Comprehensive all constraints, the block size Mb=Nb=Kb=64 is finally determined. Under this configuration, the A matrix double-buffering occupies 8192 elements, the B matrix double-buffering occupies 8192 elements, and the C matrix single-buffering occupies 4096 elements, totaling 20480 double-precision numbers, occupying a storage capacity of 160KB, with an LDM utilization rate of 62.5%, and the remaining 96KB used for local variables, DMA descriptors and stack data. This method is universal and can be extended to software management memory architectures with different LDM capacities and SIMD widths.

[0085] In terms of three-level block hierarchical optimization, the embodiment establishes a three-layer optimization system of inter-core parallelism, LDM block caching and SIMD vectorization. The first level of block division divides the tasks among the cores, divides the matrix A (MxK) by rows, the matrix B (KxN) by columns, and the result matrix C (MxN) into corresponding submatrix blocks, and uses the 8x8 core array topology to perform two-dimensional cyclic distribution on the submatrix blocks of the result matrix C. Each core is responsible for calculating the submatrix block of the result matrix C allocated to it and its corresponding A row block and B column block. The second level of block division is performed by the core, and the data block is further processed based on the block parameters determined by the multi-constraint optimization model. The submatrix block allocated to each core is further processed in the respective LDM by 64x64 submatrix small blocks, which are completely loaded into the LDM and efficiently multiplexed by the double-buffering mechanism. The third level of block division utilizes the 8-element width of the SIMD vector unit at the microkernel level for vectorization processing, ensuring that the innermost loop operates in multiples of 8 to maximize the utilization efficiency of the vector unit.

[0086] As Figure 2 shown, a three-level tiling strategy is demonstrated. The first level is inter-core parallel tiling, which divides the result matrix C into a 8x8 grid corresponding to a 64-core array. The second level is LDM capacity-adaptive tiling, which further divides each sub-matrix of the core into 64x64 blocks to organize double-buffered layout in a 256KB LDM. The third level is SIMD vectorized tiling, which processes each 64x64 block in an 8-element vector width. Figure 2 The LDM space allocation (A_buffer[2], B_buffer[2], C_buffer) and data flow path, main memory -> LDM -> SIMD, are demonstrated in the figure. The three-level hierarchical organization optimizes data locality at different granularities. The main memory (Main Memory) is the system shared main memory, which is accessed by the master core and the slave core. A_buffer[2] represents the double-buffered array of A matrix; B_buffer[2] represents the double-buffered array of B matrix; C_buffer represents the single-buffered C matrix.

[0087] The block parameter selection method has universality and scalability. For architectures with different LDM capacities, the optimal block size can be recalculated by adjusting the constraints. For example, for an architecture with a 128KB LDM capacity, the optimal block size is about 48x48x48; for an architecture with a 512KB or larger capacity, a block configuration of 80x80x80 or larger can be used. For scenarios where the three dimensions of the matrix M, N, and K are significantly unbalanced, a rectangular tiling strategy can be used, setting Mb, Nb, and Kb as independent parameters, and solving the optimal combination through multivariate constraint optimization.

[0088] In terms of master-slave core collaboration and parallel strategy, the embodiment uses two-dimensional cyclic distribution to achieve load balancing. The block of the i-th row and j-th column of the matrix C is assigned to the slave core with row coordinate i mod 8 and column coordinate j mod 8, so that the number of sub-matrix blocks allocated to each slave core is at most 1. The slave core declares a thread-local LDM buffer array (A_buffer[2][64x64], B_buffer[2][64x64], C_buffer[64x64]) through __thread_local, which is physically completely isolated to avoid data race. The slave core independently processes according to the block index list: DMA loads the matrix block into the LDM, executes the double-buffered pipeline calculation, and asynchronously writes back the result. The master core waits for all slave cores to complete through athread_join(). Among them, A_buffer[2][64x64] represents the declaration of the double-buffered array of A matrix, each group of 64x64 elements; B_buffer[2][64x64] represents the declaration of the double-buffered array of B matrix, each group of 64x64 elements; C_buffer[64x64] represents the declaration of the C matrix buffer, 64x64 elements.

[0089] In the DMA double buffering protocol, the embodiment designs a ping-pong buffer pipeline to realize the pipeline parallelism of "calculating the current block and asynchronously loading the next block". A_buffer[2][64x64] and B_buffer[2][64x64] and C_buffer[64x64] are declared in the LDM of each slave core, that is, A matrix configures two independent buffers, B matrix configures two independent buffers, and C matrix configures one buffer, and the two sets of buffers of matrix A and the two sets of buffers of matrix B are synchronously and alternately used through indexes 0 and 1. The execution process is divided into three stages: (1) the initialization stage pre-fetches the first data block to Buffer[0] through DMA_iget() and calls DMA_wait() to wait for completion; (2) the main loop stage, for iteration k, calculates the current buffer index current_buf and the next buffer index next_buf through the modulo operation k mod 2, to realize the alternate use of the two sets of buffers; if there is a subsequent data block, then asynchronously load block k+1 to Buffer[next_buf] through DMA_iget() (non-blocking call), while performing SIMD calculation on Buffer[current_buf], and calling DMA_wait() at the end of iteration to wait for the completion of transmission; (3) the end stage writes back the result matrix block through DMA. Buffer[0] represents the 0th set of buffers; Buffer[1] represents the 1st set of buffers.

[0090] There are five buffers in the LDM of one slave core: the 0th set of buffers of A matrix A_buffer[0]: 64x64 elements; the 1st set of buffers of A matrix A_buffer[1]: 64x64 elements; the 0th set of buffers of B matrix B_buffer[0]: 64x64 elements; the 1st set of buffers of B matrix B_buffer[1]: 64x64 elements; and the unique buffer of C matrix C_buffer: 64x64 elements, used for accumulation.

[0091] For example, the ping-pong buffer workflow is as follows:

[0092] Suppose the slave core needs to process 4 blocks of K dimensions, namely Block 0, Block 1, Block 2, and Block 3:

[0093] Initialization: DMA loads Block 0 to A_buffer[0] and B_buffer[0];

[0094] Iteration k=0:

[0095] SIMD computation: compute using A_buffer[0], B_buffer[0], result accumulated to C_buffer;

[0096] DMA transfer: concurrently load Block 1 to A_buffer[1], B_buffer[1] asynchronously;

[0097] Iteration k = 1:

[0098] SIMD computation: compute using A_buffer[1], B_buffer[1], result accumulated to C_buffer;

[0099] DMA transfer: concurrently load Block 2 to A_buffer[0], B_buffer[0], overwrite data used up;

[0100] Iteration k = 2:

[0101] SIMD computation: compute using A_buffer[0], B_buffer[0];

[0102] DMA transfer: concurrently load Block 3 to A_buffer[1], B_buffer[1];

[0103] Iteration k = 3:

[0104] SIMD computation: compute using A_buffer[1], B_buffer[1];

[0105] No following block, no need of DMA transfer;

[0106] This way, the computation of current block is achieved while the next block is loaded, as shown in the time axis overlap. Figure 3

[0107] To ensure the overlap efficiency of double-buffered computation and transfer, the embodiment proposes a guarantee mechanism based on time analysis. Define T_comp as the time for the slave to perform SIMD computation on a 64x64 sub-matrix block, which corresponds to the time for one iteration of computation in the double-buffered pipeline; T_dma is the time required for transferring a 64x64 A matrix block and a 64x64 B matrix block from the main memory to the LDM buffer of the slave through DMA; T_sync is the time for the slave to wait for the completion of DMA transfer, which is the time consumed by the slave to call the DMA_wait() function to wait for the completion of the asynchronous transfer started by the last DMA_iget(). When T_comp≥T_dma, the computation time is sufficient to mask the transfer time, and the waiting time T_sync is close to 0; otherwise, the slave needs to wait for the completion of the transfer.

[0108] ​This embodiment adopts four strategies: (1) early prefetch, initiating DMA transfer of the next block as soon as the core starts computing the current 64x64 sub-matrix tile using the SIMD unit, instead of waiting until the current block computation is finished, maximizing the overlap window of computation and transfer; (2) batch synchronization, using scatter-gather functionality to transfer multiple non-contiguous memory regions in a single request, reducing DMA initiation and synchronization overhead; (3) non-blocking writeback, asynchronously writing back results via DMA_iput(), concurrently with the next block computation; (4) adaptive buffer depth, dynamically switching between double buffering and triple buffering modes according to the T_comp / T_dma ratio. By placing DMA_wait() before iteration k computation instead of at the end of iteration k-1, precise synchronization timing control is achieved, maximizing pipeline efficiency.

[0109] As shown in Figure 3 the timeline is expanded horizontally, and divided vertically into DMA transfer lanes, CPE computation lanes, and synchronization event tracks. Figure 3 Buffer[0] is first loaded in the initialization phase, the ping-pong buffer pipeline (Buffer[0] computes Block k while Buffer[1] loads Block k+1) in the main loop phase, and timing comparisons between ideal (T_comp≥T_dma, T_sync=0) and non-ideal cases. Figure 3 The coordination between DMA_iget() asynchronous loading, SIMD computation, and DMA_wait() synchronization is shown.

[0110] In addition to the double buffering scheme, this embodiment can implement a triple buffering or adaptive multi-buffering scheme. The triple buffering scheme maintains three independent buffers in the LDM, which can more effectively hide latency in scenarios where transfer delay is significantly greater than computation time. The adaptive multi-buffering scheme dynamically adjusts the buffer depth by monitoring the ratio of T_comp and T_dma at runtime, switching between double buffering, triple buffering, and single buffering to adapt to different computation-to-transfer ratios. For computation-intensive scenarios, it can be simplified to a single-buffered synchronous DMA scheme, saving buffer space for larger tile sizes.

[0111] As an optional implementation, the specific implementation of the three-buffer scheme is that A_buffer[3], B_buffer[3] and C_buffer are maintained in the LDM, the first two blocks are preloaded to buffer[0] and buffer[1] in the initialization stage, and the current calculation buffer is determined by k%3 and the next loading buffer is determined by (k+2)%3 in the main loop, so as to realize the three-level pipeline of '1 calculation, 1 transmission and 1 readiness'. The capacity constraint is adjusted to 3MbKb+3KbNb+MbNb+Overhead≤LDM_SIZE, and the block size needs to be correspondingly reduced to adapt to the LDM capacity.

[0112] In the aspect of the SIMD vectorization technology, the embodiment designs a unified vectorization framework for four matrix transpose combinations, and automatically selects a strategy according to a data access mode. For the NN and NT combinations, data is arranged continuously, a continuous vectorization strategy is adopted, 8 elements are loaded by using simd_load(), and a multiplication and addition instruction is fused by using simd_vmad(). For the TN and TT combinations, data is accessed by using a stride, a horizontal reduction strategy is adopted, non-continuous data is loaded by using a vector gather instruction, and reduction is performed by using simd_reduce_add().

[0113] The boundary processing adopts a three-level strategy to cope with the case that the matrix dimension is not a multiple of 8: (1) the 8-element multiple part uses a complete SIMD vectorization kernel, and 8 elements are processed each time; (2) when there are remaining 4-7 elements, a SIMD masked vectorization technology is used, the activation state of a vector channel is controlled through a mask, and only the valid elements are processed; (3) when there are less than 4 elements, a scalar loop is used to process them one by one, so as to avoid the initialization overhead of the masked vectorization. For example, for a matrix with a dimension of 67, the first 64 elements are processed by using 8 complete vectors, and the remaining 3 elements are processed by using a scalar loop.

[0114] The memory optimization includes that the LDM array is aligned to 64 bytes through the declaration of __attribute__((aligned(64))), so that the simd_load() is completed in a single cycle; the main matrix is padded to a multiple of 64 to avoid boundary judgment; and the column-major matrix adjusts the loop order to ensure continuous access.

[0115] The core innovation of the embodiment lies in the deep integration of the three optimization technologies. The LDM block strategy reduces the memory access amount from O(MNK) to O(MNK / B) and improves the arithmetic intensity from O(1) to 16 FLOPs / Byte, so that the DGEMM is changed from memory-limited to calculation-limited; the DMA double-buffering hides the memory delay through pipeline overlap on the basis of the block; the SIMD vectorization adopts the differential strategy for different transpose combinations to maximize the calculation throughput; and the three cooperate to produce a nonlinear multiplication effect, so as to realize the performance order-of-magnitude improvement of the DGEMM.

[0116] The embodiment realizes significant performance improvement on the Shenwei SW26010pro many-core processor. For the density functional theory calculation of 1000 carbon atom system with a total of 24525 calls of DGEMM, the DGEMM execution time is reduced from 7354.458 seconds of the benchmark to 33.926 seconds, realizing 216.79 times acceleration, reducing the DGEMM time ratio from 54.6% to 0.6%, and the end-to-end software performance is improved by 2.57 times.

[0117] The technology contribution decomposition table shows the synergistic effect of the optimization strategy: the LDM-aware blocking strategy contributes 19.70 times the base speed, increases the arithmetic intensity from O(1) to 16 FLOPs / Byte, and changes DGEMM from memory-limited to computation-limited; On this basis, DMA double buffering and SIMD vectorization are introduced, which contributes an additional 11.00 times acceleration, 19.70x11.00≈216.79, the double buffering synchronization latency ratio is less than 5%, and the SIMD vectorization coverage is 92%. Comparative experiments show that the single technology has limited acceleration effect (13.8~26.1 times), the combination of two technologies improves the effect (88.7~189.5 times), and the combination of three technologies realizes the optimal performance (216.79 times), proving the multiplication effect of synergistic optimization rather than simple linear superposition.

[0118] The embodiment establishes a complete optimization methodology for software management memory architecture, including a blocking parameter determination method of a multi-constraint optimization model, a calculation and transmission overlap guarantee mechanism based on time analysis, and a unified multi-transposition combined vectorization framework, providing a systematic technical solution for high-performance scientific computing on domestic many-core processors.

[0119] Embodiment Two

[0120] The purpose of the embodiment is to provide a double-precision general matrix multiplication optimization system for memory architecture, which includes:

[0121] The main core is used for receiving a double-precision general matrix multiplication operation request, determining a unified blocking rule by using a multi-constraint optimization model, and decomposing a double-precision general matrix into a plurality of sub-matrix blocks according to a core array topology.

[0122] The slave core is configured to receive a sub-matrix block issued by the master core, divide the sub-matrix block into sub-matrix small blocks based on a uniform block rule, and load the sub-matrix small blocks into independent buffers of a local data memory based on a DMA double-buffering protocol; after the loading is completed, the sub-matrix small blocks are divided into SIMD vectors for calculation according to characteristics of a self-SIMD unit; the pipeline parallelism of the SIMD vectorization calculation and the DMA data transmission is realized by alternately switching between different independent buffers until all the sub-matrix small blocks allocated by the slave core are calculated.

[0123] The master core is configured to, after all the slave cores complete the calculation, synchronously wait for and collect results written back to the master memory by the slave cores, and output a final double-precision general matrix operation result.

[0124] In more embodiments, there are also provided:

[0125] An electronic device includes a memory and a processor, and computer instructions stored in the memory and executed on the processor, when executed by the processor, complete the method described in Embodiment One. For brevity, this will not be described here.

[0126] It should be understood that in the embodiments, the processor can be a central processing unit CPU, and the processor can also be other general-purpose processors, digital signal processors DSPs, application-specific integrated circuits ASICs, ready-to-program gate arrays FPGA or other programmable logic devices, discrete gates or transistor logic devices, discrete hardware components, etc. The general-purpose processor can be a microprocessor or the processor can also be any conventional processor, etc.

[0127] The memory can include read-only memory and random access memory, and provide instructions and data to the processor, and a part of the memory can also include non-volatile random access memory. For example, the memory can also store device type information.

[0128] A computer readable storage medium is configured to store computer instructions, when executed by a processor, complete the method described in Embodiment One.

[0129] The method in Embodiment One can be directly embodied as hardware processor execution completion, or executed by a combination of hardware and software modules in the processor. The software module can be located in a storage medium mature in the art, such as random access memory, flash memory, read-only memory, programmable read-only memory, or electrically erasable programmable memory, register, etc. The storage medium is located in the memory, and the processor reads the information in the memory, and combines the hardware to complete the steps of the above method. To avoid repetition, this will not be described in detail here.

[0130] Those skilled in the art can understand that the units and algorithm steps of each example described in combination with the present embodiment can be realized in electronic hardware or a combination of computer software and electronic hardware. Whether the functions are realized in hardware or software manner depends on the specific application and design constraints of the technical solution. The skilled person can use different methods to realize the described functions for each specific application, but such implementation should not be considered beyond the scope of the present application.

[0131] Although the specific embodiments of the present application are described above in combination with the drawings, it is not a limitation on the scope of protection of the present application. Those skilled in the art should understand that various modifications or variations made by those skilled in the art on the basis of the technical solutions of the present application without creative labor are still within the scope of protection of the present application.

Claims

1. A method for optimizing double-precision general matrix multiplication oriented towards a memory architecture, characterized in that, The method comprises the following steps: The main core receives a double-precision general matrix multiplication operation request, determines a unified blocking rule by using a multi-constraint optimization model, and decomposes the double-precision general matrix into a plurality of sub-matrix blocks according to the topology of the slave core array; The slave core receives the sub-matrix blocks issued by the main core, divides the sub-matrix blocks into sub-matrix small blocks based on the unified blocking rule, and loads the sub-matrix small blocks into independent buffer areas of the local data memory based on the DMA double-buffer protocol; after the loading is completed, the slave core divides the sub-matrix small blocks into SIMD vectors for calculation according to the characteristics of the SIMD unit of the slave core; by alternately switching between different independent buffer areas, the pipeline parallelism of SIMD vectorization calculation and DMA data transmission is realized, and the calculation of all the sub-matrix small blocks allocated by the slave core is completed; After all the slave cores complete the calculation, the main core synchronously waits and collects the results written back to the main memory by the slave cores, and obtains the final double-precision general matrix operation result; The pipeline parallelism of SIMD vectorization calculation and DMA data transmission is realized by alternately switching between different independent buffer areas, and specifically comprises the following steps: In the initialization stage, the first data block is preloaded into the 0th group of the A matrix buffer area and the 0th group of the B matrix buffer area, respectively; SIMD vector calculation is performed by using the 0th group of the A matrix buffer area and the 0th group of the B matrix buffer area, the results are accumulated into the C matrix buffer area, and the second data block is loaded into the 1st group of the A matrix buffer area and the 1st group of the B matrix buffer area at the same time; SIMD vector calculation is performed by using the 1st group of the A matrix buffer area and the 1st group of the B matrix buffer area, the results are accumulated into the C matrix buffer area, and the third data block is loaded into the 0th group of the A matrix buffer area and the 0th group of the B matrix buffer area at the same time if there are subsequent data blocks, so as to realize the pipeline parallelism of SIMD vectorization calculation and DMA data transmission; The method further comprises the following steps: the slave core initiates DMA transmission of the next block immediately when starting to calculate the current sub-matrix small block by using the SIMD unit; the dispersion-gathering function is adopted to combine the transmission requests of a plurality of discontinuous memory areas; and the double-buffer or triple-buffer mode is dynamically switched according to the ratio of the time for performing SIMD calculation of the sub-matrix small block to the time required for transferring the sub-matrix small block from the main memory to the LDM buffer area of the slave core.

2. The memory architecture oriented double precision general matrix multiplication optimization method of claim 1, wherein, The multi-constraint optimization model comprises an LDM capacity constraint, a SIMD alignment constraint, a data multiplexing efficiency constraint, a DMA transmission efficiency constraint and a load balancing constraint.

3. The memory architecture oriented double precision general matrix multiplication optimization method of claim 1, wherein, The unified blocking rule is determined by using the multi-constraint optimization model, and the double-precision general matrix is decomposed into a plurality of sub-matrix blocks according to the topology of the slave core array, and specifically comprises the following steps: The main core divides the matrix A, the matrix B and the result matrix C into sub-matrix blocks; The two-dimensional cyclic distribution of the sub-matrix blocks of the result matrix C is performed by using the topology of the slave core array, so that each slave core is responsible for calculating the distributed sub-matrix blocks of the result matrix C and the corresponding row blocks of the matrix A and the column blocks of the matrix B; The determined unified sub-matrix small block size is determined by the multi-constraint optimization model, and the blocking size of the sub-matrix small block is packaged and transmitted to the slave core.

4. The memory architecture oriented double precision general matrix multiplication optimization method of claim 1, wherein, The sub-matrix small block is divided into SIMD vectors for calculation according to the characteristics of the SIMD unit of the slave core, and specifically comprises the following steps: For NN and NT combination, continuous vectorization strategy is adopted, and element-by-element multiply-add instruction is used; For TN and TT combination, horizontal reduction strategy is adopted, and non-continuous data is loaded by vector instruction and reduction is performed.

5. The memory architecture oriented double precision general matrix multiplication optimization method of claim 1, wherein, In the process of executing SIMD vectorization calculation from the core, the matrix dimension 8 element multiple part uses full vectorization kernel, and when there are 4-7 elements left, the SIMD mask vectorization technology is used to control the activation state of the vector channel through the mask, and only the effective elements are processed; when there are less than 4 elements left, scalar loop is used to process them one by one.

6. A double precision general matrix multiplication optimization system for a memory- oriented architecture, characterized by, It comprises: The main core is used for receiving a double-precision general matrix multiplication operation request, determining a unified blocking rule by using a multi-constraint optimization model, and decomposing the double-precision general matrix into a plurality of sub-matrix blocks according to the array topology of the slave core; The slave core is used for receiving the sub-matrix blocks issued by the main core, dividing the sub-matrix blocks into sub-matrix small blocks based on the unified blocking rule, and loading the sub-matrix small blocks into independent buffer areas of the local data memory based on the DMA double-buffer protocol; after loading is completed, the sub-matrix small blocks are divided into SIMD vectors for calculation according to the characteristics of the SIMD unit itself; by alternately switching between different independent buffer areas, the pipeline parallelism of SIMD vectorization calculation and DMA data transmission is realized, and the calculation of all allocated sub-matrix small blocks is completed by the slave core; wherein, by alternately switching between different independent buffer areas, the pipeline parallelism of SIMD vectorization calculation and DMA data transmission is realized, specifically: In the initialization stage, the first data block is preloaded into the 0th group of the A matrix buffer area and the 0th group of the B matrix buffer area respectively; SIMD vector calculation is performed using the 0th group of the A matrix buffer area and the 0th group of the B matrix buffer area, and the result is accumulated in the C matrix buffer area, while the second data block is loaded into the 1st group of the A matrix buffer area and the 1st group of the B matrix buffer area asynchronously; SIMD vector calculation is performed using the 1st group of the A matrix buffer area and the 1st group of the B matrix buffer area, and the result is accumulated in the C matrix buffer area, and if there are subsequent data blocks, the third data block is loaded into the 0th group of the A matrix buffer area and the 0th group of the B matrix buffer area asynchronously, realizing the pipeline parallelism of SIMD vectorization calculation and DMA data transmission; The main core is used for synchronously waiting and collecting the results written back to the main memory by the slave core after all the slave cores complete the calculation, and outputting the final double-precision general matrix operation result; The slave core further comprises a DMA transmission initiation unit for initiating DMA transmission of the next block immediately when the current sub-matrix small block starts to be calculated using the SIMD unit; a scatter-gather function is used to combine the transmission requests of multiple non-continuous memory areas; and a double-buffer or triple-buffer mode is dynamically switched according to the ratio of the time for executing SIMD calculation of the sub-matrix small block to the time required for transferring the sub-matrix small block from the main memory to the LDM buffer area of the slave core.

7. An electronic device, comprising: The computer program product comprises a memory and a processor, and computer instructions stored on the memory and run on the processor, and when the computer instructions are run by the processor, the method of any one of claims 1-5 is completed.

8. A computer-readable storage medium, characterized in that, Computer program product for storing computer instructions which, when executed by a processor, perform the method of any one of claims 1-5.

Images