A GPU-based batch parallel LU decomposition method for small and medium-sized dense matrices for simulation systems

Abstract

The present invention implements a high-performance batch parallel LU decomposition method based on GPU for small and medium-sized dense matrices of a simulation system. The method includes: S1, calculating the Jacobian matrix of the model unit of the simulation system at each time step, and taking these matrices as sub-matrices, and splicing them into a large matrix H; S2, dividing the internal blocks of each sub-matrix, and allocating a GPU thread group to each sub-matrix according to the size of the blocks; S3, the GPU thread group sequentially completes the selection of the main element of all columns in each diagonal block, the update of the right block of the row where the diagonal block is located, the exchange of the internal column of the diagonal block with the main element column, and the update of all the block data below and below the column where the diagonal block is located in a diagonal order from the upper left to the lower right, thereby completing the LU decomposition of all sub-matrices in the large matrix H. The present invention can increase the performance in an approximately linear trend as the batch number increases, and can accelerate the solution speed when implicitly solving a model containing a large batch of small ordinary differential equations, and effectively shorten the simulation time.

Description

Technical Field

[0001] The present invention belongs to the field of cross-technology application of modeling simulation and parallel computing, and in particular to a batch parallel LU decomposition (Lower-Upper Decomposition) method for small and medium-sized dense matrices based on a GPU (Graphics Processing Unit) in complex system simulation. Background Art

[0002] LU decomposition is the most important computational step in many engineering and scientific computing problems. It has a wide range of applications in system simulation, machine learning, image processing, information encryption, and high-speed information transmission. With the continuous development of computer science, the amount of data that needs to be processed is growing exponentially. How to use computers to process experimental data efficiently and quickly has become a difficult problem. In most key applications, what needs to be solved is the parallel solution of a large number of small-scale matrices, that is, to complete a large number of matrix LU decompositions under a set of operation logic, rather than serially solving a few large linear systems.

[0003] Complex systems usually contain a large number of small and medium-sized nonlinear simulation model units, which need to be solved approximately in linear form. The most time-consuming step is the calculation of the linear equations, which requires the LU decomposition of the Jacobian matrix of the model unit. The current GPU-based LU decomposition method is mainly aimed at large matrices. It uses a block form to distribute the processing tasks to multiple threads. It is a method of dividing the internal computing tasks of a large model. However, this method is not suitable for small and medium-sized models. Since the Jacobian matrix is ​​small, the time required to read and write matrix elements is much longer than the calculation time. Matrix operations based on block form cannot increase the proportion of effective computing time, and global memory access delays constitute the main performance bottleneck. The smaller the matrix, the more serious this problem is. Summary of the invention

[0004] To solve the above technical problems, the main purpose of the present invention is to provide a batch parallel LU decomposition method for small and medium-sized dense matrices based on GPU for simulation systems, which can effectively reduce the time consumed by the decomposition process and improve computing efficiency.

[0005] The present invention proposes a batch LU decomposition method for small and medium-sized dense matrices based on GPU for a simulation system, comprising:

[0006] S1, calculate the Jacobian matrix of the model unit of the simulation system at each time step, and use these matrices as sub-matrices and concatenate them into a large matrix H;

[0007] S2, divide each sub-matrix into blocks, and allocate GPU thread groups to each sub-matrix according to the size of the blocks;

[0008] S3. The GPU thread group completes the pivot selection of all columns in each diagonal block, the update of the blocks on the right side of the row where the diagonal block is located, the exchange of the columns inside the diagonal block with the pivot columns, and the update of all block data below and to the lower right of the column where the diagonal block is located in a diagonal order from the upper left to the lower right, thereby completing the LU decomposition of all sub-matrices in the large matrix H.

[0009] Furthermore, in step S1, the simulation system has n model units, and the Jacobian matrix is ​​an m-order small and medium-sized dense matrix H i , i ranges from [0, n-1], where m is obtained by expanding the dimensions of all model units currently to be solved upward to an integer power of 2 m (i.e., the Jacobian matrix is ​​expanded to order m), and for the added rows and columns, except for the diagonal elements whose values ​​are 1, all others are filled with 0.

[0010] Furthermore, in step S2, it includes:

[0011] S21, sub-matrix H i Divide into several sub-blocks to ensure that in subsequent processing, the thread group can read the data in the entire block from the global memory at one time, reducing the number of global memory accesses;

[0012] S22, according to the size of the block for the sub-matrix H i Allocate the corresponding number of GPU threads to form a thread group Group i .

[0013] Furthermore, the submatrix H i Methods for dividing into several sub-blocks include:

[0014] If m is greater than or equal to 16, set the block size blocksize to 8 or 16;

[0015] If m is less than 16, the entire submatrix is ​​taken as a block, and the block size blocksize=m.

[0016] Furthermore, in step S3, the specific steps include:

[0017] S31, each thread group Group i Select submatrix H i The diagonal block B with inner index j j,j , determine B j,j The pivot columns corresponding to all columns in B j,j All blocks on the right row are updated; j is initialized to 0;

[0018] S32, all the obtained pivot columns and diagonal blocks B j,jInterchange the positions of the corresponding columns inside;

[0019] S33. Update all the sub - blocks below the column where the diagonal block B j,j is located;

[0020] S34. Update all the sub - blocks below and to the right of the diagonal block B j,j ;

[0021] S35. If the diagonal block B j,j is the bottom - right sub - block in the sub - matrix H i , end the batch parallel LU decomposition, and output the large matrix H. Among them, each sub - matrix H i is composed of L i and U i , which are the unit lower triangular matrix and the upper triangular matrix respectively, and are used for the subsequent fast solution of linear equations; otherwise, let j = j + 1, jump to step S31, and continue to process the next sub - block.

[0022] Furthermore, step S31 includes:

[0023] S311: Initialize a one - dimensional matrix To i with dimension m, and let To i [l]=m, 0 ≤ l < m, where

[0024] the array index l corresponds to the column number of the sub - matrix H i , and To i [l] is used to record the final pivot column number of the column with serial number l;

[0025] S312: The thread group starts from the current block B j,j , and reads all the elements in the p - th row of the block in parallel and stores them in the cache inside the thread; where p is calculated according to the row index of the sub - matrix H i , and is initialized to j×blocksize; in this step, since the row elements of the block B j,j are closely arranged in the global memory, the GPU only needs to issue one instruction to achieve the reading of all the data of the thread group, achieving the effect of memory coalesced access; where j×blocksize ≤ p < (j + 1)×blocksize);

[0026] S313: The thread group continues to read all the data in the p - th row of the next adjacent sub - block on the right side of B j,j in parallel, compare it with the data in the cache, and save the element with the larger absolute value and its column number in the internal cache; this process continues to iterate until the data of the last block of the sub - matrix is read;

[0027] S314: Compare all elements cached in the thread group, take the column q corresponding to the element with the largest absolute value as the main element column corresponding to column p, and update the element To i [p] = q;

[0028] S315: Pair matrix H i The current row p in the , starting from column index j × blocksize, traverses all subsequent elements in the row If To i [l]>q, then the following operation is performed:

[0029]

[0030] in, For elements The corresponding true pivot;

[0031] S316: For each element in the P1 region The column index y of To i [y]>q, then the following elimination operation is performed:

[0032]

[0033] Among them, the P1 region is located in the submatrix H i A matrix with [p+1, (j+1) × blocksize-1] rows and [j × blocksize, m-1] columns; p+1≤x<(j+1) × blocksize, j × blocksize≤y <m;

[0034] S317: If p is block B j,j The last row in the step S31 ends, and B j,j Contains the LU decomposition result L in the current block j,j and U j,j ; Otherwise, set p=p+1, jump to step S312, and continue processing the next line.

[0035] Further, in step S32, for B j,j Any column p in the thread group, the pth thread in the thread group executes column col p With column col To i [p] Exchange elements between .

[0036] Further, in step S33, let P2 region be located in submatrix H i For each element of the matrix between [(j+1)×blocksize,m-1] rows and [j×blocksize,(j+1)×blocksize-1] columns Perform the following operations:

[0037]

[0038] where \((j + 1)\times blocksize\leq x\lt m\) and \(j\times blocksize\leq y\lt (j + 1)\times blocksize\).

