Sparse LU decomposition solving method based on matrix calculation unit optimization and related device

Abstract

The invention discloses a matrix calculation unit optimization-based sparse LU decomposition solving method and a related device, and relates to the technical field of numerical algorithm optimization, and the method comprises the steps: calling an improved PRCM algorithm to carry out the reordering of a sparse matrix, improving the local density of data, obtaining a sparse matrix after bandwidth compression, and obtaining a sparse matrix; then converting the sparse matrix into a regular and dense Tile block structure, storing the Tile block structure in a pre-allocated continuous memory space, performing numerical decomposition on the sparse matrix after bandwidth compression by adopting a right view-based LU decomposition algorithm, and accelerating by using a matrix calculation unit in the numerical decomposition process to obtain an LU matrix; therefore, efficient solving of the to-be-solved linear equation can be completed. According to the scheme, the spatial locality of the sparse matrix is fully utilized, the strong computing power of the matrix computing unit is exerted, the computing efficiency is improved while the computing amount is reduced, and the time of sparse LU decomposition in symbol decomposition and numerical decomposition stages is remarkably shortened.

Description

Technical Field

[0001] The present application relates to the technical field of numerical algorithm optimization, and in particular to a sparse LU decomposition solution method based on matrix computing unit optimization and related devices. Background Art

[0002] Sparse LU decomposition is an important algorithm in numerical linear algebra, mainly used to solve the linear system Ax=b of sparse matrices. A is a sparse matrix with most elements being zero. LU decomposition splits the sparse matrix into the product of a lower triangular matrix (L) and an upper triangular matrix (U), so that A=LU, thereby simplifying the solution process. Compared with directly solving the matrix equation, the advantage of LU decomposition is that it reduces repeated calculations through matrix decomposition, making the subsequent solution process more efficient. Sparse LU decomposition is mainly divided into 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 matrix and the U matrix; numerical decomposition is floating-point calculation to find L and U.

[0003] However, due to its inherent complex data dependencies and irregular distribution of sparse matrices, parallelization algorithms face huge challenges. Complex data dependencies mean that it is difficult to find parts that can be parallelized, and the parallelization granularity is usually small, resulting in a weak acceleration effect. In addition, sparse matrices usually use storage formats such as CSC and CSR, and their calculations themselves require irregular memory access, resulting in additional overhead. The irregular characteristics of sparse matrices limit parallel optimization.

[0004] In summary, the existing sparse matrix LU decomposition technology has 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 dispersed and irregular, resulting in low effective computing efficiency of the matrix computing unit. (2) The data structure is non-uniformly divided, which is not suitable for the matrix computing unit structure. This conflicts with the requirement of uniform and dense input of the matrix computing unit, making it difficult to bring into play the powerful computing 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 that can utilize matrix computation unit optimization is crucial to accelerate the solution of large-scale sparse linear equations and improve its performance in practical applications. Summary of the invention

[0006] The purpose of this application is to provide a sparse LU decomposition solution method and related devices based on matrix calculation unit optimization, which can reduce the amount of calculation while improving the calculation efficiency, and significantly improve the performance of sparse LU decomposition in practical applications.

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

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

[0009] Get 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 to improve data stability and obtain a reordered sparse matrix.

[0011] The PRCM algorithm is called to reorder the sparse matrix after the first reordering to improve the local density of the data and obtain the sparse matrix after bandwidth compression.

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

[0013] A right-view based LU decomposition algorithm is used to numerically decompose the bandwidth compressed sparse matrix to obtain a corresponding LU matrix; when numerically decomposing the bandwidth compressed sparse matrix, a matrix calculation unit is used for acceleration.

[0014] The lower triangular equations and the upper triangular equations are solved based on the LU matrix to efficiently solve the linear equations to be solved.

[0015] Optionally, calling the PRCM algorithm to reorder the sparse matrix after the first reordering to improve the local density of the data and obtain a sparse matrix after bandwidth compression, specifically includes the following steps:

[0016] With the goal of minimizing the bandwidth of the sparse matrix, key nodes are found from the sparse matrix after a reordering; key nodes are nodes that have an impact on the bandwidth compression operation.

[0017] Move several key nodes to the lower right corner of the reordered sparse matrix.

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

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

[0020]

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

