Heterogeneous Parallel Helmholtz Operator Construction Method for DCU Clusters

Abstract

The present invention discloses a method for constructing a heterogeneous parallel Helmholtz operator for a DCU cluster, belonging to the technical field of computational fluid dynamics. The present invention proposes a method for constructing a heterogeneous parallel Helmholtz operator for a DCU cluster, which mainly includes two parts: block matrix multiplication memory access optimization and task decomposition algorithm design. Compared with the existing design, the present invention no longer performs fine-grained access to the main memory, but instead fully utilizes the on-chip shared memory and registers on the DCU to perform block-parallel prefetching of matrices, significantly reducing memory access overhead. In addition, the present invention enables reasonable allocation of tasks on the DCU, fully utilizing the computing power of the DCU, and the acceleration effect increases accordingly with the increase of the interpolation order. Finally, the load balancing problem between the program CPU and DCU, and between multiple DCUs is solved, thereby improving the scalability of the program. The use of multiple DCUs can further improve the acceleration ratio.

Description

Technical Field

[0001] The present invention relates to the technical field of computational fluid dynamics, and in particular to a method for constructing heterogeneous parallel Helmholtz operators for DCU clusters. Background Art

[0002] In the field of computational fluid dynamics, using the spectral element method to solve the Navier-Stokes equations without turbulence model simplification for direct numerical simulation (DNS) can capture fluid turbulence motion in complex application scenarios. However, the computational complexity of the spectral element method is relatively high. Its main hot function is to construct the Helmholtz operator Ax for the velocity and pressure equations in each solution time step. The pseudo code of the algorithm is as follows:

[0003]

[0004] The large number of small dense matrix multiplication operations generated during the Helmholtz operator construction process consumes the majority of the program's execution time, directly impacting the computational cost of DNS. Parallel optimization based on the size and number of matrices used in the Helmholtz operator construction is key to improving the program's parallel efficiency.

[0005] With the increasing demand for scientific computing, traditional Helmholtz operator solutions based on a single CPU have become unable to meet the requirements of large-scale data processing and high-precision computing. In recent years, heterogeneous computing systems have become an effective way to solve large-scale computing problems due to their high parallel computing capabilities and energy efficiency advantages.

[0006] The Deep Computing Unit (DCU) is a coprocessor product based on the Haiguang Information GPGPU architecture suitable for accelerating computationally intensive tasks. It is equipped with a large number of computing units and high-bandwidth memory. However, the application of existing Helmholtz operator solution methods in the CPU-DCU heterogeneous computing environment still faces many challenges, such as:

[0007] (1) Load balancing between CPU and DCU. Due to the large difference in computing power between CPU and DCU, tasks and data cannot be simply evenly distributed. Instead, they need to be reasonably allocated based on the application characteristics.

[0008] (2) Memory access overhead. Dense matrix multiplication itself is a memory-intensive computing task. How to make full use of the cache, reduce data transmission bottlenecks, reduce memory access time overhead, and improve overall computing efficiency is a technical problem that needs to be solved urgently.

[0009] (3) The computational complexity of a single small dense matrix multiplication task is relatively small. How to schedule tasks on the DCU and make full use of the DCU computing power is the key to improving parallel efficiency.

[0010] Based on the above content, the present invention proposes a method for generating inflow conditions for high-precision simulation of reactor thermal fluids. Summary of the Invention

[0011] The purpose of this invention is to propose a heterogeneous parallel Helmholtz operator construction method for DCU clusters to address the problems of memory-intensive, uneven task distribution, high memory overhead, and insufficient DCU computing power utilization associated with the large number of dense small matrix multiplications generated during Helmholtz operator construction in a DCU environment. This method implements multi-level prefetching of dense small matrix multiplications and rationally distributes computational tasks, reducing memory overhead while fully utilizing DCU computing performance, significantly improving program execution speed.

[0012] In order to achieve the above object, the present invention adopts the following technical solutions:

[0013] The heterogeneous parallel Helmholtz operator construction method for DCU clusters includes the following contents:

[0014] (1) Block matrix multiplication memory access optimization: Based on the block matrix multiplication, when calculating different result sub-matrices Ctile, all sub-matrices Atile and Btile are read into the shared memory space, thereby reducing the number of global memory accesses and lowering the memory access latency;

[0015] For the multiplication of sub-matrices Atile and Btile, registers are introduced to further improve the efficiency of data reuse. Each time a register is used, a row and a column of sub-matrices Atile and Btile are read to calculate a portion of sub-matrix Ctile, achieving memory access optimization.

