A sparse LU decomposition solution method and related apparatus based on matrix computation unit optimization

Abstract

This application discloses a sparse LU decomposition solution method and related apparatus based on matrix computation unit optimization, belonging to the field of numerical algorithm optimization technology. The method includes: reordering the sparse matrix using an improved PRCM algorithm to increase data local density, obtaining a bandwidth-compressed sparse matrix; then converting it into a regular, dense tile block structure and storing it in pre-allocated contiguous memory space; and finally, employing a right-look-based LU decomposition algorithm to numerically decompose the bandwidth-compressed sparse matrix. During the numerical decomposition process, a matrix computation unit is used for acceleration to obtain the LU matrix, which can then be used to efficiently solve the linear equation. This application's solution fully utilizes the spatial locality of the sparse matrix while leveraging the powerful computational capabilities of the matrix computation unit, reducing computational load while improving computational efficiency, significantly reducing the time spent on symbolic and numerical decomposition stages of sparse LU decomposition.

Description

Technical Field

[0001] This application relates to the field of numerical algorithm optimization technology, and in particular to a sparse LU decomposition solution method and related apparatus based on matrix computation unit optimization. Background Technology

[0002] Sparse LU decomposition is an important algorithm in numerical linear algebra, mainly used to solve linear systems of sparse matrices Ax = b. A is a sparse matrix in which most elements are zero. LU decomposition splits the sparse matrix into the product of a lower triangular matrix (L) and an upper triangular matrix (U), such that A = LU, thus simplifying the solution process. Compared with directly solving the matrix equation, the advantage of LU decomposition is that it reduces redundant calculations through matrix decomposition, making the subsequent solution process more efficient. Sparse LU decomposition mainly consists of three stages: (1) reordering, (2) symbolic decomposition, and (3) numerical decomposition. The purpose of reordering is to reduce the number of newly added non-zero elements and maintain numerical stability; the purpose of symbolic decomposition is to determine the structure of the L and U matrices; numerical decomposition is a floating-point calculation to find L and U.

[0003] However, the inherent complex data dependencies and irregular distribution of sparse matrices pose significant challenges to parallelization algorithms. Complex data dependencies mean it's difficult to find parallelizable components, and the parallelization granularity is typically small, resulting in minimal speedup. Furthermore, sparse matrices often use storage formats like CSC and CSR, whose computation itself requires irregular memory accesses, leading to additional overhead. This irregularity of sparse matrices limits parallel optimization.

[0004] In summary, existing sparse matrix LU decomposition techniques have the following three main problems: (1) The reordering method aims to reduce the number of newly added non-zero elements, but its side effect is that the data distribution becomes more scattered and irregular, resulting in low effective computational efficiency of the matrix computation unit. (2) The data structure is non-uniformly partitioned, which is not suitable for the matrix computation unit structure and conflicts with the matrix computation unit's requirement for uniform and dense input, making it difficult to leverage the powerful computational performance advantages of heterogeneous computing hardware. (3) During the numerical decomposition calculation process, additional format conversion and data movement are required, resulting in high indirect memory access overhead.

[0005] Therefore, developing a sparse LU decomposition solution scheme that can utilize matrix computation units for optimization is crucial for accelerating the solution of large-scale sparse linear equation systems in practical applications and improving their performance in real-world applications. Summary of the Invention

[0006] The purpose of this application is to provide a sparse LU decomposition solution method and related apparatus based on matrix computation unit optimization, which can improve computational efficiency while reducing computational load, and significantly improve the performance of sparse LU decomposition in practical applications.

[0007] To achieve the above objectives, this application provides the following solution:

[0008] In a first aspect, this application provides a sparse LU decomposition solution method based on matrix computation unit optimization, comprising the following steps:

[0009] Obtain the sparse matrix to be decomposed; the sparse matrix to be decomposed is the sparse matrix in the linear equation to be solved.

[0010] The MC64 algorithm is called to reorder the sparse matrix to be decomposed, thereby improving data stability and obtaining a sparse matrix after one reordering.

[0011] The PRCM algorithm is called to reorder the sparse matrix after the first reordering, thereby improving the local density of the data and obtaining a sparse matrix with bandwidth compression.

[0012] The bandwidth-compressed sparse matrix is ​​converted into a regular dense tile structure and stored in a pre-allocated contiguous memory space.

[0013] A right-look-based LU decomposition algorithm is used to numerically decompose the bandwidth-compressed sparse matrix to obtain the corresponding LU matrix; matrix computation units are used to accelerate the numerical decomposition of the bandwidth-compressed sparse matrix.

[0014] By using LU matrices to solve the lower and upper triangular equation systems, an efficient solution to the linear equations to be solved can be achieved.

[0015] Optionally, the PRCM algorithm is called to reorder the sparse matrix after the first reordering to improve the local density of the data and obtain a bandwidth-compressed sparse matrix. The specific steps include:

[0016] With the goal of minimizing the bandwidth of the sparse matrix, key nodes are identified from the sparse matrix after one reordering; key nodes are those that affect the bandwidth compression operation.