[0022] Optionally, a right-view based LU decomposition algorithm is used to numerically decompose the sparse matrix after bandwidth compression to obtain a corresponding LU matrix, which specifically includes the following steps:

[0023] In the sparse matrix after bandwidth compression, several tile structures are selected through a sliding window. The size of the sliding window depends on the number of tile structures that need to be calculated in the current round.

[0024] For several tile block structures in the sliding window, the LU decomposition of the current round is performed to obtain the LU matrix vector of the current round. The LU matrix vector of the current round includes the row vector of the U matrix in the current round and the column vector of the L matrix in the current round.

[0025] Move the sliding window downward layer by layer along the sparse matrix after bandwidth compression, select several Tile block structures, and jump to "for several Tile block structures in the sliding window, perform LU decomposition of the current round to obtain the LU matrix vector of the current round" until the LU decomposition of the sparse matrix after bandwidth compression is completed.

[0026] According to the LU matrix vector obtained in each round of iteration, the corresponding LU matrix is ​​obtained.

[0027] Optionally, for a plurality of Tile block structures in the sliding window, LU decomposition of the current round is performed to obtain the LU matrix vector of the current round, which specifically includes the following steps:

[0028] According to the diagonal blocks of the previous round in the sliding window, the LU matrix vector of the current round is decomposed.

[0029] According to the LU matrix vector of the current round, the diagonal block of the next round is calculated, and the Schur's complement calculation is performed on other Tile block structures in the sliding window; when performing the Schur's 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 is performing the current round of LU decomposition, the copy stream is used to store the decomposition result completed in the previous round and complete pre-fetching of the next round of data.

[0031] In a second aspect, the present application provides a sparse LU decomposition solving device based on matrix computing 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 a reordered sparse matrix.

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

[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 continuous memory space.

[0036] The sparse matrix numerical decomposition module is used to use the right-view-based LU decomposition algorithm to numerically decompose the bandwidth-compressed sparse matrix to obtain the corresponding LU matrix; when numerically decomposing the bandwidth-compressed sparse matrix, a matrix calculation unit is used for acceleration.

[0037] The linear equation solving module is used to solve the lower triangular equation group and the upper triangular equation group based on the LU matrix, so as to efficiently solve the linear equation to be solved.

[0038] In a third aspect, the present application provides a computer device comprising: 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 computing unit optimization described above.

[0039] In a fourth aspect, the present 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 computing unit optimization described above.

[0040] In a fifth aspect, the present application provides a computer program product, including a computer program, which, when executed by a processor, implements the steps of the sparse LU decomposition solution method based on matrix computing unit optimization described above.

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

[0042] The present application provides a sparse LU decomposition solution method based on matrix computing unit optimization and a related device. In the method, the MC64 algorithm is first called to reorder the sparse matrix to be decomposed to improve data stability and obtain a reordered sparse matrix. Then, the PRCM algorithm is called to reorder the sparse matrix after the reordering to improve the local density of the data and obtain a bandwidth-compressed sparse matrix. Then, the bandwidth-compressed sparse matrix is ​​converted into a regular dense Tile block structure and stored in a pre-allocated continuous memory space. The right-view-based LU decomposition algorithm is used to numerically decompose the bandwidth-compressed sparse matrix. During the numerical decomposition process, the matrix computing unit is used for acceleration to obtain the corresponding LU matrix. Finally, the lower triangular equation group and the upper triangular equation group can be solved based on the LU matrix to complete the efficient solution of the linear equation to be solved. The present application scheme adopts an improved PRCM algorithm for reordering, which can solve the problems of failure of some RCM algorithms, and further enhance the local density and compression bandwidth of data, reduce the amount of calculation and improve the calculation efficiency. The regular dense block data structure used after the reordering is completed can directly call the matrix calculation unit for acceleration in the subsequent numerical decomposition stage, and also eliminates the indirect memory access overhead caused by format conversion in traditional methods, fully utilizes the spatial locality of sparse matrices, and at the same time exerts the powerful computing power of the matrix calculation unit, significantly reducing the time of sparse LU decomposition in the symbolic decomposition and numerical decomposition stages. BRIEF DESCRIPTION OF THE DRAWINGS

[0043] In order to more clearly illustrate the embodiments of the present application or the technical solutions in the prior art, the drawings required for use in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present application. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0044] Figure 1 A flowchart of a sparse LU decomposition solution method based on matrix computing unit optimization provided in one embodiment of the present application.