[0016] Obtain fluid mechanics experimental data, and use RANS, LES and DNS methods to obtain simulation results, and store the obtained experimental data and simulation results in the data source module;

[0017] (2) Task decomposition algorithm design:

[0018] For the matrix generated during the construction of the Helmholtz operator in the DCU environment, the size of the matrix depends on the polynomial interpolation order of the spectral element method;

[0019] The shape of the matrix is ​​divided into N×N×N, N 2 ×N×N and N×N×N 2 The three situations include the following:

[0020] 4) For the N×N×N case: treat a single matrix multiplication as a block and distribute it among the thread groups;

[0021] 5) For N 2×N×N case: In this case, the size of the result matrix C is N 2 ×N; when N < 8, split the matrix into N matrix multiplications of N × N × N size, and use the method in case 1) to perform task decomposition; when N ∈ [8, 24], the minimum number of matrix rows is 64, adjust the block size to 64 * N, and the load of each thread is N elements;

[0022] 6) For N×N×N 2 Situation: Using N 2 ×N-scale matrices have a similar blocking idea. When N is less than 8, the matrix is ​​split into N N×N×N matrix multiplications, and the method in case 1) is used for task decomposition. When N∈[8,24], the number of columns of the result matrix C is much larger than the number of rows. The block size is changed to N*64, and the matrix is ​​divided into multiple blocks. The results in the blocks are distributed to threads by column. Under this blocking shape, the data of each submatrix is ​​no longer distributed continuously in the original matrix. The 64 elements in the same row are at continuous memory addresses. The 64 threads in the thread warp read 64 elements at continuous memory addresses in each round.

[0023] Preferably, the scale of the register in step (1) is tile*tile.

[0024] Preferably, for case 1), each grid generates N blocks during the Helmholtz operator construction process. Since the dimension N of the matrix may not be evenly divisible by the dimension of the thread group, in order to avoid load imbalance within the same thread group, the blocks are decomposed into four parts for calculation. The decomposition process is as follows:

[0025] ① Calculate the number of chunks elementChunkNum and the total number of chunks totalChunkNum generated by a single grid:

[0026] elementChunkNum=N

[0027] totalChunkNum=Nel*elementChunkNum

[0028] ②The global z-axis coordinate of the thread is the global number of the thread group to which the thread belongs, and the number of the block it is responsible for is obtained accordingly:

[0029] gridGroupId=blockDim.z*blockIdx.z+threadIdx.z

[0030] ③ The block is split into four parts and calculated sequentially. Each part is processed by a thread in the thread group. Each thread obtains the starting address of the calculated sub-block Ctile based on the x, y coordinates and sub-block size, which is used to write the calculation results.

[0031] Preferably, the scenario 2) specifically includes the following:

[0032] ① Calculate the number of blocks elementChunkNum and the total number of blocks totalChunkNum generated by a single grid; if there are extra rows so that N 2 If the number cannot be divided evenly by 64, it is also assigned to a thread group for calculation:

[0033]

[0034] totalChunkNum=Nel*elementChunkNum

[0035] ② Obtain the responsible block according to the global number of the thread group. Each thread in the thread group is responsible for calculating one row of elements. Each time an element of matrix A and a row of matrix B are read to achieve data reuse within the block. In terms of memory access merging optimization, the block prefetch method is used for matrix A, and 64*N blocks of size are read for access by one thread warp.

[0036] Compared with the prior art, the present invention provides a method for constructing a heterogeneous parallel Helmholtz operator for a DCU cluster, which has the following beneficial effects:

[0037] (1) This invention solves the memory access overhead problem of dense small matrix multiplication. Instead of fine-grained access to main memory, it fully utilizes the on-chip shared memory and registers on the DCU to perform parallel prefetching of matrices in blocks, significantly reducing memory access overhead.

[0038] (2) This invention rationally distributes tasks on the DCU, fully utilizing the DCU's computing power. The acceleration effect increases with the increase in the interpolation order, and the speedup ratio increases accordingly. The speedup ratio can reach a minimum of 3.61 and a maximum of 30.04.

[0039] (3) This invention solves the load balancing problem between the program CPU and DCU, and between multiple DCUs, improving the scalability of the program. Using multiple DCUs can further improve the speedup ratio; tests have shown that when the interpolation order reaches 24, the speedup ratio can be further improved by 17.9% compared to a single DCU. BRIEF DESCRIPTION OF THE DRAWINGS

[0040] Figure 1 This is a schematic diagram of the 3D thread structure mentioned in Example 1 of the present invention;