[0017] Move several key nodes to the bottom right corner of the sparse matrix after one reordering.

[0018] The RCM algorithm is used to perform bandwidth compression on the subgraph consisting of several nodes remaining after removing several key nodes, resulting in a bandwidth-compressed sparse matrix.

[0019] Optionally, the degree of the key node is greater than the node degree threshold; the node degree threshold for minimizing the bandwidth of the sparse matrix is ​​determined according to the following formula:

[0020]

[0021] Where t is the node degree threshold, and T is the set of degrees of nodes in the sparse matrix. Let Γ be the minimum estimate of the bandwidth of the sparse matrix, where a is a constant and Γ = {2}. 0 ,2 1 ,…,2 14}, where S(t) is the number of nodes with a degree greater than t.

[0022] Optionally, a right-look-based LU decomposition algorithm is used to numerically decompose the bandwidth-compressed sparse matrix to obtain the corresponding LU matrix, specifically including the following steps:

[0023] Within the bandwidth-compressed sparse matrix, a sliding window is used to select several tile structures. The size of the sliding window depends on the number of tile structures that need to be computed in the current round.

[0024] For the tile-based structure within the sliding window, perform LU decomposition for the current round to obtain the LU matrix vector for the current round. The LU matrix vector for the current round includes the row vectors of the U matrix and the column vectors of the L matrix for the current round.

[0025] The sliding window moves down layer by layer along the bandwidth-compressed sparse matrix, selects several tile block structures, and jumps to "for several tile block structures within the sliding window, perform LU decomposition for the current round to obtain the LU matrix vector for the current round", until the LU decomposition of the bandwidth-compressed sparse matrix is ​​completed.

[0026] Based on the LU matrix vectors obtained from each iteration, the corresponding LU matrix is ​​obtained.

[0027] Optionally, for several tile blocks within the sliding window, perform LU decomposition for the current round to obtain the LU matrix vector for the current round. This specifically includes the following steps:

[0028] Based on the diagonal blocks of the previous round within the sliding window, the LU matrix vector of the current round is obtained by decomposition.

[0029] Based on the LU matrix vector of the current round, the diagonal block of the next round is calculated, and the Schur complement is calculated for other tile block structures within the sliding window. When performing the Schur complement calculation, the tile block structure is used as the unit, and the matrix calculation unit is used to accelerate the 3-level BLAS operation.

[0030] Optionally, while the computation stream performs the LU decomposition for the current round, a copy stream is used to store the decomposition results completed in the previous round and to prefetch the data for the next round.

[0031] Secondly, this application provides a sparse LU decomposition solver based on matrix computation unit optimization, comprising the following modules:

[0032] The sparse matrix acquisition module is used to obtain the sparse matrix to be decomposed; the sparse matrix to be decomposed is the sparse matrix in the linear equation to be solved.

[0033] The sparse matrix reordering module is used to call the MC64 algorithm to reorder the sparse matrix to be decomposed, improve data stability, and obtain the sparse matrix after one reordering.

[0034] The sparse matrix bandwidth compression module is used to call the PRCM algorithm to reorder the sparse matrix after one reordering, improve the local density of the data, and obtain the bandwidth-compressed sparse matrix.

[0035] The sparse matrix block storage module is used to convert the bandwidth-compressed sparse matrix into a regular dense Tile block structure and store it in a pre-allocated contiguous memory space.

[0036] The sparse matrix numerical decomposition module is used to numerically decompose the bandwidth-compressed sparse matrix using a right-look-based LU decomposition algorithm to obtain the corresponding LU matrix; the matrix calculation unit is used to accelerate the numerical decomposition of the bandwidth-compressed sparse matrix.

[0037] The linear equation solving module is used to solve the lower triangular equation system and the upper triangular equation system based on the LU matrix, so as to complete the efficient solution of the linear equation to be solved.

[0038] Thirdly, this application provides a computer device, including: a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement the steps of the sparse LU decomposition solution method based on matrix computation unit optimization described above.

[0039] Fourthly, this application provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the steps of the sparse LU decomposition solution method based on matrix computation unit optimization described above.

[0040] Fifthly, this application provides a computer program product, including a computer program that, when executed by a processor, implements the steps of the sparse LU decomposition solution method based on matrix computation unit optimization described above.

[0041] According to the specific embodiments provided in this application, the following technical effects are disclosed:

[0042] This application provides a sparse LU decomposition solution method and related apparatus based on matrix computation unit optimization. In this method, the MC64 algorithm is first used to reorder the sparse matrix to be decomposed, improving data stability and obtaining a reordered sparse matrix. Then, the PRCM algorithm is used to reorder the reordered sparse matrix again, improving local data density and obtaining a bandwidth-compressed sparse matrix. Subsequently, the bandwidth-compressed sparse matrix is ​​converted into a regular dense tile block structure and stored in a pre-allocated contiguous memory space. A right-look-based LU decomposition algorithm is used to numerically decompose the bandwidth-compressed sparse matrix. During the numerical decomposition process, a matrix computation unit is used for acceleration to obtain the corresponding LU matrix. Finally, the lower triangular equation system and the upper triangular equation system can be solved based on the LU matrix to complete the efficient solution of the linear equation to be solved. This application uses an improved PRCM algorithm for reordering, which can solve some of the failures of the RCM algorithm and further enhance the local density and compression bandwidth of the data, reduce the amount of computation and improve the computational efficiency. The data structure with regular dense blocks used after reordering can directly call the matrix computation unit for acceleration in the subsequent numerical decomposition stage. It also saves the indirect memory access overhead caused by format conversion in traditional methods, makes full use of the spatial locality of sparse matrices and leverages the powerful computing power of the matrix computation unit, significantly reducing the time of sparse LU decomposition in the symbolic decomposition and numerical decomposition stages. Attached Figure Description

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

[0044] Figure 1 This is a flowchart illustrating a sparse LU decomposition solution method based on matrix computation unit optimization, provided as an embodiment of this application.

[0045] Figure 2 This is a schematic diagram comparing the RCM reordering results and the PRCM reordering results in a sparse LU decomposition solution method based on matrix computation unit optimization provided in an embodiment of this application.

[0046] Figure 3 This is a flowchart of step A3 in a sparse LU decomposition solution method based on matrix computation unit optimization provided in an embodiment of this application.

[0047] Figure 4 This is a schematic diagram of the symbolic decomposition of the GaAsH6 matrix in a sparse LU decomposition solution method based on matrix computation unit optimization provided in an embodiment of this application.

[0048] Figure 4 (a) is a schematic diagram of the original GaAsH6 matrix in a sparse LU decomposition solution method based on matrix computation unit optimization provided in an embodiment of this application.

[0049] Figure 4 (b) is a schematic diagram of the GaAsH6 matrix after reordering in a sparse LU decomposition solution method based on matrix computation unit optimization provided in an embodiment of this application.

[0050] Figure 4 (c) is a schematic diagram of the symbolic decomposition of the GaAsH6 matrix in a sparse LU decomposition solution method based on matrix computation unit optimization provided in an embodiment of this application.

[0051] Figure 5 This is a flowchart of step A5 in a sparse LU decomposition solution method based on matrix computation unit optimization provided in an embodiment of this application.

[0052] Figure 6 This is a schematic diagram illustrating the execution process of numerical decomposition in a sparse LU decomposition solution method based on matrix computation unit optimization provided in an embodiment of this application.

[0053] Figure 7 This is a schematic diagram of the execution timeline of a sparse LU decomposition solution method based on matrix computation unit optimization provided in an embodiment of this application.

[0054] Figure 8 This diagram illustrates the data dependencies during numerical decomposition in a sparse LU decomposition solution method based on matrix computation unit optimization, as provided in an embodiment of this application.

[0055] Figure 9 This diagram illustrates a comparison between a sparse LU decomposition solution method based on matrix computation unit optimization provided in an embodiment of this application and the test results of traditional SuperLU_DIST and PanguLU.

[0056] Figure 10 This is a schematic diagram of the functional modules of a sparse LU decomposition solver based on matrix computation unit optimization, provided in an embodiment of this application.

[0057] Figure 11 This is a schematic diagram of the structure of a computer device provided in an embodiment of this application. Detailed Implementation

[0058] The technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, and not all embodiments. Based on the embodiments of this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0059] To make the above-mentioned objectives, features and advantages of this application more apparent and understandable, the application will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0060] This application provides a sparse LU decomposition solution method based on matrix computation unit optimization, such as... Figure 1 As shown, it includes the following steps:

[0061] A1. Obtain the sparse matrix to be decomposed; the sparse matrix to be decomposed is the sparse matrix in the linear equation to be solved.

[0062] A2. The MC64 algorithm is called to reorder the sparse matrix to be decomposed, improving data stability and obtaining a reordered sparse matrix. In this embodiment, the MC64 algorithm is used for static pivot selection before bandwidth compression to improve numerical stability.

[0063] A3. Apply the PRCM algorithm to reorder the sparse matrix after the first reordering to improve the local density of the data, resulting in a bandwidth-compressed sparse matrix. The improved PRCM algorithm is used for this second reordering to aggregate the non-zero elements in the sparse matrix while reducing the number of newly added non-zero elements.

[0064] The PRCM (Partial Reverse Cuthill-Mckee) algorithm is based on RCM and further improves data locality. The idea behind PRCM is to move some nodes with high degrees to the bottom right corner of the matrix, and then apply the RCM algorithm to the remaining nodes. The function of RCM is to cluster the non-zero elements of the matrix near the diagonal, thus compressing the matrix bandwidth. However, RCM fails in some cases, such as matrices in circuit simulation, because these matrices often contain nodes with very high degrees. These nodes amplify the bandwidth of the entire sparse matrix, preventing non-zero elements from clustering near the diagonal.