[0039] Furthermore, in step S34, let the P3 region be the matrix between the \([(j + 1)\times blocksize, m - 1]\) - th row and the \([(j + 1)\times blocksize, m - 1]\) - th column of the sub - matrix H i ; the P1 and P3 regions are further divided into multiple sub - regions according to blocksize, and the thread group Group i successively performs the following operations on any sub - region P1 t , P3 t , \(0\leq t\lt k - j\):

[0040] Each GPU thread T i in the thread group Group y , \(0\leq y\lt blocksize\), performs the following operations:

[0041] Read the column vector A with the serial number y in the sub - region P1 t into the cache of the thread;

[0042] Calculate each element in the column vector D with the serial number y in the sub - region P3 t :

[0043] D[x]=D[x] - C x A(4)

[0044] where \((j + 1)\times blocksize\leq x\lt m\).

[0045] Furthermore, in step S35, determine whether the currently processed block B j,j is the bottom - right block of the sub - matrix H i . If so, the batch parallel LU decomposition ends; otherwise, let \(j = j + 1\), jump to step S31, and continue to process the next block.

[0046] The batch parallel LU decomposition method for medium - sized dense matrices based on GPU provided by the present invention has the following beneficial effects compared with the existing methods:

[0047] (1) The present invention provides high-performance batch parallel LU decomposition capabilities for small and medium-sized dense matrices. The performance of this method increases approximately linearly with the increase in batch size, with the highest peak approaching 450 Gflops / s. Compared with the CUBLAS library function provided by NVIDIA, the highest speedup ratio exceeds 10;

[0048] (2) The present invention is of great significance in improving the efficiency of complex system simulation. When implicitly solving a model containing a large number of small ordinary differential equations, it can speed up the solution and effectively shorten the simulation time. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0050] Figure 1 A schematic diagram of a process flow of a batch parallel LU decomposition method according to an embodiment of the present invention;

[0051] Figure 2 A schematic diagram of a matrix H and its internal submatrices according to an embodiment of the present invention;

[0052] Figure 3 A schematic diagram of internal sub-matrix block division according to an embodiment of the present invention;

[0053] Figure 4 A schematic diagram of selecting a principal element in steps S312, S313, and S314 of an embodiment of the present invention;

[0054] Figure 5 A schematic diagram of the P1, P2, and P3 regions processed in the calculation process of steps S33, S34, and S35 of an embodiment of the present invention;

[0055] Figure 6 A schematic diagram of exchanging pivot columns in step S32 of an embodiment of the present invention;

[0056] Figure 7 Schematic diagram of processing the internal elements of the P2 region in step S34 of one embodiment of the present invention;

[0057] Figure 8 It is a schematic diagram of processing the internal elements of the P3 region in step S35 of one embodiment of the present invention;

[0058] Fig. 9 Schematic diagram of the experimental test effect. DETAILED DESCRIPTION

[0059] In order to make the purpose, technical solution and advantages of the embodiments of the present invention clearer, the technical solution in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0060] In order to improve the efficiency of batch parallel LU decomposition of small and medium-sized dense matrices, the present invention proposes a batch parallel LU decomposition method of small and medium-sized dense matrices based on GPU for simulation system. The basic idea of ​​the present invention is to use memory merge access and local variables to cache part of the data of the sub-matrix, use out-of-order execution technology to update part of the key data, postpone the vector-based outer multiplication operation, and convert it into a matrix-based dense multiplication operation, thereby effectively reducing the number of global memory accesses, increasing the proportion of effective computing time, and greatly improving the decomposition efficiency.

[0061] The present invention provides a batch LU decomposition method for small and medium-sized dense matrices based on GPU for simulation systems, such as Figure 1 As shown, including:

[0062] S1, calculate the Jacobian matrix of the model unit of the simulation system at each time step, and use these matrices as sub-matrices and concatenate them into a large matrix H;

[0063] S2, divide each sub-matrix into blocks, and allocate a corresponding number of GPU threads to each sub-matrix according to the size of the blocks to form a thread group;

[0064] S3. The GPU thread group completes the pivot selection of all columns in each diagonal block, updates the blocks to the right of the row where the diagonal block is located, exchanges the columns inside the diagonal block with the pivot columns, and updates all block data below and to the right of the column where the diagonal block is located, in diagonal order from the upper left to the lower right.