[0041] Figure 2 This is a schematic diagram of the N×N matrix thread allocation strategy mentioned in Example 1 of the present invention;

[0042] Figure 3N is the N proposed in Example 1 of the present invention 2 ×N matrix thread allocation strategy diagram;

[0043] Figure 4 is the N×N proposed in Example 1 of the present invention 2 Schematic diagram of matrix thread allocation strategy. DETAILED DESCRIPTION

[0044] The technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, rather than all the embodiments.

[0045] The present invention proposes a method for constructing a heterogeneous parallel Helmholtz operator for a DCU cluster, which includes the following contents:

[0046] (1) Block matrix multiplication memory access optimization:

[0047] Matrix multiplication is a memory-intensive task, and the main optimization direction is to reduce the number of global memory accesses. Based on the block matrix multiplication, the different result sub-matrices C tile The calculation may need to access the same sub-matrix A tile or B tile , that is, there is data reuse between multiple threads. Therefore, firstly, the sub-matrix A tile ,B tile All data is read into the shared memory space. In this way, when a thread in a single thread group accesses a continuous memory address, the DCU memory controller only needs to read all the data once, greatly reducing the memory access latency.

[0048] For the submatrix A tile and B tile The multiplication of the sub-matrix A is further improved by introducing registers. tile and B tile A row and a column are used to calculate the submatrix C tile Memory access optimization can be achieved by using a portion of the register size of tile*tile.

[0049] (2) Task decomposition algorithm design:

[0050] The calculations generated during the construction of the Helmholtz operator are mainly based on a large number of small-scale matrix multiplications, and the size of the small-scale matrix depends on the polynomial interpolation order of the spectral element method. The order range used in the program is N∈[4,24], and the resulting matrix shape can be divided into N×N×N, N 2 ×N×N and N×N×N 2Three situations. The present invention performs task decomposition optimization for the three situations respectively, which are described in detail as follows.

[0051] (1) For the N×N×N case: This form of matrix-matrix multiplication has more numbers than other types, but the amount of computation involved in a single matrix multiplication is small. Therefore, consider treating a single matrix multiplication as a block and distributing it among thread groups. The threads in the block are organized into a thread bundle with 64 threads and run simultaneously. The thread bundle is composed of multiple thread groups, so the number of threads in the thread group needs to be divisible by 64. The optional thread group sizes are 2*2, 4*4, and 8*8. During the construction of the Helmholtz operator, each grid can generate N blocks. Since the dimension N of the matrix may not be divisible by the dimension of the thread group, in order to avoid load imbalance within the same thread group, consider decomposing the block into four parts for calculation. The decomposition process is as follows:

[0052] ① Calculate the number of chunks elementChunkNum and the total number of chunks totalChunkNum generated by a single grid:

[0053] elementChunkNum=N

[0054] totalChunkNum=Nel*elementChunkNum

[0055] ②The global z-axis coordinate of the thread is the global number of the thread group to which the thread belongs, which can be used to obtain the number of the block it is responsible for:

[0056] gridGroupId=blockDim.z*blockIdx.z+threadIdx.z

[0057] ③ The block is split into four parts and calculated sequentially. Each part is processed by a thread in the thread group. Each thread obtains the starting address of the calculated sub-block Ctile based on the x, y coordinates and sub-block size, and uses it to write the calculation results.

[0058] (2) For N 2 ×N×N case: This form of matrix multiplication is computationally intensive, and the resulting matrix C is of size N 2 ×N. When N<8, the matrix is ​​split into N matrix multiplications of N×N×N size, and the task decomposition is performed using the method in case (1). When N∈[8,24], the minimum number of matrix rows can reach 64, which is enough to occupy a thread warp for parallelization. Therefore, the block size is adjusted to 64*N, and the load of each thread is N elements. The specific process is as follows:

[0059] ① Calculate the number of blocks elementChunkNum and the total number of blocks totalChunkNum generated by a single grid. 2 If the number cannot be divided evenly by 64, that is, there are extra rows, then they are also assigned to a thread group for calculation:

[0060]

[0061] totalChunkNum=Nel*elementChunkNum

[0062] ② The thread group's global number determines the block it is responsible for. Each thread in the thread group is responsible for computing one row of elements. Each read of one element of matrix A and one row of matrix B reuses data within the block. To optimize memory merging, block prefetching can be used for matrix A, reading 64*N blocks for access by a single warp.

