DCU cluster-oriented heterogeneous parallel Helmholtz operator construction method

Abstract

The invention discloses a heterogeneous parallel Helmholtz operator construction method for a DCU cluster, and belongs to the technical field of computational fluid mechanics. The invention provides a heterogeneous parallel Helmholtz operator construction method oriented to a DCU cluster, which mainly comprises two parts of contents, namely block matrix multiplication memory access optimization and task decomposition algorithm design, compared with the existing design, fine-grained access is not carried out on a main memory any more, an on-chip shared memory and a register on the DCU are fully utilized, and the design efficiency is improved. Block parallel prefetching is carried out on the matrix, so that the memory access overhead is obviously reduced; in addition, tasks are reasonably distributed on the DCU, the computing power of the DCU is fully utilized, and the acceleration ratio of the acceleration effect is correspondingly increased along with the increase of the interpolation order; and finally, the problem of load balancing between the program CPU and the DCU and among multiple DCUs is solved, the expandability of the program is improved, and the acceleration ratio can be further improved by using the multiple DCUs.

Description

Technical Field

[0001] The present invention relates to the technical field of computational fluid dynamics, and in particular to a method for constructing a heterogeneous parallel Helmholtz operator for a DCU cluster. 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 the turbulent motion of fluids in complex application scenarios. However, the computational complexity of the spectral element method is relatively high. Its main hotspot function is to construct the Helmholtz operator Ax of the velocity and pressure equations in each solution time step. The algorithm pseudo code is as follows:

[0003]

[0004] The large number of small dense matrix multiplication operations generated during the Helmholtz operator construction process takes up most of the program running time, and its computational efficiency directly affects the computational cost of DNS. Parallel optimization based on the size and number characteristics of the matrix in the Helmholtz operator construction is the key to improving the parallel efficiency of the program.

[0005] With the increasing demand for scientific computing, the traditional Helmholtz operator solution method based on a single CPU has been 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] DCU (Deep Computing Unit) is a coprocessor product of Haiguang Information's GPGPU architecture suitable for accelerating computationally intensive tasks, 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. Since the computing power of CPU and DCU is very different, tasks and data cannot be simply evenly distributed when they are divided. Instead, they need to be reasonably distributed according to 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 the present invention is to propose a heterogeneous parallel Helmholtz operator construction method for DCU clusters to solve the problems of intensive memory access and uneven task distribution, high memory access overhead and insufficient utilization of DCU computing power generated in the process of Helmholtz operator construction in the DCU environment. The present invention realizes multi-level prefetching of dense small matrix multiplication and reasonable allocation of computing tasks, reduces memory access overhead and makes full use of DCU computing performance, which can greatly improve the program running 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 reducing 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 the register is used, a row and a column of sub-matrices Atile and Btile are read to calculate a part of sub-matrix Ctile, thereby optimizing memory access;

[0016] S2. Obtain fluid mechanics experimental data, and use RANS, LES and DNS methods to obtain simulation calculation 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 shapes of the matrix are 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 N2 ×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: Adopting and N 2 ×N-scale matrix, when N < 8, split the matrix into N N × N × N-scale matrix multiplications, and use the method in case 1) to perform task decomposition; when N ∈ [8, 24], the number of columns of the result matrix C is much larger than the number of rows, and the block size is changed to N * 64, the matrix is ​​divided into multiple blocks, and the results in the blocks are distributed to threads by column; under this block shape, the data of each sub-matrix is ​​no longer distributed continuously in the original matrix, and 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 S1 is tile*tile.

[0024] Preferably, for situation 1), each grid generates N blocks during the Helmholtz operator construction process. Since the dimension N of the matrix cannot divide the dimension of the thread group, in order to avoid load imbalance in the same thread group, the block is 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 for sequential calculation. Each part is processed by a thread in the thread group. Each thread obtains the starting address of the calculated sub-block Ctile according to the x, y coordinates and sub-block size to write the calculation results.

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

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

[0033]

[0034] totalChunkNum=Nel*elementChunkNum

[0035] ②According to the global number of the thread group, the responsible block is obtained. Each thread in the thread group is responsible for the calculation of one row of elements. One element of the A matrix and one row of the B matrix are read each time to realize data reuse within the block. In terms of memory access merging optimization, the A matrix is ​​pre-fetched in blocks, and blocks of 64*N size are read for access by one thread warp.

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

[0037] (1) The present invention solves the memory access overhead problem of dense small matrix multiplication. Instead of fine-grained access to the main memory, the on-chip shared memory and registers on the DCU are fully utilized to perform block-by-block parallel pre-fetching of the matrix, significantly reducing the memory access overhead.

[0038] (2) The present invention makes reasonable allocation of tasks on the DCU and fully utilizes the computing power of the DCU. The acceleration effect increases with the increase of the interpolation order. The acceleration ratio can reach a minimum of 3.61 and a maximum of 30.04.

[0039] (3) The present invention solves the problem of load balancing between the program CPU and DCU, and between multiple DCUs, and improves the scalability of the program. The use of multiple DCUs can further improve the speedup ratio; after testing, when the interpolation order reaches 24, the speedup ratio can be further improved by 17.9% compared with 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 3 is 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 Embodiment 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 described clearly and completely 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 block matrix multiplication, the calculation of different result sub-matrices Ctile may require access to the same sub-matrix Atile or Btile, that is, there is data reuse between multiple threads. Therefore, all sub-matrices Atile and Btile are first 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 once to obtain all data, greatly reducing memory access latency.