[0065] In step S1, the model can be directly described in a certain language, such as an expression in C language, and the Jacobian matrix of the model unit at each time step can be automatically calculated by symbolic differentiation. For complex system simulation, in its solution process, the Newton method is used for the n model units that constitute the system, and the Jacobian matrix of each unit is calculated when each time step arrives, and these sub-matrices are spliced ​​into a two-dimensional large matrix H. The subsequent steps perform batch parallel LU decomposition on the large matrix H.

[0066] The Jacobian matrix of each unit is an m-order small and medium-sized dense matrix. The m-order small and medium-sized dense matrix of n model units is used as the submatrix H i, spliced ​​into a large matrix H as input; the range of i is [0,n-1].

[0067] Specifically, the method for obtaining m is: expand the dimensions of all sub-simulation models to be solved to an integer power of 2 m (that is, the Jacobian matrix is ​​expanded to the mth order). For the added rows and columns, except for the diagonal elements whose values ​​are 1, all others are filled with 0. All the expanded sub-matrices are arranged in one dimension, and finally form a two-dimensional large matrix H with a dimension of m×N, where N=m×n, such as Figure 2 shown.

[0068] In step S2, for each submatrix H i The internal blocks are divided into blocks, and each sub-matrix H is divided into blocks according to the size of the blocks. i Allocate the corresponding number of GPU threads to form a thread group Group i ;

[0069] When the order m of the submatrix is ​​greater than or equal to 16, each submatrix is ​​divided into blocks. There are two options for the block size, blocksize, which can be 8 or 16. It can be determined by passing in parameters during the actual solution. If the dimension of the submatrix is ​​less than 16, the entire submatrix is ​​taken as one block, and blocksize = m. After the division, each submatrix contains k×k subblocks, where k = m÷blocksize. The block results are as follows: Figure 3 By dividing the submatrix into blocks, it can be ensured that the thread group Group i The data in the entire block can be read from the global memory at one time, effectively reducing the number of global memory accesses.

[0070] In step S3, LU decomposition is performed on all sub-matrices in the large matrix H. The specific steps include:

[0071] S31, each thread group Group i Select submatrix H i The diagonal block B with inner index j j,j , determine B j,j The pivot columns corresponding to all columns in B j,j All blocks on the right row are updated; j is initialized to 0;

[0072] S32, all the obtained pivot columns and diagonal blocks B j,j The corresponding columns inside are swapped;

[0073] S33, opposite to the diagonal block B j,j All blocks below the column are updated;

[0074] S34, opposite to the diagonal block Bj,j Update all the sub-blocks in the lower right corner;

[0075] S35. If the diagonal block B j,j is the lower right corner sub-block in the sub-matrix H i , the batch parallel LU decomposition ends, and the large matrix H is output. At this time, each sub-matrix H i in it is composed of L i and U i , which are the unit lower triangular matrix and the upper triangular matrix respectively, and are used for the subsequent fast solution of the linear equations; otherwise, let j = j + 1, jump to step S31, and continue to process the next sub-block.

[0076] In step S31, for the diagonal block B j,j to be processed currently, where 0 ≤ j < k, traverse all the diagonal elements b j,j in B p,p (j×blocksize ≤ p < (j + 1)×blocksize), and determine the column number of the corresponding pivot element. Specifically:

[0077] S311: Initialize a one-dimensional matrix To i with dimension m, and let To i [l] = m, where 0 ≤ l < m. Here,

[0078] the array index l corresponds to the column number of the sub-matrix H i , and To i [l] is used to record the column with serial number l and its final pivot column number;

[0079] S312: The thread group starts from the current block B j,j , and reads all the elements in the p-th row of the block in parallel and stores them in the cache inside the thread; where p is calculated according to the row index of the sub-matrix H i and is initialized to j×blocksize; in this step, since the row elements of the block B j,j are closely arranged in the global memory, the GPU only needs to issue one instruction to realize the reading of all the data of the thread group, achieving the effect of memory coalesced access;