[0063] (3) For N×N×N 2 Case: The matrix multiplication of this shape is similar to the previous one, so the same matrix as N is used. 2 The block idea is similar to that of a matrix of size N×N. When N<8, the matrix is ​​split into N matrix multiplications of size N×N×N, and the task decomposition is performed using the method of case (1). When N∈[8,24], since the number of columns of the result matrix C is much larger than the number of rows, the block size is changed to N×64, the matrix is ​​divided into multiple blocks, and the results in the blocks are distributed to the threads by column. The same idea as the previous matrix is ​​used to implement memory access optimization. However, the difference is that in this block shape, the data of each submatrix is ​​no longer distributed continuously in the original matrix, but the 64 elements of the same row happen to be at continuous memory addresses. The 64 threads in the thread bundle read exactly 64 elements at continuous memory addresses in each round, so the optimization of memory access merging can be achieved without shared memory.

[0064] The heterogeneous parallel Helmholtz operator construction method for DCU clusters proposed in the present invention is described below with reference to the accompanying drawings and specific examples, which specifically include the following contents.

[0065] Example 1:

[0066] See also Figure 1 In this invention, the DCU organizes threads into blocks. The thread number is a 3D vector that describes the thread's position within the entire thread bundle. This allows you to set the thread organization structure when executing kernel functions. The program can use blockDim.x, blockDim.y, and blockDim.z to obtain block size information. Regarding memory access optimization for block matrix multiplication, the optimized matrix multiplication instruction can be divided into four steps:

[0067] (1) Load data from global memory.

[0068] (2) Store the data in shared memory.

[0069] (3) Load data from shared memory to registers.

[0070] (4) Perform calculation operations.

[0071] The algorithm flow for shared memory loading and register loading is as follows:

[0072]

[0073] This embodiment uses a parallel transfer strategy. Since the on-chip shared memory space is large, M tile *K-sized A matrix blocks and K*N tile The B matrices of different sizes are moved into shared memory in blocks at the same time. Here we use A shared Represents the A matrix portion in the on-chip shared memory. Each thread calculates its own blockThreadId to determine the starting position for transfer, and transfers a piece of data every blockSize interval. This ensures that the data addresses transferred by each thread are continuous in each round of transfer. The specific algorithm flow is as follows:

[0074]

[0075] Here A reg This represents the portion of matrix A stored in the register. The register capacity is small, so before performing each vector multiplication, a row of matrix A and a column of matrix B in shared memory are moved into the register.

[0076] For the decomposition of the N×N×N matrix multiplication task, each thread calculates the starting address of the sub-block it is responsible for based on the x and y coordinates. First, the memory address offset is calculated based on the four-part division method:

[0077] rowOffset=threadIdx.y*tileN*N

[0078] colOffset = threadIdx.x*tileN

[0079] resRosOffset=(blockDim.y*tileN+threadIdx.y)*N

[0080] resColOffset=blockDim.x*tileN+threadIdx.x

[0081] The four offset combinations allow us to calculate the starting memory addresses of the four parts. We then cyclically allocate matrix blocks across thread groups until all calculations are complete. The number of blocks calculated simultaneously at any one time is equal to the number of global thread groups, i.e., the global z-axis coordinate scale is gridDim.z * blockDim.z. The specific algorithm flow is as follows:

[0082]

[0083]

[0084] The four matrix blocks generated by this example are:

[0085] (1) The upper left block is a square matrix of size tileN*blockDim.x, which can be evenly distributed among the threads in the thread group. The size of the sub-block is tileN*tileN.

[0086] (2) The upper right block is a thin and long matrix of size res*tileN*blockDim.x. The size of the sub-block is tileN*1. The number of threads in the thread group with size res*blockDim.x participates in the operation.

[0087] (3) The lower left block is a flat wide matrix of size res*tileN*blockDim.y. The size of the sub-block is 1*tileN. There are res*blockDim.y threads in the thread group participating in the operation.

[0088] (4) The lower right block is a square matrix of size res*res, which contains a small amount of calculation, so the calculation of each element is assigned to a thread in the thread group for execution.

[0089] N 2 ×N×N and N×N×N 2 The algorithm flow of the type is similar, and the only difference is the block prefetching when calculating the block multiplication at the end. 2 ×N×N type is used as an example to give the task decomposition algorithm flow when N∈[8,24]:

[0090]

[0091]

[0092] This example uses a 3D thread structure, distinguishing thread groups by z coordinate and cyclically allocating computations to blocks within a block. Since the number of matrix rows in this case is greater than 64, fully utilizing a single warp, the matrix is ​​no longer partitioned between thread groups. For each thread group, the computation flow is as follows:

[0093] (1) Get the offset chunkOffset of the responsible block according to the global number of the thread group.

[0094] chunkOffset=elementId*N 3 +chunkId*chunkSize

[0095] (2) Get the responsible block A according to the offset chunk ,B chunk ,C chunk .

[0096]

[0097] (3) Each thread uses rowOffset to offset the A it is responsible for chunk ,B chunk The blocks are matched with the corresponding shared memory addresses to form a sub-matrix A tile ,B tile ,C tile , A and B are used for pre-fetching, and C is used for pre-storing results.

[0098]

[0099] (4) Calculate the corresponding result sub-matrix C in shared memory tile .

[0100] The above description is only a preferred specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. For ordinary technicians in this technical field, several improvements and modifications can be made without departing from the principles described in the present invention. These improvements and modifications should also be regarded as the scope of protection of the present invention.

Claims

1. A method for constructing heterogeneous parallel Helmholtz operators for DCU clusters, characterized by: Includes the following: (1) Block matrix multiplication memory access optimization: Based on the block matrix multiplication, when calculating different result sub-matrices Ctile, all sub-matrices Atile and Btile are read into the shared memory space, thereby reducing the number of global memory accesses and reducing memory access latency; For the multiplication of sub-matrices Atile and Btile, registers are introduced to further improve the efficiency of data reuse. Each time a register is used, a row and a column of sub-matrices Atile and Btile are read to calculate a portion of sub-matrix Ctile, achieving memory access optimization. (2) Task decomposition algorithm design: For the matrix generated during the construction of the Helmholtz operator in the DCU environment, the size of the matrix depends on the polynomial interpolation order of the spectral element method; The shape of the matrix is ​​divided into N × N × N 、 N 2 × N × N and N × N × N 2 The three situations include the following: 1) Targeting N × N × N Case: A single matrix multiplication is treated as a block and distributed among thread groups; 2) Targeting N 2 × N × N Case: In this case, the size of the result matrix C is N 2 × N ; when N When <8, split the matrix into N indivual N × N × N Scale matrix multiplication, the method in case 1) is used to decompose the task; when N ∈[8,24], the minimum number of matrix rows is 64, and the block size is adjusted to 64* N , the load of each thread is N elements; 3) Targeting N × N × N 2 Situation: When N When <8, split the matrix into N N × N × N Scale matrix multiplication, the method in case 1) is used to decompose the task; when N ∈[8,24], the transformation block size is N *64, divide the matrix into multiple blocks and distribute the results in the blocks to threads by column; In this block shape, the data of each sub-matrix is ​​no longer distributed continuously in the original matrix. The 64 elements in the same row are at consecutive memory addresses, and the 64 threads in the thread warp read 64 elements at consecutive memory addresses in each round.

2. The method for constructing a heterogeneous parallel Helmholtz operator for a DCU cluster according to claim 1, characterized in that: The size of the register in step (1) is tile * tile.

3. The method for constructing a heterogeneous parallel Helmholtz operator for a DCU cluster according to claim 1, characterized in that: For case 1), each grid generates N blocks, due to the dimension of the matrix N Sometimes the dimension of a thread group cannot be divided evenly. To avoid load imbalance within the same thread group, the block is decomposed into four parts for calculation. The decomposition process is as follows: Calculate the number of blocks generated by a single grid elementChunkNum and the total number of blocks totalChunkNum : elementChunkNum = N totalChunkNum = Nel * elementChunkNum The global z-axis coordinate of a thread is the global number of the thread group to which the thread belongs, and the number of the block it is responsible for is obtained accordingly: gridGroup Id = blockDim.z * blockIdx.z + threadIdx.z The block is split into four parts and calculated sequentially. Each part is processed by a thread in the thread group. Each thread obtains the starting address of the calculated sub-block Ctile based on the x, y coordinates and sub-block size, and is used to write the calculation results.

4. The method for constructing a heterogeneous parallel Helmholtz operator for a DCU cluster according to claim 1, characterized in that: Scenario 2) specifically includes the following: Calculate the number of blocks generated by a single grid elementChunkNum and the total number of blocks totalChunkNum ; If there are extra rows such that N 2 If the number cannot be divided evenly by 64, it is also assigned to a thread group for calculation: The thread group obtains the block it is responsible for according to its global number. Each thread in the thread group is responsible for calculating one row of elements. Each time an element of matrix A and a row of matrix B are read to achieve data reuse within the block. In terms of memory access merging optimization, block prefetching is used for the A matrix, and blocks of size 64*N are read for access by a thread warp.

Images