Method for accelerating multi-dimensional fast fourier transform based on matrix operation

Abstract

The application relates to a multi-dimensional fast Fourier transform acceleration method based on matrix operation in the field of computing power optimization. The method comprises the following steps: S1, applying for a computing resource, copying computing data to device global memory on the host side, generating a DFT matrix and an input data matrix; S2, reading the DFT matrix data and the input data matrix from the global memory to a register and a shared memory; S3, using multi-thread parallelism to complete the calculation task of matrix element multiplication of the input data matrix and a rotation factor matrix, and obtaining intermediate calculation result data; S4, calling a matrix operation unit to perform matrix multiplication on the intermediate calculation result data obtained in the step S3 and the DFT matrix, and outputting final calculation result data; and S5, copying the final calculation result data from the device global memory to the host memory on the host side, and releasing the computing resource. The application adopts matrix operation to calculate multi-dimensional FFT, and uses a high-performance matrix operation unit of a GPU to accelerate the calculation, so that the FFT calculation efficiency on the GPU is improved.

Description

Technical Field

[0001] This invention relates to the fields of high-performance computing and computing power optimization, specifically to a method for accelerating multidimensional fast Fourier transform based on matrix operations. Background Technology

[0002] The Fast Fourier Transform (FFT) algorithm is a fast algorithm for calculating the Discrete Fourier Transform (DFT) sequence of a given sequence. FFT is one of the most important fundamental algorithms in scientific computing, widely used in important scientific research areas such as molecular dynamics calculations, electronic structure calculations, and digital signal processing. FFT also has significant applications in deep learning; compared to traditional convolution algorithms, FFT-based convolution algorithms show better performance in network training. Due to the wide range of applications of FFT, parallel optimization of the FFT algorithm has become a research hotspot in the field of high-performance computing. GPUs, due to their powerful parallel computing capabilities, have become the mainstream computing hardware in artificial intelligence and scientific computing. Therefore, optimizing and accelerating the FFT algorithm based on the characteristics of existing GPU hardware is of great significance.

[0003] To accelerate matrix multiplication, a common operation in artificial intelligence computing, NVIDIA designed the Tenor Core in its GPUs to accelerate matrix computation at the hardware level. AMD subsequently designed its counterpart, the MatrixCore. Hardware computing power needs to be accessed through software; inadequate software optimization and adaptation can waste hardware computing power and resources, increasing computational costs. To fully utilize the GPU's computing power and improve the efficiency of FFT computation on the GPU, it is necessary to optimize traditional FFT algorithms based on the characteristics of the new hardware.

[0004] Currently, there are many mature FFT computing libraries available for GPUs, such as NVIDIA's cuFFT, AMD's rocFFT, and the cross-platform Vulkan FFT. However, these libraries do not fully utilize the computing power of modern GPUs, and there is still significant room for optimization in GPU-based FFT algorithms. Summary of the Invention

[0005] To address the technical problems existing in the prior art, this invention proposes a method for accelerating FFT (Fast Fourier Transform) calculations using high-performance matrix operation units on existing GPUs. This method can improve the computational efficiency of FFT, has the characteristics of strong applicability and fast calculation speed, and has good practical value.

[0006] This invention can be achieved by adopting the following technical solutions:

[0007] The method for accelerating the multidimensional fast Fourier transform based on matrix operations includes the following steps:

[0008] S1. Request computing resources to compute an n-dimensional FFT. Copy the computation data to the device's global memory on the host side, generate the DFT matrix and the input data matrix, and initialize the variable c, setting c = n.

[0009] S2. Start the computing kernel and read the DFT matrix data and input data matrix from global memory into registers and shared memory;

[0010] S3. Use multithreading to perform the matrix element multiplication task between the input data matrix and the rotation factor matrix in parallel, and obtain intermediate calculation result data.

[0011] S4. Call the matrix operation unit to perform matrix multiplication between the intermediate calculation result data obtained in step S3 and the DFT matrix to obtain the matrix operation result; assign c to c-1; if the value of c is not 0, transpose the matrix operation result and jump to step S2; if the value of c is 0, output the final calculation result data.