[0045] Figure 2 A schematic diagram comparing RCM reordering results and PRCM reordering results in a sparse LU decomposition solution method based on matrix computing unit optimization provided in one embodiment of the present application.

[0046] Figure 3 A flowchart of step A3 in a sparse LU decomposition solution method based on matrix computing unit optimization provided in one embodiment of the present application.

[0047] Figure 4 A schematic diagram of the process of symbolic decomposition of a GaAsH6 matrix in a sparse LU decomposition solution method based on matrix computing unit optimization provided in one embodiment of the present application.

[0048] Figure 4 (a) is a schematic diagram of the original GaAsH6 matrix in a sparse LU decomposition solution method based on matrix calculation unit optimization provided in one embodiment of the present 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 computing unit optimization provided in an embodiment of the present 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 computing unit optimization provided in one embodiment of the present application.

[0051] Figure 5 A flowchart of step A5 in a sparse LU decomposition solution method based on matrix computing unit optimization provided in one embodiment of the present application.

[0052] Figure 6 A schematic diagram of the execution process of numerical decomposition in a sparse LU decomposition solution method based on matrix computing unit optimization provided in one embodiment of the present application.

[0053] Figure 7 A schematic diagram of the execution timeline of numerical decomposition in a sparse LU decomposition solution method based on matrix computing unit optimization provided in one embodiment of the present application.

[0054] Figure 8 A schematic diagram of the dependency relationship between data when performing numerical decomposition in a sparse LU decomposition solution method based on matrix computing unit optimization provided in one embodiment of the present application.

[0055] Fig. 9 A schematic diagram comparing the sparse LU decomposition solution method based on matrix computing unit optimization provided in one embodiment of the present application with the traditional SuperLU_DIST and PanguLU test results.

[0056] Fig.10 A schematic diagram of the functional modules of a sparse LU decomposition solving device based on matrix computing unit optimization provided in one embodiment of the present application.

[0057] Fig.11 A schematic diagram of the structure of a computer device provided in one embodiment of the present application. DETAILED DESCRIPTION

[0058] The following will be combined with the drawings in the embodiments of the present application to clearly and completely describe the technical solutions in the embodiments of the present application. Obviously, the described embodiments are only part of the embodiments of the present application, not all of the embodiments. Based on the embodiments in the present application, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of this application.

[0059] In order to make the above-mentioned objects, features and advantages of the present application more obvious and easy to understand, the present application is further described in detail below with reference to the accompanying drawings and specific implementation methods.

[0060] The embodiment of the present application provides a sparse LU decomposition solution method based on matrix calculation unit optimization, such as Figure 1 As shown, the following steps are included:

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

[0062] A2. Call the MC64 algorithm to reorder the sparse matrix to be decomposed to improve data stability and obtain a reordered sparse matrix. In this embodiment, before bandwidth compression, the MC64 algorithm is first used to perform static pivot selection to improve numerical stability.

[0063] A3. Call the PRCM algorithm to reorder the sparse matrix after the first reordering to improve the local density of the data and obtain the sparse matrix after bandwidth compression. Use the improved PRCM algorithm to reorder again to gather the non-zero elements in the sparse matrix and reduce 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 on the basis of the RCM algorithm. The idea of ​​the PRCM reordering algorithm is to move some high-degree nodes to the lower right corner of the matrix, and then apply the RCM algorithm to the remaining nodes. The function of the RCM algorithm is to gather the non-zero elements of the matrix close to the diagonal, which plays a role in compressing the matrix bandwidth. However, the RCM algorithm fails in some cases, such as matrices in the field of circuit simulation. This is because these matrices usually contain nodes with particularly large degrees. These nodes amplify the bandwidth of the entire sparse matrix, resulting in the inability of non-zero elements to gather near the diagonal.