[0065] The PRCM algorithm first identifies nodes with extremely high degree that significantly impact compression; these nodes are called "critical nodes." These nodes are then moved to the bottom right corner of the matrix, and RCM is applied to the remaining subgraph for bandwidth compression. Placing critical nodes in the bottom right corner reduces the addition of new non-zero elements, preventing the creation of new non-zero elements outside the diagonal and bottom right corners during Shure complement calculation. Without the influence of critical nodes, the remaining subgraph can be effectively compressed using the RCM algorithm. Figure 2 As shown, the left side represents the result of RCM reordering, and the right side represents the result of PRCM reordering. The PRCM algorithm moves the key nodes to the lower right corner, forming a downward-pointing triangular arrow structure. The bandwidth after the original RCM reordering is B = β, and the equivalent bandwidth after PRCM reordering is B' = β' + S, where S is the number of nodes with a degree greater than the node degree threshold, and β' is the bandwidth of the subgraph formed by the remaining nodes after RCM reordering. To achieve better bandwidth compression in PRCM reordering, B' needs to be... <B。

[0066] In this embodiment, as Figure 3 As shown, step A3 specifically includes the following steps:

[0067] A31. With the goal of minimizing the bandwidth of the sparse matrix, identify the key nodes in the sparse matrix after one reordering; the key nodes are those that affect the bandwidth compression operation.

[0068] A32. Move several key nodes to the bottom right corner of the sparse matrix after one reordering.

[0069] A33. Use the RCM algorithm to perform bandwidth compression on the subgraph composed of several nodes remaining after removing several key nodes, and obtain a sparse matrix after bandwidth compression.

[0070] In step A31, to find these critical nodes, a node degree threshold needs to be set; nodes with a degree greater than this threshold are considered critical nodes. How to select this threshold is a problem that needs to be solved, and its purpose is to minimize bandwidth B'.