[0012] S5. The host copies the final calculation result data from the device's global memory to the host memory, releasing computing resources.

[0013] Specifically, step S1 includes: analyzing the input sequence on the host side, generating a corresponding DFT matrix on the host side by the CPU according to the length of the input sequence, and copying the corresponding DFT matrix to the device's global memory;

[0014] Arrange the input sequence into an input data matrix in row-major order.

[0015] Specifically, step S2 includes: if the input data matrix cannot be placed into the shared memory at once, the input data matrix is ​​further decomposed and rearranged so that the row length of the data matrix is ​​less than or equal to the size of the shared memory.

[0016] Specifically, step S3 includes: distributing the task of multiplying the input data matrix and the rotation factor matrix into each thread on an equal basis, and then each thread reads the corresponding data from shared memory or global memory;

[0017] Calculate the data in the rotation factor matrix based on the position of the input sequence data in the input data matrix;

[0018] The input data matrix and the rotation factor matrix are multiplied element-wise using multiple threads in parallel. The intermediate calculation results are then written back to shared memory.

[0019] Specifically, the step of calculating the corresponding data in the rotation factor matrix based on the position of the input sequence data in the input data matrix includes: calculating the row and column numbers of the elements storing the input sequence data in the input data matrix using a GPU thread, and finding the corresponding elements in the rotation factor matrix based on the row and column numbers.

[0020] Specifically, step S4 includes:

[0021] The intermediate calculation results and DFT matrix are loaded into the registers of each thread;

[0022] The matrix operation unit in the GPU is used to perform matrix multiplication between the intermediate calculation result data and the DFT matrix, and the matrix operation result is output. The value of c is assigned to c-1.

[0023] If the value of c is not 0, transpose the matrix operation result and jump to step S2; if the value of c is 0, output the final calculation result.

[0024] Compared with the prior art, the present invention has the following advantages and beneficial effects:

[0025] This invention provides a method for accelerating multidimensional fast Fourier transform based on matrix operations. Compared with traditional methods, this method uses matrix operations to calculate multidimensional FFT, uses multithreading to perform the matrix element multiplication of the input data matrix and the rotation factor matrix in parallel, and utilizes the high-performance matrix operation unit on existing GPUs for computational acceleration, thereby improving the FFT computation efficiency on GPUs and providing computing power support for applications such as deep learning and scientific computing. Attached Figure Description

[0026] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on the structures shown in these drawings without creative effort.

[0027] Figure 1 This is a flowchart of the multidimensional fast Fourier transform acceleration method based on matrix operations in an embodiment of the present invention;

[0028] Figure 2 This is a distribution diagram of the computational data in each thread in an embodiment of the present invention;

[0029] Figure 3 This is a schematic diagram of the matrix multiplication calculation process of Matrix Core in an embodiment of the present invention;

[0030] Figure 4This is a schematic diagram of the calculation process of the three-dimensional Fourier transform in an embodiment of the present invention;

[0031] Figure 5 This is a graph showing the speedup effect of the multidimensional fast Fourier transform acceleration method in this embodiment of the invention compared to other FFT calculation libraries. Detailed Implementation

[0032] The technical solution of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments, and the implementation of the present invention is not limited thereto. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0033] The Discrete Fourier Transform (DFT) converts a finite discrete signal in the time domain to the frequency domain according to the following formula: The DFT transform formula is:

[0034]

[0035] Where x[n] represents an element of the discrete-time signal, and X[k] represents an element of the Fourier transform sequence. The time complexity of this algorithm is O(N²). In 1963, JW Cooley and JW Tukey proposed the Cooley-Tukey algorithm, which uses a divide-and-conquer strategy to reduce the time complexity of the DFT from O(N²) to O(N²). 2 The computational complexity is reduced to O(Nlog₂N), significantly improving efficiency. The Cooley-Tukey Fast Fourier Transform algorithm mainly consists of the following three steps:

[0036] • Divide the original input sequence of length N into N1 subsequences of length N2 = N / N1.

[0037] • Use the Cooley-Tukey Fast Fourier Transform algorithm to recursively compute the Fourier transform sequence of a subsequence of length N2.

[0038] • Merge the N1 Fourier transform subsequences into a Fourier transform sequence of length N to obtain the final transform result.

[0039] The matrix-form Fourier transform algorithm is developed based on the Cooley-Tukey Fast Fourier Transform algorithm. The original input sequence of length N is divided into N1 subsequences of length N2 = N / N1. These subsequences are then arranged into an N1 × N2 matrix. The correspondence between elements in the original sequence and elements in the subsequences can be expressed by the following formula:

[0040] x m[s] = x[sN1+m]

[0041] Where x m [s] represents the m-th subsequence. Assume the Fourier transform sequence of this subsequence is X. m Combining the DFT transform formula, let:

[0042]

[0043] Let N be a DFT (Discrete Fourier Transform) matrix of size N1×N1, where i is the imaginary unit.

[0044]

[0045] Let N be a rotation factor matrix of size N1×N2, where

[0046]

[0047] Among them, X in Let N be a matrix of size N1×N2 formed by arranging the original input sequence in row-major order, where N = N1×N2 and N is the length of the original input sequence. This represents the element in the N1-1th row (index starts from zero) and N2-1th column of the matrix.

[0048] The matrix form of the Fourier transform algorithm can then be expressed as follows:

[0049]

[0050] Here, ⊙ denotes matrix element multiplication. Since the Fourier transform sequence of a sequence of length 1 is the sequence itself, the Fourier transform sequence of the original sequence can be obtained by repeating the above process. This time it was obtained.

[0051] By using the matrix operation units in modern GPUs to accelerate FFT calculations, the computational efficiency of FFT can be greatly improved, thereby providing underlying computing power support for scientific computing and deep learning.

[0052] Example 1:

[0053] The multidimensional fast Fourier transform acceleration method based on matrix operations described in this invention uses matrix operations to calculate multidimensional FFT and utilizes high-performance matrix operation units (such as Tensor Cores and Matrix Cores) on existing GPUs for computational acceleration to calculate the discrete Fourier transform (DFT) sequence of a certain sequence. This improves the computational efficiency of the FFT fast Fourier transform algorithm on GPUs and provides computing power support for applications such as deep learning and scientific computing.

[0054] The proposed method for accelerating the multidimensional fast Fourier transform based on matrix operations in this invention, such as... Figure 1 As shown, it includes the following steps:

[0055] S1. Request computing resources to compute an n-dimensional FFT. The host copies the computation data to the device's global memory and generates the necessary pre-component data. Initialize the variable c, setting c = n. Specifically, the pre-component data required for computation includes the DFT matrix and the input data matrix.

[0056] S11. The input sequence is analyzed on the host side, and the corresponding DFT matrix is ​​generated by the CPU on the host side according to the length of the input sequence. For example, to calculate a sequence of length 256, a 16x16 DFT matrix is ​​generated, and to calculate a sequence of length 1024, a 32x32 DFT matrix is ​​generated. The DFT matrix is ​​generated by the CPU on the host side according to the FFT calculation principle, and then copied to the device's global memory.

[0057] Since the 16×16 DFT matrix is ​​identical for all blocks, it is pre-generated on the host machine and then copied to the device's global memory for direct reading by each block. The DFT matrix is ​​as follows:

[0058]

[0059] Among them, F 16 It is a 16×16 DFT matrix, where

[0060] S12. Arrange the input sequences into a row-matrix. For each input sequence of length 256, its Fourier transform sequence needs to be calculated, requiring the input sequences to be arranged into a 16×16 input data matrix in row-matrix. The input data matrix is ​​denoted as:

[0061]

[0062] Among them, X in Given the input data matrix, X 15

[15] represents the element located in the 15th row and 15th column of the matrix.

[0063] S2. Start the computing kernel, read the DFT matrix data and input data matrix from global memory and put them into registers and shared memory.