[0065] The PRCM algorithm first finds out the nodes with extremely large degrees that have an impact on compression, and calls these nodes "key nodes". Then, it moves these nodes to the lower right corner of the matrix, and then uses RCM to compress the bandwidth of the subgraph composed of the remaining nodes. Placing the key nodes in the lower right corner can reduce the number of new non-zero elements. During the Schur complement calculation process, no new non-zero elements will be added outside the diagonal and the lower right corner. Without the influence of the key nodes, the remaining subgraph can be well compressed by the RCM algorithm. Figure 2 As shown, the left side shows the result of RCM reordering, and the right side shows the result of PRCM reordering. The PRCM algorithm moves the key nodes to the lower right corner, forming a downward triangle 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 composed of the remaining nodes after RCM reordering. In order to achieve better bandwidth compression effect of PRCM reordering, B' <B。

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

[0067] A31. With the goal of minimizing the bandwidth of the sparse matrix, find the key nodes from the sparse matrix after a reordering; the key nodes are the nodes that have an impact on the bandwidth compression operation.

[0068] A32. Move several key nodes to the lower right corner of the reordered sparse matrix.

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

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

[0071] make T is a 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. It is necessary to select a suitable threshold from these possible thresholds to minimize the bandwidth B. The general idea is to use a traversal method to calculate each t in T, and then select a t that minimizes B. This is the threshold to minimize the bandwidth. Finally, all nodes with degrees greater than t are moved to the lower right corner of the matrix, and the remaining nodes are reordered using RCM. However, it takes too much time to calculate each t using the RCM algorithm, which is unacceptable. Therefore, in this embodiment, an estimation method is proposed, where β'(t) = at, where a is a constant, and a set of possible values ​​of a is taken as Γ = {2 0 ,2 1 ,…,2 14 The bandwidth obtained by using the estimated value at+S(t) is called the estimated bandwidth B * , we get a set of functions B for estimating bandwidth * =at+S(t),a∈Γ. Next, we need to find B from this group * Find the minimum value in the function The minimum value of the estimated bandwidth is used to approximate the minimum value of the actual bandwidth. At this time, 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 the smallest can be determined according to the following formula:

[0073]

[0074] Among them, t is the node degree threshold, T is the degree set of nodes in the sparse matrix, is the minimum estimate of the sparse matrix bandwidth, a is a constant, Γ={2 0 ,2 1 ,…,2 14}, S(t) is the number of nodes with degree greater than t. After actual testing, it is found that the difference between the estimated bandwidth and the true value does not exceed 1%, and the time complexity of the bandwidth estimation method is much better than the actual bandwidth calculation method.

[0075] By optimizing the reordering algorithm, the spatial structure of the matrix is ​​optimized, so that the non-zero elements are more concentrated around the diagonal, so that the number of newly added non-zero elements is effectively controlled during the Schur complement calculation. Compressing the bandwidth of the matrix reduces the amount of calculation and saves GPU memory, 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 structures that are easy to calculate.

[0076] A4. Convert the bandwidth compressed sparse matrix into a regular dense Tile block structure and store it in a pre-allocated continuous 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, convert the CSC format of the sparse matrix into a regular dense block structure for storage. These blocks are called Tiles. Traditional solvers determine the LU structure before dividing the blocks, while in this embodiment, the blocks are first divided and the structure is determined. Only the position of the newly added blocks needs to be determined. Compared with symbolic decomposition in element units, symbolic decomposition in blocks greatly reduces the time overhead and does not require complex structures and algorithms such as elimination trees in traditional methods.

[0078] like Figure 4 As shown, this example is the symbolic decomposition of the GaAsH6 matrix. Figure 4 (a) is the original matrix, Figure 4 (b) is the result after reordering the matrix and partitioning. Figure 4 (c) is the final structure determined after symbolic decomposition. Since sparse matrices occupy a large amount of memory, continuous storage space is opened up in the memory to store them, and the scattered data is converted into continuous data storage, which saves the CPU-GPU memory copy time.

[0079] A5. Use the right-view based LU decomposition algorithm to numerically decompose the sparse matrix after bandwidth compression to obtain the corresponding LU matrix; when numerically decomposing the sparse matrix after bandwidth compression, use the matrix calculation unit for acceleration.

[0080] The numerical decomposition process in the sparse LU decomposition algorithm is the part that takes up the most time. This step is directly calculated based on the block structure of the previous step and accelerated by the matrix calculation unit. The evenly divided dense block storage used in this embodiment not only avoids unnecessary format conversion overhead, but also can call level 3 BLAS and use the matrix calculation unit to accelerate matrix operations.