[0080] S313: The thread group continues to read all the data in the p-th row of the next adjacent sub-block on the right side of B j,j in parallel, and compares it with the data in the cache, and saves the element with the larger absolute value and its column number in the internal cache. This process continues iteratively until the data of the last block of the sub-matrix is read; this step also achieves the effect of memory coalesced access;

[0081] S314: Compare all elements cached in the thread group, take the column q corresponding to the element with the largest absolute value as the main element column corresponding to column p, and update the element To i [p] = q;

[0082] The process of steps S312, S313 and S314 is as follows: Figure 4 As shown;

[0083] S315: Pair matrix H i The current row p in the , starting from column index j × blocksize, traverses all subsequent elements in the row If To i [l]>q, then the following operation is performed:

[0084]

[0085] in For elements The corresponding true pivot is Figure 4 shown.

[0086] This step updates the p-row elements for the subsequent elimination process;

[0087] S316: For Figure 5 Eliminate the P1 region in the

[0088] Specifically, the P1 region is located in the submatrix H i A matrix with [p+1,(j+1)×blocksize-1] rows and [j×blocksize,m-1] columns, for each element The column index y of To i [y]>q, then the following elimination operation is performed:

[0089]

[0090] where p+1≤x<(j+1)×blocksize,j×blocksize≤y <m

[0091] S317: If p is block B j,j The last row in the step S31 ends, and B j,j Contains the LU decomposition result L in the current block j,j and U j,j ; Otherwise, set p=p+1, jump to step S312, and continue processing the next row.

[0092] In step S32, the thread group Group i At the same time, the current diagonal block B j,jAll columns in and their corresponding pivot columns are swapped. Figure 6 As shown, for B j,j Any column p in the thread group, the pth thread in the thread group executes column col p With column col To i [p] Exchange of elements between

[0093] Since the diagonal block B j,j All columns in the same column are swapped at the same time, and all elements in the column interval [j×blocksize, (j+1)×blocksize-1] can be read and written in memory, which effectively reduces the number of GPU memory instruction accesses.

[0094] In step S33, for Figure 5 The P2 area in the .

[0095] The P2 region is located in the submatrix H i For each element of the matrix between [(j+1)×blocksize,m-1] rows and [j×blocksize,(j+1)×blocksize-1] columns like Figure 7 As shown, perform the following operations:

[0096]

[0097] Where (j+1)×blocksize≤x <m,j×blocksize≤y<(j+1)×blocksize

[0098] From formula (3), we can see that there is a causal relationship between the columns in P2. The update of the latter column depends on the results of all the previous columns. Therefore, it cannot be completely parallelized and needs to be updated separately. In this way, it is decoupled from the subsequent P3 area, so that the update of the P3 area can be completely parallelized.

[0099] In step S34, for Figure 5 The P3 area in the .

[0100] The P3 region is located in the submatrix H i The P1 and P3 regions are further divided into multiple sub-regions according to blocksize, such as Figure 8 As shown; Thread Group Group i Process these sub-regions in turn, and for any sub-region P1 t ,P3t , when \(0\leq t < k - j\), perform the following operations:

[0101] Thread group Group i For each GPU thread T y , when \(0\leq y < blocksize\), perform the following operations:

[0102] Read in sub-region P1 t The column vector A with serial number y in it into the cache of the thread;

[0103] Calculate sub-region P3 t For each element in the column vector D with serial number y in it:

[0104] D[x] = D[x] - C x A(4)

[0105] where \((j + 1)\times blocksize\leq x < m\);

[0106] It can be seen from formula (4) that the calculation of each element in column D requires the use of vector A. Since its data has been read into the local cache of the thread, multiple reads are effectively avoided, improving the proportion of computing time; at the same time, since all the data in regions P1 and P2 have been processed before, when the current step is executed, the matrix multiplication and addition operations can be completely parallelized in a block form, which can greatly improve the computing efficiency.

[0107] In step S35, determine whether the currently processed block B j,j is the bottom-right block of sub-matrix H i , if so, the batch parallel LU decomposition ends, and sub-matrix H i contains the decomposition results L i and U i ; otherwise, let \(j = j + 1\), and jump to step S3 to continue processing the next block.