[0064] This step of the data transfer process considers that the input data matrix can be put into shared memory at once. However, shared memory is a high-speed storage, but its capacity is generally small. For example, the shared memory size of the MI250 computing card is 64KB, which can hold up to 8192 single-precision complex numbers. If the length of the calculated sequence exceeds 8192, the shared memory cannot store all the intermediate data at once, and the intermediate calculation results need to be temporarily stored in global memory.

[0065] If the input data matrix cannot fit into shared memory at once, it is further decomposed and rearranged so that the row length is less than or equal to the shared memory size. This rearrangement is considered as K1×K2. If the length K2 is still greater than the shared memory size, the decomposition continues until the length is still less than or equal to the shared memory size, allowing the shared memory to hold the data. Then, the elements of the matrix are multiplied by a factor, transposed, and a Fourier transform is performed on each row as a sequence of length k. The results are then concatenated to obtain the Fourier transform sequence of the original sequence. The detailed calculation process is described below:

[0066] Suppose we need to calculate the Fourier transform of a sequence of length N (N is very large), we need to decompose the length N into K1×K2, and then perform the following four steps:

[0067] Treat the input data as a K1×K2 matrix, perform a Fourier transform on each row to obtain K1 Fourier transform sequences of length K2.

[0068] Let each element of the matrix above be X. mn Multiply each element by

[0069] • Transpose the K1×K2 matrix to obtain a matrix of size K2×K1.

[0070] • Perform Fourier transforms on each row in parallel, and then concatenate the data to obtain a Fourier transform sequence of length N.

[0071] The algorithm above decomposes a Fourier transform sequence of length N into two batches of Fourier transforms of lengths K1 and K2, thus making the length of the Fourier transform that the kernel can solve no longer limited by shared memory.

[0072] Next, the host machine calls the device code to start the computation kernel and begin running the GPU program, reading the computation data and DFT matrix into device registers or shared memory. Since the computation data is already in global memory from the previous step, reading the computation data and DFT matrix into device registers or shared memory speeds up data retrieval efficiency and reduces memory access time. Because each block computes a 256-bit FFT, each block loads its own 256 computation data and DFT matrix into shared memory and thread registers. The loading process is parallel; for example, if each block has 64 threads, then each thread loads 4 data points into its own registers.

[0073] Assuming we are computing a 256×256×256 three-dimensional FFT, we first compute the FFT of the first dimension with a length of 256, meaning we compute 256×256 FFTs of length 256 each. This 256×256 represents the computation batch size; each block of the GPU computes one 256-length FFT, resulting in a total of 256×256 blocks, enabling parallel computation via the GPU. The corresponding DFT matrix is ​​copied to the device's global memory, and each block retrieves the relevant data from this memory for computation.

[0074] S3. Use multithreading to perform element-wise multiplication of the input data matrix and the rotation factor matrix in parallel, obtaining intermediate calculation results. These intermediate results are used to prepare for the next step of matrix multiplication with the DFT matrix. This step specifically includes:

[0075] S31. The task of multiplying the input data matrix and the rotation factor matrix by matrix elements is evenly distributed to each thread, and then each thread reads the corresponding data from shared memory or global memory.

[0076] This example computes the Fourier transform of an input sequence of length 256. Typically, existing GPUs have 64 threads per block, so each thread is responsible for computing four data points, such as... Figure 2 As shown, this illustrates the data number held by each thread. Each thread needs to read its allocated data from shared memory or global memory into its own register. The reading process can employ merged memory access to speed up the reading process.

[0077] S32. Calculate the data corresponding to the rotation factor matrix based on the position of the input sequence data in the input data matrix.

[0078] Specifically, the elements of the input sequence data storage in X are calculated using GPU threads. inGiven the row and column numbers of the input data matrix, the corresponding elements in the rotation factor matrix are located based on these row and column numbers, and element-wise multiplication is performed. This invention uses multithreading to perform element-wise multiplication of the input data matrix and the rotation factor matrix in parallel. The rotation factor matrix is ​​calculated in real-time by a GPU thread according to the following rules:

[0079]

[0080] The element in row 15 and column 15 is Each element can be calculated using its row and column numbers in the input matrix.

[0081] S33. Use multithreading to perform element-wise multiplication of the input data matrix and the rotation factor matrix in parallel, and write the intermediate calculation results back to shared memory.

[0082] Specifically, all threads execute in parallel. After each thread completes its own computation task, they need to synchronize before proceeding to the next computation to ensure data correctness. The operations performed in this step are as follows:

[0083] TEMP 16×16 =T 16×16 ⊙X in

[0084] Among them, TEMP 16×16 This represents a temporary data matrix of size 16×16, used to temporarily store calculation results.

[0085] After step S2, X in The data is stored in the thread register. If elements are read into shared memory and then multiplied element by element, unnecessary memory accesses will be added, thus reducing computational efficiency. Therefore, it is necessary to explore X. in The data distribution pattern across threads is determined by directly retrieving the corresponding rotation factor matrix elements from the thread registers and performing multiplication operations.

[0086] like Figure 2 As shown, after experimental analysis, X in The data distribution across threads: In the MI250, one block is computed by 64 threads. Therefore, when computed on a block of length 256, each thread stores 4 elements in its local register. For example, thread 0's register stores X. in The data in thread X consists of elements 0, 1, 2, and 3, while thread 1 stores elements 16, 17, 18, and 19. Based on the distribution map, the location of the element stored in thread X can be calculated using the thread ID. in The row and column numbers in the matrix can be used to find the corresponding elements in the rotation factor matrix and perform matrix element multiplication.

[0087] S4. Call the matrix operation unit to perform matrix multiplication with the intermediate calculation result data obtained in step S3 and the DFT matrix to obtain the matrix operation result; assign c to c-1. If the value of c is not 0, transpose the matrix operation result and jump to step S2; if the value of c is 0, output the final calculation result data to obtain the n-dimensional Fourier transform sequence.

[0088] Specifically, step S4 can be divided into the following steps:

[0089] S41. Load the intermediate calculation results data and DFT matrix into the registers of each thread.

[0090] First, the intermediate calculation results are loaded into Fragment A, and the DFT matrix is ​​loaded into Fragment B. This step was already completed in step S2. Fragment C is set to zero. After the calculation is complete, the calculation result is retrieved from Fragment D. Here, Fragment is a data structure in GPU programming used to store the matrix data of the input matrix operation unit.

[0091] In the MI250 GPU, shared memory consists of 32 memory banks, which can be accessed concurrently by different threads. If different threads within the same thread bundle access the same memory bank, a memory bank conflict will occur, reducing the efficiency of shared memory access.

[0092] Since the input is a sequence of complex numbers, and the intermediate data generated during the calculation is also complex, a single-precision complex number requires 8 bytes of space. However, the width of the memory bank in the MI250 is 4 bytes. Therefore, a complex number needs to be stored in two contiguous memory banks. Accessing a complex number in shared memory requires accessing two contiguous memory banks, and other threads cannot access these two banks during the access process. Thus, with 32 memory banks, a maximum of 16 threads can access the number simultaneously, resulting in a memory bank conflict. If the real and imaginary parts of the complex number are stored in the same memory bank through data rearrangement, the 32 memory banks can support simultaneous access by 32 threads, doubling the memory access efficiency.

[0093] S42. Use the matrix operation units such as Tensor Core or Matrix Core in the GPU to perform matrix multiplication operations between the intermediate calculation result data and the DFT matrix, output the matrix operation result, and assign c to c-1.

[0094] This invention utilizes matrix operation units such as Tensor Cores or Matrix Cores in the GPU to perform matrix multiplication operations between intermediate computation results and the DFT matrix. Before performing matrix multiplication, the intermediate computation results and the DFT matrix need to be loaded into the registers of each thread. Then, all threads collaborate to complete the matrix operation and output the matrix operation result. That is, the following operations are performed:

[0095] X out =F 16 ·TEMP 16×16