[0081] In this embodiment, the numerical decomposition algorithm used is a right-looking LU decomposition, which optimizes the dependency of the diagonal blocks on the L column blocks and the U row blocks in each iteration, simplifies the two layers of dependency into one layer, and uses a sliding window strategy to cover up the delays of calculation and communication, thereby enhancing data reusability and improving data-level parallelism and computing memory access ratio.

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

[0083] A51. In the sparse matrix after bandwidth compression, several tile block structures are selected through 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. In each round, all tiles to be calculated constitute a sliding window. The sliding window consists of three parts: the part of the result of the previous round of calculation that is no longer used, which must be transmitted back to the CPU; the tiles to be calculated in this round (also including some tiles used in the previous round of calculation); and the tiles to be used in the next round of calculation, data pre-fetching. These three parts together constitute a sliding window. During numerical decomposition, the sliding window method is used to complete the data return and pre-fetching while performing floating-point calculations, forming a mutual concealment of data transmission and calculation, and the reused parts are still retained inside the window, which improves the reusability of the data.

[0084] For example, in the kth iteration, the equivalent bandwidth of the matrix is ​​B k , the dimension of the Tile block structure is ts (the size of a Tile block 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 a plurality of tile block structures in the sliding window, perform LU decomposition of the current round to obtain the LU matrix vector of the current round. The step A52 specifically includes the following steps:

[0086] According to the diagonal blocks of the previous round in the sliding window, the LU matrix vector of the current round is decomposed. The LU matrix vector of the current round includes the row vector of the U matrix in the current round and the column vector of the L matrix in the current round.

[0087] According to the LU matrix vector of the current round, the diagonal block of the next round is calculated, and the Schur's complement calculation is performed on other Tile block structures in the sliding window; when performing the Schur's 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, moving the sliding window downward layer by layer along the sparse matrix after bandwidth compression, selecting several Tile block structures, and jumping to step A52 until the LU decomposition of the sparse matrix after bandwidth compression is completed.

[0089] A54. Obtain the corresponding LU matrix according to the LU matrix vector obtained in each round of iteration.

[0090] In this embodiment, in the above-mentioned numerical decomposition process, independent calculation streams and copy streams are used to achieve decoupling of calculation and data movement, reduce the dependence of calculation on data, and improve the efficiency of data transmission and calculation efficiency. That is, while the calculation stream performs the current round of LU decomposition, the copy stream is used to store the decomposition results completed in the previous round and complete the pre-fetching of the next round of data.

[0091] Specifically, the specific implementation process is as follows Figure 6 As shown, for the kth iteration, in the LU decomposition (light gray part, according to the diagonal block D of the k-1th round k-1 Decompose to get the LU matrix vectors L and U of the current round, and calculate the diagonal block D of the kth round k , and perform Schur complement calculation S on other tile block structures in the sliding window, use the copy stream to copy the decomposition result (dark gray part) completed in the previous round from the GPU back to the CPU. This part is no longer used in subsequent calculations. At the same time, the data needed in the k+1 round is copied from the CPU to the GPU to complete data prefetching. The data is transmitted bidirectionally. Figure 6 D2H means data transmission from GPU to CPU, and H2D means data transmission from CPU to GPU.

[0092] The timeline of the execution process is as follows Figure 7 As shown, the copy stream is Figure 7 In LU decomposition, the computation flow is divided into three parts: diagonal block decomposition (D), decomposition of L columns and U rows (L, U), and Schur complement calculation (S); the first computation flow solves one column of L and one row of U, that is, Figure 7 The second calculation flow decomposes the diagonal blocks in the next round of calculation, that is, Figure 7 Diag in; the third computational flow, Schur complement computation, is Figure 7 In the Schur calculation, 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 k-1 Do GETRF (diagonal block decomposition), then decompose L columns and U rows at the same time, and finally decompose D k Do GEMMTRF (Schur complement calculation + diagonal block decomposition) and S do Schur complement calculation.

[0094] After each round of iteration, the sliding window composed of copying and calculation moves once. Most of the blocks in the window are reused, with a high data reuse rate, which can significantly reduce cache misses. The calculation unit of Shul complement is in blocks, and the matrix calculation unit is used to accelerate the 3rd level BLAS operation, which improves the granularity of data parallelism and improves parallelism. Because the PRCM reordering gathers non-zero elements near the diagonal, the size of the sliding window that needs to be calculated in each round of iteration is relatively small, which saves GPU memory usage and can be calculated simultaneously in the GPU.

[0095] The PCIe bus is full-duplex, and data transmission between D2H and H2D can be performed simultaneously, with roughly the same amount of data being transmitted. According to the definition of the computational copy ratio, the computational copy ratio of the sliding window is only related to the bandwidth of the current position. When the CCR of the sliding window is greater than the CCR of the GPU, the computation can be completely covered by the communication. For the NVIDIA A800 GPU, its Theoretically, when the bandwidth B>624×8=4992, the communication time can cover the computing time. Therefore, it has better performance in larger sparse matrices.

[0096] Before each round of numerical decomposition, the L columns and U rows of the results of the previous round of decomposition are first transferred from the GPU to the CPU, and the tiles needed for the next round are pre-fetched and transferred from the CPU to the GPU. While the data is being transmitted, the L columns and U rows are decomposed and updated with Schur's complement. After the above numerical decomposition is completed, the original sparse matrix A becomes the product of the L matrix and the U matrix, that is, A = LU.

[0097] A6. Solve the lower triangular equation group Ly=b and the upper triangular equation group Ux=y based on the LU matrix to complete the efficient solution of the linear equations to be solved.

[0098] In order to evaluate the advantages of the sparse LU decomposition solution based on matrix computing unit optimization proposed in this application, 165 large-scale sparse matrices with dimensions ranging from 50K to 2.7M were selected from the SuiteSparse Matrix Collection data set. The compared solvers are the two most advanced solvers in the industry, SuperLU_DIST and PanguLU. In the most time-consuming numerical decomposition stage, compared with SuperLU_DIST, the solution proposed in this application achieved an acceleration ratio of 17.9 times and a maximum acceleration ratio of 86.2 times; compared with PanguLU, the solution proposed in this application achieved an average acceleration ratio of 8.8 times and a maximum acceleration ratio of 90.9 times. The test results are as follows Fig. 9 shown.

[0099] The sparse LU decomposition solution method based on matrix calculation unit optimization provided in the above embodiment of the present application can solve some situations where the RCM algorithm fails by improving the RCM algorithm, and further enhances the local density and compression bandwidth of the data, reduces the amount of calculation, and improves the calculation efficiency. The reordering algorithm used in this application is different from the traditional solver. It no longer aims to minimize the newly added non-zero elements, but enhances the locality of the data, so that the data reuse is higher in the sliding window, and the compression matrix bandwidth brings lower GPU video memory usage. At the same time, due to the use of a uniform dense block data structure, the 3-level BLAS interface accelerated by the matrix calculation unit can be directly called in the numerical decomposition stage; using uniform blocks, these blocks are directly operated, eliminating the indirect memory access overhead caused by format conversion in the traditional method. This application is implemented based on the NVIDIA GPU platform, making full use of the Tensor Core components, and achieving an average peak performance of 22.5%, while the traditional solver only exerts less than 1% of the peak performance.

[0100] Based on the same inventive concept, the embodiment of the present application also provides a device for implementing the above-mentioned sparse LU decomposition solution method based on matrix calculation unit optimization. The implementation scheme for solving the problem provided by the device is similar to the implementation scheme recorded in the above method, so the specific limitations in one or more embodiments of the sparse LU decomposition solution device based on matrix calculation unit optimization provided below can refer to the limitations of the sparse LU decomposition solution method based on matrix calculation unit optimization in the above text, and will not be repeated here.

[0101] In an exemplary embodiment, Fig.10 As shown, a sparse LU decomposition solution device based on matrix computing unit optimization is provided, including 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 a reordered sparse matrix.

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

[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 continuous memory space.

[0106] The sparse matrix numerical decomposition module is used to use the right-view-based LU decomposition algorithm to numerically decompose the bandwidth-compressed sparse matrix to obtain the corresponding LU matrix; when numerically decomposing the bandwidth-compressed sparse matrix, a matrix calculation unit is used for acceleration.

[0107] The linear equation solving module is used to solve the lower triangular equation group and the upper triangular equation group based on the LU matrix, so as to efficiently solve the linear equation to be solved.

[0108] certainly, Fig.10 The architecture shown is only exemplary and can be omitted according to actual needs when implementing different functions. Fig.10 One or at least two components of the system shown.

[0109] In an exemplary embodiment, a computer device is provided. The computer device may be a server or a terminal. The internal structure diagram thereof may be as follows: Fig.11 As shown. The computer device includes a processor, a memory, an input / output interface (Input / Output, referred to as I / O) and a communication interface. The processor, the memory and the input / output interface are connected through a system bus, and the communication interface is connected to the system bus through the input / output interface. The processor of the computer device is used to provide computing and control capabilities. The memory of the computer device includes a non-volatile storage medium and an internal memory. The non-volatile storage medium stores an operating system, a computer program and a database. The internal memory provides an environment for the operation of the operating system and the computer program in the non-volatile storage medium. The input / output interface of the computer device is used to exchange information between the processor and an external device. The communication interface of the computer device is used to communicate with an external terminal through a network connection. When the computer program is executed by the processor, the steps of the sparse LU decomposition solution method based on matrix calculation unit optimization provided in the above embodiment can be implemented.

[0110] Those skilled in the art will understand that Fig.11 The structure shown in the figure is only a block diagram of a part of the structure related to the solution of the present application, and does not constitute a limitation on the computer device to which the solution of the present application is applied. The specific computer device may include more or fewer components than those shown in the figure, or combine certain components, or have a different arrangement of components.

[0111] In an exemplary embodiment, a computer device is further provided, including a memory and a processor, wherein a computer program is stored in the memory, and the processor implements the steps in the above-mentioned method embodiments when executing the computer program.

[0112] In an exemplary embodiment, a computer-readable storage medium is provided, storing a computer program, and when the computer program is executed by a processor, the steps in the above method embodiments are implemented.

[0113] In an exemplary embodiment, a computer program product is provided, including a computer program, and when the computer program is executed by a processor, the steps in the above method embodiments are implemented.

[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, stored data, displayed data, 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 relevant data must comply with relevant regulations.

[0115] Those of ordinary skill in the art can understand that all or part of the processes in the above-mentioned embodiment methods can be completed by instructing the relevant hardware through a computer program, and the computer program can be stored in a non-volatile computer-readable storage medium. When the computer program is executed, it can include the processes of the embodiments of the above-mentioned methods. Among them, any reference to the memory, database or other medium used in the embodiments provided in the present 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 may be in various forms, such as static random access memory (SRAM) or dynamic random access memory (DRAM).

[0116] The database involved in each embodiment provided in this application may include at least one of a relational database and a non-relational database. The non-relational database may include a distributed database based on blockchain, etc., but is not limited thereto. The processor involved in each embodiment provided in this application may be a general-purpose processor, a central processing unit, a graphics processor, a digital signal processor, a programmable logic device, a data processing logic device based on quantum computing, etc., but is not limited thereto.

[0117] The technical features of the above embodiments may be combined arbitrarily. To make the description concise, 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 article uses specific examples to illustrate the principles and implementation methods of this application. The description of the above embodiments is only used to help understand the method and core ideas of this application. At the same time, for those skilled in the art, according to the ideas of this application, there will be changes in the specific implementation methods and application scope. In summary, the content of this specification should not be understood as limiting this application.

Claims

1. A sparse LU decomposition solution method based on matrix computing unit optimization, characterized in that: include: Get the sparse matrix to be decomposed; The sparse matrix to be decomposed is a sparse matrix in the linear equation to be solved; Call the MC64 algorithm to reorder the sparse matrix to be decomposed to improve data stability and obtain a reordered sparse matrix; Call the PRCM algorithm to reorder the sparse matrix after the first reordering to improve the local density of the data and obtain a sparse matrix after bandwidth compression; Convert the bandwidth-compressed sparse matrix into a regular dense tile structure and store it in a pre-allocated continuous memory space; The LU decomposition algorithm based on right view is used to numerically decompose the sparse matrix after bandwidth compression to obtain the corresponding LU matrix; when the sparse matrix after bandwidth compression is numerically decomposed, the matrix calculation unit is used for acceleration; The lower triangular equations and the upper triangular equations are solved based on the LU matrix, thereby completing efficient solution of the linear equations to be solved.

2. The sparse LU decomposition solution method based on matrix computing unit optimization according to claim 1 is characterized in that: The PRCM algorithm is called to reorder the sparse matrix after the first reordering to improve the local density of the data and obtain the sparse matrix after bandwidth compression, which includes: With the goal of minimizing the bandwidth of the sparse matrix, key nodes are found from the sparse matrix after a reordering; the key nodes are nodes that have an impact on the bandwidth compression operation; Moving a plurality of the key nodes to the lower right corner of the reordered sparse matrix; The RCM algorithm is used to perform bandwidth compression operation on a subgraph consisting of several nodes remaining after removing several key nodes, so as to obtain a sparse matrix after bandwidth compression.

3. The sparse LU decomposition solution method based on matrix calculation unit optimization according to claim 2 is characterized in that: The degree of the key node is greater than the node degree threshold; the node degree threshold when the sparse matrix bandwidth is minimized is determined according to the following formula: Among them, t is the node degree threshold, T is the degree set of nodes in the sparse matrix, is the minimum estimate of the sparse matrix bandwidth, a is a constant, Γ={2 0 ,2 1 ,…,2 14 }, S(t) is the number of nodes with degree greater than t.

4. The sparse LU decomposition solution method based on matrix computing unit optimization according to claim 1 is characterized in that: The LU decomposition algorithm based on right view is used to numerically decompose the sparse matrix after bandwidth compression to obtain the corresponding LU matrix, which includes: In the sparse matrix after bandwidth compression, several tile block structures are selected through 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 multiple tile block structures in the sliding window, perform LU decomposition of the current round to obtain the LU matrix vector of the current round; the LU matrix vector of the current round includes the row vector of the U matrix in the current round and the column vector of the L matrix in the current round; Move the sliding window downward layer by layer along the sparse matrix after bandwidth compression, select several Tile block structures, and jump to "for several Tile block structures in the sliding window, perform LU decomposition of the current round to obtain the LU matrix vector of the current round", until the LU decomposition of the sparse matrix after bandwidth compression is completed; According to the LU matrix vector obtained in each round of iteration, the corresponding LU matrix is ​​obtained.

5. The sparse LU decomposition solution method based on matrix computing unit optimization according to claim 4 is characterized in that: For the several Tile block structures in the sliding window, the LU decomposition of the current round is performed to obtain the LU matrix vector of the current round, which specifically includes: Decomposing the diagonal blocks of the previous round in the sliding window to obtain the LU matrix vector of the current round; According to the LU matrix vector of the current round, the diagonal block of the next round is calculated, and the Schur's complement calculation is performed on other Tile block structures in the sliding window; when performing the Schur's complement calculation, the Tile block structure is used as a unit, and the matrix calculation unit is used to accelerate the 3-level BLAS operation.

6. The sparse LU decomposition solution method based on matrix computing unit optimization according to claim 5, characterized in that: While the computation stream is performing the current round of LU decomposition, the copy stream is used to store the decomposition results completed in the previous round and to complete the pre-fetching of the next round of data.

7. A sparse LU decomposition solving device based on matrix computing unit optimization, characterized in that: include: A sparse matrix acquisition module, used to obtain a sparse matrix to be decomposed; The sparse matrix to be decomposed is a 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 a reordered sparse matrix; The sparse matrix bandwidth compression module is used to call the PRCM algorithm to re-sort the sparse matrix after the first re-sorting, so as to improve the local density of the data and obtain the sparse matrix after the bandwidth compression; 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 continuous memory space; The sparse matrix numerical decomposition module is used to numerically decompose the sparse matrix after bandwidth compression by using the right-view LU decomposition algorithm to obtain the corresponding LU matrix; when the sparse matrix after bandwidth compression is numerically decomposed, the matrix calculation unit is used for acceleration; The linear equation solving module is used to solve the lower triangular equation group and the upper triangular equation group based on the LU matrix, so as to efficiently solve the linear equation to be solved.

8. A computer device comprising: 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 sparse LU decomposition solution method based on matrix computing unit optimization described in any one of claims 1 to 6.

9. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the sparse LU decomposition solution method based on matrix computing unit optimization described in any one of claims 1 to 6 is implemented.

10. A computer program product, comprising a computer program, characterized in that When the computer program is executed by a processor, the sparse LU decomposition solution method based on matrix computing unit optimization described in any one of claims 1 to 6 is implemented.

Images