[0071] make T is the set of node degrees, S is the number of nodes greater than the threshold t, and β' is the bandwidth of the remaining subgraph after RCM reordering. We need to select a suitable threshold from these possible thresholds to minimize the bandwidth B. A common approach is to iterate through the subgraph, calculating the threshold for each t in T, and then selecting the t that minimizes B. This is the threshold that minimizes the bandwidth. Finally, we move all nodes with degrees greater than t to the bottom right corner of the matrix and reorder the remaining nodes using RCM. However, calculating the RCM algorithm for each t is too time-consuming, which is unacceptable. Therefore, this embodiment proposes an estimation method, letting β'(t) = at, where a is a constant, and taking a set of possible values ​​for a Γ = {2 0 ,2 1 ,…,2 14 The bandwidth obtained using the estimated value at + S(t) is called the estimated bandwidth B. * This yields a set of functions B for estimating bandwidth. * =at + S(t), a∈Γ. Next, we need to find B from this set. * Find the minimum value in the function The minimum estimated bandwidth is used to approximate the minimum actual bandwidth, and the corresponding t is the node degree threshold t to be found.

[0072] Specifically, in this embodiment, the degree of the key node is greater than the node degree threshold; the node degree threshold when the sparse matrix bandwidth is minimized can be determined according to the following formula:

[0073]

[0074] Where t is the node degree threshold, and T is the set of degrees of nodes in the sparse matrix. Let Γ be the minimum estimate of the bandwidth of the sparse matrix, where a is a constant and Γ = {2}. 0 ,2 1 ,…,2 14 Let S(t) be the number of nodes with a degree greater than t. Actual testing shows that the bandwidth estimate differs from the actual value by no more than 1%, and the time complexity of the bandwidth estimate calculation method is far superior to the actual bandwidth calculation method.

[0075] By optimizing the reordering algorithm, the spatial structure of the matrix is ​​improved, causing non-zero elements to cluster more closely around the diagonal. This effectively controls the number of newly added non-zero elements during Shure complement calculation. Compressing the matrix bandwidth reduces computational load and saves GPU memory usage, thereby accelerating the computational efficiency of numerical decomposition. For matrices that RCM cannot handle well in fields such as circuit simulation, PRCM can also compress their bandwidth and improve them into a computationally efficient structure.

[0076] A4. Convert the bandwidth-compressed sparse matrix into a regular dense Tile block structure and store it in a pre-allocated contiguous memory space.

[0077] After reordering, the next step is to convert the sparse matrix into a block format and determine the structure of the LU. First, the CSC format of the sparse matrix is ​​converted into a regular, dense block structure for storage; these blocks are called Tiles. Traditional solvers determine the LU structure before block decomposition, while this embodiment decomposes the blocks first and then determines the structure, only needing to determine the positions of the newly added blocks. Compared to element-wise symbolic decomposition, block-wise symbolic decomposition significantly reduces time overhead and eliminates the need for complex structures and algorithms like tree elimination found in traditional methods.

[0078] like Figure 4 As shown, this example is a symbolic decomposition of the GaAsH6 matrix. Figure 4 (a) is the original matrix. Figure 4 (b) is the reordered matrix, after being divided into blocks. Figure 4 (c) is the final structure determined after symbolic decomposition. Since sparse matrices consume a large amount of memory, contiguous storage space is allocated in memory to store them, changing the scattered data into contiguous data storage, thus saving the time of CPU-GPU memory copying.

[0079] A5. A right-look-based LU decomposition algorithm is used to numerically decompose the bandwidth-compressed sparse matrix to obtain the corresponding LU matrix; matrix computation units are used to accelerate the numerical decomposition of the bandwidth-compressed sparse matrix.

[0080] The numerical decomposition process in the sparse LU decomposition algorithm is the most time-consuming part. This step directly builds upon the block structure from the previous step and utilizes matrix computation units for acceleration. This embodiment uses uniformly partitioned dense block storage, which avoids unnecessary format conversion overhead and allows for the invocation of a 3-level BLAS and the use of matrix computation units to accelerate matrix operations.

[0081] In this embodiment, the numerical decomposition algorithm used is LU decomposition based on right-looking, which optimizes the dependency relationship between the diagonal block and the L column block and the U row block in each iteration, simplifies the two-layer dependency relationship into one layer, and uses the sliding window strategy to mask the latency of computation and communication, enhances data reusability, and improves data-level parallelism and computation-to-memory ratio.

[0082] In this embodiment, as Figure 5 As shown, step A5 specifically includes the following steps:

[0083] A51. In a sparse matrix after bandwidth compression, a sliding window is used to select several tile blocks. The size of the sliding window depends on the number of tile blocks to be computed in the current round. In each round, all the tiles to be computed constitute the sliding window. The sliding window consists of three parts: the unused portion of the result from the previous round, which is returned to the CPU; the tiles to be computed in this round (including some tiles used in the previous round); and the tiles to be used in the next round, with data prefetched. These three parts together constitute a sliding window. During numerical decomposition, the sliding window method is used to complete data backhaul and prefetching simultaneously with floating-point calculations, creating a mutual masking of data transmission and computation. Reusable portions are still retained within the window, improving data reusability.

[0084] For example, in the k-th iteration, the equivalent bandwidth of the matrix is ​​B. k If the dimension of the tile structure is ts (the size of a tile structure is ts × ts), then the dimension of this sliding window is approximately The data type is an 8-byte double-precision floating-point number.

[0085] A52. For the tile block structure within the sliding window, perform LU decomposition for the current round to obtain the LU matrix vector for the current round. Step A52 specifically includes the following steps:

[0086] Based on the diagonal blocks of the previous round within the sliding window, the LU matrix vector of the current round is obtained. The LU matrix vector of the current round includes the row vectors of the U matrix and the column vectors of the L matrix in the current round.

[0087] Based on the LU matrix vector of the current round, the diagonal block of the next round is calculated, and the Schur complement is calculated for other tile block structures within the sliding window. When performing the Schur complement calculation, the tile block structure is used as the unit, and the matrix calculation unit is used to accelerate the 3-level BLAS operation.

[0088] A53. Move the sliding window down layer by layer along the bandwidth-compressed sparse matrix, select several tile block structures, and jump to step A52 until the LU decomposition of the bandwidth-compressed sparse matrix is ​​completed.

[0089] A54. Based on the LU matrix vectors obtained from each iteration, the corresponding LU matrix is ​​obtained.

[0090] In this embodiment, the numerical decomposition process described above uses independent computation streams and copy streams to decouple computation and data movement, reducing the dependence of computation on data and improving the efficiency of data transmission and computation. That is, while the computation stream is performing the LU decomposition of the current round, the copy stream is used to store the decomposition results completed in the previous round and to prefetch the data for the next round.

[0091] Specifically, the specific execution process is as follows: Figure 6 As shown, for the k-th iteration, in the LU decomposition (light gray part, based on the diagonal block D of the (k-1)-th iteration), k-1 The LU matrix vectors L and U for the current round are obtained by decomposition, and the diagonal block D for round k is calculated. k While performing Shur complement calculations (S) on other tile structures within the sliding window, a copy stream is used to copy the decomposition results (dark gray portion) from the GPU back to the CPU. This portion is no longer used in subsequent calculations. Simultaneously, data needed in round k+1 is copied from the CPU to the GPU, completing data prefetching. Data is transferred bidirectionally. Figure 6 In Chinese, D2H means data transfer from GPU to CPU, and H2D means data transfer from CPU to GPU.

[0092] The timeline of this execution process is as follows: Figure 7 As shown, the copy stream is... Figure 7 The copy operation in LU decomposition. The computational flow is divided into three parts: diagonal block decomposition (D), decomposition of L columns and U rows (L, U), and Schur complement computation (S); the first computational flow solves for one column of L and one row of U, i.e. Figure 7 The first computational flow involves calculating the panel; the second computational flow decomposes the diagonal blocks for the next round of computation, i.e. Figure 7 Diag in; the third computational flow, Schur complement computation, i.e. Figure 7 In Schur computation, multiple steps executed in parallel are aligned within a time window through synchronization operations.

[0093] Dependencies and parallelism levels, such as Figure 8 As shown, first for D k-1 Perform GETRF (Diagonal Block Decomposition), then simultaneously decompose column L and row U, and finally simultaneously decompose column D. k Perform GEMMTRF (Shure complement calculation + diagonal block decomposition) and S-Shure complement calculation.

[0094] In each iteration, the sliding window, consisting of copying and computation, moves once. Most blocks within the window are reused, resulting in high data reuse and significantly reducing cache misses. Shure complement computation is performed in blocks, using matrix computation units to accelerate the 3-level BLAS operation, improving the granularity of data parallelism and enhancing parallelism. Because PRCM reordering clusters non-zero elements near the diagonal, the sliding window size required for each iteration is relatively small, saving GPU memory and allowing for simultaneous computation on the GPU.

[0095] The PCIe bus is full-duplex, allowing D2H and H2D data transfers to occur simultaneously with roughly the same data volume. According to the definition of the computation-copy ratio (CCR), the CCR of a sliding window depends only on the bandwidth at the current position. When the CCR of the sliding window is greater than the CCR of the GPU, the computation can be completely overridden by communication. For the NVIDIA A800 GPU, its... Theoretically, when the bandwidth B > 624 × 8 = 4992, the communication time can cover the computation time. Therefore, it performs well in large-scale sparse matrices.

[0096] Before each round of numerical decomposition, the L columns and U rows of the result from the previous round are transferred from the GPU back to the CPU. Simultaneously, the tiles needed for the next round are prefetched and transferred from the CPU to the GPU. During data transfer, L column and U row decomposition and Shur complement updates are performed. After the numerical decomposition is complete, the original sparse matrix A is transformed into the product of the L and U matrices, i.e., A = LU.

[0097] A6. Solve the lower triangular equation system Ly=b and the upper triangular equation system Ux=y based on the LU matrix to achieve an efficient solution for the linear equation to be solved.

[0098] To evaluate the advantages of the proposed sparse LU decomposition solution scheme based on matrix computation unit optimization, 165 large-scale sparse matrices with dimensions ranging from 50K to 2.7M were selected from the SuiteSparse Matrix Collection dataset. The two leading solvers in the industry, SuperLU_DIST and PanguLU, were compared. In the most time-consuming numerical decomposition stage, compared to SuperLU_DIST, the proposed scheme achieved a speedup of 17.9 times and a maximum speedup of 86.2 times; compared to PanguLU, the proposed scheme achieved an average speedup of 8.8 times and a maximum speedup of 90.9 times. Test results are as follows... Figure 9 As shown.

[0099] The sparse LU decomposition solution method based on matrix computation unit optimization provided in the above embodiments of this application can solve some cases where the RCM algorithm fails by improving the RCM algorithm, and further enhances the local density of data and compression bandwidth, reducing the amount of computation and improving computational efficiency. The reordering algorithm used in this application is different from traditional solvers. It no longer aims to minimize newly added non-zero elements, but enhances the locality of data. This makes data reuse more efficient in the sliding window, and the compression of matrix bandwidth results in lower GPU memory usage. At the same time, due to the use of a uniformly dense block data structure, the 3-level BLAS interface accelerated by the matrix computation unit can be directly called in the numerical decomposition stage; using uniform blocks and operating directly on these blocks eliminates the indirect memory access overhead caused by format conversion in traditional methods. This application is implemented on the NVIDIA GPU platform, making full use of Tensor Core components, and achieving an average peak performance of 22.5%, while traditional solvers only achieve less than 1% of the peak performance.

[0100] Based on the same inventive concept, this application also provides an apparatus for implementing the sparse LU decomposition solution method based on matrix computation unit optimization described above. The solution provided by this apparatus is similar to the implementation scheme described in the above method. Therefore, the specific limitations of one or more embodiments of the sparse LU decomposition solution apparatus based on matrix computation unit optimization provided below can be found in the limitations of the sparse LU decomposition solution method based on matrix computation unit optimization described above, and will not be repeated here.

[0101] In one exemplary embodiment, such as Figure 10 As shown, a sparse LU decomposition solver based on matrix computation unit optimization is provided, comprising the following modules:

[0102] The sparse matrix acquisition module is used to obtain the sparse matrix to be decomposed; the sparse matrix to be decomposed is the sparse matrix in the linear equation to be solved.

[0103] The sparse matrix reordering module is used to call the MC64 algorithm to reorder the sparse matrix to be decomposed, improve data stability, and obtain the sparse matrix after one reordering.

[0104] The sparse matrix bandwidth compression module is used to call the PRCM algorithm to reorder the sparse matrix after one reordering, improve the local density of the data, and obtain the bandwidth-compressed sparse matrix.

[0105] The sparse matrix block storage module is used to convert the bandwidth-compressed sparse matrix into a regular dense Tile block structure and store it in a pre-allocated contiguous memory space.

[0106] The sparse matrix numerical decomposition module is used to numerically decompose the bandwidth-compressed sparse matrix using a right-look-based LU decomposition algorithm to obtain the corresponding LU matrix; the matrix calculation unit is used to accelerate the numerical decomposition of the bandwidth-compressed sparse matrix.

[0107] The linear equation solving module is used to solve the lower triangular equation system and the upper triangular equation system based on the LU matrix, so as to complete the efficient solution of the linear equation to be solved.

[0108] certainly, Figure 10 The architecture shown is merely exemplary; it can be omitted as needed when implementing different functionalities. Figure 10 One or at least two components of the system shown.

[0109] In one exemplary embodiment, a computer device is provided, which may be a server or a terminal, and its internal structure diagram may be as follows. Figure 11 As shown, the computer device includes a processor, memory, input / output (I / O) interfaces, and a communication interface. The processor, memory, and I / O interfaces are connected via a system bus, and the communication interface is also connected to the system bus via the I / O interfaces. The processor provides computational and control capabilities. The memory includes a non-volatile storage medium and internal memory. The non-volatile storage medium stores the operating system, computer programs, and a database. The internal memory provides an environment for the operation of the operating system and computer programs stored in the non-volatile storage medium. The I / O interfaces are used for exchanging information between the processor and external devices. The communication interface is used for communicating with external terminals via a network connection. When the computer program is executed by the processor, it can implement the steps of the sparse LU decomposition solution method based on matrix computation unit optimization provided in the above embodiment.

[0110] Those skilled in the art will understand that Figure 11 The structure shown is merely a block diagram of a portion of the structure related to the present application and does not constitute a limitation on the computer device to which the present application is applied. Specific computer devices may include more or fewer components than those shown in the figure, or combine certain components, or have different component arrangements.

[0111] In one exemplary embodiment, a computer device is also provided, including a memory and a processor, wherein the memory stores a computer program, and the processor executes the computer program to implement the steps in the above-described method embodiments.

[0112] In one exemplary embodiment, a computer-readable storage medium is provided storing a computer program that, when executed by a processor, implements the steps in the above-described method embodiments.

[0113] In one exemplary embodiment, a computer program product is provided, including a computer program that, when executed by a processor, implements the steps in the above-described method embodiments.

[0114] It should be noted that the user information (including but not limited to user device information, user personal information, etc.) and data (including but not limited to data used for analysis, data stored, data displayed, etc.) involved in this application are all information and data authorized by the user or fully authorized by all parties, and the collection, use and processing of the relevant data must comply with relevant regulations.

[0115] Those skilled in the art will understand that all or part of the processes in the methods of the above embodiments can be implemented by a computer program instructing related hardware. The computer program can be stored in a non-volatile computer-readable storage medium, and when executed, it can include the processes of the embodiments of the above methods. Any references to memory, databases, or other media used in the embodiments provided in this application can include at least one of non-volatile and volatile memory. Non-volatile memory can include read-only memory (ROM), magnetic tape, floppy disk, flash memory, optical memory, high-density embedded non-volatile memory, resistive random access memory (ReRAM), magnetic random access memory (MRAM), ferroelectric random access memory (FRAM), phase change memory (PCM), graphene memory, etc. Volatile memory can include random access memory (RAM) or external cache memory, etc. By way of illustration and not limitation, RAM can take many forms, such as Static Random Access Memory (SRAM) or Dynamic Random Access Memory (DRAM).

[0116] The databases involved in the embodiments provided in this application may include at least one type of relational database and non-relational database. Non-relational databases may include, but are not limited to, blockchain-based distributed databases. The processors involved in the embodiments provided in this application may be general-purpose processors, central processing units, graphics processing units, digital signal processors, programmable logic devices, quantum computing-based data processing logic devices, etc., and are not limited to these.

[0117] The technical features of the above embodiments can be combined in any way. For the sake of brevity, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0118] This document uses specific examples to illustrate the principles and implementation methods of this application. The descriptions of the above embodiments are only for the purpose of helping to understand the methods and core ideas of this application. Furthermore, those skilled in the art will recognize that, based on the ideas of this application, there will be changes in the specific implementation methods and application scope. Therefore, the content of this specification should not be construed as a limitation of this application.