[0096] like Figure 3 As shown, the matrix multiplication process is performed on the Matrix Core of the AMD GPU. The above steps merge 16 Fourier transform sequences of length 16 into a one-dimensional Fourier transform sequence of length 256. Each Fourier transform sequence of length 16 is also obtained recursively through the above calculation steps until the sequence length becomes 1. The Fourier transform of a sequence of length 1 is itself.

[0097] S43. When the value of c is not 0, transpose the matrix operation result and jump to step S2; when the value of c is 0, output the final calculation result to obtain the n-dimensional Fourier transform sequence.

[0098] Multidimensional FFT requires transposing the computational data before proceeding to steps S2, S3, and S4. Assuming an n-dimensional FFT is being calculated, this process needs to be repeated n-1 times. When c is 0, the computational result is output, yielding an n-dimensional Fourier transform sequence. The matrix transposition operation is performed by the GPU using multi-threaded collaboration. Matrix transposition is a crucial operation, significantly impacting the overall performance of multidimensional FFT.

[0099] Specifically, transposing the matrix operation result includes: transforming the position of the elements in the matrix and loading the data at the positions to be transformed into shared memory; using a merged memory access strategy to temporarily store the data in shared memory, with threads reading and writing in global memory in a cooperative manner; and filling the data blocks in shared memory so that the data accessed by different threads is distributed in different storage bodies.

[0100] In this embodiment, an efficient matrix transpose scheme is adopted, which makes the transpose of a three-dimensional matrix equivalent to data copying. Transpose involves changing the position of elements in the matrix, loading the data at the positions to be transformed into shared memory, and then each thread retrieves data from shared memory for exchange. The key point is to optimize the access efficiency of global memory and shared memory. Therefore, this invention mainly adopts the following two optimization strategies to transpose the matrix operation result, specifically including:

[0101] First, by temporarily storing data in shared memory, threads can read and write to global memory in a cooperative manner, i.e., merging memory accesses, thus improving memory access efficiency. Second, by filling data blocks in shared memory, data accessed by different threads is distributed across different storage banks, thereby avoiding storage conflicts.

[0102] like Figure 4 As shown, Figure 4 This is a schematic diagram of the calculation process of the three-dimensional Fourier transform in this embodiment of the invention. When the dimension n=3, calculating the 3D FFT requires performing a one-dimensional FFT transform in each of the three directions. Assuming the three directions of the 3D matrix are XYZ, the initial storage order is ZYX, meaning the data is stored continuously in the X direction. After completing the one-dimensional Fourier transform in the X direction, the storage order of the data needs to be transposed to ZXY, meaning the data is continuous in the Y direction, and then a one-dimensional Fourier transform in the Y direction is performed. Otherwise, a large amount of discontinuous storage space needs to be accessed, resulting in extremely low memory access efficiency and reduced computational efficiency. Since calculating the 3D FFT here requires transposing twice, steps S2, S3, and S4 above are repeated twice. Finally, a 256×256×256 three-dimensional Fourier transform sequence is obtained.

[0103] S5. The host side copies the final calculation result data from the device's global memory to the host memory, releasing computing resources.

[0104] Finally, the FFT calculation result is stored in the device's global memory. The host needs to copy the calculation result from the global memory to the host memory for output or other uses. Since a lot of device global memory and host memory space is allocated before calculating the FFT, it needs to be released after the calculation to prevent memory leaks.

[0105] like Figure 5 The figure shows the speedup effect of the multidimensional fast Fourier transform acceleration method compared to other FFT libraries. On the AMD MI250 GPU platform, the performance comparison results are shown between the proposed method, the AMD-developed rocFFT library, and the open-source cross-platform vkFFT library. The horizontal axis represents the size of the 3D FFT, and the vertical axis represents the computational efficiency of the FFT, measured in TFLOPs (trillion floating-point operations per second). The data in the figure shows that the proposed method has higher computational efficiency, averaging 1.5 times that of rocFFT and 2.3 times that of vkFFT.