[0108] When the method ends, output the large matrix H. At this time, each sub-matrix H i is composed of L i and U i in two parts, which are the unit lower triangular matrix and the upper triangular matrix respectively, and are used for subsequent fast solution of batch linear equations, and then calculate the state variable values of all simulation units through the Newton iteration method to ensure the rapid progress of the simulation process.

[0109] According to the method of the present invention, experimental tests were carried out on the batch parallel LU decomposition of dense matrices, as Fig. 9As shown: the horizontal axis is the dimension of the large matrix, the vertical axis is the data processing capability, and the different sub-matrix dimensions are marked in the lower right corner, represented by the curves of corresponding colors in the figure. It can be seen that, whether it is divided into 8×8 blocks or 16×16 blocks, the LU decomposition processing capability basically increases linearly with the number of sub-matrices, which means that the GPU hardware is fully utilized during the execution of the method of the present invention; it can also be seen from the figure that when the dimension of the large matrix reaches 8192×8192, the peak floating-point operation processing capability reaches about 450Gflops / s. Table 1 shows some experimental comparison results of the execution efficiency of the method of the present invention and the CUBALS library function officially provided by NVIDIA. When the sub-matrix dimension is 64, as the number of sub-matrices increases, the advantage of the method of the present invention gradually increases, and the highest acceleration ratio reaches about 13 times. Similar acceleration ratio results are also obtained under other sub-matrix dimensions.

[0110] Table 1 Comparison of processing speed results between the method of the present invention and CUBLAS

[0111]

[0112] While the present invention is highly parallelized, it fully utilizes the hardware features of memory merge access and local variable cache, effectively hides the latency of memory access, and improves the proportion of computing time. Therefore, the method of the present invention has high computational efficiency, can effectively reduce the time of LU decomposition, and has the implicit real-time parallel solution capability to support millions of small simulation models.

[0113] Those skilled in the art will appreciate that the foregoing embodiments are merely intended to illustrate the technical solutions of the present invention rather than to limit the present invention. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art will appreciate that the technical solutions described in the foregoing embodiments may be modified or some or all of the technical features may be replaced by equivalents. Such modifications or replacements do not deviate the essence of the corresponding technical solutions from the scope defined by the claims of the present invention.

Claims

1. A GPU-based batch LU decomposition method for small and medium-sized dense matrices for simulation systems, characterized in that: include: S1, calculate the Jacobian matrix of the model unit of the simulation system at each time step, and use these matrices as sub-matrices and concatenate them into a large matrix H; S2, divide each sub-matrix into blocks, and allocate GPU thread groups to each sub-matrix according to the size of the blocks; S3, the GPU thread group sequentially completes the pivot selection of all columns in each diagonal block, the update of the right block of the row where the diagonal block is located, the exchange of the internal column of the diagonal block with the pivot column, and the update of all block data below and to the right of the column where the diagonal block is located, thereby completing the LU decomposition of all sub-matrices in the large matrix H; The simulation system has n model units, and the Jacobian matrix is ​​a small and medium-sized dense matrix H i , the range of i is [0, n-1]; In step S3, the specific steps include: S31, each thread group Group i Select submatrix H i The diagonal block B with inner index j j,j , determine B j,j The pivot columns corresponding to all columns in B j,j All blocks on the right row are updated; j is initialized to 0; S32, all the obtained pivot columns and diagonal blocks B j,j The corresponding columns inside are swapped; S33, opposite to the diagonal block B j,j All blocks below the column are updated; S34, opposite to the diagonal block B j,j All blocks in the lower right corner are updated; S35, if the diagonal block B j,j is the submatrix H i The lower right corner of the block, batch parallel LU decomposition is completed, and the large matrix H is output, where each submatrix H i All by L i and U i The matrix is ​​composed of a unit lower triangular matrix and an upper triangular matrix, which are used to quickly solve the subsequent linear equations; otherwise, let j=j+1 and jump to step S31 to continue processing the next block.

2. The batch LU decomposition method according to claim 1, characterized in that: In step S1, the Jacobian matrix is ​​of order m, and m is obtained by expanding the dimensions of all model units currently to be solved upward to an integer power of 2 m, and for the added rows and columns, except for the diagonal elements whose values ​​are 1, all others are filled with 0.