Claims

1. A sparse LU decomposition solution method based on matrix computation unit optimization, characterized in that, include: Obtain the sparse matrix to be decomposed; The sparse matrix to be decomposed is the sparse matrix in the linear equation to be solved; The MC64 algorithm is called to reorder the sparse matrix to be decomposed, thereby improving data stability and obtaining the sparse matrix after one reordering. The PRCM algorithm is called to reorder the sparse matrix after one reordering to improve the local density of the data and obtain a sparse matrix after bandwidth compression. The bandwidth-compressed sparse matrix is ​​converted into a regular dense Tile block structure and stored in a pre-allocated contiguous memory space; A right-look-based LU decomposition algorithm is used to numerically decompose the bandwidth-compressed sparse matrix to obtain the corresponding LU matrix; matrix computation units are used to accelerate the numerical decomposition of the bandwidth-compressed sparse matrix. Based on the LU matrix, the lower triangular equation system and the upper triangular equation system are solved to achieve an efficient solution to the linear equation to be solved. The PRCM algorithm is called to reorder the sparse matrix after the first reordering to improve the local density of the data, resulting in a bandwidth-compressed sparse matrix. Specifically, this includes: With the goal of minimizing the bandwidth of the sparse matrix, key nodes are identified from the sparse matrix after one reordering; the key nodes are those that affect the bandwidth compression operation; and the degree of the key nodes is greater than a node degree threshold. Move several of the key nodes to the bottom right corner of the sparse matrix after one reordering; The RCM algorithm is used to perform bandwidth compression on the subgraph composed of several nodes remaining after removing several key nodes, resulting in a bandwidth-compressed sparse matrix. The node degree threshold for minimizing the bandwidth of a sparse matrix is ​​determined using the following formula: ； in, t The node degree threshold. T Let be the set of degrees of the nodes in the sparse matrix. This is the minimum estimate of the bandwidth of the sparse matrix. a It is a constant. , S ( t () is a degree greater than t The number of nodes.

2. The sparse LU decomposition solution method based on matrix computation unit optimization according to claim 1, characterized in that, A right-look-based LU decomposition algorithm is used to numerically decompose the bandwidth-compressed sparse matrix to obtain the corresponding LU matrix, specifically including: In the sparse matrix after bandwidth compression, several tile block structures are selected by a sliding window; the size of the sliding window depends on the number of tile block structures that need to be calculated in the current round. For the several tile block structures within the sliding window, perform LU decomposition for the current round to obtain the LU matrix vector for the current round; the LU matrix vector for the current round includes the row vector of the U matrix and the column vector of the L matrix in the current round; The sliding window moves down layer by layer along the bandwidth-compressed sparse matrix, selects several tile block structures, and jumps to "perform LU decomposition of the current round for the several tile block structures in the sliding window to obtain the LU matrix vector of the current round", until the LU decomposition of the bandwidth-compressed sparse matrix is ​​completed; Based on the LU matrix vectors obtained from each iteration, the corresponding LU matrix is ​​obtained.

3. The sparse LU decomposition solution method based on matrix computation unit optimization according to claim 2, characterized in that, For the several tile blocks within the sliding window, perform LU decomposition for the current round to obtain the LU matrix vector for the current round, specifically including: Based on the diagonal blocks of the previous round within the sliding window, the LU matrix vector of the current round is obtained by decomposition. Based on the LU matrix vector of the current round, the diagonal block of the next round is calculated, and the Schur complement is calculated for other Tile block structures within the sliding window. When performing the Schur complement calculation, the Tile block structure is used as the unit, and the matrix calculation unit is used to accelerate the 3-level BLAS operation.

4. The sparse LU decomposition solution method based on matrix computation unit optimization according to claim 3, characterized in that, While the computation stream performs the LU decomposition for the current round, a copy stream is used to store the decomposition results from the previous round and prefetch the data for the next round.

5. A sparse LU decomposition solver based on matrix computation unit optimization, characterized in that, The sparse LU decomposition solving device based on matrix computation unit optimization is used to implement the sparse LU decomposition solving method based on matrix computation unit optimization as described in any one of claims 1-4, wherein the sparse LU decomposition solving device based on matrix computation unit optimization includes: The sparse matrix acquisition module is used to obtain the sparse matrix to be decomposed; the sparse matrix to be decomposed is the sparse matrix in the linear equation to be solved. The sparse matrix reordering module is used to call the MC64 algorithm to reorder the sparse matrix to be decomposed, improve data stability, and obtain the sparse matrix after one reordering. The sparse matrix bandwidth compression module is used to call the PRCM algorithm to reorder the sparse matrix after one reordering, improve the local density of the data, and obtain the bandwidth-compressed sparse matrix. The sparse matrix block storage module is used to convert the bandwidth-compressed sparse matrix into a regular and dense Tile block structure and store it in a pre-allocated contiguous memory space. The sparse matrix numerical decomposition module is used to numerically decompose the bandwidth-compressed sparse matrix using a right-look-based LU decomposition algorithm to obtain the corresponding LU matrix; the matrix calculation unit is used to accelerate the numerical decomposition of the bandwidth-compressed sparse matrix. The linear equation solving module is used to solve the lower triangular equation system and the upper triangular equation system based on the LU matrix, so as to complete the efficient solution of the linear equation to be solved.

6. A computer device, comprising: A memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that the processor executes the computer program to implement the sparse LU decomposition solution method based on matrix computation unit optimization as described in any one of claims 1-4.

7. A computer-readable storage medium having a computer program stored thereon, characterized in that, When executed by a processor, the computer program implements the sparse LU decomposition solution method based on matrix computation unit optimization as described in any one of claims 1-4.

8. A computer program product, comprising a computer program, characterized in that, When executed by a processor, the computer program implements the sparse LU decomposition solution method based on matrix computation unit optimization as described in any one of claims 1-4.

Images