[0048] 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 to read a row and a column of sub-matrices Atile and Btile, which are used to calculate a part of sub-matrix Ctile, memory access optimization can be achieved. The register size used is 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 generated 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 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 number of matrix-matrix multiplications than other types, but the amount of computation contained 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 in 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 according to the x, y coordinates and sub-block size, which is used to write the calculation result.

[0058] (2) For N 2 ×N×N case: This form of matrix multiplication is computationally intensive, and the size of the resulting matrix C is 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 parallel operation. Therefore, the size of the block 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 chunks elementChunkNum and the total number of chunks totalChunkNum generated by a single grid. 2 If the number cannot be divided by 64, that is, there are extra rows, they are also assigned to a thread group for calculation:

[0060]

[0061] totalChunkNum=Nel*elementChunkNum

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

[0063] (3) For N×N×N 2 Case: Matrix multiplication of this shape is similar to the previous one, so the same matrix as N is used. 2 The same idea of ​​block partitioning as that of a matrix of size N×N is used. When N<8, the matrix is ​​split into N matrix multiplications of size N×N×N, and the task is decomposed using the method in 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 threads by column, and 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 warp read the 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 by the present invention is described below with reference to the accompanying drawings and specific examples, which specifically include the following contents.

[0065] Embodiment 1:

[0066] See also Figure 1 In the present invention, DCU organizes threads in a block manner. The thread number is a 3D vector that describes the position of the thread in the entire thread bundle. The thread organization structure can be set when executing the kernel function. The program can use blockDim.x, blockDim.y, and blockDim.z to obtain the size information of the block. In terms of memory access optimization for block matrix multiplication, the matrix multiplication instruction after memory access optimization 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 into 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 It represents the A matrix in the on-chip shared memory. Each thread calculates its own blockThreadId to determine the starting position of the transfer, and transfers a piece of data at each interval of blockSize. In this way, 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 It represents the part of matrix A stored in the register. The register capacity is small, and before each vector in the matrix multiplication, one row of matrix A and one column of matrix B in the 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 division method of the four parts:

[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 starting memory addresses of the four parts can be calculated by combining the four offsets. Then the matrix blocks are allocated cyclically among the thread groups until all calculations are completed. The number of blocks calculated simultaneously each time is the number of global thread groups, that is, the global z-axis coordinate scale 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 threads of res*blockDim.x in the thread group participate in the calculation.

[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 calculation.

[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 two types is similar, and the only difference is the block prefetching when calculating the block multiplication at the end. 2 ×N×N type is taken as an example to give the task decomposition algorithm flow when N∈[8,24]:

[0090]

[0091] This example uses a 3D thread structure, uses the z coordinate to distinguish different thread groups, and cyclically allocates blocks within the block for calculation. Since the number of matrix rows in this case will be greater than 64, which is enough to fully utilize a thread warp, the matrix is ​​no longer divided into blocks between thread groups. For each thread group, the calculation process is as follows:

[0092] (1) Obtain the offset chunkOffset of the chunk for which the thread group is responsible based on its global number.

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

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

[0095]

[0096] (3) Each thread uses rowOffset to offset the A it is responsible for chunk ,B chunk The blocks and corresponding shared memory addresses A tile ,C tile Correspondingly, A is used for pre-fetching and C is used for pre-storing results.

[0097]

[0098] (4) Calculate the corresponding result C in shared memory tile .

[0099] The above description is only a preferred specific implementation manner of the present invention, but the protection scope 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, and these improvements and modifications should also be regarded as the protection scope of the present invention.

Claims

1. A method for constructing heterogeneous parallel Helmholtz operators for DCU clusters, characterized in that: 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 the register is used, a row and a column of sub-matrices Atile and Btile are read to calculate a part of sub-matrix Ctile, thereby optimizing memory access; (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 shapes of the matrix are N×N×N, N 2 ×N×N and N×N×N 2 The three situations include the following: 1) For the N×N×N case: treat a single matrix multiplication as a block and distribute it among the thread groups; 2) 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; 3) For N×N×N 2 Situation: Adopting and N 2 ×N-scale matrix has a similar block idea. When N < 8, the matrix is ​​split into N N×N×N-scale 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, the matrix is ​​divided into multiple blocks, and the results in the blocks are distributed 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 continuous memory addresses, and the 64 threads in the thread warp read 64 elements at continuous 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 S1 is tile*tile.

3. The method for constructing a heterogeneous parallel Helmholtz operator for a DCU cluster according to claim 1, characterized in that: In case 1), each grid generates N blocks during the construction of the Helmholtz operator. Since the dimension N of the matrix cannot divide the dimension of the thread group, in order to avoid load imbalance in the same thread group, the block is decomposed into four parts for calculation. The decomposition process is as follows: ① Calculate the number of chunks elementChunkNum and the total number of chunks totalChunkNum generated by a single grid: elementChunkNum=N totalChunkNum=Nel*elementChunkNum ②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: gridGroupId=blockDim.z*blockIdx.z+threadIdx.z ③ The block is split into four parts for sequential calculation. Each part is processed by a thread in the thread group. Each thread obtains the starting address of the calculated sub-block Ctile according to the x, y coordinates and sub-block size 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: The situation 2) specifically includes the following contents: ① Calculate the number of chunks elementChunkNum and the total number of chunks totalChunkNum generated by a single grid; if there are extra rows, N 2 If it cannot be divided evenly by 64, it is also assigned to a thread group for calculation: totalChunkNum=Nel*elementChunkNum ②According to the global number of the thread group, the block it is responsible for is obtained. Each thread in the thread group is responsible for calculating one row of elements. Each time an element of the A matrix and a row of the B matrix are read to achieve data reuse within the block. In terms of memory access merging optimization, block prefetching is used for the A matrix to read blocks of size 64*N for access by a thread warp.

Images