[0106] In summary, this invention first requests computing resources and generates the necessary preliminary data for computation. Then, it starts the computing kernel, uses multithreading to perform matrix element-wise multiplication in parallel, calls high-speed matrix operation units to perform classical matrix multiplication, copies the result data to host memory, completes the computation, and releases the computing resources. Because this invention uses high-speed matrix operation units on existing GPUs to perform FFT calculations, it significantly improves the computational efficiency of FFT, offering a clear advantage over traditional FFT calculation methods.

[0107] The above embodiments are preferred embodiments of the present invention, but the embodiments of the present invention are not limited to the above embodiments. Any changes, modifications, substitutions, combinations, or simplifications made without departing from the spirit and principle of the present invention shall be considered equivalent substitutions and shall be included within the protection scope of the present invention.

Claims

1. A method for accelerating multidimensional fast Fourier transform based on matrix operations, characterized in that, Includes the following steps: S1. Request computing resources to compute an n-dimensional FFT. Copy the computation data from the host to the device's global memory, generate the DFT matrix and the input data matrix, and initialize the variable c, setting c=n. S2. Start the computing kernel and read the DFT matrix data and input data matrix from global memory into registers and shared memory; S3. Use multithreading to perform the matrix element multiplication task between the input data matrix and the rotation factor matrix in parallel, and obtain intermediate calculation result data. The task of multiplying the input data matrix with the rotation factor matrix by matrix elements is evenly distributed to each thread, and then each thread reads the corresponding data from shared memory or global memory. The data in the rotation factor matrix is ​​calculated based on the position of the input sequence data in the input data matrix. The row and column numbers of the elements storing the input sequence data in the input data matrix are calculated by the GPU thread, and the corresponding elements in the rotation factor matrix are found based on the row and column numbers. The element-wise multiplication of the input data matrix and the rotation factor matrix is ​​performed in parallel using multiple threads, and the intermediate calculation results are written back to shared memory. S4. Call the matrix operation unit to perform matrix multiplication between the intermediate calculation result data obtained in step S3 and the DFT matrix to obtain the matrix operation result; assign c to c-1; if the value of c is not 0, transpose the matrix operation result and jump to step S2; if the value of c is 0, output the final calculation result data. S5. The host copies the final calculation result data from the device's global memory to the host memory, releasing computing resources.

2. The method for accelerating multidimensional fast Fourier transform based on matrix operations according to claim 1, characterized in that, Step S1 includes: The input sequence is analyzed on the host side. Based on the length of the input sequence, the CPU on the host side generates the corresponding DFT matrix and copies the corresponding DFT matrix to the device's global memory. Arrange the input sequence into an input data matrix in row-major order.

3. The method for accelerating multidimensional fast Fourier transform based on matrix operations according to claim 1, characterized in that, Step S2 includes: if the input data matrix cannot be placed into the shared memory at once, the input data matrix is ​​further decomposed and rearranged so that the row length of the data matrix is ​​less than or equal to the size of the shared memory.

4. The method for accelerating multidimensional fast Fourier transform based on matrix operations according to claim 1, characterized in that, Step S4 includes: The intermediate calculation results and DFT matrix are loaded into the registers of each thread; The matrix operation unit in the GPU is used to perform matrix multiplication between the intermediate calculation result data and the DFT matrix, and the matrix operation result is output. The value of c is assigned to c-1. If the value of c is not 0, transpose the matrix operation result and jump to step S2; if the value of c is 0, output the final calculation result.

5. The method for accelerating multidimensional fast Fourier transform based on matrix operations according to claim 4, characterized in that, The transpose of the matrix operation result includes: transforming the position of the elements in the matrix and loading the data at the positions to be transformed into shared memory; using a merged memory access strategy to temporarily store the data in shared memory, with threads reading and writing in global memory in a cooperative manner; and filling the data blocks in shared memory so that the data accessed by different threads is distributed in different storage bodies.

6. The method for accelerating multidimensional fast Fourier transform based on matrix operations according to claim 4, characterized in that, The matrix operation units in the GPU include: Tensor Cores or Matrix Cores in the GPU.

Images