3. The batch LU decomposition method according to claim 1, characterized in that: In step S2, it includes: S21, sub-matrix H i Divide into several sub-blocks to ensure that in subsequent processing, the thread group can read the data in the entire block from the global memory at one time, reducing the number of global memory accesses; S22, according to the size of the block for the sub-matrix H i Allocate the corresponding number of GPU threads to form a thread group Group i .

4. The batch LU decomposition method according to claim 3, characterized in that: Sub-matrix H i Methods for dividing into several sub-blocks include: If m is greater than or equal to 16, set the block size blocksize to 8 or 16; If m is less than 16, the entire submatrix is ​​taken as a block, and the block size blocksize=m.

5. The batch LU decomposition method according to claim 1, characterized in that: Step S31 includes: S311: Initialize a one-dimensional matrix To with dimension m i , and let To i [l]=m, 0 ≤ l < m, where Array index l corresponds to submatrix H i Column number, To i [l] is used to record the final pivot column number of the column with serial number l; S312: The thread group starts from the current block B j,j At the beginning, all elements of the pth row in the block are read in parallel and stored in the cache inside the thread; where p is based on the submatrix H i The row index calculation is initialized to j×blocksize; in this step, since block B j,j The row elements are closely arranged in the global memory, and the GPU only needs to issue one instruction to read all the data of the thread group, achieving the effect of memory merged access; S313: The thread group continues to read B in parallel j,j All the data in the next block on the right, p rows, are read and compared with the data in the cache, and the elements with larger absolute values ​​and their column numbers are saved in the internal cache; this process is iterated continuously until the last block of data in the submatrix is ​​read; S314: Compare all elements cached in the thread group, take the column q corresponding to the element with the largest absolute value as the main element column corresponding to column p, and update the element To i [p] = q; S315: Pair matrix H i The current row p in the , starting from column index j × blocksize, traverses all subsequent elements in the row If To i [l]>q, then the following operation is performed: in, For elements The corresponding true pivot; S316: For each element in the P1 region The column index y of To i [y]>q, then the following elimination operation is performed: Among them, the P1 region is located in the submatrix H i A matrix with [p+1, (j+1) × blocksize-1] rows and [j × blocksize, m-1] columns; p+1≤x<(j+1) × blocksize, j × blocksize≤y <m; S317: If p is block B j,j The last row in the step S31 ends, and B j,j Contains the LU decomposition result L in the current block j,j and U j,j ; Otherwise, set p=p+1, jump to step S312, and continue processing the next line.

6. The batch LU decomposition method according to claim 1, characterized in that: In step S32, for B j,j Any column p in the thread group, the pth thread in the thread group executes column col p With column col To i [p] Exchange elements between .

7. The batch LU decomposition method according to claim 1, characterized in that: In step S33, let P2 region be located in submatrix H i For each element of the matrix between [(j+1)×blocksize,m-1] rows and [j×blocksize,(j+1)×blocksize-1] columns Perform the following operations: Where (j+1)×blocksize≤x <m,j×blocksize≤y<(j+1)×blocksize。 8. The batch LU decomposition method according to claim 1, characterized in that: In step S34, let the P3 region be the matrix between the [(j + 1)×blocksize, m - 1] row and the [(j + 1)×blocksize, m - 1] column of the sub-matrix H i ; the P1 and P3 regions are further divided into multiple sub-regions according to blocksize, and the thread group Group i successively processes any of the sub-regions P1 t , P3 t , 0 ≤ t < k - j, and performs the following processing: Thread group Group i Each GPU thread T y , 0 ≤ y < blocksize, performs the following operations: Read in sub-area P1 t The y column vector A with sequence number in is put into the cache of the thread; Calculate sub-area P3 t Each element in the y column vector D with sequence number: D[x]=D[x]-C x A (4) Where (j+1)×blocksize≤x <m。 9. The batch LU decomposition method according to claim 1, characterized in that: In step S35, determine the current processed block B j,j Is the submatrix H i If the lower right corner block is found, the batch parallel LU decomposition ends; otherwise, let j=j+1 and jump to step S31 to continue processing the next block.

Images